Commit 3a8f8952 authored by Phil Bull's avatar Phil Bull

get_array_value testing, improved (faster) index look-up for a2 building,...

get_array_value testing, improved (faster) index look-up for a2 building, output inbound geodesic endpoint and plot it
parent 25be3eec
......@@ -20,19 +20,19 @@ int main(int argc, char* argv[]){
/// 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.7, 1.5e3); // (h, A, w)
CFLb mod(0.72, 0.7, 1.5e3); // (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, 1e3, 0.5e3); // (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
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_model_stats();
//mod.output_fnz();
//mod.output_zin();
......
......@@ -75,7 +75,6 @@ void Model::populate_timesample_vector(){
}
g_y.push_back( y );
}
cerr << "$$$ i = " << i << ", YCELLS = " << YCELLS << endl;
}
void Model::populate_grid(){
......@@ -157,7 +156,8 @@ 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, invdr;
double drm, invdr;
int xsize = g_x.size();
// Loop through x values
for(int i=0; i<XCELLS+2; i++){
......@@ -165,8 +165,10 @@ 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 = 0.5/dr;
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
......@@ -178,14 +180,14 @@ void Model::calculate_a2grid(){
for(int j=0; j<jmax+2; j++){
y = g_y[j] + xTb(x);
t = y*ts;
if(i>1){
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(i>0 and i<xsize-1){
tmp_a2[j] = a1(r, t, i) + r*(a1(r, t, i) - a1(r-drm, t, i-1))*invdr;
tmp_a2dot[j] = a1dot(r, t, i) + r*(a1dot(r, t, i) - a1dot(r-drm, t, i-1))*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{
// a2, a2dot = a1, a1dot at r=0
tmp_a2[j] = a1(r, t);
tmp_a2dot[j] = a1dot(r, t);
tmp_a2[j] = a1(r, t, 0);
tmp_a2dot[j] = a1dot(r, t, 0);
}
} // end j loop y values
......
double Model::a1(double r, double t){
return get_array_value(r, t, 0);
double Model::a1(double r, double t, int i){
return get_array_value(r, t, 0, i);
}
double Model::a1dot(double r, double t){
double Model::a1dot(double r, double t, int i){
/// Units: inverse time
return get_array_value(r, t, 1);
return get_array_value(r, t, 1, i);
}
double Model::a2(double r, double t){
return get_array_value(r, t, 2);
double Model::a2(double r, double t, int i){
return get_array_value(r, t, 2, i);
}
double Model::a2dot(double r, double t){
double Model::a2dot(double r, double t, int i){
/// Units: inverse time
return get_array_value(r, t, 3);
return get_array_value(r, t, 3, i);
}
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();
......@@ -37,8 +38,8 @@ int Model::get_x_index(double x){
// 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
......@@ -63,8 +64,8 @@ int Model::get_y_index(double y){
// Get j<=>y by searching g_y
a = 0; b = gysize - 1;
while( (b - a) > 1 ){
mid = a + floor(0.5 * (double)(b - a)); // Find the middle cell
if(y <= g_y[mid]){
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
......@@ -77,9 +78,9 @@ int Model::get_y_index(double y){
//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){
double Model::get_array_value(double r, double t, int v, int i=-1){
/// Generic interpolation routine
int i, j;
int j; int tmpi; tmpi = i; bool happened = false;
double x1, x2, y1, y2;
double xcorr = 0.0; double ycorr = 0.0;
......@@ -88,7 +89,7 @@ double Model::get_array_value(double r, double t, int v){
double gxsize = g_x.size(); double gysize = g_y.size();
i = get_x_index(x);
if(i==-1){i = get_x_index(x);} // If i=-1, it's the default value, so need to find i ourselves
j = get_y_index(y);
// Check if x, y indices are safe to do forward interpolation
......
......@@ -30,11 +30,12 @@ protected:
double z(double r);
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 a1(double r, double t, int i=-1);
double a1dot(double r, double t, int i=-1);
double a2(double r, double t, int i=-1);
double a2dot(double r, double t, int i=-1);
double get_array_value(double r, double t, int v, int i);
double xM(double x);
double xMprime(double x);
......@@ -97,9 +98,9 @@ public:
double vp_frw(double r);
//double v_dipole(double r, double ZLIM); // Use this definition if you're doing integrand output for off-centre obs.
double v_dipole(double r);
double v_dipole(double r, double* r_end);
double outgoing_redshift(double r_start, double t_start, double t_end);
double ingoing_redshift(double r_start, double t_start, double t_end);
double ingoing_redshift(double r_start, double t_start, double t_end, double* r_end);
void output_background(double t);
void output_background();
......@@ -110,7 +111,7 @@ public:
void output_velocity_mismatch();
void output_kszvpec_for_z();
void output_fnz();
void output_zin();
//void output_zin();
Model(){
// Constructor. Don't do anything here yet. Do it all in the derived class.
......
......@@ -104,15 +104,17 @@ void Model::output_velocities(){
// Get the r when z~5
//double rmax = fn_r(fn_tz(4.0));
double rmax = fn_r(fn_tz(10.0));
double r_end;
// Output things until then //rmax/50.
for(double r=0.0; r<rmax; r+=(rmax/75.)){
double tt = fn_t(r);
vel << r << "\t"
<< z(r) << "\t"
<< v_dipole(r) << "\t"
<< v_dipole(r, &r_end) << "\t"
<< (Hr(r+20.0, tt) - Hr(r, tt))/20.0 << "\t"
<< vp_frw(r)
<< vp_frw(r) << "\t"
<< r_end
<< endl;
}
vel.close();
......@@ -178,6 +180,7 @@ void Model::output_fnz(){
}
/*
void Model::output_zin(){
/// Output z_in calculated for off-centre observers
cerr << "\toutput_zin()" << endl;
......@@ -201,7 +204,7 @@ void Model::output_zin(){
}
vel.close();
cerr << "\toutput_velocities(): " << ((double)clock() - (double)tclock)/CLOCKS_PER_SEC << " seconds." << endl;
}
}*/
void Model::output_kszvpec_for_z(){
/// Output the expected kSZ velocity given a set of z values
......
......@@ -42,7 +42,7 @@ double Model::outgoing_redshift(double r_start, double t_start, double t_end){
return exp(I * 0.5 * dt) - 1.; // Redshift, z
}
double Model::ingoing_redshift(double r_start, double t_start, double t_end){
double Model::ingoing_redshift(double r_start, double t_start, double t_end, double* r_end){
/// Return the redshift for rays going out of the void
double k1, k2, k3, k4;
double dtr = dt * C; // dt in distance units
......@@ -80,11 +80,11 @@ double Model::ingoing_redshift(double r_start, double t_start, double t_end){
} // end integration loop
//gd.close();
*r_end = r;
return exp(I * 0.5 * dt) - 1.; // Redshift, z. Factor of 0.5*dt comes from simple integration
}
double Model::v_dipole(double r){
double Model::v_dipole(double r, double* r_end = NULL){
/// Find the dipole velocity (derived from the redshift) at radius r
/// (in km/sec)
......@@ -102,16 +102,14 @@ double Model::v_dipole(double r){
// Find the time at which the redshift is z=40
// (proxy for z=1100, decoupling)
// FIXME: Needs to go to z=1100
double t_dec = fn_tz(600.0); // FIXME: z=40
double t_dec = fn_tz(590.0); // FIXME: z=40
// Find the redshift for the ray going *out* of the void
// FIXME: Can be replaced by the pre-calculated geodesic
double z_out = outgoing_redshift(r, t_obs, t_dec);
// Find the redshift for the ray going *through* the void
double z_in = ingoing_redshift(r, t_obs, t_dec);
//cerr << "r=" << r << ", z_in=" << z_in << ", z_out=" << z_out << endl;
double z_in = ingoing_redshift(r, t_obs, t_dec, r_end);
// GBH give dT/T|Dipole = (z_in - z_out) / (2 + z_in + z_out)
return 3e5 * (z_in - z_out) / (2. + z_in + z_out);
......
......@@ -30,6 +30,7 @@ gdmmilne = []
rr = []
zz = []
dz = []
rend = []
# Supernova data from http://www.astro.ucla.edu/~wright/sne_cosmology.html [Kowalski]
# z, deltaDM, sigma
......@@ -122,6 +123,7 @@ for line in h.readlines():
rr.append(float(tmp[0]))
zz.append(float(tmp[1]))
dz.append(float(tmp[2]))
rend.append(float(tmp[5]))
h.close()
......@@ -172,8 +174,9 @@ title('r(z)')"""
#ylim((0.0, 3000.0))
subplot(235)
p = plot(zz, dz, 'r+', zz, dz, 'r-')
p = plot(zz, rend, 'b-')
grid(True)
title('deltaZ(z)')
title('delta_vp(z)[r], r_end[b]')
#xlim((0.0, 6000.0))
# Distance modulus
......
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