{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Adjoint analysis of a multi-layer slab\n", "\n", "In this notebook, we will show how to differentiate with respect to\n", "`DiffractionMonitor` outputs in `tidy3d` and also check the gradient\n", "values against gradients obtained using transfer matrix method (TMM),\n", "validating their accuracy for a multilayer slab problem." ], "id": "534ec54a-4d74-4c25-8627-ca59eb77475c" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from typing import List, Tuple\n", "\n", "import autograd as ag\n", "import autograd.numpy as anp\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import tidy3d as td\n", "import tidy3d.web as web\n", "import tmm" ], "id": "0e0fafb0" }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we define some global parameters describing the transmission\n", "through a multilayer slab with some spacing between each slab.\n", "\n", "The layout is diagrammed below.\n", "\n", "\"Schematic" ], "id": "126faeb9-29c5-4b14-8f0b-f39b7dc76ba4" }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# frequency we want to simulate at\n", "freq0 = 2.0e14\n", "k0 = 2 * np.pi * freq0 / td.C_0\n", "freqs = [freq0]\n", "wavelength = td.C_0 / freq0\n", "\n", "# background permittivity\n", "bck_eps = 1.0**2\n", "\n", "# space between each slab\n", "spc = 0.0\n", "\n", "# slab permittivities and thicknesses\n", "slab_eps0 = [2.0**2, 1.8**2, 1.5**2, 1.9**2]\n", "slab_ds0 = [0.5, 0.25, 0.5, 0.5]\n", "\n", "# incidence angle\n", "theta = 0 * np.pi / 8\n", "\n", "# resolution\n", "dl = 0.01" ], "id": "28b7b516" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Transfer Matrix Method (Ground Truth)\n", "\n", "Next we use the `tmm` package to write a function to return the\n", "transmission `T` of `p` polarized light given a set of slab\n", "permittivities and thicknesses. We’ll also write a function to compute\n", "the numerical gradient using TMM and will take these to be our “ground\n", "truths” when evaluating the accuracy of our values obtained through FDTD\n", "and autograd.\n", "\n", "### Transmission Calculation with TMM\n", "\n", "First, we write a function to compute transmission." ], "id": "9fb22f4a-b933-484f-8486-df11328236d3" }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def compute_T_tmm(slab_eps=slab_eps0, slab_ds=slab_ds0) -> float:\n", " \"\"\"Get transmission as a function of slab permittivities and thicknesses.\"\"\"\n", "\n", " # construct lists of permittivities and thicknesses including spaces between\n", " new_slab_eps = []\n", " new_slab_ds = []\n", " for eps, d in zip(slab_eps, slab_ds):\n", " new_slab_eps.append(eps)\n", " new_slab_eps.append(bck_eps)\n", " new_slab_ds.append(d)\n", " new_slab_ds.append(spc)\n", " slab_eps = new_slab_eps[:-1]\n", " slab_ds = new_slab_ds[:-1]\n", "\n", " # add the input and output spaces to the lists\n", " eps_list = [bck_eps] + slab_eps + [bck_eps]\n", " n_list = np.sqrt(eps_list)\n", " d_list = [np.inf] + slab_ds + [np.inf]\n", "\n", " # compute transmission with TMM\n", " return tmm.coh_tmm(\"p\", n_list, d_list, theta, wavelength)[\"T\"]" ], "id": "c99edb28" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We run this function with our starting parameters and see that we get a\n", "transmission of about 98% for the set of input parameters." ], "id": "28a760e6-4e54-4762-b8c4-034268054195" }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "T (tmm) = 0.786" ] } ], "source": [ "T_tmm = compute_T_tmm(slab_eps=slab_eps0, slab_ds=slab_ds0)\n", "print(f\"T (tmm) = {T_tmm:.3f}\")" ], "id": "1f051167" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Numerical Gradient with TMM\n", "\n", "Next, we will use our `compute_T_tmm()` function to compute the\n", "“numerical” gradient to use as a comparison against our adjoint results\n", "with FDTD.\n", "\n", "The derivative of a function $f(x)$ w.r.t. $x$ can be approximated using\n", "finite differences as\n", "\n", "$$ \\frac{df}{dx}(x) \\approx \\frac{f(x+\\Delta) - f(x-\\Delta)}{2\\Delta}$$\n", "\n", "with a small step $\\Delta$.\n", "\n", "To compute the gradient of our transmission with respect to each of the\n", "slab thicknesses and permittivities, we need to repeat this step for\n", "each of the values. Luckily, since TMM is very fast, we can compute\n", "these quantities quite quickly compared to if we were using FDTD.\n", "\n", "> Important note: We assume in our TMM numerical gradient that when the\n", "> slabs are touching (`spc=0`) and a slab thickness is modified, the\n", "> thicknesses of the neighboring slabs adjust to accommodate this\n", "> change. For example, if slab `i` increases by `dt`, slab `i-1` and\n", "> `i+1` each decrease by `dt/2`. We also account for this in our FDTD\n", "> set up by keeping the centers of all boxes constant and not tracking\n", "> the gradient through these quantities. The reason this is required is\n", "> that `tidy3d` does not recognize the space between touching `td.Box`\n", "> objects as a single interface and will instead “double count” the\n", "> gradient contribution of the interface if they are placed right next\n", "> to each other. One must therefore be careful about overlapping or\n", "> touching two `td.Box` or other geometries when computing gradients.\n", "\n", "Here we write the function to return the numerical gradient." ], "id": "b048ec11-4581-47b7-9742-9425b7018baf" }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def compute_grad_tmm(slab_eps=slab_eps0, slab_ds=slab_ds0) -> Tuple[List[float], List[float]]:\n", " \"\"\"Compute numerical gradient of transmission w.r.t. each of the slab permittivities and thicknesses using TMM.\"\"\"\n", "\n", " delta = 1e-4\n", "\n", " # set up containers to store gradient and perturbed arguments\n", " num_slabs = len(slab_eps)\n", " grad_tmm = np.zeros((2, num_slabs), dtype=float)\n", " args = np.stack((slab_eps, slab_ds), axis=0)\n", "\n", " # loop through slab index and argument index (eps, d)\n", " for arg_index in range(2):\n", " for slab_index in range(num_slabs):\n", " grad = 0.0\n", "\n", " # perturb the argument by delta in each + and - direction\n", " for pm in (-1, +1):\n", " args_num = args.copy()\n", " args_num[arg_index][slab_index] += delta * pm\n", "\n", " # NEW: for slab thickness gradient, need to modify neighboring slabs too\n", " if arg_index == 1 and spc == 0:\n", " if slab_index > 0:\n", " args_num[arg_index][slab_index - 1] -= delta * pm / 2\n", " if slab_index < num_slabs - 1:\n", " args_num[arg_index][slab_index + 1] -= delta * pm / 2\n", "\n", " # compute argument perturbed T and add to finite difference gradient contribution\n", " T_tmm = compute_T_tmm(slab_eps=args_num[0], slab_ds=args_num[1])\n", " grad += pm * T_tmm / 2 / delta\n", "\n", " grad_tmm[arg_index][slab_index] = grad\n", " grad_eps, grad_ds = grad_tmm\n", " return grad_eps, grad_ds" ], "id": "48758668" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let’s run this function and observe the gradients. These will be saved\n", "later to compare against our autograd results." ], "id": "4c612af4-94f3-4290-8df9-0152d894a6ea" }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "gradient w.r.t. eps (tmm) = [-0.2766323 0.01377339 -0.2032054 -0.28999361]\n", "gradient w.r.t. ds (tmm) = [-1.75199732 -0.21552416 1.00729645 -2.08209951]" ] } ], "source": [ "grad_eps_tmm, grad_ds_tmm = compute_grad_tmm()\n", "print(f\"gradient w.r.t. eps (tmm) = {grad_eps_tmm}\")\n", "print(f\"gradient w.r.t. ds (tmm) = {grad_ds_tmm}\")" ], "id": "4372a8c2" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## FDTD\n", "\n", "Next, we will implement the same two functions using Tidy3D.\n", "\n", "### Transmission Calculation with FDTD\n", "\n", "We first write a function to compute the transmission of a multilayer\n", "slab using Tidy3D. We use a `DiffractionMonitor` to measure our\n", "transmission amplitudes. Below, we break up the transmission calculation\n", "into a few functions to make it easier to read and reuse later." ], "id": "d10f929f-344e-429b-896a-87d41b0f9337" }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "from autograd.tracer import getval\n", "\n", "\n", "def make_sim(slab_eps=slab_eps0, slab_ds=slab_ds0) -> td.Simulation:\n", " \"\"\"Create a Simulation given the slab permittivities and thicknesses.\"\"\"\n", "\n", " # frequency setup\n", " fwidth = freq0 / 10.0\n", " freqs = [freq0]\n", "\n", " # geometry setup\n", " bck_medium = td.Medium(permittivity=bck_eps)\n", "\n", " space_above = 2\n", " space_below = 2\n", "\n", " length_x = 1.0\n", " length_y = 1.0\n", " length_z = space_below + sum(slab_ds0) + space_above + (len(slab_ds0) - 1) * spc\n", " sim_size = (length_x, length_y, length_z)\n", "\n", " # make structures\n", " slabs = []\n", " z_start = -length_z / 2 + space_below\n", " for d, eps in zip(slab_ds, slab_eps):\n", " # don't track the gradient through the center of each slab\n", " # as tidy3d doesn't have enough information to properly process the interface between touching Box objects\n", " z_center = z_start + d / 2\n", " z_center = getval(z_center)\n", "\n", " slab = td.Structure(\n", " geometry=td.Box(center=[0, 0, z_center], size=[td.inf, td.inf, d]),\n", " medium=td.Medium(permittivity=eps),\n", " )\n", " slabs.append(slab)\n", " z_start = z_start + d + spc\n", "\n", " # source setup\n", " gaussian = td.GaussianPulse(freq0=freq0, fwidth=fwidth)\n", " src_z = -length_z / 2 + space_below / 2.0\n", "\n", " source = td.PlaneWave(\n", " center=(0, 0, src_z),\n", " size=(td.inf, td.inf, 0),\n", " source_time=gaussian,\n", " direction=\"+\",\n", " angle_theta=theta,\n", " angle_phi=0,\n", " pol_angle=0,\n", " )\n", "\n", " # boundaries\n", " boundary_x = td.Boundary.bloch_from_source(\n", " source=source, domain_size=sim_size[0], axis=0, medium=bck_medium\n", " )\n", " boundary_y = td.Boundary.bloch_from_source(\n", " source=source, domain_size=sim_size[1], axis=1, medium=bck_medium\n", " )\n", " boundary_spec = td.BoundarySpec(x=boundary_x, y=boundary_y, z=td.Boundary.pml(num_layers=40))\n", "\n", " # monitors\n", " mnt_z = length_z / 2 - space_above / 2.0\n", " monitor_1 = td.DiffractionMonitor(\n", " center=[0.0, 0.0, mnt_z],\n", " size=[td.inf, td.inf, 0],\n", " freqs=freqs,\n", " name=\"diffraction\",\n", " normal_dir=\"+\",\n", " )\n", "\n", " # make simulation\n", " return td.Simulation(\n", " size=sim_size,\n", " # grid_spec=td.GridSpec.auto(min_steps_per_wvl=100),\n", " grid_spec=td.GridSpec.uniform(dl=0.01),\n", " structures=slabs,\n", " sources=[source],\n", " monitors=[monitor_1],\n", " run_time=10 / fwidth,\n", " boundary_spec=boundary_spec,\n", " medium=bck_medium,\n", " subpixel=True,\n", " shutoff=1e-8,\n", " )" ], "id": "b9533b10" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let’s generate a simulation and plot it to make sure it looks\n", "reasonable." ], "id": "fd9a46c2-3437-4d94-b1f6-e3ae935f12b0" }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgUAAAZoCAYAAAACubaaAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90\nbGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAd\nhwAAHYcBj+XxZQAAiwFJREFUeJzs3Xd4VFX+x/HPpFdSqNJBqooiKiAWUCygVAuoq1hAcO26WFDU\ntaz6W3HXdXV1FcQCWAALgigqVQWkCQiIdEKVkpDez+8P7swmZJJJmT7v1/PkccLce875jsk9n9xq\nM8YYAQCAkBfm6wEAAAD/QCgAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUA\nAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFQAUH\nDhyQzWaTzWZTz549fT0cAGVkZGSoYcOGstlsevjhh309HI/6+uuvHduiRYsWeaVPQgGC2tatWx2/\nVH369PH1cIIen7dr69at03333adTTz1VSUlJSkhIUIcOHXT77bfrxx9/DJh+fVXHE088ocOHDys5\nOVnjxo3zWD/+oF+/fo7fo3vvvVclJSWe79QAQWzLli1GkpFkevfuXa119u/f71inR48enh1gkKnN\n5x0qSktLzfjx4014eLjjM3L2NXr0aFNYWOi3/fqqDmOM2bp1q4mIiDCSzOOPP+7Wtv3Vd9995/hM\nJ0+e7PH+CAUIaoQC7yIUVO6RRx6pchIt+3XjjTea0tJSv+zXV3UYY8xNN91kJJno6Ghz8OBBt7Xr\n77p162YkmTZt2piioiKP9kUoQFBjkvIuPm/nFi9eXG6y7NKli5k1a5Y5evSoyczMNPPnzzcXXnhh\nuWWmTZvmd/36qg5jjElLS3PsJbjhhhvc0mag+O9//+v4PD/88EOP9kUoQFBjkvIuPm/nevXqVW4i\nzc7OrrBMUVGRufzyyx3LtW7d2pSUlPhVv76qwxhjxo0b52jzu+++q3N7gSQzM9PExMQYSebcc8/1\naF9+GwoKCwvNe++9Z4YPH27atm1rEhMTTWxsrGnTpo0577zzzGOPPWaWL1/udNfUU0895fjhefvt\nt40xxuzatcs88cQT5swzzzQNGjQwNpvN9O3bt8K62dnZ5pVXXjGXXHKJadq0qYmKijIpKSnmjDPO\nMA8++KDZuHFjtWuYN2+eGTFihOncubNJTEw0ERERplGjRubUU081/fr1MxMmTHDZnjva8OZ47b77\n7jtz++23m06dOpnk5GQTHR1tmjVrZgYOHGgmT55c411gixYtMnfddZfp0qWLqV+/vomKijJNmzY1\n3bp1MyNHjjRffPGFyc3NdSy/YMGCau/ivPrqq8v1VdPDB+74mVmxYoWjz8GDBxtjjv8OvP3226ZP\nnz6mcePGjpqvvvpqM3/+/Bp9flXZt2+f+eijj8yf//xnc+6555q2bduahIQEExsba5o3b24GDRpk\nJk6cWOXx4bp83q7k5eWZ1NRUI8lERESYffv2VWu99957r9Z9utPWrVvL1f/TTz9VuuyuXbscfw1L\nMgsWLPCbfn1VhzHGlJSUmObNmxtJpn79+qa4uLjSZd9++21Hv88++6zLtkeOHOlY/ttvv63wvrPf\nzdzcXPOvf/3L9OrVyzRq1MjExsaaTp06mfvvv9/s2bOnQhvffvutufbaa03btm1NTEyMadCggRkw\nYECNws2gQYMc49i8eXO116spvwwF3333nWnVqlW1NjArVqyosP6JoeCTTz4x9erVq7DuxRdfXG69\nefPmmSZNmlTZX3h4uHnooYeqTL7Z2dnl/ge6+kpPT/dIG9Xlzr727dtnLrroIpdtdO7cuVo/2Hv3\n7jUXX3xxtcb11FNPOdbzVihw18/MiRue3bt3m3POOafKdh999FGXn191NG7cuFqfU7t27cymTZuc\ntuHJUGCMMQ888IBj/RdeeKFa65Tdjf3111/XuE93mTBhgmMcp59+usvlBw8e7Fj+3nvv9Zt+fVWH\nMcf/KLC3deONN1a5rKdDwZYtW8ypp55a6c93/fr1zerVq40xxhQUFJRr39nXiy++WK3PYOLEiY51\nnnnmmWqtUxt+d0nihx9+qH79+mnXrl1uae+HH37Qn/70J2VmZlZ4r7S01PF67ty5uvLKK3XgwIEq\n2yspKdFLL72kkSNHVrrMPffco1mzZtV+0G5qw9t97dq1Sz169NCCBQtcLrtp0yadd9552r17d6XL\nbNu2TWeffbbmz59f57F5gjt/ZsrKzMxU//79tWLFiiqXe/HFFzV9+vRqj7eutm7dqt69e+uPP/7w\nWp92o0ePdryeNGmSy+W3bNmixYsXS5Jat26tSy+91GNjc2XNmjWO15dcconL5S+77DKn6/q6X1/V\nIUlfffWV47UvL3U9duyYBgwYoA0bNlS6zJEjRzRw4EDl5+dr9OjRLn9ex40bpyVLlrjsu3fv3o7X\nZT8Pt/NY3KiF9evXm9jYWEcaiouLM2PHjjXLli0z6enpJj8/3+zYscP88MMP5rnnnjPdu3c3q1at\nqtBO2T0F9q/u3bubzz77zPzxxx8V/mI7dOiQadSokWPZ+Ph48/TTT5vNmzeb/Px888cff5iPP/7Y\ndO7cuVybU6dOrdD30aNHy+02u+mmm8yiRYvM4cOHTXFxsUlPTzebNm0y8+bNMw899JDp3LmzycjI\ncHsb1eWuvoqKikz37t0d7aSkpJhnnnnGrFmzxmRnZ5vCwkKzY8cO8/rrr5f7y7pPnz5Ox1VQUGDO\nOOOMcp/34MGDzezZs83BgwdNYWGh2b9/v1m1apWZOHGiGThwoPm///u/Cu146uoDd/7MGFP+rxH7\nV7Nmzcx//vMfs2PHDlNYWGgOHz5sPvzwQ9OsWTPHMl26dKlWTVU5++yzza233mo+/PBDs2rVKpOW\nlmYKCgpMdna2Wbt2rXnmmWdMUlKSo89bb7210rY8eU5B7969HW0vXLiwymXLniH/3HPPOV2m7HFv\nd3x99tlnTvvp2rWrY5l33nnHZZ0LFy50LJ+amupy+cq4u19f1WGMMT169HC0tX79+iqX9eSeAvtX\n27ZtzXvvvWf27dtnCgsLzc6dO81TTz1lwsLCHMtccsklRjp+pcT48ePNhg0bTH5+vklPTzfTp083\nTZs2dSzbr1+/an0ODRo0MJJMWFiYycrKqtY6NeVXoWDgwIGOD+mkk05y+T+/MieGgqFDh1Z5DPvp\np58uF0SWL1/udLmsrCzTs2dPx7Lt27evcE7DsmXLHO9feeWVtRq/O9rwdl/vv/++o51OnTqZvXv3\nVrrsgQMHTNu2bR3L//jjjxWWKXu2rc1mq/X1uZ4KBe78mTGm4oanY8eO5o8//nDa5m+//WaioqIc\ny27durVaddXF77//blJSUowkExMTU2kI9WQomDZtmqPtqnYhFxUVOYJnREREpT+L3goFZUPw999/\n77LOnTt3lmu3ttf6u7tfX9VRVFTkOMkuMjLSZTueDgVdunSp9BDqE088UW7ZqKgos3jxYqfLrlmz\nxthsNiMdP8SYk5PjcqxlD6UuWrTI5fK14TeHD/bs2aPZs2c7vv/ggw902mmn1bnd+Ph4TZw4URER\nEZUu89577zleP/XUU+revbvT5RISEjRlyhRHW1u2bKlw567o6GjH60aNGtVqzO5ow9t9/ec//3G8\nfu+999S0adNKl23cuLFeeuklx/dffPFFhWXeeOMNx+v77rtPt9xyS63H5gnu/Jlx5p133lHDhg2d\nvtexY0cNGDDA8f0vv/xSg5HXTvv27XX//fdLkvLz87V06VKP93miq666Sg0aNJAkzZw5U8eOHXO6\n3Jw5cxyHdK688soqfxa9ISsry/E6MTHR5fIJCQmVru/Lfn1Vx86dO5Wfny9Jat68uSIjI2vVjrtM\nnDhRycnJTt+74447yn1///3364ILLnC6bNeuXdWjRw9Jxw8x/vrrry77btOmjeP1b7/9Vs0R14zf\nhIL58+fLGCNJOvvss9W3b1+3tDtgwAClpqZW+v6BAwe0fft2SVJERES5Y5fOnHzyyerfv7/j+xM3\n8Keddppjcp0yZYpee+21Gv8yuKMNb/aVlZXlOP7dvn37SifIssoek/z555/LvXfkyBGtXbtWkvzy\n/ubu/pk5UZs2bdSrV68qlznzzDMdrw8dOuRqyG7RtWtXx2v7/x9vio6O1s033yxJysvL09SpU50u\nN3HiRMfrqv7ffP311zLH95a65WvIkCFO+8nLyytXgysxMTHlvs/NzXW5jjf69VUdZc878nXAa9u2\nbZXbt6ZNmyolJcXx/fXXX19le126dHG8rs65OmXrd9d5dyfym1CwceNGx2t3BQKp/MbTmd9//93x\n+tRTT600AZZ1/vnnO11fOj5JvPbaawoPD1dRUZHuuece1a9fX7169dJ9992nDz74QHv27KmyfXe0\nUV3u6OvXX3913JN7y5YtioiIcHyFh4crPDxcYWFhji+bzaakpCTH+idOaps2bXIExFNOOUUnnXSS\nW2p1F3f/zJyoQ4cOLtsr+/nl5OS4XN6VrKwsTZo0SVdffbU6duyo5ORkRUREOJ5jYLPZNHjwYMfy\nR48erXOfteHqhMP9+/dr7ty5kqQWLVqoX79+XhtbZWJjYx2vCwoKXC5v/6vYLi4uzi/69VUdGRkZ\njtcn7n3wtvbt27tcpuxelHbt2lV72er8Hpddvuzn4k5+EwqOHDnieN28eXO3tVs2tTmTnp7ueF3d\nyadsWnO2cbz22mu1ePFiXXzxxbLZbCoqKtLSpUv16quvasSIEWrRooV69uypDz/80DH5eaKN6qpr\nX2X/30nHd4XZv0pLS1VaWlruL6oTnfjL4KmfBXfxxM9MWSf+heWMzWZzvK7r//9Zs2apXbt2GjVq\nlD799FP9/vvvOnbsWJUPX/HUnitXOnTooIsuukiStHr16gpntU+ePNkx7pEjRyoszPebuLIb8up8\nbtnZ2ZWu78t+fVVH2QASFRVVqzbcpaa/m66Wr+nvcdk9NCeGLnfx/W+Mh7naKNRmg1p2nbL/U8vq\n1auXvv/+e+3Zs0fTpk3Tfffdp+7duzuOLS9fvlw33HCDhg4dqqKiIo+1UV116au4uLhOfdd1UvM2\nT/3M+MLixYt11VVX1fgyw7KX83rbmDFjHK/L7i0wxuidd96RJIWHh1f7ElBPa9KkieN1dXb57ty5\n0/E6JSWl1sfQ3d2vr+ooOxEWFhbWqo1gUTYIlN1z405+Ewrq16/veO2uXePVUfZ8g/3791drnbLL\nudoT0bRpU11//fV65ZVXtHz5ch0+fFgTJ050/IX5xRdf6O9//7vH26iu2vRV9v/dJZdcUuNjsWU3\nHie2582fhery9M+MNz3yyCOOv6xbt26tl19+WUuXLtXevXuVm5urkpISx/+nsicC+9LQoUMdJ2FO\nnTrVsaFctGiRtm3bJknq37+/y71M/fr1K3d4pK5fn3/+udN+TjnlFMfr6pxMVvYa+FNPPdXl8pVx\nd7++qqPs4TlP7KGq6x813lS2/uoctqwNvwkFZX9ovv/+e6/1W/b47YYNGyo9o7mssieKVef4b1lJ\nSUkaOXKkZs6c6fi3KVOmeL0Nd/bVqVMnx1+/K1eurNbxxqqccsopjvY2btxY7YnXGU/8Ve7tnxlP\nSU9P1/LlyyUdD2IrV67Ugw8+qJ49e6pp06aKjY0tt6ctLS3NZZve2AsSFRXluBolIyPD8bNZ3RMM\nva1bt26O1999953L5efNm+d47eqcKG/266s6WrZs6Xi9b9++Gq1bnZMb67J98bay9Zf9XNzJb0KB\n/Xi2dHxi8VYwaNKkidq2bSvpeGIsu2FxZseOHeXuJlX2BLKa6Nmzp8LDwyXJcSa7L9pwR18NGzbU\nGWecIen4Rro6d5yrSmpqqmMjYoyp014QT+x69NXPjLsdOHDAcVijW7du5fbQOPPRRx+5bNNbu3pH\njx7t2F5MmjRJGRkZ+vTTTyVJzZo10xVXXOGxvmuq7Ama69atq/JyzrS0NM2ZM8fx/dChQ/2mX1/V\n0bp1a8ex+T179tToUKmrOylmZ2dr2bJltR6bt+3YscPxunPnzh7pw29CQbNmzcr90N10001V3krS\nneyXOUnSX//610p/kHJycnTjjTc6dje1b9++wqVj06dP1xNPPOH0tsplLVq0yLHbtl69em5vo7rc\n1dc999zjeP3QQw/pm2++cdn3wYMHNWbMGKef95133ul4/a9//Uvvvvuuy/acKburfvv27W47f8Gd\nPzO+Eh8f73i9Zs2aKs9mfv7557Vo0SKXbXrq8z5Ru3btdPHFF0uSFi5cqGeeecZxydxtt93mCLBV\n8dYlie3atVPPnj0d399xxx1OzzQvLi7W7bff7vhZadWqlS688MKafjQe69dXdURERDj+6CguLq7R\n9fnz5s2rcMlzWU899VS5bZ8vz5WpjnXr1kk6fq5c2T03blX7+x6534YNG0xcXJzjjk3x8fFm7Nix\nZvny5ebYsWMmPz/f7Ny50/z444/mb3/7m+nRo4fL2xzbn5JYlUOHDpmGDRs61klMTDR/+9vfzJYt\nW0xBQYE5dOiQ+eSTT8wpp5xS7m5Vzm5ZO3nyZCPJJCQkmJEjR5qZM2ea3bt3m/z8fJORkWF+/fVX\n8+yzz5a7beywYcPc3kZ1uauvE29zbLPZzLBhw8ysWbPM3r17TWFhocnOzjZbt241U6dONcOHD3fc\nlc/ZQ60KCwtNt27dyn3eQ4YMMXPmzDF//PGHKSwsNAcOHDCrV682EydONIMGDXJ6m2NjTLmHaz35\n5JNm7969VT6cqLq3OXbXz4wxzp/EVpV///vfjuVfeukll8tXpuwtk8844wwzc+ZMs3//flNUVGT+\n+OMP8+WXX5rLLruswt37xowZU2mbNf28a+uTTz6pMK6wsDCzc+dOt/dVV2Vv+SvrrnizZs0y6enp\nJisryyxYsKDcA5wkmSlTplTaXtmHDY0cOdJr/bq7veoqe9vqiRMnVrls2TsaSjJJSUnmlVdecdwu\nPDMz0yxZssRcc801FX5+Hn744QpPYKzp72bZn39XT4P9y1/+4lj2ww8/rHLZsncL7dmzp8tx1JZf\nhQJjjPn444/L3Yvf1Vd1npJYHV999ZWJjIysdr8333yz03bsk2x1v+Li4io8ec4dbVSXO/vau3ev\nOfnkk2vUXmX/D40xZvv27eXuD+7qq+xTEsu65557qlyvtk9JdNfPjDG+CwX/+Mc/qjV2m81m7rjj\nDsf3VYWCmn7etVVYWFju+ROSTP/+/d3StieMHTu22j8r119/vdPbYdtVNxS4u19PtFcdtX1KYp8+\nfaocX1hYmLn77rsr/HvZJ3H6SygI2ackDhs2THPnzlWLFi282m///v01e/ZsNW7cuMrlwsPDNXbs\nWMelTycaPny4XnjhhSrvomjXokULzZs3T506dXJ7G9Xlzr6aNm2qFStW6Nprr63WCWfNmjXT22+/\nXelusDZt2ujnn3+u81PRHn/8cY/8PLnrZ8aX7r//fpeX7sXFxWnSpEnlDu9VxVOf94kiIyN16623\nlvu322+/3eP91tbf//53Pfrooy4PbYwcOVLvvvuu207adHe/vqjj/PPPd1xNMnfu3CrvoVFW3759\nde+991b6/gsvvKBnnnmm3M3A/NWXX37peD18+HDPdeSxuFFH+fn55u233zaDBw82LVu2NLGxsSY+\nPt6cfPLJ5oILLjDjx483y5cvd5pCa7OnwC4rK8v885//NH379jVNmjQxkZGRJjk52Zx++unmgQce\nMBs3bqxWO7m5ueb99983V199tWnZsqWJiooysbGxpmXLlmbgwIHm7bffdvkADHe0UV3u7mvdunXm\n4YcfNj169DCNGzc2kZGRpl69eqZz587mlltuMZ9++qnLFF3Wt99+a0aNGmU6depkkpKSTFRUlGne\nvLk566yzzO23325mzZplcnNzK13/8OHD5oknnjDdu3c3ycnJJjw8vM57Cuzc8TPjqz0FdrNnzzaD\nBg1yjD81NdWcccYZZty4cWbbtm3GGGPmzp1brT0FxtTs866LX375xdHuSSedVKOfKV/55ZdfzN13\n3206d+5sEhMTTVxcnGnXrp257bbbzJIlS6rVRk32FLizX0+258q4ceMcNTt7cJGdswciTZ482fTo\n0cPUq1fPJCcnmz59+pjZs2c71vnyyy9N69atTUpKihk+fLjZvn274z1/2FOQlZXleIKwJw8dGOOH\nhw8AoLpef/11x0b1scce8/Vw4EFpaWmOQ8vXX399pcvV9CmJgeCtt95y1DRt2jSP9uV3hw8AoLrs\nT9O02Wx+cwdDeEbz5s11ww03SJI+/fTTGt+FM5D997//lXT88sxrr73Wo30RCgAEpBkzZjjurHfh\nhRc67h2B4PXkk08qIiJCBQUFevXVV309HK+YP3++Vq1aJel/9XsSoQBAQMnOztasWbP05z//2fFv\nd999tw9HBG85+eSTdccdd0iSXn/99XIPJwtWzz33nKTjj1keMWKEx/sjFAAICHfffbdsNpsSExM1\nePBgHT58WNLxjeVVV13l49HBW5555hk1bNhQGRkZev755309HI/65ptvtGDBAknSq6++Wq2bctWV\nZ/dDAIAHxcbGatKkSX7xiGR4R0pKSsicT3D55Zd7/Smy/CYBCDgpKSkaNGiQfvrpJ51zzjm+Hg4Q\nNGzG2zEEAAD4JfYUAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmEAgAAYOHmRaiT4uJipaWlSZLq\n1avHTWSAAFFaWqrMzExJUosWLTx+T30EBn4KUCdpaWk8iAYIcNu3b1ebNm18PQz4AUIBPC4lJUXD\nhg3zWf/h4eGKjIx0fF9UVKSSkpIaLxOoqD+06//kk09C4sFBcA9CAeqkXr16jteR4ZFqldTS8f2u\nY7tVVFKkqKgodejQwRfDcyguLlZhYaHj+6ioKMfu0qreCxbUH7r1R0VFSar4+yn973e07O8xQlvw\n/OTDJ8qeQ9AqqaW+GvGl4/sr3h+orUe3SZJOOukkr4/tRDk5OcrIyHB8b98QZmRkKCYmRpKUnJys\n+Ph4XwzP46g/tOs/8fdT+t/vKOcCwY5QgJBh39jbJ4ayE4QU3BOCRP2hXj9QHcRDhJT4+HglJydX\n+PdQmRCoP7TrB1whFAAAAEmEAoSYE48r22VkZCgnJ8f7A/Iy6g/t+gFXCAUIGSdOCMnJyeV2JQf7\nxED9oV0/UB2caIiQ4GxCKHsM+cSTz4Lt+DL1h3b9QHWxpwBBz9WEcOLJZ8H2FyP1h3b9QE0QChDU\nXE0IdsE6MVA/9RtjfD0MBBBCAYJWdScEu2CbGKif+p2dVAlUhVCAoFTTCcEuWCYG6qd+AgFqg1CA\noFPbCcEu0CcG6qf+svXbbDbfDQYBh1CAoFLXCcEuUCcG6qf+E+sHaoJQgKDhrgnBLtAmBuqnfnfW\nj9BEKEBQ8NQGMVAmBuqnfgIB3IFQgIDn6Q2iv08M1E/9BAK4C6EAAc1bG0R/nRion/oJBHAnQgEC\nlrc3iP42MVA/9RMI4G6EAgQkX20Q/WVioH7qJxDAEwgFCDi+3iD6emKgfuonEMBTCAUIKP6yQfTV\nxED91O8P9SN4EQoQMPxtg+jtiYH6qd+f6kdwIhQgIPjrBtFbEwP1U78/1o/gQyiA3/P3DaKnJwbq\np35/rh/BhVAAvxYoG0RPTQzUT/2BUD+CB6EAfivQNojunhion/oDqX4EB0IB/FKgbhDdNTFQP/UH\nYv0IfIQCeJwxpkYbxkDfINZ1YqB+6g/k+hHYCAXwiupuGINlg1jbiYH6qT8Y6kfgIhTAa1xtGINt\ng1jTiYH6qT+Y6kdgIhTA42w2m+N1ZRvGYN0gVndioH7qD8b6EXgIBfCKqjaMwb5BdDUxUD/1B3P9\nCCwRvh4AQoN9I2ff+JXdCIbCBpH6qV8K3foROAgF8JqqNoxS8G8QqZ/6pdCtH4GBwwfwqhN3pdqF\nygaR+qk/lOuH/yMUAAAASYQCeNmJJ1XZ1fbOb4GG+qk/lOuH/yMUwGucnWXtznvF+zvqp/5Qrh+B\ngRMN4RWuLrs68eSrYDu+Sv3UH8r1I3CwpwBeUdUG0dPPo/c1VxMC9VN/MNePwEIogMcZYxyvKzvL\nOlg3jNW9MQ31U38w1o/AQyiA17i67CrYNow1vVMd9VN/MNWPwEQogFdU9zrsYNkw1vbWtdRP/cFQ\nPwIXoQAeZ7PZanTiVKBvGOt6L3vqp/5Arh+BjVAAvxSoG0Z3PdyG+qk/EOtH4CMUwG8F2obR3U+7\no37qD6T6ERwIBfBrgbJh9NTjb6mf+gOhfgQPQgH8nr9vGD01IdhRP/X7c/0ILoQCBAR/3TB6ekKw\no37q98f6EXwIBQgY/rZh9NaEYEf91O9P9SM48eyDShQXF+vnn3/W/Pnz9fvvv2vfvn06duyYkpKS\n1LlzZ/Xr10/9+vVTeHi4W/rLz8/XJ598ojlz5mjPnj2Kjo5W+/btdd111+miiy7yWhv+zr7RtW+M\n7f/19r3ivT0h2FE/9Uu+rx/Bi1DgREZGhlq2bKmsrCyn73///fd67bXXdPrpp2vatGk69dRT69Tf\npk2bdPXVV2vTpk3l/n3BggV66623dO211+q9995TbGysR9sIFL7eMPpqQrCjfuqXCAbwDEKBE8XF\nxSoqKlK/fv103nnnqW3btmrevLkkadeuXZo2bZq+/vprrVu3Tn369NG6det00kkn1aqvgwcP6vLL\nL1daWpri4+M1btw49e7dWzk5Ofrggw80depUTZ8+XSUlJZo5c6bH2gg0vtow+npCsKN+6pcIBnA/\nQoETKSkpOnLkiOLi4py+f9NNN+m1117TPffco8OHD2vChAl6+eWXa9XX+PHjlZaWpvDwcM2dO1cX\nXHCB473LL79cLVq00IsvvqhPP/1Us2bN0qBBgzzSRiDy9obRXyYEO+qnfolgAPfiREMnwsPDKw0E\ndnfeeadSU1MlScuWLatVP4cOHdLkyZMlSTfeeGO5ydzumWeeUZMmTSRJ//d//+eRNgKZt06+8rcJ\nwY76qZ+TD+FOhIJaCgsLU0pKiiQpKiqqVm188cUXKikpkSSNGDHC6TKRkZG67rrrJElLly7V/v37\n3d5GoPP0htFfJwQ76qd+ggHchVBQS9u3b9eOHTskSeedd16t2vjpp58kHd8z0bNnz0qXO//88yVJ\nxpgKeyXc0UYw8NSG0d8nBDvqp36CAdyBUFBDxhh9/fXXGjBggEpLS9WmTRs99NBDtWrLfqVA8+bN\nqzxc0b59e8frjRs3ur2NYOHuDWOgTAh21E/9BAPUFScauvDTTz/pzjvvlCTl5uZq7969ys3NVXx8\nvG655RZNmDBBSUlJtWrbvhvf1ZULTZs2dbw+cOCA29twJT09vdL3ym40/YG7Tr4KtAnBjvqpX6pY\nP1BdhAIXMjMztXbt2gr/3q5dO3Xv3l316tWrddvZ2dmS5PKkxrLv29dxZxuu2E+oDBR1nRgCdUKw\no37ql/5XtzHGh6NBoCEUuHDeeedpzZo1ko7fMXD79u365JNP9MUXX+jOO+/U9OnT9dlnn9Vqb0Fx\ncbEkubwrYtn37eu4s41gVNuJIdAnBDvqp36JPQWoOUKBC4mJieratavj+549e+qGG27Q+++/r5tv\nvlkLFizQvffeq/fee6/GbSckJOjYsWPKz8+vcrm8vLxy67i7DVeOHj1a6XsZGRlq27ZtjdrzlppO\nDMEyIdhRP/UDNcWJhrU0YsQIDRgwQJI0derUKo+7V8Z+SeORI0eqXK7spFz2RCJ3teFKSkpKpV81\nbcvbqnvyVbBNCHbUT/02m83Xw0AAIRTUwSWXXCJJKikp0a+//lrj9e1XBOzatavK4347d+50vO7Q\noYPb2wh2riaGYJ0Q7Kg/tOsHaoJQUAdlb1pUmzRuPyyRk5NT4UFGZa1YsaLCOu5sIxRUNjGEyoRA\n/aFdP1BdhII6+OWXXxyvW7duXeP1Bw4c6Hg9a9asSpezv9e8efMKE7o72ggVziaGUJoQqD+06weq\ng1BQS7/99ps++OADSVK3bt0cT1GsiTPPPFPdunWTJL366qs6duxYhWXmz5/vuGvhbbfdVmGPhDva\nCCUnTgx2oTIhUH9o1w+4Qihw4oUXXtA//vEPbdmyRaWlpeXey8rK0qRJk3ThhRcqLy9PYWFheuml\nl5y2M3ToUHXt2tXpQ4rsJkyYoLCwMO3fv1+DBg3Srl27HO99//33jmcWNG3aVGPHjvVYGwAAcEmi\nE6tWrdLMmTP1l7/8ReHh4WrUqJFSU1OVk5Oj3bt3O4JCfHy83nrrLV188cVO29m0aZM2b95c5T0M\nLrroIv3zn//U/fffr8WLF+vkk09WmzZtlJOT47hbYUpKimbPnq3ExESPtREqTjyGbFfbO98FGuoP\n7foBV9hT4MT48eM1fvx4nXPOOY6/wDds2KCdO3fKGKP27dvroYce0saNG3XDDTfUub97771X8+fP\nV58+fWSM0datW7V//37Fx8frpptu0rp163TmmWd6vI1g5+ykslC6Vzz1h3b9QHWwp8CJrl27qmvX\nrnr22WdljNHhw4d16NAhRUdHq1GjRtX+a/vzzz9Xfn6+y7sNSlKfPn3Up08f5eTkaO/evYqOjlbT\npk0VGRlZ7XG7o41g5eos89reEjdQUH9o1w9UF6HABZvNpoYNG6phw4Y1XrdTp041Xic+Pr7O9xFw\nRxvBxNWE4K6H6Pgr6g/t+oGa4PABglp1r0MP1sfOUj/180Ak1AShAEGrpjemCbaJgfqp39lJlUBV\nCAUISrW9U12wTAzUT/0EAtQGoQBBp663rg30iYH6qb9s/aF8szLUHKEAQcVd97IP1ImB+qn/xPqB\nmiAUIGi4++E2gTYxUD/1u7N+hCZCAYKCpzaIgTIxUD/1EwjgDoQCBDxPbxD9fWKgfuonEMBdCAUI\naN7aIPrrxED91E8ggDsRChCwvL1B9LeJgfqpn0AAdyMUICD5aoPoLxMD9VM/gQCeQChAwPH1BtHX\nEwP1Uz+BAJ5CKEBA8ZcNoq8mBuqnfn+oH8GLUICA4W8bRG9PDNRP/f5UP4IToQABwV83iN6aGKif\n+v2xfgSfCF8PAKHh4MGDtV63uLhYBQUFju+jo6OVnZ2t7OxsdwzNbexjPHjwoKKjoxUR4Z5fL+qn\n/kCoH8GBUACPCwsLU2xsbK3WLS4uliTHBjYyMtJtG1t3i42NVVFRkeN7d4yV+qlfqlv9YWHsEEb1\n+edvF4JKZGSkmjdv7uthACEpMjLS10NAACEUwOMSEhJ05ZVX+noYQEhKSEjw9RAQQNivBAAAJBEK\nAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAA\ngIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICF\nUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAA\nAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAA\nJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQR\nCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoA\nAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACA\nhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQ\nAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAA\nACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAk\nEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEK\nAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAA\ngIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICF\nUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAA\nAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAA\nJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQR\nCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoA\nAICFUAAAACQRCgAAgIVQAAAAJBEKAACAhVAAAAAkEQoAAICFUFBN2dnZysrK8kpfWVlZKigo8Hkb\nAIDQQihwYu/evXrvvfd06623qlOnTkpMTFRiYqLq1aunhIQE9e3bV1OnTpUxxm19rl+/XjfccIOS\nkpJUr149xcTEqFWrVnrkkUd07Ngxr7UBAAhdhAInHn/8cd1yyy169913tXnzZmVnZys+Pl4JCQnK\nycnR/PnzdeONN+riiy9WZmZmnfubNm2azj77bH344YfKzMxUQkKCIiMjtXv3bv39739Xly5dtHXr\nVo+3AQAIbYSCKpx77rl69913tX//fsfhg40bN+rqq6+WJC1cuFDDhw+vUx8rVqzQrbfeqsLCQnXp\n0kUrV65UVlaWcnNz9cEHHyghIUFpaWm64oorlJ+f77E2AAAgFDhx2mmn6bvvvtNPP/2km2++WU2a\nNHG817lzZ02fPl1DhgyRJH399ddatGhRrfsaO3asCgsLlZqaqm+//VZnnXWWJCkiIkI33nij3n//\nfUnSli1b9O9//9tjbQAAQChwYuzYserbt2+l79tsNj3++OOO7+fNm1erfjZu3KjFixdLku666y41\nbty4wjJDhw5V165dJUn//e9/PdIGAAASoaDW2rZt63h94MCBWrUxa9Ysx2v7IQln7O9t27ZNGzZs\ncHsbAABIhIJa27Vrl+N1SkpKrdpYvXq1JCk2NlZdunSpdLnu3btXWMedbQAAIBEKau2DDz5wvL7g\nggtq1cbmzZslSa1atVJYWOX/K9q0aVNhHXe2AQCAJEX4egCBaOPGjXr99dclSe3atdMVV1xRq3aO\nHj0qSWrQoEGVy5V9376OO9twJT09vdL3MjIyatQWAMB/EQpqKCsrS9dcc40KCwtls9n0n//8R5GR\nkbVqKycnR5IUExNT5XKxsbGO19nZ2W5vw5XU1NQaLQ8ACEwcPqiB4uJiDR8+XJs2bZIkjR8/Xpde\nemmt27PZbJLk8s6IZd8/8RCBO9oAAEBiT0G1lZaW6pZbbtHcuXMlSbfddpueeeaZOrWZkJCgo0eP\nKi8vr8rlcnNzy63j7jZcqepwQ0ZGRrkrMQAAgYtQUE133XWXpk6dKkkaNmyY3nrrrTq32ahRI+3e\nvdvlJY1l32/UqJHb23CltldXAAACC/uRq2Hs2LF68803JUkDBw7UlClTFB4eXud2O3fuLElKS0ur\n8omG27Ztq7COO9sAAEAiFLj05JNP6uWXX5YkXXbZZZo+fXqtTyw8UY8ePSRJRUVFWrFiRaXL/fjj\nj47XZe834K42AACQCAVV+r//+z89++yzkqQ+ffro888/V3R0tNvaHzJkiONEwWnTpjldprS0VB9/\n/LEkqVu3bmrVqpXb2wAAQCIUVOr111/Xo48+Kknq1auXvvzyy3KX9VXH4cOHdeDAAf3xxx9O32/W\nrJmGDRsmSZo0aZLWrVtXYZl//OMfjrsnPvjggx5pAwAAiRMNnfrggw90zz33SJI6duyoyZMnKzs7\nu9Lr+6Ojo52ejHf++edr8+bNSkpKqvQmPy+++KK+++47HTlyRH379tWLL76oPn36KCcnRx988IHj\n0MVFF12kG264wWNtAABAKHBi0qRJjuv6N2/erI4dO1a5fN++ffXdd9/Vqq/WrVtr9uzZuvrqq7Vv\n3z6NGjWqwjK9e/fWp59+6jhM4Ik2AAAgFDiRmprq9BHEVS3vTMOGDZWRkaGkpKQq1+/Zs6d+/fVX\nvfXWW/rqq6+Ulpam6OhodejQQcOHD9d1113n8oZD7mgDABDaCAVOfPrpp25pZ8mSJdVeNiUlRY88\n8ogeeeSRWvfnjjYAAKGLPx0BAIAkQgEAALAQCgAAgCRCAQAAsBAKAACAJEIBAACwEAoAAIAkQgEA\nALAQCgAAgCRCAQAAsBAKAACAJEIBAACwEAoAAIAknpIIbyhMl7a+5etRAKGpMN3XI0AAYU8BAACQ\nRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQxH0K4A1RKVK70b4eBRCaov4lab+vR4EA\nwZ4CAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyE\nAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIA\nACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAg\niVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQ\nAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAA\nACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAs\nhAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQC\nAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAA\nIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJ\nUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAA\nAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAA\nLIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyE\nAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIA\nACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAg\niVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQ\nAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAA\nACyEAgAAIIlQAAAALIQCAAAgiVAAAAAshAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAs\nhAIAACCJUAAAACyEAgAAIIlQAAAALIQCAAAgiVAAAAAsEb4eQCA4dOiQfv75Z+3fv1+S1KlTJ51/\n/vlu76eoqEjLli3Tnj17FB0drfbt26tLly5ebwMAEJoIBZWYMmWKZs+ereXLl2vnzp3l3hs5cqRb\nQ0FpaaleeeUVPf/88zpy5Ei590455RS9/vrr6tOnj8fbAACENg4fVOLNN9/Uxx9/rJ07d6pBgwa6\n/PLLPdbX7bffrr/85S86cuSImjdvruHDh2vAgAGKjY3Vxo0b1bdvX82YMcPjbQAAQhuhoBI33XST\nPvzwQ23btk2HDh3S559/7pF+3nvvPb3zzjuSpNGjR2vHjh366KOP9OWXX2rr1q06/fTTVVpaqhEj\nRigtLc1jbQAAQCioxJgxY3Tdddepbdu2HuujpKRETzzxhCSpS5cuev311xUR8b8jOk2bNtUnn3yi\n8PBw5eXl6bnnnvNIGwAASIQCn1q0aJHjL/d77rmn3GRu17FjR11xxRWSpE8++URFRUVubwMAAIlQ\n4FNff/2143VV5yzY38vIyNCyZcvc3gYAABKhwKd+/fVXSVJqaqpatmxZ6XJnnnlmhXXc2QYAABKX\nJPrUtm3bJKnKyVySWrRoUWEdd7bhSnp6eqXvZWRk1KgtAID/IhT4UGZmpiSpXr16VS6XlJRUYR13\ntuFKampqjZYHAAQmDh/4UG5uriQpOjq6yuViYmIqrOPONgAAkNhT4FP2idzV1QCFhYUV1nFnG64c\nPXq00vcyMjI8etkmAMB7CAU+lJiYqEOHDiknJ6fK5bKzs8ut4+42XElJSanR8gCAwMThAx+ynxy4\nZ8+eKpcr+/6JJxS6ow0AACRCgU+deuqpkqQDBw5UuYt+06ZNjtennHKK29sAAEAiFPiU/amFxhgt\nXLiw0uXmz58vSYqKilKvXr3c3gYAABKhwKf69+/vOL7/3//+1+kyR48edTzdsF+/fhUuPXRHGwAA\nSIQCj5oxY4YmTpyoDz74wOn78fHx+stf/iJJmjdvnt56661y7xcVFWnUqFHKzs6WzWZzPPjI3W0A\nACBx9UGlNm/erCVLlji+Ly4uLvfexIkTHd9HRkbq5ptvrtDG+PHjtXnzZiUlJemmm25y2s8jjzyi\nb775RkuXLtWYMWM0d+5c9enTRzk5Ofroo4+0fv16SdLjjz+us88+22NtAABAKKjEkiVLdPvttzt9\n74cfftAPP/zg+D4+Pt5pKKiOmJgYzZkzR6NHj9aMGTP0+eef6/PPP3e8HxUVpWeeeUaPPPKIR9sA\nAIBQUIlOnTpp5MiR1Vq27N0Cy7r22mu1f/9+xcXFVbl+SkqKpk+fro0bN2rOnDlKS0tTdHS0OnTo\noCFDhqhhw4Yux+CONgAAoY1QUInzzz9f559/fp3aePbZZ2u0/CmnnFLnywXd0QYAIDRxoiEAAJBE\nKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgA\nAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQJEX4egAIfnlb0rT05MG+HgYQkvLS0nw9BAQQ\n9hQAAABJhAIAAGAhFAAAAEmEAgAAYCEUAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmEAgAAYCEU\nAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmEAgAAYCEUAAAASYQCAABgIRQAAABJhAIAAGAhFAAA\nAEmEAgAAYCEUAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmEAgAAYCEUAAAASYQCAABgIRQAAABJ\nhAIAAGAhFAAAAEmEAgAAYCEUAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmEAgAAYCEUAAAASYQC\nAABgIRQAAABJhAIAAGAhFAAAAEmEAgAAYCEUAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmEAgAA\nYCEUAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmEAgAAYCEUAAAASYQCAABgIRQAAABJhAIAAGAh\nFAAAAEmEAgAAYCEUAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmEAgAAYCEUAAAASYQCAABgIRQA\nAABJhAIAAGAhFAAAAEmEAgAAYCEUAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmEAgAAYCEUAAAA\nSYQCAABgIRQAAABJhAIAAGAhFAAAAEmEAgAAYCEUAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmE\nAgAAYCEUAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmEAgAAYCEUAAAASYQCAABgIRQAAABJhAIA\nAGAhFAAAAEmEAgAAYCEUAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmEAgAAYInw9QAQ/I4mJunN\nATf7ehhASDo6ZbV0NMvXw0CAYE8BAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQ\nRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQo\nAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAA\nABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAW\nQgEAAJBEKAAAABZCAQAAkOTjUFBaWurL7gEAQBkR3uhk//79WrhwoRYuXKiVK1fqjz/+0KFDh1RU\nVKTU1FQ1bNhQ7du3V+/evdW7d2+deeaZCgtjJwYAAN7ksVBQUlKiL7/8Um+88Ya+/fZbGWOcLnf4\n8GEdPnxYmzZt0qxZsyRJLVu21OjRozVq1Cg1btzYU0MEAABluD0UGGP0/vvv64knnlBaWlq59xo1\naqRTTz1Vqampql+/vqKjo3XkyBEdPXpUu3fv1m+//abS0lLt3r1b48eP19NPP62bb75Zzz33HOEA\nAAAPc2soWLJkie6//36tXr1akhQdHa0BAwZo6NChOvfcc9W2bdsq18/MzNTy5cv1/fffa9q0aUpL\nS9PEiRP18ccfa9y4cXrwwQcVHR3tziEDAACLWw/cP/DAA1q9erXat2+vN998UwcOHNCMGTP0pz/9\nyWUgkKR69erp0ksv1Ysvvqhdu3Zp/vz5GjZsmLKzs/XYY49p/fr17hwuAAAow62hoF27dpo8ebI2\nbdqkMWPGKDk5udZt2Ww2XXTRRfr444+1bt06XXXVVZx8CACAB7n18MFHH33kzuYcTjvtNM2cOdMj\nbQMAgOP40xsAAEgiFAAAAAuhAAAASPLSHQ3tcnNz9dlnn2nhwoXaunWrMjMzVVJSUuU606ZN0ymn\nnOKlEQIAELq8FgoWL16sP/3pT9qzZ0+N1svNzfXQiAAAQFleCQU7duxQ//79HRN8y5Ytdfrppysl\nJcXlZYYNGjTwxhABAAh5XgkFr776qnJzc5WUlKSpU6fqyiuv9Ea3AACgBrxyouGaNWskSX/9618J\nBAAA+CmvhILIyEhJUvfu3b3RHQAAqAWvhAL71QNHjx71RncAAKAWvBIKxowZI5vNpo8//tgb3QEA\ngFrw2p6CCRMmaOrUqXr11Ve90SUAAKghr92n4MEHH1Tjxo11//33a8qUKRo2bJg6dOiguLi4Ktfr\n0aOHEhMTvTRKAABCl1fvaJicnKzWrVtrxYoVWrFiRbXWWbFihc4++2wPjwwAAHgtFDz88MN66aWX\nvNUdAACoIa+Egnnz5jkCwRlnnKGRI0dW+46G7dq188YQAQAIeV4JBe+//74k6ZJLLtHcuXMVEeHV\noxYAAKAavHL1wa5duyRJ9913H4EAAAA/5ZVQ0LBhQ0lSSkqKN7oDAAC14JVQcOmll0qS1q9f743u\nAABALXglFIwYMUKdO3fWhAkTdOzYMW90CQAAasgroSA+Pl5z585VcnKyLrzwQv3www/e6BYAANSA\nV876e+CBB7RlyxY1aNBA3333nS644AI1bdpUHTt2dHlHw1deeYXLEgEA8AKvhIIlS5Zo1apV5f5t\n37592rdvn8t1//rXv3poVAAAoCyvhIIePXooOTm5VuvWq1fPvYOphZ9++klz5szRnj17FB0drfbt\n2+vaa69V69at/a4fb40VABB8vBIKXn/9dW9043YHDhzQiBEj9O2331Z4b9y4cXr44Yf13HPPubwr\nozf68dZYAQDBizsJVSI7O1v9+vXT2rVrFRYWphEjRqh3797KycnRtGnT9NNPP+mFF15Qbm6uXnnl\nFZ/2462xAgCCm9/+2Th//nz9+OOPPuv/ueee09q1ayVJU6dO1eTJk3XLLbforrvu0pIlS3TjjTdK\nkv71r3/VaZzu6MdbYwUABDe/DQUPPvigrrzyShUUFHi976ysLL366quSpEGDBum6664r935YWJj+\n85//KCkpSdLxSdlX/XhrrACA4OeVwwfvvfee9u/fX+3l9+7dq/Xr16u0tFTbt29X586dPTi6imbP\nnq28vDxJ0qhRo5wuk5iYqGHDhuntt9/W999/r/T09Brfxtkd/XhrrACA4OeVUPDvf/+7wiWJ1REf\nH6+mTZt6YERVW7x4sSTJZrPpwgsvrHS5Pn366O2331ZRUZF++uknXXnllV7vx1tjBQAEP68cPoiL\ni1N8fHyVX5GRkeXWOeecc/Tll186dnt704YNGyRJTZs2rbL/Tp06OV5v3LjRJ/14a6wAgODnlT0F\n9r9mXdm9e7feffddPffcczpy5Ijat2/v4ZE5t3fvXklSs2bNqlyu7Pt79uzxST/eGGt6enql72Vk\nZNSoLQCA//KrSxJbtmypJ598Ug0aNNBdd92lu+66S1988YXXx5GVlSXp+OGLqpR9Pzs72yf9eGOs\nqampNVoeABCY/PLqg9tvv13x8fGaM2eOT/4SLSwslCRFRFSdmaKiohyva3OVhDv68dZYAQDBzy9D\nQWRkpNq0aaOSkpJa7ZavK/tDmlxNnvaz/iXXf6l7qh9vjPXo0aOVfm3fvr1GbQEA/JdfHT4oy76H\nwBfPPkhKStL+/ftd7qU4duyY43VtxumOfrwxVi5fBIDQ4Jd7Cvbu3au9e/eqdevWatGihdf7P/nk\nkyUdP/GxKrt27XK8rs3jnd3Rj7fGCgAIfn4XCrKzs3XHHXdIkl566SXZbDavj+H000+XdHxvxY4d\nOypdbvXq1Y7XXbp08Uk/3horACD4eeXwwWuvveby3ABjjHbv3q1FixZp//79Ouecc7Ry5UqtXLnS\n6fL33HOPy8vwaqt///564YUXJElfffWV7rrrLqfLzZ07V5JUv3599ejRwyf9eGusAIDg55VQ8O67\n79b4joYrVqzQihUrKn3/mmuu8VgoOO+883TyySdr27ZtevXVVzVq1ChFR0eXW2bt2rWOxxTfcMMN\nCg8P90k/3horACD4eeXwQePGjdWsWTO3fpW9xM7dwsLC9Pzzz0uSfv/9d916663Kzc11vL9161YN\nGzZMpaWlSkxM1GOPPea0nfvuu0/XXHONbr75Zo/1466xAgDglT0Fc+bM8UY3bjVs2DAtXbpUr7zy\nij788EN988036t69u3JycrR06VIVFxcrIiJCn3zyiZo0aeK0jW+++UabN2+u8vbD7ujHHW0AAOB3\nJxr6k3/+85+aOHGiWrZsqaNHj+rrr7/WkiVLVFxcrHPPPVdLly5Vv379/KIfb40VABC8/PY+Bf5i\n5MiRuu222/Trr78qLS1N0dHRat++vVq2bOly3VdffVWZmZnVOtRRl37c2QYAIHQRCqrBZrOpS5cu\nNb6U77LLLvNKP+5uAwAQmtx6+CA/P9+dzZVjjOGe/QAAeJBbQ0H//v31+OOPV/mo3Zoyxuizzz7T\nGWecofXr17utXQAAUJ5bQ0FOTo6ef/55tW7dWg899FCdJvGMjAy99dZb6tatm6666iqtX79eYWGc\nFwkAgKe4dZadPn26hg8frszMTE2YMEGnn366unbtqqefflrffPNNlQ/tKS4u1po1a/Tmm2/qmmuu\nUZMmTTRmzBj98ssvateunWbOnKlu3bq5c7gAAKAMt55o2KpVK3300Ue6//779dhjj2nBggVau3at\n1q5dK+n4SXAnnXSS6tevr9TUVEVHRzsewbtv374K5yQ0a9ZMY8eO1V133aXIyEh3DhUAAJzAI1cf\n9OzZU/Pnz9emTZv0xhtvaOrUqTp69KiMMdq3b5/27dtX+YAiItS3b1/dcccdGjhwILfkBQDASzx6\nSWLnzp316quv6pVXXtHatWu1ePFirVq1SgcPHtShQ4dUWFioBg0aqGHDhmrfvr169+6t8847TwkJ\nCZ4cFgAAcMIr9ykICwvTmWeeqTPPPNMb3QEAgFrgdH4AACCJUAAAACyEAgAAIIlQAAAALIQCAAAg\niVAAAAAshAIAACCJUAAAACxeuXkRQlt9U6QHi3b5ehhASPrBFOmYrweBgMGeAgAAIIlQAAAALF4J\nBdddd51OO+00DRw4UEeOHKnROhs2bPDw6AAAgOSlULB161Zt2LBBs2fPVo8ePbR58+Zqr5OXl+eF\nEQIAAK8ePrDZbNq2bZt69uyp+fPne7NrAADggldDwRtvvKGTTjpJGRkZ6tevnyZOnOjN7gEAQBW8\nGgrOOuss/fzzz+ratauKiop0++236+GHH1Zpaak3hwEAAJzw+tUHzZs31w8//KCBAwdKkl566SVd\nffXVys3N9fZQAABAGT65JDE+Pl6ff/65HnjgAUnS559/rgsuuED79u3zxXAAAIB8eJ+CsLAw/eMf\n/9Cbb76piIgIrV69Wt27d9eaNWt8NSQAAEKaz29eNGbMGM2dO1dJSUnau3evLrjgAs2aNcvXwwIA\nIOT4PBRI0iWXXKKlS5eqTZs2ysnJ0dChQ7Vp0yZfDwsAgJDiF6FAkjp37qzly5erV69eKi0t5cRD\nAAC8zCuh4IILLtCVV16p5OTkKpdr2LCh5s+frxtuuMEbwwIAAGV45dHJ//znP6u9bHR0tKZOnapJ\nkyY5vgcAAJ7nlVBQGzExMb4eAgAAIcVvzikAAAC+RSgAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAA\nkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJII\nBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUA\nAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADA\nQigAAACSCAUAAMBCKAAAAJKkCF8PAMEvPTZGb559uq+HAYSk9C9ipHRfjwKBgj0FAABAEqEAAABY\nCAUAAEASoQAAAFgIBQAAQBKhAAAAWAgFAABAEqEAAABYCAUAAEASoQAAAFgIBQAAQBKhAAAAWAgF\nAABAEqEAAABYCAUAAEASoQAAAFgIBQAAQBKhAAAAWAgFAABAEqEAAABYCAUAAEASoQAAAFgIBQAA\nQBKhAAAAWAgFAABAEqEAAABYCAUAAEASoQAAAFgIBQAAQBKhAAAAWAgFAABAEqEAAABYCAUAAEAS\noQAAAFgIBQAAQBKhAAAAWAgFAABAEqEAAABYCAUAAEASoQAAAFgIBQAAQBKhAAAAWAgFAABAEqEA\nAABYCAUAAEASoQAAAFgIBQAAQBKhAAAAWAgFAABAEqEAAABYCAUAAEASoQAAAFgIBQAAQBKhAAAA\nWAgFAABAEqEAAABYCAUAAEASoQAAAFgIBQAAQBKhAAAAWAgFAABAEqEAAABYCAUAAEASoQAAAFgI\nBQAAQBKhAAAAWAgFAABAEqEAAABYCAUAAEASoQAAAFgIBQAAQBKhAAAAWAgFAABAEqEAAABYCAUA\nAEASoQAAAFgIBQAAQBKhAAAAWAgFAABAEqEAAABYCAUAAEASoQAAAFgIBQAAQBKhAAAAWAgFAABA\nEqEAAABYCAUAAEASoQAAAFgIBQAAQBKhAAAAWAgFAABAEqEAAABYCAUAAEASoQAAAFgIBQAAQBKh\nAAAAWAgFAABAEqEAAABYCAUAAEASoQAAAFgIBQAAQBKhAAAAWAgFAABAEqEAAABYCAUAAEASoQAA\nAFgIBQAAQBKhAAAAWAgFAABAEqEAAABYCAUAAEASoQAAAFgifD0ABL+TkuP1xm19fT0MICQtfjle\n6ft8PQoECvYUAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmEAgAAYCEUAAAASYQCAABgIRQAAABJ\nhAIAAGAhFAAAAEmEAgAAYCEUAAAASYQCAABgIRQAAABJhAIAAGAhFAAAAEmEgkrl5ubqhx9+0Msv\nv6zhw4erTZs2iomJUUxMjP785z+7vb/9+/frscce0xlnnKH69euradOm6t27t9544w0VFRV5rQ0A\nQOiK8PUA/NXll1+uH374wel77p5g58+fr2uvvVZHjx4t9+/79+/X4sWLNXHiRM2dO1eNGjXyaBsA\ngNDGnoJKxMXF6bzzztMDDzygDz/8UBs3bvRIP7///ruGDh2qo0ePqlmzZpoxY4YOHTqknTt36tln\nn1VERIRWr16tQYMGqaSkxGNtAADAnoJKfPPNN+W+z8/P90g/Dz30kDIzMxUXF6cFCxaoffv2kqQG\nDRpo/PjxatCggf785z9r+fLlmjx5skaNGuWRNgAAYE+BD+3cuVOzZs2SJI0ePdoxmZc1ZswYtWvX\nTpL0r3/9yyNtAAAgEQp86osvvnC8vu6665wuY7PZNGzYMEnSr7/+qm3btrm9DQAAJEKBT/3888+S\npKioKHXr1q3S5Xr16uV4vWLFCre3AQCARCjwqd9++02S1LJlS0VGRla6XNu2bSus4842AACQONHQ\npw4dOiRJaty4cZXLlX3fvo4723AlPT290vcyMjJq1BYAwH8RCnwoOztbkhQTE1PlcrGxsRXWcWcb\nrqSmptZoeQBAYOLwgR+w2WzVft8Y47E2AAChjT0FPpSQkKD09HTl5eVVuVzZ9xMTE93ehisn3iWx\nrIyMjHLnKwAAAlfI7Cl48803Hc8uOPHLfrmet9WvX1+S62P8f/zxR4V13NmGKykpKZV+JScn16gt\nAID/CplQUFxcrIKCAqdfhYWFPhlTp06dJEm7d++u8vbDO3bsqLCOO9sAAEAKoVDw5z//WXl5eU6/\npk+f7pMxnXXWWZKO30L5l19+qXS5pUuXVljHnW0AACCFUCgIDw+v9PBBVdf3e9KQIUMcr6sKJjNm\nzJAkdezYUR07dnR7GwAASCEUCvxRu3btdNlll0mS3njjDaWlpVVYZsqUKY4nNN59990eaQMAAIlQ\n4FFnnHGGYmJi1KRJk0qXmTBhgmJjY5WZmamLL75Y33//vYqKinTs2DG9/vrrjicann766RozZozH\n2gAAgFBQiXfffbfcIYayZ9mf+F6DBg2ctmE/kbGqxy536dJFH330keLj47V161Zdcsklio6OVnJy\nsu6++24VFBSoQ4cOmj17dqWHOdzRBgAAhIJKOLtawa6kpKTcv1c16VfHoEGDtG7dOt1xxx1q2bKl\nwsLCFBsbq9NPP11/+9vftGbNGrVo0cLjbQAAQhs3L6rErbfeqhtvvLFay1Z2N8F169aptLTU5d0G\npeMPLHrjjTckSaWlpQoLq3lec0cbAIDQRSioRHh4uMLDw+vURlRUVK3Wc8dkTiAAANQUMwcAAJBE\nKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgA\nAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAA\nFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZC\nAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEA\nAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQ\nRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQo\nAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAA\nABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAW\nQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIB\nAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAA\nkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBE\nKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgA\nAAAWQgEAAJBEKAAAABZCAQAAkEQoAAAAFkIBAACQRCgAUIWCokKl52T5ehgAvIRQAMCpwuIi3fbu\ni7rsnw9q/qZVvh4OAC8gFABwav5vq/XbgV3Kys/Rm4u+kDHG10MC4GGEAgBOzVi1QNkFeTqUlaHt\nh/dpze7ffT0kAB5GKABQwdY/9uiXtK1Kz8lSfnGh8goLNGPVQl8PC4CHEQoAVDBz1SIVFhcppzBP\nYbYwpedmaf5vq3U0J9PXQwPgQYQCAOXkFuRrzvqlSs/NUkR4hBomJiszP0cFxUWa9csPvh4eAA8i\nFAAo5+sNy5VTkKdjuTlKiU1QUmyCbLLpWF62Pl2zWKWlpb4eIgAPIRQAcDDGaMaqhcrMz1GJKVVy\nXILCw8JULzZe6blZ2pdxWEu3/+rrYQLwEEIBAIdf927X7wfTlJ6TpcSYOEWER0iSUuISj59jUJDP\nCYdAECMUAHCYuXqR8osKlVdUoJS4BMe/x0RGKTYyWhm5Wfphy3rtzzjsw1EC8BRCAQBJUmZejuZt\nXKH03CxFhUcqLiqm3PspcYnKzM9VUUmxPluzxEejBOBJhAIAkqQv1/6o/MICZeblKCU+UTabrdz7\nibFxCreFKSMvW5//skRFxcU+GikATyEUAJAxRjNXL9KxvBwZGSXFxldYJswWpqS44yccHs3O1KLf\nf/H+QAF4FKEAgFbs3KTdRw8qPTdL9WLiFR4W7nS5lLhEFZcUK7sglxMOgSBEKACgGasWKq+wQAXF\nhUqJS6x0uaiISMVHxSo9N0srd/2mnYf3e3GUADyNUACEuENZ6Vq4+Rel52YpJiJKMZFRVS6fHJeg\n7IJ8FRUXa+bqRV4aJQBvIBQAIe6LX35QUUmRMvNylOzkBMMTJcbEKSLs+PMQZq/7SflFBV4aKQBP\nIxQAIayktESfrl6sjNxs2Ww2JcVUPMHwRDabTclxicrIyz5+GeOGFV4YKQBvIBQAIeyHLev0R2a6\n0nOzlRQbr7Cw6m0SkuMSVFJaqqx8TjgEggmhAAhhM1YtVE5hnopKiqo8wfBEkeERSog+fsLhxv07\ntWn/Ts8NEoDXEAqAELU3/ZCWbt+g9JwsxUXFKNrFCYYnSolPVG5hvgqKijRzFSccAsGAUACEqJmr\nF6m4pERZBXlKrsFeArv4qBhFhkcqPTdLX29Yruz8XA+MEoA3EQqAEFRYXKRZa39Qem6WwsPClBgT\nV+M2bDabUuISdCwvW7kF+ZqzfqkHRgrAmwgFQAj6ftMqpedkKSM3S8mxCQpzcRliZZJiE2SMUWZ+\njmasWihjjJtHCsCbCAVACJq5eqGyC/JUXFqq5DKPSK6piPBwJcbGKz0nSzsO79ea3b+7cZQAvI1Q\nAISYrX/s0S9pW5Wek6X46BhFRUTWqb2UuATlFxcqr7CAyxOBAEcoAELMzFWLVFhcpJzCvBpdhliZ\n2MhoRUdEKT03S/N/W62jOZluGCUAXyAUACHEfkJgem6WIqx7DdTV8RMOE5WZn6OC4iJ98csSN4wU\ngC8QCoAQ8vWG5copyNOx3BylxCa4fM5BddWLjZdNNh3Ly9Zna5aotLTULe0C8K4IXw8AgPfMWLVQ\nmfk5KjElOpqbpYy87HLvx0fFqklSapVhIacgX39kHlWJKT/xG2OUnpulfRmHtXT7rzqv3ekeqQGA\n5xAKgBBRXFKs3w+mKcwWpkaJKRXeLzVGh7MzlByfoNjI6ErbOZJz7Pghg1jn5yOUlJZqw76dhAIg\nABEKgBARER6hBy4Zpm82/Cyj8vcTOJqTqb0Zh3U4O0OmtOp7DRhjFBMZrYaJKerYpEWF9xskJGng\n6b3cOnYA3kEoAELIn3pepj/1vKzCv3+yYr6en/tBjdpKjkvQByOfcNfQAPgBTjQEAACSCAUAAMBC\nKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSuE8BAB8wxmjrLz8qJ/Oo6p/USo1atFN0bLyvhwWE\nPEIBAK8ryMvR4X07VFxYoNzMdO3Zsl4pjZupSauOSmpwktse1ASgZggFALzOlJZIkgoL8mTychQe\nEaXDe3co/cAexcQnqlHL9mrU4mRFRsX4eKRAaCEUAPCpsPAImdIS5WVnKiw8XEWF+crLyVTa5l9U\n/6RWatyqgxJTGrL3APACQgEAnwqPiFBUTJxKiotVXJivwvxcFebnKiIySgd3b9HhfTsUl5isxi07\nqEGzNoqIjPL1kIGgRSgA4BfCIyIUHpGgKFOq4sJCFRXmqyj7mMLDI1RUWKCczHTt+m21GjZro8at\nOii+XqqvhwwEHUIBgHLyigorPFq5rJLSUp0UHasmEZHKOLy/Vn0U5uVU+p7NFqbI6BhFREWrtKRY\nRYX5KszLUWFejiKiorV/x286uHuLEpIbqHGrDmpwUiuFhbMpA9yB3yQA5fyRdbTK9y9r1EL9GrdS\nTFSUNi3/zmPjsNlsCo+IVHhEpExpqYoKC1RcmK+ignyFR0SquLBQWemHtGvjSjVsfrIat+yg2IR6\nHhsPEAoIBQDUrlEzRdjC1aFxS5fLdqt/khIjo1QvIko5x6oOEK4YYxQWHe5yOVtYmKJiYhUZHaOS\n4iIVFxYoPzdLtrAwRUZGa++2Ddq/Y5OS6jdR41YdldK4ucLCuDcbUFOEAnhcZkmmPj38qa+HgSqY\nOKMbhpyjA0cyXC7b7IhNcQVhCjc2FUZWfpihWsKkIuVIRbk1Wy9KUoRNtuISFRdkSwXZMuE2ZeQf\n1c4Dm2SiwlXSJF7FjeOlaoSOYJZZkunrISCAEAoAyGazqfVJDdX6pIYul41af0jhh/NkKyyRifTh\nX+NhNpmo8ONjKDGyFZXKllcsW5hNpqhUtl2ZikjLUmlqjIqbxKs0OVriskagSoQCAIHNZpMibDIR\nYTKlVjgoLDkeWiLCFP5HrsKO5MnERFh7D+KkyNDeewBUhlAAIHiE2WSiw2WiwqRiI1vxCXsP8osV\nsTtTRe1SVNIoztejBfwOZ+IACD42mxQZJhNjHV4oNbIVlMiWVyyVGIX/UfklkUAoY08BgOBjP4xQ\nXCpJMpFhx8NBmE0mNkJFrZJ8PEDAPxEKAAQHY/53yKDE/O9ExIjjYaCkfoxKmiSoNCmKEw6BShAK\nAAS2E/cKhIfJxIbLhNtkoiNU0iTu+KWJUZxcCLhCKIDH1Quvp6saXOXrYcBNfo3+RkfDd6vQlquI\nosg6tWULC1dkdEyNn4BojFFJcZGKCvJVUlyksLBwRcTEKyIqWmFhYUpu2ExNWnVQcqNmIf90xSfC\nn/D1EBBACAUAasRms8kWFqaw8AiVlpbWqa3SokLZbDZFRsdUb/nSUhUX5qu4sEClpaUKj4hUTHyi\nIiKiFBkdo0Yt2qlRy/aKiUuo07iAUEUoAFAjDZufrOyMw255hHFuVoaMqTpYHN8rcPyxysVFhccf\nmBQVfXyvQHi46qU2VuNWHZTapIXCwjhEANQFoQBAjTRqcbIaNG2l0pKSWreRn5ul9T/OrXIZU1qq\noqICFRcUqLS0ROEREYqJS1B4ZJQiIqPUsFlbNW7VQXGJybUeB4DyCAUAaiwsPKJOjysOL8x3+u/G\nGOtxyQUqKSqUJEVERismOlFh4eGKr5d6/HHJTVsrPKJu5zMAqIhQAMDnjDEqLixQUWGBSkuKFRYe\noaiYOEVERSs8PEL1m7ZWk1YdlJDcwNdDBYIaoQCATxUXFaqo4Pieg4jIKEXHxis8PEKxCfXUuFUH\nNWzWVhFR0T4eJRAaCAUAfMtIUTGxioiMUVh4mFKbtFTjVh1UL7VxyF9OCHgboQCA19nPB4iJS1BY\nWISiY+PVqGV7NWrRTlExsT4eHRC6CAUAvC4qJk6tOp+l3Mx01T+pFTcZAvwEoQCATzRte4qvhwDg\nBDw6GQAASCIUAAAAC6EAAABI4pwCl4qLi3Xw4EHt379fxhg1adJELVq08Fh/+/fv1549exQdHa3W\nrVurXr16PmkDABB62FPgxMqVKzVu3Dide+65io2NVfPmzXXOOeeoe/fuatmypZo3b64nnnhC2dnZ\nbuvziy++ULdu3dS0aVN1795dZ5xxhho0aKDBgwdry5YtXmsDABC6CAVOvPbaa3rxxRe1bNkyFRcX\nKyUlRWeddZbOPvtsxcbGau/evXruuefUtWtX7d69u879PfvssxoyZIjWrFmjmJgYnXXWWerUqZOK\nioo0a9Ysde3aVT/++KPH2wAAhDZCQSWio6N16623atmyZTp8+LBWrlypFStW6NChQ3rsscckSdu2\nbdOQIUNUUoenxX311Vd68sknJUn9+/dXWlqaVq5cqU2bNmnFihVq3ry5cnNzNWTIEKWnp3usDQAA\nCAVOXHrppdq4caPeeecd9ejRQ2Fh//uY4uPj9be//U133nmnJGnNmjX66quvatWPMUYPP/ywJKlV\nq1aaPn26GjT43wNfzj77bH3yySeSpMOHD+vFF1/0SBsAAEiEAqf+9Kc/qW3btlUuYw8FkrR48eJa\n9fPzzz9rw4YNkqR7771X8fHxFZY599xz1adPH0nS+++/r9LSUre3AQCARCiotcaNGzteHz16tFZt\nzJkzx/F6wIABlS535ZVXSpIOHDigVatWub0NAAAkQkGt/fbbb47XZQNCTaxbt06SlJiYqA4dOlS6\n3DnnnFNhHXe2AQCARCiotbfeesvx+pJLLqlVG/bLBFu1alXlcmXf//33393eBgAAEjcvqpWFCxdq\nypQpko7/BX7xxRfXqp2MjAxJUkpKSpXLlX3/2LFjbm/DlaquWLD3DwAIfISCGtqzZ4+uu+46GWMU\nExOjiRMn1rqtnJwcSVJMTEyVy5V9376OO9twJTU1tUbLAwACE4cPauDYsWMaMGCADh48KEl6/fXX\ndfrpp9e6vcjISEnHb6VclbLv29dxZxsAAEghtKdg7969ld7qt0GDBjrttNOqXD8vL08DBw7U2rVr\nJUkvvPCCbrvttjqNKSEhQYcPH1Zubm6Vy5V9PyEhwe1tuFLV1RUZGRkuL98EAASGkAkFn332me65\n5x6n7w0ePFiff/55pesWFhbqqquu0pIlSyRJ48aN06OPPlrnMTVr1kw7d+7Uvn37qlxu7969jtdN\nmzZ1exuuuDpfAQAQHEImFDRr1ky9e/d2+l5VewmKi4t13XXX6euvv5Yk3XfffXr++efdMqbOnTvr\nxx9/1N69e5WdnV3pX/BlrxY45ZRT3N4GAABSCIWCoUOHaujQoTVap7S0VDfffLM+++wzSdLo0aP1\nyiuvuG1M559/viZOnKjS0lL9+OOPuvzyy50uZ79jYlhYmM4991y3twEAgMSJhpUyxmj06NGaNm2a\nJOnmm2/Wm2++6dY+Bg4cqKioKEnS5MmTnS6Tl5enjz/+WJJ04YUXqmHDhm5vAwAAiVBQqQceeECT\nJk2SJF133XV65513ZLPZatTGzz//rIULF+qHH35w+n5qaqpGjx4tSfrkk08chyjKevjhh3X48GFJ\nx89l8EQbAABIIXT4oCaee+45/etf/5IkdezYUSNHjqzyoUcpKSk644wzKvz7iBEjtHnzZiUlJVV6\nk5+nn35aX3/9tbZu3arBgwfr3nvvVZ8+fZSTk6MPPvhAs2fPlnR8T8Vll13msTYAACAUOPHdd985\nXm/evFmXXnpplcv37du33Do1kZqaqnnz5mnYsGFauXKlJkyYoAkTJpRbZtSoUfrPf/7j0TYAACAU\nONG1a9caLe9sL4Ekde/eXU2aNHF5X4A2bdpo6dKlmjNnjr766iulpaUpOjpaHTp00PDhw9WtWzeX\nY3BHGwCA0EYocMJdVxi8//771V42IiJCgwcP1uDBg2vdnzvaAACELk40BAAAkggFAADAQigAAACS\nCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggF\nAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAA\nwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBC\nKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigA\nAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAA\nkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJII\nBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUA\nAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADA\nQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIo\nAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAA\nAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACS\nCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggF\nAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAA\nwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBC\nKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigA\nAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAA\nkggFAADAQigAAACSCAUAAMBCKAAAAJIIBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJII\nBQAAwEIoAAAAkggFAADAQigAAACSCAUAAMBCKAAAAJKkCF8PAMEvOztbc+bM8fUwgJCUnZ3t6yEg\ngBAK4HEFBQXauXOnIiKq9+NWXFysoqIix/eRkZHVXtdf1KUG6qd+d9Zfti3AlcD6TUNAMsaoqKio\nRhun4uJiFRQUOL6Pjo4OmInB2dilmm2cqZ/63VV/aWlprdZDaAqM3zIENJvNpsaNG9d4vZycHGVk\nZDi+T0hIUHx8vBtH5n72MdsnguTk5FqPmfqpP9DqR+DjREP4rfj4eCUnJzu+z8jIUE5Oju8G5MKJ\nG/G6TAgS9VN/YNWP4EAogF8LlA2juycEO+qn/kCoH8GDUAC/5+8bRk9NCHbUT/3+XD+CC6EAAcFf\nN4yenhDsqJ/6/bF+BB9CAQKGv20YvTUh2FE/9ftT/QhOhAIEFH/ZMHp7QrCjfur3h/oRvAgFCDi+\n3jD6akKwo37qJxjAUwgFCEi+2jD6ekKwo37qJxjAEwgFCFje3jD6y4RgR/3UTzCAuxEKENC8tWH0\ntwnBjvqpn2AAdyIUIOB5esPorxOCHfVTP8EA7kIoQFDw1IbR3ycEO+qnfoIB3IEHIiFo2DfW9o24\n/b/ueiCNv04IdtRP/VLF+oGaYE8Bgoq7/mIKtAnBjvqp/8T6gZogFCDo1HViCNQJwY76qb9s/cYY\n3w0GAYfDB1UoLi7WypUr9fvvv2vfvn06duyYkpKS1LlzZ/Xp00dJSUlu7W/Pnj2aO3eu9uzZo+jo\naLVv315XXHFFjTZI7mgjGNR2V3KgTwh21E/9EnsKUHOEAidycnJ0/fXXa9GiRcrMzHS6TExMjEaP\nHq3nn3++zhuN7OxsPfjgg5o4cWKFVJ+YmKgJEyZo9OjRHm8j2NR0YgiWCcGO+qkfqClCgRN5eXn6\n8ssvFRcXp969e6tt27Zq3ry5JGnXrl368ssvlZ6erldffVU///yzFi1apKioqFr1VVRUpCFDhuj7\n77+XJPXp00e9e/dWTk6OZsyYoZ07d2rMmDE6duyYHnroIY+1EayqOzEE24RgR/3Ub7PZfD0MBBBC\ngRMJCQn65ptvdNFFFykyMrLC+9nZ2RozZoymTZumZcuW6fXXX9cDDzxQq77++c9/Oibzv//97+Um\n7WeeeUZDhgzRvHnz9Oijj6p///467bTTPNJGMHM1MQTrhGBH/aFdP1ATnGjoRExMjC677DKngUA6\nHhr++9//KjY2VpI0Z86cWvWTn5+vF198UZJ04YUXVvgrPjY2Vu+//75iY2NVWlqqp59+2iNthILK\nTj4LlQmB+kO7fqC6CAW1lJCQ4DikUNubhHzzzTdKT0+XJN15551Ol2ncuLGGDh0q6Xj4OLEvd7QR\nKpxNDKE0IVB/aNcPVAehoJays7O1Z88eSar17nj7Ln9Juvjiiytd7pJLLpF0/FyHn376ye1thJIT\nJwa7UJkQqD+06wdcIRTUQkZGhkaNGqW8vDxFRUVp7NixtWpnw4YNko7/Jd+wYcNKlzv11FMrrOPO\nNgAAkDjR0KXNmzfr7bffliTl5uZqz549WrBggbKzs9WqVStNmTJFHTt2rFXbu3btkiTHYYjKtGjR\nosI67mzDFfvhCWcC7TroE48h21X3OvZAR/2hXT/gCqHAhR07dujll1+u8O+9evXSpEmT1KlTp1q3\nnZWVJen4fQSqkpCQUGEdd7bhSmpqao2W91fOTiqT3HevfH9H/aFdP1AdhAIXOnXqpJdeeknS8TP9\nt2/frtmzZ+unn35S586d9de//lVPPfVUrdrOz8+XJJf3OIiOjq6wjjvbCAWuzjIP9omB+kO7fqC6\nCAUutG7dusI5A3l5ebrzzjv17rvv6q9//auio6P16KOP1rjt2NhYZWZmqqCgoMrlyk7i9ssg3dmG\nK0ePHq30vYyMDLVt27ZG7XmbqwnB3U/X8zfUH9r1AzURMqHgxx9/1Geffeb0vdNOO0233HJLtduK\njY3VpEmTtGbNGq1du1YvvviiHnjggXJ/jVdHvXr1dPDgwUpvpWxXdnd/vXr13N6GKykpKTVa3p9U\n9zr0YJ0YqJ/6eSASaiJkQsGaNWucnhsgSYMHD65RKJCksLAwDRs2TGvXrtWxY8e0bt06nXPOOTVq\no02bNtqyZYvS0tKqXK7s+23atHF7G8GqpjemCbaJgfqpP9BOBIbvhUwoOP/88x3nBpyoffv2tWqz\nQYMGjtfZ2dk1Xv+0007TvHnzdPjwYe3bt09NmzZ1utzatWsdr8teWuiuNoJRbe9UFywTA/VTP4EA\ntREyoaBr167q2rWrW9vcuXOn43WTJk1qvP5ll12mf/zjH5Kkb7/9VjfffLPT5ebNmyfp+BUGvXr1\ncnsbwaaut64N9ImB+qm/bP08EAk1wc2LaiknJ0dTpkyRJJ100km1ujTx4osvdoSJf//73yotLa2w\nzM6dO/Xll19Kkq666qoK5y24o41g4q572Vd2r3x/R/3Uf2L9QE0QCpyYM2eO1q1bV+n727dvV//+\n/R3H6ceNG+c0jf/973/X2LFj9cQTTzhtJzIyUk8++aQkadWqVXrkkUfKnRSUnp6u6667TkVFRYqK\ninLajjvaCBbufrhNoE0M1E/97qwfoSlkDh/UxOTJkzVz5ky1bNlSHTt2VJMmTZSamqqcnBxt2LBB\ny5cvd/xFPmrUKN19991O23nnnXe0efNmJSUl6dlnn3W6zB133KH58+drxowZmjBhgubMmaPevXsr\nJydHs2fPdtxN8M0339TJJ5/ssTYCnac2iIGyK5n6qZ9AAHcgFDgxaNAg7d27VytWrNDu3budLnPm\nmWfq0Ucf1bBhw+rUl81m07Rp09SpUye9/PLL2rRpkzZt2uR4v3nz5nrttdc0ePBgj7YRyDy9QfT3\niYH6qZ9AAHchFDgxYsQIjRgxQunp6Vq7dq3279+vQ4cOKTo6Wo0bN1a3bt3UsmVLl+088sgjOnLk\niGJiYqpcLjIyUs8++6weeeQRLViwQGlpaYqOjlaHDh103nnnKSzM9VEed7QRiLy1QfTXiYH6qZ9A\nAHciFFQhJSVFffr0qfX6t956a42WT0hI0MCBA2vdn7vaCBTe3iD628RA/dRPIIC7Beefjwh6vtog\n+svJZ9RP/QQCeAKhAAHH1xtEX08M1E/9BAJ4CqEAAcVfNoi+mhion/r9oX4EL0IBAoa/bRC9PTFQ\nP/X7U/0IToQCBAR/3SB6a2Kgfur3x/oRfAgF8Hv+vkH09MRA/dTvz/UjuBAK4NcCZYPoqYmB+qk/\nEOpH8CAUwG8F2gbR3RMD9VN/INWP4EAogF8K1A2iuyYG6qf+QKwfgY9QAI8zxtRowxjoG8S6TgzU\nT/2BXD8CG6EAXlHdDWOwbBBrOzFQP/UHQ/0IXIQCeI2rDWOwbRBrOjFQP/UHU/0ITIQCeJzNZnO8\nrmzDGKwbxOpODNRP/cFYPwIPoQBeUdWGMdg3iK4mBuqn/mCuH4GFRyfDKyp77OyJr4N1g0j91C+F\nbv0IHIQCeE1VG0Yp+DeI1E/9UujWj8DA4QN41Ym7Uu1CZYNI/dQfyvXD/xEKAACAJEIBvOzEk6rs\nanvnt0BD/dQfyvXD/xEK4DXOzrJ2573i/R31U38o14/AwImG8ApXl12dePJVsB1fpX7qD+X6ETjY\nUwCvqGqD6Onn0fuaqwmB+qk/mOtHYCEUwOOMMY7XlZ1lHawbxuremIb6qT8Y60fgIRTAa1xddhVs\nG8aa3qmO+qk/mOpHYCIUwCuqex12sGwYa3vrWuqn/mCoH4GLUACPs9lsNTpxKtA3jHW9lz31U38g\n14/ARiiAXwrUDaO7Hm5D/dQfiPUj8BEK4LcCbcPo7qfdUT/1B1L9CA6EAvi1QNkweurxt9RP/YFQ\nP4IHoQB+z983jJ6aEOyon/r9uX4EF0IBAoK/bhg9PSHYUT/1+2P9CD6EAgQMf9swemtCsKN+6ven\n+hGcCAUIKP6yYfT2hGBH/dTvD/UjeBEKEHB8vWH01YRgR/3UTzCApxAKEJB8tWH09YRgR/3UTzCA\nJxAKELC8vWH0lwnBjvqpn2AAdyMUIKB5a8PobxOCHfVTP8EA7kQoQMDz9IbRXycEO+qnfoIB3IVQ\ngKDgqQ2jv08IdtRP/QQDuEOErwcAuIt9Y23fiNv/W9uNeKBMCHbUT/1SxfqBmmBPAYKKu/5iCrQJ\nwY76qf/E+oGaIBQg6NR1YgjUCcGO+qm/bP3GGN8NBgGHUICgVNuJIdAnBDvqp/6y9QPVxTkFCFo1\nPcYcLBOCHfVTP1BT7ClAUKvuX4zBNiHYUT/122w2Xw8DAYRQgKDnamII1gnBjvpDu36gJjh8gJBQ\n1eVaoTAhUH9o1w9UF6EAIcPVddzBPiFQf2jXD1QHoQB1UlpaWq3lCgsLPTyS6omMjFRCQoIyMzPL\n/Xu9evUUGRnpN+P0FOoP7forU93fYwQ/QgHqpOzGddex3bri/YHlvpek4uJirVmzxutjq0xkZKSi\no6PL/duhQ4dUVFTkoxF5F/WHVv3FxcWSKv5+2v9NOv57XL9+fa+PDf7HZrizBepgx44datu2ra+H\nAaAOtm/frjZt2vh6GPADhALUSXFxsdLS0iQd3wUbFnb8gpaMjAxHWNi+fTs3UglA/D8MfFX9Pywt\nLXXs6WvRooUiIthxDA4foI4iIiJc/oWRnJyslJQUL40InsD/w8Dn7P8hhwxwIu5TAAAAJBEKAACA\nhVAAAAAkEQoAAICFUAAAACQRCgAAgIX7FAAAAEnsKQAAABZCAQAAkEQoAAAAFkIBAACQRCgAAAAW\nQgEAAJBEKAAAABZCAQAAkEQoAAAAlghfDwD+Z+3atfr888+1c+dO2Ww2tW3bVkOHDtWpp55a57Zf\nfPFFrVy50uVynTt31rPPPuvTsQYCY4x+/PFHzZ07V3v27FF0dLQ6dOiga6+9Vq1atfKrfrw11kBT\nXFysefPmacGCBdq/f7/q1aunLl266JprrlHDhg3r3P4111xTreVuvPFGDRkyxKdjhR8wgCUzM9Pc\ncMMNRpLTr9tuu83k5ubWqY8rr7yy0vbLfp133nk+H6u/27dvn+nbt6/T+sPDw8348eNNaWmpX/Tj\nrbEGml9//dV06dLF6ecSHx9v3nrrrTr3UZ3fN0nmhRde8PlY4XvsKYAkqaSkRNdcc43mzZsnSRoy\nZIgGDBigkpISff7555o7d67eeecdpaena+bMmbLZbHXq75RTTtHTTz9d6fsNGjTwm7H6o8zMTF1+\n+eVav369wsPDNWLECF144YXKycnRlClTtGzZMj333HPKy8vThAkTfNqPt8YaaHbu3Km+ffvq4MGD\nio+P1+23365u3brp4MGDmjhxojZv3qzRo0fLZrNp1KhRde5vxIgRGjhwYKXvd+nSxW/GCh/ydSqB\nf3jjjTeq/Ivhsccec7w/ZcqUWvdj31PQu3dvvx+rP/vLX/7iqHH69Onl3ispKTHXX3+94/1ly5b5\ntB9vjTXQ2H8XYmJizIoVK8q9l52dbc4991wjycTFxZl9+/bVuh/7Z/vSSy/5/Vjhe4QCmNLSUtOi\nRQsjyXTt2tXpbtzi4mLTvn17I8l06tSp1n3VNRR4c6z+Kj093cTExBhJ5qqrrnK6TEZGhklMTDSS\nzIABA3zWj7fGGmjWrFnjmKwfe+wxp8ts2LDBsczYsWNr3VddQ4E3xwrf4+oDaNmyZUpLS5Mk3Xbb\nbU53t4eHh+uWW26RJP3222/69ddfvTlEh0Aaq6d8+eWXys/Pl6RKd9UmJSU5TjCbN2+eMjMzfdKP\nt8YaaKZPn+54Xdnncsopp+jcc8+VJM2YMcMr43ImkMaKuiMUQIsXL3a87tOnT6XLlX1v0aJFHhxR\n5QJprJ5i/wzCw8N1wQUXVLqc/TMoLCzU0qVLfdKPt8YaaOyfS8uWLdWmTZtKl7N/Ljt37tSuXbu8\nMbQKAmmsqDtONIQ2bNggSbLZbOrYsWOly3Xq1MnxeuPGjXXq89ChQ/rb3/6mrVu3qri4WA0aNNA5\n55yjK664QsnJyX41Vn9j/wyaN2+uhISESpc78TO4/PLLvd6Pt8YaaOw/k507d65yuRM/l7pcurly\n5Uo9/PDD2rt3r6Kjo9WsWTNdcMEF6tu3r8LDw/1qrPAdQgG0d+9eSVLDhg0VFRVV6XKpqamKjY1V\nXl6e9uzZU6c+N27cqPHjx1f494SEBD366KMaN26cwsIq7sjyxVj9jf0zaNasWZXLlX2/Np+BO/rx\n1lgDSV5eno4ePSrJu5/Lxx9/7PTf27RpozfeeMNpEPPVWOE7hAIoKytLkhQfH+9y2bi4OOXl5Sk7\nO7vW/Z188sm65JJLdPLJJ6tZs2bKysrSzz//rI8//ljZ2dkaP368Nm7cqClTplQ4Z8DbY/VH1f0M\nyr5fm8/AHf14a6yBxP6ZSN75XGw2m84991z16tVLrVq1Umpqqg4cOKB58+Zp3rx52rFjh/r3768P\nPvhAf/rTn3w6VvgeoQAqLCyUJEVEuP5xsP91XlBQUKu+/v3vfzs9LjlmzBg988wzGjRokFavXq1p\n06bpsssu08033+yzsfqr6n4GZfek1OYzcEc/3hprILF/JpJ3Ppdt27Y5/Z178MEHNX/+fA0ZMkRZ\nWVkaNWqUevfurebN/7+9ew+KqorjAP6FZWFEEHRRQBARFJVMsxpHslSYfABpGuWg+cg0ndTK8tmk\nTjqVNcrU+Aq1/lDRFE1zcFLJ55gPnLRRyiQfiagwIuADEZbdfv0h3FjZXXAfXFi/n5mduew595zf\nPTPL/vbcxwlVLVZSHy80JHh7ewOo3wf5wYMHAOr3S90caxcqhYSEICMjQ2k7JSWlVp2GjLWxqu8Y\nVB8/YNsYOKKfhoq1KakeE6BhxsXaZy4uLg6pqakAgPLycqxatcqkvKFjJfVxpsBFnTp1CosXLzZb\n1qFDByxZskT528/PDwBw+/Ztq20ajUZlWrBFixaOCfQRbdu2xfDhw5GWlobs7GwUFRVBp9M1yljV\n4ufnh8LCwjrH4M6dO8q2LWPgiH4aKtamxNfXF25ubhCRRjEuycnJeO+991BcXIxDhw6ZlDW2WMn5\nmBS4qPz8fPz4449my3r06GHyd2RkJICHj6O9ffu2xav/b9y4AYPBAADo2LGj44J9RM27CgoKCkyS\ngsYWqxoiIyNx8eLFOm/7qlluyxg4op+GirUp0Wq1CAsLQ25ubqMYF3d3d3Tq1AlZWVnIz89v1LGS\n8zEpcFHPP/+8yUNHanr0i7R79+7K9u+//47Y2Fiz+50+fVrZtvacdHvVnKZ89A6DxharGrp37469\ne/eiqKgIeXl5aNeundl69o6BI/ppqFibmu7duyM3Nxdnz56F0Wi0eEtg9bi4u7s7deXP6s+cuTt6\nGlus5GRqP1KR1Hf16tV6PaJ08uTJyqp2hYWFTounf//+AkC8vLxqrXTY2GJVw4EDB5QxSE1NtViv\nelXCNm3aiNFoVKWfhoq1qVm1apUyLocOHTJbx2g0SlhYWL1WDbVHSUmJaLVaASAJCQmNOlZyPiYF\nJCL/fxHrdDq5detWrfLr168rz6cfMmSI0+LYs2eP8g/I0nPwG0usajEYDBIeHi4ApGvXrlJRUVGr\nzqlTp8TNzU0AyIcffqhaPw0Va1Nz8+ZN8fb2FgCSmJhots6mTZuUz8Lq1audFsu0adOsJm6NKVZy\nPiYFJCIix48fF3d3dwEgsbGxUlRUpJQVFBQoq6B5eHjImTNnzLbx2WefSVJSkiQlJcn9+/drle/Y\nsUPS0tKkuLi4Vpler5c1a9aIj4+PABCtViunT592WqxN3caNG5V/wqNHjzaZUcnJyZGOHTsKAPHz\n85ObN2+abWPatGmSlJQk48ePd2o/jmjDFc2fP18ZlwULFojBYFDKjhw5Iq1atRIAEhUVJXq9vtb+\nRqNR+bzNmzfPbB9Lly6Vffv2md2/sLBQ3n33XSWGzp07m03aHBErNR1MCkjx9ddfKx98Hx8fefnl\nlyUuLk6aNWumvL927VqL+w8aNEipV1JSUqt8zpw5ypd1hw4dJCYmRoYOHSqxsbHSsmVLZd9mzZpJ\nenq6U2N1BVOnTlWOVafTyeDBg+XFF18UjUajJFaZmZkW94+MjFT2dWY/jmrD1VRWVsrgwYOVcQkN\nDZWEhAR57rnnlPdatmwp586ds7h/dT1LU/ZPPfWUAJDmzZtLly5dpH///jJ06FDp1auXeHp6KvtH\nRETIpUuXnBYrNR1MCsjEzp07pUuXLsoHvfrVrVs32bt3r9V960oKsrKy5K233pLg4OBa7QMP12If\nM2aMnD9/3umxuorVq1dLaGhorTHo06ePxZmWavVNCuztx5FtuBq9Xi+ffvqp8ku7+qXRaGTIkCFy\n5coVi/vWJylYt26dJCYmKtP/j75CQkJkwYIFcu/ePafGSk2Hm4gIiB6Rk5ODK1euAAAiIiLQqVOn\nOvf59ddfUVBQAAB49dVXodVqLdbNy8vDjRs3cPPmTXh6eiI4OBjR0dH1elKhI2J1JSKC7OxsXLt2\nDV5eXoiKirJ4lX9Ne/bsQWlpKby8vDBkyBCn9ePoNlxRZWUlzpw5g4KCAvj6+iI6OhqtW7e2uo+I\nKLcdBwQEWF01VK/XIzc3F/n5+SgpKYGvry8iIiIQHh7eILFS08GkgIiIiADwMcdERERUhUkBERER\nAWBSQERERFWYFBAREREAJgVERERUhUkBERERAWBSQERERFWYFBAREREAJgVERERUhUkBERERAWBS\nQERERFWYFBAREREAJgVERERUhUkBERERAWBSQERERFWYFBA9IdasWYOJEydi9+7daodil3379mHi\nxIlYuXKl2qEQuRw3ERG1gyAi57p06RKio6Oh0Whw+fJlBAUFqR2SzW7duoUOHTqgvLwc2dnZ6NKl\ni9ohEbkMzhQQPQFmz54NvV6PKVOmNOmEAAACAgIwbdo0GAwGzJw5U+1wiFwKZwqIXNyJEycQExMD\nT09P5OXloU2bNmqHZLeioiKEhISgoqIChw4dQr9+/dQOicglcKaAyMUtX74cAJCQkOASCQEA6HQ6\nvPLKKwD+Pz4ish+TAiIXVlhYiG3btgEAxowZo3I0jlV9PDt37sSNGzdUjobINXioHQDRk+b8+fNY\nunQpAKBfv34Wv6zLysowe/ZslJeXo23btli0aNFj95WWlga9Xg9vb28kJiZarJeZmYn09HRotVp8\n++23Fuv99NNP2LVrF/z8/JCSkmJSdvbsWSxbtgwAsHjxYrRu3RrFxcXYsmULLly4gObNm+OFF15A\nfHx8rXaPHTuGgwcPIj8/H23btsVrr71W5wWECQkJ8PHxQWlpKdatW4ePP/7Yan0iqgchogb3+uuv\nCwDRarVy8uRJs3UmTJggAMTd3V32799vUz8JCQkCQF566SWr9ZYsWSIAxMvLy2q9Tz75RABIYGBg\nrbKMjAwBIADkwoULsmXLFvH29lbeq34NGDBAysrKRETk7t27Eh8fX6uORqORlJSUOo8vNjZWAEhc\nXFyddYmobjx9QKSCtWvXIjw8HJWVlUhOTsbdu3dNytPT0/H9998DAObOnYu4uLjH7sNoNOLo0aMA\ngN69e9sf9GM4ePAgkpOT8cwzz+DLL7/E8uXLMXz4cADAL7/8gqlTp8JoNCI+Ph6HDx/GpEmTsGLF\nCixcuBDt27eH0WjEjBkzcODAAav9xMTEAACOHz8OvV7v9OMicnlqZyVET6rjx4+Lh4eHAJDk5GTl\n/cuXL4ufn58AkJiYGKmsrLSp/ezsbOWX98aNG63WdfRMga+vryxatKhWnblz5yozAe+//74EBQVJ\nTk6OSZ3CwkIJDg5WZhWs2bp1q9JnVlaW1bpEVDfOFBCppHfv3sp1Aps3b8Z3330Hg8GAkSNH4s6d\nO/D398cPP/wADw/bLv3Jz89XtgMCAhwSc3316NED8+fPr/X+nDlz4OHhAaPRiGXLlmHFihWIiooy\nqRMQEIDJkycDeDjjUFFRYbGf1q1bK9s1j5eIbMOkgEhFc+fOxYABAwAAH3zwAcaNG4esrCwADx9L\n3L59e5vbLiwsVLZ1Op19gT6mt99+2+z7/v7+6Ny5MwCgZcuWGDZsmNl6vXr1AgAYDAZcvHjRYj81\nj6vm8RKRbZgUEKnIzc0NGzZsQGBgIMrKyrBp0yYAwKRJk/DGG2/Y1XZpaamy3bx5c7vaelw9evSw\nWFY9a9GtWzdoNBqzdWrOANy+fdtiWzWP6969e48ZJRE9ikkBkcoCAwNNptqDgoLwzTff2N1uq1at\nlG1rX6zO0KJFC4tl1adD6lMHACorKy3Wq3lcNY+XiGzDpIBIZXfu3DG557+goAA///yz3e3W/LVd\nXFxsd3uNUVFRkbJd83iJyDZMCohUNmnSJPzzzz/QarXo2rUrAOCdd97B1atX7Wq35vUIBQUFdrVV\nzWAwOKQdR6l5XPZcf0FEDzEpIFLR2rVrkZ6eDgD44osvkJGRAV9fX5SUlGDUqFEwGo02tx0eHo7Q\n0FAAwG+//VavfSoqKqyearh27ZrN8TjDyZMnATw8dRAdHa1yNERNH5MCIpWcO3cO06dPBwAMHDgQ\nM2bMQGRkJFJTUwEAR48excKFC+3qo3r1wOo7Guqj+oFHj3rw4AEyMzPtisfRTpw4AQDo27cv3Nzc\nVI6GqOljUkCkgvLyciQnJ6OsrAyBgYFYv3698qU2atQojBs3DgDw+eef4/Dhwzb3k5CQAODhugT1\nvdhw+vTp+Pvvv03eKy0txdixY5Xb/u7fv2/XLIYj3L17F2fOnAHw/3ESkX2YFBCp4KOPPkJ2djbc\n3Nywfv16BAYGmpRXP9Tn33//xZtvvmlyQd3jSEpKgk6ng8FgUE5TWOPh4YH79+/j6aefRt++fTFm\nzBgMHToU7dq1w7Zt2/Dss88CeJgkDBo0CF999ZVNcTnCtm3boNfr4evri5EjR6oWB5ErYVJA1MB2\n7NihrEQ4a9YsDBw4sFYdHx8fbN68GZ6enrh+/TomTJhgU19eXl4YP348AGD9+vV11tdoNNi1axdC\nQkJw5MgRpKWlISMjA/fu3VPWIggLCwMA7N+/H/v377cpLkfYsGEDAGD06NHw8fFRLQ4iV+ImIqJ2\nEERPku3bt6O4uBgajQajR4+GVqu1WPfgwYO4dOkSAGDYsGE2Pa74ypUr6Ny5M/R6Pf7880+zF+Qt\nXboUs2bNgpeXF8rLy2E0GnHs2DHk5ubCx8cHvXv3RlBQEACgpKQEmZmZ8Pf3R2xsLDw9PQEAeXl5\n2Lt3LwBgxIgRFp9DsHv3bly/fh1hYWFmEyLg4S2U27dvBwAkJiYiODjYpPzixYuIioqCu7s7/vjj\njzqXWSai+mFSQPQEmDlzJlJSUjBixAhs2bKlVvmjSUFjN3bsWGzYsAFTpkzBypUr1Q6HyGXw9AHR\nE2DevHnQ6XTYunUrsrOz1Q7HLjk5Odi0aRP8/PzsvjuDiEzZtvwaETUp/v7+2LlzJ/766y+rqw42\nBWVlZUhNTUVUVFSDr/5I5OqYFBA9Ifr06YM+ffqoHYbdevbsiZ49e6odBpFL4ukDIiIiAsCZAiIC\nMGjQIPj7+5usTkhETx7efUBEREQAePqAiIiIqjApICIiIgBMCoiIiKgKkwIiIiICwKSAiIiIqjAp\nICIiIgBMCoiIiKgKkwIiIiICwKSAiIiIqjApICIiIgBMCoiIiKgKkwIiIiICwKSAiIiIqjApICIi\nIgBMCoiIiKgKkwIiIiICwKSAiIiIqjApICIiIgDAf30fAhjzj5SGAAAAAElFTkSuQmCC\n" } } ], "source": [ "sim = make_sim()\n", "f, ax = plt.subplots(1, 1, figsize=(10, 10))\n", "sim.plot(y=0, ax=ax)\n", "plt.show()" ], "id": "f7afb97b" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we write a function to post process some run results to get the\n", "transmission we are after." ], "id": "a18917b4-cca5-42b6-a46e-058b6807d674" }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def post_process_T(sim_data: td.SimulationData) -> float:\n", " \"\"\"Given some td.SimulationData from the run, return the transmission of \"p\" polarized light.\"\"\"\n", " amps = sim_data[\"diffraction\"].amps.sel(polarization=\"p\").values\n", " return anp.sum(abs(amps) ** 2)" ], "id": "fa387241" }, { "cell_type": "markdown", "metadata": {}, "source": [ "And finally, put everything together in a single function that relates\n", "the permittivities and thicknesses of each slab to the transmission,\n", "through a `td.Simulation` run." ], "id": "744b3e52-72f1-44be-a764-9c63de37abfd" }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def compute_T_fdtd(slab_eps=slab_eps0, slab_ds=slab_ds0) -> float:\n", " \"\"\"Given the slab permittivities and thicknesses, compute T.\"\"\"\n", " sim = make_sim(slab_eps=slab_eps, slab_ds=slab_ds)\n", " sim_data = web.run(sim, task_name=\"slab\", verbose=True)\n", " return post_process_T(sim_data)" ], "id": "94c8c21d" }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "
02:29:01 UTC Created task 'slab' with resource_id                               \n",
              "             'fdve-28fe4774-6d42-4456-9953-92313445efeb' and task_type 'FDTD'.  \n",
              "
" ] } }, { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "
             View task using web UI at                                          \n",
              "             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-28fe4774-6d4\n",
              "             2-4456-9953-92313445efeb'.                                         \n",
              "
" ] } }, { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "
             Task folder: 'default'.                                            \n",
              "
" ] } }, { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "" ] } }, { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
02:29:02 UTC Estimated FlexCredit cost: 0.107. This assumes the FDTD solver runs\n",
              "             for the full simulation time; if early shutoff is reached, the     \n",
              "             billed cost can be lower. Use 'web.real_cost(task_id)' to get the  \n",
              "             billed FlexCredit cost after a simulation run.                     \n",
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
02:29:03 UTC status = success                                                   \n",
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              ""
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
02:29:04 UTC Loading results from simulation_data.hdf5                          \n",
              "
"
            ]
          }
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0.7851578758048251"
          ]
        }
      ],
      "source": [
        "T_fdtd = compute_T_fdtd(slab_eps0, slab_ds0)\n",
        "print(T_fdtd)"
      ],
      "id": "e976bb1a"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Computing T and Gradient with FDTD\n",
        "\n",
        "Now that we have this function defined, we are ready to compute our\n",
        "transmission and gradients using Tidy3d.\n",
        "\n",
        "We first call `autograd.value_and_grad()` on our transmission\n",
        "calculation function, which returns a function that will give us both\n",
        "`T` and the gradient of `T` with respect to the input parameters in one\n",
        "shot. For more details, see the previous tutorial."
      ],
      "id": "e8b2071b-54ce-462a-81aa-d830e7bd622e"
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {},
      "outputs": [],
      "source": [
        "compute_T_and_grad_fdtd = ag.value_and_grad(compute_T_fdtd, argnum=(0, 1))"
      ],
      "id": "c9f17ffe"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we call this function on our starting parameters, which will kick\n",
        "off the original (`fwd`) T transmission simulation and then the reverse\n",
        "(`adj`) simulation, which is used in combination with `fwd` for the\n",
        "gradient calculation."
      ],
      "id": "09f4a3fa-8190-4069-9c8f-62b6343b058d"
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
             Created task 'slab' with resource_id                               \n",
              "             'fdve-69bd8f13-dac7-4fa0-915d-763ad48da03d' and task_type 'FDTD'.  \n",
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
             View task using web UI at                                          \n",
              "             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-69bd8f13-dac\n",
              "             7-4fa0-915d-763ad48da03d'.                                         \n",
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
             Task folder: 'default'.                                            \n",
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              ""
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              ""
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
02:29:07 UTC status = success                                                   \n",
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              ""
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
02:29:08 UTC Loading results from simulation_data.hdf5                          \n",
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
             Started working on Batch containing 1 tasks.                       \n",
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
02:29:11 UTC Maximum FlexCredit cost: 0.108 for the whole batch.                \n",
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
             Use 'Batch.real_cost()' to get the billed FlexCredit cost after    \n",
              "             completion.                                                        \n",
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              ""
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
             Batch complete.                                                    \n",
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
"
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              ""
            ]
          }
        },
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "
"
            ]
          }
        }
      ],
      "source": [
        "T_fdtd, (grad_eps_fdtd, grad_ds_fdtd) = compute_T_and_grad_fdtd(slab_eps0, slab_ds0)"
      ],
      "id": "20a593ed"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Checking Accuracy of TMM (Numerical) vs FDTD (Adjoint)\n",
        "\n",
        "Let’s convert these from autograd types to numpy arrays to work with\n",
        "them easier, and then display the results compared to TMM."
      ],
      "id": "29af200c-0c1d-4884-abfc-120ac46286db"
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {},
      "outputs": [],
      "source": [
        "grad_eps_fdtd = np.array(grad_eps_fdtd)\n",
        "grad_ds_fdtd = np.array(grad_ds_fdtd)"
      ],
      "id": "160fc08a"
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "T (tmm)  = 0.78581\n",
            "T (FDTD) = 0.78516"
          ]
        }
      ],
      "source": [
        "print(f\"T (tmm)  = {T_tmm:.5f}\")\n",
        "print(f\"T (FDTD) = {T_fdtd:.5f}\")"
      ],
      "id": "d9ed73b8"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We see that the transmission results match very well with TMM, giving us\n",
        "a lot of confidence that our set up is correct.\n",
        "\n",
        "Let’s look at the gradients now."
      ],
      "id": "4f77688e-cfee-49e8-a028-fc3f15291d48"
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "un-normalized:\n",
            "    grad_eps (tmm)  = [-0.2766323   0.01377339 -0.2032054  -0.28999361]\n",
            "    grad_eps (FDTD)  = [-0.28167686  0.01413911 -0.20659503 -0.29535051]\n",
            "--------------------------------------------------------------------------------\n",
            "    grad_ds  (tmm)  = [-1.75199732 -0.21552416  1.00729645 -2.08209951]\n",
            "    grad_ds  (FDTD)  = [-1.74894425 -0.21534298  1.00708371 -2.07979115]\n",
            "RMS error = 1.8039225960236371 %\n",
            "RMS error = 0.13189793532981076 %"
          ]
        }
      ],
      "source": [
        "print(\"un-normalized:\")\n",
        "print(f\"\\tgrad_eps (tmm)  = {grad_eps_tmm}\")\n",
        "print(f\"\\tgrad_eps (FDTD)  = {grad_eps_fdtd}\")\n",
        "print(80 * \"-\")\n",
        "print(f\"\\tgrad_ds  (tmm)  = {grad_ds_tmm}\")\n",
        "print(f\"\\tgrad_ds  (FDTD)  = {grad_ds_fdtd}\")\n",
        "\n",
        "\n",
        "rms_eps = np.linalg.norm(grad_eps_tmm - grad_eps_fdtd) / np.linalg.norm(grad_eps_tmm)\n",
        "rms_ds = np.linalg.norm(grad_ds_tmm - grad_ds_fdtd) / np.linalg.norm(grad_ds_tmm)\n",
        "\n",
        "print(f\"RMS error = {rms_eps * 100} %\")\n",
        "print(f\"RMS error = {rms_ds * 100} %\")"
      ],
      "id": "25764620"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "If we only care about the error in the “directions” of the gradients, we\n",
        "can compare their normalized versions to each other."
      ],
      "id": "d4efebe5-9ce6-4aa1-9f66-7521170493c0"
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {},
      "outputs": [],
      "source": [
        "def normalize(arr):\n",
        "    return arr / np.linalg.norm(arr)\n",
        "\n",
        "\n",
        "grad_eps_tmm_norm = normalize(grad_eps_tmm)\n",
        "grad_ds_tmm_norm = normalize(grad_ds_tmm)\n",
        "grad_eps_fdtd_norm = normalize(grad_eps_fdtd)\n",
        "grad_ds_fdtd_norm = normalize(grad_ds_fdtd)\n",
        "\n",
        "rms_eps = np.linalg.norm(grad_eps_tmm_norm - grad_eps_fdtd_norm) / np.linalg.norm(grad_eps_tmm_norm)\n",
        "rms_ds = np.linalg.norm(grad_ds_tmm_norm - grad_ds_fdtd_norm) / np.linalg.norm(grad_ds_tmm_norm)"
      ],
      "id": "a2991c98"
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "normalized:\n",
            "    grad_eps (tmm)  = [-0.61534061  0.03063751 -0.45200988 -0.64506151]\n",
            "    grad_eps (FDTD)  = [-0.61546818  0.03089418 -0.45141325 -0.6453453 ]\n",
            "    RMS error = 0.0720183276321135 %\n",
            "--------------------------------------------------------------------------------\n",
            "    grad_ds  (tmm)  = [-0.60214521 -0.07407365  0.34619844 -0.71559827]\n",
            "    grad_ds  (FDTD)  = [-0.60183578 -0.07410248  0.34655136 -0.71568475]\n",
            "    RMS error = 0.047813317028530644 %"
          ]
        }
      ],
      "source": [
        "print(\"normalized:\")\n",
        "print(f\"\\tgrad_eps (tmm)  = {grad_eps_tmm_norm}\")\n",
        "print(f\"\\tgrad_eps (FDTD)  = {grad_eps_fdtd_norm}\")\n",
        "print(f\"\\tRMS error = {rms_eps * 100} %\")\n",
        "print(80 * \"-\")\n",
        "print(f\"\\tgrad_ds  (tmm)  = {grad_ds_tmm_norm}\")\n",
        "print(f\"\\tgrad_ds  (FDTD)  = {grad_ds_fdtd_norm}\")\n",
        "print(f\"\\tRMS error = {rms_ds * 100} %\")"
      ],
      "id": "ed209e78"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In which case we see a very good agreement."
      ],
      "id": "9dcfade8-dc08-4a3b-bc77-f7121d7cab57"
    }
  ],
  "nbformat": 4,
  "nbformat_minor": 5,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    }
  }
}