Commit 0caa24e3 authored by Phil Bull's avatar Phil Bull

Try to use GSL to do the v_p velocity calculation. Need r(t) instead of t(r)...

Try to use GSL to do the v_p velocity calculation. Need r(t) instead of t(r) in this case. Not working at the moment.
parent fba7b7ba
......@@ -28,15 +28,17 @@ int main(int argc, char* argv[]){
//GBH mod(0.66, 0.13, 2.5e3, 0.64*2.5e3); // (h, omega_in, r0, deltar) Garcia-Bellido and Haugbolle
//GBH mod(0.66, 0.3, 1.5e3, 0.3*1.5e3); // (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_velocities();
//mod.output_velocity_mismatch();
//mod.output_kszvpec_for_z();
//mod.output_model_stats();
//mod.output_fnz();
//mod.output_zin();
cerr << "ADAPTIVE v_p = " << mod.adaptive_v_dipole(100.0) << endl;
cerr << "NORMAL v_p = " << mod.v_dipole(100.0) << endl;
// TESTING
/*
......@@ -51,10 +53,11 @@ int main(int argc, char* argv[]){
cerr << r << "\t" << mod.adaptive_z(0.0, r) << "\t" << mod.z(r) << "\t" << mod.fn_r(mod.fn_tz(mod.z(r))) << endl;
}*/
/*
cerr << "r\tz" << endl;
for(double r=0.0; r<200.0; r+=20.0){
cerr << r << "\t" << mod.z(r) << endl;
}
}*/
/*
......
......@@ -17,7 +17,7 @@ void Model::generate_geodesic(){
double zz = 0.0;
rval.push_back(r); tval.push_back(t); zval.push_back(0.0);
init_integrators();
init_integrators(); // Initialise the GSL integration routines
// Loop through t cells
while (t > t0 + 0.01 + tb(r) and zz < 1200.0){
......@@ -45,11 +45,9 @@ void Model::generate_geodesic(){
zval.push_back(zz);
}
free_integrators();
free_integrators(); // Free the memory used by the GSL integration/ODE solving routines
// TESTING
cerr << "r=" << rval[rval.size()-1] << ", t=" << tval[rval.size()-1] << ", z=" << zval[rval.size()-1] << endl;
cerr << "\tArray size was " << rval.size() << endl;
cerr << "\tgenerate_geodesic(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
}
......
double Model::a2prime(double r, double t){
/// Get the approximate radial derivative of a2(r,t)
return a2(r+1.0, t) - a2(r, t); // Since dr = 1.0
}
double Model::a1(double r, double t, bool tt){
return get_array_value(r, t, 0, tt);
}
......
......@@ -38,6 +38,7 @@ protected:
void populate_timesample_vector();
void calculate_a2grid();
double a2prime(double r, double t);
double a1(double r, double t, bool tt = TIME_GLOBAL);
double a1dot(double r, double t, bool tt = TIME_GLOBAL);
double a2(double r, double t, bool tt = TIME_GLOBAL);
......@@ -74,15 +75,22 @@ public:
double rstart; // Starting point of a calculated geodesic
double rend; // End point of a calculated geodesic
double t_rstart; // Time at which the geodesic was at rstart
bool geodesic_direction; // The direction in which the geodesic is going
Model *mod; // A pointer to a model, allowing access to all of its public functions from a static function
};
double adaptive_t_fnr(double rstart, double rend, double tstart, bool ingoing = false);
static int ft_jac(double r, const double t[], double *dfdt, double dfdr[], void *params);
static int static_ft (double r, const double t[], double dfdr[], void *params);
static int static_ft (double r, const double t[], double dtdr[], void *params);
static double static_z_integrand (double r, void * params);
double adaptive_I(double rstart, double rend, double t_rstart);
double adaptive_I(double rstart, double rend, double t_rstart, bool geodesic_direction=GEODESIC_OUTGOING);
double adaptive_v_dipole(double robs);
double adaptive_r_fnt(double tstart, double tend, double rstart, bool ingoing);
static int fr_jac(double t, const double r[], double *dfdr, double dfdt[], void *params);
static int static_fr (double t, const double r[], double drdt[], void *params);
// END TESTING
double fn_r(double t);
......
......@@ -17,13 +17,24 @@ void Model::free_integrators(){
}
int Model::static_ft (double r, const double t[], double dfdr[], void *params){
int Model::static_ft (double r, const double t[], double dtdr[], void *params){
/// Static wrapper function used for GSL integration of the radial
/// geodesic equation dt/dr = f(r, t)
ModelPointers p = *(ModelPointers *)params; // Get model parameters
// Set the system of equations to be integrated (there's only one)
dfdr[0] = p.sign * p.mod->a2(fabs(r), t[0]) / (C * sqrt(1. - p.mod->k(fabs(r))*r*r));
dtdr[0] = p.sign * p.mod->a2(fabs(r), t[0]) / (C * sqrt(1. - p.mod->k(fabs(r))*r*r));
return GSL_SUCCESS;
}
int Model::static_fr (double t, const double r[], double drdt[], void *params){
/// Static wrapper function used for GSL integration of the radial
/// geodesic equation dr/dt = f(r, t)
ModelPointers p = *(ModelPointers *)params; // Get model parameters
// Set the system of equations to be integrated (there's only one)
drdt[0] = p.sign * C * sqrt(1. - p.mod->k(fabs(r[0]))*r[0]*r[0]) / p.mod->a2(fabs(r[0]), t);
return GSL_SUCCESS;
}
......@@ -45,26 +56,49 @@ int Model::ft_jac(double r, const double t[], double *dfdt, double dfdr[], void
return GSL_SUCCESS;
}
int Model::fr_jac(double t, const double r[], double *dfdr, double dfdt[], void *params){
/// Static wrapper function used to define the Jacobian of the
/// radial geodesic equation for use in GSL ODE solver
ModelPointers p = *(ModelPointers *)params; // Get model parameters
// Construct a matrix for the Jacobian
gsl_matrix_view dfdr_mat = gsl_matrix_view_array (dfdr, 1, 1);
gsl_matrix * m = &dfdr_mat.matrix;
// Set the Jacobian matrix elements
double rr = fabs(r[0]);
double okr = sqrt(1. - p.mod->k(rr)*rr*rr);
double inva2 = 1. / p.mod->a2(rr, t);
double element_dfdr =
- C * okr * inva2 * inva2 * p.mod->a2prime(rr, t)
- ( 0.5 * C * inva2 * (p.mod->kprime(rr)*rr*rr + 2.*p.mod->k(rr)*rr) / okr );
gsl_matrix_set (m, 0, 0, p.sign * element_dfdr);
dfdt[0] = 0.0; // Set df/dr = 0 initially (FIXME: is this right?)
return GSL_SUCCESS;
}
double Model::static_z_integrand (double r, void * params){
/// Static member function for use in the GSL integration: this is
/// the integrand used when calculating z = exp(I) - 1
ModelPointers p = *(ModelPointers *)params; // Get model parameters
// Find what t should be at this point
double t = p.mod->adaptive_t_fnr(p.rstart, r, p.t_rstart, GEODESIC_OUTGOING);
double t = p.mod->adaptive_t_fnr(p.rstart, r, p.t_rstart, p.geodesic_direction);
// Return the value of the integrand at this point
return p.mod->a2dot(r, t) / (C*sqrt(1. - p.mod->k(r)*r*r));
}
double Model::adaptive_I(double rstart, double rend, double t_rstart){
/// Adaptively integrate the integrand used to find z
double Model::adaptive_I(double rstart, double rend, double t_rstart, bool geodesic_direction){
/// Adaptively integrate the integrand used to find z (as fn r)
// Prepare variables which will be used by the integrator
double result, error;
ModelPointers p; p.mod = this;
ModelPointers p; p.mod = this; p.geodesic_direction = geodesic_direction;
p.rstart = rstart; p.rend = rend; p.t_rstart = t_rstart ;
// Prepare the integrator
......@@ -100,9 +134,39 @@ double Model::adaptive_t_fnr(double rstart, double rend, double tstart, bool ing
// Integrate until we reach the desired r
while (fabs(r) < fabs(rend)){
status = gsl_odeiv_evolve_apply (ode_e, ode_c, ode_s, &sys, &r, rend, &hh, t);
if (status != GSL_SUCCESS){ cerr << "Model::adaptive_geodesic(): Integration failed." << endl; break; }
if (status != GSL_SUCCESS){ cerr << "Model::adaptive_t_fnr(): Integration failed." << endl; break; }
} // end while loop
if(isnan(t[0])){cerr << "ERROR: Model::adaptive_geodesic, t is NaN\n"; throw 111;}
if(isnan(t[0])){cerr << "ERROR: Model::adaptive_t_fnr, t is NaN\n"; throw 111;}
return *t;
}
double Model::adaptive_r_fnt(double tstart, double tend, double rstart, bool ingoing){
/// Adaptively integrate the radial geodesic ODE from
/// (rstart, tstart) until tend, and return r(tend). Set ingoing to
/// change the sign of the geodesic (going into or out of void)
int status;
// Make sure we start at the right place
ModelPointers p; p.mod = this;
double hh;
double r[1] = { rstart }; // [1] just means it's a 1D array
double t = tstart;
// Set signs depending on whether the geodesic is ingoing or outgoing
if(ingoing){ p.sign = +1.0; hh = -1.0 * 1e-16;}
else{ p.sign = -1.0; hh = +1.0 * 1e-16; }
// Define system of ODEs (ODE + Jacobian)
gsl_odeiv_system sys = {static_fr, fr_jac, 1, &p};
// Integrate until we reach the desired r
while (t > tend){
status = gsl_odeiv_evolve_apply (ode_e, ode_c, ode_s, &sys, &t, tend, &hh, r);
if (status != GSL_SUCCESS){ cerr << "Model::adaptive_r_fnt(): Integration failed." << endl; break; }
} // end while loop
if(isnan(r[0])){cerr << "ERROR: Model::adaptive_r_fnt, r is NaN\n"; throw 111;}
return *r;
}
......@@ -72,6 +72,9 @@ void Model::output_model_stats(){
<< "XCELLS = " << XCELLS << endl
<< "YCELLS = " << YCELLS << endl
<< "x_max = " << XMAX << endl
<< "Central geodesic endpoint: r=" << rval[rval.size()-1]
<< ", t=" << tval[rval.size()-1]
<< ", z=" << zval[rval.size()-1] << endl
<< "------------------------------" << endl;
}
......
......@@ -11,7 +11,7 @@ void Model::prepare_grid(){
r0 = 0.0;
t0 = y0 * ts;
XMIN = 20; // The minimum x that should be reached before we start considering whether it has reached asymp. value
//XMIN = 20; // The minimum x that should be reached before we start considering whether it has reached asymp. value
XSTEP = 0.01; // FIXME: NORMALLY 0.01 // The step to take in x when populating the grid
XMAX = 50;//50; // FIXME: XMAX = 200 normally // And then had it set to 400
XCELLS = XMAX / XSTEP; // Maximum number of x cells to calculate
......
......@@ -97,6 +97,35 @@ double Model::ingoing_redshift(double r_start, double t_start, double t_end, dou
return exp(I * 0.5 * dt) - 1.; // Redshift, z. Factor of 0.5*dt comes from simple integration
}
double Model::adaptive_v_dipole(double robs){
/// Use the adaptive algorithms to calculate the offcentre
/// velocity dipole (in km/sec)
double t_lss = fn_tz(1100.0); // Find the time at which the central observer thinks the CMB formed at
double tobs = fn_t(robs);
// The r value at which the outgoing geodesic hits the LSS
// (This geodesic is along the same path as the central geodesic)
double rlss_outgoing = fn_r(t_lss);
// The r value at which the ingoing geodesic hits the LSS
double rlss_ingoing = adaptive_r_fnt(tobs, t_lss, robs, GEODESIC_INGOING);
init_integrators(); // Initialise the GSL integration routines
// Use the adaptive integrator to find the integrand of z
double I_outgoing = adaptive_I(robs, rlss_outgoing, tobs, GEODESIC_OUTGOING);
double I_ingoing = adaptive_I(robs, rlss_ingoing, tobs, GEODESIC_INGOING);
free_integrators(); // Free the memory used by the GSL integration/ODE solving routines
// Calculate z_in, z_out
double z_in = exp(I_ingoing) - 1.0;
double z_out = exp(I_outgoing) - 1.0;
return 3e5 * (z_in - z_out) / (2. + z_in + z_out);
}
double Model::v_dipole(double r, double* r_end = NULL){
/// Find the dipole velocity (derived from the redshift) at radius r
/// (in km/sec)
......
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