C ______________________________ POWER SPECTRUM C 25 april 1990 C R. Kates (Potsdam), A. Klypin (NMSU) C C C Module FFT integer, parameter :: NGRID = 2048 PARAMETER (NARR =2*NGRID+1) ! need in FFT PARAMETER (NF67 =NGRID) end Module FFT C_________________________________________________________ Module Structures use FFT Integer*4, parameter :: Nbins = 50 ! number of bins Real*4, ALLOCATABLE, dimension(:) :: Xp,Yp,Zp,Vc ! coordinates and Vcirc Real*4, ALLOCATABLE, dimension(:) :: VVx,VVy,VVz ! velocity INTEGER*4, DIMENSION(-Nbins:0) :: Ncounts ! counts for CF Real*4, DIMENSION(-Nbins:0) :: Rcounts ! sum(radii) Real*4, DIMENSION(NGRID) :: PkSum =0, PkErr=0. Real*8, DIMENSION(NGRID) :: DPOWER,dk Integer*4, DIMENSION(NGRID) :: iVolume CHARACTER*80 CATALOG,XI Real*4 :: boxsize, ! Box in Mpch units & xmlim, ! low limit on Vcirc & xl,xr,yl,yr,zl,zr ! box limits with periodical conditions Integer*4 :: Ngal, ! number of halos inside boxsize + Ntot, ! total number of halos + ind_redshift ! real (0) or redshift (1,2,3) space REAL :: AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0,factor,Scale COMMON / GRID/ FI(NGRID,NGRID,NGRID) End Module Structures C ______________________________ ______________________________ ! PROGRAM PowerCatalog use Structures Integer*4 :: OMP_GET_THREAD_NUM,OMP_GET_MAX_THREADS NPROCS = 1 NPROCS = OMP_GET_MAX_THREADS() ! comment this line if not using OMP ! Read data and open files Call ReadCatalog(0 ) ! count halos !Ngal =1000000 ALLOCATE(Xp(Ngal),Yp(Ngal),Zp(Ngal),Vc(Ngal)) If(ind_redshift /= 0)ALLOCATE(VVx(Ngal),VVy(Ngal),VVz(Ngal)) Call ReadCatalog(1) ! read halos !CALL GalRandom !CALL GalHomogeneous ! nCatalogs = 40 ! do jj = 1,nCatalogs ! Call ReadCarmen(jj,nCatalogs) ! CALL DENSIT ! Find density ! write (*,*) ' Density is ok' ! CALL DENTES(DELR) ! Statistics ! CALL POWER(NPROCS) ! Power spectrum ! close(17) ! endDo !WRITE (18,'(6X,a,2X,2a)')'k(mpch-1)', !& 'Harmonics ',' Spectrum(h^3) rms Error ' ! DO I=1,NGRID ! IF(iVolume(I).GT.0)THEN ! wk = dk(I)/iVolume(I)*Scale ! Pk = PkSum(I)/nCatalogs ! Pkc = sqrt(max(PkErr(I)/nCatalogs-Pk**2,1.e-20)) ! WRITE (18,110) I,wk,iVolume(I),Pk,Pkc,Pkc/Pk ! ENDIF ! ENDDO !110 FORMAT(1x,I4,F9.4,I10,3X,1p,5G12.5) ! Call ReadIlian CALL DENSIT ! Find density write (*,*) ' Density is ok' CALL DENTES(DELR) ! Statistics CALL POWER(NPROCS) ! Power spectrum WRITE (18,'(6X,a,2X,2a)')'k(mpch-1)', & 'Harmonics ',' Spectrum(h^3) rms Error ' DO I=1,NGRID IF(iVolume(I).GT.0)THEN wk = dk(I)/iVolume(I)*Scale Pk = PkSum(I)/nCatalogs Pkc = sqrt(max(PkErr(I)/nCatalogs-Pk**2,1.e-20)) WRITE (18,110) I,wk,iVolume(I),Pk,Pkc,Pkc/Pk ENDIF ENDDO 110 FORMAT(1x,I4,F9.4,I10,3X,1p,5G12.5) END PROGRAM PowerCatalog C_________________________________________________________ C Compute mean density and rms SUBROUTINE ReadCatalog(iFlag ) use Structures CHARACTER*80 :: tmp_char Character*4 :: ch1,ch2,ch3 If(iFlag == 0)Then open(11,file='direct_cf_mattia.inp') read (11,'(A)') CATALOG read (11,'(A)') XI read (11,*) boxsize read (11,*) xmlim read (11,*) ind_redshift close(11) Ngal = 0 EndIf kin = 0 ! number of halos inside Box kread = 0 write(*,*) ' Index of redshift space= ',ind_redshift write(*,'(2a)') ' Counting halos in ', CATALOG OPEN(19,FILE=CATALOG,STATUS='OLD') ! read header of Catshort file read(19,'(A)') tmp_char if(iFlag == 1)write(17,'(a)')tmp_char read(19,*) ch1,ch2,aex factor = 100.*aex*sqrt(0.27/aex**3+0.73) write(*,*) ' a=',aex,factor do ic1=0,14 read(19,'(A)') tmp_char if(iFlag == 1 .and. ic1 .le.2)write(17,'(a)')tmp_char enddo If(iFlag == 0)Then if(aex >0.0)Then redshift = 1./aex - 1. else redshift = 0. End if iz1 = int(redshift) iz2 = max(min(int(100*(redshift-iz1)),99),0) write(XI,'(2a,i3.3,a,2(i2.2,a),i4.4,a)')Trim(XI),'.z',iz1,'.', & iz2,'.d',ind_redshift,'.',INT(xmlim),'.dat' Open(17,file=TRIM(XI)) End If do ! kh =1,1000 ! ----- reading BDM halo catalogs read(19,*,iostat =ierr)xct,yct,zct,vxt,vyt,vzt,xmass, xtotal, & rvir,vrms,vmax,nhal,cnfw,hasub,index_sub if(ierr /= 0)exit kread = kread +1 if(mod(kread,1000000)==0) & write(*,'(a,i10,3f10.4)') ' Reading line=',kread,xct,yct,zct ! IF ((vmax.GT.xmlim).and.(hasub.EQ.1)) THEN ! IF (vmax.GT.xmlim .and. index_sub == 0) THEN IF (vmax.GT.xmlim) THEN kin = kin + 1 If(iFlag ==1)Then Xp(kin) = xct Yp(kin) = yct Zp(kin) = zct Vc(kin) = vmax SELECT CASE (ind_redshift) CASE (1,2,3) VVx(kin) = vxt VVy(kin) = vyt VVz(kin) = vzt end SELECT End If ENDIF end do If(iFlag == 1)Then write(17,*) ' Total number of halos read = ',kread write(17,*) ' Number of halos selected = ',kin write(17,*) ' limit of halos = ',xmlim write(17,*) ' real(0)/redshift space = ',ind_redshift EndIf write(*,*) ' Flag =',iFlag,kin If(iFlag == 0)Then Ngal = kin write(*,*) NGAL,' objects used to calculate P(k)' Else ! If(kin /= Ngal)Stop 'Stop: error in counting halos' EndIf vol = boxsize**3 dens = NGAL/vol If(iFlag==1)write(17,*) ' Number-density =',dens,' (h/Mpc)^3' close(19) end SUBROUTINE ReadCatalog C----------------------------------------------------- C Compute mean density and rms SUBROUTINE DENTES(DELR) C----------------------------------------------------- use Structures Dimension x(44),amass(44),imass(44) real*8 SUM1,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 SUM1 = 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(+:SUM1,SUM2) DO K=1,NGRID DO J=1,NGRID DO I=1,NGRID SUM1 = SUM1 + FI(I,J,K) SUM2= SUM2+ FI(I,J,K)**2 ! Do l=44,1,-1 ! Dens =FI(I,J,K) +1. ! If(Dens .GT. x(l))Then ! Nn =Nn +1 ! Am =Am +Dens ! imass(l) =imass(l) +1 ! amass(l) =amass(l) +Dens ! GoTo 10 ! EndIf ! EndDo ! 10 Continue ENDDO ENDDO ENDDO Total =(FLOAT(NGRID))**3 DENM = SUM1/Total DELR = DSQRT(SUM2/Total-DENM**2) DensAv=SUM1/float(NGRID)**3 WRITE (*,150) DELR,DENM,DensAv WRITE (17,150) DELR,DENM,DensAv 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) ! &,/' bin Rho Rho Sum(rho_i)/dRho/TotMass Sum',/ ! &'
< left > Normalized density Number of cells ',/ ! &' < bin > Distribution With this density') ! Do k=1,44-1 ! drho =x(k+1)-x(k) ! rhom =x(k)+drho/2. ! Sn =drho*Total ! If(Sn .lt.1.e-7)Sn=1. ! write (*,20) k,rhom,x(k),amass(k)/Sn,imass(k) ! write (17,20) k,rhom,x(k),amass(k)/Sn,imass(k) ! 20 Format(i3,2f8.1,2x,G12.3,12X,i12) ! EndDo ! k =44 ! write (*,20) k,x(k),x(k),amass(k),imass(k) ! write (17,20) k,x(k),x(k),amass(k),imass(k) ! DO K=1,NGRID ! DO J=1,NGRID ! DO I=1,NGRID ! If( FI(I,J,K)> 10.)write(300,'(g12.5,3x,3i6)')FI(I,J,K),I,J,K ! ENDDO ! ENDDO ! ENDDO RETURN END C---------------------------------------------------- C power spectrum P(k) C for given density field FI SUBROUTINE POWER(NPROCS) C---------------------------------------------------- use Structures PARAMETER ( Pi=3.141592653, NPOWER =NGRID) DIMENSION DPW(NPOWER,NPROCS),iVol(NPOWER,NPROCS), & dkk(NPOWER,NPROCS) dimension Zf(narr) , Yf(narr) dimension si(nf67) , indx(nf67) REAL*8 SIGMA,DPW Integer :: OMP_GET_THREAD_NUM si = 0. indx = 0. Box = boxsize 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) ! write(*,'(1p,10g11.4)') si ! DO k=1,NGRID ! DO j=1,NGRID ! DO i=1,NGRID ! FI(i,j,k) = sin(i/200.) ! EndDo ! 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 !write (*,*) ' k=',k 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 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. write (*,*) ' Power .... ' C.... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (iP, k,j,i,MJ,MK,AMP,RK,MP) C$OMP+REDUCTION(+:SIGMA) Do K=1,NGRID iP = 1 iP = OMP_GET_THREAD_NUM()+1 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 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 aNyquist = Pi/(Box/NGRID) Correct = Box**3/Ngal WRITE (*,100) SQRT(SIGMA) WRITE (17,100) SQRT(SIGMA) 100 FORMAT(5x,'RMS DRho/Rho=',F7.3, & /4X,'n',2X,'k(h/Mpc)',1X, & 'Harmonics ','Power Spectrum(h^3) P(k)_noise_corrected') 110 FORMAT(1x,I4,F9.4,I10,3X,1p,5G12.5) DO I=1,NPOWER IF(iVolume(I).GT.0)THEN wk = dk(I)/iVolume(I)*Scale Pk = Box**3*DPOWER(I)/iVolume(I) Pkc = Pk - Correct*(1.-0.66667*(sin(Pi*wk/aNyquist/2.))**2) if(I.lt.32)WRITE (*,110) I,wk,iVolume(I),Pk,Pkc WRITE (17,110) I,wk,iVolume(I),Pk,Pkc PkSum(I) = PkSum(I) + Pk PkErr(I) = PkErr(I) + Pk**2 ENDIF ENDDO END SUBROUTINE POWER C--------------------------------------- C Density Field SUBROUTINE DENSIT C--------------------------------------- use Structures Real*8 :: ss,X,Y,Z,D1,D2,D3,T1,T2,T3 XN =FLOAT(NGRID)+1.-1.E-7 YN =FLOAT(NGRID) C Subtract mean density C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (M1,M2,M3) 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)**3/float(Ngal) PARTW = W write(*,*) ' Particle weight = ',PARTW xmin =1.e+5 xmax =-1.e+5 ymin =1.e+5 ymax =-1.e+5 zmin =1.e+5 zmax =-1.e+5 Do IN =1,Ngal X= Xp(IN) Y= Yp(IN) Z= Zp(IN) SELECT CASE (ind_redshift) CASE (1) X= X + VVx(IN)/factor CASE (2) Y= Y + VVy(IN)/factor CASE (3) Z= Z + VVz(IN)/factor End SELECT X= X*NGRID/boxsize +1.001 Y= Y*NGRID/boxsize +1.001 Z= Z*NGRID/boxsize +1.001 IF(X.ge.NGRID+1.)X=X-NGRID IF(Y.ge.NGRID+1.)Y=Y-NGRID IF(Z.ge.NGRID+1.)Z=Z-NGRID IF(X.lt.1.)X=X+NGRID IF(Y.lt.1.)Y=Y+NGRID IF(Z.lt.1.)Z=Z+NGRID !If(mod(IN,10000)==0) ! write (200+jfile,'(3f11.5,3x,5i12)') ! & X,Y,Z,IN,j,IROW,N,jfile 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,' Irow=',IROW,IN If(J.le.0)write (*,*) ' Y:',X,Y,Z,' Irow=',IROW,IN If(K.le.0)write (*,*) ' Z:',X,Y,Z,' Irow=',IROW,IN If(I.gt.NGRID+1)write (*,*) ' X:',X,Y,Z,' Irow=',IROW,IN If(J.gt.NGRID+1)write (*,*) ' Y:',X,Y,Z,' Irow=',IROW,IN If(K.gt.NGRID+1)write (*,*) ' Z:',X,Y,Z,' Irow=',IROW,IN D1=X-FLOAT(I) D2=Y-FLOAT(J) D3=Z-FLOAT(K) T1=1.-D1 T2=1.-D2 T3=1.-D3 T2W =T2*W D2W =D2*W 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 ss = 0. C Subtract mean density C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (M1,M2,M3) REDUCTION(+:ss) DO M3=1,NGRID DO M2=1,NGRID DO M1=1,NGRID ! FI(M1,M2,M3) = FI(M1,M2,M3)-1. ss = ss + FI(M1,M2,M3) END DO END DO END DO write (*,*) ' Number of particles read = ', Ngal write (*,*) ' Density/cell = ', ss/float(NGRID)**3 write (*,*) ' X:',xmin,xmax write (*,*) ' Y:',ymin,ymax write (*,*) ' Z:',zmin,zmax END SUBROUTINE DENSIT C-------------------------------------- Fourier Transform SUBROUTINE SETF67(IBC1,IQ1,ibc,ip,isl,l1,n2,n3,n4,n7, & si,indx,Zf,Yf) c ----------------------------------------------------- use FFT 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 SUBROUTINE SETF67 C----------------------------------------- SUBROUTINE TFOLD(IS,L,ZZZ,n2) use FFT integer n2 DIMENSION ZZZ(NARR) IH2=N2/2-1 DO I=IS,IH2 I1=I+L I2=N2-I+L A=ZZZ(I1) ZZZ(I1)=A-ZZZ(I2) ZZZ(I2)=A+ZZZ(I2) end DO RETURN END SUBROUTINE TFOLD C------------------------------------ SUBROUTINE NEG(I1,I3,I2,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) use FFT dimension Zf(narr) , Yf(narr) dimension si(nf67) , indx(nf67) integer ibc , ip, isl , l1 , n2 , n3 , n4 , n7 DO K=I1,I3,I2 Yf(K)=-Yf(K) end DO END SUBROUTINE NEG C------------------------------------- SUBROUTINE REVNEG(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) use FFT 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) use FFT dimension Zf(narr) , Yf(narr) dimension si(nf67) , indx(nf67) integer ibc , ip, isl , l1 , n2 , n3 , n4 , n7 DO I=1,N2 LI=L+I Zf(LI-1)=0.0 end DO END SUBROUTINE ZEERO C------------------------------------------ SUBROUTINE TFOLD1(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) use FFT dimension Zf(narr) , Yf(narr) dimension si(nf67) , indx(nf67) integer ibc , ip, isl , l1 , n2 , n3 , n4 , n7 II=N2-1 DO I=1,II I1=ISL+I I2=L1-I A=Zf(I1) Zf(I1)=A+Zf(I2) Zf(I2)=A-Zf(I2) end DO END SUBROUTINE TFOLD1 C------------------------------------- SUBROUTINE KFOLD(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf) use FFT 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 END SUBROUTINE KFOLD c ------------------------------------------------------ SUBROUTINE FOUR67(IBC1,IQ1,ip1,isl1,l11,n21,n71, & si,indx,Zf,Yf) c ------------------------------------------------------ use FFT 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 SUBROUTINE FOUR67 !---------------------------------------------------------------------- Subroutine GalHomogeneous use Structures ! ! Nn = 256 ii = 0 Do k=1,Nn z = (k-1)*boxsize/Nn Do j=1,Nn y = (j-1)*boxsize/Nn Do i=1,Nn x = (i-1)*boxsize/Nn ii = ii+1 Xp(ii) = x Yp(ii) = y Zp(ii) = z End Do end Do end Do Ngal =Nn**3 end Subroutine GalHomogeneous !---------------------------------------------------------------------- Subroutine GalRandom use Structures ! generating random distribution of NGAL halos instead of using real halos ! Nseed = 163863235 !Nseed = 9121071 !Nseed = 121071 Do i =1,Ngal Xp(i) = randd(Nseed)*boxsize End Do Do i =1,Ngal Yp(i) = randd(Nseed)*boxsize End Do Do i =1,Ngal Zp(i) = randd(Nseed)*boxsize Vc(i) = 0. EndDo ! Do i=1,Ngal ! read(20,*) Xp(i),Yp(i),Zp(i) ! EndDo end Subroutine GalRandom ! ------------------ FUNCTION RANDd(M) ! ------------------ ! random generator ! ------------------ DATA LC,AM,KI,K1,K2,K3,K4,L1,L2,L3,L4 + /453815927,2147483648.,2147483647,536870912,131072,256,16777216, + 4,16384,8388608,128/ ML=M/K1*K1 M1=(M-ML)*L1 ML=M/K2*K2 M2=(M-ML)*L2 ML=M/K3*K3 M3=(M-ML)*L3 ML=M/K4*K4 M4=(M-ML)*L4 M5=KI-M IF(M1.GE.M5)M1=M1-KI-1 ML=M+M1 M5=KI-ML IF(M2.GE.M5)M2=M2-KI-1 ML=ML+M2 M5=KI-ML IF(M3.GE.M5)M3=M3-KI-1 ML=ML+M3 M5=KI-ML IF(M4.GE.M5)M4=M4-KI-1 ML=ML+M4 M5=KI-ML IF(LC.GE.M5)ML=ML-KI-1 M=ML+LC RANDd=M/AM END FUNCTION RANDd C_________________________________________________________ C Compute mean density and rms SUBROUTINE ReadCarmen(iCatalog,nCatalogs ) use Structures CHARACTER*120 :: Line Character*4 :: ch1,ch2,ch3 write(Line,'(a,i2.2,a)') & 'DATA/mock_gammatweak_boss_cmass_Carmen_20', & iCatalog,'_z0p520_fof_b0p156.ff' If(iCatalog /= 1) DEALLOCATE(xp,yp,zp,vvx,vvy,vvz) write(*,*) ' Read CARMEN files = ',iCatalog redshift = 0.52 ind_redshift = 2 xmlim = iCatalog Om0 = 0.25 Open(1,file=TRIM(Line),form='unformatted') read (1) n,np,iskip,itrans,jseed write(*,*) ' Number of halos = ',Np write(*,*) ' Number of particles = ',n Allocate(xp(np),yp(np),zp(np)) Allocate(vvx(np),vvy(np),vvz(np)) read(1) rcube,h,xindx,rsm1,rth1,thnorm,pnorm,gnu,eta read(1) znow write(*,*) ' Redshift = ',znow write(*,*) ' Box = ',rcube write(*,*) ' Hubble = ',h read(1) (xp(i),i=1,np) read(1) (yp(i),i=1,np) read(1) (zp(i),i=1,np) read(1) (vvx(i),i=1,np) read(1) (vvy(i),i=1,np) read(1) (vvz(i),i=1,np) boxsize = rcube Ngal = np AEXPN = 1/redshift-1. factor = 100.*AEXPN*sqrt(Om0/AEXPN**3+1.-Om0) iz1 = int(redshift) iz2 = max(min(int(100*(redshift-iz1)),99),0) write(Line,'(2a,i1.1,a,2(i2.2,a),i4.4,a)')'Carmen','.z',iz1,'.', & iz2,'.d',ind_redshift,'.',INT(xmlim),'.dat' Open(17,file=TRIM(Line)) write(17,*)' Carmen file =',iCatalog write(17,*)' Box =',boxsize write(17,*)' Nhalos =',Ngal write(17,*)' Real/Redshift =',ind_redshift write(17,*)' Number-density=',Ngal/boxsize**3 If(iCatalog == 1)Then Open(18,file='CarmenAllSpectraY.dat') write(18,*)' Carmen files =',nCatalogs write(18,*)' Box =',boxsize write(18,*)' Nhalos =',Ngal write(18,*)' Real/Redshift =',ind_redshift write(18,*)' Number-density=',Ngal/boxsize**3 EndIf xmin = MINVAL(xp) ymin = MINVAL(yp) zmin = MINVAL(zp) write(*,'(a,1p,3g12.4)') ' Minimum of coords =',xmin,ymin,zmin xmax = MAXVAL(xp) ymax = MAXVAL(yp) zmax = MAXVAL(zp) write(*,'(a,1p,3g12.4)') ' Maximum of coords =',xmax,ymax,zmax xmin = MINVAL(vvx) ymin = MINVAL(vvy) zmin = MINVAL(vvz) write(*,'(a,1p,3g12.4)') ' Minimum of veloc =',xmin,ymin,zmin xmax = MAXVAL(vvx) ymax = MAXVAL(vvy) zmax = MAXVAL(vvz) write(*,'(a,1p,3g12.4)') ' Maximum of veloc =',xmax,ymax,zmax close (1) End SUBROUTINE ReadCarmen C_________________________________________________________ C Compute mean density and rms SUBROUTINE ReadIlian use Structures CHARACTER*120 :: Line Character*4 :: ch1,ch2,ch3 write(Line,'(a,i2.2,a)') & 'DATA/1Gpc_3546_z0_vmax250_halos' write(*,*) ' Read Ilian files ' redshift = 0. ind_redshift = 0 xmlim = 250. Om0 = 0.27 np = 812338 boxsize = 1000. Open(1,file=TRIM(Line)) read (1,*) Line write(*,*) ' Number of halos = ',Np Allocate(xp(np),yp(np),zp(np)) Allocate(vvx(np),vvy(np),vvz(np)) write(*,*) ' Redshift = ',redshift write(*,*) ' Box = ',boxsize read(1,'(a)')Line !read(1,'(a)')Line write(*,*)Line Do i=1,np read(1,*) xp(i),yp(i),zp(i) !write(*,*)i,xp(i) EndDo Ngal = np AEXPN = 1/redshift-1. factor = 100.*AEXPN*sqrt(Om0/AEXPN**3+1.-Om0) iz1 = int(redshift) iz2 = max(min(int(100*(redshift-iz1)),99),0) write(Line,'(2a,i1.1,a,2(i2.2,a),i4.4,a)')'Ilian','.z',iz1,'.', & iz2,'.d',ind_redshift,'.',INT(xmlim),'.dat' Open(17,file=TRIM(Line)) write(17,*)' Ilian file =',0 write(17,*)' Box =',boxsize write(17,*)' Nhalos =',Ngal write(17,*)' Real/Redshift =',ind_redshift write(17,*)' Number-density=',Ngal/boxsize**3 xmin = MINVAL(xp) ymin = MINVAL(yp) zmin = MINVAL(zp) write(*,'(a,1p,3g12.4)') ' Minimum of coords =',xmin,ymin,zmin xmax = MAXVAL(xp) ymax = MAXVAL(yp) zmax = MAXVAL(zp) write(*,'(a,1p,3g12.4)') ' Maximum of coords =',xmax,ymax,zmax close (1) End SUBROUTINE ReadIlian