more utilities

src/scgenerator/data/gas.toml

@@ -1,51 +1,7 @@
-[nitrogen]
-a = 0.137
-b = 1.709e-5
-
-[helium]
-a = 0.00346
-b = 2.38e-5
-
-[helium_alt]
-a = 0.00346
-b = 2.38e-5
-
-[hydrogen]
-a = 0.02453
-b = 2.651e-5
-
-[neon]
-a = 0.02135
-b = 1.709e-5
-
-[argon]
-a = 0.1355
-b = 3.201e-5
-
-[argon_alt]
-a = 0.1355
-b = 3.201e-5
-
-[argon_alt2]
-a = 0.1355
-b = 3.201e-5
-
-[krypton]
-a = 0.2349
-b = 3.978e-5
-
-[xenon]
-a = 0.425
-b = 5.105e-5
-
 [air]
 a = 0.1358
 b = 3.64e-5
 
-[vacuum]
-a = 0
-b = 0
-
 [air.sellmeier]
 B = [57921050000.0, 1679170000.0]
 C = [238018500000000.0, 57362000000000.0]
@@ -56,7 +12,12 @@ kind = 2
 [air.kerr]
 P0 = 101325
 T0 = 293.15
-n2 = 3.01e-23
+n2 = 4.0e-23
+source = "Pigeon, J. J., Tochitsky, S. Y., Welch, E. C., & Joshi, C. (2016). Measurements of the nonlinear refractive index of air, N 2, and O 2 at 10 μm using four-wave mixing. Optics letters, 41(17), 3924-3927."
+
+[nitrogen]
+a = 0.137
+b = 1.709e-5
 
 [nitrogen.sellmeier]
 B = [32431570000.0]
@@ -72,6 +33,10 @@ T0 = 273.15
 n2 = 2.2e-23
 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904."
 
+[helium]
+a = 0.00346
+b = 2.38e-5
+
 [helium.sellmeier]
 B = [2.16463842e-5, 2.10561127e-7, 4.7509272e-5]
 C = [-6.80769781e-16, 5.13251289e-15, 3.18621354e-15]
@@ -86,6 +51,10 @@ T0 = 273.15
 n2 = 3.1e-25
 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904."
 
+[helium_alt]
+a = 0.00346
+b = 2.38e-5
+
 [helium_alt.sellmeier]
 B = [14755297000.0]
 C = [426297400000000.0]
@@ -99,6 +68,10 @@ P0 = 30400.0
 T0 = 273.15
 n2 = 3.1e-25
 
+[hydrogen]
+a = 0.02453
+b = 2.651e-5
+
 [hydrogen.sellmeier]
 B = [0.0148956, 0.0049037]
 C = [1.807e-10, 9.2e-11]
@@ -113,6 +86,10 @@ T0 = 273.15
 n2 = 6.36e-24
 source = "Shelton, D. P., & Rice, J. E. (1994). Measurements and calculations of the hyperpolarizabilities of atoms and small molecules in the gas phase. Chemical Reviews, 94(1), 3-29"
 
+[neon]
+a = 0.02135
+b = 1.709e-5
+
 [neon.sellmeier]
 B = [1281450000.0, 22048600000.0]
 C = [184661000000000.0, 376840000000000.0]
@@ -126,6 +103,10 @@ T0 = 273.15
 n2 = 8.7e-25
 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904."
 
+[argon]
+a = 0.1355
+b = 3.201e-5
+
 [argon.sellmeier]
 B = [2501410000.0, 500283000.0, 52234300000.0]
 C = [91012000000000.0, 87892000000000.0, 214020000000000.0]
@@ -140,6 +121,10 @@ T0 = 273.15
 n2 = 9.7e-24
 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904."
 
+[argon_alt]
+a = 0.1355
+b = 3.201e-5
+
 [argon_alt.sellmeier]
 B = [0.0002033229, 0.0003445831]
 C = [2.0612e-16, 8.066e-15]
@@ -154,6 +139,10 @@ T0 = 273.15
 n2 = 9.7e-24
 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904."
 
+[argon_alt2]
+a = 0.1355
+b = 3.201e-5
+
 [argon_alt2.sellmeier]
 B = [0.030182943]
 C = [144000000000000.0]
@@ -169,6 +158,10 @@ T0 = 273.15
 n2 = 9.7e-24
 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904."
 
+[krypton]
+a = 0.2349
+b = 3.978e-5
+
 [krypton.sellmeier]
 B = [2536370000.0, 2736490000.0, 62080200000.0]
 C = [65474200000000.0, 73698000000000.0, 181080000000000.0]
@@ -182,6 +175,10 @@ T0 = 273.15
 n2 = 2.2e-23
 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904."
 
+[xenon]
+a = 0.425
+b = 5.105e-5
+
 [xenon.sellmeier]
 B = [3228690000.0, 3553930000.0, 60676400000.0]
 C = [46301000000000.0, 59578000000000.0, 112740000000000.0]
@@ -195,6 +192,10 @@ T0 = 273.15
 n2 = 5.8e-23
 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904."
 
+[vacuum]
+a = 0
+b = 0
+
 [vacuum.sellmeier]
 B = []
 C = []

src/scgenerator/io.py

@@ -162,7 +162,7 @@ def load_config(path: os.PathLike) -> BareConfig:
     return BareConfig(**config)
 
 
-def load_material_dico(name):
+def load_material_dico(name: str) -> dict[str, Any]:
     """loads a material dictionary
     Parameters
     ----------

src/scgenerator/physics/__init__.py

@@ -0,0 +1,97 @@
+"""
+This file contains functions that don't properly fit into the 'fiber' or 'pulse' category.
+They're also not necessarily 'pure' function as they do some io and stuff.
+"""
+
+from typing import TypeVar
+
+import numpy as np
+from scipy.optimize import minimize_scalar
+
+from ..io import load_material_dico
+from .. import math
+from . import fiber, materials, units, pulse
+
+T = TypeVar("T")
+
+
+def group_delay_to_gdd(wavelength: np.ndarray, group_delay: np.ndarray) -> np.ndarray:
+    w = units.m.inv(wavelength)
+    gdd = np.gradient(group_delay, w)
+    return gdd
+
+
+def material_dispersion(
+    wavelengths,
+    material: str,
+    pressure=None,
+    temperature=None,
+    ideal=False,
+):
+    """returns the dispersion profile (beta_2) of a bulk material.
+
+    Parameters
+    ----------
+    wavelengths : ndarray, shape (n, )
+        wavelengths over which to calculate the dispersion
+    material : str
+        material name in lower case
+    temperature : float, optional
+        Temperature of the material
+    pressure : float, optional
+        constant pressure
+    ideal : bool, optional
+        whether to use the ideal gas law instead of the van der Waals equation, by default False
+
+    Returns
+    -------
+    out : ndarray, shape (n, )
+        beta2 as function of wavelength
+    """
+
+    w = units.m(wavelengths)
+    material_dico = load_material_dico(material)
+    if ideal:
+        n_gas_2 = materials.sellmeier(wavelengths, material_dico, pressure, temperature) + 1
+    else:
+        N_1 = materials.number_density_van_der_waals(
+            pressure=pressure, temperature=temperature, material_dico=material_dico
+        )
+        N_0 = materials.number_density_van_der_waals(material_dico=material_dico)
+        n_gas_2 = materials.sellmeier(wavelengths, material_dico) * N_1 / N_0 + 1
+
+    return fiber.beta2(w, np.sqrt(n_gas_2))
+
+
+def find_optimal_depth(
+    spectrum: T, w_c: np.ndarray, w0: float, material: str = "silica", max_z: float = 1.0
+) -> tuple[T, pulse.OptimizeResult]:
+    """finds the optimal silica depth to compress a pulse
+
+    Parameters
+    ----------
+    spectrum : np.ndarray or Spectrum, shape (n, )
+        spectrum from which to remove 2nd order dispersion
+    w_c : np.ndarray, shape (n, )
+        corresponding centered angular frequencies (w-w0)
+    w0 : float
+        pump central angular frequency
+    material : str
+        material name, by default 'silica'
+    max_z : float
+        maximum propagation distance in m
+
+    Returns
+    -------
+    float
+        distance in m
+    """
+    silica_disp = material_dispersion(units.To.m(w_c + w0), material)
+
+    propagate = lambda z: spectrum * np.exp(-0.5j * silica_disp * w_c ** 2 * z)
+
+    def score(z):
+        return 1 / np.max(math.abs2(np.fft.ifft(propagate(z))))
+
+    opti = minimize_scalar(score, bracket=(0, max_z))
+    return opti.x

src/scgenerator/physics/pulse.py

@@ -13,6 +13,7 @@ import itertools
 import os
 from pathlib import Path
 from typing import Literal, Tuple, TypeVar
+from collections import namedtuple
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -48,6 +49,11 @@ P0T0_to_E0_fac = dict(
 )
 """relates the total energy (amplitue^2) to the t0 parameter of the amplitude and the peak intensity (peak_amplitude^2)"""
 
+PulseProperties = namedtuple(
+    "PulseProperties",
+    "quality mean_coherence fwhm_noise mean_fwhm peak_rin energy_rin timing_jitter",
+)
+
 
 def initial_field(t, shape, t0, peak_power):
     """returns the initial field
@@ -475,12 +481,14 @@ def g12(values):
     # Create all the possible pairs of values
     n = len(values)
     field_pairs = itertools.combinations(values, 2)
+    mean_spec = np.mean(abs2(values), axis=0)
+    mask = mean_spec > 1e-15 * mean_spec.max()
     corr = np.zeros_like(values[0])
     for pair in field_pairs:
-        corr += pair[0].conj() * pair[1]
+        corr[mask] += pair[0][mask].conj() * pair[1][mask]
-    g12_arr = corr / (n * (n - 1) / 2 * np.mean(abs2(values), axis=0))
+    corr[mask] = corr[mask] / (n * (n - 1) / 2 * mean_spec[mask])
 
-    return np.abs(g12_arr)
+    return np.abs(corr)
 
 
 def avg_g12(values):
@@ -842,39 +850,45 @@ def _detailed_find_lobe_limits(
     )
 
 
-def measure_properties(spectra, t, compress=True, debug=""):
+def measure_properties(spectra, t, compress=True, return_limits=False, debug="") -> PulseProperties:
     """measure the quality factor, the fwhm variation, the peak power variation,
 
     Parameters
     ----------
-        spectra : 2D array
+    spectra : np.ndarray, shape (n, nt)
-            set of n spectra in sc-ordering that differ only by noise
+        set of n spectra on a grid of nt angular frequency points
-        t : 1D array
+    t : np.ndarray, shape (nt, )
-            time axis of the simulation
+        time axis of the simulation
-        compress : bool, optional
+    compress : bool, optional
-            whether to perform pulse compression. Default value is True, but this
+        whether to perform pulse compression. Default value is True, but this
-            should be set to False to measure the initial pulse as output by gaussian_pulse
+        should be set to False to measure the initial pulse as output by gaussian_pulse
-            or sech_pulse because compressing it would result in glitches and wrong measurements
+        or sech_pulse because compressing it would result in glitches and wrong measurements
+    return_limits : bool, optional
+        return the time values of the limits
 
     Returns
     ----------
-        qf : float
+    PulseProperties object:
+        quality : float
             quality factor of the pulse ensemble
-        mean_g12 : float
+        mean_gmean_coherence12 : float
             mean coherence of the spectra ensemble
-        fwhm_var : float
+        fwhm_noise : float
-            relative noise in temporal width of the compressed pulse
+            relative noise in temporal width of the (compressed) pulse
-        fwhm_abs : float
+        mean_fwhm : float
-            width of the mean compressed pulse
+            width of the mean (compressed) pulse
-        int_var : flaot
+        peak_rin : float
-            relative noise in the compressed pulse peak intensity
+            relative noise in the (compressed) pulse peak intensity
-        t_jitter : float
+        energy_rin : float
+            relative noise in the (compressed) pulse total_energy
+        timing_jitter : float
             standard deviantion in absolute temporal peak position
+    all_limits : list[tuple[np.ndarray, np.ndarray]], only if return_limits = True
+        list of tuples of the form ([left_lobe_lim, right_lobe_lim, lobe_pos], [left_hm, right_hm])
     """
     if compress:
         fields = abs2(compress_pulse(spectra))
     else:
-        print("Skipping compression")
         fields = abs2(ifft(spectra))
 
     field = np.mean(fields, axis=0)
@@ -897,16 +911,28 @@ def measure_properties(spectra, t, compress=True, debug=""):
     P0 = []
     fwhm = []
     t_offset = []
+    energies = []
+    all_lims: list[tuple[np.ndarray, np.ndarray]] = []
     for f in fields:
         lobe_lim, fwhm_lim, _, big_spline = find_lobe_limits(t, f, debug)
+        all_lims.append((lobe_lim, fwhm_lim))
         P0.append(big_spline(lobe_lim[2]))
         fwhm.append(length(fwhm_lim))
         t_offset.append(lobe_lim[2])
+        energies.append(np.trapz(fields, t))
     fwhm_var = np.std(fwhm) / np.mean(fwhm)
     int_var = np.std(P0) / np.mean(P0)
+    en_var = np.std(energies) / np.mean(energies)
     t_jitter = np.std(t_offset)
 
-    return qf, mean_g12, fwhm_var, fwhm_abs, int_var, t_jitter
+    if isinstance(mean_g12, np.ndarray) and mean_g12.ndim == 0:
+        mean_g12 = mean_g12[()]
+    pp = PulseProperties(qf, mean_g12, fwhm_var, fwhm_abs, int_var, en_var, t_jitter)
+
+    if return_limits:
+        return pp, all_lims
+    else:
+        return pp
 
 
 def measure_field(t: np.ndarray, field: np.ndarray) -> Tuple[float, float, float]: