From ce3e47daa985d4c0549a5d47237dc63268eb9715 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Mon, 18 Mar 2024 10:54:51 -0700 Subject: [PATCH] comment out some eme notes --- nom-eme.py | 332 +++++++++++++++++++++++++++-------------------------- 1 file changed, 168 insertions(+), 164 deletions(-) diff --git a/nom-eme.py b/nom-eme.py index 2b89f8a..d9c7c8b 100644 --- a/nom-eme.py +++ b/nom-eme.py @@ -1,169 +1,173 @@ - -from simphony.elements import Model -from simphony.netlist import Subcircuit -from simphony.simulation import SweepSimulation - -from matplotlib import pyplot as plt +import scipy +import numpy +from numpy.typing import ArrayLike, NDarray -class PeriodicLayer(Model): - def __init__(self, left_modes, right_modes, s_params): - self.left_modes = left_modes - self.right_modes = right_modes - self.left_ports = len(self.left_modes) - self.right_ports = len(self.right_modes) - self.normalize_fields() - self.s_params = s_params - - def normalize_fields(self): - for mode in range(len(self.left_modes)): - self.left_modes[mode].normalize() - for mode in range(len(self.right_modes)): - self.right_modes[mode].normalize() - - -class PeriodicEME: - def __init__(self, layers=[], num_periods=1): - self.layers = layers - self.num_periods = num_periods - self.wavelength = wavelength - - def propagate(self): - wl = self.wavelength - if not len(self.layers): - raise Exception("Must place layers before propagating") - - num_modes = max([l.num_modes for l in self.layers]) - iface = InterfaceSingleMode if num_modes == 1 else InterfaceMultiMode - - eme = EME(layers=self.layers) - left, right = eme.propagate() - self.single_period = eme.s_matrix - - period_layer = PeriodicLayer(left.modes, right.modes, self.single_period) - current_layer = PeriodicLayer(left.modes, right.modes, self.single_period) - interface = iface(right, left) - - for _ in range(self.num_periods - 1): - current_layer.s_params = cascade(current_layer, interface, wl) - current_layer.s_params = cascade(current_layer, period_layer, wl) - - self.s_params = current_layer.s_params - - -class EME: - def __init__(self, layers=[]): - self.layers = layers - self.wavelength = None - - def propagate(self): - layers = self.layers - wl = layers[0].wavelength if self.wavelength is None else self.wavelength - if not len(layers): - raise Exception("Must place layers before propagating") - - num_modes = max([l.num_modes for l in layers]) - iface = InterfaceSingleMode if num_modes == 1 else InterfaceMultiMode - - first_layer = layers[0] - current = Current(wl, first_layer) - interface = iface(first_layer, layers[1]) - - current.s = cascade(current, interface, wl) - current.right_pins = interface.right_pins - - for index in range(1, len(layers) - 1): - layer1 = layers[index] - layer2 = layers[index + 1] - interface = iface(layer1, layer2) - - current.s = cascade(current, layer1, wl) - current.right_pins = layer1.right_pins - - current.s = cascade(current, interface, wl) - current.right_pins = interface.right_pins - - last_layer = layers[-1] - current.s = cascade(current, last_layer, wl) - current.right_pins = last_layer.right_pins - - self.s_matrix = current.s - return first_layer, last_layer - - -def stack(sa, sb): - qab = numpy.eye() - sa.r11 @ sb.r11 - qba = numpy.eye() - sa.r11 @ sb.r11 - #s.t12 = sa.t12 @ numpy.pinv(qab) @ sb.t12 - #s.r21 = sa.t12 @ numpy.pinv(qab) @ sb.r22 @ sa.t21 + sa.r22 - #s.r12 = sb.t21 @ numpy.pinv(qba) @ sa.r11 @ sb.t12 + sb.r11 - #s.t21 = sb.t21 @ numpy.pinv(qba) @ sa.t21 - s.t12 = sa.t12 @ numpy.linalg.solve(qab, sb.t12) - s.r21 = sa.t12 @ numpy.linalg.solve(qab, sb.r22 @ sa.t21) + sa.r22 - s.r12 = sb.t21 @ numpy.linalg.solve(qba, sa.r11 @ sb.t12) + sb.r11 - s.t21 = sb.t21 @ numpy.linalg.solve(qba, sa.t21) - return s - - -def cascade(first, second, wavelength): - circuit = Subcircuit("Device") - - circuit.add([(first, "first"), (second, "second")]) - for port in range(first.right_ports): - circuit.connect("first", "right" + str(port), "second", "left" + str(port)) - - simulation = SweepSimulation(circuit, wavelength, wavelength, num=1) - result = simulation.simulate() - return result.s - - -class InterfaceSingleMode(Model): - def __init__(self, layer1, layer2, num_modes=1): - self.num_modes = num_modes - self.num_ports = 2 * num_modes - self.s = self.solve(layer1, layer2, num_modes) - - def solve(self, layer1, layer2, num_modes): - nm = num_modes - s = numpy.zeros((2 * nm, 2 * nm), dtype=complex) - - for ii, left_mode in enumerate(layer1.modes): - for oo, right_mode in enumerate(layer2.modes): - r, t = get_rt(left_mode, right_mode) - s[ oo, ii] = r - s[nm + oo, ii] = t - - for ii, right_mode in enumerate(layer2.modes): - for oo, left_mode in enumerate(layer1.modes): - r, t = get_rt(right_mode, left_mode) - s[ oo, nm + ii] = t - s[nm + oo, nm + ii] = r - return s - - -class InterfaceMultiMode(Model): - def __init__(self, layer1, layer2): - self.s = self.solve(layer1, layer2) - - def solve(self, layer1, layer2): - n1p = layer1.num_modes - n2p = layer2.num_modes - num_ports = n1p + n2p - s = numpy.zeros((num_ports, num_ports), dtype=complex) - - for l1p in range(n1p): - ts = get_t(l1p, layer1, layer2) - rs = get_r(l1p, layer1, layer2, ts) - s[n1p:, l1p] = ts - s[:n1p, l1p] = rs - - for l2p in range(n2p): - ts = get_t(l2p, layer2, layer1) - rs = get_r(l2p, layer2, layer1, ts) - s[:n1p, n1p + l2p] = ts - s[n1p:, n1p + l2p] = rs - - return s +#from simphony.elements import Model +#from simphony.netlist import Subcircuit +#from simphony.simulation import SweepSimulation +# +#from matplotlib import pyplot as plt +# +# +#class PeriodicLayer(Model): +# def __init__(self, left_modes, right_modes, s_params): +# self.left_modes = left_modes +# self.right_modes = right_modes +# self.left_ports = len(self.left_modes) +# self.right_ports = len(self.right_modes) +# self.normalize_fields() +# self.s_params = s_params +# +# def normalize_fields(self): +# for mode in range(len(self.left_modes)): +# self.left_modes[mode].normalize() +# for mode in range(len(self.right_modes)): +# self.right_modes[mode].normalize() +# +# +#class PeriodicEME: +# def __init__(self, layers=[], num_periods=1): +# self.layers = layers +# self.num_periods = num_periods +# self.wavelength = wavelength +# +# def propagate(self): +# wl = self.wavelength +# if not len(self.layers): +# raise Exception("Must place layers before propagating") +# +# num_modes = max([l.num_modes for l in self.layers]) +# iface = InterfaceSingleMode if num_modes == 1 else InterfaceMultiMode +# +# eme = EME(layers=self.layers) +# left, right = eme.propagate() +# self.single_period = eme.s_matrix +# +# period_layer = PeriodicLayer(left.modes, right.modes, self.single_period) +# current_layer = PeriodicLayer(left.modes, right.modes, self.single_period) +# interface = iface(right, left) +# +# for _ in range(self.num_periods - 1): +# current_layer.s_params = cascade(current_layer, interface, wl) +# current_layer.s_params = cascade(current_layer, period_layer, wl) +# +# self.s_params = current_layer.s_params +# +# +#class EME: +# def __init__(self, layers=[]): +# self.layers = layers +# self.wavelength = None +# +# def propagate(self): +# layers = self.layers +# wl = layers[0].wavelength if self.wavelength is None else self.wavelength +# if not len(layers): +# raise Exception("Must place layers before propagating") +# +# num_modes = max([l.num_modes for l in layers]) +# iface = InterfaceSingleMode if num_modes == 1 else InterfaceMultiMode +# +# first_layer = layers[0] +# current = Current(wl, first_layer) +# interface = iface(first_layer, layers[1]) +# +# current.s = cascade(current, interface, wl) +# current.right_pins = interface.right_pins +# +# for index in range(1, len(layers) - 1): +# layer1 = layers[index] +# layer2 = layers[index + 1] +# interface = iface(layer1, layer2) +# +# current.s = cascade(current, layer1, wl) +# current.right_pins = layer1.right_pins +# +# current.s = cascade(current, interface, wl) +# current.right_pins = interface.right_pins +# +# last_layer = layers[-1] +# current.s = cascade(current, last_layer, wl) +# current.right_pins = last_layer.right_pins +# +# self.s_matrix = current.s +# return first_layer, last_layer +# +# +#def stack(sa, sb): +# qab = numpy.eye() - sa.r11 @ sb.r11 +# qba = numpy.eye() - sa.r11 @ sb.r11 +# #s.t12 = sa.t12 @ numpy.pinv(qab) @ sb.t12 +# #s.r21 = sa.t12 @ numpy.pinv(qab) @ sb.r22 @ sa.t21 + sa.r22 +# #s.r12 = sb.t21 @ numpy.pinv(qba) @ sa.r11 @ sb.t12 + sb.r11 +# #s.t21 = sb.t21 @ numpy.pinv(qba) @ sa.t21 +# s.t12 = sa.t12 @ numpy.linalg.solve(qab, sb.t12) +# s.r21 = sa.t12 @ numpy.linalg.solve(qab, sb.r22 @ sa.t21) + sa.r22 +# s.r12 = sb.t21 @ numpy.linalg.solve(qba, sa.r11 @ sb.t12) + sb.r11 +# s.t21 = sb.t21 @ numpy.linalg.solve(qba, sa.t21) +# return s +# +# +#def cascade(first, second, wavelength): +# circuit = Subcircuit("Device") +# +# circuit.add([(first, "first"), (second, "second")]) +# for port in range(first.right_ports): +# circuit.connect("first", "right" + str(port), "second", "left" + str(port)) +# +# simulation = SweepSimulation(circuit, wavelength, wavelength, num=1) +# result = simulation.simulate() +# return result.s +# +# +#class InterfaceSingleMode(Model): +# def __init__(self, layer1, layer2, num_modes=1): +# self.num_modes = num_modes +# self.num_ports = 2 * num_modes +# self.s = self.solve(layer1, layer2, num_modes) +# +# def solve(self, layer1, layer2, num_modes): +# nm = num_modes +# s = numpy.zeros((2 * nm, 2 * nm), dtype=complex) +# +# for ii, left_mode in enumerate(layer1.modes): +# for oo, right_mode in enumerate(layer2.modes): +# r, t = get_rt(left_mode, right_mode) +# s[ oo, ii] = r +# s[nm + oo, ii] = t +# +# for ii, right_mode in enumerate(layer2.modes): +# for oo, left_mode in enumerate(layer1.modes): +# r, t = get_rt(right_mode, left_mode) +# s[ oo, nm + ii] = t +# s[nm + oo, nm + ii] = r +# return s +# +# +#class InterfaceMultiMode(Model): +# def __init__(self, layer1, layer2): +# self.s = self.solve(layer1, layer2) +# +# def solve(self, layer1, layer2): +# n1p = layer1.num_modes +# n2p = layer2.num_modes +# num_ports = n1p + n2p +# s = numpy.zeros((num_ports, num_ports), dtype=complex) +# +# for l1p in range(n1p): +# ts = get_t(l1p, layer1, layer2) +# rs = get_r(l1p, layer1, layer2, ts) +# s[n1p:, l1p] = ts +# s[:n1p, l1p] = rs +# +# for l2p in range(n2p): +# ts = get_t(l2p, layer2, layer1) +# rs = get_r(l2p, layer2, layer1, ts) +# s[:n1p, n1p + l2p] = ts +# s[n1p:, n1p + l2p] = rs +# +# return s def get_t(p, left, right):