{
"cells": [
{
"cell_type": "markdown",
"id": "60244b2c",
"metadata": {},
"source": [
"## Photonic crystal slab band structure calculation\n",
"\n",
"In this notebook, we simulate a photonic crystal slab consisting of a square lattice of air holes in a dielectric slab. Our goal is to compute the band structure of this photonic crystal slab, as found in:\n",
"\n",
"> Shanhui Fan and J. D. Joannopoulos, \"Analysis of guided resonances in photonic crystal slabs,\" Phys. Rev. B 65, 235112 (2002).\n",
"\n",
"To this end, we excite the structure with several `PointDipole` sources, and we measure the response with several `FieldTimeMonitor` monitors. We excite modes with a fixed Bloch wavevector by using Bloch boundary conditions. We then use the `ResonanceFinder` to find the resonant frequencies. By sweeping the Bloch wavevector, we obtain the full band structure of the photonic crystal slab.\n",
"\n",
"See also the api reference for `ResonanceFinder` [here](https://docs.flexcompute.com/projects/tidy3d/en/v1.9.0rc2/_autosummary/tidy3d.plugins.ResonanceFinder.html)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "193cb096",
"metadata": {
"execution": {
"iopub.execute_input": "2023-02-03T05:01:47.316746Z",
"iopub.status.busy": "2023-02-03T05:01:47.316281Z",
"iopub.status.idle": "2023-02-03T05:01:48.387709Z",
"shell.execute_reply": "2023-02-03T05:01:48.387322Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
[23:01:47] WARNING This version of Tidy3D was pip installed from the 'tidy3d-beta' repository on __init__.py:103\n",
" PyPI. Future releases will be uploaded to the 'tidy3d' repository. From now on, \n",
" please use 'pip install tidy3d' instead. \n",
"
\n"
],
"text/plain": [
"\u001b[2;36m[23:01:47]\u001b[0m\u001b[2;36m \u001b[0m\u001b[31mWARNING \u001b[0m This version of Tidy3D was pip installed from the \u001b[32m'tidy3d-beta'\u001b[0m repository on \u001b]8;id=155674;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=358327;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/__init__.py#103\u001b\\\u001b[2m103\u001b[0m\u001b]8;;\u001b\\\n",
"\u001b[2;36m \u001b[0m PyPI. Future releases will be uploaded to the \u001b[32m'tidy3d'\u001b[0m repository. From now on, \u001b[2m \u001b[0m\n",
"\u001b[2;36m \u001b[0m please use \u001b[32m'pip install tidy3d'\u001b[0m instead. \u001b[2m \u001b[0m\n"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"
INFO Using client version: 1.9.0rc1 __init__.py:121\n",
"
\n"
],
"text/plain": [
"\u001b[2;36m \u001b[0m\u001b[2;36m \u001b[0m\u001b[34mINFO \u001b[0m Using client version: \u001b[1;36m1.9\u001b[0m.0rc1 \u001b]8;id=472395;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=736079;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/__init__.py#121\u001b\\\u001b[2m121\u001b[0m\u001b]8;;\u001b\\\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# standard python imports\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import xarray as xr\n",
"\n",
"import tidy3d as td\n",
"from tidy3d import web\n",
"from tidy3d.plugins import ResonanceFinder\n"
]
},
{
"cell_type": "markdown",
"id": "6022fa5a",
"metadata": {},
"source": [
"We will randomly position our sources and monitors. Here, we seed the random number generator to guarantee reproducible results."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "be1a51fc",
"metadata": {
"execution": {
"iopub.execute_input": "2023-02-03T05:01:48.389564Z",
"iopub.status.busy": "2023-02-03T05:01:48.389376Z",
"iopub.status.idle": "2023-02-03T05:01:48.391205Z",
"shell.execute_reply": "2023-02-03T05:01:48.390929Z"
}
},
"outputs": [],
"source": [
"rng = np.random.default_rng(12345)\n"
]
},
{
"cell_type": "markdown",
"id": "5b3da68c",
"metadata": {},
"source": [
"Now we define the parameters for the simulation, structure, sources, and monitors.\n",
"\n",
"We take the dipole polarization to be `Hz` and the symmetry to be `(0,0,1)` in order to excite only modes which are even with respect to the xy mirror plane."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "da1f5f85",
"metadata": {
"execution": {
"iopub.execute_input": "2023-02-03T05:01:48.392641Z",
"iopub.status.busy": "2023-02-03T05:01:48.392543Z",
"iopub.status.idle": "2023-02-03T05:01:48.395993Z",
"shell.execute_reply": "2023-02-03T05:01:48.395719Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Total runtime = 6.74 ps\n"
]
}
],
"source": [
"# Simulation parameters\n",
"runtime_fwidth = 200.0 # in units of 1/frequency bandwidth of the source\n",
"t_start_fwidth = 5.0 # time to start monitoring after source has decayed, units of 1/frequency bandwidth\n",
"dPML = 1.0 # space between PhC slabs and PML, in unit of longest wavelength of interest\n",
"\n",
"# Structure parameters (um)\n",
"a_lattice = 1 # lattice constant \"a\"\n",
"r_hole = 0.2 * a_lattice # radius of the air holes\n",
"t_slab = 0.5 * a_lattice # slab thickness\n",
"ep_slab = 12 # dielectric constant of the slab\n",
"ep_hole = 1 # dielectric constant of the holes\n",
"\n",
"# Frequency range of interest (Hz)\n",
"freq_range_unitless = np.array((0.1, 0.43)) # in units of c/a\n",
"freq_scale = (\n",
" td.constants.C_0 / a_lattice\n",
") # frequency scale determined by the lattice constant\n",
"freq_range = freq_range_unitless * freq_scale\n",
"lambda_range = (td.constants.C_0 / freq_range[1], td.constants.C_0 / freq_range[0])\n",
"\n",
"# Gaussian pulse parameters\n",
"freq0 = np.sum(freq_range) / 2 # central frequency\n",
"freqw = 0.3 * (freq_range[1] - freq_range[0]) # pulse width\n",
"\n",
"# Runtime\n",
"run_time = runtime_fwidth / freqw\n",
"print(f\"Total runtime = {(run_time*1e12):.2f} ps\")\n",
"t_start = t_start_fwidth / freqw\n",
"\n",
"# Simulation size\n",
"spacing = dPML * lambda_range[-1] # space between PhC slabs and PML\n",
"sim_size = Lx, Ly, Lz = (a_lattice, a_lattice, 2 * spacing + t_slab)\n",
"\n",
"# Number of k values to sample, per edge of the irreducible Brillouin zone\n",
"Nk = 4\n",
"\n",
"# Number of dipoles and monitors\n",
"num_dipoles = 7\n",
"num_monitors = 2\n",
"\n",
"# Dipole polarization and symmetry\n",
"polarization = \"Hz\"\n",
"symmetry = (0, 0, 1)\n"
]
},
{
"cell_type": "markdown",
"id": "27c57eee",
"metadata": {},
"source": [
"We define the materials and structures in terms of the above parameters. We take the photonic crystal slab to be lying in the xy plane. Because the photonic crystal slab is periodic in the x and y directions, we only need to simulate a single unit cell, containing the dielectric slab and a single air hole."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "cbcc09f9",
"metadata": {
"execution": {
"iopub.execute_input": "2023-02-03T05:01:48.397427Z",
"iopub.status.busy": "2023-02-03T05:01:48.397340Z",
"iopub.status.idle": "2023-02-03T05:01:48.399959Z",
"shell.execute_reply": "2023-02-03T05:01:48.399715Z"
}
},
"outputs": [],
"source": [
"mat_slab = td.Medium(permittivity=ep_slab, name=\"mat_slab\")\n",
"mat_hole = td.Medium(permittivity=ep_hole, name=\"mat_hole\")\n",
"\n",
"slab = td.Structure(\n",
" geometry=td.Box(\n",
" center=(0, 0, 0),\n",
" size=(td.inf, td.inf, t_slab),\n",
" ),\n",
" medium=mat_slab,\n",
" name=\"slab\",\n",
")\n",
"\n",
"hole = td.Structure(\n",
" geometry=td.Cylinder(\n",
" center=(0, 0, 0),\n",
" axis=2,\n",
" radius=r_hole,\n",
" length=t_slab,\n",
" ),\n",
" medium=mat_hole,\n",
" name=\"hole\",\n",
")\n",
"\n",
"structures = [slab, hole]\n"
]
},
{
"cell_type": "markdown",
"id": "dbda5f33",
"metadata": {},
"source": [
"We will excite the photonic crystal slab with several `PointDipole` sources. Each dipole will have a random position and phase."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "7cd6c4ed",
"metadata": {
"execution": {
"iopub.execute_input": "2023-02-03T05:01:48.401362Z",
"iopub.status.busy": "2023-02-03T05:01:48.401282Z",
"iopub.status.idle": "2023-02-03T05:01:48.405869Z",
"shell.execute_reply": "2023-02-03T05:01:48.405610Z"
}
},
"outputs": [],
"source": [
"dipole_positions = rng.uniform(\n",
" [-Lx / 2, -Ly / 2, 0], [Lx / 2, Ly / 2, 0], [num_dipoles, 3]\n",
")\n",
"\n",
"dipole_phases = rng.uniform(0, 2 * np.pi, num_dipoles)\n",
"\n",
"pulses = []\n",
"dipoles = []\n",
"for i in range(num_dipoles):\n",
" pulse = td.GaussianPulse(freq0=freq0, fwidth=freqw, phase=dipole_phases[i])\n",
" pulses.append(pulse)\n",
" dipoles.append(\n",
" td.PointDipole(\n",
" source_time=pulse,\n",
" center=tuple(dipole_positions[i]),\n",
" polarization=polarization,\n",
" name=\"dipole_\" + str(i),\n",
" )\n",
" )\n"
]
},
{
"cell_type": "markdown",
"id": "700880a7",
"metadata": {},
"source": [
"We create `FieldTimeMonitors` to record the field as a function of time at several random locations within the photonic crystal slab. Crucially, we start the monitors after the source pulse has decayed."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "bc6dca60",
"metadata": {
"execution": {
"iopub.execute_input": "2023-02-03T05:01:48.407272Z",
"iopub.status.busy": "2023-02-03T05:01:48.407185Z",
"iopub.status.idle": "2023-02-03T05:01:48.409629Z",
"shell.execute_reply": "2023-02-03T05:01:48.409355Z"
}
},
"outputs": [],
"source": [
"monitor_positions = rng.uniform(\n",
" [-Lx / 2, -Ly / 2, 0], [Lx / 2, Ly / 2, 0], [num_monitors, 3]\n",
")\n",
"\n",
"monitors_time = []\n",
"for i in range(num_monitors):\n",
" monitors_time.append(\n",
" td.FieldTimeMonitor(\n",
" fields=[\"Hz\"],\n",
" center=tuple(monitor_positions[i]),\n",
" size=(0, 0, 0),\n",
" start=t_start,\n",
" name=\"monitor_time_\" + str(i),\n",
" )\n",
" )\n"
]
},
{
"cell_type": "markdown",
"id": "e6939eae",
"metadata": {},
"source": [
"We will perform `3*Nk` different simulations, each with different Bloch boundary conditions, as we sweep the Bloch wavevector over the boundary of the irreducible Brillouin zone. We sweep over three lines, namely $\\Gamma X$, $XM$, and $M\\Gamma$. We use a PML in the z direction.\n",
"\n",
"Here, we simply define all of the boundary conditions we will use and put them into a single array."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "84ec0547",
"metadata": {
"execution": {
"iopub.execute_input": "2023-02-03T05:01:48.410972Z",
"iopub.status.busy": "2023-02-03T05:01:48.410887Z",
"iopub.status.idle": "2023-02-03T05:01:48.414515Z",
"shell.execute_reply": "2023-02-03T05:01:48.414247Z"
}
},
"outputs": [],
"source": [
"bspecs_gammax = []\n",
"bspecs_xm = []\n",
"bspecs_mgamma = []\n",
"for i in range(Nk):\n",
" bspecs_gammax.append(\n",
" td.BoundarySpec(\n",
" x=td.Boundary.bloch((1 / 2) * i / Nk),\n",
" y=td.Boundary.periodic(),\n",
" z=td.Boundary.pml(),\n",
" )\n",
" )\n",
" bspecs_xm.append(\n",
" td.BoundarySpec(\n",
" x=td.Boundary.bloch(1 / 2),\n",
" y=td.Boundary.bloch((1 / 2) * i / Nk),\n",
" z=td.Boundary.pml(),\n",
" )\n",
" )\n",
" bspecs_mgamma.append(\n",
" td.BoundarySpec(\n",
" x=td.Boundary.bloch((1 / 2) * (1 - i / Nk)),\n",
" y=td.Boundary.bloch((1 / 2) * (1 - i / Nk)),\n",
" z=td.Boundary.pml(),\n",
" )\n",
" )\n",
"bspecs = bspecs_gammax + bspecs_xm + bspecs_mgamma\n"
]
},
{
"cell_type": "markdown",
"id": "e80d5995",
"metadata": {},
"source": [
"Now we define the simulations we want to run. "
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "0c0a8a8a",
"metadata": {
"execution": {
"iopub.execute_input": "2023-02-03T05:01:48.415791Z",
"iopub.status.busy": "2023-02-03T05:01:48.415704Z",
"iopub.status.idle": "2023-02-03T05:01:48.423035Z",
"shell.execute_reply": "2023-02-03T05:01:48.422787Z"
}
},
"outputs": [],
"source": [
"sims = {}\n",
"for i in range(3 * Nk):\n",
" sims[f\"sim_{i}\"] = td.Simulation(\n",
" center=(0, 0, 0),\n",
" size=sim_size,\n",
" grid_spec=td.GridSpec.auto(),\n",
" structures=structures,\n",
" sources=dipoles,\n",
" monitors=monitors_time,\n",
" run_time=run_time,\n",
" shutoff=0,\n",
" boundary_spec=bspecs[i],\n",
" normalize_index=None,\n",
" symmetry=symmetry,\n",
" )\n"
]
},
{
"cell_type": "markdown",
"id": "22ff9910",
"metadata": {},
"source": [
"Let's check that the structure and source look correct. The source spectrum must fill the entire frequency range of interest."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "443a41b7",
"metadata": {
"execution": {
"iopub.execute_input": "2023-02-03T05:01:48.424363Z",
"iopub.status.busy": "2023-02-03T05:01:48.424280Z",
"iopub.status.idle": "2023-02-03T05:01:48.775003Z",
"shell.execute_reply": "2023-02-03T05:01:48.774082Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
[23:01:48] INFO Auto meshing using wavelength 3.7736 defined from sources. grid_spec.py:510\n",
"
\n"
],
"text/plain": [
"\u001b[2;36m[23:01:48]\u001b[0m\u001b[2;36m \u001b[0m\u001b[34mINFO \u001b[0m Auto meshing using wavelength \u001b[1;36m3.7736\u001b[0m defined from sources. \u001b]8;id=302674;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/components/grid/grid_spec.py\u001b\\\u001b[2mgrid_spec.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=808338;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/components/grid/grid_spec.py#510\u001b\\\u001b[2m510\u001b[0m\u001b]8;;\u001b\\\n"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAuoAAAGGCAYAAAAtuaQaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABNiElEQVR4nO3deXxU1f3/8fcsmZnsIRASwg4iCCIgCEKlFI2CUAutUqVU0VKoC1qlWuGnBcWfolaRqihCtdR9+1bqVy2KCFoFQRH6QxRFBFkDxJCFLDOZmfP7AzNlspFAkrkzeT0fj/uAuXPu5HMzcO97zpx7rs0YYwQAAADAUuyRLgAAAABAdQR1AAAAwIII6gAAAIAFEdQBAAAACyKoAwAAABZEUAcAAAAsiKAOAAAAWBBBHQAAALAggjoAAABgQQR1oBmtXr1aNptNq1evjnQpAIAoxzkl9hHUgSbw2GOPaenSpZEu44Q8//zzWrBgQaTLkCQFg0Hdf//96tq1qzwej8444wy98MIL9d6+oKBA06ZNU0ZGhhITEzVy5Eh99tlnNbZ9/fXXdeaZZ8rj8ahTp06aM2eO/H5/Y+0KAJwwzimN68knn9Rpp50mj8ejHj166JFHHqn3tl6vV7feequys7MVHx+vIUOGaMWKFTW2XbNmjc455xwlJCQoKytLN9xwg44cOdKwYg2ARtenTx8zYsSIausDgYApKyszgUCg+Yuqp7Fjx5rOnTtHugxjjDEzZ840kszUqVPN4sWLzdixY40k88ILLxx320AgYIYNG2YSExPNHXfcYR599FHTu3dvk5ycbL7++uuwtm+99Zax2Wxm5MiRZvHixeb66683drvdXH311U21awBQb5xTGs+iRYuMJHPxxRebxYsXm8svv9xIMvfee2+9tr/sssuM0+k0N998s3niiSfM0KFDjdPpNP/+97/D2m3cuNF4PB4zYMAA8/jjj5vbbrvNuN1uM3r06AbVS1BHkzly5EikS4iY2g6q0cAqB9U9e/aYuLg4c91114XWBYNBM3z4cNOhQwfj9/vr3P6ll14ykswrr7wSWnfw4EGTlpZmJk6cGNa2d+/epl+/fqaioiK07rbbbjM2m818+eWXjbRHAE4G55QRkS7jhFjlnGKMMaWlpaZ169Zm7NixYesnTZpkEhMTTX5+fp3br1u3zkgyf/7zn0PrysrKTPfu3c3QoUPD2l544YWmXbt2prCwMLRuyZIlRpJ5++23610zQR31smfPHvOb3/zGtGvXzrhcLtOlSxdz9dVXG6/Xa4wx5m9/+5uRZFavXm2uueYak5GRYdLS0kLbL1y40PTu3du4XC7Trl07c+2115rDhw+H/Yyvv/7a/OIXvzCZmZnG7Xab9u3bm0svvdQUFBSE2rzzzjvmRz/6kUlNTTWJiYnm1FNPNbNmzTpu/fXZrry83MyePdt0797duFwu06FDB3PLLbeY8vLyaq/3zDPPmLPOOsvEx8ebtLQ0M3z48NB/vM6dOxtJYUvlAXbVqlVGklm1alXY67388svmzDPPNB6Px7Ru3dpMmjTJ7NmzJ6zN5MmTTWJiotmzZ48ZN26cSUxMNG3atDF/+MMfjhtajTFm2bJlZsyYMaH3sFu3bmbu3Llh244YMaJa7XUdYCdPnlytfeUyZ86c49ZUl4ULFxpJZsuWLWHrn3/+eSOpWu9FVRMmTDCZmZnVepqmTZtmEhISQu/rli1bjCSzcOHCsHZ79+41ksxdd911UvsBoDrOKeE4pxz11FNPGUnmySefDFt/9913G0nmzTffPG5ddXnzzTdrfJ01a9YYSeaZZ56pc/tbbrnFOByOsPBtjDH33HOPkWR27dpljDGmsLDQOJ1Oc8stt4S183q9JikpyUyZMqXeNTsbNlAGLdG+ffs0ePDg0HjfXr16ae/evXr11VdVWloql8sVanvttdcqIyNDs2fPVklJiSTpjjvu0J133qmcnBxdc801+uqrr/T444/rk08+0UcffaS4uDj5fD6NGjVKXq9X119/vbKysrR371698cYbKigoUGpqqrZs2aKf/vSnOuOMMzR37ly53W598803+uijj+qsvz7bBYNB/exnP9OHH36oadOm6bTTTtPmzZv10EMP6euvv9ayZctCbe+8807dcccdGjZsmObOnSuXy6V169bpvffe0wUXXKAFCxbo+uuvV1JSkm677TZJUmZmZq31LV26VFdddZXOOusszZs3TwcOHNBf/vIXffTRR9q4caPS0tJCbQOBgEaNGqUhQ4bogQce0LvvvqsHH3xQ3bt31zXXXFPn72Hp0qVKSkrSjBkzlJSUpPfee0+zZ89WUVGR/vznP0uSbrvtNhUWFmrPnj166KGHJElJSUm1vubvfvc75eTkhK1bvny5nnvuObVt2za0Li8vr87aKiUnJ8vtdkuSNm7cqMTERJ122mlhbQYPHhx6/pxzzqn1tTZu3KgzzzxTdnv4pTiDBw/W4sWL9fXXX6tv377auHGjJGnQoEFh7bKzs9WhQ4fQ8wAaB+cUzim1ueqqq/SPf/xDM2bM0Pnnn6+OHTtq8+bNuvPOOzVlyhSNGTMm1Pbw4cMKBAJ11ihJCQkJSkhIkKRaj/cDBw6U3W7Xxo0b9etf/7rW19q4caNOPfVUpaSkhK2vPC9t2rQpVLPf76/2c1wul/r379+w80q9Iz1arCuuuMLY7XbzySefVHsuGAwaY/7b+3HOOeeEfZo+ePCgcblc5oILLgjr2Xz00UeNJPPUU08ZY46O5VKVYQpVPfTQQ0aSOXToUIPqr892zzzzjLHb7dV6aSvHsn300UfGGGO2bdtm7Ha7+fnPf16tp7byd2FM7V9TVu398Pl8pm3btub00083ZWVloXZvvPGGkWRmz54dWlfZez137tyw1xwwYIAZOHBg3b8Ec/Qrv6p+97vfhfUuG3NyX1Nu27bNpKammvPPPz/s34Fq6XWvuvztb38Lq6Nbt27VfkZJSYmRZGbOnFlnLYmJieY3v/lNtfWVPSrLly83xhjz5z//Oawn5FhnnXWWOfvss+u7+wDqgXMK55S67N+/36Snp5vzzz/feL1eM2DAANOpU6dqvdg1fdNQ03Lst7vXXXedcTgcNf7cjIwMc9lll9VZW58+fcy5555bbX3lN7OLFi0yxhjzyiuvGEnmgw8+qNZ2woQJJisr63i/hhBmfUGdgsGgli1bposuuqjaJ0NJstlsYY+nTp0qh8MRevzuu+/K5/PpxhtvDOvZnDp1qlJSUvTmm29KklJTUyVJb7/9tkpLS2uspbIX4J///KeCwWC996E+273yyis67bTT1KtXL+Xl5YWWc889V5K0atUqSdKyZcsUDAY1e/bsaj21VX8X9fHpp5/q4MGDuvbaa+XxeELrx44dq169eoV+P8e6+uqrwx4PHz5c33777XF/Vnx8fOjvxcXFysvL0/Dhw1VaWqqtW7c2uPaqSkpK9POf/1ytWrXSCy+8EPbvYMWKFfVaRo0aFdqmrKws1Lt+rMrfU1lZWZ311Hf7yj9ra3u8nwOg/jincE45nqysLC1cuFArVqzQ8OHDtWnTJj311FPVerGfe+65ep1XrrjiitA2ZWVlYd/YHKs+x/tInFcY+oI6HTp0SEVFRTr99NPr1b5r165hj7/77jtJUs+ePcPWu1wudevWLfR8165dNWPGDM2fP1/PPfechg8frp/97Gf69a9/HTrgXnrppfrrX/+q3/72t5o5c6bOO+88/eIXv9All1xS7QB3rPpst23bNn355ZfKyMio8TUOHjwoSdq+fbvsdrt69+5dr9/H8dT2+5GkXr166cMPPwxb5/F4qtXYqlUrHT58+Lg/a8uWLbr99tv13nvvqaioKOy5wsLChpZezdSpU7V9+3atWbNGrVu3Dnuu6vCY+oiPj5fX6622vry8PPR8Y2xf+WdtbY/3cwDUH+eUozin1O2yyy7Ts88+qzfffFPTpk3TeeedV63Nj370owa/bnx8vHw+X43P1ed4H4nzCkEdjepkQs2DDz6oK6+8Uv/85z/1zjvv6IYbbtC8efP08ccfq0OHDoqPj9cHH3ygVatW6c0339Ty5cv10ksv6dxzz9U777wT1utStabjbRcMBtW3b1/Nnz+/xtfo2LHjCe9XY6ptH4+noKBAI0aMUEpKiubOnavu3bvL4/Hos88+06233tqg3qSa/OUvf9ELL7ygZ599Vv3796/2fG5ubr1eJzU1NfRvqF27dlq1apWMMWE9S/v375d0dAx5Xdq1axdqe6yq27dr1y60vur7vH///tDYQwDNj3NK07LqOeX777/Xp59+Kkn64osvFAwGq314OnToUL3GqCclJYXGxbdr106BQEAHDx4Mu47K5/Pp+++/r9d5Ze/evdXW13Veqant8X5OmHoPkkGLFAgETEpKihk3blyd7SrHE1Ydc1g5Q8dbb70Vtt7r9ZrU1FRz8cUX1/qaH330kZFkbrvttlrbVF4JvmLFiuPvTB3bjRkzxrRv3z5sTGBNKsczb9y4sc52p59+er3GE1Zeaf7YY49Va3vaaaeFjROsvEK/qjlz5pjj/Vd+7bXXjCTz/vvvh61fvHhxtRkDfvrTnzZoPOEHH3xgnE6nufHGG2ttoxMYo1455rTqrC/PPfdcrWP/jnXJJZfUOOvL1KlTw8ZQfv7553XO+lJ1/CaAE8c5JRznlJpdeumlJiEhwcybN89IMg8++GC1NicyRr1yrH7VWV8q/208/fTTddZ188031zjrS+X7X3mtU0FBQZ2zvtR0/VRtGKOOOtntdo0fP17/+7//G/p0eyxjTJ3b5+TkyOVy6eGHHw5r++STT6qwsFBjx46VJBUVFVW7C2Tfvn1lt9tDXx3l5+dXe/3K3tuavl6qVJ/tfvnLX2rv3r1asmRJtbZlZWWh2QbGjx8vu92uuXPnVusxOHb/EhMTVVBQUGtNlQYNGqS2bdtq0aJFYfvwr3/9S19++WXo93OyKntNjq3R5/Ppscceq9Y2MTGx3l9b7t+/X7/85S91zjnnhK7yr8mJjFEfN26c4uLiwmo0xmjRokVq3769hg0bFlbH1q1bVVFREVp3ySWX6MCBA/rHP/4RWpeXl6dXXnlFF110UWjsYJ8+fdSrVy8tXrw4rHfm8ccfl81m0yWXXFKv3wWA4+OcwjnleF599VW99NJLuvfeezVz5kxddtlluv322/X111+HtTuRMernnnuu0tPT9fjjj4e91uOPP66EhISw309eXp62bt0ado3DJZdcokAgoMWLF4fWeb1e/e1vf9OQIUNC35SkpqYqJydHzz77rIqLi0Ntn3nmGR05ckQTJkyo9++DHnUc1549e0xWVpZJSEgwN954o3niiSfMHXfcYfr06ROat7a23g9j/vvp/IILLjCPPvqouf76643D4TBnnXWW8fl8xpijn87bt29vbrzxRvPYY4+Zhx9+2Jx11lkmLi7OrF271hhjzO9//3szYMAAc/vtt5slS5aYu+++27Rv39506NAhbF7cquqzXSAQMGPGjDE2m81cdtll5pFHHjELFiwwV199tUlPTw/brz/96U9Gkhk2bJh54IEHzCOPPGKuuOKKsFlIrr32WmOz2cxdd91lXnjhBbNy5UpjTM1z3lb+7oYMGWIWLFhgZs2aZRISEkyXLl3C5gU+md6PvLw806pVK9O5c2fz4IMPmvnz55sBAwaYfv36Vavn/vvvN5LMTTfdZJ5//nnz+uuv1/q6l1xyiXE4HGb+/PnmmWeeCVv+85//1FlTfdxyyy1Gkpk2bZpZsmRJ6M6kzz33XFi7ytkLduzYEVrn9/vN2WefbZKSksydd95pFi5caPr06WOSk5PN1q1bw7b/3//9X2Oz2cy5555rFi9ebG644QZjt9vN1KlTT3ofAITjnMI5pTYHDhwwbdq0MSNHjgx9G5GXl2cyMzPN0KFDG+UOrJX36LjkkkvMkiVLzBVXXGEkmbvvvjusXeXvoeoc9RMmTAj1lj/xxBNm2LBhxul0Vvt2YcOGDcbtdofdmdTj8ZgLLrigQfUS1FEv3333nbniiitMRkaGcbvdplu3bua6666rdnOKmg6qxhwdxtCrVy8TFxdnMjMzzTXXXBN2wPj222/Nb37zG9O9e3fj8XhMenq6GTlypHn33XdDbVauXGnGjRtnsrOzjcvlMtnZ2WbixInVbgdfVX238/l85r777jN9+vQxbrfbtGrVygwcONDceeed1b7meuqpp8yAAQNC7UaMGBH2VWlubq4ZO3asSU5ONtLxb07x0ksvhV4vPT29zptTVFWfg6oxR7/aO/vss018fLzJzs42f/zjH83bb79drZ4jR46YX/3qVyYtLc1Idd+coqabWVQuJ3vDI2OOnuzuuece07lzZ+NyuUyfPn3Ms88+W61dTUHdGGPy8/PNlClTTOvWrU1CQoIZMWJErf9GX3vtNdO/f3/jdrtNhw4dzO233x466QNoXJxTOKfU5Be/+IVJTk42O3fuDFv/z3/+00gy991333Hrqo/Fixebnj17GpfLZbp3724eeuihasOUagvqZWVl5uabbzZZWVnG7Xabs846KzTdb1X//ve/zbBhw4zH4zEZGRnmuuuuM0VFRQ2q1WbMcb5nAgAAANDsGKMOAAAAWBBBHQAAALAggjoAAABgQQR1AAAAwIII6gAAAIAFEdQBAAAAC3JGugCrCwaD2rdvn5KTk2Wz2SJdDgA0O2OMiouLlZ2dLbud/h0gGjQkv/B/3LoI6sexb98+de7UWUETrLWN3Wav83kAsLL6HsN2796tDh06NENFAE7Wvn37Qre0ry/+j1sPQf04kpOTFTRBPTD6PnVP717t+e3523Xz8lv1s1vnq03H6s8DgJXl7d6u1++bUesxTvrvcS45ObmZqwNwoir/v86ePVvt2rWrs21xcbH++Mc/8n/cggjqx1H5dVH39O7q07Z3re3i3PFyxSc2V1kA0Cji3PGSjn+Mk8TwPyCKVP5/dTgcstvt9Qrh/B+3HgYiAQAAxKjExEQVFRWpuLg40qXgBNCjDgAAEKMSExNlt9tVVFQkSQxviTIEdQAAgBhWGc4J69GHoA4AABDjCOvRiaAOAADQAhDWow9BHQAAoIUgrEcXgjoAAEALUjWsw7oI6gAAAC3MsWGd+dOti6AOAADQAlWG9dzc3AhXgtpwwyMAAIAWKjk5WYmJ3FndqgjqAAAALRhB3boI6gAAAIAFEdQBAAAACyKoAwAAABZEUAcAAAAsiKAOAAAAWBBBHQAAALAggjoAAABgQQR1AAAAwIII6gAAAIAFEdQBAABasJKSkkiXgFoQ1AEAAFqo4uJigrqFEdQBAABaoOLiYhUVFSkxMTHSpaAWURfUFy5cqC5dusjj8WjIkCFav359vbZ78cUXZbPZNH78+KYtEAAAwOIqQ3pKSgpB3cKiKqi/9NJLmjFjhubMmaPPPvtM/fr106hRo3Tw4ME6t9u5c6duvvlmDR8+vJkqBQAAsKZjQ3pycnKky0Edoiqoz58/X1OnTtVVV12l3r17a9GiRUpISNBTTz1V6zaBQECTJk3SnXfeqW7dujVjtQAAANZCSI8uURPUfT6fNmzYoJycnNA6u92unJwcrV27ttbt5s6dq7Zt22rKlCnNUSYAAIAlEdKjjzPSBdRXXl6eAoGAMjMzw9ZnZmZq69atNW7z4Ycf6sknn9SmTZvq/XO8Xq+8Xm/ocVFR0QnVCwAAYBWE9OgUNT3qDVVcXKzLL79cS5YsUZs2beq93bx585SamhpaOnbs2IRVAgAANC1CevSKmh71Nm3ayOFw6MCBA2HrDxw4oKysrGrtt2/frp07d+qiiy4KrQsGg5Ikp9Opr776St27d6+23axZszRjxozQ46KiIsI6AACISiUlJTLGENKjVNT0qLtcLg0cOFArV64MrQsGg1q5cqWGDh1arX2vXr20efNmbdq0KbT87Gc/08iRI7Vp06Zaw7fb7VZKSkrYAgAAEI1KSkoI6VEsanrUJWnGjBmaPHmyBg0apMGDB2vBggUqKSnRVVddJUm64oor1L59e82bN08ej0enn3562PZpaWmSVG09AABALHK73YT0KBZVQf3SSy/VoUOHNHv2bOXm5qp///5avnx56ALTXbt2yW6Pmi8JAAAAmpTP55PP55PL5Yp0KTgBURXUJWn69OmaPn16jc+tXr26zm2XLl3a+AUBAABYlNPpVF5entq0aUNYj0J0PwMAAMSotLQ0xcXFKS8vTz6fL9LloIEI6gAAADHKZrOpdevWhPUoRVAHAACIYXa7nbAepQjqAAAAMY6wHp0I6gAAAC0AYT36ENQBAABaCMJ6dCGoAwAAtCBVw3pFRUWkS0ItCOoAAAAtzLFhvaCgINLloBYEdQAAgGb2wQcf6KKLLlJ2drZsNpuWLVsW9rwxRrNnz1a7du0UHx+vnJwcbdu2rVFrqAzrTmfU3f+yxSCoAwAANLOSkhL169dPCxcurPH5+++/Xw8//LAWLVqkdevWKTExUaNGjVJ5eXmDfo7f75fP56t18fv9SkxMbIxdQhPgIxQAAEAzu/DCC3XhhRfW+JwxRgsWLNDtt9+ucePGSZKefvppZWZmatmyZbrsssvq/XP8fv9xw30wGKx/4WhWBHUAAAAL2bFjh3Jzc5WTkxNal5qaqiFDhmjt2rUNCurjx49XSkpKnW2Kiop00003nXC9aDoEdQAAAAvJzc2VJGVmZoatz8zMDD1XldfrldfrDT0uKiqSJKWnpx83qDNG3boYow4AABDl5s2bp9TU1NDSsWPHSJeERkBQBwAAsJCsrCxJ0oEDB8LWHzhwIPRcVbNmzVJhYWFo2b17d5PXiaZHUAcAALCQrl27KisrSytXrgytKyoq0rp16zR06NAat3G73UpJSQlbEP0YlAQAANDMjhw5om+++Sb0eMeOHdq0aZPS09PVqVMn3Xjjjfq///f/qkePHuratav+9Kc/KTs7W+PHj49c0Wh2BHUAAIBm9umnn2rkyJGhxzNmzJAkTZ48WUuXLtUf//hHlZSUaNq0aSooKNA555yj5cuXy+PxRKpkRABBHQAAoJn95Cc/kTGm1udtNpvmzp2ruXPnNmNVsBrGqAMAAAAWRFAHAAAALIigDgAAAFgQQR0AAACwIII6AAAAYEEEdQAAAMCCCOoAAACABRHUAQAAAAsiqAMAAAAWRFAHAAAALIigDgAAAFgQQR0AAACwIII6AAAAYEEEdQAAAMCCCOoAAACABRHUAQAAAAsiqAMAAAAWRFAHAAAALIigDgAAAFgQQR0AAACwIII6AAAAYEEEdQAAAMCCCOoAAACABRHUAQAAAAsiqAMAAAAWRFAHAAAALIigDgAAAFgQQR0AAACwIII6AAAAYEEEdQAAAMCCCOoAAACABUVdUF+4cKG6dOkij8ejIUOGaP369bW2XbJkiYYPH65WrVqpVatWysnJqbM9AAAAYBVRFdRfeuklzZgxQ3PmzNFnn32mfv36adSoUTp48GCN7VevXq2JEydq1apVWrt2rTp27KgLLrhAe/fubebKAQAAgIaJqqA+f/58TZ06VVdddZV69+6tRYsWKSEhQU899VSN7Z977jlde+216t+/v3r16qW//vWvCgaDWrlyZTNXDgAAADRM1AR1n8+nDRs2KCcnJ7TObrcrJydHa9eurddrlJaWqqKiQunp6U1VJgAAANAonJEuoL7y8vIUCASUmZkZtj4zM1Nbt26t12vceuutys7ODgv7VXm9Xnm93tDjoqKiEysYaCBHMCi3P6C4YFBxgYBcgaCcwaDigkEFbDZVOOyqsDt++NMun8Mhn8Mu2WyRLh0AADSBqAnqJ+vee+/Viy++qNWrV8vj8dTabt68ebrzzjubsTK0NIm+CqWXlSmt3Ku0cp9Sy71K9XqV7Kto8Gt5HQ4VeNwq8LhU6Har0OPWYY9bh+M9MgR4AACiWtQE9TZt2sjhcOjAgQNh6w8cOKCsrKw6t33ggQd077336t1339UZZ5xRZ9tZs2ZpxowZocdFRUXq2LHjiReOFi/R51N2cYmyi48ou7hEqV5fo722OxBQZkmpMktKw9Z7HQ7tS07U/qRE7U1JUj7BHQCAqBM1Qd3lcmngwIFauXKlxo8fL0mhC0OnT59e63b333+/7r77br399tsaNGjQcX+O2+2W2+1urLLREhmjNqVl6pFfoC4FRY0azOvLHQioa0GRuhYcHbrldTi0JyVJ29LTtCs1WUF71FyeAgBAixU1QV2SZsyYocmTJ2vQoEEaPHiwFixYoJKSEl111VWSpCuuuELt27fXvHnzJEn33XefZs+ereeff15dunRRbm6uJCkpKUlJSUkR2w/EpmSvTz2+P6we+QVqVe49/gbNyB0IqPvhQnU/XCivw6HtrVK1rXWa9iclMsYdAACLiqqgfumll+rQoUOaPXu2cnNz1b9/fy1fvjx0gemuXbtkP6an8PHHH5fP59Mll1wS9jpz5szRHXfc0ZylI1YZoy4FRTrjQJ6yj5REupp6cQcC6p2Xr955+Sp2xWlLRmttadtaFQ5HpEsDAADHiKqgLknTp0+vdajL6tWrwx7v3Lmz6QtCi2QzRqfkF2jA/oNKt1jveUMk+yp09t5cDcg9qC0ZbfT/MtuoPC7qDgsAEHPuuOOOapNb9OzZs94z3SE2cEYGGsAeDOq0vHz1yz2klBOYpcWq3IGgzsw9qL4HD2lrm3RtyspQicsV6bIAoEXr06eP3n333dBjp5PY1tLwjgP11LGwWOfs2huRi0ObS1zQqO/B79UrL18bs9rqP1kZCnDhKQBEhNPpPO7MdohtBHXgOJK9Pg3bvS80g0pLEBc0GrzvgHp+f1gfdsrW7tSUSJcEAC3Otm3blJ2dLY/Ho6FDh2revHnq1KlTjW25YWNsoqsMqIU9GNSZ+w7o0s+/alEh/VipXp/GbtupUd/sVHIMf5MAAFYzZMgQLV26VMuXL9fjjz+uHTt2aPjw4SouLq6x/bx585SamhpauAdMbKBHHahBarlX52//Tm3KyiNdiiV0LShS+6Ijer9LB21PT4t0OQAQ8y688MLQ38844wwNGTJEnTt31ssvv6wpU6ZUa88NG2MTQR2o4pTvD+vH3+2VKxiMdCmW4goGdf63u9S+6IjWdMyW38EXcgDQXNLS0nTqqafqm2++qfF5btgYmzjTAj9wBoIasXO3cnbsJqTXoXdevn6+9Ru14tsGAGg2R44c0fbt29WuXbtIl4JmRFAHJKWUe/WLL7fptLzDkS4lKrQuK9cvvtym7vkFkS4FAGLSzTffrPfff187d+7UmjVr9POf/1wOh0MTJ06MdGloRgx9QYvXurRMY7/eoQS/P9KlRJW4oNH53+5SQoVfmzPbRLocAIgpe/bs0cSJE/X9998rIyND55xzjj7++GNlZGREujQ0I4I6WrTsoiMa/c1OhrqchB/t3qf4igqtb58l2WyRLgcAYsKLL74Y6RJgAQR1tFjd8gt03o7dchgT6VKi3pm5h5RQ4df7XTrIENYBAGgUBHW0SL0O5WvEd3tEpGw8vb4/LHcgoBXdOito5zcLAMDJ4mJStDhdDxfqx4T0JtG1oEg/2blb4lsKAABOGkEdLUq74iM679td/MNvQqfmF+jsPfsjXQYAAFGPvIIWo3VpmUZ/s1NOenubXP8DeeqXeyjSZQAAENUI6mgRkr0+jdm2Q+4As7s0l6F79qvH98xLDwDAiSKoI+Y5gkGN+manEiuYJ725/WTnHrUpKY10GQAARCWCOmLe2Xv2qw23u48Ihzl6U6S4QCDSpQAAEHUI6ohpXQ8Xqu/B7yNdRouW6vXpx9/tYSYYAAAaiKCOmJXs9R2dKhAR1yO/UL3y8iNdBgAAUYWgjphkM0bnfbuLi0ct5Ee796kVQ5AAAKg3gjpiUp+D3yuLixgtJS5oNPy7vQyBAQCgngjqiDnxFRU6a19upMtADbKPlKhHfkGkywAAICoQ1BFzzt6Ty5AXCxu6e79cfmaBAQDgeAjqiClZxSXqyU12LC3B79egfQciXQYAAJZHUEfMsBmjc3btjXQZqIfTD+YpvbQs0mUAAGBpBHXEjG6HC7mxUZSwS/SqAwBwHAR1xAZjNGD/wUhXgQboWlCkND5YAQBQK4I6YkLHomJ606OMTVL/3EORLgMAAMsiqCMmDNhP4ItGPfIPK8nri3QZAABYEkEdUS+ruETZR0oiXQZOgMNI/Q7wIQsAgJoQ1BH1+h7Mi3QJOAm98vIVF2BedQAAqiKoI6q5/AF1KSiKdBk4CXFBo66HeQ8BAKiKoI6o1u1wgRzGRLoMnKQe+dykCgCAqgjqiGo98gsiXQIaQfuiI0rwVUS6DAAALIWgjqiV6POpfTEXkcYCu6RT+NAFAEAYgjqiVo/vCyJdAhoRQR0AgHAEdUStToXFkS4BjahtaZniKxj+AgBAJYI6opIjGFTbktJIl4FGls1QJgAAQgjqiEqZR0rlZLaXmENQBwDgvwjqiErZxUciXQKaQDveVwAAQgjqiEr0vMam9HIv49QBAPgBQR1Rxx4MKpPx6TGLD2EAABxFUEfUSfX6uBtpDGtVVh7pEgAAsASCOqJOark30iWgCaV6eX8BAJAI6ohCaQT1mJZa7ot0CQAAWAJBHVGHHtfYllbulRjaBAAAQR3Rhx7X2OYKBhXv90e6DAAAIo6gjqhDj3rs48MYAAAEdUQhtz8Q6RLQxNwB3mMAAKIuqC9cuFBdunSRx+PRkCFDtH79+jrbv/LKK+rVq5c8Ho/69u2rt956q5kqRVOwGSMn45djXhxBHQAkNTz3ILZEVVB/6aWXNGPGDM2ZM0efffaZ+vXrp1GjRungwYM1tl+zZo0mTpyoKVOmaOPGjRo/frzGjx+vzz//vJkrR2OJCwQjXQKaQVyQ9xkAGpp7EHucDd1g8uTJmjJlin784x83RT11mj9/vqZOnaqrrrpKkrRo0SK9+eabeuqppzRz5sxq7f/yl79o9OjRuuWWWyRJd911l1asWKFHH31UixYtatba0TjigvS0tgR8IAPQXM4991yNGDFCc+bMCVt/+PBhXXzxxXrvvfciVFnDc09NbnpmtVzxiZKkqz/9f2HPLRp0hiTJV8Ydoa2qwT3qhYWFysnJUY8ePXTPPfdo7969TVFXNT6fTxs2bFBOTk5ond1uV05OjtauXVvjNmvXrg1rL0mjRo2qtT2sjwDXMtCjDqC5rF69Wo8++qjGjx+vkpL/Blafz6f3338/YnWdSO5B7GlwUF+2bJn27t2ra665Ri+99JK6dOmiCy+8UK+++qoqKiqaokZJUl5engKBgDIzM8PWZ2ZmKjc3t8ZtcnNzG9Rekrxer4qKisIWWIeTANci8D4DaE7vvvuucnNzdfbZZ2vnzp2RLkdSw3NPbfnF7XTI88PiT4gPWyrXu52OZtknNNwJjVHPyMjQjBkz9J///Efr1q3TKaecossvv1zZ2dm66aabtG3btsaus9nMmzdPqampoaVjx46RLgnHCNij6rIKnKCAzRbpEgC0IO3atdP777+vvn376qyzztLq1asjXVKD1ZZfTstooz6ZGTo9K0MHfzxYeWcPUN6Q/sob0l99MjPUJzNDp2W0iXD1qM1JpZ79+/drxYoVWrFihRwOh8aMGaPNmzerd+/eeuihhxqrRklSmzZt5HA4dODAgbD1Bw4cUFZWVo3bZGVlNai9JM2aNUuFhYWhZffu3SdfPBqNz0FQbwl8Dnp3ADQP2w8dA263W88//7x+//vfa/To0XrsscciWldDc0/t+cUoGAzKGCkQ55QxRvL7Jb9fwWBQwWBQErOpWVWDU09FRYX+53/+Rz/96U/VuXNnvfLKK7rxxhu1b98+/f3vf9e7776rl19+WXPnzm3UQl0ulwYOHKiVK1eG1gWDQa1cuVJDhw6tcZuhQ4eGtZekFStW1NpeOvofNSUlJWyBdVTYCXAtQQUfyAA0E1Nlyt/bb79dzz33nB588MEIVXRUQ3NPbfmlY3qKume2UteMVMVnpiuuewc5e3aRs2cXdc9spe6ZrdQxnaxjVQ2e9aVdu3YKBoOaOHGi1q9fr/79+1drM3LkSKWlpTVCeeFmzJihyZMna9CgQRo8eLAWLFigkpKS0NXQV1xxhdq3b6958+ZJkn7/+99rxIgRevDBBzV27Fi9+OKL+vTTT7V48eJGrw3NgwDXMlQwxAlAM9mxY4cyMjLC1l188cXq1auXPv300whVddTxck99JCfEKTXRJX8gqJJyvxI9Tjmrnkv9cY1cORpLg4P6Qw89pAkTJsjj8dTaJi0tTTt27Dipwmpy6aWX6tChQ5o9e7Zyc3PVv39/LV++PHShxa5du2Q/5gQ/bNgwPf/887r99tv1f/7P/1GPHj20bNkynX766Y1eG5qHsdnkt9m46VGMq2DoC4Bm0rlz5xrX9+nTR3369GnmasIdL/fUh8NhDwvmziqPK9vAmhoc1C+//PKmqKPepk+frunTp9f4XE0Xf0yYMEETJkxo4qrQnLxOh5wV/kiXgSbEtQgAcFRduQexj7Mhok6xyxXpEtDEitzuSJcAAEDEEdQRdQo8hLhYVmG3qSSuwV/2AQAQcwjqiDqFBPWYVuh2S8yjDgAAQR3Rhx712Mb7CwDAUXy/jKhT6GaMeizjGxMAaDxFJT4Zh1fB4NHZ0gpLfLLbw7+1LC7xRaI01ANBHVGn0ONWwGaTgykaY9JhgjoANBqX0y630y5/0KjCH1Sc0y5nlaDudTLAwqoI6og6AbtdBxPj1e5IaaRLQRPYn5wU6RIAIGZ43E4leOLkDwRV6PcpwV39hkc+L3HQqvgIhai0jzAXkwrcLpW4uEMeAAASQR1Ral9yYqRLQBPgAxgAAP9FUEdUOpCYqABT+MUcPoABAPBfBHVEJb/j6Dh1xBZ61AEA+C+COqLW7pTkSJeARpQX71Ep49MBAAghqCNqbWvdKtIloBF9k54W6RIAIOaUe/0qLa9QqdcvSSqtfHzMUv7Dc7Ae5uNB1Cp2u7Q/KYFpGmOAkbStdVqkywCAmOPzB+X1B0M3PKrwBxWoMo+6zx+MRGmoB4I6otq29FYE9RiwLzlRJS7uOAsAjS0l0aW0JPfRedRLfEpNdFWbR90W4PhrVQx9QVTbnp6qAJO/RL1t6QxjAgCgKoI6oprX6dSu1JRIl4GT4LfZ9G2r1EiXAQCA5RDUEfU2t20T6RJwEr5q00o+pyPSZQAAYDkEdUS9fcmJOsCc6lEpKOk/mRmRLgMAAEsiqCP62WzamNU20lXgBGxPT1ORxx3pMgAAsCRmfUFM2JmWonyPW+nl3kiXggbYmEVvOgA0pUAgKP8Pi6TQn1XbwJoI6ogNNps2ZbXVuTt3R7oS1NN3qcnKT2DIEgA0peLSCsnpCz0uKa9+c6Pi0ormLAkNQFBHzPgmPU0Dcg+qFb3qlmckfZqdGekyACDmJSfEKTXRJX8gqJJyvxI9zmrzqMsfF5nicFyMUUfMCNpt+rBTdqTLQD18kZGuQ4kJkS4DAGKew2GX84dFUujvxy6OqsEdlsE7g5iyNyVZ3zAnt6WVOxxa3z4r0mUAAGB5BHXEnLUd26nCzj9tq1rXIUteJ6PuAAA4HtIMYk6Jy6UN7Ziu0YoOJMbryzbpkS4DAICoQFBHTPp/mW2UF++JdBk4RsBm0wedO0g2W6RLAQAgKhDUEZOCdrtWdO/MEBgLWduhnb5nOkYAAOqNgaKIWYUetz7o3F7n7WBu9UjbkZaiz9u2jnQZANDiFJX4ZBxeBYNGklRY4pPdHv7NZnGJr6ZNYQF0NyKmbWvdSltbt4p0GS1asStOq7sw5AUAIsHltMvttCvOeTTyxf3w+NjF5SQOWhU96oh5H3Zqr7YlpUrnRkjNLijp3W6dmOUFACLE43YqwRMnfyCoQr9PCe7qNzzyeTlGWxUfoRDz/A673j6ls8qcjkiX0uJ82Km9DiQlRroMAACiEkEdLUKhx6N/ndJVFXaGXzSXT9u11ReMSwcA4IQR1NFiHExK0DvduyhAVm9yWzLS9Wl2ZqTLAAAgqhHU0aLsTk3W6i4dI11GTPs2LUUfdmrPxaMAAJwkgjpanG2tW+mjjtmRLiMm7UlO0spunWQI6QAAnDQu80WLtDmzjXwOu0bs3MOn1UayvVWqVnbtqCA3mQIAyyj3+lVaXiH/D/Ool3r9cla5Xqvc649EaagHgjparK/apKvM6dT5336nuB8OYDgxn2e01kedsulJBwCL8fmD8vqDoRseVfiDClQJ6j5/MBKloR4I6mjRdqWl6I1Tu+nCbTvlCQQiXU5UWp+dqc/atWVMOgBYUEqiS2lJ7qPzqJf4lJroqjaPui3gilB1OB6+o0aLdyApUf/s1V1FLg5UDRGw2bS6c3t9lp1JSAcAoAkQ1AFJh+M9erV3D21vlRrpUqJCodul13p119YM5kkHgKbQpUsX2Wy2sOXee++NdFloZgx9AX7gczq0olsn7cnL14927ZPTMG69JtvSU/VB5w6qcHCnVwBoSnPnztXUqVNDj5OTkyNYDSKBoA4cy2bTlxmtdTAxQTnf7lKrcm+kK7IMv82mDztla2ubdIa6AEAzSE5OVlZWVqTLQAQx9AWowfcJ8fqf007RxqwM7mSqo/Ojv9q7x9GhLoR0AGgW9957r1q3bq0BAwboz3/+s/z+2qdR9Hq9KioqClsQ/ehRB2rhdzi0rkM7fdW6lc7ZtU8dio9EuqRmVxLn1JqO2UfH7hPQAaDZ3HDDDTrzzDOVnp6uNWvWaNasWdq/f7/mz59fY/t58+bpzjvvrLY+EAjK/8MiKfRn1TawJoI6cBwF8R69cWpXdT9cqKG79yupoiLSJTW5gE3a3DZDG7LbMhYdABrJzJkzdd9999XZ5ssvv1SvXr00Y8aM0LozzjhDLpdLv/vd7zRv3jy53e5q282aNStsm6KiInXs2FHFpRWS0xdaX1JevVe+uDT2z2vRiqAO1IfNpu3pafouNVmnH/xeZxzIU0IdX0FGq6Ckba1b6bN2GSr0eCJdDgDElD/84Q+68sor62zTrVu3GtcPGTJEfr9fO3fuVM+ePas973a7awzwyQlxSk10yR8IqqTcr0SPs9o86vLH1Xsf0LwI6kAD+B0ObWrXVpsz26hXXr765x5Ssi/6eyL8Npu+zEjXfzIzdMTNfPIA0BQyMjKUkZFxQttu2rRJdrtdbdu2bdB2Doc9LJg7qzyubANrIqgDJyBgt2tL2zb6sk1rdT9coH65h9SmrDzSZTVYmdOhL9u01v/LbKPyOA4HAGAFa9eu1bp16zRy5EglJydr7dq1uummm/TrX/9arVq1inR5aEZR8xEqPz9fkyZNUkpKitLS0jRlyhQdOVL7xX35+fm6/vrr1bNnT8XHx6tTp0664YYbVFhY2IxVI9YF7TZta91Kr/Y5VS/37qGNWRk6EmftrxD9Npu2pafqrVO66Jkzemt9hyxCOgBYiNvt1osvvqgRI0aoT58+uvvuu3XTTTdp8eLFkS4NzSxqzs6TJk3S/v37tWLFClVUVOiqq67StGnT9Pzzz9fYft++fdq3b58eeOAB9e7dW999952uvvpq7du3T6+++mozV4+WID8hXusS4rWufZayi0vUI/+wOhcUW2Isu99m0/7kRG1Lb6UdrVK4QBQALOzMM8/Uxx9/HOkyYAFREdS//PJLLV++XJ988okGDRokSXrkkUc0ZswYPfDAA8rOzq62zemnn67/+Z//CT3u3r277r77bv3617+W3++X0xkVu45oZLNpX0qS9qUkScYordyr7OISZRcfUXZxSbMEd7/NpgNJCdqXnKR9yYk6mJiggD1qvkADAACKkqC+du1apaWlhUK6JOXk5Mhut2vdunX6+c9/Xq/XKSwsVEpKSp0h3ev1yuv9790ouWEATorNpoJ4jwriPfqibetQcG9dVq7Ucq/Syr1K9R79030C89gGbDYVul0q9LhV4HGr0H30z0OJ8QRzAACiXFQE9dzc3GpXOTudTqWnpys3N7der5GXl6e77rpL06ZNq7NdbTcMABrFMcE9jDHy+ANyB/yKCwTlCgQVFwwqLhBQXDAov92uCrtdFY7KPx2qsNtV4oqT4UZEAIBaFJX4ZBxeBYNGklRY4pPdHn7eKC7x1bQpLCCiQb2+E/+frKKiIo0dO1a9e/fWHXfcUWfb2m4YADQpm03lcU4u6gQANCqX0y630y5/0KjCH1Sc0y5nlaDudfINrFVFNBXUd+L/rKwsHTx4MGy93+9Xfn6+srKy6ty+uLhYo0ePVnJysl577TXFHWdGjtpuGAAAABBtPG6nEjxx8geCKvT7lOCufsMjn5dOIquK6DtT34n/hw4dqoKCAm3YsEEDBw6UJL333nsKBoMaMmRIrdsVFRVp1KhRcrvdev311+XhTosAAACIElHxXcdpp52m0aNHa+rUqVq/fr0++ugjTZ8+XZdddlloxpe9e/eqV69eWr9+vaSjIf2CCy5QSUmJnnzySRUVFSk3N1e5ubkKBAKR3B0AAADguKLmu47nnntO06dP13nnnSe73a6LL75YDz/8cOj5iooKffXVVyotLZUkffbZZ1q3bp0k6ZRTTgl7rR07dqhLly7NVjsAAADQUFET1NPT02u9uZEkdenSRcaY0OOf/OQnYY8BAACAaBIVQ18AAACAliZqetQBAADQMOVev0rLK+T/YR71Uq+/2vSM5d6mv2M2TgxBHQAAIEb5/EF5/cHQDY8q/EEFqgR1n7/hd8ZG8yCoAwAAxKiURJfSktxH51Ev8Sk10VVtHnVbwBWh6nA8jFEHAAAALIigDgAAAFgQQR0AAACwIII6AAAAYEEEdQAAAMCCmPUFAAAgRgUCQfl/WCSF/qzaBtZEUAcAAIhRxaUVktMXelxSXv3mRsWlFc1ZEhqAoA4AABCjkhPilJrokj8QVEm5X4keZ7V51OWPi0xxOC6COgAAQIxyOOxhwdxZ5XFlG1gT7wwAAABgQQR1AAAAwIII6gAAAIAFEdQBAAAACyKoAwAAABbErC8AAAAxqqjEJ+PwKhg0kqTCEp/sdltYm+ISX02bwgII6gAAADHK5bTL7bTLHzSq8AcV57TLWSWoe50MsLAqgjoAAECM8ridSvDEyR8IqtDvU4K7+g2PfF7ioFXxEQoAAACwIII6AAAAYEEEdQAAAMCCCOoAAACABRHUAQAAAAviMl8AAIAYVe71q7S8Qv4f5lEv9fqrTc9Y7vVHojTUA0EdAAAgRvn8QXn9wdANjyr8QQWqBHWfPxiJ0lAPBHUAAIAYlZLoUlqS++g86iU+pSa6qs2jbgu4IlQdjocx6gAAAIAFEdQBAAAACyKoAwAANKO7775bw4YNU0JCgtLS0mpss2vXLo0dO1YJCQlq27atbrnlFvn9XPTZ0jBGHQAAoBn5fD5NmDBBQ4cO1ZNPPlnt+UAgoLFjxyorK0tr1qzR/v37dcUVVyguLk733HNPBCpGpNCjDgAA0IzuvPNO3XTTTerbt2+Nz7/zzjv64osv9Oyzz6p///668MILddddd2nhwoXy+XzNXC0iiaAOAABgIWvXrlXfvn2VmZkZWjdq1CgVFRVpy5YtDXqtQCAo/w+LpNDfj10CAaZntCqGvgAAAFhIbm5uWEiXFHqcm5tb4zZer1derzf0uKioSJJUXFohOf/bC19SXn2ce3FpxUnXjKZBjzoAAMBJmjlzpmw2W53L1q1bm+znz5s3T6mpqaGlY8eOkqTkhDilJrqU6DnaN5vocSo10RW2JCfENVldODn0qAMAAJykP/zhD7ryyivrbNOtW7d6vVZWVpbWr18ftu7AgQOh52oya9YszZgxI/S4qKhIHTt2lLf0iLxxdvkDQZWVVcgZjFOgyg2PvKVH6lUXmh9BHQAA4CRlZGQoIyOjUV5r6NChuvvuu3Xw4EG1bdtWkrRixQqlpKSod+/eNW7jdrvldrurrf/PupWKj4+XMVIgaOSw22SzHX3OBIPy+ytUXu6tth2sgaAOAADQjHbt2qX8/Hzt2rVLgUBAmzZtkiSdcsopSkpK0gUXXKDevXvr8ssv1/3336/c3Fzdfvvtuu6662oM43Wx2x1yutwyxkgBI6fj6DCcQIVfZeVHZHc4FedyNcFeojEQ1AEAAJrR7Nmz9fe//z30eMCAAZKkVatW6Sc/+YkcDofeeOMNXXPNNRo6dKgSExM1efJkzZ07t8E/y+ZwyOmMkzFGxhaU02FXwO9XWVmxnHEuJaWkqbi4uNH2DY2LoA4AANCMli5dqqVLl9bZpnPnznrrrbca/Wf7/RUqKS6Qw+FUUkqabDbmFbEygjoAAEALEPBXyFtSJIeTkB4tCOoAAAAxKhAwqvAH5PdXqOxIoewOp9wJyfIHjKRAqA2siaAOAAAQo+w2yQT8Ki8pkt3hVGJyquxVetLttggVh+MiqAMAAMQoEwyotKRQTodTrsRkOR0O2WzhydxGUrcsBicBAADEqLKSYjkcTiWmpDImPQrRow4AABCjbA6HklLSJNkkBSNcDRqKj1YAAACABRHUAQAAYpQJBHSkqEDG0JsejaImqOfn52vSpElKSUlRWlqapkyZoiNHjtRrW2OMLrzwQtlsNi1btqxpCwUAALCI+MRkBQJ+lRQVEtajUNQE9UmTJmnLli1asWKF3njjDX3wwQeaNm1avbZdsGBBtSucAQAAYp3N7lBCYqr8Ab/KjhTJHwgoEAiGLSbIPOpWFRUXk3755Zdavny5PvnkEw0aNEiS9Mgjj2jMmDF64IEHlJ2dXeu2mzZt0oMPPqhPP/1U7dq1a66SAQAAIi5oJJvDKU9iisqOFKqkuFAJSSlhM8CQ060rKnrU165dq7S0tFBIl6ScnBzZ7XatW7eu1u1KS0v1q1/9SgsXLlRWVla9fpbX61VRUVHYAgAAEI0cDpvinA553G7FJ6VKwYC8pcVy/rA+zumQw8GoA6uKiqCem5urtm3bhq1zOp1KT09Xbm5urdvddNNNGjZsmMaNG1fvnzVv3jylpqaGlo4dO55w3QAAAFbhcMYpMSVNgYCfC0yjRESD+syZM2Wz2epctm7dekKv/frrr+u9997TggULGrTdrFmzVFhYGFp27959Qj8fAADAapzOOCUltyKsR4mIjlH/wx/+oCuvvLLONt26dVNWVpYOHjwYtt7v9ys/P7/WIS3vvfeetm/frrS0tLD1F198sYYPH67Vq1fXuJ3b7Zbb7a7vLgAAAEQVZ9zRsH6k+LCOFBUoSgZYtEgRDeoZGRnKyMg4bruhQ4eqoKBAGzZs0MCBAyUdDeLBYFBDhgypcZuZM2fqt7/9bdi6vn376qGHHtJFF1108sUDAABEqWPDurfcF+lyUIuomPXltNNO0+jRozV16lQtWrRIFRUVmj59ui677LLQjC979+7Veeedp6efflqDBw9WVlZWjb3tnTp1UteuXZt7FwAAACylMqyXleyNdCmoRdR81/Hcc8+pV69eOu+88zRmzBidc845Wrx4cej5iooKffXVVyotLY1glQAAABZijMwPy9GHJmxxOJ3yJCZFuEjUJip61CUpPT1dzz//fK3Pd+nSJfSPsDbHex4AACCWBAJSReC/F4z6g0ZS1TwUNXGwxeGdAQAAiFEOhxTnsMsYI3/QyGm3Vbtbu8MRoeJwXAR1AACAWGU7Npib0PTXVdvAmqJmjDoAAADQkhDUAQAAAAsiqAMAAAAWRFAHAAAALIigDgAAAFgQs74AAADEqEDAqMIfCM2cXhEIquocL4EA95mxKoI6AABAjLLbJLvNJiMpaIzsNlu1oG6rdgMkWAVBHQAAIEbZ7DY5frjhUTBg5KhywyNjgvKWlUSwQtSFMeoAAAAtkDFBHSkqUDDgj3QpqAVBHQAAoIWpDOmBgF/xicmRLge1IKgDAAC0IMeG9KTkVrI7GQltVQR1AACAFqJqSHfGxUW6JNSBj1AAAAAtgDFBlRQVKxAkpEcLgjoAAECMMkGjQCCooAmq7EiRTDCgxKRU2ewOBQLBUBtYE0EdAAAgRgWNFAgGVHqkSMGAX/FJqbI5nAoaE9YG1kRQBwAAiFEOu+QtLZaCAcUnpcrjdofNoy5JDkfVWyDBKriYFAAAIEaVlRxRIOBXYkqaHE7GpEcbgjoAAECMClbO7kJIj0oEdQAAgBjliHMyu0sUI6gDAADEKH9FhcpLSyJdBk4QF5MCAADEKLcnQWVlR2Rk5HDFR7ocNBBB/TjMD9MXbc/fXuPzlesrvGXylfGJFUB0qfCWSar9GHfsc8YwhxsQLSr/v/oq/IpzOpT//UHFuRLkcXukKpO8lJQcCdsG1mEzvCt12rNnjzp27BjpMgAg4nbv3q0OHTpEugwA9fDtt9+qe/fuDdpm+/bt6tatWxNVhBNBj/pxZGdna/fu3UpOTq4272ikFRUVqWPHjtq9e7dSUlIiXU6zYJ/Z51hk9f01xqi4uFjZ2dmRLgVAPaWnp0uSdu3apdTU1ND6mo43hYWF6tSpU2gbWAdB/Tjsdrvle5BSUlIseXJvSuxzy9DS9tnK+3vsiR6A9dntR+cLSU1NrfG4UtPxpnIbWAfvCAAAAGBBBHUAAADAggjqUcztdmvOnDlyu92RLqXZsM8tQ0vb55a2vwCaXm3HlZrWcwyyLmZ9AQAAACyIHnUAAADAggjqAAAAgAUR1AEAAAALIqhHmfz8fE2aNEkpKSlKS0vTlClTdOTIkXpta4zRhRdeKJvNpmXLljVtoY2oofucn5+v66+/Xj179lR8fLw6deqkG264QYWFhc1YdcMsXLhQXbp0kcfj0ZAhQ7R+/fo627/yyivq1auXPB6P+vbtq7feequZKm0cDdnfJUuWaPjw4WrVqpVatWqlnJyc4/5+rKih73GlF198UTabTePHj2/aAgFEvcrzpcfjkd1ul9Pp1FlnnaX169fXegz6yU9+IpvNFracccYZ1dZ5PJ4I713LRFCPMpMmTdKWLVu0YsUKvfHGG/rggw80bdq0em27YMECy91dtT4aus/79u3Tvn379MADD+jzzz/X0qVLtXz5ck2ZMqUZq66/l156STNmzNCcOXP02WefqV+/fho1apQOHjxYY/s1a9Zo4sSJmjJlijZu3Kjx48dr/Pjx+vzzz5u58hPT0P1dvXq1Jk6cqFWrVmnt2rXq2LGjLrjgAu3du7eZKz9xDd3nSjt37tTNN9+s4cOHN1OlAKLZpEmT9OGHHyoYDOrWW29V+/btdfjwYY0cObLWY1BhYaFsNpv+9Kc/6f3339eNN96oL774QklJSdq/f39o+e677yK9ey2TQdT44osvjCTzySefhNb961//Mjabzezdu7fObTdu3Gjat29v9u/fbySZ1157rYmrbRwns8/Hevnll43L5TIVFRVNUeZJGTx4sLnuuutCjwOBgMnOzjbz5s2rsf0vf/lLM3bs2LB1Q4YMMb/73e+atM7G0tD9rcrv95vk5GTz97//valKbHQnss9+v98MGzbM/PWvfzWTJ08248aNa4ZKAUSryvNlnz59Qsebf/3rX0aSiYuLM2effXao7bHHoIyMDNOlS5ew1+rWrZtxuVzNWj9qRo96FFm7dq3S0tI0aNCg0LqcnBzZ7XatW7eu1u1KS0v1q1/9SgsXLlRWVlZzlNpoTnSfqyosLFRKSoqcTmdTlHnCfD6fNmzYoJycnNA6u92unJwcrV27tsZt1q5dG9ZekkaNGlVreys5kf2tqrS0VBUVFUpPT2+qMhvVie7z3Llz1bZtW8t+EwTAWtauXavU1FRt3bo1dLypPF9WVFQoEAiE2h57DCoqKtKhQ4fUpk0bnX766Zo1a5ZOO+00+Xw+de7cWR07dtS4ceO0ZcuWSO1ai0ZQjyK5ublq27Zt2Dqn06n09HTl5ubWut1NN92kYcOGady4cU1dYqM70X0+Vl5enu666656DxFqTnl5eQoEAsrMzAxbn5mZWev+5ebmNqi9lZzI/lZ16623Kjs7u9qHFas6kX3+8MMP9eSTT2rJkiXNUSKAGJCbm6vWrVuHHW+cTqfS0tIkqdq1XZXHIL/fr+nTp2vVqlWaNWuWnnnmGX311VdKTk7WP//5Tz377LMKBoMaNmyY9uzZ09y71eIR1C1g5syZ1S7aqLps3br1hF779ddf13vvvacFCxY0btEnqSn3+VhFRUUaO3asevfurTvuuOPkC0dE3XvvvXrxxRf12muvxeyFTcXFxbr88su1ZMkStWnTJtLlAIiwpj5f2u129evXT3379tWkSZP09NNP65tvvpHL5VL//v01YsQI/eMf/1BGRoaeeOKJRtwz1Ie1xgG0UH/4wx905ZVX1tmmW7duysrKqnbxmd/vV35+fq1DWt577z1t37499Im60sUXX6zhw4dr9erVJ1H5iWvKfa5UXFys0aNHKzk5Wa+99pri4uJOtuxG16ZNGzkcDh04cCBs/YEDB2rdv6ysrAa1t5IT2d9KDzzwgO699169++67OuOMM5qyzEbV0H3evn27du7cqYsuuii0LhgMSjraO/bVV1+pe/fuTVs0AMtoyPny+++/Dzve+P1+FRQUSJKSkpLCtqk8BlU9pwwZMkSSlJKSEloXFxenAQMG6JtvvmmEPUJDENQtICMjQxkZGcdtN3ToUBUUFGjDhg0aOHCgpKNBPBgMhv5jVTVz5kz99re/DVvXt29fPfTQQ2FBoLk15T5LR3vSR40aJbfbrddff92yva8ul0sDBw7UypUrQ9PvBYNBrVy5UtOnT69xm6FDh2rlypW68cYbQ+tWrFihoUOHNkPFJ+dE9leS7r//ft199916++23w65XiAYN3edevXpp8+bNYetuv/12FRcX6y9/+Ys6duzYHGUDsIiGnC8LCwvVp0+f0PGm8nwZFxcXdo3Wsccgj8cTdk7ZtGmTJIXOuZIUCAS0efNmjRkzplH3DfUQ6atZ0TCjR482AwYMMOvWrTMffvih6dGjh5k4cWLo+T179piePXuadevW1foaiqJZX4xp+D4XFhaaIUOGmL59+5pvvvnG7N+/P7T4/f5I7UatXnzxReN2u83SpUvNF198YaZNm2bS0tJMbm6uMcaYyy+/3MycOTPU/qOPPjJOp9M88MAD5ssvvzRz5swxcXFxZvPmzZHahQZp6P7ee++9xuVymVdffTXsvSwuLo7ULjRYQ/e5KmZ9AVAfo0ePNl26dDEul8vcdtttpnPnzqZ79+4mISHBuN1u89BDD5nExETTo0cPk5aWZj7++GPz29/+1thsNjN06FCzcOFCk5aWZiSZJ554wmzfvt1s2LDBXHbZZcbj8ZgtW7ZEehdbHIJ6lPn+++/NxIkTTVJSkklJSTFXXXVVWGDZsWOHkWRWrVpV62tEW1Bv6D6vWrXKSKpx2bFjR2R24jgeeeQR06lTJ+NyuczgwYPNxx9/HHpuxIgRZvLkyWHtX375ZXPqqacal8tl+vTpY958881mrvjkNGR/O3fuXON7OWfOnOYv/CQ09D0+FkEdQH1Uni9dLpex2WzGbrebQYMGmY8//tg88sgjJjs720gyqamp5uOPPza7du0yP/7xj43dbg8dW9u0aWPGjh0bOl5lZmaaMWPGmM8++yzSu9ci2Ywxppk78QEAAAAcB7O+AAAAABZEUAcAAAAsiKAOAAAAWBBBHQAAALAggjoAAABgQQR1AAAAwIII6gAAAIAFEdQBAAAACyKoAwAAABZEUAcAAAAsiKAOAAAAWBBBHbCgQ4cOKSsrS/fcc09o3Zo1a+RyubRy5coIVgYAiDZPP/20WrduLa/XG7Z+/PjxuvzyyyNUFerDZowxkS4CQHVvvfWWxo8frzVr1qhnz57q37+/xo0bp/nz50e6NABAFCkrK1O7du20ZMkSTZgwQZJ08OBBtW/fXu+8845GjhwZ4QpRG4I6YGHXXXed3n33XQ0aNEibN2/WJ598IrfbHemyAABR5tprr9XOnTv11ltvSZLmz5+vhQsX6ptvvpHNZotwdagNQR2wsLKyMp1++unavXu3NmzYoL59+0a6JABAFNq4caPOOussfffdd2rfvr3OOOMMTZgwQX/6058iXRrqwBh1wMK2b9+uffv2KRgMaufOnZEuBwAQpQYMGKB+/frp6aef1oYNG7RlyxZdeeWVkS4Lx0GPOmBRPp9PgwcPVv/+/dWzZ08tWLBAmzdvVtu2bSNdGgAgCj3++ONasGCBzj//fG3btk1vv/12pEvCcRDUAYu65ZZb9Oqrr+o///mPkpKSNGLECKWmpuqNN96IdGkAgChUWFio7Oxs+f1+Pf3007r00ksjXRKOg6EvgAWtXr1aCxYs0DPPPKOUlBTZ7XY988wz+ve//63HH3880uUBAKJQamqqLr74YiUlJWn8+PGRLgf1QI86AABAC3HeeeepT58+evjhhyNdCuqBoA4AABDjDh8+rNWrV+uSSy7RF198oZ49e0a6JNSDM9IFAAAAoGkNGDBAhw8f1n333UdIjyL0qAMAAAAWxMWkAAAAgAUR1AEAAAALIqgDAAAAFkRQBwAAACyIoA4AAABYEEEdAAAAsCCCOgAAAGBBBHUAAADAggjqAAAAgAX9f6t8NB4hg/O0AAAAAElFTkSuQmCC\n",
"text/plain": [
"
\n"
],
"text/plain": [
"\u001b[2;36m[23:03:16]\u001b[0m\u001b[2;36m \u001b[0mBatch complete. \u001b]8;id=75146;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/web/container.py\u001b\\\u001b[2mcontainer.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=833728;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/web/container.py#382\u001b\\\u001b[2m382\u001b[0m\u001b]8;;\u001b\\\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# initialize a batch and run them all\n",
"batch = td.web.Batch(simulations=sims)\n",
"\n",
"# run the batch and store all of the data in the `data/` dir.\n",
"batch_data = batch.run(path_dir=\"data\")\n"
]
},
{
"cell_type": "markdown",
"id": "c865c155",
"metadata": {},
"source": [
"Now that the simulations are complete, we can analyze the data. Let's first look at one of the FieldTimeMonitors to make sure the source has decayed."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "4a97628f-c2d1-4e2f-82f2-ba97d88fe48e",
"metadata": {
"execution": {
"iopub.execute_input": "2023-02-03T05:03:17.417842Z",
"iopub.status.busy": "2023-02-03T05:03:17.417763Z",
"iopub.status.idle": "2023-02-03T05:03:19.441657Z",
"shell.execute_reply": "2023-02-03T05:03:19.441405Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
[23:03:17] INFO downloading file \"output/monitor_data.hdf5\" to webapi.py:593\n",
"\"data/524d39cf-9f84-4f05-9d5c-f8d13bd6b5c2.hdf5\"\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(\n",
" batch_data[\"sim_1\"].monitor_data[\"monitor_time_0\"].Hz.t,\n",
" np.real(batch_data[\"sim_1\"].monitor_data[\"monitor_time_0\"].Hz.squeeze()),\n",
")\n",
"plt.title(\"FieldTimeMonitor data\")\n",
"plt.xlabel(\"t\")\n",
"plt.ylabel(\"Hz\")\n"
]
},
{
"cell_type": "markdown",
"id": "21bc59dc-0163-432b-9cf0-5c2d2617ef35",
"metadata": {},
"source": [
"We see that the source has mostly decayed by the time we switch on the monitors, and the remaining data shows decay and oscillation due to the resonances inside the system.\n",
"\n",
"Looking at the Fourier transform of this data, we can see resonances at the band frequencies."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "e1b95adb-1ea6-40cb-ba63-ef5ac60de470",
"metadata": {
"execution": {
"iopub.execute_input": "2023-02-03T05:03:19.451411Z",
"iopub.status.busy": "2023-02-03T05:03:19.451337Z",
"iopub.status.idle": "2023-02-03T05:03:19.532499Z",
"shell.execute_reply": "2023-02-03T05:03:19.532252Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
INFO loading SimulationData from data/524d39cf-9f84-4f05-9d5c-f8d13bd6b5c2.hdf5 webapi.py:415\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"df = 1 / np.amax(batch_data[\"sim_1\"].monitor_data[\"monitor_time_0\"].Hz.t)\n",
"minn = int(freq_range[0] / df)\n",
"maxn = int(freq_range[1] / df)\n",
"spectrum = np.fft.fft(batch_data[\"sim_1\"].monitor_data[\"monitor_time_0\"].Hz.squeeze())\n",
"plt.plot(\n",
" np.linspace(freq_range[0], freq_range[1], maxn - minn),\n",
" np.abs(spectrum[::-1][minn:maxn]),\n",
")\n",
"plt.title(\"Spectrum at single wavevector\")\n",
"plt.xlabel(\"Frequency (Hz)\")\n",
"plt.ylabel(\"Amplitude\")\n"
]
},
{
"cell_type": "markdown",
"id": "663c96ba-8b16-4056-9983-bef38e7a5c9a",
"metadata": {},
"source": [
"We use the ResonanceFinder plugin to find the band frequencies.\n",
"\n",
"We first construct a `ResonanceFinder` object storing our parameters, and then call `run()` on our list of `FieldTimeData` objects. This will add up the signals from all of the monitors before searching for resonances. The `ResonanceFinder` class has additional methods in case the signal takes another form; see the api reference [here](https://docs.flexcompute.com/projects/tidy3d/en/v1.9.0rc2/_autosummary/tidy3d.plugins.ResonanceFinder.html).\n",
"\n",
"The `run()` method returns an `xr.Dataset` containing the decay rate, Q factor, amplitude, phase, and estimation error for each resonance as a function of frequency. "
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "7496ef10-068e-41ce-8907-418bbcd60e26",
"metadata": {
"execution": {
"iopub.execute_input": "2023-02-03T05:03:19.534085Z",
"iopub.status.busy": "2023-02-03T05:03:19.534000Z",
"iopub.status.idle": "2023-02-03T05:03:20.324743Z",
"shell.execute_reply": "2023-02-03T05:03:20.324230Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
INFO loading SimulationData from data/524d39cf-9f84-4f05-9d5c-f8d13bd6b5c2.hdf5 webapi.py:415\n",
"
"
],
"text/plain": [
" decay Q amplitude phase error\n",
"freq \n",
"9.788253e+13 7.430813e+11 413.826865 0.008332 1.656388 0.000041\n",
"1.061235e+14 9.529783e+10 3498.472629 0.037352 3.054479 0.000009\n",
"1.131265e+14 2.128035e+12 167.007322 0.038965 1.695435 0.000051\n",
"1.191562e+14 1.150597e+12 325.344332 0.022350 -2.133775 0.000034\n",
"1.372927e+14 1.121520e+12 384.583068 0.031599 -0.707254 0.000053"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"resonance_finder = ResonanceFinder(freq_window=tuple(freq_range))\n",
"resonance_data = resonance_finder.run(signals=batch_data[\"sim_1\"].data)\n",
"resonance_data.to_dataframe()\n"
]
},
{
"cell_type": "markdown",
"id": "46db3d3c",
"metadata": {},
"source": [
"We see the four resonances from the previous figure. All four have reasonable Q factors, amplitudes, and errors, so they are likely to represent physical resonances. Note that in order to accurately obtain the Q factor for high-Q modes, it may be necessary to run the simulation for a longer time.\n",
"\n",
"Now we are ready to compute the band structure. We run the resonance finder at every Bloch wavevector."
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "4afe673e",
"metadata": {
"execution": {
"iopub.execute_input": "2023-02-03T05:03:20.327892Z",
"iopub.status.busy": "2023-02-03T05:03:20.327655Z",
"iopub.status.idle": "2023-02-03T05:03:50.334645Z",
"shell.execute_reply": "2023-02-03T05:03:50.334018Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
INFO downloading file \"output/monitor_data.hdf5\" to webapi.py:593\n",
"\"data/425fbe84-a204-43ca-a331-08a2bd603db7.hdf5\"\n",
"
sim_0: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_1: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_2: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_3: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_4: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_5: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_6: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_7: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_8: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_9: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_10: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_11: status = postprocess ━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸━━━━━━━━━━━ 71%0:00:23\nsim_0: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_1: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_2: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_3: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_4: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_5: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_6: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_7: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_8: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_9: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_10: status = success ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━100%0:00:00\nsim_11: status = postprocess ━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸━━━━━━━━━━━ 71%0:00:23
\n",
"text/plain": "\r\u001b[2K\u001b[1A\u001b[2K\u001b[1A\u001b[2K\u001b[1A\u001b[2K\u001b[1A\u001b[2K\u001b[1A\u001b[2K\u001b[1A\u001b[2K\u001b[1A\u001b[2K\u001b[1A\u001b[2K\u001b[1A\u001b[2K\u001b[1A\u001b[2K\u001b[1A\u001b[2Ksim_0: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_1: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_2: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_3: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_4: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_5: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_6: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_7: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_8: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_9: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_10: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_11: status = postprocess \u001b[38;2;249;38;114m━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[38;2;249;38;114m╸\u001b[0m\u001b[38;5;237m━━━━━━━━━━━\u001b[0m \u001b[35m 71%\u001b[0m \u001b[36m0:00:23\u001b[0m\nsim_0: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_1: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_2: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_3: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_4: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_5: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_6: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_7: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_8: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_9: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_10: status = success \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[35m100%\u001b[0m \u001b[36m0:00:00\u001b[0m\nsim_11: status = postprocess \u001b[38;2;249;38;114m━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[38;2;249;38;114m╸\u001b[0m\u001b[38;5;237m━━━━━━━━━━━\u001b[0m \u001b[35m 71%\u001b[0m \u001b[36m0:00:23\u001b[0m"
},
"metadata": {},
"output_type": "display_data"
}
],
"tabbable": null,
"tooltip": null
}
},
"dac5f20ef7b04279b466cfcd21d52610": {
"model_module": "@jupyter-widgets/output",
"model_module_version": "1.0.0",
"model_name": "OutputModel",
"state": {
"_dom_classes": [],
"_model_module": "@jupyter-widgets/output",
"_model_module_version": "1.0.0",
"_model_name": "OutputModel",
"_view_count": null,
"_view_module": "@jupyter-widgets/output",
"_view_module_version": "1.0.0",
"_view_name": "OutputView",
"layout": "IPY_MODEL_ac96117d4a39453fb1e6e42f45f83cc7",
"msg_id": "",
"outputs": [
{
"data": {
"text/html": "