program sM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Andrew Zenter & J.S. Bullock CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Generates sigma(M) look-up table. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CCCCC Common block with cosmological parameters CCCCCCCCCC real*8 cobe_offset common/offset/cobe_offset real*8 OmegaM, OmegaL, OmegaBh2, hubble real*8 N_m_nu, m_nus, tilt, run, Rfree common /cosmo/ OmegaM, OmegaL, OmegaBh2, hubble, & N_m_nu, m_nus, tilt, run, Rfree CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC real*8 test real*8 sigM, mass, R, z, sig8 real*8 logmin, logmax, dm, nsint, pwr real*8 Mwdm, Mf character*4 model integer nsample, lcv external sigM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC print*,' ' print*,'*****************************************' print*,'* *' print*,'* sigma(M) *' print*,'* *' print*,'*****************************************' print*,' ' print*,'***** Reading Input Parameters **********' print*,' ' open(unit=11,file='sM_in.dat',status='unknown') read(11,*) model read(11,*) logmin, logmax, nsample write(*,*) 'The model name is ',model,'.' print*,' ' open (unit=12, file = 'model.params.'//model) write(*,*) 'Compute sigma(M) at ',nsample write(*,*) 'equally spaced points in the range ' write(*,*) 'log(Mmin) = ',logmin,' to ' write(*,*) 'log(Mmax) = ',logmax,' .' print*,' ' write(12,*) 'Compute sigma(M) at ',nsample write(12,*) 'equally spaced points in the range ' write(12,*) 'log(Mmin) = ',logmin,' to ' write(12,*) 'log(Mmax) = ',logmax,' .' z = 0. write(*,*) 'Compute sigma(M) at redshift ',z print*,' ' write(12,*) 'Compute sigma(M) at redshift ',z read(11,*) OmegaM, OmegaBh2, hubble read(11,*) N_m_nu, m_nus read(11,*) tilt, run read(11,*) Mwdm ! in Kev read(11,*) sig8 if (Mwdm.gt.0.001) then Rfree = 0.2*(OmegaM*hubble**2)**.3333 Rfree = Rfree*Mwdm**(-4./3.) Rfree = Rfree*hubble !in Mpc/h else Rfree = 0. endif Mf = 3.7e14*Rfree**3*OmegaM OmegaL = 1.d0 - OmegaM print*,'The cosmological parameters are : ' print*,' ' write(*,*) '(1) Omega_Matter = ',OmegaM write(*,*) '(2) Omega_Baryon*h**2 = ',OmegaBh2 write(*,*) '(3) Hubble parameter, h = ',hubble write(*,*) '(4) The number of massive, mass' write(*,*) ' degenerate neutrinos is ',N_m_nu write(*,*) '(5) The neutrino mass is ',m_nus write(*,*) '(6) Spatial flatness is assumed so ' write(*,*) 'that Omega_Lambda = ',OmegaL print*,' ' print*,'The initial power spectrum has : ' write(*,*) '(1) n(k_cobe) = ',tilt write(*,*) '(2) dn(k)/dln(k) = ',run write(*,*) 'WDM particle mass in Kev = ', Mwdm write(*,*) 'This implies Rfree in Mpc/h = ', Rfree, Rfree/hubble write(*,*) 'With log10 Mfree Msun/h = ', log10(Mf) print*,' ' print*,'***** Input parameters specified **********' print*,' ' write(12,*) 'The cosmological parameters are : ' write(12,*) '(1) Omega_Matter = ',OmegaM write(12,*) '(2) Omega_Baryon*h**2 = ',OmegaBh2 write(12,*) '(3) Hubble parameter, h = ',hubble write(12,*) '(4) The number of massive, mass' write(12,*) ' degenerate neutrinos is ',N_m_nu write(12,*) '(5) The neutrino mass is ',m_nus write(12,*) '(6) Spatial flatness is assumed so ' write(12,*) 'that Omega_Lambda = ',OmegaL write(12,*) 'The initial power spectrum has : ' write(12,*) '(1) n(k_cobe) = ',tilt write(12,*) '(2) dn(k)/dln(k) = ',run write(12,*) 'WDM particle mass in Kev = ', Mwdm write(12,*) 'This implies Rfree in Mpc/h = ', Rfree pause open(unit=21,file='sigma_'//model,status='unknown') 211 format(f12.5,' ',f8.5,' ',I6,' ',f8.5) 212 format(f9.5,' ',f12.7) 931 format('COBE-norm Sigma_8 = ',f7.3,' *****') 932 format('****using Sigma_8 = ',f7.3,' *****') CCCCCCCCCCCC Sigma_8 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cobe_offset = 1.d0 test = sigM(8.d0,0.d0) print*,' ' write(*,931) test print*,' ' write(12,*) 'COBE-normalized sigma_8 = ', test if (sig8.gt.0.) then cobe_offset = sig8/test endif test = sigM(8.d0,0.d0) print*,' ' write(*,932) test print*,' ' write(12,*) 'using sigma_8 = ', test CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC nsint = nsample*1.d0 dm = (logmax - logmin)/nsint write(21,211) logmin, logmax, nsample, dm CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC print*,'***** Calculating sigma(M) values *****' print*,' ' do lcv=0,nsample pwr = logmin + lcv*dm mass = 10.d0**pwr R = 9.510362d-5*((mass/OmegaM)**(1.d0/3.d0)) test = sigM(R,z) write(21,212) pwr,test if (mod(lcv,40).eq.0) write(*,212) pwr, test enddo print*,' ' print*,'***** Finished *****' print*,' ' write(*,*) 'Output written to the file ' write(*,*) 'sigma_'//model print*,' ' print*,'********************' print*,' ' end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC *************************************************************************** CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Subroutines and Functions C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Calculate Quantity To Be Integrated To Get sigma(M) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC function dsig(lnk, pars, npars) integer npars real*8 pars(npars) real*8 cobe_norm, transfer, dsig external cobe_norm, transfer real*8 k_star, k, n, delta, ind, lnk real*8 wind, k_norm, T, d_star, arg, z real*8 OmegaM, OmegaL, OmegaBh2, hubble real*8 N_m_nu, m_nus, tilt, run, Rfree common /cosmo/ OmegaM, OmegaL, OmegaBh2, & hubble, N_m_nu, m_nus, tilt, run, Rfree R = pars(1) z = pars(2) k_star = 0.0023d0 k = dexp(lnk) k_norm = k/k_star n = tilt d_star = cobe_norm(tilt,OmegaM,0.d0) T = transfer(k,z,OmegaM,OmegaBh2,hubble,N_m_nu,m_nus) n = n + 0.5*run*log(k/k_star) ind = n + 3.d0 delta = (2401.d0*(d_star**2.d0))*((T**2.d0) & *(k_norm**ind)) c... Simplistic WDM approximation (see, e.g. Eke et al. 2001) if ( Rfree .gt.0.001) then aa = Rfree ! Rfree is assumed to be in units of Mpc/h delta = delta*exp(-k*aa - (k*aa)**2) endif arg = R*k c... window function if (arg.lt.1.d-5) then wind = 1.d0 elseif (arg.lt.1.d6) then wind = 3.d0*((dsin(arg) - arg*dcos(arg))/arg**3.d0) else wind = 0.d0 endif wind = wind**2.d0 dsig = delta*wind return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Linear Growth Factor in Flat, LambdaCDM Cosmology C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC function growth(z,O_m) real*8 z, O_m real*8 a, Ol, x, B, C real*8 betai, growth external betai a = 1.d0/(1.d0+z) Ol = 1.d0-O_m x = Ol*a**3.d0/(O_m + Ol*a**3.d0) B = 0.5d0/0.6d0 C = 0.2d0/0.3d0 f1 = B*betai(B,C,x)*((O_m/Ol)**(1.d0/3.d0)) f2 = dsqrt(1.d0+(O_m/(Ol*a**3.d0))) growth = f1*f2 return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Calculate COBE Normalization At COBE Pivot Point. C This is for a flat LambdaCDM cosmology. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC function cobe_norm(n,Om,r) real*8 n,Om,r, cobe_norm real*8 Ol, f, pwr, g, d0 Ol = 1.d0-Om f=0.75d0-0.13d0*(Ol**2.d0) pwr = -0.8d0 - 0.05d0*log(Om) g = 1.d0-0.18d0*(1.d0-n)*Ol - 0.03d0*r*Ol d0 = 1.91d-5*(7.d0**((n-1.d0)/2.d0)) cobe_norm = d0*(exp(1.01d0*(1.d0-n))/ & dsqrt(1.d0+r*f))*(Om**pwr)*g return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Transfer Function Fitting Formula From Eisenstein and Hu. C Omega_cb does NOT include the massive neutrinos!!!!!!!!!! C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC function transfer(kin,z,O_cb,O_bh2,h,N_nu,m_nu) real*8 kin, z, O_cb, O_bh2, h, N_nu, m_nu real*8 growth, transfer external growth real*8 k, O_nuh2, O_nu, O_b, O_c, O_m, O_L, f_c, f_b real*8 f_cb, theta, O_mh2, y_eq, b1, b2, y_d, s, q, y real*8 g2, O_m_z, O_L_z, D1, D1_0 real*8 p_c, p_cb, y_fs, alpha, rtalpha, G_eff, q_eff real*8 beta, L, C, T_sup, B, T_mas, q_nu, f_nu, f_nub k = h*kin h2 = h**2.d0 CCCCCCC First Get HDM Energy Density From Neutrino Mass O_nuh2 = N_nu*m_nu/91.5d0 O_nu = O_nuh2/h2 O_b = O_bh2/h2 O_c = O_cb-O_b O_m = O_cb + O_nu O_L = 1.d0 - O_cb - O_nu f_c = O_c/O_m f_b = O_b/O_m f_cb = (O_c + O_b)/O_m CCCCCCCCC CMB Temp in Units of 2.7K CCCCCCCCCCCCCCCCCCCCC theta = 1.0093d0 O_mh2 = O_m*h2 CCCCCCCCC Epoch of matter radiation Equality y_eq = 2.50d4*O_mh2/(theta**4) + 1.d0 b1 = (.313d0*O_mh2**-0.419d0)*(1.d0+0.607d0*O_mh2** & (0.674d0)) b2 = 0.238d0*O_mh2**0.223d0 CCCCCCCCC Epoch of Decoupling y_d = 1291.d0*(1.d0+b1*O_bh2**b2)*(O_mh2**0.251d0/ & (1.d0+0.659d0*O_mh2**0.828d0)) + 1.d0 y_d = y_eq/y_d CCCCCCCCC Sound Horizon At Decoupling s = 44.5d0*dlog(9.83d0/O_mh2)/ & dsqrt(1.d0+10.d0*(O_bh2**0.75d0)) q = theta**2*k/O_mh2 y = 1.d0+z g2 = O_m*y**3.d0 + (1.d0-O_m-O_L)*y**2.d0 + O_L O_m_z = O_m*y**3/g2 O_L_z = O_L/g2 D1 = 1.724733d0*y_eq*growth(z,O_m) D1_0 = 1.724733d0*y_eq*growth(0.d0,O_m) CCCCCCCCC First, we have a simplified case with no HDM if (m_nu.lt.1.d-4) then O_nu = 0 p_c = 0.25d0*(5.d0-dsqrt(1.d0+24.d0*f_c)) p_cb = 0.25d0*(5.d0-dsqrt(1.d0+24.d0*f_cb)) y_fs = 0.d0 alpha = (f_c/f_cb)*((5.d0-2.d0*(p_c+p_cb))/ & (5.d0-4.d0*p_cb))*(1.d0-0.553d0*f_b+0.126d0*f_b**3)* & ((1.d0+y_d)**(p_cb-p_c))*(1.d0+.5d0*(p_c-p_cb)*(1.d0+1.d0/ & ((3.d0-4.d0*p_c)*(7.d0-4.d0*p_cb)))/(1.d0+y_d)) rtalpha = dsqrt(alpha) G_eff=O_mh2*(rtalpha + (1.d0-rtalpha)/ & (1.d0+((0.43d0*s)*k)**4)) q_eff = (theta**2)*k/G_eff beta = 1.d0/(1.d0-0.949d0*f_b) L = log(2.718282d0+(1.84d0*beta*rtalpha)*q_eff) C = 14.4d0+325.d0/(1.d0+60.5d0*(q_eff**1.11d0)) T_sup = L/(L + C*(q_eff**2)) transfer = T_sup*(D1/D1_0) return else f_nu = O_nu/O_m f_nub = (O_nu+O_b)/O_m p_c = 0.25d0*(5.d0-dsqrt(1.d0+24.d0*f_c)) p_cb = 0.25d0*(5.d0-sqrt(1.d0+24.d0*f_cb)) y_fs = (17.2d0*f_nu*(1.d0+.488d0*f_nu**(-7.d0/6.d0))* & (N_nu/f_nu)**2)*(q**2) D_cb = (D1**(1.d0-p_cb))*(1.d0+(D1/(1.d0+y_fs))**0.7d0) & **(p_cb/.7d0) alpha = (f_c/f_cb)*((5.d0-2.d0*(p_c+p_cb))/ & (5.d0-4.d0*p_cb))*(1.d0-0.553d0*f_nub+0.126d0*f_nub**3)* & ((1.d0+y_d)**(p_cb-p_c))*(1.d0+0.5d0*(p_c-p_cb)*(1.d0+1.d0/ & ((3.d0-4.d0*p_c)*(7.d0-4.d0*p_cb)))/(1.d0+y_d))/ & (1.d0-0.193d0*dsqrt(f_nu*N_nu)+0.169d0*f_nu*(N_nu**0.2d0)) rtalpha = dsqrt(alpha) G_eff = O_mh2*(rtalpha+(1.d0-rtalpha)/ & (1.d0+((.43d0*s)*k)**4)) q_eff = theta**2*k/G_eff beta = 1.d0/(1.d0-0.949d0*f_nub) C = 14.4d0+325.d0/(1.d0+60.5d0*(q_eff**1.11d0)) L = dlog(2.718282d0+(1.84d0*beta*rtalpha)*q_eff) T_sup = L/(L+C*(q_eff**2)) CCCCCCCCC There are two extra corrections for HDM q_nu = (3.92d0*dsqrt(N_nu)/f_nu)*q B = 1.d0+(1.2d0*(f_nu**.64d0)* & (N_nu**(.3d0+.6d0*f_nu)))/(q_nu**(-1.6d0)+q_nu**(0.8d0)) T_mas = T_sup*B transfer = T_mas*D_cb/D1_0 return endif end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Find sigma(R) on some length scale R in h^-1 Mpc C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC function sigM(R,z) real*8 R,z,ss, sigM real*8 dsig, pars(2) real*8 INTEGRATE real*8 cobe_offset common/offset/cobe_offset external dsig, INTEGRATE pars(1) = R pars(2) = z ss = INTEGRATE(dsig,pars,2,-15.d0,30.d0,0.1d0,1.d-7) ss = dsqrt(ss) sigM = ss*cobe_offset return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Integration Stuff C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Special Functions From Numerical Recipes C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FUNCTION gammln(xx) real*8 gammln,xx INTEGER j real*8 ser,stp,tmp,x,y,cof(6) SAVE cof,stp DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, & 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, & -.5395239384953d-5,2.5066282746310005d0/ x=xx y=x tmp=x+5.5d0 tmp=(x+0.5d0)*log(tmp)-tmp ser=1.000000000190015d0 do j=1,6 y=y+1.d0 ser=ser+cof(j)/y enddo gammln=tmp+dlog(stp*ser/x) return END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C beta function C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FUNCTION betacf(a,b,x) INTEGER MAXIT real*8 betacf,a,b,x,EPS,FPMIN PARAMETER (MAXIT=100,EPS=3.d-7,FPMIN=1.d-30) INTEGER m,m2 real*8 aa,c,d,del,h,qab,qam,qap CCCCCCCCCCCC qab=a+b qap=a+1.d0 qam=a-1.d0 c=1.d0 d=1.d0-qab*x/qap if(abs(d).lt.FPMIN)d=FPMIN d=1.d0/d h=d do m=1,MAXIT m2=2*m aa=m*(b-m)*x/((qam+m2)*(a+m2)) d=1.d0+aa*d if(abs(d).lt.FPMIN)d=FPMIN c=1.d0+aa/c if(abs(c).lt.FPMIN)c=FPMIN d=1.d0/d h=h*d*c aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2)) d=1.d0+aa*d if(abs(d).lt.FPMIN)d=FPMIN c =1.d0+aa/c if(abs(c).lt.FPMIN)c=FPMIN d=1.d0/d del=d*c h=h*del if(abs(del-1.d0).lt.EPS)goto 1 enddo pause 'a or b too big, or MAXIT too small in betacf' 1 betacf=h return END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Incomplete Beta Function C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FUNCTION betai(a,b,x) real*8 betai,a,b,x real*8 bt,betacf,gammln if(x.lt.0.d0.or.x.gt.1.d0)pause 'bad argument x in betai' if(x.eq.0.d0.or.x.eq.1.d0)then bt=0.d0 else bt=exp(gammln(a+b)-gammln(a)-gammln(b) & +a*log(x)+b*log(1.d0-x)) endif if(x.lt.(a+1.d0)/(a+b+2.d0)) then betai=bt*betacf(a,b,x)/a return else betai=1.d0-bt*betacf(b,a,1.d0-x)/b return endif END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCC######################################################################### CCCCCCCCC######################################################################### CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Better Integrator. Again From Numerical Recipes C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC function INTEGRATE(FUNC,fp,np,a,b,dxinit,eps) implicit none integer np real*8 integrate real*8 FUNC real*8 a, b, eps, dxinit, fp(np) external FUNC integer maxsteps parameter(maxsteps=100000000) real*8 x, dx, dxnext, y, dydx, yscale integer Nstep x = a dx = dxinit y = 0.d0 Nstep = 0 do while ((x-b)*(b-a).lt.0.d0.and.Nstep.lt.maxsteps) Nstep = Nstep + 1 dydx = FUNC(x,fp,np) yscale = max(abs(y) + abs(dx*dydx), 1.d-12) if ((x+dx-b)*(x+dx-a).gt.0.d0) ! If stepsize overshoots, decrease it. & dx = b - x call RUNGE5VAR(y,dydx,x,dx,eps,yscale,dxnext,FUNC,fp,np) dx = dxnext enddo if (Nstep.ge.maxsteps) & write (*,*) 'WARNING: failed to converge in INTEGRATE.' INTEGRATE = y return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Fifth-order Runge-Kutta step with monitoring of local truncation error c to ensure accuracy and adjust stepsize. Input are the dependent c variable y and its derivative dydx at the starting value of the c independent variable x. Also input are the stepsize to be attempted c htry, the required accuracy eps, and the value yscale, against which the c error is scaled. On output, y and x are replaced by their new values. c hdid is the stepsize that was actually accomplished, and hnext is the c estimated next stepsize. DERIVS is the user-supplied routine that c computes right-hand-side derivatives. The argument fparm is for an c optional second argument to DERIVS (NOT integrated over). c c c by M.A.K. Gross c inherited by me from J. S. Bullock SUBROUTINE RUNGE5VAR(y,dydx,x,htry,eps,yscale,hnext,DERIVS, & fp,np) implicit none integer np real*8 fp(np) real*8 eps,hnext,htry,x,dydx,y,yscale,DERIVS external DERIVS real*8 errmax,h,hold,htemp,xnew,yerr,ytemp real*8 safety,pgrow,pshrink,errcon parameter (safety = 0.9d0) parameter (pgrow = -0.2d0) parameter (pshrink = -0.25d0) parameter (errcon = 1.89d-4) h = htry ! Set stepsize to initial accuracy. errmax = 10.d0 do while (errmax.gt.1.d0) call RUNGE(y,dydx,x,h,ytemp,yerr,DERIVS,fp,np) errmax = abs(yerr/yscale)/eps ! Scale relative to required accuracy. if (errmax.gt.1.d0) then ! Truncation error too large; reduce h htemp = safety*h*(errmax**pshrink) hold = h h = sign(max(abs(htemp),0.1d0*abs(h)),h) ! No more than factor of 10 xnew = x + h if (xnew.eq.x) then write (*,*) 'WARNING: ', & 'Stepsize underflow in RUNGE5VAR().' h = hold errmax = 0.d0 endif endif enddo if (errmax.gt.errcon) then hnext = safety*h*(errmax**pgrow) else hnext = 5.d0*h ! No more than factor of 5 increase. endif x = x + h y = ytemp return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Given values for a variable y and its derivative dydx known at x, use c the fifth-order Cash-Karp Runge-Kutta method to advance the solution c over an interval h and return the incremented variables as yout. Also c return an estimate of the local truncation error in yout using the c embedded fourth order method. The user supplies the routine c DERIVS(x,y,dydx), which returns derivatives dydx at x. c c by M.A.K. Gross c SUBROUTINE RUNGE(y,dydx,x,h,yout,yerr,DERIVS,fp,np) implicit none integer np real*8 h,x,dydx,y,yerr,yout,DERIVS,fp(np) external DERIVS real*8 ak3, ak4, ak5 ,ak6 real*8 a2,a3,a4,a5,a6 real*8 c1,c3,c4,c6,dc1,dc3,dc4,dc5,dc6 parameter(a2 = 0.2d0) parameter(a3 = 0.3d0) parameter(a4 = 0.6d0) parameter(a5 = 1.d0) parameter(a6 = 0.875d0) parameter(c1 = 37.d0/378.d0) parameter(c3 = 250.d0/621.d0) parameter(c4 = 125.d0/594.d0) parameter(c6 = 512.d0/1771.d0) parameter(dc1 = c1 - 2825.d0/27648.d0) parameter(dc3 = c3 - 18575.d0/48384.d0) parameter(dc4 = c4 - 13525.d0/55296.d0) parameter(dc5 = -277.d0/14336.d0) parameter(dc6 = c6 - 0.25d0) ak3 = DERIVS(x+a3*h,fp,np) ak4 = DERIVS(x+a4*h,fp,np) ak5 = DERIVS(x+a5*h,fp,np) ak6 = DERIVS(x+a6*h,fp,np) c c Estimate the fifth order value. c yout = y + h*(c1*dydx + c3*ak3 + c4*ak4 + c6*ak6) c c Estimate error as difference between fourth and fifth order c yerr = h*(dc1*dydx + dc3*ak3 + dc4*ak4 + dc5*ak5 + dc6*ak6) return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C c ---------------- FUNCTION erfc(x) c ---------------- c Taken from Numerical Recipes. c implicit none real*8 erfc,x real*8 gammp,gammq external gammp,gammq if(x.lt.0.d0) then erfc=1.d0+gammp(5.d-1,x**2) else erfc=gammq(5.d-1,x**2) endif return END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FUNCTION gammp(a,x) implicit none real*8 a,gammp,x real*8 gammcf,gamser,gln if (x.lt.0.d0.or.a.le.0.d0) & pause 'bad arguments in gammp' if (x.lt.a+1.d0) then call gser(gamser,a,x,gln) gammp=gamser else call gcf(gammcf,a,x,gln) gammp=1.d0-gammcf endif return END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FUNCTION gammq(a,x) implicit none real*8 a,gammq,x real*8 gammcf,gamser,gln if (x.lt.0.d0.or.a.le.0.d0) & pause 'bad arguments in gammq' if (x.lt.a+1.d0) then call gser(gamser,a,x,gln) gammq=1.d0-gamser else call gcf(gammcf,a,x,gln) gammq=gammcf endif return END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE gcf(gammcf,a,x,gln) implicit none INTEGER ITMAX real*8 a,gammcf,gln,x,EPS,FPMIN PARAMETER (ITMAX=100,EPS=3.d-7,FPMIN=1.d-30) INTEGER i real*8 an,b,c,d,del,h,gammln external gammln gln=gammln(a) b=x+1.d0-a c=1.d0/FPMIN d=1.d0/b h=d do i=1,ITMAX an=-i*(i-a) b=b+2.d0 d=an*d+b if(abs(d).lt.FPMIN) d=FPMIN c=b+an/c if(abs(c).lt.FPMIN) c=FPMIN d=1.d0/d del=d*c h=h*del if(abs(del-1.d0).lt.EPS) goto 1 enddo continue pause 'a too large, ITMAX too small in gcf' 1 gammcf=exp(-x+a*log(x)-gln)*h return END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE gser(gamser,a,x,gln) implicit none INTEGER ITMAX real*8 a,gamser,gln,x,EPS PARAMETER (ITMAX=100,EPS=3.d-7) INTEGER n real*8 ap,del,sum,gammln external gammln gln=gammln(a) if(x.le.0.d0) then if(x.lt.0.d0) pause 'x<0 in gser' gamser=0.d0 return endif ap=a sum=1.d0/a del=sum do n=1,ITMAX ap=ap+1.d0 del=del*x/ap sum=sum+del if(abs(del).lt.abs(sum)*EPS) goto 1 enddo continue pause 'a too large, ITMAX too small in gser' 1 gamser=sum*exp(-x+a*log(x)-gln) return END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C (C) Copr. 1986-92 Numerical Recipes Software .-35?421.1-9. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC