Commit f1ed288b authored by Phil Bull's avatar Phil Bull

Merge to master

parents 76a4faa0 383284bd
......@@ -5,3 +5,5 @@ twofluid
data/*
doxygen/*
figures/*
callgrind*
*.pyc
#!/usr/bin/python
class obskSZ(object):
"""A collection of kSZ data points"""
def load_ksz_data(self):
"""Load the kSZ data from a file"""
f = open(self.filename, 'r')
for line in f.readlines():
line = line.split("\t")
if len(line) > 4 and line[0]!="#":
self.z.append(float(line[1]))
self.vp.append(float(line[2]))
self.e_upper.append(float(line[3]))
self.e_lower.append(float(line[4]))
def __init__(self):
self.filename = "obsdata/GBH_kSZ_data.dat" # The file to retrieve the data from
self.z = [] # The z value taken for this measurement
self.vp = [] # The measured peculiar velocity, in km/sec
self.e_upper = [] # The upper error on vp
self.e_lower = [] # The lower error on vp
self.load_ksz_data() # Load the data from the file
def get_closest(z, zvec):
"""Get the data point from the model (zvec) with the closest z value
to some real kSZ measurement (z). Returns the index in
zvec which comes closest. Assumes that zvec is sorted."""
a = 0
b = len(zvec) - 1
while (b-a)>1:
mid = a + int(0.5*float(b-a))
if zvec[mid] > z:
b = mid
else:
a = mid
return a
def chisquared(model_z, model_vp):
"""Return the chi-squared statistic for the kSZ vp(z) data, and the no. of data points"""
# model_z and model_kSZ are vectors of data output by the model
obs = obskSZ() # Load the SNe data from the file
x2 = 0.0 # The Chi-squared statistic
count = 0
for j in range(len(obs.vp)):
count += 1
i = get_closest(obs.z[j], model_z) # Get the index of the closest of the model-output data points
# FIXME: This isn't right, because we're forgetting the systematic error, and the data clearly isn't Gaussian - lopsided error bars!
x2 += (obs.vp[j] - model_vp[i])**2.0 / (obs.e_lower[j]**2.0) # Calculate the contribution to the X^2 sum for this data point
return x2, count
#!/usr/bin/python
class obsDM(object):
"""A collection of supernova/GRB data points"""
def load_sn_data(self):
"""Load the SN/GRB data from a file"""
f = open(self.filename, 'r')
for line in f.readlines():
line = line.split(" ")
if len(line) > 5 and line[0]!="#":
self.zmin.append(float(line[0]))
self.zmax.append(float(line[1]))
self.z.append(float(line[2]))
self.DM.append(float(line[3]))
self.errDM.append(float(line[4]))
self.source.append(str(line[5][:-1]))
def __init__(self):
self.filename = "obsdata/SNe_Wright.dat" # The file to retrieve the data from
self.z = [] # The mean z value taken for this measurement
self.zmin = [] # The minimum z bound for this measuremnet
self.zmax = [] # The maximum z bound for this measuremnet
self.DM = [] # The measured distance modulus
self.errDM = [] # The error on the distance modulus
self.source = [] # A string denoting the source of the data
self.load_sn_data()
def get_closest(z, zvec):
"""Get the data point from the model (zvec) with the closest z value
to some real deltaDM measurement (z). Returns the index in
zvec which comes closest. Assumes that zvec is sorted."""
a = 0
b = len(zvec) - 1
while (b-a)>1:
mid = a + int(0.5*float(b-a))
if zvec[mid] > z:
b = mid
else:
a = mid
return a
def chisquared(model_z, model_DM):
"""Return the chi-squared statistic for the SNe data, and the no. of data points"""
# model_z and model_DM are vectors of data output by the model
obs = obsDM() # Load the SNe data from the file
x2 = 0.0 # The Chi-squared statistic
count = 0
for j in range(len(obs.DM)):
#if obs.source[j][:3] != "GRB":
count += 1
i = get_closest(obs.z[j], model_z) # Get the index of the closest of the model-output data points
x2 += (obs.DM[j] - model_DM[i])**2.0 / (obs.errDM[j]**2.0) # Calculate the contribution to the X^2 sum for this data point
return x2, count
......@@ -10,31 +10,76 @@ int main(int argc, char* argv[]){
/// Main function
clock_t tclock = clock();
//if(argc<2){cerr << "ERROR: Not enough command-line arguments." << endl; return 1;}
//double rrr = atof(argv[1]);
//if(argc<4){cerr << "ERROR: Not enough command-line arguments." << endl; return 1;}
//double hh = atof(argv[1]);
//double AA = atof(argv[2]);
//double ww = atof(argv[3]);
// Instantiate a model
/// Interesting one: CFLb mod(0.72, -5.75e-8, 1.0e-9); // (h, A, KR)
/// Another interesting one, with cusp!: CFLb mod(0.72, 0.5, 1.5e3); // (h, A, w)
//CFLa mod(0.72, 0.85, 1e3); // (h, A, w)
//CFLb mod(0.72, 0.5, 1.5e3); // (h, A, w) // Good one!
CFLb mod(0.72, 0.5, 1.5e3); // (h, A, w)
//CFLb mod(0.72, 0.7, 1.5e3); // (h, A, w) // Has been the standard one for a while
CFLb mod(0.72, 0.8, 1.0e3); // (h, A, w)
////CFLb mod(hh, AA, ww); // (h, A, w)
//CFLb mod(0.72, 0.1, 1.5e4); // (h, A, w)
//FRW mod(0.72, 0.0001, 0.9999); // (h, omega_m, omega_k) almost-Milne
//FRW mod(0.72, 0.3, 0.6); // (h, omega_m, omega_k)
///GBH mod(0.65, 0.13, 2e3, 1e3); // (h, omega_in, r0, deltar) Garcia-Bellido and Haugbolle
//GBH mod(0.65, 0.13, 1e3, 0.5e3); // (h, omega_in, r0, deltar) Garcia-Bellido and Haugbolle
///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_contour();
//mod.output_velocities();
mod.output_velocities();
//mod.output_velocity_mismatch();
mod.output_model_stats();
//mod.output_kszvpec_for_z();
//mod.output_model_stats();
//mod.output_fnz();
//mod.output_zin();
// TESTING
//mod.v_dipole(rrr);
//cerr << "ADAPTIVE v_p = " << mod.adaptive_v_dipole(100.0) << endl;
//cerr << "NORMAL v_p = " << mod.v_dipole(100.0) << endl;
// TESTING
/*
cerr << "IN t(r) = " << mod.adaptive_t_fnr(1.0, 1000.0, mod.fn_t(1.0), GEODESIC_INGOING) << endl;
cerr << endl;
cerr << "OUT t(r) = " << mod.adaptive_t_fnr(1.0, 1000.0, mod.fn_t(1.0), GEODESIC_OUTGOING) << endl;
cerr << endl;
*/
/*cerr << "r\tAD\tNORM\tr(z)" << endl;
for(double r=0.0; r<2000.0; r+=50.0){
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;
}*/
/*
cerr << "g/PI\tz" << endl;
double tend = mod.fn_tz(1.0);
double rstart = 100.0;
/*
for(double g=0.0; g<=2.*M_PI; g+=2.*M_PI/20.0){
cerr << g/M_PI << "\t" << mod.general_offcentre(rstart, mod.t_today, g, tend) << endl;
} // end loop through gamma angles
//cerr << 0 << "\t" << mod.general_offcentre(rstart, mod.t_today, 0.0, tend) << endl;
cerr << 90 << "\t" << mod.general_offcentre(rstart, mod.t_today, M_PI*0.5, tend) << endl;
//cerr << 180 << "\t" << mod.general_offcentre(rstart, mod.t_today, M_PI, tend) << endl;
*/
// Finish timing the run
cerr << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
......
......@@ -11,16 +11,23 @@
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_integration.h>
// Boolean values
const bool TIME_LOCAL = false;
const bool TIME_GLOBAL = true;
const bool GEODESIC_INGOING = true;
const bool GEODESIC_OUTGOING = false;
// 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 = 9.17e-4; // FRW scale factor at recombination
const double T_REC = 0.380; // Age of universe at recombination (Myr)
//const double A1_REC = 3.13e-4; // FRW scale factor at matter-radiation equality
//const double T_REC = 0.0766; // Age of universe at matter-radiation equality (Myr)
//const double A1_REC = 9.17e-4; // FRW scale factor at recombination
//const double T_REC = 0.380; // Age of universe at recombination (Myr)
const double A1_REC = 1e-7; // FRW scale factor at an early time
const double T_REC = 1e-7; // Age of universe at an early time (but actually, made up) (Myr)
const double A_MAX = 1.0; // The maximum a1 value to calculate
......@@ -40,6 +47,11 @@ const double rs = 100.0; // 100 MPC
#include "model_def.cpp"
#include "model_output.cpp"
// TESTING
#include "model_integrate_geodesic.cpp"
#include "offcentre_general.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
......@@ -53,23 +65,3 @@ 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();
*/
void Model::generate_geodesic(){
/// Build a geodesic by integrating along dt to find (r,t) at
/// multiple sample points
/// Also create a geodesic sample
/// USES ADAPTIVE INTEGRATOR
clock_t tclock = clock();
// Set initial values
double tnext = 0.0; double rcell = 10.0;
// Add initial values to the r,t arrays
double t = t_today;
double r = 0.0;
double I = 0.0;
double zlast = 0.0;
double zz = 0.0;
rval.push_back(r); tval.push_back(t); zval.push_back(0.0);
init_integrators(); // Initialise the GSL integration routines
// Loop through t cells
while (t > t0 + 0.01 + tb(r) and zz < 1200.0){
// Adaptive step to find t(r) and integrand of z, I
tnext = adaptive_t_fnr(r, r+rcell, t, GEODESIC_OUTGOING);
I += adaptive_I(r, r+rcell, t);
// Increment values
r += rcell;
t = tnext;
zz = exp(I) - 1.0;
// Check if integration step size could be changed
if(zz - zlast > 10.0){
rcell *= 0.1;
}else if(zz - zlast < 0.001){
rcell *= 10.0;
}
zlast = zz;
// Add values to array
rval.push_back(r);
tval.push_back(t);
zval.push_back(zz);
}
free_integrators(); // Free the memory used by the GSL integration/ODE solving routines
// TESTING
cerr << "\tgenerate_geodesic(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
}
/*
void Model::generate_geodesic(){
/// Build a geodesic by integrating along dt to find (r,t) at
/// multiple sample points
......@@ -20,7 +73,7 @@ void Model::generate_geodesic(){
// FIXME: Needs better high-z behaviour, very inaccurate ATM
// Have to stop a step early, because k4 = k4(t + dt2), so need t + dt2 > t_end
for(double t=t_today; t > t0-dt2; t+=dt2){
for(double t=t_today; t > (t0-dt2) + tb(r); t+=dt2){
// Calculate Runge-Kutta variables
k1 = -1. * sqrt(1. - k(r)*r*r) / a2(r, t);
......@@ -41,8 +94,12 @@ void Model::generate_geodesic(){
I += dI + dIlast;
dIlast = dI;
// Add the next sample, if it's time to do so
if((tlast - t) > tcell2){
//cerr << dI << endl;
if(dI>1e-5 and -dt2 > (t_today - t0)/1e10){dt2 *= 0.1; dtr = dt2 * C;} // Adaptively shrink dt2
tlast = t;
rval.push_back(r);
tval.push_back(t);
......@@ -57,5 +114,6 @@ void Model::generate_geodesic(){
// TESTING
cerr << "r=" << rval[rval.size()-1] << ", t=" << tval[rval.size()-1] << ", z=" << zval[rval.size()-1] << endl;
cerr << "\tNo. samples = " << rval.size() << endl;
cerr << "\tgenerate_geodesic(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
}
}*/
......@@ -28,13 +28,65 @@ double Model::da1dot_da(double a, double m, double k){
}
int Model::get_max_j(double x){
/// Get the j value which will ensure that the grid will be
/// populated at least up to the maximum global time for a given x
int mid;
int a = 0;
int b = g_y.size()-1;
// Target y
double y = y_today - xTb(x);
// Get i<=>y by searching g_y
while( (b - a) > 1 ){
mid = a + floor(0.5 * (double)(b - a)); // Find the middle cell
if(y <= g_y[mid]){
b = mid; // a <= y <= mid
}else{
a = mid; // mid < y <= b
} // end range check
} // end y-searching loop
return a; // Index of the maximum y value to calculate for this x
}
void Model::populate_timesample_vector(){
/// Populate the time sampling vector, which gives the local times
/// on the grid. Can use this to do variable time sampling
// FIXME: Not Bang-time-safe?
double yfinestep = 10.0/(ts*100.0); // Take 100 samples between 0 and t=50Myr
double ystep = (y_today - (10.0/ts))/((double)(YCELLS-100));
/*
// Make YCELLS no. of samples
for(int i=0; i<YCELLS; i++){
g_y.push_back( y0 + (double)i*ystep );
} // end loop through cells
*/
int i; double y;
for(i=0; i<YCELLS; i++){
if(y<10.0/ts){
y+=yfinestep;
}else{
y+=ystep;
}
g_y.push_back( y );
}
}
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 x, y, ysamp;
double x, y; int jmax;
bool reached_asymptotic = false;
// Set-up arrays etc.
populate_timesample_vector();
// Initialise the GSL adaptive integrator
//const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd; // RK8 Prince-Dormand adaptive algorithm
const gsl_odeiv_step_type * T = gsl_odeiv_step_rkf45; // Runge-Kutta-Fehlberg (4,5) adaptive algorithm
......@@ -44,39 +96,36 @@ void Model::populate_grid(){
int status; // Variable to hold status of GSL integration
// Make sure we start at the right place
double hh;
double a1[1] = { A1_REC };
double hh; double a1[1] = { A1_REC };
x = x0 - XSTEP;
ModelPointers p;
p.mod = this;
ModelPointers p; p.mod = this;
// Loop through x values
for(int i=0; i<XCELLS+2 and !reached_asymptotic; i++){
// Set-up for integration
ysamp = y0 + YSTEP;
x += XSTEP;
g_x.push_back(x); // Add the x value for this cell to the array
x += XSTEP; g_x.push_back(x); // Add the x value for this cell to the array
p.m = xM(x); p.k = xK(x);
a1[0] = A1_REC;
// Calculate the maximum j value to go up to
jmax = get_max_j(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);
vector<double> tmp_a1(jmax+2, 0.0);
vector<double> tmp_a1dot(jmax+2, 0.0);
// Define system of ODEs (ODE + Jacobian)
gsl_odeiv_system sys = {static_a1dot, a1dot_jac, 1, &p};
// Loop through sample points
y = y0; hh = 1e-16; // Make sure we reset hh
for(int j=0; j<YCELLS+2; j++){
for(int j=0; j<jmax+2; j++){
// Integrate until the next sample point
while (y < ysamp){
status = gsl_odeiv_evolve_apply (e, c, s, &sys, &y, ysamp, &hh, a1);
while (y < g_y[j]){
status = gsl_odeiv_evolve_apply (e, c, s, &sys, &y, g_y[j], &hh, a1);
if (status != GSL_SUCCESS){ cerr << "Model::populate_grid(): Integration failed." << endl; break; }
} // end while loop
if(isnan(a1[0])){throw 111;}
......@@ -84,46 +133,25 @@ void Model::populate_grid(){
// Add to mesh at this sample point
tmp_a1[j] = a1[0];
tmp_a1dot[j] = s_a1dot(a1[0], p.m, p.k)/ts;
if(i>1){
tmp_a2[j] = g_a1[i-1][j] + x*(1./XSTEP)*(a1[0] - g_a1[i-2][j]) ;
tmp_a2dot[j] = ( g_a1dot[i-1][j] + x*(1./XSTEP)*(tmp_a1dot[j] - g_a1dot[i-2][j]) );
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
}
ysamp += YSTEP;
} // end time steps
g_a1.push_back( tmp_a1 );
g_a1dot.push_back( tmp_a1dot );
g_a2.push_back( tmp_a2 );
g_a2dot.push_back( tmp_a2dot );
// FIXME: This is triggered if x>XMIN and the function is only *temporarily* too flat
// This cuts-off the integration way too early for larger voids
// Check if we're in the asymptotic region (nest for efficiency?)
// Here, a1(r) = a1(r+c), for all c > 0, so no point calculating the rest
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;
if(g_a1[i][jmax] == g_a1[i-10][jmax]){
reached_asymptotic = true;
XCELLS = g_a1.size() - 2;
XMAX = g_x[XCELLS-1];
XMAX = x; XCELLS = XMAX / XSTEP;
cerr << "\tReached asymptotic region at r=" << x*rs << endl;
}
} // end asymptote check
} // end check x region
} // 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]);
}
// Get the a2 values now the dust has settled
calculate_a2grid();
// Free all of the integration variables
gsl_odeiv_evolve_free (e);
......@@ -132,3 +160,54 @@ void Model::populate_grid(){
cerr << "\tpopulate_grid(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
} // end initial integration function
void Model::calculate_a2grid(){
/// Calculate the a2 values on the grid now that we have a1, a1dot values
double x, y, r, tl, tg; int jmax;
double drm, invdr;
int xsize = g_x.size();
// Loop through x values
for(int i=0; i<XCELLS+2; i++){
x = g_x[i]; r = x*rs;
jmax = get_max_j(x); // Calculate the maximum j value to go up to
if(i>0 and i<xsize-1){
drm = rs * (x - g_x[i-1]);
invdr = 1./drm;
}
// FIXME: Should this be 1/dr or 1/2dr? 1/dr gives *a* right-looking answer! But 2dr should be correct
// Seems that shell-crossings more likely with 1/dr
// 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+2; j++){
y = g_y[j];
tl = y*ts; // Local time
tg = tl + ts*xTb(x); // Global time
if(i>0 and i<xsize-1){
tmp_a2[j] = a1(r, tl, TIME_LOCAL) + r*invdr
* (a1(r, tl, TIME_LOCAL) - a1(r-drm, tg, TIME_GLOBAL));
tmp_a2dot[j] = a1dot(r, tl, TIME_LOCAL) + r*invdr
* (a1dot(r, tl, TIME_LOCAL) - a1dot(r-drm, tg, TIME_GLOBAL));
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{
// a2, a2dot = a1, a1dot at r=0
tmp_a2[j] = a1(r, tl, TIME_LOCAL);
tmp_a2dot[j] = a1dot(r, tl, TIME_LOCAL);
}
} // end j loop y values
g_a2.push_back( tmp_a2 );
g_a2dot.push_back( tmp_a2dot );
} // end i loop through x values
}
double Model::a1(double r, double t){
return get_array_value(r, t, 0);
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::a1dot(double r, double t){
double Model::a1(double r, double t, bool tt){
return get_array_value(r, t, 0, tt);
}
double Model::a1dot(double r, double t, bool tt){
/// Units: inverse time
return get_array_value(r, t, 1);
return get_array_value(r, t, 1, tt);
}
double Model::a2(double r, double t){
return get_array_value(r, t, 2);
double Model::a2(double r, double t, bool tt){
return get_array_value(r, t, 2, tt);
}
double Model::a2dot(double r, double t){
double Model::a2dot(double r, double t, bool tt){
/// Units: inverse time
return get_array_value(r, t, 3);
return get_array_value(r, t, 3, tt);
}
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;
int Model::get_x_index(double x){
/// Get the i index for a given x value
// FIXME: Returns (i_returned = i_actual - 1) if x exactly== g_x[i_actual]
int a, b, mid;
int gxsize = g_x.size();
// Sanity checks on x, y
if(x < g_x[a]){
if(x < g_x[0]){
if(x < 0.0){
cerr << "WARNING: x went negative, taking fabs(x) [Model::get_array_value]" << endl;
//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;
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[b-1]){
x = g_x[b-1]; // Set x to value where expn. vars. reach asymptotic values, it'll return the same results
}
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
while( (b - a) > 1 ){
mid = a + floor(0.5 * (double)(b - a)); // Find the middle cell
if(x <= g_x[mid]){
mid = a + int(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;
return a;
}
int Model::get_y_index(double y){
/// Get the j index for a given y value
// FIXME: Needs to be made safe for arbitrary tB(r)
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;
while( (b - a) > 1 ){
mid = a + int(0.5 * (double)(b - a)); // Find the middle cell
if(y < g_y[mid]){
b = mid; // a <= y <= mid
}else{
a = mid; // mid < y <= b
} // end range check
} // end y-searching loop
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
//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, bool time_type = TIME_GLOBAL){
/// Generic interpolation routine
int i, j; double y;
double x1, x2, y1, y2;
double xcorr = 0.0; double ycorr =