ccc varkap.f Time-stamp: <2004-06-03 12:06:01 hamana> implicit real (a-h,o-z) implicit integer (i-n) common /cosmo/ om,ov,crv,scrv common /dpower/rn,gamfin,sig8,anrm external dw parameter (ismin=-3000) parameter (ismax=8000) parameter (dels=0.001) dimension s(ismin:ismax),dps(ismin:ismax),pskap(ismin:ismax) dimension dpsnl(ismin:ismax),pskapnl(ismin:ismax) parameter (ia=100) dimension asf(0:ia),cd(0:ia),d(0:ia),dod(0:ia) parameter (nrkmin=-30,nrkmax=40) parameter (nrk=nrkmax-nrkmin+1) parameter (delrk=0.1) dimension rk(nrk),pnl(nrk),pnl2(nrk) parameter (nphimin=-10,nphimax=20,delphi=0.1) pi=4.0*atan(1.0) twopi=2.0*pi tps=2.0*pi**2 am2deg=pi/(180.0*60.0) cln2log=log(10.0) cv=3000.0 igrfrag=-1 write(*,*)'### varkap.f' write(*,*)'### version-1 [3/6/2004]' write(*,*)'###' write(*,*)'### Use this at your own risk' write(*,*)'### Bug report to hamanatk [at-mark] cc.nao.ac.jp' write(*,*)'###' write(*,*)'### This uses BBKS CDM P(k) with ' write(*,*)'### Sugiyama (1995) Gamma factor and' write(*,*)'### Smith et al (2003) nonlinear fitting function' write(*,*)'###' write(*,*)'### This will output a data file varkap.dat' write(*,*)'###' write(*,*)'### Hubble parameter is defined by' write(*,*)'### H_0=100h[km/sec/Mpc]' write(*,*)'#########################################' 11 write(*,*)'Cosmology 1, 2 or 3 ?' write(*,*)'1: Lambda-flat (Omega_m=0.3, Lambda=0.7)' write(*,*)'2: Eds (Omega_m=1.0, Lambda=0.0)' write(*,*)'3: something else' read(*,*)icosmo if (icosmo.eq.1) then om=0.3 ov=0.7 elseif (icosmo.eq.2) then om=1.0 ov=0.0 elseif (icosmo.eq.3) then write(*,*)'Omega_m ?' read(*,*)om write(*,*)'Lambda ?' read(*,*)ov else goto 11 endif crv=om+ov-1.0 scrv=sqrt(abs(crv)) kcv=0 if (crv.gt.0.001) kcv=1 if (crv.lt.-0.001) kcv=-1 12 write(*,*)'Hubble parameter 1, 2 or 3 ?' write(*,*)'1: h=0.7' write(*,*)'2: h=0.5' write(*,*)'3: something else' read(*,*)icosmo if (icosmo.eq.1) then hubp=0.7 elseif (icosmo.eq.2) then hubp=0.5 elseif (icosmo.eq.3) then write(*,*)'h ?' read(*,*)hubp else goto 12 endif 13 write(*,*)'Baryon density 1, 2 or 3 ?' write(*,*)'1: Omega_b=0.04' write(*,*)'2: Omega_b=0.0' write(*,*)'3: something else' read(*,*)icosmo if (icosmo.eq.1) then ob=0.04 elseif (icosmo.eq.2) then ob=0.0 elseif (icosmo.eq.3) then write(*,*)'Omega_b ?' read(*,*)ob else goto 13 endif ccc gamma factor by Sugiyama (1995) gamfin=exp(ob+sqrt(2.0*hubp)*ob/om)/(om*hubp*cv) 14 write(*,*)'sigma_8 1, 2 or 3 ?' write(*,*)'1: sigma_8=0.9' write(*,*)'2: sigma_8=0.5' write(*,*)'3: something else' read(*,*)icosmo if (icosmo.eq.1) then sig8=0.9 elseif (icosmo.eq.2) then sig8=0.5 elseif (icosmo.eq.3) then write(*,*)'sigma_8 ?' read(*,*)sig8 else goto 14 endif ccccc normalization of power spectrum call normPk coef=9.0/4.0*om*om 15 write(*,*)'source redshift 1, 2 or 3 ?' write(*,*)'1: zs=1.0' write(*,*)'2: zs=2.0' write(*,*)'3: something else' read(*,*)icosmo if (icosmo.eq.1) then zmax=1.0 elseif (icosmo.eq.2) then zmax=2.0 elseif (icosmo.eq.3) then write(*,*)'zs ?' read(*,*)zmax else goto 15 endif ccc source redshift amin=1.0/(1.0+zmax) dela=(1.0-amin)/ia dela2=0.5*dela ccc complute a-w a-r relation asf(0)=1.0 cd(0)=0.0 d(0)=0.0 do j=1,ia asf(j)=1.0-j*dela call qromb(dw,asf(j),1.0,ss) cd(j)=ss d(j)=fk(ss,kcv) enddo wi=cd(ia) fks=1.0/d(ia) do j=0,ia wij=wi-cd(j) dod1=fk(wij,kcv)*fks dod(j)=dod1*dod1 enddo do i=ismin,ismax s(i)=10.0**(i*dels) pskap(i)=0.0 dps(i)=0.0 pskapnl(i)=0.0 dpsnl(i)=0.0 enddo do i=1,ia z=1.0/asf(i)-1.0 gr2=gr(asf(i),igrfrag)**2 weight=dod(i)/sqrt(e(asf(i))) weightnl=weight/asf(i)**2 call halofit_v3(om,ov,sig8,ob,hubp,z, $ nrkmin,nrkmax,nrk,delrk,rk,pnl) yp1=(pnl(2)-pnl(1))/(rk(2)-rk(1)) ypn=(pnl(nrk)-pnl(nrk-1))/(rk(nrk)-rk(nrk-1)) call spline(rk,pnl,nrk,yp1,ypn,pnl2) pow=(log(pnl(nrk))-log(pnl(nrk-1)))/ $ (log(rk(nrk))-log(rk(nrk-1))) do ks=ismin,ismax ak=s(ks)/d(i) !!! linear dps2=weight*gr2*plin(ak) pskap(ks)=pskap(ks)+ $ (dps(ks)+dps2)*dela2 dps(ks)=dps2 !!! non-linear x=ak/cv y=spfit(rk,pnl,pnl2,nrk,pow,x)*tps/ak**3 dps2=weightnl*y pskapnl(ks)=pskapnl(ks)+ $ (dpsnl(ks)+dps2)*dela2 dpsnl(ks)=dps2 enddo enddo do ks=ismin,ismax pskap(ks)=coef*pskap(ks) pskapnl(ks)=coef*pskapnl(ks) enddo open(10,file='varkap.dat',status='unknown') write(10,'(a31)')'# theta [arcmin]; Varlin; Varnl' do n=nphimin,nphimax phi=10.0**(delphi*n) phideg=phi*am2deg var=0.0 varnl=0.0 do i=ismin,ismax x=s(i)*phideg weight=(s(i)*bessj1(x)/x)**2 var=var+pskap(i)*weight varnl=varnl+pskapnl(i)*weight enddo var=2.0*var*dels*cln2log/pi varnl=2.0*varnl*dels*cln2log/pi write(10,'(3(1pe13.5))')phi,var,varnl enddo close(10) end !!! function spfit(rk,pnl,pnl2,nrk,pow,x) implicit real (a-h,o-z) implicit integer (i-n) dimension rk(nrk),pnl(nrk),pnl2(nrk) if (x.lt.rk(1)) then y=pnl(1)*(x/rk(1))**4 elseif (x.gt.rk(nrk)) then y=pnl(nrk)*(x/rk(nrk))**pow else call splint(rk,pnl,pnl2,nrk,x,y) endif spfit=y return end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine halofit_v3(om_m0,om_v0,sig80,om_b0,hubp,z, $ nrkmin,nrkmax,ndat,delrk,rkdim,pnldim) implicit none real rn,plin,pq,ph,pnl,rk,p_cdm,aexp,z real rknl,rneff,rncur,d1,d2,rad,sig,rknl_int_pow real om_m,om_v,om_b,h,p_index,gams,sig8,amp real om_m0,om_v0,omega_m,omega_v,grow,grow0,gg real rn_pd,rn_cdm,pnl_pd,rklin,f_pd real f1a,f2a,f3a,f4a,f1b,f2b,f3b,f4b,frac,f1,f2,f3,f4 real diff,xlogr1,xlogr2,rmid real om_b0,hubp real rkmin,rkmax integer i,j,k,ndat real sig80,delrk real rkdim(ndat),plindim(ndat),pnldim(ndat) integer nrkmin,nrkmax common/cospar/om_m,om_v,h,p_index,gams,sig8,amp sig8=sig80 gams=om_m0*hubp/exp(om_b0+sqrt(2.0*hubp)*om_b0/om_m0) aexp=1./(1.+z) ! expansion factor c calculate matter density, vacuum density at desired redshift om_m=omega_m(aexp,om_m0,om_v0) om_v=omega_v(aexp,om_m0,om_v0) c calculate the amplitude of the power spectrum at desired redshift c using linear growth factors (gg given by Caroll, Press & Turner 1992, ARAA, 30, 499) grow=gg(om_m,om_v) grow0=gg(om_m0,om_v0) amp=aexp*grow/grow0 c calculate nonlinear wavenumber (rknl), effective spectral index (rneff) and c curvature (rncur) of the power spectrum at the desired redshift, using method c described in Smith et al (2002). ! write(*,*) 'computing effective spectral quantities:' xlogr1=-2.0 xlogr2=3.5 10 rmid=(xlogr2+xlogr1)/2.0 rmid=10**rmid call wint(rmid,sig,d1,d2) diff=sig-1.0 if (abs(diff).le.0.001) then rknl=1./rmid rneff=-3-d1 rncur=-d2 elseif (diff.gt.0.001) then xlogr1=log10(rmid) goto 10 elseif (diff.lt.-0.001) then xlogr2=log10(rmid) goto 10 endif ! write(*,20) 'rknl [h/Mpc] =',rknl,'rneff=',rneff, 'rncur=',rncur ! 20 format(a14,f12.6,2x,a6,f12.6,2x,a6,f12.6) c now calculate power spectra for a logarithmic range of wavenumbers (rk) ! write(*,*) 'press return to compute halofit nonlinear power' ! read(*,*) rkmin=nrkmin*delrk rkmax=nrkmax*delrk j=0 do i=nrkmin,nrkmax j=j+1 rk=10.0**(i*delrk) c linear power spectrum !! Remeber => plin = k^3 * P(k) * constant c constant = 4*pi*V/(2*pi)^3 plin=amp*amp*p_cdm(rk,gams,sig8) c calculate nonlinear power according to halofit: pnl = pq + ph, c where pq represents the quasi-linear (halo-halo) power and c where ph is represents the self-correlation halo term. call halofit(rk,rneff,rncur,rknl,plin,pnl,pq,ph) ! halo fitting formula rkdim(j)=rk plindim(j)=plin pnldim(j)=pnl enddo return end c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c halo model nonlinear fitting formula as described in c Appendix C of Smith et al. (2002) subroutine halofit(rk,rn,rncur,rknl,plin,pnl,pq,ph) implicit none real gam,a,amod,b,c,xmu,xnu,alpha,beta,f1,f2,f3,f4 real rk,rn,plin,pnl,pq,ph real om_m,om_v,h,p_index,gams,sig8,amp real rknl,y,rncur,p_cdm real f1a,f2a,f3a,f4a,f1b,f2b,f3b,f4b,frac real yy common/cospar/om_m,om_v,h,p_index,gams,sig8,amp gam=0.86485+0.2989*rn+0.1631*rncur a=1.4861+1.83693*rn+1.67618*rn*rn+0.7940*rn*rn*rn+ &0.1670756*rn*rn*rn*rn-0.620695*rncur a=10**a b=10**(0.9463+0.9466*rn+0.3084*rn*rn-0.940*rncur) c=10**(-0.2807+0.6669*rn+0.3214*rn*rn-0.0793*rncur) xmu=10**(-3.54419+0.19086*rn) xnu=10**(0.95897+1.2857*rn) alpha=1.38848+0.3701*rn-0.1452*rn*rn beta=0.8291+0.9854*rn+0.3400*rn**2 if(abs(1-om_m).gt.0.01) then ! omega evolution f1a=om_m**(-0.0732) f2a=om_m**(-0.1423) f3a=om_m**(0.0725) f1b=om_m**(-0.0307) f2b=om_m**(-0.0585) f3b=om_m**(0.0743) frac=om_v/(1.-om_m) f1=frac*f1b + (1-frac)*f1a f2=frac*f2b + (1-frac)*f2a f3=frac*f3b + (1-frac)*f3a else f1=1.0 f2=1. f3=1. endif y=(rk/rknl) ph=a*y**(f1*3)/(1+b*y**(f2)+(f3*c*y)**(3-gam)) ph=ph/(1+xmu*y**(-1)+xnu*y**(-2)) yy=-y/4.0-y**2/8.0 if (yy.gt.-88.0) then pq=plin*(1+plin)**beta/(1+plin*alpha)*exp(-y/4.0-y**2/8.0) else pq=0.0 endif pnl=pq+ph return end c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c Bond & Efstathiou (1984) approximation to the linear CDM power spectrum !!! changed to BBKS96 in v2 function p_cdm(rk,gams,sig8) implicit none real p_cdm,rk,gams,sig8,p_index,rkeff,q,q8,tk,tk8 integer flag save flag,q8,tk8 if (flag.eq.0) then rkeff=0.172+0.011*log(gams/0.36)*log(gams/0.36) q8=1.e-20 + rkeff/gams tk8=log(1.0+2.34*q8)/((2.34*q8*(1.0+3.89*q8+(16.1*q8)**2+ $ (5.46*q8)**3+(6.71*q8)**4)**0.25)) flag=1 endif q=1.e-20 + rk/gams tk=log(1.0+2.34*q)/( $ (2.34*q* $ (1.0+3.89*q+(16.1*q)**2+(5.46*q)**3+(6.71*q)**4)**0.25)) p_cdm=sig8*sig8*((q/q8)**4)*tk*tk/(tk8*tk8) return end c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c The subroutine wint, finds the effective spectral quantities c rknl, rneff & rncur. This it does by calculating the radius of c the Gaussian filter at which the variance is unity = rknl. c rneff is defined as the first derivative of the variance, calculated c at the nonlinear wavenumber and similarly the rncur is the second c derivative at the nonlinear wavenumber. subroutine wint(r,sig,d1,d2) implicit none real sum1,sum2,sum3,t,y,x,w1,w2,w3 real rk, p_cdm, r, sig, d1,d2 real om_m,om_v,h,p_index,gams,sig8,amp integer i,nint common/cospar/om_m,om_v,h,p_index,gams,sig8,amp nint=3000 sum1=0. sum2=0. sum3=0. do i=1,nint t=(float(i)-0.5)/float(nint) y=-1.+1./t rk=y d2=amp*amp*p_cdm(rk,gams,sig8) x=y*r w1=exp(-x*x) w2=2*x*x*w1 w3=4*x*x*(1-x*x)*w1 sum1=sum1+w1*d2/y/t/t sum2=sum2+w2*d2/y/t/t sum3=sum3+w3*d2/y/t/t enddo sum1=sum1/float(nint) sum2=sum2/float(nint) sum3=sum3/float(nint) sig=sqrt(sum1) d1=-sum2/sum1 d2=-sum2*sum2/sum1/sum1 - sum3/sum1 return end c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c evolution of omega matter with expansion factor function omega_m(aa,om_m0,om_v0) implicit none real omega_m,omega_t,om_m0,om_v0,aa omega_t=1.0+(om_m0+om_v0-1.0)/(1-om_m0-om_v0+om_v0*aa*aa+om_m0/aa) omega_m=omega_t*om_m0/(om_m0+om_v0*aa*aa*aa) return end c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c evolution of omega lambda with expansion factor function omega_v(aa,om_m0,om_v0) implicit none real aa,omega_v,om_m0,om_v0,omega_t omega_t=1.0+(om_m0+om_v0-1.0)/(1-om_m0-om_v0+om_v0*aa*aa+om_m0/aa) omega_v=omega_t*om_v0/(om_v0+om_m0/aa/aa/aa) return end c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c growth factor for linear fluctuations function gg(om_m,om_v) implicit none real gg,om_m,om_v gg=2.5*om_m/(om_m**(4./7.)-om_v+(1.+om_m/2.)*(1.+om_v/70.)) return end c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ccc CDM power spectrum function plin(w) implicit real (a-h,o-z) implicit integer (i-n) common /dpower/rn,gamfin,sig8,anrm q=w*gamfin plin=anrm*trnf2(q)*w return end ccc CDM transfer function by BBKS 1986 function trnf2(q) implicit real (a-h,o-z) implicit integer (i-n) qq=q*q trnf2=0.182628*log(1.0+2.34*q)**2/qq $ /sqrt(1.0+3.89*q+259.21*qq $ +162.77*qq*q+2027.17*qq*qq) return end ccc linear growth factor function gr(a,igrfrag) implicit real (a-h,o-z) implicit integer (i-n) common /cosmo/ om,ov,crv,scrv save grnorm if (igrfrag.eq.-1) then et=1.0/e(1.0) ot=om*et vt=ov*et grnorm=0.4*(ot**(0.571428571)-vt+ $ (1.0+0.5*ot)*(1.0+vt*0.014285714))/ot igrfrag=1 endif et=1.0/e(a) ot=om*a*et vt=ov*et*a**4 graw=2.5*ot/(ot**0.571428571-vt+ $ (1.0+0.5*ot)*(1.0+vt*0.014285714)) gr=graw*grnorm return end ccc linear grwoth rate function grpd(a) implicit real (a-h,o-z) implicit integer (i-n) common /cosmo/ om,ov,crv,scrv et=1.0/e(a) ot=om*a*et vt=ov*et*a**4 graw=2.5*ot/(ot**0.571428571-vt+ $ (1.0+0.5*ot)*(1.0+vt*0.014285714)) grpd=graw return end ccc function dw(a) implicit real (a-h,o-z) implicit integer (i-n) dw=1.0/sqrt(e(a)) return end ccc function e(a) implicit real (a-h,o-z) implicit integer (i-n) common /cosmo/ om,ov,crv,scrv a2=a*a e=a*om-a2*crv+a2*a2*ov return end ccc function fk(w,kcv) implicit real (a-h,o-z) implicit integer (i-n) common /cosmo/ om,ov,crv,scrv if (kcv.eq.-1) then fk=sinh(scrv*w)/scrv elseif (kcv.eq.0) then fk=w else fk=sin(scrv*w)/scrv endif return end ccccccc normalization of power spectrum subroutine normPk implicit real (a-h,o-z) implicit integer (i-n) common /dpower/rn,gamfin,sig8,anrm external delta pi=4.0*atan(1.0) coe=4.5/pi/pi rn=8.0/3000.0 wmin=-10.0 wmax=10.0 call qromb(delta,wmin,wmax,seki) y=coe*seki anrm=sig8*sig8/y return end ccccccccc k k^2 t^2 W ccccccccccccccccccccccccccccccccccc function delta(x) implicit real (a-h,o-z) implicit integer (i-n) common /dpower/rn,gamfin,sig8,anrm w=exp(x) q=w*gamfin r=w*rn filt=(sin(r)-r*cos(r))/r**3 tranf=trnf2(q) delta=w**4*filt*filt*tranf return end cccccc integration pakege cccccccccccccccccccccccccccccccc SUBROUTINE qromb(func,a,b,ss) INTEGER JMAX,JMAXP,K,KM REAL a,b,func,ss,EPS EXTERNAL func PARAMETER (EPS=1.e-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1) CU USES polint,trapzd INTEGER j REAL dss,h(JMAXP),s(JMAXP) h(1)=1. do 11 j=1,JMAX call trapzd(func,a,b,s(j),j) if (j.ge.K) then call polint(h(j-KM),s(j-KM),K,0.,ss,dss) if (abs(dss).le.EPS*abs(ss)) return endif s(j+1)=s(j) h(j+1)=0.25*h(j) 11 continue pause 'too many steps in qromb' END C (C) Copr. 1986-92 Numerical Recipes Software 41m. SUBROUTINE trapzd(func,a,b,s,n) INTEGER n REAL a,b,s,func EXTERNAL func INTEGER it,j REAL del,sum,tnm,x if (n.eq.1) then s=0.5*(b-a)*(func(a)+func(b)) else it=2**(n-2) tnm=it del=(b-a)/tnm x=a+0.5*del sum=0. do 11 j=1,it sum=sum+func(x) x=x+del 11 continue s=0.5*(s+(b-a)*sum/tnm) endif return END C (C) Copr. 1986-92 Numerical Recipes Software 41m. SUBROUTINE polint(xa,ya,n,x,y,dy) INTEGER n,NMAX REAL dy,x,y,xa(n),ya(n) PARAMETER (NMAX=10) INTEGER i,m,ns REAL den,dif,dift,ho,hp,w,c(NMAX),d(NMAX) ns=1 dif=abs(x-xa(1)) do 11 i=1,n dift=abs(x-xa(i)) if (dift.lt.dif) then ns=i dif=dift endif c(i)=ya(i) d(i)=ya(i) 11 continue y=ya(ns) ns=ns-1 do 13 m=1,n-1 do 12 i=1,n-m ho=xa(i)-x hp=xa(i+m)-x w=c(i+1)-d(i) den=ho-hp if(den.eq.0.)pause 'failure in polint' den=w/den d(i)=hp*den c(i)=ho*den 12 continue if (2*ns.lt.n-m)then dy=c(ns+1) else dy=d(ns) ns=ns-1 endif y=y+dy 13 continue return END C (C) Copr. 1986-92 Numerical Recipes Software 41m. cccccc spline fitting ccccccccccccccccccccccccccccccccc SUBROUTINE splint(xa,ya,y2a,n,x,y) INTEGER n REAL x,y,xa(n),y2a(n),ya(n) INTEGER k,khi,klo REAL a,b,h klo=1 khi=n 1 if (khi-klo.gt.1) then k=(khi+klo)/2 if(xa(k).gt.x)then khi=k else klo=k endif goto 1 endif h=xa(khi)-xa(klo) if (h.eq.0.) pause 'bad xa input in splint' a=(xa(khi)-x)/h b=(x-xa(klo))/h y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h** *2)/6. return END C (C) Copr. 1986-92 Numerical Recipes Software 41m. SUBROUTINE spline(x,y,n,yp1,ypn,y2) INTEGER n,NMAX REAL yp1,ypn,x(n),y(n),y2(n) PARAMETER (NMAX=500) INTEGER i,k REAL p,qn,sig,un,u(NMAX) if (yp1.gt..99e30) then y2(1)=0. u(1)=0. else y2(1)=-0.5 u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) endif do 11 i=2,n-1 sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) p=sig*y2(i-1)+2. y2(i)=(sig-1.)/p u(i)=(6.*((y(i+1)-y(i))/(x(i+ *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig* *u(i-1))/p 11 continue if (ypn.gt..99e30) then qn=0. un=0. else qn=0.5 un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) endif y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) do 12 k=n-1,1,-1 y2(k)=y2(k)*y2(k+1)+u(k) 12 continue return END C (C) Copr. 1986-92 Numerical Recipes Software 41m. FUNCTION bessj1(x) REAL bessj1,x REAL ax,xx,z DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6, *s1,s2,s3,s4,s5,s6,y SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4, *s5,s6 DATA r1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0, *242396853.1d0,-2972611.439d0,15704.48260d0,-30.16036606d0/,s1,s2, *s3,s4,s5,s6/144725228442.d0,2300535178.d0,18583304.74d0, *99447.43394d0,376.9991397d0,1.d0/ DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4, *.2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/.04687499995d0, *-.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/ if(abs(x).lt.8.)then y=x**2 bessj1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+ *y*(s4+y*(s5+y*s6))))) else ax=abs(x) z=8./ax y=z**2 xx=ax-2.356194491 bessj1=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y* *p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))*sign(1.,x) endif return END C (C) Copr. 1986-92 Numerical Recipes Software 41m.