{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Radio Astronomy and Data Analysis Module\n", "\n", "**Lecturer:** Poonam Chandra
\n", "**Jupyter Notebook Authors:** Poonam Chandra, David Kaplan, & Cameron Hummels\n", "\n", "This is a Jupyter notebook lesson taken from the GROWTH Winter School 2018. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-astro-school-2018-resources.html\n", "\n", "## Objective\n", "Investigate radio observations, interferometry, and synthesis imaging\n", "\n", "## Key steps\n", "- Estimate the UV tracks for East-West, North-South or inclined array at the GMRT location\n", "- Calculate the UV tracks for different duration observations\n", "- Find the flux density for a bright source or a region covering an extended source\n", "\n", "## Required dependencies\n", "\n", "See GROWTH school webpage for detailed instructions on how to install these modules and packages. Nominally, you should be able to install the python modules with `pip install `. The external package `FriendlyVRI` can be installed by following the instructions at the associated website below\n", "\n", "### Python modules\n", "* python 3\n", "* astropy\n", "* numpy\n", "* matplotlib\n", "\n", "### External packages\n", "* FriendlyVRI - https://crpurcell.github.io/friendlyVRI/" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import os\n", "import webbrowser\n", "cwd = os.getcwd()\n", "data_dir = os.path.join(cwd, 'data')" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from astropy.io import fits\n", "from astropy.wcs import WCS" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "Radio astronomy reveals objects that do not radiate in other parts of the spectrum, and they are able to pass through galactic dust clouds that obscure the view in the optical range. In continuum, radio sky is dominated by non-thermal emission, i.e. the emission which cannot be characterized by a single temperature. The non-thermal radio emission usually are sites of particle acceleration and produce synchrotron, Coherent emission, Masers emission. They also reveal new phenomenon such as Fast Radio Bursts etc. \n", "\n", "In this image below, left hand side is the optical image of a galaxy which shows big sidelobes in the radio image in the right hand side. The radio image is dominated by synchrotron emission\n", "\n", "![SegmentLocal](data/Images/RadioGal.png \"segment\")\n", "\n", "\n", "Spectral line radiation is generated at specific frequencies by atomic and molecular processes. Neutral atomic hydrogen at 1420.405 MHz, which results from the transition between two energy levels of the atom, the separation of which is related to the spin vector of the electron in the magnetic field of the nucleus. \n", "\n", "In the left image of this galaxy, the optical emission is dominated by star light, whereas the extent of Neutral H1 emission, shown in the right side image, is much larger.\n", "\n", "![SegmentLocal](data/Images/NGC.png \"segment\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# # Radio Interferometry\n", "The angular resolution of a single radio antenna is insufficient for many astronomical purposes. For example, at 21cm, to obtain 1 arcsec resolution one needs 42km dish. Thus the concept is that for high resolution, use an antenna array instead of a single antenna. Signals from different antennas are brought to focus electronically. Now resolution is better since in resolution λ/d, d is the largest separation between the antennas in the array. Imagine making an image with mirrors with holes. More antennas means lesser holes. \n", "![SegmentLocal](data/Images/aper.png \"segment\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rotation of the Earth helps in filling the gaps between mirrors. This is called aperture synthesis.\n", "\n", "In the aniation below, one can see that Y-shaped array, and making use of Earth's rotation. \n", "\n", "If you want to play the animation again, double-click in this cell and click Run.\n", "![SegmentLocal](data/Images/ApertureSynthesisAnimation.gif \"segment\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![SegmentLocal](data/Images/radioast.jpg \"segment\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![SegmentLocal](data/Images/uv2.png \"segment\")\n", "![SegmentLocal](data/Images/uv1.png \"segment\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Hands On session on UV tracks\n", "\n", "In this code below, one can estimate the UV tracks for East-West, North-South or inclined array at the GMRT location. You can select differnt numbers of hours and check the UV tracks.\n", "\n", "Try different number of hours, and different array and see how UV tracks are different for different arrays. \n", "RUn the code below and compare 2 hrs and 10 hrs tracks for East-West, North-South and Inclined array.\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This will show the uv tracks for a given baseline for two observing session of different lengths\n", "Enter the numbers of hours (e.g., 5 10): 5\n", "Enter the wavelength in metre: 1\n", "Choose the array type: East-West (press 1), North-South (press 2) or inclined (press 3): 1\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAESCAYAAADaLCNlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4FNX+x/H3SSEJCaETwFBCSWiBQCCK1FAERUC6CogVRKTYrniVCyrcn3oVBQVEEVFAAUHpRVBCR2mh9x56KGkkIeX8/pglBEwgIbs7u8n39TzzsDszO/vZSdhvzpmZM0prjRBCCJFXLmYHEEIIkT9IQRFCCGEVUlCEEEJYhRQUIYQQViEFRQghhFVIQRFCCGEVUlCEEEJYhRQUIYQQVuGUBUUpNVgpdUAptVcp9YnZeYQQQoCb2QFySykVDnQG6mmtk5VSZe71mlKlSunKlSvn6X0TEhLw9vbO0zbsTTLbh7Nldra8IJntIau827Zti9Zal87xRrTWTjUBc4A2uXlNaGiozqvVq1fneRv2Jpntw9kyO1terSWzPWSVF9iqc/Fd64xdXoFAM6XUX0qpNUqpRmYHEkIIAUo74OCQSqlVQNksFr0LjAFWA0OARsBsoIq+44MopfoD/QH8/PxCZ82aladM8fHx+Pj45Gkb9iaZ7cPZMjtbXpDM9pBV3vDw8G1a64Y53khumjOOMAHLgfBMz48Cpe/2Gunych6S2facLa/WktkerNHl5XQH5YH5QDiwWikVCBQConO7kZSUFKKiokhKSsrR+kWLFmX//v25fRtT5dfMnp6e+Pv74+7ubqdUQoiccMaCMhWYqpTaA9wA+lkqaa5ERUVRpEgRKleujFLqnuvHxcVRpEiR3Kc1UX7MrLXm8uXLREVFERAQYMdkQjgZrSEH323W5HQFRWt9A+iT1+0kJSXluJgIx6GUomTJkly6dMnsKEI4Lq2hdWuoUwcGDIDate3yts54lpfVSDFxTvJzE+IetmyB1avhyy+NotK8Ofz0k1FobMjpWihCCCHu4euvScMFF9JRAOvWwY0b8PTTNn3bAt1CcUSVK1cmOjrX5xhkiIyMZOnSpVZMZD1nz56le/fuZscQIn+7ehVmzWIH9RnPEPZg6e56+WWbv7UUlHwkNTX1vgpKamqqjRLdrnz58sydO9cu7yVEgfXjj+jERLbSkKsUR6OgWDHo2dPmby0FRakcTUV8fXO87j+mLCQkJNChQwfq1atHnTp1mD17dsayL7/8kgYNGhAcHMyBAwcAuHLlCk888QR169bloYceYteuXQCMGjWKvn370qRJE/r27ct//vMfZs+eTUhICPPmzSMhIYHnn3+esLAw6tevz4IFCwCYNm0anTp1olWrVrRu3TpH2bZt20aLFi0IDQ2lXbt2nDt3DoDx48dTq1Yt6taty5NPPgnAmjVrCAkJISQkhPr16xMXF8eJEyeoU6cOYJwU8dxzzxEcHEz9+vVZvXo1ADNnzqRr1660b9+e6tWr869//csqP2YhCgSt4euvScWN8pylGNeoyX7o1w8KF7b528sxFJMsX76c8uXLs2TJEgBiYmIylpUqVYrt27czceJEPv30U6ZMmcLIkSOpX78+8+fP588//+SZZ54hMjISgH379rF+/Xq8vLyYNm0aW7du5auvviIuLo4xY8bQqlUrpk6dyrVr1wgLC6NNmzYAbN++nV27dlGiRIl7ZktJSWHw4MEsWLCA0qVLM3v2bN59912mTp3KRx99xPHjx/Hw8ODatWsAfPrpp0yYMIEmTZoQHx+Pp6fnbe8xYcIElFLs3r2bAwcO8Mgjj3Do0CHA6LbbsWMHHh4eBAUFMXjwYCpUqGCDn4IQ+czatXDgAO5AJxaRjsIFbZzpZQfSQjFJcHAwK1eu5O2332bdunUULVo0Y1nXrl0BCA0N5cSJEwCsX7+evn37AtCqVSsuX75MbGwsAJ06dcLLyyvL9/n999/56KOPCAkJoWXLliQlJXHq1CkA2rZt+49ikl22gwcPsmfPHtq2bUtISAijR48mKioKgLp169K7d29mzJiBm5vxN0qTJk14/fXXGT9+PNeuXcuYf9P69evp08c4+7tGjRpUqlQpo6C0bt2aokWL4unpSa1atTh58mTud7AQBdHXX9/21AUNLVtCzZp2eXspKCYJDAxk+/btBAcH89577/HBBx9kLPPw8ADA1dU1R8c37jZEttaaefPmERkZSWRkJKdOnaKm5Zcru9dllU1rTe3atTO2s3v3bn7//XcAlixZwqBBg9i+fTuNGjUiNTWV4cOHM2XKFBITE2nSpElG111O3Pz8udkHQhR4Fy/CvHnsoyZ7qE0qrsZ8OxyMv0kKiknOnj1L4cKF6dOnD2+99Rbbt2+/6/rNmjVj5syZAERERFCqVCl8fX3/sV6RIkWIi4vLeN6uXTu+/PLLm+OesWPHjvvKFhQUxKVLl9i0aRNgDF2zd+9e0tPTOX36NOHh4Xz88cfExMQQHx/P0aNHCQ4O5u2336ZRo0b/KCiZP8+hQ4c4deoUQUFB98wmhMjG99+jU1L4g9bMpTtHqAalS0OXLnaLIMdQcnihj7WHMdm9ezdvvfUWLi4uuLu7M2nSpLuuP2rUKJ5//nnq1q1L4cKF+eGHH7JcLzw8PKOLa9iwYYwYMYJhw4ZRt25d0tPTCQgIYPHixbnOVqhQIebOncuQIUOIiYkhNTWVYcOGERgYSJ8+fYiJiUFrzZAhQyhWrBgjRoxg9erVuLi4ULt2bR599NGMg/gAr7zyCgMHDiQ4OBg3NzemTZt2W8tECJEL6ekweTInqMxlSlKEOAI5BC+8DYUK2S9HbkaSdNYpq9GG9+3bl/WQm9mIjY3N1fqOID9nzu3Pz5byw6iyjk4y38Py5VqDnkN3PZKRejUttFZK66NHc7yJgnqDLSGEEJlNnEg83uynJgpNA7ZDu3ZQpYpdY0hBEUIIZ3b8OCxaxHYakI4LgRzClzi7Hoy/SQqKEEI4swkTSNOKLRh3Q3+Qv8DfHzp0sHsUKShCCOGs4uNhyhQ0imasowYHCOA4DBoEbvY/50rO8hJCCGc1YwbExOAGhLGFMLaApye89JIpcaSFIoQQzkhrGD/+n/P79IGSJe2fBykopvLx8TE7gk3JcPVC2NAff8D+/SynHatoTQKWwR8HDzYtknR5FUCpqan/GFvLFmS4eiFsaNw44vFmC41Ix4WGbDXG7apb17RI0kJxAFpr3nrrLerUqUNwcHDGcPERERG0bNmS7t27U6NGDXr37p0xhMrSpUupUaMGoaGhDBkyhMcffxzgtuHqmzZt6vDD1U+bNu224epHjBhhy10tRP5w5AgsWcJWGpKGK0EcpBgxMGSIqbGkhWIxalT2yzp2hMBA4/G2bbBo0f1tJzu//vorkZGR7Ny5k+joaBo1akTz5s0BY+ytvXv3Ur58eZo0acKGDRto2LAhAwYMYO3atQQEBPDUU09lbCvzcPWnT5+mdevWTjVcfWBgIG+88YYMVy/E3VhOFd5KQ8ByqnClSsaXlYmkheIA1q9fz1NPPYWrqyt+fn60aNGCLVu2ABAWFoa/vz8uLi6EhIRw4sQJDhw4QJUqVQgICAC4raBkHq6+Q4cOTjdcfVBQkAxXL8TdxMbC1KnspTbx+FCGi1TmhGmnCmcmLRSLe7Usbg7gGxpqTPaS26HctWW4+qCgoNsGtPzrr7/uOVz90qVLee+992jdujVdunShdu3aGaMLZ7ZkyRLWrl3LokWLGDNmDLt372b48OF06NCBpUuX0qRJE1asWPGPVoq1PqMQBdqUKejYWDbRGDBaJ8rLC154weRg0kJxCM2aNWP27NmkpaVx6dIl1q5dS1hYWLbrBwUFcezYsYybb2W+fbAMVy9EPpaSAl98wRke4Bzl8CaBuuyCZ56BLHof7E1aKA6gS5cubNq0iXr16qGU4pNPPqFs2bLZ3pTKy8uLiRMn0r59e7y9vWnUqFHGsszD1aemplK1alUZrl6I/GLOHDh9mgeA5/ieeHxwV2nw+utmJzPkZmhiZ53y4/D1cXFxWmut09PT9cCBA/XYsWP/sY6jZc4JGb7e9pwtr9aSWWutdXq61iEhWhuXNN6annjCKpuX4esLsG+//ZaQkBBq165NTEwMAwYMMDuSEMKW/vgDIiOJ545joW++aU6eLEiXl5N67bXXeO2118yOIYSwl08/JQ4fxjGUqhylF7NxeehBePhhs5NlKNAtFK1zdvtf4Vjk5yYKnF27YMUK/iaMVNxwIR0XNLz1FihldroMBbageHp6cvnyZflycjJaay5fvpzjU5KFyBfGjiWZQhn3PGnCBqhaFTp3NjnY7Qpsl5e/vz9RUVFcunQpR+snJSU53ZdYfs3s6emJv7+/nRIJYbIzZ+Cnn9hBA5LwpCKn8OcMvD4BXF3NTnebAltQ3N3dM640z4mIiAjq169vw0TWJ5mFyAc+/5y0lDQ28xAAD7PRGJ7+2WfNzZWFAtvlJYQQDu/yZfj6a/ZQh2sUoySXCeKgMcxK4cJmp/sHKShCCOGoxo2DhARi8cWVNJqxzhhm5dVXzU6WpQLb5SWEEA4tNha+/BKAZqynHjvxJgH6D4bSpU0OlzUpKEII4YgmTgTL7SEAfIkDd3eHupDxTk5XUJRSs4GbowoWA65prUNMjCSEENZ1/TqMHcspKhCLL7XYZ1x38uyz4MBnODpdQdFa97r5WCn1GRBjYhwhhLC+776DS5dYxXOcoiIdWEIjl+3w9ttmJ7srpysoNymlFNATaGV2FiGEsJobN+CTTzhJRU5RES8SjSHqn3zSuJjRgSlnvVJcKdUcGKu1bpjN8v5AfwA/P7/QWbNm5en94uPj8fHxydM27E0y24ezZXa2vFCwMpddsoQan37KDHpzhGq0JIKWrGHL1Kkk5OLaudzKKm94ePi27L5js5SboYntNQGrgD1ZTJ0zrTMJeCMn28tq+PrckuGz7UMy256z5dW6AGVOSdG6WjV9hnJ6JCP1GN7R1/G02hD1d2ON4esdsstLa93mbsuVUm5AV8CON+MVQggbmzkTjhxhLcah4kZswYsk+Pe/TQ6WM856YWMb4IDWOsrsIEIIYRUpKfDBB5yjLAeogTspNGYTtG0Lme7K6sgcsoWSA08CP5sdQgghrGb6dDh2jFK40Z7lpOKGDwkwapTZyXLMKQuK1vpZszMIIYTV3LgBH34IgDupPMRfxvx27RzqBlr34qxdXkIIkX/88AOcOMEN3G+f//775uS5T1JQhBDCTDduwOjRnMafsbzORhob8x97DB580NxsuSQFRQghzDR1Kpw6RQQtScKTRLyM+U7WOgEpKEIIYZ7kZBgzhlNU4ChV8SDZuIFWp07QMOfXEzoKKShCCGGWb7+FqChWEw7AQ2w2rjtxojO7MpOCIoQQZoiPh9GjOUEljhOAJ0nGdSdduoCT3gZbCooQQphh3Dj0hQuswhgYpDGb8CTZaVsnIAVFCCHs7/Jl+OSTjAPw3iQYrZNevaBuXZPD3T+nvLBRCCGc2n//C7GxFAZe4DtiKEohNw2jR5udLE+khSKEEPZ06hR89VXGUwUUIwZefBGqVTMvlxVIQRFCCHsaNYrUG2kspgOXKWHMK1wY/vMfc3NZgRQUIYSwl3374Icf+JswttKQeXRDAwwbBuXKmZ0uz6SgCCGEvfz73ySlu7OOZgCEsxpVvDi89ZbJwaxDCooQQtjDxo2wYAEbaEIiXlTiJNU4Ytw8q1gxs9NZhRQUIYSwNa3hjTeIw4fNPARAW1ai/P1h0CCTw1mPnDYshBC2NmcObN7Mn3QiBXdqsh9/zsCoKeDlZXY6q5EWihBC2FJSErz9NvF4s5tgXEmjDaugdm3o18/sdFYlLRQhhLClcePg5El8gIFM4jQVKMkV+OwncMtfX8H569MIIYQjuXgRxozJeFqSK0YxefRR4/a++Yx0eQkhhK385z+kxSVwgCDjehMAV1f49FMzU9mMFBQhhLAB7+PH4dtv2UpDZvEk83nCWNC/P9SqZW44G5GCIoQQNlB10iQS0wsRQUsAarEPfH2d8ta+OSUFRQghrG3ZMkps2cJampOIFwEcJ5BD8N57ULq02elsRgqKEEJYU3IyDB1KNCX5mzAUmnasQAUEwODBZqezKTnLSwghrOnzz9GHD7OMPqThSgO2U5YL8PEc8PQ0O51NSQtFCCGsJSoKPvyQw1TnKFXxJInW/AHNm0P37manszlpoQghhLW8+SZcv04VjtGaP/AmAW/XZPjyS1DK7HQ2JwVFCCGsYfVqmD0bADfSaMZ6Y/6gIU59n/jckC4vIYTIq5QUGDyYOHxIJNNxktKl8/VpwneSgiKEEHk1YQLs3ctiHudLBnOcysb8jz/ON/c6yQkpKEIIkRfnz8PIkRyiOgcJIg1XShFNTK1a+W404XuRgiKEEHnx5pukxiawnPYAtCSCIiqBw0OGgEvB+ootWJ9WCCGsaeVKmDmT9TTlCiUozSXC+Bv69yc+KMjsdHYnBUUIIe5HYiIMHEg0JVlHMwAeZzGuJYrdNmR9QSIFRQgh7seHH6KPHmUxj2dcEV+JU8bQ9CVLmp3OFE5XUJRSIUqpzUqpSKXUVqVUmNmZhBAFzJ498L//oYBW/EllTtCWldCiBTz7rNnpTON0BQX4BHhfax0C/MfyXAgh7CM9HQYMgNRUACpymmf5Aa9C6fD11wXiivjsOGNB0YCv5XFR4KyJWYQQBc2338LGjVyi1O3z33kHatQwJ5ODcMahV4YBK5RSn2IUxIdNziOEKCjOnYO33+YYAfzIMzRkK4+zBAIDjYJSwCmt9b3XsjOl1CqgbBaL3gVaA2u01vOUUj2B/lrrNllsoz/QH8DPzy901qxZecoUHx+Pj49PnrZhb5LZPpwts7PlBcfJXGvUKIqv2cAkBnKFErTmD5qxnsjPP+daSMht6zpK5pzKKm94ePg2rXXDHG9Ea+1UExDDrUKogNh7vSY0NFTn1erVq/O8DXuTzPbhbJmdLa/WDpL5l1+0Br2S1nokI/UEBupUXLR+7rksV3eIzLmQVV5gq87F97MzHkM5C7SwPG4FHDYxixCiIIiOhlde4Szl2MjDKDQdWYRrqRLwv/+Znc5hOOMxlJeAcUopNyAJS7eWEELYzJAhpF66wnz6k44LjdlEBaJg3MwCe81JVnJdUJRSwUAYxjEOT+AKcAjYqLW+at14/6S1Xg+E2vp9hBACgPnz4eef2URTLlKGklymFX9C587w1FNmp3MoOSooSqkqwECgN+AHpAPXgGSgGFAYSFdKrQGmALO11uk2SSyEEPZy5Qq8/DIAjdjCFUpQnx24Fy8CkyYV6GtOsnLPYyhKqSnAXiAE+ACoD3hqrUtrrf211j5AGaAjsBvjQsP9SqmmtosthBB2MGwYXLgAgCfJdGYhFTkN48ZBuXImh3M8OWmhJAI1tNYns1tBax0NLAOWKaVeB3oAD1gnohBCmGDxYpg+nWMEUJFTuJFmzO/QAfr0MTebg7pnQdFaD87NBi1dXbPvO5EQQpjtyhUYMIDz+DGDPpTkMv35Bvei3jB5snR1ZcMZTxsWQgjb0RoGDiTt7Hnm8wTpuBDAcdxJhc8/hwek8yU7VisoSqlictxECOH0Zs6EOXNYTTjnKUtxrtKGVdCuXYEeSTgn7ue0YR+gtmWqk+nfm0OluFotnRBC2NPJkzBoECepyAaaoNB05VcKFfM2BoWUrq67umdBsZwy/BJQF6N4VLi5CDgD7AFmWv7dbZuYQghhY2lp0K8fybFJ/EYXNIrmrDUuYJz4E1SocO9tFHA5aaEsAgKB9cAKjDO4AHprrZfZKpgQQtjV2LGwZg27aMg1ilGes7RgjXHxolzAmCM5KSg1gF5a67kASql3gbHAQqXUZGC41jrehhmFEMK2du6Ed98FoCFbKcQNynMWV//yMGGCyeGcR04OyncB1tx8orWO1lo/AzyCMZT8AaVUVxvlE0II20pKMq4rSUkBjL78euyiNNHwww9QvLi5+ZzIPQuK1nqh1vpSFvNXYxxX+QaYoZRaoJTyt0FGIYSwnbffRu/Zwwoeuf0ujK+/Dq1amZfLCeXptGGtdYrW+gOMwuIF7LNKKiGEsIcFC2D8eP7iQTbRmOn0JQ0XCA6GMWPMTud0rDJ8vdb6CPCIUupJa2xPCCFs7tQpeO45zlKOlbQF4DGW4lrIDWbMAE9PkwM6n5wMDtlXKZWja0u01rMsr6mmlGqW13BCCGETqanw9NMkXb3OL/QgDVfC+JsaHISPP4a6dc1O6JRy0uX1OnBUKfWhUqpedisppUoqpXorpRYBkYAMxSmEcEwjR6I3bGAxj3OV4pTlPI/wO3TsCEOHmp3OaeVkcMj6SqlewGDgXaVUPLAfiObW/VACgIrAVWAG8LLW+ozNUgshxP1atQr+7//YTgP2UIdC3KAHv+DmXw6+/16uhs+DHB1D0VrPBmYrpapinCocijHUijdwAVgLbAAitNYpNsoqhBB5c+GCcYqw1riShjspPM5iSrpcg5/my+188yi3B+V/Bt7VWn9jizBCCGEzaWnQt2/GDbNC2EkVjuFLHLz/ITSTw755ldvThvcAS5VSEUqpJrYIJIQQNvH++7ByJbEUyZjlS5xxrck775gYLP/IVUHRWj8P1MIYFHKNUmqZUqqBTZIJIYS1LF4MH37IdurzJYPZTR1jfunSxinCrjJIujXk+sJGrfVhrXVvoB5wHdiilJqnlKpp9XRCCJFXR45Anz6cpRxLeYwU3EnD1Tj4Pn263Bveiu77Snmt9V6tdTeMA/SFgF1KqelWSyaEEHl1/Tp068b1mBvMphepuNGQrYSwE0aONG6aJazmvguKUqqwpburNrATOA48ba1gQgiRJ1rDgAGk79rNPLoRQ1Ee4AztWQ6PPQYjRpidMN/J1VleSqnxGMPZ1wAewBiYMxU4glFUfrJ2QCGEuC+TJsGMGUQQzlGqUpjr9GQObgEVja4uF6vdAV1Y5Pa04TbAXmAqxkCQe4FDcu2JEMKhbNoEw4YRjzebeQiFpjtzKep5A+bNgxIlzE6YL+WqoGita9kqiBBCWEVUFHTtCikp+JDCC3zHGR6gCsfh62lQv77ZCfMtq4w2LIQQDuH6dXjiCfT589wcQMWPi/hxEV5+Gfr1MzVefiediEKI/EFreP550rdtZw49iSTTWLZhYfDFF+ZlKyCkhSKEyB/GjIHZs/mDNuynJieoTBAH8SpX3Dhu4uFhdsJ8T1ooQgjn9+uvMGIEkdRjA01wIZ2ezMHLQ8P8+eAvdye3B2mhCCGcW2Qk9O3LafxZREfAuPNiACfguxlGd5ewC2mhCCGc14UL0KkTMdfdmMWTGXdebMg2Y8DH3r3NTligSEERQjinhATjDounT/MbXUjAmyocM66E79wZRo82O2GBIwVFCOF80tLg6adhyxYAOrCE6hymB7/gUqe2XAlvEqfb40qpekqpTUqp3UqpRUopX7MzCSHsSGsYMgQWLsyYVZpoevMTXqV8jPlFitxlA8JWnK6gAFOA4VrrYOA34C2T8wgh7OnTT2HiRDbSmL9pdGu+h4dxRldAgHnZCjhnLCiBGPewB1gJdDMxixDCnmbPhn/9i93U4XceYSmPcYEyxr1NZsyAJnIjWTM5Y0HZC3S2PO4BVDAxixDCXtatg2ee4QSVmM8TADzC78awKp99Bt27mxxQKK212Rn+QSm1CiibxaJ3gYPAeKAksBAYorUumcU2+gP9Afz8/EJnzZqVp0zx8fH4+PjkaRv2Jpntw9kyO1teAPbto8nw4VyL8+A7XiAJTx7kL9qznDNdu3Lk1VeNVooDcbb9nFXe8PDwbVrrhjneiNbaaSeM7q+/77VeaGiozqvVq1fneRv2Jpntw9kyO1teffy4TipVSsfioz9nqB7JSP0zvXQaSusuXbROTTU7YZacbT9nlRfYqnPxnex0XV5KqTKWf12A94CvzU0khLCZ8+ehbVs8oqOZzxNcoxj+RNGNebg89KBx3MTV1eyUwsLpCgrwlFLqEHAAOAt8b3IeIYQtXL1q3PP9yBEAHmcxgRziaX7CvWol4/TgwoVNDikyc7qxvLTW44BxZucQQthQQgI8/jh6166M+5oU5xpP8zOUKwcrV0Lp0qZGFP/kjC0UIUR+duMGdOtG+sZNzKEn68l0KnCJEkYxkWtNHJIUFCGE40hLgz590CtWsJBO7Kcm62lKHD7g7Q1Ll0Lt2manFNlwui4vIUQ+lZYGzz2H/uUXfucRIgnBnRR6MxNv92RYsAwefNDslOIupIUihDBfejq8+CJ6+nRWE84mGuNKGr2YTQWXs+wbMQJatzY7pbgHaaEIIcyVng4vvQTTphFBS9bSHBfS6cY8qnEUpkwlWo6ZOAVpoQghzJOeDgMGwNSpJOHBTuqh0HTlV2qxH774Ap57zuyUIoekhSKEMEd6OrzyCkyZAoAnyTzLNM5RjpocgLFjYehQk0OK3JAWihDC/rSGV1+FyZM5jX/G7GLEGMXkf/+D114zMaC4H1JQhBD2lZZmdHNNmsR6mvAdL7CWZreWf/wxvPmmefnEfZMuLyGE/aSkwLPPon/6iXU0409aodD4Emss/+9/4V//MjejuG9SUIQQ9pGUBL16oRcuZBVt2EATFJrOLCCEnTB6NLzzjtkpRR5IQRFC2F58PDzxBPqPP1jKY2yhES6k05VfqcNe+OADePdds1OKPJKCIoSwrWvXoEMH2LiRP2nFFhrhRio9mUMgh427Lb7+utkphRVIQRFC2M6lS8YQ9Dt2ANCQrRwkiEdZRoA6CV9Phv79TQ4prEUKihDCNo4ehfbtST1yHFdAAUWJ5WW+xsXVBX6cAU8/bXZKYUVy2rAQwvq2bIHGjUk4cpZpPMsaWmQscinkDvPmSTHJh6SgCCGsa+lSaNmSK5dSmcrzROFPJCEk4mncYXHxYujc2eyUwgaky0sIYT3ffQcDBnA2rQwz6U0C3pTlPL2ZiVcxT6OYNGly7+0IpyQFRQiRd1rD++/D++9zmGr8Qg9uUIiqHKUnc/CoWBaWL4eaNc1OKmxICooQIm+/Bs0FAAAY9klEQVSSk2HgQPj+ew4QxBx6ko4L9dhJJxbiWi/Y6AYrX97spMLGpKAIIe7fxYvQtSts2ACAP1EUJYY67KEVf6Jat4ZffwVfX5ODCnuQgiKEuD+7dkHHjqScOosrChc0PiQwgMl4kgx9+hjHVAoVMjupsBM5y0sIkXvz58PDDxN76ipTeZ5VtMlY5EkyDB8OP/4oxaSAkYIihMg5rY0Rgbt0ISqhGN/Qn3OU4wA1SKYQuLnB5Mnwf/8HSpmdVtiZdHkJIXImIcG49/vPP7OLYBbSiVTcCOA4PfgFj5JFjAsWW7S497ZEviQFRQhxbwcPQrdupO/dx5+0Zj1NAWjEFtqzHNfaNWHRIggIMDmoMJN0eQkh7m7uXGjUCPbuZS3NWU9TXEinA0vowFJcH38MNm6UYiKkoAghspGSYtzXvUcPiIsD4EH+wp8o+jKdRmyFt982DtDLacEC6fISQmTlzBno2RO9cSN7qEMt9uFKOl4k8QLfoby84JvpxqnBQlhIQRFC3G7FCujbl+RLMSykO3upTRT+PMpyAFT16sbB9+Bgk4MKRyNdXkIIQ1KS0cXVvj0XL8G3vMReauNBMpU4aazTtasxNL0UE5EFaaEIIWDfPnjqKfSuXfxNGCtpSypulOEivZhNSdcY+PhT41a9cn2JyIYUFCEKMq1h0iR44w1SklKZw9McpjoADdhOe5ZTqFwpmL0amjUzOaxwdFJQ7iYmxjh7Rf4iE/nRpUvwwgvG9SMYXwaupOFFIp1YSE0OQJs2MH06lC1rblbhFKSgZCctDdq2Nf4jTZxodhohrGvuXBg0iJSLV0ikCL7EoYBOLCQVN3zdk+D/PjWOqbjIoVaRM1JQsjNunHHwESAigvIvvgjNm8t/LuHcLl6k1qhRsGYNZyjPAl7CjVRe4DtcSacwiVCjBvz0E9Svb3Za4WQc8ttRKdVDKbVXKZWulGp4x7J3lFJHlFIHlVLtbBLg2DF47z0ALlOCo3GlCfz8cwgPh0OHbPKWQtiU1vDzz1CrFiXWrGcVrZnCi1ykDMl4EEcRY70BA2DbNikm4r44agtlD9AVmJx5plKqFvAkUBsoD6xSSgVqrdOs9s5aG/+pEhPRwEI6cZJKhBBJu7Ur8KpbF0aOhDffBHd3q72tEDZz7pxxR8UFCziNPwt4kmhKodA0ZhOt+BP3kkXhu/nQubPZaYUTc8gWitZ6v9b6YBaLOgOztNbJWuvjwBEgzKpvrhS8/DL4+aFRVOMIbqQSSQhf8Sp7k6ui//1vCAsz/pITwlGlpRlDydeqBQsWsJqWTOV5oilFKaJ5ge9ox++4d+4Au3dLMRF5prTWZmfIllIqAnhTa73V8vwrYLPWeobl+XfAMq313Cxe2x/oD+Dn5xc6a9asXL23W1wcVSdNotyyZURTkkV05CSVAKjBATqwBB+X60R168aJfv1I8/bOwye1jfj4eHx8fMyOkSuS2TqKHDxI9S++wPfAgYx5f9OI5bTnYTbSkgi0rzeHhw7lYni4w5/J6Ij7+F6cLXNWecPDw7dprRtm85J/0lqbMgGrMLq27pw6Z1onAmiY6flXQJ9Mz78Dut/rvUJDQ/V9W7VK6ypVdDroLYTq/zJcj2Sk/h9v6Bu4aQ1aly2r9Q8/aJ2Wdv/vYwOrV682O0KuSeY8unxZ65df1lopHYuPPkh143cUdDroC5Q2nvfoofWFC2anzTGH2sc55GyZs8oLbNW5+F43rctLa91Ga10ni2nBXV52BqiQ6bm/ZZ7ttG4Nu3YR1aMHDV12MIgJBHGQh9iMO6nGOufPQ79+0LQpbN9u0zhCZCk9Hb7/HoKCSP96Mpt1GF/xKnPpTqzlgLsCypRR7Bk1CubMgTJlTI0s8h+HPIZyFwuBJ5VSHkqpAKA68LfN39Xbm6OvvAKbN+MbXJknmUUTNmQs3k595tOZ+E27oGFD4xjM5cs2jyUEYNyLpGlTeP55TkV7MZkBLKc9yXgQwPFb6/XrB/v2ES13VBQ24pAFRSnVRSkVBTQGliilVgBorfcCc4B9wHJgkLbmGV730qgRbN2KGj0a5ekJQBouRNCSSEL4ksFs0g+SNvlbqF4dxo41BtwTwhaOHIHu3aFJE2I37WE+nZnK81zAj+Jc5Wl+4ilm4VunEqxdC9OmQcmSZqcW+ZhDFhSt9W9aa3+ttYfW2k9r3S7TsjFa66pa6yCt9TK7hytUCN59Fw4cgG7dcCWdZ5lGIIdIxoMVtONrXubI1RLoN96AoCD48UfjjBshrCE6GoYOhZo1jWHkgUV0JJIQXEmjBWt4hYkE+pyDzz4zumFlHC5hBw5ZUJxCpUrG8BUrV1KiZlme5mee5idKcIVLlGYGffiRZ0g4FW10NYSEwJIlxuFRIe5HYiJ8/DFUrUr6+C+5nnrrOqhwVlOT/QxiAuFE4N6rm/FHz+uvy/VSwm6koORVmzawcyd89hmBRc7zChNpy0q8SCQBb7xINNbbswcefxxatoRNm0yNLJxMYiKMHw9Vq6KHD+dgbFkmMZDf6JKxSnnO0Ys5lKhdHlauhFmz4IEHTAwtCiIpKNbg7m78JXjoEG7PPUMTl80MYTzdmYsLRoskliIs5VHi126Dhx82Bp5cs8bk4MKhZS4kQ4dy/JwH03iWn3mKS5TmMiVJxDiWR7ly8O23EBlp/JEjhAmkoFhT2bIwdSrs3IlXx7aU4VLGogha8jdhjGMoK3iEuFWbjdZK8+bw++/SFSZuyaaQ/EA/TlIJLxJpz3IGMQEvb1d4/304fBhefBHcHHU0JVEQyG+fLdSpAwsXwrp18PbbsGkTjdnEdQpzgBpsojFbaEQDttNk3QaKtmtnDOXy3ntGt5iDX7UsbOTKFWOolPHjjWubgES8mElvUnHDi0Qas4kH+QsPl1SjgLz/vtyrRDgMKSi21KwZbNgACxZQ+p13ePLAbM5RlrU0Zz81+ZswthHKYywl9O+/oVMnY9ylIUOgb18oXNjsTyDs4cgR+OIL+P579PXrHKMKlXHJGE6+KetxId0oJCoFevQwBiitVcvs5ELcRrq8bE0peOIJY/C9adMoF+hLL+bwChMJZjcaxQOZLva/se+wcWFkhQowfDicPm1ieGEzWsP69dC1KwQGkjrha3ZcD2QSA5lOX/ZSO2PVlqyhuVqPR88nYNcumD1biolwSFJQ7MXNLeNKZX7+mTK1y9CNX3mNzynLhYzVZtCHqTzHvit+pH/8CQQEQM+eRktHjrM4v/h4mDLF6OJs1ozrvy1nrW7KFwxjAZ25SBmKEIfG0u2plPHzv1lI6tQxN78QdyFdXvbm6gpPPml8SSxYQJHRozPG/4rDh4uUIQlPTlGRosQQlvY3DX5ZhNcvv0Dt2sY9wPv0gdKlTf4gIle2bYNvvjHuhBgfD8A6mrKW5qRgXCdSlvM0ZhN12IOrC9Cjl3FcTYqIcBLSQjGLiwt06QJbt8LSpdCiBUWI53XG8hhLKcllYijKStoyltdZSEfi954wTk9+4AHo1s14XWqq2Z9EZCcmxjjIHhoKDRty45vvSYy/9fPyIpEU3KnOYZ7hRwYwmXo+x3AdOtg4rjJrlhQT4VSkhWI2peDRR40pMpJC48cTNnMmjW5s4QjV2MxDHKUqe6hDO1YYr0lJIeXXhbj/+iuULw/PPAO9ekG9enKGmNkSE2HxYuN2u0uXQnIy5/FjG4+xi7qEso1HWAlAPXZSiZOUJtr4I2HoJ/DSS1CsmMkfQoj7IwXFkYSEGNexfPQRavJkqk+cSPXzM4imJBfww4MbAKTiyhcMoyKnqH92B9U++hiXjz4yBqTs2dM4C6huXZM/TMGhUlON4vHzzzB/PsTHE4cPewlhF3U5S/mMdaMplfHYnVRKh1aC174wfm4yRIpwclJQHFGZMjBihHENy5w5lJo8mVLr12csjsKfRLzYT032U5PCXKcm+6l1eB8BY/6Ly5gxEBhIQFgYFC1qFCppuVhXbKxxQerixTz822/Gc4vNPMgK2mUcWPckiXrspAHb8eMi+PhA795GayQ01KxPIITVSUFxZIUKGQfg+/SBgweNGyj98AOVz5/kNT5nF3XZQX2iKcU2QtlGKIW5ziAm4H3oEJUOHYIZM4wL39q3N7rV2raF4sXN/mTO6ehRWLTI6NJauxZSUriOF4epThHiqGK590h5zuJCOoEcog57COSQcTO2sDDo/1+je9KJbg0rRE5JQXEWQUHw0UcwejQsW0aRqVNpsngxD6du5BKl2Utt9lIbF9Lx5nrGy1bShrLnz1N12mwKT5tmnAzw0ENGcWnd2vgLuVAh8z6XI7twwRhvLSIC/vwTDh5EA5cozSHCOEQgp6mARlGVoxkFpQKneYv/4Umy0drs+bJxdl5IiKkfRwhbk4LibNzcoGNHY7pwATV3LmXmzKHMujW01BEk3RwsELhKMTbQBACF5gHOUD39MNU2HqH8xhGoESPA0xMefNC4qr9pU2jcGHx9zfp05jp71mh53CwiBw7ctngbDVhPU65yq4XnShqVOElN9mfMU0WL4tm1Kzz1FISHy/haosCQ33Rn5ucHgwYZ07lzqHnz8Jozx7gCW2sKcYNH+J0jVOMklYjCnyj8WU04hbnOs0yjTNIl4wv05sjHLi7G2WINGxp/UYeEGAf481sXzblzxrUhW7ca/27bZswDkvDgFBU5QVtqcICKGKMVpOPCVYrjTQLVOUwgh6jKUeNkCS8v6NjTKCKPPgoeHmZ+OiFMIQUlvyhXDl591ZjOnoV580j+8Uce3rWNh29s4gbuHCeAw1TnCNWIowjFuZrx8rl04zqFqZh+igo7TlNuxwwK862xUCmoVs0oLvXqGd1vVasakyO3ZrQ29sWhQ8YxqJtTZGRG8QC4RCnOUp6z1OM0FThHuVtXqkNGQanFPspxznKMRIO/Pzz+PLsqVqTu0KEy9poo8KSg5Efly8PgwewKDqZlo0YQEUGhZcsIWraMoGNL0EAsvsaBYkADR6lKIl4co0rGZnyJpSznaaC3U+PwQWOI9F9+uf29Spe+VVyqVDEKm5+fMZUta/xrq9ZNXBycOWMUjZvTzefHjhmFxHJVOsAN3LlMSS5RmmDOZZSMuXTnAn4Z67mQjj9RVOYEQRzMmO+tEvEOC4bHXzZGhbZc93MlIkKKiRBIQcn/vL2hQwdj0hoOH0YtW0bRiAijayw6GgUMYgKnqcApKnKaClzAj1h8icWXqhzN2Nx+ahBBS0pwhRJcofilq5S4dIESm/fjS2zGDcVuU7iwUXh8fIzH3t63Tx4eRjatCTxzBqZPh/R0Y15yslE4YmNvn+Li4MaNbD92LEXYR22iKcVlSnKZksRyqzVVgdMU5xoAgRyiBFcoz1ke4Az+RFGIFGPFunWhRddb964pVSqLdxNCgBSUgkUpCAw0pqFDjS/sQ4dg3Tp81q+n5rp11Dz2OwDpKK5SnPOUpTxnMzZxkTJcwO+2v+hv8iCZd/go4/kqWpOOCz7X4yl88jqFuIE7V3HnIiW4gi9xgHGhZhKeaBQ+KGJQGV1OabhSgisZrYnDVOMqVUjGg0S8iKNIxhTIIdph5I+jCMtpf1s+V9IowRVKc4k0XDPmt+ZP44GbGwQHQ7OOtwpIyZJ52uVCFCRSUAoypYzjIUFBxs2awOgu2rYNl8hISlomjl3LeMlDbKY6hy3tE2O6SnGuUAI3bh9XbDsNuE7WXUHhrKYFawE4RCBz6JltzOF8ZJyCC2zkYY4TkOV6VyiR8bgU0TRii6VtYkzFuHarBeXuDnXqG6dN35yCg42z3oQQ90UKirhd+fLG1LHjrXmxscbw6Tt24HHgAOWPHqX8kSNwYiOkpWWsls7tV+M/wu/E40M8PiTixQ0KkYI7Nyh02wkBbqTig3Gsw2ib3JpcSCcVN7AUlEAOUYpoPEjGk6SM9okvsRSxtHgAPLhBB5YaXWqBgRDU0CicgYFQo4YxcrOciSWEVUlBEffm62tco9K06e3zU1Ph1CljZNyjR3E5dcq4GNAyhVy4ABf2QUrKXTcfyGHe5LMcRWnM5ltPChW6VQDLhxoDLN58/sADRvEoX16GnRHCTqSgiPvn5mac2VWlSvbraA3Xrhn3S09IMKbr129/nJxsXP+iFAcPHyaoZk2jCLi4GO/h65v15OEhxUIIByIFRdiWUsbYYTkcP+xcRARBLVvaNpMQwibkBltCCCGsQgqKEEIIq5CCIoQQwiqkoAghhLAKKShCCCGsQgqKEEIIq1BaZzGYXz6jlLoEnMzjZkoB0VaIY0+S2T6cLbOz5QXJbA9Z5a2ktS6d0w0UiIJiDUqprVrrhmbnyA3JbB/OltnZ8oJktgdr5JUuLyGEEFYhBUUIIYRVSEHJuW/MDnAfJLN9OFtmZ8sLktke8pxXjqEIIYSwCmmhCCGEsAopKFlQSvVQSu1VSqUrpRresewdpdQRpdRBpVS7TPPbW+YdUUoNt3/q2zKOUkqdUUpFWqbHMi3LMr/ZHGn/3Y1S6oRSardlv261zCuhlFqplDps+TdnQyvbLuNUpdRFpdSeTPOyzKgM4y37fZdSqoEDZXbY32OlVAWl1Gql1D7Ld8VQy3yH3c93yWy9/ay1lumOCagJBAERQMNM82sBOwEPIAA4CrhapqNAFaCQZZ1aJuYfBbyZxfws8zvA/nao/XePrCeAUnfM+wQYbnk8HPjY5IzNgQbAnntlBB4DlgEKeAj4y4EyO+zvMVAOaGB5XAQ4ZMnlsPv5Lpmttp+lhZIFrfV+rfXBLBZ1BmZprZO11seBI0CYZTqitT6mtb4BzLKs62iyy282Z9l/2ekM/GB5/APwhIlZ0FqvBa7cMTu7jJ2BH7VhM1BMKVXOPklvySZzdkz/PdZan9Nab7c8jgP2Aw/gwPv5Lpmzk+v9LAUldx4ATmd6HmWZl918M71qaVpPzdQF44g5wXFzZUUDvyultiml+lvm+Wmtz1kenwf8zIl2V9lldPR97/C/x0qpykB94C+cZD/fkRmstJ8LbEFRSq1SSu3JYnKKv4zvkX8SUBUIAc5BDm/YLnKiqda6AfAoMEgp1TzzQm30FTj0qZPOkNHC4X+PlVI+wDxgmNY6NvMyR93PWWS22n4usLcA1lq3uY+XnQEqZHrub5nHXebbRE7zK6W+BRZbnt4tv5kcNdc/aK3PWP69qJT6DaML4IJSqpzW+pylG+OiqSGzll1Gh933WusLNx874u+xUsod44t5ptb6V8tsh97PWWW25n4usC2U+7QQeFIp5aGUCgCqA38DW4DqSqkApVQh4EnLuqa4o2+2C3DzzJns8pvNofZfdpRS3kqpIjcfA49g7NuFQD/Lav2ABeYkvKvsMi4EnrGchfQQEJOpy8ZUjvx7rJRSwHfAfq312EyLHHY/Z5fZqvvZ3mcaOMNk2alRQDJwAViRadm7GGc7HAQezTT/MYyzJo4C75qcfzqwG9hl+aUod6/8Zk+OtP/ukrEKxlkvO4G9N3MCJYE/gMPAKqCEyTl/xui6SLH8Hr+QXUaMs44mWPb7bjKd1egAmR329xhoitGdtQuItEyPOfJ+vktmq+1nuVJeCCGEVUiXlxBCCKuQgiKEEMIqpKAIIYSwCikoQgghrEIKihBCCKuQgiKEEMIqpKAIIYSwCikoQtiZ5WrpSKVUvzvmd1bGPXj2KaWqZfG6r5RS39kvqRC5IwVFCPvrCZQAfrpj/maMYWcCgZeyeN2nQO+sio0QjkAKihD2NwSYrrVOyTxTa31Baz0HWA/UvfNFWusTlmUD7RFSiNySgiKElSilKimltFKq8R3zv1VK/WF5XA14GJh7l00dAmpns2weRitF/u8KhyO/lEJYTz2Mwfd2ZzH/5rzWQALGAJP/oJQqDXQHKiilfLNYZSPGTZuCrRFYCGuSgiKE9dQDTmit42/OsLQkanOroIRiDB+ens02PuPW/8usWil7gTQc49bNQtxGCooQ1pO5JXJTdaBwpvllgeisXqyUCgd6WybIoqBorVOBa5btCOFQpKAIYT31MO4pkVldjG6wvZbnnhj32bmN5cZik4BvtNZLMO6MVyeb90m2bEcIhyIFRQgrUEoVxrgv9547FjUHjmmtEyzPrwDFstjEcMDX8i+W7WR3YL6YZTtCOBQpKEJYRymMu/Jl3NZVKVUK6MXt3WAHgYDML7Sc+fUOMFRrHWOZvZssWiiWg/aFMc4EE8KhSEERwjrOAfHA00qpkkqpB4HfAB9uLygbgIqWwnDTROAPrfUvmebtAcoqpUrc8T4NMbrQNlr7AwiRV1JQhLACy0WKLwIdgFPAGIz7cbtxezdYBEZ3VXsApdRTGNelDLpjkzeL0J2tlPbAGq31ZSvGF8Iq5J7yQtiZUmocUE1r3SGXr3MFTgLDtdYzbBJOiDyQFooQ9vc/IFwpFZjL1/UAEoFZ1o8kRN5JQRHCzrTWUcDzQLlcvlQBL1iuRRHC4UiXlxBCCKuQFooQQgirkIIihBDCKqSgCCGEsAopKEIIIaxCCooQQgirkIIihBDCKqSgCCGEsIr/B1nrtphHGaQIAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# This code plots the uv track for a baseline given the antenna coordinates, wavelength and the number of hours\n", "dec =2.082175*(np.pi/180.0) #in radians, declination for HD 35298\n", "\n", "gmrt_lat =19.1*(np.pi/180.0) # GMRT latitude\n", "gmrt_long =74.05*(np.pi/180.0)# GMRT longitude\n", "start_alt =20.0*(np.pi/180.0) # Starting elevation angle for GMRT\n", "print ('This will show the uv tracks for a given baseline for two observing session of different lengths')\n", "#num =int(raw_input('Enter the number of hours: (integer expected) '))\n", "num_arr =np.array([float(i) for i in input('Enter the numbers of hours (e.g., 5 10): ').split(' ')])\n", "\n", "lamb =float(input('Enter the wavelength in metre: '))\n", "arr_typ =int(input('Choose the array type: East-West (press 1), North-South (press 2) or inclined (press 3): '))\n", "R =6371.0*10**3 #Earth's radius in metres\n", "\n", "def find_u_v(arr_typ,num,dec,gmrt_long,lamb): # function to find u, v given an array of hour angles\n", " ant1_lat,ant1_long =18.97239*(np.pi/180.0),74.05211*(np.pi/180.0)\n", " ant2_lat,ant2_long =19.08447*(np.pi/180.0),74.04933*(np.pi/180.0)\n", "\n", " if arr_typ==1:\n", " ant2_lat=ant1_lat\n", " if arr_typ==2:\n", " ant2_long=ant1_long\n", " if arr_typ==3:\n", " ant2_lat=ant1_lat+abs(ant1_long-ant2_long)\n", " \n", " h_i =np.arccos((np.sin(start_alt)-np.sin(dec)*np.sin(gmrt_lat))/(np.cos(dec)*np.cos(gmrt_lat))) #starting hour angle\n", " h_f =h_i+num*15.0*(np.pi/180.0)\n", " h =np.linspace(h_i,h_f,50)\n", " X1 =R*np.cos(ant1_lat)*np.cos(ant1_long-gmrt_long)\n", " X2 =R*np.cos(ant2_lat)*np.cos(ant2_long-gmrt_long)\n", " Y1 =R*np.cos(ant1_lat)*np.sin(ant1_long-gmrt_long)\n", " Y2 =R*np.cos(ant2_lat)*np.sin(ant2_long-gmrt_long)\n", " Z1 =R*np.sin(ant1_lat)\n", " Z2 =R*np.sin(ant2_lat)\n", " Lx =X2-X1\n", " Ly =Y2-Y1\n", " Lz =Z2-Z1\n", " u =(1.0/lamb)*(Lx*np.sin(h)+Ly*np.cos(h))\n", " v =(1.0/lamb)*(-Lx*np.sin(dec)*np.cos(h)+Ly*np.sin(dec)*np.sin(h)+Lz*np.cos(dec))\n", " #print X1,Y1,Z1\n", " return u,v\n", "\n", "\n", "lat1,long1 =18.97239*(np.pi/180.0),74.05211*(np.pi/180.0)\n", "lat2,long2 =19.08447*(np.pi/180.0),74.04933*(np.pi/180.0)\n", "\n", "u1,v1 =find_u_v(arr_typ,max(num_arr),dec,gmrt_long,lamb)\n", "u2,v2 =find_u_v(arr_typ,min(num_arr),dec,gmrt_long,lamb)\n", "\n", "plt.plot(u2,v2,'r-',linewidth=5,label='shorter session')\n", "plt.plot(u1,v1,'b--',linewidth=2,alpha=0.5,label='longer session')\n", "\n", "\n", "plt.xlabel(r'$u(\\lambda)$',fontsize=15)\n", "plt.ylabel(r'$v(\\lambda)$',fontsize=15)\n", "plt.grid()\n", "plt.legend(loc='best')\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulated source and tracks for GMRT antennas\n", "In this exercise below, we have simulated a point source and a 2' Gaussian using GMRT antenna array.\n", "Here left side is the image plane and right side is the UV plane. The constant amplitude is due to the point source at the center and the structure in the amplitude at small UV baselines is due to the extended source - the fundamental property of Fourier transforms!\n", "\n", "![SegmentLocal](data/Images/source.png \"segment\")\n", "\n", "\n", "In the images below, we should UV tracks for 1, 2, 5 and 10 hours, respectively and the resulting images. One note that while white source does not get affected, the extended source quality depends sensitivity on UV coverage. Why?\n", "\n", "\n", "![SegmentLocal](data/Images/uvtracks.png \"segment\")\n", "![SegmentLocal](data/Images/images.png \"segment\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise for the UV tracks for the above simulated image\n", "In the code below, the UV tracks are for the above simulated image using GMRT array. The observations are for 10 hours in total. You need to provide the fraction of the UV tracks. The number 1 implies the full UV tracks.\n", "Compare UV tracks for fractions 0.1, 0.2, 0.5 and 1 and see how UV tracks improve with increased observations.\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Give in percentage what fraction of the data to be plotted(e.g. for 10%, say 0.1):0.1\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAEKCAYAAAClutpcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztvXmYXFW1v/+uru5OZyLdGUhCyIQEBAQCCSGADIYZfxJUwqQSES8qCHgv8QKiX72gGLzEaxy4iExBZB4jo0wRBEIIQ5ghCQkZCJnnoadavz/WObeqq6uqq6qruqq61/s89dQ5++yzz67qrvM5a++11xJVxXEcx3HKiYpid8BxHMdxssXFy3Ecxyk7XLwcx3GcssPFy3Ecxyk7XLwcx3GcssPFy3Ecxyk7XLwcx3GcssPFy3Ecxyk7XLwcx3GcsqOy2B0oR/r3768jRowodjccx3HKitdee22Nqg7IR1suXjkwYsQI5s6dW+xuOI7jlBUi8km+2vJhQ8dxHKfscPFyHMdxyg4XL8dxHKfscPFyHMdxyg4XL8dxHKfscPFyHMdxyg4XL8dxHKfscPFyHCevqMKzz4IIrF5d7N44nRUXL8dx8sq8eXD00ba9884mYuFr06bi9s3pPBRVvETkBBH5UEQWiMhlSY53E5G7g+OviMiIuGOXB+UfisjxbbUpIiODNhYEbVYnXOvrIqIiMrYwn9ZxOjc7dsB++8Hpp8OPfpS8Tp8+JmKvvtqxfXM6H0UTLxGJAH8CTgT2Bs4Ukb0Tqp0LrFfV3YH/Aa4Jzt0bOAPYBzgBuE5EIm20eQ3wP0Fb64O2w770Bi4GXinEZ3WcrsD48fD22/DRR/C736WvO26cidiSJR3Tt2zYtg26dbP+ffZZsXvjpKKYltc4YIGqfqyqDcBdwMSEOhOBGcH2fcDRIiJB+V2qWq+qi4AFQXtJ2wzOmRC0QdDmKXHXuQoTtx35/pCO01V49NHszxk+3ERiy5b2X3/Dhtjw5OLF2Z/f0ACjR0PPnrYNMHhw+/vlFIZiitcQYGnc/rKgLGkdVW0CNgL90pybqrwfsCFoo8W1RORAYKiq5vDTcxwnZPjw3M/t3RuWLcv9/GgU+vWL7Y8cmfm527ZB9+5mbc2b1/LYihW59ylbtmyxz/DSS+b04qSnSztsiEgF8Fvgkgzqnicic0Vk7mp3oXIcGhvhq181S2ennXKzduIZOjS38zZvhkjEBCxk0aL05zQ3w3/+p/W9Z0+br4snEoG1a2HQoNz6lAnRKNx0k4lmdbUJ+Lp1cPjhrUXUaU0xU6IsB+L/XXcNypLVWSYilUAfYG0b5yYrXwvUikhlYH2F5b2BLwCzbGSRQcBMETlZVVvkPFHVG4AbAMaOHevPRU6XZ/JkeOgh2968OXfxiWfNGujfP7O6qnb9r32tZfnLL0OqdHvbt5uFmOr5s7ISVq2CurqMu5wVO3bAvvvaXFpTU2vRBLjxRth//8JcPxcaG+GMM+yBIBKBgw+2Oc3KYifUUtWivDDh/BgYCVQD84B9EupcAFwfbJ8B3BNs7xPU7xac/zEQSdcmcC9wRrB9PXB+kj7NAsa21fcxY8ao43R1GhpUTzlF1WQkP6+Kisyv/8Ybrc9/66305+yxR/LrDhmiunVr+76PZOzYoTp6tGq/fqp9+6pGIsmvX1WlWlur+uKLqtFo/q7f0KB65pn2ngnRqOrs2apjxqgecYTqAQeo1tW17u9FF+XWH2Cu5klDiqadqtokIj8EnsSE52ZVfVdErgw+4EzgJuCvIrIAWIcJGEG9e4D3gCbgAlVtBkjWZnDJS4G7ROSXwBtB247j5EhVFTz4oDlKDBwYc3JI5HvfM2vmwQdbHzvrLLj9dpg1y9aGrVyZ+fX33x8eeCBmeX36adsOFm++GbO8evSwOa2ddsr8mulQhdmz4bTTbEiwvh42bjQLKxm9e9s6uLfftjm3fNDQAMcdBwsXWpuffWZWcTQKd93Vun59PXzxizZ0qmrfy/vvJ2976FD7Ox98MEyblp/+tgdRnxnMmrFjx6pnUnac1mzZYje5DRtiZfvsA++8Y9vPPAPHHANPPglvvAFTpthQVDmyY4cJ6Nattt/YaCKdjEjE1rh16wbf/z5ccUV+PveOHeYhGc73LV9uDiiJzJhh83s339zSGeSNN5IvHO/dGw44wISvqgouvNAeNCra6SUhIq+pal7W0rp45YCLl+O0zY4dtvZr9myoqSl2b/JLQ4N5BqZy8d9lF7NqqqvzJ1aNjbYA/N13zWrcts3WySWbNwuprraHiaoq+OCD5HW6d7d1d6omxJ/7nFnDVVXt628y8ilexZ5ycxynk1JTY8N0nZGJE2PC1b9/bFFzXR3MmdN+sY53klA1oVq9GtavT31OZaVZRpGIOaaAiezChbE6AwfC5z8f26+uhkcesfdyw8XLcRwnSx5+GE4+Ga66CsaONeHKheZmuPpqePppG6IDa2vp0vRBjSMROzee+Lm1vn3NqzGkogL23rtEvATzRCf5GI7jOB1HdTU88UTm9ZuabN7o/ffNktqyBXr1MqtoeeICoQwIhatfP1sXtmED1NaaZVbO1lQ2uHg5juMUmH//d7j++tTHd94ZhgTxhcLhxz33hOuus7IxY2KOMLW19r7zzoWbmyoHXLwcx3EKzPjx8Mc/miANHBizvLZtS73oNxo1N3YwR432evp1Nly8HMdpF42NttbqkUdg/nzYffdi96j0OPNME6dJkzIXoYoKO89Jjmu54zjtYvJkEy6AUaPgFU8s1IqKCree8o1/lY7TSdm40Sbt4zMZx782bszPdWbMaLk/fnx+2nWcdLh4OU4nZfBgG9JLRW1tSzELo2BkS1UVLFiQ27mOkysuXo7TiVCF554zMQoXqmbKvvvmbo197nOxSOj5ioi+dq0Ns914Y8t0J44DLl6O06mYNw8mTMj9/NAay8WSeu01mDrVAsCKmPNGLmzZYuf3729i/G//Bvfem1tbTufFYxvmgMc2dEqVrVth2DBLapgvVq60NUWZEh9tIpvby6pV5kaejOZmd3boDOQztqH/OzhOmaNqMQRVLQRQPoULTFDq6zOr+/e/Z9f2xo2xObdUwrVqVeGFa/Nmu4aIObmcemr6+UKn+Lh4OU6Z8/zzlr6iosKijBeCTIYit22zeH/ZkCr/VkWFZVVWhQEDsmuzLaJRuOkmC6ZbVWXX2mmnmJXY2Aj3329LAEqN+nqLpXjIIZk/UHRWXLwcp8w5+uj8thcujA3nmQYMgGefbfu8vfaKbXfrljpdSDwrVrTc/8lPLA5gc7PF7csH27fDrrtajqqePa1v3/2uRVxvamo9tFlVBV//euslAB1NKFQHHmgPJ3vuad/Ja69Zmpn2zG12BjzChuOUOatW2ZxUYpTxVIi0vmHX1VlQ17VrLSL5HXdYeTZzVu+/bzfYX/4SvvWtzIb6+vTJ7hrpiEZNcH78Y9tubDQh2rAh9TUqK+17q601q7VXr/z15bbb4Kc/tWtHIpZPq0eP1nWbm+FXv4KHHrK/Ta9eNoz50UexRJfxVFVZrMNMHig6My5ejlPm9O2bOtU82HDeqFHw6afwpS9ZFuMwmGtzM1x7bX4yGvfoYak8isGOHTBypKW9T0VFheXZam42a/LDD5OLSbaomiV0+ulm1e3YYX+PxL58/vPmOfn007HYhlu2mLh+/HHytquq4AtfiOX06tsXZs2y63R13NswB9zb0HFKi9GjbZlASF1dzPKKREygM7UG0xEmiVy40KzU7t1teK+tucbKShu2TJdMsrYWdtstZnl169b5hMozKTuO48Qxe7Ytjo5G4a23TFTaQ2OjWVLvvmvWWTh899lnsaSR2dDUFBOuoUPN8gstrz59Ol+iyI7AvyrHccqemhobBsyGcK7pwQdtSK5nTytXNcsqF5EK6d3bHC1CVM1xZOJEuOyy9g/ROi5ejlPWNDTYDfHhhzt/5tx8Ul9vIa3aymK8886wyy5meSVaXcOG2TybSCyLMXSdTMbFpqjiJSInANOBCHCjqk5NON4NuA0YA6wFTlfVxcGxy4FzgWbgIlV9Ml2bIjISuAvoB7wGfEtVG0TkP4DvAk3AauA7qvpJIT+34+SLE06wWIap5kU2bTIrwGnJhAkx4erTxxZIx1tejY1w6aXwjW/E5smam+Hqq20e6uij7bhbUMWjaA4bIhIBPgKOBZYBrwJnqup7cXXOB/ZT1e+LyBnAV1X1dBHZG7gTGAfsAjwN7BGclrRNEbkHeEBV7xKR64F5qvq/IvIl4BVV3SYiPwCOUtXT0/XdHTacUiE+FFMqliyxeZZ4NmyweZemJnPRnjixMP0rVerr4Ygj4KCDfK6pI+ks4aHGAQtU9WNVbcCsosSf0EQgXCp4H3C0iEhQfpeq1qvqImBB0F7SNoNzJgRtELR5CoCqPqeq24Ly2cCuBfisjpM3duywBcGZCBfY8Na2bS3LBg2Kudefckp++5cpGzbYZwijaXQk3bpZ0sw//tGFq1wppngNAeJXhSwLypLWUdUmYCM27Jfq3FTl/YANQRuprgU2DPl4Dp/FcTqM0aNtwWs29OxpQhHvNZcr779vbUUi5i6eKfX1tog5jGVYV2flqqnjGjpOKjw8VICIfBMYC/x3iuPnichcEZm7evXqju2c48SRrVddPL16WTDc2lr4+c+tLHxPx7ZtMdHZe28ri0bTi87mzS2TXdbUWNSIZKxcmd3naA/hWq2//c3zhJUzxRSv5UD8SPyuQVnSOiJSCfTBHDdSnZuqfC1QG7TR6loicgxwBXCyqiYNd6mqN6jqWFUdOyDfkUIdJw2NjRbwNhSBd99tX3thMNxf/MKsnl/8InXdaBSuuSbmzJBIMtHZvt289Hbaqe2+jB1rw6D9+7ddNxeamiyO4U47WXSKXr1s3dbdd9uiZc8TVr4Uc7T3VWBU4AW4HDgDOCuhzkxgMvAycCrwrKqqiMwE7hCR32IOG6OAOYAkazM457mgjbuCNh8GEJEDgD8DJ6jqqkJ+YMfJhcmTW6Ya2Wef9rWXGAw3Hffea+uSEunXDz75JLmojR4NqQYn6uryG0Mwnm3bYMQIm8NqaDDh3bIldWqTvfeGSZPy349cCdedPfywie3OO8Ptt8dCeTktKZp4qWqTiPwQeBJza79ZVd8VkSuBuao6E7gJ+KuILADWYWJEUO8e4D3Mxf0CVW0GSNZmcMlLgbtE5JfAG0HbYMOEvYB7za+DJaqaZWIHx8kPmzfDkCH2/oUvWATxGTPsJpxtrqyQ11+3qOS5MGkSLF4cE7Dly23dUzrefBOGDzcBO+88+NOf8ucU0dQE558Pjz1m+6omVFVVNo+Xznm6ri5Wd+xYePTR4iS4bGy0yP2hk4pq6hiHkUgsSLLTEo9tmAPuKu8Uitpam5PKli1bzAoKrY/Vq+Hxx20dWLnS0ADHHQcLFsTKNm5sO9XK4MExy0s1JlT5WjQcjcJf/2pWEthQbphqJYxFGI3anNrNN7cW1PfeS22Zgs0j7rVX57S8PLah43RSli6NWV7pSPXM2aOHpUjpDEycCP/8Z/JjFRUmUvGWV2Vl6rQjudDUBD/8IcyZE4vq3qOHzfOlGnodMMAWMH/0kYlUOvbc05YshJZXba0vfs4GFy/HKSF697aoGCtX2o0tGYsWdWyfisXDD7e2vCorYdw4s2ryYZE0NcGFF5rQhA8MoVCtWmVDedmwebMt+g4ZONBSocRTWQnnnGNDh8UYtuwsuHg5TgkycGBr62rVKisfOdIEbMSIonStw6iutmG49hCN2tDb9OmxKO5hyhGAZcvSD+FlwpFHxrbDOIcbN3qk+ELjX6vjlAnxa6pGjoxtf/SRJZvs6jQ0wIknmjW1ZYtZsatW2aLqXKmpsRBS/fqZFVZbG3vfssViHY4dm3m0Eyd/uHg5Tonz4Yc29PT88xaPL5E99kjvZddVmDgRnn02+bFevSxFSaLlVVlpcR/XrDGBUjVxqquzYcnO5jDRmXDxcpwSJ5wzSSZcISJ2o/3GNzqmT9mweXPLBcsHHQQvvJD/DMEPP9za8hLpnF57jrvK54S7yjsdxfbtuXnPTZgQs0Iee8xu6tlQX2+iuXhxrKytNV719TBmTGYRQA49FF58Mbs+OeVPZ4kq7zhdDlWYObNlzL/4V6LzwOjRuV0nfvjspJPS1922Dbp3bx2HMF64wFz4w/p9+7bue01NZsJ10EGph/fyzfbt5thyyCEmrk7nwcXLcTqQefPS587aeWcb/goF4aOPbO1QewijUSSjsdHmenbsaLudMHnjXnvFsgZnQiRiwX+bmky858zJ75Bhc7PFZ9x119av3r0tjNXs2WaNOp0Hn/NynA5g82YbcmsrOgS0zq+1ejXcdBOce27qc7ZujcUZnDIFrr3WbtZPPJF+rmfy5NTCFYmYt17fvi3L33/fhCFRwCIRGzZ8/vn8z2c1NMDxx5tr+44dZimq2vaOHW3nAxs/vuOsvbYI15bFe0FWVbnnYrb4nFcO+JyXky25hn0KGTTIojqcfbaFJornvffMGsqFxkY49VRLzLhwYfJAu6p24z/mmFhZdbUtpK6tze26ie2/8ooJaarb0fLlrRNqJmNIQpa+7t3hrbfsPV+oWrzIp56CH/84eTSMxkY46ywTqtC1fv16E6alS+27TqR7d3jppdyHisuBfM55uXjlgIuXky1vvw377Zf7+WHsQrCb57x5sP/+6Z/St21rLUbposE3NJhH4yuvZNanbt0yG25MRX29OW6sXWt9yoT+/W1uLd7y6tnTskU/8kj+4heCfR8nnBCLsiFibvarVlnUj6YmW283enTMxT5cB7ZwoYlUOsLwUNB1LC+Pbeg4ZUZ7hAtaio1I+qfzpUvtZp6MtWth991bxuarr7eo89ks5q2ubl82ZrBhzddfj+1XVrZcfA0WIWPVKht2ra21+cBDDsnPDT4+BUnv3rEYgz172sPB1q1WJx2LFqUP13XkkbGoG6Hl1b+/DQ17eKj24eLl5ERzsz0pPvqoBUjdsMF+/HV1sSfPtWtjT+fh0IqI3Sj+7d/gZz/rvAFIt283V/MlS3I7v18/+/7A0qJkQyrhCtuNjxUIJiJtCVeq+a/28OyzZnlt2WLW1Kef2nZDQ2z4sKEhNk+4YYOJwauvZja0FlpOGzfGxOm99+z/sbnZhLG9Hohjxtj3nWh5DRhgQ4qd3ZIqJi5eXRxVuznecotNhr/5pv3Yq6rsx71uXfvaTzZPsXkz/Nd/WdvTp3fOH/fo0bkL19NPW3TxXFmypKWArViROsgvmIiElteQIebhmK/I7GBCcdVVcOONsbIwGnx1tc0PZRJfcMgQS4651172/YRWUbxlE8/bb7fv/7dXLzj8cHsQqauL9XvDBujTxyzcnj1tyLN3b7Oiwj56PMPC419xJ6S52VK3v/56LKfRmjX2xN3eoZ58ctNN5oH3n//Z+YZP3nwzveW1bp2JxPjxLcu3bzcrpD0MHZpduKhu3dpO35EJjY1w2mlmGYXrvrZvN8spU2eVMBdX2P+KChPeOXPgnXdg6lQL1pupKFVU2G8gnpoaa7+pyR7SjjqqpUNHZSX86EcwbVrM+k283gsvJO/Dc8+ZlTp9emb9c3LHHTZyoFgOG9Eo3HmnpVwInzrXr7ft+fPtx9jYaGP1TU0d3r2cGTMGXn65tUt3U5PdRFTtZuBPs6XNWWfZ/2c6Qm/AeMtLpG2vwDffhMMOi1nyvXtbrEJIbXlVVNgc2syZZu2rmkXas6dlLA7jGw4YYBHg49t66aW286L17m0WaxiKqqLCI8m3hXsbFplCiJcqzJ0L//3fdtNev95+RJ98Yq+KCvuxdeUoAbvv3nKIccgQm3O75BK4557YdxMORanGtquq7Du84AL4f/+v8861FZNUllckYp517fEGDH8fU6fa7+CRR6w83hsQTHBWrbJrNjW1/yHuwANNmBKprs6/d2NXwMWryBRCvObOhYMPbj3E4eSfqVMtW22xGTHCHkyGD28diqkrE+8FGEaAD93UN2+2+bENG8wKa8/vpbraQlUNGBCz2kLLy9OdFAZ3le9kNDaa26wLV+5UVMSGnNqyvKZMKXx/PvvM5m8Seeklc/WG2NqmTz4xIQsFbMWKlgFw33ijcy5cDf/vV6+ODb2pmqdqGIoqV2pqbOiub19bD7dpU0yYNm0yh49LL3ULvJxxyysH8ml5NTTYk3cpOVKUCz16xOY1amvNQaWhwYQqjKPX2NjynH33NS+09vK731mIn0RHk40bLT5hQ0Nu7f7Hf8Bvf9u6vDP+TNuaIxs40IYbU1leFRUm+rvtZvVdlEqfTjNsKCInANOBCHCjqk5NON4NuA0YA6wFTlfVxcGxy4FzgWbgIlV9Ml2bIjISuAvoB7wGfEtVG9JdIxX5Eq/mZvNIW7DAnhJraqxs+/Z2N+10AHfdBaef3rKssrLtha3Z0tUsr+3bLXjxZZe5CHU28ileqGpRXpi4LAR2A6qBecDeCXXOB64Pts8A7g629w7qdwNGBu1E0rUJ3AOcEWxfD/wg3TXSvcaMGaP5YOpUVfu5luarokI1Eil+P0rx9bvfqTY3t/6bFuJa48er1tfn5V/OcYoKMFc1PxpSzNU144AFqvqxqjZgVlFisoiJwIxg+z7gaBGRoPwuVa1X1UXAgqC9pG0G50wI2iBo85Q2rlEQGhvtaf1rXzNPuT32sCGvQrjWhi7IIb1726LLAw80F98ZM2x466qrzMV39Ggr/9rXrLy5OZYS3TH+/nebm7z44uRr0wqxgmL27PRpVBynK1JMh40hQHzoymXAwanqqGqTiGzEhv2GALMTzg3jSSdrsx+wQVWbktRPdY02kizkxuTJ5tadC1VVNrT4uc+Z4NXW2uLjxkab/3niidxcd3/6U3sl49BDs8vd1B6GDLHPUVOTn3mpQvCVr9h7GGUhkTFjzF7K5+PP+PHmeec4Tgz3NswQETkPOA9gWLrgcW0wY4ZZNI2NsbVcYUy0igo4+WQL0/TGG7a6f8IE+MlPOm7sP4xZeMstds3EOHjtpbLSHFQqKuALX4C77269OPmqq0pXvEIGD06fouOee2zNU3tQbd/5jtOZKaZ4LQeGxu3vGpQlq7NMRCqBPphTRbpzk5WvBWpFpDKwvuLrp7pGC1T1BuAGMIeNrD5pHFVVdsMuVa691hbxJjJ4sHl7gU2of/ppS9f+3XaLhfXZfXf4xz9yX8A5Y0bbdYpNfFT2REaONLf3ePf3bEmWssRxnBjFFK9XgVGBF+ByzFnirIQ6M4HJwMvAqcCzqqoiMhO4Q0R+C+wCjALmAJKszeCc54I27grafDjdNQr0mUueKVNMhELLSwSuuAK+9a3CxR9UhRdfhOOOKw9Py1RDhiGhYLVn4bEvnXCcNsiX50cuL+Ak4CPMQ/CKoOxK4ORguwa4F3PImAPsFnfuFcF5HwInpmszKN8taGNB0Ga3tq6R6pUvb0PHeOMNVZGO8xTs3bvl/v/8T0vPwfnzY8c2bMj+84wYYeeOGKH60EPZ92/58vx9t45TSpBHb0NfpJwDnkk5v6SzvERszu/SSy0VxvLltgh48eL8pnYvNIsXt060mAz/OTqdGQ8P5XQqROCLX0zvAAHmyFKujM3w57p6tcXacxwnPZ0si5LjlBYLFpg4r23lApScnXe2+kuXtl23o2lqsgzY/frZmsGddzY3/q6c6cApHm55OU6B+PhjGDUqt3PjV2OsXGlCkS1btlhg2sT4jt26wVNPmbW7fbstXVA179D4RJBg242NdmzbNtixo2X7q1fbco4XX8y+f47THly8HKdAfO5z+Wln4MCYoESjcPPNcN55uc+P1dfDEUdA//5mEebSTq9eNue4227w7LO59SMXmposf9u0aZ7wsavjf37HySPvvQf77JPfNleujG3fe68N3eWDNXExZPr1a9vyikRMrP75T7Pe8okqzJkD558fi1ASRpGPZ9kys/aam+GPf8xvH5zywsXLcfJIe4Xrk09aDhkmMmmS3dAzsbxqauC66+CMM+Cww0wEw5QxYY6zykr44AMLy5VvGhvt2osXmxBt2ZI8uaSIvV57LfO2k+VKKxZNTXDRRbG1eZEInHKKRcwv1NpIx8XLcfLKu++2T8CGD29dFonE0qxEIpbm/rvfza7d11/PvU8h8WLUu7eVqbZMZxIK1JYtZiWtWpXdNWpqbElEKsursjKWLqXQRKNwxx3w4IMxa7B/f7P8ILb/7LOtw6jdf7/1NTFljpM/XLwcJ4/svXdri2j0aJg3L/c24/ODNTfbHFiiE0ZHMHkyPPBA9uftuaf1OZ3l1bu3DU0+8kjuYcVSEYrQAw9YPNFoFN5/377LpiYTzLo6c2BZty4mSvPmwSuvZH6do46ydkLLa9Kk/H4OpyUuXo5TYGbPtptavKdee4ifA+tIZswwZ49MLK/1621B+YQJJhyJwZfzQZhe6J13zMLr1cu+44oKE6iKChPCTZtaW3DxbNpk53/4YfLjY8bYUG4qy2v9evjxj20tX+GSKTmJuHg5Tg6sW2dODgD33Wc50MIb12efFXZOZpddYmur6upgyZJY0ORcaWiAY481IaisNDf/xODAzc02LBqNthaDTZtMOLp3t2OffWbv991n/Z0+Pf314+eNFi+2Icd4B5KaGqu3fXtMlNasabnGbOPGzD9vJGLxKdNZXuvW+dxVKePi5Tg5MHBgbPvUUzv22vE37PXrLXr9mjV2Y99jD3sPhy5FYiIQjWaem61PH/MsDBExkczWeuzVCxYtMnFPx9tvty/9To8eMes20fJqaDAxGjXKHioqK+F3v3NX+3LH/3yOkwMrV8Ysr2IT3pzD+Zx80NwM8+e3LhexlDeJqMYsLzBhqK83C+7vf8/8ukcdZRZUOssrGrXhycZG+xsMGwY77RTr3847w+23F2ao0ikdXLwcJwf69rWb6gcfmHdcPJdcYqllCjV02K0bfOMbtlgZWg/h9euXm+XVp499rqoqs2S2bo0dE7Fjs2ZlvsarsRHOOqulw0kyVK1fW7ZYUlZV89isrbVykdg2mIUWbq9dmzz0ViRic22lhKotB7h9RTfgAAAgAElEQVT1VjjkEB+ObC8ZiZeIDAdGqerTItIdqFTVNFOgjtM1+Pznk6+3KsRN6Xvfs4W5lZUmCLvuauu4qqpMnLp3N6eDQqzZyoR4sVKNZQjfsMGG9BLLQkHautVc+aPR7NZ6HXhgzHEkJLS8OjqhaWOjPVD07Wv7qjaUGwZZVoWFC+GFF+xvdd117krfXtpMiSIi/wacB/RV1c+JyCjgelU9uiM6WIp4ShSnLVavzi0eYVtcdFHbzg+FprHRrIa1a1uK0sKF7QsoPHgwHHxwa2sr0fKqrIRzzimM5dLcDL/+tYmpqr3X1JhIHnCAXS9elEKeeSbzObuDD4YLL+yalldHp0S5ABgHvAKgqvNFpAA/S8cpP3bsMMeGFSsKe52KCgsLNW1aYa+TCZMn2yLcVBx1VEvrqi3Lq6LC1scVyokiGoW//c2GWRcujC047tYt5oAS7q9ZY31LRltW4W67mccmJLe81q41a/HSS21Y02kfmfyr1KtqgwR+wCJSCXjKPMfBUoK0R7gqKuBLX4LHHsv/4txCMWOGDX0lWl4DBhR2vVM4Z/TnP8M//hGby6urM7f2UIzihUnELMUwdFO21NSYQ0hblpd7MHY8mXzV/xSRnwDdReRY4HwgC/8hx+m8zJ6dueVVWWmLYevqCt+vQlJVZeu3ciUahTvvhH/9y7bXrDFvyVAUQqtl0SJ779bNLJWtW5MnJF2yJPNr9+ljQpTK8mpstIzXo0aZOLoglS6Z/FkuA84F3ga+BzwG3FjITjlOuVBTA59+mt059fUwbpw5V0QidsMUsZtkfNinbEJAjRxpC4hDV/VSIbSWbrrJhHvAABu6e/rp9rXbvTsMGZKZ5dW9u4nWD38I3/xm15tn6qxkIl6nALep6l8K3RnH6QpMmABvvdW6vKEh9zYXLbIYiqlCHBWLefPs8yYLzzR2rM0BZWp5NTebwE+YYE4rbhF1bTL5838F+B8ReR64G3hCVZsK2y3HKQ6qFkX8619PXefSS+FXv8pu0l0VXnoJTjihZYSMfPLRRzak99BDljLl5JPhD38o7k1+//0t6nq85dXUZC7j0aiJbXW1DSOWy5yfUxq06SoPICJVwInA6cAXgadUNcukDJ0Hd5XvvLz5pk3Mt8XUqSZiyVC1+Zwjj8w923G+qK42q+bqq+Fb38p+yKypCS64wBYn9+gB27bZe0ODrVFatsxEM574YL3JypOlSjnhBHj88aw/nlNm5NNVHlXN6AVUYVbYA8CaTM9L0VZf4ClgfvBel6Le5KDOfGByXPkYbA5uAfB7YiKctF1AgnoLgLeAA4Py0cDLwLtB+emZ9H/MmDHqdE6iUdX771e1W23y16WXqjY1pW7jpZfSn1+s1/DhqqNHq1ZVqUYiqiKqlZWqNTX2HpZFIvaqqLDyQvVnzz1VjzxS9dhjVevrO+ov7BQTYK62QzviX5ksUg4trqOAWcA9wD+0HUOHIvIbYJ2qThWRywKRuTShTl9gLjAWc81/DRijqutFZA5wEbb27DHg96r6eKp2ReQk4ELgJOBgYLqqHiwiewCqtnZtl+Aae6lqipUehlteTiL19TaH8847xe5Jx7PnnjBoUGxf27C8+vY1SzDfLvXNzWYRL1nSOm1JuB9SWelpTIpBRy9SPhub6/qequZrtH4iJoYAMzBRTByEOR4bnlwHICJPASeIyCxgJ1WdHZTfhjmVPJ6m3YmY04kCs0WkVkQGq+r/DXio6qcisgoYAKQVL6fr0dwMV1wB11xT7J5kztChduOORzWW1iQatXm7ykobHlSNRWMP64bR2ocPbz3HV8hIF2GKlJUrY84coQjFC1L89vz52T08PPaYDe+OHp3fvueCBl6Zt9xiGQtWrLC/ny9oTk2b4qWqZxbgugNVNVwZ8xkwMEmdIUB8sJllQdmQYDuxPF27qdr6v9U5IjIOqAYWZvthnM7PtdeWhnD16hVzbIhGW0aD6NXL0qP06GHhh846qzhu4WGIpTfeaCk82bzPn29pUnIhzGgMbVte++/fro/aJqE1uHSp9WXwYBOmwYNtiUVY9vbbMHNma4/TSCT13GpXp03xEpHxwB+AvbCbewTYqqo7tXHe08CgJIeuiN9RVRWRvE9rZ9OuiAwG/orNqyVNKiEi52ExHhk2bFje+umUB1OmWCij9ghYdbXdjLZvT12nqqr1mq+qKth3X/jnPzOP6F5IQish9CBMdHN/+eXchSeRL33JFgxnYnlFIoVNHhk6rzzxhA19bttmCTu3brX3+nrzUl29OiZO8+bZQvZMGTcOvvzlmOU1ZUr+P0dnIZNhwz8CZwD3YvNPZwN7tHWSqh6T6piIrAyG7VYEwrEqSbXlxIYAAXbFhgGXB9vx5cuD7VTtLgeGJjtHRHYCHgWuCIciU3yeG4AbwOa8UtVzOgfNzfCTn8BvfpO+Xo8edqPZKe2jXOci3dqteBKFJ5v3QYNsSHLMmMLMSTU12aLluXPtc8SLUPgOLcvWr49ZcKmierz3XuuycePMgzWd5eXDhNmT0QoQVV0gIhFVbQZuEZE3gMvbcd2ZmCfh1OD94SR1ngSuFpEwmM5xwOWquk5ENgUW4SuYmP6hjXZnAj8Ukbswh42NgcBVAw9i82HtCHjjdDauvbZt4QJ7+u7TJ7YfHwIqGrU4gD/+sR2LT64Yn2cr8Ri0PJ7Mp6qqCg4/3HJWdXTSxcS1W4mW1/r1hXWGCMNLvfiiidzSpbZurEcPu140au+pRKhnT/jkE4vKkQs1NZYKJxPLy2MeFo5MvtJtwU3+zcCbbwXQXqN8KnCPiJwLfAKcBiAiY4Hvq+p3A5G6Cng1OOfK0HkDi694K9Adc9R4PF27mEfiSZir/DbgnKD8NOAIoJ+IfDso+7aqvtnOz+eUOVOm2M0tEwGLp6nJJtwbGuDee+E73ylM/8DiC1ZVdXzSRRETprHt9BkL54OWLYvd7MP2w/1XXzVxGDjQLKStWy3sUzbxDNMxaJBZzZlYXn36WESQMKeaU1wycZUfDqzE5rv+HegDXKeqGWav6Xy4q7wD9oR/ww3wgx/EyrqC5ZUt0SjcdZfFNFy+PCZOjz+e3XxQOurqzFklU8srnB+77DIfputI8ukqn1K8ROQU4CVVTTYf1aVx8SpfGhvhxBMteWDIV78Kd99d+iJQrtx9tzlRJLvVVFTY97/ffplbXr16WSDi/fbzuaJyo6PWeX0T+JOIbANeAl7ExKwLLsN0OguTJ7cULrBYhqni6vXubdZC4oLbYtHQAMcfb0NtO3bY/Ev8e69edlO//fbSEeNJk2yIMNHyWrmysE4ZTucmk2HDEcChwesQYBjwqqqeVOjOlSpueZUvySyvtujTJ3V23Uxoboaf/9yC+XYUvXvDU0+Zp5sLg1Mq5NPyatPxQlUXA68DbwBvYu7nJZY1yHEyo6rKckmFEfbWr09fv3dv82ZrD9de27HCBTbENn487L47/PWvNg/kOJ2JdHNeP8EsrQHAh8Ds4PVW4DLfZXHLy0lGQwMcdxwsCFyZotHMMix3BCNHwkEHFWc4UdXWU11zTWoRVTXrtq7O4w52ZjpqzutsYCvwd2zO6xVV3ZiPizpOOdLUZJ6FDz4Ycz5QtaHI6urc1w11BIsW2esf/zBr7OCDW64/ikbhb3+Dm2+OeUNGo7boNsw/Fkb+COMgxlNRYcOr++/f2o1840aLDtKcxSNvKcUddEqTtHNeQWT3cL5rPNALmIc5btzSIT0sQdzy6ppcfDH8/vfF7kX+qKw0wQndytOFrcoHI0akzpXWkZZXuMj5X/+y7TVrzJvxsMMKF1rKMTrEVT7hgpVYDq0jgO8BI1W1yzqmunh1TVJZXm3Nm5UzFRWxtVK5Wl6hd2G+I02EMRZvvtlESDUW6SMdH39sziyJiJionX56/vrotKRDhg1F5GTM4joM2AdL2PgicAk2jOg4ZUtjo60vevTR9rcViWQ3JNbRVFZa6CSw4c1o1AQIYoGAm5pg550t7FEkUjjBCVGF11+3hcrLlrVdP562IrFnyoQJNoQab3lNmpRbW+0hPtzV4MFWpmrzpbvsEttftcpDTcWT7mv4EfAs8J/Aa6qa47+I45QekyfnR7igOMIViUBtLQwZYuu7evSwoT8R266vt8C406fn52bX1GTDptFoy1BO0HJRcTypygHWroVHHrG+t5eDDrKwTZlaXhUVhRsiVIU5cywlTa9eFphYxP5OifVCcXrnHbjnnszaj0Tsb+qkF6+HsWjy38XiBd7hMf+czsKMGbBpU34ErLIyeVT5eGeOeMoptFPIJZfAddflv92jjzbrJxs6IhJ7aBk+9pjlJfvkEysPE3gmi4W4ahV07271Q557LrPrdesGX/mKLTAPr5/M8po2LX+fsdxJKV6qOh2YHsQ2PAOLJt8duBO4Mz4LseOUG1VV9uTvZMa0aWZ95cvyEim88DzxRGwubPlyS3AZjcKWLSY2W7akFqNotKUIZUvPnrZAPBPLq6ICTjjBrEdfGpA5GTls/F9lkQOAm4H93GHDHTa6KmFSwpdftiGimpr2t6kKr7xiw5mJP8lIpPRCPhWCMIDvxx/HvoNo1PKHLVliw3BgopNqO7SAPvggf/0SgT32sLbbsrwGDjRL+5FHUocc68p01Dqv8GKVwImY9XU0lhDyF/m4uOOUI5dcYtHkwYaJRo1qXSece9q6NXkbicfr69On+fjgA7tpdnT6k1xIZvmEVkay4bBw/8034YEHCtevESNg+PDMLK9evax8wgS3ikqVdN6GxwJnYnmw5gB3Aeepaoqfo+MUh2gU/vIX+P7309erqYHPPmuZPDIXpk0zR4NbbzVPt/nz29dePCKt54BCy2vGjPxdp5DMmwfHHpu/JQQisOeeZtVA25bX1q0mUvvvHxPGSMQjz3c20oWHeha4A7hfVTvxSpbs8WHD0uLuu+GMM3I7NxKxqObPP2+T5tmwY4dFgEgW8ihby0sE+vaFWbOy70epkavllWx7113d8ulMdPgiZaclLl6lRaaWVzoOPdTW2TiOUzg6dM7LcUqdigr43vfslY6NG21tVDyh5fXss9lft7HRLL5Fi8wq6NXLorknI5fjzc3w0Udm4Ym0XAwd/8wpArvtZudXVsLnPtf5nTscx8XL6TL06ZM8m2+ufPObhXUwiEc1FhUj2bGFC2P7r71WPs4djpMrHoLScbIgdGnP1VorJBUV5tRw6qml79wRpkn54Q8tmr3nG3OyxS0vp8vS2AinnQavvhorEzGvxDDCeuJ+fb2FICpFolFYudKswUGDzEOvsbF1+pNc2r3jDmt30SKLNjFwoK19GjQoFpIpDM/Uv3/q99Wr7TsVgYcfNmvyuuusbx4Q18mGojhsBKlW7gZGAIuB05J5NIrIZOCnwe4vVXVGUD4GuBXL6PwYcLGqaqp2RUSA6Zjb/zbg26r6etx1dgLeAx5S1R+21X932ChNmpvhv/4Lfvtbu2mHcz7duplI/eEPLW/gZ51lAVG7At26xWL+7dgRW1gdv719u1lvqi2Pq1oswkLwhS/AZZd5KpKuQtl7G4rIb4B1qjpVRC4D6lT10oQ6fYG5wFhAgdeAMYEYzQEuAl7BxOv3qvp4qnZF5CTgQky8Dgamq+rBcdeajmWMXufiVb5cc43dCNPRs6c99VdVmRW1bVvb7fboAb17m5VQqJt4OVFTYwt+22N5DRhgDxL5jJKeLHhwYogqEV/3VUw6g7fhROCoYHsGFrXj0oQ6xwNPqeo6ABF5CjhBRGYBO6nq7KD8NuAU4PE07U4EblNT6tkiUisig1V1RWDFDQSewITSKVOmTLF1U4mW1/btsTmVVOuu0rFtW2YiVyr07WuCm8rCytby6tHDXkccYa98W0nxKUEGDWotOvGBeJOJUrh9332WJiVTQgErJcIQWYsXw/HH+/q2dBRLvAaq6opg+zNMPBIZAiyN218WlA0JthPL07WbtC0RWQlMA74JHJPbR3FKhUgErrzSXvHU19u8z8KFsRT3VVUWHSN+u6kplvK+FKmoMKtht90smkTv3vZ5tm4tnHt8fN6t5cvtO7zggpZikuo9mfAkO/avfyVPDpkL++xjEfszsbymTMnPNVPR3GyjAdmkzIlPjzJtGjzzjC2Ed1pTMPESkaeBQUkOXRG/E8xV5X3sMsN2zwceU9Vl0sbjjYicB5wHMGzYsPx00ukQunWzuHltEXoSfvvbFrNw2zZ7X768pZOGSH5d7uOprLRXc7O9h2u7hg41J4l+/WDDBhg2zN7r6kzUGhrguONsHduGDfa+fr31Nb5s40bYe28brotEzK3+llus7WQiJGLH85F3q63PfcopsO++uVtehUjWGC9A0Si89ZYNHYdhqcL/hW3b4MtfbmmRvvBCbqJcUwPnngvnnGMhrpzkFEy8VDWlJSMiK+OG7QYDq5JUW05sCBBgV2wYcHmwHV++PNhO1e5yYGiScw4BDheR84FeQLWIbFHVVjMnqnoDcAPYnFeqz+aULyIwfnz7IpK31wmkqSm2nquxMVb+0Uf2ygfPPWfCdc45Fng21cLpeI4+2gIQZ2JxZWt5rVxpfRkzJn9DZMki1IOJ0Asv2Jq/MPaiamsxCodPQ+/KTIj3Wg059lizBDPF06NkTrGGDWcCk4GpwfvDSeo8CVwtInXB/nHA5aq6TkQ2ich4zGHjbOAPbbQ7E/ihiNyFOWxsDIYXvxFeTES+DYxNJlyOkykzZpjoNDXFrJ3497o6uykmO9aeusne01le06aZgD37bHrLq5DODdGoxaXcbTeLg/j441YuYoF1P/nEPr+qid3q1bBuXez8UHSGDrXPF8+KFfDhh/ntL5i4jB4di4iSzvJyp5DCUizxmoplZz4X+AQ4DUBExgLfV9XvBiJ1FRA+z1wZOm9gw323Yq7yjwevlO1iHoknAQswV/lzCvjZnC5MVRXce2+xe5E5Y8faqz1kMrcza5aJVW0tLF1qN/6VK/OTd+u117KrLwJHHhnbT2d5VVRYXQ8SXHp4YN4ccFd5p1xoboZf/9qcLlKh2nIuLNWcUaqhuFzndrLloINsqUNIOssrTKMSRq8P8WG54tIZXOUdp1MTOlC8+27LchHLsBvv1Rh6PDY2tj6WSLdumXtE1tTY0GGmywOee86G7q68Mrk7/L33WjzHZM+7bc3tJLO8KiqSC0xImBZlyBAbWjzjDF/I7MRwyysH3PJKT2MjnHiiufmC3ZBXrmwd0b0zc+KJNo9TrowbBxMntpyzSWV5+dyOkylueTklzeTJMeECs0Lq6lrW+drX7EbYmdJ2hDEAH3zQLKSddoJNm4rdq9yYM8ded9wBRx0VG04866zCX7u5GaZOhWXLsnOXdwuta+Hi5eSdGTNszU28gCXywANmkYXU1cGSJbGU7rnS1GSLaGfNap2pONzv1g1uvNGsi3zOe9x7L3zrW/lrrxR49117Pf64DQ1WVtpnfPJJ+44//TT2vfbs2fp9yxb7jpMdC9937IDDDrO5qxUrbCF0rnNoFRVmAXqQ386PDxvmgA8bZk5Tk1limeSW6tev/RHbL74Yfv/7zOvHLzgOt3v0sKf/MMRUGIkjjM4BLcvBRLdPH7v5psq75WROZaUNW+63X2lYXokelT4flxs+bOiUDZWVlq/pb3+LlW3ebENq8dTVWTy39jJtmj3JJ7O8tm61aBnxxD+7hdvxcQxD54jwWPieWL5pU/kOEbaHoUMt4G4+La98u6SHIa6eeCKzvGGhMMU7kqTzqHRrrzi4eDkdThiTrxBUVsKf/5z8WDRq8f+uucbWF0WjyS0vxyzL3r3NyaZPn5jQbNtm1saoUfDd7+Y3KkZIKDZXX52fJJWffQa33tr+4MrxHpWJltekSe3uZsaEi7ubm20h97BhsSggw4fb8Ht82ciRsNdecMABnWt5gA8b5oAPG3Zempvtpvn002YhhpEU4rfDBa0ffGBze4WiZ0+LrVhTE4stGG4nK8tkOxKxG+6BB1qMPpFYqpJIxOILJrrJh16GixbFboqZ3jYSo2VkQr7EJpFJkyx2Ylsks7wK4VEZDkVGo2Z1/v3vtr1pk0Xvr6hILkaffWYJPLP5GwwYYPOUxQ7y68OGTlnQ2Ahf/ao5blx2Gfz0p8Vxp25uhl/9yjL3hlZfKEDx0dkTt8OhzXCuq7bWYvyFN7G2FgDHC0OmiFiopsMOK26Cxni3+LffjkU670gyFZu2EMn/vFRoHT75pAnPkiU23PnIIzYEvnChWawbN1rE/2T7n36aOsJIOmcnsAebyy6zz5ap5dXZgvy65ZUDbnllRrogteHC3Npau0HGR07IN5kkqcyGgQPh85+37XjRy2QfzMvykUdaeluWGnffbeIZ3h5qauA737HAxYW2vAotNmFSiHDILZnw1NbG+hsKztq1tui7rs7eZ83K7zBzJBKzuNJZXuU8DFj2mZTLHRevzAgtr0cfbbtu/HovEZu7amqym/5FF8EVV+RuteVqeSUK0IYNJrT54IQTYoFoS5F4y0uksCGVEudwhg+38nhxSSxbvNgsl1Wr7O8SWjOh6FRUWFLOsKwQYpOMSAROPtkEry3Lq67OrCGR8hakbHDxKjIuXtmxdatZK/FhikTsBpNpor4BA+wmFmb4LQaJIgid1/LKlmTzYpDccgjLhg6Fl182i2jOnOL1PSReeNqyvAB23tnmD0PLstAi3xlw8SoyLl75IRq1dBznn9/yiTje8mpqinmc7b9/ZkklnY4ncZixIzjooFhQ3rYsL2gtNuDC09G4w4bTKaiosIyx556bus6OHeYhVV0Ns2d3XN+c7Jg0ySzTXCyvNWtsMXL8/Fi6cz/5xP53PJ5i18Ytrxxwy8txHCd78ml5eUATx3Ecp+zwYUPHyRBVy9p78802RLZ2rc2R5GP4qqnJvCpXrjTnFFUbTuvf397bU5b4PmCADdENHepDb0754uLlOBkybx5MmGDRNkLuv9/iNh5+eEw02hKVZMfnz7fFwB1NGDnCccoNFy+nrIlGTTxuvrljPN1GjTLrSNU83F5+Gd57z17tZe+9bfFzR1peU6a0v9+OUwxcvJyy5t574eyzi3PtzZstfNHHH9sia1VbmF1dbftjxlh8vEwsr8GDYwkfHcdpG/+pOAWlsRFOO83mioYMsSgH3brlr/1JkyyA65Qptgg6fo1YfLT4sCx+v7ExFqYq3A7LoeU54baILZTeuNHE6403UvftxRdtmHHMGJ9bcpx84+LlFJTJk+Ghh2x76VK78e+6q71v3x4Tg/hI6b162bqf229vGTYqGRUVtnB53brc+7h9e3b1M412vn69zYndf78l4zzsMAvSGx+wt73bbrE5XZWi/MuLSF/gbmAEsBg4TVXXJ6k3GfhpsPtLVZ0RlI8BbgW6A48BF6uqpmpXRASYDpwEbAO+raqvB20NA24EhgIKnKSqi/P9mbsqM2aYODz/fExgli1r+7wPPoAHH7SwUt27p6+rasKnCnvsYUFdN2yIHQ8jwq9f33r/gw9sDqvQvPuuvQpBJALTpxembccpVYqySFlEfgOsU9WpInIZUKeqlybU6QvMBcZiovIaMCYQoznARcArmHj9XlUfT9WuiJwEXIiJ18HAdFU9OLjOLOBXqvqUiPQCoqqa9tnaFylnTzQKt91maVFCayuV5bVjh80D5UJ1NbzySuZ5i+Jd1PNlDXXktlteTjlR9rENReRD4ChVXSEig4FZqrpnQp0zgzrfC/b/DMwKXs+p6ucT66VqNzxXVe+Mvz5QB9ygql/Mpv8uXoWluRl++UsbNtyxo23LK6RbN7jxRhg3zuPTOU4p0hliGw5U1RXB9mfAwCR1hgBL4/aXBWVDgu3E8nTtpmprV2CDiDwAjASeBi5T1VaxzkXkPOA8gGFhQiCnIEQi8POf28txHCcZBRMvEXkaGJTk0BXxO8FcVd7NvwzbrQQOBw4AlmDzZd8GbkrS3g3ADWCWV1476ziO42RFwcRLVY9JdUxEVorI4LjhvVVJqi3HhvZCdsWGDJcH2/Hly4PtVO0uxxwyEs+pBN5U1Y+Dfj0EjCeJeDmO4zilQ7EC884EJgfbk4GHk9R5EjhOROpEpA44DngyGBbcJCLjAy/Cs+POT9XuTOBsMcYDG4N2XgVqRWRAUG8CkIdYCY7jOE4hKZZ4TQWOFZH5wDHBPiIyVkRuBFDVdcBVmMC8ClwZlAGcj7m3LwAWAo+naxfzSPw4qP+X4HyCua0pwDMi8jYgwXHHcRynhPF8Xjng3oaO4zjZ4/m8HMdxnC6NL210OjWJi5DXrYNTToEzz7TQUo7jlCcuXk5J0twMv/61BfTdsAHq6iz8U7LtiorUSSEvuQT+939blt13H9x0k52bKtp7LttQ2LZceB0nhouXU5Jcey387GeZ17//frjhBhMxiIlbnz6Wt2r7dovAsWYN1NfDc88Vpt+F5v77LRTU6acXuyeOU1xcvJy80twMV18NTz9tKUNELEp8/PaOHXDhhXDWWaktiClTLDVJJpbX8uUwZw588om9ikFlJUycCP36Fd7ymjSpOJ/RcUoJ9zbMAfc2TM0118Bll2VWd489TMx69YItW6B3b7tJb9mSvCzcbmhomRcszKZ8yy0WLR5ailtIpmXZHqushB//GMaO9ZiKjpOOzhDb0OmkTJli4pLO8nr7bUtH8tFHuV+nrs4C8MYLG5h1EpbV1cHXv+5zRI7TGXHLKwfc8mofjY0mKKtXp7ay0lleb74JmzZldi0RuPNOnyNynFLALS+nrKmqMo+/XGlogBNPtPm1RMsLWlpe55zjc0SO0xlx8XLKjupqeOaZYvfCcZxi4jMBjuM4Ttnh4uU4juOUHS5ejuM4Ttnh4uU4juOUHS5ejuM4Ttnh3oZO2dHUBBdfbK7y6UgMt5QPRCxWYrIgwI7jdBwuXk7ZccklcN11xe3DQ/KjqCkAAAvjSURBVA/FggBD/oQyErG1aWPGeKgpx0mHi5dTdkybZtZXMSyvtWth5kx49VV7FYLbb4d//hNGjy5M+47TGXDxcsqOykr405+Kc21Vi3R/880WEDi+PJ+W1/77t68dx+nsuHg5GdHUZGlM3n/f9isr7Sbb1YLeilj0+LF5ic7mOE6uuHg5GXHJJXD99S3LnnkGrroKBg2ym3ptrUWLT7a9YYO9b90KjzxiIZ4cx3FypSjiJSJ9gbuBEcBi4DRVXZ+k3mTgp8HuL1V1RlA+BrgV6A48BlysqpqqXRERYDpwErAN+Laqvh609Rvgy9iygafCtvL+ocucadMsIG5oeS1dCh9/DB9+aK9sGDwY9tuvpaile0+WiHLjRth7b/jd78wKdByna1GUlCiBYKxT1akichlQp6qXJtTpC8wFxgIKvAaMCcRoDnAR8AomXr9X1cdTtSsiJwEXYuJ1MDBdVQ8WkUOB/waOCC77L+ByVZ2Vrv+eEiWWMTkMkJup5fXyy7ByZf76MXgwjB9vQ5cHHugu7I5TynSGlCgTgaOC7RnALODShDrHA0+p6joAEXkKOEFEZgE7qersoPw24BTg8TTtTgRuCyyq2SJSKyKDMVGsAaoBAaqAPN5aOy+RCPzsZ/bKhjCX19q17bO8li2D11+HFSvgwQet7fvvt35dmvif5DhOp6NY4jVQVVcE258BA5PUGQIsjdtfFpQNCbYTy9O1m7QtVX1ZRJ4DVmDi9UdVfT+3j+RkQntzeYWomqv65ZdDnz4xy2vKlPa37ThO6VMw8RKRp4FBSQ5dEb8TzFXlfewyk3ZFZHdgL2DXoOgpETlcVV9IUvc84DyAYcOG5bu7TpaIwLhxntfLcboqBRMvVT0m1TERWSkig1V1RTB8typJteXEhgDBBGZWUL5rQvnyYDtVu8uBoUnO+SYwW1W3BP16HDgEaCVeqnoDcAPYnFeqz+Y4juMUnmKt0JkJTA62JwMPJ6nzJHCciNSJSB1wHPBkMCy4SUTGB16EZ8edn6rdmcDZYowHNgbtLAGOFJFKEakCjgR82NBxHKfEKZZ4TQWOFZH5wDHBPiIyVkRuBAgcNa4CXg1eV4bOG8D5wI3AAmAh5qyRsl3MI/HjoP5fgvMB7gvOfxuYB8xT1b8X4gM7juM4+aMorvLljrvKO47jZE8+XeW7UGAfx3Ecp7Pg4uU4juOUHR5Yx2k3jY3wjW9A374to6vHb3sSR8dx8omLl9NuJk+Ge+/NrO6DD9pi4vj0IakEL/FYSKqyDRssF1ZVVX4+l+M4pYuLl9NuZsywlCnpLK+FC+Hpp2HuXHsVirlz4dhjbTtVjq10ubdUYd06OOWUrpfuxXHKCRcvp91UVcE996SvE43CnXfCv/5l24WwvJ59FhYsgD//uf2f6f77LVr96ae3vy3HcfKPi5fTIVRU2LzYN75RuGs0NcHFF1vEe2i/5TVpUuH66jhO+3DxcjoNlZXwpz8VuxeO43QEPqLvOI7jlB0uXo7jOE7Z4eLlOI7jlB0uXo7jOE7Z4eLlOI7jlB0uXo7jOE7Z4eLlOI7jlB2ezysHRGQz8GGx+5Ej/YE1xe5EjpRr38u13+B9Lwbl2m9ou+/DVTVJeIDs8UXKufFhvhKqdTQiMtf73rGUa7/B+14MyrXf0LF992FDx3Ecp+xw8XIcx3HKDhev3Lih2B1oB973jqdc+w3e92JQrv2GDuy7O2w4juM4ZYdbXo7jOE75oapd7gX0BZ4C5gfvdSnqTQ7qzAcmx5WPAd4GFgC/J2bBJm0XkKDeAuAt4MCgfDjwOvAm8C7w/TLq+2jg5aDfbwGnl0O/g2NPABuAR9ro8wnYkogFwGVJjncD7g6OvwKMiDt2eVD+IXB8W20CI4M2FgRtVrd1jRLv9xHY/3YTcGoWv81S6Pt/AO8F/zfPYO7d5dL372O/kzeBfwF7l0O/445/HVBgbJv9zvSfqjO9gN+EXyhwGXBNkjp9gY+D97pgO7wxzgHGYzfIx4ET07ULnBTUk+C8V4LyaqBbsN0LWAzsUiZ93wMYFWzvAqwAaku938Gxo4GvkEa8gAiwENgt+DvNI+FGAJwPXB9snwHcHWzvHdTvhv1YFwbtpWwTuAc4I9i+HvhBumuUQb9HAPsBt5GheJVQ378E9Ai2f9DWd15ifd8p7nonA0+UQ7+D/d7A88BsXLxS/sE+BAYH24OxdVuJdc4E/hy3/+egbDDwQbJ6qdoNz012/biyfsAS2havkut7UD6PQMzKod/AUaQXr0OAJ+P2LwcuT6jzJHBIsF2JLc6UxLphvVRtBuesASoTr53qGqXe77i6t5K5eJVU34PyA4AXy7TvZwKPl0u/gd8BXwZmkYF4ddU5r4GquiLY/gwYmKTOEGBp3P6yoGxIsJ1Ynq7dVG0hIkNF5K3g+DWq+mm59D1ERMZhT1gLy6nfbZDJ+f9XR1WbgI3YQ0i6z5GsvB+wIWgj8VqprlHq/c6FUuz7uZgFXzZ9F5ELRGQhNipxUTn0W0QOBIaq6qNt9Pf/6LQRNkTkaWBQkkNXxO+oqoqI5vv6mbarqkuB/URkF+AhEbkP+Btl0HcAERkM/BWbq/qHiJRFvx2nLUTkm8BY4Mhi9yUbVPVPwJ9E5Czgp9hvs2QRkQrgt8C3szmv04qXqh6T6piIrBSRwaq6Irj5rkpSbTk2tBSyK2bOLg+248uXB9up2l0ODE1xTtjfT0XkHeDwcum7iOwEPApcoaqzgbLod4Zkcn5YZ5mIVAJ9gLVtnJusfC1QKyKVwVNpfP1U1yj1fudCyfRdRI7BHrqOVNX6cup7HHcB/1sG/e4NfAGYJSJgD+4zReRkVZ2bsueZjEV3thfw37Sc5P9Nkjp9gUWY40BdsN03OJboPHBSunaxcdx454E5QfmuQPdguw74CNi3TPpejXli/aicvvO4ax1F+jmvSsxhZCSxSed9EupcQMuJ7HuC7X1oOZH9MTaJnbJN4F5aTmSfn+4apd7vuGvdSuZzXiXRd2yeayFp5nBLuO+j4q73FWBuOfQ74XqzcIeNlH+wftiNdz7wNLEb5Fjgxrh638FcOhcA58SVjwXeCf7B/0jMbTtVuwL8Kaj/dviHAY7F3HHnBe/nlVHfvwk0Yi654Wt0qfc7OPYCsBrYjo27H5+izydhDxQLMesS4Erg5GC7BvsxLsDEdbe4c68IzvuQwDMyVZtB+W5BGwuCNru1dY0033Up9Pug4Lvdij1xv5vhb7MU+v40sJLY//XMMur7dGz5ypvAcyQIUan2O6E/s8hAvDzChuM4jlN2dFVvQ8dxHKeMcfFyHMdxyg4XL8dxHKfscPFyHMdxyg4XL8dxHKfscPFynBJCREYEi9Xjy34hIj8XkTfjys4Uke0iUhXs7xuEGXOcLoGLl+OUB9uBYSLSO9g/FHgfW1Ab7r9UjI45TjFw8XKc8iAKzAUODvbHYIuwDw32DwVeBBCR/ycir4rIOyJygxifF5E5YWOBhfd2sD1GRP4pIq+JyJNBmC3HKWlcvBynfHgROFREemJiNouW4hVaXn9U1YNU9QtAd+D/U9UPgGoRGRnUOR24Oxh2/AMWwmkMcDPwqw75NI7TDly8HKe0SBXyRjFxOhQYB7yqqguB3UVkANAr2Af4koi8ElhWE7AYdGCJAE8Ptk/HMtnuiQVFfSqYU/spLYMgO05J0mmjyjtOmbIWC0ocTxiweDYWM/Aw4OXg2DIsWOrLACJSA1yHxYZbKiK/wGLTgYnVvSLyAJZBZr6I7IvFHTykcB/JcfKPW16OU0Ko6hZghYhMABCRvsAJwL9UdTOW5O8cYuL1MvAjgvkuYkK1RkR6AafGtb0QaAZ+hgkZWEDVASJySHC9KhEJLTXHKVlcvByn9Dgb+FkwjPcs8F9xQ4IvYpG4w0y1L2ORul8CUNUNwF+wCPxPAq8mtH03lhHgnqB+AyZw14jIPCwa+aE4TonjUeUdx3GcssMtL8dxHKfscPFyHMdxyg4XL8dxHKfscPFyHMdxyg4XL8dxHKfscPFyHMdxyg4XL8dxHKfscPFyHMdxyo7/H7xAjiZCMZ7CAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "infilename = os.path.join(data_dir, 'uv.fits')\n", "\n", "infile = fits.open(infilename)\n", "data = infile[0].data\n", "\n", "uarr = np.array(data.par('UU--'))\n", "varr = np.array(data.par('VV--'))\n", "\n", "fraction=float(input('Give in percentage what fraction of the data to be plotted(e.g. for 10%, say 0.1):'))\n", "plotlen = int(len(uarr)*fraction)\n", "\n", "plt.plot(uarr[:plotlen],varr[:plotlen],'b.',markersize=1)\n", "plt.xlabel('UWave')\n", "plt.ylabel('VWave')\n", "plt.show()\n", "\n", "infile.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Virtual VRI software - A Virtual Radio Interferometer application\n", "Below is the friendly virtual VRI software avaialble freely online at https://crpurcell.github.io/friendlyVRI/.\n", "\n", "Run this software and follow instructions below.\n", "https://crpurcell.github.io/friendlyVRI/HELP.html\n", "\n", "Use at least three antenna arrays (GMRT, VLA and ALMA) and using 2 model images and create the radio image.\n", "For the same observation, use different hour angles and note how image quality improves with the increased hour angle." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "os.chdir(os.path.join('data', 'friendlyVRI-master'))\n", "os.system('python vriTk.py')\n", "os.chdir(cwd)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Image Files\n", "In the code below, you can fiddle around with the image files for 1 hr, 2 hr, 5 hr and 10 hr observations. They are pre-made FITs file. You can zoom-in, zoom-out, check source structure, check the peak flux, pixel positions etc." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Enter the fits file name (these are 1hr.fits, 2hrs.fits, 5hrs.fits, 10hrs.fits): 1hr.fits\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAEKCAYAAABewe3GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztvXu0b1dVJvjNe24SHgohpGQgYBMlbTeKFlQEHJSWTQQDbRlGN1hBW1NWulOjC2wfVaMq2DVQKasUhwOUKnxkCFaklFdEiZjqNC+rqi2NCYJAwBTXoBAEQ0gIoIbce87sP/be964zz/zmnOt39nn97v7u+I27H2vNNdfaa33rm/P8HqKqWLBgwYIFe49jB+3AggULFpwtWAh3wYIFC/YJC+EuWLBgwT5hIdwFCxYs2CcshLtgwYIF+4SFcBcsWLBgn7AQ7oIFCxbsExbCXbBgwYJ9wkK4CxYsWLBPOH7QDuwFzpXz9EF46EG7sWDBWuPzuPduVf1bu7Hxbf/TQ/Uz92ym5d7z/i/epKqX7aatw4C1JNwH4aF4mly6/aLI9nPV4Vr70WZbxsKr42Eq49nz2puuteV361fUdtuubTPzwevXqh8Pr/Srt45no/LMWtvVcayC+dwznxiqddvnZuuw68n9d+j1fx43nuMz92ziD2/6irTcxqM/cuFu2zoM2LOUgoi8VkTuEpEPNtcuEJG3i8hHxv8fMV4XEXmViJwQkfeLyFOaOleO5T8iIld2OHDmBfBJZsvZsraeJZy2fmtnWkj2ZW3Zc1bOszMd2754dWy/rY+ez63NaFPwxsCOqVcu66fnczamFt5z9saL1fXaY/1jbUdEZsfQ1p3qs373bFqsDhvHSBD0bnzMPQBbhX/rgr3M4f57ADYEuAbAO1X1YgDvHM8B4DkALh5fVwP4BWAgaAA/CuBpAJ4K4Ecnkk6RLUjvPiMztiiYzRbt4mZ1PHhk1R5HfrO2PR96+uT1LSNNpoaZv9EGFpG553fUl5ZMPbLzxqzdKDKyZ6S2Sl2GaL4xm+y5MPu2rr23SygUJ3Uzfa0L9oxwVfU/A7jHXL4cwHXj8XUAntdc/1Ud8AcAzheRRwP4NgBvV9V7VPVeAG/HThLnyBaoVYReOTY5I3ue8vTatD62ZSM1nJGqtcEW+iqLJlLa2ebB1HjUTmSfPa9Kv6IybJOxBJSRP5tT1gZrszpvsg3V9sM77xUSM+JsUrj7ncN9lKp+cjz+FIBHjcePAfDxptyd4zV2fQdE5GoM6hgPwkNqxMlUkCUIFt56dlk7UXhpF7W1N7Ub1a2GfXYzsO1N7Vhlz8aAhaJMFUcRh+ez9wwirKLc7LVoDL1nZccmiziqfkebY9aH6bwl/WyTaf1nfvU+jwQKxeYekvlhw4H90UxVVURmG2lVvRbAtQDwMLlAd0wgT8l65OipwyxMO+OET6xsEjOSae219z2fq0rYs23JISpbucbuWfK25awv2cZYDYmtbetTa7uaovDmh3fdtsX6ltWN/O3d1Fo7rGwmQtrzmVbv1lyGjgD2+324fzmmCjD+f9d4/RMAHteUe+x4jV3PkYVyu7FhSc0Lf71Qn9noASOPTNl6G0uFkHvrZfDSA1k6xdZndmwZa9tuXD3E3dryrnu+ZGkT5rPXXoUYq+o+GjvWri0/m8IFNqHpa12w34R7A4Arx+MrAby1uf6947sVng7gvjH1cBOAZ4vII8Y/lj17vJbDI8L2/7acV75CLDYMzVILtp5t1yOiyEYWnttrUfjpteW1scpCq4wj26gsrL/VMfb+t9cyZRgRVbYxs/GN7EQRmGfD22ja69mmnPnVszl1YAuavtYFe5ZSEJHXA/gWABeKyJ0Y3m3wUwDeJCJXAfhzAN85Fr8RwHMBnADw1wC+DwBU9R4R+VcAbhnLvUxV7R/imAP5RGPKJFM+ESnYHGhE3F66o/UhKx/51qJK/j2o9Cur77XPNkZ23yMoRlJ2TvSmGRjZeXbaez3RDSPN1lbWj4rtyubSe28FKICTM9o77NgzwlXVF5Jbl9oLOvyw2ouIndcCeO1KTkQkE+UVM9Kyi8JTi/b/aFFam4yoIzKyC5HB87dVl9XFViWNVctYP6INzCPaiLi8Nqrqr+IrizC8euwaU9Js7lXaq6Lq/wzQNUsZZFjLT5qliPKqllTtYmYLmS1eW9barSws1odIuXoqiC3Q3eQZWdusTjROFSVuo4eKbxXi8u5H1yOCtX5FpMjUbzY2kV3veR9WFanA5iF1bS9wdhIu4KvcSOWwcLUt6y3ELJdmbWTqubXpqcDIN6+dihL07Nm6LTLy8zYaL5xnIX5GItG4ROSbqdxIXUbktopi9uxk6YbMj1Xa7C3TCQXW6F22OdabcCvpAZZOsJM7UlVZaOeRHbsWLRiPzLOFzvpk/bA2vfLMlwp6iIApcttn9myjcYuQEWVEnJXIwT4LFiF5NnqigkpqySIj2nZcV7HPjWMTnRvSEcZ6Ey5DFIZHZexC9xZMllpgypItbk9FVRZmRs7eookIJWq7Aq/dqM+MjKIIwpb1bLI+RymKipLONgF7HD0L1sdKlONFYtGmy/rD/M9sdEIBnNSFcI8+MvKL1EZFGbN2MrJloWALj/iZos0UqVffKxsRjtevyqLzFnBEZN5GZstF5OTZmspEBJXVrxBR9Px7Nng29yobgy1nfbXRQQ+R7kE6AcD4PtyFcI8+KhO+BVMS7XmkSpjCactnZJ2lAFh/vEXuqcTIRm9YWgkr51ic2fhUNgnPplcm2rAy5cpsZ+dtHbuRemWyubxq+mSy7flUWT+7wNaicNcIWVjf7votAVo1yexF5JiF74xgvbLVcD4LuSv3LKKNoEK8FduZT5UIpPWPqbgW9nlXCDzaxCpzxdvYWT/s8SpRRuSf55tX1j7vVdpmLmFRuOuPjGiyReUp4UgReAvFI3Wmkllb3nmVACthqb3vpSUiRVTdIKJ+e2TMIhEPWTrA+5+NTVXVR/2ICNSzkUUvbDP3NvvI/97yM0Eh2DyLfulr/Qm3qvKyhRmpEo/wmELOFpPXbqSKvDCQ2expv72elWO+V1FRti25V8Phim+VTS7yi/kzXY82ZxtZeXW96KY6Z1s7q0ZIUb9m4uIlpbAO6FV5URm7KKrhqmcnuxbVZeTACJuRiKdUGblmoSUra/vXM16ZwmdlrC025qxd5p/X9pzzqz1mz9fzk23qts6cKjXaHFcxB8EDujGbvcOO9dbyLalUkZGCF163r7ZdpoJsKM4WTFvW2mL2bF0PEWFbv70+2Ws2ReL1mdn3wNSY12cvrRH1x/OHXbeb0nTdEl7Ulkd67TXv2Xlqns01rw/VaKwC9rxnggLYwrH0VXNVLhOR28ef6rrGuf/NIvJHInJKRJ5v7q32U16dWF+Fu2o4zEK7ivJrJySb9B55ssUTqWKmTqOFyPqUlY/sRH57KivaXDzb1UjBHldInalA7/mwsNyLDrz22WbgoVLGtuddj+ZVVJ/Zm1MpN5jjj2YisgHg1QCeheGHCm4RkRtU9UNNsY8B+IcA/pmpO/2U1yUY9oD3jHXv3bVjButLuC2ikJuplWzSRvYjko3IpYoKwVZU3VyhYaTkK4rI28w8H6t2bL3IZ8+/iFQ9+566j8q1x1Hfowilkt6opEsq8PycCaqCTZ0l0H4qgBOqegcAiMgbMPx012nCVdU/G+/ZTxOf/imv8f70U16vn8OxFmcH4UYLhqkTO5GZcmUT0ZZltlgI6vWhB9FC3o2dFuzaqoubER8j3VZ1ZpFEpb1Mye+2v0yJR5FC5kuFlFdBtIky31bEVk3hXigitzbn146/8jLB+zmupxVdKP+U126xvoRrCS1bgFlYmOVDPbXppRo82549C0/5RMojUvX2GvMjyqNmRNGGodHiZyqv0i5T09W+R/bY2M6xITLV2573PJfdkmuljT3C8EezEg3draqX7LU/e431/aNZu3h7QlF2nYWMtr22zSwX2JtD83yx/YwWKitrbbX/e/22Cj1Sv9aW50tbzh7b/mcbH0NVobft2OfV+u49A+++HR82lmwcbH22sbPn1NqpjkU232ZUt4rZ/mi2+s9x7a5uF9aXcCuhT6asvIVvF799WTBFkqk5puI82MXt1W3PGQlWwvBMRbbHrI/eWNkxrDy/aJPr2dA84qy0xdQomwfRPGNEXPGn9ck79+Y5e/Zss91Dxbupkr4KuAXAxSJykYicC+AKDD/dVcHqP+XVifUlXE8pZIS6yuSKFE+26D3Vaa9bPxnJM5UUpTG8RZkt3IgoPeXF+s3Q1qsoQ68MU3wZMUbPyvbTtmV9a+vY4+w5eM8gUste+x5hZoTslfH6NCOmT5plr9SO6ikAL8ZAlB8G8CZVvU1EXiYi3zF0Qb5h/LmvFwD4JRG5bax7D4Dpp7xuQc9PeXVifXO4EyoThIXW0z2rBLL69jjziYVp7UKwBNO2ya5bf6uKyRJLpB4zZcxIkY0x89cjRzbezJdMlWdlI1/YvYp69srbeWfb85639SHb9LINxjveA2zN8y4FqOqNGH4fsb320ub4FgzpAq/u6j/l1YH1J9wJHnEAfMKyCVlRGNY2qzeV8Qi19dmrx9rvCYsjXxlxeP5ki9ciUnqVZ2PhjV20Cdi6nl32LLJnXZkLXvus7VWeSxVR3Wi+zggFlu9SWCvYENUL1abjtk6LyiTzFo63QCJSrZB8ZDvra3u9ve8tajsuXijN+mSR+WTbsX2ukmeLSM1mG4XdnCNf2NzoIdhonrBNl/Ul23B7N1qv7CrPg5mE4ORZ9NHe9SfcHthwzhKYh2wSeoszCscrobZHrJ6vLXFYm0zFe7YykmQbFFPg0ebWgpFIdbFHY9v6WbHhlc3636JnY6+MkRexZWmFaprAI/VoHHYBVcz1wYcjgbODcNmi6z3ObFsw4mGLyRKsd7yKj1OZSG1HitUj7qjd3jA3K5dFCcyXKFzP2sxSG9O5JT1W1x5nbTLSs/ezDSu6xjbfqJ2ZlG3TQPWDD2uBs4NwvUlbmTiV0Ckqs0oYbNV1ZeHb9pjq8dpiqYKKOmRqOFO0FRLySNGSwqrPpqps23azMj3zyitX6Qu77j2b6BnYNj1BsKo/nVAsCnf90KsGe8r2KNy2TqQore0oRcDsWx883+wCqypgW8bz1fPbnnvIFL5VZJGKjcY0y1W25diGHREUI7jKmLZ+VhSstzlkKSG2oVlEc3Om7MLyR7N1Q6/6mAsVZWTLsRSDV96GgCzsbO9FxMjasOQVLdr2fmsn23S8umxTyp5jdaNkRMUIkW1MHsllGzGLYjJ4c7mywXqI1H+2Ic0EhSxfQL4WqJKsVy6q64X5q6YOVkkRRGVYbrNty7bJ1KO1nyFKG7ByXhsVO9VnlilP1m6UTsief0XpsvYqBMzsW/Vt63ibB7sXtTWzaFEAJ2vfpbAWWN+eeirFC1MjdWfr23veOQPLq3mLNyJiFiJ65719ixSzRyrZRuClPyLfog2jvZ/1y7sWhd2r+J+hZ2Np7/Wq3ozss6igulHZNmdTubL8iORaIVOGEyqKuCfH5dXNlEy0ODLij8jaIy17va1TUbmePe++R26sXVaWbU6trYxce5WZN94RiVZCcXa/7cMcKjLakD2wdEokLkQwRw5XMd8nzY4C1ptwswnUokrM1fvZhJ3KRHYi4rOEZhdrte9Z2NizyWQhvG23l6QiVej50KZZehWZR6SrbMbsnPnInlnPmGfP3RsTb0OMfJhN4S4/k76gF1ko7C2ATM1keUTv2LvXKsQIEZmzspHqzgikYs/a8tphtqxdz0Y0Jt69SvhdGWeWqvDGwTvOoqXeSMqOj1ffVfxxVytQlUXhrg0y9ckmpreTV4ixEoaxetaPbAFVQt5MqXroCW9XIRoPdmPw2mSKOBoLu7nZNqv+rhIZMVuZP5YMe9IL0QZvy0WpELspzJHmYC4Dy0d71woegXmTmamznpRCpFA8EvDsRCQWka2HSqhqr2WqL9qUPAXWXqtsWt41Fi3YPlo/bdtVlW/70Us4UVmvX6xdT0l7Sp9tLJkflflt7VQ2oS7M9ptmRwIH0lMR+SERuU1EPigirxeRB41fHHzz+BPHbxy/RBgict54fmK8//hSI94CZQvb2+WzcnaBeyToEb21EbXp2Wj7ZP2zdexCrBBHpDQzxWzvty9vzFgfMmJk7bft2mOvjh2nyjP27nskx+x4vnjPxvO1JUlvLnnI2m/75W2Sti9Vki5CAWyppK91wb4Trog8BsD/BeASVf1aABsYvp395QBeqapPAHAvgKvGKlcBuHe8/sqxXKWh7ROJKbhMebX1I3VqFRVboCxk9Ai4Pc/8zBZ6RnxV1ewhGhvWPiNsb8No7TIFVyECL2KICMnzgZG5tem155XLbFifvL6w8WcE6Y0r65u3cc+MOb6A/KjgoHpyHMCDReQ4gIcA+CSAZwK4frx/HYDnjceXj+cY718qUmCDaGFW4E3iLFxjRBz55S2KTF1lfckImC24uVBZlGxDilRwNJbMRlSfwW5MUURh2+4Bm1/R5uiVYT4wFcs2hUjhW5szYfqk2aJw9wiq+gkAPwPgYxiI9j4A7wHw2fFnMoDtP1N8+ieMx/v3AXjkSo1XwukzjvLjKLRk9Wzo6v2fKTyrEC0hRMrQu+f1e7ewZG/Hp+JHNo7tcc/mwRRl1p+MeDJVHkUAHrFXNo2ojB273sih9XmOjSXBTD8ieSSw7380G3+k7XIAFwH4LIA3A7hsBrtXA7gaAB6Eh2xXjxEiFWOvZ2RlJyhTZJFPlUXkKRtm15Iz820usPHynkW72Hf7nLw2K3ba8WHk09bNfGbjzBRlFqazZ83K2jLsGVTHqbdeJ1SBk1vrQ6gZDuJdCt8K4KOq+mkAEJG3AHgGgPNF5PioYtufKZ5+wvjOMQXxcACfsUZV9VoA1wLAw+SCYbZ4KvBMhXjytIswsuHVa/+PFlV03fPVTnyvLrvWQ0S7DSUZudp+MCUW+buKSrVgJNaTQsjKR2W8yMT65/ls61b99kiXze02gmpts37tEkNK4ewh3IPo6ccAPF1EHjLmYi8F8CEA7wbw/LHMlQDeOh7fMJ5jvP8u1eLTbhexdz0yU1WE1r7Nk7WTuxLOR+G451fmX3aPlWf+ZqioNS/UjuywsLu9Z8fL8yUK7at98dI/9npbtt1E2jnhEbDtT+u3bdf6UInQ2nP2DGyfqvN2F9gcv08heq0L9l3hqurNInI9gD8CcArAezEo098B8AYR+Ynx2mvGKq8B8DoROQHgHgzvaKg0VDvP1CcjOm+iWvssnKso56iuh8hWZfOYQ9m2dpg/FTVoy3m+VYgzU8mM9LyyzFcvounxZbpnyTmy5fmcbVx27jE7bdmozkxQYK3+KJbhQD74oKo/CuBHzeU7ADzVKXs/ht+RX7Wx4f8oVM0WlEd4FZFdWRC2rbbN9lpl0nuLvoeA2Rj1LrhorKu2elIjTP1m16xdb1zYxmw3Fs9u5DODjYra6/baqqq2Enntth9lnF0phfX/pNkEbyJHk8hTWJXJydSop/hYuUwVRgTj1ctSKBkZMf8jROMa2Yr8jNrKyK8di1617JGf52/kN5t33hxhpBspUxYReX5VogNbvveZdGD5TbN1Q6ROvcUTERkjblum0nbkK2u3shAzhcNsZ5hL4Xj+eeOzm/ayTcbar0Q/THlGsD5U0jttuZ7n7vnJfI7GJtp0Z1a5qsDJreW7FNYDrfKx+TFLhJVd35JCj+rNfIx8j3ywdaptRuetDebfHIsvU81ZHz1EpFlV1t44MsKKUhzWNhtH1o+IqFs/WH3vuK2TRUhe/Zmx/MTOusCbsJYcs/RCu6jYwvMIvBLORb6yslGOMeqDRZaSYIozigZY+xEsEVXC4ooS9UjFa8+7HyEL/z1/KmPt+d3665Vr77NnnM0fJkTY2Eft7QJLSmFdUJ2Itqy3iDxCba95xJcpOK/9KFURkXXPAulVgdae9TOyy8p6ddp2orSDLd+DTLF7ffWetzdekSJlZOf5wa6xTYT1IUtBsHuZKJiRdBXLuxTWA9FiZGqX1WeKqL0WKeRqSsH6lIXanoJpbXgLJCPKqH+rLrRqnaxcRvCM/KKoxmuD2fHa6mnbtu/ZyaKdSAx4hMx8Yz70zO+ZcDa9S+Hs6Gk7ET311BIVUxFWzdpXe6+16/nR2swmMlvwFWKOFiDzPepfS7wVdVlVoNUx7bHHnnU2VrYtr8/R2Hubp+e/t4FFm6+t7/WPjZud196YeNdtGebLLqAqOKXH0lcFInKZiNw+fpXrNc5996teReQcEblORD4gIh8WkZfM0jkHZwfhVkPRdhFkYbB9tdcjBdpe99QOI7NMhUTEnC3U1qbXN7s5sMUZbTCRn1F71jeGzDcPTM217XtzxyNjz59MWdo2mfqOIh02ZxjZr0qUTFTMgDm+LUxENgC8GsBzADwRwAtF5ImmGPuq1xcAOE9VnwTg7wD4x+Xv3e7E+hIuU6zexPMIM1O6ng2mWiw81RaRSmWx9BCNV4eFuhHh2MVXab93s2CRBbPbq4yjzZhtqOxa5E9rv90cPH/sM2FkHc1NW6ad12wsK0S6yjwLoJjtC8ifCuCEqt6hqg8AeAOGL8lqwb7qVQE8dPyulgcDeADA52bo3g6sL+G2yBQGK5+F3l69LCxjC9WSfFW1VNR7T5+tfWazJQ1vbKJximza82gcmE9TmYoyjojSG4eIpNgm7bXp9TUSCLZ/Xp0MXpRSmVNeP2Yk3ZkI9/TXuI5ov+J1RxnzVa/XA/grDF8X+zEAP6Oq9+yuVz7W949mETxyYaTslV21Te84UneRLbv4VvHPI7JIVXlEUfErate25208LDrxSKBXcTNEpMnSBe399thuTNam13dvHKvzxhJxNMa2Huuv58MM6Hgf7oUicmtzfu34DYFz4KkANgF8OYBHAPgvIvIOVb1jJvunsb6E64X2XqrAW1iVtMBUzju2tqL7Xll7n5GfRUX5ZeW9cWnb90iVLWzWVkbiHumyftn2Wb884vYIPSPGiDzZNUbIUR/ZNdamN9+iuR6R7qpR0ooovg/3blW9JLg/fY3rhPYrXm0Z+1Wv3wXg/1HVkwDuEpHfA3AJhu93mRVnR0oB4OGqPc/CY7aoWCgYkZq3oLwQsV3kzJdMeVbBFltF6XoLvBLy2hDX9jUrP/kdtWX9Z1GO939bxqrb9r71p23Lts36ZNuL6kQbcUSeXmrCttm2y+brDFAFTm0dS18F3ALg4vHHaM/F8K2CN5gy7KteP4bhJ74gIg8F8HQAfzJD93ZgfRWuBSOJSKmwicbCSU8xVeDZqICRYEZwWehYSQtkfaz6kNn3yrLxqvht22DkwyKiyJ9sLnlkzZQ3U6MscrNgY8EUcKvuK5vWjOQ7xwcfVPWUiLwYwE0Yfpj2tap6m4i8DMCtqnoD+Fe9vhrAr4jIbQAEwK+o6vt37ZSD9SXc6uRhCioLKyOFWV0UVYKOFq9X1yMRTwVFbVlERMQIgSnBzAfPp0q4XlWB1rZ3bompQmoe+WcbYTRPsvnFbEdzoW3HG7to3ljbM2DO71JQ1RsB3GiuvbQ5dr/qVVW/4F3fC6wv4QK1BcjIwO767fXpOJp8jJjbdmxZ7157P2oz2mCq5MZUlLXBSDezUfUlIrn2vn122dh4/jFSifrPiNbz21OUbEPsjRYYIbd1qpuQNx7Z85uJe3X5aO8aoBoaedcrKoPZixSeR+A9CjsLnTM15fnI+rNqmR4FxEjCtpHVr6p322bvxpf5YKOO9nr0bLLx7iHkbAPJ5gWzsYdYvrxmneCp0ywMm8CIcrqXKSF7XgkZKwtjFVQX3YRVyvT4nBFatilF6rf1p6rmbJmMmDyitfatLU8J2zRE1KaHihKPNsXqpmbtz5RaUF2+vGZ94IVuQKxwq6qpolazyZwpamun2pdI/VZDemaHoSc8r/rjpSy8+7ZMRSlnJBeRUmafEZzd8Oz/0fixsYpSOVF5z3+P/O2cm4lom8axufxM+hogmoDZ5G2Po8UUEWy2SDy1HZGKXRRZOiIj4lVSLhlWqceiBC8899qo5hwzZcYUcHTf85nNj6qSjNR8NldaP1l6JBqH6NqklLfZ9V3oxZLDXSewhWOvVUPXTGl59ZmSjRZtlMaI2o0USsXHVcH86lHUbXlGtLtNvUThdXRsxz4a3woxen4xVcns22s7CNHZ/Ns0QzVNUWl/RSiWlMJ6ozcUj2y059UwLlNsXrjJ7HsLqtoHRoq95GvbjcjA8yELmVkdSyysbiWNZG1HdiI/vPYqmw8jaK9NVt+mA7J5ZZHZbsvNmVbYiyzFIcbZR7gTovAPyNMDnppginU69trxymcEGC0opogzJWbtVFFZqFXFlvWdbSrZJteeZ8ququ4qaZuonFenQrbVDaG95kVMdt5WNsUdvvAu9WB5l8I6IFOdU5npfiW35S0IG6JFpOuBqbzILutXlu7YbTieoZLuYIuf+e49R7aJVEg+k1NZ+ilCZWO27bBN1vOrutFE86JV5t59tmF7fs8AXf5otkbIJnp0b0K0qCsLkRGzB9ZuplRtuShsrOb/WD+y9iNE7Xv+MD97Uwk95FBRrQxsHKP6dn6wTahqr0LKzHYUjVU2hRWxpBTWBb3hHdv97QSNCMIrz9qP1Fhl4UR9yOowcmALfy41nCn+VRSovc/Gk/WJhdteu5GdTK1mpOmlEVhqwdbzSNN7rt4xU7XevM+e3wpY3qWwLugJ5Vi+iy0mb/Ix5cVIjJXbLSqhtkf2nl+RAvPgqaSIgNgzsmrWtumNn+1vRJK2bNYntgGusqlN9zx/2P0oKpqDyNu6LILbMY6+Sz0YurYQ7tEH2/EtgdgyExjRtPc8sHrRBGfnVXh9ae1VQsFKyOop+6yO10dG6j0kxtRfZQOMSNzzn82Ltk1GqPa84rf1y+sPO2d9sPPSK9uj4mdUucvbwtYJ2W5eVZVVxdRLPHMpWk8RVvyIbDE/PcKsRhMVwp/qes9olc3AnlefUbbJemMT+ZkpSAsWTVTqVxR5Ze55AuV0+3n1CmbOUBxqrDfhskXbghFmNqkjgmjPvfYzdbSb/mQo8buNAAAgAElEQVT+VpQUa7Nin13LfK2q0Mx+RkJWnUX97w3VvXuZmmQ2mIqNwJ5VTzRg77NNciaWVAi2lncprAkitRERT5RCqBAOy/l5i6hH4VYJgN2zi79CpJ79qW5GOF5brP22TuZDlSSt34w4PFKL2qgoblbP1omeQRTaM9uRvxEBR6Ije9a7xFkkcNeccIHtky0iLEbAVXKyisa277XJ7Nt6Hjwi88owv5mPmZqJFPt0v0KcmSKu1ItUfrSBWr9Z+UxZV8tVyDC7l6Vs2v8ZojHJ/OtR2j1QYPmj2TqBhVhZnSwnGLXl1Y0WTpV8enJ+Uc7Rlo0UF1NGVWXH+sYIpEfhrrohRVHKlHLINiobMVSiimoKhEVYbJOrbhyrkmW0Wc2Bs0jirj/hTvDIxrvfk4aw5SqLZypnia6CiqLyyJEtSOuzp+qjHGPPZlZRY70Kt4JVNgtGqsx2ds3zxSvLyLS9nqWrvPaq89L2O5oHIrMR5dmkcA8kWy0i54vI9SLyJyLyYRH5RhG5QETeLiIfGf9/xFhWRORVInJCRN4vIk8pN8Qm1HTM1M/0as/bOlF5a8/64y3giNS8+vZaz2Jr/e1R7dZ/qwaZSm6v9ZKl9dv6H5XznlN7bJ+x9d1Ti0wB23bt9aoatH5E42XTGra9ajQU+WrXSCZaVoAC2NqS9LUuOKg/D/4cht+B/x8AfD2ADwO4BsA7VfViAO8czwHgOQAuHl9XA/iFcitWGbSoqJ3KpGRttuUtsXsKpkKuEZExIojQkmemTNtzb+zY+DL0Llg7hlnI3557xMmI1V7znpUlIc/X6NllSne65gkG+6zZPKkSNgMTHdX6VSgAlfy1Jth3whWRhwP4Zgw/WQxVfUBVPwvgcgDXjcWuA/C88fhyAL+qA/4AwPki8ujORnkoyxQBU7Ltfe/YsxHZzkiLvZhf0YYQEUWmbJi6ZzarYfsqYCG1RxJss61umuzY2t4NuVm1Go03i5K8eRRtRNkcactF62cGRNO8J0A4CjgIhXsRgE9j+B3494rIL4vIQwE8SlU/OZb5FIBHjcePAfDxpv6d47VtEJGrReRWEbn1JL64s1VGVkPlna/sSXvKJ0LvLPJ8qqiWSHUxMvIWaEZSM4eWu0Z1nFj6gJXxrnvzJmqHtWuPvbKVDTSK4NrzSBjYcr2Rz26ghdea4CD+aHYcwFMAfL+q3iwiP4cz6QMAgKqqiHQNs6peC+BaAHiYXKBhSF4J6XonWrsYqyFlldBbGxVf2/5lbbWqqkJGts5hhNf/Ft6GasfAU64sKvHSE9ZWlLawvmXtszoWVVLsKbdtvtSqxRAcpT+aiciDAHw7gG8C8OUA/gbABwH8jqreltU/CIV7J4A7VfXm8fx6DAT8l1OqYPz/rvH+JwA8rqn/2PFajoqSjMKy3cQ2FWUyHbcEbRcaUxoVhRzVz4hkaiNSzBnk2PZXdN2WWRVs/Lwy23xNNueorenY21S98ix6YCkKpnizyGTVOcPa26sN9ogoXBH5cQC/B+AbAdwM4JcAvAnAKQA/Nf6x/+siG/uucFX1UyLycRH5alW9HcClAD40vq4E8FPj/28dq9wA4MUi8gYATwNwX5N64GB5Pk+d2InIlBFrw05y2041lLM2PRKOiMESZ0YWkeryxsTzswLd2k6s07Xp3JJxe79kn/Q9GwdP3dvjtg1bz96zz4mRbiX68WxndjNbEawyt+MRzcHdQAE9Ou9C+ENV/VFy7xUi8mUAviIycFDvw/1+AL8mIucCuAPA92FQ228SkasA/DmA7xzL3gjguQBOAPjrsWyOiBhYqF0Jv7Ny0QJmZT3b3sJn4WslReH1xVOD7X3Wp9T2se2EacnWHjMbQI14o5A6IgeWPvDSK961ykZrCbOqMlm/Ih+8vkS2WD3bl3ae9cyxMo4G4arq7wCAiDxJVT/g3L8LZyJzFwdCuKr6PgCXOLcudcoqgBd1N9LuzNGuX1GDbV22oKqKxi68qF2vXLU/rP22DlvQltRZdEDb3DrzfyFNIMfO1Nct25+mviXxiIw94qqkEjxCi8bQnjNl7ZHvKqgQqUeMbJ5Gyrata5XtnMR7SFIGHfh5ETkPwL8H8Guqel+14np/TY+dJNP5bkIij3SZuolSAO39yqKxasNTpbY+Q9R/prC8frj9M3nZ05eFvrZX33ltm+0JPSmHHXYMcbSvFiw1xIjbKuHKPPM2VEaWUX1W1/rINhTrr3e9HaeZ0wpHIYc7QVW/CcB3Y/jb0ntE5NdF5FmVuuv90V5vwrTn7bVV7LZqtaKU7HHrS2sn8jPKAdrFwFR0RtKeuqkgSANMyvU0mcoxYDq2qla3uOq1irlKvO3YMrXuESrrv0dOLHLJxtuG79Yem2NteS/t0d73ztlc8uaQFRZzKVwFjuIHG1T1IyLyLwHcCuBVAJ4sIgLgR1T1Laze+ircbKEwNdqqtqoithOwXQDWFvOVHXtKNyqTqdd2QVtVxRZtdXGxP4RNzR8T955sjOTbknHzClVvFUypVhQiu+bV9567JS3vmmeD+e/NF8/3yrP05kJ7vfWFbeC7hA0yvFcFInKZiNw+fg3ANc7980TkjeP9m0Xk8c29rxOR3xeR20TkA+Pbv1g7Xycir8TwCdlnAvj7qvo/jsevjHw8OxSu3bE90orqt7B2bDttGWur53o1VcDKVTYJVja6l6HJ21KSPCYQEeDYMWBjYyi7uQnd3DxNuiICnfq6dey06t2R4+3yzYkgmIL0FGq28tk4Rsowe57efPPUp+f/KhsnEw8zk+w2zPAuBRHZAPBqAM/C8NbTW0TkBlX9UFPsKgD3quoTROQKAC8H8A9E5DiA/wDge1T1j0XkkQBOBs39WwyflP0RVf2b6aKq/sWoeinWl3CjtAGbPNnEsuEbu8fayFRoFrJ56pyVq9ZnNtpyXsidEPG21MGEhkyxsQE558z0U2DH36oFALa2oHMl8Vh0422YTAUzm94ci1IA7XH2PNm9OZQnmwteKmLVNjIX5jH5VAAnVPUOABjfRno5hrebTrgcwI+Nx9cD+HdjGuDZAN6vqn8MAKr6maghVf17wb3XRXXXl3A94vImF8u7MZtROW8xRcfWJ6tiMvXFUgBefWufhbSsz5nqzd6NMJHtxsZAtsaGbBwDjo/TcVK728bpjMo97dqWnr4+C3qVbETG1Qgj27gneM8he+5TGTZ3rP9RpOWS787i3dCynQtF5Nbm/Nrx06UTvK8AeJqxcbqMqp4SkfsAPBLAfw9AReQmAH8LwBtU9aetAyLygcDbLwL4UwA/ORG3hxLhisi/AfDT45fMYPzqxH+qqqF8PhTIwjhGhF79bDFMZew1WzYiec9WNayv9MHatn7Z4yqit4A1KQQ5fnxII2wcG9rZ0jN1N46dKQ+MKYam3FzkasfYI59qfXbfzpUoRRSldCpzi20AbI55m+dU10spRG3vGoLiH83uVlXvraRz4DiAvwvgGzC81/+dIvIeVX2nKfftiY2vxfBWsSdHhSp4jqr+yHSiqveKyHMBHH7CjSaPhSXX9rp33KYX7MLxlGyP2vbaiYiwJ8/G/PLSJWWy3/4WsDNtbQHYGMj23HMGshUZy28BGwLI8TOkDEA3t4BjmxBsQDc3xxTEmNPdmoF0qyF8xQ6bL2xjZvOlve/5ZjdDFq2wut46YG3Y616d3g05wjwcXvkKgKnMnWPe9uEAPoNBDf9nVb0bAETkRgxfN2AJ92Oq4YT5UxH5O5GT1XcpbIxv9MXo0IMBnBeUPxzIlID3aut6517Z1nbbbqRqvHDOLkjvntcvey3rh1Vdnu89is9D+y6DMY0wnI9ke0xOq105vgE555whpSAypB+OHwfOOT4o3PFl0xCzomI7S720Y2ZVM1OelWjJizyYSm3L2TL2Za9n2Kvx3yq8ctwC4GIRuWj8BOsVGL4WoMUNGL42AACeD+BdI4HeBOBJIvKQkYj/Hrbnfie8W0S+X0S2fXxXRM4VkWeKyHUYPjlLUVW4v4ZBZv/KeP59OPPdtYcb7cS3ZOYpFG9nr6iXtmymNr22onDOs8988+p5bbflPNXVo5gZjsn2nK1uDeHjhgwEKzLkblWBza3h/2PHIKpQOQboJmTj2CCAtrZO29z1uxYyoqugHS+mEjMFzGzaa96x13Z73V6L+rAK5kotKKophdjMkJN9MQby3ADwWlW9TUReBuBWVb0BwzsLXiciJwDcg4GUp4j9FRhIWwHcOH2M1+AyAP8IwOtF5CIAnwXwoLG9/xfAz6rqeyM/S4Srqi8XkT8G8K3jpX+lqjdV6h5K9BATUxwZUWeTmZFlj4pgdTLlytqxJOIt8Mwl+1awra2h3jnHgePHB6I9ZgKr6e1hW1tN+ycBbACbm2fKAMCpU2dcmotst3XACdNZRJGleaopgvY8830qy8iclfeuV0jTppz2ADO9SwGqeiOG715pr720Ob4fwAtI3f+A4a1hkf37Afw8ho/2ngPgQgB/M/1tq4Kedyl8GMApVX3HKL2/VFU/31F/f2FDOiBeRFadZgsuImQvrPRIkakU1h8PkSKyZSIFHaUZovar2DJtTamD1q+mDREZ3g62sTH+AW0rUGoded3KRsjOe1UpI2SbPmptexFYppSz+cM2hWiuRP6f9ilutoy94fE9haqeBJB/a6FB9V0K/weG3xO7AMBXYXh7xS/C+bKZQwNGMlHIXMmNVlVgNSfW1gknNwlfvbpRyqPiW1HNulW39IzK3dIhfbCxMShcmdIHm8DGBvScc3D/4y/Anz3vOM757DE8/oYv4Phd9wGq0NbXjQ3IlkJtWgGbK/nYHUozEmzvR2A5VXvfHrf3vWPrI2vbO4/mbaTM2zILulFVuC/C8MbimwFg/Bzxl+2ZV3Mg2tGrijWzXQnlo9CvQpKR8owWp6eMmBKrqOQiTpOtbo1v/9rY7t7mFuT4+OmyLcW5996PR/3XL8Wxk1s4dv/44Z4to1YndbyxMX4ibauuaD30kkV1fLxybOPLxn63Pra2mCqOoixP6VpbM2GulMJRQJVwv6iqD0x/JR7/knf4h6kSSq1CNquQmKeKvYXQo1izxcmU0V6ok/G9tKcVrv2E2ebmqHiPARjfGra1hY17voAL/uCvBp+2toY/nh07Bjl2bPgAxKmtMx+EMGS78h/NIgU53a+kD1iqwLtn26pEV57P1n5EjLataMP3bDD/5yRdxSwf7d1PiMgTzUeGISLfoqq/m9Wtvi3sP4nIjwB48Pg1ZG8G8Nvdnh4ULMG2r/Z6tYwlQhvae9crPrCQM1sE3mJk6iqr2x53hd1kKqkOJDl9wMHcO/06tTm8JpwmV22K60DeY1u7+l6FqX3AH3f23KZrkUK1c8EjPu8ZVVIT2flks5r28PrEiDkqsxto4XW48CYR+Rcy4MEi8m8B/GSlYpVwr8HwS7sfAPCPMfwl8PB/6AHwc2/RxGLX20UT5cHYYrNgoaZH/Ozc+twutqhtr749Zj5WoSPJbo0KVRXT+3J1aws4eRL44gPD6+Sp8b23Y/pha2sor1vbPuYrIj5xr+RfEhXY1Ew1pWPv2bp287THXnmPOCP/vTSB7Q+r15aNIsEZoyTR/HXI8DQMH6D4rxjeSvYXAJ5RqVh9W9iWiPwWgN9S1U+v6uW+I1ocmWpsy0R50qmcJXRG3pEPlfCO2fMWGVsw0aKdIWQ8nVY4JmfezrW5Cd3agoy5XT1+HKe/MWxzE3qy+XImT7lubfHc7SoEHD2X6dwjQxblbPNHt88H79n0wLZZUcIZ6Xr2J7SbtpfqOl2mvyu+vzPZ2T+cxPBrvQ/G8D7cj6rWJmGocEfJ/GMicjeA2wHcLiKfFpGXRvUODbyQrqr8ogXCFHAUwkUT30sxWJteGMtCVusXu+71a4ZQsc3h6ubmoHLH99ieTi9MKYORiAEMed7NTWBrfJ1ORTTvWsBA6OGXkntg86A6N7znwMaLpSWYYvUio9bv1oa93h5HdavpABvBsXk9o8I9gimFWzAQ7jdg+Ln0F4rImysVs5n6Qxik8jeo6gWqegEGOf0MEfmhXTi898gmhA3fbD0vlGfKdjpm5BctYrZ4rVJuffAWsLeYe7BbolVHfW4psDWQ7PA9CM1vnY0pAz15EnhgfLXjN/2R7NSpgbSztuk98qwmVPrdG06zDbC9z8iRteNtEhZsHmQb81S3tWH7kfVpRVTSCYcwpXCVqr5UVU+q6idV9XLs/Bixiyyl8D0AnjV9qQMAqOodIvK/YfgoW/jt5ocGNhSMwmYvLGchF7CTAG07zEY1zKzkDldZBCxk9uyX7AU/fdPem3KyG+M7Dbb1Y1TGU853UsZj3pb+kYx98GGVfmyzS55PlaS9dAJTmd7Ye+qWkbQtG80RL93Bytlre4Ej9i4FAHfZ71MA8J8qFTPCPacl2wmq+unxo22HF27OqbBT2wnPiLK14U3SzHalPEMPMVYWK8tf9vjWfD1jS4zthyD0GCBbW6fztgCgYx73NPGOhKynTplPpx2DHNvaYZ+q28pG1kYM9plbYqwQISPT6liy8N367dXLjm1dKxqydlxy3lltFRxCBZvhd4DT35n/IAAXYUi5fk1WMSPcB1a8d/hQVSoeyUaLxVOKGexij9rKlItnu60XhZ/T/TkUDMmj6pZCNnD6Y7yqOrwrAdj2lYzY2DiT053es3vGyE51W/kbhSWJaLwjZWfrZ9FFpEy9+jZC8pSx9S275vnnqe1M6baprd2mnaiPe2N2r6CqT2rPReQpAP5JpW5GuF8vIp9zrk/MfriRqQ5GwnZiepPeKubMBzZZWfqCKR3vOCP9ntBw1UU1KVwnlTB8aAHbPxDREvCpU00aoaBeJ9vhfee5R6otIjIvymHlsjDeu8fssygkeo5sDkV1WXt2fu9FakGPpMLdBlX9IxGxvy7hIiRcVd2I7h96MIXD1Ex1UlcVcKQKmIph51WVES22im/V+q5NQ7rer0BM97a2X9/xu2XqpA8sItKN+peF7VH9DOxZ9pBVlLpo7XrztUqQnnpuUYkI5yLKI0a4IvLDzekxDF9W/heVuuv9m2b2uJeEsgWbLeosl8cUqlc22zhasE3FlmHX2nRHdewm8mt/Kt1g2xfbEKKcCLb0s+i977/1nhcjlsom3da3YXqPD1kayEttMN+t3ahu5quttxcKF4DM9JN0+4gvbY5PYcjp/kal4voSrqcSonDKO2cTjNnJJjAjYLsIvHpeeFfJUbb+eoTq+cN8jlTfjo/tGuJtvmthm3n7o5DOcdqWh95c5zb7SURRSft4UVBGqFW/KkTZlmdk2WOj2uZZAFX98VXrri/hVkIiVs6rVwkVPUK0bUX51Kr6yupki9/W9zYlb8NaNb8LnFHARvnu+vsQMmSpAdZ/73yyx+55bUeIIgg2fyIC9OZW1YdVUgxz4YjwuIj8NgJvVfU7MhvrS7hRPjYjMruosnDfQ1vPU4keKWYKl5F1tviy9EHrmz32xqCHeLvDfp6S6LZXJYmKOrX2qht1W9Yj2IjUe+yxeozUozkdKfht5Xw3u6A4Sn80+5ndGlhfwm1RDZMrubNqqqLqS0TAnq9eGe88SzUwpcTajkLsHniE2ZJr9AGKkv0g3ZPVq6ZPbLloLniRUER8Xl1bJ3ouEXmvklaIRMNcODqE+1FV/dhuDFS/LWw90S6GSMGxHBhbBNUydjGyyZ8pnuielwpoiTZTrt4Cbl9zYPpjm/cq2yDPh/WJqVjvnEU53thNbXrn1p+e52rHvX3ZTdLOuSo52j7Z471KK2jhdTjwW9OBiJT+SGaxvgqXhd1R+FiZmFHYzULPChnbcxbOVfz01Fpri517aYTKpnIYwHKRbCyyzSXqexbSs43VRhwsOmLzyvpoffDSFV4dz3evjyya643qAghwlN6l0D7Yr1zFwPoSrkUW2kWhdAXtZGS22vajNIW1Z+ux/BtbCBGxM5tssUZpioOG99w8kvPgPRdrs7LhVcbE+sRINFLAUVTU3ov6nPnByP50ed9sFxRHKYer5LiMA0spiMiGiLxXRN42nl8kIjeLyAkReaOInDteP288PzHef3yxAX+itdejXdqWiUI6W2e3qCqn6FoU+mehog1XvXtZ3YNEJV1SGZ9KPyobkDcXM+XK/PY2cy9t0fqVvap93Sto4XU48PUi8jkR+TyArxuPPycinyefyN2Bg8zh/gCGn16f8HIAr1TVJwC4F8BV4/WrANw7Xn/lWC6HzWnZCTjds+Vtmem8khaokk01ddC27eXlrIqzbbQvlkaoKjJvMUdl9xrZ5jf5YuvYe9bfaPNhfkx2vHmTHUf+WZvtJuip3N2QZvRsPWKfk6CPCOGq6oaqPkxVv1RVj4/H0/nDKjYOhHBF5LEA/mcAvzyeC4BnArh+LHIdgOeNx5eP5xjvXypSWNEREXpho7cAPWKrhmisbQ92E2CEXFm0DEwpZWqp9a+i6iMCXNX3rD4jzR4F1/PMIzXKnl3WttdG5DPbQL3oJOqTN5+965X+r4gj+H24K+OgFO7PAvjnAKZ0+SMBfFZVx6+Rwp0AHjMePwbAxwFgvH/fWL4GpnTba0yleMTBFHJr09pv61ZCS49YV5nk1ToZWXnlqmRSIemqr2xTau95dirElaUVGMFN51OZSmogU4lZH5hSZ+1m89hLL3jn1v+5cEQU7hzYd8IVkW8HcJeqvmdmu1eLyK0icutJfLG9wSdHFK4x5RaFVZaMrdJoy7T3IjKy9z2yZ4oyWxReW6xcVe146ilSiJ4KrfSnR+2zDSUq46lCr/xuNiHvnG3mTIm29T1fok0+m3utP/Z4LiggW/lrXXAQ71J4BoDvEJHnYviKx4cB+DkA54vI8VHFPhbAJ8byn8DwC5l3ishxAA8H8BlrVFWvBXAtADxMLtg5M7KJZa+3E5Wpn2nC2sUXLRrPvlXWUVusflu2JXLvf68NG5Z6920fIpUWgRGDbdsjStYmU7wZUXpkz4jIPqfK/PAQRUYeUfZEG9VNlPndEnF1PewWa6RgM+y7wlXVl6jqY1X18QCuAPAuVf1uAO8G8Pyx2JUA3joe3zCeY7z/LtUVnnhEdjudjOva4x5lyVRv6w9LS3iKywv1bJmor9kiikLrSv0q2HPw1OBUPhoLr3xrs73XOx9sOU/JV8aFhfHVqMvWszaZv9Zvz6bnY1ZvRSw53IPBvwDwwyJyAkOO9jXj9dcAeOR4/YcBXFOyZidwO3Gn++3/Lbxw3SvDFoK3eGxbkaJifWF1s357i9raYouIjZk3tqzfFTBla8c5Wuysnr1fJVfWXg/ZVMuyPmZzIxoP++w94ZCpXGZvTmjhVYCIXCYit49vId3BE9lbTEXkK0TkCyLyz3bTnQgH+sEHVf1dAL87Ht8B4KlOmfsBvGDlRlio6Tu0faJ5obidqExptvXb67ZeRX1mYbHnP/MrWoCer0z9eWFzNEbeeLGxsX2yfWckUYE3lp69Svu2v72qL5tr2fOytuxxRVzYumy+9KyjHnQQagQR2QDwagDPwvBH91tE5AZV/VBT7PRbTEXkCgxvMf0Hzf1XAPiPu/eG4zAp3HnR7tLtzhwtChu6e0Rj71mFmJGIDW2ZIozCx8j/bENgfbN9t37YsfRUUNQPb8wqiBRmpHijUDvz0Y5P+2Lk7o1ZT/+8cba+RfWi55CNoQUTD9V11AEB5kopPBXACVW9Q1UfAPAGDG8pbUHfYioizwPwUQC3zdAtivUl3AkRadoymdKppgosMdl7bFFXwvFooldDQG/xeOkCz2bU57Yvka/emHj32nYYYbD+ec8486uCKLpgiBRze53NK0vCnp1oM4jmVA857xFmItzTbx8d0b61dEeZ9i2mIvIlGFKaP77bvmQ4+75LYQILo7yQzqo5FhZ7dqMw2aomzyemnCqK2CMFptQ8IqmE0wyZcmL3KqFy1k7m28wqzY1kmD92fG1dRpzes2zLRfPM89Gbcwxzj9cO+6VSF4rIrc35teM7k+bAj2H4lOsXKp+p2g3OHsIFdk7mSojLHkAlLPYWD1tk7b2oXc+vSL2wPrL6LESNyLsKRo7MN0+xe2PZ2yaDVZTRs/f89p5DZXNu27O2K1FC9Cwr8zRrs7JOdoOa2btV9ZLg/vT20QntW0ttGfsW06cBeL6I/DSA8wFsicj9qvrvah2oY30J11OqUf4tIxhr25b1wnMWdnpEyJQn65dHzhkZVRZOtqiYj6xelfBsWsParYa63vhUFKN3zbPVlqs8p7YPXoTE2vPsWV+9Pnr9zFIK1j7zc8f4+Ca7UE8ZZLgFwMUichEGYr0CwHeZMtNbTH8f299i+k1TARH5MQBf2AuyBdaZcC3sxIwUhLdQIxJluTUWJnvXqoqxh/CYarL2PDsZaXn3vDGqqFBmbxV4ZBiRbbQJe/ervtn2I3KshOzR+HrHLIpqy2WRRetnK1jmVrozmFPVUyLyYgA3AdgA8FpVvU1EXgbgVlW9AcNbTF83vsX0HgykvK9YX8Jl6qf9P1N8UdjGVKWtHy2SiOwiZOkArzxbjJ56q6j1SPFlobZnNyLI1o7XxqrkzOpGUUc1vI7s2mOrQr02vI0wmgdzpTu8tk8f7zS5Cub66K6q3gjgRnPtpc1x+hZTVf2xebzxsb6E66kuOxG9Cd6jXjwyysgiWqi96oH5m7VRaYsp8NZG5Es0DhHZVqKKyH4V0bPuTT9ESpGNd6UND1HUZO1V0zw2CsoU9G7H3rozs2A+zFhfwo3IwaKaOmhtZ/m3jNi8SZ7l2dq61oYtG20i2cYShdOsTxVVysraTbFaj5WroHeTzWy1zzBSstO1SM1Gqr8nncR8be+zaCdS0XOM2WmbmE0pHwWsL+FWQuG2bGUSeZPey21VyKw9tgSakVmkPr32vA2lsrHY892qMnveQ9Ksfs/CzwhlN/DGuWfTztIA0X32bD1EqWi8R+oAAB/qSURBVIW2TLQRZG30YiHcswR20a1KLp7SZHW9RemFnlGqwFsgUSgY9aNXkU/HEVlG7TGfqptea9OqtEy12s2uRwW259E1FoJnm+RuUk3Vfnj3PF/b8Zw5fbDDDSwphbMHjNwq6ooRZ2vXa6+i6qKFYv2wNjM/osVXTTVk9qytXiKtRgvMB7bBRW0ym5Ev1X4xUs5UYjUCiJR01kY1HWHn+IyQrbOHcdebcHsIwpZh6QOmSrJF39a36qGijK2dan+qaZLpOAtjq/a9POV0nbVdsdtbxrbbu7l4z5ylaTyfKuPJlHF7PyoTbaKVOdQThVTEQA8US0phLcHCQ0/BtOFtFA6yyZeF1UwZW6ySM7N9aK9HaiYj/woitbtXyJRw1C9Lnp4t79lVFSqbN71Rh4fKJha1USln7e7YnH0TvVhSCusCtpCyhVVRW1kagS1Mrx222L20QaU95nsvmVbC+sjH3drvCeV7wnuPPCI/LKL50PrdlvVSRpF/9ppXLvOrSrK2LIvO9goL4a4BvAWVqVavrr3HrldCRa++DdU8W4x4vLr2fuZ3RB5eKM3sRO1XUzSeD1W7Xv22LLsflemB96wrKZksTK8o8+map0ijzamyoYYkzW/1YFG464KKQqoooypxZKFhlAushMZZDm9HyJcoQW+ReuWYX54dzwbb/JjvrD22QTFEEUjUTgURkbL+emPtPbOpvBfFeP97pJ2JAfZ8Jz+zcZ5T8S6EuybIyMuSwQQ20Vi9akhqw8qIdHrD87YN1mbUh8hWdUNi556tqsLy7Hobldd2dI9tLuzZs3wmU6g9aZXejZvds8deW5kv1kaUKpmDKBVr9au8GdabcC0qC87ea2EXVBZyV/yxi6xHxTIlsgqZeTa9+syPqCzzo2eD6q3blrOqLSsftV9pN0sRedfYWLX+e35kczLym6n/yLfo2goQLCmF9UBPPi1TBVH4591j6YKo7UqY6PXLI+oKIhXf/m/LRsrf8ysi7GyTqqYNrO1V60f32PO11yI7lc0qIt62TiWNxDaLaCPJSNqWmQNz2zvEWF/CzVICwO5JqoeoWfogCxF7fPTUEFtcreqLwMJtr27bt0yZ2eO5UbXdEwl4aZH2OlOMNmLw5gCbmxWRED1H5ofXn8hn2/aMJLko3HVGNa3Aymeh3SrqwitnyToKFXvymMyfysJeJWXAymTk3dpfJQKx/rVlo0038itSsMwmi47sMYvIPF8ixZzN6ep8Zvfnho6vswTrS7jRIo4mcYbKYsxIsmKb+R3l11pVZFV1RNjW/7ZcFqZn45aRYGSjQrRReJ7V70W0AWTKv50TTAlnfrH67BlVNrFsQ6kKhV1g+aPZOoMthgiZUqi0Y+2x8l6OkxEpa5fl+GyZzJdIoUabSaYiPaJgC7tng6yQu23L61ekVttUTKTebXtZuM6Imo0HI1pvA14FPamWXWIh3HVAhUS9a5WFvBt/vInsEaRHiIwII0R9ihZ55r+1wfrFfGLIiIvZWXUDtITH+mev2WPPLiPrXv8Y+UUpjMwWe/ZTOUby7f05UgGK+lxeA6wv4WbI8lfeorIEswoRVzcCT2llE5PZjhZXpEQ9G2zRRxFAlWB77rH7mQr2Qvi2X57f0aYS+WQVsXe9Pbc2mGKtEK+HjJhZlBW1PQOWP5qtA5hK27FLOwtiup8dr9o+W7wVgonSCpk6s7baeh6Z2vv2PFuErd0srRDZiO57ZbNr1fa8OcE2PrZBMxvWD6Y+o02DCQHbVhQZRP3I2p8LC+GuASpEY6/1qNZe9WVzf5kSWiVU99RwtX62QdhyWUja1q9g1YjB+stsZiTZ2vBUZaUt756XKmJ+eqTLxECkQpkPzA9L8JGqnplsBYvCXW9UQ2hbB1iNmL3yXujW/h8dR9eYyrL32IJlbXuKyyo9GzFkiNpkZXrDWjaGkZ2KX16dSvjN5h2bV2yusjLeeaUvbD6ysnOSruryBeRrjSzUZnWicw8sLIwUznTs2WJtRNejsD8jcY/squMQEQgrl9mJxm0VdczqeAqvogq9/217mT+t3ci/1s9s3lQ3a1Yv6u9cPHn28O1ZSLjAzknbhnQtvAnNFgMLp7NQPiKxqtrKVJDnvwcWekZtRshsMGLzymTk0wNWns2DqF5Eeiyq8TZC77l79ioqtJISyXz0+jVH2sdzYSHcswC9Cq6aTmDK1jtvy0dk6ZXN1Fm00Lz2dhO+RznA6rUs5J9zsUcbpm2fkRzbWNmGVtnwrO2KjUwsZH1l1ytzcA4ogCWlcBaA5a2yhR2pHzvJGamxcI0t+vbca5Mt+KiPLXajYm25VUN1S3K9PmTINhTbZuuTV6eyUUbPtbVjj737no9RHz1SzqI0b8NuN50dY7TThZVw9vDtmhNuL5l69XvKe2rH3q/UaxdIRT1NdaYyljDYedQP6wcjE1aX+Wfte+WYf9WQ2kOm2LLyHmna58RsRPOPbYrtfTYm1fnGNtzI5yg9wuquiLMppXBsvxsUkceJyLtF5EMicpuI/MB4/QIRebuIfGT8/xHjdRGRV4nICRF5v4g8paOxXBV4KothKhO9MuWRLU7Pz0hFTX1s+5rZz2B9jOzaOpGay/xo+1Ahkd5F7ylWq87tsd3IWr8ipWjbae15xx4ZsufAfPTsen2zc6Z9Rc/Q83cGyJamr3XBvhMugFMA/qmqPhHA0wG8SESeCOAaAO9U1YsBvHM8B4DnALh4fF0N4BdKrXgTxZvEE+yCmq6xMt7iYwvNsxmpjOl6xR6r5/lj+xLdq9S3G02k1hgxRZtFRFyZf5HP0X3bXnYtUpmer4y02T07Nowobf8y5RyNQfT8ZyRaADj9bWHZa02w74Srqp9U1T8ajz8P4MMAHgPgcgDXjcWuA/C88fhyAL+qA/4AwPki8ui0ITZJ23tt2apisgolK8NCwfa63QiysNpry5JxtNCyvlbUuyUJb6Fb8mQqOIoQotCWKbGoXxEp2XYyZW99rqhx9qwrYM8yK8PsV+Z8JEBmgAAQ1fS1LjgIhXsaIvJ4AE8GcDOAR6nqJ8dbnwLwqPH4MQA+3lS7c7xmbV0tIreKyK0n8cXti4GFQowAonuRQmEKjE3UygLy7lfCTVt+FXjhpr3X+pyReFvXex7es/HGkW1ilf54ZXuIxNucWQQVqfFIbbMNr73PfLdz0ZJ+Zdz2m+C2Cq8CROQyEbl9TD9e49w/T0TeON6/eeQfiMizROQ9IvKB8f9nztEtDwdGuCLyJQB+A8APqurn2nuq2h1IqOq1qnqJql5yDs7bOckzFctUsK3HFlWF/LY7vJNQqnWtv56vkZ1qGxUl3NqLSMUjy0gVM3Xo9SEipEp/omca9Z8Rd7apt2mIaDOO0grMP5YSYNfa/lpkxD8T5lC4IrIB4NUYUpBPBPDCMVXZ4ioA96rqEwC8EsDLx+t3A/j7qvokAFcCeN1MXduBAyFcETkHA9n+mqq+Zbz8l1OqYPz/rvH6JwA8rqn+2PFaDZ4q80gpWtTe5Lf3bJseSbf3PYXokVG1j5Gasv2xx6yst9DbvjF1xdS9JeAKITBbTG3b+0whMqLyfLbwyK8lUWa3bZ+RZjSunm+ZQmWbiW3P2+wi0p+LeLX4yvFUACdU9Q5VfQDAGzCkI1u0acvrAVwqIqKq71XVvxiv3wbgwSJy3sp9CnAQ71IQAK8B8GFVfUVz6wYMuwvG/9/aXP/e8d0KTwdwX5N64GCqx5JdRLqZuovKteQeKYv2/ioky4ik517Wv9bndvFlpDfZqfQhezFlHI29d8z8YyRsfbV9ilRxJZLwNt3WjrdxZIrTm2vRZu7N0T1Qsj7ydyiM71K4cEoZjq+rjaFK6vF0GVU9BeA+AI80Zf5XAH+kql+cq4ctDuJ9uM8A8D0APiAi7xuv/QiAnwLwJhG5CsCfA/jO8d6NAJ4L4ASAvwbwfaVWIpVgF3Bbx5uMURmmPiKbUXl7LyIB60O0eFn9jLwi2AUaqUHvvN2QPB89xcWiAY8sPHXutR35H41JRPgsCmjtso2APS/bd9aOR/hemeomyebgXJxcm3d3q+olM7XoQkS+BkOa4dl71ca+E66q/n8AGNNc6pRXAC9asbHhf494IhXiLWymULN27aK29iM/Pfts4UT+VPzOyMHaYLZZfbawK9dZ+944MvXNNtnWfja2lvhbm2wD93yN/LCRUWTPO2Z98/rAwDbAvYBirp/YqaQepzJ3ishxAA8H8BkAEJHHAvhNAN+rqn86i0cODvRdCnsOj+jYRLXhW0R+mXrxfLDtt4vWs+uRvA03WzttW6yfLAS17dtwtU0j2PLRuTcGXh0WRVQ2AHbNji8bJ+9/ZoNthN64s428GvFUSS5rxz47+2wn9EY3c8Kbd948jHELgItF5CIRORfAFRjSkS3atOXzAbxLVVVEzgfwOwCuUdXfm6lXLtabcKOwabrPHigj6QnRhGYT3C5YtlAjdcF8sccVtRstVkawvcrHI++oPY8E2T2vL9a/6tgwpe2V8TZqby5lm0tEfr3XvXNvo6tGKNmamRNaeGUmhpzsiwHchOG9/W9S1dtE5GUi8h1jsdcAeKSInADwwzjz4aoXA3gCgJeKyPvG15fN1LttWN/vUmDKNJtY3uT1lHL7v2c3Cg+9he+FrJ5f0QZi++XZzAjZC6uriDYmhkjtRSF4xceoLxVVzXyoKFFv/u1GtfZc9+Zrq9bteFbmRWu32o8iZGuenIKq3ojhbz7ttZc2x/cDeIFT7ycA/MQsTiRYX8KNJgUjtGzitmXb+4x4q/WYv4zII/+YzShEt7YYua2qbFh9RmK2ni1fHW+vjEeilbGPnl2EaFyZnZYUo029p/1VjtkGs6ovHhTlDzasA9aXcCfYHZ2Fhe3/9npEONa+VZUtWpVhEfnIVGL1GuuP148KKXv12uNImba2IsVp26xuUl4brB1P2dnIxKLaR/bsq363/1fEQWSnB5UNaEaFK1ivj+5mWH/CncAWuCWZjMQi1VWZlFkIGpG1VydTuQxM5Wf+MtKIVDZTlFm6wVNTlX6yNiKi9PqTRTiev9H9yOfIvt3Q5gzpWbqrPffExJxYCHeNEBFG5UF7Ciiq75FvWy+yHZFYxbfWr4yMesg18qdHhVfGhPniXbdjXek7q+tteJ7arfgfPdcoaqpu1t5mEG1sUWqiEtHsBcm2WAh3DeGRkne/qoi8a5lS8HyJwsWq6vL8YIuM+dWG3tWykdpnqIbVEVFbf2yd3ZJKlFaobFSVdj3iZQQdzSV2r9p/zxd2Lbq+KhRLDnetYMOwKKfXnrMJG5FQFvZn9noULSO3jCwYwXj+ZH5YWMKL+hOpw2hRM0WYjZ2nritE2eMPK9OWY2kHL3T3oirvPENbvqp091rVtk3N9C6Fo4D1Jlw7YS28kNLWXbWditqb6ljiyNR4psbYgoz665WvhKtRvxgigow2w6x9RmZsHlTUeLQxsXREWzZS0Mx/du7Z9cjU89cr64mRKKpq/ZotDaAz2jr8WG/CjdQjC8Ojuoz8IqVUsdleq6YRVlU9LF3gtVdV+dY3217PZsZI09b1nglT96uqtWysMjLNNha7ebDn39rKSNkj1bYcS7swW6yMCGb5LgVFPB/WDOtLuEx9eGrHe+CVa1k6wZZhi6Etm4V6Pb5nSpYp2F5UUggeKmTMiCYa80hxVnyqtGHvR1GIVybb1Gx5prazCIjNHzaebMP1bMyFsyejsMaEO6Gy6DJVZe14iCZ9VD8j2OoCaH3wNoIoPJ7qZOTNEIWcHslUUhK2/YzcWDqlsoExeORbSWd4KteG7p5y9dIQlUhm1T7Z61EKIhr/XWJ5H+46wFssWYiX2WD3vUUTEVymzKxtr0xEkKsok4xIIlRTBr2LNiPobDPyiMyOcRSyM+Js/V9VBXt9m9rwIgZvjLMIxvOvEo155/Z4TpJcCHcNwJRBVNa7noVsbf1sMnshoqd2Mh+ZAooI2COmyibAfGjrMSUWtVVZZFFIbdvyrlfC9oh02bl3z44za8/rX6Rw27JsQ4zUt/cMsnm96qa9ClSBzbMnp7C+hNuissBZaFdRZV7d6foERrYRAbT1IvXs2bc+VAgvI+CKMq9sHswWI3SbLvDqrKr4WXuVOt4zsKQVjSmbV948jBQ4I2bbdtR+lWy31YuLlbEo3DWCl07wFqe3qJkSqCjiyuRubVXCvIx8bL02PI3899RplTxsn5l/VYLNylTD47Z8FlZHNiqK3fPZU8keGbLoo60b2fbKs3FhStrrr9fPno2rBwvhrhE8omBE6qFSxhIVq2eJ37bDQlnWnleXESxbaO29SMEzf1m47fnCQubdPI+e58POra8REWUbj7Vn50NWl41R5HtrI9qg2RyobDh7BQWwtRDu+iCa+NVFblEhCq9MRmhVXxhxtsceqTE71lalzUqYz+5Xyq1Cxq2tqD/eBhzZqY5TdTOokmjrA7OZbQR2856DQGdVpAroksM9+ohC+2pY7tWZ6jG1ypQK86klYo/Mehan9bVtL+o3U12VfrFQOfM5G3uPjHuI1wvj2+uZf+11zw6bV5lqjTZez3bP869GRD3wlPGqQsW1j+WPZmsBRlqR6vCIitWxIVyFzJhiytRilFer9MVer4S+tq7nM/Mr2zh6+tS2VYkQojDeIwuWErC+Rc+5SmjRc2OpF28+VQkv2vw9m+w8qjsHZlXMhxvrS7gtKg/UU1TRIszSANHkbstki6GtV5n0lfAxWmCVdAAbK9tmRXH1plGq/vQqwUr53vQGi4SqYx2lPrJUR8+mEpHrKtFJLxbCXSNUFkiWx5vKtOcZvIXPcmzTPW/ReORZ8aWSUonqMp9YWc/nav2IyFZdjBUF6inJaCPyykURjddOC89GRpyVSMprN5rfPQp97pTC8uU1a4gsXOtVwMxO21Y2kb30ha3DUgBVX6PUgfXXtsn89OpGaYlM2bW22UKuKMwqqdixjjbGiPyq6QHWB3Yv8yGKXrzNIHqeFTt20992b2fVbiiA5esZ1wjZIuxZvIxcexcY848tcKZ0Ix97F2qFjNm9jDBtOxWitu1kmyNTgra+Rx6e76tEBl57nn/V9Ijnu9087HXV/rbs84hSClU13INF4a4BsjDSHu9GObGFWp1IzFdGnp7tLO9WaT9SdKwt1n5EntY2U/kVAvYQEU71GVfmRlvWhtoVgq6kFLxrma/Wr7ZM5Tl5bewZdHmXwtqgd9Gw41Y52AldmbBMbXllvXJMsWZKJtp0VtkYqqiodnteSUl45xmqm0+WZmARQDRvrA+VvmQRkfXLgzd/7b1qRONt+nPOFwV0eR/uGmCVhcoUSkQAVQKO8pft9WhxrhLWVhRTD3rJwCOvqJxnbzeq0SNRFopbpd9DsNae1zbztYWdcyw9EqVNbL0saqhEFT3zrxfLJ83WEF7oH5GyXYg9aQWvjSpReOeWDDICWlUFVuqx/vZEED0+Zio+ayeyET1jm8tsr2d+ZErXa9+z35u6YRtHVN/zsSIwZlW5C+GuByLC6128UUhVSStkypWdR9ejPrRpkMh/bwF7C9mzz671pBR6+tT62AMvirDzovJMIwVceRbsekVhR8QczRnrOysfPfdVN8sKVJd3KawFdhM2W/SmIzIbkYqJ/J4WRTXMzkLjqoJjYXK2yL1z1idWtpqyaX3PVFqvOvfqVjahjNTbZ8nSQJEi9dInHjJVWlW9OyIt31w3FoW7BshCfO+8vd7aYWogqpchUqBRGxarqOveHKNni/mbkZe3WbC6vWrWKr6euhW/vbYiP7I2WL+9+lmqwJsvlfGMNpAIs5GkQjc3Z7J1+HHsoB3YM3iTkym3ijJh16zqstdtGU+NVEJIa38VRHlJtogz9dTa6G3bIxnvudm6lbay1EUP7ObI7GXPpvI82bNoy9pnEtWx7bKxZHN1L9MJAE5/PWP2KkBELhOR20XkhIhc49w/T0TeON6/WUQe39x7yXj9dhH5trm6Z3FkCDcbzB2oEgVbJJaM2/+zCR61VUk9RERXbTvzczd9qC661odI9XkkUG0nUnvWdoW0bXlPNVqVXvErexat7wyZeGiPLRmz+63dbMz2CrqVvxKIyAaAVwN4DoAnAnihiDzRFLsKwL2q+gQArwTw8rHuEwFcAeBrAFwG4OdHe7PjSBBucTBtJV8leQsqmmge8dl8lm03AiN3z2fWl+j+KkRaye1N5Sr22zLeWGXRR5Y/9I49X1mUYW1YvyJyYna8OeT507bFVDxLubC6re9eX3o3rkqdmYhYAeiWpq8CngrghKreoaoPAHgDgMtNmcsBXDceXw/gUhGR8fobVPWLqvpRACdGe7PjSBAuaoPZDxZGVRRMZLO3DlNDbOG2ZaZrURjI/Mz6WA3LPbXF/IzqtdcqSrJnw2Ok3PrrjUllg8miktZONh72vleeteWNVYWcvQ3J2xBsnTmgOovCBfAYAB9vzu8cr7llVPUUgPsAPLJYdxYclT+aeQPytC4LVjVMYJOxEn4ypcPgte0tckv6bHJnqYnI/+qmYH1mZZgPbBwr484IJPNjjr6tSig9qayKP2zD8SIITyHb8tP9bE5n82dGFP9odqGI3NqcX6uq1+6RS3uGo0K4KUTkagBXj6dffMfWmz/oFrTzfKaNumTLb/tCAHfvuB/Z6vV5N308U/dC6OjnfrVdtbHz/pkx3Y3dOeG3NfhZ8SObt56NueaQnadn8N91WHHxedx70zv0+gsLRe9W1cuC+58A8Ljm/LHjNa/MnSJyHMDDAXymWHcWHBXCTQdk3O2uBQARuVVVL9k/91bHUfH1qPgJHB1fj4qfwN75mpBoD24BcLGIXISBG64A8F2mzA0ArgTw+wCeD+BdqqoicgOAXxeRVwD4cgAXA/jDmfzahqNCuJXBXLBgwVkKVT0lIi8GcBOADQCvVdXbRORlAG5V1RsAvAbA60TkBIB7MPAIxnJvAvAhAKcAvEhV9+TNwUeCcNlgHrBbCxYsOERQ1RsB3GiuvbQ5vh/AC0jdfw3gX++pgzgihAv4gxngKCXTj4qvR8VP4Oj4elT8BI6Wr4cWonO9vWPBggULFoQ4Ku/DXbBgwYIjj7Uj3O6PAO+tL48TkXeLyIdE5DYR+YHx+gUi8nYR+cj4/yPG6yIirxp9f7+IPGWf/d0QkfeKyNvG84vGz5yfGD+Dfu54nX4mfZ/8PF9ErheRPxGRD4vINx7iMf2h8dl/UEReLyIPOgzjKiKvFZG7ROSDzbXuMRSRK8fyHxGRK/fK37WBqq7NC8Mf1P4UwFcCOBfAHwN44gH682gATxmPvxTAf8Pw0eSfBnDNeP0aAC8fj58L4D8CEABPB3DzPvv7wwB+HcDbxvM3AbhiPP5FAP/nePxPAPzieHwFgDfus5/XAfjfx+NzAZx/GMcUwwd2Pgrgwc14/sPDMK4AvhnAUwB8sLnWNYYALgBwx/j/I8bjR+znXDhqrwN3YOZJ9I0AbmrOXwLgJQftV+PPWwE8C8DtAB49Xns0gNvH418C8MKm/Oly++DbYwG8E8AzAbxtXFx3AzhuxxbDu0W+cTw+PpaTffLz4SOJibl+GMd0+oTkBeM4vQ3Atx2WcQXweEO4XWMI4IUAfqm5vq3c8tr5WreUwr59JroXY3j4ZAA3A3iUqn5yvPUpAI8ajw/S/58F8M8BTB9cfySAz+rwmXPrC/tM+n7gIgCfBvArY/rjl0XkoTiEY6qqnwDwMwA+BuCTGMbpPTic4wr0j+GhXW+HFetGuIcSIvIlAH4DwA+q6ufaezpIgwN9q4iIfDuAu1T1PQfpRxHHMYTCv6CqTwbwVxjC39M4DGMKAGMO9HIMm8SXA3gohq//O/Q4LGO4blg3wt23z0RXISLnYCDbX1PVt4yX/1JEHj3efzSAu8brB+X/MwB8h4j8GYZvYnsmgJ8DcP74mXPry2k/zWfS9wN3ArhTVW8ez6/HQMCHbUwB4FsBfFRVP62qJwG8BcNYH8ZxBfrH8NCtt8OOdSPc0x8BHv/yewWGz08fCEREMHyc8MOq+orm1vSZboz/v7W5/r3jX4WfDuC+JsTbM6jqS1T1sar6eAxj9i5V/W4A78bwmXPPz8n/059J32s/R18/BeDjIvLV46VLMXwk81CN6YiPAXi6iDxknAuTr4duXJ32K2N4E4Bni8gjRjX/7PHaAoaDTiLP/cLwF9X/huHdCv/3AfvydzGEZe8H8L7x9VwMebl3AvgIgHcAuGAsLxi+aP1PAXwAwCUH4PO34My7FL4Sw5d4nADwZgDnjdcfNJ6fGO9/5T77+LcB3DqO629h+Av5oRxTAD8O4E8AfBDA6wCcdxjGFcDrMeSVT2KIGq5aZQwB/KPR3xMAvm+/5+tRey2fNFuwYMGCfcK6pRQWLFiw4NBiIdwFCxYs2CcshLtgwYIF+4SFcBcsWLBgn7AQ7oIFCxbsExbCXbCvEJFNEXnf+O1Zvy0i55v7Pygi94vIww/KxwUL9goL4S7Yb/yNqv5tVf1aDL8r9SJz/4UYPsDyv+y7ZwsW7DEWwl1wkPh9NF92IiJfBeBLAPxLDMS7YMFaYSHcBQcCEdnA8FHX9qPXV2D4Lof/AuCrReRRXt0FC44qFsJdsN94sIi8D2e+/u/tzb0XAniDqm5h+MIf9xdWFyw4qlg+2rtgXyEiX1DVLxGRh2D4opM3q+qrRORJGL4fYfpimXMxfNPWMw7K1wUL5sZCuAv2FRPhjsdPxvDlM18F4GUAPq+qP9mU/SiAb1HVPz8QZxcsmBlLSmHBgUFV34vhG79eiCF/+5umyG+O1xcsWAssCnfBggUL9gmLwl2wYMGCfcJCuAsWLFiwT1gId8GCBQv2CQvhLliwYME+YSHcBQsWLNgnLIS7YMGCBfuEhXAXLFiwYJ+wEO6CBQsW7BP+f3TLn1oLSOqnAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "### This code shows the 2D image from a fits file\n", "### Need to have astropy, WCSAxes packages installed\n", "imname =input('Enter the fits file name (these are 1hr.fits, 2hrs.fits, 5hrs.fits, 10hrs.fits): ')\n", "hdulist =fits.open(os.path.join(data_dir, imname))\n", "hdu =hdulist[0]\n", "data =hdu.data[0,0,:,:]\n", "hdu.header['CDELT3'] =1.0\n", "plt.subplot()\n", "plt.imshow(hdu.data[0,0,:,:],vmax=max(hdu.data[0,0,:,:].flatten()),vmin=min(hdu.data[0,0,:,:].flatten()),origin='lower')\n", "plt.xlabel('RA')\n", "plt.ylabel('Dec')\n", "cb=plt.colorbar()\n", "cb.set_label('Flux (Jy)')\n", "hdulist.close()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Image Statistics\n", "Give the pixels corresponding to the box in which you want to find out the flux and obtain the flux. The flux will be estimated for the last running file. Run the code below and give the pixel range obtained from the other image. \n", "Find the flux density for the noise in the upper right region, region covering the\n", "bright source, and region covering the extended source. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "bottom left corner of your region (in format, e.g. 502, 506): 502, 506\n", "upper right corner of your region (in format, e.g. 520, 518): 520, 518\n", "502 506\n", "520 518\n", "[[0.01005777 0.01000005 0.01001082 0.01007294 0.01015488 0.01022022\n", " 0.01023758 0.01018868 0.01007239 0.00990365 0.00970807 0.00951426\n", " 0.00934673]\n", " [0.01023271 0.01019404 0.01020637 0.01026021 0.01033302 0.01039596\n", " 0.0104222 0.01039421 0.01030773 0.01017102 0.01000033 0.00981409\n", " 0.00962825]\n", " [0.01065369 0.01070555 0.01073378 0.01073424 0.01070292 0.01063723\n", " 0.01053757 0.01040852 0.01025844 0.0100973 0.00993328 0.00977031\n", " 0.00960789]\n", " [0.01106416 0.01119297 0.01125261 0.01123405 0.0111398 0.01098323\n", " 0.01078549 0.01057081 0.01036125 0.01017184 0.0100074 0.00986225\n", " 0.00972322]\n", " [0.01130496 0.0114346 0.01152207 0.01155037 0.0115135 0.01141792\n", " 0.01127947 0.01111779 0.01095018 0.01078678 0.01062831 0.01046681\n", " 0.01028912]\n", " [0.01196422 0.01219954 0.01243916 0.01264996 0.01280291 0.01287939\n", " 0.01287272 0.01278592 0.01262754 0.01240744 0.01213369 0.01181115\n", " 0.01144209]\n", " [0.01476999 0.01581316 0.0168462 0.01774593 0.01839303 0.01869882\n", " 0.01862469 0.0181883 0.01745541 0.01652069 0.01548388 0.01442908\n", " 0.01341264]\n", " [0.02239434 0.02583338 0.02917733 0.03201284 0.03393997 0.03466692\n", " 0.03408504 0.03229723 0.02958909 0.02635303 0.0229924 0.01983726\n", " 0.01709529]\n", " [0.03624028 0.04404863 0.05166366 0.05811083 0.06244247 0.06396792\n", " 0.06243815 0.05811663 0.05170812 0.04417152 0.03648453 0.02943926\n", " 0.02352492]\n", " [0.05187708 0.06457129 0.07700859 0.08756473 0.0946558 0.09712762\n", " 0.09456699 0.08741893 0.07686369 0.06449973 0.05194939 0.04052015\n", " 0.03101415]\n", " [0.05915116 0.07409479 0.08876162 0.101214 0.10956524 0.11244636\n", " 0.10937542 0.10089201 0.08840425 0.07380734 0.05901509 0.04556262\n", " 0.03438775]\n", " [0.05212442 0.06487479 0.07736365 0.08793212 0.09497309 0.09733023\n", " 0.09460969 0.08728862 0.0765781 0.06409404 0.05145862 0.0399667\n", " 0.0304064 ]\n", " [0.03646803 0.04435004 0.05201098 0.05843535 0.06264925 0.06396261\n", " 0.06215573 0.05754007 0.05086955 0.04313643 0.03532794 0.02822574\n", " 0.02230119]\n", " [0.02242391 0.02595528 0.02932823 0.03210371 0.03387016 0.03434619\n", " 0.03345735 0.0313575 0.02838579 0.02497528 0.02154705 0.01842565\n", " 0.01579834]\n", " [0.01455818 0.01572103 0.01680105 0.01766706 0.0182039 0.01834147\n", " 0.01807346 0.01745822 0.01660094 0.01562383 0.01463523 0.01370848\n", " 0.01287605]\n", " [0.01145947 0.01179292 0.01209312 0.0123312 0.01248991 0.01256774\n", " 0.01257705 0.01253692 0.01246392 0.0123649 0.01223476 0.01205973\n", " 0.01182476]\n", " [0.01059492 0.01076944 0.01090319 0.01098621 0.01102456 0.01103797\n", " 0.0110507 0.01107992 0.01112698 0.01117585 0.01119913 0.01116871\n", " 0.01106637]\n", " [0.01042185 0.01060927 0.01072293 0.01075357 0.01071509 0.01063946\n", " 0.01056428 0.01051825 0.01051125 0.01053297 0.01055971 0.01056553\n", " 0.01053244]\n", " [0.0102014 0.01040831 0.01054295 0.01059606 0.01057942 0.01052111\n", " 0.01045479 0.01040781 0.01039292 0.01040662 0.01043344 0.01045341\n", " 0.01044915]]\n", "352.09133750513985 mJy\n" ] } ], "source": [ "bottom_x, bottom_y = map(int, input('bottom left corner of your region (in format, e.g. 502, 506): ').split(','))\n", "upper_x, upper_y = map(int, input('upper right corner of your region (in format, e.g. 520, 518): ').split(','))\n", "\n", "print (bottom_x,bottom_y)\n", "print (upper_x,upper_y)\n", "\n", "num_pix_in_beam=19 ###beam_area/pixel_size**2\n", "print (data[bottom_x:upper_x+1,bottom_y:upper_y+1])\n", "flux_density=np.sum(data[bottom_x:upper_x+1,bottom_y:upper_y+1])/num_pix_in_beam\n", "print (flux_density*1e3, \"mJy\")" ] } ], "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.7.1" } }, "nbformat": 4, "nbformat_minor": 2 }