Commit 91656a4f authored by Phil Bull's avatar Phil Bull

Refactored code into all-in-one class structure

parent a4fe8460
using namespace std;
#include "bubble.h"
#include "model.h"
#include "models_cfl.cpp"
int main(int argc, char* argv[]){
/// Main function
clock_t tclock = clock();
/*
// Use command line arguments to define model parameters
if(argc<2){cerr << "ERROR: Not enough command-line arguments." << endl; return 1;}
double h = atof(argv[1]);
double A = atof(argv[2]);
double KR = atof(argv[3]);
double M0 = atof(argv[4]);*/
// Instantiate a model
CFLb mod(0.72, -5.75e-8, 1.0e-9); // (h, A, KR)
// Define model
//CFLb v(0.72, -1e-11, 1.0e-10, 1e8); // (h, A, KR, M0)
//CFLb v(0.72, -1e-11, 1.0e-10); // (h, A, KR)
///CFLb v(0.72, -1e-9, 1.0e-10); // (h, A, KR)
///////CFLb v(0.72, -5e-8, 1.0e-10); // (h, A, KR)
//CFLb v(0.72, -5.75e-8, 1.0e-10); // (h, A, KR)
CFLb v(0.72, -5.75e-8, 1.0e-9); // (h, A, KR)
// Build a grid and geodesics
NumericGrid g(&v, 300);
Geodesic geo(&v, &g);
// Output datafiles
output_background(&v, &g, &geo);
//output_grid_tests(&g);
output_geodesic(&g, &geo);
//output_velocities(&v, &geo, &g); // Dipole velocities
output_model_stats(&g, &geo);
output_a1dot(&v);
//output_background(&gbh, &ggbh, &geo_gbh);
//output_geodesic(&ggbh, &geo_gbh);
//output_model_stats(&ggbh, &geo_gbh);
// FIXME: Output this as a model statistic in output_model_stats
//cerr << "Maximum geodesic r value: " << geo.rval[geo.rval.size()-1] << endl;
mod.output_background();
mod.output_geodesic();
mod.output_model_stats();
// Finish timing the run
cerr << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
......
using namespace std;
#include "bubble.h"
int main(int argc, char* argv[]){
/// Main function
clock_t tclock = clock();
/*
// Use command line arguments to define model parameters
if(argc<2){cerr << "ERROR: Not enough command-line arguments." << endl; return 1;}
double h = atof(argv[1]);
double A = atof(argv[2]);
double KR = atof(argv[3]);
double M0 = atof(argv[4]);*/
// Define model
//CFLb v(0.72, -1e-11, 1.0e-10, 1e8); // (h, A, KR, M0)
//CFLb v(0.72, -1e-11, 1.0e-10); // (h, A, KR)
///CFLb v(0.72, -1e-9, 1.0e-10); // (h, A, KR)
///////CFLb v(0.72, -5e-8, 1.0e-10); // (h, A, KR)
//CFLb v(0.72, -5.75e-8, 1.0e-10); // (h, A, KR)
CFLb v(0.72, -5.75e-8, 1.0e-9); // (h, A, KR)
// Build a grid and geodesics
NumericGrid g(&v, 300);
Geodesic geo(&v, &g);
// Output datafiles
output_background(&v, &g, &geo);
//output_grid_tests(&g);
output_geodesic(&g, &geo);
//output_velocities(&v, &geo, &g); // Dipole velocities
output_model_stats(&g, &geo);
output_a1dot(&v);
//output_background(&gbh, &ggbh, &geo_gbh);
//output_geodesic(&ggbh, &geo_gbh);
//output_model_stats(&ggbh, &geo_gbh);
// FIXME: Output this as a model statistic in output_model_stats
//cerr << "Maximum geodesic r value: " << geo.rval[geo.rval.size()-1] << endl;
// Finish timing the run
cerr << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
return 0;
}
// Header for main.cpp
#include <cmath>
#include <ctime>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <vector>
// Units: Mass (Msun); Time (Myr); Distance (Mpc)
const double G = 4.4986e-21; // Gravitational constant G, in units of Mpc^3 Msun^-1 Myr^-2
const double C = 0.30659458; // Speed of light, in Mpc/Myr
const double GYR = 1000.; // One Gyr in units of time (Myr)
const double fH0 = 1.02269e-4; // H0 = 100h km s^-1 Mpc^-1 = fH0*h Myr^-1
const double A1_REC = 1e-3; // FRW scale factor at recombination
const double T_REC = 0.4; // Age of universe at recombination (Myr)
const double A_MAX = 1.0; // The maximum a1 value to calculate
const double ts = 1000.0; // 1 GYR
const double rs = 100.0; // 100 MPC
#include "model_init.cpp"
#include "model_geodesic_build.cpp"
#include "model_geodesic_interp.cpp"
#include "model_grid_build.cpp"
#include "model_grid_interp.cpp"
#include "model_obs.cpp"
#include "model_params.cpp"
#include "model_def.cpp"
#include "model_output.cpp"
/*
const double G = 6.67e-11; // Gravitational constant G, in units of m^3 kg^-1 s^-2
const double C = 3e8; // Speed of light, in m/s
const double GYR = 3.16e16; // One Gyr in units of time (s)
const double fH0 = 3.24e-18; // H0 = 100h km s^-1 Mpc^-1 = fH0*h Myr^-1
const double A1_REC = 1e-3; // FRW scale factor at recombination
const double T_REC = 1.26e13; // Age of universe at recombination (400 kyr)
const double A_MAX = 1.0; // The maximum a1 value to calculate
const double ts = 1e13; // ~ 1 Myr
const double rs = 3e24; // ~ 100 MPC
*/
/*
class Model;
void Model::generate_geodesic();
double Model::z(double r);
double Model::fn_r(double t);
double Model::fn_t(double r);
double Model::fn_tz(double z);
double Model::xa1dot(double a1, double x);
double Model::xa1dot(double a1, double m, double k);
void Model::populate_grid();
double Model::a1(double r, double t);
double Model::a1dot(double r, double t);
double Model::a2(double r, double t);
double Model::a2dot(double r, double t);
double Model::get_array_value(double r, double t, int v);
void Model::build_grid();
void Model::hubble_time();
*/
double Model::m(double r){
/// Density-like function m(r)
cerr << "Model::m" << endl;
return 0.0;
}
double Model::mprime(double r){
/// Radial derivative of the density-like function,. m'(r)
return 0.0;
}
double Model::k(double r){
/// Curvature-like function k(r)
return 0.0;
}
double Model::kprime(double r){
/// Radial derivative of the curvature-like function, k'(r)
return 0.0;
}
double Model::tb(double r){
/// Bang-time function t_b(r)
return 0.0;
}
double Model::tbprime(double r){
/// Radial derivative of the bang-time function, t_b'(r)
return 0.0;
}
// Dimensionless functions to be used in the grid integration
double Model::xM(double x){
/// Dimensionless density-like function, M(x) = 8*pi*G/3 * ts^2 * m(x)
return ((8.*M_PI*G)/3.) * ts*ts * m(x*rs);
}
double Model::xMprime(double x){
/// Dimensionless radial derivative of density-like function, M'(x)
/// M'(x) = 8*pi*G/3 * ts^2 * rs * dm(r)/dr
return ((8.*M_PI*G)/3.) * ts*ts * rs * mprime(x*rs);
}
double Model::xK(double x){
/// Dimensionless curvature-like function, K(x) = k(x) * ts^2 * C^2
return k(rs*x) * ts*ts * C*C;
}
double Model::xKprime(double x){
/// Dimensionless radial derivative of curvature-like function, K'(x)
/// K'(x) = k(x) * ts^2 * C^2 * rs
return kprime(rs*x) * ts*ts * C*C * rs;
}
double Model::xTb(double x){
/// Bang-time function, t_B(x)
// FIXME
return tb(x*rs);
}
double Model::xTbprime(double x){
/// Radial derivative of Bang-time function, t_B'(x)
// FIXME
return tbprime(x*rs);
}
void Model::generate_geodesic(){
/// Build a geodesic by integrating along dt to find (r,t) at
/// multiple sample points
/// Also create a geodesic sample
clock_t tclock = clock();
double k1, k2, k3, k4;
dt = (t0 - t_today)/1e6;
double dtr = dt * 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 dI = 0.0; double dIlast = 0.0; double I = 0.0;
//double dI2 = 0.0; double dIlast2 = 0.0; double I2 = 0.0;
cerr << "Start generate_geodesic()" << endl;
// Add initial values to the r,t arrays
rval.push_back(r); tval.push_back(t_today);
// FIXME: have to stop a step early, because it isn't dealing with t+dt < t_end
cerr << "*** end time: " << t0 - dt << endl;
for(double t=t_today; t > t0 - dt; t+=dt){
// 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);
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);
k4 = -1. * sqrt(1. - k(r + dtr*k3)*(r + dtr*k3)*(r + dtr*k3))
/ a2(r + k3*dtr, t + dt);
// Add value to r
dr = (dtr/6.)*(k1 + 2.*k2 + 2.*k3 + k4);
r += dr;
// Add value to redshift integral
// FIXME: Doesn't properly deal with the first couple of values!
dI = -1. * 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){
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. );
} // end sampling check
} // end integration loop
cerr << "Geodesic::generate_geodesic(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
}
double Model::z(double r){
/// Get the redshift z(r)
// Get the index of the nearest value of r (need to search t)
int a = 0; int b = rval.size()-1;
int mid;
// Check that the r being asked for is sensible
if(!(rval[a] <= r and rval[b] >= r)){
cerr << "WARNING: r=" << r << " falls outside the range of the calculated values " << rval[a] << " to " << rval[b] << ". geodesic.cpp::z()" << endl;
return -0.0;
}
// Search for the cells bounding the correct time
while( (b - a) > 1 ){
mid = a + floor(0.5 * (double)(b - a)); // Find the middle cell
if(r <= rval[mid]){
b = mid; // a <= r <= mid
}else{
a = mid; // mid < r <= b
} // end range check
} // end searching loop
// Check if we're hitting a boundary
if(a == zval.size() - 1){
// FIXME: Use the *backwards* gradient to do the interpolation
// For now, just return the latest value
return zval[a];
}else{
// Do proper forwards linear interpolation
return zval[a] + (r - rval[a]) * (zval[a+1] - zval[a]) / (rval[a+1] - rval[a]);
} // end boundary check
}
double Model::fn_r(double t){
/// Find r given t by doing a linear interpolation
int a = 0; int b = tval.size()-1;
int mid;
// Check that the t being asked for is sensible
if(!(tval[a] >= t and tval[b] <= t)){
cerr << "WARNING: t=" << t << " falls outside the range of the calculated values " << tval[a] << " to " << tval[b] << ". geodesic.cpp::fn_r()" << endl;
return -0.0;
}
// Search for the cells bounding the correct r
while( (b - a) > 1 ){
mid = a + floor(0.5 * (double)(b - a)); // Find the middle cell
if(t >= tval[mid]){
b = mid; // a <= t <= mid
}else{
a = mid; // mid < t <= b
} // end range check
} // end searching loop
// Get the index of the nearest value of t
//int a = (int)floor((g->t_today - t)/tcell);
// Check if we're hitting a boundary
if(a == tval.size() - 1){
// FIXME: Use the *backwards* gradient to do the interpolation
// For now, just return the latest value
return rval[a];
}else{
// Do proper forwards linear interpolation
// FIXME: Check that this is interpolating the right way!
return rval[a] + (tval[a] - t) * (rval[a+1] - rval[a]) / (tval[a+1] - tval[a]);
} // end boundary check
}
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)){
cerr << "WARNING: r=" << r << " falls outside the range of the calculated values " << rval[a] << " to " << rval[b] << ". geodesic.cpp::fn_t()" << endl;
return -0.0;
}
// Search for the cells bounding the correct time
while( (b - a) > 1 ){
mid = a + floor(0.5 * (double)(b - a)); // Find the middle cell
if(r <= rval[mid]){
b = mid; // a <= r <= mid
}else{
a = mid; // mid < r <= b
} // end range check
} // end searching loop
// Check if we're hitting a boundary
if(a == rval.size() - 1){
// FIXME: Use the *backwards* gradient to do the interpolation
// For now, just return the latest value
return tval[a];
}else{
// Do proper forwards linear interpolation
return tval[a] + (r - rval[a]) * (tval[a+1] - tval[a]) / (rval[a+1] - rval[a]);
} // end boundary check
}
double Model::fn_tz(double z){
/// Find t given z by doing a linear interpolation
int a = 0; int b = zval.size()-1;
int mid;
// Check that the z being asked for is sensible
if(!(zval[a] <= z and zval[b] >= z)){
cerr << "WARNING: z=" << z << " falls outside the range of the calculated values " << zval[a] << " to " << zval[b] << ". geodesic.cpp::fn_tz()" << endl;
return -0.0;
}
// Search for the cells bounding the correct r
while( (b - a) > 1 ){
mid = a + floor(0.5 * (double)(b - a)); // Find the middle cell
if(z <= zval[mid]){
b = mid; // a <= z <= mid
}else{
a = mid; // mid < z <= b
} // end range check
} // end searching loop
// Check if we're hitting a boundary
if(a == zval.size() - 1){
// FIXME: Use the *backwards* gradient to do the interpolation
// For now, just return the latest value
return tval[a];
}else{
// Do proper forwards linear interpolation
// FIXME: Check that this is interpolating the right way!
return tval[a] + (zval[a] - z) * (tval[a+1] - tval[a]) / (zval[a+1] - zval[a]);
} // end boundary check
}
double Model::xa1dot(double a1, double x){
/// Calculate dimensionless a1dot, given a1 and r
return sqrt( (xM(x) / a1) - xK(x) );
}
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)
//cerr << (m / a1) / (-k) << endl;
return sqrt( (m / a1) - k );
}
void Model::populate_grid(){
/// Populate a mesh with values for a1 using an RK4 ODE solution
/// to the LTB Friedmann equation
clock_t tclock = clock();
double a1, x, y, ylast, m, k;
double k1, k2, k3, k4; // Runge-Kutta variables
bool reached_asymptotic = false;
int j;
x = x0 - XSTEP; // Make sure we start at the right place
// Loop through x values
for(int i=0; i<XCELLS+2 and !reached_asymptotic; i++){
// Set-up for integration
x += XSTEP;
g_x.push_back(x); // Add the x value for this cell to the array
j = 0; a1 = a1_0; ylast = 0.0;
m = xM(x); k = xK(x);
// Create new temporary time slice vectors
// Add an extra 2 elements to each slice to avoid boundary problems when interpolating
vector<double> tmp_a1(YCELLS+2, 0.0);
vector<double> tmp_a1dot(YCELLS+2, 0.0);
vector<double> tmp_a2(YCELLS+2, 0.0);
vector<double> tmp_a2dot(YCELLS+2, 0.0);
for(y=y0; j<YCELLS+2; y+=dy){
k1 = xa1dot(a1, m, k);
k2 = xa1dot(a1 + 0.5*k1*dy, m, k);
k3 = xa1dot(a1 + 0.5*k2*dy, m, k);
k4 = xa1dot(a1 + k3*dy, m, k);
//if(j==298){cerr << "\tk1 = " << k1 << endl;}
// Calculate new a1
a1 += (dy/6.) * (k1 + 2.*k2 + 2.*k3 + k4);
//cerr << "da1 = " << (dy/6.) * (k1 + 2.*k2 + 2.*k3 + k4) << endl;
// Add to mesh if we're at a sample point
if(y - ylast >= YSTEP){
tmp_a1[j] = a1;
//if(j==298){ cerr << a1 << endl; }
tmp_a1dot[j] = k1 / ts; // k1=a1dot
if(i>1){
// FIXME: This isn't accurate enough. Should probably do a different integration type!
tmp_a2[j] = g_a1[i-1][j] + x*XSTEP*(a1 - g_a1[i-2][j]) ;
tmp_a2dot[j] = ( g_a1dot[i-1][j] + x*XSTEP*(k1/ts - g_a1dot[i-2][j]) );
// Check for shell crossings
if(tmp_a2[j] <= 0.0){cerr << "WARNING: Shell crossing occurred at x=" << x << ", y=" << y << " (NumericGrid::populate_grid)" << endl; throw 99;} // end check
}
ylast = y; j++;
} // end sample-point check
} // end time integration
g_a1.push_back( tmp_a1 );
g_a1dot.push_back( tmp_a1dot );
g_a2.push_back( tmp_a2 );
g_a2dot.push_back( tmp_a2dot );
if(x>XMIN){
// Check if we've reached an asymptotic value (at y=today)
// Beware the minimum in a2! Need to set xmin to avoid this
if( (g_a2[i-1][YCELLS-1] - g_a2[i-2][YCELLS-1])
/ (g_a2[i-1][YCELLS-1] - g_a2[i-6][YCELLS-1]) < 1e-14 ){
cerr << "\tAsymptotic value reached at x=" << x << ", a2=" << g_a2[i-1][YCELLS-1] << endl;
reached_asymptotic = true;
XCELLS = g_a1.size() - 2;
XMAX = g_x[XCELLS-1];
}
} // end asymptote check
} // end loop through x values
// Tidy-up values that couldn't be calculated efficiently in the loop
for(int jj=0; jj<YCELLS+2; jj++){
// be careful with the x being used - x should be xlast, but use x0 for the first ones
g_a2[0][jj] = g_a1[0][jj] + x0*2.*(g_x[1] - g_x[0])* (g_a1[1][jj] - g_a1[0][jj]);
g_a2[1][jj] = g_a1[1][jj] + x0*2.*(g_x[2] - g_x[1])* (g_a1[2][jj] - g_a1[1][jj]);
g_a2dot[0][jj] = g_a1dot[0][jj] + x0*2.*(g_x[1] - g_x[0])* (g_a1dot[1][jj] - g_a1dot[0][jj]);
g_a2dot[1][jj] = g_a1dot[1][jj] + x0*2.*(g_x[2] - g_x[1])* (g_a1dot[2][jj] - g_a1dot[1][jj]);
}
cerr << "NumericGrid::populate_grid(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
} // end initial integration function
double Model::a1(double r, double t){
return get_array_value(r, t, 0);
}
double Model::a1dot(double r, double t){
/// Units: inverse time
return get_array_value(r, t, 1);
}
double Model::a2(double r, double t){
return get_array_value(r, t, 2);
}
double Model::a2dot(double r, double t){
/// Units: inverse time
return get_array_value(r, t, 3);
}
double Model::get_array_value(double r, double t, int v){
/// Generic interpolation
int i, j, a, b, mid;
double x1, x2, y1, y2;
double xcorr = 0.0; double ycorr = 0.0;
a = 0; b = g_x.size() - 3; // Since XCELLS = xsize - 2
double x = r/rs;
double y = t/ts;
// 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(y < y0 or y > y_today){ cerr << "ERROR: y=" << y << " outside calculated range y0=" << y0 << ", y_today=" << y_today << endl; throw 72;}
// Check if x is too big, and if so, get the asymptotic value
if(x >= g_x[b-1]){
x = g_x[b-1]; // Set x to value where expn. vars. reach asymptotic values, it'll return the same results
}
// Get i<=>x by searching g_x
while( (b - a) > 1 ){
mid = a + floor(0.5 * (double)(b - a)); // Find the middle cell
if(x <= g_x[mid]){
b = mid; // a <= x <= mid
}else{
a = mid; // mid < x <= b
} // end range check
} // end x-searching loop
i = a;
// Get j<=>y by looking it up from the index (assumes fixed grid cell size in y)
j = (int)floor( ((y - y0) / (y_today - y0)) * (double)YCELLS);
if(j >= YCELLS){j = YCELLS - 1;} // Check we're not right at the top of the range
// Check if x, y indices are safe to do forward interpolation
double gxsize = g_x.size();
if(i >= gxsize and j < YCELLS){
// We're at the x boundary. FIXME: Check for validity
xcorr = g_x[i] - g_x[i-1]; // Correct for shift in cell used to calculate gradient
i -= 1; // We're at the x boundary, interpolate using a backwards step (fix the indices)
}else if(i < gxsize and j >= YCELLS){
ycorr = YSTEP; j -= 1; // We're at the y boundary
}else if(i >= gxsize and j >= YCELLS){
xcorr = g_x[i] - g_x[i-1];
ycorr = YSTEP; i -= 1; j-=1; // We're at the x and y boundary!
}
// Get values for x1, x2, y1, y2
x1 = g_x[i]; x2 = g_x[i+1];
y1 = y0 + ((double)j * YSTEP); y2 = y1 + YSTEP;
// Do a simple (bilinear) interpolation
double arr1, arr2, arr3, arr4;
switch(v){
case 0: arr1 = g_a1[i][j]; arr2 = g_a1[i][j+1];
arr3 = g_a1[i+1][j]; arr4 = g_a1[i+1][j+1]; break;
case 1: arr1 = g_a1dot[i][j]; arr2 = g_a1dot[i][j+1];
arr3 = g_a1dot[i+1][j]; arr4 = g_a1dot[i+1][j+1]; break;
case 2: arr1 = g_a2[i][j]; arr2 = g_a2[i][j+1];
arr3 = g_a2[i+1][j]; arr4 = g_a2[i+1][j+1]; break;
case 3: arr1 = g_a2dot[i][j]; arr2 = g_a2dot[i][j+1];
arr3 = g_a2dot[i+1][j]; arr4 = g_a2dot[i+1][j+1]; break;
default: cerr << "WARNING: No case was chosen in NumericGrid::get_array_value()" << endl;
break;
}
return (1./(YSTEP*(x2 - x1))) * (
arr1 * (ycorr + y2 - y)*(xcorr + x2 - x)
+ arr2 * (ycorr + y - y1)*(xcorr + x2 - x)
+ arr3 * (ycorr + y2 - y)*(xcorr + x - x1)
+ arr4 * (ycorr + y - y1)*(xcorr + x - x1)
);
}
class Model{
// Calculated grid variables
protected:
vector<vector <double> > g_a1;
vector<vector <double> > g_a1dot;
vector<vector <double> > g_a2;
vector<vector <double> > g_a2dot;
vector<double> g_x;
vector<double> rval;
vector<double> tval;
vector<double> zval;
double tcell; // The distance between samples on the y grid
void prepare_grid();
void populate_grid();
void generate_geodesic();
void hubble_time();
double z(double r);
double fn_r(double t);
double fn_t(double r);
double fn_tz(double z);
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);
double a2dot(double r, double t);
double get_array_value(double r, double t, int v);
double xM(double x);
double xMprime(double x);
double xK(double x);
double xKprime(double x);
double xTb(double x);
double xTbprime(double x);
public:
int XCELLS, YCELLS;
double XMAX, XMIN, XSTEP, YSTEP;
double a1_0, x0, y0, dy, y_today, t_today, r0, t0;
double dt;
double h, H0;
virtual double m(double r);
virtual double mprime(double r);
virtual double k(double r);
virtual double kprime(double r);
virtual double tb(double r);
virtual double tbprime(double r);
double density(double r, double t);
double Ht(double r, double t);
double Hr(double r, double t);
double DM(double r);
double DM_milne(double r);
double dL_Milne(double r);
double dL(double r);
double DM_EdS(double r);
double dL_EdS(double r);
double EdS_fn_r(double t);
double z_EdS(double r);
double dA(double r);
void output_background();
void output_geodesic();
void output_model_stats();
//void output_velocities();