1 files changed, 140 insertions, 0 deletions
diff --git a/astro/cosmo_calc/calc_distance.cc b/astro/cosmo_calc/calc_distance.cc
new file mode 100644
index 0000000..4c96dc8
--- /dev/null
+++ b/astro/cosmo_calc/calc_distance.cc
@@ -0,0 +1,140 @@
+// calc_distance.cc
+// calculate commonly used cosmology distances
+// Ref: https://en.wikipedia.org/wiki/Distance_measures_(cosmology)
+//
+// Junhua Gu
+//
+// ChangeLogs:
+// 2012/08/12, LIweitiaNux
+//   fix a bug in `calc_angular_distance()'
+//   add `calc_transcomv_distance()' and
+//     account `omega_k >0, <0' cases
+//   add `calc_ligtrav_distance()'
+//
+
+#include <iostream>
+#include <cmath>
+#include <cstdlib>
+#include <cstddef>
+#include <cassert>
+// integrator
+#include "adapt_trapezoid.h"
+
+// extern variables in `calc_distance.h'
+extern double c;
+extern double km;
+extern double s;
+extern double Mpc;
+// cosmology parameters
+extern double H0;       // units: [km/s/Mpc]
+extern double omega_m;
+extern double omega_l;
+extern double omega_k;
+// precision to determine float number equality
+extern double eps;
+
+using namespace std;
+
+// auxiliary functions
+// for `calc_comoving_distance'
+double f_comvdist(double z);
+// for `calc_ligtrav_distance'
+double f_ligdist(double z);
+double f_age(double z);
+
+/////////////////////////////////////////////////
+// main functions
+// //////////////////////////////////////////////
+// dimensionless Hubble parameter
+double E(double z)
+{
+    return sqrt(omega_m*(1+z)*(1+z)*(1+z) + omega_k*(1+z)*(1+z) + omega_l);
+}
+
+// Hubble distance
+double calc_hubble_distance(void)
+{
+    //
+    // cerr << "*** H0 = " << H0 << " ***" << endl;
+    //
+    double _H0_ = H0 * km/s/Mpc;
+    return c / _H0_;
+}
+
+// calc `comoving distance'
+double calc_comoving_distance(double z)
+{
+    double d_H = calc_hubble_distance();
+    return d_H * adapt_trapezoid(f_comvdist, 0.0, z, 1e-4);
+}
+
+// calc `transverse comoving distance'
+double calc_transcomv_distance(double z)
+{
+    if (fabs(omega_k) <= eps) {
+        // omega_k = 0
+        // flat, d_M(z) = d_C(z)
+        return calc_comoving_distance(z);
+    }
+    else if (omega_k > 0.0) {
+        // omega_k > 0
+        // d_M(z) = d_H/sqrt(omega_k) * sinh(sqrt(omega_k) * d_C(z) / d_H)
+        //                   ^             ^
+        double d_H = calc_hubble_distance();
+        double d_C = calc_comoving_distance(z);
+        return (d_H / sqrt(omega_k) * sinh(sqrt(omega_k) * d_C / d_H));
+    }
+    else {
+        // omega_k < 0
+        // d_M(z) = d_H/sqrt(-omega_k) * sin(sqrt(-omega_k) * d_C(z) / d_H)
+        //                   ^             ^      ^
+        double d_H = calc_hubble_distance();
+        double d_C = calc_comoving_distance(z);
+        return (d_H / sqrt(-omega_k) * sin(sqrt(-omega_k) * d_C / d_H));
+    }
+}
+
+// calc `angular diameter distance'
+// d_A(z) = d_M(z)/(1+z)
+double calc_angdia_distance(double z)
+{
+  return (calc_transcomv_distance(z) / (1+z));
+}
+
+// calc `luminoisity distance'
+// d_L(z) = (1+z)*d_M(z)
+double calc_luminosity_distance(double z)
+{
+  return (calc_transcomv_distance(z) * (1+z));
+}
+
+// calc `light-travel distance'
+// d_T(z) = d_H \int_0^z 1/((1+z)*E(z)) dz
+double calc_ligtrav_distance(double z)
+{
+    double d_H = calc_hubble_distance();
+    return d_H * adapt_trapezoid(f_ligdist, 0.0, z, 1e-4);
+}
+
+
+// auxiliary functions
+// for `calc_comoving_distance'
+double f_comvdist(double z)
+{
+    return 1.0/E(z);
+}
+
+// for `calc_ligtrav_distance'
+double f_ligdist(double z)
+{
+    return 1.0 / ((1+z)*E(z));
+}
+
+double f_age(double z)
+{
+    return f_comvdist(1.0/z) / (z*z);
+}
+
+
+// EOF
+/* vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=cpp: */