C--------------------------------------------------- C Anatoly Klypin Jauary 2000 C Formulas for LCDM model: age, angular-distance C C C----------------------------------------------------- C Om = Omega_0 Omc = Omega_cdm_present C Ocurv = 0 C hsmall = H_0/100km/s/Mpc = the Hubble constant C AEXPN = expansion parameter = 1/(1+z), z =redshift C----------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / Omegas/ Om,hsmall REAL*8 INTG,INPUT EXTERNAL INTG C Constants Om =INPUT(' Enter Omega_0 = matter at z=0 =>') hsmall =INPUT(' Enter hubble constant h =>') H =100.*hsmall ! Hubble constant open(50,file='Angles.dat', access='APPEND') DO a=0.1,1.,0.1 CALL AGE(t0,GrowthDen,GrowthVel,a) ENDDO a =1. CALL AGE(t0,GrowthDen_0,GrowthVel_0,a) z = INPUT(' Enter redshift =>') c angle =INPUT(' Enter angle (arcsec) =>') D =INPUT(' Enter proper radius (Mpc/h) =>') CALL DISTANCEANG(z,angle,D) c CALL ANGDISTANCE(z,angle,D) write (*,*) ' Proper distance (Mpc/h) =',D, & ' angle (arcsec)=',angle,' z=',z write (50,100) D,angle,z 100 format(' Proper distance (Mpc/h) =',f8.4, & ' angle (arcsec)=',f7.1,' z=',f7.3) STOP END C------------------------------------------------- C C angular-distance relation : C for given proper linear size (h=1) find angle (arcsecs) C at redshift z SUBROUTINE DISTANCEANG(z,angle,D) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (one =1.d-0) COMMON / Omegas/ Om,hsmall Real*8 INTG External Fangle if(z.le.0.d+0)Then write (*,*) ' wrong redshift z=',z return EndIf angle= D*206265./( 3000./(1.+z)*INTG(Fangle,one,one+z)) Return End C------------------------------------------------- C C angular-distance relation : C for given angle (arcsecs) find proper linear size (h=1) C at redshift z SUBROUTINE ANGDISTANCE(z,angle,D) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (one =1.d-0) COMMON / Omegas/ Om,hsmall Real*8 INTG External Fangle if(z.le.0.d+0)Then write (*,*) ' wrong redshift z=',z return EndIf D =angle/206265.* 3000./(1.+z)*INTG(Fangle,one,one+z) 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 / Omegas/ Om,hsmall Common /OmegRat/ OmLOm0,OmcOm0 Real*8 INTG External Hnorm,Hage,Hgrow Oml =Max(Abs(1.-Om),zero) Ocurv =0. 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 write (*,20) t0,z,a,Hubble,GrowthDen,GrowthVel write (1,20) t0,z,a,Hubble,GrowthDen,GrowthVel 20 format(' Age=',f8.3,' z=',f6.2,' a=',f6.3, . ' Hubble/100=',f8.3,' GrowthRateDen=',f7.4, . ' GrowthRateVelocity=',f7.4) Return End C-------------------------------------------- Simpson integration REAL*8 FUNCTION INTG(FUNC,A,B) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (EPS=5.0d-4, 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---------------------------------- 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 Fangle(x) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / Omegas/ Om,hsmall Fangle =1./sqrt( (1.-Om)+Om*x**3 ) Return End