forked from PedroCarrilho/evogrowthpy
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathEvo.cpp
128 lines (99 loc) · 3.09 KB
/
Evo.cpp
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
#if HAVE_CONFIG_H
# include <config.h>
#endif
#include <cfloat>
#include <cmath>
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_min.h>
#include "Evo.h"
#include <math.h> /* pow */
#include <iostream>
#include <stdlib.h>
#include <functional>
/*cosmological parameters for evolution factors (applied to analytic solutions) */
struct param_type2 {
double omega00;
double omega_cb;
double hubble;
double par1;
double par2;
double par3;
double par4;
};
int jac (double a, const double G[], double *dfdy, double dfdt[], void *params)
{
return GSL_SUCCESS;
}
// fixed omega
// H/H0
static double HAde(double a, double omega0, double w0, double wa){
double A = -3.*(1.+w0+wa);
double omegaf = pow(a,A)*exp(3.*(-1.+a)*wa);
double omegaL= (1.-omega0)*omegaf;
return sqrt(omega0/pow(a,3)+omegaL);
}
// aH dH/da / H0^2
static double HA1de(double a, double omega0, double w0, double wa){
double A = -3.*(1.+w0+wa);
double omegaf = pow(a,A)*exp(3.*(-1.+a)*wa);
double omegaL= (1.-omega0)*omegaf;
return -3*omega0/(2.*pow(a,3)) + (A+3.*a*wa)*omegaL/2.;
}
//HA2de = -dH/dt/H^2 = -a dH/da / H
static double HA2de(double a, double omega0, double w0, double wa){
return -HA1de(a,omega0,w0,wa)/pow(HAde(a,omega0,w0,wa),2);
}
//3/(2H^2) * Omega_m_0 *a^3 = Omega_m(a)
static double HA2de2(double a, double omega0, double w0, double wa){
return 3.*omega0/(2.*pow(HAde(a,omega0,w0,wa),2)*pow(a,3));
}
// linear growth factor for wCDM CPL with interaction and nDGP
int evo_eqs_ide(double a, const double G[], double F[], void *params)
{
param_type2 p = *(param_type2 *)(params);
double omega0 = p.omega00;
double omega_cb = p.omega_cb;
double hubble=p.hubble;
double w0 = p.par1;
double wa = p.par2;
double xi = p.par3;
double Omega_rc = p.par4;
double A, omegaf, omegaL, fric,hade,hade1,hade2,f_cb,muphi;
A = -3.*(1.+w0+wa);
omegaf = pow(a,A)*exp(3.*(-1.+a)*wa);
omegaL= (1.-omega0)*omegaf;
hade = HAde(a,omega0,w0,wa);
fric = (1.+w0+(1.-a)*wa)*omegaL*xi*hubble/hade*0.0974655;
hade1 = HA2de(a,omega0,w0,wa);
hade2 = HA2de2(a,omega0,w0,wa);
if(Omega_rc == 0){
muphi = 1.0;
}
else{
muphi = 1.0 + 1.0 / (3.0 * (1.0 + hade / sqrt(Omega_rc) * (1.0 + HA1de(a,omega0,w0,wa)/(3.0*hade*hade))));
}
f_cb = omega_cb/omega0;
F[0] = -G[1]/a;
F[1] =1./a*(-(2.+fric-hade1)*G[1]-f_cb*hade2*muphi*G[0]);
return GSL_SUCCESS;
}
void get_growth(double A, double omega0, double omegacb, double h,double w0, double wa, double xi, double Om_rc, double * res, double accuracy)
{
struct param_type2 params_my = {omega0,omegacb,h,w0,wa,xi,Om_rc};
//Initial a
double a = 0.0001;
double G[2] = { a,-a};
gsl_odeiv2_system sys = {evo_eqs_ide, jac, 2, ¶ms_my};
// this gives the system, the method, the initial interval step, the absolute error and the relative error.
gsl_odeiv2_driver * d =
gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd,
1e-3, accuracy, 0.);
int status1 = gsl_odeiv2_driver_apply (d, &a, A, G);
res[0] = G[0];
res[1] = -G[1]/G[0];
gsl_odeiv2_driver_free(d);
}