Merge branch 'master' into feature/traceorder

docs/source/usage/homs/index.rst

@@ -28,3 +28,4 @@ The pages below focus on the non-simulation aspects of HOM modelling, e.g. using
     :maxdepth: 2
 
     propagating_beams
+    scattering_matrices

docs/source/usage/homs/scattering_matrices.rst

+.. include:: /defs.hrst
+
+.. _arbitrary_scatter_matrices:
+
+Computing arbitrary scattering matrices
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The internal functions for computing coupling coefficient matrices
+during a simulation are also exposed via specific methods in the
+:mod:`finesse.knm` sub-module. The most-general function for computing
+these matrices is :func:`.make_scatter_matrix` - this requires the matrix
+type as the first argument, mapping to one of the specific scattering
+matrix functions in this module.
+
+As these Python-facing functions internally call the same C functions used
+to compute coupling coefficients during simulations, you can rest assured
+that these routines are highly optimised - allowing for fast computation of
+these matrices.
+
+Calculating coupling matrices via the Bayer-Helms formalism
+```````````````````````````````````````````````````````````
+
+To compute a scattering matrix based on a system with mode-mismatches and
+tilt (angular) misalignments, one can use :func:`.make_bayerhelms_matrix` (or
+:func:`.make_scatter_matrix` with "bayerhelms" as the type, as noted above).
+
+This function requires the input and output beam parameters, in both tangential
+and sagittal planes, as well as the misalignment angles in both planes. One can
+optionally provide a medium refractive index (defaults to :math:`n_r = 1`) and /
+or a beam wavelength (defaults to :math:`\lambda = 1064\,\mathrm{nm}`) too.
+
+As an example, below, we compute a scattering matrix with only mismatches present,
+up to a `maxtem` of 5 (see :ref:`selecting_modes`):
+
+.. jupyter-execute::
+
+    import finesse
+    finesse.configure(plotting=True)
+
+    from finesse.knm import make_scatter_matrix
+
+    # Non-astigmatic beam with 20% mismatch in w0
+    # from input to output
+    qx1 = qy1 = finesse.BeamParam(w0=1e-3, z=0)
+    qx2 = qy2 = finesse.BeamParam(w0=1.2e-3, z=0)
+
+    # No angular misalignments
+    xgamma = ygamma = 0
+
+    # Compute the scattering matrix, returning a KnmMatrix object
+    kmat = make_scatter_matrix(
+        "bayerhelms",
+        qx1, qx2, qy1, qy2,
+        xgamma, ygamma,
+        maxtem=5,
+    )
+
+    # Plot the absolute value of the coupling
+    # coefficients as a colormesh
+    kmat.plot(cmap="bone");
+
+As expected, we get a plot showing non-zero couplings for even mode orders. Note
+that the ``kmat`` object here will be a :class:`.KnmMatrix` object, providing
+a convenient interface for accessing specific coupling coefficients, printing the
+matrix and plotting the matrix. The :meth:`.KnmMatrix.plot` method computes and
+plots the absolute value of the coefficients by default; this can be altered as
+shown in the documentation for this method.
+
+This function also accepts mode selection keyword arguments, allowing for custom
+mode indices to be used (see :ref:`selecting_modes` and :func:`.make_modes`). For
+example, in this case we could select just even modes as we know that no angular
+misalignments are present:
+
+.. jupyter-execute::
+
+    kmat = make_scatter_matrix(
+        "bayerhelms",
+        qx1, qx2, qy1, qy2,
+        xgamma, ygamma,
+        select="even", maxtem=8,
+    )
+    kmat.plot(cmap="bone");
+
+.. rubric:: Accessing specific couplings
+
+As implied previously, the return type of :func:`.make_scatter_matrix` is a
+:class:`.KnmMatrix` object. Along with plotting and stringifying, this object
+provides a dict-like interface for accessing coupling coefficients in the matrix.
+
+Here are a few convenient ways of accessing these coefficients:
+
+.. jupyter-execute::
+
+    # Using the "nm->n'm'" format
+    print("Coupling coeff. from 00 -> 02 = ", kmat["00->02"])
+
+    # Using the "nmn'm'" format
+    print("Coupling coeff. from 02 -> 40 = ", kmat["0240"])
+
+    # Using the (n, m, n', m') format
+    print("Coupling coeff. from 26 -> 44 = ", kmat[(2, 6, 4, 4)])

docs/source/usage/homs/selecting_modes.rst

 .. include:: /defs.hrst
 
+.. _selecting_modes:
+
 Selecting the modes to model
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

setup.py

@@ -146,6 +146,7 @@ def ext_modules():
         ("knm.pyx", open_mp_args),
         ("components/workspace.pyx", {}),
         ("components/mechanical.pyx", {}),
+        ("components/modal/workspace.pyx", {}),
         ("components/modal/mirror.pyx", {}),
         ("components/modal/beamsplitter.pyx", {}),
         ("components/modal/signal.pyx", {}),

src/finesse/__main__.py

@@ -242,6 +242,18 @@ def plot_graph(state, network, layout, graphviz=False):
 
 input_file_argument = click.argument("input_file", type=click.File("r"))
 output_file_argument = click.argument("output_file", type=click.File("w"), default="-")
+plot_option = click.option(
+    "--plot/--no-plot",
+    default=True,
+    show_default=True,
+    help="Display results as figure.",
+)
+trace_option = click.option(
+    "--trace/--no-trace",
+    default=False,
+    show_default=True,
+    help="Displays the results of a beam trace of the model.",
+)
 graphviz_option = click.option(
     "--graphviz",
     is_flag=True,
@@ -322,18 +334,8 @@ def cli(ctx):
 
 @cli.command()
 @input_file_argument
-@click.option(
-    "--plot/--no-plot",
-    default=True,
-    show_default=True,
-    help="Display results as figure.",
-)
-@click.option(
-    "--trace/--no-trace",
-    default=False,
-    show_default=True,
-    help="Displays the results of a beam trace of the model.",
-)
+@plot_option
+@trace_option
 @legacy_option
 @verbose_option
 @quiet_option
@@ -714,6 +716,37 @@ def debug(ctx):
     set_debug(ctx, None, True)
 
 
+@debug.command(name="run")
+@input_file_argument
+@plot_option
+@trace_option
+@legacy_option
+@verbose_option
+@quiet_option
+@debug_option
+@click.pass_context
+def run_debug(
+    ctx, input_file, plot, trace, legacy,
+):
+    """Run a Finesse script with Python fault handler enabled.
+
+    If the file extension is 'py', it is interpreted as Python code; otherwise it is parsed assuming
+    it to be kat script.
+    """
+    import os
+    import faulthandler
+
+    click.secho("Enabling Python fault handler", fg="yellow")
+    faulthandler.enable()
+
+    root, ext = os.path.splitext(input_file.name)
+    if ext.casefold() == ".py":
+        exec(input_file.read())
+    else:
+        # Forward arguments to normal "run" command.
+        ctx.forward(run)
+
+
 @debug.command()
 @input_file_argument
 @graph_layout_argument

src/finesse/components/beamsplitter.py

@@ -8,12 +8,10 @@ import numpy as np
 
 from finesse.exceptions import TotalReflectionError
 from finesse.parameter import model_parameter, Rebuild
-from finesse.knm import KnmMatrix
 from finesse.utilities import refractive_index
 
-from finesse.components.general import InteractionType, _conn_abcd_return
+from finesse.components.general import InteractionType
 from finesse.components.surface import Surface
-
 from finesse.components.node import NodeDirection, NodeType
 
 LOGGER = logging.getLogger(__name__)
@@ -108,13 +106,6 @@ class Beamsplitter(Surface):
         mass=np.inf,
         signal_gain=1e-9,
     ):
-        from finesse.components.modal.beamsplitter import BeamsplitterConnections
-
-        # NOTE this is needed (for now) as it is used in _check_phi which is currently called
-        #      in Surface constructor when initialising phi value
-        self._cos_alpha = 1
-        self._cos_alpha_2 = 1
-
         super().__init__(
             name,
             R,
@@ -126,8 +117,9 @@ class Beamsplitter(Surface):
             ybeta,
             mass,
             signal_gain,
-            BeamsplitterConnections,
         )
+        self._cos_alpha = 1
+        self._cos_alpha_2 = 1
 
         self._add_port("p1", NodeType.OPTICAL)
         self.p1._add_node("i", NodeDirection.INPUT)
@@ -241,7 +233,7 @@ class Beamsplitter(Surface):
 
         nr1 = refractive_index(from_node, symbolic=True) or 1
         nr2 = refractive_index(to_node, symbolic=True) or 1
-
+        
         alpha1 = np.radians(alpha)
         # we get alpha2 from Snell's law
         if from_node.port is self.p3 or from_node.port is self.p4:
@@ -613,12 +605,6 @@ class Beamsplitter(Surface):
         ws.nr2 = refractive_index(self.p3) or refractive_index(self.p4) or 1
 
         ws.set_fill_fn(beamsplitter_fill)
-        # Only allocate new memory for scattering matrices for the carrier
-        # simulation. Audio simulation will use the carrier Knm matrices.
-        if not sim.is_audio:
-            self.__allocate_knm_matrices(ws)
-        else:
-            self.__set_audio_ws_knm_matrix_refs(ws)
 
         # Initialise the ABCD matrix memory-views
         if sim.is_modal:
@@ -685,33 +671,3 @@ class Beamsplitter(Surface):
                 # TODO, can't access matrix like this!
                 sim.Mq[sim.field(self.p1.o, freq.index, 0)] = self.L / 2
                 sim.Mq[sim.field(self.p2.o, freq.index, 0)] = self.L / 2
-
-    def __allocate_knm_matrices(self, ws):
-        Is = np.eye(ws.sim.nhoms, dtype=np.complex128)
-        losses = np.ones(ws.sim.nhoms)
-
-        couplings = ["12", "21", "34", "43", "13", "31", "24", "42"]
-        for c in couplings:
-            setattr(ws, f"K{c}", KnmMatrix(Is.copy(), self, kdir=c))
-            setattr(ws, f"K{c}_loss", losses.copy())
-
-    def __set_audio_ws_knm_matrix_refs(self, ws):
-        DC = ws.sim.DC
-        DC_workspaces = list(filter(lambda x: x.owner is self, DC.workspaces))
-        if not DC_workspaces:
-            raise RuntimeError(
-                "Bug encountered! Could not find a DC simulation workspace "
-                f"associated with Beamsplitter of name {self.name}"
-            )
-        if len(DC_workspaces) > 1:
-            raise RuntimeError(
-                "Bug encountered! Found multiple workspaces in the DC simulation "
-                f"associated with Beamsplitter of name {self.name}"
-            )
-
-        DC_ws = DC_workspaces[0]
-
-        couplings = ["12", "21", "34", "43", "13", "31", "24", "42"]
-        for c in couplings:
-            setattr(ws, f"K{c}", getattr(DC_ws, f"K{c}"))
-            setattr(ws, f"K{c}_loss", getattr(DC_ws, f"K{c}_loss"))

src/finesse/components/general.py

@@ -179,12 +179,11 @@ class Connector(ModelElement):
 
     _borrows_nodes = False
 
-    def __init__(self, name, connections_type=None):
+    def __init__(self, name):
         from finesse.components.workspace import Connections
 
         super().__init__(name)
 
-        self.__connections_type = connections_type or Connections
         self.__connections = OrderedDict()
         self.__ports = OrderedDict()
         self.__interaction_types = {}
@@ -202,10 +201,6 @@ class Connector(ModelElement):
     def __nodes_of(self, node_type):
         return tuple([node for node in self.nodes.values() if node.type == node_type])
 
-    @property
-    def connections_type(self):
-        return self.__connections_type
-
     @property
     def optical_nodes(self):
         """The optical nodes stored by the connector.
@@ -455,7 +450,7 @@ class Connector(ModelElement):
             to_node = to_node.o
 
         return from_node, to_node
-    
+
     def _resymbolise_ABCDs(self):
         # By default components will not have to resymbolise optical ABCD matrices
         pass
@@ -468,8 +463,7 @@ class Connector(ModelElement):
         symbolic=False,
         copy=True,
         retboth=False,
-    ):  # This docstring is appended too in inheriting classes with more info
-        # Don't add before Parameters section
+    ):
         """
         Parameters
         ----------
@@ -529,7 +523,7 @@ class Connector(ModelElement):
         else:
             M_sym, M_num = self._abcd_matrices[key]
 
-            if M_sym is None : 
+            if M_sym is None:
                 # Symbolic M is None if an error occurred during the tracing
                 # In such a case M_num is a TotalReflectionError or
                 # some other exception instance
@@ -547,12 +541,12 @@ def _conn_abcd_return(M_sym, M_num, symbolic: bool, copy: bool, retboth: bool):
         Ms = M_sym
         Mn = M_num
 
-    if symbolic:
-        return Ms
-
     if retboth:
         return Ms, Mn
 
+    if symbolic:
+        return Ms
+
     return Mn
 
 

src/finesse/components/isolator.py

@@ -6,11 +6,8 @@ import logging
 
 import numpy as np
 
-from finesse.knm import KnmMatrix
-
 from finesse.components.modal.isolator import (
-    IsolatorWorkspace,
-    IsolatorConnections,
+    IsolatorWorkspace
 )
 from finesse.components.general import Connector, InteractionType
 from finesse.components.node import NodeDirection, NodeType
@@ -34,7 +31,7 @@ class Isolator(Connector):
     """
 
     def __init__(self, name, S=0.0):
-        super().__init__(name, IsolatorConnections)
+        super().__init__(name)
 
         self.S = S
 
@@ -69,13 +66,6 @@ class Isolator(Connector):
         ws.nr1 = refractive_index(self.p1) or 1
         ws.nr2 = refractive_index(self.p2) or 1
 
-        # Only allocate new memory for scattering matrices for the carrier
-        # simulation. Audio simulation will use the carrier Knm matrices.
-        if not sim.is_audio:
-            self.__allocate_knm_matrices(ws)
-        else:
-            self.__set_audio_ws_knm_matrix_refs(ws)
-
         ws.set_fill_fn(self._fill_matrix)
         return ws
 
@@ -93,28 +83,3 @@ class Isolator(Connector):
                 ws.owner_id, ws.connections.P2i_P1o_idx, freq.index, freq.index,
             ) as mat:
                 np.multiply(suppression, ws.K21.data, out=mat[:])
-
-    def __allocate_knm_matrices(self, ws):
-        Is = np.eye(ws.sim.nhoms, dtype=np.complex128)
-        ws.K12 = KnmMatrix(Is.copy(), self, kdir="12")
-        ws.K21 = KnmMatrix(Is.copy(), self, kdir="21")
-
-    def __set_audio_ws_knm_matrix_refs(self, ws):
-        DC = ws.sim.DC
-        DC_workspaces = list(filter(lambda x: x.owner is self, DC.workspaces))
-        if not DC_workspaces:
-            raise RuntimeError(
-                "Bug encountered! Could not find a DC simulation workspace "
-                f"associated with Isolator of name {self.name}"
-            )
-        if len(DC_workspaces) > 1:
-            raise RuntimeError(
-                "Bug encountered! Found multiple workspaces in the DC simulation "
-                f"associated with Isolator of name {self.name}"
-            )
-
-        DC_ws = DC_workspaces[0]
-
-        couplings = ["12", "21"]
-        for c in couplings:
-            setattr(ws, f"K{c}", getattr(DC_ws, f"K{c}"))

src/finesse/components/lens.py

@@ -6,13 +6,10 @@ import logging
 
 import numpy as np
 
-from finesse.knm import KnmMatrix
-
 from finesse.components.modal.lens import (
-    LensWorkspace,
-    LensConnections,
+    LensWorkspace
 )
-from finesse.components.general import Connector, InteractionType, _conn_abcd_return
+from finesse.components.general import Connector, InteractionType
 from finesse.components.node import NodeDirection, NodeType
 from finesse.parameter import model_parameter, Rebuild
 from finesse.utilities import refractive_index
@@ -48,7 +45,7 @@ class Lens(Connector):
     """
 
     def __init__(self, name, f=np.inf):
-        super().__init__(name, LensConnections)
+        super().__init__(name)
 
         self.f = f
 
@@ -90,9 +87,9 @@ class Lens(Connector):
         M_sym = np.array([[1.0, 0.0], [-1.0 / self.f.ref, 1.0]])
         M_num = np.array(M_sym, dtype=np.float64)
         key = (self.p1.i, self.p2.o, direction)
-        self._abcd_matrices[key] = ( M_sym, M_num )
+        self._abcd_matrices[key] = (M_sym, M_num)
         key = (self.p2.i, self.p1.o, direction)
-        self._abcd_matrices[key] = ( M_sym, M_num )
+        self._abcd_matrices[key] = (M_sym, M_num)
 
     @property
     def abcdx(self):
@@ -197,13 +194,6 @@ class Lens(Connector):
         ws.nr1 = refractive_index(self.p1) or 1
         ws.nr2 = refractive_index(self.p2) or 1
 
-        # Only allocate new memory for scattering matrices for the carrier
-        # simulation. Audio simulation will use the carrier Knm matrices.
-        if not sim.is_audio:
-            self.__allocate_knm_matrices(ws)
-        else:
-            self.__set_audio_ws_knm_matrix_refs(ws)
-
         ws.set_fill_fn(self._fill_matrix)
 
         if sim.is_modal:
@@ -213,28 +203,3 @@ class Lens(Connector):
             _, ws.abcd_y = self._abcd_matrices[key]
 
         return ws
-
-    def __allocate_knm_matrices(self, ws: LensWorkspace):
-        Is = np.eye(ws.sim.nhoms, dtype=np.complex128)
-        ws.K12 = KnmMatrix(Is.copy(), self, kdir="12")
-        ws.K21 = KnmMatrix(Is.copy(), self, kdir="21")
-
-    def __set_audio_ws_knm_matrix_refs(self, ws):
-        DC = ws.sim.DC
-        DC_workspaces = list(filter(lambda x: x.owner is self, DC.workspaces))
-        if not DC_workspaces:
-            raise RuntimeError(
-                "Bug encountered! Could not find a DC simulation workspace "
-                f"associated with Lens of name {self.name}"
-            )
-        if len(DC_workspaces) > 1:
-            raise RuntimeError(
-                "Bug encountered! Found multiple workspaces in the DC simulation "
-                f"associated with Lens of name {self.name}"
-            )
-
-        DC_ws = DC_workspaces[0]
-
-        couplings = ["12", "21"]
-        for c in couplings:
-            setattr(ws, f"K{c}", getattr(DC_ws, f"K{c}"))

src/finesse/components/mirror.py

@@ -7,15 +7,13 @@ import logging
 import numpy as np
 
 from finesse.parameter import model_parameter, Rebuild
-from finesse.knm import KnmMatrix
 from finesse.utilities import refractive_index
 
 from finesse.components.modal.mirror import (
     mirror_fill,
     MirrorWorkspace,
-    MirrorConnections,
 )
-from finesse.components.general import InteractionType, _conn_abcd_return
+from finesse.components.general import InteractionType
 from finesse.components.surface import Surface
 from finesse.components.node import NodeDirection, NodeType
 
@@ -99,7 +97,7 @@ class Mirror(Surface):
         signal_gain=1e-9,
     ):
         super().__init__(
-            name, R, T, L, phi, Rc, xbeta, ybeta, mass, signal_gain, MirrorConnections
+            name, R, T, L, phi, Rc, xbeta, ybeta, mass, signal_gain
         )
 
         self._add_port("p1", NodeType.OPTICAL)
@@ -336,7 +334,6 @@ class Mirror(Surface):
         """
         return super().ABCD(from_node, to_node, direction, symbolic, copy, retboth)
 
-
     def _get_workspace(self, sim):
         _, is_changing = self._eval_parameters()
 
@@ -353,12 +350,6 @@ class Mirror(Surface):
         ws.nr2 = refractive_index(self.p2) or 1
 
         ws.set_fill_fn(mirror_fill)
-        # Only allocate new memory for scattering matrices for the carrier
-        # simulation. Audio simulation will use the carrier Knm matrices.
-        if not sim.is_audio:
-            self.__allocate_knm_matrices(ws)
-        else:
-            self.__set_audio_ws_knm_matrix_refs(ws)
 
         # Initialise the ABCD matrix memory-views
         if sim.is_modal:
@@ -386,33 +377,3 @@ class Mirror(Surface):
             for freq in sim.frequencies:
                 sim.Mq[sim.field(self.p1.o, freq.index, 0)] = self.L / 2
                 sim.Mq[sim.field(self.p2.o, freq.index, 0)] = self.L / 2
-
-    def __allocate_knm_matrices(self, ws: MirrorWorkspace):
-        Is = np.eye(ws.sim.nhoms, dtype=np.complex128)
-        losses = np.ones(ws.sim.nhoms)
-
-        couplings = ["11", "12", "21", "22"]
-        for c in couplings:
-            setattr(ws, f"K{c}", KnmMatrix(Is.copy(), self, kdir=c))
-            setattr(ws, f"K{c}_loss", losses.copy())
-
-    def __set_audio_ws_knm_matrix_refs(self, ws):
-        DC = ws.sim.DC
-        DC_workspaces = list(filter(lambda x: x.owner is self, DC.workspaces))
-        if not DC_workspaces:
-            raise RuntimeError(
-                "Bug encountered! Could not find a DC simulation workspace "
-                f"associated with Mirror of name {self.name}"
-            )
-        if len(DC_workspaces) > 1:
-            raise RuntimeError(
-                "Bug encountered! Found multiple workspaces in the DC simulation "
-                f"associated with Mirror of name {self.name}"
-            )
-
-        DC_ws = DC_workspaces[0]
-
-        couplings = ["11", "12", "21", "22"]
-        for c in couplings:
-            setattr(ws, f"K{c}", getattr(DC_ws, f"K{c}"))
-            setattr(ws, f"K{c}_loss", getattr(DC_ws, f"K{c}_loss"))

src/finesse/components/modal/beamsplitter.pxd

@@ -5,6 +5,7 @@ from finesse.simulations.base cimport frequency_info_t, ModelData, NodeBeamParam
 from finesse.simulations cimport BaseSimulation
 from finesse.element cimport BaseCValues
 from finesse.components.workspace cimport ConnectorWorkspace, FillFuncWrapper
+from finesse.components.modal.workspace cimport KnmConnectorWorkspace
 
 from finesse.cyexpr cimport (
     cy_expr,
@@ -80,7 +81,7 @@ cdef class BeamsplitterValues(BaseCValues):
         double alpha
 
 
-cdef class BeamsplitterWorkspace(ConnectorWorkspace):