From a1db2d9bc41ea7232c77d610529002ebaf4f3dae Mon Sep 17 00:00:00 2001 From: = <=> Date: Mon, 8 Jun 2026 15:46:41 +0800 Subject: [PATCH] repo build --- simulation/DualPortElements.py | 734 ++++++++++++++++++ simulation/Lumerical/GDS_SIMU_DEVICE_2X2.lsf | 295 +++++++ simulation/Lumerical/mx_analysis_lib.lsf | 212 +++++ simulation/Lumerical/mx_devices_lib.lsf | 274 +++++++ simulation/Lumerical/mx_frames_lib.lsf | 56 ++ simulation/Lumerical/mx_function_lib.lsf | 105 +++ simulation/Lumerical/mx_lib_install.lsf | 27 + simulation/Lumerical/mx_poly_spiral_lib.lsf | 265 +++++++ simulation/Lumerical/mx_simulation_lib.lsf | 207 +++++ simulation/Lumerical/mx_structures_lib.lsf | 205 +++++ simulation/LumericalPATH.json | 13 + simulation/__init__.py | 6 + .../DualPortElements.cpython-310.pyc | Bin 0 -> 16024 bytes .../DualPortElements.cpython-312.pyc | Bin 0 -> 30708 bytes .../DualPortElements.cpython-38.pyc | Bin 0 -> 15929 bytes .../DualPortElements.cpython-39.pyc | Bin 0 -> 15803 bytes .../__pycache__/__init__.cpython-310.pyc | Bin 0 -> 515 bytes .../__pycache__/__init__.cpython-312.pyc | Bin 0 -> 548 bytes .../__pycache__/__init__.cpython-38.pyc | Bin 0 -> 509 bytes .../__pycache__/__init__.cpython-39.pyc | Bin 0 -> 475 bytes 20 files changed, 2399 insertions(+) create mode 100644 simulation/DualPortElements.py create mode 100644 simulation/Lumerical/GDS_SIMU_DEVICE_2X2.lsf create mode 100644 simulation/Lumerical/mx_analysis_lib.lsf create mode 100644 simulation/Lumerical/mx_devices_lib.lsf create mode 100644 simulation/Lumerical/mx_frames_lib.lsf create mode 100644 simulation/Lumerical/mx_function_lib.lsf create mode 100644 simulation/Lumerical/mx_lib_install.lsf create mode 100644 simulation/Lumerical/mx_poly_spiral_lib.lsf create mode 100644 simulation/Lumerical/mx_simulation_lib.lsf create mode 100644 simulation/Lumerical/mx_structures_lib.lsf create mode 100644 simulation/LumericalPATH.json create mode 100644 simulation/__init__.py create mode 100644 simulation/__pycache__/DualPortElements.cpython-310.pyc create mode 100644 simulation/__pycache__/DualPortElements.cpython-312.pyc create mode 100644 simulation/__pycache__/DualPortElements.cpython-38.pyc create mode 100644 simulation/__pycache__/DualPortElements.cpython-39.pyc create mode 100644 simulation/__pycache__/__init__.cpython-310.pyc create mode 100644 simulation/__pycache__/__init__.cpython-312.pyc create mode 100644 simulation/__pycache__/__init__.cpython-38.pyc create mode 100644 simulation/__pycache__/__init__.cpython-39.pyc diff --git a/simulation/DualPortElements.py b/simulation/DualPortElements.py new file mode 100644 index 0000000..4c4574a --- /dev/null +++ b/simulation/DualPortElements.py @@ -0,0 +1,734 @@ +import nazca as nd +import numpy as np +import os +import json +import time +import sys +import h5py +import matplotlib.pyplot as plt +from typing import Literal + +def __getLumericalLibPATH__(): + path = os.path.dirname(os.path.abspath(__file__)) + "\\Lumerical\\" + return path + +def __checkLumericalDIR__(): + path = os.path.dirname(os.path.abspath(__file__)) + with open(path+"\\LumericalPATH.json") as FID: + FILE_CONTEXT = json.load(FID) + + DEFAULT_PATH = FILE_CONTEXT["DEFAULT_PATH"] + for idx in range(len(DEFAULT_PATH)): + if (os.path.exists(DEFAULT_PATH[idx])): + return DEFAULT_PATH[idx] + + raise Exception("NO FDTD link found using default paths") + + +def __PathGenerate__(dev_name, path=None): + """ creating folders """ + if (path is None): + folder = f"{dev_name}_simu\\" + elif (isinstance(path,str)): + if (path.endswith("\\")): + folder = f"{path}{dev_name}_simu\\" + else: + folder = f"{path}\\{dev_name}_simu\\" + + try: + os.makedirs(folder) + except: + pass + + return folder + +def __LoggAllAttrs__(device,folder,dev_name): + attrs = dir(device) + time_curr = time.strftime('%Y%m%d-%H%M%S', time.localtime()) + attrDict = {} + attrDict["time"] = time_curr + """ recording ALL parameters inside the device """ + for attr in attrs: + if not attr.startswith("__"): + attrCur = getattr(device,attr) + if (isinstance(attrCur,float) or isinstance(attrCur,int)): + attrDict[attr] = attrCur + + elif (isinstance(attrCur,list)): + if (isinstance(attrCur[0],float) or isinstance(attrCur[0],int)): + attrDict[attr] = attrCur + + elif (isinstance(attrCur,np.ndarray)): + attrDict[attr] = attrCur.tolist() + + + """ """ + with open(folder+f"{dev_name}_mxpic.json","w") as FID: + json.dump(attrDict,FID) + + device = device.cell + +def PortParas(pin,width,height,radius=0): + port_dict = {} + port_dict["x"] = pin.x + port_dict["y"] = pin.y + port_dict["a"] = pin.a + port_dict["width"] = width + port_dict["height"] = height + port_dict["radius"] = radius + + + return port_dict + +def MonitorParas(x,y,z,dx,dy,dz): + mont_dict = {} + mont_dict["x"] = x + mont_dict["y"] = y + mont_dict["z"] = z + mont_dict["dx"] = dx + mont_dict["dy"] = dy + mont_dict["dz"] = dz + + return mont_dict + +def DEVICE_2X2_FDTD_INIT(fdtd,run=False,instrcutPATH=None,LibPATH=None): + + fdtd.eval("newproject;") + fdtd.eval("switchtolayout;") + fdtd.eval("deleteall;") + fdtd.eval("clc;") + + if (LibPATH is None): + fdtd.eval("ABSOLUTE_LIB_DIR = read(\'lib_path.txt\');") + elif (isinstance(LibPATH,str)): + if (LibPATH.endswith("\\")): + fdtd.eval("ABSOLUTE_LIB_DIR = \'"+LibPATH+"\';") + + else: + fdtd.eval("ABSOLUTE_LIB_DIR = \'"+LibPATH+"\\\\\';") + ##### Install maxwell's library ##### + + fdtd.eval("PATH_LIB = ABSOLUTE_LIB_DIR + \'mx_lib_install.lsf\';") + fdtd.eval("feval(PATH_LIB);") + + # fdtd.eval("feval(\'GDS_SIMU_DEVICE_2X2.lsf\');") + + fdtd.eval("FUNC_DEVICE_2X2(\'"+instrcutPATH+"\');") + + if (run): + + fdtd.eval("run;") + fdtd.eval("DATA_RETRIEVE_DEVICE_2X2(\'"+instrcutPATH+"\');") + +def tuple_to_complex(t): + return complex(t[0], t[1]) + +def SimuDataFigurePlot(simuPath,devName,saveFlag=True,ports=["a1","b1"]): + + if (simuPath.endswith("\\")): + pass + else: + simuPath = simuPath + "\\" + + """ data analysis """ + data = h5py.File(simuPath + devName + "_results.mat","r") + + dataDict = {} + for portName in ports: + """ Getting result of transmission """ + dataDict[portName] = {} + dataDict[portName]["trans"] = np.squeeze(data[portName]["power"]["T"][()]) + dataDict[portName]["wl"] = np.squeeze(data[portName]["power"]["lambda"][()]) + dataDict[portName]["modes"] = np.squeeze(data[portName]["modes"]["T_net"][()]) + dataDict[portName]["E"] = np.squeeze(data[portName]["E"]["E"][()]) + dataDict[portName]["H"] = np.squeeze(data[portName]["H"]["H"][()]) + + """ Calculating the propagation field """ + E_prg = np.squeeze(data["z1"]["E"]["E"][()]) + wl_prg = np.squeeze(data["z1"]["E"]["lambda"][()]) + + x = np.squeeze(data["z1"]["E"]["x"][()]) + y = np.squeeze(data["z1"]["E"]["y"][()]) + z = np.squeeze(data["z1"]["E"]["z"][()]) + + E_prg = np.vectorize(tuple_to_complex)(E_prg) + + # plt.figure(figsize=(12,6)) + fig,ax = plt.subplots(3,3,figsize=(20, 9)) + plt.subplots_adjust(wspace=0.5, hspace=0.3) + + manager = plt.get_current_fig_manager() + manager.window.wm_geometry("+100+100") + + wl_idx_plt = [0,int(np.floor(len(wl_prg)/2+1))-1,len(wl_prg)-1] + for idx in range(len(wl_idx_plt)): + # plt.subplot(3,2,2*idx+1) + Ex = np.abs(E_prg[idx,0,:].reshape(len(y),len(x))) + Ey = np.abs(E_prg[idx,1,:].reshape(len(y),len(x))) + Ez = np.abs(E_prg[idx,2,:].reshape(len(y),len(x))) + + E_mag = np.sqrt(np.square(Ex) + np.square(Ey) + np.square(Ez)) + + ax[0,idx].set_title(f"wavelength = {wl_prg[wl_idx_plt[idx]]*1e+9:.1f} nm") + ax[0,idx].pcolor(np.real(E_mag)) + + """ Plotting the port transmission """ + if ("b1" in ports and "a1" in ports): + dataDict["Ephase_11"] = np.squeeze(data["Ephase_11"][()]) + Trans_11 = dataDict["b1"]["trans"] + Tmodes = dataDict["b1"]["modes"] + + ax1 = ax[1,0] + ax2 = ax[2,0] + + ax1.set_title("a_1 to b_1 trans [Through]") + ax1.plot(dataDict["a1"]["wl"]*1e+6,Trans_11,linewidth=3) + + """ plotting the eigen mode decomposition """ + dataSZ = np.shape(Tmodes) + plt_wl = dataDict["a1"]["wl"]*1e+6 + + dataPlt = [] + + """ Plotting the Mode crosstalk for target modes """ + if (len(dataSZ) > 1): + for pltIdx in range(0,dataSZ[0]): + _data_ = np.abs(Tmodes[pltIdx,:]) + dataPlt.append(_data_) + else: + _data_ = np.abs(Tmodes) + dataPlt.append(_data_) + + for pltIdx in range(0,len(dataPlt)): + _data_ = dataPlt[pltIdx] + ax1.plot(plt_wl,_data_,linewidth=3) + Trans_dB = 10*np.log10(_data_) + ax2.plot(plt_wl,Trans_dB,label=f"mode_{pltIdx}",linewidth=3) + + + if ("b2" in ports and "a1" in ports): + dataDict["Ephase_21"] = np.squeeze(data["Ephase_11"][()]) + Trans_21 = dataDict["b2"]["trans"] + Tmodes = dataDict["b2"]["modes"] + + ax1 = ax[1,1] + ax2 = ax[2,1] + + ax1.set_title("a_1 to b_1 trans [Through]") + ax1.plot(dataDict["a1"]["wl"]*1e+6,Trans_21) + + """ plotting the eigen mode decomposition """ + dataSZ = np.shape(Tmodes) + plt_wl = dataDict["a1"]["wl"]*1e+6 + + dataPlt = [] + + """ Plotting the Mode crosstalk for target modes """ + if (len(dataSZ) > 1): + for pltIdx in range(0,dataSZ[0]): + _data_ = np.abs(Tmodes[pltIdx,:]) + dataPlt.append(_data_) + else: + _data_ = np.abs(Tmodes) + dataPlt.append(_data_) + + for pltIdx in range(0,len(dataPlt)): + _data_ = dataPlt[pltIdx] + ax1.plot(plt_wl,_data_,linewidth=3) + Trans_dB = 10*np.log10(_data_) + ax2.plot(plt_wl,Trans_dB,label=f"mode_{pltIdx}",linewidth=3) + + if ("a2" in ports and "a1" in ports): + Refl_21 = dataDict["a2"]["trans"] + # fig,ax = plt.subplots(3,2,6) + ax[1,2].set_title("a_1 to a_2 trans [Replection]") + ax[1,2].plot(dataDict["a1"]["wl"]*1e+6,Refl_21,linewidth=3) + + ax2 = ax[2,2] + Trans_dB = 10*np.log10(Refl_21) + ax2.plot(dataDict["a1"]["wl"]*1e+6,Trans_dB,linewidth=3) + + if (saveFlag): + """ in CPU mode, there will be no folder """ + try: + os.makedirs(simuPath + devName + "_simu\\") + except: + pass + plt.savefig(simuPath + devName + "_simu\\"+devName+"_results.jpg", format='jpg', dpi=300) + plt.close() + + + data.close() + +class DEVICE_PORTS: + def __init__(self,dev_name,device,simu_xs="strip", + port_width=3,path=None,wl=[1.5,1.6], + mesh_order=5,layer_heights=[0.22],FDTD_height=2, + material="Si (Silicon) - Palik", + CladMaterial = "SiO2 (Glass) - Palik", + modeIdx=[1,2,3,4], + sourceMode = 1, + ports_extend=["a1"], + SimuBox = None, + + port_radius={"a1":0}, + sample_points = 101, + Field_sample = 3, + + FDTDBuild = False, + LumericalPATH = None, + runFDTD = False, + GPUOn = True, + port_names = ["a1","b1","a2","b2"], + ) -> None: + + + self.dev_name = dev_name + + + time_curr = time.strftime('%Y%m%d-%H%M%S', time.localtime()) + """ creating folders """ + folder = __PathGenerate__(path=path,dev_name=dev_name) + + if (isinstance(device,nd.Cell)): + device = device + elif (hasattr(device,"cell")): + __LoggAllAttrs__(device=device,folder=folder,dev_name=dev_name) + device = device.cell + """ """ + else: + raise Exception("Argument type not recongized") + + if (dev_name == device.cell_name): + + self.dev_name = dev_name + "_GDS" + dev_name = self.dev_name + + """ exporting GDS """ + fname = f"{dev_name}.gds" + + with nd.Cell(name=dev_name) as CELL_INSTR: + instr = device.put() + instr.raise_pins() + for name in ports_extend: + if (device.pin[name].xs is None): + _xs_extend_ = simu_xs + else: + _xs_extend_ = device.pin[name].xs + + nd.strt(xs=_xs_extend_,length=port_width*2,width=device.pin[name].width).put(instr.pin[name]) + nd.text(layer=1001,text=time_curr,height=20).put() + + nd.export_gds(filename=folder + fname,topcells=CELL_INSTR,flat=True) + + + jsonFile = {} + + jsonFile["ports"] = {} + jsonFile["mont"] = {} + jsonFile["input"] = {} + jsonFile["layers"] = {} + + jsonFile["ports"]["names"] = port_names + + for name in jsonFile["ports"]["names"]: + jsonFile["ports"][name] = PortParas(pin=CELL_INSTR.pin[name],width=port_width,height=FDTD_height) + + # jsonFile["ports"]["a1"] = PortParas(pin=CELL_INSTR.pin['a1'],width=port_width,height=FDTD_height) + # jsonFile["ports"]["a2"] = PortParas(pin=CELL_INSTR.pin['a2'],width=port_width,height=FDTD_height) + # jsonFile["ports"]["b1"] = PortParas(pin=CELL_INSTR.pin['b1'],width=port_width,height=FDTD_height) + # jsonFile["ports"]["b2"] = PortParas(pin=CELL_INSTR.pin['b2'],width=port_width,height=FDTD_height) + + for key in port_radius: + jsonFile["ports"][key]["radius"] = port_radius[key] + + """ port Z for propagation recording """ + ports = jsonFile["ports"] + dx = abs(ports["a1"]["x"] - ports["b1"]["x"]) + cX = (ports["a1"]["x"] + ports["b1"]["x"])/2 + + if ("b2" in jsonFile["ports"]["names"]): + dy = abs(ports["b1"]["y"] - ports["b2"]["y"]) + cY = (ports["b1"]["y"] + ports["b2"]["y"])/2 + elif (SimuBox is not None): + dy = SimuBox["dy"] + cY = ports["b1"]["y"] + else: + dy = 0 + cY = ports["b1"]["y"] + + + FDTD = {} + FDTD["x"] = cX + FDTD["y"] = cY + + FDTD["dx"] = dx + FDTD["dy"] = dy + FDTD["z"] = 0 + + SimuBoxKeys = ["x","y","dx","dy"] + if (SimuBox is not None): + for key in SimuBoxKeys: + if (key in SimuBox.keys()): + # if ("dx" in SimuBox.keys()): + FDTD[key] = SimuBox[key] + # FDTD["dy"] = SimuBox["dy"] + # if ("x" in SimuBox.keys()): + # FDTD["x"] = SimuBox["x"] + # FDTD["y"] = SimuBox["y"] + + """ expansion """ + FDTD["dx"] = FDTD["dx"] + port_width + FDTD["dy"] = FDTD["dy"] + port_width + + FDTD["dz"] = FDTD_height + FDTD["wl"] = wl + FDTD["mesh_order"] = mesh_order + FDTD["sourceMode"] = sourceMode + if (GPUOn is True): + FDTD["GPUOn"] = 1 + else: + FDTD["GPUOn"] = 0 + + FDTD["Trans_sample_points"] = sample_points + FDTD["Field_sample_points"] = Field_sample + + jsonFile["mont"]["z1"] = MonitorParas(x=FDTD["x"],y=FDTD["y"],z=0,dx=FDTD["dx"],dy=FDTD["dy"],dz=0) + + """ exporting json configure files """ + + jsonFile["FDTD"] = FDTD + + jsonFile["wafer"] = {} + jsonFile["wafer"]["material"] = material + + jsonFile["clad"] = {} + jsonFile["clad"]["material"] = CladMaterial + + + jsonFile["time"] = time_curr + + jsonFile["geometry"] = {} + + jsonFile["modes"] = modeIdx + + layerNumList = [] + layerDataTypeList = [] + growxList = [] + + + for layers,growx,growy,acc in nd.layeriter(xs=simu_xs): + (a1,b1), (a2,b2),c1,c2 = growx + + layerInfo = nd.get_layer_tuple(layers) + + layerNumList.append(layerInfo.layer) + layerDataTypeList.append(layerInfo.datatype) + growxList.append(b1) + + + if (len(layer_heights) != len(layerNumList)): + raise Exception("Input wafer heigh list is not same length as layer numbers") + + jsonFile["layers"]["numbers"] = layerNumList + jsonFile["layers"]["datatype"] = layerDataTypeList + jsonFile["layers"]["growth"] = growxList + jsonFile["layers"]["heights"] = layer_heights + + + + + """ """ + jsonPath = folder+f"{dev_name}.json" + with open(jsonPath,"w") as FID: + json.dump(jsonFile,FID) + + jsonPath = os.path.abspath(jsonPath) + jsonPath = jsonPath.replace('\\',"\\\\") + self.jsonPath = jsonPath + + + + + GDSPath = folder + fname + GDSPath = os.path.abspath(GDSPath) + GDSPath = GDSPath.replace('\\',"\\\\") + self.GDSPath = GDSPath + + FolderPath = folder + FolderPath = os.path.abspath(FolderPath) + FolderPath = FolderPath.replace('\\',"\\\\") + self.FolderPath = FolderPath + + """ writing instruction to Lumerical for operation """ + instPath = f"{self.FolderPath}\\Instruction.txt" + + self.instPath = instPath + + with open(file=instPath,mode="w") as FInst: + FInst.write(f'devName = \"{self.dev_name}\";\n\r') + FInst.write(f'gdspath = \"{self.GDSPath}\";\n\r') + FInst.write(f'folder = \"{self.FolderPath}\";\n\r') + FInst.write(f'jsonPath = \"{self.jsonPath}\";\n\r') + + + + """ Internal Building FDTD """ + if (FDTDBuild): + if (LumericalPATH is None): + LuPATH = __checkLumericalDIR__() + else: + if (not os.path.exists(LumericalPATH)): + raise Exception("No Lumerical installation found in the given paths") + LuPATH = LumericalPATH + + sys.path.append(LuPATH) + # sys.path.append(os.path.dirname(__file__)) #Current directory + + # print(sys.path) + + import lumapi + + fdtd = lumapi.FDTD() + + """ constructing simulation files """ + print("Building FDTD project ... ") + DEVICE_2X2_FDTD_INIT(fdtd=fdtd,run=runFDTD,LibPATH=__getLumericalLibPATH__(),instrcutPATH=self.instPath) + print("... FDTD Project to ",self.FolderPath) + + + fdtd.close() + + if (runFDTD): + SimuDataFigurePlot(simuPath=self.FolderPath,devName=self.dev_name,ports=port_names) + self.resultPath = self.FolderPath + "\\" + self.dev_name + "_results.mat" + +class DEVICE_RING_BUS(DEVICE_PORTS): + def __init__(self, dev_name, device, r_ring,port_distance=6,Aport=None, + simu_xs="strip", port_width=3, + path=None, wl=[1.5, 1.6], mesh_order=5, layer_heights=[0.22], + FDTD_height=2, + material="Si (Silicon) - Palik", + CladMaterial="SiO2 (Glass) - Palik", + modeIdx=[1, 2, 3, 4], + sample_points=101, + FDTDBuild = False, + LumericalPATH = None, + GPUOn = True, + runFDTD = False,) -> None: + + + if (isinstance(device,nd.Cell)): + cell_dev = device + elif (hasattr(device,"cell")): + cell_dev = device.cell + else: + raise Exception("ERROR :: not recongized") + + + dx = abs(cell_dev.pin["a1"].x - cell_dev.pin["b1"].x)+port_width + + if (Aport is None): + + if (cell_dev.pin['b1'].x > r_ring): + cell_dev.pin['b2'] = nd.Pin(name="b2").put(r_ring,0, 90) + cell_dev.pin['a2'] = nd.Pin(name="a2").put(-r_ring,0, 90) + else: + x = cell_dev.pin['b1'].x + y = -np.sqrt(r_ring**2 - x**2) + a = np.arcsin(x/r_ring)/np.pi*180 + dy_ports = abs(y-cell_dev.pin['b1'].y) + + if (dy_ports > port_distance): + + y = cell_dev.pin['b1'].y + port_distance + x = np.sqrt(r_ring**2 - (abs(y))**2) + a = np.arcsin(x/r_ring)/np.pi*180 + cell_dev.pin['b2'] = nd.Pin(name="b2").put( x,y, a) + cell_dev.pin['a2'] = nd.Pin(name="a2").put(-x,y,180-a) + + else : + cell_dev.pin['b2'] = nd.Pin(name="b2").put( x,y, a) + cell_dev.pin['a2'] = nd.Pin(name="a2").put(-x,y,180-a) + + else : + cell_dev.pin['b2'] = nd.Pin(name="b2").put( r_ring*np.sin(Aport/180*np.pi),-r_ring*np.cos(Aport/180*np.pi), Aport) + cell_dev.pin['a2'] = nd.Pin(name="a2").put(-r_ring*np.sin(Aport/180*np.pi),-r_ring*np.cos(Aport/180*np.pi), Aport) + + yMax = cell_dev.pin['b2'].y + yMin = -r_ring - port_width + + dy = abs(yMax - yMin) + cy = (yMax + yMin)/2 + + if (isinstance(device,nd.Cell)): + device = cell_dev + elif (hasattr(device,"cell")): + device.cell = cell_dev + + + super().__init__(dev_name=dev_name, device=device, simu_xs=simu_xs, port_width=port_width, path=path, wl=wl, + mesh_order=mesh_order, layer_heights=layer_heights, + FDTD_height=FDTD_height, material=material, CladMaterial=CladMaterial, + modeIdx=modeIdx, ports_extend=["a1","b1"], SimuBox = {"dx":dx,"dy":dy,"y":cy}, port_radius={"a2":r_ring,"b2":r_ring}, + sample_points=sample_points,FDTDBuild=FDTDBuild,LumericalPATH=LumericalPATH,runFDTD=runFDTD, + GPUOn=GPUOn, + port_names = ["a1","b1","a2","b2"],) + +class DEVICE_COUPLER(DEVICE_PORTS): + + def __init__(self, dev_name, device, simu_xs="strip", port_width=3, path=None, + wl=[1.5, 1.6], mesh_order=5, layer_heights=[0.22], FDTD_height=2, + material="Si (Silicon) - Palik", CladMaterial="SiO2 (Glass) - Palik", + modeIdx=[1, 2, 3, 4], + sample_points=101, + sourceMode = 1, + FDTDBuild = False, + LumericalPATH = None, + runFDTD = False, + GPUOn = True, + ) -> None: + + super().__init__(dev_name, device, simu_xs, port_width, path, wl, mesh_order, + layer_heights, FDTD_height, material, CladMaterial, modeIdx, ports_extend=["a1","a2","b1","b2"],SimuBox=None,port_radius={}, + sample_points=sample_points,FDTDBuild=FDTDBuild,LumericalPATH=LumericalPATH,runFDTD=runFDTD,GPUOn=GPUOn, + port_names = ["a1","b1","a2","b2"],sourceMode=sourceMode) + +class EULER_CROW_INTER_CP(DEVICE_PORTS): + def __init__(self, dev_name, device, simu_xs="strip", + port_width=3, path=None, wl=[1.5, 1.6], + mesh_order=5, layer_heights=[0.22], + FDTD_height=2, material="Si (Silicon) - Palik", + CladMaterial="SiO2 (Glass) - Palik", + modeIdx=[1, 2, 3, 4], + SimuBox=None, + sample_points=101, + FDTDBuild = False, + LumericalPATH = None, + GPUOn=True, + runFDTD = False) -> None: + + """ Device MUST Be CROW device """ + """ The pins reconized in here is ra1,ra2,ra3,ra4 and rb1,rb2,rb3,rb4 """ + + newDev = device + newDev.cell.pin['a1'] = device.cell.pin['ra2'] + newDev.cell.pin['a2'] = device.cell.pin['rb2'] + newDev.cell.pin['b1'] = device.cell.pin['ra4'] + newDev.cell.pin['b2'] = device.cell.pin['rb4'] + + port_radius = {"a1":device.R1, "a2": device.R1, "b1":-device.R1, "b2":-device.R1} + + super().__init__(dev_name, newDev, simu_xs, port_width, path, wl, mesh_order, layer_heights, FDTD_height, material, CladMaterial, modeIdx, + ports_extend=[], + SimuBox=SimuBox, port_radius=port_radius, sample_points=sample_points, + FDTDBuild=FDTDBuild,LumericalPATH=LumericalPATH,runFDTD=runFDTD,port_names = ["a1","b1","a2","b2"], + GPUOn=GPUOn) + +class EULER_CROW_BUS(DEVICE_PORTS): + def __init__(self, dev_name, device, simu_xs="strip", + port_width=3, path=None, wl=[1.5, 1.6], + mesh_order=5, layer_heights=[0.22], + FDTD_height=2, material="Si (Silicon) - Palik", + CladMaterial="SiO2 (Glass) - Palik", + modeIdx=[1, 2, 3, 4], + SimuBox=None, + sample_points=101, + FDTDBuild = False, + LumericalPATH = None, + GPUOn = True, + runFDTD = False) -> None: + + """ Device MUST Be CROW device """ + """ The pins reconized in here is ra1,ra2,ra3,ra4 and rb1,rb2,rb3,rb4 """ + + newDev = device + # newDev.cell.pin['a1'] = device.cell.pin['ra2'] + newDev.cell.pin['a2'] = device.cell.pin['ra2'] + # newDev.cell.pin['b1'] = device.cell.pin['ra4'] + newDev.cell.pin['b2'] = device.cell.pin['ra4'] + + port_radius = {"a2": device.R1, "b2":-device.R1} + + if (SimuBox is None): + yMax = newDev.cell.pin['b2'].y + yMin = newDev.cell.pin['b1'].y - newDev.ring_cell[0].sz[1]/2 - port_width/2 + + SimuBox = {} + SimuBox["dy"] = yMax - yMin + SimuBox["y"] = (yMax+yMin)/2 + + super().__init__(dev_name, newDev, simu_xs, port_width, path, wl, mesh_order, layer_heights, FDTD_height, material, CladMaterial, modeIdx, + ports_extend=["a1","b1"], + SimuBox=SimuBox, port_radius=port_radius, sample_points=sample_points, + FDTDBuild=FDTDBuild,LumericalPATH=LumericalPATH,runFDTD=runFDTD,port_names = ["a1","b1","a2","b2"], + GPUOn=GPUOn) + + +class RESONATOR(DEVICE_PORTS): + def __init__(self, dev_name, device, + simu_xs="strip", port_width=3, + path=None, wl=[1.5, 1.6], + mesh_order=5, layer_heights=[0.22], + FDTD_height=2, material="Si (Silicon) - Palik", CladMaterial="SiO2 (Glass) - Palik", + modeIdx=[1, 2, 3, 4], + ports_extend=["a1"], + sample_points=10001, + SimuBox=None, + FDTDBuild = False, + LumericalPATH = None, + runFDTD = False) -> None: + super().__init__(dev_name, device, simu_xs, port_width, path, wl, mesh_order, layer_heights, FDTD_height, material, CladMaterial, + modeIdx, ports_extend=["a1","a2","b1","b2"], SimuBox=SimuBox, port_radius=[], sample_points=sample_points, + FDTDBuild=FDTDBuild,LumericalPATH=LumericalPATH,runFDTD=runFDTD,port_names = ["a1","b1","a2","b2"],) + +class RING_PHASE(DEVICE_PORTS): + def __init__(self, dev_name, device, simu_xs="strip", + port_width=3, path=None, wl=[1.5, 1.6], + mesh_order=5, layer_heights=[0.22], + FDTD_height=2, + material="Si (Silicon) - Palik", + CladMaterial="SiO2 (Glass) - Palik", + modeIdx=[1, 2, 3, 4], + SimuBox=None, + port_radius={ "a1": 0 }, + sample_points=101, + FDTDBuild=False, + LumericalPATH=None, + runFDTD=False, + GPUOn = True, + ) -> None: + + + if (hasattr(device,"cell")): + dev_cell = device.cell + + + elif (isinstance(device,nd.Cell)): + dev_cell = device + + dev_cell.pin['a1'] = dev_cell.pin['r1'] + dev_cell.pin['b1'] = dev_cell.pin['r3'] + + + dy = abs(dev_cell.pin['a1'].y - dev_cell.pin['b1'].y ) + + cy = (dev_cell.pin['a1'].y + dev_cell.pin['b1'].y)/2 + if (SimuBox is None): + SimuBox = {} + SimuBox["dy"] = dy + SimuBox["y"] = cy + + if (hasattr(device,"sz")): + SimuBox["dx"] = device.sz[0]/2 + SimuBox["x"] = -device.sz[0]/4 + # SimuBox["dx"] = dx + + + super().__init__(dev_name, dev_cell, simu_xs, port_width, path, wl, mesh_order, layer_heights, + FDTD_height, material, CladMaterial, modeIdx, ports_extend=[], + SimuBox=SimuBox, port_radius=port_radius, sample_points=sample_points, + FDTDBuild=FDTDBuild, LumericalPATH=LumericalPATH, runFDTD=runFDTD, port_names=["a1","b1"], + GPUOn = GPUOn) + + diff --git a/simulation/Lumerical/GDS_SIMU_DEVICE_2X2.lsf b/simulation/Lumerical/GDS_SIMU_DEVICE_2X2.lsf new file mode 100644 index 0000000..72cd5ed --- /dev/null +++ b/simulation/Lumerical/GDS_SIMU_DEVICE_2X2.lsf @@ -0,0 +1,295 @@ +function FUNC_DEVICE_2X2(instPath) +{ + inst = read(instPath); + eval(inst); + + jsonload(jsonPath); + + HwaferMax = 0; + for(idx=1;idx<=length(layers.numbers);idx=idx+1){ + + if (iscell(layers.numbers)){ + layer = num2str(layers.numbers{idx}) + ":" + num2str(layers.datatype{idx});} + + else{ + layer = num2str(layers.numbers(idx)) + ":" + num2str(layers.datatype(idx));} + + if (iscell(layers.heights)){ + Hwafer = layers.heights{idx}*1e-6;} + else{ + Hwafer = layers.heights(idx)*1e-6;} + + + gdsimport(gdspath, devName, layer,wafer.material, 0, Hwafer); + #gdsimport(folder + devName + ".gds", devName, layer,wafer.material, 0, Hwafer); + #set("name","device"); + + if (Hwafer>HwaferMax) + { + HwaferMax = Hwafer; + } + } + + ## output ports + + # HwaferMax = Hwafer; + zOffset = HwaferMax/2*1e+6; + mx_simu_area('FDTD',[FDTD.x,FDTD.y,FDTD.z+zOffset]*1e-6, + [FDTD.dx,FDTD.dy,FDTD.dz]*1e-6,FDTD.mesh_order,[40,40,40]*1e-6,'PML',10e+5); + + portList = ports.names; + + for (idx=1;idx<=length(portList);idx=idx+1) + { + portName = portList{idx}; + + eval("port_struct = ports."+portName+";"); + + if (port_struct.a < 0) { port_struct.a = port_struct.a + 360; } + + if ( 45<=port_struct.a and port_struct.a<135 ) + { + portsType = -2; + portsSZ = [port_struct.width,0,port_struct.height]*1e-6; + } + else if (225<=port_struct.a and port_struct.a<315) + { + portsType = 2; + portsSZ = [port_struct.width,0,port_struct.height]*1e-6; + } + + else if (135<=port_struct.a and port_struct.a<225) + { + portsType = 1; + portsSZ = [0,port_struct.width,port_struct.height]*1e-6; + } + + else if (port_struct.a<45 or port_struct.a>=315) + { + portsType =-1; + portsSZ = [0,port_struct.width,port_struct.height]*1e-6; + } + + + mx_power_monitor(portName,[port_struct.x,port_struct.y,zOffset]*1e-6, + portsSZ,abs(portsType)) ; + set("override global monitor settings",1); + set("frequency points",FDTD.Trans_sample_points); + + mx_mode_expansion(portName+'_modes',[port_struct.x,port_struct.y,zOffset]*1e-6, + portsSZ,abs(portsType),0, + [abs(port_struct.radius)*1e-6, + sign(port_struct.radius)*sign(portsType)*180],1,FDTD.wl*1e-6,portName); + + + if (iscell(modes)){modes = [modes{1}];} + + set("selected mode numbers",modes); + + if ( abs(port_struct.a) <= 45 or abs(port_struct.a) > 315) + { + set('theta',port_struct.a); + set('phi',0); + } + + else if ( abs(port_struct.a) > 135 and abs(port_struct.a) <= 225) + { + set('theta',port_struct.a); + set('phi',0); + } + + else{ + set('theta',90-port_struct.a); + set('phi',90); + } + + + } + + mx_power_monitor('z1',[mont.z1.x,mont.z1.y,zOffset]*1e-6, + [mont.z1.dx,mont.z1.dy,0]*1e-6,3) ; + set("override global monitor settings",1); + set("frequency points",FDTD.Field_sample_points); + + ## adding input ports + + input_struct = ports.a1; + if (input_struct.a < 0) { input_struct.a = input_struct.a + 360; } + + if ( 45<=input_struct.a and input_struct.a<135 ) + { + inputType = -2; + portsSZ = [input_struct.width,0,input_struct.height]*1e-6; + input_theta = 90-input_struct.a; + input_phi = 0; + } + else if (225<=input_struct.a and input_struct.a<315) + { + inputType = 2; + portsSZ = [input_struct.width,0,input_struct.height]*1e-6; + input_theta = 270-input_struct.a; + input_phi = 0; + } + + else if (135<=input_struct.a and input_struct.a<225) + { + inputType = -1; + portsSZ = [0,input_struct.width,input_struct.height]*1e-6; + input_theta = input_struct.a - 180; + input_phi = 0; + + } + + else if (input_struct.a<45 or input_struct.a>=315) + { + inputType = 1; + portsSZ = [0,input_struct.width,input_struct.height]*1e-6; + input_theta = input_struct.a; + input_phi = 0; + } + + mx_mode_source('a1_input',[input_struct.x,input_struct.y,zOffset]*1e-6, + portsSZ,inputType,0, + [abs(input_struct.radius)*1e-6,sign(input_struct.radius)*sign(inputType)*90],FDTD.sourceMode,FDTD.wl*1e-6); + set("theta",input_theta); + set("phi",input_phi); + + ## setting FDTD configurations + + select("FDTD"); + + set("background material",clad.material); + + if (FDTD.GPUOn == 1){ + + setnamed("FDTD", "express mode", true); + setresource("FDTD","GPU", true); + setresource("FDTD", 1, "GPU Device", "Auto"); + set('z min bc','PML'); + } + else { + setnamed("FDTD", "express mode", false); + setresource("FDTD","GPU", false); + + if (length(layers.numbers)==1){ + set('z min bc','symmetric');} + else {set('z min bc','PML');} + + } + + save(folder+"\\"+devName+"_simu.fsp"); + #save("DEVICE_2X2.fsp"); + +} + +function DATA_RETRIEVE_DEVICE_2X2(instPath) +{ + + inst = read(instPath); +eval(inst); +jsonload(jsonPath); + +portList = ports.names; + +save_cmd = ""; + +for (idx=1;idx<=length(portList);idx=idx+1) + +{ + portData = struct; + portData.name = portList{idx}; + + portData.power = getresult(portData.name,"T"); + portData.modes = getresult(portData.name+"_modes","expansion for input"); + E = getresult(portData.name,"E"); + H = getresult(portData.name,"H"); + + x = E.x; + y = E.y; + z = E.z; + + cx = floor(length(x)/2)+1; + cy = floor(length(y)/2)+1; + cz = floor(length(z)/2)+1; + + ## wavelength length + sz_wl = size(E.E,4); + + step = floor((sz_wl-1)/(FDTD.Field_sample_points-1)); + idxSect = [1:step:sz_wl+1]; + + E_save = E; + H_save = H; + + E_save = matrixdataset; + E_save.addparameter("x",E.x); + E_save.addparameter("y",E.y); + E_save.addparameter("z",E.z); + E_save.addparameter("lambda",E.lambda(1:step:end)); + E_save.addattribute("E",E.E(:,:,:,1:step:end,:)); + + H_save = matrixdataset; + H_save.addparameter("x",H.x); + H_save.addparameter("y",H.y); + H_save.addparameter("z",H.z); + H_save.addparameter("lambda",H.lambda(idxSect)); + H_save.addattribute("H",H.H(:,:,:,idxSect,:)); + + portData.E = E_save; + portData.H = H_save; + + ## Center electric field of the monitor + Ecenter = E.E(cx,cy,cz,:,:); ## (x,y,z,wl,Ex/Ex/Ez) + Hcenter = H.H(cx,cy,cz,:,:); ## (x,y,z,wl,Ex/Ex/Ez) + portData.Ecenter = Ecenter; + portData.Hcenter = Hcenter; + + eval(portData.name + " = portData;"); + save_cmd = save_cmd + portData.name + ","; + + } + +if (find(portList=="a1") and find(portList=="b1")){ + Ephase_11 = unwrap(angle(b1.Ecenter) - angle(a1.Ecenter)); + Hphase_11 = unwrap(angle(b1.Hcenter) - angle(a1.Hcenter)); + save_cmd = save_cmd + "Ephase_11" + "," + "Hphase_11" + ","; + + } + +if (find(portList=="a2") and find(portList=="b2")){ + Ephase_22 = unwrap(angle(b2.Ecenter) - angle(a2.Ecenter)); + Hphase_22 = unwrap(angle(b2.Hcenter) - angle(a2.Hcenter)); + save_cmd = save_cmd + "Ephase_22" + "," + "Hphase_22" + ","; + } + +if (find(portList=="a1") and find(portList=="b2")){ + Ephase_12 = unwrap(angle(b2.Ecenter) - angle(a1.Ecenter)); + Hphase_12 = unwrap(angle(b2.Hcenter) - angle(a1.Hcenter)); + save_cmd = save_cmd + "Ephase_12" + "," + "Hphase_12" + ","; + } + +if (find(portList=="a2") and find(portList=="b1")){ + Ephase_21 = unwrap(angle(b1.Ecenter) - angle(a2.Ecenter)); + Hphase_21 = unwrap(angle(b1.Hcenter) - angle(a2.Hcenter)); + save_cmd = save_cmd + "Ephase_21" + "," + "Hphase_21" + ","; + } + +z1= struct; +z1.E = getresult("z1","E"); +z1.H = getresult("z1","H"); +Ex = getresult("z1","Ex"); +Ey = getresult("z1","Ey"); +Ez = getresult("z1","Ez"); + +savefname = devName+"_results.mat"; + +#matlabsave(savefname,b1,b2,z1); + +eval("matlabsave(savefname,"+save_cmd+"z1);"); + + + } + + + + diff --git a/simulation/Lumerical/mx_analysis_lib.lsf b/simulation/Lumerical/mx_analysis_lib.lsf new file mode 100644 index 0000000..eba35c6 --- /dev/null +++ b/simulation/Lumerical/mx_analysis_lib.lsf @@ -0,0 +1,212 @@ +function mx_sweep(sweeps,run_on){ + + #addsweep(0); + + ## NOTICE : this can only be used in one-dimentional result calculation + ## NOTICE : such as [T .vs. lambda] + + #### @ sweeps: struct #### + #### @ sweeps: [var_names={'width','radius','wavelength'...}] + #### @ sweeps: width = [w1:dw:w2] + #### @ radius = [r1:dr:r2] + + addsweep(0); + sweep_name = 'mx_sweep'; + deletesweep(sweep_name); + setsweep("sweep", "name", sweep_name); + setsweep(sweep_name,"type","Ranges"); + + ins_temp = 'var_temp = sweeps.' + sweeps.var_names{1} + ';'; + eval(ins_temp); + setsweep(sweep_name,"number of points",length(var_temp)); + result_data = zeros(length(var_temp),length(sweeps.result_names)); + + + for (idx_vars=1;idx_vars<=length(sweeps.var_names);idx_vars=idx_vars+1){ + + para = struct; + para.Name = sweeps.var_names{idx_vars}; + para.Type = sweeps.var_types{idx_vars}; + para.Parameter = sweeps.var_select{idx_vars}; + + ins_temp = 'var_temp = sweeps.' + sweeps.var_names{idx_vars} + ';'; + eval(ins_temp); + + para.Start = var_temp(1); + para.Stop = var_temp(end); + addsweepparameter(sweep_name, para); + + } + + for (idx_rsult=1;idx_rsult<=length(sweeps.result_names);idx_rsult=idx_rsult+1){ + + result = struct; + result.Name = sweeps.result_names{idx_rsult}; + result.Result = sweeps.result_select{idx_rsult}; + addsweepresult(sweep_name, result); + + } + + if (run_on) { + + runsweep; + + for (idx_result=1;idx_result<=length(sweeps.result_names);idx_result=idx_result+1){ + + result_data(:,idx_result) = getsweepdata(sweep_name,sweeps.result_names{idx_result}); + } + + } + ################################### RESULT ################################################## + + + sweep_cur = sweeps; + sweep_cur.result = result_data; + + return sweep_cur; + + } + + +##### FWM analysis lib ##### + +##### FWM calculation pack ##### +##### @ freq_pump : pumping frequency, vector of 1*2 +##### @ freq_signal: signal frequency, single scalar +##### @ mode_idx : vector of 1*4, [p,p,l,s] denoted +##### @ mode_pol : vector of 1*4, [p,p,l,s] denoted +##### @ wg: struct of {width, bend radius} +##### @ ONLY operates in FDE !!!! ##### + +function mx_FWM_analysis(freq_pump,freq_signal,mode_idx,mode_pol,wg_bend){ + ##### Select the target modes ##### + freq_idler = freq_pump(1) + freq_pump(2) - freq_signal; + freq_range = [freq_pump,freq_signal,freq_idler]; + mode_neff = zeros(1,4); + + for (idx_mode=1;idx_mode<=4;idx_mode=idx_mode+1){ + + temp = mx_get_mode_data(c/freq_range(idx_mode),mode_pol{idx_mode},mode_idx(idx_mode),wg_bend,{'neff'}); + mode_neff(idx_mode) = temp.neff; + + } + + wavelengths = c/freq_range; + wl_neff = wavelengths/mode_neff; + k_vector = 2*pi/wl_neff; + phase_mismatch = sum(k_vector*[1,1,-1,-1]); + + return phase_mismatch; +} + + +##### Phase Matched Coupler calculation pack ##### +##### NOTICE : this is coupled by TE0/TM0 by default ##### + +function mx_DC_analysis(wafer,width_seed,wl,gap,neff,bend,outer_side,mode_idx,mode_pol){ + + cur_file_name = currentfilename; + #load('Temp_workspace.lms'); + switchtolayout; + #deleteall; + + cladding = wafer.cladding; + + max_itn = 5; + + mesh_grids = [20,20,20]*1e-9; + + ##### Adding two waveguides to the strcuture ##### + #mx_FDE_strip(wafer,width_seed); + + #mx_simu_area('FDE_y',[0,0,0],[7e-6,10e-6,3e-6],2,mesh_grids,'Metal',10000); + + width_cur = width_seed; + + neff_sweep = mx_get_mode_data(wl,mode_pol,0,bend,{'neff'}); + run; + setanalysis('use max index',0); + setanalysis('n',neff_sweep.neff); + findmodes; + + width_step = 0.002e-6; + + for (itn=1;itn<=max_itn;itn=itn+1){ + + + width_delta = [-0.01:0.002:0.01]*1e-6; + width_sweep = width_cur+width_delta; + + radius_sweep = width_sweep/2 + gap + outer_side; + + + #### Adding single sweep #### + if (bend>0){ + sweeps = struct; + sweeps.var_names = {'width','radius'}; + sweeps.radius = radius_sweep; + sweeps.width = width_sweep; + sweeps.var_select = {'::model::waveguide::x span','::model::FDE::bend radius'}; + sweeps.var_types = {'Length','Length'}; + sweeps.result_names = {'TE0_neff'}; + + sweeps.result_select = {'::model::FDE::data::mode'+num2str(mode_idx)+'::neff'}; + sweeps.bound_para_idx = [1,2]; + neff_data = mx_sweep(sweeps,1); + neff_data = abs(neff_data.result); + + mismatch = neff*bend - neff_data*radius_sweep; + + } + else { + sweeps = struct; + sweeps.var_names = {'width'}; + sweeps.width = width_sweep; + sweeps.var_select = {'::model::WG::x span'}; + sweeps.var_types = {'Length'}; + sweeps.result_names = {'TE0_neff'}; + sweeps.result_select = {'::model::FDE::data::mode'+num2str(mode_idx)+'::neff'}; + sweeps.bound_para_idx = [1]; + neff_data = mx_sweep(sweeps,1); + neff_data = abs(neff_data.result); + mismatch = neff - neff_data; + } + + if ((mismatch(1)*mismatch(end)) <=0) { + idx_match = find(abs(mismatch) == min(abs(mismatch))); + itn = max_itn + 1; + } + else { + + delta_match = mismatch(1) - mismatch(end); + cent_match = (mismatch(1) + mismatch(end))/2; + delta_width = (width_delta(1)-width_delta(end))/delta_match*cent_match; + width_cur = width_cur - delta_width; + width_cur = round(width_cur/width_step)*width_step; + print(width_cur); + + } + + + } + + #matched_width = width_sweep(idx_match); + return width_sweep(idx_match); + load(cur_file_name); + + +} + + + + + + + + + + + + + + diff --git a/simulation/Lumerical/mx_devices_lib.lsf b/simulation/Lumerical/mx_devices_lib.lsf new file mode 100644 index 0000000..b31685b --- /dev/null +++ b/simulation/Lumerical/mx_devices_lib.lsf @@ -0,0 +1,274 @@ +function mx_euler_ring(name,coord,euler_para,wafer) +{ + R1 = euler_para.R1; + R0 = euler_para.R0; + w1 = euler_para.w1; + w0 = euler_para.w0; + para = euler_para.para; + + ## setting the ring to a group ## + addstructuregroup; + set('name',name); + set('x',0); + set('y',0); + set('z',0); + + wg = mx_euler_wg2wg(euler_para,pi/2,0,[ 0,0,0],'single','linear',wafer,[0,0]); + delete; + sz = wg.sz; + wg = mx_euler_wg2wg(euler_para,pi/2,0,[0,-sz(2),0],'single','linear',wafer,[0,0]); + addtogroup(name); + wg = mx_euler_wg2wg(euler_para,pi/2,0,[0,-sz(2),0],'single','linear',wafer,[1,0]); + addtogroup(name); + wg = mx_euler_wg2wg(euler_para,pi/2,0,[0, sz(2),0],'single','linear',wafer,[0,1]); + addtogroup(name); + wg = mx_euler_wg2wg(euler_para,pi/2,0,[0, sz(2),0],'single','linear',wafer,[1,1]); + addtogroup(name); + + select(name); + set('x',coord(1)); + set('y',coord(2)); + set('z',coord(3)); + + sz = [sz(1)*2,sz(2)*2]; + + ring = struct; + ring.sz = sz; + ring.w = wg.w; + return ring; +} +## =================================================== ## +## DEVICE: ring coupler with euler bend attached +## =================================================== ## +function mx_euler_racetrack(name,coord,euler_para,dLx,dLy,wafer) +{ + + R1 = euler_para.R1; + R0 = euler_para.R0; + w1 = euler_para.w1; + w0 = euler_para.w0; + para = euler_para.para; + + ## generation of the single racetrack ## + + + ## setting the ring to a group ## + addstructuregroup; + set('name',name); + set('x',0); + set('y',0); + set('z',0); + + wg = mx_euler_wg2wg(euler_para,pi/2,0,[ dLx/2,-dLy/2,0],'dual','linear',wafer,[0,0]); + delete; + sz = wg.sz; + wg = mx_euler_wg2wg(euler_para,pi/2,0,[ dLx/2,-dLy/2-sz(2),0],'dual','linear',wafer,[0,0]); + addtogroup(name); + wg = mx_euler_wg2wg(euler_para,pi/2,0,[-dLx/2,-dLy/2-sz(2),0],'dual','linear',wafer,[1,0]); + addtogroup(name); + wg = mx_euler_wg2wg(euler_para,pi/2,0,[ dLx/2, dLy/2+sz(2),0],'dual','linear',wafer,[0,1]); + addtogroup(name); + wg = mx_euler_wg2wg(euler_para,pi/2,0,[-dLx/2, dLy/2+sz(2),0],'dual','linear',wafer,[1,1]); + + dy = sz(2); + dx = sz(1); + + addtogroup(name); + mx_rect('wg_d',[0,-dLy/2-dy,0],[dLx+1e-9,w0,H_wafer],wafer.Material,1,0); + addtogroup(name); + mx_rect('wg_u',[0, dLy/2+dy,0],[dLx+1e-9,w0,H_wafer],wafer.Material,1,0); + addtogroup(name); + mx_rect('wg_l',[-dLx/2-dx,0,0],[w1,dLy+1e-9,H_wafer],wafer.Material,1,0); + addtogroup(name); + mx_rect('wg_r',[ dLx/2+dx,0,0],[w1,dLy+1e-9,H_wafer],wafer.Material,1,0); + addtogroup(name); + + select(name); + set('x',coord(1)); + set('y',coord(2)); + set('z',coord(3)); + + sz = [dLx+2*dx,dLy+2*dy]; + + racetrack = struct; + racetrack.sz = sz; + racetrack.w = wg.w; + return racetrack; + +} +## =================================================== ## +## DEVICE: ring coupler with euler bend attached +## =================================================== ## +function mx_euler_coupler(name,coord,coupler_para,dLc,dAc,Att,coupler_type,H_wafer,mWafer){ + + w_cp = coupler_para.w_cp; + R_cp = coupler_para.R_cp; + Ratt = coupler_para.Ratt; + Rmin = coupler_para.Rmin; + w_wg = coupler_para.w_wg; + para = coupler_para.para; + + addstructuregroup; + set('name',name); + set('x',0); + set('y',0); + set('z',0); + + if (coupler_type=='straight' or coupler_type=='s' or coupler_type=='DC'){ + mx_rect('wg_coupling',[0,0,0],[dLc,w_cp,H_wafer],mWafer,1,0); + addtogroup(name); + dx_attach = dLc/2; + dy_attach = 0; + Acp = 0; + } + else if (coupler_type=='bend' or coupler_type=='b' or coupler_type=='BDC'){ + if (dAc>0){ + mx_ring('wg_coupling',[0,R_cp,0],H_wafer,[R_cp-w_cp/2,R_cp+w_cp/2],[270-dAc/2/pi*180,270+dAc/2/pi*180],mWafer,1,0); + addtogroup(name); + } + Acp = dAc/2; + dx_attach = R_cp*sin(Acp); + dy_attach = R_cp-R_cp*cos(Acp); + } + else { + Acp = 0; + dx_attach = 0; + dy_attach = 0; + } + + euler_para = struct; + euler_para.R0 = R_cp; + euler_para.R1 = Ratt; + euler_para.w0 = w_cp; + euler_para.w1 = w_cp; + euler_para.para = para; + euler_para.order = coupler_para.order; + euler_para.w_offset = coupler_para.w_offset; + + wg = mx_euler_wg2wg(euler_para,Att,Acp,[dx_attach,dy_attach,0],'single','linear',H_wafer,mWafer,[0,0]); + addtogroup(name); + wg = mx_euler_wg2wg(euler_para,Att,Acp,[-dx_attach,dy_attach,0],'single','linear',H_wafer,mWafer,[1,0]); + addtogroup(name); + + euler_para.R0 = Ratt; + euler_para.R1 = Rmin; + euler_para.R2 = Ratt; + euler_para.w0 = w_cp; + euler_para.w1 = w_wg; + + sz = wg.sz; + + mx_euler_wg2wg(euler_para,-Att-Acp,Att+Acp,[-dx_attach-sz(1),dy_attach+sz(2),0],'dual','linear',H_wafer,mWafer,[1,0]); + addtogroup(name); + wg_attach = mx_euler_wg2wg(euler_para,-Att-Acp,Att+Acp,[ dx_attach+sz(1),dy_attach+sz(2),0],'dual','linear',H_wafer,mWafer,[0,0]); + addtogroup(name); + mx_rect('patch',[0,0,0],[1e-9,w_cp,H_wafer],mWafer,1,0); + addtogroup(name); + sz_attach = wg_attach.sz; + select(name); + set('x',coord(1)); + set('y',coord(2)); + set('z',coord(3)); + + + + cp_sz = [sz_attach(1)*2+dx_attach*2+sz(1)*2,sz_attach(2)+dy_attach+sz(2)]; + return cp_sz; + +} +## =================================================== ## +## DEVICE: ring coupler with circular bend attached +## =================================================== ## +function mx_circular_coupler(name,coord,w_couple,R_couple,dAc,Att,wafer) +{ + addstructuregroup; + set('name',name); + set('x',0); + set('y',0); + set('z',0); + + + mx_ring('wg_coupling',[0,R_couple,0],H_wafer,[R_couple-w_couple/2,R_couple+w_couple/2],[270-dAc/2/pi*180,270+dAc/2/pi*180],wafer.Material,1,0); + addtogroup(name); + + mx_ring('wg_in',[R_couple*sin(dAc/2)*2,R_couple*(1-cos(dAc/2)*2),0],H_wafer,[R_couple-w_couple/2,R_couple+w_couple/2],[90,90+dAc/2/pi*180],wafer.Material,1,0); + addtogroup(name); + + mx_ring('wg_out',[-R_couple*sin(dAc/2)*2,R_couple*(1-cos(dAc/2)*2),0],H_wafer,[R_couple-w_couple/2,R_couple+w_couple/2],[90-dAc/2/pi*180,90],wafer.Material,1,0); + addtogroup(name); + + select(name); + set('x',coord(1)); + set('y',coord(2)); + set('z',coord(3)); + + cp_sz = [R_couple*(sin(dAc/2))*2*2,R_couple*(1-cos(dAc/2))*2]; + + + return cp_sz; + + } + +function mx_std_dc(name,coord,gap,w_cp,L_cp,L_attach,w_wg,R0,A,wafer) +{ + addstructuregroup; + set('name',name); + set('x',0); + set('y',0); + set('z',0); + + Lt = abs(w_wg-w_cp)/tan(5/180*pi); + taper_vtx_x = [0,0,-Lt,-Lt]; + taper_vtx_y = [w_cp/2,-w_cp/2,-w_wg/2,w_wg/2]; + + taper_vtx = [taper_vtx_x;taper_vtx_y]; + + mx_rect('cp_u',[0,w_cp/2+gap/2,0],[L_cp,w_cp,H_wafer],wafer.Material,1,0); + addtogroup(name); + mx_rect('cp_d',[0,-(w_cp/2+gap/2),0],[L_cp,w_cp,H_wafer],wafer.Material,1,0); + addtogroup(name); + + mx_ring('cp_ul',[-L_cp/2,(w_cp/2+gap/2)+R0,0],H_wafer,[R0-w_cp/2,R0+w_cp/2],[270-A,270],wafer.Material,1,0); + addtogroup(name); + mx_ring('cp_ul',[-L_cp/2-sin(A/180*pi)*R0*2,R0+(w_cp/2+gap/2)-cos(A/180*pi)*R0*2,0],H_wafer,[R0-w_cp/2,R0+w_cp/2],[90-A,90],wafer.Material,1,0); + addtogroup(name); + + y_port = R0+(w_cp/2+gap/2)-cos(A/180*pi)*R0*2+R0; + x_port = -L_cp/2-sin(A/180*pi)*R0*2-L_attach/2-Lt; + mx_poly('ul_taper',[-L_cp/2-sin(A/180*pi)*R0*2,y_port,0],taper_vtx,H_wafer,wafer.Material,1,0); + addtogroup(name); + mx_rect('ul_attach',[x_port,y_port,0],[L_attach,w_wg,H_wafer],wafer.Material,1,0); + addtogroup(name); + + mx_ring('cp_dl',[-L_cp/2,-((w_cp/2+gap/2)+R0),0],H_wafer,[R0-w_cp/2,R0+w_cp/2],[90,90+A],wafer.Material,1,0); + addtogroup(name); + mx_ring('cp_dl',[-L_cp/2-sin(A/180*pi)*R0*2,-(R0+(w_cp/2+gap/2)-cos(A/180*pi)*R0*2),0],H_wafer,[R0-w_cp/2,R0+w_cp/2],[270,270+A],wafer.Material,1,0); + addtogroup(name); + mx_poly('dl_taper',[-L_cp/2-sin(A/180*pi)*R0*2,-y_port,0],taper_vtx,H_wafer,wafer.Material,1,0); + addtogroup(name); + mx_rect('dl_attach',[x_port,-y_port,0],[L_attach,w_wg,H_wafer],wafer.Material,1,0); + addtogroup(name); + + taper_vtx(1,:) = -taper_vtx(1,:); + + mx_ring('cp_ur',[L_cp/2,(w_cp/2+gap/2)+R0,0],H_wafer,[R0-w_cp/2,R0+w_cp/2],[270,270+A],wafer.Material,1,0); + addtogroup(name); + mx_ring('cp_ur',[L_cp/2+sin(A/180*pi)*R0*2,R0+(w_cp/2+gap/2)-cos(A/180*pi)*R0*2,0],H_wafer,[R0-w_cp/2,R0+w_cp/2],[90,90+A],wafer.Material,1,0); + addtogroup(name); + mx_poly('ur_taper',[-(-L_cp/2-sin(A/180*pi)*R0*2),y_port,0],taper_vtx,H_wafer,wafer.Material,1,0); + addtogroup(name); + mx_rect('ur_attach',[-x_port,y_port,0],[L_attach,w_wg,H_wafer],wafer.Material,1,0); + addtogroup(name); + + mx_ring('cp_dr',[L_cp/2,-((w_cp/2+gap/2)+R0),0],H_wafer,[R0-w_cp/2,R0+w_cp/2],[90-A,90],wafer.Material,1,0); + addtogroup(name); + mx_ring('cp_dr',[L_cp/2+sin(A/180*pi)*R0*2,-(R0+(w_cp/2+gap/2)-cos(A/180*pi)*R0*2),0],H_wafer,[R0-w_cp/2,R0+w_cp/2],[270-A,270],wafer.Material,1,0); + addtogroup(name); + mx_poly('dr_taper',[-(-L_cp/2-sin(A/180*pi)*R0*2),-y_port,0],taper_vtx,H_wafer,wafer.Material,1,0); + addtogroup(name); + mx_rect('dr_attach',[-x_port,-y_port,0],[L_attach,w_wg,H_wafer],wafer.Material,1,0); + addtogroup(name); + + sz = [abs(x_port-L_attach/2)*2,y_port*2]; + return sz; +} \ No newline at end of file diff --git a/simulation/Lumerical/mx_frames_lib.lsf b/simulation/Lumerical/mx_frames_lib.lsf new file mode 100644 index 0000000..3ce8d9d --- /dev/null +++ b/simulation/Lumerical/mx_frames_lib.lsf @@ -0,0 +1,56 @@ +##### FUNCTION LIB ##### +function mx_FDE_dual_strip(wafer,w1,w2,gap){ + + cladding = wafer.cladding; + mx_rect('WG_inner',[-w1/2-gap/2,0,0],[w1,10e-6,wafer.Height],wafer.Material,1,0); + mx_rect('WG_outer',[w2/2+gap/2,0,0],[w2,10e-6,wafer.Height],wafer.Material,1,0); + + if (wafer.clad_on){ + mx_rect('SiO2',[0,0,0],[cladding.Size*2+w1+w2+gap,cladding.Size,cladding.Height],cladding.Material,2,0); + } + else { + mx_rect('SiO2',[0,0,-cladding.Height/2-wafer.Height/2],[cladding.Size*2+w1+w2+gap,cladding.Size,cladding.Height],cladding.Material,2,0); + } + + #inner = wl+gap; + #outer = w2+gap; + #return [inner,outer]; + } + +function mx_FDE_strip(wafer,w){ + + cladding = wafer.cladding; + mx_rect('WG',[0,0,0],[w,10e-6,wafer.Height],wafer.Material,1,0); + if (wafer.slab>=1e-10) { + mx_rect('slab',[0,0,-wafer.Height/2+wafer.slab/2],[cladding.Size*2+w,cladding.Size,wafer.slab],wafer.Material,1,0); + } + if (wafer.clad_on){ + mx_rect('SiO2',[0,0,0],[cladding.Size*2+w,cladding.Size,cladding.Height],cladding.Material,2,0); + } + else { + mx_rect('SiO2',[0,0,-cladding.Height/2-wafer.Height/2],[cladding.Size*2+w,cladding.Size,cladding.Height],cladding.Material,2,0); + } + + #inner = wl+gap; + #outer = w2+gap; + #return [inner,outer]; + } + +##### FUNCTION LIB ##### +function mx_FDE_Disk(wafer,R){ + + cladding = wafer.cladding; + mx_rect('Disk',[-R,0,0],[R*2,10e-6,wafer.Height],wafer.Material,1,0); + + if (wafer.clad_on){ + mx_rect('SiO2',[0,0,0],[cladding.Size*2+R,cladding.Size,cladding.Height],cladding.Material,2,0); + } + else { + mx_rect('SiO2',[0,0,-cladding.Height/2-wafer.Height/2],[cladding.Size*2+R,cladding.Size,cladding.Height],cladding.Material,2,0); + } + + #inner = wl+gap; + #outer = w2+gap; + #return [inner,outer]; + } + \ No newline at end of file diff --git a/simulation/Lumerical/mx_function_lib.lsf b/simulation/Lumerical/mx_function_lib.lsf new file mode 100644 index 0000000..7466da0 --- /dev/null +++ b/simulation/Lumerical/mx_function_lib.lsf @@ -0,0 +1,105 @@ +##### maxwell lib ##### +##### Function lib #### + +function get_system_time() +{ + fname="cur_time.txt"; # file name to store current time + cmd="echo %date:~0,4%%date:~5,2%%date:~8,2%_%time:~0,2%%time:~3,2%%time:~6,2%> "+fname; # system command to get current time and write to fname + system(cmd); # run command to get time and save to file + cur_time=read(fname); # read time from file + cur_time = substring(cur_time,1,15); + + return cur_time; +} + +function mode_polar_select(polar_name,current_mode_name){ + +if (polar_name=='TE') { + polar_select = 0.7; + } + +else if (polar_name=='TM') { + polar_select = -0.3; + } + + cur_pol = getresult(current_mode_name,'TE polarization fraction'); + + selected = ((cur_pol*sign(polar_select))>polar_select); + + + return selected; + } + +function mx_get_sys_time(){ + system("notepad"); + fname="cur_time.txt"; # file name to store current time + cmd="echo %time% "+fname; # system command to get current time and write to fname + rm(fname); # delete time file + system(cmd); # run command to get time and save to file + + cur_time=read(fname); # read time from file + return cur_time; +} + +#### @ result : {'neff','loss'} cell arrays ##### +function mx_get_mode_data(center_wl,mode_pol,mode_idx,wg_bend,results){ + + run; + setanalysis('use max index',1); + setanalysis('number of trial modes',20); + setanalysis('wavelength',center_wl); + if (wg_bend==0) { + setanalysis('bent waveguide',0); + + } + else if (wg_bend>0) { + setanalysis('bent waveguide',1); + setanalysis('bend radius',wg_bend); + setanalysis('bend orientation',90); + + + } + + else { + setanalysis('bent waveguide',1); + setanalysis('bend radius',abs(wg_bend)); + setanalysis('bend orientation',-90); + + + } + + n_modes = findmodes; + idx_TE = 0; + idx_TM = 0; + idx_final = 0; + for (idx_n=1;idx_n<=n_modes;idx_n=idx_n+1){ + cur_mode_name = 'FDE::data::mode'+num2str(idx_n); + + if (mode_polar_select('TE',cur_mode_name)){ + idx_TE = idx_TE+1; + if (idx_TE==(mode_idx+1) & (mode_pol == 'TE')){ + idx_final = idx_n; + idx_n = n_modes +1 ; + } + } + if (mode_polar_select('TM',cur_mode_name)){ + idx_TM = idx_TM+1; + if (idx_TM==(mode_idx+1) & (mode_pol == 'TM')){ + idx_final = idx_n; + idx_n = n_modes +1 ; + } + } + } + + final_data = struct; + final_mode_name = 'FDE::data::mode'+num2str(idx_final); + + for (idx_result=1;idx_result<=length(results);idx_result = idx_result+1){ + temp = getdata(final_mode_name,results{idx_result}); + ins = 'final_data.'+results{idx_result}+ '=temp;'; + eval(ins); + } + + return final_data; + +} \ No newline at end of file diff --git a/simulation/Lumerical/mx_lib_install.lsf b/simulation/Lumerical/mx_lib_install.lsf new file mode 100644 index 0000000..968c9c3 --- /dev/null +++ b/simulation/Lumerical/mx_lib_install.lsf @@ -0,0 +1,27 @@ +#### Lib Install #### + +PATH_LIB = ABSOLUTE_LIB_DIR; +NAME_LIB_1 = 'mx_structures_lib.lsf'; +feval(PATH_LIB+NAME_LIB_1); + +NAME_LIB_2 = 'mx_simulation_lib.lsf'; +feval(PATH_LIB+NAME_LIB_2); + +NAME_LIB_3 = 'mx_function_lib.lsf'; +feval(PATH_LIB+NAME_LIB_3); + +NAME_LIB_4 = 'mx_analysis_lib.lsf'; +feval(PATH_LIB+NAME_LIB_4); + +NAME_LIB_5 = 'mx_frames_lib.lsf'; +feval(PATH_LIB+NAME_LIB_5); + +NAME_LIB_6 = 'mx_poly_spiral_lib.lsf'; +feval(PATH_LIB+NAME_LIB_6); + +NAME_LIB_7 = 'mx_devices_lib.lsf'; +feval(PATH_LIB+NAME_LIB_7); + + +NAME_LIB = 'GDS_SIMU_DEVICE_2X2.lsf'; +feval(PATH_LIB+NAME_LIB); \ No newline at end of file diff --git a/simulation/Lumerical/mx_poly_spiral_lib.lsf b/simulation/Lumerical/mx_poly_spiral_lib.lsf new file mode 100644 index 0000000..fb7eaa7 --- /dev/null +++ b/simulation/Lumerical/mx_poly_spiral_lib.lsf @@ -0,0 +1,265 @@ +## Generate a spiral line ## +function mx_poly_spiral(r,theta,coord,order,para){ + #UNTITLED2 Summary of this function goes here + # Detailed explanation goes here + + dL = para.dL; + r_init = r(1); + r_end = r(2); + theta_init = theta(1); + theta_end = theta(2); + x_init = coord(1); + y_init = coord(2); + + K0 = 1/r_init; + K1 = 1/r_end; + + L0 = abs(theta_end - theta_init)/(K0+(K1-K0)*order/(order+1)); + + L = [0:dL:L0]; + K = K0 + (K1 - K0)/L0^order * (L0^order - abs(L-L0)^order); + + R = 1/K; + R = (R<=para.R_max)*R + (R>para.R_max)*para.R_max*ones(length(R),1); + direction = sign(theta_end-theta_init); + dt = direction*dL/R; + #theta_temp = cumsum(dt)+theta_init; + theta_temp = dt; + + x=zeros(length(L),1)+x_init; + y=zeros(length(L),1)+y_init; + + + + for (i=2;i<=length(L);i=i+1){ + theta_temp(i) = theta_temp(i)+theta_temp(i-1); + cur_theta = theta_temp(i)+theta_init; + pre_theta = theta_temp(i-1)+theta_init; + x(i) = x(i-1) + direction* R(i)*( sin( cur_theta ) - sin(pre_theta ) ); + y(i) = y(i-1) - direction* R(i)*( cos( cur_theta ) - cos( pre_theta ) ); + + } + + theta_temp = [theta_temp(1);theta_temp(2:50:end-1);theta_temp(end)]+theta_init; + x = [x(1);x(2:50:end-1);x(end)]; + y = [y(1);y(2:50:end-1);y(end)]; + + + vtx = [x,y,theta_temp]; + + return vtx; + +} + +function mx_wg_draw(vtx,width){ + #UNTITLED6 Summary of this function goes here + # Detailed explanation goes here + + z = vtx(:,1) + 1i*vtx(:,2); # complex points + #dz = diff(z); # direction of each point + dz = z(2:end) - z(1:end-1); + + dz = [transpose(dz),dz(end)]; + dir_upper = -1i*real(dz)+imag(dz); + dir_down = 1i*real(dz)-imag(dz); + + p_upper = [z + dir_upper*width/2/abs(dir_upper)]; + p_down = [z+ dir_down*width/2/abs(dir_down)]; + + wg = struct; + wg.curve_inner = [real(p_upper),imag(p_upper)]; + wg.curve_outer = [real(p_down),imag(p_down)]; + + return wg; + +} + +function mx_euler_wg(vtx,width,offset){ + #UNTITLED6 Summary of this function goes here + # Detailed explanation goes here + + z = vtx(:,1) + 1i*vtx(:,2); # complex points + #dz = diff(z); # direction of each point + dz = sin(vtx(:,3))*1i + cos(vtx(:,3)); + + dz = [transpose(dz)]; + dir_upper = -1i*real(dz)+imag(dz); + dir_down = 1i*real(dz)-imag(dz); + + p_upper = [z + dir_upper*(offset+width/2)/abs(dir_upper)]; + p_down = [z+ dir_down*(-offset+width/2)/abs(dir_down)]; + + wg = struct; + wg.curve_inner = [real(p_upper),imag(p_upper)]; + wg.curve_outer = [real(p_down),imag(p_down)]; + + return wg; + +} + + + +function mx_euler_wg2wg(euler_para,bend_angle,theta_start,coord,bend_type,width_type,Height,Material,vtx_flip) +{ + R0 = euler_para.R0; + R1 = euler_para.R1; + Win = euler_para.w0; + dW = euler_para.w1 - euler_para.w0; + order = euler_para.order; + + if (bend_type=='single'){ + vtx_start = mx_poly_spiral([R0,R1],[theta_start,bend_angle+theta_start],[0,0],order,euler_para.para); + p_start = vtx_start(1,:); + p_end = vtx_start(end,:); + vtx_euler_bend = vtx_start; + } + + else { + R2 = euler_para.R2; + + vtx_start = mx_poly_spiral([R0,R1],[theta_start,bend_angle/2+theta_start],[0,0],order,euler_para.para); + vtx_stop = mx_poly_spiral([R1,R2],[bend_angle/2+theta_start,bend_angle+theta_start],[vtx_start(end,1),vtx_start(end,2)],order,euler_para.para); + p_start = vtx_start(1,:); + p_end = vtx_stop(end,:); + vtx_euler_bend = [vtx_start;vtx_stop(2:end,:)] ; + } + + ## attaching waveguide + dx = abs(p_end(2) - p_start(1));## displacement in x direction + dL = (vtx_euler_bend(2:end,2) - vtx_euler_bend(1:end-1,2))^2 + (vtx_euler_bend(2:end,1) - vtx_euler_bend(1:end-1,1))^2; + dL = sqrt(dL); + ##L = cumsum(dL) ## L for each pieces + L = zeros(length(dL),1); + L(1) = dL(1); + for (idx=2;idx<=length(L);idx=idx+1){ + L(idx) = L(idx-1)+dL(idx); + } + L = [0;L]; + + L0 = sum(dL); + w_offset = 0; + if (width_type=='cos'){## in this situation, dW is the difference of input and output + dy = abs(p_end(2) - p_start(2)); ## displacement in y direction + vtx_euler_bend(:,2) = -dy + vtx_euler_bend(:,2); + z = vtx_euler_bend(:,3); + + #z = vtx_euler_bend(:,1) + 1i*vtx_euler_bend(:,2); + w = dW/2*cos(z*pi/abs(bend_angle)) + (Win*2+dW)/2; + } + + else if (width_type=='sin'){ ## in this situation, win = wout, dW is the middle width difference + if (abs(bend_angle-pi)<0.001){ + dy = abs(vtx_start(end,2) - vtx_start(1,2)); ## displacement in y direction + } + else { + dy = abs(p_end(2) - p_start(2)); ## displacement in y direction + } + vtx_euler_bend(:,2) = -dy + vtx_euler_bend(:,2); + z = vtx_euler_bend(:,3); + + w = dW*cos(z+pi/2)^2 + Win; ## revised 2023.03.27 + + vtx_euler_bend(:,2) = dy + vtx_euler_bend(:,2); + + } + + else if (width_type=='pumpkin'){ ## in this situation, win = wout, dW is the middle width difference + if (abs(bend_angle-pi)<0.001){ + dy = abs(vtx_start(end,2) - vtx_start(1,2)); ## displacement in y direction + } + else { + dy = abs(p_end(2) - p_start(2)); ## displacement in y direction + } + vtx_euler_bend(:,2) = -dy + vtx_euler_bend(:,2); + z = vtx_euler_bend(:,3); + + z = z; + z = z^0.5*(pi/2)^0.5; + z = sin(z)^2*pi/2; + + w = dW*sin( z )^2 + Win; ## revised 2023.05.04 + #w = dW*sin( z )^2 + Win; ## revised 2023.05.04 + #w = dW*theta/(pi/2) + Win; ## revised 2023.05.04 + vtx_euler_bend(:,2) = dy + vtx_euler_bend(:,2); + + } + + else if (width_type=='special'){ ## in this situation, win = wout, dW is the middle width difference + if (abs(bend_angle-pi)<0.001){ + dy = abs(vtx_start(end,2) - vtx_start(1,2)); ## displacement in y direction + } + else { + dy = abs(p_end(2) - p_start(2)); ## displacement in y direction + } + vtx_euler_bend(:,2) = -dy + vtx_euler_bend(:,2); + z = vtx_euler_bend(:,3); + + z = z; + z = z^0.65*(pi/2)^0.35; + z = sin(z)^2*pi/2; + + w = dW*sin( z )^2 + Win; ## revised 2023.05.04 + vtx_euler_bend(:,2) = dy + vtx_euler_bend(:,2); + + } + + else if (width_type=='sin2'){ ## in this situation, win = wout, dW is the middle width difference + if (abs(bend_angle-pi)<0.001){ + dy = abs(vtx_start(end,2) - vtx_start(1,2)); ## displacement in y direction + } + else { + dy = abs(p_end(2) - p_start(2)); ## displacement in y direction + } + vtx_euler_bend(:,2) = -dy + vtx_euler_bend(:,2); + z = vtx_euler_bend(:,1) + 1i*vtx_euler_bend(:,2); + w = dW/2*sin(angle(z)*2+abs(bend_angle))^2 + Win+dW/2; + } + + else if (width_type=='linear'){ + w = dW/L0*L + Win;} + + else if (width_type=='dual_linear'){ + w = dW/2/(L0/2)*abs(L-L0/2) + Win+dW/2;} + + else if (width_type=='linear_offset'){ + w = dW/L0*L + Win; + w_offset = euler_para.w_offset; + } + + else { ## default linear from input to output + w = dW/L0*L + Win;} + + + sz = abs([p_end(1) - p_start(1),p_end(2) - p_start(2)]); ## the size of the bending + wg = mx_euler_wg(vtx_euler_bend,w,w_offset); + + vtx = [wg.curve_outer;flip(wg.curve_inner,1)]; + + if (vtx_flip(1)==1){ + vtx(:,1) = -vtx(:,1);} + if (vtx_flip(2)==1){ + vtx(:,2) = -vtx(:,2);} + + mx_poly('euler',coord,vtx,Height,Material,1,0); + wg = struct; + wg.sz = sz; + wg.w = w; + wg.vtx = vtx_euler_bend; ## central line + + z = vtx(:,1) + 1i*vtx(:,2); # complex points + #dz = diff(z); # direction of each point + dz = z(2:end) - z(1:end-1); + + dz = [transpose(dz),dz(end)]; + dir_upper = -1i*real(dz)+imag(dz); + dir_down = 1i*real(dz)-imag(dz); + + wg.angle = -angle(dir_upper); + return wg; + +} + + + + + diff --git a/simulation/Lumerical/mx_simulation_lib.lsf b/simulation/Lumerical/mx_simulation_lib.lsf new file mode 100644 index 0000000..64ef39a --- /dev/null +++ b/simulation/Lumerical/mx_simulation_lib.lsf @@ -0,0 +1,207 @@ +## @ cord = [x,y,z] denotes a rectangle +## @ cord_span = [xs,ys,zs] denotes a rectangle + +function mx_simu_area(name,coord,span,mesh_accuracy,meshs,boundary,time) { + + if((name=='FDTD')){ + addfdtd; + set('simulation time',time*1e-15); + + FDE_on = 0; + } + + else if ((name=='FDE_x')){ + FDE_on = 1; + addfde; + set('solver type','2D X normal'); + span(1) = 0; + } + + else if (name=='FDE_y'){ + addfde; + FDE_on = 1; + set('solver type','2D Y normal'); + span(2) = 0; + } + + else if (name=='FDE_z'){ + addfde; + FDE_on = 1; + set('solver type','2D Z normal'); + span(3) = 0; + } + + + + + x = coord(1); + y = coord(2); + z = coord(3); + + set('x',x); + set('y',y); + set('z',z); + + xs = span(1); + ys = span(2); + zs = span(3); + + if (xs>0) { + set('x span',xs); + set('x min bc',boundary); + set('x max bc',boundary); + if (FDE_on) { + set('define x mesh by','maximum mesh step'); + set('dx',meshs(1)); + } + } + + if (ys>0) { + set('y span',ys); + set('y min bc',boundary); + set('y max bc',boundary); + if (FDE_on) { + set('define y mesh by','maximum mesh step'); + set('dz',meshs(2)); + } + } + + if (zs>0) { + set('z span',zs); + set('z min bc',boundary); + set('z max bc',boundary); + if (FDE_on) { + set('define z mesh by','maximum mesh step'); + set('dz',meshs(3)); + } + } + if((name=='FDTD')){ + if (mesh_accuracy==0){ + addmesh; + set('dx',meshs(1)); + set('dy',meshs(2)); + set('dz',meshs(3)); + set('x',x); + set('y',y); + set('z',z); + + xs = span(1); + ys = span(2); + zs = span(3); + + set('x span',xs); + set('y span',ys); + set('z span',zs); + + } + else { + set('mesh accuracy',mesh_accuracy); + } + } + +} + +## +function mx_mode_source(name,coord,span,s_dir,theta,bend,mode_idx,wl){ + addmode; + set('name',name); + set('injection axis',abs(s_dir)); + set('direction',1.5-sign(s_dir)/2); + + x = coord(1); + y = coord(2); + z = coord(3); + + set('x',x); + set('y',y); + set('z',z); + + xs = span(1)*((abs(s_dir))~=1); + ys = span(2)*((abs(s_dir))~=2); + zs = span(3)*((abs(s_dir))~=3); + + if (xs>0) { + set('x span',xs); + } + + if (ys>0) { + set('y span',ys); + } + + if (zs>0) { + set('z span',zs); + } + + if (bend(1)>0){ + set('bent waveguide',1); + set('bend orientation',bend(2)); + set('bend radius',bend(1)); + + } + else{ + set('bent waveguide',0); + } + + set('theta',theta); + set('theta',theta); + set('mode selection','user select'); + set('wavelength start',wl(1)); + set('wavelength stop',wl(2)); + + updatesourcemode(mode_idx); + + +} + +## +function mx_mode_expansion(name,coord,span,m_dir,theta,bend,mode_idx,wl,monitor_name){ + addmodeexpansion; + set('name',name); + set('monitor type',abs(m_dir)); + + x = coord(1); + y = coord(2); + z = coord(3); + + set('x',x); + set('y',y); + set('z',z); + + xs = span(1)*((abs(m_dir))~=1); + ys = span(2)*((abs(m_dir))~=2); + zs = span(3)*((abs(m_dir))~=3); + + if (xs>0) { + set('x span',xs); + } + + if (ys>0) { + set('y span',ys); + } + + if (zs>0) { + set('z span',zs); + } + + if (bend(1)>0){ + set('bent waveguide',1); + set('bend orientation',bend(2)); + set('bend radius',bend(1)); + + } + else{ + set('bent waveguide',0); + } + set('theta',theta); + set('mode selection','user select'); + + updatemodes(mode_idx); + + setexpansion('input',monitor_name); + + +} + + + + diff --git a/simulation/Lumerical/mx_structures_lib.lsf b/simulation/Lumerical/mx_structures_lib.lsf new file mode 100644 index 0000000..520fa18 --- /dev/null +++ b/simulation/Lumerical/mx_structures_lib.lsf @@ -0,0 +1,205 @@ +function mx_rect(name,coord,sz,material,mesh_order,refractive_index) +{ + addrect; + set('name',name); + set('x',coord(1)); + set('y',coord(2)); + set('z',coord(3)); + set('x span',sz(1)); + set('y span',sz(2)); + set('z span',sz(3)); + + set('material',material); + + set('override mesh order from material database',1); + set('mesh order',mesh_order); + + if (refractive_index>0) + { + set('material',''); + set('index',refractive_index); + } +} + +function mx_concoid(name,coord,height,R0,T0,kR,w0,res,theta,material,mesh_order,refractive_index) +{ + ## in polar axis + + dT = linspace(T0,theta+T0,res); + R = ((dT-T0)*kR+R0); + + e_theta = -1/((R0/kR)+dT-T0); + e_rou = ones(length(dT)); + e_theta(end) = 0; + e_theta(1) = 0 ; + + + + ex = cos(dT)*e_rou - sin(dT)*e_theta; + ey = sin(dT)*e_rou + cos(dT)*e_theta ; + + Lnorm = sqrt(ex^2+ey^2); + + vtx_x = R*cos(dT); + vtx_y = R*sin(dT); + + vtx_out_x = vtx_x + w0/2*ex/Lnorm; + vtx_out_y = vtx_y + w0/2*ey/Lnorm; + + vtx_in_x = vtx_x - w0/2*ex/Lnorm; + vtx_in_y = vtx_y - w0/2*ey/Lnorm ; + + vtx_in = [flip(vtx_in_x,1),flip(vtx_in_y,1)]; + vtx_out = [vtx_out_x,vtx_out_y]; + vtx = [vtx_out;vtx_in]; + + + mx_poly(name,coord,vtx,height,material,mesh_order,refractive_index); + + return vtx; + + } + +function mx_taper(name,coord,height,wa,wb,L,offset,material,mesh_order,refractive_index) +{ + vtx_x = [0,L,L,0]; + vtx_y = [wa/2,wb/2+offset,-wb/2+offset,-wa/2]; + + vtx = [vtx_x;vtx_y]; + mx_poly(name,coord,vtx,height,material,mesh_order,refractive_index); + + return vtx; + + + } + +function mx_ring(name,coord,height,radius,theta,material,mesh_order,refractive_index) +{ + addring; + set('name',name); + set('x',coord(1)); + set('y',coord(2)); + set('z',coord(3)); + set('z span',height); + set('outer radius',max(radius)); + set('inner radius',min(radius)); + set('theta start',theta(1)); + set('theta stop',theta(2)); + + set('material',material); + + set('override mesh order from material database',1); + set('mesh order',mesh_order); + + if (refractive_index>0) + { + set('material',''); + set('index',refractive_index); + } + +} + +function mx_ring_coic(name,coord,height,Ra_in,Rb_in,Ra_out,Rb_out,offset,theta,material,mesh_order,refractive_index) +{ + theta = linspace(theta(1),theta(2),1001); + xout = Ra_out*cos(theta); + yout = Rb_out*sin(theta); + + xin = Ra_in*cos(theta); + yin = Rb_in*sin(theta)+offset; + + vtx_outer = [xout,yout]; + vtx_inner = [xin,yin]; + + vtx = [vtx_outer;flip(vtx_inner,1)]; + mx_poly(name,coord,vtx,height,material,mesh_order,refractive_index); + + return vtx; + + } + + +function mx_poly(name,coord,vtx,height,material,mesh_order,refractive_index) +{ + addpoly; + set('name',name); + set('x',coord(1)); + set('y',coord(2)); + set('z',coord(3)); + set('vertices',vtx); + set('z span',height); + + set('material',material); + + set('override mesh order from material database',1); + set('mesh order',mesh_order); + + if (refractive_index>0) + { + set('index',refractive_index); + } + + } + +function mx_elipse(name,coord,height,La,Lb,wa,wb,theta,offset,material,mesh_order,refractive_index) +{ + theta = linspace(theta(1),theta(2),1001); + x = La*cos(theta); + y = Lb*sin(theta); + + ## norm direction + dX = 2*x/La^2; + dY = 2*y/Lb^2; + + L_norm = sqrt(dX^2+dY^2); + + offset = (offset(1)-offset(2))*cos(theta)^2 + offset(2); + + w = (wa-wb)*cos(theta)^2 + wb; ## width variation + + vtx_outer_x = x + dX/L_norm*(w/2+offset); + vtx_outer_y = y + dY/L_norm*(w/2+offset); + + vtx_inner_x = x + dX/L_norm*(-w/2+offset); + vtx_inner_y = y + dY/L_norm*(-w/2+offset); + + vtx_outer = [vtx_outer_x,vtx_outer_y]; + vtx_inner = [vtx_inner_x,vtx_inner_y]; + + vtx = [vtx_outer;flip(vtx_inner,1)]; + mx_poly(name,coord,vtx,height,material,mesh_order,refractive_index); + + return vtx; +} + +function mx_power_monitor(name,coord,sz,diretion) +{ + addpower; + set('name',name); + set('x',coord(1)); + set('y',coord(2)); + set('z',coord(3)); + + if (diretion==1) + { + set('monitor type','2D X-normal'); + set('y span',sz(2)); + set('z span',sz(3)); + } + + if (diretion==2) + { + set('monitor type','2D Y-normal'); + set('x span',sz(1)); + set('z span',sz(3)); + } + + if (diretion==3) + { + set('monitor type','2D Z-normal'); + set('x span',sz(1)); + set('y span',sz(2)); + } +} + + diff --git a/simulation/LumericalPATH.json b/simulation/LumericalPATH.json new file mode 100644 index 0000000..1dbbace --- /dev/null +++ b/simulation/LumericalPATH.json @@ -0,0 +1,13 @@ +{ + "DEFAULT_PATH":[ + + "C:\\Program Files\\Lumerical\\v241\\api\\python\\", + "D:\\Program Files\\Lumerical\\v241\\api\\python\\", + "F:\\Program Files\\Lumerical\\v241\\api\\python\\", + + "C:\\Program Files\\Lumerical\\v232\\api\\python\\", + "D:\\Program Files\\Lumerical\\v232\\api\\python\\", + "F:\\Program Files\\Lumerical\\v232\\api\\python\\" + + ] +} \ No newline at end of file diff --git a/simulation/__init__.py b/simulation/__init__.py new file mode 100644 index 0000000..cf7f42d --- /dev/null +++ b/simulation/__init__.py @@ -0,0 +1,6 @@ +from .DualPortElements import DEVICE_RING_BUS,DEVICE_COUPLER,EULER_CROW_INTER_CP,RESONATOR,RING_PHASE,EULER_CROW_BUS,DEVICE_PORTS + + +from .DualPortElements import DEVICE_2X2_FDTD_INIT,__getLumericalLibPATH__,__checkLumericalDIR__,SimuDataFigurePlot + + diff --git a/simulation/__pycache__/DualPortElements.cpython-310.pyc b/simulation/__pycache__/DualPortElements.cpython-310.pyc new file mode 100644 index 0000000000000000000000000000000000000000..04e9bbc8ff238dfae2f6d941eadaf2cbcb91651d GIT binary patch literal 16024 zcmch8dvIG%@cj@Z>p|#2i`2umL|c|bQqrtMsUo%BMY+8T;eZe*5TFh& zO%grGI-xshk}7taCQWL4hIG?*JWY1mPXA1sX_~fan{ArrpUjzZ(rKSZ-3g)??5`h?_fH}@6d8M-KJ9?MlrG+O-GS7h4rSE zjtQUecsi)&2B*G}^5Tl|^L{QQxq9TTpE|C&((SeeL=oDQj=@tpmEqe0C zq9J-kA8LB8X`){Yys3$SO(Wfl+@KgjZYbBgX{P&-9~L9Xk8r*p`2!+}e3J76$R8Ak zkUzxvLFAtkhmk+b`61+wh@;3K<@~VF=d{%ECxo@xl)X1}d$pRTG<0UnR@)YCt}f?fA)6^Kq5{hV15F0;db%JkEZliecFLtfc40oVk}Iz|3lpa&&YwDc z=FIr1+bd4t`WH^k6^fXy%q(1+np;?2Unykm#j?!j7Al41)ndjeluHX!tC`}RvUH}4 zx#e8Rsf?{`xV^TW&pFLW-Ynd|b9w%{ZMV@1X&QaHrAs2N7DNydXm*ToQk}J6ucWxES8R7b5rj*aQW)bE7+kGdFHd z+gEST&QE`7UJ|lhvmn-a$2}0DvTobXKFDPsHuq=h#$DTPM^z;OkkDg#L?8M%XvXz8 zQU+488x>gtpz)&>(6o847LiaO(p61{mG%y}_&uX$lm;A>I{Ky|4=WI~FoEGArN0A0 z<~AnU1Sv1fdWlstszJxxvXrrn-FprDHphENyGZGjRqae>We&6lv@Vv5A}1@!*{6P+ zNOIw+P7Dz)0O=JcQ_2Ew6-TLgpmdmC6R~EU2W7R23Ufdko z))(h3tG+lxUPjGK%Qn)M(!<;Mn@A{a3x^$vLQb(zMlxj#5g)}DbFXKLqm55f7pomJ_hzr!Q`29(adq0B z_|nAav8stVR4ut$x>y~Ux;%f`zB@gC_r~-Wr(2q(EZ3sWSSZPUOj!;9Q1PO0glnP$ zH$ub9?5aZ}lBUQ?(_XGg>rI*KubR zb+%Ks>#hDWKB|!J(5cHW;5F|WnQ_;+Ki)$72T{{J9$y5{Gi`zPrN0DG=4rq`0yaTty3-1fk|yoi{PFXp_Ue5)VA$u z5pHP_K`lMIM~Sqwh@zIB-7TUxs{64nQhF>AtM}3D2ZwxXO>>qwep_qW8$H)Tdurm% zn(&^Qj%H08YJ3i)gbF|@fm*Q?gb@WZN+cee3q~zci=r$<=9K7tj04{ZY=yVAT1>^3 z+O{I1kEDR;ccSPg=CrHerY;5sVP286KzguGtCu2`H*4)q9QpRTD!x*S(#dj=m;5|L4CCz%;Gt`d+|Pucc1FTY>u$>?Op0a z%Y*nTgy(2NQ{C9RV~=4BF>Up$9;~CkwM}gxuhj-o->&*;haOsQ>+<&7)KDFdzykUjR)bpuYM4fUTYvN~N<>*`)mlS+q4z=TGWB(+HlzmC5avFt zXvX-yF-tLu?+2;BgjO3-eOx;tQflZO?L7mnjz6||q$lD;9eZDI;_>VI;PE8m@kAYu z2e5KVH}H6>j>oaRc--2H#{;C?ekL3mvido2XgtdOf5M^h=x4<9bsT=XiNkLLL(kW- zl>~*IW?CO->f+odC@k67JG7o4Jl6LxxoZzwX>z3bh4Hp2dqU9louRGax;#JQ=@qYO zCG50kxI^MVr(tQ2T7I8})ITTAix-HtA%hNVlX-7>Z`L)*1`aB=1uFrvKU2*-z%<_GaagDX2>sv01TqW(Bw@W#EGVSWuU86d_M-zNP z8feORLe}yVGVAC_TrA`(g=#L3EPmg-lt072-@BADT>l!pF)&8`5BOaTo)|xUnjY^| zrL$392dGBZGOt4yD&?IAP*u3E+&8{)xoVD$FILTDX<6Q)wl@LNq3M+enM%$cA5X!= zqhdat4QV&h)q#vXo^;B|`}Cs4Cja9616f|pKllpw{l=SPmofJc{1rLY@}H1TCe2u| zF>hZnb3a#xt%&6U`xMTnx`cQW{75L}NNsSRJlUA!AQ8rf@e0{UA%;i?1|3 zt>j%3pl}lFA|%4~L>p;}l|m}X3Y2R-I0tu=H3JuyJVC<;DvwrkxoXZ0 zy^bNu5<9{=2W*p~1M6}1{z|b74@!F@WoN|FD(u2u*o|z^!G{HxS>9gGlrnixi+`<9 z66H0wZEYD&t@3itksC5aW07%LmXZt?N!}CJgxAdtfao5;q2^kZM-t`|CX|twPX+xu z1*ZrqSjm=)W$9YbXp8XNVM4CY`^KMHf&T|ipJF*be%f;r4SMSWVPIOWFIy~Ea;bw@ z3mhKglz`VS_xdb680laIjb@4&v?8kuUO|y@GH|g1l7*TMQZYC6P1`FnpZ2d6>Ai=( zU6V#KGV5-jj)im>4dF_%fkjDVpEjn~@oeC!rqS7QCZ7(@6N!Ny)BbtpBOqjboWKO! zzL>*YI_*aX_F56{U*zu;9WJ;b*3*8QO6*1_;z}bmF`f?G%`Fz~iSfO_Er75AOBx=T1bi``ReebJ{toMgF`(m*Va!OFG1iw)8gYHp7}ZDh z!^SA;=-VFnzaprGO%7Y!fOZWl9lgf%PSnM8;|F02+BNFqWSi>=*yuZ7^nIEypsB}< zB*vic{g}0Zp8bF*(<7$*DJ&+7r{zbi@+zLnmjH6?bx*^GB;aEcu;%lI@QXlR7eVrs z=b^a1uZxI?qQoR$c|jNLXd6ck{&LiJzOOHKhfJt;@Rz^ogZJKq*56D1O*VG!+`cYrQ#FjR6@?-_q=vR)e3-?cDjkE0WGuC-;eL3|*tq|yGHW9x|4wOb z1pRcV5h#!$D1439og)0$+={9p?`u{+=g_iCw5dL#0;;8PUW1x#;jig2e`~kL4qqg& z&H=1pP{c$#trpIf12C@+i1;L#SqDT%{oPrAcR}gurFtk538d?(ehn*%*M^;LHR#O> zE*Y2-k4>+I8io<{YbuUbky=v471ZB)4JpDssQn$_k9(^fbb8c5(Tnmvv^q$NDy90! zmwITw8Ym2*L9OnEUROT zzobg4_O0G+jc5N5MhM~Qs7g+f3cdggp+A&Z*fXzZ^xMp11QkO0LA?EVlc5p7`#HR$ zcpt{QN5obMGruTkjz=e<6~V zygbfc?4!t?s8b24=3(SdqV`n1Yz$>>D1$;yWvUP9=j(N+k&hrhPPnO!;z{8-4hrpu z<7?uvPW-Z9)lSq-)=t&N)X4=42J7>+(}3f(iQ1XkS;UD3)Cm|T=hP_7jB`svtk|YF z%|?tmj`ZoeIrHD05i(!_M>q9{CXBWuIQ$HA_$X%5g;}1(tWx##Ii!!*)8~;6*3&QW z+955#`7e^^ymsC>pw5ez*nj?lliWI}4yYH%#AxBh&^Qdm%4;qBcxdUlt;6a()^J{( z-#Q{Lh>LK3APpIF)Z^6Ji;maoMYMWRz1VDpG+G@~{k50sEnY&4m()wm7D%JTC}>Dq z8pUiUA*U{=3n)ur#uvm(;^oK2Tc>Lm)y1vj;uVr_EOoaC_A>N>~g)g-R-WR z?C4Ireb)VZyVs2wb-`GePzjQP*oi-KPHbwm$x@V8a*|AJ;ObAXKUaBw617R^6mT`D zCTTZ8i&}B4_Og0;>v>4oS0L$6i>Xc2IOFOS)$L4d(wb9ZnhZ?T&Cr}Tu_Nl0turi- zyD|P*=PZ@ufrJy+SoR|?E|830(jXbGv$Q|gaC@>Ta${B=z~KWmmk%=vA;{;cz?7#c z-3%;Lhc3%Jg1RLq>1?dzlBKc(hbKb~KcFNCF$crMv)P);#$}1G=fTp;JaIX9be3iDyukHE@ zrIl5vuiQb!^-+6r`d_4JxBB?KvzVg+d?3?V?8OD>@a1#c^PQ6yJ!S1uHd{UFavt?= z7~?#!%MyV-H(V*N%53fyj51aPT;H`j_imTcU7n^{$&eS)UMVBIQ%QHt6mmsjw-l!Q zyGGl`i7!cig8cwTYTCD!S;UwK4kG4VAk92m%!qo}sj2zjqc!OAm+-1yxIs8d@^?uN zk0A0}agv1!?;mE8OV*9KOeM*fOqN!c?;~W17*svM+BxmH>(9&b8tL<1q@_~zaR^Kv z!dRZ4+gMnr#usjopLmrGDDr~CG_q(%U2L4lRom*;0PM#@RqN2jaJ1S6vLT;27lpkD zC#NG6XC@nha}%@WWHTn5tUJRw$X}0T7fMOzK`xmuyq+t0q3r6wmDNH~P%M^kn5_Hw zlVf9JN$fN!U2=-R41sBYYJ%#x@f|N9Kn5xZ>kVr0;{=EbBd>0lV=s7Gi1=!38R8>4urj3ND=19k* z>USw0&aHFU3V7oXl&o;8hq;XkM{Qs+6JOC?4!4~QN6y;p_?8N708hMK{{Q#<80bOk zAgtOUSja>A2*st)n?<)IOm6A5M~_q1>aTY1Hj2j@*!#eczW{X3BB%$95{cQi8@6qb{wiI- zawE14S5VQbk^iZ=39Z{9e;O45+KJ;x1T}-_L|y^%30U-I9|5L@G?#wc=!*Rh{w|A8 zT7D!O=vN_e`UyXQL3iumE5N-vcn07p5N)tMBJdT!M*#moyNE4zgbdO3liH##Iz;E2 z(X{!Kuo%;HMjsKJtpw81l}S*YKwfEz+Wfu zHwgSqfSGy5UzTKA<97fFLJ!Hy$_yYvlZK-&QB)gdc&!^spb5jAt{8Q~0JczKIk-4O ziXg7JbT0)2qDV_YP9X}o6eQzL1)zxcz~46L#UmcvR?&sbM)8_X{_ z3;Ho%Z{_v1wP5Bm_fWdl^1TTx#IXl$gBsx=0bIa|VYVTMIMUv7-!tawwg2q>U`5`Z zf+pJGfz4Clef4&}i-zNWh2^g)!>qdk<|QQgY0EvIu+ zgQMke0c@o7p93Q)12CrifbA*(7*|2S4iy6IRAIm_)drYQ5x{O01?*8Vz+Tl3*as(_ zm8V`1!&Ad*L=3!ZB8NGHwuauNnqd_O98n#B2UI6uQY8QnsxH7oVCv`I$83I<%ptmh zb!3)Rq-t{d?%msWlNT-|C%v0auk0~Ab{#?_Tjbv&0O6!PAp=7S0>scy{t2bnykWcM zUsCp85wHj}6$kmRDf?RlNKKIchQPli@b3uxHULz?9kbv(m3{*N(v=6}S^Rq{{T%}T zfxzz)paWn2M*=@jfTTt$;oZ0Rhg9&72>fFLWPh=8$XcDJ1)f8hYu>^2C0+=H667kC z`7<($%Y?20wTdPKMGHAIiz~Oj%IXR($C4mV9gx3CpEMBwSELQ~E1!=|}Y;^RPZ+ zUPoP+58)V;1v+R?n>e&S=EpF80Oz^)n8rszJP|x$Jh-3o{g%gPd@tw^KCSlcIb)jV zt^CiJ+y_4S2Q=ev68J9!Xgl0!eN~K``pS5D|BIUN8vhWfr?kc-5<%QNv!EQkonnGlKMft@Hb-q#h@^KQ5)TddiJcb*yc7W!uIJo9BLyu;|s*K!YG@8RWN#@Ar|Kc-HY z(Yp39+UI2@>Uduf!W7cV3R4Df6vGFJNEoEdTB#oSsM_1QN3E18ze6K%M|&vmyQsKA zU6Qm2YEVK!d!*d$`ukHO-1DKNmd~WR?|cTTGs*4GR2OBc>tuzF{#9bPLf7SXP2U5t z*;nu0ehrqbkKYy0FJo_jM%T+%sxj(A+?KeC4c*QwC0FuMFCo^2+EQ zr&orXT=g=#uIiQHu1viw$Yln~`N6r46g0KayH+jXn4oeCL7J}ps60=!olr5RZ=&un zsJlb7J?_T28Ns<3eL9q!;w!<)TlePXlUH&{!WFLf(>38F?(t(}i5(&>Rf;oqN z_wbGfc1>Gn0Cqe1-2_gD;@@VXWM4Fay+9 z6PLq0)&TJ@nMhj#mjjo4d<6_jYab#(q`X3{?Lqu03Kg{nWl_-v7D3umZIP!-iheK# zVv_qa282f-Z|q@?7&AvZGY4I$j?1T5ru-fO7X2SliuOSMH-K61>Ix~3y8^wWig@}8 zi8D`k`xvQIaEH((VwGMaH$+O6O~<-xRI0nEzxV$!^-HMEo+c{4O%wf}1omNqRbqle zsAMKE5>SvRvVWZY4sm@1N>P&;dRwXQ)A9l0u$qHmBtjH$B?3H?&i6O_pP`w+1>l_# z&`y5nDG|CmJ$HNd^8D@3rIye>7dno?;P>qPy;?~R>iK+l6^S3Nm3TUdYKPp8?bb+y z`8ZLxbgWs~6tO5iDR5Cm#^sE{Y3B^ zE#XFQDd<4JI_8oE-%$s&kc3J=bmDF+zv;h6C=Z_zzsdLG_i7+w16~=W*i+EGUk${o zAD#ocwM%}cY+VFHl%j7cg{DGY1wsF!mBp-rbUT&4sf*}(33QhTm*J+k6M~oyLrf>0 zslRyI3$$}=12|oEi5QD$visXTyWchBxTmHV@{A{zZQ#xO-xjiRr4Ja}6}bC3DmzDj zWpmpuLz?0+a>v5v@L+|?cKZ2a%Ch$Uqmt6(5!kD!{4SM|wgRIU|I$ZigtwqSp^AN2 z?d!xMC3Y|E0FUHy*uM5}9%`N#LIo zcmUvbAdpFS3-O;V_@7+9W$5?bvPkuQ1RZgjO3V%awn>20b1%AfmQs@hE)%#$;0pwP zioj13_(K9@p-57-_#BcXO?djd3?H_B$mDi|lb$d96|~YSL?1!?+uP0`Iq2cox5Ck2 gM|dC@3AY7@J5B9tv9R{FxD&q~?}{IYk0v_)54*ZQ8UO$Q literal 0 HcmV?d00001 diff --git a/simulation/__pycache__/DualPortElements.cpython-312.pyc b/simulation/__pycache__/DualPortElements.cpython-312.pyc new file mode 100644 index 0000000000000000000000000000000000000000..5fef09c4d2c4f28c107d40bacc5abd027c597403 GIT binary patch literal 30708 zcmeHw32+-%mRL7#5C91f0Ph352}-2IOOkcb6c33KMM)%eQ5<6+G)RGlNV)-Cn1)SH zyfxrPDUmxf5w!M{pvhDO#i}x0Gs$wYn;p*XB&(_2P0$7vqOhD4H*wX*n{0}L$fd4DF^yX!4efNXc6^A|r44h@89?BMS0Xjws1nHKHPK^~J0a4UI4)ZtX?ghz|VH zcc_tUMh0#UQ_9HSp-1$L0zw9+j8Q@;mr+5;$dof`2$`5H2$>lTqh)mCk}(OB4Y7F= z>LA7ByhAbi_a#DoNAgLW0pfC<`R_{;u>~a72(hNL*g_I(hS|U(gqwC<@q`X2bxD$zn zBqK0R9E`Nc=X{||%bM_dRbt%xFr@g1Z-(|k$WL+95M`AN_%x&amo7S4m&5KJwaVj? zDNkHBZTDV?D;XC%X}{=6I;w70QkdrSNFUCr&I{inN++PVghogQs!8ohXR+T|Ga zxGr9D+r6%-$l(;d^~%p7?H*{Ka7FZPWfrCE|Ex0WQ)Gko1?|NehGsUkrF2iO}!;y zB!!UGFZCM4FT6(4$4SL6#dkOLHRxKu>^c2m79phGoRpJcNg^ZiiUsJ&s~+CrLNFc*?`ZO!d66qrU}o~WL(?b)7yEnf7pgm zzS@Cfmfr5+Zj0MB`MPCn>e3`*x#V$8j$0V#nEjI5Ye8i8ShXyUB0*#pXAp~1hE;{?d|LD zv2`6A816YU%;NaOr7q^G2TQZy!y+M$n{7771*hZnBvyC#4FX1Em*v3E^N+!K4o1&( zuc~zcCStVm+19AKVV$DowPD9UzgtB?#wweCb9MgA( zf+bONZ8WFuGfE<_2}oD9**6Z(9=tIyI}qF(mh#10A8H$8#`2JqH&%b7tzMTyw#O<; zV~UveMzi)sWcxlhAY3f)3y&76`qP}e10B*o=#Y1*l#X<=CWAyWv3>-en-oXSC^+g} zKZzIXeagMS{$qs$P?hu*p!VMC%d`UQnci{J9G4x>q=U#i>b708pF)dFFZ((e%IWYOZ_xg(wQ%HV0FbONjrZF1`yV_ z)Q~1;KhT@x7ZZpQn;wQvOGkK{p7#@B5i$HS7v-08GEUA&_ClJP;-o;UJ_K5&;1m-$ zi#I(f3d7RZ&FS=&QzW4(Q`ux;2n4lZi zho|REOkm}}7rL6WH5Y4`#+pMlM{9=SGOr658wt!SBso5fjj$oYYC>PJ8sbuztZXf$ z#Z@q!$4GXydkUB&;?07&&U(oNk1K&Sv3tENi;a!T$6&qi0>z*7#$|3;C*qRH>9}%| zv9qlGT3q3s!blcLJd5KPmob+vPRC^qr`v5+#}z;vTn;B&iRl#tAi{>l2FtZMF0nvc zFs9q(fCA75{291ikb>juLDUC2B0R!U1JtB{YJ9xY?e0Wa&}u}`M$B!3xeBX{HfK)p z_m9lk%;NYshD7f-kWP^%4rG7%YZkLDifuR2(rk^4X1HM^xQ{Ym{Zd zv8vMEP|vF8>VvJzs&ZTvEwCzPn{J$&J@OGkJb@v-D@p@*aiKkxsOb zl;Q8uH^I3{0V$FJ5`us8!^s#L`x;Ua(KbDhN{APvrbAEUf>bH_mIj}U_$1(y6CW^F z3?Va|4A^+(K%BlBr?16ndt83S#Q^uNxZre+U+~5itetUPf^~=84{+HIbY2Aq7A@(t zYm)86hy&2m=^%mp`3&Yzs= z3pm~$ct)VG46=;#vb{&VuHH|Bpz?yOWB zo{~fJ+#Fg?ST=Q=N~BrqQ($I47iLaonAc$GwCZ8QKk2+OEo{v9_;PR!b6kK`$$o9> zl6Q|Ui*dT0UZ>sd-s6)w+>SlI`pyGG$NEnW_t^UT4%lFqV%cW_ZkVaBbGy#la4*&5 zz3Q#A?qN${aD1i=+3S3gIsi913T`Xex&tTx^F-Dymb!~qZ3vgJ({6W@+cO3q>BpRx z?e2QsKheVtDVb@0j9OOtziMctjlowmWA;laM1Q#~89 zS!Hn<0wc>ci);y74h}|-FJWj02XPZ8y0EO8Ql;B)!)faq z=o@BFLY~*)=lLsefKivJzSZ}3U*O`htRyD*BiZH4GD{*5a4*Y>MgQ2c%)F+`3FH9Z zjjL58cmG3G2eS63-y2%as)%WgE7}rXTM`^z)>a~;-#6PA$!iRc^LgzHPCjqXL+xH5 zqB_Hd1Q(yzv8;U+669*EOeDK!S=O7Zjb&L@0bIYBRvFsMXK!0j@Yy?-Wjp^|rGdQ$ zSvqT%RfWllZClXtxqBB+@wr_ORo&}KO05;Bb~@$3Eia)<$VqM77%(*QnjqXD;s|lB z167s|5T5K4go6S#$JE7XxAQ7^X|Et~N-#CvOAv4KPT3MU*b(^R8}Rd>Q~|fnxpU2P zC+6FNhG6q8Q&dwHk(H4;4#y?-=D6g1^M*wo87ER%6)RF z%wcH<5 z)h$B}{mSR{{{{#d`X6>vFH?t{C;Qy5ekln03=pzj3W6a61dSve6NyR!Iyvp}BJW3@ zBpxe8AM_$_oX9hg+y#plwJxF7asbJ>}8)S33sOo7+}*yk3e_ReX&S&}M2^m%bzA`k476J(Vs@@b%*#p&#V zn3geru@CxqTo+AqkTZyV z_?&#XVm`#A=j1h}=N0{a6Q^SeIsF(}z$dmmSr?{+Q*k*Qp(BJ6l=^@nV9MT=G3DP& z-$bbs%3lmnQk14p{DLWdy&HJr11Y)zqg}`p zVr|6v3S7%J=!cTDT4bUXd0a_iwg}8#f|mJ9x#EfT47Dy5>j1i+W*IeE`#kt!{uhCm z-k!8}7mF~`+I9rn#ualVWF$+UGEPNY>BdoE>eK4G<5}&)G4+>m1sTV`h$~~PT*~YBzcU*pPig9}4@?qPg6F%vQ(}&^`U-Krp>L<8w zwn{!hR{Im&gX5-u%r)+D`JCerglDZ|yp=rB4yz=txB}KFpl4NFAa9>)OY`>a_!Hy? z!17bfeHEB4DE%$gkcUI1AWf#v+uKkf8Mo0`XbQp@B;b`r>(i! zDj{DJVBnv|EC|zAWVbb2yi=C*_`=3o&JAB+r!I|O_!@yd^G;JI)Lji`2Tme`{{(k8 zLNC~$-g3A7ywlC1{uX2pEy*HVl0{;z#UjxUVu}{AI#>^1X`%{tTZ>qQK_^H-9B72_ zHStqYWz^64QiC=zJsy{WT`fPS);u5I5@ysV(eZ!=*OV6IREiSdzqkHLeufo&_1f4#Wn1#};4A2E)+J5v57#K1EQ~` zIj$UZj=60u%^utZq6qVuk04>o9fI!Z29Uk3@k^}pICL!w3U=yE_<0Tj` zlUG%(R@L0yeP?%=S=e=Na`80Z)U$MfZ#ot^KFBu>MXQFzyluS7k^yW!K#U9y^UWut zwWq|q1a4*Z-G)02;o*h)d#^2a@r|!88TrPZrSp7aU$pYDn6;Hxm9M^F){j)>>p9rk z^;}Bcga$WaRJUWM_hT^F$ePWXYYWQecLyEw9idj!IXFX#!k zK9p_yINPw2O>B~~x6B-j$u!>@cza+)R>sT9g1%*00|=sxrtfsUC%LVZ?4_bk5g;l7BWbLOyMvQ+eKTQsNU5v7(}1JaM7gHFw! zn(GNx1?}IbgLU%*E2akC)NsEw{Koq`epIz!|6y&^^h#8>8|wFK%lxrE_oi}Q2`y3I zQim$v)!mnTpuVSGsC-`sL-l)-pQwMVUaI_E2=qkty#e`}Dl5<&7z(t{s#jFSys9{O zG2FYbd(rV>2fwvzshi&l$R62p^ubW1(#X)|{t5HiQAT1PF8iNViJ-UiHPiRp1`8TuYvxAo(>T5rh zLrySCD=HCkzFU@-zO1kPIIlR^`JN+G_pUpdx8*ZRrtYHyy)ko1@YK8wJldQahh`7W zoeK7bPXrDvYq!Pn3*W1}U4N@Sv?E&7z~?sv4y7_33ZGoqu}~k`vS-O2IP_3^XkCo} z*Rv_L?tf#~FRR-ApX_!m>GoNI@wXadxkbUY`KC}lpKC?*$_lj3Dusn9QuNBQYB$a! zh`OPj)y|!P1?sg>H*dCvt9f%%_{x$Ypj}pVrLf;<{?7TE7v?VnuSD};p48z2^=#Dt zbCrG-vY-v!q2LwXR2wejO)U#MdDHI2{k*9!GI%=jwXZ)^*)bO=_O5FLW!&Mn5C2BL zpf`g>HTj5Q=SPDOb(3-uTl^6TLo*Ua%76{XxP(zK%5j=eF={4joMJTJrx`7TrD!wn zqM2;SnFFC5NYj6x9y5Y4ESE97BOj4olUj{vGD|e4$BzvT4*|pEWnI(k1e7-L7eXz5 z+3^>@eG8wz>G%Ys6x3f>M{oSWU)TRh2jPgnFb%mZ^+PVV%K_G}mPX5QyW90T#vW_2 z)E{))J)UIjJ(<8z6SfqK_yi>xgzXkMb|(aTAr~^@0!w}seV4&;q?@dvP=hGwf?)oP zmZThIgDBy7w)*X*%?tBDD2AI6@_-_ZzAq646#o1)<3UMEU6IKYTl|KVa%k<^k}|sH zU?h{N#c3wk5-uM~lchFuoMb}82eCr%256cz^Xa5L#W!g*$eTb0bJO#BJ#=awM>7ht zeN5{+ur%0|7lVEpj4HvlmQjMGZb8bVFjG6B#z_?g82%MP4R8ap=@B=c{vxhuB54|! z$|R=a1|CF<$Dt0zY2~7>e$#^n<&2t$MkN!LO^FcXiix-&gHZB*me|9YY6|sxDl~D! z4}!QsFkuH=6rPD{(U*EN8ZM8~-q$^=tW>K}NSbb}p9VXFA^`BD07vqzLMHouj#vVX z)@UW5L@r_UDA*%DL(-R<^cm5|nK2&`v6x6W8NxUigDNrxRtS)XD|~vtfI?b*-z?Tc zgv%8p-0L+mipv3bdVp67N`X>@pUgu-OdiRZi-LPG#b4&#B9sid%b0vp(*p8snIMX0 zNC=f&;<)(BpP3WD0mf*c9u}rZEQc{6Ob|V7+fW`6b~6Nr0A>*>orE@&R{?2H&zUN( zaWK5tm((^Jx zKg6=3TS+$D^lX0AghfH_q#ruFD2m~SPeWI>#jz-maQY#}KcND&0xB!xp&1rdmo zuNz4oEyOjU50oV}@T-MiUE18-F4QF5<)rtwiTDV84@=)DEPDQWzt!L1-@;ijM3mE8 z{f*?^K?LW7=1B8NX^w<+#oQ zCEO$SKuTSm7ME0by}lFHq4c$Dy9fu?t{UiTBN?e?()YE=a0?07CBv;GtV)L4Fbv}l zt0$~k?SeVFzty`hS*uoNhq&Ux8rA0Q5K|^z6@95USHiW4xt}yn_Z|>)PIQaD)O+G! z$}jr3Rz}JU26%4`6WpCyr zpht5-4X);lGa4^iD3V8{`ufcvH2r%zT8Oq)0-_7lw_?ieO^-k7*4{+;Y?nK7!tpP-UKdl z^K<_|O4*tQw=46=iTDQ8bPNEc0E=V`lo7Z> zlQS$VjxwA;z@;g5vc*@@$&QnQ!WQqfX{Tj!$_sYIu%8?UJp&V$*$#FO`D9JwjOQL5 zmyrW=u#N(pxa2A-nuubDoL_LE4vBqTNCY*@HPAH3yw0m$a$ zB->iyyucv{T4ia<0VMpXQv@sn!FU#wbHO&n!okA0QP2Q*>}d6Bo1TJmM4q^@*X49G zwhfV1MH*BrX43%$)`g&DXqXIVLhNGz8%`m@;dB@3NgQrFlQ`d$)Waa|!pZJGgqLr3 zANGZXd}kp?N#HoG$7^wU$T&ewoR)-Y&F-;~9adzqo z>RE(SDITAkoV8++hzn-plF?D0ezXtGt}hXF78;db71;()Mg(&D|~w2ZqhJ12$n z!M>scmt1ZJ&wOGZS`sGnmZqjA3ykOqtT~Ej*kN>VmGPM{jbu(7d`ER5AkcR($9K^| zY-72Acm$AJV`z?%>z3#l`VR5wu)g?|H zyucuygBm;%v=KXNb`C%Im*}jZL-cfLi%AZ#u}D(Vj#H~3h`PIJH!5mX;!wnuFr|sU zEvuItqa}V74hQkcuaLr#{bNlCLh8AOGyT#<`?QPYF-AVk!i*qB_T=zC`)$mOS}ABh zQBWf!>t5&qr+ZB3Vwvm6(|%=LH1u`9-YwO1dhfvIFvRa#6~fL zwhN}sQE?I}5h+perX{LE8xR~+H!4wWc3|piTti9{&R@lK8zwIt5>K1aaW$4i&@HY@ zqN#B708S%dVxroFB(cb>6FLf~4GB@p>gwt5hvNrB!-Gi0VDQ;L!1|J&M}2Nw3ZpDw z9}=Lr#BqiYA@H0PI#!dY^%1~^xR!Lpz@>}*c!1GBLZ}EFhPh0l)Ht87qL=*>geNCq z;dKqvl@3QHX&#gYX7>(IeG`a%lc0qDkJy$xbUs1nFVXo=;Kb#UFQ1brC^z9Q3m2Txve|s*bWv!N7%c-mxd)FO%~0$0RAEIxc3lVhFLU9| zqw_~Y+kVvcJNrM}A8qZ5nz{q(7^oKw6(R4!nPtQNfZ}6q-rVK)&fgxtHNLE^jhTyQ zdsnl~b7R3%%h}bd+LBQ3LhVxhL+uerx_NH?94aGBALkZFN_Izc_e8XN*0iReVx^#- zFR1@WTfbIP9lFeyGzPk_AC1A;e*e5baA=LFa^@O>-4T67K#uyb)3c|8^~<`dbvVk= z5%R@MMK=%6AO3cKpeLp^2dz=zH$l(IRt4=H1v^5c;Q zGJfCDm3>3}zM*CP@QVHnuRjwx``XIc^ZeQK5yv=x_QFHGYh8*UJ{Il`0JT%!!r-!Y z*NS#8uid+>-T%1;s_+2f?=~5A?$&WuTGD2i0+`xDz@lcT(;CN$2QnxGx(NjlKITDxuG@b3(KIIz^lw;g)m zi`d=EwHE`2=1%k4(jfDZwhF2XIoI`+sT@Cgl@4^pw(VYQ<+pW+*&MvKGPLU>?N-eD zDh<_CH^=k^!F(S3cU7mKGtYVFi-SWvjC4=9`p%IBDPP^j8``3}_JAaoZJg_#eLYyk zXP1S{A@7~y@DN|o%x5<*H1oQ)fMiXp2i-uj?!2~X)nJ-C33`ijAtP_D3CqLGJ$ne!oLcs47R^q9;-o`0Wb4QSxt zQ3P%i$glj2lIoi0q%l)|aOds4xAwlf|09!it$Iti>rOvttA~DDy{?2DpyXB;#muN& z90@h>=B*<1%`t0cij_^eTjlnR=7RoR-VX{{K|3lPcD5U+l_h*}N2 zs%9;Phjp#u=S94#8hWjyM*OeFZeCRztHUlK zsk@-%vHEr~h?&>3D77JGEDW9r9u7GpMO%5})_`hNXTVd{QC&G4wN_V$T33a0X>%8X zo@H$n9JN+gLzkHHZtk4ld9wp-6?R4OJJ3Y07^`?=Rj50h7d1k`u>wGb>pCa~%4nZ& z3Mf}|O*iZ2>uzqF-xjKi<~9VBvD`uu0+JZbg&eW`GWbTjW-bn7t(x<1o}NE_(>8Al zjYQ4c1AxxuE9N@1K9Pk_Ma^vxQ&bK4bJkS*V|hh4r{pCcnU!}t*fqvdNGJ7PrHRQePzvB;0Eu4-N>;pQNp1|GVLX!&zBLzSJ_lXHx zLd@N%J5!OG9Sh@;g8j?d4l$)MEc-xrPZxSiO!yoW%0PGYl)5giHOKU&5q-^SVN`lngm@nr!=?|) z@0BlWUje*_gM?tFa=ri0_4%v%%5e7w$L<|l>i(mnKRp^bbN->;5v#1dyY0@lg>Jrb z*J9gG_WgJtUwP=ke!lYTob0vGCau-;1P@;<*Ng9YJJT$sN2+>03TDdNf(pSNI} zKPBD#g$J48Kg>Q*GEyUrSn@|&q|3RT`$jsX@f!KaHf6j;KC)LC?~tEUEAN2?8|(C5 zVkZaS;3ll)$OhSL@hqDSn6gWFM^#*Fvw^jkTS#H?q61>7NVKhl2Vj4QJOw&ZTE7~5j?1&du#>%ULW3l4$wX(`kbMQ>8qB>@& zSSziJRn)|aE%CC}^$LenK0>coQu)P?tB7AH5(`QSAiz)Q(Shd?GT-0Ir010|az+8V zUL~Uf-7aW)84bwwwV>%@H`T#XX!_8>A5ov(BN{`@;0ET!@f~^_s8#rh? zaB}Dq1(Y(t{s>&6=w<&Y`u-U@zk|*qI=I1Re~ivQN9T9Z`4`~yI*=b1pM3DZj+=S} zT2kt}^oE0-lGG(Wna)9q!V;noHBDBhY}P>`N+(lwpg_b;#HPn@0;Q)$j7(hu^v$#z zu{cvE1tlY#Q4+gc1Oe}-1IM$5-d z4PwcdF5z875ve{2N==^Hul~H0 zIy;m48rqrV1DpVbDtV`mgN6d7JYpqzCW=H~>J5q=DJuQMp7i3DgFF)@qA&F(XY4e@ zqB~9T4MOj-V7T%NH2?KtI4rnFT9lnmG?w#F4YfB87sJb{GQ~1GEb8#9x|rkKBuM6(%r@N;^Tj zfvZq*CQf3LiAd?K7Jaj29Q6;5!rB6+QykvmrUCmMr`GG??HdcqAH>Dsa&X5!L!EJ%fYC1}(dHTlNZ9d+pz(=tvWk zU55IwKSf6Z&L<=2yNV97WuGAD23i%A-x1jz_(>vBK>R)?no{yk_A|`(2p#16h@OE6 zIoZb;{A+Zg==@i7K1b)5===>jU!e2f(E0D^{10^gCpvKW5KAdV|0+7aKR)?_?dK%l1x-|<9f?NZo z3^^Eb?I@h?kX<_p%|oRR2sW?9)o`gN8b8`>o^ptLwj^>!*oENE#8i6>S->^5T+qXEu>kgr)$}oF$O{zMU&DZQ++_6-#G#ITpG^a{_kRey^;?AY!rSs96 zJ~3ZZkK(+W{qy}H6PV1dnCf{`{T!Xx#RpG?PTV>hHiSFxnfTJisJSU(YJyadBlOMp z1??fltqW089T?op@5y-+T=vGBnjyWoEciyK;BLvCk}$(pw60XN^YFjCJzBhDPKDw_ z|3{{+u+zV-y`>EegT!T}U@Kp+HSFgL_Re*#3LEW^B23>=-yd2qynkxJe(wlx*%>X^ zB_yMze?_SGPJOgs%l!i%^xW&gqZw6T?+@wLEg$6F%Uc*;w(efB?&q!hm*}O62XxeW zc5d0Ak>f)@w=NH# zS{Xje51;+{wrJM~zjq`udOmXQ{IbKj;&AZ}S7c%;GVS3fJj>q8E8c57{GYfMb@=#m zzICA8^Z^Obp4u*nIM#qfnQ`J-rSAC(N?pszzroFN!EQdMVkM`F(Wr#>%mC0IoK zAA_{a1X9-GsIC;HxVdw~AjNH-wFM-b#rUP^V*FSUNc65Bf#FeS&niK#7^+{`v9vRC z+8(*!j=bTCRW~hOjMZ#k+!Yyo4el&llr0V~d?PY^Dl+;SKVXY^u0cP&D(RH~iJPPk zN@8_wi`OEjnaD+NymrsJhEf;)HHb_03$g<2j>J3}qc< zZSlgISb1fvq(YEvted-N`TlhiRk!7H;{BAu2%&li-2s4;Pj!X zYwK7xh&6JESR;FbSVJGrP7unjV<(UI_Y9`w83QjV&@dz`OHoOI#xO(@fkrZM)&Q&l zFr|bw0Bd8sR4JW7fW@B%c??<`5>|oJLyQF20ckora8EWM8n~y{Y_i<{K9qD%^`eXK zJqcTb!-*qr&E2#-el-iAfIql#dt&uUU$F8wjze<2Y6jNDa}2)5^~xkpj$lD(M|kJL z;U(G9`3DUku?%Y#&Ma~dPJiUM3iIFq%ma80NltE_2XL<#%mZoqJYctB4>W_5j^(Ww ziOx2xKo?fv9crz4*V2gx!*CsCuxjpbxOs7Su{U!31R!{%<75KCMbX`Qjt=l!4^q3& z4dDBf0!P4Zgt`H{QSXhcqa@h^|$>HF^_ zlJ9~ZLL>*HFo{1GByFTI8=Xwol_sui%YEVelDz=*CEybo9rxn!6PKdCaDB)gJwa zCy5ylF(dAHaBl;IEDH#kCXJA3t-3EtEZyXeV$0E!L&KH>P7D4i2e?}w?=7~#^@;Wt zvZVpv`4+N6!MHYJCxvlsg0a*NoUcgy*qz`TfE)D54x$@k;*!DUG%{!Y)yW*$)nsD! z6>wA^S|jXu$ij(Rm?sDyM5NQck_>Y+NdP;M{@@M@2;h1Vr5>c$4Eb;t(o%a@b4P>9 zyahRL*%46SX|WqQvpGQ-9?Sf7Ca;xL-tN2A7e4jD*?VUf$A04av1@VU!R3E*?dR7X zOz?SU1DcOjI>N+kq|=d8j)-$Q!U`;1ACl-Jl0iwVwjLOoh1Vk|UW-+2UECKrapt43 zZ{QV>gjDkv_}W5f`9uaN03(BDiXk^PPGpH?rQVEie;1Kyp$rTTtR?+^ zGaT;GfmF&&Qgn=ltOyte075>vi$=KBQY*ro%>yb9GMJ}Pap0kqAg$U+$u~==6jnXE z7w6$VbO=-0fj(qPMadHS2LuCGKtkbG6+uOab5@XwoI_udj$2hJGR}4)JR;svqtwQR zJCwM@<4dFM`M*Be&K|(xi70Ce`taftfwx6QC=-d_NaXW3&gvw6e+_9?8vxr#-`Ac^ zz@L{~T6Oy@NIS+paNTo-NBF#50nM7$7&r#j7O27}W)p3^v3|wafCu-xGslNdtu*fA z8}}{R;mV4YM!2}5YpIiO>|J_;Z#={s4?Wl&Id%HydwJt&un5R22po)+SVC>Ljsy;Z z?b1;Um~QpKZT$NnqHDJQ%ZE{(2awH&ewjWm7y}VGNAWU3SD7kK;Roz z20dT{K97!r^AgNFp0>u#YBv6!p@&xw_6!{x=o~)wrS_^gy<7%38Pn@9Vl5L-l9$i0 zVKu(UhLwYZ;lSq+33jFpE8C7Eu?rlK_imKJeu-@+5iJw`i*RC_q{@WvGC|8SOv-e& zD>vDb$;tQ*qR9|^3X}11LFV1&W-bP#9v^M%elQ+!IDw(LOm|7rxfnU&VkGJ}pUcJc zVGp76_plOV%L_Zs4I-5dTufwuAzMzgA1SYn0XLImW%9G-Z_;p)2e2~4`zeJH0>WYm znGtpn+Wis;o4Oi+Tns=^)<9$|`VBIc+;|qT899EabExObFf(LHfIbmTM=~OA5Hf5w z9ob|6f*>|Me&b|`7o;)5!%0$Gn@o#x(wsyB8(rgD#*E7|xssC0oF@DwiOgw)&3MXp zg1JTVWcX8Z2_`9`s7!1Wyq}U^P_bpmPZZdICino8lexZMU1?AX+#!G}Iy10lT?rNj zm=40cKpaChw%$aIpgj^sog9lTuvrypE~&BLFyx|10WLrg7@=q6gw3rM>i~JlHQnB% zsDrRaR>B?`pSsCSY;n0Lk5W@a8i=JuZk^tnj%m58)V3j32kY_W)L#` zV!V+tMIe%-_`662hXD|hVhg1+&GiO%;9n^awNfA|&nXP%{XlY8c}KaNT@MLh#$H-; zd-T>Q+{97}cHZQEmSHsg%^VFNURas6AM{w+c{kj%?qD0AZCS~#5z5P?OozzFH{gNmScSr5j3wRhEb)ZvOdTDY<* zdrk?e16^8V&>uc=>n(_u&THXrqJ(Y$w$;1Brxzxcsvk%m9FClE#tN%L7Z$XO{3y_2sqxVw^U12rCGORTGrBj0SR%!3RfVB*d%K032q7nu_{5`8{g0dGW zwQxYF7JYw(4hK%$04LG#sy_HjAMp1*$)$Y?;myuqx&ls(!IpuQ;ukQ>LE(W;6*@>| z*j{wdZc;c3G=jbl(fJ`dMCCyw!bH?SxJt5I3;l#ya5qKtBVVZY3MQlb*;CLs_)mFo zKLf`AX!_5nT=HN2XOxEg&jR-sRQWHcqR(U!TK}0Ui{AgJD2LAdtg?VEeAJdjn?Ea1 z)80LH2Q%819gD+D)r;_FN4q4n;q$}rm!0Sl=vVT3+^_7W>HTx=#}vLF z_e*jhH2Ro=*W=1!+D`|!KBn;fcy|@Jn1t`g{d5lfCOxOcoNv;PD;sGqJy-JA6uwFO F{uemCIMo0E literal 0 HcmV?d00001 diff --git a/simulation/__pycache__/DualPortElements.cpython-38.pyc b/simulation/__pycache__/DualPortElements.cpython-38.pyc new file mode 100644 index 0000000000000000000000000000000000000000..5af427ef49cfe99b3871709568c674b46bb2fc05 GIT binary patch literal 15929 zcmc&*du&|US--D&dS*Pf$8X1(yghk%V<(%U7>}AwiE;vNR=o>s`ihzQmK&8l_*s!wFGEEEn6TI#P9p= zo$({tZFX5;tb5+~+;hJ3y}#od(P&t~Z|JYP@}FB!l6q{dKk`O# zKmz(OcyfLQ2PORiVx2PJmnrFmNKud zX`Jeov97MGLTy*7imhNWZM@Pjb7g5EEAyFjVJ78sv{KnI7Sq<86A*b>OfO`eK>B8x zTY{!Jn=fQdGo?y8uGbYAME2l0ae8KAu4I*p`OM5zdNEsCvStR44xTu2^w_b1Bi9zK z{H6CFxt=dzxiUR-apd|;X>snB$i6UB&Mzz#(pJ7yoEce47bZ&58ZBfOvPG-hzqsPG zn`SO+)fadre{*7Z>XKY0mv7Z3$0>H+2rDX%awD|t{a;fM(kz6Y0i@Bu8&Zd_NR+8XfPDP}S%`ogG?NRzE2_SKd(qgvg7(BUAD5Zttlj&kE z>*xiv`?JgWvQ>6Mqsy5rj-yme`FI*fAGykf@&OpF-^9iEu&G7_oLO*lKEWbP_TiF zw}LWZE3bfBU)8Ew@wkOnOI_9EB^!t%nIVnoGRsLJz7QKM6&bEq^YGcpTI15js61(B8INeDReKqb@<=K=Dvl`dug zwX!9hV74gAckyz0_Xr^fs+Iqz-diiSJt(vkyold z+pv8?8&Yf?c>2+*-}bK=1ix}0;kxZ#)v5u@ZwIWvs#c!jQV?+dupO{{wtq;u-N$Vq zw0(w30{G2iU#!`E*)qU3;!^(pPX{ITOwN)C-h2_P3#-$?l zPO8b#aU^mV`kWvRZI*|J3MKF}E{1WU(&A%u0$?U-%aXE_hWlngCl)ZhXgNj!^xh3($8VLWQARLUM(G7ZEG;ZLMkZS*q(Y7#D3H%&WiR#b;0bAN5E*^f%q&S@ z5VeftGxz|Nz)F!`Ku-(3M9gy&c>!CllybS@LSdNZG0h|BFYiX8s6hQ2emRUU z$@7#DgyjelC$vy1K6+=7s}PP>N%H2dV{CnMye)TgnjAyV*p_3YBc+Fr@hd3U${NHL zF*+_W%i|IV-D`uCM+r(3SnCk?`KYCM+EQ)Z(m<>)*!r3-{9vtu`vbNPA`P}foMBw^ zrnZpV#w}rB;pTBCrQy<1; zu52WH3NO}3dUrJuhxE4uA^HsDkN`&+DR#Np;eNl`Z%-V`!tr+ z&`fDzv5;LxQGJM_H@Q-j%uZuwve_yKl^PVtb)d2L~sgOVuK0GL|t3>-K2Zoew zBZAvTgtm+bqqkw(i13yX5%kit`74nvBckY~XY+_C1n^OJB_ig{P@e^Yn-mWl4fF7R(&R_?iWPjpMy0C>U`VR!|8(sm!vT1Af?v>R=}6tLTofRXo#A$Gp%a^a0}ZJ43m zPLSQDDleXU z??v8Ww_`Q?nB%t3cVOf$ycNQ;KcU#IIJ*OPb)1vF)@iq47oA&2RJ(FYwHuJ_q|>=& zysXRDEI{Ax6rDuypr!rQ9=zRTx7j@&3;_Ymw8QRR>#}$7&h*=MG2xEd+Cy_;_HLXo z%{5Z(vAgXatb2z|tHt}Z6^cH*-%ax+lxnZt!M(jAW%s?j&^;^FJQ^4Eck10YBsRm;aC8 z=er#gPk8YAW*vUt1azM8Ae97W9c3IJV65WWLztDUof}3E5)8fbOKv*9t+?1*e?xoO zFO!uvt7mP8$H9-e8pL+v4ZNmzV+8a1Wsa zl-V+@=W_o7^n6u1zQxj=taQ{V7Nd^Of{9#q{Dt(wO_6rAI|awLP!id)UM^Xq_H1B0%Lx0 zE?v%=0|P0I-{#3s&e9N)&XumTIgqqU$(xkXW|O}(H783;xw#MV+#i3je;8{I!S9e| zo&G-gQd}sQShuf`zL_mZ8VeacSpVc;{gbrn+9#@RT_LAl)3CV2)9~EnysC zL1~Hglrg`Rw+g_5#Z0MCl8yl#wg7J%7UcNccl_x^_;KL8DU@;pM_qRhISsHcAO@D@ z_%el3Ih)#ry})%r?ge=2vM-F2r!Gi?#tLbSA{z^yJdw82a8@8?D>O*WJk&R8F3KD> zSD<_wemgpSNlPy~ff^LT7zme;2`I`|T97sh<7gDd(P+iNWDEFcaEee2_~`hjTsqRz z%L9PQ@*u2btl|1J_Fk||5L$rh6F3rTxC}$fZ&Hg{n?#(il?Df#z+`r|U=9v!*M6I? z0|X;@Z7iQ#lGzCW@}pQ=`45qh?yG4*J)yaF^KRAs+pBk}T^b&wF&!Q&(v7vl@Wb?} z_|x_x??HPT{4){s1phiJxfgx#vg^?9JIFW5Ce6gP$*Jo=OiSh$ zrHQYP{}2g&|90+&zxf;V`}cGAeW?F1HS^L}-|YMPIp#bM64&#|K6rcc8MqmfN0Jli zLjD%DUK>pIT`Z)_<$5c#t@~tfFiUf}A@PE-DVNbLP}eOaud8q~WUVDx%n(FszjG+m z{bOBZD5$n#vT3l%JJw;rGX`Iyc2}!*!WX7l{(1PozwDx!(Hp$0zVy**x2>&3 zp&a+vJ!>%*Q!%JO<4}FH`9^z3C;VP)q<`7gtGz{-(zVuVFXm~od!ZDDpv2YMn??Ap zz81B6+;JqFZ+{jeTSSB1AtEA5z4V<8-x_)o*ijc&+bv?Ek@f`#$xawYJ4Jj5W#mm> z-t6TqP@LMS9|}PNWv|q!VD-DJJFHf_+g$h1`8#}KKwAF<;Y*Hum0aT^M( z*GH;s8~1|E+U#AT9qk<$wTskEDs_;5bN8cqAoxSwN_CIu+E@`(X-yvGndtp) zo~g?%@l1QYnYZ+KGwmxLwHt8)A5~-z|##i?h<6_6w7Qk-5y~o}!Islh_$h(p6 zM&65jFY=xtm;=anBJV&vs8{sljOu~zVf>7pr1g&b{kWvitx=^y4 z(6S$SGx7s4OgoU%8DKYn9Ic+d&0jNHsjI} zv^Ag&>NU059VqvEeNUhsLH+2E0;8i3PYTaL;Akh@RRecbqK+A(dZ>E1dZgNKAD%H_ zaz0T#igchlSUpxfjyO%1eF!GTlXf4Bg(v5GSVc{7nN1PX1*<%URi(W0ag+~wLf;;v`^MYp^Q-n z?9S>bZ^S8#IAx!zk3bnC`oMqsu&OiQJ*Vx{Xh~sZr^QJzbXR-n+3Fel%-TWmG;t}6 zI_sII)OHSS2e!36e?6 zSN?TC`#eKCQ5~|503zM?5S$}cjPev-|COXD&LaRHT+ z1rE9FzyY{L4L&TgfhJXebTV9%ND!Crzns_5mlPP30kE^3N1iun?*=e z{uE|->MCVLg85DRF>cX8szwW8DF>UZTd`TNYp5&|b>r)4akEeLZ< zW6Hlxq1@&um*k;@)((dyVhQP4d=tSnM2GXF4`&K#;RST+`t~$lK-Bq5$SS8V6Eu?i zUXnu)h?|zJWWLNNg_UHJo(YvMCmCkR;?lxRg!~ZQ@d~UH)9E_?oGjfTE!vG{lq)_C zF3COkmg`T}W@ak!nakvlU1GC^yv;C33_3eHB9W~$cse+&w%rwD_nC0C(f|Y@-!wOc z-3TF z?x%!elWyQlgMngP0IC9J=`w}1IQ+!KB*|&)gb0y2J_@@Fm>ZVqJSz1%_dcNjAO)H4$_%$$#;i8P61rmL{sSH z6(abhh4f-xevMjuixOnPj&;75AU{omU!$(D>y4n9Qs=PwY+=>~G;tpoc(7$m@Q*OV zfUR+tQ4xL0Kc_KKcf%f{It_H5rdyM5(16GV(9YEwB?SmGVOzTp$QbtmHqo{dqR%nnIng@8y3q*)ET~id1%1ff0P`q$3=Pa!VJtJ3 z3tt$$as^R}>r<1&)_?)>kLjfbY~5QMY(P2?vgsWihwozMc^ai`2x;)Dzq2bXE`2nm{I z5Z;Y4baq8)^XqXff|#x!Q5(7u0Q*&=NUvWH!)kORg>l>yh-e4ZUM-;;WS`O>TH
    ZI zv!Qq3ZU;r25p&jI{c;I4e+(^cT=V46q{PF-<3Inn3C$x`gCJOcSn^g4-j*fZ!(oKpoKXw*lDz z!I)%yP}NjT(G-RamO4*C;$~RKty3k`G7!xW51z=tY7zR^RL7ti9&PVXT=;(dX?`{?;M@2nl)>U@ixwKu>vg5BPU6}x-K?uuu( zK<1sl|CufC>wrQWr_j)?5DXH41-uD*18NB8Y;X6y#Y)}r?|dHY$30VEMI(G6n8yxu zz#9TsA_T;+M+nL8fu0_QUEgcL(}DE^_h65|Y`4Pqg7BS-eQ5ct-D+#5VKuKR^DX3y zf$bZBg%*Mh+5o#UYR6z@!z-545USc?_ln3p9Tz^}1G}fe7tsWJH(|HX2vDPsXz^|< z-4g}9HIB3PQb~?ZPF|Z#o<5x%aLOp>r?$VR59k#(u2`?}OA{JV&WW;~As~ABf}L-e z>hEIT+cwI3q}@6});V!%F3i~kqJ{ERB>2ErV6m`IeB~+}+s4XO?C}D88cSG)V;4i`TE|hxIy*fQ(>Lv8CyoAK_ zvovfb*~9XOQoj-=c?4<1t!nr|rDGe@uM5G~eG*ri6^WaT&}b>4r8!Y=(Y$qU;k^1? zG=%gF8AIs-p7MUwf+Xm834VMOW&R@p=>5H^HE3`0I~m&X=3m^alu0>yNa($_BRUFK zbPPlX!8gk;XiwxdiRLXLuf>UZyT5Sl#>AD;$qiXI{_cexsc5wyl{@k1c*<9hSPEDP zyLIWR8{G%l%us4hb}$NZZAB^CPfP^uIz+dTGi+QazuV@q5dXN#8<+*K0?AY!en*tb zsyett2b~`iRLeGkdbcF=RqXeX41!T*E6R^R5pXXaMiMmy1-^t`M5UBVsQ+iBee0ps zF;?!OE%3bCFmfL{&ey&`svqN`OBk$-)v@tkK~+DKH2yuqMh=MJ{2E2<7hq${=nWvD zd13O}2Vt*FQE}ouM@iDfI`A+`eh~%eTrJSEfR=pgs)b7$B&pj5UCVac_?{;B;zEvTIdiOX#1K_(28ZVs~eqqgJho}#wn z9N&T4;sb=CTdC*u(Bn`Ki8vN^!n%3S35)I$)Cuc*sJuhM@`y$F?rD7Q%a6OX_iCgEg;vL%ESC^LYV$hE zAPtm!eDezGYzKljq}D>=Z9~u~3gxv8ZIFhf7TayE(u+FX4D^Ex5jK1@86sdCfJ(cW z6w%^Ipe9d=rTfc`vJZzTQ2)a|j4JwlT zGzTkuA=2YaQZ*c{T-iig?f=W9C7+-*u$tRTrJa;Kg2I*v^$%h{pfD|gKk{zG)PwAe zh=bCI6!7PXG^Vunt$6qtM0ONkx85O<(K~NVWK&r2J5i|MRUn0mY@J`EQiUc*uU{J< zp1St#wJAc2eA5B=4_x=bc8#h|WtTAUaWyKt5yB<5U7vzBm3Vt?D)`v)U_U^fdJmK* zwzHV&GyC6WVm$&2UnSh&|M7|R0Ml4$+Q7O0sx;S+E`jwwU@(z$~;{u6yeFSbRT{L?-JZ=SV(cxR2F#D?}E)8EFeb+pl z84)(I)oqoa<$FKk?)2z*%7_rsHFeiX1*g?KZt51}7C5-2&l_-D?FR2p*a`4v+*{>b z2M?w?99Eo@+sW_tfQJX%HY%}`qV-Wd;NN~YE$B8boLus29*+lS#@y5j#fw)(|4#Ka ztSU%1Pw5>`26!z{#Uk8>``=awav%&jka$aV%vHw73eoBE58N3KIWZhf4!9bZCXcy1 ze;cPC1jIIl0@*`&5}=gK{TnuU2gqEK-=c)w2BcTH&H|PgJkN%y@?Ob-w{>W^|9VGK z$aSkEX)li;!!Hok&_5k1O>jEv5!Haf9@l%c9=O6?cX$|1=LG(JiiN_4VDRdovY%I4 zHqiMosP#D3@_G7hA3nRq)$CCyS=yMl_V-gQru{_HTPYR)(aT4Y(#OWfQ@s@Atz^;+ zdzA5yXDhBVhEU#(Y3-uY_b9oI#A!kVm97%vUsv$2vwW%0@8)ue`u&J@LV!!G52ht0 z9LM6Y-cwY4hLSN#$TsF%lN`8WA2|C2*c6kbF-nqXk~j@7f0LcZ#-1DD7;;_pXE91C vA4h`d!lpk&2(+LPjE0+nUExTuA-p4~H-AhC$3CW5@k{ZR_|AA=yy?FIYkDSi literal 0 HcmV?d00001 diff --git a/simulation/__pycache__/DualPortElements.cpython-39.pyc b/simulation/__pycache__/DualPortElements.cpython-39.pyc new file mode 100644 index 0000000000000000000000000000000000000000..091ec3c01b155d1e0593b0f1ff9c4852e325874e GIT binary patch literal 15803 zcmch8du&|UdEb4_`@tbOd{fj&s|QEYaz(DB-Ichjl|)h!tfd5zve#4DH={Y1fAdLIlU7*<_2k+$h%3JqKUNuwMbnr;08$27Wunx=pa0cqMe za2`sU_V<1F&hR0vqun^UgL~fhJihb2&zWK>6<6>~{USO+H79+B3H_eC5x336|^1zPYL0@fY%;I z;3&)LrpBS}7!TEls!;osrs68-Oq*>aZeFRcSh8FyR&M4)o>sFxW3}ikc@a^T)#8ff zMT)m-5cu;}a*x$Hz}TKX-Zl`sJG|YpdnbO}o5O zuN0kft$K5+UaVZJNoTrZtyopZ9$Q`a`b=}ta@xaRDc`zwZvK*KCTQ+7fTH%RlF(+y z5@~rfE-&#Xk$7^6O~nbhz@a{6Lw$R&sXF1!h>N8+wNC&q4E#29RjFpi75nw$ipsxF zB6KyO4Jb|Bi8?V?-%vNyDdkhSrokh`UEMVXl!vPPw4c)O#hh~I(@sL5?dO240I4zC z=-vV;Eof}nu2sF{)bz}`8&~E{Dry{>y_%bunxD#5%GKMsg<8ETa&^00UCaq0~E?=28FI=6SpMG;*62LvZEY^5G{UAnVJ*HV&vP!qxdo*?Vx@mTx z%FZBA)U29TN42P$!7r`J9+Z_dstF#NJx!VSY6%G!=vP-{)>Xa$`g~VwYSj}CN*#4W zlb2kes?Y)9mtFM>Kv8a^qfLhLLJvxerq)P1`ljJ(4^>t9gi=km_OPlQ$1c)7G`eyn zX#x8H)P-6_SklhTJ`S6Nls6x%gzWJ$KyEw5Y6(!Y9qC1_s<7{soh1gl7hNgdwm=4U zUiYHHx@!_W_)rY{UG4^0#x$>C<1SiNP^JS?AriEE5h#FCQjO!6(K4C_xP{fQCPz?H zA}F;Uf_odgzl#i@11ttk5^MmcP5q(rP-_}xr5SP!Hzc%4#npkQKimwv;Z1|!XO9y8 zyWtJ38F9jH#EEQZ_B^Mefb*x_h#PXllggb@E{mb;^OO?6Z*KcmJ4fxmxT}49If2qq zH?o<0sI+@ZxvHySmqPDqeTu-?Zx1)4PI|KgEus%K+NWK#_j@~NPg%}1(E@uHnSGo( z_{oDS2gUOTFCBd2;GAbTWpE(M#5wU<5QM;FW>3u*e01nofpMME^zZ$r3S9X z$vAddI=r=B1Wc#sI8ydg^Uwll#{u?N9nYwMW<0IB>P4!eC}nZo3p+I`lu5LeDf++= z^_5l6C|Q+CKIVmi7Uhy9_fS1;wilwlL3lK-S*lA}h*HT^xr6~IRfte{0X5z95pDK1 zk^=@^sVy#^t5nWWKc@LS>h1jq6qV@AP@}*k!lSGj*W&6h!larudetMEJb_jv=A6_X z5Ah7bE~dZNkO0;guClDShERbw5DrKY1}V!LQaVypVj!k`2r;;x0Is+O?~JFedFr~S z7QN8DvVfcjFInZqCC3ZPqA1rvD)J=8k*6ql2?6L#UoBVV%T(|R1*Z`dVyln_rT|oS zLm)wE$0Tc$!$lyNgB5#NfKP;Sqff9P>8M6s!v~8@{L-x z?9?PLZVqkRi{oYai__#e)XePIM%q$(c^kiqgsW^qyb<%`6th82f%t*)g284eAm@CbrA-;6L`c{t717;u~IqHFiMru zna1e3^K)0P+?by>uUtNFLK5X(&4Iy-(ZiMUEtBNsn6u^_&YzKefKVg*l;(#U+F^8e z^CrRKFyzLbk5wEy;PuiP6Ii9dl?HvR*yY7@Rscf6lx<4V`HE$!-mAl2&0}p1nN(X>?}i5Gddh{A){~tlp)g#OTrt)JY<< z67kj6L+LiUbsr8LOn}gzzH$P$k+&Ix?y3Y7yPlBDSMN9JPsEE#f;`BvDK6_EC~M zTBJ}*@Aeic2<7AcOhh{9q1_9FIVn9xyQV!y2-=-B9j$(wv0XKpc1?U&O=r6%ftrv5 z&frF%dO&?xj=~6n86>jz^_yBV*-W7&AOcrRgSwJ0<$>PO(~3-byH%%%`AWr zHhaXd8+AKez!b3ClZC1G1w-t8*XP1p?RwEepPMD?OI6-c?|fV2T>Y-{jv>=xUn_++ zw72_Py8zX0KglMMLrDZ?Nl%al?lpTci{}vcAwGb(-|fR}4l>8>Tkc28efTPd_fS@G zd$4zh@5AK57z1uERxz-%MRRabX$}Fh1GGE0jdyhUMF-G#2gCpoJZR}qa~NNDy1njj z07F0kJ?(LaHV553yfXccTh02dcGgf|=zRz~Onpr?hutA}81vra(rod4YldPJ-w#oL zS*1DR_H*rs$h*T|P~O$h>d1WqFx2n-mN?3g3(464hDhPt-7q}HFgzN-a4*(P=@ty1 z4`4X98-|;^VYru6SoHsF@H6DMKLY%WJ74*Gf?w!PRJ;(t@0)G-eG|}mA%Ij4n01_S ze4Me0b5CGauC;GyJwY%G_Aj?>|90ZyNc#)zov_R`dY$3TJpl)Q(bph$8+70`BU>$) z*Pj5dKPgU$ml&%fn|mLU`EK|zYZx!jvL5g_X-dn6MQx9*K+jjD7h0{|v!tiavl#V4 z4$NlT3s;IOw?xs??p3_dN=;aHAvAAREqpTVsh2#hF}_RhdqnzXUVD^9^oYzh5`zon zMZ4Uv7LmpC`)3zlMbFB&8wFCtA~?!p5Fp3Pp3GHHGUIbb=CjCO*~rl)ooqxv{lr%d3l?j-Y%D z>mnd7)Ff@BE>_F=9IH&8v2+sdBx457DS4EJkJxwWmesJl*j)@!lh_f~Ghm5S99W3; zTdS2Cydxd0lvxzZb=ZD=unpO0gO3W%u0?aDSS>CBTf+CsRZ+X=CGM@jX;oXX9Jwy@ zG!_|F}3%5g^@Canp;M_Xjvx13KtukzfEbvT7b;b1ww2$9wZNT0E(iGQth=-DPZXjy z8qHLSXhk*_{CT426yY32NLFYeO2ynXG;OZR#X|UAh2mZK?dddxG=!VT z1QaEqeLM4WG>CdLbq>()ZWoEYD&{dUa(1XK8Iro32})-?d~hcP$%R}hfy zt7%a^tNC{Ge%1dS&k>_U`$I(WI$pbg zz+%sM0p1`1Z<&CVzNm??h%BliO8)FcsHX3!A}LZR(aE1(Rz(NeW{`tV8}(i9sS7>G zbw5zatNm68K6V{ieJ}Yb*_63<_4@oA5Ymz5RryJbI{O_2c>etCcYglo==qDYkI237 zo&3$WfAf2zUq8#_|2R8W&W+BME9DYgj=AS^*NT<$Z7RJwksG~ODcW|s^nJr;SWNab z#3M2~{4(NgK3lYqnwdviE@H+fsFGmxz6!U4<eg9ujS-zGN)A$?t_tH;$2pnsDU>u0w3of^^Dr+ef8~s)EsiP%@ma7VRv{l&0;GJ zm1q{Kkha|6?iqj|jE(lMxO#J>3e&n(+8jYYo$d&f#2A#iR(Y3*-`6)&?y%pEg!P^G z(6U=3+1C?jg|+xEw@0 zgm^#V5yS@&4^P5qK)iQS0Zr!KQ=5m|V|eaiKqd#sOwKUM;$(adRV9uX&9`TJqOnvBRJ_~sDOM*|!o6+DjoG1NUDl#QV* zfifuGROa?0{X$T89Qh>j;{=!HDBe8YBf!rAIJ742t3)q1jpotjvF7v5G56R_14ii! z&Ep8in-k3!nhJg3=MagQK~&9xSJ)T+l9Suq z=k9f1dcqe@v}o90+Chu^m!I1_;GV(?PPwNx4~o;`3|t#XgTEZ|skHgB~)kK-c@!k=xuwu^C&yC)$R!n6taI|`>Ne;m<01~miXU0A^BD3 z=!VjqtPb)jj*$WVPRRbR*qIBwGuh^(^E}`*Mi>VEi zJLB$aZjUpuL3}wcrpbUr-3-ll14f7W-c4;(&mI5cYO^fzb?s{DCGjnkJ2Mmhd2$HyEvueF_A+2;MjOsrgBrWNfRLo!#buiU~> zBTkTl1nb?j=U#YG*6xv}?w^?0jSwGT$^971_gQNv=wfGSV2n*#bg8V5yGUF|f2dMNQ1tS!YK8K??;Taczz<^oNreiZc;$jk$F(9GAs1tYI8D**+@MM#2blq~l}tyOfVxYkd3!xbabvw7J!7Yu)D47udo? zS9HI^>mZYmvnIQ-rHi)m0R>;CfDOMSrBW1(Q2+`J4(*|Q!W@;bb;2nBBr4r;1XzfQ~!7$AQRZ3+o2J(wHpOxh6g z8+bZ4-^Wa1TAo&VlaKU}_$0%-kOI zH1gk31rsEhKE9DmUMUaPcFH(*G|ebHC^qZ>Md=Oe87+x(<1kK!=q3VeWQ|Uf!+IQ+ zsec-o!Tp4!c0@g)Wp#sWVEUmXqvj$}wO^0o?0Fa#@US|9USP|@#!l+6v-yPD5B>5? z1ouN)J2LIvw`KUrgd1@oI);~M#(YGg+p|k0`sRa zDb|y)qtiGOr{9SXtt-kaF|=d|J!7`8x-m)^)_bF8yOBHAf)9?~04qWs)AZt|38b#m zbtf-rnsDD#{2KXtZSd{{wEV}Yh!BiNk%+1qLq$Uj_+nie#zqd-zW)pVk=CvKFaD~8 z4|Y66^Yqx)5IkY;YjlqdHb3sL!NP~dFA}gHlCb$<>%-3P5b1?ZviBdTKc+0GqEmFe zl`7~z2HUSdr}*7oM*x89muD}U=WonCB0e=Me*z&J`|>9#b)SOIQ}8wgq>#vWDEKJ~ zewu=Rf?#G|{!>IH;#RF^FA{VNOPtXDuEs}*nxih$xpNl`>>em}It*;O_HzKHtbvlK zgUc+W*x~X?&vHcI%x5{uDV*3XN6AEUBT%Y);iV=8QI%ias&)0-hWtg05F_B+S2uB8 zivaXc=QcRy=ek&{E&Wr?k$`_2_kifO} zK!f9=69#hD?FI>^L4qCkd+72h$SgC$%gXbCdpUq1Se)|Jy_XaRT|z zC?G!lh}3*3C?Mi{^sf*=l_Ck{UpP61)ccy>f>Y4 z87o%4B6bpGo_-BCf&MkslAm7IsVrQSCEUaFG-x)o9_TvA=_TAm4cYZoOESmJ@0I_a zK4BBv|LX>g%vj z*w4QH2At)_`Wsl|S(>3%T8A|k!zD`cd*bpM3k(v-uLg zLO%YXNRf4K@Pa?Yr#{7W@1mHGN)x2iKa~S8cI?FIDTsQTjt`Y2X=2dA{%@mi=iaoD zQT_>a%-`=qw(p_hJoPk)1nF#P#_yKqg!s=u7Jn+K-46uObn=rkqNNzoIs!zae*kRP zoVvX3=^MZ=^TPG3pMY^OPswW^72BSMXgyd`TXFg_GEkPfp&$X-_y$%tr!zC1mM!yV~n}RaBx9OMR0!dI7#03LVu6)pMu>9cfwo z$hVNj{de*OWI#`Jmd^Clw-fW8fs}Aupg7M+O4Rm)Nr@D@XH|1x#rx6+s?Ogae`43C z=dRA4o4@*zbR0r}a_3>#?7qpqTgB-`y^hy@9l^7xIBv(H0(r1a!GU^{ZL2q1j2y`J zb@cf+7#7UP{|xQ&KLF*0$762IeEYX)_}`&`ZJ4JZ*1NRvq>3oXVIWlUuNz%$Xg>q) z;8;m{KwXfHLcXjhW1`QB12*h4VQTpV=5*4)cx2{8MkMv83#<) z(U%Q)e)fZ1XWcB=DejVTtb=V*9v&Z#$!`<(Yru9Rei^0M%I|qv4Ong%RzKa4g_lTv zEnv%Fujrdfp#br$=xpoo80&p!PFVdy)d?*NJV`HR%!ZmnmSg zn~&{m%907~_w$F8BB7Z7L2>Ty>kxn?jK9c{-0&Cle`xrpP}wgLi5x+t(h2iGPpSqC z-i$u14TE5O_i-Ff+ARLkM0yM}^6o9XeFj1d0gC^UM%j-+d`B$1UIUG&l?F-viI&#B zE2({6#eZtq>5|uHW@qywbR^m+6&c1z<4<4fzSDyEqkqY&S)uwcP5~iWC1$4hM>eGF`^Q};C^bpJISMXP@G%O0oC5MJvsKL&k0c#V zlJt>8VLbh9W)#~}WV3tGN#BM28d@nf4FRwFwjVs_TGWW9;+@gKcru!Z?}_SNpHbrJ R&nQmjQl>kzH#3^){Qn^B6S@EZ literal 0 HcmV?d00001 diff --git a/simulation/__pycache__/__init__.cpython-310.pyc b/simulation/__pycache__/__init__.cpython-310.pyc new file mode 100644 index 0000000000000000000000000000000000000000..67171b7ba59280f2da1bc77146a6d9398e872d7f GIT binary patch literal 515 zcmYk3%}(Pm5XYS+{r&*4yW+x;OAkd*Bn}86N)u?bRTCvCyRsr#A$5yITvx6mEsw+l z?ZwwU?JMjBJg8bRl7F5V|2-a!MYr2R#(KZsmdHZr4+NjBVuBm9wtxj%SOwx1ws-|s zxQ%V@V24+6mDg~M*KyqfO~q&$xWSva$y>O^+qiA=c8vPYU-Jg)*Y%$3J@{jfME>la z+|Fa&+a<=qJe+!w{^HHeM#fR_gZQ%q+)y_oFAip-B#3nD<7PM+#h&hc(n77yu2vXC zN!)kz*IhjPIV5*(;+lefq7MmKQ8g_&l`PMSDSHe@$%GL7jgb6_=Fk6+Tt71DgP8Hs z%~W>BR;8q2A=Iqj!TWBR6`_#ID=4RHwUsb3X#*TU6;K1z0S!PC0E3rpK*!(}$u9VS zf63x94WEP(YnG=;wxOa_>F{)TF*rRxKN|!a#U|eeF)I|6GE48>IOV@KEGJ7LS2W$4 iPK!)2u}%qLYo-V}*}Uk+6&H)Lpx5vTTWIbW9og^Qh>1o3 literal 0 HcmV?d00001 diff --git a/simulation/__pycache__/__init__.cpython-312.pyc b/simulation/__pycache__/__init__.cpython-312.pyc new file mode 100644 index 0000000000000000000000000000000000000000..85b6d196a7404a3d4499faad4bbdfb9c5617f9c8 GIT binary patch literal 548 zcmYk3&u`N(6vyo({oRgffVgt)fohGIG$Dj&o0LWBnkY%xWLdI8YDOS#CMS{DU&7zO z-@+9xzH;IQ+M&`e@LUzdJAC}!^Yhow&tC?EE&>~T`~4PSgud(GuNiyi_9Zx<&>F3A zjszxH7{m}JF@;4e(I5@cBu&vGEv#z{sA-e7=#Y-+lCJ2H9!59j(E0)TXxLKs9QU<1 zb6MyuU$N;bQV09Q%wGiyH&plC6)bEP`WMVw#@Y>3Cv+o!If?yHb$2I&vq|Ktga2x^ zR)4P*_+cCkE%jg@+m|+*JFx?TUaTH5cEzhjEx63mBwu82f=PVFn0m-q`i`d`{%tv4 z2pJHziepg9`i&D+|U47(LS%FKb{DKL#Q9B7ElQzHs)Bz1Z6VL*<02sXN0J;VbNcO-7 zyx{y>idUs7H+-2U*_M@+N{1K2uY-%r%P)g)tN8TiAm)W)QfBFm7pLWRwa(eUbZ5FO fGR4bHN-5uPMd|tWUbnuBGOr4D1s|}37SHTA!8(WP literal 0 HcmV?d00001 diff --git a/simulation/__pycache__/__init__.cpython-39.pyc b/simulation/__pycache__/__init__.cpython-39.pyc new file mode 100644 index 0000000000000000000000000000000000000000..f6fab24102ad4b1534887d5c2d4d4f68ce12e9b7 GIT binary patch literal 475 zcmYk2%}T>S5XU!7+N9}+;s?HfmmDlM(2Iyzo74tM6E>;J5?Erp)?l_RX#)C4K9HBa zdh!*#=#CWW!2WhV{yRHN(rTGVTKnVY2`LDDLa=KU8C=V=0nE`{$q-kF!Yic0Ribi@ zXuL|QyhdtVC%OWfiqsgy;B`{xCNX(~G-O^4k)?gfKai!1J;%NC1}??kvR@x