Working Full field

src/scgenerator/__init__.py

@@ -4,7 +4,7 @@ from .evaluator import Evaluator
 from .legacy import convert_sim_folder
 from .math import abs2, argclosest, normalized, span
 from .parameter import Configuration, Parameters
-from .physics import fiber, materials, pulse, simulate, units
+from .physics import fiber, materials, pulse, simulate, units, plasma
 from .physics.simulate import RK4IP, parallel_RK4IP, run_simulation
 from .physics.units import PlotRange
 from .plotting import (
@@ -18,5 +18,5 @@ from .plotting import (
     transform_mean_values,
 )
 from .spectra import SimulationSeries, Spectrum
-from .utils import Paths, _open_config, open_single_config
+from .utils import Paths, _open_config, open_single_config, simulations_list
 from .variationer import DescriptorDict, VariationDescriptor, Variationer, VariationSpecsError

src/scgenerator/const.py

@@ -2,6 +2,12 @@ __version__ = "0.2.4dev"
 
 from typing import Any
 
+ONE_2 = 1 / 2
+ONE_3 = 1 / 3
+ONE_FOURTH = 1 / 4
+ONE_FIFTH = 1 / 5
+ONE_6 = 1 / 6
+
 
 def pbar_format(worker_id: int) -> dict[str, Any]:
     if worker_id == 0:

src/scgenerator/math.py

@@ -1,7 +1,6 @@
 from typing import Union
 
 import numpy as np
-from scipy.integrate import cumulative_trapezoid
 from scipy.special import jn_zeros
 import numba
 
@@ -11,10 +10,9 @@ pi = np.pi
 c = 299792458.0
 
 
-def inverse_integral_exponential(y: np.ndarray, dx: float) -> np.ndarray:
-    """evaluates exp(-func(y)) from x=-inf to x"""
-    # return np.exp(-cumulative_trapezoid(y, dx=dx, initial=0))
-    return np.exp(-cumulative_simpson(y) * dx)
+def expm1_int(y: np.ndarray, dx: float) -> np.ndarray:
+    """evaluates 1 - exp( -∫func(y(x))dx ) from x=-inf to x"""
+    return -np.expm1(-cumulative_simpson(y) * dx)
 
 
 def span(*vec):

src/scgenerator/operators.py

@@ -69,6 +69,16 @@ class CurrentState:
 
 
 class Operator(ABC):
+    def values(self) -> dict[str, float]:
+        return {}
+
+    def get_values(self) -> dict[str, float]:
+        out = self.values()
+        for operator in self.__dict__.values():
+            if isinstance(operator, Operator):
+                out |= operator.get_values()
+        return out
+
     def __repr__(self) -> str:
         value_pair_list = list(self.__dict__.items())
         if len(value_pair_list) == 0:

src/scgenerator/physics/plasma.py

@@ -3,7 +3,7 @@ import numpy as np
 import scipy.special
 
 from .units import e, hbar, me
-from ..math import inverse_integral_exponential, cumulative_simpson
+from ..math import expm1_int, cumulative_simpson
 
 
 @dataclass
@@ -12,7 +12,7 @@ class PlasmaInfo:
     dn_dt: np.ndarray
     polarization: np.ndarray
     loss: np.ndarray
-    phase_effect: np.ndarray
+    debug: np.ndarray
 
 
 class IonizationRate:
@@ -64,17 +64,15 @@ class Plasma:
         field_abs = np.abs(field)
         delta = 1e-14 * field_abs.max()
         rate = self.rate(field_abs)
-        exp_int = inverse_integral_exponential(rate, self.dt)
-        electron_density = N0 * (1 - exp_int)
-        dn_dt = N0 * rate * exp_int
+        exp_int = expm1_int(rate, self.dt)
+        electron_density = N0 * exp_int
+        dn_dt = (N0 - electron_density) * rate
         out = self.dt * cumulative_simpson(
             dn_dt * self.Ip / (field + delta)
             + e ** 2 / me * self.dt * cumulative_simpson(electron_density * field)
         )
         loss = cumulative_simpson(dn_dt * self.Ip / (field + delta)) * self.dt
-        phase_effect = e ** 2 / me * self.dt * cumulative_simpson(electron_density * field)
-        phase_effect = exp_int
-        return PlasmaInfo(electron_density, dn_dt, out, loss, phase_effect)
+        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

@@ -13,6 +13,7 @@ from ..logger import get_logger
 from ..operators import CurrentState
 from ..parameter import Configuration, Parameters
 from ..pbar import PBars, ProgressBarActor, progress_worker
+from ..const import ONE_2, ONE_3, ONE_6
 
 try:
     import ray
@@ -248,14 +249,16 @@ class RK4IP:
         while not keep:
             h = h_next_step
 
-            expD = np.exp(h / 2 * self.params.linear_operator(self.state))
+            expD = np.exp(h * ONE_2 * self.params.linear_operator(self.state))
 
             A_I = expD * self.state.spectrum
             k1 = expD * (h * self.params.nonlinear_operator(self.state))
-            k2 = h * self.params.nonlinear_operator(self.state.replace(A_I + k1 / 2))
-            k3 = h * self.params.nonlinear_operator(self.state.replace(A_I + k2 / 2))
+            k2 = h * self.params.nonlinear_operator(self.state.replace(A_I + k1 * ONE_2))
+            k3 = h * self.params.nonlinear_operator(self.state.replace(A_I + k2 * ONE_2))
             k4 = h * self.params.nonlinear_operator(self.state.replace(expD * (A_I + k3)))
-            new_state = self.state.replace(expD * (A_I + k1 / 6 + k2 / 3 + k3 / 3) + k4 / 6)
+            new_state = self.state.replace(
+                expD * (A_I + k1 * ONE_6 + k2 * ONE_3 + k3 * ONE_3) + k4 * ONE_6
+            )
 
             self.cons_qty[step] = self.params.conserved_quantity(new_state)
             if self.params.adapt_step_size:
@@ -266,8 +269,8 @@ class RK4IP:
                     progress_str = f"step {step} rejected with h = {h:.4e}, doing over"
                     self.logger.debug(progress_str)
                     keep = False
-                    h_next_step = h / 2
-                elif cons_qty_change_ok < curr_p_change <= 2 * cons_qty_change_ok:
+                    h_next_step = h * ONE_2
+                elif cons_qty_change_ok < curr_p_change <= 2.0 * cons_qty_change_ok:
                     keep = True
                     h_next_step = h / self.size_fac
                 elif curr_p_change < 0.1 * cons_qty_change_ok: