From cb0d7c26642467b44f4a654141fc5ca1f3e78d97 Mon Sep 17 00:00:00 2001
From: lhausamm <loic_hausammann@hotmail.com>
Date: Mon, 4 Jun 2018 10:19:28 +0200
Subject: [PATCH] Add function giving scale factor from time
---
src/cosmology.c | 41 +++++++++++++++++++++++++++++++++++++++++
src/cosmology.h | 2 ++
2 files changed, 43 insertions(+)
diff --git a/src/cosmology.c b/src/cosmology.c
index 1f74f696a4..7da8f7594c 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 3b998c4b28..6f6f146203 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);
--
GitLab