diff options
Diffstat (limited to 'astro/cosmo_calc')
| -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 | 
10 files changed, 0 insertions, 590 deletions
| diff --git a/astro/cosmo_calc/.gitignore b/astro/cosmo_calc/.gitignore deleted file mode 100644 index 009f466..0000000 --- a/astro/cosmo_calc/.gitignore +++ /dev/null @@ -1,2 +0,0 @@ -# ignore the executable -cosmo_calc diff --git a/astro/cosmo_calc/Makefile b/astro/cosmo_calc/Makefile deleted file mode 100644 index abd12bb..0000000 --- a/astro/cosmo_calc/Makefile +++ /dev/null @@ -1,31 +0,0 @@ -## -## 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 deleted file mode 100644 index b18e120..0000000 --- a/astro/cosmo_calc/README.txt +++ /dev/null @@ -1,20 +0,0 @@ -// -// 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 deleted file mode 100644 index f99d493..0000000 --- a/astro/cosmo_calc/VERSION.txt +++ /dev/null @@ -1,10 +0,0 @@ -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 deleted file mode 100644 index 35d22b4..0000000 --- a/astro/cosmo_calc/adapt_trapezoid.h +++ /dev/null @@ -1,124 +0,0 @@ -#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 deleted file mode 100644 index 4c96dc8..0000000 --- a/astro/cosmo_calc/calc_distance.cc +++ /dev/null @@ -1,140 +0,0 @@ -// 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 deleted file mode 100644 index a1f3352..0000000 --- a/astro/cosmo_calc/calc_distance.h +++ /dev/null @@ -1,59 +0,0 @@ -// 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 deleted file mode 100644 index aa72471..0000000 --- a/astro/cosmo_calc/ddivid.cc +++ /dev/null @@ -1,56 +0,0 @@ -#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 deleted file mode 100644 index a957de5..0000000 --- a/astro/cosmo_calc/ddivid.h +++ /dev/null @@ -1,8 +0,0 @@ -#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 deleted file mode 100644 index bbd73ba..0000000 --- a/astro/cosmo_calc/main.cc +++ /dev/null @@ -1,140 +0,0 @@ -// -// 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: */ | 
