Commit d6f099d9 authored by Phil Bull's avatar Phil Bull

Add code to handle the distinction between global and local times, for when...

Add code to handle the distinction between global and local times, for when there is an inhomogeneous bang time
parent 3a8f8952
......@@ -20,16 +20,17 @@ 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
GBH mod(0.66, 0.13, 2.5e3, 0.64*2.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();
......
......@@ -13,6 +13,10 @@
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
// Boolean values
const bool TIME_LOCAL = false;
const bool TIME_GLOBAL = true;
// 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
......@@ -20,8 +24,8 @@ 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 = 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
......
......@@ -20,7 +20,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 +41,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);
......
......@@ -54,11 +54,9 @@ 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: Basic (fixed) step-size definition for now
TBMIN = 0.0; // FIXME: Should set TBMIN (minimum value of Bang-time fn) in the defined models
// FIXME: Not Bang-time-safe!
double yfinestep = ((30.0/ts) - (TBMIN/ts))/(100.0); // Take 50 samples between TBMIN and t=50Myr
double ystep = (y_today - (30.0/ts))/((double)(YCELLS-100));
// FIXME: Not Bang-time-safe?
double yfinestep = 10.0/(ts*200.0); // Take 200 samples between 0 and t=50Myr
double ystep = (y_today - (10.0/ts))/((double)(YCELLS-200));
/*
// Make YCELLS no. of samples
......@@ -68,7 +66,7 @@ void Model::populate_timesample_vector(){
*/
int i; double y;
for(i=0; i<YCELLS; i++){
if(y<50.0/ts){
if(y<10.0/ts){
y+=yfinestep;
}else{
y+=ystep;
......@@ -155,7 +153,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 x, y, r, tl, tg; int jmax;
double drm, invdr;
int xsize = g_x.size();
......@@ -178,16 +176,20 @@ void Model::calculate_a2grid(){
vector<double> tmp_a2dot(jmax+2, 0.0);
for(int j=0; j<jmax+2; j++){
y = g_y[j] + xTb(x);
t = y*ts;
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, 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;
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, t, 0);
tmp_a2dot[j] = a1dot(r, t, 0);
tmp_a2[j] = a1(r, tl, TIME_LOCAL);
tmp_a2dot[j] = a1dot(r, tl, TIME_LOCAL);
}
} // end j loop y values
......
double Model::a1(double r, double t, int i){
return get_array_value(r, t, 0, i);
double Model::a1(double r, double t, bool tt){
return get_array_value(r, t, 0, tt);
}
double Model::a1dot(double r, double t, int i){
double Model::a1dot(double r, double t, bool tt){
/// Units: inverse time
return get_array_value(r, t, 1, i);
return get_array_value(r, t, 1, tt);
}
double Model::a2(double r, double t, int i){
return get_array_value(r, t, 2, i);
double Model::a2(double r, double t, bool tt){
return get_array_value(r, t, 2, tt);
}
double Model::a2dot(double r, double t, int i){
double Model::a2dot(double r, double t, bool tt){
/// Units: inverse time
return get_array_value(r, t, 3, i);
return get_array_value(r, t, 3, tt);
}
int Model::get_x_index(double x){
......@@ -51,6 +51,7 @@ int Model::get_x_index(double x){
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();
......@@ -78,18 +79,25 @@ 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, int i=-1){
double Model::get_array_value(double r, double t, int v, bool time_type = TIME_GLOBAL){
/// Generic interpolation routine
int j; int tmpi; tmpi = i; bool happened = false;
int i, j; double y;
double x1, x2, y1, y2;
double xcorr = 0.0; double ycorr = 0.0;
double x = r/rs;
double y = t/ts - xTb(x);
if(time_type == TIME_GLOBAL){
// Evaluate this as being a global time (t in global time coordinate)
// a1, a1dot, ... arrays stored in local time, so need to correct for Bang time
y = t/ts - xTb(x);
}else{
// Evaluate this as being a local time (t since local Big Bang)
y = t/ts;
}
double gxsize = g_x.size(); double gysize = g_y.size();
if(i==-1){i = get_x_index(x);} // If i=-1, it's the default value, so need to find i ourselves
i = get_x_index(x);
j = get_y_index(y);
// Check if x, y indices are safe to do forward interpolation
......
......@@ -30,12 +30,12 @@ protected:
double z(double r);
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 a1(double r, double t, bool tt = TIME_GLOBAL);
double a1dot(double r, double t, bool tt = TIME_GLOBAL);
double a2(double r, double t, bool tt = TIME_GLOBAL);
double a2dot(double r, double t, bool tt = TIME_GLOBAL);
double get_array_value(double r, double t, int v, int i);
double get_array_value(double r, double t, int v, bool time_type);
double xM(double x);
double xMprime(double x);
......
......@@ -6,7 +6,7 @@ void Model::output_background(double t){
ofstream bk; bk.open(fname.str().c_str());
double loc_dens;
for(double r=r0; r<6000.; r += 2.){
for(double r=r0; r<20000.; r += 2.){
loc_dens = (((8.*M_PI*G/3.) * m(r) / a1(r, t)) - C*C*k(r));
bk << r << "\t"
<< density(r, t) << "\t"
......
......@@ -2,7 +2,7 @@
void Model::prepare_grid(){
/// Create new grids
YCELLS = 400; // FIXME: Normally 300
YCELLS = 500; // FIXME: Normally 300
H0 = fH0 * h;
y0 = T_REC/ts; // Initial time
......@@ -12,8 +12,8 @@ void Model::prepare_grid(){
t0 = y0 * ts;
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 = 50;//50; // FIXME: XMAX = 200 normally // And then had it set to 400
XSTEP = 0.01; // FIXME: NORMALLY 0.01 // The step to take in x when populating the grid
XMAX = 200;//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)
......
......@@ -15,7 +15,7 @@ double Model::outgoing_redshift(double r_start, double t_start, double t_end){
//ofstream gd; gd.open(TEST_FILENAME.c_str(), ios_base::app);
//double tsamp = t_start*2.; // Large starting value, so we sample in the first instance
for(double t=t_start; t > t_end - dt; t+=dt){
for(double t=t_start; t - tb(r) > (t_end - dt); t+=dt){
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);
......@@ -56,7 +56,7 @@ double Model::ingoing_redshift(double r_start, double t_start, double t_end, dou
//ofstream gd; gd.open(TEST_FILENAME.c_str(), ios_base::app);
//double tsamp = t_start*2.; // Large starting value, so we sample in the first instance
for(double t=t_start; t > t_end - dt; t+=dt){
for(double t=t_start; t - tb(r) > (t_end - dt); t+=dt){
k1 = 1. * sqrt(1. - k(fabs(r))*r*r) / a2(fabs(r), t);
k2 = 1. * sqrt(1. - k(fabs(r + 0.5*dtr*k1))*(r + 0.5*dtr*k1)*(r + 0.5*dtr*k1))
/ a2(fabs(r + 0.5*dtr*k1), t + 0.5*dt);
......@@ -88,6 +88,9 @@ double Model::v_dipole(double r, double* r_end = NULL){
/// Find the dipole velocity (derived from the redshift) at radius r
/// (in km/sec)
// FIXME: This doesn't deal with varying Bang time properly
// Should use endpoint as the local time after Big Bang?
// Find the local time at which point r is observed by a central
// observer today
double t_obs = fn_t(r);
......@@ -102,7 +105,10 @@ double Model::v_dipole(double r, double* r_end = NULL){
// 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(590.0); // FIXME: z=40
// *Local* time at decoupling
double t_global_dec = fn_tz(72.0);
double t_dec = t_global_dec - tb(fn_r(t_global_dec)); // FIXME: z=40
cerr << "t_global_dec = " << t_global_dec << ", t_dec(local) = " << t_dec << endl;
// Find the redshift for the ray going *out* of the void
// FIXME: Can be replaced by the pre-calculated geodesic
......
......@@ -46,6 +46,17 @@ public:
/// Radial derivative of curvature-like function, k'(r)
return -3. * KR * A * r*r * exp(-KR * r*r*r);
}
/*
double tb(double r){
/// Bang-time function, Myr
return 50.0*exp(-1.5*KR*(r-1500.)*(r-1500.));
}
double tbprime(double r){
/// Radial derivative of Bang-time function
return -50.0 * 2. * 1.5 * KR * (r-1500.) * exp(-1.5*KR*(r-1500.)*(r-1500.));
}*/
private:
......
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