{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Synthetic photometry" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this tutorial, we will compute the synthetic [MKO J band](http://svo2.cab.inta-csic.es/svo/theory/fps/index.php?id=MKO/NSFCam.J&&mode=browse&gname=MKO&gname2=NSFCam#filter) flux from an [IRTF spectrum of Jupiter](http://irtfweb.ifa.hawaii.edu/~spex/IRTF_Spectral_Library/References_files/Planets.html)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Getting started" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start by importing the required Python packages." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import species" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The [species](https://species.readthedocs.io/en/latest/species.html) HDF5 database is initiated by creating an instance of the [SpeciesInit](https://species.readthedocs.io/en/latest/species.core.html#species.core.init.SpeciesInit) class." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Initiating species v0.5.2... [DONE]\n", "Creating species_config.ini... [DONE]\n", "Database: /Users/tomasstolker/applications/species/docs/tutorials/species_database.hdf5\n", "Data folder: /Users/tomasstolker/applications/species/docs/tutorials/data\n", "Working folder: /Users/tomasstolker/applications/species/docs/tutorials\n", "Grid interpolation method: linear\n", "Creating species_database.hdf5... [DONE]\n", "Creating data folder... [DONE]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "species.SpeciesInit()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Jupiter spectrum" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The spectrum of Jupiter that is used as an example is now downloaded from the [IRTF website](http://irtfweb.ifa.hawaii.edu/~spex/IRTF_Spectral_Library/References_files/Planets.html)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "('data/plnt_Jupiter.txt', )" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import urllib.request\n", "urllib.request.urlretrieve('http://irtfweb.ifa.hawaii.edu/~spex/IRTF_Spectral_Library/Data/plnt_Jupiter.txt',\n", " 'data/plnt_Jupiter.txt')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The file contains the wavelength in $\\mu$m, and the flux and uncertainty in W m$^{-2}$ $\\mu$m$^-1$, which are also the units that are required by [species](https://species.readthedocs.io/en/latest/species.html). We can read the data with [numpy.loadtxt](https://numpy.org/doc/stable/reference/generated/numpy.loadtxt.html)." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "wavelength, flux, error = np.loadtxt('data/plnt_Jupiter.txt', unpack=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's create a [SpectrumBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.SpectrumBox) with the data." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "spec_box = species.create_box('spectrum',\n", " spectrum='irtf',\n", " wavelength=wavelength,\n", " flux=flux,\n", " error=error,\n", " name='jupiter')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And pass the [Box](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.Box) to the [plot_spectrum](https://species.readthedocs.io/en/latest/species.plot.html#species.plot.plot_spectrum.plot_spectrum) function together with the filter name." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Adding filter: MKO/NSFCam.J... [DONE]\n", "Plotting spectrum... [DONE]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAEJCAYAAAByjOZgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABSP0lEQVR4nO3deXhU5dn48e89W/aNBAirrLIrCCioKC644FKx1pWqrUu1tVbb/vpaW0ustVrf12prW7dq7YJb1WpxR1lUZJdFQED2ELZAyL7N8vz+mJMwCQnJTCY5k8n9ua5cZs6cec49J5i58yz3I8YYlFJKKaXilcPuAJRSSiml2pMmO0oppZSKa5rsKKWUUiquabKjlFJKqbimyY5SSiml4prL7gCiQUTcwGjrYSEQsDEcpZRSSkWHA+hufb/OGOONpJFOk+yIyETgJeA3xpgXGj09Gviiw4NSSimlVEc5CVgVyQs7RbIjIjOAbwEldseilFJKqc6lUyQ7wHJjzH9EZEEzzxfWfbNy5Upyc3Prn+jTpw8FBQVRCaIrtBXt9rQte9vTtuxtT9uytz1ty972otHWvn37GD9+fN3DwmOdeyzSmSooW8nOC42HsUSkN1AAUFBQQO/evUOfI1rvsSu0Fe32tC1729O27G1P27K3PW3L3vai0daePXvo06dP3cM+xpg9kbQT96uxZs2aZXcITYpmXNF+j7EcW7TE8nuM5diiJZbfYyzHFi2x/B5jObZoieX3GMuxtUXc9ew0NmvWLPLy8qJ1/ahm0Co8ev/tpfffXnr/7aX3v+Pk5eVx//33N/WU9uzUKSgowBhT/xWtRAdiK0vtivT+20vvv730/ttL73/HycvLq/8Mj9b8objr2Wk8Z0cp1bkZY9i0aRPr1q1jzZo1fPXVV5SWltKvXz/69OnDgAEDuPjii+nRo4fdoSqloixac3Y6y2ospVQXYoxhxYoVvPbaa6xcuZIhQ4YwYcIELr74Yn784x/jcrkoKSnhwIEDrFu3jttvv53KykrGjx/PzJkzGT58eETXXbduHf369SMjIyPK70gpZadO0bMjIuOBR4GxwD5ggzHm8pDntWdHqQ5gjGH+/PksXryYiooKxowZw4wZM0hMTGzV6wOBAMYYamtrSUpKavBcRUUFn332GXPmzGHjxo1MmDCBb33rW5x00kmISItte71e5syZw4cffsjWrVtJSEhg4sSJ5OTkMGTIEAYOHMjQoUObbCsQCPCzn/2MHTt2UF1dzVtvvYXT6WzdTVFKtZto9ex0imSnJZrsKNX+Fi9ezB133MEZZ5zB5ZdfTlpaGsuWLeOll15i5syZ3HjjjU0mCOvXr+eVV15h3rx5LFmyhDFjxpCamorH48HpdFJdXQ1At27dOP3007n44osZNmxYqxKcY6mtrWXp0qWUlpayceNGdu7cyaZNm+jfvz+XXXYZZ511FsnJyWzdupX/+Z//Yfr06Xz3u9/lgQceYNKkSUybNq1N11dKtZ0mOyE02VGq9R566CGKior41a9+RVpaWoPntm/fzpo1awgEArjdbhwOB5s2bWL+/Pn079+fWbNmHTU3xuv18uc//5m33nqLadOmMX78eMrKyvj8889Zvnw5AwcO5KabbmLKlCk4HPavidiyZQuvv/46n376KT6fj+zsbO655x7GjBkDBIeynnvuOR577DGbI1VKabITQpMdpVrnww8/5IMPPmD69Ok89NBDPP300/Tq1Yv33nuP2bNnk5GRwahRo+jWrRterxefz0evXr2YNm3aUYlRY2VlZXz22Wds2rSJQCDAKaecwsSJE/F4PB307qLDGMP06dN577337A5FqS5PJygrpcL25ptv8pOf/ITBgwczcOBA7r//fsrKypgyZQpPPvkkPXv2jLjttLQ0LrzwQi688MIoRtzxRASn04nf79d5O0rFCU12lOpCtm3bxqBBgwAYNGgQf//7322OKDYNHDiQHTt2MHjwYLtDUUpFgf0D6EqpDlFRUUFycnKbJ/52BaNGjWL9+vV2h6GUihJNdpTqIr766itGjx5tdxidgiY7SsWXqCY7IuIWkRNE5LpotquUarstW7YwdOhQu8PoFEaOHKnJjlJxpE1zdkTkROAiYAwwGvAAG4BVbQ9NKRVNW7Zs4ZxzzrE7jE4hOzuboqIiu8NQSkVJW3t2XgOKCSY4h4GLjDEzjDG/bmtgSqno2rJlC0OGDLE7jE5Dd7lWKn60dTXWHcaYDwBE5FXgSRF5D3jM6G8JpWJKYWEhOTk5dofRaeTk5HDw4EG6d+9udyhKqTZqdc+OiPzK+qr/07Au0bG+3wScB3iBD6MapVIqKnQlVusNHTqULVu22B2GUioKwhnGmgEsAJodyDZBTwA3tzEupVQUlZWVkZqaancYncrQoUP5+uuv7Q5DKRUF4SQ7xcaYT4wxzSY7IjIUwBizs82RKaWiZuvWrTpfJ0xDhgzRZEepOBFOstOaOThPRxqIUqr96OTk8OkwllLxI5wJylNFxN9ukSil2s2WLVs47bTT7A6jU0lPT6e0tNTuMJRSURBOsrMGuOsYzwvwWJuiUUq1iy1btnDDDTfYHUanZIzRid1KdXLhJDuHjTELj3WCiDzetnCUUu1h79695Obm2h1Gp9O9e3ddfq5UHIjqdhHGGN1CWakYpb0T4dMVWUrFh3CSnV4iMk9ETmm3aJRSUVe327kK35AhQ3SSslJxoNXJjjFmhDHmbGPM0vYMSCkVXVu2bGHw4MF2h9Epac+OUvEhqsNYSqnYs379ekaPHm13GJ3SkCFD2Lx5s91hKKXaSJMdpeLcunXrNNmJkC4/Vyo+aLKjVJz76quvGD58uN1hdFpZWVkcOnTI7jCUUm2gyY5SccwYQ3V1NYmJiXaH0mlNmDCBlStX2h2GUqoNNNlRKo7t2rWL4447zu4wOrXJkyezaNEiu8NQSrWBJjtKxbFPPvmEU0891e4wOrWJEyeybNkyu8NQSrWBJjtKxbG3336biy66yO4wOjWXy0VGRgaFhYV2h6KUipAmO0rFqaqqKioqKsjOzrY7lE7viiuu4NVXX7U7DKVUhDpNsiMiiSLygogsEZEVInJeR8eQl5fX0ZdUIfT+h+fll1/m4osvjlp7Xfn+X3LJJbz++uv4/X7bYujK9z8W6P3v3MQYY3cMrSIiDwO9jDE3iMjxwBJghDFmv4j0BgoACgoK6N27d3vFQGe5X/FI73/rFRcXc/nll/PBBx/gdruj0mZXv//PPPMMVVVV/OhHP7Ll+l39/ttN77899uzZQ58+feoe9jHG7ImknU7RsyMiDuBm4DkAY8xmYBUws6XXxmo2Hs24ov0eYzm2aInl99jW9nw+H3fffTf33XcfDz74YHSCirLOeP9vvvlmPvnkE+bOndvmtuzWGe+/3W1FUyy/x1iOrS06Rc+OiAwBvgZ6GmMOWMeeAHKMMdeIyABgO8DKlSvJzc2tf22fPn0oKCiIShxdoa1ot6dtdWx7+fn53HPPPcyYMYMrrrgiZt9nrLbVUnuVlZXcd999JCYmcvfdd5OTk9NhsXWFtqLdnrZlb3vRaGvfvn2MHz++7uFAY8yOSNrpLMnOacBnQKIxpsY69gBwqjHmHBG5GnjJzhiVUkop1a6uMca8HMkLO8UwVojGmZnYEoVSSimlOg2X3QG00gHrv5nNfL+k7sR33nmnwTDW+PHjo1bqvSu0Fe32tC1729O27G1P27K3PW3L3vai0da+fftCa4UtOda5x9JZhrEcQCFwmTHmU+vYx8C7xphHj7UaK5oz6LtCW9FuT9uytz1ty972tC1729O27G0vGm11qdVYxpgA8CzwXQARGQqMBWa39NpZs2a1a2yRimZc0X6PsRxbtMTye4zl2KIllt9jLMcWLbH8HmM5tmiJ5fcYy7G1Rafo2YFgUUHgKWA4weG3e40xH1rPaZ2dLkDvv730/ttL77+99P7bI1o9O51lzg7GmGrgRjtjiKUstSvS+28vvf/20vtvL73/nVun6dk5lo7q2VFKKaVUx+lSc3aUUkoppSKlyY5SSiml4lrcJTt9+vRBRGJqTw6llFJKhScvLy90CKtNdM6OUkoppWKSztlRSimllGoFTXaUUkopFdc02VFKKaVUXNNkRymllFJxTZMdpZRSSsU1TXaUUkopFdc02VFKKaVUXIu7ZCcWigr+ad7XfLK5kPIaHxU1PtviUEoppTorLSrYSKwVFfzFf75k9tJd9Y8nDerGX64bT7cUj41RKaWUUp2LFhWMUb999yu2FpY3OLZkWxF5/11vU0RKKaVU1+ay8+IiMhF4CfiNMeaFZs6ZCjwF7As5/CdjzGvtHV+4qmr9PPPJNgCmDutOr4xEPtl8kIE5KXyx67DN0SmllFJdk23JjojMAL4FlLTi9IebS4ZiyYhfvV///UOXj6HaG+DBy8ZQVuPjL/O32BiZUkop1XXZ2bOz3BjzHxFZYGMMUXP8L99r8LhXRlL99xlJbhLcTowxiEhHh6aUUkp1abbN2THG7A7j9MtEZL6IfCoivxQRW4ffGluy7RC1vgBpicGwfnPZ6KPOSXA5qPEFOjo0pZRSqsvrDBOUS4DFwDTgQuBc4He2RhTirdUFPLlgK4O6p3DygG4AzJx03FHnZSS5Ka70dnR4SimlVJcXUbIjIlNE5PsiMktEfiQi00XEHe3gAIwxq4wxvzPG+Iwx5QQTndukmfGgujo7dV/tXW/nmU+2sXBzIVW1fo41QjW0R+pRq7SUUkop1VBeXl79Z3i06uyENRwkIpOB54FyIN/6bxLwTeApEfm5MWZ2VCJr3i4gGegOHGj8ZEfX2Tnj+O789YYJZCZ5WLnzMBU1/ibP652ZxJRH5rPj4Ys6LDallFKqs8nLy6vvqGhUZydire7ZEZGhwK3AmcaYicaYy40x1xtjvmWMOQMYDgwTkUvaHFXD694pIokhh3oCtcChaF4nEgfLa3ht5W56ZSSR5HFy+tAchuWmNXlu36wkrhjft4MjVEoppVQ4w1jFwE3GmKN6UwCMMZXGmF8By9oSkIjkWBORu1mHTgKutJ5zAj8EXjTGNN2F0oFyUhP4102nNDh22pCcJs8VEXpnJBIPFauVUkqpzqTVyY4xptAYU7+cSERGNnPe/ta0JyLjrWXnY4F7ROQN66kkgr1EydbjZ4GrRWQ+sAQoBH7U2rjbW+OenGkjezZ77h/nbeHfK8NZhKaUUkqptmrLEu7HgfMifbExZiUwtYnj+QTn49Q9XgRMj/Q6seZQea3dISillFJdSovJjoiUEuxREaBuDEYI9sioMLz5g9PYvK/M7jCUUkqpLqU1PTubgSuNMcWhB0VkbrtEFMf6ZiWxdnex3WEopZRSXUpr5uxMA47qjjDGTIt+OPEt2eOkstb2edVKKaVUl9JismOMORwLK59aq66oYHsXE4xEosvJI+9vtDsMpZRSKubl5eVFrahg2BWUReRnUblyOykoKMAYE5PJjsMhBHTluVJKKdWivLw8CgoKotJWJNtFXBCVK3dRt505uEGtnapaP08v3EphWY2NUSmllFLxK5Jk5xg7QKmWPLVwKxv2ltY/PlRRw0PvbWTigx/ZGJVSSikVvyJJdnQgpo3qOnbKqr1c/1ybCk4rpZRSqgXas2ODy//yOQDnPLqQtKQjm8X7dUKPUkopFXWRJDvfb+qgtVGoaoVaf6D+v+P6ZZLgcnDLlIHc8eIXNkemlFJKxZ+wt4swxnwlIg7gLKAXRxKmmbRh+4iu4vieqRyu9DJv435SPC7Kqn3cPe14Asbw7Kfb7Q5PKaWUijuR9OwA/BfIA84lmPScBURnMXwbxXKdHYAP7jqD7BQP331hBSkJTipqfNx25mBcDh0dVEoppepEs85OpBuBZhpjTg89ICIXRiGeNisoKKB37952h9EsEWGjtT9WssdFRa0PAIdosqOUUkrVycvL49Zbb41KwhNpsrNQRAYbY7aGHBvS5mi6iOG5afTrlkxVrb9++whNdpRSSqn2Eekw1nJglYgUiMg2EdkOPBDFuOKaL2DYV1JNsic4jAXgtIaxPv260M7QlFJKqbgTac/OQ8BlwDaCdXeE4Bwe1Qr+gGFATgpuh1BuJTsOK9nZtK+MKUO72xmeUkopFVciTXa+NMbMCz0gItqz00pefwC3Q0hOcB5JdgTOHdFDh7OUUkqpKIs02dklIn8DPgfqNnXSpeet5PMbXE4hI8nNk9eNB8ApwvDcdEqqvDZHp5RSSsWXSJOda4APgVNDjsXE0vPOwBcI4HI6cDkcTB6cDQSHsTwuBz6toqyUUkpFVaTJzgPGmKdCD4jIRVGIp0s4WF6Lu1FdHYcILqdosqOUUkpFWUSrsRonOpbcNsYSFbFeVLCOy+losKOq0wEep0M3HlNKKaWIgaKCIjKvicNDgefaFk7bxXpRQYAEl4MEl4NASC+OQwR3owRIKaWU6qpioahgMfAH63s3MBZIbXM0XYQ/YEhyOym3qifDkWEspZRSSkVXpMnOTGNMZcjjj0TkcQAR6W6M0cp4x3DTlIEkup0Uh6y8cjqCPTtKKaWUiq5Ik50JcqQejIPg7udnisgZwI8JFhxUzfj5hSO448Uv2F9aXX8sOIwlOmdHKaWUirJIk50Xgc3Q4LO5GLif4Nwd1YJ5Gw/U74sFwaKC2rOjlFJKRV+kyc7PjTH/bOoJEbm9DfF0GWcP78Hba/fWP3Y6BJfpXBOU95VUk5uRaHcYSiml1DFFuvS8yUTHeu7JyMPpOh6/aiwf/+TM+sfBooLBYaxAJ6m1M+mhj/l6fxkA1z+/zOZolFJKqabpuIlNXE4Hg7sfWcDmEMHlcJDgdlDrD9gYWcvu/c+XfLh+HwDTHvuEe//zJZ9sLmTx1kM2R6aUUkodLe6Snc5SVLAxt7VdhNsRm1tGrNxZxFVPLwbgxaW7uPWfK+ufKywLbo+WX1TZ5GuVUkqpcNleVDCWdYaigk05ZVBwj6wNe0rxxWDPzpw1e1m6vYhzf78QgES3g2pvgIwkN59vOYhD0DpBSimloiaaRQVb1bMjIoki4ml0rFebr67qOR1i1doRzn50oS0xHGuukMMqNbDlQDkA1d5gQuZ2Oqio9TNjXF8Oldfy4tJd7R+oUkopFYYWkx0RuRt4F3hbRJ4UkSTrqdntGlkX5XQ4KKqoteXa0x5bSFXIcvgl2w6xOr8YAGnUaZObHlyF5fUH+J8LhpPgdvDgu19x73++7DQTrJVSSnUNrenZ+aYx5mxjzHnAG8AcEYlKLR0RmSgiW0TkxhbOu05EVorIChF5VKTxR2/8qBsK8towlLWnuJrTfjevfu5NUUUte4qr+HzrQZ77bHuDcy85sRc/nnY8Myf1RwRSE1x857QBADzywaaODl0ppZRqVmuSHaeIuACMMXOBG4G/AMe35cIiMgO4Gyhp4bzRwKPA+cDJwEnA99ty7VjmcgSTnfJqXwtnRp/XH6CoopZXV+QDUFnr51B5DaUh21qsv/98Nv/mQgIGkj1O/t/5w+uf8wcMGUluhvbQbdKUUkrFjtYkOz8FcuoeGGN2AxcBv2zjtZcbY64Fylo47ybgXWPMQWNMAHgeuK2N145ZdXNjyms6PtmpWwX2xLwt+PwBKmt9HKqo5bZ/fVF/TkqCC4/LwYCcFI7vmdbg9bW+AFdO6MvfPt/OrkO6MksppVRsaDHZMcYsMsbsa3Ss1hjzQlsubCVNrTER2BjyeAMwKmTuUFxxOASP00GZDT07OakJ9d/nH65i0ZaDPP7R1/XHQosgfnvScZxxfHcAarwBElwOan0B0hLdrCso5Y1Vrf3xKqWUUu0r7Do7IvKz9gjkGHrScKirmOCeXDlNnVxXZ6fuq7PV23EIJCc4Kav2tnxylF13Sn82/eYCAB6bu5kP1u9v8HxoEcRQyR4nmckeavwB0hKD1Qwe/+jr+pVbSimlVGvl5eXVf4ZHq85OJEUFL4jKlcPT1PKeJicpFxQUYIyp/+p8yY6QmeTmcGXzK7ICAcPKnYejet2rnl6MARJcTs4d0YP/rtkDQLcUz7FfCNw8ZSDfOXUAtb4AqQku0hKCCY8mO0oppcKVl5dX/xleUFAQlTYjSXY6eiXUASAz5HEmweSnsIPj6BAOETKTPRwsbz7ZqfT6mfnXpVG97tLtRfzNWnE1Y1zf+uNJbidz7z6DR791YrOvFREcDrGGsVz1213c9q+Vzb5GKaWU6iiRVFDu6CIqy4FhIY9HAuuNMVUdHEeHcAhkJruPWWunssZHSkL0il+/tjI4v+aOs4fUHxvUPQWf3/Dh3WeQ6HYytNFk5KYEe3bc1PhirwK0UkqprivmenZEJEdEPhWRbtahvwLTRSRbRBwEl74/1Z4x2MkhgkPkmPtj/WXBVhxR/Cm8vy44/3xM3wwgWEBwTJ8MuqV4SHQ7W92OP2AabBnxo3OiUo5JKaWUapNIkp0ma9yEW2hQRMaLyAJgLHCPiLxhPZUEDAeSAYwx6wguf/8QWAqsJljnJy45HBAw5pgZ5Quf7+BAWQ1Lt7V9l/HFWw/hDwR7YlKt3qLzR+Xy8OUnkN2K+TqhAsbgEOEPV49tc1xKKaVUtIQ9FmKM+crqYTkL6MWRhGkmcF4Y7awEpjZxPB/o3ujYbLrI9hQOEVrabeGiE3rxztq9XPXMEnY8fFHE11q7u5if/nsNAWPITU8k2RP85+B0CEkeJ9mp4SU7Lmdwf69vjO1DwBj+sXgn157Sn57W1hJKKaWUHSLp2QH4L5AHnEsw6TkLiM76sC7OIYIxx852BndPZeMDF9AnM6nFc4/l0j8toqC4Cq8/2EbjlVf9spLDau9vN57MSf0zgeAk51W7itlxsCLi+JRSSqloiHSWa6Yx5vTQAyJyYRTi6fKCPTstJzCJbienD8lh4eZCpg7r0aZr+gMBFv6/s4+a9Bw6Ybk1kjxHz+/xuCLNp5VSSqnoiPSTaKGIDG50LLxPxnZSV1Sws9XXqeOQ4ETf1rhuUv9j1uNpSd0E4oChydVdbd1v9cZTB1Dt1ZVZSimlwpeXl2drUUEILgdfJSIFIrJNRLYDD0QlojaqKyrYaZMdR8tzdupSkB5piVTVRp5M1F3msrG9I27jWC4d25tqn79d2lZKKRXf8vLyolZUMNJhrIeAy4BtBD8zheAcHhUNpnXFjBLdDqq9kSUTFTU+th+s4PtTB/OzC4a3/IIIJLqcHPBWt0vbSimlVGtFmux8aYyZF3pARGKiZ6ezM3WpYyvn7ZRURbaH1sZ9pcxZs6dNq7lakuRxUhVhMqaUUkpFS6TJzi4R+RvwOVBjHQtr6blqmiFYY6ei1s/hilqymqh1U5cGJbgcEc/Z8ThbXywwUoluB4VlNS2fqJRSSrWjSOfsXGP991R06Xl0mWAF40mDstl9+OgdMfaWVPHSsl1AcAJx3ffhenvtHp6/cUKbQm1JosvJb9/d2K7XUEoppVoSabLzgDHmO6FfwM+iGVhXFRzFEhJcDi7502dHPV9YVsM9IXNsfnTOUGoj2Ivq6U+2EWjnhVLhbDWhlFJKtZeIkh1jzFF7Uxlj3ml7OMpYPTuhe0yFOlzpZUBOSv3jRHdk82Kmj8nl5EHdWj6xDRK0xo5SSqkYoJ9GMcZgEIHJg7K5fWrjUkZww/PL6JGWUP84weXA6w+/i2ZIjzTSE91tirUljmjuVqqUUkpFKO6Snc5eVNCY4DCWiOBxNv3j6ZuVVP+92xlZstNR7tSdz5VSSkXAtqKCItI+BVmiqLMXFYTgMFZzfnTO0AaVjd1OR0Rzdv748deRhBaRtuzfpZRSqmuKZlHBcHt2nhGRSdLWfQRUs8JNCzwRDmN1VI+L2yH4Wrn9hVJKKdUewk129gBDgT+KyBMicoOI9GyHuLqsgDH1PTer84sbPLenuIqs5IbzbII9O+EnEx2VrSa6nRFXebbTDc8vwxjTYJ8yXwwPFyqllGpeuEUFrzPG+IF/Wr07E4FbReRUYJUx5t6oR9jFTDgui76ZwTk5J/bLbPDc3pIqBnZPbXDM4xJqw/wQNsYcc6gsmoJbWgRIS+yY60XLql2H+e4Ly+lt/SyuPaU/t/1rJc/fMJEhPVLbvEmqUkqpjhNWsmMlOohIOjANmA6cD3iAwqhH1wWlJbpJs1ZJBXeNONLTs3BTIWcN79HgfI/TGfYwltdvcDcz+TnaEjphz866ghJKq33M33Tkn/TspbtwCFzyp8+Y95Op9UmQUkqp2BfuBOV7ReQTYB9wF7Cd4IagPY0x10c9ui4uOAR0JJGprPUzsnd6g3PcTgl7gnKtP4C7mTo+0ZbodlLTiXY+9/oDXPzEkWKOA3NSWHrvOcyc1J/huelUewPc8eIXNkaolFIqXOH+eT8G6AXcB1xpjPmNMWaF0eU27SLJ7WhQMDA5wUWCq2FVYncEE5S9vkCzy9qjLdHlaJCwxaqSSi9XPPk54349F4BHvnkCLodwQt8MeqYnkpOawC8vHkH3tAROH9pdV5gppVQnEtYnnjHmGmAYsIjgXJ2nReRBETlDRGJib4DOXmcnVLLH1WJ1ZLfDgdcf3gdvWbUPV0clO51kGOtP879mxc7DlNf4+PmFw+nbLYnxx2Xx2JVjAbhgdC7DeqZRWePj6YVb+ePHW+pf+8/FO+wJWiml4phtdXYAjDEBY8wSY8z9BHt4DgFvWv+1XTzU2amT6HFSVes75jlul4S9SujOl1exrbCiLaG1WuOhuFj17Kfbufn0gSz8f1O59YxBeJwOav2B+irQw3PTyUz20CM9kRpfgHV7Supfu2LnYbvCVkqpuGVbnR0ReVJEThOR34jICmAXcCnwCHBmVCJS9ZLcTqpqjyQKTc2ycTkcYa/GOmNoDj+f3jH1Icuqvfz8P2s75Fpt8cOzh/CLi0ZwXHYKIsLQnmncf+moBuc4HcK8n5zJrEtGMnfDfq5/fhlj8j7grdV7bIpaKaVUa4Tbs/M94GUgF3gI6GGMmWqMedgYsybq0XVxyR4nlS307HicDnxhDmMh0mGrsU4ZlM0ZQ7t3yLUi9ffPd/CvJTsbLCfPSHJzQt/Mo84Vkfqk85PNhZRVB38+OodHKaViV7ifeE8bY/oZY242xrxujCltl6gU0Lodzd0uiem9sVITXLzxRQGvrsi3O5RmHa6spa1FnncfrsLrDzDjL4uiE5RSSqmoCbfOzu3tFYg6WrLHyYHS6vrHTX0euxyxvREoQJXXT0ml1+4wmmUMLL33nFaf37ig4OXj+jDlkfmM7pPO5n3lGGMwRnd9V0qpWBF3u57Hk9QEF6XVx04SPM7wV2PZIVYLDh8sr2HexgMkulu/mFAErjm5P/+86WQARvfJAGBdQSm1/gD3z9nAkwu3tku8SimlwqfJTgzrk5nEnuLqY57jcoY/jGVH3uGM0V6OrQfK+bKgpOUTQwjgEJgytDvPXj+BG04dwMhe6dx46gBy0xN54fMdrNQVWkopFTM02YlhDofUD129tnJ3k+e4nY6wdxW3ox/IEaNdO5HscTWufxbnjgjufzttZE+cDqHa62fmpOMYPyALgHkbD0Q1TqWUUpFrVbIjIiNFpI/1/SQRuSBWigg2Fk9FBUP99N9rmuyRiWS7CDvSjlidv+L1B0j2hPdPeXSfjKP2KHty5nj6d0tmT3EVALdMGUh5zbFX0imllGpeNIsKtjhBWUSeA3xAqogECBYPLAJuAb4ZlSiiqKCggN69e9sdRlQFrJ6btMSjf1zh9kzkF1WSf7gyKnGFIxZTnS0Hyrn++WWM7JXe8sktGJabBgR7sG48dQDTRuby8rJd3DxlUJvbVkqprigvL49bb701KglPa1ZjDTLGnAUgIuuNMaOs7+e1+eqqRcYYxv76QwC6pXja3N4/l+zknbV7+b21DUJHicUp1NVeP/6A4ZITe0WtzedumEBKgovCshqufHoxkwdnM6p3RtTaV0opFb7WDGNVisg9IvIAcFBEbheRq4DYXUscR8qqfZRaheuiMcnX6RACHVwA75lvj8fnD/D6yt0xtU9W3W7st0Sx9yUz2YPb6SAz2Q3AK8sb1hcKtLWgj1JKqbC1Jtn5JrABeBW4gOAf6dnAVW25sIgkisgLIrJERFaIyHnNnHejiKwWkQUhX6e35dqdSUnVkZzS5Wj7fHI7lqqfNiSHV5bn85N/r4mpeSxbrf3BIpmk3JIkayn7wfKaBsdP/u1HDX6mSiml2l+Ln57GmGpjzH+NMV8aY6qAdGPMX4wxxW28dh4gxphJwLXAyyLSs5lz77K2paj7+qyN1+40QpODaMzx9bgcvPq9yW1vKAwup7BxXxlATPXsPP/Zdi45sX3md4kIY/pk8N66fXywfh8Aj83dzMHyWh5696t2uaZSSqmmRdJVcEFbLyoiDuBm4DkAY8xmYBUws61txxt/wJDicXL+qObywPDU+gJMtJZHdxSP08ENk49jUPeUmNkBffP+Msb1z+SJa8a12zVuOHUAxsDq/GJqfQH+8PHXABRXepmvS9OVUqrDRJLsRKPPfxDBobCNIcc2ABOaOf97IrLQGsL6fhSu32kYY+jXLRmIXhXi9hi2ael6e0qq2VZYETM9O7OX7OTttXvb9RqJ7uD/Xq8uz+dQxZHhrMLyGu57a127XlsppdQRYe2NZYnGhI+6borQ0rXFwMgmzt0PfAi8AHQHFoiI0xjzRBTiiHmnD+3OpEHd+MNHXxONPNOu6bHPfHs8r39RQJXXT1FFLWmJLnx+Q1KYNW6i5e+LdzKoe0q7XiPBFXxvhypq+d8PNtUfL6/2kZ7obtdrK6WUOsKunp06jT97j2rbGPOeMeZvJugA8ATQbO9OXVHBuq/OXlzwptMHMqp3BobY3V+qNUSEAdnJVNb6+fGrq/ns64Oc9MBcW2OaPjp6S86bUtezA7BwUyHfGBucH+Q3BrezE/8wlVKqHeXl5dV/hkerqGAkyU40hpHqJixkhhzLDDl+LLuA45p7sqCgwNp1OvjV2ZOdeJLodnLD88tIcjuprPVTZdOQVt2S89unDm7X69T17ECwd6dfVjI5qQm4HBKzFaWVUspueXl59Z/hBQUFUWkz7GTHGPOViORa20ZMEpHcCK67lWAV5mEhx0YCyxufKCI/a3SoJ7Angmt2ascqjRPO0JSdH7F1PR0Oh+AL2DdR+c/zgzuSh7tNRLhCe3YAyqq91Pj8uJ0OnJ25m04ppTqZsJIdERkgInOBncCb1tdOEflIRFpdmc0YEwCeBb5rtTsUGAvMFpERIvJxyN5bF4rImdZ5ycD3gH+GE3e86Owfj3U9HakeF/tLj72be3v648df88gVJ7T7RO0ElxNXSA9OWY2PGl8Al1NwiHDV04vb9fpKKaWCwu3ZeRx4DEg1xuQaY3KBVOvYY2G2lQeIiCwBXgKuMcbsAzKA4UDdDM5HgV+JyHxgITAfeCjMa8WB5vtvOksS5HEF/7mlJrr4vw83A/DZ1wc7NAZjDOeO6MEFoyPpkAxPtxQPEwd0Y/qY4LV+dv5wXrttsjWMBUu3F7V7DEoppcJfjbXXGPNu6AFjjBd4R0QuDqchY0w1cGMTx5cAfUIevw28HWaccSkaPRF2blZQ18vhdjrITvGwt6Sal5fv4vShOR0WQ3mNjwkDunXIaqjuaQnMvvkU8uasByA3I5HcjEQCBjw6Z0cppTpMuD07uSJysYjU70gpIh4RuQRo/z+Vu7Ae6YmkJjSdm3aWOTvZqQncec5QElwOAsbwp2vH8fbavfj8HTd/J7+oil4ZiR12PYdDGNc/s8Exf8Dg0Dk7SinVYcJNdu4EbgMqRKRIRIqAcoLzaO6KcmwqxIOXjWby4Owmn+uMH5s1vkB90jHkF+912HX/NP9rkj2RlJeK3IxxfZk0qFv9Y2NMfS9dUUVth8bSWfl1A1WlVBuElewYY/KNMRcTXBF1jvWVa4y52Bizsz0CDFddnZ14W3Le0VWP20vdu6j1BchM9hzz3PYwPDedqcO6d/h1X771yH5kfmtJJcDWwvIOjyUWHSir5q3VzS8xnfzQx80+Z461VFEp1Wnl5eXZWmcHY0yRMWaV9RVTsyzr6uzEW7ITDev3lNTvz2QnA1TW+slJSSDZ4+RH5wztkOtWe/34Awa3s+27x7eFPwAB6wN69a5iW2OJFVN+N58nF2xt9vkDZTUEmujdOVRew4vLdrVnaEopm+Tl5dlXZ6c5InJptNpS4XE7hVpfy/NeanwBvt/OhfRa4nII/kCAZ749nrREF+/eOQVDxwxTDL/vfebFwAacgYAhEIBLTuzNg7oDOgDnj8rlvFG5LNtexKwm9g1L8Tj59dsb+Oun2xocL6v28eLSXR0670sp1fmEW2fn+Oa+gCvaKUbVgpQEF5W1vhbPq6jxcdbwHh0QUfPcLgdev+G8Ubk4HMKAnBT+vSK/wd5R0VRR42PHwYr6x2XV3na5Tjj8xuByCm5dkQXApn1lTByQhQDvrdvL3xfvxNsoeUlLdPPC5zv46Kv99cd2Hqpg9+Eq1u8pZeO+spjZZFYpFXvC7dn5AviK4G7ljb+ui25oqrVSPC7Ka1pOdsqqfaQlduzk3MY8TsdRvVBVXn+rkrVIbCus4MMN+wC49YxBzP/p1Ha5TjgCxuBxOnTLCMvq/MOceXwPCstr+NuiHcDRc5nq/t2u2HG4fjjro68OMPO5pQA8/ck2dh+u7LiglVKdSrjJzizgdGOMo/EX8Jd2iE+1QnJCcK+plhRXeslM6vhJwaHcTqG20V/twrG3w2iLy/6yiN++u5HDFbUkuZ0xMdE7EDAkuHXLCIAH39nAnDV76dctiesnH4dD4PJxfVi+4zCPvL8RCE5A7pkeXLmX5HZSXBXsnTtUXlPfzpw1e9haWKGTlZVSTYqkgvKoZp6b3bZQVKRSElrXs3O4spbM5PYvpndMTXzAX35S33aLq24u0JUxtDWD3xgSXU4c9s6Ttt3Ly3ZRWuXj6W+PR0QYnpvOmlnn8eCMMdz35jr+smArFTU+npi3hVOHZNMtxcMtZwziQFk1B0qrmb+psEF73/vnSv69crdN70YpFcvCXXruN8b8tZnnlkQnJBWuFI+LypqWe3ZqfAES3e27+WVL/P5Ag/2iAH550Yh2K7I3Y1wf3v7h6Xx9oJwdhypafkEHCARgTN8MxvXLsjsU2xhjeHVFPuW1PlJCimWmJbpJsjZoTfY4ufnvK/j93M30TEvki/umMXlwNoVlNTz83ka+2lvK9DG5fHj3GQztkcoZx3fnN29vsOstKaViWKuTHRGZaG3Yeaxz0kXkoraHpcKR4HJQ42s52fljDCw79xuOSmxEpN2WxPfvlsyo3ul857QBHbbEvSUBYzh/VC6XnNgbgL8t2m5zRB2vbrl4dkrTw6qzLhnJHWcPYfG2Q6QluBjVJx2AHmkJfPu5ZbyxqoDtD03nz9eexPE908hK9jCqdzql1T5dmaWUOko4PTurgf8Tkf8RkZNEJEdEkkQky9qp/BbgDYKTmG0Tr0UFjyXBffSk36b065bUAdEcWyBgcDYzMffnb6xl/Z6SNl9j7e5iAF63hjREhFmXjGJQ99Q2tx0NddtFJFgbo76yPN/miDre9sIKRvRK5+5zj2/y+e+cNpALR/dibL9MxvbPZHhuMNnpk3nk37CI1M/Bmjgwi9G9M7jjrCFUt+L/BaVU7LOlqKC14ecVgJPg/JwDBLeKOAjMJTiX5ypjzN6oRBahrlhU0ON0UNPCL/haX4Arx/froIia5zfNJzsHSmsoOFwVcdurdh3mntfX8tt3v2JNfjE/+feao4bMYkHAGByO4L5ZF43pddQy664g2ePkwRljyGqmZwdgYE4Kj101ll1FR1ZZuZwObpky8KipX//v/OFcdEIveqYnUFypW3AoFQ9sKypojPEaY35rjBkBJAB9gTRjTF9jzF3GmENRiUqFxeNquWenosZHqs3LziHYq9FUsnP3ucfz8cYD3PrPlUc9V+31c6i8hmueWcKLS3c1u+LmzVUFvLw8ny92FvONPy8CYEiP2OjNCeUPmPqVWL5AwPaKzh1tb0kVCa2cO9Y7M7G+2nSde6eP4NYpg5o8v7zGz50vrdK9tJRSDUT8W9ZKfPYaY7S4hc0SXM5jztkxxlBW7Wt21/SOVO31k9TEB12iu/l/igs2FXL+45+weNsh7v3PlxSGLDkO1T0tAYBaf4Azjg/ufzWiV3oUoo6u0ITP6zd4/QEWbTloc1Qdo9rrZ/JD8+ib1boh1QSXkwU/PavBMRHh59NHNHn+7VMH88WuYu56ZXVbQ1VKxZGu9SdlnPK4jj2MddEfP2PHoYoG8x3skuxx0TM94ajjrpDejbpeqrq/zgPGcLA8ODQxpEcqxZVHV0EOBAxlIcvvJw/KZvNvLmRATkpU448GY45s7Or1B6is9bOuoO1zlWJZtdfPxn2l7C+t5rQh2Vw0plerX9vcsOexzFmzR2vuKKXqabITBxJcDmr9gSY3SgQoqfJy/fPL6NctuYMjO9ptZw7iyglHzx2qi33yoGwOWj03pz4c3On6+7OPzHlP9jg5XHH0nIwtheU8vXAb8386laxkN9PH5OJxxeY/78bzlmp9Aaq98Tlv540vgpPEC8tqmL+xkIPltdx6xuAGyW20fXHfNCYOyOLsRxdyxZOft9t1lFKdR2x+Gqiw1G3BcOXTi9lbcvQE37Ot/bByMxI7OrSjhK6gCeWzkp2Zk46jrDrYQ7O/tIZl24vqz5k6rDuTB2fXV9ANVeMNcMGoXLJTPXj9huOyY69Hp86d5wxtMJRX6wtQ3YrSAZ3N8h1FfPp1cHiutNpLabWXV5fnN7vcPFq6pXiYPDiH7QcrWLHzMAPueYfZS3e26zVV/Cmp9LJ+TwkfrN/HvpJq3lm7l0PlNdpj2ElpshMHHA4hYMDlFLYeOLpwXjfrwyWWJ8JmWRWUkz1OKmt9VFnbX1z59GLuPGcoPdISeOE7J3PD5AEUV9Ye9QunuKqWm6cMJNXjYs4PT+/w+MNx25mDG/Ts1PgC1MRJz07djuWr84t5/KPNFBQHk+83vijgcEUtr6zI75AexnH9M0kPmZD/i/+sa1WVcaXqfPTVfi7642d8tGE/s/67jh+8+AXjf/MRG/eV2R2aikCbP/1EJFdEzhSRW0TkURF5MwpxqQiM6JXOV3tLm3zukStO6OBownPVxH5sf2g6SR4nVV4/b60+stxQgOzU4DyfzGQ3hyu9nPm/C+qf31ZYzoY9pWQmu3E4hIExOE+nOcYEJ1THS8/Oy8vzOf1387jsz4tYtOUQpVVejDH8Y/EOlu0o4snrTiIjqf23LDlrWA8+vPtMPE4H9186ip7pCfz1022UNNErqLqWJxdsbfGcWl+AL3Yd5tffCO6OdPLAbBZYmwhvOVCOMUaLV3YyESc7IvIPEVkILAT+AJQC/wCujlJsEemKRQUhmBCkJbopq276l3lT82RiSd3wVrLHyaItB7nnjS8B+N6Zg6jxBchJDfZOJbmdfPzVfnYVVfL7DzfhDxjOfnQhD723kb5Z9s9JilS1Nz6SncxkN7utWkkep4OCw1X8/fMdeP2G4kov/bM77mfULcVDwBhuOHUAT397Ao9/9DU7Y2TLEGWf372/sdnaVgfLa9h+sIKnFm4l0e2kb1YSB8pquOn0gQzISeH9u6bwzCfbeGlZPs9+2vUqn3c0W4oKNmGjMeZMYATwFHAtkGyMqY5KZBHqikUFAeoHdUTIL6rkF//5kkfe39jpxpeT3E72Fgf/Ca3NO4+B2SnsLakix+rZEZH6InN/nLeFlTsPA/DUzPG27/sViQcuG83QHqnsKa6ipMrLgHveYdehzlXNYfbSnQy45x1+/+Em9pfWcM3JwcR6+S/OZcrxOby1Zg8AQ7qndugkeY/LwTPXjwcg0+pNOtTE5HYV/+Zu2M/a3cUs215EZrKb0bM+OOqcfy3ZybXPLuG3737Ff9fs4b6LR5KV7CH0N+jw3HT2FFfxu/c38uyn28gv6lz/r3Y2thUVFJH+IpILYIz5rfXfgDGmLtm5UET+FpXIVEQEuOzPi5i9dBd/WbCVwrKa+o0VO4Mkj5PiKi+/v/JE0hPdjO6Twer84vqeHYB+Vg+O0yHc8/pa/veKEzhnRA+7Qm6TgTkpOB3Ckm1FXPLEZwBs2t955gTkF1Wy+3AV4/pn8sd5W/jh2UO458IRLL33HDKS3QzPTWfVrmKemjmekb3TSU9s/yGsUGcP7wkEe5ymDM1hzuo9rCsoYfmOIq5+ZjFAl6xgHc9Kq71HDTPd8o8VbCusYNZ/19MjLQFfwNT/oQTw10+3UVBcxbSRPamq9XOntY/eiX0zedZKmOss/8W5vPejKbz1g9OYv+kAv/9wU8e8MdUm4VaZex14TkTetR5XGWMKAYwxFcCvRKR3NANUrfPW6gKmj+mF2+lo8NfrNc8u4dpTjrMxsvBkpySw41AFJw/sBsDQnqn837dOZHDIvlbPXD+B7QcrcDmE17/YzdnDe8T05OuW7C0J9mTV9Vh1hiEtrz/AqyvyeeT9TaQmuHji2nFc/pfP+cFZQ0h0O+vn5QzISeGEvhmcN7Jn/apAO6QnuvnlRSN558u9XPzEZ3xrfF/yi4LDbb+fu5n/uWC4bbGp6Drn0YX0zkwiEDC8dvvk+uO/n7uZIT1SuefC4TwxbwvffPJzvjW+L/9euRuXQzh3RE+e+vZ4NuwpZVD34Lw/h0NIcDT8Y9HhEHpnJmGM4VdvrQfg7mnHN7nKVMWOcD8hVli9ON8B3gbOb3yCMWZPNAJT4fnmSX3ZW9xw2fmdZw9ha2EFlZ1oFUqSx8m+kmrSrQ/LBJeTiQO61a8og+BcjPHHZXFiv0x+/Y3R9ZOXO6u8S0c2eFzVCZKd4kovv3tvIyVVXhwOOKl/Ftsfmn7UUOKlJ/bm7mnH43CIrXWPHA5hWG4al54YLGa4eX9ZfXXt/SVHRt7nrNnDnDX6K6yzKqv2Mrp3Omvyi/myoIRb/7GSn7y6hv+94gR2FVVSWFbD8T3T+MNVY/n2pON4f/0+rhjfl7fuOI0eVrHTkb3TWzUkLiL1ldrnrLV1S0jVCuH27BgAY8z9IjLcGPOvuidEJMEY03Qdf9Xu+mQm1a80ufHUAfz0/GEUldeytfBIL0lnUVnr7/DhDjvNGNeXu19ZA8APzhpcX2colpVVeymt9nHXuUO5y9q5vLm/bM8aFjtDjEN6pLH9oekUFFfx7xVWwcOQ7UeKKmpZnV/MhaNzcTqargmlYsfhilqyUjzMemsdCW4nl57Ym6tP7s/zN05k0/4yLnj8UwblpHDXuUNZ8NOpbLcmqDscwvfOHMR3ThtAv27JuJ0O7p0e/j56j181ls+3HuSOF1dx8ZheOGJw42EV1JY/tRrPfP3ftgSi2mbGuD7cO30EAqQnukhNcNE/O5k/X3cSpwzKtju8sLxy6yS7Q+hw861lrVnJHkpjfHn0joMVzN9UyFMzx3NdJxoirSMi9M1K5rMtBxlwzzsUltXUbxzrcAj/WVXAB+v389TCbfxn1W6bo1XHctrv5vHel3vZV1rNnuIqvvPCcsb2y0REGNI9lYcuH8O8n05lSI80BuSkNEi8+2YlM6h7av0QeCQLHLqleOiZHizW+tQnLS9pV/YJt2fnchEZa30/VERCa7EPAO6MRlAqfHV/UVTW+ugfw9WDW6OzJWfRMCA7mS/um0ZGkptTH/6Yu6cdb3dIzbriqc85WF7L9oemd+qej2G5aWQmuVm87RCVVhHLundT4/NzsLyGuRv2MWNcX/uCVMc0pEcq98/ZwIQBWTxyxQkcKq+tTz5cTgfXnNy/3WOoG2J/5P1NnD8qt8H8QhU7wk12NgB/b+a5mW2MJSrq1uTPmjWryy0/h+DS7ZExuNO3OjYRqf+lGcvVAkqqvBwsr2Xrbzt3ogPw2xljeHvtHj7eeAAR+PWcDQy0Jqa+tnI3qQkuvthVjDGm07/XzmxbYTnZKQkkehz4A4Z/LN7Jwk2FpCS4mD6mF1dO6Ee110+yx0Vyt3A/0tpuUE4KX9w3jUc/3MS03y9k20MXdXgM8SovL4/7778/Km2F+y/jl8aYJnfWE5GYWH9XUFBA795dd0HY+AHdGNErze4wVBtc3QF/jUZq3sb9zLnj9Ih2Io9FdavGHCL8a+lOfnb+MAB2HqokNSH46/FwpZd5Gw9wxfiu2cNTVFFLWqKLKq+f/SXVDO3ZPr9fqr1+PE5HfS+11x9gf2k1z3wSXBa++3Bwafg/Fu+g2hvgkStOYFy/zAaLF+xQ94eK1x/o9IslYk1eXh633nprVAoLhjVnp7lEx3puSZujUW125vHdSfZ0/F83KnpcDomp2i8b95XyjT8v4tXl+fx39R5G9Y6fnsPJg7L56XnHM6p3OkluJ4+8H/yb7VBFDZXe4ETxA2XV/OqtdZQ2U538WNbvKcEfiI2uuqpaf32R0dbG9ONXV3P3K6t5beVu/vjR10x77JP655raLiEQMLy6Ip/Pvj7Y7NY1tb4Av56zgX0l1RwOKZPx8Hsbee2L3ew4WMGdL62yth6Zz9rdJfWJ0L+W7OR7Zwxmza/O48oJ/dot8YrE5v3lDM9No6LGx/5SW2vrqibop6JSMaZbioc9xVWt2rk9YH1otecqkCcXbGVNfjFJbgf/uumUuFpx4nI6uOPsoWQke3jg7Q3UWh/g1d4AaQluHr9qGAdKa+iblcSbqwroluIhOyWBnFRPqz5o/71iNzMn9SclwUVhWQ0n9M0MK75D5TVkJXuics9f+HwHJw/Momd6Ine8uIqbTh/IKYO60SMtkdJqL/M3HiA7JYEl2w6x+3Alj101lrW7S+iemsCbqwrISvZw5YS+VNT4KK/x8fB7G8m7dBRpCS7eXbeX3PREfvTyagqKq/jG2N6UVnk5e3gPLj+pL1c/s4TfffMEZv13HX2zkvnPqgLSk1ysKyhh+phe7Cmu4s3VBbzw+Q6SPU5OHZxNeqKL52+cwMZ9ZWQle5gxrg8FxVUxOyfm2esn8IePN/PC5ztYnV/Ms9dPsDskFcLWZEdEEgluNTHciuVeY8yHzZz7E4JVmgFeNsbo6i8Vl7JTPJz5vwtYm3dek0vw3/tyL0keJ+v3lLJ46yGOy07mwRljonb94spaDld6yUp28/yiHZRUednxcHzPQ8hKdlPra9hTUV7j46T+WSzbUcSpg3N4eVk+1T4/kwdlMyA7pVXJzr6Sal5als+2wnK+2FXM3B+fQY+0xAbnrNxZREaSmyE9jm7v+ueX8YerxzGkRyrvrN3LRSf0ivg9FlfVsnR7EcluJ9kpHl5atosNe0sZ0yeD++es58S+mSzZdoge6Yk4Rbj22aVcOaEvH67fz2u3nwrAO2v3kn+4kj/N28KuokrmrNlDcWUtL3y+k6paH5eN68OWA+Uc3zONc0f0ZPvBcr77wnLuu3gk/1yykxP7ZvJlQQlTh3VnX0k1Jx2XxSPvb+LOc4bym8tGM7JXOoMaJTN1VbCBmE10ALqnJTD+uCzufmUN15zcH3/AxM1wbzywu2cnDxBjzCQROR5YIiIjjDH7Q08SkQuAW4Cx1qHVIrLBGPNOh0arVAeom4Pw2NzNzLpkFJW1PvKLqujXLQmnQ3j20218sauYft2ScDkciMBTC7dy25mDw7rO/I0H2LC3lMmDs0nxuBjaI5WdRZVs3FvKYx9t5tuTjmPFjiL+cPW49nibMWXigG50T0tAgANlwbo7r90+mfRENwfKgkUuzxrenX+v2M3spbu48+whLbbpDxg+23KQcf0zSUt0Map3OvO+OkBmsocLRufWn7dgUyFz1uzhtdtPJSvZw7efW8qLt0xiW2E5WckeCoqryEn18OqK/BaTnapaP95AoEGS/P66vZRV+/D7DX6/4fNdh+iRnsDUYT14bO5mvt5fxvWTB7D1QDmXjevDlRP60S8rmS8LShjXP5Nbpgw6cp8GZvGX+VsZkJ3Cn649iev+uoTUBBcf/+RM0hNdR03kHpabxtRhPUh0O5ut9/X9qS3fy85ixri+1HgD3PPGl+wqquC5GyZ2yj374pFtyY6IOICbgcsBjDGbRWQVwVVdjzY6/XvAi3WbjIrIbOA2QJMdFXf6Zyfz/amDEYG8/65nf2k1763bR5/MJCYNyuYPV48jNyMRITgMA8Hz/rZoO985bWCrrvHmqgK+2ldKstvFc59ux+UUDpXX8mVBCeeO6MkjV5zI3xZtZ/bNp3SJlUg90xM5e1gP7v/GKEqqvMEeDqsHprLGj8Mh3HXu8byyPL/J13+5u4QxfTMIBAzPL9rOzVMG8cDbG/jVxSN55INNTBvZkwSXk3ve+BKAz+85m96ZSczbuJ+VOw+z+3AVE37zEX+4eiyb95fj9Qd46L2NDMhJ5umFW/nZBcM5UFaD1x845tYoCzYdYPP+crJS3PTLSubJhVtZtr2Iy8b25orx/Th9aA5VtX72l1bTr1sy54/KpaLGZ00GriEjyU1uRvB9nz4056j2e6Qlct0p/cmyEvJdRZX89wen10/0bkpX+7C/+uT+jOiVzjf+vIi/frqNrw+UH/UHw6ItBwE4bcjR9xiCPYLd0xJwSPPFOuvoasHWEbt2xRaRIcDXQE9jzAHr2BNAjjHmmkbn7gZ+bIx51Xr8TeBPxphe1uPeQAHoaiwVXypqfCR7nBwoq+FAaQ3r9pQ0WTukrNrLK8vzKa3yMrxXOpW1fvaVVDF5cDYZSW5KqrwUFFdTVF7DvtIaxh+XxbSRPRu0UVnrY+ehSl5fuZtfXjzyqGt0Vfe8vpaRvdO5fvIAdhysoGd6ItP/+Cm9MhLJzUiksKwGf8CQleKhutbPxn1lTB3WnS8LSnjsqrE8/N5Gkj1Orhjfl037yiit8vLWmj1U1vopLKshJ9XDN8b2YWthOQET3L7i5ikD+d37m/jWhL5Ue/2kJrjISHKz+3AVn35dyJSh3Zk4oBvnjuzBn+Zt4ZPNhZw+NIcVOw4zMCeFguIqRvZOJ9Xj4pITe9OvW3K7DKnU+PwkuLpWMtNaFTU+/vrpdgrLq8lM8pCT6kFEOFBWTa0vQFm1j9yMRG47czCJbic+f4DyGh+ZyR6ufHoxw3PT6N8tme+cNhCnQ1hXUMJXe0s5ZWA2qYku3vlyLwOyk/nboh2cNiSHS07ohS9g6J6WgDG0y/YsgYBBBNbvKaVPZhIHy2uo9QcYnpvOkm2HGNM3g7QEF6vyixnXL5OD5bV0T4t8hZoxhj179tC3b/1KyD6RbkllZ7JzGvAZkFi3zYSIPACcaow5p9G5XuBiY8wH1uNzgA+NMU7rsSY7ShHc8Xm3tcFlr4xEVu8uprzaR1qii24pHrKSPfTrlmxzlJ1LYVkNqQkukjwNP9SNMdT6A/gDhmSPC3/AcLiylvREN0UVtfU9JHXnhv717Q8YKmp9FJbV0Cczqb73o6LGR1m1j9X5xQztmUqvjER2FVVSVF7LqN4ZfLblIJnJbnYeqqSs2ktFrZ+RvdKZNKgbhWU1DOqeitMhBAKGgDH1PX/KPsYEd1h3OIRthRX0yUxize5ivL4ASR4nWwsr6J7qoazGRyBgyEgOJkZrd5dw6uBsNu0rwwAiMDA7JZhg+AIcn5vG5v3lDOuZhsspLN56iJQEJwETXCnndDhwOiBgwO10cKxUt8YXaJAc1Z1rrO+tfaLwG4PPb0jyOCmu9OILBOiVkVT/R1mV10+tL0BaopviSi8VNT56pAeTr4AJJkoux5HrBKz8wxhwiPUmQ3KS4iovgYoiHrzuzLpDESc7GGNs+QJOs+6hJ+TYA8C8Js71AueFPD4H8Ic87m21ddTX5ZdfbioqKowxxlRUVJjZs2eb/Px8U+eNN94wa9eurX+8YMECs2DBgvrHa9euNW+88Ub94/z8fDN79mxtU9vUNrVNbVPb1Dbboc1Zs2Y1+XkO9DYR5hx29uwMBTZz9DBWd2PM1Y3OLQDuNg2Hsf5sjMm1HmvPjlJKKRVn9uzZE1pUMOKeHTv7OLcCRcCwkGMjgeVNnLu8lecppZRSSjVgW7JjjAkAzwLfhfqenrHAbBEZISIfi0jdIPlTwDUikmjV5rnWOqaUUkopdUyxUGfnKRFZYsVyjTFmn4gMIFho0E1wbs77IjIKWGS97nmjNXaUUkop1Qq2JjsmWDfnxiaOLwH6NDr2KEfX31FKKaWUOiZdl6iUUkqpuKbJjlJKKaXiWtwlO3369EFEyMvLszsUpZRSSkUoLy8vdNl5m9hWZyeatM6OUkopFX/ioc6OUkoppVS702RHKaWUUnFNk50w6Dwge+n9t5fef3vp/beX3v/OTefshHcd4uF+dVZ6/+2l999eev/tpfffHjpnp5ViNRuPZlzRfo+xHFu0xPJ7jOXYoiWW32MsxxYtsfweYzm2aInl9xjLsbVFvPTs9AXyAVauXElubm79c3369KGgoCAq1+kKbUW7PW3L3va0LXvb07bsbU/bsre9aLS1b98+xo8fX/ewnzFmdyTtxEuyMw74wu44lFJKKdVuTjLGrIrkhXE/jKWUUkqpri1eenbcwGjrYSEQsDEcpZRSSkWHA+hufb/OGOONpJG4SHaUUkoppZqjw1hKKaWUimua7IQQkUQReUFElojIChE5r5nz0kTknyKyXESWWa9J7eh445GITBSRLSJyYwvnXSciK62f06MiIh0UYlxrzf0XkdNF5G0R+dj6GTwoIvq7JApa++8/5PwFIvJC+0bVdYTx+2eEiLwjIp+KyAYR+W0HhRjXWvn7xyUifxCRL0TkcxF5U0Rymzu/jv6CaiiP4NDeJOBa4GUR6dnEefcBxwGTrK/jgF92VJDxSkRmAHcDJS2cNxp4FDgfOBk4Cfh+uwcY51p7/4HfAM8YY84BzgS+CfywncOLe2Hc/7rzLwLGtWtQXUgYv38ygNnAD40xUwj+P9C//SOMb2H8+78VuAg4zRhzKlAKPN5S+5rsWKy/TG8GngMwxmwGVgEzmzh9FLDMGOM3xgSAZegvnWhYboy5Fihr4bybgHeNMQet+/88cFu7Rxf/Wnv/3wTmABhjyoG3gSZ7QVVYWnv/635f/ZTgv30VHa29/zcDc40x2wCMMYXGmKY+J1R4Wnv/RwGrjTFV1uPFtOLzV5OdIwYB2cDGkGMbgAlNnPsucLaIpIhIMnA2sLT9Q4xvYRSLmsjRP6dRIpIU/ai6jtbef2PM46bhyoZEgqsgVRuEWSzt28A7tLIXSLUsjPt/DlAlIm+IyCIReVZEstoztq4gjPv/HnCKiOSIiAu4kFZ8/mqyc0TdcFXoL49ioEfjE40xfwYWANuAHQQLGv66XaNToXpy9M9JgBxbounCRMRJcDjxSbtj6SpEJAG4BfiT3bF0UQOA2wn2rJ0OGOAlOwPqSowxbwOPEfyDd7t1+PaWXqfJztEar8U/auKriPwCOJHgXJ3+BHuFbmz3yFSopmom6CTljvdT4E1jjPZsdpw7gOeNMdV2B9JFJQBvG2O2WT2c/wecb21IrdqZiMwEvgMMJfgZXAzc09LrNNk54oD138yQY5khx0PdCTxtjKm2fuE8BTzQrtGpUAc4+udk0KGUDiUi0wkO87b4i0ZFhzU59jLg7zaH0pUdBvaHPK4bfulrQyxd0Q+BfxpjDltzNv8A3CMiKcd6kSY7R2wFioBhIcdGAsubONcDhFZx9AJp7ReaamQ5R/+c1odMWFPtTEQmAT8Cvm2M8YvIULtj6iJOBtKBj0VkAcEe5QusJeiT7QysC1lNw+kNddV993R8KF1SU5+/LiD5WC/SZMdiZYjPAt8FsH55jwVmWzUVPrbmJwB8BFwlFoLL1OfbEHaXYE1E+1REulmH/gpMF5Fsa1XKjQR711Q7aHz/RWQEwb+mbgZcVo2pX9kZYzwLvf/GmLnGmBONMVONMVOBF4D3rceL7Y00PjXx++c54GIRqUt4bgLmRbobtzq2Ju7/R8A3rW2iILhiep0x5pg9+5rsNJQHiIgsITjh7BpjzD4gAxgO1N3cHxCcH7KE4CxwF8G1/6oNRGS89dfqWILdkm9YTyURvP/JAMaYdQTninxI8P6vBv7SweHGndbef4LJ5snALoLLRMuAKR0abBwK4/7Xnf8GDXt2PB0XbfwJ4/fPYuDnwFwR+ZTgvozXdnjAcSaMf/95wHpgsYh8bp1/RYvt695YSimllIpn2rOjlFJKqbimyY5SSiml4pomO0oppZSKa5rsKKWUUiquabKjlFJKqbimyY5SSiml4pomO0rFKRF5W0T8IrJRRG4NOf6qiPwj5PGlIrLB+jqznWJ5XET2iUhee7QfRhwfi0i1iEwN83XJIvI7q4Bie8Q1SETutXZxVkpFmSY7SsUpY8zFBAsuzjHGPANgVZw+Gzgn5Lz/Am8A1xpjFrZTLHcB77dH280RkRutImWhcZwD7IuguReABcaY8iiEdhRjzDagBHioPdpXqqvTZEep+DYXODfk8ThgAdBbREaFHD8BWNOBcXUa1oan3Ywx77XzpZ4CZjT6uSilokCTHaXi21zgRBGp26zwXOBpYDswDep30i42xhgRucAa6plr7Udzl3XOQBFZLyK1IvK2dexsEdkkIl+JyBARSRCR/xORJSKyUET+JSJZzQUmIr1F5HXrOotEZJbV84SIzBaRYuvYKyKy0oopK+T1w63XLhWR/1rDTMUi8rKIXEpwN/ax1lYKHze6/FgReUNEvrSG9dw07zqg/vUicruI7BCRF6zHk0VktYjsCDmnbrjsLhH5jzWU+EcR6SciL4rIGuu69cNWxhg/wURUtx5QKso02VEqvn0GVHNk2Oo04FOCSdA069hZHNnINg242RgzDTgDuExEphhjtgPnEdwT7qcAxph5BPcmu9AYswX4HcE9s84wxpwJ7Ce4j9ZRrKRmDvC1MWYKwaG18wnuO4cx5jqCQ3BnAjONMeOtl95kvd5lvf59Y8wpwFVWG6uNMVdbQ3MPW4+nWsNXocYD3yTY03UC8I1j3MOTgS11D4wxTxIc1qp7vBi4K/QFIcNlY40xM6zrXAM8AFwPnASMbOK6m4GJx4hFKRUBTXaUimPGmBqCCc80EUmwjlUT3Dn4DKtH4xzrMQQTjN+KyCKCCdBwYJL1ugLgPaxNb61elgxjzA4rebkVeNoYU2u19XeCwzJ1uxWHOpngB/4TIXG+THAn9VBzjDFe6/uVwFDr+0nAEKxkyhhTBbwSxq153QT5rPc89Bjn5gKlYbQd6s2Q+DYRTL58Vi/OGuD4RueXWtdTSkWRzvxXKv7NBe4ATgcWWcc+JriL8GSgnzEm3zr+DsGJxNdaw1ovACkhbT0L/E1Efg58G/indTyH4O7EPxWRW6xjToI7o/cCihrF1N/672wRqTuWwtF/gB0O+b4aSLC+72P9tzDk+YON3/gxFDfTblMcQKQ7Jodex9fosbeJ6wbQP0KVijpNdpSKf3OBR4DvA78FMMYUicgq4LtAPoA1r2co8ANjTN2Hu6dRW+8CVQSHgC4FLrSOH7SOz7KGkLDazKbhB3ydXdZ/LzXGlFrnCsGkqTUKrP/24Mjqqta+Nlz7CQ7vNZYY8n1TvVeRyLCup5SKIv0LQqn4t4ZgD8gZwKqQ43MJ9s7MtR4fIpi0nAH1E5enhDZkDb88T3B+zoq6ISZjTAB4Bvhu3aRbERlGcAitKcuAL4DvhRy7Hvh9K9/TEoLzaG6xrpVEMAELVYqVpFgTp09pZduNfQEMbOL4WBFxWUnajAjbbmwQDX9GSqko0GRHqThn9dJ8BMy3kpI6cwkOmyywzgsQnET7LRH5DPgz8DVwo4jcHvK6vxKcV/Jso0v9D8F5KZ+JyDzgj8CVxhi/iDwOXGC19f+sa10CnCwin1n1cM4jONyGiDwNjAXuEZEZEiyKeCNwgYjcb821uRQ4X0SWAS8CH9JwuOkjwGe9lxOA1SLylhX749Yqql+FxBX6HkO9RHASd2NZBBOhOcByILdu1VcT12n8fu4PuW7dexZgqnU9pVQUyZHeaqWU6jxEJMcYczDk8c+AccaYa6J8HSE4fPegMeYz61geMMAYc2MUrzOT4Eq2W1s8WSkVFu3ZUUp1Vn+3hsoQkRSCvVJvR/siVs/YNcBMab/tIgYS7Pm5sz3aV6qr054dpVSnZA1tfQ8oI7iS6w1jTLtvtyAi9wC3EZyg/G9jzA/b+5pKqbbRZEcppZRScU2HsZRSSikV1zTZUUoppVRc02RHKaWUUnFNkx2llFJKxTVNdpRSSikV1zTZUUoppVRc+/9Ivphad/+4sgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "species.plot_spectrum(boxes=[spec_box, ],\n", " filters=['MKO/NSFCam.J'],\n", " xlim=(0.75, 1.8),\n", " ylim=(-2e-9, 1.8e-8),\n", " offset=(-0.12, -0.05),\n", " figsize=(7., 3.),\n", " output=None)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Synthetic flux and magnitude" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we use the [SyntheticPhotometry](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.photometry.SyntheticPhotometry) class to calculate the flux and magnitude for the [MKO/NSFCam.J](http://svo2.cab.inta-csic.es/svo/theory/fps/index.php?id=MKO/NSFCam.J&&mode=browse&gname=MKO&gname2=NSFCam#filter) filter. We first create and instance of [SyntheticPhotometry](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.photometry.SyntheticPhotometry) with the filter name from the SVO website." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Adding filter: MKO/NSFCam.J... [DONE]\n" ] } ], "source": [ "synphot = species.SyntheticPhotometry('MKO/NSFCam.J')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The average $J$-band flux is calculated with the [spectrum_to_flux](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.photometry.SyntheticPhotometry.spectrum_to_flux) method. The error on the synthetic flux is estimated with Monte Carlo sampling of the input spectrum." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Flux (W m-2 um-1) = 1.80e-09 +/- 8.58e-14\n" ] } ], "source": [ "j_flux = synphot.spectrum_to_flux(wavelength, flux, error=error)\n", "print(f'Flux (W m-2 um-1) = {j_flux[0]:.2e} +/- {j_flux[1]:.2e}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly, we calculate the synthetic magnitude with the [spectrum_to_magnitude](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.photometry.SyntheticPhotometry.spectrum_to_magnitude) method. Also the absolute magnitude can be calculated by providing the distance and uncertainty (set to `None` in the example). In [species](https://species.readthedocs.io/en/latest/species.html), the magnitude is defined relative to Vega, which is assumed to have a magnitude of 0.03 in all filters. For the selected $J$-band filter, Jupiter has a magnitude of 0.59 so the planet is comparable in brightness to Vega." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Downloading Vega spectrum (270 kB)... [DONE]\n", "Adding Vega spectrum... [DONE]\n", "Apparent magnitude = 0.59 +/- 5.03e-05\n" ] } ], "source": [ "j_mag, _ = synphot.spectrum_to_magnitude(wavelength, flux, error=error, distance=None)\n", "print(f'Apparent magnitude = {j_mag[0]:.2f} +/- {j_mag[1]:.2e}')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.0" } }, "nbformat": 4, "nbformat_minor": 2 }