Commit 11328279 authored by Phil Bull's avatar Phil Bull

Add dimensions to general offcentre geodesic solver, some bug-hunting, import...

Add dimensions to general offcentre geodesic solver, some bug-hunting, import data and write analysis tools for SNe
parent e14e0b83
......@@ -39,14 +39,22 @@ int main(int argc, char* argv[]){
// TESTING
cerr << "g/PI\tz" << endl;
double tend = mod.fn_tz(1.0);
double rstart = 1.0;
/*
for(double g=0.0; g<2.*M_PI; g+=2.*M_PI/40.0){
cerr << g/M_PI << "\t" << mod.general_offcentre(0.01, mod.t_today, g, 1e3) << endl;;
for(double g=0.0; g<=2.*M_PI; g+=2.*M_PI/20.0){
cerr << g/M_PI << "\t" << mod.general_offcentre(rstart, mod.t_today, g, tend) << endl;
} // end loop through gamma angles
*/
cerr << 0.0 << "\t" << mod.general_offcentre(0.01, mod.t_today, 0.0, 1e3) << endl;
cerr << M_PI << "\t" << mod.general_offcentre(0.01, mod.t_today, M_PI, 1e3) << endl;
cerr << 2.*M_PI << "\t" << mod.general_offcentre(0.01, mod.t_today, 2.*M_PI, 1e3) << endl;
cerr << 0.0 << "\t" << mod.general_offcentre(rstart, mod.t_today, 0.0, tend) << endl;
cerr << endl << endl << endl;
cerr << M_PI << "\t" << mod.general_offcentre(rstart, mod.t_today, M_PI, tend) << endl;
//cerr << M_PI << "\t" << mod.general_offcentre(rstart, mod.t_today, M_PI, tend) << endl;
//cerr << 2.*M_PI << "\t" << mod.general_offcentre(rstart, mod.t_today, 2.*M_PI, tend) << endl;
// Finish timing the run
cerr << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
......
......@@ -23,7 +23,7 @@ int Model::get_x_index(double x){
// Sanity checks on x, y
if(x < g_x[0]){
if(x < 0.0){
cerr << "WARNING: x went negative, taking fabs(x) [Model::get_array_value]" << endl;
//cerr << "WARNING: x went negative, taking fabs(x) [Model::get_array_value]" << endl;
x = fabs(x); // Try the +ve value of x
}else{
cerr << "ERROR: x=" << x << " outside calculated range: " << g_x[0] << " to " << g_x[gxsize-1] << endl; throw 71;
......
This diff is collapsed.
#!/usr/bin/python
# SNCAT Catalogue from http://www.sai.msu.su/sn/sncat/, downloaded 28-Jan-2011
Sname = [] # Supernova name, field 0
z = [] # Redshift, field 10
RA = [] # Supernova RA, field 25
Dec = [] # Supernova Dec, field 26
# Union2 catalogue from http://supernova.lbl.gov/Union/descriptions.html#FullTable
Uname = [] # Supernova name, field 0
Uz = [] # Supernova redshift, field 1
DM = [] # Supernova distance modulus, field 8
eDM = [] # Error on the supernova distance modulus, field 9
Ui = [] # The index in the SNCAT catalogue corresponding to this supernova (-1 == no match)
def convertRA(rastr):
"""Convert an RA string into a decimal"""
rastr = rastr.replace("- ", "-")
rastr = rastr.replace("+ ", "+")
r = rastr.split(" ")
return float(r[1])*(360./24.) + float(r[2])*(360./(24.*60.)) + float(r[2])*(360./(24.*3600.))
def convertDec(decstr):
"""Convert a Dec string into a decimal"""
decstr = decstr.replace("- ", "-")
decstr = decstr.replace("+ ", "+")
r = decstr.split(" ")
# Get the sign right
if float(r[1]) < 0.0:
sgn = -1.
else:
sgn = 1.
return sgn * ( sgn*float(r[1]) + (float(r[2])/60.) + (float(r[2])/3600.) )
# Read the SNe database file
f = open('sncat_latest_view.txt', 'r')
for line in f.readlines():
line = line.split(" |")
if len(line) > 26:
try:
# Get the variables from the file, if we can
ss = line[0].replace(" ", "")
zz = float(line[10].replace(" ", ""))
rr = convertRA(line[25])
dd = convertDec(line[26])
# Add the variables to the arrays
Sname.append(ss)
z.append(zz)
RA.append(rr)
Dec.append(dd)
except:
pass
f.close()
# Read the cosmology SNe z/deltaDM database
f = open('tableSCPUnion2.dat', 'r')
for line in f.readlines():
line = line.split("\t")
if len(line) > 9:
try:
# Get the variables from the file, if we can
uu = line[0].replace(" ", "")
zz = float(line[1])
dd = float(line[8])
ee = float(line[9])
# Add the variables to the arrays
Uname.append(uu)
Uz.append(zz)
DM.append(dd)
eDM.append(ee)
Ui.append(-1)
except:
pass
f.close()
# Search SNe and try to cross-match between Union2 and SNCAT
for i in range(len(Uname)):
# Loop through each item in Union2
candidate = -1
for j in range(len(Sname)):
# Loop through each item in SNCAT
# SNCAT names are more complicated than Union2 names - they have prefixes like "HST"
if Uname[i].lower() in Sname[j].lower():
candidate = j
break
if candidate >= 0:
Ui[i] = candidate
num = 0
for i in range(len(Uname)):
if(Ui[i]>=0):
num += 1
j = Ui[i]
print Uname[i], Sname[j], Uz[i], z[j]
print num, "matches were found, out of a possible", len(Uname), "(" + str(round(100.0*num/len(Uname), 2)) + "%)"
This diff is collapsed.
This diff is collapsed.
......@@ -13,15 +13,17 @@ double Model::general_offcentre(double robs, double tobs, double _gamma, double
double dr, dt, dth, dp;
double J = a1(r, t) * sin(g); // Angular momentum term
double J = C * r * a1(r, t) * sin(g); // Angular momentum term
double p = sqrt(1. - k(r)*r*r) * cos(g) / a2(r, t); // dr/du, change in radial coord along geodesic
double invq, inva1r, okr, a2d, a2a, a2p; // Use invq = 1/q, since multiplication is faster than division
double dz = 1e-2;
double dz = 1e-3;
double z0 = z(robs);
double z;
J = 0.0;
//cerr << "(r0, t0, th0, p0) = " << r << ", " << t << ", " << th << ", " << p << endl;
//cerr << "J = " << J << endl;
int i = 0;
for(z=z0; t>tend; z+=dz){
......@@ -31,36 +33,37 @@ double Model::general_offcentre(double robs, double tobs, double _gamma, double
okr = (1. - k(r)*r*r);
a2d = a2dot(r, t);
a2a = a2(r, t);
a2p = 0.05*(a2(r+20.0, t) - a2(r, t)); // a2', radial derivative of a2
a2p = 0.1*(a2(r+10.0, t) - a2(r, t)); // a2', radial derivative of a2
invq = 1. / (
( a2a * a2d * p*p / okr )
+ ( J*J*r*a1dot(r, t) * inva1r*inva1r*inva1r )
( a2a * a2d * p*p / (okr*C) )
+ ( J*J*r*a1dot(r, t) * inva1r*inva1r*inva1r/(C*C*C) )
);
// Do simple z-step forward
dt = dz * -1.0 * (1.0 + z) * invq;
dt = dz * -1.0 * (1.0 + z) * invq / C;
dr = dz * p * invq;
dth = dz * J * invq * inva1r*inva1r;
dth = dz * J * invq * inva1r*inva1r / C;
dp = dz * invq * (
( okr * J*J * inva1r*inva1r*inva1r / (r*a2(r,t)) )
+ ( 2.*p*(1.0 + z)*a2d / (a2a) )
( okr * J*J * inva1r*inva1r*inva1r / (a2a*C*C) )
+ ( 2.*p*(1.0 + z)*a2d / (a2a*C) )
- ( p*p * ( (a2p/a2a) + (0.5*(kprime(r)*r*r + 2.*k(r)*r)/okr )) )
);
// Increment values
if(isnan(dt)){cerr << "r=" << r << ", t=" << t << ", th=" << th << ", p=" << p << endl;}
//if(i==1){
//cerr << r << "\t" << t << "\t" << th << "\t" << p << "\t" << z << endl;// i=0;}
cout << t << "\t" << p << "\t" << dp << endl;
rr += dr; r = fabs(rr); // FIXME: Is this right?
t += dt;
th += dth;
p += dp;
i+=1;
if(i==1000000){cerr << "\t\t\t" << r << ", " << t << endl; i=0;}
if(isnan(dt)){cerr << "r=" << r << ", t=" << t << ", th=" << th << ", p=" << p << endl; throw 7;}
//i+=1;
//if(isnan(dt)){cerr << "r=" << r << ", t=" << t << ", th=" << th << ", p=" << p << endl; throw 7;}
} // end loop through z values
//cerr << "(r, t, th, p) = " << r << ", " << t << ", " << th << ", " << p << endl;
cerr << "\tgeneral_offcentre(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
//cerr << "\tgeneral_offcentre(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
return z;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment