started working on plasma

README.md

@@ -118,13 +118,13 @@ capillary_nested : int, optional
    
 
 gamma: float, optional unless beta is directly provided
-    nonlinear parameter in m^2 / W. Will overwrite any computed gamma parameter.
+    nonlinear parameter in m^-1*W^-1. Will overwrite any computed gamma parameter.
 
 effective_mode_diameter : float, optional
     effective mode field diameter in m
 
 n2 : float, optional
-    non linear refractive index
+    non linear refractive index in m^2/W
 
 A_eff : float, optional
     effective mode field area

src/scgenerator/__init__.py

@@ -10,6 +10,6 @@ from .math import abs2, argclosest, span
 from .physics import fiber, materials, pulse, simulate, units
 from .physics.simulate import RK4IP, new_simulation, resume_simulations
 from .plotting import plot_avg, plot_results_1D, plot_results_2D, plot_spectrogram
-from .spectra import Pulse
+from .spectra import Pulse, Spectrum
 from .utils.parameter import BareParams, BareConfig
 from . import utils, io, initialize, math

src/scgenerator/physics/__init__.py

@@ -26,7 +26,8 @@ def material_dispersion(
     material: str,
     pressure=None,
     temperature=None,
-    ideal=False,
+    ideal=True,
+    safe=True,
 ):
     """returns the dispersion profile (beta_2) of a bulk material.
 
@@ -41,7 +42,7 @@ def material_dispersion(
     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
+        whether to use the ideal gas law instead of the van der Waals equation, by default True
 
     Returns
     -------
@@ -50,6 +51,9 @@ def material_dispersion(
     """
 
     w = units.m(wavelengths)
+
+    order = np.argsort(w)
+
     material_dico = load_material_dico(material)
     if ideal:
         n_gas_2 = materials.sellmeier(wavelengths, material_dico, pressure, temperature) + 1
@@ -59,12 +63,19 @@ def material_dispersion(
         )
         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
-
+    if safe:
-    return fiber.beta2(w, np.sqrt(n_gas_2))
+        disp = np.zeros(len(w))
+        ind = w > 0
+        disp[ind] = material_dispersion(
+            units.To.m(w[ind]), material, pressure, temperature, ideal, False
+        )
+        return disp
+    else:
+        return fiber.beta2(w[order], np.sqrt(n_gas_2[order]))[order]
 
 
 def find_optimal_depth(
-    spectrum: T, w_c: np.ndarray, w0: float, material: str = "silica", max_z: float = 1.0
+    spectrum: T, w_c: np.ndarray, w0: float, material: str, max_z: float = 1.0
 ) -> tuple[T, pulse.OptimizeResult]:
     """finds the optimal silica depth to compress a pulse
 
@@ -86,12 +97,26 @@ def find_optimal_depth(
     float
         distance in m
     """
-    silica_disp = material_dispersion(units.To.m(w_c + w0), material)
+    w = w_c + w0
+    disp = np.zeros(len(w))
+    ind = w > (w0 / 10)
+    disp[ind] = material_dispersion(units.To.m(w[ind]), material)
 
-    propagate = lambda z: spectrum * np.exp(-0.5j * silica_disp * w_c ** 2 * z)
+    propagate = lambda z: spectrum * np.exp(-0.5j * disp * w_c ** 2 * z)
+    integrate = lambda z: math.abs2(np.fft.ifft(propagate(z)))
 
     def score(z):
-        return 1 / np.max(math.abs2(np.fft.ifft(propagate(z))))
+        return -np.nansum(integrate(z) ** 6)
 
-    opti = minimize_scalar(score, bracket=(0, max_z))
+    # import matplotlib.pyplot as plt
-    return opti.x
+
+    # to_test = np.linspace(0, max_z, 200)
+    # scores = [score(z) for z in to_test]
+    # fig, ax = plt.subplots()
+    # ax.plot(to_test, scores / np.min(scores))
+    # plt.show()
+    # plt.close(fig)
+    # ama = np.argmin(scores)
+
+    opti = minimize_scalar(score, method="bounded", bounds=(0, max_z))
+    return propagate(opti.x), opti

src/scgenerator/physics/fiber.py

@@ -647,7 +647,7 @@ def PCF_dispersion(lambda_, pitch, ratio_d, w0=None, n2=None, A_eff=None):
         if A_eff is None:
             V_eff = pi2a / lambda_ * np.sqrt(n_co ** 2 - n_FSM2)
             w_eff = a_eff * (0.65 + 1.619 / V_eff ** 1.5 + 2.879 / V_eff ** 6)
-            A_eff = interp1d(lambda_, w_eff, kind="linear")(units.m.inv(w0)) ** pipi
+            A_eff = interp1d(lambda_, w_eff, kind="linear")(units.m.inv(w0)) ** 2 * pi
 
         if n2 is None:
             n2 = 2.6e-20

src/scgenerator/physics/materials.py

@@ -1,8 +1,13 @@
+from typing import Any, Callable
 import numpy as np
+import scipy.special
+from scipy.integrate import cumulative_trapezoid
+
+from scgenerator import math
 
 from ..logger import get_logger
 from . import units
-from .units import NA, c, kB, me, e
+from .units import NA, c, kB, me, e, hbar
 
 
 def pressure_from_gradient(ratio, p0, p1):
@@ -174,3 +179,57 @@ def non_linear_refractive_index(material_dico, pressure=None, temperature=None):
 
 def adiabadicity(w: np.ndarray, I: float, field: np.ndarray) -> np.ndarray:
     return w * np.sqrt(2 * me * I) / (e * np.abs(field))
+
+
+def free_electron_density(
+    t: np.ndarray, field: np.ndarray, N0: float, rate_func: Callable[[np.ndarray], np.ndarray]
+) -> np.ndarray:
+    return N0 * (1 - np.exp(-cumulative_trapezoid(rate_func(field), t, initial=0)))
+
+
+def ionization_rate_ADK(
+    ionization_energy: float, atomic_number
+) -> Callable[[np.ndarray], np.ndarray]:
+    Z = -(atomic_number - 1) * e
+
+    omega_p = ionization_energy / hbar
+    nstar = Z * np.sqrt(2.1787e-18 / ionization_energy)
+    omega_t = lambda field: e * np.abs(field) / np.sqrt(2 * me * ionization_energy)
+    Cnstar = 2 ** (2 * nstar) / (scipy.special.gamma(nstar + 1) ** 2)
+    omega_C = omega_p / 4 * math.abs2(Cnstar)
+
+    def rate(field: np.ndarray) -> np.ndarray:
+        opt4 = 4 * omega_t(field) / om
+        return omega_C * (4 * omega_p / ot) ** (2 * nstar - 1) * np.exp(-4 * omega_p / (3 * ot))
+
+    return rate
+
+
+class Plasma:
+    def __init__(self, t: np.ndarray, ionization_energy: float, atomic_number: int):
+        self.t = t
+        self.Ip = ionization_energy
+        self.atomic_number = atomic_number
+        self.rate = ionization_rate_ADK(self.Ip, self.atomic_number)
+
+    def __call__(self, field: np.ndarray, N0: float) -> np.ndarray:
+        """returns the number density of free electrons as function of time
+
+        Parameters
+        ----------
+        field : np.ndarray
+            electric field in V/m
+        N0 : float
+            total number density of matter
+
+        Returns
+        -------
+        np.ndarray
+            number density of free electrons as function of time
+        """
+        Ne = free_electron_density(self.t, field, N0, self.rate)
+        return cumulative_trapezoid(
+            np.gradient(Ne, self.t) * self.Ip / field, self.t, initial=0
+        ) + e ** 2 / me * cumulative_trapezoid(
+            cumulative_trapezoid(Ne * field, self.t, initial=0), self.t, initial=0
+        )

src/scgenerator/physics/pulse.py

@@ -572,7 +572,7 @@ def peak_ind(values, mam=None):
         ) + 3
     except IndexError:
         right_ind = len(values) - 1
-    return left_ind, right_ind
+    return max(0, left_ind), min(len(values) - 1, right_ind)
 
 
 def setup_splines(x_axis, values, mam=None):
@@ -685,8 +685,7 @@ def find_lobe_limits(x_axis, values, debug="", already_sorted=True):
             dd_roots,
             fwhm_pos,
             peak_pos,
-            folder_name,
+            out_path,
-            file_name,
             fig,
             ax,
             color,
@@ -716,7 +715,7 @@ def find_lobe_limits(x_axis, values, debug="", already_sorted=True):
             c=color[5],
         )
         ax.legend()
-        fig.savefig(os.path.join(folder_name, file_name), bbox_inches="tight")
+        fig.savefig(out_path, bbox_inches="tight")
         plt.close(fig)
 
     else:
@@ -818,9 +817,7 @@ def _detailed_find_lobe_limits(
     # if measurement of the peak is not straightforward, we plot the situation to see
     # if the final measurement is good or not
 
-    folder_name, file_name, fig, ax = plot_setup(
+    out_path, fig, ax = plot_setup(out_path=f"measurement_errors_plots/it_{iterations}_{debug}")
-        file_name=f"it_{iterations}_{debug}", folder_name="measurements_errors_plots"
-    )
 
     new_fwhm_pos = np.array([np.max(left_pos), np.min(right_pos)])
 
@@ -842,8 +839,7 @@ def _detailed_find_lobe_limits(
         dd_roots,
         fwhm_pos,
         peak_pos,
-        folder_name,
+        out_path,
-        file_name,
         fig,
         ax,
         color,
@@ -946,7 +942,7 @@ def measure_field(t: np.ndarray, field: np.ndarray) -> Tuple[float, float, float
 
 
 def remove_2nd_order_dispersion(
-    spectrum: T, w_c: np.ndarray, beta2: float, max_z: float = -1.0
+    spectrum: T, w_c: np.ndarray, beta2: float, max_z: float = -100.0
 ) -> tuple[T, OptimizeResult]:
     """attempts to remove 2nd order dispersion from a complex spectrum
 
@@ -964,11 +960,49 @@ def remove_2nd_order_dispersion(
     np.ndarray, shape (n, )
         spectrum with 2nd order dispersion removed
     """
-    # makeshift_t = np.linspace(0, 1, len(w_c))
     propagate = lambda z: spectrum * np.exp(-0.5j * beta2 * w_c ** 2 * z)
 
     def score(z):
-        return 1 / np.max(abs2(np.fft.ifft(propagate(z))))
+        return -np.max(abs2(np.fft.ifft(propagate(z))))
 
     opti = minimize_scalar(score, bracket=(max_z, 0))
     return propagate(opti.x), opti
+
+
+def remove_2nd_order_dispersion2(
+    spectrum: T, w_c: np.ndarray, max_gdd: float = 1000e-30
+) -> tuple[T, OptimizeResult]:
+    """attempts to remove 2nd order dispersion from a complex spectrum
+
+    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)
+
+    Returns
+    -------
+    np.ndarray, shape (n, )
+        spectrum with 2nd order dispersion removed
+    """
+    propagate = lambda gdd: spectrum * np.exp(-0.5j * w_c ** 2 * 1e-30 * gdd)
+    integrate = lambda gdd: abs2(np.fft.ifft(propagate(gdd)))
+
+    def score(gdd):
+        return -np.sum(integrate(gdd) ** 6)
+
+    # def score(gdd):
+    #     return -np.max(integrate(gdd))
+
+    # to_test = np.linspace(-max_gdd, max_gdd, 200)
+    # scores = [score(g) for g in to_test]
+    # fig, ax = plt.subplots()
+    # ax.plot(to_test, scores / np.min(scores))
+    # plt.show()
+    # plt.close(fig)
+    # ama = np.argmin(scores)
+
+    opti = minimize_scalar(score, bounds=(-max_gdd * 1e30, max_gdd * 1e30))
+    opti["x"] *= 1e-30
+    return propagate(opti.x * 1e30), opti

src/scgenerator/physics/simulate.py

@@ -1,5 +1,6 @@
 import multiprocessing
 import os
+import random
 from datetime import datetime
 from pathlib import Path
 from typing import Dict, List, Tuple, Type
@@ -704,7 +705,7 @@ def new_simulation(
     if prev_sim_dir is not None:
         config_dict["prev_sim_dir"] = str(prev_sim_dir)
 
-    task_id = np.random.randint(1e9, 1e12)
+    task_id = random.randint(1e9, 1e12)
 
     if prev_sim_dir is None:
         param_seq = initialize.ParamSequence(config_dict)
@@ -718,7 +719,7 @@ def new_simulation(
 
 def resume_simulations(sim_dir: Path, method: Type[Simulations] = None) -> Simulations:
 
-    task_id = np.random.randint(1e9, 1e12)
+    task_id = random.randint(1e9, 1e12)
     config = io.load_toml(sim_dir / "initial_config.toml")
     io.set_data_folder(task_id, sim_dir)
     param_seq = initialize.RecoveryParamSequence(config, task_id)

src/scgenerator/plotting.py

@@ -30,6 +30,7 @@ def plot_setup(
     - an axis
     """
     out_path = defaults["name"] if out_path is None else out_path
+    out_path = Path(out_path)
     plot_name = out_path.name.replace(f".{file_type}", "")
     out_dir = out_path.resolve().parent