{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# Inverse design of a planewave absorber with JAX + optax\n", "\n", "Dolphindes has two ways of doing simple inverse design via its simple FDFD Maxwell solver. The first is by working directly with the polarization objectives, for which we have written out custom derivative rules. This is not always convenient: writing objectives in terms of polarization is necessary for limits but not always natural for design. Thus, the second method is with a custom jax wrapper that allows you to differentiate through the field solve.\n", "\n", "The key ingredient is `dolphindes.maxwell.jax_fdfd.build_jax_field_solver`. You hand it a\n", "solver and it hands back a `jax.grad`-differentiable function `field(source, chi) -> Ez`,\n", "where `Ez` is flattened. The heavy linear solve still runs on the CPU via scipy; JAX just\n", "wraps it so `jax.grad` works. JAX is an optional dependency (`pip install dolphindes[jax]`)\n", "and the solver classes themselves never import it.\n", "\n", "**Build it once.** `build_jax_field_solver` copies the Maxwell operator into a native-JAX\n", "sparse array, so call it outside your optimization loop and reuse the returned function.\n", "\n", "**Physics.** The time-averaged absorbed power is\n", "\n", "$$ f = \\frac{\\omega}{2}\\int \\mathrm{Im}\\,\\chi(\\mathbf r)\\,|E_z(\\mathbf r)|^2\\, d\\mathbf r . $$\n", "\n", "We parametrize the structure by a density $\\rho(\\mathbf r)\\in[0,1]$ over a design region,\n", "with local susceptibility $\\chi(\\mathbf r) = \\rho(\\mathbf r)\\,\\chi_\\mathrm{mat}$, and\n", "maximize $f$ over $\\rho$. To keep $\\rho$ in $[0,1]$ we optimize an unconstrained latent\n", "field through a sigmoid.\n", "\n", "> **64-bit JAX is required.** " ] }, { "cell_type": "code", "execution_count": null, "id": "1", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "import jax\n", "jax.config.update(\"jax_enable_x64\", True)\n", "import jax.numpy as jnp\n", "import optax\n", "\n", "from dolphindes import geometry\n", "from dolphindes.maxwell import TM_FDFD\n", "from dolphindes.maxwell.jax_fdfd import build_jax_field_solver" ] }, { "cell_type": "markdown", "id": "2", "metadata": {}, "source": [ "## 1. Geometry, source, and solver\n", "\n", "A unit-amplitude plane wave travelling in $+x$ is launched by a line current just inside\n", "the left PML. The square design region sits in the middle, padded from the PML by vacuum." ] }, { "cell_type": "code", "execution_count": null, "id": "3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "grid = 60 x 60, design pixels = 400\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAG4CAYAAADohIisAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIiRJREFUeJzt3Qt0FOX9//FvMBC5JCl4SUACIkRUEFqlRqI1qBALatHUVo0XrJ4WBKwRLRpoFawSDRrRgli8xvZQbI+i1iomp0LUE/EAiiJefliixkpMwZhEiAmX+Z/v8z+7Z3ezgdzgm919v84ZyM5OdueZncxnnmeenSfO8zxPAAAw1M3yzQEAUIQRAMAcYQQAMEcYAQDMEUYAAHOEEQDAHGEEADBHGAEAzBFGAABzhFGUmTdvnsTFxR201z/22GPlmmuukWinZdSyIrI89dRTbv//7LPP/PP4LCNDvPUKILKsXLlSkpKSJNr94Q9/kBtvvNF6NdAJ+CwjA2GENvnRj37UZbaY3lbx+++/l549e3b6aw8dOlS60tn+r371K1feSLVr1y7p1auXyXt3pc8SLaOZLoL961//kh/+8IeSkJAgQ4YMkfvuuy/scnoQe/jhh92yeuDu27evXHLJJbJ169ag5d5991254IIL5Oijj3avOWDAADn//PPlyy+/3G8z3ebNmyU7O9sdbI466iiZMWOGWzdtLlmzZo1/uXHjxsnIkSNl3bp18pOf/MQtf9xxx8k999wj+/btO2B59fVmzpwpjzzyiJx44oluHYuLi91zW7ZskdzcXP+66/NLlixp9hqtXddwTTsafPn5+W5b9+jRQ4455hj3+99++23Qcvp7uh1XrVolp5xyitvmJ5xwgjzxxBNyqP3jH/+QjIwMSU5O9m/va6+9NmiZL774Qq688sqgbXf//fcHfSa6bUK3kdLmMJ2vgRm47fr06SObNm1y2zoxMVHOPfdc91xjY6Pceeed7j0OP/xwOeKII+Tss8+W8vLyNu+vrRXus/TtS3/5y1/cuui2GT16tLz00kvNfr+1+xY6hppRhPr3v/8tkydPlrFjx8qKFStk7969UlhYKF9//XWzZadOneoOFr/97W/l3nvvlW+++cYdEDIzM+W9996TlJQU2blzp0yYMMEdaPUPTedVVVXJ6tWrpb6+vsX12LZtm2RlZUnv3r1l6dKl7g/2b3/7m/tDD0df84orrpCbb75Z7rjjDtfspwd4Db6rr776gOV+/vnn5Y033pDbb79dUlNT3ft9+OGHriyDBg1yB1Gd/+qrr7rybt++3b1Pe9Y1kB4gL7roIrfddX01TN9//3332m+99Zab9EDlo9tVy3jbbbe5bfnYY4/JddddJ8OGDZOzzjpLDgVdp0svvdRNei1RD/6ff/65vPbaa/5l/ve//7lt19TUJH/84x/dQVsPyLfccov85z//caHQHvp6P/vZz9y+p9tgz549bpo4caL7/PLy8uScc85x89auXesCUdejtftrZ9CTED0x0tfW8NS/n4svvlg++eQTF9qqtfsWOoEOIYHIk5GR4Q0YMMBraGjwz6urq/P69eunbTn+eW+99ZZ7fP/99wf9fmVlpdezZ09v9uzZ7vH69evdcs8///x+33fw4MHelClT/I9/97vfeXFxcd7mzZuDljvvvPPc661evdo/Lysry817++23g5Y96aST3PIHor+bnJzsffPNN83ea+DAgV5tbW3Q/JkzZ3qHH364f/m2rKuWUcvqs2rVKrdMYWFh0O8+88wzbv6yZcuCtpG+7+eff+6fp5+TfjZTp049YDn37Nnj7d692z89/vjj7j0C5+m0d+/e/b7Offfd537v22+/bXGZ2267Lexncv3117tt9cknn7jHum1Ct5GqqKhw85988smgbafznnjiiaBln376aTf/0UcfbXF9Wru/tkTXQ39f16ulz1LpMikpKe5vxqeqqsrr1q2bV1BQ0OZ9Cx1HM10E0lqMntHl5OS4s10fbQ658MILg5bVs1xtktBmGN/ZqU56hqfNEr5mFz1j1+aQW2+91TWD6Rlha5SVlbmmt5NOOilo/uWXXx52eX3f0047LWjeqFGj3Bl7a+jZtK5nYNOZ1lb0jFabWgLLOGnSJPe8nnm3Z10D+WoToU2Uv/jFL1xNS9chkDYx6dm0j35Oxx9/fKvKqdc4unfv7p+0RqUC5+mkZ/T78+Mf/9j9/8tf/lL+/ve/y3//+9+w5dLtEfqZaDn1mB1Yi2qrn//850GPX3nlFbcdQpsJ27O/dgZtHtS/GR+tcWlt2fcZtWXfQsfRTBeBampqXHu+/oGGCp2nzXZ6UGmpacPXHKHXFPRgfffdd8ucOXPce/Tv319+/etfy+9//3t38Atnx44drmkvVEvvp9cIQmnzVkNDg7SGrlPo++vB4U9/+pObwtHmlPasa+j7xMfHu+tMgfTAqdtcn++scv7zn/9011YCD9Dz5893JyCBtGlzf7Q5UJs1H3roIdcEqq85YsQImTt3rj+Adb3DdWH3vXZouVpLD96hvS61SVBft1u3ls+BW7u/doYDfUZt2bfQcYRRBNKagR4E9fpLqNB5Rx55pFtW2+kDr2n4BM47+eST3fUnPRjo9RBtt9ezb72IrO3+Lf1Bh7tOFW7dOkPod6h0Wxx22GFy1VVXuc4E4fgCqCPrqr+rByY9oAYGkm4r/X1fLaQz6OcQ6IMPPnD/jxkzps2vpdcVddIg0rP4goICdzFeA0ivN2q59FpaqK+++sq//yhfDTwwJPd3MA73XTfdbm+++aY7kWopkNqyvx5sbdm30HE000UgbRbSZpXnnnvONRX4aEcDPasOpL269ICpTTR6MAudQg98Sg8G2iTywAMPyA9+8AN55513WlwX7RCgB8vQZj0NtUNBz8C1uUV7AmpzX7gy+s6AO7Kuvt5gf/3rX4PmP/vss67Z1Pd8V6UHcS2/dghQur2Urrduj9DP+Omnn3b7gW5b5as96UlKoBdffLHV66CdF3R/Dex5F6o9+2tX2LfQcdSMIpT2fPrpT3/qesBpry3tTacHGg0q7X3kc8YZZ8hvfvMb9z2V9evXu6YbXUbPhvUsVf+4r7/+etcUpD2ntMeYNoXoAUHDTrst63u0RHtFaZdlPdBoLUqbV5YvXy4ff/yxe35/TTKd5cEHH5QzzzzT9XDTsuiBU4P5008/deHsu+7RkXXVbXDeeee5a2p1dXVuu/p60+l3r/TsuavRHofaLV8DZ+DAge6z1G2lTa4aTOqmm25ywaNd+HWbDB482PUy031Bt6Ve51LaFDl+/HhXs9Iagy6n11N0H2ktbRp88sknZdq0aa7Hmh7otZb09ttvu+7Sl112Wav310OltfsWOkEndIKAkRdffNEbNWqU16NHD2/QoEHePffc491xxx1Bvel8tGeT9sDr3bu365U0dOhQ7+qrr3a96NTHH3/sXX755W6+Pq+91k477TTvqaee2m9vOvXBBx9448ePd72LtMfYdddd5xUXF7v1eO+994J6040YMaLZuoXr7RSOvt6MGTPCPqe9p6699lrvmGOO8bp37+4dddRRXmZmpnfXXXe1a13DrZP2iLv11lvdfH2P/v37u15nNTU1zbbR+eef32wdtfw6tZWvh1hbvfTSS97EiRPdNtF95Oijj/YmTZrkvfHGG0HLaa+/3Nxc74gjjnDlGj58uLdw4cJmvfW2bdvmXXLJJW676f5x5ZVX+nthhvam0/0sHN2Gt99+u5eenu7WSd/znHPO8crLy9u0v3ZGb7pw+1K4/bu1+xY6Jk7/6YxQAwLp2a1+h0cvAusXRLuySFpXIFrRTIcO0+Yd7SWlzXvfffeda/LTL3lqL7yudnCPpHUFYglhhA7TaxALFy501ye0x1l6eroUFRV1yRuNRtK6ArGEZjoAgDm6dgMAzBFGAABzhBEAwBxhBAAwRxgBAMx16TDSW5LojQj1Jo2nnnqqu3liNHn99dfdkA/6vRe9D5jeYTmQfh9ZB0XT5/VmpTpSqo5UGun0ljJ6Y1G9fb/esl9vQaS3h4n2suuAfnqPM72btU56o1IdViGayxzus9d9XW/NFO3l1jJpWQOnwLvqR2u5oy6MnnnmGbfD6u3u9UaFem8ovaeYjggZLfQGm3pD0sWLF4d9Xkee1O/A6PM6fIDuyHqPtP2NvBoJdKgKvQuy3kW6tLTUfd9Hh6fW7RHNZdf7w+kQ63rPNZ10bCa9o7bvABSNZQ6kZVq2bJkL5EDRXG4dskPvq+ebdCj2WCh3u3hdlN4Xbdq0aUHzTjjhBDcyZTTSj2LlypX+x/v27fNSU1Pd/eZ8vv/+e3dPsEceecSLJtXV1a78ZWVlMVf2vn37eo899ljUl7m+vt7dj660tNTdn+/GG29086O53HqfyNGjR4d9LprL3V5dsmbU1NQkGzZscGfLgfRxeXm5xIKKigo3Tk7gNvANAxBt26C2ttb9369fv5gpu95lXYeu0NqgNtdFe5m1Jqx3Btc7fweK9nJv2bLFNcPp5Qa9K/nWrVtjotxRczsgHbBL/1hDR3vUxwdr0LauxlfOcNugtUN0RwKtFM6aNcvdpl+HBI/2smszjYaPjuvTp08fWblypRv223cAisYya+jqeEmhI9VG+2edkZHhhufQYTh0UMe77rpLMjMzXbNsNJc7qsKopdEi9cAVbgTJaBbt22DmzJluXCAdqyYWyj58+HDZuHGjG1tIB+abMmWKu4YWrWWurKx09/0rKSnxjxYbTrSVW+k1bh8dh0lPQoYOHSrFxcVy+umnR22526tLNtPp0MM63G9oLai6urrZmUS08vW6ieZtcMMNN7iRQlevXu0u7sdC2fXO4MOGDXOjhGrPMu3AogO4RWuZtbldy6C9YePj492k4fvQQw+5n31li7Zyh6ODBGooadNdtH7eURdG+gerO6/2tAqkj7WaGwu0jVl32MBtoNfS9A850reBnv1pjUhHCdWRMrWssVL2cNuisbExasuso8xq06TWBn2TBvEVV1zhftahPKKx3OHo5/zRRx9J//79o/bz7hCvi1qxYoUbVfHxxx/3PvzwQy8vL8+N+vjZZ5950UJ7GL377rtu0o+iqKjI/awjbyrtaaO9a5577jlv06ZNbiRWHV20rq7Oi2Q6OqqWa82aNW70UN+0a9cu/zLRWPb8/Hzv9ddfdyOHvv/++96cOXO8bt26eSUlJVFb5nACe9NFc7lvvvlmt49v3brVW7t2rXfBBRd4iYmJ/mNYtJa7vbpsGKklS5a4YYB1eOJTTjnF3/U3WqxevdqFUOjkG/ZYu39q91DtApqQkOCdddZZbqeNdOHKHDp0dTSWXYeu9u3POnT1ueee6w+iaC1za8IoWst96aWXunDRk+oBAwZ4OTk53ubNm6O+3O3FeEYAAHNd8poRACC2EEYAAHOEEQDAHGEEADBHGAEAzBFGAABz3SLhW8s6AJX+H0soN593LGA/j639fH+6/PeM6urqJDk52Q0zoKNjxgrKzecdC9jPY2s/j+iaEQAg+hFGAABzXW48o3379slXX30liYmJblwPrcYr3/+xgnLzeccC9vPo3s89z5P6+no32m23bgeo+3gH8Sanxx57rLsBoN7kVO9W3BqVlZUt3kiTiW3APsA+wD4gEbcN9Lh+IAelZvTMM89IXl6ePPzww3LGGWfIn//8Zzfq4YcffiiDBg3a7+9qjUidKZMkXrofjNUDABwCe2S3vCkv+4/rh7w3nY79fsopp8jSpUv980488US56KKL3OiWoV07A7tta7U9LS1NxslkiY8jjAAgUu3xdssaeaFVvaE7vQODjlaoQw1nZ2cHzdfH5eXlzZbXcNKu275JgwgAEFs6PYy2b98ue/fubTaOuz4OHe9d5efnu9T0TZWVlZ29SgCAWO1Npz3hAmlrYOg8lZCQ4CYAQOzq9JrRkUceKYcddlizWlB1dXWz2hIAAAcljHr06CGnnnqqlJaWBs3Xx5mZmWx1AMChaaabNWuWXHXVVTJmzBgZO3asLFu2TL744guZNm3awXg7AECEOyhhdOmll8qOHTvkzjvvlG3btsnIkSPl5ZdflsGDBx+MtwMARLgud9du3118+Z4RAEQ20+8ZAQDQVoQRAMAcYQQAMEcYAQDMEUYAAHOEEQDAHGEEADBHGAEAzBFGAABzhBEAwBxhBAAwRxgBAMwRRgAAc4QRAMAcYQQAMEcYAQDMEUYAAHOEEQDAHGEEADBHGAEAzBFGAABzhBEAwBxhBAAwRxgBAMwRRgAAc4QRAMAcYQQAMEcYAQDMEUYAAHOEEQDAHGEEADBHGAEAzBFGAABzhBEAwBxhBAAwRxgBAMwRRgAAc4QRAMAcYQQAMEcYAQDMEUYAAHOEEQDAHGEEADBHGAEAzBFGAABzhBEAwBxhBAAwRxgBAMwRRgAAc4QRAMAcYQQAMBff1l94/fXXZeHChbJhwwbZtm2brFy5Ui666CL/857nyfz582XZsmVSU1MjGRkZsmTJEhkxYkSb3mfl/22SpESyEgAiVV39Pul7fOuWbfPRfufOnTJ69GhZvHhx2OcLCwulqKjIPb9u3TpJTU2VCRMmSH19fVvfCgAQI9pcM5o4caKbwtFa0aJFi2Tu3LmSk5Pj5hUXF0tKSoosX75cpk6d2ux3Ghsb3eRTV1fX1lUCAES4Tm0Hq6iokKqqKsnOzvbPS0hIkKysLCkvLw/7OwUFBZKcnOyf0tLSOnOVAACxFkYaREprQoH0se+5UPn5+VJbW+ufKisrO3OVAADR2EzXGnFxcc2a70LnBdacdAIAxK5ODSPtrKC0FtS/f3///Orq6ma1pQO5+PiTJT6ue2euHgDgENrj7RaRrYe+mW7IkCEukEpLS/3zmpqapKysTDIzMzvzrQAAsVwz+u677+TTTz8N6rSwceNG6devnwwaNEjy8vJkwYIFkp6e7ib9uVevXpKbm9vZ6w4AiNUwWr9+vZx99tn+x7NmzXL/T5kyRZ566imZPXu2NDQ0yPTp0/1fei0pKZHExMTOXXMAQNSI87R3QRei3zPSLt7jZDLXjAAgwq8ZrZEXXE/ppKSk/S7L/XYAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGAuvi0LFxQUyHPPPScff/yx9OzZUzIzM+Xee++V4cOH+5fxPE/mz58vy5Ytk5qaGsnIyJAlS5bIiBEjDsb6Iwp8+sDp1quAQ2jYTWvZ3uhYzaisrExmzJgha9euldLSUtmzZ49kZ2fLzp07/csUFhZKUVGRLF68WNatWyepqakyYcIEqa+vb8tbAQBiSJtqRqtWrQp6/OSTT8rRRx8tGzZskLPOOsvVihYtWiRz586VnJwct0xxcbGkpKTI8uXLZerUqc1es7Gx0U0+dXV17S8NACD2rhnV1ta6//v16+f+r6iokKqqKldb8klISJCsrCwpLy9vsekvOTnZP6WlpXVklQAAsRRGWguaNWuWnHnmmTJy5Eg3T4NIaU0okD72PRcqPz/fhZpvqqysbO8qAQBioZku0MyZM+X999+XN998s9lzcXFxzYIrdF5gzUknAEDsalfN6IYbbpAXX3xRVq9eLQMHDvTP184KKrQWVF1d3ay2BABAu8JIazhaI9Lu3a+99poMGTIk6Hl9rIGkPe18mpqaXC887QYOAECHm+m0W7f2invhhRckMTHRXwPSjgf6vSNtisvLy5MFCxZIenq6m/TnXr16SW5ublveCgAQQ9oURkuXLnX/jxs3rlkX72uuucb9PHv2bGloaJDp06f7v/RaUlLiwgsAgA6HkTbTHYjWjubNm+cmAABag3vTAQDMEUYAAHOEEQDAHGEEADBHGAEAzBFGAABzhBEAwBxhBAAwRxgBAMwRRgAAc4QRAMAcYQQAMEcYAQDMEUYAAHOEEQDAHGEEADBHGAEAzBFGAABzhBEAwBxhBAAwRxgBAMwRRgAAc4QRAMAcYQQAMEcYAQDMEUYAAHOEEQDAHGEEADBHGAEAzBFGAABzhBEAwBxhBAAwRxgBAMwRRgAAc4QRAMAcYQQAMEcYAQDMEUYAAHOEEQDAHGEEADBHGAEAzBFGAABzhBEAwBxhBAAwRxgBAMwRRgAAc4QRAMAcYQQAMEcYAQDMEUYAAHOEEQAgssJo6dKlMmrUKElKSnLT2LFj5ZVXXvE/73mezJs3TwYMGCA9e/aUcePGyebNmw/GegMAYjWMBg4cKPfcc4+sX7/eTeecc45MnjzZHziFhYVSVFQkixcvlnXr1klqaqpMmDBB6uvrD9b6AwBiLYwuvPBCmTRpkhx//PFuuvvuu6VPnz6ydu1aVytatGiRzJ07V3JycmTkyJFSXFwsu3btkuXLl7f4mo2NjVJXVxc0AQBiS7uvGe3du1dWrFghO3fudM11FRUVUlVVJdnZ2f5lEhISJCsrS8rLy1t8nYKCAklOTvZPaWlp7V0lAECshNGmTZtcbUiDZtq0abJy5Uo56aSTXBCplJSUoOX1se+5cPLz86W2ttY/VVZWtqccAIAIFt/WXxg+fLhs3LhRvv32W3n22WdlypQpUlZW5n8+Li4uaHltvgudF0hDTScAQOxqc82oR48eMmzYMBkzZoxrYhs9erQ8+OCDrrOCCq0FVVdXN6stAQDQqd8z0pqPdkIYMmSIC6TS0lL/c01NTa7WlJmZ2dG3AQBEsTY1082ZM0cmTpzoOhlod23twLBmzRpZtWqVa4rLy8uTBQsWSHp6upv05169eklubu7BKwEAILbC6Ouvv5arrrpKtm3b5nq+6RdgNYj0u0Rq9uzZ0tDQINOnT5eamhrJyMiQkpISSUxMPFjrDwCIAnGetrN1Ifo9Iw26cTJZ4uO6W68ODoFPHzid7RxDht201noVcIjs8XbLGnnB9ZTWu/bsD/emAwCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAufiO/HJBQYHMmTNHbrzxRlm0aJGb53mezJ8/X5YtWyY1NTWSkZEhS5YskREjRnTWOiPKDLtprfUqAIjUmtG6detc4IwaNSpofmFhoRQVFcnixYvdMqmpqTJhwgSpr6/vjPUFAEShdoXRd999J1dccYU8+uij0rdvX/98rRVpDWnu3LmSk5MjI0eOlOLiYtm1a5csX7487Gs1NjZKXV1d0AQAiC3tCqMZM2bI+eefL+PHjw+aX1FRIVVVVZKdne2fl5CQIFlZWVJeXt5iU19ycrJ/SktLa88qAQBiKYxWrFgh77zzjguRUBpEKiUlJWi+PvY9Fyo/P19qa2v9U2VlZVtXCQAQSx0YNCi0s0JJSYkcfvjhLS4XFxcX9Fib70LnBdacdAIAxK421Yw2bNgg1dXVcuqpp0p8fLybysrK5KGHHnI/+2pEobUg/Z3Q2hIAAO0Ko3PPPVc2bdokGzdu9E9jxoxxnRn05+OOO871nistLfX/TlNTkwuszMzMtrwVACCGtKmZLjEx0fWQC9S7d2854ogj/PPz8vJkwYIFkp6e7ib9uVevXpKbm9u5aw4AiBod+tJrOLNnz5aGhgaZPn26/0uveo1JgwwAgHDiPO1d0IXo94y0i/c4mSzxcd2tVwcA0E57vN2yRl5wPaWTkpL2uyz3pgMAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAQGSF0bx58yQuLi5oSk1N9T/veZ5bZsCAAdKzZ08ZN26cbN68+WCsNwAglmtGI0aMkG3btvmnTZs2+Z8rLCyUoqIiWbx4saxbt84F1YQJE6S+vr6z1xsAEEXi2/wL8fFBtaHAWtGiRYtk7ty5kpOT4+YVFxdLSkqKLF++XKZOnRr29RobG93kU1dX19ZVAgDEWs1oy5YtrhluyJAhctlll8nWrVvd/IqKCqmqqpLs7Gz/sgkJCZKVlSXl5eUtvl5BQYEkJyf7p7S0tPaWBQAQC2GUkZEhTz/9tLz66qvy6KOPuvDJzMyUHTt2uJ+V1oQC6WPfc+Hk5+dLbW2tf6qsrGxvWQAAsdBMN3HiRP/PJ598sowdO1aGDh3qmuNOP/10N187NYQ234XOC6S1J50AALGrQ127e/fu7UJJm+5815FCa0HV1dXNaksAAHRaGGnHg48++kj69+/vriFpIJWWlvqfb2pqkrKyMteUBwBApzTT3XLLLXLhhRfKoEGDXI3nrrvucr3fpkyZ4pri8vLyZMGCBZKenu4m/blXr16Sm5vblrcBAMSYNoXRl19+KZdffrls375djjrqKHedaO3atTJ48GD3/OzZs6WhoUGmT58uNTU1rsNDSUmJJCYmHqz1BwBEgThPexh0IVrT0i7e42SyxMd1t14dAEA77fF2yxp5wfWUTkpK2u+y3JsOAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmCOMAADmCCMAgDnCCABgjjACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAOcIIAGCOMAIAmIuXLsbzPPf/Htkt8v9/BABEIHccDziuR1QY1dfXu//flJetVwUA0EnH9eTk5P0uE+e1JrIOoX379slXX30liYmJEhcXJ3V1dZKWliaVlZWSlJQksYJy83nHAvbz6N7PPc9zQTRgwADp1q1bZNWMdIUHDhzYbL5+YNH8obWEcscWPu/YEgufd/IBakQ+dGAAAJgjjAAA5rp8GCUkJMgdd9zh/o8llJvPOxawn8fWfh5RHRgAALGny9eMAADRjzACAJgjjAAA5ggjAIA5wggAYI4wAgCYI4wAAOYIIwCAWPt/kz0mTkMF2lcAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "wavelength = 1.0\n", "omega = 2 * np.pi / wavelength\n", "\n", "px_per_wavelength = 20\n", "dl = 1.0 / px_per_wavelength\n", "\n", "des_size = 1.0 # square design region side length (in wavelengths)\n", "pml_sep = 0.5 # vacuum padding between design region and PML\n", "pml_thick = 0.5\n", "\n", "M = int(des_size / dl)\n", "Nsep = int(pml_sep / dl)\n", "Npml = int(pml_thick / dl)\n", "Nx = Ny = M + 2 * (Nsep + Npml)\n", "\n", "# design region: where material is allowed to be placed\n", "des_mask = np.zeros((Nx, Ny), dtype=bool)\n", "des_mask[Npml + Nsep:-(Npml + Nsep), Npml + Nsep:-(Npml + Nsep)] = True\n", "des_idx = np.where(des_mask.ravel())[0]\n", "Ndes = des_idx.size\n", "\n", "# unit-amplitude planewave in +x: line current just inside the left PML\n", "ji = np.zeros((Nx, Ny), dtype=complex)\n", "ji[Npml, :] = 2.0 / dl\n", "\n", "geo = geometry.CartesianFDFDGeometry(Nx=Nx, Ny=Ny, Npmlx=Npml, Npmly=Npml, dx=dl, dy=dl)\n", "solver = TM_FDFD(omega, geo)\n", "\n", "print(f\"grid = {Nx} x {Ny}, design pixels = {Ndes}\")\n", "plt.matshow(des_mask + np.real(ji) * dl)\n", "plt.title(\"design region + source line\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "4", "metadata": {}, "source": [ "## 2. Incident field\n", "\n", "The incident (vacuum) field is just the ordinary numpy solve. A correctly normalized\n", "unit-amplitude plane wave has $|E_z|\\approx 1$ in the interior." ] }, { "cell_type": "code", "execution_count": null, "id": "5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "incident |E| near center: 1.0126723055251141\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAAGzCAYAAABZzq+8AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPM9JREFUeJzt3Ql4FGW2+P8TAwkQSQAjJFwCoiIgQYZFNgcR2RVccARFER1UGESI6AUZdFiuguAMArI4MCgwqOFxlFF/RiQugBpBQJBF9DIjSFAwgpCELYFQ/+e8c7v/6exdSSddXd/P85RJV1d3v1XBfk+ddwuzLMsSAADgKhdVdQEAAEDlIwAAAMCFCAAAAHAhAgAAAFyIAAAAABciAAAAwIUIAAAAcCECAAAAXIgAAAAAFyIAAADAhQgAAABwIQIA2LZ8+XIJCwuTAwcOVOn7T5061RwXaGlpaeazTpw4UebXrF69Wlq1aiU1a9Y0ZdyxY4ft6+bP68pzTTyf49mqVasm8fHxctddd8m+ffukoqSkpPh8TsFt2bJlFfZZAAqrVsQ+oExuvvlm+eKLL0zl4MT3txMATJs2Te6//36pU6dOqcf/8ssvMmzYMOnXr58sWrRIIiMj5aqrrpL/+q//CqrzKs4rr7wiLVq0kLNnz8rnn38uzz77rHzyySfy7bffSt26dcv9/l999ZX5+fbbb0v9+vULPX/11VeX+zMAFI8AALZdeumlZnPq+wfa//7v/8q5c+fk3nvvle7du3v316pVyxHnlZiYKB06dDC/33DDDZKXlydTpkyRf/7zn/LAAw9USAAQHR0tAwcOrJQMDgBfNAHAtoIpaU/aec+ePXL33XdLTEyMNGjQQH7/+99LZmZmodfrnaQep8fo3XHjxo3lvvvuk5ycnCLfX7333nvym9/8xhzftGlT+fOf/1xs+TRdPXToUHN3qce3bNlSFi5c6HNMWcusx/33f/+3+V0/15OmXr9+fZGfrVmC3/72t+b3IUOGmGO1Ei3uvMpS1uL4c03KwxMM/Pzzz4Wes1P+bdu2SZs2baj8gSpCBgAV7o477jCV3ogRI2TXrl0yadIks//ll1/2HvP111+bCjI2NlamT58uzZo1k8OHD8s777wjubm5phIp6KOPPpJbb71VunTpIsnJyeaOdPbs2UVWSN9884107drVBBV/+ctfJC4uTj744AMZO3asHD161NzJ+lPmBx98UH799Vd58cUX5a233vKm74tLUz/99NPSsWNHeeSRR2TGjBnSo0cPc7dbFH/LavealNf+/fvNT23GKG/5jx07JgcPHpQBAwbI+fPnCz0fHh5OYAAEmgXY9Morr1j6T2j//v3m8ZQpU8zj2bNn+xw3evRoq0aNGtaFCxe8+2688UarTp06VkZGRpnfv1OnTlbDhg2tM2fOeI/Jysqy6tWrZ47Lr2/fvlajRo2szMxMn/1jxowxZfn111/9LvPzzz/vU57SfPLJJ+b4N954o8TzKmtZC77O32tSVp7P2bRpk3Xu3DkrOzvbWrt2rRUXF2ddf/31Zl9+ZS1/fuvWrTOfUdy2e/duW2UHUHY0AaDC3XLLLT6Pr7nmGtORLCMjwzw+ffq0bNiwQQYPHlzmtvBTp07Jli1bZNCgQVKjRg3v/tq1a5s25Pz0s/TO+Pbbbzft7XqH6dluuukm8/ymTZv8KnOg2CmrnWtiR+fOnaV69erm/bQjo3b80w57OiqgvOXX9L/SbIqeQ8GNDoBA4NEEgAp3ySWX+Dz2pPPPnDljfh4/ftykqhs1alTm99TXXLhwwaSXCyq4T9PLWgFpul63omhq2p8yB4qdstq5JnasXLnStOVnZ2eb4Yx//etfTT+J999/v9zl1w6AGrRo4KXp/uJo84YGGPmDnk8//dTbvwKAfQQAqHT16tUzX/qHDh0q82v07lM7zh05cqTQcwX36bH6/joET9vgi6Kd5YJBecrqzzWxQyt/T8c/7cOgQdvf/vY3+cc//iG/+93vylV+DQB0lEFJlb/SjIPHuHHjzL8Z7e8AoPwIAFDpdFIcHRb3xhtvmLHl2hGwNFFRUaZTnaaMn3/+eW/KW+9O3333XZ9jNRWtFdb27dtNKj8iIqJCyh2IrEB5yurPNakI2rnwzTfflD/96U+m2eGiiy6yVX4dXfH999+bkRZlpZ0ydaSBDkEsLWgAUDb0AUCVmDNnjhkj36lTJ1m6dKmZYEZ7setQMq3AivI///M/5s62d+/epiLQyqhnz56mIixo3rx5ppd5t27dzLA7Ha6nleILL7wgN954o60yt27d2vveOpHP1q1biy2rP8pTVn+uSXnp3b5WxHv37pXXXnvNdvn17t+yLFNG7R9QcEtPT/c5XkeJbN682QQ6FRXMASAAQBXR8d9ffvmltG/f3lQq2sls4sSJ5i67uC95TyWXlZVlhuyNHz/eDN8r6k5SO5F50sxPPfWU9OnTxwzx0/S1VpB26Dh+LatWbtoGfe2113o7s5VHecrqzzWpCI8++qgZ7qeVsjYJ2Cm/ZwbA+fPnm3R+wW3dunU+geLatWvN8ND8HR0BlF+YDgWogPcBgAq1ePFi0+fg448/NhM0AahYBAAAgpKut6DDCPMPO9SmIs28ACg/AgAAAFyIToAAALgQAQAAAC5EAAAAgAsRAAAA4IeNGzea9TYaNmxoZuPUobil0fVPdNizDme9/PLL5aWXXip0jM7jocNqdTi0/lyzZo04cibARYsWmdnJdInXVq1aydy5c81EIaXRuc1/+uknswCJXlgAgLPo6HKdJEsrSJ0xMlB0lIguH15eOveIP/NM6JoUOpfJAw88YObdKMtS2ro41kMPPSSrVq2Szz//XEaPHm0WQ/O8XicX07k8dHIvXVxLK39dMO2zzz4zE6YFhBUAycnJVvXq1a2lS5da33zzjTVu3DgrKirK+uGHH0p9bXp6eonLhLJxDfg3wL8B/g0449+Afp8Hii6BHVdB5YyLi/NZUtsf+vo1a9aUeMyECROsFi1a+OwbOXKk1blzZ+/jwYMHW/369Su01PZdd91lBUpAMgA6e5fOBPbggw+ax3r3/8EHH5iJPWbOnFnia/XOX6UfPCjR0dGBKB4AIIB0ZsqExo293+eBoHf+uuRVeliYlKemyBKRhCNHzKqV+escTcN71v8oL7271xky8+vbt68sW7bMTImuy27rMY899lihY7T+DJRqgfij6PSoTz75pM9+Pfm0tLRCx+fk5JjNwzO3uv4hCAAAwLkqoxlXq+zo8nyO9Z/JcBMSEnx2T5kyRaZOnSoVQdfraNCggc8+faxLaWvgER8fX+wxFbGyZ6UFAHoyOkd4WU9EMwLTpk2r6GIAANxA+xiUNwDIyzOLUBXMAAQyGPLMwp9/f1HHBDKICljvjLKeiC6uosuDeraCK4EBAFBiAFDeTf7/rLNnq8gAIC4urtANcEZGhpnm+pJLLinxmII300EdAOja7rped1lPRC9ywQsPAEBlBgCBpKtcpqam+uzTVS87dOhg2v9LOqZr164BK1eFn7kOp9CxjgVPRB8H8kQAAKgMJ0+elB07dpjNM8xPfz948KA3s33fffd5jx81apT88MMPZrnuvXv3yssvv2w6AD7xxBPeY8aNG2cq/FmzZsm3335rfn744YeSlJQUsPMIyCgAPclhw4aZ6EajmiVLlpgLoxcBAICg6gPgp61bt0qPHj186jw1fPhwWb58uZn/xhMMqKZNm0pKSorp5b9w4UIzP8L8+fN95hDQG+Tk5GR56qmn5Omnn5YrrrhCVq9eHbg5AAK5GqBOBDR79mxzIRITE+WFF16Q66+/vkzDR3Tt78wTJ2gOAAAHMt/jdeqYfl2Batb11hVRUeUaBZBlWRJz6lRAyxqsgm45YAIAAHA2AgBnCNhUwAAAhGITQKggAAAAOBcBgG2sBggAgAuRAQAAOBcZANsIAAAAzqXt/+WZzOfCBXErmgAAAHAhMgAAAOeqpOl8QxEBAADAuQgAbCMAAAA4FwGAbeRNAABwITIAAADnIgNgGwEAAMC5CABsowkAAAAXIgMAAHAuMgAhGACcP/+fDQDgLJX53U0AYBtNAAAAuFDwZgAAACgNGQDbCAAAAO5dDMiyxK1oAgAAwIXIAAAA3NsEYLk3AxC8AUC1av/ZAADOUpnf3QQAtlHDAgCciwDANvoAAADgQmQAAADORQbANgIAAIBzEQDYRhMAAAAuFLQZgH3/CpOLLw6r6mIAAPx08mQlfneTAQi9AAAAgFIRANhGEwAAAC5EBgAA4FxkAGwjAAAAuHcxoAsXxK1oAgAAwIWCNgOQkiJSo0ZVlwIA4K+zZx3UBHCRe++DgzYAAACgVAQAtrk39AEAhE4AUJ7NhkWLFknTpk2lRo0a0r59e/n000+LPfb++++XsLCwQlurVq28xyxfvrzIY84GMJ1CAAAAgB9Wr14tSUlJMnnyZNm+fbt069ZN+vfvLwcPHizy+Hnz5snhw4e9W3p6utSrV0/uvPNOn+Oio6N9jtNNA4xAIQAAADhXFWQA5syZIyNGjJAHH3xQWrZsKXPnzpWEhARZvHhxkcfHxMRIXFycd9u6dascP35cHnjgAZ/j9I4//3G6BRIBAADAuSooAMjKyvLZcnJyivy43Nxc2bZtm/Tp08dnvz5OS0srU5GXLVsmvXr1kiZNmvjsP3nypNnXqFEjGTBggMkuuHYUQLWgLR0AoDjnzzvv2iQkJPg8njJlikydOrXQcUePHpW8vDxp0KCBz359fOTIkVI/R9P677//vrz22ms++1u0aGH6AbRu3doEINpscN1118nXX38tzZo1k0CgigUAiNtHAaSnp5s2eI/IyMgSX6bp+vwsyyq0ryhaydepU0duu+02n/2dO3c2m4dW/u3atZMXX3xR5s+fL4FAAAAAELcHANHR0T4BQHFiY2MlPDy80N1+RkZGoaxAQRokvPzyyzJs2DCJiIgopVgXybXXXiv79u2TQKEPAAAAZaQVtw77S01N9dmvj7t27Vriazds2CD/+te/TAfC0miwsGPHDomPj5dAIQMAAHCuKpgIaPz48eYuvkOHDtKlSxdZsmSJGQI4atQo8/ykSZPkxx9/lJUrVxbq/NepUydJTEws9J7Tpk0zTQDa3q99ADTtrwHAwoULJVAIAAAA7l0MKKz0dvuChgwZIseOHZPp06ebTn1aoaekpHh79eu+gnMCZGZmyptvvmk69xXlxIkT8vDDD5umBR022LZtW9m4caN07NhRAiXM0jxDENHIR0++du1MCQsrvT0GABBcLCtLsrNjTKVXlnb18tQVmb17S3T16vbf59w5iUlNDWhZgxUZAACAc7EWgG0EAAAA5yIAsI0AAADgXAQAtjEMEAAAFyIDAABwLjIAthEAAACciwCg8gIAHZf4/PPPm9WQdKzjmjVrfOY01lGFOqGBToygyx3qpAc6kUGrVq38+pzsbH9LBgAAAtYH4NSpU9KmTRtZsGBBkc/Pnj3brJWsz2/ZssWsZ9y7d2/JpkYHAATpcsBu5HcGoH///mYrit79z507VyZPniyDBg0y+1asWGEWSNClD0eOHFnoNbrmcv51l3VyBwAAyoQmANsqNPTZv3+/mcawT58+Pksqdu/eXdLS0op8zcyZM81sTp6t4JrMAAAgyAMAz/KIBZdE1McFl0700EUTdApGz6ZrMgMAUCY0AQTXKICwAosraNNAwX35MwS6AQDghMWAQkWFBgDa4U/p3X7+NYwzMjIKZQVKExXl6r8LADiWLjF36lRVlwKV2gTQtGlTEwSkpqZ69+Xm5sqGDRuka9euFflRAADQBFCZGYCTJ0/Kv/71L5+Ofzt27JB69epJ48aNJSkpSWbMmCHNmjUzm/5eq1YtGTp0aHnKCQBAYYwCqLwAYOvWrdKjRw/v4/Hjx5ufw4cPl+XLl8uECRPkzJkzMnr0aO9EQOvWrZPatWvbLyUAAEUhALAtzNIeekFE5wHQ4YBRUZkSFhZd1cUBAPjJsrLk1KkYM7IrOjo6oHVF5r33SnREhP33yc2VmFWrAlrWYMVaAAAA5yIDEHoBwOWXi4SHV3UpAAD+yssT2bWrkq4bAYBt7p0EGQAAFwvaDAAAAKUiA2AbAQAAwLkIAGyjCQAAABciAwAAcC4yAKEXAHTsKFKOoZ0AgCqSm1uJowBYDMg2mgAAAHChoM0AAABQKpoAbCMAAAA4FwGAbQQAAADnIgCwjT4AAAC4UNBmAHr3FqlVq6pLAQDw1+nTIsuWVdJ1IwMQegEAAAClIgCwjSYAAABciAwAAMC5yADYRgYAAOD8AKA8mw2LFi2Spk2bSo0aNaR9+/by6aefFnvs+vXrJSwsrND27bff+hz35ptvytVXXy2RkZHm55o1aySQCAAAAPDD6tWrJSkpSSZPnizbt2+Xbt26Sf/+/eXgwYMlvu67776Tw4cPe7dmzZp5n/viiy9kyJAhMmzYMPn666/Nz8GDB8vmzZslUMIsy7IkiGRlZUlMTIxk/vqrREdHV3VxAAB2vsfr1ZPMzMyAfY9764qJEyU6MtL+++TkSMysWX6VtVOnTtKuXTtZvHixd1/Lli3ltttuk5kzZxaZAejRo4ccP35c6tSpU+R7auWv5/T+++979/Xr10/q1q0rr7/+ugQCGQAAgHN5FgOyu4WFmbfRyjf/lpOTU+TH5ebmyrZt26RPnz4++/VxWlpaiUVt27atxMfHS8+ePeWTTz7xeU4zAAXfs2/fvqW+Z3kQAAAAXC8hIcFkFDxbUXfy6ujRo5KXlycNGjTw2a+Pjxw5UuRrtNJfsmSJaeN/6623pHnz5iYI2Lhxo/cYfa0/71kRGAUAABC3jwJIT0/3aQLQjngl0U58+WlresF9Hlrh6+bRpUsX83l//vOf5frrr7f1nhWBDAAAwLkqaBRAdHS0z1ZcABAbGyvh4eGF7swzMjIK3cGXpHPnzrJv3z7v47i4uHK/p78IAAAAzlXJwwAjIiLMsL/U1FSf/fq4a9euZX4fHT2gTQP5swIF33PdunV+vae/aAIAAMAP48ePN8P0OnToYCpubd/XIYCjRo0yz0+aNEl+/PFHWblypXk8d+5cueyyy6RVq1amE+GqVatMfwDdPMaNG2eaA2bNmiW33nqrvP322/Lhhx/KZ599Ju4LAFJSWA0IAJy6GlAIzwQ4ZMgQOXbsmEyfPt2M509MTJSUlBRp0qSJeV735Z8TQCv9J554wgQFNWvWNIHAe++9JzfddJP3GL3TT05OlqeeekqefvppueKKK8x8Azrk0H3zAKxaJdEsBwgAjpN1+rTE3Htv5cwDMH26RNeoYf99zp6VmD/9KaBlDVb0AQAAwIWCtwkAAIDSsBiQbQQAAADnIgCwjSYAAABcKHgzAB9+qAMuq7oUAAB/5eZW3jUjAxCCAQAAAGVdDMiusMBNtRvsaAIAAMCFyAAAAJyLJgDbCAAAAM5FAGAbAQAAwLkIAEIwANiyRSQ8vKpLAQDwV14e18wBgjcAAACgNGQAbCMAAAA4FwGAbQwDBADAhcgAAACciwyAbQQAAADnIgAIwQDg3/929RSNAOBYllXVJYCjAwAAAEpDBsA2AgAAgHOxGJBtjAIAAMCFyAAAAJyLJgDbCAAAAM5FAFA5AcDMmTPlrbfekm+//VZq1qwpXbt2lVmzZknz5s29x1iWJdOmTZMlS5bI8ePHpVOnTrJw4UJp1aqVfyU7e9a/4wEA7kMAUDl9ADZs2CCPPPKIbNq0SVJTU+X8+fPSp08fOXXqlPeY2bNny5w5c2TBggWyZcsWiYuLk969e0t2drb9UgIAgKrLAKxdu9bn8SuvvCL169eXbdu2yfXXX2/u/ufOnSuTJ0+WQYMGmWNWrFghDRo0kNdee01GjhxZsaUHALgbGYCqGQWQmZlpftarV8/83L9/vxw5csRkBTwiIyOle/fukpaWVuR75OTkSFZWls8GAIBfAUB5NpeyfeZ6tz9+/Hj57W9/K4mJiWafVv5K7/jz08ee54rqVxATE+PdEhIS7BYJAAAEOgAYM2aM7Ny5U15//fVCz4UVmMJXg4WC+zwmTZpkMgmeLT093W6RAABuQwagcocBPvroo/LOO+/Ixo0bpVGjRt792uFP6d1+fHy8d39GRkahrED+JgLdComKYi0AAHDqWgD5OocHFH0AKicDoHfyeuevQwE//vhjadq0qc/z+liDAB0h4JGbm2tGD+iQQQAA4MAMgA4B1N78b7/9ttSuXdvbrq9t9zovgKb5k5KSZMaMGdKsWTOz6e+1atWSoUOHBuocAABuRQagcgKAxYsXm5833HBDoeGA999/v/l9woQJcubMGRk9erR3IqB169aZgAEAgApFAFA5AYA2AZRGswBTp041GwAACE6sBQAAcC6WA7aNAAAA4Fw0AYRgAHDFFSLh4VVdCgCAv/LyRHbuDOkAYNGiRfL888/L4cOHzWJ3Og1+t27dijxWR85pH7odO3aY2W/1eG0m79u3r/eY5cuXywMPPFDotdqnrkaNGhII7p0DEQAAG1avXm1GvOm6N9u3bzcVf//+/eXgwYNFHq9z5uiieCkpKWbtnB49esjAgQPNa/OLjo42AUX+LVCVf3BnAAAACMIMwJw5c2TEiBHy4IMPmsd69//BBx+Yu3yd3r4gfT4/HR6vw+nfffddadu2rU8nes+EepWBDAAAQNw+FXBWgUXpNFVfFJ3cTu/i8y96p/RxcYveFXThwgXJzs72LqTncfLkSWnSpImZYXfAgAGFMgQVjQAAAOB6CQkJPgvTFXUnr44ePSp5eXl+LXpX0F/+8hc5deqUDB482LuvRYsWph+ATrOva+xo6v+6666Tffv2BexvQxMAAEDc3gSQnp5u2uA9ilyjxuaid/lp5a4dALUJoH79+t79nTt3NpuHVv7t2rWTF198UebPny/uCgD0QkREVHUpAAD+ys113CiA6OhonwCgOLGxsRIeHl7obr+kRe/ydx7UvgNvvPGG9OrVq5RiXSTXXnttQDMANAEAAFBGERER0r59e59F75Q+LmnRO73z1ynzdT2dm2++udTP0YyCDhvMv7KuezIAAAAE4SiA8ePHy7Bhw6RDhw7SpUsXWbJkiRkCOGrUKPP8pEmT5Mcff5SVK1d6K//77rtP5s2bZ9L8nuyBLqKn/Q3UtGnTzHO6iJ52QtS0vwYACxculEAhAAAAOFcVBABDhgyRY8eOyfTp081Y/cTERDPGX3vwK92Xf06Av/71r3L+/Hmzoq5uHsOHDzcd/9SJEyfk4YcfNsGBBgU6PFDnD+jYsaMESphVlhV+KpFGPnrymQ8/LNH0AQAAx8nKzZWYJUskMzOzTO3q5aorPv5Yoi++2P77nDwpMTfeGNCyBiv6AAAA4ELB2wTQs6dIrVpVXQoAgL9OnxZZsqRyrhuLAYVgAAAAQGkIAGyjCQAAABciAwAAcC4yALYRAAAAnIsAwDaaAAAAcKGgzQCc69lPzrlsTCYAhIJzWVmV92FkAEIvAAAAoFQEALbRBAAAgAuRAQAAOBcZANsIAAAAzkUAYBsBAADAucLCyrcaYFiYuFXQBgD/7/+xFAAAOJEuBYDgF7QBAAAApaIJwDYCAACAcxEA2MYwQAAAXIgMAADAucgA2EYAAABwLgKA0AsAPv5YJCKiqksBAPBXbi7XzAmCNgAAAKBUZABsIwAAADgXAYBtjAIAAMCFyAAAAJyLDIBtBAAAAOciALCNAAAA4FwsBhR6AcCWLSLh4VVdCgCAv/LyuGZOELQBAAAApaIJwDYCAACAcxEA2MYwQAAAXIgMAADAucgA2EYAAABwLgKA0AsADhz4z+gOAICzWFZVlwCODgAAACgVGQDbCAAAAM5FAGAbowAAAPDTokWLpGnTplKjRg1p3769fPrppyUev2HDBnOcHn/55ZfLSy+9VOiYN998U66++mqJjIw0P9esWSOBRAAAAHB+BqA8m59Wr14tSUlJMnnyZNm+fbt069ZN+vfvLwcPHizy+P3798tNN91kjtPj//jHP8rYsWNNhe/xxRdfyJAhQ2TYsGHy9ddfm5+DBw+WzZs3S6CEWVZwddfIysqSmJgYqVkzU8LCoqu6OAAAP1lWlpw5EyOZmZkSHR0d0Loi88SJcn1Glr5PnTqSnp7u8z56F65bUTp16iTt2rWTxYsXe/e1bNlSbrvtNpk5c2ah4ydOnCjvvPOO7N2717tv1KhRpqLXil9p5a9lef/9973H9OvXT+rWrSuvv/66BITlh0WLFlmtW7e2ateubbbOnTtbKSkp3ucvXLhgTZkyxYqPj7dq1Khhde/e3dq9e7c/H2FlZmZqQGKJ6E8NTti4Bvwb4N8A/wac9W/gP9/j+n0eKJ664sSJTOvCBcv2duKEp87x3bQuK0pOTo4VHh5uvfXWWz77x44da11//fVFvqZbt27m+fz09dWqVbNyc3PN44SEBGvOnDk+x+jjxo0bW4HiV+6jUaNG8txzz8nWrVvNduONN8qtt94qe/bsMc/Pnj1b5syZIwsWLJAtW7ZIXFyc9O7dW7KzswMTvQAAUAHS09NNxsKzTZo0qcjjjh49Knl5edKgQQOf/fr4yJEjRb5G9xd1/Pnz5837lXRMce9ZEfwKAAYOHGjaMa666iqzPfvss3LxxRfLpk2bNJMgc+fONW0igwYNksTERFmxYoWcPn1aXnvttYCdAADAvS5cKP+mNP2ffysu/e8RVmCiGq0DC+4r7fiC+/19zyrrBKgRUHJyspw6dUq6dOliOjlopNKnTx/vMXoBu3fvLmlpacW+T05Ojmn3yL8BAFCZAUBZxcbGSnh4eKE784yMjEJ38B6aDS/q+GrVqskll1xS4jHFvWeVBAC7du0yd/1auWsnBh2moMMVPAX3N4WhHSa0I4dnS0hIsHMeAAAEXEREhBnOl5qa6rNfH3ft2rXI1+hNcsHj161bJx06dJDq1auXeExx71kh/O00oB0g9u3bZ23ZssV68sknrdjYWGvPnj3W559/bjpO/PTTTz7HP/jgg1bfvn2Lfb+zZ8+azhyeLT09nU6AVd6Bh41rwL8B/g04oxNgRkamdfasZXvLyPC/rMnJyVb16tWtZcuWWd98842VlJRkRUVFWQcOHDDPa904bNgw7/Hff/+9VatWLeuxxx4zx+vr9PX/+Mc/vMdoHaqdC5977jlr79695qd2Ety0aZMVKNXsRD9XXnml+V2jF+3sN2/ePDPMQendfnx8fJlTGMUNtahdm7UAAMCJNASorL7fdtL4BV/vLx2yd+zYMZk+fbocPnzY9HlLSUmRJk2amOd1X/45AXTCIH3+sccek4ULF0rDhg1l/vz5cscdd3iP0Tt9bVZ/6qmn5Omnn5YrrrjCzDegQw6Ddh6Anj17mrT9K6+8Yk5KT3DChAnmudzcXKlfv77MmjVLRo4c6dfYztq1mQcAAJxI5wHIzq6ceQCOHCnfZ2RlZUlcXGDLGqz8ygDo7EU625FW+Dq0T6OV9evXy9q1a01PRZ0ZacaMGdKsWTOz6e+1atWSoUOHBu4MAACuVRUZAFcGAD///LOZnlDTGxp5XXPNNaby17H+Su/8z5w5I6NHj5bjx4+b1IV2Yqit+XwAACoYAUAITgVMEwAAOFNlNgEcOlT+JoBGjWgCAADAUcgA2Of3KIDKogMNwsOruhQAAH/l5Yls3145140AIAQDAAAASqON2OXpyGcFVSN45bI9FTAAAHAuMgAAAMeiCcA+AgAAgGMRANhHEwAAAC4UtBmAzp11nYCqLgUAwF85OYwCcIKgDQAAACgNTQD20QQAAIALkQEAADgWGQD7CAAAAI5FAGAfTQAAALgQGQAAgGORAQjBAKBXL5GoqKouBQDAX6dOiSxeXDnXjQAgBAMAAABKw2JA9tEHAAAAFyIDAABwLJoA7CMAAAA4FgGAfTQBAADgQkGbAeh13RmJjq5e1cUAAPgpK+tMpV0zMgAhGAAAAFAaAgD7aAIAAMCFyAAAAByLDIB9BAAAAMciALCPJgAAAFwoeDMAKSkitWpVdSkAAP46fbrSrhkZgFAMAAAAKAUBgH0EAAAAx2IxIPvoAwAAgAuRAQAAOBZNAPaRAQAAOD4AKM8WSMePH5dhw4ZJTEyM2fT3EydOFHv8uXPnZOLEidK6dWuJioqShg0byn333Sc//fSTz3E33HCDhIWF+Wx33XVXiGQAPvxQJCKiqksBAPBXbi7X7P8MHTpUDh06JGvXrjWPH374YRMEvPvuu1KU06dPy1dffSVPP/20tGnTxgQQSUlJcsstt8jWrVt9jn3ooYdk+vTp3sc1a9aU0AgAAABwcBPA3r17TcW/adMm6dSpk9m3dOlS6dKli3z33XfSvHnzQq/RLEFqaqrPvhdffFE6duwoBw8elMaNG3v316pVS+Li4myXjyYAAIC4vQkgKyvLZ8vJySl32b744gtToXsqf9W5c2ezLy0trczvk5mZaVL8derU8dn/6quvSmxsrLRq1UqeeOIJyc7O9qt8ZAAAAK6XkJDgcw2mTJkiU6dOLdd1OXLkiNSvX7/Qft2nz5XF2bNn5cknnzRNCdHR0d7999xzjzRt2tRkAHbv3i2TJk2Sr7/+ulD2oCQEAAAAcXsTQHp6uk8FGxkZWexrNDCYNm1aie+7ZcsW81Pv3AuyLKvI/UV1CNSOfRcuXJBFixYVav/3SExMlGbNmkmHDh1M/4F27dpJWRAAAADE7QFAdHS0TwBQkjFjxpTa4/6yyy6TnTt3ys8//1zouV9++UUaNGhQauU/ePBg2b9/v3z88cellk0r/erVq8u+fftCIADQ3o7h4VVdCgCAv/LyQvqaxcbGmq002tlP2++//PJL04lPbd682ezr2rVrqZW/VuaffPKJXHLJJaV+1p49e8zr4uPjy3wedAIEADhWMM8D0LJlS+nXr59J1+tIAN309wEDBviMAGjRooWsWbPG/H7+/Hn53e9+Z4b8aSe/vLw8019At9z/G17573//2wz/02MOHDggKSkpcuedd0rbtm3luuuuK3P5gjcDAACAg4cBKq3Ex44dK3369DGPdTz/ggULJD8dEqhZAaVzBrzzzjvm99/85jc+x2k2QCcAioiIkI8++kjmzZsnJ0+eNB0Yb775ZtNxMdyPzDkBAADAsYJ9MaB69erJqlWrSimD5dN3IP/jomiFv2HDhnKXjSYAAABciAwAAMCxgr0JIJgFbwDwww86gLKqSwEACLa8ej4EAPbRBAAAgAsFbwYAAIBSkAGwjwAAAOBYBAD20QQAAIALkQEAADgWGYAqygDMnDnTrGiUlJTk3acTGOhKSQ0bNpSaNWuaWYt0jmIAANw0FXDIZgB0qcMlS5bINddc47N/9uzZMmfOHFm+fLlcddVV8swzz0jv3r3NVIe1a9cu+wfk5NgtGgAACEQGQOcevueee2Tp0qVSt25dn7v/uXPnyuTJk2XQoEFmjeIVK1bI6dOn5bXXXrPzUQAAFIsMQCUHAI888ohZeKBXr14++3XdYl2xyLPogYqMjJTu3btLWlpake+Vk5MjWVlZPhsAAGVBAFCJTQDJycny1VdfmSaAgrTyVw0aNPDZr49/0Jn9iulHMG3aNH+LAQAAnQArKwOQnp4u48aNMysb1ahRo9jjtGNgfto0UHCfx6RJk8wyiJ5NPwMAAARRBmDbtm2SkZEh7du39+7Ly8uTjRs3mvWNtaOfJxMQHx/vPUZfUzArkL+JQDcAAEJtOeCQCQB69uwpu3bt8tn3wAMPSIsWLWTixIly+eWXS1xcnKSmpkrbtm3N87m5uWbd4lmzZvlXspgYFgMCACfSWjUzs1I+inkAKikA0GF82rM/v6ioKLnkkku8+3VOgBkzZkizZs3Mpr/XqlVLhg4dWo5iAgCAoJ4JcMKECXLmzBkZPXq0HD9+XDp16iTr1q3zbw4AAADKgAxAFQYA69ev93msnf10JkDdAAAIJAIA+1gMCAAAF2IxIACAY5EBCMUA4MorRcLDq7oUAAB/5eXpuPFKuW4EAPbRBAAAgAsFbwYAAIBSkAGwjwAAAOBYBAD2EQAAAByLAMA++gAAAOBCwZsB6NRJVwqq6lIAAPyVk1NpowBYDCgUAwAAAEpBE4B9NAEAAOBCZAAAAI5FBsA+AgAAgGMRANhHEwAAAC4UvBmAvn1FoqKquhQAAH+dOiWyaFGlXDcyAKEYAAAAUAoCAPtoAgAAIECOHz8uw4YNk5iYGLPp7ydOnCjxNffff7+EhYX5bJ07d/Y5JicnRx599FGJjY2VqKgoueWWW+TQoUN+lY0AAADg+AxAebZAGjp0qOzYsUPWrl1rNv1dg4DS9OvXTw4fPuzdUlJSfJ5PSkqSNWvWSHJysnz22Wdy8uRJGTBggOTpUsxlRBMAAMCxgrkJYO/evabS37Rpk3TS2W1FZOnSpdKlSxf57rvvpHnz5sW+NjIyUuLi4op8LjMzU5YtWyZ///vfpVevXmbfqlWrJCEhQT788EPpq33oyoAMAABA3J4ByMrK8tk0xV5eX3zxhUn7eyp/pal83ZeWllbia9evXy/169eXq666Sh566CHJyMjwPrdt2zY5d+6c9OnTx7uvYcOGkpiYWOr7OiIDkN3+BgmLjq7qYgAA/JSdleW4a5aQkODzeMqUKTJ16tRyveeRI0dMJV6Q7tPnitO/f3+58847pUmTJrJ//355+umn5cYbbzQVv2YG9LURERFSt25dn9c1aNCgxPd1TAAAAEBlLQaUnp4u0fluOrWiLY4GBtOmTSvxfbds2WJ+age+wp9pFbnfY8iQId7f9a6+Q4cOJhh47733ZNCgQSWcS8nvWxABAABA3N4HIDo62icAKMmYMWPkrrvuKvGYyy67THbu3Ck///xzoed++eUXc7deVvHx8SYA2Ldvn3msfQNyc3PNCIP8WQBtJujatWuZ35cAAAAAP+jQO91Ko539tMPel19+KR07djT7Nm/ebPb5U1EfO3bMZCg0EFDt27eX6tWrS2pqqgwePNjs05ECu3fvltmzZ5f5fekECABwrGAeBtiyZUsznE878elIAN30dx2ul38EQIsWLcyQPqXD+Z544gnTgfDAgQOmM+DAgQNNwHH77bebY7QT4YgRI+Txxx+Xjz76SLZv3y733nuvtG7d2jsqoCzIAAAAHCuYhwGqV199VcaOHevtsa8T9ixYsEDy0yGBmhVQ4eHhsmvXLlm5cqWZMEjv+nv06CGrV6+W2rVre1/zwgsvSLVq1UwG4MyZM9KzZ09Zvny5eX1ZhVnaayCI6PALjW4OHcosc3sMACC4vscbNYoxlVqgvsc9dcVjj2VKZKT9z8jJyZIXXghsWYNV0GYAdNKjmjWruhQAAH+dOVN51yzYMwDBLGgDAAAASkMAYB+dAAEAcCEyAAAAxyIDYB8BAADAsQgA7CMAAAA4FgFACAYAH38sEhFR1aUAAPgrN5dr5gRBGwAAAFBZiwG5EQEAAMCxaAKwj2GAAAC4EBkAAIBjkQGwjwAAAOBYBAAhGABs26arIlV1KQAA/srL45o5QdAGAAAAlIYMgH0EAAAAxyIAsI9RAAAAuBAZAACAY5EBsI8AAADgWAQAIRgAHDwoEhZW1aUAAATz9LoEAPbRBwAAABcK2gwAAAClYTEg+wgAAACORROAfTQBAADgQmQAAACORQagkgKAqVOnyrRp03z2NWjQQI4cOWJ+tyzLPL9kyRI5fvy4dOrUSRYuXCitWrXyu2A5OX6/BADgMgQAldgEoJX54cOHvduuXbu8z82ePVvmzJkjCxYskC1btkhcXJz07t1bsrOzy1FEAABQ5U0A1apVMxV7QXr3P3fuXJk8ebIMGjTI7FuxYoXJELz22msycuTIiikxAAD/hwxAJWYA9u3bJw0bNpSmTZvKXXfdJd9//73Zv3//ftMU0KdPH++xkZGR0r17d0lLSyv2/XJyciQrK8tnAwDAnwCgPJtb+RUAaJv+ypUr5YMPPpClS5eaCr9r165y7Ngxbz8AveMvro9AUWbOnCkxMTHeLSEhwe65AACAQDQB9O/f3/t769atpUuXLnLFFVeYVH/nzp3N/rAC8/dq00DBfflNmjRJxo8f732sGQCCAABAWdAEUEXDAKOiokwgoM0Ct912m9mnd/vx8fHeYzIyMgplBfLTZgLdCoqJYS0AAHDq7HyZmZXzWQQAVTQRkLbf792711T42idAOwempqZ6n8/NzZUNGzaYZgIAACoafQAqKQPwxBNPyMCBA6Vx48bmzv6ZZ54xKfvhw4ebNH9SUpLMmDFDmjVrZjb9vVatWjJ06NByFBEAAFRpAHDo0CG5++675ejRo3LppZeadv9NmzZJkyZNzPMTJkyQM2fOyOjRo70TAa1bt05q165d4QUHAIDFgCqpCSA5OVl++uknk9r/8ccf5c0335Srr77a+7xmAXS2QJ0g6OzZsyb9n5iYWI7iAQDg3CaA48ePy7Bhw7wj3fT3EydOlPgarUuL2p5//nnvMTfccEOh53Vovj9YCwAAgADRJnDNnq9du9Y8fvjhh00Q8O677xb7Gr2Jzu/999+XESNGyB133OGz/6GHHpLp06d7H9esWdOvshEAAAAcK5hHAezdu9dU/NpUrk3iSufQ0SH03333nTRv3rzI1xWcbfftt9+WHj16yOWXX+6zX/vYFTUzr+MDAL0u1YK2dACA4pw/L/Lll84KALIKzEJb3BB1f3zxxRcm7e+p/JX2ndN9OkNucQFAfj///LO89957Zr6dgl599VVZtWqVGWqv8/RMmTLFrz53VLEAANdLKDALrVam2qetPHRenPr16xfar/tKmiE3P634tVL3rLHjcc8993iH3+/evdtMqvf111/7DMUvDQEAAMCxKioDkJ6eLtHR0d79Jd39a2Awbdq0Et9XV8RVRc2EW9oMufm9/PLLprKvUaNGofZ/D+1sr0PvO3ToIF999ZW0a9euTO9NAAAAELcHANHR0T4BQEnGjBlTao/7yy67THbu3GlS+AX98ssvJc6Q6/Hpp5+avgKrV68u9Vit9KtXr25m5iUAAAAgAGJjY81WGu3sl5mZKV9++aV07NjR7Nu8ebPZV5YZcpctWybt27eXNm3alHrsnj175Ny5cz5T8Qd0KmAAAKpSMM8D0LJlS+nXr59J1+tIAN309wEDBvh0AGzRooWsWbPG57XaKfGNN96QBx98sND7/vvf/zbD/7Zu3SoHDhyQlJQUufPOO6Vt27Zy3XXXlbl8QdsE0KOHtsFUdSkAAP7KyXHeKIBA0Z76Y8eOlT59+pjHt9xyiyxYsMDnGE3za1ag4MR72ldAZ98tKCIiQj766COZN2+enDx50nRgvPnmm03HxfDw8DKXLczSTwgiGvXoEImJEzMlMrJs7TEAgOCRk5Mls2bFmEqtrO3qduuKDh0ypVo1+59x/nyWbN0a2LIGK5oAAABwoaBtAgAAoDQsBmQfAQAAwLG0Db+MQ+qLfb1b0QQAAIALBW0G4OabRS6+uKpLAQDw18mTIrNmVc51IwMQggEAAAClIQCwjyYAAABciAwAAMCxyADYRwAAAHAsAgD7aAIAAMCFgjYD0ObKUxIdXfY5jQEAwSEr61SlfRYZgBAMAAAAKA0BgH0EAAAAxyIAsI8+AAAAuBAZAACAY7EYkH0EAAAAxyrvYj4XXLwYUPAGAKdOiVxECwUAOI5+fyPoBW8AAABAKcgA2EcAAABwLAIA+8ixAwDgQmQAAACORQbAPgIAAIBjEQCEYgAQFSVy8cVVXQoAgL/cPLbOQYI3AAAAoBRkAOwjAAAAOBYBgH0EAAAAxyIAsI9hgAAAuBAZAACAY7EYkH0EAAAARzcBhIWVL4Bwq+ANAGrV+s8GAHCW8+erugRwdAAAAEApyADYRwAAAHAsAgD7GAUAAIALkQEAADgWGQD7CAAAAI5FAGAfTQAAAATIs88+K127dpVatWpJnTp1yvQay7Jk6tSp0rBhQ6lZs6bccMMNsmfPHp9jcnJy5NFHH5XY2FiJioqSW265RQ4dOuRX2QgAAACOzgCUdwuk3NxcufPOO+UPf/hDmV8ze/ZsmTNnjixYsEC2bNkicXFx0rt3b8nOzvYek5SUJGvWrJHk5GT57LPP5OTJkzJgwADJy8sr8+eEWRpqBJGsrCyJiYmRzBMnJDo6uqqLAwCw8z1ep45kZmYG7HvcU1eEh2dKWJj9z7CsLMnLiwloWdXy5ctNpX3ixIlSymOZO389duLEid67/QYNGsisWbNk5MiRpqyXXnqp/P3vf5chQ4aYY3766SdJSEiQlJQU6du3b5nKRAYAACBuzwBkZWX5bFrpVoX9+/fLkSNHpE+fPt59kZGR0r17d0lLSzOPt23bJufOnfM5RoOGxMRE7zFlQQAAAHC9hIQEk1HwbDNnzqySa6KVv9I7/vz0sec5/RkRESF169Yt9hhHjgLwtEhoBAYAcB7P93dltDBrCr98H5Nl/puenu7TBKB33cXRDnrTpk0r8V217b5Dhw62SxVWYIEDvZYF9xVUlmOCOgDwdHJIaNy4qosCACjn97neTQeC3gFr57gjRxLK/V5xcXGmN32NGjXKdPyYMWPkrrvuKvGYyy67zHZZlN7Jx8fHe/dnZGR4swJ6jHYuPH78uE8WQI/REQeODQC0HUMjsdq1a5tIRiNJTc0UjM5CmRvP2a3n7cZzVpx3aP+99U5UK3/9Pg8Uray1vVwrwooIJmqUsfJXGizoFghNmzY1FXxqaqq0bdvW7NNz3LBhg+kEqNq3by/Vq1c3xwwePNjsO3z4sOzevduMIHBsAHDRRRdJo0aNCu3X/1lC+X+YorjxnN163m48Z8V5h65A3fnnp5W2PxV3VTh48KD8+uuv5qcO0duxY4fZf+WVV8rFF19sfm/RooXpc3D77bebG18dATBjxgxp1qyZ2fR3nUdg6NCh3ms7YsQIefzxx+WSSy6RevXqyRNPPCGtW7eWXr16OTcAAAAgVPzpT3+SFStWeB977uo/+eQTM8GP+u6778zQPo8JEybImTNnZPTo0SbN36lTJ1m3bp3JjHu88MILUq1aNZMB0GN79uxphhqGh4c7dx6AYucFCPAYzWDixnN263m78ZwV5+2uvzeCU9APA9SemFOmTCmxR2aoceM5u/W83XjOivN2198bwSnoMwAAAMCFGQAAAFDxCAAAAHAhAgAAAFyIAAAAABciAAAAwIWCOgBYtGiRmRZRZ3rSqQ8//fRTCSUbN26UgQMHmukydfanf/7znz7P6wANXXRCn69Zs6aZNGLPnj3iZDrb1bXXXmsmtKhfv77cdtttZhKMUD/vxYsXyzXXXOOd+a5Lly7y/vvvh/Q5F/W398xyFsrnreej55l/88zvHqrnDGcK2gBg9erV5oti8uTJsn37dunWrZv079/fTKcYKk6dOiVt2rSRBQsWFPm8zuk8Z84c87yuLKVfIr179/YumOREOp/1I488Ips2bTLzWJ8/f96saa3XIpTPW6e3fu6552Tr1q1mu/HGG+XWW2/1fvGH4jnnp+e0ZMkSEwTlF6rn3apVKzM3u2fbtWtXyJ8zHMgKUh07drRGjRrls69FixbWk08+aYUi/VOsWbPG+/jChQtWXFyc9dxzz3n3nT171oqJibFeeuklK1RkZGSYc9+wYYOrzlvVrVvX+tvf/hby55ydnW01a9bMSk1Ntbp3726NGzfO7A/V854yZYrVpk2bIp8L1XOGMwVlBkBXPtq2bZu5M8xPH6elpYkb6CpXuhxk/mugs6d17949pK6BZ/5rXczCLeetC4IkJyebrIc2BYT6OWvG5+abby60SEkon/e+fftMil+bMHXZ2O+//z7kzxnOE5SLAR09etR8SXrWPvbQx/o/jxt4zrOoa/DDDz9IKNDEx/jx4+W3v/2tJCYmhvx5axpYK/yzZ8+aVcDWrFkjV199tfeLPxTPWQOdr776yqS6CwrVv7Uu3LJy5Uq56qqr5Oeff5ZnnnnGrNGuzT2hes5wpqAMADy080zBCqPgvlAXytdgzJgxsnPnTvnss89ccd7Nmzc3S4GeOHFC3nzzTRk+fLjpExGq55yeni7jxo0zq5iVtGRrqJ239lXy0OVZNei74oorzIpwnTt3DslzhjMFZRNAbGysWdKw4N1+RkZGocg5VHl6DYfqNXj00UflnXfeMUtiagc5N5x3RESEWQO8Q4cOpke8dgCdN29eyJ6zNuPpOegIHl22VDcNeObPn29+95xbqJ13QVFRUSYQ0GaBUP1bw5kuCtYvSv3S0F7i+eljTaW5gbYd6pdF/mugfSP0C9TJ10DvdPTO/6233pKPP/7YnKcbzru4a5GTkxOy56zrk2uzh2Y9PJsGP/fcc4/5/fLLLw/J8y5I/8Z79+6V+Pj4kP1bw6GsIJWcnGxVr17dWrZsmfXNN99YSUlJVlRUlHXgwAErVGjv6O3bt5tN/xRz5swxv//www/mee0prL2D33rrLWvXrl3W3XffbcXHx1tZWVmWU/3hD38w57R+/Xrr8OHD3u306dPeY0LxvCdNmmRt3LjR2r9/v7Vz507rj3/8o3XRRRdZ69atC9lzLkr+UQChet6PP/64+ff9/fffW5s2bbIGDBhg1a5d2/vdFYrnDGcK2gBALVy40GrSpIkVERFhtWvXzjtULFR88sknpuIvuA0fPtw7ZEiHFOmwocjISOv66683XxhOVtT56vbKK694jwnF8/7973/v/bd86aWXWj179vRW/qF6zmUJAELxvIcMGWIqdL2BadiwoTVo0CBrz549IX3OcKYw/U9VZyEAAEDlCso+AAAAILAIAAAAcCECAAAAXIgAAAAAFyIAAADAhQgAAABwIQIAAABciAAAAAAXIgAAAMCFCAAAAHAhAgAAAMR9/j+Z1FMmTejSVAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ei = solver.get_TM_field(ji) # ordinary numpy solve, no chi -> vacuum\n", "print(\"incident |E| near center:\", abs(ei[Nx // 2, Ny // 2]))\n", "plt.imshow(np.real(ei), cmap=\"bwr\")\n", "plt.title(r\"incident field $\\mathrm{Re}\\,E_z$\")\n", "plt.colorbar()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "6", "metadata": {}, "source": [ "## 3. A differentiable objective\n", "\n", "`build_jax_field_solver` gives us the differentiable solve. `density_to_chi` places the\n", "design density into a full-grid susceptibility, and `absorption` runs the solve and\n", "returns the absorbed power." ] }, { "cell_type": "code", "execution_count": null, "id": "7", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING:2026-07-18 09:52:32,516:jax._src.xla_bridge:794: An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.\n" ] } ], "source": [ "chi_material = 4 + 0.1j # lossy dielectric (Im > 0 gives absorption)\n", "im_chi = float(np.imag(chi_material))\n", "dA = dl * dl\n", "ji_j = jnp.asarray(ji)\n", "\n", "# Build the differentiable field solve ONCE, outside the optimization loop.\n", "# field(source, chi) -> flattened Ez, differentiable in both arguments.\n", "field = build_jax_field_solver(solver)\n", "\n", "def density_to_chi(rho):\n", " \"Scatter the design density into a flattened full-grid susceptibility.\"\n", " chi_flat = jnp.zeros(Nx * Ny, dtype=jnp.complex128)\n", " return chi_flat.at[des_idx].set(rho.astype(jnp.complex128) * chi_material)\n", "\n", "def absorption(latent):\n", " \"Absorbed power; rho = sigmoid(latent) keeps the density in [0, 1].\"\n", " rho = jax.nn.sigmoid(latent)\n", " chi = density_to_chi(rho)\n", " Ez = field(ji_j, chi) # differentiable solve (flat)\n", " return 0.5 * omega * jnp.sum(rho * im_chi * jnp.abs(Ez[des_idx]) ** 2) * dA\n", "\n", "# jax.grad flows through `field`: the forward + adjoint solves are handled\n", "# automatically by the custom_linear_solve rule inside dolphindes.\n", "value_and_grad = jax.value_and_grad(lambda latent: -absorption(latent))" ] }, { "cell_type": "markdown", "id": "8", "metadata": {}, "source": [ "## 4. Optimize with optax\n", "\n", "Adam on the latent field. Each step is one forward solve (the objective) plus one adjoint\n", "solve (the gradient). We start from `latent = 0.5`. This is a natural way to do optimization with jax - but a wrapper with nlopt would be relatively easy to implement. \n", "\n", "Running this cell requires optax, which is not an explicit requirement of Dolphindes." ] }, { "cell_type": "code", "execution_count": null, "id": "9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "iter 0 absorption = 0.1356\n", "iter 25 absorption = 0.5393\n", "iter 50 absorption = 0.6120\n", "iter 75 absorption = 0.6561\n", "iter 99 absorption = 0.6804\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHFCAYAAAAOmtghAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAASppJREFUeJzt3Qd4VFXex/F/ek8glZIQCDUQWkCRLoIgsiAIK1hQ14qKLuK+9oJYsIsNVnZtWJEia8EC0gULobfQSQgpJIGEENLnfc4JGQkhmDKTO+X72efu3LlzZ+bkZsj8PNXFZDKZBAAAwEG4Gl0AAAAASyLcAAAAh0K4AQAADoVwAwAAHArhBgAAOBTCDQAAcCiEGwAA4FAINwAAwKEQbgAAgEMh3ACosQ8//FBcXFzk0KFDXDUANotwAwAAHArhBgD+QnFxsZSUlHCdADtBuAFQb++//7507dpVvL29JTg4WMaMGSO7du2qdM6BAwdkwoQJ0qxZM/Hy8pKIiAgZPHiwbN682XzO8uXL5dJLL5WQkBDx8fGRFi1ayNixYyU/P/8vy/DZZ59J7969xd/fX2/dunWT9957z/x4y5Yt5eabb67yPPV+aquwcuVK3fT28ccfywMPPCDNmzfX5d2xY4c+fvZrVvj+++/1Y19//bX52N69e+W6666T8PBw/fzY2Fh55513anhFAdSHe72eDcDpzZgxQx599FG59tpr9X5WVpZMmzZNB40//vhD2rZtq6/RlVdeKaWlpfLSSy/p0JKZmSnr1q2TEydO6MdVP54RI0ZI//79dVhq1KiRpKSkyA8//CBFRUXi6+tb7bV+8skn5ZlnnpGrr75aB5KgoCDZvn27HD58uM6/n0ceeUT/DP/+97/F1dVVoqKipHv37vLBBx/IrbfeWqUvkgox6mdUdu7cKX369NE/56uvvipNmjSRH3/8Ue677z79cz/11FNO/7kBrMoEADX0wQcfmNSfjYMHD+r7x48fN/n4+JiuvPLKSuclJSWZvLy8TNddd52+n5mZqZ83c+bMal97wYIF+pzNmzfX6vdx4MABk5ubm+n666+/4HnR0dGmm266qcrxgQMH6q3CihUrdDkGDBhQ5dw333xTP5aYmGg+lp2drX/WBx54wHxs2LBhpsjISFNOTk6l50+ePNnk7e2tnwPAemiWAlBn69evl9OnT1dp7lG1HJdddpn8/PPP+r5qqmrdurW8/PLL8tprr8mmTZukrKys0nNUM5Knp6fccccd8tFHH+lmrJpYunSprhG65557LPqbVM1h57r++ut1E5Oqqanw+eefS2FhofzjH//Q9wsKCvTPrZrmVG2T6qtTsamaHfX4r7/+atGyAqiMcAOgzlQTlNK0adMqj6m+NRWPq/4o6gt/2LBhulkqPj5ewsLCdDPNyZMn9Tkq/Cxbtkw376igou6r7Y033rhgGY4dO6ZvIyMjLfqbPN/PpELaqFGjZO7cuTpQKSroXHzxxdKpUyd9X/3MKsi89dZb4uHhUWmraLZSTVMArIc+NwDqTHX8VVJTU6s8dvToUQkNDTXfj46ONnfG3bNnj3z55Ze6b47qT6P6tSiqv43aVHDYsGGDDghTpkzRnY9VZ+TzUSFJOXLkiK4xqo7q7KxqWM6lgsbZ5aygAtn5qBqa+fPn6xoj1adG9SuaPXu2+fHGjRuLm5ubTJw4sdrapFatWlVbTgD1R7gBUGeqw60a1fTJJ5/I3//+d/NxFTTUyKdx48ad93nt2rWTxx9/XBYuXCgbN26s8rgKB7169ZIOHTrIp59+qs+pLtwMHTpUn68ChipPddRoqa1bt1Y6pkJWYmLiecNNddT7qRFUqmOxCjcqNKnO1BVUU9SgQYN001uXLl10UxuAhkW4AVBnakTTE088oUdL3XjjjfpLXjXLPP300/pLv2JUkAoVkydP1gFIjZ5SX/gq/KjjDz/8sD5H1d6oY2rElAoNqm+KGjWlDBky5IKhRb2/Gi2l+v+oMqjRUmrEkqqVUWVRVE3KDTfcIHfffbfuT6NGUqkmsoqan5pSQUr9rKrvUGBgoB6hpd7vbKoprV+/froW6q677tJlVM1v+/btk2+++Ub/nACsyIqdlQE4+GipCv/9739NXbp0MXl6epqCgoJMV111lWnHjh3mx9PT000333yzqUOHDiY/Pz+Tv7+/Pv/11183lZSU6HPWr19vGjNmjB7VpEYfhYSE6FFMX3/9dY3KNnfuXNNFF12kRyOp1+/evbsub4WysjLTSy+9ZIqJidHn9OzZ07R8+fJqR0vNnz+/2vfas2ePPkdtS5cuPe856hrdcsstpubNm5s8PDxMYWFhpj59+pieffbZGv08AOrORf2fNcMTAABAQ2K0FAAAcCiEGwAA4FAINwAAwKEQbgAAgEMh3AAAAIdCuAEAAA7F6SbxU4v1qWnhAwICqp1eHQAA2BY1c42aDFOtW+fqeuG6GacLNyrYXGj9GQAAYLuSk5P/cqFcw8PNrFmz5OWXX9YL76lVdWfOnKmnLD+fm2++WT766KMqxzt27Cg7duyo0fupGpuKi6OmTgcAALYvNzdXV05UfI/bbLiZN2+eXvFXBZy+ffvKu+++K8OHD9drwqi1Zc6l1mt54YUXzPdLSkqka9eulRbs+ysVTVEq2BBuAACwLzXpUmLo8gtq1d/4+Hi9mm+F2NhYGT16tMyYMeMvn7948WK9aN3BgwclOjq6xslPLXKXk5NDuAEAwE7U5vvbsNFSRUVFkpCQIEOHDq10XN1ft25djV7jvffe06sFXyjYFBYW6gty9gYAAByXYeEmMzNTSktLJSIiotJxdT8tLe0vn6/66Hz//fdy2223XfA8VQOkkl7FRmdiAAAcm6uttZ2pVrKatKd9+OGH0qhRI92EdSGPPPKIrsKq2FRHYgAA4LgM61AcGhoqbm5uVWppMjIyqtTmnEsFoPfff18mTpwonp6eFzzXy8tLbwAAwDkYVnOjQkmPHj1k6dKllY6r+3369Lngc1etWiX79u2TW2+91cqlBAAA9sbQoeBTp07VtS89e/aU3r17y5w5cyQpKUkmTZpkblJKSUmRuXPnVulIrEZaxcXFGVRyAABgqwwNN+PHj5esrCyZPn267iCswsqSJUvMo5/UMRV2zqb6zSxcuFDPeQMAAGBT89wYgXluAACwP3Yxzw0AAIA1EG4AAIBDIdwAAACHQrgBAAAWUVZmkmMnC+XAsTxx2tFSAADAPhSVlEl6boGkqS2nfEvNKah0LONkgRSXmqRFsK+sfnCQYWUl3AAA4ORKSssk/WShpJ44LUdVaDlxWgeX1JzTOrSoY5l5hVKT8dVqBSWT+l8Nl1OyBsINAAAOzGQySW5BiaQcPy1HdXg5rfdTzgQYdUzVvpTVILh4uLlIkyBvaRroU34b5C0RgWdug7ylSaC3hAV4iYebsb1eCDcAANh5eDmeXyxHjufLkeOnzbcqwOjbE6clr7CkRsFFBZVmQT7StJEKLD46tJRv5ceCfT3F1dWY2pjaINwAAGDjcguKJSlLhZZ8Sc4uDzDJZwWZ/KLSv3yNED9PadbIR5o18i6/DfKR5o3LA0zzRj4S6u9lF8GlJgg3AADYQJ8XVcOSlJ1v3pLNt6cl53TxX75GRKCXRDb2lcjGPjqsqP3mZ/bV5uPpJs6CcAMAQAMoKC7VYeVQ5ik5nJUvh7PP3Gbl62BT+hedXlTNS2Swr0Q19pEofVseZNSmamK8PZwnvPwVwg0AABZyuqhUhxYVYA5llQeZQ1nlIUZ13r0QT3dXPYS6YisPMD7SIqQ8yPh58ZVdU1wpAABqWQOjmowOqOCSeUoOntlUgFHzvVxIgLe7tAzxk+gQ3/It2E+HF7UfEeDtMH1ejEa4AQDgPDPtqiHTKrQcOFYeXvYfy9O3qgnpQvO9BPl4SMsQX2kZ6qeDTMtQFV78pFWInzTy9TBs7hdnQrgBADitU4Ul5uCy/1j5bXmYyZOC4rJqnxfg5S6twirCi5+0CvUt3w/xk8Z+ng36M6Aqwg0AwOHngUnPLTwTYPJkf0aeblJSt2rm3eq4u7ro5qKYMH+JCfWTmDC1+esAE+rvSQ2MDSPcAAAcghptpEYj7U0/KXsz/gwyqkbmQpPYqVFIrVWA0eHF78y+v+7M627wTLuoG8INAMDu+sMkH8+XxLSTsiddbXnmMKMWdzwfN1ULE1xeC9M63E9ah/qbgwzNSI6HcAMAsOnmpEQVYNJO6lsVaPZmnKy2P4yXu6u0CfeXtuH++lZtKsCoDr1qqDWcA+EGAGA4NQNvojnA5MqetDy9X93MvCqoqADTLiJA2kb4S3t1Gx6gZ+RVtTRwboQbAECDLjOgRiftTM2VXanlQWZ32slqJ7hTQaVVqJ8OLyrItG/iL20jAnSnXkIMqkO4AQBYbbbeXWm5suNoruw8miM7j5YHmcJq+sWo9Y/aqVqYJoHSoUl5mFH9Y7zcWVYAtUO4AQBYZNVqFV62p+ToMKNuVQff8y2X5OvppsNLbNNA6dD0zyCjJr8DLIFwAwColay8wvIAc7Q8yOxIydHrKJ1PqL+XdGoWKHHNA6Vj0yDp2CxQj1pimQFYE+EGAFCt9NwCXQuzPaU8zKj96vrHqGal8iATpG87Nw+S8EBvri4aHOEGAKCdyC+SLUdyZGvyCX27LeWEHop9PqqTrwownZoFmQNNMMsOwEYQbgDASUctqaHWm5JOyMak47I56YRekuBcalS1mismToWY5kES1yxQNy0FeNM/BraLcAMATlIro0LMxsMnJOHwcdly5ITkF5VWOU+tZt0lspF0iQySrlGNdK2MrydfFbAvfGIBwAHXWFIT4m1KLg8zm5KOn7dWRq1s3a1FI+ke1Ui6t2gs3aIasRQBHALhBgDsXGZeoW5eUiFG3W49ckJOnadWRq1sHR/dWOJbNJYe0Y11cxMT4cEREW4AwM5qZdRikappqbyZ6fh5h2H7e7lL16ggHWTURq0MnAnhBgBsWHFpmR5+/fvBbL39cShbcgtKKp3jojr9hvnrENNdNTO1oFYGzo1wAwA2Nopp+9FcWb8/S9YfyJINh7KrdPxVM/yqENND1cpEq0DTmNl9gbMQbgDA4DCjFpH89UCW/HqgvHYmr7ByzUwjXw+5qGWw9GoVLBe3CpaOTQPF3c3VsDIDto5wAwANqKikTLaZm5lUzcxxOXlOmAn0dpdeMSHSW22tQ/SK2CxXANQc4QYArKiguFSPYPrtYJYONKoTcEFx5VWxA1SYaRUil8QEyyUxIXpBSUYxAXVHuAEAC9fMqCHZv+zL1M1Mm5NPSFFp5TDT2NdDNy+ppibCDGB5hBsAqKdDmadk1Z5jsmbvMd0R+Nw5ZsIDvHSIUYFG9ZtR88u4qCFOAKyCcAMAdZhrZnPycVm6M0OW7UqXfRl5lR4P8fOUvm1CpU/rEN13Ri1pQJgBGg7hBgBq2HdGNTX9uCNNft6VIVmniv78Q+rqopuYBrQLk/5tQ/VoJjoAA8Yh3ABANXLyi2XlngwdaFYmHqs034wa0TSoQ7gMiY2Qge3DJJBVsgGbQbgBgLMczjoly3ZlyLKd6fL7oWzdBFWhaZC3DO0YIcM6NZGLWgWLB3PNADaJcAPAqRWWlOoh2it2H9O1NAeOVV49u12Ev1x+JtB0bh5E3xnADhBuADidlBOnZcXuDFmZmCG/7MuS08WllfrPqFFNqrlJbS1CfA0tK4DaI9wAcHhlZSbZdGZ0kwo1ieknqwzVvrR9mAxqHy5924bSfwawc4QbAA67ZpNqbvp+e5ruEJxxstD8mKuL6BW0VYdgFWrU6CaGagOOg3ADwGGYTGr+mRPy1aYU+XZrqmSfNVxbLXEwuEO4XBYbIQPahkojX09DywrAegg3AOxeUla+DjSLN6fIwcxTlZY5GNqxiVzRuYn0bR0qnu6spA04A8INALudg+bbbUflq40psuHwcfNxHw83GdYpQsbER0rf1iHiznBtwOkQbgDYjVOFJXoyvW+3HtWzBFcsSKn60KjlDq6Ob65ravy8+NMGODPD62hnzZolrVq1Em9vb+nRo4esWbPmgucXFhbKY489JtHR0eLl5SWtW7eW999/v8HKC6BhnSwolv9tTpE7P94g8c8slXs+26g7Catg06FJgDx6ZQdZ/8hg+fjWXjKmeyTBBoCxNTfz5s2TKVOm6IDTt29feffdd2X48OGyc+dOadGixXmfc80110h6erq899570qZNG8nIyJCSkpIGLzsA63UK3p12UtfQqHloEg4fl5KzZgmODvGVK+KayFVdm0vHZoH8GgBU4WJSf0kM0qtXL4mPj5fZs2ebj8XGxsro0aNlxowZVc7/4YcfZMKECXLgwAEJDg6u03vm5uZKUFCQ5OTkSGAgfxgBW6mdUYtSlgeaY5KWW1Dp8dZhfjI8rqkM79yEYduAk8qtxfe3YTU3RUVFkpCQIA8//HCl40OHDpV169ad9zlff/219OzZU1566SX5+OOPxc/PT0aNGiXPPPOM+Pj4NFDJAVhqluAlW1Pl593psuFQ5doZbw9X6R0TIpe2L5+HJjrEj4sOoMYMCzeZmZlSWloqERERlY6r+2lpaed9jqqxWbt2re6f89VXX+nXuPvuuyU7O7vafjeqj47azk5+AIxx7GShLNmWKt9sOVpphJMSE+qnV9dWgaZXq2Dx9nDj1wSgTgwfUnDurKCqlay6mULLysr0Y59++qmumlJee+01GTdunLzzzjvnrb1RzVtPP/20lUoPoCZLH6zdlylz1x+S5bszpKKCRv0zv7hlsAyPa6JnCqZ2BoDdh5vQ0FBxc3OrUkujOgifW5tToWnTptK8eXNzsKnoo6MC0ZEjR6Rt27ZVnvPII4/I1KlTK9XcREVFWfRnAVBVzuliWZhwRD7+9XClifW6RjWSUV2byYjOTaVJkDeXDoDjhBtPT0899Hvp0qUyZswY83F1/6qrrjrvc9SIqvnz50teXp74+/vrY3v27BFXV1eJjIw873PUcHG1AWgYydn58t7ag/LlhmTJLypfbTvAy13G9oiUGy6Jljbh5f92AcAhm6VUjcrEiRN1J+HevXvLnDlzJCkpSSZNmmSudUlJSZG5c+fq+9ddd53uPPyPf/xDNzWpPjf/93//J7fccgsdigGDbUw6Lv9dc0B+2J5mbnpqF+EvN/ZuKWO6N2f+GQANxtBwM378eMnKypLp06dLamqqxMXFyZIlS/QEfYo6psJOBVVbo2p27r33Xh2IQkJC9Lw3zz77rIE/BeC8CopLdefgT39L0gtWVhjQLkxu799K+rUJZbVtAM41z40RmOcGqL/9x/Lks9+SZEHCEd23RvF0c5WrujWT2/rHSPsmAVxmAM43zw0A+3I465R8ty1VD+XenvLnlAqRjX3kul4t5O89oiQsgP5tAIxHuAFwXqpSNzH9pF6g8vvtlQONWqhyUPtw3UFYNUG5qQMAYCMINwDMSkrLZP2BLFm2M12W7crQswhXUAFGzRp8ZeemMqxThIT4U0sDwDYRbgAnp2poVK3Mok1HdOfgzLwi82Ne7q66U/CQjhEytCOBBoB9INwATiqvsEQ+/y1J5m1Iln0ZeebjIX6eMrRThAzuECF924SKjyfLIACwL4QbwMmo0U0frTsk7/9yUE7kF5traC7vGCFXxzeX/m3DxMPN1ehiAkCdEW4AJ3Eiv0j+u+agDjYnC0v0sVahfnLHgBgZ0aWpBHp7GF1EALAIwg3g4ApLSuXj9YflreX7zHPSqJmDJ1/WVq/vxEgnAI6GcAM4cEdhNS/Niz/sluTs8lFP7SMC5P7L2+nOwa4M3wbgoAg3gAOGmhWJGfLmz/vMSyKEB3jJA0PbybgeUdTUAHB4hBvAQZSVmeSHHWny9vJ9sjO1fMI9X083uXNAa7l9QCvx9eSfOwDnwF87wM4lZ+fLjzvS5Is//hzS7efpJjf0jpbb+sWwJAIAp0O4AeyMmkV4V+pJWbYrXYea3WknzY8FervLP/q2kn/0bSmNfD0NLScAGIVwA9i4rLxCvSTC5qQTug/N9qM5UlBcZn5cjXa6uGWwXBHXRM9TE8CQbgBOjnAD2HBT00870mXD4WwpM1V+PMDbXS6JCZFhnZrI4A7h0tiPWhoAqEC4AWxolNNPO9PlreV7K63ArcQ2DZSe0Y2lW1Qj6RrVSGJC/RjKDQDVINwANkA1Nz3/3S75/VB2paYmtcbT0E5NpHkjH6OLCAB2g3ADGNz89NKPiXo17oo1nm7r30pu7RcjwTQ1AUCdEG4Ag5qgFiQckWlf75BTRaXi4iJydfdIPdFeM2ppAKBeCDdAA8vJL5ZHv9qml0ZQVF+aaaM6SVzzIH4XAGABhBugAa3fnyVTv9wsqTkF4u7qotd5mjSwNUsiAIAFEW6ABnCqsERe+SlRPlx3SEwmkZYhvjJzQnc9+gkAYFmEG8DK1CKWj3+1XVJOlK/MfU3PSHlqZCfx8+KfHwBYA39dASvOLPzMtztl8ebykVBqOPfzV3eWge3CuOYAYEWEG8CCSstMsm5/ph4J9cP2NCksKRNXF9HrPU29vB21NQDQAAg3QD3kFhRLUla+HDmeL1uO5MjiTSm6s3CFTs0C5bkxnelbAwANiHAD1EBRSZnsST8p21Ny9MKVanmEg5mnJOd0cZVzg3w8ZFTXZjKuR6R0iQwSFzWJDQCgwRBugL+Yk+alH3fL/A1HpKj0z5W4zxbi5ylRwb56BNTlHZvI4Nhw8fZw47oCgEEIN0A1Mwh/szVVpn+zUzLzCs01MnHNAyWuWZB0ah4kbcP9pUWwL/1oAMDGEG6A86z39Pji7bJqzzF9v3WYnzwzOk56x4TQxAQAdoBwA5xlw6FsufH93yW/qFQ83Vxl8mVt5M6BMeLlTjMTANgLwg1wRmLaSbnlwz90sFHrPb00rovEhPlzfQDAzhBuABE9lPvG93+T3IIS6RHdWD6+tZf4eFJbAwD2yNXoAgBGyz5VpJui0nMLdSfh927qSbABADtGuIE4+4KW//jgdzlw7JQ0C/KWubdeLI18PY0uFgCgHgg3cFplZSa59/NNembhRr4eOtg0DfIxulgAgHoi3MBpzVy2R5bvzhAvd1d5/+aLpE14gNFFAgBYAOEGTumnHWny5vJ9en/G1Z0lvkVjo4sEALAQwg2czv5jeTL1yy16/+Y+LeXq+EijiwQAsCDCDZxKXmGJ3Plxgr69qGVjeWxErNFFAgBYGOEGTrVe1L++3CL7MvIkItBL3rk+Xjzc+CcAAI6Gv+xwCirYvPD9bvlhR5p4uLnIrOt7SHiAt9HFAgBYAeEGTmHWyv3y7uoDev+5MZ31LMQAAMdEuIHD+/jXw/Lyj4l6/7ErY+WanlFGFwkAYEWEGzi0/21OkSf/t13vTx7URm4fEGN0kQAAVka4gcNatjNdHvhyi5hMIjf2jpYHhrYzukgAgAbAquBwyM7D/1lzQHcgLjOJjOneXKaN7CQuLi5GFw0A0AAIN3C4hTAfXLhVvtuaqu+PjY+UF8Z2FldXgg0AOAvCDRzGwcxTcufHG2RPep64u7rIUyM7yg2XRFNjAwBOhnADh7AiMUPu+3yTnCwokbAAL5l9fbz0bBlsdLEAAAYg3MDu+9e8t/agPL9kl+5fo+avUcEmPJAJ+gDAWRk+WmrWrFnSqlUr8fb2lh49esiaNWuqPXflypW6ieHcbffu3Q1aZtiGopIyeWTRNnn2u/JgM75nlHx++yUEGwBwcobW3MybN0+mTJmiA07fvn3l3XffleHDh8vOnTulRYsW1T4vMTFRAgMDzffDwsIaqMSwFdmniuSuTxLkt4PZovoKP3plrNzarxX9awAA4mJS9foG6dWrl8THx8vs2bPNx2JjY2X06NEyY8aM89bcDBo0SI4fPy6NGjWq03vm5uZKUFCQ5OTkVApIsB9qRe+r3l4r+4+dEn8vd3nr2u4yqEO40cUCAFhRbb6/DWuWKioqkoSEBBk6dGil4+r+unXrLvjc7t27S9OmTWXw4MGyYsWKC55bWFioL8jZG+zbtK936GDTJNBbFt3dh2ADALCNcJOZmSmlpaUSERFR6bi6n5aWdt7nqEAzZ84cWbhwoSxatEjat2+vA87q1aurfR9VA6SSXsUWFcW6Qvbsmy1HZUHCEd0U9ea13aVdRIDRRQIA2HOfm5KSEnnuuefklltusVhIOHfWWNVKVt1MsirMqK1C7969JTk5WV555RUZMGDAeZ/zyCOPyNSpU833Vc0NAcc+HTmeL49+tc28TtTFrRjqDQCoZ82Nu7u7vPzyy7rGpb5CQ0PFzc2tSi1NRkZGldqcC7nkkktk79691T7u5eWl2+bO3mB/SkrL5P55m/U8Nt1bNJL7Brc1ukgAAEdplhoyZIju2Ftfnp6eeuj30qVLKx1X9/v06VPj19m0aZNuroJjm7Vyv/xx6LjuQPzG+O7i7mb4LAYAAEcZCq6Gaqumnu3bt+tw4ufnV+nxUaNG1fi1VHPRxIkTpWfPnrqJSfWnSUpKkkmTJunH1fukpKTI3Llz9f2ZM2dKy5YtpVOnTrpD8ieffKL736gNjmtT0nF54+fy2rlnRneSFiG+RhcJAOBI4eauu+7St6+99lqVx1Rfmdo0WY0fP16ysrJk+vTpkpqaKnFxcbJkyRKJjo7Wj6tjKuxUUIHmX//6lw48Pj4+OuR89913cuWVV9b2x4AdefGH3VJaZpJRXZvJmO6RRhcHAGDjDJ3nxgjMc2NfNiefkNHv/KIXwlz94CBp1sjH6CIBABx5npuCgoL6PB34S++u2q9vr+rWnGADAKiRWocb1ez0zDPPSPPmzcXf318OHDigjz/xxBPy3nvv1fblgGodzDwlP+woH01358AYrhQAwDrhRs1z8+GHH8pLL72kRzxV6Ny5s/z3v/+t7csB1Zqz+oCoRtPBHcKZrA8AYL1wo0YuqVFN119/vZ6npkKXLl1YnRsWk3GyQBZuPKL37xzYmisLALBeuFEjldq0aVPleFlZmRQXF9f25YDz+vCXQ1JUUibxLRrJRS0bc5UAANYLN2r49Zo1a6ocnz9/vl7QEqivkwXF8vGvh821NtUtxwEAgEXmuXnqqaf0xHuqBkfV1qgFLBMTE3Vz1bffflvblwOq+Pz3JL3MQkyYn1weW/OlOAAAqFPNzciRI2XevHl6sj31X9RPPvmk7Nq1S7755hu5/PLLuaqoF9UU9d7ag3r/zgEx4qqW/wYAwJo1N8qwYcP0BljaisQMSc8tlFB/LxndvTkXGABg/Zqbxx57TC9umZ+fX/t3A/7CwoTyEVJXxzcXL/c/R+MBAGC1cJOQkCBjx46Vxo0b68Uu1eKWP/zwg+Tl5dX2pYBKsk8V6ZobZWw8a0gBABoo3Kggc/z4cVm5cqVcddVVsmnTJr0AZnBwsFxyySV1LAYg8vXmFCkuNUlc80Bp3ySASwIAaLg+N2ryPlVrowKNqsEJCAiQxYsXy/795esAAXWx4MykfeOotQEANGTNzezZs2XChAnStGlT6d+/v/z000/6VjVXHTt2rD5lgRNLTDsp21NyxcPNRUZ1oyMxAKABa27uueceCQsLkwceeEAmTZr0l8uOAzVRsdTCoPbhEuz355plAABYveZGTdqn1pX64osvJDw8XHr16iUPPfSQfP/993QqRp2UlJbJoo0pen9cDzoSAwAauOZm9OjRelNycnL0UgwLFizQnYvVpH6FhYX1LBKczZq9mZKZV6hrbC5tH250cQAAztihODs7W1atWqVHTKlt+/btEhISIgMHDrR8CeE0HYmv6tZMPN1rXZkIAED9wk2XLl1k586deqTUgAED5Pbbb5dLL71U4uLiavtSgOTkF8vSHen6SjC3DQDAkHBzxx13EGZgMd9sPSpFpWXSoUmAdGpG53QAgAHhZvLkyeZ9k8mkb1VfG6AuvtqUYq614XMEALCEOnVwmDt3rnTu3Fl8fHz0ppqqPv74Y4sUCM4jOTtfEg4fF5WNR3VrZnRxAADOWnPz2muvyRNPPKFrcPr27atrb3755Rc9501mZqbcf//91ikpHM7XW47q2z6tQyQi0Nvo4gAAnDXcvPXWW3qW4htvvNF8TA0D79Spk0ybNo1wgxpRoXjxmSapq5iRGABgZLNUamqq9OnTp8pxdUw9BtTErtSTsjcjTw/9viKuCRcNAGBcuGnTpo18+eWXVY7PmzdP2rZta6lywcH9b3N5rc3gDuES6O1hdHEAAM7cLPX000/L+PHjZfXq1brPjRrhsnbtWvn555/PG3qAc5WVmcz9bdTEfQAAGFpzM3bsWPntt98kNDRUFi9erNeaUvu///67jBkzxqKFg2P6/VC2pOYUSIC3O8stAABsY/mFHj16yCeffGL50sCpmqSGxzURbw83o4sDAHAwdQo3paWl8tVXX8muXbt0s1RsbKweMeXuXqeXgxMpLCmVJdvS9P5oRkkBAKyg1mlELZKpgkxaWpq0b99eH9uzZ4+EhYXJ119/rSf3A6qzKvGY5JwulohAL+kVE8KFAgAY3+fmtttu03PaHDlyRDZu3Ki35ORkPUuxWncKuJD/nelIPLJLM3FzZdkOAIAN1Nxs2bJFNmzYII0bNzYfU/vPPfecXHTRRZYuHxzIyYJiWbazfAXw0d2bG10cAICDqnXNjWqKSk8v/4I6W0ZGhp4DB6jO99vSpLCkTGLC/FgBHABgO+Hm+eefl/vuu08WLFigm6bUpvanTJkiL774ouTm5po34GwLEo7o23E9WAEcAGA9Lia1yE8tuLr+mYfUSCml4iXOvq/21agqW6NCV1BQkOTk5EhgYKDRxXEahzJPyaWvrBTVzWbdw4OlSRALZQIArPP9Xes+NytWrKjtUwBZtLG81qZf2zCCDQDAqmodbgYOHGidksChl1tYuDHF3CQFAIBN9bkBauvXA1mScuK0Xm5haMcILiAAwKoIN2iwjsSjujZjuQUAgNURbmD1uW2WbE/V+zRJAQAaAuEGVrVkW6oUFJdJ6zA/6RbViKsNALA6wg0aaG6bKPNUAQAAGD5aqnv37jX+YlJrTQEVc9v8cei4ntvm6niWWwAA2FC4GT16tHm/oKBAZs2aJR07dpTevXvrY7/++qvs2LFD7r77buuVFHZn4Zm5bQa0C5OIQCbtAwDYULh56qmnKq0KrpZfeOaZZ6qco1YHB5Ti0jL5ckP554GOxAAAm+5zM3/+fLnxxhurHL/hhhtk4cKFlioX7Nz329MkPbdQwgK8ZGjHJkYXBwDgRGodbnx8fGTt2rVVjqtj3t40PaDcR+sO6dvre7UQT3f6rQMAbHj5BbX691133SUJCQlyySWXmPvcvP/++/Lkk09ao4ywM1uPnJCEw8fFw81FruvVwujiAACcTK3DzcMPPywxMTHyxhtvyGeffaaPxcbGyocffijXXHONNcoIO/PhmVqbv3VpJuEB1OYBAGw83CgqxBBkcD6ZeYXy7ZbyGYlv6tOSiwQAaHB16gxx4sQJ+e9//yuPPvqoZGdnm+e3SUkpX/m5NtSw8latWun+Oj169JA1a9bU6Hm//PKLuLu7S7du3Wr9nrCez39LkqLSMj0bMTMSAwDsItxs3bpV2rVrJy+++KK8/PLLOugoX331lTzyyCO1eq158+bpPjyPPfaYbNq0Sfr37y/Dhw+XpKSkCz4vJydHj9gaPHhwbYsPKw///uS3w3r/H32ptQEA2Em4mTp1qtx8882yd+/eSqOjVChZvXp1rV7rtddek1tvvVXPnaP67cycOVOioqJk9uzZF3zenXfeKdddd515EkHYhh/OGv49PK6p0cUBADipWoebP/74Q4eLczVv3lzS0tJq/DpFRUV6xNXQoUMrHVf3161bV+3zPvjgA9m/f3+liQUvpLCwUHJzcyttsG5HYoZ/AwDsKtyo2przBYTExEQJCwur8etkZmZKaWmpREREVDqu7lcXklRtkRqt9emnn+r+NjUxY8YMCQoKMm+qZgiWtyWZ4d8AADsNN1dddZVMnz5diouL9X21oKbqI6NCx9ixY2tdgHMX5DSZTOddpFMFIdUU9fTTT+s+PzWl+gGpPjoVG0tEWMebP+/VtyO7MvwbAGBnQ8FfeeUVufLKKyU8PFxOnz4tAwcO1DUtqv/Lc889V+PXCQ0NFTc3tyq1NBkZGVVqc5STJ0/Khg0bdMfjyZMn62NlZWU6DKlanJ9++kkuu+yyKs/z8vLSG6xba/Pz7gy9+vfkQW241AAA+wo3gYGBeqmF5cuX6+HfKmDEx8fLkCFDavU6np6eeuj30qVLZcyYMebj6r6qHTrf+27btq3KMHJVjgULFujh5DDG68v26Nsx3SMlJsyfXwMAwP4m8VNULcn5akpqO/Jq4sSJ0rNnT13zM2fOHN3ENWnSJHOTkpo7Z+7cueLq6ipxcXGVnq9qj1QfoHOPo+GoZRZWJh4TN1cXuW8wtTYAADsNNz///LO8/vrrsmvXLt0/pkOHDnq+mtrW3owfP16ysrJ0H57U1FQdUpYsWSLR0dH6cXXsr+a8gbFmnqm1GRvfXKJD/Ph1AAAM52JSnVZq4e2335b7779fxo0bZ55nRi2cqZqG1Lw1Ff1hbJUa6aVGTanOxaqpC3W34VC2jPv3enF3dZEV/7pUooJ9uZwAAMO/v2tdc6OGVqtam7NDzH333Sd9+/bVHYptPdzA8n1t/t4zkmADALDfoeAqOV1xxRVVjqvJ95ggz3n8eiBLftmXJR5uLnIPI6QAAPYcbkaNGqXXkTrX//73Pxk5cqSlygUbployX/upvNbmmp5REtmY5igAgO2oUbPUm2++ad5Xa0Cp5qeVK1dW6nOjVul+4IEHrFdS2IzvtqXK74eyxcvdlVobAIDNqVGH4prOIaNGTh04cEBsGR2K6ye/qESGvLpKjuYUyJQhbWXKkJrPFg0AgM10KD548GCdCwPHMnvlfh1smjfykUkDWxtdHAAA6tfnRq0nFRMTIzt37qzN0+AgkrLy5d3V5TVzT/wtVrw93IwuEgAA9Qs3Hh4eUlhYeN6FLeH4nvlupxSVlEnfNiEyrFMTo4sDAIBlRkvde++98uKLL0pJSUltnwo7tmrPMVm6M10vszBtZCcCLgDAZtV6Er/ffvtNL7+gVuHu3Lmz+PlVnnJ/0aJFliwfbICqrXn6mx16/6beLaVtRIDRRQIAwHLhplGjRjJ27NjaPg127P1fDsqBY6ck1N9Tplze1ujiAABg2XDzwQcf1PYpsGPJ2fnmxTEfuqKDBHp7GF0kAAAsvyq4cuzYMUlMTNR9L9q1aydhYWF1fSnYKDUF0lNf75CC4jLp1SpYxvWINLpIAABYvkPxqVOn5JZbbpGmTZvKgAEDpH///tKsWTO59dZbJT8/v7YvBxv2w/Y0Wb47Q68f9dyYznQiBgA4ZriZOnWqrFq1Sr755hs5ceKE3tS6UuoYyy84jpMFxTLtTCfiuwa2ljbh/kYXCQAA6zRLLVy4UBYsWCCXXnqp+diVV14pPj4+cs0118js2bNr+5KwQa/+tEfScwulZYiv3M2q3wAAR665UU1PERERVY6Hh4fTLOUgth45IR+tP6T3nx3dmZmIAQCOHW7USuBPPfWUFBQUmI+dPn1ann76afMq4bBfpWUmefSrbaKWUx3drZn0axtqdJEAALBus9Qbb7whV1xxhURGRkrXrl11J9PNmzeLt7e3/Pjjj7V9OdiYBQnJsj0lVwK83eWxER2NLg4AANYPN3FxcbJ371755JNPZPfu3Xq48IQJE+T666/X/W5g352IX/4xUe//c3BbCQvwMrpIAAA0zDw3KsTcfvvtdXkqbNg7K/ZLZl6RxIT6yY29WxpdHAAAGqbPzUcffSTfffed+f6DDz6ol2To06ePHD58uG6lgOEOZ52S99ce1PuPjYgVT/dafzQAALAJtf4Ge/75583NT+vXr5e3335bXnrpJQkNDZX777/fGmVEA5ixZLcUlZZJ/7ahclmHcK45AMB5mqWSk5OlTZs2en/x4sUybtw4ueOOO6Rv376V5r6B/Vi/P0t+2JEmri4ij4/oyEzEAADnqrnx9/eXrKwsvf/TTz/JkCFD9L4aLaWGhMP+hn5P/3an3r++V7S0bxJgdJEAAGjYmpvLL79cbrvtNunevbvs2bNHRowYoY/v2LFDWrakE6q9WZhwRHal5kqgt7vcf3k7o4sDAEDD19y88847erI+tSq4WoohJCREH09ISJBrr722/iVCgykpLZO3V+zT+/de1laC/Ty5+gAAu+diUhPVOJHc3FwJCgqSnJwcCQwMFGf2v80p8s8vNutQs/ahQeLrWaeZAQAAsKnv7zp9mx0/flzee+892bVrl+582qFDB7nlllskODi4rmVGAysrM8k7Z2ptbunbkmADAHDeZqlVq1bpvjVvvvmmDjnZ2dny1ltvSatWrfRjsA/LdqXLnvQ8CfByl4lM2AcAcCC1rrm55557ZPz48TJ79mxxc3PTx0pLS+Xuu+/Wj23fvt0a5YQFqZbId1bu1/sTe0dLkI8H1xcA4Lw1N/v375cHHnjAHGwUtT916lT9GGzfuv1ZsiX5hHi5u8ot/VoZXRwAAIwNN/Hx8bqvzbnUsW7dulmqXLCiir42117cQkL9WRwTAOCEzVJbt2417993333yz3/+U/bt2yeXXHKJPvbrr7/qIeIvvPCC9UoKi9iYdFzX3Li7usjtA2K4qgAA5xwK7urqqkdF/dWp6hzV/8aWOftQ8Ns++kOW7cqQa3pGykvjuhpdHAAAjBkKfvBg+WrRsG8HjuXpYOPiIjJpYGujiwMAgFXUKNxER0db593RoOb9kaxvB7UPl5gwf64+AMAh1XlK2p07d0pSUpIUFRVVOj5q1ChLlAsWVlRSJgsSjuj9CRdFcX0BAA6r1uHmwIEDMmbMGNm2bVulfjhqX7H1PjfOPGlf1qkiCQ/wkss6hBtdHAAAbGcouBoppWYjTk9PF19fX70a+OrVq6Vnz56ycuVK65QS9fb570n69u89I8Xdrda/dgAAHLfmZv369bJ8+XIJCwvTo6jU1q9fP5kxY4YeJr5p0ybrlBR1lpydL2v3Zer98T1bcCUBAA6t1v8Jr5qd/P3LO6OGhobK0aNHzZ2OExMTLV9C1NuXG5JFtR72axMqLUJ8uaIAAIdW65qbuLg4PalfTEyM9OrVS1566SXx9PSUOXPm6GOwLSWlZTrcKBMupiMxAMDx1TrcPP7443Lq1Cm9/+yzz8rf/vY36d+/v4SEhMi8efOsUUbUw8rEY5KeWyjBfp5yeccIriUAwOHVOtwMGzbMvK9qatSQ8OzsbGncuLF5xBRsxxd/lHckHhvfXLzc/1zsFAAAR1XneW7OFhwcbImXgYWl5RTI8t0Zen/8RXQkBgA4B8YEO7AFCclSZhK5uGWwtAlnRmIAgHMg3DgoNbni4s3lI9muYUZiAIATIdw4qMT0k7IvI0883VxlaCc6EgMAnAfhxkF9tzVV3w5sHyaB3h5GFwcAAOcJN7NmzdLLOXh7e0uPHj1kzZo11Z67du1a6du3rx527uPjIx06dJDXX3+9QctrL01S354JN3/r0tTo4gAAYH+jpepKzYszZcoUHXBUaHn33Xdl+PDhenh5ixZVR/f4+fnJ5MmTpUuXLnpfhZ0777xT799xxx2G/Ay2aGdqrhzMPCVe7q4yOJYmKQCAc3ExVSzrbQA1w3F8fLzMnj3bfCw2NlZGjx6t16qqiauvvlqHm48//rhG5+fm5kpQUJDk5ORIYGCgOKIXf9gts1fulys6NZF/T+xhdHEAAKi32nx/G9YsVVRUJAkJCTJ06NBKx9X9devW1eg11CKd6tyBAwdaqZT2R2XViv42f+tKkxQAwPkY1iyVmZmpF+GMiKjcbKLup6WlXfC5kZGRcuzYMSkpKZFp06bJbbfdVu25hYWFejs7+TmybSk5kpSdLz4ebnJZh3CjiwMAgPN1KD53yQZV8/BXyzioTscbNmyQf//73zJz5kz5/PPPqz1XNW+paqyKLSrKsRePrKi1uSw2XHw9De1SBQCAIQz79gsNDRU3N7cqtTQZGRlVanPOpUZXKZ07d5b09HRde3Pttdee99xHHnlEpk6dWqnmxlEDTqVRUp1pkgIAOCfDam48PT310O+lS5dWOq7u9+nTp1Zf6Gc3O53Ly8tLdzw6e3NUm5NPSMqJ0+Lr6SaDaJICADgpQ9stVI3KxIkTpWfPntK7d2+ZM2eOJCUlyaRJk8y1LikpKTJ37lx9/5133tFDxNX8NooaCv7KK6/Ivffea+SPYTMqam2GxEaItwcrgAMAnJOh4Wb8+PGSlZUl06dPl9TUVImLi5MlS5ZIdHS0flwdU2GnQllZmQ48Bw8eFHd3d2ndurW88MILeq4bZ1dWZpIl25i4DwAAQ+e5MYKjznOTcPi4jJ29Tvy93GXD40OouQEAOBS7mOcGlrV0Z7q+VX1taJICADgzwo2DWLarPNxc3pHlFgAAzo1w4wDUOlL7MvLE3dVFBrYLM7o4AAAYinDjAH4+U2vTKyZYgnw8jC4OAACGItw4gJ/O9Le5nBXAAQAg3Ni746eKZMOhbL0/mHADAADhxt6tSMyQMpNIhyYBEhXsa3RxAAAwHM1Sdo5RUgAAVEa4sWOFJaWyKvGY3mcIOAAA5Qg3dmz9/iw5VVQqEYFeEtcsyOjiAABgEwg3DtAkpToSu7q6GF0cAABsAuHGTqklwZbtzND7NEkBAPAnwo2d2p6SK2m5BeLr6Sa9Y0KMLg4AADaDcGOnlp5pkhrQNoyFMgEAOAvhxs6XXBjCQpkAAFRCuLFD6bkFsuNorri4iAxqz0KZAACcjXBjh1Ymlnck7hrZSEL8vYwuDgAANoVwY4eW7y4PN4PahxtdFAAAbA7hxs4UlZTJ2r2Zev+yDoQbAADORbixM38cytazEof6e0mnZoFGFwcAAJtDuLEzK8xNUmHMSgwAwHkQbuzM8jOdiQfRJAUAwHkRbuzI4axTcuDYKXF3dZF+bUONLg4AADaJcGOHTVI9WzaWQG8Po4sDAIBNItzYkeWJx/Qto6QAAKge4cZO5BeVyK8HsvQ+89sAAFA9wo2dWLcvS89xE9nYR9qE+xtdHAAAbBbhxk6sqBgl1T5cXNSiUgAA4LwIN3bAZDKZOxPT3wYAgAsj3NiBPel5cjSnQLzcXeWSmBCjiwMAgE0j3NjRQpm9W4eIj6eb0cUBAMCmEW7sAE1SAADUHOHGxuXkF0tC0nG9zxBwAAD+GuHGxq3ae0xKy0zSNtxfooJ9jS4OAAA2j3Bj42iSAgCgdgg3NkzV2KzaU77kAquAAwBQM4QbG7blyAnJPlUkAd7u0iO6sdHFAQDALhBu7KBJakC7MPFw41cFAEBN8I1pB/PbXNY+3OiiAABgNwg3Nio9t0B2HM0VtYzUwPZhRhcHAAC7Qbix8SapLpGNJNTfy+jiAABgNwg3NoomKQAA6oZwY4MKS0pl7b5Mvc8q4AAA1A7hxgb9fjBb8otKJSzASzo1CzS6OAAA2BXCjQ1asfvMxH3tw8TV1cXo4gAAYFcINzbGZDLJ8t3pep8mKQAAao9wY2P2ZeTJoax88XRzlX5tGQIOAEBtEW5szE87y2tt+rQJEX8vd6OLAwCA3SHc2JilZ8LN5R0jjC4KAAB2iXBjY7MSb04+ofcvjyXcAABQF4QbG7JsV3mtTbeoRhIe6G10cQAAsEuGh5tZs2ZJq1atxNvbW3r06CFr1qyp9txFixbJ5ZdfLmFhYRIYGCi9e/eWH3/8URwFTVIAANh5uJk3b55MmTJFHnvsMdm0aZP0799fhg8fLklJSec9f/Xq1TrcLFmyRBISEmTQoEEycuRI/Vx7l1dYIuv2Zen9YZ1okgIAoK5cTGpiFYP06tVL4uPjZfbs2eZjsbGxMnr0aJkxY0aNXqNTp04yfvx4efLJJ2t0fm5urgQFBUlOTo6u/bEV321NlXs+2yitQv1k+QMDxUUtBw4AAGr9/W1YzU1RUZGufRk6dGil4+r+unXravQaZWVlcvLkSQkODq72nMLCQn1Bzt5s0dKdaeZRUgQbAADqzrBwk5mZKaWlpRIRUbkJRt1PSyv/ov8rr776qpw6dUquueaaas9RNUAq6VVsUVFRYmuKS8vMq4APZQg4AAD23aH43FoK1UpWk5qLzz//XKZNm6b77YSHh1d73iOPPKKrsCq25ORkscWFMnMLSiTEz1O6t2hsdHEAALBrhk2BGxoaKm5ublVqaTIyMqrU5pxLBZpbb71V5s+fL0OGDLnguV5eXnqzh1FSg2PDxY2FMgEAsM+aG09PTz30e+nSpZWOq/t9+vS5YI3NzTffLJ999pmMGDFC7J2qqfppR3nAG9qxidHFAQDA7hm6eNHUqVNl4sSJ0rNnTz1nzZw5c/Qw8EmTJpmblFJSUmTu3LnmYHPjjTfKG2+8IZdccom51sfHx0f3p7FHO47mytGcAvHxcJN+bUONLg4AAHbP0HCjhnBnZWXJ9OnTJTU1VeLi4vQcNtHR0fpxdezsOW/effddKSkpkXvuuUdvFW666Sb58MMPxR79eKbWpn/bUPH2cDO6OAAA2D1D57kxgi3Nc6Mu/aBXVsqhrHx5Y0I3uapbc0PLAwCArbKLeW4gsj0lVwcbbw9XGcJCmQAAWAThxkDfbD2qbwfHRoifl6EthAAAOAzCjUHKykzy7ZbycDOySzOjigEAgMMh3BhkY9JxPUrK38tdLm0fZlQxAABwOIQbg3xzptZmaKcIRkkBAGBBhBsDlJSWyXfbUvX+yK40SQEAYEmEGwP8djBbMvOKpJGvh/Rrw8R9AABYEuHGwCap4XFNxcONXwEAAJbEN2sDKyopk++3l89KPLJr04Z+ewAAHB7hpoGt3XdMck4XS1iAl/RqFdLQbw8AgMMj3DSwb7aUdyQe0bmpuLm6NPTbAwDg8Ag3DehUYYn8dGahTEZJAQBgHYSbBvTt1qNyqqhUWob4SnyLRg351gAAOA3CTQP6/PdkfTv+ohbi4kKTFAAA1kC4aSC703Jlc/IJcXd1kXE9IhvqbQEAcDqEmwbyxZlamyGxEXqkFAAAsA7CTQMoKC6VRRuP6P0JF0c1xFsCAOC0CDcN4PvtqZJbUCLNG/lI/7asAA4AgDURbhq0I3EUc9sAAGBlhBsr25eRJ78fzBY1X9/fe9KRGAAAayPcWNm8P5L07aD24dI0yMfabwcAgNMj3FhRYUmpLNyYovcnXNzC6T9sAAA0BMKNFS3dmS7Zp4okItBLBrWnIzEAAA2BcGMlJpNJ/rPmoN6/pmeUuLtxqQEAaAh841rJrweyZUvyCfF0d5Ube7e01tsAAIBzEG6sZPaq/fr2mp6RzEgMAEADItxYwfaUHFm955ie0+bOAa2t8RYAAKAahBsr1tr8rUtTiQr2tcZbAACAahBuLOxg5in5fluq3p80kFobAAAaGuHGwuasPiBlJjVpX5jENg209MsDAIC/QLixoIzcAlmYUL76992D2ljypQEAQA0RbizovbUHpai0THpGN5aLWgZb8qUBAEANEW4sJCe/WD759bDev+tS+toAAGAUd8Pe2cEkH8/X89lEurvpRTIBAIAxCDcWEtc8SH5+4FJJyy0QV1cXS70sAACoJZqlLEhN2te8kY8lXxIAANQS4QYAADgUwg0AAHAohBsAAOBQCDcAAMChEG4AAIBDIdwAAACHQrgBAAAOhXADAAAcCuEGAAA4FMINAABwKIQbAADgUAg3AADAoRBuAACAQ3EXJ2MymfRtbm6u0UUBAAA1VPG9XfE9fiFOF25Onjypb6OioowuCgAAqMP3eFBQ0AXPcTHVJAI5kLKyMjl69KgEBASIi4uLxVOlCk3JyckSGBho0dcG19oofK651o6Iz7X9XWsVV1Swadasmbi6XrhXjdPV3KgLEhkZadX3UL88wk3D4Fo3HK4119oR8bm2r2v9VzU2FehQDAAAHArhBgAAOBTCjQV5eXnJU089pW9hXVzrhsO15lo7Ij7Xjn2tna5DMQAAcGzU3AAAAIdCuAEAAA6FcAMAABwK4QYAADgUwo2FzJo1S1q1aiXe3t7So0cPWbNmjaVe2mnNmDFDLrroIj2bdHh4uIwePVoSExMrnaP6w0+bNk3PWOnj4yOXXnqp7Nixw7AyO9K1VzN4T5kyxXyMa205KSkpcsMNN0hISIj4+vpKt27dJCEhgWttYSUlJfL444/rv83q70NMTIxMnz5dz1Rfgc913a1evVpGjhyp//6qvxeLFy+u9HhNrm1hYaHce++9EhoaKn5+fjJq1Cg5cuRIPUr155ujnr744guTh4eH6T//+Y9p586dpn/+858mPz8/0+HDh7m29TBs2DDTBx98YNq+fbtp8+bNphEjRphatGhhysvLM5/zwgsvmAICAkwLFy40bdu2zTR+/HhT06ZNTbm5uVz7Ovr9999NLVu2NHXp0kV/lrnWlpWdnW2Kjo423XzzzabffvvNdPDgQdOyZctM+/bt41pb2LPPPmsKCQkxffvtt/o6z58/3+Tv72+aOXMm19oClixZYnrsscf0318VJ7766qtKj9fk7/OkSZNMzZs3Ny1dutS0ceNG06BBg0xdu3Y1lZSU1KtshBsLuPjii/Uv6GwdOnQwPfzww5Z4eZyRkZGh/wGtWrVK3y8rKzM1adJE/wOqUFBQYAoKCjL9+9//5rrVwcmTJ01t27bVf2gGDhxoDjdca8t56KGHTP369av2ca615aj/ILrlllsqHbv66qtNN9xwA9faws4NNzX5HJ84cUJXDKgKggopKSkmV1dX0w8//FCv8tAsVU9FRUW6Onno0KGVjqv769atq+/L4yw5OTn6Njg4WN8ePHhQ0tLSKl17NUnUwIEDufZ1dM8998iIESNkyJAhlY5zrS3n66+/lp49e8rf//533dzavXt3+c9//sO1toJ+/frJzz//LHv27NH3t2zZImvXrpUrr7xS3+dzbT01ubbqu7O4uLjSOaoJKy4urt5/w51u4UxLy8zMlNLSUomIiKh0XN1Xv1hYhvoPg6lTp+o/VuqDr1Rc3/Nd+8OHD3Ppa+mLL76QjRs3yh9//FHlMa615Rw4cEBmz56tP8+PPvqo/P7773LffffpP/w33ngj19qCHnroIf0fRR06dBA3Nzf9t/q5556Ta6+9Vj/O59p6anJt1Tmenp7SuHFji39/Em4sRHWmOvfL+NxjqLvJkyfL1q1b9X91ce0tLzk5Wf75z3/KTz/9pDvFV4fPef2pzqyq5ub555/X91XNjepkqQKPCjdca8uZN2+efPLJJ/LZZ59Jp06dZPPmzbqTvKoduOmmm7jWDaAufzMs8f1Js1Q9qR7e6r8Izk2ZGRkZVRIr6kb1pFdV+StWrJDIyEjz8SZNmuhbrn39qeph9ZlVI/3c3d31tmrVKnnzzTf1fsVnmWtdf02bNpWOHTtWOhYbGytJSUl6n8+15fzf//2fPPzwwzJhwgTp3LmzTJw4Ue6//349GpBrbV01+Ryrc1TXjuPHj1d7Tl0RbupJVampL4SlS5dWOq7u9+nTp74v79RUelc1NosWLZLly5fr4ZxnU/fVP46zr736h6K+lLn2tTN48GDZtm2b/i/bik3VLlx//fV6Xw2h5VpbRt++fatMaaD6hERHR+t9PteWk5+fL66ulb/m1H+MVgwF51pbT02urfru9PDwqHROamqqbN++vf5/w+vVHRmVhoK/9957eij4lClT9FDwQ4cOcYXq4a677tI961euXGlKTU01b/n5+eZzVE98dc6iRYv0UMNrr72WoeAWcvZoKa61ZYfau7u7m5577jnT3r17TZ9++qnJ19fX9Mknn3CtLeymm27Sw4wrhoKrvxOhoaGmBx98kGttodGVmzZt0puKE6+99prer5gGpSZ/n9VI48jISD0dghoKftlllzEU3Ja88847eu4KT09PU3x8vHm4MupO/WM536bmvjl7uOFTTz2lhxx6eXmZBgwYoP8RwfLhhmttOd98840pLi5Of2bVtBFz5syp9DjX2jLUl6j6DKv5sby9vU0xMTF6XpbCwkKutQWsWLHivH+jVais6ef49OnTpsmTJ5uCg4NNPj4+pr/97W+mpKSkepfNRf1f/ep+AAAAbAd9bgAAgEMh3AAAAIdCuAEAAA6FcAMAABwK4QYAADgUwg0AAHAohBsAAOBQCDcArOrSSy/VixXaErUo3+LFi40uBgArYRI/AFaVnZ2t148JCAiQli1b6qDTUGFn2rRpOsSo9bHOphbza9y4sXh5eTVIOQA0LPcGfj8ATiY4ONjir6kW4FOL1tZ3xWIAjolmKQAN0iylbg8fPiz333+/bhZSW4V169bJgAEDxMfHR6KiouS+++6TU6dOmR9XNT7PPvus3HzzzRIUFCS33367Pv7QQw9Ju3btxNfXV69c/sQTT0hxcbF+7MMPP5Snn35atmzZYn4/dex8zVJqRfTLLrtMv39ISIjccccdkpeXZ35cve/o0aPllVdekaZNm+pz7rnnHvN7AbAthBsADWLRokUSGRkp06dPl9TUVL1VBIthw4bJ1VdfLVu3bpV58+bJ2rVrZfLkyZWe//LLL0tcXJwkJCToEKOopi4VWHbu3ClvvPGG/Oc//5HXX39dPzZ+/Hh54IEHpFOnTub3U8fOlZ+fL1dccYVupvrjjz9k/vz5smzZsirvv2LFCtm/f7++/eijj/T7VoQlALaFZikADdY85ebmpgPJ2c1CKrRcd9115n44bdu2lTfffFMGDhwos2fPFm9vb31c1az861//qvSajz/+eKXaHRVmVDh68MEHdS2Mv7+/uLu7X7AZ6tNPP5XTp0/L3Llzxc/PTx97++23ZeTIkfLiiy9KRESEPqbCjzqufoYOHTrIiBEj5OeffzbXIgGwHYQbAIZSNTH79u3TIaOCyWSSsrIyOXjwoMTGxupjPXv2rPLcBQsWyMyZM/XzVTNSSUmJBAYG1ur9d+3aJV27djUHG6Vv3776/RMTE83hRtUAqWBTQTVPqVonALaHcAPAUCpE3HnnnbqfzblatGhh3j87fCi//vqrTJgwQferUc1aqi/OF198Ia+++mqt3l8FqbP7/5zt7ONqxNe5j6myA7A9hBsADUaNcCotLa10LD4+Xnbs2CFt2rSp1Wv98ssvEh0dLY899pj5mOqw/Ffvd66OHTvqPjSqA3NFgFKv7erqqjsrA7A/dCgG0GBUv5jVq1dLSkqKZGZmmkc8rV+/Xo8+UvPR7N27V77++mu59957L/haKgwlJSXp2hrV0Vf10/nqq6+qvJ9q2lKvq96vsLCwyutcf/31ul/PTTfdJNu3b9cdhtV7T5w40dwkBcC+EG4ANBg1UurQoUPSunVrCQsL08e6dOkiq1at0qGmf//+0r17dz0aSvVpuZCrrrpKDytXo5q6deumh5NXjKKqMHbsWD0SatCgQfr9Pv/88yqvo4aR//jjj3qywYsuukjGjRsngwcP1p2HAdgnZigGAAAOhZobAADgUAg3AADAoRBuAACAQyHcAAAAh0K4AQAADoVwAwAAHArhBgAAOBTCDQAAcCiEGwAA4FAINwAAwKEQbgAAgEMh3AAAAHEk/w991KjJch28uAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "latent = jnp.zeros(Ndes) + 0.5\n", "optimizer = optax.adam(learning_rate=0.1)\n", "opt_state = optimizer.init(latent)\n", "\n", "history = []\n", "n_iters = 100\n", "for i in range(n_iters):\n", " neg_abs, grad = value_and_grad(latent)\n", " updates, opt_state = optimizer.update(grad, opt_state)\n", " latent = optax.apply_updates(latent, updates)\n", " history.append(-float(neg_abs))\n", " if i % 25 == 0 or i == n_iters - 1:\n", " print(f\"iter {i:4d} absorption = {history[-1]:.4f}\")\n", "\n", "plt.plot(history)\n", "plt.xlabel(\"iteration\")\n", "plt.ylabel(\"absorbed power\")\n", "plt.title(\"loss curve\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "10", "metadata": {}, "source": [ "## 5. Optimized structure and field" ] }, { "cell_type": "code", "execution_count": null, "id": "11", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAywAAAFfCAYAAACsp0jbAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZiZJREFUeJzt3QmcFMX5P/6nZ/bk2FVArnCICoKiqKCIF6CCQcQrRjy+iFd+EBBFvECTiEhETOSPSsQQ78SDb6Io8YsKiXIokgBCRDQeAWFBFgLCLix7zXT9X0/hrNNdtTvdPcf29H7eeXWwa7t7e3pmu6a66qnHEEIIAgAAAAAA8KFQY58AAAAAAABAfdBgAQAAAAAA30KDBQAAAAAAfAsNFgAAAAAA8C00WAAAAAAAwLfQYAEAAAAAAN9CgwUAAAAAAHwLDRYAAAAAAPAtNFgAAAAAAMC30GABAAAAAADfQoMFAAAAAAB8Cw0WAAAAyGrV1dV0ww03UOfOnamoqIhOP/10WrlyZWOfFgCkCBosAAAAkNUikQh169aNPvzwQ9q3bx/9/Oc/p4svvpgOHjzY2KcGAClgCCFEKg4EAAAA4BetWrWi999/n/r06dPYpwIASUIPSxP18ssv0+zZs7U/MwyDpk6dmtLfl45jJrJ06VL5e/lfL/h8ef/G9vzzz8vz+Oabbxy9fwAAfnHkkUcq9/5FixbJe1p9yzPPPOPpuPH+/e9/U2VlJR199NGOznPatGl03HHHkWmadWU8pIx/B/fYeJHs/vXd/+szf/58Ov7446mwsFDus379elf7e/3durqS38Mf/ehHVFFR4fr3AuigwdJENfSF96OPPqKbb745pb8vHcdsKoYPHy6vX4cOHerK0GABgGz18ccfy3/ffPNNeW+zLz/96U+TOj4PAxs1ahT94he/oBYtWiTc/ttvv6VHHnlENlpCoZClwfHAAw8k1WBJZn83/vvf/8rXzA20d955R17HHj16aOuPTBg9ejQ1b95cXleAVMhJyVEgUDhYMRuO2VQcccQRcgEACEqDhQPjR4wYkfJe7NraWrryyitlb8m9997raJ/HHnuMDjvsMLr88sspW3355Zfytf/P//wPDRw4sK68WbNmjVJ/5OTk0JgxY+jBBx+ke+65R54HQDLQw5IlPvjgAzrvvPOoZcuW8g//jDPOoP/7v//Tdt8uWbJEzpbC43f5CQdXCps2barbbtCgQXLfLVu2WLrh6xu+FTvue++9Rz/72c+odevWsrK57rrrZHdvaWmprCD4hs9Pce68805544xnPyZ359c3HCB+CNdXX31F11xzDbVt25by8/OpV69e9Lvf/U7b/f/jH/9YXps2bdrQ2LFjaf/+/Y6vL1+Pk046Sf4ODtz87W9/q92OQ76efPJJuS13ux9++OF0xRVXWK5v7Br37t2bVq9eTWeffbY8r6OOOooefvhhy5AD/u/p06fTscceK4/H1/DEE0+UFWh93fL1vX98bt27d6cLLrhAOe8DBw5QcXExjR8/3vE1AQBIh7Vr18q4klQ3Vvh+yvVSOByWQ5KcHL+mpkZuy/VMfO8K11d33XWX/G+uE+z1U6I6OdH+X3/9tayn+Z7N+/PwKa6rN2zY4Pp1X3/99XTWWWfJ/x45cqT8PVxPNDSsy2nd6rWuZNdeey2Vl5fTq6++6vo1AdihhyULLFu2jIYMGSK/yPKNlW8U/KWZb26vvPKKvEHFu+mmm+T2PGyopKREdovzzeuTTz6RX4h53//3//4f/ec//6EFCxY4Pg8e0sVPoPjms27dOvn0imdm+eKLL2Q5H/Nvf/sbzZw5kzp27EiTJk2q91j8e3kayviKhhsZ/MW/S5cusuyzzz6TlQCvP/roo9S+fXt699136dZbb6Xdu3fT/fffL7fbuXOnfKKUm5srX1u7du3opZdeoltuucXR6/r73/9Ol1xyCQ0YMEC+tmg0Krux+bh2/MSIKwA+B36d3333nRxGwOf5r3/9S/7uGG7I8Q37jjvukOfKr3nKlCny2nClyvj3cMXG79E555wjG3rc+GpoCEF97x9XShMmTKCJEyfKyogrwpgXX3xRVhxosABAY9qzZw9t3bqVLrroIll/2HFjw2tDhu/PO3bskEOi+Am/E//4xz/kOQ0ePFip7/j+/sQTT9Drr79eN6SKe26c1MkN7R8bhsYP//ghFveA8LYvvPAC9e/fX9av/BDLqV/+8pd02mmnyfv7Qw89JF8LP1Ssj9O6NZm6kvFxe/bsKRs5N954o+PXA6DFs4SBv51++umibdu2Yv/+/XVlkUhE9O7dW3Tq1EmYpinLnnvuOZ7xTVx22WWW/T/88ENZPn369Lqy4cOHi65du2p/H297//33163HjjthwgTLdpdeeqksnzVrlqX8pJNOEqecckqDx7S75ZZbRE5Ojli0aFFd2QUXXCBfX1lZmbJtQUGB+O677+T6PffcIwzDEOvXr7dsN2TIEPl733//fdGQ/v37i44dO4rKysq6svLyctGqVSu5f8xHH30k1x999FHL/iUlJaKwsFDcfffddWUDBw6U2/7jH/+wbHvcccfJ1xVz0UUXyevVkNj137x5c8L3j8+7ZcuW4rbbblN+7+DBgxv8PQAAqcb3qfh7/+LFi+X9rL7l008/9XTcb775Ru7PdUPz5s3rluXLlzd4nJkzZ8r9SktLlZ/95je/Ue69burk+vbX4f1rampE9+7dxe23397g/V+H6zne7s9//rOlXLe/07rVvq/TujLetddeK9q1a5fw9QMkgiFhPsdDrvgJEA87ig8e5KdQHGC3bds22cMRj5/qx+MnKV27dpXTOyaDn4jF4y5kxkF99nIeruQUP2GaM2cOPfXUUzRs2DBZVlVVJZ/mXHbZZbK7nJ/ExZYLL7xQ/nzVqlVyW35dPDOKfepK7u52cn152Bb3EBUUFNSVczc/Py2L99Zbb8knfzxGOP58+CkS/277bGRczk+94vETufhrwz/nnplx48bJJ1zcC5IMPm8eZsC9QLHZWXgoHz9Rc9rjBO7wZ5HfNy8L7wvQ1IaDMe514HuvfYn1QLjFdRw/G+OZwXgIbGzhIbkN4Z4Ovq/zUOJ01ck6XHdwbwi/3ry8PNkjxP9y7/jnn39O6eKmbvVaV8bjIWe7du3S9qZBelUFrG7CkDCf27t3r7wJ62b44KFFjLuz7V+U7bjMvp1bHBMTj2+u9ZU7/bD/6U9/kkPLfvWrX8mhbDF8rnyD4+50XnS46zq2LY+ltdNdB9315eFo9V2zeNztze9F/LCveByjEo+7++146ABXqDE8RIzjjPg6cIONKz0eGsbDzfr160de8LAwbgDysDgeOsb/3alTJ9mVD6nFn/Nu3X5EpaXfedqfP2ObN2+2fAEACHrAPX/eOakj3+8aG9+PeTix03PxUifr8JBpjhnhgHQe0szxkBxDw0PJ4uuIVHNTt3qtK+Pxe83Xi++VTmZsg1TWTe2ptLQsMHUTGiw+F7uJ8bhc3ZMhZn8yxLETdlx2zDHHkJ/w5AA8rpUDBnnqR/vrjj2xqi/uItZI4YZBfa85Ef49/HTNyf58nXnbFStWyIaHna4sEX6qxhUXLxy3wjFA3IDjwHmOP/Iyswq/z9xTxZUh/7tw4UJ5ff3w5SBoOGCXGyvfbPlfKipy916Vlx+kI7teKY/hp0oBIN0NFp6QxMn9iB+ycG9A/FN+vv/GAsxTge/r/DfIx+aHR+mok3X4IRXHMnIvi72xwLGm6eKmbvVaV8bj2ByuG9FYaYy6qYy+KXmMiooKXe1bXl5JR3a+zXd1ExosPsc3UA7C4+5zno2DZ5Ji/KSDb3j85JznWo/HT9Z/8pOfWOaC52FI8XlQ7E/6M40TWvE5nnvuuTRv3jzl5/xFnQMHOfiQh1HFenN0eDsO/OOhVfHDwnjSASfXl4dl8fX9zW9+U/fHyTOM/fWvf1WGxPHwte3bt8tZ0VKNKykeZsDH58B5ntWlvuERid6/2267jYYOHSrnwufKiWd3g/QpalFARS3cVQoUN1scQFNQVlYmJ1ZxGoDNeVri72k83IoDvlOJg8IZT2LCdY3uIVT8vdZNnazbP4a//NsfcnFwOt//0/lw0U3d6rWujMfvt9dhfpC8Fi24sejuYWr8TKZ+ggZLFpgxY4ackYRvMjxlMN9geEaSTz/9VM5IYp9RZc2aNbJxwsm3+Cn9fffdJ6dM5DiJmBNOOEHeeObOnUt9+/aVT4y8DkFyi8dH8lhZvtHz6+Hzjcc3N57hhKf25SdpPAb55z//uZwKmW+OPB0k3yA5NoPxl/tnn31WxtLwFMGxWcJ4ti0neJ54nhKZrzHP6MUzn/CQLL5B89OhmDPPPFMOseIYET5nHrrF2/CTNp7ikq8pn6cbPPaXnzbyteeZYrhhyQk9eTx2/CxfdoneP34tfB05vodjbngcMaQR3+Dd3uR9WikApLN3hYcH8X1TFyfB9VTnzp2Vch46y7Edb7zxRsp7imPT//L52BssfJ9lXBfxwx8eOsazdzmtk+vbn+M++AEYxxpyg4l/L8f2cEOAGzzp5rRu9VpXxn/x/ec//2kZ7g2ZJURELm738SM0WLIAj2/lGwhPNcjDp/gmwD0JPNTHHgjPeJrFP/7xj3TVVVfJqYP5pso3qPhYE35atXHjRjn8iJ96cSVyaDKv9OMbW6w7nW98dvwlmysR/sLNFRzfJHnaXw7c414I/iLPDZ74sZY8zSS/Jr758hMkDijk2A0ncRt8DlwR8u/g6Sj5eNy446di9qFqv//972USTP6XKyh+L3jcMjdm7AH2TvB789prr9HTTz8tG3L8u/l8eJpKrtzq4+T9414gnjIZwfYZgAYLgOMM948//rhc7Pg+aP9yy9PGc5D7okWLHPcGuMENJP7izr05/EAqHtdD3Fji6Yb/8Ic/yPt9rH5yUic3tD/XyXyP58YPTw5wyimnyIdQXA+lm9O6NZm6kvFENFw/2ScCgswRIioXt/v4kcFThTX2SUBq8NMafvrPM3lkqrcE/Is/A/ykjz8PkB7cyOSEnN+VvkZFRc1d7ltBrdr/RFboDeVMAMhm/PSev9THJw52atasWfSXv/yFFi9erMRAJHNcO35oxF/AuYebe3kgNThOhoeEffjhh7ikjVQ3le5+1FMMS/s2d/iubkIPC0DAblI8LIGnYOYhBm4Sg0ISTOFhSBieFQHUh4e78tBe7slId8A2T9V76qmnyt4O7pmH5HFM0Pz58+sdXgaZITAkDAD8iLv5eZgZz5zGwxUuvfTSxj4lAADXeCgVT80a3+PBw6m4YZFq3BPNQ7Z4SBcP2+KYQEjO1q1bZeMvlTO6QdOGHpYA4e5xXqDp4rHRGOXZCBDDApBSPM17JvHkJ7xAavCDM17ADzEsEdf7+BEaLAAAyUKDBUCLp2fPpuMCBIkwI3Jxu48focECAJAsNFgAAMBvROTQ4nYfH0rbQE2e8pWzpXJyIc4TwdlpAQACSXyfh8XNwvtARqFeAoCmGHQvXC5NpoeFZ4bgZH5cOXB+Cs5ZMWzYMPrss8+oS5cuDe7LAW/ffvutTKxkT4gIAJAKHOfDidI4h04qAmwNYcrF7T6QOcnUSwx1EwBkW91EPLzLrHW/T1PJw9K/f3+ZBImnJYzp1auXnLGIpw1syLZt27SZbgEAUq2kpCSpzNKxue73bnqRilo2c7fv/oN0+FHX+W6u+6BKpl5iqJsAINvqpm3b7qaionyX+1ZTp06P+K5uSnkPS01Njcz/MHnyZEv50KFDaeXKlQn3556VQ7h3BT0sAJAO/JxGxN1vkoQYFl9Ltl5iqJsAINvqJiGHeIVd7+NHKW+w7N69m6LRKLVr185SzuulpaXK9tXV1XKJ4a6wQ9BgAYB0EqkbdioTR7rsrEbiyIxxWy8x1E0AkP11U4TIDAdiSFjagu7tF5tHnuneAO6K526r2ILhYACQddwG3HvpkYGM1UsMdRMAZD0z4m1pCg2WNm3aUDgcVp5a7dq1S3m6Fctmy+PkYguP2wMAyCposPia23qJoW4CgOwX/WFqY6cL79MUGix5eXlyGuMlS5ZYynn9jDPOULbPz8+XQT3xCwBAVuG5S3jWL1dLyuc7gRTVSwx1EwBkO8OMeFqazLTGkyZNolGjRlG/fv1owIABNG/ePNq6dSuNHTs2Hb8OAAAA9RIAQEClpcEycuRI2rNnD02bNo127NhBvXv3pkWLFlHXrl3T8esAABoXZgnzPdRLANDkmByT4nIwVVPqYWHjxo2TCwBA4GGWsKyAegkAmhQTDRYAAKirFDzM+oVZwgAAII0MESFDhFzv06R6WAAAmgwOonfbAOF9AAAA0vowLep+Hx9CgwUAIEmGacrF7T4AAADpYshZv9wloTSaWuJIAAAAAACAZKGHBQAgJXlYXOZVQR4WAABIJzPqYZYwfyaORIMFACBZCLoHAABfzhJmuN/Hh9BgAQBIFhosAADgM4YZJcNlDwvv40dosAAAJAt5WAAAwG+EhyFhvI8PocECAJAs9LAAAIAvZ7CMut7HjzBLGAAAAAAA+BYaLAAAKRkSZrpcXM4qBgAA4HqWMA+LC8uXL6cRI0ZQx44dyTAMeuONN+rddsyYMXKb2bNnk1tosAAApGpaY7cLAABAWoPuo64XNyoqKqhPnz40Z86cBrfjhsw//vEP2bDxAjEsAADJQgwLAAD4jRn1MK2xuwbLsGHD5NKQ7du30y233ELvvvsuDR8+nLxAgwUAIFncW+J2iBd6WAAAII0M2WNiNOq0xqZp0qhRo+iuu+6i448/3vNx0GABAEj+jnxocbsPNGEun3pmZCR3Kj+TIouvcabfi3TfC4L2XmQ6msHM3vfG9N7DUl5ebinOz8+Xi1szZ86knJwcuvXWWykZiGEBAAAAAIA6nTt3puLi4rplxowZ5NbatWvpscceo+eff14G2ycDPSwAAMlCDwsAAPiMYQrXeVV4H1ZSUkJFRUV15V56V1asWEG7du2iLl261JVFo1G644475Exh33zzjeNjocECAJAsZLoHAAC/MXlImId9iGRjJb7B4gXHrpx//vmWsgsuuECW33DDDa6OhQYLAECyhHlocbsPAABAuggPDRbex4UDBw7Q119/Xbe+efNmWr9+PbVq1Ur2rLRu3dqyfW5uLrVv356OPfZYV78HDRYAgGShh6UJyUt6LLbTEFLDyFXLNPsZRmrCUYXDRrSwfwPS7CdI96XHTF9D3vE1sG5nUNjRsXTXXdnG4Tk4uc7KNZaFZmquseZYjjl6jbrrF07JNc6uz7tOw79TyNkjqyhVDGGSIVzOEubys7FmzRoaPHhw3fqkSZPkv6NHj5axK6mCBgsAQLIQwwIAAAEaEubUoEGDvm9oOeMmbiUeZgkDAAAAAADfQg8LAECyMCQMAAB82ftvuN/Hh9BgAQBISYPFbeJIPyRzAwCAwDLRYAEAAEuDxWUDBA2WrBQOtUhJ0K/9GLrg71BIfaYYMtSysC0432nwshO64G/TFpRrilpn+5kRzXZRB4HQXp/4Jp6gIKSZ2EB3jXXvuf09Czn8XNivny5gW3cdTBFJfN0zfP2cfm7tn0n9dfceiO+E/TNpfx/qey/s1113jbVlDt9X+8+jZgqD7k2TDJdvv9u8LZmCHhYAgKR5mNbY85cIAAAApz0sHvbxITRYAACShR4WAADwGzM4DRbMEgYAAAAAAL6FHhYAgGShhwUAAPzGDE4PCxosAADJQuLIJiMnXJiaoHsHQchhTfC3frvcBtd1v88pbfC8g0D5qCYg3AypZVEngeSac3CSpVwfKJ++6+410732Wgkn1yono9dKdx1Sea1CmgD+dH5u7cH09b0XQkQbfB/q4/SzrAbdU+qIqPsJYVzHY2YGGiwAAMlCDwsAAPiMgVnCAACgDhosAADgN2ZwhoQh6B4AAAAAAHwLQ8IAAJKFGJYmg5MKpiaGJZwwFsAwNEn5NGVqMj9v8QG6Y+uYtjH9QpOZLmroYjDUspCtTBeXoYsZEIaThJPkIHFkKGXXXXeNdUyDEu6ni3UI2cqE7pmzcBif4iABpJP4FF28lLaMEic39fz508arqOy/UbefNmbG9n7pzlIfm5LrKJlk4uMkIUA9LGiwAAAkS4hDi9t9AAAA0jpc2XS/jw9hSBgAQKpiWNwuHjz55JPUrVs3KigooL59+9KKFSsa3P6ll16iPn36ULNmzahDhw50ww030J49ezy+UAAAyBpm5uom3zVYli9fTiNGjKCOHTuSYRj0xhtvWH4uhKCpU6fKnxcWFtKgQYNo48aNqTxnAIAmWSnMnz+fJk6cSPfddx+tW7eOzj77bBo2bBht3bpVu/0HH3xA1113Hd10003yPvznP/+ZVq9eTTfffDMFCeolAIAGhiu7XYLQYKmoqJBP6+bMmaP9+SOPPEKzZs2SP+eKsX379jRkyBDav39/Ks4XAMB/hIcKwcNc93xv5cYHNzh69epFs2fPps6dO9PcuXO1269atYqOPPJIuvXWW2WvzFlnnUVjxoyhNWvWUJCgXgIACHaDxXUMCz/N40WHe1e4AuWnf5dffrkse+GFF6hdu3b08ssvy4oSAAB+UF5ebrkc+fn5crGrqamhtWvX0uTJky3lQ4cOpZUrV2ov6RlnnCHvx4sWLZL37V27dtFf/vIXGj58eKDegiDUS/oA4HDCgOND24USBsHrQ4W9UYOjwykL/o44/K4UFQ0HRstTcPBQwNRsY+gCozWB5Pbj24Pp6+PkvOzJCuXxPSb00yeFDHtKAJkTym8wmF7upw3gV7fL9N+TfdIHbZJIbeLNxEH+kGUxLJs3b6bS0lJZgcZwxTtw4MB6K1QAgKY8JIx7SIqLi+uWGTNmaH/F7t27KRqNyi/a8Xid77v1NVg4hmXkyJGUl5cne7wPO+wweuKJJ6ipQL0EAE2WGZwYlpTOEharNHUV6pYtW7T7VFdXy6W+p40AAL7HD9zc3uS/f0hXUlJCRUVFdcW63pV4HDto70Gwl8V89tlncjjYr371K7rgggtox44ddNddd9HYsWPpmWeeoabAS73EUDcBQNYTPPzYcLlPE2iweKlQ+WniAw88kI7TAADwfaZ7bqzEN1jq06ZNGwqHw0pvCg/zsn8Zj7+/nnnmmbKRwk488URq3ry5DNafPn26nDWsqXBTLzHUTQCQ9QTXTR72CfqQMB5uwNxUqFOmTKGysrK6hZ82AgBkE2EKT4sbPKSLpzFesmSJpZzXeeiXzsGDBykUst7mudEjz9mnlVKqeamXGOomAMh6JoaEafEsNFw5cAV68skn1wWKLlu2jGbOnKndp74AUwCArJGhxJGTJk2iUaNGUb9+/WjAgAE0b948OaUxD/GKfcnevn07vfjii3Kdp6D/2c9+JmcRiw0J42mRTzvtNDn1fFPgpV7KRN1kz3gddRhkrcvMbg8C1maU12Vqtz+z1HwktVnEnWQa15yDNsjZHriuyWqvO5YTukkLQiHrwJKwNqu9JtjcUD8LOZTvKbA8StbXE6EfhsXHmIb6eYjaro32WpkR19nV66O7DiFbcL4uoF93HezHcprpXpt53sNnVLedbmID/bG8Bdl7ve6pb7CQ+32CMCTswIED9PXXX1sCGtevX0+tWrWiLl26yMrwoYceou7du8uF/5sTll1zzTWpPncAgCaFg+c56eO0adNk46N3795yBrCuXbvKn3NZfE6W66+/Xk4pz9PM33HHHTLg/txzz23wi3o2Qr0EABBsrhssPH//4MGDLU/82OjRo+n555+nu+++myorK2ncuHG0d+9e6t+/Py1evJhatmyZ2jMHAAhADItbfG/lRYfvwXYTJkyQS5ChXgIAqCfm3mUPi8cOJf81WDhzfUNjnzmIkTPd8wIA0CRksMECKtRLAAAaTXlIGAAA2KDBAi7ZYzf0yeh08Ry65Iu1CWM3vI651yU5tMdO6LfRxavo4its+2qO5SQWQB+vkqduaDt8OJzvKAbDHq/C8oxm1m2Es69UEfvxNd8Po7b3lNVSZcJ4FVNo4mHMmpRdU7LFrOhiTAwHSSh1sS/2+Jj6tkvl51vZRvv5S/yN30msTX2vJ/1T7pP7fXwIDRYAgGShwQIAAH5jBqfBkuGmHgBA8PAwWdfTGjeRaYUBAKCRCI+LC8uXL5czUvLMkxwW8sYbb9T9rLa2lu655x464YQTZA4w3ua6666jb7/91vVLQYMFAAAAAABcq6iooD59+sjZKHW5wD7++GP65S9/Kf99/fXX6csvv6SLL77Y9e/BkDAAgGRhSBgAAPiMMA25uNvH3e8YNmyYXHSKi4uVZMdPPPGEzAXGU/BzOhSn0GABAEgWGizgkprEztm3BF0AsBrwG3F2LCVpY+LEjofKIg6C6TXJHrUBzd6GRhpk/RImNN/JtNfU3Xe3BpMa2oPsc0kT5K9je8kRzbE9XpZ6gsY1web291DzC/XvYcPvw6FCdfCOSfage12AvbOgeycB/F7pP+8Ogu51kwNoJwzInqD78vLytCTTLSsrk0PHOC+YGxgSBgCQqgaL2wUAACBdhEFkuly+b/137txZ9pDElhkzZiR9OlVVVTR58mSZTL6oqMjVvuhhAQBIFgfQuw2iR9A9AAD4dEhYSUmJpVGRbO8KB+BfddVVZJomPfnkk673R4MFACBJQcomDAAAAWF+32viap9D/3BjxW0vSEONlSuvvJI2b95M7733nqfjosECAAAAAAApF2usfPXVV/T+++9T69atPR0HDRYAgGQh6B58zkkwsZMAe+2EAQ4yiNdHG7StbKQLt7VlXNdkp3caCO2EKdTXaNqug3293mPZ99Mc2wnngd7qtbEPYDV0557CbmD1M6J5HzSfNd3XVEeZ5zOdUd6vxA8xKc73cbf5gQMH6Ouvv65b516U9evXU6tWrWTelSuuuEJOafzWW29RNBql0tJSuR3/PC/P4UQVaLAAAKQAGiwAANAEpzVes2YNDR48uG590qRJ8t/Ro0fT1KlTaeHChXL9pJNOsuzHvS2DBg1y/HvQwwIAkCTEsAAAgO+YIQ8xLO66WLjRIRqYRKahn7mBBgsAQLKEh2mKMUsYAAD4NOjeb9BgAQBoxORc0DQpCeR0+ffIaYyCJvGgjTA0MRi2mAFDaBL+acqUxJGa89TGGTh4jU6TB9pfs9N4jpDt+CHD2blHNUkUo4atzOEzC/t+UV3CTs052M81HMpz9JodxS9p4pAcJQ11eCNT32fd+5WTMEmkbl/7e+r0b8Lr34nuuvgiSaSGEIZcgvAsrfGvJgAAAAAAQD3QwwIAkCRhCrm43QcAACCbY1gyBQ0WAIBkYUgYAAD4ckIYw+U+aLAAAAQT39/d3uP9WScAAEBQCA9B927ztmQIelgAAJKEIWHgli6YWN0m5CwI2UE4qi44WtgSFkY1iftMkThxny4IWcdJEHyOka9skxMqSLifju687AHu2gB7TRC8E2FNgkYd+/FrRaWjc7C/92HNtQo7TKCZKLCcRcwqtUxUJ9zPyefBaZC67vMetr33hhFOWcC77vNgCgeTJDgMxM+OoHuD/AgNFgCAZGFIGAAA+DKGJeRyH/Klxm/+AQAAAAAA1AM9LAAASUKmewAA8BthGh6C7jEkDAAgmDAkDAAAfEYghgUAAKxTR7qsSHw6Thgahz2QuL5gc11QupMAY12mdnugcEibcT0/YRCyjtMA6rxQM+u6YV1nBaK5UpYv1GtjV22oQeMHjXLLepWwrrNa86BSFjGr1TLb8XXB3zpOJjvQTbiQa7tWBUaRsk0zKkrZtarKqVDKaoT12tRorpXu86ELZvc6oYR9YgFdcLvTCRDsnEzCYJ94oL7XbGpu8oKs733amcGJYcGQMACAZAkPN3lMawwAAGkkMCQMAADqKgVxaHFVkaDBAgAAaSQCNCQMs4QBAAAAAIBvYUgYAECSEMMCAAC+YyKGBQAA6ioFDzEsPg1shMywB1XrAo61Wd81QfAhW7C3aQvqltuQGhAeNRIHzzvJuK6jC5bOCannbg+ybykOV7Y5zNQElxt5lvUwqcNYqkz1OuwxrAHo/zXUAOpKc69SFjHVbPR2huYa6zgJvM4JFSplhWHrtTlctFW2aS0OU8oKNJMBRG1BdAfNGmWbfSHrNWb7bZfZDKmvJaK5tzkJuncSYK8r0wXYOz2Wk/PUTVhhp06bwH+H6mcr0wRiWAAA4IdKAbOEAQCAv4gAxbBgSBgAQJIQdA8AAL4jPAwJ8+mEMGiwAAAkizMDu80O7NNswgAAEAwCQ8IAAAAgnfTJF53FSdiZmrgJe5kuXkU3pt9r4khThBOeQ45Qn6O21MTytMm3blcQ1sSwRNVjhauscSDlITVepUyTyDES3U+JGA5iJJhwcP10cUL5RgvLeltTjffpUKDGneivjfUx+u5qdZsDmvfC/n6ZmmuVysSRus+kPVGkNl7FwSS4Ie02mjLbpdHFzOhjvdRjRTXxZeAMelgAAJKEGBYAAPDncGXD9T5+hAYLAECSghTYCAAAAWEacliY2338yFUkzowZM+jUU0+lli1bUtu2benSSy+lL774wrKNEIKmTp1KHTt2pMLCQho0aBBt3Lgx1ecNAOC7Hha3C6QG6iYAAJUQIU+LH7k6q2XLltH48eNp1apVtGTJEopEIjR06FCqqKio2+aRRx6hWbNm0Zw5c2j16tXUvn17GjJkCO3fn3j8JwBA1na7u22w+LTbPRuhbgIAaGBCGLdLtg8Je+eddyzrzz33nOxpWbt2LZ1zzjmyd2X27Nl033330eWXXy63eeGFF6hdu3b08ssv05gxY1J79pBSkVnXK2Xl66zJpIpPtyb5koacqhSJdu2s6wXWZF3MiKjBeoXF6jlEIvusx3I4555hi5TLyVETalWWPa+UiRzrn4VRVaUee+dO9RcuWa0Ula06aFkvOlkNiMyZpJ4DZBcMCWtc2Vg3mbYutpDhMGGdg9ufbj9dcHTEtN7bopptoprkd8J27vb1+oRDavLFaNh6rgUha2A5Kwy1Ucra5Fvv763z1WDmqGbYZa4tgH+XLQifleqCxoWmHrB9hQppkj3qRG1JKAU5C1wvsgXZdylUJyPo2kJ9L8Kaz9aeausECBURdUKECvOApmyPZb0mqm4T1SShTBQ4X19ZWDPhgv3a5IQKHCVFzRUFrhNJHvp92RsoLwI0XDmpfp+ysjL5b6tWreS/mzdvptLSUtnrEpOfn08DBw6klStXJnuuAAAAqJsAAJoYz0H3/MRq0qRJdNZZZ1Hv3r1lGTdWGD+1isfrW7Zs0R6nurpaLjHl5eVeTwkAoHEEKLAx26FuAgAIXh4Wzz0st9xyC33yySf0yiuvKD8zDEOpQOxl8cGSxcXFdUvnzp29nhIAQKNmune7QOqhbgIAyFzQ/fLly2nEiBFysi3+rv/GG2+kZTIuTw2WCRMm0MKFC+n999+nTp061ZVzgH18T0vMrl27lF6XmClTpsihZbGlpKTEyykBADT6OGG3C6QW6iYAALWHxe3iBk+81adPHznZlk6qJuNyNSSMW0lcISxYsICWLl1K3bp1s/yc1/lEeAaxk08+WZbV1NTIGVxmzpypPSbHuPACPvCzQ8Go8Yp3Whufxueb1P227VDLQrYP/GFqYCPVqIF5kag6JLDmkWut53DjReSEePYty3rePWpvIH0fh2U5fp41MN7Yp2ZCNnSvuYsaGFp8wVHWc2p3qFFvgaD7rBekbvdsFIS6SWgCeyOmGvBuGroAbXv2cTXgOKIJnrcHR0c1v09oAvHtZcKW/bw+kaga2F1ryyAfyVHPITdHnaykvOwIy3pRjhpAbWq6MXfZft8O42vNear1gpLunOuU3LaW9dZ5R5MTe2r+Y1mvrt3h6Bx2kPVc/13ZXNnmuxp1YpyQZoRLecR6nb8N/VfZZmfkS6WsOvKdZd3UBNg7+TwYpH4WDEP9Shox1IkawqH8hJ/tHE2wfsSwbpdjO059WewFmQmz2tu30U2sEdSg+2HDhslFf6zUTXjiqsHCUxrzL3jzzTdlLpZYTwoP5eJuHu4KmjhxIj300EPUvXt3ufB/N2vWjK655ho3vwoAIGt4GeKFIWGpg7oJACC1DRZ7TLmXhziJJuNKW4Nl7ty58l8ef2afQvL66w9NR3v33XdTZWUljRs3jvbu3Uv9+/enxYsXywYOAABAqqFuAgBILXtM+f333y9jUdzwMhlXyoaEJcK9LPyC3L4oAIBsFaS57rMR6iYAgHrqJtNb3cQx5UVFRXXlyQyRdTMZV8qnNYYAipteOqb9KX+wrFdU73J0KPsY1tqILrmUNami3O7A6+qxcq3jlx2PvJlojX2pHf9TZZPcFmrcTjhkHQOcm9PC0fhbneb51jHOO76+09F+kF1M05CL232g6bLfIyOaGBZDEz+iY49Z0SXu08W1mMK6ndBsQ7oYFocxK+p+KjNqS15pVijb/McWN8E22eIdjFpdTIQ6r5A9Tkf3+3TXISesJh5um9/Tst41eiQ5Ec63xknssCWSZJGoNWEy21v5hWV9bWibemxNXIYusaf9PdTFKumvjf29N1P2WTCEpl7VxLXYP7dRU41xihhqos9wyLpdVKjXKqRJJmnY5qfSx6vo4loaP+Gk8DDrV+wBEDdW4hssXsRPxtWhQwdHk3GlJXEkAABgWmMAAGias4Q1JH7Ck5jYhCdnnHEGuYEGCwBAFk1r/OSTT8pKoKCggPr27UsrVqxocHtOzMsztHTt2lV26R999NH07LPPenylAACQLUQG6qYDBw7Q+vXr5RILtOf/3rp1q2UyLp7F8dNPP5Ux714m48KQMACALIlhmT9/vrz5c6PlzDPPpN///vdyOsnPPvuMunTpot3nyiuvpJ07d9IzzzxDxxxzjOyKj0ScDTUCAIDsJTJQN61Zs4YGDx5ctz5p0iT57+jRo+n5559P2WRcaLAAAGQJTr5100030c033yzXeX77d999V86SNWPGDGX7d955R3a9b9q0iVq1aiXLjjzS2Th7AACARHjm4IYmPknVZFxosECd3COsQeq6ADtDkzyLNIGN9tGGuoRQhXkd1N8XcjBK0ck2zDQTHlt3DlW1e6zrNbt1B1eLNIGNB6utQZG57UY1dMaQpUxhyMXtPm7wuN+1a9fS5MmTLeU8vz3PZ6+zcOFC6tevn8w0/Mc//pGaN29OF198MT344IMydxY0Hl0gtJ1JmiSRZiRhwK8ugNrUJNezB9QLfVg8eaGdmERTD4RD1s9hSJPwzx4srQt8jkTVAPFotNLTazY0gdc5YTVJY8hWz9Vq3i8d+366YzuZDCCqSS4ZNRNfY93vNEjdRvde2D9HUc2EAV4natBuo/sibDt+1NAEvGvP3bqdqTnPkOYzqgvE92OAvQ7fZtwnNSZfQoMFAKARM907Tc61e/duikaj2vnsY3Pd23HPygcffCDjXXj8MB+Du+W/++47xLEAAAScCNCU+wi6BwBIUaZ7t0ssOVdxcXHdohva5XU+e9M05c9eeuklOu200+jCCy+Uw8p4XDGPKQYAgOAS309r7HbxI/SwAAAkySQPQ8LIXXKuNm3aUDgcVnpTGprPnue9/9GPfiQbQjG9evWSjZxt27ZR9+7dXZ0zAABkDzMDw5UzxZ/NKACAJjJ1ZCw5V2ypr8GSl5cnpzGOn8+e8Xp989nzTGLffvutnHYy5ssvv6RQKESdOnVK6TUAAACfMT3kYPFpUmP0sECd/Nw2ytUoKuhsWc+1ZYFnNaYmi70tuC2sCVr7UegEpSy0a5dSZn6fKfWHAm+RY7pj9yg8Tynbnr+hwddS3+vJC7VQymrNg5b18qoSZZuqmm+VMgAdni5y1KhRMpB+wIABNG/ePDnX/dixY+XPp0yZQtu3b6cXX3xRrvM89xxgf8MNN9ADDzwgY1juuusuuvHGGxF07zPa7NkOAuzlvvZAck0wcWqD5+2TqugC5QuUssKcQzPVxSvKsU580iFqrXPYkXnq9KcFYeuXqm+rrNnP2efGv5Wy/1ZbyyKawHVd1HFNRN2uVHxmWd+X4+xeXlVtzWJfG93v6BwMw/pe5IR/6DmNOSK/p1LWS6hlHQusExlURdXg9m9q1PPaEbbWYeWRHco2lZHvlLKoac08L7STQKivWT8RhJ1mP+3fgP2cNKdgaI4VMhvMfF8fQzshEXiFBgsAQJK4t8TMQGDjyJEjac+ePTRt2jTasWMH9e7dmxYtWiSTQjIu4wZMTIsWLWQPzIQJE2Qjp3Xr1jIvy/Tp013/bgAAyC4iQEH3aLAAAGRRpcCzfPGiw8H0dj179lSGkQEAQPAJNFgAACDG9DDgxqdT3QMAQEAINFggiL4bd45Sln/zWQ0mY2TGl98oZR/Nto47ffprNTHWmqotSlmb42cpZTkh6/joCe2GkRNP7Hzbsh4x1TGznXL6KGUXNbvAsn7zMWoCrwET1afjoocmg7gtWWX10x8omzSb/aq6H2SVIFUKEACasfOGUOPuhGFPDKw5lCZez55cT5f4sCB8mFLWJnRo6GK8o8wfWdZPP0KNmbnulDVKWdvu1vrj83+epGzz+40nKmWvm3st62UR6zoTmgSQwlTjOaptZdW16Y1HNIT1ujfPVWcH/HGh+prHHP+1UtbrtPWW9V1fqe/Nix/3VcpW7e5hWd+kee93G2rdXhXdlzDRpy6RI2lithzFWTliOkoAaf/aEwqpg5OcxrVkmilCcnG7jx/586wAAAAAAAAQwwIAkDxTuJ+7nvcBAABIa++/GYzefwTdAwAkCUPCAADAb0SAhiujwQIAkJIeFvf7AAAApItAgwWCKG/yJUqZKMhP/C2ruEgpKgivtax/W2lNGsXWP67u9+YTFyplN/77Lcv6jK3PkhP5udagz2d7XqRsc8kENZDywlut51oQ1gRgnqoGI4qWanIzChkJrzEh6D7rBalSgMxQksoJZ4nnQpSbMHTY0M1B5/HjZtgC7FnYNhFKXkgNvC4MqUkNWwj1nt8qz3r8o1uqiYi7Xr5OKYseb50wpfNONSlvzmfHJAz+1gXYhwz19YTDzRwk7HQ295/9fdVd42jUmnSYmaKiwdfCcjTvc+eu6rXJu8qaoLPrRvUaH/2VmnDyy3JrguRdNep7WqF573UJT+2imolxhG2igeTYr3vIUQB/NieAND3kCHO7faaghwUAIEmm/IroslLw+g0SAACgiT1My95mIwAAAAAABB56WAAAkiTEocXtPgAAAOkiAtTDggYLAECSgjROGAAAgsEMUN2EBgv8oLkaaKikeK0oV7dZ/JFS9Kt/WYPu3vrv+eqhNdd+xFVq6e6q023noGbIdfJ6REGBeg6aYLq3rrauX9r2b8o2CzWvmX58plpWVJT4GkPWEx5iWHgfgIYD7NUqWhjqPdJwkJnaaUC4cmxtYHLi3xclNcj6oKEG1O+psU5W8mW5eo/c8vrJSlnoTevr+d+Pz1a2ebvyc6WsqmaHZT0vp72yTZd8dVKVIlMNJK82rEHilYYaKK9TKKwB/PnCNrkNEZXnlyllW6vXNvha6nvNR61Qr82VFdZjmVH1Guveiz01kYTvqe69d/IZChnqhBLCYxZ7J4HyTj/b2Rx0L2Tvv9seFvIlNFgAAJKEIWEAAOA3AkPCAAAgiN3uAAAQnAaLGZAYluzt5wIAAAAAgMDDkDAAgCRxPIrbmBTEsAAAQDoJDAmDQKqsVIpCW7Za1o8a+r6yzday5ep+oTzLulF9lrKNyFHby4Y9yF9uaCsLOewYtO1nRDQZ6zXHsm/37oGnlW1yx9QoZV2K1WzCmxYPtqybXbskOGnIRqY4tLjdB6DhDNvOAo6FEU0YYK8L1ncSiK/bxrRlh68x1YlQIkLNWl4dUgO094e+s6x/u0cNgv/kb6coZbtsmeA/ib6tbFNZs10py8s5wrLeJ3eIss2ROdaJANjBqHodtpvWwPhyYy85USQOs6x3CKvZ4o8OW7dhh5utLOv/oiXKNiUHlillD9R8rZQtKB1kWW8btk4EwDaZu5Sy74xSy/pBU33NtaY6+YApIok/o5qydAbP6/7mQh4D7E3d6yHr32VjEGiwAABADGJYAADAb8wAxVdiSBgAQJIwJAwAAPxGoIcFAABiMCQMAAD8xkQPCwRRQYcxSllUMxbVLhxSx772ajHcsm7YYmGkNq3VsqgaZxLaZk2OJdZ8mfCc5O/s18OybnbqoG4T1sTR7N7T4Gthn+//q1K2ZZ+aYDJ82t8SXisAAB2n4+nt8VBCm5aXEsYMmEKX8E8zNt8Wa6j9fZox/VWao5fTJst6qeZYn2njA6wvOhxqoWzTrvmpSlkv83jLevcWakLhPdXq79tIm5Wy7TUfW9YjEWcxLKU5h1vW9xWoMTqniKOUspNaWGNdmh/Q1E35G5Wy3Zpkkv+oecmybugmDdF+/hLHgWgp8SMOEzSKxAkmnSZ2tJ9r2MhxFDdmP1f931etozhFP8S1pFokEqGpU6fSSy+9RKWlpdShQwe6/vrr6Re/+AWFnMYcO4AhYQAAScKQMAAAaIp108yZM+mpp56iF154gY4//nhas2YN3XDDDVRcXEy33XYbpQoaLAAAScKQMAAAaIoxLB999BFdcsklNHz4oR6/I488kl555RXZcEklV301c+fOpRNPPJGKiorkMmDAAHr77R+mEhRCyG6hjh07UmFhIQ0aNIg2blS7JwEAgjhO2O0CqYG6CQAgtXVTeXm5ZamuVqcqZ2eddRb9/e9/py+/PDRc/1//+hd98MEHdOGFF1KjNVg6depEDz/8sGw18XLuuefKVlWsUfLII4/QrFmzaM6cObR69Wpq3749DRkyhPbv35/SkwYA8BPhcYHUQN0EAFB/D4vbhXXu3FkO64otM2bM0PwGonvuuYeuvvpq6tmzJ+Xm5tLJJ59MEydOlGWpZAjuFklCq1at6De/+Q3deOONsmeFT5JPnnFrrF27dnJ825gxakC3Drfi+MIcakvhCWQmGZoRguGwNZCxRWEnZZtrD7tMKZt2tjVAsXkXNdDsw+U/Usoe/VxtQ1fZAkF/3EFN6qXzzg5rQ7nAFqjH7uilBs+deY412VjFVjUI71cruillL+1boJQdqNxmWY9G1cRpwpaEDTKBb3smlZWVyd5ir2L3qz/3vYOahfNd7XswWk0/Xfto0ucAma2bCvO6OQ7ydSOZJHZR2z3SnqRPbmOqT0eFbTvdfkIXiG9PAqgNJHb61cJazxuG+ndUmKcmkzw6bA2o71vYVtnm+GL13JvlWM91wz416P6tg+rIkG0HViQM0G6Zr9YLOvurrfWjEOr169TibKXsombWCQNOOEydxuBgRP0cbSxT6761ldakkP+Jrla2qayxJok8dK7VqXmfdcH62iB49dxDtu0MzX7hUH7C/XI0nzV7QL+OLuheN2GFl2SSPBFGZc3mlNVNC/peSs3DiV9TvIpoLV229g0qKSmxnEN+fr5c7F599VW666675P2WY1jWr18v77fcgTF69GhKFc933Wg0Kk+yoqJCDg3bvHmznB1g6NChddvwCxs4cCCtXLkyVecLAACAugkAII1i4R+xRddYYdxYmTx5Ml111VV0wgkn0KhRo+j222+vt0cmY0H3GzZskA2UqqoqatGiBS1YsICOO+64ukYJP7WKx+tbtmyp93j8pCt+XBy3CgEAsgk/RzM97AOpg7oJACDzQfcHDx5Upi8Oh8Nk2qY/z3iD5dhjj5XdPfv27aPXXntNdvcsW7as7ueGYX2hPOLMXhaPW2APPPCA29MAAPCNIGUTzlaomwAArEwy5OKG2+1HjBhBv/71r6lLly5ySNi6devkcDAejptKroeE5eXl0THHHEP9+vWTjY0+ffrQY489JgPsGQ8Li7dr1y6l1yXelClT5Fi92MJj5gAAsrGHxe0CqYO6CQDAxkvAvcuHaU888QRdccUVNG7cOOrVqxfdeeedMjbwwQcfpFRKOg8L96DwkK5u3brJRsuSJUvkDAGspqZG9r5wYGN96gvigcw7s3icUrbsCdvXqlN6qjvmaDLDfmste3vyQWWbUZ8vVsp27n5YPb498NRh5tTb7d2RmgC4dm0mK2V/jP4Qh8WGPaxmp3/8frUR/nhkrHoSH//bsjpwgnruH5Q9ru4HWQV5WPwn2+omXTCu00zZUbPGdYD9oePbttMF3ScVUG9lGHlKWU6YJ9n5QbsCa2A5G5ynlo3rvcmyftKQdxydwz/fHmRZ/+v2jso22w+qcbd5OUcoZaMP/4llfVSP+oe/x/vjl+db1l/Y+5p6DhXqOXwjjrWsX32MNXCenTZsqaNzWL/kLMv6k58eyqER7/1cdfKBnVXWski0TNlGCOvn8fvSxJPNaCYf0H0mo7bg+ZAmeJ6D1xMG4oecfSl2EojvV6aHKfTdbt+yZUuaPXu2XNLJVYPl3nvvpWHDhsmpzniqYg66X7p0Kb3zzjty2BfPCvDQQw9R9+7d5cL/3axZM7rmmmvS9woAABoZMt03LtRNAADBHq7sqsGyc+dOGf2/Y8cOOV0aJ5HkxgrnWmF33303VVZWym6hvXv3Uv/+/Wnx4sWy9QUAAJAOqJsAAILNVYPlmWeeafDn3MvCme55AQBoKjAkrHGhbgIACPYMlknHsEBwPHu6moSKTuptXW+uxnNQ3LTUMc9PtCZRunaENYEi27lCE6/iRAqnytPFzJjTXrSsPz9RTQY2+n81B9Ndm5N6WFafPf1TZZMe7zo7V/AvDAmDpD9DmjH3pmacv2nqYlHsiSNrHCT80yUsND3Gq6hDSHQJIHPDhyllh+V3saz3Mq33THZBRzXdwcnDfpidlIX6tVG2qf6bWqd9uMMas/IvUu/JukSO5+ZfopQ9eNHfLettnplCTvS4yZqfouR/1WO/U/G0UmY/1w932OpnDjM9UKiU5Z+vJsc8Ocd6/S74Tr1+pZvU96Im3xqPuq96q7JNbXSfx4STwlkMlS3fuan53BqkiT9z8tVBE9cStu1oaDbykiQyE0RTHRIGAAAq9LAAAIA/6ybD9T5+hAYLAECS0GABAAC/EWTIxe0+foQGCwBAkoJUKQAAQDCYGZjWOFNcJ44EAAAAAADIFPSwQJ0Xv26vXI3xU9da1qf90xogyV7Zt0Apyw9bp7Ie9YtfKtsYNZrkUlFNMqkaa0ApRTTb6OTYPt55muRPYfVPIPSL/7GsT5mjZmud1G+/Unb1YZcpZb86zRqQ+OLXXRs8ZchOHAPqdtyvLW4UmniQvdBEBOsD7NXgedOWOFJQraNAcjUKWfehNBI+6zQ0ifXCoeZKWX6OmuKg0CiyHktzCtsq1AlN/vnWuZb1itfUIOsVpWqC37e+22NZ31n1ibJNbk4rpeyoFmpdUVB0wLJeWbWdnLDvpzt2brV6Djsrref66s5OyjYHXrpCKTv77zuVsub51QmvMfcbJ3q/KjXvqX0SCBa153HWbON80gfTwb202kFgvCZQXhc8HzITJpL0Q4B90Icro8ECAJCkIE0dCQAAwSACNFwZDRYAgCQFaepIAAAIBjNAMSxosAAAJAk9LAAA4DdCuB9+7NfhymiwAAAkKUjjhAEAIBgEGWRiSBgEzdO731PKZr+9w7JeoQsqFGpgaItm1uzwxt69yjahrzcpZf0u3qKUfV2zwrIeiarZi3Vywtbsvsfkna1ss2ahGgRvHnOUZf1g7W5lmwMHNytlT1X8Rin749s/sqy3yPtPgrMGgKbIa4C9LsjeWYC9TlgpMXRfdgzrs86QJqu9LjBZp0pYA9BLjP8q27xb2lopW77Lem/da6r1wqbQx0rZnuovLOumac3czvJzDlfK1pVbz5M99JeLLesn/O3f5MSGfRcnPHZOSM1YXxv5zrL+WfViZZude49Vyt7fo5YdHrJmtq801c9Macg6QYHu/dLRBqXbPiPaT6Pmu4RwEHSvD57X7Gb7KJumw7lzzYaD8A8dWt3RMDARbyqhhwUAIElcN7rtMEEHCwAApJMIUHwlGiwAACkZEuY2sBGXHQAA0sdE0D0AAMSghwUAAPxGBKj3Hz0sUGd/tRqfUlG1NfF4Zo21Zw+2rEf+v7eUbVo9vlIpq4mqCRlJk8jJiZpImWX9k5rXlW1aDFGTXn136xkNvhZ27LubHY2/PVD5jXUTj68F/A1B9+CWPVGkLvGcNomdLkGdo/uKZoy9/X6uGXNvkBqPYNhiWOzr9ak1K5WyA8Ias3LQUOMmtmmOVVNTYT227X7PTJE43lGX9LI2aj02+zSk1lf/2WeNdckva0FOVNviQCrFXkfnYD/XSLRc2WZnxWql7L/GBqUsN6fYsp4XVhN96j5q9uSLUU2clY79M6KL7hCaGCrSJEElW2+2Ps5Fs5s9tkvzdcZRXItmm1BI8zegOa1Mx7WY6GEBAIC6SgGJIwEAwGfMANVNmMIAAAAAAAB8Cw0WAIAUJedyu3jx5JNPUrdu3aigoID69u1LK1ZYp/2uz4cffkg5OTl00kknefvFAACQlbOECZeLH6HBAgCQouRcbhbex6358+fTxIkT6b777qN169bR2WefTcOGDaOtW62xZnZlZWV03XXX0XnnnZfEqwQAgGxifh/D4nbxIwTdQ52DmqD7UMia7Km6Ug1c1wVq2j3eR30KXFX7pqMg03C4mWX9xBaXkxOfHLCeazSqJgir0iQb+91bx1jWb/2XmnCyhq5Qf6Hm3PMLL094jSH7eekx8dLDMmvWLLrpppvo5ptvluuzZ8+md999l+bOnUszZsyod78xY8bQNddcQ+FwmN544w33vxh8EZjvdIS5Ptmj7j4dShiA7jUhninUYGkzqgZo1whrsLzQ7Ce0gd3WPyBD83UmN6dVwoSMEc1EABHN5C/7q/apZbqodE80CTs170VOuGXC5JLa1xNRg/qra7+1raufGUOTENR+XiHthAuJPx/6z5CzY9kTpfIjIHUjB5EYukktDPU9NZUJdTR/E8KfiSNFgGYJa/yrCQAQkMBGt4sbNTU1tHbtWho6dKilnNdXrlRnMIp57rnn6D//+Q/df//9Hl8dAABkIxM9LAAA8EOl4D4RZGz78nLr1KT5+flysdu9ezdFo1Fq166dpZzXS0tLtb/jq6++osmTJ8s4F45fAQCApsPELGEAAJAKnTt3puLi4rqloaFdzDBsOQiEUMoYN254GNgDDzxAPXr0wJsFAABZC4/cAAAacZxwSUkJFRUV1ZXreldYmzZtZAyKvTdl165dSq8L279/P61Zs0YG599yyy2yzDRN2cDh3pbFixfTueee6/KsAQAgWwgPs375dZYwNFigTusWvZWrUVo6zbIuwt4+Mret6a+U3dm8RimrqVmk7hzyGmp1VcIUtnl5FyY8V+H496vb1ez/X8t6+/a/UrbZvf9jh8eHIA4J48ZKfIOlPnl5eXIa4yVLltBll11WV87rl1xyibI9H3PDhg3KlMjvvfce/eUvf5FTI0M20tyPbMG9hi5juDaLvT3o3tu9VihByfoJA3TB80r2cW0gu/oFKjfnCMt6+4ITlG3aRTsqZRUhawb5LZF1yja1kd2aM1CvabP8oy3rbXKd/U3trt1sWT9YvcXRtcoNt7esd805WdmmualmrN9ZYA2wZ6VVGxK+ZiHUSWnskyIII+wsWN/BZ01bpg26D3manEIony2Hk1rYg/MN3Wdbs58PvvcLD0PC/Bp0jwYLAECWzBI2adIkGjVqFPXr148GDBhA8+bNk1Majx07Vv58ypQptH37dnrxxRcpFApR797WhxBt27aV+Vvs5QAAEDyCPPSweGhpcb1zzz330Ntvv02VlZVyGPIzzzwjH7KlChosAABZEtg4cuRI2rNnD02bNo127NghGx6LFi2irl27yp9zWaKcLAAA0DSYSfT+O7V3714688wzafDgwbLBwg/GeGbKww47jFIJDRYAgCyoFGLGjRsnF53nn3++wX2nTp0qFwAACD6RgTwsM2fOlJPH8BT6MUceeSSlGvKwAAAAAABAHZ5yP36prtYlcCVauHChHKb805/+VPaunHzyyfSHP/yBUg09LFBnX+U36gdx4AeW9RErNbMKOQhKF/kFSlntvv91FrjmZeyM9iRMZ+egOVdHNEH99uunu8aQ/YKUTRj8QxdwTJoM6E720wfUpybIXh90r8tYX+vgL0EdPx8OqZNStLIFvJ8guivbHF6gfsXZUpVnWf9GmxFd/essKlSnBh+aN8iyfnob9TrorNptPdbi0FJlm7LKz9Szsp1ra/NwZZuuBYVKWdsaNRC/Jv+gZX23qU6CEzWtOaK+L23wnOpn/9yq741hOPtKav8s6zLd6wPxbZ93x+fuje4ctH/TGUgc6UZse+41icfJh3U99Js2baK5c+fKGMt7772X/vnPf9Ktt94qZ7y87rrrKFXQYAEAyKIhYQAAAOmOr3Q65T5Pl889LA899JBc5x6WjRs3ykYMGiwAAH6biYXSPxMLAABAJvKwOJ1yv0OHDnTcccdZynr16kWvvfYapRJ6WAAAUjHXvdtpjXHVAQAgy2ewPPPMM+mLL76wlH355Zd1s1emChosUCcatY5pZVd8PMeyXrujp7KN2aFj4rgWTXyH8ZImSeS1FyaOKXGayNH2O43qKkfnIH52ZeLfp3k9oR3fJrx+urHekP0wJAzSwWkiPUpZAkinCfhMB/EqUU/NdEMTo5MTVmMw8oxmCQ9dGVXPfU9or2W9prpMcxZqMsTWIfXLV48i6y89sc1/yYnvaqwJINd+px67jKxfAOW5RqznuifX+lpY26h+2E6i66e7xqaoVMqEsMe6CGfvvfKQX/fZdhB3ov1867Zx8Pl2mgAyi4kM5Ai7/fbb6YwzzpBDwq688koZw8I5wnhJJcwSBgAAAAAArp166qm0YMECeuWVV2RusAcffJBmz55N1157LaUSelgAAJKEWcIAAMBvTNlv5XKWMA/xlRdddJFc0impHpYZM2aQYRg0ceLEujIhhJz2rGPHjlRYWEiDBg2SswUAAAR9SJjbBVIP9RIAQPDqJs8NltWrV8vxaSeeeKKl/JFHHqFZs2bRnDlz5Dbt27enIUOG0P79+1NxvgAAviM8/g9SC/USAECc72NY3Cx+rZo8DQk7cOCAHJvGmSynT59u6V3hcWv33XcfXX755bLshRdeoHbt2tHLL79MY8aMSd2ZQ8oJUZUw+Crc9fr0XvnxT6f3+NlyDpBVEHTf+LK9XjI0gd5OR2Zo93VA2JIA6pi6pJC2BJD6APvExz7Eeu6GkmBQP4lApbAGoP+btirb1JpqZu49kU0JkyPqrmcFqQHuG/cdaT2W6EBO/LtMJDy2oXnz7ee6JbJO2eaAcZRSlkv5Ca+ffoKHXAefGWfvvRrIrU7UoA25N6yJPnV075fu9aiTRejOPVhB92aGhoT5todl/PjxNHz4cDr//PMt5Zs3b6bS0lIaOnSoJdHMwIEDaeXKlcmfLQCADwmPC6QO6iUAACu3vSteZhXzbQ/Lq6++Sh9//LHserfjxgrjJ1fxeH3Lli3a41VXV8slprxcfdoBAACQqXoJdRMAQBb3sJSUlNBtt91Gf/rTn6igwJYbIw4H4sfjLnl7WXyAZHFxcd3SuXNnN6cEANDoghTYmG3SUS8x1E0AEJTEkabLJesbLGvXrqVdu3ZR3759KScnRy7Lli2jxx9/XP537AlW7IlWDO9jf7oVM2XKFCorK6tbuPIBAMgmQep2zzbpqJcY6iYAyHZmgB6muRoSdt5559GGDRssZTfccAP17NmT7rnnHjrqqKPkrGBLliyhk08+Wf68pqZGVh4zZ87UHpNjXHgBAMhWXp5K+fUpVrZJR73kl7rJaTCxvUyXrdtJBm/hIMD+UGHE46fZSBhcbhg5js69OmodPr7HULOy10QPqPvV7rEfXXOa6jWuNtVjbbeXlbUgJ7ZHrPtV0wFH52C/7lU1/1U22a25Vnlh9bxMZeIEdT/de2EIW0+lNkBb943XdnzdZ017JN3nPXEgvk4oZH89zj5r9jJdsL7u79IPhId4SZ+2V9w1WFq2bCmzWMZr3rw5tW7duq6cc7I89NBD1L17d7nwfzdr1oyuueaa1J45AIBPYJawxoN6CQCgobrJ5SxhPm2xpDzT/d13302VlZU0btw42rt3L/Xv358WL14sKxUAAIBMQ70EAJDdkm6wLF261LLOQYyc6Z4XAIAmwUtMik+fYgUB6iUAAG/xkqKp9LAAADQ1iGGBZDkdA5/KsfL2mBVdDIsuPkUdw5/ENxwHr8ceb8FqTWvMihAVyjaRqFpmahIkq9RzigpNEsrwbtsvdHBo3i9k3S8arfY0J5LutdREyjTbRRJ+jnTX2NH7pU0aqiMSx4FoXrM+rsq6nWE4jPVycE21CScNs8Hf72dmgOIr0WABAEhSkJ5iAQBAMIgA1U1osAAAJClIT7EAACAYzADVTWiwAAAkiZMQ8uJ2HwAAgHQRHvKq+LVqyp6BeAAAAAAA0OSghwUAIEnIwwLpkMoAe11QtT0YWxcILbRB1V4fwepeT+LXqEvmF4lWJgzONoWaTJI0r9GJiKkGxleE91oLHL5dFbQ34bGdiTp6zRHNS9YlhXTG/iJDHq+xcPZZ0/wNqJMIaALsNYHx9r8nx39f9lN1mNbESZB/uommmjgSAACCXSkAAEAwmB6GhDWZxJEAAE1NkCoFAAAIBoFZwgAAIAYNFgAA8BsTs4QBAIB1SJjLWcJw+QAAII3MAPX+Y0gYAABAgOizlmuyK9iD2R1mundGjUw2HEQrC+15Jk4hb2oy0esnDFDPygldtvhq84C1wGGMtX0/3bGdnauzwHWT1Gvj5FS174VyRur1E9pr6uRbsNP33pbZ3paJvr6/ASXoXhesbz/2oULbOWlO07fZS4IDDRYAgCQF6SkWAAAEgwjQhDBosAAAJClIgY0AABAMZoAepqHBAgCQJI5fMV3HsPi0VgAAgEAQAXqYhgYLAECSglQpgH/oEibqkt2ZZiThfro4CUG1DhrRjf9B1ce12Mt0MQSpiyvQxUTUmrrElInZ99PHHHnl7DooiTZTmKTUO01MjvbzZ71epiZJZEgTn2L/u9B9rgwjnDCuJZviVUzMEgYAAEGsFAAAIEB1k3C/jx/5oUkNAAAAAABZbMaMGWQYBk2cODHlx8aQMACAJAkh5OJ2HwAAgCDMErZ69WqaN28enXjiiZQO6GEBAEjRTCxuFwAAgHQRHuolL8/SDhw4QNdeey394Q9/oMMPPzwdLwU9LAAAyTI9zBLmdnsAXYA9ExRtOKBaltU6SDKYuiSR2uehTgK7NRMGaDezvWan+3mlu6ZRe7JKh6dg30937JTSXBthWK+f4fV2pHtPtd947efg9BfqJgywl2g+25rPnz04XzeBhfK50gTda/dL8+cvqQlhPOzj1vjx42n48OF0/vnn0/Tp0ykdMCQMACAV3e5uZwnDVQcAAJ9OCFNeXm4pz8/Pl4vdq6++Sh9//LEcEpZOGBIGAJCiHha3CwAAQLqYcpiXcLkc2rdz585UXFxct3BAvV1JSQnddttt9Kc//YkKCgrS+kaihwUAAAAAACyNkaKiorp1Xe/K2rVradeuXdS3b9+6smg0SsuXL6c5c+ZQdXU1hcNqbhsv0GABAMiSccIAAACZmCWMGyvxDRad8847jzZs2GApu+GGG6hnz550zz33pKyxwtBgAQBIEoLuIVlOg3Z1QcH2oG1tELc2sNtrIHRihjYQ3xvda7afuz4jutfXowv0NhNPgOBwkL19P/17b6Y1W7yhXD9KGGyezHvv/ZPlIIBfN8mE5tyFLejeNDVZ7UO6ySK8/q2q2xkZjsQweUiYh32catmyJfXu3dtS1rx5c2rdurVSniw0WAAAkiTH/bqdJQxdLAAAkEbi+/+53ceP0GABAEhSkCoFAAAIBjPNPSw6S5cupXRAgwUAIEl8f3dbKaC5AgAAfp3W2G8wrTEAAAAAAPgWelgAAJKEoHtwy0ngrpMAe/n5swfdazN/pzIo3fD47DPkMcBet6H9+qU5073m+KawXmdhRj2997pjp5bm+ML2HjrM+u7sPTUdbGcm8XkUDj7btQkz3evO3VC24a2y96uyEB6GK/s0vjJ73wUAAJ8IUqUAAADBYAZoSBgaLAAASUIPCwAA+I0I0MM0xLAAAKSoweJ28eLJJ5+kbt26UUFBgcwuvGLFinq3ff3112nIkCF0xBFHyARgAwYMoHfffTeJVwoAANk2IYzpYvFncwUNFgCApHlrrrjveJ8/fz5NnDiR7rvvPlq3bh2dffbZNGzYMNq6dat2++XLl8sGy6JFi2jt2rU0ePBgGjFihNwXModjFuyLsg1F1UWzn36ptSwyvsO+aL+apIah+Z+MibAvnqXv3LV010+zcDyRdXH6fln3c/r70iuF11jz3iufj3Sfu6PrqS7a9ysN9/KM5ggT7hc/Qg8LAECWmDVrFt1000108803U69evWj27NnUuXNnmjt3rnZ7/vndd99Np556KnXv3p0eeugh+e9f//rXjJ87AACAV2iwAABkwZCwmpoa2UsydOhQSzmvr1y50tl5mibt37+fWrVq5ep3AwBA9iY1Fi7/l/UNlqlTp5JhGJalffv2lkAd3qZjx45UWFhIgwYNoo0bN6bjvAEAAtFgKS8vtyzV1dXa37F7926KRqPUrl07Szmvl5aWOjrPRx99lCoqKujKK6+kIEHdBACgchu/koEBl5nrYTn++ONpx44ddcuGDRvqfvbII4/IIQtz5syh1atXy8YMj5/mJ3oAAEFlevwf4yFdxcXFdcuMGTMa/F38oCgePyiyl+m88sor8os9x8G0bduWggZ1EwBA400I47tpjXNyciy9KvGVJo+X5mDQyy+/XJa98MIL8unfyy+/TGPGjEnNGQMA+IwwBAnD3XOpWLd7SUmJnMErJj8/X7t9mzZtKBwOK70pu3btUnpd7LiRwrEvf/7zn+n888+nIApm3WQ6Shxp3y69SSL9MZo808NW9EkUbddB+97ojmV6S5aZxuuX6jD41CWcdPI+q9s4SSapJpKs57xsp2U4nEBCt519wg2nx/JKBtG7/FsJTND9V199JYd88bSaV111FW3atEmWb968WVak8eOrueIdOHCg4/HVAADZSHh4ghWrULmxEr/U12DJy8uT0xgvWbLEUs7rZ5xxRoM9K9dff738cj58+HAKKtRNAADBjWFx1cPSv39/evHFF6lHjx60c+dOmj59uqwoOU4l9tRPN756y5Yt9R6Tx2vHj9nmMdwAAKCaNGkSjRo1ivr16ydzqsybN09OaTx27Fj58ylTptD27dvlfTrWWLnuuuvoscceo9NPP73uPs0xhjz8LChQNwEABJurBgvP9x9zwgknyArz6KOPlt3rXBl6GV/N47UfeOAB92cOAOATHI9iuAxVjMWwuDFy5Ejas2cPTZs2TcYQ9u7dW+ZY6dq1q/w5l8XnZPn9739PkUiExo8fL5eY0aNH0/PPP09BgboJAEDlJSbFrzEsSQ2ea968uWy4cFd8bOyw2/HV/ESwrKysbuHx3AAA2SRTiSPZuHHj6JtvvpE90zzN8TnnnFP3M26ELF26tG6d/5sfGtmXIDVWdFA3AQBQoILuk2qwcIX5+eefU4cOHWRMCzda4sdXc96AZcuWNTi+msdr28dwAwBkE9MwPS2QHtlYNznOkk5RZclkVnvmLKu9bklh1vKMS3xejh9UZPj9csRhJnhnNO+97fNh/wzJz1Gm3y9Rq1ki6mL/m9P8XfqVaKoxLHfeeSeNGDGCunTpIntOOIaFY054eAEP+5o4cWJdJuVYVuVmzZrRNddck75XAADQRIaEgR7qJgCA+ieEcSMQDZZt27bR1VdfLROYHXHEETJuZdWqVXXjp++++26qrKyUQxb27t0rAyEXL15MLVu2TNf5AwA0OjRYGhfqJgAAFffkG0YwHqYZggc0+wj32ByavSbUKDODA0BTwLc9U8bNJTPUJ3a/Ou6wayls5LnaNypq6LN9LyV9DpAZsfe6MK+b69wJuiEjpi1nR9Ss1uyn5vUwhbqdPf+HPq9H6vKwGBS2FTjMZ+GI5suS5vqpT4G95vBw8Ppkofp6DMpNuI2W7fUIW34Q3Tb691WkLJeOdkhWOt9XR6/vUGn6PrfqM/uQoU4rHwrlJfx86O4J2jL7dbdtw/eKyprNKaubTi682lPdtK7yFd/VTa4TRwIAgJWXIHqvQfcQDGojxnT4GXHy5Z0ClSQy03TXU3dl7F+wDYdvg5NEkX4dlpM9ySSdJMt0lphV2BNM6j8MjqQ7UaQd10xGQGYJQ4MFACBJJkXJcJmtmvcBAABIFxGgh2losAAAJOnQvCpuKwV/PsUCAIBgMGWHkNseFn9CgwUAIElBCmwEAIBgMANUN6HBAgCQkiFh7sYmY0gYAACkkxmgKffRYAEAAHCIZ/dy2zh1Enitm0mME9mphU4C8b3OHqWb1cjp7FGUxtmjRApfs0e66254PYPErzm1El8/oXvvHbxmz5My6GbQ0r73TgLxhbfXrNtE83E3bUH3ujPSDQkO2WeR4+1svR32v3s/J6FsbGiwAAAkzX1go39HCgMAQBCY6GEBAIC6SkHw0/KQh30AAADSQwRolrCmN8E6AECaKgW3CwAAQDqD7k0PixszZsygU089lVq2bElt27alSy+9lL744ouUvxY0WAAAUhCP4GUBAABIFyEHhbn7n9uHacuWLaPx48fTqlWraMmSJRSJRGjo0KFUUVGR0teCGBYAgCQdmlUlGDOxQMMi0UoP2apNbfB+4gB7TdZtbUNXZDbruz0wWHs9HH6+7UHHvshP5Cz4Ww0S9/Y37WxSgcxz9t57PbjTyRVSyX583UQXuv2qLWtRzd+qYWgC7DV/v6Zo+Gt3qoPuhXw45u5+5fZh2jvvvGNZf+6552RPy9q1a+mcc86hVEEPCwAAAAAAJK2srEz+26pVK0ol9LAAACQJme4BACBIvf/l5eWW8vz8fLk0RAhBkyZNorPOOot69+5NqYQeFgCAJAkR9bQAAACki+kxioV17tyZiouL6xYOrk/klltuoU8++YReeeWVlL8W9LAAACQJMSxNR9Q8QIbhdeC+29gN3ZPRdI7zd/YkVhkTn1Tj2+uY/VRdB6fH0cVcpFPjx3M42aueoA+PMh23ozv3qIO4kpDDeDONBPFv3EOR+hgWw/U+rKSkhIqKiurKE/WuTJgwgRYuXEjLly+nTp06UaqhwQIA0IiVAgAAgN8ephUVFVkaLA01srixsmDBAlq6dCl169aN0gENFgCAJPETONfJuVI8GwwAAECmE0fylMYvv/wyvfnmmzIXS2lpqSznYWSFhYWUKohhAQAAAAAA1+bOnStnBhs0aBB16NChbpk/fz6lEnpYAACShBgWAADwG1MOPTY87ONcquNu6oMGCwBAkuSsX25jWDBLWJaqISFSEHTvC16/aDTF+Cs/JLRMJ3wWnF8bTWB+yi6zyLohYZmCBgsAQJKQhwUAAPzGFB56WHz6MA0NFgCAlATdu+1h8edTLAAACAaBHhYAAPgBDwlzy59PsQAAIEgNlqjrffwIs4QBAAAAAIBvYUgYAECSDg3vwpAwAADwV91kBqRuQoMFACBJaLAAAIDfCDm8ywjEkDA0WAAAknToGVYwKgUAAAgG4WHGL79OuY8GCwBAktDDAgAAfmMG6GEagu4BAAAAAMC30MMCAJCkIHW7AwBAMIgATQiDBgsAQAoy3XPnu/t9AAAA0kN4yPflZZ9MQIMFACBJXp5I+fUpFgAABIMQHh6myX38Bw0WAIAkocECAAB+IzwE0Ps16B4NFgCAJAWpUgAAgGAQMlZSBKL3H7OEAQAAAACAb6GHBQAgSRgSBgAAfiMCFF+JBgsAQJKCVCkAAEAwiAANV/Zdg+WH2Qn8OUsBAASBSPFsKF5u8P6sFEAPdRMAZFvdJAL0MM13DZb9+/d//1/8ZqHRAgDpvd8UFxcnfZwgVQqgh7oJALKubiL0sKRNx44dqaSkhFq2bCnfsM6dO8v1oqIiyibl5eU4d1x3fGZ8ip9e8f2F7zcpOV6AKgXQQ93UuFCn4ro3hc9MyusmEZxZwnzXwxIKhahTp07yvw3DkP/yBy1bPmx2OHdcd3xm/CkVT6+g6UDd5A+oU3Hdg/6ZQd2UJQ0WAIBsE6RswgAAEBTCQ7ykP+smNFgAAJLG3e6HeoSzvVIAAIBgEHJ4lxGIh2m+brDk5+fT/fffL//NNjh3XHd8ZpqOIFUKkBju75mHa944cN2zmyAPdZNPH6YZArUmAIDnoE4eb2wYzepi7pziW68QB6msrCyrxlcDAEC21E3NPdZNFb6rm3zdwwIAkBU89LAQelgAACCdRHDqplBjnwAAAAAAAEB90GABAEiS8Pg/L5588knq1q0bFRQUUN++fWnFihUNbr9s2TK5HW9/1FFH0VNPPeXxVQIAQLbFsAgPSybqJrfQYAEASJrpcXFn/vz5NHHiRLrvvvto3bp1dPbZZ9OwYcNo69at2u03b95MF154odyOt7/33nvp1ltvpddeew3vOQBA4Jm+rJsC1WBJd0stFZYvX04jRoyQGUk5qOmNN95QApemTp0qf15YWEiDBg2ijRs3kh/MmDGDTj31VGrZsiW1bduWLr30Uvriiy+y4vznzp1LJ554Yl0yqAEDBtDbb7/t+/Ou733gzw7/ofv9/Pmc+Fzjl/bt2/v+vDNDHBr362bx0MMya9Ysuummm+jmm2+mXr160ezZs2UmZ/6b0OHelC5dusjteHve78Ybb6Tf/va3KXjNTRPqpvRC3dT4sqleYqibsq9uCkyDJRMttVSoqKigPn360Jw5c7Q/f+SRR+SbyD9fvXq1/HI3ZMgQ2r9/PzU2HiYyfvx4WrVqFS1ZsoQikQgNHTpUvia/n3+nTp3o4YcfpjVr1sjl3HPPpUsuuaTu5unX87bjc5s3b55sfMXz8/kff/zxtGPHjrplw4YNWXHe6Sc8dLmLutlc4pfq6mrtb6ipqaG1a9fKv9N4vL5y5UrtPh999JGy/QUXXCD/bmpra1P26psK1E3ph7qpcWVjvcRQN2VX3eSJ8KHTTjtNjB071lLWs2dPMXnyZOFXfCkXLFhQt26apmjfvr14+OGH68qqqqpEcXGxeOqpp4Tf7Nq1S76GZcuWZeX5H3744eLpp5/OmvPev3+/6N69u1iyZIkYOHCguO2222S5n8///vvvF3369NH+zM/nnU6VlZXydX9/h3e9tGjRQinj66yzfft2+fMPP/zQUv7rX/9a9OjRQ7sPf8b45/F4fz7Ot99+m8Ir0TSgbso81E2Zk431EkPdlH11kxe+62HJWEstzXjseGlpqeV1cAKmgQMH+vJ18HzbrFWrVll1/tFolF599VXZM8RDw7LlvLl3a/jw4XT++edbyv1+/l999ZUcEsDDNa+66iratGlTVpx3uvCQVX7t/PfjZdm2bZtSNmXKlAZ/p31OfX5e0tA8+7rtdeXQMNRNjQN1U+Zka73EUDdlX92U9XlYdu/eLb+EtmvXzlLO6/wHky1i56p7HVu2bCE/4Q/VpEmT6KyzzqLevXtnxfnzUCRuoFRVVVGLFi1owYIFdNxxx9XdPP163owbWB9//LHsVrfz83Xv378/vfjii9SjRw/auXMnTZ8+nc444ww5FM/P552JioGXdGvTpg2Fw2HlPrhr1y7lusfwsA3d9jk5OdS6deu0nm/QoG7KPNRNmZOt9RJD3ZR9dVMgGiyZaqllSja8jltuuYU++eQT+uCDD7Lm/I899lhav3497du3T854NHr0aDn22e/nXVJSQrfddhstXry4wRuJH8+f48hiTjjhBNlgPProo+mFF16g008/3bfnHRR5eXlyAhKOObvsssvqynmdY7h0+D3661//ainjz16/fv0oNzc37eccREH5jGfD60DdlBnZXC8x1E3ZVzd54bshYZlqqaVbbPYkv7+OCRMm0MKFC+n999+XwezZcv78B3LMMcfIL148owlPfvDYY4/5/rx5uCOfC/9x81NuXrih9fjjj8v/jp2jX88/XvPmzWXDhbvi/X7dg4J7Qp9++ml69tln6fPPP6fbb79dTkYyduxY+XPusr/uuuvqtudyfgLK+/H2vN8zzzxDd955ZyO+iuyEuimzUDdlTpDqJYa6yX91UyAbLPEttXi8zsNPsgWP8ecvcfGvg8dA803AD6+Dn4zw06vXX3+d3nvvPXm+2XT+utfDM1j4/bzPO+88OZyNe4diCze6rr32WvnfnNjPz+cfj68335g6dOjg++seFCNHjpTTRU6bNo1OOukkObX6okWLqGvXrvLnPHNb/GyK/L7wz5cuXSq3f/DBB+WXkJ/85CeN+CqyE+qmzEDdlHlBqpcY6ib/1U0pIXzo1VdfFbm5ueKZZ54Rn332mZg4caJo3ry5+Oabb4Sf8Iwa69atkwtfylmzZsn/3rJli/w5z6jBs2i8/vrrYsOGDeLqq68WHTp0EOXl5Y196uLnP/+5PLelS5eKHTt21C0HDx6s28av5z9lyhSxfPlysXnzZvHJJ5+Ie++9V4RCIbF48WJfn3d94mdj8fP533HHHfLzsmnTJrFq1Spx0UUXiZYtW9b9Xfr1vAFSBXVT+qFu8odsqZcY6qamwZcNFva73/1OdO3aVeTl5YlTTjmlbrpdP3n//fe108GNHj26bipAngaOp5bLz88X55xzjvxD94P6prJ77rnn6rbx6/nfeOONdZ+NI444Qpx33nl1jRU/n7fTisGv5z9y5EhZQfHDhI4dO4rLL79cbNy40ffnDZBKqJvSC3WTP2RLvcRQNzUNBv9f6vprAAAAAAAAKLgxLAAAAAAAADFosAAAAAAAgG+hwQIAAAAAAL6FBgsAAAAAAPgWGiwAAAAAAOBbaLAAAAAAAIBvocECAAAAAAC+hQYLAAAAAAD4FhosAAAAAADgW2iwAAAAAACAb6HBAgAAAAAAvoUGCwAAAAAAkF/9/5DR4tH82BMWAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "rho_opt = np.asarray(jax.nn.sigmoid(latent))\n", "rho_grid = np.zeros(Nx * Ny)\n", "rho_grid[des_idx] = rho_opt\n", "rho_grid = rho_grid.reshape(Nx, Ny)\n", "\n", "chi_opt = np.asarray(density_to_chi(jnp.asarray(rho_opt))).reshape(Nx, Ny)\n", "Ez_opt = solver.get_TM_field(ji, chi_opt) # numpy solve for plotting\n", "\n", "fig, ax = plt.subplots(1, 2, figsize=(10, 4))\n", "im0 = ax[0].imshow(rho_grid, cmap=\"magma\")\n", "ax[0].set_title(\"optimized density\")\n", "plt.colorbar(im0, ax=ax[0])\n", "im1 = ax[1].imshow(np.abs(Ez_opt) ** 2, cmap=\"inferno\")\n", "ax[1].set_title(r\"$|E_z|^2$ (total field)\")\n", "plt.colorbar(im1, ax=ax[1])\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "12", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "keep_output": true, "kernelspec": { "display_name": "dolphindes", "language": "python", "name": "python3" }, "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.10.19" } }, "nbformat": 4, "nbformat_minor": 5 }