astropy_related.ipynb

{
 "cells": [
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# usage of astropy functions\n\n## 1. astropy.cosmoloy\nastropy.cosmology provides different cosmological models for conversion between redshift(z) and related parameters (comoving distance, luminosity distance etc.). User can also modify existing models' parameters(Hubble constant H0...). The demonstration below use Plack15 as the cosmological model."
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "import matplotlib.pyplot as plt\n## note that Planck15 is a class FLRW object\nfrom astropy.cosmology import Planck15 as cosmo\nimport numpy as np\nimport seaborn as sns\nsns.set_style('whitegrid')\nsns.set_context('notebook')\nimport pandas as pd",
   "execution_count": 16,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "## Suppose you have an array of redsift data\nz = np.arange(1.0e-8, 0.2, 0.001)\n\n## obtain parameters using cosmo\ndl = cosmo.luminosity_distance(z)\nage = cosmo.age(z)\nd_co = cosmo.comoving_distance(z)\n\nprint(dl)\nprint(age)\nprint(d_co)",
   "execution_count": 24,
   "outputs": [
    {
     "output_type": "stream",
     "text": "[4.42563419e-05 4.42907719e+00 8.86490260e+00 1.33075131e+01\n 1.77569012e+01 2.22130596e+01 2.66759809e+01 3.11456577e+01\n 3.56220825e+01 4.01052480e+01 4.45951469e+01 4.90917716e+01\n 5.35951148e+01 5.81051692e+01 6.26219273e+01 6.71453817e+01\n 7.16755252e+01 7.62123502e+01 8.07558494e+01 8.53060155e+01\n 8.98628411e+01 9.44263188e+01 9.89964412e+01 1.03573201e+02\n 1.08156591e+02 1.12746603e+02 1.17343231e+02 1.21946466e+02\n 1.26556302e+02 1.31172732e+02 1.35795747e+02 1.40425340e+02\n 1.45061505e+02 1.49704234e+02 1.54353519e+02 1.59009353e+02\n 1.63671729e+02 1.68340639e+02 1.73016077e+02 1.77698034e+02\n 1.82386504e+02 1.87081478e+02 1.91782951e+02 1.96490914e+02\n 2.01205360e+02 2.05926282e+02 2.10653673e+02 2.15387525e+02\n 2.20127831e+02 2.24874583e+02 2.29627775e+02 2.34387398e+02\n 2.39153446e+02 2.43925912e+02 2.48704788e+02 2.53490067e+02\n 2.58281741e+02 2.63079804e+02 2.67884247e+02 2.72695064e+02\n 2.77512248e+02 2.82335791e+02 2.87165685e+02 2.92001925e+02\n 2.96844502e+02 3.01693408e+02 3.06548638e+02 3.11410184e+02\n 3.16278037e+02 3.21152192e+02 3.26032641e+02 3.30919376e+02\n 3.35812391e+02 3.40711678e+02 3.45617230e+02 3.50529040e+02\n 3.55447100e+02 3.60371403e+02 3.65301943e+02 3.70238711e+02\n 3.75181702e+02 3.80130906e+02 3.85086318e+02 3.90047930e+02\n 3.95015735e+02 3.99989725e+02 4.04969894e+02 4.09956235e+02\n 4.14948739e+02 4.19947401e+02 4.24952212e+02 4.29963166e+02\n 4.34980255e+02 4.40003473e+02 4.45032812e+02 4.50068265e+02\n 4.55109825e+02 4.60157484e+02 4.65211236e+02 4.70271074e+02\n 4.75336990e+02 4.80408977e+02 4.85487028e+02 4.90571136e+02\n 4.95661294e+02 5.00757495e+02 5.05859732e+02 5.10967996e+02\n 5.16082283e+02 5.21202584e+02 5.26328892e+02 5.31461200e+02\n 5.36599501e+02 5.41743788e+02 5.46894055e+02 5.52050293e+02\n 5.57212496e+02 5.62380657e+02 5.67554768e+02 5.72734823e+02\n 5.77920815e+02 5.83112736e+02 5.88310580e+02 5.93514340e+02\n 5.98724007e+02 6.03939577e+02 6.09161041e+02 6.14388392e+02\n 6.19621623e+02 6.24860728e+02 6.30105700e+02 6.35356531e+02\n 6.40613214e+02 6.45875743e+02 6.51144110e+02 6.56418309e+02\n 6.61698332e+02 6.66984173e+02 6.72275825e+02 6.77573280e+02\n 6.82876532e+02 6.88185573e+02 6.93500397e+02 6.98820998e+02\n 7.04147367e+02 7.09479498e+02 7.14817384e+02 7.20161019e+02\n 7.25510394e+02 7.30865505e+02 7.36226342e+02 7.41592900e+02\n 7.46965172e+02 7.52343150e+02 7.57726829e+02 7.63116200e+02\n 7.68511258e+02 7.73911995e+02 7.79318404e+02 7.84730479e+02\n 7.90148213e+02 7.95571598e+02 8.01000628e+02 8.06435297e+02\n 8.11875597e+02 8.17321522e+02 8.22773064e+02 8.28230217e+02\n 8.33692975e+02 8.39161329e+02 8.44635274e+02 8.50114803e+02\n 8.55599909e+02 8.61090585e+02 8.66586824e+02 8.72088620e+02\n 8.77595966e+02 8.83108855e+02 8.88627280e+02 8.94151234e+02\n 8.99680712e+02 9.05215706e+02 9.10756208e+02 9.16302214e+02\n 9.21853715e+02 9.27410706e+02 9.32973179e+02 9.38541128e+02\n 9.44114546e+02 9.49693427e+02 9.55277763e+02 9.60867548e+02\n 9.66462775e+02 9.72063438e+02 9.77669531e+02 9.83281045e+02\n 9.88897975e+02 9.94520314e+02 1.00014806e+03 1.00578119e+03] Mpc\n[13.79761652 13.78319259 13.76878974 13.75440795 13.74004716 13.72570735\n 13.71138848 13.6970905  13.68281338 13.66855709 13.65432158 13.64010682\n 13.62591277 13.61173939 13.59758665 13.58345451 13.56934294 13.55525189\n 13.54118133 13.52713122 13.51310153 13.49909222 13.48510325 13.47113459\n 13.4571862  13.44325805 13.4293501  13.4154623  13.40159464 13.38774707\n 13.37391956 13.36011206 13.34632455 13.33255699 13.31880935 13.30508158\n 13.29137366 13.27768555 13.26401721 13.25036861 13.23673972 13.22313049\n 13.2095409  13.19597091 13.18242049 13.16888959 13.1553782  13.14188626\n 13.12841376 13.11496065 13.1015269  13.08811248 13.07471735 13.06134147\n 13.04798483 13.03464737 13.02132908 13.0080299  12.99474982 12.9814888\n 12.9682468  12.95502379 12.94181974 12.92863461 12.91546838 12.90232101\n 12.88919246 12.8760827  12.86299171 12.84991945 12.83686588 12.82383097\n 12.8108147  12.79781702 12.78483791 12.77187734 12.75893526 12.74601166\n 12.7331065  12.72021974 12.70735135 12.69450131 12.68166958 12.66885613\n 12.65606093 12.64328394 12.63052514 12.61778449 12.60506197 12.59235753\n 12.57967116 12.56700282 12.55435247 12.54172009 12.52910565 12.51650911\n 12.50393044 12.49136963 12.47882662 12.4663014  12.45379393 12.44130418\n 12.42883212 12.41637773 12.40394097 12.3915218  12.37912021 12.36673616\n 12.35436963 12.34202057 12.32968896 12.31737478 12.30507799 12.29279856\n 12.28053646 12.26829167 12.25606415 12.24385387 12.23166081 12.21948494\n 12.20732622 12.19518463 12.18306014 12.17095271 12.15886233 12.14678896\n 12.13473257 12.12269314 12.11067063 12.09866502 12.08667628 12.07470438\n 12.06274928 12.05081097 12.03888942 12.02698459 12.01509646 12.00322499\n 11.99137017 11.97953196 11.96771034 11.95590528 11.94411674 11.93234471\n 11.92058915 11.90885003 11.89712734 11.88542103 11.87373109 11.86205749\n 11.85040019 11.83875918 11.82713441 11.81552588 11.80393354 11.79235738\n 11.78079736 11.76925346 11.75772564 11.7462139  11.73471819 11.72323849\n 11.71177477 11.70032701 11.68889518 11.67747925 11.6660792  11.654695\n 11.64332663 11.63197405 11.62063724 11.60931618 11.59801083 11.58672118\n 11.5754472  11.56418885 11.55294612 11.54171897 11.53050739 11.51931135\n 11.50813081 11.49696576 11.48581617 11.47468201 11.46356326 11.45245989\n 11.44137187 11.43029919 11.41924181 11.40819971 11.39717287 11.38616126\n 11.37516485 11.36418362 11.35321754 11.34226659 11.33133075 11.32040998\n 11.30950427 11.29861359] Gyr\n[4.42563415e-05 4.42465249e+00 8.84720809e+00 1.32677098e+01\n 1.76861564e+01 2.21025467e+01 2.65168794e+01 3.09291533e+01\n 3.53393672e+01 3.97475200e+01 4.41536103e+01 4.85576371e+01\n 5.29595991e+01 5.73594951e+01 6.17573241e+01 6.61530848e+01\n 7.05467760e+01 7.49383967e+01 7.93279456e+01 8.37154217e+01\n 8.81008237e+01 9.24841507e+01 9.68654014e+01 1.01244575e+02\n 1.05621670e+02 1.09996685e+02 1.14369619e+02 1.18740472e+02\n 1.23109242e+02 1.27475928e+02 1.31840530e+02 1.36203045e+02\n 1.40563473e+02 1.44921813e+02 1.49278063e+02 1.53632224e+02\n 1.57984293e+02 1.62334270e+02 1.66682153e+02 1.71027942e+02\n 1.75371636e+02 1.79713234e+02 1.84052734e+02 1.88390136e+02\n 1.92725439e+02 1.97058642e+02 2.01389743e+02 2.05718742e+02\n 2.10045638e+02 2.14370430e+02 2.18693117e+02 2.23013697e+02\n 2.27332171e+02 2.31648537e+02 2.35962795e+02 2.40274943e+02\n 2.44584980e+02 2.48892905e+02 2.53198719e+02 2.57502419e+02\n 2.61804005e+02 2.66103476e+02 2.70400831e+02 2.74696070e+02\n 2.78989191e+02 2.83280193e+02 2.87569076e+02 2.91855839e+02\n 2.96140482e+02 3.00423002e+02 3.04703400e+02 3.08981674e+02\n 3.13257825e+02 3.17531850e+02 3.21803749e+02 3.26073522e+02\n 3.30341168e+02 3.34606685e+02 3.38870074e+02 3.43131333e+02\n 3.47390461e+02 3.51647458e+02 3.55902324e+02 3.60155057e+02\n 3.64405656e+02 3.68654121e+02 3.72900452e+02 3.77144647e+02\n 3.81386705e+02 3.85626627e+02 3.89864411e+02 3.94100057e+02\n 3.98333563e+02 4.02564930e+02 4.06794157e+02 4.11021243e+02\n 4.15246186e+02 4.19468988e+02 4.23689647e+02 4.27908162e+02\n 4.32124532e+02 4.36338758e+02 4.40550838e+02 4.44760772e+02\n 4.48968560e+02 4.53174200e+02 4.57377692e+02 4.61579035e+02\n 4.65778229e+02 4.69975274e+02 4.74170168e+02 4.78362912e+02\n 4.82553504e+02 4.86741944e+02 4.90928231e+02 4.95112366e+02\n 4.99294347e+02 5.03474173e+02 5.07651845e+02 5.11827362e+02\n 5.16000723e+02 5.20171928e+02 5.24340976e+02 5.28507867e+02\n 5.32672600e+02 5.36835175e+02 5.40995591e+02 5.45153848e+02\n 5.49309945e+02 5.53463882e+02 5.57615659e+02 5.61765274e+02\n 5.65912728e+02 5.70058021e+02 5.74201150e+02 5.78342117e+02\n 5.82480921e+02 5.86617561e+02 5.90752038e+02 5.94884349e+02\n 5.99014496e+02 6.03142478e+02 6.07268294e+02 6.11391944e+02\n 6.15513427e+02 6.19632744e+02 6.23749894e+02 6.27864876e+02\n 6.31977690e+02 6.36088336e+02 6.40196814e+02 6.44303122e+02\n 6.48407262e+02 6.52509231e+02 6.56609031e+02 6.60706661e+02\n 6.64802120e+02 6.68895409e+02 6.72986526e+02 6.77075472e+02\n 6.81162246e+02 6.85246849e+02 6.89329278e+02 6.93409536e+02\n 6.97487620e+02 7.01563532e+02 7.05637270e+02 7.09708835e+02\n 7.13778226e+02 7.17845442e+02 7.21910485e+02 7.25973353e+02\n 7.30034046e+02 7.34092564e+02 7.38148907e+02 7.42203075e+02\n 7.46255067e+02 7.50304883e+02 7.54352523e+02 7.58397987e+02\n 7.62441275e+02 7.66482386e+02 7.70521320e+02 7.74558078e+02\n 7.78592659e+02 7.82625062e+02 7.86655288e+02 7.90683336e+02\n 7.94709207e+02 7.98732900e+02 8.02754416e+02 8.06773753e+02\n 8.10790912e+02 8.14805893e+02 8.18818695e+02 8.22829319e+02\n 8.26837765e+02 8.30844032e+02 8.34848120e+02 8.38850029e+02] Mpc\n",
     "name": "stdout"
    }
   ]
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "## More useful parameters:\nx4 = cosmo.H(z)\nx5 = cosmo.angular_diameter_distance(z)\nx6 = cosmo.comoving_volume(z)\n\n## note that all functions assumed the observer positioned\n## at z = 0, change the setting by including initial\n## redshift z0 in the bracket(z0, z)\n\n## for the function distmod, z cannot include 0 (it will\n## lead to a log(0) error)",
   "execution_count": 3,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "fig, ax = plt.subplots(3, figsize=(6, 12), squeeze=True, sharex=True)\nax[0].plot(z, dl, \"-\")\nax[0].set_ylabel(\"d_L(Mpc)\")\nax[0].set_title(\"redshift against luminosity distance\")\nax[0].set_xlim(left=0)\nax[0].set_ylim(bottom=0)\n",
   "execution_count": 25,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 25,
     "data": {
      "text/plain": "(0, 1056.0702509599225)"
     },
     "metadata": {}
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x864 with 3 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAK8CAYAAAAqDlIfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XtAzXfjB/B395tIdETItaIQYRjFXAo11mwuIzZjV8azmUbkOmY9chvD49meeXgws8JDuY2hbGQbpZRLueV0KlKp0+mcz+8PP+dZQ6pv9T3V+/XP+l7O+b472Xmf7/l+v5+vkRBCgIiISAJjuQMQEVHNxzIhIiLJWCZERCQZy4SIiCRjmRARkWQsEyIikoxlQpItWrQIa9euLXWdl156CRcvXnxi/sWLFzF9+nQAQHp6Ovz9/TFixAj89ttveOutt5CdnV0lmf9s9erViIiIqNBjc3NzERQU9NRlwcHB2LJli5Roen9+nSrT0aNHsWTJEgDA8ePHsXr16go/V3Z2NlxdXZ943meRuj0yLKZyB6C6rVOnTlizZg0A4JdffkHjxo3x7bffAgBOnz5dLRk++uijCj82JyfnqSVZ2f78OlWmgQMHYuDAgQAeFVZOTk6lP++zVOb2SH4sEwLw6I186dKlsLa2Rn5+Pn744QecOnUKGzZsgEajgaWlJWbPno2uXbsiLy8Pc+fORVJSEhQKBUxMTODl5QUA2L59O3bs2AEzMzNYWFhg0aJFaNeuHQBg586dCA0NRXZ2NkaMGIGZM2fil19+weLFixESEoJVq1YhNzcXEyZMQPPmzQEAEydOxKZNm9C0aVN91szMTMyfPx9ZWVlQqVRwcnLCqlWr0KhRI1y4cAELFiyARqNBy5YtcefOHQQHB6NHjx74/PPP8ccffyA/Px9CCCxZsgReXl4IDg5G+/btMXnyZHTq1AlTp07F6dOnkZGRgbfffhvjxo2DSqXC7Nmzce/ePQCAj48PZsyYgc8++wyFhYUYMWIE9uzZAxMTk6e+vq6uroiNjYW9vX2J6ZSUFKxcuRJNmzbF9evXYWVlhalTp2Lr1q24fv06hgwZgjlz5uhfp/379yM4OBj16tXD5cuXcffuXbi6uuKLL76AjY0Nzp07hxUrVqCgoABmZmaYMWMGvL29n5l/z549iI6Oxvvvv48dO3ZAq9XC1tYWFy5cwNChQ/H6668DANavX4/79+9jzpw5JX6vQ4cOITw8HFZWVvDw8NDPf/y8GzduxKFDh7BhwwYYGRnBxMQEn376KczNzUts75133sGCBQuQlpaG+/fvw8bGBmFhYWjTpg0mTJgAT09PnD9/Hunp6ejduzcWL14MY2Nj/PTTT1i1ahV0Oh2sra2xcOFCuLm54fz58wgLC0NBQQGMjY3x4YcfYsCAAZL/P6FSCCIhxJkzZ4Sbm5u4deuWEEKI69evC39/f5GdnS2EECI5OVm8+OKLIj8/XyxdulR8+umnQqfTiaysLOHt7S3WrFkjiouLhbu7u1AqlUIIIX788UexY8cOIYQQAwYMEIsWLRJCCJGRkSE8PDzEnTt3xJkzZ8Tw4cOFEEL88MMPYurUqfpMLi4uIisr64ms3377rdi4caMQQgidTifefvttsWXLFqHRaIS3t7c4fvy4EEKI2NhY4erqKs6cOSPOnz8vpk2bJrRarRBCiI0bN4p33nlHCCHE7NmzxT/+8Q/9Nrdu3SqEEOLixYvCw8NDFBYWinXr1ol58+YJIYTIz88XM2bMEA8ePBA3b94Unp6eT31N//q8f/5dHk+fOXNGdOjQQSQkJAghhJg8ebIYPXq0UKvVIisrS7i7u4u7d++WeJ1mz56tX6eoqEiMHDlS7N69W2RnZ4vevXuL33//Xf8369mzp7hx48Yz8//5NV+zZo1YuHChEEKIw4cPi1dffVUIIYRWqxUDBgwQV69eLfH7qVQq4eXlJVJSUoQQQnz99dfCxcXlib/lwIEDxW+//SaEEOLkyZNi7dq1T2zv4MGDYvHixfrnnjdvnv7fy/jx48X06dOFVqsVubm5om/fviI2Nla//cevXXR0tJg8ebK4f/++GDJkiLh586YQQoi7d+8Kb29vcfv27af+nahycM+E9Jo2bQonJycA0H8ynzRpkn65kZERbty4gdjYWMyZMwdGRkawt7fH4MGDAQAmJibw8/PDmDFj0L9/f/Tt2xc+Pj76x/v7+wMAHBwc0LhxY2RlZVUo58SJE3Hu3Dl88803SE1NRUpKCrp06YLk5GQA0G+zV69eaN++PQCga9euaNCgAXbs2IGbN2/il19+gY2NzVOf//HXM+7u7igqKsLDhw/Rr18/TJ06Fenp6ejTpw8+/vhj2NraVsrXNM2bN0fHjh0BAC1btoStrS3Mzc1hb28PGxubp26jX79+MDc3BwC4uLggJycHFy5cQMuWLdGlSxcAQPv27dGtWzf8+uuvz8z/LAMGDMDSpUuRlJQEpVKJ5s2bo02bNiXWiYuLg4uLi37Pc/To0Vi5cuUTzzV8+HB8+OGH8PHxwYsvvogpU6Y8sY6fnx9atGiBrVu3Ii0tDb/++iu6du1aIo+xsTHq1asHZ2dn5OTk4Pz582jfvr3+tRsyZAiGDBmCEydOQKVS4YMPPtA/3sjICJcvX0azZs2e+TuTNCwT0rO2ttb/rNPp0Lt3b6xatUo/Lz09HQqFAgAg/jSk25+/2gkLC0NycjJiYmKwadMmREZG6g+ympr+75+bkZFRiecojy+//BIXLlzAq6++ihdeeAHFxcUQQsDExOSJ53yc7fjx41i6dCnefPNNDBw4EG3atMHevXuf+vwWFhb6jI9/186dO+Po0aOIjY3FmTNn8Nprr2Hz5s2ws7Mrd/6ioqIS049L4bE/v07PYmlpqf/58Wup1Wr1mR8TQqC4uPiZ+Z/FxMQEo0ePxu7du5GRkYExY8Y8db0/v97Pyj1z5ky8+uqrOH36NPbs2YN//vOf2L17d4l1tm/fjl27duGNN95AQEAA7OzscOvWrVJ/X1NT0xK/rxACly9fhlarRdu2bfH999/rlymVSv1XjFQ1eDYXPVXv3r1x+vRpXL16FQBw4sQJvPzyyygsLES/fv2we/du6HQ65OTk4OjRowAenc3j4+MDOzs7TJo0CTNmzJB0cNrExATFxcVPzD916hQmTpyIkSNHolGjRoiJidG/gZibm+Pnn38GAFy4cAHJyckwMjLC6dOnMWDAAIwbNw4eHh44cuQItFptmbOEhYVh/fr1GDRoEObOnYt27dohJSUFpqam0Gq1zy1Ge3t7/Wuxf//+crwKZefp6Ylr167hwoULAICUlBScPXsWPXv2fGb+P/vr6/3aa6/hyJEjSEhI0O99/lmPHj1w5coVJCUlAXh0nOSviouL8dJLL6GgoABjx45FaGgoLl++jKKiohLbO3XqFF555RW89tpraN26NY4dO/bcv0+XLl1w9epV/e9x9OhRzJo1C56enkhLS8PZs2cBAImJifD19YVSqSzrS0kVwD0Teqp27dph0aJF+Nvf/qb/FLhhwwbY2Nhg2rRpCA0NxdChQ2Fvbw8XFxcAj94w33vvPUyaNAmWlpYwMTF57umhpfHz88OECROwdu1a/TYA4IMPPsCKFSuwevVqmJmZoVu3brhx4wZMTU2xdu1ahIaGYuXKlWjVqhUaN24MS0tLjBkzBh9//DECAgJQXFyMF198EYcOHYJOpytTlokTJyI4OBj+/v4wNzeHq6srhg8fDhMTE3Tu3BnDhw/Htm3b0LBhw6c+PiQkBIsWLUL9+vXRp08fODg4VPh1eRZ7e3usXr0aixcvRmFhIYyMjLBs2TK0bt36mfn/XGy9evXCJ598gsWLF2PevHlo1KgRPDw80LZtW5iZmT11e2FhYfjkk09gZmaGHj16PLGOqakp5syZg08++US/J/H555/D3Ny8xPbeeustzJ8/X7/H4unpqf/a8lkaN26MsLAwzJ49G1qtFvXq1UN4eDjs7e2xZs0arFixAmq1GkIIrFixQn9SB1UNI1HR7xqIDNAXX3yByZMno3HjxkhPT8eIESNw5MgR1K9fX+5oNU52djZGjRqFbdu2lTibjuhpuGdCtYqTkxMmTZoEU1NT/em/LJLy27VrF1auXIlp06axSKhMuGdCRESS8QA8ERFJxjIhIiLJau0xE51Oh/z8fJiZmT1x7j0RET2dEAIajQY2NjYwNi77/katLZP8/PznnlpIRERP5+LiUuooCX9Va8vk8XnxLi4uT1xhbCji4+NLDI5naJhPGuaTztAz1sZ8RUVFSE5Ofuq1RaWptWXy+Kstc3Nz/fAYhsiQswHMJxXzSWfoGWtrvvIeHuABeCIikoxlQkREkrFMiIhIsiotk7y8PPj7++uHko6JiUFAQACGDBmC8PBw/XqJiYkIDAyEr68v5s6dqx9J9M6dO3jjjTfg5+eH9957D/n5+VUZl4iIKqjKyuSPP/7A2LFjkZqaCgAoLCzEnDlzsH79ehw4cADx8fE4ceIEAGDWrFmYP38+oqOjIYTArl27AAALFy7EuHHjEBUVBQ8PD6xfv76q4hIRkQRVVia7du1CaGio/mZKFy5cgLOzM1q0aAFTU1MEBAQgKioKt2/fRmFhITw9PQEAgYGBiIqKgkajwdmzZ+Hr61tiPhERPZ1WJ7Dz7A1Exd+t9m1X2anBS5cuLTGdkZFR4h4OCoUCSqXyifkODg5QKpW4d+8e6tWrp7972+P5RET0pGRlLoJ/uIDzN+5jct/W8PNwrNbtV9t1Jjqd7olbbBoZGT1z/uP//llFhkWJj4+veOhqEBcXJ3eEUjGfNMwnnaFnlDtfkVbgh8Q8RCTlw8rMCNN7NoC340N9rurKV21l4ujoCJVKpZ9WqVRQKBRPzM/MzIRCoYC9vT1yc3Oh1WphYmKiX7+8PDw8DPaiori4OHh5eckd45mYTxrmk87QM8qdL/ZqFub8eBHXM/MR2M0JIcM7wt7mfyN+VCSfWq2u0Ifwajs1uEuXLrh+/TrS0tKg1Wqxf/9+eHt7w8nJCRYWFvr2jIyMhLe3N8zMzNC9e3ccOHAAABAREQFvb+/qiktEZLDuPyzCp7v/wNjNZ6DVCfx78gtY+bpniSKpbtW2Z2JhYYHly5dj2rRpUKvV8PHxgZ+fHwAgLCwMISEhyMvLg7u7O4KCggAAoaGhCA4OxoYNG9C0aVOsXLmyuuISERkcIQT2/nEHi/dfwr2HGrzr0xYfDWwPK3MTuaNVfZkcO3ZM/3Pv3r2xd+/eJ9Zxc3PD7t27n5jv5OSErVu3Vmk+IqKa4Gb2Q8yLjMfxyyp0ad4A3731Ajo2M5xbUtfagR6JiGqDYq0O38ak4u+HkmFkBMz374iJfVrBxNiw7tPEMiEiMlDxt3MQvOcC4m8/wEtuCiwe6QEnOyu5Yz0Vy4SIyMA8LCrGqiMp2HLqOhpam+Orcd0wrJOjQd81lmVCRGRATiSrMPfHi7h1rwBje7ZEsJ8bGliX70ZVcmCZEBEZgMw8NRbvv4TI3++grYMNdr3TGz1b28sdq8xYJkREMhJC4Pu4W/j8QCLy1cX4aGB7vD+gLSxM5T/dtzxYJkREMrmemY85ey4i9loWerRqiGWBndBOYSt3rAphmRARVbOiYh02n7yG1UdTYGFqjM9f6YQxPVrA2MBO9y0PlgkRUTU6f+MePvvhIi4rczGskyMWBLhDUd9S7liSsUyIiKpBbqEGX0ZfxtYzaXCsb4nNQd0xuGMTuWNVGpYJEVEVO5RwF/MjE6DMLcTE3q3wia8r6lnUrrff2vXbEBEZkDv3C7BgbwIOXVLCzdEWG8Z3Q9eWDeWOVSVYJkRElaxYq8O/YtOw8tBlaIVA8FA3TO7bGmYm1XbXj2rHMiEiqkQXbt3HZ3suIuHOAwxwdcCiER5oYW8td6wqxzIhIqoEuYUa/P1QMr6LTUXjehZY/0Y3DPUw7PG0KhPLhIhIAiEEohPuInRvAjJy1ZjQyxmf+LqivqXhj6dVmVgmREQVlJGvxdv/OoejSRno0LQ+Nk7oDs8WdnLHkgXLhIionIq1OnxzOhVh0ZkwNjZGyPAOmNSnFUxr8QH252GZEBGVw2837mHOj/FITH+A7k0tsCqoD5o3rP0H2J+HZUJEVAYPCjX4Muoy/v1LGprYWuLr8V5oXHiLRfL/WCZERKUQQuDAxbtYuC8BmXlqTOrTCh8PeXQFe1zcbbnjGQyWCRHRM9zMfoh5kfE4flkFD6f6+MfE7ujcvG4eYH8elgkR0V9otDr84+R1rD6aDBMjI8z374ig3s51+gD787BMiIj+JC7tHub+eBFJd3Ph694EC152R9MGVnLHMngsEyIiADkPNfgiOgnbf7mBZg1q3xDxVY1lQkR1mhACe/+4g8X7E5Gdr8bbfVtj5mAX2NSyIeKrGl8tIqqz0rLyERIRj5MpmejSvAG+fbMHPJwayB2rRmKZEFGdoy7W4uvj1/DV8SswNzHGwpfdMb6XM0xq8D3Y5SZLmURGRmLTpk0AAG9vb8yePRuJiYmYO3cu8vPz0b17dyxcuBCmpqa4c+cOZs2ahaysLLRu3RphYWGwsbGRIzYR1QKnr2RiXkQ8rmXmw79zU8zz74gmteAe7HKr9vPcCgoKsHTpUmzduhWRkZE4d+4cYmJiMGvWLMyfPx/R0dEQQmDXrl0AgIULF2LcuHGIioqCh4cH1q9fX92RiagWyHhQiOn/+Q1v/OMX6ITAd2/1xLpx3VgklaTay0Sr1UKn06GgoADFxcUoLi6GqakpCgsL4enpCQAIDAxEVFQUNBoNzp49C19f3xLziYjKSqsT+FdMKgb+/QSiEu5ixqD2iJrhDW8XB7mj1SrV/jVXvXr18NFHH2Ho0KGwsrJCjx49YGZmBgeH//1hHRwcoFQqce/ePdSrVw+mpqYl5hMRlcUfN+8jJCIeF2/noF/7xlg0wgOtG/Nr8qpQ7WWSlJSEH374AT/99BNsbW3xySef4PTp0yXuRiaEgJGRkf6/f1beu5bFx8dXSu6qEhcXJ3eEUjGfNMwnXUUy5hfpsD0+D9FXH8LO0hh/69UAfZqbIDstCdlp8uerTtWVr9rL5NSpU+jduzcaNWoE4NFXV1u2bIFKpdKvk5mZCYVCAXt7e+Tm5kKr1cLExAQqlQoKhaJc2/Pw8ICFhUWl/g6VJS4uDl5eXnLHeCbmk4b5pCtvRiEEIn+/gyVHH10zMunFVvjbYBfYVtFdDw39NaxIPrVaXaEP4dV+zMTNzQ0xMTF4+PAhhBA4duwYevbsCQsLC32DRkZGwtvbG2ZmZujevTsOHDgAAIiIiIC3t3d1RyaiGuBKRh7e+McvmLHzdzg1tMLeD/siNMC9yoqESqr2PZO+ffvi0qVLCAwMhJmZGTp16oSpU6di8ODBCAkJQV5eHtzd3REUFAQACA0NRXBwMDZs2ICmTZti5cqV1R2ZiAxYoUaLdceuYOPPV2FlZoKlr3hgTI+WvGakmslyncnUqVMxderUEvPc3Nywe/fuJ9Z1cnLC1q1bqysaEdUgPyVlYP7eeNzMLkBgNyfMGdYBjesZ5tfatR2vgCeiGufO/QIs2ncJUQl30U5RD/+Z0gu92zaSO1adxjIhohpDo9Xh29OpCD+SDJ0QmOXriin92sDclPcZkRvLhIhqhHOp2QiJiEfS3VwMdFNgwcvuaGHP+68bCpYJERm0e/lFWH4wCTvP3USzBpbYOMELQzo2Kfc1Z1S1WCZEZJB0OoGj1x/i7f8eR25hMd7xboPpA9vzPiMGin8VIjI4iekPMC8iHufSHqBHq4ZYMrITXB1t5Y5FpWCZEJHBeFCowarDKfhXbCoaWJnhg+718XFgbxjzmhGDxzIhItk9HgZl6YFEZOapMa5nS8zydcXVxIsskhqCZUJEskpW5mJeRDx+uZ6NLs0bYMvE7ujc3E7uWFROLBMikkWeuhirjyTjm9OpqGdpis9f6YQxPVpwT6SGYpkQUbUSQmD/hXQs+e8lKB+oMbZnC8zydYO9jbnc0UgClgkRVZsrGXkI3RuP01ey4OFUH1+P90LXlg3ljkWVgGVCRFUuX12MtceuYMupa7AyM8HiEe4Y94IzR/atRVgmRFRlhBCIir+LRfsvIT2nEKO8miN4qBtH9q2FWCZEVCWuqfIQujcBJ1My0aFpfawd2xXdW9nLHYuqCMuEiCpVQZEWX/10BZt+vgYLU2MsCOiI8b2cYWrCkX1rM5YJEVUKIQQOXVJi0b5LuH2/AIFdnRA8zA0KW0u5o1E1YJkQkWRpWflYsDcBP11WwbWJLXZO7YUX2vBmVXUJy4SIKqxQo8X641fx9YmrMDM2QsjwDpjYpxXM+JVWncMyIaIKOZqoxIJ9CbiZXYCXuzTD3OEd0KQ+v9Kqq8pVJkVFRTA2NoapKTuIqK66mf0QC/cl4EhiBtop6mH7lBfQp21juWORzJ7bCllZWdi0aRMOHz6M9PR0GBkZoXnz5vDz88OkSZNgb89T/YjqgkKNFhv+/ystE2MjfDbUDW++2Jr3XycAzymTiIgIfPfdd/Dz88O6devQvHlzmJmZ4ebNmzh58iTefPNNTJo0Ca+88kp15SWiaiaEQHSCEkv+ewm37hXAv3NTzB3eAU0bWMkdjQxIqWWSk5OD3bt3w9i45CcPFxcXuLi4YNKkSdi6dWuVBiQi+VxV5WHB/1946NrEFv+Z0gu92/IsLXpSqWUyceJE/c83b95EixYtkJeXhxs3bqBjx44wMTHBpEmTqjojEVWzPHUx1h5NwT9PX4elqQlCAzpiAi88pFKU6V/G1q1b8f777wMA7t27h2nTpuH777+v0mBEVP0e3fHwNl4KO46NP1/DSE8nHPukP958sTWLhEpVptOydu7ciR07dgAAWrRogYiICIwbNw6vvfZalYYjouqTmP4AoZEJ+DU1G52cGuDrCV7oxuHhqYzK9FFDq9WiXr16+mlbW1sYGVV86Ohjx44hMDAQQ4cOxZIlSwAAMTExCAgIwJAhQxAeHq5fNzExEYGBgfD19cXcuXNRXFxc4e0S0ZNyHmoQGhmP4WtOIiUjF8sCOyHigxdZJFQuZSqTNm3aICwsDDdv3sTNmzexevVqtGrVqkIbvHnzJkJDQ7F+/Xrs3bsXly5dwokTJzBnzhysX78eBw4cQHx8PE6cOAEAmDVrFubPn4/o6GgIIbBr164KbZeIStIJgR2/3sCAvx/H1jNpGN/LGT990h9je7bkfUao3MpUJgsXLkRqaipGjhyJUaNGITU1FQsWLKjQBg8fPoxhw4bB0dERZmZmCA8Ph5WVFZydndGiRQuYmpoiICAAUVFRuH37NgoLC+Hp6QkACAwMRFRUVIW2S0T/8/vN+/jsaDaC91xEWwcb7J/WD4tGeMDOmrfOpYop0zGTxo0bY926dZWywbS0NJiZmeHdd99Feno6+vfvj/bt28PBwUG/jkKhgFKpREZGRon5Dg4OUCqVlZKDqC7KzFNjRVQSdp27hYaWxlg12hMjPJtJ+tqaCChjmahUKixZsgSnTp2CiYkJXnrpJXz22Wdo0KBBuTeo1Wpx7tw5bN26FdbW1njvvfdgaWlZ4h+zEAJGRkbQ6XRPnV8e8fHx5c5YneLi4uSOUCrmk8ZQ8ml1AlFXH2JHQh7UxQIjXG3wWgcbWOnu4vz5u3LHK5WhvIbPwnyPlKlMgoOD0bFjR0RERECr1WLnzp2YN28e1qxZU+4NNm7cGL1799YPwzJo0CBERUXBxMREv45KpYJCoYCjoyNUKpV+fmZmJhQKRbm25+HhAQsLw7xFaFxcHLy8vOSO8UzMJ42h5DtzLQsL9iYg6W4u+rVvjNAAd7RT1DOYfKUx9Iy1MZ9ara7Qh/AyHTO5e/cuPv74Y7Ro0QKtWrXC7NmzceXKlXJvDAAGDBiAU6dO4cGDB9BqtTh58iT8/Pxw/fp1pKWlQavVYv/+/fD29oaTkxMsLCz0zRoZGQlvb+8KbZeorknPKcC0//yGMZvOILewGF+P98J3b/VEO0W95z+YqJzKtGfSrFkz3LhxAy1btgQAZGRklHsP4bEuXbrg7bffxrhx46DRaPDiiy9i7NixaNOmDaZNmwa1Wg0fHx/4+fkBAMLCwhASEoK8vDy4u7sjKCioQtslqivUxVpsOXUd645dQbFO4KOB7fGuT1tYmZs8/8FEFVSmMjE2NsbIkSPRt29fmJiYIDY2Fo6Ojnj33XcBAF9//XW5Njpq1CiMGjWqxLzevXtj7969T6zr5uaG3bt3l+v5ieqq45czsHDfJVzPzMfgjk0wb3hHtGxkLXcsqgPKVCZ+fn76PQUA/KqJyMCkZuZjyX8v4UhiBlo3tsG3b/ZAf9eKfXtAVBFlKhMOMU9kmPLUxVh37Ar+eeo6zEyMMNvPDW/1bQULU36lRdWr1DLp2rVrqafinj9/vtIDEdHz6XQCP/52G19EJSEjV43Abk4I9nODgrfNJZmUWiYeHh5ITU1FQEAAAgICUL9+/erKRUTP8PvN+1iwNwG/37yPLi3ssHGCF7pyHC2SWallsnXrVty5cwcRERH49NNP0aZNGwQGBqJfv35P3DCLiKpWRm4hVkRdxu64W3CwtUDYa10Q2NUJxhxHiwzAc4+ZNGvWDO+//z7ef/99nD9/HhEREfjiiy8wYMAAzJo1qzoyEtVpRcU6fHP6OtYeuwJ1sRbv+LTBhwPawdbSTO5oRHrl2r1o1aoV2rVrBwsLCxw9erSqMhHR/zuWpITvqp+x7GASXmhtj0MzffDZ0A4sEjI4z90zUavVOHLkCCIiIpCQkABfX1+EhobqR/Ilosp3VZWHxfsv4fhlFdo0tsE3b/bAAJ7qSwas1DL57LPPcOzYMXTv3h2vv/46+vfvDzMzfiIiqioPCjVYezQF35xOhZWZCeYO64CJfVrB3JTHKMmwlVomP/74IxwcHHDjxg2sWbPmiYEd9+3bV6XhiOoKnU5gd9wtrIhOQlZ+EV73aoFPfF3hYGuYg5QS/VWpZfLdd99VVw6iOisu7R4W7kvAhVs56NbSDv+c1AOdm9vJHYuoXEotE3t7e7Rr167UJ0hJSUH79u0rNRRAUfI8AAAgAElEQVRRXaB8UIjlB5Pw42+30aS+BW9URTVaqWWyadMm1K9fH2PGjHmiVK5du4atW7fiwYMH+Pvf/16lIYlqk0LNo1F9v/rpCoq1Ah8MaIv3+7eDjUWZRjciMkil/utdsWIFDh48iOnTp0OtVsPZ2Rk6nQ43btyApaUlPvjgAwwfPry6shLVaEIIHL6kxJL/JuJG9kMM7tgEIcM7wLmRjdzRiCR77kehoUOHYujQoUhJScG1a9dgZGSE1q1b86stonK4fDcXS/57CSdTMtFOUQ9bJ/dEv/YOcsciqjRl3q9u3759iQL5z3/+g7Fjx1ZJKKLaIitPjfAjydj+yw3YWpphvn9HTOjtDDMTnupLtUuFv6T98ssvWSZEz1BUrMPe5Hzs2XccD4u0mNDLGTMGuaChjbnc0YiqRIXLRAhRmTmIagUhBI4mZmDpgURcz8yHt4sD5g3vgPZNbOWORlSlKlwmPH2RqKQ/Hxdp62CDuX0bYop/T7ljEVWLUsskISHhqfOFENwzIfp/2flFWHn4sv64SGhAR4zv5YwLv/8mdzSialNqmUybNu2Zyxo25M14qG4rKtbhu9hUrD6awuMiVOeVWibHjh177hPs378f/v7+lRaIyNDxuAjRkyRfcrtlyxaWCdUZfz0uwqHhiR6RXCY8dkJ1wbOOi/B6EaJHJJcJz+qi2ozHRYjKhiPLET0Fj4sQlQ/LhOgveFyEqPxkO2byxRdf4N69e1i+fDkSExMxd+5c5Ofno3v37li4cCFMTU1x584dzJo1C1lZWWjdujXCwsJgY8MRVqlq8LgIUcWVWiYTJkwo9ZjId999h4CAgHJvNDY2Fj/++CP69+8PAJg1axaWLFkCT09PzJkzB7t27cK4ceOwcOFCjBs3DsOHD8dXX32F9evXY9asWeXeHlFp1MVafHs6Fet+usLjIkQVVOpHrvHjx+ONN96AQqGAtbU1JkyYgEmTJqFhw4ZwdnYGAEyePLlcG7x//z7Cw8Px7rvvAgBu376NwsJCeHp6AgACAwMRFRUFjUaDs2fPwtfXt8R8osoihMB/L6Rj0MoTWHYwCT1a2SN6Rj8sHOHBIiEqp1L3TB6/kW/ZsgU7duyAsfGj7unfvz9Gjx5doQ3Onz8fM2fORHp6OgAgIyMDDg7/u6+Dg4MDlEol7t27h3r16sHU1LTEfKLK8NuNe1jy30TEpd2Dm6Mt/j35BfRt31juWEQ1VpmOmdy7dw9qtRpWVlYAgPz8fOTk5JR7Y99//z2aNm2K3r17Y8+ePQAAnU5X4qs0IQSMjIz0//2zipyGHB8fX+7HVKe4uDi5I5SqtuXLyNdi28VcnLpZCDtLY7zXvT4GtLKCyYM0xMWlyZ6vuhl6PsDwMzLfI2UqE39/f7z++usYPHgwhBCIiorC66+/Xu6NHThwACqVCiNGjEBOTg4ePnwIIyMjqFQq/TqZmZlQKBSwt7dHbm4utFotTExMoFKpoFCU/4waDw8PWFhYlPtx1SEuLg5eXl5yx3im2pQvt1CD9cevYsup6zA2Aqa/1A7v+LSt0vuu16bXTy6GnrE25lOr1RX6EF6m/5M++ugjuLu748yZMwCA4OBg+Pj4lHtj33zzjf7nPXv24Ndff8WyZcvg7++v/6UjIyPh7e0NMzMzdO/eHQcOHEBAQAAiIiLg7e1d7m1S3Vas1WHnuZtYeSgZWflFCOzqhE98XdHMzkruaES1Spk/lg0aNAiDBg2qkhBhYWEICQlBXl4e3N3dERQUBAAIDQ1FcHAwNmzYgKZNm2LlypVVsn2qnY5fzsDS/yYiJSMPPVvb45vhHdC5uZ3csYhqJdkuWgwMDERgYCAAwM3NDbt3735iHScnJ2zdurW6o1ENl3T3AZb+NxEnUzLh3MgaX4/3gq97Ew79Q1SFeAU81RqqXDVWHk7GzrM3UM/CFPP8O2JCL2eYm/KiQ6KqxjKhGq9Qo8WWU9ex/qcrUBfrMLFPK0x/qT2vFSGqRiwTqrF0QiDy99v44mAS7uQUYkjHJgge6oY2DvXkjkZU57BMqEY6m5qNz45m48o9JTyc6mPlaE/0atNI7lhEdRbLhGqUtKx8fBGVhAMX78Leyhh/f60LXunqBGNjHlwnkhPLhGqE+w+LsO7YFfwrNhWmxsb422AXdK+Xgz5ezeWORkRgmZCBK9RosTU2DWuPpSBXXYzXvJrj4yGuaFLf0uCHsSCqS1gmZJB0OoF9F+7gy+jLuHWvAD4uDvhsmBvcHOvLHY2InoJlQgYn9moWlh1MxIVbOejYtD7+PbkzR/QlMnAsEzIYKcpcLD+YhKNJGWjWwJIH14lqEJYJyS4jtxDhh1Ow8+wN2JibYrafG958sRUszUzkjkZEZcQyIdnkq4ux6edr2HzyGoqKdQjq3QrTB7aHPa9cJ6pxWCZU7Yq1Ouw6dwvhR5KhylVjWCdHfOrrhlaNbeSORkQVxDKhaiOEwLGkDCw/mISUjDx4OTfE1+O94OXcUO5oRCQRy4SqxcVbOVh64BLOXMtG68Y2+Hp8N/i6O3JYeKJagmVCVepm9kOEHbqMyN/vwN7GHAtfdse4F1rCzITDwhPVJiwTqhI5DzX46vgVfHs6FUZGwPv92+Ld/m1R39JM7mhEVAVYJlSp1MWPhz+5ggeFGgR2bY6Ph7jwnutEtRzLhCqFTiew/2I6voxOws3sAvRr3xifDe2Ajs04/AlRXcAyIclOpWRieVQi4m8/gJujLf71Vk/4uDjIHYuIqhHLhCos/nYOvohKwsmUTDjZWWHl610wwtMJJhz+hKjOYZlQud3Ieoi/H350hpadtRlChnfA+F7OHP6EqA5jmVCZZeWpsfbYFWz7JQ0mxkZ4v39bvOPTFg2seIYWUV3HMqHnelhUjC0nr2Pjz9fwsKgYo3u0wEcDXeDYwFLuaERkIFgm9EzFOoF/n0nDqiMpyMxTY0jHJvjUzxXtFLZyRyMiA8MyoScIIXAw/i4WR2ciPU+J7s4NsXFCN3g528sdjYgMFMuESjhzLQvLDibhj5v30by+Kf4R1B0DOyg4hhYRlYplQgCAxPQHWBGVhJ8uq9C0gSVWjOqMVlCiZ8cmckcjohpAljJZt24dDh48CADw8fHBp59+ipiYGCxbtgxqtRpDhw7FzJkzAQCJiYmYO3cu8vPz0b17dyxcuBCmpuzAynLr3kOsPJyMH3+7DVsLUwQPdcOkPo/uchgXlyF3PCKqIap96NaYmBicOnUKP/74IyIiIpCQkID9+/djzpw5WL9+PQ4cOID4+HicOHECADBr1izMnz8f0dHREEJg165d1R25VrqXX4Ql+y/hpbAT2H8hHVP7tcHPnw7Auz5teb0IEZVbtX/Ed3BwQHBwMMzNH92atW3btkhNTYWzszNatGgBAAgICEBUVBTatWuHwsJCeHp6AgACAwOxZs0ajBs3rrpj1xoFRVp8E3MdG45fRZ66GKO6NcfMwRyIkYikqfYyad++vf7n1NRUHDx4EOPHj4eDw//GclIoFFAqlcjIyCgx38HBAUqlslrz1hbFWh2+j7uFVUeSoXygxqAOCszydYOrI0/zJSLpZDv4kJKSgnfeeQeffvopTExMkJqaql8mhICRkRF0Ol2Js4gezy+P+Pj4yopcJeLi4qr0+XVCIPZWIXbE5+FOnhYu9mb4oL893B2MkXc7GXG35c0nFfNJY+j5AMPPyHyPyFImcXFxmD59OubMmYPhw4fj119/hUql0i9XqVRQKBRwdHQsMT8zMxMKhaJc2/Lw8ICFhUWlZa9McXFx8PLyqpLnFkLgRLIKX0ZfRsKdB3BpUg+bXnHF4I5NylzIVZmvMjCfNIaeDzD8jLUxn1qtrtCH8Govk/T0dHzwwQcIDw9H7969AQBdunTB9evXkZaWhubNm2P//v149dVX4eTkBAsLC/0LEhkZCW9v7+qOXOPEpd3Diqgk/HI9G80bcjRfIqp61V4mW7ZsgVqtxvLly/XzxowZg+XLl2PatGlQq9Xw8fGBn58fACAsLAwhISHIy8uDu7s7goKCqjtyjXH5bi6+jL6MI4lKNK5ngYUvu2Nsz5YwN+X91omoalV7mYSEhCAkJOSpy/bu3fvEPDc3N+zevbuqY9VoN7IeIvxIMiJ+v416FqaY5euKN19sBWtzXo9DRNWD7zY1WEZuIdYdu4L//HoDxkZGmOrdBu/5tIWdtbnc0YiojmGZ1EA5BRpsPHEV35xORZFWh9E9WmD6S+05JDwRyYZlUoMUFGnxbUwqvj5xFTkFGrzcpRn+NtgFrRrbyB2NiOo4lkkNoNHqsPPsTaw5moKMXDUGuDrgE19XuDdrIHc0IiIALBODptMJ7LtwBysPJyMt6yF6tGqIdeO6oWdr3leEiAwLy8QACSHw0+UMfBmdjMT0B3BztMU3k3qgv6sD7ytCRAaJZWJgfr2ejS+jk3A29R6cG1lj9RhPBHRuBmNecEhEBoxlYiDib+fg74cu46fLKihsLbBkpAdG92gBMxNecEhEho9lIrNkZS7CDyfjYPxdNLAyw2y/RzensjLnPUWIqOZgmcgkNTMfq3+5j5O7f4aNuSk+Gtgek/u1Rn1LM7mjERGVG8ukmt2+X4B1x1Kw69wtmBgJTPVug3e926KhDa9aJ6Kai2VSTTJyC7H+p6vY/ssNAMCEXs7o2ygfg17sIHMyIiLpWCZV7F5+ETb+fA3/ink09MlrXs0xbWB7ONlZGfxNdYiIyoplUkVyCzXYcuo6tpy8jryiYozo0gwzBnHoEyKqnVgmlexhUTG+i03D1yeu4v5DDfzcHTFzsAvvtU5EtRrLpJKoi7X4zy83sO6nq8jMU6O/qwM+HuyKTs05fhYR1X4sE4k0Wh1+iLuFNUdTcCenEC+0tseG8d3QoxXHzyKiuoNlUkFancC+P+5g1ZFkpGY9RJcWdlgxqgtebNeI42cRUZ3DMiknIQSiE+5i5eFkJCvz4OZoi38EdcfADgqWCBHVWSyTMhJC4EhiBlYdSUbCnQdo42CDdeO6YphHUw7CSER1HsvkOYQQOH5ZhfAjybhwKwfOjawR9loXjPRsBlMOwkhEBIBl8kxCCJxMycTKw8n4/eZ9NG9ohRWjOuOVrk4cyZeI6C9YJn8hhEDM1SyEH07GubR7cLKzwrLATni1W3OYm7JEiIiehmXyJ2euZWHl4WT8ej0bjvUtsXikB17v3hwWphwOnoioNCwTAOdSs7HycDJirmZBYWuBBQEdMaZnS1iasUSIiMqiTpfJ+Rv3EH44GSdTMtG4njnm+XfEGy+wRIiIyqtOlskfN+8j/Egyjl9Wwd7GHHOGuWF8L2dYm9fJl4OISLIa8e65b98+bNiwAcXFxZg4cSLeeOONCj1P/O0crDqSjCOJGbCzfnSL3KDezrCxqBEvAxGRwTL4d1GlUonw8HDs2bMH5ubmGDNmDF544QW0a9euzM+RmP4Aq44kIzpBifqWpvhkiAsm9mkFW94il4ioUhh8mcTExKBXr16ws7MDAPj6+iIqKgoffvhhmR4/PzIeO8/fha2FKWYMao+3+vI+60RElc3gyyQjIwMODg76aYVCgQsXLjz3cUIIAMDl2/cwa1AbjO7eErZWpgB0UKvVVRW33Awpy9MwnzTMJ52hZ6xt+YqKigD87z20rAy+THQ6XYkBFIUQZRpQUaPRAADmeTcE8BBpV5OqKqIk8fHxckcoFfNJw3zSGXrG2ppPo9HA0tKyzOsbfJk4Ojri3Llz+mmVSgWFQvHcx9nY2MDFxQVmZmYczZeIqIyEENBoNLCxKd8txg2+TPr06YO1a9ciOzsbVlZWOHToEBYvXvzcxxkbG8PWlrfKJSIqr/LskTxm8GXSpEkTzJw5E0FBQdBoNBg1ahQ6d+4sdywiIvoTI1HeoyxERER/wWFwiYhIMpYJERFJxjIhIiLJWCZERCQZy4SIiCRjmRARkWQsEyIikoxlQkREkrFMiIhIMpYJERFJxjIhIiLJWCZERCQZy4SIiCRjmRARkWQsEyIikoxlQkREkrFMiIhIMpYJERFJxjIhIiLJWCZERCQZy4SIiCRjmRARkWQsEyIikoxlQkREkrFMiIhIMpYJERFJxjIhIiLJWCZERCQZy4SIiCQziDLJy8uDv78/bt269cSyxMREBAYGwtfXF3PnzkVxcbEMCYmIqDSyl8kff/yBsWPHIjU19anLZ82ahfnz5yM6OhpCCOzatat6AxIR0XPJXia7du1CaGgoFArFE8tu376NwsJCeHp6AgACAwMRFRVV3RGJiOg5TOUOsHTp0mcuy8jIgIODg37awcEBSqWyTM+r0+mQn58PMzMzGBkZSc5JRFQXCCGg0WhgY2MDY+Oy72/IXial0el0JYpACFHmYsjPz0dycnJVRSMiqtVcXFxga2tb5vUNukwcHR2hUqn005mZmU/9OuxpzMzMADx6QczNzaskn1Tx8fHw8PCQO8YzMZ80zCedoWesjfmKioqQnJysfw8tK4MuEycnJ1hYWCAuLg5eXl6IjIyEt7d3mR77eA/G3NwcFhYWVRlTEkPOBjCfVMwnnaFnrK35ynt4QPYD8E8zZcoUXLx4EQAQFhaGZcuWwc/PDw8fPkRQUJDM6YiI6K8MZs/k2LFj+p83b96s/9nNzQ27d++WIxIREZWRQe6ZEBFRzcIyISIiyVgmREQkGcuEiIgkY5kQEZFkLBMiIpKMZUJERJKxTIiISDKWCRERScYyISIiyVgmREQkGcuEiIgkY5kQEZFkLBMiIpKMZUJERJKxTIiISDKWCRERScYyISIiyVgmREQkGcuEiIgkY5kQEZFkLBMiIpKMZUJERJKxTIiISDKWCRERScYyISIiyVgmREQkGcuEiIgkY5kQEZFkLBMiIpJM9jLZt28fhg0bhiFDhmDbtm1PLE9ISMCrr76Kl19+Ge+88w4ePHggQ0oiIiqNrGWiVCoRHh6O7du3IyIiAjt37sSVK1dKrLN06VJMnz4de/fuRevWrbFlyxaZ0hIR0bPIWiYxMTHo1asX7OzsYG1tDV9fX0RFRZVYR6fTIT8/HwBQUFAAS0tLOaISEVEpZC2TjIwMODg46KcVCgWUSmWJdYKDgxESEoK+ffsiJiYGY8aMqe6YRET0HKZyblyn08HIyEg/LYQoMV1YWIi5c+fi22+/RefOnfHNN99g9uzZ2LRpU5m3ER8fX6mZK1tcXJzcEUrFfNIwn3SGnpH5HpG1TBwdHXHu3Dn9tEqlgkKh0E8nJyfDwsICnTt3BgCMHj0aq1evLtc2PDw8YGFhUTmBK1lcXBy8vLzkjvFMzCcN80ln6BlrYz61Wl2hD+Gyfs3Vp08fxMbGIjs7GwUFBTh06BC8vb31y52dnXH37l1cu3YNAHD06FF06tRJrrhERPQMsu6ZNGnSBDNnzkRQUBA0Gg1GjRqFzp07Y8qUKZg+fTo6deqEZcuWYcaMGRBCoFGjRvj888/ljExERE8ha5kAQEBAAAICAkrM27x5s/5nHx8f+Pj4VHcsIiIqB9kvWiQiopqPZUJERJKxTIiISDKWCRERScYyISIiyVgmREQkGcuEiIgkY5kQEZFkLBMiIpKMZUJERJKxTIiISDKWCRERScYyISIiyVgmREQkGcuEiIgkY5kQEZFkLBMiIpKMZUJERJKxTIiISDKWCRERScYyISIiyVgmREQkGcuEiIgkY5kQEZFkLBMiIpKMZUJERJKxTIiISDKWCRERScYyISIiyWQvk3379mHYsGEYMmQItm3b9sTya9euYcKECXj55ZcxefJk5OTkyJCSiIhKI2uZKJVKhIeHY/v27YiIiMDOnTtx5coV/XIhBN577z1MmTIFe/fuRYcOHbBp0yYZExMR0dPIWiYxMTHo1asX7OzsYG1tDV9fX0RFRemXJyQkwNraGt7e3gCAd999F2+88YZccYmI6BlkLZOMjAw4ODjopxUKBZRKpX76xo0baNy4MebMmYNXXnkFoaGhsLa2liMqERGVwlTOjet0OhgZGemnhRAlpouLi/Hrr7/i3//+Nzp16oRVq1Zh+fLlWL58eZm3ER8fX6mZK1tcXJzcEUrFfNIwn3SGnpH5HpG1TBwdHXHu3Dn9tEqlgkKh0E87ODjA2dkZnTp1AgD4+/tj+vTp5dqGh4cHLCwsKidwJYuLi4OXl5fcMZ6J+aRhPukMPWNtzKdWqyv0IVzWr7n69OmD2NhYZGdno6CgAIcOHdIfHwGArl27Ijs7G0lJSQCAY8eOwd3dXa64RET0DLLumTRp0gQzZ85EUFAQNBoNRo0ahc6dO2PKlCmYPn06OnXqhK+++gohISEoKCiAo6MjVqxYIWdkIiJ6ClnLBAACAgIQEBBQYt7mzZv1P3fp0gW7d++u7lhERFQOsl+0SERENR/LhIiIJGOZEBGRZCwTIiKSjGVCRESSsUyIiEgylgkREUnGMiEiIslYJkREJBnLhIiIJGOZEBGRZCwTIiKSjGVCRESSsUyIiEgylgkREUnGMiEiIslYJkREJBnLhIiIJGOZEBGRZCwTIiKSjGVCRESSsUyIiEgylgkREUnGMiEiIslYJkREJBnLhIiIJGOZEBGRZCwTIiKSjGVCRESSsUyIiEgy2ctk3759GDZsGIYMGYJt27Y9c73jx4/jpZdeqsZkRERUVqZyblypVCI8PBx79uyBubk5xowZgxdeeAHt2rUrsV5mZia++OILmVISEdHzyLpnEhMTg169esHOzg7W1tbw9fVFVFTUE+uFhITgww8/lCEhERGVhax7JhkZGXBwcNBPKxQKXLhwocQ63333HTp27IguXbpUaBvx8fGSMla1uLg4uSOUivmkYT7pDD0j8z0ia5nodDoYGRnpp4UQJaaTk5Nx6NAhfPvtt7h7926FtuHh4QELCwvJWatCXFwcvLy85I7xTMwnDfNJZ+gZa2M+tVpdoQ/hsn7N5ejoCJVKpZ9WqVRQKBT66aioKKhUKrz66quYOnUqMjIyMG7cODmiEhFRKWQtkz59+iA2NhbZ2dkoKCjAoUOH4O3trV8+ffp0REdHIzIyEps2bYJCocD27dtlTExERE8ja5k0adIEM2fORFBQEEaOHAl/f3907twZU6ZMwcWLF+WMRkRE5SDrMRMACAgIQEBAQIl5mzdvfmK95s2b49ixY9UVi4iIykH2ixaJiKjmY5kQEZFkLBMiIpKMZUJERJKxTIiISDKWCRERScYyISIiyVgmREQkGcuEiIgkY5kQEZFkLBMiIpKMZUJERJKxTIiISDKWCRERScYyISIiyVgmREQkGcuEiIgkY5kQEZFkLBMiIpKMZUJERJKxTIiISDKWCRERScYyISIiyVgmREQkGcuEiIgkY5kQEZFkLBMiIpKMZUJERJKxTIiISDLZy2Tfvn0YNmwYhgwZgm3btj2x/MiRIxgxYgRefvllvP/++8jJyZEhJRERlUbWMlEqlQgPD8f27dsRERGBnTt34sqVK/rleXl5WLBgATZt2oS9e/fC1dUVa9eulTExERE9jaxlEhMTg169esHOzg7W1tbw9fVFVFSUfrlGo0FoaCiaNGkCAHB1dUV6erpccYmI6BlkLZOMjAw4ODjopxUKBZRKpX66YcOGGDx4MACgsLAQmzZtwqBBg6o9JxERlc5Uzo3rdDoYGRnpp4UQJaYfy83NxQcffAA3Nze88sor5dpGfHy85JxVKS4uTu4IpWI+aZhPOkPPyHyPyFomjo6OOHfunH5apVJBoVCUWCcjIwOTJ09Gr169MGfOnHJvw8PDAxYWFpKzVoW4uDh4eXnJHeOZmE8a5pPO0DPWxnxqtbpCH8Jl/ZqrT58+iI2NRXZ2NgoKCnDo0CF4e3vrl2u1Wrz77rsYOnQo5s6d+9S9FiIikp+seyZNmjTBzJkzERQUBI1Gg1GjRqFz586YMmUKpk+fjrt37+LSpUvQarWIjo4G8GhPY+nSpXLGJiKiv5C1TAAgICAAAQEBJeZt3rwZANCpUyckJSXJEYuIiMpB9osWiYio5mOZEBGRZCwTIiKSjGVCRESSsUyIiEgylgkREUnGMiEiIslYJkREJBnLhIiIJGOZEBGRZCwTIiKSjGVCRESSsUyIiEgylgkREUnGMiEiIslYJkREJBnLhIiIJGOZEBGRZCwTIiKSjGVCRESSsUyIiEgylgkREUnGMiEiIslYJkREJBnLhIiIJGOZEBGRZCwTIiKSjGVCRESSsUyIiEgy2ctk3759GDZsGIYMGYJt27Y9sTwxMRGBgYHw9fXF3LlzUVxcLENKIiIqjaxlolQqER4eju3btyMiIgI7d+7ElStXSqwza9YszJ8/H9HR0RBCYNeuXTKlJSKiZ5G1TGJiYtCrVy/Y2dnB2toavr6+iIqK0i+/ffs2CgsL4enpCQAIDAwssZyIiAyDqZwbz8jIgIODg35aoVDgwoULz1zu4OAApVJZpucWQgAAioqKKilt1VCr1XJHKBXzScN80hl6xtqW7/F75uP30LKStUx0Oh2MjIz000KIEtPPW14ajUYDAEhOTq6ktFUjPj5e7gilYj5pmE86Q89YW/NpNBpYWlqWeX1Zy8TR0RHnzp3TT6tUKigUihLLVSqVfjozM7PE8tLY2NjAxcUFZmZmZS4gIqK6TggBjUYDGxubcj1O1jLp06cP1q5di+zsbFhZWeHQoUNYvHixfrmTkxMsLCwQFxcHLy8vREZGwtvbu0zPbWxsDFtb26qKTkRUa5Vnj+QxI1HeL8Yq2b59+7Bx40ZoNBqMGjUKU6ZMwZQpUzB9+nR06tQJSUlJCAkJQV5eHtzd3bFs2TKYm5vLGZmIiP5C9jIhIqKaT/aLFomIqOZjmRARkWQsEyIikoxlQkREkrFMiIhIMpYJERFJxjIhIiLJWCZERCQZy4SIiCRjmRARkWQsEyIikoxlQkREkrFMiIhIMpYJERFJxq2k+DAAABQCSURBVDIhIiLJWCZERCQZy4SIiCRjmRARkWQsEyIikoxlQkREkrFMiIhIMpYJERFJxjIhIiLJWCZERCQZy4SIiCRjmRARkWQsEyIikoxlQkREkrFMiIhIMpYJERFJZhBlkvd/7d1/TNT34cfxF1agO2tCu3GwdIlbsohrAdtgUsIaTDsFuXpQlcVWLUtWUOoaOrb+WMXUzcXZNUtwXZu20GZmKSS6tiIsDTBnmnWBpOWyteI0jDTLVlvucDhFBnjA+/uH4fPdqXeCb7zP4Z6Pv3h/3p+7zysf3uF1nzvv4/nzWrt2rT799NPL5k6cOKH169eruLhYdXV1mpiYcCEhACAW18vko48+0sMPP6y///3vV5x/6qmn9Nxzz6mjo0PGGB08eDC+AQEAV7XQ7QAHDx7Url279PTTT182d+rUKY2Njemuu+6SJK1fv14vvviiNm3adNXnnZqa0sjIiJKTk5WUlDTnuQHgRmSMUTgc1qJFi7RgwcyvN1wvkz179kSdC4VCSk9Pd8bp6ekKBoMzet6RkRH19fVZ5wOA/0VLly7V4sWLZ7y/62USy9TUVMRVhTFmxlcZycnJki6ekJSUlOuSz1Zvb6+ys7PdjhEV+eyQz16iZ7wR8124cEF9fX3O39CZSugyyczM1ODgoDM+ffq0vF7vjB47XTopKSlKTU29LvnmQiJnk8hni3z2Ej3jjZpvth8PuP4BfCy33367UlNTFQgEJEmHDx9WYWGhy6kAAJdKyDKpqqrSsWPHJEm/+MUvtHfvXq1Zs0b/+c9/VFFR4XI6AMClEuZtrqNHjzo/NzY2Oj8vW7ZMb731lhuRAAAzlJBXJgCA+YUyAQBYo0wAANYoEwCANcoEAGCNMgEAWKNMAADWKBMAgDXKBABgjTIBAFijTAAA1igTAIA1ygQAYI0yAQBYo0wAANYoEwCANcoEAGCNMgEAWKNMAADWKBMAgDXKBABgjTIBAFijTAAA1igTAIA1ygQAYI0yAQBYo0wAANYoEwCANcoEAGDN9TJpa2uTz+dTUVGRmpqaLps/fvy4NmzYoNLSUm3btk3nzp1zISUAIBZXyyQYDKq+vl7Nzc1qaWnRgQMH1N/fH7HPnj17VFNTo9bWVn3ta1/TG2+84VJaAEA0rpZJV1eX8vPzlZaWJo/Ho+LiYrW3t0fsMzU1pZGREUnS6Oiobr75ZjeiAgBicLVMQqGQ0tPTnbHX61UwGIzY50c/+pF27type++9V11dXXrooYfiHRMAcBUL3Tz41NSUkpKSnLExJmI8Njamuro67d+/X7m5ufr1r3+tZ555Rg0NDTM+Rm9v75xmnmuBQMDtCDGRzw757CV6RvJd5GqZZGZmqqenxxkPDg7K6/U6476+PqWmpio3N1eStHHjRv3yl7+c1TGys7OVmpo6N4HnWCAQUF5entsxoiKfHfLZS/SMN2K+8fHxa3oR7urbXAUFBeru7tbQ0JBGR0fV2dmpwsJCZ37JkiUaGBjQJ598Ikn6wx/+oJycHLfiAgCicPXKJCMjQ7W1taqoqFA4HFZ5eblyc3NVVVWlmpoa5eTkaO/evfr+978vY4y++MUv6mc/+5mbkQEAV+BqmUiS3++X3++P2NbY2Oj8vHLlSq1cuTLesQAAs+D6lxYBAPMfZQIAsEaZAACsUSYAAGuUCQDAGmUCALBGmQAArFEmAABrlAkAwBplAgCwRpkAAKxRJgAAa5QJAMAaZQIAsEaZAACsUSYAAGuUCQDAGmUCALBGmQAArFEmAABrlAkAwBplAgCwRpkAAKxRJgAAa5QJAMAaZQIAsEaZAACsUSYAAGuUCQDAmutl0tbWJp/Pp6KiIjU1NV02/8knn+iRRx5RaWmpHn30UZ09e9aFlACAWFwtk2AwqPr6ejU3N6ulpUUHDhxQf3+/M2+M0WOPPaaqqiq1trbqG9/4hhoaGlxMDAC4ElfLpKurS/n5+UpLS5PH41FxcbHa29ud+ePHj8vj8aiwsFCSVF1drc2bN7sVFwAQhatlEgqFlJ6e7oy9Xq+CwaAz/sc//qEvfelL2rFjh9atW6ddu3bJ4/G4ERUAEMNCNw8+NTWlpKQkZ2yMiRhPTEzogw8+0JtvvqmcnBzt27dPzz//vJ5//vkZH6O3t3dOM8+1QCDgdoSYyGeHfPYSPSP5LnK1TDIzM9XT0+OMBwcH5fV6nXF6erqWLFminJwcSdLatWtVU1Mzq2NkZ2crNTV1bgLPsUAgoLy8PLdjREU+O+Szl+gZb8R84+Pj1/Qi3NW3uQoKCtTd3a2hoSGNjo6qs7PT+XxEku6++24NDQ3p5MmTkqSjR4/qzjvvdCsuACAKV69MMjIyVFtbq4qKCoXDYZWXlys3N1dVVVWqqalRTk6OXn75Ze3cuVOjo6PKzMzUCy+84GZkAMAVuFomkuT3++X3+yO2NTY2Oj8vX75cb731VrxjAQBmwfUvLQIA5j/KBABgjTIBAFijTAAA1igTAIA1ygQAYI0yAQBYo0wAANYoEwCANcoEAGCNMgEAWKNMAADWKBMAgDXKBABgjTIBAFijTAAA1igTAIA1ygQAYI0yAQBYo0wAANYoEwCANcoEAGCNMgEAWKNMAADWKBMAgDXKBABgjTIBAFijTAAA1igTAIA1ygQAYM31Mmlra5PP51NRUZGampqi7vfee+/p/vvvj2MyAMBMLXTz4MFgUPX19XrnnXeUkpKihx56SPfcc4++/vWvR+x3+vRp/fznP3cpJQDgaly9Munq6lJ+fr7S0tLk8XhUXFys9vb2y/bbuXOnHn/8cRcSAgBmwtUyCYVCSk9Pd8Zer1fBYDBin9/85je64447tHz58njHAwDMkKtvc01NTSkpKckZG2Mixn19fers7NT+/fs1MDBwTcfo7e21znk9BQIBtyPERD475LOX6BnJd5GrZZKZmamenh5nPDg4KK/X64zb29s1ODioDRs2KBwOKxQKadOmTWpubp7xMbKzs5WamjqnuedKIBBQXl6e2zGiIp8d8tlL9Iw3Yr7x8fFrehHu6ttcBQUF6u7u1tDQkEZHR9XZ2anCwkJnvqamRh0dHTp8+LAaGhrk9XpnVSQAgPhwtUwyMjJUW1uriooKPfjgg1q7dq1yc3NVVVWlY8eOuRkNADALrr7NJUl+v19+vz9iW2Nj42X7feUrX9HRo0fjFQsAMAuuf2kRADD/USYAAGuUCQDAGmUCALBGmQAArFEmAABrlAkAwBplAgCwRpkAAKxRJgAAa5QJAMAaZQIAsEaZAACsUSYAAGuUCQDAGmUCALBGmQAArFEmAABrlAkAwBplAgCwRpkAAKxRJgAAa5QJAMAaZQIAsEaZAACsUSYAAGuUCQDAGmUCALBGmQAArLleJm1tbfL5fCoqKlJTU9Nl80eOHFFZWZlKS0u1fft2nT171oWUAIBYXC2TYDCo+vp6NTc3q6WlRQcOHFB/f78zf/78ef34xz9WQ0ODWltblZWVpV/96lcuJgYAXImrZdLV1aX8/HylpaXJ4/GouLhY7e3tznw4HNauXbuUkZEhScrKytLnn3/uVlwAQBSulkkoFFJ6eroz9nq9CgaDzvjWW2/V6tWrJUljY2NqaGjQqlWr4p4TABDbQjcPPjU1paSkJGdsjIkYTxseHtb3vvc9LVu2TOvWrZvVMXp7e61zXk+BQMDtCDGRzw757CV6RvJd5GqZZGZmqqenxxkPDg7K6/VG7BMKhfToo48qPz9fO3bsmPUxsrOzlZqaap31eggEAsrLy3M7RlTks0M+e4me8UbMNz4+fk0vwl19m6ugoEDd3d0aGhrS6OioOjs7VVhY6MxPTk6qurpaJSUlqquru+JVCwDAfa5emWRkZKi2tlYVFRUKh8MqLy9Xbm6uqqqqVFNTo4GBAf31r3/V5OSkOjo6JF280tizZ4+bsQEAl3C1TCTJ7/fL7/dHbGtsbJQk5eTk6OTJk27EAgDMgutfWgQAzH+UCQDAGmUCALBGmQAArFEmAABrlAkAwBplAgCwRpkAAKxRJgAAa5QJAMAaZQIAsEaZAACsUSYAAGuUCQDAGmUCALBGmQAArFEmAABrlAkAwBplAgCwRpkAAKxRJgAAa5QJAMAaZQIAsEaZAACsUSYAAGuUCQDAGmUCALBGmQAArFEmAABrlAkAwJrrZdLW1iafz6eioiI1NTVdNn/ixAmtX79excXFqqur08TEhAspAQCxuFomwWBQ9fX1am5uVktLiw4cOKD+/v6IfZ566ik999xz6ujokDFGBw8edCktACCahW4evKurS/n5+UpLS5MkFRcXq729XY8//rgk6dSpUxobG9Ndd90lSVq/fr1efPFFbdq06arPbYyRJF24cOE6pZ8b4+PjbkeIiXx2yGcv0TPeaPmm/2ZO/w2dKVfLJBQKKT093Rl7vV59/PHHUefT09MVDAZn9NzhcFiS1NfXN0dpr4/e3l63I8REPjvks5foGW/UfOFwWDfffPOM93e1TKamppSUlOSMjTER46vNx7Jo0SItXbpUycnJM34MAPyvM8YoHA5r0aJFs3qcq2WSmZmpnp4eZzw4OCiv1xsxPzg46IxPnz4dMR/LggULtHjx4rkLCwD/I2ZzRTLN1Q/gCwoK1N3draGhIY2Ojqqzs1OFhYXO/O23367U1FQFAgFJ0uHDhyPmAQCJIcnM9lOWOdbW1qbXXntN4XBY5eXlqqqqUlVVlWpqapSTk6OTJ09q586dOn/+vO68807t3btXKSkpbkYGAFzC9TIBAMx/rn9pEQAw/1EmAABrlAkAwBplAgCwNq/K5FpvCvnZZ59p8+bNWrNmjR577DGNjIxIks6dO6etW7eqpKREmzdvjvhOSzzzBQIBlZeXq6ysTN/5znd06tQpSdIHH3yge+65R2VlZSorK9Ozzz7rSr5Dhw7p3nvvdXLU19dLin5e45nvX//6l5OrrKxM999/v+6++25J8T9/055++mm98847zjhR1l+0fImy/qLlS5T1d6V8ibT+jhw5orKyMpWWlmr79u06e/aspPitP5l5YmBgwNx3333mzJkzZmRkxPj9fvO3v/0tYp8HHnjA/PnPfzbGGPPss8+apqYmY4wxW7duNb/73e+MMca89NJL5oUXXjDGGPOTn/zEvPbaa8YYYw4dOmSeeOIJV/Ldd9995sSJE8YYY37729+a6upqY4wxb7zxhnn11VevOdNc5du9e7dpa2u77Dmjndd455s2OTlptmzZYlpbW40x8T9/AwMDZtu2bSY3N9e8/fbbzvZEWX/R8iXK+ouWL1HWX7R809xcf8PDw+ab3/ymGRgYMMYYs2/fPvPTn/7UGBOf9WeMMfPmyuS/bwrp8Xicm0JOu9JNIdvb2xUOh/Xhhx+quLg4Yrskvffee/L7/ZKktWvX6o9//KNzT6945btw4YKeeOIJLVu2TJKUlZWlzz//XJJ07Ngx/elPf5Lf71d1dbWzPZ75pnMcOnRIfr9fTz75pM6ePRvzvMY737S3335bX/jCF5zfaTzPn3TxleO3vvUtlZSUONsSZf1Fy5co6y9avukcbq+/WPmmubn+wuGwdu3apYyMDEn//3uM1/qT5tHbXFe6KeR/3/Qx2k0hz5w5o1tuuUULFy6M2H7pYxYuXKhbbrlFQ0NDcc2XkpKisrIySRfvRfbSSy9p1apVkqTFixfrkUceUVtbm1auXKna2tprymaTb/rn7du3q7W1VV/+8pe1e/fumOc13vkkaXJyUq+++qp++MMfOtvief4kqbKyUt/+9rcjtiXK+ouWL1HWX7R8UmKsv1j5JPfX36233qrVq1dLksbGxtTQ0KBVq1bFbf1J86hMrvWmkJfuJynqjR+NMVqw4NpOie1NKy9cuKAnn3xSExMT2rZtmyRp9+7dKioqkiQ9/PDD6u/v1/DwcNzzvfzyy8rLy1NSUpIqKyv1/vvvz+q8Xu98kvT+++/rq1/9qrKyspxt8Tx/0STK+rsat9dfLImw/q4mUdbf8PCwtm7dqmXLlmndunVxW3/SPCqTS2/6ONObQt52220aHh7W5OTkZY/zer06ffq0JGliYkIjIyPO/60Sr3ySNDIyosrKSk1MTOiVV15RcnKypqam9Morrzi5p910001xzTc8PKz9+/c7240xuummm2Ke13jmm3bkyBH5fD5nHO/zF02irL9YEmH9RZMo6+9qEmH9hUIhbdq0SVlZWdqzZ4+k+K0/aR6VybXeFDI5OVkrVqzQu+++K0lqaWlxHrdy5Uq1tLRIkt59912tWLFCycnJcc0nXfzfJJcsWaJ9+/Y59x1bsGCBfv/736ujo8PJvXz5cnk8nrjm83g8ev311/XRRx9Jkt58802tXr065nmNZ75pf/nLX7RixQpnHO/zF02irL9YEmH9RZMo6+9q3F5/k5OTqq6uVklJierq6pyrj3itP0nz519zGWNMa2ureeCBB0xRUZFpaGgwxhhTWVlpPv74Y2OMMSdOnDAbNmwwxcXF5gc/+IEZHx83xhjz6aefmi1btpiSkhLz3e9+1/z73/82xhhz5swZs23bNuPz+czGjRvNP//5z7jnO378uFm6dKnx+XymtLTUlJaWmsrKSmOMMX19fWbjxo3G5/OZLVu2mM8++yzu+Ywx5sMPPzQPPvigWbNmjamurjbnzp0zxkQ/r/HOZ4wxubm5ZmxsLOL54n3+pj3zzDMR/9onUdbflfIl0vq7Uj5jEmf9RctnjPvrr7Oz02RlZTm/w9LSUrNjxw5jTPzWHzd6BABYmzdvcwEAEhdlAgCwRpkAAKxRJgAAa5QJAMAaZQIAsEaZAACsUSYAAGv/B2oS18n6YPTEAAAAAElFTkSuQmCC\n"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "## try to plot the results of x1, x2, x3\nplt.figure()\nplt.subplot(311)\nplt.title(\"redshift against luminosity distance\")\nplt.xlabel(\"d_L(Mpc)\")\nplt.ylabel(\"z\")\nplt.plot(x1, z, \"-\")\nplt.grid(True)\n\n##repeat for x2, x3\nplt.subplot(312)\nplt.title(\"redshift against age\")\nplt.xlabel(\"Age(Gyr)\")\nplt.ylabel(\"z\")\nplt.plot(x2, z, \"-\")\nplt.grid(True)\n\nplt.subplot(313)\nplt.title(\"redshift against comoving distance\")\nplt.xlabel(\"D_A(Mpc)\")\nplt.ylabel(\"z\")\nplt.plot(x3, z, \"-\")\nplt.grid(True)\n\nplt.tight_layout()\nplt.show()\n\n",
   "execution_count": 4,
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOyde3hcR333P1+ttKu7JUu24lsix3biXIE4N0igzgXIG6ChLe0LLS0U2rxAS8vbUm5teQIFSgstkJK3kIYQWi4mDYQmFCgBYiBASOKEXJx7bCe2Y1u2JEtaXVYr6ff+MbPSSl7Jsqz1SvLv8zznOXNm5sz5zZzLd+eyMzIzHMdxHGeuUVZqAxzHcRynEC5QjuM4zpzEBcpxHMeZk7hAOY7jOHMSFyjHcRxnTuIC5TiO48xJXKCOYyRtlLRrGvGukfTlKcK3StoY3ZL0RUmdku6ZRXNnhKS0pJOPwXVM0toipPtSSU/MdroTrvE5SX9bpLTHPWP5z4rjHI7yUhvgzH/M7Iy8w4uBlwMrzaw3foy+bGYrS2Rb7dGmIekmYJeZ/c3RW3RkmNlPgVOLfI235dzFvl8TnpWCSGoFtgMVZjZUDDuc+YHXoBYIkubKj42TgB1m1ltqQxzHmd+4QM1jJO2Q9F5JDwG9ksolLZf0DUn7JW2X9Gd58ask3RSb3x4FzpuQ3nsl7ZbUI+kJSZflBScl/XsM2yrp3Al2XC7prcANwItj09ongO8Cy+NxWtLyAvl4laQHJHVL2inpmgnhfyDpWUntkv42d70Ydr6kX0g6KGmPpM9KSuadO9r0FvN+naT/jvn4paQ1MUySPiWpLdrxsKQzJV0N/B7wnmj/7dO4L5sl/VHe8Zsl3TXBpndIeira8XeS1kj6ebz2zbk8FGgi2yHp3ZIektQl6euSKvPC/1jS05I6JN2WK+/J8pdXLh+RVFPofknqk9SUd41z4vNVUSDvh3vGJt67+6I9+yT9c4z2k7g/GG14cSyfH8Vn4ICkr0hqOIJyuUrSr+K1npF0RfRfJOkL8dnZHcshcbh77BwjzMy3eboBO4BfAauAKsIPji3AB4EkcDKwDXhljP9x4KfA4njOI4SmKwjNSDuB5fG4FVgT3dcAA8CVQAL4e+DuCXZcHt1vBu7KC9uYu8YU+dgInBXtPxvYB7w2hp0OpAlNh0ngk0A273obgAsJzdWtwGPAu/LSNmBtdN8EtAPnx/hfATbFsFfGsmsABJwGLMs77yOHyUP+dTYDf5QXNrFMDPgvoB44A8gAP4z3axHwKPCmQuUXy/oeYHm8j48Bb4thlwIHgHOAFPAvwE+OJH+F7hfwHeDtecefAv5lknKY9Bkr8Kz8Avj96K4FLsx79gwozztvLaHpOAUsIYjYp6dZLucDXfH8MmAFsD6G3Qp8HqgBlsY0/k+p323fwuY1qPnPtWa208z6Cb9Wl5jZh81s0My2Af8GvD7G/R3go2bWYWY7gWvz0hkmvPynS6owsx1m9kxe+F1m9h0zGwb+A3jBbGXAzDab2cNmNmJmDwFfA34tBr8OuN3M7jKzQYL4Wt65W8zsbjMbMrMdhI/NrzE5t5rZPRb6Nr4CvDD6Z4E6YD0gM3vMzPbMVh4L8I9m1m1mWwkf8e+b2TYz6yLUYl40xbnXmtnzZtYB3J6Xh98DbjSz+80sA7yfUJtt5ejy9yXgjQCxdvEGwjNQiKmesYlkgbWSms0sbWZ3TxbRzJ42szvMLGNm+4F/5tD7PFm5vJVQLnfEZ2y3mT0uqYXwo+tdZtZrZm0E8X09zpzABWr+szPPfRKheeZgbgM+ALTE8OUT4j+bc5jZ08C7CLWlNkmbNL45bm+euw+o1Cz1e0m6QNKdsdmoC3gb0FzIZjPrI9SCcueeIunbkvZK6gY+lnduISbmozam+yPgs8B1hPxfL6l+FrI3Gfvy3P0Fjqca3FEwD4Syyr+naUJZrTjK/P0X4YfLakItpMvMJhuhOekzVoC3AqcAj0u6V9KrJ4soqSU+k7vjff4yh97nycplFfAMh3ISUAHsyXtfPk+oSTlzABeo+U/+dPQ7ge1m1pC31ZnZlTF8D+FlzXHiuITMvmpmFxNeXAP+YZbtm4yvArcBq8xsEfA5QjMUBJtHR5RJqgKa8s79V+BxYJ2Z1RMEWcwAM7vWzDYQmhVPAf7qCPKQTy9QnXd8wkzsmQHPE+4dALFPqQnYDVPmL59D8mpmA8DNhFrU7zN57QkO84xNSPcpM3sDQRD+Abgl2lyovD8W/c+K9/mNTP8+7wTWTOKfAZrz3pd6m8ZIQ+fY4AK1sLgH6FEY7FAlKRE7+nMd1TcD75fUKGkl8M7ciZJOlXSppBShv6kfGJkFm/YBTZIWTRGnDugwswFJ5wO/mxd2C/AaSS+JAweuYfyHqQ7oBtKS1gNvn4mRks6LNbkKgsAMMJb/fYT+oenyK+A3JVUrDNB460xsmgFfA/5Q0gvjffwY8Esz23GY/OUz2f36d0Jf2q8ztUBN+oxNRNIbJS0xsxHgYPQeAfbHfX6Z1xH6IrskraCwuE7GFwjlcpmkMkkrJK2PTZzfB/5JUn0MWyNpqiZi5xjiArWAiP1Drya0vW8ndJjfQOh4B/gQocllO+HFzP/QpAgd3AcITSVLCX0YR2vT44QP57bYjHLIKD7gHcCHJfUQ+phuzjt/K+Ejt4nw6zwNtBF++QK8myBoPYT+tq/P0NT6eH4noYzagU/EsC8QmrgOSvrWNNL6FDBI+Nh/idDXVXTM7AfA3wLfIJTVGsb6U6bKX34aBe+Xmf2MIBr3m9lUzXZTPWMTuQLYKikNfAZ4vZn1x2bcjwI/izZcGNM9hzDY4b+Bbx6mOPLzdA/wh4T70gX8mLGa5h8QBt88SiibW4Bl003bKS4y8wULnfmDpFrCr+11Zra91PYcT0j6EfBVM7uh1LY4xwdeg3LmPJJeE5vLagjDzB8mDCt2jhGxmfgcZl5DdZwjxgXKmQ9cRRgA8DywjtAU5FX/Y4SkLwE/IAzH7im1Pc7xgzfxOY7jOHOSos7fJmkHofN6GBgys3OnPsNxHMdxAsdigtFLzOzAdCI2Nzdba2vrjC/U29tLTU3NjM9fqHi5TI6XTWG8XArj5VKY3t5eHn/88QNmtmQ2050rM2AD0Nrayn333Tfj8zdv3szGjRtnz6AFgpfL5HjZFMbLpTDzqVyGR4yDfYN09g3Sng77jt7shOOw33jqUv7i5afM+FqbN2/mkksumervBzOiqH1QkrYT/ltgwOfN7PoCca4GrgZoaWnZsGnTphlfL51OU1t71Mv/LDi8XCbHy6YwXi6FKVW5mBkDw9AzaKQHje5BI501egYhPWj0ZG00LOfuy04+BUplAmqToq5C1CXF2UsSXH7SIZPTT5t0Os1rXvOaLbPdjVPsGtTFZrZb0lLgDkmPm9lP8iNE0boe4Nxzz7Wj+XUyn37dHEu8XCbHy6YwXi6Fma1yyQwN09mbHa3BtPcO0tk7OOlxZ2+WweHCE7tUJERjdZLFNUmaFyc5pTbJ4uokjTVJFldXsLg2FY8rWFyTpLE6SWXF7K4osnnz5llNL0dRBcrMcnOAtUm6lTDt/U+mPstxHGf+MDJiHOwfE5uO3rGts3eQjr7x7s7eLOnM5AsFN1RXjArMysZqXrCyIYhNTQWN1UmaapOjgtRYk6QuVY40o+kn5zxFE6j4p8oyM+uJ7lcAHy7W9RzHcY4WM6NvcHhMZPrGajIdvYNsfTrD13beNyZAfVkO9g0yMklbWlVFgsU1Y2KyurmGxTWpIDY1oaaTH95QVUF5wv+emqOYNagW4Nao7OWEKVK+V8TrOY7jjGNwaISDfbEWk84XnCwdvRk6+rLjBKijb5DBocJNaYkyUVMOy7J9NNZUcOoJdaFGE8Ul13y2OM9dlfTFeY+GoglUXCxv1ha1cxzH6R8cpr03Q0dv6KfpSI/12XSk4z4vvGdg8qa0+sryUTFZ3lDJGcvrx9Vmcs1sOQGqryznxz/+MRs3vuwY5vj4Zk4NM3cc5/jBzOjJDOUJSxCXnNh05PXntKfDvj87XDCt8jKNiktTbZKzGhvCAIGaFItro8jk1W4aqiuo8Ka0OY8LlOM4s8LYYIHMqKC05zWfjdVusqO1nOxw4c6byooymmpSo4KydkltcNeO9ds01SZjf06o3SzUgQLHMy5QjuMUZHBoZNyotFCzyYwTnnwBmmqwQF2qPIhLTZIVDZWctaKexTUpmmrGajX5NaDqpH+aHBcoxzlumNh/0zlOeMb33+zr6qP/e98tmI4EDVXhPzVNNalQu1mdHCc2udpPU21oTkuV+2AB58hxgXKceUrf4BDtUVja05lR92g/zgz6bxbXJDlzxSLW1Axy9iknj/bfLK4Z2zdUJ0mUeXOaU3xcoBxnjjA4NEJH7yAH0pkxoUkPciAdBKijd5ADeWI0meAU6r9pnCAyh+u/CTMmrDsW2XacSXGBcpwiMTxioxNztkexaY/iM7HWcyCdmXRIdEVC45rMTm6uoakmSVNtKu7Hu73/xlko+JPsONPEzOjuHwpiEwXmQG60WjozWrsZbVbrG6TQXMxlYnTKmqaaFGcsr6e5dkyAmmpSNNfmajopH6HmHLe4QDnHNZkhY2dHX2hWi2JzYGJtJ9aAphoWXV9ZTnNtiqbaMJ3Nua2Lac7VcqLYNMdajvfhOM70mJZASfoh8E9m9p08v+vN7OqiWeY4M2BkxOjqz3IgnWF/rOEc6MnQ3pvhQE9oSjsQ/dt7MwxkR+AHdx6STnUyMdpPs2xRJWeuqB/fpFYThKe5NkVjdZJkuf/p03Fmm+nWoFYD75V0npl9KPr58u3OMWFoeISOvsFRgckXm3wROhCb14YK/BknEUeqNdeG5rM18Y+fXW27OO/s02iuHT802vtxHKf0TPctPAhcBlwr6XbgjcUzyTkeyAwNx2a0fKHJFBShyfpykokymmuTNNelOCHWcoIAhdrNktoUzXXhuKGqgrICzWqbN7ex8dxVxyDHjuMcKdMVKJnZEPAOSW8G7gIai2aVMy/pHxweE5uezOjw6FyT2pgIZeieZMRadTIxWss5qamaDa2No8c58cmJ0kJeB8dxnOkL1OdyDjO7SdLDwJ8UxyRnLjE8YnT0DrK/JwhPW/cA+9OZcNyToa0nCM7+ngw9kyzCVl9ZPlqTOe2EeprXhsED+WKzJNZ6vGnNcZwc0/oamNnnJxxvAd5SFIucY0JvZmhUdPb3jBeetp4xAWrvHWS4QJ9OXaqcJXWhCe205fX8WhSgJVFsmmtTNNeF0Ws+zY3jODPBf64uIIZHjPZ0FJh0hv3dYf/AYxn+c/f9tPUMjApP7+ChsxAkysSSKDIt9ZWcuXwRS+vHRGdpfYoltZU013lNx3Gc4uNfmXmAmXGwL8ve7gH2dQ/Q1p1hX/cA+3oG2NuVoa0n+O/vyRScTbqqHJb3d7OkLsVZKxvyxCaIT87dWJ0sOJDAcRynFLhAlZiegSz7ukMT276eAfblxKd7zN3WnWFw+NBlqBurK2ipr6SlvpL1J9TRUl/J0roUS+oqg/DUBQG6+2c/ZePGjcc+c47jOEeBC1SRMDM6+7Ls6epnz8EB9nT183zXAHsO9rM3rxZUqKmtNlVOS31oZjuvdTFL61O01AUhOmFRiqVRgCorvG/HcZyFiwvUDDAzugeGRsXn+VERCkKU2w9kx9d6ystES30lyxZVctryejaeunRUiMKWYml9JbUpvy2O4zj+JSxArs9nZ2cfuzr72dkR9/F4z8H+Q2o+ZWJUfE5fXs/lpy3lhEVVLF9UybKGsG+qTfkcbI7jONOkqAIl6QrgM0ACuMHMPl7M6x0JmaFhdnb0sW1/L89FAdqVJ0gTBWhRVQWrFlexdkktL1u3hGWLKlnWUMmyRVUsb6hkSW2K8oTPx+Y4jjNbFE2gJCWA64CXA7uAeyXdZmaPFuuaEzEzdh/s55n9vWzfn2ZHex/bDvSy/UCa3Z3940a81abKWdlYxcrGal68pomVjdWsiscrF1dRX1lxrMx2HMdxKG4N6nzgaTPbBiBpE3AVUBSB6h7I8kTHMM/9YgeP7+3hib09PLm3Z9zsBrWpclY31/CiVY38xotWcnJzDaubazhxcTUN1RU+bY7jOM4copgCtQLYmXe8C7igWBf7j188yyfuGQC2Ul9Zzvpl9fzGOSs49YQ61i2to7W5miW1KRchx3GceULJB0lIuhq4GqClpYXNmzfPKJ3mvhHefrqxbmk1jSkhZYAM9B+g71l49NnZs3m+kU6nZ1yuCx0vm8J4uRTGy6Uw6XS6KOkWU6B2A/nrGKyMfuMws+uB6wEk7b/kkkuORkqagQNHcf5CxctlcrxsCuPlUhgvl8I0AyfNdqKyQgvtzEbCUjnwJGEdqd3AvcDvmtnWolwwXPM+M/OFFCfg5TI5XjaF8XIpjJdLYYpVLkWrQZnZkKQ/Bf6HMMz8xmKKk+M4jrOwKGoflJl9B/hOMa/hOI7jLEwW2j9Lry+1AXMUL5fJ8bIpjJdLYbxcClOUcilaH5TjOI7jHA0LrQblOI7jLBBcoBzHcZw5yYIRKElXSHpC0tOS3ldqe4qBpBsltUl6JM9vsaQ7JD0V943RX5KujeXxkKRz8s55U4z/lKQ35flvkPRwPOdazZNpNyStknSnpEclbZX059H/uC4bSZWS7pH0YCyXD0X/1ZJ+GfPydUnJ6J+Kx0/H8Na8tN4f/Z+Q9Mo8/3n73klKSHpA0rfj8XFfLpJ2xOf8V5Lui36le4/MbN5vhGHszwAnA0ngQeD0UttVhHy+DDgHeCTP7x+B90X3+4B/iO4rge8CAi4Efhn9FwPb4r4xuhtj2D0xruK5/6vUeZ5muSwDzonuOsL/704/3ssm2lob3RXAL2MebgZeH/0/B7w9ut8BfC66Xw98PbpPj+9UClgd37XEfH/vgL8Avgp8Ox4f9+UC7ACaJ/iV7D1aKDWo0YlpzWwQyE1Mu6Aws58AHRO8rwK+FN1fAl6b5//vFrgbaJC0DHglcIeZdZhZJ3AHcEUMqzezuy08Sf+el9acxsz2mNn90d0DPEaYC/K4LpuYv9wcNBVxM+BS4JboP7FccuV1C3BZ/IV7FbDJzDJmth14mvDOzdv3TtJK4FXADfFYeLlMRsneo4UiUIUmpl1RIluONS1mtie69wIt0T1ZmUzlv6uA/7wiNr+8iFBbOO7LJjZj/QpoI3wongEOmllumv/8vIzmP4Z3AU0ceXnNBz4NvAfILXvdhJcLhB8w35e0RWGeVCjhe1TyyWKd2cPMTNJx+78BSbXAN4B3mVl3fvP28Vo2ZjYMvFBSA3ArsL7EJpUcSa8G2sxsi6SNpbZnjnGxme2WtBS4Q9Lj+YHH+j1aKDWoaU1Mu0DZF6vOxH1b9J+sTKbyX1nAf14gqYIgTl8xs29Gby+biJkdBO4EXkxoisn9OM3Py2j+Y/gioJ0jL6+5zkXAr0vaQWh+u5Sw8vfxXi6Y2e64byP8oDmfUr5Hpe6Um42NUBPcRuiozHVKnlFqu4qU11bGD5L4BOM7MP8xul/F+A7Me2ysA3M7ofOyMboXW+EOzCtLnd9plokI7dmfnuB/XJcNsARoiO4q4KfAq4H/ZPxggHdE958wfjDAzdF9BuMHA2wjDASY9+8dsJGxQRLHdbkANUBdnvvnwBWlfI9KXiizWLhXEkZvPQP8dantKVIevwbsAbKE9tu3EtrCfwg8Bfwg70EQcF0sj4eBc/PSeQuhQ/dp4A/z/M8FHonnfJY408hc34CLCW3nDwG/ituVx3vZAGcDD8RyeQT4YPQ/OX4ono4f5VT0r4zHT8fwk/PS+uuY9yfIG3k13987xgvUcV0uMf8Pxm1rzu5Svkc+1ZHjOI4zJ1kofVCO4zjOAsMFynEcx5mTuEA5juM4cxIXKMdxHGdO4gLlOI7jzElcoBzHcZw5iQuU4xwGSddIevckYTdJel0B/2V5yzhslGSS/igv/IXRr2C607DpB7llDxxnoeIC5TjF4S+Af8s7fgT4nbzjNxD+EDlT/oOwDITjLFhcoBynAJL+WtKTku4CTp1BEr8FfC/v+FmgUlJLXKrhCsJUL7nrbZb0mbhQ3COSzo/+tZK+GBd5e0jSb8VTbiOInOMsWHw2c8eZgKQNhDnXXkh4R+4HthzB+auBTjPLTAi6BfhtwvRD9wMTw6vN7IWSXgbcCJwJ/C3QZWZnxbQbAcysM6702mRm7UeaR8eZD3gNynEO5aXArWbWZ2bdhNrKkbAM2F/A/2aCQL2BMK/iRL4GowtT1sclMi4nzHdGDOvMi98GLD9C2xxn3uAC5TizTz9hgtFxmNlewkS/LydMvnlIlMMcT6QyXstxFiQuUI5zKD8BXiupSlId8JojPP9JwrIohfgg8F4LCwlO5H8DSLqY0KzXRVgF909yEXJNfLEf6wRgxxHa5jjzBhcox5mAmd0PfJ0wyu67wL2HOeXzknbF7Rdm1gs8I2ltgbR/bmbfmiSdAUkPENYiemv0+wjQGAdOPAhcEv03AHfb2BLljrPg8OU2HKcISPoNYIOZ/c00428G3m1m900z/meA28ysUFOh4ywIfBSf4xQBM7tVUlMRL/GIi5Oz0PEalONMA0nXARdN8P6MmX2xFPY4zvGAC5TjOI4zJ/FBEo7jOM6cxAXKcRzHmZO4QDmO4zhzEhcox3EcZ07iAuU4juPMSVygHMdxnDmJC5TjOI4zJ3GBchzHceYkLlCO4zjOnMQFypm3SNooadc04l0j6ctThG+VtDG6FZdY75R0zyyaOyMkpSWdXGo7HKcU+GSxznGPmZ2Rd3gxYUHBlWbWG4Xry2a2skS21R5tGpJuAnZNd2Z1x5kreA3KmRNImis/lk4CdsQ1nRzHKSEuUE7JkLRD0nslPQT0SiqXtFzSNyTtl7Rd0p/lxa+SdFNsfnsUOG9Ceu+VtFtSj6QnJF2WF5yU9O8xbKukcyfYcbmktwI3AC+OTWufICxYuDwepyUtL5CPV0l6QFK3pJ2SrpkQ/geSnpXULulvc9eLYedL+oWkg5L2SPqspGTeuZZb+DDm/TpJ/x3z8UtJa2KYJH1KUlu042FJZ0q6Gvg94D3R/tsnuRefibZ3S9oi6aUTyv1Lsdwfk/Se/KbVqe6Z4xwVZuabbyXZCMuV/wpYBVQRfjBtISyLngROBrYBr4zxPw78FFgcz3mE0HQFcCqwE1gej1uBNdF9DTAAXAkkgL8nrEabb8fl0f1m4K68sI25a0yRj43AWdH+s4F9wGtj2OlAmtB0mAQ+CWTzrrcBuJDQ3N4KPAa8Ky9tA9ZG901AO3B+jP8VYFMMe2UsuwZAwGnAsrzzPnKYPLwRaIrp/iWwF6jMK/cfA43ASuChvHKf8p755tvRbF6DckrNtWa208z6CTWiJWb2YTMbNLNtwL8Br49xfwf4qJl1mNlO4Nq8dIaBFHC6pAoz22Fmz+SF32Vm3zGzYeA/gBfMVgbMbLOZPWxmI2b2EPA14Ndi8OuA283sLjMbJHzILe/cLWZ2t5kNmdkO4PN55xbiVjO7x8JS718BXhj9s0AdsJ6wjM5jZrbnCPLwZTNrj3b8E6EsT43BvwN8zMw6zWwX48v9cPfMcWaMC5RTanbmuU8iNKcdzG3AB4CWGL58Qvxncw4zexp4F6G21CZp04TmuL157j6gcrb6vSRdIOnO2MTVBbwNaC5ks5n1EWpBuXNPkfRtSXsldQMfyzu3EBPzURvT/RHwWeA6Qv6vl1R/BHl4d2y+64rlvmiyPHBk98xxZowLlFNq8lfM3AlsN7OGvK3OzK6M4XsITXs5ThyXkNlXzexiwkfTgH+YZfsm46vAbcAqM1sEfI7QzAbB5tERgJKqCE1pOf4VeBxYZ2b1hI+7mAFmdq2ZbSA0K54C/NV08hD7m95DqCk1mlkD0DVZHhh/Dw53zxxnxrhAOXOJe4CeONihSlIidvTnBkPcDLxfUqOklcA7cydKOlXSpZJShP6mfmBkFmzaBzRJWjRFnDqgw8wGJJ0P/G5e2C3AayS9JA5+uIbxAlQHdANpSeuBt8/ESEnnxZpcBdBLKINc/vcR+oamsn8I2A+US/ogkF/7yi/3FcCf5oUd7p45zoxxgXLmDLF/6NWEfpXtwAHCqLqcOHyI0Ky3Hfg+oS8pR4rQmX+A0Ay2FHj/LNj0OKFPaVtswjpkFB/wDuDDknoIfUw3552/lSCkmwg1kTTQBmRilHcTBK2H0Hfz9RmaWh/P7ySUUTvwiRj2BULf3EFJ3ypw7v8A3wOejOcOML4Z78PALkK5/4AgupmYv8PdM8eZMTKbTguG4zizgaRa4CChSW97qe2ZCZLeDrzezKYazOE4R43XoBynyEh6jaRqSTWEYeYPE4a2zwskLZN0kaQySacShqHfWmq7nIWPC5TjFJ+rgOfjto5Q+5hPTRdJwvD3HuBHwH8B/6+kFjnHBUVt4pO0g/BQDwNDZnbu1Gc4juM4TuBYzH92iZkdOAbXcRzHcRYQc2WCTgCam5uttbV1Rudmh42u/kGaa1OzaxTQ29tLTU3NrKdbSjxP8wPP0/zA8wRbtmw5YGZLZtOGYguUAd+XZMDnzez6iRHiZJZXA7S0tPDJT35yRhe67ZlBvvlUluXNCf747BR1yRn917Eg6XSa2tqjXvVgTuF5mh94nuYHnie45JJLnj18rCOkmBP9ASvifinwIPCyqeJv2LDBZsrIyIh96efbbd0HvmMXfPQHtvmJthmnNZE777xz1tKaK3ie5geep/mB58kMuM/m02SxZrY77tsIw1LPL9a1JPEHL27lm+94CbWV5bzpxnt4zy0P0tWfLdYlHcdxnCJSNIGSVCOpLucGXkFYHqGonLliEd9+58W8feMavnH/bl7xqR9z6wO7GBmZT6N6HcdxnGLWoFqAuyQ9SJiv67/N7HtFvN4olRUJ3nvFem59x0toqa/k/379QX7zX3/O/c91HovLO47jOLNA0QZJWFgXZtbW3JkJZ69s4FvvuIhvPrCbf/ze4/zm//s5Lz+9hT+/bB1nrvCpwhzHceYyc2qYeTEoKxOv27CS/3XmCfzbT7dx46ozTO8AACAASURBVF3befWj+7j8tKW889J1vGBVQ6lNdBzHcQpw3Ex1VJMq512Xn8Jd77uUv3j5Kdy7o5OrrvsZv/WvP+f2B58nOzwbKzM4juM4s8WCr0FNpL6ygj+7bB1/eFEr/3nfLr70ix2882sP0FKf4g3nn8hvnbOSVYurS22m4zjOcc9xJ1A56ioreMvFq3nzS1rZ/GQbX/zZDj79g6f49A+e4oLVi0Oz4FnLqE0dt0XkOI5TUo77r29Zmbh0fQuXrm9hV2cft96/m2/cv4u/uuUhPvhfW9l46hJOLBvinIEs9ZUVpTbXcRznuOG4F6h8VjZW887L1vGnl67l/uc6ufWB3fzP1n18tyfDjVvv4CVrmnnlGSew8dQlLG+oKrW5juM4CxoXqAJIYsNJi9lw0mI+/OtncuN//Yi21HK+98hePnDrwwCsWVLDS9ct4WWnNHPB6iZqvCnQcRxnVvGv6mEoKxNrGxP80cbTeP//Ws9TbWl+8uR+fvrUATbd+xw3/XwHFQnxolWNnNvayHmtiznnxEYWVXtzoOM4ztHgAnUESOKUljpOaanjj156MgPZYbY828lPntrP3ds6uP4n2/h/m58B4NSWOja0NnLuSY2cvXIRq5trSZTN3gzrjuM4C51pCZSkHwL/ZGbfyfO73syuLppl84DKigQXrW3morXNAPQPDvOrnQe5b0cH9z7bye2/ep6v/vI5AKqTCc5YXs+ZKxZx1opFnLliEWuWuGg5juNMxnRrUKuB90o6z8w+FP18+fYJVCUTvHhNEy9e0wTA8IjxVFsPD+/q4pHdXTzyfDeb7tnJF7M7QvyKBOtaalm3tI5TWmo5paWOdS21LF9URZkLl+M4xznTFaiDwGXAtZJuB95YPJMWDokysf6EetafUM9vn7sKCKK1bX+ah3d38cjubp7c18NPn9rPN+7fNXpeTTLB2qW1rGupY82SWlqbqjmpqYaTmqp9MIbjOMcN0/3aycyGgHdIejNwF9BYNKsWMIkysa6ljnUtdfzmOWP+XX1Znmzr4cl9PTy1L82T+3rY/MR+btmya9z5S+pStDZV09pUQ2tzEK2TFteworGKxuoKJK95OY6zMJiuQH0u5zCzmyQ9DPxJcUw6PllUXcF5rYs5r3XxOP+egSzPtvfxbHsfO9p72XGgl2fb+/jxk/v5zwniVVWRYHlDJcsbqljRUMXy0a2SlQ3VnLCokmT5cTP9ouM485xpCZSZfX7C8RbgLUWxyBlHXWUFZ8ZBFRPpGxwaFa/nD/aHrauf3Z39PLanhwPpzLj4EjTXpmipT5EYHOB77Q+xtC7FkvpKltalwlZfyZLalAuZ4zglxzs05jHVyXJOW1bPacvqC4YPZIfZ2zXA8wf72ZUTsIP9tPVk2PZ8mh881kZ7bwYrsNjw4ppkEK+4NdUkWVyTYnFNRdwng19tkrpUuTctOo4z67hALWAqKxK0Noe+qols3ryZjRs3MjQ8QnvvIG3dGdp6BtgX9209Gdq6M+zvGeCZtjTtvYNkhgovSVKREI3VySBatUkaq5PjBK2+qoJFBbbyhNfSHMeZHBeo45zyRBkt9ZW01FcCU68y3Dc4RHt6kI7eQTr6BumI7vbeQTrjvqM3w/MHu2lPZ+geGJoyvdpUOYuqcgJWXlDE6vP2dalyaivLqU2VU5Ms96H4jrPAcYFypk11spzqxeXTXi8rOzxCZ98g3f1ZuvK3vixd/UPj/Lr7s+w40EdXf5aD/YMMZA+/gGRtKohVTrSyff18fdcWaqJ/XeX48HBcMequTiaoSZWTKi/zJkrHmYO4QDlFoyJRxtK6SpbWVR7xuZmh4VHh6urP0j0wRHpgiHQm7HsyQ/Rmxvx6MkM83wPP7E+PhqczQwX71yZSpjACsjoVRKuqIghXvrsqmaC6IkF18tB4ubBRdzIRxDyZoMKbMR1nxrhAOXOSVHmCpXWJIxK30K/2a6PHZkbf4HAQsDxxS2ey9AwEgevLDtM/OExvZpj+7BB9g8NxC/H392ToHRyif9R/+IjyUZEQleUJUhUJKivKqMztyxPRPZl/8HvuuSz779tJZUUQxInxUxPiuyA6CwkXKGfBIomaVDk1qXJaCg90PGJGRoyBoSBUOdEaL2BjItcf3QPZEfqzw2SywwwMheOBbIjb0TsY/AaHGRgK/gPZYUbya36PPjRt+xJlorK8bFT8kuVlJBNlJMvLSJWH/Zg7MS4sVT4hXiLEOfS8uE8kSFUUTj+Z8GZT5+gpqkBJugL4DJAAbjCzjxfzeo5TbMrKFJvvivfqmBnZ4SCEP/rxT9lw3oUMZIfpz46J20B2TNAyef6jcYaGyWRHGBweYXBomMGhETJDIwwOjZDODOWFBf9MjDM4PDKtZtHpkCwvI1VAvDL9/TRs/RnJhKhIlI1uyXJRXjbmzg+rGBdXJMuDu7xszD0aliijPO+csXAdkl4yUeaDbeYwRXvLJCWA64CXA7uAeyXdZmaPFuuajrMQkESyPHx4G1Jl0x6UMhvkxDEnXoN54pXJE7kxcRvOizcWlpksbGiEvW391FdVkB0aITs8Qu/g8Kg7bDbOPRjdsyWcE0mUKQhdooyK8vHiVR7d5YkyKspEoiwch30Q1ERCdOwf4Pa2B4Nf9C8vE4mEqCgri36iPIrqOPfofuy8UXciXCcxwT/45dkS4+b8ysSCqMEWswZ1PvC0mW0DkLQJuApwgXKcOUq+OJIqzjVCX+H5R3ze8EgQrsHhkShodoioDQ6PMDR8+HiF4wb30MgIg0P554Tzctfvzw4zNGIMDY+M+vX0jvBcfztDIyHNXPjQSHAPjxRJXaegYoLgJaKwjRPQKGivOL2Fd1627pjbeDhkRfpZIul1wBVm9kfx+PeBC8zsTyfEuxq4GqClpWXDpk2bimLP0ZBOp6mtrS21GbOK52l+4HmaHxwuTyNmjBgMG2E/AkM5v5HgP2xBhEcMhqJ/OMdi2Ph4ubSGRuON+Y+PX/g6o26DM5oSvLJ1/CrgR3qfLrnkki1mNqvLMJVcoCacsx94tigGHR3NwIFSGzHLeJ7mB56n+YHnCU4ysyWzaUAxm/h2A6vyjldGv0mZ7czNFpLum+1fBqXG8zQ/8DzNDzxPxaGYf5q4F1gnabWkJPB64LYiXs9xHMdZQBStBmVmQ5L+FPgfwjDzG81sa7Gu5ziO4ywsivo/KDP7DvCdYl7jGHF9qQ0oAp6n+YHnaX7geSoCRRsk4TiO4zhHg0/c5TiO48xJXKAcx3GcOYkL1BRIWiXpTkmPStoq6c9LbdNsISkh6QFJ3y61LbOBpAZJt0h6XNJjkl5capuOFkn/Nz53j0j6mqQjX7ekxEi6UVKbpEfy/BZLukPSU3HfWEobj5RJ8vSJ+Ow9JOlWSQ2ltPFIKZSnvLC/lGSSmo+1XS5QUzME/KWZnQ5cCPyJpNNLbNNs8efAY6U2Yhb5DPA9M1sPvIB5njdJK4A/A841szMJI2FfX1qrZsRNwBUT/N4H/NDM1gE/jMfziZs4NE93AGea2dnAk8D7j7VRR8lNHJonJK0CXgE8d6wNAheoKTGzPWZ2f3T3ED56K0pr1dEjaSXwKuCGUtsyG0haBLwM+AKAmQ2a2cHSWjUrlANVksqBauD5EttzxJjZT4COCd5XAV+K7i8Brz2mRh0lhfJkZt83s6F4eDdhYoJ5wyT3CeBTwHuAkoymc4GaJpJagRcBvyytJbPCpwkP3eHXVZ8frAb2A1+MzZY3SKoptVFHg5ntBj5J+OW6B+gys++X1qpZo8XM9kT3XqCllMYUgbcA3y21EUeLpKuA3Wb2YKlscIGaBpJqgW8A7zKz7lLbczRIejXQZmZbSm3LLFIOnAP8q5m9COhl/jUbjSP2y1xFEN/lQI2kN5bWqtnHwv9cFsx/XST9NaFr4CultuVokFQNfAD4YCntcIE6DJIqCOL0FTP7ZqntmQUuAn5d0g5gE3CppC+X1qSjZhewy8xytdtbCII1n7kc2G5m+80sC3wTeEmJbZot9klaBhD3bSW2Z1aQ9Gbg1cDv2fz/g+kawo+jB+O3YiVwv6QTjqURLlBToLDi1xeAx8zsn0ttz2xgZu83s5Vm1krodP+Rmc3rX+ZmthfYKenU6HUZ83/dseeACyVVx+fwMub5wI88bgPeFN1vAv6rhLbMCnH18PcAv25mfaW252gxs4fNbKmZtcZvxS7gnPiuHTNcoKbmIuD3CbWMX8XtylIb5RTkncBXJD0EvBD4WIntOSpibfAW4H7gYcK7WvKpZ44USV8DfgGcKmmXpLcCHwdeLukpQk3x46W08UiZJE+fBeqAO+J34nMlNfIImSRPJcenOnIcx3HmJF6DchzHceYkLlCO4zjOnMQFynEcx5mTuEA5juM4cxIXKMdxHGdO4gLlOAWQ9No4g/P6o0zn05JeFt3lkj4WZ/HO/W3hr48wvaSkn8T5+RxnQeMC5TiFeQNwV9zPCElNwIVxIk6AjxCmLTrLzF4IvBSoOIL0ys1skDAD+P+eqV2OM1/w/0E5zgTi3ItPAJcAt5vZqZLKCH/GvBTYCWSBG83sFkkbgH8GaoEDwJvNbI+kq4HlZnZNnNtsJ9AaZ8afeM0PAx1m9ul4/FHCFEAPAn8HdALrzewUSS8A/t7M/E/jzoLGa1COcyhXEdaWehJojwL0m0ArcDphdpEXw+hcjf8CvM7MNgA3Ah+N6VwE5CblXQs8V0icIjcCfxDTLCNMQ5WbI/Ec4M/N7JR4/Ahw3tFn03HmNt6O7TiH8gbCAogQJtR9A+Fd+U8zGwH2Srozhp8KnEmY4gbCwoK5pSSWEZYBOQRJf0hYNLIJeImZ7ZDULulFhOUnHjCz9pjmPWa2PXeumQ1LGpRUN4XgOc68xwXKcfKQtJjQjHeWJCMIjgG3TnYKsNXMCi0x3w/klml/GjgxJypm9kXC+lWPxGtAWEDyzcAJhBpVjt4CaaeAgWlnzHHmId7E5zjjeR3wH2Z2UpzJeRWwnbDa6G9JKpPUAmyM8Z8AlkgabfKTdEYMe4zQtEec4foLwGclVca4CSCZd+1bCctunwf8z2QGxsEXB+IyHI6zYHGBcpzxvIFDa0vfINRqdhGW8fgyYZbxrjiq7nXAP0h6EPgVY+s2/TdjQgbw14Tmv0ckPQD8lLDk+fMQlqoH7gRuNrPhKWy8JKbtOAsaH8XnONNEUq2ZpWMN5h7gosOtjyPpLuDVZnZwGumXEYTvt83sqSnifRN4XxzE4TgLFu+Dcpzp821JDYRmub+b5uJtfwmcCEwpUJJOB74N3HoYcUoC33Jxco4HvAblOI7jzEm8D8pxHMeZk7hAOY7jOHMSFyjHcRxnTuIC5TiO48xJXKAcx3GcOYkLlOM4jjMncYFyHMdx5iQuUI7jOM6cxAXKcRzHmZO4QDmO4zhzEhcox3EcZ07iAuUgaaOkXdOId42kL08RvlXSxuiWpC9K6pR0zyyaOyMkpSWdXGo7SomkE2M5JA4f+6iv9eY4k3vu+Lgvf+fIcYFyZg0zO8PMNsfDi4GXAyvN7PzpimARbas1s21Hk4akmyR9ZLZsOtaY2XOxHKZaa6pY1z5s+Zf6GXHmHi5QCwxJc2UJlZOAHWZWaLlyx3Gcw+ICtQCQtEPSeyU9BPRKKpe0XNI3JO2XtF3Sn+XFr4q1gU5JjxKWGM9P772SdkvqkfSEpMvygpOS/j2GbZV07gQ7Lpf0VuAG4MWxaecTwHeB5fE4LWl5gXy8StIDkrol7ZR0zYTwP5D0rKR2SX+bu14MO1/SLyQdlLRH0mfj2km5c03S2ui+SdJ1kv475uOXktbEMEn6lKS2aMfDks6UdDXwe8B7ov23T3IvzpB0h6QOSfskfSD6pyR9WtLzcfu0pFQM2yhpl6T3xOvukfRaSVdKejKm9YG8a0yV1mOSXp0Xtzw+A+dIao3lUB7DNkv6O0k/i+XwfUnN0ynvAvluknRbLLN7gDUTwvPL/0pJj8Zr7pb0bkk1hZ6Rad7Xt0l6Ksa5TpLywv84lklPvOY50X/S98OZQ5iZb/N8A3YQlhpfBVQRfnhsAT5IWFzvZGAb8MoY/+OE5cYXx3MeAXbFsFOBncDyeNwKrInua4AB4EogAfw9cPcEOy6P7jcDd+WFbcxdY4p8bATOivafDewDXhvDTgfShKbDJPBJIJt3vQ3AhYRFOFuBx4B35aVtwNrovgloB86P8b8CbIphr4xl1wAIOA1YlnfeR6awv46wpPtfApXx+IIY9mHgbmApsAT4OWHRw1y+h+L9qgD+GNgPfDWmcQbQD6yeRlofBL6SZ9OrgMfy7qUB5fF4M/AMcArhudkMfHw65V0g75uAm4Ea4Exg94T7n1/+e4CXRncjcM5kz8g07+u34/06MZbbFTHst6Md58V7uZZQs5/y/fBt7mwlN8C3WbiJQRjeknd8AfDchDjvB74Y3dtyL3E8vpoxgVoLtAGXAxUT0rgG+EHe8elA/wQ7ZixQBfL1aeBT0f1B4Gt5YdXA4BQfzHcRVqfNHU8UqBvywq4EHo/uS4En40exbEKaNzG1QL0BeGCSsGeAK/OOX0loAs2VTT+QiMd10d4L8uJvYUysp0prLdADVMfjrwAfjO5WDhWov8lL5x3A9460vAk/VrLA+jy/jzG5QD0H/B+gfkI6h31GJrmvF+cd3wy8L7r/B/jzAmlM+X74Nnc2b+JbOOzMc59EaCo5mNuADwAtMXz5hPjP5hxm9jThI3AN0CZpk8Y3x+Uvc94HVGqW+r0kXSDpztjs0gW8Dcg1OY2z2cz6CLWg3LmnSPq2pL2SugkfyGYmZ2I+amO6PwI+C1xHyP/1kuqnmYVVBPEoxHLyyjm688u13cYGL/TH/b688P6cjVOlFe/fY8BrJFUDv06oiU1GwXLgMOU9gSWEGk7BZ6oAv0X4UfCspB9LevFkEad5XyfLw2T343DvhzNHcIFaOFieeyew3cwa8rY6M7syhu8hvLw5ThyXkNlXzexiwotswD/Msn2T8VXgNmCVmS0CPkdomoFg88pcRElVQFPeuf8KPA6sM7N6wgdHzAAzu9bMNhBqiKcAfzXNPOwkNBcV4nlCeeY4MfrNhMOl9TVCbe4q4NEoWkfK4co7n/2EJspJn6l8zOxeM7uK0ET5LUKtBwqX79Hc151M6AvL85/q/XDmCC5QC5N7gB6FwQ5VkhKxoz83GOJm4P2SGiWtBN6ZO1HSqZIujZ3uA4Rf7iOzYNM+oEnSoini1AEdZjYg6Xzgd/PCbiHUCl4SO8mvYfyHqg7oBtKS1gNvn4mRks6LNbkKoJdQBrn872NyAYLQF7JM0rviQIY6SRfEsK8BfyNpSRyI8EFg0v+UHYbDpbUJeAWhDKaqPU3F4cp7lFjz+yZwjaRqSacDbyoUV1JS0u9JWmRmWcI9yy/fic/I0dzXG4B3S9qgwFpJJ3H498OZI7hALUDiB+PVwAuB7cABwsuae/E/RGiC2Q58H/iPvNNThEEUBwhNJ0sJ7fNHa9PjhA/rttiscsgoPkIfyIcl9RA+ujfnnb+VIKSbCL/u04S+skyM8m6CoPUA/wZ8fYam1sfzOwll1A58IoZ9ATg92v+tAnnsIfz36zWEsnsKuCQGfwS4D3gIeBi4P/rNhCnTMrM9wC+AlzDDcphGeU/kTwlNa3sJfXVfnCL53wd2xCa7txFGR072jMz4vprZfwIfJYh0D6G2tnga74czR5DZdFpeHGduIakWOEho+tleansWOl7eTinwGpQzb5D0mtiEVEMY9vwwYeSgUwS8vJ1S4wLlzCeuIgwGeB5YB7zevAmgmHh5OyXFm/gcx3GcOUlR522TtIPQOTkMDJnZuVOf4TiO4ziBYzGx6CVmdmA6EZubm621tXXGF+rt7aWmpmbG5x8r3M7ZZb7YCfPHVrdzdpnLdhowYkZCOio7t2zZcsDMlsymbXNl5msAWltbue+++2Z8/ubNm9m4cePsGVQk3M7ZZb7YCfPHVrdzdplLdpoZz+xP87On27nr6QPcva2d33jRCj581ZlHZaekqWYPmVmaxeyDkrSd8H8SAz5vZtcXiHM1YS44WlpaNmzatGnG10un09TW1h4+YolxO2eX+WInzB9b3c7ZpZR2mhl7eo0nO4d5vGOYxztGOJgJ3/0lVeL0pgQbWhKcvaT8qOy85JJLtsx2N06xBWqFme2WtBS4A3inmf1ksvjnnnuueQ1q7uB2zj7zxVa3c3Y5lnaOjBhPtvXwy20d3LO9g19u7+BAOvy/ekldigtWL+aitc1ctKaZE5uqZ81OSbMuUEVt4jOz3XHfJulWwvIGkwqU4ziOc2QMZIfZ+nw3DzzXyT3bO7h3RwedfVkAli+q5KXrmrlg9WLOX72Y1c015C2XNecpmkDFP/eVmVlPdL+CsI6N4ziOMwPMjOe7Brj/2U4eeO4g9z/XyaPPdzM4HKYzPHFxNZef1sL5qxdz4clNrGysmleCNJFi1qBagFtj4ZQDXzWz7xXxeo7jOAuK3swQj+4JtaOcIO3rDs11qfIyXrCygT+8qJUXndjIOSc2sLS+ssQWzy5FEygz2wa8oFjpO47jLCR6M0Nsfb6bh3d38cjuLh7e3cUz+9PkhgmsWlzFhSc3cc6JjbzoxAZOW1ZPRWJhTwY0p4aZO47jHA+kM0NsjSKUE6NtB3pHxWhpXYqzViziVWct46wVizh71SKW1i2s2tF0cIFyHMcpEsMjxo72Xu7ZO8SW7z/B43t7eHxvNzs7+kfjtNQHMXrNC5Zz1opFnLVi0YJrqpspLlCO4zizQEfvII/v6eaxvT08vqebJ/b18MTeHjJDYQBDmZ5mdXMNZ69s4H+fu4rTl9dz5orjs2Y0XVygHMdxpomZsadrgKfb0mHbn+aZtjTP7E9zID04Gq+pJslpy+p544Unsf6EOnp3P8nrr9xIZUWihNbPP1ygHMdxJjA0PMKzHX2jQvRMnhj1Dg6PxquvLGft0louXb+UdUvrWL+sjvUn1LOkLjUuvc3pZ1ycZoALlOM4xyVDwyPsPtjP9gO9PNvex472XnZE987OPrLDY7PsnFBfydqltfz2uatYs7SWtUtqWbu0luba5Lz+n9FcxwXKcZwFS3Z4hN2d/Wxv7+XZA73siEL0bHsfOzv6GBoZE6HqZILWphrWL6vjijNPYE0UoTVLa6lN+aeyFHipO44zbxkeMXZ19rGrs5+dHWG/q7OfnZ197O7sZ09XP3kaRE0yQWtzDacvq+fKs07gpKYaVjfXcFJTNUtqU14bmmNMS6Ak/RD4JzP7Tp7f9WZ2ddEscxznuCc7PEJbT4bdnf2HCtHBPp7v7Gf4+3eOxpdCc9zKxirOX72YlY1VnLi4mtbmGlqbarxJbp4x3RrUauC9ks4zsw9FP18d13GcGZMTnz0H+9nTNcDergGe7+qP+wH2dvWzvyczrgYE4X9DKxurOefERl7QMMTFL1zPysZqVjZWsayhklS5D0ZYKExXoA4ClwHXSrodeGPxTHIcZ76TzgzR1j1AW0+Gtp4M+6YpPjXJBMsaqli2qJJTW5ZwwqIqli+qZFlDFasaq1jeUDVuNNzmzZvZeP6Jxzh3zrFiugIlMxsC3iHpzcBdQGPRrHIcZ85hZhzsy0bRGaCtOzPm7smwv3vM3Zc3FDtHdTLBskWVLG+o4pSlS0aFKGyh9lOXKvcmOGeU6QrU53IOM7tJ0sPAnxTHJMdxjhUjI0ZXf5b23gwH0oO0pwdp782w5alBfnDw4VER2h+33LIO+dQkEyytr2RJXWp0ZoSl9SmW1qVG3S31ldRXuvg4R8a0BMrMPj/heAvwlqJY5DjOjDEz+gaHaU8PcqA3EwQnnaG9d5AD6cyoAIX9IB29gwxPbGeLNO7ZQ3NtiqX1KVY3L2ZpXYoldSmW1lfSEvdL61LU+BBsp0j4k+U4c5SREaN7IEtnX5bOvkEO9g3S2Ztzj9939mU52BcEJzf320RqU+U01SZpqkmyanE1L1zVEI9TNNUmaa5NjR4/dO/PuezSS45xjh1nPC5QjlNkhkeM9MAQXf1ZdnQN87OnD9DVnycwvWMCky86Xf3ZQwYR5CgTNFQnaaiuoLE6yYqGSs5YXk9jdUUUmhRNNckgONF9JFPtJMq8Kc4pPS5QjnMYzIyB7AjdA1m6+rN098f9QJauvizdUXy6c379Wbr7h0bj9AwMjU/wF78cd1hVkaCxuoKG6iSNNRUsa6iiMQpPQ3Uyzx32jdVJ6irLKXMRcRY4LlDOgmZkxOgdHCKdGSI9MERPbj8wRDqTjfvgl86E8J6BoVGx6Y5iU2hwQD7VyQSLqiqor6xgUVUFyxsqWb+sbvS4virsn3v6cS4670XUV42Jjk8i6jiFcYFy5hxmxuDwCH2ZYdKZIfoGw75nIDteTOL+ye0Zbt695RCxSQ8MkR4cGl2ldCpqkglqK8upTZVTW1lBfWU5KxqrxolOfVX5hOOwr6ssn/bS25t7nuaCk5uOsoQc5/jABco5KvLFpHdwiN64zx33DQ6RzgzTlxmid3Bs35sJYb2Z4bCPfr1RkIYm63yZQE0yQVIjNGXT1KbKqass54T6yuiuoLaynLpU+aj41FWGrTZVMSZIqXLvc3GcOYgL1HGAmZEZGqF/cJi+7DD9g8MMZIfpGxymPztM/+AQ/bnj3JYd5qntGb7X/hB9g8OHiElfZmi0djNdMYHQFFadLKcmlaAm7huqk6xoDP61qXKqkwlqcvtkOdWpRJ64VMRaTjk1ySAsmzdvZuPGXytiCTqOUwqKKlCSrgA+AySAG8zs48W83nzEzMgOGwNDQTRy4tA3OMxAnrt/Ylg2iEX/4Aj92aHRsEPjhOPpNHPlkygTyTKjrqONquT0xaQmFQSlJjkWVhP3VRUJ79h3HGfaFE2gJCWA64CXA7uAeyXdZmaPFuuas0GutjGQHWYgG/dDT0fFHAAACe9JREFUee7onxkaHh8nOzIqMgPZETKTnDeQHeZgug9+eseo/xFUQEaprCijOhk++vnuhuokyyoSVCcTVCYTVFckqEomqIx+VfG4qiKITFWyjKqKcqqS8ZyKEJYsL4s1k42zXsaO4zjToZg1qPOBp81sG4CkTcBVQFEE6vG93WzemWXbXdtHhSGTLwwTBWVo8vCZUpEQleUJUlE0KnP78iAGi2vKSFUkOHggQ+uqE8aFV1aMF5TxQpInMMkEleX/v72zj9WyrOP453sOHF4dL0pE6gLypVErXpxCmkkZL1aKG22gI3phbBUrqs1gbpZttVnNtGliK1sr81gGZKSiIMtZG8pBXkOEhBSGAk0liATk1x/X74H7PHDgnMPznHOdh99nu3eu+3df93V9n+e+7vN7ruu+7t8VPZEgCGofWVvHflpbsDQVmGRms3x/BnCFmc0pyzcbmA0wePDgMY2Nje2q74lth2ncfKiZrXsdNNRD9zrRUA8NddBQL7fLj5XZ6qB7PTTUtfK42xrqoa6Vccb2799P37592/U5O5LQWXm6itbQWVnOBp3jx49vMrPKLsNkZlXZgKmk506l/RnAPac6Z8yYMdZe3jp4yBY+vtzeOPC2HTx0xI4ePdrusqrNihUrOltCqwidlaeraA2dleVs0Amssgr7kWoO8e0ELizsX+C2Fmlqator6V9nUOd5wN4zOL+jCJ2VpavohK6jNXRWlrNB53srKQSqO8TXDXiJtNDhTuB54CYz21iVClOdq6zSXcwqEDorS1fRCV1Ha+isLKGzfVStB2VmRyTNAZaSppk/UE3nFARBENQWVX0PysweAx6rZh1BEARBbdK6AGJdh593toBWEjorS1fRCV1Ha+isLKGzHVTtGVQQBEEQnAm11oMKgiAIaoRwUEEQBEGW1IyDkjRJ0mZJWyXN64T6H5C0W9KGgm2gpKckbfG/A9wuST91reskjS6cM9Pzb5E0swo6L5S0QtI/JG2U9PUctUrqKek5SWtd5+1uHyZppet5WFKD23v4/lY/PrRQ1ny3b5Y0sZI6vfx6SS9IWpKrRq9ju6T1ktZIWuW2rK67l99f0iOSXpS0SdK43HRKutS/x9K2T9Lc3HR6+d/we2iDpIf83sqyjZ5Apd/87YyNNI39n8BwoAFYC4zoYA1XA6OBDQXbD4F5np4H3OHp64DHAQFjgZVuHwi87H8HeHpAhXUOAUZ7+hzSu2ojctPq9fX1dHdgpdf/e2Ca2xcAX/b0V4AFnp4GPOzpEd4eegDDvJ3UV/g7/SbwO2CJ72en0evZDpxXZsvqunsdvwZmeboB6J+jzoLeeuA10ouqWekEzge2Ab0KbfPzubbRE/RXu4KO2IBxwNLC/nxgfifoGEpzB7UZGOLpIcBmT98PTC/PB0wH7i/Ym+WrkuY/kSLOZ6sV6A2sBq4gveXerfy6k963G+fpbp5P5W2hmK9C2i4AlgMfB5Z4nVlpLJS7nRMdVFbXHehH+oeqnHWWaZsA/C1HnSQH9SrJAXbzNjox1zZavtXKEF/pIpTY4bbOZrCZ7fL0a8BgT7ekt0M/h3ffR5F6J9lp9aGzNcBu4CnSr7Y3zezISeo8psePvwWc2wE67wJuAUph8M/NUGMJA56U1KQUpBnyu+7DgD3Ar3zY9BeS+mSos8g04CFPZ6XTzHYCPwZeAXaR2lwT+bbRZtSKg8oeSz87spnTL6kv8EdgrpntKx7LRauZvWNmI0m9lMuB93eypGZI+jSw28yaOltLK7nKzEYDk4GvSrq6eDCT696NNFR+n5mNAg6QhsqOkYlOAPzZzfXAH8qP5aDTn4HdQHL87wH6AJM6U1NbqBUH1ebAtB3E65KGAPjf3W5vSW+HfA5J3UnO6UEzW5izVgAzexNYQRqK6K8U57G8zmN6/Hg/4N9V1nklcL2k7UAjaZjv7sw0HsN/TWNmu4FFJKef23XfAewws5W+/wjJYeWms8RkYLWZve77uem8FthmZnvM7DCwkNRus2yj5dSKg3oeuNhnpjSQutyPdrImSBpKs3Jmkp73lOyf85k9Y4G3fFhgKTBB0gD/5TPBbRVDkoBfApvM7M5ctUoaJKm/p3uRnpNtIjmqqS3oLOmfCjztv2AfBab57KRhwMXAc5XQaGbzzewCMxtKanNPm9nNOWksIamPpHNKadL12kBm193MXgNelXSpmz5BWuQ0K50FpnN8eK+kJyedrwBjJfX2e7/0fWbXRk9KtR9yddRGmiXzEuk5xa2dUP9DpDHew6RfgV8ijd0uB7YAy4CBnlfAva51PXBZoZwvAlt9+0IVdF5FGnZYB6zx7brctAIfAl5wnRuA29w+nHRjbCUNq/Rwe0/f3+rHhxfKutX1bwYmV+n6X8PxWXzZaXRNa33bWLpHcrvuXv5IYJVf+8Wk2W056uxD6l30K9hy1Hk78KLfR78hzcTLro2ebItQR0EQBEGW1MoQXxAEQVBjhIMKgiAIsiQcVBAEQZAl4aCCIAiCLAkHFQRBEGRJOKggCIIgS8JBBUEBSe/48gkblZb6+Jak094nku6StLM8r6Qpkm7z9HclmaSLCsfnuu2ydmgdJOmJtp4XBF2FcFBB0JyDZjbSzD5Ail4xGfjOqU5wp3QjKZjmx8oO3wL8rLC/nhR1osRnSS/Othkz2wPsknRle84PgtwJBxUELWApZt1sYI6HiWmJa0hO5j5S6BsAJF0CvG1mewt5F5OCdyLpfaRo0XsL5+yX9BPvwS2XNMjtF0la5r261X5uqbybz/SzBkGOhIMKglNgZi+TFqR71ymyleKxLQI+5cF4IQXlXF2Wdx8p1twH8QXhyo73AVZ5D+6vHO+9PQjca2YfBj5CCqsFKSTQR9v6uYKgKxAOKgjOAA9OfB2w2NKyJStJC8JBWpBuz0lOayQ5pykkp1bkKMed1m+BqzzI6/lmtgjAzP5nZv/1PLtJyygEQc3R7fRZguDsRdJw4B2OL5tQzkTSkuTrfRSwN3CQtHLpQdJyBeUsAX5E6intO/Xo4WnXE+rp9QRBzRE9qCBoAX/+swC4x1qOqjwdmGVmQy0tuzEM+KSk3qTlQS4qP8F7P98Gvn+S8uo4vgzCTcCzZvYfYIekKa6rh5cPcAkpSnUQ1BzhoIKgOb1K08xJyyU8SVqu4ATcSUwC/lKymdkB4FngM8AzwKiTTbAws0YzK38+BWkF2cslbSAtgPg9t88AviZpHfB34N1uH1+sPwhqiVhuIwiqiKS7gT+b2bJW5t9vZn3bUP4zwA1m9kZ7NQZBrkQPKgiqyw9Iz6Uqjg9B3hnOKahVogcVBK1A0kTgjjLzNjO7sTP0BMHZQDioIAiCIEtiiC8IgiDIknBQQRAEQZaEgwqCIAiyJBxUEARBkCX/Bx5UWM3HFQz9AAAAAElFTkSuQmCC\n",
      "text/plain": "<Figure size 432x288 with 3 Axes>"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "For most of the time we do not know the redshift directly and one has to infer it from the other paramerters. One  of the common parameters is the magnitude. Knowing the absoulte magnitude of a certain object will give relative magnitude of nearby objects by comparing their magnitudes. As a result, the distance modulus could be calculate from the magnitude."
  },
  {
   "metadata": {
    "scrolled": true,
    "trusted": false
   },
   "cell_type": "code",
   "source": "## provide unit for the distance modulus\nimport astropy.units as u\nfrom astropy.cosmology import z_at_value\n\n## suppose we have a distance modulus\nDval = 24*u.mag\n\n## infer redshift from the distance modulus\nzval = z_at_value(cosmo.distmod, Dval)\nprint(zval)\n\n## if no. of object is large (up to 10^6), then we use \n## interpolation instead\n\n## an array of uniformly random distmod between 24 and 59\nDvals_1 = (24 + np.random.rand(pow(10,6)) * 36) * u.mag\n\n## produce a relation from 24-30 using cosmo\nzmin = z_at_value(cosmo.distmod, Dvals_1.min())\nzmax = z_at_value(cosmo.distmod, Dvals_1.max())\n## generate uniform log space between zmin, zmax\nzgrid = np.logspace(np.log10(zmin), np.log10(zmax), 100)\nDgrid = cosmo.distmod(zgrid)\n\nzvals_1 = np.interp(Dvals_1.value, zgrid, Dgrid.value)\nprint(zvals_1)\n\n",
   "execution_count": 5,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "0.00014255181662361678\n[53.78794129 52.87823414 53.68137639 ... 53.79377029 52.67955284\n 52.46171127]\n"
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## 2. astropy.table\nOne can import table from other sources(sncosmo) and turn it into astropy table. The astropy.table allows user to manipulate the data in a way that numpy does."
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "import numpy as np\nfrom astropy.table import Table\n\na = [1,2,3]\nb = ['Happy', 'New', 'Year']\nc = [3.14, 2.71, 5.00001]\n\ntb = Table([a,b,c], names =('number', 'words', 'number again'))\nprint (tb)",
   "execution_count": 17,
   "outputs": [
    {
     "output_type": "stream",
     "text": "number words number again\n------ ----- ------------\n     1 Happy         3.14\n     2   New         2.71\n     3  Year      5.00001\n",
     "name": "stdout"
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df = tb.to_pandas()",
   "execution_count": 18,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df",
   "execution_count": 19,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 19,
     "data": {
      "text/plain": "   number  words  number again\n0       1  Happy       3.14000\n1       2    New       2.71000\n2       3   Year       5.00001",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>number</th>\n      <th>words</th>\n      <th>number again</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>1</td>\n      <td>Happy</td>\n      <td>3.14000</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>2</td>\n      <td>New</td>\n      <td>2.71000</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>3</td>\n      <td>Year</td>\n      <td>5.00001</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Apart from putting data in columns (which is default), user can put data in rows (which is more common in terms of observation), user can also define data type in each column:"
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "import astropy.units as u\n\nwavelength = []\nfor i in range(0,6):\n    m = np.random.rand()*100 + 900\n    m = m*u.AA\n    wavelength.append((m.value, m.to(u.nm).value, m.to(u.m).value))\n\nprint(wavelength)",
   "execution_count": 7,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "[(982.3229901460558, 98.23229901460559, 9.823229901460559e-08), (996.2501099798594, 99.62501099798595, 9.962501099798596e-08), (920.7264727547558, 92.07264727547559, 9.20726472754756e-08), (908.0093369228149, 90.80093369228149, 9.08009336922815e-08), (989.8302339492756, 98.98302339492756, 9.898302339492757e-08), (929.294410621016, 92.9294410621016, 9.292944106210162e-08)]\n"
    }
   ]
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "tb_1 = Table(rows = wavelength, names = ('angstrom','nm','mm'), dtype = ('f8','f4','i4'))\nprint(tb_1)",
   "execution_count": 8,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "     angstrom         nm     mm\n----------------- --------- ---\n982.3229901460558   98.2323   0\n996.2501099798594  99.62501   0\n920.7264727547558  92.07265   0\n908.0093369228149 90.800934   0\n989.8302339492756 98.983025   0\n 929.294410621016  92.92944   0\n"
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "user can obtain data from table as in numpy arrays"
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "tb_1[0]",
   "execution_count": 9,
   "outputs": [
    {
     "data": {
      "text/html": "<i>Row index=0</i>\n<table id=\"table139791717962192\">\n<thead><tr><th>angstrom</th><th>nm</th><th>mm</th></tr></thead>\n<thead><tr><th>float64</th><th>float32</th><th>int32</th></tr></thead>\n<tr><td>982.3229901460558</td><td>98.2323</td><td>0</td></tr>\n</table>",
      "text/plain": "<Row index=0>\n     angstrom        nm     mm \n     float64      float32 int32\n----------------- ------- -----\n982.3229901460558 98.2323     0"
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ]
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "tb_1[0][0]",
   "execution_count": 10,
   "outputs": [
    {
     "data": {
      "text/plain": "982.3229901460558"
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ]
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "tb_1['angstrom']",
   "execution_count": 11,
   "outputs": [
    {
     "data": {
      "text/html": "&lt;Column name=&apos;angstrom&apos; dtype=&apos;float64&apos; length=6&gt;\n<table>\n<tr><td>982.3229901460558</td></tr>\n<tr><td>996.2501099798594</td></tr>\n<tr><td>920.7264727547558</td></tr>\n<tr><td>908.0093369228149</td></tr>\n<tr><td>989.8302339492756</td></tr>\n<tr><td>929.294410621016</td></tr>\n</table>",
      "text/plain": "<Column name='angstrom' dtype='float64' length=6>\n982.3229901460558\n996.2501099798594\n920.7264727547558\n908.0093369228149\n989.8302339492756\n 929.294410621016"
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ]
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "tb_2 = tb_1[0:2]\nprint(tb_2)",
   "execution_count": 12,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "     angstrom        nm     mm\n----------------- -------- ---\n982.3229901460558  98.2323   0\n996.2501099798594 99.62501   0\n"
    }
   ]
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "tb_3 = tb_2['angstrom', 'mm']\nprint(tb_3)",
   "execution_count": 13,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "     angstrom      mm\n----------------- ---\n982.3229901460558   0\n996.2501099798594   0\n"
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "The QTable accepts quantities while creating the table, and their units are taken account in numeriacal operations"
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "from astropy.table import QTable\nfrom astropy.coordinates import SkyCoord\nfrom astropy.time import Time\nimport astropy.units as u\n\ntb_4 = QTable()\ntb_4['time'] = Time(['2000:002', '2002:345'], format = 'yday', scale = 'utc')\ntb_4['SkyCoord'] = SkyCoord([10,20],[-45, 40], unit = 'deg')\ntb_4['velocity'] = [300, 400]*u.km/u.s\ntb_4",
   "execution_count": 17,
   "outputs": [
    {
     "data": {
      "text/html": "<i>QTable length=2</i>\n<table id=\"table139791712338440\" class=\"table-striped table-bordered table-condensed\">\n<thead><tr><th>time</th><th>SkyCoord</th><th>velocity</th></tr></thead>\n<thead><tr><th></th><th>deg,deg</th><th>km / s</th></tr></thead>\n<thead><tr><th>object</th><th>object</th><th>float64</th></tr></thead>\n<tr><td>2000:002:00:00:00.000</td><td>10.0,-45.0</td><td>300.0</td></tr>\n<tr><td>2002:345:00:00:00.000</td><td>20.0,40.0</td><td>400.0</td></tr>\n</table>",
      "text/plain": "<QTable length=2>\n         time          SkyCoord  velocity\n                       deg,deg    km / s \n        object          object   float64 \n--------------------- ---------- --------\n2000:002:00:00:00.000 10.0,-45.0    300.0\n2002:345:00:00:00.000  20.0,40.0    400.0"
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## 3. astropy.SkyCoord\nThe astropy.SkyCoord allows user to handle the coordinates of objects(or events) on sky, and support transformation between coordinate system. Moreover, Skycoord provides distance calculation between two objects:"
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "from astropy import units as u\nfrom astropy.coordinates import SkyCoord\n\n#create coord objects in different ways\ncoord_1 = SkyCoord(10.625*u.degree, 41.2*u.degree, frame = 'icrs')\ncoord_2 = SkyCoord('00 42 30 +41 12 00', unit=(u.hourangle, u.deg))\n\n#SkyCoord also accept array input\nra = [10, 11, 12, 13]*u.degree\ndec = [41, -5, 42, 0]*u.degree\ncoord_3 = SkyCoord(ra, dec, frame = 'icrs')\n\nprint(coord_1)\nprint(coord_2)\nprint(coord_3)",
   "execution_count": 19,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "<SkyCoord (ICRS): (ra, dec) in deg\n    (10.625, 41.2)>\n<SkyCoord (ICRS): (ra, dec) in deg\n    (10.625, 41.2)>\n<SkyCoord (ICRS): (ra, dec) in deg\n    [(10., 41.), (11., -5.), (12., 42.), (13.,  0.)]>\n"
    }
   ]
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "#access the data in coordinate objects as np.array\nprint(coord_1.ra)\nprint(coord_2.ra.hour)\nprint(coord_3[0].ra.hms)\nprint(coord_3[2].dec.radian)",
   "execution_count": 21,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "10d37m30s\n0.7083333333333334\nhms_tuple(h=0.0, m=40.0, s=4.263256414560601e-13)\n0.7330382858376184\n"
    }
   ]
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "#the coordinate are default in observer = origin, but one\n#could change to heliocentric or galactocentic\ncoord_4 = SkyCoord(10.68458, 41.46917, 10, unit = (u.deg, u.deg, u.kpc))\nprint(coord_4)\nprint(coord_4.galactic)\n\nfrom astropy.coordinates import Galactocentric as galcen\nprint(coord_4.transform_to(galcen))",
   "execution_count": 26,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "<SkyCoord (ICRS): (ra, dec, distance) in (deg, deg, kpc)\n    (10.68458, 41.46917, 10.)>\n<SkyCoord (Galactic): (l, b, distance) in (deg, deg, kpc)\n    (121.18204208, -21.37301734, 10.)>\n<SkyCoord (Galactocentric: galcen_coord=<ICRS Coordinate: (ra, dec) in deg\n    (266.4051, -28.936175)>, galcen_distance=8.3 kpc, galcen_v_sun=(11.1, 232.24, 7.25) km / s, z_sun=27.0 pc, roll=0.0 deg): (x, y, z) in kpc\n    (-13.13328657, 7.96690902, -3.60167428)>\n"
    }
   ]
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "#most of the time we perform calculation in cartensian\n#coordinates, for a ra, dec, dist input\n\ncoord_4.representation_type = 'cartesian'\ncoord_4\n\n#in case you would like to have the distnce between two objects\n",
   "execution_count": 31,
   "outputs": [
    {
     "data": {
      "text/plain": "<SkyCoord (ICRS): (x, y, z) in kpc\n    (7.36321131, 1.38924091, 6.6221695)>"
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ]
  },
  {
   "metadata": {
    "trusted": false
   },
   "cell_type": "code",
   "source": "c = []\nra = [10, 11, 12, 13]*u.degree\ndec = [11, 12, 13, 14]*u.degree\ndist = [5, 6, 7, 8]*u.kpc\n\nfor i in range(0,4):\n    c.append(SkyCoord(ra[i], dec[i], dist[i]))\n    \nprint(c)\nprint(c[0].separation(c[1]),\nc[1].separation_3d(c[2]))",
   "execution_count": 36,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "[<SkyCoord (ICRS): (ra, dec, distance) in (deg, deg, kpc)\n    (10., 11., 5.)>, <SkyCoord (ICRS): (ra, dec, distance) in (deg, deg, kpc)\n    (11., 12., 6.)>, <SkyCoord (ICRS): (ra, dec, distance) in (deg, deg, kpc)\n    (12., 13., 7.)>, <SkyCoord (ICRS): (ra, dec, distance) in (deg, deg, kpc)\n    (13., 14., 8.)>]\n1d24m00.2881s 1.0124163712587475 kpc\n"
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# 4. astropy.time"
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python",
   "version": "3.7.1",
   "mimetype": "text/x-python",
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "pygments_lexer": "ipython3",
   "nbconvert_exporter": "python",
   "file_extension": ".py"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}