      program shearap

c The code solves the Sachs equations
c in a similar fashion as the Dyer-Roeder approximation.
c It includes the smoothness parameter c and the Weyl shearing (c, if c=0 then the Dyer-Roeder approximation is recovered)

c INPUT:  param.ini /cosmological parameters/
c OUTPUT: disred /z,D_A/

c author: Krzysztof Bolejko, The University of Sydney, bolejko@gmail.com
c last modification: 11 March 2013.

c The code can be used freely and/or modified.
c However, do not use the code if it could be harmful or endanger others.
c If you use this code, please cite: 


      implicit none
      integer I,n
      character ch
      double precision mu,lu,tu,gcons,hp,cs,kap,kapc2,pi
      double precision Om,Ho,gkr,Ol,Ok,lb,Hz,zlim
      double precision rho,l,dl,c,al,a,t
      double precision ks1,ks2,ks3,ks4,si,vs,s
      double precision kz1,kz2,kz3,kz4,zi,vz,z,zz
      double precision kd1,kd2,kd3,kd4,di,vd,d
      double precision kv1,kv2,kv3,kv4,vi,v

c Input
     	open(11,file='param.ini')
c Output
     	open(12,file='disred')
c reading the parameter file:
     	rewind(11)
	read(11,*) ch
	read(11,*) hp
	read(11,*) ch
	read(11,*) Om
	read(11,*) ch
	read(11,*) Ol
	read(11,*) ch
	read(11,*) al
	read(11,*) ch
	read(11,*) c
	read(11,*) ch
	read(11,*) zlim
c physical constants
	Ok = 1d0-Om-Ol
	mu=1.989d45
	lu=3.085678d19
	tu=31556926*1d6
	pi=4d0*datan(1d0)
	gcons= 6.67259d-11*((mu*(tu**2))/(lu**3))
	cs=299792458*(tu/lu)
	kap=8d0*pi*gcons*(1d0/(cs**4))
	kapc2=8d0*pi*gcons*(1d0/(cs**2))
	Ho=(tu/(lu))*hp
	gkr=3d0*(((Ho)**2)/(8d0*pi*gcons))
	gkr=gkr*kapc2*Om
	lb=3d0*Ol*(((Ho)**2)/(cs*cs))	
	Ho = Ho/cs
        n = 100000
        dl = -5d1
        
c initial data
      zi = 0d0
      si = 0d0
      di = 0d0
      vi = -1d0
      v = vi
      s = si
      d = di
      z = zi
      t = 0d0



      do I=1,n
      l = I*dl

      zz = 1d0+z
      Hz = Ho*dsqrt(Om*(zz**3)+Ok*(zz**2)+Ol)
      rho = gkr*(zz**3)
      vz = -1d0*zz*zz*Hz
      vs = -2d0*s*t + c*5d-1*rho*zz*zz*(al-1d0)
      a = -d*s*s-5d-1*rho*zz*zz*al*d
      kz1 = vz*dl
      kv1 = a*dl
      ks1 = vs*dl
      kd1 = v*dl
      s= si+5d-1*ks1
      v= vi+5d-1*kv1
      z= zi+5d-1*kz1
      d= di+5d-1*kd1

      zz = 1d0+z
      Hz = Ho*dsqrt(Om*(zz**3)+Ok*(zz**2)+Ol)
      rho = gkr*(zz**3)
      t = v/d
      vz = -1d0*zz*zz*Hz
      vs = -2d0*s*t + c*5d-1*rho*zz*zz*(al-1d0)
      a = -d*s*s-5d-1*rho*zz*zz*al*d
      kz2 = vz*dl
      kv2 = a*dl
      ks2 = vs*dl
      kd2 = v*dl
      s= si+5d-1*ks2
      v= vi+5d-1*kv2
      z= zi+5d-1*kz2
      d= di+5d-1*kd2

      zz = 1d0+z
      Hz = Ho*dsqrt(Om*(zz**3)+Ok*(zz**2)+Ol)
      rho = gkr*(zz**3)
      t = v/d
      vz = -1d0*zz*zz*Hz
      vs = -2d0*s*t + c*5d-1*rho*zz*zz*(al-1d0)
      a = -d*s*s-5d-1*rho*zz*zz*al*d
      kz3 = vz*dl
      kv3 = a*dl
      ks3 = vs*dl
      kd3 = v*dl
      s= si+ks3
      v= vi+kv3
      z= zi+kz3
      d= di+kd3

      zz = 1d0+z
      Hz = Ho*dsqrt(Om*(zz**3)+Ok*(zz**2)+Ol)
      rho = gkr*(zz**3)
      t = v/d
      vz = -1d0*zz*zz*Hz
      vs = -2d0*s*t + c*5d-1*rho*zz*zz*(al-1d0)
      a = -d*s*s-5d-1*rho*zz*zz*al*d
      kz4 = vz*dl
      kv4 = a*dl
      ks4 = vs*dl
      kd4 = v*dl

      s= si+(ks1+2d0*(ks2+ks3)+ks4)/6d0
      v= vi+(kv1+2d0*(kv2+kv3)+kv4)/6d0
      z= zi+(kz1+2d0*(kz2+kz3)+kz4)/6d0
      d= di+(kd1+2d0*(kd2+kd3)+kd4)/6d0

      si = s
      zi = z
      vi = v
      di = d
      t = v/d

	if(z.ge.1000d0 .or. z.ge.zlim) goto 137
	if(d.le.0d0) then
		print *, 'approximation breaks down /caustics/'
		stop
	endif

c output: z, D_A in Gpc
      if(mod(I,5).eq.0) write(12,*) z,d*1d-6
      enddo

 137	continue
      end
