c======================================================================= c cvir.f c----------------------------------------------------------------------- c by James Bullock c using routines made available by A. Kravtsov, J. Navarro and M. Gross c --------------------------------------------------------------------- c Implements toy model outlined in Bullock et al. (2001). c Returns cvir = Rvir/Rs given a halo mass Mvir and redshift z c --------------------------------------------------------------------- c The program is set up to return cvir for a given Mvir and z c --------------------------------------------------------------------- c --------------------------------------------------------------------- c The two free parameters of the model, F & K, have been set to match c the results of N-body simulations for the LCDM and SCDM cosmologies. c ---------------------------------------------------------------------- c Here we use F=0.01 and K=3.75. These are the values appropriate for the c correct power spectrum we adopt here. c ---------------------------------------------------------------------- c============================================================================ implicit none real dum real omega0,lambda real mcall real cvir,zcoll, mvir, z, F, K,zl,zh real cvir_low, cvir_high, Klow, Khigh !+/- 1 sigma values real sMmin, sMmax, dsM, sigm, om0,Par, rf real sigma_8, ns, Rfree, h character*4 model integer nsample common / SIGM / sMmin , ! see sigma_init_read & sMmax , & dsM , & sigm(10000),sigma_8, ns, Rfree,h, & nsample, om0, Par(6), rf open(unit=57,file='sM_in.dat',status='old') read(57,*) model read(57,*) read(57,*) omega0, dum, h om0 = omega0 lambda = 1.0 - omega0 read(57,*) read(57,*) ns, dum read(57,*) read(57,*) sigma_8 write(*,*) 'Mass of halo in Msun/h?' read(*,*) mvir write(*,*) 'redshift?' read(*,*) z K = 3.75 F = 0.01 open (unit=21, file='sigma_'//model, status='unknown') call sigma_init_read ( ) ! reads in mass power spectra c... determine median cvir mcall = mvir call compute_cvir (cvir, zcoll, omega0, lambda, & mcall, z, F, K) call compute_Vmax (cvir,mvir) c... determine +/- 1-sigma values for cvir Klow = 10**(log10(K) - .14) Khigh = 10**(log10(K) + .14) call compute_cvir (cvir_low, zcoll, omega0, lambda, & mcall, z, F, Klow) zl = (1. + zcoll)*cvir_low/cvir - 1. call compute_cvir (cvir_high, zcoll, omega0, lambda, & mcall, z, F, Khigh) zh = (1. + zcoll)*cvir_high/cvir - 1. c... write output write(*,*) write(*,11) write(*,5) omega0, lambda, h,sigma_8, ns write(*,6) write(*,7) K,F write(*,*) write(*,8) mcall, z write(*,*) write(*,13) zl, zh write(*,12) zcoll write(*,*) write(*,9) cvir_low, cvir_high write(*,10) cvir write(*,11) write(*,*) 5 format ( 1x,'Om_m = ', f3.1,1x, & ' Om_l = ',f3.1,1x,' h = ',f4.2,1x, & 'sig_8 = ', f4.2,1x,' n_s = ',f4.2) 6 format (1x, & '---------------------------------------------------------') 7 format (1x, 'Using parameters: K = ', f4.2,1x,'F = ', f5.3) 8 format (1x, 'Mvir = ', g20.10,1x,'Msun/h at z = ', f4.2,1x) 9 format (1x, 'cvir = ', f10.5,' - ', f10.5,' (68% range)') 10 format (1x, 'median cvir = ',f10.5) 12 format (1x, 'median z_c =', f6.2) 13 format (1x, 'z_c =', f6.2, ' -- ', f6.2, ' (68% range)') 11 format (1x, & '=========================================================') STOP END c c================================================================= subroutine compute_cvir(c,zcoll,omega0,lambda,mvir,z0,F,K) c------------------------------------------------------------------ c cvir is supplied for mvir(z) in units of (Msun h^-1) c================================================================== implicit none real sMmin, sMmax, dsM, sigm, om0,Par, rf integer nsample real sigma_8, ns, Rfree,h common / SIGM / sMmin , ! see sigma_init_read & sMmax , & dsM , & sigm(10000), sigma_8, ns, Rfree, h, & nsample, om0, Par(6), rf real lambda, mvir, mstar real c, zcoll, omega0,z0 real a, z, astar real diff, msg, dcz real F, K c.... find expansion factor, astar, when F*mvir was M_* a = 1./(1.+z0) call findmstar(z0, omega0, lambda, dcz, mstar) astar = 10**(log10(a) - (mstar-log10(mvir))/15.) if (astar.gt.1.) astar = 1. 111 continue z = 1./astar - 1. call deltacrit(z,omega0,lambda,dcz) dcz = dcz if (sigm(1).lt.dcz) then astar = 1.2*astar goto 111 endif msg = log10(F*mvir) mstar = 1. diff = mstar/5. do while((abs(diff)/abs(mstar).gt.0.002)) 112 continue z = 1./astar - 1. call deltacrit(z,omega0,lambda,dcz) dcz = dcz if (sigm(1).lt.dcz) then write(*,*) z, astar, dcz, sigm(1) astar = 1.2*astar goto 112 endif call findmstar(z, omega0, lambda, dcz, mstar) diff = mstar - msg astar = 10**(log10(astar) - diff/80.) enddo c... zcoll: typical collapse redshift for halos of mass F*mvir zcoll = 1/astar - 1. c... This determines c c = K*((1.0+zcoll)/(1.0+z0)) return end c c========================================================= subroutine findmstar(z, omega0, lambda, dcz, mstar) c ---------------------------------------------------- c returns the mstar, typical collapsing mass at c redshift z, in Msun h^-1 c========================================================= implicit none real sMmin, sMmax, dsM, sigm, om0,Par, rf integer nsample real sigma_8, ns, Rfree,h common / SIGM / sMmin , ! see sigma_init_read & sMmax , & dsM , & sigm(10000), sigma_8, ns, Rfree, h, & nsample, om0, Par(6), rf real z, omega0, lambda, dcz, mstar, mstar_g, sig real sigma, diff real mguess call deltacrit(z,omega0,lambda,dcz) dcz = dcz mstar_g = 15. ! first guess mguess = 10**mstar_g if (sigm(1).lt.dcz) then write(*,*) 'Mass is out of range for this power spectrum' pause endif sig=sigma(mguess) diff = log10(sig) - log10(dcz) do while (abs(diff).gt.0.002) mstar_g = mstar_g +5.*diff mguess = 10**mstar_g sig = sigma(mguess) diff = log10(sig) - log10(dcz) enddo mstar = mstar_g return end c c=============================================================== subroutine sigma_init_read ( ) c --------------------------------------------------------- c reads rms fluctuations sigma(M) c where Mass M in Msun/h c c=============================================================== real sMmin, sMmax, dsM, sigm, om0,Par,rf integer i, nsample real sigma_8, ns, Rfree,h real rho_c parameter (rho_c = 2.7755e11 ) ! critical density in h**2 Msun common / SIGM / sMmin , ! see sigma_init_read & sMmax , & dsM , & sigm(10000), sigma_8, ns, Rfree, h, & nsample, om0,Par(6),rf read (21,*) sMmin, sMmax, nsample, dsM do i = 1, nsample read(21,*) dum, sigm(i) c write(*,*) i, sigm(i) enddo return end c============================================================== real function sigma ( sM ) c --------------------------- ---------------------------- c modified from a function supplied by A. Kravtsov c -------------------------------------------------------- c get rms fluctuation on a mass scale sM c sMmin , sMmax , dsM , nsample are read from file 57 c [sM] = [Msun] - physical mass (not multiplied by h_100) c=============================================================== real sM, sMmin, sMmax, dsM, sigm, om0, Par, rf integer nsample real sigma_8, ns, Rfree,h common / SIGM / sMmin , ! see sigma_init_read & sMmax , & dsM , & sigm(10000), sigma_8, ns, Rfree, h, & nsample, om0, Par(6), rf ibin1 = (log10(sM) - sMmin) / dsM + 1 IF ( (ibin1.gt.1).and.(ibin1.lt.nsample) ) THEN sM1 = 10.0**(sMmin + dsM*ibin1 - 0.5*dsM) if ( sM .lt. sM1 ) then ibin2 = ibin1 - 1 sM2 = 10.d0**(sMmin + dsM*ibin2 - 0.5d0*dsM) a = (sigm(ibin2) - sigm(ibin1)) / (sM2-sM1) b = sigm(ibin1) - a * sM1 sigma = a * sM + b else ibin2 = ibin1 + 1 sM2 = 10.d0**(sMmin + dsM*ibin2 - 0.5d0*dsM) a = (sigm(ibin2) - sigm(ibin1)) / (sM2-sM1) b = sigm(ibin1) - a * sM1 sigma = a * sM + b endif ELSE write(*,*) ibin1, sM, sMmin, sMmax write(*,*) & 'mass out of range, change sMmin or sMmax in sigma_init routine' pause ENDIF return end c========================================================== subroutine tofzl(z,omega0,lambda,tz,xxz) c ----------------------------------------------------- c Modified from a routine made available by Julio Navarro c ----------------------------------------------------- c returns the time at any z in units of H0^-1 for a model with c Omega=1, Omega < 1, or Omega+Lambda=1 c=========================================================== real lambda c--- Omega=1 if (omega0.eq.1.0) then tz=(2./3)*(1.0+z)**(-1.5) return end if c c-- Omega<1, Lambda=0 if (omega0.lt.1.0 .and. lambda.eq.0) then step=1.e-3 z1=1.0+z a=1.0/z1 xxmax=1.e4 do xx=step,xxmax,step atmp=0.5*omega0*(cosh(xx)-1.0)/(1.0-omega0) if (atmp.gt.a) then xxz=xx go to 1 end if end do 1 continue if (atmp.lt.a) stop 'SUB. TOFZL:t_o: Increase xxmax..' tz=0.5*omega0*(sinh(xxz)-xxz)/(1.0-omega0)**1.5 return end if c c-- Lambda+Omega=1, Peebles 1994 eq.13.20 if (lambda+omega0.eq.1.0) then z1=1.0+z omega1=1.0-omega0 const0=2./(3*sqrt(omega1)) param=sqrt(omega1/omega0)*(z1**(-3./2)) tz=const0*log(param+sqrt(1.0+param*param)) return end if end c c=============================================================== subroutine lingrowth(z,omega0,lambda,omega,b,hubble) c ---------------------------------------------------------- c By Julio Navarro c ----------------------------------------------------------- c c This subroutine returns the value of the growth factor of c linear density perturbations at a given redshift. c c The Hubble constant at that redshift is also returned, in c units such that: H_0 = H(z=0) = 1 c c The value of Omega(z) is also returned. c b(z) is from Peebles (1980), eq.11.15. c=============================================================== real lambda double precision func2,x,int,aofx,aofxn double precision xn,w, zero common /bzero/ b0 external func2 c z1 = 1.0+z c c--Compute omega and Hubble's constant at this redshift. if(omega0.lt.1 .and. lambda.eq.0.0) then omega = omega0*z1/(1.0+omega0*z) hubble=z1*sqrt(1.0+omega0*z) ! in units of h0=1.0 else if (omega0.eq.1.0) then omega = omega0 hubble=z1**1.5 ! in units of h0=1.0 else if (omega0+lambda.eq.1.0) then omega=omega0*z1**3/(1.0-omega0+omega0*z1**3) hubble=sqrt(1.0-omega0+omega0*z1**3) else stop 'LINGROWTH: Omega+Lambda .ne. 1' end if c c-- Compute b(z) c--a) Omega<1 if (omega0.lt.1 .and. lambda.eq.0.0) then call tofzl(0.0,omega0,lambda,t0,xx0) b0=(3.0*sinh(xx0)*(sinh(xx0)-xx0)/(cosh(xx0)-1.)**2.)-2.0 call tofzl(z,omega0,lambda,tz,xx) b=(3.0*sinh(xx)*(sinh(xx)-xx)/(cosh(xx)-1.)**2.)-2.0 c c--b) Omega=1 else if (omega0.eq.1.0 .and. lambda.eq.0.0) then t0=(2./3) b0=(3*t0/2)**(2./3) tz=(2./3)*z1**(-1.5) b=(3*tz/2)**(2./3) c c--c) Omega+Lambda=1. In this case b(z) is normalized so that b(0)=1. else if (omega0+lambda.eq.1.0) then a=1./z1 w = omega0**(-1.0) - 1.0 if (abs(omega0-1) .lt. 1.e-5) then b = a else xn = (2.0*w)**(1.0/3) int = 0. zero = 0.0 call simp(func2,zero,xn,int) aofxn = ((xn**3.0+2.0)**0.5)*(int/xn**1.5) b0=aofxn x = a*xn call simp(func2,zero,x,int) aofx = ((x**3+2)**0.5)*(int/x**1.5) b = aofx/aofxn end if c end if c end c c================================================================ subroutine deltacrit(z,omega0,lambda,deltac) c ----------------------------------------------------------- c by Julio Navarro c ----------------------------------------------------------- c c This subroutine computes the critical density threshold delta_c c ---see Lacey& Cole 1993, MNRAS 262, 627 for a definition -- as a c function of redshift. This version works for C C a) Omega=1 Lambda=0. C b) Omega<1 Lambda=0. C c) Omega+ Lambda=1. C================================================================= real dum real lambda common /bzero/ b0 c pi=4*atan(1.0) c c-- Omega=1, Lambda=0 if (omega0.eq.1.0 .and. lambda.eq.0.0) then deltac=(3./20)*(12*pi)**(2./3)*(1.0+z) return end if c c-- Omega < 1, Lambda=0. if (omega0.lt.1.0 .and. lambda.eq.0) then tomega=pi*omega0*(1.0-omega0)**(-3./2) call tofzl(z,omega0,lambda,tauz,dum) call lingrowth(z,omega0,lambda,omega,b,hubble) deltac=(3./2)*b0*(1.0+(tomega/tauz)**(2./3)) return end if c c-- Omega+Lambda = 1 if (omega0+lambda.eq.1.0) then call lingrowth(z,omega0,lambda,omega,b,hubble) dcz = (3./20)*(12*pi)**(2./3)*(omega**0.0055) deltac=dcz/b return end if c stop 'DELTACRIT: Unknown option' c end c cccccccccccccccccccccccc function func2(x) c c-- Auxiliary function double precision x,func2 if (x.lt.1.e-5) then func2 = 0. else func2 = (x/(x**3+2))**1.5 endif return end c ccccccc subroutine simp(func,a,b,s) c c-- This subroutine from Numerical Recipes. c implicit none double precision func,a,b,eps,ost,os,s,st integer jmax,j external func ost = -1.e-30 os = -1.e-30 eps = 1.e-5 jmax = 40 do j = 1,jmax c write(*,*) a,b,s,j, ost call trap(func,a,b,st,j) c write(*,*) j, st, 'simp' s = (4.*st - ost)/3. if (abs(s-os) .lt. eps*abs(os)) goto 3 os = s ost = st enddo pause 'Too many steps.' 3 end c c cccccccccccccccccccccccccccccccccccc subroutine trap(func,a,b,s,n) c c-- This subroutine from Numerical Recipes. c implicit none double precision func,a,b,s,del,sum,x integer j,n,it,tnm external func if (n .eq. 1) then s = 0.5*(b-a)*(func(a)+func(b)) it = 1 else tnm = it del = (b-a)/tnm x = a + 0.5*del sum = 0.0 do j = 1,it sum = sum + func(x) x = x + del enddo s = 0.5*(s+(b-a)*sum/tnm) it = 2*it endif end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Quadrature using fifth order Runge-Kutta with adaptive step size. c Based on Press et al, Numerical Recipes in C, 2nd ed, pp 719-722. c c Runge-Kutta driver with adaptive stepsize control. Integrate starting c value y from a to b with accuracy eps, storing intermediate results in c global variables. dxinit should be set as a guessed first stepsize. c c Pass a second parameter to FUNC in fparm. c c by M.A.K. Gross c c double precision function INTEGRATE(FUNC,fp,np,a,b,dxinit,eps) implicit none integer np double precision a, b, eps, dxinit, FUNC, fp(np) external FUNC integer maxsteps parameter(maxsteps=100000000) double precision 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) c c yscale is the scaling used to monitor accuracy. This general-purpose c choice can be modified if need be. c 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. 1 dx = b - x call RUNGE5VAR(y,dydx,x,dx,eps,yscale,dxnext,FUNC,fp,np) dx = dxnext end do if (Nstep.ge.maxsteps) 1 write (*,*) 'WARNING: failed to converge in INTEGRATE.' INTEGRATE = y return end 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 SUBROUTINE RUNGE5VAR(y,dydx,x,htry,eps,yscale,hnext,DERIVS, 1 fp,np) implicit none integer np double precision fp(np) double precision eps,hnext,htry,x,dydx,y,yscale,DERIVS external DERIVS double precision errmax,h,hold,htemp,xnew,yerr,ytemp double precision 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: ', 1 'Stepsize underflow in RUNGE5VAR().' h = hold errmax = 0.d0 end if end if end do c c Step succeeded. Compute estimated size of next step. c if (errmax.gt.errcon) then hnext = safety*h*(errmax**pgrow) else hnext = 5.d0 * h ! No more than factor of 5 increase. end if x = x + h y = ytemp return end 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 double precision h,x,dydx,y,yerr,yout,DERIVS,fp(np) external DERIVS double precision ak3, ak4, ak5 ,ak6 double precision a2,a3,a4,a5,a6 double precision 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 c c ---------------- FUNCTION erfc(x) c ---------------- c Taken from Numerical Recipes. c REAL erfc,x c USES gammp,gammq REAL gammp,gammq external gammp,gammq if(x.lt.0.)then erfc=1.+gammp(.5,x**2) else erfc=gammq(.5,x**2) endif return END c c ------------------- FUNCTION gammln(xx) c ------------------- REAL gammln,xx INTEGER j DOUBLE PRECISION 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 11 j=1,6 y=y+1.d0 ser=ser+cof(j)/y 11 continue gammln=tmp+log(stp*ser/x) return END c c ------------------- FUNCTION gammp(a,x) c ------------------- REAL a,gammp,x c USES gcf,gser REAL gammcf,gamser,gln if(x.lt.0..or.a.le.0.)pause 'bad arguments in gammp' if(x.lt.a+1.)then call gser(gamser,a,x,gln) gammp=gamser else call gcf(gammcf,a,x,gln) gammp=1.-gammcf endif return END c c ------------------- FUNCTION gammq(a,x) c ------------------- REAL a,gammq,x c USES gcf,gser REAL gammcf,gamser,gln if(x.lt.0..or.a.le.0.)pause 'bad arguments in gammq' if(x.lt.a+1.)then call gser(gamser,a,x,gln) gammq=1.-gamser else call gcf(gammcf,a,x,gln) gammq=gammcf endif return END c c ------------------------------ SUBROUTINE gcf(gammcf,a,x,gln) c ------------------------------ INTEGER ITMAX REAL a,gammcf,gln,x,EPS,FPMIN PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30) c USES gammln INTEGER i REAL an,b,c,d,del,h,gammln external gammln gln=gammln(a) b=x+1.-a c=1./FPMIN d=1./b h=d do 11 i=1,ITMAX an=-i*(i-a) b=b+2. d=an*d+b if(abs(d).lt.FPMIN)d=FPMIN c=b+an/c if(abs(c).lt.FPMIN)c=FPMIN d=1./d del=d*c h=h*del if(abs(del-1.).lt.EPS)goto 1 11 continue pause 'a too large, ITMAX too small in gcf' 1 gammcf=exp(-x+a*log(x)-gln)*h return END c c ------------------------------- SUBROUTINE gser(gamser,a,x,gln) c ------------------------------- INTEGER ITMAX REAL a,gamser,gln,x,EPS PARAMETER (ITMAX=100,EPS=3.e-7) c USES gammln INTEGER n REAL ap,del,sum,gammln external gammln gln=gammln(a) if(x.le.0.)then if(x.lt.0.)pause 'x < 0 in gser' gamser=0. return endif ap=a sum=1./a del=sum do 11 n=1,ITMAX ap=ap+1. del=del*x/ap sum=sum+del if(abs(del).lt.abs(sum)*EPS)goto 1 11 continue pause 'a too large, ITMAX too small in gser' 1 gamser=sum*exp(-x+a*log(x)-gln) return END C (C) Copr. 1986-92 Numerical Recipes Software .-35?421.1-9. c c c -------------------------- REAL FUNCTION gasdev(idum) c -------------------------- INTEGER idum CU USES randd ! changed to use our random numer generator INTEGER iset REAL fac,gset,rsq,v1,v2,randd SAVE iset,gset DATA iset/0/ if (iset.eq.0) then 1 v1=2.*randd(idum)-1. v2=2.*randd(idum)-1. rsq=v1**2+v2**2 if(rsq.ge.1..or.rsq.eq.0.)goto 1 fac=sqrt(-2.*log(rsq)/rsq) gset=v1*fac gasdev=v2*fac iset=1 else gasdev=gset iset=0 endif return END C (C) Copr. 1986-92 Numerical Recipes Software .-35?421.1-9. c c ------------------------------------- FUNCTION zbrent(fp,np,func,x1,x2,tol) c ------------------------------------- INTEGER ITMAX REAL zbrent,tol,x1,x2,EPS integer np real fp(np) real*8 func EXTERNAL func PARAMETER (ITMAX=100,EPS=3.e-8) INTEGER iter REAL a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm a=x1 b=x2 fa=func(fp,a) fb=func(fp,b) if((fa.gt.0..and.fb.gt.0.).or.(fa.lt.0..and.fb.lt.0.)) then write(*,*) 'root must be bracketed for zbrent' zbrent = 0. return endif c=b fc=fb do 11 iter=1,ITMAX if((fb.gt.0..and.fc.gt.0.).or.(fb.lt.0..and.fc.lt.0.))then c=a fc=fa d=b-a e=d endif if(abs(fc).lt.abs(fb)) then a=b b=c c=a fa=fb fb=fc fc=fa endif tol1=2.*EPS*abs(b)+0.5*tol xm=.5*(c-b) if(abs(xm).le.tol1 .or. fb.eq.0.)then zbrent=b return endif if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then s=fb/fa if(a.eq.c) then p=2.*xm*s q=1.-s else q=fa/fc r=fb/fc p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.)) q=(q-1.)*(r-1.)*(s-1.) endif if(p.gt.0.) q=-q p=abs(p) if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then e=d d=p/q else d=xm e=d endif else d=xm e=d endif a=b fa=fb if(abs(d) .gt. tol1) then b=b+d else b=b+sign(tol1,xm) endif fb=func(fp,b) 11 continue pause 'zbrent exceeding maximum iterations' zbrent=b return END c c ------------------------------------- FUNCTION zbrent2(fp,np,func,x1,x2,tol) c ------------------------------------- INTEGER ITMAX REAL zbrent2,tol,x1,x2,EPS integer np real fp(np) real*8 func EXTERNAL func PARAMETER (ITMAX=100,EPS=3.e-8) INTEGER iter REAL a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm a=x1 b=x2 fa=func(fp,a) fb=func(fp,b) if((fa.gt.0..and.fb.gt.0.).or.(fa.lt.0..and.fb.lt.0.))pause *'root must be bracketed for zbrent2' c=b fc=fb do 11 iter=1,ITMAX if((fb.gt.0..and.fc.gt.0.).or.(fb.lt.0..and.fc.lt.0.))then c=a fc=fa d=b-a e=d endif if(abs(fc).lt.abs(fb)) then a=b b=c c=a fa=fb fb=fc fc=fa endif tol1=2.*EPS*abs(b)+0.5*tol xm=.5*(c-b) if(abs(xm).le.tol1 .or. fb.eq.0.)then zbrent2=b return endif if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then s=fb/fa if(a.eq.c) then p=2.*xm*s q=1.-s else q=fa/fc r=fb/fc p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.)) q=(q-1.)*(r-1.)*(s-1.) endif if(p.gt.0.) q=-q p=abs(p) if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then e=d d=p/q else d=xm e=d endif else d=xm e=d endif a=b fa=fb if(abs(d) .gt. tol1) then b=b+d else b=b+sign(tol1,xm) endif fb=func(fp,b) 11 continue pause 'zbrent2 exceeding maximum iterations' zbrent2=b return END C (C) Copr. 1986-92 Numerical Recipes Software .-35?421.1-9. C********************************************************************* c c ---------------------- real function RANDd(M) c ---------------------- c c random number generator c c to be initialized with a positive integer M, which must c never be changed by user after that c 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 RETURN END