mirror/spectra_python_files
1
0
Fork 0
mirror of https://github.com/glatterf42/spectra_python_files.git synced 2024-09-19 15:53:44 +02:00
Code Releases Activity

Adapted power spectra evaluation

Browse source
This commit is contained in:
glatterf42 2022-05-12 13:04:28 +02:00
parent ad04f7b888
commit efb271e3d4
3 changed files with 57 additions and 24 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

56
agora_plot_power_spectra.py
View file

@ -6,6 +6,7 @@ Created on Wed Mar 16 15:12:21 2022
@author: ben @author: ben
""" """
from matplotlib import scale
import pandas as pd import pandas as pd
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from pathlib import Path from pathlib import Path
@ -15,21 +16,30 @@ import numpy as np
def Dplus( lambda0, a ): 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 ) 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/") # basedir = Path("/home/ben/sims/gadget4/examples/agora_test/output/spectra/")
#choose waveform and Lbox: #choose waveform and Lbox:
Lbox = 85.47 # 60 #only option as of now Lbox = 85.47 # 60 #only option as of now
Nres1 = 128 Nres1 = 512
# Nres2 = 256 # Nres2 = 256
k0 = 2 * 3.14159265358979323846264338327950 / Lbox k0 = 2 * 3.14159265358979323846264338327950 / Lbox
knyquist1 = Nres1 * k0 knyquist1 = Nres1 * k0
# knyquist2 = Nres2 * 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, 0.04, 0.08, 0.1]
# times = [0.01, 0.02] # times = [0.01, 0.02]
# scale_factor = 0 # give index of a list above # 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_m = 0.272
Omega_L = 0.728 Omega_L = 0.728
lambda_0 = Omega_L / Omega_m 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 zstart = 100 # 99.09174098173423
a_ics = 1 / (1 + zstart) 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' # 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) df_ics = pd.read_csv(f'{filename_ics}.txt', sep=' ', skipinitialspace=True, header=None, names=columns, skiprows=1)
#only consider rows above resolution limit #only consider rows above resolution limit
@ -53,17 +69,26 @@ Dplus0 = Dplus(lambda0=lambda_0, a=a_ics)
D_squared_ics = Dplus0 ** 2 D_squared_ics = Dplus0 ** 2
p1_ics_noramlised = p1_ics / D_squared_ics p1_ics_noramlised = p1_ics / D_squared_ics
print(D_squared_ics)
# p1_ics_ic_normalised = p1_ics_noramlised / p1_ics_noramlised # p1_ics_ic_normalised = p1_ics_noramlised / p1_ics_noramlised
# plt.loglog(k_ics, p1_ics_ic_normalised, label='ICs') # plt.loglog(k_ics, p1_ics_ic_normalised, label='ICs')
for scale_factor in range(len(times)): # #This is for 128 particles with fewer output times:
filename = basedir / f"{Lbox}mpc/agora_a{scale_factor}_cross_spectrum" # 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}_small_dt_cross_spectrum"
# filename = basedir / f'agora_a{scale_factor}_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 # 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 # 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 # 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) 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]"] k = df["k [Mpc]"]
p1 = df["P1"] p1 = df["P1"]
p1_error = df["err. P1"] p1_error = df["err. P1"]
# p2 = df["P2"] # p2 = df["P2"]
# p2_error = df["err. P2"] # p2_error = df["err. P2"]
# pcross = df["Pcross"] # 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_normalised = p1 / D_squared
p1_ic_normalised = p1_normalised / p1_ics_noramlised p1_ic_normalised = p1_normalised / p1_ics_noramlised
print(D_squared)
# Plot the power spectra: # 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.loglog(k, p2, label="P2")
plt.title(f"Power Spectra Agora 128") plt.title(f"Power Spectra Agora {Nres1}")
savedir = Path(f"/home/ben/Pictures/swift/agora/spectra/{Lbox}mpc/") savedir = Path(f"/home/ben/Pictures/swift/agora/spectra/{Lbox}mpc_{Nres1}")
plt.xlabel(r"k [$\mathrm{Mpc}^{-1}$]") plt.xlabel(r"k [$\mathrm{Mpc}^{-1}$]")
plt.ylabel("P") 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.vlines(knyquist2, ymin=min(p2), ymax=max(p2), color="black", linestyles="dashed", label=f"{Nres2}")
plt.legend() plt.legend()
# plt.ylim(1, 3) # plt.ylim(1, 3)
plt.savefig(f"{savedir}/agora_{Nres1}_power.png")
plt.show() plt.show()
# plt.savefig(f"{savedir}_2.png")

17
plot_cross_spectra.py
View file

@ -15,16 +15,17 @@ basedir = Path("/home/ben/sims/data_swift/monofonic_tests/spectra/")
#choose Nres and Lbox: #choose Nres and Lbox:
waveforms = ['DB2', "DB4", "DB8", "shannon"] #DB2, DB4, DB8, shannon are all we have right now waveforms = ['DB2', "DB4", "DB8", "shannon"] #DB2, DB4, DB8, shannon are all we have right now
Lbox = 100.0 #only option as of 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 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] a = [0.166666, 0.333333, 0.5, 0.666666, 1.0]
scale_factor = 4 # give index of a list above scale_factor = 4 # give index of a list above
for wave in waveforms: 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}_a{scale_factor}_{Nres1}_{Nres2}_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}_ics_{Nres1}_{Nres2}_cross_spectrum" # for ICs
#find columns in file manually #find columns in file manually
#is k really in Mpc? Swift doesn't use /h internally at least. #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: # Plot the Cross Correlation:
plt.plot(k, pcross, label=f"{wave}") 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 savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/cross_{Nres1}_{Nres2}_{Lbox:.0f}_ics") # for ICs
plt.title(f"Cross correlation N={Nres} L={Lbox:.0f} a=0.02") # 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}") # savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/cross_{Nres1}_{Nres2}_{Lbox:.0f}_a{scale_factor}")
# plt.title(f"Cross correlation N={Nres} L={Lbox:.0f} a={a[scale_factor]}") # plt.title(f"Cross correlation N=({Nres1}, {Nres2}) L={Lbox:.0f} a={a[scale_factor]}")
plt.xscale("log") plt.xscale("log")
plt.xlabel(r"k [$\mathrm{Mpc}^{-1}$]") plt.xlabel(r"k [$\mathrm{Mpc}^{-1}$]")

8
plot_power_spectra.py
View file

@ -10,20 +10,20 @@ import pandas as pd
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from pathlib import Path 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: #choose waveform and Lbox:
waveform = "shannon" #DB2, DB4, DB8 or shannon waveform = "shannon" #DB2, DB4, DB8 or shannon
Lbox = 100.0 #only option as of now Lbox = 100.0 #only option as of now
Nres1 = 128 Nres1 = 128
Nres2 = 256 Nres2 = 512
k0 = 2 * 3.14159265358979323846264338327950 / Lbox k0 = 2 * 3.14159265358979323846264338327950 / Lbox
knyquist1 = Nres1 * k0 knyquist1 = Nres1 * k0
knyquist2 = Nres2 * k0 knyquist2 = Nres2 * k0
a = [0.166666, 0.333333, 0.5, 0.666666, 1.0] a = [0.166666, 0.333333, 0.5, 0.666666, 1.0]
scale_factor = 4 # give index of a list above 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 # 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 # 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 # 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.loglog(k, p2, label="P2")
plt.title(f"Power Spectra {waveform} L={Lbox:.0f} a={a[scale_factor]}") 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.xlabel("k [Mpc]")
plt.ylabel("P") plt.ylabel("P")