import numpy as np
from more_itertools import pairwise
from scipy.sparse import diags
from finesse.components import Beamsplitter, Mirror
from finesse.exceptions import FinesseException
from ...solutions import BaseSolution
from . import Operator
from .base import Action
[docs]class PseudoLockCavitySolution(BaseSolution):
"""Solution to the :class:`.PseudoLockCavity` action.
Attributes
----------
operator : :class:`.OperatorSolution`
Solution of the operator action that is used to compute the roundtrip
operator in the cavity
w : array_like
2D Eigenvector array
v : array_like
1D Eigenvalue array
v_idx : int
Index of the eigenmode being locked to
phase : float
Phase in degrees of the eigenmode that is being locked to. This is what
is used to set the operating point of the cavity
cavity : str
Which cavity is being locked
lowest_loss : bool
If False, the target `mode` attribute is used. Otherwise
mode : array_like
[n, m] target mode to try and lock to
"""
# IMPORTANT: renaming this class impacts the katscript spec and should be avoided!
[docs]class PseudoLockCavity(Action):
"""An action that locks a cavity defined by a `Cavity` element to a specific mode
without using any radio-frequency sensing scheme. This will only work on simple
cavities that are not coupled in any way. You can specify whether to try and lock to
a particular HG mode with the `mode=[n,m]` keyword argument, or just pick the lowest
loss mode, `lowest_loss=True`.
Parameters
----------
cavity : Cavity
Cavity element describing some Fabry-Perot like optical cavity
mode : (n, m), optional
HG mode to try and lock to, default is [0,0]
lowest_loss : bool, optional
Select the eigenmode which has the lowest loss, most likely the fundamental
mode of the cavity. Using lowest loss will override the `mode` selection.
feedback : Parameter optional
If `None` the required cavity tuning to lock to the calculated mode will be
determined from the cavity objects source node element, and the relevant
`phi` parameter will be used. Alternatively you can specify which tuning
parameter is used instead. Which should be a `phi` of some mirror in the
cavity or a `DegreeOfFreedom` which controls the cavity length.
name : str, optional
Name of the solution generated by this action
Examples
--------
A Fabry-Perot based on aLIGO cavities:
>>> import finesse
>>> model = finesse.Model()
>>> model.parse('''
... l l1 P=1
... m ITM T=0.014 L=0 Rc=-1945
... s sARM ITM.p2 ETM.p1 L=3994
... m ETM R=1 L=0 Rc=2145
... link(l1, ITM)
... cav arm ETM.p1.o
... modes(maxtem=1)
... ''')
The cavity source node is on the ETM so the pseudo-lock will use that node and
component `phi` parameter to feedback to. In this case the ETM.phi will be
corrected to match the ITM.phi which will make the HG00 resonant.
>>> model.ITM.phi = 0
>>> model.ETM.phi = 10
>>> model.run("pseudo_lock_cavity(arm)")
>>> print(model.ITM.phi, model.ETM.phi)
0.0 degrees 0.0 degrees
>>> model.ITM.phi = 10
>>> model.ETM.phi = 0
>>> model.run("pseudo_lock_cavity(arm)")
>>> print(model.ITM.phi, model.ETM.phi)
10.0 degrees 10.0 degrees
This lock will also handle any misalignments, mismatches, or maps applied to
the cavity.
>>> model.ITM.xbeta = 5e-8
>>> model.ITM.phi = 11
>>> model.ETM.phi = 0
>>> model.run("pseudo_lock_cavity(arm)")
>>> print(model.ITM.phi, model.ETM.phi)
11.0 degrees 10.984168676865762 degrees
We can also lock to other HG modes:
>>> sol = model.run("series(pseudo_lock_cavity(arm, mode=[1,0]), noxaxis())")
>>> print(model.ITM.phi, model.ETM.phi)
>>> print(model.homs)
>>> print(sol['noxaxis']['Ecirc'])
11.0 degrees 28.663091273523214 degrees
[[0 0]
[1 0]
[0 1]]
[-0.16465261-0.11468501j 0.41262682-0.27112167j 0. +0.j ]
"""
def __init__(
self,
cavity,
*,
mode=None,
lowest_loss=False,
feedback=None,
name="pseudo_lock_cavity",
):
super().__init__(name)
self.cavity = cavity
self.mode = [0, 0] if mode is None else mode
self.lowest_loss = lowest_loss
self.feedback = feedback
def _requests(self, model, memo, first=True):
cav = (
model.elements[self.cavity] if isinstance(self.cavity, str) else self.cavity
)
memo["output_nodes"].extend((_.full_name for _ in cav.path.nodes))
if not any(np.all(model.homs == self.mode, axis=1)):
raise FinesseException(
f"HG mode {self.mode} was not found in the simulation"
)
# If a user provided feedback is given then use that to send locking
# signal towards
if self.feedback is not None:
if isinstance(self.feedback, str):
memo["changing_parameters"].append(model.get(self.feedback))
elif hasattr(self.feedback, "parameter"):
memo["changing_parameters"].append(self.feedback.parameter)
else:
memo["changing_parameters"].append(self.feedback)
else:
# if not use the cavity source component
memo["changing_parameters"].append(cav.source.component.phi)
def _do(self, state):
model = state.model
key = id(self)
sol = PseudoLockCavitySolution(self.name)
sol.mode = self.mode
sol.lowest_loss = self.lowest_loss
sol.cavity = str(self.cavity)
# Try and get data already computed for this action
ws = state.action_workspaces.get(key, None)
if ws is None:
state.action_workspaces[key] = ws = {}
cav = ws["cavity"] = (
model.elements[self.cavity]
if isinstance(self.cavity, str)
else self.cavity
)
# if the cavity is using a via definition for it's path use that
# if not use the opposite to ensure we don't select a dumb option
via = cav.source.opposite if cav.via is None else cav.via
ws["operator"] = Operator(cav.source, cav.source, via=via)
ws["hg_idx"] = int(np.argmax(np.all(model.homs == self.mode, axis=1)))
# Get the phi parameter of the cavity source component
component = cav.source.component
if self.feedback is None:
ws["phi"] = component.phi
elif isinstance(self.feedback, str):
ws["phi"] = model.get(self.feedback)
elif hasattr(self.feedback, "parameter"):
ws["phi"] = self.feedback.parameter
else:
ws["phi"] = self.feedback
# Decide on the direction of phase application, depending on side
# of mirror/bs the cav source node is on
if isinstance(component, Mirror):
if cav.source in component.p1:
ws["scale"] = -1
else:
ws["scale"] = +1
elif isinstance(component, Beamsplitter):
# TODO need angle of incidence here
if cav.source in component.p1 or cav.source in component.p2:
ws["scale"] = -1
else:
ws["scale"] = +1
else:
raise FinesseException(
f"From the cavity node {cav.source}, its owner {cav.source.component} is neither a mirror nor a beamsplitter."
)
# Compute eigenmodes and vectors
sol.operator = state.apply(ws["operator"])
sol.w, sol.v = np.linalg.eig(sol.operator.operator)
if not self.lowest_loss:
# select mode with largest self.mode component
# which is probably the best guess for making that
# mode resonant
sol.v_idx = np.argmax(abs(sol.v[:, ws["hg_idx"]]))
else:
# select eigenmode which has the lowest loss, likely the fundamental
# minimize -> roundtip loss = 1 - abs(sol.w)**2
sol.v_idx = np.argmax(abs(sol.w))
sol.phase = np.angle(sol.w[sol.v_idx], deg=True)
ws["phi"].value += ws["scale"] * sol.phase / 2
return sol
[docs]class PseudoLockDRFPMISolution(BaseSolution):
"""TODO: needs updating for all the various attributes this stores."""
[docs] def plot_PRC_SRC_eigenvectors(self):
import matplotlib.pyplot as plt
"""Plots eigenvector maxtrix for PRC and SRC"""
plt.imshow(abs(self.M_evec), origin="upper")
plt.text(0, 0, "PRC", c="b", horizontalalignment="center")
plt.text(0, 3, "SRC", c="b", horizontalalignment="center")
plt.hlines(2.5, -0.5, 5.5)
plt.colorbar()
plt.title("Eigenvectors")
# IMPORTANT: renaming this class impacts the katscript spec and should be avoided!
[docs]class PseudoLockDRFPMI(Action):
"""Pseudo-locking is attempting to find an operating point for a LIGO like model
without needing to use RF sidebands and readouts. Although it is not physically
accurate it does provide a useful tool for analysing detectors from a more
theoretical basis. This generates a :class:`.PseudoLockDRFPMISolution` solution
containing various operators and results.
This action is hardcoded to work with a LIGO like model. Mirrors should be named
with ITMX, ETMX, PRM, etc.
This currently action only really works for finding the lock points for the PRC,
SRC, XARM, and YARM – which all have defined cavity roundtrips which allow
eigendecomposition of roundtrip operators. The eigenvectors describe the HOM mix
for each resonant mode in a cavity and the eigenvalues the roundtrip phase and
loss of the mode. This code looks for eigenvectors with the largest HG00 content
and then uses the eigenvalues to compute what cavity tunings need to make this
mode resonant.
The corner is the most complicated here as the PRC and SRC are coupled via the
beamsplitter. The eigendecomposition is performed on the 2x2 operator matrix for
the PRC, SRC, and the coupling matrices between them. When the coupling is small
the results are the same as performing the decomposition on each SRC and PRC
separately.
Note that the only degree of freedom that this does not handle currently is MICH.
MICH is awkward because it is not a cavity. It is essentially the beamsplitter
position that makes the anti-symmetric port dark.
Parameters
----------
frequency : float, optional
Frequency to use for calculating the operators
apply_tunings : bool, optional
When True the action will modify the model tunings
name : str, optional
Name of the solution generated by this action
"""
def __init__(self, frequency=0, *, apply_tunings=True, name="operator_lock"):
super().__init__(name)
self.frequency = frequency
self.apply_tunings = apply_tunings
# connections to compute operators for in a LIGO like model
# name : (from, to, via)
self.operators = {
"BS13": ("BS.p1.i", "BS.p3.o", None),
"BS31": ("BS.p3.i", "BS.p1.o", None),
"BS12": ("BS.p1.i", "BS.p2.o", None),
"BS21": ("BS.p2.i", "BS.p1.o", None),
"BS24": ("BS.p2.i", "BS.p4.o", None),
"BS42": ("BS.p4.i", "BS.p2.o", None),
"BS34": ("BS.p3.i", "BS.p4.o", None),
"BS43": ("BS.p4.i", "BS.p3.o", None),
"MICHX1": ("BS.p3.o", "ITMX.p1.i", None),
"MICHY1": ("BS.p2.o", "ITMY.p1.i", None),
"MICHX2": ("ITMX.p1.o", "BS.p3.i", None),
"MICHY2": ("ITMY.p1.o", "BS.p2.i", None),
"R1X": ("ITMX.p1.i", "ITMX.p1.o", None),
"R2X": ("ITMX.p2.i", "ITMX.p2.o", None),
"T1X": ("ITMX.p1.i", "ITMX.p2.o", None),
"T2X": ("ITMX.p2.i", "ITMX.p1.o", None),
"GX": ("ITMX.p2.o", "ITMX.p2.i", "ETMX.p1.i"),
"GY": ("ITMY.p2.o", "ITMY.p2.i", "ETMY.p1.i"),
"R1Y": ("ITMY.p1.i", "ITMY.p1.o", None),
"R2Y": ("ITMY.p2.i", "ITMY.p2.o", None),
"T1Y": ("ITMY.p1.i", "ITMY.p2.o", None),
"T2Y": ("ITMY.p2.i", "ITMY.p1.o", None),
"SRC": ("BS.p4.o", "BS.p4.i", "SRM.p1.o"),
"PRC": ("BS.p1.o", "BS.p1.i", "PRM.p2.o"),
"SQZ": ("sqz.p1.o", "SRM.p2.o", None),
}
def _requests(self, model, memo, first=True):
for _from, _to, via in self.operators.values():
memo["output_nodes"].append(_from)
memo["output_nodes"].append(_to)
if via:
memo["output_nodes"].append(via)
if self.apply_tunings:
memo["changing_parameters"].append("SRCL.DC")
memo["changing_parameters"].append("PRCL.DC")
def _do(self, state):
sim = state.sim
model = state.model
key = id(self)
f_idx = None
N = sim.model_settings.num_HOMs
sol = PseudoLockDRFPMISolution(self.name)
# Try and get data already computed for this action
ws = state.action_workspaces.get(key, None)
if ws is None:
try:
frequency = float(self.frequency)
except TypeError:
frequency = float(self.frequency.value)
# find the right frequency index
for freq in state.sim.carrier.optical_frequencies.frequencies:
if freq.f == frequency:
f_idx = freq.index
break
if f_idx is None:
raise RuntimeError(
f"Could not find an optical carrier frequency with a value of {self.frequency}Hz"
)
state.action_workspaces[key] = ws = {}
# loop through each operator we need to get
for name, (start, end, via) in self.operators.items():
nodes = model.path(start, end, via_node=via).nodes
connections = tuple(
(n1.full_name, n2.full_name) for n1, n2 in pairwise(nodes)
)
operators = []
# Get all the matrix views to eventually multiply together later
for _ in pairwise(nodes):
# if we have a 1D array it's just a diagonal matrix
# so convert it to a sparse array for easy multiplying
# later
# TODO can have some option to select between different
# frequencies or signal/carrier at some point
Mview = sim.carrier.connections[_]
if Mview.ndim == 2:
# We have a frequency scattering matrix
# TODO : not sure on user interface for getting
# different frequency couplings yet, for now it's
# just same freq in and out
M = Mview[f_idx, f_idx]
else:
# no frequency coupling
M = Mview[f_idx]
operators.append(M)
ws[name] = (nodes, connections, operators)
# Loop through all the opeators are multiply them together to get
# each reduced operator
sol.operators = {}
for name, (_, _, operators) in ws.items():
M = np.eye(sim.model_settings.num_HOMs)
for op in operators:
if op.view.ndim == 1:
M = diags(-op.view) @ M
else:
M = -op.view @ M
sol.operators[name] = M
# Give all the opertors variable names
GX = sol.operators["GX"]
GY = sol.operators["GY"]
SRC = sol.operators["SRC"]
PRC = sol.operators["PRC"]
R1X = sol.operators["R1X"]
R1Y = sol.operators["R1Y"]
R2X = sol.operators["R2X"]
R2Y = sol.operators["R2Y"]
# T1X = sol.operators["T1X"]
# T1Y = sol.operators["T1Y"]
# T2X = sol.operators["T2X"]
# T2Y = sol.operators["T2Y"]
BS12 = sol.operators["BS12"]
BS21 = sol.operators["BS21"]
BS13 = sol.operators["BS13"]
BS31 = sol.operators["BS31"]
BS24 = sol.operators["BS24"]
BS42 = sol.operators["BS42"]
BS34 = sol.operators["BS34"]
BS43 = sol.operators["BS43"]
MICHX1 = sol.operators["MICHX1"]
MICHX2 = sol.operators["MICHX2"]
MICHY1 = sol.operators["MICHY1"]
MICHY2 = sol.operators["MICHY2"]
# Now finally multiply the various operators together to compute
# the different eigendecompositions
# I = np.eye(N, dtype=np.complex128)
GRTX = GX @ R2X
GRTY = GY @ R2Y
RX = R1X # + T2X @ np.linalg.inv(I - GRTX) @ GX @ T1X
RY = R1Y # + T2Y @ np.linalg.inv(I - GRTY) @ GY @ T1Y
MX = MICHX2 @ RX @ MICHX1
MY = MICHY2 @ RY @ MICHY1
MICHSS = BS34 @ MX @ BS43 + BS24 @ MY @ BS42
MICHSP = BS31 @ MX @ BS43 + BS21 @ MY @ BS42
MICHPP = BS31 @ MX @ BS13 + BS21 @ MY @ BS12
MICHPS = BS34 @ MX @ BS13 + BS24 @ MY @ BS12
sol.SRC = MICHSS @ SRC
sol.PRC = MICHPP @ PRC
sol.SRC_PRC = MICHSP @ SRC
sol.PRC_SRC = MICHPS @ PRC
# eigenvalues and vectors of the SRC roundtrip operator gives
# the various resonant mode structure regardless of HOM content
# and defects
sol.X_eval, sol.X_evec = np.linalg.eig(GRTX)
sol.Y_eval, sol.Y_evec = np.linalg.eig(GRTY)
sol.X_HG00_idx, sol.Y_HG00_idx = (
abs(sol.X_evec[0, :]).argmax(),
abs(sol.Y_evec[0, :]).argmax(),
)
sol.xarm_phi = -np.angle(sol.X_eval[sol.X_HG00_idx], deg=True) / 2
sol.yarm_phi = -np.angle(sol.Y_eval[sol.Y_HG00_idx], deg=True) / 2
# Make a 2x2 block matrix for corner coupling between SRC and PRC
sol.M = np.block([[sol.PRC, sol.PRC_SRC], [sol.SRC_PRC, sol.SRC]])
sol.M_eval, sol.M_evec = np.linalg.eig(sol.M)
# Compute which eigenvectors are the most HG00 like
# Grab the HG00 row in and select the column with the
# max value, will be the mode with the most HG00 content
sol.PRC_HG00_idx, sol.SRC_HG00_idx = (
abs(sol.M_evec[0, :]).argmax(),
abs(sol.M_evec[N, :]).argmax(),
)
# Compute phases for each cavity
sol.corner_phases = -np.angle(sol.M_eval, deg=True) / 2
# Then select the cavity detuning to pick the most HG00 mode
# Resonant sideband extraction
sol.src_phi = sol.corner_phases[sol.SRC_HG00_idx]
sol.prc_phi = sol.corner_phases[sol.PRC_HG00_idx]
sol.PRC_evec = sol.M_evec[:N, sol.PRC_HG00_idx]
sol.SRC_evec = sol.M_evec[N:, sol.SRC_HG00_idx]
if self.apply_tunings:
state.model.SRCL.DC += sol.src_phi
state.model.PRCL.DC += 90 - sol.prc_phi
return sol