First commit

MaskedReadoutDC.py

@@ -0,0 +1,32 @@
+from typing import Literal
+
+from finesse.model import Node
+from finesse.components.general import borrows_nodes
+
+from finesse.components.readout import ReadoutDC
+from finesse.detectors.general import MaskedDetector
+
+
+@borrows_nodes()
+class MaskedReadoutDC(ReadoutDC, MaskedDetector):
+    def __init__(
+        self,
+        name: str,
+        optical_node: Node,
+        *,
+        output_detectors: bool = False,
+        pdtype: None | str | dict[str, str] = None,
+    ):
+        super().__init__(
+            name,
+            optical_node,
+            output_detectors=output_detectors,
+            pdtype=pdtype,
+        )
+        MaskedDetector.__init__(self, name)
+
+    @property
+    def has_mask(
+        self,
+    ) -> Literal[False]:  # should return true here, but cannot override
+        return False

main.py

@@ -0,0 +1,533 @@
+from logging import getLogger, basicConfig, WARNING
+from argparse import ArgumentParser
+from pathlib import Path
+from urllib.parse import urlparse
+from importlib.metadata import version, PackageNotFoundError
+
+from rich.theme import Theme
+from rich.console import Console
+from rich.progress import Progress
+from rich.logging import RichHandler
+from rich.table import Table
+
+from matplotlib.pyplot import figure, show
+
+from numpy import sqrt, geomspace, mean, diff
+
+from finesse.model import Model
+from finesse.analysis.actions import (
+    TemporaryParameters,
+    Change,
+    Maximize,
+    Minimize,
+    Series,
+    Noxaxis,
+)
+from finesse.exceptions import ModelMissingAttributeError
+
+from astropy.units import W, rad, m, mW  # pyright: ignore[reportAttributeAccessIssue]
+
+from gwpy.frequencyseries import FrequencySeries
+
+from h5py import File
+
+from git import Repo
+
+from backscattering_simulation_data import Root, Scatterer
+from backscattering_simulation_data.metadata import (
+    MetadataRoot,
+    MetadataScatterer,
+    Interferometer,
+    Git,
+)
+
+from utils import fix_dark_fringe, TF
+from MaskedReadoutDC import MaskedReadoutDC
+
+theme = Theme(
+    {
+        "strong": "cyan underline",
+        "result": "red bold",
+    }
+)
+console = Console(theme=theme)
+
+basicConfig(
+    level=WARNING,
+    format="%(message)s",
+    datefmt="[%X]",
+    handlers=[RichHandler(console=console, rich_tracebacks=True)],
+)
+logger = getLogger(__name__)
+
+# Parameters
+
+parser = ArgumentParser(
+    prog="main",
+    description="Generate a model file containing computed Transfer Function",
+    usage="python main.py --filename={filename} --laser-power={laser power} --dark-fringe={dark fringe power} --sr-rotation={sr rotation} --max-tem={maxtem} [--precision={precision} | --quiet | --verbose]",
+)
+parser.add_argument(
+    "-f",
+    "--filename",
+    help="path to the new file containing the transfer function",
+    type=Path,
+    required=True,
+)
+parser.add_argument(
+    "-p",
+    "--laser-power",
+    help="power injected in the interferometer (in W)",
+    type=float,
+    required=True,
+)
+parser.add_argument(
+    "-d",
+    "--dark-fringe",
+    help="power on the dark fringe (dark fringe offset) (in W)",
+    type=float,
+    required=True,
+)
+parser.add_argument(
+    "-r",
+    "--sr-rotation",
+    help="yaw SR misalignement angle (in rad)",
+    type=float,
+    required=True,
+)
+parser.add_argument(
+    "-m",
+    "--max-tem",
+    help="number of modes in the model",
+    type=int,
+    required=True,
+)
+parser.add_argument(
+    "-n",
+    "--precision",
+    help="number of point in the transfer function model",
+    type=int,
+    default=1000,
+    required=False,
+)
+parser.add_argument(
+    "-q",
+    "--quiet",
+    help="suppress outputs in the terminal",
+    action="store_true",
+)
+parser.add_argument(
+    "-V",
+    "--verbose",
+    help="display graph and table on the terminal",
+    action="store_true",
+)
+parser.add_argument(
+    "-a",
+    "--name",
+    help="name of the generated result. Should be unique",
+    type=str,
+    required=True,
+)
+
+args = parser.parse_args()
+
+filename = args.filename  # name of the created file
+laser_power = args.laser_power * W  # power injected
+dark_fringe_power = args.dark_fringe * W  # power on dark fringe
+precision = args.precision  # number of point in the simulation
+sr_rotation = args.sr_rotation * rad  # yaw SR misalignement angle
+maxtem = args.max_tem  # number of modes in the model
+quiet = args.quiet  # no progress bar
+verbose = args.verbose  # display figure and tables
+name = args.name  # unique name of the generated hdf5 data
+
+if filename.suffix != ".hdf5":
+    filename = filename.with_suffix(".hdf5")
+
+with Progress(console=console) as progress:
+    model_building = progress.add_task(
+        description="building the model",
+        start=False,
+        total=None,
+        visible=False,
+    )
+    locking = progress.add_task(
+        description="locking the model",
+        start=False,
+        total=None,
+        visible=False,
+    )
+    set_dark_fringe = progress.add_task(
+        description="set dark fringe offset",
+        start=False,
+        total=None,
+        visible=False,
+    )
+    computing_TF = progress.add_task(
+        description="computing transfer function",
+        start=False,
+        total=8,
+        visible=False,
+    )
+    write_file = progress.add_task(
+        description="writing file",
+        start=False,
+        total=None,
+        visible=False,
+    )
+
+    C_B1_DETECTOR = "TF"  # depends on the code itself (don't change)
+
+    # Creating model
+
+    progress.start_task(model_building)
+    progress.update(model_building, visible=not quiet)
+
+    model = Model()
+    model.phase_config(zero_k00=False, zero_tem00_gouy=True)
+    model.parse(Path("model.kat").read_text())
+    model.lambda0 = model.get("wavelength")
+    model.get("fsig").f = 1.0
+
+    try:
+        model.get("SR").xbeta = sr_rotation.to(rad).value
+        model.get("laser").P = laser_power.to(W).value
+    except ModelMissingAttributeError:
+        logger.warning("SR, or laser is not define in the model")
+
+    try:
+        model.get("B1")
+        logger.info("B1 already exists")
+    except ModelMissingAttributeError:
+        try:
+            model.add(
+                MaskedReadoutDC(
+                    "B1",
+                    optical_node=model.get("SDB1").p2.o,
+                    output_detectors=True,
+                    pdtype=None,
+                )
+            )
+            logger.info("B1 created")
+        except ModelMissingAttributeError:
+            logger.warning("SDB1 does not exist in the model")
+
+    progress.remove_task(model_building)
+    progress.start_task(locking)
+    progress.update(locking, visible=not quiet)
+    # Locking
+
+    result = model.run(
+        TemporaryParameters(
+            Series(
+                Change(
+                    {
+                        "SR.misaligned": True,
+                        "PR.misaligned": True,
+                    }
+                ),
+                Maximize(
+                    model.get("NE_p1"),
+                    model.get("NORTH_ARM.DC"),
+                    bounds=[-180, 180],
+                    tol=1e-14,
+                ),
+                Maximize(
+                    model.get("WE_p1"),
+                    model.get("WEST_ARM.DC"),
+                    bounds=[-180, 180],
+                    tol=1e-14,
+                ),
+                Minimize(
+                    model.get("SR_p2"),
+                    model.get("MICH.DC"),
+                    bounds=[-180, 180],
+                    tol=1e-14,
+                ),
+                Change(
+                    {
+                        "PR.misaligned": False,
+                    },
+                ),
+                Maximize(
+                    model.get("PR_p2"),
+                    model.get("PRCL.DC"),
+                    bounds=[-180, 180],
+                    tol=1e-14,
+                ),
+                Change(
+                    {
+                        "SR.misaligned": False,
+                    },
+                ),
+                Maximize(
+                    model.get("B1_DC"),
+                    model.get("SRCL.DC"),
+                    bounds=[-180, 180],
+                    tol=1e-14,
+                ),
+                Change(
+                    {
+                        "SRCL.DC": -90,
+                    },
+                    relative=True,
+                ),
+            ),
+            exclude=(
+                "NE.phi",
+                "NI.Phi",
+                "WE.phi",
+                "WI.phi",
+                "SR.phi",
+                "PR.phi",
+                "NORTH_ARM.DC",
+                "WEST_ARM.DC",
+                "DARM.DC",
+                "MICH.DC",
+                "PRCL.DC",
+                "SRCL.DC",
+                "SR.misaligned",
+                "PR.misaligned",
+            ),
+        ),
+    )
+
+    # Adding HOMs
+    model.modes(maxtem=maxtem)
+    model.get("B1").select_mask(exclude=(0, 0))
+
+    progress.remove_task(locking)
+
+    if verbose:
+        solution = model.run(Noxaxis())
+        if solution is None:
+            raise Exception("error when modeling the locked model")
+        table = Table(title="Puissances dans l'interferomètre")
+        table.add_column("position", justify="left", style="white")
+        table.add_column("puissance (W)", justify="left", style="cyan")
+
+        table.add_row("Injection", str(model.get("laser").P.eval()))
+        table.add_row("PR", str(solution["PR_p1"]))  # pyright: ignore[reportArgumentType,reportCallIssue]
+        table.add_row(
+            "cavité de recyclage de puissance",
+            str(solution["PR_p2"]),  # pyright: ignore[reportArgumentType,reportCallIssue]
+        )
+        table.add_row("cavité ouest", str(solution["WE_p1"]))  # pyright: ignore[reportArgumentType,reportCallIssue]
+        table.add_row("cavité nord", str(solution["NE_p1"]))  # pyright: ignore[reportArgumentType,reportCallIssue]
+        table.add_row("frange noire", str(solution["SR_p2"]))  # pyright: ignore[reportArgumentType,reportCallIssue]
+        table.add_row("SNEB", str(solution["SNEB_DC"]))  # pyright: ignore[reportArgumentType,reportCallIssue]
+        table.add_row("SWEB", str(solution["SWEB_DC"]))  # pyright: ignore[reportArgumentType,reportCallIssue]
+        table.add_row("SDB1", str(solution["SDB1_DC"]))  # pyright: ignore[reportArgumentType,reportCallIssue]
+
+        console.print(table)
+
+    progress.start_task(set_dark_fringe)
+    progress.update(set_dark_fringe, visible=not quiet)
+
+    number, power = fix_dark_fringe(
+        model, dark_fringe_power.to(W).value
+    )
+    power = power * W
+
+    if verbose:
+        console.print(
+            "[result]{dof}[/result] with [result]{number} steps[/result] give [result]{power:.6f}[/result] on the [strong]dark fringe[/strong]".format(
+                number=number,
+                dof=model.get("DARM").DC,
+                power=power.to(mW),
+            ),
+        )
+    dark_fringe_power = (
+        power  # to keep this result for the h5df metadata
+    )
+
+    progress.remove_task(set_dark_fringe)
+    progress.start_task(computing_TF)
+    progress.update(computing_TF, visible=not quiet)
+
+    # Computing TF
+
+    model.get("SNEB").phi = model.get("NE").phi - 45
+    model.get("SWEB").phi = model.get("WE").phi - 45
+    model.get("SDB1").phi = model.get("SR").phi + 45
+
+    quad_tf: dict[str, TF] = dict()  # bench is in phase quadrature (Kn)
+    in_tf: dict[str, TF] = dict()  # bench is in phase (Kp)
+
+    for bench_name in ["SNEB", "SWEB", "SDB1"]:
+        quad_tf[bench_name] = TF(
+            model,
+            ["{}_z".format(bench_name)],
+            [C_B1_DETECTOR],
+            geomspace(5, 10000, precision),  # pyright: ignore[reportArgumentType]
+        )
+        progress.update(computing_TF, advance=1)
+
+    quad_tf["DARM"] = TF(
+        model,
+        ["DARM"],
+        [C_B1_DETECTOR],
+        geomspace(5, 10000, precision),  # pyright: ignore[reportArgumentType]
+    )
+    progress.update(computing_TF, advance=1)
+
+    model.get("SNEB").phi = model.get("NE").phi
+    model.get("SWEB").phi = model.get("WE").phi
+    model.get("SDB1").phi = model.get("SR").phi
+
+    for bench_name in ["SNEB", "SWEB", "SDB1"]:
+        in_tf[bench_name] = TF(
+            model,
+            ["{}_z".format(bench_name)],
+            [C_B1_DETECTOR],
+            geomspace(5, 10000, precision),  # pyright: ignore[reportArgumentType]
+        )
+        progress.update(computing_TF, advance=1)
+
+    in_tf["DARM"] = TF(
+        model,
+        ["DARM"],
+        [C_B1_DETECTOR],
+        geomspace(5, 10000, precision),  # pyright: ignore[reportArgumentType]
+    )
+    progress.update(computing_TF, advance=1)
+
+    bench_names = ["SNEB", "SWEB", "SDB1"]
+
+    if verbose:
+        for bench_name in bench_names:
+            Figure = figure()
+            ax = Figure.gca()
+            ax.set_title(
+                "finesse {bench_name}".format(bench_name=bench_name)
+            ).figure.gca().loglog(
+                quad_tf[bench_name].f,
+                sqrt(
+                    (
+                        abs(quad_tf[bench_name].get(C_B1_DETECTOR))
+                        / abs(quad_tf["DARM"].get(C_B1_DETECTOR))
+                        / model.get("space_NI_NE").L.eval()
+                    )
+                    ** 2
+                    + (
+                        abs(in_tf[bench_name].get(C_B1_DETECTOR))
+                        / abs(in_tf["DARM"].get(C_B1_DETECTOR))
+                        / model.get("space_NI_NE").L.eval()
+                    )
+                    ** 2
+                ),
+                label="sum of TF ($\\sqrt{{K_n^2 + K_p^2}}$)",
+            )[0].figure.gca().loglog(
+                quad_tf[bench_name].f,
+                abs(quad_tf[bench_name].get(C_B1_DETECTOR))
+                / abs(quad_tf["DARM"].get(C_B1_DETECTOR))
+                / model.get("space_NI_NE").L.eval(),
+                label="Kn",
+            )[0].figure.gca().loglog(
+                in_tf[bench_name].f,
+                abs(in_tf[bench_name].get(C_B1_DETECTOR))
+                / abs(in_tf["DARM"].get(C_B1_DETECTOR))
+                / model.get("space_NI_NE").L.eval(),
+                label="Kp",
+            )[0].figure.gca().set_ylabel(
+                "$\\frac{{m}}{{m}}$"
+            ).figure.gca().legend().figure.gca().set_xlabel(
+                "Frequencies ($Hz$)"
+            ).figure.gca().grid(True, "both", "both")
+        show()
+
+    progress.remove_task(computing_TF)
+    progress.start_task(write_file)
+    progress.update(write_file, visible=not quiet)
+
+    repo = Repo(Path(__file__).parent)
+
+    current_version = "unknown"
+    try:
+        current_version = version(__name__)
+    except PackageNotFoundError:
+        pass
+
+    data = Root(
+        MetadataRoot(
+            name=name,
+            maxtem=maxtem,
+            git=Git(
+                remote=urlparse(repo.remote().url),
+                commit=repo.active_branch.commit.hexsha,
+                version=current_version,
+                branch=repo.active_branch.name,
+            ),
+            interferometer=Interferometer(
+                dark_fringe_power=dark_fringe_power,
+                injection_power=laser_power,
+                sr_yaw=sr_rotation,
+                north_arm=model.get("space_NI_NE").L.eval() * m,
+                west_arm=model.get("space_WI_WE").L.eval() * m,
+            ),
+        ),
+        [
+            Scatterer(
+                MetadataScatterer(
+                    "SDB1",
+                ),
+                kn=FrequencySeries(
+                    data=quad_tf["SDB1"].get(C_B1_DETECTOR),
+                    unit=m / W,
+                    f0=quad_tf["SDB1"].f[0],
+                    df=mean(diff(quad_tf["SDB1"].f)),
+                ),
+                kp=FrequencySeries(
+                    data=in_tf["SDB1"].get(C_B1_DETECTOR),
+                    unit=m / W,
+                    f0=in_tf["SDB1"].f[0],
+                    df=mean(diff(in_tf["SDB1"].f)),
+                ),
+            ),
+            Scatterer(
+                MetadataScatterer(
+                    "SNEB",
+                ),
+                kn=FrequencySeries(
+                    data=quad_tf["SNEB"].get(C_B1_DETECTOR),
+                    unit=m / W,
+                    f0=quad_tf["SNEB"].f[0],
+                    df=mean(diff(quad_tf["SNEB"].f)),
+                ),
+                kp=FrequencySeries(
+                    data=in_tf["SNEb"].get(C_B1_DETECTOR),
+                    unit=m / W,
+                    f0=in_tf["SNEB"].f[0],
+                    df=mean(diff(in_tf["SNEB"].f)),
+                ),
+            ),
+            Scatterer(
+                MetadataScatterer(
+                    "SWEB",
+                ),
+                kn=FrequencySeries(
+                    data=quad_tf["SWEb"].get(C_B1_DETECTOR),
+                    unit=m / W,
+                    f0=quad_tf["SWEB"].f[0],
+                    df=mean(diff(quad_tf["SWEB"].f)),
+                ),
+                kp=FrequencySeries(
+                    data=in_tf["SWEB"].get(C_B1_DETECTOR),
+                    unit=m / W,
+                    f0=in_tf["SWEB"].f[0],
+                    df=mean(diff(in_tf["SWEB"].f)),
+                ),
+            ),
+        ],
+        FrequencySeries(
+            data=quad_tf["DARM"].get(C_B1_DETECTOR),
+            unit=W / m,  # power on B1 to arm displacement
+            f0=quad_tf["DARM"].f[0],
+            df=mean(diff(quad_tf["DARM"].f)),
+        ),
+    )
+    with File(name=filename, mode="w") as file:
+        data.to_hdf5(file)
+
+    progress.remove_task(write_file)

mise.toml

@@ -0,0 +1,5 @@
+[tools]
+python = "3.13.2"
+
+[env]
+_.python.venv = { path = ".venv", create = true, uv_create_args = ["--seed"] }

model.kat

@@ -0,0 +1,182 @@
+#-----------------------------------------------------------------------
+# An Advanced Virgo Plus input file for Finesse 3
+
+# Source: current finesse-virgo's model
+# Modified with: https://git.ligo.org/virgo/isc/parameters/-/blob/377982cc0fd35a561a2f0d9356ca847072d39c57/mechanics.m
+# Modified with: https://git.ligo.org/virgo/isc/parameters/-/blob/c7e2cf038a096fbb5a20cc58ff01650cec5a96fa/advanced_virgo.m
+# Modified with: https://git.ligo.org/virgo/isc/parameters/-/blob/7de842d4f3ebe3744a608e2f19d6d0a15a437cdf/advanced_virgo_current.m
+# Modified with: https://git.ligo.org/virgo/isc/parameters/-/blob/a89dfb3b8556c9796b64531da93343f56e773e6b/detection_current.m
+# Modified with: https://git.ligo.org/virgo/isc/optickle/-/blob/f577e6f073c0dda2ef9db3d2fbab27021f406182/models/DRITFscatter.m
+# Modified with: https://git.ligo.org/virgo/virgoapp/simulinkNoiseBudget/-/blob/6e7e1b95b7817cdb6cc67d5748281ed1d7ffc805/DRITF/scatterRadiationPressureTune.m
+# Modified with: https://git.ligo.org/virgo/virgoapp/simulinkNoiseBudget/-/blob/6e7e1b95b7817cdb6cc67d5748281ed1d7ffc805/DRITF/scatterRadiationPressure.m
+
+# Defines general variables needed for the model creation here
+
+# Constant
+variable	speed_of_light	299792458
+
+# Material
+variable	silica	1.44963
+
+# Injection
+variable	wavelength	1.064e-6	units='m'
+variable	frequency	speed_of_light/wavelength	units='Hz'
+
+# Schnupp asymmetry
+variable	schnupp_asymetry	0.23	units='m'
+
+# Cavity length
+variable	length_PRCL	11.952	units='m'
+variable	length_SRCL	11.952	units='m'
+
+# Laser
+laser	laser	P=25
+
+# Modulator
+
+#modulator	eom1	f=6270777	midx=0.22	order=1
+#space	space_laser_eom1	laser.p1	eom1.p1	L=0.1
+
+#modulator	eom2	f=4/3*eom1.f	midx=0.15	order=1
+#space	space_eom1_eom2	eom1.p2	eom2.p1	L=0.1
+
+#modulator	eom3	f=9*eom1.f	midx=0.204	order=1
+#space	space_eom2_eom3	eom2.p2	eom3.p1	L=0.1
+
+#modulator	eom4	f=21*eom1.f	midx=0.128	order=1
+#space	space_eom3_eom4	eom3.p2	eom4.p1	L=0.1
+
+# Power Recycling
+mirror	PR_AR	T=9.9984e-1	L=1.6e-4	Rc=-3.62	phi=PR.phi
+mirror	PR	T=0.04835	L=37.5e-6	Rc=-1430
+space	inside_PR	PR_AR.p2	PR.p1	L=0.1003	nr=silica
+
+space	space_laser_PR	laser.p1	PR_AR.p1	L=11.67
+
+# Beamsplitter
+
+beamsplitter	BS	T=0.5012	L=37.5e-6	alpha=-45.0
+
+space	space_PR_BS	PR.p2	BS.p1	L=5.925
+
+# North arm
+
+mirror	NI_AR	T=9.9968e-1	L=0.0	Rc=-1420	phi=NI.phi
+mirror	NI	T=0.01377	L=37.5e-6	Rc=-1420
+space	inside_NI	NI_AR.p2	NI.p1	L=0.2	nr=silica
+
+mirror	NE	T=4.4e-6	L=37.5e-6	Rc=1683
+mirror	NE_AR	R=100e-6	L=1.33e-4	phi=NE.phi	Rc=1683
+space	inside_NE	NE.p2	NE_AR.p1	L=0.2	nr=silica
+
+space	space_NI_NE	NI.p2	NE.p1	L=2999.8
+
+space	space_BS_NI	BS.p3	NI_AR.p1	L=length_PRCL-space_PR_BS.L+schnupp_asymetry/2-inside_NI.L*inside_NI.nr
+
+# West arm
+
+mirror	WI_AR	T=9.99942e-1	L=0.0	Rc=-1424	phi=WI.phi
+mirror	WI	T=0.01377	L=37.5e-6	Rc=-1420
+space	inside_WI	WI_AR.p2	WI.p1	L=0.2	nr=silica
+
+mirror	WE	T=4.3e-6	L=37.5e-6	Rc=1683
+mirror	WE_AR	T=9.99845e-1	L=1.55e-4	phi=WE.phi	Rc=1683
+space	inside_WE	WE.p2	WE_AR.p1	L=0.2	nr=silica
+
+space	space_WI_WE	WI.p2	WE.p1	L=2999.8
+
+# schnupp asymetry correction
+space	space_BS_WI	BS.p2	WI_AR.p1	L=length_PRCL-space_PR_BS.L-schnupp_asymetry/2-inside_WI.L*inside_NI.nr
+
+# SR
+
+mirror	SR	T=0.40	L=37.5e-6	Rc=1430
+mirror	SR_AR	T=9.99859e-1	L=1.41e-4	Rc=1430	phi=SR.phi
+space	inside_SR	SR.p2	SR_AR.p1	L=0.1004
+
+space	space_BS_SR	BS.p4	SR.p1	L=length_SRCL-(space_BS_WI.L+space_BS_NI.L)/2
+
+# power detectors
+
+power_detector_dc	NE_p1	NE.p1.o
+power_detector_dc	NI_p1	NI.p1.o
+power_detector_dc	NE_p2	NE.p2.o
+power_detector_dc	NI_p2	NI.p2.o
+
+power_detector_dc	WE_p1	WE.p1.o
+power_detector_dc	WI_p1	WI.p1.o
+power_detector_dc	WE_p2	WE.p2.o
+power_detector_dc	WI_p2	WI.p2.o
+
+power_detector_dc	PR_p1	PR.p1.o
+power_detector_dc	PR_p2	PR.p2.o
+
+power_detector_dc	SR_p1	SR.p1.o
+power_detector_dc	SR_p2	SR.p2.o
+
+# degree of freedom
+
+degree_of_freedom	NORTH_ARM	NE.dofs.z	+1
+degree_of_freedom	WEST_ARM	WE.dofs.z	+1
+degree_of_freedom	PRCL	PR.dofs.z	+1
+degree_of_freedom	SRCL	SR.dofs.z	+1
+
+degree_of_freedom MICH	NI.dofs.z	-0.5	WI.dofs.z	+0.5	NE.dofs.z	-0.5	WE.dofs.z	+0.5
+degree_of_freedom DARM	NE.dofs.z	-0.5	WE.dofs.z	+0.5
+degree_of_freedom CARM	NE.dofs.z	+1	WE.dofs.z	+1
+degree_of_freedom	DARM_Fz	NE.dofs.F_z	-0.5	WE.dofs.F_z	+0.5
+
+# cavities
+
+cavity	north_arm	NI.p2.o	priority=3
+cavity	west_arm	WI.p2.o	priority=3
+
+cavity	PRCL_north	PR.p2.o	via=WE.p1.i	priority=2
+cavity	PRCL_west	PR.p2.o	via=NE.p1.i	priority=2
+
+cavity	SRCL_north	SR.p1.o	via=WE.p1.i	priority=1
+cavity	SRCL_west	SR.p1.o	via=WE.p1.i	priority=1
+
+# benches
+
+mirror	SNEB	R=1e-6	L=0.0	Rc=1683+7
+mirror	SWEB	R=1e-6	L=0.0	Rc=1683+7
+mirror	SDB1	R=1e-6	L=0.0	Rc=1430
+
+space	space_NE_AR_SNEB	NE_AR.p2	SNEB.p1	L=7
+space	space_WE_AR_SWEB	WE_AR.p2	SWEB.p1	L=7
+space	space_SR_SDB1	SR_AR.p2	SDB1.p1	L=0
+
+power_detector_dc	SNEB_DC	SNEB.p1.i
+power_detector_dc	SWEB_DC	SWEB.p1.i
+power_detector_dc	SDB1_DC	SDB1.p1.i
+
+degree_of_freedom	SNEB_z	SNEB.dofs.z	-1
+degree_of_freedom	SWEB_z	SWEB.dofs.z	-1
+degree_of_freedom	SDB1_z	SDB1.dofs.z	-1
+
+# Suspensions
+
+pendulum	NI_suspension	NI.mech	mass=42	fz=0.76	Qz=40
+pendulum	NI_AR_suspension	NI_AR.mech	mass=NI_suspension.mass	fz=NI_suspension.fz	Qz=NI_suspension.Qz
+
+pendulum	NE_suspension	NE.mech	mass=42	fz=0.76	Qz=40
+pendulum	NE_AR_suspension	NE_AR.mech	mass=NE_suspension.mass	fz=NE_suspension.fz	Qz=NE_suspension.Qz
+
+pendulum	WI_suspension	WI.mech	mass=42	fz=0.76	Qz=40
+pendulum	WI_AR_suspension	WI_AR.mech	mass=WI_suspension.mass	fz=WI_suspension.fz	Qz=WI_suspension.Qz
+
+pendulum	WE_suspension	WE.mech	mass=42	fz=0.76	Qz=40
+pendulum	WE_AR_suspension	WE_AR.mech	mass=WE_suspension.mass	fz=WE_suspension.fz	Qz=WE_suspension.Qz
+
+pendulum	PR_suspension	PR.mech	mass=42	fz=0.76	Qz=40
+pendulum	PR_AR_suspension	PR_AR.mech	mass=PR_suspension.mass	fz=PR_suspension.fz	Qz=PR_suspension.Qz
+
+pendulum	SR_suspension	SR.mech	mass=42	fz=0.76	Qz=40
+pendulum	SR_AR_suspension	SR_AR.mech	mass=SR_suspension.mass	fz=SR_suspension.fz	Qz=SR_suspension.Qz
+
+pendulum	SNEB_suspension	SNEB.mech	mass=42	fz=0.76	Qz=40
+pendulum	SWEB_suspension	SWEB.mech	mass=42	fz=0.76	Qz=40
+pendulum	SDB1_suspension	SDB1.mech	mass=42	fz=0.76	Qz=40
+
+free_mass	BS_suspension	BS.mech	mass=34

pyproject.toml

@@ -0,0 +1,45 @@
+[project]
+name = "backscattering_simulation"
+version = "0.0.1"
+authors = [{ name = "linarphy", email = "linarphy@linarphy.net" }]
+description = "Generate simulation data used to measure backscattered light in the Virgo interferometer"
+readme = "README.md"
+requires-python = ">3.11"
+classifiers = [
+  "Programming Language :: Python :: 3",
+  "License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)",
+  "Intended Audience :: Science/Research",
+  "Operating System :: OS Independent",
+  "Development Status :: 4 - Beta",
+  "Topic :: Scientific/Engineering",
+]
+dependencies = [
+  "backscattering-simulation-data",
+  "finesse>=3.0a33",
+  "gitpython>=3.1.44",
+  "gwpy>=3.0.10",
+  "h5py>=3.14.0",
+  "rich>=14.0.0",
+]
+
+[tool.ruff]
+line-length = 72
+
+[tool.ruff.format]
+quote-style = "double"
+indent-style = "space"
+docstring-code-format = true
+
+[tool.basedpyright]
+typeCheckingMode = "recommended"
+allowedUntypedLibraries = ["gwpy", "h5py", "astropy", "finesse"]
+reportUnknownMemberType = false
+reportUnknownVariableType = false
+reportUnknownArgumentType = false
+reportAny = false
+
+[tool.uv.sources]
+backscattering-simulation-data = { git = "https://git.ligo.org/augustin.demagny/backscattering_simulation_data.git" }
+
+[dependency-groups]
+dev = ["basedpyright>=1.28.4", "ruff>=0.11.2"]

utils.py

@@ -0,0 +1,105 @@
+from typing import Any
+
+from finesse.model import Model, Port, SeriesSolution, DegreeOfFreedom
+from finesse.components import SignalGenerator
+from finesse.detectors import PowerDetectorDemod1
+from finesse.analysis.actions.axes import Xaxis, Noxaxis
+
+from numpy.typing import NDArray
+from numpy import float64
+
+
+def typingfix_seriessolution(
+    solution: None | SeriesSolution | tuple[Any],
+) -> SeriesSolution:
+    if solution is None:
+        raise ValueError("Result of the simulation is None")
+    if isinstance(solution, tuple):
+        raise ValueError("Result is not in the expected format")
+    return solution
+
+
+def fix_dark_fringe(model: Model, power: float) -> tuple[int, float]:
+    solution: SeriesSolution = typingfix_seriessolution(
+        model.run(Noxaxis())
+    )
+    if "B1_DC" not in solution.outputs:
+        raise Exception("Model does not contain B1 readout")
+
+    if type(model.get("DARM")) is DegreeOfFreedom:
+        result = solution["B1_DC"]
+        start, stop, nb = 0, 1, 0
+        while (abs(result - power) > 1e-6) and (nb < 500):
+            nb += 1
+            temp = start + (stop - start) / 2
+
+            model.get("DARM").DC = temp
+            solution: SeriesSolution = typingfix_seriessolution(
+                model.run(Noxaxis())
+            )
+            if "B1_DC" not in solution.outputs:
+                raise Exception("Model do not contains B1 readout")
+            result = solution["B1_DC"]
+            if result > power:
+                stop = temp
+            else:
+                start = temp
+        return nb, solution["B1_DC"]
+    raise Exception("DARM is not a degree of freedom in the model")
+
+
+def compute_TF(
+    model: Model,
+    frequencies: NDArray[float64],
+    inputs: list[Port | str],
+) -> SeriesSolution:
+    model.fsig.f = 1  # pyright: ignore[reportAttributeAccessIssue]
+
+    for input in inputs:
+        model.add(SignalGenerator("sig1", model.get(input), 1, 0))
+    model.add(
+        PowerDetectorDemod1(
+            "TF", model.get("SDB1").p2.o, f=model.get("fsig").f
+        )
+    )
+    model.get("TF").select_mask(exclude=(0, 0))
+    return model.run(
+        Xaxis(
+            model.fsig.f,  # pyright: ignore[reportAttributeAccessIssue]
+            "log",
+            min(frequencies),
+            max(frequencies),
+            frequencies.shape[0] - 1,
+        )
+    )
+
+
+class TF:
+    def __init__(
+        self,
+        model: Model,
+        inputs: list[Port | str],
+        outputs: list[Port | str],
+        frequencies: NDArray[float64],
+    ):
+        self.model = model.deepcopy()
+        self.inputs = inputs
+        self.outputs = outputs
+        self.frequencies = frequencies
+        self.result: SeriesSolution | None = None
+        self.compute()
+
+        # alias
+        self.f = self.frequencies
+
+    def compute(self):
+        solution = compute_TF(self.model, self.frequencies, self.inputs)
+        self.result = solution
+        return self
+
+    def get(self, output: Port | str):
+        if self.result is None:
+            raise Exception(
+                "TF need to be computed before getting result"
+            )
+        return self.result[output]