C--------------------------------------------------- C Spectrum of perturbations: P(k), sigma_8,V(50h-1Mpc) C Anatoly Klypin December 1997 C Holtzman appoximation C or BBKS-style + Hu&Sugiyama approx C Needs cdm.fit file -- parameters of models and fitting C----------------------------------------------------- C Om = Omega_0 Omc = Omega_cdm_present C Omb = Omega_baryon Omnu =Omega_neutrino_present C Ocurv = Omega_curvature_present C hsmall = H_0/100km/s/Mpc = the Hubble constant C ns = slope of the power spectum (ns=1 for Harr-Zeld) C AEXPN = expansion parameter = 1/(1+z), z =redshift C All calculations are done in real Mpc, C Only final outputs are scaled to h^{-1} C----------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NtabM = 100000) PARAMETER (Nm = 100) PARAMETER (Na = 10) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup,GrowthDen_0 COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),StepK,alog0,iFlag,Ntab CHARACTER Name*16,Answer*1,Header*45,StartFile*16 REAL*8 INTG,ns,INPUT,P2,P,Ptophat,Pgauss,Vtophat DIMENSION aNdens(Nm,Na),aMasses(Nm),aRedshifts(Na) EXTERNAL INTG,P2,P,Ptophat,Pgauss,Vtophat DATA StartFile/'InStart.dat'/ WRITE (*,*) '===Enter File Name for output file [lcdm.dat]: ' READ (*,'(a)',IOSTAT=iStat) Name if(iStat.ne.0)Name='lcdm.dat' if(TRIM(Name)=='')Name='lcdm.dat' OPEN(1,file=Name) OPEN(10,file='number'//Name) OPEN(2,file='cdm.fit') ! parameters of models C Constants PI =3.1415926535 AEXPN =1. WRITE (*,'(A,$)') ' Enter sigma_8 and slope n: ' READ (*,*) Sigma8,ns CALL MODEL(Line) H =100.*hsmall ! Hubble constant DO ia=0,50 a =0.1*10**(ia/50.) z =1./a-1. CALL AGE(t0,GrowthDen,GrowthVel,a) write (*,20) t0,z,a,GrowthDen,GrowthVel write (1,20) t0,z,a,GrowthDen,GrowthVel write (10,20) t0,z,a,GrowthDen,GrowthVel 20 format(' Age=',f8.3,' z=',f6.2,' a=',f6.3, . ' GrowthRateDen=',f7.4, . ' GrowthRateVelocity=',f7.4) ENDDO AEXPN =1. a =AEXPN z =1./a-1. write (*,*) ' a=',a,' Aexpn=',AEXPN CALL AGE(t0,GrowthDen_0,GrowthVel_0,a) write (*,20) t0,z,a,GrowthDen_0,GrowthVel_0 write (1,20) t0,z,a,GrowthDen_0,GrowthVel_0 write (10,20) t0,z,a,GrowthDen_0,GrowthVel_0 rf =8./hsmall ! top-hat for sigma_8 Sn = (Sigma8)**2 / & INTG(Ptophat,1.d-5,5.d+1) ! normalization of P(k) c bias parameter Sigma8 = sqrt(Sn*( & INTG(Ptophat,1.d-5,5.d+1))) write (*,10) hsmall,Sigma8,ns ! check if all is self-consistent write (1,10) hsmall,Sigma8,ns write (10,10) hsmall,Sigma8,ns 10 Format(' Hubble=',f7.3,' Sigma_8 =',G11.3, & ' Slope n=',f7.2) c Evolution of Sigma(mass,a),Fraction_collapsed(a) write(1,15)' Z a' 15 format(a,T17,3(3x,'Mass',7x,'Sigma', & ' FracMass')) Do i=0,30 z =i/3. AEXPN = 1./(1.+z) aMass1 = 1.e11 ss1 = SigMass(aMass1) fr_mass1 = 0.5*derfc(1.68d0/sqrt(2.)/ss1) aMass2 = 1.e12 ss2 = SigMass(aMass2) fr_mass2 = 0.5*derfc(1.68d0/sqrt(2.)/ss2) aMass3 = 1.e13 ss3 = SigMass(aMass3) fr_mass3 = 0.5*derfc(1.68d0/sqrt(2.)/ss3) write (1,'(2f8.3,3(g11.3,f9.3,g11.3))')z,AEXPN, & aMass1,ss1,fr_mass1, & aMass2,ss2,fr_mass2,aMass3,ss3,fr_mass3 EndDo !AEXPN =1. !write(*,*) ' Mass sigma(Mass)' !Do i=0,150 ! aM = 1.e8*10.**(i/15.) ! ss1 = SigMass(aM) ! write (*,'(1P,2g11.4)')aM*hsmall,ss1 !EndDo C Mass function at different redshifts Do i=1,Nm aMasses(i) = 1.d10*10.**(i/float(Nm)*5.) EndDo Do i=1,Na aRedshifts(i) = (i/float(Na)*10.) EndDo aRedshifts(1) =0. Do j=1,Na z = aRedshifts(j) AEXPN = 1./(1.+z) Do i=1,Nm aMass = aMasses(i) aNdens(i,j) = aNumber(aMass) EndDo EndDo write (10,'(2a)')' M(sun/h)',' Redshifts' write (10,'(11x,16g11.3)') aRedshifts Do i=1,Nm write (10,'(g11.3,16g12.4)') & aMasses(i)*hsmall,(aNdens(i,j)/hsmall**3,j=1,Na) EndDo WRITE (*,'(A,$)') & 'Enter expansion parameter for mass function:' READ (*,*) AEXPN z = 1./AEXPN -1. ! mass function for a specific moment CALL AGE(t0,GrowthDen,GrowthVel,AEXPN) gr0 = GrowthDen/GrowthDen_0 Do i=1,Nm aMass = aMasses(i) aNdens(i,1) = aNumber(aMass) EndDo write (10,'(2a,f8.4,a,f8.5)')' M(sun/h) Mdn/dM (Mpch**-3)', & ' sigma Redshift=', & z,' Growth Rate=',gr0 Do i=1,Nm write (10,'(g11.3,16g12.4)') & aMasses(i)*hsmall,(aNdens(i,1)/hsmall**3), & sigMass(aMasses(i)) EndDo !z =0. !AEXPN =1. !rf =50./hsmall ! top-hat for bulk velocity !Veloc = GrowthVel*H* !!& sqrt(Sn*(INTG(Vtophat,1.d-5,1.d+0))) !write (*,30) Veloc,rf*hsmall !write (1,30) Veloc,rf*hsmall 30 Format(' Bulk Velocity =',F8.2, . 'km/s for radius top-hat=',f8.1,'h-1Mpc') !write (1,*) ' k/h Pk*h^3 Power Spectrum at z=',1/AEXPN-1. !DO i=1,110 ! w =1.e-4*10.**(i/20.) ! Pk0 =P(w)*Sn*hsmall**3*(2.*pi**2) ! write (1,14) w/hsmall,Pk0 !ENDDO 14 Format(2G11.4,5G12.4) C..........................................................Power spectrum and C Amplitude of fluctuations in a Box at redshift z WRITE (*,'(A,$)') & 'Enter Box(h^-1Mpc), Nparticles(1D) and redshift:' READ (*,*) Box, NROW, z Box =Box/hsmall ! scale it to real Mpc Uklow =2.*pi/Box ! frequency range for integrals Ukup =2.*pi/Box*(NROW/2.) AEXPN =1./(1.+z) ! expansion parameter a =AEXPN CALL AGE(t0,GrowthDen,GrowthVel,a) sigma =(GrowthDen/GrowthDen_0)* & sqrt(Sn*INTG(P2,Uklow/sqrt(2.),Ukup)) WRITE (*,40) z,sigma,NROW,Box*hsmall 40 format(' z=',f8.3,' delta\rho/rho in box=',f9.5,/ . 5X,'Particles=',i4,' Box=',f8.2) write (1,*) ' k/h Pk*h^3 Power Spectrum at z=',1/AEXPN-1. DO i=1,130 w =1.e-4*10.**(i/20.) Pk0 =(GrowthDen/GrowthDen_0)**2*P(w)*Sn*hsmall**3*(2.*pi**2) write (1,14) w/hsmall,Pk0 ENDDO C............................................................................................ WRITE (*,'(A,$)') ' Do you need a file to run PMstart (Yes/No) ' READ (*,'(a)') Answer IF(Answer.eq.'Y'.or.Answer.eq.'y')Then ! prepare input file OPEN(30,file=StartFile) ! parameters of a run HEADER='N=128x256 L=20h-1CDMsig=0.239 z=30-----------' write (*,*) '------ Enter Header up to 45 characters' write (*,*) ' Example:' write (*,'(A)') HEADER read (*,'(A)') HEADER write (30,'(A)') HEADER write (30,*) AEXPN,' Expansion Parameter' ASTEP =INPUT(' Step of the expansion parameter =') write (30,*) ASTEP,' Step in dAEXPN ' write (30,*) sigma,' DRho/rho in box ' write (30,*) ns, ' Slope of P(k) ' write (30,*) Box*hsmall,' Box in Mpc/h ' Nseed =INPUT(' Enter random seed integer 1-2^31 =') write (30,*) Nseed,' Random seed ' write (30,*) Line,' Line number in cdm.fit' Write (*,'(A,$)') ' Enter Min and Max number of mass species =' read (*,*) Nmin_lev,Nmax_lev write (30,*) Nmin_lev,Nmax_lev,' Min and Max number of species' write (*,*) ' Results were written to file: ',StartFile CLOSE (30) ENDIF STOP END c----------------------------------- Pcold(k)*k^2 REAL*8 FUNCTION P2(WK) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup,GrowthDen_0 Real*8 ns P2=WK**2*P(WK) RETURN END C--------------------------------------------- C Power spectrum C wk = k= in real Mpc REAL*8 FUNCTION P(wk) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NtabM = 100000) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / FACT/ qqscaleb COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),StepK,alog0,iFlag,Ntab Real*8 ns If(iFlag.gt.0)Then ! use approximations If(Par(6).ne.0.)Then ! Holtzman approx sk= sqrt(wk) P = wk**ns / . (1.+sk*(Par(2) . +sk*(Par(3) . +sk*(Par(4) . +sk*(Par(5) )))) )**(2.*Par(6)) Else ! BBKS + Sugiyama approx c Gammaeff =Om*hsmall/exp(Omb*(1.+sqrt(hsmall/0.5)/Om)) c Q = wk/hsmall/Gammaeff Q = wk*qqscaleb P = wk**ns*( LOG(1.+Par(1)*Q)/(Par(1)*Q) )**2/ . sqrt(1.+Q*(Par(2)+ . Q*(Par(3)**2+ . Q*(Par(4)**3+ . Q*(Par(5)**4) )))) EndIf Else ! use table for P(k) P = Ppk(wk) EndIf RETURN END C------------------------------------------------- C Age of the Universe: t0 (z=0) C Expansion parameter: a =1/(1+z) C Growth_Rate_Density at a: GrowthDen C GrowthDen =\delta\rho/\rho C normalized to GrowthDen =a for Omega=1 C GrowthDen < 1 at a=1 for Lcdm or OCDM C Growth_Rate_velocity at a: GrowthVel C GrowthVel = V_peculiar= (a/delta)(d delta/d a) SUBROUTINE AGE(t0,GrowthDen,GrowthVel,a) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (zero =1.d-12) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall Common /OmegRat/ OmLOm0,OmcOm0 Real*8 INTG,ns External Hnorm,Hage,Hgrow Oml =Max(Abs(1.-Om-Ocurv),zero) OmLOm0 =Oml/Om OmcOm0 =Ocurv/Om t0 =9.766/hsmall/sqrt(Om)*INTG(Hage,zero,a) z =1./a-1. Hubble = hsmall*sqrt(Om/a**3)*Hnorm(a) ww = INTG(Hgrow,zero,a) GrowthVel = -(1.5+OmcOm0*a)/Hnorm(a)**2+ & sqrt(a**5)/Hnorm(a)**3/ww GrowthDen =2.5*Hnorm(a)/sqrt(a**3)*ww Return End c----------------------------------- P(k)*k^2*gauss(rf) REAL*8 FUNCTION Pgauss(WK) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup,GrowthDen_0 Real*8 ns Pgauss=WK**2*P(WK)*EXP(-(rf*WK)**2) RETURN END c----------------------------------- P(k)*k^4*gauss(rf) REAL*8 FUNCTION P4gauss(WK) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup,GrowthDen_0 Real*8 ns P4gauss=WK**4*P(WK)*EXP(-(rf*WK)**2) RETURN END C-------------------------------------- P*k^2*Top-Hat Filter REAL*8 FUNCTION Ptophat(wk) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup,GrowthDen_0 Real*8 ns IF (wk.lt.1.d-4) THEN Ptophat =wk**ns*wk**2 ELSE X =wk*rf TOPHAT =( (SIN(X)-x*COS(X))*3./X**3 )**2 Ptophat=P(wk)*wk**2*TOPHAT ENDIF RETURN END C-------------------------------------- Number-density of halos C Gives M*dN/dM in units Mpc**-3 C aMass is in real Msun units c Calls function Fn(sigma) REAL*8 FUNCTION aNumber(aMass) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (pi = 3.1415926535d0) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup,GrowthDen_0 REAL*8 INTG,ns,INPUT,Ptophat,kmax EXTERNAL INTG,Ptophat kmax = 200. !2.*pi/250.*1024.*0.7 ! rf is in comoving real Mpc rf = (aMass/(1.154e12*Om*hsmall**2))**0.33333333 a = AEXPN CALL AGE(t0,GrowthDen,GrowthVel,a) Sig =(GrowthDen/GrowthDen_0)* & sqrt(Sn*INTG(Ptophat,1.d-4,kmax)) aM1 = aMass*1.03 rf = (aM1/(1.154e12*Om*hsmall**2))**0.33333333 Sig1=(GrowthDen/GrowthDen_0)* & sqrt(Sn*INTG(Ptophat,1.d-4,kmax)) aM0 = aMass/1.03 rf = (aM0/(1.154e12*Om*hsmall**2))**0.33333333 Sig0=(GrowthDen/GrowthDen_0)* & sqrt(Sn*INTG(Ptophat,1.d-4,kmax)) derivative = (Sig0-Sig1)/Sig/(aM1-aM0) factor = Fst(Sig) aNumber = 2.75e11*(Om*hsmall**2)* & derivative * factor c write (*,'(9g12.4)')aMass,Sig,Sig0,Sig1, c & derivative,factor,aNumber if(aNumber.lt.0..or.derivative.lt.0.) & write (*,'(a,12g12.4)')' Error PS:',aMass,Sig,Sig1,Sig0, & derivative,factor,aNumber,aM0,aM1 RETURN END C-------------------------------------- Auxiliary for aNumber() C PS approximation C REAL*8 FUNCTION Fps(sigma) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (pi = 3.1415926535d0) PARAMETER (Deltac = 1.686) x = Deltac/sigma Fps =sqrt(2./pi)*x*exp(-0.5*x**2) RETURN END C-------------------------------------- Auxiliary for aNumber() C ST approximation C REAL*8 FUNCTION Fst(sigma) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (pi = 3.1415926535d0) PARAMETER (Deltac = 1.686) PARAMETER (A = 0.3222) PARAMETER (As = 0.707) PARAMETER (p = 0.3) x = Deltac/sigma Fst =A*sqrt(2.*As/pi)*(1.d0+(1./(As*x**2))**p)*x*exp(-0.5*As*x**2) RETURN END C-------------------------------------- Auxiliary for aNumber() C Tinker approximation C REAL*8 FUNCTION Ftn(sigma) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (A = 0.186) PARAMETER (b = 2.57) PARAMETER (c = 1.19) PARAMETER (ac = 1.47) Ftn =A*(1.d0+(b/sigma)**ac)*exp(-c/sigma**2) RETURN END C-------------------------------------- rms fluctuations C for mass aMass C aMass is in real Msun units REAL*8 FUNCTION SigMass(aMass) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (pi = 3.1415926535d0) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup,GrowthDen_0 REAL*8 INTG,ns,INPUT,Ptophat EXTERNAL INTG,Ptophat rf= (aMass/(1.154e12*Om*hsmall**2))**0.33333333 a =AEXPN CALL AGE(t0,GrowthDen,GrowthVel,a) SigMass =(GrowthDen/GrowthDen_0)* & sqrt(Sn*INTG(Ptophat,1.d-5,2.d+2)) RETURN END C-------------------------------------- P*Top-Hat Filter REAL*8 FUNCTION Vtophat(wk) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup,GrowthDen_0 Real*8 ns IF (wk.lt.1.d-4) THEN Vtophat =wk**ns ELSE X=wk*rf TOPHAT=( (SIN(X)-x*COS(X))*3./X**3 )**2 Vtophat=P(wk)*TOPHAT ENDIF RETURN END C-------------------------------------- Linear Pair-wise V REAL*8 FUNCTION Vpairwise(wk) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup,GrowthDen_0 Real*8 ns X =wk*rf Vpairwise=P(wk)*(1.-sin(X)/X) RETURN END C________________________________________Read parameters of the model C Om = Omega_0 C Omb = Omega_baryon C Omnu = Omega_neutrino C hsmall = hubble constant C Par = fitting parameters C Par(6) = 0 --- bbks-style+Hu&Sugiyama C Par(6) ne 0 --- holtzman-style SUBROUTINE MODEL(Line) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NtabM = 100000) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / FACT/ qqscaleb COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),StepK,alog0,iFlag,Ntab INCLUDE 'PMinitialp.h' Real*8 ns,INPUT Character Header*79 Line1 = INPUT (' Enter Line Number in cdm.fit =') Line = Line1 If (Line1.gt.0) Then ! use analytical approximation iFlag = Line1 OPEN (2, file = 'cdm.fit') Read (2, '(A)') Header1 If (Line1.gt.2) Then Do i = 1, Line1 - 2 Read (2, * ) a EndDo EndIf Read (2, * ) Om, Omb, Omc, Omnu, hsmall, Par CLOSE (2) ! = T_cmb/2.7 theta = 2.726 / 2.7 ! Hu&Sugiyama fitting Ob0 = Omb / Om Omh2 = Om * hsmall**2 a1 = (46.9 * Omh2) **0.670 * (1. + (32.1 * Omh2)**( - 0.532) ) a2 = (12.0 * Omh2) **0.424 * (1. + (45. * Omh2)** ( - 0.582) ) alph = 1. / a1**Ob0 / a2** (Ob0**3) qqscaleb = theta**2 / (1. - Ob0) **0.60 / sqrt (alph) / Omh2 Write (*,20) Om, Omb, Omc, Omnu, hsmall, Par ! Write (16,20)Om,Omb,Omc,Omnu,hsmall,Par 20 Format (' Model: Om_0=', F6.3, ' Om_baryon=', F6.3, ' Om_cold=', & F6.3, ' O_nu=', F6.3, ' hsmall=', F6.2, / 8x, 'Parameters=', & (6G12.3) ) Else ! read spectrum from table iFlag0 = -Line If(iFlag0 == 1)Then Call ReadPkTable Om0 = Omc Else hsmall = ParSet(1,iFlag0) Omb = ParSet(2,iFlag0) Om = ParSet(3,iFlag0) Omc = Om Open(50,file=TableNames(iFlag0)) CALL ReadTable Close (50) EndIf Write ( *, 22) Om, Omb, hsmall,TableNames(iFlag0) Write (1, 22) Om, Omb, hsmall,TableNames(iFlag0) 22 Format (' Model: Om_0=', F6.3, ' Om_baryon=', F6.3, ' hsmall=', & F6.2, ' Using table for P(k)=>',a) EndIf Return END SUBROUTINE MODEL !--------------------------------------- SUBROUTINE ReadPkTable ! interpolate table with p(k) !--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NtabM = 100000) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),StepK,alog0,iFlag,Ntab Real*8 ns Character*80 pkTable,Line Logical :: exst Inquire(file='PkTable.dat',exist=exst) If(exst)Then Open(50,file='PkTable.dat') Else write(*,*) ' Table with the power', & ' spectrum PkTable.dat is not found' stop End If Call ReadHeader ! read header, get Om0, sigma8 Ntab = 0 12 READ (50, *, end = 32, err = 32) xx, pp Ntab = Ntab + 1 xkt (Ntab) = xx*hsmall Pkt (Ntab) = pp ! Pk GOTO 12 32 Write (*,*) Write (*,*) ' Read ', Ntab, ' lines from P(k) table ' If (Ntab.le.1) stop 'wrong table for p(k)' StepK = log10 (xkt (Ntab) / xkt (1) )/(Ntab-1) alog0 = log10 (xkt (1) ) ! test that spacing is the same Do k = 2, Ntab - 1 ss = log10 (xkt (k + 1) / xkt (k) ) If (abs (ss / StepK - 1.) .gt.2.e-2) Then Write (*,*) ' error in K spacing. k=', k, xkt (k + 1),xkt (k) STOP EndIf EndDo END !------------------- Subroutine ReadHeader IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NtabM = 100000) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),StepK,alog0,iFlag,Ntab Real*8 ns Logical :: exst Omb = ParseLine() Omc = ParseLine() OmL = ParseLine() Om0 = ParseLine() sigma8 = ParseLine() hsmall =0.678 Ocurv =0. Om =Om0 ns =0.96 write(*,'(3(a,ES12.3))') & ' Model: Ombar =',Omb,' Om_matter =',Om0,' Sigma8 =',sigma8 end Subroutine ReadHeader !------------------- ! read line from Header of Pk ! Function ParseLine() IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NtabM = 100000) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall Real*8 ns Character :: Line*120,Line2*120,Line3(120) Read(50,'(a)')Line Ieq =INDEX(Line,'=',BACK=.TRUE.) !write(*,*) ' Ieq =',Ieq backspace (50) write(Line2,'(a1,i2,a)') '(',Ieq,'a1,g12.5)' !write(*,'(a)') Line2 Read(50,Line2)(Line3(i),i=1,Ieq),dummy ParseLine = dummy !write(*,'(a,ES12.3)') ' Result =',ParseLine end Function ParseLine C-------------------------------------------- Simpson integration REAL*8 FUNCTION INTG(FUNC,A,B) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (EPS=3.0d-5, JMAX=22) EXTERNAL FUNC OST=-1.D30 OS= -1.D30 ST =0. DO 11 J=1,JMAX CALL TRAPZD(FUNC,A,B,ST,J) INTG=(4.0d0*ST-OST)/3.0d0 IF (ABS(INTG-OS).Le.EPS*ABS(OS)) RETURN OS=INTG OST=ST 11 CONTINUE WRITE (*,*)'Integration did not converge' RETURN END C---------------------------------------------- SUBROUTINE TRAPZD(FUNCC,A,B,S,N) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) SAVE IT EXTERNAL FUNCC IF (N.EQ.1) THEN S=0.5d0*(B-A)*(FUNCC(A)+FUNCC(B)) IT=1 ELSE TNM=IT DEL=(B-A)/TNM X=A+0.5D0*DEL SUM=0.0D0 DO 11 J=1,IT SUM=SUM+FUNCC(X) X=X+DEL 11 CONTINUE S=0.5D0*(S+(B-A)*SUM/TNM) IT=2*IT ENDIF RETURN END C-------------------------------------------- Simpson integration REAL*8 FUNCTION INTG1(FUNC,A,B) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (EPS=1.0d-3, JMAX=22) EXTERNAL FUNC OST=-1.d30 OS= -1.d30 ST =0. DO 11 J=1,JMAX CALL TRAPZD1(FUNC,A,B,ST,J) INTG1=(4.0D0*ST-OST)/3.0D0 IF (ABS(INTG1-OS).Le.EPS*ABS(OS)) RETURN OS=INTG1 OST=ST 11 CONTINUE c WRITE (1,*)'Integration did not converge' RETURN END C---------------------------------------------- SUBROUTINE TRAPZD1(FUNCC,A,B,S,N) C--------------------------------------- IMPLICIT Real*8 (A-H,O-Z) SAVE IT EXTERNAL FUNCC IF (N.EQ.1) THEN S=0.5D0*(B-A)*(FUNCC(A)+FUNCC(B)) IT=1 ELSE TNM=IT DEL=(B-A)/TNM X=A+0.5D0*DEL SUM=0.0D0 DO 11 J=1,IT SUM=SUM+FUNCC(X) X=X+DEL 11 CONTINUE S=0.5D0*(S+(B-A)*SUM/TNM) IT=2*IT ENDIF RETURN END C---------------------------------- Read in variables REAL*8 FUNCTION INPUT(text) C----------------------------------------- Character text*(*) write (*,'(A,$)')text read (*,*) x INPUT =x Return End C----------------------------------------- REAL*8 FUNCTION sinhh(x) C----------------------------------------- IMPLICIT REAL*8 (A-H,P-Z) sinhh =0.5*(exp(MIN(x,33.d+0))-exp(-MIN(x,33.d+0))) RETURN END C--------------------------------------- REAL*8 FUNCTION Hh(x) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) Hh =(sqrt(x/(1.d+0 +x**3)))**3 Return End C--------------------------------------- REAL*8 FUNCTION Hnorm(x) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) Common /OmegRat/ OmLOm0,OmcOm0 Hnorm =sqrt(1.d+0 & +OmLOm0*x**3 & +OmcOm0*x ) Return End C--------------------------------------- REAL*8 FUNCTION Hage(x) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) Hage =sqrt(x)/Hnorm(x) Return End C--------------------------------------- REAL*8 FUNCTION Hgrow(x) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) Hgrow =(sqrt(x)/Hnorm(x))**3 Return End C--------------------------------------- REAL*8 FUNCTION Ppk(x) c interpolate table with p(k) c x is in real 1/Mpc C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NtabM = 100000) COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),StepK,alog0,iFlag,Ntab If(x.ge.xkt(Ntab))Then ! slope is ns =-3 Ppk =Pkt(Ntab)/(x/xkt(Ntab))**3 Return EndIf If(x.lt.xkt(1))Then ! slope is ns=1 Ppk =Pkt(1)*(x/xkt(1)) Return EndIf ind = INT((log10(x)-alog0)/StepK) +1 dk = xkt(ind+1)-xkt(ind) Ppk = (Pkt(ind)*(xkt(ind+1)-x)+Pkt(ind+1)*(x-xkt(ind)))/dk Return End C--------------------------------------- SUBROUTINE ReadTable c interpolate table with p(k) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NtabM = 100000) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),StepK,alog0,iFlag,Ntab REAL*8 ns Character*80 pkTable If (iFlag.eq. -1) Then ! W.Hu table twop = 2. * 3.14145926**2 Ntab = 0 10 READ (50, *, end = 30, err = 30) xx, pp Ntab = Ntab + 1 xkt (Ntab) = xx *hsmall Pkt (Ntab) = pp / xx**3 * twop ! k^3Pk -> Pk GOTO 10 30 Write (*,*) Write (*,*) ' Read ', Ntab, ' lines from WHUpk.dat table' Else ! read header Do i = 1, 9 READ (50, '(a)', end = 100, err = 100) Line EndDo Ntab = 0 12 READ (50, *, end = 32, err = 32) xx, pp Ntab = Ntab + 1 xkt (Ntab) = xx *hsmall Pkt (Ntab) = pp ! Pk GOTO 12 32 Write (*,*) Write (*,*) ' Read ', Ntab, ' lines from P(k) table ' EndIf ! end case iFlag If (Ntab.le.1) stop 'wrong table for p(k)' StepK = log10 (xkt (Ntab) / xkt (1) )/(Ntab-1) alog0 = log10 (xkt (1) ) ! test that spacing is the same Do k = 2, Ntab - 1 ss = log10 (xkt (k + 1) / xkt (k) ) If (abs (ss / StepK - 1.) .gt.2.e-2) Then Write (*,*) ' error in K spacing. k=', k, xkt (k + 1),xkt (k) STOP EndIf EndDo Return 100 Write (*,*) ' Error reading table of P(k)' STOP END SUBROUTINE ReadTable