diff --git a/docs/source/api/profiles/index.rst b/docs/source/api/profiles/index.rst index a93251fc..3ee152a6 100644 --- a/docs/source/api/profiles/index.rst +++ b/docs/source/api/profiles/index.rst @@ -10,4 +10,4 @@ Laser Profiles combined_profile longitudinal/index transverse/index - speckled + speckled/index diff --git a/docs/source/api/profiles/speckled.rst b/docs/source/api/profiles/speckled.rst deleted file mode 100644 index 45057300..00000000 --- a/docs/source/api/profiles/speckled.rst +++ /dev/null @@ -1,4 +0,0 @@ -Speckled Laser Profile -====================== - -.. autoclass:: lasy.profiles.SpeckleProfile diff --git a/docs/source/api/profiles/speckled/fm_ssd.rst b/docs/source/api/profiles/speckled/fm_ssd.rst new file mode 100644 index 00000000..40d2a5c6 --- /dev/null +++ b/docs/source/api/profiles/speckled/fm_ssd.rst @@ -0,0 +1,5 @@ +Frequency-modulated Smoothing by Spectral Dispersion (FM-SSD) Laser Profile +=========================================================================== + +.. autoclass:: lasy.profiles.speckled.FM_SSD_Profile + :members: evaluate diff --git a/docs/source/api/profiles/speckled/gp_isi.rst b/docs/source/api/profiles/speckled/gp_isi.rst new file mode 100644 index 00000000..972da5d6 --- /dev/null +++ b/docs/source/api/profiles/speckled/gp_isi.rst @@ -0,0 +1,5 @@ +Smoothing by Induced Spatial Incoherence (ISI) Laser Profile +============================================================ + +.. autoclass:: lasy.profiles.speckled.GP_ISI_Profile + :members: evaluate diff --git a/docs/source/api/profiles/speckled/gp_rpm_ssd.rst b/docs/source/api/profiles/speckled/gp_rpm_ssd.rst new file mode 100644 index 00000000..4f7d17d7 --- /dev/null +++ b/docs/source/api/profiles/speckled/gp_rpm_ssd.rst @@ -0,0 +1,5 @@ +Random Phase Modulated (RPM) Smoothing by Spectral Dispersion (RPM-SSD) Laser Profile +===================================================================================== + +.. autoclass:: lasy.profiles.speckled.GP_RPM_SSD_Profile + :members: evaluate diff --git a/docs/source/api/profiles/speckled/index.rst b/docs/source/api/profiles/speckled/index.rst new file mode 100644 index 00000000..7571ed82 --- /dev/null +++ b/docs/source/api/profiles/speckled/index.rst @@ -0,0 +1,11 @@ +Speckled Laser Profiles +======================= + +.. toctree:: + :maxdepth: 4 + :hidden: + + rpp_cpp + fm_ssd + gp_rpm_ssd + gp_isi diff --git a/docs/source/api/profiles/speckled/rpp_cpp.rst b/docs/source/api/profiles/speckled/rpp_cpp.rst new file mode 100644 index 00000000..3f0cb292 --- /dev/null +++ b/docs/source/api/profiles/speckled/rpp_cpp.rst @@ -0,0 +1,5 @@ +RPP/CPP only Laser Profile +========================== + +.. autoclass:: lasy.profiles.speckled.PhasePlateProfile + :members: evaluate diff --git a/lasy/profiles/__init__.py b/lasy/profiles/__init__.py index d2e8ee20..034b1a9f 100644 --- a/lasy/profiles/__init__.py +++ b/lasy/profiles/__init__.py @@ -3,7 +3,6 @@ from .from_array_profile import FromArrayProfile from .from_openpmd_profile import FromOpenPMDProfile from .from_insight_file import FromInsightFile -from .speckle_profile import SpeckleProfile __all__ = [ "CombinedLongitudinalTransverseProfile", @@ -11,5 +10,4 @@ "FromArrayProfile", "FromOpenPMDProfile", "FromInsightFile", - "SpeckleProfile", ] diff --git a/lasy/profiles/speckle_profile.py b/lasy/profiles/speckle_profile.py deleted file mode 100644 index 5f808ea0..00000000 --- a/lasy/profiles/speckle_profile.py +++ /dev/null @@ -1,540 +0,0 @@ -import numpy as np -from scipy.constants import c - -from .profile import Profile - - -def gen_gaussian_time_series(t_num, dt, fwhm, rms_mean): - """Generate a discrete time series that has gaussian power spectrum. - - Parameters - ---------- - t_num: number of grid points in time - fwhm: full width half maximum of the power spectrum - rms_mean: root-mean-square average of the spectrum - - Returns - ------- - temporal_amplitude: a time series array of complex numbers with shape [t_num] - """ - if fwhm == 0.0: - temporal_amplitude = np.zeros(t_num, dtype=np.complex128) - else: - omega = np.fft.fftshift(np.fft.fftfreq(t_num, d=dt)) - psd = np.exp(-np.log(2) * 0.5 * np.square(omega / fwhm * 2 * np.pi)) - spectral_amplitude = np.array(psd) * ( - np.random.normal(size=t_num) + 1j * np.random.normal(size=t_num) - ) - temporal_amplitude = np.fft.ifftshift( - np.fft.fft(np.fft.fftshift(spectral_amplitude)) - ) - temporal_amplitude *= rms_mean / np.sqrt( - np.mean(np.square(np.abs(temporal_amplitude))) - ) - return temporal_amplitude - - -class SpeckleProfile(Profile): - r""" - Class for the profile of a speckled laser pulse. - - Speckled lasers are used to mitigate laser-plasma interactions in fusion and ion acceleration contexts. - More on the subject can be found in chapter 9 of `P. Michel, Introduction to Laser-Plasma Interactions `__. - A speckled laser beam is a laser that is deliberately divided transversely into :math:`N_{bx}\times N_{by}` beamlets in the near-field. - The phase plate provides a different phase to each beamlet, with index :math:`ml`, which then propagate to the far field and combine incoherently. - - The electric field in the focal plane, as a function of time :math:`t` and the coordinates - :math:`\boldsymbol{x}_\perp=(x,y)` transverse to the direction of propagation, is: - - .. math:: - - \begin{aligned} - E_u(\boldsymbol{x}_\perp,t) &= Re\left[ E_0 - {\rm sinc}\left(\frac{\pi x}{\Delta x}\right) - {\rm sinc}\left(\frac{\pi y}{\Delta y}\right)\times p_u - \right. - \\ - & \times\sum_{m,l=1}^{N_{bx}, N_{by}} A_{ml} - \exp\left(i\boldsymbol{k}_{\perp ml}\cdot\boldsymbol{x}_\perp - + i\phi_{{\rm RPP/CPP},ml}+i\psi_{{\rm SSD},ml}(t)\right) - \Bigg] - \end{aligned} - - where :math:`u` is either :math:`x` or :math:`y`, :math:`p_u` is - the polarization vector, and :math:`Re` represent the real part [Michel, Eqns. 9.11, 87, 94]. - Several quantities are computed internally to the code depending on the - method of smoothing chosen, including the beamlet amplitude :math:`A_{ml}`, - the beamlet wavenumber :math:`k_{\perp ml}`, - the relative phase contribution :math:`\phi_{{\rm RPP/CPP},ml}` of beamlet :math:`ml` induced by the phase plate, - and the phase contribution to beamlet :math:`ml`, :math:`\psi_{{\rm SSD},ml}(t)`, from the temporal smoothing. - The beam widths are :math:`\Delta x=\frac{\lambda_0fN_{bx}}{D_{x}}`, - :math:`\Delta y=\frac{\lambda_0fN_{by}}{D_{y}}`. - The other parameters in these formulas are defined below. - - This profile admits several options for calculating the amplitudes and phases of the beamlets: - - * Random phase plates (RPP), ``'RPP'``: - Here the phase plate contribution to the phase of beamlet :math:`ml` is :math:`\phi_{{\rm RPP},ml}\in\{0,\pi\}` (drawn randomly for each :math:`m,l`), :math:`\psi_{{\rm SSD},ml}(t)=0`, and :math:`A_{ml}=1` - * Continuous phase plates (CPP), ``'CPP'``: - :math:`\phi_{{\rm CPP},ml}\in[0,\pi]` (drawn randomly with uniform probability between :math:`0` and :math:`\pi`, for each :math:`m,l`), :math:`\psi_{{\rm SSD},ml}(t)=0`, and :math:`A_{ml}=1` - * CPP + Smoothing by spectral dispersion (SSD), ``'FM SSD'``: - :math:`\phi_{{\rm CPP},ml}\in[0,\pi]` (drawn randomly with uniform probability between :math:`0` and :math:`\pi`, for each :math:`m,l`), - - .. math:: - - \begin{aligned} - \psi_{{\rm SSD},ml}(t)&=\delta_{x} \sin\left(\omega_{x} t + 2\pi\frac{mN_{cc,x}}{N_{bx}}\right)\\ - &+\delta_{y} \sin\left(\omega_{y} t + 2\pi\frac{lN_{cc,y}}{N_{by}}\right), - \end{aligned} - - and :math:`A_{ml}=1`. The modulation frequencies :math:`\omega_x,\omega_y` are determined by the - laser bandwidth and modulation amplitudes by the relation - - .. math:: - - \omega_x = \frac{\Delta_\nu r }{2\delta_x} - - where :math:`\Delta_\nu` is the relative bandwidth of the laser pulse and :math:`r` is an additional rotation factor - supplied by the user that determines how much of the modulation is in x and how much is in y. [Michel, Eqn. 9.69] - - * Gaussian Process Randomly phase-modulated SSD, ``'GP RPM SSD'``: - CPP + a generalization of SSD that has temporal stochastic variation in the beamlet phases; that is, :math:`\phi_{{\rm CPP},ml}\in[0,\pi]`, :math:`\psi_{{\rm SSD},ml}(t)` is sampled from a Gaussian stochastic process, and :math:`A_{ml}=1` - * Induced spatial incoherence (ISI), ``'GP ISI'``: - a smoothing technique with temporal stochastic variation in the beamlet phases and amplitudes; that is, :math:`\phi_{{\rm CPP},ml}=0`, :math:`\psi_{{\rm SSD},ml}(t)=0`, and the :math:`A_{ml}` are complex amplitudes sampled from a Gaussian stochastic process to simulate the random phase differences and amplitudes the beamlets pick up when passing through the ISI echelons - - This is an adapation of work by `Han Wen `__ to LASY. - - - Notes - ----- - This assumes a flat-top rectangular laser and so a rectangular arrangement of beamlets in the near-field. - - Parameters - ---------- - wavelength : float (in meters) - The main laser wavelength :math:`\lambda_0` of the laser, which - defines :math:`\omega_0` in the above formula, according to - :math:`\omega_0 = 2\pi c/\lambda_0`. - - pol : list of 2 complex numbers (dimensionless) - Polarization vector. It corresponds to :math:`p_u` in the above - formula ; :math:`p_x` is the first element of the list and - :math:`p_y` is the second element of the list. Using complex - numbers enables elliptical polarizations. - - laser_energy : float (in Joules) - The total energy of the laser pulse. The amplitude of the laser - field (:math:`E_0` in the above formula) is automatically - calculated so that the pulse has the prescribed energy. - - focal_length : float (in meters) - Focal length of lens :math:`f` just after the RPP/CPP. - - beam_aperture : list of 2 floats (in meters) - Widths :math:`D_x,D_y` of the rectangular beam in the near-field, i.e., size of the illuminated region of the RPP/CPP. - - n_beamlets : list of 2 integers - Number of RPP/CPP elements :math:`N_{bx},N_{by}` in each direction, in the near field. - - temporal_smoothing_type : string - Which method for beamlet production and evolution is used. - Can be ``'RPP'``, ``'CPP'``, ``'FM SSD'``, ``'GP RPM SSD'``, or ``'GP ISI'``. - (See the above bullet list for an explanation of each of these methods.) - - relative_laser_bandwidth : float - Bandwidth :math:`\Delta_\nu` of the laser pulse, relative to central frequency. - Only used if ``temporal_smoothing_type`` is ``'FM SSD'``, ``'GP RPM SSD'`` or ``'GP ISI'``. - - ssd_phase_modulation_amplitude :list of 2 floats - Amplitudes :math:`\delta_{x},\delta_{y}` of phase modulation in each transverse direction. - Only used if ``temporal_smoothing_type`` is ``'FM SSD'`` or ``'GP RPM SSD'``. - - ssd_number_color_cycles : list of 2 floats - Number of color cycles :math:`N_{cc,x},N_{cc,y}` of SSD spectrum to include in modulation - Only used if ``temporal_smoothing_type`` is ``'FM SSD'`` or ``'GP RPM SSD'``. - - ssd_transverse_bandwidth_distribution: list of 2 floats - Determines how much SSD is distributed in the :math:`x` and :math:`y` directions. - if `ssd_transverse_bandwidth_distribution=[a,b]`, then the SSD frequency modulation is :math:`a/\sqrt{a^2+b^2}` in :math:`x` and :math:`b/\sqrt{a^2+b^2}` in :math:`y`. - Only used if ``temporal_smoothing_type`` is ``'FM SSD'`` or ``'GP RPM SSD'``. - - do_include_transverse_envelope : boolean, (optional, default False) - Whether to include the transverse sinc envelope or not. - I.e. whether it is assumed to be close enough to the laser axis to neglect the transverse field decay. - """ - - supported_smoothing = "RPP", "CPP", "FM SSD", "GP RPM SSD", "GP ISI" - - def __init__( - self, - wavelength, - pol, - laser_energy, - focal_length, - beam_aperture, - n_beamlets, - temporal_smoothing_type, - relative_laser_bandwidth, - ssd_phase_modulation_amplitude=None, - ssd_number_color_cycles=None, - ssd_transverse_bandwidth_distribution=None, - do_include_transverse_envelope=False, - ): - super().__init__(wavelength, pol) - self.laser_energy = laser_energy - self.focal_length = focal_length - self.beam_aperture = np.array(beam_aperture, dtype="float") - self.n_beamlets = np.array(n_beamlets, dtype="int") - self.temporal_smoothing_type = temporal_smoothing_type - self.laser_bandwidth = relative_laser_bandwidth - - # time interval to update the speckle pattern, roughly update 50 times every bandwidth cycle - self.dt_update = 1 / self.laser_bandwidth / 50 - self.do_include_transverse_envelope = do_include_transverse_envelope - - self.x_lens_list = np.linspace( - -0.5 * (self.n_beamlets[0] - 1), - 0.5 * (self.n_beamlets[0] - 1), - num=self.n_beamlets[0], - ) - self.y_lens_list = np.linspace( - -0.5 * (self.n_beamlets[1] - 1), - 0.5 * (self.n_beamlets[1] - 1), - num=self.n_beamlets[1], - ) - self.Y_lens_matrix, self.X_lens_matrix = np.meshgrid( - self.y_lens_list, self.x_lens_list - ) - self.Y_lens_index_matrix, self.X_lens_index_matrix = np.meshgrid( - np.arange(self.n_beamlets[1], dtype=float), - np.arange(self.n_beamlets[0], dtype=float), - ) - - if "SSD" in self.temporal_smoothing_type.upper(): - # Initialize SSD parameters - # the amplitude of phase along each direction - self.ssd_phase_modulation_amplitude = ssd_phase_modulation_amplitude - # number of color cycles - self.ssd_number_color_cycles = ssd_number_color_cycles - # bandwidth distributed with respect to the two transverse direction - self.ssd_transverse_bandwidth_distribution = ( - ssd_transverse_bandwidth_distribution - ) - ssd_normalization = np.sqrt( - self.ssd_transverse_bandwidth_distribution[0] ** 2 - + self.ssd_transverse_bandwidth_distribution[1] ** 2 - ) - ssd_frac = [ - self.ssd_transverse_bandwidth_distribution[0] / ssd_normalization, - self.ssd_transverse_bandwidth_distribution[1] / ssd_normalization, - ] - self.ssd_phase_modulation_frequency = [ - self.laser_bandwidth * sf * 0.5 / pma - for sf, pma in zip(ssd_frac, self.ssd_phase_modulation_amplitude) - ] - self.ssd_time_delay = ( - ( - self.ssd_number_color_cycles[0] - / self.ssd_phase_modulation_frequency[0] - if self.ssd_phase_modulation_frequency[0] > 0 - else 0 - ), - ( - self.ssd_number_color_cycles[1] - / self.ssd_phase_modulation_frequency[1] - if self.ssd_phase_modulation_frequency[1] > 0 - else 0 - ), - ) - - # Check user input - assert ( - temporal_smoothing_type.upper() in SpeckleProfile.supported_smoothing - ), "Only support one of the following: " + ", ".join( - SpeckleProfile.supported_smoothing - ) - assert relative_laser_bandwidth > 0, "laser_bandwidth must be greater than 0" - assert np.size(n_beamlets) == 2, "has to be a size 2 array" - if "SSD" in self.temporal_smoothing_type.upper(): - assert ( - ssd_number_color_cycles is not None - ), "must supply `ssd_number_color_cycles` to use SSD" - assert ( - ssd_transverse_bandwidth_distribution is not None - ), "must supply `ssd_transverse_bandwidth_distribution` to use SSD" - assert ( - ssd_phase_modulation_amplitude is not None - ), "must supply `ssd_phase_modulation_amplitude` to use SSD" - for q in ( - ssd_number_color_cycles, - ssd_transverse_bandwidth_distribution, - ssd_phase_modulation_amplitude, - ): - assert np.size(q) == 2, "has to be a size 2 array" - assert q[0] > 0 or q[1] > 0, "cannot be all zeros" - - def init_gaussian_time_series( - self, - series_time, - ): - r"""Initialize a time series sampled from a Gaussian process. - - At every time specified by the input `series_time`, calculate the random phase and/or amplitudes as determined by the smoothing type. - - * If the smoothing type is "SSD", then this function returns a time series with random phase offsets in x and y at each time. - The phase offsets are real-valued and centered around the user supplied ``ssd_phase_modulation_amplitude`` - :math:`\delta_{x},\delta_{y}`, with distribution FWHM ``ssd_phase_modulation_frequency``. - - If the smoothing type is "ISI", this function returns a time series with complex numbers defining beamlet phase and amplitude. - Each beamlet has a complex number whose magnitude is the amplitude of the beamlet and - whose phase gives the beamlet phase offset. - The complex numbers have mean 1 and FWHM given by twice the laser bandwidth - - Parameters - ---------- - series_time: array of times at which to sample from Gaussian process - ssd_time_delay: only required for "SSD" type smoothing - ssd_phase_modulation_frequency: only required for "SSD" type smoothing - - Returns - ------- - array-like, either the supplied `series_time` if "ISI" smoothing or `series_time` with some padding at the end for "SSD" smoothing - array-like, either with 2 (for "SSD" smoothing) or `n_beamlets[0] x n_beamlets[1]` ("ISI" smoothing) random numbers at every time - """ - if "SSD" in self.temporal_smoothing_type.upper(): - pm_phase0 = gen_gaussian_time_series( - series_time.size - + int(np.sum(self.ssd_time_delay) / self.dt_update) - + 2, - self.dt_update, - 2 * np.pi * self.ssd_phase_modulation_frequency[0], - self.ssd_phase_modulation_amplitude[0], - ) - pm_phase1 = gen_gaussian_time_series( - series_time.size - + int(np.sum(self.ssd_time_delay) / self.dt_update) - + 2, - self.dt_update, - 2 * np.pi * self.ssd_phase_modulation_frequency[1], - self.ssd_phase_modulation_amplitude[1], - ) - time_interp = np.arange( - start=0, - stop=series_time[-1] + np.sum(self.ssd_time_delay) + 3 * self.dt_update, - step=self.dt_update, - )[: pm_phase0.size] - return ( - time_interp, - [ - (np.real(pm_phase0) + np.imag(pm_phase0)) / np.sqrt(2), - (np.real(pm_phase1) + np.imag(pm_phase1)) / np.sqrt(2), - ], - ) - elif "ISI" in self.temporal_smoothing_type.upper(): - complex_amp = np.stack( - [ - np.stack( - [ - gen_gaussian_time_series( - series_time.size, - self.dt_update, - 2 * self.laser_bandwidth, - 1, - ) - for _i in range(self.n_beamlets[1]) - ] - ) - for _j in range(self.n_beamlets[0]) - ] - ) - return series_time, complex_amp - else: - return series_time, None - - def beamlets_complex_amplitude( - self, t_now, series_time, time_series, temporal_smoothing_type="FM SSD" - ): - """Calculate complex amplitude of the beamlets in the near-field, before propagating to the focal plane. - - If the temporal smoothing type is "RPP" or "CPP", this returns a matrix of ones, giving no modification to the amplitude. - If the temporal smoothing type is "FM SSD", this returns the complex phases as calculated in, for example, Introduction to Laser-Plasma Interactions eqn. 9.87. - If the temporal smoothing type is "GP RPM FM ", this returns complex phases modeled as random variables. - If the temporal smoothing type is "ISI", this returns an array of random complex numbers that gives both amplitude and phase of the beamlets. - - Parameters - ---------- - t_now: float, time at which to calculate the complex amplitude of the beamlets - series_time: 1d array of times at which the stochastic process was sampled to generate the time series - time_series: array of random phase and/or amplitudes as determined by the smoothing type - temporal_smoothing_type: string, what type of temporal smoothing to perform. - - Returns - ------- - array of complex numbers giving beamlet amplitude and phases in the near-field - """ - if any( - rpp_type in temporal_smoothing_type.upper() for rpp_type in ["RPP", "CPP"] - ): - return np.ones_like(self.X_lens_matrix) - if temporal_smoothing_type.upper() == "FM SSD": - phase_t = self.ssd_phase_modulation_amplitude[0] * np.sin( - self.ssd_x_y_dephasing[0] - + 2 - * np.pi - * self.ssd_phase_modulation_frequency[0] - * ( - t_now - - self.X_lens_matrix * self.ssd_time_delay[0] / self.n_beamlets[0] - ) - ) + self.ssd_phase_modulation_amplitude[1] * np.sin( - self.ssd_x_y_dephasing[1] - + 2 - * np.pi - * self.ssd_phase_modulation_frequency[1] - * ( - t_now - - self.Y_lens_matrix * self.ssd_time_delay[1] / self.n_beamlets[1] - ) - ) - return np.exp(1j * phase_t) - elif temporal_smoothing_type.upper() == "GP RPM SSD": - phase_t = np.interp( - t_now - + self.X_lens_index_matrix - * self.ssd_time_delay[0] - / self.n_beamlets[0], - series_time, - time_series[0], - ) + np.interp( - t_now - + self.Y_lens_index_matrix - * self.ssd_time_delay[1] - / self.n_beamlets[1], - series_time, - time_series[1], - ) - return np.exp(1j * phase_t) - elif temporal_smoothing_type.upper() == "GP ISI": - return time_series[:, :, int(round(t_now / self.dt_update))] - else: - raise NotImplementedError - - def generate_speckle_pattern( - self, t_now, exp_phase_plate, x, y, series_time, time_series - ): - """Calculate the speckle pattern in the focal plane. - - Calculates the complex envelope defining the laser pulse in the focal plane at time `t=t_now`. - This function first gets the beamlet complex amplitudes and phases with the function `beamlets_complex_amplitude` - then propagates the the beamlets to the focal plane. - - Parameters - ---------- - t_now: float, time at which to calculate the speckle pattern - exp_phase_plate: 2d array of complex numbers giving the RPP / CPP phase contributions to the beamlets - x: 3d array of x-positions in focal plane - y: 3d array of y-positions in focal plane - series_time: 1d array of times at which the stochastic process was sampled to generate the time series - time_series: array of random phase and/or amplitudes as determined by the smoothing type - - Returns - ------- - speckle_amp: 2D array of complex numbers defining the laser envelope at focus at time `t_now` - """ - lambda_fnum = self.lambda0 * self.focal_length / self.beam_aperture - X_focus_matrix = x[:, :, 0] / lambda_fnum[0] - Y_focus_matrix = y[:, :, 0] / lambda_fnum[1] - x_focus_list = X_focus_matrix[:, 0] - y_focus_list = Y_focus_matrix[0, :] - x_phase_focus_matrix = np.exp( - -2 - * np.pi - * 1j - / self.n_beamlets[0] - * self.x_lens_list[:, np.newaxis] - * x_focus_list[np.newaxis, :] - ) - y_phase_focus_matrix = np.exp( - -2 - * np.pi - * 1j - / self.n_beamlets[1] - * self.y_lens_list[:, np.newaxis] - * y_focus_list[np.newaxis, :] - ) - - bca = self.beamlets_complex_amplitude( - t_now, - series_time=series_time, - time_series=time_series, - temporal_smoothing_type=self.temporal_smoothing_type, - ) - speckle_amp = np.einsum( - "jk,jl->kl", - np.einsum("ij,ik->jk", bca * exp_phase_plate, x_phase_focus_matrix), - y_phase_focus_matrix, - ) - if self.do_include_transverse_envelope: - speckle_amp = ( - np.sinc(X_focus_matrix / self.n_beamlets[0]) - * np.sinc(Y_focus_matrix / self.n_beamlets[1]) - * speckle_amp - ) - return speckle_amp - - def evaluate(self, x, y, t): - """ - Return the envelope field of the laser. - - Parameters - ---------- - x, y, t: ndarrays of floats - Define points on which to evaluate the envelope - These arrays need to all have the same shape. - - Returns - ------- - envelope: ndarray of complex numbers - Contains the value of the envelope at the specified points - This array has the same shape as the arrays x, y, t - """ - # General parameters - t_norm = t[0, 0, :] * c / self.lambda0 - t_max = t_norm[-1] - - # Calculate auxiliary parameters - if "RPP" == self.temporal_smoothing_type.upper(): - phase_plate = np.random.choice([0, np.pi], self.n_beamlets) - elif any( - cpp_smoothing_type in self.temporal_smoothing_type.upper() - for cpp_smoothing_type in ["CPP", "SSD"] - ): - phase_plate = np.random.uniform( - -np.pi, np.pi, size=self.n_beamlets[0] * self.n_beamlets[1] - ).reshape(self.n_beamlets) - elif "ISI" in self.temporal_smoothing_type.upper(): - phase_plate = np.zeros(self.n_beamlets) # ISI does not require phase plates - else: - raise NotImplementedError - exp_phase_plate = np.exp(1j * phase_plate) - if self.temporal_smoothing_type.upper() == "FM SSD": - self.ssd_x_y_dephasing = np.random.standard_normal(2) * np.pi - - series_time = np.arange(0, t_max + self.dt_update, self.dt_update) - - if "GP" in self.temporal_smoothing_type.upper(): - new_series_time, time_series = self.init_gaussian_time_series(series_time) - else: - new_series_time, time_series = series_time, None - - envelope = np.zeros(x.shape, dtype=complex) - for i, t_i in enumerate(t_norm): - envelope[:, :, i] = self.generate_speckle_pattern( - t_i, - exp_phase_plate=exp_phase_plate, - x=x, - y=y, - series_time=new_series_time, - time_series=time_series, - ) - return envelope diff --git a/lasy/profiles/speckled/__init__.py b/lasy/profiles/speckled/__init__.py new file mode 100644 index 00000000..5d015ad1 --- /dev/null +++ b/lasy/profiles/speckled/__init__.py @@ -0,0 +1,13 @@ +from .speckle_profile import SpeckleProfile +from .rpp_cpp import PhasePlateProfile +from .fm_ssd import FM_SSD_Profile +from .gp_rpm_ssd import GP_RPM_SSD_Profile +from .gp_isi import GP_ISI_Profile + +__all__ = [ + "SpeckleProfile", + "PhasePlateProfile", + "FM_SSD_Profile", + "GP_RPM_SSD_Profile", + "GP_ISI_Profile", +] diff --git a/lasy/profiles/speckled/documentation_splice.py b/lasy/profiles/speckled/documentation_splice.py new file mode 100644 index 00000000..b0983f2f --- /dev/null +++ b/lasy/profiles/speckled/documentation_splice.py @@ -0,0 +1,41 @@ +class _DocumentedMetaClass(type): + """Metaclass that combines the __doc__ of the SpeckleProfile base and of the implementation.""" + + def __new__(cls, name, bases, attrs): + # "if bases" skips this for the _ClassWithInit (which has no bases) + # "if bases[0].__doc__ is not None" skips this for the picmistandard classes since their bases[0] (i.e. _ClassWithInit) + # has no __doc__. + if bases and bases[0].__doc__ is not None: + implementation_doc = attrs.get("__doc__", "") + # print('implementation doc',implementation_doc) + base_doc = bases[0].__doc__ + param_delimiter = "Parameters\n ----------\n" + opt_param_delimiter = " do_include_transverse_envelope" + + if implementation_doc: + # The format of the added string is intentional. + # The double return "\n\n" is needed to start a new section in the documentation. + # Then the four spaces matches the standard level of indentation for doc strings + # (assuming PEP8 formatting). + # The final return "\n" assumes that the implementation doc string begins with a return, + # i.e. a line with only three quotes, """. + implementation_notes, implementation_params = implementation_doc.split( + param_delimiter + ) + base_doc_notes, base_doc_params = base_doc.split(param_delimiter) + base_doc_needed_params, base_doc_opt_params = base_doc_params.split( + opt_param_delimiter + ) + attrs["__doc__"] = ( + base_doc_notes + + implementation_notes + + param_delimiter + + base_doc_needed_params + + implementation_params[:-5] + + "\n\n" + + opt_param_delimiter + + base_doc_opt_params + ) + else: + attrs["__doc__"] = base_doc + return super(_DocumentedMetaClass, cls).__new__(cls, name, bases, attrs) diff --git a/lasy/profiles/speckled/fm_ssd.py b/lasy/profiles/speckled/fm_ssd.py new file mode 100644 index 00000000..6687387f --- /dev/null +++ b/lasy/profiles/speckled/fm_ssd.py @@ -0,0 +1,144 @@ +import numpy as np +from .speckle_profile import SpeckleProfile +from .documentation_splice import _DocumentedMetaClass + + +class FM_SSD_Profile(SpeckleProfile, metaclass=_DocumentedMetaClass): + r"""Speckled laser profile information specific to smoothing by frequency modulated (FM) spectral dispersion (SSD). + + In frequency-modulated smoothing by spectral dispersion, or FM-SSD, the amplitude of the beamlets is always :math:`A_{ml}(t)=1`. + There are two contributions to the phase :math:`\phi_{ml}` of each beamlet: + + .. math:: + + \phi_{ml}(t)=\phi_{PP,ml}+\phi_{SSD,ml}. + + The phase plate part :math:`\phi_{PP,ml}` is the initial phase delay from the randomly sized phase plate sections, + drawn from uniform distribution on the interval :math:`[0,2\pi]`. + The temporal smoothing is from the SSD term: + + .. math:: + + \begin{aligned} + \phi_{SSD,ml}(t)&=\delta_{x} \sin\left(\omega_{x} t + 2\pi\frac{mN_{cc,x}}{N_{bx}}\right)\\ + &+\delta_{y} \sin\left(\omega_{y} t + 2\pi\frac{lN_{cc,y}}{N_{by}}\right). + \end{aligned} + + The modulation frequencies :math:`\omega_x,\omega_y` are determined by the + laser bandwidth and modulation amplitudes according to the relation + + .. math:: + + \omega_x = \frac{\Delta_\nu r_x }{2\delta_x}, + \omega_y = \frac{\Delta_\nu r_y }{2\delta_y}, + + where :math:`\Delta_\nu` is the relative bandwidth of the laser pulse + and :math:`r_x, r_y` are additional rotation factors supplied by the user + in the `transverse_bandwidth_distribution` parameter that determine + how much of the modulation is in x and how much is in y. [Michel, Eqn. 9.69] + + Parameters + ---------- + relative_laser_bandwidth : float + Resulting bandwidth :math:`\Delta_\nu` of the laser pulse, relative to central frequency, due to the frequency modulation. + + phase_modulation_amplitude : list of 2 floats + Amplitudes :math:`\delta_{x},\delta_{y}` of phase modulation in each transverse direction. + + number_color_cycles : list of 2 floats + Number of color cycles :math:`N_{cc,x},N_{cc,y}` of SSD spectrum to include in modulation + + transverse_bandwidth_distribution : list of 2 floats + Determines how much SSD is distributed in the :math:`x` and :math:`y` directions. + if `transverse_bandwidth_distribution=[a,b]`, then the SSD frequency modulation is :math:`r_x=a/\sqrt{a^2+b^2}` in :math:`x` and :math:`r_y=b/\sqrt{a^2+b^2}` in :math:`y`. + """ + + def __init__( + self, + wavelength, + pol, + laser_energy, + focal_length, + beam_aperture, + n_beamlets, + relative_laser_bandwidth, + phase_modulation_amplitude, + number_color_cycles, + transverse_bandwidth_distribution, + do_include_transverse_envelope=True, + long_profile=None, + ): + super().__init__( + wavelength, + pol, + laser_energy, + focal_length, + beam_aperture, + n_beamlets, + do_include_transverse_envelope, + long_profile, + ) + self.laser_bandwidth = relative_laser_bandwidth + # the amplitude of phase along each direction + self.phase_modulation_amplitude = phase_modulation_amplitude + # number of color cycles + self.number_color_cycles = number_color_cycles + # bandwidth distributed with respect to the two transverse direction + self.transverse_bandwidth_distribution = transverse_bandwidth_distribution + normalization = np.sqrt( + self.transverse_bandwidth_distribution[0] ** 2 + + self.transverse_bandwidth_distribution[1] ** 2 + ) + frac = [ + self.transverse_bandwidth_distribution[0] / normalization, + self.transverse_bandwidth_distribution[1] / normalization, + ] + self.phase_modulation_frequency = [ + self.laser_bandwidth * sf * 0.5 / pma + for sf, pma in zip(frac, self.phase_modulation_amplitude) + ] + self.time_delay = ( + ( + self.number_color_cycles[0] / self.phase_modulation_frequency[0] + if self.phase_modulation_frequency[0] > 0 + else 0 + ), + ( + self.number_color_cycles[1] / self.phase_modulation_frequency[1] + if self.phase_modulation_frequency[1] > 0 + else 0 + ), + ) + self.x_y_dephasing = np.random.standard_normal(2) * np.pi + self.phase_plate = np.random.uniform( + -np.pi, np.pi, size=self.n_beamlets[0] * self.n_beamlets[1] + ).reshape(self.n_beamlets) + + def beamlets_complex_amplitude( + self, + t_now, + ): + """Calculate complex amplitude of the beamlets in the near-field, before propagating to the focal plane. + + Parameters + ---------- + t_now: float, time at which to evaluate complex amplitude + + Returns + ------- + array of complex numbers giving beamlet amplitude and phases in the near-field + """ + phase_t = self.phase_modulation_amplitude[0] * np.sin( + self.x_y_dephasing[0] + + 2 + * np.pi + * self.phase_modulation_frequency[0] + * (t_now - self.X_lens_matrix * self.time_delay[0] / self.n_beamlets[0]) + ) + self.phase_modulation_amplitude[1] * np.sin( + self.x_y_dephasing[1] + + 2 + * np.pi + * self.phase_modulation_frequency[1] + * (t_now - self.Y_lens_matrix * self.time_delay[1] / self.n_beamlets[1]) + ) + return np.exp(1j * (self.phase_plate + phase_t)) diff --git a/lasy/profiles/speckled/gp_isi.py b/lasy/profiles/speckled/gp_isi.py new file mode 100644 index 00000000..20126167 --- /dev/null +++ b/lasy/profiles/speckled/gp_isi.py @@ -0,0 +1,115 @@ +import numpy as np +from .speckle_profile import SpeckleProfile +from .stochastic_process_utilities import gen_gaussian_time_series +from .documentation_splice import _DocumentedMetaClass + + +class GP_ISI_Profile(SpeckleProfile, metaclass=_DocumentedMetaClass): + r"""Speckled laser profile information specific to smoothing inspired by Induced Spatial Incoherence (ISI). + + This is a smoothing technique with temporal stochastic variation in the beamlet phases and amplitudes + to simulate the random phase differences and amplitudes the beamlets pick up when passing through ISI echelons. + In this case, :math:`\phi_{ml}(t)` and :math:`A_{ml}(t)` are chosen randomly. + Practically, this is done by drawing the complex amplitudes :math:\tilde A_{ml}(t)` + from a stochastic process with Guassian power spectral density having mean of 1 and FWHM of twice the laser bandwidth. + + Parameters + ---------- + relative_laser_bandwidth : float + Bandwidth :math:`\Delta_\nu` of the incoming laser pulse, relative to the central frequency. + + """ + + def __init__( + self, + wavelength, + pol, + laser_energy, + focal_length, + beam_aperture, + n_beamlets, + relative_laser_bandwidth, + do_include_transverse_envelope=True, + long_profile=None, + ): + super().__init__( + wavelength, + pol, + laser_energy, + focal_length, + beam_aperture, + n_beamlets, + do_include_transverse_envelope, + long_profile, + ) + self.laser_bandwidth = relative_laser_bandwidth + self.dt_update = 1 / self.laser_bandwidth / 50 + return + + def init_gaussian_time_series( + self, + series_time, + ): + r"""Initialize a time series sampled from a Gaussian process. + + At every time specified by the input `series_time`, calculate the random phases and/or amplitudes. + + * This function returns a time series with random phase offsets in x and y at each time. + The phase offsets are real-valued and centered around the user supplied ``phase_modulation_amplitude`` + :math:`\delta_{x},\delta_{y}`, with distribution FWHM ``phase_modulation_frequency``. + + Parameters + ---------- + series_time: array of times at which to sample from Gaussian process + + Returns + ------- + array-like, the supplied `series_time` + array-like, either with 2 random numbers at every time + """ + complex_amp = np.stack( + [ + np.stack( + [ + gen_gaussian_time_series( + series_time.size, + self.dt_update, + 2 * self.laser_bandwidth, + 1, + ) + for _i in range(self.n_beamlets[1]) + ] + ) + for _j in range(self.n_beamlets[0]) + ] + ) + return complex_amp + + def setup_for_evaluation(self, t_norm): + """Create or update data used in evaluation.""" + self.x_y_dephasing = np.random.standard_normal(2) * np.pi + self.phase_plate = np.random.uniform( + -np.pi, np.pi, size=self.n_beamlets[0] * self.n_beamlets[1] + ).reshape(self.n_beamlets) + + t_max = t_norm[-1] + series_time = np.arange(0, t_max + self.dt_update, self.dt_update) + + self.time_series = self.init_gaussian_time_series(series_time) + return + + def beamlets_complex_amplitude( + self, + t_now, + ): + """Calculate complex amplitude of the beamlets in the near-field, before propagating to the focal plane. + + Parameters + ---------- + t_now: float, time at which to evaluate complex amplitude + + Returns + ------- + array of complex numbers giving beamlet amplitude and phases in the near-field + """ + return self.time_series[:, :, int(round(t_now / self.dt_update))] diff --git a/lasy/profiles/speckled/gp_rpm_ssd.py b/lasy/profiles/speckled/gp_rpm_ssd.py new file mode 100644 index 00000000..3d981112 --- /dev/null +++ b/lasy/profiles/speckled/gp_rpm_ssd.py @@ -0,0 +1,200 @@ +import numpy as np +from .speckle_profile import SpeckleProfile +from .stochastic_process_utilities import gen_gaussian_time_series +from .documentation_splice import _DocumentedMetaClass + + +class GP_RPM_SSD_Profile(SpeckleProfile, metaclass=_DocumentedMetaClass): + r"""Speckled laser profile information specific to smoothing by a random phase modulated (RPM) spectral dispersion (SSD). + + This provides a version of smoothing by spectral dispersion (SSD) where the phases are randomly modulated. + Here the amplitude of the beamlets is always :math:`A_{ml}(t)=1`. + There are two contributions to the phase :math:`\phi_{ml}` of each beamlet: + + ..math:: + + \phi_{ml}(t) = \phi_{PP,ml} + \phi_{SSD,ml}. + + The phase plate part :math:`\phi_{PP,ml}` is the initial phase delay from the randomly sized phase plate sections, + drawn from uniform distribution on the interval :math:`[0,2\pi]`. + The phases :math:`\phi_{SSD,ml}(t)` are drawn from a stochastic process + with Gaussian power spectrum with means :math:`\delta_x,\delta_y` given by the `phase_modulation_amplitude` argument + and FWHM given by the modulation frequencies :math:`\omega_x,\omega_y`. + The modulation frequencies :math:`\omega_x,\omega_y` are determined by the + laser bandwidth and modulation amplitudes according to the relation + + .. math:: + + \omega_x = \frac{\Delta_\nu r_x }{2\delta_x}, + \omega_y = \frac{\Delta_\nu r_y }{2\delta_y}, + + where :math:`\Delta_\nu` is the resulting relative bandwidth of the laser pulse + and :math:`r_x, r_y` are additional rotation factors supplied by the user + in the `transverse_bandwidth_distribution` parameter that determine + how much of the modulation is in x and how much is in y. [Michel, Eqn. 9.69] + + Parameters + ---------- + relative_laser_bandwidth : float + Bandwidth :math:`\Delta_\nu` of the laser pulse, relative to central frequency. + Only used if ``temporal_smoothing_type`` is ``'FM SSD'``, ``'GP RPM SSD'`` or ``'GP ISI'``. + + phase_modulation_amplitude : list of 2 floats + Amplitudes :math:`\delta_{x},\delta_{y}` of phase modulation in each transverse direction. + Only used if ``temporal_smoothing_type`` is ``'FM SSD'`` or ``'GP RPM SSD'``. + + number_color_cycles : list of 2 floats + Number of color cycles :math:`N_{cc,x},N_{cc,y}` of SSD spectrum to include in modulation + Only used if ``temporal_smoothing_type`` is ``'FM SSD'`` or ``'GP RPM SSD'``. + + transverse_bandwidth_distribution : list of 2 floats + Determines how much SSD is distributed in the :math:`x` and :math:`y` directions. + if `transverse_bandwidth_distribution=[a,b]`, then the SSD frequency modulation is :math:`a/\sqrt{a^2+b^2}` in :math:`x` and :math:`b/\sqrt{a^2+b^2}` in :math:`y`. + Only used if ``temporal_smoothing_type`` is ``'FM SSD'`` or ``'GP RPM SSD'``. + """ + + def __init__( + self, + wavelength, + pol, + laser_energy, + focal_length, + beam_aperture, + n_beamlets, + relative_laser_bandwidth, + phase_modulation_amplitude, + number_color_cycles, + transverse_bandwidth_distribution, + do_include_transverse_envelope=True, + long_profile=None, + ): + super().__init__( + wavelength, + pol, + laser_energy, + focal_length, + beam_aperture, + n_beamlets, + do_include_transverse_envelope, + long_profile, + ) + self.laser_bandwidth = relative_laser_bandwidth + # the amplitude of phase along each direction + self.phase_modulation_amplitude = phase_modulation_amplitude + # number of color cycles + self.number_color_cycles = number_color_cycles + # bandwidth distributed with respect to the two transverse direction + self.transverse_bandwidth_distribution = transverse_bandwidth_distribution + normalization = np.sqrt( + self.transverse_bandwidth_distribution[0] ** 2 + + self.transverse_bandwidth_distribution[1] ** 2 + ) + frac = [ + self.transverse_bandwidth_distribution[0] / normalization, + self.transverse_bandwidth_distribution[1] / normalization, + ] + self.phase_modulation_frequency = [ + self.laser_bandwidth * sf * 0.5 / pma + for sf, pma in zip(frac, self.phase_modulation_amplitude) + ] + self.time_delay = ( + ( + self.number_color_cycles[0] / self.phase_modulation_frequency[0] + if self.phase_modulation_frequency[0] > 0 + else 0 + ), + ( + self.number_color_cycles[1] / self.phase_modulation_frequency[1] + if self.phase_modulation_frequency[1] > 0 + else 0 + ), + ) + self.dt_update = 1 / self.laser_bandwidth / 50 + return + + def init_gaussian_time_series( + self, + series_time, + ): + r"""Initialize a time series sampled from a Gaussian process with the correct power spectral density. + + At every time specified by the input `series_time`, calculate the random phases and/or amplitudes. + + This function returns a time series with random phase offsets in x and y at each time. + The phase offsets are real-valued and centered around the user supplied ``phase_modulation_amplitude`` + :math:`\delta_{x},\delta_{y}`. + The time series has Gaussian power spectral density, or autocorrelation, with a FWHM ``phase_modulation_frequency``. + + Parameters + ---------- + series_time: array of times at which to sample from Gaussian process + time_delay: only required for "SSD" type smoothing + phase_modulation_frequency: only required for "SSD" type smoothing + + Returns + ------- + array-like, the supplied `series_time` with some padding at the end for "SSD" smoothing + array-like, 2 random numbers at every time + """ + pm_phase0 = gen_gaussian_time_series( + series_time.size + int(np.sum(self.time_delay) / self.dt_update) + 2, + self.dt_update, + 2 * np.pi * self.phase_modulation_frequency[0], + self.phase_modulation_amplitude[0], + ) + pm_phase1 = gen_gaussian_time_series( + series_time.size + int(np.sum(self.time_delay) / self.dt_update) + 2, + self.dt_update, + 2 * np.pi * self.phase_modulation_frequency[1], + self.phase_modulation_amplitude[1], + ) + time_interp = np.arange( + start=0, + stop=series_time[-1] + np.sum(self.time_delay) + 3 * self.dt_update, + step=self.dt_update, + )[: pm_phase0.size] + return ( + time_interp, + [ + (np.real(pm_phase0) + np.imag(pm_phase0)) / np.sqrt(2), + (np.real(pm_phase1) + np.imag(pm_phase1)) / np.sqrt(2), + ], + ) + + def setup_for_evaluation(self, t_norm): + """Create or update data used in evaluation.""" + self.x_y_dephasing = np.random.standard_normal(2) * np.pi + self.phase_plate = np.random.uniform( + -np.pi, np.pi, size=self.n_beamlets[0] * self.n_beamlets[1] + ).reshape(self.n_beamlets) + + t_max = t_norm[-1] + series_time = np.arange(0, t_max + self.dt_update, self.dt_update) + + self.series_time, self.time_series = self.init_gaussian_time_series(series_time) + return + + def beamlets_complex_amplitude( + self, + t_now, + ): + """Calculate complex amplitude of the beamlets in the near-field, before propagating to the focal plane. + + Parameters + ---------- + t_now: float, time at which to evaluate complex amplitude + + Returns + ------- + array of complex numbers giving beamlet amplitude and phases in the near-field + """ + phase_t = np.interp( + t_now + self.X_lens_index_matrix * self.time_delay[0] / self.n_beamlets[0], + self.series_time, + self.time_series[0], + ) + np.interp( + t_now + self.Y_lens_index_matrix * self.time_delay[1] / self.n_beamlets[1], + self.series_time, + self.time_series[1], + ) + return np.exp(1j * (self.phase_plate + phase_t)) diff --git a/lasy/profiles/speckled/rpp_cpp.py b/lasy/profiles/speckled/rpp_cpp.py new file mode 100644 index 00000000..a0c4816c --- /dev/null +++ b/lasy/profiles/speckled/rpp_cpp.py @@ -0,0 +1,67 @@ +import numpy as np +from .speckle_profile import SpeckleProfile +from .documentation_splice import _DocumentedMetaClass + + +class PhasePlateProfile(SpeckleProfile, metaclass=_DocumentedMetaClass): + r"""Laser profile information for the random phase plate / continuous phase plate class of speckled lasers. + + This has no temporal smoothing. + The amplitude of the beamlets is always :math:`A_{ml}(t)=1` and + the relative phases of the beamlets, resulting from the randomly sized phase plate sections, + are assigned randomly. + If the user specifies Random Phase Plate (RPP: `rpp`), the beamlet phases :math:`\phi_{ml}(t)` are drawn with equal probabilities from the set :math:`{0,2\pi}`. + If the user specifies Continuous Phase Plate (CPP: `cpp`), the beamlet phases :math:`\phi_{ml}(t)` are drawn from a uniform distribution on the interval :math:`[0,2\pi]`. + + Parameters + ---------- + rpp_cpp: string, can be 'rpp' or 'cpp' + Whether to assign beamlet phases according to RPP or CPP scheme + """ + + def __init__( + self, + wavelength, + pol, + laser_energy, + focal_length, + beam_aperture, + n_beamlets, + rpp_cpp, + do_include_transverse_envelope=True, + long_profile=None, + ): + super().__init__( + wavelength, + pol, + laser_energy, + focal_length, + beam_aperture, + n_beamlets, + do_include_transverse_envelope, + long_profile, + ) + self.rpp_cpp = rpp_cpp + + def beamlets_complex_amplitude( + self, + t_now, + ): + """Calculate complex amplitude of the beamlets in the near-field, before propagating to the focal plane. + + Parameters + ---------- + t_now: float, time at which to evaluate complex amplitude + + Returns + ------- + array of complex numbers giving beamlet amplitude and phases in the near-field + """ + if self.rpp_cpp.upper() == "RPP": + phase_plate = np.random.choice([0, np.pi], self.n_beamlets) + else: + phase_plate = np.random.uniform( + -np.pi, np.pi, size=self.n_beamlets[0] * self.n_beamlets[1] + ).reshape(self.n_beamlets) + exp_phase_plate = np.exp(1j * phase_plate) + return exp_phase_plate diff --git a/lasy/profiles/speckled/speckle_profile.py b/lasy/profiles/speckled/speckle_profile.py new file mode 100644 index 00000000..4a631580 --- /dev/null +++ b/lasy/profiles/speckled/speckle_profile.py @@ -0,0 +1,230 @@ +import numpy as np +from scipy.constants import c +from ..profile import Profile + + +class SpeckleProfile(Profile): + r"""Profile of a speckled laser pulse. + + Speckled lasers are used to mitigate laser-plasma interactions in fusion and ion acceleration contexts. + More on the subject can be found in chapter 9 of `P. Michel, Introduction to Laser-Plasma Interactions `__. + A speckled laser beam is a laser that is deliberately divided transversely into :math:`N_{bx}\times N_{by}` beamlets in the near-field. + The phase plate provides a different phase to each beamlet, with index :math:`ml`, which then propagate to the far field and combine incoherently. + + The electric field in the focal plane, as a function of time :math:`t` and the coordinates + :math:`\boldsymbol{x}_\perp=(x,y)` transverse to the direction of propagation, is: + + .. math:: + + \begin{aligned} + E_u(\boldsymbol{x}_\perp,t) &= Re\left[ E_0 + {\rm sinc}\left(\frac{\pi x}{\Delta x}\right) + {\rm sinc}\left(\frac{\pi y}{\Delta y}\right)\times p_u + \right. + \\ + & \times\sum_{m,l=1}^{N_{bx}, N_{by}} A_{ml}(t) + \exp\left(i\boldsymbol{k}_{\perp ml}\cdot\boldsymbol{x}_\perp + + i\phi_{ml}(t)\right) + \Bigg] + \end{aligned} + + where :math:`u` is either :math:`x` or :math:`y`, :math:`p_u` is + the polarization vector, and :math:`Re` represent the real part [Michel, Eqns. 9.11, 87, 94]. + Several quantities are computed internally to the code depending on the + method of smoothing chosen, including the beamlet amplitude :math:`A_{ml}(t)`, + the beamlet wavenumber :math:`k_{\perp ml}`, + the relative phase contribution :math:`\phi_{ml}(t)` of beamlet :math:`ml` induced by the phase plate and temporal smoothing. + The beam widths are :math:`\Delta x=\frac{\lambda_0fN_{bx}}{D_{x}}`, + :math:`\Delta y=\frac{\lambda_0fN_{by}}{D_{y}}`. + The other parameters in these formulas are defined below. + + Parameters + ---------- + wavelength : float (in meters) + The main laser wavelength :math:`\lambda_0` of the laser, which + defines :math:`\omega_0` in the above formula, according to + :math:`\omega_0 = 2\pi c/\lambda_0`. + + pol : list of 2 complex numbers (dimensionless) + Polarization vector. It corresponds to :math:`p_u` in the above + formula ; :math:`p_x` is the first element of the list and + :math:`p_y` is the second element of the list. Using complex + numbers enables elliptical polarizations. + + laser_energy : float (in Joules) + The total energy of the laser pulse. The amplitude of the laser + field (:math:`E_0` in the above formula) is automatically + calculated so that the pulse has the prescribed energy. + + focal_length : float (in meters) + Focal length of lens :math:`f` just after the RPP/CPP. + + beam_aperture : list of 2 floats (in meters) + Widths :math:`D_x,D_y` of the rectangular beam in the near-field, i.e., size of the illuminated region of the RPP/CPP. + + n_beamlets : list of 2 integers + Number of RPP/CPP elements :math:`N_{bx},N_{by}` in each direction, in the near field. + + do_include_transverse_envelope : boolean (optional, default: False) + Whether to include the transverse sinc envelope or not. + I.e. whether it is assumed to be close enough to the laser axis to neglect the transverse field decay. + + long_profile : Lasy Longitudinal laser object (optional, default: None). + If this is not None, the longitudinal profile is applied individually to the beamlets in the near-field. + + Notes + ----- + This is an adaptation of work by `Han Wen `__ to LASY. + + This assumes a flat-top rectangular laser and so a rectangular arrangement of beamlets in the near-field. + The longitudinal profile is currently applied to the beamlets + individually in the near-field before they are propagated to the focal plane. + """ + + def __init__( + self, + wavelength, + pol, + laser_energy, + focal_length, + beam_aperture, + n_beamlets, + do_include_transverse_envelope, + long_profile, + ): + super().__init__(wavelength, pol) + self.laser_energy = laser_energy + self.focal_length = focal_length + self.beam_aperture = np.array(beam_aperture, dtype="float") + self.n_beamlets = np.array(n_beamlets, dtype="int") + self.do_include_transverse_envelope = do_include_transverse_envelope + self.long_profile = long_profile + + self.x_lens_list = np.linspace( + -0.5 * (self.n_beamlets[0] - 1), + 0.5 * (self.n_beamlets[0] - 1), + num=self.n_beamlets[0], + ) + self.y_lens_list = np.linspace( + -0.5 * (self.n_beamlets[1] - 1), + 0.5 * (self.n_beamlets[1] - 1), + num=self.n_beamlets[1], + ) + self.Y_lens_matrix, self.X_lens_matrix = np.meshgrid( + self.y_lens_list, self.x_lens_list + ) + self.Y_lens_index_matrix, self.X_lens_index_matrix = np.meshgrid( + np.arange(self.n_beamlets[1], dtype=float), + np.arange(self.n_beamlets[0], dtype=float), + ) + return + + def beamlets_complex_amplitude( + self, + t_now, + ): + """Calculate complex amplitude of the beamlets in the near-field, before propagating to the focal plane. + + This function can be overwritten to define custom speckled laser objects. + + Parameters + ---------- + t_now: float, time at which to evaluate complex amplitude + + Returns + ------- + array of complex numbers giving beamlet amplitude and phases in the near-field + """ + return np.ones_like(self.X_lens_matrix) + + def setup_for_evaluation(self, t_norm): + """Create or update data used in evaluation.""" + pass + + def generate_speckle_pattern(self, t_now, x, y): + """Calculate the speckle pattern in the focal plane. + + Calculates the complex envelope defining the laser pulse in the focal plane at time `t=t_now`. + This function first gets the beamlet complex amplitudes and phases with the function `beamlets_complex_amplitude` + then propagates the the beamlets to the focal plane. + + Parameters + ---------- + t_now: float, time at which to calculate the speckle pattern + exp_phase_plate: 2d array of complex numbers giving the RPP / CPP phase contributions to the beamlets + x: 3d array of x-positions in focal plane + y: 3d array of y-positions in focal plane + series_time: 1d array of times at which the stochastic process was sampled to generate the time series + time_series: array of random phase and/or amplitudes as determined by the smoothing type + + Returns + ------- + speckle_amp: 2D array of complex numbers defining the laser envelope at focus at time `t_now` + """ + lambda_fnum = self.lambda0 * self.focal_length / self.beam_aperture + X_focus_matrix = x[:, :, 0] / lambda_fnum[0] + Y_focus_matrix = y[:, :, 0] / lambda_fnum[1] + x_focus_list = X_focus_matrix[:, 0] + y_focus_list = Y_focus_matrix[0, :] + x_phase_focus_matrix = np.exp( + -2 + * np.pi + * 1j + / self.n_beamlets[0] + * self.x_lens_list[:, np.newaxis] + * x_focus_list[np.newaxis, :] + ) + y_phase_focus_matrix = np.exp( + -2 + * np.pi + * 1j + / self.n_beamlets[1] + * self.y_lens_list[:, np.newaxis] + * y_focus_list[np.newaxis, :] + ) + bca = self.beamlets_complex_amplitude(t_now) + if self.long_profile is not None: + # have to unnormalize t_now to evaluate in longitudinal profile + bca = bca * self.long_profile.evaluate(t_now / c * self.lambda0) + # propagate from near-field to focus + speckle_amp = np.einsum( + "jk,jl->kl", + np.einsum("ij,ik->jk", bca, x_phase_focus_matrix), + y_phase_focus_matrix, + ) + if self.do_include_transverse_envelope: + speckle_amp = ( + np.sinc(X_focus_matrix / self.n_beamlets[0]) + * np.sinc(Y_focus_matrix / self.n_beamlets[1]) + * speckle_amp + ) + return speckle_amp + + def evaluate(self, x, y, t): + """ + Return the envelope field of the laser. + + Parameters + ---------- + x, y, t: ndarrays of floats + Define points on which to evaluate the envelope + These arrays need to all have the same shape. + + Returns + ------- + envelope: ndarray of complex numbers + Contains the value of the envelope at the specified points + This array has the same shape as the arrays x, y, t + """ + # General parameters + t_norm = t[0, 0, :] * c / self.lambda0 + self.setup_for_evaluation(t_norm) + + envelope = np.zeros(x.shape, dtype=complex) + for i, t_i in enumerate(t_norm): + envelope[:, :, i] = self.generate_speckle_pattern( + t_i, + x=x, + y=y, + ) + return envelope diff --git a/lasy/profiles/speckled/stochastic_process_utilities.py b/lasy/profiles/speckled/stochastic_process_utilities.py new file mode 100644 index 00000000..93890987 --- /dev/null +++ b/lasy/profiles/speckled/stochastic_process_utilities.py @@ -0,0 +1,33 @@ +import numpy as np + + +def gen_gaussian_time_series(t_num, dt, fwhm, rms_mean): + """Generate a discrete time series that has gaussian power spectrum. + + Credit Han Wen, possibly NRL + + Parameters + ---------- + t_num: number of grid points in time + fwhm: full width half maximum of the power spectrum + rms_mean: root-mean-square average of the spectrum + + Returns + ------- + temporal_amplitude: a time series array of complex numbers with shape [t_num] + """ + if fwhm == 0.0: + temporal_amplitude = np.zeros(t_num, dtype=np.complex128) + else: + omega = np.fft.fftshift(np.fft.fftfreq(t_num, d=dt)) + psd = np.exp(-np.log(2) * 0.5 * np.square(omega / fwhm * 2 * np.pi)) + spectral_amplitude = np.array(psd) * ( + np.random.normal(size=t_num) + 1j * np.random.normal(size=t_num) + ) + temporal_amplitude = np.fft.ifftshift( + np.fft.fft(np.fft.fftshift(spectral_amplitude)) + ) + temporal_amplitude *= rms_mean / np.sqrt( + np.mean(np.square(np.abs(temporal_amplitude))) + ) + return temporal_amplitude diff --git a/tests/test_laser_profiles.py b/tests/test_laser_profiles.py index eca27bce..b3e33801 100644 --- a/tests/test_laser_profiles.py +++ b/tests/test_laser_profiles.py @@ -7,7 +7,7 @@ from lasy.laser import Laser from lasy.profiles.profile import Profile, SummedProfile, ScaledProfile -from lasy.profiles import GaussianProfile, FromArrayProfile, SpeckleProfile +from lasy.profiles import GaussianProfile, FromArrayProfile from lasy.profiles.longitudinal import ( GaussianLongitudinalProfile, SuperGaussianLongitudinalProfile, @@ -24,6 +24,7 @@ SummedTransverseProfile, ScaledTransverseProfile, ) +from lasy.profiles.speckled import GP_ISI_Profile from lasy.utils.exp_data_utils import find_center_of_mass @@ -314,18 +315,20 @@ def test_speckle_profile(): focal_length = 3.5 # m beam_aperture = [0.35, 0.5] # m n_beamlets = [24, 32] - temporal_smoothing_type = "GP ISI" + do_sinc_profile = False relative_laser_bandwidth = 0.005 + long_profile = None - profile = SpeckleProfile( + profile = GP_ISI_Profile( wavelength, polarization, laser_energy, focal_length, beam_aperture, n_beamlets, - temporal_smoothing_type=temporal_smoothing_type, - relative_laser_bandwidth=relative_laser_bandwidth, + relative_laser_bandwidth, + do_sinc_profile, + long_profile, ) dimensions = "xyt" dx = wavelength * focal_length / beam_aperture[0] diff --git a/tests/test_speckles.py b/tests/test_speckles.py index 63dac7b8..4ec91542 100644 --- a/tests/test_speckles.py +++ b/tests/test_speckles.py @@ -1,11 +1,56 @@ import numpy as np - from lasy.laser import Laser -from lasy.profiles.speckle_profile import SpeckleProfile +from lasy.profiles.speckled import ( + FM_SSD_Profile, + GP_ISI_Profile, + GP_RPM_SSD_Profile, + PhasePlateProfile, +) import pytest from scipy.constants import c +def _get_arg_string( + temporal_smoothing_type, + speckle_args, + ssd_args=None, + isi_args=None, +): + if temporal_smoothing_type.upper() in ["RPP", "CPP"]: + args = [*speckle_args, temporal_smoothing_type] + elif temporal_smoothing_type.upper() in ["FM SSD", "GP RPM SSD"]: + if ssd_args is None: + raise ValueError(f"require ssd_args for SSD smoothing") + else: + args = [*speckle_args, *ssd_args] + elif temporal_smoothing_type.upper() == "GP ISI": + if isi_args is None: + raise ValueError(f"require isi_args for ISI smoothing") + else: + args = [*speckle_args, *isi_args] + else: + raise ValueError(f"Invalid smoothing type provided: {temporal_smoothing_type}") + return args + + +def _get_laser_profile( + temporal_smoothing_type, + *args, + **kw_args, +): + if temporal_smoothing_type.upper() in ["RPP", "CPP"]: + profile = PhasePlateProfile(*args, **kw_args) + elif temporal_smoothing_type.upper() in "FM SSD": + profile = FM_SSD_Profile(*args, **kw_args) + elif temporal_smoothing_type.upper() == "GP RPM SSD": + profile = GP_RPM_SSD_Profile(*args, **kw_args) + elif temporal_smoothing_type.upper() == "GP ISI": + profile = GP_ISI_Profile(*args, **kw_args) + else: + raise ValueError(f"Invalid smoothing type provided: {temporal_smoothing_type}") + return profile + + @pytest.mark.parametrize( "temporal_smoothing_type", ["RPP", "CPP", "FM SSD", "GP RPM SSD", "GP ISI"] ) @@ -18,29 +63,40 @@ def test_intensity_distribution(temporal_smoothing_type): wavelength = 0.351e-6 # Laser wavelength in meters polarization = (1, 0) # Linearly polarized in the x direction + laser_energy = 1.0 # J (this is the laser energy stored in the box defined by `lo` and `hi` below) focal_length = 3.5 # m beam_aperture = [0.35, 0.5] # m n_beamlets = [24, 32] - relative_laser_bandwidth = 0.005 - laser_energy = 1.0 # J (this is the laser energy stored in the box defined by `lo` and `hi` below) - - ssd_phase_modulation_amplitude = (4.1, 4.5) - ssd_number_color_cycles = [1.4, 1.0] - ssd_transverse_bandwidth_distribution = [1.8, 1.0] - - profile = SpeckleProfile( + do_sinc_profile = False + long_profile = None + speckle_args = ( wavelength, polarization, laser_energy, focal_length, beam_aperture, n_beamlets, - temporal_smoothing_type=temporal_smoothing_type, - relative_laser_bandwidth=relative_laser_bandwidth, - ssd_phase_modulation_amplitude=ssd_phase_modulation_amplitude, - ssd_number_color_cycles=ssd_number_color_cycles, - ssd_transverse_bandwidth_distribution=ssd_transverse_bandwidth_distribution, ) + opt_args = { + "do_include_transverse_envelope": do_sinc_profile, + "long_profile": long_profile, + } + + relative_laser_bandwidth = 0.005 + phase_modulation_amplitude = (4.1, 4.5) + number_color_cycles = [1.4, 1.0] + transverse_bandwidth_distribution = [1.8, 1.0] + ssd_args = ( + relative_laser_bandwidth, + phase_modulation_amplitude, + number_color_cycles, + transverse_bandwidth_distribution, + ) + isi_args = (relative_laser_bandwidth,) + + args = _get_arg_string(temporal_smoothing_type, speckle_args, ssd_args, isi_args) + profile = _get_laser_profile(temporal_smoothing_type, *args, **opt_args) + dimensions = "xyt" dx = wavelength * focal_length / beam_aperture[0] dy = wavelength * focal_length / beam_aperture[1] @@ -92,25 +148,36 @@ def test_spatial_correlation(temporal_smoothing_type): focal_length = 3.5 # m beam_aperture = [0.35, 0.35] # m n_beamlets = [24, 32] - relative_laser_bandwidth = 0.005 - - ssd_phase_modulation_amplitude = (4.1, 4.1) - ssd_number_color_cycles = [1.4, 1.0] - ssd_transverse_bandwidth_distribution = [1.0, 1.0] - - profile = SpeckleProfile( + do_sinc_profile = False + long_profile = None + speckle_args = ( wavelength, polarization, laser_energy, focal_length, beam_aperture, n_beamlets, - temporal_smoothing_type=temporal_smoothing_type, - relative_laser_bandwidth=relative_laser_bandwidth, # 0.005 - ssd_phase_modulation_amplitude=ssd_phase_modulation_amplitude, - ssd_number_color_cycles=ssd_number_color_cycles, - ssd_transverse_bandwidth_distribution=ssd_transverse_bandwidth_distribution, ) + opt_args = { + "do_include_transverse_envelope": do_sinc_profile, + "long_profile": long_profile, + } + + relative_laser_bandwidth = 0.005 + phase_modulation_amplitude = (4.1, 4.1) + number_color_cycles = [1.4, 1.0] + transverse_bandwidth_distribution = [1.0, 1.0] + ssd_args = ( + relative_laser_bandwidth, + phase_modulation_amplitude, + number_color_cycles, + transverse_bandwidth_distribution, + ) + isi_args = (relative_laser_bandwidth,) + + args = _get_arg_string(temporal_smoothing_type, speckle_args, ssd_args, isi_args) + profile = _get_laser_profile(temporal_smoothing_type, *args, **opt_args) + dimensions = "xyt" dx = wavelength * focal_length / beam_aperture[0] dy = wavelength * focal_length / beam_aperture[1] @@ -175,25 +242,35 @@ def test_sinc_zeros(temporal_smoothing_type): focal_length = 3.5 # m beam_aperture = [0.35, 0.35] # m n_beamlets = [24, 48] - relative_laser_bandwidth = 0.005 - ssd_phase_modulation_amplitude = (4.1, 4.1) - ssd_number_color_cycles = [1.4, 1.0] - ssd_transverse_bandwidth_distribution = [1.0, 1.0] - - profile = SpeckleProfile( + do_sinc_profile = True + long_profile = None + speckle_args = ( wavelength, polarization, laser_energy, focal_length, beam_aperture, n_beamlets, - temporal_smoothing_type=temporal_smoothing_type, - relative_laser_bandwidth=relative_laser_bandwidth, - ssd_phase_modulation_amplitude=ssd_phase_modulation_amplitude, - ssd_number_color_cycles=ssd_number_color_cycles, - ssd_transverse_bandwidth_distribution=ssd_transverse_bandwidth_distribution, - do_include_transverse_envelope=True, ) + opt_args = { + "do_include_transverse_envelope": do_sinc_profile, + "long_profile": long_profile, + } + + relative_laser_bandwidth = 0.005 + phase_modulation_amplitude = (4.1, 4.1) + number_color_cycles = [1.4, 1.0] + transverse_bandwidth_distribution = [1.0, 1.0] + ssd_args = ( + relative_laser_bandwidth, + phase_modulation_amplitude, + number_color_cycles, + transverse_bandwidth_distribution, + ) + isi_args = (relative_laser_bandwidth,) + + args = _get_arg_string(temporal_smoothing_type, speckle_args, ssd_args, isi_args) + profile = _get_laser_profile(temporal_smoothing_type, *args, **opt_args) dimensions = "xyt" dx = wavelength * focal_length / beam_aperture[0] dy = wavelength * focal_length / beam_aperture[1] @@ -215,46 +292,55 @@ def test_sinc_zeros(temporal_smoothing_type): assert abs(F[:, -1, :]).max() / abs(F).max() < 1.0e-8 -def test_FM_SSD_periodicity(): +def test_FM_periodicity(): """Test that the frequency modulated Smoothing by spectral dispersion (SSD) has the correct temporal frequency.""" + temporal_smoothing_type = "FM SSD" wavelength = 0.351e-6 # Laser wavelength in meters polarization = (1, 0) # Linearly polarized in the x direction laser_energy = 1.0 # J (this is the laser energy stored in the box defined by `lo` and `hi` below) focal_length = 3.5 # m beam_aperture = [0.35, 0.35] # m n_beamlets = [24, 32] - temporal_smoothing_type = "FM SSD" - relative_laser_bandwidth = 0.005 - - ssd_phase_modulation_amplitude = [4.1, 4.1] - ssd_number_color_cycles = [1.4, 1.0] - ssd_transverse_bandwidth_distribution = [1.0, 1.0] - - laser_profile = SpeckleProfile( + do_sinc_profile = False + long_profile = None + speckle_args = ( wavelength, polarization, laser_energy, focal_length, beam_aperture, n_beamlets, - temporal_smoothing_type=temporal_smoothing_type, - relative_laser_bandwidth=relative_laser_bandwidth, - ssd_phase_modulation_amplitude=ssd_phase_modulation_amplitude, - ssd_number_color_cycles=ssd_number_color_cycles, - ssd_transverse_bandwidth_distribution=ssd_transverse_bandwidth_distribution, ) + opt_args = { + "do_include_transverse_envelope": do_sinc_profile, + "long_profile": long_profile, + } + + relative_laser_bandwidth = 0.005 + phase_modulation_amplitude = [4.1, 4.1] + number_color_cycles = [1.4, 1.0] + transverse_bandwidth_distribution = [1.0, 1.0] + ssd_args = ( + relative_laser_bandwidth, + phase_modulation_amplitude, + number_color_cycles, + transverse_bandwidth_distribution, + ) + args = _get_arg_string(temporal_smoothing_type, speckle_args, ssd_args=ssd_args) + laser_profile = _get_laser_profile(temporal_smoothing_type, *args, **opt_args) + nu_laser = c / wavelength - ssd_frac = np.sqrt( - ssd_transverse_bandwidth_distribution[0] ** 2 - + ssd_transverse_bandwidth_distribution[1] ** 2 + frac = np.sqrt( + transverse_bandwidth_distribution[0] ** 2 + + transverse_bandwidth_distribution[1] ** 2 ) - ssd_frac = ( - ssd_transverse_bandwidth_distribution[0] / ssd_frac, - ssd_transverse_bandwidth_distribution[1] / ssd_frac, + frac = ( + transverse_bandwidth_distribution[0] / frac, + transverse_bandwidth_distribution[1] / frac, ) phase_mod_freq = [ relative_laser_bandwidth * sf * 0.5 / pma - for sf, pma in zip(ssd_frac, ssd_phase_modulation_amplitude) + for sf, pma in zip(frac, phase_modulation_amplitude) ] t_max = 1.0 / phase_mod_freq[0] / nu_laser