-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcal_a_vir.py
70 lines (50 loc) · 2.09 KB
/
cal_a_vir.py
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
import numpy as np
import scipy as sp
# ==================== PARAMETERS ====================
max_iter = int(1e6) # maximum iteration
# ==================== SCALE FACTOR AT THE VIRIALIZATION ====================
# --------------------- LCDM
def cal_a_vir_LCDM(a_ta, omega):
def integrand_1(y):
return np.sqrt(y / (a_ta**(-3) + omega*y**3))
def integrand_2(y):
return np.sqrt(y / (a_ta**(-3) + omega*y**3))
def equation(y_vir):
result_1, _ = sp.integrate.quad(integrand_1, 0, y_vir)
result_2, _ = sp.integrate.quad(integrand_2, 0, 1)
return 2*result_2 - result_1
y_vir = sp.optimize.ridder(equation, 1, 100, maxiter=max_iter)
a_vir = y_vir * a_ta
return a_vir
# --------------------- LsCDM (Pre-Turnaround)
def cal_pre_a_vir_LsCDM(a_ta, omega, y_dag):
def integrand_1(y):
return np.sqrt(y / (a_ta**(-3) + omega*y**3))
def integrand_2(y):
return np.sqrt(y / (a_ta**(-3) - omega*y**3))
def integrand_3(y):
return np.sqrt(y / (a_ta**(-3) + omega*y**3))
def equation(y_vir):
result_1, _ = sp.integrate.quad(integrand_1, 1, y_vir)
result_2, _ = sp.integrate.quad(integrand_2, 0, y_dag)
result_3, _ = sp.integrate.quad(integrand_3, y_dag, 1)
return result_1 - result_2 - result_3
y_vir = sp.optimize.ridder(equation, 1, 100, maxiter=max_iter)
a_vir = y_vir * a_ta
return a_vir
# --------------------- LsCDM (Post-Turnaround)
def cal_post_a_vir_LsCDM(a_ta, omega, y_dag):
def integrand_1(y):
return np.sqrt(y / (a_ta**(-3) - omega*y**3))
def integrand_2(y):
return np.sqrt(y / (a_ta**(-3) + omega*y**3))
def integrand_3(y):
return np.sqrt(y / (a_ta**(-3) - omega*y**3))
def equation(y_vir):
result_1, _ = sp.integrate.quad(integrand_1, 1, y_dag)
result_2, _ = sp.integrate.quad(integrand_2, y_dag, y_vir)
result_3, _ = sp.integrate.quad(integrand_3, 0, 1)
return result_1 + result_2 - result_3
y_vir = sp.optimize.ridder(equation, 1, 100, maxiter=max_iter)
a_vir = y_vir * a_ta
return a_vir