<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<!-- saved from url=(0044)http://www.atlasoftheuniverse.com/cosmodis.c -->
<HTML><HEAD>
<META content="text/html; charset=windows-1252" http-equiv=Content-Type>
<META name=GENERATOR content="MSHTML 8.00.6001.18904"></HEAD>
<BODY><PRE>/*
cosmodis - A Cosmological Distances Program - version 1.1
by Richard Powell - http://www.atlasoftheuniverse.com/
This is a simple piece of code to provide comoving, angular diameter,
luminosity, and light travel distances for any given redshift. The
Hubble constant (H0), Omega_matter (OM) and Omega_lambda (OL) are defined
within the body of the program and can be adjusted to your favourite
values.
For a summary of the formulae used in this program, see:
David Hogg, Distance Measures in Cosmology, (2000), astro-ph/9905116.
http://arxiv.org/abs/astro-ph/9905116
The standard command for compiling this program in Linux is:
cc cosmodis.c -lm -ocosmodis
This program is in the Public Domain. There is no copyright.
*/
#include <stdio.h>
#include <math.h>
#define c 299792.458
int main()
{
float H0 = 70; // Hubble constant (km/s/Mpc) - adjust according to taste
float OM = 0.27; // Omega(matter) - adjust according to taste
float OL = 0.73; // Omega(lambda) - adjust according to taste
float OR = 0.42/(H0*H0); // Omega(radiation) - this is the usual textbook value
long i;
long n=10000; // Number of steps in the integral
float OK = 1-OM-OR-OL; // Omega(k) defined as 1-OM-OR-OL
float HD = 3.2616*c/H0/1000; // Hubble distance (billions of light years). See section 2 of Hogg
float z, a, adot; // Redshift "z", Scale Factor "a", and its derivative "adot"
float DC, DCC=0, DT, DTT=0, DA, DL, DM;
float age, size; // The age and size of the universe
printf("Enter Redshift: "); // Ask for a redshift z
scanf("%f",&z);
for(i=n; i>=1; i--) { // This loop is the numerical integration
a = (i-0.5)/n; // Steadily decrease the scale factor
// Comoving formula (See section 4 of Hogg, but I've added a radiation term too):
adot = a*sqrt(OM*pow(1/a,3)+OK*pow(1/a,2)+OL+OR*pow(1/a,4)); // Note that "a" is equivalent to 1/(1+z)
DCC = DCC + 1/(a*adot)/n; // Running total of the comoving distance
DTT = DTT + 1/adot/n; // Running total of the light travel time (see section 10 of Hogg)
if (a>=1/(1+z)) { // Collect DC and DT until the correct scale is reached
DC = DCC; // Comoving distance DC
DT = DTT; // Light travel time DT
}
}
// Transverse comoving distance DM from section 5 of Hogg:
if (OK>0.0001) DM=(1/sqrt(OK))*sinh(sqrt(OK)*DC);
else if (OK<-0.0001) DM=(1/sqrt(fabs(OK)))*sin(sqrt(fabs(OK))*DC);
else DM=DC;
age = HD*DTT; // Age of the universe (billions of years)
size = HD*DCC; // Comoving radius of the observable universe
DC = HD*DC; // Comoving distance
DA = HD*DM/(1+z); // Angular diameter distance (section 6 of Hogg)
DL = HD*DM*(1+z); // Luminosity distance (section 7 of Hogg)
DT = HD*DT; // Light travel distance
printf("-------------------------------------------------------------------\n");
printf("For Redshift %.3f, (Ho=%.1fkm/s/Mpc, Omega_M=%.2f, Omega_L=%.2f):\n",z,H0,OM,OL);
printf("-------------------------------------------------------------------\n");
printf("* Age of the universe now = %.3f Gyr\n",age);
printf("* Age of the universe then = %.6f Gyr\n",age-DT);
printf("* Comoving horizon of the universe now = %.3f Gyr\n",size);
printf("* Comoving horizon of the universe then = %.3f Gyr\n",size/(1+z));
printf("* Comoving distance of the source now = %.3f Gly\n",DC);
printf("* Comoving distance of the source then = %.3f Gly\n",DC/(1+z)); // In a flat universe, this is equal to DA
printf("* Angular Diameter distance = %.3f Gly\n",DA);
printf("* Luminosity distance = %.3f Gly\n",DL);
printf("* Light Travel Time Distance = %.3f Gly\n",DT);
printf("-------------------------------------------------------------------\n");
}
</PRE></BODY></HTML>