model_grid_interp.cpp 5.01 KB
Newer Older
1

2 3 4 5 6
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
}

7 8
double Model::a1(double r, double t, bool tt){
	return get_array_value(r, t, 0, tt);
9
}
10
double Model::a1dot(double r, double t, bool tt){
11
	/// Units: inverse time
12
	return get_array_value(r, t, 1, tt);
13
}
14 15
double Model::a2(double r, double t, bool tt){
	return get_array_value(r, t, 2, tt);
16
}
17
double Model::a2dot(double r, double t, bool tt){
18
	/// Units: inverse time
19
	return get_array_value(r, t, 3, tt);
20 21
}

22 23
int Model::get_x_index(double x){
	/// Get the i index for a given x value
24
	// FIXME: Returns (i_returned = i_actual - 1) if x exactly== g_x[i_actual]
25 26
	int a, b, mid;
	int gxsize = g_x.size();
27 28
	
	// Sanity checks on x, y
29
	if(x < g_x[0]){
30
		if(x < 0.0){
31
			//cerr << "WARNING: x went negative, taking fabs(x) [Model::get_array_value]" << endl;
32 33
			x = fabs(x); // Try the +ve value of x
		}else{
34
			cerr << "ERROR: x=" << x << " outside calculated range: " << g_x[0] << " to " << g_x[gxsize-1] << endl; throw 71;
35
		} // end check whether x is +ve
36
	} // end sanity check
37 38
	
	// Check if x is too big, and if so, get the asymptotic value
39
	if(x >= g_x[gxsize-1]){
40
		 x = g_x[gxsize-1]; // Set x to value where expn. vars. reach asymptotic values, it'll return the same results
41
	}
42 43
	
	// Get i<=>x by searching g_x
44
	a = 0; b = gxsize - 3; // Since XCELLS = xsize - 2
45
	while( (b - a) > 1 ){
46 47
		mid = a + int(0.5 * (double)(b - a)); // Find the middle cell
		if(x < g_x[mid]){
48 49 50 51 52
			b = mid; // a <= x <= mid
		}else{
			a = mid; // mid < x <= b
		} // end range check
	} // end x-searching loop
53 54 55 56 57 58
	
	return a;
}

int Model::get_y_index(double y){
	/// Get the j index for a given y value
59
	// FIXME: Needs to be made safe for arbitrary tB(r)
60 61 62 63 64 65 66 67 68
	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
69
	
70 71 72
	// Get j<=>y by searching g_y
	a = 0; b = gysize - 1;
	while( (b - a) > 1 ){
73 74
		mid = a + int(0.5 * (double)(b - a)); // Find the middle cell
		if(y < g_y[mid]){
75 76 77 78 79
			b = mid; // a <= y <= mid
		}else{
			a = mid; // mid < y <= b
		} // end range check
	} // end y-searching loop
80
	return a;
81
	
82
	// Get j<=>y by looking it up from the index (assumes fixed grid cell size in y)
83 84
	//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
85 86
}

87
double Model::get_array_value(double r, double t, int v, bool time_type = TIME_GLOBAL){
88
	/// Generic interpolation routine
89
	int i, j; double y;
90 91 92 93
	double x1, x2, y1, y2;
	double xcorr = 0.0; double ycorr = 0.0;
	
	double x = r/rs;
94 95 96 97 98 99 100 101
	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;
	}
102 103 104
	
	double gxsize = g_x.size();	double gysize = g_y.size();
	
105
	i = get_x_index(x);
106
	j = get_y_index(y);
107 108
	
	// Check if x, y indices are safe to do forward interpolation
109
	if(i >= gxsize and j < gysize){
110 111 112
		// 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)
113 114 115
	}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){
116
		xcorr = g_x[i] - g_x[i-1];
117
		ycorr = g_y[j] - g_y[j-1]; i -= 1; j-=1; // We're at the x and y boundary!
118 119 120 121
	}
	
	// Get values for x1, x2, y1, y2
	x1 = g_x[i]; x2 = g_x[i+1];
122 123
	y1 = g_y[j]; y2 = g_y[j+1];
	
124
	// Do a simple (bilinear) interpolation
125
	// FIXME: This is wrong! May not be [i, j+1] for a shifted grid
126 127 128 129 130 131 132 133 134 135 136 137 138 139
	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;
	}
	
140
	// Return linear (grid) interpolation value
141
	return (1./((y2 - y1)*(x2 - x1))) * (
142 143 144 145 146 147
				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)
		);
}