Commit 3118af93 authored by Daniel Brown's avatar Daniel Brown
Browse files

adding in longitudinal zpk suspension

parent c9dd962e
Pipeline #239816 passed with stages
in 27 minutes and 26 seconds
......@@ -10,7 +10,7 @@ from finesse.components.workspace cimport ConnectorWorkspace, FillFuncWrapper
from finesse.components.workspace import Connections
from finesse.cymath cimport complex_t
from finesse.components import Connector, NodeType, NodeDirection
from finesse.parameter import float_parameter
from finesse.parameter import float_parameter, info_parameter, bool_parameter
from finesse.components.general import DOFDefinition
from cpython.ref cimport PyObject, Py_XINCREF, Py_XDECREF
......@@ -440,3 +440,142 @@ class Pendulum(Connector):
) as mat:
omega0 = 2 * PI * ws.values.fpitch
mat[:] = 1 / ws.values.I_pitch * 1/(s**2 + s * omega0/ws.values.Qpitch + omega0**2)
class SuspensionZPKWorkspace(ConnectorWorkspace):
pass
@info_parameter("z", "Zeros")
@info_parameter("p", "Poles")
@float_parameter("k", "Gain")
@bool_parameter("enabled", "Enabled")
class SuspensionZPK(Connector):
"""A suspension that models multiple poles and zeros for the z motion of an optic.
ZPKs should describe the force to displacement transfer function. The user must
ensure that minus signs are correct for this transfer function as well as defining
complex conjugae pairs for physically correct behaviour.
ZPK terms are in units of radians.
The `k` gain parameter is a model parameter, so can be tuned and set to a symbol
during a simulation. Zeros and poles however are fixed during a simulation.
Parameters
----------
name : str
Element name
mech_port : Element or mechanical port
Mechanical port or element to attach this suspension to
z, p : array[float] or None
Zeros and poles of the transfer function
k : float
Gain of the transfer function
Example
-------
Free mass can be implemented using two poles at 0Hz and a gain of 1/Mass:
import matplotlib.pyplot as plt
model = finesse.Model()
model.parse('''
l l1
m m1 R=1 T=0
link(l1, m1)
fsig(1)
frequency_response(
geomspace(1, 100, 100),
l1.amp,
m1.mech.z
)
''')
model.add(SuspensionZPK('m1_sus', model.m1.mech, [], [0, 0], 1/10))
# Generate a bode plot for looking at the suspension transfer function
plt.figure()
model.m1_sus.bode_plot(label='10kg')
model.m1_sus.k = 1/100
model.m1_sus.bode_plot(label='100kg')
plt.legend()
# Compute the transfer function from laser amplitude to mirror motion
plt.figure()
out = model.run()
out.plot();
"""
def __init__(self, name, mech_port, z, p, k, enabled=True):
super().__init__(name)
self.z = z if z is not None else []
self.p = p if p is not None else []
self.k = k
self.enabled = enabled
# Add motion and force nodes to mech port.
# Here we duplicate the already created mechanical
# nodes in some other connector element
self._add_port("mech", NodeType.MECHANICAL)
self.mech._add_node("z", None, mech_port.z)
self.mech._add_node("F_z", None, mech_port.F_z)
# We just have direct coupling between DOF, no cross-couplings
self._register_node_coupling(
"F_to_Z", self.mech.F_z, self.mech.z
)
def _get_workspace(self, sim):
if sim.signal and (self.enabled or self.enabled.is_changing):
refill = (
sim.model.fsig.f.is_changing or self.enabled.is_changing
)
ws = SuspensionZPKWorkspace(self, sim)
ws.signal.add_fill_function(self.fill, refill)
ws.z = np.array(self.z, copy=True)
ws.p = np.array(self.p, copy=True)
ws.k = self.k
return ws
else:
return None
def bode_plot(self, f=None, Hz=True, *args, **kwargs):
"""Plot bode function for this suspensions force to displacement
transfer function. When Hz is True the plot will be shown in units
of Hertz. The input frequency f is always in units of Hertz.
*args and **kwargs are passed to the Maptlotlib plotting function.
"""
from scipy import signal
import matplotlib.pyplot as plt
fig = plt.gcf()
Nax = len(fig.axes)
axs = fig.axes
if Nax == 0:
fig.subplots(2, 1)
axs = fig.axes
elif Nax != 2:
raise RuntimeError("Figure does not have two axes to plot magnitude and phase on")
w, mag, phase = signal.bode((self.z, self.p, self.k), w=2*np.pi*f if f is not None else None)
if Hz:
w /= 2*np.pi
axs[0].semilogx(w, mag, *args, **kwargs)
axs[1].semilogx(w, phase, *args, **kwargs)
if Hz:
axs[1].set_xlabel("Frequency [Hz]")
else:
axs[1].set_xlabel("Frequency [Radians]")
axs[1].set_ylabel("Phase [Hz]")
axs[0].set_ylabel("Magnitude [m/N (dB)]")
def fill(self, ws):
s = 2j* np.pi * ws.sim.model_data.fsig
if ws.values.enabled:
H = ws.k * np.prod(s - ws.z)/np.prod(s - ws.p)
else:
H = 0
if ws.signal.connections.F_to_Z_idx >= 0:
with ws.sim.signal.component_edge_fill3(
ws.owner_id, ws.signal.connections.F_to_Z_idx, 0, 0,
) as mat:
mat[:] = H
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment