ccc cosmodist.f Time-stamp: <2004-06-01 08:21:12 hamana> implicit real (a-h,o-z) implicit integer (i-n) common /omegas/ om,ov,ok common /curv/ sqok,iok external dlda external dtda parameter (ch0=2997.9) parameter (h0=9.78437) pi=4.0*atan(1.0) deg=pi/180.0 amin=pi/(180.0*60.0) asec=pi/(180.0*60.0*60.0) write(*,*)'### cosmodist.f' write(*,*)'### version-1 [20/06/2003]' write(*,*)'###' write(*,*)'### Use this at your own risk' write(*,*)'### Bug report to Takashi Hamana' write(*,*)'###' write(*,*)'### Hubble parameter is defined by' write(*,*)'### H_0=100h[km/sec/Mpc]' write(*,*)'#########################################' 11 write(*,*)'Cosmology [1], [2] [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.2) then om=1.0 ov=0.0 ok=0.0 sqok=0.0 iok=0 elseif (icosmo.eq.1) then om=0.3 ov=0.7 ok=0.0 sqok=0.0 iok=0 elseif (icosmo.eq.3) then write(*,*)'Omega_m ?' read(*,*)om write(*,*)'Lambda ?' read(*,*)ov ok=om+ov-1.0 sqok=sqrt(abs(ok)) iok=0 if (ok.lt.-1.0e-3) iok=-1 if (ok.gt.1.0e-3) iok=1 else goto 11 endif 22 write(*,*)'Redshift ?' read(*,*)z if (z.le.0.0) then write(*,*)'redshift should be larger than zero' goto 22 endif a=1.0/(1.0+z) !scale factor a0=1.0 call qromb(dlda,a,a0,ss) comov=ch0*ss cangd=ch0*angd(ss) ddeg=cangd*tan(deg) dmin=cangd*tan(amin) dsec=cangd*tan(asec)*1000.0 puti=1.0e-9 call qromb(dtda,puti,a0,ss) t0=h0*ss call qromb(dtda,puti,a,ss) t=h0*ss write(*,*)'######################################################' write(*,100)' Radial distance =',comov,'[comoving Mph/h]' write(*,100)'Angular diameter distance =',cangd,'[comoving Mph/h]' write(*,100)' 1 degree corresponds to ',ddeg,'[comoving Mpc/h]' write(*,100)' 1 arcmin corresponds to ',dmin,'[comoving Mpc/h]' write(*,100)' 1 arcsec corresponds to ',dsec,'[comoving Mpc/h]' write(*,101)' Age t_0 =',t0,'[Gyr/h]' write(*,101)' Age t(z) =',t,'[Gyr/h]' 100 format(a27,1pe12.5,a16) 101 format(a27,1pe12.5,a7) end !!! comoving distance -> angular diameter distance function angd(x) implicit real (a-h,o-z) implicit integer (i-n) common /omegas/ om,ov,ok common /curv/ sqok,iok if (iok.eq.0) then angd=x elseif (iok.eq.-1) then angd=sinh(sqok*x)/sqok else angd=sin(sqok*x)/sqok endif return end !!! dt/da=1.0/sqrt(a/om-ok+a2*ov) function dtda(a) implicit real (a-h,o-z) implicit integer (i-n) common /omegas/ om,ov,ok a2=a*a dtda=1.0/sqrt(om/a-ok+a2*ov) return end !!! dl/da=1.0/sqrt(a*om-a2*ok+a2*a2*ov) function dlda(a) implicit real (a-h,o-z) implicit integer (i-n) common /omegas/ om,ov,ok a2=a*a dlda=1.0/sqrt(a*om-a2*ok+a2*a2*ov) return end 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.