Commit 84d3d8cd authored by Phil Bull's avatar Phil Bull

WORKING - Tidy-up code, fix problems due to off-by-one errors, question about derivatives though

parent 973588a0
......@@ -29,7 +29,7 @@ int main(int argc, char* argv[]){
mod.output_background();
mod.output_geodesic();
//mod.output_contour();
//mod.output_velocities();
mod.output_velocities();
//mod.output_velocity_mismatch();
mod.output_model_stats();
//mod.output_fnz();
......
......@@ -108,7 +108,7 @@ void Model::populate_grid(){
// Loop through sample points
y = y0; hh = 1e-16; // Make sure we reset hh
for(int j=0; j<jmax; j++){
for(int j=0; j<jmax+2; j++){
// Integrate until the next sample point
while (y < g_y[j]){
......@@ -143,10 +143,7 @@ void Model::populate_grid(){
void Model::calculate_a2grid(){
/// Calculate the a2 values on the grid now that we have a1, a1dot values
double x, y, r, t; int jmax;
double dr = XSTEP*rs; double invdr = 1./(2.*dr);
// TESTING
double a11, a12, a13, a14, a15, a16, a17, a18;
double dr, invdr;
// Loop through x values
for(int i=0; i<XCELLS+2; i++){
......@@ -154,24 +151,20 @@ void Model::calculate_a2grid(){
x = g_x[i]; r = x*rs;
jmax = get_max_j(x); // Calculate the maximum j value to go up to
dr = rs * (x - g_x[i-1]);
invdr = 1./dr;
// FIXME: SHould this be 1/dr or 1/2dr? 1/dr gives the right-looking answer! But 2dr should be correct
// Create new temporary time slice vectors
// Add an extra 2 elements to each slice to avoid boundary problems when interpolating
vector<double> tmp_a2(jmax+2, 0.0);
vector<double> tmp_a2dot(jmax+2, 0.0);
for(int j=0; j<jmax; j++){
for(int j=0; j<jmax+2; j++){
y = g_y[j] + xTb(x);
t = y*ts;
if(i>1){
a11 = a1(r, t);
a12 = a1(r+dr, t);
a13 = a1(r-dr, t);
a15 = a1(r-2.*dr, t);
a16 = a1(r+2.*dr, t);
a17 = a1(r+3.*dr, t);
a14 = r*(a1(r+dr, t) - a1(r-dr, t))*invdr;
tmp_a2[j] = a11 + a14;
//tmp_a2[j] = a1(r, t) + r*(a1(r+dr, t) - a1(r-dr, t))*invdr;
tmp_a2[j] = a1(r, t) + r*(a1(r+dr, t) - a1(r-dr, t))*invdr;
tmp_a2dot[j] = a1dot(r, t) + r*(a1dot(r+dr, t) - a1dot(r-dr, t))*invdr;
if(tmp_a2[j] <= 0.0){cerr << "WARNING: Shell crossing occurred at r=" << x*rs << ", y/y_today=" << y/y_today << " (NumericGrid::populate_grid)" << endl; throw 99;} // end check
}else{
......
......@@ -14,18 +14,10 @@ double Model::a2dot(double r, double t){
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;
double x = r/rs;
double y = t/ts - xTb(x);
double gxsize = g_x.size();
double gysize = g_y.size();
int Model::get_x_index(double x){
/// Get the i index for a given x value
int a, b, mid;
int gxsize = g_x.size();
// Sanity checks on x, y
if(x < g_x[0]){
......@@ -35,18 +27,12 @@ double Model::get_array_value(double r, double t, int v){
}else{
cerr << "ERROR: x=" << x << " outside calculated range: " << g_x[0] << " to " << g_x[gxsize-1] << 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: t=" << y*ts << " outside calculated range t0=" << ts*y0 << ", t_today=" << ts*y_today << ", y-y_today=" << y-y_today << endl; throw 72;
}
} // end sanity check
// Check if x is too big, and if so, get the asymptotic value
if(x >= g_x[gxsize-1]){
if(x >= g_x[gxsize-1]){
x = g_x[gxsize-1]; // Set x to value where expn. vars. reach asymptotic values, it'll return the same results
}
}
// Get i<=>x by searching g_x
a = 0; b = gxsize - 3; // Since XCELLS = xsize - 2
......@@ -58,7 +44,21 @@ double Model::get_array_value(double r, double t, int v){
a = mid; // mid < x <= b
} // end range check
} // end x-searching loop
i = a;
return a;
}
int Model::get_y_index(double y){
/// Get the j index for a given y value
int a, b, mid;
int gysize = g_y.size();
// 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: t=" << y*ts << " outside calculated range t0=" << ts*y0 << ", t_today=" << ts*y_today << ", y-y_today=" << y-y_today << endl; throw 72;
} // end sanity check
// Get j<=>y by searching g_y
a = 0; b = gysize - 1;
......@@ -70,11 +70,26 @@ double Model::get_array_value(double r, double t, int v){
a = mid; // mid < y <= b
} // end range check
} // end y-searching loop
j = a;
return 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
}
double Model::get_array_value(double r, double t, int v){
/// Generic interpolation routine
int i, j;
double x1, x2, y1, y2;
double xcorr = 0.0; double ycorr = 0.0;
double x = r/rs;
double y = t/ts - xTb(x);
double gxsize = g_x.size(); double gysize = g_y.size();
i = get_x_index(x);
j = get_y_index(y);
// Check if x, y indices are safe to do forward interpolation
if(i >= gxsize and j < gysize){
......@@ -92,11 +107,8 @@ double Model::get_array_value(double r, double t, int v){
x1 = g_x[i]; x2 = g_x[i+1];
y1 = g_y[j]; y2 = g_y[j+1];
//cerr << "x1=" << x1 << ", x2=" << x2 << ", y1=" << y1 << ", y2=" << y2 << endl;
//cerr << "x=" << x << ", y=" << y << endl;
// Do a simple (bilinear) interpolation
// FIXME: This is wrong! May not be [i, j+1]
// FIXME: This is wrong! May not be [i, j+1] for a shifted grid
double arr1, arr2, arr3, arr4;
switch(v){
case 0: arr1 = g_a1[i][j]; arr2 = g_a1[i][j+1];
......@@ -111,24 +123,7 @@ double Model::get_array_value(double r, double t, int v){
break;
}
/*
cerr << "-------------\nV = " << v << endl;
cerr << " ca[i,j]=" << arr1 * (ycorr + y2 - y)*(xcorr + x2 - x)
<< ", ca[i,j+1]=" << arr2 * (ycorr + y - y1)*(xcorr + x2 - x)
<< ", ca[i+1,j]=" << arr3 * (ycorr + y2 - y)*(xcorr + x - x1)
<< ", ca[i+1,j+1]=" << arr4 * (ycorr + y - y1)*(xcorr + x - x1)
<< endl;
cerr << " a[i,j]=" << arr1
<< ", a[i,j+1]=" << arr2
<< ", a[i+1,j]=" << arr3
<< ", a[i+1,j+1]=" << arr4
<< endl;
cerr << "(ycorr + y2 - y)=" << (ycorr + y2 - y)
<< ", (xcorr + x2 - x)=" << (xcorr + x2 - x) << endl;
cerr << " ycorr=" << ycorr << ", xcorr=" << xcorr << endl;
*/
// Return linear (grid) interpolation value
return (1./((y2 - y1)*(x2 - x1))) * (
arr1 * (ycorr + y2 - y)*(xcorr + x2 - x)
+ arr2 * (ycorr + y - y1)*(xcorr + x2 - x)
......
......@@ -22,6 +22,8 @@ protected:
void generate_geodesic();
void hubble_time();
int get_x_index(double x);
int get_y_index(double y);
int get_max_j(double x);
void populate_timesample_vector();
void calculate_a2grid();
......
......@@ -103,7 +103,7 @@ void Model::output_velocities(){
// Get the r when z~5
//double rmax = fn_r(fn_tz(4.0));
double rmax = fn_r(fn_tz(0.2));
double rmax = fn_r(fn_tz(4.0));
// Output things until then //rmax/50.
for(double r=0.0; r<rmax; r+=(rmax/50.)){
......
......@@ -2,7 +2,7 @@
void Model::prepare_grid(){
/// Create new grids
YCELLS = 300;
YCELLS = 300; // FIXME: Normally 300
H0 = fH0 * h;
y0 = T_REC/ts; // Initial time
......@@ -13,7 +13,7 @@ void Model::prepare_grid(){
XMIN = 20; // The minimum x that should be reached before we start considering whether it has reached asymp. value
XSTEP = 0.01; // The step to take in x when populating the grid
XMAX = 400; // FIXME: XMAX = 200 normally
XMAX = 50; // FIXME: XMAX = 200 normally // And then had it set to 400
XCELLS = XMAX / XSTEP; // Maximum number of x cells to calculate
// Get the time when the Hubble constant is reached locally (x=0, y=now)
......
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