Adapted power spectra plot so that k has correct units.

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")
 
 
 

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)

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
         )
 

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():