Hello
I need help in translating a C++ code into python..
Can you help please?
Please find attached the two codes.
The python result is wrong because I may have misread the C++ code
Thanks
Vickram
from mpmath import *
mp.dps=5
c = mpf('299792.458')
H0 = float(70) # Hubble constant (km/s/Mpc) - adjust according to taste
OM = float(0.27) # Omega(matter) - adjust according to taste
OL = float(0.73) # Omega(lambda) - adjust according to taste
OR = float(0.42/(H0*H0)) #Omega(radiation) - this is the usual textbook value
n=10000 # Number of steps in the integral
OK = float(1-OM-OR-OL) # Omega(k) defined as 1-OM-OR-OL
HD = float(3.2616*c/H0/1000) # Hubble distance (billions of light years). See
section 2 of Hogg
# Redshift "z", Scale Factor "a", and its derivative "adot"
DC =float(0)
DCC=float(0)
DT = float(0)
DTT=float(0)
a = float(0)
adot = float(0)
# The age and size of the universe
z = input("Enter Redshift: ") #Ask for a redshift z
for i in range(n,0,-1): # This loop is the numerical integration
a = (float(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 = float(a)*sqrt(OM*pow(1/float(a)
,3)+OK*pow(1/float(a),2)+OL+OR*pow(1/float(a),4))# Note that "a" is equivalent
to 1/(1+z)
DCC = float(DCC) + 1/(float(a)*float(adot))/long(n) # Running total of the
comoving distance
DTT = float(DTT) + 1/float(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)
elif OK<-0.0001:
DM=(1/sqrt(fabs(OK)))*sin(sqrt(fabs(OK))*DC)
else:
DM=DC
print a
print adot
print DT
print DTT
print DC
print DCC
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
print"-------------------------------------------------------------------"
print "For Redshift",z,", Ho= ",H0,"km/s/Mpc, Omega_M=",OM,", Omega_L=",OL
print "-------------------------------------------------------------------"
print "* Age of the universe now = ",age,"Gyr"
print "* Age of the universe then = ",age-DT,"Gyr"
print "* Comoving horizon of the universe now = ",size,"Gyr"
print "* Comoving horizon of the universe then = ",size/(1+z),"Gyr"
print "* Comoving distance of the source now = ",DC,"Gly"
print "* Comoving distance of the source then = ",DC/(1+z),"Gly"# In a flat
universe, this is equal to DA
print "* Angular Diameter distance = ",DA,"Gly"
print "* Luminosity distance = ",DL,"Gly"
print "* Light Travel Time Distance = ",DT,"Gly"
print "-------------------------------------------------------------------"
/*
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");
}
_______________________________________________
Tutor maillist - [email protected]
To unsubscribe or change subscription options:
http://mail.python.org/mailman/listinfo/tutor
