Commit 46fabc99 authored by Phil Bull's avatar Phil Bull

Replace obsolete a1dot functions, clean-up code

parent 18727e72
int Model::static_a1dot (double t, const double a[], double adot[], void *params) {
ModelPointers p = *(ModelPointers *)params;
adot[0] = p.mod->s_a1dot(a[0], p.m, p.k);
return GSL_SUCCESS;
}
int Model::a1dot_jac(double t, const double a[], double *dfda, double dfdt[], void *params){
ModelPointers p = *(ModelPointers *)params;
gsl_matrix_view dfda_mat = gsl_matrix_view_array (dfda, 1, 1);
gsl_matrix * m = &dfda_mat.matrix;
gsl_matrix_set (m, 0, 0, p.mod->da1dot_da(a[0], p.m, p.k) );
dfdt[0] = 0.0;
dfdt[1] = 0.0;
return GSL_SUCCESS;
}
double Model::s_a1dot(double a1, double m, double k){
return sqrt( (m / a1) - k );
}
double Model::da1dot_da(double a, double m, double k){
-0.5 * m / ( a*a*sqrt( m/a - k ) );
}
//////////////////////
......
......@@ -12,7 +12,6 @@ void Model::generate_geodesic(){
double dr = 0.0; double r = r0;
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;
......@@ -50,17 +49,7 @@ void Model::generate_geodesic(){
// 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.
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-1;
dtr = dt2 * C;
tcell2 *= 1e-1;
switched = true;
}*/
zval.push_back( exp(I * 0.5) - 1. );
} // end sampling check
......@@ -68,7 +57,5 @@ void Model::generate_geodesic(){
// TESTING
cerr << "r=" << rval[rval.size()-1] << ", t=" << tval[rval.size()-1] << ", z=" << zval[rval.size()-1] << endl;
cerr << "\tgenerate_geodesic(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
}
int Model::static_a1dot (double t, const double a[], double adot[], void *params) {
/// Static wrapper function used for GSL integration of the LTB
/// Friedmann-like equation
ModelPointers p = *(ModelPointers *)params;
adot[0] = p.mod->s_a1dot(a[0], p.m, p.k);
return GSL_SUCCESS;
}
double Model::xa1dot(double a1, double x){
/// Calculate dimensionless a1dot, given a1 and r
return sqrt( (xM(x) / a1) - xK(x) );
int Model::a1dot_jac(double t, const double a[], double *dfda, double dfdt[], void *params){
/// Static wrapper function used to define the LTB Friedmann-like
/// equation's Jacobian for GSL integration
ModelPointers p = *(ModelPointers *)params;
gsl_matrix_view dfda_mat = gsl_matrix_view_array (dfda, 1, 1);
gsl_matrix * m = &dfda_mat.matrix;
gsl_matrix_set (m, 0, 0, p.mod->da1dot_da(a[0], p.m, p.k) );
dfdt[0] = 0.0;
dfdt[1] = 0.0;
return GSL_SUCCESS;
}
double Model::xa1dot(double a1, double m, double k){
/// Calculate dimensionless a1dot, given a1 and x
/// Alternative function definition which prevents multiple
/// evaluations of M(x) and K(x)
double Model::s_a1dot(double a1, double m, double k){
return sqrt( (m / a1) - k );
}
double Model::da1dot_da(double a, double m, double k){
-0.5 * m / ( a*a*sqrt( m/a - k ) );
}
......@@ -63,7 +77,7 @@ void Model::populate_grid(){
// Integrate until the next sample point
while (y < ysamp){
status = gsl_odeiv_evolve_apply (e, c, s, &sys, &y, ysamp, &hh, a1);
if (status != GSL_SUCCESS){ cerr << "Yikes!" << endl; break; }
if (status != GSL_SUCCESS){ cerr << "Model::populate_grid(): Integration failed." << endl; break; }
} // end while loop
if(isnan(a1[0])){throw 111;}
......@@ -116,10 +130,5 @@ void Model::populate_grid(){
gsl_odeiv_control_free (c);
gsl_odeiv_step_free (s);
cerr << "Density(r=0, today) = " << density(0.0, t_today) << endl;
cerr << "Density(r=500, today) = " << density(500.0, t_today) << endl;
cerr << "Density(r=0, 0.9*today) = " << density(0.0, t_today*0.9) << endl;
cerr << "Density(r=500, 0.9*today) = " << density(500.0, t_today*0.9) << endl;
cerr << "\tpopulate_grid(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
} // end initial integration function
......@@ -23,8 +23,6 @@ protected:
double z(double r);
double xa1dot(double a1, double x);
double xa1dot(double a1, double m, double k);
double a1(double r, double t);
double a1dot(double r, double t);
double a2(double r, double t);
......
......@@ -9,7 +9,6 @@ void Model::prepare_grid(){
x0 = 0.0; // Radial coordinate to start at (x=0 is the void centre)
a1_0 = A1_REC; // Initial a1 at x=0
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
......@@ -36,10 +35,10 @@ void Model::hubble_time(){
for(y=0.0; a1>A1_REC; y+=dy2){
k1 = xa1dot(a1, m, k);
k2 = xa1dot(a1 + 0.5*k1*dy2, m, k);
k3 = xa1dot(a1 + 0.5*k2*dy2, m, k);
k4 = xa1dot(a1 + k3*dy2, m, k);
k1 = s_a1dot(a1, m, k);
k2 = s_a1dot(a1 + 0.5*k1*dy2, m, k);
k3 = s_a1dot(a1 + 0.5*k2*dy2, m, k);
k4 = s_a1dot(a1 + k3*dy2, m, k);
// Calculate new a1
a1 += (dy2/6.) * (k1 + 2.*k2 + 2.*k3 + k4);
......
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