1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
|
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstddef>
#include <cassert>
#include "adapt_trapezoid.h"
//calc_distance
//usage:
//calc_distance z
using namespace std;
static double cm=1;
static double s=1;
static double km=1000*100;
static double Mpc=3.08568e+24*cm;
static double kpc=3.08568e+21*cm;
static double yr=365.*24.*3600.;
static double Gyr=1e9*yr;
static double H=71.*km/s/Mpc;
static const double c=299792458.*100.*cm;
//const double c=3e8*100*cm;
static const double omega_m=0.27;
static const double omega_l=0.73;
static const double arcsec2arc_ratio=1./60/60/180*3.1415926;
double E(double z)
{
double omega_k=1-omega_m-omega_l;
return sqrt(omega_m*(1+z)*(1+z)*(1+z)+omega_k*(1+z)*(1+z)+omega_l);
}
double f_dist(double z)
{
return 1/E(z);
}
double f_age(double z)
{
return f_dist(1/z)/(z*z);
}
double calc_angular_distance(double z)
{
//return c/H*integer(f_dist,0,z)/(1+z);
//return c/H*adapt_trapezoid(f_dist,0.,z,1e-4)/(1+z);
return adapt_trapezoid(f_dist,0.,z,1e-4)/(1+z);
}
double calc_luminosity_distance(double z)
{
//return c/H*integer(f_dist,0,z)/(1+z);
return c/H*adapt_trapezoid(f_dist,0.,z,1e-4)*(1+z);
}
//EOF
|
