Remove old code
This commit is contained in:
parent
6e5a0d66ff
commit
2bb202a622
GPG key ID:
434A18DA831BC610
1 changed files with 0 additions and 461 deletions
|
|
@ -1,461 +0,0 @@
|
||||||
from pathlib import Path
|
|
||||||
from tomllib import load
|
|
||||||
|
|
||||||
from numpy import argmin, pi, log10, logspace, array
|
|
||||||
from science_signal import interpolate
|
|
||||||
from science_signal.signal import Signal
|
|
||||||
from backscattering_analyzer import logger
|
|
||||||
from backscattering_analyzer.acquisition import AcquisitionType
|
|
||||||
from backscattering_analyzer.data import (
|
|
||||||
Data,
|
|
||||||
load_coupling,
|
|
||||||
)
|
|
||||||
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.debug(
|
|
||||||
_("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)
|
|
||||||
)
|
|
||||||
self.calibrations: dict[str, float] = {
|
|
||||||
"bench": 1.0,
|
|
||||||
"mirror": 1.0,
|
|
||||||
}
|
|
||||||
try:
|
|
||||||
self.calibrations = {
|
|
||||||
"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)
|
|
||||||
)
|
|
||||||
logger.debug(
|
|
||||||
_("value of calibrations (bench, mirror): ({bench}, " + "{mirror})").format(
|
|
||||||
bench=self.calibrations["bench"],
|
|
||||||
mirror=self.calibrations["mirror"],
|
|
||||||
)
|
|
||||||
)
|
|
||||||
self.factors: dict[str, float | None] = {
|
|
||||||
"pre": None,
|
|
||||||
"true": None,
|
|
||||||
}
|
|
||||||
try:
|
|
||||||
self.factors = {
|
|
||||||
"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)
|
|
||||||
)
|
|
||||||
logger.debug(
|
|
||||||
_(
|
|
||||||
"values of scattering factor (outscale, real): ({pre}, " + "{true})"
|
|
||||||
).format(
|
|
||||||
pre=self.factors["pre"],
|
|
||||||
true=self.factors["true"],
|
|
||||||
)
|
|
||||||
)
|
|
||||||
self.powers: dict[str, float] = {
|
|
||||||
"laser": 23,
|
|
||||||
"dark_fringe": 8e-3,
|
|
||||||
}
|
|
||||||
try:
|
|
||||||
self.powers = {
|
|
||||||
"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)
|
|
||||||
)
|
|
||||||
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,
|
|
||||||
)
|
|
||||||
)
|
|
||||||
data_folder_name: str = "/home/demagny/data"
|
|
||||||
try:
|
|
||||||
data_folder_name = 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)
|
|
||||||
)
|
|
||||||
logger.debug(
|
|
||||||
_("data folder containing signal values: {folder}").format(
|
|
||||||
folder=data_folder_name,
|
|
||||||
)
|
|
||||||
)
|
|
||||||
|
|
||||||
data_folder = Path(data_folder_name)
|
|
||||||
|
|
||||||
self.data = Data(data_folder, bench, date)
|
|
||||||
|
|
||||||
self.reference_movement: None | Signal = process_movement(
|
|
||||||
self.data.reference.movements.bench,
|
|
||||||
self.data.reference.movements.mirror,
|
|
||||||
self.calibrations,
|
|
||||||
vibration_frequency,
|
|
||||||
)
|
|
||||||
self.excited_movement: None | Signal = None
|
|
||||||
try:
|
|
||||||
self.excited_movement = process_movement(
|
|
||||||
self.data.excited.movements.bench,
|
|
||||||
self.data.excited.movements.mirror,
|
|
||||||
self.calibrations,
|
|
||||||
vibration_frequency,
|
|
||||||
)
|
|
||||||
except AttributeError:
|
|
||||||
pass
|
|
||||||
logger.debug(_("movement processed"))
|
|
||||||
|
|
||||||
model_powers: dict[str, float] = {
|
|
||||||
"laser": 0,
|
|
||||||
"dark_fringe": 0,
|
|
||||||
}
|
|
||||||
correct_power: bool = False
|
|
||||||
try:
|
|
||||||
correct_power = values["dates"][date]["correct-power"]
|
|
||||||
if correct_power:
|
|
||||||
model_powers = {
|
|
||||||
"laser": values["dates"][date]["coupling"]["laser"],
|
|
||||||
"dark_fringe": values["dates"][date]["coupling"]["dark_fringe"],
|
|
||||||
}
|
|
||||||
except KeyError as error:
|
|
||||||
print(error)
|
|
||||||
logger.warning(
|
|
||||||
_(
|
|
||||||
"cannot find if coupling values should be "
|
|
||||||
+ "corrected due to the diff between power used in "
|
|
||||||
+ "model and the real one. No correction will be "
|
|
||||||
+ "applied"
|
|
||||||
)
|
|
||||||
)
|
|
||||||
|
|
||||||
try:
|
|
||||||
self.coupling = load_coupling(
|
|
||||||
data_folder,
|
|
||||||
bench,
|
|
||||||
date,
|
|
||||||
self.modelisation_file,
|
|
||||||
self.powers,
|
|
||||||
model_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_reference = self.compute_projection(
|
|
||||||
AcquisitionType.REFERENCE,
|
|
||||||
)
|
|
||||||
try:
|
|
||||||
self.projection_excited = self.compute_projection(
|
|
||||||
AcquisitionType.EXCITED,
|
|
||||||
)
|
|
||||||
except Exception:
|
|
||||||
pass
|
|
||||||
logger.debug(_("end of experiment initialisation"))
|
|
||||||
|
|
||||||
def get_factors(
|
|
||||||
self,
|
|
||||||
phase: float,
|
|
||||||
*signals: Signal,
|
|
||||||
start: float | None = None,
|
|
||||||
end: float | None = None,
|
|
||||||
type: AcquisitionType = AcquisitionType.EXCITED,
|
|
||||||
) -> tuple[Signal, ...]:
|
|
||||||
"""
|
|
||||||
Get factor to compute the projection for a given scatter factor
|
|
||||||
"""
|
|
||||||
coupling_n = self.coupling[0]
|
|
||||||
coupling_d = self.coupling[1]
|
|
||||||
|
|
||||||
if type is AcquisitionType.EXCITED:
|
|
||||||
if self.excited_movement is not None:
|
|
||||||
factor_n = (self.excited_movement * phase).sin().psd().sqrt()
|
|
||||||
factor_d = (self.excited_movement * phase).cos().psd().sqrt()
|
|
||||||
else:
|
|
||||||
logger.critical(
|
|
||||||
_(
|
|
||||||
"no excited signal given, cannot get factors "
|
|
||||||
+ "of the excited signal"
|
|
||||||
)
|
|
||||||
)
|
|
||||||
raise Exception(
|
|
||||||
_(
|
|
||||||
"no excited signal given, cannot get factors "
|
|
||||||
+ "of the excited signal"
|
|
||||||
)
|
|
||||||
)
|
|
||||||
elif type is AcquisitionType.REFERENCE:
|
|
||||||
if self.reference_movement is not None:
|
|
||||||
factor_n = (self.reference_movement * phase).sin().psd().sqrt()
|
|
||||||
factor_d = (self.reference_movement * phase).cos().psd().sqrt()
|
|
||||||
else:
|
|
||||||
logger.critical(
|
|
||||||
_(
|
|
||||||
"no reference signal given, cannot get factors "
|
|
||||||
+ "of the reference signal"
|
|
||||||
)
|
|
||||||
)
|
|
||||||
raise Exception(
|
|
||||||
_(
|
|
||||||
"no reference signal given, cannot get factors "
|
|
||||||
+ "of the reference signal"
|
|
||||||
)
|
|
||||||
)
|
|
||||||
else:
|
|
||||||
logger.critical(_("unknown acquisition type"))
|
|
||||||
raise Exception(_("unknown acquisition type"))
|
|
||||||
|
|
||||||
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)
|
|
||||||
|
|
||||||
return interpolate(
|
|
||||||
factor_n,
|
|
||||||
coupling_n,
|
|
||||||
factor_d,
|
|
||||||
coupling_d,
|
|
||||||
*signals,
|
|
||||||
)
|
|
||||||
|
|
||||||
def fit_factors(
|
|
||||||
self,
|
|
||||||
scatter_min: float = 1e-20,
|
|
||||||
scatter_max: float = 1,
|
|
||||||
start_frequency: float = 20,
|
|
||||||
end_frequency: float = 35,
|
|
||||||
precision: int = 3,
|
|
||||||
) -> dict[str, None | float]:
|
|
||||||
"""
|
|
||||||
Find the best factor (first order only) to get the projection
|
|
||||||
that best match the excited signal
|
|
||||||
"""
|
|
||||||
phase = 4 * pi / self.wavelength
|
|
||||||
(
|
|
||||||
factor_n,
|
|
||||||
coupling_n,
|
|
||||||
factor_d,
|
|
||||||
coupling_d,
|
|
||||||
excited,
|
|
||||||
reference,
|
|
||||||
) = self.get_factors(
|
|
||||||
phase,
|
|
||||||
self.data.excited.sensibility.psd().sqrt(),
|
|
||||||
self.data.reference.sensibility.psd().sqrt(),
|
|
||||||
start=start_frequency,
|
|
||||||
end=end_frequency,
|
|
||||||
type=AcquisitionType.EXCITED,
|
|
||||||
)
|
|
||||||
|
|
||||||
scatter_max = 1e10
|
|
||||||
scatter_min = 1e-10
|
|
||||||
|
|
||||||
for index in range(precision):
|
|
||||||
logger.debug(_("search for a local minimum"))
|
|
||||||
|
|
||||||
x = logspace(log10(scatter_min), log10(scatter_max), 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))
|
|
||||||
]
|
|
||||||
)
|
|
||||||
from matplotlib.pyplot import figure, show
|
|
||||||
|
|
||||||
diff_1 = abs(
|
|
||||||
Signal(
|
|
||||||
factor_n.x,
|
|
||||||
compute_light(
|
|
||||||
pre_scatter_factor=x[100],
|
|
||||||
factor_n=factor_n,
|
|
||||||
coupling_n=coupling_n,
|
|
||||||
factor_d=factor_d,
|
|
||||||
coupling_d=coupling_d,
|
|
||||||
phase=phase,
|
|
||||||
),
|
|
||||||
)
|
|
||||||
+ reference
|
|
||||||
- excited
|
|
||||||
)
|
|
||||||
diff_2 = abs(
|
|
||||||
Signal(
|
|
||||||
factor_n.x,
|
|
||||||
compute_light(
|
|
||||||
pre_scatter_factor=x[900],
|
|
||||||
factor_n=factor_n,
|
|
||||||
coupling_n=coupling_n,
|
|
||||||
factor_d=factor_d,
|
|
||||||
coupling_d=coupling_d,
|
|
||||||
phase=phase,
|
|
||||||
),
|
|
||||||
)
|
|
||||||
+ reference
|
|
||||||
- excited
|
|
||||||
)
|
|
||||||
|
|
||||||
fig = figure()
|
|
||||||
fig.gca().loglog(diff_1.x, diff_1.y, label=str(x[100]))
|
|
||||||
fig.gca().loglog(diff_2.x, diff_2.y, label=str(x[900]))
|
|
||||||
fig.gca().legend()
|
|
||||||
show()
|
|
||||||
|
|
||||||
if argmin(y) == 0:
|
|
||||||
logger.warning(_("smaller than current range allows"))
|
|
||||||
scatter_max = x[1]
|
|
||||||
elif argmin(y) == len(x) - 1:
|
|
||||||
logger.warning(_("bigger than current range allows"))
|
|
||||||
scatter_min = x[-2]
|
|
||||||
else:
|
|
||||||
scatter_min = x[argmin(y) - 1]
|
|
||||||
scatter_max = x[argmin(y) + 1]
|
|
||||||
|
|
||||||
logger.debug(_("local minimum found"))
|
|
||||||
|
|
||||||
pre_scatter_factor = scatter_min / 2 + scatter_max / 2
|
|
||||||
|
|
||||||
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, type: AcquisitionType = AcquisitionType.EXCITED
|
|
||||||
) -> Signal:
|
|
||||||
"""
|
|
||||||
Compute the projection of the noise from current values
|
|
||||||
"""
|
|
||||||
phase = 4 * pi / self.wavelength
|
|
||||||
(
|
|
||||||
factor_n,
|
|
||||||
coupling_n,
|
|
||||||
factor_d,
|
|
||||||
coupling_d,
|
|
||||||
) = self.get_factors(phase, type=type)
|
|
||||||
|
|
||||||
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,
|
|
||||||
),
|
|
||||||
)
|
|
||||||
Loading…
Reference in a new issue