attoworld.data

13def read_dwc(file_path):
14    """
15    Reads files in the .dwc format produced by many FROG scanners
16
17    Args:
18        file_path: path to the .dwc file
19
20    Returns:
21        Spectrogram: the loaded data
22    """
23    with open(file_path, "r") as file:
24        lines = file.readlines()
25        delay_increment = float(lines[2].strip().split("=")[1])
26
27    wavelength_vector = pd.read_csv(
28        file_path, skiprows=5, nrows=1, delimiter="\t", header=None
29    ).values[0]
30
31    data_array = pd.read_csv(
32        file_path, skiprows=8, delimiter="\t", header=None, dtype=float
33    ).values
34
35    freqs, first_spec = spectrum.wavelength_to_frequency(
36        wavelength_vector, data_array[:, 0]
37    )
38    if freqs is not None:
39        data_array_freq = np.zeros((len(freqs), data_array.shape[1]))
40        data_array_freq[:, 0] = first_spec
41        for _i in range(data_array.shape[1]):
42            _f, _dat = spectrum.wavelength_to_frequency(
43                wavelength_vector, data_array[:, _i], freqs
44            )
45            data_array_freq[:, _i] = _dat
46        dt = 1e-15 * delay_increment
47        delays = dt * np.array(range(data_array.shape[1]))
48        delays -= np.mean(delays)
49        return data_structures.Spectrogram(
50            data=data_array_freq, time=delays, freq=freqs
51        )
52    else:
53        raise Exception("Interpolation failure reading dwc file")

56def load_mean_spectrum_from_scarab(filename: str):
57    data = np.loadtxt(filename)
58    return data_structures.IntensitySpectrum(
59        spectrum=np.mean(data[:, 1::], axis=1),
60        wavelength=1e-9 * data[:, 0],
61        freq=1e9 * constants.speed_of_light / data[:, 0],
62        is_frequency_scaled=False,
63        phase=None,
64    )

128def load_waves_from_matfile(filename: str, phase: Optional[float] = None):
129    """Load the contents of an attolab scanner file in .mat format
130
131    Args:
132        phase (float): phase to use when interpreting the lock-in data
133        filename (str): path to the mat file
134    Returns:
135        time_delay: array of time delay values
136        signal: signals corresponding to the time delays
137    """
138
139    datablob = sio.loadmat(filename)
140    stage_position = datablob["xdata"][0, :]
141    time_delay = -2e-3 * stage_position / 2.9979e8
142    lia_x = datablob["x0"]
143    lia_y = datablob["y0"]
144    if phase is None:
145        optimized_phase = np.atan2(np.sum(lia_y[:] ** 2), np.sum(lia_x[:] ** 2))
146        signal = np.fliplr(
147            lia_x * np.cos(optimized_phase) + lia_y * np.sin(optimized_phase)
148        )
149        return time_delay, signal
150    else:
151        signal = np.fliplr(lia_x * np.cos(phase) - lia_y * np.sin(phase))
152        return time_delay, signal

 97def load_waveform_from_text(
 98    filename: str,
 99    time_multiplier: float = 1e-15,
100    time_field: str = "delay (fs)",
101    wave_field: str = "field (a.u.)",
102    sep="\t",
103) -> data_structures.Waveform:
104    """Loads a waveform from a text file
105
106    Args:
107        filename (str): path to the file
108        time_multiplier (float): multiplier needed to convert the time unit in the file to seconds
109        time_field (str): name (header) of the column containing the times
110        wave_field (str): name (header) of the column containing the waveform
111        sep (str): separator used in the file format (tab is default)
112
113    Returns:
114        Waveform: the waveform
115    """
116
117    data = pd.read_csv(filename, sep=sep)
118    time = time_multiplier * data[time_field].to_numpy()
119    wave = data[wave_field].to_numpy()
120    dt = time[1] - time[0]
121    diff_time = np.diff(time)
122    uniform = bool(np.all(np.isclose(diff_time, diff_time[0])))
123    return data_structures.Waveform(
124        wave=wave, time=time, dt=dt, is_uniformly_spaced=uniform
125    )

67def load_spectrum_from_text(
68    filename: str,
69    wavelength_multiplier: float = 1e-9,
70    wavelength_field: str = "wavelength (nm)",
71    spectrum_field: str = "intensity (a.u.)",
72    sep: str = "\t",
73):
74    """
75    Load a spectrum contained in a text file.
76
77    Args:
78        filename (str): path to the file
79        wavelength_multiplier (float): multiplier to convert wavelength to m
80        wavelength_field (str): name of the field in the data corresponding to wavelength
81        spectrum_field (str): name of the field in the data corresponding to spectral intensity
82        sep (str): column separator
83    Returns:
84        IntensitySpectrum: the intensity spectrum"""
85    data = pd.read_csv(filename, sep=sep)
86    wavelength = wavelength_multiplier * np.array(data[wavelength_field])
87    freq = constants.speed_of_light / wavelength
88    spectrum = np.array(data[spectrum_field])
89    return data_structures.IntensitySpectrum(
90        spectrum=spectrum,
91        wavelength=wavelength,
92        freq=freq,
93        phase=np.zeros(spectrum.shape, dtype=float),
94    )

1115@yaml_io
1116@dataclass(slots=True)
1117class FrogData:
1118    """
1119    Stores data from a FROG measurement
1120
1121    Attributes:
1122        spectrum (ComplexSpectrum): the reconstructed complex spectrum
1123        pulse (Waveform): time-domain reconstructed field
1124        measured_spectrogram (Spectrogram): measured (binned) data
1125        reconstructed_spectrogram (Spectrogram): spectrogram resulting from reconstructed field
1126        f0: the central frequency of the spectrum
1127        dt: the time step
1128    """
1129
1130    spectrum: ComplexSpectrum
1131    pulse: Waveform
1132    measured_spectrogram: Spectrogram
1133    reconstructed_spectrogram: Spectrogram
1134    raw_reconstruction: np.ndarray
1135    f0: float
1136    dt: float
1137
1138    def lock(self):
1139        """
1140        Make the data immutable
1141        """
1142        self.raw_reconstruction.setflags(write=False)
1143        self.spectrum.lock()
1144        self.pulse.lock()
1145        self.measured_spectrogram.lock()
1146        self.reconstructed_spectrogram.lock()
1147
1148    def save(self, base_filename):
1149        """
1150        Save in the Trebino FROG format
1151
1152        Args:
1153            base_filename: base of the file path; 4 files will be made from it: .A.dat, .Arecon.dat, .Ek.dat, and .Speck.dat
1154        """
1155        self.measured_spectrogram.save(base_filename + ".A.dat")
1156        self.reconstructed_spectrogram.save(base_filename + ".Arecon.dat")
1157
1158        t = 1e15 * self.dt * np.array(range(len(self.raw_reconstruction)))
1159        t -= np.mean(t)
1160        f = np.fft.fftshift(np.fft.fftfreq(len(t), d=self.dt)) + self.f0
1161        lam = constants.speed_of_light / f
1162        raw_spec = np.fft.fftshift(np.fft.fft(self.raw_reconstruction))
1163
1164        with open(base_filename + ".Ek.dat", "w") as time_file:
1165            for _i in range(len(self.raw_reconstruction)):
1166                time_file.write(
1167                    f"{t[_i]:.15g}\t{np.abs(self.raw_reconstruction[_i]) ** 2:.15g}\t{np.angle(self.raw_reconstruction[_i]):.15g}\t{np.real(self.raw_reconstruction[_i]):.15g}\t{np.imag(self.raw_reconstruction[_i]):.15g}\n"
1168                )
1169
1170        with open(base_filename + ".Speck.dat", "w") as spec_file:
1171            for _i in range(len(self.raw_reconstruction)):
1172                spec_file.write(
1173                    f"{lam[_i]:.15g}\t{np.abs(raw_spec[_i]) ** 2:.15g}\t{np.angle(raw_spec[_i]):.15g}\t{np.real(raw_spec[_i]):.15g}\t{np.imag(raw_spec[_i]):.15g}\n"
1174                )
1175
1176    def plot_measured_spectrogram(self, ax: Optional[Axes] = None, log: bool = False):
1177        """
1178        Plot the measured spectrogram.
1179
1180        Args:
1181            ax: optionally plot onto a pre-existing matplotlib Axes
1182        """
1183        if ax is None:
1184            fig, ax = plt.subplots()
1185        else:
1186            fig = ax.get_figure()
1187        if log:
1188            self.measured_spectrogram.plot_log(ax)
1189        else:
1190            self.measured_spectrogram.plot(ax)
1191        ax.set_title("Measured")
1192        return fig
1193
1194    def plot_reconstructed_spectrogram(
1195        self, ax: Optional[Axes] = None, log: bool = False
1196    ):
1197        """
1198        Plot the reconstructed spectrogram.
1199
1200        Args:
1201            ax: optionally plot onto a pre-existing matplotlib Axes
1202        """
1203        if ax is None:
1204            fig, ax = plt.subplots()
1205        else:
1206            fig = ax.get_figure()
1207        if log:
1208            self.reconstructed_spectrogram.plot_log(ax)
1209        else:
1210            self.reconstructed_spectrogram.plot(ax)
1211        ax.set_title(
1212            f"Retrieved (G': {self.get_error():0.2e}; G: {self.get_G_error():0.2e})"
1213        )
1214        return fig
1215
1216    def plot_pulse(
1217        self, ax: Optional[Axes] = None, phase_blanking: float = 0.05, xlim=None
1218    ):
1219        """
1220        Plot the reconstructed pulse.
1221
1222        Args:
1223            ax: optionally plot onto a pre-existing matplotlib Axes
1224            phase_blanking: only show phase information (instantaneous frequency) above this level relative to max intensity
1225            xlim: pass arguments to set_xlim() to constrain the x-axis
1226        """
1227        return self.pulse.to_complex_envelope().plot(ax, phase_blanking, xlim)
1228
1229    def plot_spectrum(
1230        self, ax: Optional[Axes] = None, phase_blanking: float = 0.05, xlim=None
1231    ):
1232        """
1233        Plot the reconstructed spectrum and group delay curve.
1234
1235        Args:
1236            ax: optionally plot onto a pre-existing matplotlib Axes
1237            phase_blanking: only show phase information (group delay) above this level relative to max intensity
1238            xlim: pass arguments to set_xlim() to constrain the x-axis
1239        """
1240        return self.spectrum.to_intensity_spectrum().plot_with_group_delay(
1241            ax, phase_blanking=phase_blanking, shift_from_centered=True, xlim=xlim
1242        )
1243
1244    def plot_all(
1245        self,
1246        phase_blanking=0.05,
1247        time_xlims=None,
1248        wavelength_xlims=None,
1249        figsize=None,
1250        log: bool = False,
1251    ):
1252        """
1253        Produce a 4-panel plot of the FROG results, combining calls to plot_measured_spectrogram(),
1254        plot_reconstructed_spectrogram(), plot_pulse() and plot_spectrum() as subplots, with letter labels.
1255
1256        Args:
1257            phase_blanking: relative intensity at which to show phase information
1258            time_xlim: x-axis limits to pass to the plot of the pulse
1259            wavelength_xlim: x-axis limits to pass to the plot of the spectrum
1260            figsize: custom figure size
1261        """
1262        if figsize is None:
1263            default_figsize = plt.rcParams["figure.figsize"]
1264            figsize = (default_figsize[0] * 2, default_figsize[1] * 2)
1265        fig, ax = plt.subplots(2, 2, figsize=figsize)
1266        self.plot_measured_spectrogram(ax[0, 0], log=log)
1267        label_letter("a", ax[0, 0])
1268        self.plot_reconstructed_spectrogram(ax[1, 0], log=log)
1269        label_letter("b", ax[1, 0])
1270        self.plot_pulse(ax[0, 1], xlim=time_xlims, phase_blanking=phase_blanking)
1271        label_letter("c", ax[0, 1])
1272        self.plot_spectrum(
1273            ax[1, 1], xlim=wavelength_xlims, phase_blanking=phase_blanking
1274        )
1275        label_letter("d", ax[1, 1])
1276        return fig
1277
1278    def get_error(self) -> float:
1279        """
1280        Get the G' error of the reconstruction
1281        """
1282        norm_measured = np.linalg.norm(self.measured_spectrogram.data)
1283        norm_retrieved = np.linalg.norm(self.reconstructed_spectrogram.data)
1284        return np.sqrt(
1285            np.sum(
1286                (
1287                    self.measured_spectrogram.data[:] / norm_measured
1288                    - self.reconstructed_spectrogram.data[:] / norm_retrieved
1289                )
1290                ** 2
1291            )
1292            / np.sum((self.measured_spectrogram.data[:] / norm_measured) ** 2)
1293        )
1294
1295    def get_G_error(self) -> float:
1296        """
1297        Get the G (note: no apostrophe) error. This one doesn't mean much, but is useful
1298        for comparing reconstructions of the same spectrogram between different programs.
1299        """
1300        return np.sqrt(
1301            (1.0 / float(len(self.measured_spectrogram.data[:]) ** 2))
1302            * np.sum(
1303                (
1304                    self.measured_spectrogram.data[:]
1305                    / np.max(self.measured_spectrogram.data[:])
1306                    - self.reconstructed_spectrogram.data[:]
1307                    / np.max(self.reconstructed_spectrogram.data[:])
1308                )
1309                ** 2
1310            )
1311        )
1312
1313    def get_fwhm(self) -> float:
1314        """
1315        Get the full-width-at-half-max value of the reconstructed pulse
1316        """
1317        return self.pulse.get_envelope_fwhm()

155def read_Trebino_FROG_matrix(filename: Path | str) -> data_structures.Spectrogram:
156    """
157    Read a spectrogram file made by the Trebino FROG code
158
159    Args:
160        filename (Path | str): the name (path) of the file
161    """
162    with open(filename, "r") as f:
163        line = str(f.readline())
164        line = line.split()
165        n1 = int(line[0])
166        n2 = int(line[1])
167        line = str(f.readline())
168        line = line.split()
169    measured_data = pd.read_csv(filename, sep="\t", header=None, skiprows=2)
170    measure = []
171    raw_freq = (
172        1e9 * constants.speed_of_light / np.array(measured_data[0][0:n2]).squeeze()
173    )
174    df = np.mean(np.diff(raw_freq))
175    freq = raw_freq[0] + df * np.array(range(raw_freq.shape[0]))
176    time = 1e-15 * np.array(measured_data[0][n2 : (n2 + n1)]).squeeze()
177    for i in range(n1):
178        measure.append(measured_data[0][(i + 2) * n2 : (i + 3) * n2])
179    data = np.array(measure)
180    return data_structures.Spectrogram(data=data, time=time, freq=freq)

183def read_Trebino_FROG_speck(filename: Path | str) -> data_structures.ComplexSpectrum:
184    """
185    Read a .Speck file made by the Trebino FROG code
186
187    Args:
188        filename (Path | str): the name (path) of the file
189    """
190    data = np.array(pd.read_csv(filename, sep="\t", header=None), dtype=float)
191    raw_freq = 1e9 * constants.speed_of_light / data[:, 0]
192    df = np.mean(np.diff(raw_freq))
193    freq = np.linspace(0.0, raw_freq[-1], int(np.ceil(raw_freq[-1] / df)))
194    spectrum = interpolate(freq, raw_freq, data[:, 3]) + 1j * interpolate(
195        freq, raw_freq, data[:, 4]
196    )
197    return data_structures.ComplexSpectrum(spectrum=spectrum, freq=freq)

200def read_Trebino_FROG_data(filename: str) -> data_structures.FrogData:
201    """
202    Read a set of data produced by the Trebino FROG reconstruction code.
203
204    Args:
205        filename: Base filename of the .bin file; e.g. if the data are mydata.bin.Speck.dat etc., this will be "mydata.bin" """
206    spectrum = read_Trebino_FROG_speck(filename + ".Speck.dat")
207    pulse = spectrum.to_centered_waveform()
208    measured_spectrogram = read_Trebino_FROG_matrix(filename + ".A.dat")
209    reconstructed_spectrogram = read_Trebino_FROG_matrix(filename + ".Arecon.dat")
210    raw_speck = np.array(
211        pd.read_csv(filename + ".Speck.dat", sep="\t", header=None), dtype=float
212    )
213    raw_ek = np.array(
214        pd.read_csv(filename + ".Ek.dat", sep="\t", header=None), dtype=float
215    )
216    f0 = 1e9 * np.mean(constants.speed_of_light / raw_speck[:, 0])
217    dt = 1e-15 * (raw_ek[1, 0] - raw_ek[0, 0])
218    raw_reconstruction = raw_ek[:, 3] + 1.0j * raw_ek[:, 4]
219    return data_structures.FrogData(
220        spectrum=spectrum,
221        pulse=pulse,
222        measured_spectrogram=measured_spectrogram,
223        reconstructed_spectrogram=reconstructed_spectrogram,
224        raw_reconstruction=raw_reconstruction,
225        dt=dt,
226        f0=float(f0),
227    )

 39class LunaResult:
 40    """Loads and handles the Luna simulation result.
 41
 42    The result must be in the HDF5 format using the saving option in the Luna.Interface.prop_capillary() function [filepath="..."].
 43    As opposed to most of the analysis routines here, units are in SI!
 44
 45    Attributes:
 46        z: position along the fiber (for the field data)
 47        fieldFT: complex FFT of the field (positive freq axis); shape: (n_z, n_modes, n_freq) [or (n_z, n_freq) after mode selection/averaging]
 48        omega: angular frequency axis for the FFT field; shape: (n_z, n_modes, n_freq) [or (n_z, n_freq) after mode selection/averaging]
 49
 50        stats_z: position along the fiber (for stats data)
 51        stats_density: particle density along the fiber (m^-3)
 52        stats_electrondensity: electron density along the fiber (m^-3)
 53        stats_pressure: gas pressure profile along the fiber (bar)
 54        stats_energy: pulse energy along the fiber (J)
 55        stats_peakpower: peak power of the pulse along the fiber (W)
 56        stats_peakintensity: peak intensity of the pulse along the fiber
 57        stats_peak_ionization_rate: peak ionization rate along the fiber
 58    """
 59
 60    def __init__(self, filename):
 61        """Constructor of class LunaResult.
 62
 63        Args:
 64            filename: saved result, file path
 65        """
 66
 67        self.filename = filename
 68        self.fieldFT = None
 69        self.omega = None
 70        self.z = None
 71
 72        self.stats_z = None
 73        self.stats_energy = None
 74        self.stats_electrondensity = None
 75        self.stats_density = None
 76        self.stats_pressure = None
 77        self.stats_peakpower = None
 78        self.stats_peakintensity = None
 79        self.stats_zdw = None
 80        self.stats_peak_ionization_rate = None
 81
 82        self.open_Luna_result(filename)
 83
 84    def open_Luna_result(self, filename):
 85        """Opens the Luna result file and loads the data"""
 86        with h5py.File(filename, "r") as data:
 87            # FIELD DATA
 88            self.fieldFT = np.array(data["Eω"])
 89            self.omega = np.array(data["grid"]["ω"])
 90            self.z = np.array(data["z"])
 91
 92            # STATS
 93            self.stats_z = np.array(data["stats"]["z"])
 94            self.stats_energy = np.array(data["stats"]["energy"])
 95            self.stats_electrondensity = np.array(data["stats"]["electrondensity"])
 96            self.stats_density = np.array(data["stats"]["density"])
 97            self.stats_pressure = np.array(data["stats"]["pressure"])
 98            self.stats_peakpower = np.array(data["stats"]["peakpower"])
 99            self.stats_peakintensity = np.array(data["stats"]["peakintensity"])
100            self.stats_zdw = np.array(data["stats"]["zdw"])
101            self.stats_peak_ionization_rate = np.array(
102                data["stats"]["peak_ionisation_rate"]
103            )
104
105    def average_modes(self):
106        """Averages the propagation modes in the Luna result file"""
107        if len(self.fieldFT.shape) == 3:
108            self.fieldFT = np.mean(self.fieldFT, axis=1)
109            self.stats_zdw = None
110            self.stats_peakpower = None
111            self.stats_energy = np.sum(self.stats_energy, axis=1)
112
113    def select_mode(self, mode: int):
114        if len(self.fieldFT.shape) < 3:
115            print("WARNING: No mode to select")
116        elif mode >= self.fieldFT.shape[1] or mode < 0:
117            print("WARNING: mode ", mode, " is out of range")
118        else:
119            self.fieldFT = self.fieldFT[:, mode, :]
120            self.stats_zdw = self.stats_zdw[:, mode]
121            self.stats_peakpower = self.stats_peakpower[:, mode]
122            self.stats_energy = self.stats_energy[:, mode]
123
124    def get_time_field(self, position=None):
125        """Get the electric field in time from the Luna result file. If no mode was previously selected, the method computes the average of all modes.
126        Therefore, after calling get_time_field(), mode selection is not possible any more.
127
128        Args:
129            position (float): position along the fiber in m. If None, the end of the fiber is used.
130
131        Returns:
132            timeV (numpy.ndarray): time axis in seconds
133            fieldV (numpy.ndarray): electric field in V/m
134        """
135        self.average_modes()
136        if position is None:
137            position = self.z[-1]
138        index = np.argmin(np.abs(self.z - position))
139        if position > np.max(self.z) or position < np.min(self.z):
140            print("WARNING: position ", position, "m is out of range")
141        check_equal_length(self.fieldFT[index], self.omega)
142        fieldFFT = np.concatenate(
143            (self.fieldFT[index, :], np.conjugate(self.fieldFT[index, :][::-1]) * 0)
144        )
145        freq = np.concatenate((self.omega, -self.omega[::-1])) / 2 / np.pi
146        timeV, fieldV = inverse_fourier_transform(freq, fieldFFT)
147        return timeV, np.real(fieldV)
148
149    def get_wavelength_spectrum(self, position=None):
150        """Get the spectrum from the Luna result file (|FFT|^2 * (2 * pi * c / λ^2))
151
152        Args:
153            position (float): position along the fiber in m. If None, the end of the fiber is used.
154
155        Returns:
156            wvl (numpy.ndarray): wavelength axis in m
157            wvlSpectrum (numpy.ndarray): electric field spectrum in V/m
158        """
159        self.average_modes()
160        if position is None:
161            position = self.z[-1]
162        index = np.argmin(np.abs(self.z - position))
163        if position > np.max(self.z) or position < np.min(self.z):
164            print("WARNING: position ", position, "m is out of range")
165        wvl = 2 * np.pi * constants.speed_of_light / self.omega[::-1]
166        wvlSpectrum = np.abs(self.fieldFT[index, ::-1]) ** 2 * (
167            2 * np.pi * constants.speed_of_light / wvl**2
168        )
169        return wvl, wvlSpectrum
170
171    def get_spectral_phase(self, position=None):
172        """Get the spectral phase from the Luna result file
173
174        Args:
175            position (float): position along the fiber in m. If None, the end of the fiber is used.
176
177        Returns:
178            wvl (numpy.ndarray): wavelength axis in m
179            phase (numpy.ndarray): spectral phase in rad
180        """
181        self.average_modes()
182        if position is None:
183            position = self.z[-1]
184        index = np.argmin(np.abs(self.z - position))
185        if position > np.max(self.z) or position < np.min(self.z):
186            print("WARNING: position ", position, "m is out of range")
187        wvl = 2 * np.pi * constants.speed_of_light / self.omega[::-1]
188        phase = np.angle(self.fieldFT[index, ::-1])
189        return wvl, phase
190
191    def plot_stats(self):
192        """Plots the 'stats_' attribute of the simulation stored in the present object."""
193
194        fig, axs = plt.subplots(2, 2, figsize=[7, 5])
195        axs[0, 0].plot(self.stats_z, self.stats_pressure)
196        axs[0, 0].set_xlabel("z (m)")
197        axs[0, 0].set_ylabel("gas pressure (bar)")
198        ax2 = axs[0, 0].twinx()
199        ax2.plot(self.stats_z, self.stats_density, color="r")
200        ax2.set_ylabel("gas particle density ($m^{-3}$)", color="r")
201        ax2.tick_params(axis="y", colors="r")
202        ax2.yaxis.label.set_color("r")
203        axs[0, 1].plot(self.stats_z, self.stats_electrondensity)
204        axs[0, 1].set_xlabel("z (m)")
205        axs[0, 1].set_ylabel("electron density ($m^{-3}$)")
206        ax2 = axs[0, 1].twinx()
207        if len(self.stats_energy.shape) == 2:
208            ax2.plot(self.stats_z, np.sum(self.stats_energy, axis=1) * 1e6, color="r")
209        else:
210            ax2.plot(self.stats_z, self.stats_energy * 1e6, color="r")
211        ax2.set_ylabel(f"pulse energy ({Char.mu}J)", color="r")
212        ax2.tick_params(axis="y", colors="r")
213        ax2.yaxis.label.set_color("r")
214        if self.stats_peakpower is not None:
215            if len(self.stats_peakpower.shape) == 2:
216                axs[1, 0].plot(self.stats_z, self.stats_peakpower[:, 0])
217                axs[1, 0].set_xlabel("z (m)")
218                axs[1, 0].set_ylabel("peak power (W)")
219            elif len(self.stats_peakpower.shape) == 1:
220                axs[1, 0].plot(self.stats_z, self.stats_peakpower)
221                axs[1, 0].set_xlabel("z (m)")
222                axs[1, 0].set_ylabel("peak power (W)")
223        axs[1, 1].plot(self.stats_z, self.stats_peak_ionization_rate)
224        axs[1, 1].set_xlabel("z (m)")
225        axs[1, 1].set_ylabel("peak ionization rate")
226        plt.tight_layout()
227
228        return fig