{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "**This is a fixed-text formatted version of a Jupyter notebook**\n", "\n", "- Try online [![Binder](https://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.16?urlpath=lab/tree/models.ipynb)\n", "- You can contribute with your own notebooks in this\n", "[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\n", "- **Source files:**\n", "[models.ipynb](../_static/notebooks/models.ipynb) |\n", "[models.py](../_static/notebooks/models.py)\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Gammapy Models\n", "\n", "\n", "This is an introduction and overview on how to work with models in Gammapy. \n", "\n", "The sub-package `~gammapy.modeling` contains all the functionality related to modeling and fitting\n", "data. This includes spectral, spatial and temporal model classes, as well as the fit\n", "and parameter API. We will cover the follwing topics in order:\n", "\n", "1. [Spectral Models](#Spectral-Models)\n", "1. [Spatial Models](#Spatial-Models)\n", "1. [SkyModel and SkyDiffuseCube](#SkyModel-and-SkyDiffuseCube)\n", "1. [Model Lists and Serialisation](#Model-Lists-and-Serialisation)\n", "1. [Implementing as Custom Model](#Implementing-a-Custom-Model)\n", "\n", "The models follow a naming scheme which contains the category as a suffix to the class name. An overview of all the available models can be found in the :ref:`model-gallery`.\n", "\n", "Note that there is a separate tutorial [modeling](modeling.ipynb) that explains about `~gammapy.modeling`,\n", "the Gammapy modeling and fitting framework. You have to read that to learn how to work with models in order to analyse data.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Setup" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from astropy import units as u\n", "from gammapy.maps import Map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Spectral Models\n", "\n", "All models are imported from the `gammapy.modeling.models` namespace. Let's start with a `PowerLawSpectralModel`:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from gammapy.modeling.models import PowerLawSpectralModel" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PowerLawSpectralModel\n", "\n", " name value error unit min max frozen\n", "--------- --------- ----- -------------- --- --- ------\n", " index 2.000e+00 nan nan nan False\n", "amplitude 1.000e-12 nan cm-2 s-1 TeV-1 nan nan False\n", "reference 1.000e+00 nan TeV nan nan True\n" ] } ], "source": [ "pwl = PowerLawSpectralModel()\n", "print(pwl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To get a list of all available spectral models you can import and print the spectral model registry or take a look at the [model gallery](https://docs.gammapy.org/stable/modeling/gallery/index.html#spectral-models):" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Registry\n", "--------\n", "\n", " ConstantSpectralModel\n", " CompoundSpectralModel\n", " PowerLawSpectralModel\n", " PowerLaw2SpectralModel\n", " SmoothBrokenPowerLawSpectralModel\n", " ExpCutoffPowerLawSpectralModel\n", " ExpCutoffPowerLaw3FGLSpectralModel\n", " SuperExpCutoffPowerLaw3FGLSpectralModel\n", " SuperExpCutoffPowerLaw4FGLSpectralModel\n", " LogParabolaSpectralModel\n", " TemplateSpectralModel\n", " GaussianSpectralModel\n", " AbsorbedSpectralModel\n", " NaimaSpectralModel\n", " ScaleSpectralModel\n", "\n" ] } ], "source": [ "from gammapy.modeling.models import SPECTRAL_MODELS\n", "\n", "print(SPECTRAL_MODELS)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Spectral models all come with default parameters. Different parameter\n", "values can be passed on creation of the model, either as a string defining\n", "the value and unit or as an `~astropy.units.Quantity` object directly:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "amplitude = 1e-12 * u.Unit(\"TeV-1 cm-2 s-1\")\n", "pwl = PowerLawSpectralModel(amplitude=amplitude, index=2.2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For convenience a `str` specifying the value and unit can be passed as well:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PowerLawSpectralModel\n", "\n", " name value error unit min max frozen\n", "--------- --------- ----- -------------- --- --- ------\n", " index 2.200e+00 nan nan nan False\n", "amplitude 2.700e-12 nan cm-2 s-1 TeV-1 nan nan False\n", "reference 1.000e+00 nan TeV nan nan True\n" ] } ], "source": [ "pwl = PowerLawSpectralModel(amplitude=\"2.7e-12 TeV-1 cm-2 s-1\", index=2.2)\n", "print(pwl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The model can be evaluated at given energies by calling the model instance:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[2.70000000e-12 2.40822469e-13 1.70358483e-14 1.51948705e-15] 1 / (cm2 s TeV)\n" ] } ], "source": [ "energy = [1, 3, 10, 30] * u.TeV\n", "dnde = pwl(energy)\n", "print(dnde)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The returned quantity is a differential photon flux. \n", "\n", "For spectral models you can computed in addition the integrated and energy flux\n", "in a given energy range:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.108034597491956e-12 1 / (cm2 s)\n", "4.982075849517389e-12 TeV / (cm2 s)\n" ] } ], "source": [ "flux = pwl.integral(emin=1 * u.TeV, emax=10 * u.TeV)\n", "print(flux)\n", "\n", "eflux = pwl.energy_flux(emin=1 * u.TeV, emax=10 * u.TeV)\n", "print(eflux)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This also works for a list or an array of integration boundaries:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1.64794383e-12 4.60090769e-13 1.03978226e-13] 1 / (cm2 s)\n" ] } ], "source": [ "energy = [1, 3, 10, 30] * u.TeV\n", "flux = pwl.integral(emin=energy[:-1], emax=energy[1:])\n", "print(flux)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In some cases it can be useful to find use the inverse of a spectral model, to find the energy at which a given flux is reached:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0 TeV\n" ] } ], "source": [ "dnde = 2.7e-12 * u.Unit(\"TeV-1 cm-2 s-1\")\n", "energy = pwl.inverse(dnde)\n", "print(energy)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a convenience you can also plot any spectral model in a given energy range:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3xV9f3H8dcnYa+wQVkBUTAQkK1MqUxlWMC9pSJWFKHOqq1WLY4WFUURq+KoOBjKEAEXS0SWsrcoW4aA7JHP74/EX9M0O/fmJDfv5+NxH3q/595z3mmP+eR7xueYuyMiIpIVUUEHEBGR/EfFQ0REskzFQ0REskzFQ0REskzFQ0REskzFQ0REsqxQ0AFyQ8WKFT02NjboGCIi+crixYv3uHul1JYViOIRGxvLokWLgo4hIpKvmNmPaS2L6MNWZtbTzEYfOHAg6CgiIhEloouHu0929wExMTFBRxERiSgRXTxERCQ8VDxERCTLVDxERCTLIrp46IS5iEh4RHTxyOkJc3dnwpKtHD91OsTJRETyt4guHjm1+MdfGPrB9/QYMZclP/0SdBwRkTxDxSMdzWPL88aNLTh0/BR9X/6ax6as4siJU0HHEhEJnIpHBjrWr8yMIe25plVNXpv7A12fm828DXuCjiUiEigVj0woXawwj18az/sDzqdQVBTX/GsB941bxoGjJ4OOJiISCBWPLGhVpwLTBrdjYIezGLdkK52Hz2LGyp1BxxIRyXUqHllUrHA093evz0d/bEOFUkUZ8PZibn93CXsOHQ86mohIrono4hHO+zziq8cwaVAb7u5yDjNX7qLT8FlMXLoVdw/5tkRE8horCL/smjdv7uFsyb7h51+5d9wylvy0nwvrVeKJ38dTrWzxsG1PRCQ3mNlid2+e2rKInnnklrqVS/PhwNb8tWccCzbto8vwWbw9fzMJCZFfmEWkYFLxCJHoKOOmNrWZMaQ9TWqW4+GPV3Ll6G/YtPtQ0NFEREJOxSPEapQvwdv9W/J0v0as2XmQbs/P4eWvNnLqdELQ0UREQkbFIwzMjMub1+CzoR3oWK8ST326hktfmsfK7WrQKCKRQcUjjCqXKcaoa5vx0jVN2XngGL1enMcz09dw7KQaLYpI/qbiEWZmxsXxZzBzSAcuPa8aI7/cyCUj5rD4x31BRxMRyTYVj1xSrmQR/nl5Y968uSXHTibQb9R8Hpm0ksPH1WhRRPIfFY9c1uGcSkwf0p7rz6/Fm/M30+XZ2cxetzvoWCIiWZIvioeZ1TGz18xsXLKxS83sVTP72My6BJkvq0oVLcSjvRvywa0XULRwFNe//i13f/g9+4+cCDqaiEimhL14mNnrZvazma1IMd7NzNaa2QYzuz+9dbj7Jnfvn2LsI3e/BbgRuCLkwXNBi9jyfHJnO27veBYTl26j87Oz+XTFjqBjiYhkKDdmHmOAbskHzCwaGAl0B+KAq8wszszizWxKilflDNb/UNK68qVihaO5p2t9Pr69DZVKFWXgO0u47Z3F/PzrsaCjiYikqVC4N+Dus80sNsVwS2CDu28CMLP3gN7uPgzokZn1mpkBTwLT3H1J6BIHo2G1GD4e1IbRszfx/Ofr+XrjXh7uEUffptVI/FFFRPKOoM55VAO2JHu/NWksVWZWwcxGAU3M7IGk4TuATkA/MxuYyncGmNkiM1u0e3f+OCFdODqK2zvWZdrgdpxTpRR3f/g917/+LVv2HQk6mojIf8mVrrpJM48p7t4w6f1lQFd3/0PS++uAlu5+Rzi2H+6uuuGQkOD8e8GPPDltDQ7c27Ue118QS1SUZiEikjvyYlfdrUCNZO+rA9tDvZFwPs8j3KKijOsuiGX6kPa0iC3PI5NXcfkr89nwsxotikjwgioeC4Gzzay2mRUBrgQmhXoj7j7Z3QfExMSEetW5pnq5Eoy5qQX/vKwxG3Yf4uLn5zDyyw2cVKNFEQlQblyqOxaYD9Qzs61m1t/dTwGDgOnAauADd18Zhm3n25lHcmZG32bVmTmkA53iKvPM9LX0fnEeK7bl759LRPIvPUkwH/p0xU4e/ngF+w6fYED7Ogy+6GyKFY4OOpaIRJi8eM5DcqBbw6p8NqQDfZtW4+WvNnLx83NYuFmNFkUk90R08YiUw1apiSlRmKf7Nead/q04cTqBy0bN5y8fr+CQGi2KSC7QYasIcOTEKZ6ZvpYxX2/mzJjiPPH7hlxYL6Mb80VE0qfDVhGuRJFC/LVnA8YNbE3xItHc+MZChn7wHb8cVqNFEQmPiC4ekXzYKjXNapVj6p1tGdSxLpO+207nZ2cxddkOCsLsUkRylw5bRahV2w9y3/hlLN92gC5xVXj80oZULlMs6Fgiko/osFUBFHdmGSb+sTX3d6/PrHW7uWj4LD5YuEWzEBEJiYguHgXtsFVKhaKjGNjhLKYNbse5Z5Th3vHLuO41NVoUkZzTYasCIiHBeffbnxj2yWoSHO7pWo8bWscSrUaLIpIGHbYSoqKMa8+vxYyhHWhVpzx/m7KKfqO+Zv2uX4OOJiL5kIpHAVOtbHHeuLEFz11xHpv3HOaSEXMZ8fl6TpxSo0URyTwVjwLIzLi0STVmDu1A14ZVGT5zHb1enMuyrfuDjiYi+UREF4+CfsI8IxVLFeWFq5rw6vXN+eXICS4dOY9hn6zm2MnTQUcTkTxOJ8wFgANHT/LktNWM/XYLsRVK8GTfRpxfp0LQsUQkQDphLhmKKV6YYX0a8e4fWpHgcOXob3hw4nJ+PXYy6GgikgepeMh/aV23Ip/e1Y4/tK3N2G9/osuzs/liza6gY4lIHqPiIf+jRJFCPNQjjvG3taZ0sULcPGYRd723lH1qtCgiSSK6eOiEec40qVmOKXe0Y/BFZzNl2Q46DZ/FpO+3q8WJiOiEuWTOmp0HuW/cMr7feoBO5yY2Wqwao0aLIpFMJ8wlx+pXLcP421rz4MXnMnfDbjo/O4v3vv1JsxCRAkrFQzKtUHQUt7Svw6eD29PgzDLcP2E5V7+6gB/3Hg46mojkMhUPybLYiiUZe8v5DOsTz4ptB+j63Gz+NWcTpxM0CxEpKAqlt9DMlmRiHbvdvWuI8kg+YWZc1bImHetV5sGJy3l86momL9vB030bUa9q6aDjiUiYpVs8gKJAr3SWGzAhdHEkv6kaU4x/3dCcyct28MiklfR4YQ63d6zLHy+sS5FCmtiKRKqMisft7r4xvQ+Y2Z0hzCP5kJnRq/GZtK1bkUcnr+S5z9YzbflOnurXiPNqlA06noiEQbp/Grr7VxmtIDOfCYru88hd5UsW4fkrm/DaDc05cPQkfV6axxNTV3H0hBotikSaDO/zMLMWwLVAO+AM4CiwApgKvOvuef5pQrrPI/f9euwkT05bw78X/ETN8iV4sm88rc+qGHQsEcmCbN/nYWZTgEHALOBSoDbQFHgcKAtMNbMeoY0rkaB0scI88ft4xt5yPmZw9asLeGDCcg6q0aJIREh35mFmVdw93a54ZlbZ3X8OebIQ0swjWEdPnObZz9bxrzmbqFS6KI9fGk/nuCpBxxKRDOTkDvMHzKxleh/I64VDgle8SDR/vvhcJv6xDeVKFOGWtxZxx9il7D10POhoIpJNGRWPLcBIM9toZk+YWcPcCCWRqXGNskwa1JYhnc7h0xWJjRY//m6bWpyI5EMZXW31T3dvAXQBjgBjzWyFmf3ZzOrkSkKJKEUKRTG409lMvbMdtSqUZPB739H/zUVs33806GgikgVZ7qprZs2AfwGN3D06LKlCTOc88qbTCc6Yrzfzj+lriY4yHri4Ple1qElUlAUdTUQIQVddM4s2s+5m9iaJl+huAq4IYUYpgKKjjP5tazP9rvY0rhHDgxNXcNWr3/DDHjVaFMnrMrpUt6OZjQa2AXcCXwBnu3tfdx+XGwEl8tWsUIJ3+rfiqb7xrNpxkG7PzeaVWRs5dToh6GgikoaMZh5/A5YC8e7e3d3fzO2bAs2sjpm9Zmbjko2da2ajzGycmd2Wm3kkPMyMK1rU5LOhHWh/TiWGTVtDn5e/ZvWOg0FHE5FUZHTCvJ27v+zuu83sfDO7HsDMKphZzYxWbmavm9nPZrYixXg3M1trZhvM7P4MMmxy9/4pxla7+0DgciDV43GSP1UpU4zR1zXjxaubsH3/UXq+MJfhM9Zy/JRanIjkJZk95/EQ8FfgoaShYsC7mfjqGKBbinVFAyOB7kAccJWZxZlZvJlNSfGqnE6mXsBc4PPM/AySf5gZPRqdycwhHejV+ExGfLGBS0bMZfGPvwQdTUSSZLZndj/gYuAwgLtvA8pk9CV3nw3sSzHcEtiQNKM4AbwH9Hb35e7eI8UrzRsQ3X2Su7cGrsnkzyD5TLmSRRh+xXm8cVMLjhw/Rb9RX/Po5JUcOXEq6GgiBV5mi8dxT7ym1wHMrEQOtlmNxJsPf7M1aSxVSYfIRgFNzOyBpLELzWyEmb0CfJLG9waY2SIzW7R79+4cxJWgdaxXmelD2nNtq1q8MW8zXZ6dzdz1e4KOJVKgZfQ8j99MMLORQIyZ3QT0B17P5jZTu4g/zZtN3H0vMDDF2FfAV+ltxN1HA6Mh8T6PrIaUvKV0scI8dmlDejQ6g/snLOfa1xZwefPqPHhJHDHFCwcdT6TAydTMw92fAqYAk4DGwBPu/lw2t7kVqJHsfXVgezbXlS49zyPytKpTgWmD2zGww1mMX7KNzsNnMX3lzqBjiRQ4GXXVneHuXXK0AbNYYIq7N0x6XwhYB1xE4v0jC4Gr3X1lTraTHt1hHpmWbz3AveOXsXrHQS6JP4NHejWgUumiQccSiRg5ucO8Ug43PBaYD9Qzs61m1t/dT5H4jJDpwGrgg3AVDs08Ilt89RgmDWrDPV3rMXPVLjoNn8X4xVvVaFEkF2Q089gE3J3WcnefEI5QoaaZR+Tb8POv3DtuGUt+2k+Hcyrx9z7xVCtbPOhYIvlaejOPjIrHXuBj0jjJ7e43hyZieKl4FAynE5y352/m6elrMeC+7vW5tlUtNVoUyaacFI8l7t40bMnCzMx6Aj3r1q17y/r164OOI7lky74j/Hnicuas30OL2HI82bcRZ1UqFXQskXwnJ+c88vWfbO4+2d0HxMTEBB1FclGN8iV46+aWPNOvEWt3/kr35+fw0lcb1GhRJIQyKh7X5UoKkRAzMy5rXoPP/tSB39WrzNOfruXSl+axcrsunhAJhYwaI65Ib3lep6utpHLpYoy6rhkvX9OUnQeO0+vFeTwzfQ3HTqrRokhOZPlJgvmRTpgLwP4jJ3h86mrGLd7KWZVK8lTfRjSPLR90LJE8K8dPEkyxsnJm1ijnsURyV9kSRfjHZY156+aWHDuZwGWvzOeRSSs5fFyNFkWyKrMt2b8yszJmVh74HnjDzIaHN1rO6bCVpKb9OZWYMaQ9N1wQy5vzExstzl6n5pkiWZHZmUeMux8E+gBvuHszoFP4YoWGrraStJQsWohHejXgw1svoGjhKK5//Vvu/vB79h85EXQ0kXwhs8WjkJmdQeKT+6aEMY9IrmoeW55P7mzH7R3PYuLSbXQaPptpy3cEHUskz8ts8fgbib2oNrj7QjOrA+iuO4kIxQpHc0/X+kwa1IYqZYpy27+XMPDtxfx88FjQ0UTyLF1tJZLMydMJ/GvODzz72TqKFYri4R5x9GtWHbN8fb+sSLaE9Gqr/EQnzCWrCkdHcduFZzFtcDvqVy3DPeOWcf3r37Jl35Ggo4nkKZp5iKQhIcH594IfeXLaGhy4t2s9rr8gVo0WpcAosDMPkZyIijKuuyCWGUM70CK2PI9MXsVlr8xnw8+/Bh1NJHAZFg8zq29mF5lZqRTj3cIXSyTvqFa2OGNuasHwyxuzcfchLn5+Li9+sZ6TarQoBVi6xcPM7iTxeR53ACvMrHeyxX8PZzCRvMTM6NO0OjOHdKBzgyr8Y8Y6er04jxXbdD5NCqaMZh63AM3c/VLgQuBhMxuctCzPH/jVCXMJtUqlizLy6qa8cl0z9h46Tu+R83hymhotSsGT0cOgVrl7XLL3pYBxwCrgd+5+Xvgj5pxOmEs4HDhykr9/spr3F22hTsWSPNm3ES1rq9GiRI6cnDDfaWb/XyDc/RDQA6gIxIcuokj+E1OiME/1a8Q7/VtxMiGBy1+Zz8MfreCQGi1KAZBR8bge2Jl8wN1Pufv1QPuwpRLJR9qeXZHpd7XnpjaxvLPgR7oMn8WXa38OOpZIWGX0MKit7r4T/tOK3cyamllT4GiuJBTJB0oUKcRfezZg3MDWlChaiJveWMjQ97/jl8NqtCiRqVBmPmRmjwE3AhuB306SOPC78MQSyZ+a1SrH1DvbMvKLDbz01UZmr9/No70acnF8VbU4kYiSqTvMzWwtEO/u+fLPKJ0wlyCs3nGQe8ctY/m2A3SJq8JjlzakSpliQccSybRQ3GG+Aigbukgike/cM8ow8Y+teaB7fWat202n4bP4YOEWCkJLIIl8mZ15NCfxZsEVwPHfxt29V/ii5ZyZ9QR61q1b95b169VBXoLzw57D3Dd+Gd/+sI+2dSsyrE88NcqXCDqWSLrSm3lktnisBF4BlgP/35PB3WeFKmQ46bCV5AUJCc673/7Ek9PWcDrBuadrPW5oHUu0Gi1KHhWK4jHL3TuEPFkuUfGQvGT7/qM8OHE5X67dTZOaZXm6byPOrlI66Fgi/yMU5zwWm9kwM7vgt0t1ky7XFZEsOrNscV6/sQXPX3kem/cc5pIRcxnx+XpOnFKjRck/MnWpLtAk6Z/nJxvTpboi2WRm9D6vGm3qVuTRyasYPnMdnyzfwVN9G9G4hq5NkbxPD4MSyQNmrtrFQx8tZ/evx7mlXR3u6nQOxYtEBx1LCrgcH7Yys7+bWdlk78uZ2eOhCihS0HWOq8LMoR24okUNXpm9ie7Pz+abTXuDjiWSpsye8+ju7vt/e+PuvwAXhyeSSMFUplhhhvVpxLt/aEWCw5Wjv+HBicv59djJoKOJ/I/MFo9oMyv62xszKw4UTefzIpJNretW5NO72vGHtrUZ++1PdHl2Nl+s2RV0LJH/ktni8Q7wuZn1N7ObgZnAm+GLJVKwlShSiId6xDH+ttaUKlqIm8cs4q73lrJPjRYlj8hU8XD3p4HHgXOBBsBjSWNhZ2Z1zOw1MxuXYrykmS02sx65kUMkCE1qlmPqne0YfNHZTF2+g07DZzHp++1qcSKBy+gZ5v9/66u7f+rud7v7n9x9emqfSeX7r5vZz2a2IsV4NzNba2YbzOz+9DK4+yZ375/KovuAD9L7rkgkKFIoiiGdz2HyHW2pUa44d45dyi1vLWLngWNBR5MCLKOZx5dmdoeZ1Uw+aGZFzOx3ZvYmcEM63x8DdEvx3WhgJNAdiAOuMrM4M4s3sykpXpVTW6mZdSLxUbg6ECwFRv2qZZjwxzY8ePG5zN2wh87DZzH22580C5FAZHSTYDfgZmCsmdUG9gPFSSw6M4Bn3f27tL7s7rPNLDbFcEtgg7tvAjCz94De7j6MxEfcZkZHoCSJxeeomX3i7ro9VyJedJRxS/s6dGlQhfvHL+eBCcuZ9N12nuwbT60KJYOOJwVIRk8SPObuL7l7G6AWcBHQxN1rufst6RWOdFQDtiR7vzVpLFVmVsHMRgFNzOyBpFwPuvtdwLvAq6kVDjMbYGaLzGzR7t27sxFTJO+qVaEk797SimF94lmx7QBdn5vNv+Zs4nSCZiGSOzLbngR3PwnsCME2UztHkuYe7+57gYFpLBuTzvdGA6Mh8Q7zrEUUyfvMjKta1qRjvco89NFyHp+6msnLdvB030bUq6pGixJemb1UN5S2AjWSva8ObA/Hhsysp5mNPnDgQDhWL5InVI0pxqvXN2fEVU3Ysu8IPV6Yw7Mz16nRooRVEMVjIXC2mdU2syLAlcCkcGzI3Se7+4CYmJhwrF4kzzAzejU+k8+GduDi+DN4/vP19HxhLt9t2Z/xl0WyIaNLdaeb2RAzq5+dlZvZWGA+UM/MtppZf3c/BQwCpgOrgQ/cfWV21p+J7WvmIQVK+ZJFeP7KJrx+Y3MOHjtJn5fm8fiUVRw5cSroaBJh0u2qa2ZVSbziqhtwDrAA+BT43N0P5UrCEFBXXSmIfj12kmHT1vDugp+oWb4ET/aJp3XdikHHknwkx08STFpJFNCKxPszLgKOAjNy607znFDxkILsm017uX/8MjbvPcKVLWrw50vOpUyxwkHHknwgJMUjlZVWBLq6+79zEi6czKwn0LNu3bq3rF+/Pug4IoE5euI0z322jlfnbKJS6aI8cWk8neKqBB1L8riwFI/8RDMPkUTfb9nPfeOXsWbnr/RsfCZ/7RlHxVJqkC2pC8UzzEUkAjSuUZZJg9oytPM5fLpiB52Hz2Li0q1qcSJZFtHFQ1dbifyvIoWiuPOis5l6ZztiK5ZkyPvfc/OYhWzffzToaJKPZLt4mNlNoQwSDrrPQyRt51QpzbiBrflLjzi+2bSPLs/O5p1vfiRBLU4kE3Iy83g0ZClEJBDRUcbNbWsz/a72NK4Rw0MfreDKV7/hhz2Hg44meVxG93ksS2sRcI675+kzbbraSiTz3J0PFm3h8amrOXEqgaGdz6F/29oUio7oo9uSjmxfbWVmu4CuwC8pFwFfu/uZIUsZRrraSiTzdh08xsMfrWDGql00qh7Dk30aEXdmmaBjSQBycrXVFKCUu/+Y4rUZ+CrEOUUkD6hSphivXNeMkVc3Zfv+o/R6cS7/nLGW46dOBx1N8hDd5yEiafrl8Akem7qKCUu2UbdyKZ7q24hmtcoFHUtyie7zEJFsKVeyCMMvP483bmrBkeOn6Dfqax6dvJLDx9VosaDLqKvukoxWkJnPBEX3eYiERsd6lZkxtAPXnV+LN+Ztputzs5mzXk/oLMgyOmF+FEjvMiUDYty9ZqiDhZIOW4mEzrc/7OP+8cvYtOcwlzevzoMXxxFTQo0WI1F6h60yegxtZp7jobNoIgVIy9rl+WRwO0Z8vp5XZm/iy7W7eax3Q7o1rBp0NMlFOmEuItm2YtsB7h23jFU7DnJJ/Bk80qsBlUrn6du/JAt0wlxEwqJhtRg+HtSGe7rWY+aqXXQaPovxi9VosSCI6OKhE+Yi4Vc4OorbO9blk8FtqVu5FH/68HtufGMh29RoMaJlqniYWVwqYxeGPE2IqTGiSO6pW7k0H956AY/0jGPh5n10GT6Lt+ZvVqPFCJXZmccHZnafJSpuZi8Aw8IZTETyn6go48Y2iY0Wm9Yqx18+XskVo+ezcfehoKNJiGW2eLQCagBfAwuB7UCbcIUSkfytRvkSvHVzS/5xWWPW7TpE9+fnMPLLDZw8nRB0NAmRzBaPk8BRoDhQDPjB3bUXiEiazIx+zaozc2h7LqpfmWemr+XSkfNYsU3nICNBZovHQhKLRwugLXCVmY0LWyoRiRiVSxfj5Wub8fI1Tdl18Di9R87jmelrOHZSt4jlZ5ktHv3d/S/uftLdd7p7b+DjcAYTkcjSPf4MPhvant83qcbILzdy8Yg5LNq8L+hYkk2ZLR4/m1nN5C9gVjiDiUjkKVuiCP+4rDFv3dyS4ycTuOyV+TwySY0W86NM3WFuZssBJ7GXVTGgNrDW3RuEN17O6EmCInnX4eOneGb6Wt6cv5kzY4ozrE887c+pFHQsSSbbTxJMZ4VNgVvd/dachssNak8iknct2ryPe8cvY9Puw/RrVp2HLjmXsiWKBB1LCEN7EndfQuLJcxGRHGkeW55P7mzHoI51mbh0G52Gz2ba8h1Bx5IMZNRVFwAzG5rsbRTQFFAzfxEJiWKFo7m7az26x1fl3nHLuO3fS+jWoCp/692AymWKBR1PUpHZmUfpZK+iwFSgd7hCiUjB1ODMGD6+vQ33davPF2t/ptPwWXy4aIsaLeZBaskuInnSxt2HuH/8MhZu/oV2Z1fk77+Pp0b5EkHHKlCyfcLczCaTeJVVqty9V87jhZ+Kh0j+lJDgvLPgR56atgYH7u1aj+sviCUqyoKOViDkpHh0SG/F7p4v7vVQ8RDJ37btP8qfJyxn1rrdNKtVjqf6xlO3cumgY0W8nBSPmu7+U9iS5RIVD5H8z92ZuHQbf5uyiiPHT3PnRXW5tcNZFI6O6McSBSonl+p+lGwl40OaSkQkC8yMPk2rM3NIBzo3qMI/Zqyj14tqtBiUjIpH8gOLdcIZJM0AZnXM7LXkjRjN7EIzm2Nmo/LDQ6lEJHQqlS7KyKubMvq6Zuw9lNho8clparSY2zIqHp7Gv2eKmb1uZj+b2YoU493MbK2ZbTCz+9MN4L7J3funkusQia1StmY1l4jkf10aVGXm0A70a1qdUbM20v35OSzYtDfoWAVGRuc8TgOHSZyBFAeO/LYIcHcvk+7KzdqT+Ev+LXdvmDQWDawDOpP4i38hcBUQzf8+nfBmd/856Xvj3L1f0r9HuXuCmVUBhrv7Nenl0DkPkcg2b8Me7p+wjC37jnLd+bW4t1s9ShcrHHSsfC+9cx7p3mHu7tE52bC7zzaz2BTDLYEN7r4pKdx7QG93Hwb0yOR6f3sQ1S8k3rQoIgVYm7oVmX5Xe/45Yx2vz/uBz1fv4onfx9OxfuWgo0WsIC5TqAZsSfZ+a9JYqsysgpmNApqY2QNJY33M7BXgbeDFNL43wMwWmdmi3bvVSUUk0pUoUoiHe8Qx/rbWlCxaiJvGLGTI+9+x7/CJoKNFpEz1tgqx1O7uSe9GxL3AwBRjE4AJ6W3E3UcDoyHxsFXWY4pIftS0Zjmm3NmWkV9u5KUvNzB73W4e7d2AS+LPwEw3F4ZKEDOPrUCNZO+rA9vDsSEz62lmow8c0KV8IgVJ0ULRDO18DpPvaEu1csUZ9O5Sbn17MbsOHgs6WsQIongsBM42s9pmVgS4EpgUjg25+2R3HxATExOO1YtIHnfuGWWYcFtr/nxxfWat202n4bN4f+FParQYAmEtHmY2FpgP1DOzrWbW391PAYOA6cBq4AN3Xxmm7WvmIVLAFYqOYkD7s/j0rvbEnVGG+8Yv59rXFvDT3iMZf1nSpK66IlJgJCQ4Yxf+xLBP1nA6wbm7az1ubB1LtBotpirkTxIUEcmPoqKMa1rVYubQ9lxwVriHatgAAAvISURBVAUem7KKvi9/zbpdvwYdLd+J6OKhw1YikpozYorz2g3Nef7K8/hx72EuGTGH5z9bz4lTCRl/WQAdthKRAm7voeM8MnkVk7/fTv2qpXmqbyMa1ygbdKw8QYetRETSUKFUUV64qgmvXt+cX46c4PcvzePvn6zm6Ak1WkxPRBcPHbYSkczqHFeFmUM7cEWLmoyevYnuz89m/kY1WkxLRBcP3echIllRplhhhvWJ591bWuHAVa9+w58nLufgsZNBR8tzIrp4iIhkR+uzKvLp4PYMaF+H9779iS7DZ/P56l1Bx8pTIrp46LCViGRX8SLR/Pnic5nwxzbEFC9M/zcXMfi9pew9dDzoaHmCrrYSEcnAiVMJvPzVRl78cj2lixXmrz3j6NX4zIhvtKirrUREcqBIoSgGdzqbKXe0o0b5Egx+7ztueWsROw8U3EaLKh4iIplUr2ppJtzWmocuOZe5G/bQefgsxn5bMBstRnTx0DkPEQm16CjjD+3qMP2u9jSsFsMDE5Zz9asL2LzncNDRcpXOeYiIZJO7897CLfx96mpOJiTwp871uLlt7YhptKhzHiIiYWBmXNWyJjOHdqBt3Uo88clq+rw0jzU7DwYdLexUPEREcqhqTDFevb4ZL1zVhK2/HKXHiLkMn7mO46cit8WJioeISAiYGT0bn8nMoR3o2fhMRny+np4vzGXpT78EHS0sVDxEREKofMkiPHvFebxxYwt+PXaKPi9/zWNTVnHkxKmgo4VURBcPXW0lIkHpWL8yM4a05+qWNXlt7g90e24OX2/YE3SskIno4qHGiCISpNLFCvPE7+N5b8D5REcZV/9rAfePX8aBo/m/0WJEFw8Rkbzg/DoVmDa4Hbd2qMMHi7bQefgsZqzcGXSsHFHxEBHJBcUKR/NA93P56PY2lC9ZhAFvL2bQu0vYk08bLap4iIjkokbVyzL5jrbc3eUcZqzcRafhs5i4dGu+a3Gi4iEikssKR0cx6HdnM/XOttSuWJIh73/PTWMWsm3/0aCjZZqKh4hIQM6uUppxA1vzlx5xLNi0jy7DZ/H2Nz+SkJD3ZyEqHiIiAYqOMm5uW5vpd7XnvJplefijFVz56jf8kMcbLUZ08dB9HiKSX9SsUIJ3+rfi6b6NWL3jIN2em82oWRs5dToh6GipUlddEZE8ZtfBYzz80QpmrNpFfLUYnu7XiHPPKJPrOdRVV0QkH6lSphivXNeMl65pyo4DR+n5wlz+OWNtnmq0qOIhIpIHmRkXx5/BzCEd6HXembzwxQYuGTGXxT/mjUaLKh4iInlYuZJFGH75eYy5qQVHT5ym36iveXTySg4fD7bRooqHiEg+cGG9ykwf0p7rzq/FG/M20/W52cxdH1yjRRUPEZF8olTRQvytd0M+uPUCikRHce1rC7h33PccOJL7jRZVPERE8pmWtcvzyeB23HbhWYxfso1Oz87i0xW522hRxUNEJB8qVjia+7rV5+Pb21CpVFEGvrOY2/+9hN2/5k6jRRUPEZF8rGG1GD4e1IZ7utZj5qrERovjF4e/0WKeLx5mVsfMXjOzccnGoszsCTN7wcxuCDKfiEjQCkdHcXvHunwyuB11K5fiTx9+zw1vLGTrL0fCts2wFg8ze93MfjazFSnGu5nZWjPbYGb3p7cOd9/k7v1TDPcGqgEnga2hTS0ikj/VrVyKD2+9gEd7NWDR5n10eXY2b83fHJZthXvmMQbolnzAzKKBkUB3IA64yszizCzezKakeFVOY731gPnuPhS4LYz5RUTylago44bWscwY0p7mseVZu/PXsGynUFjWmsTdZ5tZbIrhlsAGd98EYGbvAb3dfRjQI5Or3gqcSPr3vHO/vohIHlG9XAnevKkFJ8LUWDGIcx7VgC3J3m9NGkuVmVUws1FAEzN7IGl4AtDVzF4AZqfxvQFmtsjMFu3evTtE0UVE8g8zo2ih6LCsO6wzjzRYKmNpXhbg7nuBgSnGjgApz4Ok/N5oYDQkdtXNekwREUlLEDOPrUCNZO+rA9vDsSE9z0NEJDyCKB4LgbPNrLaZFQGuBCaFY0PuPtndB8TExIRj9SIiBVa4L9UdC8wH6pnZVjPr7+6ngEHAdGA18IG7rwzT9jXzEBEJAz1JUEREUlVgnySomYeISHhEdPHQOQ8RkfAoEIetzGw3sB9IbwoSk87yikBwT13JnvR+nry8rZysK6vfzeznM/O5jD4TafsX5N4+pv0ruP2rlrtXSnWJuxeIFzA6u8uBRUHnD/XPm1e3lZN1ZfW7mf18Zj5X0PavUP//nlvb0f4VuldEH7ZKYXIOl+c3ufnzhHJbOVlXVr+b2c9n5nMFbf+C3PuZtH/lwf2rQBy2yikzW+RpXHEgklPavyScwrV/FaSZR06MDjqARDTtXxJOYdm/NPMQEZEs08xDRESyTMVDRESyTMVDRESyTMUji8yspJm9aWavmtk1QeeRyGJmdczsNTMbF3QWiUxmdmnS76+PzaxLdtej4gGY2etm9rOZrUgx3s3M1prZBjO7P2m4DzDO3W8BeuV6WMl3srJ/ufsmd0/3QWciKWVxH/so6ffXjcAV2d2mikeiMUC35ANmFg2MBLoDccBVZhZH4sOrfnuMrp6fLpkxhszvXyLZMYas72MPJS3PFhUPwN1nA/tSDLcENiT9JXgCeA/oTeKTEKsnfUb/+0mGsrh/iWRZVvYxS/QUMM3dl2R3m/rll7Zq/GeGAYlFoxowAehrZi8TmS0nJHekun+ZWQUzGwU0MbMHgokmESKt32F3AJ2AfmY2MLsrL5SzbBHNUhlzdz8M3JTbYSTipLV/7QWy/R+0SDJp7WMjgBE5XblmHmnbCtRI9r46sD2gLBJ5tH9JuIV1H1PxSNtC4Gwzq21mRYArgUkBZ5LIof1Lwi2s+5iKB2BmY4H5QD0z22pm/d39FDAImA6sBj5w95VB5pT8SfuXhFsQ+5gaI4qISJZp5iEiIlmm4iEiIlmm4iEiIlmm4iEiIlmm4iEiIlmm4iEiIlmm4iEFnpmdNrPvkr3uz/hb4Wdmm81suZk1N7OJSdk2mNmBZFlbp/HdP5jZ2ynGqiS17S5sZu+b2T4zuzR3fhqJNLrPQwo8Mzvk7qVCvM5CSTdp5WQdm4Hm7r4n2diFwN3u3iOD75YD1gPV3f1Y0tggIN7db016/w6Jz6b5KCc5pWDSzEMkDUl/+T9qZkuSZgD1k8ZLJj18Z6GZLTWz3knjN5rZh2Y2GZhhZlFm9pKZrTSzKWb2iZn1M7OLzGxisu10NrMJOcjZwsxmmdliM5tmZlXc/Rfga+CSZB+9Ehib3e2IJKfiIQLFUxy2Sv50tT3u3hR4Gbg7aexB4At3bwF0BJ4xs5JJyy4AbnD335H41MlYIB74Q9IygC+Ac82sUtL7m4A3shPczIoCzwN93b0Z8A7wWNLisSQWDMysRlKW2dnZjkhKaskuAkfd/bw0lv02I1hMYjEA6AL0MrPfikkxoGbSv890998eytMW+NDdE4CdZvYlJPbETjofca2ZvUFiUbk+m9nPBRoAn5kZQDSJ3VQhsQneCDMrReLjRj9IyiKSYyoeIuk7nvTP0/znvxcj8S/9tck/aGatgMPJh9JZ7xskPkzsGIkFJrvnRwxY5u7tUi5w98Nm9hmJTyi8Ergtm9sQ+R86bCWSddOBOyzpT30za5LG5+aS+NTJKDOrAlz42wJ3307isxUeIvH509m1isQnELZMylLEzBokWz4WuAco6+4Lc7Adkf+i4iHyv+c8nszg848BhYFlZraC/5xjSGk8iYeQVgCvAAuAA8mW/xvY4u6rshvc3Y8D/YDhZvY9sBRolewjn5J4SO297G5DJDW6VFckjMyslLsfMrMKwLdAG3ffmbTsRWCpu7+Wxnc3k+JS3RBn06W6km2aeYiE1xQz+w6YAzyWrHAsBhqReHVUWnYDn5tZ81CHMrP3gTYknnMRyTLNPEREJMs08xARkSxT8RARkSxT8RARkSxT8RARkSxT8RARkSxT8RARkSz7P7WVvRzyjSY2AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pwl.plot(energy_range=[1, 100] * u.TeV)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Spatial Models" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Spatial models are imported from the same `gammapy.modeling.models` namespace, let's start with a `GaussianSpatialModel`:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "from gammapy.modeling.models import GaussianSpatialModel" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "GaussianSpatialModel\n", "\n", " name value error unit min max frozen\n", "----- --------- ----- ---- ---------- --------- ------\n", "lon_0 0.000e+00 nan deg nan nan False\n", "lat_0 0.000e+00 nan deg -9.000e+01 9.000e+01 False\n", "sigma 2.000e-01 nan deg 0.000e+00 nan False\n", " e 0.000e+00 nan 0.000e+00 1.000e+00 True\n", " phi 0.000e+00 nan deg nan nan True\n" ] } ], "source": [ "gauss = GaussianSpatialModel(lon_0=\"0 deg\", lat_0=\"0 deg\", sigma=\"0.2 deg\")\n", "print(gauss)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again you can check the `SPATIAL_MODELS` registry to see which models are available or take a look at the [model gallery](https://docs.gammapy.org/stable/modeling/gallery/index.html#spatial-models)." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Registry\n", "--------\n", "\n", " ConstantSpatialModel\n", " TemplateSpatialModel\n", " DiskSpatialModel\n", " GaussianSpatialModel\n", " PointSpatialModel\n", " ShellSpatialModel\n", "\n" ] } ], "source": [ "from gammapy.modeling.models import SPATIAL_MODELS\n", "\n", "print(SPATIAL_MODELS)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The default coordinate frame for all spatial models is ``\"icrs\"``, but the frame can be modified using the\n", "``frame`` argument:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "gauss = GaussianSpatialModel(\n", " lon_0=\"0 deg\", lat_0=\"0 deg\", sigma=\"0.2 deg\", frame=\"galactic\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can specify any valid `~astropy.coordinates` frame. The center position of the model can be retrieved as a `~astropy.coordinates.SkyCoord` object using `SpatialModel.position`: " ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "print(gauss.position)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Spatial models can be evaluated again by calling the instance:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[13061.88470839 10172.60603928] 1 / sr\n" ] } ], "source": [ "lon = [0, 0.1] * u.deg\n", "lat = [0, 0.1] * u.deg\n", "\n", "flux_per_omega = gauss(lon, lat)\n", "print(flux_per_omega)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The returned quantity corresponds to a surface brightness. Spatial model\n", "can be also evaluated using `gammapy.maps.Map` and `gammapy.maps.Geom` objects:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAEMCAYAAAAPqefdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2de7xkRXXvv7+ZYV7AyGMAEZAhkcgFoiiI+Li5fERkUK8QFUQjDkokckEIJlch8QYf4XMJisYHoojISDSASAQNiAQuGiMPB4IgqGEiKCMThuE1zIzMMDPr/rGrmT7V1b2r99ndZ5/T6/v59Od0VVfVrv04q6vXqrWWzAzHcRxncjBtoifgOI7j5ONC23EcZxLhQttxHGcS4ULbcRxnEuFC23EcZxLhQttxHGcSMVShLWmhpF9KWirp9FD3PEk3SrpK0lah7n2S7pZ0p6QfSdq7bYxFku4Lr0Vt9TdJWjDM83Ecxxk2QxPakqYD5wGHA3sDbw/C+BTg/cCFwDtD82+Y2R+a2X7AOcCnwhjbAWcCLwcOBM6UtG2f87ikhtNpHFP1vGDqnpufl1OFYa60DwSWmtmvzGw9cClwBDAd2BReAjCzVW39tgRaHkCHAdeb2WNm9jhwPbAwfPYYsDFjHm8e74k0lKl6XjB1z83Py+mbGUM81i7Ag23lZRQr5k8AlwBPAu9ofSjpJOADwEzgNT3G2AXAzPxBcRxnyjNMoa1EnZnZr4E/SnxwHnCepHcAHwYWdRuj9MDFz7WWUJ+7/fbbN9J3P3VyucycOZP5Xc5rPOP2YlDjxsyaOZMdh3DPBnWAbuP2umfjGXeikYSkNW1VV5rZsf2MsXDhQlu5cmVW29tvv/06M1tY3nJqMEyhvQzYra28K/BQRr9LgfPbxjg4GuOmsgHCA3MsgCR77LHHMg7bHzl6prI2VcZolWcAv+tyXt369NumrE+VNjljbAtQcs82lYxR9nnVNqk+OW1g7D3L7dNPm7rOuV923313HnjggS3HM8bKlStZsmRJVltJ88dzrMnGMHXaPwH2lLSHpJnAMcDVqYaS9mwrvgG4L7y/DnidpG2DAfJ1oc5xnCnHpszXaDG0lbaZbZB0MoWQnQ5cZGb3dGl+sqTXAs8Aj1OoRjCzxyR9nOILAOBjZlb/stlxnAnGgA0TPYlGMkz1CGZ2DXBNRrtTe3x2EXBRnfPql0GpDnLUFjmfl407rD4phvXTbliqjrr6xMTXqcr5pK51lTYTR3Nm0iSGKrQdx3HyMFxop3Gh7ThOQ3GhncKFtuM4DcRX2t1woe04TkNxoZ3ChXYGdRkVq/Spw0AY3+RhGiKrXJc6qGs/dR2GyHgPxHj2dtdN6jhVjKD147tHuuFC23GchuIr7RQutB3HaSCu0+6GC23HcRqKC+0ULrQT1OHQktOvTPdcV58cnXZZmypzy20zCOrSaZfpo1Na1zIHlirHyelTNo/UXHL6TYyO21fa3fB0Y47jNJQNma/eSLpI0gpJP2ur+4SkX0i6S9I/Sdqm7bMzQnatX0o6rK1+/5BRa6mkz0pSqJ8l6bJQf+ugM2i50HYcp4G0Vtq1BIy6mM3JUlpcD+xrZi8C/gM4AyBk0zoG2Cf0+ULIugVFtNETgD3DqzXm8cDjZvYC4NPA3/V1qn3iQttxnIZSj9A2sx9SZLZqr/u+mbWW6bdQhHmGIpvWpWa2zszuB5YCB0raGZhnZjebmQFfA45s67M4vL8COKS1Ch8EI6/THtZ+6pw2ObrmKvrpsnJdx6nrnKswiP3UOW1SfeK6+BxTfcrapPrE17/KOecwMUGmhqrTfg9wWXi/C4UQb9HKjvVMeB/Xt/o8CM9GM30S2B7Iy+LQJyMvtB3HaSrZQnu+pPaMCReY2QU5HSX9NcV34tdbVYlm1qO+V5+B4ELbcZyGki20V5rZAf2OLmkR8EbgkKDygO4ZtpaxWYXSXt/eZ5mkGcBziNQxdeI6bcdxGkjLjX38u0dSSFoIfAh4k5mtbfvoauCYsCNkDwqD421mthx4StJBQV/9LuCqtj6Lwvu3Aje2fQnUjq+0HcdpIPXptCX9I0Vu2fmSlgFnUuwWmQVcH2yGt5jZ+8zsHkmXA/dSfCOcZGYbw1AnUuxEmQNcG14AXwEukbSUYoV9TC0T78JICu1+nWcGZYgclFFx5pD65LQpuy51mdjjZU0V55Qco2JcXj+gPjnkGFIHQVWnnXqO1D9m9vZE9Vd6tD8LOCtRvwTYN1H/NHDUeObYDyMptB3HmQy4R2QKF9qO4zQQd2Pvhgttx3EaiMfT7sbIC+1hZUCHepxeynTPqTbD6jM90WaLqDxROu1nEn02RuUcXXPcJnVdyvqkno06+vT7eTdynGnK+tSDr7RTjLzQdhynqbjQTuFC23GcBuI67W640HYcp6G40E7hQttxnAbiK+1ujKTQ7tcRoA7HmVRdmZExp02qz+wa+sQGxFmJPlUMkXGbVJ8qxEbFuFzVELmuZJynE33KjJWp5yd17LI+VcgxTlbJVDMY5xrfPZJiJIW24zhNx1fa3XCh7ThOQ3GhncKFtuM4DcRX2t0YeaE9rAwzqbocXXNZm1gXDZ3657hNqk9Zm5ROO26TOudYp13mbAMwrURBuinxv1zmTJPSacca05R+OtZpx23i80m1ifX2VZ65Mp13VXKy6OTgzjXDY+SFtuM4TcWFdgoX2o7jNBCPPdINF9qO4zQQ12l3Y+SEthirs6srWUHOPu0y/XTO/ukc/fTckjbx51Ck4ujVZ8tEn5nRSc9OTGZmdFLTowszPXExq+i0N0Z1G6NF2vqEUvjpSPm8PjHumrhPVP5dYn5levzU3vRBZamPyQkGFT+7E7fedaGdYuSEtuM4kwUX2ilcaDuO00BcPdINF9qO4zQQN0R2w4W24zgNxVfaKUZeaFdxdKiambzMuSbH6SU2IuYYFbfO6LNV3Caa3NxEp9jwmGOIjMszEheqiiFyQ4nhMccQGZcB5q4dW14bjbM6Mb/4vj4VlatkS8ohJwN9jhiM17cTk7nG1SPdGHmh7ThOU3GhncKFtuM4DcWFdgoX2o7jNBBXj3RjJIV2v9nWq+i06wr+VKbDjnXRqbq4PC9xAbbaqr8ydOq5UzrtuG6idNopfXVct3ZtZ5vVq3uXZyaU2jOi+cWnU1fihxwddlmfnDY5SRHqdwaqb/eIpIuANwIrzGzfULcdcBmwAHgAONrMHg+fnQEcT5En4xQzuy7U7w9cTGE2ugY41cxM0izga8D+wKPA28zsgVomn2BQjleO4zjjZFPmq5SLgYVR3enADWa2J3BDKCNpb+AYYJ/Q5wuSWt+z5wMnAHuGV2vM44HHzewFwKeBv+vvPPvDhbbjOA2lHqFtZj8EHouqjwAWh/eLgSPb6i81s3Vmdj+wFDhQ0s7APDO72cyMYmV9ZGKsK4BDJKm/c82nsUJb0kWSVkj6WVvdRyT9VtKd4fX6UH+wpIsnbLKO49RMS6edJbTnS1rS9joh4wA7mdlygPB3x1C/C/BgW7tloW6X8D6uH9PHzDYATwLb93O2/dBknfbFwOcpvtHa+bSZfbLqoIMKGJWzT7tMh52j045Vy/EebIDnROV50YG3npfoE9XNi8opnXZcl9rLHdfl6LRTde3E+utUXazTTumr47pYXw3levstEsaLGauicjSXKkuwnD3XOTru+NJWaVO20qtviZltiFxpZgfUdNDU9K1Hfa8+A6GxK+0uP2m6sZ7i281xnClBXyvtKjwcVB6EvytC/TJgt7Z2uwIPhfpdE/Vj+kiaQbFuypVdfdNYod2DkyXdFdQn2wKY2Y/N7NSJnpjjOHWyIfNViauBReH9IuCqtvpjJM2StAeFwfG2oEJ5StJBQV/9rqhPa6y3AjcGvfdAmGxC+3zg94H9gOXAuTmdJF0iaY2kNb7z03EGyyMrV9L6fwuvS/ofpb6VtqR/BG4GXihpmaTjgbOBQyXdBxwaypjZPcDlwL3A94CTzGxjGOpE4EIK4+R/AteG+q8A20taCnyAsBNlUDRZp92BmT3cei/py8B3M/sdCxwLMF0a2Deg4ziww/z5rF69OpU3oz9SG/IrYGZv7/LRIV3anwWclahfAuybqH8aOGo8c+yHSSW0Je3csvgCfwz8rFf7HOoKGFUlC02c4Twn+FNsD4yNjtBpeNxmm97lnDY5xsscB5xZHSncB2OJXFfBcebJVZ1t5kZ1OYbUjmw8T0TlRPCqeCWxMSqnxFdZmyrGy1RdTlad2n/BmtUmtKcajRXa4SfNwRTbeZYBZwIHS9qP4hl/APizCZug4ziDxYV2ksYK7S4/ab4y9Ik4jjN8zOAZT4KQorFC23GcEcd8pZ1iJIV2mXNNmaNMXQGjcpxrYueZjuBPiQPF+ujttutdzmmT0oPHDjjT5yVsT3MjrfzcqM2MRPikSjrtsRreWWvXROXOvOnzVo1tk9LJP1HBOags4FWHjhvYEOm5Y3EV668Bnilpk1qnxnWpK13Faaf2bWiu0+7KSAptx3EmAS60k7jQdhyneRgutLvgQttxnAbi6pFujKTQ7jdgVF2Jfct02ql92nFdnMAgtX+6TKc9f35nn7guLs/ZLqFx32bbaHKJ8FUdUaUinfbMLTr7VNFpr480vJFOOxUNavqqsSl3t33i8Y42s2eP3fAdT61Uf02n7NmYSuIQRarYELVJbO3u0GnH5VSfnOe0LLHvoBITj8F3j3RlJIW24ziTAF9pJ3Gh7ThO8zB8y18Xsn7VSNpd0mvD+zmSUmGcHcdxaiLotHNeI0ap0Jb0XooUOl8KVbsC3x7kpBzHcVxop8lRj5wEHAjcCmBm90nasXeX5lJH5pq6AkbF5Tg4FCQyq0cVcdAmqGaI3DG6o7PmRz+m5ieyJ+V44GwVjROfQOytAmmHm3Y2JFxN4lQ1HWnUxxodAXgi8nKZHYfwgjkzHh1T3nFaYpyIWI7EdtOUHfWZaPrro0BV6xLHietyDJHxoVOmvvj5LjNMdqsbF+5c05Ucob3OzNa38lSGzAwe3tRxnMFhpL/dnCyh/QNJfwXMkXQo8L+A7wx2Wo7jjDaW3hvpZP2qOR14BLibIhTqNcCHBzkpx3FGnJZHpOu0OyhdaZvZJuDL4TUlGIZzTUqnHbuRZGVjjwaKVcJx0CaoptPu0GHvuENU3qmzUzxQTnaFOIDU7IQmXyWPpSV+Nj8dBYSKA0TF+mtIOPp06rRj75m4xfxNnTrusszwcRng6ad7l9ck+sTPS9Ql+QzGw6SudNnzPxSdNuZb/rrQ9b9D0t300F2b2YsGMiPHcRwYyVV0Dr2WNG8Mf08Kf1vJOf8ESCRvchzHqQkzN0R2oavQNrNfA0h6lZm9qu2j0yX9G/CxQU/OcZwRxlfaSXJ2j2wp6dVm9iMASa8Exp9puSFU2addJXECdOpD43LqosaJcWOddipwfxxEKlYrJ4M/xfuwYx12vJEbOvXeqewKc2I9d+xMm9LklzyWSqzA5kQa3TmRrjnWpUPnvuwK0Z/mbIh3R8M2kUJ6TfS7NJVkOK6Ly1smdNpxWof4eepM+1DtWc7RadeO+e6RbuQI7eOBiyS1En8/AbxncFNyHMfBV9pdyNk9cjvwYknzAJnZk4OfluM4I417RHalVGhL+puoDICZuU7bcZzBUeOWP0mnAX9KsSPubuDdFOHqLwMWAA8AR5vZ46H9GRRaho3AKWZ2XajfH7iYIurENcCpZjZUD/Ec9dSattdG4HCKk3QcxxkMrd0jOa8SJO0CnAIcYGb7AtOBYygcB28wsz2BG0IZSXuHz/cBFgJfkNQKiHM+cAKwZ3gtrPO0c8hRj5zbXpb0SeDqgc1owFQJGFWXc01Z5pqZicnMjhqVGSahM4hUhwNOnHEGyj1wYqMjdBort0h47RAbIuMJp8JkJbLZjKHT+NdpdouOMydxnB3LjkPnT/T1UZimpztDOc1bvXxMOY5dtaYziU5Hm/i+J5+NaGo5QctyDObNcK6hbvXIDIpQHM9QrLAfAs4ADg6fLwZuAj4EHAFcambrgPslLQUOlPQAMM/MbgaQ9DXgSODaOidaRpVrPRf4vbon4jiO8yzWVzzt+ZKWtL1OGDuU/Rb4JPAbYDnwpJl9H9jJzJaHNsuB1hapXYAH24ZYFup2Ce/j+qGSo9Nu94ycDuwAfHyQk3Icx+ljpb3SzA7o9qGkbSlWz3tQ7H77pqR39hhPiTrrUT9Ucrb8vbHt/QbgYbNU8AfHcZyaMOrcp/1a4H4zewRA0pXAK4GHJe1sZssl7QysCO2XAbu19d+VQp2yLLyP64dKjtD+WzM7tr1C0iVx3WSlik47xyEhFcY/rov7xHrMVF2Zjhs69dzT50VuO6ms6bEHTlnUKUjosFM67Vh/Hh87lYO+7LFMrRlij5VYh50YM1Zpb5fQlcc665ws7/PGZjDYatXYPql7VnafU8/GjOiU4+cr9QzmONdUef7rx2BjbWvD3wAHSZpLYfw4BFhCsbliEXB2+HtVaH818A1JnwKeR2FwvM3MNkp6StJBFElh3gV8rq5J5pIjtPdpL4QkCPsPZjqO4zjUuk/bzG6VdAVwB8U3/r8DF1BYqy+XdDyFYD8qtL9H0uXAvaH9SWbWSpd0Ipu3/F3LkI2Q0DvK3xlAK/lBa+kgiuiOFwxhbo7jjDI17h4xszOBM6PqdRSr7lT7s4CzEvVLgH1rm1gFuv7SMbP/a2ZbA58ws3nhtbWZbW9mZwxxjo7jjCKeBCFJr5X2Xmb2CwpL60vjz83sjoHOzHGc0cXd2LvSS6f9AQrPn3MTnxnwmoHMaAiUGVLqcC5IuW3EdXE5lZg8rssxRHbUxRHuUh45cdb0jtCAiaw0HY4zCaedjjbPicopQ2TKLaSdVJ7x2FJXktEd6DBozknExdsmqlsVpUlPXcvoes+d278hMr7vqWdji8gQWfZ8QbVnOcfw2N4mtS+uEh7lL0mveNqtDeqHm9mYWJOSUvE0Hcdx6sGTIHQl50v0x5l1juM49WGb8l4jRi+d9nMpXDTnSHoJm3/1zCP9m9ZxHKceXKfdlV467cOA4yi8fj7VVv8UxVbAKUFdzgU5bWIta6xznJ64G2W6zZSuc1asvIqzjsdl6NTNdmR6STjkdAR/SrWJddixjjulaaui0y7rszFRF+uwE/OfG2Vx77hOiWsZ1c2avXJMOcd2EZdTz0bH8xOVqz6nzXCuwYV2F3rptBcDiyW9xcy+NcQ5OY4z6vhKuys5oVm/JekNFJ6Rs9vqPQmC4ziDYwoKbUnTgLea2eVVxyj9pSPpi8DbgPdT6LWPAnavekDHcZxSzOCZDXmvSYSZbQJOHs8YObFHXmlmL5J0l5l9VNK5wJXjOehEM16dXKp/vDc1tVe1TOc4PTHwjBn9ldONoiPPTOzg7dgQHuu0U7rnuE3KPh3XxeOkkiDEecVjcvZgx8d9OtEm4xzj69ChfE5cy/h6R/djxoxOQVN2X1PPRsfzE5VTz2BcV5d+unY9tzElV9qB6yX9JUWqs2c38ZvZYzmdc4R2y1qzVtLzgEcp4tI6juMMCJvK2/neE/6e1FZnZCaXyRHa35W0DfAJiihZBny5nxk6juP0zRRdaZvZuBa9pb9qzOzjZvZE2EGyO7AX8M/jOajjOE5PWuqRKRgwStJRkrYO7z8s6crgC5NFX6ooM1tnZk8C3+xzno7jOH1QXzb2BvJ/zOwpSa+m8IdZDHwxt3OOeiRFbTFheh5EWgh8hsLGcqGZnS1pOwoF/gLgAeBoM3tc0sHAcWZ23DDmVgfTpvUu57TJM0RmWC9j45ly8nbHRrhUm7Ic4SmjY5XHMnaeiY9Tcf7xdSgxMibrMgyRZfc559lg8i06u2M2lQNGtR7WNwDnm9lVkj6S27mq0XfgySwlTQfOAw4H9gbeLmlv4HTgBjPbE7ghlB3HmWpMUfUI8FtJXwKOBq6RNIs+ZHGv2CPfIS2cBWzf7ywrcCCw1Mx+FeZzKUVG5SOAg0ObxcBNwIcofJufHMK8HMcZNFN7y9/RwELgk2b2REgq/L9zO/f6HfrJip/VxS7Ag23lZcDLgZ3MbDlAyKK8Y3j/Yzz6oONMEabulj8zW0ubr0uQZ8tz+/eKPfKD8U1t3KT05pXUMpIuAd7cbVDHcepjxcqVSGrP/HClmR3b90BTd6U9LqoaIofBMmC3tvKuwEPAw5J2DqvsnYEVZQOFB+ZYgC2kgevjHWeU2XH+fH61enUi/GEfGGzc4P+qKYYWZbECPwH2lLSHpJnAMcDV4bUotFkEXDVB83McZ0CYwaaNea/JgqTrJJ0maa/xjNNYoW1mGygCq1wH/By43MzuAc4GDpV0H3BoKDuOM8XYZHmvScQi4HHgI5LukHS+pCMkJRKNdqdUPSLpeuAoM3silLcFLjWzw6rMuh/M7BrgmqjuUeCQQR97GMQqu5QKr6xN0rcgriwrA2yIliwWtVHqQM/EgyTaxAkL4nJO8KeYdRWOU3H+8XWIr1PyWva+3qkuZfc559mYSkzFzSNm9l/AxcDFIUTryym2NH9Q0u+A75vZOWXj5Oi057cEdjjw460dG47jOAPBJpfqo19CiNabw+tvJM2n8I4sJUc9sknS81sFSbszBOcax3FGmzp9ayRtI+kKSb+Q9HNJr5C0naTrJd0X/m7b1v4MSUsl/VLSYW31+0u6O3z2WUm1bEgzs5Vm9vWctjlC+6+BH0m6JGyd+yFwxngm6DiO0wuzQhOV88rkM8D3zGwv4MUUdrKkd3XwvD6GIlvXQuALwUMb4HzgBGDP8FpYywn3QU6Uv+8BL6WI93E5sL+ZXTfoiTmOM7rUGeRP0jzgj4CvAJjZ+qDyPYLCq5rw98jw/ggKu906M7sfWAocGLYYzzOzm83MgK+19RkavdzY9zKzX0h6aah6KPx9vqTnm9kdg5/eYBivfSPVP9YXpfRH8aIgHicVH6eKTbGzUXTk9bEBDlgfGe6ejjKVz0llfomzma9NtImzwZRlTYd05vR2UtnY4/nFc0nNLZ5/4hzj6xBfp9S1LDFWVrBdJp+NjucnKqeewbiuLltf7TbD/nTa8yUtaStfYGYXtJV/D3gE+KqkFwO3A6fSxbuawhv7lrb+y0LdM+F9XD9uJL3bzL6a07aXIfIDFD8Dzk18ZsBrKszNcRynFKOv7XwrzeyAHp/PoNAWvN/MbpX0GXoHmuvmjV2bl3aCjwLjE9pmdkJ4e7iZjVmCSEolDHQcx6kHq3XL3zJgmZndGspXUAjtbt7V3byxl4X3cX0Wku7q9hGwU+44OVv+fkzxLVVW5ziOUwtmsCGheao2lv2XpAclvdDMfknh53FveC2icNBr966+GviGpE8Bz6MwON5mZhslPSXpIOBW4F3A5/qYyk4U2/oej+pFH8Hueum0n0uhr5kTUuG0fhrMI512e1KS+jKP68rKuW1iFV38TG5M6Do7VKglZYB1kWp21to1YyviMsDq1VGbWKf9VGcfYkeuVGb1MueZ1KNUpvdO6bRjnXUcpTc1/9XlbeLr0HGdEtcyqovvR+qeld3X1LPR8fxE5arPaZXnfxDU7O34fuDrISTGr4B3U2zEuFzS8cBvgKMAzOweSZdTCPUNwElm1rq8J1I4yMwBrg2vXL4LbGVmd8YfSLopd5BeK+3DgOMofgKcy2ahvQr4q9wDOI7j9EvdHpFBUKb03knvajM7CzgrUb8E2LfiHI7v8dk7csfppdNeDCyW9JaQ1NdxHGc41KvTnlLkONfsL2mbVkHStpL+doBzchzHmXJR/uoixxB5uJk9qw4JsUdeD3x4cNMaLGVf4HXotFM2lLguLufoOp+O9KNrE9uP47pZZXpZgNWRPveJJ8aW5yb01XPiupzHKf4vS+3/LhsntTk9vhCxfjq2/QBE5/i7JxJNorr4OqWuZXS94/uRumfxfc2xXZQ9T6lnsC77TK8+daii+9zyN1Lk/JdNlzTLzNYBSJpDOoW24zhOLdS5e2SqkSO0/wG4QdJXKb4A38Nm10/HcZz6cZ12V0qFtpmdI+luCiurgI977BHHcQaNq0fSZOWINLN+9yM6juNUxnyl3ZWczDUHUXj9/DcKr4fpwBozmzfguQ2FupwLYtNYyqgd18V9YmNUqi7HEBnbxuatGuvsMX1Vwomkw/AY5WWdnTBj7LjF2PIWnU06zzIO0pRyyEkO1EZK2RmPGxsIE0bGZ1aOLT/2WGeb+LrE5cS13Bhd7w5/nAxDZFkZyp+51DMY92m0c40L7SQ5K+3PU8SW/SbF5vR3AS8Y5KQcxxltrL8ofyNFrnpkqaTpwZXzq5Ky/eQdx3Gq0EeCg5EiR2ivDf76d0o6B1gObFnSx3EcpzK+0u5OjtA+lkKPfTJwGkXIwrcMclKDxBirk6ui00vmHYjKqZBGcV2splyfmEyZDjvl2/HkqrHlraK4Tts+kXA0iXXWM6PytAzn2e0SuuY5sa556/jAiYGqONfEVzPSNaccZ2Id9opHOtusLNF7J67lquj6x/cj6Y8T3dcOZ5vUsxGVy/LRQ+eVS13JOpzL6sB12mlytvz9Orz9HUWgbsdxnIHiHpHd6RWa9W56eKSa2YsGMiPHcRzf8teVXivtNw5tFo7jOG24G3t3eoVm/XW3zxzHcQaNq0fSuHNNRl2OQ0KOgWddSTmRA4W5JYbHuYnEL3Mjw9cTUZvZszs9NebMeHRsRY7hMf79+nR8RsA2kSFybmQQnJ1wrlHJY2mJqxtnTY8jG8ZOMam62OgIsOLhqM3Y6/S7xzqvZTzsUxmGyDIHnNSzUfY85RjMc57liXCuqTsJwlTCnWscx2kevuWvK+5c4zhOI/GVdhp3rnEcp3GYuU67G7nONdOYIs41MHHONbExPNaGpnK4rI0GinWfsxO+KTNn9i7PSNz1HaeNdUbpCA+VWvasj7SoqczksadJ7OkTTw5gRkkG95R/c5zaJb5QccYZ6FQ+pwJGRTrsdSvHjpNSg3f438QxpqJLAgmddpyxqLNLR11HJqREn8niXGPAhtTknL6ca57GnWscxxkGrtPuStctApKOkHRSW/lWSb8Kr7cOZ3qO44wqmyzvNWr0Wml/kGLXSIlhcMIAABUESURBVItZwMso9NlfBa4Y4LwcxxlhPAlCd3oJ7Zlm9mBb+Udm9ijwqKRJa4isEjAqRw8Y6w9TFzbWQcah/uPQStAZyn9mVLFFSiU8o3c5Zwv2/E1jdbdzUu5p8b7s1AbkWIcdJ1eYmUh4kFK6t5NSdq6P5hfr11NzixMYJII/xfuwYx12FZ12HEAqNb14tqlno8wmkhO0rMpe7qEFjKpRPSJpOrAE+K2ZvVHSdsBlwALgAeBoM3s8tD0DOJ4ij8QprfSKkvYHLqbI3HENcKqZDX2t3+vfd9v2gpmd3FbcYTDTcRzHCW7sG/NemZwK/LytfDpwg5ntCdwQykjam0LDsA+wEPhCEPgA5wMnAHuG18JxnmYlegntWyW9N66U9GfAbYObkuM4o07LIzLnVYakXYE3ABe2VR8BLA7vFwNHttVfambrzOx+YClwoKSdgXlmdnNYXX+trc9Q6fU79DTg25LeAdwR6van0G1PyGQdxxkR+tNpz5e0pK18gZld0Fb+ewobXXsw953MbDmAmS2XtGOo3wW4pa3dslD3THgf1w+dXgGjVgCvlPQaip8KAP9sZjcOZWaO44w0fei0V5rZAakPJL0RWGFmt0s6OGMsJeqsR/3QydmnfSMwpQT1MJxrcoyVZYZJKCJ0tTMjmsyMhFFreqT0qhL7Kbb1bZNIBz5v9fKxx52XmMzcKCBUbIhMOdJUMUTGys3YEBkHkKIza3rK6aXM/yblj1NmiIwDSAGsiq5/nLA9kcC91PBY9TlthHNNfbtHXgW8SdLrKdIkzZP0D8DDknYOq+ydgRWh/TIKB8IWuwIPhfpdE/VDJ+Pf2XEcZ/jUsU/bzM4ws13NbAGFgfFGM3sncDWwKDRbBFwV3l8NHCNplqQ9KAyOtwVVylOSDpIkisB5VzEBZAWMchzHGSZDSIJwNnC5pOOB3wBHFce1eyRdDtxL8UPkpBAoD+BENm/5uza8ho4LbcdxGkndzjVmdhNwU3j/KHBIl3ZnAWcl6pcA+9Y7q/4ZSaFdptOO9XyxDimlU8oJGBVf7Fibm9Jpx2GO4mPPSB0oEe+/ndQ/Q6wmjuMvrUkoVWOHkK1WdQaMmjt3bN2s2ZE3Skp/XUmnPbZuXUkWe+icf8rpJdY/x/rpnNwKHQGjEvcsdqaJ73tOwKi6srGXtRlWwKhRdFHPYSSFtuM4Dcfd2LsyNEOkCj4raamkuyS9tO2zYyTdIenP2+rOkvSgpNXROMdJekTSneH1p6F+gaSbhnU+juMMDqPY8pfzGjWGuXvkcDa7f55A4RLa4hiKYFQHSWoFq/gOcGCXsS4zs/3C68IubRzHmaxYfR6RU41hqkeOAL4WXEBvkbRNa58kmzeuP7uJ3cxuASh212SxEUjsmu1NXfu0Y/1h6tuwTDeeCv1f1iZ5deLJRDrVjRV02imd8JqcJMNRXWdChs6rmaprJ0OlnTX/uC4n4W68lztHDx7rsJ/s7FIaICq1TztO5FtXwKgq+7TrxsyTIHRjmCvtXYD2qIHtbqBXUkTgWmJmiRQjHbwlqFiukLQbgJk9aGZvrnXGjuNMGL7STjNMod3VDdTMFpvZS8zs3IxxvgMsMLMXAf/C5qAv3Q8sXSJpjaQ1I3iPHWeoPLJyJa3/t/C6pN8x6gwYNdUYqNCWdFLLYEjh8plyD+0LM3vUzFq/DL9MEcSqrM+xZralmW3pLqCOM1h2mD+f1v9beB1bZZxNma9RY6AyzMzOaxkMgW8D7wq7SA4CnmxF2eqHECegxZsYGyPXcZwpggvtNMM0RF4DvJ4iPu1a4N29Gks6B3gHMFfSMuBCM/sIcIqkN1HYUB4DjutnEnHmmhRlhscc55oqxsoc42UOsU/ChujAGxLm2mfi7N8VnFNSmeHjupzM8GUBrqo4ByXiXVU6x7IydAZ/ipskunTUxVOp4lxTJStNTpsy430dPjGWmIdTMDShHXaNnFTacHP7D1LEwI3rzwDOqHFqjuM0jJzF1ajiHpGO4zQSF9ppXGg7jtNIXGinGUmh3f4wpNSnVZxrYlKODTE5Ou0qxJ69HfNP/DesjxxCcvS9sb46pdOOddgTpdOOy9B5jim9d3zea6NxUvrp+FLFjgepzOplSQ9SOu2yAFFVA0ZVcaYZSMComsecKoyk0HYcp/m40E7jQttxnMbhu0e640LbcZxG4ivtNCMvtKvo63JWADl7rnP03mXkBLyKddyp48bBh9ZEjbZMdJoZnVCOTnt69MTFSYihmk47DoK1sUTHDZ067PWJceO0DrFuOUc/HffJSdJbVobyfdnDTOxbN67T7s7IC23HcZqJC+00LrQdx2kcvtLujgttx3EaiQvtNC60HcdpHL57pDsjJ7Sr/OyKH57URatirKxCjlEoNjw+U1KGTkNkbFNMGdxmRwefkbCwbRHVxRnnkwbbCobIuCrnnON7lDL2lWWHyQnkFJfjMVNtchxlytpUdZwZb8CouvCVdpqRE9qO4zQf12l3x4W24ziNxIV2Gk/k4jhOI6krCYKk3ST9P0k/l3SPpFND/XaSrpd0X/i7bVufMyQtlfRLSYe11e8v6e7w2WfVR+bxuhj5lXaVb/McPWBdVHF0KHOuSel347pYxzor0Sfym0lmk4912HGbVJ8qNyU+x5xzruJ0VHadUuPE5Sp9cpJqVHGuyXl+Jsq5psb/qQ3AX5jZHZK2Bm6XdD1FApUbzOxsSacDpwMfkrQ3cAywD/A84F8k/YGZbQTOB04AbqFI7LIQuLa+qZbjK23HcRpHS6ddx0rbzJab2R3h/VMUKQp3AY5gc2LwxcCR4f0RwKVmts7M7qfItnVgSHU4z8xuDkldvtbWZ2iM/ErbcZxm0seKfr6kJW3lC8zsglRDSQuAlwC3Aju18tSa2XJJO4Zmu1CspFssC3XPhPdx/VBxoe04TuPoc/fISjM7oKyRpK2AbwF/bmareqijUx9Yj/qhMpJCuywJQq/2qXKKuvRx/SYhTtXFc8nZ8xvrq1P7tKvotOPrXZcVJ/7PKdu3DXk67bJrN6w+VYI/Vd2nXZedZ7zUOaakLSgE9tfN7MpQ/bCkncMqe2dgRahfBuzW1n1X4KFQv2uifqi4TttxnEZS4+4RAV8Bfm5mn2r76GpgUXi/CLiqrf4YSbMk7QHsCdwWVClPSToojPmutj5DYyRX2o7jNJuad4+8CjgWuFvSnaHur4CzgcslHQ/8BjgKwMzukXQ5cG+Yxklh5wjAicDFwByKXSND3TkCLrQdx2kgdXpEmtmP6K6JO6RLn7OAsxL1S4B9a5paJVxoO47TSNwjMs3IC+2cByNW/Fd9mMp+7lVxdKgSvCrVp8wQluoT16XalGWcH5YhMschqoqxr0rG8yp9coyKdTjOpJgI55phHmeyMfJC23Gc5uEBo7rjQttxnEbiQjuNC23HcRqHJ0HozkgK7TLnmmF9w8cPZZVN8zl6yxzdc5neOzW3nDZlOu26HAXqCKxVRW88KKeXQWVNr6tNWZ868JV2mpEU2o7jNBvXaXfHhbbjOI3EhXYaF9qO4zQOX2l3x4W24ziNxA2RaUZeaKe+zetwpskZNy7nGAirONfkGNzKjIpVHGdy2wyCugxuZde/LqeXKhnQq/SJqRLlb1iZa3ylnWbkhbbjOM3EhXYaF9qO4zQOX2l3x4W24ziNxIV2GhfaCeKHpYqOu4rTTko/WqYTrqKTT82tzNGnqr56suu0B+G0M6jjDCr400QJTxfaaVxoO47TONyNvTsutB3HaRyu0+6OC23HcRqJC+00LrQzyNEJT1Rg+Bw9cpU+g9JXT5ROO6fNsHTNde2NHpR+uinCsinzaBoutB3HaRyuHumOC23HcRqJC+00LrQdx2kcvnukOy60HcdpJL7STtNYoS1pNvBDYBbFPK8wszMlbQdcBiwAHgCONrPHJR0MHGdmxw16blWDTFVx0skxIvbbZ5hGxbI2g8pc0+/nVds02ellMhkdY1yn3Z1hGfOrsA54jZm9GNgPWCjpIOB04AYz2xO4IZQdx5libMp8jRqNFdpWsDoUtwgvA44AFof6xcCR4f164MmhTtJxnIHQWmm70O6kseoRAEnTgduBFwDnmdmtknYys+UAZrZc0o7h/Y+BH5eNOXfuXPbaa6/a55rz7acaxojblI2Z26eKeqTKsauMkYOVfJ7zz102Rmqc1LjxOFWOXYdqo8r51MGTT9azdnJDZJpGC20z2wjsJ2kb4J8k7VtlHEmXAG8GmDZtGlJdomIzOf8gdQiWbtx///3sscce4xihuUzVc2vCedX/n1Ccl6Q1bVVXmtmxfQ5zHTA/s+3KPsee1MgsR9xMPJLOBNYA7wUODqvsnYGbzOyFfYyzxsy2HNQ8J4qpel4wdc/Nz8upQmN12pJ2CCtsJM0BXgv8ArgaWBSaLQKumpgZOo7jDJ8mq0d2BhYHvfY04HIz+66km4HLJR0P/AY4aiIn6TiOM0waK7TN7C7gJYn6R4FDxjH0lePo22Sm6nnB1D03Py+nbyaNTttxHMdpsE7bcRzH6cSFtuM4ziRi0gttSQsl/VLSUkmnh7rnSbpR0lWStgp175N0t6Q7Jf1I0t5tYyySdF94LWqrv0nSgmGfUy8kXSRphaSftdV9RNJvw7ndKen1of5gSRdP2GT7oMt93E7S9eG+XC9p21A/Yeclabak2yT9VNI9kj4a6rvdg5mSvhqevZ+GGDmtsQ6WtETSOYnjfE7S6qjtk23j/03bZw8M8pxT9LgOjbtnU41JLbTDzpLzgMOBvYG3B2F8CvB+4ELgnaH5N8zsD81sP+Ac4FNhjO2AM4GXAwcCZ7YetIZyMbAwUf9pM9svvK4Z8pzGRY/72MQ4M91i4kD6HrwXwMz+EDgUOFdS6//uROC/A9MlPeumK+kAYJvEsf+1bfyP1X9qfeGxgSaISS20KYTsUjP7lZmtBy6liE0ync2hCQRgZqva+m3JZgfFw4DrzewxM3scuJ7NQvExYOPAz6IPzOyHFPPKYbLEY+l2HxsXZ6ZHTJxu7E0hvDCzFcATwAHhs2lsDrMhePYL7BPAB/uY1iN9tK0Fjw00cUx2ob0L8GBbeVmo+zzwJeB9wD+0PpR0kqT/pFhpn1IyBmb2ZjNr/6zJnCzprqA+2RaKeCxmdupETyyDbvdgTJwZ4Nk4MxN5XpKmS7oTWEHxhX9r+KjjHgA/BY6QNEPSHsD+wG7hswsp4uVMM7Oft8YArm6dd8QrgjriWkn7tCrN7GU1n2IWXa5DI+/ZVGKyC+1U6AQzs1+b2R+Z2f80s6faPjjPzH4f+BDw4V5jDGCug+R84PcpfqYuB86d2On0zaS6B2a2MajZdgUODDFxut2Diyi+hJYAf08hpDeEca4zs5ea2V9AYYuhcBb7XOKwdwC7B3XE54BvD+j0sulyHZwBM9mF9jI2r1qgeHgeyuh3KZt/tlUdozGY2cPhH2gT8GUKdcNkots9eDjElyH8XTEBc+uKmT0B3AQs7HYPzGyDmZ0W9NBHUOiq7+sy5EsoIlouDcbFuZKWhnFWtdQRQV++haTcgEoDpf060PB7NhWY7EL7J8CekvaQNBM4hiI2SQeS9mwrvoHN/zjXAa+TtG34Sfu6UDdpaP2TBP4Y+Fm3tg2l231sXJwZdYmJ0+0eSJoracvw/lBgg5ndmxrbzP7ZzJ5rZgvMbAGw1sxeEPo+V1JL730gxf/uowM5yQy6XQcaeM+mGo11Y8/BzDZIOplCyE4HLjKze7o0P1nSa4FngMcJD5aZPSbp4xSCA+BjZpZr6Bs6kv4ROBiYL2kZxc6XgyXtR6FSeAD4swmbYAW63UdJZ9O8ODPdYuJc0uUe7AhcJ2kT8Fug3xClLd4KnChpA/A74BibWHdmjw00Qbgbu+M4ziRisqtHHMdxRgoX2o7jOJMIF9qO4ziTCBfajuM4kwgX2o7jOJMIF9qO4ziTCBfaI4SknSR9Q9KvJN0u6WZJf1zSZ4HawsD2ebzjgmt2q3yh2kLilvQ9WNJ3qxw3F0k/Dn8XSHpHhf7HSfp8/TNznO640B4Rgjfdt4Efmtnvmdn+FJ6Huw7wsMcBzwptM/vTbt6AE4GZvTK8XQD0LbQdZyJwoT06vAZYb2ZfbFWEwFqfg2dXm/8q6Y7wemU8QK82kj6ozYH+z5b0VooQpF9XEbR/joqkEgeE9gvDGD+VdEPuSUg6RNK/h2NdJGlWqH9A0kfDmHe34lMHd+vrQ/2XJP26FbNDm5MMnA389zDP0+IVtKTvKiQvkPRuSf8h6QfAq9ra7CDpW5J+El7PfuY4deJCe3TYhyJSXDdWAIea2UuBtwGfzW0j6XCKAFwvD1HozjGzKygi2/1JCJb0u9YgknagCKr0ltA+y9VZ0myKJBBvC0kFZlAkEmixMsztfOAvQ92ZwI2h/p+A5yeGPp3NCQY+3eP4OwMfpRDWh1LEym7xGYokCC8D3kIRdtVxamdSxx5xqiPpPODVFKvvl1EEsf98iJ+xEfiDRLdubV4LfNXM1kIRz6Xk8AdRqGnuz2zf4oXA/Wb2H6G8GDiJIuQpwJXh7+3Am8P7V1MEcMLMvifp8cxjpXg5cJOZPQIg6TLGXoO9Q0wngHmStm4PDew4deBCe3S4h2IFCICZnRTUBEtC1WnAw8CLKX6BPZ0Yo1sb0V/8637bt/frxbrwdyObn+2yPik2MPZX6Oy2993mPQ14RfsvCscZBK4eGR1uBGZLalcnzG17/xxgeYgHfSxFtL2Ybm2+D7xH0lx4Nu8mwFPA1olxbgb+h4pMLu3ty/gFsEDSC0L5WOAHJX1+BBwdjvM6IJX/M57nA8B+kqZJ2o3N8clvpYiouL2kLRir1vk+RdYZwrH2yzojx+kTF9ojQgjjeSSFsLxf0m0U6oUPhSZfABZJuoXiJ/+axDDJNmb2PYo4yktUpJ9q6ZMvBr7YMkS2zeUR4ATgSkk/BS7rMu1DJC1rvSiSBLwb+KakuylyK36xS98WH6WIl34HReLg5RRCup27gA3BKHoa8G/A/cDdwCcJtoCQPusjFF86/8JYG8EpwAEq0o3dS5HqznFqx0OzOlOasLtkY4jZ/Qrg/JAiy3EmJa7TdqY6z6cIyj+NIiP4eyd4Po4zLnyl7TiOM4lwnbbjOM4kwoW24zjOJMKFtuM4ziTChbbjOM4kwoW24zjOJOL/AwgDlr/yCeV9AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "m = Map.create(skydir=(0, 0), width=(1, 1), binsz=0.02, frame=\"galactic\")\n", "m.quantity = gauss.evaluate_geom(m.geom)\n", "m.plot(add_cbar=True);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again for convenience the model can be plotted directly:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAEMCAYAAAAPqefdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2dfbxdVXnnv7+8J4QYIIBIFLCmZYDxDUSqHYePVAlqC1WxqS1EpTI68WXsdCq0TtVaZtD6Uq0I0ooE1CIiUzIdESmUtlbeAlURX0oqKJEUDAkQCAm5yTN/7LXuWXdl7332ufecc8895/l+Pvuz9l577b3X2ufc5z7nWc/zLJkZjuM4zsxg1nR3wHEcx2mOC23HcZwZhAttx3GcGYQLbcdxnBmEC23HcZwZhAttx3GcGURfhbakSyQ9JOl7Sd0zJN0o6RpJi0PdfElflrRB0q2SDk/ar5Z0T9hWJ/U3pe0cx3GGkX5r2pcCK7O6dwHvBP4K+J1Qdxaw1cyeA3wC+DCApP2B9wMvBo4H3i9pv046IOnyyXZ+kBnWccHwjs3H5UyGvgptM/tHYEtWPRvYEzaFulOBtWH/KuAkSQJOBq43sy1mthW4ntY/gS3A7gbdeO3kRzDQDOu4YHjH5uNyOmbOdHcA+DRwOfAo8MZQdyhwP4CZjUl6FDggrQ9sDHWYmX9RHMcZeqZdaJvZT4CXZdUqa1pTX0v4uRaF+qIDDjhg6GL3582bxzCOC4Z3bMM6LklIeiKputrMzujkHitXrrTNmzc3anvHHXdcZ2a52XVomXahXcFG4JnARklzgKdRmD82Aicm7ZYDN7W7WfjCnAEgybZsyS00w8GwjguGd2zDOK7DDjuM++67b5+p3GPz5s2sX7++UVtJy6byrJnGoLr8rQOiZ8jrgRutyGx1HfBKSfuFCchXhjrHcYaOPQ230aKvmrakv6bQlJdJ2gi838w+V9L0c8DlkjZQaNirAMxsi6QPAbeHdn9iZsOnqjjOyGPA2HR3YiDpq9A2s99q2G4HcHrFuUuAS7rZL8dxBpHR06KbMKg2bcdxRhrDhXY5LrQdxxlQXGiX4ULbcZwBxDXtKlxoO44zoLjQLsOFtuM4A4h7j1ThQttxnAHFNe0yXGg7jjOAuE27ChfajuMMKC60y3Ch7TjOAOKadhWDmnvEcZyRZ6zhVk/Fill/JumHkr4r6f9IWpqcOzesmvUjSScn9cdKuiuc+1TI8V+70lYvcKHtOM4AEjXtriSMupS9V8y6HjjGzJ4L/CtwLoCkoyhyHR0drvmMpNnhmguBs4EVYYv3LF1pq1e40HYcZ0DpjtAuWzHLzL5hZlFNv4UizTMUq2ZdYWY7zexeYANwvKRDgCVmdnPIOHoZcFpyTdlKWz3BhbbjOANIVzXtdrwFuDbsV62OdWjYz+snXBP+EcSVtnqCT0Q6jjOgNBbIyySlKyZcbGYXN7lQ0h9RGMa/GKtKmrVbNWtSK2pNFhfajuMMKI2F9mYzO67Tu0taDbwGOCmYPKC1alZkOfBAqF9eUp9ek6+01RPcPOI4zgASw9in7j1ShqSVwHuBXzez7cmpdcCq4BFyBMWE421mtgnYJumEYK8+E7gmuaZspa2e4Jq24zgDSPf8tMtWzKLwFpkPXB/mDG8xs7eZ2d2SrgS+T/EfYY2Z7Q63ejuFJ8pCCht4tIOXrrTVK1xoO44zoHRHaFesmFW2zGFsfx5wXkn9euCYkvrKlbZ6gQttx3EGFI+ILMOFtuM4A4iHsVfhQttxnAHE82lX4ULbcZwBxTXtMlxoO44zoLjQLsOFtuM4A4jbtKtwoe04zoDiQrsMF9qO4wwgrmlX4ULbcZwBxb1HynCh7TjOAOKadhUutB3HGVBcaJfhQttxnAHENe0qXGg7jjOguNAuw4W24zgDigvtMlxoO31hulbb8D/7mYrnHqnChbbjOAOI27SrcKHtTJlBXrOuSd9cNAwq/smU4ULbcZwBxYV2GS60nY6ZimbdL628kz/3vE8uKgYBN49U4ULbcZwBxCciq3Ch7TjOgOKadhkutJ1aOjFnTMb00S1zSf7nPZUJyLJrXXz0GzePVOFC23GcAcWFdhkutJ0JNNFQ27XptXY+VaIoqHp2majwycrpwN9yGS60HccZQNw8UsUgx0U4fWBWtrU7P6vNuXbbnGybzD16dd8mY2763pypEr1Hmmz1SLpE0kOSvpfU7S/pekn3hHK/5Ny5kjZI+pGkk5P6YyXdFc59SpJC/XxJXw71t0o6vAsvoBL/vjmOM6Dsabi15VJgZVZ3DnCDma0AbgjHSDoKWAUcHa75jKTZ4ZoLgbOBFWGL9zwL2GpmzwE+AXy4s3F2hgvtEaSJhtiJ9tmJJtxOQ57s1vR5k9Wwm74v/4PqJt0R2mb2j8CWrPpUYG3YXwucltRfYWY7zexeYANwvKRDgCVmdrOZGXBZdk2811XASVEL7wV9/Y5JWhl+cmyQFP+zPUPSjZKukbQ41L0t/Az5tqRvhv9+8R6rw0+aeyStTupv6vXPEsdx+kW0aTcS2sskrU+2sxs84GAz2wQQyoNC/aHA/Um7jaHu0LCf10+4xszGgEeBAzoZbSf0bSIy/MS4AHgFxYBvl7QOOBN4J/Bs4HeAi4AvmdlF4bpfBz4OrJS0P/B+4DiKT/UOSevMbGu/xjGTaacpNq2fVXGuk7ZN+tQJVfrWnpLzs7K6/LhJ27rn5tc4k6XxG9xsZsd16aFlGrLV1Ndd0xP6qWkfD2wwsx+b2VPAFRQ/K2bT+pcpADN7LLluH1ov4GTgejPbEgT19bTsSluA3T0fheM4faAjTXsyPBhMHoTyoVC/EXhm0m458ECoX15SP+EaSXOAp7G3OaZr9FNoV/3s+DTwWeBtwBfiSUlrJP0b8BHgXW3ugZm91szSc47jzGi64z1SwTogmldXA9ck9auCR8gRFBOOtwUTyjZJJwR79ZnZNfFerwduDHbvntBPP+3SnxBm9hPgZSUnLgAukPRG4H0UL2VSP0MkXQ68trPuznya/EduZ+KoM2/k5yZjSum2eSTXu6rq687NKmlTdW2dKaSJKWUY2bx5M5KeSKquNrMzOrtL9/y0Jf01cCKF7XsjhYn1fOBKSWcBPwVOBzCzuyVdCXyf4j/CGjOLv+DfTuGJshC4NmwAnwMul7SBQsNe1ZWOV9BPoV31s6MdV1C42sR7nJjd46Z2NwhfmDMAJPXsP6DjOLBs2TIef/zxfaZ8oz3dEdpm9lsVp06qaH8ecF5J/XrgmJL6HQSh3w/6KbRvB1aEnxw/o/hv9MayhpJWmNk94fDVQNy/DvhfiSP8K4Fze9flmUk77XWqGnGVht2kbbv6qrqUTrTmsvq8rqpsQtkEZbtEVKOicU8Js64J7WGjb0LbzMYkvYNC8M4GLjGzuyuav0PSrwK7gK0Ee5GZbZH0IYp/AAB/YmY9M/g7jjONuNAupa+5R8zsa8DXGrR7d825S4BLutmvYaATd75O7NOdaNpzOmjba027TnuO+2MdtJ2MFp7fI+IadwPMYJcvglCGJ4xyHGcwMf+3VoYL7SGmmxp2+kWp0qzLNO28rhNNvorJ2KdTna3K8yPXvPPrmvapXfCO0wC3aVfiQttxnMHEhXYpLrRnOFP1xW5nly7TtOdUlE3atntOWX9z6jTtXFsey8p0f052PKukbbvn1PWvnYbdxONkZDFcaFfgQttxnAHEzSNVuNCeoTTxFunE57pKE87Lsrp5DdrGNlX3Lwt1nZ0dlyWWiZFSuUYd/9yfyurL2sY2de8nvybX0svwJFNTwL1HKnGh7TjOYOKadikutB3HGTwMd/mroJHQlnQYsMLM/k7SQmCOmW3rbdecMiaTl7rOva7dRGFu+iirm1dRpm3mZuXs7Djtk7I2kWgeSZPHxD/rXVmbeBz7kv7QfiorcxPIU0nbqveai5MmZpL8nm4mqcNt2lW0dT6Q9FaKJXQ+G6qWA3/Ty045juOwZ0+zbcRoommvoVjA4FYAM7tH0kH1lzj9YjIBM03c96q06LK6WC4oaRs16fkV15RNREYNu2oCL52QzCcicy16Zyh30SKe25H1Ide80zqyc3Wadbvps06STI0sHlxTSROhvdPMnorrVIaVGTy9qeM4vcOAMfceKaOJ0P4HSX8ILJT0CuC/Av+3t91ycjpJtzqZhE7tNOwFtIj7c7PjvCyryzXu3LYNMKfNYMcSBSy3Zeca9o6sTPdnZ2XZ+2lHmdbcxB0wvbYMt20b7B7d0dfR5Pt5DvBz4C7gv1Bk6XtfLzvlOM6IEyMi3aa9F201bTPbA/xl2JwBoZPw9U5s2lX26UXJNbn2HM8tzOqhWJU5rZsfHjQ/3HheKGcnnZob9mdlg4x/n2nMxe6w/1RQsXfGMtRHrTpd+yrWPRmfl9V3ki62TGTkYfF1bfP7j54IqsLc5a+CSqEt6S5qbNdm9tye9MhxHAdGUotuQp2m/ZpQrgnl5aH8bWB7z3rkTKBdetXJLAdW5nNdVS7KSmhp1Ptk5/bJSoBF4UaLFk0sFwTVO2rac5JOzWmjaafzU2OZpr0jqMvbt08sFyRuIFHrjhp2fHSV18qEPtScy9s0SY1bdc3Ip3U184nICiqFdlglHUkvNbOXJqfOkfTPwJ/0unOO44wwrmmX0sR7ZB9Jv2Jm3wSQ9BImKlPONDCZdKtNEjrldurcXg2wbyj3yY4XhwcuXtxqG/dj2U7jhs407XYa9uOPTywB5oX9ueF+TTTsKuqWKIs0sXGPrEZdhbn3SBVNhPZZwCWSnhaOHwHe0rsuOY7j4Jp2BU28R+4AnidpCSAze7T33XIcZ6TxiMhK2gptSX+cHQNgZm7T7hGduJx1EsZeZh6pCqKJZTSL7EuLaBaJP72WhIv3XRLql7TaLlkysYxmkoXBPLKwxDwS96vMI08lk4px/8lgHnkyM4s89lgYT+KHODdOfoZzs/JY9RIms9p7PplYNcmYX99p26Gliy5/kt4D/C6FR9xdwJsprH9fBg4H7gPeYGZbQ/tzKawMu4F3mdl1of5Y4FKKP42vAe82s75GiDcx4z2RbLuBUygG6TiO0xui90iTrQ2SDgXeBRxnZsdQTGOsoggcvMHMVgA3hGMkHRXOHw2sBD4jKU59XAicDawI28puDrsJTcwjH0uPJX0UWNezHjm1TMblr5OJyDw0vcyNL2rYS0OjpUvLy3Q/17Tzicn5aUROnImck3094x9o8oe6M5uAzCce8wlPaGny40E8j4QTaax7IE9SVZXyNe1plRZetaJQWZuR0abr6K55ZA5FKo5dFBr2A8C5wInh/FrgJuC9wKnAFWa2E7hX0gbgeEn3AUvM7GYASZcBpwHXdrOj7ZjMhPki4Nnd7ojjOM440abdLIx9maT1yXb2xFvZz4CPAj8FNgGPmtk3gIPNbFNoswmI2UsPBe5PbrEx1B0a9vP6vtLEpp1GRs4GDgQ+1MtOjSqd/AftxOWvamGDdL+dy19q04427KhF779/eZm2yUstDrp7LBckToXzgr5fadRuJVqdv6MIRp//eBEyszSUjwTtuSOXwi2hTGzcUcOOGnWeoCpfkzLdz139qmzbVXVVjIw23lzT3mxmx1WdlLQfhfZ8BIX321ck/U7N/cqWLLWa+r7SxOXvNcn+GPCgmXmokuM4vcPopp/2rwL3mtnPASRdDbwEeFDSIWa2SdIhwEOh/Ubgmcn1yynMKRvDfl7fV5oI7T81szPSCkmX53VOb2iakrVJutW6MPYqr5Foy16cPCB6ieSa9rJlE4/TuoX7RwP4fkW5JOju48btxGoeNe0qm3aiabM9BKUHI7YeK1bB22/R1mI8C3bsdat2iajGtrTO7YrnwnFM+Zpr2qkWk7/nJrbtKvE0Mlr1XlgrG9jU+SlwgqRFFHnCTgLWUzhXrAbOD+U1of064EuSPg48g2LC8TYz2y1pm6QTKBaFORP4i251silNhPbR6UFYBOHY3nTHcRyHrvppm9mtkq4C7qT4//ovwMXAYuBKSWdRCPbTQ/u7JV0JfD+0X2Nm8f/022m5/F1LnychoT7L37lAXPzgsVhNkWf+4j70zUmoW5y3advJ+GmPa9pJaPrTKjTtWB6ULEY3f1nQqJcdMLHRuFtJ8EVZlKSkGs8ilS3tOxb+blJH7eg28liI+Ro3ZhfLLSyc83DRp1l7r0Odh8XH2+5K/cDDNz86lsT3knuRpK7e+SIInXxme/Wx5pqhp4veI2b2fuD9WfVOCq27rP15wHkl9euBY7rWsUlQ+T0ys/9tZvsCf2ZmS8K2r5kdYGbn9rGPjuOMIr4IQil1mvaRZvZDipnWF+bnzezOnvbMcZzRxcPYK6mzaf8eReTPx0rOGfDynvTI6YhuBddUBdXEfNipeWRJZh6JZZx0HDeJADz94HDywImNopkk3mxR4vIX3f+UfT2j09KOJ1t128P+Y/E+oefzwmqUYQZyfusKlu0pTCVVubh3JEE243WxDRPLunfazkwykqHpneBZ/kqpy6cdHdRPMbMJsWKSFpRc4jiO0x18EYRKmniPfAvIzSNldc4kmUxYap3W1klwTdSwoya61/qPYX6wTNPONe5xt7446QgtDfvpT594LmraC2PMexq+E5+efT0V/ogXJjrEwjDBGDX1du6CwMKxYhpxaVCjq3JxT6h7amLP4vuKPUnfaZyUbLeCUNln1oluOfTugL5GZCl1Nu2nU4RoLpT0AlrRQEuYuPqU4zhOd3GbdiV1mvbJwJsoon4+ntRvo3AFdHrEZFzD6trWaXhR046aYtQg46rp+UozsHfSp/EEUTFwpiy6JmrYBwUb99xl8aJ416RX0b49l4lER7vEph2vWxiuOSi7ZtxwvbNVt6PYX7p9E9BKLpUnsUr354c288Pt8gRbTX7pNKFKex5J+7cL7VLqbNprgbWSXmdmX+1jnxzHGXVc066kSWrWr0p6NUVk5IKk3hdBGCA6sWmnH3oMX8mDbOaHnZhwKdU+F2Z27vHkTzE0Pc3NmkfejGvYsQza+QSbdnxY/vWMdunE6DyulYe2UfXdP2jlMeQ9hrtDEvL+WBjHExPGlY41jj++j3mZph3fX5n3SJNfOk4NQyi0Jc0CXm9mV072Hm2/P5IuAn4TeCeFXft04LDJPtBxHKctZkUymCbbDMLM9gDvmMo9mniPvMTMnivpu2b2QUkfA66eykOd5lRpZ5PxHilLGDU3qxvXvOdNLNNFBBbm2vd4mtWgesfQdGi5mIx7icQyatjR/l2maac+GdDyy0g9TnO7d/QwCXbvJbFM+rT4kQn9XrQoaNrhtmULJsRydlDyc//29J029Rqp8x6pW7l9JDCGUtMOXC/p9ymWOhv/CWhmW6ovadFEaMdZn+2SngE8TJGX1nEcp0fYMLv8vSWUa5I6o+HiMk2E9t9KWgr8GUWWLAP+spMeOtNLrtmp5Fysi5rj7GzFr7KFd8eXCIsRjDEacYL7Rb408OLsOJaJHXxck67StPN62NuzJNx30SMlfZq48ML8bKGEsgUTZmcm8/i+6t6p266nyJBq2mY2JaW37ffKzD5kZo8ED5LDgCOB/zeVhzqO49QSzSNDmDBK0umS9g3775N0dYiFaURHyoCZ7TSzR4GvdNhPx3GcDujeauwDyP80s22SfoUiHmYtcFHTi5uYR8ooWyut60haCXySYn7sr8zsfEn7UxjwDwfuA95gZlslnQi8ycze1I++DTKduJrFD3LcdS00mpuZR9Ko8HHzwbjtZG5WJvaF8bUf8zVx8lUo08nF2CZN85T2MiVeHwPKs+fF50+w72T9DeOYF/z50rHG/fH3Ed7P7KDglf0huKtfFzAb5oRRcUGFVwMXmtk1kj7Q9OLJfo96vpilpNnABcApwFHAb0k6CjgHuMHMVgA3hGPHcYaNITWPAD+T9FngDcDXJM2nA1lcl3vk/1IunAUcUFLfbY4HNpjZj0N/rqBYUflU4MTQZi1wE/BeilmqR/vQrxnP7JL9vVwKZ5WXE/bzZc3H1fLkCcpDenInuTysB1oadtXXc3eyP672l99fJX2q6PesWWMTqtP9fF3JeDg7K50uMdwuf28AVgIfNbNHwqLC/6PpxXXmkY9O8ly3OBS4PzneCLwYONjMNgGEVZQPCvvfosg+6DjOjGd4Xf7MbDtJrEuQZ5uaXl+Xe+Qfpta1KVNmLpyUWUbS5cBrp9Ydx3GasHnzZiQleQO42szO6PhGw6tpT4nJTkT2g43AM5Pj5cADwIOSDgla9iHAQ+1uFL4wZwBI6rk93nFGmWXLlvH444/v075lDQa7x/xPtYxBFtq3AyskHQH8DFgFvJEi7nk1cH4or5m2Hs5Qdpfs5zpNVHLycsJ+dLfKlzUfS54QlwiLCxiMB8HE46eyEqotxDG9ato27lfc30r6VNHvurHmSl883J2VTncwgz1D9lIlXQd8Hbg2rL87KQbWC8nMxigSq1wH/AC40szuphDWr5B0D/CKcOw4zpCxx5ptM4jVwFbgA5LulHShpFMlLW53YUpbTVvS9cDpZvZION4PuMLMTp5MrzvBzL4GfC2rexg4qdfPnsnsaVOmxO98VGrGQqOYPG0sK6G1EG5rgYEsDepTiSYcF+EdXyJsPJVNKKtC1tNejT853jSp256VT05sE5+f9invbxhHbJKONe6Pv489E3tWJjM6ef9OOcPoPGJm/w5cClwaUrS+mMKl+Q8kPQl8w8w+0u4+Tcwjy6LADg/eGj02HMdxesIQmkdSQorWm8P2x5KWUURHtqWJeWSPpGfFA0mH0YfgGsdxRptuxtZIWirpKkk/lPQDSb8saX9J10u6J5T7Je3PlbRB0o8knZzUHyvprnDuU5K6Eh1uZpvN7ItN2jbRtP8I+Kak6AL4MuDsyXbO6T/5z3MrORfr4jTe7swsMsG6EPZ3BgvE/GiCiKvDTFjOPJpHwqrp41n+qtaBhPb5tNOVa7Zl5eMTj+PzJ/Qp9DP0O44jjisdaxx/fB/x/cT3VfdOh+zXfV8xmzh33AU+CXzdzF4vaR7Fl+wPKaKrz5d0DkV09XtD5PUqitW6ngH8naRfNLPdwIUU8u8WCtPtSuDarva0DU2y/H0deCFFvo8rgWPN7Lped8xxnNGlm0n+JC2hUDY/B2BmTwWT76kUUdWE8rSwfyrFvN1OM7sX2AAcH1yMl5jZzWZmwGXJNX2jLoz9SDP7oaQXhqoHQvksSc8yszt73z2nalXuukmudhNhaV603AEvKje51rkjmft7MuxH5XX+40FzjcuaP5ZkE3gsrFwT82ovzNZ03Ksn0JpobLJGZNSwt4YyTL88GcqwDuSEPsV+hn7HccRxpWPN38P4hG3W6/Sdtnv/dZ9Z1T1Gjs5s2sskrU+OLzazi5PjZwM/Bz4v6XnAHcC7qYiupojGviW5fmOo2xX28/opI+nNZvb5Jm3rzCO/R/Ez4GMl5wx4+ST65jiO0xajI3e+zWZ2XM35ORTWgnea2a2SPkl9ormqaOyuRWmX8EFgakLbzKLd+hQzS32skLSg5BJnGmmiaY9lJSSadVbuzDTs1CT8ZNiPCuvSoLHqsaD1PvJIq3FcJSamQT0olOOm7NibJ1vXVNq789VpoGXDDs/ctbkot2yZWKZ9Cv200O84jjiudKxx/PF95O8p17zTfXf9mwLWVZe/jcBGM7s1HF9FIbSroqurorE3hv28vhGSvlt1Cji46X2aTER+i+K/VLs6x3GcrmAGY7vat2t2L/t3SfdL+iUz+xFFnMf3w1YWXb0O+JKkj1NMRK4AbjOz3ZK2SToBuBU4E/iLDrpyMIVb39asXnSQ7K7Opv10CnvNwrAUTvxpsITW9L7TA1IFo2pV7ka+mm1KaOmu4xp2LIO6uL1E+xw3CYcyKrH7LQrfxQXJ4gXzwv6c7Ku2f3jywmxNx+IGoayyaac//IJ2H23YUbPe/HAoM80b4JGtE/qdj2eCo0nYj+8jD6SP76/JL50mVLUdRe28y9GO7wS+GDxHfgy8meLP6EpJZwE/BU4HMLO7JV1JIdTHgDXBcwTg7RQBMgspvEY68Rz5W2CxmX07PyHppqY3qdO0TwbeRPET4GO0hPZjFK4yjuM4PaHbEZFBUJbZvUujq83sPOC8kvr1wDGT7MNZNefe2PQ+dTbttcBaSa8Li/o6PaIT7Tm/phPvkTxFE7Q8oaMGuSMrt2f2a2g5ZMQFzhfElb0WFFctnPNwq3GuYeeh70uCpr0osTnHJcKUXRuTP+1IbNrRDzt2alzTDhr25p+HstWnJ7cU/XwkdzAJZTrWOP78vcT3lf9SgeY27TrvkSYMtfbdXZv2UNFEThwraWk8kLSfpD/tYZ8cx3HYs7vZNmo0mYg8xczGzSEh98irgPf1rltOU5po2rnXSJmnQ9QYoyYZM9gvCCpkqn22NOuijGvmRqX6oFnbxtuOW7fHNeygo8aoxCVPK8pFyTTJ+A2zFK0xRC4NWYyqcPTDfqTctr1zc6tPUQmPTXONe4KmHR4V30d8P7l/dtk7dS+SydOhy99I0URoz5Y038x2AkhayN7LZDuO43SNbnqPDBtNhPYXgBskfZ7iH+BbaIV+Oo7jdB+3aVfSVmib2Uck3UUxyyrgQ557pP9UfX/LJiU6Ca7Jg0XyCbdoFpiXmAzmButFbhbJVywHWLanMEssjGrTjmAeiTaIxcE2sShZnSoG4rSbxISWmWU8hD4G+BRufXHSMZpEoDru5tES80huFoll/t46Ca4pw00ne+PmkXIaLTdmZp36IzqO40wac027kiYr15xAEfXzHyg8xGYDT5jZkh73zaG9O2DdRGTU+OK1c7J62FvDjlN/0akuBpLPTR4wJ2ikcys07PSPLSrHS0M8+NLtmwBQnPVbHDTs6OYHLU276sapph3d/0JIegxNzycZ0yj2Kk17W9S0k/5HTTtbD2cvjbssYVTTCck6RlluudAup4mm/WmK3LJfoXBOPxN4Ti875TjOaDOMC/t2i6bmkQ2SZodQzs9Lahwn7zRnMiHqearWdL9JcE2uhUdNe25FPcCscINZifYKLc2obD3JmHhp3JS9uNBhFy0qyvlpCrJoy+Z0T+IAABZxSURBVK6yaScP2JkltMpD0qNCn2raVVr4Y6GvLefAvVegrLJpTya4JqUboe7DRpcXQRgamgjt7SFe/9uSPgJsAvZpc43jOM6kcU27miZC+wwKResdwHsoUha+rpedcqqp07DzNlW27VSTz7Xv+IXIbdyl2n9otCfYhuOK5RMWY8804cVhtbGFIZZmYRagU+wXN5o1K7UUJybtkqXP4gIGedrYPEQdWl4i0YYdNey4TMITraaVQTW5LbtbwTWjokU3wW3a5TRx+ftJ2H2SIlG34zhOT/GIyGrqUrPeRc2qDGb23J70yClNzVpl7y6rb2fTTts+ldVVlbX9DTcZixp3jaYdo9XzZFOppl3l993EZp6nks1t3BPqwv2iDTtq1alNu53XyFT8tJv4bTdpO5S4y18ldZr2a/rWC8dxnAQPY6+mLjXrT6rOOY7j9Bo3j5TjwTUzlE4mJMey+lkl+0/ROfH+cZJ/V5woTCb9qswjVRkC0/0m5pGxbPKzykyyPRngE1m5PTtOV6Bs5/I3lYnIFLcETKTbiyAMEx5c4zjO4OEuf5V4cM2AU6VR12nY+bX5ca55l1E3AZlr2FU5uQF2BJU0arrzwyTg/Czp1Ozkm9guPH5XMoDdmaYdV02Pazrmia/SuqhR51p02v94riqYpi6MvZOVa/Jrq45HCde0y/HgGsdxBg4zt2lX0TS4ZhYeXDNQNAl570TTbpeQqqxuV1bGdRPTiPQdWd388PAQP8PsoMrOTa6Z08bPcCzpVHx21PrzVeXLtOe8Lg+cKWtbpVmX2bQ9MdTUMSbOXTgtOgmu2YEH1ziO0w/cpl1JXXDNqcByM7sgHN8KHBhO/4GZXdWH/jmBpila64iKS/qhd2Lnzq+J8TC7szJ1r42aalyfbl5Wxr4ouWZ2eEBVIFH6txx/Qefh+LnGnfYpT0dbl/ypnYZdpk13Y21I17zdPFJFnab9BxReI5H5wIso7NmfB1xoO47TE3wRhGrqhPY8M7s/Of6mmT0MPCzJJyIHhLow9irqtOomIdftVnlPItL3SkAVbdd5Cti0/8raRKKGnSpguX091/bzVdOhvQdIWZrVvE3dEm5NNWyXSfV00zwiaTawHviZmb1G0v7Al4HDgfuAN5jZ1tD2XOAsiq/Tu+LyipKOBS4FFgJfA95tZn3/PVA35bNfemBm70gOD8RxHKdHmBX5tJtsDXk38IPk+BzgBjNbAdwQjpF0FIWF4WhgJfCZIPABLgTOBlaEbeUUhzkp6jTtWyW91cz+Mq2U9F+A23rbLaeKdomj2p3LqYuWrLpH3J+THZdpqrFN1L7zRRfKbNqRKk07Jbdp55pwXUKnKi26rm2Vpl23CEVOJ771o0o3IyIlLQdeDZwH/F6oPhU4MeyvBW4C3hvqrzCzncC9kjYAx0u6D1hiZjeHe14GnMY0rJ1bJ7TfA/yNpDcCd4a6Yyls26f1umOO44wwndm0l0lanxxfbGYXJ8d/TjFHt29Sd7CZbQIws02SDgr1hwK3JO02hrpdYT+v7zt1CaMeAl4i6eUUPxUA/p+Z3diXnjmOM9J0YNPebGbHlZ2Q9BrgITO7Q9KJDe5V9qPPaur7ThM/7RsBF9QDRp0LYK6gNAl5z90B61zaqibYctNHup/n7Z6THZclsaqiri9NJgrbtWmSG7vJ+6nqYxluFplIF71HXgr8uqRXUcR4LZH0BeBBSYcELfsQ4KHQfiNFAGFkOfBAqF9eUt93muS4dxzH6Tt7rNlWh5mda2bLzexwignGG83sd4B1wOrQbDVwTdhfB6ySNF/SERQTjrcFU8o2SSdIEkXivGuYBholjHIGl07C2TvRuOtWrqnSOnMtPa1rtzJOtzTtdqv2lNXVTSo2DZSZTLpV166r6cMiCOcDV0o6C/gpcHrxXLtb0pXA9yk+/jUhUR7A22m5/F3LNExCggttx3EGlG4H15jZTRReIoSYk5Mq2p1H4WmS168HjulurzrHhfYQ0y6ta1rXjlRTbZcmtmy1907WoOxE087rmmjC7ezSk9Hk6/rnGnXn+MK+1bjQdhxn8PAw9kr6JrSD8f6TwKsocsu/yczuDOdWUfhRXmZmfx7qvg4cEvr4TwTbkqT5wGUUPuMPA79pZvdJOhy41MxO7NeYBonJBte0u1+dzTkvyzTtKo2615p2Xf1ktPJOPELclj11jO6GsQ8T/fQeOYVW+OfZFCGhkVUUyahOkLQ41L3BzJ5HYUM6kDBRQJETYKuZPQf4BPDhPvTdcZx+EjTtJtuo0U/zyKkUmrQBt0haGv0kaTmujzuxm1lcHnYORRS0Jff5QNi/Cvh00OJ3A1t6PooZQDuPkjKtvJPvfpUWXqc1d6JpT4apaMSdtMnrq47LGEH5MmnMfBGEKvqpaR8KpFkD0zDQqykycK03s22xgaTrKJzet9FKBTt+HzMbAx4FDjCz+83stT0dgeM4fcM17XL6KbQrw0DNbK2ZvcDMPjbhpNnJFHbt+cDL292n8sHS5ZKekPREXTvHcabO5s2biX9vYbu803vEhFEutPemp+YRSWuAt4bD2ykPD63FzHZIWkdhFrmeVpjpRklzgKfRxixiZmdQrHWJpJFxJMq/z3WBMu3MJE0mOjsxj9T1aTK063cnrnmdTC56aPreLFu2jMcff3zKOfdH9f21o6eatpldYGbPN7PnA38DnKmCE4BHY5atHEmLQz4AgmB+FfDDcDoNP309RVjqyAhixxkV9jTcRo1+TkR+jUL4bqBw+XtzTdt9gHXBvW82RcKqi8K5zwGXhzy3W5i4JJpTQ51GXKWV1wXkdBK8M92adl1dJ9qyTzj2B6OzdUtHib4J7aANr2nY9kEKF8Cycztouf85jjOEGP7PrwqPiBxB6uzTeZs6W3cn7oJVGnev/jAnY2uein3aBUz38Xdajgttx3EGEhfa5bjQHnHaeZhMVcPu9HmTZTLBQVO5lwuU3uLmkWpcaDuOM5C40C7HhbYzgU58u3OmooFPlW5o/d2+vzN53HukGhfajuMMJP6PshwX2k4tdTbtJm3bXTsd3iO9vNbpDm7TrsaFtuM4A4kL7XJcaDuOM3C4pl2NC22nY6bitjeIf4iD2CfHP5cqXGg7jjNwuPdINS60nSkzmQnIfuHa2szFP7tyXGg7jjNwuE27mulSgJwRo2lu5G5vzsylW5+xpGdK+ntJP5B0t6R3h/r9JV0v6Z5Q7pdcc66kDZJ+JOnkpP5YSXeFc58K69P2FRfajuMMJF38xzwG/Hcz+w/ACcAaSUcB5wA3mNkK4IZwTDi3CjgaWAl8RtLscK8LgbOBFWFbOcVhdowLbcdxBo44Edlka3svs01mdmfY3wb8gGKB8FOBtaHZWuC0sH8qcIWZ7TSzeykWbjk+rKa1xMxuDusDXJZc0zfcpu04zsDRoU17maT1yfHFZnZxWUNJhwMvAG4FDo5LHprZJkkHhWaHArckl20MdbvCfl7fV1xoO44zkHQgtDeb2XHtGklaDHwV+G9m9liNObrshNXU9xU3jziOM3BETbtbk82S5lII7C+a2dWh+sFkAfFDgIdC/Ubgmcnly4EHQv3ykvq+4kLbcZyBpIveI6JYEPwHZvbx5NQ6YHXYXw1ck9SvkjRf0hEUE463BVPKNkknhHuemVzTN9w84jjOQNJFl82XAmcAd0n6dqj7Q+B84EpJZwE/JSwYbmZ3S7oS+D7FXOcaM9sdrns7cCmwELg2bH1FxSTo6CBptAbsOH3msMMO47777puS//ISyV7UsO2NcEcTm/aw4Jq24zgDh0dEVuNC23GcgcSFdjkutB3HGUhcaJfjQttxnIHDzSPVuNB2HGcgcaFdjgttx3EGDl8EoRoX2o7jDCSuaZfjQttxnIHDbdrVuNB2HGcgcaFdjgttx3EGDte0q3Gh7TjOQOITkeW40HYcZ+BwTbsaF9qO4wwkLrTLcaHtOM7A4Zp2NS60HccZSFxol+NC23GcgcSFdjkutB3HGTg8jL0aF9qO4wwcbtOuxoW24zgDiQvtclxoO44zkLjQLseFtuM4A4ebR6pxoe04zkDiQrscF9qO4wwc7j1SjQttx3EGEte0y5nVrwdJOlLSzZJ2Svr97NwqSXdK+m9J3bGS7pK0QdKnJCnUz5f05VB/q6TDQ/3hkm7q13gcx+kd0abdZBs1+ia0gS3Au4CPlpxbBbwIOEHS4lB3IXA2sCJsK0P9WcBWM3sO8Angw73stOM400M3hbaklZJ+FJS9c3rR337RN6FtZg+Z2e3ArpLTis0ASToEWGJmN5uZAZcBp4U2pwJrw/5VwElBC99N8Y/BcZwZTjc1bUmzgQuAU4CjgN+SdFQPut0XBsWmfTWwHviCmW2T9EvAxuT8RuDQsH8ocD+AmY1JehQ4wMzuB17b7kGLFi3iyCOP7GrnHcdp8eijj3blPl2ciDwe2GBmPwaQdAWF8vf97j2ifwyE0DaztbS0Z2hp3hOaNThXiqTLCQJ91qxZBPP4UHHvvfdyxBFHTHc3esKwjm2YxyXpiaTqajM7o8PbXAcsa9h2gaT1yfHFZnZxcjyu6AU2Ai/usD8DQ0+FtqQ1wFvD4avM7IGGl24ElifHy4EHknPPBDZKmgM8jTZmkfCFOSP06Yn169fv07AfMwZJTzz88MNDNy4Y3rEN87jMbErjMrOV7Vs1pmNFb5DpqU3bzC4ws+eHranAxsw2AdsknRDs1WcC14TT64DVYf/1wI3B7u04jlNGVPQiqRI441C/5J2kp1PYrZdQzB88DhxlZo9VtD8OuBRYCFwLvNPMTNIC4HLgBRQa9qpoq2rYjylrAYPIsI4LhndsPq7+EH6R/ytwEvAz4HbgjWZ297R2bJL0zaZtZv/ORJNHu/brgWNK6ncAp0+hK1dP4dpBZljHBcM7Nh9XHwgOC++gsJPPBi6ZqQIb+qhpO47jOFOnn8E1juM4zhRxoe04jjODmPFCW9Ilkh6S9L2k7hmSbpR0TQyLr8pZEs6tlnRP2FYn9Tel7fpJWdhtxbjeFnK0fFvSN9NIr0EcVxUV491f0vWh/9dL2i/Unyjp0mnq5wJJt0n6jqS7JX0w1H9A0s/C5/BtSa8K9fMkfT58Rt+RdGJyrxMlrZf0kaTuc6HddyVdlXzOCjl4NoRzL0yuua9f40+ZKZ/Z0GFmM3oDXga8EPheUnc+cDTwa8DbQt1/BS4K+6uAL4f9/YEfh3K/sL9fOHcTcPg0jGk28G/As4F5wHcowm/LxrUkue7Xga8P6rgmMd6PAOeENucAHw77JwKXTlNfBSwO+3OBW4ETgA8Av1/Sfg3w+bB/EHAHMCscf5nCO+pjwJEln+fHk/G/isKLSuF5tybt7vPPbHS2Ga9pm9k/sndwzWxaqQmiY31VzpKTgevNbIuZbQWup5WcagtFTpN+Mx52a2ZPATHsdq9x2USXyX1oBQ0M4riqqBpv+pmtpZV/5imgO7HSHWIFj4fDuWGrm80/CrghXPsQ8AhwXDg3i1aajQmfZ/huLkzufSpwWXj+LcDSkKMH4OddGFqnzJjPbNiY8UK7gk8DnwXeBnwh1E3IWULxBTqA8hDXQ0O711qR06TfVPWpbFxIWiPp3yi0nHe1ucd0jquKqr4ebEWgFaE8KOx/y8ze3fdeBiTNlvRt4CGKf4y3hlPvCKaLS6JZgEIDPVXSHElHAMfSCvT4K+BbFJr3D5L7fx74d+BI4C9Cdd3n+aKuD7I9M+ozGyaGUmib2U/M7GVm9mtmti1UV4WyDmKIa2mfKsaFFZGnvwC8F3hf3T160NduMJP6ipntNrPnU8QdHC/pGIpUwr8APB/YRGHyALiEQqCtB/6cQkiPhftcZ2YvNLP/nt3/zcAzgB8AvxmqB+0dDVp/RoahFNoVjIeyZjlLBjHEdbJ9uoLWz9FBHFcVVX19MJoAQvnQNPStEjN7hGJ+YKWZPRiE+R7gLynMB5jZmJm9x4pUDqcCS4F7Gtx7N4XN+3WhatA+zxn5mQ0DoyS0q3KWXAe8UtJ+4SftK0PddHI7sELSEZLmUUycritrKGlFcvhqWgJhEMdVRdV4089sNa38M9OGpAMlLQ37C4FfBX6Y2JcBfgP4XmizSNI+Yf8VwJiZlaYEDR4iz4n7FBPOPwyn1wFnhjYnAI9GM8Q0MWM+s6FjumdCp7oBf03xc3QXxX//syraLQC+AmwAbgOenZx7S6jfALx5uscU+vQqinwJ/wb8UU27TwJ3A98G/h44epDH1cl4KeYcbqD4R3QDsP8A9PO5wL8A36UQzH8c6i8H7gr164BDQv3hwI8oTB1/BxxWc+9ZwD+H+3wP+CLBm4TCHHFBeD93AccNwLuYEZ/ZsG0exu44jjODGCXziOM4zozHhbbjOM4MwoW24zjODMKFtuM4zgzChbbjOM4MwoW24zjODMKF9ggh6WBJX5L0Y0l3SLpZ0m+0ueZwJWlvO3zemyQ9Izn+qzR1bJtrT5T0t5N5blMkfSuUh0t64ySuf5OkT3e/Z45TjQvtESFE2P0N8I9m9mwzO5Yiiq3xup2T4E0UOTQAMLPftYpowOnAzF4Sdg8HOhbajjMduNAeHV4OPGVmF8UKKxJQ/QWMa5v/JOnOsL0kv0FdG0l/kCT6P1/S6ylSkH4xLAqwUMXiC8eF9ivDPb4j6Yamg5B0kqR/Cc+6RNL8UH+fpA+Ge94l6chQf2BIxn+npM9K+omkZeFcTLF6PvCfQj/fk2vQkv5WYfECSW+W9K+S/gF4adLmQElflXR72MbPOU43caE9OhwN3Flz/iHgFWb2QorMcp9q2kbSKRSJql5sZs8DPmJmV1FktvttK5IlPRlvIulAiqRKrwvtT28yAEkLgEuB3zSz/wjMAd6eNNkc+nYh8Puh7v0UeWZeCPwf4Fkltz4H+KfQz0/UPP8Q4IMUwvoVFLmyI58EPmFFmtTXUaRddZyuM2e6O+BMD5IuAH6FQvt+EUUy/09Lej7FAgm/WHJZVZtfpVidZTuAmeWLUuScQGGmubdh+8gvAfea2b+G47UUK8P8eTi+OpR3AK8N+79CkcAJM/u6pK0Nn1XGi4GbzOznAJK+zMR3cFRhhQJgiaR9LUmh6zjdwIX26HA3rTSfmNmaYCZYH6reAzwIPI/iF9iOkntUtRGd5VLutH16XR07Q7mb1ne73TVljDHxV+iCZL+q37OAX05/UThOL3DzyOhwI7BAUmpOWJTsPw3YZEU+6DMoljbLqWrzDeAtkhZBsbhrqN8G7Ftyn5uB/6xiJZe0fTt+CBwe05eGPvxDm2u+CbwhPOeVFOtl5uT9vA94vqRZkp5JyI1NsR7kiZIOkDSXiWadbwDviAfh14jjdB0X2iOCFekcT6MQlvdKuo3CvPDe0OQzwGpJt1D85H+i5Dalbczs6xTpSNerWIYr2pMvBS6KE5FJX34OnA1cLek7FMn+yzhJ0sa4AS8A3gx8RdJdFGsrXlRxbeSDFHnF7wROoUjjm5ssvguMhUnR91CkR72XIgXqRwlzAVbkr/4AxT+dv2PiHMG7gONULDf2fYol4Ryn63hqVmeoCd4lu81sTNIvAxdasVSY48xI3KbtDDvPAq6UNItiRfC3TnN/HGdKuKbtOI4zg3CbtuM4zgzChbbjOM4MwoW24zjODMKFtuM4zgzChbbjOM4M4v8Dp/f+ksYlBfEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "gauss.plot(add_cbar=True);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All spatial models have an associated sky region to it e.g. to illustrate the extend of the model on a sky image. The returned object is an `~regions.SkyRegion` object:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Region: EllipseSkyRegion\n", "center: \n", "width: 0.4 deg\n", "height: 0.4 deg\n", "angle: 0.0 deg\n" ] } ], "source": [ "print(gauss.to_region())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can plot the region on an sky image:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAEMCAYAAAAPqefdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2deZwdVZn3v7/ubN3ZISwhQYgSQEBFYTCvKKK4RJ0RVMCgbIqivjjqqKPgOzMuDPO6jYzb4DiCLCqLisL4ioggMAsgAZFFQAIESYjEJiF7Z+l+3j/qnO7qStW9dTv33lR3P9/Ppz5169SpqnO6733uc5/zLDIzHMdxnJFBx84egOM4jlMeF9qO4zgjCBfajuM4IwgX2o7jOCMIF9qO4zgjCBfajuM4I4i2Cm1JF0laKen+VNtekm6SdI2kKaFtoqQrJS2RdIekfVP9T5P0SNhOS7XfnO7nOI4zGmm3pn0xsDDT9iHgr4HvACeHtjOA1Wa2H3A+8AUASbsAnwZeChwBfFrSzEYGIOmy4Q6+yozWecHonZvPyxkObRXaZnYrsCrT3An0h02h7VjgkvD6R8AxkgS8HrjBzFaZ2WrgBga/BFYBfSWG8dbhz6DSjNZ5weidm8/LaZhxO3sAwDeAy4A1wDtC2xzgSQAz2yZpDbBruj2wLLRhZv5GcRxn1LPThbaZPQEclWlWXtca7TUJP9eiUO/eddddR13s/oQJExiN84LRO7fROi9JSNqQarrazE5p5B4LFy60np6eUn3vuuuu680sa3Ydtex0oV3AMmBvYJmkccB0EvPHMuDoVL+5wM31bhbeMKcASLJVq7IWmtHBaJ0XjN65jcZ57bPPPixdunTyjtyjp6eHxYsXl+oradaOPGukUVWXv2uB6BlyPHCTJZmtrgdeJ2lmWIB8XWhzHGfU0V9yG1u0VdOWdDmJpjxL0jLg02Z2YU7XC4HLJC0h0bAXAZjZKknnAneGfp8zs9GnqjjOmMeAbTt7EJWkrULbzE4q2a8XOKHg3EXARc0cl+M4VWTsadFlqKpN23GcMY3hQjufqtq0HccZ8zTHpl0QiX2lpHvCtlTSPaF9X0mbUue+lbrmMEn3hUjtr4XYkZoR3K3AhbbjOBUkatpNWYi8mEwktpm93cwONbNDgR8DV6dOPxrPmdn7U+0XAGcC88MW75kbwd0qXGg7jlNRmiO0CyKxAQja8onA5bXuIWk2MM3MbguebJcCx4XTRRHcLcGFtuM4FSR6j5TZmCVpcWo7s4EHvQJ42sweSbXNk/RbSbdIekVom0MSJxIZiMYmE8FNEt29awNjaAhfiHQcp6KUXojsMbPDh/mQkxiqZa8AnmNmz0g6DPippIOpHY09rEjt4eJC23GcCtJ675EQbf1W4LCBp5ptBjaH13dJehTYn0Sznpu6fC7wVHhdFMHdEtw84jhORWl5RORrgIfMbMDsIWk3SZ3h9XNJFhwfM7MVwDpJC4K9+lTgmnBZUQR3S3Ch7ThOBWme90iIxL4NOEDSMklnhFOL2H4B8ijgXkm/I1lUfH8q6voDJHn/lwCPAteF9guBXUME90eBsxudbSO4ecRxnIrSnDD2okhsMzs9p+3HJC6Aef0XA4fktBdGcLcCF9qO41QQj4gswoW24zgVxYV2Hi60HcepIK5pF+FC23GciuJCOw8X2o7jVBQX2nm40HYcp4J4EYQiXGg7jlNB3KZdhAttx3EqigvtPFxoO45TUVxo5+FC23GcCuLmkSJcaDuOU0F8IbIIF9qO41QU17TzcKHtOE5FcaGdhwttx3EqiNu0i3Ch7ThORXGhnYcLbcdxKohr2kW40HYcp6K490geLrQdx6kgrmkX4ULbcZyK4kI7DxfajuNUENe0i3Ch7ThORXGhnUfHzh6A4zhOPv0lt9pIukjSSkn3p9o+I2m5pHvC9sbUuXMkLZH0sKTXp9oPk3RfOPc1SQrtEyVdGdrvkLRvM2ZfhAttx3EqSMw9Umary8XAwpz2883s0LD9HEDSQcAi4OBwzb9K6gz9LwDOBOaHLd7zDGC1me0HnA98oZGZNooLbcdxKki0ae+4pm1mtwKrSj74WOAKM9tsZo8DS4AjJM0GppnZbWZmwKXAcalrLgmvfwQcE7XwVuBC23GcilJaaM+StDi1nVnyAR+UdG8wn8wMbXOAJ1N9loW2OeF1tn3INWa2DVgD7NrITBvBhbbjOBWltNDuMbPDU9u3S9z8AuB5wKHACuCfQ3uehmw12mtd0xJcaDuOU0GaZx7JvbvZ02bWZ2b9wL8DR4RTy4C9U13nAk+F9rk57UOukTQOmE55c0zDuNB2HKeCNHUhcjuCjTryFiB6llwLLAoeIfNIFhx/Y2YrgHWSFgR79anANalrTguvjwduCnbvluB+2o7jVJTm+GlLuhw4msT2vQz4NHC0pENJvh2WAu8DMLMHJF0F/J7kG+EsM+sLt/oAiSdKF3Bd2AAuBC6TtIREw17UlIEX4ELbcZwK0ryISDM7Kaf5whr9zwPOy2lfDByS094LnLAjY2wEF9qO41QUj4jMw4W2U4pWLX74x9Ipxt8debjQdhyngnjCqCJcaDtDaLc7UZnn+Ud3LBK9R5wsLrQdx6ko/nWdhwvtMchwtOlWa+C1Pp7ZZ/tHeazg/+k82vprWNLCkO5wiaSzQ9tekm6SdI2kKaHt/SEF4j2S/itk3or3OE3SI2E7LdV+c6tTIjqO0y5aGxE5kmmbph3SG34TeC1J2Oedkq4liSz6a+C5wMnAt4AfmNm3wnVvBr4CLJS0C4lj/OEk/9W7JF1rZqvbNY+RTJlv6Gb1KUv8yNW6Z/Zj2UhfZyTj/8082qlpHwEsMbPHzGwLcAVJSsNOBr8yBWBma1PXTWYw+crrgRvMbFUQ1DcwmNN2FdCH4zijANe0i2inTTsv5eFLgS8Bl5GkM3xHPCnpLOCjwATg1TXuMQfAzN7aqoE7jrMzcO+RPNqpaeemLzSzJ8zsKDP7KzNblzrxTTN7HvBJ4O9q3aPug6XLJG2QtGFYIx+hdGS2Rvpk21u1jQvbcK8fzpyd1tLT00P8vIXtssbv4pp2Ee18XxelPKzHFQxWiBjWPczsFDObbGaTS47VcZxhMmvWLOLnLWynDOtG/f3ltjFGO80jdwLzQ7rD5SSZsN6R11HSfDN7JBy+CYivrwf+KVVl4nXAOa0b8sik3jdxkYbayHEjz6tFrYXI7LkyH8+iRcux99Ee4ZiNSYFchrYJbTPbJumDJIK3E7jIzB4o6P5BSa8BtgKrCblqzWyVpHNJvgAAPmdmLUs27jjOTsSFdi5tDa4JFY9/XqLfh2ucuwi4qJnjGg0MVxPuKOhT1F6r746Q9/Hsz+yHozWX0eSdCmIGW30hMg+PiHQcp5qYf63m4UJ7FFNGI85q1PU077J9ylJGw8726a/Rt5Fnur27wrhNuxAX2o7jVBMX2rm40B7hDCfsvJbWnN2PyxzX6tvImCK1NOVtBX1qadpltPOshu0JqSqI4UK7ABfajuNUEDePFOFCe4QyHG+RPM24SKMuOs67T542XmacUE57zmrc23L61tPK8+5fdNxR45zTJtx7pBCP9HUcp5o0KSJS0kWSVkq6P9X2JUkPSbpX0k8kzQjt+0raFNJC3yPpW6lrDgspo5dI+pokhfaJkq4M7Xe0OkW0C23HcaqHkbj8ldnqczGD2UAjNwCHmNkLgT8wNLL6UTM7NGzvT7VfAJwJzA9bvOcZwGoz2w84H/hCg7NtCDePjBKG6843LrPPmjyy+zLX5JlfisgzXxSZQ7L79OuOzHGeKaUstRYtnXbRPJu2md2a1X7N7Jepw9uB42vdQ9JsYJqZ3RaOLyXJiXQdSYrpz4SuPwK+IUlmVjeZ3XBwTdtxnGpS3jwyS9Li1HZmg096N4nwjcyT9FtJt0h6RWibQ5KwLjKQFppUymgz20aSZnrXBsdQGte0RyFl3fjSr+vtJ+RcM6Ggb56mnZdTFwbz6uZp2kX7LTl9txT0zWrg6WdltfA8F0APxNlJNBZc02Nmhw/nMZL+D8lb4fuhaQXwHDN7RtJhwE8lHUzttNDDShk9XFxoO45TPQzY1lrvkVBj9i+BY6Ipw8w2A5vD67skPQrsT6JZz01dnk4LHVNGL5M0DphOUkmrJbjQHmEMJ5Cllk27yIY9oWBf69z4zD59/6iKdGbG1pfZQ5LaMd22pWCffj2uoE/enIs07FoiwjXudmPQ17q/sqSFJAVWXmlmG1PtuwGrzKxP0nNJFhwfCxlG10laANxBUtv26+Gya0kykd5GYhu/qVX2bHCh7ThOFWliRKSky4GjSWzfy0iKg58DTARuCJ57twdPkaOAz0naRqI3vD+V/vkDJJ4oXSQ28GgHvxC4TNISEg17UVMGXoAL7VHCcELT06+zGvakguP067gfX3AMg5p1bCsKGd+aastq2JszfXpTfXsz+y0Fx2ntPEvUsOPfoJaocI27XVjTsvyZ2Uk5zRcW9P0x8OOCc4uBQ3Lae4ETdmSMjeBC23GcauJh7Lm40B7FZLXaPO+Oerbs7rBPa9qxrStzLruH5Pdn0bPT5GnaWW15Y6Y9ry0ej8u05/0Sydq98/zAi3CNu8WYtXwhcqTiQttxnGrimnYuLrRHCDtSaLeMTTurYWe15nQZ+6lhHzXuKZnjSSlXk0nhBuND2/iCd1w6N1BfeN3bO3QfteiNg11Zn9lH23kZ3/Ey1NP1avl2OzuAtdZ7ZCTjQttxnGrimnYuLrQdx6keXm6sEBfaI5wyVdMbSRiVNY9Ek8dUBommkulhPy1cNCXYSbq7B/vG19FMMm7c0H38BZz27tqyZeh+Y7CHrA82kHXrB/t2B/tFdlE0mknWhX0Zl8j0AmeWbCKq7L08yVQL8MK+ubjQdhynerj3SCEutMcAZTTtooXIqMGmlOcBDXtG6DRjRrKfNi3ZT5022HdqRvueEB7QmXnnpZWq+FndEDTsTRlNe+3awb7PPjt0H10Mo6Yd55O6pO5CZJ7GHYcXh12keTtNxM0juZQS2pKOJMkXu0+4RoCZ2XNbNzTHccYsbtMupKymfSHwN8BdDM3t41SQelXS8+o9ZjXubGj6FAaJNuyoYe+yS/4xDGrf0d6ddQHsDANIe3dlXf6yGnbUqtP3HbCdhywRE3qHzmM4hRny2rIpX2td77btHcSFdi5lhfYaM7uufjfHcZwmYLifdgFlhfavJX0JuJrB/D2Y2d0tGZXTVGoF5mS9R7IadtaLBAa126hFRw17992TfVrTnhleT5kWnjol+KFMCtbnjmyyVmBbEtQ+vTd5q+2xPvEBWZOxX6fHkNXkJ/Qk+86gnTcSUFNL095uqCXu4xr3cLDBn1zOEMoK7ZeGfbo6hAGvbu5wHMdxcJt2DUoJbTN7VasH4rSeMr7c2UIGAxp3KjQ92o+jl0jUrOM+atwAE2dNHXoyqsbdXWzq7+ehbX1s7O+jt7+fzf3G7uPH87yJE5k5bhz0bkr6BqP29KBiT+/pGbj/tGn9Q8Y0oGmH8XaEiXWmtPMs2eLC/TXOZdvz/qauYTcJF9q5lPUemU6SOPyo0HQL8DkzW9OqgTmjk2Vbt/K9Z1fzq02buGP9euZNmsSUzk66OjqY0NHB01u2sKS3lxmdnby0u5sTZs7kr8J5Z4zhQjuXsuaRi4D7gRPD8SnAd4G3tmJQzuhjfX8/X1i7ln9dvoxF06fzkT324Kj585nW1b1d3/6+Ph7r7eU/Vz3Dd3p6eN+GDRw7ZQrv7ejgyAkTcu7ujDrcPFJIWaH9PDN7W+r4s5LuacWAnNZTJqQ7Lg/GYJVJqSTZA+aRzIJkXHQcMIkA7LkHv+nt5a0rn+ZV06Zxz4tfzt5dXYMXTQrhOylNuqO/n/2A/Xo38S5gRc8zXLliBe987FFeKPH5vffmoInJyGZOWxnGl9g/smaRPAW9L5hKsqaPrZnjvD61TCkD4y845+aSBnHvkVzK/ubcJOnl8SAE22xqzZCc0cR/bdrEm556igvmzeOy/eYnArtBZk+axEfmzeOhF7+Eo6dP5+hly3jPn/7EM30eMjBqiWHsZbYxRllN+wPAJcG2LZLilae3alBO+1Fmn81LPT4nR3bUuKO73YBbX1h0XLJ1K29buZIf7H8Arz3wwOTcrF2T/YCvXtDKlfNW7AsZo8JC5KRnn+WjBx7Iu5cv59OP/IGXLH+Ky5+zNy+bMoWuSSsA2HdCsh+XuV36l3Z83RfcAaOGHSvZpL8KstVssuHreZXcizQhX6BsEE8YlUspTdvM7jGzFwEvBF5gZi82s9+1dmjOSOfjz67mo3vO5rXRkbtJzBg/nq8edDDfOPhg3vLYY3zp6afpN2vqM5ydTLRpl9nqIOkiSSsl3Z9q20XSDZIeCfuZqXPnSFoi6WFJr0+1HybpvnDuawpl3CVNlHRlaL9D0r5N/VtkqKlpSzrZzL4n6aOZdgDM7CstHJuzE4i27O1cAVPvlKjFRvvxgL07as3TpvHr9eu5d9s2rjjwAOjsHNSwd98j3DBG4MTsUunKknEwQX+dviHsVw/cH+Cv9tuP3+y+G4vuuIPfAN/fbz4TghF7LsuBwV/P6c92NvXrlhDyHjXudL3Kepp2Xli7a9JNonkLkRcD3wAuTbWdDdxoZp+XdHY4/qSkg4BFwMHAXsCvJO1vZn3ABcCZwO3Az4GFwHXAGcBqM9tP0iLgC8DbmzX4LPU07Zg6eWrONqXoIsf57MqVnLfHnkzqzIl4bCL7TJ7Mza98JX0Gx/3hYTa5xj06aKKmbWa3kph00xwLXBJeXwIcl2q/wsw2m9njwBLgCEmzgWlmdpuZGckXwHE59/oRcEzUwltBTU3bzP4tvPyVmf13+lxYjHTGIFHTjulVB+zdITR97cQJ3NW7ibfsuceg7TruBzTsGIETTSdp179sAtRNQ/vOTOkLkyYyEbhq/AROvv02Tli7lqv3P4AJ4cO815bExh21ati+9mTvU8k+5mdIp2bN2rvjVOPI4kjziiCUtW1nr3cCrXX528PMVgCY2QpJ8Q05h0STjiwLbVvD62x7vObJcK9tktYAuwI9ZJDUARxvZlcNd+BlvUe+XrLNcbhl7VoWTJnCpDYGxIzr6OCyBf+LTon3PPYo5hr3yMYsqfZcZoNZkhantjN34Ml5GrLVaK91zfaNZv3AB4c3tIR6Nu3/BbwM2C1j157GoPnTaQNFYdNltLYyYdrxHdaX0ydL1n22c0C1TN4SN2/YyKtmzoQJEwf9sKO9e8CGHTXsYOseUtAs6rNxNFszfVJug13hLfycxOPl8le+kpfedBMXzujmPbNm0RVC4ffsHYxjz5Yvi8cbQ5e0L2vUuqPFPWrcRTbu9Gu3be8ARiOado+ZHV6/2xCeljQ7aNmzgZWhfRmwd6rfXOCp0D43pz19zTJJ40jqhGTNMWlukPRx4EpgQ2w0s1rXDFBPFZpAYrsex1B79lrg+DIPcMYez2zbyl4TJtbv2AK6x43jhwsWcM5TT3FvlMbOCMQSl78y2/C4FjgtvD4NuCbVvih4hMwD5gO/CaaUdZIWBHv1qZlr4r2OB26y2j/13g2cBdxKUqPgLmBx2YHXs2nfAtwi6WIze6LsTZ3WU0ubK9LwyiREyuq2A/uUI3L8nMR9VvPe1N/PpI7wizGaSAb8sItKBqfdAmOfOJPNmfbxbEdXGOBemzlwrzmc//TTnLhkCXfNns3kjg5mpmqUzVqbDDg2ZQsGb0zNNYr93sw+atxZ63t61GVt204BTbJpS7ocOJrEjLKMJI/S54GrJJ0B/BE4AcDMHpB0FfB7kn/rWcFzBJJ4lYtJfupdFzZIisRcJmkJiYa9qNZ4zGzejsynbHDNxpBP+2BSvllm5qlZne3o7e+nKy9Pdhs5ec5c/mPlSv559Wr+Yddd61/gVIvGzCO1b2V2UsGpYwr6nwecl9O+GDgkp72XIPTLIOkE4Bdmtk7S3wEvAc41s9+Wub7sStH3gYeAecBngaXAnWUH6YwtJqiDTf07P8T88wccyNdWr2blGAx1HvmM6jD2vw8C++XA60ncBb9V9uKymvauZnahpA+nTCa3DGOwDSFpIfBVkkXP7wRH+F1IDPj7knx5nGhmqyUdDZxuZqe3elwjjSIzCWy/oBbNIX3ZfeqzEd3n4udl4FyoOPPciRN5bNOmOn602aqU6eCauNA4vuA4TRxxWD6cGWwcu69lHnDyww/xuY0b+casWQNXzFiVrDnFQM1YCScmwFqXyr2dreCTrVqfNZPA9gE3ZRYkfdEyg9loThgVP1ZvAi4ws2skfabsxWU17fjJWCHpTZJezNCV1KYjqRP4JvAG4CDgpBCtFCOZ5gM3hmOnQjy/u4v7K7II+Hdz9+aKnh6e9ORSI48mBddUkOWS/o0k1fXPJU2kgYp4ZTXtfwzJoj5G4p89jaQ6eys5AlhiZo8BSLqCJPLoWJJFBUh+VtwMfJJE6fGiDBS7+NVyT4vnomiLGuTAAlwq4iRq2hs2Dj0XazoeNn48561bl1SeidVnYvKnzqxuH5+Yfs9Gjbood3Y60DwObF3YhxQSQY2etdtuvHX2bC7ftIlPhLbpMxJNO1tfcqCye0rTjsul8e5xRNnammWqvbs23QBNtGlXkBNJQuC/bGbPBpfDvy17cdlyYz8LL9cA7So9NhBlFFhGUqsyN5LJzP4H+J82jc2pwcGTJrG5v5/7N23aftVmJ/DOvebwod/ezSfSFYedimOjNsufmW0kKZIej1cAK8peX7bc2G7Ae0nsyAPXmNm7yz5oGJSOMqp7I+kyRkmVnTKBG9lztdKIZvXeqGFvzhynNe1o+diUCVKJVdM7NmzgpClT+P6f/sT/nffcoZ1i8qeBEJaoNccnwpDgmSFkSzPA9nbvkC5nQH2ezCsmd7Pq7ru4v7OTQ7q6BgJ9pkxJxhuL52TrTAJM6B36xCJNO/1BKnL5Gyt1JXt6epC0IdV0tZmd0vCNRq+mvUOUNY9cA/wn8CuGphtuJUWRSUWRTIWEN8wpAJI8vrkNvGPaNN68fDnnmdHRutw5peiQOH7mTK5dsyYR2k5LmTVrFuvXr59cv2cNDPq2+Uc1j7JCu9vMPtnSkWzPncD8EJW0nMRh/R3ALiTRR59naCSTE6gX8l7LeyRq1lH/3ZjZw6DSHPcxSGVNsAVPf/ZZXmjGdDN+uXQpC2fOHHTRiOlVB4JpYnBN2nsk6yUSNey8gmDxXEbn7Yx5YxMh/dIZM7iqpyc5DomtJk1KNO3JQcPeLtUsMClo2lHDzj4tL4AmG3DTiBeJk2AGFfAabSqSrgd+AVxnZg8N9z5lVyx/JumNw33IcDCzbSSJVa4HHgSuMrMHSIT1ayU9Arw2HDsVQxL/MHkyn/rjE5UoUPCi7smV8WhxytFv5bYRxGnAauAzku6WdIGkYyU1lOa6rKb9YeBTkjaTqDsCzMym1b5sxzCzn5MkG0+3PUNBJNNYopZ9NNunjE17S2afDdtOi7t1GQ07KtEDynRPkpHybWZ8eeJELn/8Md45M3h1DFQBju/TMj7Y0aIcZ5BOnppN4ZQhhNHP6+7mj5s30yfROS55Vkwtmy3q0Jn6VHRm9tkybHneI9Roc8oxGp1HzOxPJGHwF4cUrS8lcWn+hKRNwC/N7Iv17lPWe2Rq/V6OMxRJfHHWLE770594W19fywsi1CKmid3S31+4zOlUiFFoHkkTUrTeFrZ/kDSLJDqyLvVSsx5oZg9JeknBg+9udLDO2OKo7m6OmDSJv33oQb5+8M5zADQztpi1Nce3s2OMNk27FmbWQ5IupC71NO2Pkbj6/XPecwBPGFURar2/a+XTLlqIzJpF1qeu6Q6dozkketcNFqkZfMLMaSv5t44ODluxgld2dHB8d9Bzw2LgQD7sIUSzSBxFNEpETT1VhmYg7CWTj8+GFofs7e9ngjRQ3xQG84Ars0/Xwxyf2RfV0KwVXFN0XCYz41jFDLaNYk17R6iXmvW9Yd+ugBpnFDKjo4Mr583jjY8+yosPOIDnTWl/edFntm1j2riySzjOzmY02rSbRT3zSM2AFDO7utZ5pxrUWoiMGl428VF2ITKtaUebcNS0Y1DKQBh4qtzjpElJp8MndfEPXd0ce8st3HLQwewaFgN5TrxpelQx8CZq0fGJedmro4YdY89jkuxwbQijv/WZHl4+ZQps2TxgLC2swJNDPBU1bWXa8zTtehp2rWpDY17zHuU27SyS3mVm3y3Tt57q8Vc1zhmpUEzHqcdZ3d38saODhQ89yK/22YfpE4pyizSfG9au5TXTWurs5DQRY8S58+0onwV2XGib2buaMhynZTRS2bsRm3bcRz027ZAXY0+iI96kUNkuBqUMCQMPcnnfCUlqhS/svRcfAV59/S+4/sDnMyukc2WvVBh7TK8akz/F0PTct2v8LRA07K1hMNEfce1aNvX18Ytnn+WcmbvAxk0QElttDZOMqWVHZmrmUYqNPvOIpHuLTgF7lL1P2dwj/wR80cyeDcczgY+Z2d+VfZDjQOIG+C/TpvN/JI584H5+OHMGL5wxo/6FO8BXn3iCl3V3s//EnVO30mkcs4H07KOJPUjc+lZn2kUDye7Krsy8wcw+FQ9C0YE3Ai60RxC1akQWeY/kJUTKelTEpEoTesI+ZfWIHnZxDXAuywH4p7mzef7EiRxz883837324oyDDxn07Ng9aMlRmEdjeWeOOSV6iUQbdtSwe54B4OmVf+bLjz3G7bvtnioKmfSNSbCyRR3S9TCLamg6rWcUmkd+Bkwxs3uyJyTdXPYmZZ1WO0Oi7viALoamWnOchjmlq5tb58/nqytXcsrv7qFny5b6FzXA5r4+3rHkEd67++7sNz4v4tKpKtF7ZDTVQDCzM8zsvwrOvaPsfcpq2t8DbpT0XZK/57tJChA4FaERn988m3bW46FI44btCwIM+DAHRTYdv5KNZYna7F5bEhv38/faxB0TJ3L2xo0ceNONfGyXXfjw8/aju7MzpWkHm3ZI/jTkpvFTG4stBG16wzOrOOGRP7BLXx//2NUFTw8mg1wfqrHHhFcbM8Uctqa+O7Jl2Cyzd1rEKLRpN4tSmsJ45VoAABhqSURBVHaIh/9H4PkkFdnPLRMj7zhl6Jb42u67c9tznsPdvb3s/9u7+eaKFawahuZtZvy/lSs5/L772GP8eH6w55507uTUsM7w6O8rt9VD0gGS7kltayV9RNJnJC1Ptb8xdc05kpZIeljS61Pth0m6L5z7mtT+N1cj0QYPAtvM7FeSuiVNNbN1da9y2k4jqVkjWY17S+Y4zw95bUGfzlS5roFn9w/dR3m8Z2/SeebatcwHfgjcPn06X3lqOZ96YilHd3Xxzt135/XTpjF9crBt54Sir964gZvXruWCFX/iyW1b+fK0abxxUhcKGvbmnsG36uqMg8m6Ao0bBkszRNlQFF3qSmFzaabLn5k9DBwKA7VnlwM/Ad4FnG9mX073D7VoF5EoqHsBv5K0v5n1ARcAZwK3kySzWwhc15yRlqOs98h7SQa6C/A8klJg38Kz7TktYMH48Vw1fTprp07lJ+vXc2FPD+9+4gn2mjCBA7u62DeYSTb199Pb388fNm3iwY0bOHLqVI6fPp13zZzJ+HWuT4xkWug9cgzwqJk9UUNJPha4wsw2A49LWgIcIWkpMM3MbgOQdClwHFUU2sBZJIV27wAws0dibUbHaRXTOjs5bfp0Tpszh21mPGzwh02beGLLFjqUZO7r6ujgPXvuwYKJE5nY0ZH4YTsjn8Zs2rMkLU4df9vMvl3QdxFweer4g5JOBRaTuDGvJlFKb0/1WRbatobX2fa2UlZobzazLfGbSdI4fC1mxFCm/mBeiHv6mlph2rUWRvqCqSRrFokmiGiSmLV2cFQzViUmjVg1nSlTGQccPGkiBwOMy3qCrBs0bobAmejWFxcdo0kEYNWqoft1a4eOpTdlSo+Wkq2ZfS2ziJtKmkMD5pEeMzu8XidJE4A3A+eEpguAc0lk2bkkifHeTXF92qbVrd0RygrtWyR9CuiS9FrgfwP/0bphOY4zlrHWeI+8AbjbzJ5OnpHsAST9O4kfNRTXp10WXmfb20pZoX02cAZwH/A+EgP8d1o1KKc5NFLhuyiCO89/o5GM1ANBKUGbzWra2TqTMOjpF1OFxKrpsabjkMoyYTAx+dPWEvePGnZMeDUQcxMXJFPjzybOipp2NvzfNe7m0wKhfRIp00gsEB4O3wLcH15fC/xA0ldIFiLnA78xsz5J6yQtIDEVnwp8vemjrEPZyjX9kn4K/NTM/tziMTmOM8ZpdmFfSd0kNWXfl2r+oqRDSUwcS+M5M3tA0lXA70m+l88KniMAHyApGdZFsgDZ1kVIqJ+aVcCnSQrsKjT1AV83s8+1YXxOE2hE487atodT+zAvTD5qqFuCytobflRGO/L6VO7XbHGFrpDqNVZNT6fFVmYwMflT1Ojj/del7r9uMJfUkOetDdek09BmNexsRcq8tYDs39fdAodHM4sgmNlGYNdM2yk1+p8HnJfTvhjYeSWYqP/5+whwJPAXZrarme1CUozySEl/0/LROY4zJomadjOCa0Yb9cwjpwKvDfXLADCzxySdDPwSOL+Vg3PaR1ZjzCs3EOnNacu7F2yvocbjGLSyMWi5aU14atCws0UVYiKqdEKqqGlH2/a2TJrVrJdK+vWAvTsMbk08nxp/dCCMc86mrs2zabtm3Rw8jD2fekJ7fFpgR8zsz5I8A4/jOC1hDBZBKE09oV0r+UNzU7I5LadWwYSikPdadQGi9lkrpDub0jRbsjdqshtTD1oXy5iFfba4Qtp7ZHzBOzimV816k8CgH3a2aHE8TsdSZsoFb6dhl9G03cY9DDxhVCH1hPaLJK3NaReDBUwcx3GayigtgtAU6pUb66x13nEcp1W4eSSfRrL8OaOERnJvRxopn1irMk42P3U0O6QX/+JPuFjUPVbGmRT2aU2iaGEl+9xUBcrtAma2M9XU6Fu0EJnn8ldkFqk13jJ9xwKxCIKzPS60HcepHk0OrhlNuNAe4zQSeBPJLr6VWYgsqvoeteq0G2FsiwuCsa5d9PTL07SLFlKjRp/+/Ge1/OxxeixlNey8akDZsTiN4Zp2Pi60HcepHGZu0y7ChbYDDC/UvVYATrZvVjONWnO2FmX6XHYfNey0Hbszs49EzTp+7vMCfrKBP9l9rXPZ+dSyaXuVm8YxBgOknKG40HYcp3q4TbsQF9rOEIaTzrXInpx+nd0X2bhhULMeV2effnZR4ag8TbveGNIKXj3NupYdv55G7elca+PmkXxcaDuOUzlaVARhVOBC28mllsZdpIXnad5FGnbWHp5+I27JtMV9R07fMiXPsmPNjqVI885rK3ON27Kbg5tH8nGh7ThO5TBrbj7t0YQLbacmtZJMNUJWCy/SwKFYsy5TZLiIMsWM8zTiIs266DjvPh7t2DgeEVmMC23HcaqH27QLcaHtOE4lcZt2Pi60nVLkLUTmnSvqW7TPM79kTSl5C5BF9y9DvfzWZRYty1xT77llz41F3HukGBfajuNUEvfTzseFttMwWQWoltZcpNUWad7Z1+njHakQn0eZijJltfLhBMq4IllMs4sgSFpKkoOsD9hmZodL2gW4EtgXWAqcaGarQ/9zgDNC/w+Z2fWh/TDgYqAL+DnwYTNr69fLcN7rjuM4Lae/v9zWAK8ys0PN7PBwfDZwo5nNB24Mx0g6CFgEHAwsBP5VUkxvcwFwJjA/bAt3dJ6N4pq2s8MMJ71r0bXZ1+njZmnYRWMoOk63Dafeo2vUjdOmwr7HAkeH15cANwOfDO1XmNlm4HFJS4AjgrY+zcxuA5B0KXAccF3LR5rCNW3HcaqHNaRpz5K0OLWdmX9HfinprtT5PcxsBUDY7x7a5wBPpq5dFtrmhNfZ9rbSNk1bkoCvAm8kqeh0upndHc4tAj4BXGpm/xLafgHMDmP8T+AsM+uTNBG4FDgMeAZ4u5ktlbQvcLGZHd2uOTlDydMom6F952njRc8ZDo1oyY1o5cN5npNgNOTy15MyeRRxpJk9JWl34AZJD9Xom5d/zGq0t5V2atpvYNAOdCaJbSiyCPgLYIGkKaHtRDN7EXAIsBtwQmg/A1htZvsB5wNfaMPYHcdpJ41p2vVvZ/ZU2K8EfgIcATwtaTZA2K8M3ZcBe6cunws8Fdrn5rS3lXYK7WNJNGkzs9uBGfEPxuA32MC3mZmtDW3jSLJ1Wuo+l4TXPwKOCVp8H7CqtVNwGqWfYm203rYtbNnjbSXPld2Knpt3LrvlzafRv4WzPWZJEYQyWz0kTZY0Nb4GXgfcD1wLnBa6nQZcE15fCyySNFHSPBJF8zfBhLJO0oIgc05NXdM22rkQWWQnWgFcDSwGvmdmsTQgkq4n+Ua8jkRAD7mPmW2TtAbY1cyeBN7a6kk4jtMemhhcswfwk0TOMg74gZn9QtKdwFWSzgD+SPg1b2YPSLoK+D3Jd/dZZhaNNR9g0OXvOtq8CAntFdqF9iAzu4RB7XnwpNnrJU0Cvg+8Grih1n0KHyxdhgt0x2kLPT09SNqQarrazE5p5B7NTBhlZo8BL8ppfwY4puCa84DzctoXk5hsdxotFdqSzgLeGw7vJN9OVBMz65V0LYlZ5AYG7U3LJI0DplPHLBLeMKeEMXmc1U6gKCCn6Hyj922mna/MWJrVZzQya9Ys1q9fP3lH7zNW/371aKlN28y+GZzZDwV+CpyqhAXAmuhuk0XSlNQCwTgSj5O42pu2Qx0P3NTuiCTHcVpPmXWPsSjY22ke+TmJ8F1C4vL3rhp9JwPXBve+TuAm4Fvh3IXAZcHhfRWJ54kzwqj1YWtECx+OS2EjDOe+Y1GQNBtjaI51Z5C2Ce2gDZ9Vsu/TJC6Aeed6GXT/cxxnFGL4l18RHsbuVI569u9afduFC5TW43/jfFxoO45TSVxo5+NC26k8ZT68rVpRd8Gxc3DzSDEutB3HqSQutPNxoe2MCvwDPrpw75FiXGg7jlNJ/Is4HxfajuNUDrdpF+NC23GcSuJCOx8X2o7jVA7XtItxoe04TiVxoZ2PC23HcSqHe48U40LbcZxK4pp2Pi60HcepHG7TLsaFtuM4lcSFdj4utB3HqSQutPNxoe04TuXwhchiXGg7jlM53KZdTEtrRDqO4wyXZtWIlLS3pF9LelDSA5I+HNo/I2m5pHvC9sbUNedIWiLpYUmvT7UfJum+cO5rktS8GZfDNW3HcSpHkzXtbcDHzOxuSVOBuyTdEM6db2ZfTneWdBBJ7dmDgb2AX0na38z6gAuAM4HbSereLgSua95Q6+OatuM4laRZmraZrTCzu8PrdcCDwJwalxwLXGFmm83scZJi5EdImg1MM7PbQs3bS4HjhjW5HcCFtuM4laQBoT1L0uLUdmbRPSXtC7wYuCM0fVDSvZIukjQztM0Bnkxdtiy0zQmvs+1txc0jjuNUjga9R3rM7PB6nSRNAX4MfMTM1kq6ADg3PO5c4J+BdwN5dmqr0d5WXGg7jlM5mu09Imk8icD+vpldDWBmT6fO/zvws3C4DNg7dflc4KnQPjenva24ecRxnErSRO8RARcCD5rZV1Lts1Pd3gLcH15fCyySNFHSPGA+8BszWwGsk7Qg3PNU4JodmOKwcE3bcZxK0kRN+0jgFOA+SfeEtk8BJ0k6lESxXwq8D8DMHpB0FfB7EivNWcFzBOADwMVAF4nXSFs9RwCULIKOHSSNrQk7TpvZZ599WLp06Q75L0+R7EUl+/4P3FXGpj1acE3bcZxK4hGR+bjQdhyncnjukWJcaDuOU0lc087HhbbjOJXDE0YV40LbcZxK4kI7HxfajuNUDte0i3Gh7ThOJfGFyHxcaDuOUzlc0y7GhbbjOJXEhXY+LrQdx6kcrmkX40LbcZxK4kI7HxfajuNUEhfa+bjQdhyncngYezEutB3HqRxu0y7GhbbjOJXEhXY+LrQdx6kkLrTzcaHtOE7lcPNIMS60HcepJC6083Gh7ThO5XDvkWJcaDuOU0lc086no10PknSgpNskbZb08cy5RZLulvSRVNthku6TtETS10LJekJZ+ytD+x2S9g3t+0q6uV3zcRyndUSbdpltrNE2oQ2sAj4EfDnn3CLgL4AFkqaEtguAM4H5YVsY2s8AVpvZfsD5wBdaOWjHcXYOzRTakhZKejgoe2e3Yrztom1C28xWmtmdwNac04rdAEmaDUwzs9vMzIBLgeNCn2OBS8LrHwHHBC28j+SLwXGcEU4zNW1JncA3gTcABwEnSTqoBcNuC1WxaV8NLAa+Z2brJB0ALEudXwbMCa/nAE8CmNk2SWuAXc3sSeCt9R7U3d3NgQce2NTBO44zyJo1a5pynyYuRB4BLDGzxwAkXUGi/P2+eY9oH5UQ2mZ2CYPaMwxq3kO6lTiXi6TLCAK9o6ODYB4fVTz++OPMmzdvZw+jJYzWuY3meUnakGq62sxOafA21wOzSvadJGlx6vjbZvbt1PGAohdYBry0wfFUhpYKbUlnAe8Nh280s6dKXroMmJs6ngs8lTq3N7BM0jhgOnXMIuENc0oY04bFixdPLjmOEYOkDc8888yomxeM3rmN5nmZ2Q7Ny8wW1u9VmoYVvSrTUpu2mX3TzA4NW1mBjZmtANZJWhDs1acC14TT1wKnhdfHAzcFu7fjOE4eUdGLpJXAEYfaJe8k7Ulit55Gsn6wHjjIzNYW9D8cuBjoAq4D/trMTNIk4DLgxSQa9qJoqyo5jh3WAqrIaJ0XjN65+bzaQ/hF/gfgGGA5cCfwDjN7YKcObJi0zaZtZn9iqMmjXv/FwCE57b3ACTswlKt34NoqM1rnBaN3bj6vNhAcFj5IYifvBC4aqQIb2qhpO47jODtOO4NrWo6k0yX95c4eh+Ok8fel00wq4fI3HCQtBL5K8nPnO2b2+XDqxHDuaTM7t07f6Hi/GFhuZn8Z2k4HXgVsAlYA40lMNSea2ZY2TG878uYgaW+SwKM9SdYJvm1mX63i+Iso+t8UzPd0KjCvGn/3ScCtwESSz9aPzOzT4bK89+VSYB1JYNg2Mzs8tP8N8B4SD4f7gHeZWW+F5l807hnAd8KYDHg3cEAVxjyqMLMRt5F8kB8FngtMAH5HEul0OvDO0OfKWn1T9/oo8APgZ6m200kWKgBuDPtPAS+u2HxnAy8JfaaSLLYcVLXxD2Netf6/O31eNf7uAqaE9vHAHcCCvPdleL0UmJW59xzgcaArHF8FnF6l92XeuEP7JcB7wusJwIyqjHk0bSPVPDIQ4WTJt3WMcAKI4VhWr6+kucCbSLSDLNGr5c9hv4VEg9oZ5M7BzFaY2d0AZrYOeJDByNEqjb+Iov9Nrf/vTp9X0d/dEtaHbuPDFt+H2fdlLcYBXcHroZuh7mk7ff55SJoGHAVcCGBmW8zs2XC6kmMeqYxUoZ0X4TRnGH3/BfgE1U8WVne+Idvhi0m0u5FC0bwa+f/uVLJ/d0mdku4BVgI3mFmt/4cBv5R0l6QzAcxsOUlStT+SmBPWmNkvWzeDYbHduEl+Ff0Z+K6k30r6jqTKuP2NJkaqTTs3wsnMLk4dLKrVNywMrTSzuyQdnbnRdvcxs7zshO2iZkRXyIz4Y+Ajlvi9XzzQqRrjL6JoXqX/vztzXjl/d8ysDzg02Hd/IumQgvclwJFm9pSk3YEbJD1EYsM+FpgHPAv8UNLJZva9Cs0/b9wbgZeQxFPcIemrwNlm9vcVGfOoYaRq2o1EOBX1PRJ4c1hUuQJ4taTvNX+oTaFwvpLGkwiO75tZpfxjS1A0r8pHsNX7uwfTwM0MphTeDgtRwma2EvgJiVnoNcDjZvZnM9tK4vP8sqZPYAcoGPcyYFnql8WPSIS402RGqtC+E5gvaZ6kCST5uK9tpK+ZnWNmc81s39B2k5md3I7BD4PcOYQQ/wuBB83sKzt1hMOj6P/YyP+37RT93SXtFjRsJHWRCOCHCu4xWdLU+Bp4HXA/iVlkgaTu8JxjSGzmlaBo3JYEzz0ZMnRCMu4RmUWv6oxI84g1EOHUSN+qUjQHSS8nSYR1X7CjAnzKzH6+s8baCLX+NxX/nx1Jzt+dRNu8JLiRdgBXmdnPCu6xB4n5BJLP4Q/M7BcAkn4E3E2SnfS3wLcL7rEzKBw38NfA98MX7WPAu3bOEEc3HhHpOI4zghip5hHHcZwxiQttx3GcEYQLbcdxnBGEC23HcZwRhAttx3GcEYQLbcdxnBGEC23HcZwRhAttBwBJfZLukXS/pP9IRfbtFYI96l2/vqD9OEkH1bn2d5IuH97Im0PZeTrOzsaFthPZZGaHmtkhJAWTz4Ikz4SZHb8D9z2OJNd0LpKeT/I+PGpnZoVrwjwdpy240HbyuI2QClXSvpLuD6+7JV0l6V5JV0q6Q9Lh8SJJ5wWt+XZJe0h6GfBm4EtBi39ezrPeAVwG/DL0jff6kKTfh2ddEdqmSPqupPtC+9tC++sk3Sbpbkk/DNn3kLRU0mdD+32SDgztrwzjuSekEZ2ameek1HN+K+lVof10SVdL+oWkRyR9scl/d8epiwttZwghb8Yx5Cdo+t/AajN7IXAucFjq3GTgdjN7EUnJrfea2f+E+/xt0OIfzbnn24ErgcuBk1LtZ5NUN3kh8P7Q9vck+aVfENpvkjQL+DvgNWb2EpLScR9N3acntF8AfDy0fRw4y8wOBV5BUgorTfyV8YIwpkuUlBIDODSM+QXA25WUHnOctuFC24l0heRHzwC7ADfk9Hk5SRpbzOx+4N7UuS1ATI50F7BvvQdK+gvgz2b2BHAj8BJJM8Ppe0mSD51MkjgJkqx534zXm9lqknJeBwH/HcZ/GrBP6jExbWp6TP8NfEXSh4AZZraNobycRPvHzB4CngD2D+duNLM1ZtZLksVuHxynjbjQdiKbgua5D0l9v7Ny+uQVJ4hstcHsY32UyyB5EnCgkpzmjwLTgLeFc28iEdCHAXcpKb0lti/XJZIKMYeG7SAzOyN1fnN2TJYUD34P0AXcHs0mJee5OfW67Dwdp2m40HaGYGZrgA8BH1eS6D/NfwEnAgSPkBeUuOU6kuK3Q5DUAZwAvNDM9g15zY8FTgrn9jazX5OUg5sBTCGxe38wdY+ZwO3AkZL2C23dkvanBpKeZ2b3mdkXSMwpWaF9K/DO0Hd/4DnAwyXm6jgtx4W2sx1m9luSCuiLMqf+FdhN0r3AJ0lMGGuozRXA34YFvfRC5FHA8lATMXIrialjDvA9SfeR5JM+P1SC+UdgZnBL/B3wKjP7M0nF78vDuG5neyGc5SOpe2wCrsuZZ2d4/pUk1dA3Z2/iODsDz6ftlCYsUo43s94ggG8E9g8V0x3HaQNuj3MaoRv4dTCbCPiAC2zHaS+uaTuO44wg3KbtOI4zgnCh7TiOM4Jwoe04jjOCcKHtOI4zgnCh7TiOM4L4/xVYkvpQ06n9AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# create and plot the model\n", "gauss_elongated = GaussianSpatialModel(\n", " lon_0=\"0 deg\", lat_0=\"0 deg\", sigma=\"0.2 deg\", e=0.7, phi=\"45 deg\"\n", ")\n", "ax = gauss_elongated.plot(add_cbar=True)\n", "\n", "# add region illustration\n", "region = gauss_elongated.to_region()\n", "region_pix = region.to_pixel(ax.wcs)\n", "ax.add_artist(region_pix.as_artist());" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `.to_region()` method can also be useful to write e.g. ds9 region files using `write_ds9` from the `regions` package:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "from regions import write_ds9\n", "\n", "regions = [gauss.to_region(), gauss_elongated.to_region()]\n", "\n", "filename = \"regions.reg\"\n", "write_ds9(regions, filename, coordsys=\"galactic\", fmt=\".4f\", radunit=\"deg\")" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "# Region file format: DS9 astropy/regions\r\n", "galactic\r\n", "ellipse(0.0000,0.0000,0.2000,0.2000,0.0000)\r\n", "ellipse(96.3373,-60.1886,0.1428,0.2000,45.0000)\r\n" ] } ], "source": [ "!cat regions.reg" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# SkyModel and SkyDiffuseCube" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `~gammapy.modeling.models.SkyModel` class combines a spectral and a spatial model. It can be created\n", "from existing spatial and spectral model components:" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SkyModel\n", "\n", " Name : my-source\n", " Spectral model type : PowerLawSpectralModel\n", " Spatial model type : GaussianSpatialModel\n", " Temporal model type : \n", " Parameters:\n", " lon_0 : 0.000 deg \n", " lat_0 : 0.000 deg \n", " sigma : 0.200 deg \n", " e (frozen) : 0.000 \n", " phi (frozen) : 0.000 deg \n", " index : 2.200 \n", " amplitude : 2.70e-12 1 / (cm2 s TeV)\n", " reference (frozen) : 1.000 TeV \n", "\n", "\n" ] } ], "source": [ "from gammapy.modeling.models import SkyModel\n", "\n", "model = SkyModel(spectral_model=pwl, spatial_model=gauss, name=\"my-source\")\n", "print(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is good practice to specify a name for your sky model, so that you can access it later by name and have meaningful identifier you serilisation. If you don't define a name, a unique random name is generated:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "uxWqU-Zd\n" ] } ], "source": [ "model_without_name = SkyModel(spectral_model=pwl, spatial_model=gauss)\n", "print(model_without_name.name)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The spectral and spatial component of the source model can be accessed using `.spectral_model` and `.spatial_model`:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.spectral_model" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.spatial_model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And can be used as you have seen already seen above:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hUdfr+8feTkNAJHem9k4gQOgR16VUBFV11VQR1RVF21+5vXSvqCiJiARU7FlAEFCmWBBCVJlXpvTfpPZ/fHxO/m2UhySQzOTOT+3VduWDOzJy5s7hzz5lzznPMOYeIiIg/orwOICIi4UflISIiflN5iIiI31QeIiLiN5WHiIj4TeUhIiJ+y+d1gNxQunRpV61aNa9jiIiElYULF+51zpU53315ojyqVavGggULvI4hIhJWzGzThe7T11YiIuI3lYeIiPhN5SEiIn6L6PIws55mNubgwYNeRxERiSgRXR7OuSnOuUFxcXFeRxERiSgRXR4iIhIcKo8MnE11TFy4lbOpGlsvIpKeyiMDM1bs5G+fLqHPqz+waudhr+OIiIQMlUcGujS6iJeuvYQt+4/RY9RsXpy1mlNnUr2OJSLiOZVHBsyMXhdXYNbQ9nSPL8+Ls9bQc9Qclmz53etoIiKeUnlkQcnCsbzY/xLe/EsiB4+f5spX5vLUlys5fuqs19FERDyh8vDDn+qXY8bQJK5tXoWxszfQZWQKP6zb63UsEZFcp/LwU7ECMTx1ZTzjB7bEgOvG/sSDny3l0InTXkcTEck1Ko9salWzFNOGJHFbUg0+nr+FjsOTmbVyl9exRERyhcojBwrGRvNgt/pMurMNJQrFcuu7C7hr/GL2HTnpdTQRkaBSeQRAQqXiTB7clqEd6/D18h10GJ7MF79swzmdXCgikUnlESCx+aK4+0+1+fLudlQtVZghH/3CgHcWsP33415HExEJOJVHgNUpV5SJd7Tm0R4NmLduH51GpPDBT5tI1YgTEYkgKo8giI4yBrStzvR7kri4chwPf76ca8f+yIa9R72OJiISEBFdHl5fz6NKqUK8P6AFz/aNZ+WOQ3R5MYXXk9dx5qxGnIhIeLO8sFM3MTHRLViwwNMMuw6d4JFJy5m5chcJleJ4tm8C9csX8zSTiEhGzGyhcy7xfPdF9JZHKClXrABjbmjK6OuasP334/QcNYfhM1Zx8oxGnIhI+FF55CIzo3tCeWbe255eF1fgpW/X0v2lOSzafMDraCIiflF5eKBE4ViGX9OYcTc349jJM/R99Qcen7KSY6fOeB1NRCRLVB4euqxuWabfm8T1Lary1twNdH4xhblrNWhRREKfysNjRQvE8MQVjfh4UEtioqL48xs/cf+EpRw8rkGLIhK6VB4hokWNUnw1pB13XFqTCYu20nF4MtNX7PQ6lojIeak8QkiBmGju71KPSX9tQ6ki+bntvYXc+cEi9hzWoEURCS0qjxAUXymOyYPb8I/OdZm5chcdRyTz2aKtGrQoIiFD5RGiYqKjuPOyWnw1pC01Shdm6CdLuPnt+WzToEURCQEqjxBXq2xRPr29NY/1bMDPG/bTaXgy783bqEGLIuIplUcYiI4ybmrjG7TYpGoJHv1iBf3H/Mj6PUe8jiYieZTKI4xULlmId29pzvP9Evht5yG6jJzNq99r0KKI5D6VR5gxM65KrMysoe25rG4Znv36N654ZS4rtnszOVhE8iaVR5gqW6wAr9+QyKt/bsLOgyfp9fJcnp/+GydOa9CiiASfyiPMdY0vz6yhSVzRuCKjv1tH95dms3DTfq9jiUiEU3lEgOKFYnnh6ot595bmnDidSr/X5vHY5BUcPalBiyISHCqPCJJUpwwz7k3iL62q8c68jXQakULK6j1exxKRCKTyiDCF8+fjsV4N+fS2VuSPieLGt37m758u4fdjp7yOJiIRROURoRKrleSru9tx52U1+XzxNjoMT2Hash1exxKRCKHyiGAFYqL5R+d6TB7chnLF8nPHB4u4/b2F7D50wutoIhLmVB55QMMKcXxxZxvu61KXb1ftpsPwZD5dsEWDFkUk21QeeUS+6Cj+emktpg1pR92LivKPCUu58a2f2bL/mNfRRCQMhXx5mFkNM3vTzCakW3aFmY01sy/MrJOX+cJNzTJF+HhQKx7v3ZBFmw7Q+cUU3p67QYMWRcQvQS0PM3vLzHab2fJzlncxs1VmttbMHshoHc659c65Aecsm+ScGwjcBFwT8OARLirKuLFVNabfm0SzaiV5bMpKrnp9Hmt3H/Y6moiEiWBvebwNdEm/wMyigdFAV6ABcK2ZNTCzeDObes5P2UzW/0jauiQbKpUoxNs3N+OFqy5m3Z4jdBs5h9HfreW0Bi2KSCbyBXPlzrkUM6t2zuLmwFrn3HoAM/sI6O2cewbokZX1mpkBw4BpzrlFF3jMIGAQQJUqVbKVPy8wM/o2rURSnTI8NmUFz09fxdSlO3i+XwKNKsZ5HU9EQpQX+zwqAlvS3d6atuy8zKyUmb0GXGJmD6YtvgvoAPQzs9vP9zzn3BjnXKJzLrFMmTIBih65yhTNz+jrmvD6DU3Zd+QkvUfPZdg0DVoUkfML6pbHBdh5ll1wb61zbh9w+znLXgJeCnAuATo3vIiW1Uvx9Fe/8lryOmas2Mmwvgk0r17S62giEkK82PLYClROd7sSsN2DHHIBcYVieLZfAu8PaMHp1FSufn0ej05azhENWhSRNF6Ux3ygtplVN7NYoD8wORgvZGY9zWzMwYO6UFJ2tK1dmun3JHFzm2q8/9MmOg1P5rtVu72OJSIhINiH6o4H5gF1zWyrmQ1wzp0BBgPTgV+BT5xzK4Lx+s65Kc65QXFx2vGbXYVi8/HPng2ZcHtrCuXPx83j5jP04184cFSDFkXyMssLIyoSExPdggULvI4R9k6eOcvob9fyyvfrKF4ohn/1akS3+IvwHfwmIpHGzBY65xLPd1/In2EuoSN/vmiGdqrLlLvaUj6uIHd+uIjb39egRZG8SOUhfqtfvhif/7U1D3atx/er9tBheDKfzNegRZG8JMOvrczsvCfgnWOPc65z4CIFjpn1BHrWqlVr4Jo1a7yOE5E27D3K/ROX8vOG/bStVZqnr4ynSqlCXscSkQDI6GurzMpjBdAro3UDnznnEnIWMbi0zyO4UlMdH/68mWHTfuNsquPvnetyU+tqREdpX4hIOMuoPDI7SfBO59y6TFZ+d7aTSUSIijKub1mVy+uV5eHPl/HE1JV8uXQ7z/ZNoHa5ol7HE5EgyHCfh3Pu+8xWkJXHSN5QoXhB3rqpGS9e05gNe4/S/aU5jPpmjQYtikSgTHeYm1kzMxtpZovMbIeZrTezyWZ2m5npY6X8FzPjiksqMnNoezo1LMcLM1fTc9Qclm793etoIhJAGZaHmU3Fd0JfMnAFUB1oAjwJFAe+NLMsTcL1gs4w907pIvl5+bomjL0xkQPHTnHF6Lk889WvGrQoEiEy22Fezjm3K8MVmJV1zoX0zArtMPfWweOneearX/lo/haqlSrEsL4JtKxRyutYIpKJnJwk+KCZNc/oAaFeHOK9uIIxDOubwIe3tiDVQf8xP/Lw58s4fOK019FEJJsyK48twGgzW2dmT5lZo9wIJZGpdS3foMVb21Zn/M+b6TQihW9/y3DDVkRCVGZHW73gnGsGdAKOAePNbLmZPWRmNXIloUSUgrHRPNKjARPvaE3RAvm45e0FDPloMfs1aFEkrPg9GNHMmgJvAAnOueigpAow7fMITafOpDL6u7W88v1aihaI4bFeDemZUF6DFkVCRI4HI5pZtJl1NbN3gC+B9cA1AcwYFDraKrTF5ovi3o51mHJXWyqXKMjd4xcz8N2F7DyoQYsioS6zo60uA67FN6JkMfARvnEkh3MnXmBoyyP0nU11vDVnAy/MXEVMVBQPda9P/2aVtRUi4qGcbHk8jq804p1zXZ1z74RbcUh4iI4yBibV4OshSTSsWIwHP1vGdWN/YtO+o15HE5HzyGyHeTvn3KvOuT1m1tLMbgQws1JmViV3IkpeUq10YcYPbMnTV8azfNtBOr+Ywhuz13M2VePeRUJJVvd5PAL8E3gkbVEB4MNghZK8zcy4rkUVZg5tT9tapXnyy1/p8+oPrNqpjV6RUJHVi0H1A7oBRwGcc9uAYsEKJQJwUVwBxt6YyEvXXsLW/cfoMWo2L85azakzGrQo4rWslsdJ59uz7gDMTFf7kVxhZvS6uAIzh7anW3x5Xpy1hh6jZvPLFg1aFPFSVsvjMzMbDcSZ2c3ADOCt4MUKDB2qGzlKFo5lZP9LeOumRA6fOEOfV+by1JcrOX5KgxZFvJDlkwTNrCu+M80NmO6cmxbMYIGkQ3Ujy+ETpxk27Tc++GkzVUoWYljfeFrXLO11LJGIk5PL0M5wznUKWrJcovKITD+u38cDE5eycd8xrm1emQe61ieuYIzXsUQiRk7O8ygThDwiAdGyRimmDUnitqQafDx/C51GJDNzpQYtiuSGzK5hHmdmfS50p3PuswDnEfFLwdhoHuxWn27x5bl/4lIGvruAHgnleaxXQ0oXye91PJGIlWl5AD3w7ec4lwNUHhISLq5cnMmD2/Ja8jpGfbuGuWv38s+eDenduIJGnIgEQWb7PBY555rkYp6g0D6PvGX1rsPcP3Epizf/zuX1yvLkFY2oULyg17FEwk5O9nnoI5uEnTrlijLh9tb8vx4NmLduH51GpPD+j5tI1YgTkYDJrDxuyJUUIgEWHWXc0rY6M+5NonHl4jwyaTnXjv2RDXs1aFEkEDIbjLg8t4IEg04SlMolC/HegOY81zeBlTsO0eXFFF5PXseZsxpxIpITfl9JMBxpn4cA7Dp0gkcnLWfGyl3EV4zjuX4J1C+vEW0iF5LjKwmes7ISZpaQ81giuatcsQK8fkNTRl/XhB0Hj9Nz1ByGz1jFyTMacSLir6yOZP/ezIqZWUlgCTDOzIYHN5pI4JkZ3RPKM/Pe9vRqXIGXvl1L95fmsHDTAa+jiYSVrG55xDnnDgF9gHHOuaZAh+DFEgmuEoVjGX51Y8bd3IxjJ8/Q77UfeHzKSo6dOuN1NJGwkNXyyGdm5YGrgalBzCOSqy6rW5YZQ9tzQ8uqvDV3A51GpDBnzV6vY4mEvKyWx+PAdGCtc26+mdUA1gQvlkjuKZI/H4/3bsQnt7UiNjqK69/8ifsmLOHgsdNeRxMJWTraSiSdE6fPMvKbNYxJWU/JwrE80bsRXRpd5HUsEU8E9GgrkUhWICaa+7vU44s721CmSH5uf38hd36wiD2HT3odTSSkqDxEzqNRxTi+GNyGf3Suy8yVu+g4IpmJC7eSF7bURbJC5SFyATHRUdx5WS2+GtKOmmWK8LdPl3DTuPls+/2419FEPJdpeZhZPTP7k5kVOWd5l+DFCgyNJ5FAqFW2CJ/c1orHejZg/sb9dBqezLvzNmrQouRpmY1kvxu4E/gVaAwMcc59kXZf2Ixr1w5zCZQt+4/x0OfLmL1mL82qlWBY3wRqlimS+RNFwlBOdpgPBJo6564ALgUeNbMhf6w3cBFFwkPlkoV495bmPN8vgVU7D9N15Gxe+X4tpzVoUfKYzMoj2jl3BMA5txFfgXRNG02i8pA8ycy4KrEys/7WnsvrluW5r1dxxei5LN+mr0cl78isPHaaWeM/bqQVSQ+gNBAfzGAioa5s0QK8dkNTXv1zE3YdOknv0XN57uvfOHFagxYl8mVWHjcCO9MvcM6dcc7dCCQFLZVIGOkaX55ZQ5O48pKKvPL9Orq9NJsFG/d7HUskqDK7GNRW59xO+M8odjNrYmZNAB2vKJKmeKFY/n3VxbxzS3NOnk7lqtfn8djkFRw9qUGLEpmyNJ7EzJ4AbgLWAX88wTnnLg9etMDR0VaSm46ePMPz01fxzryNVIgryDN94kmqU8brWCJ+y+hoq6yWxyog3jl3KtDhcoPKQ7ywYON+7pu4lPV7jtKvaSUe6V6f4oVivY4lkmWBmG21HCgeuEgikS+xWkm+ursdd15Wk88Xb6PD8BSmLdvhdSyRgMjqlkci8AW+Evm/CXHOuV7BixY42vIQr63YfpD7Jy5l+bZDdGl4EY/3bkjZYgW8jiWSoUB8bbUCeB1YBvzf2VDOueRAhQwmlYeEgjNnUxk7ewMjZq2mQL4oHu3RgH5NK2GmU6YkNAWiPJKdc+0DniyXqDwklKzbc4QHJi5l/sYDtKtdmqevjKdyyUJexxL5H4HY57HQzJ4xs1Z/HKqbdriuiPipZpkifDyoFU/0bsiiTQfo/GIKb8/doEGLElayuuXx3XkW61BdkRza9vtxHvpsGcmr99C0agme7RtPrbJFvY4lAgTga6twp/KQUOac4/PF23h86kqOnTzLkA61GZRUg5hoXW5HvJXjr63M7GkzK57udgkzezJQAYNF1/OQcGBm9GlSiZn3tqdjw3I8P30VvV7WoEUJbVn9aNPVOff7HzeccweAbsGJFDjOuSnOuUFxcXFeRxHJVJmi+Rl9XRNev6Epe4/4Bi0Om6ZBixKasloe0WaW/48bZlYQyJ/B40Ukmzo3vIhZ97anX5NKvJa8jm4jZ/PzBg1alNCS1fJ4H/jGzAaY2S3ATOCd4MUSydviCsXwbL8EPri1BadTU7n69Xk8Omk5RzRoUUJElneYp12zvAO+i0DNcM5ND2awQNIOcwlnx06d4d/TVzPuhw2UL1aAp/rEc1ndsl7Hkjwg20dbmZm5TNolK4/xmspDIsGizQe4b8JS1u4+wpWXVOTRHg0oWViDFiV4cnK01XdmdpeZVTlnhbFmdrmZvQP8JVBBReTCmlQpwZd3t+Xuy2sxZcl2Og5PZurS7YT4ZzeJUJmVRxfgLDDezLab2Uoz2wCsAa4FRjjn3g5yRhFJkz9fNEM71WXKXW2pWKIggz9czG3vLWTXoRNeR5M8xp99HjH4rl1+PP1hu+FAX1tJJDpzNpW35m7ghRmric0XxSPd63N1YmUNWpSACcRsK5xzp51zO8KtOEQiVb7oKAYl1WT6PUk0KF+M+ycu4/o3f2LzvmNeR5M8QPMPRMJctdKFGT+wJU9d2YglWw7S+cUU3pyzgbMatChBpPIQiQBRUcafW1Rl5tAkWtUsxRNTV9LvtR9Ys+uw19EkQmVYHmY23czuNbN6uRVIRLKvfFxB3vxLIiP7N2bj3qN0e2k2I2et4dSZ1MyfLOKHzLY8/gIcAB4zs0Vm9qqZ9TazIrmQTUSywczo3bgis4a2p0uj8oyYtZpeL89hyRbtrpTA8edoqyigBdAV+BNwHN+Z5s8FL15g6GgryctmrtzFI5OWsefwSQa2q8E9HepQMDba61gSBoJyPQ8zKw10ds59kJNwuUHlIXndoROneearXxn/8xaqlSrEsL4JtKxRyutYEuICcqjuuZxze8OhOEQEihWI4Zk+CXx4awtSHfQf8yMPfb6MQydOex1NwpSOthLJQ1rXKs30e5IY2K46H/28mU7DU/j2t11ex5IwpPIQyWMKxkbzcPcGfPbXNsQVjOGWtxcw5KPF7Dty0utoEkayXR5mdnMgg4hI7mpcuThT7mrLPR1q89WyHXQckcLkJRq0KFmTky2PfwUshYh4IjZfFPd0qMPUu9pRuWQh7h6/mIHvLmDHweNeR5MQl9n1PJZe6C6gjnMuLC5Fq6OtRDJ3NtUxbu4G/j1jFTFRUTzYrT79m1UmKkqDFvOqjI62ypfJc8sBnfGdKPhf6wR+CEA2EQkR0VHGre1q0LFBOR6YuIyHPl/G5CXbGNYngWqlC3sdT0JMZl9bTQWKOOc2nfOzEfg+6OlEJNdVLVWYDwe24Jk+8azYdoguI1MYm7KeM2c14kT+I9snCeYWM6sBPAzEOef6pS2rDwzBd32Rb5xzr2a0Dn1tJZI9Ow+e4JFJy5j1624urhTHs/0SqHdRMa9jSS4JykmCWXzht8xst5ktP2d5FzNbZWZrzeyBjNbhnFvvnBtwzrJfnXO3A1cD5/3FRCTnLoorwNgbE3np2kvYeuA4PUfNYcTM1Rq0KJlO1V2U2Qoyeczb+C5lm/7x0cBofDOyGgDXmlkDM4s3s6nn/JTN4HV7AXOAbzLLKCLZZ2b0urgCM4e2p3t8eUZ+s4Yeo2bziwYt5mmZHW11HN/1yi/4EHxfJ1XJYB3VgKnOuUZpt1sBjznnOqfdfhDAOfdMhkHNJvzxtdU5y790znU/z/JBwCCAKlWqNN20aVNGqxeRLPr2t108/Plydh06wS1tqjO0Ux0KxWZ27I2Eo5wcbZWV63ic9TNPRWBLuttb8U3rPS8zKwU8BVxiZg86554xs0uBPkB+4KvzPc85NwYYA759Hn5mFJELuLxeOWbcW5Jh037jjTkbmLFyF8P6xNO6Vmmvo0kuyrA8nHPB+Lh+voPGL/jm7pzbB9x+zrLv0dFeIp4pWiCGp66Mp+fFFXhg4lKue+Mn+jerzIPd6hNXMMbreJILvJhttRWonO52JWC7BzlEJIda1ijF1/ckcVv7GnyyYAudRiQzY8VOr2NJLvCiPOYDtc2supnFAv2BycF4ITPraWZjDh48GIzViwhQICaaB7vWZ9KdbShRKJZB7y1k8IeL2KtBixEtS+VhZg3Os+zSLDxvPDAPqGtmW81sgHPuDDAYmA78CnzinFvhV+oscs5Ncc4NiouLC8bqRSSdhErFmTy4LX/rWIcZK3bRcXgykxZv06DFCJWlkwTTztN4D3gOKJD2Z6JzrlVw4wWGThIUyV1rdh3mvolLWbz5dy6vV5Ynr2hEheIFvY4lfgrESYIt8O2n+AHf107bgTaBiScikaZ2uaJMuL01/69HA+at20enESm8/+MmUlO1FRIpsloep4HjQEF8Wx4bnHM6xVRELig6yrilbXVm3JtE48rFeWTScvqP/ZENe496HU0CIKvlMR9feTQD2uI7K3xC0FIFiHaYi3ivcslCvDegOc/1S+C3HYfo8mIKryWv06DFMJfVfR6JzrkF5yy7wTn3XtCSBZD2eYiEht2HTvDoF8uZvmIX8RXjeLZvAg0qaNBiqArEPo/dZlYl/Q+QHLiIIpIXlC1WgNeub8ro65qw4+Bxer08hxdmrOLkGX8HVYjXsjqQ5kt8Z4Ebvn0e1YFVQMMg5RKRCGVmdE8oT+uapXhi6kpGfbuWact38mzfBJpWLeF1PMmiLG15OOfinXMJaX/WBprjm2grIpItJQrHMvyaxoy7uRnHTp6h32s/8K8pKzh26ozX0SQLsnWGuXNuEb6d5yFNO8xFQt9ldcsyY2h7rm9RlXFzN9JpRApz1uz1OpZkIqs7zIemuxkFNAFK/TFWPdRph7lIePh5w34emLiU9XuPcnViJR7u3kCDFj0UiB3mRdP95Me3D6R3YOKJiPg0r16Sr4a0445LazJx0TY6Dk9mugYthqSQv4Z5IGjLQyT8LN92kPsmLGXljkN0jy/PY70aUqZofq9j5SnZvhiUmU0h42tt9MphNhGR82pUMY4vBrdhTMp6Rs5aw5y1e3m0RwP6NqmI2fkuCyS5KbPL0LbP6MnOubA410NbHiLhbe3uI9w/cSkLNx0gqU4Znr6yEZVKFPI6VsTLaMsjs/Ko4pzbHLRkQWZmPYGetWrVGrhmTUaXYheRUJea6nh33kaem74KgPu71OOGllWJitJWSLDkZIf5pHQrmRjQVLlA1/MQiRxRUcZNbaoz/Z4kmlYtwT8nr+Dq1+exbs8Rr6PlSZmVR/pKrxHMICIiWVG5ZCHevaU5z/dLYM3uI3QdOZtXvl/LaQ1azFWZlYe7wN9FRDxjZlyVWJmZQ5PoUL8sz329iitGz2X5Np0QnFsyK4+LzeyQmR0GEtL+fsjMDpvZodwIKCJyIWWLFuCVPzflteubsPvwSXqPnstzX//GidMatBhsGR6q65yLzq0gIiLZ1aVReVrVKM2TX67kle/X8fWKnTzXN4HEaiW9jhaxsjXbSkQk1MQViuH5qy7mvQHNOXk6laten8c/v1jOkZMatBgMEV0eGowokve0q12GGfcm8ZdW1Xj3x010HpFC8uo9XseKOBpPIiIRa+Gm/dw3YSnr9hylb5NKPNqjPsULxXodK2wEYjCiiEjYaVq1JF/e3Y7Bl9Vi0i/b6DA8hWnLdngdKyKoPEQkohWIiebvnesyeXAbLorLzx0fLOL29xay+9AJr6OFNZWHiOQJDSvEMemvbbi/Sz2+XbWbDsOT+WTBFvLCV/fBoPIQkTwjX3QUd1xak6+HtKPeRcW4b8JSbnzrZ7bsP+Z1tLCj8hCRPKdGmSJ8NKglT/RuyKJNB+g0IoVxczdwNlVbIVml8hCRPCkqyrihVTVmDG1Pixol+deUlVz12g+s3X3Y62hhQeUhInlaxeIFGXdTM0ZcczHr9x6l28g5vPztGg1azEREl4dOEhSRrDAzrrykErOGtqdjw3L8e8Zqeo6aw7Kteu+4EJ0kKCJyjukrdvLopOXsO3qKge1qcE+H2hSIyXuj/nSSoIiIHzo3vIiZQ9tzVdNKvJa8jq4jZ/PT+n1exwopKg8RkfOIKxjDsL4JfHBrC86kpnLNmB95dNJyDp847XW0kKDyEBHJQJtapZl+TxID2lbn/Z98gxa/W7Xb61ieU3mIiGSiUGw+Hu3RgIl3tKZw/nzcPG4+Qz/+hQNHT3kdzTMqDxGRLGpSpQRT727L3ZfXYvKS7XQYnszUpdvz5IgTlYeIiB/y54tmaKe6TLmrLRVLFGTwh4sZ9N5CduWxQYsqDxGRbKhfvhif3dGah7vVJ2X1HjoMT+ajnzfnma0QlYeISDbli45iYFINpt+TRMMKxXjgs2X8+Y2f2Lwv8gctqjxERHKoWunCfHhrS56+Mp6lWw/S6cVk3pi9PqIHLUZ0eWg8iYjklqgo47oWVZg5NInWNUvz5Je/0vfVH1i9KzIHLWo8iYhIgDnnmLxkO/+aspLDJ04z+LLa3HFpTWLzhdfndY0nERHJRWZG78YVmXlvEl0alWfErNX0enkOS7b87nW0gFF5iIgESaki+Rl17SW8cWMivx87zZWvzOXpr37l+KmzXkfLMZWHiEiQdWhQjhlDk7imWWXGpK2hIZcAAAeoSURBVKyny8gU5q0L70GLKg8RkVxQrEAMz/RJ4MOBLQC4duyPPPjZMg6F6aBFlYeISC5qXbM0Xw9JYlBSDT6ev5lOw1P45tddXsfym8pDRCSXFYyN5qFu9fnsr22IKxjDgHcWcPf4xew7ctLraFmm8hAR8UjjysWZcldb7u1Qh2nLd9BxRApf/LItLEacqDxERDwUmy+KIR1qM/WudlQuWYghH/3Cre8sYMfB415Hy5DKQ0QkBNS9qCif3dGaR7rXZ+66vXQansKHP20mNURHnKg8RERCRHSUcWs736DFRhXjeOjzZVz3xo9s3HvU62j/Q+UhIhJiqpYqzIcDWzCsTzwrth2iy8gUxqaE1qBFlYeISAgyM/o3r8LMoe1pW6sMT331K31emctvOw95HQ1QeYiIhLSL4gow9samjLr2ErYeOE7PUXMYMXM1J894O+JE5SEiEuLMjJ4XV2Dm0Pb0SKjAyG/W0HPUHBZvPuBZpoguD13PQ0QiScnCsYy4pjHjbmrG4RNn6PPqDzwxdSXHTp3J9Sy6noeISBg6fOI0z379G+//uJkqJQsxrE88rWuVDuhr6HoeIiIRpmiBGJ68Ip6PB7UkOsq47o2feGDiUg4ez51BiyoPEZEw1qJGKaYNacdt7WvwyYItdByezIwVO4P+uioPEZEwVyAmmge71mfSnW0oWTiWQe8tZPCHi9gbxEGLKg8RkQiRUMk3aPFvHeswY8UuOgxP5otftgXltVQeIiIRJCY6irv+VJsv725LjdKF2RCk0Sb5grJWERHxVO1yRfn09takBumIWpWHiEiEio4yorGgrFtfW4mIiN9UHiIi4jeVh4iI+E3lISIiflN5iIiI31QeIiLiN5WHiIj4LU+MZDezPcAmIA7w9+IepYG9AQ8lF5Kdf6NwEKq/lxe5gv2awVh/INaZ03Vk9/k5eQ+r6pwrc7478kR5/MHMxjjnBvn5nAUXmmcvgZedf6NwEKq/lxe5gv2awVh/INaZ03Vk9/nBeg/La19bTfE6gGQqUv+NQvX38iJXsF8zGOsPxDpzuo6Q+m8oT215ZIe2PEQknGnLwztjvA4gIpIDQXkP05aHiIj4TVseIiLiN5WHiIj4TeUhIiJ+U3n4ycwKm9k7ZjbWzP7sdR4REX+YWQ0ze9PMJuRkPSoPwMzeMrPdZrb8nOVdzGyVma01swfSFvcBJjjnBgK9cj2siMg5/HkPc86td84NyOlrqjx83ga6pF9gZtHAaKAr0AC41swaAJWALWkPO5uLGUVELuRtsv4eFhAqD8A5lwLsP2dxc2BtWkufAj4CegNb8RUI6H8/EQkBfr6HBYTe/C6sIv/ZwgBfaVQEPgP6mtmrhNi4ABGRdM77HmZmpczsNeASM3swuyvPl9N0EczOs8w5544CN+d2GBERP13oPWwfcHtOV64tjwvbClROd7sSsN2jLCIi/grqe5jK48LmA7XNrLqZxQL9gckeZxIRyaqgvoepPAAzGw/MA+qa2VYzG+CcOwMMBqYDvwKfOOdWeJlTROR8vHgP02BEERHxm7Y8RETEbyoPERHxm8pDRET8pvIQERG/qTxERMRvKg8REfGbykPyPDM7a2a/pPt5IPNnBZ+ZbTSzZWaWaGafp2Vba2YH02VtfYHn3mpm752zrFza2O4YM/vYzPab2RW589tIpNF5HpLnmdkR51yRAK8zX9pJWjlZx0Yg0Tm3N92yS4G/O+d6ZPLcEsAaoJJz7kTassFAvHPutrTb7+O7Ns2knOSUvElbHiIXkPbJ/19mtihtC6Be2vLCaRffmW9mi82sd9rym8zsUzObAswwsygze8XMVpjZVDP7ysz6mdmfzOzzdK/T0cw+y0HOZmaWbGYLzWyamZVzzh0AfgC6p3tof2B8dl9HJD2VhwgUPOdrq2vS3bfXOdcEeBX4e9qyh4FvnXPNgMuA582scNp9rYC/OOcux3fVyWpAPHBr2n0A3wL1zaxM2u2bgXHZCW5m+YGRQF/nXFPgfeCJtLvH4ysMzKxyWpaU7LyOyLk0kl0EjjvnGl/gvj+2CBbiKwOATkAvM/ujTAoAVdL+PtM598dFedoCnzrnUoGdZvYd+GZip+2PuN7MxuErlRuzmb0+0BCYZWYA0fimqYJvCN5LZlYEuAbfbKPUbL6OyH9ReYhk7GTan2f5z/9fDN8n/VXpH2hmLYCj6RdlsN5x+C4mdgJfwWR3/4gBS51z7c69wzl31Mxm4bt6XH/gjmy+hsj/0NdWIv6bDtxlaR/1zeySCzxuDr6rTkaZWTng0j/ucM5tx3dthUfwXX86u1biuzpc87QssWbWMN3944F/AMWdc/Nz8Doi/0XlIfK/+zyGZfL4J4AYYKmZLec/+xjONRHfV0jLgdeBn4CD6e7/ANjinFuZ3eDOuZNAP2C4mS0BFgMt0j3ka3xfqX2U3dcQOR8dqisSRGZWxDl3xMxKAT8DbZxzO9PuexlY7Jx78wLP3cg5h+oGOJsO1ZVs05aHSHBNNbNfgNnAE+mKYyGQgO/oqAvZA3xjZomBDmVmHwNt8O1zEfGbtjxERMRv2vIQERG/qTxERMRvKg8REfGbykNERPym8hAREb+pPERExG//Hy8OD9eerTNEAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "model.spectral_model.plot(energy_range=[1, 10] * u.TeV);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In some cases (e.g. when doing a spectral analysis) there is only a spectral model associated with the source. So the spatial model is optional:" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SkyModel\n", "\n", " Name : source-spectrum\n", " Spectral model type : PowerLawSpectralModel\n", " Spatial model type : None\n", " Temporal model type : \n", " Parameters:\n", " index : 2.200 \n", " amplitude : 2.70e-12 1 / (cm2 s TeV)\n", " reference (frozen) : 1.000 TeV \n", "\n", "\n" ] } ], "source": [ "model_spectrum = SkyModel(spectral_model=pwl, name=\"source-spectrum\")\n", "print(model_spectrum)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Additionally the `gammapy.modeling.models.SkyDiffuseCube` can be used to represent source models based on templates, where the spatial and energy axes are correlated. It can be created e.g. from an existing FITS file:\n", "\n" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "from gammapy.modeling.models import SkyDiffuseCube" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SkyDiffuseCube\n", "\n", " Name : gll_iem_v06_gc.fits\n", " Parameters:\n", " norm : 1.000 \n", " tilt (frozen) : 0.000 \n", " reference (frozen) : 1.000 TeV \n", "\n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "WARNING: FITSFixedWarning: 'datfix' made the change 'Set DATE-REF to '1858-11-17' from MJD-REF'. [astropy.wcs.wcs]\n" ] } ], "source": [ "diffuse = SkyDiffuseCube.read(\n", " \"$GAMMAPY_DATA/fermi-3fhl-gc/gll_iem_v06_gc.fits.gz\"\n", ")\n", "print(diffuse)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Model Lists and Serialisation\n", "\n", "In a typical analysis scenario a model consists of mutiple model components, or a \"catalog\" or \"source library\". To handle this list of multiple model components, Gammapy has a `Models` class:" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "from gammapy.modeling.models import Models" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Models\n", "\n", "Component 0: SkyModel\n", "\n", " Name : my-source\n", " Spectral model type : PowerLawSpectralModel\n", " Spatial model type : GaussianSpatialModel\n", " Temporal model type : \n", " Parameters:\n", " lon_0 : 0.000 deg \n", " lat_0 : 0.000 deg \n", " sigma : 0.200 deg \n", " e (frozen) : 0.000 \n", " phi (frozen) : 0.000 deg \n", " index : 2.200 \n", " amplitude : 2.70e-12 1 / (cm2 s TeV)\n", " reference (frozen) : 1.000 TeV \n", "\n", "Component 1: SkyDiffuseCube\n", "\n", " Name : gll_iem_v06_gc.fits\n", " Parameters:\n", " norm : 1.000 \n", " tilt (frozen) : 0.000 \n", " reference (frozen) : 1.000 TeV \n", "\n", "\n" ] } ], "source": [ "models = Models([model, diffuse])\n", "print(models)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Individual model components in the list can be accessed by their name:" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SkyModel\n", "\n", " Name : my-source\n", " Spectral model type : PowerLawSpectralModel\n", " Spatial model type : GaussianSpatialModel\n", " Temporal model type : \n", " Parameters:\n", " lon_0 : 0.000 deg \n", " lat_0 : 0.000 deg \n", " sigma : 0.200 deg \n", " e (frozen) : 0.000 \n", " phi (frozen) : 0.000 deg \n", " index : 2.200 \n", " amplitude : 2.70e-12 1 / (cm2 s TeV)\n", " reference (frozen) : 1.000 TeV \n", "\n", "\n" ] } ], "source": [ "print(models[\"my-source\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note:**To make the access by name unambiguous, models are required to have a unique name, otherwise an error will be thrown.\n", "\n", "To see which models are available you can use the `.names` attribute:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['my-source', 'gll_iem_v06_gc.fits']\n" ] } ], "source": [ "print(models.names)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that a `SkyModel` object can be evaluated for a given longitude, latitude, and energy, but the `Models` object cannot. This `Models` container object will be assigned to `Dataset` or `Datasets` together with the data to be fitted as explained in other analysis tutorials (see for example the [modeling](modeling.ipynb) notebook).\n", "\n", "The `Models` class also has in place `.append()` and `.extend()` methods:" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [], "source": [ "model_copy = model.copy(name=\"my-source-copy\")\n", "models.append(model_copy)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This list of models can be also serialised toa custom YAML based format: " ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "components:\n", "- name: my-source\n", " type: SkyModel\n", " spectral:\n", " type: PowerLawSpectralModel\n", " parameters:\n", " - {name: index, value: 2.2, unit: '', min: .nan, max: .nan, frozen: false}\n", " - {name: amplitude, value: 2.7e-12, unit: cm-2 s-1 TeV-1, min: .nan, max: .nan,\n", " frozen: false}\n", " - {name: reference, value: 1.0, unit: TeV, min: .nan, max: .nan, frozen: true}\n", " spatial:\n", " type: GaussianSpatialModel\n", " frame: galactic\n", " parameters:\n", " - {name: lon_0, value: 0.0, unit: deg, min: .nan, max: .nan, frozen: false}\n", " - {name: lat_0, value: 0.0, unit: deg, min: -90.0, max: 90.0, frozen: false}\n", " - {name: sigma, value: 0.2, unit: deg, min: 0.0, max: .nan, frozen: false}\n", " - {name: e, value: 0.0, unit: '', min: 0.0, max: 1.0, frozen: true}\n", " - {name: phi, value: 0.0, unit: deg, min: .nan, max: .nan, frozen: true}\n", "- name: gll_iem_v06_gc.fits\n", " type: SkyDiffuseCube\n", " filename: $GAMMAPY_DATA/fermi-3fhl-gc/gll_iem_v06_gc.fits.gz\n", " parameters:\n", " - {name: norm, value: 1.0, unit: '', min: .nan, max: .nan, frozen: false}\n", " - {name: tilt, value: 0.0, unit: '', min: .nan, max: .nan, frozen: true}\n", " - {name: reference, value: 1.0, unit: TeV, min: .nan, max: .nan, frozen: true}\n", "- name: my-source-copy\n", " type: SkyModel\n", " spectral:\n", " type: PowerLawSpectralModel\n", " parameters:\n", " - {name: index, value: 2.2, unit: '', min: .nan, max: .nan, frozen: false}\n", " - {name: amplitude, value: 2.7e-12, unit: cm-2 s-1 TeV-1, min: .nan, max: .nan,\n", " frozen: false}\n", " - {name: reference, value: 1.0, unit: TeV, min: .nan, max: .nan, frozen: true}\n", " spatial:\n", " type: GaussianSpatialModel\n", " frame: galactic\n", " parameters:\n", " - {name: lon_0, value: 0.0, unit: deg, min: .nan, max: .nan, frozen: false}\n", " - {name: lat_0, value: 0.0, unit: deg, min: -90.0, max: 90.0, frozen: false}\n", " - {name: sigma, value: 0.2, unit: deg, min: 0.0, max: .nan, frozen: false}\n", " - {name: e, value: 0.0, unit: '', min: 0.0, max: 1.0, frozen: true}\n", " - {name: phi, value: 0.0, unit: deg, min: .nan, max: .nan, frozen: true}\n", "\n" ] } ], "source": [ "models_yaml = models.to_yaml()\n", "print(models_yaml)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The structure of the yaml files follows the structure of the python objects.\n", "The `components` listed correspond to the `SkyModel` and `SkyDiffuseCube` components of the `Models`. \n", "For each `SkyModel` we have informations about its `name`, `type` (corresponding to the tag attribute) and sub-mobels (i.e `spectral` model and eventually `spatial` model). Then the spatial and spectral models are defiend by their type and parameters. The `parameters` keys name/value/unit are mandatory, while the keys min/max/frozen are optionnals (so you can prepare shorter files).\n", "\n", "If you want to write this list of models to disk and read it back later you can use:" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [], "source": [ "models.write(\"models.yaml\", overwrite=True)" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [], "source": [ "models_read = Models.read(\"models.yaml\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Additionally the models can exported and imported togeter with the data using the `Datasets.read()` and `Datasets.write()` methods as shown in the [analysis_mwl](analysis_mwl.ipynb) notebook." ] }, { "cell_type": "markdown", "metadata": { "jupyter": { "outputs_hidden": false } }, "source": [ "# Implementing a Custom Model\n", "\n", "In order to add a user defined spectral model you have to create a SpectralModel subclass.\n", "This new model class should include:\n", "\n", "- a tag used for serialization (it can be the same as the class name)\n", "- an instantiation of each Parameter with their unit, default values and frozen status\n", "- the evaluate function where the mathematical expression for the model is defined.\n", "\n", "As an example we will use a PowerLawSpectralModel plus a Gaussian (with fixed width).\n", "First we define the new custom model class that we name `MyCustomSpectralModel`:" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [], "source": [ "from gammapy.modeling.models import SpectralModel, Parameter\n", "\n", "\n", "class MyCustomSpectralModel(SpectralModel):\n", " \"\"\"My custom spectral model, parametrising a power law plus a Gaussian spectral line.\n", " \n", " Parameters\n", " ----------\n", " amplitude : `~astropy.units.Quantity`\n", " Amplitude of the spectra model.\n", " index : `~astropy.units.Quantity`\n", " Spectral index of the model.\n", " reference : `~astropy.units.Quantity`\n", " Reference energy of the power law.\n", " mean : `~astropy.units.Quantity`\n", " Mean value of the Gaussian.\n", " width : `~astropy.units.Quantity`\n", " Sigma width of the Gaussian line.\n", " \n", " \"\"\"\n", "\n", " tag = \"MyCustomSpectralModel\"\n", " amplitude = Parameter(\"amplitude\", \"1e-12 cm-2 s-1 TeV-1\", min=0)\n", " index = Parameter(\"index\", 2, min=0)\n", " reference = Parameter(\"reference\", \"1 TeV\", frozen=True)\n", " mean = Parameter(\"mean\", \"1 TeV\", min=0)\n", " width = Parameter(\"width\", \"0.1 TeV\", min=0, frozen=True)\n", "\n", " @staticmethod\n", " def evaluate(energy, index, amplitude, reference, mean, width):\n", " pwl = PowerLawSpectralModel.evaluate(\n", " energy=energy,\n", " index=index,\n", " amplitude=amplitude,\n", " reference=reference,\n", " )\n", " gauss = amplitude * np.exp(-((energy - mean) ** 2) / (2 * width ** 2))\n", " return pwl + gauss" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is good practice to also implement a docstring for the model, defining the parameters and also definig a `tag`, which specifies the name of the model for serialisation. Also note that gammapy assumes that all SpectralModel evaluate functions return a flux in unit of `\"cm-2 s-1 TeV-1\"` (or equivalent dimensions).\n", "\n", "\n", "\n", "This model can now be used as any other spectral model in Gammapy:" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "MyCustomSpectralModel\n", "\n", " name value error unit min max frozen\n", "--------- --------- ----- -------------- --------- --- ------\n", "amplitude 1.000e-12 nan cm-2 s-1 TeV-1 0.000e+00 nan False\n", " index 2.000e+00 nan 0.000e+00 nan False\n", "reference 1.000e+00 nan TeV nan nan True\n", " mean 3.000e+00 nan TeV 0.000e+00 nan False\n", " width 1.000e-01 nan TeV 0.000e+00 nan True\n" ] } ], "source": [ "my_custom_model = MyCustomSpectralModel(mean=\"3 TeV\")\n", "print(my_custom_model)" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$1.1443958 \\times 10^{-12} \\; \\mathrm{\\frac{1}{s\\,cm^{2}}}$$" ], "text/plain": [ "" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "my_custom_model.integral(1 * u.TeV, 10 * u.TeV)" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hVVdbH8e9KhzQgJIBUCUgHkSZSoiNVBRSxDzYGxY44xTIzOpZhHN9BURHFsY5jxQZYKJYAYgEUBEQ6SqihGFoCIdnvHzdxIqbn3tyS3+d58sg959xzVoTclb332nubcw4REZGKCPN3ACIiEnyUPEREpMKUPEREpMKUPEREpMKUPEREpMKUPEREpMIi/B1Adahfv75r0aKFv8MQEQkqS5cu3e2cSy7uXI1IHi1atGDJkiX+DkNEJKiY2Q8lnVO3lYiIVJiSh4iIVJiSh4iIVJiSh4iIVJiSh4iIVJiSh4iIVJiSh0g1WLUti5zcPH+HIeI1Sh4iPjZ9aQZnP7qQl7/80d+hiHhNUCQPM2tpZs+Y2fQix841s6fN7F0zG+TP+ERKsmj9bm5/81sAlv64z8/RiHiPz5OHmT1rZrvMbOVxx4eY2RozW29mt5d2D+fcRufcmOOOveOcGwtcCVzk9cBFqmjdzgNc+9JSWibH0q91fb7N+MnfIYl4TXW0PJ4HhhQ9YGbhwBRgKNAeuMTM2ptZJzObddxXShn3/3PBvbxuf04ulz79BV9u3OOL20sIy8nN46rnFxMTGc6zV/agb6v6bNmbzd5DR/0dmohX+Dx5OOfmA3uPO9wTWF/QojgKvAqMcM6tcM6dc9zXruLuax4PAh84574u5vw1ZrbEzJZkZmZWKvZtP2Xz497DXDTtC+54awVZ2bmVuo/UPKu37ydjXzZ/Pac9TerWpnOTOgAsV+tDQoS/xjwaA1uKvM4oOFYsM0sysyeBrmZ2R8Hhm4ABwCgzG3f8e5xz05xz3Z1z3ZOTi10UskxtGyYw59b+jO13Iq8t/pGBk9L5cOX2St1LapaNmYcA6HBCAgCdmiRiBt9uyfJnWCJe469Vda2YY66ki51ze4Bxxx17FHjUy3H9Su2oCO46uz3DupzA7W+uYNxLXzOofQPuHdGRhokxvn68BKkNmQeJCDOa1qsNQFx0BKnJcRr3kJDhr5ZHBtC0yOsmwDY/xVIunZvU4d0b+3D70Lakr81k4KR0/vPFD+Tnl5jzpAbbkHmQ5km1iQz/349Y5yaJLM/Iwjn9m5Hg56/ksRhobWYnmlkUcDEww0+xlFtkeBjj0lKZc2t/OjdN5C/vrOTCpz5n3c4D/g5NAsyGzEOkJsf94tjJTeuw++ARtmfl+CkqEe+pjlLdV4DPgTZmlmFmY5xzx4AbgdnAauB159wqHzx7mJlNy8rybj9z86RYXhrTi4dGdWbdroOc9egCHpm3liPHNINYIDcvnx/2HCI15ZfJ4+dB8y3qupLgVx3VVpc45xo55yKdc02cc88UHH/fOXeScy7VOfeAj5490zl3TWJiotfvbWZc0L0pH92WxtCOjXhk3jrOfnQhSzYfX1gmNc2WvYfJzXO/anm0axRPZLixPEOD5hL8gmKGeSCrHxfNo5d05bkre5B9NI9RT37On99ZwYEclfXWVBsKKq1Sk2N/cTw6Ipy2DRM0aC4hQcnDS85om8KcW/tzVZ8W/PfLHxk4aT5zVu3wd1jiBxsyDwLQ8riWB3gGzVdkZKnQQoJeSCcPX415lCQ2OoK7h3Xg7ev7UKd2JNf8ZynXvbSUXfs1QFqTbNh1kOT4aBJrRf7qXJcmdThw5Bib9hzyQ2Qi3hPSycOXYx6lOblpHWbe1Jc/DG7DR9/v4sxJ6bz85Y/6bbOG2JB58FddVoW6NNWguYSGkE4e/hQZHsYNZ7Ri9vj+dDghgTvfXsHFT3/xc5eGhCbnXLFluoVSk2MJs//NQBcJVkoePnZi/VheGXsqD57fie+372fo5AU8/vE6jh7L93do4gN7Dx0lKzu32PEOgIjwMOrHRbPrgLoyJbgpeVQDM+OiHs2Yd1saA9s34P/mrGXYYwv5Wvs7hJySKq2KSkmIZteBI9UVkohPKHlUo5T4GKZcegr/vrw7+3NyOX/qIu6ZsYqDR475OzTxksJuyZK6rcDz72DXfiUPCW4hnTyqu9qqvAa0b8CcW/tz+anNeeHzzQyalM5Hq3f6Oyzxgg27DhIdEUbjOrVKvCYlXi0PCX4hnTz8VW1VHvExkfxtREemjzuNuJgIxrywhBtf/ppMfagEtQ2ZB2mZHEdYWHELR3ukxEez59ARjuVp3EuCV0gnj2DQrXldZt3UjwkDT2LOqp0MmJTO64u3aOXVIOWptCp5vAMgOSEG52CPdhWUIKbkEQCiIsK4+czWvH9LP9o0iOePb37LpU9/yabdKucMJjm5eWzZd7jU8Q7wtDwAjXtIUFPyCCCtUuJ49ZpTeeC8jqzcmsWQR+Yz5ZP15Kp7Iyj8uPcwzkHLMloePycPletKEFPyCDBhYcZlvZoz77Y0zmiTwkOz1zDssYWakRwECvfpaJRY8mA5QEqCZwdKDZpLMFPyCFANEmJ4cnQ3nhrdjX2Hj3LeE59x78zvOKSy3oC18+fkUfr2xMlxnpbHTq15JkEspJNHoJbqVsTgDg2ZOyGNS3s149nPNjHo4fl8smaXv8OSYuwoSAYpCdGlXhcVEUbd2pFqeUhQC+nkEciluhWREBPJ/ed2Yvq43tSKCueq5xZzy6vfsPugPnwCyY79OdSLjSI6IrzMazVRUIJdSCePUNO9RT3eu7kv4we05v0V2xkwKZ3pSzNU1hsgdmbl0CCh9C6rQikJ0WRqwFyCmJJHkImOCGf8gJN4/+Z+tEqO4/dvLGf0M1/xg/aH8Lsd+3NoWEaXVaFkzTKXIKfkEaRaN4jn9Wt7c9+5HVm25ScGPzKfp9I3aNayH+3cn0PDMgbLC6XEx5B54Ij2eJGgpeQRxMLCjNGnNmfuhP70a53MxA++Z8SUz1i5NXgLBILV0WP57D54tPzdVvHRHMt37DusWeYSnJQ8QkCjxFpMG92NqZedwq4DRxgx5TP+/v5qso/m+Tu0GqNwwl/DCox5eN6nrisJTiGdPEKhVLe8zIyhnRox79Y0LuzehGnzNzLokXQWrMv0d2g1QuGcjQYV6LYCJQ8JXiGdPEKlVLciEmtHMnFkZ1695lQiw8IY/cxXTHhtGXu1CJ9P7cjyJIFytzx+Xt9KFVcSnEI6edRkp7ZM4v1b+nHjGa2YsXwbAyal8843W1XW6yOFEwTVbSU1hZJHCIuJDOf3g9sw6+a+NKtXm/GvLeOK5xazZe9hf4cWcnbuzyEqIow6tSPLdX3tqAjioiO0f4sELSWPGqBtwwTevO407hnWnqWb9zLo4fn8e8FGlfV60Y6sHBomxGBW8iZQx/PsZa5uKwlOSh41RHiYcWWfE5kzIY3eqUnc/95qRk5dxKptoV9MUB08EwTL12VVKCU+WkuUSNBS8qhhGtepxTNXdOexS7qy7adshj/+GRM/UFlvVe3cn1PuSqtCKfExGvOQoKXkUQOZGcO6nMC8CWmcf0pjnkrfyJDJ8/ls/W5/hxaUnHMF3VblW5qkUEq8p9tKRQwSjJQ8arA6taP456guvDy2FwZc9u8v+f0by9mnst4KycrO5cix/HLPLi+UkhBNTm4+B7RHiwQhJQ/htNT6fDi+P9ednsrb32xlwKR03l2mst7y+rlMtxLdVqC9zCU4hXTyqEkzzKsqJjKcPw1py8wb+9K4bi1ueXUZVz+/mK0/Zfs7tIC3I6ticzwKaS9zCWYhnTxq4gzzqmp/QgJvX9+Hv5zTni827mXgpHSeXbiJPK3+WqKflyapRLcVoLkeEpRCOnlI5YSHGWP6nsicW/vTo0U97p31HSOnLuL7Hfv9HVpAKlyapKLJI1ndVhLElDykRE3r1eb5q3ow+eKTydh7mHMeXchDs78nJ1dlvUXt2J9DUmwUUREV+3FKiIkgOiJM3VYSlJQ8pFRmxoiTGzNvQhojTm7MlE82MHTyAr7YuMffoQWMnfvLv/1sUWZGUmwUe1TdJkEoorSTZvZ1Oe6R6Zwb7KV4JEDVjY3iXxd24dyuJ3DX2yu5eNoXXNyjKXcMbUdiOddzClU7ssq/g+DxkuKiteKxBKVSkwcQDQwv5bwBb3kvHAl0/VonM3t8fx6Zt5Z/L9zER9/v4m/DOzC0Y8MKresUSnbuz6FL0zqVem9SXBR7Dip5SPApq9vqBufchlK+1gM3V0egEjhqRYVzx1ntePeGPqTER3P9f79m7ItL2Z5V88p6jxzLY8+hoxUu0y2UFBvNnoMaMJfgU2rycM59WtYNynONhKaOjRN594Y+3HlWWxauz2TgpPm8+PnmGlXWW1gp1TCxYkuTFEqKi2L3oaOakClBp8wBczPrYWaTzexrM9tuZhvNbIaZXWtm8dURpASuiPAwrumfypzxaXRtVoe/vruKC55cxNqdB/wdWrX43+zyWpV6f1JsFEeP5XNIC1NKkCk1eZjZLOBGIB04FzgROAW4H6gDvGdm5/g6SAl8zZJq8+LVPZl0YRc27T7E2Y8uYNKcNRw5FtofioUbazWpW8nkEedpsajrSoJNWQPmY5xzO487lgN8VfD1oJml+CQyCTpmxshTmpB2UjL3v7eaRz9ez6wV2/nHyM70PLGev8PziYx9nnGexnUqmzyiANh98CjNk2K9FpeIr5XVbXWHmfUs7QLn3C4vxuNVWtvKP5Lionn4opN54eqeHMnN58KnPufOt1ewPyfX36F53Za9h2mQEE1MZHil3l8/Vi0PCU5lJY8twBQz22BmD5hZx+oIylu0tpV/pZ2UzJxb+zOm74m8+tWPDPhXOh+u3OHvsLxqy77DNKlbu9Lvr1fQ8tBcDwk2ZVVb/cs51wMYBBwGXjGzlWZ2p5m1rJYIJajFRkfwl3Pa884NfUiKi2bcS0u59j9Lfl5MMNhl7MumaSXHO8AzYA5olrkEnXItT1Iwp+MB51wn4ArgAmCdTyOTkNK5SR1m3NiH24e25dM1mQz4VzovffED+UFc1nssL5/tWTk0rVf5lkdMZDhx0RHsVreVBJlyJQ8zCzezoWb2AvAesBG4yKeRSciJDA9jXFoqs8f3p1OTRP78zkoumvY563cFZ1nv9qwc8vJdpSutCmmWuQSjskp1zzCzacBWPDPJPwZaO+fOd85Nr44AJfS0qB/Lf3/Xi4dGdWbtzoOcNXkhk+et4+ixfH+HViFb9nnKdJtWYcwDoF5sFHsOqeUhwaWslse9wDdAJ+fcUOfcC8654Pw1UQKKmXFB96bMm5DG4I4NeXjeWs5+dAFLf9jr79DKLWOvp0y3Kt1WULhEiVoeElzKGjDv55yb6pzLNLNTzexyADNLMrNm1ROihLLk+Ggeu6Qrz13Zg8NH8xj15Of85Z2VHAiCst4t+w4TZhXfu/x49eO0LLsEn/KOefwZuBv4c8GhGOBlXwUlNc8ZbVOYfWt/rujdgpe+/IGBk+Yz97vj56cGlox92TRKrEVkeNW2xUmKi2LvoaNBXTwgNU95/9WPAs4CDgE457YCCb4KSmqmuOgI7hnegbeuO406tSMZ++ISrv/vUnYFaFnvlr2HaVqvaoPlAPVio8nLd2RlB35rS6RQeZPHEedZ9tMBmFnVOnlFStG1WV1m3tSXPwxuw7zVuzhzUjqvfPVjwP1mXtUJgoXqx2muhwSf8iaPt8xsCpBoZlcBc4BnfReW1HSR4WHccEYrPrylH+0bJXDHWyu4+Okv2JB50N+hAZ59PHbuP1LlSivwDJiDliiR4FLeSYIPArOAGUAX4AHn3CO+DEwEoGVyHK+MPZV/jOzE6u37GTp5AY9/7P+y3q37Ciutqt5tlaSWhwShsvYwn+OcGwTgnPsA+KBaohIpIizMuLhnM37TNoV7Zq7i/+asZeby7fzj/E50bVbXLzFtKUge3ui2+jl5qOUhQaSslkdytUQhUg4pCTE8cVk3nr68O1nZuYycuoh7Zqzi4JFj1R5LRuEEQS+0POrWVstDgk9Z+3kkmtnIkk46597ycjwiZRrYvgGntqzHQ7PX8MLnm5mzagf3n9eR37RtUG0xbNmbTWS40SC+anM8wDO+U6d2pCYKSlApM3kA5wBWzDkHKHmIX8THRHLviI6MOLkxt7/5LVc/v4RzOjfi7mEdSI6v3H7iFZGx7zCN69QiLKy4H42KS9ISJRJkykoePzjnrq6WSHzAzIYBw1q1auXvUMRHujWvy3s39+PJ9A08/vF6FqzbzV1nt+OCbk0w884He3G27Muu8rIkRSXFRbNbLQ8JImWNefjup68aaDOomiEqIoybz2zN+7f046QGcfxx+rdc9u8v2bz7kM+embH3cJVX0y0qKTZKA+YSVMpKHqOrJQoRL2iVEsdr1/TmgfM6siIji8GPzGfqpxvIzfNuWe/ho8fYc+ioVyqtChUuUSISLMpaGHFldQUi4g1hYcZlvZozd0Iap7dJ5sEPv2f445/xbcZPXntGxj7vrKZbVFJsNPsO53LMy4lOxFeqtqKbSIBqmBjDU6O78+Rvu7H30BHOnfIZ9836jkNeKOtdvNmzbHybBvFVvlehwiVK9h5W60OCQ4WTh5nVNbPOvghGxNuGdGzI3AlpXNKzGc8s3MSgh+fz6ZpdVbrnjGXbSE2O5aQGcV6K0rM4IqByXQka5V2S/VMzSzCzesBy4Dkzm+Tb0ES8IyEmkgfO68Qb43oTExnGlc8tZvyr31RqgHpHVg5fbd7L8C6NvVrNVTjLXOMeEizK2/JIdM7tB0YCzznnugEDfBeWiPf1aFGP92/px81ntua9FdsZMCmdN5dm4FkwunxmfbsN52BYl0Zeja2w22q3Kq4kSJQ3eUSYWSPgQjwLJIoEpeiIcCYMPIn3bu5Hy+Q4bntjOZc/+xU/7jlcrvfPXL6Njo0TaJnsvS4rKLqyrloeEhzKmzzuBWYD651zi82sJbDOd2GJ+NZJDeJ549re3DeiA9/8+BODHkkvs6x38+5DLM/IYniXE7weT2KtSMLDTC0PCRrlXZL9DedcZ+fc9QWvNzrnzvdtaCK+FRZmjO7dgrkT+pN2kqesd9hjC/nmx33FXj9z+TYAzu7s/eQRFmY0TIhh60/ZXr+3iC+oVFdqvEaJtXhqdHeeGt2Nnw7nct4Ti7j1tWU/r5wL4JxjxvJt9GhRl8Z1vDezvKjUlLiA2exKpCxlrW0lUmMM7tCQ01KTeOLTDTy7cBPvfbud87o2JizM2JGVzbpdB7lvRAefPT81OZbFm/aSn++8tuCiiK8oeYgUER8TyZ+GtOXy3s2ZNGctby/bSkJMJMnx0Qxq34DhXRr77NmpyXFk5+axY38OJ/iodSPiLWUmDzNrCzQGvnTOHSxyfIhz7kNfBifiL40Sa/HQBV3456jOPl2dt6jUggquDZkHlTwk4JU65mFmNwPvAjcBK81sRJHTf/dlYCKBoLoSB0BqSiwAG3Zp3EMCX1ktj7FAN+fcQTNrAUw3sxbOuckE+XLtIoEmOS6a+JgINmT6bil5EW8pK3mEF3ZVOec2m9npeBJIc5Q8RLzKzEhNVsWVBIeySnV3mNnJhS8KEsk5QH2gky8DE6mJlDwkWJSVPC4HdhQ94Jw75py7HOjvs6hEaqjUlFh27j/CgZxcf4ciUqpSu62ccxmFfzazukDTIu/RVFgRLyusuNqYeYguTev4ORqRkpVrnoeZ3QdcCWwACpcgdcBvfBOWSM1UtFxXyUMCWXknCV4IpDrntOSniA81q1eb8DBjoyquJMCVd22rlYB+DRLxsaiIMJrXq61Bcwl45W15TAS+MbOVwM9rRjvnhvskKpEarKUqriQIlDd5vAA8CKwASt7wQESqLDUllvlrMzmWl09EuBa+lsBU3uSx2zn3qE8jKUHBxlN34dkKd1TBsXbALXjmm3zknJvqj9hEfCE1OY6jeflk7MumRf1Yf4cjUqzy/lqz1MwmmllvMzul8KusN5nZs2a2q6C7q+jxIWa2xszWm9ntpd2jYOOpMccdW+2cG4dnIL97Ob8HkaBQtOJKJFCVt+XRteC/pxY5Vp5S3eeBx4EXCw+YWTgwBRgIZACLzWwGEI5nbKWoq51zu4q7sZkNB24vuL9IyEhNLlggMfMgZ7Zr4OdoRIpXruThnDujMjd3zs0vWFCxqJ549kLfCGBmrwIjnHMT8Sx9Ut57zwBmmNl7wMvHnzeza4BrAJo1a1aZ8EX8ok7tKBokRPPpmkzG9mtZrSv7ipRXubqtzOzvZlanyOu6ZnZ/JZ/ZGNhS5HVGwbGSnp1kZk8CXc3sjoJjp5vZo2b2FPB+ce9zzk1zznV3znVPTk6uZKgi/jEuLZVFG/Ywe9VOf4ciUqzyjnkMdc79VPjCObcPOKuSzyzu1yhXzLHCZ+1xzo1zzqUWtE5wzn3qnLvZOXetc25KJeMQCVijT21Omwbx3P/ed+Tk5vk7HJFfKW/yCDez6MIXZlYLiC7l+tJk4Fkjq1ATYFsl7yUSkiLCw7h7eHsy9mXzVPpGf4cj8ivlTR4vAR+Z2RgzuxqYi2fuR2UsBlqb2YlmFgVcDMyo5L1EQtZpqfU5u1Mjnvh0PRn7Dvs7HJFfKFfycM79E7gfaAd0AO4rOFYqM3sF+BxoY2YZZjbGOXcMuBGYDawGXnfOrarsN1DG84eZ2bSsrCxf3F7E5+48ux1mMPbFpSrdlYBizpU43ICZmSvtgnJe42/du3d3S5Ys8XcYIpXyyZpdTHhtGTm5+dw9rD0X9WiqCiypFma21DlX7Fy6sloen5jZTWb2i1pXM4sys9+Y2QvAFd4KVER+7Yw2KXw4vj+nNK/D7W+tYOyLS9mRlePvsKSGKyt5DAHygFfMbJuZfWdmm4B1wCXAw865530co0iN1yAhhv9c3Yu7zmrHgnWZDJyUzktf/EB+fkA3+iWEldpt9YsLzSLxrCWVXbRsN5CZ2TBgWKtWrcauW7fO3+GIeMXm3Ye48+0VLNqwhx4t6jJxZGdapcT5OywJQaV1W5U7eQQzjXlIqHHO8caSDB54fzXZR/O48TetGJeWSlSEVuEV76nKmIeIBCAz48IeTZk3IY1BHRowae5aznlsAUt/2Ofv0KSGUPIQCWLJ8dE8fukpPHNFdw7mHGPUk4u4+92VHDxyzN+hSYgrNXmY2Wwzu9XM2lZXQCJScWe2a8CcCWlc0bsFL37xAwMnpTPvO62LJb5TVsvjCmAfcI+ZfW1mU81shJlpdE4kwMRFR3DP8A68dd1pJNaK5HcvLuGG/37NrgMq6xXvq0i1VRjQCxgKnAlkA3PKM9PcX1RtJTVVbl4+0+ZvZPJH64iJCOPOs9ppcqFUmE+qrcysPjDYOfffqgRXHVRtJTXVxsyD3PHWCr7ctJdeJ9Zj4shOtExWx4GUj0+qrZxzu4MhcYjUZC2T43hl7KlMHNmJ77bvZ8jkBUz5ZD25efn+Dk2CnKqtREJcWJhxSc9mfDQhjQHtUnho9hqGPbaQb35UWa9UnpKHSA2RkhDDE5d1Y9robvx0OJeRUxfxt5mrOKSyXqmESicPM7vKm4GISPUY1KEhcyf0Z/SpzXl+0WYGPTyfT77f5e+wJMhUpeXxN69FISLVKj4mkntHdGT6uN7UjgrnqucXc9Mr37D74BF/hyZBoqz9PL4t6RRwknOuslvRVguV6oqU7eixfKZ+uoEpn6ynVlQ4d53djgu6NVFZr1S+VNfMdgKD8UwU/MUpYJFz7gSvRelDKtUVKdv6XQe4460VLN68j9NSk/j7eZ1oUT/W32GJH1WlVHcWEOec++G4r83Ap16OU0T8qFVKPK9d05sHzuvIiowsBj8yn6mfblBZrxRLS7KLyK/syMrh7hkrmb1qJ+0aJfDg+Z3o3KSOv8OSaqYl2UWkQhomxvDU6O48+dtu7Dl4hHOnfMZ9s77j8FGV9YpHWavqfl3WDcpzjYgEpyEdGzJ3QhoX92zGMws3MXDSfNLXZvo7LAkAZQ2YZ+PZr7zES4BE51wzbwfmTeq2Eqm6rzbt5Y63vmVD5iHOPfkE/nJOe5LiArrgUqqoKtVWzctx/zznXEZlg/MlleqKeNeRY3lM+WQDUz9dT1x0BH85pz3ndW2sst4QpT3M1fIQ8aq1Ow9w+5vf8vWPP9GvdX3+fl4nmtar7e+wxMs0YC4iXnVSg3imjzuNe0d04Osf9jHw4XSmzd/AMZX11hhKHiJSKWFhxuW9WzB3Qhp9W9Xn7+9/z7lPfMbKrVn+Dk2qQbmSh5m1L+bY6V6PRkSCzgl1avH05d2Zcukp7Mg6wogpnzHxg9VkH83zd2jiQ+VtebxuZn8yj1pm9hgw0ZeBiUjwMDPO7tyIjyakMeqUJjyVvpHBj8xn4brd/g5NfKS8yaMX0BRYBCwGtgF9fBWUiASnxNqRPDiqM6+MPZXwMOO3z3zJba8vZ9+ho/4OTbysvMkjF8gGagExwCbnnEbGRKRYvVOT+OCWftxwRirvLtvKgEnpvLtsKzWhurOmKG/yWIwnefQA+gKXmNl0n0UlIkEvJjKcPwxuy8yb+tKkXm1ueXUZVz63mIx9h/0dmnhBeZPHGOfcX51zuc65Hc65EcC7vgxMREJDu0YJvHXdafz1nPYs3ryXQQ/P55mFm8jLVyskmJVrkqCZFbv8iHPuR69H5EWaYS4SWDL2HebP76zk0zWZdGmSyMSRnWl/QoK/w5ISVHmGuZmtAByetaxigBOBNc65Dt4M1Fc0w1wkcDjnmLF8G/fO/I6s7FzG9m/JLWe2JiYy3N+hyXGqPMPcOdfJOde54L+tgZ7AQm8GKSI1g5kx4uTGfHRbGud2bczUTzcw5JH5LFqvst5gUqkZ5s65r/EMnouIVEqd2lH83wVd+O/veuGAS//9JX+cvpyfDqusNxhElOciM5tQ5GUYcAqgRf1FpMr6tKrP7PH9eWTeOp5esJGPv9/FPcM7cHanRlqtN4CVt+URX+QrGngPGOGroESkZomJDOf2oW2ZcWMfGiXW4saXv+F3Lyxh60/Z/g5NSqAl2UUkoBzLy+f5RZv511gnWukAAA3VSURBVJy1hBn8YXAbRvduQXiYWiHVrSqbQc3EU2VVLOfc8KqH53tKHiLBZ8vew9z1zkrmr83k5KZ1ePD8zrRpGO/vsGqUqiSPtNJu7JxLr2Js1ULJQyQ4Oed4Z9lW7pu1mv3ZuYxLS+XG37RSWW81KS15lDVgvinQJwKKSOgyM87r2oS0k1K4f9Z3PP7Jet5fuZ2J53WiV8skf4dXo5U1YP5O4R/M7E0fxyIiUqx6sVFMuuhkXry6J7l5+Vw07QvueGsFWdm5/g6txioreRQdoWrpy0BERMrS/6RkZo/vzzX9W/La4h8ZOCmdD1du93dYNVJZycOV8OegYGbDzGxaVpa2xRQJFbWjIrjzrHa8e0NfkuOjGffS11zz4hJ2ZOX4O7QapawB8zzgEJ4WSC2gcC1lA5xzLihWNNOAuUhoOpaXzzMLN/HwvLVEhoXxx6FtuaxnM8JU1usVlV7byjkX7pxLcM7FO+ciCv5c+DooEoeIhK6I8DCuTUtlzvg0ujStw1/eWcmFT33Oup0H/B1ayKvU2lYiIoGkWVJt/jOmJ/+6oAvrMw9y1qMLeHjuWo4cy/N3aCFLyUNEQoKZcX63JsybkMZZnRox+aN1nP3oQpZs3uvv0EKSkoeIhJT6cdFMvrgrz13Vg+yjeYx68nP+/M4K9ueorNeblDxEJCSd0SaFObf25+o+J/Lyl56y3jmrdvg7rJCh5CEiISs2OoK/DmvP29f3oW7tKK75z1Kue2kpO/errLeqlDxEJOR1aVqHmTf15Y9D2vDx97sYMCmdl7/8kfz8oJu+FjCUPESkRogMD+P601vx4fj+dDwhkTvfXsHF075g/a6D/g4tKCl5iEiNcmL9WF4e24t/nt+ZNTsPcNbkBTz60TqOHsv3d2hBRclDRGocM+PCHk2ZNyGNQR0aMGnuWs55bAFLf9jn79CChpKHiNRYyfHRPH7pKTxzRXcO5Bxj1JOLuPvdlRw8cszfoQU8JQ8RqfHObNeAuRPSuKJ3C1784gcGTkrno9U7/R1WQFPyEBEB4qIjuGd4B9687jTiYyIY88ISbnj5azIPHPF3aAFJyUNEpIhTmtVl1k39uG3gScxdtZMBk9J5ffEWSluBvCZS8hAROU5URBg3ndmaD8b3o03DeP745rdc8vQXbNp9yN+hBYyQTh7aDEpEqiI1OY5Xx57K38/rxKpt+xn8yHymfLKe3DyV9Za6GVSo0GZQIlJVu/bncPeMVXywcgdtG8bzj/M7c3LTOv4Oy6cqvRmUiIh4pCTEMPW33XhqdDf2HT7KyCc+428zV3Gohpb1KnmIiFTA4A4NmTshjUt7NeO5zzYz6OH5fLJml7/DqnZKHiIiFZQQE8n953Zi+rje1IoK56rnFnPzK9+w+2DNKetV8hARqaTuLerx3s19GT+gNR+s3M6ASelMX5pRI8p6lTxERKogOiKc8QNO4oNb+tEqOY7fv7Gc0c98xQ97QrusV8lDRMQLWqXE8/q1vbn/3I4s3/ITgx+Zz5PpGzgWomW9Sh4iIl4SFmb89tTmzJ2QRv/Wyfzjg+8ZMeUzVm4NvblmSh4iIl7WMDGGaZd3Z+plp7DrwBGGP76QB977jsNHQ6esV8lDRMRHhnZqxLwJaVzUoxlPL9jEoIfnM39tpr/D8golDxERH0qsFcnEkZ147ZpTiYoI4/Jnv+LW15ax99BRf4dWJUoeIiLVoFfLJN6/uR83/aYVM5dv48x/fcrb3wRvWa+Sh4hINYmJDOe2QW147+Z+tKgfy62vLeeK5xazZe9hf4dWYUoeIiLVrE3DeKaPO42/De/A0s17GfTwfP69YGNQlfUqeYiI+EF4mHHFaS2YOyGNPq2SuP+91Zz3xCJWbQuOsl4lDxERPzqhTi2evrw7Uy49he1ZOQx//DMmfrCa7KN5/g6tVEoeIiJ+Zmac3bkRH01IY9QpTXgqfSNDJs/ns/W7/R1aiZQ8REQCRGLtSB4c1ZmXx/bCgMv+/SW3vb6cfQFY1qvkISISYE5Lrc+H4/tz/empvLtsKwMmpfPusq0BVdar5CEiEoBiIsP545C2zLypL03q1uKWV5dx1fOLydgXGGW9Sh4iIgGsXaME3rq+D385pz1fbfKU9T67cBN5+f5thSh5iIgEuPAwY0zfE5lza396nliPe2d9x8ipi1i9fb/fYlLyEBEJEk3q1ua5K3sw+eKTydh7mGGPLeSh2d+Tk1v9Zb1KHiIiQcTMGHFyY+ZNSOPcro2Z8skGhk5ewOcb9lRrHEoeIiJBqG5sFP93QRdeGtOLvHzHJU9/wZ+mf0vW4dxqeX7AJw8za2lmz5jZ9OOOx5rZUjM7x1+xiYj4W9/W9Zk9vj/X9m/J9K8zOHNSOrO+3ebzsl6fJg8ze9bMdpnZyuOODzGzNWa23sxuL+0ezrmNzrkxxZz6E/C6N+MVEQlGtaLCueOsdrx7Qx8aJcZw48vf8LsXlrDtp2yfPdPXLY/ngSFFD5hZODAFGAq0By4xs/Zm1snMZh33lVLcTc1sAPAdsNO34YuIBI+OjRN5+/rTuOusdizasIeBk9J58fPNPnlWhE/uWsA5N9/MWhx3uCew3jm3EcDMXgVGOOcmAuXtgjoDiMWTfLLN7H3n3C/WMjaza4BrAJo1a1bp70FEJJhEhIcxtn9LhnRsyJ1vr2DNjgO+eY5P7lq6xsCWIq8zgF4lXWxmScADQFczu8M5N9E5d1fBuSuB3ccnDgDn3DRgGkD37t0DZ06/iEg1aFqvNi9e3ZOjPtojxB/Jw4o5VuKHu3NuDzCuhHPPeykmEZGQY2ZER4T75N7+qLbKAJoWed0E2OaHOEREpJL8kTwWA63N7EQziwIuBmb4IQ4REakkX5fqvgJ8DrQxswwzG+OcOwbcCMwGVgOvO+dW+ej5w8xsWlZWcGzrKCISLCyQ1of3le7du7slS5b4OwwRkaBiZkudc92LOxfwM8xFRCTwKHmIiEiFKXmIiEiF+WOeR7Uxs2HAMGC/ma0DEoHKjJ7XB3Z7MzYpUWX/jgJdoH5f/ojL18/0xf29cc+q3sMfn1/NSzpRIwbMC5nZNOfcNZV435KSBo3Euyr7dxToAvX78kdcvn6mL+7vjXtW9R6B9vlV07qtZvo7AClTqP4dBer35Y+4fP1MX9zfG/es6j0C6t9QjWp5VJZaHiISrNTy8K9p/g5ARKSSfPL5pZaHiIhUmFoeIiJSYUoeIiJSYUoeIiJSYUoeFWRmsWb2gpk9bWaX+TseEZGKMLOWZvaMmU2vyn2UPAAze9bMdpnZyuOODzGzNWa23sxuLzg8EpjunBsLDK/2YEVEjlORzzDn3Ebn3JiqPlPJw+N5YEjRA2YWDkwBhgLtgUvMrD2enQ8L92DPq8YYRURK8jzl/wzzCiUPwDk3H9h73OGewPqCLH0UeBUYgWcb3SYF1+j/n4j4XQU/w7xCH34la8z/WhjgSRqNgbeA881sKgG2XICISBHFfoaZWZKZPQl0NbM7KnvzkF5Vt4qsmGPOOXcIuKq6gxERqaCSPsP2AOOqenO1PEqWATQt8roJsM1PsYiIVJRPP8OUPEq2GGhtZieaWRRwMTDDzzGJiJSXTz/DlDwAM3sF+BxoY2YZZjbGOXcMuBGYDawGXnfOrfJnnCIixfHHZ5gWRhQRkQpTy0NERCpMyUNERCpMyUNERCpMyUNERCpMyUNERCpMyUNERCpMyUNqPDPLM7NlRb5uL/tdvmdmm81shZl1N7O3C2Jbb2ZZRWI9rYT3/s7M/nPcsQYFy3ZHmtlrZrbXzM6tnu9GQo3meUiNZ2YHnXNxXr5nRMEkrarcYzPQ3Tm3u8ix04HfO+fOKeO9dYF1QBPnXE7BsRuBTs65awtev4Rnb5p3qhKn1ExqeYiUoOA3/7+Z2dcFLYC2BcdjCzbfWWxm35jZiILjV5rZG2Y2E5hjZmFm9oSZrTKzWWb2vpmNMrMzzeztIs8ZaGZvVSHOHmaWbmZLzewDM2vgnNsHLALOLnLpxcArlX2OSFFKHiJQ67huq4uKnNvtnDsFmAr8vuDYXcDHzrkewBnAQ2YWW3CuN3CFc+43eHadbAF0An5XcA7gY6CdmSUXvL4KeK4ygZtZNDAZON851w14Cbiv4PQreBIGZta0IJb5lXmOyPG0JLsIZDvnTi7hXGGLYCmeZAAwCBhuZoXJJAZoVvDnuc65wk15+gJvOOfygR1m9gl41sQuGI/4rZk9hyepXF7J2NsBHYB5ZgYQjmc1VfAsgveomcUBF+FZ2yi/ks8R+QUlD5HSHSn4bx7/+3kxPL/pryl6oZn1Ag4VPVTKfZ/Ds5lYDp4EU9nxEQO+dc71O/6Ec+6Qmc3Ds3vcxcB1lXyGyK+o20qk4mYDN1nBr/pm1rWE6xbi2XUyzMwaAKcXnnDObcOzt8Kf8ew/XVnf4dkdrmdBLFFm1qHI+VeAPwB1nHOLq/AckV9Q8hD59ZjHP8q4/j4gEvjWzFbyvzGG472JpwtpJfAU8CWQVeT8f4EtzrnvKhu4c+4IMAqYZGbLgW+AXkUu+RBPl9qrlX2GSHFUqiviQ2YW55w7aGZJwFdAH+fcjoJzjwPfOOeeKeG9mzmuVNfLsalUVypNLQ8R35plZsuABcB9RRLHUqAznuqokmQCH5lZd28HZWavAX3wjLmIVJhaHiIiUmFqeYiISIUpeYiISIUpeYiISIUpeYiISIUpeYiISIUpeYiISIX9P9HYLlnmGBHCAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "my_custom_model.plot(energy_range=[1, 10] * u.TeV)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a next step we can also register the custom model in the `SPECTRAL_MODELS` registry, so that it becomes available for serilisation:" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [], "source": [ "SPECTRAL_MODELS.append(MyCustomSpectralModel)" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [], "source": [ "model = SkyModel(spectral_model=my_custom_model, name=\"my-source\")\n", "models = Models([model])\n", "models.write(\"my-custom-models.yaml\", overwrite=True)" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "components:\r\n", "- name: my-source\r\n", " type: SkyModel\r\n", " spectral:\r\n", " type: MyCustomSpectralModel\r\n", " parameters:\r\n", " - {name: amplitude, value: 1.0e-12, unit: cm-2 s-1 TeV-1, min: 0.0, max: .nan,\r\n", " frozen: false}\r\n", " - {name: index, value: 2.0, unit: '', min: 0.0, max: .nan, frozen: false}\r\n", " - {name: reference, value: 1.0, unit: TeV, min: .nan, max: .nan, frozen: true}\r\n", " - {name: mean, value: 3.0, unit: TeV, min: 0.0, max: .nan, frozen: false}\r\n", " - {name: width, value: 0.1, unit: TeV, min: 0.0, max: .nan, frozen: true}\r\n" ] } ], "source": [ "!cat my-custom-models.yaml" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly you can also create custom spatial models and add them to the `SPATIAL_MODELS` registry. In that case gammapy assumes that the evaluate function return a normalized quantity in \"sr-1\" such as the model integral over the whole sky is one." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.0" }, "nbsphinx": { "orphan": true } }, "nbformat": 4, "nbformat_minor": 4 }