formatting; added pbh box length file forfuture reference; spectra_plot can now plot power spectra of z=1 INSTEAD of end, just comment/uncomment the lines with _z1 and _end

find_box_length_for_pbh.py

@@ -0,0 +1,131 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jan 13 17:19:55 2022
+
+@author: Ben Melville
+"""
+
+import numpy as np
+import scipy.integrate as integrate
+import matplotlib.pyplot as plt
+
+# Cosmological parameters:
+n_s = 0.97
+Omega_M_0 = 0.3099
+Omega_Lambda_0 = 0.690021
+Omega_R_0 = 0.0
+h = 0.67742
+R = 40  # for R=40, sigma is only 0.113, R>40 yields nan in one integration
+z = 0
+prop_factor = 289.98204077151246  # found from running R=8, z=0, requiring that sigma=0.8 should hold
+
+# L_box / N_res = 0.009235859375        #for 10^5 M_sun particle mass -> 4096 almost enough to go to z=0
+# L_box / N_res = 0.019897945879081797  #for 10^6 M_sun particle mass -> 2048 just enough to go to z=0
+
+
+def E_squared(z, Omega_Lambda_0, Omega_M_0, Omega_R_0):
+    Omega_0 = Omega_Lambda_0 + Omega_M_0 + Omega_R_0
+    return (
+        Omega_Lambda_0
+        + (1 - Omega_0) * (1 + z) ** 2
+        + Omega_M_0 * (1 + z) ** 3
+        + Omega_R_0 * (1 + z) ** 4
+    )
+
+
+def Omega_M(z, Omega_Lambda_0, Omega_M_0, Omega_R_0):
+    return Omega_M_0 * (1 + z) ** 3 / E_squared(z, Omega_Lambda_0, Omega_M_0, Omega_R_0)
+
+
+def Omega_Lambda(z, Omega_Lambda_0, Omega_M_0, Omega_R_0):
+    return Omega_Lambda_0 / E_squared(z, Omega_Lambda_0, Omega_M_0, Omega_R_0)
+
+
+def g(z, Omega_Lambda_0, Omega_M_0, Omega_R_0):
+    Omega_M_local = Omega_M(z, Omega_Lambda_0, Omega_M_0, Omega_R_0)
+    Omega_Lambda_local = Omega_Lambda(z, Omega_Lambda_0, Omega_M_0, Omega_R_0)
+    return (
+        5
+        * Omega_M_local
+        / (
+            2
+            * (
+                Omega_M_local ** (4 / 7)
+                - Omega_Lambda_local
+                + (1 + Omega_M_local / 2) * (1 + Omega_Lambda_local / 70)
+            )
+        )
+    )
+
+
+def D(z, Omega_Lambda_0, Omega_M_0, Omega_R_0):
+    return g(z, Omega_Lambda_0, Omega_M_0, Omega_R_0) / (1 + z)
+
+
+def keq(Omega_M_0, h):
+    return Omega_M_0 * h ** 2
+
+
+def T(k, k_eq):
+    if k < k_eq:
+        return 1
+    elif k >= k_eq:
+        return np.log(k / k_eq) / (k / k_eq) ** 2
+
+
+def P_i(k, n_s):
+    return k ** n_s
+
+
+def P(k, z, n_s, Omega_M_0, Omega_Lambda_0, Omega_R_0, h):
+    k_eq = keq(Omega_M_0, h)
+    return (
+        P_i(k, n_s) * T(k, k_eq) ** 2 * D(z, Omega_Lambda_0, Omega_M_0, Omega_R_0) ** 2
+    )
+
+
+def W_R(k, R):
+    argument = k * R
+    return 3 / argument ** 2 * (np.sin(argument) - argument * np.cos(argument))
+
+
+def integrand(k, z, R, n_s, Omega_M_0, Omega_Lambda_0, Omega_R_0, h):
+    return P(k, z, n_s, Omega_M_0, Omega_Lambda_0, Omega_R_0, h) * W_R(k, R) * k ** 2
+
+
+k_lower = 0.0
+k_eq = keq(Omega_M_0, h)
+k_higher = 100 * k_eq
+
+integral, error = integrate.quad(
+    integrand,
+    k_lower,
+    k_eq,
+    args=(z, R, n_s, Omega_M_0, Omega_Lambda_0, Omega_R_0, h),
+    limit=100,
+)
+integral_2, error_2 = integrate.quad(
+    integrand,
+    k_eq,
+    k_higher,
+    args=(z, R, n_s, Omega_M_0, Omega_Lambda_0, Omega_R_0, h),
+    limit=100,
+)
+
+sigma_squared = (integral + integral_2) / (2 * np.pi ** 2)
+sigma = np.sqrt(sigma_squared) * prop_factor
+
+print(sigma)
+
+
+# k_plot = np.logspace(np.log10(k_lower), np.log10(k_eq), 1000)
+# T_plot = [T(k, k_eq) for k in k_plot]
+
+# k_plot_2 = np.logspace(np.log10(k_eq), np.log10(k_higher), 1000)
+# T_plot_2 = [T(k, k_eq) for k in k_plot_2]
+
+# plt.title("Linear transfer function")
+# plt.xlabel("log(k)")
+# plt.ylabel("log(T(k))")
+# plt.loglog(k_plot, T_plot, ".", markersize=5)
+# plt.loglog(k_plot_2, T_plot_2, ".", markersize=5)

music_shift_analysis_v3.py

@@ -53,5 +53,5 @@ print(
     f"Highest centre values for same Music shift: {(centre_music_ics + maximum_positive_boundary_for_same_shift) * 100}"
 )
 print(
-    f'Difference between these limits: {(maximum_positive_boundary_for_same_shift + maximum_negative_boundary_for_same_shift) * 100}'
+    f"Difference between these limits: {(maximum_positive_boundary_for_same_shift + maximum_negative_boundary_for_same_shift) * 100}"
 )

spectra_plot.py

@@ -47,6 +47,15 @@ def spectra_data(
             names=columns,
             skiprows=1,
         )
+    elif time == "z=1":
+        spectra_data = pd.read_csv(
+            f"{dir}/{waveform}_{Lbox}_a2_{resolution_1}_{resolution_2}_cross_spectrum.txt",
+            sep=" ",
+            skipinitialspace=True,
+            header=None,
+            names=columns,
+            skiprows=1,
+        )
     elif time == "end":
         spectra_data = pd.read_csv(
             f"{dir}/{waveform}_{Lbox}_a4_{resolution_1}_{resolution_2}_cross_spectrum.txt",
@@ -57,7 +66,7 @@ def spectra_data(
             skiprows=1,
         )
     else:
-        raise ValueError(f"invalid time ({time}) should be (ics|end)")
+        raise ValueError(f"invalid time ({time}) should be (ics|z=1|end)")
 
     # only consider rows above resolution limit
     spectra_data = spectra_data[spectra_data["k [Mpc]"] >= k0]
@@ -74,8 +83,10 @@ def create_plot(mode):
     )
     for i, waveform in enumerate(waveforms):
         ax_ics: Axes = axes[i][0]
+        # ax_z1:  Axes = axes[i][1]
         ax_end: Axes = axes[i][1]
         bottom_row = i == len(waveforms) - 1
+        # for is_end, ax in enumerate([ax_ics, ax_z1]):
         for is_end, ax in enumerate([ax_ics, ax_end]):
             if bottom_row:
                 ax.set_xlabel("k [Mpc$^{-1}$]")
@@ -91,6 +102,7 @@ def create_plot(mode):
             ax.text(
                 0.98 if mode == "cross" else 0.93,
                 0.85,
+                # "z=1" if is_end else "ics",
                 "end" if is_end else "ics",
                 size=13,
                 horizontalalignment="right",
@@ -124,9 +136,18 @@ def create_plot(mode):
                 comp_p1 = comp_data["P1"]
                 end_p1 /= comp_p1
 
+                # z1_data = spectra_data(waveform, resolution, resolution, Lbox, "z=1")
+                # z1_k = z1_data["k [Mpc]"]
+                # z1_p1 = z1_data["P1"]
+                # comp_data = spectra_data(waveform, 512, 512, Lbox, 'z=1')
+                # comp_p1 = comp_data["P1"]
+                # z1_p1 /= comp_p1
+                # ax_z1.semilogx(z1_k, z1_p1, color=colors[j])
+
                 ax_ics.semilogx(ics_k, ics_p1, color=colors[j])
                 ax_end.semilogx(end_k, end_p1, color=colors[j])
                 for ax in [ax_ics, ax_end]:
+                # for ax in [ax_ics, ax_z1]:
                     ax.set_ylim(0.9, 1.10)
 
 
@@ -152,6 +173,7 @@ def create_plot(mode):
             ax_end.set_xlim(right=k0 * resolutions[-1])
             ax_end.set_ylim(0.8, 1.02)
         if bottom_row:
+            # ax_z1.legend()
             ax_end.legend()
 
         # fig.suptitle(f"Cross Spectra {time}") #Not needed for paper