Skip to content

Commit

Permalink
Merge 943e027 into 6cc7788
Browse files Browse the repository at this point in the history
  • Loading branch information
@vitenti
vitenti committed Mar 22, 2024
2 parents 6cc7788 + 943e027 commit f557100
Show file tree
Hide file tree
Showing 44 changed files with 8,574 additions and 3,641 deletions.
18 changes: 10 additions & 8 deletions docs/numcosmo-docs.sgml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,9 +77,10 @@
<section>
<title>Power spectrum functions</title>
<xi:include href="xml/ncm_powspec.xml"/>
<xi:include href="xml/ncm_powspec_corr3d.xml"/>
<xi:include href="xml/ncm_powspec_filter.xml"/>
<xi:include href="xml/ncm_powspec_sphere_proj.xml"/>
<xi:include href="xml/ncm_powspec_corr3d.xml"/>
<xi:include href="xml/ncm_powspec_spline2d.xml"/>
</section>
<section>
<title>Harmonic Oscillator</title>
Expand Down Expand Up @@ -271,14 +272,15 @@
<chapter>
<title>Perturbations</title>
<xi:include href="xml/nc_hipert.xml"/>
<xi:include href="xml/nc_hipert_bg_var.xml"/>
<xi:include href="xml/nc_hipert_wkb.xml"/>
<xi:include href="xml/nc_hipert_adiab.xml"/>
<xi:include href="xml/nc_hipert_gw.xml"/>
<xi:include href="xml/nc_hipert_two_fluids.xml"/>
<xi:include href="xml/nc_hipert_boltzmann.xml"/>
<xi:include href="xml/nc_hipert_bg_var.xml"/>
<xi:include href="xml/nc_hipert_boltzmann_std.xml"/>
<xi:include href="xml/nc_hipert_boltzmann.xml"/>
<xi:include href="xml/nc_hipert_em.xml"/>
<xi:include href="xml/nc_hipert_first_order.xml"/>
<xi:include href="xml/nc_hipert_gw.xml"/>
<xi:include href="xml/nc_hipert_two_fluids.xml"/>
<xi:include href="xml/nc_hipert_wkb.xml"/>
<section>
<title>Gravitational theories</title>
<xi:include href="xml/nc_hipert_grav.xml"/>
Expand Down Expand Up @@ -320,9 +322,9 @@
<section>
<title>Matter Power Spectrum</title>
<xi:include href="xml/nc_powspec_ml.xml"/>
<xi:include href="xml/nc_powspec_ml_fix_spline.xml"/>
<xi:include href="xml/nc_powspec_ml_transfer.xml"/>
<xi:include href="xml/nc_powspec_ml_cbe.xml"/>
<xi:include href="xml/nc_powspec_ml_spline.xml"/>
<xi:include href="xml/nc_powspec_ml_transfer.xml"/>
<xi:include href="xml/nc_powspec_mnl.xml"/>
<xi:include href="xml/nc_powspec_mnl_halofit.xml"/>
</section>
Expand Down
150 changes: 27 additions & 123 deletions examples/example_zeta_exp_potential.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,16 +55,16 @@ def test_zeta_exp_potential(new_config: bool = True) -> None:
Vexp.props.Omegac = 1.0
Vexp.props.H0 = 67.8
else:
Vexp.props.alphab = +1.0e-20
Vexp.props.sigmaphi = +8.0e-1
Vexp.props.dphi = +5.0e-1
Vexp.props.xb = 2.0e38
Vexp.props.alphab = +8.3163e-2
Vexp.props.sigmaphi = +9.0
Vexp.props.dphi = -9.0e-4
Vexp.props.xb = 2.0e36
Vexp.props.OmegaL = 1.0
Vexp.props.Omegac = 1.0
Vexp.props.H0 = 67.8

k = 1.0e0
tc = Vexp.tau_xe(1.0e15)
k = 0.05
tc = +Vexp.tau_xe(1.0e15)
reltol = 1.0e-14

adiab.set_ti(Vexp.tau_min())
Expand All @@ -84,151 +84,55 @@ def test_zeta_exp_potential(new_config: bool = True) -> None:
# exit ()

print("# Preparing ADIAB")
found, adiab_ti = adiab.find_adiab_time_limit(Vexp, -20.0, -1.0e-2, 1.0e-6)
assert found, "Cannot find initial time"
adiab.set_init_cond_adiab(Vexp, adiab_ti)
adiab.prepare(Vexp)

# exit ()

print("# Preparing GW")
found, gw_ti = gw.find_adiab_time_limit(Vexp, -20.0, -1.0e-2, 1.0e-8)
assert found, "Cannot find initial time"
gw.set_init_cond_adiab(Vexp, gw_ti)
gw.prepare(Vexp)

# (t0, t1) = adiab.get_t0_t1 (Vexp)
(t0, t1) = gw.get_t0_t1(Vexp)
t0 = gw.get_ti()
t1 = gw.get_tf()

print(f"# BACKG (t0, t1) = ({Vexp.tau_min(): 21.15e}, {Vexp.tau_max(): 21.15e})")
print(f"# ADIAB (t0, t1) = ({t0: 21.15e}, {t1: 21.15e})")

(Delta_zeta_c, _Delta_Pzeta_c) = adiab.eval_Delta(Vexp, tc)
(Delta_h_c, _Delta_Ph_c) = gw.eval_Delta(Vexp, tc)
(Delta_zeta_c, _Delta_Pzeta_c) = adiab.eval_powspec_at(Vexp, tc)
(Delta_h_c, _Delta_Ph_c) = gw.eval_powspec_at(Vexp, tc)

print(f"# Time of x_e = 10^15: tau_c = {tc: 21.15f}")
print(f"# Power spectrum at tau_c: PS_ADIAB = {Delta_zeta_c: 21.15e}")
print(f"# Power spectrum at tau_c: PS_GW = {Delta_h_c: 21.15e}")
print(
f"# Power spectrum at tau_c: r = {2.0 * Delta_h_c / Delta_zeta_c: 21.15e}"
)

tau_a = np.linspace(-5.0, t1, 100000)
Delta_zeta = np.array([adiab.eval_powspec_at(Vexp, t) for t in tau_a])
Delta_h = np.array([gw.eval_powspec_at(Vexp, t) for t in tau_a])
mylw = 1.0

t_a = np.linspace(t0, t1, 100000)
alpha_a = []

plt.plot(tau_a, np.abs(Delta_zeta[:,0]), lw=mylw, label=r"$\Delta_{\zeta}a$")
plt.plot(tau_a, Delta_h[:,0], lw=mylw, label=r"$\Delta_h$")

Delta_zeta = []
# Delta_Pzeta = []

Delta_h = []
# Delta_Ph = []

zeta_a_a = []
zeta_b_a = []
Pzeta_a_a = []
Pzeta_b_a = []

h_a_a = []
h_b_a = []
Ph_a_a = []
Ph_b_a = []

nu_a = []
m_a = []
mnu_a = []
dlnmnu_a = []

epsilon_a = []
gamma_a = []
sin_thetab_a = []
cos_thetab_a = []

for t in t_a:
(Delta_zeta_v, _Delta_Pzeta_v) = adiab.eval_Delta(Vexp, t)
(Delta_h_v, _Delta_Ph_v) = gw.eval_Delta(Vexp, t)

(zeta_a_v, zeta_b_v, Pzeta_a_v, Pzeta_b_v) = adiab.eval_QV(Vexp, t)
(h_a_v, h_b_v, Ph_a_v, Ph_b_v) = gw.eval_QV(Vexp, t)

(epsilon_v, gamma_v, sin_thetab_v, cos_thetab_v) = gw.eval_AA(Vexp, t)

epsilon_a.append(epsilon_v)
gamma_a.append(gamma_v)
sin_thetab_a.append(sin_thetab_v)
cos_thetab_a.append(cos_thetab_v)

nu_adiab = adiab.eval_nu(Vexp, t, k)
mnu_adiab = adiab.eval_mnu(Vexp, t, k)
m_adiab = mnu_adiab / nu_adiab
dlnmnu_adiab = math.fabs(adiab.eval_dlnmnu(Vexp, t, k))

alpha = 0
if t > 0.0:
alpha = +0.5 * t**2 / math.log(10.0)
else:
alpha = -0.5 * t**2 / math.log(10.0)

alpha_a.append(alpha)

nu_a.append(nu_adiab)
m_a.append(m_adiab)
mnu_a.append(mnu_adiab)
dlnmnu_a.append(dlnmnu_adiab)

Delta_zeta.append(Delta_zeta_v)
Delta_h.append(Delta_h_v)

zeta_a_a.append(zeta_a_v)
zeta_b_a.append(zeta_b_v)

Pzeta_a_a.append(Pzeta_a_v)
Pzeta_b_a.append(Pzeta_b_v)

h_a_a.append(h_a_v)
h_b_a.append(h_b_v)

Ph_a_a.append(Ph_a_v)
Ph_b_a.append(Ph_b_v)

mylw = 1

h_a_npa = np.array(h_a_a)
h_b_npa = np.array(h_b_a)

# Ph_a_npa = np.array(Ph_a_a)
# Ph_b_npa = np.array(Ph_b_a)

zeta_a_npa = np.array(zeta_a_a)
zeta_b_npa = np.array(zeta_b_a)

# Pzeta_a_npa = np.array(Pzeta_a_a)
# Pzeta_b_npa = np.array(Pzeta_b_a)

# plt.plot (alpha_a, Delta_zeta, lw=mylw, label = r'$\Delta_{\zeta}$')
# plt.plot (alpha_a, Delta_h, lw=mylw, label = r'$\Delta_{h}$')

plt.plot(alpha_a, zeta_a_npa, lw=mylw, label=r"$\tilde{\zeta}^a$")
plt.plot(alpha_a, zeta_b_npa, lw=mylw, label=r"$\tilde{\zeta}^b$")
plt.plot(alpha_a, h_a_npa, lw=mylw, label=r"$\tilde{h}^a$")
plt.plot(alpha_a, h_b_npa, lw=mylw, label=r"$\tilde{h}^b$")

# plt.plot (alpha_a, epsilon_a, lw=mylw, label = r'$\epsilon$')
# plt.plot (alpha_a, gamma_a, lw=mylw, label = r'$\gamma$')
# plt.plot (alpha_a, sin_thetab_a, lw=mylw, label = r'$\sin(\theta_b)$')
# plt.plot (alpha_a, cos_thetab_a, lw=mylw, label = r'$\cos(\theta_b)$')

# plt.plot (alpha_a, nu_a, lw=mylw, label = r'$\nu_h$')
# plt.plot (alpha_a, m_a, lw=mylw, label = r'$m_h$')
# plt.plot (alpha_a, mnu_a, lw=mylw, label = r'$m_h\nu_h$')
# plt.plot (alpha_a, 1.0 / (np.array (mnu_a)), lw=mylw, label = r'$(m_h\nu_h)^{-1}$')
# plt.plot (alpha_a, dlnmnu_a, lw=mylw, label = r'$d\ln(m_h\nu_h)$')

# plt.plot (alpha_a, (Ph_a_a * h_b_a), lw=mylw, label = r't_1')
# plt.plot (alpha_a, (Ph_b_a * h_a_a), lw=mylw, label = r't_2')

plt.grid(b=True, which="both", linestyle=":", color="0.75", linewidth=0.5)
plt.grid(which="both", linestyle=":", color="0.75", linewidth=0.5)
plt.legend(loc="upper left")

# plt.xscale('symlog', linthreshx=1.0e-5)
plt.yscale("symlog", linthreshy=1.0e23)
plt.yscale("log")
# plt.yscale('log', linthreshy=1.0e-20)

plt.show()
plt.clf()


if __name__ == "__main__":
test_zeta_exp_potential()
test_zeta_exp_potential(new_config = False)
30 changes: 15 additions & 15 deletions notebooks/primordial_perturbations/magnetic_dust_bounce.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,21 +76,13 @@
"metadata": {},
"outputs": [],
"source": [
"try:\n",
" import gi\n",
"\n",
" gi.require_version(\"NumCosmo\", \"1.0\")\n",
" gi.require_version(\"NumCosmoMath\", \"1.0\")\n",
"except:\n",
" pass\n",
"\n",
"import sys\n",
"import math\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from gi.repository import NumCosmo as Nc\n",
"from gi.repository import NumCosmoMath as Ncm"
"\n",
"from numcosmo_py import Nc, Ncm"
]
},
{
Expand Down Expand Up @@ -153,7 +145,7 @@
"outputs": [],
"source": [
"class PyCSQ1DMagDust(Ncm.CSQ1D):\n",
" def __init__(self, Omega_m=0.3, xb=1.0e25, cc=1.0, h=0.7, a_ns=1.0e-12):\n",
" def __init__(self, Omega_m=0.3, xb=1.0e25, cc=1.0e-50, h=0.7, a_ns=1.0e-12):\n",
" Ncm.CSQ1D.__init__(self)\n",
"\n",
" self.h = h\n",
Expand Down Expand Up @@ -312,9 +304,10 @@
"kf = 4.0e3\n",
"k_a = np.geomspace(ki, kf, 3)\n",
"\n",
"csq1d.set_k(kf)\n",
"(Found2, tf) = csq1d.find_adiab_time_limit(None, +1.0e-25, +1.0e15, 1.0e0)\n",
"csq1d.set_max_order_2(False)\n",
"\n",
"csq1d.set_k(ki)\n",
"(Found2, tf) = csq1d.find_adiab_time_limit(None, +1.0e-25, +1.0e15, 1.0e0)\n",
"print(Found2, tf)"
]
},
Expand All @@ -340,7 +333,7 @@
" # print (csq1d.eval_adiab_at (None, -5))\n",
"\n",
" (Found1, ti) = csq1d.find_adiab_time_limit(None, -1.0e15, -1.0e-25, 1.0e-3)\n",
" # print (Found1, ti)\n",
" print (Found1, ti)\n",
"\n",
" (Found2, tfa) = csq1d.find_adiab_time_limit(None, +1.0e-25, +1.0e15, 1.0e0)\n",
" tf = tfa * 20\n",
Expand Down Expand Up @@ -647,6 +640,13 @@
"print(\"x_eq = % 22.15e\" % (xeq))\n",
"print(\"t_eq = % 22.15e\" % (teq))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand All @@ -668,5 +668,5 @@
}
},
"nbformat": 4,
"nbformat_minor": 2
"nbformat_minor": 4
}
Loading

0 comments on commit f557100

Please sign in to comment.