plasma is working but overflows everytime

src/scgenerator/data/materials.toml

@@ -78,6 +78,8 @@ source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field
 [helium]
 a = 0.00346
 b = 2.38e-5
+ionization_energy = 3.9393356074281e-18
+atomic_number = 2
 
 [helium.sellmeier]
 B = [2.16463842e-5, 2.10561127e-7, 4.7509272e-5]
@@ -131,6 +133,8 @@ source = "Shelton, D. P., & Rice, J. E. (1994). Measurements and calculations of
 [neon]
 a = 0.02135
 b = 1.709e-5
+ionization_energy = 3.45501365359425e-18
+atomic_number = 10
 
 [neon.sellmeier]
 B = [1281450000.0, 22048600000.0]
@@ -148,6 +152,8 @@ source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field
 [argon]
 a = 0.1355
 b = 3.201e-5
+ionization_energy = 2.5249661793774e-18
+atomic_number = 18
 
 [argon.sellmeier]
 B = [2501410000.0, 500283000.0, 52234300000.0]
@@ -203,6 +209,8 @@ source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field
 [krypton]
 a = 0.2349
 b = 3.978e-5
+ionization_energy = 2.2429831039374e-18
+atomic_number = 36
 
 [krypton.sellmeier]
 B = [2536370000.0, 2736490000.0, 62080200000.0]
@@ -220,6 +228,8 @@ source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field
 [xenon]
 a = 0.425
 b = 5.105e-5
+ionization_energy = 1.94342415157935
+atomic_number = 54
 
 [xenon.sellmeier]
 B = [3228690000.0, 3553930000.0, 60676400000.0]

src/scgenerator/evaluator.py

@@ -375,10 +375,10 @@ default_rules: list[Rule] = [
     Rule("n_op", operators.MarcatiliRefractiveIndex),
     Rule("n_op", operators.MarcatiliAdjustedRefractiveIndex),
     Rule("n_op", operators.HasanRefractiveIndex),
+    Rule("gas_op", operators.ConstantGas),
     Rule("raman_op", operators.NoRaman, priorities=-1),
     Rule("loss_op", operators.NoLoss, priorities=-1),
     Rule("conserved_quantity", operators.NoConservedQuantity, priorities=-1),
-    Rule("conserved_quantity", operators.conserved_quantity),
 ]
 
 envelope_rules = default_rules + [
@@ -419,6 +419,7 @@ envelope_rules = default_rules + [
     Rule("dispersion_op", operators.ConstantPolyDispersion),
     Rule("dispersion_op", operators.DirectDispersion),
     Rule("linear_operator", operators.EnvelopeLinearOperator),
+    Rule("conserved_quantity", operators.conserved_quantity),
 ]
 
 full_field_rules = default_rules + [

src/scgenerator/operators.py

@@ -801,6 +801,7 @@ class ConstantGamma(AbstractGamma):
 class Plasma(Operator):
     mat_plasma: plasma.Plasma
     gas_op: AbstractGas
+    ionization_fraction = 0.0
 
     def __init__(self, dt: float, gas_op: AbstractGas):
         self.gas_op = gas_op
@@ -812,7 +813,12 @@ class Plasma(Operator):
 
     def __call__(self, state: CurrentState) -> np.ndarray:
         N0 = self.gas_op.number_density(state)
-        return self.mat_plasma(state.field, N0)
+        plasma_info = self.mat_plasma(state.field, N0)
+        self.ionization_fraction = plasma_info.electron_density[-1] / N0
+        return plasma_info.polarization
+
+    def values(self) -> dict[str, float]:
+        return dict(ionization_fraction=self.ionization_fraction)
 
 
 class NoPlasma(NoOpTime, Plasma):

src/scgenerator/physics/plasma.py

@@ -9,10 +9,7 @@ from ..math import expm1_int, cumulative_simpson
 @dataclass
 class PlasmaInfo:
     electron_density: np.ndarray
-    dn_dt: np.ndarray
     polarization: np.ndarray
-    loss: np.ndarray
-    debug: np.ndarray
 
 
 class IonizationRate:
@@ -66,12 +63,11 @@ class Plasma:
         rate = self.rate(field_abs)
         electron_density = free_electron_density(rate, self.dt, N0)
         dn_dt = (N0 - electron_density) * rate
-        out = self.dt * cumulative_simpson(
+        polarization = self.dt * cumulative_simpson(
-            dn_dt * self.Ip / (field + delta)
+            dn_dt * self.Ip / np.where(field_abs < delta, 1e-14, field)
             + e ** 2 / me * self.dt * cumulative_simpson(electron_density * field)
         )
-        loss = cumulative_simpson(dn_dt * self.Ip / (field + delta)) * self.dt
+        return PlasmaInfo(electron_density, polarization)
-        return PlasmaInfo(electron_density, dn_dt, out, loss, loss)
 
 
 def adiabadicity(w: np.ndarray, I: float, field: np.ndarray) -> np.ndarray:

src/scgenerator/physics/simulate.py

@@ -164,12 +164,19 @@ class RK4IP:
         state = self.init_state.copy()
         yield len(self.stored_spectra) - 1, state
 
-        integrator_args = [
+        if self.params.adapt_step_size:
-            self.params.compute(a) for a in solver.Integrator.factory_args() if a != "init_state"
+            integrator_args = [
-        ]
+                self.params.compute(a)
-        integrator = solver.Integrator.create(
+                for a in solver.Integrator.factory_args()
-            self.params.integration_scheme, state, *integrator_args
+                if a != "init_state"
-        )
+            ]
+            integrator = solver.Integrator.create(
+                self.params.integration_scheme, state, *integrator_args
+            )
+        else:
+            integrator = solver.ConstantStepIntegrator(
+                state, self.params.linear_operator, self.params.nonlinear_operator
+            )
 
         for state in integrator:
 

src/scgenerator/scripts/__init__.py

@@ -261,7 +261,7 @@ def finish_plot(fig: plt.Figure, legend_axes: plt.Axes, all_labels: list[str], p
 
 def plot_helper(config_path: Path) -> Iterable[tuple[dict, list[str], Parameters]]:
     cc = cycler(color=[f"C{i}" for i in range(10)]) * cycler(ls=["-", "--"])
-    for style, (descriptor, params) in zip(cc, Configuration(config_path)):
+    for style, (descriptor, params), _ in zip(cc, Configuration(config_path), range(20)):
         yield style, descriptor.branch.formatted_descriptor(), params
 
 

src/scgenerator/solver.py

@@ -130,7 +130,7 @@ class ConstantStepIntegrator(Integrator, scheme="constant"):
         linear_operator: LinearOperator,
         nonlinear_operator: NonLinearOperator,
     ):
-        super().__init__(init_state, linear_operator, nonlinear_operator, 0.0)
+        super().__init__(init_state, linear_operator, nonlinear_operator)
 
     def __iter__(self) -> Iterator[CurrentState]:
         while True:
@@ -156,6 +156,7 @@ class AdaptiveIntegrator(Integrator):
     target_error: float
     current_error: float = 0.0
     steps_rejected: float = 0
+    min_step_size = 0.0
 
     def __init__(
         self,
@@ -190,7 +191,8 @@ class AdaptiveIntegrator(Integrator):
         else:
             next_h_factor = 2.0
         h_next_step = h * next_h_factor
-        accepted = self.current_error <= 2 * self.target_error
+        accepted = self.current_error <= 2 * self.target_error or h_next_step < self.min_step_size
+        h_next_step = max(h_next_step, self.min_step_size)
         if not accepted:
             self.steps_rejected += 1
             self.logger.info(

src/scgenerator/spectra.py

@@ -17,6 +17,7 @@ from .plotting import (
     mean_values_plot,
     propagation_plot,
     single_position_plot,
+    transform_1D_values,
     transform_2D_propagation,
 )
 from .utils import load_spectrum, simulations_list, load_toml
@@ -246,6 +247,19 @@ class SimulationSeries:
         vals = self.retrieve_plot_values(plot_range, z_pos, sim_ind)
         return single_position_plot(vals, plot_range, self.params, ax, **kwargs)
 
+    def plot_values_1D(
+        self,
+        left: float,
+        right: float,
+        unit: Union[Callable[[float], float], str],
+        z_pos: int,
+        sim_ind: int = 0,
+        **kwargs,
+    ) -> tuple[np.ndarray, np.ndarray]:
+        plot_range = PlotRange(left, right, unit)
+        vals = self.retrieve_plot_values(plot_range, z_pos, sim_ind)
+        return transform_1D_values(vals, plot_range, self.params, **kwargs)
+
     def plot_mean(
         self,
         left: float,

src/scgenerator/utils.py

@@ -11,7 +11,7 @@ import os
 import re
 from collections import defaultdict
 from dataclasses import dataclass
-from functools import cache
+from functools import cache, lru_cache
 from pathlib import Path
 from string import printable as str_printable
 from typing import Any, Callable, MutableMapping, Sequence, TypeVar, Set, Union
@@ -21,7 +21,7 @@ import pkg_resources as pkg
 import tomli
 import tomli_w
 
-from .const import PARAM_FN, PARAM_SEPARATOR, SPEC1_FN, Z_FN, ROOT_PARAMETERS
+from .const import PARAM_FN, PARAM_SEPARATOR, SPEC1_FN, SPECN_FN1, Z_FN, ROOT_PARAMETERS
 from .logger import get_logger
 from .errors import DuplicateParameterError
 
@@ -134,7 +134,7 @@ def load_previous_spectrum(prev_data_dir: str) -> np.ndarray:
     return load_spectrum(prev_data_dir / SPEC1_FN.format(num))
 
 
-@cache
+@lru_cache(20000)
 def load_spectrum(file: os.PathLike) -> np.ndarray:
     return np.load(file)
 
@@ -535,7 +535,7 @@ def simulations_list(path: os.PathLike) -> list[Path]:
     """
     paths: list[Path] = []
     for pwd, _, files in os.walk(path):
-        if PARAM_FN in files:
+        if PARAM_FN in files and SPEC1_FN.format(0) in files:
             paths.append(Path(pwd))
     paths.sort(key=branch_id)
     return [p for p in paths if p.parent.name == paths[-1].parent.name]