import numpy as np
from finesse.element import ModelElement
from finesse.parameter import float_parameter
from finesse.solutions.base import BaseSolution
[docs]@float_parameter("ASD", "Amplitude spectral density", units="ASD")
# IMPORTANT: renaming this class impacts the katscript spec and should be avoided!
class ClassicalNoise(ModelElement):
def __init__(self, name, node, ASD):
super().__init__(name)
self._add_to_model_namespace = True
self._namespace = (".noises",)
self.node = node
self.ASD = ASD
def _on_add(self, model):
if model is not self.node._model:
raise Exception(
f"{repr(self)} is using a node {self.node} from a different model"
)
[docs]class NoiseSolution(BaseSolution):
[docs] def __init__(self, name, f, output_nodes, noises, parents=None):
super().__init__(name, parents)
self.f = f
self.dtype = np.dtype([(n.name, float) for n in noises])
self.noises = {
onode: np.zeros(f.shape, dtype=self.dtype) for onode in output_nodes
}
[docs] def total_ASD(self, output_node):
return np.sqrt(sum(v**2 for v in self.noises[output_node].values()))
[docs] def plot(
self,
output_node,
noises=None,
show_total=True,
figsize_scale=1,
ax=None,
):
import matplotlib.pyplot as plt
if ax is None:
fig, ax = plt.subplots()
if not isinstance(output_node, str):
output_node = output_node.full_name
if output_node not in self.noises:
raise Exception(
f"The noise at {output_node} was not calculated, please recalculate."
)
for k, v in self.noises[output_node].items():
if noises is not None and k not in noises:
continue
if hasattr(v, "shape"):
ax.loglog(self.f, v, label=k)
else:
ax.loglog(self.f, v * np.ones_like(self.f), label=k)
if (noises is not None and len(noises) > 1) or (show_total and noises is None):
ax.loglog(
self.f, self.total_ASD(output_node), c="k", ls=":", lw=3, label="Total"
)
ax.legend()
ax.set_title(f"Noise at {output_node}")
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("ASD [UNIT/rtHz]")
plt.gcf().tight_layout()
return ax
[docs]def get_loop_network(model):
import networkx
import sympy
from finesse.components import NodeType, Wire
net = networkx.digraph.DiGraph()
def is_electric(x):
return all(p.type == NodeType.ELECTRICAL for p in x.ports)
optics_in = set() # nodes going into the optics
optics_out = set() # nodes coming out of the optics
Omega = sympy.var("Omega")
noise_nodes = tuple(n.node.full_name for n in model.noises)
for i, o, d in model.network.edges(data=True):
# Keep any edge going into an electronic element
if (
d["out_ref"]().type == NodeType.ELECTRICAL
or d["in_ref"]().type == NodeType.ELECTRICAL
):
if (
len(d["in_ref"]().connections) > 0
or i in noise_nodes
or o in noise_nodes
):
# Store the nodes that inject into an optical node
# and those optical nodes that feed into the electronics
if d["in_ref"]().type == NodeType.OPTICAL:
optics_out.add(d["out_ref"]())
elif d["out_ref"]().type == NodeType.OPTICAL:
optics_in.add(d["in_ref"]())
else:
el = d["owner"]()
sym = sympy.var(el.name)
if type(el) is Wire:
if not el.delay.is_changing:
sym = sympy.exp(-sympy.I * Omega * float(el.delay))
elif el.delay == 0:
sym = 1
net.add_edge(i, o, symbol=sym)
elif d["in_ref"]().type == NodeType.MECHANICAL and i in noise_nodes:
# Mechanical noise inputs need to be included so we can project them
# These are really just inputs into the optical system.
# I guess we could also have mechanical to electrical sensors...
optics_in.add(d["in_ref"]())
optic_couplings = {}
# Add in placeholds for any opto-mechanic coupling
for i in optics_in:
for o in optics_out:
sym = sympy.var(
f"{i.full_name.replace('.','')}__{o.full_name.replace('.','')}"
)
net.add_edge(i.full_name, o.full_name, symbol=sym)
optic_couplings[sym] = (i, o)
return net, optic_couplings
[docs]def nx_to_coo_sparse(G, nodelist=None):
from itertools import chain
if nodelist is None:
nodelist = list(G)
nlen = len(nodelist)
if nlen == 0:
raise Exception("Graph has no nodes or edges")
if len(nodelist) != len(set(nodelist)):
raise Exception("Ambiguous ordering: `nodelist` contained duplicates.")
index = dict(zip(nodelist, range(nlen)))
# build I-M sparse matrix
coefficients = zip(
*chain(
# Add in I matrix
((i, i, 1) for i in range(nlen)),
# negative off-diagonal elements
(
(index[u], index[v], -d["symbol"])
for u, v, d in G.edges(nodelist, data=True)
if u in index and v in index
),
)
)
try:
row, col, data = coefficients
except ValueError:
# there is no edge in the subgraph
row, col, data = [], [], []
return row, col, data
# class NoiseAnalysis(Action):
# def __init__(
# self,
# name,
# mode,
# start,
# stop,
# steps,
# output_node,
# *opt_output_nodes,
# noises=None,
# ):
# super().__init__(name)
# if mode not in ("lin", "log"):
# raise Exception("mode must be lin or log")
# if mode == "lin":
# self.f = np.linspace(start, stop, steps)
# else:
# self.f = np.geomspace(start, stop, steps)
# if noises is None:
# self.noises = tuple(n.name for n in output_node.component._model.noises)
# else:
# self.noises = tuple(n.name for n in noises)
# self.output_nodes = tuple(
# (output_node.full_name, *(o.full_name for o in opt_output_nodes))
# )
# self._info.parameters_changing = tuple(("fsig.f",),)
# self._info.makes_solution = True
# def fill_info(self, p_info):
# info = self.copy_info()
# p_info.add(info)
# def setup(self, s_prev, model: finesse.model.Model):
# import sympy
# from sympy.matrices import SparseMatrix
# from sympy import Matrix
# ws = NoiseAnalysisWorkspace(s_prev, model)
# ws.s_prev = s_prev
# ws.model = model
# ws.info = self.copy_info()
# ws.fn_do = self.do
# ws.params = tuple(
# finesse.analysis.actions.get_param(model, p)
# for p in self._info.parameters_changing
# )
# for p in ws.params:
# if not p.is_tunable:
# raise ParameterLocked(
# f"{repr(p)} must set as tunable " "before building the simulation"
# )
# # Setup the symbolic sparse matrix for solving noise propagations
# ws.loop_network, oc = get_loop_network(model)
# ws.noises = tuple(n for n in model.noises if n.name in self.noises)
# assert len(ws.noises) > 0
# # returns a I-M matrix for solving
# rcv = nx_to_coo_sparse(ws.loop_network)
# ws.M_nodes = nodes = tuple(ws.loop_network.nodes)
# M = SparseMatrix(
# None, {k: v for k, v in zip(tuple(zip(rcv[0], rcv[1])), rcv[2])}
# )
# # Always ensure that the Signal frequency variable
# # is present as we always fill that
# ws.Omega = sympy.var("Omega")
# ws.syms = M.free_symbols | set((ws.Omega,))
# ws.Minv = Matrix(M).inv()
# # Fill up a dict of how each noise
# # couples into each output node
# ws.TF = tuple([] for i in range(len(self.output_nodes)))
# for i, onode in enumerate(self.output_nodes):
# for j, noise in enumerate(ws.noises):
# noise_tf = ws.Minv[
# nodes.index(noise.node.full_name), nodes.index(onode)
# ].expand()
# ws.TF[i].append(
# sympy.lambdify(ws.syms, abs(noise_tf) ** 2, "numpy", dummify=False)
# )
# ws.values = {s: 0 for s in ws.syms}
# ws.OMEGA = ws.values[ws.Omega] = 2 * np.pi * self.f
# ws.lambdas = tuple(noise.ASD.lambdify(ws.model.fsig.f) for noise in ws.noises)
# # Get the simulations to do some calculations with in the do step
# ws.carrier = model.carrier_simulation
# try:
# ws.signal = model.signal_simulation
# ws.optical_couplings = {}
# ins = []
# outs = []
# # Store optical coulpings node indices for injecting later
# for key, (inode, onode) in oc.items():
# idx = ws.signal.field(inode, 0, 0)
# odx = ws.signal.field(onode, 0, 0)
# if idx in ins:
# x = ins.index(idx)
# else:
# ins.append(idx)
# x = len(ins) - 1
# if odx in outs:
# y = outs.index(odx)
# else:
# outs.append(odx)
# y = len(outs) - 1
# ws.optical_couplings[key] = (x, y)
# ws.input_rhs_indices = np.array((tuple(ins)), dtype=int)
# ws.output_rhs_indices = np.array((tuple(outs)), dtype=int)
# ws.out_fsig_sweep = np.zeros(
# (len(ins), len(self.f), len(outs)), dtype=np.complex128
# )
# except AttributeError:
# raise Exception(
# "Noise analysis requires a signal simulation to work, please set the model fsig.f frequency"
# )
# return ws
# def update_symbolic_values(self, ws):
# run_fsig_sweep(
# ws.carrier,
# ws.signal,
# self.f,
# ws.input_rhs_indices,
# ws.output_rhs_indices,
# ws.out_fsig_sweep,
# True,
# )
# for sym in ws.syms:
# if sym is ws.Omega:
# continue
# el = ws.model.elements.get(str(sym))
# if el:
# ws.values[sym] = el.eval(ws.OMEGA)
# else:
# # Now we loop over the actual simulation and run each point
# i, o = ws.optical_couplings[sym]
# ws.values[sym] = ws.out_fsig_sweep[i, :, o]
# def do(self, ws):
# ws.sol = NoiseSolution(
# self.name, self.f, self.output_nodes, ws.noises, ws.s_prev
# )
# ws.sol.Minv = ws.Minv
# ws.sol.Minv_nodes = ws.M_nodes
# self.update_symbolic_values(ws)
# ws.sol.optical_tfs = ws.out_fsig_sweep
# v = ws.values.values()
# for i, onode in enumerate(self.output_nodes):
# ws.sol.noises[onode] = {}
# for j, noise in enumerate(ws.noises):
# y_noise = ws.lambdas[j](self.f) * ws.TF[i][j](*v)
# ws.sol.noises[onode][noise.name] = y_noise