diff --git a/docs/source/index.rst b/docs/source/index.rst index 688308baa..5dfa8a232 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -72,6 +72,7 @@ new ideas. compressible-rt-compare.ipynb adding_a_problem_jupyter.ipynb + rt_average.ipynb advection-error.ipynb compressible-convergence.ipynb diff --git a/docs/source/rt_average.ipynb b/docs/source/rt_average.ipynb new file mode 100644 index 000000000..14fdaa201 --- /dev/null +++ b/docs/source/rt_average.ipynb @@ -0,0 +1,155 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e405b84e-09c0-4c8f-b178-02102b2782be", + "metadata": {}, + "source": [ + "# Horizontal Averages of Rayleigh-Taylor" + ] + }, + { + "cell_type": "markdown", + "id": "6395d3eb-9572-454b-823a-b5c2f8c4267d", + "metadata": {}, + "source": [ + "Here we show how to do a lateral / horizonal average of the compressible Rayleigh-Taylor\n", + "simulation to see what the average vertical profile of the mixing looks like." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "571afb49-3a5d-49d1-a7e3-8e51bb6c7d11", + "metadata": {}, + "outputs": [], + "source": [ + "from pyro import Pyro" + ] + }, + { + "cell_type": "markdown", + "id": "97989792-3858-43f4-8e59-dca0a2686777", + "metadata": {}, + "source": [ + "First we'll run the `rt` problem using the problem defaults." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7e5a5084-c403-41d4-aec7-0bd48f0be310", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/zingale/development/pyro2/pyro/compressible/problems/rt.py:71: RuntimeWarning: invalid value encountered in divide\n", + " 0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :]\n" + ] + } + ], + "source": [ + "p = Pyro(\"compressible\")\n", + "p.initialize_problem(\"rt\")\n", + "p.run_sim()" + ] + }, + { + "cell_type": "markdown", + "id": "9ab0b191-8a7a-409f-8773-ba0f1b1dd49e", + "metadata": {}, + "source": [ + "Now we'll get the density variable and use `np.average()` to compute the average in the x-direction.\n", + "Note that we operate only on the valid data (`dens.v()`)." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "bdbb2ac1-54ae-4f5a-b96a-fdc3eb602221", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "1a0909a9-5a61-4606-9e87-83d4559213c6", + "metadata": {}, + "outputs": [], + "source": [ + "dens = p.get_var(\"density\")\n", + "dens_avg = np.average(dens.v(), axis=0)" + ] + }, + { + "cell_type": "markdown", + "id": "81ae06fb-6125-4a67-a4f1-fbd157f79b03", + "metadata": {}, + "source": [ + "Finally, we can plot the profile." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "950658f0-5c24-4392-94a9-17dbfc769a63", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'average density')" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMe0lEQVR4nO3deXwU9f0/8NfsbnZzX4RcJJBAuI8QzoIHiBwCxaIVVGyB4lEVK15Y6E+l1CPiVyu21Vq0hWK9L7QeSAQJilAECVfkDiSE3NdmN8me8/tjDwgkkE12d3ZnXs/HIw+yk5ndd4awefE5BVEURRARERHJhErqAoiIiIi8ieGGiIiIZIXhhoiIiGSF4YaIiIhkheGGiIiIZIXhhoiIiGSF4YaIiIhkRSN1Af5mt9tx9uxZREVFQRAEqcshIiKiDhBFEY2NjUhNTYVKdem2GcWFm7NnzyI9PV3qMoiIiKgTSkpKkJaWdslzFBduoqKiAOfNiY6OlrocIiIi6gC9Xo/09HT37/FLUVy4cXVFRUdHM9wQEREFmY4MKeGAYiIiIpIVhhsiIiKSFYYbIiIikhWGGyIiIpIVhhsiIiKSFYYbIiIikhWGGyIiIpIVhhsiIiKSFYYbIiIikhWGGyIiIpIVhhsiIiKSFYYbIiIikhXFbZxJRESkFKIowmYXYbGJMNvssNjssNpEWO12qFWC40MQoFGpoFIBapXjc60muNs+GG6IiCjgiKIIk9UOk8UOs83uPi4IgOD+XIBKALQaFbRqFTTq4P6F7CKKIprMNhhMVjS2WGEwWWFoscJgsqCx5bxjrb5ucT92HWs221rdO09oVALCtWqEazUI16kRodUgTKtGhFaN6LAQRIeGICYsBNFhGkSHhrR5LC5C6/V70+H6JXtlIiKSPVEUUW0wo7S+GZX6FtQ1mVFjNKPOaEat0YJaowm1TRY0NlvQbLGhxWJDi8WOFqsNoujZa6kEQKdRO8KORgWd80+tWgVdiBo6tQq6kHPHdRr1eZ+r3NfqXB8haue1KoSoVRCcgQo4P2A5PhzHBFhsdlhsIiw2O8xWu7u1xGy1o9lig9FkhdHs/NNkQ5PZ2uqYweR4bPfwe/eERiVApRJgt4uwiWKb99lqF6FvsULfYu3Ua8SGh6DgialdL7aTGG6IiKhLRFFEWUMLTlQZcLzSgBNVBhRVG3G2vgWl9c0wWzvXeuApuwg0W2xottj88nq+plYJiNRpEKnTICr0vD9DQ1odi9RpEBmqQZTzT9fXwrQaaNUqhKgFhKgdAc0VbM7nCjk2uwi7KMJiFdFksaLJbEOTM4A1mW0wOoNYY4sV+maLI/w0W9DQbIG+xQJ9s9X5pwXRoSGS3Tcw3BARkacq9C3YW1yPgpJ67Cupx4HSBhhM7f8PXyUASdGhSIoORbcILeIitIiP0CIuXIv4iBDER+gQHero9ggNUSNUo0ZoiAqhWsfnIWrB3WICZ5hysdkdY0nMVseHyf1hcz9u/eeFx22tznE9bnXMYoPFZofofn3nn+c9EJ2fapxBQqdROQOFAK1GjRCVgDCtGpE6DcK1GkTo1IjQaRB+3rFInaMLKCpUgyhdCEJDVK2+b19RqQSoICBE7TygBWLQtXBi82XTUwcw3BAR0SWZrDZsPVKFrUeqsO1oFUrrmy86R6MS0KtbOLISI9GneyR6d49EWlwYesSGITkmFCFeHA9z/i98jVqARq1CuHTDO6gNapXvQ9mlSBpucnNz8dFHH+Hw4cMICwvD+PHjsWrVKvTv3/+S173//vt4/PHHcerUKfTt2xerVq3CjBkz/FY3EZESlNQ24a1dxXjvhxLUGM3u4yoB6JcUhZyesRieHovs9Fj06R7p1QBD1BWShpv8/HwsXrwYo0ePhtVqxR/+8AdMnToVhYWFiIiIaPOa77//Hrfeeityc3Px85//HG+99RZmz56NH3/8EUOGDPH790BEJCc2u4htR6vwn52nseVIpbsLJjk6FNcNScaEft0xJjMeETo2/FPgEkTR0/HovlNVVYXExETk5+fj6quvbvOcm2++GUajEZ999pn72M9+9jMMHz4cr7766kXnm0wmmEwm92O9Xo/09HQ0NDQgOjraR98JEVFwaWiy4O0fivHm/06jpPZct9NVfRNw29hemDwwUTZTrSk46fV6xMTEdOj3d0BF74aGBgBAfHx8u+fs2LEDDz30UKtj06ZNw4YNG9o8Pzc3FytXrvRypURE8tDQZMHr353Euu2n0OgcFBwdqsGcUem4bWxP9O4eKXWJRB4LmHBjt9vxwAMP4Iorrrhk91J5eTmSkpJaHUtKSkJ5eXmb5y9fvrxVGHK13BAREfCbdbvwY3E9AKB/UhRuvyoTs4alIkyrvuy1RIEqYMLN4sWLcfDgQXz33XdefV6dTgedTufV5yQikoPTNUb8WFwPjUrAX2/NwbTByRetgUIUjAIi3Nx333347LPPsG3bNqSlpV3y3OTkZFRUVLQ6VlFRgeTkZB9XSUQkL9uOVQMARvaKw/ShKVKXQ+Q1ko4OE0UR9913Hz7++GNs2bIFmZmZl71m3Lhx2Lx5c6tjeXl5GDdunA8rJSKSn21HqwAAV/frLnUpRF4lacvN4sWL8dZbb+GTTz5BVFSUe9xMTEwMwsLCAADz589Hjx49kJubCwBYsmQJJkyYgBdeeAEzZ87EO++8g927d2PNmjVSfitEREHFYrNjx4kaAMDVfRluSF4kbbn5+9//joaGBkycOBEpKSnuj3fffdd9TnFxMcrKytyPx48fj7feegtr1qxBdnY2PvjgA2zYsIFr3BAReWBvcT0MJiviI7QYnMplMUheJG256cgSO1u3br3o2Jw5czBnzhwfVUVEJH/fHnN0SV2ZlcBBxCQ7XJGJiEiBXONtruqbIHUpRF7HcENEpDB1RjP2lzoWTb2K421IhhhuiIgUZn9pA0QR6N09AskxoVKXQ+R1DDdERApT59zhOymKwYbkieGGiEhhGpotAICYsBCpSyHyCYYbIiKF0TPckMwx3BARKYy75Sac4YbkieGGiEhh2C1FcsdwQ0SkMAw3JHcMN0RECsNwQ3LHcENEpDAMNyR3DDdERArD2VIkdww3REQKw5YbkjuGGyIiBbHY7DCabQDDDckYww0RkYK4uqQAIJrhhmSK4YaISEFcXVJROg3UKkHqcoh8guGGiEhBXOGGrTYkZww3REQKwsHEpAQMN0RECsJwQ0rAcENEpCBc44aUgOGGiEhB2HJDSsBwQ0SkIO5wE85wQ/LFcENEpCBsuSElYLghIlIQTgUnJWC4ISJSELbckBIw3BARKUhDsxVguCGZY7ghIlIQTgUnJWC4ISJSEHZLkRIw3BARKYTVZofBxG4pkj+GGyIihdC3WN2fR4dqJK2FyJcYboiIFMLVJRWp00Cj5ts/yRd/uomIFILjbUgpGG6IiBSCC/iRUjDcEBEpxLmWG463IXljuCEiUgh2S5FSMNwQESkEF/AjpWC4ISJSCLbckFIw3BARKURDE8MNKQPDDRGRQuhbOFuKlIHhhohIIRqdKxRHcXVikjmGGyIihWh0ttxE6dhyQ/LGcENEpBBsuSGlYLghIlIIvTvcsOWG5I3hhohIIdzdUmy5IZljuCEiUgCz1Q6T1Q4AiGbLDckcww0RkQK4Wm0AIJItNyRzDDdERArgGkwcoVVDrRKkLofIpxhuiIgUoJGDiUlBGG6IiBSg0eTolmKXFCkBww0RkQJwjRtSEoYbIiIFYLcUKQnDDRGRAnCNG1IShhsiIgVwtdxEM9yQAjDcEBEpwLmWG3ZLkfwx3BARKYB7zI2OLTckf5KGm23btmHWrFlITU2FIAjYsGHDZa958803kZ2djfDwcKSkpGDRokWoqanxS71ERMGKs6VISSQNN0ajEdnZ2Xj55Zc7dP727dsxf/583H777Th06BDef/997Nq1C3feeafPayUiCmZ6dkuRgkga4adPn47p06d3+PwdO3YgIyMD999/PwAgMzMTv/3tb7Fq1ap2rzGZTDCZTO7Her2+i1UTEQUfttyQkgTVmJtx48ahpKQEX3zxBURRREVFBT744APMmDGj3Wtyc3MRExPj/khPT/drzUREgYADiklJgircXHHFFXjzzTdx8803Q6vVIjk5GTExMZfs1lq+fDkaGhrcHyUlJX6tmYgoELDlhpQkqMJNYWEhlixZgieeeAJ79uzBxo0bcerUKdx9993tXqPT6RAdHd3qg4hIac6tc8OWG5K/oIrwubm5uOKKK7B06VIAwLBhwxAREYGrrroKTz31FFJSUqQukYgo4FhsdjRbbABbbkghgqrlpqmpCSpV65LVajUAQBRFiaoiIgpsBmerDbgrOCmEpOHGYDCgoKAABQUFAICioiIUFBSguLgYcI6XmT9/vvv8WbNm4aOPPsLf//53nDx5Etu3b8f999+PMWPGIDU1VbLvg4gokLm6pMJC1AhRB9X/aYk6RdIIv3v3blxzzTXuxw899BAAYMGCBVi3bh3KysrcQQcAFi5ciMbGRvztb3/Dww8/jNjYWEyaNOmSU8GJiJROz00zSWEEUWH9OXq9HjExMWhoaODgYiJShB0nanDrazvRu3sEtjw8UepyiDrFk9/fbJ8kIpI5rnFDSsNwQ0Qkc+emgbNbipSB4YaISOYMJi7gR8rCcENEJHPubikdu6VIGRhuiIhkjlsvkNIw3BARyZzeHW7YckPKwHBDRCRzjVznhhSG4YaISObYLUVKw3BDRCRzXOeGlIbhhohI5rjODSkNww0Rkcw1ckAxKQzDDRGRzHFAMSkNww0RkYxZbHYYzTYAQHQYW25IGRhuiIhkrMZgBgCoVQJiGW5IIRhuiIhkrNpgAgDER2ihUglSl0PkFww3REQyVuUMNwmROqlLIfIbhhsiIhmrbnSEm+5RDDekHAw3REQyVu0cc5MQqZW6FCK/YbghIpIx15ib7uyWIgVhuCEikrFqjrkhBWK4ISKSMXe4iWK3FCkHww0RkYxVNbLlhpSH4YaISMbODShmuCHlYLghIpIpq82OuiaGG1IehhsiIpmqNZohioBKcKxQTKQUDDdERDJV5d56QQc1t14gBWG4ISKSKS7gR0rFcEMkExv2lmLkk3nYdKhc6lIoQHDrBVIqhhsiGThdY8QfPj6AGqMZT3xyCE1mq9QlUQDgAn6kVAw3REHOZhex9P39aDLbAADl+ha8tq1I6rIoAJwLN+yWImVhuCEKcmu3F2HXqVpEaNX4/XUDAACv5p9Ahb5F6tJIYlzAj5SK4YYoiFXqW/DCpqMAgMd+Pgh3T+iNET1j0Wyx4c/O46RcXMCPlIrhhiiIvfj1MTRbbBjRMxa3jE6HIAh4ZFp/AMCmQg4sVrpz+0ox3JCyeBxuJkyYgPXr16O5udk3FRFRhxyvNOC93SUAgOUzBkIQHOuYDE6JAQDUNVnQ7ByHQ8rkCjfd2XJDCuNxuMnJycEjjzyC5ORk3Hnnndi5c6dvKiOiS3pu42HY7CImD0zC6Ix49/HoMA3CtWoAQFkD/xOiVDa7iFqjs1uKO4KTwngcblavXo2zZ89i7dq1qKysxNVXX41Bgwbh+eefR0VFhW+qJKJWfiyuw6bCCqgE4PfX9W/1NUEQkBobBgAoa+CgYqWqNZphFwFBAOLDGW5IWTo15kaj0eDGG2/EJ598gjNnzmDevHl4/PHHkZ6ejtmzZ2PLli3er5SI3P66+RgA4Jcj0tA3Keqir6fEhAIASuvZcqNUri6p+HAtNGoOryRl6dJP/K5du7BixQq88MILSExMxPLly5GQkICf//zneOSRR7xXJRG5HSxtwDdHqqASgHuvyWrznNQYZ8tNPVtulIoL+JGSaTy9oLKyEm+88QbWrl2LY8eOYdasWXj77bcxbdo094DGhQsX4rrrrsPzzz/vi5qJFO1vW44DAGZlpyIzIaLNc1JiHS03HHOjXJV610wpdkmR8ngcbtLS0tCnTx8sWrQICxcuRPfu3S86Z9iwYRg9erS3aiQip6MVjdjo3DvqvnZabXBey81ZjrlRrJK6JgBAely41KUQ+Z3H4Wbz5s246qqrLnlOdHQ0vvnmm67URURteHXrCQDA9CHJbY61cXEPKOaYG8UqrnGGm3iGG1Iej8fcrFixAvX19Rcd1+v1mDRpkrfqIqILVDWa8N/9ZwEAv53Q55LnurqlztY3QxRFv9RHgeV0rSPc9OrGcEPK43G4yc/Ph9lsvuh4S0sLvv32W2/VRUQXeHtXMSw2EcPTYzE8PfaS57q6pYxmG/Qt3CFciYqd4aYnW25IgTrcLbV//34AgCiKKCwsRHn5uaXdbTYbNm7ciB49evimSiKFs9jsePN/pwEAC8dnXPb8MK0aseEhqG+yoKyhGTFhIX6okgJFk9nq3jSzV3zbg86J5KzD4Wb48OEQBAGCILTZ/RQWFoa//vWv3q6PiABsPFiOCr0J3aN0mDE0pUPXpMSEOcJNfQsGJEf7ukQKICW1jrFW0aEaxIQz2JLydDjcFBUVQRRF9O7dG7t27Wo1S0qr1SIxMRFqtdpXdRIp2r+/PwUAmDemJ7SajvUm94gNxU9lepzldHDFOV1jBAD06sZWG1KmDoebXr16AQDsdrsv6yGiCxyvNGD36TpoVAJuG9uzw9eluKaDc8aU4nC8DSldh8LNp59+iunTpyMkJASffvrpJc+9/vrrvVUbEQH4pKAUADChX3ckRod2+Dr3Qn5cpVhx3OGGM6VIoToUbmbPno3y8nIkJiZi9uzZ7Z4nCAJsNps36yNSNFEU8fFeR7iZnePZgP1zC/mx5UZp2HJDStehcHN+VxS7pYj8Z/fpOpypa0akToPJA5M8upY7gyuXawG/Xgw3pFBe2Sq2rUX9iKjrXK021w1JRpjWswH7rp3ByxpaYLdzIT+lsNlFnKlztNZxdWJSKo/DzapVq/Duu++6H8+ZMwfx8fHo0aMH9u3b5+36iBTLZLXh8/1lAIAbPeySAoDkmFAIAmC22lFjvHjhTZKncn0LzDY7NCrB3XpHpDQeh5tXX30V6enpAIC8vDx8/fXX2LhxI6ZPn46lS5f6okYiRco/UoWGZguSo0Mxtnc3j68PUavQPVIHcHdwRXF1SaXFhUGtEqQuh0gSHm+cWV5e7g43n332GebOnYupU6ciIyMDY8eO9UWNRIr01aEKAMCMoSmd/iWVGK1DZaMJNQa23ChFca1jjZueXOOGFMzjlpu4uDiUlJQAADZu3IjJkycDzlkdnClF5B1Wmx2bDzvCzbTBng0kPl9cuBYAUMtuKcU4N1OKXVKkXB6HmxtvvBHz5s3DlClTUFNTg+nTpwMA9u7di6ysLI+ea9u2bZg1axZSU1MhCAI2bNhw2WtMJhP+3//7f+jVqxd0Oh0yMjLwr3/9y9Nvgyig/XCqDvVNFsSFh2Bkr7hOP098BMON0px2z5Riyw0pl8fdUi+++CIyMjJQUlKC5557DpGRkQCAsrIy3HvvvR49l9FoRHZ2NhYtWoQbb7yxQ9fMnTsXFRUV+Oc//4msrCyUlZVxejrJTl6ho9Xm2oFJ0Kg7P6nRHW6aGG6U4kSVa+sFzpQi5fI43ISEhOCRRx656PiDDz7o8YtPnz7d3fLTERs3bkR+fj5OnjyJ+Ph4AEBGxqV3SDaZTDCZTO7Her3e4zqJ/EkURWwqLAcATB3U+S4pAIh3dkvVseVGESw2O05UGgCAm6WSonkcbgDg2LFj+Oabb1BZWXlRq8kTTzzhrdou8umnn2LUqFF47rnn8MYbbyAiIgLXX389nnzySYSFtd2/nJubi5UrV/qsJiJv+6msEWfqmhEaosJVfbt34Ir2xbFbSlFO1xhhttkRrlUjLY5jbki5PA43r732Gu655x4kJCQgOTkZgnBuFocgCD4NNydPnsR3332H0NBQfPzxx6iursa9996LmpoarF27ts1rli9fjoceesj9WK/Xu2d7EQUiV6vNVX27e7xw34Vc3VJ17JZShCPljlabvklRUHEaOCmYx+HmqaeewtNPP43f//73vqnoEux2OwRBwJtvvomYmBgAwJ///GfcdNNNeOWVV9psvdHpdNDpdH6vlaizvjlcCQCY0sUuKXC2lOIcKXd0u/dPipS6FCJJeTxSsa6uDnPmzPFNNZeRkpKCHj16uIMNAAwcOBCiKOLMmTOS1ETkTXVGM/aXNgDOXcC76lzLjaXLz0WB70hFIwCgP8fbkMJ5HG7mzJmDTZs2+aaay7jiiitw9uxZGAwG97GjR49CpVIhLS1NkpqIvOn7EzUQRaBfUiSSokO7/HxxESEAgPomM2zcX0r2jlY43hv7J0VJXQqRpDzulsrKysLjjz+OnTt3YujQoQgJCWn19fvvv7/Dz2UwGHD8+HH346KiIhQUFCA+Ph49e/bE8uXLUVpaivXr1wMA5s2bhyeffBK/+c1vsHLlSlRXV2Pp0qVYtGhRuwOKiYLJd8erAABXZnW91QbndUvZRUDfbHEPMCb5abHYcKrGMQ28fzLDDSmbx+FmzZo1iIyMRH5+PvLz81t9TRAEj8LN7t27cc0117gfuwb+LliwAOvWrUNZWRmKi4vdX4+MjEReXh5+97vfYdSoUejWrRvmzp2Lp556ytNvgyjgiKKIb49VAwCu6pvglecMUasQFapBY4sVtU1mhhsZO1ZhgCg6uiITIvn3TMrmcbgpKiry2otPnDgRoth+U/m6desuOjZgwADk5eV5rQaiQHG6pgln6poRohYwtne81563W4QWjS1Wx1o33mkQogDkGm/TLymy1SxWIiXq9NKnZrMZR44cgdVq9W5FRAr17XFHq82InnEI13ZqCao2uVprajhjStaOOsMNF+8j6kS4aWpqwu23347w8HAMHjzY3W30u9/9Ds8++6wvaiRShG+POsbbeKtLyoWrFCvD4XJXyw3H2xB5HG6WL1+Offv2YevWrQgNPTebY/LkyXj33Xe9XR+RIlhtduw4WQMAuLKLqxJfKI77SynC0XLXNHCucUPkcdv3hg0b8O677+JnP/tZq37dwYMH48SJE96uj0gRfiprRGOLFVGhGgztEdOBKzrOvdYNW25kq77JjHJ9C+BcnZhI6TxuuamqqkJiYuJFx41GIwexEXXSrlO1AIBRveKg9vKy+edWKeZCfnJ1wLnwY69u4YgODbns+URy53G4GTVqFD7//HP3Y1egef311zFu3DjvVkekELud4WZ0pvdmSbnEOxfy4/5S8rX/jCPcDEuLlboUooDgcbfUM888g+nTp6OwsBBWqxUvvfQSCgsL8f3331+07g0RXZ4oivjBFW4yvB9uuL+U/B1whRsvd2kSBSuPW26uvPJKFBQUwGq1YujQodi0aRMSExOxY8cOjBw50jdVEsnYqZomVBvM0GpUGJbm/V9O3Blc/vafqQcADPXBzw9RMOrUYhp9+vTBa6+95v1qiBTohyJHq012Wgx0GrXXn98VbthyI09VjSacbWiBIABD2HJDBHQ03Oj1+g4/YXQ0F5Ai8oQvu6RwXrhpbLHCYrMjRN3ptTspAB10Dibu0z0SkTrvLf5IFMw69C8hNja2wzOhbDZbV2siUhRfh5vo0BCoBMfmmXVNZiRGdX23cQoc+5xdUhxvQ3ROh8LNN9984/781KlTWLZsGRYuXOieHbVjxw78+9//Rm5uru8qJZKhysYWnKppgiAAI3rF+eQ1VCoBceFa1BjNqDUy3MiNazAxx9sQndOhcDNhwgT353/605/w5z//Gbfeeqv72PXXX4+hQ4dizZo1WLBggW8qJZKhPafqAAD9k6IQE+a79UniIs6FG5IPURSxv5TTwIku5HHn+44dOzBq1KiLjo8aNQq7du3yVl1EirC3xNGlMNJHrTYu5/aX4kJ+clKub0FVowlqlYBBKRzvSOTicbhJT09vc6bU66+/jvT0dG/VRaQIrim82em+/V93nHMhP+4vJS+uxfv6JkYiTOv9mXZEwcrjofUvvvgifvnLX+LLL7/E2LFjAQC7du3CsWPH8OGHH/qiRiJZsttFHCx1zET0xfo25+P+UvK0t9gZjtklRdSKxy03M2bMwLFjx3D99dejtrYWtbW1mDVrFo4ePYoZM2b4pkoiGTpZbYTBZEVoiApZ3X27kzPXupGnPacdM+1GZvi2W5Mo2HRqUYS0tDQ8/fTT3q+GSEEOlDr+1z0kNQYaH689wy0Y5MdktWGfs1vKV8sIEAUrruZFJJF9Jf6bwpsQqQMA1BhNPn8t8o+DpQ0wW+3oFqFFRrdwqcshCigMN0QSOeCewuv7cOPqlqoxsOVGLn5wLiMwKiOuw4usEikFww2RBKw2Ow6d9d/6JN0ineGG3VKysdsVbnqxS4roQgw3RBI4VmlAi8WOSJ0Gmd0ifP563SIc3VJ1RjPsdtHnr0e+JYqiezDxKA4mJrpIp8KN1WrF119/jX/84x9obGwEAJw9exYGg8Hb9RHJkmvJ/CE9oqFS+b5LwdUtZbWL0LdwIb9gd6LKiLomC3QaFQanctsFogt5PFvq9OnTuO6661BcXAyTyYQpU6YgKioKq1atgslkwquvvuqbSolkxLXZob/WJ9FqVIgK1aCxxYpqgxmxztlTFJx2OzdbHZ4eC62GDfBEF/L4X8WSJUswatQo1NXVISwszH38hhtuwObNm71dH5EsHSx1tdz473/drhlTnA4e/HafPjeYmIgu5nHLzbfffovvv/8eWm3r//llZGSgtLTUm7URyZLNLuJIhaM7d3Cq//YDio/QoqjaiBoDp4MHu11FrvE2HExM1BaPW27sdjtsNttFx8+cOYOoqChv1UUkWyW1TWix2KHTqNDLD4OJXbpFcMaUHJTUNqG4tglqlcDF+4ja4XG4mTp1KlavXu1+LAgCDAYDVqxYwe0XiDrgcLmj1aZvUiTUfhhM7OKeDs61boLajhM1AIDstBhE6jq1yDyR7Hn8L+OFF17AtGnTMGjQILS0tGDevHk4duwYEhIS8Pbbb/umSiIZOeIMN/2T/NclhfOmg9dyleKg9v2JagDAFVkJUpdCFLA8DjdpaWnYt28f3nnnHezfvx8GgwG33347brvttlYDjImobUed4236J/t2s8wLuaaDV7NbKmiJoojtzpabcX26SV0OUcDqVJumRqPBr371K+9XQ6QAh8v1AID+yX5uuXF2S9WyWyponagyoKrRBJ1GhRE9OVOKqD0eh5tPP/20zeOCICA0NBRZWVnIzMz0Rm1EstNiseFUTRMAYECyfwfgu7qluHlm8Np+3NFqMyojDqEhaqnLIQpYHoeb2bNnQxAEiGLrJdxdxwRBwJVXXokNGzYgLo7/syA634kqA2x2ETFhIUiM0vn1tTmgOPi5xtuM78PxNkSX4vFsqby8PIwePRp5eXloaGhAQ0MD8vLyMHbsWHz22WfYtm0bampq8Mgjj/imYqIg5h5MnBzl952cXeGmrskMG/eXCjo2u4idJx3r24zneBuiS/K45WbJkiVYs2YNxo8f7z527bXXIjQ0FHfddRcOHTqE1atXY9GiRd6ulSjoucKNv7ukACDOueWCXQTqm8zoFunfliPqmoOlDWhotiBSp8FQP65sTRSMPG65OXHiBKKjLx4IGR0djZMnTwIA+vbti+rqau9USCQjrjVu+iX5P9yEqFWIDQ8BuAVDUNp6pAoAcGVWAjRq7idFdCke/wsZOXIkli5diqqqKvexqqoqPProoxg9ejQA4NixY0hPT/dupUQy4JoGLkXLDc6fDs5xN0HnmyOVAICJ/btLXQpRwPM43Pzzn/9EUVER0tLSkJWVhaysLKSlpeHUqVN4/fXXAQAGgwGPPfaYL+olCloNTRaUNbQAAPpJFG4SIrh5ZjCqNZrdO8lP7J8odTlEAc/jMTf9+/dHYWEhNm3ahKNHj7qPTZkyBSqVIyvNnj3b+5USBbnjVQYAQHJ0KKJDQySpId69vxSngweTbUerIIqOFr/kmFCpyyEKeJ1axE+lUuG6667Ddddd5/2KiGTqVLURAJCZ4L/NMi/kmjHFbqngstXZJXXNALbaEHVEp8KN0WhEfn4+iouLYTa3fpO8//77vVUbkaycrnGEmwwpw42z5Yb7SwUPm11E/lHHGMeJ/TjehqgjPA43e/fuxYwZM9DU1ASj0Yj4+HhUV1cjPDwciYmJDDdE7Shyrkyc0S1cshpc07+5kF/w2HemHnVNFkSFajCiFxdGJeoIjwcUP/jgg5g1axbq6uoQFhaGnTt34vTp0xg5ciSef/5531RJJAOubikpW27OjblhuAkWWw87uqSu6puAEE4BJ+oQj/+lFBQU4OGHH4ZKpYJarYbJZEJ6ejqee+45/OEPf/BNlURBThRFnKoJnDE3NQZ2SwWLTYUVAIBrByRJXQpR0PA43ISEhLhnRSUmJqK4uBgAEBMTg5KSEu9XSCQDtUYzGlusAICe8RJ2S7k3z2TLTTA4XWPE4fJGqFUCrh3IwcREHeXxmJucnBz88MMP6Nu3LyZMmIAnnngC1dXVeOONNzBkyBDfVEkU5Fw7gafGhEq6m7Nrs876JgtaLDbuLB3gNh1ytNqMzYxHrHP7DCK6PI9bbp555hmkpKQAAJ5++mnExcXhnnvuQVVVFdasWeOLGomCXiCMtwGA2PAQRGgdgaa0vlnSWujyNhWWAwCmDU6WuhSioOJRy40oikhMTHS30CQmJmLjxo2+qo1INlzjbXp1kzbcCIKAHnFhOFphQGldM/p0j5S0HmpftcGE3afrAABTBnG8DZEnPGq5EUURWVlZHFtD5CFXt1RmgnTjbVx6xIYBAM7UseUmkH1dWAFRBIb2iEGq8++MiDrGo3CjUqnQt29f1NTU+K4iIhlyd0tJ3HIDAD3iHL8oS+ubpC6FLsE1S2raYLbaEHnK4zE3zz77LJYuXYqDBw/6piIimRFFMWDG3ABAWpyj9aiULTcBq6HZgu+OVQMApnK8DZHHPJ4tNX/+fDQ1NSE7OxtarRZhYa2bS2tra71ZH1HQqzWa0WiyQhCknQbu4uqW4oDiwJVXWAGzzY6+iZHolyTNDvJEwczjcLN69WrfVEIkU67BxKkxYQEx9drVLcUxN4Hrs/1nAQA/H5YqdSlEQcnjcLNgwQLfVEIkU6eqHWNbekm4p9T50pzhpkLfAovNziX9A0yd0ezukpo5LEXqcoiCUqfe1U6cOIHHHnsMt956KyorHfuefPnllzh06JBHz7Nt2zbMmjULqampEAQBGzZs6PC127dvh0ajwfDhwz2un8ifTgfINHCXhAgdtBoV7CJQ3tAidTl0gU2F5bDaRQxIjkJWIqfqE3WGx+EmPz8fQ4cOxf/+9z989NFHMBgMAIB9+/ZhxYoVHj2X0WhEdnY2Xn75ZY+uq6+vx/z583Httdd6dB2RFM46A4SrxURqKpXgHndTUscZU4Hms/1lAIBZ2eySIuosj8PNsmXL8NRTTyEvLw9a7bnlwCdNmoSdO3d69FzTp0/HU089hRtuuMGj6+6++27MmzcP48aN8+g6Iim4WkdSYkKlLsXNPaiY424CSo3BhO9POJbamDmUXVJEneVxuDlw4ECbYSQxMRHV1dXeqqtda9euxcmTJzvcSmQymaDX61t9EPlTWYMjQCQHULhJi+OMqUD0331nYbOLGNojJiCWDSAKVh6Hm9jYWJSVlV10fO/evejRo4e36mrTsWPHsGzZMvznP/+BRtOxsdC5ubmIiYlxf6Snp/u0RqLziaKIMnfLTWB0S4GrFAesD38sBQDcOMK376VEcudxuLnlllvw+9//HuXl5RAEAXa7Hdu3b8cjjzyC+fPn+6ZKADabDfPmzcPKlSvRr1+/Dl+3fPlyNDQ0uD+4dQT5U6PJiiazDQCQHB04LTfuVYoZbgLG0YpGHChtgEYl4HqOtyHqEo+ngj/zzDNYvHgx0tPTYbPZMGjQIHfweOyxx3xTJYDGxkbs3r0be/fuxX333QcAsNvtEEURGo0GmzZtwqRJky66TqfTQafT+awuoktxjbeJDQ9BmFb6NW5c3KsUs1sqYHy45wwA4JoBiegWyfcsoq7wONxotVq89tprePzxx3Hw4EEYDAbk5OSgb9++vqnQKTo6GgcOHGh17JVXXsGWLVvwwQcfIDMz06evT9QZri6pQGq1wXktN2UNzbDZRahVgtQlKZrVZsfHex1dUr8ckSZ1OURBz+Nw89133+HKK69Ez5490bNnzy69uMFgwPHjx92Pi4qKUFBQgPj4ePTs2RPLly9HaWkp1q9fD5VKhSFDhrS6PjExEaGhoRcdJwoU5QE4mBgAkqJ0UKsEWGwiKhtbEEjjgZTo2+PVqGw0IS48BJMGJEpdDlHQ83jMzaRJk5CZmYk//OEPKCws7NKL7969Gzk5OcjJyQEAPPTQQ8jJycETTzwBACgrK0NxcXGXXoNISmUBOA0cADRqlbs1ieNupPfBbkeX1PXZqdBquGI0UVd5/K/o7NmzePjhh5Gfn48hQ4Zg+PDh+L//+z+cOXPG4xefOHEiRFG86GPdunUAgHXr1mHr1q3tXv/HP/4RBQUFHr8ukb9U6F3dUoHXMuLaxPNUDRfyk1K1wYRNheUAgJtHd601nIgcPA43CQkJuO+++7B9+3acOHECc+bMwb///W9kZGS0OaCXSMkCteUGAPonO3abPlzGtZ+k9OGeM7DYRGSnx2JQarTU5RDJQpfaPzMzM7Fs2TI8++yzGDp0KPLz871XGZEMuGZLBdqYGwAY4Ao35Y1Sl6JYoiji7V2Orvd5Y7gGF5G3dDrcbN++Hffeey9SUlIwb948DBkyBJ9//rl3qyMKcoHccjMwxdFKcLicLTdS2XGyBqdqmhCp0+Dnw7i2DZG3eDxbavny5XjnnXdw9uxZTJkyBS+99BJ+8YtfIDw83DcVEgWpJrMVDc0WIEBbbvolRUEQgGqDGVWNJnSP4toq/vb2Lseior8YnooIncdvx0TUDo//NW3btg1Lly7F3LlzkZCQ4JuqiGTA1SUVqdMgKjRE6nIuEqZVI7NbBE5WG/FTmR7do7pLXZKiVDWasPGgYyubW8dwIDGRN3kcbrZv3+6bSohkxhVukqIDt0VkQEoUTlYbcbhcj6v7Mdz401v/K4bFJiKnZyyG9IiRuhwiWel0O2hhYSGKi4thNptbHb/++uu9URdR0AvEDTMvNDA5Gl8cKMfhMg4q9iez1Y43/3caALBwfIbU5RDJjsfh5uTJk7jhhhtw4MABCIIAURQBAILgWL7dZrN5v0qiIFSuD9yZUi4DnIOKf+KMKb/aeKgclc5xTtOHpEhdDpHseDxbasmSJcjMzERlZSXCw8Nx6NAhbNu2DaNGjbrkgntESlMewDOlXFzTwY9XNsJis0tdjmL8+/tTAIDbxvbkisREPuDxv6odO3bgT3/6ExISEqBSqaBSqXDllVciNzcX999/v2+qJApCZQG8xo1LWlwYonQaWGwiTlQZpC5HEfafqcee03UIUQuYN5YDiYl8weNwY7PZEBXl+N9eQkICzp49CwDo1asXjhw54v0KiYJUud6xZ1Mgt9wIgoABKa6Vitvumjpa0Yh/f38KD7+3Dyv/ewgmK7ueu+If204CAGYOTUFiVOD+bBAFM4/H3AwZMgT79u1DZmYmxo4di+eeew5arRZr1qxB7969fVMlURCq0JsAIOB/gQ1IjsYPp+rwU7kes9Gj1df+u+8s7n9nL5xD6wAAJbVNeOW2kexO6YSiaiO+POCY/n33xD5Sl0MkWx6/Oz322GOw2x1983/6059QVFSEq666Cl988QX+8pe/+KJGoqAjiiLqjI6ZhPERWqnLuaTBzv2M8o9UuScIAMDpGiOWf3QAogiMyYjHnVdlQqdR4eufKrHknb2wcoyOx9ZsOwm7CEwakIgBydxHishXPG65mTZtmvvzrKwsHD58GLW1tYiLi3PPmCJSOoPJCqvdERTiwgM73Fw3JBkr/1uIw+WN2HmyFuP6dIPZasf9b++FwWTF6Iw4vHXnWGjUKlyRlYC71u/BlwfL8Wr+Cdw3qa/U5QeNSn0LPtxzBgBwL1ttiHzKK+3K8fHxDDZE56lvcmy7oNOoEKZVS13OJcWGa3HjCEd31NrtRQCAZ788jH1nGhATFoLVt+RAo3a8VUzsn4inbhgCAHj9uyIYTFYJKw8u//yuCGabHaMz4jAqI17qcohkjZ3mRD5Q1+Tokgr0VhuX31zhWEgu76cKPP/VEfzLGXKeu2kYesS2XoTwlyPS0DshAvVNFryx47Qk9QabWqMZb+x03Kt72GpD5HMMN0Q+UOdsuYkND7w9pdqSlRiFq/t1hygCf/vmOADgwcn9MG1w8kXnqlUCFl+TBQB47duTaDJ3rfWmzmjG+h2ncOf63dhyuKJLzxWoXv/2JJrMNgztEYNr+idKXQ6R7DHcEPlAfZC13OC81hsAmJWdivuvzWr33F8MT0XP+HDUGs1463/FnX7N93eXYMwzX+OJTw4hr7ACL2w62unnClR1RrN70b77r+3LLnwiP2C4IfIB10ypuIjgaLkBgAl9u2P28FRMH5KM/7tp2CV/CWvUKtznbL15Nf8EjJ0Ye2Ox2bFq42FYbCKyEiMB55o6Zqu8ZmH9a3sRjGYbBqZEY/JAttoQ+QPDDZEPnOuWCp6WG5VKwOpbcvD3X41EaMjlB0HfMKIHMrqFo9pgdg9E9sSWw5WoNpjRPUqHL5dchahQx0rJxyvls1JyrdGMddsdrTZLrs1iqw2RnzDcEPnAuW6p4Gm58VSIWoUHp/QDAPwj/6S7taqj3vuhBABw44geCFGrMMi5ieehsw0+qFYaL319FI0mKwalRGPqoIvHLxGRbzDcEPmAq+UmmMbcdMasYakYmBKNRpMVr+af6PB1FfoWfHOkEgAwZ2Q6AGBwagwA4NBZvY+q9a8TVQa86RyP9NjMgVCp2GpD5C8MN0Q+4JoKHkzdUp2hUgl4dFp/AMC670+hUt/Soes++rEUdhEY1SvOPd7GtVJyoUzCzbNfHobVLmLSgESMz0qQuhwiRWG4IfKBenfLjXy7pVwm9u+O4emxMFnt+ODHM5c9XxRFvL/b0SU1d1S6+/jgHs5wU6aH3S62e30w2HmyBnmFFVCrBCyfPkDqcogUh+GGyAeU0nID587i88b0BAC8v/tMq/2p2vLVoQqcrDYiUqfBzGEp7uN9ukdCq1HBYLKiuLbJ53X7it0u4pkvfgIA3DI6HX2ToqQuiUhxGG6IfEBJLTcAMHNYCsK1ahRVG7H7dF2759ntIlZ/7VjLZuH4DETozm1vF6JWYUCyIwgE87ibT/edxf4zDYjQqvHA5H5Sl0OkSAw3RF5mttrdey7JfUCxS4ROg587W2Hedc6CastXh8pxuLwRkToN7rgq86Kvu8bdBOuMqRaLDf/31RHAuc1C9yid1CURKRLDDZGX1Tc7uqQEAYgOU0bLDQDcPNoxfubz/WVtbqjpaLU5BjhXQ26ry25QkM+Y+tf2IpTWNyMlJhS3X9lb6nKIFIvhhsjLXF1SMWEhUCto+u+InnHo3T0CzRYbNuwtvejrH+0txZGKRkTpNLijnV/851pugi/clDU0429bHPtyPTy1f8DvBk8kZww3RF7m3npBIV1SLucPLH5h0xFUNZrcXztdY8QfPz0EALh7Yh/EtDMWaWByNAQBqDaYUG0wtXlOoHrys0I0mW0Y0TMWN+b0kLocIkVjuCHysmDbEdyb5o/LwMCUaNQ1WfDYhgMQRRFmqx33v70XBpMVozPi8Nur2++uCdOq0T3SMU6lrL5ja+YEgvyjVfjiQDlUAvDU7KFcsI9IYgw3RF4WjDuCe4tWo8Lzc4ZBoxLw1aEKPPrBftz62k7sO9OAmLAQrL4lBxr1pd92kqJDAecqxsGgxWLDik8OAgAWjs/EIGfXGhFJh+GGyMuU3HID5zYKv5vUFwDw/p4z2HO6DioBWPXLYegRG3bZ613hpjxIws3qr4/hVE0TEqN0eHBKX6nLISIAmg6cQ0QeUHLLjcu91/RBXZMZ+hYLctJjMa5PN2Qldmwxu+QYR7dUR7dykNK+knqs2ebYU+vpG4YiKlSZgZYo0DDcEHlZnQJ2BL+cELUKf7x+cKeuTYoKjpYbk9WGRz/YD7sIXJ+diimDkqQuiYic2C1F5GXnuqWU23LTFUkxrnAT2LOl/rL5GI5UNKJbhLbTQY6IfIPhhsjLlDoV3FuSnWNuArlbaldRLV7Z6uiOemr2EMRH8O+aKJAw3BB5GbuluibQBxTrWyx48N0CiCJw08g0TB+a0oGriMifGG6IvKye3VJd4mq5qW+yoMVik7qci6z45BBK65vRMz6c3VFEAYrhhsiLRFFEfbNzR/AIttx0RnSYBjqN462pMsDG3Xy67yw+3lsKlQC8eHM2InWck0EUiBhuiLxI32KFzS4CHHPTaYIgIDkm8LqmSuub8f8+PgAAuG9SX4zsFS91SUTUDoYbIi9yrXETGqJCaAg3TuysQBt3Y7OLePi9AjS2WDE8PRb3T8qSuiQiugSGGyIvck0DZ6tN1yQF2Iypv2w+hp0naxGuVWP1zcMvu4UEEUmL/0KJvEjvHG8TE8bxNl2RHO1Ypbi8Qfpws/VIJf6y5RgA4OkbhiAjIULqkojoMhhuiLzIYLICAKJCOdC0KwKlW6q0vtk97fu2sT1xQ06apPUQUccw3BB5kaHFEW44i6ZrznVLSTdbqtlsw91v7EFdkwVDe8Tg8Z8PkqwWIvIMww2RFzU6W24iuYFil0g9W0oURTz64X4cKG1AfIQWr9w2ggPEiYIIww2RFxld4UbHX4RdkXxet5Qoin5//Ve2nsB/952FRiXg77eNQHp8uN9rIKLOY7gh8iKDid1S3tA9yjGg2Gy1o8E5SNtf8gor8PymIwCAlb8YjLG9u/n19Ymo6xhuiLyo0T3mht1SXREaonbvzeXPrqmjFY144J29EEXg1z/rhdvG9vLbaxOR9zDcEHmRu+WGs6W6zD1jyk/TweuMZty5fjeMZht+1jseT8ziAGKiYMVwQ+RFhhZHF0oUu6W6zJ8zpkxWG+56YzdO1zQhLS4Mr9w2EiFcqI8oaPFfL5EXGU2OXazZctN1Kc4ZU2cbmn36OqIo4tEP9uOHU3WI0mnwr4WjER/BFaaJghnDDZEXuaaCR7DlpsvS4sIAACW1vg03L359DJ8UOGdG/Wok+iVF+fT1iMj3GG6IvMhgcnRLcbZU17mmX5fUNfnsNT7ccwZ/2ezYWuGp2UNwZd8En70WEfmPpOFm27ZtmDVrFlJTUyEIAjZs2HDJ8z/66CNMmTIF3bt3R3R0NMaNG4evvvrKb/USXY5rhWJuv9B17nBT65tws+NEDZZ9tB8AcM/EPrhlTE+fvA4R+Z+k4cZoNCI7Oxsvv/xyh87ftm0bpkyZgi+++AJ79uzBNddcg1mzZmHv3r0+r5WoI7jOjfekxznCTbm+BSarzavPfbzSgN++sRsWm4iZQ1OwdGp/rz4/EUlL0nfg6dOnY/r06R0+f/Xq1a0eP/PMM/jkk0/w3//+Fzk5OT6okKjjTFYbLDbHarocUNx1CZFahIWo0Wyx4Wx9CzK9tBt3jcGERet+gL7FipyesXhhbjZUKsErz01EgSGox9zY7XY0NjYiPj6+3XNMJhP0en2rDyJfcHVJAUCEluGmqwRBOG9QsXe6plosNtz1xh4U1zYhPT4Mr80fxT2jiGQoqMPN888/D4PBgLlz57Z7Tm5uLmJiYtwf6enpfq2RlMPVJRWuVUPNlgCv8OagYteU7z2n6xAdqsHahaOREKnzQpVEFGiCNty89dZbWLlyJd577z0kJia2e97y5cvR0NDg/igpKfFrnaQc57ZeYKuNt/R0Dyru+nTwlzYfw6fOzTBf/dVIZCVyyjeRXAXlu/A777yDO+64A++//z4mT558yXN1Oh10Ov7vjHzPyK0XvM5b3VKfFJRi9dfnpnyPz+KUbyI5C7qWm7fffhu/+c1v8Pbbb2PmzJlSl0Pk5uqW4tYL3uONbqkfi+uw9APHlO+7ru7NKd9ECiDpu7DBYMDx48fdj4uKilBQUID4+Hj07NkTy5cvR2lpKdavXw84u6IWLFiAl156CWPHjkV5eTkAICwsDDExMZJ9H0Tgppk+4ZoO3tmWmzN1Tbhr/W6YrXZMHpiE3183wMsVElEgkrTlZvfu3cjJyXFP437ooYeQk5ODJ554AgBQVlaG4uJi9/lr1qyB1WrF4sWLkZKS4v5YsmSJZN8DkYtrzA1nSnlPeryjW6quyeIOjx3V2GLB7et2o9pgxsCUaLx0y3AO9CZSCEnfhSdOnAhRFNv9+rp161o93rp1qx+qIuocttx4X1RoCGLDQ1DfZEFJbRMGpkR36DqrzY7fvb0XRyoakRilwz8XjOJ+X0QKEnRjbogClZFjbnyiZye2YXj6i5+w9UgVQkNUeH3BKKTGhvmwQiIKNAw3RF7ingrOlhuvco+7qevYdPA3dp7G2u2nAAB/njscw9JifVofEQUehhsiLzm3r1SI1KXISlp8x6eDbztahT9+eggAsHRaf8wYmuLz+ogo8DDcEHmJwb2IH5fz96aOzpg6VtGIxW/+CJtdxI0jeuDeiX38VCERBRqGGyIv4YBi3+jVzRFuDpc3tjsBocZgwqJ//4BGkxWjM+KQe+NQCAJnRhEpFcMNkZewW8o3RvaKg1ajQml9M45WGC76umszzJLaZvSMD8ervxoJnYatZ0RKxnBD5CXnwg1bbrwpXKvBlc7tEvIKy1t9zW4XsfS8zTD/tXA0unEzTCLFY7gh8hLXmJsodkt53ZRBSQCAvJ8qWx1f/fVR/LfVZpiRElVIRIGE4YbIS1wtN1wszvuuHZgIANhXUo8KfQsA4MM9Z/CXLY7tW565cSg3wyQiN4YbIi+w20V2S/lQYlQohqc71qvZ/FMlvjlSiWUfOTbDvGdiH8wdlS5xhUQUSPguTOQFTRab+3N2S/nGlEFJKCipx1+3HEOFvgV2EZgxNBlLp/aXujQiCjBsuSHyAtd4G41KgE7Df1a+4Bp3U9bgCDZzR6XhxZuHQ8XNMInoAvwvJpEXGEwWwLnGDddX8Y2+iZEYlBKNIxWNeGzmQCwcn8F7TURtYrgh8gL3vlIcb+MzgiDg3d/+DI0tVm6ESUSXxHdiIi/gYGL/iAoNQVQoF0kkokvj4AAiLzAy3BARBQyGGyIvcHdLcaYUEZHkGG6IvIDdUkREgYPhhsgLuPUCEVHgYLgh8gL31gtahhsiIqkx3BB5gbtbii03RESSY7gh8oL6ZsciftGcpkxEJDmGGyIvqG8yAwDiI7RSl0JEpHgMN0ReUGd0tNzEhrPlhohIagw3RF7garmJC2fLDRGR1BhuiLygluGGiChgMNwQdVGLxYYWix0AEBvBbikiIqkx3BB1UZ2z1UajEhDFFYqJiCTHcEPURecGE2shCILU5RARKR7DDVEXnRtMzC4pIqJAwHBD1EUcTExEFFgYboi6qK6Ja9wQEQUShhuiLqo3suWGiCiQMNwQdZGr5SaOWy8QEQUEhhuiLuKAYiKiwMJwQ9RFdRxQTEQUUBhuiLqolgOKiYgCCsMNURe5u6U45oaIKCAw3BB1UZ2RY26IiAIJww1RF1htduhbrADH3BARBQyGG6IuaGi2uD+PCWPLDRFRIGC4IeoC10yp6FANNGr+cyIiCgR8NybqAi7gR0QUeBhuiLrANZg4luNtiIgCBsMNURfUu1puOFOKiChgMNwQdYFrzE08W26IiAIGww1RF9Q2sVuKiCjQMNwQdUG9kd1SRESBhuGGqAtc3VKxnC1FRBQwGG6IuoADiomIAg/DDVEXcEAxEVHgYbgh6oI6DigmIgo4DDdEnXS80oAa5yJ+idE6qcshIiInhhuiTnox7yhEEZgyKAkJkQw3RESBguGGqBMOljbg8wNlEATg4an9pC6HiIjOI2m42bZtG2bNmoXU1FQIgoANGzZc9pqtW7dixIgR0Ol0yMrKwrp16/xSK5GLxWbH85uOAACuz07FgORoqUsiIqLzaKR8caPRiOzsbCxatAg33njjZc8vKirCzJkzcffdd+PNN9/E5s2bcccddyAlJQXTpk3zS83tsdlFlDU0QxQdj0URECE6/wREUXT+CaDV8XPn2UWxc9ef9zX3Ne1cj4vOafu57c7zztXiOM9FvPDrznNwwWs5n7FVvRcecz1Jq69f+Lid12/99dbH3M/Txjnnf19mq93xYbPBZLHDbHM8NlntMJqsMDg/jCYrGlusMFntAAC1SsCDk9lqQ0QUaCQNN9OnT8f06dM7fP6rr76KzMxMvPDCCwCAgQMH4rvvvsOLL74oebipMZhw5apvJK2B/EcQgLuu7o2MhAipSyEiogtIGm48tWPHDkyePLnVsWnTpuGBBx5o9xqTyQSTyeR+rNfrfVKbIAjQaVQQBECA4PzTcVwAgPMfX/A1Qbjg8wuuh/ucyz+36rxzcNFrtL4erV7z3Lkq5wHX1xwvIbif4/zH53//gvvzC2p3P49wwdfbf27X61/4PZx7vnPPDef31ZHnxvnfv/M8rVoNrUYFnUbl/tP1eYROgwidBlE6DSJDNYjQahAV6jgWouaQNSKiQBRU4aa8vBxJSUmtjiUlJUGv16O5uRlhYWEXXZObm4uVK1f6vLbuUTocearjrVBERETkG7L/r+fy5cvR0NDg/igpKZG6JCIiIvKhoGq5SU5ORkVFRatjFRUViI6ObrPVBgB0Oh10Oq5BQkREpBRB1XIzbtw4bN68udWxvLw8jBs3TrKaiIiIKLBIGm4MBgMKCgpQUFAAOKd6FxQUoLi4GHB2Kc2fP999/t13342TJ0/i0UcfxeHDh/HKK6/gvffew4MPPijZ90BERESBRdJws3v3buTk5CAnJwcA8NBDDyEnJwdPPPEEAKCsrMwddAAgMzMTn3/+OfLy8pCdnY0XXngBr7/+uuTTwImIiChwCOL5K6MpgF6vR0xMDBoaGhAdzZVliYiIgoEnv7+DaswNERER0eUw3BAREZGsMNwQERGRrDDcEBERkaww3BAREZGsMNwQERGRrDDcEBERkaww3BAREZGsBNXGmd7gWrNQr9dLXQoRERF1kOv3dkfWHlZcuGlsbAQApKenS10KEREReaixsRExMTGXPEdx2y/Y7XacPXsWUVFREATBa8+r1+uRnp6OkpISbutwGbxXnuH96jjeK8/wfnUc75VnfHG/RFFEY2MjUlNToVJdelSN4lpuVCoV0tLSfPb80dHR/MHvIN4rz/B+dRzvlWd4vzqO98oz3r5fl2uxceGAYiIiIpIVhhsiIiKSFYYbL9HpdFixYgV0Op3UpQQ83ivP8H51HO+VZ3i/Oo73yjNS3y/FDSgmIiIieWPLDREREckKww0RERHJCsMNERERyQrDDREREckKw40HXn75ZWRkZCA0NBRjx47Frl27Lnn++++/jwEDBiA0NBRDhw7FF1984bdapebJvVq3bh0EQWj1ERoa6td6pbJt2zbMmjULqampEAQBGzZsuOw1W7duxYgRI6DT6ZCVlYV169b5pdZA4On92rp160U/W4IgoLy83G81SyU3NxejR49GVFQUEhMTMXv2bBw5cuSy1yn1fasz90up711///vfMWzYMPcCfePGjcOXX355yWv8/XPFcNNB7777Lh566CGsWLECP/74I7KzszFt2jRUVla2ef7333+PW2+9Fbfffjv27t2L2bNnY/bs2Th48KDfa/c3T+8VnKtYlpWVuT9Onz7t15qlYjQakZ2djZdffrlD5xcVFWHmzJm45pprUFBQgAceeAB33HEHvvrqK5/XGgg8vV8uR44cafXzlZiY6LMaA0V+fj4WL16MnTt3Ii8vDxaLBVOnToXRaGz3GiW/b3XmfkGh711paWl49tlnsWfPHuzevRuTJk3CL37xCxw6dKjN8yX5uRKpQ8aMGSMuXrzY/dhms4mpqalibm5um+fPnTtXnDlzZqtjY8eOFX/729/6vFapeXqv1q5dK8bExPixwsAEQPz4448vec6jjz4qDh48uNWxm2++WZw2bZqPqws8Hblf33zzjQhArKur81tdgaqyslIEIObn57d7jpLfty7UkfvF965z4uLixNdff73Nr0nxc8WWmw4wm83Ys2cPJk+e7D6mUqkwefJk7Nixo81rduzY0ep8AJg2bVq758tFZ+4VABgMBvTq1Qvp6emX/B+A0in156qrhg8fjpSUFEyZMgXbt2+XuhxJNDQ0AADi4+PbPYc/X+d05H6B712w2Wx45513YDQaMW7cuDbPkeLniuGmA6qrq2Gz2ZCUlNTqeFJSUrt99+Xl5R6dLxeduVf9+/fHv/71L3zyySf4z3/+A7vdjvHjx+PMmTN+qjp4tPdzpdfr0dzcLFldgSolJQWvvvoqPvzwQ3z44YdIT0/HxIkT8eOPP0pdml/Z7XY88MADuOKKKzBkyJB2z1Pq+9aFOnq/lPzedeDAAURGRkKn0+Huu+/Gxx9/jEGDBrV5rhQ/V4rbFZwCz7hx41ol/vHjx2PgwIH4xz/+gSeffFLS2ii49e/fH/3793c/Hj9+PE6cOIEXX3wRb7zxhqS1+dPixYtx8OBBfPfdd1KXEhQ6er+U/N7Vv39/FBQUoKGhAR988AEWLFiA/Pz8dgOOv7HlpgMSEhKgVqtRUVHR6nhFRQWSk5PbvCY5Odmj8+WiM/fqQiEhIcjJycHx48d9VGXwau/nKjo6GmFhYZLVFUzGjBmjqJ+t++67D5999hm++eYbpKWlXfJcpb5vnc+T+3UhJb13abVaZGVlYeTIkcjNzUV2djZeeumlNs+V4ueK4aYDtFotRo4cic2bN7uP2e12bN68ud0+xnHjxrU6HwDy8vLaPV8uOnOvLmSz2XDgwAGkpKT4sNLgpNSfK28qKChQxM+WKIq477778PHHH2PLli3IzMy87DVK/vnqzP26kJLfu+x2O0wmU5tfk+TnymdDlWXmnXfeEXU6nbhu3TqxsLBQvOuuu8TY2FixvLxcFEVR/PWvfy0uW7bMff727dtFjUYjPv/88+JPP/0krlixQgwJCREPHDgg4XfhH57eq5UrV4pfffWVeOLECXHPnj3iLbfcIoaGhoqHDh2S8Lvwj8bGRnHv3r3i3r17RQDin//8Z3Hv3r3i6dOnRVEUxWXLlom//vWv3eefPHlSDA8PF5cuXSr+9NNP4ssvvyyq1Wpx48aNEn4X/uPp/XrxxRfFDRs2iMeOHRMPHDggLlmyRFSpVOLXX38t4XfhH/fcc48YExMjbt26VSwrK3N/NDU1uc/h+9Y5nblfSn3vWrZsmZifny8WFRWJ+/fvF5ctWyYKgiBu2rRJFAPk54rhxgN//etfxZ49e4parVYcM2aMuHPnTvfXJkyYIC5YsKDV+e+9957Yr18/UavVioMHDxY///xzCaqWhif36oEHHnCfm5SUJM6YMUP88ccfJarcv1xTlS/8cN2fBQsWiBMmTLjomuHDh4tarVbs3bu3uHbtWomq9z9P79eqVavEPn36iKGhoWJ8fLw4ceJEccuWLRJ+B/7T1n0C0Ornhe9b53Tmfin1vWvRokVir169RK1WK3bv3l289tpr3cFGDJCfK0F0/KUSERERyQLH3BAREZGsMNwQERGRrDDcEBERkaww3BAREZGsMNwQERGRrDDcEBERkaww3BAREZGsMNwQERGRrDDcEBERkaww3BAREZGsMNwQERGRrDDcEFHQW79+Pbp16waTydTq+OzZs/HrX/9asrqISBoMN0QU9ObMmQObzYZPP/3UfayyshKff/45Fi1aJGltROR/DDdEFPTCwsIwb948rF271n3sP//5D3r27ImJEydKWhsR+R/DDRHJwp133olNmzahtLQUALBu3TosXLgQgiBIXRoR+ZkgiqIodRFERN4wcuRI3HTTTZg6dSrGjBmDU6dOIT09XeqyiMjPNFIXQETkLXfccQdWr16N0tJSTJ48mcGGSKHYckNEstHQ0IDU1FRYrVasX78eN998s9QlEZEEOOaGiGQjJiYGv/zlLxEZGYnZs2dLXQ4RSYThhohkpbS0FLfddht0Op3UpRCRRNgtRUSyUFdXh61bt+Kmm25CYWEh+vfvL3VJRCQRDigmIlnIyclBXV0dVq1axWBDpHBsuSEiIiJZ4ZgbIiIikhWGGyIiIpIVhhsiIiKSFYYbIiIikhWGGyIiIpIVhhsiIiKSFYYbIiIikhWGGyIiIpKV/w/5qXuLOBy97QAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "g = dens.g\n", + "ax.plot(g.y[g.jlo:g.jhi+1], dens_avg)\n", + "ax.set_xlabel(\"y\")\n", + "ax.set_ylabel(\"average density\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.13.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}