generate_error_cube.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import fits\n",
    "from random import randint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ******** Make a mini-cube of the central 3x3 region ********"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Filename: /Users/orion/phd_research/3C297/0final_3C297_medsub.fits\n",
      "No.    Name      Ver    Type      Cards   Dimensions   Format\n",
      "  0  PRIMARY       1 PrimaryHDU     779   (66, 68, 2210)   float32   \n"
     ]
    }
   ],
   "source": [
    "fits_filename = ('/Users/orion/phd_research/3C297/finalcube_3C297_mediansubtracted.fits')\n",
    "hdul = fits.open(fits_filename)\n",
    "hdul.info() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cube = hdul[0].data\n",
    "cube_header = hdul[0].header\n",
    "cw = cube_header['CRVAL3']         #central wavelength\n",
    "mp = cube_header['CDELT3']         #microns per pixel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error_cube = data_cube[:,32:35,32:35]                 #3x3 cube: (32,32) to (34,34)\n",
    "#print (np.shape(error_cube))\n",
    "stderror_cube = np.std(error_cube,axis=(1,2))\n",
    "#print (np.shape(stderror_cube))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3x3 error cube\n",
    "summed = np.arange(len(data_cube), dtype=float)\n",
    "summed = 0\n",
    "for i in range(32,35):\n",
    "    for j in range(32,35):\n",
    "        print (i,j)\n",
    "        summed += data_cube[:,i,j]\n",
    "\n",
    "#print (np.shape(summed))\n",
    "\n",
    "y = summed             \n",
    "x = np.arange(len(y))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.64999997615814 0.000195000000530854\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([1.43452498, 1.43471998, 1.43491498, ..., 1.86488998, 1.86508498,\n",
       "       1.86527998])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print (cw,mp)\n",
    "start_w = cw - cube_header['CRPIX3']*mp \n",
    "x = start_w + x*mp\n",
    "\n",
    "\n",
    "# line-free region in spectrum => continuum \n",
    "select_cont = (x>1.70) & (x<1.80)\n",
    "\n",
    "#x[select_cont]     #check\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc8c2728d10>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2deZgU1dX/v6e32ZkZGHYGBxRQEBQBF8S4oeIeNSrmfaPx1RiTmJjVaEyMMSGa6C8u2dwSNUYlLknADRHcUFEEFWQEZN+XgZlhlp7p9f7+qLrVt6qrp2d6me7pOp/n4aGnqrrq3q6633vuOefeIiEEGIZhmMLElesCMAzDMNmDRZ5hGKaAYZFnGIYpYFjkGYZhChgWeYZhmALGk+sCqNTU1Ii6urpcF4NhGKZPsWLFiv1CiIF2+/JK5Ovq6rB8+fJcF4NhGKZPQURbE+1jdw3DMEwBwyLPMAxTwLDIMwzDFDAs8gzDMAVM2iJPRLVE9CYRrSGieiK6Ud/en4heJ6L1+v/V6ReXYRiG6QmZsOTDAH4khDgCwPEAvkNE4wHcDGCxEGIMgMX63wzDMEwvkrbICyF2CyE+1j+3AlgDYDiACwE8oR/2BIAvp3sthmEYpmdk1CdPRHUAJgP4EMBgIcRuQOsIAAxK8J3riGg5ES1vaGjIZHEYhmGyzr6WTiys35PrYiQkYyJPROUAXgDwfSFES3e/J4R4WAgxVQgxdeBA2wlbDMMwecsVj3yA655cgVAkmuui2JIRkSciLzSBf0oI8W99814iGqrvHwpgXyauxTAMk09sa/QDAKJ5+gKmTGTXEIC/AVgjhPiDsms+gKv0z1cBmJfutRiGYfINAgEAovlpyGdk7ZoTAXwNwGdE9Km+7WcA7gLwLBFdA2AbgEszcC2GYZi8JJKnlnzaIi+EeBfQu7J4Tk/3/AzDMH2BgnXXMAzDOBrdxI1GWeQZhmEKDunGiLDIMwzDFB6kq3y++uRZ5BmGYdIg37NrWOQZhmHSgC15hmEYB8CBV4ZhmAKEA68MwzAFDOn+GnbXMAzDFCDSkmd3DcMwTCHCgVeGYZjChX3yDJMGm/e3o9kfzHUxGCYh0ifPefIMkwKn3vMWzrl/Sa6LwTBJ4QXKGCZFdh3szHURGCYp7JNnGIYpQIhXoWQYhilcOPDKMAxTwPBkKIZhmAImNhkqp8VICIs8wzBMBmBLnmEYpgDhwCvDMExBo/vkWeQZpmeIPB3+MowKvzSEYVIkTw0jhjHBq1AyTIrk6zRxhlFxhCVPRH8non1EtFrZ1p+IXiei9fr/1Zm4FuMc8tXHyTAq5BCf/OMAZlm23QxgsRBiDIDF+t8M023y1DBiGFvydeSZEZEXQrwDoNGy+UIAT+ifnwDw5UxcK9c0+4N4a92+XBfDEeTr8Jdh7HDiZKjBQojdAKD/P8juICK6joiWE9HyhoaGLBYnM3zjH8vx9cc+wsGOUK6LUvDkq2XEMCqO8MmngxDiYSHEVCHE1IEDB+a6OEnZsK8NQP763wqJrrIV1u1p5Y6WyQucnF2zl4iGAoD+f0H5OOxyuMORKJZttnqtmFTpqs2cdd87+OojH/ReYRgmAU5eoGw+gKv0z1cBmJfFa/WYW/79GW54+uMef0/eRjsBemLpVlz20FK8yT77jJBotCQtpvpdLb1ZHIbpkoK25InoGQBLAYwjoh1EdA2AuwCcQUTrAZyh/503PLNsG15atbvH35OdtZ2/uKldexfpJ9ua0yobo5FoxmsgnKcRLsbR5KsL15OJkwghrkiw6/RMnD8fCdvc0KFVxQCAHY3+3i5OQZJo+BsIR3q5JLnji72tCIajOHJ4Za6LwiQgFnjNbTkSkRGRdxLSuozY3FE5KWIbi3xGSGQYdYacY8mfee87AIAtd52b45IwieBVKAsMeRvDNkmxoYi2raEt0IslKlwSNRonWfJM/mPMeHVg4LWgsfPJS5Hn/O7MkOh3dJIlz+Q/hruGLfkCQb+Pdj75kO7C6c7MtyeXbsEP/vVpBgtWeCRqNGzJM/mEk/PkCxLDXWPjk5eWfHfWQf/FvHr855OdmSxawcE+eaYv4OQ8+YKmK3eNnZXP9JxE7hq25Jl8Qhp19y1aj00NbTkuTTyOF/m2QNgQ555gJ+RB/TzZzuO+7MGluOLhwp/tyT55pi+gSsGqHQdzV5AEOF7kj/zla7juH8u7fbyRQmnnkw9r27JtaS7b0oilmw5g3qc7se1A4aZrJoptsCXP5BOqMVJRnH9Z6Y4XeQB4c133V7/sjk8+GI4m9Muv3dOCtkC4x2W048a5n+Lyh5dm5Fz5CFvyTF9AfUxLffkn8vlXoj5C1ymUmjvH6ybzd6ICs+5bguNH989YOVoKeCXGZD55ItvdDNOrqM+pQP7F49iS7yGiGymUgL1fXvrsP9rSlLHyVJX6MnaufCNRCqW05FnjmXxAfU7zMVeeRb6HyJ460sWMV0Bz2ViRIu9WTNB0c2ury7xpfT+fSfTTSEveVeCmfN3NLxuftx3wo6Wzd0dtzf4gWnv5mn0R9TnNx8w6FvkUsUvIUUXeLjgohd+l/OrpPhTVBWzJO9knb43pfOnuN/HlP7/Xq2U4+o7XceycxVk7fyQqujWnJF0C4QgeeWcTQpEoVu88iLqbX8bWA+0ZO79ah3ycEFWwIr+ruQP/WLol4+eV9zOZJR+wEaKQjSXf0+GdtVHk2l2zr7Uza+dOtnZNOCryslFFowL3vv4F9qexhpE/GG8kbGqICVMkKnrFNdARyk4mUzgSxaE/ewW/W7AuK+dXeXTJZsx5ZQ3mLtuG/+oTEF9dvSdj51eNkUA42usjrmQUrMh/9ZEPcNu8+rRfEXewI4RdzR1x2+3z5GPbTrnnrfj9hiUfE/lQD9/+a/X1F3lydwsXrN6NY+csxoebDqR1HiHsLbqESw0rHWgwhTkO2aZ+VwvuX7we35+b+rIVjfq7CRIxbc4inHjXG906VzActX2Gc0mn/hz//b3Ncfv2tXRmtAOTotsWiKCyRHNvJktY2HqgHTuakqcnf++ZT9Dkj53rxrmfYNLtC3tlhNJdClbkt+j54+n+2OfcvwTTlcYkz2afJ28WHOu1pch7FJG3W7J4V3NHF0FHs2WVS0tWBpBX7kj9JSmdoQhG3fIK/vzmhrh9iW6d6grLR5F36/d37Z7WlM+RTOQb24PY09K9UdRPX1iF6Xe9Effs5JKAXhbr89vYHsSxv12M3y9Ym/QcQogexQwEBPpJke8MdRnnOPnutzDjd28mPef8lbsAACeNqQEQS75INCFy7Z6WXp8VW7AiL0nX573TagEZ7prEKZQS642Wf7sVkbeWr6E1gOl3vYF7FtoPY63+6EytlyGE6LEIyM4qnd9YWlQPvbMpbl/CBcpCiQPcC1bvxpG/fC2ngiY7nnTcNclEvrv84r+rjTWSevKb2PmZWztD2HMwM+45aclb4y6N7dpv9lp9cnfKP5ZuxcTbF8a3USvKJaSn9GBHWItz/CkzcQ6f2yyl7Qnmwsy6bwlO+39vZ+Sa3aXgRd5OKFKx7sMWAbdNoYwKlHjdxt+qX3X93lac98d3AZizQqzr0ssH9qVVu7DnYCe2W15AIhvqzCMGo7Z/SbcF9qC/a4vn5c924/BfLMDGHlgZHn0egN1opLvIiWF2lk/CwKtqyVu+99tX1qItEMbuHohRZyjSLVHpLnaZVT0lUyL/5Adbjc9WA+G1+j3wB+3FSPXFt+r36PjfLsbxd6YeiF2weg/+9MZ6vSy6Ja/cYiEE3l2/H4C26NePn1vZZbD55c+013fazfp+9qPteMvmXcvSQGj2a7/vpv1dB2Bnd3OyoTdO5LvuUHvTmi94kbdblyYVf58Uo1gKpb27pqwoNr9MbUDqmhYmSz4iUL/roNHx7NZFPhwRuH1+PX703ErTNaTAXTR5OEq9HkQiWvCx7uaX8bd34/2bADDv05046o6F+Nzy4uudzR341Yv1CIQjeG+D5ld/Y033X0IuA8ihNCx5+bvaCWPCyVBdWPLyO3aB8UTct2g9vvnkCry/YX+3v9MVqjsp1SBcpkReRbXkV2xtwjefXIHfJwh8qrOyD/pDaPYH0W4TDO4OGxvacPFf3sP1/1yBexZ+gXAkapuYMH/lLtz+4ucAtDkQz6/YgU+3N3fbKOsIRnD7fC0Od9MLq/D1xz6KO0bem61Kx9DV+T/Y1Bi3bfXOg3Hf8VliY3uTJCSs38ci3yPW723FJX99Hyu2NuGjLY34Ym/MF2onxqm4OFo6dJHvcjJUFGVFMUu+Q2kUHmX2q2rJz1+5C+c+8C5e+UyzJKUlH4pEsb8tgLZOs6UlrbFirwtuFyEcFWjXO5Nfv/S5bdnvX7zedG7JxX95D4+9twUfb23GoIoiAD3zI8shd6KhaXew1g8A/vD6F1iwek/CtWtUS946EpL3xy47JRENrZqLYGuGXtuoCpiaEaPy0qpdeH+j1qlsO+DH1N+8js2KVZmNDA3VOt+wT7vPiToTv2KJNvmD+GJv6qJ049xP8LHycvvtTR2meyjZmED49rd13eFJ4X7qw614/P0tePidjbbHRSLCGDGqr+i0nr8r0X/ls90474/vGqMIiccyu/3SB5di2eb4DkJityxKtigIkQ+Eo1ixtQkH2gK49MGlxnsxAXsxthP+P7+5ARv2teGe19YZPkg1l9ba6OwCnqFI1HTzVKFRh3OqJf+J/vBvOdAOIQR2NWsWQEtHGG2BcJw1K62xYq8bHjchEo2ixUYoJUIIQ2isPtm9LZq4bWhoM+q3dk/M2t/e6MeKrYkfVGntNSdxBXVFq9JBhCNRXP/kCjyweD2u/+eKbmXXhCIC2w74samhDXU3v2x0ZMmGyyrVpVowrsnftZjsa+3E9+d+gqYkVrYaDLYK10F/COv2tOKGpz/BVx/5EADwn092Yn9bEP/6aLtxXHfXN+qJ61G9//KZkHW38vP/rjY+N3eETKOTno6E5TMdu3abbXzA6vKQbEhk9erFkL/VAf2++Nxu02EyGOoPRWyvu3xLo+EmAro2EOp3aSPydRZjyG5i3l/eik8mkNi9PjRbFMTaNfLhCNn0jgc7QmhqD6K6LJZPbif8d7+2Dne/pg1dz5owBBNHVOLku98y9ltTrhIta3D0yCpDaFTLSRV29XNbQDvvS6t24+7X1mHGYVqUPhiJoskfRL9icyOMibwLLtIsebt0sGhU4PU1ezGqpszY1qp0BkIIEGmW74a9rYbYbmxoQzQq4HIRTvq9ll2w5a5z0RGMoC0QxkDd4gdiFnxXaaq3z6/HzCMGw+dxYdWOZlx70mjTftWSb2gLYIHiG39WET3Tb6AIzqUPLkVbIBy3TlBPRhcypXWv4sdvaA1gX2snJgyrNLZ9658fY8XWJswcPxjnTRqW8HyqIK7ZbXaRXf7w0rjRkrQC1biPv5udVGcoihKfO/mBMD+P6/TRbquNgRCJCryruK6a/UGUKQtvdYQiKFfckgfaAhhQXmQ6xzWPf4SzjhyCy6bWxo0WNja0YcygirjrJspI2dbYjnsWrkOzP4jFPzolbr98huRzqI6mZXkBbWQdFQLVpV5T2uO3nvoYAPDGj07G6IHlxsjYDmnQdIYipvts51rsKm2VLfkeIv1hwUh8w7ht3mpc/bjZN5cs7TAUjeKm582+cGnpxlIo7Zct6FfsxX++PR2A2V2jToBSNN6wQqQYqI2rqT0UZzXJhlDkccPjIkSFvcj/a/l2fPPJFfjJ86vi6iDPI5/LDQ1txjk6Q9E4t05nKIKj71iI8/XAMaBljsz7VEsfO9hhb9nua+nE4+9vwdWPL8NlDy3Fb15eE2d5qhbrCXea874Xr7WPDwRCUZTpwia/b+3gZUNtC4Tx3PLteGHFDttzATDS8LY3xep93h+X4NwH3jWVd8VWLWW0MxRFOBLFc8u32z5LMk4wbnAFXv5sN77Y24oFq3djxdbGOIEPRaLGyKDJHzIC5F0Jjcrj72/pdhBPHQHJgH6jzejFau02+0MmAfYr9+zVz3Zjym8W4eNt2m/T2qnVYfHafbjp+VVxMbESrxv7WgK2M8IPKJ2BOsLrCEawYmsTNiZwfclnQD7D1vYg69MZiiAQiqLIY98pypHUxn3m68h7PO/TnXjqw20ANHePOlq0G1CpHejBjpDpPrEl30MMkbexBLbs98ctZJUsI2V/awDPLjeLgnSJxNaT17av3dOCJ97fgh+eMQ7N/iCqSr2GZaUO+9SbqlryXbkVgpGo4bKIRAX+sXSLkbZY7HVrPvmIMLlrZPlkI165vRk+twvhaNQQs3Akaipbsz+EgDcKj+7j37CvDbX9S439H29tQiAcxZ6WTrQHwigr8pgm+mxr9OsjA/MvvVIPNvcv8xmuoXV7WzGkX7ExUzeVZZc7wxFUFHu7DATK+l324FJ8rnegl0wZYTrmoD+EP7y+Ds8s267/DjGRkeX94bMr8fuvTDLNbWj2B/Hou5tx16taLvelU2tN55WCePExw3Hnq2tN7kMrP35upWHVvfDxDrzw8Q6s/fWshC6DjmDEFOT73YK1eHTJJqz4xRkJr2F8Vxc7IQR26B2anU/eOsu12R8yrZOulu3FVVpHf/+i9agpL8KLq3aZ2qFVcKtKvTjYEbJdmkJ1g6nZYNa5EC2dIcx+6APjvrZZRpT7WmOpq52hiGFsyXoVeeNtW7eL8JkeUL3iEfMLeULRKIpcbnygTPpbv7fNNFq0s+Qb24NGu7j4L++ZOik7r0O2KAhLXg7VgzY/XFsgjNaA5t9+/L3NuPu1tUkteXWYde2MUQDiH1Zpyc+6bwmeWbYdD729EaGIwMljB6LUqzUINbtGLZvqv9tmE+zzWLJvAOD1z/fgVy9+jl/Orweg+VI1n7zZkv/e3E9x/J2LTSOAgRVF6FfiRWtnGG+s3YvDbn3VNIGpPRDGwY4Qjh89AETxb7fZofweMitALffelgB22aQsfqZfY2hlibFt1n1LcMo9b+FAWwDPfrTd1l2QjEAoivIkL2doD4QRCEcMIXAR8KsX6/G8YtEfdcdCPLE0lmLoD0bQHggbU98BzV/+Wv0eU4fS5A8aWVB28RApcqqrLBHzPt1lCrgCwLVPLMcba/fhuFH98djXp5n2HXHbArzwsdkAUS3gnc0dOKDn51tHTR3BCP6xdAvW7G41OiIp8p/vasHTupXaYelgmvxBk3Cro4wt+7Xn4O0vGvDCxzviDC3VlVfmc6OyxIvmjpCtb1yWZWBFkUnY1RGIEAJ7DnYa9xWIuWtkAF0V+ZaOEPyKuyYQiqLY48aYQeWma888YhDW7G6x7VylILcFIhhQ5sMdF07Apv3tplntdnGKQDhqdCzWUYg1JTubZF3kiWgWEa0jog1EdHM2rlGkB1q6yk/ec7ATt7/4Of785saklrwqWMOrS+Bzu7C/LYh9LZ1GXu+yLU1YvGavcdy7G/bD6yZMPaTasORVi0idDUtJ1q4ZUllsfJYWwtKNB/S/tU6tutQHt8ul+eQVN8yLK3dhb0vANLweWFGEimIPWjpC+PfHmoBJq6Si2IO2QAQHO0IYUV2CMYPK8dGWRqzeGRP6HYqgy8ylKkvA7tv/XBEnKnJGpt3w//v/+hQ3vbAKD769EV43mTq2o0ZUois6LT5hO+59/QujDocPqUBUAI+9twU/fm4lbpu32vZ37whF8KsX6/H9f5mXI3hp5W5DQADNrSLTRq2xACBmyQ+rKonbJ3n9B1/CpfrI4nOL31667MqKPHFZGwDiUj1lGXY0+XHiXW/ghqc/MZVDsmxzI26bV49rntDcl6Nqygzr+ZwHluBn//kMQgjjuZ09rRYjqktwsMPirtGFcMv+9riyW5Eif//so/HhrTNRWeLF65/vNVxfQKwN7G8P4JyJQ/DkNceazqEKflsgHOfPtrpr1Ht16j1voVOx5APhCIq8Lsy97njMve5447jph9agyR8yZeZJpCDvau7A2MEVuEwfuanPUKL4984me798b65WmVWRJyI3gD8DOBvAeABXENH4TF9HDl+7CrbtVaaAJ8sOUH3SVaVeDKksxq7mDpx1X2zY/c4XDbjmidhrA3c1d6BfsRcetwulusjfNq/e6HhU32SyjIgh/WIiLx8G1Zc7qKIYLpcmjJolH1/vdabji1BRpFnysrHK4fDAiiK0B8Jo7gihssSLo2ur8O6G/cbELSDmqyaKTTw5oKed3Xv5URg7uBwrdxyMG+bL4JZdCpw6WanE68ZHt86MlVepvx2BcDTpa9bagxFc8ldtIsu0OvNLWv6xdCs+2hKfNeQPRmzzlzftbzPNXn36w21Ysl57m9j9i9bHTcaRv7HaWat8fXodxgyuwE/PPtzYNv3QAXHHlfrcJtee5L96LEQiZ1surNeMjmV63aQY//CMsQBg1Fn+9pNHVqE9GDHNYm0LhA1L/ozxg1Fd6kOzP2jyoe9rCeDOV9fgRX1K/5PXHIt/f3u6kYarIkV+RHUJyos8xtoxzykjKn8wjEhUcyGNqC7FuMHmoKya9dTsj49TSf+9HBWqi+a1ByNYrWfEdAQjCISjKPK4MKC8CMePjv3mowdqoy67d7S2dobx5zc3YGNDG4ZWFaPY6zbauCSRj/2Me9/BQpuJdlaRX7B6j+0zmQmybckfC2CDEGKTECIIYC6ACzN9ESnyXa1joT7Ir1hyXK2o7prKEi+GVhZj98EOU0TeSktn2HAhqLNepdiGutHry3Q2NVtCupbUYeSgflpjMrJrbOqtWkqau8aD1s6w4UKaqweZBpYXoSMUQTAcRWWpF2MHx2c9bGv0o6bch8EVxdjT0olQJIrdBzvwvdMOw0WTR+DaGVrGzP2L1ptGU9LHbTe9X5WuB782xTSJbHC/eLGQtHZqVmUyS15l2qj4N3E9uiR+GYWOYMR0b774zdm4dsYobGv0Y1+LuQ7bG7Vn5EB7EOf9cQl+9WI9fv3S56i7+WUEQprfPNEy0LLDrykvwiEDtNjH2ROHxh0XFSJhWqGKz+NCNCoMQXW7CEII417LzDJ1dmdNuQ9TDqkGANMs1mZ/yOisS7xuVOmZKKrL5IHF6/HQ25vwwBvrUeRx4aQxA3HMyGocNzq+o5JlkuJeURyfsjnx9oVYsr4BwXAUo2vKQEQmQ0e1zJv9oThB3XqgHZGoMMRexlNke5JtZ2dzB97feADFSvu87/Kj8ftLJmG4PupSR7CSh97ZiLtfW4dmfwjDdNejHHkeObwfALOP/bGrp+Gpa48z/r7uyRVx57S6a377yho8qbgOM0m2RX44ADUPboe+zYCIriOi5US0vKGh++9aVXG7CC7qOoinLuZ056tdL35kFfnhVSVxub52SOtSXWVSNhjVXZMo/1sO79V1MOSxqh9UXegsGrVfpEkd4g6qKEZliRdN/mDc8FFNiaws8RoWjcr2Rj8GlBVhSGUxtjX68Zc3NyIqgMG6pSpdNw+9swkPvh2biGLtFH/95SNx3iRNzFRf/KQRVaZg4uCKxJb8lF8vAgCTyJ+giMtdF080PtcNKMWvLpiAkUoQ+bGrp6GiyIMl6+Nnt3aEIqZ77/O4MHJAKTpDUSN7xI6WzjAee2+LMeO4oS2AIo/L1goHgHFDYh3plJGa0B4zsgrPX3+CuTzBSMJzqDT5Qxj9s1fwnu7GCYajaGgLGG4yKbAqx40egAFl8Z3Qd5/5BPfoqcTFug/9YEfI9DzJeEwoIkznPs6mM5UiLxcGs5sEBQB3vqK1SRnHGDkgds/UkWCTP2hY8vdefhR+ef54rN7ZYmu4Pf+t6fjTVycb55QjGHXE/+XJw3HZtFqj7ckRxlUnHGIco05qUp8lrbyab181boZWFscdZ0XtFGScYWiCkV+6ZFvk7Z5Qk8IJIR4WQkwVQkwdOHBgyhfyeVxdBvF6stTqXsVqK/K4MbSqOPkiSAAqimIP/DdP1qxbKcDJ3hoFALeecwRmT6vFV48baWyT68L4AxFjpTtp3bndhHBUy5RRRc86lDx0UBmGVpZgz8FOk1UE2Ih8jTkgBWiBrAHlPgzpV4xlmxtx76IvAAADyrTv9lfEQl2nu6k9aCrL0SOqcOaEIQDMna7VKh9s466RLi4pNmrgdVpdtfH50qm1OGP8YFx5wiF488en4KrpdSYxG1xRjEm1lSYf8wmjB+BHZ4xFJCqwrzWA6lIvfn/JJACxRv23dzdjwrB++PNXj4krm5UdTR0J0/Tm33AivnZ8TEDOnDAYw6tKMLqmHEcON8ciOkIReF2Jm+g3Tx6NEdUxv/9yZfS27YDfsGBLFcv1jPGDAQD/d2Kd7Ujj0+3NxnlKfW7FXRP7vVS3nBqbueSYEThlnLkNy3aTbInfTfs1N9ko3chQ3VfmeEjQcHUMrijG7Gkj9aB6/GzvMp8H500ahjd/fAr++j+x+2Y3e1e17gGzy1A9vk7vMGRcbaA+P0DtBL1uF2r7l+KRK6fa1hXQXMZPf7gNC1bvQWN7EMFINKF7L12yLfI7AKj5ZSMA7EpwbFr43K64qcYqqa4VMaSy2BjWAsDJYxN3RKqf+HI9ONNiI/J2OcIAcEhNGe66ZJLpRSCqJX/YoHI8euVU/PGKyQBg+OT9wYjJcrBaBCeMHoBhVcVoDYTRHoyYXoCtinxVic8kGipHDO0XZxEOKNfKqU40269kdjT5g8bkLgA4bFA5amysRyt2AUtrGEUNvqkTcdwuwiNXTsUdFx5pNET19xhYURQnpl6Py+QuuvH0Mbhsmnb/jhjaz9h+6ZQROHfSUPziPHNYqX+ZzzSa2NnUYazzf+mUEZild2yANmpRA++zjhyK924+DSU+d5zQJLPkpx3SH18+erjtvv1twZjIKx3t7y+ZhH9ddzymHNLf1DnbId01BztC6AhGUOZzx708vaokdo4SnxvfPe0w0/6t+/0o9rqMTu/QgTEjQnVrhiLa5Lwa3XD47mlj8OiVUzG6pswk8tsO+A0XpttFKPG5ERWx507+XhXFHtOrMacfGnsOE4341fabyNVWV2O20GX7Udu37JitHZ5KKBrFz/7zGa7/5wpjhKFmoWWSbIv8RwDGENEoIvIBmA1gfjYuJJDY113mc3e5jkgRohgAAB5ESURBVIRKsZJDu/znM1FTXoTTDh9sDEXVtDjpl5ONRfU3yuHpbf+tx+Q7FposoUQz+6r076gvAgnrr0jzByMo83kwc/xgI4ddrl3jD4ZNlq1Hsf5q+5dgQHkRhlfFHs47LjzS+Dyw3GzJe9wuLP/5TEwYFhM3ALhx5pg4S0NayP2VBtHQGsDMP7yNUbe8gqgwBz1LfG4cPbLK1n2gMtymo7HmIauB9GRi5XG78KsLJuDwIRXoX+aLc1MIIUxCqHay6qhihj6SOlO3hiXzvnMiTj9ikPH3zuaYyN996VH4i25FJiunlSmH9LfN3pH4PC5cPq0W3zhplOHukp19S2fIGNnK57KyxIvqMp/hO6/uhshXlngNES0v9hijN0mlJcuqssR8zgX1ezBIcb/dfPbhhkuwpsJ8bEWRx3B1ul2EmeO1mdLqyGH9vjbDkpeZR8ePjj1jh+uusO+dNsY0mqos9WLLXefi5LEDMeei2POv8tr3v2R8TvQyHtleZGc3yE7kPdpOr9uFey8/yvY86kxv6WXIlrsmq5OhhBBhIroBwGsA3AD+LoSoz8a1ukqf7F/uQ3tj99w1Q/oVY8sBP4jM4iWHb6prYdEPT8Y76xvwzLLtaGwPmiwB+VkGg97fGJtIkWhdbyk0qn86GtUWVYpEBUot07VVS75/mQ9/vGIyqkt9uG2etu7I41dPwynjNPEZVhV7gNT1SqzuGkALCNbVlKFeWbWyX7EX1598KGYeMRjn/0nLvJENvp9FtNW1Rop9bnz1uJGo0H+3Up8HL94wA1+6W1sy4fbz45Ot7LI0rCJ/6dRaLPxcyyaRI4quuGp6Ha6aXgcAKC+K72TUYLdVuI4c3g+rd7YYVmht/1JcNHm4sU57ZanXNHkMMN9Dl4vwt6umYryl47Tj5e/NQP8yH9oDERwyoDRuqWnrNWr7l+LWc8dj8dp92NTQjtr+JdjW6EdLR8gwQvqVePDezafFucWqlPt24dHDjBnMkmLdXQNoLrsij9t4dg8dWIaNDe1xHXad4kufVleNj7Y0mYyDYq8bJx5ag00N7aalErRyxt8XVWzHD+2H9fvaDJ+8Wzdm/vo/UzD5168DAOZcNBFD+hUndH088X/H2m4HtBHk0984Diu2NMFqL86eVhs3CgNiBkl/pfNTjayLJo/AD/5lnj0PwJg5CwAf6gZoX/XJQwjxihBirBDiUCHEnGxdp6s3BA3tF7MMH7t6WsLjgFjaW5nPYwqgymn0qpDX1ZThyhPqUF4Uv8/qk12xtQklXjfOmTgkoSUvHyA18BrWRVyWSSVmyWtW/vlHDcOMMTXGb6G6XlTrWB2KqpaZKm4jbFwmJT43Jio57P1KPEY57Dj/qGG45Jjh+O1FE3HLOUcY22v7l+CbJ4/GS9+dga+fOCrue2U2mTNCxNxcPzlrnMkNZLUuk2GXflmq/LZVFrF55hvH48OfnW5q4PdefrTxudznibPSrZbg6UcM7tZwfMKwSgytLMFhg8rhdbu6zK5RryEtzKGVJSDSfN/S/11RrCUPWAXZo5/75LEDcf/sybji2Fqcf1RsTR7prgG07DSfx2WI36m68WD9rTxuF1bdfiZevfEkY90fq3jJDrW8yIPlP5+JK47VYlB2Izy1szxyeD9sPdAes+T15666zGdkx5T53Gn5tqcfWoPvnj4mLs35gqOGmWJlPzlrHADgmJHVuOfSo3Cf8jxYXyCSjLe/aIDHRXHr/2SKgljWAEjsqgGA284fj22Nfry9rsE2A0BFpm71swiBbFB2giYXEUuWB+F1E9wuV1xZz500FB8q61Zbp13LoZ01oOp2ERpaA2hoDeDo2ipj+wNXTMYj72xC3YCYa0l1y1SVevHSd2fg/Y37Mbgytr1CEVe1g7C6biSq6L14wwx43ISz71+Cn5w1DoFQBDecNiZunW35vVvOPiJue1dEhUCHvgRERbHH5MZQfa/dwTpbVgjzb2v1x1YUe21T/yQuF2FybRWuPrEOc5dtR0coYmuVpoKrC5+8+hjJEVl1qRcVRR60dIaNDqKrOQUrf3mm4Ru/8+JJeH/DfiP/3et2GSK/r7UTYwdX4ITRAzB/5S6cMm4QHn13c9ykOEBrD/2GxjK1rDUo1p+J0iIPasqLMNoIZsaXT31+qst8CEWixmxztS3+4bKj8NtX1sSNqFLFKifWuRv/c9wh+J/jtAD6VyzLZVgnsN0/+2i4iLDnYCfmvLIm7lob9rVheFVJtzKpUqFgRL4rZLDtHJtcZCtDdGvL2khlUCwYiWL2tFocpkyLvv6UQ7F47b6ED9ghA0qx9YAfPo/LNLNTYs3YsFoCMnhrtXDVYWGZIlLHjKzGX/93iulYVZCrSjXLRwYgf3HeeHy0udEkKNLy//r0Ovxg5ljTua6dMcpYxVAiLfwNc842LMSesuSmUxOuox6JCsPHXF7kMTUI60qdyaiw/I7WoKedcCXD43bhl+dPwK7mDrxWvxfjhyZ3zXQH9XH56azDcc/CdYa7QnX7XXzMcCxasxdHDq/Egvo9eOeLBmza3w6fp+vRQJy7xbIUg4xPhCICRR4XfnfJJFwzYxTGDanA+KH9cJRiXFiR37VOkivSf2ufLobSNx8Kx1tqsi0UebTgbTgq4ix5QEsJnXfDjIRl6SnHjDTXyy7dNBHW3/tCJTh+1fQ6HPWrhdrs6gsm4M9vbsC+1kDWMmsAh4i8dehc7HXZLpAExIaWVqtZniMQiuIuPb1OMq2uP5b97HTUJBhuDe5XjK0H/PC6XbbrTluxWr9S+OwseUmJr/u30rqG+DUzRuGaGWa3ycThVTh8SAWuml4X56P++XmJJy2nKvAATJ3kY1dPw9XKm32iIpYVUV7kMXVa1qyUZKhW+Q2nHoarptfhQHvAdn9XPH71tLjZvLLjtYplqpBiBxMBr954Ev70xgbMX7nLNDP0tMMHY+2vzwYAPPzOJiOe0tNXEQ6rKsGD/zvFeA2k6o7xeVwo8bkNYX/lxpO6PNcp4wbimJFV+NGZ40zb5f2Sz2+NTRqiRLo9S31uuIkgRCzI2dUoJ10OG1SBl747w5j5nWyGtUqyYLnspIu9LoyqKWORzwRWERhQVpQw710G/axW89UnjsIXe9uM4F3c92xyu+/+itYZLNLXuPFY1mgBzJaaJE7k9WUL4i352JetHYAdL94wA6+u3m1KXUvEwIoiLFCyDXob6fOVCCHiskVSRXXX/Fj3rVaVejGwoghuom4Pm0+xlBEAzps0FC9/ttuUdpsOalEIwNjBFXjgisl4QE+jtaOnIxsrs46MpXyqI9pEuf9dlePf3z4xbrvMYJMdoiHyNh2SbAulvtg6PupkwGyitsOeGC/W4KyVUFS+3U3LNvtwc2Na70lOhiNE3ur+qCn3JRR5mY9tFc3qMh8e/NoUm28kRi5BK5cY8LpdhvVRU16E/W0BW5+1tbxyQpVVnFUxslvIysrEEZWmwGlfIipiv0NPrCo77JZE8Lpd+OCW0xO+2Lq7nD1xKNbcMavbL/JIiiry3dS0ZCmqPcHr1mbuRqIiYVphTyn2mC15GbS2uy+yfZQVxdbxkYkL2fJhS7LViciYXJHHjeu/dChe/WwPLjrGfr5DJnCEyFuHdYmi2LecfThq+2u+6JPHxltpqSLzkb2umE9+Wl01QhGBq0+sizueiOBzu4zhq5yGbRUO1fVjXR620IiK2Lr5diJ/xNB+GNm/e5NJEq1743ZR2qMEIP4+pYN6j7vj6gO6Xv0yFXxuFzqiEcOXni7SZpVtYVBFEX4663CcM3FI3LGyYynxeYzj5ZLEni5mA2eC7qwblA7FXheqy3x456ZTs3qdglhPXuWV73XtJwSA608+1Hb75Xou7JKbTsUVx9baHpMK0gcejEQN66O8yINHr5qKE5VUQBXVwr9PfxG31cJXJ2D05J2mfYWJyszUqBDGsrh2k4pevfEkPPS1xNPIVbJtAWaSVEp64mGx2beJZjD3BPks9jQ1MBHW7BgiwrdOORSHDIiPYxiWvM9tdHL3LVpv+n62yL7IZ84Y6IqCE/nuTDg5dlR/bLnr3Ljt0uKv7V+a1K/WE2Qe94G2gPFg2uWCq6giL98raU2tVPPtxwyOX3Omr/PPa4/DbH15ASG0V9V53dSjFSgTcePpY7qcGJMvqNZ7d59JOaP1O6cemjQ42h3ks2j3RqVUCCvLEnT32qpPXpJtn3xXAdRM0FsiXzDumiH9ik2LXqVCth6aCfpypC2dYeMayQKl50wcgo372rFUeeVYkeUt9HJy0JUnHGKIYSEh17ef+9F2w5KvLvVlpAP+wRljkx+UB6hV7W6ty4s82PjbczJm6UqXSaZ88rKTtpvZbEUGkVs7Q3H1cWdZhL0Zqm8iijPUaSajYER+wfdPSulVcird9Xn2FPXN9DJ9KpnI/+bLE/Hc8u1mkbda8noaqN1060JB3pOo0HyxPV3/pa9DJku++9/LpCvDsOR7mF2TiPMnDYM/GMHF3Qg2XnDUMNz92jo0tgfjjDB3lp/5rlYAzQTFGfo9k1EwIl9V6jMtLJUK2bLk3S7ClEOqMbqmzHCxlHYjr906PLX6ROW5MmVh5SOyHUej2qqWiVYHLFRUHcuWEZIM+dzZZYKlgstFxlIGyajtX4oHrpiMcYMrsEp5LzHQGz55dtcUHNl8aF741nQAwE3Pa4sVdSev3dqorVlCcsZjIYu8vCdRIdDYHsThQ2Ixl7oBpagscNFX73iuB2u5es4u0NfT+Xy3+a1N2fbJu12E40f3x1Un1GXl/OyuSYOvTBmBDzYdwI4EL9FNRG+4PAzruxs3OFmnI8/VWxZBLlDdNU3+kGmdmrd+kt3Us3zAFHjNYTmA3BsTbov7JNuWPBFh7nUnJD9Q55JjRpiWwE4GW/JpcM+l2hrOdTe/nOOSxBOzvpPf4GSWSpm++mWyTJ2+jNS4cCSKZn/QtPyzE6A8MuUzlSefKtb2kG9xqP93mf3a8YnorU6zcNUhT4lZ38lvcDIf7J0XT8K0uh1xiykVEvI3aO4IISqSv+ii0MgrSz7LeePJyFVMIlv0VidV8CKvvsQgH5Drbvjc3bDkkwR++pf5cO1JozNSrnxFNuwD+uvdnJZdo5IrkTOm4feSDzkR2fbBFyoFLfKrbj8zY7P0MoXhrsmAJe8EZLuWqz06LbvGlWIKZSYR+kIEOffJZznbpbe48fQxeK1+T/IDM0RBi3y6q/FlA9lou7MSZLbX5ugLyCHtfoda8qlMhsoWmUqhTJVCseR/cMbYXp2MV9Ain4/cN/toPLNsW7deKsEaHxM5uSiVo33yOU+hzG3gNduTnwoVFvleZkR1KX5y1uHdOpYfasUnLxcnc5i7xpxck2OffK7dNQViyfc2bCvmMd1ZI77Qkb9Asz8In9uV0WV8+wL54K6RSwPn2pLn9pAaLPJ5DAdeYyIXDEdz7hPOBea1a3L7PHg9ub2+dTIU0z0c/au9/oMv4XeXTMx1MRLCgddYRxcMRx1vyeW69rk2Oth9mRqOVpExgytw+bTuLZSUC1jjYShbIBzN+ksc8p2cpVDqTvlcSyz75FMjrVZDRJcSUT0RRYloqmXfLUS0gYjWEdFZ6RXTmfBDrVjykSi8Dv89cmVJn33kUAC5T191+kguVdLNrlkN4GIAD6kbiWg8gNkAJgAYBmAREY0VQhTeO+p6iZtmjct1EXKCbNbBcDTrL3HId3Jlyf/wjLH4vxmjci7yuXYX9VXSajVCiDVCiHU2uy4EMFcIERBCbAawAUD+v2stzwhHtGHy4UMq8O1TDstxaXKDbNehSLRgJsP0NVwuyrnAA4UzGaq3yZZpNBzAduXvHfq2OIjoOiJaTkTLGxoaslScvklU94U6eZiqBl7ZJ+/c5wBg92WqJHXXENEiAENsdt0qhJiX6Gs224TNNgghHgbwMABMnTrV9hinEjFeeOxccTO5a5wu8rkuQI5xsrGTDklFXggxM4Xz7gCgvll6BIBdKZzH0UwYVomzjxzSZ146nQ1ICbw6vZE73SetplB+46RROSxJ3yJbyxrMB/A0Ef0BWuB1DIBlWbpW2rz70/x8w5DP48Jf/3dKrouRU2I+eZH1FyvnOw7XeJO75tZzx+ewJH2LtESeiC4C8EcAAwG8TESfCiHOEkLUE9GzAD4HEAbwnXzOrBlRnV9rzjMxVF3L9YzLXOPs2vPkwFRJS+SFEP8B8J8E++YAmJPO+RlGfXm50xu54y15h7vrUsXZrYbJe0yWvMMDr0635XlZg9Rweqth8hw1bdDrcEvO6RmEnEKZGizyTF6jGm8eh1vyTs+T58lQqeHsVsPkPWZ3jbMbubNrb47PMN2HRZ7Ja9TccE6hzHUJmL6Is1sNk/eowub4FEpnV59JERZ5Jq9RLXlOoWSVZ3oOv8ib6TOwT56ZVleNC462XeuQSQCLPJPXmHzynF2T6yLknOeun57rIvQ5nN1qmLyHUyhjsMQzqeDsVsPkParI+5zurnF29ZkUYZFn8hpT4NXhlrzTlxpmUsPZrYbJe1RZc/qMR2fXnkkVFnkmr1GDjT6Hv8ibVZ5JBYe3GibfMQVenZ4nzyrPpICzWw2T95jcNRx4ZZgewyLP5DVqsNHp64lz4JVJBRZ5Jq9Rdc3h3hq25JmUcHizYfId1Xp1uiXr7NozqcIiz/QZnD6t3+HVZ1KERZ7Ja0zuGseLnON/ACYFWOSZvIbdNTG4k2NSgUWeyWvYko/hdHcVkxos8kxeo1rvThc5Z9eeSZW0RJ6I7iaitUS0ioj+Q0RVyr5biGgDEa0jorPSLyrjRFRhc7q7xuHVZ1IkXUv+dQBHCiEmAfgCwC0AQETjAcwGMAHALAB/ISJ3mtdinAi7awx4WQMmFdISeSHEQiFEWP/zAwAj9M8XApgrhAgIITYD2ADg2HSuxTgTDrzGcHj1mRTJpE/+/wC8qn8eDmC7sm+Hvi0OIrqOiJYT0fKGhoYMFocpBFRdc7rIOb3+TGokfccrES0CMMRm161CiHn6MbcCCAN4Sn7N5nhhd34hxMMAHgaAqVOn2h7DOBe25GOwu4ZJhaQiL4SY2dV+IroKwHkAThdCSJHeAaBWOWwEgF2pFpJxLuYUSmeLnMOrz6RIutk1swD8FMAFQgi/sms+gNlEVEREowCMAbAsnWsxzkS1Xh0feHV4/ZnUSGrJJ+FPAIoAvK7nMH8ghLheCFFPRM8C+ByaG+c7QohImtdiHAgpZojT8+SdPpJhUiMtkRdCHNbFvjkA5qRzfoYx58nnrBh5gcOrz6QIz3hl8hpT4NXhKs+GPJMKLPJMXsNr16g4/gdgUoBFnslr1MCr033yDq8+kyIs8kxewymUMZxefyY1WOSZvEbVNadLnNPrz6QGizyT1/CM1xgOrz6TIizyTF7Da9fE4GUNmFRgkWfyGrbkYzi8+kyKsMgzeY0p8Orwp5VFnkkFhzcbJt8htuQNnJ5CyqQGizzTZ3D6ZCiHV59JERZ5ps/gdEvW4dVnUoRFnukzON1d4/T6M6nBIs/0GdhdwzA9h0We6TM43pJ1ePWZ1GCRZ/oMrPEO/wGYlGCRZ/oMTrfkHV59JkVY5Jk+g9NF3un1Z1KDRZ7pM3DglWF6Dos802fgPPlcl4Dpi7DIM30GtuQd/gMwKcEiz/QZHO+Tdnj1mdRgkWf6DE4XeaePZJjUYJFn+gzk8KfV6TEJJjXSajZE9GsiWkVEnxLRQiIapuy7hYg2ENE6Ijor/aIyTsfplryza8+kSrq20d1CiElCiKMBvATgNgAgovEAZgOYAGAWgL8QkTvNazEOx+nuCof3cUyKpCXyQogW5c8yAEL/fCGAuUKIgBBiM4ANAI5N51oMw5a8s+vPpIYn3RMQ0RwAVwI4COBUffNwAB8oh+3Qt9l9/zoA1wHAyJEj0y0OU8A4XOMdX38mNZJa8kS0iIhW2/y7EACEELcKIWoBPAXgBvk1m1MJm20QQjwshJgqhJg6cODAVOvBOADHW/LOrj6TIkkteSHEzG6e62kALwP4JTTLvVbZNwLArh6XjmEUHC/y7K5hUiDd7Joxyp8XAFirf54PYDYRFRHRKABjACxL51oMw4HXXJeA6Yuk65O/i4jGAYgC2ArgegAQQtQT0bMAPgcQBvAdIUQkzWsxDsfpeeLOrj2TKmmJvBDiki72zQEwJ53zMwwTw+nuKiY1HD6HkGH6DqzxTCqwyDNMH8Hp7iomNVjkGYZhChgWeYZhmAKGRZ5hGKaAYZFnmDxnWGVxrovA9GHSXruGYZjs8t/vnIg1e1pzXQymj8IizzB5zqB+xRjUj615JjXYXcMwDFPAsMgzDMMUMCzyDMMwBQyLPMMwTAHDIs8wDFPAsMgzDMMUMCzyDMMwBQyLPMMwTAHDIs8wDFPAsMgzDMMUMCzyDMMwBQyLPMMwTAHDIs8wDFPAsMgzDMMUMCzyDMMwBQyLPMMwTAGTEZEnoh8TkSCiGmXbLUS0gYjWEdFZmbgOwzAM0zPSfjMUEdUCOAPANmXbeACzAUwAMAzAIiIaK4SIpHs9hmEYpvtkwpK/F8BNAISy7UIAc4UQASHEZgAbABybgWsxDMMwPSAtS56ILgCwUwixkojUXcMBfKD8vUPfxjA95rGvT0NniAeBDJMKSUWeiBYBGGKz61YAPwNwpt3XbLYJm20gousAXAcAI0eOTFYcxoGcevigXBeBYfosSUVeCDHTbjsRTQQwCoC04kcA+JiIjoVmudcqh48AsCvB+R8G8DAATJ061bYjYBiGYVIjZZ+8EOIzIcQgIUSdEKIOmrAfI4TYA2A+gNlEVEREowCMAbAsIyVmGIZhuk3a2TV2CCHqiehZAJ8DCAP4DmfWMAzD9D4ZE3ndmlf/ngNgTqbOzzAMw/QcnvHKMAxTwLDIMwzDFDAs8gzDMAUMizzDMEwBQ0LkT2o6ETUA2JrrcqRADYD9uS5EL8N1dgZc577BIUKIgXY78krk+ypEtFwIMTXX5ehNuM7OgOvc92F3DcMwTAHDIs8wDFPAsMhnhodzXYAcwHV2BlznPg775BmGYQoYtuQZhmEKGBZ5hmGYAoZFvguI6O9EtI+IVndxzClE9CkR1RPR25Z9biL6hIheyn5pM0M6dSaiKiJ6nojWEtEaIjqhd0qdHmnW+Qf6ttVE9AwRFfdOqdMjWZ2J6Cd6fT/V6xYhov76vllEtI6INhDRzb1b8tRJtc5EVEtEb+rPdD0R3djbZU8LIQT/S/APwJcAHANgdYL9VdCWUx6p/z3Isv+HAJ4G8FKu69IbdQbwBIBr9c8+AFW5rk826wztlZabAZTofz8L4Ou5rk8m6mw59nwAb+if3QA2Ahit3+OVAMbnuj5ZrvNQaO/KAIAKAF/0lToLIdiS7wohxDsAGrs45KsA/i2E2KYfv0/uIKIRAM4F8GhWC5lhUq0zEfWD1oj+pm8PCiGas1zcjJDOfYa2XHcJEXkAlCLBG9DyjW7UWeUKAM/on48FsEEIsUkIEQQwF8CFWShixkm1zkKI3UKIj/XPrQDWoA+9s5pFPj3GAqgmoreIaAURXansuw/ATQCiuSla1khU59EAGgA8pruoHiWistwVM6PY1lkIsRPAPQC2AdgN4KAQYmEOy5lxiKgUwCwAL+ibhgPYrhyyA31I8LqDTZ3VfXUAJgP4sHdLlTos8unhATAFmsV+FoBfENFYIjoPwD4hxIqcli472NZZ334MgL8KISYDaAfQZ/y1SUh0n6uhWbGjAAwDUEZE/5u7YmaF8wG8J4SQFjDZHFNoedjWOgMAiKgcmvB/XwjRkpOSpUBWXv/nIHYA2C+EaAfQTkTvADgKmthdQETnACgG0I+I/imEKAQBSFTnJQB2CCGkhfM8CkfkE9UZADYLIRoAgIj+DWA6gH/mpphZYTZirhpA+y1qlb9HoI+4qHqAtc4gIi80gX9KCPHvnJQqRdiST495AE4iIo8+xDsOwBohxC1CiBFCeyXibGgBnEIQeCBxnfcA2E5E4/TjTocWrCwEbOsMzU1zPBGVEhFBq/OaHJYzoxBRJYCTodVf8hGAMUQ0ioh80J7v+bkoXzawq7N+b/8G7Tn/Q67KlipsyXcBET0D4BQANUS0A8AvAXgBQAjxoBBiDREtALAKmu/9USFEwjS8vkCadf4ugKf0xr8JwNW9Xf5USKfORPQ8gI+hvbD+E/SRKfHJ6qwfdhGAhfoIBvq+MBHdAOA1aJk2fxdC1Pdm2VMl1ToDOBHA1wB8RkSf6tt+JoR4pVcKnia8rAHDMEwBw+4ahmGYAoZFnmEYpoBhkWcYhilgWOQZhmEKGBZ5hmGYAoZFnmEYpoBhkWcYhilg/j+ZnWOsjjh/bwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cont_wave = x[select_cont]\n",
    "\n",
    "cont_data = y[select_cont]              #-- for 3x3 \n",
    "#cont_data = central_pix[select_cont]  # -- for per spaxel fitting\n",
    "\n",
    "plt.plot(cont_wave,cont_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/orion/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in sqrt\n",
      "  \n",
      "/Users/orion/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in sqrt\n",
      "  import sys\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(2210, 68, 66)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s = np.mean(cont_data)  # mean continuum => signal\n",
    "n = np.std(cont_data)   # std. dev. in continuum => noise\n",
    "\n",
    "snr = (s/n)\n",
    "\n",
    "sn_cube = np.sqrt(data_cube/s)*snr   # (datacube / mean continuum)^1/2  =>  error cube \n",
    "err = np.sqrt(data_cube*s)/snr   \n",
    "\n",
    "#check size & dimensions\n",
    "np.shape(sn_cube)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "hdu = fits.PrimaryHDU(sn_cube)\n",
    "hdul = fits.HDUList(fits.PrimaryHDU(sn_cube))\n",
    "hdul[0].header = cube_header\n",
    "hdul.writeto('3x3spaxel_error_cube.fits') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}