From 727f12dda2441ef57c8d73a435c75eb8df7b33c3 Mon Sep 17 00:00:00 2001 From: Igor Andriyash Date: Thu, 9 Oct 2025 11:57:00 +0200 Subject: [PATCH 1/3] revised use of phto_num and lambda_0 arguments in utils replaced numpy trapz with trapezoid --- synchrad/utils.py | 31 +++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/synchrad/utils.py b/synchrad/utils.py index 1e71202..def8932 100644 --- a/synchrad/utils.py +++ b/synchrad/utils.py @@ -61,13 +61,6 @@ def get_full_spectrum(self, spect_filter=None, \ if normalize_to_weights: val /= self.total_weight - if phot_num: - ax = self.Args['omega'] - val /= ax[:,None,None] - else: - if lambda0_um is not None: - val *= J_in_um / lambda0_um - return val def get_energy_spectrum(self, spect_filter=None, \ @@ -79,13 +72,13 @@ def get_energy_spectrum(self, spect_filter=None, \ if self.Args['mode'] == 'far': theta_loc = 0.5 * (self.Args['theta'][1:] + self.Args['theta'][:-1]) val_loc = 0.5 * (val[:,1:,:] + val[:,:-1,:]) - int_theta = np.trapz( val_loc * np.sin(theta_loc)[None,:,None], + int_theta = np.trapezoid( val_loc * np.sin(theta_loc)[None,:,None], theta_loc, axis=1) val = self.Args['dph'] * int_theta.sum(-1) elif self.Args['mode'] == 'near': r_loc = self.Args['radius'] - int_r = np.trapz( val * r_loc[None,:,None], r_loc, axis=1) + int_r = np.trapezoid( val * r_loc[None,:,None], r_loc, axis=1) val = self.Args['dph'] * int_r.sum(-1) return val @@ -96,7 +89,13 @@ def get_energy(self, spect_filter=None, \ val = self.get_energy_spectrum(spect_filter=spect_filter, \ phot_num=phot_num, lambda0_um=lambda0_um, **kw_args) - val = np.trapz(val, self.Args['omega']) + if phot_num: + val /= self.Args['omega'] + else: + if lambda0_um is not None: + val *= J_in_um / lambda0_um + + val = np.trapezoid(val, self.Args['omega']) return val def get_spot(self, k0=None, spect_filter=None, \ @@ -106,13 +105,17 @@ def get_spot(self, k0=None, spect_filter=None, \ phot_num=phot_num, lambda0_um=lambda0_um, **kw_args) if k0 is None: + if phot_num: + ax = self.Args['omega'] + val /= ax[:,None,None] + else: + if lambda0_um is not None: + val *= J_in_um / lambda0_um + if val.shape[0]>1: - val = np.trapz(val, self.Args['omega'], axis=0) + val = np.trapezoid(val, self.Args['omega'], axis=0) else: val = val[0] * self.Args['dw'] - #if val.shape[0]>1: - # val = 0.5*(val[1:] + val[:-1]) - #val = (val*self.Args['dw'][:, None, None]).sum(0) else: ax = self.Args['omega'] indx = (ax Date: Thu, 9 Oct 2025 12:46:19 +0200 Subject: [PATCH 2/3] re-added tracksFromOPMD method to utils for compatibility --- synchrad/utils.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/synchrad/utils.py b/synchrad/utils.py index def8932..5bbe435 100644 --- a/synchrad/utils.py +++ b/synchrad/utils.py @@ -5,6 +5,8 @@ from scipy.interpolate import griddata from scipy.ndimage import gaussian_filter +from .converters import tracksFromOPMD + try: from tvtk.api import tvtk, write_data tvtk_installed = True From b7dc6b9bc1239161471501dd468e26605a4f1efa Mon Sep 17 00:00:00 2001 From: Igor Andriyash Date: Thu, 9 Oct 2025 13:08:07 +0200 Subject: [PATCH 3/3] revised arguments in utils in PIC tutorial plotter --- tutorials/PIC/plot_spectrum.py | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/tutorials/PIC/plot_spectrum.py b/tutorials/PIC/plot_spectrum.py index 5fc1170..115dd41 100644 --- a/tutorials/PIC/plot_spectrum.py +++ b/tutorials/PIC/plot_spectrum.py @@ -53,11 +53,6 @@ calc = SynchRad(file_spectrum="spectrum.h5") -# with h5py.File("spectrum.h5", "r") as f: -# calc.Data["radiation"] = {} -# for key in f["radiation"].keys(): -# calc.Data["radiation"][key] = f[f"radiation/{key}"][...] / total_particle_weight - # compute total emitted energy E_cutoff = 1.0 # keV eph_keV_m = 1.24e-9 @@ -68,21 +63,21 @@ spect_filter = (energy_axis_full > E_cutoff)[:, np.newaxis, np.newaxis] energy_tot = calc.get_energy( - lambda0_um=1e6, phot_num=False, spect_filter=spect_filter + lambda0_um=1e6, spect_filter=spect_filter ) energy_per_C = energy_tot/(e*calc.total_weight) print(f"Total energy emitted in >{E_cutoff:g} keV: {energy_per_C:g} J/C") # plot energy spectrum pyplot.figure() - energy_spectrum1D = calc.get_energy_spectrum(lambda0_um=1e6) - pyplot.semilogx(energy_axis, energy_spectrum1D) + energy_spectrum1D = calc.get_energy_spectrum() + pyplot.semilogx(energy_axis, energy_spectrum1D * 1e-3) # 1e-3 is for 0.1% units pyplot.xlabel("Photon energy (keV)") pyplot.ylabel("Brightness (ph./0.1%b.w./e$^-$)") pyplot.savefig("energy_spectrum_1D.png") # plot real-space emission - spotXY_far, ext_far = calc.get_spot_cartesian(bins=(512, 512), lambda0_um=1.0) + spotXY_far, ext_far = calc.get_spot_cartesian(bins=(512, 512), lambda0_um=1e6) fig = pyplot.figure(figsize=(8, 8)) Plot2D(