From efb271e3d41dab74d9f010983f42840b3f9ce765 Mon Sep 17 00:00:00 2001 From: glatterf42 Date: Thu, 12 May 2022 13:04:28 +0200 Subject: [PATCH] Adapted power spectra evaluation --- agora_plot_power_spectra.py | 56 +++++++++++++++++++++++++++++-------- plot_cross_spectra.py | 17 +++++------ plot_power_spectra.py | 8 +++--- 3 files changed, 57 insertions(+), 24 deletions(-) diff --git a/agora_plot_power_spectra.py b/agora_plot_power_spectra.py index 1da7bed..972f117 100644 --- a/agora_plot_power_spectra.py +++ b/agora_plot_power_spectra.py @@ -6,6 +6,7 @@ Created on Wed Mar 16 15:12:21 2022 @author: ben """ +from matplotlib import scale import pandas as pd import matplotlib.pyplot as plt from pathlib import Path @@ -15,21 +16,30 @@ import numpy as np def Dplus( lambda0, a ): return a * np.sqrt(1.0+lambda0*a**3) * sf.hyp2f1(3.0/2.0, 5.0/6.0, 11.0/6.0, -lambda0 * a**3 ) -basedir = Path("/home/ben/sims/swiftsim/examples/agora/spectra/") +basedir = Path("/home/ben/sims/swiftsim/examples/agora/") # basedir = Path("/home/ben/sims/gadget4/examples/agora_test/output/spectra/") #choose waveform and Lbox: Lbox = 85.47 # 60 #only option as of now -Nres1 = 128 +Nres1 = 512 # Nres2 = 256 k0 = 2 * 3.14159265358979323846264338327950 / Lbox knyquist1 = Nres1 * k0 # knyquist2 = Nres2 * k0 -times = [0.01, 0.02, 0.04, 0.08, 0.166666, 0.333333, 0.5, 0.666666, 1.0] +# #This is for 128 particles with fewer output times +# times = [0.01, 0.02, 0.04, 0.08, 0.166666, 0.333333, 0.5, 0.666666, 1.0] # times = [0.01, 0.02, 0.04, 0.08, 0.1] # times = [0.01, 0.02] # scale_factor = 0 # give index of a list above +#This is for 512 particles with more output times: +time_file = basedir / 'snap_times_agora1024' +df_times = pd.read_csv(f'{time_file}.txt', sep=' ', skipinitialspace=True, header=0, engine='python') + +output_numbers = [0, 1, 2, 4, 8, 16, 32] +redshifts = df_times.loc[output_numbers] + + Omega_m = 0.272 Omega_L = 0.728 lambda_0 = Omega_L / Omega_m @@ -42,7 +52,13 @@ columns = ["k [Mpc]", "Pcross", "P1", "err. P1", "P2", "err. P2", "P2-1", "err. zstart = 100 # 99.09174098173423 a_ics = 1 / (1 + zstart) -filename_ics = basedir / f'{Lbox}mpc/agora_ics_cross_spectrum' +filename_ics = basedir / f'spectra/{Lbox}mpc_{Nres1}/agora_ics_cross_spectrum' + +# #just as a test for normalising wrt to the first snapshot +# zstart = 12.0 +# a_ics = 1 / (1 + zstart) +# filename_ics = basedir / f'spectra/{Lbox}mpc_{Nres1}/agora_0000_cross_spectrum' + # filename_ics = basedir / 'agora_a0_cross_spectrum' df_ics = pd.read_csv(f'{filename_ics}.txt', sep=' ', skipinitialspace=True, header=None, names=columns, skiprows=1) #only consider rows above resolution limit @@ -53,17 +69,26 @@ Dplus0 = Dplus(lambda0=lambda_0, a=a_ics) D_squared_ics = Dplus0 ** 2 p1_ics_noramlised = p1_ics / D_squared_ics +print(D_squared_ics) + # p1_ics_ic_normalised = p1_ics_noramlised / p1_ics_noramlised # plt.loglog(k_ics, p1_ics_ic_normalised, label='ICs') -for scale_factor in range(len(times)): - filename = basedir / f"{Lbox}mpc/agora_a{scale_factor}_cross_spectrum" +# #This is for 128 particles with fewer output times: +# for scale_factor in range(len(times)): +# filename = basedir / f"{Lbox}mpc/agora_a{scale_factor}_cross_spectrum" # filename = basedir / f"agora_a{scale_factor}_small_dt_cross_spectrum" # filename = basedir / f'agora_a{scale_factor}_cross_spectrum' # filename = basedir / f"{waveform}_{Lbox:.0f}/{waveform}_{Lbox:.0f}_ics_vsc_cross_spectrum" # for ICs # savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/power_{waveform}_{Lbox:.0f}_ics_vsc") # for ICs # plt.title(f"Power Spectra {waveform} L={Lbox:.0f} a=0.02 vsc") # for ICs +#This is for 512 particles with more output times: +for k in range(len(output_numbers)): + output_number = output_numbers[k] + redshift = float(redshifts.loc[output_number]) + scale_factor = 1 / (1 + redshift) + filename = basedir / f'spectra/{Lbox}mpc_{Nres1}/agora_{output_number:04d}_cross_spectrum' df = pd.read_csv(f"{filename}.txt", sep=" ", skipinitialspace=True, header=None, names=columns, skiprows=1) @@ -73,20 +98,26 @@ for scale_factor in range(len(times)): k = df["k [Mpc]"] p1 = df["P1"] p1_error = df["err. P1"] + # p2 = df["P2"] # p2_error = df["err. P2"] # pcross = df["Pcross"] - D_squared = Dplus(lambda0=lambda_0, a=times[scale_factor]) ** 2 + # D_squared = Dplus(lambda0=lambda_0, a=times[scale_factor]) ** 2 + D_squared = Dplus(lambda0=lambda_0, a=scale_factor) ** 2 + p1_normalised = p1 / D_squared p1_ic_normalised = p1_normalised / p1_ics_noramlised + print(D_squared) + # Plot the power spectra: - plt.loglog(k, p1_ic_normalised, label=f"{times[scale_factor]}") + # plt.loglog(k, p1_ic_normalised, label=f"{times[scale_factor]}") + plt.loglog(k, p1_ic_normalised, label=f"{scale_factor}") # plt.loglog(k, p2, label="P2") -plt.title(f"Power Spectra Agora 128") -savedir = Path(f"/home/ben/Pictures/swift/agora/spectra/{Lbox}mpc/") +plt.title(f"Power Spectra Agora {Nres1}") +savedir = Path(f"/home/ben/Pictures/swift/agora/spectra/{Lbox}mpc_{Nres1}") plt.xlabel(r"k [$\mathrm{Mpc}^{-1}$]") plt.ylabel("P") @@ -94,10 +125,11 @@ plt.vlines(knyquist1, ymin=min(p1_ic_normalised), ymax=max(p1_ic_normalised), co # plt.vlines(knyquist2, ymin=min(p2), ymax=max(p2), color="black", linestyles="dashed", label=f"{Nres2}") plt.legend() # plt.ylim(1, 3) + +plt.savefig(f"{savedir}/agora_{Nres1}_power.png") + plt.show() -# plt.savefig(f"{savedir}_2.png") - diff --git a/plot_cross_spectra.py b/plot_cross_spectra.py index dad9726..cddf29b 100644 --- a/plot_cross_spectra.py +++ b/plot_cross_spectra.py @@ -15,16 +15,17 @@ basedir = Path("/home/ben/sims/data_swift/monofonic_tests/spectra/") #choose Nres and Lbox: waveforms = ['DB2', "DB4", "DB8", "shannon"] #DB2, DB4, DB8, shannon are all we have right now Lbox = 100.0 #only option as of now -Nres = 256 #128 and 256 exist for now +Nres1 = 128 #Nres1 should always be smaller than Nres2 +Nres2 = 256 #128, 256 and 512 exist for now k0 = 2 * 3.14159265358979323846264338327950 / Lbox -knyquist = Nres * k0 +knyquist = Nres2 * k0 #Not used at the moment anyway except for upper limit, for which we need the larger Nres a = [0.166666, 0.333333, 0.5, 0.666666, 1.0] scale_factor = 4 # give index of a list above for wave in waveforms: - # filename = basedir / f"{wave}_{Lbox:.0f}/{wave}_{Lbox:.0f}_a{scale_factor}_cross_spectrum" - filename = basedir / f"{wave}_{Lbox:.0f}/{wave}_{Lbox:.0f}_ics_local_cross_spectrum" # for ICs + # filename = basedir / f"{wave}_{Lbox:.0f}/{wave}_{Lbox:.0f}_a{scale_factor}_{Nres1}_{Nres2}_cross_spectrum" + filename = basedir / f"{wave}_{Lbox:.0f}/{wave}_{Lbox:.0f}_ics_{Nres1}_{Nres2}_cross_spectrum" # for ICs #find columns in file manually #is k really in Mpc? Swift doesn't use /h internally at least. @@ -46,11 +47,11 @@ for wave in waveforms: # Plot the Cross Correlation: plt.plot(k, pcross, label=f"{wave}") -savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/cross_{Nres}_{Lbox:.0f}_ics_local") # for ICs -plt.title(f"Cross correlation N={Nres} L={Lbox:.0f} a=0.02") # for ICs +savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/cross_{Nres1}_{Nres2}_{Lbox:.0f}_ics") # for ICs +plt.title(f"Cross correlation N=({Nres1}, {Nres2}) L={Lbox:.0f} a=0.02") # for ICs -# savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/cross_{Nres}_{Lbox:.0f}_a{scale_factor}") -# plt.title(f"Cross correlation N={Nres} L={Lbox:.0f} a={a[scale_factor]}") +# savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/cross_{Nres1}_{Nres2}_{Lbox:.0f}_a{scale_factor}") +# plt.title(f"Cross correlation N=({Nres1}, {Nres2}) L={Lbox:.0f} a={a[scale_factor]}") plt.xscale("log") plt.xlabel(r"k [$\mathrm{Mpc}^{-1}$]") diff --git a/plot_power_spectra.py b/plot_power_spectra.py index 3a46db8..1a412d3 100644 --- a/plot_power_spectra.py +++ b/plot_power_spectra.py @@ -10,20 +10,20 @@ import pandas as pd import matplotlib.pyplot as plt from pathlib import Path -basedir = Path("/home/ben/sims/swift/monofonic_tests/spectra/") +basedir = Path("/home/ben/sims/data_swift/monofonic_tests/spectra/") #choose waveform and Lbox: waveform = "shannon" #DB2, DB4, DB8 or shannon Lbox = 100.0 #only option as of now Nres1 = 128 -Nres2 = 256 +Nres2 = 512 k0 = 2 * 3.14159265358979323846264338327950 / Lbox knyquist1 = Nres1 * k0 knyquist2 = Nres2 * k0 a = [0.166666, 0.333333, 0.5, 0.666666, 1.0] scale_factor = 4 # give index of a list above -filename = basedir / f"{waveform}_{Lbox:.0f}/{waveform}_{Lbox:.0f}_a{scale_factor}_cross_spectrum" +filename = basedir / f"{waveform}_{Lbox:.0f}/{waveform}_{Lbox:.0f}_a{scale_factor}_{Nres1}_{Nres2}_cross_spectrum" # filename = basedir / f"{waveform}_{Lbox:.0f}/{waveform}_{Lbox:.0f}_ics_vsc_cross_spectrum" # for ICs # savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/power_{waveform}_{Lbox:.0f}_ics_vsc") # for ICs # plt.title(f"Power Spectra {waveform} L={Lbox:.0f} a=0.02 vsc") # for ICs @@ -50,7 +50,7 @@ plt.loglog(k, p1, label="P1") plt.loglog(k, p2, label="P2") plt.title(f"Power Spectra {waveform} L={Lbox:.0f} a={a[scale_factor]}") -savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/power_{waveform}_{Lbox:.0f}_a{scale_factor}") +savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/power_{waveform}_{Lbox:.0f}_{Nres1}_{Nres2}_a{scale_factor}") plt.xlabel("k [Mpc]") plt.ylabel("P")