Commit c5c1c14d authored by Phil Bull's avatar Phil Bull

Make (basic) geodesic code more flexible, play with v_p plotting again

parent 31ab48e1
......@@ -23,16 +23,16 @@ int main(int argc, char* argv[]){
//FRW mod(0.72, 0.3, 0.6); // (h, omega_m, omega_k) almost-Milne
///GBH mod(0.65, 0.13, 2e3, 1e3); // (h, omega_in, r0, deltar) Garcia-Bellido and Haugbolle
//mod.output_background();
//mod.output_geodesic();
mod.output_background();
mod.output_geodesic();
//mod.output_contour();
//mod.output_velocities();
//mod.output_velocity_mismatch();
//mod.output_model_stats();
mod.output_model_stats();
//mod.output_fnz();
// TESTING
mod.v_dipole(2000.0, 40.0);
//mod.v_dipole(2000.0, 40.0);
// Finish timing the run
cerr << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
......
......@@ -27,8 +27,8 @@ const double rs = 100.0; // 100 MPC
#include "model_geodesic_interp.cpp"
#include "model_grid_build.cpp"
#include "model_grid_interp.cpp"
//#include "model_velocity_profile.cpp"
#include "model_velocity_profile_testing.cpp"
#include "model_velocity_profile.cpp"
//#include "model_velocity_profile_testing.cpp"
#include "model_obs.cpp"
#include "model_params.cpp"
#include "model_def.cpp"
......
......@@ -7,27 +7,31 @@ void Model::generate_geodesic(){
clock_t tclock = clock();
double k1, k2, k3, k4;
dt = (t0 - t_today)/1e6;
double dtr = dt * C; // dt in distance units
double dt2 = (t0 - t_today)/1e6;
double dtr = dt2 * C; // dt in distance units
double dr = 0.0; double r = r0;
double tlast = t_today; // How far we've travelled in t since the last sample
double tlast = 2.*t_today; // How far we've travelled in t since the last sample. Set to large value initially so we get the first sample.
double dI = 0.0; double dIlast = 0.0; double I = 0.0;
double zz = 0.0; // Temporary redshift variable
double tcell2 = tcell;
bool switched = false;
// Add initial values to the r,t arrays
rval.push_back(r); tval.push_back(t_today);
rval.push_back(r); tval.push_back(t_today); zval.push_back(0.0);
cerr << "Integration must stop before t0=" << t0 << endl;
// FIXME: Needs better high-z behaviour, very inaccurate ATM
// FIXME: have to stop a step early, because it isn't dealing with t+dt < t_end
for(double t=t_today; t > t0 - dt; t+=dt){
for(double t=t_today; t > t0-dt2; t+=dt2){
// Calculate Runge-Kutta variables
k1 = -1. * sqrt(1. - k(r)*r*r) / a2(r, t);
k2 = -1. * sqrt(1. - k(r + 0.5*dtr*k1)*(r + 0.5*dtr*k1)*(r + 0.5*dtr*k1))
/ a2(r + 0.5*dtr*k1, t + 0.5*dt);
/ a2(r + 0.5*dtr*k1, t + 0.5*dt2);
k3 = -1. * sqrt(1. - k(r + 0.5*dtr*k2)*(r + 0.5*dtr*k2)*(r + 0.5*dtr*k2))
/ a2(r + 0.5*dtr*k2, t + 0.5*dt);
/ a2(r + 0.5*dtr*k2, t + 0.5*dt2);
k4 = -1. * sqrt(1. - k(r + dtr*k3)*(r + dtr*k3)*(r + dtr*k3))
/ a2(r + k3*dtr, t + dt);
/ a2(r + k3*dtr, t + dt2);
// Add value to r
dr = (dtr/6.)*(k1 + 2.*k2 + 2.*k3 + k4);
......@@ -35,22 +39,35 @@ void Model::generate_geodesic(){
// Add value to redshift integral
// FIXME: Doesn't properly deal with the first couple of values!
dI = -1. * a2dot(r, t) / a2(r, t);
dI = -1. * dt2 * a2dot(r, t) / a2(r, t);
I += dI + dIlast;
dIlast = dI;
// Add the next sample, if it's time to do so
if((tlast - t) > tcell){
if((tlast - t) > tcell2){
tlast = t;
rval.push_back(r);
tval.push_back(t);
// This is OK, the dimensions work out; C is cancelled in
// the calculation, and we're only multiplying by the
// constant factor (dt/2) once for speed.
zval.push_back( exp(I * 0.5 * dt) - 1. );
zz = exp(I * 0.5) - 1.;
zval.push_back( zz );
// Check if we need to change the integration step size yet (getting to high z)
if(zz > 20. and !switched){
dt2 *= 1e-2;
dtr = dt2 * C;
tcell2 *= 1e-2;
switched = true;
cerr << "Switched! Now, t=" << t << ", dt=" << dt2 << ", tcell=" << tcell2 << endl;
}
} // end sampling check
} // end integration loop
cerr << "\tgenerate_geodesic(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
}
......@@ -79,13 +79,9 @@ double Model::fn_r(double t){
double Model::fn_t(double r){
/// Find t given r by doing a linear interpolation
// Get the index of the nearest value of r
//int i = (int)floor((r - g->r0)/tcell);
// Get the index of the nearest value of r (need to search t)
int a = 0; int b = rval.size()-1;
int mid;
//bool found = false;
// Check that the r being asked for is sensible
if(!(rval[a] <= r and rval[b] >= r)){
......
......@@ -27,13 +27,18 @@ double Model::get_array_value(double r, double t, int v){
// Sanity checks on x, y
if(x < g_x[a]){
cerr << "ERROR: x=" << x << " outside calculated range: " << g_x[a] << " to " << g_x[b] << endl; throw 71;
if(x < 0.0){
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[a] << " to " << g_x[b] << endl; throw 71;
} // end check whether x is +ve
}
// FIXME: Sometimes this gets thrown when compiled with -O3 (but seemingly doesn't happen when DEBUG flags used)
// Perhaps should tweak it just to warn, or to be *less* sensitive?
// Yes, y-y_today is *very* small, so must be a precision thing?
if( y < y0 or (y > y_today and (y - y_today > 1e-10)) ){
cerr << "ERROR: y=" << y << " outside calculated range y0=" << y0 << ", y_today=" << y_today << ", y-y_today=" << y-y_today << endl; throw 72;
cerr << "ERROR: t=" << y*ts << " outside calculated range t0=" << ts*y0 << ", t_today=" << ts*y_today << ", y-y_today=" << y-y_today << endl; throw 72;
}
// Check if x is too big, and if so, get the asymptotic value
......
......@@ -75,8 +75,8 @@ public:
double z_EdS(double r);
double dA(double r);
double v_dipole(double r, double ZLIM); // Use this definition if you're doing integrand output for off-centre obs.
//double v_dipole(double r);
//double v_dipole(double r, double ZLIM); // Use this definition if you're doing integrand output for off-centre obs.
double v_dipole(double r);
double outgoing_redshift(double r_start, double t_start, double t_end);
double ingoing_redshift(double r_start, double t_start, double t_end);
......
......@@ -91,7 +91,7 @@ void Model::output_velocities(){
clock_t tclock = clock();
stringstream fname(stringstream::in | stringstream::out);
fname << "data/velocity-" << 40;
fname << "data/velocity";
ofstream vel; vel.open(fname.str().c_str());
// Get the r when z~5
......@@ -102,7 +102,7 @@ void Model::output_velocities(){
double tt = fn_t(r);
vel << r << "\t"
<< z(r) << "\t"
<< v_dipole(r, 40.0) << "\t"
<< v_dipole(r) << "\t"
<< (Hr(r+20.0, tt) - Hr(r, tt))/20.0
<< endl;
}
......@@ -127,7 +127,7 @@ void Model::output_velocity_mismatch(){
vel << r << "\t"
<< zz << "\t"
<< C*zz - (fH0*h*dL(r)) << "\t"
<< -1.*v_dipole(r, 40.0) // FIXME: Remove this def.
<< -1.*v_dipole(r)
<< endl;
}
vel.close();
......
......@@ -56,16 +56,16 @@ for s in rad:
print "r_observer =", s, "-- Inbound (in):", len(gi1[i]), "-- Inbound (out):", len(gi2[i]), "-- Outbound:", len(go[i])
subplot(331+i)
plot(ro[i], go[i], 'kx', ro[i], go[i], 'k-')
plot(ri1[i], gi1[i], 'b+', ri1[i], gi1[i], 'b-')
plot(ri2[i], gi2[i], 'r1', ri2[i], gi2[i], 'r-')
plot(to[i], go[i], 'kx', to[i], go[i], 'k-')
plot(ti1[i], gi1[i], 'b+', ti1[i], gi1[i], 'b-')
plot(ti2[i], gi2[i], 'r1', ti2[i], gi2[i], 'r-')
title("r = " + str(rad[i]))
yscale('log')
xscale('log')
xlim((5e0, 5e3))
ylim((5e-5, 5e-4))
#xlim((5e0, 5e3))
#ylim((5e-5, 5e-4))
grid(True)
xlim((5e0, 1e4)) # For the last plot
ylim((1e-4, 1e-1))
#xlim((5e0, 1e4)) # For the last plot
#ylim((1e-4, 1e-1))
show()
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