Add fit

src/backscattering_analyzer/__init__.py

@@ -1,7 +1,6 @@
 from __future__ import annotations
-from typing import Any
 from numpy.typing import NDArray
-from numpy import arange, float64
+from numpy import arange, float64, zeros, pi, sqrt
 
 
 def interpolate_abciss(signals: tuple[Signal, ...]) -> NDArray[float64]:
@@ -32,4 +31,71 @@ def interpolate(signals: tuple[Signal, ...]) -> tuple[Signal, ...]:
     return tuple(new_signals)
 
 
+def compute_light(
+    scatter_factor: list[float],
+    coupling: list[Signal],
+    movement: Signal,
+    wavelength: float,
+    power_in: float,
+    power_out: float,
+) -> Signal:
+    """
+    Compute the projection from a given coupling and movement
+    """
+    frequencies = interpolate_abciss(
+        (movement.sin().psd().sqrt(), coupling[0].abs())
+    )
+    parts = zeros(
+        (
+            len(coupling),
+            len(frequencies),
+        )
+    )
+
+    for index in range(2):
+        phase = (index + 1) * 4 * pi / wavelength
+
+        factor_n = (movement * phase).sin().psd().sqrt()
+        coupling_n = coupling[0].abs()
+        factor_d = (movement * phase).cos().psd().sqrt()
+        coupling_d = coupling[1].abs()
+
+        factor_n, coupling_n, factor_d, coupling_d = interpolate(
+            (factor_n, coupling_n, factor_d, coupling_d)
+        )
+
+        parts[index] = (
+            sqrt(scatter_factor[index])
+            * power_in
+            / power_out
+            * (coupling_n * factor_n + coupling_d * factor_d).y
+        )
+
+    return Signal(
+        frequencies,
+        sum(parts),
+        movement.settings,
+    )
+
+
+def opt_compute_light(
+    scatter_factor: float,
+    factor_n: Signal,
+    coupling_n: Signal,
+    factor_d: Signal,
+    coupling_d: Signal,
+    power_in: float,
+    power_out: float,
+) -> NDArray[float64]:
+    """
+    Optimized computing of light with pre-computed factor
+    """
+    return (
+        sqrt(scatter_factor)
+        * power_in
+        / power_out
+        * (coupling_n * factor_n + coupling_d * factor_d).y
+    )
+
+
 from backscattering_analyzer.signal import Signal  # no circular import

src/backscattering_analyzer/analyzer.py

@@ -2,20 +2,23 @@
 from sys import argv
 from backscattering_analyzer.settings import Settings
 from backscattering_analyzer.signal import Signal
-from backscattering_analyzer import interpolate, interpolate_abciss
+from backscattering_analyzer import (
-from numpy import loadtxt
+    compute_light,
+    opt_compute_light,
+    interpolate,
+)
+from numpy import loadtxt, logspace, where, zeros, argmin, intp, pi
 from scipy.io.matlab import loadmat
 
-# maths
-from numpy import zeros, pi, sqrt
-
 
 class Analyzer:
     """
     Utility class to study backscattering light in VIRGO
     """
 
-    def __init__(self, arguments: None | list[str] | str = None) -> None:
+    def __init__(
+        self, arguments: None | list[str] | str = None
+    ) -> None:
         if arguments is None:
             options = []
         elif isinstance(arguments, str):
@@ -135,47 +138,68 @@ class Analyzer:
 
     def compute_light(self) -> None:
         """
-        Compute psd of the computed projection with current bench
+        Compute projection with current bench excitation
-        excitation
         """
-        results = zeros(
+        self.projection = compute_light(
-            (
+            scatter_factor=self.settings.scattering_factor,
-                len(self.coupling),
+            coupling=self.coupling,
-                len(
+            movement=self.movement,
-                    interpolate_abciss(
+            wavelength=self.settings.wavelength,
-                        (self.movement.psd(), self.coupling[0].abs())
+            power_in=self.settings.power_in,
-                    )
+            power_out=self.settings.power_out,
-                ),
-            )
         )
 
-        # frequencies depends of psd result, which we do not have yet
+    def fit_scatter_factor(
-        frequencies = zeros(results.shape[1])
+        self, start: int, stop: int, number: int
-
+    ) -> tuple[intp, float]:
-        for index in range(len(self.coupling)):
+        """
-            phase = (index + 1) * 4 * pi / self.settings.wavelength
+        Find the best scatter factor (first order only) in the given
+        range
+        """
+        import matplotlib.pyplot as plt
+        factors = logspace(start, stop, number)
+        sums = zeros(number)
 
+        phase = 4 * pi / self.settings.wavelength
         factor_n = (self.movement * phase).sin().psd().sqrt()
         coupling_n = self.coupling[0].abs()
         factor_d = (self.movement * phase).cos().psd().sqrt()
         coupling_d = self.coupling[1].abs()
 
-            factor_n, coupling_n, factor_d, coupling_d = interpolate(
+        coupling_d.cut(5, 40) # cut signal between 5 and 40 Hz
-                (factor_n, coupling_n, factor_d, coupling_d)
+
+        factor_n, coupling_n, factor_d, coupling_d, data, reference = (
+            interpolate(
+                (
+                    factor_n,
+                    coupling_n,
+                    factor_d,
+                    coupling_d,
+                    self.data_signal.psd().sqrt(),
+                    self.reference_signal.psd().sqrt(),
+                )
+            )
         )
 
-            # no need to redefine it each time but simpler here
+        reference = reference.y
-            frequencies = factor_n.x
+        data = data.y
 
-            results[index] = (
+        for index in range(number):
-                sqrt(self.settings.scattering_factor[index])
+            self.settings.log("{}".format(index))
-                * self.settings.power_in
+            projection = opt_compute_light(
-                / self.settings.power_out
+                scatter_factor=factors[index],
-                * (coupling_n * factor_n + coupling_d * factor_d).y
+                factor_n=factor_n,
+                coupling_n=coupling_n,
+                factor_d=factor_d,
+                coupling_d=coupling_d,
+                power_in=self.settings.power_in,
+                power_out=self.settings.power_out,
             )
+            diff = abs(projection + reference - data)
+            _ = plt.loglog(projection + reference)
+            _ = plt.loglog(data)
+            _ = plt.show()
+            sums[index] = sum(diff)
+        min_index = argmin(sums)
 
-        self.projection = Signal(
+        return min_index, factors[min_index]
-            frequencies,
-            sum(results),
-            self.settings,
-        )

src/backscattering_analyzer/signal.py

@@ -15,7 +15,10 @@ class Signal:
     """
 
     def __init__(
-        self, x: NDArray[float64], value: NDArray[float64], settings: Settings
+        self,
+        x: NDArray[float64],
+        value: NDArray[float64],
+        settings: Settings,
     ) -> None:
         if x.shape != value.shape:
             raise Exception("x and y does not have the same dimension")