Add astro/cosmo_calc from "chandra-acis-analysis"

11 files changed, 593 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
index c3edc65..91fcee7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -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: */