Added tests for LongitudinalProfileFromData and fix for issue #185

lasy/profiles/longitudinal/longitudinal_profile_from_data.py

@@ -29,7 +29,7 @@ class LongitudinalProfileFromData(LongitudinalProfile):
 
         axis : ndarrays of floats
             The horizontal axis of the pulse duration measurement.
-            The array must be monotonously increasing.
+            The array must be monotonically increasing or decreasing.
             When datatype is 'spectral' axis is wavelength in meters.
             When datatype is 'temporal' axis is time in seconds.
 
@@ -62,14 +62,21 @@ def __init__(self, data, lo, hi):
         if data["datatype"] == "spectral":
             # First find central frequency
             wavelength = data["axis"]
-            assert np.all(
-                np.diff(wavelength) > 0
-            ), 'data["axis"] must be in monotonously increasing order.'
-            spectral_intensity = data["intensity"]
+            assert np.all(np.diff(wavelength) > 0) or np.all(
+                np.diff(wavelength) < 0
+            ), 'data["axis"] must be in monotonically increasing or decreasing order.'
+            if np.all(np.diff(wavelength) < 0):
+                wavelength = wavelength[::-1]
+                spectral_intensity = data["intensity"][::-1]
+            else:
+                spectral_intensity = data["intensity"]
             if data.get("phase") is None:
                 spectral_phase = np.zeros_like(wavelength)
             else:
-                spectral_phase = data["phase"]
+                if np.all(np.diff(wavelength) < 0):
+                    spectral_phase = data["phase"][::-1]
+                else:
+                    spectral_phase = data["phase"]
             dt = data["dt"]
             cwl = np.sum(spectral_intensity * wavelength) / np.sum(spectral_intensity)
             cfreq = c / cwl

tests/test_laser_profiles.py

@@ -9,6 +9,7 @@
 from lasy.profiles.longitudinal import (
     CosineLongitudinalProfile,
     GaussianLongitudinalProfile,
+    LongitudinalProfileFromData,
     SuperGaussianLongitudinalProfile,
 )
 from lasy.profiles.profile import Profile, ScaledProfile, SummedProfile
@@ -153,11 +154,14 @@
 
     wavelength = 800e-9
     tau_fwhm = 30.0e-15
+    omega_fwhm = 4 * np.log(2) / tau_fwhm  # Assumes fully-compressed
     t_peak = 1.0 * tau_fwhm
     cep_phase = 0.5 * np.pi
-    omega_0 = 2.0 * np.pi * c / wavelength
+    omega0 = 2.0 * np.pi * c / wavelength
 
     t = np.linspace(t_peak - 4 * tau_fwhm, t_peak + 4 * tau_fwhm, npoints)
+    omega = np.linspace(omega0 - 4 * omega_fwhm, omega0 + 4 * omega_fwhm, npoints)
+    wavelength_axis = 2.0 * np.pi * c / omega
 
     # GaussianLongitudinalProfile
     print("GaussianLongitudinalProfile")
@@ -234,6 +238,51 @@
     print("cep_phase = ", cep_phase_cos)
     assert np.abs(cep_phase_cos - cep_phase) / cep_phase < 0.02
 
+    # LongitudinalProfileFromData - initially assuming zero phase
+    print("LongitudinalProfileFromData")
+    data = {}
+    data["datatype"] = "spectral"
+    data["dt"] = np.abs(
+        t[1] - t[0]
+    )  # Generate spectral data assuming unchirped Gaussian
+    profile = np.exp(
+        tau**2 * ((omega - omega0) ** 2) / 4.0 + +1.0j * (cep_phase + omega * t_peak)
+    )
+
+    print("Case 1: monotonically increasing spectral data")
+    data["axis"] = wavelength_axis
+    data["intensity"] = spectral_intensity
+    profile_data = LongitudinalProfileFromData(data, np.min(t), np.max(t))
+    field_data = profile_data.evaluate(t)
+
+    std_gauss_data = np.sqrt(np.average((t - t_peak) ** 2, weights=np.abs(field_data)))
+    std_gauss_th = tau / np.sqrt(2.0)
+    print("std_th = ", std_gauss_th)
+    print("std = ", std_gauss_data)
+    assert np.abs(std_gauss_data - std_gauss_th) / std_gauss_th < 0.01
+
+    t_peak_gaussian_data = t[np.argmax(np.abs(field_data))]
+    print("t_peak_th = ", t_peak)
+    print("t_peak = ", t_peak_gaussian_data)
+    assert np.abs(t_peak_gaussian_data - t_peak) / t_peak < 0.01
+
+    print("Case 2: monotonically decreasing spectral data")
+    data["axis"] = wavelength_axis[::-1]
+    data["intensity"] = spectral_intensity[::-1]
+    profile_data = LongitudinalProfileFromData(data, np.min(t), np.max(t))
+    field_data = profile_data.evaluate(t)
+
+    std_gauss_data = np.sqrt(np.average((t - t_peak) ** 2, weights=np.abs(field_data)))
+    std_gauss_th = tau / np.sqrt(2.0)
+    print("std_th = ", std_gauss_th)
+    print("std = ", std_gauss_data)
+    assert np.abs(std_gauss_data - std_gauss_th) / std_gauss_th < 0.01
+
+    t_peak_gaussian_data = t[np.argmax(np.abs(field_data))]
+    print("t_peak_th = ", t_peak)
+    print("t_peak = ", t_peak_gaussian_data)
+    assert np.abs(t_peak_gaussian_data - t_peak) / t_peak < 0.01
+
 
 def test_profile_gaussian_3d_cartesian(gaussian):
     # - 3D Cartesian case