diff options
author | Weitian LI <liweitianux@gmail.com> | 2014-06-18 22:20:59 +0800 |
---|---|---|
committer | Weitian LI <liweitianux@gmail.com> | 2014-06-18 22:20:59 +0800 |
commit | e3923265d0d6949407a83726e9a9bd5d97079221 (patch) | |
tree | 9afd8520595f4cf80815b9bccfc3dcf2879ebe47 /tools/cosmo_calc/calc_distance.cc | |
download | chandra-acis-analysis-e3923265d0d6949407a83726e9a9bd5d97079221.tar.bz2 |
Initial commit
Added files:
* mass_profile: developed by Junhua GU, modified by Weitian LI
* opt_utilities: developed by Junhua GU
* tools/cosmo_calc: originated from 'calc_distance', modified
* scripts: scripts used to process Chandra ACIS data
* files: useful files used in processing
* HOWTO_chandra_acis_process.txt
* README.md
Diffstat (limited to 'tools/cosmo_calc/calc_distance.cc')
-rw-r--r-- | tools/cosmo_calc/calc_distance.cc | 140 |
1 files changed, 140 insertions, 0 deletions
diff --git a/tools/cosmo_calc/calc_distance.cc b/tools/cosmo_calc/calc_distance.cc new file mode 100644 index 0000000..4c96dc8 --- /dev/null +++ b/tools/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: */ |