{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b35870ee",
   "metadata": {},
   "source": [
    "<h1>An Image Colourization Demonstration Script (Using False Colour)</h1>\n",
    "<p><b>George J. Bendo</b><br/>\n",
    "<b>02 January 2022</b></p>\n",
    "\n",
    "<p>This script converts a FITS image into greyscale and false colour png images.  It is intended to be a demonstration of how professional astronomy images in general can be colourized.</p>\n",
    "\n",
    "<p>The example used below is based on the mid-infrared (24 micron) image of M81 from <a href=\"http://ned.ipac.caltech.edu/uri/NED::Image/fits/2012MNRAS.423..197B/NGC_3031:I:MIPS24:bgm2012\">http://ned.ipac.caltech.edu/uri/NED::Image/fits/2012MNRAS.423..197B/NGC_3031:I:MIPS24:bgm2012</a>.  This file should be downloaded to the same directory as the Jupyter Notebook before running the script.</p>\n",
    "    \n",
    "<p>The script will also work with other FITS files, although some of the lines would need to be modified.  The locations of these lines are indicated in the comments within the code below using the words <i>MODIFICATION OPTION</i>.</p>\n",
    "\n",
    "<p>The following packages need to be installed to use this script: astropy, matplotlib, numpy, and PIL.  The matplotlib and numpy packages are standard python utilities.  The astropy package contains the tools needed to import FITS files.  PIL is used to export images as png files.</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d69637e7",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99635f0f",
   "metadata": {},
   "source": [
    "Perform a series of prepratory steps first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "392a84d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import packages.\n",
    "import numpy\n",
    "import matplotlib.pyplot as pp\n",
    "from astropy.io import fits\n",
    "from PIL import Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a008df5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set functions that describe the colours (r, g, and b).  The input image will\n",
    "# be rescaled to values between 0 and 1000 before these conversions are applied.\n",
    "# These functions will then scale those pixel values into rgb colours that\n",
    "# range from 0 to 255.  These conversion will produce a colours that change\n",
    "# from black through red and yellow to white, producing something like a \"heat\"\n",
    "# colour scale.\n",
    "#\n",
    "# [MODIFICATION OPTION: Change these lines to change the colours.  This can be\n",
    "# quite complicated, but the plots in the next step can be used check how the\n",
    "# science valeus are mapped into red, green, and blue pixel values.  Keep in\n",
    "# mind how these three colours combine to produce other colours.]\n",
    "x=numpy.asarray(range(1001))\n",
    "r=x*0.6\n",
    "g=(x-350)*0.52\n",
    "b=(x-650)*(255./350)\n",
    "\n",
    "# Perform some additional steps to adjust the ranges of the r, g, and b arrays\n",
    "# so that they stay between 0 and 255.\n",
    "r[r>255]=255\n",
    "g[g<0]=0\n",
    "g[g>255]=255\n",
    "b[b<0]=0\n",
    "b[b>255]=255"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9d68a59e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Output PNG Pixel Value')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAFzCAYAAACQKhUCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABN80lEQVR4nO3deZzNZf/H8dcne5I2rZLSSgv3LdUvlVZka1FpsURFSNKifd9TWpRCWUqK7FmjCIVQtsgWRURUKMss1++P68x9zy1mzow5c53l/bwf8/ie851zzryde5r5zHVd389lzjlEREREJJx9QgcQERERSXUqyEREREQCU0EmIiIiEpgKMhEREZHAVJCJiIiIBKaCTERERCSwoqED7I1DDjnEVaxYMXQMERERkVzNnj37N+dcud19LqELsooVKzJr1qzQMURERERyZWar9vQ5TVmKiIiIBKaCTERERCQwFWQiIiIigakgExEREQlMBZmIiIhIYCrIRERERAJTQSYiIiISmAoyERERkcBUkImIiIgEpoJMREREJDAVZCIiIiKBJfRelpIP33wDv/0WOoVISku3DBbst5K1JTaFjiIiwA8Lj+S4w46jYd2LgmVQQZZKFi6EGjVCpxBJWQ4Y3gjufhlWVAqdRiTFZRp8Wh+63ANTzueESz+lYd1wcVSQpZKPPoJ99oHx42G//UKnEUkpC0otp+MxXZlYdhZV/j6OfsubcsL2ozEXOplIatm+vQhjRx7PgH6nsmrlARx+xFaa3Dedli0ODJpLBVmqcA4GDoRateDii0OnEUkZm9jEozxKd7pTlrJ0oxut921N0Ur68StSmH77Dd56C7p1gw0b4F//guefg8aN96No0bNDx4vdon4zO9rMvjCz781soZndGTn/uJmtMbPvIh+XZ3vOA2a2zMx+MLPascqWkubNgyVL4NprQycRSQmZZNKLXpzESXSnO7dzO0tZSjvaUVR/C4sUmqVLoW1bqFABHnvMr9z54guYNQuaNIGicfKfYyxjpAN3O+fmmFkZYLaZfRb5XFfnXJfsDzazykAToApwJDDBzE50zmXEMGPq+PhjKFIErroqdBKRpDeTmbSnPd/wDTWpSTe6cQZnhI4lklKmTYMuXWD4cChWDJo2hU6doHLl0Ml2L2YjZM65tc65OZHbW4BFwFE5PKUR8JFzbodz7kdgGaAV6AXBORg0CC68EMqVC51GJGltYAO3citnczY/8zMf8AFf8qWKMZFCkpkJI0bAuedCzZrw5Zfw0EOwahX06hW/xRgUUh8yM6sIVANmRE61N7N5ZvaemWWtojsK+Dnb01azmwLOzG4zs1lmNmvDhg2xjJ085s+HZcugcePQSUSSUgYZvMmbnMiJ9KEPnejED/zAjdyIYaHjiSS9nTuhb1847TRo1AjWrIHXX4effoKnnoLDDw+dMHcxL8jMbD9gMNDRObcZ6A5UAqoCa4GX8/J6zrkezrnqzrnq5TTaE50hQ8AMrrgidBKRpDONaVSnOu1pz7/4F3OZSxe6sD/7h44mkvS2boWuXaFSJWjRwq/M+eADv27sjjugdOnQCaMX06VsZlYMX4z1d84NAXDO/Zrt8z2BTyN31wBHZ3t6+cg52VuDB/ux28MOC51EJGmsYx2d6Uw/+lGe8gxkII1prBExkUKwYYMfAXvzTfj9dzj/fOjRA+rU8eMPiSiWV1ka8C6wyDn3SrbzR2R72JXAgsjtEUATMythZscCJwAzY5UvZSxZAgsWwNVXh04ikhQyyKAb3TiJkxjAAB7gARazmGu4RsWYSIz9+CO0b++vmHz6abjgAvj6a5g8GerWTdxiDGI7QnYu0BSYb2bfRc49CFxvZlXxTatXAq0BnHMLzWwg8D3+Cs12usKyAAwZ4o+6ulJkr81hDq1pzSxmcSmX0o1unMiJoWOJJL25c+GFF3w7zX328VdM3nsvnHxy6GQFx5xL3DbR1atXd7NmzQodI76deab/7p0xI/fHishubWYzj/AI3ehGOcrxKq9yHddpREwkhpyDSZN8ITZunN9gpk0b6NgRjsqpZ0McM7PZzrnqu/tcoVxlKYGsWuU732m6UiRfHI5P+IRTOIU3eIM2tGExi2lCExVjIjGSmQnDhsHZZ8NFF8G338Kzz8LPP8NLLyVuMZabOOlPKzGh6UqRfPuRH2lPe0YzmqpUZQhDOIuzQscSSVrp6X5K8tlnYeFCOO446N4dmjeHUqVCp4s9jZAlsyFD4PTT4fjjQycRSRhppPE8z1OFKkxmMq/wCt/wjYoxkRjZudM3bT35ZLjxRn+uf3/44Qc/RZkKxRhohCx5rVvn9414/PHQSUQSxjSm0ZrWLGQhV3Ilr/EaR/9PNx4RKSh//+0LsZdegtWroXp1GDoUGjb0S59TTQr+k1PE0KF+RaTWj4nkajObaUtbalKTLWxhBCMYwhAVYyIxsHmzX6h/7LFw553+OHYszJzp+5enYjEGGiFLXsOGwYknxvfGXSJx4FM+5XZuZw1r6EhHnuIp9mO/0LFEks7Gjb6Z6+uvwx9/QO3afp/J884LnSw+pGgdmuQ2b4YvvvDjvoncJU8khtaznuu5ngY04AAO4Gu+pitdVYyJFLB163zPsGOOgSefhFq14Jtv/KiYirH/0ghZMho/HtLSfEEmIv/D4fiAD+hIR7awhSd4gvu5n+IUDx1NJKn89BO8+KJfJ5aWBk2awAMPwKmnhk4Wn1SQJaMRI+Dgg+Gcc0InEYkrq1hFa1ozjnGcwzn0oheV0bS+SEH68UffuqJPHz9J07w5dO6sC/5zo4Is2aSnw6hRUL8+FNX/vSLw3/0nH+IhAF7nddrSliIUCZxMJHmsWAHPPAP9+kGRIr5lxX33wdG6NiYq+o2dbL76CjZt0nSlSMRiFnMzNzOd6dSlLm/zNhWoEDqWSNJYvvy/hVjRotC2rR8RO/LI0MkSixb1J5sRI6B4cbjsstBJRILKIIMudKEqVVnCEt7nfUYxSsWYSAFZtgxatICTToIBA6B9ez9K9tprKsbyQyNkyWbkSLjwQihTJnQSkWCWsIQWtOBrvqYRjXibtzmcw0PHEkkKS5fC00/DBx9AiRLQoYO/ivKII0InS2waIUsmP/wAS5ZoulJSVgYZdKUrZ3AGi1nMB3zAUIaqGBMpAD/8AE2b+i2OBg2Cjh39iNgrr6gYKwgaIUsmI0b4Y4MGYXOIBLCMZdzMzUxlKvWpTw96cAT6LSGytxYv9iNiAwb4EbFOneCee+Cww0InSy4qyJLJiBFQtaouaZGUkkkmb/ImnelMcYrTl740pSmGmiKL7I3Fi+Gpp3whVqoU3H23L8QOPTR0suSkKctk8dtv/gpLTVdKClnBCi7iIjrQgVrUYiELaUYzFWMie2H5cmjWDKpUgeHD/fqwlSt9k1cVY7GjEbJkMXo0ZGZqulJSgsPRi17cxV0UoQjv8i43c7MKMZG98NNPfmryvff8xfqdOvk+YuXKhU6WGlSQJYvRo+Hww+Hf/w6dRCSmfuVXbuVWRjKSi7mY3vTmaDRNL5Jfa9f6zvo9evj7bdv6LY60UL9wqSBLBunpfv/KRo20mbgktZGMpBWt2MxmutKVDnRgH628EMmXDRvghRfgzTf9r5GWLeHhh7UMORQVZMlg5kz4/XeoWzd0EpGY2MpWOtGJnvTkDM7gcz7nVLRDsUh+/P47dOniG7hu2+ZbWTz6KBx3XOhkqU0FWTIYMwb22QcuvTR0EpEC9zVf05SmrGAFnenMEzxBCUqEjiWScDZvhldf9X3D/vwTmjSBxx7zfcUkPI31J4MxY+Ccc+DAA0MnESkwaaTxKI9Sk5qkk84kJvE8z6sYE8mjv/7yV0gee6wvwC66CObN8+0sVIzFDxVkie7XX2H2bE1XSlJZwhLO5Vye4ima0pS5zOV8zg8dSySh7NgBb7zhpyI7d4azzoJvvoEhQ+C000Knk11pyjLRjRvnjyrIJAk4HH3oQ3vaU5KSfMInXM3VoWOJJJSMDPjwQ78ubOVKqFULhg6F//u/0MkkJxohS3Rjxvj9K6pWDZ1EZK/8yZ/cwA20pCU1qME85qkYE8kD52DkSP/roFkzOOggfwH+55+rGEsEKsgSWUaG/6+tTh2/qF8kQU1nOlWpyiAG8QzPMIEJHMVRoWOJJIwpU6BmTb9Zy44d8PHHfnry0kvVDSlR6Ld4Ips5EzZt0nSlJKxMMnmO56hJTRyOKUzhQR6kCEVCRxNJCHPnQr16cP75fnqyRw9YuBCuvVZ/pycarSFLZGp3IQnsF36hKU35nM+5jut4m7c5gANCxxJJCMuX+zViH34IBxzgG7y2bw/77hs6meSXCrJENmYMnH22XyggkkA+5VNa0IJtbNM+lCJ5sG4dPPWUHwkrVsxvcXTvvep6lAw0oJmofvvNt7uoUyd0EpGo7WQnHelIAxpQnvLMZjYtaaliTCQXf/wBDz0ElSr5YuzWW/0o2bPPqhhLFhohS1QTJ/pLai67LHQSkaisZCXXcR0zmUkHOvACL1CSkqFjicS1HTuge3c/KrZpE1x/PTz5JBx/fOhkUtBUkCWqCROgbFmoXj10EpFcjWQkzWlOBhnqLSYSBedg4EA/Jfnjj36p8AsvQLVqoZNJrGjKMhE5B5995ve/KKKr0SR+pZHGfdxHQxpSkYrMYY6KMZFcfPml76rfpAmUKeP7f48fr2Is2akgS0TLl8OqVbq6UuLaalZTi1q8xEvczu18xVdUolLoWCJxa9EiaNQILrgA1q6FPn1gzhytTEkVmrJMRJ995o+XXBI2h8gejGUsN3ETO9jBAAbQhCahI4nErXXr4PHHoVcv37bi2WehY0coVSp0MilMGiFLRBMmQIUKWtUpcSeddB7mYepSlyM5klnMUjEmsgdbt8ITT/gf5e++C+3a+QmQBx5QMZaKNEKWaDIy/MZkV1+t/TAkrmxgA01owud8Tita8Tqvsy/qUimyq/R0eO89eOwxPzrWuDE895z+xk51KsgSzezZviGNpisljsxkJldzNRvYwHu8x83cHDqSSNxxDkaP9o1cFy2Cc8+FoUN9f28RTVkmmqz1YxdfHDaHSERPenIe51GEInzFVyrGRHZj/ny/OL9+fT/RMWSI3xBcxZhkUUGWaCZMgKpVoVy50EkkxW1nO7dyK7dxG7WoxWxm8y/+FTqWSFxZvx7atPE/tmfPhtdegwUL4MortepE/pcKskTy118wbZraXUhwP/ET53EevejFgzzIaEZzMAeHjiUSN3bsgBdfhBNO8Av227eHZcugQwe/B6XIrrSGLJFMmQJpaVo/JkFNZCJNaMIOdjCUoVzBFaEjicQN5/x05H33wYoVfoqySxc46aTQySTeaYQskUyYAMWLQ82aoZNICnI4XuIlLuMyDuVQvuEbFWMi2cyeDbVq+asm993Xd9cfOVLFmERHBVkimTQJzjnH/5cuUoi2sY2buIn7uI+ruIrpTOck9FtGBOCXX+Dmm+HMM/3Vk2+/Dd9+q9UlkjcqyBLFH3/4/8Jr1QqdRFLMalZzHucxgAE8y7MMZCBlKBM6lkhwf/8NTz0FJ54IH34I99wDS5dC69ZQVAuCJI/0LZMopk6FzEwVZFKopjOdK7mSrWxlOMNpQIPQkUSCcw4GDfIF2M8/w1VX+QX8lbRVq+wFjZAlikmToEQJNa2RQtOHPlzABZSmNNOZrmJMBJg3Dy68EK67Dg46yP9oHjxYxZjsPRVkiWLSJF+MlSwZOokkuXTS6UQnbuZmzuM8ZjKTKlQJHUskqE2b4I47oFo13+S1e3e/iP+CC0Ink2QRs4LMzI42sy/M7HszW2hmd0bOH2Rmn5nZ0sjxwMh5M7PXzWyZmc0zM3WYzKL1Y1JIfud36lGPrnSlAx0Yy1gO4qDQsUSCyciAd97x68Teegtuv92vE2vTBooUCZ1OkkksR8jSgbudc5WBs4F2ZlYZuB+Y6Jw7AZgYuQ9QFzgh8nEb0D2G2RKL1o9JIVjCEs7iLL7gC3rRi9d4jaJaZiopbNo0f+VkmzZw6qn+7+Ju3fxUpUhBi1lB5pxb65ybE7m9BVgEHAU0AvpGHtYX/tPIqBHQz3nTgQPM7IhY5UsoWj8mMTaZyZzN2fzBH3zO57SiVehIIsGsWQM33eRbPm7YAB9/DF98AaefHjqZJLNCWUNmZhWBasAM4DDn3NrIp9YBh0VuHwX8nO1pqyPndn2t28xslpnN2rBhQ+xCx5Os/mNaPyYx0Je+XMqlHM7hTGc6NVHjYUlNO3bACy/4Rq6ffAIPPwyLF8O112rfSYm9mBdkZrYfMBjo6JzbnP1zzjkHuLy8nnOuh3OuunOuerlU2GBb68ckRjLJ5GEepgUtOJ/z+YqvOI7jQscSCWLUKD8tef/9vqHr99/7HmOlS4dOJqkipgWZmRXDF2P9nXNDIqd/zZqKjBzXR86vAY7O9vTykXOpTevHJAa2s50buIFneIZbuIUxjOEADggdS6TQrVwJjRr5PSeLFoVx42DoUDhOf5tIIYvlVZYGvAsscs69ku1TI4DmkdvNgeHZzjeLXG15NvBntqnN1JW1fuyss0InkSSxnvVcxEUMZCAv8iI96EExioWOJVKoduyAZ56BypVh4kTf2HXuXLjsstDJJFXF8hKqc4GmwHwz+y5y7kHgeWCgmbUCVgHXRj43GrgcWAb8Ddwcw2yJQ+vHpAB9z/fUpz7rWMcnfMJVXBU6kkihmzAB2rWDJUv8RuBdu0L58qFTSaqLWUHmnJsK7GkZ5MW7ebwD2sUqT0LassWvH3voodBJJAlMZCJXczWlKMVkJnMmZ4aOJFKo1qyBTp1g4EA4/ngYOxZq1w6dSsRTp/54Nn26Xz9WU1e9yd7pT3/qUIejOZoZzFAxJiklLQ1efhlOPhlGjIAnn/Td9lWMSTxRQRbPpk6FffZR/zHJN4fjRV7kJm7iPM5jKlOpQIXQsUQKzZQp8K9/+Y3AL7gAFi6ERx7RKhCJPyrI4tm0aXDGGbD//qGTSALKIIM7uZPOdKYJTRjDGMpSNnQskULx66/QvDmcf75f/TFsGIwcqasnJX6pIItXaWl+ylLTlZIP29lOE5rwBm/QiU70pz8lKBE6lkjMZWb6jb9POgkGDIAHH/Q9xRo1UnNXiW/aqC5ezZ0Lf/2lgkzy7Hd+pxGNmMIUXuEV7uKu0JFECsW8eXDbbTBjBlx0Ebz5pl83JpIINEIWr6ZO9cdzzw2bQxLKz/xMTWoygxl8xEcqxiQl/P2377D/73/D8uXw/vu+tYWKMUkkGiGLV1OnQsWKcNQ/tvMU2a35zKcuddnCFsYylgu5MHQkkZgbOxZuv9133G/Z0jd4Pfjg0KlE8k4jZPHIOV+QabpSovQVX3E+5+NwTGWqijFJeuvWwfXXQ926fjOTSZPg3XdVjEniUkEWj5Yv95cIqSCTKIxjHJdwCeUox1d8xWmcFjqSSMxkZkKPHnDKKTBkCDzxhF9ye8EFoZOJ7B1NWcajadP8UQWZ5GIgA7mJm6hCFcYxjkM5NHQkkZhZsABat4avvoJateDtt/3VlCLJQCNk8WjqVDjgAP8noMge9KQnTWjCWZzFF3yhYkyS1rZtfge5atXghx+gTx/4/HMVY5JcNEIWj6ZO9VdX7qN6WXbvRV6kM52pS10+4RP2Zd/QkURiYtIkuPVWWLbMN3rt0gUOOSR0KpGCp9/48ea332DxYrW7kN1yOO7nfjrTmeu4jmEMUzEmSenPP/305IUX+uucJk70I2MqxiRZqSCLNzNm+KMKMtlFBhm0pjUv8AJtaEN/+lOc4qFjiRS4kSOhcmXo1cvvQTlvnm/0KpLMVJDFm+nToUgR3+FQJCKNNG7iJnrSkwd4gLd4iyIUCR1LpECtX+9bWTRs6NtXzJgBL70E+2oQWFKACrJ4M306nHYalC4dOonEiZ3spAlN+IiPeJ7neZZnMbQpnyQP56B/fz8qNmQIPPUUzJoF1auHTiZSeFSQxZPMTJg5E84+O3QSiRPb2c7VXM0QhtCVrnSmc+hIIgXq55+hfn246SY48UT49lt4+GEortl4STEqyOLJ4sWweTOcdVboJBIHtrGNRjTiUz7lLd6iIx1DRxIpMJmZ8NZbflRs0iR47TWYMsXfF0lFansRT6ZP90eNkKW8v/iLBjRgEpN4l3dpScvQkUQKzJIlcMstvgC79FLfeb9ixdCpRMLSCFk8mTHDN4Q98cTQSSSgzWymDnWYzGT60U/FmCSNjAx45RU44wzfdb9PHxg3TsWYCGiELL5Mnw41aqghbAr7gz+oQx1mM5uP+IhruCZ0JJECsXQp3Hyz3xmuYUO/7dERR4ROJRI/9Js/Xmzd6v9k1HRlyvqd37mES5jDHD7hExVjkhQyM+HVV/2o2MKF8P77MGyYijGRXUVVkJlZTTO7OXK7nJkdG9tYKWjWLP+TSwVZSvqDP7iMy5jPfIYxjEY0Ch1JZK8tW+Y3Ab/rLrj4Yl+Q3XQTmLq2iPxDrgWZmT0GdAYeiJwqBnwQy1ApKWtBf40aYXNIofuTP6lNbeYylyEM4XIuDx1JZK9kZsIbb8Dpp/su+336wIgRcOSRoZOJxK9o1pBdCVQD5gA4534xszIxTZWKZsyAE07w7aklZWxhC3WpyxzmMJjB1KNe6Egie2XFCmjZEiZPhrp1oWdPOOqo0KlE4l80U5Y7nXMOcABmphbyBc05P0Km6cqUspWt1KUuM5nJQAbSkIahI4nkW2YmvPmm32jk22/hvfdg1CgVYyLRimaEbKCZvQMcYGa3Ai2BnrGNlWJ++gnWrVND2BTyF39Rj3pMZzoDGMCVXBk6kki+rVrlr6D84guoXduPih19dOhUIokl14LMOdfFzC4FNgMnAY865z6LebJUMnOmP6ogSwl/8zf1qc9UptKf/rqaUhKWc9CvH3To4EfIevTwDV+1aF8k76LqQxYpwFSExcqsWX7jttNOC51EYmwb22hIQ77kS/rRjyY0CR1JJF82bIDWrWHoUDjvPOjbF47V9fci+RbNVZZbzGxz5GO7mWWY2ebCCJcyZs3ylyOVKBE6icTQTnZyDdfwOZ/Tm97cyI2hI4nky6ef+r8fR42CF1/0U5UqxkT2TjRTlv+5otLMDGgEaPV5QXEOZs+GJhopSWYZZNCMZoxiFN3pTjOahY4kkmdbtkCnTtCrl/8bcvx4fxSRvZenTv3OGwbUjk2cFLR8Ofz5J1SvHjqJxIjD0YY2fMzHvMiLtKFN6EgieTZtGlStCu++C507+6WvKsZECk6uI2RmdlW2u/sA1YHtMUuUambN8kcVZEnJ4bibu+lFLx7mYe7l3tCRRPJk50547DE/NXnMMfDll1CzZuhUIsknmkX9DbLdTgdWgvZ1KTCzZ/u1Y1WqhE4iMfAkT9KVrtzBHTzJk6HjiOTJggXQtCl89x20agVdu0IZtQUXiYlo1pDdXBhBUtasWX7X3WLFQieRAtaVrjzO47SgBa/yKoZ6AUhicA5ef91PTe6/PwwfDg3Vt1gkpvZYkJnZG0S68++Oc65DTBKlksxMP0LWtGnoJFLAetGLTnSiMY3pSU/2ydtyTZFgfv0VWrSAsWOhfn2/ZuzQQ0OnEkl+OY2QzSq0FKlq2TJ/2dK//x06iRSgIQyhNa2pQx3605+i0bX7Ewlu1CjfcX/LFr8N0u23q8mrSGHZ428K51zfwgySkrSgP+lMZjI3cANncRaDGUxxioeOJJKrbdvgvvugWzd/5eQXX2hZq0hhi+Yqy3JAZ6AyUDLrvHPuohjmSg2zZkHJklC5cugkUgDmMY+GNOQ4jmMkI9mXfUNHEsnV/Plwww1+AX/HjvDcc/7HkogUrmgWtvQHFgHHAk/gr7L8JoaZUsfs2b6xT1FNaSW6laykDnUoQxnGMY6DOTh0JJEcZS3cP/NMvw3SmDH+KkoVYyJhRFOQHeycexdIc85Nds61BDQ6trcyM2HOHE1XJoENbKA2tdnGNsYxjqM5OnQkkRz9+ivUqwd33gmXXALz5kGdOqFTiaS2aIZm0iLHtWZWD/gFOCh2kVLEkiWwdasW9Ce4rWylPvX5iZ+YwASqoIU3Et/GjPFXUW7e7NeMtW2rhfsi8SCnthfFnHNpwNNmVha4G3gD2B+4q5DyJa+sBf0qyBJWGmk0pjGzmMVQhnIu54aOJLJHO3fCQw9Bly5+Y/CJE+HUU0OnEpEsOY2QrTGzEcAAYLNzbgFwYeHESgHffec79J9ySugkkg8Ox63cyjjG0ZOeNERdMyV+/fgjNGni95+8/XZ4+WUoVSp0KhHJLqc1ZKfgF+8/DPxsZq+Z2dmFEysFfPed/zNVC/oT0lM8RV/68jiPcwu3hI4jskeffALVqsHixTBoELz1looxkXi0x4LMObfROfeOc+5CoAawAuhqZsvN7JlCS5iMnPMFWdWqoZNIPrzP+zzGYzSnOY/yaOg4Iru1fbtfH3bNNXDSSf5HTuPGoVOJyJ5EtZ+Lc+4X4F2gO7AFNCSwV9asgY0bVZAloElMohWtuJAL6UEP7U8pcWnxYjjrLOjeHe65B6ZMgWOPDZ1KRHKSY0FmZiXN7BozGwIsw7e7uB84sjDCJa25c/1RBVlCWcQiruRKjud4hjBEXfglLvXt67vprFnjt0J66SUorm9VkbiX01WWHwKXAJPxzWFvcM5tL6xgSe277/zx9NODxpDorWc99ahHcYozmtEcwAGhI4n8j61boV076NcPLrgA+veHo44KnUpEopXTCNlYoJJz7hrn3OC8FmNm9p6ZrTezBdnOPW5ma8zsu8jH5dk+94CZLTOzH8ysdt7/KQnku++gUiUoUyZ0EonCNrbRkIasYx0jGUlFKoaOJPI/Fizwo2Lvvw+PPeZbWqgYE0ksOW0u3m8vX7sP0A3Y9XW6Oue6ZD9hZpWBJkAV/HToBDM70TmXsZcZ4pMW9CeMTDK5iZuYyUyGMIQa1AgdSeR/9OsHbdpA2bK+ELtQzYlEElJUi/rzwzn3JbApyoc3Aj5yzu1wzv2IX6+WnL/5tmyBZctUkCWIB3mQIQzhZV7mCq4IHUfkP7Zvh9tug+bN/QL+b79VMSaSyGJWkOWgvZnNi0xpHhg5dxTwc7bHrI6cSz7z5/ujCrK49wEf8AIv0JrWdKRj6Dgi/7F8Ofzf/0HPnvDAA/DZZ3D44aFTicjeyGlR/1U5PdE5NyQfX6878BTgIseXgZZ5eQEzuw24DaBChQr5iBBY1oL+M84IGkNyNoMZ3MIt1KIWb/CG2ltI3Bg+3I+KmcHIkVC/fuhEIlIQcmoT3yCHzzkgzwWZc+7XrNtm1hP4NHJ3DXB0toeWj5zb3Wv0AHoAVK9e3eU1Q3DffQcHHQTly4dOInuwmtVcwRUcyZEMYhDFKBY6kghpaX4vypde8lvgDhqk3mIiySSnRf03F/QXM7MjnHNrI3evBLKuwBwBfGhmr+AX9Z8AzCzorx8Xshb0m0Zc4tHf/M0VXMFWtjKBCRzCIaEjifDLL34vyilT/F6UXbv6rXBFJHnkupGimR0GPAsc6ZyrG7ki8hzn3Lu5PG8AUAs4xMxWA48BtcysKn6EbSXQGsA5t9DMBgLfA+lAu6S8wjI93a8ha9s2dBLZDYejJS2ZwxyGM5wqVAkdSYTPP4frr/d9xvr3hxtuCJ1IRGIhmp2t+wC9gYci95cAH+O3Utoj59z1uzm9x+c4554BknuPzKVL/aVRWtAfl57hGT7mY57neRrkOGMvEnuZmfD88/DII3DiifDFF1C5cuhUIhIr0VxleYhzbiCQCeCcSweSb/SqMMyb54/q0B93hjKUR3iEm7iJ+7gvdBxJcX/+CVde6deMXXcdfPONijGRZBfNCNlfZnYwfpoRMzsb+DOmqZLVggVQpAicfHLoJJLN93xPU5pSgxr0pKeuqJSgFi2CK67wrS1efx3at9eSU5FUEE1Bdjd+0X0lM5sGlAOuiWmqZDV/vp970GrcuPEnf3IlV7If+zGEIZSkZOhIksKGDPEtLfbd168dO//80IlEpLDkWpA552ab2QXASYABP0SOklcLFsC//hU6hURkkklzmrOc5XzO5xyVpL2IJf5lZPi1Ys8957vuf/KJOuOIpJpc15CZ2SSgvHNuoXNuAVAV+CbGuZLPX3/BihVw6qmhk0jE8zzPcIbzMi9zPhqKkDA2bYJ69XwxduutMHmyijGRVBTNlOVzwFgzex2/ndHlQIH3KEt6ixaBc3DaaaGTCDCe8TzMw9zADXSgQ+g4kqLmzvWL99esgR49fEEmIqkpminLcWbWBvgM+A2o5pxbF/NkyWZBpAeuRsiCW8lKrud6TuVUetBDi/gliA8/hFtu8Rt3TJ4MZ58dOpGIhBTNlOUjwBvA+cDjwCQzqxfjXMlnwQIoWRKOOy50kpS2jW1cxVVkkslQhlKa0qEjSYpJT4dOneDGG6F6dZg9W8WYiEQ3ZXkwUMM5tw342szGAr2AUTFNlmzmz/eNhIoUCZ0kZTkct3M73/Itn/IplagUOpKkmI0b4ZprfJPXO+6Al1+GYtoqVUSIbsqy4y73VwGXxipQ0lqwAC65JHSKlPYu79KXvjzGY9RDg7xSuBYuhAYN/L6UfftCs2ahE4lIPNljQWZmrzrnOprZSCJNYbNzzjWMabJksmmT/ymsBf3BzGMed3AHl3Ipj/BI6DiSYkaO9HtQ7refXy921lmhE4lIvMlphOz9yLFLYQRJagsX+qMW9AexhS1cwzUcyIF8wAcUQdPGUjicgxdegAcf9C0Ihw1TSwsR2b2cCrKFZtYROB6YD7wb2cdS8kpXWAbjcLSmNctYxud8zqEcGjqSpIht2/xVlB9+CE2awHvvQalSoVOJSLzK6SrLvkB1fDFWF3i5UBIlowULoGxZOEqd4AtbT3oygAE8yZNcwAWh40iKWLMGLrjAF2PPPOOPKsZEJCc5jZBVds6dBmBm7wIzCydSEpo/34+OaYfgQvUd39GBDlzGZTzAA6HjSIqYOdNvDr5li5+ibNQodCIRSQQ5jZClZd3QVOVecM6PkGlBf6HazGau5VoO5mA+4AP2yb3lnshe69/fbwhesiR8/bWKMRGJXk4jZGeY2ebIbQNKRe4b4Jxz+8c8XTJYtw5+/x2qVAmdJGVkrRtbznImMYlylAsdSZJcZiY89BA8/zzUqgWDBsEhh4ROJSKJZI8FmXNOl6IVhMWL/fGUU8LmSCG96c1HfMQzPMN5nBc6jiS5v/+Gm26CoUOhTRt4/XU1exWRvIumU7/sjUWL/PHkk8PmSBFLWEIHOnAhF9KZzqHjSJJbuxYaNoQ5c+C116CD9qkXkXxSQRZrixZBmTJw5JGhkyS9nezkBm6gBCXoRz/1G5OYmjcP6tf3fZ+HD/e3RUTySwVZrC1e7EfHdIVlzD3CI8xmNkMYQnnUfVNiZ/RouO46381m6lSoWjV0IhFJdLr0LNYWLdL6sUIwgQm8yIvcxm1cyZWh40gSe/NNvyflCSfAjBkqxkSkYOyxIDOzVmZ2b7b7a8xss5ltMbM2hRMvwW3Z4jtEqiCLqd/4jWY042RO5hVeCR1HklRGBtx5J7Rv76cnv/xSvZ5FpODkNELWBngv2/31kVYX5YDrY5oqWWRdYakF/THjcLSiFRvZyId8SGlKh44kSWjLFt9T7PXXoVMnGDLEbxQuIlJQclpDZs65jdnuDwJwzm03M20CEo2sKyw1QhYz7/AOIxjBK7xCNaqFjiNJaPVqPyK2YAG89RbcfnvoRCKSjHIqyA7Ifsc59yyAme0DqOVhNBYv9g2JjjsudJKktIxl3M3dXMql3MmdoeNIEvruO6hXz4+QjRoFtWuHTiQiySqnKcvxZvb0bs4/CYyPUZ7ksmgRHH+8ukTGQAYZNKc5xSjGe7ynrZGkwH32md8GaZ994KuvVIyJSGzlNEJ2L9DLzJYBcyPnzgBmAbfEOlhSWLRIWybFSBe68BVf8T7vq8WFFLj334eWLaFyZd/iQov3RSTWcto66S/gejM7DsiqKr53zi0vlGSJLi0Nli+Hxo1DJ0k685jHozzK1VzNjdwYOo4kEef8fpQPPggXXwyDB/teYyIisbbHgszMKkRupvPfEbL/nHfO/RTbaAlu2TJIT9cVlgVsJztpRjMO4AC60x1DDXelYGRkwB13QPfucOON8N57ULx46FQikipymrIcBTj4n994Dt/24lDQvjQ50hWWMfEETzCXuQxnOOUoFzqOJIm//4brr4cRI+D+++GZZ/zaMRGRwpLTlOVp2e+bWUWgM3AJ8GxsYyWBrB5kJ50UNkcSmc50nud5WtCChjQMHUeSxIYNvvP+zJnQrRu0axc6kYikolz3sjSzE4CHgLOAl4EOzrm0WAdLeIsWwdFHq3tkAfmbv2lOc8pTnld5NXQcSRLLl0OdOr7X2JAhcMUVoROJSKrKaQ3ZqfhCrArwItDKOZdRWMES3pIlGh0rQA/zMEtYwkQmUhatspa9N2sWXH65Xzs2cSL83/+FTiQiqSynEbK5wM/4tWQ1gBpm/11O5pzrENtoCcw5X5Bdrx2mCsJ0pvMqr3I7t3MRF4WOI0lgwgQ/GlauHIwdq7+dRCS8nAqyVvhF/JJXGzfCH3/AiSeGTpLwdrCDVrSiPOV5nudDx5EkMGiQv4ry5JNh3Dg44ojQiUREcl7U38fMygHHAMucc38UWqpEt3SpP55wQtgcSeAZnuF7vmc0o9mf/UPHkQT39tvQtq2fnhw5Eg48MHQiERFvjxd2m9ktwELgDWCxmemytmipICsQ85jHczxHU5pSl7qh40gCcw6eftpvDH755TB+vIoxEYkvOU1ZdgSqOOc2RLr19wdGFEqqRLdkCRQpAsceGzpJwkonnZa05CAOoitdQ8eRBJaZCXfdBa+/Ds2aQa9e2l5WROJPTgXZTufcBgDn3AozK1FImRLf0qW+GNNP/XzrSldmM5uBDORgDg4dRxJUWhrcfDP07w+dOsFLL6nhq4jEp5wKsvJm9vqe7usqyxwsXarpyr2whCU8yqNcwRU0RnuBSv789Rdccw2MGQPPPQedO4Nppy0RiVM5FWT37nJ/diyDJA3nfEF23nmhkySkTDK5lVspQQne4i3tVSn5smkT1K8PM2ZAz55wyy2hE4mI5Cynqyz7FmaQpLFuHWzdqhGyfOpDH77kS3rSkyNQPwLJu3Xr4NJL/VLOQYPgqqtCJxIRyV2uWydJHukKy3zbwAbu5V5qUpOWtAwdRxLQTz/BxRfD2rV+qvIi9REWkQShgqygZRVkagqbZ/dxH5vZzNu8zT577sgisltLl8Ill8Cff8Jnn8E554ROJCISvVx/65nZudGck4ilS6F4cahQIXSShDKZyfShD/dwD1WoEjqOJJgFC+D88+Hvv+GLL1SMiUjiiWYY4o0ozwn4hSvHHef7kElUdrCDNrThWI7lER4JHUcSzOzZcMEF/grKyZOhWrXQiURE8m6PU5Zmdg7wf0A5M+uU7VP7A6o29kQtL/KsC11YzGJGMYp92Td0HEkgU6dCvXq+6/7EiVCpUuhEIiL5k9MIWXFgP3zRVibbx2ZQc6jdysyEZcu0fiwPlrOcp3maxjTmci4PHUcSyIQJULs2HH44TJmiYkxEEltObS8mA5PNrI9zblVeX9jM3gPqA+udc6dGzh0EfAxUBFYC1zrnfjczA14DLgf+Blo45+bk9WsGt2YNbN+uEbIoORxtaUsxivEqr4aOIwlkxAjf9PXkk/2+lIcdFjqRiMjeiWYNWR8z+3zXj2ieB9TZ5dz9wETn3AnAxMh9gLrACZGP24DuUaWPN8uW+ePxx4fNkSAGM5jxjOdpnuYojgodRxLExx/73mJVq/oF/CrGRCQZRNP24p5st0sCVwPpuT3JOfelmVXc5XQjoFbkdl9gEtA5cr6fc84B083sADM7wjm3Nop88WPFCn/U3Emu/uZv7uZuTud02tI2dBxJEB98AM2bw7nnwqefwv77h04kIlIwci3InHO7bpk0zcxm5vPrHZatyFoHZP1texTwc7bHrY6cS7yCrGhRKF8+dJK49wIv8BM/8T7vU1Tt8CQK/fpBixZQqxaMHAmlS4dOJCJScHL9TRhZ95VlH+DfQNm9/cLOOWdmLq/PM7Pb8NOaVIi3Xl8rVsAxx/iiTPboR37kBV6gCU04n/NDx5EE0Ls3tGrlu/APHw776mJcEUky0VQOswEHGH6q8kegVT6/3q9ZU5FmdgSwPnJ+DXB0tseVj5z7B+dcD6AHQPXq1fNc0MXUihW+B5nkqBOdKEIRXuKl0FEkAfTqBbfd5venHDYMSpUKnUhEpODluqjfOXesc+64yPEE59xlzrmp+fx6I4DmkdvNgeHZzjcz72zgz4RbPwYqyKIwnvEMYxgP8zDl0dSu5Oydd+DWW317i+HDVYyJSPKKZsqyJNAWqIkfKZsCvO2c257L8wbgF/AfYmargceA54GBZtYKWAVcG3n4aHzLi2X4thc35+cfE9TmzfDbbyrIcrCTnXSgA8dzPJ3olPsTJKW99Ra0a+cbvw4eDCVKhE4kIhI70UxZ9gO28N/tkm4A3geuyelJzrnr9/Cpi3fzWAe0iyJL/PrxR39UQbZHb/AGP/ADn/IpJdBvV9mzN96ADh2gQQMYNEjFmIgkv2gKslOdc5Wz3f/CzL6PVaCEldXyQgXZbv3KrzzBE1zO5dSjXug4EsdefRXuuguuuML3HCtePHQiEZHYi6Yx7JzIui4AzOwsYFbsIiUoFWQ5epRH2cY2utI1dBSJY6+95ouxq6+GgQNVjIlI6ohmhOzfwFdm9lPkfgXgBzObj59tPD1m6RLJihV+h+MDDgidJO4sYAG96MUd3MGJaJ9P2b0334SOHX0X/gEDoFix0IlERApPNAXZrtsfye7oCss9uod72J/9eZRHQ0eRONWzJ7RvDw0bqhgTkdQUTUH2tHOuafYTZvb+rudS3ooVcMYZoVPEnbGMZRzjeIVXOIiDcn+CpJw+faB1a6hbV9OUIpK6ollDViX7HTMrip/GlCwZGbBypUbIdpFOOvdwD5WoRLsEv4hWYqN/f2jZEi65BIYM0dWUIpK69jhCZmYPAA8CpcxsM75TP8BOIp3yJeKXX2DnThVku3iP91jIQj7hE4qjYQ/5XwMHQrNmfm/KYcOgZMnQiUREwtnjCJlz7jnnXBngJefc/s65MpGPg51zDxRixvinKyz/YQtbeIRHqElNruKq0HEkzgwdCjfcAOee6zcK196UIpLqollDNsbM/rEDtHPuyxjkSUwqyP7hBV5gPesZyUjsP4OrIr4Au+46qFEDRo2C0qVDJxIRCS+aguzebLdLAjXwG45fFJNEiWjFCihSBI4+OvfHpoCf+ZmXeZkbuIEa1AgdR+LIuHHQuLG//mXMGChTJnQiEZH4kGtB5pxrkP2+mR0NvBqrQAlpxQpfjOlafQCe4AkyyeRZng0dReLIlClw5ZVQuTKMHw9ly4ZOJCISP6K5ynJXq4FTCjpIQlu1CipWDJ0iLixmMb3pTVvacgzHhI4jcWL2bL9J+DHH+GLswANDJxIRiS+5jpCZ2RuAi9zdB6gKzIlhpsSzahVc/I8901PSwzzMvuzLgzwYOorEie+/h9q14aCD4LPPoFy50IlEROJPNGvIsu9bmQ4McM5Ni1GexJOW5tteVKgQOklw3/ANgxnM4zxOOfRbV+DHH+HSS/1s/oQJUL586EQiIvEpmoLsY+D4yO1lzrntMcyTeFavhsxMPxeT4h7gAQ7hEDrRKXQUiQO//OIbvm7fDpMnw/HH5/4cEZFUlVNj2KLAs0BLYBW+MezRZtYbeMg5l1Y4EePcqlX+mOIF2QQmMJGJvMqrlEGXzqW6337zI2MbNsDEiXDqqaETiYjEt5wW9b8EHAQc65z7t3PuX0Al4ACgSyFkSww//eSPKVyQORwP8AAVqEAb2oSOI4Ft3gx16viLj0eOhDPPDJ1IRCT+5TRlWR840TmXtaAf59xmM7sdWAzcGetwCSFrhCyFe5ANYQizmEVvelMCbUaYyv7+Gxo0gLlzYfhwuOCC0IlERBJDTgWZy16MZTuZYWb/OJ+yVq2Cww9P2Y340knnIR6iMpVpStPQcSSgtDS49lqYOhUGDIDLLw+dSEQkceQ0Zfm9mTXb9aSZ3YQfIRPwBVkKT1d+yIf8wA88xVMUoUjoOBKIc3DrrX4rpO7dfWEmIiLRy2mErB0wxMxa4rdKAqgOlAKujHWwhLFqFVStGjpFEOmk8xRPUZWqXKlviZR2//3Qty88+STcdlvoNCIiiWePBZlzbg1wlpldBFSJnB7tnJtYKMkSQWamX9TfqFHoJEF8wAcsYxnDGKYNxFPYK6/Aiy9Cu3bw8MOh04iIJKZo9rL8HPi8ELIkng0bYMeOlJyyTCedp3maalSjIQ1Dx5FA3n8f7r4brrkGXnsNTHW5iEi+RNMYVvYkhXuQvc/7LGc5wxmu0bEUNWYMtGwJF13kC7MiWkIoIpJv+dlcXLKkaEGWRhpP8RT/5t80oEHoOBLAjBnQuDGcfjoMHQol1O1ERGSvaIRsb6RoQdaPfvzIj7zO6xodS0GLFvmWFkccAaNHw/77h04kIpL4NEK2N1atgrJl/UeKSCONp3ma6lSnHvVCx5FCtmYN1K7tNwsfPx4OOyx0IhGR5KARsr2xahVUqBA6RaHqS19WspJudNPoWIrZvNmPjP3xh98s/LjjQicSEUkeKsj2xk8/pdR0ZTrpPMdzVKc6l6M27KkkLc2vGfv+e9/8tVq10IlERJKLpiz3xk8/pdQI2UAGsoIVPMRDGh1LIc75Zq+ffQY9e8Jll4VOJCKSfFSQ5ddff8Hvv6fMpuKZZPIcz1GZyuo7lmKeeAL69IHHH4cWLQKHERFJUpqyzK81a/yxfPmwOQrJp3zKAhbwPu+zj+r4lNG7ty/IWrSARx8NnUZEJHnpN2t+rV7tjylQkDkcz/AMFalIE5qEjiOFZPx4P1V56aXQo4e68IuIxJJGyPIrhQqyL/iCmcykO90pqm+ZlDB3rl/EX7kyfPKJb3MhIiKxoxGy/MqasjzqqLA5CsGzPMvhHE4LWoSOIoXg5599e4uyZdX4VUSksGi4I79Wr4aDDoJSpUIniakZzGAiE3mJlyhJydBxJMa2bIF69WDrVpg6NSX+3hARiQsqyPJr9eqUmK58juc4kANpTevQUSTGMjLg+ut9r7HRo+G000InEhFJHZqyzK8UKMgWspDhDKcDHShDmdBxJMbuvdc3fX3jDfUaExEpbCrI8isFCrKXeZlSlKI97UNHkRh75x3o2hU6dIDbbw+dRkQk9aggy48dO2D9+qReYLOOdfSnPy1owSEcEjqOxNCECdCuHdStCy+/HDqNiEhqUkGWH2vX+mMSj5B1oxtppHEXd4WOIjG0eLFvb3HyyfDRR1BUq0pFRIJQQZYfSd6D7C/+ojvdaUQjTuCE0HEkRjZuhPr1oXhx+PRTtbcQEQlJfw/nR5IXZH3owyY2cTd3h44iMbJzJ1x1lf9W/vxzqFgxdCIRkdSmgiw/krggyyCDrnTlLM7iXM4NHUdiwDlo0wa+/BL694f/+7/QiURERFOW+bF6Ney3X1LO8YxgBMtZzt3cjaHNC5NRly5+0/BHH4UbbgidRkREQAVZ/iRxy4uXeZmKVORKrgwdRWJg7Fjo3BmuuQYefzx0GhERyaIpy/xYsyYpC7LpTGca03iN17SJeBJasgSaNIHTT/cjZKYBUBGRuKERsvxI0hGyrnSlLGVpScvQUaSAbd4MjRpBsWIwbBiULh06kYiIZKdhkLzKyPB9yJKsKewa1jCYwdzJnezHfqHjSAHKzIQbb4SlS30TWF1RKSISf4IUZGa2EtgCZADpzrnqZnYQ8DFQEVgJXOuc+z1EvhytX++LsiQryN7mbTLJpB3tQkeRAvboo77PWLduUKtW6DQiIrI7IacsL3TOVXXOVY/cvx+Y6Jw7AZgYuR9/srr0H3542BwFaAc76EEP6lGP4zgudBwpQIMGwTPPQKtW0LZt6DQiIrIn8bSGrBHQN3K7L3BFuCg5yCrIjjgibI4CNJCBrGc9d3BH6ChSgObOhRYt4Jxz4M03tYhfRCSehSrIHDDezGab2W2Rc4c55yLVDuuAw8JEy0USFmTd6MZJnMQlXBI6ihSQ337zi/gPPBAGD4YSJUInEhGRnIRa1F/TObfGzA4FPjOzxdk/6ZxzZuZ298RIAXcbQIUKFWKfdFfr1vljkkxZzoz87w3eYJ+4GjCV/EpP933G1q2DKVOS6m8HEZGkFeQ3sHNuTeS4HhgK1AB+NbMjACLH9Xt4bg/nXHXnXPVy5coVVuT/WrsWDjooaYYc3uANylCG5jQPHUUKyP33w6RJ0KMHnHlm6DQiIhKNQi/IzKy0mZXJug1cBiwARsB/qoLmwPDCzhaVtWuTZsjhV37lYz6mBS0oQ5nQcaQADBoEL78M7dpBs2ah04iISLRCTFkeBgw1v8K4KPChc26smX0DDDSzVsAq4NoA2XK3dm3STFf2oAdppNGe9qGjSAFYtAhatoSzz4ZXXgmdRkRE8qLQCzLn3ArgjN2c3whcXNh58mztWjjvvNAp9lo66fSgB5dyKSdyYug4spe2bIGrr4ZSpfwoWfHioROJiEheqFN/XjiXNFOWYxjDalbzGq+FjiJ7yTnfZ+yHH+Czz5JyVy8RkaSngiwv/vgDdu5MioLsHd7hcA6nAQ1CR5G99OqrflTshRfgootCpxERkfxQn4O8SJIeZD/xE2MYQytaUYxioePIXvjyS7j3XrjySn8UEZHEpIIsL5KkIHuXd3E4buGW0FFkL6xdC9ddB5UqQe/e6sQvIpLINGWZF0mwj2U66fSiF7WpTUUqho4j+ZSWBtdeC5s3w4QJULZs6EQiIrI3VJDlRRKMkI1iFL/wC2/yZugoshcefBCmToUPP4QqVUKnERGRvaUpy7xYtw723RfKJG4T1Xd4hyM5kvrUDx1F8mnkSOjSxTd/vf760GlERKQgqCDLi6yWFwm6WGcVqxjLWFrRiqIaHE1IP/0EzZtDtWq+KBMRkeSggiwvErwHWS96AWgxf4JKS/OL+NPTYeBAKFkydCIRESkoKsjyIoELsgwyeI/3qEtdKlAhdBzJhwcfhOnToVcvOP740GlERKQgqSDLi3XrEvYKy/GM5xd+oSUtQ0eRfPj0Uz9F2batv7pSRESSiwqyaO3YAX/+CYcdFjpJvvSmNwdzsDrzJ6CsdWNVq8LLL4dOIyIisaCCLFobNvjjoYeGzZEPm9jEcIZzIzdSHO06nUjS0qBJE3/UujERkeSlS+2itX69PyZgQfYhH7KTndzMzaGjSB499BB8/TV89BGccELoNCIiEisaIYtWAhdkfehD1cj/JHGMHg0vvQS33+6vrhQRkeSlgixaWQVZuXJhc+TRfOYzm9kaHUsw69ZBixZw+unwyiuh04iISKypIItWgo6Q9aY3xSjGDdwQOopEKTPTL+LfuhUGDNC6MRGRVKA1ZNFavx5KlEiobZPSSOMDPqABDTiEQ0LHkSi9+iqMHw9vvw2VK4dOIyIihUEjZNFav96PjiXQtkmjGMUGNmi6MoF8+y3cfz9ccQXcdlvoNCIiUlhUkEUrqyBLIH3ow+EcTh3qhI4iUfjrL79Z+KGH+m78CVT7i4jIXtKUZbQSrCDbxCZGM5o7uEMbiSeIjh1hyRKYOBEOPjh0GhERKUwaIYvWhg0JVZANYhBppHEjN4aOIlH45BM/Knb//XDhhaHTiIhIYVNBFg3nEm6ErD/9OZmTqUa10FEkFz//DLfeCjVqwBNPhE4jIiIhqCCLxtatsH17whRkq1jFFKZwIzdiaCFSPMvIgJtu8scPP4RixUInEhGRELS4KBoJ1oNsAAMA1HssAXTpAl9+Cf36QaVKodOIiEgoGiGLRoIVZP3pzzmcw3EcFzqK5GDuXHjkEWjc2I+SiYhI6lJBFo0EKsjmMY8FLNBi/ji3Ywc0beqvpnz7bbW4EBFJdZqyjEYCFWT96U8RinAt14aOIjl45BGYP99vIK4WFyIiohGyaCTIxuKZZDKAAdSmNuWI76yp7Msv/dqxNm2gbt3QaUREJB6oIIvG+vV+D8sSJUInydEUpvAzP2u6Mo5t3gzNmvkF/F26hE4jIiLxQlOW0di4EQ6J/825BzCA0pSmEY1CR5E96NjR9x2bNg1Klw6dRkRE4oVGyKKxcSMcdFDoFDlKJ53BDKYBDSiNftPHo+HDoXdveOABOPvs0GlERCSeqCCLxsaNcb/yehKT+I3fuIZrQkeR3Vi/3nfjr1YNHn00dBoREYk3KsiisWlT3BdkgxhEaUpTF60SjzfO+QX8mzfDBx9A8eKhE4mISLzRGrJoxPkIWTrpDGEIDWhAKUqFjiO7+PhjGDoUXnoJKlcOnUZEROKRRshyk54Of/wR1wVZ1nSleo/Fn19/hfbt/Zqxu+4KnUZEROKVCrLc/P67P8ZxQTaQgezHftShTugoko1zcPvtfm/63r2hSJHQiUREJF5pyjI3Gzf6Y5wWZOmkM5Shmq6MQwMH+qnKF1+Ek08OnUZEROKZRshyE+cFma6ujE+//grt2sFZZ0GnTqHTiIhIvFNBlps4L8g0XRl/nIO2bTVVKSIi0dOUZW7iuCDT1ZXxadAgGDIEnn8eTjkldBoREUkEGiHLTRwXZJOZzEY2aroyjqxf76cqa9SAu+8OnUZERBKFCrLcbNwIRYv6zcXjzDCGUYpS1KZ26CgS0b69bwDbu7f/thEREYmGfmXkJqsprFnoJP/D4RjGMGpTm33ZN3QcAUaM8NOVTz+tBrAiIpI3GiHLTZx26Z/NbFazmiu4InQUwY+KtW0Lp50G990XOo2IiCQajZDl5vff4aCDQqf4h2EMowhFqE/90FEEeOAB+OUXv5i/WLHQaUREJNFohCw3f/wBBxwQOsU/DGMY53M+BxN/o3epZto06N4dOnTwi/lFRETySgVZbv78E8qWDZ3ifyxlKQtZqOnKOLBjB9x6Kxx9tF87JiIikh+asszNH3/EXUE2nOEANKJR4CTy/POwaBGMHg377Rc6jYiIJCqNkOXEubgcIRvGMKpRjWM4JnSUlPb99/DMM3D99VC3bug0IiKSyOKuIDOzOmb2g5ktM7P7g4b5+2/IyIirguxXfuUrvtJ0ZWCZmX6qskwZePXV0GlERCTRxdWUpZkVAd4ELgVWA9+Y2Qjn3PdBAv35pz/G0aL+EYzA4VSQBfb22/DVV9CnDxx6aOg0IiKS6OKqIANqAMuccysAzOwjoBEQpCAbt3Apjat8Bb1KQp/5ISL8w3aqUYLZtOa00FFS2rx5cMkl0KxZ6CQiIpIM4q0gOwr4Odv91cBZ2R9gZrcBtwFUqFAhpmGKFStKiTLbYF+Dohkx/VrRKgEcyeHsT3ztHJBqLr8cunSJuw0cREQkQcVbQZYr51wPoAdA9erVXSy/1kW1zuW3r2P5FURERETib1H/GuDobPfLR86JiIiIJK14K8i+AU4ws2PNrDjQBBgROJOIiIhITMXVlKVzLt3M2gPjgCLAe865hYFjiYiIiMRUXBVkAM650cDo0DlERERECku8TVmKiIiIpBwVZCIiIiKBqSATERERCUwFmYiIiEhgKshEREREAlNBJiIiIhKYCjIRERGRwFSQiYiIiASmgkxEREQkMHPOhc6Qb2a2AVgV4y9zCPBbjL+G/JPe9zD0voeh9z0Mve9hpPL7foxzrtzuPpHQBVlhMLNZzrnqoXOkGr3vYeh9D0Pvexh638PQ+757mrIUERERCUwFmYiIiEhgKshy1yN0gBSl9z0Mve9h6H0PQ+97GHrfd0NryEREREQC0wiZiIiISGAqyHJgZnXM7AczW2Zm94fOkwrM7D0zW29mC0JnSSVmdrSZfWFm35vZQjO7M3SmVGBmJc1sppnNjbzvT4TOlCrMrIiZfWtmn4bOkkrMbKWZzTez78xsVug88URTlntgZkWAJcClwGrgG+B659z3QYMlOTM7H9gK9HPOnRo6T6owsyOAI5xzc8ysDDAbuELf77FlZgaUds5tNbNiwFTgTufc9MDRkp6ZdQKqA/s75+qHzpMqzGwlUN05l6p9yPZII2R7VgNY5pxb4ZzbCXwENAqcKek5574ENoXOkWqcc2udc3Mit7cAi4CjwqZKfs7bGrlbLPKhv5JjzMzKA/WAXqGziGRRQbZnRwE/Z7u/Gv2CkhRgZhWBasCMwFFSQmTq7DtgPfCZc07ve+y9CtwHZAbOkYocMN7MZpvZbaHDxBMVZCLyH2a2HzAY6Oic2xw6TypwzmU456oC5YEaZqap+hgys/rAeufc7NBZUlRN59y/gLpAu8gyFUEFWU7WAEdnu18+ck4kKUXWMA0G+jvnhoTOk2qcc38AXwB1AkdJducCDSNrmT4CLjKzD8JGSh3OuTWR43pgKH55kKCCLCffACeY2bFmVhxoAowInEkkJiKLy98FFjnnXgmdJ1WYWTkzOyByuxT+IqLFQUMlOefcA8658s65ivif6587524KHCslmFnpyEVDmFlp4DJAV9RHqCDbA+dcOtAeGIdf4DzQObcwbKrkZ2YDgK+Bk8xstZm1Cp0pRZwLNMWPFnwX+bg8dKgUcATwhZnNw/8R+JlzTm0YJFkdBkw1s7nATGCUc25s4ExxQ20vRERERALTCJmIiIhIYCrIRERERAJTQSYiIiISmAoyERERkcBUkImIiIgEpoJMJMWY2dbcH5Xn16xoZjfk8Llt2dppfGdmxc2shZl1M7OHsp3PyHa7g5mdZGaTIvcXmVmPXF7/ezN728z2MbOGZnZ/Pv89k8ys+i7nHjOz53Y5V9XMFuXwOo+b2T35ybCb1yplZpMjWy1VNDNnZk9n+/whZpZmZt3y8rXNrImZPRRlhuJm9qWZFc3/v0REdkcFmYgUhIrAbguyiOXOuarZPnZmfcI590zWeWBbtse8DrwOdI3cPwV4I6fXB04HKgNXOOdGOOeeL4B/W5YBwHW7nGsSOV8YWgJDnHMZkfs/4jfIznINkJ9eiXWBXHtBmVnRyP9vE/nn+yAie0kFmUiKMrNakZGgT8xssZn1j3Tsx8xWmtmLZjbfzGaa2fGR833MrHG218gabXseOC8ySnVXAcY8Aliddcc5Nz+nB0caOn8FHJ81AhfJOdzMmkVutzaz/pHbl5nZ12Y2x8wGRfby3NNrLwF+N7Ozsp2+FhhgZrea2TdmNtfMBpvZvrs+P/uoW2Q0a2XkdhEzeyny/Hlm1noPEW4Ehme7/zewKNtI3nXAwN183UpmNifb/ROy7kf+/64KzDGzC7KNTn5rZmUi3yNTzGwE8H3kJYZFsohIAVJBJpLaqgEd8aNKx+E79mf50zl3GtANeDWX17kfmBIZyeq6m89XyvbL/s085OsKfG5mY8zsLotsM7QnkULoYmDXwu024FEzOw+4G7jDzA4BHgYuiWx2PAvolEueAfhRMczsbGCTc24pfuTqTOfcGfidPfKyw0Qr/Ht9JnAmcKuZHbvLv6s4cJxzbuUuz/0IaGJmRwMZwC+7vrhzbjnwp5lVjZy6GegduV0NmOt8h/B7gHaRkcbzgG2Rx/wLuNM5d2Lk/oJIThEpQCrIRFLbTOfcaudcJvAdfuoxy4Bsx3P28utkn7JsF+2TnHO9gVOAQUAtYLqZldjNQyuZ2XfANPx2LGN2eZ1fgUfxm3ff7ZzbBJyNL0SnRZ7bHDgml0gfA43NbB/+d7ry1MhI0nz86FGVaP+N+P38mkUyzAAOBk7Y5TGHAH/s5rlj8ftfNolk25NewM1mVgQ/kvZh5HwdIOu9mga8YmYdgAMio43gv0d+zHqhyJTpTovsSSgiBUMLM0VS245stzP4358Jbje304n8IRcpSorHNB3gnPsFeA94z8wWAKcCs3d5WNYaspycBmwEjozcN/zekdfnIcvPZvYjcAFwNf8tVPvg163NNbMW+OJxV/9574CS2c4bcIdzblwOX3rbLs/JyrPTzGbjR/0qAw338PzBwGPA58Bs59zGyPnLIv8OnHPPm9ko4HJ8kVo78pi/dvN6JYDtOeQVkTzSCJmI7Ml12Y5fR26vBP4dud0QKBa5vQUo8BETM6tjZsUitw/Hjx6tycfr1MAvXq8G3BOZEpwOnJttfVxpMzsxh5fJMgA/lbrCOZe1vq0MsDaSdU/rq1by3/eucbbz44Dbs/07TzSz0tmf6Jz7HShiZv8oyoCXgc6RUb/dcs5tj3yd7kSmK82sLFA0qzgzs0rOufnOuRfwG52fvLvXMrODgd+cc2l7+noikncqyERkTw40s3nAnUDWQv2ewAVmNhc/OpQ1ejIPyIgsai/IRf2XAQsiX28ccK9zbl1eXiAyxdkTaBkZbbsbP+L2G9ACvyh/Hr7o3G0RsotB+CnJ7FdXPoKfbpwGLN7D87rgC69v8VOQWXrhF8zPiYwAvsPuZy/GAzV3PemcW+ic67ubxxflf0dA+wOZkdcBP9U5IdvnO5rZgsh7kcZ/pzJ3dSEwag+fE5F8Mr+WU0TkvyJXAFZ3zv0WOot4ZvYv4C7nXNMoHz8U6OmcGx25fw9Q1jn3SOR+L6CXc256HnMMAe6PXHUqIgVEa8hERBKAc26OmX1hZkWy9SLbrcjFBUuIjIZFirNKwEXZXu+WvGaIXO05TMWYSMHTCJmIiIhIYFpDJiIiIhKYCjIRERGRwFSQiYiIiASmgkxEREQkMBVkIiIiIoGpIBMREREJ7P8B63WkFcibQQwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# For display purposes, convert from x to the scientific units of the specific\n",
    "# input image used in this demo.  This includes accounting for adjustable\n",
    "# parameter values.  \n",
    "#\n",
    "# [WARNING: This will not produce accurate results if the parameters for \n",
    "# converting between the scientific values and rgb values (imgoffset, levelmin,\n",
    "# levelmax, and levelpower) are not updated to match the values used below.\n",
    "# in the script below.]\n",
    "imgoffset=0.22\n",
    "levelmin=-0.65\n",
    "levelmax=0.75\n",
    "levelpower=0.9\n",
    "scival=10**((levelmax-levelmin)*(x/1000)**(1/levelpower)+levelmin)-imgoffset\n",
    "\n",
    "# Plot the conversion curves.\n",
    "pp.figure(figsize=[10,6])\n",
    "pp.plot(scival,r,color='#FF0000')\n",
    "pp.plot(scival,g,color='#00FF00')\n",
    "pp.plot(scival,b,color='#0000FF')\n",
    "pp.xlabel('Input FITS Pixel Value (MJy/sr)')\n",
    "pp.ylabel('Output PNG Pixel Value')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abb5e11b",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a866898",
   "metadata": {},
   "source": [
    "Read in and prepare the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d3973fda",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the data.  Since this is a simple colourization exercise, the header is \n",
    "# not read into the python environment.  \n",
    "#\n",
    "# [MODIFICATION OPTION: Use a different FITS file as the input in the first line \n",
    "# below.]\n",
    "file=fits.open('NGC_3031_I_MIPS24_bgm2012.fits')\n",
    "img=file[0].data\n",
    "file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "22b3f6b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Crop the image.\n",
    "imgcrop=img[480:1600,180:900]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbc57d76",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3e101e9",
   "metadata": {},
   "source": [
    "Create and export a greyscale version of the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c471c4ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert the image to log values, which will reveal more detail when the image \n",
    "# is converted to png.  \n",
    "#\n",
    "# This specific image is background-subtracted, which means that, in areas that \n",
    "# do not cover the galaxy, the mean pixel value is 0.  Consequently, many of \n",
    "# the pixels are less than 0.  To make the colour scaling look more natural, \n",
    "# it is appropriate to add an offset to the data.  In this case, the offset is \n",
    "# 5 times the noise levels (or, in statistical terms, 5 times the standard\n",
    "# deviation in the background).  \n",
    "#\n",
    "# [MODIFICATION OPTION: Change the value of imgoffset to adjust the initial \n",
    "# background offset applied to the data.  For images that are not background-\n",
    "# subtracted, it would be appropriate to place this close to the mean \n",
    "# background level.  Otherwise, it may be appropriate to apply an offset\n",
    "# equal to 3-5 times the noise levels or to set it to 0.]\n",
    "imgoffset=0.22\n",
    "imglog=numpy.log10(imgcrop+imgoffset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "abbf9c40",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Range of imglog values:  -2.4564111  to  2.389815\n"
     ]
    }
   ],
   "source": [
    "# Print the range of values in imglog, which is useful to know\n",
    "# before rescaling the image.\n",
    "print('Range of imglog values: ',numpy.min(imglog),' to ',numpy.max(imglog))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7429d253",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\georg\\AppData\\Local\\Temp/ipykernel_14140/2561334731.py:15: RuntimeWarning: invalid value encountered in power\n",
      "  imglog=((imglog-levelmin)/(levelmax-levelmin))**levelpower*255.\n"
     ]
    }
   ],
   "source": [
    "# Rescale imglog so that the pixel values run from 0 to 256 between\n",
    "# minumum (levelmin) and maximum (levelmax) values.  The levelpower\n",
    "# value changes the exponent in the function describing this conversion.\n",
    "#\n",
    "# [MODIFICATION OPTION: Change the levelmin, levelmax, and levelpower\n",
    "# parameters to adjust the colour levels.  Usually, levelpower=1 works\n",
    "# for most images.  The levelmin and levelmax values should be set to\n",
    "# include at least part of the range of values in imglog given by the\n",
    "# step above.]\n",
    "#\n",
    "# [NOTE: Ignore the RuntimeWarning.]\n",
    "levelmax=0.75\n",
    "levelmin=-0.65\n",
    "levelpower=0.9\n",
    "imglog=((imglog-levelmin)/(levelmax-levelmin))**levelpower*255.\n",
    "imglog[numpy.isnan(imglog)]=0\n",
    "imglog[numpy.where(imglog<0)]=0\n",
    "imglog[numpy.where(imglog>255)]=255"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "18fca9dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write a greyscale version of the image to disk.\n",
    "#\n",
    "# [MODIFICATION OPTION: Use a different filename if desired.]\n",
    "imgsav=Image.fromarray(numpy.flip(imglog.astype('uint8'),0),mode='L')\n",
    "imgsav.save('ngc3031_spitzer24micron_bw.png','PNG')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c75e9895",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "201d4b87",
   "metadata": {},
   "source": [
    "Create and export a false colour version of the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4623c9bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Convert the image to log values (again).\n",
    "#\n",
    "# [MODIFICATION OPTION: Change the value of imgoffset to adjust the initial \n",
    "# background offset applied to the data.]\n",
    "imgoffset=0.22\n",
    "imglog=numpy.log10(imgcrop+imgoffset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "92dda614",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Range of imglog values:  -2.4564111  to  2.389815\n"
     ]
    }
   ],
   "source": [
    "# Print the range of values in imglog.\n",
    "print('Range of imglog values: ',numpy.min(imglog),' to ',numpy.max(imglog))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "cebf7635",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\georg\\AppData\\Local\\Temp/ipykernel_14140/2964144245.py:11: RuntimeWarning: invalid value encountered in power\n",
      "  imglog=((imglog-levelmin)/(levelmax-levelmin))**levelpower*1000.\n"
     ]
    }
   ],
   "source": [
    "# Rescale imglog so that the pixel values run from 0 to 1000 between\n",
    "# minumum (levelmin) and maximum (levelmax) values.\n",
    "#\n",
    "# [MODIFICATION OPTION: Change the levelmin, levelmax, and levelpower\n",
    "# parameters to adjust the colour levels.]\n",
    "#\n",
    "# [NOTE: Ignore the RuntimeWarning.]\n",
    "levelmax=0.75\n",
    "levelmin=-0.65\n",
    "levelpower=0.9\n",
    "imglog=((imglog-levelmin)/(levelmax-levelmin))**levelpower*1000.\n",
    "imglog[numpy.isnan(imglog)]=0\n",
    "imglog[numpy.where(imglog<0)]=0\n",
    "imglog[numpy.where(imglog>1000)]=1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "62ea8403",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create r, g, and b version of the image.\n",
    "imgr=r[imglog.astype(int)]\n",
    "imgg=g[imglog.astype(int)]\n",
    "imgb=b[imglog.astype(int)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "a606b4ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Insert the imgr, imgg, and imgb arrays into an image cube named imgrgb.\n",
    "imgrgb=numpy.zeros((1120,720,3))\n",
    "imgrgb[:,:,0]=numpy.flip(imgr,0)\n",
    "imgrgb[:,:,1]=numpy.flip(imgg,0)\n",
    "imgrgb[:,:,2]=numpy.flip(imgb,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8b859eef",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write the rgb image to disk.\n",
    "#\n",
    "# [MODIFICATION OPTION: Use a different filename if desired.]\n",
    "imgsav=Image.fromarray(imgrgb.astype('uint8'),mode='RGB')\n",
    "imgsav.save('ngc3031_spitzer24micron_falsecolour.png','PNG')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f720f0ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write the separate colour frames to disk.  This is primarily\n",
    "# for demonstration purposes and can be skipped.\n",
    "#\n",
    "# [MODIFICATION OPTION: Use different filenames if desired.]\n",
    "imgrgb=numpy.zeros((1120,720,3))\n",
    "imgrgb[:,:,0]=numpy.flip(imgr,0)\n",
    "imgsav=Image.fromarray(imgrgb.astype('uint8'),mode='RGB')\n",
    "imgsav.save('ngc3031_spitzer24micron_r.png','PNG')\n",
    "\n",
    "imgrgb=numpy.zeros((1120,720,3))\n",
    "imgrgb[:,:,1]=numpy.flip(imgg,0)\n",
    "imgsav=Image.fromarray(imgrgb.astype('uint8'),mode='RGB')\n",
    "imgsav.save('ngc3031_spitzer24micron_g.png','PNG')\n",
    "\n",
    "imgrgb=numpy.zeros((1120,720,3))\n",
    "imgrgb[:,:,2]=numpy.flip(imgb,0)\n",
    "imgsav=Image.fromarray(imgrgb.astype('uint8'),mode='RGB')\n",
    "imgsav.save('ngc3031_spitzer24micron_b.png','PNG')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae1e5ed1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}