C ______________________________ POWER SPECTRUM C 25 april 1990 C R. Kates (Potsdam), A. Klypin (NMSU) C NROW = number of particles per row C NGRID= number of cells on a side (normal grid) C C Warning: comment out lines with declarations of C common blocks /FOURAR/ and /F67COM/ C in files PMparameters.h C INCLUDE 'PMparameters.h' INCLUDE 'PMmesh.h' REAL INPUT C................................................................... C Read data and open files OPEN(17,FILE='Spectrum.DAT',STATUS='UNKNOWN') CALL RDTAPE write (*,*) ' RDTAPE is done' WRITE (17,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,wspecies(1), + EKIN, + NROWC,NGRID,NRECL,Om0,Oml0,hubble, + Ocurv,lspecies(1),Nspecies 100 FORMAT(1X,'Header=>',A45,/ + 1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/ + 1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.3,/ + 1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I9,/ + 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/ + 1x,' Omega_curvature=',F7.3,/ + 1x,' Number of part.1specie=',i10,' Nspecies=',i3) Box = extras(100) If(Box.le.1.e-3)Then Box =INPUT(' Enter box size in Mpc/h =') EndIf C CALL DENSIT ! Find density write (*,*) ' Density is ok' c Call DensCR CALL DENTES(DELR) ! Statistics CALL POWER(Box) ! Power spectrum STOP END C----------------------------------------------------- C Compute mean density and rms SUBROUTINE DENTES(DELR) C----------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMmesh.h' Dimension x(44),amass(44),imass(44) real*8 SUM,SUM2 Data x/ 0., 0.2, 0.4, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, . 1.3, 1.5, 1.75, 2., 2.25, 2.5, 3.0, 3.5, 4., 5., . 7., 10., 12., 15., 20., 25., 30., 40., 50., 70., . 100., 120.,150.,200.,300., 500., 700.,1000.,1200.,1500., . 2.e3, 3.e3,5.e3,1.e4/ Do i=1,44 amass(i) =0. imass(i) =0 EndDo SUM = 0. SUM2= 0. Nn = 0 Am = 0. C.... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE ( k,j,i) C$OMP+REDUCTION(+:SUM,SUM2) DO K=1,NGRID DO J=1,NGRID DO I=1,NGRID SUM = SUM + FI(I,J,K) SUM2= SUM2+ FI(I,J,K)**2 c Do l=44,1,-1 c Dens =FI(I,J,K) +1. c If(Dens .GT. x(l))Then c Nn =Nn +1 c Am =Am +Dens c imass(l) =imass(l) +1 c amass(l) =amass(l) +Dens c GoTo 10 c EndIf c EndDo c 10 Continue ENDDO ENDDO ENDDO Total =(FLOAT(NGRID))**3 DENM = SUM/Total DELR = DSQRT(SUM2/Total-DENM**2) DensAv=PARTW*float(NROW)**3/(float(NGRID))**3 DensOne= float(NGRID)**3/float(NROW)**3 WRITE (*,150) DELR,DENM,DensAv,DensOne WRITE (17,150) DELR,DENM,DensAv,DensOne 150 format(20('-'),' Density Distribution',20('-'),/ & ' Density is in units of average density', & ' in the Box',/20x, & ' RMS Delta Rho/Rho =',G11.4,/20x, & ' Mean Delta Rho/Rho =',G11.4,/20x, & ' Average density =',G11.4,/20x, & ' One particle density=',G11.4) c &,/' bin Rho Rho Sum(rho_i)/dRho/TotMass Sum',/ c &'
< left > Normalized density Number of cells ',/ c &' < bin > Distribution With this density') c Do k=1,44-1 c drho =x(k+1)-x(k) c rhom =x(k)+drho/2. c Sn =drho*Total c If(Sn .lt.1.e-7)Sn=1. c write (*,20) k,rhom,x(k),amass(k)/Sn,imass(k) c write (17,20) k,rhom,x(k),amass(k)/Sn,imass(k) c 20 Format(i3,2f8.1,2x,G12.3,12X,i12) c EndDo c k =44 c write (*,20) k,x(k),x(k),amass(k),imass(k) c write (17,20) k,x(k),x(k),amass(k),imass(k) RETURN END C---------------------------------------------------- C power spectrum P(k) C for given density field FI SUBROUTINE POWER(Box) C---------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMmesh.h' PARAMETER ( Pi=3.141592653, NPOWER =NGRID) PARAMETER (NPROCS =64) DIMENSION DPOWER(NPOWER),iVolume(NPOWER),dk(NPOWER) DIMENSION DPW(NPOWER,NPROCS),iVol(NPOWER,NPROCS), & dkk(NPOWER,NPROCS) dimension Zf(narr) , Yf(narr) dimension si(nf67) , indx(nf67) REAL*8 SIGMA,DPOWER,DPW NSPEC = NGRID/2 If(Box.le.0..or.Box.gt.1.e+4)Then write (*,*) ' Wrong Box size=',Box,' STOP' Return EndIf SIGMA =0. Scale = 2.*Pi/Box DO M =1,NPOWER iVolume(M) =0 DPOWER(M) =0. dk(M) =0. Do i=1,NPROCS iVol(M,i) =0 DPW(M,i) =0. dkk(M,i) =0. EndDo ENDDO ib1 = 3 iq = int(alog(float(NGRID))/alog(2.)+0.5) FT = (2.*NGRID)**3*float(NGRID)**3 C FFT of density field write (*,*) ' Do fft of rho along x-axis' call setf67(ib1,iq,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) C.... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE ( k,j,i,zf,yf,ip,isl,l1,n2,n4 ) DO K=1,NGRID DO J=1,NGRID DO I=1,NGRID Zf(I) =FI(I,J,K) ENDDO call four67(ib1,iq,ip,isl,l1,n2,n7, & si,indx,Zf,Yf) DO I=1,NGRID FI(I,J,K) = Yf(I) ENDDO c If(J.lt.12.and.K.eq.1)write (*,'(2i4,16g12.3)')J,K,(Zf(I),I=1,16) ENDDO ENDDO C.... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE ( k,j,i,zf,yf,ip,isl,l1,n2,n4 ) DO K=1,NGRID DO I=1,NGRID DO J=1,NGRID Zf(J) =FI(I,J,K) ENDDO call four67(ib1,iq,ip,isl,l1,n2,n7,si,indx, & Zf,Yf) DO J=1,NGRID FI(I,J,K) = Yf(J) ENDDO ENDDO ENDDO write (*,*) ' Swap k<->i' C.... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE ( k,j,i,aa) DO J=1,NGRID DO K=1,NGRID-1 DO I=K+1,NGRID aa = FI(I,J,K) FI(I,J,K) =FI(K,J,I) FI(K,J,I) =aa ENDDO ENDDO ENDDO write (*,*) ' ....z' C.... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE ( k,j,i,zf,yf,ip,isl,l1,n2,n4 ) DO K=1,NGRID DO J=1,NGRID DO I=1,NGRID Zf(I) =FI(I,J,K) ENDDO call four67(ib1,iq,ip,isl,l1,n2,n7,si,indx, & Zf,Yf) DO I=1,NGRID FI(I,J,K) = Yf(I)**2/FT ENDDO ENDDO ENDDO FI(1,1,1) =0. kStep =NGRID/NPROCS If(kstep*NPROCS.ne.NGRID)STOP '-- Wrong number of NPROCS' write (*,*) ' Power .... kStep=',kStep C.... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (iP,kStart,kEnd, k,j,i,MJ,MK,AMP,RK,MP) C$OMP+REDUCTION(+:SIGMA) Do iP=1,NPROCS kStart =(iP-1)*kStep+1 kEnd =iP*kStep Do K=kStart,kEnd IF(K.LE.(NSPEC+1)) THEN MK =(K-1)**2 ELSE MK =(K-1-NSPEC)**2 ENDIF DO J=1,NGRID IF(J.LE.(NSPEC+1)) THEN MJ =(J-1)**2 ELSE MJ =(J-1-NSPEC)**2 ENDIF DO I=1,NGRID AMP = FI(I,J,K) SIGMA = SIGMA +AMP C Find Abs(K) IF(I.LE.(NSPEC+1)) THEN RK =SQRT( FLOAT((I-1)**2 +MK+MJ) ) ELSE RK =SQRT( FLOAT((I-1-NSPEC)**2+MK+MJ) ) ENDIF C Power Spectrum MP =INT(RK +0.5) IF(MP.LE.NPOWER.and.MP.gt.0)THEN DPW(MP,iP) =DPW(MP,iP) +AMP iVol(MP,iP) =iVol(MP,iP)+1 dkk(MP,iP) =dkk(MP,iP) +RK ENDIF ENDDO ENDDO ENDDO ENDDO Do MP=1,NPOWER Do iP=1,NPROCS DPOWER(MP) =DPOWER(MP) +DPW(MP,iP) iVolume(MP) =iVolume(MP)+iVol(MP,iP) dk(MP) =dk(MP) +dkk(MP,iP) EndDo EndDo WRITE (*,100) SQRT(SIGMA),Scale*NROW/2. WRITE (17,100) SQRT(SIGMA),Scale*NROW/2. 100 FORMAT(5x,'RMS DRho/Rho=',F7.3,' k_Nyquist for Particles=',f6.2, & /4X,'n',2X,'k(h/Mpc)',1X, & 'Harmonics ','Power Spectrum(h^3)') 110 FORMAT(1x,I4,F8.3,I10,3X,5G12.4) DO I=1,NPOWER IF(iVolume(I).GT.0)THEN wk = dk(I)/iVolume(I)*Scale Pk = Box**3*DPOWER(I)/iVolume(I) WRITE (*,110) I,wk,iVolume(I),Pk WRITE (17,110) I,wk,iVolume(I),Pk ENDIF ENDDO RETURN END C********************************************************************* C Corr.function for given R FUNCTION FCORR(R) INCLUDE 'PMparameters.h' INCLUDE 'PMmesh.h' PARAMETER ( MMAX=NMAX-2) PARAMETER (PI=3.141592653) IF(R.LE.1.E-5)THEN WRITE (*,*) ' Fcorr(R): R can not be zero' STOP ENDIF QF =2.*PI*R/FLOAT(NGRID) QSHF =(FLOAT(NGRID)/2.)**2 M =NMAX +1 SIGMA =0. FK =SQRT(3.)*NMAX SIGMA =SIGMA + SIN(QF*FK) * FI(M,M,M) FK =SQRT(2.)*NMAX SIGMA =SIGMA + SIN(QF*FK) *(FI(M,M,1)+FI(M,1,M)+FI(1,M,M)) FK = NMAX SIGMA =SIGMA + SIN(QF*FK) *(FI(M,1,1)+FI(1,M,1)+FI(1,1,M)) DO K1 =1,MMAX M1 =K1 +1 L1 =K1 +1+NMAX K12 =K1**2 SIGMA =SIGMA + SIN(QF*K1)*( FI(M1,1,1)+FI(L1,1,1)+ . FI(1,M1,1)+FI(1,L1,1)+ . FI(1,1,M1)+FI(1,1,L1)) FK =SQRT(K12+2.*QSHF) SIGMA =SIGMA + SIN(QF*FK)*( FI(M1,M,M)+FI(L1,M,M)+ . FI(M,M1,M)+FI(M,L1,M)+ . FI(M,M,M1)+FI(M,M,L1)) FK =SQRT(K12+ QSHF) SIGMA =SIGMA + SIN(QF*FK)*( FI(M1,M,1)+FI(L1,M,1)+ . FI(M,M1,1)+FI(M,L1,1)+ . FI(M1,1,M)+FI(L1,1,M)+ . FI(M,1,M1)+FI(M,1,L1)+ . FI(1,M1,M)+FI(1,L1,M)+ . FI(1,M,M1)+FI(1,M,L1)) DO K2 =1,MMAX M2 =K2 +1 L2 =K2 +1+NMAX K22 =K2**2 FK =SQRT(Float(K12+ K22)) SIGMA =SIGMA + SIN(QF*FK)*( FI(M1,M2,1)+FI(M1,L2,1)+ . FI(L1,M2,1)+FI(L1,L2,1)+ . FI(1,M1,M2)+FI(1,M1,L2)+ . FI(1,L1,M2)+FI(1,L1,L2)+ . FI(M1,1,M2)+FI(M1,1,L2)+ . FI(L1,1,M2)+FI(L1,1,L2)) FK =SQRT(K12+ K22 +QSHF) SIGMA =SIGMA + SIN(QF*FK)*( FI(M1,M2,M)+FI(M1,L2,M)+ . FI(L1,M2,M)+FI(L1,L2,M)+ . FI(M,M1,M2)+FI(M,M1,L2)+ . FI(M,L1,M2)+FI(M,L1,L2)+ . FI(M1,M,M2)+FI(M1,M,L2)+ . FI(L1,M,M2)+FI(L1,M,L2)) DO K3 =1,MMAX M3 =K3 +1 L3 =K3 +1+NMAX K32 =K3**2 FK =SQRT(Float(K12+ K22 +K32)) SIGMA =SIGMA + SIN(QF*FK)*( FI(M1,M2,M3)+FI(M1,L2,M3)+ . FI(L1,M2,M3)+FI(L1,L2,M3)+ . FI(M1,M2,L3)+FI(M1,L2,L3)+ . FI(L1,M2,L3)+FI(L1,L2,L3)) ENDDO ENDDO ENDDO FCORR = SIGMA/QF RETURN END C--------------------------------------- C Density Field SUBROUTINE DENSIT C--------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMmesh.h' XN =FLOAT(NGRID)+1.-1.E-7 YN =FLOAT(NGRID) C Subtract mean density DO M3=1,NGRID DO M2=1,NGRID DO M1=1,NGRID FI(M1,M2,M3) = -1. END DO END DO END DO W = (FLOAT(NGRID)/FLOAT(NROW))**3 PARTW = W WPAR = wspecies(1) N =0 xmin =1.e+5 xmax =-1.e+5 ymin =1.e+5 ymax =-1.e+5 zmin =1.e+5 zmax =-1.e+5 C Loop over particles If(Nspecies.eq.0)Then ! old case: 1 complete set Npages =NROW ! Number of pages N_in_last=NPAGE ! Number of particles in last page Else ! multiple masses N_particles =lspecies(Nspecies) ! Total number of particles Npages = (N_particles -1)/NPAGE +1 N_in_last=N_particles -NPAGE*(Npages-1) EndIf write (*,*) Nspecies, + lspecies(Nspecies),wspecies(Nspecies) Nsp_current =1 ff = float(NGRIDC) ind = INT(log(float(NGRID)/ff)/log(2.)+10.)-10 WPAR = wspecies(1) *8.**ind Lspec_next=lspecies(1)+1 DO IROW=1,Npages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last If(IROW/50*50.eq.IROW) + write (*,*) ' Read page=',IROW,' N=',N, + ' Npages=',Npages,In_page iL = NPAGE*(IROW-1) C Loop over particles READ (21,REC=IROW) RECDAT DO IN=1,In_page N = N +1 If(N.ge.Lspec_next)Then ! next specie write (*,*) ' Another specie. Current=',Nsp_current write (*,*) ' Particle=',N,WPAR, & ' lspecies=',lspecies(Nsp_current), & ' next=',lspecies(Nsp_current+1) Nsp_current=Nsp_current+1 Lspec_next =lspecies(Nsp_current)+1 WPAR = wspecies(Nsp_current) EndIf c X=XPAR(IN) c Y=YPAR(IN) c Z=ZPAR(IN) X=(XPAR(IN)-1.)*NGRID/ff +1. Y=(YPAR(IN)-1.)*NGRID/ff +1. Z=(ZPAR(IN)-1.)*NGRID/ff +1. c X= X+0.499 IF(X.ge.NGRID+1.)X=X-NGRID c Y= Y+0.499 IF(Y.ge.NGRID+1.)Y=Y-NGRID c Z= Z+0.499 IF(Z.ge.NGRID+1.)Z=Z-NGRID xmin =MIN(xmin,X) ymin =MIN(ymin,Y) zmin =MIN(zmin,Z) xmax =MAX(xmax,X) ymax =MAX(ymax,Y) zmax =MAX(zmax,Z) I=INT(X) J=INT(Y) K=INT(Z) If(I.le.0)write (*,*) ' X:',X,Y,Z,I,' Irow=',IROW,IN,ifile If(J.le.0)write (*,*) ' Y:',X,Y,Z,I,' Irow=',IROW,IN,ifile If(K.le.0)write (*,*) ' Z:',X,Y,Z,I,' Irow=',IROW,IN,ifile D1=X-FLOAT(I) D2=Y-FLOAT(J) D3=Z-FLOAT(K) T1=1.-D1 T2=1.-D2 T3=1.-D3 T2W =T2*WPAR D2W =D2*WPAR I1=I+1 IF(I1.GT.NGRID)I1=1 J1=J+1 IF(J1.GT.NGRID)J1=1 K1=K+1 IF(K1.GT.NGRID)K1=1 FI(I ,J ,K ) =FI(I ,J ,K ) +T3*T1*T2W FI(I1,J ,K ) =FI(I1,J ,K ) +T3*D1*T2W FI(I ,J1,K ) =FI(I ,J1,K ) +T3*T1*D2W FI(I1,J1,K ) =FI(I1,J1,K ) +T3*D1*D2W FI(I ,J ,K1) =FI(I ,J ,K1) +D3*T1*T2W FI(I1,J ,K1) =FI(I1,J ,K1) +D3*D1*T2W FI(I ,J1,K1) =FI(I ,J1,K1) +D3*T1*D2W FI(I1,J1,K1) =FI(I1,J1,K1) +D3*D1*D2W ENDDO ENDDO write (*,*) ' Number of particles read=', N write (*,*) ' X:',xmin,xmax write (*,*) ' Y:',ymin,ymax write (*,*) ' Z:',zmin,zmax RETURN END C------------------------------------ c Subroutine densCR c INCLUDE 'PMparameters.h' c INCLUDE 'PMmesh.h' c XN =FLOAT(NGRID)+1.-1.E-7 c YN =FLOAT(NGRID) c SUMM = 0. c WD =0. c open(40,file='dens_cr.dat',status='unknown') c read(40,*) FI c Do k=1,NROW c Do j=1,NROW c Do i=1,NROW c SUMM =SUMM +FI(i,j,k) c WD =WD +FI(i,j,k)**2 c EndDo c EndDo c EndDo c sigma = sqrt(MAX(1.d-10,WD/NROW**3)) c write (*,*) ' Average density=',SUMM/NROW**3 c write (*,*) ' RMS density=',sigma c IFOUR = INT(ALOG(FLOAT(NROW))/ALOG(2.)+0.5) c Write (*,*) ' Ifour =',IFOUR,' Nrow=',NROW c Return c End C--------------------------------------------------- C Read current data from disk/tape, C Open files C Nrecl is the number of values in a record C Npage is the number of particles in a page SUBROUTINE RDTAPE C--------------------------------------------------- INCLUDE 'PMparameters.h' C Open file on a tape OPEN(UNIT=9,FILE='PMcrd.DAT', + FORM='UNFORMATTED', STATUS='UNKNOWN') C Read control information C and see whether it has proper C structure READ (9,err=20,end=20) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, + NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5 + ,Ocurv,extras WRITE (*,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN,EKIN1,EKIN2, + NROWC,NGRID,NRECL,Om0,Oml0,hubble, + Ocurv WRITE (16,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN,EKIN1,EKIN2, + NROWC,NGRID,NRECL,Om0,Oml0,hubble, + Ocurv 100 FORMAT(1X,'Header=>',A45,/ + 1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/ + 1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',3E12.3,/ + 1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I9,/ + 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/ + 1x,' Omega_curvature=',F7.3) IF(NROW.NE.NROWC) THEN WRITE (*,*) + ' NROW in PARAMETER and in TAPE-FILE are different' ENDIF IF(NGRID.NE.NGRIDC) THEN WRITE (*,*) + ' NGRID in PARAMETER and in TAPE-FILE are different:' write (*,*) ' Ngrid=',NGRID,' NgridC=',NGRIDC ENDIF C Open work file on disk write (*,*) ' start openning PMcrs0' 10 NBYTE = NRECL*4 NACCES= NBYTE / nbyteword OPEN(UNIT=21,FILE='PMcrs0.DAT',ACCESS='DIRECT', + STATUS='UNKNOWN',RECL=NACCES) REWIND 9 write (*,*) ' done openning PMcrs0' RETURN 20 write (*,*) ' Error reading the header file: Abort' stop END C---------------------------------- Read in variables REAL FUNCTION INPUT(text) C------------------------------------------------ Character text*(*) write (*,'(A,$)')text read (*,*) x INPUT =x Return End C-------------------------------------- Fourier Transform SUBROUTINE SETF67(IBC1,IQ1,ibc,ip,isl,l1,n2,n3,n4,n7, & si,indx,Zf,Yf) c ----------------------------------------------------- INCLUDE 'PMparameters.h' include 'PMmesh.h' dimension Zf(narr) , Yf(narr) dimension si(nf67) , indx(nf67) integer ibc , ip, isl , l1 , n2 , n3 , n4 , n7 c INCLUDE 'Control.param' PI=DATAN(1.D0)*4. IBC=IBC1 IQ=IQ1 IF(IBC.LT.3) GO TO 101 IQ=IQ-1 101 N3=2**IQ N7=N3/2 N5=N3/4 I=1 INDX(I)=N5 SI(I)=0.5*SQRT(2.) K=I I=I+1 102 IL=I IF(I.EQ.N7) GO TO 104 103 K1=INDX(K)/2 INDX(I)=K1 SI(I)=SIN(PI*FLOAT(K1)/FLOAT(N3)) K1=N7-K1 I=I+1 INDX(I)=K1 SI(I)=SIN(PI*FLOAT(K1)/FLOAT(N3)) K=K+1 I=I+1 IF(K.EQ.IL) GO TO 102 GO TO 103 104 RETURN END C----------------------------------------- SUBROUTINE TFOLD(IS,L,ZZZ,n2) INCLUDE 'PMparameters.h' include 'PMmesh.h' integer n2 DIMENSION ZZZ(NARR) IH2=N2/2-1 DO 100 I=IS,IH2 I1=I+L I2=N2-I+L A=ZZZ(I1) ZZZ(I1)=A-ZZZ(I2) ZZZ(I2)=A+ZZZ(I2) 100 CONTINUE RETURN END C------------------------------------ SUBROUTINE NEG(I1,I3,I2,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) INCLUDE 'PMparameters.h' include 'PMmesh.h' dimension Zf(narr) , Yf(narr) dimension si(nf67) , indx(nf67) integer ibc , ip, isl , l1 , n2 , n3 , n4 , n7 DO 100 K=I1,I3,I2 Yf(K)=-Yf(K) 100 CONTINUE RETURN END C------------------------------------- SUBROUTINE REVNEG(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) INCLUDE 'PMparameters.h' include 'PMmesh.h' dimension Zf(narr) , Yf(narr) dimension si(nf67) , indx(nf67) integer ibc , ip, isl , l1 , n2 , n3 , n4 , n7 DO 100 I=1,N7 J=N3+1+I K=N4+1-I A=Yf(J) Yf(J)=-Yf(K) Yf(K)=-A 100 CONTINUE RETURN END C------------------------------------- SUBROUTINE ZEERO(L,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) INCLUDE 'PMparameters.h' include 'PMmesh.h' dimension Zf(narr) , Yf(narr) dimension si(nf67) , indx(nf67) integer ibc , ip, isl , l1 , n2 , n3 , n4 , n7 DO 100 I=1,N2 LI=L+I Zf(LI-1)=0.0 100 CONTINUE RETURN END C------------------------------------------ SUBROUTINE TFOLD1(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) INCLUDE 'PMparameters.h' include 'PMmesh.h' dimension Zf(narr) , Yf(narr) dimension si(nf67) , indx(nf67) integer ibc , ip, isl , l1 , n2 , n3 , n4 , n7 II=N2-1 DO 100 I=1,II I1=ISL+I I2=L1-I A=Zf(I1) Zf(I1)=A+Zf(I2) Zf(I2)=A-Zf(I2) 100 CONTINUE RETURN END C------------------------------------- SUBROUTINE KFOLD(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) INCLUDE 'PMparameters.h' include 'PMmesh.h' dimension Zf(narr) , Yf(narr) dimension si(nf67) , indx(nf67) integer ibc , ip, isl , l1 , n2 , n3 , n4 , n7 JS1=N2 I=1 J5=ISL+N2 IS1=ISL IC1=L1 JS1=JS1/2 IF(JS1.NE.1) GO TO 200 K1=INDX(I) SN=SI(I) IS0=IS1 IS1=IS1+JS1 IC0=IC1 IC1=IC1+JS1 ODD1=SN*(Zf(IC1)-Zf(IS1)) Yf(K1+1)=Zf(IC0)+ODD1 N3MK1=N3-K1 Yf(N3MK1+1)=Zf(IC0)-ODD1 IF(IBC.LT.3) GO TO 110 ODD2=SN*(Zf(IC1)+Zf(IS1)) N3PK1=N3+K1 Yf(N3PK1+1)=Zf(IS0)+ODD2 N4MK1=N4-K1 Yf(N4MK1+1)=-Zf(IS0)+ODD2 110 RETURN 200 SN=SI(I) IS1=IS1+JS1 IC1=IC1+JS1 J3=IS1+JS1 210 IS0=IS1-JS1 IC0=IC1-JS1 ODD1=SN*(Zf(IC1)-Zf(IS1)) ODD2=SN*(Zf(IC1)+Zf(IS1)) Zf(IC1)=Zf(IC0)-ODD1 Zf(IS1)=-Zf(IS0)+ODD2 Zf(IC0)=Zf(IC0)+ODD1 Zf(IS0)=Zf(IS0)+ODD2 IS1=IS1+1 IC1=IC1+1 IF(IS1.NE.J3) GO TO 210 I=I+1 300 IS1=ISL IC1=L1 JS1=JS1/2 IF(JS1.EQ.1) GO TO 400 310 SN=SI(I) I=I+1 CS=SI(I) IS1=IS1+JS1 IC1=IC1+JS1 J3=IS1+JS1 320 IS0=IS1-JS1 IC0=IC1-JS1 ODD1=CS*Zf(IC1)-SN*Zf(IS1) ODD2=SN*Zf(IC1)+CS*Zf(IS1) Zf(IC1)=Zf(IC0)-ODD1 Zf(IC0)=Zf(IC0)+ODD1 Zf(IS1)=-Zf(IS0)+ODD2 Zf(IS0)=Zf(IS0)+ODD2 IS1=IS1+1 IC1=IC1+1 IF(IS1.NE.J3) GO TO 320 IS1=IS1+JS1 IC1=IC1+JS1 J3=IS1+JS1 330 IS0=IS1-JS1 IC0=IC1-JS1 ODD1=SN*Zf(IC1)-CS*Zf(IS1) ODD2=CS*Zf(IC1)+SN*Zf(IS1) Zf(IC1)=Zf(IC0)-ODD1 Zf(IC0)=Zf(IC0)+ODD1 Zf(IS1)=-Zf(IS0)+ODD2 Zf(IS0)=Zf(IS0)+ODD2 IS1=IS1+1 IC1=IC1+1 IF(IS1.NE.J3) GO TO 330 I=I+1 IF(IS1.EQ.J5) GO TO 300 GO TO 310 400 K1=INDX(I) SN=SI(I) I=I+1 CS=SI(I) IS0=IS1 IS1=IS1+JS1 IC0=IC1 IC1=IC1+JS1 ODD1=CS*Zf(IC1)-SN*Zf(IS1) Yf(K1+1)=Zf(IC0)+ODD1 N3MK1=N3-K1 Yf(N3MK1+1)=Zf(IC0)-ODD1 IF(IBC.LT.3) GO TO 410 ODD2=SN*Zf(IC1)+CS*Zf(IS1) N3PK1=N3+K1 Yf(N3PK1+1)=Zf(IS0)+ODD2 N4MK1=N4-K1 Yf(N4MK1+1)=-Zf(IS0)+ODD2 410 IS1=IS1+1 IC1=IC1+1 K1=INDX(I) IS0=IS1 IS1=IS1+JS1 IC0=IC1 IC1=IC1+JS1 ODD1=SN*Zf(IC1)-CS*Zf(IS1) Yf(K1+1)=Zf(IC0)+ODD1 N3MK1=N3-K1 Yf(N3MK1+1)=Zf(IC0)-ODD1 IF(IBC.LT.3) GO TO 420 ODD2=CS*Zf(IC1)+SN*Zf(IS1) N3PK1=N3+K1 Yf(N3PK1+1)=Zf(IS0)+ODD2 N4MK1=N4-K1 Yf(N4MK1+1)=-Zf(IS0)+ODD2 420 IS1=IS1+1 IC1=IC1+1 I=I+1 IF(IS1.NE.J5) GO TO 400 RETURN END c ------------------------------------------------------ SUBROUTINE FOUR67(IBC1,IQ1,ip1,isl1,l11,n21,n71, & si,indx,Zf,Yf) c ------------------------------------------------------ INCLUDE 'PMparameters.h' include 'PMmesh.h' dimension Zf(narr) , Yf(narr) dimension si(nf67) , indx(nf67) integer ibc , ip, isl , l1 , n2 , n3 , n4 , n7 c INCLUDE 'Control.param' IBC=IBC1 IQ=IQ1 ip = ip1 isl = isl1 l1 = l11 n2 = n21 n7 = n71 A5=0.5*SQRT(2.0) N4=2**IQ N3=N4 GO TO (103,103,101,102),IBC 101 Zf(1)=Zf(1)/2.0 Zf(N3+1)=Zf(1) N2=N3 CALL TFOLD(0,1,Zf,n2) N3=N3/2 IQ=IQ-1 GO TO 103 102 N3=N3/2 IQ=IQ-1 103 N5=N3/4 N7=N3/2 N11=3*N7 N31=N3+1 GO TO(300,400,500,600),IBC 300 Zf(N31)=0.0 Zf(1)=0.0 N2=N3 DO 301 I=2,IQ CALL TFOLD(1,1,Zf,n2) 301 N2=N2/2 Yf(N7+1)=Zf(2) JF=N5 ISL=1 DO 302 IP=2,IQ L1=N2+1 CALL ZEERO(1,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) CALL KFOLD(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) I1=3*JF+1 I2=4*JF I3=I1+(N2/2-1)*I2 CALL NEG(I1,I3,I2,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) N2=N2+N2 302 JF=JF/2 RETURN 400 Zf(1)=Zf(1)/2.0 Zf(N31)=Zf(N31)/2.0 N2=N3 DO 401 I=2,IQ CALL TFOLD(0,N31-N2,Zf,n2) 401 N2=N2/2 L1=N31-N2 A=Zf(L1)+Zf(L1+2) Yf(1)=-A-Zf(L1+1) Yf(N31)=-A+Zf(L1+1) Yf(N7+1)=Zf(L1)-Zf(L1+2) DO 402 IP=2,IQ ISL=N31-N2 L1=ISL-N2 idum = isl CALL ZEERO(idum,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) CALL KFOLD(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) 402 N2=N2+N2 CALL NEG(1,N31,2,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) RETURN 500 N2=N3 L2=N4 DO 501 IP=2,IQ CALL TFOLD(1,1,Zf,n2) CALL TFOLD(0,L2-N2+1,Zf,n2) 501 N2=N2/2 L1=L2-N2+1 A=Zf(L1)+Zf(L1+2) Yf(N7+1)=2.0*(-Zf(L1)+Zf(L1+2)) Yf(1)=2.0*(A+Zf(L1+1)) Yf(N31)=2.0*(A-Zf(L1+1)) Yf(N11+1)=2.0*Zf(2) DO 502 IP=2,IQ Zf(N2+1)=2.0*Zf(N2+1) ISL=N2+1 CALL TFOLD1(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) L1=L1-N2 Zf(L1)=-2.0*Zf(L1) CALL KFOLD(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) 502 N2=N2+N2 Yf(1)=Yf(1)*A5 Yf(N31)=Yf(N31)*A5 RETURN 600 Zf(1)=Zf(1)*A5 Zf(N31)=Zf(N31)*A5 N2=N3 DO 601 IP=2,IQ CALL TFOLD(0,N31-N2,Zf,n2) CALL TFOLD(1,N31,Zf,n2) 601 N2=N2/2 L1=N31-N2 A=Zf(L1)+Zf(L1+2) Yf(1)=2.0*(A+Zf(L1+1)) Yf(N31)=2.0*(A-Zf(L1+1)) Yf(N7+1)=2.0*(-Zf(L1)+Zf(L1+2)) Yf(N11+1)=2.0*Zf(N3+2) DO 602 IP=2,IQ ISL=N31+N2 Zf(ISL)=2.0*Zf(ISL) CALL TFOLD1(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) L1=L1-N2 Zf(L1)=-2.0*Zf(L1) CALL KFOLD(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) 602 N2=N2+N2 CALL REVNEG(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) idum = n3 CALL NEG(2,idum,2,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) N2=N4 CALL TFOLD(1,1,Yf,n2) Yf(N4+1)=0.0 RETURN END