Added tests for LongitudinalProfileFromData and fix for issue #185

lasy/profiles/longitudinal/longitudinal_profile_from_data.py

@@ -25,12 +25,15 @@ class LongitudinalProfileFromData(LongitudinalProfile):
 
         datatype : string
             The domain in which the data has been passed. Options
-            are 'spectral' and 'temporal'
+            are 'spectral' and 'temporal'. Spectral data may either be,
+            a) Wavelength dependent (lambda): data["axis_is_wavelength"] = True or
+            b) Angular freuqency dependent (omega): data["axis_is_wavelength"] = False
 
         axis : ndarrays of floats
             The horizontal axis of the pulse duration measurement.
-            The array must be monotonously increasing.
-            When datatype is 'spectral' axis is wavelength in meters.
+            The array must be monotonically increasing or decreasing.
+            When datatype is 'spectral' axis is wavelength in meters OR
+            angular frequency in 1/seconds.
             When datatype is 'temporal' axis is time in seconds.
 
         intensity : ndarrays of floats
@@ -60,16 +63,27 @@ class LongitudinalProfileFromData(LongitudinalProfile):
 
     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"]
+            if "axis_is_wavelength" not in data:  # Set to wavelength data by default
+                data["axis_is_wavelength"] = True
+            if data["axis_is_wavelength"]:
+                wavelength = data["axis"]  # Accept as wavelength
+                spectral_intensity = data["intensity"]
+            else:
+                wavelength = 2.0 * np.pi * c / data["axis"]  # Convert to wavelength
+                spectral_intensity = (
+                    data["intensity"] * 2.0 * np.pi * c / wavelength**2
+                )  # Convert spectral data
+            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 data.get("phase") is None:
                 spectral_phase = np.zeros_like(wavelength)
             else:
                 spectral_phase = data["phase"]
+            if np.all(np.diff(wavelength) < 0):  # Flip arrays
+                wavelength = wavelength[::-1]
+                spectral_intensity = spectral_intensity[::-1]
+                spectral_phase = spectral_phase[::-1]
             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 @@ def test_longitudinal_profiles():
 
     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
 
     t = np.linspace(t_peak - 4 * tau_fwhm, t_peak + 4 * tau_fwhm, npoints)
+    omega = np.linspace(omega_0 - 4 * omega_fwhm, omega_0 + 4 * omega_fwhm, npoints)
+    wavelength_axis = 2.0 * np.pi * c / omega  # Note: monotonically decreasing
 
     # GaussianLongitudinalProfile
     print("GaussianLongitudinalProfile")
@@ -234,6 +238,91 @@ def test_longitudinal_profiles():
     print("cep_phase = ", cep_phase_cos)
     assert np.abs(cep_phase_cos - cep_phase) / cep_phase < 0.02
 
+    # LongitudinalProfileFromData
+    print("LongitudinalProfileFromData")
+    data = {}  # Generate spectral data assuming analytic Fourier transform of GaussianLongitudinalProfile
+    data["datatype"] = "spectral"
+    data["dt"] = 1e-16
+    profile = np.exp(
+        -(tau**2) * ((omega - omega_0) ** 2) / 4.0 + 1.0j * (cep_phase + omega * t_peak)
+    )
+    spectral_intensity = np.abs(profile) ** 2 / np.max(np.abs(profile) ** 2)
+    spectral_phase = np.unwrap(np.angle(profile))
+
+    print("Case 1: monotonically decreasing data on wavelength axis")
+    data["axis"] = wavelength_axis
+    data["intensity"] = spectral_intensity
+    data["phase"] = spectral_phase
+    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 increasing data on wavelength axis")
+    data["axis"] = wavelength_axis[::-1]
+    data["intensity"] = spectral_intensity[::-1]
+    data["phase"] = spectral_phase[::-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
+
+    print("Case 3: monotonically increasing data on angular frequency axis")
+    data["axis"] = omega
+    data["intensity"] = spectral_intensity
+    data["phase"] = spectral_phase
+    data["axis_is_wavelength"] = False
+    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 4: monotonically decreasing data on angular frequency axis")
+    data["axis"] = omega[::-1]
+    data["intensity"] = spectral_intensity[::-1]
+    data["phase"] = spectral_phase[::-1]
+    data["axis_is_wavelength"] = False
+    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