{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# In this notebook we prepare building a SIMPUT of Jupiter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In particular, we will create an extended source representing Jupiter based on a Chandra observation, which sees both emission from the planet as a whole as well as auroral emission."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# imports\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import fits\n",
    "import numpy as np\n",
    "import scipy.ndimage as ndimage\n",
    "\n",
    "from matplotlib.colors import LogNorm\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7efc167d2128>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAADGCAYAAAAzHSAUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztfX3sZkd13nPqLxqnXu9S23LWVteY\nJbACAlsLdksqpXYNmKCYVpCCltpNLVkRpiUliBj5jyitVglSFUiklsSJ09j1ho8aUiwHQhxD1SDZ\njtcOcfE6ZNfONt7Y9RK8XlOqGJuc/vHeu8zOzseZmTNz576/+0g//d733rkzZz7umXOec+59iZmx\nYMGCBQvWG39nagEWLFiwYEF9LMp+wYIFCzYAFmW/YMGCBRsAi7JfsGDBgg2ARdkvWLBgwQbAouwX\nLFiwYAOgirInorcS0TeI6BAR3VijjQULFixYIAdp59kT0WkA/hzAlQCOAHgAwHuY+YBqQwsWLFiw\nQIwalv0bABxi5seZ+bsAPgXg6grtLFiwYMECIU6vUOdWAE8Y348AeKNdiIiuB3A9AJyG0/7hD+Cc\nYKXf23I2TnvmO4pixut2HU+Vwy5fsx8ufHfr2Tjzr+qNG4Cm/ekZsbnNmfua81eC3PvAvK71vbCO\n+DaO/TUznycpW0PZk+PYKVwRM98M4GYAOIe28BvpCgDA8T27AACb9t138gXHPDVrwFe363iqHHb5\nmv1w4cmK7R0b/lv1H9+z69T52wgY5ta1hk+MSepcVJg/7z3mKLdp332n/AeQfx+Y1xnjpbFecuuR\nXqddTgN/yHf8b2nZGjTOEQAXG98vwmrJirBp332n3CTrjlp9PL5n14m6S9pIvXYKRT/KeHjv7ir1\npsBew+MxzTZKELvH7M3A/l9Dntr1hMZY2r52OaDt3NdQ9g8A2E5ElxDRmQDeDeDO3MpcizJlgLRv\nfk1ILaySerVu1F42XXMDMzH2b9tN9ybXF0KLjatWG66+hcZOcz3WXi+x+u3zvXqaLY0AdWXPzC8C\neD+ALwF4FMBnmPkRjbpzFFfqzR+C9uD7+uFrR9K+7UKa1+TKP24aUrlq3uguq1mKkKKrDanidZ0v\nmTfXMV995rnSOaw9rnb9LZS7bWjU3tBCc5UD9dTLHJicvQY0ODOzjim5ulidpW32wq/3IkdLtOiz\nqSxajm/NvqXUXZvHL4FGG3/IdzzIzJdJynbxBO2Y1SGBZKfTcI3MOmpwdVK46gzdwNrWRk59uTy3\nZp2t5C5po3STliLVG9Iah1rjmdKXVO/Zd52LPuzd+7HRhbKPpV9pWiep10snNDc2oK0YUxd3bDxy\nNk5NmiW3zhwZpN4bcOp85xgRJdBUFDUoENcG0zLuI1HmpRumdJxSdEPNMeqaxqntSvnqT6FtgH6D\nPxJM6a6Wjn8Mh/fuVonZrCPFpDX2ofK1Kch1Q854zI7GcSG0GLV2P5/lId31S4KFZnta5XLg66vv\nWE6g0TdGdiBZe/MMKXozEDllRk7OeKbU5YOdgWMfT60n9VwI5tyYVrHWvd9LZpmN2pk5XVn2LQIv\nWtfXRo6FBcQXzNR9rm3Nt4a23HMdBxdK+tJ6HHLa05TRlUUnqXu2lr1G4KXV9bXhsrxClo3Ey/BZ\ncymIpeZpWPmtIRmPWHxCOk+xemty8a2vKQkKt14PvnRl6TUlcCl2jXvVRleWfQw1d1JJGa00rpy0\nzty25gwtzn3d0cuct5KjtR5ogdz4xmwt+xhKHzLQzOpJwSi3i2pJCQRL24pdm2OVt05LPL5nVzNF\nrxkHKoV2dtbUqDGumv3NTc1MgTT+V3Jegg1r2dduy3V9DSse6PtmD2FKqyrHs1sXjP1K7Z/E0+pp\nzHL7OQeMfZqdZS99qKrmgyGlnKHk+lwKyIfSbCATB2/bKW5Xcl6CFK+mJHPJdUzS9ropiBFjv1L7\nJ/G0ajxjESvrWyO5/Qy14fuuCQ1PwIUulH2Nd1prppBpX5+qwGpilGH7NQ+ddFziVraQ32WZxTZA\nybGpMdXct157JcFOn4L1BTNbBEzN77Voz1rrtQtlv84IWZUp1mXpzdkLd5pjwYVuvJT4g6TtVm9J\nlVBIOeek7bbaAFvw6xKUxL1yyvVoYMxK2Ze6UnYAMFa/BrQmvdSS9gVup0xx0ygbS+U0LVlJfSlB\n4dQ1pKVwatKZ9rWtA/MmNIOntdd5D156DN0GaNcxqDIVxsBaaspnzQ2hx/ktkUmrP1OMy5L+eypy\n+lCz32udeqkR6Gm522rRLCFrMbeN0VqVpnyOFp0dANYcT3PT6YU/LrlRfQ9XxeSx+98y28zVpobH\nm0M9lgTgtdsA8rynmnO34VIvteDbJVOP10Ivcsyp3bmnodZGLA1xHazxuUFjzGdn2ae8z16K1IyN\nnOM5kFjwteVoYamVtGe3Ky1bY3ym5mIlMQCJ5+cKzJZ6FLXGRqPeVjG5HDlGlK7X1P50oexTUi9D\nbn/K4p3it2nNndylQEtuuLE/EuXQc1C2NiTUipTuKmlHipRAbYyasz9rpRtrUCvaVFaLFFyJUVBT\noafW3YWyNyHhVXMyS+x6Rx67Voqbq46aC9AMwNZaeDWVtn3j1Eo1jWU1aWZPpaKmB5QSmHdBOmaS\nzdR1bU4wNBc1NuKaSQxaca1ZcfY1eEVpnfaG0jPHGZJNU267rlZj0tpDSe1Xb+Vz60i5N+aQsdXz\nPZuL2XH2Uvgseq06Qzuo7VH4glypyOGHS7JKail6s257HFOtPUm7khzsFlSK3VZqUDlHAdWiOWKQ\neIyacZ3U4HwM66boU9G9ss9xuWPUjC94Uxrky7k2xxVswSOXbiiSfuVY6LGy5vljO0hUp8ZmFMp0\n8dXTetNKRY7XpqH4U9GKXpxj/Sa6V/Y5Cj1WXy53XqIUJBkTmiixlOyyKZkNKf1LpS5i51xxGami\nDrUj8ShdytDHT0tQmw8ekWMVS5MbYvERX/0pkNRfcs9pexfSNVFDT3Sv7H0I0Sgpi8y+VtpmqjJt\nYcH5gsGarrWv/to8umRsXcpWwkub/3Ot15gBMZV1HkIuh53ySona1MlYf2gDKpFBe95qe+8hRJU9\nEf0WER0loq8bx7YQ0d1EdHD4v3k4TkT0q0R0iIgeJqKd/pr1UWINpVr2uTdyTT7d5Mx9x+1zEoUe\noyYkdMQUtESK4k6hX3zrwFdeKkuLdOCY7LXbrWXc+DagUupySuWsDYll/9sA3moduxHAPcy8HcA9\nw3cAuArA9uHvegCf0BGzD9TMBJLQIqUusUspa20wNmyLV5OWqO2626htOIxo8ctctbJmpO3mJjmY\nhkquFxaC5npJMahS6ilFVNkz8/8E8Ix1+GoAtw6fbwXwDuP4bbzCfQDOJaILJYJM7eZKeenQ7p/K\n20kUYKmSdN0kJZDGLWpSGJo3uQQ53Htsk24RZC+5JgWavLPEUJFSc6ntaq2ZmEGVcw+FykmRy9lf\nwMxPAcDw//zh+FYATxjljgzHTgERXU9E+4lo/wt4vokbFLoBpUHbkBKr4fKVur2li9huP7axScex\nRkylFrSzrEo8k1JZctrMRQo11qLdqQ3KESkZZRrlRogeqiKibQDuYuZXD9+fZeZzjfPHmHkzEf0e\ngF9k5q8Ox+8B8GFmfjBUf+sXoU2JUWHmBsd6gJbsLcdgruMtlbu0f3MZH5+cc5FfGy0eqnp6pGeG\n/0eH40cAXGyUuwjAk7HKarwIrQQptEeqtRByRbUsDw1XPxa4lbqiIZj1aNITuemTsfNTWIZSOizH\newhdX4paYyWhO12yuJITNhpylf2dAK4dPl8L4PPG8WuGrJxdAI6PdE8INX6DdkQJb5hLt+QiN7CX\nU49dlzRwWxok1io/XpPLh8c2Wx8dFaLVaiuUlDGa2sptEQiWUpwjlTnFmPS0wUhSLz8J4F4AP0xE\nR4joOgC/BOBKIjoI4MrhOwB8AcDjAA4B+A0A76sidQK0JjjGoZYEJGMKxNWeBkr431rBrNzrNbyw\n1E3CtUHWCB62qKsUNdNGJQHL3jyTFokKqZBk47yHmS9k5jOY+SJmvoWZv8XMVzDz9uH/M0NZZuYb\nmPlSZn4NM++v34UVQm53K6s7VFbqgtvKq5Y1MiVl5Lo+NZPJhkaGjvbN6Nq4tTfrqRXICGna6JQp\nsS4ZXJ9d86RBi6Vm7Wljdk/Q5mSB1Nr1U3jd1Im2LYIUZahB9UiVsNZNUZKBoHXzaClP30Y9NbVS\ngtLNuMR4KR23kHfs+uyi63qy0HPR5SuO5xBZT5Uxd7GH2kmpM2dMJdfMYa6mgHZ2TEl9U2eAzSGD\nZmpZctuf9SuOU7jXlAySlDIS5AQgta2alDpjnoXLI5BY7b3crDZaWl8a4yKNzeRkiml4uGNdqdx8\nSIn1tHbmqOhT0Z2yDykw6eJNSVGTXOND6JqczIwSBVWaJSPZOHq6OWNoyc9rcMehdelSBtIYRYxW\nkWKsK/WVDqU03zpAIzFAA90p+xCkyiiWc+u7ViMIY55LTYHMkVsiS4oMGpjDDa2V5ppSNiezyl5H\nNVJaUw2FkvnVVGy115lWKm3puGklmnSj7FOCnTn1mRZSSUaHljyp7dYMTEplSK1Le05DdcXqzWm3\nNAPLblfTqyyFvYGkKJTahoVPyeZsliXIpV5D9Y0Y6bAW4z2iqwDt1EGSKZHS91DAC0h7qKr2eGsG\nGmsiR67De3dj8wE+5brDe3efQnfE5qb2uPjqd8kaq6PXOTRRkv0jqVsz+F6C2QZoNQJIueenhoaC\njlkiPstIK93QF6g0g3oaizyHXimlwQ7v3X1SH5/74qXYsvPoSefH/z5FX9J+6Rz56pcoes1ArxQa\n/a0lZ2m9WnKljlE3lv0r33tT99ZCDqay2DSu6dWSa+2RjDfV+Ju2ozXvkkPiyfQwnj3IMAXm0O8U\nGWdp2adMgOaj2bUs/hI3sjZ/mpquWZJtohVcagGbXzeV/LEdhM0HGNtuuveE0h/LmFa/hLLSyN4p\nhS1nS0y5FlqNfS9BbBPdWPZTvOK45S7fg0WhLYNWfT2MzYhRlue+eCkAgG8/D4D/Bhx5+2M7KPuX\npkzDoLU3FYr/9OaRlrZle2u98O4lmKVlD7RJpbL5Y1fWTg15SjhZqSfTOp83J13VV08IqVk9oTn1\nYbTOR5xz1WM456rHTshmnx/rHSmdkKJ3xTPM+kx+2fzvk1szPhWK/9RE7RhFrK0c6lNyLLWOlpi9\nZe9znddl57aRyrdr1xu6rtcxk8CVIhnKYDGzcFz8vl1XrC27/ha/R+vCXOZwikyyUJmpvNzZWvYm\npLugPTBaGQMlXHtNpPLtufXmjL9mdk+qHL7yEo/pHz/8NyeO2f1xXb/tpntPOucag1A2yHjOx5tv\nPsCnHA/FPjRjWKkGgutzC2hY5hrXa2co1dzAulP2JYFNbaQukJKgly91UQOSesf2Ncbfd21O/2J1\nxSzp0PWbDzAO792NP3rtS7wyhsbh2A46qXyp0rW5e1f7Lnmm8gJcm1wMNdd5CDmGSE7/NFBrfGZP\n40hQm2poQWVM7V5PFTCs0a5Jw9hKNXcufZtEbPOMUTYuGUJy9U6rrWNQPwdaRu1a0DiaKNmhbUvE\nteum1p/q/oY4YUkbGmjBj7raSR3P0LERJk0y0iku+ifFQ3EFWe3jLmy76V6vR5ATaLZl0oLWegrJ\nlGN116QMa95rsXVRAxvCstdCbaukptcxJT3WoxVmB9ZaemQ2XZMScJ/Kw0o9XwJpgDrFwwH0175G\nUkNJXcf37MIDt39osexroCRdLJaBIanfV2+o3VDAsMQqyuU+tZEih8nzm+Ni15HLvYe8A9/8++I8\ndkzCp+ila07DK4itYV8bqW1LYxChdZXiXY3IHaPYdZI4UM49knpNt5b9lOlnpSixekKcq9SSKWkP\nCN/I2tRAbxZ/DofuKgP4xzFm1YU8P7NeU9baY2nXP+f7Mxc53k7qvKR6nGvB2ZcupFK+s6Rd01rU\n5Flzz0nbC1lA2hyjPU6u89IxlMY9XN/H/6M1b667WHzGxliHOVa2dT32OadPdtumrCmWdkoZs0++\ntkvbyEWre3psK7b+NSx2O+6jFeMAOrbsa2EqazJ3h49dX5Pnt48BfaTEpuLgbTux/ZqHouVSeXHX\nmPjGTvI6hZI5Hi3tKdZ3atxhCoTmtnYsTlI2V4YUy75bZa+9MKQ3ck/UgpQyyAns9HJz9jLervVh\nK/Px+/hSNKncMbpF0qZ5zFfPgmkxxaahSuMQ0cVE9BUiepSIHiGiDwzHtxDR3UR0cPi/eThORPSr\nRHSIiB4mop0iqS2UcN4jXO9Rz3HFxnpru40ulz+GmoEdaQCqBCmWUA5iQUpbudrymFTTSM+MT9AC\np/7ikGudxCz68S2aZhv291aKPmedt6RTeobWnITozRJIOPsXAfwsM78KwC4ANxDRDgA3AriHmbcD\nuGf4DgBXAdg+/F0P4BPqUgdgDrjvJssZSG3O2iWHL0tjamjyhjFo17dp333Ol5hJ2krhwV3Wuo93\nT90MJHXEIC2fs85brNUpN5TStnP1jTaiyp6Zn2Lmh4bP3wbwKICtAK4GcOtQ7FYA7xg+Xw3gNl7h\nPgDnEtGF6pIXYFQAI1ImQ7Osba2lIhTcbIXUOIS0vtBGmNqeaYnbtIgdME3dAEx6BlhZ6S5rfMTB\n23aqjFlJ/GeOmFL+Upahl7FPysYhom0AXg/gfgAXMPNTwGpDAHD+UGwrgCeMy44Mx7qCL5Mhx9rL\nLeviY6W0w+G9u4OZMz0ixSJNCZDG6rBpGNdG4uPrQ8Fqu+5QOuJYzgwUhygTOwsmxzMZkRPTqYWU\nfsyVHupNyY8QK3si+kEAnwXwM8z8XKio49gpUWAiup6I9hPR/hfwvFSMItiLx2XN5d4YqRa/z9qK\nUTnjsR5ynHOoBM0bIDX2YH+35y42H65jEtrQVY8vpdFXV2oGSUmMQyNe5oN0jEPHQ+222CA0DUJX\nfbX6IFL2RHQGVop+HzN/bjj89EjPDP/HX18+AuBi4/KLADxp18nMNzPzZcx82Rk4K1f+k5A6CdKg\nbaxOM5NDek2tIIwJrfo1qASzfAvLMUXxuGi9UKA+dCzUrnnORy1JEVs/uWta6mHmtKkx7744iFQG\n7fZSrpNisgAtERGAWwA8ysy/bJy6E8C1w+drAXzeOH7NkJWzC8Dxke6pjancJlNppLikudxtbCGU\nupElfLnESimZJ+mGnqJ4ju/Z5aT1fHRNrP3UzbHEQq9lFdZQoBr1lMa5WqGU7p0kQAvgTQD+JYDL\niehrw9/bAPwSgCuJ6CCAK4fvAPAFAI8DOATgNwC8T11qRaQumNh7U2ptOCk3n9Ty0LLYQ9dqj0ft\nDV2DmpJsNLkBYddG7EoxjhkdJTRPayVbY51K0dt6K4EkG+erzEzM/Fpmft3w9wVm/hYzX8HM24f/\nzwzlmZlvYOZLmfk1zLy/fjdWMIOXUqTygrH3pkjbacXT+drISbGbE0yl5sqF922e9mfTonfNmWte\nbYXoWo8uD8SniH0xJTt+I4k95VCXknWTu54lHuucFW6sfy3vwW7fjROC76ZwBS9zAqcuSFztEtct\nZdPJXZypXGBqW7241bYSHtdDKgXgWl8uhRpTsqlvcXStDQk96FPCJfMiDQTnrOfQptsSNdp1jUvO\n/Vdy3ka3yl5iJdemCKTWSutFmpOlkHrtWLaUe8yBNMgZkiO0AZsB8tgG71PwdjnXdRK5U4/b53wb\nmNn/0kBwyfla16ZCc/xLoHk/pY5ft+/GqYGYlTIVJNZTb3KPMkllm6IPZpt2+76Xo0nnwvwPhFMj\nU/qeel3NcdWuu/UaSG2vlnw1+70WrziWIJUXrJkFUgKfXLlyt+AgU70rqaVbIpO0/eN7dp3ycFOo\nPh/XP34PKXqzbddT26GAayp9UuoR2ci5X0IB4NbrWDJ+kvI15WiJDWPZayj8ViixBHxBw1YWWmlb\nuVaw69z4WmHzs1Z7pXJqzomvrt48wt7kqY0W/V07y17DCkzln0vbL5FZku0QutYXuJYg1UqWtJ/S\njganOd5kpqLffICDWTEp7bksfJcMZlnpnPjGxZbd5u01LOkYcn+ucUTPir6GN1zbk0mdj64t+yn4\n4F456I34M3A5KLVyfVRNaj12+ZBcqe/H781CbuHN9dZnKWrLvTaWvdRKNMtp8muSjCBNhNqTKnqN\ncXBldaTWGbJ2c66VXpeTBROykmOKPjbeZmzDVXb0QEyvwfyzMb4Az5Y/dE0KWnD7OZ5nKKOoZvxn\nCr6+Vsyta8s+B6kWnBZH3CtyZO65nzHZUmU/vHf3CavatyZSLM9YbMhF60jiAKZnJ/lR9JAMIbli\nx1pDOx6Rc9+nyJVbLnQ94J/LtfhZwlxMuUB7uDlcyFWQGv2ZckzsHxFPCSxLxgxIf3hIY6xLqaV1\nw5z6r3lPHd+zCw/c/qF50zjawU2gPLjkQkkgNLe9Gm62j5pIUT65bZdAQjf5frfAhn3ORbvY1riP\nWvCtC/u47xe0fPWY82KuBSndmQMterQmchIapupPTrKEbz2l3lvdWPavfO9Ns9mdU5FiUU4NLXl6\n6pdr8zKPuWiRVKrE1Ubsx0w0xqdm4D5Vxp7mXIKWFGeq5wbIlPksA7Q5u3NtaLWfYlFODS0rsVa/\nQrL5LGWXLKZlbr9I7PieXdh8wG0EjRa1eb1Z31hHTAlLA72+7+N/sw2p9yqdX2mw21V+6ntXgpw1\nmqqwc+In9loK1Z+CbpS9C1oKI9Wla+Ee94zeNiAToTlx/SBI7IZxBUrtm81VJiQfgOBmEbrORSfF\nrpWu0RzFM16XatH2vIZawJ5P1xylbJ72NTnj27WyHyEdqNiN5EPoBgvRLyH0tEnEuP6ccyX9qzE2\nNmdvZ1z4lGZo8zCvNRV/jrK064yV88kTUv4xOifVUpe2uyAOn3fpgy8GVbKJdsPZa74uoTV3mJvt\nMgWmlkXafomctmdmKmtJ2p2dPTNep9UP12Yk6Y/9ygfJ9VPPt4meZAHqytOqr7Pk7HNgWqxau18O\nUj2HKZEbXCpFifuZKoNtiZrH7JeRuSi8mKIfj0lkcrnxLo9B0h/bcj+2g06p24Rvo5kKvdGiqWtR\n+0eRYvEZbXSr7GPUA3DyTV2iUCVtxa5fZ+RSDyXUR4hak8AXdzF5fbte0+L2yWz36fieXTh4205v\nuZD8ttchWUdmWdePssTazKUlSyA1xErSo320h6S8FNpZT6VrPBXdKnuXdRaCrbBTglZmWykcdcxa\nHWVyvdpWEyUWQiknHxo3jY04F7ZCD23o9jz6jAhfH13vxZdw3HYbknEKWeola6vmHMXurRG5vzBn\ntmF/tlHqYZp1pJ6bGl1y9imWoF3WdsN7olBcmIOMPvT+craUGztkfbp4f9f5VJjjF6oj91xN9LJu\nW8kR0iupMZhQ3allN/TrEuaK1JhDLzdbS9Toc05wXeLR5dy8rnjB2M5GnO8YJGOiYZBMsaFIsZYB\nWs0goaRcbnu5FEoqbdXqxi+NZ9h1lUBKB/jaNv9K6jTH3qboxrdShuSw63JtHuN6CMUNtJESO6jZ\nvu+7Dck9EHqCWYrSey2VEtWs00S3ln2pS7SO7m1O3T1ZhCHKTbsN8//4vvjxPxAOmKb0IVUurXI1\ny9ZEL3L0AI2xWAvLPjf7InZ9bEfMtRpj7eZA2qcQerqxbCs6RVFJz7mCnWMGzvhfouhDQdXxnOQ3\nZSXWvpnNkzMmkjVdcx1I2tcIjJbIIC3TErnrP9cDiyp7InoJEf0xEf0pET1CRL8wHL+EiO4nooNE\n9GkiOnM4ftbw/dBwfluSRJmQZjSUnI9lVYSQQ4eU0BZSaGTupCJ1kYfKp8ynfZOYtM6Ys26e910z\n1mm+EsEnx6Z9951UzlwH42bhyuaxYcvsyxjyySCtPwexdWpvsDXWc04mkwsSKqk13WXLnZK9ZUJi\n2T8P4HJm/hEArwPwViLaBeCjAD7GzNsBHANw3VD+OgDHmPnlAD42lJsE2pNRwp+mcvJm3b4AoUb/\neowTaLZnKlab3jPnxMftuq6xz9mwLf5QJo9ZLmStm3JoGgGaXpZvnYbKad+jWnEh33fzmCaX32Lj\niCp7XuH/Dl/PGP4YwOUA7hiO3wrgHcPnq4fvGM5fQUTfN5uE0HjAopQKkpTP2QBKAjY+xZPSXo13\n+6e0r3VdyOoaYStzU1H61om9MW/ad+qTt3b95nfzF6XMuTLbO75n9WbN8dixHXSSd+FSiBJFlIrU\nOnzB6VjZEq84BS0D2lJI+p4bB0qBiLMnotOI6GsAjgK4G8BjAJ5l5heHIkcAbB0+bwXwBAAM548D\neKmjzuuJaD8R7X8Bz5/Sps/SsheYa3fMGbgYr+pCyHW3ZZNcV9Km9Jrje3aJUtF8sqdmTNjtx9pM\niX+U3iA+y9xF3ZiIvaZghD3O5roYNxRzM4rNizk+pQpMYwNOSWmckqufIm7lov20UYPGATN/j5lf\nB+AiAG8A8CpXseG/y4o/JeWHmW9m5suY+bIzcJZUXu8NZH9ORY2HgzTdvNCxFEhl8rnkEusyJaBq\n11U6ZjEq4dgOOolaMRW+yYePdY3nzPXhWishl9yMDbhkdJW1ZX7ui5eedG3tcQpd1yIYWuL9atUt\npVlc52pQPaVIysZh5mcB/A8AuwCcS0SnD6cuAvDk8PkIgIsBYDi/CcAzoXq/t+XsFDFOQSs3rZVV\nFVIKLa2UFCtbWjbEWafKJIHd3vjj4uM5Hzdufrc9NXt+zA3CZ8kd20Fivn6kjUYv9oWLnwfffl7V\nudeM3cSsWdeG5qPffOdjx0PIMXhy1rQGNOuO5tkT0XkAXmDmZ4no7wL4A6yCrtcC+Cwzf4qIfg3A\nw8z8n4noBgCvYeafJqJ3A/jnzPyToTZqPUErCYqlnqvplk0FzT7NYXykdJ9EaeXUYR9zKbrxyc/n\nvngp+PbzcPQt3xVn7WiOv2QjmwNayV7STs612nn2FwL4ChE9DOABAHcz810Afg7AB4noEFac/C1D\n+VsAvHQ4/kEANyZJrwhf5oLkZo/x8b1Aw6vR7NOU45Pi+kusbJ916bJMXdePdcTiAmMMYGzv8N7d\n2LLzKI7v2YVnHjofm/Z9/2UTlM1gAAAUpklEQVRrsT5KaLWUNVNT0ftokhqeumsuQ+3nypDiAade\nK6kjBEk2zsPM/Hpmfi0zv5qZ//1w/HFmfgMzv5yZ38XMzw/H/2b4/vLh/ONZkgmQs2jHzxK3UGvR\ntaKZStqeUkYbtW40ad2+ALtJ/9gbRojGMpW46/yWnUdPOr/5AJ9Q8vbPG5YG6LXq0ICPJmkVzPXp\nBI3YkQ8a9ebW0e0TtCZ8N2iNrBbXomsVFJUgxm3mtp0io/bGYGdYteRAXTnxqXAFdc0xGp+Q3XyA\ncfC2nSfKjxlgzzx0PoBV4He08mPvql9wMiTGW2xD0aQyQ3JNZVh1+26cEHK49rljXfvVC0JvRxyp\nFvN1xIBbORzeu/uk9/CEvEjznT0vXPw8tl/zUPE816Jb5rb2WsYaxvprtuOrey3ejZOLnhdlScpa\nzRQuaZZRbjZSatraFAil3m7at3pQysXnm2NyeO9uvHDxyc+MmP0brfZjOwhH3/Jd/Id3/Q5e8f4D\nJxS9Wa/reglyPRMNLlnSTiu0oIXs+kNBd602SjBLZd9KoU8xYSFLcIR2uplrsZaUy2k7FyljkTNu\nrqdE7Tk6vHf3Cev8/C+deeLz+FTspn33gd77Tbzx8kdA7/0mNh9g/Oj2Q7jlFZfg6d3PBbNsSgJ+\nUmhy1HOJEdVuv9Q7q4FZ0ji50HKz5ujWtkYPY1RbhoO37cQF5x3HOVc9doK+ecX7D+Dp3c+d4Okv\nOO84AODHtz6Cuz76YyfRQVqYcqx7u6dq1RP7rt2eFLOmcWq6/FMvptJMH80dv2fLJoaUFMuS62Pv\nftl+zUM456rHAKwefNq07z48vfu5E+fG8+dc9Rj+6LUvOUEHxZBK97m8Qddas70RDWjNs3Y92rpC\nI6Mppb0a6E7Z13T5NVCi6EszfUJtS6L+ofTAFEztgqeOg+t6SbmQYravT3lNcQxSuk/63c4UAuSv\nB5l6rntHixRurTnoQtn7XpfQeqGVPEwRgyTVU7v+kAIo6WduELAFpqQzYmjBi8fkSA0Au9ZMiSJr\n5Z3WXgehjCzzuAbdo9WXWXH2LXi92CROBWnfe+DKUyCR11dGmoLbakzmNvZaWJd+zzGle9acfQga\nXHmsHq3MBLvNUo60xsNRUtTMsMihLHzHUx6gqYEeUxRj9F5uW+Z6Tu13Cy+vRsZaLblbeb2zUvYh\nhPhp1/fasBeORgZGifuu6fqXlquFnimcFLTKCy9pq2Q9S9tsTTVqQisW50Ku4TgrGicVNV0vbZqg\nV6qhV/d1QT5K57SEelugi7WlcYA0C7Y0a2OE78EaSTtStKQaXIG2UPk5Y8kmORXaAUPXGNekPXpB\nbv9S7j9NzE7ZS5/yDCHV6hgflQf0f7+15nMFPqQ+CVtCEWlAkt7mQw6fPPUzEFPUX4LamWYuSFKN\na9ef+xxH6v2XWr8Ps1P2PpRyZNIJsl85W4qSoFBtBeB7Ja+NmKVX62G42M2WG6TLuQmlRkYuQvVL\nDZDUQHvKw11Tx8TsY7XWXO1rXUhJMAmhK2Wfsrhs1OYgtdpJhTTjpAZSg3C+tFUNlz4lg6TUctKA\nj9ooqcvXb+k8SY2c2Pi51mDrsY5Z3jU235Z9dPWl9B7qQtmPD1WVpOFpWk0aVorkZs+xinrm0H3K\nRGM8fTGSnAd9NDyPmMelOU8157zEi5mSWkqRe3xPUUkdIdQYh5jnkoO1zcYxrcycB5Kmzibo9eEu\nKULvh29Zhwl7TrXmuKd1sxFRa8ynnsuYDji+ZxceuP1D65uN40KM2rFv8BZZKLalmWrhljzcVUKH\nSa+PIUVJl9IToXpsd7g2Ddaroi8NaJZw85JrtTzz1ASK0Lj0MJexmF4KZmPZj7tsyW7bSy57TTla\n9LGWhbwR4Ru71vOocX9NjRTZU7yxnDFpNY5rmWdfEghqvVNLsldqp1WO0M7UCLXp48J7Thu00UrW\n2Jq010hqJo0ErjhIjfhVK0iycVzjHopB5HojPW6Ys7Hse8JUVlCN9nLrtMegVLa5WpRTeFI9QEOm\ndbh/psZaWfap/HOqJeLj0kMce6qXoWUdlXg1mnWO1x3eu1st9W7cMLTGKjUjJ6XekpeA5UCaXNDS\nCtfot3ZmXa21Lql7Dlg7yz6Xt0u9dkpMyfOWImXMe4mx1Lp27pijdR9qt7esKokMVSx7IjqNiP6E\niO4avl9CRPcT0UEi+jQRnTkcP2v4fmg4v03ahgZSJsguK+XSNfK0Y3WGyoR43tI2tODzllxj7kMv\nMRagnrWpgVSrfkpP04fa42TXH8tjb2VkhL6XZBm5ILbsieiDAC4DcA4zv52IPgPgc8z8KSL6NQB/\nysyfIKL3AXgtM/80Eb0bwD9j5n8Rqnsqzr4HqzK17Snk0pAnx+PqwbqqiXXqX8996Vk2CULyq1v2\nRHQRgB8H8JvDdwJwOYA7hiK3AnjH8Pnq4TuG81cM5btH7RxsKUqtXW0rKZShIR2nHI+rtxtUKyNm\nqjzumtZzb3NlwvTY58i9a42tlMb5OIAPA/jb4ftLATzLzC8O348A2Dp83grgCQAYzh8fyp8EIrqe\niPYT0f4X8Hym+GVIoRS00GKx1d4QfPGCOdxI2m8tdUGSehtDjbFspZC1KK/UVMgQYgZEDWo2BMnm\noy1DVNkT0dsBHGXmB83DjqIsOPf9A8w3M/NlzHzZGThLJOyI3F26ZPB8baZmfdS44XL6pS3Hpn06\nP+eogVB2VunrF1ophRrPYrTajEsyvGL11FpjrQ0/ifeqHYOTWPZvAvATRHQYwKewom8+DuBcIjp9\nKHMRgCeHz0cAXAwAw/lNAJ6RCiRR5CluvhY142sz9EBGaZtSzPlBmFxIH2jRTLdzbWg1x77GhhyC\nxMrseW1NLVtOckcJUtdHVNkz80eY+SJm3gbg3QC+zMx7AHwFwDuHYtcC+Pzw+c7hO4bzX+aE/E4t\nvjaVFy3J1Ze005Lm8NEsvcKlSCQ507X7ZOb+a/LwWnJrvAfGPC6xMlPvy9a5/7ntldJ7ZlJBCHba\ncUuUPFT1cwA+SESHsOLkbxmO3wLgpcPxDwK4MbViDSsidVGaE6X9sMdYZ65lpcHraWwAtSmLmEWe\nSolpyDvOW8yCT1HyWha7xsvmUrOkUhGqe4q4hK/NcSxz74mczbDmuLuwdg9VhWAubN/n3tCzbKmI\npVTa1lEpBVN73NZpbnKxjMG0WKvXJWgi92GkXGhZljXqlaCU2rIRs37M86Vz0oJWi8VrNgJSrFPN\n7BpJeynHc8vllp1izcxe2edOauvA6RSBGS1lFmtTewG3Sk/NCa66yrTekFPqbpFqKsE43lNnp0nv\nyRp0jFl2Cm9o9sreFxSpHcBLVXCacYCWPL3keu0F3CL46mvXRq5y0uSqSxSki9dvObYtlLvdTsoc\ntUycmNoD3FCcvQu1OMcpYgIa7fTAwcY4fa36SupsjRI559JHDcx1rmOy+c6ncPYbXtnb6HFB9ChT\nLWj2VauuWGC5JXqQIQVzkxcok9n83eQWfV+U/YJJ4VrkkoU/R8VQG9rpmqlt1/Z6e8YcPOVZZuNo\n8Fk1+e8FcrgWtyQYqn1ThPjYucx7SmZSreypWFup53MePJvimRCN9djTptaNsp9qYLXc/Bpl1wG+\n/rbIKAopSs3gtUsxtU4ptFFbycTGTzt1NuchybkiNMcl66obZR9CSQdbpfHVKFuKddhYWo5XqsVp\nUiymgpdY4zXmpoZnZP63P4faL1VY5nhqeywlcmlcE0NoHkvmuDtlr+3e18p1nwPmbN3UguazB/Zx\nCVU1JQdvyyBFyvMWof6l9N3My5e2n4JYurZ0s5rigczcNdydstdyr1twwqlomdObgpBS8p2X1leS\nLlgDmmsgVpckdmEiNg+SNiX1SmVIfTYk5tFo0axSPj8nbmBvMFp6pGSD1Wgf6FDZ2wjlXIfQ8sEc\niTUAyINtNeROcb9LrSnt+MuUG6R0bnPqtD/nZDCF6vbVG7pWanVPMScxz0BjY7Sv863DEuU91Xpe\nUi8z0PrBjRapamZ+cEm7uWmXqShtZ0o5pXDNiRZK5NIep7mkYuagdE3GMMvUSwl8bm6q9RWqR7J7\nx/hI7UBPygLItRpcSkXL7dZWDBrt1JLTXosltGRI0adQFzXjYBqY6hkC12dfmVykWPQlqa4SzErZ\n+27Q8bgkQCapZ7w2ZfAlrl8IPbh5vSAU2+jRAkxVpjXTjF1ccy1jQXudtsy0kdCU2sHg0vpKx2dW\nyt4H35v9ci3TWIDGrlvbKtfMPootkJKbu+SNirExq6HUS28W33iGFEcOz6tlUWoGLnPLpmSV1NzI\nW3jHvro02jazk3LRFWffE3c3ytILt6nF39obWW5mgY+eSF3YU8+3LUMPMo3Ika0n+TWgsV5rQyu2\nlYPZcvYtqIySelOsFCCPg/Md1wrUSW6cXG4x1auoZb2nzLEtQ0owrTaksuVSNtooyVbx1VeSFaYp\nS6g+rdhW7TXVlbKvCUnaVqqL7isnkSO1npY8vhanXFPx+BRLjXZTNnK7rNa8lVIvU2xOvdQXo2Sl\n8OkQrbGtvVF3ReP4IAmW1kqFrJGiN/X1oXoBufLo0aVe4EbNNE4f5rpGat5fmpTU8T278MDtH1pe\ncbwO0F50JZumRnkNzFWBpGAj9NEHad97NJimmLfZcvZzgUaWi8T1a+0Wp1JIrWga87uWS14LId46\nJw1X2lZK/amQ1Ht47261LCJJ26Vrr8ba7X2DXjtlr5025YKEF/YtWledWlxiSpmUtMmcTIPQ/5w2\nY997VP61s3xaKRdJMHHbTfeKMoU02m6JnHurV4hoHCI6DODbAL4H4EVmvoyItgD4NIBtAA4D+Elm\nPkZEBOBXALwNwP8D8K+Y+aFQ/RIaZw4UQ49oOQ5abfU0d6H0x57kzEEv8s/p3k5Nh60paypnn2LZ\n/xNmfp1R8Y0A7mHm7QDuGb4DwFUAtg9/1wP4RKzi7205+5RjpRZwD4t4jiixUrSzJ2LoKf1RC5KU\n3NT0UhujV1fDEk9BSkKACa05SO2fS3FrU6MpSB2HEhrnagC3Dp9vBfAO4/htvMJ9AM4logtDFZ32\nzHdOOaZ9k/XiarWWQ/NGypU956aSXDt13rWUSorRWJI8eTsFuOR+kGTl1FinpfnpOTEKzX5obvTH\ndpBaXVJIlT0D+AMiepCIrh+OXcDMTwHA8P/84fhWAE8Y1x4ZjhUj1dUzP+dem4LcoGvtDaDktQYm\nJBtBSk66pB1NKy7XkpTCN7fjcYki7wUtn5OQIicGkvJsjNTIkJyPXet6w2yt9kZIOfsfYuYnieh8\nAHcD+DcA7mTmc40yx5h5MxH9HoBfZOavDsfvAfBhZn7QqvN6rGgeAHg1gK8X92b++PsA/npqISbG\nMgbLGADLGIyIjcM/YObzJBWdLinEzE8O/48S0e8CeAOAp4noQmZ+aqBpjg7FjwC42Lj8IgBPOuq8\nGcDNAEBE+6VBhnXGMg7LGADLGADLGIzQHIcojUNEZxPR3xs/A3gzVlb4nQCuHYpdC+Dzw+c7AVxD\nK+wCcHykexYsWLBgwTSQWPYXAPjdVUYlTgfwO8z8+0T0AIDPENF1AP4SwLuG8l/AKu3yEFaplz+l\nLvWCBQsWLEhCVNkz8+MAfsRx/FsATkmO51UQ4IZEOW5OLL+uWMZhGQNgGQNgGYMRauPQxbtxFixY\nsGBBXazd6xIWLFiwYMGpWJT9ggULFmwATK7sieitRPQNIjpERDfGr5gniOhiIvoKET1KRI8Q0QeG\n41uI6G4iOjj83zwcJyL61WFcHiaindP2QA9EdBoR/QkR3TV8v4SI7h/G4NNEdOZw/Kzh+6Hh/LYp\n5dYCEZ1LRHcQ0Z8N62H3Bl0H/264F75ORJ8kopes+1ogot8ioqNE9HXjWPLcE9G1Q/mDRHStqy0b\nkyp7IjoNwH/C6n06OwC8h4h2TClTRbwI4GeZ+VUAdgG4Yeir2juGZoQPAHjU+P5RAB8bxuAYgOuG\n49cBOMbMLwfwsaHcOuBXAPw+M78Sq+SHR7HB1gERbQXwbwFcxsyvBnAagHdj/dfCbwN4q3Usae6H\nl1D+PIA3YvXM08+PG0QQzDzZH4DdAL5kfP8IgI9MKVPDvn8ewJUAvgHgwuHYhQC+MXz+dQDvMcqf\nKDfnP6wesrsHwOUA7gJAWD0heLq9JgB8CcDu4fPpQzmaug+F/T8HwF/Y/diA62B8rcqWYW7vAvCW\njbAWsHpT8Ndz5x7AewD8unH8pHK+v6lpnGrv0ekZgwv6egD3Y4J3DE2MjwP4MIC/Hb6/FMCzzPzi\n8N3s54kxGM4fH8rPGS8D8E0A/2Wgsn5zeFhxQ60DZv4rAP8Rq2d0nsJqbh/ExloLI1LnPmtNTK3s\nXa9+W+tcUCL6QQCfBfAzzPxcqKjj2KzHhojeDuAon/yepFA/124MsLJKdwL4BDO/HsB38H233YV1\nHAMMtMPVAC4B8EMAzsaKtrCxzmshBl+fs8ZiamUveo/OuoCIzsBK0e9j5s8Nh58e3i2EnHcMzQxv\nAvATtPoxnE9hReV8HKvXYI8P+Jn9PDEGw/lNAJ5pKXAFHAFwhJnvH77fgZXy30jrAAD+KYC/YOZv\nMvMLAD4H4B9hY62FEalzn7Umplb2DwDYPkTgz8QqQHPnxDJVAa3eN3ELgEeZ+ZeNUxvmHUPM/BFm\nvoiZt2E1119m5j0AvgLgnUMxewzGsXnnUH7W1hwz/x8ATxDRDw+HrgBwABtoHQz4SwC7iOgHhntj\nHIcNsxYMpM79lwC8mYg2Dx7Sm4djYXQQrHgbgD8H8BiAm6aWp2I/fxQrV+thAF8b/t6GFe94D4CD\nw/8tQ3nCKlPpMQD/C6ushcn7oTgePwbgruHzywD8MVbvU/pvAM4ajr9k+H5oOP+yqeVW6vvrAOwf\n1sJ/B7B5I64DAL8A4M+werHifwVw1rqvBQCfxCpG8QJWFvp1OXMP4F8PY3EIwE9J2l5el7BgwYIF\nGwBT0zgLFixYsKABFmW/YMGCBRsAi7JfsGDBgg2ARdkvWLBgwQbAouwXLFiwYANgUfYLFixYsAGw\nKPsFCxYs2AD4/7WRdZtcDONDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# load and plot the data\n",
    "imgfile = \"JupiterACIS2011Oct2.fits\"\n",
    "\n",
    "raw_img = fits.open(imgfile)[0].data\n",
    "\n",
    "plt.imshow(raw_img+1, norm=LogNorm())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to separate background, the planet and the auroral region.\n",
    "\n",
    "One way to do it would be by counts, for instance finding the region within which the countrate is greater than some value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Planet by eye, aurorae by convolution & filtering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the planet itself, though, we can probably fit by eye, using an ellipsis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "center = [236, 550]\n",
    "radius = 107 # both in pixels\n",
    "\n",
    "axis_angle = 115 * np.pi/180\n",
    "polar_ratio = 66.854 / 71.492\n",
    "planet_mask = np.zeros_like(raw_img)\n",
    "\n",
    "for ii in range(raw_img.shape[0]):\n",
    "    for jj in range(raw_img.shape[1]):\n",
    "        dx = ii-center[0]\n",
    "        dy = jj-center[1]\n",
    "        # get distance along equatorial and polar axis\n",
    "        d_equat = dx*np.cos(axis_angle) + dy*np.sin(axis_angle)\n",
    "        d_polar = -dx*np.sin(axis_angle) + dy*np.cos(axis_angle)\n",
    "        if (np.sqrt((d_equat)**2 + (d_polar/polar_ratio)**2) < radius):\n",
    "            planet_mask[ii,jj] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(350, 100)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAD8CAYAAAD9uIjPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnX+wXVd13z8rli0BxsiWwDzLHj8U\nu6EkMxUWxgbNZBg5TcHJxGQGp4bWocSJEwpJodPWQDrFbeoZ3GnihJY6I2IHGQLGGCgexpnUIFKm\nBAy2Y2yMySDEI5YtZCQjgkMkI3v1j3vO1Xnn7bPPPj/uffee9/3MaN69++yz97rnPe373Wuvtbe5\nO0IIMSR+YrUNEEKIvtHAJoQYHBrYhBCDQwObEGJwaGATQgwODWxCiMFRO7CZ2c1m9riZfa1QdoaZ\n3WVm38x+np6Vm5m918z2mtkDZnbBJI0XQogQKYrtA8CrS2XvAD7r7ucDn83eA7wGOD/7dzVwYz9m\nCiFEOrUDm7t/HniiVHwZsDt7vRt4baH8Fh/xJWCjmS30ZawQQqSwruV9Z7r7AQB3P2BmL8jKtwCP\nFOrtz8oOlBsws6sZqTpO4qTtz+a0lqYIIYbGD/n+IXd/ftv72w5sVVigLJiz5e67gF0Ap9kZfpFd\n0rMpQoh55TN++3e63N92VfRgPsXMfj6ele8HzinUOxt4rL15QgjRnLYD2x3AG7PXbwQ+VSj/1Wx1\n9GLgB/mUVQghpkXtVNTMPgK8CthsZvuBdwPvAW4zs6uAvwUuz6rfCVwK7AV+BLxpAjYLIUSU2oHN\n3V9fcWmFU8xHeyC9patRQgjRBWUeCCEGhwY2IcTg0MAmhBgcGtiEEINDA5sQYnBoYBNCDA4NbEKI\nwaGBTQgxODSwCSEGhwY2IcTg0MAmhBgcGtiEEINDA5sQYnBoYBNCDA4NbEKIwaGBTQgxODSwCSEG\nhwY2IcTg0MAmhBgcGtiEEINDA5sQYnBoYBNCDA4NbEKIwaGBTQgxOGoPTBZilli3dXH8+vi+pVWz\nQ8w2UmxCiMEhxSZ65fjO7QBsWDp8oqyhssrbWLfn3uo6gTaLaq5Nv2I4SLEJIQaHFJtIJldEuRLK\n3x9d3DSuc/DC9aMXF54FwJbr/2qFkipSVFV5vXW52tu6OG67qABhpOrKik4KTeRIsYkkYoOTELOG\nFJsIUlZj6wqKqazcDl5+1or7z/zKsXHdvF7uO4Pl/rNxecCnti6/N1BepSDH7Zb6jfnsqij3IeYD\nKTYhxOCQYhPLVFX5/YasrOjrGr/Orp37sceAkaqJrUxWKabcf1ZWZU1sLve1rG7Wb/GeVCUmpTaf\nSLGtceQ7E0NEik2s4NCOBTZnr8cKp+jrKqyClgn5u+rUURtVdDzgY0u5R6wNagc2M7sZ+EXgcXf/\nmazsWuA3gO9l1d7l7ndm194JXAU8DfyOu//FBOwWLSg70ddtXeTQjgUANhb+05/66FMrHO/j6eK+\npfH9RyNO+eIAV566lq/HpoexqW3qQBVqV4PcsEmZin4AeHWg/AZ335b9ywe1lwBXAD+d3fO/zOyk\nvowVQogUahWbu3/ezBYT27sMuNXdjwHfNrO9wMuBL7a2UHQi5GR/9JpXjl8XwzLG9ywdHjvyczU2\nfl9UWAnhE0cXN43VXjGQdzy1rVFR43CTFqEaTUkJJK5KwldYyGzRZfHgrWb2gJndbGanZ2VbgEcK\ndfZnZSsws6vN7B4zu+fHHOtghhBCLKft4sGNwO8Bnv38feDXAAvU9VAD7r4L2AVwmp0RrCPaU1QW\nR658BQAbPzgSzlv2jMqP79y+zHc2ETv23HtC7bXoI6bUUpPeQ/68GLEE+/I1KbXZpJVic/eD7v60\nuz8DvJ/RdBNGCu2cQtWzgce6mSiEEM1opdjMbMHdD2Rvfxn4Wvb6DuDDZvYHwFnA+cCXO1spaimv\nXhYDanOlVqaoppaVB9KgYsquuNoaSnNqooDakBp0e3zfUnBlOHZP01AVKbjZICXc4yPAq4DNZrYf\neDfwKjPbxmiauQT8JoC7P2RmtwFfZ+Rvfou7Pz0Z00VOcUqZU34P1btzFMvKr2OLB6FBIjTI1IVb\nNHXa5+VNswjWbV0c57yGBvTgoB1pMzSVTw1jEZMlZVX09YHimyL1rwOu62KUEEJ0QZkHc0hZATy5\n5RTW7VkKXju+c/uJTIGaQNpyWaxPCDv2U4NhQztuxEI7osqpQhmteBY9KKZlyrFiKl/uK/UZi/5Q\nrqgQYnBIsc0h+bf9vg9vA2DrG75Y6adqem5Al7qh3T2qKO7oUc5HjfmkYtdSd/uoaidl4aNqF5Sq\n9ov3xJ5PaiiKSEOKbQ44vnP7sk0ahRBxpNhmmNAZAHDCF7X1DctDFqraaKpc2tp5fN/SCT9ZTbhH\n0QeYEoTb92dougJa5Tus3LW34fOUWusXDWwzyvGd20+EJpT+6Me7ZHR0Rve1F1tsQaEurGPd0uHx\n52gy2PS9j1zd4TBttlxququIFhf6Q1PRGUTTTiG6oYFNCDE4zH31889PszP8Irtktc2YOar8N6FY\ntWls65P33Xaa1HSq1TZ6v63fMDUdTEyez/jt97r7y9reL8UmhBgcUmwzSKpSC90TSsqOJWp32Y4n\nddugOiatxprkkkK3zS1Dzzan7vemxYMTdFVsGthmmGB6FOGdNlYz0Xo1+m4zTYwNHNOadjb5AlrL\naCoqhBAlFMe2yqSkCI2VRj41CsRwNZlKhmLA6u6P7bkWI7QAUGdXavpWG2LxaE1PyYpRF5icWi7a\nIcUmhBgc8rGtMiFllHoOQcoiQN4mpG8zNE1SfUtdUqtSn1OszZy6TIimfs9JLobMM/KxCSFECSm2\nVaLPb+qY6plE8G7dho5t2ghda9pmkz5C5Kd5bf7CgaT7+ljJTP2sxTS7aQVjryYK91hD9BFnVdXm\nak5rYlPlMm2nfHXtVMX71Q3ebQf5WXjus4ymokIIUULhHqtA6vSjas/+fNui4n77KcfKzWIuZNup\nctH21FCX2HOvU4B1U+amz1YBupNFik0IMTik2KZI7pyuOsC4TKoSKaqe4GaOFZsy9qUQYuEUtTv4\nFgJ+U3yHseDZoxWncTX1Z6UGDzcNH4ku8nT8XTTxU64FpNiEEINDq6IzwqRWybqGe/R5JkKRUJpV\n+Vpdm5P2SaWsilYppTafK+QDXKurp11XRTUVnRH6iocq/4dYt3S48lDfdVsXObRjATgRu1W2JTYA\n5VO/PFMiVP/4vqXOeaEhju9bqp1+tdnmqepa5XOo6LtqWl5ur3xPjLU6yLVBU1EhxOCQYpsg08wV\nbOuUzhcyiqouNM0LLQbkR+yFFGEx9zW0sFEXhpJi+6Qc5cUjBKtCblJoG7TbpD+puDBSbEKIwSHF\ntsq0+aZtek+l4ztB5XVVBCmqKtWuFW13sC2mpoo2lw+AbmNDH4shVSlfdfatVaTYhBCDQ+EeU6R4\nZkGVYgiFZ/Sxm0YVValXffkFU8JNUlYP265wluv17e+cNHU2zaLNfaAkeCGEKCEfW89ElUWmXEKr\niGOVElA3sRXKurp1QaZt/HWNVnALn6dqVbTsJ+rLx9d0f7W+VU8famraNg+FWsVmZueY2efM7GEz\ne8jM/k1WfoaZ3WVm38x+np6Vm5m918z2mtkDZnbBpD+EEEIUqfWxmdkCsODu95nZc4F7gdcC/wp4\nwt3fY2bvAE5392vM7FLgt4FLgYuAP3L3i2J9DMnH1nZ//Ta7xrZVBH37Zbq0F1OAdWqzT99jmxSt\ntvdAXFlXUT5XNtTOUJi4j83dD7j7fdnrHwIPA1uAy4DdWbXdjAY7svJbfMSXgI3Z4CiEEFOhkY/N\nzBaBlwJ3A2e6+wEYDX5m9oKs2hbgkcJt+7OyA4UyzOxq4GqADTy7hemzScx/FVMUbfxKwWyAjiqu\nyb2pmRWx60cXN1XGYNX5AvtOzG/KMvs6bhuUlGEQ8dGK5SSviprZqcDHgbe5+9/FqgbKVsx33X2X\nu7/M3V92MutTzRBCiFqSFJuZncxoUPszd/9EVnzQzBYytbYAPJ6V7wfOKdx+NvBYXwbPC7F8yzqq\nlF75esp2QG1zDZtE1Kcq0WD/DRVOm76W7TSSn/YU2UK9DamfY7X9oWuF2oHNzAy4CXjY3f+gcOkO\n4I3Ae7KfnyqUv9XMbmW0ePCDfMq6FqjbTbb4vur+1IT21DSeWCJ3H9PO1M/TZrrbB1XpUsVrR658\nRXRn46aBssvSnRIHvabT49TFi7U4OKYoth3AlcCDZnZ/VvYuRgPabWZ2FfC3wOXZtTsZrYjuBX4E\nvKlXi4UQogalVLWg/E1ZlQI0VVUScV5PIgWrSXvTTgvqnLg/BYXTdLGh6rkP9awDpVQJIUQJKbYp\nkBJcWiyv8sk13Wo6xa62iwtN6jRVQLOomNrQaVulwDmxfS0MzQNdFZsGtgkyqcyAqlXSpnSZVnbp\nd5Lt9UmbBYCqdqBbvF8ffeR0PeBnGmgqKoQQJbS7R0v62jk1dl9orzZgRQxWjL6VQF/5lMnZFT1O\nGevi/kLhNH3sTNuXsl1hX0F5NeljQ+TksqEgxSaEGBzysbUkFIibvy9eL5aF1EeThYVJUufrigbo\nFj5XG9tjyjOWEdEmZKLPEJDGe9NNIQRoln2WTZCPTQghSkixtSQlxzLJD7cKAZa95Uf20E5M7TZV\nsSmryX2lt9XZUNd2H30NGYV7rBIpf8B9bCHU5v7UvtvmcTbNgW0bx9bknkm11+fgnfpl14dN8z5Q\naioqhBAlpNh6oE+FEfqmjWUc9NFnig1t2oDVWQBJmYqmXusjyr9qoalcPxaOktpvH1P2WUCKTQgh\nSmhgE0IMDk1FW1K3Khq61qQOnDgXM7YBYp19dX20ISWBv86OlNXguinxNLIXmm48MI1pXt+LEbOI\npqJCCFFCuaKJpKqUQzsWgHqVlbIQkKrUglkOU/gG75IvOz4bM6LKapVYomOdGlVY1VZKLm7MptDv\nNvXvKJrpUVS+gfxRCGeBDFXdhdBUtGcmGafUZQUMSN6CJ2Vltq6vaW3BE6LNbsJtVoK7bP8zqT3V\nUldnZx1NRYUQooQUWw+kqqK83tHFTeOpWJdk8T6j8vtQDk1UT9esjCaO83l0tq+F7IIYUmxCCFFC\niwcNqY1uj2xJNPaN7VuCkvrq4oiP0SbaPmRP38qgrbKrWnRJbWuWMiFipPjImirQeVGrfSDFJoQY\nHFJsDemympav1m1YOnziW7e03XcbYja0VYBNtv5Juaeuv7q+ym2F2qvrp49A3tSV36rfSdvzRI/v\nW4qGfoRWzNeCMqtCiwcJNN3Hq0l7Te9r2ndqX122LQox6X3mmgywTR3xoS+g2H1VfTb9m5lktsi8\nocUDIYQooaloAl3DBFIi0Ce1eBDrq03bqfmbXabWKaRMy8ahN3mWQ4FoIGt+8lOg39TpdxvHfpeg\n2qY5rUNHik0IMTik2CZIk/SdkNM3ZePE2PUqJunfS/E1dbGlkQoqKa9lfrSK59+Wpp+t7IM8urip\n8RmmdalvbUJihoIUmxBicGhVtCGzEORYXLUr2xLz5xXLYrQN6q2jyu8zqV0oQu209UX1ZVNstXgc\n+lNxvS1Vym2WfW86pWrG6LoLRt/xR6ntlf9Tlf8zhHJgJx07NY1YrNRti6pi0jaUFib6CA/pwix8\n8faBwj2EEKKEFFvPFKcaTSLZu+T6NV3a7/tbvY9FgS799rW4Ms3o/dg0uUib7IYhIMUmhBAlasM9\nzOwc4BbghcAzwC53/yMzuxb4DeB7WdV3ufud2T3vBK4CngZ+x93/YgK2T4W6wEmIKIZMqcWUQ5Nv\n3LY5kFX1Y878NmlJTZRTim+raWhJ+f4U6tKnuuR2RvsNLJSEnk/s76xuUWSoai6F2qmomS0AC+5+\nn5k9F7gXeC3wK8CT7v7fS/VfAnwEeDlwFvAZ4B+5+9NVfczDVLTpymLxnqaO4rrTqSaxslfXRuh6\naHU2tY/iYkSXJPnQfSsyDwqD0qRzWLvQ9PdYrJeS3TJPTHwq6u4H3P2+7PUPgYeBLZFbLgNudfdj\n7v5tYC+jQU4IIaZCIx+bmS0CLwXuzoreamYPmNnNZnZ6VrYFeKRw234CA6GZXW1m95jZPT/mWGPD\nhRCiiuRVUTM7Ffi/wHXu/gkzOxM4BDjwe4ymq79mZu8DvujuH8ruuwm4090/XtX2PExF+yBl1avJ\ndjxtp1WTWq1d1kfidLZNnB20n0p2XSlt0k+XKWDULzvD0+m+mMqqqJmdDHwc+DN3/wSAux9096fd\n/Rng/ZyYbu4HzincfjbwWFsDhRCiKSmLBwbsBp5w97cVyhfc/UD2+u3ARe5+hZn9NPBhTiwefBY4\nf94XD2JpTHXMSnxUHylLk4hZq7Kny7mdfRBbHYb0LIOU+l23LRr3MRA111WxpezusQO4EnjQzO7P\nyt4FvN7MtjGaii4Bvwng7g+Z2W3A1xltrPCW2KAmhBB9o8yDBNps89y1v2IfTfNPi7TNXexKUeHm\nYRfFvMpJ+7iq2uiqgJqEZKT4x7q013SzgnkK/5iGYlvzFP8QyjFc5YTxlNiyuqlom8Go7bSyavrT\nNTVrWWxb/sx6sq/J/WX7msb5jT9Hxf2V0+iaoO062sRNAtHDgeZhQOsLpVQJIQaHpqIdKauKVLnf\nJKG9ru5qTItT7JtEX6H6sTqpfVWVpVxrQ52TPzU0qOravKMkeCGEKCEfWw8sCx4tOIWbOHbH3+BF\nZ3t+rSb5vY/wkarFh3G+Ze4nq0nkbxraUgzpaJr7mlonpb3jO7cH/VJ9bI0UVLk1GyR0/dxrPQle\nik0IMTjkY2tJSnhFyu4ebXfqSLUvNcShi7+q2Ebfq5htd0ips6+Pz9vUvtDfTGgHklC7qeEjQ0Hh\nHqtEF0d1LHyEwHkDXTIHmk6J83bbDCht4sTqqJuGl21Lvac8oLQdlLuGxVQduVeuGzqeb93WxYkf\nTD2vaCoqhBgcUmwTJBaesUKFRb552ziWc4rf9CnKser9LEatxz5PqF6xzjgLIjL9heqpYpUNIVLV\nZTSoO3A6mKhGik0IMTik2HqgaSBpVb28blO1FF2MCKnDBLUQU0NNHeaHdiwA1VudF/utay9Wv0m9\nqkWEcp085Cak7Jr2GbveZHePNgHhaw0pNiHE4FC4xypQF+6Rqtj63vkitV4fIRh9hCn0FbIxCWK/\n49T7i6xW6tdqoXCPGaBLpHzIaZ/qFG87JSnWbTO1SpnGFfsK1m84oMU+X5NBuSl1g2dsmt/HIJOa\nXTCUAa0vNBUVQgwOTUV7pu9v7Kb9hugSoV9uv83n6kMxdVI9gWlv3m4snKPvfMu+FNZaUGfa3UMI\nIUrIx9YzMf9ZV2L+nq591fkJYz622nzGhk7uptfqFEx5F9xlNkdsb9JH6HcTK2vjF4ztjiuWI8Um\nhBgcUmxTIHSAR9XKYsq3ettl/kklqEf9exVH6BXvLR/2ElOi5f3hQnVymnzeOnXWV8BtjNpAbym1\nZKTYhBCDQ6uiU6TvFdO2bTRJAYspkb4VYKrPrkv7s9ZeU59bSt0hoFVRIYQoIR9bD1R9m5a/jYtn\nGYyzCgoHC+d1U7/FQ8qm7TY65etVfrOQfX0pt9iqZIzU1Ko+ksfbKKc2K7kh1oJS6wspNiHE4JCP\nbYJM2gfVtK8ucVWp9sSS2+uUbeharO+6TSDbcuTKVwD12yzF6Mufl+rvHBpdfWwa2NYAsWllH+12\n2WWiKhykT0KLNuP+O9je58BS15YWD5qhqagQYnBo8WAKNFE2oW/mWJ28PBQEXNVnqP9yezHbU7bq\nSWVd4VSu1Cloua/UNKdg2lRg6pyqiPpOYI+psrWg0vpEik0IMTjkY5sgVd/OIb9S111Wm6RSpSid\nafmOQvak7K5bt7VS3T1Nty2qarvNwkcTm0O2rAXkYxNCiBJSbKtMF6WW6rNLUXO5coHCmZs1dnXZ\ngie1flLAcYXCaxoYG1NxqWc0NFkNLrZX/j2thc0kY0w83MPMNgCfB9YzWmy43d3fbWYvAm4FzgDu\nA65096fMbD1wC7AdOAz8c3dfivWxFge2tnFn07YjJS6tbyd6UzuatF0evGKDeNFl0OXwmbUWqtEH\n05iKHgN2uvs/AbYBrzazi4HrgRvc/Xzg+8BVWf2rgO+7+3nADVk9IYSYGrXhHj6SdE9mb0/O/jmw\nE3hDVr4buBa4Ebgsew1wO/A/zcx8Fua8M0TMCV211N9EwaUEpVaFZiyzpbKH8N5oqUqs6b11CwnF\n+1Ke7YbS+2C7xTMSOgQRh9RhXZhNsZ5UXnOSFg/M7CQzux94HLgL+BZwxN3zv/v9wJbs9RbgEYDs\n+g+ATZQws6vN7B4zu+fHHOv2KYQYKBrU2pEUoOvuTwPbzGwj8EngH4eqZT8tcq3Y5i5gF4x8bEnW\nDpCYSknZpaKKZfmbkbCIse+oov5Y2QT6qHLWl8uaKrDQvW2DeIu0CcVIDUGJKbBlfrzSuQXj9yz/\nHYhuNAr3cPcjwF8CFwMbzSwfGM8GHste7wfOAciuPw94og9jhRAihVrFZmbPB37s7kfM7FnAzzFa\nEPgc8DpGK6NvBD6V3XJH9v6L2fU98q9VE1IEqQG1qelOsXbKKVgr7Gqwo0TRnxezLdRu2eZynabp\naKF268pCvrqqVLVyO8EV1aysGEpz8ML1AGzZM3ofOudBdCdlKroA7DazkxgpvNvc/dNm9nXgVjP7\nr8BfAzdl9W8CPmhmexkptSsmYPfgKB4RVzeFSsmBLN7bx+aTqTFhoUGmaU5n2//kXbcxCg6QkbaW\nfe7ClBJGv8/jhdcAh3YscO7HsonNGo9TmzQpq6IPAC8NlO8DXh4oPwpc3ot1QgjRAqVUCSEGh7Yt\nmhG6rPbVtdOk7ap7mk6ZmtRvm7a1wieXcE+M6BZSFaewj32K2XTzO5efBcBztizAjgUANn/hAACn\nPvpUrT9Q9IMUmxBicCgJfgaZRo5o3k/fbfedvJ2yEWeTtprk5tZRXJgJLVw8es0rAXjOY6P/Y5u/\ncIBDmYrrcp7CWqBrrqimojNIH4NCSspO1z3cUsMt6oj1FbO5zRdASnpZMQWqPGCVp6f5FDQP4zh3\naXT9O5efxZlfObbs3uPARq2CTgVNRYUQg0OKbY5pk46VSkoqU6qaqptGloOUQ7FjdfuXFYkpwKrP\nc3Rx07LUp3GdXLEF9q1bt+deDmVH9eXTzfz6mV85NvHTt0Q1UmxCiMEhxTbHpPqUYsomVfVMwrZy\nX7Eo//KGj7G65S2RQlskldmwdHhs88EsZOPcpUWObDklq7Fy6yEKCvTUR59aZue6rYugLYdWDSk2\nISaEBrXVQ4ptDknZrntcd5X+czU5tyC2qroiT7bwmY9k/q08ALa8TVDu/xqrqYCCKirAvP98NRNW\nhmUcLeSEblg6vEKpVdktposGtjkk5pSeZAZD8d5inVgy/PGd28fTxqp7j5emjMv6qxoMd24fDzrl\nveLy9laEVmxdXDl4BhYjnsxizdbtWWLvDRcD8OL/8d0VfR3ftzS2TzvdzhaaigrRAxrUZgtlHgyE\nLqcolVm3dXGFekqZ9kJamEj5em77k5mjPp9aho4EDIWUHCrlZJaJnUiV971h6fA4zzOfij655ZQV\n01ydODUddGCyEEKUkI9tILRValX+sbJvK+XUq1BZaAdfKKixxU1jpfb3Z42Oy9hcaKccNFu2BwoO\n/siZoXl53tc//cS3APjK4YMAfOOBsznv7X8FnFiU2PjBLy4L6aj6rGL2kGITQgwOKbY1TJWaaqPA\nQu8hS1UKtZupqoMXrucfXvgMAM/67kixFX1d5ST0sS+xoMbKYR/FvvM+8vbuvv5GAC665s3L7tm8\n48ThamO/WtlmMTdoYFvDFP/TLsuBDMSexeLSQoNi0SlfnHbm5FPC5zzmnPux7y6zK69fdN4fLS0w\nwKbxYHdqVj8fEE8MlKMJyficAeAnP/pbI1Oy+LNQaIgGs/lHU1EhxOCQYhNAfeR83VF4eXk+LXzi\nNf8AwPqvnsWW60dO+eJJTnnE/oalw3zjt18IwHlv/9LIlrz9rC04odTyhYJHr3klXDjayDFv/9iv\nbwPgxf/xCHBC+RXtPe/to9fl7cTFsJBiE2sWTTmHiwJ0xQrqDllObSMndt7no9e8chwQm4dnFH1s\n5VzNYrt122yn5J6K2URbgwuRgAa0tYUUm5gIoRSvIwWfWTlVSYgiSqkSQogSmoqKiRDypxV9YVqN\nFJNEik0IMTg0sAkhBocGNiHE4NDAJoQYHBrYhBCDQwObEGJwaGATQgwODWxCiMFRO7CZ2QYz+7KZ\nfdXMHjKz/5yVf8DMvm1m92f/tmXlZmbvNbO9ZvaAmV0w6Q8hhBBFUjIPjgE73f1JMzsZ+H9m9ufZ\ntX/v7reX6r8GOD/7dxFwY/ZTCCGmQq1i8xFPZm9Pzv7FMucvA27J7vsSsNHMFrqbKoQQaSTliprZ\nScC9wHnA+9z9bjN7M3Cdmf0n4LPAO9z9GLAFeKRw+/6s7ECpzauBq7O3xz7jt3+t0yeZLpuBQ6tt\nRCLzZCvMl73zZCvMl70/1eXmpIHN3Z8GtpnZRuCTZvYzwDuB7wKnALuAa4D/AlioiUCbu7L7MLN7\numxRMm3myd55shXmy955shXmy14zu6fL/Y1WRd39CPCXwKvd/UA23TwG/Cnw8qzafuCcwm1nA48h\nhBBTImVV9PmZUsPMngX8HPCN3G9mZga8FsinkncAv5qtjl4M/MDdDwSaFkKIiZAyFV0Admd+tp8A\nbnP3T5vZHjN7PqOp5/3Ab2X17wQuBfYCPwLelNDHrsaWry7zZO882QrzZe882QrzZW8nW2dia3Ah\nhOgTZR4IIQaHBjYhxOCY2sBmZieZ2V+b2aez9y8ys7vN7Jtm9lEzOyUrX5+935tdX5yWjTX2zmwK\nmZktmdmDmV33ZGVnmNld2fO9y8xOnwV7K2y91sweLTzbSwv135nZ+jdm9s+mbOtGM7vdzL5hZg+b\n2Stm9blG7J3VZ/tTBZvuN7O/M7O39fZ83X0q/4B/C3wY+HT2/jbgiuz1HwNvzl7/a+CPs9dXAB+d\nlo019n4AeF2g3qXAnzNaRLntbJ0OAAADMElEQVQYuHsVbF0CNpfK/hujoGmAdwDXz4K9FbZeC/y7\nQN2XAF8F1gMvAr4FnDRFW3cDv569PgXYOKvPNWLvTD7bki0nMYqJPbev5zsVxWZmZwO/APxJ9t6A\nnUCeZ7qbUcgIjFKydmevbwcuyepPjbK9NcxqClnxOZaf7yzaG+Iy4FZ3P+bu32a00v7ymnt6wcxO\nA34WuAnA3Z/yURznTD7XiL1VrNqzDXAJ8C13/w49Pd9pTUX/EPgPwDPZ+03AEXfPT2HL066gkJKV\nXf9BVn+alO3NuS6TwTeY2fqsrCqFbJo48H/M7F4bpaoBnOlZ/GD28wVZ+WrbG7IV4K3Zs705n36w\nurZuBb4H/GnmkvgTM3sOs/tcq+yF2Xu2Za4APpK97uX5TnxgM7NfBB539+JBk7G0q6SUrElRYS+M\nUsheDFwInMEohQxW2d6MHe5+AaOdVd5iZj8bqbva9oZsvRH4SWAbo5zi38/qrqat64ALgBvd/aXA\n3zOaGlWx2s+1yt5ZfLZjbORb/yXgY3VVA2WV9k5Dse0AfsnMloBbGU1B/5CRlMwDhItpV+OUrOz6\n84AnpmBnpb1m9iGf4RQyd38s+/k48MnMtoN2IjtkAXg8q76q9oZsdfeD7v60uz8DvJ/ZeLb7gf3u\nfnf2/nZGA8dMPlcq7J3RZ1vkNcB97n4we9/L8534wObu73T3s919kZHk3OPu/wL4HPC6rNobgU9l\nr+/I3pNd3+OZ93AaVNj7L21GU8jM7Dlm9tz8NfDzmW3F51h+vqtib5WtJV/JL7P82V5ho5XyFzHa\n4+/L07DV3b8LPGJm+S4TlwBfZwafa8zeWXy2JV7PiWloblf35zvl1Y9XcWKVcSujB7mXkQxdn5Vv\nyN7vza5vnaaNEXv3AA8y+sP4EHBqVm7A+xitKj0IvGzKNm5ltLr1VeAh4Hez8k2MtpP6ZvbzjNW2\nN2LrBzNbHsj+gBcK9/xuZuvfAK+Z8rPdBtyT2fW/gdNn8bnW2DuTzzbr/9nAYeB5hbJenq9SqoQQ\ng0OZB0KIwaGBTQgxODSwCSEGhwY2IcTg0MAmhBgcGtiEEINDA5sQYnD8f5rBZfa63EzxAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(planet_mask*raw_img)\n",
    "#plt.imshow(planet_mask)\n",
    "plt.xlim(400, 700)\n",
    "plt.ylim(350, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's say this mask gives us roughly the planet. Now, we need to extract the auroral regions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To extract the auroral regions, we could so some form of blurring, follwed by thresholding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# convolution: do a Gaussian filter\n",
    "convimg = ndimage.gaussian_filter(raw_img, 4, mode=\"wrap\")\n",
    "\n",
    "# aurorae have to be inside the planet!\n",
    "convimg *= planet_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7efc168f0ba8>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAADGCAYAAAAzHSAUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADnVJREFUeJzt3XGs3WV9x/H3xxbK1CEUhZS2GRia\nTbPEwhqsY384qoLMWP6ABGJmw5r0H5bhNHGw/WFM9ocmizgTQ2zEWY1TGeJoCJFhwSz7Q6RMhmBF\nrujoXRmVUdDNyGB+98d5Lh7LLffce8/h2Pu8X8nJ+T3P77nnPL/n/Pq5z33OOb+mqpAkrWyvmHYH\nJEmTZ9hLUgcMe0nqgGEvSR0w7CWpA4a9JHVgImGf5OIkDyeZSXLtJJ5DkjS6jPtz9klWAd8H3g7M\nAvcCV1bVd8f6RJKkkU1iZn8+MFNVj1bV/wJfArZP4HkkSSNaPYHHXA8cHCrPAm8+ulGSXcAugFWs\n+r1XcvIEuiJJK9dPOfJkVb1ulLaTCPvMU/eitaKq2g3sBjg5a+vN2TaBrkjSyvX1uvnfR207iWWc\nWWDjUHkDcGgCzyNJGtEkwv5eYFOSs5OcCFwB7J3A80iSRjT2ZZyqej7JnwJ3AKuAz1TVQ+N+HknS\n6CaxZk9V3Q7cPonHliQtnt+glaQOGPaS1AHDXpI6YNhLUgcMe0nqgGEvSR0w7CWpA4a9JHXAsJek\nDhj2ktQBw16SOmDYS1IHDHtJ6oBhL0kdMOwlqQOGvSR1wLCXpA4Y9pLUAcNekjpg2EtSBwx7SeqA\nYS9JHTDsJakDhr0kdWDBsE/ymSSHkzw4VLc2yZ1JHmn3p7b6JPlEkpkkDyQ5b5KdlySNZpSZ/WeB\ni4+quxbYV1WbgH2tDPBOYFO77QJuGE83JUnLsWDYV9U/A08dVb0d2NO29wCXDtV/rga+CZySZN24\nOitJWpqlrtmfUVWPA7T701v9euDgULvZVvciSXYl2Z9k/3M8u8RuSJJGMe43aDNPXc3XsKp2V9WW\nqtpyAmvG3A1J0rClhv0Tc8sz7f5wq58FNg612wAcWnr3JEnjsNSw3wvsaNs7gFuH6t/bPpWzFXhm\nbrlHkjQ9qxdqkOSLwFuB1yaZBT4EfAS4KclO4DHg8tb8duASYAb4GXDVBPosSVqkBcO+qq48xq5t\n87Qt4OrldkqSNF5+g1aSOmDYS1IHDHtJ6oBhL0kdMOwlqQOGvSR1YMGPXkp6sTsO3f/C9kVnbp5i\nT6TRGPbSiIYD/uh6A1+/7gx76SUcK+Cl441hL83DkNdKY9hLQ44V8i7T6Hhn2EsY8lr5DHt1b76g\nN+S10hj26trRQW/Ia6XyS1XqlkGvnhj2Ega9Vj6XcdQtA149cWYvSR0w7CWpA4a9JHXAsJekDhj2\nktQBw16SOrBg2CfZmOTuJAeSPJTkmla/NsmdSR5p96e2+iT5RJKZJA8kOW/SByFJemmjzOyfBz5Q\nVW8AtgJXJ3kjcC2wr6o2AftaGeCdwKZ22wXcMPZeS5IWZcGwr6rHq+pf2/ZPgQPAemA7sKc12wNc\n2ra3A5+rgW8CpyRZN/aeS5JGtqg1+yRnAecC9wBnVNXjMPiFAJzemq0HDg792GyrkyRNychhn+TV\nwFeA91XVT16q6Tx1Nc/j7UqyP8n+53h21G5IkpZgpLBPcgKDoP9CVd3Sqp+YW55p94db/SywcejH\nNwCHjn7MqtpdVVuqassJrFlq/yVJIxjl0zgBbgQOVNXHhnbtBXa07R3ArUP1722fytkKPDO33CNJ\nmo5Rrnp5AfDHwHeSzF0A/C+BjwA3JdkJPAZc3vbdDlwCzAA/A64aa48lSYu2YNhX1b8w/zo8wLZ5\n2hdw9TL7JUkaI79BK0kdMOwlqQOGvSR1wLCXpA4Y9pLUAcNekjpg2EtSBwx7SeqAYS9JHTDsJakD\nhr0kdcCwl6QOGPaS1AHDXpI6YNhLUgcMe0nqgGEvSR0w7CWpA4a9JHXAsJekDhj2ktQBw16SOmDY\nS1IHDHtJ6sCCYZ/kpCTfSvJvSR5K8uFWf3aSe5I8kuTLSU5s9WtaeabtP2uyhyBJWsgoM/tngQur\n6k3AZuDiJFuBjwLXV9Um4Aiws7XfCRypqnOA61s7SdIULRj2NfDfrXhCuxVwIXBzq98DXNq2t7cy\nbf+2JBlbjyVJi7Z6lEZJVgH3AecAnwR+ADxdVc+3JrPA+ra9HjgIUFXPJ3kGOA148qjH3AXsAjiJ\nVy7vKKQV6o5D9/9K+aIzN0+pJzrejRT2VfV/wOYkpwBfBd4wX7N2P98svl5UUbUb2A1wcta+aL/U\ns6NDXlqukcJ+TlU9neQbwFbglCSr2+x+A3CoNZsFNgKzSVYDrwGeGl+XpZVnlHB3Vq/lGOXTOK9r\nM3qS/AbwNuAAcDdwWWu2A7i1be9tZdr+u6rKmbt0DM7i9XIYZWa/DtjT1u1fAdxUVbcl+S7wpSR/\nDXwbuLG1vxH4fJIZBjP6KybQb2lFGDXondVruRYM+6p6ADh3nvpHgfPnqf85cPlYeietMIudxRvy\nGpdFrdlLOrZxL8cY9Bonw15aJkNexwPDXlqiSbyxatBrUgx7aQmczet4Y9hLU2LA6+Vk2EtLcNGZ\nm5c0uzfgNS2GvbREc8E9F/oGuX6dGfbSMhnyOh74P1VJUgcMe0nqgGEvSR0w7CWpA4a9JHXAsJek\nDhj2ktQBw16SOmDYS1IHDHtJ6oBhL0kdMOwlqQOGvSR1wLCXpA4Y9pLUAcNekjowctgnWZXk20lu\na+Wzk9yT5JEkX05yYqtf08ozbf9Zk+m6JGlUi5nZXwMcGCp/FLi+qjYBR4CdrX4ncKSqzgGub+0k\nSVM0Utgn2QD8EfDpVg5wIXBza7IHuLRtb29l2v5trb0kaUpGndl/HPgg8ItWPg14uqqeb+VZYH3b\nXg8cBGj7n2ntf0WSXUn2J9n/HM8usfuSpFEsGPZJ3gUcrqr7hqvnaVoj7PtlRdXuqtpSVVtOYM1I\nnZUkLc3qEdpcALw7ySXAScDJDGb6pyRZ3WbvG4BDrf0ssBGYTbIaeA3w1Nh7Lkka2YIz+6q6rqo2\nVNVZwBXAXVX1HuBu4LLWbAdwa9ve28q0/XdV1Ytm9pKkl89yPmf/F8D7k8wwWJO/sdXfCJzW6t8P\nXLu8LkqSlmuUZZwXVNU3gG+07UeB8+dp83Pg8jH0TZI0Jn6DVpI6YNhLUgcMe0nqgGEvSR0w7CWp\nA4a9JHXAsJekDhj2ktQBw16SOmDYS1IHDHtJ6oBhL0kdMOwlqQOGvSR1wLCXpA4Y9pLUAcNekjpg\n2EtSBwx7SeqAYS9JHTDsJakDhr0kdcCwl6QOjBT2SX6U5DtJ7k+yv9WtTXJnkkfa/amtPkk+kWQm\nyQNJzpvkAUiSFraYmf0fVtXmqtrSytcC+6pqE7CvlQHeCWxqt13ADePqrCRpaZazjLMd2NO29wCX\nDtV/rga+CZySZN0ynkeStEyjhn0B/5TkviS7Wt0ZVfU4QLs/vdWvBw4O/exsq5MkTcnqEdtdUFWH\nkpwO3Jnkey/RNvPU1YsaDX5pzP3iePbrdfODI/ZlJXst8OS0OzFljoFjAI7BnIXG4bdGfaCRwr6q\nDrX7w0m+CpwPPJFkXVU93pZpDrfms8DGoR/fABya5zF3A7sBkuwfei+gW46DYwCOATgGc8Y5Dgsu\n4yR5VZLfnNsG3gE8COwFdrRmO4Bb2/Ze4L3tUzlbgWfmlnskSdMxysz+DOCrSeba/31VfS3JvcBN\nSXYCjwGXt/a3A5cAM8DPgKvG3mtJ0qIsGPZV9Sjwpnnq/wvYNk99AVcvsh+7F9l+pXIcHANwDMAx\nmDO2ccggmyVJK5mXS5CkDhj2ktSBqYd9kouTPNyupXPtwj9xfEqyMcndSQ4keSjJNa2+u2sMJVmV\n5NtJbmvls5Pc08bgy0lObPVrWnmm7T9rmv0elySnJLk5yffa+fCWTs+DP2//Fh5M8sUkJ630cyHJ\nZ5IcTvLgUN2iX/skO1r7R5LsmO+5jjbVsE+yCvgkg+vpvBG4Mskbp9mnCXoe+EBVvQHYClzdjrXH\nawxdAxwYKn8UuL6NwRFgZ6vfCRypqnOA61u7leBvga9V1e8w+PDDATo7D5KsB/4M2FJVvwusAq5g\n5Z8LnwUuPqpuUa99krXAh4A3M/jO04fmfkG8pKqa2g14C3DHUPk64Lpp9ullPPZbgbcDDwPrWt06\n4OG2/SngyqH2L7Q7nm8MvmS3D7gQuI3BN66fBFYffU4AdwBvadurW7tM+xiWefwnAz88+jg6PA/m\nLquytr22twEX9XAuAGcBDy71tQeuBD41VP8r7Y51m/YyTpfX0Wl/gp4L3EN/1xj6OPBB4BetfBrw\ndFU938rDx/nCGLT9z7T2x7PXAz8G/q4tZX26fVmxq/Ogqv4D+BsG39F5nMFrex99nQtzFvvaL+mc\nmHbYj3QdnZUkyauBrwDvq6qfvFTTeeqO67FJ8i7gcFXdN1w9T9MaYd/xajVwHnBDVZ0L/A+//LN9\nPitxDGjLDtuBs4EzgVcxWLY42ko+FxZyrGNe0lhMO+xHuo7OSpHkBAZB/4WquqVVP9GuLcRSrjF0\nnLkAeHeSHwFfYrCU83EGl8Ge+4Lf8HG+MAZt/2uAp17ODk/ALDBbVfe08s0Mwr+n8wDgbcAPq+rH\nVfUccAvw+/R1LsxZ7Gu/pHNi2mF/L7CpvQN/IoM3aPZOuU8TkcH1Jm4EDlTVx4Z2dXONoaq6rqo2\nVNVZDF7ru6rqPcDdwGWt2dFjMDc2l7X2x/Vsrqr+EziY5Ldb1Tbgu3R0HjSPAVuTvLL925gbh27O\nhSGLfe3vAN6R5NT2F9I7Wt1L+zV4s+IS4PvAD4C/mnZ/Jnicf8DgT60HgPvb7RIG6477gEfa/drW\nPgw+qfQD4DsMPrUw9eMY43i8Fbitbb8e+BaD6yn9A7Cm1Z/UyjNt/+un3e8xHftmYH87F/4ROLXH\n8wD4MPA9BhdW/DywZqWfC8AXGbxH8RyDGfrOpbz2wJ+0sZgBrhrlub1cgiR1YNrLOJKkl4FhL0kd\nMOwlqQOGvSR1wLCXpA4Y9pLUAcNekjrw/xgcZtIjpZdMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# now do a simple thresholding - where the countrate exceeds our threshold, we have aurorae\n",
    "aurora_mask = convimg > 0.6\n",
    "\n",
    "plt.imshow(aurora_mask)\n",
    "#plt.xlim(400,700)\n",
    "#plt.ylim(100, 400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A rough estimate of surface brightness is the sum of the pixel values in a given region, divided by the nr of pixel in that region.\n",
    "\n",
    "i.e. the mean of all pixels in that region:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rough surface brightnesses (cts per pixel):\n",
      "Full image: 0.071611166\n",
      "Background only: 0.05515322\n",
      "\n",
      "Full Planet: 0.30665636\n",
      "Auroral Region: 1.2284216\n",
      "Planet without aurorae: 0.24537095\n"
     ]
    }
   ],
   "source": [
    "print(\"Rough surface brightnesses (cts per pixel):\")\n",
    "\n",
    "print(\"Full image:\", np.mean(raw_img))\n",
    "print(\"Background only:\", np.mean(raw_img[planet_mask != True]))\n",
    "print(\"\")\n",
    "print(\"Full Planet:\", np.mean(raw_img[planet_mask == True]))\n",
    "print(\"Auroral Region:\", np.mean(raw_img[aurora_mask == True]))\n",
    "print(\"Planet without aurorae:\", np.mean(raw_img[planet_mask * np.logical_not(aurora_mask) == True]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(350, 100)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAD8CAYAAAD9uIjPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHeZJREFUeJzt3XmQXeV55/Hvo0ZSq1tLa18ZCbAM\nlk0hY5kw5RnHBeNgGE+wq0gZzcSmHGeU8VJje7LYJFNeZiZTIbEhceIikYEgxzGYwc6YYpxKMDgh\niw2WMIuRTNQWkrVvTQtJTUtN9zN/PO/RubR6ud1913N/n6pb955zz7l6+3D43fd9z/uea+6OiEiR\nTKt3AUREKk3BJiKFo2ATkcJRsIlI4SjYRKRwFGwiUjjjBpuZ3WNmR8zsxyXrFpjZI2a2Mz3PT+vN\nzL5kZt1m9qyZXVnNwouIjKScGtu9wLuGrfs08Ki7rwUeTcsA1wNr02MTcGdliikiUr5xg83dHwd6\nhq2+EdiSXm8B3lOy/qsefgB0mdnyShVWRKQcF0xyv6XufhDA3Q+a2ZK0fiWwt2S7fWndweEfYGab\niFodwFvU2ScimSE45u6LJ7v/ZINtNDbCuhHnbLn7ZmAzQJuZt1e4ICLSvPpgz1T2n2xF6XDWxEzP\nR9L6fcCFJdutAg5MvngiIhM32WB7CLglvb4F+HbJ+g+kq6NXAyeyJquISK2M2xQ1s/uAdwCLzGwf\n8Fng94AHzOxDwM+AX0qbfwe4AegG+oAPVqHMIiJjska4bZH62ESkVB9sc/cNk91fFyNFpHAUbCJS\nOAo2ESkcBZuIFI6CTUQKR8EmIoWjYBORwlGwiUjhKNhEpHAUbCJSOAo2ESkcBZuIFI6CTUQKR8Em\nIoWjYBORwlGwiUjhKNhEpHAUbCJSOAo2ESkcBZuIFI6CTUQKR8EmIoWjYBORwlGwiUjhjPtL8CLV\ntHEK+95XsVJI0SjYpOqmEl5T+VwFX+tSsElFnQJsETAAr5yAF4A/qVNZhgefgq51KNhkyv4YWAdc\nfUV6sQA4DbN2wfqniLRrACPV8BR2xWTuXu8y0Gbm7fUuhIwrC4brgDnAGuCy2cDVwIa0wQqgDTgE\n/Ah4BI7cC4eBP6pxeSdLYVd/fbDN3TdMdn/V2GRMpbWc24GTwADQASyZB7wV+DdEsF0O2BKgHwZf\nhnlAZ4TgUE1LPTWlf7NCrjkp2OQ8w5tsd80G5gLTYO5JYBDoBC4ELiZqaYsB60h79EMfkYInoT+t\nvQ34VJXLXmnZsVDANRcFm7xGaaj9ATB/NhFaKdjoSm+2p/XtRNCdBPb3wbS+eL2LuHKwC44CZ4Dp\nNfkLqkMB11wUbDJip/ptpAybkVYMEsHWQYRZ9jgN7CHCbE/a9iSwH3gGeC76104TXW/NTs3U5qBg\na2HDA+0mYBlwBWAz08qsHdlHVLkWkofdWaI6dpToeOsnAnAA6AX2wpFTEWr9RK2tSDaicGtUCrYW\n9f+ILrJOIp/6iArZMsBWE+H0crwx2Jfv19ZPNEEHiSpYP3AcOAA/OxMXCbImZ5Z1C9Pzlmr/UXWg\nJmpjGjfYzOwe4N3AEXd/U1r3OeA/E9/VAL/t7t9J790KfIg49f+ru/9NFcotE7QR+DOgbTURTF1E\np/8A8V9xN9BDpNt04DT4qXjrJBF8p4H2Y7DgWITV9PTebuDHwGPp31pMDAW5ND0uWwS0wZbD1f4r\n60eDgRtLOTW2e4nB418dtv4Od/9C6QozWwfcDLyR+N/mu2b2encfrEBZZRI2Ah8ErgTa3gBcRlTV\nFhMDabOO/v70OAkcgFf6ItROkDcl+4hvq57YhLuG/VtZTa0XeDo9ADhWhT9MZAzjBpu7P25ma8r8\nvBuB+939DPCimXUDVwHfn3QJZdI2Ah8G1gKzVhCBtgJYnZ47yfvIphHDz05Ev/9R4M56FLog1ESt\nr6nctuhjZvasmd1jZvPTupXA3pJt9qV15zGzTWa21cy21n/uQ/Fk/2OtB+YuImpnbeSd+1kH2Emi\nj+wonD0WzcrtKNQqZSPVuwmAjG6yFw/uBP4n4On5i8CvADbCtiPmlrtvBjZDTKmaZDlkmOx/ot8g\nzW5aQQzRgKiGZWF2lJgSsBfYDr4jhp3tBb5V0xK3Bl1Bra1JBZu7n+sGNrOvAA+nxX1EgyeziuiO\nkSorDbTFwMLpRC1tDnHZ8xB5J9mcbKNYf7Y7hpztJL8AIJWncKudSQWbmS1394Np8b3ERTGAh4Cv\nm9ntRIVhLfDklEspY9oIfJ7IsVkzieDqIgbQ9gM9cDr1nZ0App2Bhcci23qJ5mc38A91KHurUbjV\nRjnDPe4D3gEsMrN9wGeBd5jZeqKZuRv4NQB3f97MHiC6aV4FPqorotWzEfgk8MbVRJiVyob598Z9\n0XYTFz8PkIZtAE/VqJzyWpq9UH26bVGT2gjcClwymxjLMUT0nWVjMzJ9cLYnQu0L532KNAKF2/l0\n26IWs5EYU/NWYNkKYgRsF3Fls59oW2b3FmoDpkW3mkKtcal5WnkKtiayEfiPxE1ql60mamoriDA7\nStTWTsBLA2le5hD89zqVVSZG4VZZCrYm8yZgySIi3S4mams9nJsa8MpALPYSlTaRVqRgawJZZ/Mf\nAZ0ziZnqXcQcpkEixVITNN3bkaPAg7UvqkyBZitUjoKtgZVePbuNGI7WOUAk117yatnT8PKxWHUA\n+EatCyoVpYCbOgVbgxo+DScb5TwwBEsOpYUOYCCmQh0g7vOoWQPFoX63yZvKXFGpkpHmFh4nBtge\nAI6cIapne4FD0ad2FIVaEWme6eRoHFuDGetEvo6YBbU0Pc8hutk+UoNySX21Ws1N49haTJohde4G\nt79f3+KINCQFWwMZr9mhWxGLlEd9bA1CfSkyFp0fE6MaW53phJVyaRhI+VRjqyOFmkyGzpvxKdjq\nRCenTIXOn7Ep2OpAJ6VIdSnYRJqUviBHp2ATkcJRsNWYvmWlknQ+jUzBViP6fUmpFp1X51Ow1YBO\nPKk2nWOvpWATkcJRsFWZvkmlVnSu5RRsVaQTTaQ+FGwiBaIv06BgEykYhZuCrWp0ckk9tfr5p2Cr\nglY/qUTqTcFWYQo1kfpTsIlI4SjYKki1NWkkrXw+KthEpHAUbBXSyt+OIo1GwVYBCjVpVK16birY\nRAquFcNt3GAzswvN7HtmtsPMnjezj6f1C8zsETPbmZ7np/VmZl8ys24ze9bMrqz2H1FPrXjSiDS6\ncmpsrwK/7u5vAK4GPmpm64BPA4+6+1rg0bQMcD2wNj02AXdWvNQiImMYN9jc/aC7P5VenwR2ACuB\nG4EtabMtwHvS6xuBr3r4AdBlZssrXnIRKVurtSwm1MdmZmuANwNPAEvd/SBE+AFL0mYrgb0lu+1L\n64Z/1iYz22pmW33i5RYRGVXZwWZms4FvAp9w95fH2nSEdedll7tvdvcN7r5hpB2aQat9C0pza6Xz\ntaxgM7PpRKj9pbt/K60+nDUx0/ORtH4fcGHJ7quAA5UpbuNopZNEpNmUc1XUgLuBHe5+e8lbDwG3\npNe3AN8uWf+BdHX0auBE1mQVEamFcmpsbwPeD1xjZk+nxw3A7wHvNLOdwDvTMsB3gF1AN/AV4COV\nL7aITEartDTMvf5d921m3l7vQkxAq5wcUkz31bsAZeiDbe6+YbL7a+aBiBSOgm2CVFuTZtcK57CC\nTUQKR8E2Aa3wTSdSBAo2ESkcBZuIFI6CTUQKR8FWJvWviTQPBZuIFI6CrQyqrYk0FwWbSAsq+pe1\ngk1ECkfBJiKFo2ATkcJRsIlI4SjYRKRwFGwiUjgKNhEpHAWbiBSOgm0cRR/IKFJECjYRKRwFm4gU\njoJNRApHwSYihaNgE5HCUbCJSOEo2ESkcBRsIlI4CjYRKRwFm4gUjoJNRApHwSYihaNgE2lB99W7\nAFWmYBORwlGwiUjhjBtsZnahmX3PzHaY2fNm9vG0/nNmtt/Mnk6PG0r2udXMus3sBTO7rpp/gIjI\ncBeUsc2rwK+7+1NmNgfYZmaPpPfucPcvlG5sZuuAm4E3AiuA75rZ6919sJIFFxEZzbg1Nnc/6O5P\npdcngR3AyjF2uRG4393PuPuLQDdwVSUKKyJSjgn1sZnZGuDNwBNp1cfM7Fkzu8fM5qd1K4G9Jbvt\nY4QgNLNNZrbVzLb6hIstIjK6soPNzGYD3wQ+4e4vA3cClwDrgYPAF7NNR9j9vOxy983uvsHdN4y0\ng4jIZJUVbGY2nQi1v3T3bwG4+2F3H3T3IeAr5M3NfcCFJbuvAg5UrsgiImMr56qoAXcDO9z99pL1\ny0s2ey/w4/T6IeBmM5tpZhcBa4EnK1fk2ir6QEaRIirnqujbgPcDz5nZ02ndbwMbzWw90czcDfwa\ngLs/b2YPANuJK6of1RVRkcbRCl/W5l7/rvs2M2+vdyHGoN8WlSJphmDrg23uvmGy+2vmgUgLaYZQ\nqwQFm4gUjoJNRApHwSYihaNgE2kRrdK/Bgo2ESkgBVsZWumbToqp1c5hBZuIFI6CTaTgWq22Bgq2\nsrXiySHSrBRsIlI4CjaRAmvVloaCTUQKR8E2Aa367SfNqZXPVwXbBLXyySLSLBRsIlI4CjaRAmr1\nloWCbRJa/aQRaXQKNpGC0Revgm3SdPKINC4Fm0iB6As3KNimQCeRNBKdjzkFm4gUjoJNpABUW3st\nBZtIk1OonU/BNkU6qUQaj4KtAhRuUi8690amYKsQnWAijUPBJtKk9GU6OgVbBelEk1rRuTY2BZuI\nFI6CrcL0TSrVpnNsfAq2KtCJJ9Wic6s8CrYq0QkolXQfOqcmYtxgM7N2M3vSzJ4xs+fN7PNp/UVm\n9oSZ7TSzb5jZjLR+ZlruTu+vqe6fICLyWuXU2M4A17j7FcB64F1mdjVwG3CHu68FXgI+lLb/EPCS\nu78OuCNt15L0DSulPg98EHhfekj1jBtsHk6lxenp4cA1wINp/RbgPen1jWmZ9P61ZmYVK7FIk/nf\nwF0zYeVsWAesBS4EbgU+U8b++oKcuAvK2cjM2oBtwOuALwM/BXrd/dW0yT5gZXq9EtgL4O6vmtkJ\nYCFwbNhnbgI2ARQ59e4DNta7EFIX7wYuA5ZMJ5JsDszvgfnHgT7wIbDpwMDon6FQm5yyLh64+6C7\nrwdWAVcBbxhps/Q8Uk75eSvcN7v7BnffUORgA52creiXgTXAUoDFwAqiqrYBuBK4FGxVvHfX0pE/\nQ+fN5JVVY8u4e6+Z/R1wNdBlZhekWtsq4EDabB/x/bTPzC4A5gE9lStyc1LNrTV8mGieTAfagfnT\niWBbCaxOKw8Bu4j/Y3qB0+d/jkJtasYNNjNbDAykUJsF/DvigsD3gJuA+4FbgG+nXR5Ky99P7z/m\n7ufV2FpRdrIq4Irj54jMuhS4ZDrxFT+XCKwB4mv94rTBpeSJdzS2eWkf/GbJ5ynQKqOcGttyYEvq\nZ5sGPODuD5vZduB+M/tfwI+Au9P2dwN/YWbdRE3t5iqUW6Tufh6YQ1TELukA3kS0VTqJXuaTRIjN\nSc9txBiD08AJ4GhU3DIKtcoZN9jc/VngzSOs30X0tw1f3w/8UkVKJ1Jnvwr0p9ftRIVrgMisHqLZ\nuQaiNnZ5WphJhFkPURVYkHYeBI4Dh4FDcPxMusomFTehPjapDPW3NYffIgIsa6q0pfX96TGdyKzO\neUSgXZyeO4hgO0nU0GamRz/nQo1D0Rp9OH2mamuVpWCrE4Vb4/o40TU2g6hoDZEPH8hCrY1ocS4G\nWJZedKUdu9Ij23ggPU4TNbajQE+8BIVaNSjY6kjh1ng+SWTUHCLYSp0F+oigG0zbLZlJDOWYR96P\n1kYE2zSixtZLhNkA0TzthZeHoptNoVYdCrY6U7g1jk8SLcnODqI5OYNIMOJ5Vj+0n4q8mgOsnEY+\njWAheT9aP3m/Whvwclp3lGiGHo6X6oiuHmuEkRhtZt5e70I0CIVcfdwOzF1BXOJcQF5dy5qSp0te\nQ3SwrSOmFryOqLVl7dZBItD6iasD24GtwHPQebgGf0wB9ME2d98w2f1VY5OW9z5g7jxiDNoKohk5\nnQipLMzmEOGWdbi1E9W71el5ARFmp4mLBkeB/cAO4Acw+CQ8Vas/SBRsjUZN09p6J5FZ5zr8u8hr\nbINEqGXNy+E1tgXpkY1TG0jb9AI/IQapbYVXtsZo9f9Qk79IQMHWkBRulfMZ8lblcaIy1Uv04fcT\nla/pkN+3ppMItyyosmDLQqufvKmZ9aENkg+63ZsePwJ2Ak/Bour+iTICBVuDUrhNzR8Ds2YTQQVw\nGvpPRbBlWdVGZFkb5LWxrJnZQTQ7s+ZnW8mHZ7W2LNCOpv0PEGGWamq+B2ZX6e+TsSnYGlgrzS29\nkRgO1k4MqzgBPE30u0/ETcTNM2YtIiZxziECqBc6duRBNoPIsLmkGltpUzPboOSKKGdLlqcRYXci\nLR8mqoA7geeAXbD/BLx+gmWXylGwNYEi195WEgP2FxOjJhYTuZE1E+cQuTFA9MWX+gNimFjpzTKy\nYWTMSDvPy7fvIpqhfeQtz660GZBf0RzM92GIvEmaVfWGiKDbmwqV1da2x6T2XcDbJ3U0pFI03KPJ\nNHvA/VsiE7KLi+uJueOzsrtfZMMm+oHvw+AzcTVxO9Hia0+bXEcabzYNBk/lXVuD5Lc+s2zoRmd6\nI7tNUNaU7EjvdRCpeiGRsivIm7BZx1waWEsf0fzsAf4RXhqIl4eAX6jgcWp1Gu7RYpq19nYNkS0z\nyLNmHvCWDuAKYgL564iBrnBu2EQbcMUzsWoXeU3uXKWqHdpmwOqeyKceonkJ5LUs0j+4lvxqQSdR\nVetMy3NSwRaSz6c6Sx6E2eyBE+n1XvjhQITp+6d8dKTSFGxNqNn63t5JPml8HtGXtoI0Hek6ItiG\n19Z6ONcGndEP616IzzqaPmcQIqQ6Yx+bB0tOwJJeIuHmEu3abEhGF1Ej60rrs3lTWVOhNz1ns92z\nq6A9RHrtAvYQNbge2D8Ef4amRDUqBVsTa+SA+3mi8pPVoiByZBmRYbNWpxfXE7W11WmHs2mHduK+\n2qeBXujshdWH89sGtUN+r7M5RBgtLnlzTvq8BeRBdiXR1Fwyl9dWD4/DqafyGtnL5INstxNj0rbD\n8VPR5Owhfq1Ioda4FGxSUT9HVKKyPEkVqnOD9y8GZl1EhMyVxP2YL5pBzE1qhxm9MLsbGIoPOJE+\naBksOQntffFZszrIe/5L52VmtbXswsHCtP9S4KoVROS+hZgu0E5U1bph9m7o6YlQO5Aee4G/h9N7\nIt92AY+gQGsGCrYCaJSa27vJ+9CyO/l0ToPTQ9Gqawc6lxIZth54K3DRvwLemPa4gKgTAcv+5bW3\nAVoKDMLcXuLKQxZaWYJm/WVLyQMveywkmrm8D7g2/XtdwCmgO/5N7zl3ZZNd6dEN/5JCrRv4bIWP\nl1SPgq1A6h1wK4mcWUEaR7YMmBPNyM4+ov/qQvL7/18KETJriKGsWfx1wYwO6OrLpzhlVyqzqwbZ\nyP8sRRcRAXcJed9Z6WyCGUtK/tFVwKvAbuB54J/gh8A/E5PVfwLshP2n4Lu89jcJpDko2AqotKlU\nq5D7OHD5bKK/7AqiVpYNSusm75xfTLRH15FqUf1pgxPED5wdgldSFa8jbZ/d/CwLs+z2QIPktbS1\nwPIFRFOT9A/uj88+NgQDR2D514mfxO1Mz9+Ex/viZ4k2w/EDeYXtI5U+QFJTCraCq1XIDUIeWpcT\ng9Oyzvw+4tbYWUdbdlmzD5i7B+gHPxBNwaPEBYTsymTpvMzSYRrZlKdpRFNzIURNbBkRamno7v6h\nmBFwFFj3OKx+PNrLe4ka2t8A34V/OAzvqtKxkdpTsLWQajZV/wS4CyJ0lhLt0sXEVc0O8vFgWdBl\nVz7n7MqnFWTTB7IAzPr2e8jDbiitzwarZVOdTgMzjhG1v/3g2yPQtgM/Jlqd/0xe8zsA/BBe+b7u\nvFFECrYWVLWAy+ZMniQPtL60nK3rI8LoLPnQjNJg20t+caCTqGkdTfu2EX1qJ9P7benz0m90svIA\nzDsQy3uIYHuBvH3Zy7n7rA32lAzklcJRsLWwSjdTf9oHlzxH9J3NIGpH2ZzKHvK7Z2QBl11QyOZa\nvpCeB8mbrNmwi9I+ujXkMwiy6U6DRFJNI4KtdMjGXvjZUGTbv6/A3ymNT8EmQGVCbivQtQMWthOB\nk/Wx7U4bzCDva+sjn6Z0mgihXWkZItjSrbWPnMoH+S54EZZsT58/PX3GYWJGQJqc7qdidVZR7EH9\nZ61Gk+BlQsoJvQ8Soy46iAzrAGatIEZaZGG3H9gDuw/Hy8PktzVrI/JviAionSWfnd3kdmlaziYI\n9ABfm8ofJg1lqpPgFWwyJWMF3fvIf7DpcmBudtEy1eJ+dgL+irylOm2UzxmJRv8Xm4JNGtpGYubU\nOiLcFk6Ds0PROn2GGG0xGoVX61KwiUjhTDXYJlL7FxFpCgo2ESkcBZuIFI6CTUQKR8EmIoWjYBOR\nwlGwiUjhjBtsZtZuZk+a2TNm9ryZfT6tv9fMXjSzp9NjfVpvZvYlM+s2s2fN7Mpq/xEiIqXKmQR/\nBrjG3U+Z2XTgH83sr9N7v+nuDw7b/nrifqZrid/2uDM9i4jUxLg1Ng+n0mJ2F/mxpivcCHw17fcD\noMvMlk+9qCIi5SnrtkVm1gZsI36M8cvu/oSZfRj4XTP7DPAo8Gl3P0PcO3Vvye770rqDwz5zE7Ap\nLZ7pi/ucNotFwLF6F6JMzVRWaK7yNlNZobnKe+lUdi4r2Nx9EFhvZl3AX5nZm4Bbid9KmwFsBj4F\n/A/ARvqIET5zc9oPM9s6lXlhtdZM5W2mskJzlbeZygrNVV4z2zqV/Sd0VdTde4G/A97l7gdTc/MM\n8OfAVWmzfcSPrGVWEbcRFBGpiXKuii5ONTXMbBbx290/yfrNzMyA95A3JR8CPpCujl4NnHD3gyN8\ntIhIVZTTFF0ObEn9bNOAB9z9YTN7zMwWE03Pp4H/krb/DnAD8XNBfcQNVcezecIlr69mKm8zlRWa\nq7zNVFZorvJOqawNcT82EZFK0swDESkcBZuIFE7Ngs3M2szsR2b2cFq+yMyeMLOdZvYNM5uR1s9M\ny93p/TW1KuM45W3YKWRmttvMnkvl2prWLTCzR9LxfcTM5jdCeUcp6+fMbH/Jsb2hZPtbU1lfMLPr\nalzWLjN70Mx+YmY7zOxfN+pxHaO8jXpsLy0p09Nm9rKZfaJix9fda/IA/hvwdeDhtPwAcHN6/afA\nh9PrjwB/ml7fDHyjVmUcp7z3AjeNsN0NwF8TF1GuBp6oQ1l3A4uGrft9YtA0wKeB2xqhvKOU9XPA\nb4yw7TriN19mAhcBPwXaaljWLcCvptcziF/+a8jjOkZ5G/LYDitLGzEmdnWljm9Namxmtor4Ee67\n0rIB1wDZPNMtxJARiClZW9LrB4Fr0/Y1M7y842jUKWSlx3H48W3E8o7kRuB+dz/j7i8SV9qvGmef\nijCzucDbgbsB3P2sxzjOhjyuY5R3NHU7tiO4Fvipu++hQse3Vk3RPwR+i/gNXICFQK+7v5qWs2lX\nUDIlK71/Im1fS8PLm/ndVA2+w8xmpnWjTSGrJQf+1sy2WUxVA1jqafxgel6S1te7vCOVFeBj6dje\nkzU/qG9ZLyZ+i/nPU5fEXWbWSeMe19HKC413bIe7mfzXFityfKsebGb2buCIu28rXT3Cpl7Ge1U3\nSnkhppBdBryV+A3gT2W7jPAxtR5D8zZ3v5K4s8pHzeztY2xb7/KOVNY7iR+PX0/MKf5i2raeZb2A\n+EnUO939zcBpomk0mnof19HK24jH9hyLvvVfBP7PeJuOsG7U8taixvY24BfNbDdwP9EE/UOiKpkN\nEC6ddnVuSlZ6fx7xQ+G1cl55zexr3sBTyNz9QHo+Qvy4+lXAYctnhywHjqTN61rekcrq7ofdfdDd\nh4Cv0BjHdh+wz92fSMsPEsHRkMeVUcrboMe21PXAU+5+OC1X5PhWPdjc/VZ3X+Xua4gq52Pu/p+A\n7wE3pc1uAb6dXj+UlknvP+ap97AWRinvL1uDTiEzs04zm5O9Bn4hla30OA4/vnUp72hlHdZX8l5e\ne2xvtrhSfhFxj78na1FWdz8E7DWz7C4T1wLbacDjOlZ5G/HYDrORvBmalWvqx7fGVz/eQX6V8WLi\nQHYT1dCZaX17Wu5O719cyzKOUd7HgOeIE+NrwOy03oAvE1eVngM21LiMFxNXt54Bngd+J61fSNxO\namd6XlDv8o5R1r9IZXk2ncDLS/b5nVTWF4Dra3xs1wNbU7n+LzC/EY/rOOVtyGOb/v0O4Dgwr2Rd\nRY6vplSJSOFo5oGIFI6CTUQKR8EmIoWjYBORwlGwiUjhKNhEpHAUbCJSOP8fFb4/B+Uq1AIAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# \"Reconstructed\" image: Make the disk have a flat brightness, then add the aurorae on top of it\n",
    "\n",
    "# intensity of the flat disk\n",
    "flat_intens = np.mean(raw_img[(planet_mask * np.logical_not(aurora_mask)) == True])\n",
    "\n",
    "# blur the auroral contribution\n",
    "recons = (planet_mask) * flat_intens + ndimage.gaussian_filter((raw_img-flat_intens), 2)*aurora_mask*planet_mask\n",
    "\n",
    "plt.imshow(recons+1, vmax=5, cmap=\"hot\", norm=LogNorm())\n",
    "plt.xlim(400, 700)\n",
    "plt.ylim(350, 100)\n",
    "#plt.savefig(\"Recons_1.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's try and create images for the simput\n",
    "\n",
    "We'll create two image extensions - one for the aurorae, and one for the planet.\n",
    "\n",
    "Both will be subsets of the original image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x0=100; x1=350\n",
    "y0=400; y1=700\n",
    "disk_img = (planet_mask * flat_intens)[x0:x1, y0:y1]\n",
    "aurora_img = (ndimage.gaussian_filter((raw_img-flat_intens), 2)*aurora_mask)[x0:x1, y0:y1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# need WCS keys from original image\n",
    "orig_header = fits.open(imgfile)[0].header"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Disk image\n",
    "disk_hdu = fits.ImageHDU(data=disk_img)\n",
    "disk_hdu.name = \"DISK\"\n",
    "\n",
    "# copy the header, but adjust the crpix's!\n",
    "disk_hdu.header[\"RADESYS\"] = orig_header[\"RADESYS\"]\n",
    "disk_hdu.header[\"EQUINOX\"] = orig_header[\"EQUINOX\"]\n",
    "disk_hdu.header[\"CTYPE1\"] = orig_header[\"CTYPE1\"]\n",
    "disk_hdu.header[\"CTYPE2\"] = orig_header[\"CTYPE2\"]\n",
    "disk_hdu.header[\"CRVAL1\"] = orig_header[\"CRVAL1\"]\n",
    "disk_hdu.header[\"CRVAL2\"] = orig_header[\"CRVAL2\"]\n",
    "disk_hdu.header[\"CUNIT1\"] = orig_header[\"CUNIT1\"]\n",
    "disk_hdu.header[\"CUNIT2\"] = orig_header[\"CUNIT2\"]\n",
    "disk_hdu.header[\"CRPIX1\"] = orig_header[\"CRPIX1\"] - y0\n",
    "disk_hdu.header[\"CRPIX2\"] = orig_header[\"CRPIX2\"] - x0\n",
    "disk_hdu.header[\"CDELT1\"] = orig_header[\"CD1_1\"]\n",
    "disk_hdu.header[\"CDELT2\"] = orig_header[\"CD2_2\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Aurora image\n",
    "aurora_hdu = fits.ImageHDU(data=aurora_img)\n",
    "aurora_hdu.name = \"AURORAE\"\n",
    "\n",
    "# copy the header, but adjust the crpix's!\n",
    "aurora_hdu.header[\"RADESYS\"] = orig_header[\"RADESYS\"]\n",
    "aurora_hdu.header[\"EQUINOX\"] = orig_header[\"EQUINOX\"]\n",
    "aurora_hdu.header[\"CTYPE1\"] = orig_header[\"CTYPE1\"]\n",
    "aurora_hdu.header[\"CTYPE2\"] = orig_header[\"CTYPE2\"]\n",
    "aurora_hdu.header[\"CRVAL1\"] = orig_header[\"CRVAL1\"]\n",
    "aurora_hdu.header[\"CRVAL2\"] = orig_header[\"CRVAL2\"]\n",
    "aurora_hdu.header[\"CUNIT1\"] = orig_header[\"CUNIT1\"]\n",
    "aurora_hdu.header[\"CUNIT2\"] = orig_header[\"CUNIT2\"]\n",
    "aurora_hdu.header[\"CRPIX1\"] = orig_header[\"CRPIX1\"] - y0\n",
    "aurora_hdu.header[\"CRPIX2\"] = orig_header[\"CRPIX2\"] - x0\n",
    "aurora_hdu.header[\"CDELT1\"] = orig_header[\"CD1_1\"]\n",
    "aurora_hdu.header[\"CDELT2\"] = orig_header[\"CD2_2\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "hdulist = fits.HDUList(fits.PrimaryHDU())\n",
    "hdulist.append(disk_hdu)\n",
    "hdulist.append(aurora_hdu)\n",
    "\n",
    "hdulist.writeto(\"emission_regions.img\", overwrite=True)"
   ]
  }
 ],
 "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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}