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::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, tt); } double Model::a2(double r, double t, bool tt){ return get_array_value(r, t, 2, tt); } double Model::a2dot(double r, double t, bool tt){ /// Units: inverse time return get_array_value(r, t, 3, tt); } 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[0]){ if(x < 0.0){ //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[0] << " to " << g_x[gxsize-1] << endl; throw 71; } // end check whether x is +ve } // end sanity check // Check if x is too big, and if so, get the asymptotic value 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 + 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 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 } 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 = 0.0; double x = r/rs; 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(); 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){ // 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 >= gysize){ ycorr = g_y[j] - g_y[j-1]; j -= 1; // We're at the y boundary }else if(i >= gxsize and j >= gysize){ xcorr = g_x[i] - g_x[i-1]; ycorr = g_y[j] - g_y[j-1]; 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 = g_y[j]; y2 = g_y[j+1]; // Do a simple (bilinear) interpolation // 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]; 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 linear (grid) interpolation value return (1./((y2 - y1)*(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) ); }