From 3a79c4f231ce965273576642cc94fb4803971471 Mon Sep 17 00:00:00 2001 From: glatterf42 Date: Tue, 5 Jul 2022 10:10:28 +0200 Subject: [PATCH] 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 | 131 +++++++++++++++++++++++++++++++++++++ music_shift_analysis_v3.py | 2 +- spectra_plot.py | 24 ++++++- 3 files changed, 155 insertions(+), 2 deletions(-) create mode 100644 find_box_length_for_pbh.py diff --git a/find_box_length_for_pbh.py b/find_box_length_for_pbh.py new file mode 100644 index 0000000..f4bd43f --- /dev/null +++ b/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) diff --git a/music_shift_analysis_v3.py b/music_shift_analysis_v3.py index 8a60d39..aed193e 100644 --- a/music_shift_analysis_v3.py +++ b/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}" ) diff --git a/spectra_plot.py b/spectra_plot.py index e03b9b4..c9e5cc0 100644 --- a/spectra_plot.py +++ b/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