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
