diff --git a/tutorials/units-and-integration/units-and-integration.ipynb b/tutorials/units-and-integration/units-and-integration.ipynb new file mode 100644 index 0000000..eb9e433 --- /dev/null +++ b/tutorials/units-and-integration/units-and-integration.ipynb @@ -0,0 +1,550 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "WIVPH0GmXNbD" + }, + "source": [ + "# Usando `scipy.integrate`\n", + "\n", + "## Autores\n", + "Zach Pace, Lia Corrales, Stephanie T. Douglas\n", + "\n", + "## Tradução\n", + "\n", + "Luciano Cordeiro\n", + "\n", + "## Objetivo do aprendizado\n", + "* realizar integração numérica no contexto `astropy` e científico python\n", + " * aproximação trapezoidal\n", + " * quadratura gaussiana\n", + "* usar as curvas de corpo negro embutidas do `astropy`\n", + "* entender como as unidades do `astropy` interagem umas com as outras\n", + "* definir uma classe Python\n", + " * como funciona o método `__call__`\n", + "* adicionar rótulos $\\LaTeX$ em figuras `matplotlib` usuando o formatador `latex_inline` \n", + "\n", + "## Palavras-chaves \n", + "modelagem, unidades, synphot, OOP, LaTeX, astroestatística, matplotlib, física\n", + "\n", + "## Conteúdo complementar\n", + "* http://synphot.readthedocs.io/en/latest/\n", + "* [Using Astropy Quantities for astrophysical calculations](http://www.astropy.org/astropy-tutorials/rst-tutorials/quantities.html)\n", + "\n", + "## Sumário\n", + "Neste tutorial, usaremos os exemplos da função de Planck e da função de massa inicial estelar (IMF) para ilustrar como integrar numericamente, usando a aproximação trapezoidal e a quadratura gaussiana. Também exploraremos a criação de uma classe personalizada, cuja instância pode ser chamada da mesma maneira que uma função. Além disso, encontraremos as unidades integradas do astropy e teremos uma primeira amostra de como converter entre elas. Finalmente, usaremos o $\\LaTeX$ para facilitar a leitura dos rótulos dos eixos das figuras." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "id": "RAdawwnRXNbH" + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "from scipy import integrate\n", + "from astropy.modeling.models import BlackBody\n", + "from astropy import units as u, constants as c\n", + "import matplotlib.pyplot as plt\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "OhPAWSAlXNbJ" + }, + "source": [ + "## A função de Planck\n", + "\n", + "A função de Planck descreve como um corpo negro irradia energia. Vamos explorar como encontrar a luminosidade bolométrica usando a função de Planck tanto no espaço de frequência quanto no espaço de comprimento de onda. " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "aqXb3VB3XNbJ" + }, + "source": [ + "Digamos que temos um corpo negro a 5000 Kelvin. Podemos descobrir a intensidade total (bolométrica) deste objeto, integrando a função de Planck. \n", + "A maneira mais simples de fazer isso é aproximar a integral usando a regra do trapézio.\n", + "\n", + "Vamos fazer isso primeiro usando a definição de frequência da função de Planck. Vamos definir uma grade de frequência de fótons e avaliar a função de Planck nessas frequências. \n", + " \n", + "Esses serão usados ​​para integrar numericamente usando a regra trapezoidal. Ao multiplicar uma matriz `numpy` por uma unidade `astropy`, obtemos uma `quantidade`, que é efetivamente uma combinação de um ou mais números e uma unidade." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "j07CMpw6XNbK" + }, + "source": [ + "
\n", + "**Nota sobre unidades de impressão**: \n", + "\n", + "Quantidades e unidades podem ser impressas em strings usando [Format String Syntax](https://docs.python.org/3/library/string.html#format-string-syntax). Esta demonstração usa o formato `latex_inline` que está embutido no pacote 'astropy.units`. Para ver formas adicionais de formatar quantidades, consulte [Getting Started](http://docs.astropy.org/en/stable/units/#getting-started) nas páginas de documentação do astropy.units.\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "id": "yJD8dBdYXNbL", + "outputId": "3f8f436e-d17e-4315-bcdf-cb743d3e6718", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 299 + } + }, + "outputs": [ + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEaCAYAAAAL7cBuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxc1Xn/8c+j3ZZlS7bkBduybDCLAwbbAkxYwg+yAFlI05CFJCwhpbRJk/TXpFmbkLRpm/7aNCE7IQTIvpKQhCwkIRDCamNjzGLjXbItS7a1W4slPb8/7pUZlJE0I83MnRl936/XvHTvnTv3PndGmkfnnHvOMXdHRERkpIKoAxARkeykBCEiInEpQYiISFxKECIiEpcShIiIxKUEISIicSlBSFLM7I9m9o5sOYeZnWRmG82s08zenc64Rpy31sy6zKwwDcfuMrNlE3ztX5lZQ3iMVamOTaYWJQj5C2a2y8x6wi+ZA2Z2m5nNiDquUfwzcK+7V7j7Tek6SfievHR43d33uPsMdx9M9bnC4+6Y4Mv/G3hXeIwNqYxLph4lCBnNq919BrAaqAc+GnE8o1kCPBV1EFlk1PfDzIoyHIvkOCUIGZO77wV+BZw68jkzO97M/mBmh8zsoJl928wqY57fZWbvM7NNZtZuZt83s7KY5y8Pq4c6zGy7mV0S5xwLwte/P85zfwD+D/CFsLRz4sjqKTO7xsweiFl3M7vBzJ4zszYz+6KZWczzf2Nmz4RVVk+b2Woz+yZQC/w8PM8/m1ldeKyi8HXHmdldZnbYzLaZ2d/EHPNGM/uBmd0RHvcpM6sf7T0Pj3tCuHxbGOMvw9c+YmbHx3lNqZl1AYXAE2a2PeYz+ICZbQK6zazIzNaa2YPh9T9hZhfGHGepmd0XnuseM/uCmX0rfO5CM2sccd5jJSszKzCzD4af5aHwmmeHzw2/X1eb2Z7w9+UjMccpNLMPh6/tNLP1ZrY4vPb/GXHOu8zsH0d7/ySF3F0PPV7wAHYBLw2XFxP8R/qv4fofgXeEyycALwNKgRrgfuCzI47zKHAcMBt4BrghfO4soD18fQGwEDg59hzAUmArcP0YsR6LZ5T1a4AHYtYd+AVQSfCl3wJcEj53BbAXOBOw8PqWjHxPwvW68FhF4fr9wJeAMuCM8LgXhc/dCPQClxF8gf8H8PAY1+TACeHybcCh8P0qAr4NfC+R18bEvTH8HKeF7/OhMJaC8P0/BNSE+z8EfCb8TC8AOoFvhc9dCDSO8bvyHuBhYFH4+q8C3x3xfn0tjON0oA84JXz+/cCTwEnhe386MCe87n1AQbhfNXAEmBf138lUeORVCcLMbjWzZjPbnKLjDYb/4W40s7tSccwc8lMzawMeAO4D/n3kDu6+zd3vcfc+d28h+GJ5yYjdbnL3fe5+GPg5wZcnwHXAreHrh9x9r7s/G/O6FcC9wMfd/eYUX9t/unubu+8JzzEc0zuA/3L3xzywzd13j3cwM1sMnAt8wN173X0jcAtwVcxuD7j73R60WXyT4AswUXe6+6PuPkCQIM4Y7wUj3OTuDe7eA7wVuDuMZcjd7wHWAZeZWS1BcvyX8DO9n+AzS9QNwEfcvdHd+wgS4+tHVG19wt173P0J4Amefx/eAXzU3beE7/0T7n7I3R8l+Efi4nC/NwF/dPcDSb4HMgH5Vid5G/AF4I4UHa/H3ZP9Y8wXr3X33421g5nNAz4HnA9UEPxH2jpit6aY5SMEpQkI/qO9e4zDvwXYBvwoiZgTNTKm4Qb4xcD2CRzvOOCwu3fGbNtN0HYz2jnLzKwo/NKfaLyJaohZXgJcYWavjtlWTJAojwNa3b075rndBO9LIpYAd5rZUMy2QWBezPpE3vvbCRLbPeHPzyUYj0xSXpUgwv94DsduC+vJfx3Waf7JzE6OKLx89O8E1QanuftMgj9eG/slxzQAf1GXHuNG4CDwHUvuVtJuYHrM+vwkXjtWTGMNe7wPmG1mFTHbagmqq7JBbOwNwDfdvTLmUe7u/wnsB6rMrDxm/9qY5Re8t+HnUjPi2JeOOHaZB+1Y4xnrvf8WcLmZnQ6cAvw0geNJCuRVghjFzcA/uPsa4H0E9cSJKjOzdWb2sJm9Nj3h5bQKoAtoN7OFBPXIifo6cK2ZXRw2bi4ckbyPErQJlAN3mFmiv6sbgdeZ2fSwofe6JGK6BXifma2xwAlmtiR87gAQt2+CuzcADwL/YWZlZrYyPO+3kjh3pnwLeLWZvSJsGC4LG58XhdVp64BPmFmJmZ0HxJY0thL8TbzSzIoJ7mwrjXn+K8Cnht8zM6sxs8sTjOsW4F/NbHn43q80szkA7t4IPEZQNffjsKpMMiCvE4QF9+6/GPihmW0kaDRbED73OjPbHOfxm5hDLHH3euBK4LPx7h6Z4j5BcBtsO/BL4CeJvjCsW74W+N/w9fcRVFHE7tMPvI6giuLWBJPE/wL9BF/otxPU2Sca0w+BTwHfIWic/SlB4zoEDcsfDe/8eV+cl7+ZoCF2H3AnQdvJmFV0UQiT2eXAhwka0hsIEvvwe3slcDZBSfzjxFTXuns78PcEX+Z7CUoUsXc1fQ64C/itmXUSNFifnWBonwF+APwW6CD4B2JazPO3A6cRJAnJEHPPrwmDzKwO+IW7n2pmM4Et7r4gBce9LTxuOurERbKSmd1IcFfUWyOO4wKC0s8Sz7cvrSyW1yUId+8AdprZFQBh0TWhu0fMrMrMSsPlaoK7VJ5OW7AiEldYnfUe4BYlh8zKqwRhZt8luI/7JDNrNLPrCO6Guc7MniC4nz/ROtFTgHXh6+4luDVSCUIkg8zsFKCNoGr4sxGHM+XkXRWTiIikRl6VIEREJHWUIEREJK686UldXV3tdXV1UYchIpJT1q9ff9Dda+I9lzcJoq6ujnXr1kUdhohITjGzUccbUxWTiIjEpQQhIiJxKUGIiEhcShAiIhKXEoSIiMSlBCEiInHlzW2uueTo4BCP7TpMcWEBq2urKCxIdI4dEZHMUYLIsO0tXfzN7evYcTCY1XHFgpl87ep6FlZOG+eVIiKZpSqmDGrvOcrbb3uM9p6jfPHK1XzmDafTcPgIb/zqQ7Qd6Y86PBGRF1CCyKD/95tnaWzt4atvW8MrVy7gdasXcft1Z3Ggo5cP3/lk1OGJiLyAEkSG7Gjp4ruPNnDlWbXU180+tn11bRXvfemJ3P1kEw88dzDCCEVEXkgJIkNue3AXhWa8++Llf/HcdectZfHsafznr59B83OISLZQgsiArr4BfvL4Xl61cgE1FaV/8XxZcSHvvPAENu/t4MHthyKIUETkLylBZMA9TzfR1TfAlWfXjrrPa1ctpHpGKTffvyODkYmIjE4JIgN+vbmJ+TPLWF1bNeo+ZcWFXHXOEu7b2sKu8BZYEZEoKUGkWU//IPdtbeHlL5pHwTgd4q6oX0SBwY8fb8xQdCIio1OCSLP7trbQe3SIS140f9x9F8yaxnnLa/jx+kaGhtRYLSLRUoJIs/u2tlBRWsRZS2ePvzNwxZpF7Gvv5eGdaqwWkWgpQaTZIzsOcdbS2RQVJvZWX3zKXMqKC/j15qY0RyYiMjYliDQ60NHLjoPdrF02J+HXTC8p4iUn1vDrzU2qZhKRSClBpNFDYZ+Gc45PPEEAXHrqApo7+9jQ0JqOsEREEqIEkUYP7zjEzLIiTlkwM6nXXXTKXIoLjV89qWomEYmOEkQabdjTxqoJzPcws6yYc0+o5nfPHEhTZCIi41OCSJPuvgGea+7k9MWVE3r9hSfWsOvQEXYfUqc5EYmGEkSabN7bzpDD6YtmTej1LzlpLhDcJisiEgUliDTZ1NgOwMpFEytBLK0uZ8mc6dy3RQlCRKKhBJEmTzS2sbByWtzRWxP1khNreHD7IXqPDqYwMhGRxChBpMmTe9tZOcHqpWEXnlRDz9FB1u3S7a4iknlKEGnQ3TfA7kNHWJHk7a0jrV02h5LCAu7b2pyiyEREEqcEkQbPNXcBcNL8ikkdZ3pJEauXVPLQDo3LJCKZpwSRBluaOoDJJwgIShFP7eug/cjRSR9LRCQZShBpsKWpi7LiAhZXTZ/0sc5ZNgd3eHTX4RREJiKSOCWINNh6oJMT51WMO0FQIs6oraS0qODYuE4iIpmiBJEGWw50ctK8yVcvAZQWFbJmSRUPqx1CRDJMCSLFDnf309LZl5L2h2HnLJvDM00dtB3pT9kxRUTGk/EEYWaLzexeM3vazJ4ys/fE2cfM7CYz22Zmm8xsdabjnKgtTZ0AnJiiEgTA2uODdohHdqodQkQyJ4oSxADwT+6+AlgLvNPMVozY51Jgefi4HvhyZkOcuK0HggSRyhLE6YsqKStWO4SIZFbGE4S773f3x8PlTuAZYOGI3S4H7vDAw0ClmS3IcKgTsvNgN+UlhcydxBAbI5UUFbBmSRWPqgQhIhkUaRuEmdUBq4BHRjy1EGiIWW/kL5MIZna9ma0zs3UtLdkxqN3Og90srSnHbPJ3MMVas2Q2zzZ10NU3kNLjioiMJrIEYWYzgB8D73X3jokcw91vdvd6d6+vqalJbYATtPNgN3VzylN+3PolVQw5bNijcZlEJDMiSRBmVkyQHL7t7j+Js8teYHHM+qJwW1brHxiisfUIS6tTnyBW1VZSYGjgPhHJmCjuYjLg68Az7v6ZUXa7C7gqvJtpLdDu7vszFuQE7Tl8hCEnLQmioqyYk+bPZP1uJQgRyYyiCM55LvA24Ekz2xhu+zBQC+DuXwHuBi4DtgFHgGsjiDNpuw4G04PWpSFBAJxZV8WP1zcyMDhEUaG6sIhIeiWUIMxsdgK7Dbl723g7ufsDwJgtuO7uwDsTiS2b7AwTxLI0JYg1S6q446HdPNvUyakLJzfXhIjIeBItQewLH2N9sRcSlgKmqp2HuqmcXkzl9JK0HL++LsjT63e3KkGISNolmiCecfdVY+1gZhtSEE9O29nSnZb2h2ELK6exYFYZ63a3cvWL69J2HhERSLyR+pwU7ZPXdh3qZmkabnGNtWZJFes19LeIZMC4CcLMXgZ83szOCNevj7efu/emOLac0jcwyP72XmrnTH4OiLHUL6liX3sv+9p60noeEZFEShBvB94PvNXMLgLOSG9IuWlfW5AfUzFJ0FiG2yHW6XZXEUmzRBJEp7u3ufv7gJcDZ6Y5ppzUcPgIAIuqpqX1PCfPr2BacSGPK0GISJolkiB+Obzg7h8E7khfOLmrsTWo8lk0O70liKLCAlYumsWGhnHvKBYRmZRxE4S7/wzAzKrD9c+nO6hc1Nh6hKICY/7MsrSf64zaSp7Z10HfwGDazyUiU1cy3XFvTVsUeaCxtYfjKqdRmIJ5qMezanEV/YNDPLVvQmMciogkJJkEkf5vvhzW2Hok7e0Pw1bVVgKwYY+qmUQkfZJJEJ62KPJAQ2tPxhLEvJllHDerjI1qhxCRNFIJIgV6jw7S0tmX9ltcY62qrdLcECKSVskkiA+lLYoct7dt+A6mzJQgIKhmamztoaWzL2PnFJGpJeEE4e6b0xlILjt2i2sGSxBnLA7aIVTNJCLpktSkAmZWb2Z3mtnjZrbJzJ40s03pCi5XNLZmppNcrFMXzqKowFTNJCJpk+yEQd8mGHbjSWAo9eHkpobDPRQXGvMq0t8HYlhZcSGnLJipO5lEJG2STRAt7n5XWiLJYY2tR1hYOY2CDPSBiLWqtpIfr29kcMgz0v9CRKaWZOet/LiZ3WJmbzaz1w0/0hJZDmls7WFhBquXhq2qraS7f5Dnmjszfm4RyX/JliCuBU4Ginm+ismBn6QyqFzT1N7L+curM37eMxZXAUGHuZPnz8z4+UUkvyWbIM5095PSEkmOGhgcormzlwWzMtf+MKxuznQqpxezYU8rbz5rSs/2KiJpkGwV04NmtiItkeSolq4+hhzmz8p8FZOZccbiSt3qKiJpkWwJYi2w0cx2An0Evavd3VemPLIcsb89mCgoihIEBAP33be1hY7eo8wsK44kBhHJT8kmiEvSEkUOawoTxPyoEkRtJe6wqaGd8yJoBxGR/JVUgnD33ekKJFdFXYI4/ViP6lYlCBFJqWR7Ut9uZpUx61VmNqXniWhq76GsuIBZ06Kp3pk1rZjja8rVYU5EUi7ZRuqV7n7sm8jdW4FVqQ0pt+xv72XBrGmYRddRbVVtFRsa2nDXiOwikjrJJogCM6saXjGz2STfjpFXmtp7MzLN6FjOWFzJ4e5+Gg73RBqHiOSXZL/c/wd4yMx+GK5fAXwqtSHllv3tvZy9dHakMayuDTvMNbRSOydzI8qKSH5LqgTh7ncArwMOhI/Xufs30xFYLhgccg509EZ2B9OwE+fNYHpJodohRCSlkq4ecvengafTEEvOOdTVx8CQR3YH07CiwgJWLpqlob9FJKWSbYOQGPuP9YHIfC/qkVbVVvHUvg56jw5GHYqI5ImkE4SZfSAdgeSiqPtAxFq1uJKBIeepfe1RhyIieWLcKiYz+0HsKnAG8Om0RZRDmtqDu4ayIUGcURt0T9mwp401S6JtNBeR/JBICaLD3d8QPq4AfjeZE5rZrWbWbGZx57g2swvNrN3MNoaPj03mfOm0v6OXksICZpeXRB0KcyvKWFQ1TQ3VIpIyiTRSj7yN9SOTPOdtwBeAO8bY50/u/qpJniftmtqDO5ii7CQXa1VtFet3HY46DBHJE+OWINx9J4CZVYfrk/oGcvf7gbz4FtvfHv0trrFWLa5kX3vvsQEERUQmI5lG6kyOuXSOmT1hZr8ysxeNtpOZXW9m68xsXUtLSwbDCzS1RzNR0GhW1T4/cJ+IyGQlkyAyVY/yOLDE3U8HPg/8dLQd3f1md6939/qampoMhXfs3Bzo6GVuRWlGzzuWFcfNpKSwQO0QIpISySSIjIwE5+4d7t4VLt8NFA9Xb2WTzr4B+gaGmFuRPSWI0qJCXrRwphKEiKRE1pUgzGy+ha2+ZnYWQYyHMnHuZDR39AEwd2b2lCAgmGFu0942jg4ORR2KiOS4ZBLER8xs8WRPaGbfBR4CTjKzRjO7zsxuMLMbwl1eD2w2syeAm4A3eRaOY93SGSSImhlZliBqK+k9OsSWps6oQxGRHJfwWEzuvsnMngROm8wJ3f3N4zz/BYLbYLNaS1eYILKoDQKeb6jesKeVUxfOijgaEcllyQ618biZnZmWSHLMsRJEliWIhZXTqKkoVTuEiExasqO5ng28xcx2A90E7RLu7itTHlmWa+4MelFHNdXoaMyMVYsr2dCgBCEik5NsgnhFWqLIQS2dfdRUlGZNL+pYq2qr+O3TB2jt7qcqC4YBEZHclGwV01nAYXffDbwN+F9gSo4M19LZR3WWVS8Ne77DnEoRIjJxySaIf3H3TjM7D3gp8HXgK6kPK/u1dPZlVSe5WCsXzaKwwDSBkIhMSrIJYng2mlcCN7v7L4EpWYcxXMWUjaaXFHHy/Aq1Q4jIpCSbIPaa2VeBNwJ3m1npBI6R844ODnH4SH/W9YGItaq2ko172hgayrouJCKSI5L9cn8D8BvgFe7eRtD+8P6UR5XlDnX14559vahjrVpcRWffANtbuqIORURyVFJ3Mbn7EeAnMev7gf2pDirbZWsv6ljDDdWP72ll+byKiKMRkVw05aqHUqGlK5hvIVvbIACWVpdTOb2Y9bvVUC0iE6MEMQHZ2os6lplRv6SKdUoQIjJBSScIM/tAOgLJJcMjuWZzggCor5vNjpZuDoXjRomIJGPcNggz+0HsKnAG8Om0RZQDWrr6mDWtmNKiwqhDGdOZdVUArNvdyiteND/iaEQk1yRSguhw9zeEjyuA36U7qGyXzX0gYp26cBYlRQWs25UXU4CLSIYlkiA+BRAzq9tH0hdObmjO4l7UsUqLCjl90Sy1Q4jIhIybINx9Z7h4a7g+5f8dzZUSBATtEJv3ttPTPzj+ziIiMbJuytFs5+5BgsjiPhCx6pdUcXTQeaJRw26ISHKSSRAaswHo7h+k5+hgzpQg1iwJG6rVDiEiSVIJIknNHUEnuWweZiNW5fQSTpw3g8d2qR1CRJKTTIL4UNqiyCHPD7NRFnEkiVuzZDaP72llUAP3iUgSEk4Q7r45nYHkipau3OgkF+vMuio6ewfYeqAz6lBEJIck1ZPazK4ws4pw+aNm9hMzW52e0LLTcC/qXLjNddiZdcGkf2qHEJFkpGJGuS+nPqzs1dLVR3GhMWtacdShJGxR1TTmzSxVO4SIJEUzyiWppbOP6hmlFBTkTpu9mVFfN5vHdh3GXe0QIpIYzSiXpGyei3osa5fOZn97L3sOH4k6FBHJEZpRLknNOdSLOtbaZXMAeHjHoYgjEZFckVSCcPcj7v4Td38uXN/v7r9NT2jZKZeG2Yh1wtwZVM8o4eEdaqgWkcRMqeqhyRoccg53584wG7HMjLOXzeGh7YfUDiEiCdGEQUk41NXHkEPNzNzpJBdr7bI5NHX0svuQ2iFEZHyaMCgJzcd6UedeCQLgnJh2iLrq8oijEZFspwmDkpCLvahjHV9TTvWMUh5SQ7WIJCDhCYNiTNkJg1pysBd1LDNj7bLZPLxD7RAiMr6EJwwanlFushMGmdmtZtZsZnHHdrLATWa2zcw2ZdNQHrleggA45/g5HOjoY5faIURkHMk0Ut+aonPeBlwyxvOXAsvDx/Vk0VAeLZ19VJQVUVZcGHUoEzbcH+Kh7apmEpGxZXw+CHe/HxirFHI5cIcHHgYqzWxBKs49WbnaizrWsupyaipK1WFORMaVjTPKLQQaYtYbw22Ra+7szenqJQjaIc5ZNoeH1A4hIuPI6RnlzOx6M1tnZutaWlrSfr6gF3Vu9oGIdd4J1bR09rFF80OIyBiycUa5vcDimPVF4ba/4O43u3u9u9fX1NSkPbCWztzsRT3SecurAXjguYMRRyIi2SypGeXM7PdmdlnsdjO7OcUx3QVcFd7NtBZod/f9KT5H0rr7BujuH8yZuajHclzlNE6YO4P7lSBEZAzj9qQeYSnwATM7090/EW6rT+YAZvZd4EKg2swagY8DxQDu/hXgbuAyYBtwBLg2yRjToiXHe1GPdP7yar7zyB56jw7m9F1ZIpI+ySaINuBi4CYz+znw1mRP6O5vHud5B96Z7HHTLR/6QMQ6f3k13/jzLtbtaj1W5SQiEivZwfrM3Qfc/e+BHwMPAHNTH1b2GS5B5EMVE8DZS+dQXGj86bn0N+6LSG5KNkF8ZXjB3W8DrgGmxHwQzR29QP5UMZWXFrFmSRV/UjuEiIwi2QmDvjpifb27vz21IWWnlq4+igqMqun5MwX3+ctreHp/x7HSkYhIrITaIMzs84zRUc7d352yiLJUS2cf1TNKKSjIuu4gE3bB8hr+32+28OdtB3ntqqzoiygiWSTREsQ6YH34eE3M8vAj7+XqXNRjedFxM6maXsz9W9UOISJ/KaEShLvfPrxsZu+NXZ8qWjr7mJejM8mNpqDAuPCkufxxawuDQ05hHpWORGTyJjIn9ZQcwCdfelGPdNHJcznc3c/GhtaoQxGRLDORBDHlDA45h7r78+YW11gXnFhDUYHxu2eaow5FRLJMQgnCzDrNrMPMOoCVw8vD29McY+QOd/czOOR51wYBMGtaMWfWzeYPShAiMkJCCcLdK9x9ZvgoilmucPeZ6Q4yavk2zMZIF58yly0HOmk4rFnmROR5qmJKwPAwG/lYxQRw8SnzAPjDsypFiMjzlCAS8Hwv6vy6i2nY0upyltWU83slCBGJkdRgfWb2f+NsbgfWu/vG1ISUffJtoL54Lj55Lrc/uJuuvgFmlCY7hqOI5KNkSxD1wA0EU4AuBP4WuAT4mpn9c4pjyxotnX1UlBYxrSR/h8W++JR59A8O8YAG7xORULIJYhGw2t3/yd3/CVhDMJrrBQQD9+WlfOxFPVL9kiqqphfz681NUYciIlki2QQxF4gd2e0oMM/de0ZszystnX1U53mCKCos4OUr5vO7Z5rpGxiMOhwRyQLJJohvA4+Y2cfN7OPAn4HvmFk58HTKo8sSB6dACQLg0tPm09U3oLmqRQRIIkGYmQG3AdcTzCzXBtzg7p909253f0t6QoxeS2cfc6dAgnjx8dXMLCvi7idVzSQiSdzF5O5uZne7+2kEo7tOCT39g3T2DUyJEkRJUQEvWzGfe55uon/gNEqKdBe0yFSW7DfA42Z2ZloiyVLNnUEfiLkV+dkHYqRLT51PR+8AD25XNZPIVJdsgjgbeNjMtpvZJjN70sw2pSOwbNE8PBf1FChBAJy3vJoZpUW6m0lEkusoB7wiLVFkseaO/B5mY6Sy4kJeespcfrW5iU9efqqqmUSmsGT/+vcA5wNXu/tugrkh5qU8qiwy1aqYAC5ftZD2nqPcu0VDb4hMZckmiC8B5wBvDtc7gS+mNKIs09zZR1GBUTmtOOpQMub8E6qpnlHCTzfsjToUEYlQ0m0Q7v5OoBfA3VuBkpRHlUWaO4I+EAVTaDrOosICXn36cfz+mWbae45GHY6IRCTZBHHUzAoJpx01sxpgKOVRZZHmzt4p00Ad669WLaR/cIi7n9wfdSgiEpFkE8RNwJ3AXDP7FPAA8O8pjyqLtHT2UTOF2h+GnbZwFstqyrlT1UwiU1ZSCcLdvw38M/AfwH7gte7+w3QEli2aO/umzB1MscyM161ayKM7D2umOZEpKul7GN39WXf/ort/wd2fSUdQ2aJ/YIjD3f1TsooJ4LWrFmIGP1zfGHUoIhKBhBKEmT2ein1yzcHhqUanYBUTwKKq6VywvIYfPNbAwGBeNzWJSByJdpQ7ZZwe0wbMSkE8WaVlivWijufKs2v522+u594tLbxsRV53eRGRERJNECcnsE/eTSJwbJiNKdgGMeyik+cyt6KU7z66RwlCZIpJKEGEvaannKnYi3qk4sIC3lC/mC/9cRt723pYWDkt6pBEJEMiGWjHzC4xsy1mts3MPhjn+WvMrMXMNoaPd0QRZ3NHH2ZQPSOv+wKO641nLsaB7z/WEHUoIpJBGU8QYUe7LwKXAiuAN5vZiji7ft/dzwgft2Q0yFBzZx9zyksoKpzaA9Ytnj2dC0+s4TuP7NF0pCJTyKS++cxsgZklW0F/FrDN3Xe4ez/wPfY2bEAAABCLSURBVODyycSRLi2dvVTPmLrtD7Heft5SDnb1cdfGfVGHIiIZMtl/jb8JPGtm/53EaxYCsXUVjeG2kf46nHPiR2a2ON6BzOx6M1tnZutaWlqSCCExQSe5qdv+EOu8E6o5aV4FX39gJ+4edTgikgGTShDu/lJgGfCN1IRzzM+BOndfCdwD3D7K+W9293p3r6+pqUlxCEEbxFS+xTWWmXHd+Ut5tqmTB7cfijocEcmASVeue+CpJF6yF4gtESwKt8Ue85C794WrtwBrJhdl8oaGnINdShCxXnP6cVTPKOGWP+2IOhQRyYDJtkH8t5l92cxWJvGyx4DlZrbUzEqANwF3jTjugpjV1wAZH9Lj8JF+BoZcCSJGWXEhb1tbx71bWtjS1Bl1OCKSZpMtQdwI/CPwFjN7WSIvcPcB4F3Abwi++H/g7k+Z2SfN7DXhbu82s6fM7Ang3cA1k4wzaQc6gj4Q89QG8QJXnbOE8pJCPv+H56IORUTSLNk5qUd6NUED8zSCqqAlibzI3e8G7h6x7WMxyx8CPjTJ2CalqT1IEAvUMewFqspLuObcOr70x+2850Any+dVRB2SiKRJooP1fdnM/s7MzjOzmTFP9QL3A58maKzOG/uHE8QslSBGesd5y5heXMhNf9gWdSgikkaJVjFtIOjU9ilgp5ntNrOfA2cCy9x9r7vnVQ+qpvZeCgtM/SDiqCov4aoX1/GLTfvY1qy2CJF8lVCCCG8n/Qd3f4m7zwHOBb4EtAOvTGeAUdnf3su8ilIKp9Bc1Mn4m/OXMa24kM/cszXqUEQkTSbUBuHujQQd3H6V2nCyR1NHD/NVvTSq2eUlXH/BMj77u+dYt+sw9XWzow5JRFJsag8yNIb97b0smKUG6rFcf8Ey5s0s5d9++Yx6V4vkISWIONydpvZelSDGMb2kiH96+UlsbGjj55v2Rx2OiKSYEkQcHb0DHOkf1B1MCfjr1Ys4ZcFMPv2rZ+npz6v7FESmPCWIOIb7QKgEMb7CAuPjr17B3rYePvt7NViL5BMliDj2t/cA6gORqLXL5vCG+kXc8qedPL2vI+pwRCRFlCDieL4EoUbqRH34slOonFbMh+58ksEhNViL5AMliDj2t/dihgbqS0Ll9BI+9uoVPNHQxjf+vDPqcEQkBZQg4tjf3kPNjFKKp/hUo8l6zenH8bIV8/ivX2/hqX3tUYcjIpOkb8A49rb1cJwG6UuamfHpv15J5fRi3vO9jbqrSSTHKUHE0XC4h8Wzp0cdRk6aXV7C/7zhdLY1d/Fvv3w66nBEZBKUIEYYHHL2tfWwuEoliIk6f3kNf3vBMr79yB5+sK5h/BeISFZSghihqaOXgSFXCWKS3v+KkzjvhGo+eudm1u9ujTocEZkAJYgRGg4fAWCRShCTUlRYwBeuXMX8WWXc8K317GvriTokEUmSEsQIwwlicZVKEJNVOb2EW66up7d/kLd9/REOd/dHHZKIJEEJYoSG1h7M0F1MKXLivApuubqextYervnGo3T2Ho06JBFJkBLECI2Hj7BgZhklRXprUuXsZXP48ltX8/S+Dt5+22NKEiI5Qt+CIzS29rBIDdQpd9HJ8/jcm1axYU8bV35N1U0iuUAJYoQ9h4+o/SFNXrlyATdftYatBzp5w1cfOtbeIyLZSQkiRlffAE0dvSyrKY86lLx10cnzuP3tZ3Ggo5fLv/hnHt5xKOqQRGQUShAxdrR0AXB8zYyII8lva5fN4WfvPJfK6cW89ZZH+Mafd2rKUpEspAQRY0dLNwAnzFUJIt2W1czgp+88l5ecWMMnfv401972GC2dfVGHJSIxlCBibG/porDAqJ2tBJEJM8uKueXqej55+Yt4aPshXvHZ+/nR+kaVJkSyhBJEjO0tXSyZPV23uGaQmXHVOXX84h/Oo27OdN73wyd441cf5tkmzUwnEjV9E8bY0tTJ8XPV/hCF5fMq+NENL+bTf30aW5s7ufRzf+K939vAzoPdUYcmMmUVRR1AtujqG2DHwW5ec/rCqEOZsgoKjDeeWcvLV8znq/fv4LYHd/LzTft51coFXHvuUs5YXBl1iCJTihJE6Km97bjDykWzog5lyqsqL+GDl57M28+r46v37eD7jzXws437WFVbyZVn1XLJqfOpKCuOOkyRvKcqptCTe4MpMk9dqASRLeZWlPEvr1rBQx+6iI+/egVtR47y/h9tov7ffsc7v/04v968X8N2iKSRShChTY3tzJ9ZRk1FadShyAgVZcVce+5SrnlxHRsa2vjZhr38YtN+fvnkfooKjDVLqrjwpLmcc/wcViyYqZsMRFJECQIYGnIe3H6ItctmRx2KjMHMWF1bxeraKj76qhWs393KfVtb+OOWFj7962cBKCkq4LSFs1hdW8mLjpvF8nkzOL5mBmXFhRFHL5J7IkkQZnYJ8DmgELjF3f9zxPOlwB3AGuAQ8EZ335WueJ7e38HBrj5ecmJNuk4hKVZcWMDaZXNYu2wOH7jkZJo7elm3u5XHd7eyoaGN2x/aTf/AEAAFBkvmlHN8zQwWVU1jUdU0FlZOY2H4s2p6CQUFFvEViWSfjCcIMysEvgi8DGgEHjOzu9w9dob764BWdz/BzN4EfBp4Y7pi+tH6RooLjf9z8tx0nULSbO7MMi47bQGXnbYAgP6BIXYd6mbrgU62HujiuQOdbGvu4sHtBznSP/iC1xYWGFXTi5ldXsLs8hLmlJcyu7yEGWVFzCgNHuWlRcwoLaQ8XC4vKaKkqIDSogJKhh+FwUPJRvJFFCWIs4Bt7r4DwMy+B1wOxCaIy4Ebw+UfAV8wM/M0dLH91ZP7+d5je3jVyuOonqH2h3xRUlTAifMqOHFexQu2uzttR46yt62HvW097Gvr4VBXP4e6+znc3cfh7n6e2d/B4SP9dPYOMDiU/K9ccaEFyaKogOLCAgoLjAIzCgssXObY+rGfBUZhuD1YDrZbTK6xcMXg2HaLee755WOveMF+zy/bsf1esC3OMfNRPl7VtefWsaq2KuXHjSJBLAQaYtYbgbNH28fdB8ysHZgDHIzdycyuB64HqK2tnVAwB7v7OWleBR+69OQJvV5yi5lRVV5CVXnJuHesuTt9A0N09Q3Q3TcQ/hyku2+A7v4B+geGgsdg8LNvxHr/wBBHB4cYHHIG3RkacgY9aPMacmcw5uegB+cbHAoeA0ND9A0EycmB4X+NnOdX/Fic4OGa+4h9w+sgZt/gOX/Bfj7imOTpaCd5elm09aTnbr6cbqR295uBmwHq6+sn9Nm/5axarjyrlkJVC8gIZkZZcSFlxYUqXcqUFMX9gHuBxTHri8JtcfcxsyJgFkFjdcoVhMV+ERF5oSgSxGPAcjNbamYlwJuAu0bscxdwdbj8euAP6Wh/EBGR0WW8iilsU3gX8BuC21xvdfenzOyTwDp3vwv4OvBNM9sGHCZIIiIikkGRtEG4+93A3SO2fSxmuRe4ItNxiYjI8zQmgYiIxKUEISIicSlBiIhIXEoQIiISl+XL3aNm1gLsnuDLqxnRSzuH6VqyT75cB+hastVkrmWJu8cdqTRvEsRkmNk6d6+POo5U0LVkn3y5DtC1ZKt0XYuqmEREJC4lCBERiUsJInBz1AGkkK4l++TLdYCuJVul5VrUBiEiInGpBCEiInEpQYiISFxTPkGY2SVmtsXMtpnZB6OOZzxmtsvMnjSzjWa2Ltw228zuMbPnwp9V4XYzs5vCa9tkZqsjjv1WM2s2s80x25KO3cyuDvd/zsyujneuiK7lRjPbG342G83sspjnPhReyxYze0XM9kh//8xssZnda2ZPm9lTZvaecHvOfS5jXEsufi5lZvaomT0RXssnwu1LzeyRMK7vh1MmYGal4fq28Pm68a4xIe4+ZR8Ew41vB5YBJcATwIqo4xon5l1A9Yht/wV8MFz+IPDpcPky4FcE0/CuBR6JOPYLgNXA5onGDswGdoQ/q8Llqiy5lhuB98XZd0X4u1UKLA1/5wqz4fcPWACsDpcrgK1hvDn3uYxxLbn4uRgwI1wuBh4J3+8fAG8Kt38F+Ltw+e+Br4TLbwK+P9Y1JhrHVC9BnAVsc/cd7t4PfA+4POKYJuJy4PZw+XbgtTHb7/DAw0ClmS2IIkAAd7+fYH6PWMnG/grgHnc/7O6twD3AJemP/oVGuZbRXA58z9373H0nsI3gdy/y3z933+/uj4fLncAzBHPC59znMsa1jCabPxd3965wtTh8OHAR8KNw+8jPZfjz+hFwsZkZo19jQqZ6glgINMSsNzL2L1Q2cOC3ZrbezK4Pt81z9/3hchMwL1zOhetLNvZsv6Z3hVUvtw5Xy5Aj1xJWS6wi+G81pz+XEdcCOfi5mFmhmW0EmgkS7nagzd0H4sR1LObw+XZgDpO8lqmeIHLRee6+GrgUeKeZXRD7pAflypy8dzmXYw99GTgeOAPYD/xPtOEkzsxmAD8G3uvuHbHP5drnEudacvJzcfdBdz8DWETwX//JmY5hqieIvcDimPVF4bas5e57w5/NwJ0EvzgHhquOwp/N4e65cH3Jxp611+TuB8I/6iHgazxflM/qazGzYoIv1G+7+0/CzTn5ucS7llz9XIa5extwL3AOQZXe8EygsXEdizl8fhZwiEley1RPEI8By8M7A0oIGnfuijimUZlZuZlVDC8DLwc2E8Q8fNfI1cDPwuW7gKvCO0/WAu0x1QbZItnYfwO83MyqwqqCl4fbIjeifeevCD4bCK7lTeGdJkuB5cCjZMHvX1hP/XXgGXf/TMxTOfe5jHYtOfq51JhZZbg8DXgZQZvKvcDrw91Gfi7Dn9frgT+EJb/RrjExmWyZz8YHwV0ZWwnq9z4SdTzjxLqM4I6EJ4CnhuMlqGv8PfAc8Dtgtj9/J8QXw2t7EqiPOP7vEhTxjxLUhV43kdiBtxM0tm0Drs2ia/lmGOum8A9zQcz+HwmvZQtwabb8/gHnEVQfbQI2ho/LcvFzGeNacvFzWQlsCGPeDHws3L6M4At+G/BDoDTcXhaubwufXzbeNSby0FAbIiIS11SvYhIRkVEoQYiISFxKECIiEpcShIiIxKUEISIicSlBiIhIXEoQIilgZnVm1hMOJz0nZmjpphFDTZeYWdeI115jZl8Y49jTwtf2m1l1+q9GJFA0/i4ikqDtHoydA8G4P5jZjUCXu//38E5Bh9/EuXsPcIaZ7UpNmCKJUQlCZBRmNsvMDsSsrzezWWk+5w0xpY2dZnZvOs8nMhaVIERG4e7tZjbdzIo8GEL5CYIhEP40yUNPC4dxHjabcKwfd/8K8JVw0Lk/AJ+J83qRjFCCEBlbE8FMZQ0Ewy03peCYPTFVUZjZNUD9iH0+RzDg2s9TcD6RCVGCEBnbPuA4MzsbOOjuz6X7hGHCWAK8K93nEhmLEoTI2PYRjOx5afgTM/s9cJWHc3OkkpmtAd4HnO/B/AUikVGCEBnbPuBK4CJ3P2hmBcAJJD4fdbLeRdAmcW94t9M6d39Hms4lMiYN9y2SBDM7FXi7u//fEdvrgF+4+6lpPPcugvkXDqbrHCKxdJurSBLcffPI5BAaBGaNuDspJYY7ygHFgKqdJGNUghARkbhUghARkbiUIEREJC4lCBERiUsJQkRE4lKCEBGRuJQgREQkLiUIERGJSwlCRETi+v8p6IrtHJop7wAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + } + } + ], + "source": [ + "bb = BlackBody(5000. * u.Kelvin)\n", + "\n", + "nu = np.linspace(1., 3000., 1000) * u.THz\n", + "bb5000K_nu = bb(nu)\n", + "plt.plot(nu, bb5000K_nu)\n", + "plt.xlabel(r'$\\nu$, [{0:latex_inline}]'.format(nu.unit))\n", + "plt.ylabel(r'$I_{\\nu}$, ' + '[{0:latex_inline}]'.format(bb5000K_nu.unit))\n", + "plt.title('Planck function in frequency')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "pw9HsyXWXNbL" + }, + "source": [ + "### Usando $LaTeX$ para rotular eixos\n", + "\n", + "Aqui, usamos a marcação $LaTeX$ para adicionar rótulos nos eixos com melhor aparência. Para fazer isso, colocamos o texto de marcação LaTeX em cifrões, dentro de uma string `r'\\$ ... \\$'`. O `r` antes das aspas indica que a string é \"raw\" e as barras invertidas são tratadas literalmente. Este é o formato sugerido para o texto do rótulo do eixo que inclui marcação. \n", + "\n", + " Agora integramos numericamente usando a regra do trapézio.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "id": "anjDzqDwXNbM", + "outputId": "d90bf5d4-3ae8-41d9-d788-e7700b63ce86", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 38 + } + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/latex": "$1.1280849 \\times 10^{10} \\; \\mathrm{\\frac{erg}{s\\,sr\\,cm^{2}}}$", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "execution_count": 3 + } + ], + "source": [ + "np.trapz(x=nu, y=bb5000K_nu).to('erg s-1 cm-2 sr-1')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "M9xbbIV6XNbM" + }, + "source": [ + "Agora podemos fazer algo semelhante, mas para uma grade de comprimento de onda. Queremos integrar em uma faixa de comprimento de onda equivalente à faixa de frequência que fizemos anteriormente. Podemos transformar a frequência máxima no comprimento de onda (mínimo) correspondente usando o método `.to() `, com a adição de uma *equivalência*." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "id": "dkNoQ4E0XNbN", + "outputId": "17627d50-cd87-4cfb-d282-95bcdea27e97", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 303 + } + }, + "outputs": [ + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEeCAYAAABYEGiuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5xUd73/8dd7tgNL2V0glIUlQAqphIUUEo0xmmJMLGnEbjS/eM1VrxqjV6/12u/Va2I0xnKtId2IisaouTGRECCBBEIaZRcWCGXou2z//P44Z2BYli2zM3tmZz/Px2Mee9qc85mz5bPfcr5fmRnOOedcJsSiDsA551zu8iTjnHMuYzzJOOecyxhPMs455zLGk4xzzrmM8STjnHMuYzzJuG5J+j9JH8iWa0g6XtIKSfskfSSTcXW47iRJ+yXlZeDc+yUdm+7zZpqk8yXVRXTtL0r6dRTXdj3nScYBIKlG0oHwj91WST+XNCzquI7iU8CjZlZqZrdm6iLhPbkwsW5mG8xsmJm1pfta4XnXpfu8uSLKZOb6xpOMS/ZmMxsGnAFUA5+LOJ6jmQw8H3UQzrnueZJxRzCzTcCfgJM77pM0VdLfJcUl7ZD0G0kjk/bXSPqkpOck7ZF0j6TipP1XhFVdeyWtlXRxJ9cYF77/5k72/R14HfD9sNR1XMeqNknvlfRE0rpJulHSK5J2S7pdkpL2f1DSC2H122pJZ0j6FTAJ+H14nU9JqgrPlR++b7ykBZJ2Sloj6YNJ5/yipHsl/TI87/OSqo92z8PzTguXfx7G+MfwvU9JmnqU9/1C0ifC5QnheT6c9L3aKSkmaZSkP0jaLmlXuDwxPO4aScs6nPffJC0Il4sk/ZekDWEp9w5JJUeJZ7ykB8LrrE+uzuzunoT3fXm4777wZ+c/JQ0l+HkcH34v9ksaH76tsKf32EXDk4w7gqRK4FJgeWe7ga8D44ETgUrgix2OuRq4GJgCnAq8NzzvHOCXwM3ASOA1QE2Ha08BHgO+b2bf7nhxM7sAeBy4KaxiermHH+syYHYYz9XAReH1rgrjfzcwHLgciJvZu4ANhKU7M/tWJ+e8G6gL78WVwNckXZC0//LwmJHAAuD7PYwV4FrgS8AoYA3w1aMc9xhwfrj8WmAdwX1NrD9uZu0Ev+v/S1AKnAQcSIrn98DxkqYnnfc64K5w+RvAccDpwDRgAvD5joFIioXnejY85vXAxyRdlHRYp/dEUiHwW+DnQBkwH3grgJnVA5cAm8PvxTAz29zV+Vz2yMkkI+lnkrZJWtXD468O/4N9XtJd3b8jZz0kaTfwBMEfr691PMDM1pjZI2bWZGbbge8Q/DFLdquZbTaznQR/dE4Pt18P/Cx8f7uZbTKzF5PeNwN4FPiCmd2Z5s/2DTPbbWYbwmskYvoA8C0zW2qBNWZW293JwkQ8F7jFzBrNbAXwE4JklfCEmS0M23B+BZzWi3h/a2ZLzKwV+E1SvB09Bpwb/oF/DfCtMC4Ivi+PAZhZ3MweMLMGM9tHkLReG+5rAH4HzAs/23TgBGBBWOK7Afg3M9sZvvdrBEmwo9nAaDP7spk1h21MP+5w7NHuyVlAPsHPTouZPQgs6cF96ss9dv0gJ5MMwX9DR1TDdCb8hfoMMNfMTgI+lsG4st1bzGykmU02s38xswMdD5A0VtLdkjZJ2gv8GqjocNirScsNQKIDQSWwtovrvwPYBNyf+kc4qlRjOprxQOKPbkItwX/wR7tmcaKqrQeOFu9hzGwtUE+QhM4D/gBslnQ8SUlG0hBJP5JUG37f/gGM1KGecncRJhmCUsxDYfIZDQwBng6rGncDfw63dzSZoEprd9Kx/w6M7eJzJe7JeGCTHT5i78aj3ZwenM9liZxMMmb2D2Bn8rawfvrPkp6W9LikE8JdHwRuN7Nd4Xu39XO4A83XAANOMbPhwDsJqtB6YiPQadtC6IvADuAu9a6bcD3BH8KEY3rx3q5i6mqI8s1AmaTSpG2TCJJkf3uMoLquMGxPewx4D0FV24rwmE8AxwNnht+3RJVa4nv3CDBa0ukEySZRot9BULV2UvgPyEgzGxF2EOloI7A+6biRYQ/AS3vwGbYAE8KSU0Jl0rIPFz9A5WSSOYo7gX81s1nAJ4EfhNuPA46T9E9Ji9VJQ7Q7TCmwH9gjaQJB+0pP/RR4n6TXh43RE5KSPUALcBUwFPhlWAXUEyuAt4X/rU8jqJbrqZ8An5Q0S4FpkiaH+7YCnT67YmYbgUXA1yUVSzo1vG4Uz208BtxEUDoB+L9w/Ymk7talBMlit6Qy4AvJJzCzFuA+4NsEbSKPhNvbCaq8vitpDBzsYJDczpKwBNgn6RZJJZLyJJ0saXYPPsOTQBtwk6R8SVcAc5L2bwXKJY3owblcFhkUSUbB8x7nAPdJWgH8CBgX7s4HphM0ns4Dfqyk3lLuCF8i6OK8B/gj8GBP32hmS4D3Ad8N3/8YQRVL8jHNwNsIqlh+1sNE812gmeAP0S8I2jB6GtN9BO0TdwH7gIcI/shC0MHhc2HVzyc7efs8oIqgVPNbgrakv/b02mn0GEESSSSZJwhKdv9IOuZ/gBKCksligiqvju4CLgTuC9uCEm4h6HywOKxq+ytBqegwYUK7jKDqbn14rZ8A3SaGpO/79cBughLyH4CmcP+LBJ0B1oXfj/FHO5fLLsrVScskVQF/MLOTJQ0HXjKzcZ0cdwfwlJn9b7j+N+DTZra0P+N1zh1O0lPAHYnfTTcwDYqSjJntBdaH3VUJq0USvVAeIuwCKqmCoPrMn7x2rp9Jeq2kY8LqsvcQdDfvrMTlBpCcTDKS5hPU8R4vqU7S9QQ9l66X9CzB0+JXhIc/DMQlrSbo2nqzmcWjiNu5Qe54gmdsdhN0VLjSzLZEG5Lrq5ytLnPOORe9nCzJOOecyw6eZJxzzmVMzj0ZW1FRYVVVVVGH4ZxzA8bTTz+9w8w6G8Whz3IuyVRVVbFs2bLuD3TOOQeApG7H60uVV5c555zLmMiSTHcjJUt6h4I5RVZKWpT0XItzzrkBIsqSzM/peqTk9cBrzewU4CsEY48555wbQCJrkzGzf4RDvxxt/6Kk1cXAxEzH5JxzLr0GSpvM9QTTrzrnnBtAsr53maTXESSZc7s45gaC2fuYNGlSP0XmnHOuO1ldkgnn6PgJcEVX44mZ2Z1mVm1m1aNHZ6Srt3POuRRkbZKRNIlgrpJ3mdnLUceTDvVNrdxy/3Ns3NkQdSjOOdcvIqsuC0dKPh+okFRHMFNfAYCZ3QF8HigHfhDOyNpqZtXRRJseD63YxD3LNjJuZDEfu/C4qMNxzrmMi7J32bxu9n8A+EA/hdMv5i/ZAMCiNXE+dmHEwTjnXD/I2uqyXLOybg+rNu1l7PAilm/cRUNza/dvcs65Ac6TTD+5a0ktxQUx/uOyGbS0GUtrdkUdknPOZZwnmX6wv6mV363YzJtPHc8FJ4yhIE8sWrsj6rCccy7jPMn0gwUrNtPQ3Ma8MycxpDCfmZWjWLTGZ3h2zuU+TzL9YP6SDZxwTCkzK0cCcPbUclZt3sOehpaII3POuczyJJNhK+v2sHLTHubNmUTYFZu50yowgyfXeWnGOZfbPMlk2PylGyjKj/GWmRMObju9ciQlBXneLuOcy3meZDKovqmV3y3fxGWnjmdEScHB7YX5MWZPKWPRWi/JOOdymyeZDPr9s5upb27jujMrj9g3d2o5a7btZ+vexggic865/uFJJoPmL9nAcWOHccakUUfsO2dqBQBPemnGOZfDPMlkyKpNe3i27vAG/2Qzxg9nREkB/1zj7TLOudzlSSZD7g4b/N+a1OCfLC8mzj62nEVr45hZP0fnnHP9w5NMBjQ0t/LQ8s286ZRxjBxSeNTjzplWzqbdB9jgQ/8753KUJ5kM+MOzW9jf1Mp1Z3Y9S2eiXeaf/vS/cy5HeZLJgLuWbGD6mGHMmnxkg3+yqaOHMqa0yJ+Xcc7lLE8yabZ6815WbNx91Ab/ZJKYO62CJ9fGaW/3dhnnXO7xJJNmdy/dQGF+jLed0XmDf0dnTy0nXt/My9v2ZTgy55zrf55k0uhAcxu/fWZTtw3+yeZO83YZ51zu8iSTRn94bjP7mlqZN6frBv9kE0aWUFU+hCe9XcY5l4M8yaTR/CUbmDp6KLOrum7w7+jsqRU8tW4nrW3tGYrMOeei4UkmTV58dS/PbOhZg39Hc6eVs6+plec27clQdM45Fw1PMmly95KNFObFePsZE3v93rOPLQd8HDPnXO7xJJMGB5rbePCZOi455RhGDe1Zg3+y8mFFnHBMqY9j5pzLOZ5k0mDhyi3sbexdg39H50ytYFntLhpb2tIYmXPORSuyJCPpZ5K2SVp1lP2SdKukNZKek3RGf8fYU/OXbODYiqGcOaUs5XPMnVZOc2s7z9TuSmNkzjkXrShLMj8HLu5i/yXA9PB1A/DDfoip117euo9ltbtSavBPNmdKGXkx+WyZzrmcElmSMbN/ADu7OOQK4JcWWAyMlDSuf6LruflLNgQN/rN63+CfrLS4gFMnjuCf/ryMcy6HZHObzARgY9J6XbgtazS2tPHgM5u46ORjKEuhwb+juVMreK5uD/saW9IQnXPORS+bk0yPSbpB0jJJy7Zv395v1/3Tqi3sOdDCvDmVaTnfOVPLaWs3lqzvqoDnnHMDRzYnmU1A8l/vieG2I5jZnWZWbWbVo0eP7pfgAOY/tZGq8iEHn3PpqzMmj6IoP+bjmDnnckY2J5kFwLvDXmZnAXvMbEvUQSWs2baPJTU7+9zgn6y4II/qqlE+v4xzLmdE2YV5PvAkcLykOknXS7pR0o3hIQuBdcAa4MfAv0QUaqfmL9lIQZ763ODf0TlTK3jx1X3s2N+U1vM651wU8qO6sJnN62a/AR/up3B6pbGljQeeqeONJx1DxbCitJ77nKlB1dvidXEuO3V8Ws/tnHP9LZury7LWn1e9yu6GFq7rwxP+R3PKhBGUFuV7u4xzLid4kknBXUs2MDmNDf7J8vNinHlsmbfLOOdygieZXlqzbT9L1u/k2tmTiMXS0+Df0dlTK6iNN1C3qyEj53fOuf7SqzYZST0ZnKvdzHanGE/Wu3vJBvJj4so0N/gnmzstKCEtWhvn6uohGbuOc85lWm8b/jeHr67+hc8D0t9YkQUONfiPZXRpehv8kx0/tpTyoYU8uTbO1dXpedDTOeei0Nsk84KZzezqAEnL+xBPVnv4+VfZ1dDSpyH9e0ISZ08t559rdmBmaXsOxznn+ltv22TOTtMxA9L8JRuoLCth7tSKjF9r7rQKtu1rYu32+oxfyznnMqVXScbMGtNxzEC0bvt+Fq/LbIN/ssTzMt7LzDk3kPU4yUh6g6QfSzo9XL8hc2Fln7uXbiQ/Jq6qzlyDf7JJZUOYMLLEp2R2zg1ovSnJvB+4GXinpAuA0zMTUvZpam3j/qfruPDEsYwpLe6Xa0rinKnlLF63k7Z265drOudcuvUmyewzs91m9kngjcDsDMWUdf7y/FZ21jdz3Zn922lu7rQK9hxoYfXmvf16XeecS5feJJk/JhbM7NPAL9MfTnaav2QDE0eVcO60zDf4J/N2GefcQNfjJGNmvwOQVBGu35apoLLJ+h31LFobZ96c/mnwTzZmeDHTxgzjn2t9HDPn3MCUyrAyP0t7FFns7qUbyIuJqzL4hH9X5k4tZ+n6nTS3tkdyfeec64tUksygeTKwubWd+5fVceGJYxgzvH8a/Ds6e2oFB1raWLExZ0fqcc7lsFSSzKDp6vTI6q3E65sz/oR/V84+thwJ78rsnBuQvCTThYWrtjB2eBHnTR8dWQwjhhRw8vgRPOntMs65ASiVJPOZtEeRpdZvr2fGuOHk9XODf0fnTCtn+cZdNDS3RhqHc871Vq+TjJmtykQg2cbMqI3XM7l8aNShMHdqBS1txtKaXVGH4pxzvZLSpGWSqiX9VtIzkp6TtFLSc+kOLko79jdT39xGVXn087lUV42iIE8s8nYZ59wA09uh/hN+QzDEzEogJ/vW1saD0Y8nV0RfkhlSmM/MSaNY5O0yzrkBJtXpl7eb2QIzW29mtYlXWiOLWG08mPq4KguqyyB4+n/V5j3sbmiOOhTnnOuxVJPMFyT9RNI8SW9LvNIaWcRq4/XkxcSEkSVRhwIE45iZweJ1Xppxzg0cqVaXvQ84ASjgUHWZAQ+mI6hsUBNvYPzIYgrzU83D6XXaxJEMKcxj0do4F588LupwnHOuR1JNMrPN7Pi0RpJlauP1WVNVBlCYH2N2VZk/lOmcG1BS/Td9kaQZfbmwpIslvSRpjaRPd7J/kqRHJS0Pe7Bd2pfr9VZNvIHJWdCzLNncaeWs3V7P1r05Ofmocy4HpZpkzgJWhEmi112YJeUBtwOXADOAeZ0krc8B95rZTOBa4Acpxtpruxua2XOgJatKMgDnTA2mGvCh/51zA0Wq1WUX9/G6c4A1ZrYOQNLdwBXA6qRjDBgeLo8ANvfxmj1WE/Ysy4YHMZPNGDecESUFLFoT560zoxkV2jnneiPVksyXgT1JXZf3Al/oxfsnABuT1uvCbcm+SDDVcx2wEPjXFGPttcQzMtnwIGayWEycfWw5i9bGMRs045Q65wawVJPMqWZ2cOx5M9sFzExPSAfNA35uZhOBS4FfSeo0Xkk3SFomadn27dv7fOHaeAMSVJZlV5KBoF1m0+4DB5/jcc65bJZqkolJGpVYkVRG76reNgGVSesTw23JrgfuBTCzJ4FioNP5j83sTjOrNrPq0aP7PmJyTbyeccOLKS7I6/O50u2caYl2GX9exjmX/VJNMv8NPCnpK5K+AiwCvtWL9y8FpkuaIqmQoGF/QYdjNgCvB5B0IkGS6XsxpQdq4w1MyrKqsoRjK4YydngR//TGf+fcAJBSkjGzXwJvA7aGr7eZ2a968f5W4CbgYeAFgl5kz0v6sqTLw8M+AXxQ0rPAfOC91k8NEdn2jEwyScydWsHitXHa271dxjmX3VLtXYaZrebw3mC9ff9Cggb95G2f73D+uameP1X7GlvYsb8563qWJTt7ajkPLt/ES1v3ceK44d2/wTnnIpIdY6ZkkUMDY2ZndRl4u4xzbuDwJNNBbZY+I5NswsgSqsqH+Pwyzrmsl3KSkXRLOgPJFjWJeWSyuCQDQWnmqfU7aW3Lyel8nHM5osdJRtK9Sa/7gA9kMK7IbIg3MLq0iKFFKTdX9Yu5UyvY39TKc5v2RB2Kc84dVW/+ku41s4OJRdIPMxBP5Gri9VndHpNw1rFlACxas4MzJo3q5mjnnItGb6rLvgogKfFA5GfTH070auMNWd0ek1A+rIgTxw33xn/nXFbrcZIxs/Xh4s/C9Z0ZiShCB5rbeHVvI5OzcDiZzpwztZxltbtobGmLOhTnnOtUKg3/SnsUWWLDzrBnWUX2l2QgGMesubWdZ2p3RR2Kc851KpUkk7OPmddk6ejLRzNnSjl5MfkQM865rOUlmSSJIf4nlw2MksywonxOmziCJ17xJOOcy06pJJnPpD2KLFETb2DUkAJGDCmIOpQeu+ikY3i2bg9rtu2POhTnnDtCKknmeUmV3R828NTG6wdEz7JkbztjIvkxcd+yjd0f7Jxz/azXSSYcCXlhtwcOQLXxhgHTHpMwurSIC04YwwPP1NHiT/8757JMqsPKPCNpdlojiVhTaxubdx8YcCUZgGtmV7JjfzN/f3Fb1KE459xhUk0yZxJMWrZW0nOSVkp6Lp2B9be6XQdot+wfs6wzrz1uNGNKi7h3qVeZOeeyS6oDdF2U1iiywMGeZQOwJJOfF+PKWRO547G1vLqnkWNGFEcdknPOAamXZOYAO82sFngX8F2gLG1RRaBmR/bPI9OVq6sraTd44Jm6qENxzrmDUk0y/2Fm+ySdC1wI/BS4I31h9b/aeD2lRfmUDS2MOpSUVFUM5cwpZdy7bKNPy+ycyxqpJpnEYFlvAu40sz8CA/Ovc6gm3sDkiiFIA/dZ02tmV1Ibb+Cp9Tk3rJxzboBKNclskvQj4BpgoaSiPpwrKwzEZ2Q6uuTkcZQW5XOvPzPjnMsSqSaGq4GHgYvMbDdBe8zNaYuqn7W2tVO368CAbY9JKCnM4/LTx7Nw5Rb2HGiJOhznnEstyZhZg5k9aGavhOtbzOwv6Q2t/2ze3Uhruw34kgwEVWZNre0seHZz1KE459zAruJKl0OjLw/8JHPKhBGcOG64PzPjnMsKKScZSbekM5AoHXpGZmBXlwFI4prqiazctIfVm/dGHY5zbpDrcZKRdG/S6z7gA325sKSLJb0kaY2kTx/lmKslrZb0vKS7+nK9rtTEGyguiDGmtChTl+hXb5k5gcL8mHcAcM5Frjclmb1mdnX4ugr4a6oXlZQH3A5cAswA5kma0eGY6QTTCsw1s5OAj6V6ve7UxuupKh86oLsvJxs5pJCLTjqG3y7f5FMzO+ci1Zsk89UO65/tw3XnAGvMbJ2ZNQN3A1d0OOaDwO1mtgvAzDI2+mNNvCEnqsqSXVNdyZ4DLfxl9daoQ3HODWI9TjJmth5AUkW43pcn/iYAyXU5deG2ZMcBx0n6p6TFki7uw/WOqr3d2LCzISca/ZOdM7WciaNKvAOAcy5SqTT8/yztUXQuH5gOnA/MA34saWRnB0q6QdIyScu2b9/eq4u8ureR5tb2nOi+nCwWE1fNquSJNTvYuLMh6nCcc4NUKkkmHQ0Xm4Dk2TUnhtuS1QELzKwlLEW9TJB0jmBmd5pZtZlVjx49uleBHOq+nFvVZQBXVk9EwmfNdM5FJpUkk47RF5cC0yVNkVQIXAss6HDMQwSlmEQV3XHAujRc+zC18eC//MkVuVWSAZgwsoTzpo/mvqfraPNBM51zEYikJGNmrcBNBEPTvADca2bPS/qypMvDwx4G4pJWA48CN5tZvK/X7qgmXk9hXoxjhufmHCzXVFeyZU8jj7/Su2pE55xLh1QmLftMOi5sZguBhR22fT5p2YCPh6+Mqd3RQGVZCXmx3Oi+3NGFM8YwakgB9y7byPnHj4k6HOfcINPrkoyZrZJ0laRSAEmfk/SgpDPSH17m1YTPyOSqovw83jpzIo+s3kp8f1PU4TjnBpl0Tlr2w/SF1T/MjNp4Q871LOvomtmVtLQZv13esW+Fc85l1qCetGz7viYOtLRRVZF7PcuSHX9MKadVjuTeZRsJaiGdc65/DOpJy2oSPctyvCQDQQeAl7fuZ8XG3VGH4pwbRAb1pGW1OfyMTEdvPm0cJQV5Pmimc65fDepJy2rjDeTHxISRJVGHknGlxQW86dRx/P7ZLTQ0t0YdjnNukBhwVVzpVBOvZ+KoEvLzBsdtuGZ2JfubWvnjc1uiDsU5N0gMjr+uR1Ebb2DSIGiPSaiePIpjRw/1KjPnXL8ZtDNjmln4jEzut8ckSOLq6kqW1uxi7fb9UYfjnBsEIpsZM2q7GlrY19g6KHqWJXvbGRPIi8lLM865fhHJzJjZIJdHX+7KmNJiLjhhDA88vYmWtvaow3HO5bioZsaMXKL78mAryUDwzMyO/U38/cWMTTbqnHNAdDNjRq423oAElWW53325o/OPH82Y0iKfNdM5l3HZPDNmRtXGGxg/ooSi/LyoQ+l3+Xkx3j5rIo++tI2texujDsc5l8OimhkzcjXx+pwfs6wrV1dX0m5w/9N1UYfinMthUc2MGbnaeAOTygZfe0zClIqhzJlSxn0+aKZzLoMGZUlmz4EWdtY3D7qeZR1dU11JTbyBp9YP6OY151wWSyXJpGVmzChtGESjL3fl0lPGUVqU7x0AnHMZk9LMmB23STpX0u3pCSnzDj4jM4jbZABKCvN48+njWbhqC3sbW6IOxzmXg7pMMpLyJc2SdEQXLEkzJX1bUg3wFeDFDMWYdolnZCaVDe4kA0GVWWNLOwtWbI46FOdcDuquJHMzMB24BUDScZK+IOlF4DZgAyAze52Z3ZbZUNOnNt7A2OFFDCnMjzqUyJ06cQQnHFPqw8w45zKiuyTzLeB54Nvh+ovABcBlZnZumFjajvbmbFUbbxj07TEJkrhmdiXP1e3hhS17ow7HOZdjukwyZtZmZivNLFFh/zZgC/CopB9Lej0DsLfZYBt9uTtvOX0ChXkx7vEOAM65NOtVw7+ZPWRm1wInAI8CHwPGSPqhpDdmIsB0a2huZdu+Ji/JJBk1tJA3njSWh1Zsoql1wBVMnXNZLNXpl+vN7C4zezMwEVhO2G6T7WoPdl/2kkyya2ZXsruhhb88vzXqUJxzOaTPM2Oa2S4zu9PMXt+b90m6WNJLktZI+nQXx71dkkmq7muscKhnWZWXZA4zd2oFE0aWeAcA51xaRTL9ctgl+nbgEmAGME/SjE6OKwU+CjyVrmvXhCWZSV6SOUwsJq6qnsjjr+xg486GqMNxzuWIlJKMpI938rpe0uk9PMUcYI2ZrTOzZuBu4IpOjvsK8E0gbUMF18brKR9ayPDignSdMmdcVV2JBPf5oJnOuTRJtSRTDdwITAhf/w+4GPixpE/14P0TgOR6mbpw20GSzgAqzeyPKcbYqaD7spdiOjNhZAnnTqvg/mUbaWv3QTOdc32XapKZCJxhZp8ws08As4AxwGuA9/Y1KEkx4DvAJ3p4/A2Slklatn379i6PrY03eHtMF66ZXcnmPY38edWrUYfinMsBqSaZMUBT0noLMNbMDnTYfjSbgMqk9YnhtoRS4GTg/8Jha84CFhyt8T/seFBtZtWjR48+6kUbW9rYvOeAd1/uwhtnHMOMccP53EMreXWPT2jmnOubVJPMb4CnwiFmvgD8E7hL0lBgdQ/evxSYLmmKpELgWmBBYqeZ7TGzCjOrMrMqYDFwuZktSzFeAOp2NWDmA2N2pTA/xm3XzaSxpZ1/u2eFV5s55/qk10lGkoCfAzcAu8PXjWb25fD5mXd0dw4zawVuAh4GXgDuNbPnJX1Z0uW9jamnanaEPct8YMwuTR09jC9dcRJProtzx2Nrow7HOTeA9XqESDMzSQvN7BQg5ZKFmS0EFnbY9vmjHHt+qtdJVuPPyPTYVbOC7szfeeRlzjq2nFmTR0UdknNuAEq1uuwZSbPTGkk/qI03MLw4n5FDvHWMjTUAABV3SURBVPtydyTx1beezPiRxXxk/nL2HPD5ZpxzvZdqkjkTWCxpraTnJK2U9Fw6A8uEmng9VRVDCWr8XHeGFxdw67Uz2bq3kX9/cCVm3j7jnOudVCdUuSitUfST2ngDp1WOjDqMAWXmpFF8/I3H8a0/v8RrllVwzexJUYfknBtAUi3JbADOA95jZrWAAWPTFlUGtLS1s2n3AR/iPwU3vmYqc6eV88UFq1mzbV/U4TjnBpBUk8wPgLOBeeH6PoKxyLLWpl0HaGs3f0YmBbGY+M7Vp1NSmMdNdy2nscWnA3DO9UzKbTJm9mHCMcXMbBdQmLaoMuBQzzIvyaRi7PBi/vuq03jx1X18408vRh2Oc26ASDXJtIQjKRuApNFAe9qiyoBD88h4SSZVrzthDO+fO4WfL6rhr6t93hnnXPdSTTK3Ar8lmBXzq8ATwNfSFlUG1MTrGVKYR8WwrC5wZb1bLjmeGeOGc/P9z/qwM865bqU6M+ZvgE8BXwe2AG8xs/vSGVi6BaMve/flvirKz+O262bS1OrDzjjnupfypGVm9qKZ3W5m3zezF9IZVCbUxOu9PSZNpo4expcuD4ad+eH/rYk6HOdcFutVkpH0TDqO6W9t7cbGnQ3eHpNGV86ayOWnjee7f32Fp2t3Rh2Ocy5L9fZhzBO7ebJfwIg+xJMRW/YcoKXNvCSTRpL4z7eezPKNu/jI/BUs/Oh5jCjx4Xqcc4frbZI5oQfHZN1DFN6zLDMSw85cdceT/PuDK/n+dTO9zcs5d5heJZnw6f4B5+AzMj6PTNolDztz3tIKrp3jw8445w5JueF/IKmNN1CUH2NsaXHUoeSkG18zlXOnVfDF3z/vw8445w4zKJJMzY56JpUNIRbzqpxMCIadOY0hhfk+7Ixz7jBpSTKSxkkqSse5MiHxjIzLnDE+7IxzrhPpKsn8CnhR0n+l6Xxp095u1O70Z2T6gw8745zrKC1JxswuBI4F/jcd50unbfuaaGxpZ3KFl2T6wy2XHM9J433YGedcIG1tMhZ4Pl3nS5daH325XxXl53HbvGDYmY/ds9yHnXFukEtrw7+kKklvTec5+yrxjEyVt8n0m2PDYWcWr9vpw844N8j1OclIukXSI5IWAx8CsqoDQE28noI8MW6Ed1/uTz7sjHMO0lOSaQPeDPwDWG9md6fhnGlTG2+gctQQ8vMGRW/trJEYdmb8yGI+Mn8Few60RB2Scy4CvR0g84eSPiTpXEnDAczsv8ys0cw+BayV9KOMRJqimng9k709JhKJYWe27m3ko3cvZ1+jJxrnBpve/nu/HJgBfBVYL6lW0u8lfU3StWb2CHBLT04k6WJJL0laI+nTnez/uKTVkp6T9DdJk3sZK2bmz8hEbOakUXzpipN4/JUdXHbbEzy7cXfUITnn+lGvkoyZ3Wlm/2pmrzWzcmAu8ANgD/Cm8Jhu/4qEUzffDlxCkLTmSZrR4bDlQLWZnQrcD3yrN7ECxOub2d/U6iWZiL3jzMncc8NZtLYZb//hIn78j3W0e68z5waFPjVUmFmdmf3JzL5pZu/qxVvnAGvMbJ2ZNQN3A1d0OPejZtYQri4GJvY2vkPdl70kE7XqqjIWfuQ8Xn/iGL668AXe/4ul7NjfFHVYzrkMi6o1fAKwMWm9Ltx2NNcDf+rtRWp2JIb495JMNhgxpIA73jmLr7zlZBatjXPp9x5n0ZodUYflnMugrO9yJemdQDXw7S6OuUHSMknLtm/ffnB77c4GYoKJozzJZAtJvOusyfzuw3MZXlLAO376FN9++EVa29qjDs05lwFRJZlNQGXS+sRw22EkXQh8FrjczI5atxK2FVWbWfXo0aMPbq+N1zNhVAmF+VmfSwedE8cNZ8FNc7l6ViW3P7qWa+5cTN2uhu7f6JwbUKL667sUmC5piqRC4FpgQfIBkmYCPyJIMNtSuUhNvMHbY7LYkMJ8vnnlqdw6byYvvbqPS7/3OH9etSXqsJxzaRRJkjGzVuAm4GHgBeBeM3te0pclXR4e9m1gGHCfpBWSFhzldEdV68/IDAiXnzaehR85jykVQ7nx18/wuYdW+pw0zuWIXk2/nE5mthBY2GHb55OWL+zL+Xc3NLO7oYXJZV6SGQgmlQ/hvhvP4b//8hI/+sc6ltXs4vvXzWTamNKoQ3PO9UHONlYkBsb0kszAUZgf4zOXnsjP3zeb7fuauOy2J7hn6QbM/Jka5waqnE0yNYlnZHwemQHn/OPH8KePnsesyaO45YGV/Ov85ez1IWmcG5ByNskkSjKTyrwkMxCNGV7Mr95/JjdfdDx/WvUqb7r1cVb4kDTODTg5nWTGjSimuCAv6lBcimIx8eHXTePe/3c27e1w5Q8X8aPH1vqQNM4NIDmcZLxnWa6YNXkUCz96Hm88aSxf/9OLvPfnPiSNcwNFziYZf0Ymt4woKeD2687ga289hafWxbnke4/z0yfWs6fB22qcy2Y5mWT2N7WyY3+TD/GfYyRx3ZmTWHDTuUwuG8JX/rCaM7/+V26+71lWbNztvdCcy0KRPSeTSYnRl726LDcdf0wp93/oHFZv3stvnqrloeWbuO/pOk4aP5x3njWZy08bz9CinPzRdm7AycmSjD8jMzjMGD+cr771FJ767IX851tOpq3d+MyDKznza3/jPx5axUuv7os6ROcGvZz8d6/mYEnGq8sGg2FF+bzzrMm848xJPLNhN79ZXMs9yzbyq8W1VE8exTvPmszFJx/jPQ2di0BOJpnaHQ1UDCtimFeZDCqSmDV5FLMmj+I/LpvBA8/U8ZunNvCxe1Yw6vcFXF1dybw5k/wBXef6kXKtsbS6utqm3XAbrW3G/R86J+pwXMTa240n18X59eJa/rJ6K23txnnTK3jHmZO58MQx5OflZI2xc70i6Wkzq87EuXPyX/3aeAPnTK2IOgyXBWIxMXdaBXOnVbB1byP3LN3I/CUbuPHXTzN2eBHXzp7EtXMqGTeiJOpQnctJOZdkzGDLnkaqvNHfdTB2eDEfef10/uX8qTz60nZ+81Qtt/79FW77+yu8/sSxvP2MiZw5pYxRQwujDtW5nJFzSaapNZjGd7LXu7ujyM+L8YYZY3nDjLFsiDcwf+kG7l26kUdWbwVg+phhzJ5SxpyqMmZPKWPCSC/lOJeqnEsyza3BZFdeknE9Mal8CLdcfAL/duFxPFu3myXrd7K0Zie/X7GZu57aAMCEkSXMrhp1MPFMGzMMSRFH7tzAkHtJpi0syfhkZa4XCvNjzK4qY3ZVGQBt7caLr+5l6fqdLKnZyRNr4jy0YjMAo4YUUF11qKRz0vjhFHgHAuc6lXNJpqm1nbFDChgxpCDqUNwAlhcTJ40fwUnjR/DeuVMwM2riDQeTztKanQer14YU5jFz0khmh4ln5qRRlBT6MznOQQ4mmebWdn8I06WdJKZUDGVKxVCunl0JwNa9jSyt2Rkmnl1872+vYAb5MXHyhBHMmVLG9DHDOHb0UKrKh1I2tNCr2dygk3NJpqm13dtjXL8YO7yYy04dz2Wnjgdgz4EWnq7dyZL1u1has5P//ed6WtoOPYdWWpzPsRVDqaoIkk4i+VRVDGVEiZe8XW7KuSTT0uYlGReNESUFXHDCWC44YSwQ/CzW7TrA+h37Wb+jgZod9azfUc+yml0seHYzyc9Blw8tpCosKU0Jk9CUiqFUVQxhSGHO/Zq6QSQnf3q9JOOyQUFe7GDS6KixpY0NOxtYv6P+YPJZv6Oex1/Zzv1P1x127NjhRQfPM3HUEEaXFgWvYcHX8qGFPnKBy1o5mWS8JOOyXXFBHseNLeW4saVH7KtvaqUmXk/NjoZDpaB4PQ8/v5Wd9c1HHC9B2ZBCRpcWUREmnkQSqigtZPSw4oPbRpYUEIt5u5DrPzmZZLwk4wayoUX5B3u2dXSguY0d+5vYvr+J7fsOvXYk1vc3UVNTz/Z9TQcfTE6WFxMVww4lpIphQeIZXlJAaXE+w4s7LudTWlxAaVG+JyeXkpxLMjGJMh8WxOWoksI8KsuGUFnW9T9SZsb+ptZDiWh/Ezv2HZ6cduxv5sUt+9jb2EJDc1uX55OCKRWGF4cJqKSA4UdJSkOL8hlSmEdJQT4lhXnhct7B5eL8PE9Yg0hkSUbSxcD3gDzgJ2b2jQ77i4BfArOAOHCNmdV0d97C/Jh3E3WDnqSgBFJcwLGjh3V7fEtbO/sbW9nb2MLeA8HXfUnLextb2XugJdweLG/a3cgLB/axr7GFfU2t9GZA9+KCGEMK8w9PPgXB18OX8ykuyKMoP3bwVZgfoyg/j8L8GIV5MYoKEl/zDl9PPi4/Rp4ntkhEkmQk5QG3A28A6oClkhaY2eqkw64HdpnZNEnXAt8Erunu3EX53gDqXG8V5MUYNbQw5cFB29uN/c2t7GtsZX9jKwda2mhobuVAc1u43HbYcmO4/9By8NpZ30zdruRjW2lsObLaLxX5MYUJKkZBXuIl8pOXYzq4Lz9PSduTjw/WC/Nj5MeC9xfmiVgseH9MIi+x3GFb4tXZtryYyOtsW4ftMQUxHLYt6Zz54bZs+Wc7qpLMHGCNma0DkHQ3cAWQnGSuAL4YLt8PfF+SrJsJcAo9yTjX72IxBdVlxel/3qe93Whua6e5rZ2mlsTXtmBbaztNrYmvbQfXD20Lvnbc39zaTkt7O61tRmt7Oy1tRktbsN7S1k5Dcyut7Za0PemY9uBr4vjW9uyck0visMQTC5djCpoVlLScSVElmQnAxqT1OuDMox1jZq2S9gDlwI6uTlzkXTmdyymxmCiOBVVoFEcdzZHMgmTUbkZbe5B02hNfw22JV2fb2jquJ7a1Hb6v3YzWo2zr6bUB2s3CVxB7ezs8lcH7kxMN/5JuAG4AmFA5OeJonHODiSQK87OjaipV38rguaP6t38TUJm0PjHc1ukxkvKBEQQdAI5gZneaWbWZVR8zxmfEdM65bBFVklkKTJc0RVIhcC2woMMxC4D3hMtXAn/vrj3GOedcdomkuixsY7kJeJigC/PPzOx5SV8GlpnZAuCnwK8krQF2EiQi55xzA0hkbTJmthBY2GHb55OWG4Gr+jsu55xz6eNdsZxzzmWMJxnnnHMZ40nGOedcxniScc45lzGeZJxzzmWMcu3RE0nbgdqo4+hnFXQz3M4g4/fjSH5PDuf343DHm9mRM+ilQU4MK5PMzEZHHUN/k7TMzKqjjiNb+P04kt+Tw/n9OJykZZk6t1eXOeecyxhPMs455zLGk0xuuDPqALKM348j+T05nN+Pw2XsfuRcw79zzrns4SUZ55xzGeNJJotI+pmkbZJWJW0rk/SIpFfCr6PC7ZJ0q6Q1kp6TdEbSe94THv+KpPckbZ8laWX4nluVLZOAd0JSpaRHJa2W9Lykj4bbB+X9AJBULGmJpGfDe/KlcPsUSU+Fn+OecPoMJBWF62vC/VVJ5/pMuP0lSRclbb843LZG0qf7+zOmQlKepOWS/hCuD/b7URP+XK9I9BqL9PfGzPyVJS/gNcAZwKqkbd8CPh0ufxr4Zrh8KfAnQMBZwFPh9jJgXfh1VLg8Kty3JDxW4Xsvifozd3EvxgFnhMulwMvAjMF6P8J4BQwLlwsIZs09C7gXuDbcfgfwoXD5X4A7wuVrgXvC5RnAs0ARMAVYSzDlRl64fCxQGB4zI+rP3YP78nHgLuAP4fpgvx81QEWHbZH93nhJJouY2T8I5s5JdgXwi3D5F8Bbkrb/0gKLgZGSxgEXAY+Y2U4z2wU8Alwc7htuZost+En5ZdK5so6ZbTGzZ8LlfcALwAQG6f0ACD/b/nC1IHwZcAFwf7i94z1J3Kv7gdeH/3VeAdxtZk1mth5YA8wJX2vMbJ2ZNQN3h8dmLUkTgTcBPwnXxSC+H12I7PfGk0z2G2tmW8LlV4Gx4fIEYGPScXXhtq6213WyPeuF1RozCf5zH9T3I6waWgFsI/jFXwvsNrPW8JDkz3Hws4f79wDl9P5eZbP/AT4FtIfr5Qzu+wHBPx5/kfS0pBvCbZH93uTcE/+5zMxM0qDqDihpGPAA8DEz25tc/TsY74eZtQGnSxoJ/BY4IeKQIiPpMmCbmT0t6fyo48ki55rZJkljgEckvZi8s79/b7wkk/22hkVUwq/bwu2bgMqk4yaG27raPrGT7VlLUgFBgvmNmT0Ybh609yOZme0GHgXOJqjiSPzDmPw5Dn72cP8IIE7v71W2mgtcLqmGoCrrAuB7DN77AYCZbQq/biP4R2QOEf7eeJLJfguARM+O9wC/S9r+7rB3yFnAnrA4/DDwRkmjwh4kbwQeDvftlXRWWA/97qRzZZ0wxp8CL5jZd5J2Dcr7ASBpdFiCQVIJ8AaCtqpHgSvDwzrek8S9uhL4e1iPvgC4NuxtNQWYTtCYuxSYHvbOKiRoHF+Q+U+WGjP7jJlNNLMqglj/bmbvYJDeDwBJQyWVJpYJft5XEeXvTZS9IPx1RK+Q+cAWoIWgrvN6gjrjvwGvAH8FysJjBdxOUCe/EqhOOs/7CRov1wDvS9peHf7ArQW+T/gwbja+gHMJ6pafA1aEr0sH6/0I4z0VWB7ek1XA58PtxxL8UVwD3AcUhduLw/U14f5jk8712fBzv0RS76DwHr8c7vts1J+5F/fmfA71Lhu09yP87M+Gr+cTMUf5e+NP/DvnnMsYry5zzjmXMZ5knHPOZYwnGeeccxnjScY551zGeJJxzjmXMZ5knIuYpHxJ/x4+i+FcTvEk41z0ZgDrgdOiDsS5dPMk41z0XgSmEjxw6lxO8YcxnXPOZYyXZJxLkaRTJNVK+lAPjq2SdCAcpj95+1skmaQTkraVhLMaNkuqyETszvUXTzLOpcjMVhIMmvjuHr5lrZmd3mHbPOCJ8GvivAfC4zanJVDnIuRJxrm+2QaclMobw7lyziUYCPXadAblXLbwJONc33wDKJI0OYX3XgH82cxeBuKSZqU3NOei50nGuRRJugQYCvyR1Eoz8wgm2yL8Oq+LY50bkHz6ZedSIKkY+CZwOfA+4GRgYS/eX0Ywk+Mp4VS4eYBJutm8y6fLIV6ScS41nwN+aWY1BJM9nZzYIelvkiZ08/4rgV+Z2WQzqzKzSoIHMs/LVMDORcGTjHO9JOl4gqmP/yfcdDDJSIoB04Cd3ZxmHsH868kewKvMXI7x6jLnesnMXgLO7LB+Rrg6A3jAzA50c47XdbLt1nTG6Vw28JKMc2lkZqvM7OOd7GoDRnR8GLMziYcxgQKgPd0xOteffFgZ55xzGeMlGeeccxnjScY551zGeJJxzjmXMZ5knHPOZYwnGeeccxnjScY551zGeJJxzjmXMZ5knHPOZYwnGeeccxnz/wFK/8LF9QZA9gAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + } + } + ], + "source": [ + "lam = np.linspace(nu.max().to(u.AA, equivalencies=u.spectral()),\n", + " nu.min().to(u.AA, equivalencies=u.spectral()), 1000)\n", + "bb_lam = BlackBody(bb.temperature, \n", + " scale=1.0 * u.erg / (u.cm ** 2 * u.AA * u.s * u.sr))\n", + "bb5000K_lam = bb_lam(lam)\n", + "plt.plot(lam, bb5000K_lam)\n", + "plt.xlim([1.0e3, 5.0e4])\n", + "plt.xlabel(r'$\\lambda$, [{0:latex_inline}]'.format(lam.unit))\n", + "plt.ylabel(r'$I_{\\lambda}$, ' + '[{0:latex_inline}]'.format(bb5000K_lam.unit))\n", + "plt.title('Planck function in wavelength')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "id": "A4VLHFtHXNbN", + "outputId": "01a83a99-4944-4d30-968a-0375a0877477", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 38 + } + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/latex": "$1.146776 \\times 10^{10} \\; \\mathrm{\\frac{erg}{s\\,sr\\,cm^{2}}}$", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "execution_count": 5 + } + ], + "source": [ + "np.trapz(x=lam, y=bb5000K_lam).to('erg s-1 cm-2 sr-1')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "-BC40k7BXNbN" + }, + "source": [ + "Observe que isso está dentro de alguns por cento da resposta que obtivemos no espaço de frequência, apesar de nossa amostragem ruim em pequenos comprimentos de onda!\n", + "\n", + "Muitas funções `astropy` usam unidades e quantidades diretamente. À medida que você ganha confiança ao trabalhar com eles, considere incorporá-los ao seu fluxo de trabalho regular. Leia mais [aqui](http://docs.astropy.org/en/stable/units/) sobre como usar unidades." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "rlwI6oPGXNbO" + }, + "source": [ + "###Como simular observações reais\n", + "\n", + "Desde outono de 2017, o `astropy` não suporta explicitamente a construção de observações sintéticas de modelos como curvas de corpo negro. A [biblioteca synphot](https://synphot.readthedocs.io/en/latest/) permite isso. Você pode usar `synphot` para realizar tarefas como transformar espectros em magnitudes visuais por convolução com uma curva de filtro." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "88l25iyEXNbO" + }, + "source": [ + "## A função de massa inicial estelar (IMF)\n", + "\n", + "A função de massa inicial estelar nos diz quantas estrelas de cada massa são formadas. Em particular, as estrelas de baixa massa são muito mais abundantes do que as estrelas de alta massa. Vamos explorar mais a funcionalidade do `astropy` usando esta ideia.\n", + "\n", + "As pessoas geralmente pensam no FMI como uma função de densidade de probabilidade de lei de potência. Em outras palavras, se você contar as estrelas que nasceram recentemente de uma nuvem de gás, sua distribuição de massas seguirá o FMI. Vamos escrever uma pequena classe para nos ajudar a acompanhar isso:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "id": "NToG17WYXNbO" + }, + "outputs": [], + "source": [ + "class PowerLawPDF(object):\n", + " def __init__(self, gamma, B=1.):\n", + " self.gamma = gamma\n", + " self.B = B\n", + " def __call__(self, x):\n", + " return x**self.gamma / self.B" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "gqNFgJqNXNbP" + }, + "source": [ + "### O método `__call__` \n", + "\n", + "Ao definir o método `__call__`, estamos dizendo ao interpretador Python que uma instância da classe pode ser chamada como uma função. Quando chamada, uma instância desta classe recebe um único argumento, `x`, mas usa outros atributos da instância, como `gamma` e `B`.\n", + "\n", + "### Mais sobre as aulas\n", + "\n", + "As classes são estruturas de dados mais avançadas, que podem ajudá-lo a acompanhar a funcionalidade em seu código que funciona em conjunto. Você pode aprender mais sobre aulas neste [tutorial](https://www.codecademy.com/ja/courses/learn-python/lessons/introduction-to-classes/exercises/why-use-classes)." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "SPGk79FOXNbP" + }, + "source": [ + "## Integração usando quadratura gaussiana\n", + "\n", + "Nesta seção, exploraremos um método de integração numérica que não requer já ter sua grade de amostragem configurada. `scipy.integrate.quad` com referência [aqui](https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.integrate.quad.html) recebe uma função e um valor menor e limite superior, e nossa classe `PowerLawPDF` cuida disso muito bem.\n", + "\n", + "Agora podemos usar nossa nova classe para normalizar nosso IMF dados os limites de massa. Isso equivale a normalizar uma função de densidade de probabilidade. Usaremos a quadratura gaussiana (`quad`) para encontrar a integral. `quad` retorna o valor numérico da integral e sua incerteza. Nós só nos importamos com o valor numérico, então vamos empacotar a incerteza em `_` (uma variável de espaço reservado). Nós imediatamente jogamos a integral em nosso objeto IMF e a usamos para normalizar!\n", + "\n", + "Para ler mais sobre *empacotamento e desempacotamento generalizado* em Python, veja a proposta original, [PEP 448](https://www.python.org/dev/peps/pep-0448/), que foi aceita em 2015." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "id": "8_PJmxbYXNbP", + "outputId": "35681e7d-0473-4b2d-dc24-5deba60a4037", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 287 + } + }, + "outputs": [ + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEOCAYAAACXX1DeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUddrG8e+ThN6kBKUKSBMJKkSkJha6BUUs2F0VQbHAFvXVXcuuq65rQBAUC6KuiqygIhaKhUiTJtJ7EbAAotJrnvePDO/mZZNxQjKZmcz9ua65kvObmXNuf4Z55pwz8xxzd0RERPKSEOkAIiIS3VQoREQkKBUKEREJSoVCRESCUqEQEZGgVChERCSopEgHCIdq1ap5vXr1Ih1DRCSmzJ8/f7u7Jx87XiwLRb169Zg3b16kY4iIxBQz25jbuA49iYhIUCoUIiISlAqFiIgEpUIhIiJBqVCIiEhQKhQiIhKUCkUO23YdYOuu/ZGOISISVaK+UJjZJWb2opm9bWZdwrmthz9YSueMTP49bxO6ToeISLaIFAozG2VmW81syTHj3cxspZmtMbP7ANz9PXe/FegHXBnOXIM6N6bxieX54zuLuH7UHDbt2BvOzYmIxIRI7VGMBrrlHDCzRGA40B1oBvQxs2Y5HvJg4P6wOSW5PG/3bcujPU9jwcaf6Tokk9Ez1pOVpb0LEYlfESkU7p4J7DhmuDWwxt3XuftBYAzQ07I9CXzs7gvyWqeZ9TWzeWY2b9u2bcedLSHBuL5tPSYNTCO1XhUe/mAZl4+cxZqtu457nSIisSyazlHUAjblWN4cGLsT6AT0NrN+eT3Z3V9w91R3T01O/q+eVvlWu3JZXr3pLJ6+/HTWbN1Nj2emM/zzNRw6klXgdYuIxJKobwro7kOBoZHYtplxWavapDVO5uEJS3lq0komLvqep3q3oHmtSpGIJCJS5KJpj2ILUCfHcu3AWMQlVyjF8Gta8vy1rdi++wA9h8/giY9XsP/QkUhHExEJu2gqFHOBRmZW38xKAlcBEyKc6f/p1vwkpg5Mp3fL2jw/bS09nvmSOeuPPdUiIlK8ROrjsW8Bs4AmZrbZzG5298PAAGASsBwY6+5LI5EvmEplS/Bk7xb86+azOXgkiytGzuLP7y1h94HDkY4mIhIWVhy/WJaamupFceGivQcP889Jq3hl5npqVCzNY71SOLdJ9bBvV0QkHMxsvrunHjseTYeeYk7Zkkn85aJmjOvfjnKlkrjplbkMenshP+85GOloIiKFRoWiELSsW5mJd3XgrvMaMuGb7+iUMY2Ji75TGxARKRZUKApJqaREBnVpwgd3dqDmCWUY8ObX3Pb6fH7cqSaDIhLbVCgK2ak1KvLu7e24r3tTpq3aRqeMabw991vtXYhIzFKhCIOkxAT6pZ/CJ/ekcWqNitw7bjHXvawmgyISm1Qowqh+tXKMubUNf7ukOQs3/UKXwZm8PH09R9RkUERiiApFmCUkGNe2OZnJA9No06AKf524jN7Pz2T1j2oyKCKxQYWiiNQ8oQyjbjyLIVeewYbte7hg6HSGfrqag4fVZFBEopsKRREyMy45sxZTBqXTtflJZExZxcXPTmfR5l8iHU1EJE8qFBFQrXwphvU5kxevT+XnvQe5ZPgMHv9oOfsOqsmgiEQfFYoI6tzsRKYMSufKs+owMnMd3Z/JZPa6nyIdS0Tk/1GhiLCKpUvweK8WvHnL2WQ5XPXCbB54dzG79h+KdDQREUCFImq0a1iNSfekcWvH+rw151u6DM7ksxU/RjqWiIgKRTQpUzKRBy5oxvjb21OxdAl+N3oed4/5mp92H4h0NBGJYyoUUeiMOifwwZ0duKdTIz5a/D2dB2cy4Rs1GRSRyFChiFIlkxK4p1NjJt7ZkTpVynLXW19z62vz+eFXNRkUkaKlQhHlmpxUgfH92/FAj1OZvmYbnTOm8dYcNRkUkaKjQhEDEhOMW9MaMOmeNJrXqsT94xdz9YtfsfGnPZGOJiJxQIUihpxctRxv3no2j/dKYcmWX+k6JJOXvlynJoMiElYqFDHGzOjTui5TBqXToWE1/vbhcno9N5OVP6jJoIiEhwpFjDqpUmlevD6VoX3OZNOOvVw47EuGTF2lJoMiUuhUKGKYmXHx6TWZOiidC1JqMGTqai4aNp2Fm9RkUEQKjwpFMVClXEmGXHUmo25MZef+Q/QaMYO/TVymJoMiUihUKIqR85qeyOSBafRpXZeXpq+n65BMZq7ZHulYIhLjVCiKmQqlS/DYpSmM6duGBIOrX/qK+8cvYqeaDIrIcVKhKKbaNKjKx3encVt6A96eu4nOGdOYskxNBkUk/1QoirEyJRO5v/upvHdHeyqXLcmtr81jwJsL2K4mgyKSDyoUcaBF7ROYMKADgzo3ZvLSH+mcMY33vt6iNiAiEhIVijhRMimBu85vxId3deDkquW45+2F3PzqPL77ZV+ko4lIlFOhiDONTqzAuP7t+POFzZi19ie6DM7kX7M3kqU2ICKSBxWKOJSYYNzcoT6T7knj9DqVePC9JfR5cTbrt6vJoIj8NxWKOFa3aln+dfPZPHlZCsu+30m3IZmMnLaWw0fUBkRE/iPqC4WZlTOzV83sRTO7JtJ5ihsz48qz6jJ1UDppjZN5/OMV9HpuJsu/3xnpaCISJSJSKMxslJltNbMlx4x3M7OVZrbGzO4LDPcC3nH3W4GLizxsnDixYmleuK4Vz159Jt/9so+Lhk0nY/JKDhxWGxCReBepPYrRQLecA2aWCAwHugPNgD5m1gyoDWwKPEyvWmFkZlzYoiZTBqZz8ek1GfrZGi4cOp0F3/4c6WgiEkERKRTungnsOGa4NbDG3de5+0FgDNAT2Ex2sYAgec2sr5nNM7N527ZtC0fsuFG5XEkyrjyDV246iz0HDnPZczN59INl7D14ONLRRCQCoukcRS3+s+cA2QWiFjAeuMzMngM+yOvJ7v6Cu6e6e2pycnJ4k8aJc5tUZ/KgdK5rczKjZqyny+BMpq9Wk0GReBNNhSJX7r7H3W9y9/7u/kak88Sb8qWSeLRnc8be1pYSiQlc+/JX/Omdb/h1r5oMisSLaCoUW4A6OZZrB8YkCrSuX4WP7+5I/3NOYdyCLXQaPI1PlvwQ6VgiUgSiqVDMBRqZWX0zKwlcBUyIcCbJoXSJRO7t1pT372hPcvlS9PvXfG5/Yz5bd+2PdDQRCaNIfTz2LWAW0MTMNpvZze5+GBgATAKWA2PdfWkk8klwzWtV4v0B7flj1yZMXbaVzhmZjJu/WU0GRYopK47/uFNTU33evHmRjhEX1mzdzb3jFjF/48+kN07m771SqHVCmUjHEpHjYGbz3T312PFoOvQkMahh9fKMva0tD1/UjLkbdtAlYxqvzdqgJoMixYgKhRRYYoJxY/vsJoMtT67MX95fypUvzGLttt2RjiYiheA3C4WZpRRFEIl9daqU5bXfteap3i1Y+cMuuj/zJSO+WMMhNRkUiWmh7FGMMLM5Zna7mVUKeyKJaWbG5al1mDoonfOaVOcfn6zkkuEzWLLl10hHE5Hj9JuFwt07AteQ/R2H+Wb2ppl1DnsyiWnVK5bm+eta8dw1Lflx5wF6Dp/BU5NWsP+Q2nWJxJqQzlG4+2rgQeBeIB0YamYrzKxXOMNJ7OueUoOpg9LodWYthn++lh5Dv2TehmPbfIlINAvlHEULMxtM9ncbzgMucvdTA78PDnM+KQZOKFuSpy4/ndd+15oDh7K4fOQsHp6wlD0H1GRQJBaEskcxDFgAnO7ud7j7AgB3/47svQyRkKQ1TmbywDRuaFuPV2dtoMvgTKatUqdfkWgXSqF4191fd/d9RwfM7G4Ad389bMmkWCpXKomHLz6Nf9/WltIlErhh1Bx+P/Ybftl7MNLRRCQPoRSK63MZu7GQc0icSa1XhQ/v6siAcxvy/sItdMrI5KPF30c6lojkItiFgPqY2QdAfTObkOP2Of990SGRfCtdIpE/dG3C+wPac1KlUtz+xgL6vT6frTvVZFAkmiQFuW8m8D1QDXg6x/guYFE4Q0l8Oa1mJd67vT0vfrmewVNXMTNjOw9e2IzLW9XGzCIdTyTuqSmgRJV127KbDM7d8DMdG1Xj75emUKdK2UjHEokL+W4KaGbTAz93mdnOHLddZrYznGElfjVILs/bfdvyaM/TWLDxZ7oOyeSVGevVZFAkgvIsFO7eIfCzgrtXzHGr4O4Viy6ixJuEBOP6tvWYNDCNs+pV4ZEPlnH5yFms2bor0tFE4lIoX7g7xcxKBX4/x8zuMrMTwh9N4l3tymUZfdNZZFxxOmu37abHM9N59rPVajIoUsRC+XjsOOCImTUEXiC759ObYU0lEmBm9GpZmykD0+l82on8c/IqLn5WTQZFilIohSIrcJnSS4Fh7v5HoEZ4Y4n8f8kVSjH86paMvK4VP+3ObjL4xMdqMihSFEIpFIfMrA9wAzAxMFYifJFE8tb1tJOYMjCd3i1r8/y0tfR45kvmrNfXekTCKZRCcRPQFnjM3debWX1ArTskYiqVLcGTvVvwxi1ncygriytGzuLP7y1h1/5DkY4mUizpexQS0/YePMw/J63ilZnrqVGxNI/1SuHcJtUjHUskJuX7exQ5ntjezKaY2SozW2dm681sXXhiiuRP2ZJJ/OWiZozr345ypZK46ZW5DHp7IT/vUZNBkcLym3sUZrYCGAjMB/7vzKG7/xTeaMdPexTx6cDhIwz/fC0jPl9DpTIleKTnaVyQUkNtQERCdNx7FMCv7v6xu29195+O3sKQUaRASiUlMqhzYz64swO1KpdhwJtf0/f1+fyoJoMiBRJKofjczJ4ys7Zm1vLoLezJRI7TqTUqMr5/O/6nR1MyV22jU8Y03p77LcXxfJxIUQjl0NPnuQy7u58XnkgFp0NPctSG7Xu4d9wivlq/g/YNq/L4pS2oW1VNBkVyk9ehJ33qSYq9rCznrbnf8vhHKziS5fyhaxNubFePxASduxDJqSCfejrRzF42s48Dy83M7OZwhBQJh4QE45qzT2bywDTaNKjCXycuo/fzM1n9o5oMioQilHMUo4FJQM3A8irgnnAFEgmXmieUYdSNZ/HMVWewYfseLhg6naGfrubgYTUZFAkmlEJRzd3HAlkAgb5ParAjMcnM6HlGLaYOSqdr85PImLKKi5+dzjebfol0NJGoFUqh2GNmVQEHMLM2gFp3SkyrWr4Uw/qcyYvXp/Lz3oNcOmIGj3+0nH0H9R5I5FjBrpl91CBgAnCKmc0AkoHeYU0lUkQ6NzuRsxtU4fGPljMycx2Tlv7A471a0PaUqpGOJhI1QvrUk5klAU0AA1a6e5F2XzOzS4ALgIrAy+4+Odjj9aknOR4z12znvvGL+XbHXq4+uy73dW9KxdJqlCzxI98fjzWzXsFW6O7jQ9zwKOBCYKu7N88x3g14BkgEXnL3J0JYV2Xgn+4e9FNXKhRyvPYdPELGlJW8PH091SuU5rFLm3P+qSdGOpZIkTieQvFK4NfqQDvgs8DyucBMd78wxA2nAbuB144WCjNLJPvTU52BzcBcoA/ZRePxY1bxO3ffGnje08Ab7r4g2DZVKKSgFm76hXvfWcTKH3fR84ya/OXCZlQtXyrSsUTC6ri/cGdmk4Eb3P37wHINYLS7d83HxusBE3MUirbAw0fXYWb3A7j7sUXi6PMNeAKY4u5T83hMX6AvQN26dVtt3Lgx1HgiuTp4OIvnvljLs5+vpkLpEjx0UTMuPr2mmgxKsVWQpoB1jhaJgB+BugXMUwvYlGN5c2AsL3cCnYDeZtYvtwe4+wvunuruqcnJyQWMJwIlkxK4u1MjJt7ZkTpVynL3mIXc+to8vv91X6SjiRSpUArFp2Y2ycxuNLMbgQ+BXN/Vh4u7D3X3Vu7ez92fL8ptizQ5qQLj+7fjwQtOZfqa7XTJyOTNr74lK6v4tb8Ryc1vFgp3HwA8D5weuL3g7ncWcLtbgDo5lmsHxkSiUmKCcUvHBky6J43mtSrxP+8u5uqXZrPxpz2RjiYSdkXSFDCXcxRJZJ/MPp/sAjEXuNrdlxbG9nQyW8LJ3RkzdxN//3A5h7Ky+EOXJtzUvr6aDErMK8g5ioJu+C1gFtDEzDab2c2BNiADyO4htRwYW1hFQiTczIw+resyZVA6HRpW428fLqfXczNZ+YOaDErxpDbjIgXg7nyw6HsembCUnfsPcfs5Dbnj3IaUTAr7ezCRQleQNuMXmZn+6kVyYWZcfHpNpgxK54KUGjzz6WouHPYlC9VkUIqRUArAlcBqM/uHmTUNdyCRWFSlXEmGXHUmo25MZdf+w/QaMYO/TVymJoNSLITyqadrgTOBtcBoM5tlZn3NrELY04nEmPOansjkgWn0aV2Xl6avp+uQTGau2R7pWCIFEtIhJXffCbwDjAFqAJcCC8ysoB+TFSl2KpQuwWOXpjCmbxsSDK5+6SvuG7eIX/cVaS9NkUITyjmKnmb2LvAFUAJo7e7dyf5Oxe/DG08kdrVpUJVP7knjtvQGjJ23ic4Z05i89IdIxxLJt1D2KHoBg909xd2fOtqgz933Arp2tkgQpUskcn/3U3nvjvZUKVeSvq/PZ8CbC9i++0Cko4mELJRC8YO7Z+YcMLMnAdz907CkEilmWtQ+gQkDOvD7zo2ZvPRHOmVM492vN1McP54uxU8ohaJzLmPdCzuISHFXMimBO89vxId3daB+tXIMfPsbbho9ly2/qMmgRLc8C4WZ9TezxUBTM1uU47YeWFR0EUWKl0YnVuCdfu34y4XN+GrdDrpkTOP12RvVZFCiVrALF1UCKpN9IaH7cty1y913FEG246ZvZkus2LRjL/ePX8z0NdtpXb8KT17WgvrVykU6lsSp4/lmtrv7BuAOYFeOG2ZWJRwhReJNnSplef3m1vzjshYs/34n3YZk8vy0tRw+khXpaCL/J1iheDPwcz4wL/Bzfo5lESkEZsYVZ9Xh00HppDdO5omPV3DpiJks+25npKOJAGoKKBJV3J2Pl/zAX95fwi97D9H/nFMYcF5DSiUlRjqaxIG8Dj0lBXlCy2ArdPcFhRFMRP7DzOiRUoO2Dary1w+XMeyzNXy85AeevKwFrU6uHOl4EqeCncz+PMjz3N3PC0+kgtMehRQXX6zcygPvLuG7X/dxY7t6/KFLE8qVyvP9nUiB5LVHoUNPIlFu94HD/OOTFbw2ayO1K5fh8V4pdGyUHOlYUgzlu1CY2Xnu/pmZ9crtfncfX8gZC40KhRRHc9bv4L5xi1i3fQ9XpNbmgR7NqFS2RKRjSTGS73MUQDrwGXBRLvc5ELWFQqQ4al2/Ch/d3ZGhn65mZOY6Pl+5jb/2bE635idFOpoUczr0JBKDlmz5lT+9s4hl3++kR8pJPHJxc5IrlIp0LIlxBbkUalUzG2pmC8xsvpk9Y2ZVwxNTRELRvFYl3h/Qnj92bcLU5VvplDGNcfPVZFDCI5SmgGOAbcBlQO/A72+HM5SI/LYSiQnccW5DPrqrIw2rl+f3//6GG16Zy+af90Y6mhQzv3noycyWuHvzY8YWu3tKWJMVgA49SbzJynJen72RJz9ZAcC93ZpyXZuTSUiwCCeTWHLch56AyWZ2lZklBG5XAJMKP6KIHK+EBOOGdvWYdE8arU6uzEMTlnLFyFms3bY70tGkGAj28dhdZH+6yYBywNEuZQnAbnevWCQJj4P2KCSeuTvjFmzhrxOXse/QEe4+vxF90xpQIjGU94USz/K9R+HuFdy9YuBngrsnBW4J0VwkROKdmdG7VW2mDErj/KbVeWrSSi4ZPoMlW36NdDSJUSG9xTCzymbW2szSjt7CHUxECqZ6hdI8d20rnrumJT/uPEDP4TP4xycr2H/oSKSjSYwJ5eOxtwCZZJ+XeCTw8+HwxhKRwtI9pQafDkrn0jNrMeKLtfQY+iVzN0T1tcckyoSyR3E3cBaw0d3PBc4EfglrKhEpVJXKluCfl5/Oa79rzYFDWVwxchYPvb+E3QcORzqaxIBQCsV+d98PYGal3H0F0CS8sUQkHNIaJzN5YBo3tK3Ha7M30nVwJtNWbYt0LIlyoRSKzWZ2AvAeMMXM3gc2hjeWiIRLuVJJPHzxabzTry2lSyRww6g5/H7sN/yy92Cko0mUylevJzNLByoBn7h71P5V6eOxIqHZf+gIz362huenreWEsiV4tGdzeqTUiHQsiZCCfOEOM2tpZncBLYDN0VwkRCR0pUsk8oeuTXh/QHtOqlSa299YQL/X57N15/5IR5MoEsqnnv4CvApUBaoBr5jZg+EOdkyGcmY2z8wuLMrtisSL02pW4r3b23Nvt6Z8tjK7yeDYeZvUZFCA0PYorgHOcveH3P0hoA1wXSgrN7NRZrbVzJYcM97NzFaa2Rozuy+EVd0LjA1lmyJyfJISE+h/zil8cndHmp5UkT+9s4jrR81h0w41GYx3oRSK74DSOZZLAVtCXP9ooFvOATNLBIYD3YFmQB8za2ZmKWY28ZhbdTPrDCwDtoa4TREpgAbJ5RnTtw1/7XkaCzb+TNchmbwyYz1HsrR3Ea/yvMKdmQ0ju9fTr8BSM5sSWO4MzAll5e6eaWb1jhluDaxx93WB7YwBerr748B/HVoys3PI7jXVDNhnZh+5e9axjxORwpOQYFzXth7nnXoi/zN+MY98sIyJi77nyctSaFi9QqTjSRELdinUox8bmg+8m2P8iwJusxawKcfyZuDsvB7s7g8AmNmNwPa8ioSZ9QX6AtStW7eAEUUEoNYJZRh901m8+/UWHp24jB7PTOeu8xtyW/opajIYR/IsFO7+6tHfzawk0DiwuNLdD4U7WC55Rv/G/S8AL0D2x2OLIpNIPDAzerWsTcdGyTz8wVL+OXkVHy7+gad6t6B5rUqRjidFIJRPPZ0DrCb7vMIIYFUBmwJuAerkWK5N6Oc8RCRCkiuUYvjVLRl5XSt+2p3dZPCJj9VkMB6Esu/4NNDF3dPdPQ3oCgwuwDbnAo3MrH5gT+UqYEIB1iciRajraScxZVA6vVvW5vlpa+nxzJfMWa8mg8VZKIWihLuvPLrg7quAEqGs3MzeAmYBTcxss5nd7O6HgQFkd6FdDox196X5jy4ikVKpTAme7N2CN245m0NZ2U0G//zeEnbtL/Kj0lIEQrlm9ivAEeBfgaFrgER3/12Ysx03tfAQKTp7Dx7m6cmrGDVjPTUqluaxXimc26R6pGPJcShIC49+ZH+P4a7AbRnQv3DjiUisKlsyiT9f2Ixx/dtRrlQSN70yl4FvL2THHnX6KS6C7lEEvhy31N2bFl2kgtMehUhkHDh8hOGfr2XE52uoVKYEj/Q8jQtSamBmkY4mITiuPQp3PwKsNDN9MUFEflOppEQGdW7MB3d2oFblMgx482v6vj6fH9VkMKaFcuipMtnfzP7UzCYcvYU7mIjErlNrVGR8/3b8T4+mZK7aRqeMaYyZ862aDMaoUE5mp+c27u7TwpKoEOjQk0j02LB9D/eNX8TsdTtod0pVHu+VwslVy0U6luQir0NPeRYKMytN9onshsBi4OXAR1ujngqFSHTJynLGzN3E3z9azuGsLP7QpQk3ta9PYoLOXUST4zlH8SqQSnaR6E72F+9ERPItIcG4+uy6TBmURrtTqvG3D5dz2XMzWfXjrkhHkxAEKxTN3P1adx8J9AY6FlEmESmmalQqw8s3pPLMVWew8ac9XDD0S4Z+upqDh9UQOpoFKxT/9xXLWDnkJCLRz8zoeUYtpg5Kp1vzGmRMWcXFz07nm02/RDqa5CFYoTjdzHYGbruAFkd/N7OdRRVQRIqnquVLMazPmbx0fSq/7D3EpSNm8PePlrPvoJoMRptgbcYTizKIiMSnTs1OpHWDKjz+0QpeyFzHpKU/8ESvFrQ9pWqko0mArjwiIhFXsXQJHu+Vwpu3Zl/DrM+Ls7l//GJ2qslgVFChEJGo0e6Uanxydxq3dqzP23O/pUtGJp8u/zHSseKeCoWIRJUyJRN54IJmjL+9PZXKlODmV+dx11tf89PuA5GOFrdUKEQkKp1R5wQ+uLMDAzs15uMl39N5cCbvL9yiNiARoEIhIlGrZFICd3dqxMQ7O1KnSlnuHrOQW16dx/e/7ot0tLiiQiEiUa/JSRUY378dD15wKjPWbqdzRiZvfLWRrCztXRQFFQoRiQmJCcYtHRsw+Z50WtSuxAPvLuHql2azYfueSEcr9lQoRCSm1K1aljduOZsneqWwdMtOug7J5IXMtRw+ojYg4aJCISIxx8y4qnVdpgxKp2OjZP7+0Qoue24mK35Q04hwUKEQkZh1UqXSvHh9K4b1OZPNP+/jwqHTyZiyigOH1QakMKlQiEhMMzMuOr0mUwalc9HpNRn66WouGjadr7/9OdLRig0VChEpFqqUK8ngK89g1I2p7Np/mF7PzeSvE5ex96CaXxeUCoWIFCvnNT2RyQPTuLp1XV6evp5uQ75k5prtkY4V01QoRKTYqVC6BI9dmsLbfduQmGBc/dJX3DduEb/uU5PB46FCISLF1tkNqvLx3R25Lb0BY+dtonPGNCYv/SHSsWKOCoWIFGulSyRyf/dTee+O9lQpV5K+r89nwJsL2K4mgyFToRCRuNCidnaTwT90aczkpT/SKWMa7369WU0GQ6BCISJxo0RiAgPOa8RHd3egQbVyDHz7G24aPZctv6jJYDAqFCISdxpWr8C/+7XjoYua8dW6HXTJmMbrs9VkMC8qFCISlxITjJva12fywDRanlyZP7+3hKtemM26bbsjHS3qqFCISFyrU6Usr/2uNf/o3YIVP+yk2zNf8twXajKYU9QXCjNLMLPHzGyYmd0Q6TwiUvyYGVek1mHqoHTObZLMk5+s4JIRM1j2nZoMQpgLhZmNMrOtZrbkmPFuZrbSzNaY2X2/sZqeQG3gELA5XFlFRKpXLM3z17ZixDUt+eHX/Vz87HSenrwy7psMhnuPYjTQLeeAmSUCw4HuQDOgj5k1M7MUM5t4zK060ASY6e6DgP5hzisicc7M6JFSgykD07n4jJoM+2wNFwydzvyN8dtkMKyFwt0zgR3HDLcG1rj7Onc/CIwBerr7Yne/8JjbVrRtRLIAAAsuSURBVLL3Io7+H4rvsi4iRaZyuZJkXHEGo286i30Hj9D7+Zk88sFS9hyIvyaDkThHUQvYlGN5c2AsL+OBrmY2DMjM60Fm1tfM5pnZvG3bthVOUhGJe+c0qc6kgWlc3+ZkXpmxga5DMvlydXy9xkT9yWx33+vuN7v7ne4+PMjjXnD3VHdPTU5OLsqIIlLMlS+VxCM9m/Pvfm0pmZTAdS/P4Y///oZf98ZHk8FIFIotQJ0cy7UDYyIiUe2selX46K6O3H7OKYz/egudBk/jkyXFv8lgJArFXKCRmdU3s5LAVcCECOQQEcm30iUS+VO3prx/R3uSy5ei37/mc/sb89m6a3+ko4VNuD8e+xYwC2hiZpvN7GZ3PwwMACYBy4Gx7r40nDlERApb81qVeH9Ae/7UrQlTl2+lc0Ym78wvnk0GrTj+R6Wmpvq8efMiHUNE4sSarbu5b9wi5m38mbTGyfz90ubUrlw20rHyzczmu3vqseNRfzJbRCTaNaxenrG3teXRnqcxb8MOugzO5NWZG4pNk0EVChGRQpCQYFzfth6TB6aRWq8KD01YyhUjZ7G2GDQZVKEQESlEtSuX5dWbzuLpy09n9dbddH/mS4Z/voZDMdxkUIVCRKSQmRmXtarN1EHpdDq1Ok9NWsklw2ewZMuvkY52XFQoRETCJLlCKUZc04rnr23J1l0H6Dl8Bv/4ZAX7D8VWNyIVChGRMOvWvAZTB6bT68xajPhiLT2GfsncDce2wYteKhQiIkWgUtkSPHX56bx+c2sOHs7iipGzeOj9JeyOgSaDKhQiIkWoY6NkJt2Txo3t6vHa7I10HZzJtFXR3WRQhUJEpIiVK5XEQxedxjv92lGmZCI3jJrDoLEL+WXvwUhHy5UKhYhIhLQ6uTIf3tWBO89ryISF39EpYxofLf4+0rH+iwqFiEgElUpK5PddmjBhQAdqVCrD7W8s4LbX57F1Z/Q0GVShEBGJAs1qVuTd29txf/emfLFyG50ypjF27qaoaDKoQiEiEiWSEhO4Lf0UPrknjaY1KvKncYu47uU5bNqxN6K5VChERKJM/WrlGHNrG/52SXMWbvqFLoMzGTV9PUci1GRQhUJEJAolJBjXtjmZyQPTaNOgCo9OXMblz89k9Y+7ij5LkW9RRERCVvOEMoy68SwGX3k667bv4YKh0xn66WoOHi66JoMqFCIiUc7MuPTM7CaDnU87kYwpq7j42eks2vxLkWxfhUJEJEZUK1+K4Ve3ZOR1rdix5yCXDJ/B4x8tD3uTQRUKEZEY0/W0k5gyKJ0rUuswMnMd3YZkMnvdT2HbngqFiEgMqlSmBE9c1oI3bjmbI+5c9cJsHnh3Mbv2Hyr0balQiIjEsPYNqzHpnjRu7lCft+Z8y+IwXBwpqdDXKCIiRapsyST+fGEzbmxXjzpVyhb6+rVHISJSTISjSIAKhYiI/AYVChERCUqFQkREglKhEBGRoFQoREQkKBUKEREJSoVCRESCsmi4zF5hM7NtwMbAYiXg2K8qHjuWc7kasD1M0XLLUljPCfa4vO4LZW5yG4vm+Qr1eYU1X7mNa76C3xfv8xXs/kjP18nunvxfo+5erG/AC781lnMZmFeUWQrrOcEel9d9ocxNrM1XqM8rrPn6rfmJ5/nK6754n69g90frfMXDoacPQhjL7THhcDzbCfU5wR6X132hzE1uY9E8X6E+r7DmK7dxzVfw++J9voLdH5XzVSwPPRWEmc1z99RI54gVmq/80Xzlj+Yrf8I1X/GwR5FfL0Q6QIzRfOWP5it/NF/5E5b50h6FiIgEpT0KEREJSoVCRESCUqEQEZGgVCjywcwuMbMXzextM+sS6TzRzswamNnLZvZOpLNEKzMrZ2avBv6urol0nminv6n8KazXrLgpFGY2ysy2mtmSY8a7mdlKM1tjZvcFW4e7v+futwL9gCvDmTfSCmm+1rn7zeFNGn3yOXe9gHcCf1cXF3nYKJCf+YrXv6mc8jlfhfKaFTeFAhgNdMs5YGaJwHCgO9AM6GNmzcwsxcwmHnOrnuOpDwaeV5yNpvDmK96MJsS5A2oDmwIPO1KEGaPJaEKfLzm++SrQa1bS8T4x1rh7ppnVO2a4NbDG3dcBmNkYoKe7Pw5ceOw6zMyAJ4CP3X1BeBNHVmHMV7zKz9wBm8kuFguJrzdu/yef87WsaNNFn/zMl5ktpxBes+LyDzOHWvzn3Rxk/6OtFeTxdwKdgN5m1i+cwaJUvubLzKqa2fPAmWZ2f7jDRbm85m48cJmZPUfRta6IBbnOl/6m8pTX31ehvGbFzR5FYXD3ocDQSOeIFe7+E9nHRiUP7r4HuCnSOWKF/qbyp7Bes+J9j2ILUCfHcu3AmORO83X8NHf5o/nKn7DOV7wXirlAIzOrb2YlgauACRHOFM00X8dPc5c/mq/8Cet8xU2hMLO3gFlAEzPbbGY3u/thYAAwCVgOjHX3pZHMGS00X8dPc5c/mq/8icR8qSmgiIgEFTd7FCIicnxUKEREJCgVChERCUqFQkREglKhEBGRoFQoREQkKBUKEREJSoVCRESCUqGQYsHMHjCzpWa2yMwWmtnZZnaCmd0e4vN35/Z7LDCzema2z8wW5hi7zczczM7JMXZHYKxzkHWVCczfQTOrFuboEiNUKCTmmVlbsq+H0dLdW5DdVnkTcAIQUqEowLbNzKLh39Fadz8jx3IK8A3QFMDMygK3ANuARXmtxN33BdbzXRizSoyJhj9wkYKqAWx39wMA7r7d3b8j+4ItpwTeIT8FYGbXmtmcwNjIwJXBcmVm75nZ/MCeSt8c4/UCl5x8DVhCjq6dgftWmNloM1tlZm+YWSczm2Fmq82sdbD1W/Y1tD80s2/MbImZXZnbWAhz0gIYQ6BQAHcB/way3P3HkGZV5Ch31023mL4B5cm+QtwqYASQHhivByzJ8bhTyb44UInA8gjg+sDvu3M8bnfgZ5XAzzJkF4SqOdabBbTJJUs94DDZ7+gTgPnAKMDIvkLbezke+1/rBy4DXszxmEq5jeWyzSXHjG0DGgCfkL1n9TVwDjA1xDndAFSL9P9b3aLjpj0KiXnuvhtoBfQl+wXybTO7MZeHnh943NzA8fzzyX4xzctdZvYNMJvsvYZGOe7b6O6z83jeendf7O5ZwFLgU3d3YDHZL+rB1r8Y6GxmT5pZR3f/NY+xPJlZHeAnz74sZnXgj8AwoHFgXUcf18TMRpjZ02ZWM9g6Jb7pCndSLLj7EeAL4AszWwzcEFjOyYBX3f03L6EZOAncCWjr7nvN7AugdI6H7Any9AM5fs/KsZxF4N9cXut391Vm1hLoAfzNzD5190dzGwuy/RT+UxB2Ad3IvqbyEGBBYPvVgPuBPwFlgX+a2Y3ufjDIeiVOaY9CYl7gnXHOd/tnABvJfpGskGP8U7KvHVw98LwqZnZyHqutBPwceBFvCrQp5Ni5rj/wzn6vu/8LeApomdvYb6y7Bf8pFE8BAwKFNGcBSQOeAS4FTgTGAacV1n+cFC/ao5DioDwwzMxOIPv8wBqgr7v/FDiJvAT42N3/aGYPApMDn1Q6BNxBdlE51idAPzNbDqwk+/BQYcpr/SnAU2aWFcjXP4+xYFLIfuHH3SfmGG9G9qEwgNVAF3d/GsDMrgCmFfQ/SoonXbhIJMaZWT1gors3z+fz+gOdyT4k9y93H5fjvg1AqrtvL7ykEqtUKERiXODk9UyyT2Cf8VuP/411lSH7MpvJQIq77yiEiBLjVChERCQoncwWEZGgVChERCQoFQoREQlKhUJERIJSoRARkaBUKEREJCgVChERCUqFQkREgvpfRxS+q0XN8uQAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + } + } + ], + "source": [ + "salpeter = PowerLawPDF(gamma=-2.35)\n", + "salpeter.B, _ = integrate.quad(salpeter, a=0.01, b=100.)\n", + "\n", + "m_grid = np.logspace(-2., 2., 100)\n", + "plt.loglog(m_grid, salpeter(m_grid))\n", + "plt.xlabel(r'Stellar mass [$M_{\\odot}$]')\n", + "plt.ylabel('Probability density')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "RotHN5xIXNbP" + }, + "source": [ + "### Quantas estrelas M existem a mais do que estrelas O?\n", + "\n", + "Vamos comparar o número de estrelas anãs M (massa inferior a 60% solar) criadas pelo FMI, com o número de estrelas O (massa superior a 15 vezes solar)." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "id": "OCmyGKTSXNbQ", + "outputId": "7c90dc41-8192-471f-b9fd-f8ae017acc27", + "colab": { + "base_uri": "https://localhost:8080/" + } + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "20936.017868337385\n" + ] + } + ], + "source": [ + "n_m, _ = integrate.quad(salpeter, a=.01, b=.6)\n", + "n_o, _ = integrate.quad(salpeter, a=15., b=100.)\n", + "print(n_m / n_o)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "_yOghSWOXNbQ" + }, + "source": [ + "Existem quase 21.000 estrelas de baixa massa nascidas do que estrelas de alta massa!" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "rdzV9nGgXNbQ" + }, + "source": [ + "### Onde está toda a massa?\n", + "\n", + "Agora vamos calcular as massas totais relativas para todas as estrelas O e todas as estrelas M nascidas. Para fazer isso, pondere o Salpeter IMF em massa (ou seja, adicione um fator extra de massa à integral). Para fazer isso, definimos uma nova função que toma a velha lei de potência IMF como um de seus argumentos. Como esse argumento não é alterado em toda a integral, ele é passado para a tupla `args` dentro de `quad`. É importante que haja apenas *um* argumento que mude sobre a integral, e que seja o *primeiro* argumento que a função que está sendo integrada aceita.\n", + "\n", + "Matematicamente, a integral para as estrelas M é\n", + "\n", + "$$ m^M = \\int_{.01 \\, M_{\\odot}}^{.6 \\, M_{\\odot}} m \\, {\\rm IMF}(m) \\, dm $$\n", + "\n", + "e equivale a ponderar a função densidade de probabilidade (o FMI) em massa. De maneira mais geral, você encontra o valor de alguma propriedade $\\rho$ que depende de $m$ calculando\n", + "\n", + "$$ \\rho(m)^M = \\int_{.01 \\, M_{\\odot}}^{.6 \\, M_{\\odot}} \\rho(m) \\, {\\rm IMF}(m) \\, dm $$" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "id": "QgrkLAWxXNbQ", + "outputId": "40a8267a-0e6c-4edc-886c-17f2240b17e0", + "colab": { + "base_uri": "https://localhost:8080/" + } + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "20.29197629920483" + ] + }, + "metadata": {}, + "execution_count": 9 + } + ], + "source": [ + "def IMF_m(m, imf):\n", + " return imf(m) * m\n", + "\n", + "m_m, _ = integrate.quad(IMF_m, a=.01, b=.6, args=(salpeter, ))\n", + "m_o, _ = integrate.quad(IMF_m, a=15., b=100., args=(salpeter, ))\n", + "\n", + "m_m / m_o" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "AeEYuH2WXNbR" + }, + "source": [ + "\n", + "\n", + "Assim, cerca de 20 vezes mais massa está ligada em estrelas M do que em estrelas O." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "36L_Y2BPXNbR" + }, + "source": [ + "### Extras\n", + "\n", + "* Agora compare a luminosidade total de todas as estrelas O com a luminosidade total de todas as estrelas M. Isso requer uma relação massa-luminosidade, como esta que você usará como $\\rho(m)$:\n", + "\n", + "$$\n", + " \\frac{L}{L_{\\odot}} (M) =\n", + " \\begin{cases} \n", + " \\hfill .23 \\left( \\frac{M}{M_{\\odot}} \\right)^{2.3} \\hfill , \\hfill & .1 < \\frac{M}{M_{\\odot}} < .43 \\\\\n", + " \\hfill \\left( \\frac{M}{M_{\\odot}} \\right)^{4} \\hfill , \\hfill & .43 < \\frac{M}{M_{\\odot}} < 2 \\\\\n", + " \\hfill 1.5 \\left( \\frac{M}{M_{\\odot}} \\right)^{3.5} \\hfill , \\hfill & 2 < \\frac{M}{M_{\\odot}} < 20 \\\\\n", + " \\hfill 3200 \\left( \\frac{M}{M_{\\odot}} \\right) \\hfill , \\hfill & 20 < \\frac{M}{M_{\\odot}} < 100 \\\\\n", + " \\end{cases},\n", + "$$\n", + "\n", + "* Pense em quais estrelas estão produzindo a maior parte da luz e quais estrelas têm a maior parte da massa. Como isso pode resultar em dificuldade em inferir massas estelares a partir da luz que elas produzem? Se estiver interessado em saber mais, consulte este [artigo de revisão](https://ned.ipac.caltech.edu/level5/Sept14/Courteau/Courteau_contents.html)." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "vjo7iJ5IXNbR" + }, + "source": [ + "## Problema desafiador\n", + "\n", + "* No momento, não estamos preocupados com os limites da lei de potência, mas o FMI deve cair para zero de probabilidade em massas abaixo de 0,01 massas solares e acima de 100 massas solares. Modifique `PowerLawPDF` de uma forma que permita entradas `float` e `numpy.ndarray`.\n", + "* Modifique a classe `PowerLawPDF` para usar explicitamente as construções `units` do `astropy`.\n", + "* Deduza uma relação entre a taxa de formação de estrelas recente e a luminosidade de $H\\alpha$. Em outras palavras, encontre um valor de $C$ para a função\n", + "\n", + "$${\\rm SFR \\, [\\frac{M_{\\odot}}{yr}]} = {\\rm C \\, L_{H\\alpha} \\, [\\frac{erg}{s}]} \\, .$$\n", + "\n", + "* Como isso depende da inclinação e dos pontos finais do FMI?\n", + "* Dê uma olhada no Apêndice B de [Hunter & Elmegreen 2004, AJ, 128, 2170](http://adsabs.harvard.edu/cgi-bin/bib_query?arXiv:astro-ph/0408229)\n", + "* Que efeito a alteração do índice da lei de potência ou do limite de massa superior do FMI tem sobre o valor de $C$?\n", + "* Preveja o efeito sobre o valor de $C$ de usar uma forma diferente do FMI, como Kroupa ou Chabrier (ambos são mais leves na extremidade de baixa massa).\n", + "* Se você ainda não está cansado de IMFs, tente definir uma nova classe que implemente uma IMF de lei de potência quebrada (Kroupa) ou log-parabola (Chabrier). Faça os mesmos cálculos acima." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "Sbov-qwiXNbR" + }, + "outputs": [], + "source": [ + "" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython" + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python" + }, + "colab": { + "name": "units-and-integration.ipynb", + "provenance": [] + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file