maubin/sbmy_control
0
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

new: crop fields for power spectra

Browse source
This commit is contained in:
Mayeul Aubin 2025-04-30 11:45:59 +02:00
parent f63a20bf5b
commit faeed99b09
1 changed files with 42 additions and 7 deletions
Download patch file Download diff file Expand all files Collapse all files

49
analysis/power_spectrum.py
View file

@ -5,6 +5,24 @@ kmax = 2e0
Nk = 50
AliasingCorr=False
def crop_field(field, Ncrop):
if Ncrop is None or Ncrop == 0:
return
elif Ncrop > 0:
field.data = field.data[Ncrop:-Ncrop, Ncrop:-Ncrop, Ncrop:-Ncrop]
d0 = field.L0/field.N0
d1 = field.L1/field.N1
d2 = field.L2/field.N2
field.N0 -= 2*Ncrop
field.N1 -= 2*Ncrop
field.N2 -= 2*Ncrop
field.L0 = field.N0*d0
field.L1 = field.N1*d1
field.L2 = field.N2*d2
def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None):
from pysbmy.power import PowerSpectrum
from pysbmy.fft import FourierGrid
@ -91,7 +109,8 @@ def plot_power_spectra(filenames,
figsize=(8,4),
dpi=300,
ax=None,
fig=None,):
fig=None,
Ncrop=None,):
import matplotlib.pyplot as plt
from pysbmy.field import read_field
@ -110,6 +129,7 @@ def plot_power_spectra(filenames,
for i, filename in enumerate(filenames):
field = read_field(filename)
crop_field(field, Ncrop)
_, G, k, _ = add_power_spectrum_to_plot(ax=ax,
field=field,
Pk_ref=Pk_ref,
@ -163,7 +183,9 @@ def plot_cross_correlations(filenames_A,
figsize=(8,4),
dpi=300,
ax=None,
fig=None,):
fig=None,
Ncrop=None,
):
import matplotlib.pyplot as plt
from pysbmy.field import read_field
@ -181,9 +203,11 @@ def plot_cross_correlations(filenames_A,
markers = [None for f in filenames_A]
field_B = read_field(filename_B)
crop_field(field_B, Ncrop)
for i, filename_A in enumerate(filenames_A):
field_A = read_field(filename_A)
crop_field(field_A, Ncrop)
_, G, k, _ = add_cross_correlations_to_plot(ax=ax,
field_A=field_A,
field_B=field_B,
@ -255,6 +279,7 @@ if __name__ == "__main__":
parser.add_argument('-t','--title', type=str, default=None, help='Title of the plot.')
parser.add_argument('-yrp', '--ylim_power', type=float, nargs=2, default=[0.9,1.1], help='Y-axis limits.')
parser.add_argument('-yrc', '--ylim_corr', type=float, nargs=2, default=[0.99,1.001], help='Y-axis limits.')
parser.add_argument('--crop', type=int, default=None, help='Remove the outter N pixels of the fields.')
args = parser.parse_args()
@ -264,7 +289,9 @@ if __name__ == "__main__":
if args.reference is not None:
from pysbmy.field import read_field
G, _, Pk_ref = get_power_spectrum(read_field(args.directory+args.reference), kmin=kmin, kmax=kmax, Nk=Nk)
F_ref = read_field(args.directory+args.reference)
crop_field(F_ref, args.crop)
G, _, Pk_ref = get_power_spectrum(F_ref, kmin=kmin, kmax=kmax, Nk=Nk)
else:
Pk_ref = None
G = None
@ -294,7 +321,9 @@ if __name__ == "__main__":
kmax=kmax,
Nk=Nk,
ax=axes[0],
fig=fig)
fig=fig,
Ncrop=args.crop,
)
plot_cross_correlations(filenames_A=filenames,
filename_B=args.directory+args.reference,
@ -311,7 +340,9 @@ if __name__ == "__main__":
kmax=kmax,
Nk=Nk,
ax=axes[1],
fig=fig)
fig=fig,
Ncrop=args.crop,
)
axes[1].legend(loc='lower left')
axes[0].set_title("Power Spectrum")
@ -334,7 +365,9 @@ if __name__ == "__main__":
bound2=0.02,
kmin=kmin,
kmax=kmax,
Nk=Nk)
Nk=Nk,
Ncrop=args.crop,
)
ax.legend()
if args.title is not None:
ax.set_title(args.title)
@ -353,7 +386,9 @@ if __name__ == "__main__":
bound2=0.002,
kmin=kmin,
kmax=kmax,
Nk=Nk)
Nk=Nk,
Ncrop=args.crop,
)
ax.legend(loc='lower left')
if args.title is not None:
ax.set_title(args.title)