Skip to content

Instantly share code, notes, and snippets.

@will-henney
Created January 22, 2020 05:05
Show Gist options
  • Save will-henney/6982617a9e6bb6d14774338db454122b to your computer and use it in GitHub Desktop.
Save will-henney/6982617a9e6bb6d14774338db454122b to your computer and use it in GitHub Desktop.
Demo of contours for Tere
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from matplotlib import pyplot as plt\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Generate grid of parameter values (`R_coma` and `L_crit` in Tere's example)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"X, Y = np.meshgrid(\n",
" np.linspace(0.0, 0.5),\n",
" np.linspace(0.0, 0.5)\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Fake data generated on this grid. Not exactly the same as Tere's, but close enough"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"data = 1.0/(1.0 + 0.1*(X**2 + Y**2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now generate the figure with image, labelled contours, and color key. "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAAEUCAYAAABqA/5qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2dd5hU5fX4P2dmC71LB0FFpNgRxK6oAaMSLLFjLDEaa4yJJibRdE35Gb9qLLHFHhuKCiL2AiggSAfpLH0pCyzLlpnz+2MG2X3vC/PuMDuzM/t+nuc+MHfOfe+5d2bPvPec95wjqorH4/F46o5QphXweDyeXMcbWo/H46ljvKH1eDyeOsYbWo/H46ljvKH1eDyeOsYbWo/H46lj0mpoRWSoiMwXkYUicofl/ZNEpEREpse336VTP4/H46kL8tJ1IhEJAw8BpwFFwGQRGa2qcwzRz1T1zHTp5fF4PHVNOme0A4GFqrpYVSuAl4DhaTy/x+PxZIR0GtouwIpqr4vi+0wGi8g3IjJWRPqlRzWPx5ONiMiTIrJORGbt5n0Rkf+LuytniMgR1d7boyszlaTNdQCIZZ+Z//s1sK+qbhORM4A3gF6BgUSuAa4BkIKCI/Pbt0+1rnaSyFa2XXS6zp3scVKH52rVpBGdWrcgJMLm0h0Ub9lGRVWkVmPYZJx0NmSSOSZ2nCaUCZzLOrYhFHU5l+Vk5i6rjDmO5WTRmjLW9Pxq+3ZQSoWWp+wrniRPAw8Cz+zm/WHEbEgvYBDwMDCoFq7MlJBOQ1sEdKv2uiuwqrqAqm6p9v8xIvJvEWmnqsWG3GPAYwCF3bppl1t+lnJlrX+EyfwxO8g4GQ7rH6HxOmo5W0DGok80sUzwXInHCUXsMi2aNeGKUwZwwbGHkhcO8c7kuTw67ktWbiixXmcoUvPkYhnX3Gc7txjjhKqSO1eo0hwn+AGGKmteiHlum0yoIngyqYzs8TUAlTUvRCoqE8poeXlARMsrar7eYZGp3CXzpX4QPE+aUdVPRaTHHkSGA89o7Fdjkoi0EpFOQA/irkwAEdnpyqwTQ5tO18FkoJeI9BSRAuBCYHR1ARHpKCIS///AuH4b0qijJw1s2Ladf4z+lKF/eoIXP53O0CN78+ZvLudX559MuxZNM62eJ7fYncvS1ZWZEtI2o1XVKhG5ARgHhIEnVXW2iFwbf/8R4DzgOhGpAsqACzVReTGBaMEukaQfewP6ughZ1HF40kNrzjxtOgeOc5n1Wk5mzjyts1VjJmx7qnQZJ6Cj7We82nHrd2znb29/wlOfTuEnpw3i3GMOZvigfrzw2XSe/GAyW8riM6qQcb9Clus0zqWWc4eMcWw3Xg2ZkO0azAcH24OEoZA5U7Ydp5ZxAipLUMi2L5FMSp73K5Mb5XsnN9UNG22PO0GmziifDeyotuux+BOtK7tzWbq4MlNGOl0HqOoYYIyx75Fq/3+QmL/F04BYv6WUP732IU9/PJXrTx/MFScP4PzBB/Po+C958bNviETc/ig92UHxxghfjuvqJJvfadEOVR2wF6fbncuyYDf76wSfGeapNxRtKOFXz7/L+f94jhnLVvOL4Sfy5h0jOe2wQDzUk9UoEY06bSlgNDAyvvrgaKBEVVfj4MpMJWmd0Xo8LixYXcx1j73B4N7due3sE/jHlWcyfckq/jHqE2YuXZNp9Tx7iQLRFD2li8iLwElAOxEpAu4C8uG7p+UxwBnAQmA7cEX8PasrMyVKWfCG1lNvmTh/Oef/43lGHNWP6884huduvYi3vprD/73xGetLSjOtnidJFKVSU+MOUtWLEryvwPW7eS/gyqwrst/QiqL5ux4xbL+TYosyJCLJYFgwiGUJXhhC9oCZeYxFJOoQVDODWLaYjLF+0rqUzDyXQ8DMttosYcDMHAPltSmzGDNjPj8eMpDLTzqCIYcewGPvfcWzH39NZdx/awbDbAEz81xqCSKZwS+1jGMeZw1imddlWUpmBszCQZHARxEI6IGbA9AhGFbbvxKpSj6klqoZbbbgfbSerKCsopL/G/sFw+99hkkLlnPLWcfx2u2XMbh390yr5qklCkRQpy1X8IbWk1UUbSzhlife4tpHXgfg0evO5W+Xn0H7ln79bTYRRZ22XMEbWk9WMmHeMs6991keHDOBk/vvz6g7L+eSkw4nbHu09tQrFIioOm25Qg74aIG8PX8gSbVUT8avC4555+ZK9cTj2IYxfaJOMhbHaSAZwerrNVJTXZIaLEYvoR/Xdi8C/tfYv+VEeOSjL3lnxjzuHH4yvzz3JM4a1Je7Xx7P3GXrgtdgjBOyJiyYMpbvQSAf2SaT4DXBWU7EMu8JnN/mVw4OnRQOKtdk+974aBsWfkbryXpWbCjhuv+8wa1Pv03b5k144ZaL+NnZx9MoP/vnEbmIqlLhuOUK3tB6cobxM77lB/c+w6gvZ3PFkAG8cvulDDjALQPJkz5i62jdtlzBG1pPTrF1Rzm/f+V9rn7wVUIiPHnj+dx5/ik0LsjPtGqe7xAijluu4A2tJyf56tsVnHvPszzz0VTOP+YQXr/jMgb08rPb+oASK33rsuUK2e/EcgiGJYWLf8ihMpc9qcGMdDn8crtUy7LImKeyBbECSQ3JVvgyz2WrspUgYGar/2oGzMzkBNs+CcF2qvjbO58yfvZC/nTB6Txx8/m88Ok0/vXW5+yorCJqq4Rl6myrzBVKnLAQ+BpYZYxKYdZsktp/t23VvJKZUSVMarB9EI7k0mzVBT+j9eQ805au4tz/9xzPfzKNi084nJd/cQl9u3XItFoNFgUqNeS05Qq5cyUezx7YUVnFvaM+5qoHX6EwP49nb7mAq04/yr58y1OnxDLDGpaPNutdB20Km3Bhr0PiC5yjRKJRKqIRdlRVUR6pYkekirKqSrZWlrO1ooJtleU5tRDaUzsmLyzivL89x2/OH8JNZx3HsX16cOcz77J609ZMq9ZgUMS6ZjiXyXpD26VZC+49dlitjtleVcGm8rL4tp0NO7azcUcZa8u2srZsG2u3b2Nd2VZWlW6lLFKt/5JToRnDf+fkx3VouGdJNHDy9Tr4cc199qSGxAkLRBL7LoM+WmNcW5JDZM+vIbiO35YgtvOySqrK+cWLY/h07hLuPPdkXv7VpfzxlQ94d/qCwAw3avGRmmNblyEFEg1sMgleW7E5qBNKEDX+1J3MnK1jR/V9e/E0EE02IShLyXpDO3fzOo4Z9QAhCREWIU9CFITDNArnUxjOozCcR5O8fJrlF9Isv5Dm+YW0KCikVUFjWjdqTJvCJnRr1op2jZrSLL8wMH7xjlKKtpVQtG0zK7aVsHjLRhZv2cCSLRvZWF6WgSv2pIq3p85l2pKV3HPJMP4+8vsMOrA7f3/lY3ZUWkpteVLGTtdBQyLrDW1VNMLq7Xv/2KcqNM0roH3jZnRo0owOjZvTpWkLujZrSdemrejbpiOnd+tNQXhXMbuS8jIWlmxgfsl6FmwuZv7m9SzYVMyGHdv3Wh9Peli5cQs/euhlrh96DD8+dSCHdu/EL556hyVrN2ZatRxGiORQoMuFrDe0qaS0qoIlWzeyZKv9jyyM0KVpS/Zr0YaeLdqwX4u29GrZlmHdD+LiXo2/k1tduoVZG9cye+NaZm5Yw4ziNawr25auy/DUkkhU+b8xXzBlURF/vWQoL952MX96+QPenjw306rlJApUWqvv5i7e0NaCiCrLt21m+bbNfLxqcY339mnUlANb7cNBrdrTr00H+rftwCld9iccrwC9clsJ04pXMW39KqatW83MDWuoiPqmg/WJCfOXccG9z/HXy8/gz5cN5fD9uvC3lz+kosp/TqlE1c9osw4RCOfVLitarY74xJGugES1cTZUbWNi8TYmFi/5bl/jcD4HtWrPIW07c3jbzhzergtn9ugDwI6qSqYVr+KrdcuZtHYFX69fSXmkyjp2TGeLfmbQykHGeunm5MKWkmMGumwBMzOw5RAwCwTiHFqJW7sMGLbQ1t0hkNRgkVmzvZSrH32VG4Yew9VDBtK3e3tuffJtVm3cskvHQOAtcStxt+7WDn5Lq0hio2VKqKZgRrk3wTDvo/WkirJIJV8Xr+Tr4pU8Hd/XrlFTjmjbhaPad2NQ++7c0P9Ybj4kRHmkisnrVvD56qV8tnoJczasy6Gyx9lFJKrcP+YLvlm2mj9f+D1euu1ibv/vGCbOX55p1XKCWDDMz2g9dUjxjlLeK1rAe0ULAGieX8gRbbtxTMd9Ob5zT+444mTu4GQ27NjOp6sW81HRIj5ZtYTN5TsyrHnD4+PZi7nony9w31Vn8fC15/DAO1/wxPuTM61WDuBdB540s7WynI9XLeLjVYvg65iv99hOPTih036c2GU/RuzXn0g0ypT1K/lgxULGLV/A0q2bMq12g2FFcQmX3fcSd194GjefdRwHdtmH3z/3HmUVfglYssTKJHpDm12Ikpe/y0Fn97/WHqeaMpZzuXTBNY/T8K6DNkS2MrpoJqNXzCIkwqFtOnNy5wM4pfMB/HrAyfx6wMnM37yecSvmM275AmZvWrtrIJekBpvzMlC1JSgSmIDY/LimP9j0x2Lx4waSERJ3K4g6+HFthVUC53JIIhCBUq3iFy+OZc7q9fzsjOPo0b41tzwx+rtsMltxGrMYjEtRGbUWkEnm+2wxYg7f58BR1sJB1RMWaqNTdVWEilT4iLOI7De0OUxUlWkbVjJtw0r++c0ndGnagtO69mZot95c3+8Ybjr4OJZs2cjopXMYvXQOizb7tZ91yVMfT2HhmmLuvWQYL956MT978i2mLVmVabWykmgDcx00rKvNclaWbuHp+ZO58P3nGPj6/3H7pDGsLC3hhv7H8MHZ1/DOmVdwTb+BtG/cLNOq5iyfzVvKpfe9xJaych6//jzOHNAn0yplHTuDYS5bruBntFnKxvLt/G/hN/xv4Tfs07gpZ3bvw9k9+vLrI0/h9sNP4rPVS3h10UzeW/Yt5X69bkpZum4Tl9z3IvddcSZ/uXQo3du24uExEzOtVtagCBFf68CTbawvK+Wp+VN4au5UejRvzTn79+fc/fvz4Ak/oKRiB6MXz+GFBdOZsynYGdaTHFvLyrn20VH89odDuHbY0XRr15K7Xxjvkxsc8cGwLEOAvLy9/3K7BNECQaykA2ZGhShrEMvYYVmkH6iApMLyHRv41+xPuH/2Jwzu0INzex7M+QcczGUHHcE3xat44dvpjF46h7KqyhrH1RjGFjAzb7FLwMwStAokLJgBIVugyzjGpXqX2KpumffQsnDAjP/sqQpYBVF+++p4lq3fzC1nHkenti24+cnRbC7dgXmDQpZoVDDWlKKkBuv30sGwGV+6kCXgWSPRIcmEBVWhsoEFwxrWz0oDQoEJa5dy64S3GDTqAe6e/B6N8vK5d/AZfHnuDfx2wBC6N2uVaTVzgic+mMxt/32Hft068OxNF9KlTYtMq1SvUSCiIactV8idK/Hsli0VO/jvgikMfec/nPvuM3y0chEjex/Jxz+4lv+cdB7HdNw30ypmPeOmL+Dqh1+jddPGPHvzhRzYpV2mVarXNLRgWO5ciceJqetXcvPnozn29X/z4MwJHL5PZ144/SLGnnUl5+zXn/yQ/0oky/Qlq7j8wZeJRKI8efMPOcp33bWiCFF123IFsRYrySKa9Oqkvf91ZcrHdfGt2jC/HNYGCwFfb2IfbTQaNIBBmcTtWE39CkJhzuren6sOGshBrdqzevsWnpw7mZcWfsPWyvJqB7r4cR1aDxjHmd11A75gHH20poxtsb1DckSoypRJPI55TMeWzXjk6nPo3q4lv3ruXd775tuATOw4I3nDKrPnYwDClYllQhWG/7UyeIPCFTX3iU2mfNfFT5r1KFtKV9baGnbv30J//upAJ9lb+nwwVVUH1PYc9Q0/fWngVEQjvLp4BsPGPM6PPnqJJVs2cueRQ/hixE/5+aEn0KqgceJBPDVYU7KNyx/4H7OWr+XvI7/P+ccckmmV6hVKLGHBZcsV0nolIjJUROaLyEIRuWMPckeJSEREzkunfg2dT1Yv5pIPX+CsMU/x+eql3HjwsXw+4jp+ediJtCn0Brc2bNlezjWPvMZnc5fwu/OHMPKUIzOtUr1Bia06cNlyhbQZWhEJAw8Bw4C+wEUi0nc3cvcC49Klm6cmMzeu4aefjeL0t/7DhysXcW2/wXx2znXcccRJtCpslGn1sobyygi3PPUW706bz89HnMB1ZwzOtEr1hobWbjydM9qBwEJVXayqFcBLwHCL3I3Aa4BfXZ9hFpQUc9Pnb3LaW//hvRULuKbfID4bcR23HHoczfILMq1eVlAViXL7s2MZNXEW1w47mttGnJBplTKOqjQ410E6Exa6ACuqvS4CBlUXEJEuwAjgFOCo3Q0kItcA1wAUtG9Bo/w9l6xLJnqZTODLdpytyJVLMCwQVAsHAxNmgMzaYMGQCQSfCFaRihqr9JdsX8+tk0bz7zkTuPWQE7jl0OO4/KAjeXT2JJ6aP2VXZwjjOFv1rkBHBYdEiGBb7sTVsgKBOes4tiQCI6hmGSZQyCwoUkMmgnLXK+MprazkslOOpHHjAv7wyvuBcyWb1CCBDz5x23mxfOfM5VS2B/ca3+e9sIO5tEbWhXRerUvay7+A21UDf341D1J9TFUHqOqA/Jbed5guFm4p5qefv85ZY59kevEq7jjiZD46+yeM6Nk/hx7y6gZVuHfUxzw2/kvOG3wwd//wtL3pBJPVxOrRitPmQqLYj4i0FpFRIjJDRL4Skf7V3rtZRGaJyGwRuSV1V1mTdM5oi4Bu1V53BcwacwOAl+K1RNsBZ4hIlaq+kR4VPS7M2rSGKz56mUHtu3HnkUO479izuOKgAfx5yod8uXZF4gEaMA+MmUBVJMpPhw5GFH7/4ninVO5cQhEqo6kJdFWL/ZxGzMZMFpHRqjqnmtivgemqOkJEDorLD4kb3B8Tc2tWAO+KyDuq+m1KlKtGOme0k4FeItJTRAqAC4HR1QVUtaeq9lDVHsCrwE+9ka2/fLluBcPHPs0tX4ymbaMm/O97l/DYSef41N4EPDxuEv9+dyIjBvfnrosa5sw2hZlhLrGfvsAHAKo6D+ghIh2APsAkVd2uqlXAJ8RclyknbTNaVa0SkRuIrSYIA0+q6mwRuTb+/iPJjCsCBeHaFZVJtguD6Te1+mgTHGM7f7K+3qjhY7SNY3YjsCU17K0fd/SKmYxbNZcrew3iun6DGT/8ah6f+yUPzZzI9qqaq+kDiQ6mlbH4VoPFagIiARmrG9dMagiKBP2kVjez6Vu1yCS4zIc+nATAT4cORgXufnm89VHZ1MfmDw76W22BgUTHWHzPNj9q9XGSLSpDrbK+2onIlGqvH1PVx6q9Thj7Ab4BzgE+F5GBwL7EnqhnAX8WkbZAGXAGMIU6IK3Vu1R1DDDG2Gc1sKr6o3To5EkN5ZEq/j1nAq8umcHth53M9f2P5ZyeB/OXrz/kraVzM61eveThcbuMbVSVP77wfoY1Sh+1KJNYnCAzzCX2cw9wv4hMB2YC04AqVZ0rIvcC44FtxAxynTSDa1ihP0+ds65sGz+f+BbnvPdfineU8sDxP+DZUy9k3+atM61aveThcZO+C5A1lKVfqhBRcdocSBj7UdUtqnqFqh4GjAT2AZbE33tCVY9Q1ROAjUDK/bPgDa2njphWvJLhY//Lb78cx6FtO/HeWVdzff9jfNEaCw+MmcBzn07jslOObDBJDSksKpMw9iMireLvAVwNfKqqW+LvtY//252Ye+HFFF1iDbK+8Len/hJV5dkFXzNuxQJ+N+BUfnH4iQzv2ZdfTRzHlHVFmVavXvG3Nz6mWX4+1w47mtIdFTzz4dRMq1Rn7EzBTclYbrGfPsAzIhIB5gBXVRvitbiPthK4XlU3pUQxg6w3tIJSGK6dWyXpYBiJg1guv8Lm+atslbkcxo0YETLbdUWMqEzE7MtNHQbM4ucqrtzKTRNH8caymfx+wPd4bdilPLtgKvdO+4ht5VXGGIFhg0kOFhGnJZeBhAXbqVwSFswAVeIqYLa/tOrLuhS465X3adQon5+POIHSqkpemTAjoKQ1qSGwyyVhISgS/NIFRSS8a+xkV6XFisqkbqlFotiPqk4Eeu3m2ONTpsgeyHpD68kePly1kAlvLePnh57IFQcdxUmd9+eXX4xh4trlmVatXhBV5VfPv0ujgnx+c+4plJSW8f7UOnEZZhjJqfRaFxrW1XoyTlmkkj99/T4/HP8sVdEIL37vYu4eeBqN8/IzrVq9oCoS5bb/vs30pav466VDc7Z4eCozw7IBb2g9GWHq+iLOGPMET8yZzMjeRzDmzCs4rF2nTKtVLyivjHDjE2+yvLiE+358ds61xUnxqoOsIOtdByFRGtXSR5vsL6VLwoKL78ksqBEOWQrGOBSeqTIcbTafWcTwrZp+3di+uvHjErF0Ua02dgWV/Gn6eMavms8/B5/Fq0Mv48GZX/DgzAlU6a57ohEzG8FypVVGgZ08i4xLcRrzLloSKAKdcq3FYAzfqq3jrvHXV/0KSirKufbx13nu+gv590/P4bL7X2Llxi3Y/a/G+a3+V8P37OCjtXWpqPHIvxd20LsOPJ408+W65Qwb8zijl83mlkOP59Whl/l1t8Q6Nfzk0dcpCId55Cfn0KppbtQCVoQqDTltuULuXIknq9laWc7PJ77F9Z+Ookfz1rzz/SsY3rNfptXKOIvXbuSGJ96kU+vm3H/l2RTkZX/XgZ2rDhpSc0ZvaD31ineWzWPY208wZ+M67j/ubO45ehiNwlnv4dorpi9ZxZ0vjOOI/bpw98WnZVqdlNDQCn/nzpV4cobV27dy0fjneWjmBC7sdRhvnjmSA1q2zbRaGWXc9AXc/87nfP+oPtmfPeY4m82lGW3WTxVCojTKq0wsWA2XDzBZGfNX2BZ4i2rUeO3QYcEiEzbSWa1JDUbUI2IJvNVVwExsMsbYYgS6diY5KMo/Z33MV8XL+Ofgs3nrrJH85stxvL5klr1zgxmislUBc/q7daiEZSauWCSckhqM14k6Nfzn48ns26411w47muUbNvP2lFixHvNzD3ZcCI5tD4bVMmCWpB3cWfi7IeFntJ56zWdrlnDGO08yY8Nq/t+xZ/GHo05v0PUS/vDy+3y1YDm/v/A0Du/ZOdPqJIUSy4Z02XKF3LkST86yrmwbl7z/Io/OmcTI3kfy4ukXsU+jpplWKyNURaL87Km3WblxC/ddeRadWjfPtEpJ0dBcB97QerKCiCp//fojbvzsDfq37cjbZ/2Iw9tl54xub9laVs5Nj79Jfl6Yf111No3ys8sDuLPwd0MytNn1CVkIoTQO79lHm0yGiUvE0+pbTaLwjIuv17amMGz4ek1fK0DU8MnaHsdS5cetMvytYnHyRQIypoAt0WDXMe+snM2iccU8csK5/G/oxfzmq3d5dfGMGokQsBt/rNGIw604jUtd6RQlNdjOZBy287Yv3riJXz4/hoeu/AF3X3Y6t/93TECmxtjmOJaEE/Nkti4M1S/DLCJUG7yP1uOp58zdvI6z332KyetW8PfBZ3LnEUMINcDGW5/NW8q/xn7OsCN6c+WpR2VaHXfUuw48nqygpGIHP/roJf47fwpX9xnEoyedQ5MGWJjmyY+mMHbqPG76/rEM7t090+o44RMWPJ4sIqLK3VPe467J4zily/68MvRSOjbJzuDQ3nD3S+NZtGYD915+Bp3btMi0OglRxK868HiyjWcWTOWqD1+le7NWvHHGSPq0bp9pldJKWUUVtzwxms2lO2jfslmm1XFCVZy2XCHrg2EiSoGtNFItCQa/gi3MnZIYDCe/LRAXSGqwVeYyZPIcgmpVoWAevDkrsAZlHAJmZrjRFugy95nJCBAMGwXzCizHGGWkbB1UP1m3kPPHP8NTJ1/Ay0Mv4dqPX+eL1ctqyLh0BHDrGuDS3ttIOLGVKDCLbtlayhv7zKDWzuOWby5h+N/+S1SVUNj2XUlc4csMfpnJJWBU9Nqb6l0+GObxZCfzS9Zzzrj/srK0hKeH/JDhPftmWqW0EjCmxFZ1dGzdnJGnHslPzzyGoQN6Z0CzmqgPhnk82c2asq38cPyzTF23kvuPP5sr+wzItEoZo2mjAr4/oA8/GTqIw/fvwvyi9dwy4gSO798z06p514HHk+1srSzn8vf/x33Hn8XvjjqVto2a8Pdpn2ZarbQz4uh+HNh5HybMW8a4r+YB0KpZIw7u0Ykvpi3JoGa5NVt1IesNrUvCgolbwRiHhAWLn8n0ydrGCRSVsYxj+mRt45h+3JDFp5YnqfHjhsykBpsf15AxjwGojOxZxiXJwU5NmQoquXHCKDZXDOX6g4+hRWEhv/vyvRpXpg4PdFafrUvXWaeiMg4yZqKBiz4KPzltEMf37cmfX/2QuSvXEQoLB3Xdh+8P6suzH061+oyjxuDWDgu2RIdaoti/P7lM1htaj2d3RFW5c/JYtlbu4Cd9B9M0r5BfTHibiC3qlEPkhUJ0bNWcXz0/lhXFJRTkhenTeR+OOrAbk+Yt45NZizOroNoDf7mMN7SenOee6R+xtbKc2w49iUbhPG7+/E0qzelbDlEVjdKhVTNGDOzP+zMW0q9be3q0bY0Ab0yaTSTqMpevW/yqA48nB3lo9gT+MPl9ztj3IB4+8RwKLC6UXOL2Z8fSsVVzbj3rONo2b8Li1Rt46J2JLFy94TuZwgwVo1F8MMzjyVmenDeZimgVfxo0lIdPPIdrP3iDimhwvXQusHVHOXe/PJ7KSARVCBlhjNOOPJBBB3VnW1kFqzaU8PIn36RROx8MyzpEoDAlCQup6boQIXEygksXhspozRmXGUADyDeOqxRLoCtFATOJ1NxnC3SZgSxbiNKUqXKwc4FECMu5g98ASxASeH7R10RF+cvAYTwy5Adc98nrNYxt0gGygJDxmdoSDQI7LO3ZXcYx/op3DlOuEfp2a88BHdvy9qS5373fp2t7vndUb96b9i2VkQhXnnYUSzdsZvLMmgke1u97NR33xlamIqiWTWS9ofV4asuLC6cB8JeBw/j3iSO47pPXc9ZnO2flOtq3bEZeOISqEokqXdu2ZF3JNiL7I4UAACAASURBVD74ZiGVVRHyw2EG9OoaMLR1haq9NVMu4320ngbJiwunceekdzm1ay/uP2444Rwus/jxnMWg0LZ5EwBaNm1E5zYtqayKUJAX5ti+PVhQVJxWnXxmmMfTQHj+22nfBcj+fsyZOR0HP7RHJ/508VAAXp0wk7kr1vGXkUMZ/+cfUxWJ8N60BQB0aJ2eojSqbluukPWugxBRmoQq9iiTzFKSSJIdFtx8tEbhGcvvXZ7UdF5WWVaYm8kRtkSDVPlxzbErosFxTBl7wsKeo/3WYjVOCQs1MbtGAFBl+L2Bp779iiYFedx26EmUVpXz20nja8i4JCy4FIyx+lbN74YticCwNiGLbzNQL8YyzuTlK7lSI9x+7kk89/k0GhXm0bFNc577eBqPj/8KDUPb5k3553XD+eXjb1NUXGLNjpDq59+LXybvOqhDRGSoiMwXkYUicofl/eEiMkNEpovIFBE5Lp36eRomD82ewKNzJnLZgUfyi8NPyLQ6dcbtz4+lMhLhgmMOpUPr5vzh5Q94+sMpRFUJh4T1JaWM+Wout55zYp3qobgt7colY5y2Ga2IhIGHgNOAImCyiIxW1TnVxD4ARquqisghwMvAQenS0dNwuWf6R7QoaMT1Bx/Dxh1lPDF3cqZVSjlbysr5x9ufAcHlXjuTGJ778GuO69eTa4YdzeOjJ9WNIuq2gieXSKfrYCCwUFUXA4jIS8Bw4DtDq6rbqsk3xbU8qMeTAn4z+V1a5Tfht0cNoXhHKW8umZP4oCwk1l8t9qclEnM99Onanr6d9qFj6+Yc1K09VZEobZo3YePW7XWjRAP7y06noe0CrKj2uggYZAqJyAjgr0B74PvpUc3jiflDb/lsNK0Kf8ithx3P2GXzKY/k3rKvqMZmr6cd1otOrZuzf8e2HNunB2O/mocIXPH//kdxSSmlW8vrTIdccgu4kE5D69K3GVUdBYwSkROAPwKnBgYSuQa4BqBVp0Y0Mp+DEuAU6HJY7mMbJ9hhwVa9K3EwzJTJ1+DK/koj6hG2Vb5yCJiFzGCYJSpjBrZCltJOFVLz6yS2vtcJxjUrgKWS4B938L5XUMlPPnuVJnn5VIYqkXBQRhVaFBSypSJuiCxGQ/MSX0cgqm4zPsZHoQ5JDS7txru0bskvzzmJlz6fzrdrivnXmM/ZvKms5jEOFb6SJZdWFLiQTkNbBHSr9rorsGp3wqr6qYjsLyLtVLXYeO8x4DGArv1aNrCPzFPXbK0sZ2tlcDbXurAxx3balxv6H8vCkg2ERbj2wzcyoOHeU7ShhFN//58a+8KSHgO4s9ZBqhCRocD9xH6WHlfVe4z3WwNPAvsDO4ArVXVW/L2fAVfH1ZoJXKGqO1KmXJx0rjqYDPQSkZ4iUgBcCIyuLiAiB4jEppIicgRQAGwIjOTxpJmezVszsvcRnNxlfx6aMZHrP3qTpvkFXNDrkEyrtldUf3BL2yxTic3eXbYEVAuyDwP6AheJiNnD6NfAdFU9BBhJzCgjIl2Am4ABqtqfmKG+MFWXWZ20zWhVtUpEbgDGEbugJ1V1tohcG3//EeBcYKSIVAJlwAVqe1byeNJIWITLeh8BwAMzJrBk02YA1pRuJRzKbl9jpv66LKU7kiVhkJ2YAf4rgKrOE5EeItIh/l4e0Dhuc5qwh6fsvSGtCQuqOgYYY+x7pNr/7wXurc2YIkphLX20Lt0TbEQwfWHBcQI+WcvfoelbtXdqSJz4EDL+SpL144aitZexJSMEjrEkRwRlMvc7ai9FZNz3cJS7BpxOl6YtuemLN9heVUlhQQEndO5JfjjEuKIFkGexGlXG98ByKwLGxmb1HLrpBsaxqGMeZzN0lxx/OL06t+OuF8fvVqb6vuSf/lO6RtYlyP4NcA7wuYgMBPYFuqrqVBH5B7Cc2MTuPVV9L1WKVcen4Ho8eyAsQkE4zN1T3mN7VSUdGzfnlC77c2ynfZlevJrN5WU5k7rbtFEBIwb1Z3Dv7nV/MnXcoF08eWnndo0xkkuQ/R6gtYhMB24EpgFVcd/tcKAn0BloKiKXpuDqAmR9Cq7HU5dEVCmrquIvg4bxftG3dG3aksJQPvM3F/P8/Gk7V6NmWMvU8OT7U/j+kX341bmncO69zxKp3Pvyo1ZqV72rWFX31Mo4YZBdVbcAVwDEY0BL4tv3gCWquj7+3uvAMcBzrsq54me0Hk8C/jB1PF+sWUr3Zq2ZWrySUYtn89x3RjbG8P36MKTb/pzcdb+M6bm3VEYi3Pv6R/Ro35rLTjqibk/mPqNNhEuQvVX8PYitMPg0bnyXA0eLSJO4AR4CzKUO8DNaj8eBR+fE0lHbFDahuDS2+ickQn4ozO0DTqJFQSHrtm/j+C49yA+FeW/xokyqmzRfzFvGhzMXcs3pA3ln4mzWl5TWzYlS5KN1DLL3AZ4RkQixINlV8fe+FJFXga+JueynEV82mmqy3tCGUBpJ7YJhtic9W9KAiRn8iliSGlwCbWEjymBNajDGtsmYwTCXgJn5GixVt2zncpAJHOOQsFDlIFNXuDy+Vu+40LKgEZf2Opxn5k1jU3kZitK4oIA2jRrzn9lfMXvjOsYuX8BvjjqZT4qWUh6pqjaOhUAwzNbRwNTZMkw4NRW+dn4t//7Wp7z5y5HcOPx4fvv8OEOm2kB7YytT+LE7BNknAr12c+xdwF2p08aOdx14PI6UVOzgjaWz2F5Z8V1zx57N29ChSTNmb1wHgKLMLF5Tw8hmG0UbSnjmk685e1BfDu7RMfUnSOE62mzBG1qPpxasKN1MRJXTux0IwPQNq1izfSvXH3w0txx2LA+eeDazNq7NsJZ7z38++Ir1JaXcNqJuSiY2tMLf3tB6PLWkSqOct9/B/PqIUwD42WfvMKhjd87dvz8PfDORNxbPoXPT5px7QD8GduyaYW2TY3t5JQ+NmcBhPTtz6mHWp+69I3XBsKwg6320gqakqEw+iduxRiRxEkFQxuLvNAqyWBMfHJIjXHy9Lj7asOEsDFkceGYSg5nAYMNWnCYZmVRh+mQ1HFyRH5CxrdoPCzdNGMUTJ17AX44eSufGLYhohJs/f5NpxavYr1UbXh16GR8VLaJfmw78beqnfLh0cfBc5qXbDIvpW7V1j3VIRjD3WT8+Q2bU5NlccuLh3HL2cXw0ZzGVkQikqnttDrkFXMh6Q+vxZIJtVRXcNGEUjcL5dGrUkglrltGhcTOO69SD7s1a8c6Sefx20ngObdeJuwYNYcqqlbuqfWUJUVX+MfpTHr32HC467lCe+eTr1AysYCn+Vu8RkaOBB4itYigg9hNXqqotEh2blOtARM6o9v+zkxnD48l21pZtY9m2TUxYs4zOTVrws0OP54f7H0q7Rk254MBD6NK0BQe13oeirSVURBI/MdVHJsxfxudzl/KT0wfRonFhikZ1DITVv1nvg8BFwLdAY2Jrch9wObDWM1oR+T6xCjkQW+BxMcYCYY+noXHO/v1pll/Izye8TXmkihZ5jbnn2KFs2LGdz1cvzepVCPe99Rmv3HYpV5wygAfe+CI1g2ap/1VVF4pIWFUjwFMiMsHluGRcB+2IFWBoR8yr85ckxvB4coqmeQUs2rKB8kgVbRs1oXlBIf/7dgZvL5lHQSicrXYFgAWrixk7bR6XnHA4L344neItKUhiyM4bsj2eYTZdRP4GrCbWcishyRjaRUAxsVslxHKLZyQxTkoISe0TFsyAlVXG8tiSb7yOWsYxA1K2cwU6LFhlau6rtHY9qOnoqowGP86wIWMG0GxjW5MaklhrY6sCli5sgUoNOyQoBGQs3xXjHmo4xGPzJvH8kItplBemPBIhKlGWbN2AhJVKqmJVtQPjJG43bgaoNBy8pxII8lm6MBhBLKuP1LzUaq8fGjeR0w87kB8PG8hfX/koPohlDFey09BeRuyu3AD8jFiNhXNcDkzGR9saaAvsQ2xWu08SY3g8OcXmijJ+PnE0ZZEq1u/YxrjlC5izaV2m1UoZKzaUMOrL2Zx7zMF0adty7wbL3oSFH6jqDlXdoqq/V9VbgTNdDkzG0L4NtFLV/+7ckhjD48k55m1ez/0zP+O5b7/m45XB5VzZzqPjJxGJKNd8b+BejyVRt62ecbll349cDqy160BVVUSOEpGLgJL4vjEJDvN4GhzHdNyXymiEyeuKMq1KSli3pZRXJ8zgwuMP4/H3vuLbTCuUJuK27mKgp4hUD/w3x7HVVrLraN8nto5sHzLsbREgX2q3dMYtOcGlYExwHNPfavP1mr5dW6JBwG9q+Xk39QmbrU4JdliwdsENJU5qCBRAsTwLBca2+IwJ1z76bnagsPlfo0ZLmXyHx05bPZtwIGEhKBPoVmA5l+Qpfzr6dPJDYYa98zhb1XLdxuBWN3ggqcHie3b4bEy/rS3xIXA/LDPKJz6awnnHHMKPhx3NBw8F33clg+77ZJhALPDVDvhntf1bcYxPJZuC+201t8GCJMfweHKWiCq/nPQOXZq25JeHnZRpdVJG8ZZS/vfFN5x55EF7N1AW+WhVdZmqfqyqg1X1k2rb16q2X9AgyRraEdX+7xMWPB4LU9cX8fT8yYzsPYAB+3TJtDop46kPp1BRtRcJGK51DurZrFdEjhaRySKyTUQqRCQiIltcjk3W0HYQkf1FZD9ivXY8Ho+Ff0z/hJWlJfz1mKHkh3KjhtOGbdv5w8vv79UYWRoMS19mWJzfANfH/393kmOkhBBRGoUq9noc009qrpkFN9+q6f918fVWWBZQmutfk/XjhgwHnumzjSkU3BUcyBw3yemG6bdNwmebLKZvN9+lqIxlHFMmGgqOI/H7vEMruWvKOB4/8Ydc038g/541cdc4pk/d6ls1dtiKepvrZm0yoT2/Bgh8VW2/C/F9b0+b1xDX0SadGZbUT6yqLlfV24HjVHVpMmN4PA2FD1ct5J1lc7n5kOPYt3nrTKtTP8hC1wFGZpiI/AzHzLC9fZbptJfHezwNgt9Pfp+KSIQ/DfpeplXJOKLuWz2jemZYKbHMsHNdDkxoaEXkARG5RkQGi0jzvVLT42mgrCvbxt+mfczxnXpy5r59Mq1O5smiVQc7UdVlxNbOFu7MDFPVhS7HusxoZwKHAPcAS0VkiYiMFpE/x0/q8XgceP7baczcsJo7jzyFJnm2KEADIotcBxLjbhEpBuYBC0RkvYj8znWMhMEwVa3RfldEuhIzvAcTa/GbUcTogmsLGrlgJj1Erb9BiQNd5vnDlnEiZkdZh2QEl4BZhSaObdrqSJmJDrbCM8kEzJwwgmORUHDZUDBhwSJjJCzYkgii4ZrH2f6O88KJxzE7M1i76VoKzyjK3VPf47XTL+f6gwfz968/Mw5J3AUXl4IxDkE1tWRrmONYGxSnaCVAPVxRsCduAY4FjlLVJQDxFVcPi8jPVPW+RAPU+k9DVYtUdYyq3quql9ZaZY+nAfN18UpGLZnJ1X0G0b1Zq0yrkxmyz0c7Erhop5EFUNXFwKXx9xKSGwv7PJ4s4p7pH1EVjfKbAadkWpXMkUWuAyBfVYvNnaq6HvtK0ADe0Ho8aWZd2Tb+PXsCp3c/kEEdumVancyQXYZ2Twv1nRbxZ31zRpGa/lWXAjO2jrJBLAVjDF9hyNaZ1njeiVicWiFS48etdPLJGgkLyXadNVR08uM63GbT/5oqbIVn8nDw4waKaCeWCVkSFkyfsRgJHk99+yWX9DqCOwecwvCxT8c+bYv/NVD425pEkLjwt5nEYBvH3GcJQdT8TPfio6tHbgEXDt1Nqq0AjVwG8DNajycD7IhU8ffpn3BI204M79kv0+qknyya0apqWFVbWLbmqupdBx5PfeaNJbOYsWE1vzjsRArDWf9w6U72BcP2Gm9oPZ4MocCfp35Al6Yt+VHvIzOtTnrJohltKvCG1uPJIF+uW8HHKxdxbb/BNM8vyLQ66aOBGdqsf14RlAKHjgk1D0oc6LIRDGLZvglGUoNl3IgxjrXDrUPAzMQMfIE90aGucOkuHDUjLC4BM7Nali3BwwgIFdiSGgJVt4KfTThaUyG1BLrCIXNhf3Acs2uFNUEgLvPPGZ/w1rAruarfAO6f8YUhY4xt+8qZ57Jcl/nRuATDXGSSQcgtt4ALaZ3RishQEZkvIgtF5A7L+5eIyIz4NkFEDk2nfh5PJpi1aQ1jl8/j6r4DaVXoFMTOfhrYjDZthlZEwsBDwDCgL3CRiPQ1xJYAJ6rqIcAfgcfweBoA9838lKb5BVzb7+hMq1L3aNYW/k6adM5oBwILVXWxqlYALwHDqwuo6gRV3RR/OQnomkb9PJ6M8W1JMaOXzOGy3kfQurBxptWpexrYjDadPtouwIpqr4uAQXuQvwoYm2jQZLrg2gg5+GjDZlcBB1+v6WuFoG/X5lutND4ae+KDmYyQYZe7eam2GUkCmYAPF4iYPlHLZxX0vwY/mwKXhAWHwjNhoyVNxOJ/NZMYbP5gM4nhgVlfcHbPvlzd7yj+Pv2T+MmMsS331KkLg+HHFYs+pozN15uqyoXeR1t32D4i6+0WkZOJGdrbd/P+NSIyRUSmbN6490bW46kPLNqygXeWzeXy3kfSqiDHZ7UNbEabTkNbRKwi+U66AqtMIRE5BHgcGK6qG2wDqepjqjpAVQe0apO+qLrHU9c8MPMLmuUXcmWfozKtSt3hamS9oU2KyUAvEekZ77tzITC6uoCIdAdeBy5T1QVp1M3jqRcsKClmzLJ5/Kj3kTTPL8y0OnVGKjPDHFYztRaRUfHVTF+JSP/4/t4iMr3atkVEbkntlcZIm6FV1SpivXbGAXOBl1V1tohcKyLXxsV+B7QF/h2/8Cnp0s/jqS88PHsiLQoaceEBubu6MVWrDhxXM/0amB5fzTQSuB9AVeer6mGqehhwJLAdGJWyi6xGWqMnqjoGGGPse6Ta/68m1ivdGUHJT8E6ELMyhG0ResThJ9YMqgUCaEDE0NcWMAuM61AprL6RTAKDLRhmtnm3JSyYnRnsAbOaMlWhoH5mha9I1HLfjUCXmcAQO5dZ4cuSsGAE0XYGo2ZtXs3Etcu4ss9RPD13KpXRaDUZhy4MtnOZFb4swwQ+rgQJC3sVGEvdV/e71UwAIrJzNdOcajJ9gb8CqOo8EekhIh1UdW01mSHAonhfsJTjU3A9nnrIY3Mn0alJC87skYONHFPro7WtZupiyHwDnAMgIgOBfQkuHb0QeNH9ImqHN7QeTz3k41WLmL95Pdf029MKyOxEarEB7XauMIpv11iGMzFN9D1AaxGZDtwITAOqvhsgFjM6G3hlry5sD2R9rQOPJ1d5fN6X/P3oMxnccV8mrqmTJ9rM4e46KFbVAXt4P+FqJlXdAlwBsY62xDJQl1QTGQZ8bbgSUkrWG1oB8lPQmtP06dn8n25JDYaPz2Ecmx/XPJWtgE2FsaugbpoVJI2ZaAAEExQMGbPgTmyfi/81ceGZKsMJmWcpGGP6ZE1/LAQ/YzOBITZOzQ8nmS4Mby2fxR2HncKVfY5k4tql8YEs35WIQ6GepIrKBM9lS3RIhhSm1363mglYScwFcHGNc4m0ArbHM1KvBj6NG9+dXEQdug3Auw48nnpLRTTCC99OY0jXXrnXMTdFPlrH1Ux9gNkiMo/Y7PXmnceLSBPgNGLLSusMb2g9nnrMc/O/JqJRRuZSYXDHNbSui2pUdYyqHqiq+6vqn+P7Htm5oklVJ6pqL1U9SFXPqVZPBVXdrqptVbWkLi51J97Qejz1mLVl2xi7bD4/POAQGuc5tafKDnxmmMfjqU88M38qLQoacVYOLfVqaD3Dsj4YFkIpSIFn3anDghlwsRxjBtFsATQzmFJhOXWB8SUzkxxs2Cp8pRMzkJVv6XxhJh+YMrYkB5ekBpcKX1VGtMdMYIBgEoOZwACOAbMUdmGYumEF8zev5+Jeh/G/BTMDMoEOCzYDZbYStwXDjPtaVx0WYgOlaJwswc9oPZ4s4KWF0zisXRf6tm6faVX2HvWFvz0eTz1k1JJZ7Kiq5OIDD8u0KqnB+2g9Hk99o6RiB2NXzOPs/fpSGMru0qA7mzN6H20WEeuwULtjIpYPMOTw82n6X60dbl0qbRgiBWpZOO/g6w2glo9TqoL7UkDU8htt88kGjjOcgxHD6WfrluHio62SmsbH9nnmGc+iZgIDBJMYnIrKWD7zZLowRCzXVV3F15bMZETPgzl93168vXTed/sDiQWWDgvm18fma3XplFtjnPpRVCYr8DNajydLmLh2KStLSzh3v4MzrcpeI6pOW67gDa3HkyUo8PqiWZzQuSftGzfLtDrJk9rqXVmBN7QeTxbx2qJZhEMhfrCfWds6u/CrDjweT71l6dZNTF+/iuE9+2Valb3CB8OyjFj1rl249MR1KUBk+zE1Y70hiw/JJWAWSHywBEECwS9r9MLhajVx23KXql/mdTgFvixBIjPYFUg0sARyQmZg0PIXmJdgXAgmMZgJDODWhcHsdmG7feKgs/mx27owRM0gWkgZvWwOvxtwKvu3as2iLRuDA9kslHk/bI0aTBGbjK16WDLkkBF1wc9oPZ4s461lc4iqcnaPLJ3VprioTDbgDa3Hk2WsLytl4pplnN0ji/20Phjm8XjqO28tm0PPFm2yMiVXAImq05YrZL2PFhHC1ZxJyebMRAx/q22cgFfSoVuRNRHCPM7l+2QLwZo+RhefrQWzU0PY4pwrcPDJmskatuQDs/CM2V3CekygC66tU67hf5XgJ2h+FmYCQ+y4muey+4ONRANbpwaz04bFZ2z66sWSQGGeX+Ov3y9aQCQ6lKHdezO7eH3Ngyz3R83rcElYsI4TPC4Zcskt4IKf0Xo8WciG8u1MWV/E6d0OzLQqtcevo/V4PNnCuBXzOahVe3o0b51pVWqNX0fr8XiygveKFgDwve45PKvNEbyh9XiylJWlJczZuJZTuu6faVVqTUNb3pX1wTAhNb8WIbNClE3IIWBmxjycxrF8oypcStmbz1a2BfjGtMDWqSFsROeilnHMRIew5a4XiNkW3FL5yqhUFkhgsHatqHmMLUBl7jMTGAAqzQpflnsRHCdxwMxMTrCNY0uSMXfZWpKbFb3MgNVHqxfykz6DaVlYyJaKcmA3SQUpSlio8ZEmGxhTcmpFgQt+RuvxZDEfrVpIXijEiZ33y7QqtcO7DjweT7YwfcMqNu7Yzsldssd94At/ezyerCKqyierFnNil/0QsmQSqLqbDpK5S9Yb2liHhdpNzKNOH7LN71Z7P64Vh4QFs7Ovi8/W9MdC0Cdrdte16mPB9J1GLVdvyoQsnSOCPtk9JzBAMNHA9NmCW+GZ/JBRMMZyT4NdGJI7l3lLXfy4Np9oYJ9lnE9XLWHEfv3p26Y9szet3Y1zNYGCtn0OftxkyaXZqgvedeDxZDmfr1kCwHGdemRWkdrgfbQejyebWF9WyrxN6ziuc89Mq+KGgkTUacsVvKH1eHKAz1YvYWD7bhSGs8Qb6Ge0Ho8n25iwZhmF4TyOaNc506o44VcdZCHmgvuEOIjb2nsHg2i2ylxmZ4QkgmOWoW1JDYFKYRYZs5V5heVcZrtsl1biZhALgkkMZgKDbWwzgcGWRBBIWLB1tggkVFiCWMYNs1VWcwl0OXV8MJIPqixty80AmTVgZiQfRC1F1CSkTC1eQSQa5eiO3ZmwqsgiZIxtjbzVfGkLfLnE2ZxoYKsO0jqjFZGhIjJfRBaKyB2W9w8SkYkiUi4it6VTN48nm9laWc6cTWsZ2L57plVxoqHNaNNmaEUkDDwEDAP6AheJiFkifiNwE/CPdOnl8eQKX61bweHtOpNv6XNWnxBteIW/0/mJDAQWqupiVa0AXgKGVxdQ1XWqOhmoTKNeHk9OMHn9Chrl5XNw246ZViUxUcctR0inj7YLsKLa6yJgUDIDicg1wDUA3bqECdX69yLJT9CpM4KDLyyppAabP7jmS2sihiFjSwhIpvCMzQdqJjHYOwDvuUCMTb9AkoMlQSVYMMbmW3UpTlP7Aja2BAqXwjMuPloXmZ3fsa+LY77Zw9t35uvilVaZnQQ6Llhk9nSuvUW8j7bOcAj5uKGqj6nqAFUd0K5tss1rPJ7cYv2OUoq2bebw+r7ywHVpVw7Z4nTOaIuAbtVedwVWpfH8Hk/OM23DKg7fp0um1UhAw6t1kM4Z7WSgl4j0FJEC4EJgdBrP7/HkPNOKV9K1WUvaN26WaVX2iF91UEeoahVwAzAOmAu8rKqzReRaEbkWQEQ6ikgRcCvwGxEpEpEW6dLR48l2vtkQe0is1wGxFKfgOiwbbS0io0Rkhoh8JSL9q73XSkReFZF5IjJXRAan8Eq/I60JC6o6Bhhj7Huk2v/XEHMp1Apbe+w9K5Ls74sR9HDyOqcmqcHW7DtsBrEcnFouFb7MBIbYcUaigW0chy4MZuKFS8ICmtgPn0wXBntyROJ249YW8gmw61PztbWgViDRwBYw2/X/eZvXEVXl4LYd+HDlwu/2Bz5S68kSB3J1N/+vNSlyHVRbNnoaMffkZBEZrapzqon9GpiuqiNE5KC4/JD4e/cD76rqefEn7SYpUcygfi+483g8taIsUsmikg30r88zWkhlMCzhslFi6/Y/AFDVeUAPEekQf1o+AXgi/l6Fqm7euwuz4w2tx5NjzNq4hn5tOmRajT0iqk4b0E5EplTbrjGGsi0bNaOB3wDnAIjIQGBfYk/O+wHrgadEZJqIPC4iTevgcr2h9Xhyjdkb19G5aQvaFDbOtCq7Z2eXhUQbFO9cyhnfHjNGcnHg3QO0FpHpwI3ANKCKmOv0COBhVT0cKAUCPt5UkPVFZWJdcGv5e2HzAxpErJX+zfNYxqmjpIZokoVnooZzzuYrNJMEbH5cM1UvZLl201deabn2QPGXgH/Y1tFgz8dAsNCMLaEieO7E/ld74Zna+4Nd/K+2cYLHJBRh3ua1APRu3Y6Ja5fHD6w/RWVEU1prGLWiJAAACndJREFUNuGyUVXdAlwBICICLIlvTYAiVf0yLvoqdWRo/YzW48kx5m9eD8CBrfbJsCZ7wH1Gm4iEy0bjKwsK4i+vBj5V1S3x4PsKEekdf28IUD2IljKyfkbr8Xhqsr6slE3lZfSu74Y2JcNolYjsXDYaBp7cuWw0/v4jQB/gGRGJEDOkV1Ub4kbg+bghXkx85ptqvKH1eHKQBZvX119Dq6S0YIzDstGJQK/dHDsdGJA6bex414HHk4N8u7mYA1q2zbQau6UWqw5yggY5o3UKnlkCLmaAzD5OapIazAQFq8YO1bvMJAZbRS0XAgEpawcKI6nBcg8rjeQDM6jmEqxz6sLgEFhySTxItsNCMjI2zOCTvaKWuU9YvGUjLQsb07qwMZvKy5y6iri0G08ZOWREXfAzWo8nB1m8ZSMA+7Vok2FNLKhCNOq25Qje0Ho8OcjiLRuAempowRf+9ng82c/K0hIqIhF61lNDm0v+VxdywtDaOtbWFnORvtX/avgBU5XUYBsnUCjH2hnBPHcQc1G+zY/r1GHBoVOuS1JDwL9q+Gyt3WsDSQ3J/ZE6+W2TqM3nUnjGzR9rKxjjcJxlRwRlZWkJ3Zq1Sq2vVXbz/9riDa3H48kFikpL6Nq0ZabVCKJADjVedMH7aD2eHGXFts10bdYq02pYcMwKy6FZr5/Rejw5StG2EvZp3JRG4TzKKutZZCmHVhS44Ge0Hk+OsrJ0CwBdmtazJiU7XQcuW46Q9TNaQQJV8U0ilopQJmZAzVbBKhDockhqcMG2+N9MYrB1WDCv2ibjQiBQk+T3Oxj8CnZGMDszmME6W8Uvt3Ob4zp85knKJBMws+HUSjxwjPv468q2AtC+cTMWilHP2qWVuK2oXEoCawoOf5O5RNYbWo/HY2dt2TYA2jdunmFNLOSQ/9UFb2g9nhxl7faYoe3QpJ51xG2Aqw68ofV4cpTSqgpKKyvoUB9bjzewYFiDMLSmDzcZny3Y/bbpwtbp18Uf7NIp1yw0Y+twa167mZxgPbe1qIzDgeY4KcrFdEtYSOxnTicufts9HVO8o5S2jeqksetekFtLt1xoEIbW42mobCrfTpvCemZoFT+j9Xg8ucPG8jL2aVQnjV33jgY2o/XraD2eHGZT+Xba1DvXAT4zzOPx5A6bystoVVDf2o7nVjKCCw3S0NoSHFwCZCYuFb7Q5LowJJX4YBmnrr7P1sX+WvvAUaL242DpyuDQYcElbpls4kGge4LlZMkkIyRLcOxdX4StleU0zS8gJGKt3JYRFDSSbHpNduJdBx5PDrO1ohyAZvkFCSTTjHcdeDyeXGFbZczQNs8vZEvc6Gacna1sGhDe0Ho8Ocy2ygoAmhcUQmmGlalODs1WXfCGNk4ySQ2e+otLh1sXzKI3qSRVxWn2xPZIzNA2zsuv+Ya1YIzhV67DNrjqZ7QejydX2FFVBUCjcH36U88t/6sL9enuezyeFFMeiRnawvpkaBVoYKsO6tHd93g8qWanoa1PM1oF1K+j9Xg8ucKuGW1mi+PUQBte4e+0rqMVkaEiMl9EForIHZb3RUT+L/7+DBE5Ip36ZRthY8s0YYnW2DypR0QD257YmfgSStCFxIZKcEsVGlWnLVdI24xWRMLAQ8BpQBEwWURGq+qcamLDgF7xbRDwcPxfj8eTBDtXz4RDdbeCoLZsZdO496Mvt3MUL65TZdJEOl0HA4GFqroYQEReAoYD1Q3tcOAZVVVgkoi0EpFOqro6jXp6PDnDzrTbRH310omqDs20DukmnXe/C7Ci2uui+L7ayng8Hkd2zmjrck2sJzHpnNEmbvXqJoOIXANcE39ZHu60cNZe6pZu2pFdj0TZpi94nQFYCgg3p3LI6vSuq4FzjXQa2iKgW7XXXYFVScigqo8BjwGIyBRVHZBaVeuWbNM52/QFr3M6EJEpmdYhW0in62Ay0EtEeopIAXAhMNqQGQ2MjK8+OBoo8f5Zj8eT7aRtRquqVSJyAzCO2GqkJ1V1tohcG3//EWAMcAawENgOXJEu/Twej6euSGvCgqqOIWZMq+97pNr/Fbi+lsM+lgLV0k226Zxt+oLXOR1km74ZQ7SBFXfweDyedFN/Ftd5PB5PjpI1hjbb0ncd9D1IRCaKSLmI3JYJHU0cdL4kfm9niMgEETk0E3oaOiXSeXhc3+kiMkVEjsuEntX02aO+1eSOEpGIiJyXTv12o0uie3ySiJTE7/F0EfldJvSs16hqvd+IBc8WAfsBBcA3QF9D5gxgLLG1uEcDX9ZzfdsDRwF/Bm7Lknt8DNA6/v9hmbzHtdC5GbtcZIcA8+qzvtXkPiQWzzgvC+7xScDbmdSzvm/ZMqP9Ln1XVSuAnem71fkufVdVJwGtRKRTuhWNk1BfVV2nqpOBykwoaMFF5wmquin+chKxdc6ZxEXnbRq3BkBTnPrj1hku32OAG4HXgHXpVG43uOrs2QPZYmizLX23PuniSm11vorYE0QmcdJZREaIyDzgHeDKNOlmI6G+ItIFGAE8Qv3A9XsxWES+EZGxItIvPaplD9liaFOWvpsm6pMurjjrLCInEzO0t9epRolx0llVR6nqQcAPgD/WuVa7x0XffwG3q2p9aUHgovPXwL6qeijwAPBGnWuVZWSLoU1Z+m6aqE+6uOKks4gcAjwODFfVDWnSbXfU6j6r6qfA/iLiWqIv1bjoOwB4SUSWAucB/xaRH6RHPSsJdVbVLaq6Lf7/MUB+Bu9x/STTTmKXjVhixWKgJ7sc8v0Mme9TMxj2VX3Wt5rs3dSPYJjLPe5OLGvvmEzrWwudD2BXMOwIYOXO1/VRX0P+aTIfDHO5xx2r3eOBwPJM3eP6umVFKxvNsvRdF31FpCMwBWgBREXkFmLR3C31VWfgd0BbYrMsgCrNYBEUR53PJVY/oxIoAy7QuEWop/rWKxx1Pg+4TkSqiN3jCzN1j+srPjPM4/F46phs8dF6PB5P1uINrcfj8dQx3tB6PB5PHeMNrcfj8dQx3tB6PB5PHeMNrcfj8dQx3tB6PB5PHeMNrSeliMhPRGRNvMDIIhEZmWmdPJ5M4xMWPClFRB4CZsaz3wYCY1TV5717GjR+RutJNQcD8+P/XwJUZFAXj6de4A2tJ9UcDMyXWDGEG4A7M6yPx5NxvOvAkzJEpBuxWewsYsWhZwCn+gIjnoZOVlTv8mQNhwCfquopItKamMEdLCIrgT8Ra80yClhDrDzkDuAtVX1TRH4EnEys+tNqIB/oD9wFXAB0AB5X1YlpvSKPJwV4Q+tJJQcD0wBUdZOIvECsTnA+8AdV/RZARP4B/FZVl4jIK8Cb8ePHqeoLIvKBqg4RkV8Tq4HaCFgLXAZ4Q+vJOryP1pNKvjO0cd4iViNYgGi1/cKudijV3Qo7a/Guj/9bAfyFWHuXR4EmKdbX40kLfkbrSRmqeonx+lPgcBHpCdwtIquB0cSM5h9FZDvwYoJhPyLWm2xtHajs8aQFHwzzeDyeOsa7Djwej6eO8YbW4/F46hhvaD0ej6eO8YbW4/F46hhvaD0ej6eO8YbW4/F46hhvaD0ej6eO8YbW4/F46hhvaD0ej6eO+f8OBi3qBRPNewAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 360x360 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(figsize=(5, 5))\n",
"\n",
"# False color image of data\n",
"im = ax.imshow(data, \n",
" extent=[X.min(), X.max(), Y.min(), Y.max()],\n",
" origin=\"lower\")\n",
"\n",
"# Contours of two selected values\n",
"cs = ax.contour(X, Y, data, [0.97, 0.98], colors=\"w\")\n",
"# Label the contours\n",
"ax.clabel(cs, inline=True, fmt=\"%.2f\")\n",
"\n",
"ax.set(\n",
" xlabel=r\"$R_\\mathrm{coma}$\",\n",
" ylabel=r\"$L_\\mathrm{crit}$\",\n",
")\n",
"\n",
"# Add a key to the colors\n",
"fig.colorbar(im, ax=ax, shrink=0.8).set_label(\"Data\")\n",
"None"
]
}
],
"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.8.1"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
#!/usr/bin/env python
# coding: utf-8
# In[1]:
import numpy as np
from matplotlib import pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
# Generate grid of parameter values (`R_coma` and `L_crit` in Tere's example)
# In[2]:
X, Y = np.meshgrid(
np.linspace(0.0, 0.5),
np.linspace(0.0, 0.5)
)
# Fake data generated on this grid. Not exactly the same as Tere's, but close enough
# In[3]:
data = 1.0/(1.0 + 0.1*(X**2 + Y**2))
# Now generate the figure with image, labelled contours, and color key.
# In[4]:
fig, ax = plt.subplots(figsize=(5, 5))
# False color image of data
im = ax.imshow(data,
extent=[X.min(), X.max(), Y.min(), Y.max()],
origin="lower")
# Contours of two selected values
cs = ax.contour(X, Y, data, [0.97, 0.98], colors="w")
# Label the contours
ax.clabel(cs, inline=True, fmt="%.2f")
ax.set(
xlabel=r"$R_\mathrm{coma}$",
ylabel=r"$L_\mathrm{crit}$",
)
# Add a key to the colors
fig.colorbar(im, ax=ax, shrink=0.8).set_label("Data")
None
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment