Re-implemented missing knm_loss calculations

src/finesse/components/modal/isolator.pxd

@@ -35,6 +35,10 @@ cdef class IsolatorWorkspace(KnmConnectorWorkspace):
         KnmMatrix K12
         KnmMatrix K21
 
+        # Arrays of scattering losses for each mode coupling
+        double[::1] K12_loss
+        double[::1] K21_loss
+
         # Refractive indices of adjacent spaces
         double nr1
         double nr2

src/finesse/components/modal/lens.pxd

@@ -45,6 +45,10 @@ cdef class LensWorkspace(KnmConnectorWorkspace):
         KnmMatrix K12
         KnmMatrix K21
 
+        # Arrays of scattering losses for each mode coupling
+        double[::1] K12_loss
+        double[::1] K21_loss
+
         # Refractive indices of adjacent spaces
         double nr1
         double nr2

src/finesse/components/modal/modulator.pxd

@@ -75,6 +75,11 @@ cdef class ModulatorWorkspace(KnmConnectorWorkspace):
         # Complete scattering matrices for each propagation direction
         KnmMatrix K12
         KnmMatrix K21
+
+        # Arrays of scattering losses for each mode coupling
+        double[::1] K12_loss
+        double[::1] K21_loss
+
         # Refractive indices of adjacent spaces
         double nr1
         double nr2

src/finesse/components/modal/space.pyx

@@ -448,37 +448,29 @@ cdef int c_set_gouy_phase(ConnectorWorkspace cws) except -1:
     # components, so if a mismatch is present then report a bug
     # with an informative error message.
     if not ceq(q_p2i.qx, qx_p1o_propagated):
-        node1_name = list(ws.sim.carrier.nodes.keys())[ws.P1i_id]
-        node2_name = list(ws.sim.carrier.nodes.keys())[ws.P2o_id]
-        z_diff = fabs(creal(q_p2i.qx) - creal(qx_p1o_propagated))
-        zr_diff = fabs(cimag(q_p2i.qx) - cimag(qx_p1o_propagated))
-        mag_diff = cabs(q_p2i.qx - qx_p1o_propagated)
         raise RuntimeError(
-            f"Bug encountered! Mismatch in space {ws.owner.name} "
-            f"from {node1_name} -> {node2_name}:"
-            f"\n    STORED qx @ {node2_name} = {q_p2i.qx}"
-            f"\n    EXPECTED qx = {qx_p1o_propagated}"
-            f"\n        given qx @ {node1_name} = {q_p1o.qx},"
-            f"\n        with ABCD = {ws.abcd.base.tolist()}"
-            f"\n    DIFFERENCES: |z2 - z1| = {z_diff}, |zr2 - zr1| = {zr_diff};"
-            f"\n                 |q1 - q2| = {mag_diff}"
+            _space_mismatch_bug_report(
+                ws,
+                ws.P1i_id,
+                ws.P2o_id,
+                q_p1o.qx,
+                q_p2i.qx,
+                qx_p1o_propagated,
+                "x"
+            )
         )
 
     if not ceq(q_p2i.qy, qy_p1o_propagated):
-        node1_name = list(ws.sim.carrier.nodes.keys())[ws.P1i_id]
-        node2_name = list(ws.sim.carrier.nodes.keys())[ws.P2o_id]
-        z_diff = fabs(creal(q_p2i.qy) - creal(qy_p1o_propagated))
-        zr_diff = fabs(cimag(q_p2i.qy) - cimag(qy_p1o_propagated))
-        mag_diff = cabs(q_p2i.qy - qy_p1o_propagated)
         raise RuntimeError(
-            f"Bug encountered! Mismatch in space: {ws.owner.name} "
-            f"from {node1_name} -> {node2_name}:"
-            f"\n\tSTORED qy @ {node2_name} = {q_p2i.qy}"
-            f"\n\tEXPECTED qy = {qy_p1o_propagated}"
-            f"\n\t\tgiven qy @ {node1_name} = {q_p1o.qy},"
-            f"\n\t\twith ABCD = {ws.abcd.base.tolist()}"
-            f"\n    DIFFERENCES: |z2 - z1| = {z_diff}, |zr2 - zr1| = {zr_diff};"
-            f"\n                 |q1 - q2| = {mag_diff}"
+            _space_mismatch_bug_report(
+                ws,
+                ws.P1i_id,
+                ws.P2o_id,
+                q_p1o.qy,
+                q_p2i.qy,
+                qy_p1o_propagated,
+                "y"
+            )
         )
 
     ws.owner.gouy_x = degrees(
@@ -489,3 +481,31 @@ cdef int c_set_gouy_phase(ConnectorWorkspace cws) except -1:
     )
 
     return 0
+
+
+cdef _space_mismatch_bug_report(
+    SpaceWorkspace ws,
+    Py_ssize_t node1_id,
+    Py_ssize_t node2_id,
+    complex_t q1,
+    complex_t q2,
+    complex_t q2e, # expected q
+    direction, # "x" or "y"
+):
+    node1_name = list(ws.sim.carrier.nodes.keys())[node1_id]
+    node2_name = list(ws.sim.carrier.nodes.keys())[node2_id]
+
+    z_diff = fabs(creal(q2e) - creal(q2))
+    zr_diff = fabs(cimag(q2e) - cimag(q2))
+    mag_diff = cabs(q2e - q2)
+
+    return (
+        f"Bug encountered! Mismatch in space {ws.owner.name} "
+        f"from {node1_name} -> {node2_name}:"
+        f"\n    STORED q{direction} @ {node2_name} = {q2}"
+        f"\n    EXPECTED q{direction} = {q2e}"
+        f"\n        given q{direction} @ {node1_name} = {q1},"
+        f"\n        with ABCD = {ws.abcd.base.tolist()}"
+        f"\n    DIFFERENCES: |z2 - z1| = {z_diff}, |zr2 - zr1| = {zr_diff};"
+        f"\n                 |q1 - q2| = {mag_diff}"
+    )

src/finesse/components/modal/workspace.pxd

@@ -36,6 +36,8 @@ cdef struct KnmInfo:
     # Pointer to bayerhelms scatter matrix function to use
     # -> see end of set_knm_info definition for details
     void (*bhelms_mtx_func)(KnmWorkspace*, KnmWorkspace*, const int[:, ::1], complex_t*)
+    # Pointer to buffer of knm loss data
+    double* loss
 
 
 cdef class KnmConnectorWorkspace(ConnectorWorkspace):
@@ -52,6 +54,7 @@ cdef class KnmConnectorWorkspace(ConnectorWorkspace):
     cdef void flag_changing_knm_workspaces(KnmConnectorWorkspace self)
     cdef void update_changing_knm_workspaces(KnmConnectorWorkspace self) nogil
     cdef void compute_scattering_matrices(KnmConnectorWorkspace self)
+    cdef void compute_knm_losses(KnmConnectorWorkspace self) nogil
 
     cpdef set_knm_info(
         self,

src/finesse/components/modal/workspace.pyx

@@ -22,6 +22,7 @@ from finesse.knm.bayerhelms cimport (
     fast_compute_refl_knm_matrix_bh,
     fast_compute_trns_knm_matrix_bh,
 )
+from finesse.knm.matrix cimport knm_loss
 from finesse.cymath.gaussbeam cimport c_transform_q
 from finesse.components.workspace cimport ConnectorWorkspace
 from finesse.components.general import NodeType
@@ -72,6 +73,7 @@ cdef class KnmConnectorWorkspace(ConnectorWorkspace):
             self.onode_ids[2*i] = sim.carrier.node_id(p.i)
             self.onode_ids[2*i+1] = sim.carrier.node_id(p.o)
 
+        cdef np.ndarray[double, ndim=1, mode='c'] scatter_loss
         for i, (conn, (f, t)) in enumerate(self.o2o.items()):
             # TODO ddb should probably use some fixed index rather than a list of the order
             # they are defined in the element definition
@@ -91,6 +93,12 @@ cdef class KnmConnectorWorkspace(ConnectorWorkspace):
             # port coupling indices
             setattr(self, f"K{coupling}", knm)
 
+            # Make the buffer of knm loss data
+            scatter_loss = np.zeros(self.sim.model_data.num_HOMs)
+            setattr(self, f"K{coupling}_loss", scatter_loss)
+            # Keep ptr to the buffer for use in compute_knm_losses
+            self.oconn_info[i].loss = &scatter_loss[0]
+
     def __dealloc__(self):
         if self.onode_ids:
             free(self.onode_ids)
@@ -297,3 +305,16 @@ cdef class KnmConnectorWorkspace(ConnectorWorkspace):
                     self.sim.model_data.homs_view,
                     info.mtx.ptr
                 )
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    @cython.initializedcheck(False)
+    cdef void compute_knm_losses(KnmConnectorWorkspace self) nogil:
+        cdef:
+            KnmInfo *info
+            Py_ssize_t i
+
+        for i in range(self.N_opt_conns):
+            info = &self.oconn_info[i]
+            if knm_ws_is_changing(info.K_ws_x) or knm_ws_is_changing(info.K_ws_y):
+                knm_loss(info.mtx.ptr, info.loss, self.sim.model_data.num_HOMs)

src/finesse/knm/bayerhelms.pyx

@@ -736,21 +736,6 @@ cdef complex_t rev_gouy(
 
     return crotate(k, phase2 - phase1)
 
-# TODO (sjr) Move this elsewhere
-cdef void knm_loss(
-    const complex_t[:, ::1] knm_mat,
-    double[::1] out
-) nogil:
-    """Computes total scattering losses from each mode (n, m)
-    to all mode couplings."""
-    cdef:
-        Py_ssize_t i, j
-        Py_ssize_t N = knm_mat.shape[0]
-
-    for i in range(N):
-        for j in range(N):
-            out[i] -= cnorm(knm_mat[j][i])
-
 
 ### The functions for computing the coefficients themselves ###
 

src/finesse/knm/matrix.pxd

@@ -18,3 +18,8 @@ cdef class KnmMatrix:
 
 cpdef make_unscaled_X_scatter_knm_matrix(int[:,::1] modes)
 cpdef make_unscaled_Y_scatter_knm_matrix(int[:,::1] modes)
+
+
+# Compute loss from scattering for each coupling,
+# required for  quantum noise calculations
+cdef void knm_loss(const complex_t* knm_mat, double* out, Py_ssize_t N) nogil

src/finesse/knm/matrix.pyx

@@ -3,6 +3,7 @@
 cimport numpy as np
 import numpy as np
 
+from finesse.cymath.complex cimport cnorm
 from finesse.cymath.homs cimport field_index
 
 import copy
@@ -417,3 +418,17 @@ cpdef make_unscaled_Y_scatter_knm_matrix(int[:, ::1] modes):
                     k.data[i, j] = np.sqrt(m+1)/2
 
     return k
+
+
+cdef void knm_loss(const complex_t* knm_mat, double* out, Py_ssize_t N) nogil:
+    """Compute per mode losses due to scattering.
+
+    Each entry in out is then 1 - \sum_{n'm'} |k_{n'm'nm}|^2 for the mode nm.
+    """
+    cdef:
+        Py_ssize_t i, j
+
+    for i in range(N):
+        out[i] = 1.0
+        for j in range(N):
+            out[i] -= cnorm(knm_mat[j*N + i])

src/finesse/simulations/basematrix.pyx

@@ -971,6 +971,9 @@ cdef class CarrierSignalMatrixSimulation:
         for ws in self.to_scatter_matrix_compute:
             ws.update_changing_knm_workspaces()
             ws.compute_scattering_matrices()
+            # TODO (sjr) Probably want a flag to check if quantum noise calcs
+            #            being performed and to only do this call if so
+            ws.compute_knm_losses()
 
     cdef int set_gouy_phases(self) except -1:
         cdef ConnectorWorkspace ws