{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tidy3D FDTD accuracy benchmark\n",
    "\n",
    "Tidy3D offers a lightning-fast FDTD solver capable of completing simulations in a fraction of the time required by traditional FDTD solvers, reducing runtimes from hours or days to mere minutes. This is achieved through cutting-edge integration between our hardware and state-of-the-art numerical methods, which involves no approximations. Therefore, this enhanced speed does not compromise accuracy, as reliable results remain paramount to the utility of any simulation tool.\n",
    "\n",
    "In order to investigate Tidy3D speed advantage and confirm its robustness, we will compare the results of a MMI power divider, available in our [application library](https://tidy3d.simulation.cloud/workbench?taskId=pa-5e58d0d7-c3dc-446c-9309-fab8b94da423), with those from the open-source software Meep."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import tidy3d as td\n",
    "from tidy3d import web"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FDTD Simulation\n",
    "\n",
    "The model is based on our web GUI [example](https://tidy3d.simulation.cloud/workbench?taskId=pa-5e58d0d7-c3dc-446c-9309-fab8b94da423). Below is the code exported directly from the GUI with minor modifications, including added comments and changes to variable names for improved readability:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# defining parameters\n",
    "si_index = 3.48\n",
    "sio2_index = 1.46\n",
    "central_wavelength = 1.55\n",
    "waveguide_width = 0.5\n",
    "waveguide_thickness = 0.22\n",
    "waveguide_length = 2\n",
    "mmi_width = 3.6\n",
    "mmi_length = 11.5\n",
    "output_separation = 0.6\n",
    "taper_width = 1.2\n",
    "taper_length = 20.05\n",
    "SiO2_medium = td.Medium(\n",
    "    name=\"SiO2\",\n",
    "    permittivity=2.1316,\n",
    ")\n",
    "\n",
    "# operating frequency\n",
    "freq0 = 193.414e12\n",
    "fwidth = 6.24568e12\n",
    "\n",
    "freqs = np.linspace(freq0 - fwidth, freq0 + fwidth, 50)\n",
    "\n",
    "# mode source\n",
    "Mode_source = td.ModeSource(\n",
    "    name=\"Mode source\",\n",
    "    num_freqs=10,\n",
    "    center=[-mmi_length / 2 - taper_length - waveguide_length / 2, 0, 0],\n",
    "    size=[0, 3 * waveguide_width, 4 * waveguide_thickness],\n",
    "    source_time=td.GaussianPulse(\n",
    "        freq0=193414489032258.06,\n",
    "        fwidth=6245676208333.328,\n",
    "    ),\n",
    "    direction=\"+\",\n",
    "    mode_spec=td.ModeSpec(\n",
    "        target_neff=3.48,\n",
    "    ),\n",
    ")\n",
    "\n",
    "# defining monitors for calculating transmittance and field visualization\n",
    "\n",
    "Flux_monitor = td.FluxMonitor(\n",
    "    name=\"Flux monitor\",\n",
    "    center=[\n",
    "        mmi_length / 2 + taper_length + waveguide_length / 2,\n",
    "        (output_separation + taper_width) / 2,\n",
    "        0,\n",
    "    ],\n",
    "    size=[0, 3 * waveguide_width, 4 * waveguide_thickness],\n",
    "    freqs=freqs,\n",
    "    normal_dir=\"+\",\n",
    ")\n",
    "\n",
    "Flux_monitorRef = td.FluxMonitor(\n",
    "    name=\"Ref monitor\",\n",
    "    center=(Mode_source.center[0] + 0.5, Mode_source.center[1], Mode_source.center[2]),\n",
    "    size=Mode_source.size,\n",
    "    freqs=freqs,\n",
    "    normal_dir=\"+\",\n",
    ")\n",
    "\n",
    "\n",
    "Mode_monitor = td.ModeMonitor(\n",
    "    name=\"Mode monitor\",\n",
    "    center=[\n",
    "        mmi_length / 2 + taper_length + waveguide_length / 2,\n",
    "        (output_separation + taper_width) / 2,\n",
    "        0,\n",
    "    ],\n",
    "    size=[0, 3 * waveguide_width, 4 * waveguide_thickness],\n",
    "    freqs=freqs,\n",
    "    mode_spec=td.ModeSpec(\n",
    "        num_modes=3,\n",
    "        target_neff=3.48,\n",
    "    ),\n",
    ")\n",
    "\n",
    "Field_monitor = td.FieldMonitor(\n",
    "    name=\"Field monitor\", size=[td.inf, td.inf, 0], freqs=[freq0], colocate=False\n",
    ")\n",
    "\n",
    "# defining mediums\n",
    "Si_medium = td.Medium(\n",
    "    name=\"Si\",\n",
    "    permittivity=12.1104,\n",
    ")\n",
    "\n",
    "# defining structures\n",
    "MMI_section = td.Structure(\n",
    "    name=\"MMI section\",\n",
    "    geometry=td.Box(size=[mmi_length, mmi_width, waveguide_thickness]),\n",
    "    medium=Si_medium,\n",
    ")\n",
    "\n",
    "Input_taper = td.Structure(\n",
    "    name=\"Input taper\",\n",
    "    geometry=td.PolySlab(\n",
    "        slab_bounds=[-waveguide_thickness / 2, waveguide_thickness / 2],\n",
    "        vertices=[\n",
    "            [-mmi_length / 2 - taper_length, waveguide_width / 2],\n",
    "            [-mmi_length / 2, taper_width / 2],\n",
    "            [-mmi_length / 2, -taper_width / 2],\n",
    "            [-mmi_length / 2 - taper_length, -waveguide_width / 2],\n",
    "        ],\n",
    "    ),\n",
    "    medium=Si_medium,\n",
    ")\n",
    "\n",
    "Output_taper__top_ = td.Structure(\n",
    "    name=\"Output taper (top)\",\n",
    "    geometry=td.PolySlab(\n",
    "        slab_bounds=[-0.11, 0.11],\n",
    "        vertices=[\n",
    "            [mmi_length / 2, output_separation / 2],\n",
    "            [mmi_length / 2, output_separation / 2 + taper_width],\n",
    "            [\n",
    "                mmi_length / 2 + taper_length,\n",
    "                (output_separation + waveguide_width + taper_width) / 2,\n",
    "            ],\n",
    "            [\n",
    "                mmi_length / 2 + taper_length,\n",
    "                (output_separation - waveguide_width + taper_width) / 2,\n",
    "            ],\n",
    "        ],\n",
    "    ),\n",
    "    medium=Si_medium,\n",
    ")\n",
    "\n",
    "Output_taper__bottom_ = td.Structure(\n",
    "    name=\"Output taper (bottom)\",\n",
    "    geometry=td.PolySlab(\n",
    "        slab_bounds=[-0.11, 0.11],\n",
    "        vertices=[\n",
    "            [mmi_length / 2, -output_separation / 2],\n",
    "            [mmi_length / 2, -output_separation / 2 - taper_width],\n",
    "            [\n",
    "                mmi_length / 2 + taper_length,\n",
    "                -(output_separation + waveguide_width + taper_width) / 2,\n",
    "            ],\n",
    "            [\n",
    "                mmi_length / 2 + taper_length,\n",
    "                -(output_separation - waveguide_width + taper_width) / 2,\n",
    "            ],\n",
    "        ],\n",
    "    ),\n",
    "    medium=Si_medium,\n",
    ")\n",
    "\n",
    "Input_straight_waveguide = td.Structure(\n",
    "    name=\"Input straight waveguide\",\n",
    "    geometry=td.Box(center=[-32.9, 0, 0], size=[14.2, 0.5, 0.22]),\n",
    "    medium=Si_medium,\n",
    ")\n",
    "\n",
    "Output_straight_waveguide__top_ = td.Structure(\n",
    "    name=\"Output straight waveguide (top)\",\n",
    "    geometry=td.Box(center=[32.9, 0.9, 0], size=[14.2, 0.5, 0.22]),\n",
    "    medium=Si_medium,\n",
    ")\n",
    "\n",
    "Output_straight_waveguide__bottom_ = td.Structure(\n",
    "    name=\"Output straight waveguide (bottom)\",\n",
    "    geometry=td.Box(center=[32.9, -0.9, 0], size=[14.2, 0.5, 0.22]),\n",
    "    medium=Si_medium,\n",
    ")\n",
    "\n",
    "grid_spec = td.GridSpec.auto(min_steps_per_wvl=11)\n",
    "\n",
    "sim_base = td.Simulation(\n",
    "    size=[55.6, 6.925, 2.1],\n",
    "    symmetry=[0, -1, 0],\n",
    "    run_time=3e-12,\n",
    "    grid_spec=grid_spec,\n",
    "    medium=SiO2_medium,\n",
    "    sources=[Mode_source],\n",
    "    monitors=[Flux_monitor, Mode_monitor, Field_monitor, Flux_monitorRef],\n",
    "    structures=[\n",
    "        MMI_section,\n",
    "        Input_taper,\n",
    "        Output_taper__top_,\n",
    "        Output_taper__bottom_,\n",
    "        Input_straight_waveguide,\n",
    "        Output_straight_waveguide__top_,\n",
    "        Output_straight_waveguide__bottom_,\n",
    "    ],\n",
    "    boundary_spec=td.BoundarySpec(x=td.Boundary.pml(), y=td.Boundary.pml(), z=td.Boundary.pml()),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAACuCAYAAAC7pkCSAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAARAhJREFUeJzt3Xd8FGX+B/DPM7MtdUMJCSGUSK+CdAFByQ88I8WGclRFUIoecp6IUgQLKKdyKuW8U7hDFMSj+APkREAFBBSI/mgiQugkoaSXLTPf3x+bnd3NbpLd9ITv+/WKZp955pln57uzfPPMzDOCiAiMMcYYY6zGk6q6A4wxxhhjrHxwYscYY4wxVktwYscYY4wxVktwYscYY4wxVktwYscYY4wxVktwYscYY4wxVktwYscYY4wxVktwYscYY4wxVktwYscYY4wxVktwYscYq7W+/fZbCCHw7bffVnVXWDn7/PPPUbduXWRnZ1f6tk+cOAGdTodjx45V+rYZKwkndoyxGm/ZsmVYtWpVVXejVD799FMsWbKkqrsBAFBVFW+99Rbi4uJgMpnQqVMnfPbZZ36vn56ejkmTJiEyMhIhISG4++67ceTIEZ91v/zyS9xxxx0wmUxo0qQJ5s2bB7vd7td2FEXBvHnz8MwzzyA0NNTv/pWXdu3aISEhAXPnzq30bTNWEsHPimWM1XQdOnRA/fr1vUbmVFWF1WqFwWCAJFXPv2Pvv/9+HDt2DOfOnavqrmDWrFlYtGgRJk6ciO7du2Pz5s3YunUrPvvsMzz22GPFrquqKvr164dffvkFf/nLX1C/fn0sW7YMFy9exOHDh9GyZUut7ldffYWEhAQMGDAAI0eOxNGjR7F06VJMmjQJy5cvL7GfmzZtwoMPPoiLFy+iUaNGZX7fpfHVV1/hvvvuw++//47mzZtXSR8Y84kYY7VCdnZ2VXehyrRv35769+9f1d0olYSEBGratGlVd4MuXbpEer2epk6dqpWpqkr9+vWj2NhYstvtxa6/bt06AkDr16/XylJTUykiIoJGjhzpUbddu3Z0++23k81m08pefvllEkLQyZMnS+zr0KFDqW/fvv6+tQphtVqpTp06NGfOnCrtB2OFVc8/YRm7xV2+fBkTJkxATEwMjEYj4uLiMHnyZFitVgDAqlWrIITAd999hylTpqBBgwaIjY3V1l+2bBnat28Po9GImJgYTJ06Fenp6R7bOH36NB566CFER0fDZDIhNjYWjz32GDIyMrQ6O3bsQN++fREREYHQ0FC0bt0aL730Uon992c9i8WCefPmoUWLFjAajWjcuDFeeOEFWCwWr/Y++eQT9OjRA8HBwahTpw7uuusufP311wCAZs2a4fjx4/juu+8ghIAQAgMGDABQ9DV269evR9euXREUFIT69etj9OjRuHz5sked8ePHIzQ0FJcvX8bw4cMRGhqKyMhIPP/881AUpcR9sHnzZiQkJGgxbN68OV599VWPdQcMGICtW7fi/PnzWt+bNWtWZJvjx4/X6hX+eeWVV0rsU0n9tdlsmDJlilYmhMDkyZNx6dIl7N+/v9j1v/jiC0RFReHBBx/UyiIjIzFixAhs3rxZi+uJEydw4sQJTJo0CTqdTqs7ZcoUEBG++OKLYreTn5+P7du3Iz4+3qP83LlzEEL4PCVfeP+88sorEELgt99+w+jRo2E2mxEZGYk5c+aAiHDx4kUMGzYM4eHhiI6Oxttvv+3Vpl6vx4ABA7B58+Zi+8tYZdOVXIUxVpmuXLmCHj16aNcrtWnTBpcvX8YXX3yB3NxcGAwGre6UKVMQGRmJuXPnIicnB4DjH6358+cjPj4ekydPxqlTp7B8+XL89NNP2LdvH/R6PaxWKwYPHgyLxYJnnnkG0dHRuHz5MrZs2YL09HSYzWYcP34c999/Pzp16oQFCxbAaDTi999/x759+4rtvz/rqaqKoUOHYu/evZg0aRLatm2Lo0eP4t1338Vvv/2GTZs2aXXnz5+PV155BXfeeScWLFgAg8GAgwcPYteuXRg0aBCWLFmiXWv18ssvAwCioqKK7N+qVavw+OOPo3v37li4cCFSUlLwt7/9Dfv27UNiYiIiIiK0uoqiYPDgwejZsyf++te/4ptvvsHbb7+N5s2bY/LkycXuh1WrViE0NBQzZsxAaGgodu3ahblz5yIzMxOLFy8GALz88svIyMjApUuX8O677wJAsdeMPfXUU14Jzfbt27FmzRo0aNBAK7t+/XqxfXMKCwuD0WgEACQmJiIkJARt27b1qNOjRw9ted++fYtsKzExEXfccYfXKe8ePXrgww8/xG+//YaOHTsiMTERANCtWzePejExMYiNjdWWF+Xw4cOwWq244447/HqPxXn00UfRtm1bLFq0CFu3bsVrr72GunXr4u9//zvuuecevPnmm1izZg2ef/55dO/eHXfddZfH+l27dsXmzZuRmZmJ8PDwMveHsXJR1UOGjDFPY8eOJUmS6KeffvJapqoqERGtXLmSAFDfvn09TpGlpqaSwWCgQYMGkaIoWvkHH3xAAOjjjz8mIqLExESv02aFvfvuuwSArl27FlD//Vlv9erVJEkS7dmzx6N8xYoVBID27dtHRESnT58mSZLogQce8Hg/RK59QVT0qdjdu3cTANq9ezcROU6fNWjQgDp06EB5eXlavS1bthAAmjt3rlY2btw4AkALFizwaLNLly7UtWvX4ncCEeXm5nqVPfXUUxQcHEz5+flaWVlOxZ4+fZrMZjP9z//8j8fnAIBfPytXrvTox2233ea1jZycHAJAL774YrF9CQkJoSeeeMKrfOvWrQSAtm/fTkREixcvJgB04cIFr7rdu3enXr16Fbudf/7znwSAjh496lGelJTk9Z6cANC8efO01/PmzSMANGnSJK3MbrdTbGwsCSFo0aJFWnlaWhoFBQXRuHHjvNr99NNPCQAdPHiw2D4zVpn4VCxj1Yiqqti0aROGDBniNaIBOE4puZs4cSJkWdZef/PNN7BarZg+fbrHyMnEiRMRHh6OrVu3AgDMZjMA4L///S9yc3N99sU5crV582aoqur3e/BnvfXr16Nt27Zo06YNrl+/rv3cc889AIDdu3cDcFwkr6oq5s6d6zUSVHhf+OPQoUNITU3FlClTYDKZtPKEhAS0adNG2z/unn76aY/X/fr1w9mzZ0vcVlBQkPZ7VlYWrl+/jn79+iE3Nxe//vprwH0vLCcnBw888ADq1KmDzz77zONzsGPHDr9+Bg8erK2Tl5enjd65c+6nvLy8Yvvj7/rO/xdVt6Tt3LhxAwBQp06dYuv548knn9R+l2UZ3bp1AxFhwoQJWnlERARat27tM+bOPvg7QspYZeBTsYxVI9euXUNmZiY6dOjgV/24uDiP1+fPnwcAtG7d2qPcYDDgtttu05bHxcVhxowZeOedd7BmzRr069cPQ4cO1a43Ahynqf75z3/iySefxIsvvoiBAwfiwQcfxMMPP1zsHab+rHf69GmcPHkSkZGRPttITU0FAJw5cwaSJKFdu3Z+7Y+SFLV/AKBNmzbYu3evR5nJZPLqY506dZCWllbito4fP47Zs2dj165dyMzM9Fjmfh1jaU2cOBFnzpzBDz/8gHr16nksK3y61h9BQUE+r2/Mz8/XlpfH+s7/F1W3pO04UTlM6NCkSROP12azGSaTCfXr1/cqdyaUvvpQmj8yGKsonNgxVoP5+4+gL2+//TbGjx+PzZs34+uvv8azzz6LhQsX4sCBA4iNjUVQUBC+//577N69G1u3bsX27duxbt063HPPPfj66689RogK96mk9VRVRceOHfHOO+/4bKNx48alfl/lqaj3WJL09HT0798f4eHhWLBgAZo3bw6TyYQjR45g5syZAY2A+vK3v/0Nn332GT755BN07tzZa3lycrJf7ZjNZu0z1LBhQ+zevRtE5JGoXL16FYDjGrjiNGzYUKvrrvD6DRs21MoLx/nq1avaNX1FcSaxaWlpHjcMFaW4BNBXfIuKua92nAl+4USQsarEp2IZq0YiIyMRHh5e6hntmzZtCgA4deqUR7nVakVSUpK23Kljx46YPXs2vv/+e+zZsweXL1/GihUrtOWSJGHgwIF45513cOLECbz++uvYtWuXdqq0KCWt17x5c9y8eRMDBw5EfHy8149zRK158+ZQVRUnTpwodnv+jpgUtX+cZYX3T2l9++23uHHjBlatWoU//elPuP/++xEfH+/z9GGgoz179uzB888/j+nTp2PUqFE+6zRs2NCvn3Xr1mnrdO7cGbm5uTh58qRHWwcPHtSWF6dz5844cuSIV9J68OBBBAcHo1WrVh7tHDp0yKPelStXcOnSpRK306ZNGwBAUlKSz+VZWVker1NSUoptryySkpIgSZL23hirDjixY6wakSQJw4cPx//+7/96/cMHlHz6KT4+HgaDAe+9955H3Y8++ggZGRlISEgAAGRmZnrN8t+xY0dIkqSdIrt586ZX+85/dH2dRnPyZ70RI0bg8uXL+Mc//uFVNy8vT7vDd/jw4ZAkCQsWLPBKGNzfX0hIiNd0Lr5069YNDRo0wIoVKzzew1dffYWTJ09q+6esnKM+7n20Wq1YtmyZV92QkBC/T81evXoVI0aMQN++fbU7a30pzTV2w4YNg16v9+gjEWHFihVo1KgR7rzzTo9+/Prrr7DZbFrZww8/jJSUFGzYsEEru379OtavX48hQ4Zo19S1b98ebdq0wYcffugx9cvy5cshhMDDDz9c7D7o2rUrDAaDz+MDgNcfHRs3btTeS3k7fPgw2rdvr12+wFh1wKdiGatm3njjDXz99dfo37+/NhXI1atXsX79euzdu9djOo7CIiMjMWvWLMyfPx/33nsvhg4dilOnTmHZsmXo3r07Ro8eDQDYtWsXpk2bhkceeQStWrWC3W7H6tWrIcsyHnroIQDAggUL8P333yMhIQFNmzZFamoqli1bhtjY2GKnvfBnvTFjxuDzzz/H008/jd27d6NPnz5QFAW//vorPv/8c/z3v/9Ft27d0KJFC7z88st49dVX0a9fPzz44IMwGo346aefEBMTg4ULFwJw/GO/fPlyvPbaa2jRogUaNGig3YjhTq/X480338Tjjz+O/v37Y+TIkdp0J82aNcNzzz1X2rB5uPPOO1GnTh2MGzcOzz77LIQQWL16tc/komvXrli3bh1mzJiB7t27IzQ0FEOGDPHZ7rPPPotr167hhRdewNq1az2WderUCZ06dQJQumvsYmNjMX36dCxevBg2mw3du3fHpk2bsGfPHqxZs8bjFOWsWbPwr3/9C0lJSdq8ew8//DB69eqFxx9/HCdOnNCePKEoCubPn++xrcWLF2Po0KEYNGgQHnvsMRw7dgwffPABnnzySa/pVgozmUwYNGgQvvnmGyxYsMBr+fbt2zFq1Cjcdddd+O233/Dhhx8iODgYX3/9Nbp37477778/4H3ji81m0+aRZKxaqZJ7cRljxTp//jyNHTuWIiMjyWg00m233UZTp04li8VCRK7pTnxNiULkmN6kTZs2pNfrKSoqiiZPnkxpaWna8rNnz9ITTzxBzZs3J5PJRHXr1qW7776bvvnmG63Ozp07adiwYRQTE0MGg4FiYmJo5MiR9NtvvxXbd3/Xs1qt9Oabb1L79u3JaDRSnTp1qGvXrjR//nzKyMjwqPvxxx9Tly5dtHr9+/enHTt2aMuTk5MpISGBwsLCCIA29Unh6U6c1q1bp7VXt25dGjVqFF26dMmjzrhx4ygkJMTr/TmnyijJvn37qFevXhQUFEQxMTH0wgsv0H//+1+v/mRnZ9Mf//hHioiIIADFTn3Sv3//IqcucZ/Oo7QURaE33niDmjZtSgaDgdq3b0+ffPKJVz3nVDBJSUke5Tdv3qQJEyZQvXr1KDg4mPr371/kZ3Tjxo3UuXNnMhqNFBsbS7Nnzyar1epXPzds2EBCCI8pU5zTnbzxxhsUHx9PRqOR4uLi6IsvvqCXXnqJgoODaf78+UTkimHhKXmKinn//v2pffv2HmVfffUVAaDTp0/71WfGKgs/K5YxxliNoigK2rVrhxEjRuDVV18F4HjyRFxcHFauXInx48dXeB+GDx8OIYR2qpex6oKvsWOMMVajyLKMBQsWYOnSpcjOzq707Z88eRJbtmzRkkrGqhNO7BhjjNU4jz76KG7evFnsI9gqStu2bWG32/2eb5KxysSJHWOMMcZYLcHX2DHGGGOM1RI8YscYY4wxVktwYscYY4wxVkvwBMV+UlUVV65cQVhYGD/wmTHGGGOVhoiQlZWFmJgYSFLxY3Kc2PnpypUreP/996HTuXYZESEtLa1cHlVjMBi0YKmqCqvVWuY2A6XT6Tzen8ViqZDH8BRHCKE9eggA7Ha716OvKgPHw4Hj4cLxcOF4OHA8XDgeLmWNhxACderU8RpEeuONN3Dx4kXExsYWuz4ndn4KCwuDTqdDsMkAZ7JMBDSLjoEslX0Ez2onZOQ6noVpDpZg0FX+qKBKhBtZjj4YdALm4Ko5U5+Rq8Jqd3wh1AuTIFXBCCnHw4Xj4cDxcOF4uHA8HDgeLmWNh6ISrmeqcO+681HZYWFhJa7PiZ2fnJmzJAF6nQRVJagEGA26Mh9EVpuK9FwFQQbHgZBrBYKNEgz6yjswVCLcyFSg10kw6gTybQRFFQgLlkteuRxl5SpQSSDUJMFiJ+RZBeqFy5V6cHI8XDgeDhwPF46HC8fDgePhUh7xsNoJkmSDJACpYODIZndkdv5cCsY3T1Qxq03F9SwFelmgvlmH+mYd9LLA9SwFVptaKX1wHpQ2hVA/TEa9cB3CgyRk5qnIylUqpQ+A46DMzFMRHiShXrgO9cNk2BRH39RKGtLneLhwPBw4Hi4cDxeOhwPHw6U6xAPgxK5KuX8InH9VSMLxe2V9GAoflM6/KsKC5Uo9ON0PSudfeQa9VKkHJ8fDhePhwPFw4Xi4cDwcOB4u1SEeTpzYVRFfHwKnyvowFHVQOlXWwenroHSqrIOT4+HC8XDgeLhwPFw4Hg4cD5fqEA93nNhVgeI+BE4V/WEo6aB0quiDs7iD0qmiD06OhwvHw4Hj4cLxcOF4OHA8XKpDPLy2V6GtMy/+fAicKurD4O9B6VRRB6c/B6VTRR2cHA8XjocDx8OF4+HC8XDgeLhUh3j43FaFtcy8BPIhcCrvD0OgB6VTeR+cgRyUTuV9cHI8XDgeDhwPF46HC8fDgePhUh3iUeR2KqTVWkqWZZT2s1CaD4FTeX0YSntQOpXXwVmag9KpvA5OjocLx8OB4+HC8XDheDhwPFyqQzyKw/PYBcimAJLkCAQRYLOX/MGw2VSk5SjQyQLhwRLsCgAE/oEKD5aQlq0gNcOOOiEy9AEcWESEtGwFdoVQJ0QGhNAmcQyE0SAhWCGk5zjaCgkK7MDKyVOQna8i1CTBaJBK1QcIAXOQhLQcBanphDqhckCPeeN4uHA8HDgeLhwPF46HA8fDpTLiYbMTiAgqABTkfYHkoJzYBcBut0OWZdjsgCw5QmmxqVDUoj8UqgpY7SqMesfs0xYboTQfAqeggg9zjkWFQQFKeGQcAMcHwmonSAIINclQCFCspf8rQZYFQkwSrAqB8hwfcH/YFYJNIYSYJMiyQH4Z+gA43ovVriIzV4VBJ+DPscnxcOF4OHA8XDgeLhwPB46HS2XFw644HoBQ2mmVObELQHp6Otq3iEauFZAEYNRLuK+7GeaQyp1dmzHGGGO1U0aOgm0/ZcCol6AveLJVvsX/Z95yYhcAVVVhMkgIDZKQmmGHSoA5REbdMN6NjDHGGCsfOllArxPaI0utNv/X5ZsnSsGgl1AnRAZV0mNKGGOMMXZrUomQmev/qWceaiolvV6CofIeg8cYY4yxW4RakMc570ZWAriksFqN2C1cuBDdu3dHWFgYGjRogOHDh+PUqVPFrrNq1SoIITx+TCaTRx0iwty5c9GwYUMEBQUhPj4ep0+fLnN//bnwlDHGGGMsEFa7CqtN1aaYMQf7n3BUq9Tku+++w9SpU3HgwAHs2LEDNpsNgwYNQk5OTrHrhYeH4+rVq9rP+fPnPZa/9dZbeO+997BixQocPHgQISEhGDx4MPLz8yvy7TDGGGOMBUwIgbRsBVa7Y95Af+8eBqrZqdjt27d7vF61ahUaNGiAw4cP46677ipyPSEEoqOjfS4jIixZsgSzZ8/GsGHDAAD//ve/ERUVhU2bNuGxxx4rvzfAGGOMMVaFqtWIXWEZGRkAgLp16xZbLzs7G02bNkXjxo0xbNgwHD9+XFuWlJSE5ORkxMfHa2Vmsxk9e/bE/v37i2zTYrEgMzPT44cxxhhjrKIROSZPNugcT6iwK/7frFltEztVVTF9+nT06dMHHTp0KLJe69at8fHHH2Pz5s345JNPoKoq7rzzTly6dAkAkJycDACIioryWC8qKkpb5svChQthNpu1n8aNG/voY2neGWOMMcZY0Qw6CQa9pD1+LCOAu2KrbWI3depUHDt2DGvXri22Xu/evTF27Fh07twZ/fv3x4YNGxAZGYm///3vZdr+rFmzkJGRof1cvHjRY7nNpsJq58yOMcYYY+XLeXOm89mycgDZWrVM7KZNm4YtW7Zg9+7diI2NDWhdvV6PLl264PfffwcA7dq7lJQUj3opKSlFXpcHAEajEeHh4R4/TtaCZ8UF8nw5xhhjjLFAScLxXFq/61dgXwJGRJg2bRo2btyIXbt2IS4uLuA2FEXB0aNH0bBhQwBAXFwcoqOjsXPnTq1OZmYmDh48iN69ewfcvl0hXM9yPN/OOSM0Y4wxxlhFkQIYSAr4rtikpCTs2bMH58+fR25uLiIjI9GlSxf07t3ba/64QE2dOhWffvopNm/ejLCwMO0aOLPZjKCgIADA2LFj0ahRIyxcuBAAsGDBAvTq1QstWrRAeno6Fi9ejPPnz+PJJ58E4Lhjdvr06XjttdfQsmVLxMXFYc6cOYiJicHw4cMD6p8QjvPcJoOM8GCp4AHAjDHGGGPVg9+J3Zo1a/C3v/0Nhw4dQlRUFGJiYhAUFISbN2/izJkzMJlMGDVqFGbOnImmTZuWqjPLly8HAAwYMMCjfOXKlRg/fjwA4MKFC5DcZgZOS0vDxIkTkZycjDp16qBr16744Ycf0K5dO63OCy+8gJycHEyaNAnp6eno27cvtm/fHnAiajQaIUtAvXAZdgUAOLFjjDHGWPUhyI8Hnnbp0gUGgwHjxo3DkCFDvO4QtVgs2L9/P9auXYv//Oc/WLZsGR555JEK63RVyMzMxOzZs9GxZSRCg/Sw2gn5VhWP9KuDumHVajpAxhhjjNVQN7PsWL8nDSaDpF3ylZtvx5PPzEFGRobHNf+++JWRLFq0CIMHDy5yudFoxIABAzBgwAC8/vrrOHfunP/voAaxWq0BnedmjAVmzb6TVd0FVkqj+rSt6i4wxuBnYldcUldYvXr1UK9evVJ3iDHGGGOsuiNFAVLTgLrhEEZD8ZXzrwOXN3uWNRoGmOoXvw2LDQHNdYIyPFIsNTUVqampUAvN0tupU6fSNskYY6yGupqeA0kAkpAgSYCAgCQAIQlIEJCEgJAACQJCCLc6AkKAp49iNQZl54LOXQUlXQHyLRBdWkO08H6IQZm3QwBuZgIgSLoSEkc3ASd2hw8fxrhx43Dy5Ek4L88TQoCIIISAoiiBNskYY6yG23X8QpnWl4RwJHqS498UR8JXkBwW/F8SEoSAlgy66jiXF1pHEp4JpihIJqWC9gFHHSEgFfzuTz9c2yzcj0J9ck9wC8pZzUSqCiTfhJp0Gbh6HbDZHbMIqyqgVuSNlASoBGG1+r1GwIndE088gVatWuGjjz5CVFQUf1AZY4yVmUoEgFCbxwacI5RawicBQkiQ4JZ4FkomvRNLz4TSPWmUROE2CtbXRk2d23Elr+7lHv0TjnWE5J1Au/ojaX2srckr5VlA56+Czl4GsnIdw2g6GQg2AUIAufkV3wlJIJBphwNO7M6ePYv//Oc/aNGiRaCrMsYYY7csAkEht+RV0f5Tq0gFCZ/zlLvr9HvBKKnzd1GwTCqUiPozGlvEKGlRbQhRMCJb6BIBzz4CkiRBgCDfyIB8IRnylWuATXHkVQY9hCx7vlki0PGzUE4mFbtPRItISG4zrJFNQN15HLB711VUPSi/IQDFMatagLlywIndwIED8csvv3BixxhjjDEvKqmus5M1KG/VKwoa5eSgaWYWwq02qETIEwI2SQJIAHk2ADbXCgLQqQQpJ8/xoojZ4/Sqigu/22FokKmV6ax61LmZDZCAWmiUMxcGqEKFAhV2qFCFgI38fzZ9wIndP//5T4wbNw7Hjh1Dhw4doNfrPZYPHTo00CYZY4wxxqqE2WJBk6xsxGZnw6CoIAA2SYLiOCdd9IoE2IVU4oiaXMx0wVZJgip5nma1kQyCgApADXS4DqVI7Pbv3499+/bhq6++8lrGN08wxhhjrKYwKAp6Jacg2OY4J5qr04Fq+HWCgU2OAuCZZ57B6NGjcfXqVaiq6vHDSR1jjDHGagqrLOPHqAZICg+DRZZhUhSY7HZIqv+nPqubgEfsbty4geeeew5RUVEV0Z9aI0/Jw8k8z1n02+qbIijjuGfFul0BXYj2cuvvl/B/19LwQKsmaFPPXBldZYwxxm5ZaSYT0kwmnLTbEZudgyZZWQiz2SDsKuzOa+yKGMWTVRUSCAJFXmIHqZhTsXoiqIUGxXRQIAQ5puABeV2DV5KAE7sHH3wQu3fvRvPmzQNdlflhx7mr+NexM/jHL7+jV0x9jGzXDIOaxcCok0tembFSUomgEoEIUFUCoeC1CqhwlquO/4OgqgAROZapBFUrL6hL5FhO7q89y1WCW51S3PrFGKvWSn13LOAqr4q7Y4kgbmRAd+EqjCk3Iex2kACEQe+Y6gRu31Y5eYDBUMLXl0Cr5g0gBZ3TSsgmoGYEA3bvpM+m6iFZBWQI6OC8hq8Cpztp1aoVZs2ahb1796Jjx45eN088++yzgTbJCjHIEvSSwK7zyfj2Qgpiw4LxaNumeLh1UzQ1h1Z1924p5J6QeCU0BQmRV0LjuiuMfCQwnq/dkh2VHBfLFiRUzuTKIwlSCaTVKWgDBcmYW11n264y17Z8JVnVgSwCvjKEsRqp8Hx2jrniyjafna8kKND57JxJT1nmsxNu263RwoKAZtGg7DzHPHZJl4HsPMBqBfQFCZ7jzUJ0bAGpVZPi28u/Dlx2vRR6gjyoo89HislZdog9aYBeAq7fCLjrpborNjQ0FN999x2+++47j2VCiHJJ7JYuXYrFixcjOTkZt99+O95//3306NGjyPrr16/HnDlzcO7cObRs2RJvvvkm7rvvPm05EWHevHn4xz/+gfT0dPTp0wfLly9Hy5Yty9zXQFzKysP1VItnYf4NQM7VXqblWwEIhBn1CDPqYVVUJOfkYfHBE1ie+BvubhKNx9o2Q/8mUdBJVfsPIbknKgVJT9EjMu5JBXkmPaoKFc7ExS2hKUh0Cic0WhIDaKNFVEQipRY7elR4HbfXbv1gjFUefgIFq05EaBBE+9tAbZoBV645njyRkuaYmFiWij7/Wp5UAlCB050kJRU/CV9ZrVu3DjNmzMCKFSvQs2dPLFmyBIMHD8apU6fQoEEDr/o//PADRo4ciYULF+L+++/Hp59+iuHDh+PIkSPo0KEDAOCtt97Ce++9h3/961+Ii4vDnDlzMHjwYJw4cQImk8mrzYoyfddx/JJSKPuW0gEIEAoSJThGXW12xy3XACHCoIdChBybHZt+u4D//f0SYkODkNCsEe6JjUI9k8HHKEzJCY2vdbxGlrRRJFUboXIfwWKMMQDo37axV/LFz4xltYWQJaBxFOTGUaCMbMezYs9dASy24qdEKfuWHVPk6fQlV3WuQVQ+6ebVq1exevVqvPDCC2Vqp2fPnujevTs++OADAI7reho3boxnnnkGL774olf9Rx99FDk5OdiyZYtW1qtXL3Tu3BkrVqwAESEmJgZ//vOf8fzzzwMAMjIyEBUVhVWrVuGxxx7zq1+ZmZl46aWX0CgmHDqd5BidUgUaNsiFXu+dSSvCitzgyx5lHxzIRWpeHsxu6bQVRhQ+dy6V8OVmU1Xk2BWoRAjSybg9woyekXXRIiyk5g9/s1uWLCQoAUzCyaqXUX3aVnUXGKtUZLcDKTeBehEQJkPxlfOvA5c3e5Y1GubzVOzNLDvW70mDySBBr9gASUKeAjz5zBxkZGQgPDy82E2V6lmxvpw/fx4//vhjmRI7q9WKw4cPY9asWVqZJEmIj4/H/v37fa6zf/9+zJgxw6Ns8ODB2LRpEwDHCGNycjLi4+O15WazGT179sT+/fuLTOwsFgssFtdp08zMTBgMOqgqIEGBY1xNQI88GHxMra3ABoU8H9pbX2+DpNph1ruSLxvsoMBnnQHgGJ21KDYkZaXiYtY1NAgyok9kPdxel++mZTWPBAlyAKcbWPVis2RUdRcYq3z1jQDyAEte8fWEHoh92Lvcx3Fjt6ogUkCqAtILAASyV+Cp2LS0NI/XiqLg7NmzOHnyJJYtWxZocx6uX78ORVG8plKJiorCr7/+6nOd5ORkn/WTk5O15c6your4snDhQsyfP9+jbP7sZyFUQFLyYZDssAs9GsknEeJ2jZwTQYVq93w48OQmKmyqAp3bqBpBKvMJTYUczx+UhUCoLgWxclAZW2Ss8jm+vlhNlZp0uqq7wFitkGkxwW5pB6vdDlV2JHSK4v/ZuIATu40bN/osf/3117Fp0yY89dRTgTZZLc2aNctjJDAzMxP//nARDLIFefYgGGQrJEFQSA87eQ/BklCguD9TDkCWHbAoBIPbAJ1rBpzAEBFsjtseoZNk1DHoUd9oQIheBzvx6VhW84iiH7XIagBJx39QMlYeJLsREBJQcJMQENhkUAEndkUZOXIkXnvttTK1Ub9+fciyjJSUFI/ylJQUREdH+1wnOjq62PrO/6ekpKBhw4YedTp37lxkX4xGI4xGo0eZJAFG2QJZUpFnC4JBZ4UMK3Sw+GhBLThl6xKuU5AvCj/5ze1fsoIFshCQfVwrR3CMzqnkSAWNsoQ6RiPMeh1kSSqoYdOmuQD4H0pWc/CIXc2m2L3PXDDGAqcoKhxzXimgguuOKYABm3JL7H755Rd06dKlTG0YDAZ07doVO3fuxPDhwwE4bp7YuXMnpk2b5nOd3r17Y+fOnZg+fbpWtmPHDvTu3RsAEBcXh+joaOzcuVNL5DIzM3Hw4EFMnjw5oP4JAehlIFhnhRAy7KqMeiESwo3e18gRAfmqZyAMZj0stkJJoKSDey6ea1eQZ7PDILvaJALsBcGVhEC4QY8IkwGhel2hu3GK+h1ahkfeRQXl5PzFRx1yq+f5i/udsZ7t+d6mq323/xIK1SPvtjx+J891OBuoHXiguUa7kWmFKBREj28k4SoRbhWEW233rzP3cvd6bs14tuZrZEMI729FrZ5rBZ/fnHwjGqsiil2ASAWBtMQu1+b/DB4BJ3aFb1QAHKNfmzdvRkJCgsfyd955J9DmMWPGDIwbNw7dunVDjx49sGTJEuTk5ODxxx8HAIwdOxaNGjXCwoULAQB/+tOf0L9/f7z99ttISEjA2rVrcejQIXz44YcAHHeXTp8+Ha+99hpatmypTXcSExOjJY/+UhRAkvSQZYEQI5BrFTAGR8JksnrXJTty7OkeZdEhoZCtntcowhQJCFcY9l++ht9zshBi0MGmqLAqKoQAwgx6NI8Iw20RoQjWl1s+XrsUTBfjTPa0CVkK8kD3UUxnQkrkSkJJa6bQsoJE06ONguWe23Rvo2A9n/UK+uaWnBa1Tfckt7aPvhb3SB5WA9TSR4UXzEELV+LpSAm1U2TClSRqvwvn7XVwJZcFy1y/u5Z5/K616/7acz3PdgutB9esCsItsS2qb6z6scgGCEkPSZIgyypyrHrkK/7/ux9whpCYmOizvHv37khNTUVqaiqA0s9F9Oijj+LatWuYO3cukpOT0blzZ2zfvl27+eHChQuQ3CbmvfPOO/Hpp59i9uzZeOmll9CyZUts2rRJm8MOAF544QXk5ORg0qRJSE9PR9++fbF9+/Yyz2GnkyvmXyGFCDlWO/SShMbhwWheJwyNQoN5KpOSeHxJav+pXbwS1kLJa8EvXoklCielRSev3kmoexm5JcrOLnkmob6SV+f2C6+n9YEIQgiU0+xLjJUbcjueCkqqsDcVwyOZdBvR9J0sFiz1WOaeNLonuSiU0BZOQgvW0UZNhdvvhRPTQsmtW/Lq7Ie2LR/9rmlsigSDrCDHqkeOVY9gnfcAUlHKbR672i4zMxMf/PUl1DfrYdAL2BQJVkXG/7Q+B3MRI3ZphUbs6kglj9gdTr6B85k5uC0iFLeZQxFuLGFuHMZqkd+T06u6C4yxWqio5NW5zCNZLPiPxwinW0JbVPKanuPrevvA5dpMOHL1dtgVGUIAJl0+dMKGP816v2LmsWP+kYUO9fXeEw8iqKF3mZs7ouqiS1RdHp1jjDHGyknhMxcFr6qkL/4Qbn0zyFaoqv85gV8z49577704cOBAifWysrLw5ptvYunSpX53gHmqFQ9PZowxxlipEQR0sh1CEHJtwQFdf+zXiN0jjzyChx56CGazGUOGDEG3bt0QExMDk8mEtLQ0nDhxAnv37sW2bduQkJCAxYsXl/a91Ch8Epsxxhhj5U0SKoyy49Ruji0EeYr/80T6ldhNmDABo0ePxvr167Fu3Tp8+OGHyMjIAOAYYWrXrh0GDx6Mn376CW3b1uLnBQrnBeEEVQVsih6k2KAq/l/UyBgrmgR7VXeBMcaqlAQFslAAEGRJRYg+B9nWEL/X9/saO6PRiNGjR2P06NEAgIyMDOTl5aFevXrQ6/UBd7wmkiWAVBsUu4rM/DAIocCSdx0WlSfmZKw8hIjyufiYMVYLuU2AWPiKJeHrN1FoucdrUeQy9+llCi8qVM2jLZ0sCtcolUyLgF4iGGQVOkkBZMXjmruSlPrmCbPZDLPZXNrVayRFBSD0yLQFQYWEYJ0dxqD6MJr4HyPGykNOTmZVd4GxGqm4uz4dv3vOvedrypLCd3265ttzm7Kk8Dx+PqYs0dYDfE9fUuguVLc1vacz8ah3azBKRghJByEBUsGzYg3lfY0dK0BAtjUIBAnhJitUkiBkPaQKms+OsVuNyl9JrAJU+0mGnQmVj/nXeJLhW4+Q9QUJuvuzYithxO5WpAo97CRQL9gxQmetpTOtM8ZYoMKC9PAeBXIb4/EY1fGd9LgnRM42fE1qWzhxciU6nPQwxoldQATCDBboZYJN8WumGMYYuyVEmf2/uJsxVnE4OwmAIBv0klrV3WCMMcYY8yngxG7cuHH4/vvvK6Iv1V4g57gZY4wxxipbwKdiMzIyEB8fj6ZNm+Lxxx/HuHHj0KhRo4roG2PsFtMiOqKqu8AYYzVawIndpk2bcO3aNaxevRr/+te/MG/ePMTHx2PChAkYNmxYrZ7TTpYBVbVBUQiqKkNVdLDkXkO+wvPYMcYYY6zsLJZgqEojqMIOBY67NAN5Vmypbp6IjIzEjBkzMGPGDBw5cgQrV67EmDFjEBoaitGjR2PKlClo2bJlQG3abDbMnj0b27Ztw9mzZ2E2mxEfH49FixYhJiamyPVeeeUVzJ8/36OsdevW+PXXX7XX+fn5+POf/4y1a9fCYrFg8ODBWLZsGaKiogLqIxEKbj+mgjuzJMg6I2Qdn6JljDHGWNnJitGRaxTkGUBgc/mV6a7Yq1evYseOHdixYwdkWcZ9992Ho0ePol27dnjrrbfw3HPP+d1Wbm4ujhw5gjlz5uD2229HWloa/vSnP2Ho0KE4dOhQseu2b98e33zzjfZap/N8W8899xy2bt2K9evXw2w2Y9q0aXjwwQexb9++gN6vqgJCyJAkAUESBMnQGcKhN5oCaocxxhhjzBcdGSAk2TFBccENm0JU4Dx2NpsNX375JVauXImvv/4anTp1wvTp0/HHP/4R4eHhAICNGzfiiSeeCCixM5vN2LFjh0fZBx98gB49euDChQto0qRJ0W9Cp0N0dLTPZRkZGfjoo4/w6aef4p577gEArFy5Em3btsWBAwfQq1cvv/vIGGOMMVadBZzYNWzYEKqqYuTIkfjxxx/RuXNnrzp33303IiIiyty5jIwMCCFKbOv06dOIiYmByWRC7969sXDhQi0RPHz4MGw2G+Lj47X6bdq0QZMmTbB///4iEzuLxQKLxfWosMxMftQRY4wxxqq3gKc7effdd3HlyhUsXbrUZ1IHABEREUhKSipTx/Lz8zFz5kyMHDlSGwn0pWfPnli1ahW2b9+O5cuXIykpCf369UNWVhYAIDk5GQaDwSs5jIqKQnJycpHtLly4UHsertlsRuPGjcv0fhhjjDHGKlrAid2YMWNgMpX9mrI1a9YgNDRU+9mzZ4+2zGazYcSIESAiLF++vNh2/vCHP+CRRx5Bp06dMHjwYGzbtg3p6en4/PPPy9S/WbNmISMjQ/u5ePFimdpjjDHGGKtoVfZIsaFDh6Jnz57aa+dceM6k7vz589i1a1exo3W+REREoFWrVvj9998BANHR0bBarUhPT/cYtUtJSSnyujwAMBqNMBqNAW2bMcYYY6wqVVliFxYWhrCwMI8yZ1J3+vRp7N69G/Xq1Qu43ezsbJw5cwZjxowBAHTt2hV6vR47d+7EQw89BAA4deoULly4gN69ewfUtiQBRApUlUCqBFIBuzUTNpEfcD8ZY4wxxgqzW00gVQEJG1ThuCuWqJKmOylPNpsNDz/8MI4cOYItW7ZAURTtGri6devCYDAAAAYOHIgHHngA06ZNAwA8//zzGDJkCJo2bYorV65g3rx5kGUZI0eOBOC423bChAmYMWMG6tati/DwcDzzzDPo3bt3wHfECgEQqSAiEASIVCh2CxQ5rxz3BGOMMcZuVYrdkV8QCEQ1OLG7fPkyvvzySwDwuilj9+7dGDBgAADgzJkzuH79urbs0qVLGDlyJG7cuIHIyEj07dsXBw4cQGRkpFbn3XffhSRJeOihhzwmKA6UogCSpIcsC6iQIMkyjMGRMJmsgb9hxhhjjLFCLLIBkqyHJEmQZUdip6gVOI9dRWnWrBmISu74uXPnPF6vXbu2xHVMJhOWLl2KpUuXlrZ7jDHGGGPVXsB3xd7KCP4PhTLGGGOMVTZO7AJAQg+byruMMcYYY9VTtTkVWxMQgEyLEZJkAwAoqkCWRV+1nWKMMcZYrZFl0UNRBexCwDn+Zq+J19jVBKoIgk3R4Vq2AZJQAQHsOxsLWSp6h9sVAbsqQSep0Mn+B6Y4gbapkoBNkSBA0MsqRDmcUQ60TSLApkggCOhlFVIADzQuzzY5Hg4cDxeOhwvHw4Hj4cLxcKmseCiqQKbFAAFo+8afexCcOLELiIAkFCgkQyEJMqnQy3bIBfPMFJZn0yHXpkOw3oYgvb3ceqGXnG3rEYzi27arEnKtRsgSIcxgKZeD0kkWErIsRigkEGbIL7JtIiDLZoKiCoQZLdBJvvdXaeglIMtqQq5VLrFtjocDx8OF4+HC8XDgeLhwPFwqMx6SkCBIX5BsOhO6GjjdSU0gBCBBhU5SYFENUCGgKjYYdN7TneTaTMi3yQjV5yJYX/4TGIfqAQkm5NpMkGDzuQ2bKiPbEgqdsMJsyC7XgxIAZAmQjPnItIQi2yLDbPTeBhGQYQmFSirMxmzoJaV8OyGACEMeMiyhyLboEF7ENjgeDhwPF46HC8fDgePhwvFwqex4GOR8CGGCECpEQWIXyPX9nNgFgAgQkoCQBAxQHB98Wyj0ujzoZVemn2PVI9duQIjRihADAaiYR5OFygQh2ZFjDYaQdAgx2LRlNkVCljUIOlmF2WSDJCqmD0YZiJCsyMgPQqYtAmZTHqSCg1MlIDM/CCokRATlQS/rUFEfuQjZhox8HbKsZphNHA+OB8fDHcfDgePhwvFwqW7xIGGHY4ROAgRBUQUUyH63xYldwASEcCR4MhFkCcjID0JEkAV6WXV8CGx6hBptCDEoqOgbj0ONCoSwI8dqhBASQgw22BQJGflG6CSCOcgCSVRsHww6ICLIgvQ8IzLzg2EOcvxFk5nvGD537BugIveFLBx9yMgzcTw4HhqOhwPHw4Xj4cLxcKhu8ciyGKCQBElSoZCATZXhOiVbMk7sykBIQLjRghyrAel5RuhlFVZFRojB5vHXT0VzbivHqodNkWBTpIKDMl/766ei6WVVOzjT80wA4HZQlt81EcWRBGAOykdGnonjwfHQcDwcOB4uHA8XjodDdYqHogqk2U2w2h03dwgQZOH/aWielK2MRMGHgSBgVWQYZKVSPwROIQYbDLICqyKDICr1oHTSyyrMJgvsqgS7KsFsqryD0knieGg4Hi4cDweOhwvHw4Xj4VBd4uG4gYK0hyIYZCWgeHBiVw7ybK657Jx/8VS2wtt171NlUQnIsRq01zlWAwKYeqfccDwcOB4uHA8XjocDx8OF4+FSHeJhVwTcZzexq/5fXwdwYldmeTYdcqx6hBhsqB+SC51ESM8zVuqHwaZISM9zXBNRPyQXIQYbcqx65Fgr7+BUCcjIM8GuCtQJykedoHzYVYGMPFOlHpzO983x4Hg4cTxcOB4OHA8XjodLdYlHpsVxg4ZBUiALKrjOzv8+cGJXBooqkGvTa+fgncO4lflhcD8oncPnzv5U1sHpflA6r4lwXjNRmQen+0HJ8eB4ABwPdxwPB46HC8fDpTrFQ5YIsiAICTDoHMkdBZCu8c0TfiIiCAEoBJBdhY0cu04v5UGmfORbXHVNkhXZSjBu5OgRos+BvhwnVHRnUyXk2EIgCxtMUi6sbtPpybBCL4zIzDPCarMjSGcpuqEyIAKybcFQiBCiz4ZiV6G4zd1okq3IsYXgRrYOofrccp8LySnPbkS+XYZJlweZLBwPjgfHww3Hw4Hj4cLxcKlu8TDK+bAgHIpCUEGQoEIpSNf8eQKFoECeU3ELu3TpEho3blzV3WCMMcbYLerixYuIjY0ttg4ndn5SVRVXrlxBWFgYREX92VCLZWZmonHjxrh48SLCw8Oruju3HN7/VYv3f9XjGFQt3v9lQ0TIyspCTEwMJKn407J8KtZPkiSVmCWzkoWHh/NBXYV4/1ct3v9Vj2NQtXj/l57ZbParHt88wRhjjDFWS3BixxhjjDFWS3BixyqF0WjEvHnzYDRWzAOUWfF4/1ct3v9Vj2NQtXj/Vx6+eYIxxhhjrJbgETvGGGOMsVqCEzvGGGOMsVqCEzvGGGOMsVqCEzvGGGOMsVqCEztWIpvNhpkzZ6Jjx44ICQlBTEwMxo4diytXrnjUu3nzJkaNGoXw8HBERERgwoQJyM7OLrbt/Px8TJ06FfXq1UNoaCgeeughpKSkeNS5cOECEhISEBwcjAYNGuAvf/kL7HZ7ES3WThs2bMCgQYNQr149CCHw888/e9XxZ18WRkSYO3cuGjZsiKCgIMTHx+P06dMedUoT11vB0qVL0axZM5hMJvTs2RM//vhjsfXXr1+PNm3awGQyoWPHjti2bZvHcn9icav6/vvvMWTIEMTExEAIgU2bNnksL+2+KymGpTmmapuFCxeie/fuCAsLQ4MGDTB8+HCcOnXKow5/91QzxFgJ0tPTKT4+ntatW0e//vor7d+/n3r06EFdu3b1qHfvvffS7bffTgcOHKA9e/ZQixYtaOTIkcW2/fTTT1Pjxo1p586ddOjQIerVqxfdeeed2nK73U4dOnSg+Ph4SkxMpG3btlH9+vVp1qxZFfJeq6t///vfNH/+fPrHP/5BACgxMdGrTkn70pdFixaR2WymTZs20S+//EJDhw6luLg4ysvL0+qUJq613dq1a8lgMNDHH39Mx48fp4kTJ1JERASlpKT4rL9v3z6SZZneeustOnHiBM2ePZv0ej0dPXpUq+NPLG5V27Zto5dffpk2bNhAAGjjxo0ey0uz7/yJYWmOqdpm8ODBtHLlSjp27Bj9/PPPdN9991GTJk0oOztbq8PfPdULJ3asVH788UcCQOfPnyciohMnThAA+umnn7Q6X331FQkh6PLlyz7bSE9PJ71eT+vXr9fKTp48SQBo//79ROT4QpckiZKTk7U6y5cvp/DwcLJYLBXx1qq1pKQkn4mdP/uyMFVVKTo6mhYvXuzRjtFopM8++4yIShfXW0GPHj1o6tSp2mtFUSgmJoYWLlzos/6IESMoISHBo6xnz5701FNPEZF/sWAOhRO70u67kmJYmmPqVpCamkoA6LvvviMi/u6pjvhULCuVjIwMCCEQEREBANi/fz8iIiLQrVs3rU58fDwkScLBgwd9tnH48GHYbDbEx8drZW3atEGTJk2wf/9+rd2OHTsiKipKqzN48GBkZmbi+PHjFfDOaiZ/9mVhSUlJSE5O9ljHbDajZ8+eHvs/0LjWdlarFYcPH/bYb5IkIT4+vsh9vX//fo/6gONz7KzvTyyYb6XZd/7EsDTH1K0gIyMDAFC3bl0A/N1THXFixwKWn5+PmTNnYuTIkdrDnJOTk9GgQQOPejqdDnXr1kVycrLPdpKTk2EwGLTk0CkqKkpbJzk52SOpcy53LmMO/uxLX+s46xS1TmniWttdv34diqIUu98KK+pz7L6fnWX+tskcSrPv/IlhaY6p2k5VVUyfPh19+vRBhw4dAPB3T3XEiR3zsmbNGoSGhmo/e/bs0ZbZbDaMGDECRITly5dXYS9rr+L2P2OMVZWpU6fi2LFjWLt2bVV3hRWDEzvmZejQofj555+1H+dQuDOpO3/+PHbs2KGN1gFAdHQ0UlNTPdqx2+24efMmoqOjfW4nOjoaVqsV6enpHuUpKSnaOtHR0V53VzlfF9VuTVfU/i+OP/vS1zrOOkWtU5q41nb169eHLMvF7rfCivocu+9nZ5m/bTKH0uw7f2JYmmOqNps2bRq2bNmC3bt3IzY2Vivn757qhxM75iUsLAwtWrTQfoKCgrSk7vTp0/jmm29Qr149j3V69+6N9PR0HD58WCvbtWsXVFVFz549fW6na9eu0Ov12Llzp1Z26tQpXLhwAb1799baPXr0qMcB7kwq27VrV55vu9rwtf9L4s++LCwuLg7R0dEe62RmZuLgwYMe+z/QuNZ2BoMBXbt29dhvqqpi586dRe7r3r17e9QHHJ9jZ31/YsF8K82+8yeGpTmmaiMiwrRp07Bx40bs2rULcXFxHsv5u6caquq7N1j1Z7VaaejQoRQbG0s///wzXb16VftxvzP13nvvpS5dutDBgwdp79691LJlS49b0y9dukStW7emgwcPamVPP/00NWnShHbt2kWHDh2i3r17U+/evbXlzulOBg0aRD///DNt376dIiMjb7npTm7cuEGJiYm0detWAkBr166lxMREunr1qlanpH1JRNS6dWvasGGD9nrRokUUERFBmzdvpv/7v/+jYcOG+ZxyoLi43orWrl1LRqORVq1aRSdOnKBJkyZRRESEdvf2mDFj6MUXX9Tq79u3j3Q6Hf31r3+lkydP0rx583xOd1JSLG5VWVlZlJiYSImJiQSA3nnnHUpMTNTuyvdn391zzz30/vvva69LiiGRf8dUbTd58mQym8307bffenz35+bmanX4u6d64cSOlcg5xYavn927d2v1bty4QSNHjqTQ0FAKDw+nxx9/nLKysrzacV8nLy+PpkyZQnXq1KHg4GB64IEHPJIVIqJz587RH/7wBwoKCqL69evTn//8Z7LZbBX9tquVlStX+tz/8+bN0+r4sy8B0MqVK7XXqqrSnDlzKCoqioxGIw0cOJBOnTrlsU5Jcb1Vvf/++9SkSRMyGAzUo0cPOnDggLasf//+NG7cOI/6n3/+ObVq1YoMBgO1b9+etm7d6rHcn1jcqnbv3u3z8+/cx/7su6ZNm3ocL0TFx5DIv2Oqtivqu9/9e4S/e6oXQURUSYODjDHGGGOsAvE1dowxxhhjtQQndowxxhhjtQQndowxxhhjtQQndowxxhhjtQQndowxxhhjtQQndowxxhhjtQQndowxxhhjtQQndowxxhhjtQQndowxVk4++ugjDBo0qMK3s337dnTu3Bmqqlb4thhjNQsndowxVg7y8/MxZ84czJs3r8K3de+990Kv12PNmjUVvi3GWM3CiR1jjJWDL774AuHh4ejTp0+lbG/8+PF47733KmVbjLGagxM7xhhzc+3aNURHR+ONN97Qyn744QcYDAbs3LmzyPXWrl2LIUOGeJQNGDAA06dP9ygbPnw4xo8fr71u1qwZXnvtNYwdOxahoaFo2rQpvvzyS1y7dg3Dhg1DaGgoOnXqhEOHDnm0M2TIEBw6dAhnzpwp/ZtljNU6nNgxxpibyMhIfPzxx3jllVdw6NAhZGVlYcyYMZg2bRoGDhxY5Hp79+5Ft27dSrXNd999F3369EFiYiISEhIwZswYjB07FqNHj8aRI0fQvHlzjB07FkSkrdOkSRNERUVhz549pdomY6x24sSOMcYKue+++zBx4kSMGjUKTz/9NEJCQrBw4cIi66enpyMjIwMxMTGl3t5TTz2Fli1bYu7cucjMzET37t3xyCOPoFWrVpg5cyZOnjyJlJQUj/ViYmJw/vz5Um2TMVY7cWLHGGM+/PWvf4Xdbsf69euxZs0aGI3GIuvm5eUBAEwmU6m21alTJ+33qKgoAEDHjh29ylJTUz3WCwoKQm5ubqm2yRirnTixY4wxH86cOYMrV65AVVWcO3eu2Lr16tWDEAJpaWkltqsoileZXq/XfhdCFFlWeHqTmzdvIjIyssRtMsZuHZzYMcZYIVarFaNHj8ajjz6KV199FU8++aTXaJk7g8GAdu3a4cSJE17LCp8+PXv2bLn0MT8/H2fOnEGXLl3KpT3GWO3AiR1jjBXy8ssvIyMjA++99x5mzpyJVq1a4Yknnih2ncGDB2Pv3r1e5Zs3b8aGDRtw5swZvP766zhx4gTOnz+Py5cvl6mPBw4cgNFoRO/evcvUDmOsduHEjjHG3Hz77bdYsmQJVq9ejfDwcEiShNWrV2PPnj1Yvnx5ketNmDAB27ZtQ0ZGhkd5QkIC3nrrLbRr1w7ff/89li1bhh9//BGrV68uUz8/++wzjBo1CsHBwWVqhzFWuwhyv3+eMcZYqT3yyCO44447MGvWLACOeew6d+6MJUuWlOt2rl+/jtatW+PQoUOIi4sr17YZYzUbj9gxxlg5Wbx4MUJDQyt8O+fOncOyZcs4qWOMeeERO8YYqyAVNWLHGGNF4cSOMcYYY6yW4FOxjDHGGGO1BCd2jDHGGGO1BCd2jDHGGGO1BCd2jDHGGGO1BCd2jDHGGGO1BCd2jDHGGGO1BCd2jDHGGGO1BCd2jDHGGGO1xP8D6n3PUXdKJxoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plotting the simulation\n",
    "sim_base.plot(z=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Meep simulation was built to be equivalent to the Tidy3D one. The geometry was imported from the GDS file generated by the Tidy3D model, and the source and output monitor were positioned exactly as they are in Tidy3D. Since Meep sources are not normalized like those in Tidy3D, an additional monitor was placed immediately after the source in the Meep simulation to normalize the output. A schematic of the Meep model is shown in the figure below, where the red line represents the mode source, the blue lines represent the monitors, and the green area indicates the PML.\n",
    "\n",
    "<img src=\"img/MeepMMIBenchmark.png\" width=\"600\" alt=\"Schematic of the experiment\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Meep resolution was set to 30 steps per $\\mu\\text{m}$, which corresponds to approximately 13 steps per wavelength in silicon. Therefore, we will create and run a new Tidy3D simulation object with this grid specification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">08:29:32 CEST </span>Created task <span style=\"color: #008000; text-decoration-color: #008000\">'MMI'</span> with task_id                                   \n",
       "<span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">              </span><span style=\"color: #008000; text-decoration-color: #008000\">'fdve-771b7520-abc1-42d2-a778-e26c111992bc'</span> and task_type <span style=\"color: #008000; text-decoration-color: #008000\">'FDTD'</span>. \n",
       "</pre>\n"
      ],
      "text/plain": [
       "\u001b[2;36m08:29:32 CEST\u001b[0m\u001b[2;36m \u001b[0mCreated task \u001b[32m'MMI'\u001b[0m with task_id                                   \n",
       "\u001b[2;36m              \u001b[0m\u001b[32m'fdve-771b7520-abc1-42d2-a778-e26c111992bc'\u001b[0m and task_type \u001b[32m'FDTD'\u001b[0m. \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">              </span>View task using web UI at                                         \n",
       "<span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">              </span><a href=\"https://tidy3d.simulation.cloud/workbench?taskId=fdve-771b7520-abc1-42d2-a778-e26c111992bc\" target=\"_blank\"><span style=\"color: #008000; text-decoration-color: #008000\">'https://tidy3d.simulation.cloud/workbench?taskId=fdve-771b7520-ab</span></a>\n",
       "<span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">              </span><a href=\"https://tidy3d.simulation.cloud/workbench?taskId=fdve-771b7520-abc1-42d2-a778-e26c111992bc\" target=\"_blank\"><span style=\"color: #008000; text-decoration-color: #008000\">c1-42d2-a778-e26c111992bc'</span></a>.                                       \n",
       "</pre>\n"
      ],
      "text/plain": [
       "\u001b[2;36m             \u001b[0m\u001b[2;36m \u001b[0mView task using web UI at                                         \n",
       "\u001b[2;36m              \u001b[0m\u001b]8;id=731777;https://tidy3d.simulation.cloud/workbench?taskId=fdve-771b7520-abc1-42d2-a778-e26c111992bc\u001b\\\u001b[32m'https://tidy3d.simulation.cloud/workbench?\u001b[0m\u001b]8;;\u001b\\\u001b]8;id=54548;https://tidy3d.simulation.cloud/workbench?taskId=fdve-771b7520-abc1-42d2-a778-e26c111992bc\u001b\\\u001b[32mtaskId\u001b[0m\u001b]8;;\u001b\\\u001b]8;id=731777;https://tidy3d.simulation.cloud/workbench?taskId=fdve-771b7520-abc1-42d2-a778-e26c111992bc\u001b\\\u001b[32m=\u001b[0m\u001b]8;;\u001b\\\u001b]8;id=414616;https://tidy3d.simulation.cloud/workbench?taskId=fdve-771b7520-abc1-42d2-a778-e26c111992bc\u001b\\\u001b[32mfdve\u001b[0m\u001b]8;;\u001b\\\u001b]8;id=731777;https://tidy3d.simulation.cloud/workbench?taskId=fdve-771b7520-abc1-42d2-a778-e26c111992bc\u001b\\\u001b[32m-771b7520-ab\u001b[0m\u001b]8;;\u001b\\\n",
       "\u001b[2;36m              \u001b[0m\u001b]8;id=731777;https://tidy3d.simulation.cloud/workbench?taskId=fdve-771b7520-abc1-42d2-a778-e26c111992bc\u001b\\\u001b[32mc1-42d2-a778-e26c111992bc'\u001b[0m\u001b]8;;\u001b\\.                                       \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">              </span>Task folder: <a href=\"https://tidy3d.simulation.cloud/folders/9b36e144-ddb6-41f8-8dd8-30b62b26a870\" target=\"_blank\"><span style=\"color: #008000; text-decoration-color: #008000\">'default'</span></a>.                                           \n",
       "</pre>\n"
      ],
      "text/plain": [
       "\u001b[2;36m             \u001b[0m\u001b[2;36m \u001b[0mTask folder: \u001b]8;id=336990;https://tidy3d.simulation.cloud/folders/9b36e144-ddb6-41f8-8dd8-30b62b26a870\u001b\\\u001b[32m'default'\u001b[0m\u001b]8;;\u001b\\.                                           \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9d1e9c29c17d441f80a13461ba7b34f9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">08:29:35 CEST </span>Maximum FlexCredit cost: <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">0.158</span>. Minimum cost depends on task      \n",
       "<span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">              </span>execution details. Use <span style=\"color: #008000; text-decoration-color: #008000\">'web.real_cost(task_id)'</span> to get the billed \n",
       "<span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">              </span>FlexCredit cost after a simulation run.                           \n",
       "</pre>\n"
      ],
      "text/plain": [
       "\u001b[2;36m08:29:35 CEST\u001b[0m\u001b[2;36m \u001b[0mMaximum FlexCredit cost: \u001b[1;36m0.158\u001b[0m. Minimum cost depends on task      \n",
       "\u001b[2;36m              \u001b[0mexecution details. Use \u001b[32m'web.real_cost\u001b[0m\u001b[32m(\u001b[0m\u001b[32mtask_id\u001b[0m\u001b[32m)\u001b[0m\u001b[32m'\u001b[0m to get the billed \n",
       "\u001b[2;36m              \u001b[0mFlexCredit cost after a simulation run.                           \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">08:29:36 CEST </span>status = success                                                  \n",
       "</pre>\n"
      ],
      "text/plain": [
       "\u001b[2;36m08:29:36 CEST\u001b[0m\u001b[2;36m \u001b[0mstatus = success                                                  \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d6b3f436c68d4d74899028be9aa34a57",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">08:29:40 CEST </span>loading simulation from simulation_data.hdf5                      \n",
       "</pre>\n"
      ],
      "text/plain": [
       "\u001b[2;36m08:29:40 CEST\u001b[0m\u001b[2;36m \u001b[0mloading simulation from simulation_data.hdf5                      \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "grid_spec = td.GridSpec.auto(min_steps_per_wvl=13)\n",
    "sim = sim_base.updated_copy(grid_spec=grid_spec)\n",
    "sim_data = web.run(simulation=sim, task_name=\"MMI\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we load the results obtained from running the Meep model on a local machine and compare them with those from Tidy3D. As we can see, the results are very similar. Tidy3D completes the simulation in about 2 minutes, while Meep takes approximately 690 minutes to run on a laptop using 10 i7-1365U processors — a time difference of over 300 times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAG0CAYAAADacZikAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWFFJREFUeJzt3Xd8VHW+//HXzKR30gtJCC0QegtVAi6Kigi6K6gICIrXtcvqqtdVV1bNdfen4qpXFHtZ0V1BXVGuLk2QTuhC6IQWQoD0PnN+fwwMjBQTkskkzPv5eJzHTL6nzOecROfNOd/zPSbDMAxEREREPIjZ3QWIiIiINDYFIBEREfE4CkAiIiLicRSARERExOMoAImIiIjHUQASERERj6MAJCIiIh7Hy90FNEU2m41Dhw4RHByMyWRydzkiIiJSC4ZhUFxcTHx8PGbzhc/xKACdw6FDh0hMTHR3GSIiInIR9u/fT8uWLS+4jALQOQQHBwP2AxgSEuLmakRERKQ2ioqKSExMdHyPX4gC0DmcuuwVEhKiACQiItLM1Kb7ijpBi4iIiMdRABIRERGPowAkIiIiHkd9gEREROrAarVSXV3t7jI8kre3NxaLpUG2pQAkIiJSC4ZhkJubS0FBgbtL8WhhYWHExsbWe5w+BSAREZFaOBV+oqOjCQgI0EC5jcwwDMrKysjLywMgLi6uXttTABIREfkVVqvVEX4iIiLcXY7H8vf3ByAvL4/o6Oh6XQ5TJ2gREZFfcarPT0BAgJsrkVO/g/r2w1IAEhERqSVd9nK/hvodKACJiIiIx1EAEhEREY+jACQiIiIep0kEoNdff51WrVrh5+dH3759WbVq1XmXff/99zGZTE6Tn5+fY351dTWPPvooXbp0ITAwkPj4eCZMmMChQ4caY1dERESalNtuuw2TycRdd9111rx77rkHk8nEbbfd1viFuZnbA9Bnn33G1KlTefrpp8nKyqJbt24MHz7ccZ//uYSEhHD48GHHtG/fPse8srIysrKyePLJJ8nKymL27NlkZ2dz3XXXNcbuiIiINDmJiYnMmjWL8vJyR1tFRQX/+Mc/SEpKcmNl7uP2APTSSy8xZcoUJk2aRFpaGjNmzCAgIIB33333vOuYTCZiY2MdU0xMjGNeaGgoP/zwA2PGjCE1NZV+/frx2muvsXbtWnJychpjl0RExAMYhkFpZblbJsMw6lRrz549SUxMZPbs2Y622bNnk5SURI8ePRxtNpuNzMxMUlJS8Pf3p1u3bvzrX/9y2tbmzZu5+uqrCQoKIiYmhvHjx5Ofn++YP2TIEO69917uvfdeQkNDiYyM5Mknn6xzza7m1oEQq6qqWLt2LY8//rijzWw2M2zYMJYvX37e9UpKSkhOTsZms9GzZ0+ef/55OnXqdN7lCwsLMZlMhIWFnXN+ZWUllZWVjp+LiorqvjMiIuJRyqoqaPPAULd89q5XFhLo61+ndSZPnsx7773HuHHjAHj33XeZNGkSixYtciyTmZnJxx9/zIwZM2jXrh0//vgjt956K1FRUWRkZFBQUMDll1/OHXfcwcsvv0x5eTmPPvooY8aMYcGCBY7tfPDBB9x+++2sWrWKNWvWcOedd5KUlMSUKVMaZP8bglsDUH5+Plar1ekMDkBMTAzbtm075zqpqam8++67dO3alcLCQv7f//t/DBgwgC1bttCyZcuzlq+oqODRRx/l5ptvJiQk5JzbzMzM5Jlnnqn/DomIiDRRt956K48//rij28hPP/3ErFmzHAGosrKS559/nv/85z/0798fgNatW7N06VLefPNNMjIyeO211+jRowfPP/+8Y7vvvvsuiYmJbN++nfbt2wP2S24vv/wyJpOJ1NRUNm3axMsvv6wAVB/9+/d3/GIABgwYQMeOHXnzzTf5y1/+4rRsdXU1Y8aMwTAM3njjjfNu8/HHH2fq1KmOn4uKikhMTGz44kVE5JIR4OPHrlcWuu2z6yoqKooRI0bw/vvvYxgGI0aMIDIy0jF/586dlJWVccUVVzitV1VV5bhMtmHDBhYuXEhQUNBZ29+1a5cjAPXr189pwML+/fvz4osvYrVaG+xp7vXl1gAUGRmJxWLhyJEjTu1HjhwhNja2Vtvw9vamR48e7Ny506n9VPjZt28fCxYsOO/ZHwBfX198fX3rvgMiIuKxTCZTnS9DudvkyZO59957Afsd2GcqKSkBYO7cuSQkJDjNO/UdWVJSwsiRI3nhhRfO2nZ9H07a2NwagHx8fOjVqxfz589n9OjRgL0D1vz58x2/oF9jtVrZtGkT11xzjaPtVPjZsWMHCxcu1IPrREREgKuuuoqqqipMJhPDhw93mpeWloavry85OTlkZGScc/2ePXvyxRdf0KpVK7y8zh8hVq5c6fTzihUraNeuXZM5+wNN4BLY1KlTmThxIr179yY9PZ3p06dTWlrKpEmTAJgwYQIJCQlkZmYCMG3aNPr160fbtm0pKCjgb3/7G/v27eOOO+4A7OHnd7/7HVlZWXzzzTdYrVZyc3MBCA8Px8fHxz07KiIi4mYWi4WtW7c63p8pODiYhx9+mIceegibzcagQYMoLCzkp59+IiQkhIkTJ3LPPfcwc+ZMbr75Zv74xz8SHh7Ozp07mTVrFm+//bZjmzk5OUydOpX/+q//Iisri1dffZUXX3yx0ff3QtwegMaOHcvRo0d56qmnyM3NpXv37sybN8/RMTonJwez+fTd+idOnGDKlCnk5ubSokULevXqxbJly0hLSwPg4MGDfP311wB0797d6bMWLlzIkCFDGmW/REREmqILdQn5y1/+QlRUFJmZmezevZuwsDB69uzJf//3fwMQHx/PTz/9xKOPPsqVV15JZWUlycnJXHXVVU7f1RMmTKC8vJz09HQsFgsPPPAAd955p8v3rS5MRlO7Mb8JKCoqIjQ0lMLCwgv+oYiIiGeoqKhgz549pKSkOD19QM42ZMgQunfvzvTp012y/Qv9Lury/e32gRBFREREGpsCkIiIiHgct/cBEhERkUvHmSNLN2U6AyQiIiIeRwFIREREPI4CkIiIiHgcBSARERHxOApAIiIi4nEUgERERDzYbbfd5nge5/kMGTKEBx98sFHqaSwKQCIiIpcok8l0wenPf/4zr7zyCu+//36Dfu7s2bPp3bs3YWFhBAYG0r17dz766COnZYYMGeKow9fXl4SEBEaOHMns2bMbtJbz0ThAIiIil6jDhw873n/22Wc89dRTZGdnO9qCgoIICgpq8M8NDw/niSeeoEOHDvj4+PDNN98wadIkoqOjnZ5CP2XKFKZNm0ZNTQ0HDhxgzpw53HTTTdx222289dZbDV7XmXQGSERE5BIVGxvrmEJDQzGZTE5tQUFBZ10CKy0tZcKECQQFBREXF3fWU9ynTZtG586dz/qs7t278+STTwL2szvXX389HTt2pE2bNjzwwAN07dqVpUuXOq0TEBBAbGwsLVu2pF+/frzwwgu8+eabzJw5k//85z8Nf0DOoAAkIiJysaw17plc6JFHHmHx4sV89dVXfP/99yxatIisrCzH/MmTJ7N161ZWr17taFu3bh0bN25k0qRJZ23PMAzmz59PdnY2gwcP/tXPnzhxIi1atHD5pTBdAhMREbkY1hrI+ujXl3OFnuPB0vBf4SUlJbzzzjt8/PHH/OY3vwHggw8+oGXLlo5lWrZsyfDhw3nvvffo06cPAO+99x4ZGRm0bt3asVxhYSEJCQlUVlZisVj43//9X6644opfrcFsNtO+fXv27t3bsDv3y89x6dZFRESk2di1axdVVVX07dvX0RYeHk5qaqrTclOmTOHTTz+loqKCqqoq/vGPfzB58mSnZYKDg1m/fj2rV6/mueeeY+rUqbV+TphhGJhMpnrvz4XoDJCIiMjFsHjZz8S467PdaOTIkfj6+jJnzhx8fHyorq7md7/7ndMyZrOZtm3bAvb+QVu3biUzM5MhQ4ZccNtWq5UdO3Y4zi65igKQiIjIxXJzEGlobdq0wdvbm5UrV5KUlATAiRMn2L59OxkZGY7lvLy8mDhxIu+99x4+Pj7cdNNN+Pv7X3DbNpuNysrKX63hgw8+4MSJE/z2t7+t3878ikvrNyciIiIXLSgoiNtvv51HHnmEiIgIoqOjeeKJJzCbz+4xc8cdd9CxY0cAfvrpJ6d5mZmZ9O7dmzZt2lBZWcm3337LRx99xBtvvOG0XFlZGbm5uU63wb/88sv8/ve/Z+jQoa7bURSARERE5Ax/+9vfKCkpYeTIkQQHB/OHP/yBwsLCs5Zr164dAwYM4Pjx4059hsB+K/3dd9/NgQMH8Pf3p0OHDnz88ceMHTvWabmZM2cyc+ZMfHx8iIiIoFevXnz22Wdcf/31Lt1HAJNhGIbLP6WZKSoqIjQ0lMLCQkJCQtxdjoiIuFlFRQV79uwhJSUFPz8/d5fTJBiGQbt27bj77ruZOnVqo33uhX4Xdfn+1hkgERERqZOjR48ya9YscnNzzzn2T3OgACQiIiJ1Eh0dTWRkJG+99RYtWrRwdzkXRQFIRERE6uRS6D2jgRBFRETE4ygAiYiI1NKlcOajuWuo34ECkIiIyK/w9vYG7OPWiHud+h2c+p1cLPUBEhER+RUWi4WwsDDy8vIACAgIcPmzqsSZYRiUlZWRl5dHWFgYFoulXttTABIREamF2NhYAEcIEvcICwtz/C7qQwFIRESkFkwmE3FxcURHR1NdXe3ucjySt7d3vc/8nKIAJCIiUgcWi6XBvoTFfdQJWkRERDyOApCIiIh4HAUgERER8TgKQCIiIuJxFIBERETE4ygAiYiIiMdRABIRERGPowAkIiIiHkcBSERERDyOApCIiIh4HAUgERER8TgKQCIiIuJxFIBERETE4ygAiYiIiMdRABIRERGPowAkIiIiHkcBSERERDyOApCIiIh4HAUgERER8TgKQCIiIuJxFIBERETE4ygAiYiIiMdRABIRERGPowAkIiIiHkcBSERERDyOApCIiIh4HAUgERER8TgKQCIiIuJxFIBERETE4ygAiYiIiMdRABIRERGPowAkIiIiHkcBSERERDxOkwhAr7/+Oq1atcLPz4++ffuyatWq8y77/vvvYzKZnCY/Pz+nZQzD4KmnniIuLg5/f3+GDRvGjh07XL0bIiIi0ky4PQB99tlnTJ06laeffpqsrCy6devG8OHDycvLO+86ISEhHD582DHt27fPaf5f//pX/v73vzNjxgxWrlxJYGAgw4cPp6KiwtW7IyIiIs2A2wPQSy+9xJQpU5g0aRJpaWnMmDGDgIAA3n333fOuYzKZiI2NdUwxMTGOeYZhMH36dP70pz8xatQounbtyocffsihQ4f48ssvG2GPREREpKlzawCqqqpi7dq1DBs2zNFmNpsZNmwYy5cvP+96JSUlJCcnk5iYyKhRo9iyZYtj3p49e8jNzXXaZmhoKH379j3vNisrKykqKnKaRERE5NLl1gCUn5+P1Wp1OoMDEBMTQ25u7jnXSU1N5d133+Wrr77i448/xmazMWDAAA4cOADgWK8u28zMzCQ0NNQxJSYm1nfXREREpAlz+yWwuurfvz8TJkyge/fuZGRkMHv2bKKionjzzTcvepuPP/44hYWFjmn//v0NWLGIiIg0NW4NQJGRkVgsFo4cOeLUfuTIEWJjY2u1DW9vb3r06MHOnTsBHOvVZZu+vr6EhIQ4TSIiInLpcmsA8vHxoVevXsyfP9/RZrPZmD9/Pv3796/VNqxWK5s2bSIuLg6AlJQUYmNjnbZZVFTEypUra71NERERubR5ubuAqVOnMnHiRHr37k16ejrTp0+ntLSUSZMmATBhwgQSEhLIzMwEYNq0afTr14+2bdtSUFDA3/72N/bt28cdd9wB2O8Qe/DBB3n22Wdp164dKSkpPPnkk8THxzN69Gh37aaIiIg0IW4PQGPHjuXo0aM89dRT5Obm0r17d+bNm+foxJyTk4PZfPpE1YkTJ5gyZQq5ubm0aNGCXr16sWzZMtLS0hzL/PGPf6S0tJQ777yTgoICBg0axLx5884aMFFEREQ8k8kwDMPdRTQ1RUVFhIaGUlhYqP5AIiIizURdvr+b3V1gIiIiIvWlACQiIiIeRwFIREREPI4CkIiIiHgcBSARERHxOApAIiIi4nEUgERERMTjuH0gRPEwhgHWKqipPDlV2F9tNYAJTGYwYX91/Hzy1ewF3v7gEwgWbzfviIiINGcKQNKwDAOqSqC8ACoKT74WQU35ycBTBTTA2Jtmb/DxB+8A+3TqvU8Q+IeBbwiYdYJTRETOTQFILl51BZTkQfmJ04GnovDk2ZxfYfYGLx/w8rO/Wnzs4cmwAcbp94Zx+mdrFVSX219t1VBRbQ9X52Iyg1+oPQz5h4FfmIKRiIg4KABJ7VWWQMkRKM6Dklx76DmXX4YP3xDwCbCHHC8/8PIFs+Xi67BWQ3WZPQxVlZ18X2Z/X1lsr8tWczKYnTi7Nv8WEBh5evILUygSEfEwCkByfhVFUJx7cjpiv7T1S36hEBBxOuz4h4FPsGsDhcUbLKH2zz6XMy/DnZoqCk4Ho7Jj9ulotn15s5d9HwIjzghF59m2iIhcEhSA5DSbDUqPQkEOFOy3X85yYrKHhKBoCI6FoBjw9mv0Mg3DoKyqgtKKMkoqSik++VpaWUaN1YqXxYLFbJ+8zBbM5kC8fELw8muFt7WKSIuVaIuBpfw4lB6zX04rOWKfTrH4QEC485ki3+BG31eXs1lPnkGrsL/WVNoDpMmEvRP6qQVNZ7SZ7MfCL1Sd0UWk2VIA8nQ1VVB4wB54Cg/Y+9ecYjJDYBSEnAw7QdEu/8IrrShj//Fc9h87zP5jh8nJP+R4n198nOLyUkory7EZtnp9jo+XN4kRcSRHxJMaHkn7kEBSAgJI8DMT7WXCl6rTZ79OsficDkMBERDQwn55z2Q6/we5k80GVcVQUQyVRfazYlXlpy8fnupPVR8+gWf0s2px8jXU3q9LRKQJMxmG0QC35FxaioqKCA0NpbCwkJCQEHeX0/Cqy+H4Xjix196J+cwwYfGBsEQIbWmfXPRFZhgG+48dZuXODazatYFNOdvJOXaI4yUFtd6GyWQiyDeAIL8AgvwCCfILwMtswWrYqLFasdqsWG02rLYa+8+Gjeqaao4WHafGZj3/djHoER1PRlJr+sQn0LlFOFFe9vazF7aAX8jpL/5Tr36h9evnVBs22xl9oUrs/Z9OhZ3KYnufqNrccWey2M/keQecPKNnwtHxHM74+zijc3pFkX0Ig/PxCbSH5vBWEJIAFv1bS0Rcry7f3wpA53BJBqCaSjixD47vhqJcnL4Y/ULtoScsyX7GxwX9d6w2K9mHdrNi5wZW7dzAyp3rOVxw9JzLhgWEkBgRS2JEnH2KjCcpIo7o0EiC/QIJ9gsk0C+AAB8/TBdx9sVqs3Lw+BFy8g+xL/8ge48edLzm5B+ioOzsO8viQyMZktKeAS2T6RUdS1KAL5bK4gvc8WYCv2Dw8rd3+j51t5uX3y/en3H3m2GzX5I689Wwnr5MVVUO1aVndPwu//WdNXvZL1f5Btvr8Q6wj6V05uvFhtzqipNDHZw4/VpecHZdZi8ITYAWyRCaqLNDIuIyCkD1dMkEoJoqe3+e43ug6JDzmZ7AKPu/zsOS7GcwXCC/+ATfrlvE95uWsmrnBorKnTtRe5ktdE3uQN+23emV0omU6EQSI+II8Q9yST21daykgLW7N7Fix3pW7drIhn1bqbY6B51gv0C6J3cgPbEtPeMS6BwZTbS3BdOpoQDqe2mptkzm04ND+gTZQ45vsP3SnG+w/e67RlZZXkzuoWyMgn14FR3EXF2OYRgY2PsWVQXEUBkcR1VQLDERCcSERl5UkBUR+SUFoHpq1gHIWgOF++2hp+CA/QzCKQHhEJ4CLVq5PPT8O2s+P2VnOfXVCfQNoE+bLqS36Ubftt3okdKJAJ/G70RdV2VVFazbs4VVuzawcucGVu/aRGll2VnLRQa3oFtyR7ond6RHQgrdY+KJ8PPDZK06PeK1Y/Trk23WKsBkv1xmMp8c8dpivyxlMtvPxpksJ0NOAHgHnnwNsL96+TV6HyTDMDhadJx9+QdPTodOnk07RE7+QQ4XHOXM/62Em2pIslSTZK4ixHT678EADtm8yfWJJLxlJ9IS2pLWsh1pCW1pH9cK/2bwtyEiTYsCUD01uwBks9o7MB/fY+/MfOZlGb9Qe+gJT7H3T3GBC4WebskdubbnUDI6ppOW0BavS6AvSI21hq0Hd7F+31Y27NvK+n1b2XpwF9Zz9Cvy8/YlJjSS2LBI+2toFDFhkcSGRhITFkVMaAThQWGE+gc1uWNjGAaHC46y7dAuth3cTfbh3Ww7tJsduXspq7zw5bcAX3/HJUoTJsdriKmalqZKEqgg1FRNUXkpNpuNIsPMNqsfu60+1GDCbDLTJiaRtJbt6NOmK9f1/A3RoRGNtOci0lwpANVTswhANhsUHbR3Zi7Icb7k4hMEESdDT4BrvjQqqiv5dt0iZi37hqXZa88KPSN7Xs7InpeTHJXgks9vasqrKvj5wE7WnwxEG/ZtZUfuXuryn1doQDBhASGEBYYQHhhKWGCI4+fI4BZEBrcgKiT85PtwwgKCMdezv1ZVTTVHCvPJKzxGbuFRDh3PY3vuHrYd3EX24T1nXbY8xWwyE98imuSoBJIj40mKiD/9PjKBiKCwX7+sVVFI5aGN5O1aQe6xQxw6cZT9BfksPVbEmtJqSozTncjNJjODUnsxus+VjOgxhNCAS3BIAhGpNwWgemqyAchmtfflOZFjv4PrzNDjHWAPPS1SICjKZSVsPbiLT5Z+xb9WznPqLOyJoefXVFRXcqTwGLkFeeQW5HOkMJ/cgqP218J8jpxsK64ovajtW8wWIoLCiAxuQXhQGL7ePvh6eePt5Y23xdvx3sfijbeXF2aTmaPFxzly8rNzC/N/9a47i9lCm5gkUuNa0yGhNR3iW5Ma15qkyHh8vBpoSARrNeTvgLytUFGEYRgUV5Syu8pEVqmN2dmbWLN7k2NxHy9vftNpANenX8kVXQbqUpmIOCgA1VOTCkA1VfY+PSdyoPCgfdC+U7z8Tl/eCop2WV+Q0ooyvlr7Hz5e+hVZe7Y42hPCY7l5wLXc2PdqhZ56qLbWUFBaREFZEQWlRZz4xfvjJQUcKykgv/gE+UXHyS8+cc471S6Wt8Xr5GW6KGJCI2kXm0yH+DakxqfQOjoJX+9GumvLMOyXco/8bA/6p/i34KBPFP/amc3sNf8h+9Bux6xA3wCu7p7B7/pexeAOfep9RkxEmjcFoHpyewCqKrP35SnYB0WHne/e8g6w307cItk+GrOLQo9hGKzft5VPln7FnNU/ODr9epktDO82mHGDriOjYzoWV491I+dUVVPN8VOhqPgEx0sKqKyppqqmiuqaaqpqaqiqqaLKWnPy5yqsNhuRwS0cfY9O9UMKDwxtendhlRfAka1wbOfpPm3eARDbia1VFuZkLWL2qv/jwPHTA1W2iUli8pAbGdv/GoL8At1Tt4i4lQJQPTV6AKqpOvmQ0cP2MXrKjjnP9wuFFkkQlmwfhdiFX1alleV8ufp7PvhxNhtzsh3traMTuWXgdYztfw1RIeqMKo2kptL+zLYjP58eX8jiA9GpGFEdWXtwN1+s/D++WDXP0V8pyC+AmweMZPKQ35ESnejG4kWksSkA1ZPLA9BZgec4Z43YGxhpP8sTluSyu7fOtPXgLj5aMod/rvjO0SfF18uHa3sOZdygUfRv16PpnSUQz2GzwrFdkLv59DPqTGaIaAOxnSg1+fL5im95Z+E/2Xlkn322ycSwzgO4Y+gYBndM19+viAdQAKonlwWg/B320/rnCjx+IRAUa3/uVnBcowxgV1ldxdx1C/ngx9ms3LnB0Z4S1ZIJg69nTP8RRASFubwOkVozDPvl4dzNzg+vDW0JUanYguNZnL2GmQs+Y8GW5Y7Z7WJbcfvQG7mx3zUE+vq7oXARaQwKQPXksgB05GfIWWl/7xR4Yu0j+TaCU3175qz6nn+tmue4C8hitnBVt8uYOPgGBqX2VmdSafpK8uDwJvswEKd4B0BUO4hsx66CE7y76J/MWjbX0YetRWAIEy67nklDfkdsmOvulhQR91AAqieXBaDKEvu/Whsx8Jyy/fAe5qz+gS9Xf8+eowcc7fEtohk3aBTjBl6nLwRpnioK7f2E8nfa+wydEhIPUakU+0Uwa/m3vL3wc/blHwTsd76N6n0Fdw27mc6J7d1UuIg0NAWgenL7XWAN5ODxI3y5xh56Nu3f7mj39/FjeNfLuCF9OJd36tfkRiAWuSg2q/2Bv/nb7XdPnuLlB5FtsUa04f+yN/Lm/E+dLvkOTO3Ff/3mZoZ1HqAznyLNnAJQPTXXAFRjrWFjzjaWZq9lweblrNi53jHPy2xhSKd+3NDnSoZ3vYxAv8Z/SKZIo6kosgeh/J3OT6f3C4HgeLYUlzJj5QJmZy12PMKkTUwSd/7mJm7sd02zeEadiJxNAaiemksAstlsbDu0i6XZa1mybTXLd6yjpML5IZ392vXghj5Xcm3PywkPCnVTpSJuYrPZBxI9ut0+kKjTzQcmjtkszMnexPvrV7C7rBIbJsICQhh/2WgmDfkd8S2i3VW5iFwEBaB6aqoBqKqmmp8P7iRrzxaWb1/HT9vXnvUog9CAYAa278mA1F5c3S2DhPAY9xQr0tTUVJ0ceuKwfaTpU7fTY39syardm/ly20a2FxRQapipMHkzqMtgbrn8Znq07eXS8bdEpGEoANVTUwhANpuNnUf2sW7vz6zfu5V1e3/m54M7qKqpdlouwNeffm27Myi1N4M69KJTy3YanVmkNqpK7UHoVCCqLsdm2Pj5wC6WbFvNriOn7y5LjExgQJdB9GifjpdfCJjP8xy0M0OSYdhHsTas9v5Jthr7qO6O91b7GSqzGUwW+7hG5pOvpyazxT7PNxB8Q8E/1D4wqpeviw+OSPOkAFRPjRmADMMgv/g4O3L3sTN3Hzty9/LzwZ1s2Lf1rMtZYL+Nt3urNPq07splHXrTvVUa3urELFJ/5SegOA+qiqGylH2HtrN4/QK27NmA1WrvJxQWGMyA9j3p27a7e8cT8vIFvzB7GPILsb8GRjT63aUiTY0CUD25KgAdPpHHpv3Z7DgZdHaeDD3ne7Clv48fXZNS6d4qjR7JafRolUZSZLxGtBVpREcLj/KPRZ8ze+lsKktPEGiy4Wsxk96mO8O6DCC9TVd8vLw5a3BTTGD2On0Wx2y2/2yy2NvMXvazPIbt9FmiU2eIDNvpyVYDFcVQWQjlhVB99j+MHAIiICzR/uicAD2yRjyPAlA9uSoAPf/lG/x93gdntZtMJpIi4mkbm0zb2GRS41LontyR9nEpukVdpImorK7iyzU/8M7Cf7IxZ5ujvUVgCNf1GsaYftfQM6WT6/+BYq223+VWUWDvx1RRZH94bPkJ5+V8Ak8+TifRPuiqbvEXD6AAVE+uCkBfrvmBV+d9SLvYVo6w0z62FSnRifjrtluRZmPrwV38a+V3fLHy/8gtPOpobx2dyI39rua36VeRFBnfuEVVl0PBASjYZ+/XZKs5Pc/iY39cSIsk+/MF1U9QLlEKQPXUFDpBi0jTZ7VZ+Sl7LZ+v+I656xZSXlXhmNc1qQPdkjuQltCWzontSUto23jjb1lr7B27C3Ls05kjZHv5QXQqRHVolGcOijQmlwegjz76iBkzZrBnzx6WL19OcnIy06dPJyUlhVGjRl104U2FApCI1FVpRRnfrl/MP1d8y5LsNfzyf60mk4lWkQl0SmxP55btSGvZjnaxyVjMFnvvoZPLGxgYxskJ+40SZVXllFSUUVJRRmml8+up9yZMBPsHEhoQTLDfyVf/IEL9AgmjihbVBYSU5+FrO3knqckM4SkQkwaBkY16rERcxaUB6I033uCpp57iwQcf5LnnnmPz5s20bt2a999/nw8++ICFCxfWq/imQAFIROrj8Ik8Vu3ayM8HdrD5wA627N/hdKnMXcwmGBwbwxVRLejcIpSkyARiQiMxB8fYg1BYsvoKSbPm0gCUlpbG888/z+jRowkODmbDhg20bt2azZs3M2TIEPLz8+tVfFOgACQiDS2/+IQ9EO3fwZYD29lyYCf7Tj2Y+GTHaZPJhAmTYzgh+3sTAT7+BPkFEOgXQJBfAEG+Aad/9g0gyC8Qm2GjuLyEwvIS+2tZCUXlZ0xlxdScfOwHQISphg5elaT6mkiOjCM5MoGWMa1JTLuc0OReGmtImqW6fH/X+RajPXv20KNHj7PafX19KS0trevmREQ8QmRwCwZ3TGdwx3S3fP6pMcey9mxh7Z4trNm9iXX7tpJVUkr7ilzaH96H76af4D8fERcRT1Lab+h/2ThaxqS4pV4RV6tzAEpJSWH9+vUkJyc7tc+bN4+OHTs2WGEiItJwTCYTUSERDO82mOHdBgP2ByhvO7SbtXs2k7VrIwf2rCaoYA8cO8ThJR+xbMlHmGI6ktZrFFf3uUaP1pFLSp0D0NSpU7nnnnuoqKjAMAxWrVrFp59+SmZmJm+//bYrahQRERfwsnjRObE9nRPbM3HwDQAcLTrG4uVfsGvd1+Tn7sA4spWfv93K1/9+Gb/EnlzR+2pG9rqc2LAoN1cvUj8XdRfYJ598wp///Gd27doFQHx8PM888wy33357gxfoDuoDJCIC+Qc2s275p+zYsZI9efuxGrDX5sNmqz+d2vRi4uAbuKbHkJMjYYu4X6ONA1RWVkZJSQnR0dEXu4kmSQFIROQMxbkc3/4jW35ezMacrezOO8g+mzc/1/hhCY7m1kGjuHXQaF0iE7dzaQDas2cPNTU1tGvXzql9x44deHt706pVqzoX3NQoAImInEPJUTi0nhOHfmb1ro2s2LGeHaUVbLH6cQQ/hncdxKQhv2NQam89s1DcwqUBKCMjg8mTJzNx4kSn9o8//pi3336bRYsW1bngpkYBSETkAkqPQe5mavJ3smX/dpZtz2Jt7kF+rvFlr82H1jGtuC3jBsb0H0GIf5C7qxUP4tIAFBISQlZWFm3btnVq37lzJ71796agoKDOBTc1CkAiIrVQWQJ5P0Pedg4fP8Ty7etYuvtnNlSY2GH1weITwPhBo7ln+K3EhGq0aXE9l44DZDKZKC4uPqu9sLAQq9V6jjVEROSS5BsEiekQ1524o9ncEJXENT0yWLt7M0t2rGfpsePMWvgJHy6Zw/jLRnPv8PEKQtJk1PkM0MiRI/H39+fTTz/FYrE/UdhqtTJ27FhKS0v57rvvXFJoY9IZIBGRi2CzwfFdkLsFo+w423P38v3GpazLy2W31YdcczC/G3Q99w4fr9voxSVcegns559/ZvDgwYSFhXHZZZcBsGTJEoqKiliwYAGdO3e++MqbCAUgEZF6KjwIR37GKDzIjsO7+X7jUvYcPUiuzYsDpiAyBvyOu6+6TUFIGpTLb4M/dOgQr732Ghs2bMDf35+uXbty7733Eh4eftFFNyUKQCIiDaS6Ao7vwcjfyY5dq/lh01L25B3ECuSaAujRayTjR9xNXHisuyuVS0CjjQN0qVIAEhFxgYoijPydbNv8HxasnceevIP2dos3vbsNY+SQWwmMSgH/cMcDYkXqwuUBqKCggFWrVpGXl4fNZnOaN2HChLpurslRABIRcS2jOI8N679l4Y+fcDAvB4AgvwCu7HYZfVP7YAmJh+BYCIlVIJJac2kA+ve//824ceMoKSkhJCTEabArk8nE8ePHL67qJkQBSESkcRhWKwtWf8MHc/+XmoIDRJutxIe2YESPy0lLaGP/jrH4QHAMBEaCXxgEtADfEIUiOYtLA1D79u255ppreP755wkICKhXoU2VApCISOOqttbw4Y9zeOmbmZjKjhNtruE3SSlM7DmAlmHnuHXeZAH/MPBv4fzqq4EXPZlLA1BgYCCbNm2idevW9SqyKVMAEhFxj8KyYv4+7wPeXvA5lTVVmE0wsVcGD152DTHeZig/ARWFYKs59wYsPvYg5Bd2MhiF2cORz6X5D3Zx5tIAdMMNN3DTTTcxZsyYehXZlCkAiYi4V07+ITK/msGc1d8D4Ofty20Zv+XuK8YRHRIOlcX2MFR+AspOQEUBlBcC5/lKs/iccbYoDAKjICACzOZG2iNpDC4NQO+88w7Tpk1j0qRJdOnSBW9vb6f51113Xd0rbmIUgEREmoasPVuY9sWrrNi5HrAHoVsHjeKeK28lrkW088I2q/3sUHmBPRCdCkYVxZwzGJm97P2KgmPsHa4Do8DiffZy0my4NACZL5CWTSbTJfE4DAUgEZGmwzAM5m9ezsvfvsvaPZsB8PHy5paB13Hv8PG0/LUxhGzWk2eITk5lx6EkD6xVv1jQBIER9kAUdHLy9nPBHomraBygelIAEhFpegzDYMm21bw49x1W7twAgLfFi7H9R3DfVRNJjoyvy8bsoag4D4pzoeQIVJX+YiEThMZDRFsISwJLnR+fKY1MAaieFIBERJouwzBYtj2Ll759l5+y1wJgMVv4Xd+ruOfKW2kfl3JxG64ssQeh4jwoybWfLTrF4gMtWkFkG/uZId2C3yS5PACVlpayePFicnJyqKpyPoV4//3312lbr7/+On/729/Izc2lW7duvPrqq6Snp//qerNmzeLmm29m1KhRfPnll472kpISHnvsMb788kuOHTtGSkoK999/P3fddVeta1IAEhFpHlbuXM/L377Hop9XOtraxCQxrMtArug8kPS23fDxush+PRWFcGwX5O+CqpLT7T6BENkWItqAX2g990AakksD0Lp167jmmmsoKyujtLSU8PBw8vPzCQgIIDo6mt27d9d6W5999hkTJkxgxowZ9O3bl+nTp/PPf/6T7OxsoqOjz7ve3r17GTRoEK1btyY8PNwpAN15550sWLCAt99+m1atWvH9999z9913M3v27Fp30FYAEhFpXrL2bGb6d++zYPNyamyn+6IG+QUwJK0fwzoP4Ded+xMVElH3jRuG/cxQ/i44vgds1afnBUZBdCqEt9EdZU2ASwPQkCFDaN++PTNmzCA0NJQNGzbg7e3NrbfeygMPPMANN9xQ62317duXPn368NprrwFgs9lITEzkvvvu47HHHjvnOlarlcGDBzN58mSWLFlCQUGBUwDq3LkzY8eO5cknn3S09erVi6uvvppnn332nNusrKyksrLS8XNRURGJiYkKQCIizUxReQmLf17Ffzb/xPzNy8gvPuE0v3tyGsO6DKBLYnuSo1qSHBmPv08dOjpba6Agx35mqPAgjrvLfIIgrqv9zJDZ0nA7JHXi0gAUFhbGypUrSU1NJSwsjOXLl9OxY0dWrlzJxIkT2bZtW622U1VVRUBAAP/6178YPXq0o33ixIkUFBTw1VdfnXO9p59+mo0bNzJnzhxuu+22swLQnXfeybp16/jyyy+Jj49n0aJFXHfddcydO5fBgwefc5t//vOfeeaZZ85qVwASEWm+bDYb6/dt5T+bl/GfTT+xMefc30+xoVG0ikogOSqBVlEJtIpqSXJkAq2jEwkLvMB3QFUZHNsJuVugpsLe5h1gD0JR7RWE3KAuAajOXdq9vb0dt8JHR0eTk5NDx44dCQ0NZf/+/bXeTn5+PlarlZiYGKf2mJiY84aopUuX8s4777B+/frzbvfVV1/lzjvvpGXLlnh5eWE2m5k5c+Z5ww/A448/ztSpUx0/nzoDJCIizZfZbKZnSid6pnTijyOncKQwn/mbl/Hj1tXsOrKPvUcPUlxRSm7hUXILjzrGGjpTbGgUHRJakxqXQmp8GzrE298H+gXYR5eO6wrRaZC/HQ5vhOoyyFkBhzdAbBeIStXdY01UnX8rPXr0YPXq1bRr146MjAyeeuop8vPz+eijj+jcubMragSguLiY8ePHM3PmTCIjz/FcmJNeffVVVqxYwddff01ycjI//vgj99xzD/Hx8QwbNuyc6/j6+uLr6+uq0kVEpAmICY3kloHXcctAe39QwzA4XlrIvqMH2Hv0IHuPHmRf/kH2nXx/KhjlFh516mQNkBgRZw9D8a3p364Hl3Xog09UKuTvsAehqlLYv+pkEOoM0R01yGITU+dLYGvWrKG4uJihQ4eSl5fHhAkTWLZsGe3ateOdd96he/futdpOXS+BrV+/nh49emCxnD6laLPZAHvKz87OJj4+ntDQUObMmcOIESMcy91xxx0cOHCAefPm1ao2dYIWEZHi8lK2H97DtkO72HZoN9mHdrPt0G7yio6dtWyIfxBXdh3ENd2HMCQtnYCiA3Bow+m7xyw+EN/dHoTUWdplXHoJrHfv3o730dHRtQ4Vv+Tj40OvXr2YP3++IwDZbDbmz5/Pvffee9byHTp0YNOmTU5tf/rTnyguLuaVV14hMTGRiooKqqurzxqt2mKxOMKSiIhIbQT7B9KrdWd6tXa+unGspIDsk4Fo8/7t/LDpJ/KKjvGvlfP418p5+Pv48ZvOA7i2ewbD4lMJOrEDKorsZ4SOboOWfaBFkpv2Sk6pcwC6/PLLmT17NmFhYU7tRUVFjB49mgULFtR6W1OnTmXixIn07t2b9PR0pk+fTmlpKZMmTQJgwoQJJCQkkJmZiZ+f31mX2E7VcKrdx8eHjIwMHnnkEfz9/UlOTmbx4sV8+OGHvPTSS3XdVRERkbNEBIUxoH1PBrTvCdj/8b5m9ybmrlvIN+sWcfB4Lt9kLeCbrAX4eHmT0aEPN6d24DdhfvhWFMHO+RASB4npEBDu5r3xXHUOQIsWLTpr8EOAiooKlixZUqdtjR07lqNHj/LUU0+Rm5tL9+7dmTdvnqNjdE5OzgWfPXYus2bN4vHHH2fcuHEcP36c5ORknnvuuToNhCgiIlJbZrOZ9LbdSG/bjT//7gE25GxjbtZC5q5byO68/fyweRk/bF5GVEAQD/boxfVJCYQDbPnafrdYQg/w9nf3bnicWvcB2rhxIwDdu3dnwYIFhIefTq1Wq5V58+bx5ptvsnfvXpcU2pjUB0hEROrLMAy2HdrN3HUL+eeK79iXfxCAEDPc1aEd17dpT6uoBExevva7yWLSdOt8PblkHCCz2Yzp5LNPzrWKv78/r776KpMnT76IkpsWBSAREWlIVpuVHzb+xMyFnzmeXxZlquHG+HBGpHahW1JHvALCoGVvCG/l1lqbM5cEoH379mEYBq1bt2bVqlVERUU55vn4+BAdHe10h1ZzpgAkIiKu8vOBHcxc8DmzV/0flTWVtDZXkRHixeXtu9GvbXeCo9vY+wcFXsRjOzycS0eCXrx4MQMHDsTLy7n7kNVq5aeffrrggIPNhQKQiIi4Wn7xCT5aMof3F8/mWOFROnlV0N3HyoB23chI60tYUk9I6GkfcFFqxaUByGKxcPjw4bMeVnrs2DGio6OxWq3nWbP5UAASEZHGUlVTzTdZC3hz/ix27NtCD+9y2nlb6dOmK0O7DCK87WUQ00kjSteCS8cBMgzD0RfoTMeOHSMwMLCumxMREfFoPl7e3JA+nOv7XMnirat4+dt3+WbXWnKzN7Jq5wZ6t/6RId1/Q1TaFRDR2t3lXjJqHYBOPeXdZDJx2223OT06wmq1snHjRgYMGNDwFYqIiHgAk8nEkLS+ZHRMZ9n2LF6e+y7bdi6neOcmVu/eRI+s7xnS51riul4LQVG/vkG5oFoHoNDQUMB+Big4OBh//9NjFvj4+NCvXz+mTJnS8BWKiIh4EJPJxMDUXgxM7cWqnRt4Ze477Ni+mOrdm8nas4XuK78mY+AYWna+EgLUUfpi1bkP0DPPPMPDDz98SV/uUh8gERFpSrL2bOF/577F8ewFtDZXYTKZ6NOmC5cPGENkagYERf/6RjyASztBewIFIBERaYo279/OG1+/xpFtC2hlrsbLYmFgai+G9h1NUOv+EBLv7hLdqsEDUM+ePZk/fz4tWrSgR48e5+wEfUpWVlbdK25iFIBERKQpW7N7Ey/Nnk753lWkWKoI9PFlSKd+DOwxHL8kz33YaoPfBTZq1ChHp+dTT24XERER9+jduguf/OFt5m9ezktz/o4l72dK1i1m6bY1XNF1EOmdM/BK6AEtWkEdn6npKXQJ7Bx0BkhERJoLm83G7NX/x/Sv3yC4cB+plkriQsIY3m0wXVO6Yo5obb99PigGLnAF51LQaH2ASkpKsNlsTm2XQmBQABIRkeamsrqKj5Z+yWvfvkNUWS7tvSppFxHDiB5DaRebDN4B9iAU0fqSvXvMpQFoz5493HvvvSxatIiKigpH+6kBEjUStIiIiPuUVJQy4z+f8sYPnxBcXUSKuYorE5MY1WMw8S1O3i3mFwoRbSA8Bfwune85lwaggQMHYhgGDzzwADExMWd1iM7IyKh7xU2MApCIiDR3R4uO8/K37/Lhj3Ow2WpIsNQwIa0rt3TqQURg8OkFAyIgOAYCo+wDLPoGn3+jTZxLA1BQUBBr164lNTW1XkU2ZQpAIiJyqdh79AD/89WbfLnmBwACvSw80G8ot3XpSUhVEfCLGODtfzIMRdsDUUBks3kOmUsD0NChQ3niiScYNmxYvYpsyhSARETkUrNh3zaenfMaS7atASDYL5CHrriJSb0G4l9VDCV5UHYcDOe+vZjM4N/CfqbIPxT8wuyvPkFNrlO1SwPQrl27uOuuu7j11lvp3Lkz3t7eTvO7du1a94qbGAUgERG5FBmGweKtq3h2zuts3r8dgOiQCH5/xS1MuOx6Ar19oOwYlB61B6KSPKguP/fGzF72vkR+oSeDUSj4h4FvqNtuvXdpAFqxYgW33HILe/fuPb0Rk0mdoEVERJoJm83Gl2t+IPOrGew/dhiAsIAQ7rh8DLcPvZEWgaGnF64ssQei8hNQXggVBVBRdPaZolPMXvazRUGR9stnjdivyKUBKC0tjY4dO/LHP/7xnJ2gk5OT615xE6MAJCIinqCqppovVv0fr/3fh+w6kgNAgK8/Ey67nruG3Uxs2HmeOm+zQVXx6UBUXggVhVBeALbqs5f38rX3KwqMPP3q7dfg++PSABQYGMiGDRto27ZtvYpsyhSARETEk1htVuauW8Tf533guDTm4+XNTf2v5Z4rbyU5KqF2GzIMexAqPQqlJy+lnatfEUBMGiT1bcC9cHEAGjlyJLfddhu//e1v61VkU6YAJCIinsgwDBZsWc4r373Pql0bATCbzFzf5wruGnYLXZIu4g5wm9UegkrzTwajfHtISkyH2E4NWr9LA9Bbb73Fs88+y+TJk+nSpctZnaCvu+66ulfcxCgAiYiIp1uxYz2vzHufhVtWONo6tWzHmH7XcEP6lUSF1GM06Zoq+6uXTz2rdObSAGS+QM9udYIWERG5tGzM2cbr33/Md+sXU1Vj799jMVu4vFM/xvYfwRVdBuHr3bBB5mI12rPALlUKQCIiIs5OlBby1Zr/8PmKb8nas8XRHhYQwug+VzC2/wi6J3c86+aoxtToAaigoICwsLD6bqbJUAASERE5vx25e/l8+bf8a+V3HC446mhvF9uKYV0G0qd1F3q37kJ0aOM+dNWlAeiFF16gVatWjB07FoAbb7yRL774gri4OL799lu6det28ZU3EQpAIiIiv85qs7J02xo+Wz6X79Yvpry60ml+UmS8Iwz1adOVDvGt8XLhYzVcGoBSUlL45JNPGDBgAD/88ANjxozhs88+4/PPPycnJ4fvv/++XsU3BQpAIiIidVNcXsr/bfyRVbs2smbXJrYe2sUvI0aArz89W3WiT5suDOs8kF6tOzdoDXX5/q5zDMvNzSUxMRGAb775hjFjxnDllVfSqlUr+vZt2Pv5RUREpHkI9g/kd32v5nd9rwbsgShrz2ZW797E2t2bWbN7E8UVpSzNXsPS7DVUVlc1eACqizoHoBYtWrB//34SExOZN28ezz77LGAfO+BSuANMRERE6i/YP5CMtL5kpNlPjthsNrIP72HN7k2s2b2JoZ36ubW+OgegG264gVtuuYV27dpx7Ngxrr7anvTWrVt3SY8OLSIiIhfPbDbTMaENHRPaMP6y0e4up+4B6OWXX6ZVq1bs37+fv/71rwQFBQFw+PBh7r777gYvUERERKShaRygc1AnaBERkebHpZ2gAXbs2MHChQvJy8vDZnN+wNlTTz11MZsUERERaTR1DkAzZ87k97//PZGRkcTGxjqN+GgymRSAREREpMmrcwB69tlnee6553j00UddUY+IiIiIy53/yabnceLECW688UZX1CIiIiLSKOocgG688cZLYrRnERER8Vx1vgTWtm1bnnzySVasWEGXLl3w9vZ2mn///fc3WHEiIiIirnBRzwI778ZMJnbv3l3votxNt8GLiIg0Py69DX7Pnj0XXZiIiIhIU1DnPkAiIiIizd1FDYR44MABvv76a3JycqiqqnKa99JLLzVIYSIiIiKuUucANH/+fK677jpat27Ntm3b6Ny5M3v37sUwDHr27OmKGkVEREQaVJ0vgT3++OM8/PDDbNq0CT8/P7744gv2799PRkaGxgcSERGRZqHOAWjr1q1MmDABAC8vL8rLywkKCmLatGm88MILDV6giIiISEOrcwAKDAx09PuJi4tj165djnn5+fkNV5mIiIiIi9S5D1C/fv1YunQpHTt25JprruEPf/gDmzZtYvbs2fTr188VNYqIiIg0qDoHoJdeeomSkhIAnnnmGUpKSvjss89o166d7gATERGRZqFOAchqtXLgwAG6du0K2C+HzZgxwyWFiYiIiLhKnfoAWSwWrrzySk6cOOGqekRERERcrs6doDt37nxJPO9LREREPFedA9Czzz7Lww8/zDfffMPhw4cpKipymkRERESaulo/DX7atGn84Q9/IDg4+PTKJpPjvWEYmEwmrFZrw1fZyPQ0eBERkeanLt/ftQ5AFouFw4cPs3Xr1gsul5GRUftKmygFIBERkeanLt/ftb4L7FROuhQCjoiIiHi2OvUBOvOSl4iIiEhzVadxgNq3b/+rIej48eP1KkhERETE1eoUgJ555hlCQ0NdVYuIiIhIo6hTALrpppuIjo52VS0iIiIijaLWfYDU/0dEREQuFbUOQLW8W15ERESkyav1JTCbzebKOkREREQaTZ0fhdHQXn/9dVq1aoWfnx99+/Zl1apVtVpv1qxZmEwmRo8efda8rVu3ct111xEaGkpgYCB9+vQhJyengSsXERGR5sqtAeizzz5j6tSpPP3002RlZdGtWzeGDx9OXl7eBdfbu3cvDz/8MJdddtlZ83bt2sWgQYPo0KEDixYtYuPGjTz55JP4+fm5ajdERESkman1ozBcoW/fvvTp04fXXnsNsF9mS0xM5L777uOxxx475zpWq5XBgwczefJklixZQkFBAV9++aVj/k033YS3tzcfffRRreuorKyksrLS8XNRURGJiYl6FIaIiEgzUpdHYbjtDFBVVRVr165l2LBhp4sxmxk2bBjLly8/73rTpk0jOjqa22+//ax5NpuNuXPn0r59e4YPH050dDR9+/Z1CkjnkpmZSWhoqGNKTEy86P0SERGRps9tASg/Px+r1UpMTIxTe0xMDLm5uedcZ+nSpbzzzjvMnDnznPPz8vIoKSnhf/7nf7jqqqv4/vvvuf7667nhhhtYvHjxeWt5/PHHKSwsdEz79++/+B0TERGRJq9OAyG6U3FxMePHj2fmzJlERkaec5lTd6qNGjWKhx56CIDu3buzbNkyZsyYcd4Hufr6+uLr6+uawkVERKTJcVsAioyMxGKxcOTIEaf2I0eOEBsbe9byu3btYu/evYwcOdLRdirweHl5kZ2dTWJiIl5eXqSlpTmt27FjR5YuXeqCvRAREZHmyG2XwHx8fOjVqxfz5893tNlsNubPn0///v3PWr5Dhw5s2rSJ9evXO6brrruOoUOHsn79ehITE/Hx8aFPnz5kZ2c7rbt9+3aSk5Ndvk8iIiLSPLj1EtjUqVOZOHEivXv3Jj09nenTp1NaWsqkSZMAmDBhAgkJCWRmZuLn50fnzp2d1g8LCwNwan/kkUcYO3YsgwcPZujQocybN49///vfLFq0qLF2S0RERJo4twagsWPHcvToUZ566ilyc3Pp3r078+bNc3SMzsnJwWyu20mq66+/nhkzZpCZmcn9999PamoqX3zxBYMGDXLFLoiIiEgz5NZxgJqquowjICIiIk1DsxgHSERERMRdFIBERETE4ygAiYiIiMdRABIRERGPowAkIiIiHkcBSERERDyOApCIiIh4HAUgERER8TgKQCIiIuJxFIBERETE4ygAiYiIiMdRABIRERGPowAkIiIiHkcBSERERDyOApCIiIh4HAUgERER8TgKQCIiIuJxFIBERETE4ygAiYiIiMdRABIRERGPowAkIiIiHkcBSERERDyOApCIiIh4HAUgERER8TgKQCIiIuJxFIBERETE4ygAiYiIiMdRABIRERGPowAkIiIiHkcBSERERDyOApCIiIh4HAUgERER8TgKQCIiIuJxFIBERETE4ygAiYiIiMdRABIRERGPowAkIiIiHkcBSERERDyOApCIiIh4HAUgERER8TgKQCIiIuJxFIBERETE4ygAiYiIiMdRABIRERGPowAkIiIiHkcBSERERDyOApCIiIh4HAUgERER8TgKQCIiIuJxFIBERETE4ygAiYiIiMdRABIRERGPowAkIiIiHkcBSERERDyOApCIiIh4HAUgERER8TgKQCIiIuJxFIBERETE4ygAiYiIiMdRABIRERGPowAkIiIiHkcBSERERDyOApCIiIh4nCYRgF5//XVatWqFn58fffv2ZdWqVbVab9asWZhMJkaPHn3eZe666y5MJhPTp09vmGJFRESk2XN7APrss8+YOnUqTz/9NFlZWXTr1o3hw4eTl5d3wfX27t3Lww8/zGWXXXbeZebMmcOKFSuIj49v6LJFRESkGXN7AHrppZeYMmUKkyZNIi0tjRkzZhAQEMC777573nWsVivjxo3jmWeeoXXr1udc5uDBg9x333188skneHt7X7CGyspKioqKnCYRERG5dLk1AFVVVbF27VqGDRvmaDObzQwbNozly5efd71p06YRHR3N7bfffs75NpuN8ePH88gjj9CpU6dfrSMzM5PQ0FDHlJiYWPedERERkWbDrQEoPz8fq9VKTEyMU3tMTAy5ubnnXGfp0qW88847zJw587zbfeGFF/Dy8uL++++vVR2PP/44hYWFjmn//v213wkRERFpdrzcXUBdFBcXM378eGbOnElkZOQ5l1m7di2vvPIKWVlZmEymWm3X19cXX1/fhixVREREmjC3BqDIyEgsFgtHjhxxaj9y5AixsbFnLb9r1y727t3LyJEjHW02mw0ALy8vsrOzWbJkCXl5eSQlJTmWsVqt/OEPf2D69Ons3bvXNTsjIiIizYZbA5CPjw+9evVi/vz5jlvZbTYb8+fP59577z1r+Q4dOrBp0yantj/96U8UFxfzyiuvkJiYyPjx4536FAEMHz6c8ePHM2nSJJfti4iIiDQfbr8ENnXqVCZOnEjv3r1JT09n+vTplJaWOsLKhAkTSEhIIDMzEz8/Pzp37uy0flhYGICjPSIigoiICKdlvL29iY2NJTU11fU7JCIiIk2e2wPQ2LFjOXr0KE899RS5ubl0796defPmOTpG5+TkYDa7/W59ERERuYSYDMMw3F1EU1NUVERoaCiFhYWEhIS4uxwRERGphbp8f+vUioiIiHgcBSARERHxOApAIiIi4nEUgERERMTjKACJiIiIx1EAEhEREY+jACQiIiIeRwFIREREPI4CkIiIiHgcBSARERHxOApAIiIi4nEUgERERMTjKACJiIiIx1EAEhEREY+jACQiIiIeRwFIREREPI4CkIiIiHgcBSARERHxOApAIiIi4nEUgERERMTjKACJiIiIx1EAEhEREY+jACQiIiIeRwFIREREPI4CkIiIiHgcBSARERHxOApAIiIi4nEUgERERMTjKACJiIiIx1EAEhEREY+jACQiIiIeRwFIREREPI4CkIiIiHgcBSARERHxOApAIiIi4nEUgERERMTjKACJiIiIx1EAEhEREY+jACQiIiIeRwFIREREPI4CkIiIiHgcBSARERHxOApAIiIi4nEUgERERMTjKACJiIiIx1EAEhEREY+jACQiIiIex8vdBTRFhmEAUFRU5OZKREREpLZOfW+f+h6/EAWgcyguLgYgMTHRzZWIiIhIXRUXFxMaGnrBZUxGbWKSh7HZbBw6dIjg4GBMJlOd1y8qKiIxMZH9+/cTEhLiggoFdJwbi45z49Bxbhw6zo3DXcfZMAyKi4uJj4/HbL5wLx+dAToHs9lMy5Yt672dkJAQ/QfWCHScG4eOc+PQcW4cOs6Nwx3H+dfO/JyiTtAiIiLicRSARERExOMoALmAr68vTz/9NL6+vu4u5ZKm49w4dJwbh45z49BxbhzN4TirE7SIiIh4HJ0BEhEREY+jACQiIiIeRwFIREREPI4CkIiIiHgcBaBf8eOPPzJy5Eji4+MxmUx8+eWXF1x+0aJFmEyms6bc3Fyn5V5//XVatWqFn58fffv2ZdWqVS7ci6bPFcc5MzOTPn36EBwcTHR0NKNHjyY7O9vFe9K0uerv+ZT/+Z//wWQy8eCDDzZ88c2Iq47zwYMHufXWW4mIiMDf358uXbqwZs0aF+5J0+aK42y1WnnyySdJSUnB39+fNm3a8Je//KVWz5a6VNX1OANUVlbyxBNPkJycjK+vL61ateLdd991Wuaf//wnHTp0wM/Pjy5duvDtt9+6aA/OTQHoV5SWltKtWzdef/31Oq2XnZ3N4cOHHVN0dLRj3meffcbUqVN5+umnycrKolu3bgwfPpy8vLyGLr/ZcMVxXrx4Mffccw8rVqzghx9+oLq6miuvvJLS0tKGLr/ZcMVxPmX16tW8+eabdO3ataHKbbZccZxPnDjBwIED8fb25rvvvuPnn3/mxRdfpEWLFg1dfrPhiuP8wgsv8MYbb/Daa6+xdetWXnjhBf7617/y6quvNnT5zcbFHOcxY8Ywf/583nnnHbKzs/n0009JTU11zF+2bBk333wzt99+O+vWrWP06NGMHj2azZs3u2IXzs2QWgOMOXPmXHCZhQsXGoBx4sSJ8y6Tnp5u3HPPPY6frVarER8fb2RmZjZQpc1bQx3nX8rLyzMAY/HixfUr8BLRkMe5uLjYaNeunfHDDz8YGRkZxgMPPNBgdTZ3DXWcH330UWPQoEENW9wlpKGO84gRI4zJkyc7td1www3GuHHjGqDK5q82x/m7774zQkNDjWPHjp13mTFjxhgjRoxwauvbt6/xX//1Xw1RZq3oDJCLdO/enbi4OK644gp++uknR3tVVRVr165l2LBhjjaz2cywYcNYvny5O0pt1s53nM+lsLAQgPDw8MYo7ZLya8f5nnvuYcSIEU5/11J3FzrOX3/9Nb179+bGG28kOjqaHj16MHPmTDdV2rxd6DgPGDCA+fPns337dgA2bNjA0qVLufrqq91RarN06m/1r3/9KwkJCbRv356HH36Y8vJyxzLLly8/6/8Xw4cPb9TvQT0MtYHFxcUxY8YMevfuTWVlJW+//TZDhgxh5cqV9OzZk/z8fKxWKzExMU7rxcTEsG3bNjdV3fz82nH+JZvNxoMPPsjAgQPp3LmzGypunmpznGfNmkVWVharV692c7XNV22O8+7du3njjTeYOnUq//3f/83q1au5//778fHxYeLEiW7eg+ahNsf5scceo6ioiA4dOmCxWLBarTz33HOMGzfOzdU3H7t372bp0qX4+fkxZ84c8vPzufvuuzl27BjvvfceALm5uef8Hjxf/0KXaLRzTZcAanHq71wGDx5s3HrrrYZhGMbBgwcNwFi2bJnTMo888oiRnp7eEGU2ew1xnH/prrvuMpKTk439+/fXs7pLR0Mc55ycHCM6OtrYsGGDY74ugTlrqL9nb29vo3///k7L3HfffUa/fv3qW+IloaGO86effmq0bNnS+PTTT42NGzcaH374oREeHm68//77DVht81Wb43zFFVcYfn5+RkFBgaPtiy++MEwmk1FWVmYYhv3v+R//+IfTeq+//roRHR3d4DWfjy6BNYL09HR27twJQGRkJBaLhSNHjjgtc+TIEWJjY91R3iXjzON8pnvvvZdvvvmGhQsX0rJlSzdUdmk58zivXbuWvLw8evbsiZeXF15eXixevJi///3veHl5YbVa3Vxt8/XLv+e4uDjS0tKclunYsSM5OTmNXdol5ZfH+ZFHHuGxxx7jpptuokuXLowfP56HHnqIzMxMN1bZvMTFxZGQkEBoaKijrWPHjhiGwYEDBwCIjY11+/egAlAjWL9+PXFxcQD4+PjQq1cv5s+f75hvs9mYP38+/fv3d1eJl4QzjzOAYRjce++9zJkzhwULFpCSkuLG6i4dZx7n3/zmN2zatIn169c7pt69ezNu3DjWr1+PxWJxc7XN1y//ngcOHHjWMA7bt28nOTm5sUu7pPzyOJeVlWE2O381WiwWbDZbY5fWbA0cOJBDhw5RUlLiaNu+fTtms9nxj9D+/fs7fQ8C/PDDD436Pag+QL+ipKTE6V8He/bsYf369YSHh5OUlMTjjz/OwYMH+fDDDwGYPn06KSkpdOrUiYqKCt5++20WLFjA999/79jG1KlTmThxIr179yY9PZ3p06dTWlrKpEmTGn3/mgpXHOd77rmHf/zjH3z11VcEBwc7ri2Hhobi7+/fuDvYRDT0cQ4ODj6rT1VgYCAREREe3dfKFX/PDz30EAMGDOD5559nzJgxrFq1irfeeou33nqr0fevqXDFcR45ciTPPfccSUlJdOrUiXXr1vHSSy8xefLkRt+/pqKux/mWW27hL3/5C5MmTeKZZ54hPz+fRx55hMmTJzv+3/vAAw+QkZHBiy++yIgRI5g1axZr1qxp3L/nRrvY1kydum3yl9PEiRMNwzCMiRMnGhkZGY7lX3jhBaNNmzaGn5+fER4ebgwZMsRYsGDBWdt99dVXjaSkJMPHx8dIT083VqxY0Uh71DS54jifa3uA8d577zXejjUxrvp7PpP6ALnuOP/73/82OnfubPj6+hodOnQw3nrrrUbao6bJFce5qKjIeOCBB4ykpCTDz8/PaN26tfHEE08YlZWVjbhnTUtdj7NhGMbWrVuNYcOGGf7+/kbLli2NqVOnOvr/nPL5558b7du3N3x8fIxOnToZc+fObaQ9sjMZhgcPbykiIiIeSX2ARERExOMoAImIiIjHUQASERERj6MAJCIiIh5HAUhEREQ8jgKQiIiIeBwFIBEREfE4CkAiIiLicRSARERExOMoAImIiIjHUQASkSZtyJAhPPjgg+4uw+Fi6zl27BjR0dHs3bu3wWv6pZtuuokXX3zR5Z8j0pwpAIkIM2bMIDg4mJqaGkdbSUkJ3t7eDBkyxGnZRYsWYTKZ2LVrVyNX2bgaOng999xzjBo1ilatWjXYNs/nT3/6E8899xyFhYUu/yyR5koBSEQYOnQoJSUlrFmzxtG2ZMkSYmNjWblyJRUVFY72hQsXkpSURJs2bdxRarNUVlbGO++8w+23394on9e5c2fatGnDxx9/3CifJ9IcKQCJCKmpqcTFxbFo0SJH26JFixg1ahQpKSmsWLHCqX3o0KEAzJs3j0GDBhEWFkZERATXXnut05mht956i/j4eGw2m9PnjRo1ismTJwNgs9nIzMwkJSUFf39/unXrxr/+9a/z1lqb5YcMGcL999/PH//4R8LDw4mNjeXPf/6z0zLFxcWMGzeOwMBA4uLiePnllx1nfW677TYWL17MK6+8gslkwmQyOV26stlsF9z2L3377bf4+vrSr18/p/alS5fi7e3tFDD37t2LyWRi3759DBkyhPvuu48HH3yQFi1aEBMTw8yZMyktLWXSpEkEBwfTtm1bvvvuu7M+c+TIkcyaNeuCdYl4MgUgEQHsZ4EWLlzo+HnhwoUMGTKEjIwMR3t5eTkrV650BKDS0lKmTp3KmjVrmD9/Pmazmeuvv94ReG688UaOHTvmtN3jx48zb948xo0bB0BmZiYffvghM2bMYMuWLTz00EPceuutLF68+Jx11nb5Dz74gMDAQFauXMlf//pXpk2bxg8//OCYP3XqVH766Se+/vprfvjhB5YsWUJWVhYAr7zyCv3792fKlCkcPnyYw4cPk5iYWOtt/9KSJUvo1avXWe3r16+nY8eO+Pn5OdrWrVtHixYtSE5OdnxWZGQkq1at4r777uP3v/89N954IwMGDCArK4srr7yS8ePHU1ZW5rTt9PR0Vq1aRWVl5XnrEvFohoiIYRgzZ840AgMDjerqaqOoqMjw8vIy8vLyjH/84x/G4MGDDcMwjPnz5xuAsW/fvnNu4+jRowZgbNq0ydE2atQoY/LkyY6f33zzTSM+Pt6wWq1GRUWFERAQYCxbtsxpO7fffrtx8803G4ZhGBkZGcYDDzxgGIZRq+VPrTNo0CCnZfr06WM8+uijhmEYRlFRkeHt7W3885//dMwvKCgwAgICHJ915uee6de2fS6/PAan3HHHHcaECROc2p566iljyJAh5/ysmpoaIzAw0Bg/fryj7fDhwwZgLF++3Gk7GzZsMABj7969561LxJPpDJCIAPbLRqWlpaxevZolS5bQvn17oqKiyMjIcPQDWrRoEa1btyYpKQmAHTt2cPPNN9O6dWtCQkIcHXxzcnIc2x03bhxffPGF40zEJ598wk033YTZbGbnzp2UlZVxxRVXEBQU5Jg+/PDDc3ayrsvyXbt2dfo5Li6OvLw8AHbv3k11dTXp6emO+aGhoaSmptbqWF1o2+dSXl7udJbnlPXr19O9e3entnXr1jm1nflZFouFiIgIunTp4miLiYkBOOvz/f39Ac46MyQidl7uLkBEmoa2bdvSsmVLFi5cyIkTJ8jIyAAgPj6exMREli1bxsKFC7n88ssd64wcOZLk5GRmzpzp6OvTuXNnqqqqnJYxDIO5c+fSp08flixZwssvvwzY7zQDmDt3LgkJCU71+Pr6nlVjXZb39vZ2+tlkMp3VF+li1XXbkZGRnDhxwqnNarWyefNmevTo4dSelZXFb3/72wt+1pltJpMJ4KzPP378OABRUVG/tjsiHkkBSEQchg4dyqJFizhx4gSPPPKIo33w4MF89913rFq1it///veAfVyb7OxsZs6cyWWXXQbYO/X+kp+fHzfccAOffPIJO3fuJDU1lZ49ewKQlpaGr68vOTk5jsB1IXVd/nxat26Nt7c3q1evdpzNKiwsZPv27QwePBgAHx8frFbrRX/GmXr06HHWHVnZ2dlUVFQQHx/vaFu+fDkHDx4866zQxdi8eTMtW7YkMjKy3tsSuRQpAImIw9ChQ7nnnnuorq52ChgZGRnce++9VFVVOTpAt2jRgoiICN566y3i4uLIycnhscceO+d2x40bx7XXXsuWLVu49dZbHe3BwcE8/PDDPPTQQ9hsNgYNGkRhYSE//fQTISEhTJw40Wk7dV3+fIKDg5k4cSKPPPII4eHhREdH8/TTT2M2mx1nVFq1asXKlSvZu3cvQUFBhIeHYzZfXK+B4cOH8/jjj3PixAlatGgB2C9/Abz66qvcf//97Ny5k/vvvx/A6QzaxVqyZAlXXnllvbcjcqlSHyARcRg6dCjl5eW0bdvW0bcE7AGouLjYcbs8gNlsZtasWaxdu5bOnTvz0EMP8be//e2c27388ssJDw8nOzubW265xWneX/7yF5588kkyMzPp2LEjV111FXPnziUlJeWc26rr8ufz0ksv0b9/f6699lqGDRvGwIEDne7Ievjhh7FYLKSlpREVFeXUr6muunTpQs+ePfn8888dbevXr2f48OHs3r2bLl268MQTT/DMM88QEhLC3//+94v+LICKigq+/PJLpkyZUq/tiFzKTIZhGO4uQkTE3UpLS0lISODFF190yYCFc+fO5ZFHHmHz5s2YzWaGDx9Onz59ePbZZxv8s9544w3mzJnD999/3+DbFrlU6BKYiHikdevWsW3bNtLT0yksLGTatGmAfZBGVxgxYgQ7duzg4MGDJCYmsmHDBsdgkA3N29ubV1991SXbFrlU6AyQiHikdevWcccdd5CdnY2Pjw+9evXipZdecrrF3FVyc3OJi4tjy5YtpKWlufzzRORsCkAiIiLicdQJWkRERDyOApCIiIh4HAUgERER8TgKQCIiIuJxFIBERETE4ygAiYiIiMdRABIRERGPowAkIiIiHkcBSERERDyOApCIiIh4nP8Pys4opNmtYkoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "# loading Meep data\n",
    "meep_data = np.loadtxt(\"misc/MeepMMI.txt\")\n",
    "freqs = meep_data.T[0]\n",
    "flux = meep_data.T[1]\n",
    "\n",
    "ax.plot(td.C_0 / freqs, flux, label=\"Meep\")\n",
    "\n",
    "ax.plot(\n",
    "    td.C_0 / sim_data[\"Flux monitor\"].flux.f,\n",
    "    sim_data[\"Flux monitor\"].flux,\n",
    "    alpha=0.4,\n",
    "    label=\"Tidy3D\",\n",
    ")\n",
    "ax.legend()\n",
    "ax.set_ylim(0.425, 0.525)\n",
    "\n",
    "ax.set_ylabel(\"Transmittance\")\n",
    "ax.set_xlabel(r\"Wavelength ($\\mu$m)\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final Remarks\n",
    "\n",
    "As we can see, the results deviate only due to numerical errors, likely caused by differences in meshing and subpixel averaging strategies.\n",
    "\n",
    "Additionally, one could argue that the speed comparison may not be entirely fair, as the Tidy3D simulation was run on more powerful hardware. However, a direct comparison is not possible since Meep does not support GPU acceleration.\n",
    "\n",
    "From a user’s perspective, both simulations were built and submitted on the same laptop. Therefore, Tidy3D offers the advantage of cloud computing, meaning that users do not need powerful hardware on hand to benefit from the speed advantage, nor do they need to dedicate a workstation for hours to a single simulation."
   ]
  }
 ],
 "metadata": {
  "description": "This notebook presents a benchmark of an MMI interference power divider compared to Meep.",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "keywords": "MMI, Meep, mode solver, modes, Tidy3D, FDTD",
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.13"
  },
  "title": "Benchmark of an MMI Interference Power Divider Compared to Meep | Flexcompute"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
