Skip to content

Commit

Permalink
fix RECFAST in 1.10 by removing struct and using NamedTuple
Browse files Browse the repository at this point in the history
  • Loading branch information
xzackli committed Aug 5, 2024
1 parent 6a42247 commit 6e5d29d
Show file tree
Hide file tree
Showing 2 changed files with 15 additions and 25 deletions.
38 changes: 14 additions & 24 deletions src/ionization/recfast.jl
Original file line number Diff line number Diff line change
Expand Up @@ -409,18 +409,8 @@ end_He_evo_condition(z, 𝕣) = (x_H0_H_Saha(𝕣, z) - 0.985)



struct RECFASTHistory{T, R, HES, HS}
𝕣::R
zinitial::T
zfinal::T
z_He_evo_start::T
z_H_He_evo_start::T
sol_He::HES
sol_H_He::HS
end


function Xe(rhist::RECFASTHistory, z)
function Xe_RECFAST(rhist, z)
𝕣 = rhist.𝕣
if (z > 8000.)
return 1 + 2*𝕣.fHe
Expand All @@ -442,7 +432,7 @@ function Xe(rhist::RECFASTHistory, z)
end
end

function Tmat(rhist::RECFASTHistory, z)
function Tmat_RECFAST(rhist, z)
𝕣 = rhist.𝕣
if (z > rhist.z_He_evo_start)
return Tmat_early(𝕣, z)
Expand Down Expand Up @@ -480,16 +470,16 @@ function recfastsolve(𝕣, alg=Tsit5(), zinitial=10000., zfinal=0.)
prob3 = ODEProblem{false}(ion_recfast, y3, (z3, zfinal), 𝕣)
sol_H_He = solve(prob3, alg, reltol=𝕣.tol)

return RECFASTHistory(𝕣, zinitial, zfinal,
return (; 𝕣, zinitial, zfinal,
z_He_evo_start, z_H_He_evo_start, sol_He, sol_H_He)
end


function reionization_Xe(rh::RECFASTHistory, z)
function reionization_Xe(rh, z)
𝕣 = rh.𝕣
X_fin = 1 + 𝕣.Yp / ( 𝕣.not4*(1-𝕣.Yp) ) #ionization frac today
zre,α,ΔH,zHe,ΔHe,fHe = 7.6711,1.5,0.5,3.5,0.5,X_fin-1 #reion params, TO REPLACE
x_orig = Xe(rh, z)
x_orig = Xe_RECFAST(rh, z)
x_reio_H = (X_fin - x_orig) / 2 * (
1 + tanh(( (1+zre)^α - (1+z)^α ) / ( α*(1+zre)^-1) ) / ΔH)) + x_orig
x_reio_He = fHe / 2 * ( 1 + tanh( (zHe - z) / ΔHe) )
Expand All @@ -498,7 +488,7 @@ function reionization_Xe(rh::RECFASTHistory, z)
end


function reionization_Tmat_ode(Tm, rh::RECFASTHistory, z)
function reionization_Tmat_ode(Tm, rh, z)
𝕣 = rh.𝕣
x_reio = reionization_Xe(rh, z)

Expand All @@ -519,17 +509,17 @@ struct TanhReionizationHistory{T, IH, TS}
end


function Xe(trhist::TanhReionizationHistory, z)
function Xe_TanhReio(trhist::TanhReionizationHistory, z)
if (z > trhist.zre_ini)
return Xe(trhist.ionization_history, z)
return Xe_RECFAST(trhist.ionization_history, z)
else
return reionization_Xe(trhist.ionization_history, z)
end
end

function Tmat(trhist::TanhReionizationHistory, z)
function Tmat_TanhReio(trhist::TanhReionizationHistory, z)
if (z > trhist.zre_ini)
return Tmat(trhist.ionization_history, z)
return Tmat_RECFAST(trhist.ionization_history, z)
else
return trhist.sol_reionization_Tmat(z)
end
Expand All @@ -538,7 +528,7 @@ end
function tanh_reio_solve(ion_hist, zre_ini=50.0)
𝕣 = ion_hist.𝕣
reio_prob = ODEProblem(reionization_Tmat_ode,
Tmat(ion_hist, zre_ini), (zre_ini, ion_hist.zfinal), ion_hist)
Tmat_RECFAST(ion_hist, zre_ini), (zre_ini, ion_hist.zfinal), ion_hist)
sol_reio_Tmat = solve(reio_prob, Tsit5(), reltol=𝕣.tol)
trh = TanhReionizationHistory(zre_ini, ion_hist, sol_reio_Tmat);

Expand Down Expand Up @@ -697,13 +687,13 @@ function IonizationHistory(𝕣::RECFAST{T}, par::AbstractCosmoParams{T}, bg::Ab


xinitial_RECFAST = z2x(rhist.zinitial)
Xe_initial = Xe(rhist, rhist.zinitial)
Xe_initial = Xe_RECFAST(rhist, rhist.zinitial)

Xₑ_function = x -> (x < xinitial_RECFAST) ?
Xe_initial : Xe(trhist, x2z(x))
Xe_initial : Xe_TanhReio(trhist, x2z(x))
Trad_function = x -> 𝕣.Tnow * (1 + x2z(x))
Tmat_function = x -> (x < xinitial_RECFAST) ?
Trad_function(x) : Tmat(trhist, x2z(x))
Trad_function(x) : Tmat_TanhReio(trhist, x2z(x))

# =====================================================
#j - do we really need bg to be passed to IonizationHistory separately from 𝕣.bg?
Expand Down
2 changes: 1 addition & 1 deletion test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ end
rhist = Bolt.recfastsolve(𝕣)
# xe_bespoke, Tmat_bespoke = Bolt.recfast_xe(𝕣; Nz=1000, zinitial=10000., zfinal=0.)
#change to only test pre-reion (z≧50)
@test all(abs.(Xe_fort .- Bolt.Xe.((rhist,), z⃗)) .< 1e-4)
@test all(abs.(Xe_fort .- Bolt.Xe_RECFAST.((rhist,), z⃗)) .< 1e-4)
end

##
Expand Down

0 comments on commit 6e5d29d

Please sign in to comment.