Polynomial Spectral Phase element (#263)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Remi Lehe <[email protected]>

docs/source/api/optical_elements/index.rst

@@ -5,3 +5,4 @@ Optical elements
    :maxdepth: 1
 
    parabolic_mirror
+   polynomial_spectral_phase

docs/source/api/optical_elements/polynomial_spectral_phase.rst

@@ -0,0 +1,5 @@
+Polynomial spectral phase
+=========================
+
+.. autoclass:: lasy.optical_elements.PolynomialSpectralPhase
+    :members:

lasy/optical_elements/__init__.py

@@ -1,3 +1,4 @@
 from .parabolic_mirror import ParabolicMirror
+from .polynomial_spectral_phase import PolynomialSpectralPhase
 
-__all__ = ["ParabolicMirror"]
+__all__ = ["ParabolicMirror", "PolynomialSpectralPhase"]

lasy/optical_elements/optical_element.py

@@ -27,19 +27,24 @@ def amplitude_multiplier(self, x, y, omega, omega0):
 
             \tilde{\mathcal{E}}_{out}(x, y, \omega) = T(x, y, \omega)\tilde{\mathcal{E}}_{in}(x, y, \omega)
 
+        Parameters
+        ----------
+        Return the amplitude multiplier.
+
         Parameters
         ----------
         x, y, omega : ndarrays of floats
             Define points on which to evaluate the multiplier.
             These arrays need to all have the same shape.
-        omega0: float
-            Central angular frequency of the laser.
+        omega0 : float (in rad/s)
+            Central angular frequency, as used for the definition
+            of the laser envelope.
 
         Returns
         -------
         multiplier : ndarray of complex numbers
-            Contains the value of the multiplier at the specified points
-            This array has the same shape as the arrays x, y, omega
+            Contains the value of the multiplier at the specified points.
+            This array has the same shape as the array omega.
         """
         # The base class only defines dummy multiplier
         # (This should be replaced by any class that inherits from this one.)

lasy/optical_elements/parabolic_mirror.py

@@ -34,14 +34,17 @@ def amplitude_multiplier(self, x, y, omega, omega0):
 
         Parameters
         ----------
-        x, y, omega: ndarrays of floats
+        x, y, omega : ndarrays of floats
             Define points on which to evaluate the multiplier.
             These arrays need to all have the same shape.
+        omega0 : float (in rad/s)
+            Central angular frequency, as used for the definition
+            of the laser envelope.
 
         Returns
         -------
-        multiplier: ndarray of complex numbers
-            Contains the value of the multiplier at the specified points
-            This array has the same shape as the arrays x, y, omega
+        multiplier : ndarray of complex numbers
+            Contains the value of the multiplier at the specified points.
+            This array has the same shape as the array omega.
         """
         return np.exp(-1j * omega * (x**2 + y**2) / (2 * c * self.f))

lasy/optical_elements/polynomial_spectral_phase.py

@@ -0,0 +1,68 @@
+import numpy as np
+
+from .optical_element import OpticalElement
+
+
+class PolynomialSpectralPhase(OpticalElement):
+    r"""
+    Class for an optical element that adds spectral phase (e.g. a dazzler).
+
+    The amplitude multiplier corresponds to:
+
+    .. math::
+
+        T(\omega) = \exp(i(\phi(\omega)))
+
+    where :math:`\phi(\omega)` is the spectral phase given by:
+
+    .. math::
+
+        \phi(\omega) = \frac{\text{GDD}}{2!} (\omega - \omega_0)^2 + \frac{\text{TOD}}{3!} (\omega - \omega_0)^3 + \frac{\text{FOD}}{4!} (\omega - \omega_0)^4
+
+    The other parameters in this formula are defined below.
+
+    Parameters
+    ----------
+    gdd : float (in s^2), optional
+        Group Delay Dispersion (by default: ``gdd=0``). ``gdd > 0`` corresponds to a positive
+        chirp, i.e. the low-frequency part of the spectrum arriving earlier than the
+        high-frequency part of the spectrum.
+    tod : float (in s^3), optional
+        Third-order Dispersion (by default: ``tod=0``). For a Gaussian pulse, adding a positive
+        TOD (``tod > 0``) results in the apparition of post-pulses, i.e. lower intensity pulses
+        arriving after the main pulse.
+    fod : float (in s^4), optional
+        Fourth-order Dispersion (by default: ``fod=0``).
+    """
+
+    def __init__(self, gdd=0, tod=0, fod=0):
+        self.gdd = gdd
+        self.tod = tod
+        self.fod = fod
+
+    def amplitude_multiplier(self, x, y, omega, omega0):
+        """
+        Return the amplitude multiplier.
+
+        Parameters
+        ----------
+        x, y, omega : ndarrays of floats
+            Define points on which to evaluate the multiplier.
+            These arrays need to all have the same shape.
+        omega0 : float (in rad/s)
+            Central angular frequency, as used for the definition
+            of the laser envelope.
+
+        Returns
+        -------
+        multiplier : ndarray of complex numbers
+            Contains the value of the multiplier at the specified points.
+            This array has the same shape as the array omega.
+        """
+        spectral_phase = (
+            self.gdd / 2 * (omega - omega0) ** 2
+            + self.tod / 6 * (omega - omega0) ** 3
+            + self.fod / 24 * (omega - omega0) ** 4
+        )
+
+        return np.exp(1j * spectral_phase)

tests/test_polynomial_spectral_phase.py

@@ -0,0 +1,110 @@
+# -*- coding: utf-8 -*-
+"""
+This test checks the implementation of the polynomial spectral phase
+by initializing a Gaussian pulse (with flat spectral phase),
+adding spectral phase to it, and checking the corresponding
+temporal shape of the laser pulse again analytical formulas.
+"""
+
+import numpy as np
+
+from lasy.laser import Laser
+from lasy.optical_elements import PolynomialSpectralPhase
+from lasy.profiles.gaussian_profile import GaussianProfile
+
+# Laser parameters
+wavelength = 0.8e-6
+w0 = 5.0e-6  # m
+pol = (1, 0)
+laser_energy = 1.0  # J
+t_peak = 0.0e-15  # s
+tau = 15.0e-15  # s
+gaussian_profile = GaussianProfile(wavelength, pol, laser_energy, w0, tau, t_peak)
+
+# Grid parameters
+dim = "xyt"
+lo = (-12e-6, -12e-6, -100e-15)
+hi = (+12e-6, +12e-6, +100e-15)
+npoints = (100, 100, 200)
+
+
+def test_gdd():
+    """
+    Add GDD to the laser and compare the on-axis field with the
+    analytical formula for a Gaussian pulse with GDD.
+    """
+    gdd = 200e-30
+    dazzler = PolynomialSpectralPhase(gdd=gdd)
+
+    # Initialize the laser
+    dim = "xyt"
+    lo = (-12e-6, -12e-6, -100e-15)
+    hi = (+12e-6, +12e-6, +100e-15)
+    npoints = (100, 100, 200)
+    laser = Laser(dim, lo, hi, npoints, gaussian_profile)
+
+    # Get field before and after dazzler
+    E_before = laser.grid.get_temporal_field()
+    laser.apply_optics(dazzler)
+    E_after = laser.grid.get_temporal_field()
+    # Extract peak field before dazzler
+    E0 = abs(E_before[50, 50]).max()
+
+    # Compute the analtical expression in real space for a Gaussian
+    t = np.linspace(laser.grid.lo[-1], laser.grid.hi[-1], laser.grid.npoints[-1])
+    stretch_factor = 1 - 2j * gdd / tau**2
+    E_analytical = (
+        E0 * np.exp(-1.0 / stretch_factor * (t / tau) ** 2) / stretch_factor**0.5
+    )
+
+    # Compare the on-axis field with the analytical formula
+    tol = 1.2e-3
+    assert np.all(
+        abs(E_after[50, 50, :] - E_analytical) / abs(E_analytical).max() < tol
+    )
+
+
+def test_tod():
+    """
+    Add TOD to the laser and compare the on-axis field with the
+    analytical formula from the stationary phase approximation.
+    """
+    tod = 5e-42
+    dazzler = PolynomialSpectralPhase(tod=tod)
+
+    # Initialize the laser
+    dim = "xyt"
+    lo = (-12e-6, -12e-6, -50e-15)
+    hi = (+12e-6, +12e-6, +250e-15)
+    npoints = (100, 100, 400)
+    laser = Laser(dim, lo, hi, npoints, gaussian_profile)
+    t = np.linspace(laser.grid.lo[-1], laser.grid.hi[-1], laser.grid.npoints[-1])
+
+    # Get field before and after dazzler
+    E_before = laser.grid.get_temporal_field()
+    laser.apply_optics(dazzler)
+    E_after = laser.grid.get_temporal_field()
+    # Extract peak field before dazzler
+    E0 = abs(E_before[50, 50]).max()
+
+    # Compare data in the post-pulse region to the stationary phase approximation.
+    # The stationary phase approximation result was obtained under the additional
+    # assumption t >> tau^4/tod, so we only compare the data in this region
+    E_compare = abs(E_after[50, 50, t > t_peak + 2 * tau**4 / tod])  # On-axis field
+    t = t[t > t_peak + 2 * tau**4 / tod]
+    prediction = abs(
+        2
+        * E0
+        * tau
+        / (8 * tod * t) ** 0.25
+        * np.exp(-(tau**2) * t / (2 * tod))
+        * np.cos(
+            2 * t / 3 * (2 * t / tod) ** 0.5
+            - tau**4 / (8 * tod) * (2 * t / tod) ** 0.5
+            - np.pi / 4
+        )
+    )
+
+    # Compare the on-axis field with the analytical formula
+    tol = 2.4e-2
+    assert np.all(abs(E_compare - prediction) / abs(E0) < tol)