Add `nsr` flag to the various quantum detectors.

This allows `qnoisedS` and `qshotS` to parse and work properly, however
only 0 or 1 RF demodulations are currently supported (as the necessary
powerdetector functions are missing)

src/finesse/detectors/compute/power.pxd

+from finesse.cymath cimport complex_t
+
+from finesse.detectors.workspace cimport (
+        DetectorWorkspace,
+        MaskedDetectorWorkspace,
+)
+from finesse.simulations cimport BaseSimulation
+from finesse.element cimport BaseCValues
+
+cdef class PD1Workspace(MaskedDetectorWorkspace):
+    cdef public:
+        bint output_real
+        BaseSimulation DC
+        BaseSimulation AC
+        int dc_node_id
+        int ac_node_id
+        int[:,::1] homs
+        bint is_f_changing
+        bint is_phase_changing
+        bint is_audio_mixing
+        complex_t Aij
+        PD1Values v
+    cdef:
+        Py_ssize_t num_pre_set_beats
+        Py_ssize_t* pre_set_beats[2] # 2 x num_pre_set_beats
+
+    cpdef update_beats(self)
+
+cdef class PD1Values(BaseCValues):
+    cdef public:
+        double f
+        double phase
+
+cdef object c_pd1_AC_output(DetectorWorkspace dws)
+
+cdef class PD2Workspace(MaskedDetectorWorkspace):
+    cdef public:
+        bint output_real
+        BaseSimulation DC
+        BaseSimulation AC
+        int dc_node_id
+        int ac_node_id
+        int[:,::1] homs
+        bint is_f1_changing
+        bint is_f2_changing
+        bint is_phase1_changing
+        bint is_phase2_changing
+        bint is_audio_mixing
+        complex_t z1
+        complex_t z2
+        PD2Values v
+
+cdef class PD2Values(BaseCValues):
+    cdef public:
+        double f1
+        double phase1
+        double f2
+        double phase2
+
+cdef object c_pd2_AC_output(DetectorWorkspace dws)

src/finesse/detectors/compute/power.pyx

@@ -91,10 +91,6 @@ cdef c_pd0_DC_output_masked(DetectorWorkspace dws):
 ### PD1 workspace & output funcs ###
 
 cdef class PD1Values(BaseCValues):
-    cdef public:
-        double f
-        double phase
-
     def __init__(self):
         cdef ptr_tuple_2 ptr = (&self.f, &self.phase)
         cdef tuple params = ("f", "phase")
@@ -102,22 +98,6 @@ cdef class PD1Values(BaseCValues):
 
 
 cdef class PD1Workspace(MaskedDetectorWorkspace):
-    cdef public:
-        bint output_real
-        BaseSimulation DC
-        BaseSimulation AC
-        int dc_node_id
-        int ac_node_id
-        int[:,::1] homs
-        bint is_f_changing
-        bint is_phase_changing
-        bint is_audio_mixing
-        complex_t Aij
-        PD1Values v
-    cdef:
-        Py_ssize_t num_pre_set_beats
-        Py_ssize_t* pre_set_beats[2] # 2 x num_pre_set_beats
-
     def __cinit__(self, *args):
         self.num_pre_set_beats = 0
         self.pre_set_beats[0] = NULL
@@ -264,12 +244,6 @@ cdef object c_pd1_DC_output(DetectorWorkspace dws):
 ### PD2 workspace & output funcs ###
 
 cdef class PD2Values(BaseCValues):
-    cdef public:
-        double f1
-        double phase1
-        double f2
-        double phase2
-
     def __init__(self):
         cdef ptr_tuple_4 ptr = (&self.f1, &self.phase1, &self.f2, &self.phase2)
         cdef tuple params = ("f1", "phase1", "f2", "phase2")
@@ -277,22 +251,6 @@ cdef class PD2Values(BaseCValues):
 
 
 cdef class PD2Workspace(MaskedDetectorWorkspace):
-    cdef public:
-        bint output_real
-        BaseSimulation DC
-        BaseSimulation AC
-        int dc_node_id
-        int ac_node_id
-        int[:,::1] homs
-        bint is_f1_changing
-        bint is_f2_changing
-        bint is_phase1_changing
-        bint is_phase2_changing
-        bint is_audio_mixing
-        complex_t z1
-        complex_t z2
-        PD2Values v
-
     def __init__(self, owner, sims):
         super().__init__(owner, sims, PD2Values())
         self.v = <PD2Values>self.values

src/finesse/detectors/compute/quantum.pyx

@@ -2,6 +2,7 @@
 
 cimport numpy as np
 import numpy as np
+import cython
 
 from libc.stdlib cimport calloc, free
 from finesse.cymath cimport complex_t
@@ -11,6 +12,12 @@ from finesse.detectors.workspace cimport (
     DetectorWorkspace,
     OutputFuncWrapper,
 )
+from finesse.detectors.compute.power cimport (
+    PD1Workspace,
+    PD2Workspace,
+    c_pd1_AC_output,
+    c_pd2_AC_output,
+)
 from finesse.simulations cimport BaseSimulation
 from finesse.simulations.base cimport frequency_info_t
 from finesse.element cimport BaseCValues
@@ -53,11 +60,13 @@ cdef class QND0Workspace(DetectorWorkspace):
         complex_t[::1] cov
         BaseSimulation DC
         BaseSimulation AC
+        PD1Workspace pd_ws
         int dc_node_id
         int ac_node_id
         int neq
+        bint nsr
 
-    def __init__(self, owner, sims):
+    def __init__(self, owner, sims, nsr):
         super().__init__(owner, sims)
         for sim in sims:
             if sim.is_audio:
@@ -71,6 +80,14 @@ cdef class QND0Workspace(DetectorWorkspace):
         self.s = np.zeros(self.neq, dtype=np.complex128)
         self.set_output_fn(qnd0_output)
 
+        self.nsr = nsr
+        if self.nsr:
+            self.pd_ws = PD1Workspace(owner, sims)
+            self.pd_ws.DC = self.DC
+            self.pd_ws.AC = self.AC
+            self.pd_ws.dc_node_id = self.dc_node_id
+            self.pd_ws.ac_node_id = self.ac_node_id
+
     cpdef fill_selection_vector(self):
         cdef:
             Py_ssize_t i, j
@@ -78,6 +95,7 @@ cdef class QND0Workspace(DetectorWorkspace):
             frequency_info_t fc, fs
             complex_t tmp
 
+        self.update_parameter_values()
         for i in range(self.AC.num_frequencies):
             fs = self.AC.frequency_info[i]
             fc = self.DC.frequency_info[fs.audio_carrier_index]
@@ -97,13 +115,17 @@ cdef class QND0Workspace(DetectorWorkspace):
         self.cov = v
 
 qnd0_output = OutputFuncWrapper.make_from_ptr(c_qnd0_output)
-cdef c_qnd0_output(DetectorWorkspace self):
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.initializedcheck(False)
+cdef object c_qnd0_output(DetectorWorkspace self):
     cdef:
         QND0Workspace ws = <QND0Workspace> self
         double rtn = 0
         double f0 = ws.AC.model_data.f0
         double unit_vacuum = ws.AC.model_data.UNIT_VACUUM
         double nf_factor = 0.25**1  # factor 0.25 per demodulation
+        complex_t pdo
 
     for i in range(ws.neq):
         rtn += creal(ws.cov[i] * conj(ws.s[i]))
@@ -112,6 +134,10 @@ cdef c_qnd0_output(DetectorWorkspace self):
     # frequency as this is what a network analyser would do.
     rtn *= 2
 
+    if ws.nsr:
+        pdo = c_pd1_AC_output(ws.pd_ws)
+        rtn /= pdo.real * pdo.real + pdo.imag * pdo.imag
+
     return sqrt(
         unit_vacuum
         * rtn
@@ -126,6 +152,7 @@ cdef class QNDNWorkspace(DetectorWorkspace):
         complex_t[::1] cov
         BaseSimulation DC
         BaseSimulation AC
+        DetectorWorkspace pd_ws
         int num_demod
         int Ndm
         int Nf
@@ -137,12 +164,13 @@ cdef class QNDNWorkspace(DetectorWorkspace):
         double[::1] demod_f
         double[::1] demod_phi
         QNDNValues v
+        bint nsr
 
         double[::1] freqs
         double[::1] phases
         int neq
 
-    def __init__(self, owner, sims, num_demod):
+    def __init__(self, owner, sims, num_demod, nsr):
         if num_demod == 1:
             values = QND1Values()
         else:
@@ -175,6 +203,18 @@ cdef class QNDNWorkspace(DetectorWorkspace):
         self.freqs = np.empty(self.num_demod, dtype=np.float64)
         self.phases = np.empty(self.num_demod, dtype=np.float64)
 
+        self.nsr = nsr
+        if self.nsr:
+            if self.num_demod == 2:
+                self.pd_ws = PD2Workspace(owner, sims)
+                (<PD2Workspace>self.pd_ws).DC = self.DC
+                (<PD2Workspace>self.pd_ws).AC = self.AC
+                (<PD2Workspace>self.pd_ws).dc_node_id = self.dc_node_id
+                (<PD2Workspace>self.pd_ws).ac_node_id = self.ac_node_id
+                (<PD2Workspace>self.pd_ws).setup_mask(self.DC)
+            else:
+                raise ValueError(f"Can't calculate NSR for {self.num_demod-1} RF demodulations")
+
     cdef fill_carrier_qnoise_contributions(self):
         cdef:
             Py_ssize_t i, j, k, nf
@@ -250,6 +290,7 @@ cdef class QNDNWorkspace(DetectorWorkspace):
             frequency_info_t fc, fs
             complex_t tmp
 
+        self.update_parameter_values()
         self.fill_carrier_qnoise_contributions()
 
         self.s[:] = 0
@@ -282,7 +323,10 @@ cdef class QNDNWorkspace(DetectorWorkspace):
 
 
 qndN_output = OutputFuncWrapper.make_from_ptr(c_qndN_output)
-cdef c_qndN_output(DetectorWorkspace self):
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.initializedcheck(False)
+cdef object c_qndN_output(DetectorWorkspace self):
     cdef:
         QNDNWorkspace ws = <QNDNWorkspace> self
         Py_ssize_t i, j, k
@@ -292,6 +336,7 @@ cdef c_qndN_output(DetectorWorkspace self):
         double f0 = ws.AC.model_data.f0
         double unit_vacuum = ws.AC.model_data.UNIT_VACUUM
         double nf_factor = 0.25**ws.num_demod  # factor 0.25 per demodulation
+        complex_t pdo
         complex_t car_sum
 
     for i in range(ws.neq):
@@ -312,6 +357,17 @@ cdef c_qndN_output(DetectorWorkspace self):
                         )
             rtn += cabs(car_sum)**2 * (1 + f / f0)
 
+    if ws.nsr:
+        if ws.num_demod == 2:
+            (<PD2Workspace>ws.pd_ws).v.f1 = ws.v.f1
+            (<PD2Workspace>ws.pd_ws).v.phase1 = ws.v.phase1
+            (<PD2Workspace>ws.pd_ws).v.f2 = 0
+            (<PD2Workspace>ws.pd_ws).v.phase2 = 0
+        else:
+            raise ValueError(f"Can't calculate NSR for {ws.num_demod-1} RF demodulations")
+        pdo = c_pd2_AC_output(ws.pd_ws)
+        rtn /= pdo.real * pdo.real + pdo.imag * pdo.imag
+
     # Compensate for demod 0.5 factor when demodulating a signal
     # frequency as this is what a network analyser would do.
     rtn *= 2
@@ -328,12 +384,15 @@ cdef class QShot0Workspace(DetectorWorkspace):
     cdef:
         BaseSimulation DC
         BaseSimulation AC
+        PD1Workspace pd_ws
         int Nf
         int dc_node_id
+        int ac_node_id
         long[:, ::1] demod_vac_contri
         double[::1] demod_f
+        bint nsr
 
-    def __init__(self, owner, sims):
+    def __init__(self, owner, sims, nsr):
         super().__init__(owner, sims)
         for sim in sims:
             if sim.is_audio:
@@ -341,6 +400,7 @@ cdef class QShot0Workspace(DetectorWorkspace):
                 self.AC = sim
 
         self.dc_node_id = self.DC.node_id(owner._node)
+        self.ac_node_id = self.AC.node_id(owner._node)
         self.set_output_fn(qshot0_output)
 
         # As we're only doing 1 demodulation, we have a maximum of 2 signal frequencies for each
@@ -350,6 +410,14 @@ cdef class QShot0Workspace(DetectorWorkspace):
         self.demod_vac_contri = np.empty((self.Nf, 2), dtype=np.int)
         self.demod_f = np.full(self.Nf, np.nan, dtype=np.float64)
 
+        self.nsr = nsr
+        if self.nsr:
+            self.pd_ws = PD1Workspace(owner, sims)
+            self.pd_ws.DC = self.DC
+            self.pd_ws.AC = self.AC
+            self.pd_ws.dc_node_id = self.dc_node_id
+            self.pd_ws.ac_node_id = self.ac_node_id
+
     cdef fill_carrier_qnoise_contributions(self):
         cdef:
             Py_ssize_t i
@@ -386,7 +454,10 @@ cdef class QShot0Workspace(DetectorWorkspace):
 
 
 qshot0_output = OutputFuncWrapper.make_from_ptr(c_qshot0_output)
-cdef c_qshot0_output(DetectorWorkspace self):
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.initializedcheck(False)
+cdef object c_qshot0_output(DetectorWorkspace self):
     cdef:
         QShot0Workspace ws = <QShot0Workspace> self
         Py_ssize_t i, j
@@ -397,6 +468,7 @@ cdef c_qshot0_output(DetectorWorkspace self):
         double unit_vacuum = ws.AC.model_data.UNIT_VACUUM
         double nf_factor = 0.25**1  # factor 0.25 per demodulation
         complex_t c1, c2
+        complex_t pdo
 
     ws.fill_carrier_qnoise_contributions()
 
@@ -419,6 +491,10 @@ cdef c_qshot0_output(DetectorWorkspace self):
     # frequency as this is what a network analyser would do.
     rtn *= 2
 
+    if ws.nsr:
+        pdo = c_pd1_AC_output(ws.pd_ws)
+        rtn /= pdo.real * pdo.real + pdo.imag * pdo.imag
+
     return sqrt(
         unit_vacuum
         * rtn
@@ -432,21 +508,24 @@ cdef class QShotNWorkspace(DetectorWorkspace):
     cdef:
         BaseSimulation DC
         BaseSimulation AC
+        DetectorWorkspace pd_ws
         int num_demod
         int Ndm
         int Nf
         int dc_node_id
+        int ac_node_id
         long[:, ::1] demod_vac_contri
         long[:, ::1] demod_vac_contri_phi
         double[::1] demod_f
         double[::1] demod_phi
         QNDNValues v
+        bint nsr
 
         double[::1] freqs
         double[::1] phases
         double[::1] tmp_fs
 
-    def __init__(self, owner, sims, num_demod):
+    def __init__(self, owner, sims, num_demod, nsr):
         if num_demod == 1:
             values = QND1Values()
         else:
@@ -459,6 +538,7 @@ cdef class QShotNWorkspace(DetectorWorkspace):
                 self.AC = sim
 
         self.dc_node_id = self.DC.node_id(owner._node)
+        self.ac_node_id = self.AC.node_id(owner._node)
         self.set_output_fn(qshotN_output)
 
         # N demodulations, so we have a maximum of 2**N signal frequencies for each carrier
@@ -476,6 +556,18 @@ cdef class QShotNWorkspace(DetectorWorkspace):
         self.phases = np.empty(self.num_demod, dtype=np.float64)
         self.tmp_fs = np.empty(self.Ndm, dtype=np.float64)
 
+        self.nsr = nsr
+        if self.nsr:
+            if self.num_demod == 2:
+                self.pd_ws = PD2Workspace(owner, sims)
+                (<PD2Workspace>self.pd_ws).DC = self.DC
+                (<PD2Workspace>self.pd_ws).AC = self.AC
+                (<PD2Workspace>self.pd_ws).dc_node_id = self.dc_node_id
+                (<PD2Workspace>self.pd_ws).ac_node_id = self.ac_node_id
+                (<PD2Workspace>self.pd_ws).setup_mask(self.DC)
+            else:
+                raise ValueError(f"Can't calculate NSR for {self.num_demod-1} RF demodulations")
+
     cdef fill_carrier_qnoise_contributions(self):
         cdef:
             Py_ssize_t i, j, k
@@ -534,7 +626,10 @@ cdef class QShotNWorkspace(DetectorWorkspace):
 
 
 qshotN_output = OutputFuncWrapper.make_from_ptr(c_qshotN_output)
-cdef c_qshotN_output(DetectorWorkspace self):
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.initializedcheck(False)
+cdef object c_qshotN_output(DetectorWorkspace self):
     cdef:
         QShotNWorkspace ws = <QShotNWorkspace> self
         Py_ssize_t i, j, k
@@ -545,6 +640,7 @@ cdef c_qshotN_output(DetectorWorkspace self):
         double unit_vacuum = ws.AC.model_data.UNIT_VACUUM
         double nf_factor = 0.25**ws.num_demod  # factor 0.25 per demodulation
         complex_t car_sum
+        complex_t pdo
 
     ws.fill_carrier_qnoise_contributions()
 
@@ -569,6 +665,17 @@ cdef c_qshotN_output(DetectorWorkspace self):
     # frequency as this is what a network analyser would do.
     rtn *= 2
 
+    if ws.nsr:
+        if ws.num_demod == 2:
+            (<PD2Workspace>ws.pd_ws).v.f1 = ws.v.f1
+            (<PD2Workspace>ws.pd_ws).v.phase1 = ws.v.phase1
+            (<PD2Workspace>ws.pd_ws).v.f2 = 0
+            (<PD2Workspace>ws.pd_ws).v.phase2 = 0
+        else:
+            raise ValueError(f"Can't calculate NSR for {ws.num_demod-1} RF demodulations")
+        pdo = c_pd2_AC_output(ws.pd_ws)
+        rtn /= pdo.real * pdo.real + pdo.imag * pdo.imag
+
     return sqrt(
         unit_vacuum
         * rtn

src/finesse/detectors/quantum_noise_detector.py

@@ -29,19 +29,24 @@ class QuantumNoiseDetector(Detector, NoiseDetector):
 
     node : :class:`.Node`
         Node to read output from.
+
+    nsr : bool, optional
+        If true, calculate the noise-to-signal ratio rather than the absolute
+        noise value.
     """
 
-    def __init__(self, name, node):
+    def __init__(self, name, node, nsr=False):
         Detector.__init__(
             self, name, node, dtype=np.float64, unit="1/sqrt(Hz)", label="ASD"
         )
         NoiseDetector.__init__(self, NoiseType.QUANTUM)
 
+        self.__nsr = nsr
+
         self._request_selection_vector(name)
 
     def _get_workspace(self, sims):
-        ws = QND0Workspace(self, sims)
-        ws.update_parameter_values()
+        ws = QND0Workspace(self, sims, self.__nsr)
         return ws
 
 
@@ -61,9 +66,19 @@ class QuantumNoiseDetectorDemod1(Detector, NoiseDetector):
 
     node : :class:`.Node`
         Node to read output from.
+
+    f : float
+        Demodulation frequency in Hz
+
+    phase : float
+        Demodulation phase in degrees
+
+    nsr : bool, optional
+        If true, calculate the noise-to-signal ratio rather than the absolute
+        noise value.
     """
 
-    def __init__(self, name, node, f, phase):
+    def __init__(self, name, node, f, phase, nsr=False):
         Detector.__init__(
             self, name, node, dtype=np.float64, unit="1/sqrt(Hz)", label="ASD"
         )
@@ -71,12 +86,12 @@ class QuantumNoiseDetectorDemod1(Detector, NoiseDetector):
 
         self.f = f
         self.phase = phase
+        self.__nsr = nsr
 
         self._request_selection_vector(name)
 
     def _get_workspace(self, sims):
-        ws = QNDNWorkspace(self, sims, 1)
-        ws.update_parameter_values()
+        ws = QNDNWorkspace(self, sims, 1, self.__nsr)
         return ws
 
 
@@ -98,9 +113,25 @@ class QuantumNoiseDetectorDemod2(Detector, NoiseDetector):
 
     node : :class:`.Node`
         Node to read output from.
+
+    f1 : float
+        First demodulation frequency in Hz
+
+    phase1 : float
+        First demodulation phase in degrees
+
+    f2 : float
+        Second demodulation frequency in Hz
+
+    phase2 : float
+        Second demodulation phase in degrees
+
+    nsr : bool, optional
+        If true, calculate the noise-to-signal ratio rather than the absolute
+        noise value.
     """
 
-    def __init__(self, name, node, f1, phase1, f2, phase2):
+    def __init__(self, name, node, f1, phase1, f2, phase2, nsr=False):
         Detector.__init__(
             self, name, node, dtype=np.float64, unit="1/sqrt(Hz)", label="ASD"
         )
@@ -110,12 +141,12 @@ class QuantumNoiseDetectorDemod2(Detector, NoiseDetector):
         self.f2 = f2
         self.phase1 = phase1
         self.phase2 = phase2
+        self.__nsr = nsr
 
         self._request_selection_vector(name)
 
     def _get_workspace(self, sims):
-        ws = QNDNWorkspace(self, sims, 2)
-        ws.update_parameter_values()
+        ws = QNDNWorkspace(self, sims, 2, self.__nsr)
         return ws
 
 
@@ -135,16 +166,21 @@ class QuantumShotNoiseDetector(Detector):
 
     node : :class:`.Node`
         Node to read output from.
+
+    nsr : bool, optional
+        If true, calculate the noise-to-signal ratio rather than the absolute
+        noise value.
     """
 
-    def __init__(self, name, node):
+    def __init__(self, name, node, nsr=False):
+        self.__nsr = nsr
+
         Detector.__init__(
             self, name, node, dtype=np.float64, unit="1/sqrt(Hz)", label="ASD"
         )
 
     def _get_workspace(self, sims):
-        ws = QShot0Workspace(self, sims)
-        ws.update_parameter_values()
+        ws = QShot0Workspace(self, sims, self.__nsr)
         return ws
 
 
@@ -166,19 +202,29 @@ class QuantumShotNoiseDetectorDemod1(Detector):
 
     node : :class:`.Node`
         Node to read output from.
+
+    f : float
+        Demodulation frequency in Hz
+
+    phase : float
+        Demodulation phase in degrees
+
+    nsr : bool, optional
+        If true, calculate the noise-to-signal ratio rather than the absolute
+        noise value.
     """
 
-    def __init__(self, name, node, f, phase=None):
+    def __init__(self, name, node, f, phase, nsr=False):
         Detector.__init__(
             self, name, node, dtype=np.float64, unit="1/sqrt(Hz)", label="ASD"
         )
 
         self.f = f
         self.phase = phase
+        self.__nsr = nsr
 
     def _get_workspace(self, sims):
-        ws = QShotNWorkspace(self, sims, 1)
-        ws.update_parameter_values()