Refactorisation
This commit is contained in:
parent
b8d54cdd21
commit
c60e8a0e18
9 changed files with 624 additions and 540 deletions
|
|
@ -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]
|
||||
|
|
|
|||
|
|
@ -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"))
|
||||
|
|
|
|||
|
|
@ -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
|
||||
161
src/backscattering_analyzer/data.py
Normal file
161
src/backscattering_analyzer/data.py
Normal file
|
|
@ -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
|
||||
379
src/backscattering_analyzer/experiment.py
Normal file
379
src/backscattering_analyzer/experiment.py
Normal file
|
|
@ -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,
|
||||
),
|
||||
)
|
||||
|
|
@ -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)
|
||||
45
src/backscattering_analyzer/utils.py
Normal file
45
src/backscattering_analyzer/utils.py
Normal file
|
|
@ -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
|
||||
)
|
||||
1
typings/__builtins__.pyi
Normal file
1
typings/__builtins__.pyi
Normal file
|
|
@ -0,0 +1 @@
|
|||
def _(str: str) -> str: ...
|
||||
14
values.toml
14
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"
|
||||
|
|
|
|||
Loading…
Reference in a new issue