diff --git a/pyproject.toml b/pyproject.toml index 3646bec..9a1bd1a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "backscattering_analyzer" -version = "0.0.1" +version = "0.2.0" authors = [ { name="linarphy", email="linarphy@linarphy.net" }, ] @@ -12,13 +12,12 @@ classifiers = [ "License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)", "Intended Audience :: Science/Research", "Operating System :: OS Independent", - "Development Status :: 2 - Pre-Alpha", + "Development Status :: 4 - Beta", "Topic :: Scientific/Engineering", ] dependencies = [ "numpy", "scipy", - "rich", ] [project.urls] diff --git a/src/backscattering_analyzer/__init__.py b/src/backscattering_analyzer/__init__.py index c27c8ae..2a72377 100644 --- a/src/backscattering_analyzer/__init__.py +++ b/src/backscattering_analyzer/__init__.py @@ -1,112 +1,34 @@ -from science_signal.signal import Signal -from science_signal import interpolate, interpolate_abciss -from numpy.typing import NDArray -from numpy import float64, zeros, pi, sqrt +from logging import getLogger + +from gettext import install + +install(__name__) + +logger = getLogger(__name__) -def compute_light( - scatter_factor: list[float], - coupling: list[Signal], - movement: Signal, - wavelength: float, - power_in: float, - power_out: float, -) -> Signal: +def show_projection(experiment: "Experiment") -> None: """ - Compute the projection from a given coupling and movement + Show projection data with matplotlib """ - frequencies = interpolate_abciss( - movement.sin().psd().sqrt(), abs(coupling[0]) - ) - parts = zeros( - ( - len(coupling), - len(frequencies), - ) - ) + logger.debug(_("showing experiment result")) + from matplotlib.pyplot import loglog, show, legend, close - for index in range(2): - phase = (index + 1) * 4 * pi / wavelength - - factor_n = (movement * phase).sin().psd().sqrt() - coupling_n = abs(coupling[0]) - factor_d = (movement * phase).cos().psd().sqrt() - coupling_d = abs(coupling[1]) - - factor_n, coupling_n, factor_d, coupling_d = interpolate( - factor_n, coupling_n, factor_d, coupling_d - ) - - parts[index] = opt_compute_light( - scatter_factor[index], - factor_n, - coupling_n, - factor_d, - coupling_d, - power_in, - power_out, - phase, - ) - - return Signal( - frequencies, - sum(parts), - ) - - -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, - phase: float, -) -> NDArray[float64]: - """ - Optimized computing of light with pre-computed factor - """ - return ( - sqrt(scatter_factor) - / phase - * power_in - / power_out - * (coupling_n * factor_n + coupling_d * factor_d).y - ) - - -def fit_compute_light( - scatter_factor: float, - factor_n: Signal, - coupling_n: Signal, - factor_d: Signal, - coupling_d: Signal, - power_in: float, - power_out: float, - data: Signal, - reference: Signal, - phase: float, -) -> float: - """ - Scalar function used to find the right scattering factor - """ - return sum( - abs( - Signal( - factor_n.x, - opt_compute_light( - scatter_factor=scatter_factor, - factor_n=factor_n, - coupling_n=coupling_n, - factor_d=factor_d, - coupling_d=coupling_d, - power_in=power_in, - power_out=power_out, - phase=phase, - ), - ) - + reference - - data - ).y - ) + excited = experiment.signals["excited"].psd().sqrt() + reference = experiment.signals["reference"].psd().sqrt() + loglog( + experiment.projection.x, + experiment.projection.y, + label="projection", + ) # type: ignore[ReportUnusedCallResult] + loglog(excited.x, excited.y, label="excited bench") # type: ignore[ReportUnusedCallResult] + loglog(reference.x, reference.y, label="reference bench") # type: ignore[ReportUnusedCallResult] + loglog( + (experiment.projection + reference).x, + (experiment.projection + reference).y, + label="sum reference + excited", + ) # type: ignore[ReportUnusedCallResult] + legend() # type: ignore[ReportUnusedCallResult] + show() + close() + logger.debug(_("experiment result showed")) diff --git a/src/backscattering_analyzer/analyzer.py b/src/backscattering_analyzer/analyzer.py deleted file mode 100644 index 031b8d0..0000000 --- a/src/backscattering_analyzer/analyzer.py +++ /dev/null @@ -1,225 +0,0 @@ -from sys import argv - -from backscattering_analyzer.settings import Settings -from science_signal.signal import Signal -from science_signal import interpolate -from backscattering_analyzer import ( - compute_light, - fit_compute_light, - opt_compute_light, -) -from numpy import argmin, loadtxt, logspace, pi, array -from scipy.io.matlab import loadmat - - -class Analyzer: - """ - Utility class to study backscattering light in VIRGO - """ - - def __init__( - self, arguments: None | list[str] | str = None - ) -> None: - if arguments is None: - options = [] - elif isinstance(arguments, str): - options = arguments.split(" ") - else: - options = arguments - self.settings = Settings(argv[1:] + options) - self.load() - self.process_movement() - self.compute_light() - - def load(self): - """ - Load all data - """ - self.settings.log("Load all data") - try: - self.load_bench() - self.load_mirror() - self.load_data() - self.load_reference() - self.load_coupling() - except Exception: - raise Exception("An error occured during data loading") - - def load_bench(self): - """ - Load bench movement - """ - file = self.settings.bench_file() - self.settings.log("loading bench movement: {}".format(self.settings.bench_file())) - try: - data = loadtxt(file).T - self.bench_movement = Signal(data[0], data[1]) * 1e-6 # um - except OSError: - raise Exception("{file} does not exist".format(file=file)) - - def load_mirror(self): - """ - Load mirror movement - """ - file = self.settings.mirror_file() - self.settings.log("loading mirror movement: {}".format(self.settings.mirror_file())) - try: - data = loadtxt(file).T - self.mirror_movement = Signal(data[0], data[1]) * 1e-6 # um - except OSError: - raise Exception("{file} does not exist".format(file=file)) - - def load_data(self): - """ - Load excited h(t) - """ - file = self.settings.data_file() - self.settings.log("loading excited h(t): {}".format(self.settings.data_file())) - try: - data = loadtxt(file).T - self.data_signal = Signal(data[0], data[1]) - except OSError: - raise Exception("{file} does not exist".format(file=file)) - - def load_reference(self): - """ - Load reference h(t) - """ - file = self.settings.reference_file() - self.settings.log("loading reference h(t): {}".format(self.settings.reference_file())) - try: - data = loadtxt(file).T - self.reference_signal = Signal(data[0], data[1]) - except OSError: - raise Exception("{file} does not exist".format(file=file)) - - def load_coupling(self): - """ - Load modelisation coupling data - """ - file = self.settings.modelisation_file() - self.settings.log("loading matlab modelisation data: {}".format(self.settings.modelisation_file())) - try: - self.modelisation = loadmat(file) - coupling_values = self.modelisation[ - self.settings.coupling_name() - ] - self.coupling = [ - Signal( - self.modelisation["freq"][0], - coupling, - ) - for coupling in coupling_values - ] - except OSError: - raise Exception("{file} does not exist".format(file=file)) - - def process_movement(self): - """ - Clean and compute relative movement between the bench and the - miror - """ - self.settings.log("computing relative movement") - - self.movement = ( - ( - self.bench_movement * self.settings.calib_bench - - self.mirror_movement * self.settings.calib_mirror - ) - .detrend("linear") - .filter(end=5 * self.settings.vibration_frequency) - ) - - def compute_light(self) -> None: - """ - Compute projection with current bench excitation - """ - self.settings.log("computing backscatterd light projection") - self.projection = compute_light( - scatter_factor=self.settings.scattering_factor, - coupling=self.coupling, - movement=self.movement, - wavelength=self.settings.wavelength, - power_in=self.settings.power_in, - power_out=self.settings.power_out, - ) - - def fit_scatter_factor(self, guess: None | float = None) -> float: - """ - Find the best scatter factor (first order only) in the given - range - """ - if guess is None: - guess = self.settings.scattering_factor[0] - - import matplotlib.pyplot as plt - - phase = 4 * pi / self.settings.wavelength - factor_n = (self.movement * phase).sin().psd().sqrt() - coupling_n = abs(self.coupling[0]) - factor_d = (self.movement * phase).cos().psd().sqrt() - coupling_d = abs(self.coupling[1]) - - coupling_d = coupling_d.cut( - 10, - 200, - ) # cut signal between 10 and 200 Hz (updated to 11-15 to test) - - 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(), - ) - ) - - x = logspace(-15, -3, 1000) - y = array( - [ - fit_compute_light( - x[i], - factor_n, - coupling_n, - factor_d, - coupling_d, - self.settings.power_in, - self.settings.power_out, - data, - reference, - phase, - ) - for i in range(len(x)) - ] - ) - - factor: float = x[argmin(y)] - - projection = Signal( - factor_n.x, - opt_compute_light( - scatter_factor=factor, - 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, - phase=phase, - ), - ) - _ = plt.loglog(x, y) - _ = plt.show() - - _ = plt.loglog(projection.x, projection.y, label="projection") - _ = plt.loglog(reference.x, reference.y, label="référence") - _ = plt.loglog(data.x, data.y, label="data") - _ = plt.loglog( - reference.x, projection.y + reference.y, label="somme" - ) - _ = plt.legend() - _ = plt.show() - - return factor diff --git a/src/backscattering_analyzer/data.py b/src/backscattering_analyzer/data.py new file mode 100644 index 0000000..045fc55 --- /dev/null +++ b/src/backscattering_analyzer/data.py @@ -0,0 +1,161 @@ +from pathlib import Path +from numpy.lib.npyio import loadtxt +from science_signal.signal import Signal +from scipy.io.matlab import loadmat +from backscattering_analyzer import logger + + +def load_signals( + folder: Path, bench: str, date: str +) -> dict[str, Signal]: + """ + Load excited and reference signal of an experiment + """ + logger.info( + _("loading signals for {bench} in {date}").format( + bench=bench, date=date + ) + ) + excited_file = folder / "{bench}_{date}_dat.csv".format( + bench=bench, + date=date, + ) + reference_file = folder / "{bench}_{date}_ref.csv".format( + bench=bench, + date=date, + ) + try: + excited_data = loadtxt(excited_file).T + reference_data = loadtxt(reference_file).T + logger.info( + _( + "files successfully loaded:\n{excited}\n{reference}" + ).format(excited=excited_file, reference=reference_file) + ) + except OSError as e: + logger.critical( + _("some file cannot be read: {error}").format(error=e) + ) + raise e + return { + "excited": Signal(excited_data[0], excited_data[1]), + "reference": Signal(reference_data[0], reference_data[1]), + } + + +def load_movements( + folder: Path, bench: str, date: str +) -> dict[str, Signal]: + """ + Load excited movement of the bench and the nearest mirror of an + experiment + """ + logger.info( + _("loading movements for {bench} in {date}").format( + bench=bench, date=date + ) + ) + bench_file = folder / "{bench}_{date}_ben.csv".format( + bench=bench, + date=date, + ) + mirror_file = folder / "{bench}_{date}_mir.csv".format( + bench=bench, + date=date, + ) + try: + bench_movement = loadtxt(bench_file).T + mirror_movement = loadtxt(mirror_file).T + bench_movement[1] = bench_movement[1] * 1e-6 # um -> m + mirror_movement[1] = mirror_movement[1] * 1e-6 # um -> m + logger.info( + _("files successfully loaded:\n{bench}\n{mirror}").format( + bench=bench_file, mirror=mirror_file + ) + ) + except OSError as e: + logger.critical( + _("some file cannot be read: {error}").format(error=e) + ) + raise e + return { + "bench": Signal(bench_movement[0], bench_movement[1]), + "mirror": Signal(mirror_movement[0], mirror_movement[1]), + } + + +def load_coupling( + folder: Path, + bench: str, + date: str, + file: str, + powers: dict[str, float], +) -> list[Signal]: + """ + Load and correct coupling date from modelisation + """ + logger.info( + _("loading coupling for {bench} in {date}").format( + bench=bench, date=date + ) + ) + modelisation_file = folder / file + try: + modelisation_data = loadmat(modelisation_file) + logger.info( + _("file successfully loaded:\n{file}").format( + file=modelisation_file + ) + ) + except OSError: + logger.critical( + _("cannot read {file}").format(file=modelisation_file) + ) + raise Exception( + _("cannot read {file}").format(file=modelisation_file) + ) + if bench == "SDB2": + modelisation_name = "SDB1coupling" + elif bench == "SIB2": + modelisation_name = "SIB1coupling" + else: + modelisation_name = "{bench}coupling".format(bench=bench) + + try: + modelisation_signal = [ + Signal( + modelisation_data["freq"][0], + abs(coupling), + ) + for coupling in modelisation_data[modelisation_name] + ] + except KeyError: + logger.critical( + _("{name} is not defined in the modelisation data").format( + name=modelisation_name + ) + ) + raise Exception( + _("{name} is not defined in the modelisation data").format( + name=modelisation_name + ) + ) + + if bench in ["SDB1", "SDB2"]: + modelisation_signal[0] = ( + modelisation_signal[0] + * (powers["laser"] / 40) + * (powers["dark_fringe"] / 5e-6) ** (1 / 2) + ) # radiation + modelisation_signal[1] = modelisation_signal[1] * ( + powers["dark_fringe"] / 5e-6 + ) ** (1 / 2) # phase + elif bench in ["SWEB", "SNEB"]: + modelisation_signal[1] = ( + modelisation_signal[1] * powers["laser"] / 40 + ) # radiation + elif bench in ["SIB2", "SPRB"]: + logger.info( + _("cannot fix projection power for SIB2 and SPRB coupling") + ) + return modelisation_signal diff --git a/src/backscattering_analyzer/experiment.py b/src/backscattering_analyzer/experiment.py new file mode 100644 index 0000000..799b9f9 --- /dev/null +++ b/src/backscattering_analyzer/experiment.py @@ -0,0 +1,379 @@ +from pathlib import Path +from tomllib import load + +from numpy import argmin, float64, pi +from numpy.core.function_base import logspace +from numpy.core.multiarray import array +from numpy.typing import NDArray +from science_signal import interpolate +from science_signal.signal import Signal +from backscattering_analyzer import logger +from backscattering_analyzer.data import ( + load_coupling, + load_movements, + load_signals, +) +from backscattering_analyzer.utils import ( + compute_light, + process_movement, +) + + +class Experiment: + """ + An experiment in Virgo consisting of exciting (oscillating) a bench + """ + + def __init__( + self, + bench: str, + date: str, + values_file: str, + vibration_frequency: float = 0.2, + ) -> None: + logger.info( + _( + "initializing the experiment on the bench {bench} in {date}" + ).format(bench=bench, date=date) + ) + self.bench = bench + self.date = date + self.vibration_frequency = vibration_frequency + try: + with open(values_file, "rb") as file: + values = load(file) + logger.info( + _("successfully loaded values from {file}").format( + file=values_file + ) + ) + except OSError: + logger.critical( + _("{file} cannot be read").format(file=values_file) + ) + raise Exception( + _("{file] cannot be read").format(file=values_file) + ) + + try: + self.wavelength: float = values["interferometer"][ + "wavelength" + ] + except KeyError: + logger.warning( + _( + "cannot find the wavelength in {file}, default" + + "value will be used" + ).format(file=values_file) + ) + self.wavelength = 1.064e-6 + logger.debug( + _("value of wavelength: {wavelength}").format( + wavelength=self.wavelength + ) + ) + try: + self.calibrations: dict[str, float] = { + "bench": values["benches"][bench]["calib"]["bench"], + "mirror": values["benches"][bench]["calib"]["mirror"], + } + except KeyError as e: + logger.error( + _( + "cannot find the calibration in {file}, default " + + "value will be used, the result will likely be " + + "wrong: {error}" + ).format(file=values_file, error=e) + ) + self.calibrations: dict[str, float] = { + "bench": 1.0, + "mirror": 1.0, + } + logger.debug( + _( + "value of calibrations (bench, mirror): ({bench}, " + + "{mirror})" + ).format( + bench=self.calibrations["bench"], + mirror=self.calibrations["mirror"], + ) + ) + try: + self.factors: dict[str, float | None] = { + "pre": 1e6 + * values["benches"][bench]["scatter_factor"][0], + "true": values["benches"][bench]["scatter_factor"][0], + } + except KeyError: + logger.warning( + _( + "cannot find the scatter factor in {file}, no " + + "value will be set and the scatter factor will " + + "be searched with a fit" + ).format(file=values_file) + ) + self.factors: dict[str, float | None] = { + "pre": None, + "true": None, + } + logger.debug( + _( + "values of scattering factor (outscale, real): ({pre}, " + + "{true})" + ).format( + pre=self.factors["pre"], + true=self.factors["true"], + ) + ) + try: + self.powers: dict[str, float] = { + "laser": values["dates"][date]["laser"], + "dark_fringe": values["dates"][date]["dark_fringe"], + } + except KeyError as e: + logger.error( + _( + "cannot find the power values in {file}, default " + + "one will be used, the result will likely be " + + "wrong: {error}" + ).format(file=values_file, error=e) + ) + self.powers: dict[str, float] = { + "laser": 23, + "dark_fringe": 8e-3, + } + logger.debug( + _( + "values of powers (laser, dark_fringe): ({laser}, " + + "{dark_fringe})" + ).format( + laser=self.powers["laser"], + dark_fringe=self.powers["dark_fringe"], + ) + ) + try: + self.modelisation_file: str = values["dates"][date][ + "modelisation" + ] + except KeyError: + logger.error( + _( + "cannot find the modelisation file linked to this " + + "date in {file}, default one will be used, the " + + "result will likely be wrong" + ).format(file=values_file) + ) + self.modelisation_file = "scatterCouplingO4.mat" + logger.debug( + _("modelisation file: {file}").format( + file=self.modelisation_file, + ) + ) + try: + data_folder_name: str = values["general"]["data_folder"] + except KeyError: + logger.error( + _( + "cannot find the data folder containing the " + + "signals values in {file}, default one will be " + + "set, if its doesn't exist the program will crash" + ).format(file=values_file) + ) + data_folder_name = "/home/demagny/data" + logger.debug( + _("data folder containing signal values: {folder}").format( + folder=data_folder_name, + ) + ) + + data_folder = Path(data_folder_name) + + try: + self.signals: dict[str, Signal] = load_signals( + data_folder, bench, date + ) + logger.debug( + _("excited and reference signal successfully loaded") + ) + except OSError: + logger.critical(_("cannot load signals")) + raise Exception(_("cannot load signals")) + + try: + self.movements: dict[str, Signal] = load_movements( + data_folder, bench, date + ) + logger.debug( + _( + "movements of the bench and the mirror " + + "successfully loaded" + ) + ) + except OSError: + logger.critical(_("cannot load movements")) + raise Exception(_("cannot load movements")) + + self.movements["true"] = process_movement( + self.movements["bench"], + self.movements["mirror"], + self.calibrations, + vibration_frequency, + ) + logger.debug(_("movement processed")) + + try: + self.coupling = load_coupling( + data_folder, + bench, + date, + self.modelisation_file, + self.powers, + ) + logger.debug(_("coupling signals successfully loaded")) + except OSError: + logger.critical(_("cannot load coupling data")) + raise Exception(_("cannot load coupling data")) + + if self.factors["pre"] is None: + logger.debug( + _( + "scattering factor is not defined, the script will " + + "try to find the right one by fitting" + ) + ) + self.factors = self.fit_factors() + + self.projection = self.compute_projection() + logger.debug(_("end of experiment initialisation")) + + def get_factors( + self, + start: float | None = None, + end: float | None = None, + ) -> tuple[Signal, Signal, Signal, Signal, Signal, Signal, float]: + """ + Get factor to compute the projection for a given scatter factor + """ + phase = 4 * pi / self.wavelength + factor_n = (self.movements["true"] * phase).sin().psd().sqrt() + coupling_n = self.coupling[0] + factor_d = (self.movements["true"] * phase).cos().psd().sqrt() + coupling_d = self.coupling[1] + + if start is None: + start = coupling_n.x[0] + if end is None: + end = coupling_n.x[-1] + + coupling_n = coupling_n.cut(start, end) + + ( + factor_n, + coupling_n, + factor_d, + coupling_d, + excited, + reference, + ) = interpolate( + factor_n, + coupling_n, + factor_d, + coupling_d, + self.signals["excited"].psd().sqrt(), + self.signals["reference"].psd().sqrt(), + ) + + return ( + factor_n, + coupling_n, + factor_d, + coupling_d, + excited, + reference, + phase, + ) + + def fit_factors( + self, + start: float = 15, + end: float = 100, + x: NDArray[float64] | None = None, + ) -> dict[str, float]: + """ + Find the best factor (first order only) to get the projection + that best match the excited signal + """ + ( + factor_n, + coupling_n, + factor_d, + coupling_d, + excited, + reference, + phase, + ) = self.get_factors(start=start, end=end) + + if x is None: + x = logspace(-10, 0, 1000) + + y = array( + [ + sum( + abs( + Signal( + factor_n.x, + compute_light( + pre_scatter_factor=x[i], + factor_n=factor_n, + coupling_n=coupling_n, + factor_d=factor_d, + coupling_d=coupling_d, + phase=phase, + ), + ) + + reference + - excited + ).y + ) + for i in range(len(x)) + ] + ) + + pre_scatter_factor: float = x[argmin(y)] + + logger.info( + _("found a scattering factor of {factor}").format( + factor=pre_scatter_factor * 1e-6 + ) + ) + + return { + "pre": pre_scatter_factor, + "true": pre_scatter_factor * 1e-6, + } + + def compute_projection(self) -> Signal: + """ + Compute the projection of the noise from current values + """ + ( + factor_n, + coupling_n, + factor_d, + coupling_d, + _, + _, + phase, + ) = self.get_factors() + + return Signal( + factor_n.x, + compute_light( + pre_scatter_factor=self.factors["pre"], + factor_n=factor_n, + coupling_n=coupling_n, + factor_d=factor_d, + coupling_d=coupling_d, + phase=phase, + ), + ) diff --git a/src/backscattering_analyzer/settings.py b/src/backscattering_analyzer/settings.py deleted file mode 100644 index 2cf8a7a..0000000 --- a/src/backscattering_analyzer/settings.py +++ /dev/null @@ -1,198 +0,0 @@ -from pathlib import Path -from tomllib import load -from rich.console import Console -from rich.theme import Theme - -from rich.traceback import install as traceback_install - - -class Settings: - def __init__(self, options: list[str]) -> None: - self.theme = Theme( - { - "main": "white bold", - "option": "grey50 italic", - "argument": "red", - "description": "green italic", - "warning": "bold red", - } - ) - self.console = Console(theme=self.theme) - _ = traceback_install(console=self.console, show_locals=True) - self.version = "0.1.0" - self.help = ( - "[main]display[/main] [option]\\[options][/option]" - "\n [argument]-b --bench[/argument] [description]bench of the experiment[/description]" - "\n [argument]-d --date[/argument] [description]date of the experiment[/description]" - "\n [argument]-h --help[/argument] [description]print this help and exit[/description]" - "\n [argument]-v --verbose[/argument] [description]be verbose[/description]" - "\n [argument]-V --version[/argument] [description]print version number and exit[/description]" - ) - self.verbose = False - self.bench = "SDB1" - self.date = "2023_03_24" - - with open("values.toml", "rb") as file: - values = load(file) - - self.folder = Path(values["general"]["data_folder"]) - self.wavelength: float = values["interferometer"]["wavelength"] - - index = 0 - while index < len(options): - option = options[index] - try: - if option[0] != "-": - if option == "h": - option = "--help" - if option == "V": - option = "--version" - if "v" in option: - options[index] = option.replace("v", "") - option = "--verbose" - index -= 1 - if option == "b": - index += 1 - try: - options[index] = "--bench={}".format( - options[index] - ) - except IndexError: - raise Exception( - "bench name needed after b option" - ) - if option == "d": - index += 1 - try: - options[index] = "--date={}".format( - options[index] - ) - except IndexError: - raise Exception( - "date needed after d option" - ) - if option[0] == "-": - if option[1] != "-": - if option == "-h": - option = "--help" - if option == "-V": - option = "--version" - if "v" in option: - options[index] = option.replace("v", "") - option = "--verbose" - index -= 1 - if option == "-b": - index += 1 - try: - options[index] = "--bench={}".format( - options[index] - ) - except IndexError: - raise Exception( - "bench name needed after -b option" - ) - if option == "-d": - index += 1 - try: - options[index] = "--date={}".format( - options[index] - ) - except IndexError: - raise Exception( - "date needed after -d option" - ) - if option[:2] == "--": - if option == "--help": - self.console.print(self.help) - exit(0) - if option == "--version": - self.console.print(self.version) - exit(0) - if option == "--verbose": - self.verbose = True - try: - if option[:8] == "--bench=": - self.bench = str(option[8:]).upper() - except IndexError: - continue - try: - if option[:7] == "--date=": - self.date = option[7:] - except IndexError: - continue - else: - raise Exception("unknown option {}".format(option)) - except IndexError: - pass - index += 1 - - if self.bench not in values["benches"].keys(): - raise ValueError("unknow bench {}".format(self.bench)) - - self.calib_mirror: float = values["benches"][self.bench][ - "calib" - ]["mirror"] - self.calib_bench: float = values["benches"][self.bench][ - "calib" - ]["bench"] - self.scattering_factor: list[float] = values["benches"][ - self.bench - ]["scatter_factor"] - - if self.date not in values["dates"].keys(): - raise ValueError("unknow date {}".format(self.date)) - - self.power_in: float = values["dates"][self.date]["power_in"] - self.power_out: float = values["dates"][self.date]["power_out"] - self.modelisation: str = values["dates"][self.date][ - "modelisation" - ] - - self.vibration_frequency = 0.2 - - def bench_file(self) -> Path: - return self.folder / ( - "{bench}_{date}_ben.csv".format( - bench=self.bench, - date=self.date, - ) - ) - - def mirror_file(self) -> Path: - return self.folder / ( - "{bench}_{date}_mir.csv".format( - bench=self.bench, - date=self.date, - ) - ) - - def data_file(self) -> Path: - return self.folder / ( - "{bench}_{date}_dat.csv".format( - bench=self.bench, - date=self.date, - ) - ) - - def reference_file(self) -> Path: - return self.folder / ( - "{bench}_{date}_ref.csv".format( - bench=self.bench, - date=self.date, - ) - ) - - def modelisation_file(self) -> Path: - return self.folder / (self.modelisation) - - def coupling_name(self) -> str: - bench = self.bench - if self.bench == 'SDB2': - bench = 'SDB1' - elif self.bench == 'SIB2': - bench = 'SIB1' - return "{bench}coupling".format(bench=bench) - - def log(self, message: str) -> None: - if self.verbose: - self.console.log(message, _stack_offset=2) diff --git a/src/backscattering_analyzer/utils.py b/src/backscattering_analyzer/utils.py new file mode 100644 index 0000000..43fcd5d --- /dev/null +++ b/src/backscattering_analyzer/utils.py @@ -0,0 +1,45 @@ +from science_signal.signal import Signal +from backscattering_analyzer import logger + +from numpy.typing import NDArray +from numpy import float64, sqrt + + +def process_movement( + bench: Signal, + mirror: Signal, + calibrations: dict[str, float], + vibration_frequency: float, +) -> Signal: + """ + Clean and compute relative movement between the bench and the mirror + """ + logger.info(_("computing relative movement")) + + return ( + ( + bench * calibrations["bench"] + - mirror * calibrations["mirror"] + ) + .detrend("constant") + .filter(end=5 * vibration_frequency) + .detrend("constant") + ) + + +def compute_light( + pre_scatter_factor: float, + factor_n: Signal, + coupling_n: Signal, + factor_d: Signal, + coupling_d: Signal, + phase: float, +) -> NDArray[float64]: + """ + Optimized computing of light with pre-computed factor + """ + return ( + sqrt(pre_scatter_factor) + / phase + * (coupling_n * factor_n + coupling_d * factor_d).y + ) diff --git a/typings/__builtins__.pyi b/typings/__builtins__.pyi new file mode 100644 index 0000000..f247b54 --- /dev/null +++ b/typings/__builtins__.pyi @@ -0,0 +1 @@ +def _(str: str) -> str: ... diff --git a/values.toml b/values.toml index b4921a6..58da462 100644 --- a/values.toml +++ b/values.toml @@ -23,7 +23,7 @@ scatter_factor = [9.5e-9, 0] [benches.SNEB.calib] bench = 1.01 - mirror = 1.101 + mirror = 1.01 [benches.SPRB] scatter_factor = [0, 0] [benches.SPRB.calib] @@ -36,14 +36,14 @@ mirror = 0.98 [dates] [dates.2023_03_24] - power_in = 23 - power_out = 8e-3 + laser = 23 + dark_fringe = 8e-3 modelisation = "scatterCouplingO4.mat" [dates.2023_12_21] - power_in = 12.1 - power_out = 4.445e-3 + laser = 12.1 + dark_fringe = 4.445e-3 modelisation = "scatterCouplingMisalignSR.mat" [dates.2024_01_21] - power_in = 12.1 - power_out = 4.445e-3 + laser = 12.1 + dark_fringe = 4.445e-3 modelisation = "scatterCouplingMisalignSR.mat"