{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Constructing Scalar Models\n", "\n", "Once the modes have been found, we can then make use of them to construct scalar models of our meta-atoms." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# setup 2D plotting \n", "%matplotlib inline\n", "from openmodes.ipython import matplotlib_defaults\n", "matplotlib_defaults()\n", "import matplotlib.pyplot as plt\n", "\n", "# the numpy library contains useful mathematical functions\n", "import numpy as np\n", "\n", "# import useful python libraries\n", "import os.path as osp\n", "\n", "# import the openmodes package\n", "import openmodes\n", "from openmodes.sources import PlaneWaveSource" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialising simulation" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Pruned cell types: vertex, line\n" ] } ], "source": [ "sim = openmodes.Simulation(notebook=True)\n", "mesh = sim.load_mesh(osp.join(openmodes.geometry_dir, \"SRR.geo\"), parameters={'inner_radius': 2.5e-3}, mesh_tol=0.8e-3)\n", "ring = sim.place_part(mesh)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Searching for modes" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "26ba383ee27e41f1a40c88fd073ef7a3", "version_major": 2, "version_minor": 0 }, "text/plain": [ "HBox(children=(Label(value='Refining modes'), FloatProgress(value=0.0, max=3.0)))" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "start_freq = 2e9\n", "start_s = 2j*np.pi*start_freq\n", "\n", "num_modes = 4\n", "estimates = sim.estimate_poles(start_s, modes=num_modes, cauchy_integral=False)\n", "refined = sim.refine_poles(estimates, rel_tol=1e-12)\n", "modes = refined.add_conjugates()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZwcdZ3/8debkEgOAoQMAoEEOUQDrIIDhEVdRIUQQbKsB6wrh2gEPHHxh7KKioi6q+LBEcLCj6AYDg8MbDi8uH5LCEk4IgGEQEIGAplMSMJkAiHh8/ujvhM6nZ6Z7tAzVTPzfj4e/Ziq+n676tPVPfXpqvr296uIwMzMrGi2yDsAMzOzSpygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygLFeSFkr6QA9v8xZJJ/XAdt4s6S5JL0n6cXdvrydI+oSk23t4m/8sabGkVkn79/C275D06Z7cpr1uy7wDsHxI+lfgK8DbgJeAB4HvRcQ9uQbWAyLiqB7a1CRgGTA8+sgPDiPiGuCaHt7sj4DPR8Qfeni7ljOfQfVDkr4C/BS4AHgzMBq4BDg2z7j6oDHA/I6SkyR/QazOGOCRvIOwHESEH/3oAWwDtAIf7aTOm8gS2HPp8VPgTansMKAJ+D/AUmAJMBGYAPwdWA6cU7KubwO/Aa4jO1ObC7yjpHwh8IE0vQXwNWAB0AJcD4xIZR8HniI7GwE4CngeaOjgNYwD/hdYATwEHFZSdgfw6TQ9APgx2ZnO08DngQC2LNlfV6TX+SxwPjAglZ0M3EP2Df/F9PyjUtlVwKvA2rS/P1CyL34FrAI+3dlrTuv5JLAolf1H2f66Cji/pO5hQFPJ/M7Ab4HmFNsXy96X64Gr0/vyCNBYUr4r8Lv03BbgotLXXFLvbcAf0/v+OPCxkrIJwPy0/meBszp4r7YAvpFe59IU0zZkn8PW9H6sBhZ08PwAzgCeSNv6LrAHcG/az9cDg0rqfwZ4MsU8Hdi5pOyDwGPASuAi4E7SZyWVfwp4NL3ftwFj0nIBF6b4VwIPA/vm/f/e2x+5B+BHD7/hMB5YRzoAd1DnPGAmsAPQQHag/24qOyw9/1xgYPpnbwZ+DWwN7AO8DOye6n+b7ED9kVT/rHSwHJjKF/L6AffLabu7pIPTZcC0kriuITsob0+WOI/uIP5R6aA6IR38PpjmG1L5HbyeoE5LB9FdgO2AP7FxgroxxTE07Y9ZwGdT2cnptX2GLNGdnuJSKr+KjRNI+76YmOIa3NlrBsaSHaDfm8p+kvZ9lwkqrX9Oep8GAbuTJfgjS2J5Oe2jAcD3gZmpbABZUr8wve6tgHeXvOZ70vRQYDFwCtntggPIEv0+qXwJ8J40vR1wQAfv16fIEsbuwDCyxPjLkvIA9uzk8xpkiWY42efvFeDPaX3bpPf3pFT38BTjAWmf/gK4K5WNJEto7Z/VM9P+bv+sTExxvj293m8A/5vKjkz7e1uyZPV2YKe8/997+yP3AOr+guBKsm8xf6ui7nvJvtGvAz5SVnYS2TeyJ9o/3H3hAXwCeL6LOguACSXzRwIL0/RhwBpeP4vYOh0gDi6pPweYmKa/3X7gS/NblB24FvL6AfdR4P0ldXciO6C3J4ttgWeAecBlncR/dukBLi27reQgdUfJQecvpIST5j+QXs+WZJc/XwEGl5SfAPw1TZ8MPFlSNiQ9d8c0fxWbJqi7yuLq8DWTJZdrS8qGkp2RVZOgDgaeKdvW14H/WxLLn0rKxgJr0vQhZF86NvkSw8YJ6uPA3WXllwHfStPPAJ8lnfV28n79GTijZH7vsve9mgR1aNnn7+yS+R8DP03TVwD/WVI2LG1rN+BENv6siuxqQftn5Rbg1LLPchvZJcjDya4gjAO26Mn/6b786Iv3oK4iO0uoxjNk/3C/Ll0oaQTwLbJ/8oOAb0narn4h5qoFGNnF/Y+dyS63tFuUlm1YR0SsT9Nr0t8XSsrXkP3jt1vcPhERr5H905eur90Y4PeSVkhaQXbwXk+WKIiIFcANwL5kB52OjAE+2r6etK53kx38K73WxSXzpdNjyL5JLylZz2VkZ1Ltni95bW1psvS1l1tcNt/Za94otohYTfb+VWMMsHPZPjgnrXeT2MkOtFulz8WuwKKIWFfFNg4u28YngB1T+b+QnaEtknSnpEM6WE+lz1v7F4RqlX/+Ovo8brStiGgl26ej2HR/B5t+Hn5W8lqXkyWxURHxF7JLghcDL0iaIml4DfFbBX0uQUXEXWQfnA0k7SHpVklzJN0t6W2p7sKIeBh4rWw1RwJ/jIjlEfEi2TX2apNe0d1LdmlnYid1niP7Z2w3Oi3bXLu2T0jaguxyVqX1LSa7h7NtyWOriHg2PfedZJeDpgE/72R7i8nOoErXMzQiflCh7pIUzyaxpvW8AowsWc/wiNinitfckagQa0eveQkb77shZJc3260mO2trt2PJ9GLg6bL1bh0RE6qIcTEwuopGHIuBO8u2MSwiTgeIiPsj4liyhH4j2b2gSip93taxcZKpl422JWko2T6ttL/Fpp+Hz5a93sER8b8AEfHziHgX2WXGtwJf7Yb4+5U+l6A6MAX4QvrwnEXWYq0zo9j4m1NTWtbrRcRKsktHF0uaKGmIpIGSjpL0n6naNOAbkhokjUz1f/UGNvsuScelA96XyQ76MyvUmwx8T9IYgLT9Y9P0VimGc8jueYySdEYH2/sVcIykIyUNkLSVpMMk7VKh7vXAlySNkrQt2eVBACJiCXA78GNJwyVtkb7s/NPm7IQOdPiayRpUHC3p3ZIGkd0bLP2ffRCYIGmEpB3J9m27WcAqSWdLGpz2w76SDqwipllkB+sfSBqa9t+hFerdDLxV0ifTZ2igpAMlvV3SoPSbqW0i4lWyezvrK6wDss/bmZLeImkYWevS66o4g9scvwZOkfROSW9K27ovIhYC/wPsU/JZ/SIbJ/3JwNcl7QMgaRtJH03TB0o6WNJAsi8OL9Px67Uq9fkElT7w/wjcIOlBsks0lS71bPS0CsvKv/n2WhHxE7LfQH2D7F7DYrLWazemKucDs8laIs0ju093/hvY5B/I7le8SNYq7bh00Cr3M7Kb3bdLeoksiR2cyr5Pdn/l0oh4Bfg34HxJe1V4fYvJmsyfU/L6vkrlz/vlZEnoYeABYAbZt/f2g8uJZI0M5qf4f0PXn59adPiaI+IR4HNkB9UlaftNJc/9JVljhoXpNVzXXpAuwR4DvJOsUcoy4L/JGg10quS5e5JdBm8ie//K670EHAEcT3Zm8jzwQ7LGB5C91wslrSJrjPJvHWzyyvRa7kqxvgx8oas4N0dE/Bn4JlnrxiVkrf2OT2XLgI8CPyC77LcX8P9Knvt7std3bXpNfyNrTQpZA43Lyd6j9laXP+qO19CftLc26lMk7QbcHBH7puvAj0dEhwcVSVel+r9J8yeQNUv+bJq/DLgjIqZ1d+x9jaRvk93g7ujgVCiSjgImR8SYLivnQNJCspv2f8o7FrPu1ufPoCJiFfB0yam4JL2ji6fdBhwhabvUOOKItMz6mHT5a4KkLSWNImsc8/u84zKzPpigJE0jawiwt6QmSaeStSw6VdJDZD9IbL+vcaCkJrLT+sskPQIQEcvJfux3f3qcl5ZZ3yPgO2SXZh4ga0V3bq4RmRnQRy/xmZlZ79fnzqDMzKxv6FOdVY4cOTJ22223vMMwM7MazJkzZ1lENJQvzzVBSRpA1pz52Yg4uqxMZE1wJ5D9yv3kiJjb2fp22203Zs+e3V3hmplZN5C0qNLyvC/xfYnspnQlR5H9DmEvsnF1Lu2poMzMLH+5Jaj0q/4Pkf14sJJjgasjMxPYVlI9fyBpZmYFlucZ1E/JxhQq7wevXVXdDUmaJGm2pNnNzc31j9LMzHKRS4KSdDSwNCLmdFatwrJN2sRHxJSIaIyIxoaGTe6xmZlZL5XXGdShwIdTty3XAodLKu+MtImNexLuqAfsN2zWrFm0trZumG9tbWXWrFndsSkzs15v+TXXsK7l9ZFf1rW0sPyaa+q+nVwSVER8PSJ2iYjdyDpq/EuFvtqmAyemronGAStT79J1NWvWLGbMmMHUqVNpbW2ltbWVqVOnMmPGDCcpM7Myy6+5hhe+ez6LTjqJdS0trGtpYdFJJ/HCd8+ve5Iq1O+gJJ0GEBGTyXqVnkA2xHIb2RALdTd27Fjuv/9+mpubueSSbBSOtrY2GhoaGDt2bHds0sys1xo+fjwvTpvG2icX8NQxHwZg/fLlDNpzD4aPr++weX2qq6PGxsbYnN9Btba2cskll9DWlg2IOmTIEM444wyGDetsYFQzs/5pXUsLTx3zYdYvz7ooHTBiBLvfNJ0tt9++i2dWJmlORDSWL8/7d1BmZmYV9fsE1X7Pqa2tjSFDhjBkyBDa2to23JMyM7PXtd9zWr98OQNGjGDAiBGsX758wz2peur3CWr+/Pk0NzfT0NDAGWecwRlnnEFDQwPNzc3Mnz8/7/DMzApl1a23svbJBQzacw92v2k6u980nUF77sHaJxew6tZb67ot34Mia8k3duzYDfecWltbmT9/PgcddFC9QzQz6/WWX3MNw8eP33DPaV1LC6tuvZURn/jEZq2vo3tQTlBmZpYrN5IwM7NexQnKzMwKyQnKzMwKyQnKzMwKyQnKzMwKyQnKzMwKyQnKzMwKyQnKzMwKyQnKzMwKyQnKzMwKyQnKzMwKKZcEJWkrSbMkPSTpEUnfqVDnMEkrJT2YHufmEauZmeUjryHfXwEOj4hWSQOBeyTdEhEzy+rdHRFH5xCfmZnlLJcEFVkX6u2jAQ5Mj77TrbqZmb1hud2DkjRA0oPAUuCPEXFfhWqHpMuAt0jap4dDNDOzHOWWoCJifUS8E9gFOEjSvmVV5gJjIuIdwC+AGyutR9IkSbMlzW5ubu7eoM3MrMfk3oovIlYAdwDjy5aviojWND0DGChpZIXnT4mIxohobGho6ImQzcysB+TViq9B0rZpejDwAeCxsjo7SlKaPogs1paejtXMzPKRVyu+nYCpkgaQJZ7rI+JmSacBRMRk4CPA6ZLWAWuA46MvjU9vZmadyqsV38PA/hWWTy6Zvgi4qCfjMjOz4sj9HpSZmVklTlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZIuSQoSVtJmiXpIUmPSPpOhTqS9HNJT0p6WNIBecRqZmb5yGXId+AV4PCIaJU0ELhH0i0RMbOkzlHAXulxMHBp+mtmZv1ALmdQkWlNswPTI8qqHQtcnerOBLaVtFNPxtnTWu99jvWtazfMr29dS+u9z+UYkZkVybTHptGypmXDfMuaFqY9Ni3HiLpXXmdQSBoAzAH2BC6OiPvKqowCFpfMN6VlS8rWMwmYBDB69Ohui7e7td77HCv+sIDWe5fQMGk/AJqnzGPd0jYAhh2yc57hmVnOpj02jQvuu4DrHruOK468AoBTbzuVBSsXAHDC207IM7xukVuCioj1wDslbQv8XtK+EfG3kiqq9LQK65kCTAFobGzcpLy3GLzfSFrvXcK6pW28cOFcAF5b/Spb7jCEwfuNzDk6M8vbEWOO4LrHrmPBygUcN/04AJa/vJw9ttmDI8YckXN03aPLBCXp3BrXeUdE3FVt5YhYIekOYDxQmqCagF1L5ncB+uz1rgHDBtEwaT9euHAur61+FYAthg6kYdJ+DBg2KOfozCxv2w/eniuOvILjph/H8peXAzBiqxFcceQVbD94+5yj6x7V3INSjY+uVyg1pDMnJA0GPgA8VlZtOnBias03DlgZEUswM7N+ocszqIjYpAl4HewETE33obYAro+ImyWdlrY5GZgBTACeBNqAU7ohjsJY37qW5inzeG31q2wxdCCQXeJrnjLPZ1FmRsuaFk697VSWv7ycEVuNALJLfKfedmqfPYuq6R6UpA8B7wOGkCWOaRGxRNL/RMSHql1PRDwM7F9h+eSS6QA+V0t8vdmaectYt7SNLXcYskkjiTXzlrmRhFk/d/ui21mwcgF7bLPHJo0kbl90e59sJKEsD3RRSdqK7JLbQcB9wApgH2AM8FlgckQM78Y4q9LY2BizZ8/OO4zN1nrvcwzeb+SGs6X1rWudnMxsg2mPTeOIMUdsOFtqWdPSJ5KTpDkR0bjJ8ioT1H+RJad/iYhlJcsPB64Cdo6I3FoEtuvtCcrMrD/qKEFV+0PdjwOfKU1OABHxF+AYoOpWe2ZmZtWo9qxnBPBEpYKIeAg4vG4RmZmZUf0Z1CLgXZUKJE2UtK5+IZmZmVWfoH4KXCVpbPsCSYMknU7Wi8Nr3RGcmZn1X1Vd4ouIyyU1ALMkLQFWkvWh10T2W6W/dl+IZmbWH1Xd8i4iLpB0KXAoMAx4KiJmAUi6oJviMzOzfqqmpuER8SJwc4Xl369bRGZmZhSgs1gzM7NKqjmDqqoDWDMzs3rKq7NYMzOzTuUy5LuZmVlXnKDMzKyQnKDMzKyQnKDMzKyQak5Qkj4o6QpJN6X5xjTshpmZWd3UlKAkfQG4lKxn8/emxWuA82tcz66S/irpUUmPSPpShTqHSVop6cH0qPX3WGZm1ovVOsjgl4H3R8RCSWenZY8Be9e4nnXAv0fEXElbA3Mk/TEi5pfVuzsijq5x3WZm1gfUeolva2Bxmm4fincgsLaWlUTEkoiYm6ZfAh4FRtUYi5mZ9WG1Jqi7gK+VLfsib6A3c0m7AfsD91UoPkTSQ5JukbRPB8+fJGm2pNnNzc2bG4aZmRWMIqLrWu2VpZ2Am4CRZGc8TwGrgGMi4vmaNy4NA+4EvhcRvysrGw68FhGtkiYAP4uIvTpbX2NjY8yePbvWMMzMLEeS5kREY/nyWnszXyLpQOBAYAzZ5b5ZEVHzgIWSBgK/Ba4pT05pW6tKpmdIukTSyIhYVuu2zMys96m1Fd9ZkZkVETdExMyIeE3SV2pcj4ArgEcj4icd1Nkx1UPSQSnWllq2Y93vgdtupm3lig3zbStX8MBtm4zIYlY8sy6H1pLbAq3N2TIrjFpb8Z0L/KjC8m8AFRNNBw4FPgnMk/RgWnYOMBogIiYDHwFOl7SOrCn78VHL9Ujrdg/cdjN/uXIyD90+g4+dm41Zef1559DS9AwA+x/pBphWULMuhxlnwf3/DSelL1RTj4bmx7Lpgz6TX2y2QVUJquSHuAMkvY+Nh+DYHXiplo1GxD10MYxHRFwEXFTLeq1n7T3u3Tx0+wxamp7hqrM+B8CaVSvZfpfR7D3u3TlHZ9aJsROz5NT8GFwyLlvWtgwa3paVWSFU1UhC0tNpcjTwTElRAC8A34+I6fUPrzZuJNHz2lau4KqzPseaVSsBGDx8G07+0cUM2WbbnCMz60Jrc5ac2tJt7SEj4YyZMKwh37j6oTfUSCIi3pJWcnVEnFjv4MzMzMrV1EgiIk6U9GZJx0g6RdKn2h/dFaAVV9vKFVx/3jmsWbWSwcO3YfDwbVizaiXXn3fORg0nzAqntTm759S2LDtzGjIym5569MYNJyxXtbbimwgsAM4DLgO+kP5+sv6hWdE9PvMeWpqeYftdRnPyjy7m5B9dzPa7jKal6Rken3lP3uGZdWz+jdn9p4a3ZZf1zpiZTTc/lpVZIdTaiu984JSIuEHSixGxv6RTgIq9PFjf1t5Kb+9x795wz+lj517A4zPvcQs+K7b2VnpjJ75+z+mkm7Pk5BZ8hVFrTxKrImJ4mn4xIraTtAXwfETs0F1BVsuNJMzMep+OGknU2hffUklvTtMLJR0C7AEMeKMBmpmZlao1QV0OtP/A5UKyTmIfAi6pZ1BmZma13oP6r/Z+9yLiakl3AEMj4tG6R2ZmZv1a1QlK0gCgVdK2EfEKQEQ808XTzMzMNkvVl/giYj3wd2D77gvHzMwsU+slvmuAmyX9DGji9VF1iYi/1DMwMzPr32pNUKenv98uWx5kncaamZnVRa0DFr6luwIxMzMrVWszczMzsx7hBGVmZoWUS4KStKukv0p6VNIjkr5UoY4k/VzSk5IelnRAHrGamVk+am0kUS/rgH+PiLmStgbmSPpjRMwvqXMUsFd6HAxcmv6amVk/kMsZVEQsiYi5afol4FFgVFm1Y4GrIzMT2FbSTj0cqpmZ5aSmMyhJ53VQ9ArZ76JujYgXalznbsD+wH1lRaOAxSXzTWnZkrLnTwImAYwePbqWTZuZWYHVegb1VuBs4H3Anunv2WQJ5nTgKUnjq12ZpGHAb4EvR8Sq8uIKT9lkbJCImBIRjRHR2NDQUO2mzcys4GpNUFsAx0fEeyLiXyPiPcDHgPURMQ44A/hBNSuSNJAsOV0TEb+rUKUJ2LVkfhfguRrjNTOzXqrWBHUkML1s2c1kDRoAfkU2PlSnJAm4Ang0In7SQbXpwImpNd84YGVELOmgrpmZ9TG1tuJbQHYp76KSZael5QAjgdVVrOdQ4JPAPEkPpmXnAKMBImIyMAOYADwJtAGn1BirmZn1YrUmqE8Dv5N0NvAs2WW3dcBxqXxv4JtdrSQi7qHyPabSOgF8rsb4zMysj6i1L765kvYCxgE7k7WouzciXk3ldwF31T1KMzPrd2ptZj4IOBl4JzAsLf60JCLixDrHZmZm/Vitl/imAu8AbgJq+r2TmZlZLWpNUOOBt0TEiu4IxszMrF2tzcyfAd7UHYGYmZmVqvUM6mrgD2nI940u8XnIdzMzq6daE9Tn098LypZ7yHczM6srD/luZmaF5BF1zcyskLo8g5L03vQDXCQd3lE934MyM7N6quYS3yXAvmn6ig7q+B6UmZnVVZcJKiL2LZn2PSgzM+sRNd2DkjRI0nmSnpC0Ov39rqStuitAMzPrn2ptZj6ZbFTdLwKLgDHA18mGYv9UfUMzM7P+rNYEdSywR0lXR/Ml3Uc2ZpMTlJmZ1U2tzcyfB4aULRtMNuyGmZlZ3dR6BvVL4FZJvwCagF3JBhW8urQJupucm5nZG1Vrgvps+ntO2fLT0gOqaHIu6UrgaGBpaSvBkvLDgD8AT6dFv4uI82qM1czMerG8ujq6CriIrPPZjtwdEUfXaXtmZtbL5NLVUeqZYnke2zYzs96h1iHftyFrYr4/rw/5DkBEHFHHuAAOkfQQ8BxwVkQ80kFMk4BJAKNHj65zCGZmlpda70HdAAwAfg+sqX84G8wFxkREq6QJwI3AXpUqRsQUYApAY2NjdGNMZmbWg2pNUOOA7SPi1e4Ipl1ErCqZniHpEkkjI2JZd27XzMyKo9Z7UPcAb++OQEpJ2lGS0vRBZHG2dPd2zcysOGo9gzoZmJF6jygf8r3qZuCSpgGHASMlNQHfAgam9UwGPgKcLmkd2aXE4yPCl+/MzPqRWhPU98h+nLsQGF6yvKbkEREndFF+EVkzdDMz66dqTVDHA2+NCHdtZGZm3arWe1BPAd3aQMLMzAw2ry++6akvvvJ7UO5/z8zM6qbWBPW59PeCsuUe8t3MzOoqr774zMzMOlXrGRSS3gwcBIwE1L48Iq6sY1xmZtbP1doX30TgV8ATwD7AI8C+ZD/gdYIyM7O6qbUV3/nAKRGxP7A6/Z0EzKl7ZGZm1q/VmqBGR8QNZcumAifWKR4zMzOg9gS1NN2DAlgo6RBgD7Iezs3MzOqm1gR1OfDuNH0h8FfgIeCSegZlZmZWazPzH5ZMXy3pDmBoRDxa78DMzKx/q7UV31jgPcAIsiHb746I+d0RmJmZ9W9VJag0NtMVwElAE9kw7KOAnSX9EviUh8MwM7N6qvYe1CSy8ZvGRcSYiDgkIkYDh5CdUX22m+IzM7N+qtoE9UngixFxf+nCNP/lVG5mZlY31SaoscCdHZTdmcqrJulKSUsl/a2Dckn6uaQnJT0s6YBa1m9mZr1ftQlqQES8VKkgLa+1ufpVwPhOyo8C9kqPScClNa7fzMx6uWpb8Q2U9D5KOofdzPUAEBF3SdqtkyrHAlenhhczJW0raSeP5Gtm1n9Um1iW0nlnsEvrEEupUcDikvmmtGyTBCVpEtlZFqNHj65zGGZmlpeqElRE7NbNcZSrdKZWsRl7REwBpgA0Nja6qbuZWR9R672jntIE7FoyvwvZb6/MzKyfKGqCmg6cmFrzjQNW+v6TmVn/UvOIuvUgaRrZD39HSmoCvgUMBIiIycAMYALwJNAGnJJHnGZmlp9cElREnNBFeQCf66FwzMysgIp6ic/MzPo5JygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyuk3BKUpPGSHpf0pKSvVSg/TNJKSQ+mx7l5xGlmZvnIZch3SQOAi4EPAk3A/ZKmR8T8sqp3R8TRPR6gmZnlLq8zqIOAJyPiqYhYC1wLHJtTLGZmVkB5JahRwOKS+aa0rNwhkh6SdIukfSqtSNIkSbMlzW5ubu6OWM3MLAd5JShVWBZl83OBMRHxDuAXwI2VVhQRUyKiMSIaGxoa6hymmZnlJa8E1QTsWjK/C/BcaYWIWBURrWl6BjBQ0sieC9HMzPKUV4K6H9hL0lskDQKOB6aXVpC0oySl6YPIYm3p8UjNzCwXubTii4h1kj4P3AYMAK6MiEcknZbKJwMfAU6XtA5YAxwfEeWXAc3MrI9SXzrmNzY2xuzZs/MOw8zMaiBpTkQ0li93TxJmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlDWL827o4m2VWs3zLetWsu8O5pyjKh3uvrehSxrfWXD/LLWV7j63oV5hWN9TC6dxZrlad4dTdx17d+Zd+ezTDxzfwBuvPABXlyyGoD9Dtslz/B6javvXci5f3iEX967iGmTxgFwwpSZPLG0FYATD9ktv+CsT3CCsn5njwN2YN6dz/LiktVc+937AFjz0qtst9NQ9jhgh5yj6z0m7LcTv7x3EU8sbeXIC+8CoCSDxo4AAAmpSURBVGX1WvbaYRgT9tsp5+isL/AlPut3hgwfxMQz92fw1gNZ89KrrHnpVQZvPZCJZ+7PkOGD8g6v1xg57E1MmzSO7YcOomX1WlpWr2X7oYOYNmkcI4e9Ke/wrA9wgjIzs0JygrJ+p23VWm688IENZ07tZ1I3XvjARg0nrHPLWl/hhCkzN5w5tZ9JnTBl5kYNJ8w2V24JStJ4SY9LelLS1yqUS9LPU/nDkg7II07rexbMXcqLS1az3U5DOf6bB3P8Nw9mu52G8uKS1SyYuzTv8HqNGfOW8MTSVvbaYRi3nflebjvzvey1wzCeWNrKjHlL8g7P+oBcGklIGgBcDHwQaALulzQ9IuaXVDsK2Cs9DgYuTX/N3pD2Vnp7HLDDhntOE8/cnwVzl7oFXw3aW+lN2G+nDfecpk0ax4x5S9yCz+oilyHfJR0CfDsijkzzXweIiO+X1LkMuCMipqX5x4HDIqLDr2Ye8t3MrPcp2pDvo4DFJfNNaVmtdczMrI/KK0GpwrLyU7lq6iBpkqTZkmY3NzfXJTgzM8tfXgmqCdi1ZH4X4LnNqENETImIxohobGhoqHugZmaWj7wS1P3AXpLeImkQcDwwvazOdODE1JpvHLCys/tPZmbWt+TSii8i1kn6PHAbMAC4MiIekXRaKp8MzAAmAE8CbcApecRqZmb5yKUVX3eR1AwsyjuOGo0EluUdxGZw3D2rt8YNvTd2x91zxkTEJvdo+lSC6o0kza7UvLLoHHfP6q1xQ++N3XHnz10dmZlZITlBmZlZITlB5W9K3gFsJsfds3pr3NB7Y3fcOfM9KDMzKySfQZmZWSE5QZmZWSE5QfUwSR+V9Iik1yR12BRU0kJJ8yQ9KCn3LtpriLvTcb56mqQRkv4o6Yn0d7sO6hVif/fWcdKqiPswSSvT/n1Q0rl5xFlO0pWSlkr6WwflRd3fXcVdyP1ds4jwowcfwNuBvYE7gMZO6i0ERuYdby1xk/UKsgDYHRgEPASMzTnu/wS+lqa/BvywqPu7mv1H1rvKLWSdKY8D7ivAZ6OauA8Dbs471gqxvxc4APhbB+WF299Vxl3I/V3rw2dQPSwiHo2Ix/OOo1ZVxn0Q8GREPBURa4FrgWO7P7pOHQtMTdNTgYk5xtKVavbfscDVkZkJbCtpp54OtEwR3/eqRMRdwPJOqhRxf1cTd5/gBFVcAdwuaY6kSXkHU6UijuH15kidDKe/O3RQrwj7u7eOk1ZtTIdIekjSLZL26ZnQ3rAi7u9q9cb9vZFcOovt6yT9CdixQtF/RMQfqlzNoRHxnKQdgD9Keix9a+o2dYi7qjG86q2zuGtYTY/v7wrqNk5aD6smprlk/a21SpoA3Ajs1e2RvXFF3N/V6K37eyNOUN0gIj5Qh3U8l/4ulfR7ssso3XrArEPcVY3hVW+dxS3pBUk7RcSSdGlmaQfr6PH9XUHdxknrYV3GFBGrSqZnSLpE0siIKHqnpkXc313qxft7I77EV0CShkraun0aOAKo2FqnYKoZ56unTQdOStMnAZucCRZof/fWcdK6jFvSjpKUpg8iO/a09HiktSvi/u5SL97fG8u7lUZ/ewD/TPat7BXgBeC2tHxnYEaa3p2sJdRDwCNkl9gKH3eanwD8naxVVxHi3h74M/BE+juiyPu70v4DTgNOS9MCLk7l8+ikJWjB4v582rcPATOBf8w75hTXNGAJ8Gr6fJ/aS/Z3V3EXcn/X+nBXR2ZmVki+xGdmZoXkBGVmZoXkBGVmZoXkBGVmZoXkBGVmZh3qqmPasrrvlTRX0jpJHykru1XSCkk3V7ttJygzM+vMVcD4Kus+A5wM/LpC2X8Bn6xlw05QZmbWoajQMa2kPdIZ0RxJd0t6W6q7MCIeBl6rsJ4/Ay/Vsm0nKDOrSNKjkt6RdxxWSFOAL0TEu4CzgEu6YyPui8+shKSFwJuB9SWL3xqpr77+JCLenncMVjyShgH/CNyQelMCeFN3bMsJymxTx0TEnzoqlLRlRKzryYDMCmQLYEVEvLMnNmRmXUhDwp8t6WFgtaQtJe0s6beSmiU9LemLZc/ZP7VoeknSdZKulXR+KgtJe5bUvaq9LM13uO4Uy1lpCPKVad1blZTvKul36bktki5Ky78q6bdlMf5C0k8rvN5PSLqtZH5rSc9K+qeyerum17L9ZuxW64Ui6yn9aUkfBUgd6XbLpWAnKLPqnQB8CNiW7CbwTWSdcY4C3g98WdKRAKlX7xuBXwIjgBuAf6lmI5K26GzdycfIWla9BfgHspZTSBoA3AwsAnZLz782PedXwHhJ26a6WwIfTzGW25eNe3T/d7Lhxe8srRQRi4HVwH7VvDbrfSRNA+4F9pbUJOlU4BPAqZLaO1g+NtU9UFIT8FHgMkmPlKznbrL/g/en9RxZvq1yvsRntqkbJbVfwrsjItqHif95OiAj6WCgISLOS2VPSbqcbKiJ24BxwEDgp5H1yPwbSV+pcvsHdrHu9lieS7HcBLRfbjmIrKf2r5ZchrwHshGFJd1FdvC4nCzBLYuIORVi2Bf4bVr/AOB0sh6ykdQAbB0RT6W664DBVb4262Ui4oQOijZpeh4R95ONmVVpPe+pddtOUGabmtjBPajSob/HADtLWlGybABwd5reGXg2Nh4uYFGV2+9q3QDPl0y3pe1BNrjeok7ukU0lSzaXA/9G5bMnyBLUd0qmd+D15PgVsmEoviZpMLA1HQwEafZG+BKfWfVKk81i4OmI2LbksXVETEjlS4BR7YPGJaNLptuAISXzpUPWd7XuziwGRqfLd5XcCPyDpH2Bo4FryiukVlq7AvPTolHAi/H6KK3jeT0h/RPwIvBAFbGZ1cQJymzzzAJWpYYTgyUNkLSvpANT+b1kl76+mBpUHEd2+a3dg8C/pueNJzvQV7vuruJaAvxA2UjBW0k6tL0wIl4GfkP2S/9ZEfFMhXXsQ5Yg29L8cmC4shFzTwAGAWPTvaxvk13G3OSHmWZvlBOU2WaIiPXAMWT3fp4GlgH/DWyTytcCx5E1XniRrDHC70pW8aX0/BVkN5xvrHbdVca1J1m3M01p26WmkjVq6OzyXmkDifvJGlo8SDZy64fJfgfzBHAf8MOu4jLbHB5R16yHSLoKaIqIb+Qcx2jgMWDHkst2peU/A5ZHxHc2ebJZD/IZlFk/kpqwfwW4toPkNIysKf1fezo2s3JuxWfWT0gaCrxA1ppwkybCkg4jazQxlY1bDJrlwpf4zMyskHyJz8zMCskJyszMCskJyszMCskJyszMCskJyszMCskJyszMCskJyszMCun/Awex9a8XUMngAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.plot(modes.s.imag, np.abs(modes.s.real), 'x')\n", "plt.xlabel('Frequency $j\\omega$')\n", "plt.ylabel('Damping rate $|\\Omega|$')\n", "plt.title('Complex eigenfrequencies of modes')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Having found these modes, we can now use them as a basis for constructing models of our resonator which can describe all the dynamics in a relatively simple form." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exact solution to extinction\n", "\n", "Now we utilise this model to calculate the extinction cross section of each mode. First we setup the incident wave polarisation and direction and the frequency range of interest." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "num_freqs = 101\n", "freqs = np.linspace(1e9, 30e9, num_freqs)\n", "\n", "e_inc = np.array([1.0, 1.0, 0])\n", "k_dir = np.array([0, 0, 1])\n", "plane_wave = PlaneWaveSource(e_inc, k_dir, p_inc =1.0)\n", "\n", "area = np.pi*(0.5*mesh.max_distance)**2\n", "\n", "extinction = np.empty(num_freqs, np.complex128)\n", "extinction_modes = np.empty((num_freqs, len(modes)), np.complex128)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For comparison purposes we calculate the extinction directly, and time how long this calculation takes. Note that if you increase the mesh density or model a more complex structure, the direct calculation can become very slow, since a large impedance matrix must be filled and solved at every frequency." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "60a230d282af42c39b76e11541763fba", "version_major": 2, "version_minor": 0 }, "text/plain": [ "HBox(children=(Label(value='Frequency Sweep'), FloatProgress(value=0.0)))" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "2.85 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n" ] } ], "source": [ "%%timeit -n 1 -r 1\n", "\n", "for freq_count, s in sim.iter_freqs(freqs):\n", " Z = sim.impedance(s)\n", " V = sim.source_vector(plane_wave, s)\n", " extinction[freq_count] = np.vdot(V, Z.solve(V))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution based on modes\n", "\n", "Knowing the location of the mode frequencies (or poles) $s_n$, together with the the mode current distributions $\\mathbf{v}_n$, these form a compact but accurate model of the currents $\\mathbf{j}$ excited on the meta-atom by an incident field $\\mathbf{E}$. \n", "\n", "$$\\mathbf{j}(\\mathbf{r}, s) = \\sum_{n} \\left(\\frac{1}{s - s_n} + \\frac{1}{s_n}\\right)\\mathbf{v}_{n}(\\mathbf{r}) \\frac{\\int \\mathbf{v}_{n}(\\mathbf{r}')\\cdot\\mathbf{E}(\\mathbf{r}', s) \\mathrm{d} \\mathbf{r}'}\n", "{ \\int \\mathbf{v}_{n}(\\mathbf{r}')\\cdot \\frac{d}{ds}\\mathbf{Z}(s)\\big|_{s=s_n} \\cdot \\mathbf{v}_{n}(\\mathbf{r}') \\mathrm{d} \\mathbf{r}'} $$\n", "\n", "Note that the normalisation integration in the denominator to the derivative of the impedance matrix $\\mathbf{Z}$ has already been applied to the mode currents. Now we calculate the extinction based on this model" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "eb173388b363450e8e35f0ae771aa525", "version_major": 2, "version_minor": 0 }, "text/plain": [ "HBox(children=(Label(value='Frequency Sweep'), FloatProgress(value=0.0)))" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "161 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n" ] } ], "source": [ "%%timeit -n 1 -r 1\n", "for freq_count, s in sim.iter_freqs(freqs):\n", " V = sim.source_vector(plane_wave, s) \n", " I_modes = (1/(s - modes.s) + 1/modes.s)*modes.vl.dot(V)\n", " extinction_modes[freq_count] = V.vdot(modes.vr*I_modes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Timings will vary depending on the exact structure to be solved and the performance of the computer. However in all cases, once the scalar model has been constructed, the system can be solved for different frequencies, input polarisation etc with very low computational cost. This becomes much more important for complex structures with many mesh elements.\n", "\n", "We can now see that the extinction cross-section contribution from each of the modes." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXwU9f348dd7jyQcgQQI932qKIcigogiHngCVvGsZ5V6VW1/9Wi1te23VqtWW63WelUsXiiKR0GtB95FQUU5BbnvJBBICMlen98fn9mwiUlIQrI72Xk/H4957O7M7Mx7dpN57+eYz4gxBqWUUsptfKkOQCmllKqOJiillFKupAlKKaWUK2mCUkop5UqaoJRSSrmSJiillFKupAlKNYiI/FpEHm+C7Y4VkeWNvd067LeTiHwoIsUi8hex/iUiO0Tk87rGJSIXiMjbyYi5MYjIVSKyVURKRKR9kvY5V0QuT8a+9kVEjIj0T3Ucqnqi10F5h4isAToB0YTZTxljrt3H+8YB040x3ZsgJgMMMMasbOxt1zOO3wDDgTONMUZExgLPAYOMMbtTGVtTEZEgsAsYZYxZmMT9zsX+PTX6D5wGxOKKvz9VvUCqA1BJd7ox5p1UB+FCvYAlZu8vtl7AmnRNTo5OQBawONWBKFUdreJTAIjIP0TkpYTXfxaRd0WkFTAH6OpUA5WISFcR+Z2ITHfW7e1UlVwsIutEpEBEbk3Ylt+pEvzeqUJbICI9RORDZ5WFznbPEZFxIrIh4b0HOlVCRSKyWEQmJix7SkQeEpH/ONudJyL9ajnGUSLyqbOthU7JEBF5CrgYuMmJ46fA48Bo5/Xvq4mrh4i8LCL5IlIoIn935l8iIh8nrHeAiPxXRLaLyHIRObuu8YvI4IT3bnU+w84iUppYHScihzlxBKs55kwR+auIbHKmvzrzBgLxKssiEXmvPp+Zs+xSEVnqxL7K+dwS3ztJRL4WkV3Od39SwuJeIvKJ8963RaRDDfsfJyIbROQmEdkmIptFZLKInCIi3zmfza/3dbwJy290trFJRC6r5rO61/kb3ioij4hIi+riUklijNHJIxOwBji+hmUtge+AS4CxQAHQ3Vk2DthQZf3fYatpAHoDBngMaAEMBcqBA53lNwLfAoMAcZa3d5YZoH/Cdiv2BQSBlcCvgQxgPFCMrXYDeArYDozE1gY8Azxfw/F1AwqBU7A/zE5wXuclbOuPCetfAnxcQ1x+YCFwP9AKWwo5qur7nGXrgUud+A51PtfB+4ofyAY2A//P2X42cISzbDZwVUJs9wMP1nDcfwD+B3QE8oBPgf+r8r0FGviZnQr0c77TY4BS4FBn2Uhgp/Men7OtA5xlc4HvgYHYv5e5wF01xDAOiAC/df4ergDygWedz2QwUAb0rcPxngRsBQ52vptnSfj7A/4KvAa0c7b9OnBnqv9vvTylPACdkvhl2wRVAhQlTFckLB/pnDDXAuclzB9H3RJU94TlnwPnOs+XA5NqiKm2BDUW2AL4EpY/B/zOef4U8HjCslOAZTXs52bg31XmvQVcnLCtuiao0c5J8gcndionqHOAj6os/ydw+77iB84DvqrhWM4BPnGe+53PaGQN634PnJLwegK26jLxe6spQdX6mVWz/izg+oTjvL+G9eYCtyW8vhp4s4Z1xwF7AL/zOtuJ+YiEdRYAk+twvE+SkAixCdIA/bFJdjfQL2H5aGB1U/5P6lT7pG1Q3jPZ1NAGZYz5XERWYX99zmjAtrckPC8FWjvPe2BPHPXVFVhvjIklzFuL/TW+r31W1QuYIiKnJ8wLAu83IK4ewFpjTGQf6/UCjhCRooR5AeDfCa8b8pm9CjwiIn2xJ9mdxpjPa1i3K/Yzi1vrzKuLWj8zETkZuN2JwYcthX+bEP/sWrZd1+8NoNAYE+/Ys8d53JqwfE/C+2s73q7YZJa4LC7PiX+BiMTnCfYHgEoRTVCqgohcA2QCm4CbgDudRfvb1XM9tipoUT3ftwnoISK+hCTVE1sV2ZAY/m2MuaIB761uWz1FJLCPJLUe+MAYc0ID93FedQuMMWUiMgO4ADiAygmvqk3YRBPvCNHTmVfXGKr9zJx2nZnARcCrxpiwiMzCntTj762xPbAJ1Xa8m7GJk4RlcQXYRDfYGLOxqYNUdaOdJBQATqP5H4EfAxdiOwwMcxZvBdqLSNsGbv5x4P9EZIBYQxIa+bcCfWt43zxstctNIhJ0GuhPB55vQAzTgdNFZILYThtZTgN8Q7rOf4492d0lIq2cbY2pZr03gIEicqETf1BEDheRA+uwjzeAziJyg9N4ny0iRyQsfxpbnTjRObaaPAfcJiJ5TkeE3+5j/US1fWYZ2B8z+UDEKU2dmPDeJ4BLReQ4EfGJSDcROaCO+90ftR3vDOASETlIRFpiS38AOD+AHgPuF5GOAE7ME5IQs6qBJijveV329sYrEZFXRCSA/Sf+szFmoTFmBbZjwr9FJNMYswz7j7/K6c1V1yqiuPuwJ4e3sdfdPIFtHAfbljXN2e7ZiW8yxoSwJ+CTsb9wHwYucuKpF2PMemCSc1z52F/4N9KA/wGnuul0bNvFOmADtl2o6nrF2JP2udhf8VuAP2NP7PvaRzG2g8HpzvtWAMcmLP8EiAFfGmPW1LKpPwLzgW+w1W9fOvP2qbbPzInvOuz3ugM4H9vBIP7ez7GdQ+7Hdpb4AFuyaWo1Hq8xZg62I8R72M43VXsu3uzM/5+I7ALewXbsUSmiF+oq1Uw5XcOfNS644FWppqAJSqlmSEQOB/4L9HBKM0qlHa3iU6qZEZFp2OqnGzQ5qXSmJSillFKupCUopZRSruT666A6dOhgevfuneowlFJK7YcFCxYUGGPy6vMe1yeo3r17M3/+/FSHoZRSaj+IyNp9r1WZVvEppZRyJU1QSimlXEkTlFJKKVdKeoISkRwReUlEljk3Oxud7BiUUkq5Xyo6SfwNe++Xs0QkAzvEvVJKKVVJUhOUiLQBjsaOwhwfDDSUzBiUUko1D8mu4uuLHRX5XyLylYg8LiKtqq4kIlNFZL6IzM/Pz09yiEoppdwg2QkqABwK/MMYMxx7r59bqq5kjHnUGDPCGDMiL69e13Wpeli/vTTVISilVI2SnaA2ABuMMfOc1y9hE5ZKslAkxvH3fZDqMJRSqkZJTVDGmC3AehGJ3wTsOGBJMmNQVigaozwSIxrTwYKVUu6Uil58PwOecXrwrcLedVMlWSgSq3hskeFPcTRKKfVDSU9QxpivgRHJ3q+qLBx1ElQ0Rgs0QSml3EdHkvCoeAkqnqiUUsptNEF5VCiqCUop5W6aoDyqogQV0U4SSil30gTlUYltUEop5UaaoDxK26CUUm6nCcqjEruZK6WUG2mC8ijtJKGUcjtNUB5VUYLSBKWUcilNUB4VjppKj0op5TaaoDwqFI0CENY2KKWUS2mC8qj49U/aBqWUcitNUB5VrtdBKaVcThOUR4UrroPSNiillDtpgvKoeMlJr4NSSrmVJiiP0pEklFJupwnKo8J6oa5SyuU0QXlUKBJDRDtJKKXcSxOUR4WiMVoG/Xq7DaWUa2mC8qhQJEarzIBW8SmlXEsTlEeFozFaa4JSSrmYJiiPipegyrWbuVLKpTRBeVQ4amiV6dcSlFLKtTRBeVR5JEarDK3iU0q5lyYojwpF450ktBefUsqdAsneoYisAYqBKBAxxoxIdgzKjsXXKtPPrrJIqkNRSqlqJT1BOY41xhSkaN8KpwSVEaCwJJTqUJRSqlpaxedR4aheB6WUcrdUJCgDvC0iC0RkanUriMhUEZkvIvPz8/OTHJ43hCL2Oigd6kgp5VapSFBjjDGHAicD14jI0VVXMMY8aowZYYwZkZeXl/wIPaCik4QOdaSUcqmkJyhjzCbncRvwCjAy2TGo+IW6fi1BKaVcK6kJSkRaiUh2/DlwIrAomTEoK6TXQSmlXC7Zvfg6Aa+ISHzfzxpj3kxyDArbSaKljiShlHKxpCYoY8wqYGgy96mqF+8koRfqKqXcSruZe5Qdiy9Qcet3pZRyG01QHmSMqbhQVztJKKXcShOUB4WjhqBfyAz4tA1KKeVamqA8KBSNEfT7CAZ8hLWKTynlUpqgPCgciZER8BH0i3aSUEq5liYoD6ooQfl8hKIxjNEkpZRyH01QHhSKxMjw+/D5hIBPiMQ0QSml3EcTlAeForaKDyDo144SSil30gTlQfESFEBGwKcDxiqlXEkTlAeFq5SgyqPRFEeklFI/pAnKg0KRGEG/AJChPfmUUi6lCcqDKrVB6bVQSimX0gTlQbYEpZ0klFLupgnKg8JRQ2ZCG5SOx6eUciNNUB6UWILSNiillFtpgvKgUDRa0QaVoQPGKqVcShOUB4UjplIblN4TSinlRpqgPKi8ynVQ2gallHIjTVAeFE4YSSLo127mSil30gTlQYnXQWUEtJOEUsqdNEF50A9KUFrFp5RyIU1QHhS/HxRoG5RSyr00QXlQ5So+LUEppdxJE5QHVR4sVruZK6XcKSUJSkT8IvKViLyRiv17XSgSSxjqSLQEpZRypVSVoK4HlqZo355X9X5Q2otPKeVGSU9QItIdOBV4PNn7VlbV0cy1ik8p5UapKEH9FbgJqPGsKCJTRWS+iMzPz89PXmQeEY4a7SShlHK9pCYoETkN2GaMWVDbesaYR40xI4wxI/Ly8pIUnXeUVypBaRuUUsqdkl2CGgNMFJE1wPPAeBGZnuQYPC+xDSpD26CUUi6V1ARljPmVMaa7MaY3cC7wnjHmx8mMQdk2qIqRJAJ6oa5Syp30OigPqtqLTztJKKXcKJCqHRtj5gJzU7V/L0sc6ihDx+JTSrmUlqA8KKSDxSqlmgFNUB4UqlTFJ4Qi2klCKeU+mqA8qGonCS1BKaXcSBOUByV2ksjUKj6llEtpgvKgxNHMtQSllHIrTVAelDjUkXYzV0q5lSYoDwpVGeoopCNJKKVcSBOUxxhjbC8+vQ5KKeVymqA8Jhw1BP2Cz+e0QWmCUkq5lCYoj0kcRQKcThLaBqWUciFNUB4TjuztYg62ik/boJRSbqQJymOqlqC0DUop5VaaoDwmcRQJgGBAtJu5UsqVNEF5TOI4fKCdJJRS7qUJymPC0colqIBPiMQMsZi2Qyml3EUTlMeEqnSSEBHbDhXTUpRSyl00QXlMOLp3HL64oF8Ia08+pZTLaILymPIqJSiADL0WSinlQpqgPCZxHL447SihlHKjBiUoEflxYweikiMcNWQGfpigyrUEpZRymXolKBHpKSK9gcuqzL+vEWNSTai6ElSG3hNKKeVCgXqu3xebnIaLyLvAcmc6obEDU00jHP1hG5R2klBKuVG9EpQxZi4wV0ROB94A+gMHAWc3fmiqKWgblFKquWhoJ4mTAYwxK4wxrwJrGy8k1ZSqjiQBzl11NUEppVymoQnqTeAlEekiIjcC79blTSKSJSKfi8hCEVksIr9v4P5VA1Udiw+0m7lSyp0amqDeAgqAdYAAR9XxfeXAeGPMUGAYcJKIjGpgDKoBqmuDsiOaaxuUUspdGpqgPgDmAUOxHSQG1eVNxipxXgadSc+MSVRdCcp2ktASlFLKXerbiy9urDEmDCAi5wLTgNPq8kYR8QMLsB0sHjLGzKtmnanAVICePXs2MERVnar3gwK9Dkop5U4NLUEdKSJviMh0YAdwdV3faIyJGmOGAd2BkSJycDXrPGqMGWGMGZGXl9fAEFV1qu0koddBKaVcqKEJ6h7gPKCrMSYGPFrfDRhjioC5wEkNjEE1gO1mXnmwWL2rrlLKjRqaoIqNMcUJrzPq8iYRyRORHOd5C+B4YFkDY1ANEI7GqhnqSNuglFLu09A2qH+JyDNAjogcie3sUBddgGlOO5QPmGGMeaOBMagGqHo/KLDdzEPai08p5TINSlDGmOki8hVwhjNdto+3xN/3DTC8IftUjSMcNdWPJKGdJJRSLrPPBCUilwLXAL2BDcB7wL3GmMXA4iaNTjW6aktQ2gallHKhWtugnO7edwEvAlcC/wIOBZaLyISmD081tpq6mYe0BKWUcpl9laCuA84wxnyaMO9vInIC8IyIHA1sAQ42xnzcVEGqxlNdCUoHi1VKudG+evH1qJKcADDG/Bf4BTAdW803uAliU02g2pEkAqKdJJRSrrOvBJUvIj1qWDYD2+HhemPMPxs3LNVUah6LT0tQSil32VeCeg74Uw3L/ECBMealxg1JNaXq2qD0jrpKKTfaV4K6ExgsIrNFZGiVZbdiB4xVzUj1g8VqglJKuU+tnSSMMaUiMh54AFggIuuwXc17Oasc18TxqUZW4w0LI9oGpZRyl31eB+WMmXeRiNyGHZooD1gN/McYs7uJ41ONLByt/nYbekddpZTb1HkkCWPMOuDJJoxFJUGNF+rqdVDNU0k+bP4aCr4D8UMgA/yZ0LYbdDsMMrNTHaFSDdbQsfhUM2WHOqo8mrm2QTUzmxfC/x6B1R9CqBi6DIO8AwADkXI77VgNW76F9v2gxyg4+EfQczSI7HPzSrmFJiiPqfZC3YBPq/jczhhY+Q58+gAUrIQjfgrH3Ai5fWpOOpFy2PwNrPkIXr8BYhEY/mMYdgFkd0pu/Eo1gCYoj7H3g9LroJqVkm3w6jWwcwOMuQEGn2Gr8vYlkAk9DrfTUT+HDfPhq6fhoZEw9Dw7TxOVcrGG3g9KNUPGGNuL7wfXQQlhHUnCnb57Gx4ZC50PgZ9+CEPPqVtyqkrEJqqJD8I1n9vXDx8Bb90KpdsbP26lGoEmKA+Jtz/5fNoG5XqxKMy5Bf7zCzjrSTjut+Cv623X9iG7E5x0J1z1KYT32BLV/CftPpVyEU1QHhKuZhQJ0NHMXScagVd+CtsWw5UfQe8xTbOfNl3htPvgwlfgmxnw2HhbDaiUS2iC8pDqOkiAk6C0BOUOkRC8dAnsKYLzZ0CL3KbfZ+dD4NI5MOpqeP4CmHMzhPQSR5V6mqA8pKYSlHaScIlwGbxwge2xd+4zEGyRvH2L2Patqz+zyfHh0bBqbvL2r1Q1NEF5SHk14/CBvd1GWIc6Sq1YDGb+BDJawZSnbA+8VGjZDn70TzjlXph1Dbx+PZSXpCYW5XmaoDykulttgJagXOHd39vedGf8s/E6Q+yPgSfC1Z9CNAyPjIG1n6U6IuVBmqA8pLou5qAX6qbcV8/AkllwzvTUlZyqk9UWJj8ME/4EL14M//2tvfhXqSTRBOUhNXWS0BJUCq35xJ74z58BrdqnOprqHXAqXPkJFH4Pjx0HW5ekOiLlEZqgPMR2kvjhsDj2Oihtg0q6nRvgxUvgzMcgb1Cqo6ld6zxbwjtiKkw7DT572LabKdWENEF5SHkNJSi/TzDGEI1pkkqaaARmXg6jr4Z+41MdTd2IwKEXweXvwOKXYfoZsHNjqqNSaSypCUpEeojI+yKyVEQWi8j1ydy/19mRJKr/ynU0iST78G4IZMGRzfBfoF1fuPRN6DUGHj0GFr2c6ohUmkp2CSoC/D9jzIHAKOAaETkoyTF4VigSI7OaEhTYdijtKJEkqz+CBdNsjz1fM63E8AfgmJvg/Bfg/Tvg5alQtjPVUak0k9T/DmPMZmPMl87zYmAp0C2ZMXhZTRfqgu3JpzctTILdhXYYo8kPpcdI4t0Os4PYZrSGf4zRi3tVo0rZzzcR6Q0MB+ZVs2yqiMwXkfn5+fnJDi1t1dSLD+I9+bQNqkkZA69fZ28e2P/4VEfTeDJa2TH9Tv8rzLraGSqpNNVRqTSQkgQlIq2BmcANxphdVZcbYx41xowwxozIy8tLfoBpKlRrCUq0DaqpffsibF8F43+T6kiaRv/j4apPoLQQ/jkW1n+e6ohUM5f0BCUiQWxyesYYo62rSVRbCUoHjG1ixVvgzV/ZC1/ddDFuY2uRC2c+bm8P8sKP4e3b7C09lGqAZPfiE+AJYKkx5r5k7ls5CaqGElSG3nKj6Rhjb7k+4lLoOjzV0STHQZPs/aZ2brA3XFz3g5p8pfYp2SWoMcCFwHgR+dqZTklyDJ5V01h8oN3Mm9Q3L0DROjj6plRHklytOtiBb4/7Dcy4CP7zSygvTnVUqhlJdi++j40xYowZYowZ5kyzkxmDl9VWggr6tQ2qSezabG+rPvnhht2qPR0cNMnexiOyBx4aBd+9leqIVDPRTC/CUA1Razdzv4+Q3nKj8c3+JRx2CXQdlupIUqtlO5j0kO1eP+dmeOFCHYVC7ZMmKA8pr6WKLyOgVXyNbslrkL8cjr4x1ZG4R99xtjTV8UB45Cj49O/2lh5KVUMTlIeEI6bawWJBRzRvdHuKYM5NMPEBCGalOhp3CbaAY38NP/kvrPwv/PNoWPVBqqNSLqQJykNC0WiNQx1pJ4lG9t/fwqCTodeRqY7EvTr0hwtnwbhb4NVrbbXfjrWpjkq5iCYoDykLx8gM+qtdFgz4KNdu5o1j9Uew8h04/nepjsT9RGwnims/h04H28Fn3/k9lP3g+n3lQZqgPKSwpJz2rarvSWZ78Wknif0WLoPXr4dT7rF3pFV1E2wB426GKz+2FzU/eBh8/pi2T3lcINUBqOQp3B2ifevqRzHQNqhG8uHd0GmwvQutqr+23eGMf8Dmb+woFPMegWNvhYMmN+nI7yWhEpbvWM6y7cvYWLKRneU7KSovYlf5LnziIyuQRaY/kzYZbeiR3YOebXrSM7sn/XP7k+lP45FBUkwTlIcUFJfToXVNJShNUPttyyJY8JQdQUHtny5D4KJX4fv34L3/g4/usxf8DjjRVgvup/JoOV9s+YIP1n/Ap5s+JX9PPgNyBjCo3SB6ZvdkYO5AcjJzaJPRhqiJUh4tpzxSzs7QTtYXr+edte+wdtda1hWvY2DuQA7teCiHdTqMI7ocQVZAO8U0Fk1QHmGMoWB3iA41laACOtTRfolF7Ujl438D2Z1THU16EIH+x9k7Di97w3Y8+eBuOOZmGHBCvRNVzMSYt3keM1fM5JONnzAwdyBHdz+a+8bdR7+cfgR89T8dloZLWVSwiC+3fcm0JdO45aNbGNVlFON7jmdcj3FkZ2TXe5tqL01QHlFcHiHD7yOrpk4SeruN/fPF4+DPhEMvTnUk6UcEDjwdBp0CS16Fd263N0k85mYYeNI+q/6Kyop4acVLzPxuJq2CrZgycAq3HnEruVm5+x1ay2BLRnYZycguI7ly6JUUlRXxwYYPeHvt29w5707Gdh/LGQPOYGTnkfhEm/zrSxOURxSWhGhfQ/UeQIYOddRwOzfA3Lvgsrea7x1ymwOf395L66DJtkQ190/w7h9gzHVw8Fk/GEqqcE8h05ZM4+UVLzOu+zjuPvpuDu5wMNIIVYQ1ycnKYVL/SUzqP4kdZTuYvXo2935xLyXhEs4edDZnDjiTtpnaeaauNEF5REFJeY3Ve2BLUNrNvAHiI5WPugryBqY6Gm/w+eCgibZU9f178Mlf4b07YNSVcOhF7BR47JvHeGXlK5zc52RePO1FurTukvQwc7NyueDAC7jgwAtYXLCYZ5c9y8kvn8yE3hO48MAL6ZvTN+kxNTeaoDyiti7mYK+DKimPJDGiNPHNDCjeDGNuSHUk3hNvo+p/HGxcQPjTv/PC/Ad4LKcNx/ccz8sTX6ZTq06pjhKAwR0Gc8dRd1Cwp4AXl7/IZW9dxpC8IfzkkJ8wNG9oqsNzLU1QHpFfEqJDtpagGlVJPrx9K5w/w7sjlbvEh2Y3dwfy6TFwLE+YHPp/PhM2rIIRP4GBE2z1oAt0aNGBq4ZdxaUHX8qslbO4+cOb6dKqC1cMuYLRXUY3afVjc6QJyiMKS8rpUEsJqm2LIN+WhpIYURqYcyMMOx+6HZrqSDxry+4t3P3F3SzfvpxfHfErjup2lF1w/B2w+BX46C8w+0Y7ovyw86Ftt5TGG5cVyOLcA87lrIFnMWf1HO76/C7aZLThqqFXcWTXIzVRObRF1yMKSsprLUF1y2nBpqKyJEbUzC19w15MOu5XqY7Ek6KxKNOXTGfK61Pol9OPlye9vDc5gR2ZYtj5cMW7cO4zsGsj/ONImH6WHWU+4o4fYwFfgNP7nc4rE1/hggMv4J4v7uHHs3/Mxxs/xhjtVaslKI8oLAnRvm/tCWpj0Z4kRtSMlW63v8rPfNyeCFVSfV/0Pb/95LcE/UGmnTyNvm330dmg6zDo+leY8CfbTX3eI/DGDTD4RzD0PFsCTnGJxe/zc3If24Hi7TVvc88X9/BIxiNcM+waRnUZ5dkSlSYoj7C9+Gqu4uvcNottxWVEojECNdzUUDlm32h7kfUek+pIPCUcCzNt8TSeXvw01wy7himDptTv2qKMljDsPDvtWGM7uLx8OYgPDj7TTnmDmiz+uvCJj5P6nMQJvU7gzTVv8qd5f6JdVjuuHX4th3c+PKWxpYImKI+w10HVXILKCPho1yqDrcXldMvRUkGNFr8Cm7+Gn36U6kg8ZeWOldz6ya3kZObw/GnP07V11/3bYG5vOOYmezPJjQtg0cvw9CRo0Q4Gn2F/gKQwWfl9fk7teyoTek9g9urZ3P7p7XRt1ZVrhl/D8I7DUxZXsmmC8oj8knLyaklQ4FTz7dijCaomxVts6em85+2vcdXkorEo/17yb55c9CTXHXodZw44s3Gru0Sg+wg7nfhHWPeZrQZ8ejJkZttEdcCp0GVYSqoBA74AE/tN5OQ+J/PG929wy4e30Lttb64edrUnuqdrgvKA8kiUsnCUNi1q/7q75bZkk7ZDVc8YeO06O5RR9xGpjsYT1u9az22f3IZPfDx76rN0z+7etDv0+Wy1be8xcNJdtmS19DWYeQWESuywSoNOht5jk/4DJegLcsaAMzit72nM+n4Wv/zgl/TP6c/VQ6/mkLxDkhpLMmmC8oDCkhDtW2Xu85dn15ws7ShRk6/+DcWb4JzpqY4k7RljePG7F/n7V3/niiFXcMGBFyR/HDufD3ocbqcT/w8KVsDyOfDJ3+Cly6DnKOh/AvQ/Htr3S1rpKugPMmXgFCb1m8QrK17h53N/Tv/c/lw19Kq0LFFpgvKAgpLyWsfhi+ue04Ilm4uTEFEzk/8dvPM7uPgNvSC3iW3dvZXbP7udorIinjrpKfcMB9RhgFC2F7QAABsVSURBVJ3GXAdlO2HVXFjxNnz6gO1k0Xcc9DsWeh8NrfOaPJwMfwbnHHAOZww4g1krZ3HjBzfSp20fpg6ZymGdDmvy/SeLJigPKCyp+TYbibrltuCdpduSEFEzEi6Dly61t9HodFCqo0lbxhhmr57N3V/czbmDzuXyIZcT9AVTHVb1stra29QfNMlW/RassGMCfjMDXv+5vRi4z9HQ+yjoORpadWiyUDL8GZw96GzO6H8Gr37/Krd9fBsdW3Zk6pCpaXHBb1ITlIg8CZwGbDPGHJzMfXtZfh1LUF31Wqgfevs2aN/fjkSgmsT2su388X9/ZFXRKh4+/mEGtx+c6pDqTsQOEpw30A5WG43AloWw+kNYMA1mXQ1tutpE1XM09DwCcno1epVg0B/krIFnMbn/ZN5c8yb3fHEPGf4MLjvkMk7oeQJ+lwz1VF/JLkE9BfwdeDrJ+/W0wpLQPnvwwd5efMaYZv/Lq1EsfR1WvGW7lOvn0STeXfsuf5z3R07rexp3jr2z+d8+3R+AbofZ6aif24S19VtY+6lz08XfAAI9Rjq9Bw+3PQQbqdNFwBfgtL6ncUqfU5i7fi5PLnqSB758gEsGX8LEfhOb3d1+k5qgjDEfikjvZO5T2TaoTm32/Y+fnRUk6BeKSsPk1jJunyfsWAtv/BzOfQ5a5KQ6mrSzo2wHd867kyXbl/CXY/7CoZ3SdDxDfwC6DrfT6GtslWDRWlj/OWyYD4tnwbaltpTebTh0PdSObJF34H61d/rEx/ie4xnfczxfbfuKJ799koe+foizBp7FeQecR4cWTVft2Jhc2QYlIlOBqQA9e/ZMcTTNX2FJOYO7tqnTuvFqPk8nqFApvHABHPUL24tLNap31r7DHfPu4JQ+p/DimBdpEfDQdXci9iLh3N4w5Gw7L1wGWxfBpq9s4pr3iP2BlDcQugyFzkOg8yHQabC9NquehncczoPHPcjqnat5ZukzTJw1kWN7HMv5B57v+upUVyYoY8yjwKMAI0aM0BET91PBPkaRSNQ9twUbduzh4G4eveunMfDaz6DjQfYmhKrRbCvdxp3z7mRl0UruH3c/wzoOS3VI7hDM2nuxcFyoFLYutqOWbPkWFj5nS1qtO0Kng+3fZ6eDoONgaNfXltT2oU/bPtw26jZ+NvxnvPTdS/z8/Z+T1yKPcw84lwm9J5Dhd9+PUlcmKNW49jUOXyI7qrmHO0p8+iAUrrC3b9d2p0YRMzFmrpjJg18+yJRBU7jr6Luaf1tTU8toufc6rLhYFApX2sS1bQl88yJs+50d4aRdP+h4AOQd4HSJH2Svzwr88HNum9mWnxzyEy4ZfAkfbPiA55Y9x73z72Viv4mcNfAserXplbzj3AdNUB5QUMdu5uDxnnzfvwef/R0uf1dHKW8k3+34jjv+dwfhWJjHJzzOwNyBqQ6p+fL57fiAeYOAH+2dHyqFgu8gfxnkL4dvX7KPRetsD8IOA6D9AJuw2vez7V3ZXfH7/BXtVGt3rWXmiplcNOci+uf0Z3L/yRzf6/iUV78mu5v5c8A4oIOIbABuN8Y8kcwYvCYWMxSVhmhXxzalbrkt+Hp9URNH5UJbl9ghbc6eBjk9Uh1Ns1caLuUfC//Bqytf5Zph13DWwLOabVdn18to6dxSpEqVaSRkR20vXGGv1dr0FSyaCYXf24uNc3vZ6sHcPvRq14df5B7Cz449ifeKV/PK6te56/O7OKHXCUzqP4lhecNS0rM32b34zkvm/hTsKA3ROitAsI630PBkCWrnBnhmCpx0p724UjWYMYY317zJfQvuY2Tnkbw86eVm02Ms7QQy9l6jVVV5iU1e21fBjtW2k8ay/xDcsYYJuzYyoWUHtuZ25/W1X3D76rcJi49TOh/BqX1Pp2+3UbbdLBmHkJS9qJQp3F336j2wwx1t3OGhBLWnyN5l9Yipe3tVqQZZXLiYP3/+Z8oiZdw19q60GnIn7WS2hs4H26mqaASKN9GpaB2X71jLT3asZemOZfxn3WdcvvZt2kfCTAjBhGAHerTpBW262dEz2nSzVYrZXexjNe1f9aUJKs0VFJfTvh5dxju0zqS4PEJZOEpWMM2rZMJl8PwF0PcYOPK6VEfTbG3ZvYWHvn6Ijzd+zM+G/4xJ/SZpdV5z5g9ATk879T4KAQ5ypl/Eony5dQFvrXyVCzfMpaO/gONjmRxXVEbfDQuQ4s1QvNl23MhqA9ldIbuznRpAE1SaK6hnCcrnE7q0taOa98tr3YSRpVikHF68xI6TNuFP2mOvAXaW7+TJRU8yc8VMpgycwuuTX6d1Rhr/zSj8Pj+HdxnJ4V1G8qtYlAVbF/Duune5cv17ZGVkceywkzm6+9EM6zCEwJ4dULLFJqtdmxq0P01Qaa6guO5dzOPiQx6lbYKKlMOMi22vqB89Zh9VnZWGS3l22bM8vfhpxvccz8zTZ9KpVadUh6WSzO/zM7LLSEZ2GcktI29hSeES3l//Pnd/cTebdm/iyK5HMrbbWEb3GE2HFhOAy+q9D01Qaa5wd3mdL9KNS+troRKT01n/0ttn1ENpuJTnlz/PtMXTGNl5pLtuh6FSSkQY3GEwgzsM5trh17Jl9xY+2vgR7617jzvn3dngm01qgkpzBcUhhvao31hyaduTL7zHVuv5Apqc6mFn+U5eWP4Czy59lhGdR/DEiU/QP7d/qsNSLta5VWemDJzClIFTCMfCfJv/LS/yYr23owkqzdkSVD2r+HJb8L/vC5soohQp3Q7PnmMbfif/Q5NTHWzdvZXpS6fzyspXGNd9HE9MeIJ+Of1SHZZqZoK+YIMHA9YEleby6zGKRFz3dCtB7VgD08+EA06D4263t/NWNVqYv5BnljzDJ5s+YWK/ibx0+kt0btWwXlhK7Q9NUGmusB7j8MWlVRXfxgXw3Plw9C9h5BWpjsa19kT28Naat5ixfAY7ynZw/oHn85vRvyE7o/6jZyvVWDRBpTFjjDNQbP1KUF1ysti2q5xozOD3NdPu18bAl0/Du3+AiQ/AAaemOiJXWr59OTNXzGT26tkMzRvK1CFTGdttrF7HpFxBE1Qa2x2KAtAyo34nm8yAn7Ytg2zdVUbXnGY4aGqoFGb/EjZ+CZe9aQfLVBUK9hQwe9VsXvv+NXaGdjK5/2RePO1FurTukurQlKpEE1QaK3RKTw0Z5HFkn3bMWbSFnxzVpwkia0LblsHMy+2tB654FzJapToiV9hZvpN3173LW2ve4tv8bzm257HcePiNHN75cHyibXLKnTRBpbFvNuykT4eGnaAvP6oP1z77FReP7kWgjgPNplQsam+V8cnfYPxv4LBLPD86RMGeAuaun8u7697lq21fMbrLaM4YcAb3j7uflsGWqQ5PqX3SBJXGpv9vLReN7t2g9w7vmUu3nBbMXrSFiUO7Nm5gja1gBcy62g5OecV79nbaHmSMYfmO5Xy88WPmrp/LqqJVjOk2hon9JnLvMffSKqilSdW8aIJKU99tLWZ1wW5OHNzwIWiuOLovf3v3O04f0iUl94LZp/IS+OheWDANxt0Ch1/huS7kBXsK+Hzz53y2+TM+2fgJWYEsjup2FFcOvZKRnUe68jbeStWVJqg0Nf1/azl3ZM863weqOscd0JE75yzls1WFHNnPRff0MQa+mQHv/A76HA1XfQptvNHAX7CngC+3fsmX275k3uZ5bC3dyohOIziiyxFcfsjlrrpdt1L7SxNUGtpdHuHVrzfx5g1j92s7Pp9wxdi+PPbhKnckKGNg+WyYe6cdrujsadBjZKqjajLRWJSVRSv5puAbFm5byML8hRSWFTK843AO7XgofzjyDxzY/kACPv03VulJ/7LT0KyvNzKqbzu6tN3/LuJnDO/GX97+ju+2FjOwU4ou2ozF4Ls5MPcuwMC4X8GgU9KqE0Q0FmXtrrUs2b6EJYVLWFywmOU7lpPXIo8heUMYmjeUCw+6kAG5A7TXnfIMTVBpxhjDvz9by22nHtQo28sK+rl4dC/+9s4KHjxvOL5kXrhbXgxfPwvzHoHMbDjmJhh0arNuZzLGUFhWyMqilazcsZKVRStZvn053+/8nrwWeQxqN4jB7Qdz1bCrOLDdgbTNbJvqkJVKGU1QaWbB2h2UR2Ic2a99o23zkjG9ufRfX3D9C19z75QhZAaacJQBY2DTl/D1c7DoJdvGNOlh6DmqWZWYQtEQG4o3sHbXWtbuWsvqXatZvdNOMRNjQO4A+uf054B2BzC5/2QG5A7QXnZKVSHGmFTHUKsRI0aY+fPnpzqMZiEWM1w5fQEj+7Tj8rGNe5+esnCUX8z4moKSEI9eeBg5LRu5d9iONbB4Fix8DiJlMPQ8GHYB5PRo3P00kpiJUbCngE0lm+y0exMbijewvng9G4o3ULCngC6tu9Azuye92vSiT9s+FVP7rPbu7BWpVBMSkQXGmBH1eo8mqPSwuzzCL2Z8TWFJiH9dejjZWcFG30csZrhzzlLeW7aN+88ZxpDu9bvPVCXGwNZFsGw2LHvd3hZ60Ckw9FzoOTqlpaVQNETBngLy9+SzrXRbxbS1dCtbdm9hy+4t5Jfmk52RTdfWXe3Uqivds7vbqXV3urTuQtDX+N+BUm5njAFjMBhMzFS8DmZm1jtBaRVfGthYtIfLp81ncNc2PHDe8CargvP5hFtPPYgBHbO5avqXdGmbxWVH9eHEgzrVbbSJovWw5mNY9T58/74dhmjQyXDy3dDjiCa79Xo0FmVnaCdF5UUUlRWxo2wH28u328ey7RTuKaSwrJDCPYUU7CmgNFJK+6z2dGjRgY4tO1ZMR3Y9kk4tO9G5VWc6texEViCrSeJt7kwshjGGWCyGMTH7Omac+bGK5ZUfY3b9mKl1nfgJr9J2TT3XiW/TWZeE53vXq7pdZ72EY6PK8srPnW1TdVv2OTW8D5MQezXrJH6+djtUG3vifvdu01T/3oTlxGKYhH0kxlNpXxio2Cd7l9uV7B+CCD6fD5AG/97UElQzVhqK8OaiLdw1ZxlXjO3L5WP7JK3qKBKN8faSrTz58WrW7yhl7IA8RvVtz6i+7eie2xJCu2HLItj0FayfZ6dIOfQ6EvqNh37H1nnEB2MMeyJ7KAmXUBIqqXjcFd5FSaiE4lAxu0K77GP5LnaFdrGzfCc7QzvZWb6T3eHdZGdkk5OZQ9vMtuRm5tA+ox05mTnkBtuSm5lLbkYOucEccjNzae1vaf/ZYjFisajzGCMWjVac4OzJNEosGj8JRomZxGXVrV/5vZXW/cFyZ5tV5xvzg/fGTzx7t7v35Jz4/srv2fvequvET0KxKutVXqfmBAQg4kN8PnuC8tkTVXyeiNjH+PP4fJ84y/xI/OQWX9dZL749EV/FOuITqPTa/liqcTvxuOL7FpxlzrZI2C5VY5Ufxu3MAxL2L3uPR3x7T9bizE/cLnb/1cUTf5+I2OXx44rHVs26CBWfD/GYK+JOPK69+6s4ThEQ8Dk/Fvd+lgnHHz+miuOIxyKV911Fs6jiE5GTgL8BfuBxY8xdta2vCaqy3eURFq4v4pWvNvLW4i0c2iuXK8b2ZUz/ul2nZGIxotEoJholGo0Qi0YrJhOLEo1EnROvMzkn4ZizrolGK07WsWiEWDjEts3r2LphFaWF66B0K23NDrLYw65ADqWZ7Qi1aE+sRQ6+zBb4iBKJholEQoQjYaLRMJFomGgkQiQaIRqJEItFKscWi+FHCBAggB8/PjsZP34En/HhQ/AZQWIgBsQYMNhfedG9J+dYLArG2JOU3++crPwVJ5b4vPhJzOf37z3ZVqwfP6H68fnjrysv8/kqb+eHz/3Vzq88r/ptxE/UlZZVnPT9lZOCz4dP4vFJlfdUThaV9pGwn4qTZOKJOnG5TyrHpO1rqhoNSVBJreITET/wEHACsAH4QkReM8YsSWYcdWWMsSfhSNSeOKPOiTPhdSwaJRaJEHVO2BXzI85rZ3nFCd056UarzItGIphYlHA4Qll5mLLyELvLQuwpD1FaFmL3nnJK9pQTjoRpFRA6t/BxXZYP39IYSxdFWVxl+xW/3qPRihO0idliOH7nV5RPwOfD+ACfYARM/FHA+AwxMcSIESNKTGJEiRIjRoQoUYkREUNUhIjP2F9eQR8i2eDPJSZ+oiZCtLSA6O5CwjE/EeNDfAH8vgwC/iCBQAZBfzYZLbPIDGaSEcwiK9iCFhl2apnRmlaZLcnKyCAYCBIM+skMBggGAgSDfoLBIBkBP4FggGDATzAYIOh3HgN+fP4APn+V5FPDLzyllLskuw1qJLDSGLMKQESeByYBjZKgXn/wn+xesr6iHjRmnLpYW1+ztz6WeMNdzKkvjdexxpznpmJ+DMA5lxlbfrUvxJ7rESFeBjVOUTkmBkicb+x7nWhiFduriA4EYhi7fTEV2zPOepkCmUEfvgwQEWIGisoMPnwgwb1VJFlOVYf4naK4H5/4sTlI7AT4jHGmGD4TxUcMXyyKLxrBF43ii4XxRUP4ohH8/kx8gQz8gZb4Ai3xZ7TAl5mNP7MN/ozW+HxB/D4/NZ7yExbEDJRHYoQiUcojhnAsRigSI1JmCJfEiMQMkWiMqDFEY4bdZje7osVEYhBz5sWMcb46Q9T5nmMx+/kZY9czznNb62CrQnzxKhHiVRF7l4HN1xXL7LOE9ZxDSUhsPqlycInrxdet+HORxNUq1q3uM5OE7daURyttL3GnVN7uvr4TqWatal9Jdctq2HQDcn9NEfioTw2P1Cm+2vdd+7Lat7831vrFUfkYxdT87pr+ZqrbTuWXVfdRZUO1rVvb0VTZpdTr+9q3ZCeobsD6hNcbgCOqriQiU4GpAD179qzzxjPXF7Ogk/b7qF3iH5Bga1rr2zlhD4T3QHhb44VFQigN6Pzmc6aGiH8i0Qa+XynVNJLaBiUiU4AJxpjLndcXAiONMT+r6T2N0gYVi0EsArEwRMP2efwxPkXDYKLO66qP4b3bSFzHxJzHqs+dElj8ebxni1NCq+jlQpWfMeL8Poo/+vwgfvva57fjz/kCdp4/AP4M8AXBH7S3mvBn2ufBFhDIgmBLCDqPegtvpVQKub4NCltiSrzysjuwqcn36vOBLwPQWw8opVRzkexBzb4ABohIHxHJAM4FXktyDEoppZqBpJagjDEREbkWeAvb2vCkMWZxMmNQSinVPCS9R4ExZjYwO9n7VUop1bw03/sWKKWUSmuaoJRSSrmSJiillFKupAlKKaWUK2mCUkop5Uquv92GiBQDy1MdRwp0AApSHUQKePG4vXjMoMftNYOMMdn1eUNzGLhueX2Hx0gHIjJfj9sbvHjMoMed6jiSTUTqPWadVvEppZRyJU1QSimlXKk5JKhHUx1Aiuhxe4cXjxn0uL2m3sft+k4SSimlvKk5lKCUUkp5kCYopZRSruTaBCUiJ4nIchFZKSK3pDqeZBGRNSLyrYh83ZBumc2FiDwpIttEZFHCvHYi8l8RWeE85qYyxqZQw3H/TkQ2Ot/51yJySipjbAoi0kNE3heRpSKyWESud+an9Xdey3Gn9XcuIlki8rmILHSO+/fO/Hp9365sgxIRP/AdcAL2LrxfAOcZY5akNLAkEJE1wAhjTFpfyCciRwMlwNPGmIOdeXcD240xdzk/SnKNMTenMs7GVsNx/w4oMcbcm8rYmpKIdAG6GGO+FJFsYAEwGbiENP7Oaznus0nj71xEBGhljCkRkSDwMXA98CPq8X27tQQ1ElhpjFlljAkBzwOTUhyTakTGmA+B7VVmTwKmOc+nYf+R00oNx532jDGbjTFfOs+LgaVAN9L8O6/luNOasUqcl0FnMtTz+3ZrguoGrE94vQEPfKkOA7wtIgtEZGqqg0myTsaYzWD/sYGOKY4nma4VkW+cKsC0quaqSkR6A8OBeXjoO69y3JDm37mI+EXka2Ab8F9jTL2/b7cmKKlmnvvqIpvGGGPMocDJwDVOlZBKb/8A+gHDgM3AX1IbTtMRkdbATOAGY8yuVMeTLNUcd9p/58aYqDFmGNAdGCkiB9d3G25NUBuAHgmvuwObUhRLUhljNjmP24BXsNWdXrHVqbOP191vS3E8SWGM2er8M8eAx0jT79xpi5gJPGOMedmZnfbfeXXH7ZXvHMAYUwTMBU6int+3WxPUF8AAEekjIhnAucBrKY6pyYlIK6chFRFpBZwILKr9XWnlNeBi5/nFwKspjCVp4v+wjjNIw+/caTR/AlhqjLkvYVFaf+c1HXe6f+cikiciOc7zFsDxwDLq+X27shcfgNPt8q+AH3jSGHNHikNqciLSF1tqAjvS/LPpetwi8hwwDnvrga3A7cAsYAbQE1gHTDHGpFWHghqOexy2qscAa4Cfxuvp04WIHAV8BHwLxJzZv8a2x6Ttd17LcZ9HGn/nIjIE2wnCjy0IzTDG/EFE2lOP79u1CUoppZS3ubWKTymllMdpglJKKeVKmqCUUkq5kiYopZRSrqQJSimllCtpglKqAURkkIh8JSLFInJdDev8VET+msSY7hORK5O1P6WamiYopRrmJmCuMSbbGPNA1YXOBea3AfckzhOR3zq3kdnt3G5hjoicmLDOGhE5vsq2LhGRj+sQ0z3Arc6+lWr2NEEp1TC9gMW1LJ8ELDPGbEyY95Iz/yIgF+gD/A04tTECci70XAZMbIztKZVqmqCUqicReQ84Fvi7iJSIyMBqVjsZ+CDhPcdj7282yRgzzxgTcqY3jTHX12Pf5zj7jE/lIjI3YZW5NFLCUyrVNEEpVU/GmPHY4WuuNca0NsZ8V81qhwDLE14fD8wzxmzYz32/4OyzNdAVWAU8l7DKUmDo/uxDKbcIpDoApdJUDlCc8LoDsCX+QkTaYZOLAJnGmKyEdWeJSCThdQbwZeLGRcQHPIttB/tnwqJiZ99KNXtaglKqaewAshNeFwIVI1gbY7YbY3KAw4DMKu+dbIzJiU/A1dVs/w5n+1V7EGYDRfsbvFJuoAlKqabxDZDYNvUucLiIdN/fDYvIudjRsM8yxoSrLD4QWLi/+1DKDTRBKdU0ZgPHxF8YY94G3sdW3x3hdDkPAqPqs1ERGQ48iC1l5VezyjHAnIaHrZR7aIJSqmm8DhwgIl0T5v0IeAOYjq2GWw1cgL3TaF1NwnZR/zihJ98cqLgJ3kHY+2op1ezp/aCUaiIiMhU4yBhzQ5L29xfge2PMw8nYn1JNTROUUkopV9IqPqWUUq6kCUoppZQraYJSSinlSpqglFJKuZImKKWUUq6kCUoppZQraYJSSinlSv8fjrqAUgNEVAoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.plot(freqs*1e-9, extinction_modes.real/area);\n", "plt.xlim(0, freqs[-1]*1e-9)\n", "plt.xlabel('f (GHz)')\n", "plt.ylabel('$Q_{ext}$')\n", "plt.title('Extinction efficiency of each mode')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In addition, we obtain the reactive stored enerrgy of each mode, given by $W_{m}-W_{e}$, the difference between magnetic and electric and stored energies. At resonance these terms go through zero, following from the imaginary part of the admittance.\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAEZCAYAAAAg+KppAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydd3hcxdm37+ecbeqyins3NraxHRPcqKa30NNMQgshDpCQ8BJqAoEUAiG8H70ESCgvCSUkASe0gIlpoZlgcMHGDfemLq209cz3x5yVVkJ1LWlX0tzXNdeZMzPPzLNreX475cwRpRQGg8FgMPQlrHQ7YDAYDAZDVzHiZTAYDIY+hxEvg8FgMPQ5jHgZDAaDoc9hxMtgMBgMfQ4jXgaDwWDocxjxMhh6CBH5qYg8lG4/ehIRWSIiF6TbD8PAw4iXod8hIp+LSIOI1InIThF5RERye7jNw0Vka3KaUuo3Sqle69hF5DwReau32jMY0okRL0N/5WSlVC4wE9gfuCbN/mQ8IuJJtw+p0Ff9NuwdRrwM/Rql1E7gZbSIASAifhG5VUQ2i8guEblfRLLcvEEi8k8R2SMilW58ZJJtkYg8LCLb3fxnRSQHeBEY7o726kRkuIjcICKPu3YvicgPk30TkY9F5Aw3PllEXhGRChFZIyLfaOszuSOsDSJSKyIbReTbIjIFuB840G2/yi1bICKPuZ9nk4hcKyJWUj1vi8htIlIB3OCmny8in7qf72URGZPU9jEislpEqkXkbkDa8dMSkatFZL2IlIvI0yJS5OaNFRElIue6/w5lIvKzLtp+V0Q2A6+56ee4n7FcRK5zR+BHi8hQEakXkeKk+g9wvxNvW/4bMhsjXoZ+jSs8JwDrkpJ/C0xCC9o+wAjg526eBTwMjAFGAw3A3Um2/wdkA/sBg4HblFJBt43tSqlcN2xv4cqfgTOT/JrqtvG8K36vuGUGu+XuFZH9Wvk8OcCdwAlKqTzgIGCZUupT4ELgHbf9QtfkLqAAGA/MB84BvpNU5Vxgg9vujSJyGvBT4AygFHgTeMJtuwT4K3AtUAKsBw5u6WMSPwJOc9sdDlQC97QocwiwL3AU8HNXhDtrOx+YAhznfp/3At8GhrmfeQQ0/oBZAiT/IDgLeFIpFW3Hf0Mmo5QywYR+FYDPgTqgFlDAYqDQzRMgCExIKn8gsLGNumYClW58GOAAg1opdziwtUXaDcDjbjzPbXeMe38j8Ec3/k3gzRa2vweub6WdHKAK+CqQ1SLvPOCtpHsbCANTk9K+DyxJKr+5RR0vAt9NureAerTQngO8m5QnwFbggja+u0+Bo5LuhwFRwAOMdf9tRiblvw8s6ILt+KT8nwNPJN1nAxHg6KTv+O2k72UnMCfdf6smpB7MyMvQXzlN6ZHJ4cBk9EgB9GgiG/hQRKrc6bWX3HREJFtEfu9OP9UAbwCFImIDo4AKpVRlV51RStUCzwML3KQFwJ/c+BhgbsIf16dvA0NbqSeI7ogvBHaIyPMiMrmNZksAH7ApKW0T7ojEZUsLmzHAHUl+VKBFagR6BNRYXmklaGnfsq6/J9X1KRAHhiSV2ZkUrwdyu2Cb3HZL3+qB8qT854CpIjIeOAaoVkq9347vhgzHiJehX6OUeh14BLjVTSpDTwXup5QqdEOB0ps7AH6Cnsaaq5TKBw5z0wXdORaJSCFfpDOvZ3gCOFNEDgSygH+76VuA15P8KVR66u+iNj7Ty0qpY9CjkdXAg234UIYerYxJShsNbGvH7y3A91v4kqWU+g+wAy3gAIiIJN+3whb09GZyXQGl1LZ2bLpim+z7DiB5bTILaFzjUkqFgKfRPwrORk//GvowRrwMA4HbgWNEZKZSykF39reJyGAAERkhIse5ZfPQ4lblbhC4PlGJUmoHelrtXndjh1dEEuK2CygWkYJ2/HgBLSS/BJ5yfQH4JzBJRM526/SKyOyk9Z9GRGSIiJzirn2F0dOj8SQfRoqIz/U3ju6wbxSRPHfjxWXA4+34eD9wTWK9zd3w8XU373lgPxE5Q/QOvx/RyuiwRV03JjZ8iEipiJzaTvm9sX0GOFlEDnI//y/44maSx9BTpafQ/ndg6AMY8TL0e5RSe9Ad13Vu0lXoDRzvulODr6JHW6CFLgs9ankXPaWYzNno0cxqYDdwqdvGavTIaoM71TW8FT/CwN+Ao9GbMxLptcCx6KnE7eiptN8C/lY+joUeHW5HT+nNBy52814DVgI7RaTMTbsEvda2AXjLbfePrX5R2pe/u20/6X43K9CbUVBKlQFfB25GT8lNBN5uqy7gDmAR8C8RqUV/n3PbKZ+yrVJqJfqzPokehdWi/33CSWXeRq9Z/lcp9Xkn/TBkKKKnrQ0Gg6H/IPqh9CpgolJqY1L6a8CflVL9+uSTgYAZeRkMhn6BiJzsbrjJQa9xLkfvPE3kzwa+DDyVHg8N3YkRL4PB0F84FT2duh09pbnA3RGJiDyKnh6+1J2mNfRxzLShwWAwGPocZuRlMBgMhj6HOdCyGykpKVFjx45NtxsGg8HQp/jwww/LlFKlXbEx4tWNjB07lqVLl6bbDYPBYOhTiMimjks1x0wbGgwGg6HPYcTLYDAYDH0OI14Gg8Fg6HMY8TIYDAZDn8OIl8FgMBj6HEa8DAaDwdDnMOJlMBgMhj6HES9DWgiGY1QGI+l2w2Aw9FGMeBnSwhPvb+buf69LtxsGg6GPYsTLkBb21IUJReMdFzQYDIZWMMdDGdKCmTI0GAx7gxl5GdJCZX2UaNy8jsdgMKSGES9DWqgMRojEnXS7YTAY+ihGvAxpobI+QjRmxMtgMKSGES9DWtDThka8DAZDahjxMvQ6jqOoqo+YNS+DwZAyRrwMvU5NKIqjMGteBoMhZYx4GXqdCnebvJk2NBgMqWLEy9DrVNZH8dpixMtgMKSMES9Dr1MZjDA4L0A0Zta8DAZDahjxMvQ6FfURBuf7iTpm5GUwGFLDiJeh16mqjzA4z2+mDQ0GQ8oY8TL0OhXBqJk2NBgMe4URL0OvY0ZeBoNhbzHiZeh1KoIRhuQHzHNeBoMhZYx4GXqdyvoIpflm5GUwGFLHiJeh16msjzIkL2COhzIYDCljxMvQ61QGI5Tk+XCUIu4YATMYDF3HiJehV3EcRXVDlEHZPry2ZaYODQZDShjxMvQqtaEYWT4br23hM+JlMBhSxIiXoVepqI8wKNsH4J5vaKYNDQZD1zHiZehVKusjDMpJiJcZeRkMhtTIOPESkeNFZI2IrBORq1vJFxG5083/RES+3JGtiBSJyCsista9DnLTi0Xk3yJSJyJ3t2jnABFZ7tZ1p4hIT37ugUJlMMKgbC+gxSsSM+JlMBi6TkaJl4jYwD3ACcBU4EwRmdqi2AnARDcsBO7rhO3VwGKl1ERgsXsPEAKuAy5vxZ373PoTbR3fDR9xwFNZH6XInTb0eczIy2AwpEZGiRcwB1inlNqglIoATwKntihzKvCY0rwLFIrIsA5sTwUedeOPAqcBKKWCSqm30CLWiFtfvlLqHaWUAh5L2Bj2jspg8rShWfMyGAypkWniNQLYknS/1U3rTJn2bIcopXYAuNfBnfBjawd+ACAiC0VkqYgs3bNnTwfVGvSGjaZpQzPyMhgMqZBp4tXaulLLn+ZtlemMbXf6oROVekApNUspNau0tDTF5gYOVS02bJjzDQ0GQypkmnhtBUYl3Y8EtneyTHu2u9ypwMSU4O5O+DGyAz8MKVARbNoq77MtYmba0GAwpECmidcHwEQRGSciPmABsKhFmUXAOe6uw3lAtTsV2J7tIuBcN34u8Fx7Trj11YrIPHeX4Tkd2Rg6R2V9tFG8PLaYaUODwZASnnQ7kIxSKiYiPwReBmzgj0qplSJyoZt/P/ACcCKwDqgHvtOerVv1zcDTIvJdYDPw9USbIvI5kA/4ROQ04Fil1CrgIuARIAt40Q2GvaQyGKHITBsaDIa9JKPEC0Ap9QJaoJLT7k+KK+AHnbV108uBo9qwGdtG+lJgWmf9NnSOypYbNsxzXgaDIQUybdrQ0I9RSlFVH6Ww8Tkvs1XeYDCkhhEvQ69RE4oR8Nr4PPrPzmyVNxgMqWLEy9Br6G3y3sZ7s+ZlMBhSxYiXoddI3iYPZuRlMBhSx4iXodeorG8uXj5bzIYNg8GQEka8DL1GZTDauE0eEiMvs2HDYDB0HSNehl6jsj5CYXbSmpfHrHkZDIbUMOJl6DUq6yONr0MBs+ZlMBhSx4iXodeoCEYpzGmx5mXEy2AwpIARL0OvUV4XptiseRkMhm7AiJeh19hQFmRscU7jvde2iJjdhgaDIQUy7mxDQ/8kGnfYXFHP+NIk8fKYNa9m1FfA1g9gy3tQthbiUXBiYNkwdAaMngcjZ0MgP92eGgxpx4iXoVfYXFHP0PwAAa/dmGbWvIBQDaz4K3z0f1qwRnwZRs6BaWeA7QfLA/EwbPsvvPm/sH0ZTP4KHH41FE9It/cGQ9ow4mXoFdbtrmOfwbnN0gb0mlfVZnjjVlj5LIw/DOZfDROOBLuN/5JTTtbXcC28ex88dDRMOQkO/ynkD+s9vw2GDMGseRl6hdbEyzMQt8rXbIfnfwK/PwxySuGSpfDNx2HSsW0LVzL+PJh/JVzyIWQNggePgC0f9LzfBkOGYcTL0Cus313HhKT1Lhhg04axMLzxO7jvYPBmww8/hKOug9zBqdWXXQTH/BJOuh2eWADLnuhefw2GDMdMGxp6hXV76vj2vNHN0gbMtOG6V+GFK2HwFFi4BAaN6b669z0ezvsnPHEmVKyHI6/tvroNhgzGiJehx1FKsX53HfuU5jVL7/cnbNRXwItX6d2DJ/4OJh3XM+0MngIXLIaHT4DsYph3Uc+0YzBkEGba0NDj7KwJkeXzUJB0riH08+e8Pv0n3HugFpOL3+k54UqQUwxnPQNv3wGrnuvZtgyGDMCMvAw9jt6skfOFdJ+nH655NVTBi1fC1qXw9UdgzIG913bhaDjzSXj8DMgdCqPn9l7bBkMvY0Zehh6ntZ2G0A/XvDa+CfcfoncEXvhW7wpXguEz4fQH4KmzoHZn77dvMPQSRrwMPc663XXsU9qWePWDkVcsAv+6Dv72PTjpNvjK/4IvO33+TDwaZn0HnvsBqH7048BgSMKIl6HH0SOvvC+ke+1+8D6vsrXw0FFQvk6PtiYek26PNIddoTeMfPBQuj0xGHoEI16GHmf9nmCr04a+vjzyUgo+fBT+eBwccB4s+DPklKTbqyZsL5zxACy5SQuswdDPMBs2DD1KdX2UUDTOkHz/F/K8HiEa64PTWvUV8I8fQcVGOO8FGDw53R61TslEOOKnejrzu69oQTMY+glm5GXoUdbtqWVCaQ4i8oW8PrnmtfFNuP9QKBiln63KVOFKMOu7ECiE936fbk8Mhm7FjLwMPcq63XVMaGXKEPrYmlcsAkt+o49hOvXuzFnb6ggR/YD0H46F6V+HvCHp9shg6BbMyMvQo7S1TR760JpX2Vr4wzGwa1VmbcroLCUTYf+z4NUb0u2JwdBtZJx4icjxIrJGRNaJyNWt5IuI3OnmfyIiX+7IVkSKROQVEVnrXgcl5V3jll8jIsclpS9x05a5IcUTVAc2bW2TB/DaktnPeSkFS/+oN2V8+Wz41lOQW5pur1Jj/pWw4d+w5f10e2IwdAsZJV4iYgP3ACcAU4EzRWRqi2InABPdsBC4rxO2VwOLlVITgcXuPW7+AmA/4HjgXreeBN9WSs10w+7u/rwDgbXtjLxsS3CUIu5koIDV7oQ/f0PvKPzOizD7Aj0F11fx58HRv4AXrgAnnm5vDIa9JqPEC5gDrFNKbVBKRYAngVNblDkVeExp3gUKRWRYB7anAo+68UeB05LSn1RKhZVSG4F1bj2GbmBzeT2haJyxxV88GgpARDJz08bKZ/VJGcNmwgWvQum+6faoe5jxDfBmwbI/pdsTg2GvyTTxGgFsSbrf6qZ1pkx7tkOUUjsA3GtiCrCj9h52pwyvk9a2ywEislBElorI0j179nT0+QYUr366iyP2HYxltT1i8VpCLFNGXvUV8JfvwGu/ggVPwJE/61/by0XgmF/Bkpsh2pBubwyGvSLTxKu1Xq5lz9ZWmc7YdqW9byulpgOHuuHs1ipQSj2glJqllJpVWtpH10N6iNdW7+aoKe3vbvN6LKKZcLJ84hT4/OF6U8ao2en2qGcYNVuPKM3JG4Y+TqaJ11ZgVNL9SGB7J8u0Z7vLnVrEvSbWr9q0UUptc6+1wJ8x04ldojYUZdmWKg6d2P6pE2mfNqzbDX85D175uT4F/rgb9dRaf+bIa+Gt2yFUk25PDIaUyTTx+gCYKCLjRMSH3kyxqEWZRcA57q7DeUC1OxXYnu0i4Fw3fi7wXFL6AhHxi8g49CaQ90XEIyIlACLiBU4CVvTEB+6vvPFZGQeMGUSOv/1HCX3petZLKf3M1n0HQeEYuOjt9JwCnw6GTNXb/d+5O92eGAwpk1EPKSulYiLyQ+BlwAb+qJRaKSIXuvn3Ay8AJ6I3V9QD32nP1q36ZuBpEfkusBn4umuzUkSeBlYBMeAHSqm4iOQAL7vCZQOvAg/2/DfQf1i8ehdHTen46YK0bJcvWwvPX6bfvfXtZ/RrRAYah18DD8yHOQsz60xGg6GTZJR4ASilXkALVHLa/UlxBfygs7ZuejlwVBs2NwI3tkgLAgd01XeDJu4olqzZw2XHTOqwbK9OG0ZD8Nb/g/cf1Keuz1kIdsb9F+gdBo3RJ268+f/g+N+k2xuDoctk2rShoR/w0eZKBuf5GTmo43daeW2LSE9v2FAKVr8A98yB3e4pGQdePHCFK8GhP9Hb5mt3pdsTg6HLDPD/vYaeYPHq3RzdwS7DBF5PD4+8ytbBS1dB1WY4+XaYcGTPtdXXyBsKM74J/7lTb1QxGPoQZuRl6HYWf7qLIzux3gXg66k1r/oKePEqfSbhuPlw4dtGuFrj4B/DR49DnXlG0dC3MOJl6FbW7a6jIhhh5sjCTpXv9jWvWBjeuQfung1ODH74ARz8I/D4uq+N/kTBCJj+NbPz0NDnMOJl6FbuW7Kes+eNbfdUjWS67bUoTlxvfb9rFmx8A857Hr7yv2YnXWc4+FL476MQLE+3JwZDpzFrXoZuY3N5PYtX7+L1K47otI3X3ssTNpSC1f+Ef/9GHz57xgMD53mt7qJwFEw9Fd69F466rsebizpRdgZ3srV2K9vrtrO7YTdl9WWUNZRRE6khGA1SG6klEo+gUDjKwRabHF8Oud5c8nx5DMsZxojcEQzPHc6+g/ZlXME4bMvuuHFDvyEl8RKR+cAVQBVwDlBiTl033Pf6Os6aO4aCrM6fB+jzpLjmpRR89pIWLYCjrodJx/Xtk9/TySGX6ee+DvohZA3quHwnCMfDrK9az2eVn7GhagMbqzeysWYj2+u2U5pVyoi8EQzPGc7g7MHsM2gf5g2fR4GvgFxfLrneXHy2D0ssBCGu4gSjQeqiddSEa9gR3MG2um0s3ryYe5bdQ0WogslFk5k9dDaHjzycKcVTsMRMLPVnUh15/Q793NRzSilHRB5Dv1LEMEDZVtXAC8t3suTyw7tk1+U1LycOny6CN/9Xn0J5xDWw74lGtPaWQWNg0gn6Gbj5V3bZvCZSw6ryVawuX82nFZ+yumI12+q2MSpvFJMGTWKfwn04bZ/TGFcwjlF5o/B284HH1eFqVpSt4J3t73D1m1cTjAY5avRRfG3S19i3qJ+8FcDQjFTFq1YpVZt00LpZDR/g3L9kPQvmjGJQTtf+FDq95hWLwPK/wNu3gz8fjrjWjLS6m0P+Bx4+AeZdDP7W38EG0BBr4NPyT1lRtoIV5StYWbaSPQ17mFI0hSnFUzhw+IGcP+18xheM73aRaosCfwEHjziYg0cczOWzL+fz6s95YeMLXLz4YobmDOWb+36TE8adgNfqR28JGOC0K14ikq+Uau30zodF5E/od2kdBJi/iAHMzuoQiz7ezuKfzO+ybYcjr4Yq+PBheO/3UDoZTrgFxh9uRKsnKJ0EYw+GDx/R04eAoxw2Vm/kkz2fsLxsOcvLlvN59edMKJzAtJJpHDT8IBZOX5hxa05jC8Zy8cyLWThjIW9ufZPHP32c+z++n4u+dBEnjjsxo3w1pEZHI6+tIvIIcJdSam0iUSn1uIh8BJzuhvN7zkVDJqOU4ufPreBbc0dTkuvvsr3XltY3bJSt1YK1/C8w8Vj41tMwbEY3eGxoj4q532P5cwv5xBflk4pVrCxbSYG/gOml05lRMoNT9zmVyUWT8dtd/7dOBx7LwxGjj+CI0Ufw/o73uXvZ3Ty0/CGumnMVBw0/KN3uGfaCjsTrQOCHwIci8hZwu1LqX6APtQVWtmds6P88/u4mtlU1cNe39k/JXo+83A0bThzW/gs++APsWAZfPhcufke/Y8vQ7UTiEVZXrGZ52XI+2fMJn+z5hKpwFdMK85m+awVnzTiH6aXTKQoUpdvVbmHOsDk8OvRRXt/6Or/4zy/Yf8j+XDHrCoqzitPtmiEFRJ9z20EhkUL06OpiIArcBTyilKrvWff6FrNmzVJLly5Ntxu9xqc7avj2Q+/xzIUHMr607TWS9rj5xdUMtSo5L+st+PBRyB0Ms74L074K3kA3ezxwcZTD5prNLC9bzoqyFSwvW87ayrWMzh/NjNIZzCiZwfSS6YwvHI+15QP42wVwyX/715ukk6iP1nPfx/exaP0irpp9FSeOPzHdLg1oRORDpdSsLtl0RrySGhDgK8AlwCzgYaXU5V3ysh8zkMSrPhLjlLvf5qL5E/jqASO7XkE8Cmv/xfqX72NEzTICXzoDZp0/MF9P0s0opdhVv4uVZStZUb6CFWUrWFm+kjxvHvuV7Mf0kulML5nO1OKpZHvbODz5kZNg5rdh5pm963wvs6p8FVe+cSWzhsziqjlXkeXp5y8izVC6XbxE5CygwA35SfECYCYwVCllVj5dBop4xR3FpU8twxa4fUEXpguVgp2fwMdP6rWsogm8knUcqwqP5McnGtFKhYRQrSpf1Sw4ymG/kv2YVjKNacXT2K9kP0qyunDayIYl8Pzl8IP3oJ9vbghGg/zq3V+xpmINt86/lQmFE9Lt0oAjFfHqaM3rMfSDyA8DW4DPgJoWwTCAiMUdLnv6YyqDER48p5N/axUbYcUzsPyvEA3CjAVw/stQPIH1r6+nPhjpWaf7CY5y2FK7RT9HVb6a1RX6mSqAKcVTmFo0la9O/CrXzbuOoTlDkb3ZkTluPmQVwqrnYNoZ3fQJMpMcbw43HXITz657lvNfPp8bD7mRQ0Yckm63DB3QkXjtD/wYOBt4kha7Dg0Di2jc4dInl1EXjvHQubMIeNv5RV61WXd8K/8OlZtgv9PgpNtg1Fywmk4+6LazDfsZwWiQtZVr+azyMz6r/IzVFatZW7mWQn8hk4smM7l4MgsmL2By0WSGZA/ZO6FqDRH9ws7Fv4T9Tu/3jyaICKdPPJ2xBWO5bMllLJyxkDMn9+8p01RRSqGUgxN3UPE4jhPHcdx43I27aU48jkq+OnFt5yTSdVoqtCteSqmPgfNFpAT4PvCaiCwD7lBKvZpSi4Y+SW0oyv889TGOUvz+7ANaF66ytfqcwVWLoGoTTP4KHPFT/Su+jYV//UqUgSte0XiUjTUbWVe5jnVV61hbuZa1VWupCFUwoWACk4omMWnQJI4fezyTiiaR78vvPecmHguv/Uofw7XvCb3XbhrZf/D+PHbCY/xg8Q/YXLOZK2Zf0aljppRSuuOOx3SHHI/hxOPEY7Gk9HhSiOHE3I4/FiOeEIFYK2UT6U68WRvK0fUrVxCa7JLiSWLRrL5k4UkWnbbSmolPHESwbRuxbSzLwrKa4mJZWLat0xrjFmLZWLbV3MYtnwpd3bDhAb6GHo3lobfOP5RSy/2Q/rrm9cnWKi554iMO3qeE60+eit/jCpcTh60fwJoXdQjXaMGafBKMPaRTO9WefH8zH22u4rdf69/PcIViITbVbGJD9QbWV61vvG6r28bw3OHsU7gP+xTuw8RBE5lYOJFReaMy40Halc/ql1VesLhXRl9KKZTjEI9GicdixGP66sRijfdN8RhOLEo8nnzf4hpvWVaLRSK9UThiMeLxWGNeJBphffk6/OJlSGAwKh7X+YlOvbH+JrERy8K2PVge3XFbHo/uuG0Plm1h2R5s2023bF3OTpRpXi6RZtse3dknp1st762mOu2mILadJDBuObcNsZrstd9JQpTUfnPxcfP2QnDaotvXvETkNJo2aySHjcChwO8BI179FMdR/PHtjdy3ZD2/PHUaX5kxDIJlsHIxrHtVh4IRMOl4OO0+GL5/synBztDt7/NKI0op9jTsYVPNJjZWb+Tzms/1tfpz9jTsYWTuSMYXjmd8wXiOGXMM35/xfcYVjMNnp+90NeU4xGJR4pEosWiEeLTpGo9GianxxHdGiP3zQeIlU7WYJPIS5RJpCbFJ5CcEwy0Ti8aIR6M4sZZ5TUISi0UREWyPF9vjwfZ6sTy600+kWR5Pi7i+t2zbTfNiu8Jge91yto0vK7uxLqtZ+aY02xUR2+PhIIlx20d3UO6Dy2Zfid/n1523x25qL0msurtDN7RPR2teN6A3bFQnhXJgA7DEvTf0Q97bUM4v/7mKfG+cF08VBu95EB74N5Svg3GHwT5HwdHXQ0EK2+ST8Hr61pqXUoryUDmbazazpXYLm2o2sbl2M5trNrOpZhMBT4Cx+WMZkz+GcQXjmD1kNmMLxjIyb2S75+opxyEWiRCNhIlFwsQiETe0EY9GmsXjraVFo8QiEeLRCLFoVF8bRUrHHSeOx+PVnbzXi8fnw/b68Hg82D6fFonwFDwvPYs9vszN9+p018b2eBuFwZOUlhAO2+vFtr3YXk9TXsv8JCGyMmHE6XLrpN9z5RtXcsPaW7j98NsJeMyzh5lCl6YNDe3TH6YN1++q5i//fJ68He/w9aINlFZ9jJRMhAlHwPgj9IaLbnwr8UsrdvD3j7bx+7O7NGPQo4TjYbbVbWNb7Ta21m1la80WtlVuYUfVVnZX7yRb+RkeGMJQbymDvSWUeAoZZBeQb+XiiQvRcIhoOItug0sAACAASURBVEwsrIUoGk4IUrgxPRpJukYixGMxPF4fHr8fj8+H1+fT9z59r0NT3PYmpXm9bp6+2l6va6vraCZKyXGfF8v2dLzZw4nDvfPg+Jv1j5YBRsyJcc2b1xCMBrn9iNvTOlLur/TEVnlDfycWhu3L2PrxYipX/ZuxDSv4XvZQCmYejWfChTDmYMjuueOBmh0P1Q04TpxIQwORhgaioQYiIR2PhBqIhkJEGhoIBmuoqi2jqq6c2mAV9fW1NNTXEQ2FiEXCSNTB73jwxm3sGFgxh9FeDxP8fvyByfgCWXgDAbw+P96A4PGHwVdNQyCEx+fH6w+QU1CI1x9oEiN/AK/f3yg4Ou5vzPd4fd2/Y7C7sGyYfxUsuQkmHNnvdx62xGN5+M2hv+GK16/gitev4NbDbzWn02cARrwGGsEy2PI+bH2f2OfvoHZ8zGaG85FMoWTaN9ln/uPkDRrSa+4k1rzisSjh+nrC9UEi9fU63qDjkQZ9H2lIhAbCSfGmawPxaBRPwIf4vOCziXsg4nEIS4x6K0wt9YSsGFlZOWRn55GfN4iC4cMYmldKaeFQhhaMoDR/KIHsHLx+P75AFh6fz6xn7Hc6vH4LrF8M+xydbm96Ha/l5XeH/Y5Ll1zKT9/8KTcfenNmbKgZwJhpw24k46YNoyHYuRy2LYVtH8LWpahgGbsLpvNBbAKLKsdQMPEgTpo9iUP2KcG2Uv9F7ThxwsEgobpaQsE6wnV1hOqDhIN1hIJBwm48HAwSbqhvjNfW1hEK1uEVhS8rG392Nr7sHH3Nysafpe99WVnEvBCyo4SsKLXSQDVBqpxayuKV7I5XsDO6hxonSGl2KUNzhjIke4i+5ujr0JyhDM0eSlGgKHNHOZnMir/CO/fCBa/229GXUkq/5NRRoBTKUeCgr0oRjoa59s1rGZY9lEv3v9Qtp8vrsvo+UR5HoRrrI6lMU92opDKO0j646c3tFcpp4ZtqYZ/wRbX0JclOJdWdaCdRT7M6WvjQWp2ubXJ7zco31pu4/2KbKBj128N69mxDQ/ukVbyiDbBrlT6Nfccy2P4RlK1DlexDReEMPlETeLFyBC/tKmDuhFKOmTqE4/YbSkFW8+kP5TiEgnU01NYSqqtxr7U01NZoYaqrpaGurjEeqqsjHKwjEmrAn5WNPzeXQE4u/pxcAtk5TffZOfhzcvHn5Oj0nBx8Wdl8tKeSu/7zKdd+dQIVoQrKGsooD5VT1lCm4w3l7GnYQ1Woinx/PqVZpQzOHszg7MGUZJUwOHswQ7KHUJpdypDsIQwKDDKvf+8iiQ618RpXTR1pPCkvFkM9dS7M/SFqxOzmNo6CuPpimpPUSbdSvlknn5yWLBpOUrl4U0fdVK6p0/1CWqMg4Nom59NCCFw/ASwBC8QSLdSWIBYgghKoiFTg8/goCBTosiK6rKVtxbVBEnXQvFwi3Q0iNLbT0qaxLotW6nXbFLdcow1JdSe139gWzXxp/IyNfibVmfyZmn0O105oO48W/iSXT7K3LKvnxUtE3lBKHdYlowFCr4lX3W7YtVKPqnat0NeKjVCyD7HB09kamMTS2FgWl5fy7uZahvsd5gz1Mq3IZkIexOvraKitpr6mmoaaGhpqE9caQsE6fIEssvLyycrLJ5CXR1ZuHoG8fAK5uWTl6msgJ9dNyyOQm4s/KxtEqIvWURmqpDJcqa+hSipCFY0hcV8eKqcyVInX8hMJZzNj2EiKs4opySqhOFBMaXapjmcVU5pVSlGgCI+VGbPciQ5QxRXEnKZON65Qcacxrss4Op7olFsrn6grrp9xasxrLNeyTAthSbTZUhiSxSTuNBeHZJFQNHVsdqKjFrClsYNtTI/UQHAXMmRfsC3dUSXbJMoKiG0172gT9SXZNLOzAMtyr807bOV2gMrtUFUivbHTpLEjVNLUjrJ0tkp0+tBYBrcTV9CUZyV69sTowf33hqZ7FNXhan725s84buxxnDT+5GblE38jyeWb0hv/iBJJbl5y+aT85DpIzm/yq3l6kq+NdTf50MxEqQ7LtPRTJd18oY0kw5Z1tExLrgfgS0eO7hXxipvDeFun28UrVAN71sDuVbBntRas3auIx6IEC6dS4R/P5lgpG4M5bKy2qKioJh6spsiOkq9CeKNBnEiI7Lx8svILyM7PJyu/kOyCArLzCsjKLyArX4tUdn6BFqvcPCzbJubEqInUUB2ubgxV4SqqwlXN4lXhKipDlY1xv+2n0F9IUaCIQYFBFPoLKQ4UN94XBYooyiqiOFDMoMAg1u8Kc9nTy3jp0qbfQ42desxBxRQq5mgRiLmdcCw5rsvGo3FU1MGJJtlEE3ZNNs3qSXTwSQJDi46exLXxFzuIo1BJvyRV4he5Jc06RmXpX+mJX+uN6eh447UxrUWeG5zGPEGhUAiO23k4tEhXyk1LmsUBHLczTaQ32rudoEIap5SadbyJjsdpuqrqrahAIXhzGss3dmxJ7STX1yyfpLpVYvao9TKNuL/WxY1/8V4aBxB6xJBkR+LXflM5N6mpbKOuJZWlaUSRmGZuqleIOVG21m2lJKuEfH9eY6XSRvmmOps+VLKfje01fo7mZZrVRVN+cnlpVrD57G7jZ0n6Thvrl+SkFvaJ7zD5e0sUbuZLcp60WmeybXLZ+Wfu2/fFS0SOB+4AbOAhpdTNLfLFzT8RqAfOU0r9tz1bESkCngLGAp8D31BKVbp51wDfBeLAj5RSL7vpBwCPAFnAC8CPVQdfVkripRQE90DZZzi711C/ZRV129dRt2sbtcEItZ7BVKlCKiN+akIQaohANEzIChDy5mDn5JNdUEhRSQkjhg9m7Mgh5A0aRHZBIb7cXKJ+RW20jtpILbWRWmoiNU0hXENNfR11oTrqQvUEQ/U0hBqINkSRCBRZRRRY+eRbBeSTS57kkkM2WZJFtgrgJ4BfefEpL17HgwcbK07TaCGuELfzF/dXvijd+Yty05TCiSt8liBK/4AGt8NNDqqpE07cx5Vq0hahsdNOCIKTEIFGEWn6BZ4Qk4Tw6DRpnCLBtlwxEsRuSkNwRxB6hGC5tpY7uhD3V7803tN6WrNrUwdpJUYqSXZA8/REO4lOMlE2qV5a5FvuZ9NZ7ijH7UEa60anN3agbpnGDizRAW95H1nyG+Tsv4LtbWwXkspIi3hLAUnqrJuVsZqLSqavTa6tXMsF/7qAWw67hbnD5qbbnT5Lj7/Py22kx8RLRGz0yfXHAFuBD4AzlVKrksqciH6f2InAXPQ5i3PbsxWRW4AKpdTNInI1MEgpdZWITAWeAOYAw4FXgUlKqbiIvI8+ButdtHjdqZR6sT3/Oyte9VvXsOr/bqCmvJy6mnpqo15qY1nUR/UxL3FfNrV2HuWSQ8QTIJCTS252LrlZWeT4c8jx+/A5gh1XqIiDRBUSFay4Dh5HsB0Lr7Lx4cGLjRcbj7KwsfAg2ICN7nRt0Wqv+3FpEoyECCQ6e3cqR7mBxNSPbbmVCHgtPV1kW4hHEI8FHkE8NuIVLI+FeC3Ea2N5LCrCUe58Yz2//toMLJ+N7bOxPOIeU6M7U8u2sNzpJcsSnW4Lljv91Bc6uX7HY6fBlJNg9gXp9iTtfLDzAy5//XIePPZBJg2alG53+iT94TmvOcA6pdQGABF5EjgVWJVU5lTgMXcU9K6IFIrIMPSoqi3bU4HDXftH0aeDXOWmP6mUCgMbRWQdMEdEPgfylVLvuHU9BpwGtCteneX95dsJ7TwK2wqQn+OnSPx4xIdHbGwRbFFYgisyWjhiYQcnrIhbirgoHAscCz2n7wECFpZXH3fj8Xnw+L14/V4sr4XtsxGfYHttLL++t/0Wts/SYuOxtAi4107RTRt9YjWwnSoKSmJArEUb6PFwvPkUWseudY9vGcPefp6esJ/2A/jHZVB8GPhaeaFlJ9vs9L9Vp4p1tq52yrWV147NvqqIK4eez8+f/gE3HnIjpYHiDptpUXlnC3aiqhTrcn/8dfjv0eaPxLb7jWYmzecxO+dbG2SaeI1AvzcswVb06KqjMiM6sB2ilNoBoJTaISKDk+p6t5W6om68ZfoXEJGFwEKA0aNHt/PRmhjst/l7/laarb42iydo5w8pMTwC6OOvwzpUKe6468Nea096TNx6cPSXiQNLdTI89Mhed0L9henM5ok/PaunaA09Tiri1ZP/Mq3V3bKnaatMZ2w7216n61JKPQA8AHrasIP2AJh65GFMPdJs2ASoro9y6C2v8ckNx6XbFUNXqdgADx4FF78DeUPT7U3aUUpx43s38nnNZu45+h5zCkcX+NlvftNlm1QeiHk9BZvOshUYlXQ/EtjeyTLt2e5ypxZxr7s7UdfIVtIN3YzXI916PJShFykaD18+B169Id2eZAQiwtVzrsZn+/jlO7/sf9PXGUaXxUspdURPOOLyATBRRMaJiA9YACxqUWYRcI5o5gHV7pRge7aLgHPd+LnAc0npC0TELyLjgInA+259tSIyz93deE6SjaEb6U+vRBmQHHY5bFgCWz5ItycZgcfycMtht/BZ5Wfc/8n96XanX5NRRxEopWLAD4GXgU+Bp5VSK0XkQhG50C32AvqVLOuAB4GL27N1bW4GjhGRtejdiDe7NiuBp9GbOl4CfqCUSryT+iL0u8rWAevpps0ahuZ4LCHmKBzH/Ertk/jz4Ogb4MUrwDE/QgCyvdncc9Q9PLfuOZ5d92y63em3mOOhupGMO9uwjzDpZy+y/BfHNr2h2dC3cBz443F6CvHLZ6fbm4xhQ9UGvvPyd7jpkJs4aMRB6XYno0llq3xGjbwMAxOvbda9+jSWBSf8Fl77FTRUptubjGF84XhuO/w2rnnrGlZXrE63OxmHUoq4o1JeNkhpq7yInAAcC2wGPgE+UUrtSckDw4DH67GIxhzwp9sTQ8qM+DJMOQVevhZOuyfd3nQapRQxRxGLK2KOQ9zR941XN91JKpeIx1vcO0m2TWlDOKr0Qs5/8ULOGXsruXYpjttpJ8rpoytbpin3/OCmdKVU46kyjtO8TKKcSrJJHP2VsFGqqf7EfaN9s7oSeU1tOUllVUvbFuUTx4u1tHHcY78S95A4YS21DeypPud1G3AeUALMAy4AzkyxLsMAx2za6CccfT3ceyDxta8RGTOfcCxOOOYQjjpE4nFCUYdI3CESSwrx5tdoXAedpojFE2mKSNwhGmu6T5SNub/eY3FF1NE2Oq6vMbdMolxCmGJxB0fpdVfbEry2hSX679F20zy2YLvHdnktfdKLbYHHsvC4J77YostZIk12yXnWFEbICTyy/qccnH09WZ78xjpt1yYRtyzBY1n4PbpTt5PSbfcoL0sSbZEUF/eUM7c+VxS0rU5PvheSyzXVYyWVt6R5nQKNr00SwbXX6VZSXKT5fcI3ocmvhE0Cuanrf26pitcS4GOlVAPwzxTrMBgA8NkWESNevYLjKBqicYLhGHXhGPURHa+PxqkPxwlGYjRE4tRH4jRE4zREYu7VIRSNE4rqdB13CMXihKM6LxxzmBP/Fr94fCGnxH+H483B57Hweyz3auurbeH3Wvhsne7zWHhtHfweC68tjfc5fg8+28LjpiXHdRA8toXX0lePrUXG69EC4rF0WuLqtSxsO5GX6PR746Himdz53yze2X43fzjuD2R7WzmVxNAlUhWvt4APRORp4EPgI6WUeQ7KkBJmzavzxB1FTUOUqoYolfURquujVDfoUNMQpTYco6YhSk0oSm0oRm1Ii1RdKCFWMQJemxy/h1y/hyyvTa7fQ7bfJttnk+X16KvPJstrU5LrJ8tnE/Dq+4DXJuC1msX9Hh33ey0CnmPw/mMTHwXe1etghkYu2f8SKkIVXPrvS7nnqHvw2uYh5r0hVfG6BjgLPW34JeCb6GehDIYuM9CnDSMxh921IXbXhtldE2ZPXZg9tWHK68KU10WoCEYoD4apdIUqx2dTmO1jULaX/Cwvhdk+CrI8FGR5Kc7xMa4kh7yAh7yAV1/9HnIDWqyyfZ69emN2pzjuJrj3QJh6Kowxu+wSiAjXzruWnyz5CVe/eTW/Pey3GfOOur5Iqt/cv4A1SqmPgFe60R/DAMRrW0Ri/Ve8qhuibKmo16Gynu1VIbZVNbC9qoGd1SFqQlFKcv0MzvNTmhegNM9PaZ6fyUPzKM71U5TjozjHx6AcH4VZXjx2hm8Szi6CU+6Evy2E77+h7w2A+xDz/Fu4ZPEl/Pztn/PrQ35t3vydIu2Kl4jkK6VqWsk6DPhYRJ4APgKWKaU29oSDhv6Pnjbs2+IVdxSbK+pZu6uWtbvrWL+njo1lQT4vCxKJOYwqytZhkL7OHVfEsMIshhcEKM719/xoqLeZdBxseB0WXQLffNwc3puE3/Zzx5F3cPGrF/PLd37J9Qdeb17pkwIdjby2isgjwF1KqbWJRKXUASJSDMxww4nA93rMS0O/Rk8b9p01r1jcYfXOWj7eWsWKbTWs2lHDZztrKcrxMWlILhOH5DF3XBFnzhnN2OIcSnJ9A7NzOvp6+MOx8MFDMMd0D8lkebK4+6i7+f4r3+em92/imjnXDMy/kb2gI/E6EH3k0oci8hZwu1LqXwBKqXLg324wGFIm09e8QtE4/91UyTsbynlvYwUrt1UzrDCLmaMKmT6igDO+PILJQ/PIC5gF+GZ4/PC1P8IfjoFRc2HYjHR7lFHkeHO49+h7ueiVi/jVu7/i2nnXminELtCp46FEpBA4H32OYBS4C3hEKVXfs+71LczxUKlxzh/f5/yDx3L4voM7LtxLbC6v59VPd/Ha6t38d3Mlk4fmMW98MXPHFzNzVCEFWUaoOs3yZ2DxL+CC1yC3NN3eZBzBaJCLX72YkXkj+eVBv8S2Bt4xaakcD9Wlsw3dE9a/AlwCzAIeVkpd3iUv+zFGvFLjgkc/4JuzR3PM1CFp9WNbVQPPLdvGomXbKQ9GOHLfwRw5ZTAH71NCrt/sCtsrXvs1bHwDzlkE3kC6vck46qP1/PjfP6bAX8BNh9w04LbRd7t4ichZQIEb8pPiBcBMYKhSauD9TGgDI16pcdHjH3Lyl4Zz4vRhvd52LO7wyqpdPPrO56zZWcvx04Zx+v4jmDVmULMTAAx7iePAM+eB7YczHjAbOFohHA9zxetXUB+t57YjbiPPl5dul3qNVMSro5+TjwFVwMPAFuAzoKZFMBj2inSseQXDMf703iYe/c8mhhUEOO/gsRw7dSg+j1lz6BEsC067Hx45Ed74Hcy/Mt0eZRx+289th9/GTe/fxLkvncu9R93L0Bzzhuq26Ei89gd+DJwNPEmLXYcGQ3fQm895ReMOT76/mTtfW8fccUXcf9YBTB9Z0CttD3h82XDmk/DH4/V7wOZdlG6PMg7bsvnZ3J/x8MqHOfvFs7njiDuYWjw13W5lJO2Kl1LqY+B8ESkBvg+8JiLLgDuUUq/2hoOG/o/P0zvHQ73+2R6uf24Fo4qyefi82UwbYUSr18kbCucugke+ApbHbKFvBRHh/GnnMyJ3BBe+ciGXz76cUyackm63Mo5OrUIrpcqAG0Xkt8DXgF+JyO3orfMP9aSDhv5PT08b1kdi3Pj8pyxZs4cbT5+WUbsaBySFo/XGjUdOAtsHB5ybbo8ykuPGHsf4gvFc+u9LWb5nOVfOvnLAbeRoj3Yn+EXkNBE5V0QuEZGfATcC84GN6E0bv+8FHw39nJ4Ur2VbqjjxjjcJRR1evPRQI1yZQtE4PQJ7/RZ4+07owq7ngcTEQRN54qQn2Bncydkvns3GanOQUYKORl43oDdsVCeFcmAD+rUo1T3om2GA4O2hV6K8umoXV/71E248bRonpGEno6EDiifAd1+GP30dqrfA8TfDAHzGqSPyffnceeSdPLXmKc558RwunnkxC/ZdMOBP5OhozWtmbzliGLj4bCEa695f3k8v3cItL63hD+fOYv/Rg7q1bkM3UjASzn8Jnjpbh68+CL6cdHuVcYgICyYvYO6wufz0zZ/y2ubXuHbetYzJH5Nu19KG2RdsSDvdPW34wBvruXPxWp76/jwjXH2BQAF8+xl9+vwDh8POFen2KGMZVzCOx058jENGHMJZL5zFXR/dRUOsId1upQUjXoa04/V0n3i9uHwHj/5nE89ceBATSnO7pU5DL+Dxwal3wyGXwWOnwNKHzTpYG3gtL+fudy5/OfkvbKrZxGnPnsZz654j7sTT7VqvYsTLkHa6a81rzc5afvbsCu4/6wCGFpgjiPokM8+E77ykT6J/8ttQvTXdHmUsQ3OGcuv8W/nNob/hb2v/xumLTuflz1/GUZl7yHV3YsTLkHa6431e1fVRFv7fUq47aYp56LivUzoJvveaPoX+/kPh3ftggI0qusIBQw7gkeMf4arZV/Hoykc5+e8n86dP/0QwGky3az2KES9D2vHa1l5t2HAcxY+e/Iijpwzh9P1HdqNnhrTh8cPhV8N3/wWrn4ffz4fPXjZTiW0gIhw84mD+dOKfuPGQG/nvrv9y3F+P49fv/pplu5fRlQPY+wrmqGxD2tnbDRv/XL6D8mCYP5zQpXM9DX2Bkolw7j+0gL1yPbxxKxx5LYw7zBzu2woiwszBM5k5eCY76nbwjw3/4Lq3ryOu4hw/9njmj5rPtOJp/eK1K0a8DGnHa0vKa17hWJzfvbyaW776JTy2mUjol4jAlJNg3xNgxV/h+cvAm63PRtzvDPOKlTYYljuMhTMW8r3p32Nl+Upe2fQKN/znBipCFRw0/CBmDZnF/kP2Z1z+uD75zJgRL0Pa8e3FyOv/3tnEpMF5HDihuJu9MmQclg0zvgHTvgbrF8N79+vR2PSvwfSvw/D9zWisFUSEaSXTmFYyjf854H/YVreNt7e9zdJdS3ngkwdoiDUwuWgy+xbty6RBkxhfMJ6ReSMp8Gf22nHGiJeIFAFPAWOBz4FvKKUqWyl3PHAHYAMPKaVu7sheRK4BvgvEgR8ppV520w8AHgGygBeAHyullIicB/wO2OY2e7c5w7Hn8NoWsRQO5q2uj3LfkvU8uXBeD3hlyFgsCyYeo0P5evjkKXjmfC1uU0+FicfCiFlgZ0z3llGMyB3BN/b9Bt/Y9xsA7AruYk3lGtZUrOH1ra/z2KrH2FK7BY/lYXjOcEqySijNLqU4UEyBv4A8Xx55vjyyPFn4bT8BO4DH8mCJhS02IoKjnMYQiYeJhmuJhKsJhaoIhasJhWuoD9cQjNSkvLEkk/51rwYWK6VuFpGr3furkguIiA3cAxwDbAU+EJFFSqlVbdmLyFRgAbAfMBx4VUQmKaXiwH3AQuBdtHgdD7zoNveUUuqHPfuRDaCf80pl2vCeJes4dr8hTBwycF7aZ2hB8QQ44qdw+DWw7UNY8wK8cDlUbdHrYmMOgtHzYMh0I2ZtMCRnCENyhnDYyMMa05RSVDaUsaNyPXtqt7KnbjvldbvZVfYZayPV1EbqCMVDhOJhwvEIMSdGXMVxVBylHGzHQZSDpRx8ThwvgldsssQmy/Lit7zk2AGyPQEKPFkp+Z1J/5qnAoe78UfRZyde1aLMHGCdUmoDgIg86dqtasf+VOBJpVQY2Cgi64A5IvI5kK+Ueset6zHgNJrEy9BLpLJVfntVA08v3cK/Lj2s48KG/o8IjJylw1E/h5rtsOF12PIufPgoVG2GwVNg6HQYOg1KJkHRBMgf3v+mGuNRCNfqEKlriodrkuIt0kI1zdIkXEtRtJ4ibw4E8vX71/x54E/Eh0N2gc7z5TYv48trSvfnapsONohcQNf/DTJJvIYopXYAKKV2iEhrx3+PQL/ROcFWYG4H9iPQI6tkmxFA1I23TE/wVRE5DP326P9RSiW3a+hG9JpX16YN//rhVk750nAG55vFekMr5A/XDzzPPFPfh6ph10rYuRy2L4Plz+gpx0gdFIyCghHaJm8YZJdATok+ripQAH63k/Zmgzdr7w8PVgqcGMTCbgjpEG3Q10hQx6NBiNRDtF6nNYY6V5Tq3PtaHU+IVTzSJCL+5OAKSaBAC0vByBailN9Uxp+ny1iZuwmqV8VLRF4FWnuv9c86W0UraR31em3ZtFfXP4AnlFJhEbkQPZI7stXKRRaipx4ZPXp0B64YWqOrb1JWSrHo4+3cdMb0HvTK0K8IFOgpxDEHNU8P1ehTPGq26VC7E8rX6RFbsEyPTEI1+hpt0MH26efQbC9YXv1STbH0CE5Ei5NSgNIPVzsxcKIQj2lhiUd0OU+WPhbLE9DBm6Wvvhwd92ZpAfFm67dQ+/Igu1jn+3KaRja+3Cax8bvl+9toshV6VbyUUke3lSciu0RkmDtqGgbsbqXYVmBU0v1IYLsbb8u+LZutbvwLdSmlypPSHwR+285negB4AGDWrFn970nAXiDHbxOMxDpdfs2uWuojcb5sDt017C2BfAhMhSFTO1deKXekFNaiFI/oq1KgHB0SQoZoYUsE2+sGn3n1SzeQSdOGi4BzgZvd63OtlPkAmCgi49A7ARcA3+rAfhHwZxH5f+gNGxOB95VScRGpFZF5wHvAOcBdAAkRdO1PAT7tzg9qaE6u30tdqPPitWjZdk760jAsq///ujRkGCJNoyJDWsmkCc2bgWNEZC16N2FiC/xwEXkBQCkVA34IvIwWlKeVUivbs3fzn0Zv6ngJ+IG70xDgIuAhYB2wnqbNGj8SkZUi8jHwI+C8nvrQBsgNeKgLd068lFL845PtnDxjeA97ZTAYMpmMGXm5U3VHtZK+HTgx6f4F9Lb2Ttm7eTcCN7aSvhSY1kr6NcA1XXDfsBfk+GxC0ThxR2F3MJr6aEsVXttiv+H5veSdwWDIRDJp5GUYoIgIOX5Pp6YO//Hxdk750vA+eZyNwWDoPox4GTKCPL+H2nC03TJxR/H8Jzs4+UtmytBgGOgY8TJkBJ1Z93pvQzmleX7zhmSDwWDEy5AZ5AU63nG45LM9HLdfa48JGgyGgYYRL0NGkOv39a9unAAAC4FJREFUUNuBeK3YVs2XRhX2kkcGgyGTMeJlyAhyAx5q25k2VEqxcnuN2WVoMBgAI16GDCGvg92GWysbCHj/f3t3H2NHVYdx/Pt0d7tLuxdbWop9QSlYkVIVtBajRlAKFkwsEtASE0uiASNEiH8gilE0aYIgiG8YMZDUIC8NKBRTEASq8oflTV5aSqGWRrat3ZZCu1volrY//7hz6+3l3tnu7t29M5fnk2z2ztw5M+fkNPt0Zs7MGcXEzvYRrJWZZZXDyzKhs72V3pTRhqs2bmfWlGxPjmdmI8fhZZlQ6GhLveflS4ZmVs7hZZnQ2ZE+YGPlhu0cP9VnXmZW5PCyTCi0pz/n5TMvMyvn8LJM6OyoPWCje8cudu/dx9RxfpO3mRU5vCwTCh21Xw9VOuvy+wzNrMThZZnQmTJU3iMNzaySw8syoZDykPLKDTuY6ftdZlbG4WWZkDab8qpN25nlkYZmVsbhZZlQqPFW+e1vvMW23t1MnzC2AbUys6xyeFkmjElmU96zd98B61dt2s5xkw9lVD8zLJvZO4vDyzKhNJvyzr69B6xftcHPd5nZ2zm8LDOqzaa8aqPfrGFmb+fwssyoNpvy+lff4OiJvt9lZgdyeFlmVHs575aePiYVOhpUIzPLKoeXZUblg8oRwZbePg4veA4vMzuQw8syo3I25R1v7qG9ZRSHjG5pYK3MLIscXpYZlbMpd/fs8lmXmVXl8LLMKHS00rPr/6MNt/T4kqGZVefwsszobG87YLSh73eZWS0OL8uMytmUu3d4pKGZVefwssyonE3ZZ15mVktmwkvSYZIelPRS8nt8je3mSVojaa2kyw+mvKTvJtuvkfS5svWLJL0iqbfiGO2S7kjKrJB0VP1bbJUq73l179jFJIeXmVWRmfACLgceiogZwEPJ8gEktQC/Bs4AZgLnSZqZVj75fgFwPDAPuCHZD8C9wJwqdfka8FpEvA/4GfCTurTQUlW+YcNnXmZWS5bCaz6wOPm8GDiryjZzgLURsS4idgO3J+XSys8Hbo+Ivoh4GVib7IeI+GdEbOqnLncCp8pz0A+7yoeUPdrQzGrJUngdUQqS5PekKttMBV4pW+5K1qWVTytTy/4yEbEH2A5MqLahpAskPSHpiS1btvSzW0tTOZtyd0+fLxuaWVWtI3kwSX8F3l3lqysOdhdV1kUjy0TEjcCNALNnz+5vv5aifDblvj172dm3h/FjRje4VmaWRSMaXhExt9Z3kjZLmhwRmyRNBrqrbNYFHFm2PA3YmHyuVT6tTC2lMl2SWoF3Adv6KWNDVCgbKr+1dzcTxrZ7EkozqypLlw2XAguTzwuBe6ps8zgwQ9J0SaMpDsRY2k/5pcCCZAThdGAG8NgA6nIO8HBE+KxqmI0Z3ULfnuJsyr7fZWZpshReVwGnSXoJOC1ZRtIUSctg//2ni4G/AKuBJRGxKq188v0S4HngfuCiiNib7PtqSV3AGEldkq5M9nUTMEHSWuDbVBn5aPVXPpuyh8mbWZoRvWyYJiJeBU6tsn4jcGbZ8jJg2cGWT75bBCyqsv4y4LIq63cB5w6g+lYnpdmUPUzezNJk6czLbP+ElL5saGZpHF6WKaUHlT1M3szSOLwsU0oPKvvMy8zSOLwsU0qzKXf39HG43yhvZjU4vCxTDk1ezrvVlw3NLIXDyzKls73VAzbMrF8OL8uUzvY2Nrz2Ju1to+hoa+m/gJm9Izm8LFM6O1p5eetOn3WZWSqHl2VKob2VdVt6fb/LzFI5vCxTCh2tbNy+yyMNzSyVw8sypbOj+MYyn3mZWRqHl2VKZ3sxvHzPy8zSOLwsUwo+8zKzg+DwskwpdLQBPvMys3QOL8sUXzY0s4Ph8LJMGTO6hVGCSR5taGYpHF6WKZK45esnMX5MW6OrYmYZlpmZlM1KPnHMxEZXwcwyzmdeZmaWOw4vMzPLHYeXmZnljsPLzMxyx+FlZma54/AyM7PccXiZmVnuOLzMzCx3FBGNrkPTkNQDrGl0PYbRRGBroysxjNy+/GrmtkHzt+/YiCgMpIDfsFFfayJidqMrMVwkPeH25Vczt6+Z2wbvjPYNtIwvG5qZWe44vMzMLHccXvV1Y6MrMMzcvnxr5vY1c9vA7XsbD9gwM7Pc8ZmXmZnljsPLzMxyx+FVJ5LmSVojaa2kyxtdn3qTtF7Sc5KeHsyw1qyRdLOkbkkry9YdJulBSS8lv8c3so6DVaNtV0rakPTf05LObGQdh0LSkZIekbRa0ipJlyTrm6X/arUv930oqUPSY5KeSdr2o2T9gPvO97zqQFIL8CJwGtAFPA6cFxHPN7RidSRpPTA7IpriQUlJnwZ6gd9HxKxk3dXAtoi4KvkPyPiI+E4j6zkYNdp2JdAbET9tZN3qQdJkYHJEPCWpADwJnAWcT3P0X632fYmc96EkAWMjoldSG/AocAlwNgPsO5951cccYG1ErIuI3cDtwPwG18lSRMTfgW0Vq+cDi5PPiyn+wcidGm1rGhGxKSKeSj73AKuBqTRP/9VqX+5FUW+y2Jb8BIPoO4dXfUwFXilb7qJJ/rGVCeABSU9KuqDRlRkmR0TEJij+AQEmNbg+9XaxpGeTy4q5vKRWSdJRwInACpqw/yraB03Qh5JaJD0NdAMPRsSg+s7hVR+qsq7Zrsd+MiI+ApwBXJRcmrL8+A1wDHACsAm4trHVGTpJncBdwKURsaPR9am3Ku1rij6MiL0RcQIwDZgjadZg9uPwqo8u4Miy5WnAxgbVZVhExMbkdzfwJ4qXSpvN5uR+Q+m+Q3eD61M3EbE5+aOxD/gdOe+/5H7JXcAfIuKPyeqm6b9q7Wu2PoyI14HlwDwG0XcOr/p4HJghabqk0cACYGmD61Q3ksYmN46RNBY4HViZXiqXlgILk88LgXsaWJe6Kv1hSHyRHPdfctP/JmB1RFxX9lVT9F+t9jVDH0o6XNK45PMhwFzgBQbRdx5tWCfJsNXrgRbg5ohY1OAq1Y2koymebUFxJoJb894+SbcBp1CcamIz8EPgbmAJ8B7gP8C5EZG7gQ812nYKxctNAawHLizdY8gbSZ8C/gE8B+xLVn+P4n2hZui/Wu07j5z3oaQPURyQ0ULx5GlJRPxY0gQG2HcOLzMzyx1fNjQzs9xxeJmZWe44vMzMLHccXmZmljsOLzMzyx2Hl1lGSTpW0r8k9Uj6Vo1tLpR0/QjW6TpJ3xip45nV4vAyy67LgOURUYiIX1R+mTwQ/33gmvJ1kn6QTM+zM5lC4z5Jp5dts17S3Ip9nS/p0YOo0zXAFcmxzRrG4WWWXe8FVqV8Px94ISI2lK27M1n/VWA8MB34OfD5elQoeSj2BeAL9dif2WA5vMwySNLDwGeAX0nqlfT+KpudAfytrMxcinPKzY+IFRGxO/m5PyIuGcCxv5wcs/TTJ2l52SbLqVMYmg2Ww8ssgyLisxRfEXRxRHRGxItVNvsgsKZseS6wIiK6hnjsO5JjdgJTgHXAbWWbrAY+PJRjmA1Va6MrYGaDNg7oKVueCPy3tCDpMIrBI6A9IjrKtr1b0p6y5dHAU+U7lzQKuJXifbffln3VkxzbrGF85mWWX68BhbLlV4H9bx6PiG0RMQ74KNBeUfasiBhX+gG+WWX/i5L9V450LACvD7XyZkPh8DLLr2eB8nthDwEfkzRtqDuWtIDiW8zPiYi3Kr4+DnhmqMcwGwqHl1l+LQNOLi1ExAPAIxQvCZ6UDJtvAz4+kJ1KOhH4JcWzsy1VNjkZuG/w1TYbOoeXWX7dC3xA0pSydWcDfwZuoXhp72XgKxRnqz1Y8ykOs3+0bMThfbB/QsSZFOc+M2sYz+dllmOSLgBmRsSlI3S8a4F/R8QNI3E8s1ocXmZmlju+bGhmZrnj8DIzs9xxeJmZWe44vMzMLHccXmZmljsOLzMzyx2Hl5mZ5c7/AFj+djE+HFOgAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.plot(freqs*1e-9, extinction_modes.imag)\n", "plt.xlim(0, freqs[-1]*1e-9)\n", "plt.xlabel('f (GHz)')\n", "plt.ylabel('$W_m-W_e$')\n", "plt.title('Reactive stored energy')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compare model with exact solution \n", " \n", "The accuracy of this modal decomposition can be checked by comparison with the direct numerical calculation" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEZCAYAAACXRVJOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3zV9fX48de5NwkZhB32CKiAEMIwoAZkCAIuxFWL1NVStELVX+vWKtpvW2tb2yqtSq2jlrqw4F4oiDiQgGgVUFFBNgkjJKzc3JzfH/eTcAkZ997c3Jt8cp6PRx7c8Rnnkxs+5763qCrGGGOMJ94BGGOMaRgsIRhjjAEsIRhjjHFYQjDGGANYQjDGGOOwhGCMMQawhGBiTESSRURFpGsMznWXiMyuh+OOE5FPo33cEM7bRUQ+EJEiEfmNiHhEZK6I7BGRJaHGJSI/EZGXYhGzaVzExiEYESkOepoKHAL8zvMrVXVuDftOBGar6rEhnisZOAB0U9VNEYZc5zjCOG69xBthLL8Beqrqxc7z04CHgf6qeiCesRl3SIh3ACb+VLV5+WMRWQ9MU9WF8YvIVKMHsLrS828tGZhosSojUysRSRGRv4nIVhHZJCJ/EJFEEWkLzAd6iUix89NWRIaLyDIRKRSRLSLyZxEJ6cuHiLQRkX+JyDYR2Sgid4qIx3nvMRGZG7TtX0XklRriuEdEHnG27SsipSJyhXMN+SJyQ9CxEpxzfSsie0VkuYh0BJY4m3zpHHeyiEwUkXVB+w4QkfecqpvPROT0oPeeFpG/iMgbTlXP+yLSo4brP8X53e0RkZUiMtx5/SngIuBXThx3ArOB0c7zW6uIK1NEXhCRAufnT87rV4nIwqDtskTkHRHZLSJrRGRyqPGLyMCgfbeJyC9FpLuI7BORFkHbDXf+Fuye05Cpqv3YT8UPsB4YV+m1e4H3gHZAB2A5cJvz3kRgXaXthwFDAS9wDLAOuMp5LxlQoGs1538NeIBA1VUn4BPgMue9dOA74IfAWGAH0LGGOO4BHnEe93XO+zcnhqFACdDLef9XzrmOJfBFaTDQqqp4g8/lvL8B+CWQCEwAiglU7QA87cQ5xHl/HvB4NdeeCewExjkxnAHkA62DjnV70PZXAQuriSsRWOP8DlKBFCC38n5AC2ArMNX5vIYCu4Bja4sfaO3ENxNo5hxrqPPeO8AVQbE9CPwh3n/f9lPzj2VrE4qpwJ2qWqCq24H/Ay6pbmNV/VhVl6uqX1W/AR4BRtV2Eueb50jgF6q6X1W3AvcTSACoahFwKYFvxk8QSDLbwryWO1X1oKouB9YC2c7r04CbVXWdqpap6iequieE453i/HufqvpU9Q3gLQLf5ss9q6orVdUH/AcYVM2xLgP+q6oLnRheJVBFND7MawQYQeAGfavzuzygqh9Usd25wOeqOtf5vJYDLwHnhxD/ZAIJaLaqHlLVvc7+EPh8fgQgIknAD4AnI7gOE0PWhmBqJCICdCTwLbjcBqBLDfv0A/5E4FtlCoG/s/dDOF0PAt+48wOnBQLflNcFbbOUwDfaFALVROHwq2pB0PP9QHPnGrsA34R5PIDOwPeqGtw7o/LvJzhp7QeaU7UewBQRuTDotUTnHOHqBnynqmW1bNcDGCkiwckvAdgd9Ly6+LtR/e/seeABEelCoNSxSVU/CzV4Ex9WQjA1cm502wjcOMp1BzaXb1LFbv8AVgLHqGoL4G5Aqtiuso0Eqltaq2or56eFqg4J2uYXgA/YC1wXHGoo11MV5xo3E6jeOurtWnbfQuD3ESz49xOOjQSquFoF/aSp6p8jPFZmCHX2G4E3K52zuapeV8t+5ftW9TtDVYsJJOyLCZQmrXTQCFhCMKF4CrjTaahtD9wG/Nt5bzvQXkSCv/WmA4WqWiwi/YGfhnISVf0O+Ai4V0TSJdDP/jgRGQGBxk/gdgJVET8C7nBKI9XFEY5HgN+KSC8JGCwirVT1EFAI9Kpmv/cAj4hc5zRMn0agiue5CGJ4ArhQRMaKiFcCjfljncbtcC0FioBfi0iqc6zcKrZbAAwWkYsk0FEgSUROEpHeIZxjAXCsiPzM2a+FiAwNev9fBKriJgLVdl02DYclBBOKOwjUZX8BrCJQ/XOv896nwIvABqdnTBvg/wHTJDC+4W/AM2GcawqBxty1BBo3nwE6OPXQ/wbuUtXVqrqaQMnjSRFJrCaOcNwDvEKgMXQv8BCBhtLy63/OOe6k4J1U9SBwFnABgQbh+4CLnLaTsKjqtwTq7u8CCghUPV1LBP9Pnfr+M4CBwCbge+C8KrbbTaAh/AoCVXFbCLQRJYZwjt3AaQTaeHYAXxJouyi3iEDV3lKnPcg0cDYwzRhTb0TkA+DvqvrvWjc2cWclBGNMvXDGUPQm0MBsGgHrZWSMiToReZpAVdQMtZHUjYZVGRljjAGsysgYY4yjwVcZtWvXTjMzM+MdhjHGNBorVqwoUNWMcPdr8AkhMzOTvLy8eIdhjDGNhohsqH2ro1mVkTHGGMASgjHGGIclBGOMMUAjaEMwpqnw+Xxs2rSJgwcPxjsU00gkJyfTtWtXEhNrnWkkJJYQjGkgNm3aRHp6OpmZmQRN/21MlVSVnTt3smnTJnr27BmVY8a8ykhEWonIPBFZ6yzXd3KsY2gSZrUM/JhG4+DBg7Rt29aSgQmJiNC2bduolijjUUL4K/C6ql7gzGCZGocYjGmQLBmYcET77yWmCcFZdHskcDmAqpYQWNfWREvlUkH581mFsY/FGNOoxLrKqBeBRbkfE5FPROQREUmrvJGITBeRPBHJy8/Pj3GIxjRdIsIllxxeLru0tJSMjAzOOuussI6TmZlJQUFBnbd5/PHHmTlzZljnbkimTZvG6tWrAfjtb397xHu5uVWtVxRfsU4ICQTW2X1QVQcD+4CbK2+kqnNUNUdVczIywh593bTNKjyyNFD5uTE1SEtL4/PPP+fAgcAEpW+99RZdulS7fLapxSOPPEK/foFF/SonhA8++CAeIdUo1glhE4HFtpc5z+cRSBDGmAbi9NNP55VXXgHgqaeeYsqUKRXv7dq1i8mTJ5Odnc1JJ53EZ599BsDOnTsZP348gwcP5sorryR4FuXJkydzwgkn0L9/f+bMmVPr+R977DF69+7NqFGjeP/99ytez8/P5/zzz2fo0KEMHTq04r3i4mKuuOIKBgwYQHZ2Ns8//3xF7AMGDCArK4ubbrqp4jjNmzfnpptu4oQTTmDcuHF8/PHHjB49ml69evHiiy8CgZLJOeecw8SJE+nTpw933XVXxf733XcfWVlZZGVl8Ze//AWAffv2ceaZZzJw4ECysrJ45pnAIoGjR48mLy+Pm2++mQMHDjBo0CCmTp1aEQcEegvdcMMNZGVlMWDAgIp9Fy9ezOjRo7ngggvo27cvU6dOpd5np1bVmP4QWIO2j/N4FvCHmrY/4YQT1ITvooc/0BdXbY53GCYMq1evjncImpaWpp9++qmef/75euDAAR04cKAuWrRIzzzzTFVVnTlzps6aNUtVVd9++20dOHCgqqr+/Oc/17vuuktVVV9++WUFND8/X1VVd+7cqaqq+/fv1/79+2tBQYGqqvbo0aNim3JbtmzRbt266Y4dO/TQoUOam5urM2bMUFXVKVOm6Hvvvaeqqhs2bNC+ffuqquqNN96o1157bcUxdu3apZs3b644js/n0zFjxuj8+fNVVRXQV199VVVVJ0+erKeddpqWlJToqlWrKq7nscce044dO2pBQUFF3MuXL9e8vDzNysrS4uJiLSoq0n79+unKlSt13rx5Om3atIoY9uzZo6qqo0aN0uXLl1f8biv/rlVV582bp+PGjdPS0lLdtm2bduvWTbds2aKLFi3SFi1a6MaNG9Xv9+tJJ51Ucf3Bqvq7AfI0gvtzPHoZ/RyY6/Qw+pbAWq4mynx+xecvi3cYpg4yb34l6sdcf8+ZtW6TnZ3N+vXreeqppzjjjDOOeG/p0qUV38BPPfVUdu7cSWFhIUuWLOG///0vAGeeeSatW7eu2Of+++9n/vz5AGzcuJGvv/6atm3bVnnuZcuWMXr0aMqrii+66CK++uorABYuXFhRHw+wd+9eioqKWLhwIU8//XTF661bt2bJkiVHHGfq1KksWbKEyZMnk5SUxMSJEwEYMGAAzZo1IzExkQEDBrB+/fqK45x22mkVcZ533nksXboUEeHcc88lLS2t4vX33nuPiRMncv3113PTTTdx1llnccopp9T6ew7+nU6ZMgWv10uHDh0YNWoUy5cvp0WLFgwbNoyuXbsCMGjQINavX8+IESNqOWLkYp4QVHUVkBPr8zY1Pn8ZpX5b/KgxC+XmXV8mTZrE9ddfz+LFi9m5c2fF61pFlUV518equkAuXryYhQsX8uGHH5Kamsro0aNr7TdfXVfKsrIyPvzwQ1JSUo54XVWP2qeqOMslJiZWbO/xeGjWrFnF49LS0mrjEJFqj9u7d29WrFjBq6++yi233ML48eO54447qo0h1FjLYwPwer1HxFcfbC4jlyopLaPESggmQj/+8Y+54447GDBgwBGvjxw5krlz5wKBm327du1o0aLFEa+/9tpr7N69G4DCwkJat25Namoqa9eu5aOPPqrxvCeeeGJFEvL5fDz33HMV740fP57Zs2dXPF+1alWVr+/evZsTTzyRd999l4KCAvx+P0899RSjRo0K63fw1ltvsWvXLg4cOMCCBQsYPnw4I0eOZMGCBezfv599+/Yxf/58TjnlFLZs2UJqaio/+tGPuP7661m5cuVRx0tMTMTn8x31+siRI3nmmWfw+/3k5+ezZMkShg0bFlas0WJTV7iUz19mVUYmYl27duXaa6896vVZs2ZxxRVXkJ2dTWpqKk888QQAd955J1OmTGHIkCGMGjWK7t27AzBx4kQeeughsrOz6dOnDyeddFKN5+3UqROzZs3i5JNPplOnTgwZMgS/3w8Eqp5mzJhBdnY2paWljBw5koceeojbb7+dGTNmkJWVhdfr5c477+S8887jd7/7HWPGjEFVOeOMMzjnnHPC+h2MGDGCSy65hHXr1nHxxReTkxOo2Lj88ssrbtjTpk1j8ODBvPHGG9xwww14PB4SExN58MEHjzre9OnTyc7OZsiQIRXJE+Dcc8/lww8/ZODAgYgI9957Lx07dmTt2rVhxRsNDX5N5ZycHLUFcsI36g+L+NGJPfjpyF7xDsWEaM2aNRx//PHxDsMQ6GWUl5d3RMmjoarq70ZEVqhq2FXzVmXkUj6rMjLGhMmqjFyqxK/WqGxMhC6//HIuv/zyeIcRc1ZCcClrQzDGhMsSgkv5/GX4yiwhGGNCZwnBpXz+MnylVmVkjAmdJQQXUlUbqWyMCZslBBfyOY3JpVZlZMLk9XoZNGgQ/fv3Z+DAgdx3332UOX9HeXl5XHPNNVE5z+OPP86WLVuicqzySeLCtXjx4lqn9V61ahWvvvpqxfMXX3yRe+65J6LzNQbWy8iFyksGJVZl5H5RXgApJSWlYgTwjh07uPjiiyksLOSuu+4iJyenYnBWsNLSUhISwruVPP7442RlZdG5c+eoxF1fVq1aRV5eXsWcTpMmTWLSpElxjqr+WAnBhcoTgpUQTF20b9+eOXPmMHv2bFT1iG/Us2bNYvr06YwfP55LL70Uv9/PDTfcwNChQ8nOzubhhx+uOM69997LgAEDGDhwIDfffDPz5s0jLy+PqVOnMmjQoIq1F8qtW7eOcePGMXDgQIYMGcI333xDcXExY8eOZciQIQwYMIAXXnihypgrnwsOT0ENUFBQQGZm5lH7ffzxx+Tm5jJ48GByc3P58ssvKSkp4Y477uCZZ55h0KBBPPPMM0cs2LNhwwbGjh1LdnY2Y8eO5fvvvwcCXVavueYacnNz6dWrF/PmzavbBxFDVkJwofIBadaG4GIxWiq1V69elJWVsWPHjqPeW7FiBUuXLiUlJYU5c+bQsmVLli9fzqFDhxg+fDjjx49n7dq1LFiwgGXLlpGamsquXbto06YNs2fP5o9//GOVJY6pU6dy8803c+6553Lw4EHKyspISkpi/vz5tGjRgoKCAk466SQmTZp0xAR0r7322lHnClXfvn1ZsmQJCQkJLFy4kFtvvZXnn3+eu++++4gRy48//njFPjNnzuTSSy/lsssu49FHH+Waa65hwYIFAGzdupWlS5eydu1aJk2axAUXXBByLPFkCcGFytsQfDYwzURBddPbTJo0qWLm0TfffJPPPvus4ttwYWEhX3/9NQsXLuSKK64gNTUVgDZt2tR4rqKiIjZv3sy5554LQHJyMgA+n49bb72VJUuW4PF42Lx5M9u3b6djx44V+4Z7rmCFhYVcdtllfP3114hIlZPQVfbhhx9WTPl9ySWXcOONN1a8N3nyZDweD/369WP79u0hxxFvlhBcyFdqJQTXKy8J1FPJoNy3336L1+ulffv2rFmz5oj3ytcEgEDSeOCBB5gwYcIR27z++uvVTmddleqSz9y5c8nPz2fFihUkJiaSmZl51DTaVU2DDZCQkFDRMF7d1Nu/+tWvGDNmDPPnz2f9+vWMHj065JjLBZ87eNrqhj5fXDBrQ3Ch8rYDSwimLvLz87nqqquYOXNmrTf1CRMm8OCDD1Z8s/7qq6/Yt28f48eP59FHH2X//v0AFdU46enpFBUVHXWcFi1a0LVr14qql0OHDrF//34KCwtp3749iYmJLFq0iA0bNhy1b3XnyszMZMWKFQDV1ucXFhZWrB0dXC1UXZwAubm5FQvzzJ07t14XrokVSwguVN67yKqMmoBZhVEtHZSv+9u/f3/GjRvH+PHjufPOO2vdb9q0afTr148hQ4aQlZXFlVdeSWlpKRMnTmTSpEnk5OQwaNAg/vjHPwKBhterrrqqykblJ598kvvvv5/s7Gxyc3PZtm0bU6dOJS8vj5ycHObOnUvfvn2PiqG6c11//fU8+OCD5ObmUlBQUGX8N954I7fccgvDhw+vmG4bYMyYMaxevbqiUTnY/fffz2OPPUZ2djZPPvkkf/3rX2v9PTV0Nv21C326cQ/n/O19Bndvxfyrh8c7HBMim/7aRMKmvzY18vnLSPJ6bLZTY0xYLCG4UIm/jJQkr7UhGGPCYgnBhXx+Jc0SQqPU0KtwTcMS7b8XSwgu5CstLyHYzaUxSU5OZufOnZYUTEhUlZ07d1aM1YgGG4fgQj5/GWnNEigoOhTvUEwYunbtyqZNm8jPz493KKaRSE5OpmvXrlE7niUEFyrxl5GS6KXESgiNSmJiIj179ox3GKYJi3lCEJH1QBHgB0oj6RplaubzK2nNEqwNwRgTlniVEMaoatUjREyd+ZxeRqWWEIwxYbBGZRfy+ctITbRGZWNMeOKREBR4U0RWiMj0qjYQkekikiciedbAFj6fX0lN8uIrK7MeK8aYkMUjIQxX1SHA6cAMERlZeQNVnaOqOaqak5GREfsIGzmfv4zkRC8C+MssIRhjQhPzhKCqW5x/dwDzgWGxjsHtfKVlJHo9JHo9lFpCMMaEKKYJQUTSRCS9/DEwHvg8ljE0BT7/4YRQYg3LxpgQxbqXUQdgvjO3egLwH1V9PcYxuF6JX2mRJCR6pWKxHGOMqU1ME4KqfgsMjOU5m6Ly2U6tysgYEw7rdupCR1QZWQnBGBMiSwgudDghiJUQjDEhs4TgQiWlSqJXSPB6bPoKY0zILCG4kM9fRlJCoMrIEoIxJlSWEFzI5y8jwROoMrLpK4wxobKE4EKBNgSxEoIxJiyWEFzI51cSEzwkeMQSgjEmZJYQXKh8HEJSgseqjIwxIbOE4ELB4xBsTQRjTKgsIbhQid/pdmpVRsaYMFhCcKGK2U6tysgYEwZLCC5UMQ7BSgjGmDBYQnChI9sQrIRgjAmNJQQX8vkPT11h6yEYY0JlCcGFSsq7nXqtysgYEzpLCC5kVUbGmEhYQnAhX2kZCVZlZIwJkyUEFwq0IQSqjKyEYIwJlSUEl1FVfGWBKiNbD8EYEw5LCC7jL1M8Ing9zmynZZYQjDGhsYTgMuVdToHAegilVmVkjAmNJQSXKXF6GAG2HoIxJiyWEFymfOprgASvUGpVRsaYEFlCcBlfpRJCiVUZGWNCFJeEICJeEflERF6Ox/ndzFeqJCYE2hCSvB4rIRhjQhavEsK1wJo4ndvVgtsQEmzqCmNMGGKeEESkK3Am8Eisz90UBLchBBqVrcrIGBOaeJQQ/gLcCFT71VVEpotInojk5efnxy4yF/D5A9NWgNPt1EoIxpgQxTQhiMhZwA5VXVHTdqo6R1VzVDUnIyMjRtG5Q+VGZZu6whgTqliXEIYDk0RkPfA0cKqI/DvGMbhaSakebkPw2OR2xpjQxTQhqOotqtpVVTOBHwLvqOqPYhmD25WWHW5DSEqwKiNjTOhsHILLBKqMytsQrMrIGBO6hHidWFUXA4vjdX63qlxlZCUEY0yorITgMj5/GYkJVmVkjAmfJQSXOWIuI4+NQzDGhM4Sgssc0YaQ4KHUSgjGmBBZQnCZEv/hNoREj1BiJQRjTIgsIbiMr9TWQzDGRMYSgsv4/GUkJQSth2AJwRgTIksILuPzl5HgOTwOwRqVjTGhsoTgMke0IXg9+MrKULWkYIypnSUElwmuMvJ6BAH8ZZYQjDG1s4TgMoFGZal4nuj1UGoJwRgTAksILlNadrjKCJx1la1h2RgTAksILhO8hCY4i+SUWkIwxtTOEoLL+EoPT10BkGBVRsaYEFlCcJnA5HaH2xCSvB5KrIRgjAmBJQSX8fkrtyGIlRCMMSGxhOAyldsQEmz6CmNMiCwhuEzw9Ndg8xkZY0IXUUIQEVsHuYHyVdXLyKavMMaEIKyEICLdRSQT+HGl1++LYkymDnylSkLlgWlWQjDGhCDcNZV7EUgGg0XkbeBL5+e0aAdmInNUG4JHbGCaMSYkYSUEVV0MLBaRs4GXgWOBfsAPoh+aiUTlNoSkBJvx1BgTmkgblU8HUNWvVfUFYEP0QjJ1UXkcQoLH1kQwxoQm0oTwOjBPRDqJyA3A21GMydRB6VHjEKyXkTEmNJEmhDeAAuB7QIARUYvI1ElJ5W6nVmXUuMxqGfip/NiYGIg0IbwLLAMGEmhQ7hPKTiKSLCIfi8inIvKFiNwV4flNNY7qduoRKyE0dKHc+N2QHKpLdpWvzQ3X2kiF28uo3Cmq6gMQkR8CTwBnhbDfIeBUVS0WkURgqYi8pqofRRiHqSQwdUXlbqdWQmhQKm6EhVW/XtPjqo5R3fEagrrG1pCvzYUiTQi5TtvBHuBS4OpQdtLAWo7FztNE58fuVlHkKy0jMeHIqSus22kDUNWNLZJvwaGUJCqfp76FkpxCSXaVr62655Yc6k2kVUZ/AKYAnVW1DJgT6o4i4hWRVcAO4C1VXVbFNtNFJE9E8vLz8yMMsWmq3IaQ5LUqowYl1CQQ6k2vthtqNNsjwjletKt9aro2EzWRJoQiVS0Kep4U6o6q6lfVQUBXYJiIZFWxzRxVzVHVnIyMjAhDbJoqtyEkWJVRfIVz46ouCUTyjTicNolQ6/ZrO08k11n5cSjXWl0VmiWIOou0yugxEZkLtBKRXAJVP2FR1T0ishiYCHweYRwmiN+Z5trrObINwaqM4qTWm2nh0dtVd7Osbp/gappQ4ohGNVVExwghzkj2sSQQVRElBFX9t4h8Apzr/Py4ll0AEJEMwOckgxRgHPD7SGIwR/P5y0jwHlnoS/KKlRDirXLd+lE32FC+FYe4Tdg33TAassM5fzjJ7qj2hmoSYSgxVXc+E5JaE4KIXAHMADKBTcA7wB9V9QvgizDP1wl4QkS8BKqrnlXVl8M8hqlG5fYDcBqVbcW02Kq1cTQKN6pQShLhlCLCPXeox67rtdYlOZiw1ZgQRGQ68GvgPuAboAuBEsGXInKBqr4RzslU9TNgcISxmlr4SsuO6HIKgSqjfSWlcYqoCarPm2NI56/lHDVVOdX0OHifms5Zn9cYShXaEe9ZaSFctZUQrgHOVdUPgl77q4icBswVkZHANiBLVZfWV5AmNJWXzwRnPYRSqzKKuVrGBxwq9fP55r2s3LCbzzYXsveAj4M+P4dKy2iZkkjnVil0aZXMMRnNGdKjNR1aJEcWQ1WPIxWrG38o4n1+l6otIXSrlAwAUNW3ROQXwL8JVAPdDVhCiLPKPYzAGZhWZlVG9S7Ehtd1O4p55L1veenTLWS2S+OEHq0Z3TuDNmlJNEv00CzBy94DPjbtOcDm3QeYt2ITt8z/H2lJCQzNbM3oPu0Z6Wxft3hDqM9vTDfd2sZ3WGkhJLUlhHwR6aaqG6t471kCI5QvUtV50Q/NhMvnLyMpoXIbgo1DiCvn5vPV9iLuff1LVm3czY9O6sGSG8fQtnmzkA6hqnxXsI+Pvt3FK//byq8WfM5xHZozoX9HJvTvSGa7tPq8AtOE1JYQngJ+C1xSxXteoMCSQcNRedoKcLqdWpVR/amunt15rqo8+dEG/rLwa2aOOZYHpgwmJckb1ilEhF4ZzemV0ZyLT+zOoVI/y77dxRtfbOOChz6kXfMkJmZ1ZGJWR/p0SEdEaj+om4VaWrCSwlFqSwi/IzDf0KvALar6adB7txGY4M40EFVXGYlVGcXJzuJD3DjvM/KLD/H8z3LpGaVv8s0SvIzsncHI3hncfU4WKzbs5o0vtvGTx/NISvAwvn8HxvbtwJDurY7qhmyCzGppSaGSGhOCqu4XkVOB+4EVIvI9ga6nPZxNxtZzfCYMlZfPBFsPod7UMs/O5j0HmDLnI07P6siDPzrhqKq8aPF6hGE92zCsZxtuP/N4/re5kIWrtzPrxS/YWniAkb0zGH5MO3KPbUvX1qn1EkODVlVpwbqtVqvWcQiquge4VERuJzCQLAP4DnhFVffVc3wmDL7So8chBBKCVRnF0qbd+5nyj4+47ORMpp3SK2bnFRGyu7Yiu2srfjG+D1v2HODdr/J5b10B976xlpQkL4O7tSa7a0sGdmtF7w7ptEwJe5KBaqkq+0r85BcdYte+Q+wsLmH3/hKKDpay75Cf/SWllDqj6QVoluihdWoSrVKTyEhvRq92aXRplYLHU89VXjXO/9S0S8byqngAABqHSURBVAwhj1RW1e+BR+sxFlNHPr8esXwmON1OrYQQfdW0GWzcFUgGPx7ekx+P6Bmn4AI6t0phyrDuTBnWHVVl3Y5iVm3cw2ebCnnp0y2s21FMcqKXnu3S6NI6hXbNm9G2eRKtU5NIdno8JXo9+MvK8PkVn7+MfYdK2XuwlKKDpezZX0JBcQk7nZt/ftEhADLSm9EmLYm2aUm0TksiPTmB5s0SaJWaRKJXUOf7yUGfny17DvLFlr1s33uQb/KLKTpYyrHtm5PTow3Dj23LsJ5tSE+OUtKykkGtIp3LyDRA1XY7tRJCdFXTKLlrXwkXP/IRPz2lF5flZsY+rhqICMd1SOe4DulcmNMNCHyjzy86xDf5+9haeICC4kMUFJewvmAfJaVlHCoto6S0DK9HSPR6SPAKac0SSE9OoEVyIt3bpNImLYl2zZNo27wZ7dObkdasbreUvQd9fLWtiGXf7eLR97/jmqc+YUiP1pw7uAsT+nes8/FrH2TXtBucLSG4SIm/jARPpW6nHpvcrl45N45Sfxkz/7OSMwd0bnDJoDoiQvsWybSPZNBbPWmRnEhOZhtyMtswY8yxHCjxs3DNduZ/spk7X/yC07M6Mn1kL45tn16/gTTRBmdLCC4SGIdwZJVRUoJVGUVNDQ3Jv3ttLQleDzdMCGk1WROilCQvZw/szNkDO1NQfIj/LPueH875iEHdWnP1mGMY0r11ZAeubgGfJl6tZAnBRaqqMkrwWJVRffvvyk0sXLOdF2YMP2LqcRNd7Zo345qxx/HTU3oxb8VGZsxdyaBurbj59L70aBulwXnVzv7aNEoL1knZRXylVc1lZN1Oo6byAi6zCvn6Z5v4v1fWMOeSHFql1nE6CROSlCQvl5ycyTu/HE1Wl5ZM/tv7/Prl1RQd9EV+0CZyw6+NJQQX8ZUdXUKwKqP64/OX8cvnPuX68X3o07Ge67TNUVKSvMwYcyxv/r9RFB30MeHPS3h7zfbID1jTim1NZEU2SwguEhiHcGSVRYLHxiFExRFLSwZuHA8u/oZWqUlMGdYtvrE1cRnpzbj3goH88cKB3P3yamb8ZyU7iw/Vz8lcnhQsIbhIldNfJ3gotRJC1H2xpZAnPljP788fYHMHNRC5x7bjjetG0qVVCqf/9T3e/So/sgOVlxQa68yvdWCNyi5S4i8jsdIUCYkeocRKCJGrpmfRL1u9xC1nHE+nlilxCMpUJznRy61nHM/o3hn88rlPmZjVkZsm9iU5MbwJBY/SREY3WwnBRaobmGZtCNHXpVUK5w/pEu8wTDVyj23Ha9eewrbCg5z/4Ads2BnhLDsuutmHwhKCi/j8VbQheMWqjOqiUtXBxmu2MsjzHLMm9beqogauVWoSf586hB/kdOO8v3/AG19si+xANTU2w+H2JRe0L1iVkYv4/EqLZJvcrj7d9dJqpo3oSbc2TXDm0EZIRLgsN5MBXVvy8/98wsoNu7lhQh+bFrwa9ltxkZLSaqqMyspQtaQQtko9i96+6Cu+yS/mpyNjN4OpiY4h3Vvz0s9H8PmWQq54fDl79peEf5DgxmaXdk+1hOAigTaEI6sxvB5BAH+ZJYS6OOjzM+ulL5g1qT/NEurYQGniok1aEk9cMYy+HdOZNPt91m7bW38na6RJwaqMXKSqRmVwZjwtU+w+FqIqehYlA/2PeYdRvTPiEpKJjgSvh9vO7Ef/zi25+B/L+N15A5jQv2NkB3Ph4juWEFykqnEI4Kyr7C+re9e7Ju62M4+PdwgmSiYP7kLPdmlc+eQK1u0o5urRx0Snk0AjnwspplVGItJNRBaJyBoR+UJEro3l+d2uqnEI4KyrbA3LoatUR3x9/yXcc+Iya0h2mYHdWrFgxnDe+GIb1z2zioM+f+QHayQ3/NrEug2hFPilqh4PnATMEJF+MY7BtaqaugICxWQbixC5xV/mc/WYY+IdhqkHHVsm8+yVJ+MvU3445yN2FB2M/GChNDY38EbnmCYEVd2qqiudx0XAGsBG90RJaVnVVUZJXg8lpZYQalXpP6veuYeLu7zBteOOo0W0lnE0DU5yopcHpgxmTJ/2nPu3D1i9pR4bmxu4uLUhiEgmMBhYVsV704HpAN27d49pXI1ZdY3KCV6pWNzchO7tNTvYUXSIKUNt8jq3ExGuHXccx7RP45J/Bhqbx0ezsfmobRpm+0JcEoKINAeeB65T1aPSsarOAeYA5OTk2J0sRFWNQwCbvqJW1cxX9Lv0F7jtzONtEFMTclZ2Z7q1TuXKJ1fwTf4+rhrVq0mNSI95QhCRRALJYK6q/jfW53ezqpbQBEsIkWqfnsyYPu3jHYaJsYHdWjF/Ri7Tnshj3Y5ifnteVuRjT6pbqrO61+JcWoh1LyMB/gmsUdX7YnnupqD6bqdi01fUpFJj4P5bd3Ji0vPcfHrfJvXt0BzWqWUKz111MvtLSpn6j2UU1Nf6ClWJY8NzrMvCw4FLgFNFZJXzc0aMY3CtGgemWQkhZP987ztyMtswsFureIdi4ig1KYG/XTyE3GPacs7s9/liSxS+vYcy9UUcxbTKSFWXAvaVq568vudseJij/tgSPEKJJYSqVVoFbWfxIR69713mXz08vnGZBsHjEX4xvg/HdUjnkn9+zG8mZ3H6gE6xOXkc1l2wkcpNQFKCzXgaqgfeWcekgZ3JbJcW71BMA3L2wM5ktk3jyifzWLN1L9eN643HU8fvtqH0RooxSwhuUE0vmfI/uASPrYlwlGp+Zy94nuOtX4yKQ0CmoRvQtSUvzBzB1XNX8MW/8vjzDwdFf3xKTXMhxaDx2frTNQHWyyh0VwzvSbvmzeIdhmmgMtKbMXfaSXRulcLkv73Puh1F0TlwbYvwVLlP9BufrYTgBuV/SJVKBuVskZwqVPqdffqTDfz0X3ksPqVnHIMyjUFSgodfT87imeXf84OHP+L/JmdxRjTbFcLtqhpFlhBcotRfVu2HGeh2aiWEmtzz2lquHXccqUn2X8KE5qKh3enXqSU/m7uCT77fzU0T+8ZnEGMUq5Ksysgl9hzwMcQ7r8o/iEC3UyshVKi0EtqiKV+zveggF+XYFBUmPAO6tuSlmSP4ansxU/7xEVsLD0T/JLV1VT1q+5ac0MlzQiSnsoTgErv3ldA6teoGrgRnPQRzNH+Z8vvX1nLjhDh9uzONXuu0JB67fCij+7Tn7AfeZ9GXO+r/pPU0jsHKxy6xc18JbdKSqnwvyWu9jIAqexZ5gdT2rzKhf4e4hGTcweMRZow5lpwerbnumVVMGtSZX57Wh6Qq1iepk3oei2BfiVxidw0J4a5PhnP5W4NiHFHjcftZ/WyKChMVJ/Zqy8s/H8G67cWc/+AHfJtfXH8nC7cqKQSWEFxi1/7qE4JxVPrP85cRy5nZexFDureOY1DGbdo2b8Yjl+VwYU5XLnjoQ55dvhHVGLXh1TE5WJWRS+wqLqF1aqWEUMuAtabu8Q/W89LMEfEOw7iQiHDpyZmc1Kst1z69ijdXb+O35w2gfXpy/ZwwSv+nrYTgElZCqEWlnkU3Zi3hoqHdbJ1kU696d0jnhRnD6duxBWf8dSmvfLa1/k86q5AVW8tWRLKrJQSXqLINoXL9YpTrGxur1Vv28s7afGaMOTbeoZgmICnBw/UT+vCPS0/gT299yVVPrmDH3jqs3VyPLCG4xM59JbS2EsLRKg/vn9WSfnO6ce3YY22dZBNTg7u35tVrTuGY9mmc/tf3Ytu2ECJLCC6xe38JbSq3ITiWXfotF3Z8LcYRNWxThtla3Sb2khO93DChL//6yTD+9dF6pvzjo+jNhxQFlhBcYvc+X7VtCB1bJrO1sGEWUetdpWqyk5v9l48v+84GoZm46t+5JQuuHs6E/h258KEP+f3ra9lfUhrvsCwhuMWuGsYhdGiRzI69hygra1jF03pXxWyQw3q2YVjPNnEKyJjDErwerhjekzeuG8mm3QcY96d3efHTLXGtRrKE4AIHSvz4VUlNqnoh8OREL82TE9i5ryTGkTUc62duYbDnOW45/fh4h2LMEdq3SOaBKYP580WDeGjxN1z40If8b1N8On/YOAQX2LW/hLZpSTWOtu3YIplthQfJSG8Cc/1XMf4iE5g+eiUdW9ZTP3Bj6ujEXm156ecjeDZvIz9+Yjm5x7Tl+vF9Yto12koILhCY2K7mHkadWibXz0yMjchPRthaB6Zh83qEKcO6s/j60fRsl8bZs5dy90urKSg+FJPzW0JwgZomtivXsWUy2xpo3+eoq9SQPDTxeVb+eH30Jxozpp6kNUvgunG9eev/jcJfVsa4+97lntfWsrueq33tf4gL1DSxXblOTaGnUTVLCp45oJPNV2QapYz0Ztx1ThavXHMKhQd8jPnTYn7/+tp6KzFYQnCBmnoYlevYMoVtbk8IlXzwo28Ynjyf6yf0iXcoxtRJl1Yp/O68Abw0cwRFB32M/dO7zHrxC7bsiW41sDUqu8Cupt6GUM0kfremLeDuc/rTvJn9mRt36NYmlf+bPIBrTj2OR5Z+xxn3v8eo3hn89JReZHWp+zrLMS0hiMijIrJDRD6P5Xndbtf+Eto0D6ENoYmVEIb0aM3Y423hG+M+7Vskc+sZx7PkxjH079yCn/4rjx88/CGv/m9rndZPj3WV0ePAxBif0/V276t+2opyHVsE2hAa2twpdVbFlN7v/PArRqTMZ9ak/nEKypjYaJGcyPSRx7DkxjFcclIPHn9/Paf8flHEx4tpWVpVl4hIZizP2RTs2ldC67SaJ2pLa5ZAswQPe/b7XD8J3s3P/4/ZFw+xyetMk5Ho9XD2wM6cPbAzq7fspf9tkR2nQVauish0YDpA9+42CVltQmlUBujUMoWthQfdkRCqaTf4aa+3OS+juU1PYZqsfp1bRLxvg+xlpKpzVDVHVXMyMjLiHU6DtzvExXECYxFc2rDs2LLnAL84rXe8wzCmUWqQCcGErqxMA9VAtbQhgIvGIlTRbrDiiu/ISZjH36cOsQFoxkTI/uc0cnsP+khJ8pIYwnTObu5pNPM/n/D787Pp0TYt3qEY02jFtA1BRJ4CRgPtRGQTcKeq/jOWMbjNrn2Bie1C0allMsvX767niOpRNe0GF3V6nQt6trEupsbUUax7GU2J5fmagt37Q186s0OLRlxCqGJKinJJCR6uG2ftBsbUlVUZNXI7i2sfg1CuU8uUxj/BXVC7wSNjP2F8yxeZffEQvJ7qp/42xoSmQXY7NaELtYcRNNI2hGqqiQD+ufQ75v0sl5YpNt7AmGiwhNDI7aphLeXKWiQnUKZK0UEf6Q190FYNVUQAJ3jn8a/LcujSKiVGARnjflZl1Mjt2nco5DYEEWl8pYRKaxus+skGchLm8eeLBtG/c90n8zLGHGYlhEZu1z4fx3VID3n7d4rOgQc54ibboNRQRQQw7Ynl3HtBNiN724BFY6LNSgiN3O79oTcqN3g1VBOVlwzuvSCbU/ta91Jj6oOVEBq5wMR2ISSE6r55N8SSwqzCI+J77+t8rnt8OX+40JKBMfXJEkJjNqslC4Dv0rbEO5LIVVUqCHrt2eUbufeNL3nokhMYmmkT1hlTnywhuEBIvYzKSwINuWRQyZ9yP+aFRet45sqTOCajebzDMcb1LCE0RpW+Vbf8fTvn9dBv8nHtelpLl9I9N+bzy2c/Zfe6Av57dS7tmjeLUWDGNG3WqNzUzCpkRu9FzFuxKd6RVOvM+5fSs10aT08/2ZKBMTFkJYTGqI7VP1fkZnL9c59y2cmZeGIx5UMtJQKAktt38/C73/DEh+v57dn9GN+/Y/3HZYw5giWEJuiEHq1pnpzAu1/lM+bp4wIvRrtNIcxkNWn2Ujq1TOaFmSNs9LExcWJVRo3UoVI/gz3PsfGarWHvKyJcntuTxz5YH/3AKquhdLDrhnxmDfmAnITn+dnoY3j08qGWDIyJIyshNFKL1ubTu0M63dqkRrT/BS/154LgF4Jv3OGUFoJLAiFUDQUb+6fFnD2wM29cdwptra3AmLizhNBIPb9yE+cP6Vp/J6jqRl/5cZj23lzAU8u+57H31zO4eyuen9CHXtad1JgGwxJCY+PckD/iWe77wcA6HKeQ7XsP0uG+Kkb+hvJNv5Y5h4LPU/7eyHsXMap3Bo9clkNWF5uYzpiGxhJCI3Vq3/Z1HkfQoUVy7RsdUZUU/k38woc+4LuEeUwe1IWXh2fStXVkVVzGmPonqhrvGGqUk5OjeXl58Q4j/mr6Bl5H3+YXM+2JPN4pPifiY5TdsQfP3a0AmNL5Df63uZDcY9pywQldGdO3PYle679gTKyIyApVzQl3P/tf2hhE8M08HL0ymjP/6uEVz2cct6ji8cHbdgXFcTj5/G/a97x47pqK54N//VbF45+M6Mny28Yx59IcxvfvaMnAmEbCqowakV035NPmD846AFEeN9AyNRFmFbKz+BAj12yHrwOvZ896k6+cmqmsO9/gc2cc243Pf0bPdql8efIyBnRpxVs9WkF6IKZxUY3MGBMrVmXUUNVWKojR5HSqyt6DpRwo8ZPazEtaUoItaG9MAxdplZGVEOIp0u6cMZypVERomZJoC9kb0wRYQqiLUPvqh3Pjr6lk0AimrDbGNF4xb+0TkYki8qWIrBORm6N24Fktg2621TwOdbtQ969rvFU9NsaYOIlpCUFEvMDfgNOATcByEXlRVVfHMo6oqO6GXtONPpIbv5UKjDExEtNGZRE5GZilqhOc57cAqOrvqtun1kZlF3y7fvuirxj7TO/Ak0a4spkxpmFpLOMQugAbg55vcl47gohMF5E8EcnLz8+PWXCRWn75dxWP1119eOGZwpsK8P1q9+ENg2/uQY/HHl/V9BGFlgyMMTEV64RQVX/Fo4ooqjpHVXNUNScjI6PmI1a+cdb0ONTtwtw/ePH3Y9unVzxumZJY86Cs4GNaAjDGxFmsexltAroFPe8KbIlxDNETSiKp7T1jjGkgYt2GkAB8BYwFNgPLgYtV9Yvq9mmyA9OMMSZCjWJgmqqWishM4A3ACzxaUzIwxhgTOzEfmKaqrwKvxvq8xhhjambTUBpjjAEsIRhjjHFYQjDGGANYQjDGGONo8OshiEgR8GW846gn7YCCeAdRj+z6Gje7vsarj6qm177ZkRrD9NdfRtKftjEQkTy3XhvY9TV2dn2Nl4hENHjLqoyMMcYAlhCMMcY4GkNCmBPvAOqRm68N7PoaO7u+xiuia2vwjcrGGGNiozGUEIwxxsSAJQRjjDFAA04IIjJRRL4UkXUicnO844k2EVkvIv8TkVWRdhFrSETkURHZISKfB73WRkTeEpGvnX9bxzPGuqjm+maJyGbnM1wlImfEM8ZIiUg3EVkkImtE5AsRudZ53RWfXw3X55bPL1lEPhaRT53ru8t5PezPr0G2IYiIl8C6CacRWFRnOTBFVVfHNbAoEpH1QI6qumJgjIiMBIqBf6lqlvPavcAuVb3HSeqtVfWmeMYZqWqubxZQrKp/jGdsdSUinYBOqrpSRNKBFcBk4HJc8PnVcH0/wB2fnwBpqlosIonAUuBa4DzC/PwaaglhGLBOVb9V1RLgaeCcOMdkaqCqS4BdlV4+B3jCefwEgf+EjVI11+cKqrpVVVc6j4uANQTWOnfF51fD9bmCBhQ7TxOdHyWCz6+hJoQuwMag55tw0QfoUOBNEVkhItPjHUw96aCqWyHwnxJoH+d46sNMEfnMqVJqlFUqwUQkExgMLMOFn1+l6wOXfH4i4hWRVcAO4C1Vjejza6gJQap4reHVbdXNcFUdApwOzHCqJEzj8iBwDDAI2Ar8Kb7h1I2INAeeB65T1b3xjifaqrg+13x+qupX1UEE1qkfJiJZkRynoSaETUC3oOddgS1xiqVeqOoW598dwHwC1WRus92pvy2vx90R53iiSlW3O/8Ry4B/0Ig/Q6fu+Xlgrqr+13nZNZ9fVdfnps+vnKruARYDE4ng82uoCWE5cJyI9BSRJOCHwItxjilqRCTNadxCRNKA8cDnNe/VKL0IXOY8vgx4IY6xRF35fzbHuTTSz9BplPwnsEZV7wt6yxWfX3XX56LPL0NEWjmPU4BxwFoi+PwaZC8jAKcL2F8AL/Coqv4mziFFjYj0IlAqgMCMs/9p7NcnIk8BowlMKbwduBNYADwLdAe+By5U1UbZMFvN9Y0mUN2gwHrgyvI628ZEREYA7wH/A8qcl28lUM/e6D+/Gq5vCu74/LIJNBp7CXzJf1ZV7xaRtoT5+TXYhGCMMSa2GmqVkTHGmBizhGCMMQawhGCMMcZhCcEYYwxgCcEYY4zDEoJpMkSkj4h8IiJFInJNNdtcKSJ/iWFM94nIVbE6nzE1sYRgmpIbgcWqmq6q91d+0xkEeTvwh+DXROQOZyr2fc50ya+JyPigbdaLyLhKx7pcRJaGENMfgNuccxsTV5YQTFPSA/iihvfPAdaq6uag1+Y5r18KtAZ6An8FzoxGQM5AqLXApGgcz5i6sIRgmgQReQcYA8wWkWIR6V3FZqcD7wbtM47AmhznqOoyVS1xfl5X1WvDOPdFzjnLfw6JyOKgTRYTpQRjTF1YQjBNgqqeSmD6gpmq2lxVv6piswHAl0HPxwHLVHVTHc/9jHPO5kBn4FvgqaBN1gAD63IOY6IhId4BGNOAtAKKgp63A7aVPxGRNgRu5gI0U9XkoG0XiEhp0PMkYGXwwUXEA/yHQDvGw0FvFTnnNiaurIRgzGG7gfSg5zuBihkxVXWXqrYCTgCaVdp3sqq2Kv8Brq7i+L9xjl+5h1M6sKeuwRtTV5YQjDnsMyC4beFtYKiIdK3rgUXkhwRm17xAVX2V3j4e+LSu5zCmriwhGHPYq8Co8ieq+iawiEB10IlOF9RE4KRwDioig4EHCJQi8qvYZBTwWuRhGxMdlhCMOewloK+IdA567TzgZeDfBKp1vgOmEliRKlTnEOiyujSop9FrULFISz8Ca0cYE1e2HoIxQURkOtBPVa+L0fn+BHyjqn+PxfmMqYklBGOMMYBVGRljjHFYQjDGGANYQjDGGOOwhGCMMQawhGCMMcZhCcEYYwxgCcEYY4zj/wNcA9WWswGurwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.plot(freqs*1e-9, extinction.real/area)\n", "plt.plot(freqs*1e-9, np.sum(extinction_modes.real, axis=1)/area, '+')\n", "plt.legend(['Modal decomposition', 'Direct calculation'])\n", "plt.xlim(0, freqs[-1]*1e-9)\n", "plt.xlabel('f (GHz)')\n", "plt.ylabel('$Q_{ext}$')\n", "plt.title('Total extinction efficiency')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At low frequencies agreement is good, but the extinction is over-estimated at high frequencies. This discrepency is mostly due to some higher order modes which have been missed. \n", "\n", "In a later example, we will show how to obtain higher accuracy from the models. However, if we are using this is SRR in a metamaterial, we most likely only care about the frequency region of the fundamental mode, which is very accurately modelled." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.7.3" } }, "nbformat": 4, "nbformat_minor": 4 }