removed phase term from photoionization

src/scgenerator/evaluator.py

@@ -386,7 +386,6 @@ default_rules: list[Rule] = [
     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),
 ]
@@ -422,6 +421,7 @@ envelope_rules = default_rules + [
     Rule("spm_op", operators.NoEnvelopeSPM, priorities=-1),
     Rule("spm_op", operators.EnvelopeSPM),
     Rule("raman_op", operators.EnvelopeRaman),
+    Rule("raman_op", operators.NoEnvelopeRaman, priorities=-1),
     Rule("nonlinear_operator", operators.EnvelopeNonLinearOperator),
     Rule("loss_op", operators.CustomLoss, priorities=3),
     Rule("loss_op", operators.CapillaryLoss, priorities=2, conditions=dict(loss="capillary")),
@@ -454,5 +454,6 @@ full_field_rules = default_rules + [
     ),
     Rule("plasma_op", operators.Plasma, conditions=dict(photoionization=True)),
     Rule("plasma_op", operators.NoPlasma, priorities=-1),
+    Rule("raman_op", operators.NoFullFieldRaman, priorities=-1),
     Rule("nonlinear_operator", operators.FullFieldNonLinearOperator),
 ]

src/scgenerator/operators.py

@@ -630,7 +630,7 @@ class CustomLoss(ConstantLoss):
 
 
 class AbstractRaman(Operator):
-    f_r: float = 0.0
+    fraction: float = 0.0
 
     @abstractmethod
     def __call__(self, state: CurrentState) -> np.ndarray:
@@ -647,27 +647,34 @@ class AbstractRaman(Operator):
         """
 
 
-class NoRaman(NoOpTime, AbstractRaman):
+class NoEnvelopeRaman(NoOpTime, AbstractRaman):
+    pass
+
+
+class NoFullFieldRaman(NoOpFreq, AbstractRaman):
     pass
 
 
 class EnvelopeRaman(AbstractRaman):
     def __init__(self, raman_type: str, t: np.ndarray):
         self.hr_w = fiber.delayed_raman_w(t, raman_type)
-        self.f_r = 0.245 if raman_type == "agrawal" else 0.18
+        self.fraction = 0.245 if raman_type == "agrawal" else 0.18
 
     def __call__(self, state: CurrentState) -> np.ndarray:
-        return self.f_r * np.fft.ifft(self.hr_w * np.fft.fft(state.field2))
+        return self.fraction * np.fft.ifft(self.hr_w * np.fft.fft(state.field2))
 
 
 class FullFieldRaman(AbstractRaman):
-    def __init__(self, raman_type: str, t: np.ndarray, chi3: float):
+    def __init__(self, raman_type: str, t: np.ndarray, w: np.ndarray, chi3: float):
         self.hr_w = fiber.delayed_raman_w(t, raman_type)
-        self.f_r = 0.245 if raman_type == "agrawal" else 0.18
+        self.fraction = 0.245 if raman_type == "agrawal" else 0.18
-        self.multiplier = units.epsilon0 * chi3 * self.f_r
+        self.factor_in = units.epsilon0 * chi3 * self.fraction
+        self.factor_out = 1j * w ** 2 / (2.0 * units.c ** 2 * units.epsilon0)
 
     def __call__(self, state: CurrentState) -> np.ndarray:
-        return self.multiplier * np.fft.ifft(np.fft.fft(state.field2) * self.hr_w)
+        return self.factor_out * np.fft.rfft(
+            self.factor_in * state.field * np.fft.irfft(self.hr_w * np.fft.rfft(state.field2))
+        )
 
 
 ##################################################
@@ -703,19 +710,20 @@ class NoFullFieldSPM(NoOpTime, AbstractSPM):
 
 class EnvelopeSPM(AbstractSPM):
     def __init__(self, raman_op: AbstractRaman):
-        self.fraction = 1 - raman_op.f_r
+        self.fraction = 1 - raman_op.fraction
 
     def __call__(self, state: CurrentState) -> np.ndarray:
         return self.fraction * state.field2
 
 
 class FullFieldSPM(AbstractSPM):
-    def __init__(self, raman_op: AbstractRaman, chi3: float):
+    def __init__(self, raman_op: AbstractRaman, w: np.ndarray, chi3: float):
-        self.fraction = 1 - raman_op.f_r
+        self.fraction = 1 - raman_op.fraction
-        self.factor = self.fraction * chi3 * units.epsilon0
+        self.factor_out = 1j * w ** 2 / (2.0 * units.c ** 2 * units.epsilon0)
+        self.factor_in = self.fraction * chi3 * units.epsilon0
 
     def __call__(self, state: CurrentState) -> np.ndarray:
-        return self.factor * state.field2 * state.field
+        return self.factor_out * np.fft.rfft(self.factor_in * state.field2 * state.field)
 
 
 ##################################################
@@ -804,19 +812,20 @@ class Plasma(Operator):
     gas_op: AbstractGas
     ionization_fraction = 0.0
 
-    def __init__(self, dt: float, gas_op: AbstractGas):
+    def __init__(self, dt: float, w: np.ndarray, gas_op: AbstractGas):
         self.gas_op = gas_op
         self.mat_plasma = plasma.Plasma(
             dt,
             self.gas_op.material_dico["ionization_energy"],
             plasma.IonizationRateADK(self.gas_op.material_dico["ionization_energy"]),
         )
+        self.factor_out = -w / (2.0 * units.c ** 2 * units.epsilon0)
 
     def __call__(self, state: CurrentState) -> np.ndarray:
         N0 = self.gas_op.number_density(state)
         plasma_info = self.mat_plasma(state.field, N0)
         self.ionization_fraction = plasma_info.electron_density[-1] / N0
-        return plasma_info.polarization
+        return self.factor_out * np.fft.rfft(plasma_info.dp_dt)
 
     def values(self) -> dict[str, float]:
         return dict(ionization_fraction=self.ionization_fraction)
@@ -902,7 +911,7 @@ def conserved_quantity(
     if not adapt_step_size:
         return NoConservedQuantity()
     logger = get_logger(__name__)
-    raman = not isinstance(raman_op, NoRaman)
+    raman = not isinstance(raman_op, NoEnvelopeRaman)
     loss = not isinstance(loss_op, NoLoss)
     if raman and loss:
         logger.debug("Conserved quantity : photon number with loss")
@@ -1023,4 +1032,4 @@ class FullFieldNonLinearOperator(NonLinearOperator):
 
     def __call__(self, state: CurrentState) -> np.ndarray:
         total_nonlinear = self.spm_op(state) + self.raman_op(state) + self.plasma_op(state)
-        return self.factor / self.beta_op(state) * np.fft.rfft(total_nonlinear)
+        return total_nonlinear / self.beta_op(state)

src/scgenerator/physics/plasma.py

@@ -12,7 +12,7 @@ T_a = TypeVar("T_a", np.floating, np.ndarray)
 @dataclass
 class PlasmaInfo:
     electron_density: np.ndarray
-    polarization: np.ndarray
+    dp_dt: np.ndarray
     rate: np.ndarray
     dn_dt: np.ndarray
     debug: np.ndarray
@@ -79,18 +79,15 @@ class Plasma:
         rate = self.rate(field_abs)
         electron_density = free_electron_density(rate, self.dt, N0)
         dn_dt: np.ndarray = (N0 - electron_density) * rate
-        integrand = np.zeros_like(field)
-        integrand[nzi] = dn_dt[nzi] * self.ionization_energy / field[nzi]
 
-        energy_loss = self.dt * cumulative_simpson(integrand)
+        loss_term = np.zeros_like(field)
-        added_phase = (
+        loss_term[nzi] = dn_dt[nzi] * self.ionization_energy / field[nzi]
-            self.dt ** 2
+
-            * e ** 2
+        phase_term = self.dt * e ** 2 / me * cumulative_simpson(electron_density * field)
-            / me
+
-            * cumulative_simpson(cumulative_simpson(electron_density * field))
+        # polarization = -loss_integrated + phase_integrated
-        )
+        dp_dt = loss_term
-        polarization = energy_loss + added_phase
+        return PlasmaInfo(electron_density, dp_dt, rate, dn_dt, phase_term)
-        return PlasmaInfo(electron_density, polarization, rate, dn_dt, (energy_loss, added_phase))
 
 
 def adiabadicity(w: np.ndarray, I: float, field: np.ndarray) -> np.ndarray:

src/scgenerator/solver.py

@@ -19,6 +19,11 @@ from .operators import (
 from .utils import get_arg_names
 import warnings
 
+
+class IntegratorError(Exception):
+    pass
+
+
 # warnings.filterwarnings("error")