diff --git a/docs/source/usage/homs/index.rst b/docs/source/usage/homs/index.rst index 668c3a178783dddb0c02c8ec0a44c2c07af7305d..33c6103515ea254fa6a0b23ffff954f4ee61d711 100644 --- a/docs/source/usage/homs/index.rst +++ b/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 diff --git a/docs/source/usage/homs/scattering_matrices.rst b/docs/source/usage/homs/scattering_matrices.rst new file mode 100644 index 0000000000000000000000000000000000000000..aadece2352919ca857042d7c963ac73a1ca33a76 --- /dev/null +++ b/docs/source/usage/homs/scattering_matrices.rst @@ -0,0 +1,101 @@ +.. 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)]) diff --git a/docs/source/usage/homs/selecting_modes.rst b/docs/source/usage/homs/selecting_modes.rst index 63222930e8e791adc880dd7c8271f2bc590f76d4..20ae98545346f192333ac660a3795e1d7081474a 100644 --- a/docs/source/usage/homs/selecting_modes.rst +++ b/docs/source/usage/homs/selecting_modes.rst @@ -1,5 +1,7 @@ .. include:: /defs.hrst +.. _selecting_modes: + Selecting the modes to model ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/setup.py b/setup.py index efcc91b8bfbc0518bef991ebf99c8dcc83ab1168..26ec70474c3534f21bec932b2e4c84767f72101e 100644 --- a/setup.py +++ b/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", {}), diff --git a/src/finesse/__main__.py b/src/finesse/__main__.py index ac8abe4eaf76ee86e7df6394f6224bbadeff2db8..50414ecf0d64d4c4fbcf955f34e872733070b8c0 100644 --- a/src/finesse/__main__.py +++ b/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 diff --git a/src/finesse/components/beamsplitter.py b/src/finesse/components/beamsplitter.py index e0710ec6253ca08e9e900a6b262a35d0977b0644..de147553e729a878424e1d7f60dfc8e9aa8ba815 100644 --- a/src/finesse/components/beamsplitter.py +++ b/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")) diff --git a/src/finesse/components/general.py b/src/finesse/components/general.py index 1b87d67bf8cf674bf39e1b628875267a8bdcbc2f..bdcc2d4c96f1a58a219c4df08956411626018f42 100644 --- a/src/finesse/components/general.py +++ b/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 diff --git a/src/finesse/components/isolator.py b/src/finesse/components/isolator.py index 1050bd154d45188451b551c2c7b853a90b5cc19c..af664f5546bacd27a3bf616e2c9c26c4205d4b45 100644 --- a/src/finesse/components/isolator.py +++ b/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}")) diff --git a/src/finesse/components/lens.py b/src/finesse/components/lens.py index ccda45028e551e3bfce957b7f510339b2f6bf981..57cd2258ac479f85baf4bf25acffb320144e8b18 100644 --- a/src/finesse/components/lens.py +++ b/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}")) diff --git a/src/finesse/components/mirror.py b/src/finesse/components/mirror.py index 16a9c4532bb8c15a498585ad196d614c44ed2818..202aaa66958c1f3acc0cacc8ab1b312b47e531ff 100644 --- a/src/finesse/components/mirror.py +++ b/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")) diff --git a/src/finesse/components/modal/beamsplitter.pxd b/src/finesse/components/modal/beamsplitter.pxd index 3f330844b1f9056d455d3fe598b6f2935c8da77b..74e3259884a99cb97729b874045d6862cce6f8fa 100644 --- a/src/finesse/components/modal/beamsplitter.pxd +++ b/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): +cdef class BeamsplitterWorkspace(KnmConnectorWorkspace): cdef public: # Complete scattering matrices for each propagation direction KnmMatrix K12 @@ -185,12 +186,11 @@ cdef class BeamsplitterWorkspace(ConnectorWorkspace): cpdef compile_abcd_cy_exprs(self) cpdef update_parameter_values(self) + cdef void initialise_knm_workspaces(self) nogil + cdef void free_knm_workspaces(self) nogil + cdef void flag_changing_knm_workspaces(self) + cdef void update_changing_knm_workspaces(self) nogil + cdef void compute_scattering_matrices(self) cdef object c_beamsplitter_fill(ConnectorWorkspace cws) - -cdef void initialise_knm_workspaces(BeamsplitterWorkspace ws) nogil -cdef void free_knm_workspaces(BeamsplitterWorkspace ws) nogil -cdef void flag_changing_knm_workspaces(BeamsplitterWorkspace ws) -cdef void update_changing_knm_workspaces(BeamsplitterWorkspace ws) nogil -cdef void compute_scattering_matrices(BeamsplitterWorkspace ws) diff --git a/src/finesse/components/modal/beamsplitter.pyx b/src/finesse/components/modal/beamsplitter.pyx index a7faf7686cea925993e8e147c96a174cb1a04576..e4c1fff783cd1ae78b888fbf42efab6a03e4e424 100644 --- a/src/finesse/components/modal/beamsplitter.pyx +++ b/src/finesse/components/modal/beamsplitter.pyx @@ -80,8 +80,10 @@ cdef class BeamsplitterValues(BaseCValues): cdef tuple params = ("R","T","L","phi","Rcx","Rcy","xbeta","ybeta","mass","signal_gain","alpha") self.setup(params, sizeof(ptr), &ptr) -cdef class BeamsplitterWorkspace(ConnectorWorkspace): +cdef class BeamsplitterWorkspace(KnmConnectorWorkspace): def __init__(self, owner, BaseSimulation sim, refill): + self.str_couplings = ("12", "21", "34", "43", "13", "31", "24", "42") + super().__init__( owner, sim, @@ -123,9 +125,9 @@ cdef class BeamsplitterWorkspace(ConnectorWorkspace): calloc(4, sizeof(cy_expr*)), calloc(4, sizeof(cy_expr*)), ] self.abcd_elements[:] = [ - NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, - NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL - ] + NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, + NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL + ] def __dealloc__(self): cdef Py_ssize_t k, i @@ -194,7 +196,6 @@ cdef class BeamsplitterWorkspace(ConnectorWorkspace): cy_expr_init( self.sym_abcd_elements[k][2 * i + j], ch_sym, - self.sim ) @cython.boundscheck(False) @@ -214,6 +215,487 @@ cdef class BeamsplitterWorkspace(ConnectorWorkspace): self.sym_abcd_elements[k][i] ) + cdef void initialise_knm_workspaces(self) nogil: + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p1o = self.sim.trace[self.P1o_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + NodeBeamParam q_p2o = self.sim.trace[self.P2o_id] + NodeBeamParam q_p3i = self.sim.trace[self.P3i_id] + NodeBeamParam q_p3o = self.sim.trace[self.P3o_id] + NodeBeamParam q_p4i = self.sim.trace[self.P4i_id] + NodeBeamParam q_p4o = self.sim.trace[self.P4o_id] + + # Beam parameters for matched propagation of inputs + complex_t qx_p1o_matched_refl, qy_p1o_matched_refl + complex_t qx_p2o_matched_refl, qy_p2o_matched_refl + complex_t qx_p3o_matched_refl, qy_p3o_matched_refl + complex_t qx_p4o_matched_refl, qy_p4o_matched_refl + complex_t qx_p1o_matched_trns, qy_p1o_matched_trns + complex_t qx_p2o_matched_trns, qy_p2o_matched_trns + complex_t qx_p3o_matched_trns, qy_p3o_matched_trns + complex_t qx_p4o_matched_trns, qy_p4o_matched_trns + + double xbeta = self.bv.xbeta + double ybeta = self.bv.ybeta + double cos_alpha = self.cos_alpha + double cos_alpha_2 = self.cos_alpha_2 + + # Geometric factors for beam splitter beta angles giving + # gamma angle for each coupling + double beta_factor_12x = 2.0 + double beta_factor_12y = 2.0 * cos_alpha + double beta_factor_21x = beta_factor_12x + double beta_factor_21y = beta_factor_12y + double beta_factor_34x = beta_factor_12x + double beta_factor_34y = -2 * cos_alpha_2 + double beta_factor_43x = beta_factor_34x + double beta_factor_43y = beta_factor_34y + double beta_factor_13 = -(1 - self.nr1 / self.nr2) + double beta_factor_31 = 1 - self.nr2 / self.nr1 + double beta_factor_24 = -(1 - self.nr1 / self.nr2) + double beta_factor_42 = 1 - self.nr2 / self.nr1 + + # Flags for which direction is the "back-side" of a beam splitter + # giving total reflection + bint p1p2_total_refl = self.abcd_p1p2_x is None or self.abcd_p1p2_y is None + bint p2p1_total_refl = self.abcd_p2p1_x is None or self.abcd_p2p1_y is None + bint p3p4_total_refl = self.abcd_p3p4_x is None or self.abcd_p3p4_y is None + bint p4p3_total_refl = self.abcd_p4p3_x is None or self.abcd_p4p3_y is None + + double lambda0 = self.sim.model_data.lambda0 + int maxtem = self.sim.model_data.maxtem + + # Reflection at port p1 -> p2 + if p1p2_total_refl: + knm_ws_init_totalrefl(&self.K12_ws_x) + knm_ws_init_totalrefl(&self.K12_ws_y) + else: + qx_p2o_matched_refl = transform_q(self.abcd_p1p2_x, q_p1i.qx, self.nr1, self.nr1) + knm_ws_init( + &self.K12_ws_x, qx_p2o_matched_refl, q_p2o.qx, xbeta, beta_factor_12x, self.nr1, lambda0, maxtem, + ) + qy_p2o_matched_refl = transform_q(self.abcd_p1p2_y, q_p1i.qy, self.nr1, self.nr1) + knm_ws_init( + &self.K12_ws_y, qy_p2o_matched_refl, q_p2o.qy, ybeta, beta_factor_12y, self.nr1, lambda0, maxtem, + ) + + # Reflection at port p2 -> p1 + if p2p1_total_refl: + knm_ws_init_totalrefl(&self.K21_ws_x) + knm_ws_init_totalrefl(&self.K21_ws_y) + else: + qx_p1o_matched_refl = transform_q(self.abcd_p2p1_x, q_p2i.qx, self.nr1, self.nr1) + knm_ws_init( + &self.K21_ws_x, qx_p1o_matched_refl, q_p1o.qx, xbeta, beta_factor_21x, self.nr1, lambda0, maxtem, + ) + qy_p1o_matched_refl = transform_q(self.abcd_p2p1_y, q_p2i.qy, self.nr1, self.nr1) + knm_ws_init( + &self.K21_ws_y, qy_p1o_matched_refl, q_p1o.qy, ybeta, beta_factor_21y, self.nr1, lambda0, maxtem, + ) + + # Reflection at port p3 -> p4 + if p3p4_total_refl: + knm_ws_init_totalrefl(&self.K34_ws_x) + knm_ws_init_totalrefl(&self.K34_ws_y) + else: + qx_p4o_matched_refl = transform_q(self.abcd_p3p4_x, q_p3i.qx, self.nr2, self.nr2) + knm_ws_init( + &self.K34_ws_x, qx_p4o_matched_refl, q_p4o.qx, xbeta, beta_factor_34x, self.nr2, lambda0, maxtem, + ) + qy_p4o_matched_refl = transform_q(self.abcd_p3p4_y, q_p3i.qy, self.nr2, self.nr2) + knm_ws_init( + &self.K34_ws_y, qy_p4o_matched_refl, q_p4o.qy, ybeta, beta_factor_34y, self.nr2, lambda0, maxtem, + ) + + # Reflection at port p4 -> p3 + if p4p3_total_refl: + knm_ws_init_totalrefl(&self.K43_ws_x) + knm_ws_init_totalrefl(&self.K43_ws_y) + else: + qx_p3o_matched_refl = transform_q(self.abcd_p4p3_x, q_p4i.qx, self.nr2, self.nr2) + knm_ws_init( + &self.K43_ws_x, qx_p3o_matched_refl, q_p3o.qx, xbeta, beta_factor_43x, self.nr2, lambda0, maxtem, + ) + qy_p3o_matched_refl = transform_q(self.abcd_p4p3_y, q_p4i.qy, self.nr2, self.nr2) + knm_ws_init( + &self.K43_ws_y, qy_p3o_matched_refl, q_p3o.qy, ybeta, beta_factor_43y, self.nr2, lambda0, maxtem, + ) + + # Transmission at port p1 -> p3 + qx_p3o_matched_trns = transform_q(self.abcd_p1p3_x, q_p1i.qx, self.nr1, self.nr2) + knm_ws_init( + &self.K13_ws_x, qx_p3o_matched_trns, q_p3o.qx, xbeta, beta_factor_13, self.nr2, lambda0, maxtem, + ) + qy_p3o_matched_trns = transform_q(self.abcd_p1p3_y, q_p1i.qy, self.nr1, self.nr2) + knm_ws_init( + &self.K13_ws_y, qy_p3o_matched_trns, q_p3o.qy, ybeta, beta_factor_13, self.nr2, lambda0, maxtem, + ) + + # Transmission at port p3 -> p1 + qx_p1o_matched_trns = transform_q(self.abcd_p3p1_x, q_p3i.qx, self.nr2, self.nr1) + knm_ws_init( + &self.K31_ws_x, qx_p1o_matched_trns, q_p1o.qx, xbeta, -beta_factor_31, self.nr1, lambda0, maxtem, + ) + qy_p1o_matched_trns = transform_q(self.abcd_p3p1_y, q_p3i.qy, self.nr2, self.nr1) + knm_ws_init( + &self.K31_ws_y, qy_p1o_matched_trns, q_p1o.qy, ybeta, beta_factor_31, self.nr1, lambda0, maxtem, + ) + + # Transmission at port p2 -> p4 + qx_p4o_matched_trns = transform_q(self.abcd_p2p4_x, q_p2i.qx, self.nr1, self.nr2) + knm_ws_init( + &self.K24_ws_x, qx_p4o_matched_trns, q_p4o.qx, xbeta, beta_factor_24, self.nr2, lambda0, maxtem, + ) + qy_p4o_matched_trns = transform_q(self.abcd_p2p4_y, q_p2i.qy, self.nr1, self.nr2) + knm_ws_init( + &self.K24_ws_y, qy_p4o_matched_trns, q_p4o.qy, ybeta, beta_factor_24, self.nr2, lambda0, maxtem, + ) + + # Transmission at port p4 -> p2 + qx_p2o_matched_trns = transform_q(self.abcd_p4p2_x, q_p4i.qx, self.nr2, self.nr1) + knm_ws_init( + &self.K42_ws_x, qx_p2o_matched_trns, q_p2o.qx, xbeta, -beta_factor_42, self.nr1, lambda0, maxtem, + ) + qy_p2o_matched_trns = transform_q(self.abcd_p4p2_y, q_p4i.qy, self.nr2, self.nr1) + knm_ws_init( + &self.K42_ws_y, qy_p2o_matched_trns, q_p2o.qy, ybeta, beta_factor_42, self.nr1, lambda0, maxtem, + ) + + cdef void free_knm_workspaces(self) nogil: + knm_ws_free(&self.K12_ws_x) + knm_ws_free(&self.K12_ws_y) + knm_ws_free(&self.K21_ws_x) + knm_ws_free(&self.K21_ws_y) + knm_ws_free(&self.K34_ws_x) + knm_ws_free(&self.K34_ws_y) + knm_ws_free(&self.K43_ws_x) + knm_ws_free(&self.K43_ws_y) + knm_ws_free(&self.K13_ws_x) + knm_ws_free(&self.K13_ws_y) + knm_ws_free(&self.K31_ws_x) + knm_ws_free(&self.K31_ws_y) + knm_ws_free(&self.K24_ws_x) + knm_ws_free(&self.K24_ws_y) + knm_ws_free(&self.K42_ws_x) + knm_ws_free(&self.K42_ws_y) + + cdef void flag_changing_knm_workspaces(self): + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + NodeBeamParam q_p3i = self.sim.trace[self.P3i_id] + NodeBeamParam q_p4i = self.sim.trace[self.P4i_id] + + bint q_p1_or_p2_changing = q_p1i.is_changing or q_p2i.is_changing + bint q_p3_or_p4_changing = q_p3i.is_changing or q_p4i.is_changing + bint q_p1_or_p3_changing = q_p1i.is_changing or q_p3i.is_changing + bint q_p2_or_p4_changing = q_p2i.is_changing or q_p4i.is_changing + + bint xbeta_changing = self.owner.xbeta.is_changing + bint ybeta_changing = self.owner.ybeta.is_changing + + # Reflection at port p1 -> p2 + self.K12_ws_x.is_mm_changing = q_p1_or_p2_changing + self.K12_ws_x.is_alignment_changing = xbeta_changing + self.K12_ws_y.is_mm_changing = q_p1_or_p2_changing + self.K12_ws_y.is_alignment_changing = ybeta_changing + + if knm_ws_is_changing(&self.K12_ws_x) or knm_ws_is_changing(&self.K12_ws_y): + LOGGER.info("%s.K12 is changing", self.owner.name) + + # Reflection at port p2 -> p1 + self.K21_ws_x.is_mm_changing = q_p1_or_p2_changing + self.K21_ws_x.is_alignment_changing = xbeta_changing + self.K21_ws_y.is_mm_changing = q_p1_or_p2_changing + self.K21_ws_y.is_alignment_changing = ybeta_changing + + if knm_ws_is_changing(&self.K21_ws_x) or knm_ws_is_changing(&self.K21_ws_y): + LOGGER.info("%s.K21 is changing", self.owner.name) + + # Reflection at port p3 -> p4 + self.K34_ws_x.is_mm_changing = q_p3_or_p4_changing + self.K34_ws_x.is_alignment_changing = xbeta_changing + self.K34_ws_y.is_mm_changing = q_p3_or_p4_changing + self.K34_ws_y.is_alignment_changing = ybeta_changing + + if knm_ws_is_changing(&self.K34_ws_x) or knm_ws_is_changing(&self.K34_ws_y): + LOGGER.info("%s.K34 is changing", self.owner.name) + + # Reflection at port p4 -> p3 + self.K43_ws_x.is_mm_changing = q_p3_or_p4_changing + self.K43_ws_x.is_alignment_changing = xbeta_changing + self.K43_ws_y.is_mm_changing = q_p3_or_p4_changing + self.K43_ws_y.is_alignment_changing = ybeta_changing + + if knm_ws_is_changing(&self.K43_ws_x) or knm_ws_is_changing(&self.K43_ws_y): + LOGGER.info("%s.K43 is changing", self.owner.name) + + # Transmission at port p1 -> p3 + self.K13_ws_x.is_mm_changing = q_p1_or_p3_changing + self.K13_ws_x.is_alignment_changing = xbeta_changing + self.K13_ws_y.is_mm_changing = q_p1_or_p3_changing + self.K13_ws_y.is_alignment_changing = ybeta_changing + + if knm_ws_is_changing(&self.K13_ws_x) or knm_ws_is_changing(&self.K13_ws_y): + LOGGER.info("%s.K13 is changing", self.owner.name) + + # Transmission at port p3 -> p1 + self.K31_ws_x.is_mm_changing = q_p1_or_p3_changing + self.K31_ws_x.is_alignment_changing = xbeta_changing + self.K31_ws_y.is_mm_changing = q_p1_or_p3_changing + self.K31_ws_y.is_alignment_changing = ybeta_changing + + if knm_ws_is_changing(&self.K31_ws_x) or knm_ws_is_changing(&self.K31_ws_y): + LOGGER.info("%s.K31 is changing", self.owner.name) + + # Transmission at port p2 -> p4 + self.K24_ws_x.is_mm_changing = q_p2_or_p4_changing + self.K24_ws_x.is_alignment_changing = xbeta_changing + self.K24_ws_y.is_mm_changing = q_p2_or_p4_changing + self.K24_ws_y.is_alignment_changing = ybeta_changing + + if knm_ws_is_changing(&self.K24_ws_x) or knm_ws_is_changing(&self.K24_ws_y): + LOGGER.info("%s.K24 is changing", self.owner.name) + + # Transmission at port p4 -> p2 + self.K42_ws_x.is_mm_changing = q_p2_or_p4_changing + self.K42_ws_x.is_alignment_changing = xbeta_changing + self.K42_ws_y.is_mm_changing = q_p2_or_p4_changing + self.K42_ws_y.is_alignment_changing = ybeta_changing + + if knm_ws_is_changing(&self.K42_ws_x) or knm_ws_is_changing(&self.K42_ws_y): + LOGGER.info("%s.K42 is changing", self.owner.name) + + cdef void update_changing_knm_workspaces(self) nogil: + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p1o = self.sim.trace[self.P1o_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + NodeBeamParam q_p2o = self.sim.trace[self.P2o_id] + NodeBeamParam q_p3i = self.sim.trace[self.P3i_id] + NodeBeamParam q_p3o = self.sim.trace[self.P3o_id] + NodeBeamParam q_p4i = self.sim.trace[self.P4i_id] + NodeBeamParam q_p4o = self.sim.trace[self.P4o_id] + + complex_t qx_p1o_matched_refl, qy_p1o_matched_refl + complex_t qx_p2o_matched_refl, qy_p2o_matched_refl + complex_t qx_p3o_matched_refl, qy_p3o_matched_refl + complex_t qx_p4o_matched_refl, qy_p4o_matched_refl + complex_t qx_p1o_matched_trns, qy_p1o_matched_trns + complex_t qx_p2o_matched_trns, qy_p2o_matched_trns + complex_t qx_p3o_matched_trns, qy_p3o_matched_trns + complex_t qx_p4o_matched_trns, qy_p4o_matched_trns + + double xbeta = self.bv.xbeta + double ybeta = self.bv.ybeta + + # Reflection at port p1 -> p2 + if not self.K12_ws_x.total_reflection: + if self.K12_ws_x.is_mm_changing: + qx_p2o_matched_refl = transform_q(self.abcd_p1p2_x, q_p1i.qx, self.nr1, self.nr1) + knm_ws_recompute(&self.K12_ws_x, qx_p2o_matched_refl, q_p2o.qx, xbeta) + elif self.K12_ws_x.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K12_ws_x, xbeta) + + if not self.K12_ws_y.total_reflection: + if self.K12_ws_y.is_mm_changing: + qy_p2o_matched_refl = transform_q(self.abcd_p1p2_y, q_p1i.qy, self.nr1, self.nr1) + knm_ws_recompute(&self.K12_ws_y, qy_p2o_matched_refl, q_p2o.qy, ybeta) + elif self.K12_ws_y.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K12_ws_y, ybeta) + + # Reflection at port p2 -> p1 + if not self.K21_ws_x.total_reflection: + if self.K21_ws_x.is_mm_changing: + qx_p1o_matched_refl = transform_q(self.abcd_p2p1_x, q_p2i.qx, self.nr1, self.nr1) + knm_ws_recompute(&self.K21_ws_x, qx_p1o_matched_refl, q_p1o.qx, xbeta) + elif self.K21_ws_x.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K21_ws_x, xbeta) + + if not self.K21_ws_y.total_reflection: + if self.K21_ws_y.is_mm_changing: + qy_p1o_matched_refl = transform_q(self.abcd_p2p1_y, q_p2i.qy, self.nr1, self.nr1) + knm_ws_recompute(&self.K21_ws_y, qy_p1o_matched_refl, q_p1o.qy, ybeta) + elif self.K21_ws_y.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K21_ws_y, ybeta) + + # Reflection at port p3 -> p4 + if not self.K34_ws_x.total_reflection: + if self.K34_ws_x.is_mm_changing: + qx_p4o_matched_refl = transform_q(self.abcd_p3p4_x, q_p3i.qx, self.nr2, self.nr2) + knm_ws_recompute(&self.K34_ws_x, qx_p4o_matched_refl, q_p4o.qx, xbeta) + elif self.K34_ws_x.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K34_ws_x, xbeta) + + if not self.K34_ws_y.total_reflection: + if self.K34_ws_y.is_mm_changing: + qy_p4o_matched_refl = transform_q(self.abcd_p3p4_y, q_p3i.qy, self.nr2, self.nr2) + knm_ws_recompute(&self.K34_ws_y, qy_p4o_matched_refl, q_p4o.qy, ybeta) + elif self.K34_ws_y.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K34_ws_y, ybeta) + + # Reflection at port p4 -> p3 + if not self.K43_ws_x.total_reflection: + if self.K43_ws_x.is_mm_changing: + qx_p3o_matched_refl = transform_q(self.abcd_p4p3_x, q_p4i.qx, self.nr2, self.nr2) + knm_ws_recompute(&self.K43_ws_x, qx_p3o_matched_refl, q_p3o.qx, xbeta) + elif self.K43_ws_x.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K43_ws_x, xbeta) + + if not self.K43_ws_y.total_reflection: + if self.K43_ws_y.is_mm_changing: + qy_p3o_matched_refl = transform_q(self.abcd_p4p3_y, q_p4i.qy, self.nr2, self.nr2) + knm_ws_recompute(&self.K43_ws_y, qy_p3o_matched_refl, q_p3o.qy, ybeta) + elif self.K43_ws_y.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K43_ws_y, ybeta) + + # Transmission at port p1 -> p3 + if self.K13_ws_x.is_mm_changing: + qx_p3o_matched_trns = transform_q(self.abcd_p1p3_x, q_p1i.qx, self.nr1, self.nr2) + knm_ws_recompute(&self.K13_ws_x, qx_p3o_matched_trns, q_p3o.qx, xbeta) + elif self.K13_ws_x.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K13_ws_x, xbeta) + + if self.K13_ws_y.is_mm_changing: + qy_p3o_matched_trns = transform_q(self.abcd_p1p3_y, q_p1i.qy, self.nr1, self.nr2) + knm_ws_recompute(&self.K13_ws_y, qy_p3o_matched_trns, q_p3o.qy, ybeta) + elif self.K13_ws_y.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K13_ws_y, ybeta) + + # Transmission at port p3 -> p1 + if self.K31_ws_x.is_mm_changing: + qx_p1o_matched_trns = transform_q(self.abcd_p3p1_x, q_p3i.qx, self.nr2, self.nr1) + knm_ws_recompute(&self.K31_ws_x, qx_p1o_matched_trns, q_p1o.qx, xbeta) + elif self.K31_ws_x.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K31_ws_x, xbeta) + + if self.K31_ws_y.is_mm_changing: + qy_p1o_matched_trns = transform_q(self.abcd_p3p1_y, q_p3i.qy, self.nr2, self.nr1) + knm_ws_recompute(&self.K31_ws_y, qy_p1o_matched_trns, q_p1o.qy, ybeta) + elif self.K31_ws_y.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K31_ws_y, ybeta) + + # Transmission at port p2 -> p4 + if self.K24_ws_x.is_mm_changing: + qx_p4o_matched_trns = transform_q(self.abcd_p2p4_x, q_p2i.qx, self.nr1, self.nr2) + knm_ws_recompute(&self.K24_ws_x, qx_p4o_matched_trns, q_p4o.qx, xbeta) + elif self.K24_ws_x.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K24_ws_x, xbeta) + + if self.K24_ws_y.is_mm_changing: + qy_p4o_matched_trns = transform_q(self.abcd_p2p4_y, q_p2i.qy, self.nr1, self.nr2) + knm_ws_recompute(&self.K24_ws_y, qy_p4o_matched_trns, q_p4o.qy, ybeta) + elif self.K24_ws_y.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K24_ws_y, ybeta) + + # Transmission at port p4 -> p2 + if self.K42_ws_x.is_mm_changing: + qx_p2o_matched_trns = transform_q(self.abcd_p4p2_x, q_p4i.qx, self.nr2, self.nr1) + knm_ws_recompute(&self.K42_ws_x, qx_p2o_matched_trns, q_p2o.qx, xbeta) + elif self.K42_ws_x.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K42_ws_x, xbeta) + + if self.K42_ws_y.is_mm_changing: + qy_p2o_matched_trns = transform_q(self.abcd_p4p2_y, q_p4i.qy, self.nr2, self.nr1) + knm_ws_recompute(&self.K42_ws_y, qy_p2o_matched_trns, q_p2o.qy, ybeta) + elif self.K42_ws_y.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K42_ws_y, ybeta) + + + cdef void compute_scattering_matrices(self): + # Reflection p1 -> p2 + if not self.K12_ws_x.total_reflection and not self.K12_ws_y.total_reflection: + if knm_ws_is_changing(&self.K12_ws_x) or knm_ws_is_changing(&self.K12_ws_y): + compute_knm_matrix_bh(&self.K12_ws_x, &self.K12_ws_y, self.sim.homs_view, self.K12.data_view) + reverse_gouy_phases(self.K12.data_view, self.sim.homs_view, &self.K12_ws_x, &self.K12_ws_y, self.K12.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K12.data_view, self.K12.data_view) + + flip_odd_horizontal(self.K12.data_view, self.sim.homs_view, self.K12.data_view) + + knm_loss(self.K12.data_view, self.K12_loss) + + # Reflection p2 -> p1 + if not self.K21_ws_x.total_reflection and not self.K21_ws_y.total_reflection: + if knm_ws_is_changing(&self.K21_ws_x) or knm_ws_is_changing(&self.K21_ws_y): + compute_knm_matrix_bh(&self.K21_ws_x, &self.K21_ws_y, self.sim.homs_view, self.K21.data_view) + reverse_gouy_phases(self.K21.data_view, self.sim.homs_view, &self.K21_ws_x, &self.K21_ws_y, self.K21.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K21.data_view, self.K21.data_view) + + flip_odd_horizontal(self.K21.data_view, self.sim.homs_view, self.K21.data_view) + + knm_loss(self.K21.data_view, self.K21_loss) + + # Reflection p3 -> p4 + if not self.K34_ws_x.total_reflection and not self.K34_ws_y.total_reflection: + if knm_ws_is_changing(&self.K34_ws_x) or knm_ws_is_changing(&self.K34_ws_y): + compute_knm_matrix_bh(&self.K34_ws_x, &self.K34_ws_y, self.sim.homs_view, self.K34.data_view) + reverse_gouy_phases(self.K34.data_view, self.sim.homs_view, &self.K34_ws_x, &self.K34_ws_y, self.K34.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K34.data_view, self.K34.data_view) + + flip_odd_horizontal(self.K34.data_view, self.sim.homs_view, self.K34.data_view) + + knm_loss(self.K34.data_view, self.K34_loss) + + # Reflection p4 -> p3 + if not self.K43_ws_x.total_reflection and not self.K43_ws_y.total_reflection: + if knm_ws_is_changing(&self.K43_ws_x) or knm_ws_is_changing(&self.K43_ws_y): + compute_knm_matrix_bh(&self.K43_ws_x, &self.K43_ws_y, self.sim.homs_view, self.K43.data_view) + reverse_gouy_phases(self.K43.data_view, self.sim.homs_view, &self.K43_ws_x, &self.K43_ws_y, self.K43.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K43.data_view, self.K43.data_view) + + flip_odd_horizontal(self.K43.data_view, self.sim.homs_view, self.K43.data_view) + + knm_loss(self.K43.data_view, self.K43_loss) + + # Transmission p1 -> p3 + if knm_ws_is_changing(&self.K13_ws_x) or knm_ws_is_changing(&self.K13_ws_y): + compute_knm_matrix_bh(&self.K13_ws_x, &self.K13_ws_y, self.sim.homs_view, self.K13.data_view) + reverse_gouy_phases(self.K13.data_view, self.sim.homs_view, &self.K13_ws_x, &self.K13_ws_y, self.K13.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K13.data_view, self.K13.data_view) + + knm_loss(self.K13.data_view, self.K13_loss) + + # Transmission p3 -> p1 + if knm_ws_is_changing(&self.K31_ws_x) or knm_ws_is_changing(&self.K31_ws_y): + compute_knm_matrix_bh(&self.K31_ws_x, &self.K31_ws_y, self.sim.homs_view, self.K31.data_view) + reverse_gouy_phases(self.K31.data_view, self.sim.homs_view, &self.K31_ws_x, &self.K31_ws_y, self.K31.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K31.data_view, self.K31.data_view) + + knm_loss(self.K31.data_view, self.K31_loss) + + # Transmission p2 -> p4 + if knm_ws_is_changing(&self.K24_ws_x) or knm_ws_is_changing(&self.K24_ws_y): + compute_knm_matrix_bh(&self.K24_ws_x, &self.K24_ws_y, self.sim.homs_view, self.K24.data_view) + reverse_gouy_phases(self.K24.data_view, self.sim.homs_view, &self.K24_ws_x, &self.K24_ws_y, self.K24.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K24.data_view, self.K24.data_view) + + knm_loss(self.K24.data_view, self.K24_loss) + + # Transmission p4 -> p2 + if knm_ws_is_changing(&self.K42_ws_x) or knm_ws_is_changing(&self.K42_ws_y): + compute_knm_matrix_bh(&self.K42_ws_x, &self.K42_ws_y, self.sim.homs_view, self.K42.data_view) + reverse_gouy_phases(self.K42.data_view, self.sim.homs_view, &self.K42_ws_x, &self.K42_ws_y, self.K42.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K42.data_view, self.K42.data_view) + + knm_loss(self.K42.data_view, self.K42_loss) + beamsplitter_fill = FillFuncWrapper.make_from_ptr(c_beamsplitter_fill) cdef object c_beamsplitter_fill(ConnectorWorkspace cws): @@ -287,14 +769,13 @@ cdef object c_beamsplitter_fill(ConnectorWorkspace cws): double f0 = ws.sim.model_data.f0 double k0 = ws.sim.model_data.k0 - complex_t _it = 1.0j * sqrt(ws.bv.T) double _r = sqrt(ws.bv.R) + double _t = sqrt(ws.bv.T) + complex_t _r1, _r2, _t0 - double phi = radians(ws.bv.phi) - double phase_shift - - complex_t _tuning - complex_t _ctuning + double phi = DEG2RAD * ws.bv.phi + double alpha = DEG2RAD * ws.bv.alpha + double phase_shift_scaling double fval # current value of the frequency object double f_car # current value of the frequency object's carrier @@ -319,27 +800,38 @@ cdef object c_beamsplitter_fill(ConnectorWorkspace cws): bint is_lower_sb frequency_info_t freq - ws.cos_alpha = cos(DEG2RAD*ws.bv.alpha) - ws.cos_alpha_2 = cos(asin(ws.nr1 / ws.nr2 * sin(DEG2RAD * ws.bv.alpha))) + if ws.bv.alpha == 0.0: + ws.cos_alpha = 1 + ws.cos_alpha_2 = 1 + else: + ws.cos_alpha = cos(alpha) + ws.cos_alpha_2 = cos(asin(ws.nr1 / ws.nr2 * sin(alpha))) + + # Phase on reflection is not equal if nr1 != nr2 and AoI != 0 + # so the usual i on transmission phase no longer works. + # Uses N=-1, Eq.2.25 in Living Rev Relativ (2016) 19:3 DOI 10.1007/s41114-016-0002-8 + phi_r1 = 2 * phi * ws.cos_alpha * ws.nr1 + phi_r2 = -2 * phi * ws.cos_alpha_2 * ws.nr2 + phi_t = 0.5 * (phi_r1 + phi_r2) + (np.pi/2) for i in range(ws.sim.num_frequencies): freq = ws.sim.frequency_info[i] - phase_shift = 2.0 * phi * (1 + freq.f / f0) + phase_shift_scaling = (1 + freq.f / f0) + _r1 = _r * cexp(1j * phi_r1 * phase_shift_scaling) + _r2 = _r * cexp(1j * phi_r2 * phase_shift_scaling) + _t0 = _t * cexp(1j * phi_t * phase_shift_scaling) - _tuning = cexp(1.0j * phase_shift * ws.cos_alpha) - _ctuning = cexp(-1.0j * phase_shift * ws.cos_alpha_2) + (ws.bc.P1i_P2o.views[freq.index]).fill_za_zm_2(_r1, ws.K12.ptr, ws.K12.stride1, ws.K12.stride2) + (ws.bc.P2i_P1o.views[freq.index]).fill_za_zm_2(_r1, ws.K21.ptr, ws.K21.stride1, ws.K21.stride2) - (ws.bc.P1i_P2o.views[freq.index]).fill_za_zm_2(_r * _tuning, ws.K12.ptr, ws.K12.stride1, ws.K12.stride2) - (ws.bc.P2i_P1o.views[freq.index]).fill_za_zm_2(_r * _tuning, ws.K21.ptr, ws.K21.stride1, ws.K21.stride2) + (ws.bc.P3i_P4o.views[freq.index]).fill_za_zm_2(_r2, ws.K34.ptr, ws.K34.stride1, ws.K34.stride2) + (ws.bc.P4i_P3o.views[freq.index]).fill_za_zm_2(_r2, ws.K43.ptr, ws.K43.stride1, ws.K43.stride2) - (ws.bc.P3i_P4o.views[freq.index]).fill_za_zm_2(_r * _ctuning, ws.K34.ptr, ws.K34.stride1, ws.K34.stride2) - (ws.bc.P4i_P3o.views[freq.index]).fill_za_zm_2(_r * _ctuning, ws.K43.ptr, ws.K43.stride1, ws.K43.stride2) + (ws.bc.P1i_P3o.views[freq.index]).fill_za_zm_2(_t0, ws.K13.ptr, ws.K13.stride1, ws.K13.stride2) + (ws.bc.P3i_P1o.views[freq.index]).fill_za_zm_2(_t0, ws.K31.ptr, ws.K31.stride1, ws.K31.stride2) - (ws.bc.P1i_P3o.views[freq.index]).fill_za_zm_2(_it, ws.K13.ptr, ws.K13.stride1, ws.K13.stride2) - (ws.bc.P3i_P1o.views[freq.index]).fill_za_zm_2(_it, ws.K31.ptr, ws.K31.stride1, ws.K31.stride2) - - (ws.bc.P2i_P4o.views[freq.index]).fill_za_zm_2(_it, ws.K24.ptr, ws.K24.stride1, ws.K24.stride2) - (ws.bc.P4i_P2o.views[freq.index]).fill_za_zm_2(_it, ws.K42.ptr, ws.K42.stride1, ws.K42.stride2) + (ws.bc.P2i_P4o.views[freq.index]).fill_za_zm_2(_t0, ws.K24.ptr, ws.K24.stride1, ws.K24.stride2) + (ws.bc.P4i_P2o.views[freq.index]).fill_za_zm_2(_t0, ws.K42.ptr, ws.K42.stride1, ws.K42.stride2) if ws.sim.is_audio: # carrier = sim.DC @@ -526,485 +1018,3 @@ cdef object c_beamsplitter_fill(ConnectorWorkspace cws): ws.K34.ptr, ws.K34.stride1, ws.K34.stride2, c_p3_i, 1 # 1D contiguous array from outview ) - - -cdef void initialise_knm_workspaces(BeamsplitterWorkspace ws) nogil: - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p1o = ws.sim.trace[ws.P1o_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - NodeBeamParam q_p2o = ws.sim.trace[ws.P2o_id] - NodeBeamParam q_p3i = ws.sim.trace[ws.P3i_id] - NodeBeamParam q_p3o = ws.sim.trace[ws.P3o_id] - NodeBeamParam q_p4i = ws.sim.trace[ws.P4i_id] - NodeBeamParam q_p4o = ws.sim.trace[ws.P4o_id] - - # Beam parameters for matched propagation of inputs - complex_t qx_p1o_matched_refl, qy_p1o_matched_refl - complex_t qx_p2o_matched_refl, qy_p2o_matched_refl - complex_t qx_p3o_matched_refl, qy_p3o_matched_refl - complex_t qx_p4o_matched_refl, qy_p4o_matched_refl - complex_t qx_p1o_matched_trns, qy_p1o_matched_trns - complex_t qx_p2o_matched_trns, qy_p2o_matched_trns - complex_t qx_p3o_matched_trns, qy_p3o_matched_trns - complex_t qx_p4o_matched_trns, qy_p4o_matched_trns - - double xbeta = ws.bv.xbeta - double ybeta = ws.bv.ybeta - double cos_alpha = ws.cos_alpha - double cos_alpha_2 = ws.cos_alpha_2 - - # Geometric factors for beam splitter beta angles giving - # gamma angle for each coupling - double beta_factor_12x = 2.0 - double beta_factor_12y = 2.0 * cos_alpha - double beta_factor_21x = beta_factor_12x - double beta_factor_21y = beta_factor_12y - double beta_factor_34x = beta_factor_12x - double beta_factor_34y = -2 * cos_alpha_2 - double beta_factor_43x = beta_factor_34x - double beta_factor_43y = beta_factor_34y - double beta_factor_13 = -(1 - ws.nr1 / ws.nr2) - double beta_factor_31 = 1 - ws.nr2 / ws.nr1 - double beta_factor_24 = -(1 - ws.nr1 / ws.nr2) - double beta_factor_42 = 1 - ws.nr2 / ws.nr1 - - # Flags for which direction is the "back-side" of a beam splitter - # giving total reflection - bint p1p2_total_refl = ws.abcd_p1p2_x is None or ws.abcd_p1p2_y is None - bint p2p1_total_refl = ws.abcd_p2p1_x is None or ws.abcd_p2p1_y is None - bint p3p4_total_refl = ws.abcd_p3p4_x is None or ws.abcd_p3p4_y is None - bint p4p3_total_refl = ws.abcd_p4p3_x is None or ws.abcd_p4p3_y is None - - double lambda0 = ws.sim.model_data.lambda0 - int maxtem = ws.sim.model_data.maxtem - - # Reflection at port p1 -> p2 - if p1p2_total_refl: - knm_ws_init_totalrefl(&ws.K12_ws_x) - knm_ws_init_totalrefl(&ws.K12_ws_y) - else: - qx_p2o_matched_refl = transform_q(ws.abcd_p1p2_x, q_p1i.qx, ws.nr1, ws.nr1) - knm_ws_init( - &ws.K12_ws_x, qx_p2o_matched_refl, q_p2o.qx, xbeta, beta_factor_12x, ws.nr1, lambda0, maxtem, - ) - qy_p2o_matched_refl = transform_q(ws.abcd_p1p2_y, q_p1i.qy, ws.nr1, ws.nr1) - knm_ws_init( - &ws.K12_ws_y, qy_p2o_matched_refl, q_p2o.qy, ybeta, beta_factor_12y, ws.nr1, lambda0, maxtem, - ) - - # Reflection at port p2 -> p1 - if p2p1_total_refl: - knm_ws_init_totalrefl(&ws.K21_ws_x) - knm_ws_init_totalrefl(&ws.K21_ws_y) - else: - qx_p1o_matched_refl = transform_q(ws.abcd_p2p1_x, q_p2i.qx, ws.nr1, ws.nr1) - knm_ws_init( - &ws.K21_ws_x, qx_p1o_matched_refl, q_p1o.qx, xbeta, beta_factor_21x, ws.nr1, lambda0, maxtem, - ) - qy_p1o_matched_refl = transform_q(ws.abcd_p2p1_y, q_p2i.qy, ws.nr1, ws.nr1) - knm_ws_init( - &ws.K21_ws_y, qy_p1o_matched_refl, q_p1o.qy, ybeta, beta_factor_21y, ws.nr1, lambda0, maxtem, - ) - - # Reflection at port p3 -> p4 - if p3p4_total_refl: - knm_ws_init_totalrefl(&ws.K34_ws_x) - knm_ws_init_totalrefl(&ws.K34_ws_y) - else: - qx_p4o_matched_refl = transform_q(ws.abcd_p3p4_x, q_p3i.qx, ws.nr2, ws.nr2) - knm_ws_init( - &ws.K34_ws_x, qx_p4o_matched_refl, q_p4o.qx, xbeta, beta_factor_34x, ws.nr2, lambda0, maxtem, - ) - qy_p4o_matched_refl = transform_q(ws.abcd_p3p4_y, q_p3i.qy, ws.nr2, ws.nr2) - knm_ws_init( - &ws.K34_ws_y, qy_p4o_matched_refl, q_p4o.qy, ybeta, beta_factor_34y, ws.nr2, lambda0, maxtem, - ) - - # Reflection at port p4 -> p3 - if p4p3_total_refl: - knm_ws_init_totalrefl(&ws.K43_ws_x) - knm_ws_init_totalrefl(&ws.K43_ws_y) - else: - qx_p3o_matched_refl = transform_q(ws.abcd_p4p3_x, q_p4i.qx, ws.nr2, ws.nr2) - knm_ws_init( - &ws.K43_ws_x, qx_p3o_matched_refl, q_p3o.qx, xbeta, beta_factor_43x, ws.nr2, lambda0, maxtem, - ) - qy_p3o_matched_refl = transform_q(ws.abcd_p4p3_y, q_p4i.qy, ws.nr2, ws.nr2) - knm_ws_init( - &ws.K43_ws_y, qy_p3o_matched_refl, q_p3o.qy, ybeta, beta_factor_43y, ws.nr2, lambda0, maxtem, - ) - - # Transmission at port p1 -> p3 - qx_p3o_matched_trns = transform_q(ws.abcd_p1p3_x, q_p1i.qx, ws.nr1, ws.nr2) - knm_ws_init( - &ws.K13_ws_x, qx_p3o_matched_trns, q_p3o.qx, xbeta, beta_factor_13, ws.nr2, lambda0, maxtem, - ) - qy_p3o_matched_trns = transform_q(ws.abcd_p1p3_y, q_p1i.qy, ws.nr1, ws.nr2) - knm_ws_init( - &ws.K13_ws_y, qy_p3o_matched_trns, q_p3o.qy, ybeta, beta_factor_13, ws.nr2, lambda0, maxtem, - ) - - # Transmission at port p3 -> p1 - qx_p1o_matched_trns = transform_q(ws.abcd_p3p1_x, q_p3i.qx, ws.nr2, ws.nr1) - knm_ws_init( - &ws.K31_ws_x, qx_p1o_matched_trns, q_p1o.qx, xbeta, -beta_factor_31, ws.nr1, lambda0, maxtem, - ) - qy_p1o_matched_trns = transform_q(ws.abcd_p3p1_y, q_p3i.qy, ws.nr2, ws.nr1) - knm_ws_init( - &ws.K31_ws_y, qy_p1o_matched_trns, q_p1o.qy, ybeta, beta_factor_31, ws.nr1, lambda0, maxtem, - ) - - # Transmission at port p2 -> p4 - qx_p4o_matched_trns = transform_q(ws.abcd_p2p4_x, q_p2i.qx, ws.nr1, ws.nr2) - knm_ws_init( - &ws.K24_ws_x, qx_p4o_matched_trns, q_p4o.qx, xbeta, beta_factor_24, ws.nr2, lambda0, maxtem, - ) - qy_p4o_matched_trns = transform_q(ws.abcd_p2p4_y, q_p2i.qy, ws.nr1, ws.nr2) - knm_ws_init( - &ws.K24_ws_y, qy_p4o_matched_trns, q_p4o.qy, ybeta, beta_factor_24, ws.nr2, lambda0, maxtem, - ) - - # Transmission at port p4 -> p2 - qx_p2o_matched_trns = transform_q(ws.abcd_p4p2_x, q_p4i.qx, ws.nr2, ws.nr1) - knm_ws_init( - &ws.K42_ws_x, qx_p2o_matched_trns, q_p2o.qx, xbeta, -beta_factor_42, ws.nr1, lambda0, maxtem, - ) - qy_p2o_matched_trns = transform_q(ws.abcd_p4p2_y, q_p4i.qy, ws.nr2, ws.nr1) - knm_ws_init( - &ws.K42_ws_y, qy_p2o_matched_trns, q_p2o.qy, ybeta, beta_factor_42, ws.nr1, lambda0, maxtem, - ) - -cdef void free_knm_workspaces(BeamsplitterWorkspace ws) nogil: - knm_ws_free(&ws.K12_ws_x) - knm_ws_free(&ws.K12_ws_y) - knm_ws_free(&ws.K21_ws_x) - knm_ws_free(&ws.K21_ws_y) - knm_ws_free(&ws.K34_ws_x) - knm_ws_free(&ws.K34_ws_y) - knm_ws_free(&ws.K43_ws_x) - knm_ws_free(&ws.K43_ws_y) - knm_ws_free(&ws.K13_ws_x) - knm_ws_free(&ws.K13_ws_y) - knm_ws_free(&ws.K31_ws_x) - knm_ws_free(&ws.K31_ws_y) - knm_ws_free(&ws.K24_ws_x) - knm_ws_free(&ws.K24_ws_y) - knm_ws_free(&ws.K42_ws_x) - knm_ws_free(&ws.K42_ws_y) - -cdef void flag_changing_knm_workspaces(BeamsplitterWorkspace ws): - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - NodeBeamParam q_p3i = ws.sim.trace[ws.P3i_id] - NodeBeamParam q_p4i = ws.sim.trace[ws.P4i_id] - - bint q_p1_or_p2_changing = q_p1i.is_changing or q_p2i.is_changing - bint q_p3_or_p4_changing = q_p3i.is_changing or q_p4i.is_changing - bint q_p1_or_p3_changing = q_p1i.is_changing or q_p3i.is_changing - bint q_p2_or_p4_changing = q_p2i.is_changing or q_p4i.is_changing - - bint xbeta_changing = ws.owner.xbeta.is_changing - bint ybeta_changing = ws.owner.ybeta.is_changing - - # Reflection at port p1 -> p2 - ws.K12_ws_x.is_mm_changing = q_p1_or_p2_changing - ws.K12_ws_x.is_alignment_changing = xbeta_changing - ws.K12_ws_y.is_mm_changing = q_p1_or_p2_changing - ws.K12_ws_y.is_alignment_changing = ybeta_changing - - if knm_ws_is_changing(&ws.K12_ws_x) or knm_ws_is_changing(&ws.K12_ws_y): - LOGGER.info("%s.K12 is changing", ws.owner.name) - - # Reflection at port p2 -> p1 - ws.K21_ws_x.is_mm_changing = q_p1_or_p2_changing - ws.K21_ws_x.is_alignment_changing = xbeta_changing - ws.K21_ws_y.is_mm_changing = q_p1_or_p2_changing - ws.K21_ws_y.is_alignment_changing = ybeta_changing - - if knm_ws_is_changing(&ws.K21_ws_x) or knm_ws_is_changing(&ws.K21_ws_y): - LOGGER.info("%s.K21 is changing", ws.owner.name) - - # Reflection at port p3 -> p4 - ws.K34_ws_x.is_mm_changing = q_p3_or_p4_changing - ws.K34_ws_x.is_alignment_changing = xbeta_changing - ws.K34_ws_y.is_mm_changing = q_p3_or_p4_changing - ws.K34_ws_y.is_alignment_changing = ybeta_changing - - if knm_ws_is_changing(&ws.K34_ws_x) or knm_ws_is_changing(&ws.K34_ws_y): - LOGGER.info("%s.K34 is changing", ws.owner.name) - - # Reflection at port p4 -> p3 - ws.K43_ws_x.is_mm_changing = q_p3_or_p4_changing - ws.K43_ws_x.is_alignment_changing = xbeta_changing - ws.K43_ws_y.is_mm_changing = q_p3_or_p4_changing - ws.K43_ws_y.is_alignment_changing = ybeta_changing - - if knm_ws_is_changing(&ws.K43_ws_x) or knm_ws_is_changing(&ws.K43_ws_y): - LOGGER.info("%s.K43 is changing", ws.owner.name) - - # Transmission at port p1 -> p3 - ws.K13_ws_x.is_mm_changing = q_p1_or_p3_changing - ws.K13_ws_x.is_alignment_changing = xbeta_changing - ws.K13_ws_y.is_mm_changing = q_p1_or_p3_changing - ws.K13_ws_y.is_alignment_changing = ybeta_changing - - if knm_ws_is_changing(&ws.K13_ws_x) or knm_ws_is_changing(&ws.K13_ws_y): - LOGGER.info("%s.K13 is changing", ws.owner.name) - - # Transmission at port p3 -> p1 - ws.K31_ws_x.is_mm_changing = q_p1_or_p3_changing - ws.K31_ws_x.is_alignment_changing = xbeta_changing - ws.K31_ws_y.is_mm_changing = q_p1_or_p3_changing - ws.K31_ws_y.is_alignment_changing = ybeta_changing - - if knm_ws_is_changing(&ws.K31_ws_x) or knm_ws_is_changing(&ws.K31_ws_y): - LOGGER.info("%s.K31 is changing", ws.owner.name) - - # Transmission at port p2 -> p4 - ws.K24_ws_x.is_mm_changing = q_p2_or_p4_changing - ws.K24_ws_x.is_alignment_changing = xbeta_changing - ws.K24_ws_y.is_mm_changing = q_p2_or_p4_changing - ws.K24_ws_y.is_alignment_changing = ybeta_changing - - if knm_ws_is_changing(&ws.K24_ws_x) or knm_ws_is_changing(&ws.K24_ws_y): - LOGGER.info("%s.K24 is changing", ws.owner.name) - - # Transmission at port p4 -> p2 - ws.K42_ws_x.is_mm_changing = q_p2_or_p4_changing - ws.K42_ws_x.is_alignment_changing = xbeta_changing - ws.K42_ws_y.is_mm_changing = q_p2_or_p4_changing - ws.K42_ws_y.is_alignment_changing = ybeta_changing - - if knm_ws_is_changing(&ws.K42_ws_x) or knm_ws_is_changing(&ws.K42_ws_y): - LOGGER.info("%s.K42 is changing", ws.owner.name) - -cdef void update_changing_knm_workspaces(BeamsplitterWorkspace ws) nogil: - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p1o = ws.sim.trace[ws.P1o_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - NodeBeamParam q_p2o = ws.sim.trace[ws.P2o_id] - NodeBeamParam q_p3i = ws.sim.trace[ws.P3i_id] - NodeBeamParam q_p3o = ws.sim.trace[ws.P3o_id] - NodeBeamParam q_p4i = ws.sim.trace[ws.P4i_id] - NodeBeamParam q_p4o = ws.sim.trace[ws.P4o_id] - - complex_t qx_p1o_matched_refl, qy_p1o_matched_refl - complex_t qx_p2o_matched_refl, qy_p2o_matched_refl - complex_t qx_p3o_matched_refl, qy_p3o_matched_refl - complex_t qx_p4o_matched_refl, qy_p4o_matched_refl - complex_t qx_p1o_matched_trns, qy_p1o_matched_trns - complex_t qx_p2o_matched_trns, qy_p2o_matched_trns - complex_t qx_p3o_matched_trns, qy_p3o_matched_trns - complex_t qx_p4o_matched_trns, qy_p4o_matched_trns - - double xbeta = ws.bv.xbeta - double ybeta = ws.bv.ybeta - - # Reflection at port p1 -> p2 - if not ws.K12_ws_x.total_reflection: - if ws.K12_ws_x.is_mm_changing: - qx_p2o_matched_refl = transform_q(ws.abcd_p1p2_x, q_p1i.qx, ws.nr1, ws.nr1) - knm_ws_recompute(&ws.K12_ws_x, qx_p2o_matched_refl, q_p2o.qx, xbeta) - elif ws.K12_ws_x.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K12_ws_x, xbeta) - - if not ws.K12_ws_y.total_reflection: - if ws.K12_ws_y.is_mm_changing: - qy_p2o_matched_refl = transform_q(ws.abcd_p1p2_y, q_p1i.qy, ws.nr1, ws.nr1) - knm_ws_recompute(&ws.K12_ws_y, qy_p2o_matched_refl, q_p2o.qy, ybeta) - elif ws.K12_ws_y.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K12_ws_y, ybeta) - - # Reflection at port p2 -> p1 - if not ws.K21_ws_x.total_reflection: - if ws.K21_ws_x.is_mm_changing: - qx_p1o_matched_refl = transform_q(ws.abcd_p2p1_x, q_p2i.qx, ws.nr1, ws.nr1) - knm_ws_recompute(&ws.K21_ws_x, qx_p1o_matched_refl, q_p1o.qx, xbeta) - elif ws.K21_ws_x.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K21_ws_x, xbeta) - - if not ws.K21_ws_y.total_reflection: - if ws.K21_ws_y.is_mm_changing: - qy_p1o_matched_refl = transform_q(ws.abcd_p2p1_y, q_p2i.qy, ws.nr1, ws.nr1) - knm_ws_recompute(&ws.K21_ws_y, qy_p1o_matched_refl, q_p1o.qy, ybeta) - elif ws.K21_ws_y.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K21_ws_y, ybeta) - - # Reflection at port p3 -> p4 - if not ws.K34_ws_x.total_reflection: - if ws.K34_ws_x.is_mm_changing: - qx_p4o_matched_refl = transform_q(ws.abcd_p3p4_x, q_p3i.qx, ws.nr2, ws.nr2) - knm_ws_recompute(&ws.K34_ws_x, qx_p4o_matched_refl, q_p4o.qx, xbeta) - elif ws.K34_ws_x.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K34_ws_x, xbeta) - - if not ws.K34_ws_y.total_reflection: - if ws.K34_ws_y.is_mm_changing: - qy_p4o_matched_refl = transform_q(ws.abcd_p3p4_y, q_p3i.qy, ws.nr2, ws.nr2) - knm_ws_recompute(&ws.K34_ws_y, qy_p4o_matched_refl, q_p4o.qy, ybeta) - elif ws.K34_ws_y.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K34_ws_y, ybeta) - - # Reflection at port p4 -> p3 - if not ws.K43_ws_x.total_reflection: - if ws.K43_ws_x.is_mm_changing: - qx_p3o_matched_refl = transform_q(ws.abcd_p4p3_x, q_p4i.qx, ws.nr2, ws.nr2) - knm_ws_recompute(&ws.K43_ws_x, qx_p3o_matched_refl, q_p3o.qx, xbeta) - elif ws.K43_ws_x.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K43_ws_x, xbeta) - - if not ws.K43_ws_y.total_reflection: - if ws.K43_ws_y.is_mm_changing: - qy_p3o_matched_refl = transform_q(ws.abcd_p4p3_y, q_p4i.qy, ws.nr2, ws.nr2) - knm_ws_recompute(&ws.K43_ws_y, qy_p3o_matched_refl, q_p3o.qy, ybeta) - elif ws.K43_ws_y.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K43_ws_y, ybeta) - - # Transmission at port p1 -> p3 - if ws.K13_ws_x.is_mm_changing: - qx_p3o_matched_trns = transform_q(ws.abcd_p1p3_x, q_p1i.qx, ws.nr1, ws.nr2) - knm_ws_recompute(&ws.K13_ws_x, qx_p3o_matched_trns, q_p3o.qx, xbeta) - elif ws.K13_ws_x.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K13_ws_x, xbeta) - - if ws.K13_ws_y.is_mm_changing: - qy_p3o_matched_trns = transform_q(ws.abcd_p1p3_y, q_p1i.qy, ws.nr1, ws.nr2) - knm_ws_recompute(&ws.K13_ws_y, qy_p3o_matched_trns, q_p3o.qy, ybeta) - elif ws.K13_ws_y.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K13_ws_y, ybeta) - - # Transmission at port p3 -> p1 - if ws.K31_ws_x.is_mm_changing: - qx_p1o_matched_trns = transform_q(ws.abcd_p3p1_x, q_p3i.qx, ws.nr2, ws.nr1) - knm_ws_recompute(&ws.K31_ws_x, qx_p1o_matched_trns, q_p1o.qx, xbeta) - elif ws.K31_ws_x.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K31_ws_x, xbeta) - - if ws.K31_ws_y.is_mm_changing: - qy_p1o_matched_trns = transform_q(ws.abcd_p3p1_y, q_p3i.qy, ws.nr2, ws.nr1) - knm_ws_recompute(&ws.K31_ws_y, qy_p1o_matched_trns, q_p1o.qy, ybeta) - elif ws.K31_ws_y.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K31_ws_y, ybeta) - - # Transmission at port p2 -> p4 - if ws.K24_ws_x.is_mm_changing: - qx_p4o_matched_trns = transform_q(ws.abcd_p2p4_x, q_p2i.qx, ws.nr1, ws.nr2) - knm_ws_recompute(&ws.K24_ws_x, qx_p4o_matched_trns, q_p4o.qx, xbeta) - elif ws.K24_ws_x.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K24_ws_x, xbeta) - - if ws.K24_ws_y.is_mm_changing: - qy_p4o_matched_trns = transform_q(ws.abcd_p2p4_y, q_p2i.qy, ws.nr1, ws.nr2) - knm_ws_recompute(&ws.K24_ws_y, qy_p4o_matched_trns, q_p4o.qy, ybeta) - elif ws.K24_ws_y.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K24_ws_y, ybeta) - - # Transmission at port p4 -> p2 - if ws.K42_ws_x.is_mm_changing: - qx_p2o_matched_trns = transform_q(ws.abcd_p4p2_x, q_p4i.qx, ws.nr2, ws.nr1) - knm_ws_recompute(&ws.K42_ws_x, qx_p2o_matched_trns, q_p2o.qx, xbeta) - elif ws.K42_ws_x.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K42_ws_x, xbeta) - - if ws.K42_ws_y.is_mm_changing: - qy_p2o_matched_trns = transform_q(ws.abcd_p4p2_y, q_p4i.qy, ws.nr2, ws.nr1) - knm_ws_recompute(&ws.K42_ws_y, qy_p2o_matched_trns, q_p2o.qy, ybeta) - elif ws.K42_ws_y.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K42_ws_y, ybeta) - - -cdef void compute_scattering_matrices(BeamsplitterWorkspace ws): - # Reflection p1 -> p2 - if not ws.K12_ws_x.total_reflection and not ws.K12_ws_y.total_reflection: - if knm_ws_is_changing(&ws.K12_ws_x) or knm_ws_is_changing(&ws.K12_ws_y): - compute_knm_matrix_bh(&ws.K12_ws_x, &ws.K12_ws_y, ws.sim.homs_view, ws.K12.data_view) - reverse_gouy_phases(ws.K12.data_view, ws.sim.homs_view, &ws.K12_ws_x, &ws.K12_ws_y, ws.K12.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K12.data_view, ws.K12.data_view) - - flip_odd_horizontal(ws.K12.data_view, ws.sim.homs_view, ws.K12.data_view) - - knm_loss(ws.K12.data_view, ws.K12_loss) - - # Reflection p2 -> p1 - if not ws.K21_ws_x.total_reflection and not ws.K21_ws_y.total_reflection: - if knm_ws_is_changing(&ws.K21_ws_x) or knm_ws_is_changing(&ws.K21_ws_y): - compute_knm_matrix_bh(&ws.K21_ws_x, &ws.K21_ws_y, ws.sim.homs_view, ws.K21.data_view) - reverse_gouy_phases(ws.K21.data_view, ws.sim.homs_view, &ws.K21_ws_x, &ws.K21_ws_y, ws.K21.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K21.data_view, ws.K21.data_view) - - flip_odd_horizontal(ws.K21.data_view, ws.sim.homs_view, ws.K21.data_view) - - knm_loss(ws.K21.data_view, ws.K21_loss) - - # Reflection p3 -> p4 - if not ws.K34_ws_x.total_reflection and not ws.K34_ws_y.total_reflection: - if knm_ws_is_changing(&ws.K34_ws_x) or knm_ws_is_changing(&ws.K34_ws_y): - compute_knm_matrix_bh(&ws.K34_ws_x, &ws.K34_ws_y, ws.sim.homs_view, ws.K34.data_view) - reverse_gouy_phases(ws.K34.data_view, ws.sim.homs_view, &ws.K34_ws_x, &ws.K34_ws_y, ws.K34.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K34.data_view, ws.K34.data_view) - - flip_odd_horizontal(ws.K34.data_view, ws.sim.homs_view, ws.K34.data_view) - - knm_loss(ws.K34.data_view, ws.K34_loss) - - # Reflection p4 -> p3 - if not ws.K43_ws_x.total_reflection and not ws.K43_ws_y.total_reflection: - if knm_ws_is_changing(&ws.K43_ws_x) or knm_ws_is_changing(&ws.K43_ws_y): - compute_knm_matrix_bh(&ws.K43_ws_x, &ws.K43_ws_y, ws.sim.homs_view, ws.K43.data_view) - reverse_gouy_phases(ws.K43.data_view, ws.sim.homs_view, &ws.K43_ws_x, &ws.K43_ws_y, ws.K43.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K43.data_view, ws.K43.data_view) - - flip_odd_horizontal(ws.K43.data_view, ws.sim.homs_view, ws.K43.data_view) - - knm_loss(ws.K43.data_view, ws.K43_loss) - - # Transmission p1 -> p3 - if knm_ws_is_changing(&ws.K13_ws_x) or knm_ws_is_changing(&ws.K13_ws_y): - compute_knm_matrix_bh(&ws.K13_ws_x, &ws.K13_ws_y, ws.sim.homs_view, ws.K13.data_view) - reverse_gouy_phases(ws.K13.data_view, ws.sim.homs_view, &ws.K13_ws_x, &ws.K13_ws_y, ws.K13.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K13.data_view, ws.K13.data_view) - - knm_loss(ws.K13.data_view, ws.K13_loss) - - # Transmission p3 -> p1 - if knm_ws_is_changing(&ws.K31_ws_x) or knm_ws_is_changing(&ws.K31_ws_y): - compute_knm_matrix_bh(&ws.K31_ws_x, &ws.K31_ws_y, ws.sim.homs_view, ws.K31.data_view) - reverse_gouy_phases(ws.K31.data_view, ws.sim.homs_view, &ws.K31_ws_x, &ws.K31_ws_y, ws.K31.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K31.data_view, ws.K31.data_view) - - knm_loss(ws.K31.data_view, ws.K31_loss) - - # Transmission p2 -> p4 - if knm_ws_is_changing(&ws.K24_ws_x) or knm_ws_is_changing(&ws.K24_ws_y): - compute_knm_matrix_bh(&ws.K24_ws_x, &ws.K24_ws_y, ws.sim.homs_view, ws.K24.data_view) - reverse_gouy_phases(ws.K24.data_view, ws.sim.homs_view, &ws.K24_ws_x, &ws.K24_ws_y, ws.K24.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K24.data_view, ws.K24.data_view) - - knm_loss(ws.K24.data_view, ws.K24_loss) - - # Transmission p4 -> p2 - if knm_ws_is_changing(&ws.K42_ws_x) or knm_ws_is_changing(&ws.K42_ws_y): - compute_knm_matrix_bh(&ws.K42_ws_x, &ws.K42_ws_y, ws.sim.homs_view, ws.K42.data_view) - reverse_gouy_phases(ws.K42.data_view, ws.sim.homs_view, &ws.K42_ws_x, &ws.K42_ws_y, ws.K42.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K42.data_view, ws.K42.data_view) - - knm_loss(ws.K42.data_view, ws.K42_loss) diff --git a/src/finesse/components/modal/cavity.pxd b/src/finesse/components/modal/cavity.pxd index b7274f038a6dbbb498c9ffd6c8c4cd080f98fb5e..0fa80690dd42d1e27e635b6726b01526c2eadb5a 100644 --- a/src/finesse/components/modal/cavity.pxd +++ b/src/finesse/components/modal/cavity.pxd @@ -1,4 +1,11 @@ from finesse.cymath cimport complex_t +from finesse.cyexpr cimport ( + cy_expr, + cy_expr_new, + cy_expr_init, + cy_expr_free, + cy_expr_eval, +) from finesse.element cimport ElementWorkspace from finesse.tracing.ctracer cimport TraceTree @@ -9,46 +16,43 @@ cdef enum cavity_plane: cdef class CavityWorkspace(ElementWorkspace): - cdef readonly: + cdef: # Round-trip ABCDs in both planes double[:, ::1] abcd_x double[:, ::1] abcd_y + # The internal cavity tree (grabbed from the + # simulations' changing forest) + TraceTree tree + # Flag for whether the internal cavity tree is # in the changing forest of the simulation bint is_changing - bint is_stable # is 0 < gx, gy < 1 + bint is_stable # is 0 < {gx, gy} < 1 # Store both symbolic and numeric expressions for # the relevant cavity parameters - the symbolic # fields consist of only the changing symbols if any - object sym_length + cy_expr* sym_length double length - bint length_is_changing - object sym_loss + cy_expr* sym_loss double loss - bint loss_is_changing - object sym_finesse + cy_expr* sym_finesse double finesse - bint finesse_is_changing - object sym_fsr + cy_expr* sym_fsr double fsr - bint fsr_is_changing - object sym_fwhm + cy_expr* sym_fwhm double fwhm - bint fwhm_is_changing - object sym_tau + cy_expr* sym_tau double tau - bint tau_is_changing - object sym_pole + cy_expr* sym_pole double pole - bint pole_is_changing # These parameters are dependent upon the round-trip ABCD # so are computed differently (via the cavity trace tree) @@ -74,8 +78,6 @@ cdef class CavityWorkspace(ElementWorkspace): double gx double gy - cdef: - TraceTree tree cdef void __update_geometric(self) diff --git a/src/finesse/components/modal/cavity.pyx b/src/finesse/components/modal/cavity.pyx index 6fb3d35d18d17172a70673e81796abd1b333f8a3..382993cdb69540c34479facf638deacdf3606d78 100644 --- a/src/finesse/components/modal/cavity.pyx +++ b/src/finesse/components/modal/cavity.pyx @@ -40,57 +40,73 @@ cdef class CavityWorkspace(ElementWorkspace): self.abcd_x = owner._ABCDx self.abcd_y = owner._ABCDy - self.sym_length = owner._length.eval(keep_changing_symbols=True) - self.length_is_changing = isinstance(self.sym_length, Symbol) - if self.length_is_changing: - self.length = self.sym_length.eval() + ch_sym_length = owner._length.eval(keep_changing_symbols=True) + self.length = float(ch_sym_length) + if isinstance(ch_sym_length, Symbol): + self.sym_length = cy_expr_new() + cy_expr_init(self.sym_length, ch_sym_length) else: - self.length = self.sym_length + self.sym_length = NULL - self.sym_loss = owner._loss.eval(keep_changing_symbols=True) - self.loss_is_changing = isinstance(self.sym_loss, Symbol) - if self.loss_is_changing: - self.loss = self.sym_loss.eval() + ch_sym_loss = owner._loss.eval(keep_changing_symbols=True) + self.loss = float(ch_sym_loss) + if isinstance(ch_sym_loss, Symbol): + self.sym_loss = cy_expr_new() + cy_expr_init(self.sym_loss, ch_sym_loss) else: - self.loss = self.sym_loss + self.sym_loss = NULL - self.sym_finesse = owner._finesse.eval(keep_changing_symbols=True) - self.finesse_is_changing = isinstance(self.sym_finesse, Symbol) - if self.finesse_is_changing: - self.finesse = self.sym_finesse.eval() + ch_sym_finesse = owner._finesse.eval(keep_changing_symbols=True) + self.finesse = float(ch_sym_finesse) + if isinstance(ch_sym_finesse, Symbol): + self.sym_finesse = cy_expr_new() + cy_expr_init(self.sym_finesse, ch_sym_finesse) else: - self.finesse = self.sym_finesse + self.sym_finesse = NULL - self.sym_fsr = owner._FSR.eval(keep_changing_symbols=True) - self.fsr_is_changing = isinstance(self.sym_fsr, Symbol) - if self.fsr_is_changing: - self.fsr = self.sym_fsr.eval() + ch_sym_fsr = owner._FSR.eval(keep_changing_symbols=True) + self.fsr = float(ch_sym_fsr) + if isinstance(ch_sym_fsr, Symbol): + self.sym_fsr = cy_expr_new() + cy_expr_init(self.sym_fsr, ch_sym_fsr) else: - self.fsr = self.sym_fsr + self.sym_fsr = NULL - self.sym_fwhm = owner._FWHM.eval(keep_changing_symbols=True) - self.fwhm_is_changing = isinstance(self.sym_fwhm, Symbol) - if self.fwhm_is_changing: - self.fwhm = self.sym_fwhm.eval() + ch_sym_fwhm = owner._FWHM.eval(keep_changing_symbols=True) + self.fwhm = float(ch_sym_fwhm) + if isinstance(ch_sym_fwhm, Symbol): + self.sym_fwhm = cy_expr_new() + cy_expr_init(self.sym_fwhm, ch_sym_fwhm) else: - self.fwhm = self.sym_fwhm + self.sym_fwhm = NULL - self.sym_tau = owner._tau.eval(keep_changing_symbols=True) - self.tau_is_changing = isinstance(self.sym_tau, Symbol) - if self.tau_is_changing: - self.tau = self.sym_tau.eval() + ch_sym_tau = owner._tau.eval(keep_changing_symbols=True) + self.tau = float(ch_sym_tau) + if isinstance(ch_sym_tau, Symbol): + self.sym_tau = cy_expr_new() + cy_expr_init(self.sym_tau, ch_sym_tau) else: - self.tau = self.sym_tau + self.sym_tau = NULL - self.sym_pole = owner._pole.eval(keep_changing_symbols=True) - self.pole_is_changing = isinstance(self.sym_pole, Symbol) - if self.pole_is_changing: - self.pole = self.sym_pole.eval() + ch_sym_pole = owner._pole.eval(keep_changing_symbols=True) + self.pole = float(ch_sym_tau) + if isinstance(ch_sym_pole, Symbol): + self.sym_pole = cy_expr_new() + cy_expr_init(self.sym_pole, ch_sym_pole) else: - self.pole = self.sym_pole + self.sym_pole = NULL self.__update_geometric() + def __dealloc__(self): + cy_expr_free(self.sym_length) + cy_expr_free(self.sym_loss) + cy_expr_free(self.sym_finesse) + cy_expr_free(self.sym_fsr) + cy_expr_free(self.sym_fwhm) + cy_expr_free(self.sym_tau) + cy_expr_free(self.sym_pole) + cdef void __update_geometric(self): cdef: bint qx_ok, qy_ok @@ -127,20 +143,20 @@ cdef class CavityWorkspace(ElementWorkspace): cdef void update(self): # Update all the changing symbolic parameters - if self.length_is_changing: - self.length = self.sym_length.eval() - if self.loss_is_changing: - self.loss = self.sym_loss.eval() - if self.finesse_is_changing: - self.finesse = self.sym_finesse.eval() - if self.fsr_is_changing: - self.fsr = self.sym_fsr.eval() - if self.fwhm_is_changing: - self.fwhm = self.sym_fwhm.eval() - if self.tau_is_changing: - self.tau = self.sym_tau.eval() - if self.pole_is_changing: - self.pole = self.sym_pole.eval() + if self.sym_length != NULL: + self.length = cy_expr_eval(self.sym_length) + if self.sym_loss != NULL: + self.loss = cy_expr_eval(self.sym_loss) + if self.sym_finesse != NULL: + self.finesse = cy_expr_eval(self.sym_finesse) + if self.sym_fsr != NULL: + self.fsr = cy_expr_eval(self.sym_fsr) + if self.sym_fwhm != NULL: + self.fwhm = cy_expr_eval(self.sym_fwhm) + if self.sym_tau != NULL: + self.tau = cy_expr_eval(self.sym_tau) + if self.sym_pole != NULL: + self.pole = cy_expr_eval(self.sym_pole) # Update round-trip ABCD and dependent params if necessary if self.is_changing: diff --git a/src/finesse/components/modal/isolator.pxd b/src/finesse/components/modal/isolator.pxd index 3d153605dcc93403b4a82e2864ba37fffba79127..30384b46b4400c017aba3810c2de97ca55876ba1 100644 --- a/src/finesse/components/modal/isolator.pxd +++ b/src/finesse/components/modal/isolator.pxd @@ -3,7 +3,7 @@ from finesse.knm cimport KnmMatrix, KnmWorkspace from finesse.components.workspace cimport ConnectorWorkspace from finesse.element cimport BaseCValues from finesse.simulations.base cimport BaseSimulation - +from finesse.components.modal.workspace cimport KnmConnectorWorkspace from cpython.ref cimport PyObject @@ -29,7 +29,7 @@ cdef class IsolatorConnections: isolator_connections conn_ptrs -cdef class IsolatorWorkspace(ConnectorWorkspace): +cdef class IsolatorWorkspace(KnmConnectorWorkspace): cdef public: # Complete scattering matrices for each propagation direction KnmMatrix K12 @@ -55,8 +55,8 @@ cdef class IsolatorWorkspace(ConnectorWorkspace): KnmWorkspace K21_ws_y -cdef void initialise_knm_workspaces(IsolatorWorkspace ws) nogil -cdef void free_knm_workspaces(IsolatorWorkspace ws) nogil -cdef void flag_changing_knm_workspaces(IsolatorWorkspace ws) -cdef void update_changing_knm_workspaces(IsolatorWorkspace ws) nogil -cdef void compute_scattering_matrices(IsolatorWorkspace ws) + cdef void initialise_knm_workspaces(self) nogil + cdef void free_knm_workspaces(self) nogil + cdef void flag_changing_knm_workspaces(self) + cdef void update_changing_knm_workspaces(self) nogil + cdef void compute_scattering_matrices(self) diff --git a/src/finesse/components/modal/isolator.pyx b/src/finesse/components/modal/isolator.pyx index 9d6ab9634eb1f07f5822c9a5f5f5c5c6bb96279f..07a9f856e90f1586dee37ee714c95efca7a46a66 100644 --- a/src/finesse/components/modal/isolator.pyx +++ b/src/finesse/components/modal/isolator.pyx @@ -24,6 +24,7 @@ ctypedef (double*, ) ptr_tuple_1 LOGGER = logging.getLogger(__name__) +# TODO (sjr) make c_isolator_fill function? cdef class IsolatorValues(BaseCValues): def __init__(IsolatorValues self): @@ -41,8 +42,10 @@ cdef class IsolatorConnections: self.conn_ptrs.P2i_P1o = self.P2i_P1o.views -cdef class IsolatorWorkspace(ConnectorWorkspace): +cdef class IsolatorWorkspace(KnmConnectorWorkspace): def __init__(self, owner, BaseSimulation sim, refill): + self.str_couplings = ("12", "21") + super().__init__( owner, sim, @@ -55,109 +58,106 @@ cdef class IsolatorWorkspace(ConnectorWorkspace): self.P2i_id = sim.node_id(owner.p2.i) self.P2o_id = sim.node_id(owner.p2.o) -# TODO (sjr) make c_isolator_fill function? + cdef void initialise_knm_workspaces(self) nogil: + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p1o = self.sim.trace[self.P1o_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + NodeBeamParam q_p2o = self.sim.trace[self.P2o_id] -cdef void initialise_knm_workspaces(IsolatorWorkspace ws) nogil: - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p1o = ws.sim.trace[ws.P1o_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - NodeBeamParam q_p2o = ws.sim.trace[ws.P2o_id] - - # Beam parameters for matched propagation of inputs - complex_t qx_p1o_matched_trns, qy_p1o_matched_trns - complex_t qx_p2o_matched_trns, qy_p2o_matched_trns - - double lambda0 = ws.sim.model_data.lambda0 - int maxtem = ws.sim.model_data.maxtem - - # Transmission at port p1 -> p2 - qx_p2o_matched_trns = q_p1i.qx - knm_ws_init( - &ws.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx, 0, 0, ws.nr2, lambda0, maxtem, - ) - qy_p2o_matched_trns = q_p1i.qy - knm_ws_init( - &ws.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy, 0, 0, ws.nr2, lambda0, maxtem, - ) - - # Transmission at port p2 -> p1 - qx_p1o_matched_trns = q_p2i.qx - knm_ws_init( - &ws.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx, 0, 0, ws.nr1, lambda0, maxtem, - ) - qy_p1o_matched_trns = q_p2i.qy - knm_ws_init( - &ws.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy, 0, 0, ws.nr1, lambda0, maxtem, - ) - -cdef void free_knm_workspaces(IsolatorWorkspace ws) nogil: - knm_ws_free(&ws.K12_ws_x) - knm_ws_free(&ws.K12_ws_y) - knm_ws_free(&ws.K21_ws_x) - knm_ws_free(&ws.K21_ws_y) - -cdef void flag_changing_knm_workspaces(IsolatorWorkspace ws): - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - - bint is_changing = q_p1i.is_changing or q_p2i.is_changing - - # Transmission at port p1 -> p2 - ws.K12_ws_x.is_mm_changing = is_changing - ws.K12_ws_y.is_mm_changing = is_changing - - if knm_ws_is_changing(&ws.K12_ws_x) or knm_ws_is_changing(&ws.K12_ws_y): - LOGGER.info("%s.K12 is changing", ws.owner.name) - - # Transmission at port p2 -> p1 - ws.K21_ws_x.is_mm_changing = is_changing - ws.K21_ws_y.is_mm_changing = is_changing - - if knm_ws_is_changing(&ws.K21_ws_x) or knm_ws_is_changing(&ws.K21_ws_y): - LOGGER.info("%s.K21 is changing", ws.owner.name) - -cdef void update_changing_knm_workspaces(IsolatorWorkspace ws) nogil: - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p1o = ws.sim.trace[ws.P1o_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - NodeBeamParam q_p2o = ws.sim.trace[ws.P2o_id] - - complex_t qx_p1o_matched_trns, qy_p1o_matched_trns - complex_t qx_p2o_matched_trns, qy_p2o_matched_trns - - if ws.K12_ws_x.is_mm_changing: - qx_p2o_matched_trns = q_p1i.qx - knm_ws_recompute_mismatch(&ws.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx) + # Beam parameters for matched propagation of inputs + complex_t qx_p1o_matched_trns, qy_p1o_matched_trns + complex_t qx_p2o_matched_trns, qy_p2o_matched_trns + + double lambda0 = self.sim.model_data.lambda0 + int maxtem = self.sim.model_data.maxtem - if ws.K12_ws_y.is_mm_changing: + # Transmission at port p1 -> p2 + qx_p2o_matched_trns = q_p1i.qx + knm_ws_init( + &self.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx, 0, 0, self.nr2, lambda0, maxtem, + ) qy_p2o_matched_trns = q_p1i.qy - knm_ws_recompute_mismatch(&ws.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy) + knm_ws_init( + &self.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy, 0, 0, self.nr2, lambda0, maxtem, + ) - if ws.K21_ws_x.is_mm_changing: + # Transmission at port p2 -> p1 qx_p1o_matched_trns = q_p2i.qx - knm_ws_recompute_mismatch(&ws.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx) - - if ws.K21_ws_y.is_mm_changing: + knm_ws_init( + &self.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx, 0, 0, self.nr1, lambda0, maxtem, + ) qy_p1o_matched_trns = q_p2i.qy - knm_ws_recompute_mismatch(&ws.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy) + knm_ws_init( + &self.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy, 0, 0, self.nr1, lambda0, maxtem, + ) + + cdef void free_knm_workspaces(self) nogil: + knm_ws_free(&self.K12_ws_x) + knm_ws_free(&self.K12_ws_y) + knm_ws_free(&self.K21_ws_x) + knm_ws_free(&self.K21_ws_y) + + cdef void flag_changing_knm_workspaces(self): + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + + bint is_changing = q_p1i.is_changing or q_p2i.is_changing + + # Transmission at port p1 -> p2 + self.K12_ws_x.is_mm_changing = is_changing + self.K12_ws_y.is_mm_changing = is_changing + + if knm_ws_is_changing(&self.K12_ws_x) or knm_ws_is_changing(&self.K12_ws_y): + LOGGER.info("%s.K12 is changing", self.owner.name) + + # Transmission at port p2 -> p1 + self.K21_ws_x.is_mm_changing = is_changing + self.K21_ws_y.is_mm_changing = is_changing + + if knm_ws_is_changing(&self.K21_ws_x) or knm_ws_is_changing(&self.K21_ws_y): + LOGGER.info("%s.K21 is changing", self.owner.name) + + cdef void update_changing_knm_workspaces(self) nogil: + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p1o = self.sim.trace[self.P1o_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + NodeBeamParam q_p2o = self.sim.trace[self.P2o_id] + + complex_t qx_p1o_matched_trns, qy_p1o_matched_trns + complex_t qx_p2o_matched_trns, qy_p2o_matched_trns + + if self.K12_ws_x.is_mm_changing: + qx_p2o_matched_trns = q_p1i.qx + knm_ws_recompute_mismatch(&self.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx) + + if self.K12_ws_y.is_mm_changing: + qy_p2o_matched_trns = q_p1i.qy + knm_ws_recompute_mismatch(&self.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy) + + if self.K21_ws_x.is_mm_changing: + qx_p1o_matched_trns = q_p2i.qx + knm_ws_recompute_mismatch(&self.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx) + if self.K21_ws_y.is_mm_changing: + qy_p1o_matched_trns = q_p2i.qy + knm_ws_recompute_mismatch(&self.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy) -cdef void compute_scattering_matrices(IsolatorWorkspace ws): - # Transmission p1 -> p2 - if knm_ws_is_changing(&ws.K12_ws_x) or knm_ws_is_changing(&ws.K12_ws_y): - compute_knm_matrix_bh(&ws.K12_ws_x, &ws.K12_ws_y, ws.sim.homs_view, ws.K12.data_view) - reverse_gouy_phases(ws.K12.data_view, ws.sim.homs_view, &ws.K12_ws_x, &ws.K12_ws_y, ws.K12.data_view) + cdef void compute_scattering_matrices(self): + # Transmission p1 -> p2 + if knm_ws_is_changing(&self.K12_ws_x) or knm_ws_is_changing(&self.K12_ws_y): + compute_knm_matrix_bh(&self.K12_ws_x, &self.K12_ws_y, self.sim.homs_view, self.K12.data_view) + reverse_gouy_phases(self.K12.data_view, self.sim.homs_view, &self.K12_ws_x, &self.K12_ws_y, self.K12.data_view) - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K12.data_view, ws.K12.data_view) + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K12.data_view, self.K12.data_view) - # Transmission p2 -> p1 - if knm_ws_is_changing(&ws.K21_ws_x) or knm_ws_is_changing(&ws.K21_ws_y): - compute_knm_matrix_bh(&ws.K21_ws_x, &ws.K21_ws_y, ws.sim.homs_view, ws.K21.data_view) - reverse_gouy_phases(ws.K21.data_view, ws.sim.homs_view, &ws.K21_ws_x, &ws.K21_ws_y, ws.K21.data_view) + # Transmission p2 -> p1 + if knm_ws_is_changing(&self.K21_ws_x) or knm_ws_is_changing(&self.K21_ws_y): + compute_knm_matrix_bh(&self.K21_ws_x, &self.K21_ws_y, self.sim.homs_view, self.K21.data_view) + reverse_gouy_phases(self.K21.data_view, self.sim.homs_view, &self.K21_ws_x, &self.K21_ws_y, self.K21.data_view) - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K21.data_view, ws.K21.data_view) + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K21.data_view, self.K21.data_view) diff --git a/src/finesse/components/modal/lens.pxd b/src/finesse/components/modal/lens.pxd index 7c752a45669c594140cb2e5acfdbb9a73e69d5e7..84c7848f66d8d237022068ad14cad978419963ea 100644 --- a/src/finesse/components/modal/lens.pxd +++ b/src/finesse/components/modal/lens.pxd @@ -5,6 +5,7 @@ from finesse.knm cimport KnmMatrix, KnmWorkspace from finesse.components.workspace cimport ConnectorWorkspace from finesse.element cimport BaseCValues from finesse.simulations.base cimport BaseSimulation +from finesse.components.modal.workspace cimport KnmConnectorWorkspace from finesse.cyexpr cimport ( cy_expr, @@ -38,7 +39,7 @@ cdef class LensConnections: lens_connections conn_ptrs -cdef class LensWorkspace(ConnectorWorkspace): +cdef class LensWorkspace(KnmConnectorWorkspace): cdef public: # Complete scattering matrices for each propagation direction KnmMatrix K12 @@ -73,10 +74,8 @@ cdef class LensWorkspace(ConnectorWorkspace): cpdef compile_abcd_cy_exprs(self) cpdef update_parameter_values(self) - - -cdef void initialise_knm_workspaces(LensWorkspace ws) nogil -cdef void free_knm_workspaces(LensWorkspace ws) nogil -cdef void flag_changing_knm_workspaces(LensWorkspace ws) -cdef void update_changing_knm_workspaces(LensWorkspace ws) nogil -cdef void compute_scattering_matrices(LensWorkspace ws) + cdef void initialise_knm_workspaces(self) nogil + cdef void free_knm_workspaces(self) nogil + cdef void flag_changing_knm_workspaces(self) + cdef void update_changing_knm_workspaces(self) nogil + cdef void compute_scattering_matrices(self) diff --git a/src/finesse/components/modal/lens.pyx b/src/finesse/components/modal/lens.pyx index 454a7f6b42ec91de8ece8963939afea759015b14..4a242b04660dae4a47fa2fa200d0a0d793135ca9 100644 --- a/src/finesse/components/modal/lens.pyx +++ b/src/finesse/components/modal/lens.pyx @@ -45,8 +45,10 @@ cdef class LensConnections: self.conn_ptrs.P2i_P1o = self.P2i_P1o.views -cdef class LensWorkspace(ConnectorWorkspace): +cdef class LensWorkspace(KnmConnectorWorkspace): def __init__(self, owner, BaseSimulation sim, refill): + self.str_couplings = ("12", "21") + super().__init__( owner, sim, @@ -72,8 +74,10 @@ cdef class LensWorkspace(ConnectorWorkspace): cpdef compile_abcd_cy_exprs(self): cdef: dict abcd_handles = self.owner._abcd_matrices - cdef object[:, ::1] Mx_sym = abcd_handles["x"][0] - cdef object[:, ::1] My_sym = abcd_handles["y"][0] + tuple keyx = (self.owner.p1.i, self.owner.p2.o, "x") + tuple keyy = (self.owner.p1.i, self.owner.p2.o, "y") + cdef object[:, ::1] Mx_sym = abcd_handles[keyx][0] + cdef object[:, ::1] My_sym = abcd_handles[keyy][0] # NOTE (sjr) Only element C of a Lens ABCD matrix can possibly change @@ -81,13 +85,13 @@ cdef class LensWorkspace(ConnectorWorkspace): ch_sym = Mx_sym[1][0].eval(keep_changing_symbols=True) if isinstance(ch_sym, Symbol): self.sym_abcd_Cx = cy_expr_new() - cy_expr_init(self.sym_abcd_Cx, ch_sym, self.sim) + cy_expr_init(self.sym_abcd_Cx, ch_sym) if isinstance(My_sym[1][0], Symbol): ch_sym = My_sym[1][0].eval(keep_changing_symbols=True) if isinstance(ch_sym, Symbol): self.sym_abcd_Cy = cy_expr_new() - cy_expr_init(self.sym_abcd_Cy, ch_sym, self.sim) + cy_expr_init(self.sym_abcd_Cy, ch_sym) @cython.boundscheck(False) @cython.wraparound(False) @@ -101,110 +105,110 @@ cdef class LensWorkspace(ConnectorWorkspace): if self.sym_abcd_Cy != NULL: self.abcd_y[1][0] = cy_expr_eval(self.sym_abcd_Cy) + cdef void initialise_knm_workspaces(self) nogil: + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p1o = self.sim.trace[self.P1o_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + NodeBeamParam q_p2o = self.sim.trace[self.P2o_id] + + # Beam parameters for matched propagation of inputs + complex_t qx_p1o_matched_trns, qy_p1o_matched_trns + complex_t qx_p2o_matched_trns, qy_p2o_matched_trns + + double lambda0 = self.sim.model_data.lambda0 + int maxtem = self.sim.model_data.maxtem + + # Transmission at port p1 -> p2 + qx_p2o_matched_trns = transform_q(self.abcd_x, q_p1i.qx, self.nr1, self.nr2) + knm_ws_init( + &self.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx, 0, 0, self.nr2, lambda0, maxtem, + ) + qy_p2o_matched_trns = transform_q(self.abcd_y, q_p1i.qy, self.nr1, self.nr2) + knm_ws_init( + &self.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy, 0, 0, self.nr2, lambda0, maxtem, + ) + + # Transmission at port p2 -> p1 + qx_p1o_matched_trns = transform_q(self.abcd_x, q_p2i.qx, self.nr2, self.nr1) + knm_ws_init( + &self.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx, 0, 0, self.nr1, lambda0, maxtem, + ) + qy_p1o_matched_trns = transform_q(self.abcd_y, q_p2i.qy, self.nr2, self.nr1) + knm_ws_init( + &self.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy, 0, 0, self.nr1, lambda0, maxtem, + ) -# TODO (sjr) make c_lens_fill function? + cdef void free_knm_workspaces(self) nogil: + knm_ws_free(&self.K12_ws_x) + knm_ws_free(&self.K12_ws_y) + knm_ws_free(&self.K21_ws_x) + knm_ws_free(&self.K21_ws_y) + + cdef void flag_changing_knm_workspaces(self): + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + + bint is_changing = q_p1i.is_changing or q_p2i.is_changing + + # Transmission at port p1 -> p2 + self.K12_ws_x.is_mm_changing = is_changing + self.K12_ws_y.is_mm_changing = is_changing + + if knm_ws_is_changing(&self.K12_ws_x) or knm_ws_is_changing(&self.K12_ws_y): + LOGGER.info("%s.K12 is changing", self.owner.name) + + # Transmission at port p2 -> p1 + self.K21_ws_x.is_mm_changing = is_changing + self.K21_ws_y.is_mm_changing = is_changing + + if knm_ws_is_changing(&self.K21_ws_x) or knm_ws_is_changing(&self.K21_ws_y): + LOGGER.info("%s.K21 is changing", self.owner.name) -cdef void initialise_knm_workspaces(LensWorkspace ws) nogil: - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p1o = ws.sim.trace[ws.P1o_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - NodeBeamParam q_p2o = ws.sim.trace[ws.P2o_id] - - # Beam parameters for matched propagation of inputs - complex_t qx_p1o_matched_trns, qy_p1o_matched_trns - complex_t qx_p2o_matched_trns, qy_p2o_matched_trns - - double lambda0 = ws.sim.model_data.lambda0 - int maxtem = ws.sim.model_data.maxtem - - # Transmission at port p1 -> p2 - qx_p2o_matched_trns = transform_q(ws.abcd_x, q_p1i.qx, ws.nr1, ws.nr2) - knm_ws_init( - &ws.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx, 0, 0, ws.nr2, lambda0, maxtem, - ) - qy_p2o_matched_trns = transform_q(ws.abcd_y, q_p1i.qy, ws.nr1, ws.nr2) - knm_ws_init( - &ws.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy, 0, 0, ws.nr2, lambda0, maxtem, - ) - - # Transmission at port p2 -> p1 - qx_p1o_matched_trns = transform_q(ws.abcd_x, q_p2i.qx, ws.nr2, ws.nr1) - knm_ws_init( - &ws.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx, 0, 0, ws.nr1, lambda0, maxtem, - ) - qy_p1o_matched_trns = transform_q(ws.abcd_y, q_p2i.qy, ws.nr2, ws.nr1) - knm_ws_init( - &ws.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy, 0, 0, ws.nr1, lambda0, maxtem, - ) - -cdef void free_knm_workspaces(LensWorkspace ws) nogil: - knm_ws_free(&ws.K12_ws_x) - knm_ws_free(&ws.K12_ws_y) - knm_ws_free(&ws.K21_ws_x) - knm_ws_free(&ws.K21_ws_y) - -cdef void flag_changing_knm_workspaces(LensWorkspace ws): - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - - bint is_changing = q_p1i.is_changing or q_p2i.is_changing - - # Transmission at port p1 -> p2 - ws.K12_ws_x.is_mm_changing = is_changing - ws.K12_ws_y.is_mm_changing = is_changing - - if knm_ws_is_changing(&ws.K12_ws_x) or knm_ws_is_changing(&ws.K12_ws_y): - LOGGER.info("%s.K12 is changing", ws.owner.name) - - # Transmission at port p2 -> p1 - ws.K21_ws_x.is_mm_changing = is_changing - ws.K21_ws_y.is_mm_changing = is_changing - - if knm_ws_is_changing(&ws.K21_ws_x) or knm_ws_is_changing(&ws.K21_ws_y): - LOGGER.info("%s.K21 is changing", ws.owner.name) - -cdef void update_changing_knm_workspaces(LensWorkspace ws) nogil: - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p1o = ws.sim.trace[ws.P1o_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - NodeBeamParam q_p2o = ws.sim.trace[ws.P2o_id] - - complex_t qx_p1o_matched_trns, qy_p1o_matched_trns - complex_t qx_p2o_matched_trns, qy_p2o_matched_trns - - if ws.K12_ws_x.is_mm_changing: - qx_p2o_matched_trns = transform_q(ws.abcd_x, q_p1i.qx, ws.nr1, ws.nr2) - knm_ws_recompute_mismatch(&ws.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx) - - if ws.K12_ws_y.is_mm_changing: - qy_p2o_matched_trns = transform_q(ws.abcd_y, q_p1i.qy, ws.nr1, ws.nr2) - knm_ws_recompute_mismatch(&ws.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy) - - if ws.K21_ws_x.is_mm_changing: - qx_p1o_matched_trns = transform_q(ws.abcd_x, q_p2i.qx, ws.nr2, ws.nr1) - knm_ws_recompute_mismatch(&ws.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx) - - if ws.K21_ws_y.is_mm_changing: - qy_p1o_matched_trns = transform_q(ws.abcd_y, q_p2i.qy, ws.nr2, ws.nr1) - knm_ws_recompute_mismatch(&ws.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy) - - -cdef void compute_scattering_matrices(LensWorkspace ws): - # Transmission p1 -> p2 - if knm_ws_is_changing(&ws.K12_ws_x) or knm_ws_is_changing(&ws.K12_ws_y): - compute_knm_matrix_bh(&ws.K12_ws_x, &ws.K12_ws_y, ws.sim.homs_view, ws.K12.data_view) - reverse_gouy_phases(ws.K12.data_view, ws.sim.homs_view, &ws.K12_ws_x, &ws.K12_ws_y, ws.K12.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K12.data_view, ws.K12.data_view) - - # Transmission p2 -> p1 - if knm_ws_is_changing(&ws.K21_ws_x) or knm_ws_is_changing(&ws.K21_ws_y): - compute_knm_matrix_bh(&ws.K21_ws_x, &ws.K21_ws_y, ws.sim.homs_view, ws.K21.data_view) - reverse_gouy_phases(ws.K21.data_view, ws.sim.homs_view, &ws.K21_ws_x, &ws.K21_ws_y, ws.K21.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K21.data_view, ws.K21.data_view) + cdef void update_changing_knm_workspaces(self) nogil: + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p1o = self.sim.trace[self.P1o_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + NodeBeamParam q_p2o = self.sim.trace[self.P2o_id] + + complex_t qx_p1o_matched_trns, qy_p1o_matched_trns + complex_t qx_p2o_matched_trns, qy_p2o_matched_trns + + if self.K12_ws_x.is_mm_changing: + qx_p2o_matched_trns = transform_q(self.abcd_x, q_p1i.qx, self.nr1, self.nr2) + knm_ws_recompute_mismatch(&self.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx) + + if self.K12_ws_y.is_mm_changing: + qy_p2o_matched_trns = transform_q(self.abcd_y, q_p1i.qy, self.nr1, self.nr2) + knm_ws_recompute_mismatch(&self.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy) + + if self.K21_ws_x.is_mm_changing: + qx_p1o_matched_trns = transform_q(self.abcd_x, q_p2i.qx, self.nr2, self.nr1) + knm_ws_recompute_mismatch(&self.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx) + + if self.K21_ws_y.is_mm_changing: + qy_p1o_matched_trns = transform_q(self.abcd_y, q_p2i.qy, self.nr2, self.nr1) + knm_ws_recompute_mismatch(&self.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy) + + + cdef void compute_scattering_matrices(self): + # Transmission p1 -> p2 + if knm_ws_is_changing(&self.K12_ws_x) or knm_ws_is_changing(&self.K12_ws_y): + compute_knm_matrix_bh(&self.K12_ws_x, &self.K12_ws_y, self.sim.homs_view, self.K12.data_view) + reverse_gouy_phases(self.K12.data_view, self.sim.homs_view, &self.K12_ws_x, &self.K12_ws_y, self.K12.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K12.data_view, self.K12.data_view) + + # Transmission p2 -> p1 + if knm_ws_is_changing(&self.K21_ws_x) or knm_ws_is_changing(&self.K21_ws_y): + compute_knm_matrix_bh(&self.K21_ws_x, &self.K21_ws_y, self.sim.homs_view, self.K21.data_view) + reverse_gouy_phases(self.K21.data_view, self.sim.homs_view, &self.K21_ws_x, &self.K21_ws_y, self.K21.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K21.data_view, self.K21.data_view) + + +# TODO (sjr) make c_lens_fill function? diff --git a/src/finesse/components/modal/mirror.pxd b/src/finesse/components/modal/mirror.pxd index f696a239c31f2688382f042a8aa0c968502c8554..407706ab18885d61b59261e2d6f8e6722d84244e 100644 --- a/src/finesse/components/modal/mirror.pxd +++ b/src/finesse/components/modal/mirror.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, @@ -66,7 +67,7 @@ cdef class MirrorConnections: cdef: mirror_connections conn_ptrs -cdef class MirrorWorkspace(ConnectorWorkspace): +cdef class MirrorWorkspace(KnmConnectorWorkspace): cdef public: complex_t field_to_Fz complex_t z_to_field @@ -141,9 +142,8 @@ cdef class MirrorWorkspace(ConnectorWorkspace): cpdef compile_abcd_cy_exprs(self) cpdef update_parameter_values(self) - -cdef void initialise_knm_workspaces(MirrorWorkspace ws) nogil -cdef void free_knm_workspaces(MirrorWorkspace ws) nogil -cdef void flag_changing_knm_workspaces(MirrorWorkspace ws) -cdef void update_changing_knm_workspaces(MirrorWorkspace ws) nogil -cdef void compute_scattering_matrices(MirrorWorkspace ws) + cdef void free_knm_workspaces(self) nogil + cdef void initialise_knm_workspaces(self) nogil + cdef void flag_changing_knm_workspaces(self) + cdef void update_changing_knm_workspaces(self) nogil + cdef void compute_scattering_matrices(self) diff --git a/src/finesse/components/modal/mirror.pyx b/src/finesse/components/modal/mirror.pyx index 2c6b3ff9862a7db9acbc33c127a79f4d3c539b46..c80f55b3a1898126d5004001e879dee851bca5b1 100644 --- a/src/finesse/components/modal/mirror.pyx +++ b/src/finesse/components/modal/mirror.pyx @@ -81,8 +81,10 @@ cdef class MirrorConnections: self.conn_ptrs.Z_P2o = self.Z_P2o.views -cdef class MirrorWorkspace(ConnectorWorkspace): +cdef class MirrorWorkspace(KnmConnectorWorkspace): def __init__(self, owner, BaseSimulation sim, refill): + self.str_couplings = ("11", "12", "21", "22") + super().__init__( owner, sim, @@ -142,7 +144,7 @@ cdef class MirrorWorkspace(ConnectorWorkspace): ch_sym = M_sym[1][0].eval(keep_changing_symbols=True) if isinstance(ch_sym, Symbol): self.sym_abcd_Cs[i] = cy_expr_new() - cy_expr_init(self.sym_abcd_Cs[i], ch_sym, self.sim) + cy_expr_init(self.sym_abcd_Cs[i], ch_sym) @cython.boundscheck(False) @cython.wraparound(False) @@ -155,6 +157,236 @@ cdef class MirrorWorkspace(ConnectorWorkspace): if self.sym_abcd_Cs[i] != NULL and self.abcd_Cs[i] != NULL: self.abcd_Cs[i][0] = cy_expr_eval(self.sym_abcd_Cs[i]) + cdef void free_knm_workspaces(self) nogil: + knm_ws_free(&self.K11_ws_x) + knm_ws_free(&self.K11_ws_y) + knm_ws_free(&self.K22_ws_x) + knm_ws_free(&self.K22_ws_y) + knm_ws_free(&self.K12_ws_x) + knm_ws_free(&self.K12_ws_y) + knm_ws_free(&self.K21_ws_x) + knm_ws_free(&self.K21_ws_y) + + cdef void flag_changing_knm_workspaces(self): + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + + bint q_p1_or_p2_changing = q_p1i.is_changing or q_p2i.is_changing + bint xbeta_changing = self.owner.xbeta.is_changing + bint ybeta_changing = self.owner.ybeta.is_changing + + # Reflection at port p1 -> p1 + self.K11_ws_x.is_mm_changing = q_p1i.is_changing + self.K11_ws_x.is_alignment_changing = xbeta_changing + self.K11_ws_y.is_mm_changing = q_p1i.is_changing + self.K11_ws_y.is_alignment_changing = ybeta_changing + + if knm_ws_is_changing(&self.K11_ws_x) or knm_ws_is_changing(&self.K11_ws_y): + LOGGER.info("%s.K11 is changing", self.owner.name) + + # Reflection at port p2 -> p2 + self.K22_ws_x.is_mm_changing = q_p2i.is_changing + self.K22_ws_x.is_alignment_changing = xbeta_changing + self.K22_ws_y.is_mm_changing = q_p2i.is_changing + self.K22_ws_y.is_alignment_changing = ybeta_changing + + if knm_ws_is_changing(&self.K22_ws_x) or knm_ws_is_changing(&self.K22_ws_y): + LOGGER.info("%s.K22 is changing", self.owner.name) + + # Transmission at port p1 -> p2 + self.K12_ws_x.is_mm_changing = q_p1_or_p2_changing + self.K12_ws_x.is_alignment_changing = xbeta_changing + self.K12_ws_y.is_mm_changing = q_p1_or_p2_changing + self.K12_ws_y.is_alignment_changing = ybeta_changing + + if knm_ws_is_changing(&self.K12_ws_x) or knm_ws_is_changing(&self.K12_ws_y): + LOGGER.info("%s.K12 is changing", self.owner.name) + + # Transmission at port p2 -> p1 + self.K21_ws_x.is_mm_changing = q_p1_or_p2_changing + self.K21_ws_x.is_alignment_changing = xbeta_changing + self.K21_ws_y.is_mm_changing = q_p1_or_p2_changing + self.K21_ws_y.is_alignment_changing = ybeta_changing + + if knm_ws_is_changing(&self.K21_ws_x) or knm_ws_is_changing(&self.K21_ws_y): + LOGGER.info("%s.K21 is changing", self.owner.name) + + cdef void update_changing_knm_workspaces(self) nogil: + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p1o = self.sim.trace[self.P1o_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + NodeBeamParam q_p2o = self.sim.trace[self.P2o_id] + + complex_t qx_p1o_matched_refl, qy_p1o_matched_refl + complex_t qx_p2o_matched_refl, qy_p2o_matched_refl + complex_t qx_p1o_matched_trns, qy_p1o_matched_trns + complex_t qx_p2o_matched_trns, qy_p2o_matched_trns + + double xbeta = self.mv.xbeta + double ybeta = self.mv.ybeta + + if self.K11_ws_x.is_mm_changing: + qx_p1o_matched_refl = transform_q(self.abcd_p1p1_x, q_p1i.qx, self.nr1, self.nr1) + knm_ws_recompute(&self.K11_ws_x, qx_p1o_matched_refl, q_p1o.qx, xbeta) + elif self.K11_ws_x.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K11_ws_x, xbeta) + + if self.K11_ws_y.is_mm_changing: + qy_p1o_matched_refl = transform_q(self.abcd_p1p1_y, q_p1i.qy, self.nr1, self.nr1) + knm_ws_recompute(&self.K11_ws_y, qy_p1o_matched_refl, q_p1o.qy, ybeta) + elif self.K11_ws_y.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K11_ws_y, ybeta) + + if self.K22_ws_x.is_mm_changing: + qx_p2o_matched_refl = transform_q(self.abcd_p2p2_x, q_p2i.qx, self.nr2, self.nr2) + knm_ws_recompute(&self.K22_ws_x, qx_p2o_matched_refl, q_p2o.qx, xbeta) + elif self.K22_ws_x.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K22_ws_x, xbeta) + + if self.K22_ws_y.is_mm_changing: + qy_p2o_matched_refl = transform_q(self.abcd_p2p2_y, q_p2i.qy, self.nr2, self.nr2) + knm_ws_recompute(&self.K22_ws_y, qy_p2o_matched_refl, q_p2o.qy, ybeta) + elif self.K22_ws_y.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K22_ws_y, ybeta) + + if self.K12_ws_x.is_mm_changing: + qx_p2o_matched_trns = transform_q(self.abcd_p1p2_x, q_p1i.qx, self.nr1, self.nr2) + knm_ws_recompute(&self.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx, xbeta) + elif self.K12_ws_x.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K12_ws_x, xbeta) + + if self.K12_ws_y.is_mm_changing: + qy_p2o_matched_trns = transform_q(self.abcd_p1p2_y, q_p1i.qy, self.nr1, self.nr2) + knm_ws_recompute(&self.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy, ybeta) + elif self.K12_ws_y.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K12_ws_y, ybeta) + + if self.K21_ws_x.is_mm_changing: + qx_p1o_matched_trns = transform_q(self.abcd_p2p1_x, q_p2i.qx, self.nr2, self.nr1) + knm_ws_recompute(&self.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx, xbeta) + elif self.K21_ws_x.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K21_ws_x, xbeta) + + if self.K21_ws_y.is_mm_changing: + qy_p1o_matched_trns = transform_q(self.abcd_p2p1_y, q_p2i.qy, self.nr2, self.nr1) + knm_ws_recompute(&self.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy, ybeta) + elif self.K21_ws_y.is_alignment_changing: + knm_ws_recompute_misalignment(&self.K21_ws_y, ybeta) + + + cdef void compute_scattering_matrices(self): + # Reflection p1 -> p1 + if knm_ws_is_changing(&self.K11_ws_x) or knm_ws_is_changing(&self.K11_ws_y): + compute_knm_matrix_bh(&self.K11_ws_x, &self.K11_ws_y, self.sim.homs_view, self.K11.data_view) + reverse_gouy_phases(self.K11.data_view, self.sim.homs_view, &self.K11_ws_x, &self.K11_ws_y, self.K11.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K11.data_view, self.K11.data_view) + + flip_odd_horizontal(self.K11.data_view, self.sim.homs_view, self.K11.data_view) + + knm_loss(self.K11.data_view, self.K11_loss) + + # Reflection p2 -> p2 + if knm_ws_is_changing(&self.K22_ws_x) or knm_ws_is_changing(&self.K22_ws_y): + compute_knm_matrix_bh(&self.K22_ws_x, &self.K22_ws_y, self.sim.homs_view, self.K22.data_view) + reverse_gouy_phases(self.K22.data_view, self.sim.homs_view, &self.K22_ws_x, &self.K22_ws_y, self.K22.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K22.data_view, self.K22.data_view) + + flip_odd_horizontal(self.K22.data_view, self.sim.homs_view, self.K22.data_view) + + knm_loss(self.K22.data_view, self.K22_loss) + + # Transmission p1 -> p2 + if knm_ws_is_changing(&self.K12_ws_x) or knm_ws_is_changing(&self.K12_ws_y): + compute_knm_matrix_bh(&self.K12_ws_x, &self.K12_ws_y, self.sim.homs_view, self.K12.data_view) + reverse_gouy_phases(self.K12.data_view, self.sim.homs_view, &self.K12_ws_x, &self.K12_ws_y, self.K12.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K12.data_view, self.K12.data_view) + + knm_loss(self.K12.data_view, self.K12_loss) + + # Transmission p2 -> p1 + if knm_ws_is_changing(&self.K21_ws_x) or knm_ws_is_changing(&self.K21_ws_y): + compute_knm_matrix_bh(&self.K21_ws_x, &self.K21_ws_y, self.sim.homs_view, self.K21.data_view) + reverse_gouy_phases(self.K21.data_view, self.sim.homs_view, &self.K21_ws_x, &self.K21_ws_y, self.K21.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K21.data_view, self.K21.data_view) + + knm_loss(self.K21.data_view, self.K21_loss) + + cdef void initialise_knm_workspaces(self) nogil: + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p1o = self.sim.trace[self.P1o_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + NodeBeamParam q_p2o = self.sim.trace[self.P2o_id] + + # Beam parameters for matched propagation of inputs + complex_t qx_p1o_matched_refl, qy_p1o_matched_refl + complex_t qx_p2o_matched_refl, qy_p2o_matched_refl + complex_t qx_p1o_matched_trns, qy_p1o_matched_trns + complex_t qx_p2o_matched_trns, qy_p2o_matched_trns + + double xbeta = self.mv.xbeta + double ybeta = self.mv.ybeta + + # Geometric factors for mirror beta angles giving + # gamma angle for each coupling + double beta_factor_11 = 2.0 + double beta_factor_22 = 2.0 + double beta_factor_12 = -(1 - self.nr1 / self.nr2) + double beta_factor_21 = 1 - self.nr2 / self.nr1 + + double lambda0 = self.sim.model_data.lambda0 + int maxtem = self.sim.model_data.maxtem + + # Reflection at port p1 -> p1 + qx_p1o_matched_refl = transform_q(self.abcd_p1p1_x, q_p1i.qx, self.nr1, self.nr1) + knm_ws_init( + &self.K11_ws_x, qx_p1o_matched_refl, q_p1o.qx, xbeta, beta_factor_11, self.nr1, lambda0, maxtem, + ) + qy_p1o_matched_refl = transform_q(self.abcd_p1p1_y, q_p1i.qy, self.nr1, self.nr1) + knm_ws_init( + &self.K11_ws_y, qy_p1o_matched_refl, q_p1o.qy, ybeta, beta_factor_11, self.nr1, lambda0, maxtem, + ) + + # Reflection at port p2 -> p2 + qx_p2o_matched_refl = transform_q(self.abcd_p2p2_x, q_p2i.qx, self.nr2, self.nr2) + knm_ws_init( + &self.K22_ws_x, qx_p2o_matched_refl, q_p2o.qx, xbeta, beta_factor_22, self.nr2, lambda0, maxtem, + ) + qy_p2o_matched_refl = transform_q(self.abcd_p2p2_y, q_p2i.qy, self.nr2, self.nr2) + knm_ws_init( + &self.K22_ws_y, qy_p2o_matched_refl, q_p2o.qy, ybeta, -beta_factor_22, self.nr2, lambda0, maxtem, + ) + + # Transmission at port p1 -> p2 + qx_p2o_matched_trns = transform_q(self.abcd_p1p2_x, q_p1i.qx, self.nr1, self.nr2) + knm_ws_init( + &self.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx, xbeta, beta_factor_12, self.nr2, lambda0, maxtem, + ) + qy_p2o_matched_trns = transform_q(self.abcd_p1p2_y, q_p1i.qy, self.nr1, self.nr2) + knm_ws_init( + &self.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy, ybeta, beta_factor_12, self.nr2, lambda0, maxtem, + ) + + # Transmission at port p2 -> p1 + qx_p1o_matched_trns = transform_q(self.abcd_p2p1_x, q_p2i.qx, self.nr2, self.nr1) + knm_ws_init( + &self.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx, xbeta, -beta_factor_21, self.nr1, lambda0, maxtem, + ) + qy_p1o_matched_trns = transform_q(self.abcd_p2p1_y, q_p2i.qy, self.nr2, self.nr1) + knm_ws_init( + &self.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy, ybeta, beta_factor_21, self.nr1, lambda0, maxtem, + ) + mirror_fill = FillFuncWrapper.make_from_ptr(c_mirror_fill) cdef object c_mirror_fill(ConnectorWorkspace cws): @@ -403,236 +635,3 @@ cdef void single_z_mechanical_frequency_signal_calc ( ws.K22.ptr, ws.K22.stride1, ws.K22.stride2, c_p2_i, 1 # 1D contiguous array from outview ) - - - - -cdef void initialise_knm_workspaces(MirrorWorkspace ws) nogil: - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p1o = ws.sim.trace[ws.P1o_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - NodeBeamParam q_p2o = ws.sim.trace[ws.P2o_id] - - # Beam parameters for matched propagation of inputs - complex_t qx_p1o_matched_refl, qy_p1o_matched_refl - complex_t qx_p2o_matched_refl, qy_p2o_matched_refl - complex_t qx_p1o_matched_trns, qy_p1o_matched_trns - complex_t qx_p2o_matched_trns, qy_p2o_matched_trns - - double xbeta = ws.mv.xbeta - double ybeta = ws.mv.ybeta - - # Geometric factors for mirror beta angles giving - # gamma angle for each coupling - double beta_factor_11 = 2.0 - double beta_factor_22 = 2.0 - double beta_factor_12 = -(1 - ws.nr1 / ws.nr2) - double beta_factor_21 = 1 - ws.nr2 / ws.nr1 - - double lambda0 = ws.sim.model_data.lambda0 - int maxtem = ws.sim.model_data.maxtem - - # Reflection at port p1 -> p1 - qx_p1o_matched_refl = transform_q(ws.abcd_p1p1_x, q_p1i.qx, ws.nr1, ws.nr1) - knm_ws_init( - &ws.K11_ws_x, qx_p1o_matched_refl, q_p1o.qx, xbeta, beta_factor_11, ws.nr1, lambda0, maxtem, - ) - qy_p1o_matched_refl = transform_q(ws.abcd_p1p1_y, q_p1i.qy, ws.nr1, ws.nr1) - knm_ws_init( - &ws.K11_ws_y, qy_p1o_matched_refl, q_p1o.qy, ybeta, beta_factor_11, ws.nr1, lambda0, maxtem, - ) - - # Reflection at port p2 -> p2 - qx_p2o_matched_refl = transform_q(ws.abcd_p2p2_x, q_p2i.qx, ws.nr2, ws.nr2) - knm_ws_init( - &ws.K22_ws_x, qx_p2o_matched_refl, q_p2o.qx, xbeta, beta_factor_22, ws.nr2, lambda0, maxtem, - ) - qy_p2o_matched_refl = transform_q(ws.abcd_p2p2_y, q_p2i.qy, ws.nr2, ws.nr2) - knm_ws_init( - &ws.K22_ws_y, qy_p2o_matched_refl, q_p2o.qy, ybeta, -beta_factor_22, ws.nr2, lambda0, maxtem, - ) - - # Transmission at port p1 -> p2 - qx_p2o_matched_trns = transform_q(ws.abcd_p1p2_x, q_p1i.qx, ws.nr1, ws.nr2) - knm_ws_init( - &ws.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx, xbeta, beta_factor_12, ws.nr2, lambda0, maxtem, - ) - qy_p2o_matched_trns = transform_q(ws.abcd_p1p2_y, q_p1i.qy, ws.nr1, ws.nr2) - knm_ws_init( - &ws.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy, ybeta, beta_factor_12, ws.nr2, lambda0, maxtem, - ) - - # Transmission at port p2 -> p1 - qx_p1o_matched_trns = transform_q(ws.abcd_p2p1_x, q_p2i.qx, ws.nr2, ws.nr1) - knm_ws_init( - &ws.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx, xbeta, -beta_factor_21, ws.nr1, lambda0, maxtem, - ) - qy_p1o_matched_trns = transform_q(ws.abcd_p2p1_y, q_p2i.qy, ws.nr2, ws.nr1) - knm_ws_init( - &ws.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy, ybeta, beta_factor_21, ws.nr1, lambda0, maxtem, - ) - -cdef void free_knm_workspaces(MirrorWorkspace ws) nogil: - knm_ws_free(&ws.K11_ws_x) - knm_ws_free(&ws.K11_ws_y) - knm_ws_free(&ws.K22_ws_x) - knm_ws_free(&ws.K22_ws_y) - knm_ws_free(&ws.K12_ws_x) - knm_ws_free(&ws.K12_ws_y) - knm_ws_free(&ws.K21_ws_x) - knm_ws_free(&ws.K21_ws_y) - -cdef void flag_changing_knm_workspaces(MirrorWorkspace ws): - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - - bint q_p1_or_p2_changing = q_p1i.is_changing or q_p2i.is_changing - bint xbeta_changing = ws.owner.xbeta.is_changing - bint ybeta_changing = ws.owner.ybeta.is_changing - - # Reflection at port p1 -> p1 - ws.K11_ws_x.is_mm_changing = q_p1i.is_changing - ws.K11_ws_x.is_alignment_changing = xbeta_changing - ws.K11_ws_y.is_mm_changing = q_p1i.is_changing - ws.K11_ws_y.is_alignment_changing = ybeta_changing - - if knm_ws_is_changing(&ws.K11_ws_x) or knm_ws_is_changing(&ws.K11_ws_y): - LOGGER.info("%s.K11 is changing", ws.owner.name) - - # Reflection at port p2 -> p2 - ws.K22_ws_x.is_mm_changing = q_p2i.is_changing - ws.K22_ws_x.is_alignment_changing = xbeta_changing - ws.K22_ws_y.is_mm_changing = q_p2i.is_changing - ws.K22_ws_y.is_alignment_changing = ybeta_changing - - if knm_ws_is_changing(&ws.K22_ws_x) or knm_ws_is_changing(&ws.K22_ws_y): - LOGGER.info("%s.K22 is changing", ws.owner.name) - - # Transmission at port p1 -> p2 - ws.K12_ws_x.is_mm_changing = q_p1_or_p2_changing - ws.K12_ws_x.is_alignment_changing = xbeta_changing - ws.K12_ws_y.is_mm_changing = q_p1_or_p2_changing - ws.K12_ws_y.is_alignment_changing = ybeta_changing - - if knm_ws_is_changing(&ws.K12_ws_x) or knm_ws_is_changing(&ws.K12_ws_y): - LOGGER.info("%s.K12 is changing", ws.owner.name) - - # Transmission at port p2 -> p1 - ws.K21_ws_x.is_mm_changing = q_p1_or_p2_changing - ws.K21_ws_x.is_alignment_changing = xbeta_changing - ws.K21_ws_y.is_mm_changing = q_p1_or_p2_changing - ws.K21_ws_y.is_alignment_changing = ybeta_changing - - if knm_ws_is_changing(&ws.K21_ws_x) or knm_ws_is_changing(&ws.K21_ws_y): - LOGGER.info("%s.K21 is changing", ws.owner.name) - -cdef void update_changing_knm_workspaces(MirrorWorkspace ws) nogil: - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p1o = ws.sim.trace[ws.P1o_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - NodeBeamParam q_p2o = ws.sim.trace[ws.P2o_id] - - complex_t qx_p1o_matched_refl, qy_p1o_matched_refl - complex_t qx_p2o_matched_refl, qy_p2o_matched_refl - complex_t qx_p1o_matched_trns, qy_p1o_matched_trns - complex_t qx_p2o_matched_trns, qy_p2o_matched_trns - - double xbeta = ws.mv.xbeta - double ybeta = ws.mv.ybeta - - if ws.K11_ws_x.is_mm_changing: - qx_p1o_matched_refl = transform_q(ws.abcd_p1p1_x, q_p1i.qx, ws.nr1, ws.nr1) - knm_ws_recompute(&ws.K11_ws_x, qx_p1o_matched_refl, q_p1o.qx, xbeta) - elif ws.K11_ws_x.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K11_ws_x, xbeta) - - if ws.K11_ws_y.is_mm_changing: - qy_p1o_matched_refl = transform_q(ws.abcd_p1p1_y, q_p1i.qy, ws.nr1, ws.nr1) - knm_ws_recompute(&ws.K11_ws_y, qy_p1o_matched_refl, q_p1o.qy, ybeta) - elif ws.K11_ws_y.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K11_ws_y, ybeta) - - if ws.K22_ws_x.is_mm_changing: - qx_p2o_matched_refl = transform_q(ws.abcd_p2p2_x, q_p2i.qx, ws.nr2, ws.nr2) - knm_ws_recompute(&ws.K22_ws_x, qx_p2o_matched_refl, q_p2o.qx, xbeta) - elif ws.K22_ws_x.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K22_ws_x, xbeta) - - if ws.K22_ws_y.is_mm_changing: - qy_p2o_matched_refl = transform_q(ws.abcd_p2p2_y, q_p2i.qy, ws.nr2, ws.nr2) - knm_ws_recompute(&ws.K22_ws_y, qy_p2o_matched_refl, q_p2o.qy, ybeta) - elif ws.K22_ws_y.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K22_ws_y, ybeta) - - if ws.K12_ws_x.is_mm_changing: - qx_p2o_matched_trns = transform_q(ws.abcd_p1p2_x, q_p1i.qx, ws.nr1, ws.nr2) - knm_ws_recompute(&ws.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx, xbeta) - elif ws.K12_ws_x.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K12_ws_x, xbeta) - - if ws.K12_ws_y.is_mm_changing: - qy_p2o_matched_trns = transform_q(ws.abcd_p1p2_y, q_p1i.qy, ws.nr1, ws.nr2) - knm_ws_recompute(&ws.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy, ybeta) - elif ws.K12_ws_y.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K12_ws_y, ybeta) - - if ws.K21_ws_x.is_mm_changing: - qx_p1o_matched_trns = transform_q(ws.abcd_p2p1_x, q_p2i.qx, ws.nr2, ws.nr1) - knm_ws_recompute(&ws.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx, xbeta) - elif ws.K21_ws_x.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K21_ws_x, xbeta) - - if ws.K21_ws_y.is_mm_changing: - qy_p1o_matched_trns = transform_q(ws.abcd_p2p1_y, q_p2i.qy, ws.nr2, ws.nr1) - knm_ws_recompute(&ws.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy, ybeta) - elif ws.K21_ws_y.is_alignment_changing: - knm_ws_recompute_misalignment(&ws.K21_ws_y, ybeta) - - -cdef void compute_scattering_matrices(MirrorWorkspace ws): - # Reflection p1 -> p1 - if knm_ws_is_changing(&ws.K11_ws_x) or knm_ws_is_changing(&ws.K11_ws_y): - compute_knm_matrix_bh(&ws.K11_ws_x, &ws.K11_ws_y, ws.sim.homs_view, ws.K11.data_view) - reverse_gouy_phases(ws.K11.data_view, ws.sim.homs_view, &ws.K11_ws_x, &ws.K11_ws_y, ws.K11.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K11.data_view, ws.K11.data_view) - - flip_odd_horizontal(ws.K11.data_view, ws.sim.homs_view, ws.K11.data_view) - - knm_loss(ws.K11.data_view, ws.K11_loss) - - # Reflection p2 -> p2 - if knm_ws_is_changing(&ws.K22_ws_x) or knm_ws_is_changing(&ws.K22_ws_y): - compute_knm_matrix_bh(&ws.K22_ws_x, &ws.K22_ws_y, ws.sim.homs_view, ws.K22.data_view) - reverse_gouy_phases(ws.K22.data_view, ws.sim.homs_view, &ws.K22_ws_x, &ws.K22_ws_y, ws.K22.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K22.data_view, ws.K22.data_view) - - flip_odd_horizontal(ws.K22.data_view, ws.sim.homs_view, ws.K22.data_view) - - knm_loss(ws.K22.data_view, ws.K22_loss) - - # Transmission p1 -> p2 - if knm_ws_is_changing(&ws.K12_ws_x) or knm_ws_is_changing(&ws.K12_ws_y): - compute_knm_matrix_bh(&ws.K12_ws_x, &ws.K12_ws_y, ws.sim.homs_view, ws.K12.data_view) - reverse_gouy_phases(ws.K12.data_view, ws.sim.homs_view, &ws.K12_ws_x, &ws.K12_ws_y, ws.K12.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K12.data_view, ws.K12.data_view) - - knm_loss(ws.K12.data_view, ws.K12_loss) - - # Transmission p2 -> p1 - if knm_ws_is_changing(&ws.K21_ws_x) or knm_ws_is_changing(&ws.K21_ws_y): - compute_knm_matrix_bh(&ws.K21_ws_x, &ws.K21_ws_y, ws.sim.homs_view, ws.K21.data_view) - reverse_gouy_phases(ws.K21.data_view, ws.sim.homs_view, &ws.K21_ws_x, &ws.K21_ws_y, ws.K21.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K21.data_view, ws.K21.data_view) - - knm_loss(ws.K21.data_view, ws.K21_loss) diff --git a/src/finesse/components/modal/modulator.pxd b/src/finesse/components/modal/modulator.pxd index bb02b40fa7984e925b0e5506a3eb515f4191a95a..a37f83b0988cacb1e8c70e2377c0f13f2519377c 100644 --- a/src/finesse/components/modal/modulator.pxd +++ b/src/finesse/components/modal/modulator.pxd @@ -3,12 +3,13 @@ from finesse.element cimport BaseCValues from finesse.cmatrix cimport SubCCSView, SubCCSView1DArray from finesse.knm cimport KnmMatrix, KnmWorkspace from finesse.simulations.base cimport frequency_info_t, ModelData, BaseSimulation +from finesse.components.modal.workspace cimport KnmConnectorWorkspace import numpy as np cimport numpy as np -cdef class ModulatorWorkspace(ConnectorWorkspace): +cdef class ModulatorWorkspace(KnmConnectorWorkspace): cdef public: # Complete scattering matrices for each propagation direction KnmMatrix K12 @@ -32,9 +33,8 @@ cdef class ModulatorWorkspace(ConnectorWorkspace): KnmWorkspace K21_ws_x KnmWorkspace K21_ws_y - -cdef void initialise_knm_workspaces(ModulatorWorkspace ws) nogil -cdef void free_knm_workspaces(ModulatorWorkspace ws) nogil -cdef void flag_changing_knm_workspaces(ModulatorWorkspace ws) -cdef void update_changing_knm_workspaces(ModulatorWorkspace ws) nogil -cdef void compute_scattering_matrices(ModulatorWorkspace ws) + cdef void initialise_knm_workspaces(self) nogil + cdef void free_knm_workspaces(self) nogil + cdef void flag_changing_knm_workspaces(self) + cdef void update_changing_knm_workspaces(self) nogil + cdef void compute_scattering_matrices(self) diff --git a/src/finesse/components/modal/modulator.pyx b/src/finesse/components/modal/modulator.pyx index 65901eced605028d6414cf3278de367371dba72d..5cf3923808100d58aa3d497b5243017a86b28c49 100644 --- a/src/finesse/components/modal/modulator.pyx +++ b/src/finesse/components/modal/modulator.pyx @@ -31,9 +31,12 @@ cdef extern from "constants.h": LOGGER = logging.getLogger(__name__) -cdef class ModulatorWorkspace(ConnectorWorkspace): +cdef class ModulatorWorkspace(KnmConnectorWorkspace): def __init__(self, object owner, BaseSimulation sim, bint refill): from finesse.components.workspace import Connections + + self.str_couplings = ("12", "21") + super().__init__(owner, sim, refill, Connections(sim)) self.P1i_id = sim.node_id(owner.p1.i) @@ -41,6 +44,110 @@ cdef class ModulatorWorkspace(ConnectorWorkspace): self.P2i_id = sim.node_id(owner.p2.i) self.P2o_id = sim.node_id(owner.p2.o) + cdef void initialise_knm_workspaces(self) nogil: + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p1o = self.sim.trace[self.P1o_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + NodeBeamParam q_p2o = self.sim.trace[self.P2o_id] + + # Beam parameters for matched propagation of inputs + complex_t qx_p1o_matched_trns, qy_p1o_matched_trns + complex_t qx_p2o_matched_trns, qy_p2o_matched_trns + + double lambda0 = self.sim.model_data.lambda0 + int maxtem = self.sim.model_data.maxtem + + # Transmission at port p1 -> p2 + qx_p2o_matched_trns = q_p1i.qx + knm_ws_init( + &self.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx, 0, 0, self.nr2, lambda0, maxtem, + ) + qy_p2o_matched_trns = q_p1i.qy + knm_ws_init( + &self.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy, 0, 0, self.nr2, lambda0, maxtem, + ) + + # Transmission at port p2 -> p1 + qx_p1o_matched_trns = q_p2i.qx + knm_ws_init( + &self.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx, 0, 0, self.nr1, lambda0, maxtem, + ) + qy_p1o_matched_trns = q_p2i.qy + knm_ws_init( + &self.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy, 0, 0, self.nr1, lambda0, maxtem, + ) + + cdef void free_knm_workspaces(self) nogil: + knm_ws_free(&self.K12_ws_x) + knm_ws_free(&self.K12_ws_y) + knm_ws_free(&self.K21_ws_x) + knm_ws_free(&self.K21_ws_y) + + cdef void flag_changing_knm_workspaces(self): + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + + bint is_changing = q_p1i.is_changing or q_p2i.is_changing + + # Transmission at port p1 -> p2 + self.K12_ws_x.is_mm_changing = is_changing + self.K12_ws_y.is_mm_changing = is_changing + + if knm_ws_is_changing(&self.K12_ws_x) or knm_ws_is_changing(&self.K12_ws_y): + LOGGER.info("%s.K12 is changing", self.owner.name) + + # Transmission at port p2 -> p1 + self.K21_ws_x.is_mm_changing = is_changing + self.K21_ws_y.is_mm_changing = is_changing + + if knm_ws_is_changing(&self.K21_ws_x) or knm_ws_is_changing(&self.K21_ws_y): + LOGGER.info("%s.K21 is changing", self.owner.name) + + cdef void update_changing_knm_workspaces(self) nogil: + cdef: + NodeBeamParam q_p1i = self.sim.trace[self.P1i_id] + NodeBeamParam q_p1o = self.sim.trace[self.P1o_id] + NodeBeamParam q_p2i = self.sim.trace[self.P2i_id] + NodeBeamParam q_p2o = self.sim.trace[self.P2o_id] + + complex_t qx_p1o_matched_trns, qy_p1o_matched_trns + complex_t qx_p2o_matched_trns, qy_p2o_matched_trns + + if self.K12_ws_x.is_mm_changing: + qx_p2o_matched_trns = q_p1i.qx + knm_ws_recompute_mismatch(&self.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx) + + if self.K12_ws_y.is_mm_changing: + qy_p2o_matched_trns = q_p1i.qy + knm_ws_recompute_mismatch(&self.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy) + + if self.K21_ws_x.is_mm_changing: + qx_p1o_matched_trns = q_p2i.qx + knm_ws_recompute_mismatch(&self.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx) + + if self.K21_ws_y.is_mm_changing: + qy_p1o_matched_trns = q_p2i.qy + knm_ws_recompute_mismatch(&self.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy) + + cdef void compute_scattering_matrices(self): + # Transmission p1 -> p2 + if knm_ws_is_changing(&self.K12_ws_x) or knm_ws_is_changing(&self.K12_ws_y): + compute_knm_matrix_bh(&self.K12_ws_x, &self.K12_ws_y, self.sim.homs_view, self.K12.data_view) + reverse_gouy_phases(self.K12.data_view, self.sim.homs_view, &self.K12_ws_x, &self.K12_ws_y, self.K12.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K12.data_view, self.K12.data_view) + + # Transmission p2 -> p1 + if knm_ws_is_changing(&self.K21_ws_x) or knm_ws_is_changing(&self.K21_ws_y): + compute_knm_matrix_bh(&self.K21_ws_x, &self.K21_ws_y, self.sim.homs_view, self.K21.data_view) + reverse_gouy_phases(self.K21.data_view, self.sim.homs_view, &self.K21_ws_x, &self.K21_ws_y, self.K21.data_view) + + if self.sim.model_data.zero_K00: + zero_tem00_phase(self.K21.data_view, self.K21.data_view) + def modulator_fill(object modulator, ModulatorWorkspace ws, dict values, dict cporders, unicode mod_type): cdef: @@ -92,109 +199,3 @@ def modulator_fill(object modulator, ModulatorWorkspace ws, dict values, dict cp mat[:] = factors_12[_order][:] with sim.component_edge_fill3(ws.owner_id, ws.connections.P2i_P1o_idx, f1.index, f2.index) as mat: mat[:] = factors_21[_order][:] - - -cdef void initialise_knm_workspaces(ModulatorWorkspace ws) nogil: - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p1o = ws.sim.trace[ws.P1o_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - NodeBeamParam q_p2o = ws.sim.trace[ws.P2o_id] - - # Beam parameters for matched propagation of inputs - complex_t qx_p1o_matched_trns, qy_p1o_matched_trns - complex_t qx_p2o_matched_trns, qy_p2o_matched_trns - - double lambda0 = ws.sim.model_data.lambda0 - int maxtem = ws.sim.model_data.maxtem - - # Transmission at port p1 -> p2 - qx_p2o_matched_trns = q_p1i.qx - knm_ws_init( - &ws.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx, 0, 0, ws.nr2, lambda0, maxtem, - ) - qy_p2o_matched_trns = q_p1i.qy - knm_ws_init( - &ws.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy, 0, 0, ws.nr2, lambda0, maxtem, - ) - - # Transmission at port p2 -> p1 - qx_p1o_matched_trns = q_p2i.qx - knm_ws_init( - &ws.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx, 0, 0, ws.nr1, lambda0, maxtem, - ) - qy_p1o_matched_trns = q_p2i.qy - knm_ws_init( - &ws.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy, 0, 0, ws.nr1, lambda0, maxtem, - ) - -cdef void free_knm_workspaces(ModulatorWorkspace ws) nogil: - knm_ws_free(&ws.K12_ws_x) - knm_ws_free(&ws.K12_ws_y) - knm_ws_free(&ws.K21_ws_x) - knm_ws_free(&ws.K21_ws_y) - -cdef void flag_changing_knm_workspaces(ModulatorWorkspace ws): - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - - bint is_changing = q_p1i.is_changing or q_p2i.is_changing - - # Transmission at port p1 -> p2 - ws.K12_ws_x.is_mm_changing = is_changing - ws.K12_ws_y.is_mm_changing = is_changing - - if knm_ws_is_changing(&ws.K12_ws_x) or knm_ws_is_changing(&ws.K12_ws_y): - LOGGER.info("%s.K12 is changing", ws.owner.name) - - # Transmission at port p2 -> p1 - ws.K21_ws_x.is_mm_changing = is_changing - ws.K21_ws_y.is_mm_changing = is_changing - - if knm_ws_is_changing(&ws.K21_ws_x) or knm_ws_is_changing(&ws.K21_ws_y): - LOGGER.info("%s.K21 is changing", ws.owner.name) - -cdef void update_changing_knm_workspaces(ModulatorWorkspace ws) nogil: - cdef: - NodeBeamParam q_p1i = ws.sim.trace[ws.P1i_id] - NodeBeamParam q_p1o = ws.sim.trace[ws.P1o_id] - NodeBeamParam q_p2i = ws.sim.trace[ws.P2i_id] - NodeBeamParam q_p2o = ws.sim.trace[ws.P2o_id] - - complex_t qx_p1o_matched_trns, qy_p1o_matched_trns - complex_t qx_p2o_matched_trns, qy_p2o_matched_trns - - if ws.K12_ws_x.is_mm_changing: - qx_p2o_matched_trns = q_p1i.qx - knm_ws_recompute_mismatch(&ws.K12_ws_x, qx_p2o_matched_trns, q_p2o.qx) - - if ws.K12_ws_y.is_mm_changing: - qy_p2o_matched_trns = q_p1i.qy - knm_ws_recompute_mismatch(&ws.K12_ws_y, qy_p2o_matched_trns, q_p2o.qy) - - if ws.K21_ws_x.is_mm_changing: - qx_p1o_matched_trns = q_p2i.qx - knm_ws_recompute_mismatch(&ws.K21_ws_x, qx_p1o_matched_trns, q_p1o.qx) - - if ws.K21_ws_y.is_mm_changing: - qy_p1o_matched_trns = q_p2i.qy - knm_ws_recompute_mismatch(&ws.K21_ws_y, qy_p1o_matched_trns, q_p1o.qy) - - -cdef void compute_scattering_matrices(ModulatorWorkspace ws): - # Transmission p1 -> p2 - if knm_ws_is_changing(&ws.K12_ws_x) or knm_ws_is_changing(&ws.K12_ws_y): - compute_knm_matrix_bh(&ws.K12_ws_x, &ws.K12_ws_y, ws.sim.homs_view, ws.K12.data_view) - reverse_gouy_phases(ws.K12.data_view, ws.sim.homs_view, &ws.K12_ws_x, &ws.K12_ws_y, ws.K12.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K12.data_view, ws.K12.data_view) - - # Transmission p2 -> p1 - if knm_ws_is_changing(&ws.K21_ws_x) or knm_ws_is_changing(&ws.K21_ws_y): - compute_knm_matrix_bh(&ws.K21_ws_x, &ws.K21_ws_y, ws.sim.homs_view, ws.K21.data_view) - reverse_gouy_phases(ws.K21.data_view, ws.sim.homs_view, &ws.K21_ws_x, &ws.K21_ws_y, ws.K21.data_view) - - if ws.sim.model_data.zero_K00: - zero_tem00_phase(ws.K21.data_view, ws.K21.data_view) diff --git a/src/finesse/components/modal/space.pyx b/src/finesse/components/modal/space.pyx index 811f97341bfc7d97d53189a090bedf8150735b69..90c47c7a924d18e10179e90be8f861968e44c87a 100644 --- a/src/finesse/components/modal/space.pyx +++ b/src/finesse/components/modal/space.pyx @@ -75,7 +75,7 @@ cdef class SpaceWorkspace(ConnectorWorkspace): ch_sym = M_sym[0][1].eval(keep_changing_symbols=True) if isinstance(ch_sym, Symbol): self.sym_abcd_B = cy_expr_new() - cy_expr_init(self.sym_abcd_B, ch_sym, self.sim) + cy_expr_init(self.sym_abcd_B, ch_sym) @cython.boundscheck(False) @cython.wraparound(False) diff --git a/src/finesse/components/modal/workspace.pxd b/src/finesse/components/modal/workspace.pxd new file mode 100644 index 0000000000000000000000000000000000000000..473596c525932efbf05cc26d1b7be847bde4641c --- /dev/null +++ b/src/finesse/components/modal/workspace.pxd @@ -0,0 +1,13 @@ +from finesse.components.workspace cimport ConnectorWorkspace + +cdef class KnmConnectorWorkspace(ConnectorWorkspace): + # Sequence of string reprs of the optical couplings + # used by initialise_knm_matrices for convenience + cdef readonly tuple str_couplings + + cdef int initialise_knm_matrices(self) except -1 + cdef void initialise_knm_workspaces(KnmConnectorWorkspace ws) nogil + cdef void free_knm_workspaces(KnmConnectorWorkspace ws) nogil + cdef void flag_changing_knm_workspaces(KnmConnectorWorkspace ws) + cdef void update_changing_knm_workspaces(KnmConnectorWorkspace ws) nogil + cdef void compute_scattering_matrices(KnmConnectorWorkspace ws) diff --git a/src/finesse/components/modal/workspace.pyx b/src/finesse/components/modal/workspace.pyx new file mode 100644 index 0000000000000000000000000000000000000000..a1ebb2c136e3afcbe047198b2c5052fcb78603a6 --- /dev/null +++ b/src/finesse/components/modal/workspace.pyx @@ -0,0 +1,66 @@ +from finesse.knm cimport KnmMatrix +from finesse.components.workspace cimport ConnectorWorkspace + +import numpy as np + +import finesse.components as components + + +cdef class KnmConnectorWorkspace(ConnectorWorkspace): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + self.initialise_knm_matrices() + + cdef int initialise_knm_matrices(self) except -1: + # K_loss arrays should also be allocated for Surfaces + # -> these give scattering losses from each (n,m) -> all (n',m') + cdef bint make_knm_losses = isinstance(self.owner, components.Surface) + + # Only allocate new memory for scattering matrices for the carrier simulation... + if not self.sim.is_audio: + for coupling in self.str_couplings: + kmat = KnmMatrix(self.sim.homs_view, self.owner.name, coupling) + setattr(self, f"K{coupling}", kmat) + + if make_knm_losses: + setattr(self, f"K{coupling}_loss", np.ones(self.sim.nhoms)) + + # ... audio simulation will use the carrier Knm matrices. + else: + DC = self.sim.DC + DC_workspaces = list(filter(lambda x: x.owner is self.owner, DC.workspaces)) + if not DC_workspaces: + raise RuntimeError( + "Bug encountered! Could not find a DC simulation workspace " + f"associated with component of name {self.owner.name}" + ) + if len(DC_workspaces) > 1: + raise RuntimeError( + "Bug encountered! Found multiple workspaces in the DC simulation " + f"associated with component of name {self.owner.name}" + ) + + DC_ws = DC_workspaces[0] + + for coupling in self.str_couplings: + setattr(self, f"K{coupling}", getattr(DC_ws, f"K{coupling}")) + if make_knm_losses: + setattr(self, f"K{coupling}_loss", getattr(DC_ws, f"K{coupling}_loss")) + + return 0 + + cdef void initialise_knm_workspaces(KnmConnectorWorkspace self) nogil: + pass + + cdef void free_knm_workspaces(KnmConnectorWorkspace self) nogil: + pass + + cdef void flag_changing_knm_workspaces(KnmConnectorWorkspace self): + pass + + cdef void update_changing_knm_workspaces(KnmConnectorWorkspace self) nogil: + pass + + cdef void compute_scattering_matrices(KnmConnectorWorkspace self): + pass diff --git a/src/finesse/components/modulator.py b/src/finesse/components/modulator.py index 6bddddaa40199583236d687ae533ea6b5d580de6..e2b7bb7f6ebb35dfb6d9fe0325601aaeaededff0 100644 --- a/src/finesse/components/modulator.py +++ b/src/finesse/components/modulator.py @@ -9,7 +9,6 @@ from finesse.components.general import InteractionType, Connector, FrequencyGene from finesse.components.node import NodeType, NodeDirection from finesse.components.modal.modulator import modulator_fill, ModulatorWorkspace from finesse.parameter import model_parameter, info_parameter, Rebuild -from finesse.knm import KnmMatrix from finesse.utilities import refractive_index LOGGER = logging.getLogger(__name__) @@ -160,13 +159,6 @@ class Modulator(Connector, FrequencyGenerator): # rather than having to use `in` checks all the time. ws.coupling_orders = {} - # 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) - for f1 in sim.frequencies: ws.coupling_orders[f1, f1] = 0 for f2 in sim.frequencies: @@ -215,28 +207,3 @@ class Modulator(Connector, FrequencyGenerator): # TODO ddb : move all this over to new C verions values, is_changing = self._eval_parameters() modulator_fill(self, ws, values, ws.coupling_orders, self.mod_type) - - 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 Modulator of name {self.name}" - ) - if len(DC_workspaces) > 1: - raise RuntimeError( - "Bug encountered! Found multiple workspaces in the DC simulation " - f"associated with Modulator 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}")) diff --git a/src/finesse/components/space.py b/src/finesse/components/space.py index 66e5560e7f18a46f4702a99b8b30d073b2b60acd..bb7b7e5c23909a5eb5cc2e4333970ce788f477b7 100644 --- a/src/finesse/components/space.py +++ b/src/finesse/components/space.py @@ -13,7 +13,7 @@ from finesse.components.general import ( _conn_abcd_return, ) from finesse.components.workspace import ConnectionSetting -from finesse.components.modal.space import space_fill, SpaceWorkspace, SpaceConnections +from finesse.components.modal.space import space_fill, SpaceWorkspace from finesse.components.node import NodeDirection, NodeType, Port @@ -91,7 +91,7 @@ class Space(Connector): "Cannot create an unconnected space without providing a name" ) - super().__init__(name, SpaceConnections) + super().__init__(name) self.__portA = portA self.__portB = portB @@ -197,8 +197,8 @@ class Space(Connector): return value def _check_nr(self, value): - if value <= 0: - raise ValueError("Index of refraction must be > 0") + if value < 1: + raise ValueError("Index of refraction must be >= 1") return value diff --git a/src/finesse/components/surface.py b/src/finesse/components/surface.py index 0746e183f49deda82e37316bf644eb1864159f8d..84accb5dbaa35344325507a299ae701f19b57f59 100644 --- a/src/finesse/components/surface.py +++ b/src/finesse/components/surface.py @@ -80,10 +80,9 @@ class Surface(ABC, Connector): xbeta, ybeta, mass, - signal_gain, - connections_type=None, + signal_gain ): - super().__init__(name, connections_type=connections_type) + super().__init__(name) # only in constructor allow setting of none of R, T, L # -> default to equally reflective and tranmissive with no loss diff --git a/src/finesse/components/workspace.pxd b/src/finesse/components/workspace.pxd index cfb475d1dc2a3afe450ad6e0218d92eab8b96f4d..de4b0c294e218a283107e2a84ed04bef61148f7a 100644 --- a/src/finesse/components/workspace.pxd +++ b/src/finesse/components/workspace.pxd @@ -25,11 +25,11 @@ cdef class ConnectorWorkspace(ElementWorkspace): readonly object connections FillFuncWrapper fn_c # C fill function object fn_py # Python fill function - FillFuncWrapper fn_rhs_c # C fill function - object fn_rhs_py # Python fill function + FillFuncWrapper fn_rhs_c # C RHS fill function + object fn_rhs_py # Python RHS fill function readonly bint refill # Should this matrix be refilled continuously readonly BaseSimulation sim - readonly ConnectorCallbacks callback_flag + readonly ConnectorCallbacks callback_flag # flag stating which fill methods are to be called readonly dict connection_settings cpdef compile_abcd_cy_exprs(self) diff --git a/src/finesse/components/workspace.pyx b/src/finesse/components/workspace.pyx index a797a2836754ba6f240572bf5ec18147f212cc76..c8b9ed751426d664bcecfee6cbbef257319a704c 100644 --- a/src/finesse/components/workspace.pyx +++ b/src/finesse/components/workspace.pyx @@ -1,3 +1,4 @@ + cdef class FillFuncWrapper: """ Helper class for wrapping a C fill function that diff --git a/src/finesse/cyexpr.pxd b/src/finesse/cyexpr.pxd index d35bc45cb064c4168a503da2b6ff9a559182b173..dcdf8c74935a533298a68a664b326337003c37fa 100644 --- a/src/finesse/cyexpr.pxd +++ b/src/finesse/cyexpr.pxd @@ -5,6 +5,8 @@ cdef extern from "tinyexpr.h" nogil: ctypedef struct te_variable: const char* name const void* address + int type + void* context cdef te_expr* te_compile( const char* expression, @@ -31,11 +33,7 @@ cdef struct cy_expr: cdef cy_expr* cy_expr_new() nogil -cdef int cy_expr_init( - cy_expr* ex, - object operation, - object sim, -) except -1 +cdef int cy_expr_init(cy_expr* ex, object operation) except -1 cdef void cy_expr_free(cy_expr* ex) nogil diff --git a/src/finesse/cyexpr.pyx b/src/finesse/cyexpr.pyx index 95d97c41efe24c0e8e00da043bcf52135325b226..8895bdf4bdc662ee841327ff4aca91d90bee06da 100644 --- a/src/finesse/cyexpr.pyx +++ b/src/finesse/cyexpr.pyx @@ -15,7 +15,6 @@ based math parsing and evaluation engine, `tinyexpr malloc(sizeof(cy_expr)) @@ -23,11 +22,9 @@ cdef cy_expr* cy_expr_new() nogil: ce_p.variables = NULL return ce_p -cdef int cy_expr_init( - cy_expr* ex, - object operation, - object sim -) except -1: +cdef int cy_expr_init(cy_expr* ex, object operation) except -1: + """Initialise the te_expr and te_variable type fields using the + operation object (should be an instance of Operation).""" cdef str op_str = str(operation) cdef list params = operation.parameters() # get the dependent parameter-refs cdef int Nparams = len(params) @@ -36,60 +33,23 @@ cdef int cy_expr_init( if not ex.variables: raise MemoryError() - cdef list conn_workspaces = sim.workspaces - cdef list var_workspaces = sim.variable_workspaces - cdef Py_ssize_t i - cdef int pindex - cdef str param_name, pname cdef Parameter p - cdef ElementWorkspace ws for i in range(Nparams): pref = params[i] # the ParameterRef object - p = pref.parameter - # Get name into a form corresponding to tinyexpr variable grammar, rules - # are *start with lowercase letter* followed by any combination of: - # - lower case letters a through z, - # - digits 0 through 9, - # - and underscore, - param_name = pref.name - pname = param_name.replace(".", "_").lower() # i.e. replace "comp.PaRaM" with "comp_param" - # replace the corresponding parameter name in the operation string - op_str = op_str.replace(param_name, pname) + # Replace the corresponding parameter name in the operation string with + # the tinyexpr compatible name version (see ParameterRef.cyexpr_name) + op_str = op_str.replace(pref.name, pref.cyexpr_name.decode()) # Also need to replace "quantity**n" with "quantity^n" as te expects powers in this form op_str = op_str.replace("**", "^") - # NOTE (sjr) Need to get workspace corresponding to the parameter currently - # as there is no double value field in Parameter that we can take - # an address of at the moment - # TODO (sjr) *If we do end up storing a double value in Parameter* then - # all this can be replaced and second arg to te_variable - # below can then just be the address of this double field - ws = get_associated_workspace(conn_workspaces, var_workspaces, p.owner) - if ws is None: - raise ValueError( - "Bug encountered! Could not find a workspace corresponding " - f"to the connector of name {p.owner.name}." - ) - - # Get the index of the changing parameter corresponding to its name - # in the chprm field of ElementWorkspace - # NOTE (sjr) Again (as above), don't need this if we implement a double value field in Parameter - pindex = get_chprm_index(ws, p.name) - if pindex == -1: - raise ValueError( - "Bug encountered! Could not find index corresponding " - f"to the changing parameter of name {p.name}." - ) - - # Need to make a byte string to correspond to char* in C - byte_param_str = pname.encode("UTF-8") - ex.variables[i] = te_variable(byte_param_str, ws.chprm_target[pindex]) + p = pref.parameter + ex.variables[i] = te_variable(pref.cyexpr_name, &p.__cvalue, 0, NULL) cdef int err # position of parsing error in operation expression - # Again make byte str of operation expression, and store it - # in cy_expr instance for info / debugging purposes later + # Make byte str of operation expression, and store it in + # cy_expr instance for info / debugging purposes later byte_op_str = op_str.encode("UTF-8") ex.expression = byte_op_str ex.expr = te_compile(byte_op_str, ex.variables, Nparams, &err) @@ -116,37 +76,5 @@ cdef void cy_expr_free(cy_expr* ex) nogil: ex = NULL cdef double cy_expr_eval(cy_expr* ex) nogil: + """Evaluate the cythonised symbolic expression.""" return te_eval(ex.expr) - - -cdef int get_chprm_index(ElementWorkspace ws, pname): - cdef: - int i - int N = ws.num_changing_parameters - - for i in range(N): - if (ws.chprm[i]).name == pname: - return i - - return -1 - - -cdef ElementWorkspace get_associated_workspace( - list connector_workspaces, - list variable_workspaces, - object connector, -): - cdef: - Py_ssize_t Nc = len(connector_workspaces) - Py_ssize_t Nv = len(variable_workspaces) - Py_ssize_t i - - for i in range(Nc): - if connector_workspaces[i].owner is connector: - return connector_workspaces[i] - - for i in range(Nv): - if variable_workspaces[i].owner is connector: - return variable_workspaces[i] - - return None diff --git a/src/finesse/detectors/camera.py b/src/finesse/detectors/camera.py index 6a00e1c9e4ff7e3d902a9acd1770e33d9824a045..2b70b599592663d4911b148e61a4268abea672c2 100644 --- a/src/finesse/detectors/camera.py +++ b/src/finesse/detectors/camera.py @@ -30,7 +30,7 @@ from finesse.detectors.compute import ( field_camera_output, ccd_output, ) -from finesse.utilities.misc import find_nearest +from finesse.utilities.misc import find_nearest, is_iterable LOGGER = logging.getLogger(__name__) @@ -682,7 +682,7 @@ def _check_limits(value): elif isinstance(value, np.ndarray): value = [value.min(), value.max()] - elif hasattr(value, "__getitem__"): + elif is_iterable(value): if len(value) != 2: raise TypeError( "Expected limit to be a single number or sequence of size 2 " diff --git a/src/finesse/detectors/cavity_detector.py b/src/finesse/detectors/cavity_detector.py index 1575f3dbadb81ded588d3c34086d8c5909b4485b..1fbcc8806c1003b7638f00789e9b7fdc61bec48f 100644 --- a/src/finesse/detectors/cavity_detector.py +++ b/src/finesse/detectors/cavity_detector.py @@ -38,8 +38,11 @@ class CavityPropertyDetector(Detector): name : str Name of newly created cavity property detector. - cavity : str - Name of the cavity to probe. + cavity : str or :class:`.Cavity` + The cavity to probe. If the name is provided then the + :attr:`.CavityPropertyDetector.cavity` attribute will point + to the corresponding :class:`.Cavity` object when adding + this detector to a :class:`.Model` instance. prop : str or :class:`.CavityProperty` Property of the cavity to probe. diff --git a/src/finesse/detectors/general.py b/src/finesse/detectors/general.py index f94e5e3d0ebe07dee76a7db3ce5f68d3e6d4040d..f905bde22608f3863f65c1b9e0965ded01cad847 100644 --- a/src/finesse/detectors/general.py +++ b/src/finesse/detectors/general.py @@ -131,10 +131,10 @@ class Detector(ABC, ModelElement): Identifier for the mode indices to generate. This can be: - An iterable of mode indices, where each element in the iterable must unpack to two - integer convertible values. + integer convertible values. - A string identifying the type of modes to include, must be one of "even", "odd", - "tangential" (or "x") or "sagittal" (or "y"). + "tangential" (or "x") or "sagittal" (or "y"). maxtem : int, optional; default: None Optional maximum mode order, applicable only for when `select` is a string. This is diff --git a/src/finesse/element.pxd b/src/finesse/element.pxd index b316ce28ad51c4bd1c54efa8767146188bd148d5..5f0a82ef54d6c112181f8bfb0d4890a3744bcf5f 100644 --- a/src/finesse/element.pxd +++ b/src/finesse/element.pxd @@ -33,5 +33,5 @@ cdef class ElementWorkspace: cy_expr** chprm_expr # Compiled changing symbolic expressions array bint first - cpdef int compile_cy_exprs(self, object sim) except -1 + cpdef int compile_cy_exprs(self) except -1 cpdef update_parameter_values(self) diff --git a/src/finesse/element.pyx b/src/finesse/element.pyx index 1f0b960cf1ef3c571085b6aee78223947d6bb1bc..64498128c73463c69ea538d36389e332b46cb12e 100644 --- a/src/finesse/element.pyx +++ b/src/finesse/element.pyx @@ -437,8 +437,9 @@ cdef class ElementWorkspace: symbolic_params.append(p) else: numeric_params.append(p) - else: - numeric_params.append(p) + # ddb : I don't see why a parameter that isn't changing should be added here? + # else: + # numeric_params.append(p) params = numeric_params + symbolic_params if self.num_changing_parameters > 0: @@ -450,19 +451,22 @@ cdef class ElementWorkspace: if p.is_changing: if p.state == ParameterState.Symbolic: self.chprm_expr[i] = cy_expr_new() - + self.chprm[i] = p Py_XINCREF(self.chprm[i]) (self.values).get_param_ptr(p.name, &self.chprm_target[i]) i += 1 - cpdef int compile_cy_exprs(self, object sim) except -1: + cpdef int compile_cy_exprs(self) except -1: if not self.num_changing_parameters: return 0 for i in range(self.num_changing_parameters): if (self.chprm[i]).state == ParameterState.Symbolic: - if cy_expr_init(self.chprm_expr[i], (self.chprm[i]).value, sim) == -1: + if cy_expr_init( + self.chprm_expr[i], + (self.chprm[i]).value + ) == -1: return -1 return 0 @@ -491,9 +495,9 @@ cdef class ElementWorkspace: for i in range(self.num_changing_parameters): if self.chprm_target[i] == NULL or self.chprm[i] == NULL: raise MemoryError() - + if (self.chprm[i]).state == ParameterState.Numeric: - self.chprm_target[i][0] = (self.chprm[i]).__value + self.chprm_target[i][0] = (self.chprm[i]).__cvalue elif (self.chprm[i]).state == ParameterState.NONE: self.chprm_target[i][0] = 0 # maybe should be NaN? elif (self.chprm[i]).state == ParameterState.Symbolic: @@ -501,7 +505,6 @@ cdef class ElementWorkspace: else: raise ValueError("Parameter state was unexpected") - def __dealloc__(self): if self.num_changing_parameters > 0: for i in range(self.num_changing_parameters): diff --git a/src/finesse/knm.pxd b/src/finesse/knm.pxd index 0693fe27668c392500dd011599a1bc745e264d8a..436eec05e952df308c1c58648a4a9529f2b328f7 100644 --- a/src/finesse/knm.pxd +++ b/src/finesse/knm.pxd @@ -8,16 +8,16 @@ cdef extern from "constants.h": cdef class KnmMatrix: + cdef readonly: + np.ndarray data + str name + cdef: - readonly np.ndarray data complex_t[:, ::1] data_view complex_t* ptr int stride1 int stride2 - int[:, ::1] modes_view - - object comp - unicode kdir + const int[:, ::1] modes_view cdef (Py_ssize_t, Py_ssize_t) field_indices_from(self, key) cdef complex_t coupling(self, int n1, int m1, int n2, int m2) nogil diff --git a/src/finesse/knm.pyx b/src/finesse/knm.pyx index c933868a988839da8012ec118f075581a872f473..fa03239f30a738d1ca1f91aae0e7178888529192 100644 --- a/src/finesse/knm.pyx +++ b/src/finesse/knm.pyx @@ -21,55 +21,99 @@ from finesse.cymath cimport complex_t from finesse.cymath.complex cimport carg, conj, creal, cimag, csqrt # Standard complex.h functions from finesse.cymath.complex cimport ceq, inverse_unsafe, crotate, crotate2, cpow_re, cnorm from finesse.cymath.math cimport exp, sqrt, cos, sin # Standard math.h functions -from finesse.cymath.math cimport factorial, msign, nmin +from finesse.cymath.math cimport factorial, msign, nmin, nmax from finesse.cymath.gaussbeam cimport gouy, waist_size, z_q, z_R from finesse.cymath.homs cimport field_index +import copy + cdef class KnmMatrix: """ - Higher-Order-Mode ( scattering matrix container. Essentially a wrapper + Higher-Order-Mode scattering matrix container. Essentially a wrapper around a 2D NumPy array with methods for conveniently accessing specific - couplings and plotting the matrix as a colormesh. + couplings and plotting the matrix as a color-mesh. + + The underlying :class:`numpy.ndarray` object can be accessed via the ``data`` + attribute. - All component coupling coefficient matrices are stored as objects of this - type - e.g. the scattering matrix on reflection from the first surface of some - mirror ``M`` that is part of a model ``ifo`` can be accessed with ``ifo.M.K11``. - The `data` attribute of this object then allows direct access to the underlying - NumPy array, whilst the `modes` attribute gives access to a copy of the mode - indices used to produce this scattering matrix (i.e. the - :attr:`finesse.model.Model.homs` array of the associated model). + Construction of KnmMatrix objects should, generally, not be performed manually. Objects + of this type are the return type of the general scattering matrix computing routines, + see :func:`.make_scatter_matrix`. To make a KnmMatrix object from a pre-existing + matrix of complex coupling coefficients, use :meth:`.KnmMatrix.from_buffer`. Parameters ---------- - matrix : :class:`numpy.ndarray` - A 2D array of complex values corresponding to the underlying - buffer of scattering matrix data. + modes : array-like + A 2D array, or memory-view to the array, of the mode indices + associated with the scattering matrix. - component : :class:`.Connector` - The component that the scattering matrix belongs to. + comp_name : str, optional; default: "" + Name of the component associated with the matrix (if any). - kdir : str + kdir : str, optional; default: "" A string representing the coupling direction of the matrix, e.g. - "11" for a reflection at the first surface of an optic. + "11" for a reflection at the first surface of a mirror-like optic. """ - def __init__( - self, - np.ndarray[complex_t, ndim=2, mode='c'] matrix, - object component, - unicode kdir="", - ): - self.data = matrix - self.comp = component - self.kdir = kdir + def __init__(self, const int[:, ::1] modes, comp_name="", kdir=""): + if comp_name: + self.name = f"{comp_name}.K{kdir}" + else: + self.name = f"K{kdir}" + self.modes_view = modes + matrix = np.eye(self.modes_view.shape[0], dtype=np.complex128) + self.data = matrix self.data_view = self.data - self.modes_view = component._model.homs self.ptr = &self.data_view[0,0] self.stride1 = matrix[:].strides[0]//sizeof(complex_t) self.stride2 = matrix[:].strides[1]//sizeof(complex_t) + @staticmethod + def from_buffer( + np.ndarray[complex_t, ndim=2, mode='c'] buffer, + const int[:, ::1] modes, + **kwargs, + ): + """Construct a KnmMatrix object from a pre-existing 2D array buffer. + + This method is useful for creating KnmMatrix objects from scattering matrix + data which have been computed externally (i.e. not via Finesse). + + Parameters + ---------- + buffer : :class:`numpy.ndarray` + A 2D array of complex values corresponding to the underlying + buffer of scattering matrix data. + + modes : array-like + A 2D array, or memory-view to the array, of the mode indices + associated with the scattering matrix. + + kwargs : keyword arguments, optional + Additional args to pass to constructor of KnmMatrix. + + Returns + ------- + kmat : :class:`.KnmMatrix` + The KnmMatrix wrapper around the array buffer. + """ + cdef Py_ssize_t N = modes.shape[0] + if buffer[:].shape[0] != N or buffer[:].shape[1] != N: + bshape = (buffer[:].shape[0], buffer[:].shape[1]) + raise ValueError( + f"Error in KnmMatrix.from_buffer:\n" + f" Shape of matrix buffer {bshape} inconsistent " + f"with number of modes ({N})" + ) + + cdef KnmMatrix kmat = KnmMatrix(modes, **kwargs) + kmat.data = buffer + kmat.data_view = buffer + + return kmat + def __str__(self): cdef: Py_ssize_t i, j @@ -77,8 +121,6 @@ cdef class KnmMatrix: Py_ssize_t N = self.modes_view.shape[0] s = "" - comp_name = self.comp.name - for i in range(N): n1 = self.modes_view[i][0] m1 = self.modes_view[i][1] @@ -86,27 +128,70 @@ cdef class KnmMatrix: n2 = self.modes_view[j][0] m2 = self.modes_view[j][1] - s += (f"{comp_name}.K{self.kdir}[{i}][{j}]: coupling {n1} {m1} " + s += (f"{self.name}[{i}][{j}]: coupling {n1} {m1} " f"-> {n2} {m2} = {self.data_view[i][j]}\n") return s cdef (Py_ssize_t, Py_ssize_t) field_indices_from(self, key): cdef: - Py_ssize_t i, j + Py_ssize_t i = 0 + Py_ssize_t j = 0 int n1, m1, n2, m2 Py_ssize_t N = self.modes_view.shape[0] - Py_ssize_t key_size = len(key) - if key_size == 2: # getting [i, j] element directly - i, j = key - elif key_size == 4: # getting nm -> n'm' coupling - n1, m1, n2, m2 = key + Py_ssize_t key_size + str skey + str mode1, mode2 + list smodes + + if isinstance(key, str): + skey = "".join(key.split()) # strip all whitespace from the key + key_size = len(skey) + smodes = skey.split("->") + if not smodes or len(smodes) > 2: + i = j = N + else: + # Key is a str of format "nmn'm'" + if len(smodes) == 1: + if key_size != 4: + i = j = N + else: + mode1, mode2 = smodes[0][:2], smodes[0][2:] + + # Key is a str of format "nm->n'm'" + elif len(smodes) == 2: + if key_size != 6: + i = j = N + else: + mode1, mode2 = smodes + + if not i and not j: + # Have to convert to python int class + # first before setting C ints + nn1 = int(mode1[0]) + mm1 = int(mode1[1]) + nn2 = int(mode2[0]) + mm2 = int(mode2[1]) + n1 = nn1 + m1 = mm1 + n2 = nn2 + m2 = mm2 + + i = field_index(n1, m1, self.modes_view) + j = field_index(n2, m2, self.modes_view) - i = field_index(n1, m1, self.modes_view) - j = field_index(n2, m2, self.modes_view) else: - i = j = N + key_size = len(key) + if key_size == 2: # getting [i, j] element directly + i, j = key + elif key_size == 4: # getting nm -> n'm' coupling + n1, m1, n2, m2 = key + + i = field_index(n1, m1, self.modes_view) + j = field_index(n2, m2, self.modes_view) + else: + i = j = N return i, j @@ -118,14 +203,7 @@ cdef class KnmMatrix: i, j = self.field_indices_from(key) if i >= N or j >= N: - msg = """Key must be an integer sequence of: - - - length 2 -- i.e. an i, j pair corresponding to the field indices; where both i and j - must be less than the number of modes present. - - or length 4 -- i.e. a n1, m1, n2, m2 coupling; where each mode index must be present - in the mode indices of the model. - """ - raise ValueError(msg) + raise ValueError(f"Invalid or non-existent key {key}") return self.data_view[j][i] @@ -135,8 +213,15 @@ cdef class KnmMatrix: Parameters ---------- - key : integer sequence - An integer sequence of: + key : str or integer sequence + A string of the format: + + - "nm->n'm'" or, + - "nmn'm'", + + where n, m, n' and m' are all convertible to integers + + Or an integer sequence of: - length 2 -- i.e. an i, j pair corresponding to the field indices. - or length 4 -- i.e. a n1, m1, n2, m2 coupling. @@ -167,13 +252,8 @@ cdef class KnmMatrix: return self.data_view[j][i] - @property - def component(self): - """The component that the scattering matrix belongs to.""" - return self.comp - - def plot(self, mode="amplitude", cmap=None, show=True, filename=None): - """Plots the coupling coefficient matrix. + def plot(self, mode="amplitude", log=False, cmap=None, show=True, filename=None): + """Plots the coupling coefficient matrix as a color-mesh. Parameters ---------- @@ -199,6 +279,7 @@ cdef class KnmMatrix: A handle to the figure. """ import matplotlib.pyplot as plt + import matplotlib.colors as colors from finesse.plotting.tools import add_colorbar knm = self.data[:] @@ -223,14 +304,28 @@ cdef class KnmMatrix: raise ValueError(f"Invalid mode={mode} argument. Must be one of " "'amplitude', 'phase', 'real' or 'imag'.") - cax = ax.pcolormesh(func(knm), cmap=cmap) + if log: + if mode != "amplitude": + raise ValueError( + f"Logarithmic scaling only supported for amplitude mode." + ) + + norm = colors.LogNorm() + # Make sure any zero couplings get displayed as black rather than white + cmap_handle = copy.copy(plt.cm.get_cmap(cmap)) + cmap_handle.set_bad((0,0,0)) + else: + norm = None + cmap_handle = cmap + + cax = ax.pcolormesh(func(knm), cmap=cmap_handle, norm=norm) add_colorbar(cax, label=label) numrows, numcols = knm.shape hom_ticks = [] for i in range(self.modes_view.shape[0]): - hom_ticks.append(f"[{self.modes_view[i][0]}, {self.modes_view[i][1]}]") + hom_ticks.append(f"{self.modes_view[i][0]}{self.modes_view[i][1]}") A = np.arange(1, self.modes_view.shape[0] + 1) - 0.5 ax.set_xticks(A) @@ -248,7 +343,7 @@ cdef class KnmMatrix: if col >= 0 and col < numcols and row >= 0 and row < numrows: z = knm[row, col] return ( - f"x={hom_ticks[col]}, y={hom_ticks[row]}, z={z:.4g}" + f"nm={hom_ticks[col]}, n'm'={hom_ticks[row]}, z={z:.4g}" f" -> {label}(z) = {func(z):.4g}" ) @@ -606,11 +701,16 @@ cdef void _bayer_helms_outer_loop( ) +##### C functions for use by connector workspaces in ##### +##### computing and modifying scattering matrices ##### + + cdef void compute_knm_matrix_bh( KnmWorkspace* kws_x, KnmWorkspace* kws_y, const int[:, ::1] homs, complex_t[:, ::1] out ): + """Compute a Bayer-Helms scattering matrix.""" cdef: Py_ssize_t N = homs.shape[0] Py_ssize_t i @@ -630,19 +730,13 @@ cdef void zero_tem00_phase( const complex_t[:, ::1] knm_mat, complex_t[:, ::1] out ) nogil: - """ - Rotates all coupling coefficients in the matrix `knm_mat` by the phase of - the :math:`k_{0000}` coupling coefficient. + """Rotates all coupling coefficients in the matrix `knm_mat` + by the phase of the zeroth (00 -> 00) coupling coefficient. - Parameters - ---------- - knm_mat : contiguous array - The matrix of coupling coefficients. - - out : contiguous array, optional - Optional matrix in which store the computed results. If not specified - then a new :class:`numpy.ndarray` is returned. - """ + This has the effect of zeroing the phase of the 00 -> 00 coupling + coefficient, resulting in cavities being resonant for the 00 mode + by default if this function (along with zeroing of space Gouy phase + for 00 mode) is applied.""" cdef: Py_ssize_t i, j Py_ssize_t N = knm_mat.shape[0] @@ -660,22 +754,7 @@ cdef void flip_odd_horizontal( const int[:, ::1] homs, complex_t[:, ::1] out ) nogil: - """ - Flips the sign of all odd couplings in the sagittal plane. - - Parameters - ---------- - knm_mat : contiguous array - The matrix of coupling coefficients. - - homs : contiguous array - Array of HG modes to compute couplings between, formatted as - `[(n0, m0), (n1, m1), ...]`. - - out : contiguous array, optional - Optional matrix in which store the computed results. If not specified - then a new :class:`numpy.ndarray` is returned. - """ + """Flips the sign of all odd couplings in the sagittal plane.""" cdef: Py_ssize_t i, j int n2 @@ -703,31 +782,6 @@ cdef void reverse_gouy_phases( is added explicitly to the amplitude coefficients in a :class:`.Space` whereas the coupling coefficients are derived using a formula in which the Gouy phase resides in the equation for the spatial profile. - - Parameters - ---------- - knm_mat : contiguous array - The matrix of coupling coefficients. - - homs : contiguous array - Array of HG modes that the couplings were computed between, formatted as - `[(n0, m0), (n1, m1), ...]`. - - qx1 : np.complex128 - Input beam parameter in tangential plane. - - qy1 : np.complex128 - Input beam parameter in sagittal plane. - - qx2 : np.complex128 - Output beam parameter in tangential plane. - - qy2 : np.complex128 - Output beam parameter in sagittal plane. - - out : contiguous array, optional - Optional matrix in which store the computed results. If not specified - then a new :class:`numpy.ndarray` is returned. """ cdef: Py_ssize_t i, j @@ -752,27 +806,205 @@ cdef void knm_loss( const complex_t[:, ::1] knm_mat, double[::1] out ) nogil: - """ - Computes total losses from each mode (n, m) to all - other mode couplings (including itself). + """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]) + + +##### Python functions for computing generic scattering matrices ##### +##### not associated with any component or workspace ##### + +def compute_bayerhelms_coeff( + qx1, qx2, qy1, qy2, + double xgamma, double ygamma, + int n1, int m1, int n2, int m2, + double nr=1.0, double wavelength=1064e-9, +): + """Computes the coupling coefficient, via Bayer-Helms :cite:`BayerHelms`, + from mode (n1, m1) -> mode (n2, m2). .. note:: - If `out` is specified then this function expects it - to be initialised with all ones - e.g. via `np.ones`. + + Use :func:`.make_bayerhelms_matrix` to compute a scattering matrix of these + coefficients rather than calling this function in a loop. Parameters ---------- - knm_mat : contiguous array - The matrix of coupling coefficients. + qx1 : :class:`.BeamParam` or complex + Input beam parameter in the tangential plane. + + qx2 : :class:`.BeamParam` or complex + Output beam parameter in the tangential plane. + + qy1 : :class:`.BeamParam` or complex + Input beam parameter in the sagittal plane. + + qy2 : :class:`.BeamParam` or complex + Output beam parameter in the sagittal plane. + + xgamma : float + Misalignment angle in the tangential plane (in radians). - out : contiguous array, optional - Optional array in which store the computed results. If not specified - then a new :class:`numpy.ndarray` is returned. + ygamma : float + Misalignment angle in the sagittal plane (in radians). + + n1 : int + Input tangential mode index. + + m1 : int + Input sagittal mode index. + + n2 : int + Output tangential mode index. + + m2 : int + Output sagittal mode index. + + nr : float, optional; default: 1.0 + Refractive index of the associated medium. + + wavelength : float, optional; default: 1064 nm + Wavelength of the beam (in metres). + + Returns + ------- + coeff : complex + The complex coupling coefficient for the specified mode coupling. """ cdef: - Py_ssize_t i, j - Py_ssize_t N = knm_mat.shape[0] + KnmWorkspace kws_x + KnmWorkspace kws_y + int maxtem = nmax(n1 + m1, n2 + m2) - for i in range(N): - for j in range(N): - out[i] -= cnorm(knm_mat[j][i]) + knm_ws_init(&kws_x, qx1, qx2, xgamma, 1.0, nr, wavelength, maxtem) + knm_ws_init(&kws_y, qy1, qy2, ygamma, 1.0, nr, wavelength, maxtem) + + cdef complex_t coeff = k_nmnm(n1, m1, n2, m2, &kws_x, &kws_y) + + knm_ws_free(&kws_x) + knm_ws_free(&kws_y) + + return coeff + + +def make_bayerhelms_matrix( + qx1, qx2, qy1, qy2, + double xgamma, double ygamma, + double nr=1.0, double wavelength=1064e-9, + **mode_selection_kwargs, +): + """Constructs and computes a coupling coefficient scattering + matrix using the Bayer-Helms :cite:`BayerHelms` analytic + formalism. + + See the function :func:`.make_modes` for the arguments which + should be passed via `mode_selection_kwargs`. + + Parameters + ---------- + qx1 : :class:`.BeamParam` or complex + Input beam parameter in the tangential plane. + + qx2 : :class:`.BeamParam` or complex + Output beam parameter in the tangential plane. + + qy1 : :class:`.BeamParam` or complex + Input beam parameter in the sagittal plane. + + qy2 : :class:`.BeamParam` or complex + Output beam parameter in the sagittal plane. + + xgamma : float + Misalignment angle in the tangential plane (in radians). + + ygamma : float + Misalignment angle in the sagittal plane (in radians). + + nr : float, optional; default: 1.0 + Refractive index of the associated medium. + + wavelength : float, optional; default: 1064 nm + Wavelength of the beam (in metres). + + mode_selection_kwargs : keyword arguments + See :func:`.make_modes`. + + Returns + ------- + kmat : :class:`.KnmMatrix` + The resulting scattering matrix as a :class:`.KnmMatrix` object. + + Examples + -------- + See :ref:`arbitrary_scatter_matrices`. + """ + from finesse.utilities.homs import make_modes + + cdef: + KnmWorkspace kws_x + KnmWorkspace kws_y + + modes = make_modes(**mode_selection_kwargs) + cdef int[:, ::1] modes_view = modes + cdef int maxtem = np.max(modes[:, 0] + modes[:, 1]) + + cdef KnmMatrix kmat = KnmMatrix(modes_view) + + knm_ws_init(&kws_x, qx1, qx2, xgamma, 1.0, nr, wavelength, maxtem) + knm_ws_init(&kws_y, qy1, qy2, ygamma, 1.0, nr, wavelength, maxtem) + + compute_knm_matrix_bh(&kws_x, &kws_y, modes_view, kmat.data_view) + + knm_ws_free(&kws_x) + knm_ws_free(&kws_y) + + return kmat + + +# TODO (sjr) Create functions corresponding to different scattering +# matrix types as they are implemented in the code (i.e. +# things like aperture and map knm matrix computations) +def make_scatter_matrix(mtype, *args, **kwargs): + """Constructs and computes a coupling coefficient scattering + matrix of the specified type (via `mtype`). + + Parameters + ---------- + mtype : str + Type of scattering matrix to compute. + + The valid options are: + + - "bayerhelms" --- computes a Bayer-Helms scattering matrix, + see :func:`.make_bayerhelms_matrix`. + + *args : positional arguments + The positional arguments to pass to the relevant scatter matrix function. + + **kwargs : keyword arguments + The keyword arguments to pass to the relevant scatter matrix function. + + Returns + ------- + kmat : :class:`.KnmMatrix` + The resulting scattering matrix as a :class:`.KnmMatrix` object. + + Examples + -------- + See :ref:`arbitrary_scatter_matrices`. + """ + mtype_to_func = {"bayerhelms": make_bayerhelms_matrix} + + func = mtype_to_func.get(mtype.casefold()) + if func is None: + raise ValueError( + f"Unrecognised / not-yet-implemented mtype argument: {mtype}" + ) + + return func(*args, **kwargs) diff --git a/src/finesse/model.py b/src/finesse/model.py index 9662d0588492caae34247fd3f3fbb067b79baa23..6e5cd3ee5472784ee8538180e518cec602668235 100644 --- a/src/finesse/model.py +++ b/src/finesse/model.py @@ -1165,7 +1165,7 @@ class Model: If the matrix has already been built, the component has already been added to the model or `obj` is not of a valid type. """ - if hasattr(obj, "__getitem__"): + if is_iterable(obj): for o in obj: self.add(o) return @@ -2288,7 +2288,7 @@ class Model: return OpticalPath(fullpath) - if not hasattr(via_node, "__getitem__"): + if not is_iterable(via_node): via_node = [via_node] via_node = [via.full_name for via in via_node] @@ -2594,14 +2594,14 @@ class Model: qx_w0, qy_w0 = qx.w0, qy.w0 qx_z, qy_z = qx.z, qy.z - if not hasattr(w0_mm, "__getitem__"): + if not is_iterable(w0_mm): qx_w0 *= 1 + (w0_mm / 100) qy_w0 *= 1 + (w0_mm / 100) else: qx_w0 *= 1 + (w0_mm[0] / 100) qy_w0 *= 1 + (w0_mm[1] / 100) - if not hasattr(z_mm, "__getitem__"): + if not is_iterable(z_mm): qx_z *= 1 + (z_mm / 100) qy_z *= 1 + (z_mm / 100) else: diff --git a/src/finesse/parameter.pxd b/src/finesse/parameter.pxd index 51a8d883c0d808e208c72c0ae42ad7c1a639ce11..aa6b178467f968ebe346f6b94d282fdfc158ff5a 100644 --- a/src/finesse/parameter.pxd +++ b/src/finesse/parameter.pxd @@ -14,6 +14,7 @@ cpdef enum ParameterState: cdef class Parameter: cdef: object __value + double __cvalue str __units str __name str __full_name @@ -31,6 +32,7 @@ cdef class Parameter: cdef __cdeepcopy__(self, Parameter new, dict memo) cdef __update_state(self) + cdef bint _is_changing(self) cdef _get_value(self) cdef void _set_value(self, value) cdef void set_double_value(self, double value) @@ -40,9 +42,5 @@ cdef class GeometricParameter(Parameter): cdef: bint is_nr - # An array of pointers to changing ABCD matrix - # elements dependent upon this parameter - #cy_expr** dep_abcd_elements - cdef void update_abcd_matrices(self) cdef void set_double_value(self, double value) diff --git a/src/finesse/parameter.pyx b/src/finesse/parameter.pyx index 7669839bfeabeaac4cea57a67c3d1aba7a589957..49557e615a46c041d376047c66de1aa116eb66fa 100644 --- a/src/finesse/parameter.pyx +++ b/src/finesse/parameter.pyx @@ -25,6 +25,7 @@ class ParameterRef(Symbol): def __init__(self, param): self.__param = weakref.ref(param) self.__name = f"{self.owner.name}.{self.parameter.name}" + self.__cyexpr_name = self.__name.replace(".", "_").lower().encode("UTF-8") @property def parameter(self): @@ -34,6 +35,28 @@ class ParameterRef(Symbol): def name(self): return self.__name + @property + def cyexpr_name(self): + """Name of the parameter reference in ``cyexpr`` compatibility + format. + + This is equivalent to :attr:`.ParameterRef.name` but with ``"."`` + replaced with ``"_"``, converted to lower case and encoded + in UTF-8 format as a ``bytes`` object. + + The above format makes this compatible with passing to the underlying + math evaluator engine (``tinyexpr``) used via the :mod:`.cyexpr` sub-module. + + .. note:: + + This should, typically, never need to be used outside of internal usage. It + exists primarily to act as the owner for the parameter name strings (avoiding + dangling pointers in the expression code). + + :getter: Returns the ``cyexpr`` compatible name format of the pref name (read-only). + """ + return self.__cyexpr_name + @property def owner(self): return self.parameter.component @@ -75,12 +98,6 @@ class ParameterRef(Symbol): return result def eval(self, keep_changing_symbols=False, subs=None, **kwargs): - #print(self) - kw_name = self.name.replace(".", "_") - kw_value = kwargs.get(kw_name) - if kw_value is not None: - return kw_value - p = self.parameter if keep_changing_symbols and p.is_changing: @@ -97,6 +114,11 @@ cdef class Parameter: from finesse.element import ModelElement self.__value = value + if value is not None: + self.__cvalue = float(value) + else: + self.__cvalue = 0.0 + self.__units = units self.__name = name self.__full_name = f"{component.name}.{name}" @@ -121,6 +143,7 @@ cdef class Parameter: new.state = self.state new.__rebuild_flag = self.__rebuild_flag new.__value = deepcopy(self.__value, memo) + new.__cvalue = self.__cvalue new.__component_type = self.__component_type new.__validator = self.__validator new.__post_setter = self.__post_setter @@ -180,12 +203,15 @@ cdef class Parameter: def ref(self): return ParameterRef(self) - @property - def is_changing(self): + cdef bint _is_changing(self): if self.state == ParameterState.Symbolic: return self.value.is_changing - else: - return self.__is_tunable + + return self.__is_tunable + + @property + def is_changing(self): + return self._is_changing() @property def is_tunable(self): @@ -300,6 +326,9 @@ cdef class Parameter: else: self.__value = self.__validator(self.owner, value) + if isinstance(value, numbers.Number): + self.__cvalue = value + self.__update_state() if self.__post_setter is not None: @@ -314,11 +343,16 @@ cdef class Parameter: UNSAFE! Only use when you know the parameter should be changing. """ + # TODO (sjr) Not sure setting of self.__value here is reqd + # anymore as just setting __cvalue should suffice + # here now (as this method should only be called + # during simulations anyway) if self.__validator is None: self.__value = value else: self.__value = self.__validator(self.owner, value) + self.__cvalue = value if self.__post_setter is not None: self.__post_setter(self.owner)() @@ -531,8 +565,10 @@ cdef class GeometricParameter(Parameter): Parameter.set_double_value(self, value) # If the parameter is changing that means we're in a simulation so # don't call update_abcd_matrices as ABCD matrix elements are updated - # in an optimised way with cy_expr pointers during a simulation - if not self.is_changing: + # in an optimised way with cy_expr pointers during a simulation. Note + # we also use _is_changing (C call) for speed here instead of is_changing + # property. + if not self._is_changing(): self.update_abcd_matrices() # NOTE (sjr) Seems that Cython doesn't currently support overriding @@ -546,7 +582,7 @@ cdef class GeometricParameter(Parameter): def value(self, value): Parameter._set_value(self, value) # NOTE (sjr) See above comment in set_double_value for logic behind condition here - if not self.is_changing: + if not self._is_changing(): self.update_abcd_matrices() diff --git a/src/finesse/simulations/base.pxd b/src/finesse/simulations/base.pxd index c183477501e10d57effcb17541c9850c9ce0d057..695384ba5ad2a7f95002ec3155949d046ec6288d 100644 --- a/src/finesse/simulations/base.pxd +++ b/src/finesse/simulations/base.pxd @@ -182,8 +182,6 @@ cdef class BaseSimulation: cpdef _initialise_frequencies(self) - cdef bint _initialise_knm_workspaces(self, ConnectorWorkspace ws) - cdef void _determine_changing_beam_params(self) # TODO (sjr) Would prefer to have the below methods in finesse.tracing.ctracer # (see commented out blocks there) but circular import nonsense diff --git a/src/finesse/simulations/base.pyx b/src/finesse/simulations/base.pyx index fdb72b7ef4c344e0f10a1177e556a7500e167312..05c5751f6660777488baf3ef4eeada3cbb5fee15 100644 --- a/src/finesse/simulations/base.pyx +++ b/src/finesse/simulations/base.pyx @@ -18,47 +18,9 @@ from finesse.element cimport ElementWorkspace from finesse.frequency cimport Frequency from finesse.components.workspace cimport ConnectorCallbacks from finesse.components.node import NodeType + +from finesse.components.modal.workspace cimport KnmConnectorWorkspace from finesse.components.modal.cavity cimport CavityWorkspace -from finesse.components.modal.mirror cimport ( - MirrorWorkspace, - initialise_knm_workspaces as init_mirror_knm_ws, - free_knm_workspaces as free_mirror_knm_ws, - flag_changing_knm_workspaces as flag_changing_mirror_knm, - update_changing_knm_workspaces as update_mirror_knm_ws, - compute_scattering_matrices as compute_mirror_knm -) -from finesse.components.modal.beamsplitter cimport ( - BeamsplitterWorkspace, - initialise_knm_workspaces as init_bs_knm_ws, - free_knm_workspaces as free_bs_knm_ws, - flag_changing_knm_workspaces as flag_changing_bs_knm, - update_changing_knm_workspaces as update_bs_knm_ws, - compute_scattering_matrices as compute_bs_knm -) -from finesse.components.modal.lens cimport ( - LensWorkspace, - initialise_knm_workspaces as init_lens_knm_ws, - free_knm_workspaces as free_lens_knm_ws, - flag_changing_knm_workspaces as flag_changing_lens_knm, - update_changing_knm_workspaces as update_lens_knm_ws, - compute_scattering_matrices as compute_lens_knm -) -from finesse.components.modal.modulator cimport ( - ModulatorWorkspace, - initialise_knm_workspaces as init_mod_knm_ws, - free_knm_workspaces as free_mod_knm_ws, - flag_changing_knm_workspaces as flag_changing_mod_knm, - update_changing_knm_workspaces as update_mod_knm_ws, - compute_scattering_matrices as compute_mod_knm -) -from finesse.components.modal.isolator cimport ( - IsolatorWorkspace, - initialise_knm_workspaces as init_isol_knm_ws, - free_knm_workspaces as free_isol_knm_ws, - flag_changing_knm_workspaces as flag_changing_isol_knm, - update_changing_knm_workspaces as update_isol_knm_ws, - compute_scattering_matrices as compute_isol_knm -) from finesse.components.modal.space cimport ( SpaceWorkspace, set_gouy_phase @@ -520,35 +482,10 @@ cdef class BaseSimulation: ws.update() cdef void compute_knm_matrices(self): - cdef: - ConnectorWorkspace ws - MirrorWorkspace mws - BeamsplitterWorkspace bws - LensWorkspace lws - ModulatorWorkspace mdws - IsolatorWorkspace isws - + cdef KnmConnectorWorkspace ws for ws in self.to_scatter_matrix_compute: - if isinstance(ws, MirrorWorkspace): - mws = ws - update_mirror_knm_ws(mws) - compute_mirror_knm(mws) - elif isinstance(ws, BeamsplitterWorkspace): - bws = ws - update_bs_knm_ws(bws) - compute_bs_knm(bws) - elif isinstance(ws, LensWorkspace): - lws = ws - update_lens_knm_ws(lws) - compute_lens_knm(lws) - elif isinstance(ws, ModulatorWorkspace): - mdws = ws - update_mod_knm_ws(mdws) - compute_mod_knm(mdws) - elif isinstance(ws, IsolatorWorkspace): - isws = ws - update_isol_knm_ws(isws) - compute_isol_knm(isws) + ws.update_changing_knm_workspaces() + ws.compute_scattering_matrices() cdef int set_gouy_phases(self) except -1: cdef ConnectorWorkspace ws @@ -696,11 +633,14 @@ cdef class BaseSimulation: # the next few calls as trace data is required for these cdef TraceForest model_trace_forest = self.model._trace_forest if self.is_modal and not self.is_audio: + # Make sure the model trace forest gets re-planted + # when building a new simulation + self.model._rebuild_trace_forest = True self.model.beam_trace(**self.model.beam_trace_args) self.retrace = self.model.retrace if self.retrace and self.is_modal and not self.is_audio: - # construct the forest of changing trace trees - need to do this + # Construct the forest of changing trace trees - need to do this # here as the _get_workspace method of Connectors requires the # simulation trace_forest field to be initialised when making # the refill flags @@ -832,28 +772,12 @@ cdef class BaseSimulation: set(self.electrical_frequencies.values()) ) - - cdef bint _initialise_knm_workspaces(self, ConnectorWorkspace ws): - if isinstance(ws, MirrorWorkspace): - init_mirror_knm_ws(ws) - elif isinstance(ws, BeamsplitterWorkspace): - init_bs_knm_ws(ws) - elif isinstance(ws, LensWorkspace): - init_lens_knm_ws(ws) - elif isinstance(ws, ModulatorWorkspace): - init_mod_knm_ws(ws) - elif isinstance(ws, IsolatorWorkspace): - init_isol_knm_ws(ws) - else: - return False - - return True - cdef void _determine_changing_beam_params(self): cdef: Py_ssize_t i Py_ssize_t num_nodes = len(self.nodes) ConnectorWorkspace ws + KnmConnectorWorkspace kws # re-set all beam parameter changing flags to false initially for i in range(num_nodes): @@ -868,17 +792,8 @@ cdef class BaseSimulation: # now tell each knm workspace whether it is changing or not # so that only changing scattering matrices get recomputed # from here on - for ws in self.to_scatter_matrix_compute: - if isinstance(ws, MirrorWorkspace): - flag_changing_mirror_knm(ws) - elif isinstance(ws, BeamsplitterWorkspace): - flag_changing_bs_knm(ws) - elif isinstance(ws, LensWorkspace): - flag_changing_lens_knm(ws) - elif isinstance(ws, ModulatorWorkspace): - flag_changing_mod_knm(ws) - elif isinstance(ws, IsolatorWorkspace): - flag_changing_isol_knm(ws) + for kws in self.to_scatter_matrix_compute: + kws.flag_changing_knm_workspaces() cdef void _setup_trace_forest(self): cdef: @@ -889,13 +804,6 @@ cdef class BaseSimulation: tree = self.trace_forest.forest[tree_idx] self._setup_single_trace_tree(tree) - # Also need to set the parent node_id as this will be used - # in the beam tracing propagation routine - if tree.parent is not None: - self.trace_forest.forest[tree_idx] = tree.parent - tree.parent.node_id = self.node_id(tree.parent.node) - tree.parent.opp_node_id = self.node_id(tree.parent.node.opposite) - if self.trace_forest.N_trees: LOGGER.info("Determined changing trace trees:%s", self.trace_forest.draw()) @@ -1116,6 +1024,7 @@ cdef class BaseSimulation: cpdef _initialise_elements(self): from finesse.components import Connector, Cavity, Variable + # Store all the edge owners we'll need to loop over and call fill on self._edge_owners = [] # use set for unique edges @@ -1155,8 +1064,10 @@ cdef class BaseSimulation: # if the connector scatters modes then initialise the # knm workspaces here and store the connector workspace # in to_scatter_matrix_compute for future use - if self._initialise_knm_workspaces(ws): + if isinstance(ws, KnmConnectorWorkspace): + ( ws).initialise_knm_workspaces() self.to_scatter_matrix_compute.append(ws) + self.workspaces.append(ws) # store all workspaces here if isinstance(ws, ConnectorWorkspace): @@ -1184,14 +1095,15 @@ cdef class BaseSimulation: # in ElementWorkspace.chprm_expr which is used for fast evaluating # of the changing symbolic expressions for ws in self.workspaces: - ws.compile_cy_exprs(self) + ws.compile_cy_exprs() # Also compile the changing ABCD matrix elements, these are # stored in the relevant cy_expr** field of the associated # workspace -> note that the cy_expr* element is NULL for non # changing elements - ws.compile_abcd_cy_exprs() + if self.is_modal: + ws.compile_abcd_cy_exprs() for vws in self.variable_workspaces: - ws.compile_cy_exprs(self) + ws.compile_cy_exprs() # refil gets hammered hard so lets stick it in a type array rather than a list self.num_to_refill = len(self.to_refill) @@ -1303,16 +1215,8 @@ cdef class BaseSimulation: cpdef _clear_workspace(self): if self.to_scatter_matrix_compute is not None: for ws in self.to_scatter_matrix_compute: - if isinstance(ws, MirrorWorkspace): - free_mirror_knm_ws(ws) - elif isinstance(ws, BeamsplitterWorkspace): - free_bs_knm_ws(ws) - elif isinstance(ws, LensWorkspace): - free_lens_knm_ws(ws) - elif isinstance(ws, ModulatorWorkspace): - free_mod_knm_ws(ws) - elif isinstance(ws, IsolatorWorkspace): - free_isol_knm_ws(ws) + (ws).free_knm_workspaces() + # empty out any workspaces left self.to_scatter_matrix_compute.clear() if self.trace_forest: diff --git a/src/finesse/solutions/beamtrace.py b/src/finesse/solutions/beamtrace.py index 6e88ca93d5e5f23549655840c82ba994a686ef9e..e6fe4fda4c55ac4e8fa589b74ca9f4f581913a35 100644 --- a/src/finesse/solutions/beamtrace.py +++ b/src/finesse/solutions/beamtrace.py @@ -1527,6 +1527,17 @@ class AstigmaticPropagationSolution(BaseSolution): at a given location of the path. See :meth:`.PropagationSolution.q` for parameters description. + + Returns + ------- + O : float or :class:`.Operation` + If `at` is an :class:`.OpticalNode` instance, or a :class:`.Connector` + object *and* `which` was specified, then the plane overlap corresponding to + this node is returned. + + Os : dict + Otherwise, a dictionary of the input and output node mappings to their + respective plane overlaps is returned. """ qxs = self.qx(at, which) qys = self.qy(at, which) diff --git a/src/finesse/symbols.pyx b/src/finesse/symbols.pyx index 66e9e1f51428b50419a0f0b9af2eeccf7f3d22ec..71c2ba187a90d82f8f6a7c503e3e0d0aaadd4d7a 100644 --- a/src/finesse/symbols.pyx +++ b/src/finesse/symbols.pyx @@ -7,6 +7,7 @@ import warnings import numpy as np from finesse import constants +from finesse.utilities import is_iterable import operator import logging @@ -119,7 +120,7 @@ def evaluate(x): A single value for the evaluated expression if `x` is not array-like, otherwise an array of the evaluated expressions. """ - if hasattr(x, "__getitem__"): + if is_iterable(x): y = np.array(x, dtype=np.complex128) if not np.any(y.imag): # purely real symbols in array with warnings.catch_warnings(): @@ -134,76 +135,6 @@ def evaluate(x): return x.eval() -cdef construct_grids(dict kwargs): - """Constructs N-dimensional meshgrids from the arguments - provided. - - If a combination of parameters in the keyword-arguments are - all of type :class:`numpy.ndarray`, then a meshgrid is formed - from these arguments. The dictionary returned from this function - contains the parameter names with updated values (i.e. grids). - - Note that for each combination, if any of the parameters in the - combo are not :class:`numpy.ndarray` instances then the entry in - the returned dictionary for these parameters is unchanged from - the value passed in. - - Parameters - ---------- - kwargs : keyword-arguments - Parameters with their value(s). - - Returns - ------- - formed : dict - A new dictionary with updated values from `kwargs`, or just - `kwargs` itself if the number of array-like values in `kwargs` - is less than two. - - changed : bool - Flag indicating whether any meshgrids were constructed. - """ - - cdef: - bint changed = False # were any meshgrids constructed? - bint all_arrays - int n_arrays = 0 - list combos - dict values = {} - - for v in kwargs.values(): - if isinstance(v, np.ndarray): - n_arrays += 1 - - if n_arrays < 2: # cannot form any meshgrids with fewer than 2 arrays - return kwargs, changed - - combos = list(combinations(kwargs.items(), n_arrays)) - - for kv_pairs in combos: - ks = [kvp[0] for kvp in kv_pairs] # parameter names - vs = [kvp[1] for kvp in kv_pairs] # parameter values - - all_arrays = True - for v in vs: - if not isinstance(v, np.ndarray): - all_arrays = False - break - - if all_arrays: - vs = np.meshgrid(*vs) - for k, v in zip(ks, vs): - values[k] = v - - changed = True - else: - for k, v in zip(ks, vs): - if k not in values: - values[k] = v - - return values, changed - - class Symbol: __add__ = OPERATORS["__add__"] __sub__ = OPERATORS["__sub__"] @@ -374,24 +305,21 @@ class Operation(Symbol): self.args = tuple(as_symbol(_) for _ in args) self.op = operation - def eval(self, bint grids=False, **kwargs): - """Evaluates the operation, substituting parameter values for those - in `kwargs` if present. - - Parameters - ---------- - grids : bool, optional; default: False - Whether to construct meshgrids of multiple array-like parameter - substitutions present in `kwargs`. + def eval(self, **kwargs): + """Evaluates the operation. + + Parameter substitutions can be given via an optional + ``subs`` dict (mapping parameters to substituted values). + + Returns + ------- + result : number or array-like + The single-valued result of evaluation of the operation (if no + substitutions given, or all substitutions are scalar-valued). Otherwise, + if any parameter substitution was a :class:`numpy.ndarray`, then a corresponding array + of results. """ - if grids: - if "__grids_constructed" not in kwargs: - formed, changed = construct_grids(kwargs) - if changed: - kwargs.update(formed) - kwargs["__grids_constructed"] = True - - return self.op(*(_.eval(grids=grids, **kwargs) for _ in self.args)) + return self.op(*(_.eval(**kwargs) for _ in self.args)) def __eq__(self, obj): if isinstance(obj, Operation): diff --git a/src/finesse/tracing/ctracer.pyx b/src/finesse/tracing/ctracer.pyx index 33e91fcdf20452cbd5bf61b7a78c122bd120917f..ffc71b4683f7419371d5bc36e2c320aa44e835e0 100644 --- a/src/finesse/tracing/ctracer.pyx +++ b/src/finesse/tracing/ctracer.pyx @@ -407,6 +407,12 @@ cdef class TraceForest: return self.model.network.nodes[name]["weakref"]() cdef TraceForest make_changing_forest(self): + """Constructs a new TraceForest from this forest, consisting + of only the trees which will have changing beam parameters. + + This method is called in BaseSimulation._initialise for setting + up the simulation trace forest used for efficient beam tracing. + """ cdef: Py_ssize_t tree_idx TraceTree tree @@ -416,8 +422,31 @@ cdef class TraceForest: for tree_idx in range(self.N_trees): tree = self.forest[tree_idx] + # From this tree, obtain the broadest sub-trees which + # will have changing beam parameters changing_trees.extend(tree.get_broadest_changing_subtrees()) + cdef Py_ssize_t Nchanging_nominal = len(changing_trees) + cdef set roots = set() # Parents of changing trees + cdef set trees_to_remove = set() + for tree_idx in range(Nchanging_nominal): + tree = changing_trees[tree_idx] + # Now set each changing tree to the parent tree so that + # the root is used in beam tracing... + if tree.parent is not None: + # ... but only do this for parents not yet added + # otherwise could get duplicate trees + if tree.parent not in roots: + changing_trees[tree_idx] = tree.parent + roots.add(tree.parent) + else: + # This tree will already be encapsulated by the previous + # parent so just mark it to be removed + trees_to_remove.add(tree) + + for tree in trees_to_remove: + changing_trees.remove(tree) + changing_forest.forest = changing_trees changing_forest.N_trees = len(changing_trees) changing_forest.symmetric = self.symmetric diff --git a/src/finesse/utilities/homs.py b/src/finesse/utilities/homs.py index 9a7fc75fd044aed99f9d0a5afaedf8361ff7f009..be66e4848212386a31bef7cd277bef517ece7459 100644 --- a/src/finesse/utilities/homs.py +++ b/src/finesse/utilities/homs.py @@ -71,7 +71,10 @@ def make_modes(select=None, maxtem=None, unique=True): """ if select is None and maxtem is None: - raise ValueError("""Arguments 'select' and 'maxtem' cannot both be None.""") + raise ValueError( + f"Error in make_modes:\n" + f" arguments select and maxtem cannot both be None" + ) if select is None: _check_maxtem(maxtem) diff --git a/src/finesse/utilities/misc.py b/src/finesse/utilities/misc.py index 602c23198623231fb21f3ac44b147e511646cd19..824277b2c642893d6d13ff7e425d024bcec9e98a 100644 --- a/src/finesse/utilities/misc.py +++ b/src/finesse/utilities/misc.py @@ -105,6 +105,10 @@ def find_nearest(x, value, index=False): def is_iterable(obj): """Reliable check for whether an object is iterable. + Note that strings are treated as non-iterable objects + when performing this check. This will only return true + for iterable non-str objects. + Returns ------- flag : bool @@ -115,7 +119,7 @@ def is_iterable(obj): except Exception: return False else: - return True + return not isinstance(obj, str) @contextmanager diff --git a/tests/validation/plane_wave/test_beamsplitter_power_conservation.py b/tests/validation/plane_wave/test_beamsplitter_power_conservation.py new file mode 100644 index 0000000000000000000000000000000000000000..390ce67273606f5ba4cadfb0f61f94b4d6f8ac3e --- /dev/null +++ b/tests/validation/plane_wave/test_beamsplitter_power_conservation.py @@ -0,0 +1,45 @@ +import pytest +import finesse +import numpy as np +from numpy.testing import assert_allclose + +@pytest.fixture +def model(): + kat = finesse.Model() + kat.parse_legacy(""" + l i1 1 0 0 nlaser + s sin 0 nlaser nBS1 + + bs BS 0.5 0.5 0 0 nBS1 nBS2 nBS3 nAS + s sX 0 1 nBS3 nX + s sY 0 1 nBS2 nY + + m mX 1 0 0 nX dump + m mY 1 0 0 nY dump + + pd PIN nlaser* + pd PREFL nlaser + pd PAS nAS + + xaxis BS phi lin 0 1000 10 + yaxis lin abs:deg + + maxtem off + """) + return kat + +def test_vacuum_normal(model): + out = model.run() + assert_allclose(out['PREFL'] + out['PAS'], out['PIN'], rtol=1e-14, atol=1e-14) + +def test_vacuum_AoI(model): + model.BS.alpha = 45 + out = model.run() + assert_allclose(out['PREFL'] + out['PAS'], out['PIN'], rtol=1e-14, atol=1e-14) + +def test_glass_AoI(model): + model.BS.alpha = 45 + model.spaces.sX.nr = 1.45 + model.spaces.sY.nr = 1.45 + out = model.run() + assert_allclose(out['PREFL'] + out['PAS'], out['PIN'], rtol=1e-14, atol=1e-14) \ No newline at end of file