diff options
| -rw-r--r-- | .gitignore | 3 | ||||
| -rw-r--r-- | astro/cosmo_calc/.gitignore | 2 | ||||
| -rw-r--r-- | astro/cosmo_calc/Makefile | 31 | ||||
| -rw-r--r-- | astro/cosmo_calc/README.txt | 20 | ||||
| -rw-r--r-- | astro/cosmo_calc/VERSION.txt | 10 | ||||
| -rw-r--r-- | astro/cosmo_calc/adapt_trapezoid.h | 124 | ||||
| -rw-r--r-- | astro/cosmo_calc/calc_distance.cc | 140 | ||||
| -rw-r--r-- | astro/cosmo_calc/calc_distance.h | 59 | ||||
| -rw-r--r-- | astro/cosmo_calc/ddivid.cc | 56 | ||||
| -rw-r--r-- | astro/cosmo_calc/ddivid.h | 8 | ||||
| -rw-r--r-- | astro/cosmo_calc/main.cc | 140 |
11 files changed, 593 insertions, 0 deletions
@@ -11,3 +11,6 @@ __init__.py *.pyc **/__pycache__ +# C/C++ +*.o + diff --git a/astro/cosmo_calc/.gitignore b/astro/cosmo_calc/.gitignore new file mode 100644 index 0000000..009f466 --- /dev/null +++ b/astro/cosmo_calc/.gitignore @@ -0,0 +1,2 @@ +# ignore the executable +cosmo_calc diff --git a/astro/cosmo_calc/Makefile b/astro/cosmo_calc/Makefile new file mode 100644 index 0000000..abd12bb --- /dev/null +++ b/astro/cosmo_calc/Makefile @@ -0,0 +1,31 @@ +## +## Makefile for `calc_distance' +## +## Junhua Gu +## last modified: August 12, 2012 +## + +CPP= g++ +CPPFLAGS= -Wall -g + +TARGET= cosmo_calc + +all: $(TARGET) + +cosmo_calc: calc_distance.o ddivid.o main.o + $(CPP) $(CPPFLAGS) $^ -o $@ + +calc_distance.o: calc_distance.cc calc_distance.h + $(CPP) $(CPPFLAGS) -c $< + +ddivid.o: ddivid.cc ddivid.h + $(CPP) $(CPPFLAGS) -c $< + +main.o: main.cc ddivid.h calc_distance.h + $(CPP) $(CPPFLAGS) -c $< + +clean: + rm -f *.o + rm -f $(TARGET) + + diff --git a/astro/cosmo_calc/README.txt b/astro/cosmo_calc/README.txt new file mode 100644 index 0000000..b18e120 --- /dev/null +++ b/astro/cosmo_calc/README.txt @@ -0,0 +1,20 @@ +// +// simple `cosmology calculator' +// can calc commonly used cosmology distances: +// hubble distance +// comoving distance +// transverse comoving distance +// angular diameter distance +// luminoisity distance +// light-travel distance +// in addition, this calculator also calc some other +// useful infomation related with Chandra +// +// Junhua Gu +// LIweitiaNux <liweitianux@gmail.com> +// +// Version 2.0 +// +// Last Modified: +// August 12, 2012 +// diff --git a/astro/cosmo_calc/VERSION.txt b/astro/cosmo_calc/VERSION.txt new file mode 100644 index 0000000..f99d493 --- /dev/null +++ b/astro/cosmo_calc/VERSION.txt @@ -0,0 +1,10 @@ +Version 2.2 +2013/02/09 + +ChangeLogs: +v2.2: + 2013/02/09, LIweitiaNux + add 'hubble_parameter E(z)' + modify output format + + diff --git a/astro/cosmo_calc/adapt_trapezoid.h b/astro/cosmo_calc/adapt_trapezoid.h new file mode 100644 index 0000000..35d22b4 --- /dev/null +++ b/astro/cosmo_calc/adapt_trapezoid.h @@ -0,0 +1,124 @@ +#ifndef ADAPT_TRAPEZOID_H +#define ADAPT_TRAPEZOID_H +#include <list> +#include <utility> +#include <cassert> +//#include <iostream> + +//using namespace std; + +template<typename T1,typename T2,typename T3> +class triple +{ + public: + T1 first; + T2 second; + T3 third; + triple(T1 x1,T2 x2,T3 x3) + :first(x1),second(x2),third(x3) + { + } +}; + + +template<typename T1,typename T2,typename T3> +triple<T1,T2,T3> make_triple(T1 x1,T2 x2,T3 x3) +{ + return triple<T1,T2,T3>(x1,x2,x3); +} + +template <typename T> +T trapezoid(T (*fun)(T),T x1,T x2,T err_limit) +{ + int n=256; + T result; + const int max_division=24; + T old_value=0; + for(int i=1;i<max_division;i++) + { + result=0.; + n*=2; + T step=(x2-x1)/n; + for(int j=0;j<n;j++) + { + result+=(fun(x1+(j+1)*step)+fun(x1+j*step))*step/T(2.); + } + old_value-=result; + old_value=old_value<0?-old_value:old_value; + if(old_value<err_limit) + { + return result; + } + old_value=result; + } +} + + +template <typename T> +T adapt_trapezoid(T (*fun)(T),T x1,T x2,T err_limit) +{ + // const err_limit=.001; + typedef triple<T,T,bool> interval; + /*T for interval type, + bool for state trur for still to be updated, + false for do not need to be updated + */ + std::list<interval> interval_list; + T current_sum=((fun(x1)+fun(x2))/2.*(x2-x1)); + interval_list.push_back(make_triple(x1,current_sum,true)); + interval_list.push_back(make_triple(x2,(T)0.,true)); + bool int_state=1; + int n_intervals=1; + while(int_state) + { + //std::cout<<n_intervals<<std::endl; + int_state=0; + typename std::list<interval>::iterator i1=interval_list.begin(); + typename std::list<interval>::iterator i2=interval_list.begin(); + i2++; + for(;i2!=interval_list.end();i1++,i2=i1,i2++) + { + //cout<<i1->first<<"\t"<<i2->first<<endl; + //assert(i2->first>i1->first); + if(i1->third) + { + interval new_interval((i1->first+i2->first)/2,0,true); + + T sum1,sum2; + sum1=(fun(new_interval.first)+fun(i1->first))/2*(new_interval.first-i1->first); + sum2=(fun(new_interval.first)+fun(i2->first))/2*(i2->first-new_interval.first); + new_interval.second=sum2; + T err; + err=i1->second-sum1-sum2; + err=err<0?-err:err; + + if(err>err_limit/n_intervals) + { + i1->second=sum1; + interval_list.insert(i2,new_interval); + n_intervals++; + if(n_intervals>10e6) + { + + break; + } + } + else + { + i1->third=false; + } + int_state=1; + } + } + + } + T result=0; + for(typename std::list<interval>::iterator i=interval_list.begin();i!=interval_list.end();i++) + { + result+=i->second; + } + return result; +} + +#endif +//end of the file 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: */ diff --git a/astro/cosmo_calc/calc_distance.h b/astro/cosmo_calc/calc_distance.h new file mode 100644 index 0000000..a1f3352 --- /dev/null +++ b/astro/cosmo_calc/calc_distance.h @@ -0,0 +1,59 @@ +// calc_distance.h +// 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 +// + +#ifndef CALC_DISTANCE_H +#define CALC_DISTANCE_H + +double cm=1.0; +double s=1.0; +double km=1000*100*cm; +double Mpc=3.08568e+24*cm; +double kpc=3.08568e+21*cm; +double yr=365.*24.*3600.*s; +double Gyr=1e9*yr; +double arcsec2arc_ratio=1./60/60/180*3.1415926; +double c=299792458.*100.*cm; + +// cosmology parameters +double H0=71.0; // units: [km/s/Mpc] +double omega_m=0.27; +double omega_l=0.73; +double omega_k=1.0-omega_m-omega_l; + +// precision to determine float number equality +double eps=1.0e-7; + +// open functions +// +// dimensionless Hubble parameter +extern double E(double z); +// Hubble distance +extern double calc_hubble_distance(void); +// calc `comoving distance' +extern double calc_comoving_distance(double z); +// calc `transverse comoving distance' +extern double calc_transcomv_distance(double z); +// calc `angular diameter distance' +// d_A(z) = d_M(z)/(1+z) +extern double calc_angdia_distance(double z); +// calc `luminoisity distance' +// d_L(z) = (1+z)*d_M(z) +extern double calc_luminosity_distance(double z); +// calc `light-travel distance' +// d_T(z) = d_H \int_0^z 1/((1+z)*E(z)) dz +extern double calc_ligtrav_distance(double z); + +#endif /* CALC_DISTANCE_H */ + +//EOF +/* vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=cpp: */ diff --git a/astro/cosmo_calc/ddivid.cc b/astro/cosmo_calc/ddivid.cc new file mode 100644 index 0000000..aa72471 --- /dev/null +++ b/astro/cosmo_calc/ddivid.cc @@ -0,0 +1,56 @@ +#include <cassert> + +double ddivid(double (*foo)(double),double x1,double x2,double err) +{ + assert(x2>x1); + assert(foo(x2)*foo(x1)<0); + while(x2-x1>err) + { + double x=(x2+x1)/2.; + double y=foo(x); + double y1=foo(x1); + double y2=foo(x2); + if(y1*y<0) + { + x2=x; + } + else if(y2*y<0) + { + x1=x; + } + else + { + assert(0); + } + + } + return (x1+x2)/2.; +} + + +double ddivid(double (*foo)(double),double z,double x1,double x2,double err) +{ + assert(x2>x1); + assert((foo(x2)-z)*(foo(x1)-z)<0); + while(x2-x1>err) + { + double x=(x2+x1)/2.; + double y=foo(x)-z; + double y1=foo(x1)-z; + double y2=foo(x2)-z; + if(y1*y<0) + { + x2=x; + } + else if(y2*y<0) + { + x1=x; + } + else + { + assert(0); + } + + } + return (x1+x2)/2.; +} diff --git a/astro/cosmo_calc/ddivid.h b/astro/cosmo_calc/ddivid.h new file mode 100644 index 0000000..a957de5 --- /dev/null +++ b/astro/cosmo_calc/ddivid.h @@ -0,0 +1,8 @@ +#ifndef DDIVID +#define DDIVID + +extern double ddivid(double (*foo)(double),double x1,double x2,double err); +extern double ddivid(double (*foo)(double),double z,double x1,double x2,double err); + +#endif +//EOF diff --git a/astro/cosmo_calc/main.cc b/astro/cosmo_calc/main.cc new file mode 100644 index 0000000..bbd73ba --- /dev/null +++ b/astro/cosmo_calc/main.cc @@ -0,0 +1,140 @@ +// +// simple `cosmology calculator' +// can calc commonly used cosmology distances: +// hubble distance +// comoving distance +// transverse comoving distance +// angular diameter distance +// luminoisity distance +// light-travel distance +// in addition, this calculator also calc some other +// useful infomation related with Chandra +// +// Junhua Gu +// +// Modified by: LIweitiaNux +// +// ChangeLogs: +// v2.1, 2012/08/12, LIweitiaNux +// improve cmdline parameters +// v2.2, 2013/02/09, LIweitiaNux +// add 'hubble_parameter E(z)' +// modify output format +// + + +#include "calc_distance.h" +#include "ddivid.h" +#include <cstdio> +#include <cstdlib> +#include <cmath> +#include <iostream> + +using namespace std; + +// extern variables in `calc_distance.h' +extern double cm; +extern double s; +extern double km; +extern double Mpc; +extern double kpc; +extern double yr; +extern double Gyr; +extern double H0; // units: [km/s/Mpc] +extern double c; +extern double omega_m; +extern double omega_l; +extern double omega_k; +extern double arcsec2arc_ratio; + +// README, ABOUT +static char DESC[] = "simple cosmology calculator"; +static char VER[] = "v2.2, 2013-02-09"; + +// setting parameters +static double PI = 4*atan(1.0); +// chandra related +static double arcsec_per_pixel = 0.492; +// above cosmology paramters can also be changed +// through cmdline paramters + +// functions +void usage(const char *name); + +////////////////////////////////////////////// +// main part +////////////////////////////////////////////// +int main(int argc,char* argv[]) +{ + double z; // given redshift + + // cmdline parameters + if (argc == 2) { + // get input redshift + z = atof(argv[1]); + } + else if (argc == 3) { + z = atof(argv[1]); + // use specified `H0' + H0 = atof(argv[2]); + } + else if (argc == 4) { + z = atof(argv[1]); + H0 = atof(argv[2]); + // get specified `Omega_M' + omega_m = atof(argv[3]); + omega_l = 1.0-omega_m; + } + else { + usage(argv[0]); + exit(-1); + } + + // calc `Hubble parameter E(z)' + double E_z = E(z); + // calc `comoving distance' + double d_C = calc_comoving_distance(z); + // calc `angular diameter distance' + double d_A = calc_angdia_distance(z); + // calc `luminoisity distance' + double d_L = calc_luminosity_distance(z); + + // output results + // parameters + printf("Parameters:\n"); + printf(" z= %lf, H0= %lf, Omega_M= %lf, Omega_L= %lf\n", + z, H0, omega_m, omega_l); + printf("Distances:\n"); + printf(" Comoving_distance: D_C(%lf)= %lg [cm], %lf [Mpc]\n", + z, d_C, d_C/Mpc); + printf(" Angular_diameter_distance: D_A(%lf)= %lg [cm], %lf [Mpc]\n", + z, d_A, d_A/Mpc); + printf(" Luminoisity_distance: D_L(%lf)= %lg [cm], %lf [Mpc]\n", + z, d_L, d_L/Mpc); + printf("Chandra_related:\n"); + printf(" kpc/pixel (D_A): %lf\n", + (d_A / kpc * arcsec_per_pixel* arcsec2arc_ratio)); + printf(" cm/pixel (D_A): %lg\n", + (d_A * arcsec_per_pixel* arcsec2arc_ratio)); + printf("Other_data:\n"); + printf(" Hubble_parameter: E(%lf)= %lf\n", z, E_z); + printf(" kpc/arcsec (D_A): %lf\n", (d_A / kpc * arcsec2arc_ratio)); + printf(" norm (cooling_function): %lg\n", + (1e-14 / (4.0 * PI * pow(d_A*(1+z), 2)))); + + //cout<<ddivid(calc_distance,d,0,1,.0001)<<endl; + + return 0; +} + +// other auxiliary functions +void usage(const char *name) +{ + cerr << "Usage: " << endl; + cerr << " " << name << " z [H0] [Omega_M]" << endl; + // description + cout << endl << "About:" << endl; + cout << DESC << endl << VER << endl; +} + +/* vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=cpp: */ |