diff --git a/src/cosmology.c b/src/cosmology.c index 6983c033f65e5ae2147685ceed35bf8de672a743..0eb60965065956c13dc3c3040691e73a0ef7e8e2 100644 --- a/src/cosmology.c +++ b/src/cosmology.c @@ -82,14 +82,17 @@ static INLINE double interp_table(const double *restrict y_table, #ifdef SWIFT_DEBUG_CHECKS if (ii <= 1) error("Wrong ii - too small!"); if (ii >= cosmology_table_length - 1) error("Wrong ii - too large!"); + if (x < x_table[ii - 2] || x > x_table[ii + 1]) + error("Extrapolating not interpolating!"); #endif /* Initialise the interpolation range with 4 data points * The point of interest is in the range [ii - 1, ii] */ - gsl_interp_init(poly, &x_table[ii - 2], &y_table[ii - 2], 4); + // gsl_interp_init(poly, &x_table[ii - 2], &y_table[ii - 2], 4); /* Interpolate! */ - return gsl_interp_eval(poly, &x_table[ii - 2], &y_table[ii - 2], x, acc); + // return gsl_interp_eval(poly, &x_table[ii - 2], &y_table[ii - 2], x, acc); + return y_table[ii - 1] + (y_table[ii] - y_table[ii - 1]) * (xx - ii); } /** @@ -351,7 +354,7 @@ void cosmology_init_tables(struct cosmology *c) { /* We want to extend the tables on both sides so that we never have to worry about edges when interpolating */ - const double fac = 1. + 10. / ((double)cosmology_table_length); + const double fac = 1. + 16. / ((double)cosmology_table_length); const double a_begin = c->a_begin / fac; const double a_end = c->a_end * fac; c->log_a_table_begin = log(a_begin); @@ -641,6 +644,8 @@ double cosmology_get_drift_factor(const struct cosmology *c, if (ti_end < ti_start) error("ti_end must be >= ti_start"); #endif + // if (ti_start == ti_end) return 0.; + const double log_a_start = c->log_a_begin + ti_start * c->time_base; const double log_a_end = c->log_a_begin + ti_end * c->time_base; @@ -671,6 +676,8 @@ double cosmology_get_grav_kick_factor(const struct cosmology *c, if (ti_end < ti_start) error("ti_end must be >= ti_start"); #endif + // if (ti_start == ti_end) return 0.; + const double log_a_start = c->log_a_begin + ti_start * c->time_base; const double log_a_end = c->log_a_begin + ti_end * c->time_base; @@ -702,6 +709,8 @@ double cosmology_get_hydro_kick_factor(const struct cosmology *c, if (ti_end < ti_start) error("ti_end must be >= ti_start"); #endif + // if (ti_start == ti_end) return 0.; + const double log_a_start = c->log_a_begin + ti_start * c->time_base; const double log_a_end = c->log_a_begin + ti_end * c->time_base; @@ -733,6 +742,8 @@ double cosmology_get_corr_kick_factor(const struct cosmology *c, if (ti_end < ti_start) error("ti_end must be >= ti_start"); #endif + // if (ti_start == ti_end) return 0.; + const double log_a_start = c->log_a_begin + ti_start * c->time_base; const double log_a_end = c->log_a_begin + ti_end * c->time_base; @@ -764,6 +775,8 @@ double cosmology_get_therm_kick_factor(const struct cosmology *c, if (ti_end < ti_start) error("ti_end must be >= ti_start"); #endif + // if (ti_start == ti_end) return 0.; + const double log_a_start = c->log_a_begin + ti_start * c->time_base; const double log_a_end = c->log_a_begin + ti_end * c->time_base; @@ -814,6 +827,8 @@ double cosmology_get_delta_time(const struct cosmology *c, if (ti_end < ti_start) error("ti_end must be >= ti_start"); #endif + // if (ti_start == ti_end) return 0.; + const double log_a_start = c->log_a_begin + ti_start * c->time_base; const double log_a_end = c->log_a_begin + ti_end * c->time_base; diff --git a/tests/testIntegration.c b/tests/testIntegration.c index 27f50ba822735c2f0b09c19ae983b43fa203b913..8de857849dd60ae3fff2f46e8f48146fd3902dae 100644 --- a/tests/testIntegration.c +++ b/tests/testIntegration.c @@ -160,7 +160,7 @@ int main(int argc, char *argv[]) { message("Sum error: %14e", sum_err); message("Mean error: %14e", sum_err / num_steps); - if (max_err > 1e-8 || min_err < -1e-8) + if (max_err > 1e-4 || min_err < -1e-4) error("Error too large to be acceptable"); }