diff --git a/__pycache__/write_spectra_jobs.cpython-38.pyc b/__pycache__/write_spectra_jobs.cpython-38.pyc index 9fca6ad..4a08f64 100644 Binary files a/__pycache__/write_spectra_jobs.cpython-38.pyc and b/__pycache__/write_spectra_jobs.cpython-38.pyc differ diff --git a/agora_plot_power_spectra.py b/agora_plot_power_spectra.py index 6959d09..1da7bed 100644 --- a/agora_plot_power_spectra.py +++ b/agora_plot_power_spectra.py @@ -16,37 +16,50 @@ 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/gadget4/examples/agora_test/output/spectra/") #choose waveform and Lbox: -Lbox = 100.0 #only option as of now +Lbox = 85.47 # 60 #only option as of now Nres1 = 128 # Nres2 = 256 k0 = 2 * 3.14159265358979323846264338327950 / Lbox knyquist1 = Nres1 * k0 # knyquist2 = Nres2 * k0 -times = [0.166666, 0.333333, 0.5, 0.666666, 1.0] +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 Omega_m = 0.272 Omega_L = 0.728 lambda_0 = Omega_L / Omega_m +# print(Dplus(lambda_0, 1.0)) + #find columns in file manually #is k really in Mpc? Swift doesn't use /h internally at least. columns = ["k [Mpc]", "Pcross", "P1", "err. P1", "P2", "err. P2", "P2-1", "err. P2-1", "modes in bin"] -zstart = 100 +zstart = 100 # 99.09174098173423 a_ics = 1 / (1 + zstart) -filename_ics = basedir / 'agora_ics_cross_spectrum' +filename_ics = basedir / f'{Lbox}mpc/agora_ics_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 df_ics = df_ics[df_ics['k [Mpc]'] >= k0] +k_ics = df_ics['k [Mpc]'] p1_ics = df_ics['P1'] -D_squared_ics = Dplus(lambda0=lambda_0, a=a_ics) ** 2 +Dplus0 = Dplus(lambda0=lambda_0, a=a_ics) +D_squared_ics = Dplus0 ** 2 p1_ics_noramlised = p1_ics / 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"agora_a{scale_factor}_cross_spectrum" + 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 @@ -73,15 +86,17 @@ for scale_factor in range(len(times)): # plt.loglog(k, p2, label="P2") plt.title(f"Power Spectra Agora 128") -savedir = Path(f"/home/ben/Pictures/swift/agora/spectra/power_spectra") +savedir = Path(f"/home/ben/Pictures/swift/agora/spectra/{Lbox}mpc/") -plt.xlabel("k [Mpc]") +plt.xlabel(r"k [$\mathrm{Mpc}^{-1}$]") plt.ylabel("P") -plt.vlines(knyquist1, ymin=min(p1), ymax=max(p1), color="black", linestyles="dashed", label=f"k_ny {Nres1}") +plt.vlines(knyquist1, ymin=min(p1_ic_normalised), ymax=max(p1_ic_normalised), color="black", linestyles="dashed", label=f"k_ny {Nres1}") # plt.vlines(knyquist2, ymin=min(p2), ymax=max(p2), color="black", linestyles="dashed", label=f"{Nres2}") plt.legend() +# plt.ylim(1, 3) +plt.show() -plt.savefig(f"{savedir}.png") +# plt.savefig(f"{savedir}_2.png") diff --git a/plot_cross_spectra.py b/plot_cross_spectra.py index b49be87..dad9726 100644 --- a/plot_cross_spectra.py +++ b/plot_cross_spectra.py @@ -10,7 +10,7 @@ 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 Nres and Lbox: waveforms = ['DB2', "DB4", "DB8", "shannon"] #DB2, DB4, DB8, shannon are all we have right now @@ -23,8 +23,8 @@ 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_vsc_cross_spectrum" # for ICs + # 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 #find columns in file manually #is k really in Mpc? Swift doesn't use /h internally at least. @@ -46,14 +46,14 @@ 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_vsc") # for ICs -# plt.title(f"Cross correlation N={Nres} L={Lbox:.0f} a=0.02 vsc") # for ICs +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_{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_{Nres}_{Lbox:.0f}_a{scale_factor}") +# plt.title(f"Cross correlation N={Nres} L={Lbox:.0f} a={a[scale_factor]}") plt.xscale("log") -plt.xlabel("k [Mpc]") +plt.xlabel(r"k [$\mathrm{Mpc}^{-1}$]") plt.ylabel("C = Pcross") plt.ylim(0.8, 1.0) plt.xlim(k[0], knyquist) diff --git a/prepare_spectral_evaluation_agora.py b/prepare_spectral_evaluation_agora.py index bd40582..d478b7d 100644 --- a/prepare_spectral_evaluation_agora.py +++ b/prepare_spectral_evaluation_agora.py @@ -11,12 +11,12 @@ from pathlib import Path from write_spectra_jobs import write_spectra_jobs_agora from directories import monofonic_tests_basedir, agora_test_basedir -def main(scale_factor: int, Nres1: int, Nres2: int): - a = [0.166666, 0.333333, 0.5, 0.666666, 1.0] +def main(scale_factor: int, Nres1: int, Nres2: int, Lbox: float): + a = [0.01, 0.02, 0.04, 0.08, 0.166666, 0.333333, 0.5, 0.666666, 1.0] print(f"Chose scale factor a = {a[scale_factor]} (index {scale_factor} / {len(a) - 1})") - write_spectra_jobs_agora(scale_factor, Nres1, Nres2, - output_basedir=agora_test_basedir / "spectra/") + write_spectra_jobs_agora(scale_factor, Nres1, Nres2, Lbox, + output_basedir=agora_test_basedir / f"spectra/{Lbox}mpc") @@ -25,10 +25,12 @@ if __name__ == "__main__": scale_factor = int(argv[1]) Nres1 = int(argv[2]) Nres2 = int(argv[3]) + Lbox = float(argv[4]) main( scale_factor=scale_factor, Nres1=Nres1, - Nres2=Nres2 + Nres2=Nres2, + Lbox=Lbox ) diff --git a/write_spectra_jobs.py b/write_spectra_jobs.py index 4b3f6c2..75284f0 100644 --- a/write_spectra_jobs.py +++ b/write_spectra_jobs.py @@ -42,9 +42,9 @@ set -u -def write_spectra_jobs_agora(scale_factor: int, Nres1: int, Nres2: int, output_basedir: Path): +def write_spectra_jobs_agora(scale_factor: int, Nres1: int, Nres2: int, Lbox: float, output_basedir: Path): Nres_max = max(Nres1, Nres2) - input_dir = str(agora_test_basedir) + input_dir = str(agora_test_basedir / f'{Lbox}mpc/') output_dir = output_basedir output_dir_string = str(output_dir) if output_dir.exists():