diff --git a/book/pages/scatparamtest.py b/book/pages/scatparamtest.py new file mode 100644 index 0000000..ae7c7ea --- /dev/null +++ b/book/pages/scatparamtest.py @@ -0,0 +1,151 @@ +# Imports +import meep as mp +import numpy as np +import matplotlib.pyplot as plt +import os +from pathlib import Path +from gplugins.gmeep.get_meep_geometry import get_meep_geometry_from_component +from gdsfactory.read import import_gds +import gdsfactory as gf +import pathlib +mp.verbosity(3) + +res = 20 # the resolution of the simulation in pixels/um +sim_is_3D = False # Turn this to false for a 2D simulation + +pwd = pathlib.Path(__file__).parent.resolve() +gds_file = pwd.parent / "files/mmi2x2.gds" # The name of our gds file + +# The Parameters for the frequencies we'll be using +lcen = 1.55 # Center wavelength +fcen = 1 / lcen # Center frequency +df = 0.05*fcen # Frequency Diameter + +# The thickness of each material in our simulation (only used in 3D simulations) +t_oxide = 1.0 +t_Si = 0.22 +t_air = 0.78 + +dpml = 1 # Diameter of perfectly matched layers +cell_thickness = dpml + t_oxide + t_Si + t_air + dpml # Cell thickness + +# Materials used in the simulation +Si = mp.Medium(index=3.45) +SiO2 = mp.Medium(index=1.444) + +# Sets the min and max values for the cell and the silicon. Our simulation will be centered at y=0 +cell_zmax = 0.5*cell_thickness if sim_is_3D else 0 +cell_zmin = -0.5 * cell_thickness if sim_is_3D else 0 + +# Create a 2D array to hold the S-Parameters for the device +n_ports = 4 # The number of ports, also the size of our array +s_params = np.zeros((n_ports, n_ports)) + +mode_parity = mp.NO_PARITY if sim_is_3D else mp.EVEN_Y + mp.ODD_Z + +from gdsfactory.technology import LayerLevel, LayerStack +layers = dict(core=LayerLevel( + layer=(1,0), + thickness=t_Si, + zmin=-t_Si/2, + material="si", + mesh_order=2, + sidewall_angle=0, + width_to_z=0.5, + orientation="100",) + ) +layer_stack = LayerStack(layers=layers) + +mmi_comp = import_gds(gds_file) +geometry = get_meep_geometry_from_component(mmi_comp, is_3d=sim_is_3D, wavelength=lcen, layer_stack=layer_stack) +# Use this to modify the material of the loaded geometry if needed. +# geometry = [mp.Prism(geom.vertices, geom.height, geom.axis, geom.center, material=mp.Medium(index=3.45)) for geom in geometry] + +# ################################################### +# Now we actually import the geometry +cell_x = 32 +cell_y = 6 +cell_z = 3 + +port_xsize = 0 +port_ysize = 1 +port_zsize = 0.5 + +port_xdisp = 13 +port_ydisp = (0.75+0.5)/2 + +port_size = mp.Vector3(port_xsize, port_ysize, port_zsize) if sim_is_3D else mp.Vector3(port_xsize, port_ysize, 0) +cell = mp.Vector3(cell_x, cell_y, cell_z) if sim_is_3D else mp.Vector3(cell_x, cell_y, 0) + +port1 = mp.Volume(center=mp.Vector3(-port_xdisp,-port_ydisp,0), size=port_size) +port2 = mp.Volume(center=mp.Vector3(-port_xdisp,port_ydisp,0), size=port_size) +port3 = mp.Volume(center=mp.Vector3(port_xdisp,port_ydisp,0), size=port_size) +port4 = mp.Volume(center=mp.Vector3(port_xdisp,-port_ydisp,0), size=port_size) +source1 = mp.Volume(center=port1.center-mp.Vector3(x=0.5),size=port_size) +source2 = mp.Volume(center=port4.center+mp.Vector3(x=0.5), size=port_size) + +subtraction_geom = [mp.Block(size=mp.Vector3(mp.inf, 0.5, t_Si), material=Si)] +subtraction_cell = mp.Vector3(5, 4, cell_z if sim_is_3D else 0) +subtraction_sources = [ + mp.EigenModeSource( + src = mp.GaussianSource(fcen, fwidth=df), + volume=mp.Volume(center=mp.Vector3(-0.5,0,0), size=port_size), + eig_band=1, + eig_parity = mode_parity, + eig_match_freq = True, + ) +] + +subtraction_sim = mp.Simulation( + resolution=res, + cell_size=subtraction_cell, + boundary_layers=[mp.PML(dpml)], + sources=subtraction_sources, + geometry=subtraction_geom, + default_material=SiO2 +) + +subtraction_monitor_region = mp.FluxRegion(center=mp.Vector3(0), size=port_size) +subtraction_monitor = subtraction_sim.add_flux(fcen, 0, 1, subtraction_monitor_region) + +plot_plane = mp.Volume(center=mp.Vector3(z=0), size=mp.Vector3(cell.x, cell.y, 0)) +subtraction_sim.plot2D(output_plane=plot_plane if sim_is_3D else None) +subtraction_sim.run(until_after_sources=mp.stop_when_dft_decayed) + +subtraction_data = subtraction_sim.get_flux_data(subtraction_monitor) +norm_mode_coeff = subtraction_sim.get_eigenmode_coefficients(subtraction_monitor, [1], mode_parity).alpha[0,0,0] + +# Set up the first source for the simulation. I'll start with port1 (the lower left) +sources = [ + mp.EigenModeSource( + src = mp.GaussianSource(fcen, fwidth=df), + volume=source1, + eig_band=1, + eig_parity = mode_parity, + eig_match_freq = True, + ) +] + +# Create Simulation +sim = mp.Simulation( + resolution=res, # The resolution, defined further up + cell_size=cell, # The cell size, taken from the gds + boundary_layers=[mp.PML(dpml)], # the perfectly matched layers, with a diameter as defined above + sources = sources, # The source(s) we just defined + geometry = geometry, # The geometry, from above + default_material=SiO2 +) + +# Adds mode monitors at each of the ports to track the energy that goes in or out + +mode_monitor_1 = sim.add_mode_monitor(fcen, 0, 1, mp.ModeRegion(volume=port1)) +mode_monitor_2 = sim.add_mode_monitor(fcen, 0, 1, mp.ModeRegion(volume=port2)) +mode_monitor_3 = sim.add_mode_monitor(fcen, 0, 1, mp.ModeRegion(volume=port3)) +mode_monitor_4 = sim.add_mode_monitor(fcen, 0, 1, mp.ModeRegion(volume=port4)) + +sim.load_minus_flux_data(mode_monitor_1, subtraction_data) + +# Plot the simulation +plot_plane = mp.Volume(center=mp.Vector3(z=0), size=mp.Vector3(cell.x, cell.y, 0)) +# sim.plot2D(output_plane=plot_plane if sim_is_3D else None) # No parameters are needed for a 2D simulation. + diff --git a/book/pages/scattering_parameters.ipynb b/book/pages/scattering_parameters.ipynb index 49bdaec..41c17b4 100644 --- a/book/pages/scattering_parameters.ipynb +++ b/book/pages/scattering_parameters.ipynb @@ -76,15 +76,18 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 6, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "\u001b[32m2024-02-14 15:02:04.135\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mgplugins.gmeep\u001b[0m:\u001b[36m\u001b[0m:\u001b[36m39\u001b[0m - \u001b[1mMeep '1.28.0' installed at ['/home/andeloth/miniconda3/envs/photonics/lib/python3.11/site-packages/meep']\u001b[0m\n" - ] + "data": { + "text/plain": [ + "3" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ @@ -96,7 +99,8 @@ "from pathlib import Path\n", "from gplugins.gmeep.get_meep_geometry import get_meep_geometry_from_component\n", "from gdsfactory.read import import_gds\n", - "import gdsfactory as gf" + "import gdsfactory as gf\n", + "mp.verbosity(1)" ] }, { @@ -109,7 +113,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -157,17 +161,9 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 8, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\u001b[32m2024-02-14 15:02:04.266\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mgdsfactory.pdk\u001b[0m:\u001b[36mactivate\u001b[0m:\u001b[36m309\u001b[0m - \u001b[1m'generic' PDK is now active\u001b[0m\n" - ] - } - ], + "outputs": [], "source": [ "from gdsfactory.technology import LayerLevel, LayerStack\n", "layers = dict(core=LayerLevel(\n", @@ -190,7 +186,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ @@ -227,67 +223,16 @@ "\n", "In order to find all the S-Parameters for a device, we have to run a number of simulations where we set the source to be in one of the four ports and then record the output field in all ports. Normally we would have to do this as many times as number of ports, setting a different one to be the input in each time, and we will demonstrate how we can do this in this notebook.\n", "\n", - "Note, however, that the MMI we are using has mirror symmetry about both the x and y axes, which means we could just do the simulation once and then assign the s-parameters using symmetry arguments. This is a useful trick when working with devices with symmetries about the ports that can save a lot of time.\n", - "\n", - "Before we run the simulations, we will do a different simulation of a straight waveguide and record the fields going through a monitor. This simulation will simply include the section between the source and the reflection monitor. We will need to provide this field for the reflection monitors ($s_{11}, s_{22}$ etc.) so we can get the reflected field accurately." + "Note, however, that the MMI we are using has mirror symmetry about both the x and y axes, which means we could just do the simulation once and then assign the s-parameters using symmetry arguments. This is a useful trick when working with devices with symmetries about the ports that can save a lot of time." ] }, { "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiEAAAG2CAYAAACpnFbhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABT90lEQVR4nO3de3yU1bkv8N/kNhAMUa5JFCV1V2BDiwKioICIRsRr6/F4K8ZWVBB0a+p2gxcIIqLU2h43iqKAeKl6ToGtFUSpQKKCghAqKCDK1ZBIAE0iSgLJOn+8TNYkmcv7zqz3/vt+PvlAwkxmMc9az/PLZJIJCCEEiIiIiCyWYvcCiIiIyJ8YQoiIiMgWDCFERERkC4YQIiIisgVDCBEREdmCIYSIiIhswRBCREREtmAIISIiIlswhBAREZEtGEKIiIjIFq4JITNmzMDZZ5+NrKwsdOnSBVdffTW2bdsW93olJSXo378/2rRpg1/84hd47rnnLFgtERERxeOaEFJSUoLx48fjk08+wfLly3Hs2DEUFBTg8OHDUa+zc+dOjBo1CkOGDEFZWRkeeOAB3H333Vi4cKGFKyciIqJIAm59Abuqqip06dIFJSUlGDp0aMTL/Nd//RfefvttbNmypeljY8eOxb/+9S+sWbPGqqUSERFRBGl2LyBR1dXVAIAOHTpEvcyaNWtQUFDQ7GOXXHIJ5s6di6NHjyI9Pb3Vderq6lBXV9f0fmNjIw4dOoSOHTsiEAgoWj0REZE3CSFQW1uLvLw8pKTE/oaLK0OIEAJFRUU4//zz0adPn6iXq6ysRNeuXZt9rGvXrjh27BgOHDiA3NzcVteZMWMGpk6dqnzNREREfrJ3716ccsopMS/jyhAyYcIEfP755/joo4/iXrbloxeh7z5Fe1Rj0qRJKCoqanq/uroap556Kka/NBpDzhgS9Xae+vgpbD24FT079kTReUVRL2emj/d8jAUbFyCYFsS04dNwYtsTLV/DDz//gIdXPoy6Y3UoPLMQ5516nqHr7/5hN6aXTsecK+egR8ceCa/j5X+9jLkb5uLWfrfi5r43J/x5krGlaguK3itC/kn5ePLiJ5GZkWn5Gn6q/wn3Lb8PO7/fiacueQq9OveyfA0A6xHCekhW1mPbwW24/e3b8eDQB3Haiac1fTzZfqWK1+bHh199iFdueQVZWVnxLyxcZsKECeKUU04RO3bsiHvZIUOGiLvvvrvZxxYtWiTS0tJEfX29rturrq4WAMSc1XOiXmbEghECxRAjFozQ9TnNMG/DPIFiiKzHskR5dbktayivLhdZj2UJFEPM2zAvoc+xft96gWKI9fvWJ7yOaSXTBIohppVMS/hzJOvTbz8V7We0F4PnDhY1R2psWUPNkRoxeO5g0X5Ge/Hpt5/asgYhWI8Q1kOyuh6R+oqKfqWCF+fHnNVzBABRXV0d97KuCSGNjY1i/PjxIi8vT3z11Ve6rnP//feLXr16NfvY2LFjxbnnnqv7duOFEC9uoESoOtDJhhA/NthIOPAk1kPyaz1a9hUGEMmM+eHJEDJu3DiRnZ0tVq1aJSoqKprefvrpp6bLTJw4UYwePbrp/R07dojMzExx7733ii+//FLMnTtXpKeni7///e+6bzdWCPHqBjJK5YFOJoT4tcG2xIEnsR6Sn+sR3lcYQCSz5ocnQwiAiG/z589vukxhYaEYNmxYs+utWrVKnHXWWSIjI0N0795dzJ4929DtRgshXt5ARqg+0ImGED832HAceBLrIfm9HqG+smz7MgaQ48ycH54MIXaJFEK8voH0MuMrikRCiN8bbAgHnsR6SKyH7CuZ0zMZQIT584MhRKGWIcQPG0gPsx7SNBpC2GA1HHgS6yGxHprSXaUCxWAAEdbMD4YQhcJDiF82UDxmfk/VSAhhg9Vw4Emsh8R6aGqO1Ii+s/sKFENMXjnZljUI4Z8AIgRDiFKhENLzTz19s4FiMftJXXpDCBushgNPYj0k1kMTqke76e2S/tH/ZPgpgAjBEKJUKIRgon82UDRWPKtcTwhhg9Vw4Emsh8R6aMLrsWDjAttCiN8CiBAMIUqFPxJiF78EECHihxA2WA0HnsR6SKyHpmU9VPwSxET4MYAIwRCilJ7fmGomPwUQIWKHEDZYDQeexHpIrIcmUj3sCCF+DSBCMIQoZWcI8VsAESJ6s2CD1XDgSayHxHpootXD6hDi5wAiBEOIUnaFED8GECEiNws2WA0HnsR6SKyHJlY9rAwhfg8gQjCEKGVHCLF7Awlh32srtGwWbLAaDjyJ9ZBYD028elgVQhhANDM/mMkQoorVIcQJG8jO11YIbxZssBoOPIn1kFgPjZ56WBFCGEA05dXlIvhQkCFEFStDiFM2kJ2vrRBqFuPeGccGKzjwwrEeEuuh0VsPs0MIA4imaX5MBEOIKlaFEEdtIBt/tXGoWbDBcuCFYz0k1kNjpB5mhhAGEE34/Ch8o5AhRBUrQojTNpCdr60wcM7ApkdC7OK2BmsmDjwN6yG5sR5mhRAGEE3L+cEnpipkdghx4gayS+hA2/nrld3YYM3CgadhPSS31sOMEMIAook0P4yEkBSQbeaXzccf3v4DsjKysHX8VuS1z7N8Dftq9qHnMz1RW1+LeVfOw+/P+r3lawCAi16+CB/s/AAD8wbacvsAsLZ8LS5+5WL06dIHy25ahqxgluVrqK2rxcjXRmLz/s1YPno5Bp5sz/3xaOmjeHjlw5g2fBoeGvqQLWtgPSTWQ+OUeoT61Yj8Efjnzf+0ZQ2emR9mpyS3M+uREKcmWDuEf0Vh169XdutXeGbgV9wa1kNyez1U9hU+AqKJNT/47RiFzAghTt9AVmp5oO0IIW5vsCpx4GlYD8kL9VDVVxhANPHmB0OIQqpDiBs2kFUiHWirQ4gXGqwqHHga1kPySj1U9BUGEI2e+cEQopDKEOKWDWSFaAfayhDilQarAgeehvWQvFSPZPsKA4hG7/xgCFFIVQhx0wYyW6wDbVUI8VKDTRYHnob1kLxWj2T6CgOIxsj8YAhRSEUIcdsGMlO8A21FCPFag00GB56G9ZC8WI9E+woDiMbo/GAIUSjZEOLGDWQWPQfa7BDixQabKA48DeshebUeifQVBhBNIvODIUShZEKIWzeQGfQeaDNDiFcbbCI48DSsh+TlehjtKwwgmkTnB0OIQomGEDdvINWMHGizQoiXG6xRHHga1kPyej2M9BUGEE0y84MhRKFEQojbN5BKRg+0GSHE6w3WCA48Desh+aEeevsKA4gm2fnBEKKQ0RDihQ2kSiIHWnUI8UOD1YsDT8N6SH6ph56+wgCiUTE/GEIUMhJCvLKBVEj0QKsMIX5psHpw4GlYD8lP9YjXVxhANKrmB0OIQnpDiJc2ULKSOdCqQoifGmw8HHga1kPyWz1i9RUGEI3K+cEQopCeEOK1DZSMZA+0ihDitwYbCweehvWQ/FiPaH2FAUSjen4whCgUL4R4cQMlSsWBTjaE+LHBRsOBp2E9JL/WI1JfYQDRmDE/GEIUihVCvLqBEqHqQCcTQvzaYCPhwNOwHpKf69GyrzCAaMyaHwwhCkULIV7eQEapPNCJhhA/N9iWOPA0rIfk93qE9xUGEI2Z84MhRKFIIcTrG8gI1Qc6kRDi9wYbjgNPw3pIrIfsKwPnDGQAEebPD4YQhVqGED9sIL3M+IrCaAhhg5U48DSsh8R6aEJ9hQHEmvnBEKJQeAjxywbSw6yHNI2EEDZYiQNPw3pIrIc07p1xTY+E2MVP88OzIaSkpERcfvnlIjc3VwAQixcvjnn5lStXCgCt3rZs2aL7NkMhpPCNQt9soHjM/J6q3hDCBitx4GlYD4n1kEL1MOM1qfTyUwARwlgISYGLHD58GH379sWsWbMMXW/btm2oqKhoevvlL39p+LYXbFyArIwsbB2/FXnt8wxfP1n7avah5zM9UVtfi3lXzsPvz/q95WsAgItevggf7PwAI/JH4J83/9OWNawtX4uLX7kYfbr0wbKbliErmGX5GmrrajHytZHYvH8zlo9ejoEnD7R8DQDwaOmjeHjlw5g2fBoeGvqQLWtgPSTWQ+O0eowbMM6W2weA+WXz8Ye3/8D5EY2pcchEMPBIyPfff5/w7YQeCQk+FPRFgo3FimeVx3skhF/hSfyKW8N6SKyHFF4Ps16dOx6/PQIS4tlvx4QzEkK6d+8ucnJyxIUXXihWrFgR8zpHjhwR1dXVTW979+4VAMTMD2YqXL1+fgogQsQOIWywEgeehvWQWA+pZT3sCCF+DSBCMIQ02bp1q5gzZ45Yv369WL16tRg3bpwIBAKipKQk6nWmTJkS8Xkkel9FVyW/BRAhojcLNliJA0/DekishxSpHlaHED8HECEYQmK6/PLLxRVXXBH136M9EmJ1CPFjABEicrNgg5U48DSsh8R6SNHqYWUI8XsAEcLDT0xV4dxzz8X27duj/nswGET79u2bvVnNKU8i4pNQNU57kh2f9Mh6hLAekhPqwSehJsCCUGQKJPhIyDXXXCOGDx+u+/J6XkVXJbsTbIhdv9o4/CsWfoUn8StuDeshsR5SvHpY8UgIHwGRev6ppze/HVNbWyvKyspEWVmZACCeeuopUVZWJnbv3i2EEGLixIli9OjRTZf/y1/+IhYvXiy++uorsXnzZjFx4kQBQCxcuFD3bVoZQpyygex8bYVQs1iwcQEb7HEceBrWQ2I9JD31MDuEMIBIIxaMEJgIb4aQaL98rLCwUAghRGFhoRg2bFjT5Z944glx+umnizZt2oiTTjpJnH/++WLJkiWGbtOqEOKoDWTjrzYONYt209uxwQoOvBDWQ2I9JL31MDOEMIBIofnh2UdC7GBFCHHaBrLztRUmr5wsUAzRd3ZfNlgOPCEE6xGO9ZCM1MOsEMIAIoXPD1/8dIxVzA4hTtxAdgkdaBRDlO4qtWUNbmywZuHAk1gPjVvrYUYIYQCRWs4PhhCFzAwhTt1Adggd6Mzpmba9xoNbG6wZOPAk1kPj5nqoDiEMIFKk+cEQopBZIcTJG8hq4Qd62fZltoQQNzdY1TjwJNZD4/Z6qAwhDCBStPnBEKKQGSHE6RvISi0PtB2/XtntDVYlDjyJ9dB4oR6q+goDiBRrfjCEKKQ6hLhhA1kl0oG2OoR4ocGqwoEnsR4ar9RDRV9hAJHizQ+GEIVUhhC3bCArRDvQVoYQrzRYFTjwJNZD46V6JNtXGEAkPfODIUQhVSHETRvIbLEOtFUhxEsNNlkceBLrofFaPZLpKwwgkt75wRCikIoQ4rYNZKZ4B9qKEOK1BpsMDjyJ9dB4sR6J9hUGEMnI/GAIUSjZEOLGDWQWPQfa7BDixQabKA48ifXQeLUeifQVBhDJ6PxgCFEomRDi1g1kBr0H2swQ4tUGmwgOPIn10Hi5Hkb7CgOIlMj8YAhRKNEQ4uYNpJqRA21WCPFygzWKA09iPTRer4eRvsIAIiU6PxhCFEokhLh9A6lk9ECbEUK83mCN4MCTWA+NH+qht68wgEjJzA+GEIWMhhAvbCBVEjnQqkOIHxqsXhx4Euuh8Us99PQVBhAp2fnBEKKQkRDilQ2kQqIHWmUI8UuD1YMDT2I9NH6qR7y+wgAiqZgfDCEK6Q0hXtpAyUrmQKsKIX5qsPFw4Emsh8Zv9YjVVxhAJFXzgyFEIT0hxGsbKBnJHmgVIcRvDTYWDjyJ9dD4sR7R+goDiKRyfjCEKBQvhHhxAyVKxYFONoT4scFGw4EnsR4av9YjUl9hAJFUzw+GEIVihRCvbqBEqDrQyYQQvzbYSDjwJNZD4+d6tOwrDCCSGfODIUShaCHEyxvIKJUHOtEQ4ucG2xIHnsR6aPxej/C+wgAimTU/GEIUihRCvL6BjFB9oBMJIX5vsOE48CTWQ8N6yL4yeeVkBpDjzJwfDCEKtQwhfthAepnxFYXREMIGK3HgSayHhvXQhPoKA4jG7PnBEKJQeAjxywbSw6yHNI2EEDZYiQNPYj00rIe0YOMCgWKIzOmZDCAWzA+GEIVCIWTmBzN9s4HiMfN7qnpDCBusxIEnsR4a1kP69NtPRbvp7QSKIZZtX2bLGvwUQIRgCFEqFEKCDwV9s4FiMftJXXpCCBusxIEnsR4a1kMK1aPv7L5KXw7CCL8FECEYQpQKhRBM9M8GisaKZ5XHCyFssBIHnsR6aFgPKbwepbtKbQkhfgwgQjCEKBUKIYVvFNq2Br8EECFihxA2WIkDT2I9NKyH1LIeZrw6dzx+DSBCMIQoZfRVdFXzUwARInoIYYOVOPAk1kPDekiR6mF1CPFzABGCIUQpO0OI3wKIEJGbBRusxIEnsR4a1kOKVg8rQ4jfA4gQDCFK2RVC/BhAhGjdLNhgJQ48ifXQsB5SrHpYFUIYQDQMIQrZEULs3kBC2PfaCuHNgg1W4sCTWA8N6yHFq4cVIYQBRCp8o5AhRBWrQ4gTNpCdr60Qahalu0rZYI/jwJNYDw3rIemph9khhAFEmrdhnsBEMISoYmUIccwGsvFXG4eaRd/ZfdlgBQdeONZDw3pIeuthZghhAJFC8yP4UJAhRBWrQoiTNpCdr62wbPsygWKIdtPbscFy4DVhPTSsh2SkHmaFEAYQKXx+zPxgJkOIKlaEEKdtIDtfWyFzeqZAMcSCjQtsWYMQ7muwZuHAk1gPyY31MCOEMIBILecHn5iqkNkhxIkbyA7hB9quX68shDsbrBk48CTWQ3JrPVSHEAYQKdL88GwIKSkpEZdffrnIzc0VAMTixYvjXmfVqlWiX79+IhgMivz8fDF79mxDt2lmCHHqBrJa+IGevHKybSHErQ1WNQ48ifWQ3FwPlSGEAUSKNj88G0KWLl0qHnzwQbFw4UJdIWTHjh0iMzNT/Md//If48ssvxQsvvCDS09PF3//+d923aVYIcfIGslLLA23Hr1cWwt0NViUOPIn1kNxeD1V9hQFEijU/PBtCwukJIffff7/o2bNns4/dcccd4txzz9V9O2aEEKdvIKtEOtB2hBC3N1hVOPAk1kPyQj1U9BUGECne/GAIOW7IkCHi7rvvbvaxRYsWibS0NFFfX6/rdlSHEDdsICtEO9BWhxAvNFgVOPAk1kPySj2S7SsMIJKe+WEkhKTBwyorK9G1a9dmH+vatSuOHTuGAwcOIDc3t9V16urqUFdX1/R+TU0NAGD3D7uxoWJDUusZ+85YrNu3DmfnnY2ZF89M+vMl4u2tb2Nq6VRkpmXi//2v/4fKw5WoPFxp6RqqfqzCb//vb/HTsZ8wZegU9M3p23RfbDmwpdmfZnpxw4uY/dlsjBswDqN+OcqWemzevxl3LrkTp3c4HY+PeBzbD223fA2H6w9jwrsT8M2hb/DsZc8iLSXNlvuC9dCwHpKqeiTTV2L1Kyu5aX7s/mG37s8ZEEIIlYu0SiAQwOLFi3H11VdHvcwZZ5yB3//+95g0aVLTxz7++GOcf/75qKioQE5OTqvrFBcXY+rUqa0/2UQAbRQsnIiIyMuOAHgcqK6uRvv27WNe1NOPhOTk5KCysnlK279/P9LS0tCxY8eI15k0aRKKioqa3q+pqUG3bt0w58o56N+9f0LrCP+KYky/MQl9jmSFf0Ux69JZaJfRzvI1tPwKr0+XPq0us+XAFvxu0e/w6m9fRa9OvSxfIxF5TyJ9RU+/soIb58f6Xetx++O36/rcng4hgwYNwj/+8Y9mH3v//fcxYMAApKenR7xOMBhEMBhs9fEeHXugX24/w2t4tPRRzP5sNqYNn4aHhj5k+PoqrC1fi7vevQt9c/pi2U3LkBXMsnwNtXW1GPnaSOz6YRdWFK7AwJMHxrx8r069Erq/iYii0dtXjPYrs7h1fvxY+6Puz5+SzOKs9uOPP2Ljxo3YuHEjAGDnzp3YuHEj9uzZA0B7FOPmm29uuvzYsWOxe/duFBUVYcuWLZg3bx7mzp2L++67z5L1Plr6KB5e+bDtG+jiVy5Gny59bA8gm/dvxvLRy2070ERE8TilX/llfrjqkZDPPvsMw4cPb3o/9G2TwsJCvPTSS6ioqGgKJACQn5+PpUuX4t5778UzzzyDvLw8PP3007jmmmtMX6tfNlA8TjnQRETxOKVf+Wl+uCqEXHDBBYj1PNqXXnqp1ceGDRuGDRusfRaxnzZQLE450ERE8TilX/ltfrjq2zFu4LcNFI1TDjQRUTxO6Vd+nB8MIQr5cQNF4pQDTUQUj1P6lV/nB0OIIn7dQC055UATEcXjlH7l5/nBEKKAnzdQOKccaCKieJzSr/w+PxhCkuT3DRTilANNRBSPU/oV5wdDSFK4gTROOdBERPE4pV9xfmgYQhLEDaRxyoEmIorncP1hR/Qrzg/JVb8nxCm4gTQMIETkJhPenYBdP+xiAHHA/AhhCDGIG0jDAEJEbvPNoW9sfy0Yzo/m+O0YA7iBNAwgRORGz172LAOIgwIIwEdCdHv5Xy9j7pdzfb+BGECIyK36dOljy+0ygETHR0J0mruBAYQBhIjIGAaQ2BhCdLq1362+3kAMIERExjCAxMcQotPNfW+25XadsIEYQIiIjGEA0YchxMGcsIEYQIiIjGEA0Y8hxKGcsIEYQIiIjGEAMYYhxIGcsIEYQIiIjGEA0fxU/5PuyzKEOIwTNhADCBGRMQwgmtq6Wty3/D7dl2cIcRCnbCAGECIi/RhANKH5sfP7nbqvwxDiEE7aQAwgRET6MIBowufHU5c8pft6/I2pDuC0DWRnAHlxw4u23C4RkVEMIJqW8+PI4SO6r8tHQmzmxA1k52srzP5sti23TURkBAOIJtn5wRBiIy9sIFVCB3rcgHG23D4RkV4MIBoV84MhxCZe2UAqhB/oMf3G2LIGIiI9GEA0quYHQ4gNvLSBkuWEA01EpIcT+pXX5gdDiMW8toGS4YQDTUSkhxP6lRfnB0OIhby4gRLlhANNRKSHE/qVV+cHQ4hFvLqBEuGEA01EpIcT+pWX5wdDiAW8vIGMcsKBJiLSwwn9yuvzgyHEZF7fQEY44UATEenhhH7lh/nBEGIiP2wgvZxwoImI9HBCv/LL/GAIMYlfNpAeTjjQRER6OKFf+Wl+MISYwE8bKB4nHGgiIj2c0K/8Nj8YQhTz2waKxQkHmohIDyf0Kz/OD4YQhfy4gaJxwoEmItLDCf3Kr/ODIUQRv26gSJxwoImI9HBCv/Lz/HBdCHn22WeRn5+PNm3aoH///vjwww+jXnbVqlUIBAKt3rZu3ap0TX7eQC054UATEenhhH7l9/nhqhDy5ptv4p577sGDDz6IsrIyDBkyBJdeein27NkT83rbtm1DRUVF09svf/lLZWvy+wYK54QDTUSkhxP6FeeHy0LIU089hVtvvRVjxoxBr1698Ne//hXdunXD7NmzY16vS5cuyMnJaXpLTU1Vsh5uIMkJB5qISA8n9CvOD41rQkh9fT3Wr1+PgoKCZh8vKCjA6tWrY173rLPOQm5uLkaMGIGVK1fGvGxdXR1qamqavUXCDSQ54UBTEgYMAE45RfvTRVy6bLKZE/oV54fkmhBy4MABNDQ0oGvXrs0+3rVrV1RWVka8Tm5uLubMmYOFCxdi0aJF6NGjB0aMGIHS0tKotzNjxgxkZ2c3vXXr1q3VZbiBJCccaEpSZSVQXq796SIuXTbZ6MUNL9rerzg/mkuz7ZYTFAgEmr0vhGj1sZAePXqgR48eTe8PGjQIe/fuxZNPPomhQ4dGvM6kSZNQVFTU9H5NTU2zIMINJDGAEJGbzP5sNgOIQ+ZHiGseCenUqRNSU1NbPeqxf//+Vo+OxHLuuedi+/btUf89GAyiffv2zd5CuIEkBhAicptxA8YxgDhgfoRzTQjJyMhA//79sXz58mYfX758OQYPHqz785SVlSE3N9fw7W+p2sINdBwDCBG50Zh+Y2y5XQaQ6Fz17ZiioiKMHj0aAwYMwKBBgzBnzhzs2bMHY8eOBaB9K6W8vBwvv/wyAOCvf/0runfvjt69e6O+vh6vvvoqFi5ciIULFxq/7feKcGb3M32/gRhAiIj0YwCJzVUh5LrrrsPBgwfxyCOPoKKiAn369MHSpUtx2mmnAQAqKiqa/c6Q+vp63HfffSgvL0fbtm3Ru3dvLFmyBKNGjTJ82/kn5ft+AzGAEBHpxwASn6tCCADceeeduPPOOyP+20svvdTs/fvvvx/333+/ktt98uInfb2BGECIiPRjANHHNc8JsVtmRqblt+mUDcQAQkSkHwOIfgwhDuWUDcQAQkSkHwOIMQwhDuSUDcQAQkSkHwOIca57TohdGhsb0djYaPrt1NbVYtTfRmFz1Wa8d9N7GJA7wJLbbWn6h9MxedVkPHLBI3jg/AcsW0Podqy6vwkIHH8TAISr7nO58sZGYfdiyMGs6Ctry9fiktcuQZ/OfbD0hqVol97O8h7mlPmxYOMC3ZdlCNHJiqFYW1eLy964DF9UfYF3b3jXtg302EePYUrJFEwdNhWTzptk6RoYQqwX/nKO7rrP5crdtW6ymtl9ZW35Wlz6+qXo3bk33rn+HdsCiFPmx7yyebovzxBigBDmfbVVW1eLy9+8vGkDnZ13tqm3F81jHz2G4tJiFA8txqTzJtmyhhA7b9uv3Hqfu3XdZD3Ve2XdvnUygFz3Dk7IOMHy/ei0+XFL31vwEl7SdR2GEJ1SUlKQmpoa/4IJqK2rxRVvXoEvqr7Aeze9Z9v38KZ/OB3FpcV45IJH8OCQB21Zw5YDWwCYe39TdG69z926brJGSkpK058q90roEZA+nftg6Y1LbXsOiNPmx5CuQxhCVEtJSWnayCrV1tVi1Ova9/DsfhLq5FWTbX9xp/Hvjgdg3v1N0QUABFx5nweQkhL5RSyJgOYhRFVfaXoOiM1PQnXi/Cj9Kvor1bfkxo7jGU55FrMTfgom9Kzy0zucbsvtExHpxZ+CkZKdHwwhNvHKBlIh/EDPunSWLWsgItKDAURSMT8YQmzgpQ2UrJYHul1GO1vWQUQUDwOIpGp+MIRYzGsbKBlOONBERHo4oV95cX7wiak67dixAx3qOyT1OQ4fPYyxH4/F19VfY86QOcj8PhObv9+saIX6Pb/lecz6chYm/PsEXN3hamzebP0aNh3ahNs/vB3/lv1veOqsp7B7+24AwNfff639+fXXyDiYYfm6/KjH0aNIB3D06FFss2EvJOro0R4A0nH06FFs3rzN7uWQgyXbV6L1Kyu5aX7s+HaH7s/HEKLTps834ciuIwlf/0jjEfz3D/+Nfcf24T9O+g9Uf1mNj/CRwhXqs/THpfjH4X/ginZXoPfB3vjoI+vXsOvoLvyf7/8P8tLyUJhSiLJPy5r+bc/RPQCAjRs34lD6IcvX5kfd6+uRDqCuvt6W/ZCo+vruANJRX1/nqnWT9ZLpK7H6lVXcNj+212zX/TkZQnQqKyvD3oa9CAQCTb8EJvzvIaGPhf/b0ZSj+PC0D1EdrMaw3cOw4+cd2AH9SbGlSLehx5edvsQXXb9A7+96o82BNliBFa0+ZzKfX8+6D7Y5iNLTStG+rj1+vfvXWN24utm/HwoeAk4H1q1bh2/qvol5+4mur+X1on2elvcJgJjXS+b+Uvm59HxuQP5frq2rwwkA6urqsGJF5D2RzJriXc9IPcLXXVd3LYATWq070fXpPdsqPr+Ky+q9XiJnO5H/a6QzEu3zqdoTeq/3fZvvm/rK10e+1v35DrU9hJLTSpBdl92qX+mlZy/FWk+s+RHvbCd6tiKJND+ifb6qlKqon6fVbQv+qsGYampqkJ2djY69O6LND20MX78xvREHRx7E0Q5H0WlpJ2RU2fMthpqzalA7oBZZn2WhfVl7W9ZQ37keB0YdQPqhdHRc1hEpR1s/Jam+Yz2qfluFzos689sxFllXUYHcxkZUpKTg7Nxcu5ejW0XFOjQ25iIlpQK5uWfbvRxysET6ip5+ZTa3zo8jJx7BwS8Oorq6Gu3bx748HwnR6eCBg8B3Bq+UAeB3AE4EsACoKtefDpUaCmAAgBVAbWktalFr/RpOBjASQCVQ/2o9KuorIl/u+EsdVFVVAVEuQmo1hP5sbER5ebmtazFGW3ljY4PL1k2WM9pX9PYrM7l5fhzT/+n50zFmCW2gLgBeAWBXjxwK4EIAKwDo/yV2ap0MYDSA/QBeBVBv0zqIiOJxQr/y0fxgCDGDjzZQXE440EREejihX/lsfjCEqOazDRSTEw40EZEeTuhXPpwfDCEq+XADReWEA01EpIcT+pVP5wefmKpXADFffVFkCIibBNAZCLwaQKAiYEvEE0MExHCBwMoAAh/ZtIY8AfE7AVQBgb8FEDimfx0iICCg/QhhgK+Mao3Gxqa/uumVi8OW7ap1k/Vi9ZVk+pWy9XlsfjQGGuNf6DiGEJ1SU1ORkhq5IiJDoOH6BqAzkPp6KlK+SwFSLV4ggIbzGiAuEEhZlYLU1am2rKExrxGNNzQiUBVA6hupCDQEDK1DpAocwzGkpqYikMoQYomwaZ6aasOmSVB4CHHTusl60fpKsv1Kydo8OD8CqQE0NP3cXWwMIToJISL+UpZmG+hvqQjsC0DA+l+90nh+IxovaETKqhSkfJRiyxpEnkDDDQ1AFZDyegpQD8PrCN3HQgjY8F/wPbf+2iC3rpusEamvqOhXSa/Lo/PDyHlkCNGpsaERjcdaPMSUAeB6AJ0BvAI0lOtLfsoNBXABgBVAY2kjGqH/oTBlTgZwI5q+p9pQn+B9cfxqDQ0Nhn7WnNQ4dsydd7pb100WadlXVPWrZHh5fhj4r/AbqYny6ZOIInLCk7qIiPRwQr/i/GjCR0ISwQ0kOeFAExHp0RnAKDCAAM6YH+AjIcZxA0kMIETkJpeBAQRwxvw4jiHECG4giQGEiNzmEBhAnDA/wjCE6JUGbqAQBhAicqOlYACxe360wBCiVwG4gQAGECJyLzt+iIoBJCaGEL1OAjcQAwgRkX4MIHExhOj1Hvy9gRhAiIj0YwDRhSFErwM23a4TNhADCBGRfgwguukOId9++62Z66BInLCBGECIiPRjADFEdwjp06cPXnnlFTPXosuzzz6L/Px8tGnTBv3798eHH34Y8/IlJSXo378/2rRpg1/84hd47rnnLFppkpywgRhAiIj0YwAxTHcIeeyxxzB+/Hhcc801OHjwoJlriurNN9/EPffcgwcffBBlZWUYMmQILr30UuzZsyfi5Xfu3IlRo0ZhyJAhKCsrwwMPPIC7774bCxcutHjlBjlhAzGAEBHpxwAidTJwWWHAjh07xPDhw0XXrl3FW2+9ZeSqSgwcOFCMHTu22cd69uwpJk6cGPHy999/v+jZs2ezj91xxx3i3HPP1X2b1dXVAoBAV2h/mv02FALFx/+04vYivZ0MgYkQ+AMEMiy+7dzj//9cG///PnvbCwhx/E+712Lsba/Qlr7XAWvhm6PfzO4rGdD65URo/dOu/6dT5keR9vfq6uq4M9bQa8fk5+djxYoVmDVrFq655hr06tULaWnNP8WGDRuMfErd6uvrsX79ekycOLHZxwsKCrB69eqI11mzZg0KCgqafeySSy7B3LlzcfToUaSnp7e6Tl1dHerq6prer6mpUbB6nZyQYPkICBGRfnwERArNj+/0X8XwC9jt3r0bCxcuRIcOHXDVVVe1CiFmOXDgABoaGtC1a9dmH+/atSsqKysjXqeysjLi5Y8dO4YDBw4gNze31XVmzJiBqVOnqlu4Xk7aQHYGEL6kIhG5BQOIFD4/3td/NUMt/4UXXsAf//hHXHTRRdi8eTM6d+5saI0qBAKBZu8LIVp9LN7lI308ZNKkSSgqKmp6v6amBt26dUt0ufo4bQPZ+doKo2y4XSIioxhApJbz4yT9V9UdQkaOHIm1a9di1qxZuPnmm40uMWmdOnVCampqq0c99u/f3+rRjpCcnJyIl09LS0PHjh0jXicYDCIYDKpZtB5O3EB2vrZCBxtum4jICAYQKcn5ofunYxoaGvD555/bEkAAICMjA/3798fy5cubfXz58uUYPHhwxOsMGjSo1eXff/99DBgwIOLzQSzngQ2kRPiBXmLD7RMR6cUAIqmYH7p/TMQB3njjDZGeni7mzp0rvvzyS3HPPfeIdu3aiV27dgkhhJg4caIYPXp00+V37NghMjMzxb333iu+/PJLMXfuXJGeni7+/ve/675N0346xinPYrbrp2BCby2fVc6fjrH8jT8dwzfPv6nqK/wpGPkWa34cn5fKfzrGbtdddx0OHjyIRx55BBUVFejTpw+WLl2K0047DQBQUVHR7HeG5OfnY+nSpbj33nvxzDPPIC8vD08//TSuueYau/4LGq8k2GRF+oqi9XOFiYjsx0dAJIXzIyDE8WdqUkQ1NTXIzs4GusLQjx1F5bENlLBoBzoXwB0AngdQYcO6fGgvgFMAfAvA5KdgK+belZPFku0rDCCSnvlxfF5WV1ejffv2MT8dX8DOSm7ZQGZzyoEmIorHKf3Ko/ODIcQqHt1AhjnlQBMRxeOUfuXh+cEQYgUPbyBDnHKgiYjicUq/8vj8YAgxm8c3kG5OOdBERPE4pV/5YH4whJjJBxtIF6ccaCKieJzSr3wyPxhCzOKTDRSXUw40EVE8TulXPpofDCFm8NEGiskpB5qIKB6n9CufzQ+GENV8toGicsqBJiKKxyn9yofzw1W/MdVOKakpSEmLndkaz29E4wWNSFmVgpTVKbbcuyJPoOHGBqAKSH0jFYHGgOXrEBkCDTc0AJ2B1L+lIvCd/jWIVIEGNCA1NRWBtOivjkwKHTvW9Ne0NPe0hLBlu2rdZL1YfSWZfqWSl+ZHY2ojGtGo67I8uToFAgEEAtGHYsN5DU0bKPXjVMCG+dmY14iGGxoQqApoG+howPJ1iAyBxusbtQP9eipSKlKMreH4ZePd32QOt97nbl03WSRKX0m6Xynitflh5DwyhOjU0NAA0RD5N9yLIQLiAoHAygDwIdCABotXpyVYcYMAqgC8BjTW60uhSteQISCuF0BnIPBqAGKfMHxfhO7jhoYGBBo4WKzW0GD93lXBresma0TqKyr6lZK1eXB+NDbovz5DiF4CaGyMcMcOBTAcwApAlAoI2PBSPCdD+37m8e/hiXob1pEB4EYAnQG8AojyBNdw/CpCCIhGvqyR1SLucRdw67rJIi37iqp+lSyvzg8DV+cTU5PhwycRReSUJ3UREcXjlH7F+QGAISRx3EAapxxoIqJ40uCMfsX50YTfjkkEN5CGAYSI3GQUgJPAAOKE+XEcHwkxihtIwwBCRG7TAQwgTpgfYfhIiE6du3RG3cA61PSvQfv17ZG9MxvoZv066jrVoerSKqR/n47Oqzojpav1ObIxvRFVl1Th6ElH0fndzgimBJXdF/Ud6/EdvkPXrl2RkZah5pNSTKn79gEN2u9Q6JaXZ/dydNu3LxUNDUBqairy8mw4jOQaob5y0pqTcELKCbb07uozq30zP46ceARV31XpuixDiE5Z52WhKqcKv6r6FXq36w2cZ/0aDrY5iFXdVqFjXUcMqx6G9IHplq/haMpRlJxSAhEUuGjvRejYoyPQQ93nPxQ8hPfxPvr27YsOdR3UfWKKKuOtt4Cff0ZGRgbOO8+GjZ2gt97KwM8/w3XrJuuF+srZp52NDjnW95UvOn6BvZ33+mZ+7A/sx4pNK3RdliFEpx0n7sCVJ1yJy3Ivs+X2d9bvxP8c+h+cmnYq7s67G21Ob2P5Go40HsHTh57Gj8d+xB87/BH5p+Qrv409R/fg/QPv4+yzz8ap6acq//zUWnDZMuDnnxEMBjF8+HC7l6PbsmVB/PwzXLdusp6dfWVJ7RJs+nGTr+bH9prtWPE6Q4hS1/ziGkw7f5ott/35wc8xq3QWenXohTlD56BdejvL13D46GHcXno7vsN3eOnCl/Drjr825Xa+/P5LYDlw1lln4d9P+ndTboOaCwaDTX8OGTLE5tXo59Z1k/Xs6ivPffkc3q54G3f3uRtj/32sZbcbzo750fbbtrovyxCi093D7kavM3pZfrtry9fi9rdux69zfo1lNy1DVjDL8jXU1tVi5Gsj8c2P3+CDwg8w8OSBpt3WzxU/AwDy8/PRK9f6+9uXjr/uSnpaGnr1cs99Hnq5mLS0dFetm6xnR195tPRRPL35aUwbPg0PDX3Ikttsya75UZWq7/kgAH86xtHWlq/Fxa9cjD5d+tgeQDbv34zlo5ebGkCIiLzg0dJH8fDKh20PIHbPDz0YQhzKCRuIAYSIyBgGEGMYQhzICRuIAYSIyBgGEM1P9T/pvixDiMM4YQMxgBARGcMAoqmtq8V9y+/TfXk+MVWnhoYG018ufG35Woz820j07twbS65fgsy0zKi3GQgEIIT6V1ysravFqNdH4YuqL7DsxmXon9Pf0pdJD91Wy/vbrP9vsiKtS89ajf5/zPz/pwAIQHvhy8bj93m820t2PWr+P3LlDQZeOtwORv6/Tt3rZon2/9V7tvTcX9H6it7rxzP9w+mYUjIFU4dNxQPnPxCzZyZytvSs0cj8MEttXS0ue/0y7Px+p+7rMIToJIQw9eXC1+1bh0tfvxS9O/fGO9e9g3bp7Sx/efLaulpc/ubl+KLqC7x7w7sYkDvA8jWEDprZ9zdJ4Q+Huus+lyt317rJamb2lcc+egzFpcUoHlqMSedNsnzwA86bH3+66E8YP228rusxhBgQCARM+bxry9c2baAl1y+x7SG08ABi17dg5m2c1/R3s+5vis6t97lb103WU7lXQgEk9AiIHZw4P44cPqL7ugwhOqWkpCAlRf1TaEIbqE/nPlh641L7NtAb2gZ676b3bAsg0z+cjuc3PA/AvPubYnPrfe7WdZM1QvtDZV8JfQvmkQsewYNDHlTyOY1y6vwo/Ur/q/MxhOhkxlBcW74Wl7x2ie1PIhr1+ihsrrL3SaiPlj6KyasmY9yAcZj92WyGEBsEAARceZ8HkJLCR0IoOtUhJNSv7H4SqlPnh5H72I0dxxOc8ixmJ/wUTPizysf0G2PLGoiI9OBPwWhUzQ+GEBt4aQMlywkHmohIDyf0K6/ND4YQi3ltAyXDCQeaiEgPJ/QrL84PhhALeXEDJcoJB5qISA8n9Cuvzg+GEIt4dQMlwgkHmohIDyf0Ky/PD9eEkO+//x6jR49GdnY2srOzMXr0aPzwww8xr3PLLbcgEAg0ezv33HOtWXAYL28go5xwoImI9HBCv/L6/HDNj+jeeOON+Pbbb7Fs2TIAwO23347Ro0fjH//4R8zrjRw5EvPnz296PyMjw9R1tuT1DWSEEw40EZEeTuhXfpgfrgghW7ZswbJly/DJJ5/gnHPOAQC88MILGDRoELZt24YePXpEvW4wGEROTo5VS23GDxtILyccaCIiPZzQr/wyP1zx7Zg1a9YgOzu7KYAAwLnnnovs7GysXr065nVXrVqFLl264IwzzsBtt92G/fv3x7x8XV0dampqmr0lwi8bSA8nHGgiIj2c0K/8ND9cEUIqKyvRpUuXVh/v0qULKisro17v0ksvxWuvvYYVK1bgz3/+M9atW4cLL7wQdXV1Ua8zY8aMpuedZGdno1u3bobX66cNFI8TDjQRkR5O6Fd+mx+2hpDi4uJWTxxt+fbZZ58BiPyiQ0KImC9GdN111+Gyyy5Dnz59cMUVV+Ddd9/FV199hSVLlkS9zqRJk1BdXd30tnfvXkP/J79toFiccKBJh5wc4OSTtT9dxKXLJodyQr/y4/yw9TkhEyZMwPXXXx/zMt27d8fnn3+O7777rtW/VVVVoWvXrrpvLzc3F6eddhq2b98e9TLBYBDBYFD35wznxw0UjRMONOl0POi7jUuXTQ7khH7l1/lhawjp1KkTOnXqFPdygwYNQnV1NdauXYuBA7U75dNPP0V1dTUGDx6s+/YOHjyIvXv3Ijc3N+E1R+PXDRSJEw40EZEeTuhXfp4frnhOSK9evTBy5Ejcdttt+OSTT/DJJ5/gtttuw+WXX97sJ2N69uyJxYsXAwB+/PFH3HfffVizZg127dqFVatW4YorrkCnTp3wm9/8Run6/LyBWnLCgSYi0sMJ/crv88MVIQQAXnvtNfzqV79CQUEBCgoK8Otf/xqvvPJKs8ts27YN1dXVAIDU1FRs2rQJV111Fc444wwUFhbijDPOwJo1a5CVpa7Ift9A4ZxwoImI9HBCv+L8cMnvCQGADh064NVXX415GSFE09/btm2L9957z9Q1cQNJTjjQRER6OKFfcX5oXPNIiNNwA0lOONBERHo4oV9xfkiueSTESbiBJCccaCIiPV7c8CJmfzabAcQh8wPgIyGGcQNJDCBE5CYMIM6ZHyEMIQZwA0kMIETkNuMGjGMAccD8CMcQotOWqi3cQMcxgBCRG43pN8aW22UAiY4hRKei94q4gcAAQkRkBANIbAwhOuWflO/7DcQAQkSkHwNIfAwhOj158ZO+3kAMIERE+jGA6MMQolNmRqblt+mUDcQAQkSkHwOIfgwhDuWUDcQAQkSkHwOIMQwhDuSUDcQAQkSkHwOIcQwhDuOUDcQAQkSkHwOI9PK/XtZ9WYYQB3HKBmIAISLSjwFEerT0UczdMFf35RlCHMJJG4gBhIhIHwYQKTQ/bu13q+7rMIQ4gNM2kJ0BZPP+zbbcLhGRUQwgUvj8uLnvzbqvxxBiMyduIDtfW+HOJXfacttEREYwgEjJzA+GEBt5YQOpEjrQp3c43ZbbJyLSiwFESnZ+MITYxCsbSIXwAz3r0lm2rIGISA8GEEnF/GAIsYGXNlCyWh7odhntbFkHEVE8DCCSqvnBEGIxr22gZDjhQBMR6eGEfuXF+cEQYiEvbqBEOeFAExHp4YR+5dX5wRBiEa9uoEQ44UATEenhhH7l5fnBEGIBL28go5xwoImI9HBCv/L6/GAIMZnXN5ARTjjQRER6OKFf+WF+MISYyA8bSC8nHGgiIj2c0K/8Mj8YQkzilw2khxMONBGRHk7oV36aHwwhJvDTBorHCQeaiEgPJ/Qrv80PhhDF/LaBYnHCgSYi0sMJ/cqP84MhRCE/bqBonHCgiYj0cEK/8uv8YAhRxK8bKBInHGgiIj2c0K/8PD8YQhTw8wZqyQkHmohIDyf0K7/PD4aQJPl9A4VzwoEmItLDCf2K84MhJCncQJITDjQRkR5O6FecHxqGkARxA0lOONBERHo4oV9xfkgMIQngBpKccKCJiPTYvH+z7f2K86M514SQ6dOnY/DgwcjMzMSJJ56o6zpCCBQXFyMvLw9t27bFBRdcgC+++CKpdXADSQwgROQmdy65kwEEzpgfIa4JIfX19bj22msxbtw43deZOXMmnnrqKcyaNQvr1q1DTk4OLr74YtTW1ia0Bm4giQGEiNzm9A6nM4A4YH40I1xm/vz5Ijs7O+7lGhsbRU5Ojnj88cebPnbkyBGRnZ0tnnvuOd23V11dLQCIdze9KwbPHSzaz2gvPv3200SWrsS0kmkCxRDTSqbZtoZPv/1UtJ/RXgyeO1jUHKlR+rnX71svUAyxft96pZ+XiPwr1FdKd5Xacvs1R2p8NT9KtpUIAKK6ujruZdPsDkFm2blzJyorK1FQUND0sWAwiGHDhmH16tW44447Il6vrq4OdXV1Te/X1NQAAMYvHY/vjn2HZy97FmkpadhQscHc/0AEL254EbM/m41xA8Zh1C9H2bKGzfs3484ld+L0Dqfj8RGPY/uh7Uo//5YDW5r9SUSUrFA/2VOzx/K+ebj+MCa8OwHfHPrGN/Nj28Ftui8bEEII01Zigpdeegn33HMPfvjhh5iXW716Nc477zyUl5cjLy+v6eO33347du/ejffeey/i9YqLizF16tTW/zARQJskFk5EROQHRwA8DlRXV6N9+/YxL2rrIyFRB36YdevWYcCAAQnfRiAQaPa+EKLVx8JNmjQJRUVFTe/X1NSgW7duuPyMyzF1ZOy1mmXsO2Oxbt86nJ13Np67/Dlb1vD21rcxtXQqMtMyseh/L0LnEzqbcjtbDmzB7xb9Dq/+9lX06tSr2b+1/IqiT5c+pqwhnvCvKMb0G2PLGsIfkZp16Sy0y2hn+RpYD4n1kJxYjz01e6L2FbNU/ViF3/7f3+KnYz9hytApuLLnlZbcbkt2zI9FGxdh+uPT9V3Y1G8MxVFVVSW2bNkS8+3nn39udh29zwn55ptvBACxYcOGZh+/8sorxc0336x7jaHnhMxZPUf3dVQasWCEQDHEiAUjbLl9IYSYt2GeQDFE1mNZory63NTbivacEL99TzUWM5+ToxfrIbEeklPrYfVzzcqry0XWY1kCxRDzNsyz5DYjsWt+zFk9R/dzQjz/xNQnnnii6WN1dXUJPzHVjhDitwAiRORmwQYrceBJrIeG9ZCi1cPKEOL3ACKER0PI7t27RVlZmZg6dao44YQTRFlZmSgrKxO1tbVNl+nRo4dYtGhR0/uPP/64yM7OFosWLRKbNm0SN9xwg8jNzRU1NfqbhV0hxI8BRIjWzYINVuLAk1gPDeshxaqHVSGEAUTjyRBSWFgoALR6W7lyZdNlAIj58+c3vd/Y2CimTJkicnJyRDAYFEOHDhWbNm0ydLt2hBC7N5AQ9gQQIZo3CzZYiQNPYj00rIcUrx5WhBAGEKnwjULvhRC7WB1CnLCB7AogQjT/eX42WA0HnsR6aFgPSU89zA4hDCDSvA3zBCaCIUQVK0OIYzaQTQFECNks+s7uywYrOPDCsR4a1kPSWw8zQwgDiBSaH8GHggwhqlgVQpy0gewKIEIIsWz7MoFiiHbT27HBcuA1YT00rIdkpB5mhRAGECl8fsz8YCZDiCpWhBCnbSC7Akh5dbnInJ4pUAyxYOMCW9YghPsarFk48CTWQ3JjPcwIIQwgUsv54cknptrF7BDixA1kh/ADbedrx7ixwZqBA09iPSS31kN1CGEAkSLND4YQhcwMIU7dQFYLP9CTV062LYS4tcGqxoEnsR6Sm+uhMoQwgEjR5gdDiEJmhRAnbyArtTzQdr2KrpsbrEoceBLrIbm9Hqr6CgOIFGt+MIQoZEYIcfoGskqkA21HCHF7g1WFA09iPSQv1ENFX2EAkeLND4YQhVSHEDdsICtEO9BWhxAvNFgVOPAk1kPySj2S7SsMIJKe+cEQopDKEOKWDWS2WAfayhDilQabLA48ifWQvFSPZPoKA4ikd34whCikKoS4aQOZKd6BtiqEeKnBJoMDT2I9JK/VI9G+wgAiGZkfDCEKqQghbttAZtFzoK0IIV5rsIniwJNYD8mL9UikrzCASEbnB0OIQsmGEDduIDPoPdBmhxAvNthEcOBJrIfk1XoY7SsMIFIi84MhRKFkQohbN5BqRg60mSHEqw3WKA48ifWQvFwPI32FAURKdH4whCiUaAhx8wZSyeiBNiuEeLnBGsGBJ7EektfrobevMIBIycwPhhCFEgkhbt9AqiRyoM0IIV5vsHpx4Emsh+SHeujpKwwgUrLzgyFEIaMhxAsbSIVED7TqEOKHBqsHB57Eekh+qUe8vsIAIqmYHwwhChkJIV7ZQMlK5kCrDCF+abDxcOBJrIfkp3rE6isMIJKq+cEQopDeEOKlDZSMZA+0qhDipwYbCweexHpIfqtHtL7CACKpnB8MIQrpCSFe20CJUnGgVYQQvzXYaDjwJNZD8mM9IvUVBhBJ9fxgCFEoXgjx4gZKhKoDnWwI8WODjYQDT2I9JL/Wo2VfYQCRzJgfDCEKxQohXt1ARqk80MmEEL822JY48CTWQ/JzPcL7CgOIZNb8YAhRKFoI8fIGMkL1gU40hPi5wYbjwJNYD8nv9Qj1lWXblzGAHGfm/GAIUShSCPH6BtLLjK8oEgkhfm+wIRx4EushsR6yr2ROz2QAEebPD4YQhVqGED9sID3MekjTaAhhg9Vw4Emsh8R6aEp3lQoUgwFEWDM/GEIUCg8hftlA8Zj5PVUjIYQNVsOBJ7EeEuuhqTlSI/rO7itQDDF55WRb1iCEfwKIEAwhSoVCSM8/9fTNBorF7Cd16Q0hbLAaDjyJ9ZBYD02oHu2mt1P+chBG+CmACMEQolQohGCifzZQNFY8q1xPCGGD1XDgSayHxHpowuuxYOMC20KI3wKIEAwhSoU/EmIXvwQQIeKHEDZYDQeexHpIrIemZT3MenXuePwYQIRgCFEqkVfRVclPAUSI2CGEDVbDgSexHhLroYlUDztCiF8DiBAMIUrZGUL8FkCEiN4s2GA1HHgS6yGxHppo9bA6hPg5gAjBEKKUXSHEjwFEiMjNgg1Ww4EnsR4S66GJVQ8rQ4jfA4gQDCFK2RFC7N5AQtj32gotmwUbrIYDT2I9JNZDE68eVoUQBhDNzA9mMoSoYnUIccIGsvO1FcKbBRushgNPYj0k1kOjpx5WhBAGEE15dbkIPhRkCFHFyhDilA1k52srhJrFuHfGscEKDrxwrIfEemj01sPsEMIAommaHxPBEKKKVSHEURvIxl9tHGoWbLAceOFYD4n10Biph5khhAFEEz4/Ct8o9F4IefTRR8WgQYNE27ZtRXZ2tq7rFBZqd0T42znnnGPodq0IIU7bQHa+tsLAOQObHgmxi9sarJk48DSsh+TGepgVQhhANC3nh5EnpqbAJerr63Httddi3Lhxhq43cuRIVFRUNL0tXbrUpBUmZn7ZfPzh7T8gKyMLW8dvRV77PMvXsK9mH3o+0xO19bWYd+U8/P6s31u+BgC46OWLsHbfWgDAmH5jbFnD2vK1uPiVi9GnSx8su2kZsoJZlq+htq4WI18bic37N2P56OUYePJAy9cAAI+WPoqHVz6MacOn4aGhD9myBtZDYj00TqnHRS9fhA92foAR+SPwz5v/acsaPDE/LAhJSs2fP9/QIyFXXXVVUrdn5iMhTkywdgl9RRF6JMSOX6/sxq/wzMKvuDWsh+Tmeqh+JISPgGiizQ9PPhKSqFWrVqFLly4444wzcNttt2H//v12LwmARxKsIuFfUcy+fLYta+BXeBK/4tawHhLrIfEREI2y+WFmSjKDkUdC3njjDfHOO++ITZs2ibffflv07dtX9O7dWxw5ciTqdY4cOSKqq6ub3vbu3av8kRAnJ1irtfyKwo5fr+zmr/BU41fcGtZD8kI9VPUVPgKiiTc/XPPLyqZMmdLqiaMt39atW9fsOkZCSEv79u0T6enpYuHChYbXpCqEuGEDWSXSgbY6hHihwarCgadhPSSv1ENFX2EA0eiZH0ZCSFpyD8gkZ8KECbj++utjXqZ79+7Kbi83NxennXYatm/fHvUykyZNQlFRUdP7NTU16Natm5Lb99RDaElywkOafIhZ4kP+GtZDYj0kJ/Qrr84PW0NIp06d0KlTJ8tu7+DBg9i7dy9yc3OjXiYYDCIYDCq/ba9uoEQ44UCzwUoceBrWQ2I9JCf0K0/PD9UP1Zhl9+7doqysTEydOlWccMIJoqysTJSVlYna2tqmy/To0UMsWrRICCFEbW2t+OMf/yhWr14tdu7cKVauXCkGDRokTj75ZFFTo/9hRRU/HeOWh9CsEO8hTSu+HeOVh5hV4EP+GtZD8mI9Eu0r/BaMxuj8cM1zQoyI9IvHAIiVK1c2XQaAmD9/vhBCiJ9++kkUFBSIzp07i/T0dHHqqaeKwsJCsWfPHkO3m2wIceMGMoueA212CPFig00UB56G9ZC8Wo9E+goDiCaR+eHJEGKXZEKIWzeQGfQeaDNDiFcbbCI48DSsh+TlehjtKwwgmkTnB0OIQomGEDdvINWMHGizQoiXG6xRHHga1kPyej2M9BUGEE0y84MhRKFEQojbN5BKRg+0GSHE6w3WCA48Desh+aEeevsKA4gm2fnBEKKQ0RDihQ2kSiIHWnUI8UOD1YsDT8N6SH6ph56+wgCiUTE/GEIUMhJCvLKBVEj0QKsMIX5psHpw4GlYD8lP9YjXVxhANKrmB0OIQnpDiJc2ULKSOdCqQoifGmw8HHga1kPyWz1i9RUGEI3K+cEQopCeEOK1DZSMZA+0ihDitwYbCweehvWQ/FiPaH2FAUSjen4whCgUL4R4cQMlSsWBTjaE+LHBRsOBp2E9JL/WI1JfYQDRmDE/GEIUihVCvLqBEqHqQCcTQvzaYCPhwNOwHpKf69GyrzCAaMyaHwwhCkULIV7eQEapPNCJhhA/N9iWOPA0rIfk93qE9xUGEI2Z84MhRKFIIcTrG8gI1Qc6kRDi9wYbjgNPw3pIrIfsKwPnDGQAEebPD4YQhVqGED9sIL3M+IrCaAhhg5U48DSsh8R6aEJ9hQHEmvnBEKJQeAjxywbSw6yHNI2EEDZYiQNPw3pIrIc07p1xTY+E2MVP84MhRKFQCCl8o9A3GygeM7+nqjeEsMFKHHga1kNiPaRQPcx4TSq9/BRAhGAIUSoUQjDRPxsoFrOf1KUnhLDBShx4GtZDYj2kUD1Cj4TYEUL8FkCEYAhRKhRCgg8FfbOBorHiWeXxQggbrMSBp2E9JNZDCq+HWa/OHY8fA4gQDCFKhULIzA9m2nL7fgogQsQOIWywEgeehvWQWA+pZT3sCCF+DSBCMIQoZfRVdFXyWwARInqzYIOVOPA0rIfEekiR6mF1CPFzABGCIUQpu0KIHwOIEJGbBRusxIGnYT0k1kOKVg8rQ4jfA4gQDCFK2RFC7N5AIXb8ZsGWzYINVuLA07AeEushxaqHVSGEAUTDEKKQ1SHECRtICPteWyG8WbDBShx4GtZDYj2kePWwIoQwgEg9/9STIUQVK0OIUzaQna+tEGoWCzYuYIM9jgNPw3pIrIekpx5mhxAGEGnEghECE8EQoopVIcRRG8jGX20cahbtprdjgxUceCGsh8R6SHrrYWYIYQCRQvODj4QoZEUIcdoGsvO1FSavnCxQDNF3dl82WA48IQTrEY71kIzUw6wQwgAihc8PPidEIbNDiBM3kF1CBxrFEKW7Sm1ZgxsbrFk48CTWQ+PWepgRQhhApJbzgyFEITNDiFM3kB1CBzpzeqZtv17ZrQ3WDBx4EuuhcXM9VIcQBhAp0vxgCFHIrBDi5A1ktfADvWz7MltCiJsbrGoceBLroXF7PVSGEAYQKdr8YAhRyIwQ4vQNZKWWB9qOX6/s9garEgeexHpovFAPVX2FAUSKNT8YQhRSHULcsIGsEulAWx1CvNBgVeHAk1gPjVfqoaKvMIBI8eYHQ4hCKkOIWzaQFaIdaCtDiFcarAoceBLrofFSPZLtKwwgkp75wRCikKoQ4qYNZLZYB9qqEOKlBpssDjyJ9dB4rR7J9BUGEEnv/GAIUUhFCHHbBjJTvANtRQjxWoNNBgeexHpovFiPRPsKA4hkZH4whCiUbAhx4wYyi54DbXYI8WKDTRQHnsR6aLxaj0T6CgOIZHR+MIQolEwIcesGMoPeA21mCPFqg00EB57Eemi8XA+jfYUBREpkfjCEKJRoCHHzBlLNyIE2K4R4ucEaxYEnsR4ar9fDSF9hAJESnR8MIQolEkLcvoFUMnqgzQghXm+wRnDgSayHxg/10NtXGECkZOYHQ4hCRkOIFzaQKokcaNUhxA8NVi8OPIn10PilHnr6CgOIlOz88FwI2blzp/jDH/4gunfvLtq0aSN+8YtfiMmTJ4u6urqY12tsbBRTpkwRubm5ok2bNmLYsGFi8+bNhm7bSAjxygZSIdEDrTKE+KXB6sGBJ7EeGj/VI15fYQCRVMwPz4WQd999V9xyyy3ivffeE99884146623RJcuXcQf//jHmNd7/PHHRVZWlli4cKHYtGmTuO6660Rubq6oqdF/6PWGEC9toGQlc6BVhRA/Ndh4OPAk1kPjt3rE6isMIJKq+eG5EBLJzJkzRX5+ftR/b2xsFDk5OeLxxx9v+tiRI0dEdna2eO6553Tfjp4Q4rUNlIxkD7SKEOK3BhsLB57Eemj8WI9ofYUBRFI5P4yEkDS4VHV1NTp06BD133fu3InKykoUFBQ0fSwYDGLYsGFYvXo17rjjjojXq6urQ11dXbPbAYAPv/ow4uV/+PkHPLzyYdQdq0PhmYU4duQYXljzQiL/paQ89fFT2HpwK3p27InrfnmdLWv4eM/HWLBxAYJpQTw85GEs+WKJ4c+x+4fdwBFg/a71+LH2R8PX/6n+J9y3/D7s/H4nnrrkKRw5fASlX5Ua/jzJevlfL2Puhrm4td+tGJoz1JY1bKnagqL3ipB/Uj4eHvgwynaXWb4G1kNiPSSr67Ht4DbgCLBo4yKs37UegJp+lSyvzo/QvBRCxL9w0pHHBl9//bVo3769eOGFF6Je5uOPPxYARHl583R72223iYKCgqjXmzJligDAN77xjW984xvfknjbu3dv3Hlu6yMhxcXFmDp1aszLrFu3DgMGDGh6f9++fRg5ciSuvfZajBkzJu5tBAKBZu8LIVp9LNykSZNQVFTU9H5jYyMOHTqEjh07xrye09TU1KBbt27Yu3cv2rdvb/dyfIH3ubV4f1uL97f13HqfCyFQW1uLvLy8uJe1NYRMmDAB119/fczLdO/evenv+/btw/DhwzFo0CDMmTMn5vVycnIAAJWVlcjNzW36+P79+9G1a9eo1wsGgwgGg80+duKJJ8a8LSdr3769qzavF/A+txbvb2vx/raeG+/z7OxsXZezNYR06tQJnTp10nXZ8vJyDB8+HP3798f8+fORkpIS8/L5+fnIycnB8uXLcdZZZwEA6uvrUVJSgieeeCLptRMREVFyYk9yh9i3bx8uuOACdOvWDU8++SSqqqpQWVmJysrKZpfr2bMnFi9eDED7Nsw999yDxx57DIsXL8bmzZtxyy23IDMzEzfeeKMd/w0iIiIK44qfjnn//ffx9ddf4+uvv8Ypp5zS7N9E2LNvt23b1vTTLABw//334+eff8add96J77//Hueccw7ef/99ZGVlWbZ2uwSDQUyZMqXVt5bIPLzPrcX721q8v63nh/s8IISen6EhIiIiUssV344hIiIi72EIISIiIlswhBAREZEtGEKIiIjIFgwhHrdr1y7ceuutyM/PR9u2bXH66adjypQpqK+vt3tpnjZ9+nQMHjwYmZmZrv5ld0717LPPIj8/H23atEH//v3x4YeRX9uJ1CgtLcUVV1yBvLw8BAIB/M///I/dS/KsGTNm4Oyzz0ZWVha6dOmCq6++Gtu2bbN7WaZhCPG4rVu3orGxEc8//zy++OIL/OUvf8Fzzz2HBx54wO6leVp9fT2uvfZajBs3zu6leM6bb76Je+65Bw8++CDKysowZMgQXHrppdizZ4/dS/Osw4cPo2/fvpg1a5bdS/G8kpISjB8/Hp988gmWL1+OY8eOoaCgAIcPH7Z7aabgj+j60J/+9CfMnj0bO3bssHspnvfSSy/hnnvuwQ8//GD3UjzjnHPOQb9+/TB79uymj/Xq1QtXX301ZsyYYePK/CEQCGDx4sW4+uqr7V6KL1RVVaFLly4oKSnB0KFD7V6OcnwkxIeqq6vRoUMHu5dBZFh9fT3Wr1+PgoKCZh8vKCjA6tWrbVoVkXlCv4DTqz2bIcRnvvnmG/z3f/83xo4da/dSiAw7cOAAGhoaWr0IZdeuXVu9jAOR2wkhUFRUhPPPPx99+vSxezmmYAhxqeLiYgQCgZhvn332WbPr7Nu3DyNHjsS1116LMWPG2LRy90rkPidzBAKBZu8LIVp9jMjtJkyYgM8//xyvv/663UsxjSteO4ZamzBhAq6//vqYl+nevXvT3/ft24fhw4dj0KBBmDNnjsmr8yaj9zmp16lTJ6SmprZ61GP//v2tHh0hcrO77roLb7/9NkpLS1u9ZpqXMIS4VKdOndCpUyddly0vL8fw4cPRv39/zJ8/HykpfAAsEUbuczJHRkYG+vfvj+XLl+M3v/lN08eXL1+Oq666ysaVEakhhMBdd92FxYsXY9WqVcjPz7d7SaZiCPG4ffv24YILLsCpp56KJ598ElVVVU3/lpOTY+PKvG3Pnj04dOgQ9uzZg4aGBmzcuBEA8G//9m844YQT7F2cyxUVFWH06NEYMGBA0yN7e/bs4fOcTPTjjz/i66+/bnp/586d2LhxIzp06IBTTz3VxpV5z/jx4/G3v/0Nb731FrKyspoe9cvOzkbbtm1tXp0JBHna/PnzBYCIb2SewsLCiPf5ypUr7V6aJzzzzDPitNNOExkZGaJfv36ipKTE7iV52sqVKyPu58LCQruX5jnR+vX8+fPtXpop+HtCiIiIyBZ8cgARERHZgiGEiIiIbMEQQkRERLZgCCEiIiJbMIQQERGRLRhCiIiIyBYMIURERGQLhhAiIiKyBUMIEblGQ0MDBg8ejGuuuabZx6urq9GtWzc89NBDNq2MiBLB35hKRK6yfft2nHnmmZgzZw5uuukmAMDNN9+Mf/3rX1i3bh0yMjJsXiER6cUQQkSu8/TTT6O4uBibN2/GunXrcO2112Lt2rU488wz7V4aERnAEEJEriOEwIUXXojU1FRs2rQJd911F78VQ+RCDCFE5Epbt25Fr1698Ktf/QobNmxAWlqa3UsiIoP4xFQicqV58+YhMzMTO3fuxLfffmv3cogoAXwkhIhcZ82aNRg6dCjeffddzJw5Ew0NDfjnP/+JQCBg99KIyAA+EkJErvLzzz+jsLAQd9xxBy666CK8+OKLWLduHZ5//nm7l0ZEBjGEEJGrTJw4EY2NjXjiiScAAKeeeir+/Oc/4z//8z+xa9cuexdHRIbw2zFE5BolJSUYMWIEVq1ahfPPP7/Zv11yySU4duwYvy1D5CIMIURERGQLfjuGiIiIbMEQQkRERLZgCCEiIiJbMIQQERGRLRhCiIiIyBYMIURERGQLhhAiIiKyBUMIERER2YIhhIiIiGzBEEJERES2YAghIiIiWzCEEBERkS3+P/u+HkSlL4XdAAAAAElFTkSuQmCC", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "subtraction_geom = [mp.Block(size=mp.Vector3(mp.inf, 0.5, t_Si), material=Si)]\n", - "subtraction_cell = mp.Vector3(5, 4, cell_z if sim_is_3D else 0)\n", - "subtraction_sources = [\n", - " mp.EigenModeSource(\n", - " src = mp.GaussianSource(fcen, fwidth=df),\n", - " volume=mp.Volume(center=mp.Vector3(-0.5,0,0), size=port_size),\n", - " eig_band=1,\n", - " eig_parity = mode_parity,\n", - " eig_match_freq = True,\n", - " )\n", - "]\n", - "\n", - "subtraction_sim = mp.Simulation(\n", - " resolution=res,\n", - " cell_size=subtraction_cell,\n", - " boundary_layers=[mp.PML(dpml)],\n", - " sources=subtraction_sources,\n", - " geometry=subtraction_geom,\n", - " default_material=SiO2\n", - ")\n", - "\n", - "subtraction_monitor_region = mp.FluxRegion(center=mp.Vector3(0), size=port_size)\n", - "subtraction_monitor = subtraction_sim.add_flux(fcen, 0, 1, subtraction_monitor_region)\n", - "\n", - "plot_plane = mp.Volume(center=mp.Vector3(z=0), size=mp.Vector3(cell.x, cell.y, 0))\n", - "subtraction_sim.plot2D(output_plane=plot_plane if sim_is_3D else None)\n", - "subtraction_sim.run(until_after_sources=mp.stop_when_dft_decayed)\n", - "\n", - "subtraction_data = subtraction_sim.get_flux_data(subtraction_monitor)\n", - "norm_mode_coeff = subtraction_sim.get_eigenmode_coefficients(subtraction_monitor, [1], mode_parity).alpha[0,0,0]" - ] - }, - { - "cell_type": "code", - "execution_count": 6, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ - "# Set up the first source for the simulation. I'll start with port1 (the lower left)\n", + "# Set up the first source for the simulation in port 1 (lower left).\n", "sources = [\n", " mp.EigenModeSource(\n", " src = mp.GaussianSource(fcen, fwidth=df),\n", @@ -301,12 +246,13 @@ "# Create Simulation\n", "sim = mp.Simulation(\n", " resolution=res, # The resolution, defined further up\n", - " cell_size=cell, # The cell size, taken from the gds\n", - " boundary_layers=[mp.PML(dpml)], # the perfectly matched layers, with a diameter as defined above\n", - " sources = sources, # The source(s) we just defined\n", - " geometry = geometry, # The geometry, from above\n", + " cell_size=cell, # The simulation size, taken from the gds\n", + " boundary_layers=[mp.PML(dpml)], # the boundary layers to absorb fields that leave the simulation\n", + " sources = sources, # The sources\n", + " geometry = geometry, # The geometry\n", " default_material=SiO2\n", - ")\n" + ")\n", + "sim.init_sim()" ] }, { @@ -314,20 +260,110 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "After creating the simulation object, we can add mode monitors at each of the ports. These will read how much light comes in/out of each port. \n", + "After creating the simulation, we can add mode monitors at each of the ports. These will collect the fields on each port and accumulate the fourier transform of the time-domain fields, giving us frequency-domain field information at the end.\n", "\n", - "Now that our simulation is completely set up for the first source, I can plot it to make sure everything looks right before running the simulation. We do that using the plot2D() method of the simulation object." + "After our simulation is completely set up for the first source, we can plot it to make sure everything looks right before running the simulation. We do that using the plot2D() method." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "metadata": { "tags": [ "hide-output" ] }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " prism, center = (-7.75,-0.625,5e+19)\n", + " height 1e+20, axis (0,0,1), sidewall angle: 0 radians, 4 vertices:\n", + " (-2.75,-1.125,0)\n", + " (-12.75,-0.875,0)\n", + " (-12.75,-0.375,0)\n", + " (-2.75,-0.125,0)\n", + " dielectric constant epsilon diagonal = (12.0946,12.0946,12.0946)\n", + " prism, center = (-7.75,0.625,5e+19)\n", + " height 1e+20, axis (0,0,1), sidewall angle: 0 radians, 4 vertices:\n", + " (-2.75,0.125,0)\n", + " (-12.75,0.375,0)\n", + " (-12.75,0.875,0)\n", + " (-2.75,1.125,0)\n", + " dielectric constant epsilon diagonal = (12.0946,12.0946,12.0946)\n", + " prism, center = (7.75,0.625,5e+19)\n", + " height 1e+20, axis (0,0,1), sidewall angle: 0 radians, 4 vertices:\n", + " (2.75,0.125,0)\n", + " (2.75,1.125,0)\n", + " (12.75,0.875,0)\n", + " (12.75,0.375,0)\n", + " dielectric constant epsilon diagonal = (12.0946,12.0946,12.0946)\n", + " prism, center = (7.75,-0.625,5e+19)\n", + " height 1e+20, axis (0,0,1), sidewall angle: 0 radians, 4 vertices:\n", + " (2.75,-1.125,0)\n", + " (2.75,-0.125,0)\n", + " (12.75,-0.375,0)\n", + " (12.75,-0.875,0)\n", + " dielectric constant epsilon diagonal = (12.0946,12.0946,12.0946)\n", + " prism, center = (0,0,5e+19)\n", + " height 1e+20, axis (0,0,1), sidewall angle: 0 radians, 4 vertices:\n", + " (-2.75,-1.25,0)\n", + " (-2.75,1.25,0)\n", + " (2.75,1.25,0)\n", + " (2.75,-1.25,0)\n", + " dielectric constant epsilon diagonal = (12.0946,12.0946,12.0946)\n", + " prism, center = (15.25,-0.625,5e+19)\n", + " height 1e+20, axis (0,0,1), sidewall angle: 0 radians, 4 vertices:\n", + " (12.75,-0.875,0)\n", + " (12.75,-0.375,0)\n", + " (17.75,-0.375,0)\n", + " (17.75,-0.875,0)\n", + " dielectric constant epsilon diagonal = (12.0946,12.0946,12.0946)\n", + " prism, center = (15.25,0.625,5e+19)\n", + " height 1e+20, axis (0,0,1), sidewall angle: 0 radians, 4 vertices:\n", + " (12.75,0.375,0)\n", + " (12.75,0.875,0)\n", + " (17.75,0.875,0)\n", + " (17.75,0.375,0)\n", + " dielectric constant epsilon diagonal = (12.0946,12.0946,12.0946)\n", + " prism, center = (-15.25,0.625,5e+19)\n", + " height 1e+20, axis (0,0,1), sidewall angle: 0 radians, 4 vertices:\n", + " (-12.75,0.875,0)\n", + " (-12.75,0.375,0)\n", + " (-17.75,0.375,0)\n", + " (-17.75,0.875,0)\n", + " dielectric constant epsilon diagonal = (12.0946,12.0946,12.0946)\n", + " prism, center = (-15.25,-0.625,5e+19)\n", + " height 1e+20, axis (0,0,1), sidewall angle: 0 radians, 4 vertices:\n", + " (-12.75,-0.375,0)\n", + " (-12.75,-0.875,0)\n", + " (-17.75,-0.875,0)\n", + " (-17.75,-0.375,0)\n", + " dielectric constant epsilon diagonal = (12.0946,12.0946,12.0946)\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAACcCAYAAACHmVqXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkoElEQVR4nO3de3gU5d0+8Ht2QzbnQBIgRAhGBaXEEkOqEgUU5VSKvEoVrUKsgBWNFcG+FyppguHgkaogYITw01KsVK21VcC0IoEflEBCUGhEtIGQQIycNiSQ0+7z/rFuJgmbZA+zM7O79+e6cimbOTyZ77Pz3DuzMyMJIQSIiIiIfJxB6wYQERERKYGhhoiIiPwCQw0RERH5BYYaIiIi8gsMNUREROQXGGqIiIjILzDUEBERkV9gqCEiIiK/EKR1A9RktVpx4sQJREZGQpIkrZtDREREThBC4Pz580hISIDB0PnxmIAKNSdOnMCAAQO0bgYRERG54fjx4+jfv3+nvw+oUBMZGQkAmP7/pmPk4JEuzXvs3DEsKVyCvDvycHXs1U7N886Bd7CuZB1mps7EjGEzXG6vEsp+KMO8rfOQ1CsJL499GWHBYaq34ULTBTxV8BTKz5Zj+fjlGNJ7iOptAFgPO9ZDxnrIWA8b1kOmZD0Onz6Mhz9+GM+OehYDew50er5zF88ha1sWGusbgT/I43inRAAxm80CgMjblefyvMUnigVyIIpPFDs1fe72XIEciNztuS6vSyl7KveIqGVRIn1duqhtqNWkDbUNtSJ9XbqIWhYl9lTu0aQNQrAedqyHjPWQsR42rIdM6Xq4OoYKIUSVuUpELo0UyIHI+HOGACDMZnOX8zDUOMmVgvhjh3QHdxAy1kPGetiwHjLWQ+av9XA11LQNNPkl+SJvVx5DTUdqhBp/7ZCu4g5CxnrIWA8b1kPGesj8uR6uhJqOgUYIwVDjiLdDjT93SFdwByFjPWSshw3rIWM9ZP5eD2dDjaNAIwRDjUPeDDX+3iGdxR2EjPWQsR42rIeM9ZAFQj2cCTWdBRohGGoc8laoCYQO6QzuIGSsh4z1sGE9ZKyHLFDq0V2o6SrQCMFQ45A3Qk2gdMjucAchYz1krIcN6yFjPWSBVI+uQk13gUYIhhqHlA41gdQhu8IdhIz1kLEeNqyHjPWQBVo9Ogs1zgQaIRhqHFIy1ARah+wMdxAy1kPGetiwHjLWQxaI9XAUapwNNEL4YahZunSpSEtLExEREaJ3795iypQp4uuvv3ZpGUqFmkDskI5wByFjPWSshw3rIWM9ZIFaj46hxpVAI4Qfhprx48eL9evXi4MHD4rS0lIxadIkkZiYKOrq6pxehhKhZs4/5gRkh+yIOwgZ6yFjPWxYDxnrIQvkerQNNa4GGiH8MNR0VFNTIwCI7du3Oz2PEqEmUDtkW9xByFgPGethw3rIWA9ZoNfDPoZuObLF5UAjhPOhxmcfaGk2mwEAMTExnU7T2NiIxsbG1n/X1tYCsD2csuRkiUvrm/Gh7YFivxzyS/x80M9dnl8JB2sO4tFPHsWVMVfi+duex5EzR1RvQ31TPTI3Z+K7M99h1aRVCDIEabIt1pasxep9qzEnbQ7rwXoAYD3aYj1sWA+Z1vUoO1UGAPif9/4HDS0NyB6VjWHxw5zeFsfOHXNqOkkIIdxupUaEEJgyZQrOnj2LHTt2dDpdTk4OFi1adOkvFgAI8V77iIiISEENAJ63HdCIiorqdDKfDDWPPfYYPvnkE+zcuRP9+/fvdDpHR2oGDBiAvF15GH75cKfWZU/YvxzyS7xf9j423LUBQ+K0eRQ9ERGRLyo7VYYHPnwAi25ZhF8M/oXL8xcfLcbD6Q93G2p87vTT448/jo8//hiFhYVdBhoAMJlMMJlMl7x+dezVSO2X2u26Fhcuxup9q5F7ay5+PujneL/sfQyJG+LUvERERNTeLwb/wq0xtO58nVPTGVxeskaEEMjMzMSHH36Izz//HElJSV5d3+LCxcjaloXcW3OxcNRCr66LiIiIPOczR2oee+wxbNy4EX/7298QGRmJ6upqAEB0dDRCQ0MVXRcDDRERke/xmSM1q1evhtlsxi233IJ+/fq1/rz33nuKroeBhoiIyDf5zJEaNb7PzEBDRETku3zmSI23MdAQERH5NoYaMNAQERH5A585/eQtDDREXfPBW1n5BEmStG4Ckd8J6FDDQEPUPQ6+ROQrAjbUMNAQOae6uhoWi4XhRiFCCBiNRsTHx2vdFCK/E5Ch5p0D72Ddf9Yx0JBX6OF0jSRJsFqtEEI4DCOdvd7x90IIXHHFFbh48aI3mxtwIiMjce7cOQDta9FdXdpOY/+vwWDQTZ8j0lpAhpp1JeuQO9F3Ak1aGlBdDcTHA/v2ad0a/RNCwGq1qrpO+yBjNBp1s3M3GJS5DiA0NBQXL15s/RvJfQaDAVarFSEhIYrVB9BPoLAf0VO7nxgMBt1sAz0LhLEkIEPNzNSZPhNoAFsnrKrSuhXusVgsqu3g7EcnevToAaPRqMo6HamtrVV0wHKVxWJBdHQ0HnroIezatQuRkZGwWCwuLaPt0YC2RxTIM/awffr0aaSm2p5/48zRmY6MRiPOnz+P9PR05Ofnw2w2a9rnrVYroqKiNG1Dc3OzqketJEnS9O91hy+PJc4KyFAzY9gMrZugGfspCbWo/aY3Go2oqalBaWkpQkJCVPtbLRYLYmNjsXjxYrz//vsIDg5GS0uLKut2xGAwaLp+6prVasX+/fs9Xs7hw4fxxz/+UfUjk20FBQWhqakJ99xzD5555hmcPn1atfe9JEm4ePEirrvuOvTp00eVdbbl6ocFT9hPNVLXAjLUBDK13xSvvvoqzp07h6CgIK8HDKvViujoaGzcuBF79+716rq609TUpOn67YOcEofkeYTGO5Sqjdbh1d7XN23ahE2bNmnShuuvvx733XcfzGaz1/dxkiShpaUFPXv2xNy5c726LnIdQ02AsFqtMBgMyMzMRGlpKcLDw73+6U6SJBQUFHh1HZ3p0aMHWlpaVD/PruUnZkcYSPTLH2uj9ocmIQSCgoJQVFSEoqIiVdcNAJs3b/b6e95gMKC+vh6pqal4/fXXYbFYfO60l5ok4Y/vrE7U1tYiOjoa4VeEo8eZHi7N29K7BXX31yHiTxEI+kHdLGg2H4QQl0GSqhAdnez2ciRJwtmzZxVsmXPrVPsNqPUnVyJSlxb7GTW/L2jXq1cvj9ap1FjiDk/H0Pqb69H8j2aYzWZERUV1Ol1AHqmpr68Hzrk4U6jtP3V1da7P6zFbJ277pU1PqPlpymq1MmQQkVdpdRpO7X2p5x9KlR1LXOLJGDoKQDKAf3Q/aUCGmkCnt1MkRES+iPtSFYwCMAbA/3ducoYaIiIi0h97oPkcwGHnZuH1YURERKQvbQNNofOzBeaRGsn1c6FCEhCw3SRLMqh9RY38/56ew+XhUiIi3+T5/l+5ZbnKlTFUjBQQtwpI2yRIOyXAAFgl58auwAw1wo3B/ccvnAshIKzaXTCmRCjp0aMHLBaLKp2aXxImIn8VFOT9IdRqtcJoNKK5uVnRD6Wqf8B1dgwdBeBWAJ8DotAWhNrO352ADDXXXXcdejb2dGme8xHnsQ/7kJaWhsi6SO80rBO7dpnQ2AiYTCakp9/q1jKEEDCZTNi6dSuam5sBqNep1f5EwKNRRIFHzf2M/Y7dan1os+/Txo8fj8bGRrfvv6XEWOIuZ8bQo4lHUZ5UjqTyJFxuvNwWbn50znQO+7d0fxfugLxPzfbD2zFq8CiX5i05WYLhecNR/HAxUvuleqmFjvXvb3tex2WXAZWVni0rPz8fR44cQWhoqCo3jcrOzvbqOvRILw/WC6C3dkBjf9POokWLWh/i6S0GgwEXLlzANddcgwcffNCjZSk5lriquzF0ceFiZG3LQu6tjh82XfhNIUZfPVq5+9RUVlaif//+zk7uNatWrcJLL72EkydPYujQoXj11VcxcuRIl5ZhsVhcfmaHfXp35vWcAYAEQMBicT+ICCHw0EMPKdYqZ4wZMwYNDQ2qfIqyP39p+fLl+NOf/oSgoCBVT3/ZHwCpp517UFBQ692k3cVTiN7hyakLe01bWlp01d/UfkK3/T0+ffp0zJ07V7XnTlksFoSHhyM9Pd3r62rL87ukKzOWuKOrMXTJjiXI3p6NRaMX4embnnY4xjo77jr9rkpOTsaKFSswffp0Z2dR3HvvvYe5c+di1apVuOmmm/Dmm29i4sSJ+M9//oPExESnl2M0Gl3u+Pbp3ZlXOZ7fNVPtu2DefPPNqq3LbsOGDVixYoWqdbJYLOjVqxdyc3Oxbt06xMTEaHYazH7r+JKSktZAwlNy+uNpWLTXNDU1VZNHgtgZDAacOXMGs2bNwsKFC3H27FlN3ntaUPtDk3Lf4VH/DsydjaGLCxcje3t2p0doOs7fHae30NKlS/HYY4/ho48+Ql5eHmJjY52dVTHLly/HzJkzMWvWLAC2hyVu3boVq1evxrJly1Rvjy9SuyOr/VRwwPY3arWTy8rKQlZWlibr7mjlypUoKSlBRESEW6HG/ql77dq1aGxs9EILA1dISAhmzZrl9nvDYDCgrq4OqampyMzMVLh1ntHqvaf2EXRlQ0Zg6u6Ukzucrsijjz6KiRMnYubMmRg6dCjy8vJwxx13KNIIZzQ1NaG4uBgLFixo9/q4ceOwa9cuh/M0Nja22xnX1tZ6tY3eEh/f/r++RO0vCdtpdUjefvpJ6+85WK1WxQa7jRs3MtQoLCIiAitWrFBkWS0tLZq9z+zsfV6rfq/Fs598kZ7GEm8EGsDFq5+SkpLw+eefY+XKlZg6dSqGDBlySVItKSlRrHFtnTp1ChaLBX379m33et++fVFdXe1wnmXLlmHRokVeaY+a9u3TugW+R8udqx4YDAY0NTXBarW6/T0HPX5HyF9YLBY0Nja6vW3ttTEYDAgODla4deSv9DKWeCvQAG5c0n3s2DF88MEHiImJwZQpU1Q//NZx0OjqU/HTTz+NefPmtf67trYWAwYMcGu9B2sOujWf09LSgOpqW4Tu0PO6+BV1QssjNXpgtVoVG+z0EtT8idFohMlkUmRZejhSYxfoHyb0Tg9jydqStVi9b7VXAg3gYqh56623MH/+fNx+++04ePAgevfurXiDOhMXFwej0XjJUZmamppLjt7YmUwmRXYcRVVFePSTRz1eTpeqq23X2rn2K93T6js1gb5zNRgMWLFiBfbv34/w8HCPvm9w4cIFBVtGgO1JxZmZmW5/gdtoNKK+vl6X36nRihbfqdFLmHSWHsYSbwYawIVQM2HCBBQVFWHlypWYMWOGVxrTleDgYAwfPhwFBQW48847W18vKCjAlClTvLbeoqoijP3jWFwZcyW+/P5Lr61HLWpf/aTVF+nsV2Co9bdaLBbExMQgNzcXa9eu1cXVT/v3d3+jKtJGQ0MD3njjDY+Xs379euTn5+vi6qfZs2dj4cKFOHPmjGrfb5EkqfXqJy2+U6P21U/+8L2hOWlzvBZoABdCjcViwZdffqnpvWrmzZuH6dOnIy0tDSNGjEBeXh4qKirwyCOPuLQcZ+81U1RVhAkbJ2Bo76FYNmYZbnnnFq/dp0a+ewBgvWT5yt2nRu2QsXPnTly8eFHV+9S88sorePfddzW7Tw0AVFRUqLbervA+Nfql1H1q9BJes7Ky8Pvf/16T+9Q88MADqt+nJiwsTPVbVvjDfWp+PezXbo2hit+npqCgwOVGKG3atGk4ffo0nnvuOZw8eRLJycn49NNPMXDgQJeW48y9ZuyBJrlPMrbcvwVHzhxxel5PSOjqsmvPk3p+fj6++eYb1e4onJOT49V1dEXtwdi+M9fLKSghBO9To2NK3adGT/1N7VPN9m24YcMGbNiwQdV1A0BOTo6qdxT+9a9/rdBS9XOfGlfn705APiYhZXxKl89+qo2sxYFhBxBeH46ffvlTBFmCcD7iPIrTijF833CvPPtp065d6N3UhB+Cg3FPh7tU7tq1CU1NvREc/APS0+9xa/n2Zz999tlnSjTXJfajBWrhAE4UeLR49pPaxo0b5+GznzwfS9zl6Rh6znQOpVtLu31MQkCGGvQF8H0nE10GYDqAGgAbADT9+Ho/AL8B8CaAk8q37TiA/gAqAVx6fVbXv3UVn9JNROQ5tZ/S7TllxxKXeDqG/jhuK/bsJ78iOU71IkFAPCCAHwBpowSpRbKdggQgJNsj0CVJgmTwwqHGNkcXOrat7YEHT4OI1WpV/SndRET+SO2ndHu+/5f/X+0rtzwdQ62Sc+NVYIYa4WBAvwzAA2g9QiOabAVoOw/w43ljq3cPbnUVNhhEiIgCk5L7f9XHEk/HUCdn8a2L7L2ls1NORERE5DMYahhoiIiI/EJgnn6yC9BAo+a5VJ4uIyJ/xX2p/gRkqAkPD4d0tYS6u+pgPG1ExN8jIIVJQFjn87REtKAOdYiIiEBQT+U3m2Q2Az8+x6pndHS735nNEoSw3Y8iOrqn2+uw3/lTzTeHFnfB5BVXRIFFi/2MxWJRPWh4eqdypcYSd3g6hjaHN6Me9d1OF5CXdK/ZuQb/u+N/W2+sF2nq/pr5kpMlGJ43HMUPFyO1X6ryjevf3/ZQjssuAyornf2V0+x3IH388cdRWlqK8PBwr78hJUnS5L44gO2ydS1uHc9PUxTIVL+i5se7pCtzubPrxo8f7/X3vNFoRF1dHVJTU/Haa695dIdwJcYSd3k6hhZ+U4jRV4/mJd2OzNs6DymXpzgdaPyB/U2wYsUKVdf72muv4dy5c6p8irJarYiOjsbGjRtRVFQEQD9PzdaKEqEu0Leht/hjbbQI9c3Nzbjhhhtw3333wWw2q/ZIlp49e+KJJ57w+ro68rWHaKotIENNUq+kgAo0ban91Gwt3vT33XcfDhw4gJCQEFUfaBkbG4vFixfjL3/5C0wmk+pPDW5LkiQ0NzfrbtAjmVK1UfNmmo4YjUY0NjZi2rRpeOaZZ1R7/hJg6+cNDQ1ISUlB7969VVlnW2q+x33xqeBaCMhQ8/LYlwMy0ADqp3w1nwouSRKsViv69OmDsWPHqrLOjjZt2oTa2lpNn6ZrsVgQFRWFhx56CLt370ZkZKTDOrR9AKf933bix+93CSFw4MABnlZTmMFgQEpKisPt33bbd2Sfxmg0ora2FjfddBPWrVunmz6npebmZhgMBlX3N/7w1Gx/E5ChJiy4i28Ek6K0eGiaEEL1Qdg+CBmNRs137nb5+fmKLCc2NhZnzpzpdKAl59m3YWxsLIqLixVbrl76nP3hjmr3E4PBgB49eqi6TtKngAw1uhQf3/6/zv2KHND6E5QeBn77UauuPvU7swwhBBobGwHo4+/ydfZt2NDQ4NGpYHtt7Kck9FAbrd931L1AGEsYavRi3z53fkU6pPYVV51R6lTjt99+C6vVqpu/y9fZj+gpeSqYtSFn+PJYUvZDmVPTMdQQUZfi/fljHRHpXlFVEeZtnefUtAw1RNQlPZza8Ec8ukLUvaKqIoz941gk9UrCIRzqdnqGGiLqEgdfItKCPdAk90lG1vVZmIiJ3c7Di96JiIhIV9oGmi33b3H6qmWGGiIiItKNjoHGlfvKMdQQERGRLngSaAAfCTVHjx7FzJkzkZSUhNDQUFx55ZXIzs5GU1OT1k0jIiIiBXgaaAAf+aLw119/DavVijfffBNXXXUVDh48iNmzZ6O+vh4vv/yy1s0jIiIiDygRaAAfCTUTJkzAhAkTWv99xRVX4PDhw1i9ejVDDRERkQ9TKtAAPhJqHDGbzYiJielymsbGxtZbvANAbW0tAODw6cOIOBnh0vrKTpW1+y8RERE5p7Mx9GDNQTz6yaO4MuZKPH/b8zhy5ojD+Q+fPuzUeiThg3fW+u6775CamopXXnkFs2bN6nS6nJwcLFq06NJfLAAQ4r32ERERkYIaADxvO6DR1QNcNQ01nYaONvbu3Yu0tLTWf584cQKjR4/G6NGjsXbt2i7ndXSkZsCAAXj202dxV8pdLrW17FQZHvjwAbw1+S2sKFqB7858h1WTViG5T7JLy1HK2pK1WL1vNeakzcGs1M6DnTe1TdgrJ65EeHC46m2ob6pH5uZM1gOsR1ushw3rIWM9ZFrUwz6GbrhrA4bEDcHHX3+MRYWLEBYUhg/v+RC9I3p3Of+HpR9iyc+XdBtqIDT0ww8/iLKysi5/Ll682Dp9VVWVGDx4sJg+fbqwWCwur89sNgsAIm9XnsvzFp8oFsiBGLZ6mIhaFiX2VO5xeRlKyd2eK5ADkbs9V7M27KncI6KWRYn0demitqFWkzbUNtSK9HXprIdgPdpiPWxYDxnrIdOqHvYxtPhEscgvyRfIgYhcGimqzFVOzZ+3K08AEGazucvpNA01rqisrBSDBg0S9957r2hpaXFrGZ6Emi1HtgjkQIQvCQ/IDtkWdxAy1sOG9ZCxHjLWw4b1kEPN77f93uVAI4SfhZqqqipx1VVXiTFjxojKykpx8uTJ1h9XuBtqqsxVImxJmEAOxNulb7s0r5K4g7DhDkLGeshYDxvWQ8Z6yLSuhz3UuBNohPCzULN+/XoBwOGPK9wJNVXmKhG5NLK1GMUnil1tviK07pBCcAfRFuthw3rIWA8Z62HDesjeLn1bIAcibEmYy4FGCD8LNUpxNdS0DTT2Q2ZahBo9dEjuIGSshw3rIWM9ZKyHDesh21O5R4QvCRfIgdhyZItby2CoccCVUNM20OSX5Lf7kpOa9NIhuYOwYT1sWA8Z6yFjPWxYD5m9HsNWD/NoDGWoccDZUNMx0AghNAk1euqQ3EGwHnash4z1kLEeNqyHrG09Co8WMtQozZlQ4yjQCKF+qNFbh+QOgvUQgvVoi/WQsR42rIesYz08HUMZahzoLtR0FmiEUDfU6LFDaoE7CBnrIWM9bFgPGesh02s9GGq8oKtQ01WgEUK9UKPXDqk27iBkrIeM9bBhPWSsh0zP9WCo8YLOQk13gUYIdUKNnjukmriDkLEeMtbDhvWQsR4yvdeDocYLHIUaZwKNEN4PNXrvkGrhDkLGeshYDxvWQ8Z6yHyhHgw1XtAx1DgbaITwbqjxhQ6pBu4gZKyHjPWwYT1krIfMV+qhVqgJcvFBmz5N/PhA8h3f7MC5i+eQtS0LjS2NyEjJQEtDC97a/Van8x47dwxoAIqPFqPufJ1ibXrnwDtYV7IOM1NnYlT8KBR+U6jYsp1V9kMZ5m2dh6ReSci6Pgv7j+1XvQ0Xmi7gqYKnUH62HMvHL0dDfYMm24L1sGE9ZKyHjPWwYT1kztbj8OnDQIPtadvFR4tdXs+Ob3YAkMfxzkiiuyn8SGVlJQYMGKB1M4iIiMgNx48fR//+/Tv9fUCFGqvVihMnTiAyMhKSJLk0b21tLQYMGIDjx48jKirKSy30P9xuruM2cw+3m+u4zdzD7eY6T7eZEALnz59HQkICDAZDp9MF1Okng8HQZcJzRlRUFDuxG7jdXMdt5h5uN9dxm7mH2811nmyz6OjobqfpPO4QERER+RCGGiIiIvILDDVOMplMyM7Ohslk0ropPoXbzXXcZu7hdnMdt5l7uN1cp9Y2C6gvChMREZH/4pEaIiIi8gsMNUREROQXGGqIiIjILzDUEBERkV9gqHHCkiVLkJ6ejrCwMPTs2dPhNJIkXfKzZs0adRuqM85st4qKCkyePBnh4eGIi4vDb3/7WzQ1NanbUJ27/PLLL+lbCxYs0LpZurJq1SokJSUhJCQEw4cPx44dO7Rukq7l5ORc0qfi4+O1bpauFBYWYvLkyUhISIAkSfjoo4/a/V4IgZycHCQkJCA0NBS33HILDh06pE1jdaS77fbggw9e0vduvPFGxdbPUOOEpqYm3H333ZgzZ06X061fvx4nT55s/cnIyFCphfrU3XazWCyYNGkS6uvrsXPnTvz5z3/GBx98gPnz56vcUv177rnn2vWthQsXat0k3Xjvvfcwd+5cPPvss9i/fz9GjhyJiRMnoqKiQuum6drQoUPb9amvvvpK6ybpSn19PYYNG4aVK1c6/P2LL76I5cuXY+XKldi7dy/i4+MxduxYnD9/XuWW6kt32w0AJkyY0K7vffrpp8o1wK1ngAeo9evXi+joaIe/AyD++te/qtoeX9HZdvv000+FwWAQVVVVra+9++67wmQydft4+UAycOBA8Yc//EHrZujW9ddfLx555JF2r11zzTViwYIFGrVI/7Kzs8WwYcO0bobP6Lh/t1qtIj4+Xjz//POtrzU0NIjo6GixZs0aDVqoT47GxYyMDDFlyhSvrZNHahSUmZmJuLg4/OxnP8OaNWtgtVq1bpKu7d69G8nJyUhISGh9bfz48WhsbERxseuPpvdnL7zwAmJjY5GSkoIlS5bwFN2PmpqaUFxcjHHjxrV7fdy4cdi1a5dGrfINR44cQUJCApKSknDvvffiv//9r9ZN8hnl5eWorq5u1+9MJhNGjx7NfueEL774An369MHgwYMxe/Zs1NTUKLbsgHqgpTfl5ubitttuQ2hoKP71r39h/vz5OHXqFE8TdKG6uhp9+/Zt91qvXr0QHByM6upqjVqlP0888QRSU1PRq1cvFBUV4emnn0Z5eTnWrl2rddM0d+rUKVgslkv6Ud++fdmHunDDDTfgnXfeweDBg/H9999j8eLFSE9Px6FDhxAbG6t183TP3rcc9btjx45p0SSfMXHiRNx9990YOHAgysvLkZWVhTFjxqC4uFiRuw0H7JEaR1+U6/izb98+p5e3cOFCjBgxAikpKZg/fz6ee+45vPTSS178C7Sh9HaTJOmS14QQDl/3J65sxyeffBKjR4/GT3/6U8yaNQtr1qzBunXrcPr0aY3/Cv3o2F8CoQ95YuLEiZg6dSquvfZa3H777fjkk08AAG+//bbGLfMt7HeumzZtGiZNmoTk5GRMnjwZmzdvxjfffNPaBz0VsEdqMjMzce+993Y5zeWXX+728m+88UbU1tbi+++/vyTN+zIlt1t8fDz27NnT7rWzZ8+iubnZr7aZI55sR/uVAt9++23Af6qOi4uD0Wi85KhMTU2N3/chJYWHh+Paa6/FkSNHtG6KT7BfKVZdXY1+/fq1vs5+57p+/fph4MCBivW9gA01cXFxiIuL89ry9+/fj5CQkE4vZfZVSm63ESNGYMmSJTh58mTrjuGzzz6DyWTC8OHDFVmHXnmyHffv3w8A7XamgSo4OBjDhw9HQUEB7rzzztbXCwoKMGXKFA1b5lsaGxtRVlaGkSNHat0Un5CUlIT4+HgUFBTguuuuA2D7ftf27dvxwgsvaNw633L69GkcP35csf1ZwIYaV1RUVODMmTOoqKiAxWJBaWkpAOCqq65CREQE/v73v6O6uhojRoxAaGgotm3bhmeffRYPP/xwQD/FtbvtNm7cOPzkJz/B9OnT8dJLL+HMmTN46qmnMHv2bERFRWnbeJ3YvXs3/v3vf+PWW29FdHQ09u7diyeffBJ33HEHEhMTtW6eLsybNw/Tp09HWloaRowYgby8PFRUVOCRRx7Rumm69dRTT2Hy5MlITExETU0NFi9ejNra2oC/DUVbdXV1+Pbbb1v/XV5ejtLSUsTExCAxMRFz587F0qVLMWjQIAwaNAhLly5FWFgYfvWrX2nYau11td1iYmKQk5ODqVOnol+/fjh69CieeeYZxMXFtftQ4hGvXVflRzIyMgSAS362bdsmhBBi8+bNIiUlRURERIiwsDCRnJwsXn31VdHc3KxtwzXW3XYTQohjx46JSZMmidDQUBETEyMyMzNFQ0ODdo3WmeLiYnHDDTeI6OhoERISIq6++mqRnZ0t6uvrtW6arrzxxhti4MCBIjg4WKSmport27dr3SRdmzZtmujXr5/o0aOHSEhIEHfddZc4dOiQ1s3SlW3btjncf2VkZAghbJd1Z2dni/j4eGEymcSoUaPEV199pW2jdaCr7XbhwgUxbtw40bt3b9GjRw+RmJgoMjIyREVFhWLrl4QQQpl4RERERKSdgL36iYiIiPwLQw0RERH5BYYaIiIi8gsMNUREROQXGGqIiIjILzDUEBERkV9gqCEiIiK/wFBDREREfoGhhoh8ksViQXp6OqZOndrudbPZjAEDBmDhwoUatYyItMI7ChORzzpy5AhSUlKQl5eH+++/HwAwY8YMHDhwAHv37kVwcLDGLSQiNTHUEJFPe/3115GTk4ODBw9i7969uPvuu1FUVISUlBStm0ZEKmOoISKfJoTAmDFjYDQa8dVXX+Hxxx/nqSeiAMVQQ0Q+7+uvv8aQIUNw7bXXoqSkBEFBQVo3iYg0wC8KE5HPy8/PR1hYGMrLy1FZWal1c4hIIzxSQ0Q+bffu3Rg1ahQ2b96MF198ERaLBf/85z8hSZLWTSMilfFIDRH5rIsXLyIjIwO/+c1vcPvtt2Pt2rXYu3cv3nzzTa2bRkQaYKghIp+1YMECWK1WvPDCCwCAxMREvPLKK/jd736Ho0ePats4IlIdTz8RkU/avn07brvtNnzxxRe4+eab2/1u/PjxaGlp4WkoogDDUENERER+gaefiIiIyC8w1BAREZFfYKghIiIiv8BQQ0RERH6BoYaIiIj8AkMNERER+QWGGiIiIvILDDVERETkFxhqiIiIyC8w1BAREZFfYKghIiIiv8BQQ0RERH7h/wDksr9ueIjkAwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "# Adds mode monitors at each of the ports to track the energy that goes in or out\n", "\n", @@ -336,8 +372,6 @@ "mode_monitor_3 = sim.add_mode_monitor(fcen, 0, 1, mp.ModeRegion(volume=port3))\n", "mode_monitor_4 = sim.add_mode_monitor(fcen, 0, 1, mp.ModeRegion(volume=port4))\n", "\n", - "sim.load_minus_flux_data(mode_monitor_1, subtraction_data)\n", - "\n", "# Plot the simulation\n", "plot_plane = mp.Volume(center=mp.Vector3(z=0), size=mp.Vector3(cell.x, cell.y, 0))\n", "sim.plot2D(output_plane=plot_plane if sim_is_3D else None) # No parameters are needed for a 2D simulation. \n" @@ -348,7 +382,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "As we see from the output of sim.plot2D, our simulation is set up correctly. The red is our source, the blue are our 4 mode monitors, an the black is the geometry. We are ready to run the simulation! Actually running the simulation is the most computationally intense part of this, so it may take some time. The until_after_sources parameter for sim.run() means the run the simulation until 100 meep time units after the sources have turned off. This makes sure the all of the light has time to propagate through the mmi." + "As we see from the output of sim.plot2D, our simulation is set up correctly. The red line is our source, the blue are our 4 monitors. We are ready to run the simulation! Actually running the simulation is the most computationally intense part of this, so it may take some time if you have a high resolution and are running a 3D simulation. The `until_after_sources=mp.stop_when_dft_decayed` parameter for sim.run() means meep will start checking how much the cumulative Fourier Transform of all the fields in the simulation changes only after the sources have turned off. Then, it will wait until the DFT's change over one time step is below a certain threshold (1e-11 relative change by default), at which point it will stop the simulation" ] }, { @@ -370,7 +404,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now that our simulation has been run, we'll quickly compute the s-parameters, store them in our array, and print them out. The function sim.get_eigenmode_coefficients() returns an array that holds the parameters for every mode, every frequency, and in the forward and backward direction. In our simulation we only simulated the fundamental TE mode and at only our center frequency. For most applications, only the fundamental TE mode matters, but often we do want to see how the S-parameters change with frequency. To do that, simply change the df and nfreq parameters for the mode monitors. However, in our simulation we do care about the direction. In our simulation, light will go in port 1 going forward (towards positive X) and exit ports 3 and 4 in the same direction, but any light exiting port 1 and 2 will be going backwards (towards negative X). So, we change the alpha to be (0,0,1) for port 1 and 2. The 1 in the third spot indicates the backward propagating light. We'll also see how much light is inserted at port 1 going forward." + "Now that our simulation has been run, the monitors contain the s-parameters which we can retrieve. The function sim.get_eigenmode_coefficients() returns an array that holds the parameters for every mode, every frequency, and in the forward and backward direction. In our simulation we only simulated the fundamental TE mode and at only our center frequency. For most applications, only the fundamental TE mode matters, but often we do want to see how the S-parameters change with frequency. To do that, simply change the df and nfreq parameters for the mode monitors. However, in our simulation we do care about the direction. In our simulation, light will go in port 1 going forward (towards positive X) and exit ports 3 and 4 in the same direction, but any light exiting port 1 and 2 will be going backwards (towards negative X). Therefore, we retrieve the index [0,0,1] for port 1 and 2 which is the index for backward propagating mode." ] }, { @@ -380,10 +414,11 @@ "outputs": [], "source": [ "# Finds the S parameters\n", - "port1_coeff = sim.get_eigenmode_coefficients(mode_monitors[0], [1], eig_parity=mode_parity, direction=-mp.X, ).alpha[0, 0, 1] / norm_mode_coeff\n", - "port2_coeff = sim.get_eigenmode_coefficients(mode_monitors[1], [1], eig_parity=mode_parity, direction=-mp.X).alpha[0, 0, 1] / norm_mode_coeff\n", - "port3_coeff = sim.get_eigenmode_coefficients(mode_monitors[2], [1], eig_parity=mode_parity).alpha[0, 0, 0] / norm_mode_coeff\n", - "port4_coeff = sim.get_eigenmode_coefficients(mode_monitors[3], [1], eig_parity=mode_parity).alpha[0, 0, 0] / norm_mode_coeff\n", + "norm_mode_coeff = sim.get_eigenmode_coefficients(mode_monitor_1, [1], eig_parity=mode_parity).alpha[0, 0, 0]\n", + "port1_coeff = sim.get_eigenmode_coefficients(mode_monitor_1, [1], eig_parity=mode_parity).alpha[0, 0, 1] / norm_mode_coeff\n", + "port2_coeff = sim.get_eigenmode_coefficients(mode_monitor_2, [1], eig_parity=mode_parity).alpha[0, 0, 1] / norm_mode_coeff\n", + "port3_coeff = sim.get_eigenmode_coefficients(mode_monitor_3, [1], eig_parity=mode_parity).alpha[0, 0, 0] / norm_mode_coeff\n", + "port4_coeff = sim.get_eigenmode_coefficients(mode_monitor_4, [1], eig_parity=mode_parity).alpha[0, 0, 0] / norm_mode_coeff\n", "# Store the S parameters in s_params\n", "s_params[0] = [port1_coeff, port2_coeff, port3_coeff, port4_coeff]\n", "\n", @@ -400,9 +435,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Just for fun, we'll also compute and print the transmittance through each port and the total insertion loss. We won't do this with the other ports, but it's a good exercise here. \n", - "\n", - "We find that this component is actually terrible and would never be used in an actual photonic design. Almost a third of the light is lost, and the light that isn't lost is not split in any specific ratio. Fortunately, it was never meant to be used, and just exists as an example of the basic shape of an mmi." + "Just for fun, we'll also compute and print the transmittance through each port and the total insertion loss. We won't do this with the other ports as inputs, but it's a good exercise here. " ] }, { @@ -434,6 +467,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "We find that this component is actually terrible and would never be used in an actual photonic design. Almost a third of the light is lost, and the light that isn't lost is not split in any specific ratio. Fortunately, it was never meant to be used, and just exists as an example of the basic shape of an mmi.\n", + "\n", "Finally, before we move on to the next step, we'll run the simulation again and create a plot of the steady state of the mmi. " ] },