{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Evaluate the structure function of the phase screen" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "from astropy.table import Table\n", "import matplotlib.pyplot as plt\n", "import matplotlib\n", "import MegaScreen\n", "import numpy as np\n", "from theory import sf_integrated\n", "from MegaScreen import VonKarmanSpectrum, NestedSpectra\n", "\n", "%matplotlib inline\n", "matplotlib.rcParams[\"savefig.dpi\"] = 200\n", "matplotlib.rcParams[\"text.usetex\"] = False\n", "\n", "\n", "def multiscreen_sf(\n", " r0=5,\n", " L0=10000,\n", " diameter=128,\n", " decimate=20,\n", " numIter=10000,\n", " nfftWoofer=512,\n", " nfftTweeter=256,\n", "):\n", " \"\"\"Return the x and y structure functions of the woofer, tweeter, and summed screen\"\"\"\n", " args = locals()\n", " args[\"MegaScreenVersion\"] = MegaScreen.__version__\n", " sf = []\n", " # Use debug=True to get the woofer and tweeter screens\n", " for screens in MegaScreen.MegaScreen(\n", " r0=r0,\n", " L0=L0,\n", " windowShape=[diameter, diameter],\n", " dx=diameter,\n", " nfftWoofer=nfftWoofer,\n", " nfftTweeter=nfftTweeter,\n", " debug=True,\n", " numIter=numIter,\n", " ):\n", " sf.append([average_sf_xy(screens[i], decimate) for i in range(3)])\n", " sf = np.mean(sf, axis=0)\n", " r = np.arange(sf.shape[2]) + 1\n", " return r, sf, args\n", "\n", "\n", "def structure_function_brute_force(sig):\n", " \"\"\"Return the structure function over the first dimension of an ndarray\n", " \n", " This could be done faster using an FFT but it is simpler to code the brute force\n", " version.\n", " \"\"\"\n", " l = len(sig)\n", " a = [\n", " np.sum((sig[:-lag] - sig[lag:]) ** 2, axis=0) / (l - lag) for lag in range(1, l)\n", " ]\n", " return np.array(a)\n", "\n", "\n", "def average_sf(sig):\n", " \"\"\"Return the structure function in the first dimension, averaged over the last dimension\"\"\"\n", " a = structure_function_brute_force(sig)\n", " return np.mean(a, axis=-1)\n", "\n", "\n", "def average_sf_xy(screen, decimate):\n", " \"\"\"Return the average structure function of a 2-d screen in two dimensions\n", " \n", " Reduce the computation time by only averaging over every `decimate`th line\n", " \"\"\"\n", " return (\n", " average_sf(screen[:, ::decimate]),\n", " average_sf(screen.transpose()[:, ::decimate]),\n", " )" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This may take some time...\n", "..done\n" ] } ], "source": [ "print(\"This may take some time...\")\n", "r,sf,args=multiscreen_sf(numIter=1000)\n", "print(\"..done\")" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "lines_to_next_cell": 2 }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAADlCAYAAABEdm2gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd0FNX7x/H3TYeEGnqvIqAU6b0TEEFRsYBSUiiCiPIVEBQEBEQFaQImJBTpvWgkCb1KL1Kl95JQAgmpu/f3R4I/RJANKbO7eV7neI47OzP7AR+TZ+/MvaO01gghhBBC2DIHowMIIYQQQqSWNDRCCCGEsHnS0AghhBDC5klDI4QQQgibJw2NEEIIIWyeNDRCCCGEsHnS0AgAVJKZSqk7SqndRucRQgghUkIaGiunlPpCKfX7Y9tOPWXbe6n4qPpAC6CI1rpmKs4jBEqpqEf+MSulYh553SmDMlxXStXPiM8SAkApVV8ptUMpFamUuq2U2q6UqmF0rsxCGhrrtwWoq5RyBFBKFQScgaqPbSuTvO/zKg6c11pHp/RApZRTKj5X2CGttcfDf4CLQNtHts0zOt+zKKUclFLy81FYTCmVHfgVmAzkBgoDw4E4I3NlJvI/rPXbQ1IDUyX5dQNgI3DysW1ntNZXlVJ1lVJ7kr8h7FFK1X14IqVUIaXU6uRvDqeVUn7J232AGUCd5G/Qw5O3v6aUOqiUupv8raPSI+c6r5QaqJQ6DERLUyMspZTyUErFJv8CQCk1UikVp5TKkvz6e6XUt8n/nkUpNUEpdSl5xGWyUsr1kXO1V0odTq7RrUqpCsnblwD5gNDkmu6bvL2BUmpX8v77lVL1HjnXH0qpEUqpXcADoFBG/Z0Iu/ACgNZ6gdbapLWO0VqHaq0PK6W+VkrNfbijUqqEUko//LmplNqklPom+edslFJqjVLKUyk1Tyl1L/lneQlj/li2QxoaK6e1jgd2AQ2TNzUEtgLbHtu2RSmVG/gNmAR4AuOB35RSnsn7LQQuk/SD+m1gtFKqqdY6EOgJ7Ez+Bj1MKVUVCAJ6JJ/rZ2D1o79MgPeBNkBOrXVi2v/phT3SWkcBh0lqxAEakVSXtR95vTn538cDRYCXgXIk/dIYBKCUqg1MBbqRVKO/ACuVUk5a6w7ATaBlck1PSv6FsBIYQtI36C+T98/1SLwPgM5ANuB6mv7Bhb37CzAppWYrpVo/VleWeA/4kKSRndLATmAmSbV6HBiWlmHtkTQ0tmEz/9+8NCCpodn62LbNJDUXp7TWv2itE7XWC4ATQFulVFGgHjBQax2rtT5I0qhM56d8ZnfgZ631ruRvG7NJGjqt/cg+k7TWl7TWMWn3RxWZxGagUXKDXBaYlvw6G1AJ2J787dUH+ERrfVdrHQl8S9IPfkhqtqdorfcl16g/4ApUe8pndgGWa63Xaa3NWutg4BjQ8pF9ZmitT2qtE6RJFymhtb5H0r2IGggAwpNHxPNbeIqZWuszyXX+O0mj7uuS63AJUDVdgtsRaWhswxagfvIITF6t9SlgB0n31uQGXkrepxBw4bFjL5DU8RcCbmut7z/hvScpDvRPHpq/q5S6CxTln8Pwl1L55xKZ12agMVAL2AtsIGlkph7wZ/Ivh0IkXW49+kgNriTpUhIk1ejgx2o0L/9d0x88tn91pKZFGtFaH9dad9VaFyHp53IhYIKFh9945N9jnvDaI21S2i+578E27ARyAH7Adkj6NqCUupq87arW+lzy6+KPHVsMWAtcBXIrpbI90tQUA6485TMvAaO01qP+I5c8ql08r21AZZJGFTcDB4EXSRoteXi56RqQCJTWWt96wjkuAb9prcc95TMer89LJI3AfPwfuaSmRZrQWp9QSs0iaSRxP5D1kbcLGBLKzskIjQ1IvqSzF/iMpEtND21L3vZwdlMw8IJSqqNSykkp9S5QAfhVa32JpFGdMUopt+QbfH2AuTxZANBTKVVLJXFXSrVJviQgRKpore8CR4FewGattZmkGvcluaHRWieQdB/XRKVUnuQ6LKqUapF8Gn/gY6VU9eT3PJRS7ZRSD39x3ABKPfKxs4EOSqlmSinH5BuOmyml5JeLSDWl1ItKqf5KqSLJr4uSdJ/hHyQ17A2VUsWUUjmALwyMarekobEdm0kaat/2yLatydu2ACR/i30N6A/cAgYAr2mtI5L3fx8oQdJozQpgmNZ63ZM+TGu9l6TRnynAHeA00DUt/0Ai09sMKJK+vT587c4/a7wfSfW6F4gkabSxDIDWejvQl6Qb1u+SdFNmR/5/lGUUMCr58lIfrfVZ4C2SptJGkHTJ9RPk56BIG/dJuoS6SykVTVIjcwTor7UOAxaRdDP8PpKmd4s0prSWEVYhhBBC2Db5ZiKEEEIImycNjRBCCCFsnjQ0QgghhLB50tAIIYQQwuZJQyOEEEIIm2cXC+sppdoCbbNly+b3wgsvGB1H2Il9+/ZFaK3zZtTnSR2LtCY1LOyBpXVsV9O2q1evrvfu3Wt0DGEnlFL7tNbVM/pzpY5FWpEaFvbA0jqWS05CCCGEsHnS0AghhBDC5klDI4QQQgibJw2NEEIIIWyeNDRCCCGEsHnS0AghhBDC5ll1Q6OUKq+Umq6UWqqU6mV0HiFSSmpY2DqpYWErMryhUUoFKaVuKqWOPLa9lVLqpFLqtFJqEIDW+rjWuifwDlAvo7MK8SRSw8LWSQ0Le2TECM0soNWjG5RSjsBPQGugAvC+UqpC8nvtgN+A4IyNKTKD+7EJz3PYLKSGhW2bhdSwsBKrf/+N7bv3pPo8Gd7QaK23ALcf21wTOK21Pqu1jgcWAq8n779aa90a6JSxSYW9C9m5i29/asGx89dSdJzUsLAWWmuWn1rO/fj7KT1OalhYhVUr5uNxqAfR2/qn+lzW8iynwsClR15fBmoppRoDbwKuPOWbgVKqO9AdoFixYumbUtiNtdu3Mu9wDw55OvDavc3Ae6k95XPXMEgdi5TTWjNx/0QCjwRyO/Y2vi/7pvaUUsMiQwUv+hm38yP4X4Hc1MhfipapPJ+1NDRPpLXeBGx6xj7+gD8kPT8k/VMJWxeydR2zjvXleBYHhr7oS51KqW5mnsqSGk7eT+pYWMyszYzdPZb5J+bT4YUOeL/knW6fJTUs0kPYvHGYr4zji/yeVPSswHdNJ6b6nNYyy+kKUPSR10WSt1lEKdVWKeUfGRmZ5sGEfVm7bS3+xz/mhKsD37z8CW/X7pdWp05VDYPUsbCMyWxi+M7hzD8xn84VOvNV7a9wUGnyo1xqWGSILbO/JuraOAbn86RqvqoEtAoih2uOVJ/XWhqaPUBZpVRJpZQLSeP/qy09WGu9RmvdPUeO1P+FCPsVsuVXph37jPMujnxXZSBtq/ml5elTVcMgdSyeLcGcwBfbvmD5qeX0qNSD/1X/H0qptDq91LBIV9psZkfQAC5HBDA0ryd1CtZmWsufcXd2T5PzGzFtewGwEyinlLqslPLRWicCfYAQ4DiwWGt9NAXnlG8F4j+FbF3FpJMDuOrswA/VhtGyyofPfa70qOHk80odi6dKMCXw+ebP+f3c73zyyif0qdrnuZsZqWGR0bTZzO7ATzkZOZ8xeXLTuEhjJjf/iSxOWdLsM5TW9nOps3r16nrv3r1GxxBWJnTbSn44MZhIRwfG1fyG+hXfsOg4pdQ+rXX1dI73L1LH4nGxibF8tukztl7ZyqCag+hU3rLJRlLDwhpos5ld/r04EBvM1Fw58SruxZiGY3B2cLboeEvr2KpvChYitUJ2rOC7k1/ywMGBCTW/pU7F14yOJESKPEh4QN8Nfdl9fTdD6wylwwsdjI4khMW02cyu6d35IyGUwFw5aVeqHSPqjcDRwTHNP8ta7qFJFRnmFE8SumslY45/SZxSTKj9vdU3M1LH4nFR8VH0WteLPTf2MKr+KKtvZqSGxaPMJjO7pnqzOTGMwJw56PBCB0bWH5kuzQzYSUMjN6KJx4XsWc3Io0PQwITaP1CrQmujIz2T1LF4VGRcJD3CenA4/DBjG46lbem2Rkd6Jqlh8ZDZZOKPn7oQpjcxN0d2Or3YKS1n5D2RXTQ08q1APCpk/68M/3MwzhrG1/qB6hVaPfsgKyB1LB66G3sXv1A/jt0+xrjG42hVQmpY2A6TycQfUzoT7LCNxdmz4VPRm4E1B6bljLwnsouGRr4ViIdCDwcz7OAg3M1mvq/2HdVesv6RmYekjgXArZhbeId6c+buGSY1mUTTYk2NjmQxqWFhMpnYObkzq5x2siqbBz0r9eSTav3SvZkBuSlY2JGwY2v5ct8APM0mvqk6hmpV2hgdSYgUCX8Qjm+oL1ejrjKl2RTqFKpjdCQhLGYymdgxuQsrXP4gzN2dT6r2xbdSmq739Z/sYoRGhjnFupPr+GLX5+RNNDGi4nCqVbNsarY1kTrO3K5HX6dbSDeuRV9javOpNtnMSA1nXiaTia2TO7PUZSdh7ln5X/X/ZWgzA3bS0MgwZ+a27vQ6Bu74lMKJCQwrO5Aadd4xOtJzkTrOvK5GXaXb2m5ExETg38KfGgVqGB3puUgNZ04mk5ktU7qy2HUXG9yz8kXNL+hSsUuG55BLTsKmrT+7ns+3fUqphHg+L9qXmo27Gh1JiBS5fP8yPiE+3I+/T0CLAF7O+7LRkYSwmMlkZtOUrix22cmOrFkYWnsoHcoZs7yANDTCZq07t57Pt3xK2fh4+uXpSu3WvY2OJESKXLx3Ee8Qb2JNsQR4BVDRs6LRkYSwmMlkZsOUrixy3clutyyMqDuC9mXbG5bHLi45yXXbzCfsXBj/2/wpL8bH0i/r69R98wujI6Wa1HHmci7yHN3WdiPeFE9gy0C7aGakhjMPk8nM+inezHPdyW43N0Y1GG1oMwN20tDIddvMJeRcCP/b3J+KcbH0VQ2o22ms0ZHShNRx5nH27lm8Q7xJ1IkEegVSLnc5oyOlCanhzMFkMhP6kw9z3bZzwM3NahZ+tIuGRmQewWeDGbDlcyrHxdA3tiK1ffwhA9Y3ECKtnLpzim4h3QAI8gqibK6yBicSwnImk5mQqd2Z67qdw65Z+K7RD7Qu9arRsQC5h0bYkDVn1jBk2xCqxsTS734RKvdbiEqnZ4IIkR5O3j6JX6gfzg7OzPCaQckcJY2OJITFzCYzwdN6Md91C8dd3BjXZDzNijc3OtbfpKERNmHV6VV8tf0rXomJZeBtD8r0W4mDs6vRsYSw2InbJ/AN9cXV0ZUgryCKZy9udCQhLGY2mVkz7SPmu2zkpIsbPzaZQJPi1rWKtTQ0wuqtOLWCYTuG8cqDeIZFKPL1/hXnrHKNXtiOY7eO4RfqR1bnrAS1DKJo9qJGRxLCYmazZvXPHzPfZSOnXNyY1HQSDYs1NjrWv9jFPTRyZ739Wn5qOcN2DKPqAxNjb0aRo+tK3HMXMjpWupA6tk9HI47iG+qLh7MHM71m2nUzIzVsf8xmzYqfP2Gu0zpOu7gyqdkUq2xmwE4aGrmz3j4t+WsJw3YMo3KMYtLNm5jfmk/uouWNjpVupI7tz5GII/iF+pHdJTtBrYIokq2I0ZHSldSwfTGbNcv8P2O+YwhnXVyZ3OwnGhRtaHSsp7KLhkbYn8UnFzNi5wheinXF//oFbjedQsGXrPd/JCEedzj8cFIz45qdmV4zKexR2OhIQljMbNYsmTGA+Y7BnHNx5afmU6lXtIHRsf6TNDTC6iw8sZCRf4ykQnwOZl87xYUqgynZ4F2jYwlhsUPhh+ge1p1cbrmY1WoWBT0KGh1JCItprVkU+AUL1WouOrsytcU06hSpb3SsZ5KGRliVBScWMGrXKCqYCzL3yp+cKdaRF98YYHQsISx28OZBeoT1wNPNkyCvIAq4FzA6khAW01qzIOgrFrGCS86uTG0+jdqF6xkdyyLS0AirMe/4PEbvGk1FpzLMPr+bszkbUL7rFKNjCWGxAzcPSDMjbJbWmnmzhrPIvJTLySMztYrYRjMD0tAIKzHv+Dy+3f0tVdwrMe2vHdx0LU3ZngtAFs4TNuLAzQP0DOtJvqz5CPIKIr97fqMjCWExrTVz53zD4sSFXHV2YVqLadS0kZGZh+yioZGpgrZt7rG5fLv7W2rmrsU3f+7E5OhO3h4rcXDLZnS0DCV1bLv239j/dzMT6BWYaZsZqWHbpLXmlzljWBI/j2vOzkxtPo0aNtbMgJ00NDJV0HbNPTaXsXvGUr9AQ/ru309eFY1Dp8Vk8bTftTqeRurYNu2/sZ+e6/6/mcmXNZ/RkQwjNWx7tNbMmfs9S+PncM3ZmWnNplHDhi4zPcouGhphmx42M42LNOWdved52XyWO14/kbt0daOjCWGRh81M/qz5CfIKytTNjLA9Wmtmz/+RpbFBXHNyZlqzqVQvav2zmZ5GGhphiIfNTLNizWj6ZwJNErZxruoACtd52+hoQljk8WYmb9a8RkcSwmJaa+YsnMzSB/5cd3JmWrOfqG7l68w8izQ0IsM9bGaaF2tO/YvFaH93HqcKvU7p178wOpoQFpFmRti62YunsiRqKjecnJnWdDLVi9n+wqXycEqRoR5tZrwSmtL4bDfOe1SmrPcMUMroeEI8kzQzwtbNWjydZfemJDUzTSZRvXhjoyOlCWloRIaZd3ze35eZ3svTlRJL23HfKTdFeiwFJxej4wnxTAduHqDXul7SzAibNWvJDJZFTuS6kxPTmkykeokmRkdKM9LQiAzxcJ2ZZsWa0bfcQGL825BNxWLqshqn7HIjpbB+B28e/Mc6M9LMCFsza/kslt4dxw0nJ35q/CPVSzQ1OlKasup7aJRSbyilApRSi5RSLY3OI57P/OPz+Xb3tzQt2pShNUdxPsiP8pzlfpvpZCtW2eh46Upq2D48fJxB3qx5CfQKzFTNjNSwfZi94heW3Rqb1Mw0GkfNks2NjpTmMryhUUoFKaVuKqWOPLa9lVLqpFLqtFJqEIDWeqXW2g/oCcjTCW3QghMLGLN7DE2LNmVsg+/ZEDCEpolbuVT1f+Sv0d7oeM9FajhzOXjzID3X9SRv1rx2MzVbajhzmb1qPksjRnPdyZEpDb6jZin77EuNGKGZBbR6dINSyhH4CWgNVADeV0pVeGSXL5PfFzZk4YmFjN41msZFG/NDox9YvWgmb96ZxfmCrSnebojR8VJjFlLDmcKh8EP0XNcTTzdPAlva1aJ5s5AazhTm/LqEZTdHcs3Jkcn1xlCrTGujI6WbDG9otNZbgNuPba4JnNZan9VaxwMLgddVkrHA71rr/RmdVTy/xScXM2rXKBoXacz4RuPZsGUbr54ayjX3cpTwnmnTM5qkhjOHP8P/pGdYT3K75ba7xxlIDWcOc4OXs/TaUK46OTK57mhqv9DW6EjpylruoSkMXHrk9eXkbR8DzYG3lVI9n3SgUqq7UmqvUmpveHh4+icVz7TkryWM/GMkDYs0ZFzjcRw9fZHym7qT4JiVfH7LwDmL0RHTw3PXMEgdW5sjEUfoEdaDnK45M9NTs6WG7cjc31ey5PKXXHVyZGLtkdQp187oSOnOqmc5aa0nAZOesY8/4A9QvXp1nRG5xNMt+2sZI3aOoEHhBvzY+EciIuOIW9iFAuo2cR3X4JyriNERM5QlNZy8n9SxlTgacZTuod3J7pqdma1mZpZm5qmkhm3PgpBfWXppMJedHZhY62vqlbfN+xVTylpGaK4Ajz6NsEjyNovIE16tw4pTK/h659fUK1yPH5v8iMnkyG7/3tTSf3K7yXdkK1PX6IjpKVU1DFLH1uDYrWP4hfklNTNema6ZkRq2A/PD1rLowgAuuTjwY42vqF8h8zxOxloamj1AWaVUSaWUC/AesNrSg+UJr8ZbdXoVw3YMo26hukxsMhEXBxeWBY3ljdhVXCzbhQKNfIyOmN5SVcMgdWy0E7dP4Bfqh4ezB4FegRT0KGh0pIwmNWzjFq0PZcm5z7jk4sC4al/Q8KXMNSnNiGnbC4CdQDml1GWllI/WOhHoA4QAx4HFWuujKTinfCsw0Joza/hq+1fUKliLiU0m4uroyrI1K+lwbRyXc9ag2HvjjY6YptKjhpPPK3VskJO3T+Ib6ktW56wEeQVR2KOw0ZHSldSw/VmyYSOLTvfjgosDP1QdQOOXOxkdKcMpre3nUmf16tX13r17jY6Rqfx29jcGbxtMjfw1mNxsMlmcsrDtwJ+UWdkWJxc3PPttR7l7Gh3zuSil9mmtq2f050odZ6xTd07hHeKNi6MLs7xmUTR70WcfZCOkhjOHZZs3M//ER5x1VXxf6TOav+JtdKQ0ZWkdW8slp1SRbwXGWHtuLYO3DaZa/mp/NzOnr93CY6U3OdUDPLosttlmxghSxxnv9J3T+Ib64uLgwkyvmXbVzBhBajjjrdi6nQXJzczYl/vaXTOTEnbR0Mh124wXcj6EQVsHUSVvFaY0nUIWpyzci4nnZGAPqqi/iH51Mm5FKhkd06ZIHWess3fP4hPqg6NyJNArkGLZixkdyeZJDWeslTt2Mv9od067Kr6t8BEtq3U3OpKh7KKhERlr3YV1DNwykEp5KzG1+VSyOmfFZNasDBhJm8Qwrr78EZ41M9fNaMK2nIs8h0+oDwrFDK8ZlMhRwuhIQqTImj/2MP+wH6fcFKPLd8er5kdGRzKcXTQ0MsyZcTZc3MDnmz+nYp6KTG02FXdndwAWLV3I+7d+4nLehhRq/43BKW2T1HHGuHDvAj4hPpi1mUCvQErlKGV0JLshNZwxgvfsZ+7BbvzlBiNe8ObVWn2NjmQV7KKhkWHOjLH50mb6b+5Pec/yTG8+HQ8XDwDW/3GAFkcHEOlakMLec8DB0eCktknqOP1duncJ7xBvEs2JzGg5g9I5Sxsdya5IDae/tfsO8cveLpxwg6/LdKZd3c+MjmQ1rHqlYGE9tl3ZxqebPuWFXC8wvcV0srlkA+DEpZvk/d0HD4cEHL2XoLLkMjipEE92JeoKPqE+xJniCGwZSNlcZY2OJESKrDtwlDm7P+BoFs2wUh15o/4AoyNZFbsYoZFhzvS14+oOPtnwCWVylsG/hT/ZXbIDcCcqjrOzelBJnSGu7TRcClR4xpnEf5E6Tj/Xoq7hE+JDVEIUAS0CKJe7nNGR7JLUcPrZcPgEM3d25EgWzZASHXir4WCjI1kdu2hoZJgz/ey+tpu+G/pSIkcJ/Fv4k8M16e840WRm9YzhvGrawLUqfcn5SuZ4Vkh6kjpOH9ejr+Md4s29uHsEtAigvGd5oyPZLanh9LHlyGmCtr3LYXczg4q9wbuNhxkdySrJJSfxVHuv76XPhj4UzVaUgJYB5HTL+fd7C5csoNOdaVzJ35jC7YYbmFKIp7v54Ca+ob7cibuDfwt/KuapaHQkIVJk+/HzzNjyNofczXxeqA0dm8qki6exixEakfYO3DzAR+s/ooB7AQJaBpDbLfff74Xs3IfX8UHccSuSfBOwlJGwPhExEfiE+BD+IJzpzadTKa+siyRsy+5Tl/Df0J4D7iY+LdCSzi2+NTqSVbOL30Ry3TZtHQo/RK91vciXNR+BLQPJkyXP3+8dv3iTgmv9cHdIIEe3xeAmQ8tpReo47dyOvY1viC83HtxgavOpVMlXxehImYLUcNrZe/YK00LfYL9HIh/nbYK31zijI1k9u2ho5Lpt2jkacZReYb3I7ZabwJaB5M2a9+/37kTFcWZ2TyqpMyS0nYZLAbkXIS1JHaeNu7F38Qv140rUFX5q9hPV8lczOlKmITWcNg6fv87U4NfZ6xHPR54N6P7qJKMj2YQU30OjlHIHYrXWpnTIIwx0/NZx/ML8yO6ancCWgeR3z//3eyazZlXgSLqa1nO98scUkJuAhRWKjIuke1h3zkeeZ0qzKdQoUCNDPz82wcSNu9HcvXWN2MgIEqLvYo6PAVMcCjMODo44ODrh6OSMo7Mrjs5uOLm44uTsgrOjA84k4mKOwSXxPi4J93FOuIdjXCSOiTEoNDg6Y3Zw4UGiJjImgQfxGkdTLB6VXiXfi/Uy9M8q0seRS+FM/LUde7LF4ZurFr1em2p0JJvxzIZGKeUAvAd0AmoAcYCrUioC+A34WWt9Ol1TinT3152/6B7WHXdndwK9AinoUfAf789fupiOt6dyNV8DCr0+wqCUQjzd/fj79Ajrwem7p5nUdBJ1CtVJt8/SWnP5djSnju4j6uxuXCOOkCv6LAVN1yisblFcmdP08+K1I6BwwoSD0ngAHo+8/wc5pKGxAyeu3mLiqrbszhZD1+yv8Em7GUZHsimWjNBsBNYBXwBHtNZmAKVUbqAJMFYptUJrPTf9Yor0dObuGfxC/XBxdCGoZRCFPQr/4/11uw7idXQA91wLUsh7rtwELKxOVHwUPdf15OSdk0xoPIH6heun+WfEJpj44/Bxbu1fRe5rW6li+pOmKirpPVy56VaC6GzVOJWjCI45CuGcLQ/OWXPi5OaBg7MrGgdMJhNmUyKJCfGYEuJIjI/FlBiHOTGeRJMmASfilCsPHDy4rzyIUh48cPAgxuRAbKKZRJPG3clMgWzOFMmZBU93Z0yOrpTK4ZLmf16RsU5fv8v4ZW35I3s0ndxf5rM3ZhkdyeZY0tA0B0zAIK314Ycbtda3gWXAMqWUczrls4hSqi3QtkyZMkbGsEnnIs/hE5L8xOGWgRTNXvQf75+8Ek6eYF+yOcTi1C0YsuR8yplEakkdP58HCQ/ovb43xyKO8UPjH2hUtFGandtk1uw8fp6LW+dT6tpvNOQYDkpz2zEvEYWbEV22IfnL18Mt3wsUk0d+SA0/p3M37/H9ktfYmf0+72Z5kYFvzUMpZXQsm6O01pbtqNQOrXXddM6TKtWrV9d79+41OobNuHjvIt3WdiNRJzLTayalcv7zIX2RDxLYNL4TryeGEPlaADmqv2NQUmMopfZpratn9OdKHVsuJjGG3ut7s+/GPr5r+B1eJbzS5Ly3ouII27ieLAeDaJ64BXcVR4RLER6Ua0/B2h1wLlQJbOAXjtSw9bsYEcXo+W3YnuM2b7mWZti7K6SZeYyldZySm4IPK6WGASMfXnYStuvhc23izfEEeQX9q5kxmzXLZozGOzGE6y/1oEAma2aE9YszxdF3Q1/2Xt/Phnb5AAAgAElEQVTLmAZj0qSZuRARxYbfFvLimSDeczhKnHLlZonXcG7iR57itW2iiRG249rdB3w7vx3bc9ymnXMxhr2zXJqZVEhJQ5MbaAT0UkrtAg4Dh7XWS9IlmUg316Ov4xPiw4OEBwR6PfkhfQtXrKDTrUlczVOHQm+OMSClEE8Xb4qn38Z+7Lq2i5H1RtKmVJtUne98eBTrVs2m1sUAujmcI9IlLxHVB5OnoR9Fs+Z+9gmESKGbkTF8M6c9W3OE86pjYUa+txol9yemisUNjdb6HQCllCtQEXgZqAVIQ2NDbkTfwDvEm8i4SGa0nMGLuV/81z6b9h2l8eH+RLvkoaD3PJB7A4QVSTAl0H9Tf7Zd2cawOsN4vczrz32um/diWb1yIdVPT8LX4Qy3sxTmXqMfyVHzA3CSG21F+rgVFceI2W+xJcdVWjjkY0zH33CQn7OpZsm0baUfudFGax0H7E/+54n7COsUEROBb6gvt2Ju4d/yyc+1OXPjLu6rffFUUdA5FOXuaUBSIZ4s0ZzIwK0D2XR5E0NqDeHtF95+rvPExJtY+nsYxfePwVcd5K5rfu41Hk/u2l3AUR5xJ9LP3QfxDA96h805LtFEefL9+2ulmUkjFk3bVkotA1ZprS8+3KiUcgHqA11Imto9K10SijRxO/Y2fqF+3Hhwg2nNp1E5b+V/7RMVl8iBwL68rY5x22sKuYtWNSCpEE9mMpsYvHUwYRfCGFBjAO+9+F6Kz6G1JmzfCe4FD6ejaS1xju7cqfMVuRr3AWe3dEgtxP+7F5vA0MCObMxxlgY6J+M7heDoZOgkYbtiSUPTCvAGFiilSgJ3gSwkPTYhFJigtT6QfhGfTaYK/rfIuEi6h3bn0v1LTG029YlLwWutWRw0Hu/4VVwr15mCdT40IGnmJnX8dGZtZuiOofx+/nc+rfYpH1ZIeX2eD48idP543rrtT04VTcSLncj/+kiyyj0yaUZq+Omi4xIZOqMzG7OfpI45GxM/CMXJydXoWHbF4mnbAMnrzeQBYrTWd9Mt1XOSqYL/dj/+Pn6hfvx15y+mNJ1C3cJPnnm/+NffabunM3dzVqRg3zBwlG8NMuXVOpi1mRE7R7Ds1DJ6V+lNz8o9U3R8gsnM4rXrKbv7K2qq49zIWQXPDpNwKvzvUUp7IzVsHWLiTQwO8Ga9+z6qmz2Y/sF6XFzcjY5lM9Jj2jZa6wTgmlLKQylVFzhmjY2NSBKdEP336qkTm0x8ajOz88gpau/5mDin7BTwXSjNjLAaWmvG7BrDslPL8HvZL8XNzNHLEeyZO5z3Y+aT6JiFyGbjyV+nm6x2LTJMbIKJLwN7ssF9H1XNWZjWKUyamXTyPA+nLAlMJOlxCF2VUpu01vPTPJlIlQcJD/ho3UccjTjKuEbjaFik4RP3uxRxH/NSHwqqOyR2/BWVrUAGJxXiybTWfL/3exaeXEjXil35uOrHFh+bYDKz8Ne1vLLvC7o6nOda0dYUfG8SeORLx8RC/FN8opmvAvuyPstOXjK5Mr1jGK6u2YyOZbee53b+D4DSQDhwGfBWSoVprcPTNJl4brGJsfTd0JeD4QcZ23AszYo3e/J+CSa2B/bnPQ4R0XgseUqn38P8hEgJrTUT90/kl2O/0PHFjnxW7TOLFxw7fT2S7XO+4v3oucQ5Zye67UwKVnkznRML8U+JJjNDg/qzzm0zL5qcCXgvhCzy6Jh0leKGRms9UilVEJgDRJD0gMpPgcFpnE08h4cLju2+vptR9UfRqkSrJ+6ntWb+nGl4xyziasm3KdSoRwYnFeLpph+aTuCRQN5+4W0G1RxkUTOjtWblxp0U3dyPLuok14p4UbDjNJClB0QGM5k1w2YNJtQljNImJwLeCSarex6jY9k9ixsapVQRIBdwFvgM+AioAkzSWi9Nn3giJR4uOLb96naG1x1O29Jtn7rv6vWb6XDxG65nK0+hjj/Jku7Casz4cwZTD03l9dKv81XtryxqZu5Ex7N09gTevTEeJwdFpNdPFKzVSepaZDizWTNiznBCHH6lmMmRwA5ryJatoNGxMgVLFtYrASwH8gOxQAFgA/A/rfX49AwnLPf4gmNvln36EPuB05eosPUjcHQhn89iWX9DWI05R+cwcf9EWpdszfC6w3FQz755d+9fl7ix8GP8zBu5kbMyebvMIWvuEukfVojHaK0ZNe9bgvUSCpodmfnmCnJkL2p0rEzDklv9xwI/a60La61LAzmANUCwUkoWG7ACJrOJIduGEHYhjM+rf/6fC47dvBfD7Xl+lFLXUB1m4pCrWAYmFeLpFp5YyPd7v6dF8RaMrj8ax2esnmo2axasDibX3Ba0Mm/m5iufkL/vBhykmREG0FozduEEfk2cS16zAzNfX0yuXKWefaBIM5ZccnpBa/3uwxda60TAXyl1HhgGpNsKbEqpUsAQIIfW+vnWOLdzZm3m651fE3wumE9e+YTOFTs/dd8Ek5mQgCF8qHdyo/YQ8pd/8s3CIu1IDVtmxakVjNo1isZFGjO2wVicHP77R9PtqDhWBI2h060pxDlnI/7dFeR7oXHGhM1kpIYtM37ZdFbFzCCnVgS9Np88ef79nDyRviwZoXniynta61CgfEo/UCkVpJS6qZQ68tj2Vkqpk0qp00qpQcmfcVZr7ZPSz8gstNaM+mMUK0+vpFflXvi+7Puf+y9c+Asd7wVxpZAX+b0+z6CU9kdqOG39evZXhu0YRr1C9RjXeBzOz1gH6cDZq+wc/y4+t3/kTp5qZP90F1mkmUkRqeG0NXllEMvvTcFDK4JazaJA/peNjpQpWdLQFFBK+SilaimlPB5773keSDmLpMcp/E0p5Qj8BLQGKgDvK6UqPMe5Mw2tNd/t+Y7Ffy3G+yVvelXu9Z/7h+7YTZu/BnMrSwkKdwmSmyVTZxZSw2ki9HwoQ7YNoUaBGkxoMgEXx6c/4Vprzcr1W8kyuyWtzZu4UbUfBXsHo2RtmecxC6nhNPFz8AKW3BqHq4aAFv4ULpzhCzOLZJZccvqapNlMnYGXlFL3gSPJ/6R4FTat9ZbkG40fVRM4rbU+C6CUWgi8DhxL6fkzA601E/ZPYO7xuXxQ/gP6vdLvP2eCnLh0k0IhPXBzMOPsvRhcH+9LRUpIDaeNjRc3MnDLQCrnrczkppNxc3r6zemxCSbmzfmZDhdH4uDoSMxbi8hf0SsD09oXqeG0MSt0BfOvfYNS8HOjyZQo9uTV2EXGeOYIjdbaX2v9sda6kdbak6QnbP8E3AE2p1GOwsClR15fBgorpTyVUtOBqkqpL550oFKqu1Jqr1Jqb3h45ljbb/qh6QQdCeLdcu8yoMaA/2xmIh/Ec2Z2T15SZ4lvNw3nfC9kYNJM5blrGDJfHW+/sp3+m/tT3rM8U5tNJatz1qfue+VONMvH98bn0hfEeBQja59tuEszkx6khlNgwaZgfrn4JWYF0+v/QNnSTY2OlOk9z8J6l0kq9N/TPs6/PusW8J8Pb9Fa+wP+kPRAtPTOZLSHa3S8UeYNBtca/J/NjNmsWRU4is6J67laqTeFqr6RgUkFWFbDyftlmjrefW03n2z8hNI5SzOt+TQ8XJ4+Yrj3xHmiFvrQkb1cLfEmhTpNBecsGZhWSA3/2/LtGwg6PYAYB5he6xvKv9Da6EgCy+6hyQhXgEcn6xdJ3mYRpVRbpZR/ZGRkmgezJr8c+4WJ+yfyaslX+brO189co2PJyhW8FzGZy571KPTGyAxKmWmlqoYhc9TxgZsH6LOhD0WzFcW/hT85XHM8dd+V6zeTc0FrGrCf8AbfUKhLkDQz6Utq2AK/7tnJz8f6ct9RM7nal1Sq2N7oSCKZtTQ0e4CySqmSSikX4D1gtaUHa63XaK2758jx9B+Otm7xycV8t+c7mhdrzqj6o565Rse2g8doeKg/913yUdhnLjxjf5FqqaphsP86/jP8T3qt60X+rPkJaBlALrdcT9wvwWRm5i9BNNnyHvkdo4l9fzl5m30sN7KnP6nhZ1h3cD9TD3bnthOMf7k/1Sq/b3Qk8YgMb2iUUguAnUA5pdRlpZRP8to2fYAQ4DiwWGt9NAXntOtvBStOrWDkHyNpVKQR3zX87plrdFwMv4frSh9yqWjcOy9CZc2dQUkzh/So4eTz2m0dn7h9gh7repDLNRczWs4gT5YnP9fmTlQc8yYO5sPT/YnLWpCsfbbgXq5JBqe1f1LDKbf12FEm7O7CTSfN9y9+RN3q3kZHEo9RWtvPpc7q1avrvXv3Gh0jTQWfDWbQ1kHULlibyc0m4+ro+p/7x8SbCB7nzVtxK4loMZk89Z6+0J74b0qpfVrrDJ+DaW91fPrOabqFdMPNyY3ZrWZTyKPQE/c7c/0OR2d0p11iKFcLNKVQt19kRl4qSQ2njV2nTjF6w9tccjUxplRXvBr+z+hImYqldWwtl5xSxV6/Fay7sI7B2wZTLX81Jjad+MxmRmvNktkTeStuJZdf+FCaGRtjj3V8PvI8vqG+ODs4E9gy8KnNzK5jZ4iY/hrtEkO5Vqk3hbovk2bGBtljDR+6cIFvN3TgoquJ4UXflWbGitlFQ2OP1223XN7C51s+56U8LzGl2RSyOD37Zsg1Yet4+/K3XM1emSLvyHNDbY291fGl+5fwCfVBo5nRcgbFsj/5uWG/btxG3kVteIUT3GoxkYJvjgYHu/jRlOnYWw0fv3KNkb+356xrIl8WaEvbZl8ZHUn8hxRP2xbpb+fVnXy68VPK5izL1OZTcXd2f+YxB/86T6XtfYh3cqeAzyJwevqKq0Kkt+vR1/EL9SPOFEdgy0BK5fz3Q/rMZs38JQtoc+xznBwdiH9/BZ5lGxqQVoh/O3sjgmFr2vKXWzwD8zTnrVZjjI4knsEuvgbZ0zDn3ut76buhL8VzFMe/hT/ZXbI/85ib9x5wb4EPRVQ4Tu/NwSFHwQxIKtKavdRx+INwfEJ8iIyL5OfmP1Mud7l/7RObYGLOz9/R4djHJLrlJkuvTbhLM2Pz7KWGr9y5x+DlbTieJY5+2evR6bUJRkcSFrCLhsZehjkPhx+m9/reFPQoiH8Lf3K65XzmMQkmMxv8B9JQ7yWi3jA8yjbIgKQiPdhDHd+OvY1fqB/hMeFMaz6Ninkq/nufqDiWT+hH1xujuZWrCnk+2YJT3tIGpBVpzR5q+Oa9aAYsbM3RrA/4KEtVvN/82ehIwkJyyclKHL91nJ7repLbLTcBLQKeOq31cUsWzOS9+79wqWhbijbvm84phXi6yLhIeoT14HLUZaY1n0aVfFX+tc/5G3c4HuBNx8QNXC72OkU6z5DLo8Jq3ImOpf+8Vzmc9R4+zi/S6505RkcSKWAXIzS2Psx5+s5puod1x93ZnUCvQPK757fouHU7dvHqqa+4mbUMRTv7y8JjNs6W6zgqPope63px5u4ZJjaZSI0CNf61z6HTF7k+rR2tEzdwtUo/inSbLc2MnbHpGo5N4NM5r3Ew620+cCxJv/cXGx1JpJBdNDS2PMxp6bTWx528dINCId1xclB4ei8Cl6c/3E/YBlut4wcJD+i9vjfHbx1nXKNx1Ctc71/7bN1/GJdfXqM6RwlvNoFCbwyXBtwO2WoNx8Qn8sms19mX9QYdKMSAjiulPm2QXHIy0OX7l/EN9X3mtNbH3YuJ5/zs7rRQF7j3xjw85P4DYZDYxFj6buzLwfCDjG04libF/r2q72/rN1Jliy+5HKKJfms+eV9qZUBSIZ4sPtFMv5lvszvLJdqZ8vBVl99QsmyATZL/aga5Hn0d31BfYhJj8G/h/8RprU9iNmvWzBiBV+ImrlXpR87KbdI5qRBPlmBK4LNNn7H72m5G1htJqxL/bFS01ixauoj6Wzri7mRGdQ0mhzQzwookmsz0m9mJHW5naJWYg286h6Ic5Xu+rbKL/3JKqbZA2zJlyhgdxSIRMRH4hfoRGRfJjJYznjit9WmWr1rGOxFTuZi3IcXaDU1VjoSEBC5fvkxsbGyqzmPr3NzcKFKkCM7OzobmsKU6TjQnMmDLALZe2crQOkNpV7rdP983mVkwawrvXBzBXbeC5O6+BmfPEsaEFRnGlmrYbNZ8NtuXrS5HaBqflbFd16GcjP0ZIFJHnuWUwe7E3sE7xJsrUVf4ucXPVM1X1eJjdxw6RunlrXFwyUqez3agsjz5acWWOnfuHNmyZcPT0xOVSa8Xa625desW9+/fp2TJkv94T56D82Qms4nB2wYTfC6YATUG8GGFD//xfmyCiaU/D6dj+CSuZX+ZQr1WyQNSDSI1/GRaaz6f04cQtlA/zoXJXTbh5JrN6FjiKTLVs5xsxb34e/QI68HFexeZ3HRyipqZS+GRuK7wJod6gEeXhaluZgBiY2MzdTMDoJTC09Mz049SWcqszYz4YwTB54L55JVP/tXMRD6IZ82Ej/kgYiJX8jWkcN9QaWaE1RkyfyAhbKFWnCMTO4VKM2MnpKHJINEJ0fRa14tTd08xockEahWsZfGxsQkm9gf2oRrHud9iPFmKVE6zXJm5mXlI/g4so7Xm293fsvzUcnpU6oHvy77/eP/m3Wi2TfiQDtHzuFj8TYr2XA7Oz34GmRAZafjiYfyWEMwrcYop7/+Oi7un0ZFEGrGLe2isXUxiDH3W9+FoxFHGNRpHgyKWr+artWbZ7Al0il3NxbJdKFbvw2cfJEQa01rz474fWXBiAV0qdKF3ld7/eP/89VucD+hIG9MfXKzYk2JvfyvTXoXVGbvye1Y8WMZL8TD17TW4ZZPHxNgTuxihsebFnOJN8fTb2I99N/Yxuv5omhVvlqLjg9evp/2lsVzOXoVi741Lp5TCGlhzHU87NI2ZR2fybrl36V+9/z9GtY6du0TE9NdobPqDK7WGUazDWGlmMilrruGJwT+x6O5sXojX/PTGUtxzFzc6kkhjdtHQWOtiTgnmBPpv7s+OqzsYXnc4r5Z6NUXH/3n6Ai9t/Yg4J3cK+i4CR/u6A3/Pnj1UqlSJ2NhYoqOjqVixIkeOHDE6lmGstY6DjgQx7dA03ijzBoNrDf5HM7Pv6AnUrNeowklutJhC4dafGZhUGM1aazhg/Wzm3JhGiQQzP7WZR858LxodSaQDueSUTkxmE4O3DmbTpU0MrjWY9mXbp+j4W/djiJzvTXkVwYN3VuGYvUA6JTVOjRo1aNeuHV9++SUxMTF88MEHvPTSS0bHEo+Yd3weP+77kdYlWvN1na9xUP//HWjHvoMUXP0uBdQdItvPJX/llDXsQmSEuVuX4n/xewolmpnSfAZ5C//7GWPCPkhDkw7M2sywHcNYe34t/av15/0X30/R8YkmMxsCBtLBvJerdYZTqFz6P0F7+JqjHLt6L03PWaFQdoa1/ffTlh81dOhQatSogZubG5MmTUrTzxeps/zUcr7d/S1NizZlVINRODo4/v3ehm3beTHsQ7I7xBL/3jI8M6BGhUip5XvWMuXU1+Qxm5nccBKFStY1OpJIR9LQpDGtNaN3jWbVmVV8VPkjur7UNcXnWLZ4Nh0i53ChSBuKt/wk7UNakVu3bhEVFUVCQgKxsbG4u7sbHUkAv539ja93fE29wvX4vtH3ODv8/+XOtWEhVN/mi5OjA3T5lRzFXzEwqRBPtvbwZn7483OyaTMTao6mRLnmRkcS6UwamjSktWbc3nEsOrmIbi91o2flnik+x+Y/duN1Ygg3spSieJcZGXZz5bNGUtJLjx49GDlyJOfOnWPgwIFMmTLFkBzi/62/sJ4h24ZQvUB1fmz8Iy6O//9E7DVrltNob28SnNxx8/0Nt4KWr3ItREbZ/NdevtnTB1dMjKs0mHKV3jA6ksgA0tCkoamHpjL72Gzef/F9Pn3l0xSvb3LmSjj5fvfDyUGR23ux3T9Be86cOTg7O9OxY0dMJhN169Zlw4YNNG3a1OhomdbWy1v535b/UTFPRSY3nUwWp6R1ZLTWrFg8m9bHPueeS35y9fwNF0+ZJSKsz97zR/lqqzcOysR3ZT6mUo0PjI4kMohdNDTW8PyQoCNBTD80nfZl2jOo5qAUNzPRsQmcmeVHc3WBO+1+wSOf9T8LJbU6d+5M586dAXB0dGTXrl0GJzKW0XW85/oePt30KWVzlmVa82m4Oydd/jObNUt/mcIbZ4cRnqUUBXoH45gtnyEZhXUzuoaPXj3L5+s6kehg4ruiXanRsJchOYQxZNp2Gph/fH7STJCSrRlWZ9g/ZoJYQmvNmsCRtEzYyOVKffGs2jadkgprZmQdH7x5kN7re1M0W1F+bvEz2V2yA2Aya5bMGMNbZ7/iRraKFPokTJoZ8VRG1vCZ8Cv0DX6bGMcERuZtT/0WAzI8gzCWXYzQGGnFqRWM2T0maSZI/X/OBLHUmt9W8ubNKVzIU5/ib3yd9iGF+A/Hbh2j17pe5MuaD/8W/uRyS3pOWILJzMqfh/HuzUmcz1Wb4h8tR7nITdvC+ly5G06vVa9zzzGOr7M3p1nbb4yOJAwgDU0qBJ8NZtiOYU+cCWKpvUdOUHPPp0Q656OYz1xwsItBM2EjTt05RY+wHmR3yc6MljPImzUvAPGJZtZMHUiH2/6cy9OEkj0XgZOrwWmF+LfwqLv4Lm1LhFMsX7nUos1bE42OJAwiDc1zWn9hPYO3DX7iTBBLXb99H7W0GzlVNKbOK1FZU/8EbSEsdT7yPH6hfrg4uDCj5QwKuCct3hgbn8jvU/vz1t1ZnM3filLd59rdKtXCPkTGRtF14Wtcd4pikKpA+/cDjY4kDCTDAc9h25VtT5wJkhLxiWZ2B/ShGse42+x73IvJ6pUi41yJuoJvqC8aTYBXAEWzFwUgJi6R0Mm9aX93FmcLt6NUj/nSzAir9CAhlg/nteWy010+SyzJux8ulGeIZXLS0KTQnut76Lex379mgqTUyrmTaBezknOlP6BAg65pG1KI/3Aj+ga+Ib7EJMbg38KfUjlKARAVm8CGSb60u7+QM8XfoZTPbHiOe8KESG/xifF8MPd1zjlF0DcuPx92WyGX64VcckqJQ+GH6L2+N0U8ivxjJkhKhW3aRJtzo7mU7WVKvv9jGqcU4uluxdzCL8yPO3F3CGgRQLncSQvj3YuJY8ekrrSJCeZMqQ8p/eFk+bYrrFKiOZHO89/llMNVPorOjk/3YHCUX2VCGhqLHb91nF5hvcibJS8BLQP+ngmS4vOcv0SZjT2Id8ya9ARtp5TfeyPE84iMi6RHWA+uRV1jeovpvJz3ZQDuRsWwd/IHtIpbx9lyfpR+73tpZoRVMmszPou7clSfxve+Cz39QuRmdfE3qx6jU0q5K6VmK6UClFKdjMpx5u4ZeoT1wMPF4x8zQVIqMjqOiF+8KarCUR1m45SzcBontS1Dhw5lwoQJf78eMmQIEyfa1wwFa6nh6IRoeq3rxdnIs0xsOpFq+asBcCsyigMT36V53DrOvtSXUtLMiMdYSw1rrem9ojf74w7RORL6dA1FuXoYFUdYoQxvaJRSQUqpm0qpI49tb6WUOqmUOq2UGpS8+U1gqdbaD2iX0VkBLt67iF+oH04OTsxoOYOCHgWf6zxmsyYsYBANTLu5VnMIOcs3SuOktsfb25s5c+YAYDabWbhwIR98YP3LlNtaDcckxtB7fW+O3zrOuEbjqFso6YnDN+/c4/jkt2iSsJmzVQZQ6u2R0sxkErZWw1prBgQPYtv9bbwfmUCf99fg6OFpRBRhxYy45DQLmALMebhBKeUI/AS0AC4De5RSq4EiwJ/Ju5kyNiZci7qGb6gvCeYEZrWaRbHsxZ77XKuWzaX9nZmcK9iKkq0/S8OUaeT3QXD9z2fvlxIFXobW3z717RIlSuDp6cmBAwe4ceMGVatWxdPTJn5IzcJGajjeFM+nGz9l/439fNfwO5oUawLAtVt3uDD1Leqb9nGu+lBKvdY/o6MJY83CRmoYYMTGMayNCOate7H0eX0ZWfKWMCKGsHIZ3tBorbcopUo8trkmcFprfRZAKbUQeJ2k/6mKAAfJ4NGk8Afh+Ib6EhUfRaBXIKVzln7uc+3Yf4BGRwZx060EJboFyrfgR/j6+jJr1iyuX7+Ot7e30XEsYis1nGBO4PPNn7P96nZG1B1Bq5KtALgWfovL09pT03SYC/VGUbJln4yMJayArdQwwPgdU1h6aQFt7z+gT7NAshevlNERhI2wlpuCCwOXHnl9GagFTAKmKKXaAGuedKBSqjvQHaBYsecfQXnUndg7dA/rTnhMOP4t/CnvWf65z3Xxxm1yrPbGVZnJ2m2x9V7z/Y+RlPTUvn17hg4dSkJCAvPnzzckQxp57hqGtK9jk9nEl9u+ZMOlDXxR8wval20PwM07kVya/hbVTYe53OgHijf1TfVnCbthVTUMEHBgFjNP/UzLqAf0rjGWPBUap8l5hX2ylobmibTW0UC3Z+zjr5S6BrR1cXGpltrPvB9/nx5hPbh0/xJTm02lSr7nX/AuNsHE8cAeeHGW8DazyFvghdTGszsuLi40adKEnDlz4uhof2ueWFLDyfulWR1rrRn5x0iCzwXT75V+dCzfEYDb96I5/VMH6poOcL7+WEpIMyMsYEQNAyw8voxJh8fROPoBvcv9j8K13k7tKYWds5ZZTleAoo+8LpK8zSJp9YTXBwkP+GjdR5y6e4rxjcdTs2DN5z6X1prVQaPxig/lfIVe5K3RPlXZ7JXZbOaPP/7Ax8fH6CiplaoahrSpY601Y/eMZdmpZfSo1AOfl5P+XiOjYzky5V3qJu7ibI1hlGjR87k/Q9gtq6hhgDWngxm962vqxMTQu8CHlGreK1XnE5mDtTQ0e4CySqmSSikX4D1gtaUHK6XaKqX8IyMjnztAbGIsfTf05XDEYb5r+B0NizR87nMBrA35jdevTuB8ztqUeHtUqs5lr44dO0aZMmVo1qwZZcuWNTpOaqWqhiFt6njygcnMOz6PDyt8SO8qvQGIjo1n3+RONHbHjw0AAA6YSURBVIzfyukqAyjVxgpvShfWwCpqeOPFTXy5bRBV42Lpm6UlL7Yf9tznEpmLEdO2FwA7gXJKqctKKR+tdSLQBwgBjgOLtdZHLT1nar8VJJgS6L+5P7uv7+abet/QoniL5zrPQ4dPnqbyzr7cd/KkmO98WT7+KSpUqMDZs2cZN26c0VFSJD1qGFJfxwGHAwj4M4AOL3Tg8+qfo5QiNj6R7ZO8aRq7jlMV+lDmjSHPdW5hX6y1hnde2clnGz6hfHws/c2VeenDKTKJQljMiFlO7z9lezAQnMFxSDQnMnDrQLZc3sLQOkNpW7ptqs4XcS+auIXd8FT3iP+/9u49OKr67uP4+7chcQklyCWhSLgESB4khoIFCghIoaTyKAK1ts3DPNzC1cIUB31qKzbgIDzMoB1vQKEqVAKFCo9yMQihFQZa5SKaFoncSiQUDEm4yDWGnOcP0GaYXDa72T17zn5eM5nJnj1z8v3mfDb55mTPOaNy8OhaCa4TbhkGWPnpSl468BLDOgxjVu9ZGGO4/lU5f3llMkOvbORwp0xSHp1rR2kShsIxwweKDvDz3MdIKrvGry63pev0Vbo/k9SJK9Li72HOCquC3+z+DdsKtvFkjyd5NOXRgOoov1HB7qUz6GnlcXbAfBp36BnQ9iSy+JvjdYfXsWDvAoa0G8Kz9z2Lx3gov1HBtkWPM/TiWxxul0HKqOf1l64Enb8ZPlh8kElbJtKq7BpZ55uQNvX/dFsYqTNXDDT+HOa0LIvnPniOjcc3Mq3bNEanjg64jndWL2H4pbUcb/cTEgdNDHh7Eln8yfHm45uZ87c59GvdjwX9F9DA04AbFRY5S/6Hh879gc/uGknKmEUaZiQk/MnwkXNHyMzJpGnZFeYXe+gyZRMeb+MgViluFdanbQeLZVk8v+951h5ey/h7xjOp66SAt7lj927SjzzLyUZ30+G/X6mHKkVqtr1gO0/vepqe3+7Jbwf+luioaCzLYvOyLB4+u4zDCQ/wHxNe02F7CVsnLpxg3LvjaFj2JS8UXaPjhFyi4xLsLkscyhU/6ep6mHPxJ4tZ8ekKMjpnMOPeGZgA/3o9VniaxK0TqfDE0DJzre7+apP27dtTXFwc8Dp2qUuOd53axRM7n+CeFvfw8qCX8TbwYlkWG9/4Xx4+/SKHmw0kZXK23pAuIVWXDJ+6dIqx746D6+d5+cwF2maso2FChxBUKW7lioGmLoc53/jHGyz+ZDEjOo3gqV5PBTzMXLr2FYXLx9PenKb8R68R07x+rpApkcfXHO89s5cZf5lB8p3JLPrBImKjYwHYnP0iDxUs4EhcH5KnroGoiDwAKzbyNcNFV4oY++54rl0tZvGZYlqOeJO4pO4hqlLcyhUDja9W56/mhf0vMLT9UGb3mY3HBNa+ZVm8t+xp7i//KyfvfYLmaYGd7h2JTpw4QefOnRk7diwpKSmMGjWK3Nxc7rvvPpKTk9mzZw+lpaWMGDGCrl270rt3b/Ly8gAoKSkhPT2d1NRUJkyYgGVZ32x35cqV9OrVi27dujF58mRu3LDlnnr17pOznzBt+zQSv5XI74b8jriYOABy1izhgSNz+Gfj7nSath4T7bW5UpGqlV4rZfyWTM5dPs2SM1/Q7PsvE5/6fbvLEhdwxUDjy2HOt4++zbwP5zGwzUCe6/8cUfVwKH7ThjWMKF7K8fjBtB/264C3F6mOHj3KzJkzyc/PJz8/n1WrVrFr1y4WLlzIvHnzyMrKonv37uTl5TFv3jxGj775Bu45c+bQr18/Dh48yMiRI/n8888BOHToEGvWrGH37t18/PHHREVFkZ2dbWeLPqktx/ml+UzNnUrzhs1Zlr6Mpt6mAGx9ewWDP/01J2O7kDRtAyYmNpRli3yjtgxfuH6Bie9N4l8XCljyxRmafHc2rfsEdnapyNdccUzasqyNwMYePXpUeWrRln9uIeuvWfRp1YeF9y8k2hMd8Nfcl/d3+nz0JEUxbUjKXO74s0gW7FlAfml+vW6zc7PO/LLXL2tdLykpibS0NABSU1MZPHgwxhjS0tI4ceIEBQUFrFu3DoBBgwZRUlLCxYsX2blzJ+vXrwfgwQcfpGnTm7/gt2/fzv79++nZ8+Zp81evXiUhIfzfaFhTjo+fP86krZNoFN2I36f/nvjYeAByN6/h/gMzOe3tSJtpOjtE7FVThi9/dZkpuVM5du4wi4q+oEmnqST98DEbqhS3csVAU5sj54/QLb4bLw56kTuiAn/D7umS89yxfhyxpgxr7B8x3rh6qDJy3XHHv/eJx+P55rHH46G8vJzo6LoNoJZlMWbMGObPn1+vddqp8FIhDRs0ZGn6Uu761l0A7Nj2Dn33TKcoJpFWP8+hQaOmNlcpUr2SK6WcKjrKC2eLaJEwkuRHZttdkrhMRAw007tPp+xGGTFRgV+o6Xr5DT5eNpWhHOH0D5fQqnVqPVRoP1+OpNilf//+ZGdn88wzz/D+++/TokUL4uLiGDBgAKtWrWLWrFnk5ORw7tw5AAYPHszw4cN5/PHHSUhIoLS0lC+//JJ27drZ3In/BiQOoPfI3t9kePeOrdy7azLno+OJfyyHmLgWNlcoUrMbhf8it+AzCuKHkKxrI0kQuGKgMcYMA4Z16tSp2nXqY5gB2PSH53nk2rscS8mkY58qrx4u9Wz27NmMHz+erl27Ehsby4oVKwDIysoiIyOD1NRU+vbtS9u2N88w69KlC3PnziU9PZ2Kigqio6N59dVXw36gqS3HX2f4ww92kPrnsVxu0IQ7p+TgbdoqhFWKVK+mDCelfo+iy6/RqftDujaSBIWpfGaI0/Xo0cPat29f0Laf++dt9NuRwZm4rrSfsdXxp8UeOnSIu+++2+4ywkJV3wtjzH7LsnqEupaacvzR/g9pt+HHVETF4J20lcbf7hji6sRJwjHDInXla441Jvvo0LECOu+YyuWoJiROXO34YUacJy/vAK03/BTj8RAzbpOGGRGRSjTQ+OD8patcyB5DS3OOqIw3aRDX0u6SJMIcyj9I83WP4jXlMPodmrTRkTURkco00NSiosJix7KZ9K44wOm+s7kzua/dJUkEijq1n8aea5T91zqaJXWzuxwRkbDjioHG31vW+2LDn15n+IVsjrYeTtsh0+p9+3Zz03uo/BUu34OacpwyeDSNnvw78Snfs6EyEd8E82exSG1cMdD4c8t6X/xt7x4GfTqLQm8yHccscd1phl6vl5KSkrD5hW4Hy7IoKSnB67X/VgG15TgqVteZkfAWrJ/FIr7QO1urcfJMMc03Z2I8HuIz17rycvKJiYkUFhZy9uxZu0uxldfrJTEx0e4yREQkABpoqnD1ejlHXs9kICcpHpZN43h33tI+OjqapKQku8sQEREJmCv+5VSfLMsi5/VnGVT2PifSfkHCvQ/aXZKIiIjUQgPNbd7b8g7DzrzC8Wb96TAyy+5yRERExAcaaCrJy/+M7h/8gtLolrSf8KYuzy0iIuIQrngPzdf3DwEuGmOKgK/PGWxS6fMWQLHPG32mWU3PVt6uL2pbv7rnq1p++7KaHvvff83q0r+/vVf3XCj7D+nNnyrl+Iox5tCtxcHah7dvO9B13ZxhX9YP1/7DIcMQHvsxmPvw9mXV7VO9hv3r37ccW5blqg9gaTWf7wvG16iP9at7vqrlty+r6XE49O9v7+Hcf7A/QtVDKPZjuO7DSH8NB/ujlp4c9bPIl31Y036r/Fiv4eBm2I3/U9lYzefB+hr1sX51z1e1/PZlNT0Oh/797b2658Kh/2ALVQ+h2I/hug8j/TUcbLX1GKyvE8i6gezD25dFWv9hkWFX3W27JsaYfZYNd50NF+rf+f27oYdAqH939O+WPvwRyb1D8Pt34xGa6iy1uwCbqX/nc0MPgVD/7uCWPvwRyb1DkPuPmCM0IiIi4l6RdIRGREREXEoDjYiIiDieBhoRERFxvIgdaIwxjYwxK4wxy4wxo+yuJ9SMMR2MMa8ZY96yuxY7GGNG3Nr3a4wx6XbX4w9lWBlWhp1NGa7fDLtqoDHGvG6MKTLG/OO25Q8YYz4zxhw1xjx1a/GPgLcsy5oIPBzyYoOgLv1blnXcsqxMeyoNjjr2//atfT8F+Kkd9VZFGVaGlWFnU4bty7CrBhpgOfBA5QXGmCjgVWAo0AXIMMZ0ARKBk7dWuxHCGoNpOb7370bLqXv/s249Hy6Wowwrw5Uow46zHGXYlgy7aqCxLGsnUHrb4l7A0VuTcBnwR2A4UMjNFxO45PtQx/5dpy79m5sWADmWZX0U6lqrowwrwyjDjqYM25dhVwSoFq35918AcPMF1BpYDzxijFmMcy8v7osq+zfGNDfGLAG6G2N+ZU9pIVHd/p8O/AD4sTFmih2F1YEyrAwrw86mDIcgw66427Y/LMu6DIyzuw67WJZVws3/W0Yky7JeAl6yu45AKMPKMMqwoynD9ZvhSDhCcwpoU+lx4q1lkUL9O79/N/QQCPXv/P7d0EMg1H8I+o+EgWYvkGyMSTLGxAA/AzbYXFMoqX/n9++GHgKh/p3fvxt6CIT6D0X/lmW55gNYDZwGvuLm/+gyby3/T+AwcAx42u461b/6d3MP6j+y+3dDD+rfmf3r5pQiIiLieJHwLycRERFxOQ00IiIi4ngaaERERMTxNNCIiIiI42mgEREREcfTQCMiIiKOp4FGREREHE8DjYiIiDieBpoIZoyJsrsGkUAow+J0ynD9idi7bUcqY8yfgFLgO8AmYK69FYnUjTIsTqcMB4cGmsiTBqy1LKu33YWI+EkZFqdThoNA93KKIMYYL/A5cJdlWeV21yNSV8qwOJ0yHDx6D01kSQU+1ItIHEwZFqdThoNEA01kSQPy7C5CJADKsDidMhwkGmgii15I4nTKsDidMhwkeg+NiIiIOJ6O0IiIiIjjaaARERERx9NAIyIiIo6ngUZEREQcTwONiIiIOJ4GGhEREXE8DTQiIiLieBpoRERExPH+HxemIkB9dBAhAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def plot_all(r, results, meta):\n", " plt.figure(figsize=(9, 3))\n", " for i, title in enumerate([\"Woofer\", \"Tweeter\", \"Sum\"]):\n", " plt.subplot(1, 3, i + 1)\n", " sf_x, sf_y = results[i]\n", " plt.loglog(r, sf_x, label=\"x\")\n", " plt.loglog(r, sf_y, label=\"y\")\n", " plot_theory(r, meta)\n", " if i == 0:\n", " plt.ylabel(r\"$D_\\phi(r)$\")\n", " plt.legend(loc=\"lower right\")\n", " plt.xlabel(r\"$r$\")\n", " plt.title(title)\n", " plt.ylim(2e-1, 2e3)\n", "\n", "\n", "def plot_theory(r, meta):\n", " model = sf_integrated(r, r0=meta[\"r0\"], L0=meta[\"L0\"])\n", " plt.loglog(r, model, label=\"model\")\n", "\n", "\n", "plot_all(r, sf, args)\n", "# plt.savefig(\"component_sf.png\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "jupytext": { "formats": "ipynb,py:percent" }, "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": 2 }