diff --git a/src/cosmology.c b/src/cosmology.c index 1f74f696a48c14ec3a56357dbb6b1466f505f244..7da8f7594c95fc2a2a1133e46d9b36704a40fa27 100644 --- a/src/cosmology.c +++ b/src/cosmology.c @@ -650,6 +650,47 @@ double cosmology_get_a_from_z(const struct cosmology *c, double redshift) { return 1. / (1. + redshift); } +/** + * @brief Compute scale factor from time since big bang. + * + * WARNING: This function is slow, therefore it should not be called more than + * once per time step. + * + * Find the scale factor from the time since big bang using the bisection method + * on $ t - f(a) = 0 $ where $ f(a) $ is #cosmology_get_time_since_big_bang + * + * By knowing that this function is a decreasing monotonic function, we can + * avoid a few checks. + * + * @param cosmo The current #cosmology. + * @param t since the big bang + * @return The scale factor. + */ +double cosmology_get_scale_factor(const struct cosmology *cosmo, double t) { + /* convergence limit and current error */ + const double limit = 1e-8; + double eps = 1.; + + /* limits in scale factor */ + double a_low = 1e-9; + double a_up = 1.; + + /* do the bisection method */ + while (eps > limit) { + double a_tmp = 0.5 * (a_low + a_up); + double f_tmp = t - cosmology_get_time_since_big_bang(cosmo, a_tmp); + + if (f_tmp > 0) + a_low = a_tmp; + else + a_up = a_tmp; + + eps = a_up - a_low; + } + + return 0.5 * (a_low + a_up); +} + /** * @brief Prints the #cosmology model to stdout. */ diff --git a/src/cosmology.h b/src/cosmology.h index 3b998c4b287e4bdcc2ce1c4d72a59100b62b5b0f..6f6f1462038ce3c018e670ccada931f49aedc285 100644 --- a/src/cosmology.h +++ b/src/cosmology.h @@ -186,6 +186,8 @@ double cosmology_get_delta_time(const struct cosmology *c, double a1, double cosmology_get_a_from_z(const struct cosmology *c, double redshift); +double cosmology_get_scale_factor(const struct cosmology *cosmo, double t); + double cosmology_get_time_since_big_bang(const struct cosmology *c, double a); void cosmology_init(struct swift_params *params, const struct unit_system *us, const struct phys_const *phys_const, struct cosmology *c);