11# Define w0waCDMCosmology type
2- @kwdef mutable struct w0waCDMCosmology <: AbstractCosmology
3- ln10Aₛ:: Number = 3.0
4- nₛ:: Number = 0.96
5- h:: Number = 0.67
6- ωb:: Number = 0.022
7- ωc:: Number = 0.12
8- ωk:: Number = 0.0
9- mν:: Number = 0.0
10- w0:: Number = - 1.0
11- wa:: Number = 0.0
2+ @kwdef mutable struct w0waCDMCosmology{T1, T2, T3, T4, T5, T6, T7, T8, T9} <: AbstractCosmology
3+ ln10Aₛ:: T1 = 3.0
4+ nₛ:: T2 = 0.96
5+ h:: T3 = 0.67
6+ ωb:: T4 = 0.022
7+ ωc:: T5 = 0.12
8+ ωk:: T6 = 0.0
9+ mν:: T7 = 0.0
10+ w0:: T8 = - 1.0
11+ wa:: T9 = 0.0
1212end
1313
14+ _call_interpolant (interp:: Ref , y:: T ) where {T} = interp[](y):: T
15+
1416function _F (y)
1517 f (x, y) = x^ 2 * √ (x^ 2 + y^ 2 ) / (1 + exp (x))
1618 domain = (zero (eltype (Inf )), Inf ) # (lb, ub)
3335
3436function _ΩνE2 (a, Ωγ0, mν; kB= 8.617342e-5 , Tν= 0.71611 * 2.7255 , Neff= 3.044 )
3537 Γν = (4 / 11 )^ (1 / 3 ) * (Neff / 3 )^ (1 / 4 )
36- return 15 / π^ 4 * Γν^ 4 * Ωγ0 / a^ 4 * F_interpolant[](_get_y (mν, a))
38+ y = _get_y (mν, a)
39+ val = _call_interpolant (F_interpolant, y)
40+ return 15 / π^ 4 * Γν^ 4 * Ωγ0 / a^ 4 * val
3741end
3842
3943function _ΩνE2 (a, Ωγ0, mν:: AbstractVector ; kB= 8.617342e-5 , Tν= 0.71611 * 2.7255 , Neff= 3.044 )
4044 Γν = (4 / 11 )^ (1 / 3 ) * (Neff / 3 )^ (1 / 4 )
41- sum_interpolant = sum (mymν -> F_interpolant[](_get_y (mymν, a)), mν)
45+ sum_interpolant = sum (mymν -> begin
46+ y = _get_y (mymν, a)
47+ _call_interpolant (F_interpolant, y)
48+ end , mν)
4249 return 15 / π^ 4 * Γν^ 4 * Ωγ0 / a^ 4 * sum_interpolant
4350end
4451
4552function _dΩνE2da (a, Ωγ0, mν; kB= 8.617342e-5 , Tν= 0.71611 * 2.7255 , Neff= 3.044 )
4653 Γν = (4 / 11 )^ (1 / 3 ) * (Neff / 3 )^ (1 / 4 )
47- return 15 / π^ 4 * Γν^ 4 * Ωγ0 * (- 4 * F_interpolant[](_get_y (mν, a)) / a^ 5 +
48- dFdy_interpolant[](_get_y (mν, a)) / a^ 4 * (mν / kB / Tν))
54+ y = _get_y (mν, a)
55+ val_F = _call_interpolant (F_interpolant, y)
56+ val_dFdy = _call_interpolant (dFdy_interpolant, y)
57+ return 15 / π^ 4 * Γν^ 4 * Ωγ0 * (- 4 * val_F / a^ 5 +
58+ val_dFdy / a^ 4 * (mν / kB / Tν))
4959end
5060
5161function _dΩνE2da (a, Ωγ0, mν:: AbstractVector ; kB= 8.617342e-5 , Tν= 0.71611 * 2.7255 , Neff= 3.044 )
5262 Γν = (4 / 11 )^ (1 / 3 ) * (Neff / 3 )^ (1 / 4 )
53- sum_interpolant = 0.0
54- for mymν in mν
55- sum_interpolant += - 4 * F_interpolant[](_get_y (mymν, a)) / a^ 5 +
56- dFdy_interpolant[](_get_y (mymν, a)) / a^ 4 * (mymν / kB / Tν)
57- end
63+ sum_interpolant = sum (mymν -> begin
64+ y = _get_y (mymν, a)
65+ val_F = _call_interpolant (F_interpolant, y)
66+ val_dFdy = _call_interpolant (dFdy_interpolant, y)
67+ - 4 * val_F / a^ 5 + val_dFdy / a^ 4 * (mymν / kB / Tν)
68+ end , mν)
5869 return 15 / π^ 4 * Γν^ 4 * Ωγ0 * sum_interpolant
5970end
6071
@@ -74,11 +85,19 @@ function _dρDEda(a, w0, wa)
7485 return 3 * (- (1 + w0 + wa) / a + wa) * _ρDE_a (a, w0, wa)
7586end
7687
88+ function _E_a_scalar (a:: Number , Ωcb0, Ωγ0, Ων0, ΩΛ0, Ωk0, h, mν, w0, wa)
89+ return sqrt (Ωγ0 * a^- 4 + Ωcb0 * a^- 3 + Ωk0 * a^- 2 + ΩΛ0 * _ρDE_a (a, w0, wa) + _ΩνE2 (a, Ωγ0, mν))
90+ end
91+
7792function E_a (a, Ωcb0, h; mν= 0.0 , w0= - 1.0 , wa= 0.0 , Ωk0= 0.0 )
7893 Ωγ0 = 2.469e-5 / h^ 2
79- Ων0 = _ΩνE2 (1.0 , Ωγ0, mν)
94+ Ων0 = _ΩνE2 (1.0 , Ωγ0, mν):: promote_type (Float64, typeof (Ωγ0), typeof (mν) <: AbstractVector ? eltype (mν) : typeof (mν))
8095 ΩΛ0 = 1.0 - (Ωγ0 + Ωcb0 + Ων0 + Ωk0)
81- return @. sqrt (Ωγ0 * a^- 4 + Ωcb0 * a^- 3 + Ωk0 * a^- 2 + ΩΛ0 * _ρDE_a (a, w0, wa) + _ΩνE2 (a, Ωγ0, mν))
96+ if a isa AbstractArray
97+ return _E_a_scalar .(a, Ωcb0, Ωγ0, Ων0, ΩΛ0, Ωk0, h, Ref (mν), w0, wa)
98+ else
99+ return _E_a_scalar (a, Ωcb0, Ωγ0, Ων0, ΩΛ0, Ωk0, h, mν, w0, wa)
100+ end
82101end
83102
84103function E_a (a, cosmo:: w0waCDMCosmology )
@@ -98,16 +117,24 @@ function E_z(z, cosmo::w0waCDMCosmology)
98117 return E_z (z, Ωcb0, cosmo. h; mν= cosmo. mν, w0= cosmo. w0, wa= cosmo. wa, Ωk0= Ωk0)
99118end
100119
120+ function _dlogEdloga_scalar (a:: Number , Ωcb0, Ωγ0, Ων0, ΩΛ0, Ωk0, h, mν, w0, wa)
121+ return a * 0.5 / (_E_a_scalar (a, Ωcb0, Ωγ0, Ων0, ΩΛ0, Ωk0, h, mν, w0, wa)^ 2 ) *
122+ (- 3 (Ωcb0)a^- 4 - 4 Ωγ0 * a^- 5 - 2 Ωk0 * a^- 3 + ΩΛ0 * _dρDEda (a, w0, wa) + _dΩνE2da (a, Ωγ0, mν))
123+ end
124+
101125function _dlogEdloga (a, Ωcb0, h; mν= 0.0 , w0= - 1.0 , wa= 0.0 , Ωk0= 0.0 )
102126 Ωγ0 = 2.469e-5 / h^ 2
103- Ων0 = _ΩνE2 (1.0 , Ωγ0, mν)
127+ Ων0 = _ΩνE2 (1.0 , Ωγ0, mν):: promote_type (Float64, typeof (Ωγ0), typeof (mν) <: AbstractVector ? eltype (mν) : typeof (mν))
104128 ΩΛ0 = 1.0 - (Ωγ0 + Ωcb0 + Ων0 + Ωk0)
105- return a * 0.5 / (E_a (a, Ωcb0, h; mν= mν, w0= w0, wa= wa, Ωk0= Ωk0)^ 2 ) *
106- (- 3 (Ωcb0)a^- 4 - 4 Ωγ0 * a^- 5 - 2 Ωk0 * a^- 3 + ΩΛ0 * _dρDEda (a, w0, wa) + _dΩνE2da (a, Ωγ0, mν))
129+ if a isa AbstractArray
130+ return _dlogEdloga_scalar .(a, Ωcb0, Ωγ0, Ων0, ΩΛ0, Ωk0, h, Ref (mν), w0, wa)
131+ else
132+ return _dlogEdloga_scalar (a, Ωcb0, Ωγ0, Ων0, ΩΛ0, Ωk0, h, mν, w0, wa)
133+ end
107134end
108135
109136function _Ωma (a, Ωcb0, h; mν= 0.0 , w0= - 1.0 , wa= 0.0 , Ωk0= 0.0 )
110- return Ωcb0 * a^- 3 / (E_a (a, Ωcb0, h; mν= mν, w0= w0, wa= wa, Ωk0= Ωk0))^ 2
137+ return Ωcb0 . *. ^-3 . /E_a (a, Ωcb0, h; mν= mν, w0= w0, wa= wa, Ωk0= Ωk0)). ^ 2
111138end
112139
113140function _Ωma (a, cosmo:: w0waCDMCosmology )
@@ -222,16 +249,18 @@ end
222249function _growth_solver (z, Ωcb0, h; mν= 0.0 , w0= - 1.0 , wa= 0.0 , Ωk0= 0.0 )
223250 amin = 1 / 139
224251 loga = vcat (log .(_a_z .(z)))
225- u₀ = [amin, amin]
252+
253+ T = promote_type (eltype (z), typeof (Ωcb0), typeof (h), typeof (mν), typeof (w0), typeof (wa), typeof (Ωk0))
254+ u₀ = T[amin, amin]
226255
227- logaspan = (log (amin), log (1.01 ))# to ensure we cover the relevant range
256+ logaspan = (T ( log (amin)), T ( log (1.01 ) ))# to ensure we cover the relevant range
228257
229- p = [Ωcb0, mν, h, w0, wa, Ωk0]
258+ p = T [Ωcb0, mν, h, w0, wa, Ωk0]
230259
231- prob = ODEProblem (_growth!, u₀, logaspan, p)
260+ prob = ODEProblem {true} (_growth!, u₀, logaspan, p)
232261
233- sol = solve (prob, Tsit5 (), reltol= 1e-5 ; saveat= loga)[ 1 : 2 , :]
234- return sol
262+ sol = solve (prob, Tsit5 (), reltol= 1e-5 ; saveat= loga)
263+ return Array ( sol)[ 1 : 2 , :] :: Matrix{T}
235264end
236265
237266function D_z (z, Ωcb0, h; mν= 0.0 , w0= - 1.0 , wa= 0.0 , Ωk0= 0.0 )
0 commit comments