- #1
PetradR
- 10
- 1
- TL;DR Summary
- Hi, I am new to the forum &New to Fortran: HELP!
Trying to put together a program to calculate on a single graph, a plot the following:
*the luminosity distance,
*the present proper distance,
*the angular diameter distance,
*the relativistic Hubble law distance estimate and
*the non-relativistic Hubble law distance for values of z between 0 and 4. Expressed these in units of c/H_0
What I have done so far is below: still cant get it to just give me just me 5 rows of values.
Fortran:
PROGRAM Test
IMPLICIT NONE
REAL, PARAMETER :: Omega_m0 = 0.27, Omega_rel0 = 8.24E-5, Omega_L0 = 0.73, Omega_0 = 1.02
REAL :: z_max, z, dz
REAL :: I, g_z, g_zplus, S
REAL :: d_p0, d_L, d_A, d_nr, d_r
INTEGER :: n = 100000
OPEN(UNIT = 10)
z_max = 0
DO WHILE (z_max <= 4.001)
dz = z_max/n
!integrand of Eq. (29.168) z=0 g_z = 1/SQRT(Omega_m0*(1 + z)**3 + Omega_rel0*(1 + z)**4 + Omega_L0 + (1 - Omega_0)*(1 + z)**2)
I = 0
DO WHILE (z < z_max)
z = z + dz
g_zplus = 1/SQRT(Omega_m0*(1 + z)**3 + Omega_rel0*(1 + z)**4 + Omega_L0 + (1 - Omega_0)*(1 + z)**2)
I = (I + g_z)*(dz/2)
END DO
S = SIN(I*SQRT(Omega_0 - 1))/sqrt(Omega_0 - 1) !Eq. (29.174) for Omega_0 > 1
d_L = S*(1 + z_max) !Eq. (29.184) in units of c/H_0
d_p0 = I !Eq. (29.169) in units of c/H_0
d_A = d_L/(1 + z_max)**2 !Eq. (29.192) in units of c/H_0
d_r = ((z + 1)**2 - 1)/((z + 1)**2 + 1) !Eq. (27.7) in units of c/H_0
d_nr = z !Eq. (27.8) in units of c/H_0
WRITE(10,'(7E13.6)') z, d_L, d_p0, d_A, d_r, d_nr, ABS((d_L - d_r)/d_L)
z_max = z_max + 0.01
END DO
PAUSE
STOP
END
Last edited by a moderator: