{
 "cells": [
  {
   "cell_type": "raw",
   "id": "eaf29ef7",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/restructuredtext",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    ".. _cookbook-iso:\n",
    "\n",
    "Isochrone Color-Magnitude Diagram\n",
    "=================================\n",
    "\n",
    "This :ref:`cookbook` recipe demonstrates how to create a color-magnitude diagram (CMD) for an isochrone, using the\n",
    "MSG Python interface.\n",
    "\n",
    "Preparation\n",
    "-----------\n",
    "\n",
    "Download :download:`isochrone-cmd.tar.bz2 <../downloads/isochrone-cmd.tar.bz2>` and unpack it into your working directory using the :command:`tar` utility. This archive contains the following items:\n",
    "\n",
    "* :file:`MIST.iso` --- isochrone data file\\ [#isochrone]_\n",
    "* :file:`read_mist_models.py` --- Python module for reading :file:`MIST.iso`\n",
    "\n",
    "Initialization\n",
    "--------------\n",
    "\n",
    "To start, import modules and set other parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2eb4244c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import standard modules\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy.constants as con\n",
    "\n",
    "# Import pymsg and read_mist_models\n",
    "\n",
    "import pymsg\n",
    "import read_mist_models\n",
    "\n",
    "# Set constants\n",
    "\n",
    "R_SUN = con.R_sun.cgs.value\n",
    "PC = con.pc.cgs.value\n",
    "\n",
    "# Set plot parameters\n",
    "\n",
    "plt.rcParams.update({'font.size': 12})"
   ]
  },
  {
   "cell_type": "raw",
   "id": "07a01f15",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "Create PhotGrid Objects\n",
    "-----------------------\n",
    "\n",
    "Next, create a pair of :py:class:`pymsg.PhotGrid` objects for interpolating fluxes in the Johnson B and V bands:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9c7a8544",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create PhotGrid objects using the demo specgrid\n",
    "\n",
    "MSG_DIR = os.environ['MSG_DIR']\n",
    "GRID_DIR = os.path.join(MSG_DIR, 'data', 'grids')\n",
    "PASS_DIR = os.path.join(MSG_DIR, 'data', 'passbands')\n",
    "\n",
    "photgrid_B = pymsg.PhotGrid(os.path.join(GRID_DIR, 'sg-demo.h5'),\n",
    "                            os.path.join(PASS_DIR, 'pb-Generic-Johnson.B-Vega.h5'))\n",
    "photgrid_V = pymsg.PhotGrid(os.path.join(GRID_DIR, 'sg-demo.h5'),\n",
    "                            os.path.join(PASS_DIR, 'pb-Generic-Johnson.V-Vega.h5'))"
   ]
  },
  {
   "cell_type": "raw",
   "id": "5b4d9d5f",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "Read Isochrone Data\n",
    "-------------------\n",
    "\n",
    "Now read in the isochrone data file and extract the stellar parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2d008331-53c4-468b-983a-8682d3ab0882",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading in: MIST.iso\n"
     ]
    }
   ],
   "source": [
    "# Read isochrone data file\n",
    "\n",
    "iso = read_mist_models.ISO('MIST.iso')\n",
    "\n",
    "# Extract stellar parameters\n",
    "\n",
    "Teff = 10**iso.isos[0]['log_Teff']\n",
    "logg = iso.isos[0]['log_g']\n",
    "\n",
    "R = 10**iso.isos[0]['log_R']*R_SUN"
   ]
  },
  {
   "cell_type": "raw",
   "id": "7e98770e-7c62-43a0-9ea8-27b1d3e50736",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "Evaluate Irradiances\n",
    "--------------------\n",
    "\n",
    "Using these parameters, evaluate the photometric irradiance in each band:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "163c7134-4bf4-433c-9e89-af740f708f45",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Set the distance to the standard 10 parsecs used to define\n",
    "# absolute magnitudes\n",
    "\n",
    "d = 10*PC\n",
    "\n",
    "# Evaluate irradiances\n",
    "\n",
    "n = len(Teff)\n",
    "\n",
    "F_obs_B = np.empty(n)\n",
    "F_obs_V = np.empty(n)\n",
    "\n",
    "for i in range(n):\n",
    "\n",
    "    # Set up photospheric parameters dict\n",
    "\n",
    "    x = {'Teff': Teff[i],\n",
    "         'log(g)': logg[i]}\n",
    "    \n",
    "    # Evaluate irradiances. Use try/execpt clause to deal with \n",
    "    # points that fall outside the grid\n",
    "\n",
    "    try:\n",
    "        F_obs_B[i] = (R[i]/d)**2 * photgrid_B.flux(x)\n",
    "        F_obs_V[i] = (R[i]/d)**2 * photgrid_V.flux(x)\n",
    "    except (ValueError, LookupError):\n",
    "        F_obs_B[i] = np.NAN\n",
    "        F_obs_V[i] = np.NAN"
   ]
  },
  {
   "cell_type": "raw",
   "id": "fc36f462-c4f3-4c76-ab16-ae76a67b6b87",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "\n",
    "Plot the CMD\n",
    "------------\n",
    "\n",
    "As a final step, convert fluxes to absolute magnitudes and plot the CMD:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d30546b7-42d4-4d3d-a1c6-4a1ec93676df",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAGzCAYAAAAPGELKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABNYUlEQVR4nO3dd3xT5f4H8E+S7pVuulu6KC0tpWUUVJC9FBAHoOJAvQ4cyPWnICrFBd7rQq+ggOIAVBBwFmQUZCO0VFpmoaUbSlfSmTbJ+f1RGym0paTjZHzer1de3p6ck3yfmz7Nh+c85zkSQRAEEBEREdENkYpdABEREZExYogiIiIi0gNDFBEREZEeGKKIiIiI9MAQRURERKQHhigiIiIiPTBEEREREenBQuwCDI1Wq0VhYSEcHR0hkUjELoeIiIjaQRAEVFZWwsfHB1Jp94wRMURdpbCwEP7+/mKXQURERHrIy8uDn59ft7wXQ9RVHB0dATR+CE5OTiJXQ0RERO2hVCrh7++v+x7vDgxRV2k6hefk5MQQRUREZGS6cyoOJ5YTERER6YEhioiIiEgPDFFEREREemCIIiIiItKDSYSo5ORkzJo1CxEREbC3t4evry8mT56MlJQUsUsjIiIiE2USIWr58uW4cOECnnvuOSQlJWHp0qUoLi5GQkICkpOTxS6PiIiITJBEEARB7CI6qri4GJ6ens22VVVVITQ0FH369MGOHTva/VpKpRJyuRwKhYJLHBARERkJMb6/TWIk6uoABQAODg6IjIxEXl6eCBURERGRqTPZxTYVCgVSU1MxYsSINvdTqVRQqVS6n5VKZVeXRkRERCbAJEaiWjJ79mxUV1djwYIFbe63ePFiyOVy3YP3zSMiIqL2MLgQtXv3bkgkknY90tLSWnyNV199FWvXrsUHH3yA+Pj4Nt9v/vz5UCgUugdP/xEREVF7GNzpvF69emHlypXt2jcgIOCabYsWLcKbb76Jt956C08//fR1X8Pa2hrW1tY3XCcRERGZN5O4Oq/JokWLkJiYiMTERCxcuFCv1+DVeURERMaHV+d1wBtvvIHExES88soregcoIiIiovYyuNN5+njvvffw2muvYdy4cZg4cSIOHTrU7PmEhASRKiMiIiJTZRIh6pdffgEAbN26FVu3br3meRM6Y0lEREQGwiRC1O7du8UugYiIiMyMycyJIiIiIupODFFEREREemCIIiIiItIDQxQRERGRHhiiiIiIiPTAEEVERESkB4YoIiIiIj0wRBERERHpgSGKiIiISA8MUURERER6MInbvhAREQmCgNLqeuSX10JZ24Caeg1Uag2cbCzham+F8B6OsLWSiV0mmRCGKCIiMjqCICC3rAYHzpfieH4FThVVIvNSJarrNa0eYyWTYt74CMy6uWc3VkqmjCGKiIiMyjtbT+PntEIUVNS2uZ+znSVc7KyQXVINAKjXaPHBjrMMUdRpGKKIiMioXK5UoaCiFpYyCWL9naGobcDZS1XX7FdR04CKmoZm26YP8O+uMskMMEQREZFRaDqF52JniTBPB8htLXFRWYe8srZHpK709IiwLqyQzA1DFBERGYVpnx3CnxfK9DrWw9Eavzx9M+S2lp1cFZkzLnFARERGoae7PaxkUtjpcYXdxGhveMltuqAqMmcMUUREZPAulFTD3toCvi62qGnjCryWBLvb419Dg7uoMjJnPJ1HREQGqaCiFknHi/DzX4VIL1C0up+tpQwqtQZaofn2Pr5OuKe/P+7p7w8bS64PRZ2PIYqIiAyCIAg4UajEtpOXsOPkJZwsUra6r52VDKGeDlDWNqBIUacLUHJbS9zRzxd39/dDlI+8myonc8UQRUREolGpNTh4vhQ7Tl3CjpPFuKisa3VfqQS4KdQdznZWKKyoRUpOue65Xj0c8fBNQZgc68tVyanbMEQREVG3Kq+ux64zxdh+8hL2nL3c5irjANDb2wljIntApdZiS0YR9maWAGgMVaMje+ChIT2REOwKiUTSHeUT6TBEERFRl7tQUo0dpy5h28lLOHqh7Jr5S1cL83TA6MgeiPFzxqGsUny+LxtVKjUAwMnGAtMHBmBmQiD8Xe26oXqiljFEERFRp9NqBaTlV2D7yUvYfvISzhVfu6L4lSQSID7ABWOiemBEhCeyS2qw7nAOlv9xHsLfgSvEwx4P3dQTd8b5ws6KX18kPv4WEhFRpxAEAcfyKvDb8SIkpRehSNH6/CYAsLaQ4pYwd4yO7IGRvXugtl6D9UfzcN+qw7ikVOn2G97LAw/d1BO3hLpDKuUpOzIcDFFERKQ3QRDwV74Cvx0vRFL6xXbdFHhEhCfGRHphaLg7pBIJkk8X49/r/8KezMu6USdXeyvcHe+HaQP8Eezh0A0tIbpxDFFERHRDBEFARoESv6YX4rfjRcgvbzs4+bnYYkykF8ZE9UD/QBfIpBIczSnHG7+ewm/HC6GsU+v2vTnUHTMGBmB0ZA9YWXA9aDJsDFFERHRdTWs4/ZZehN+OFyG3rKbN/SO8HDG+jzfGRPVAhJcjJBIJskuq8VHyOWw+lt/spsHechvc0c8X0wb4I9DNvqubQtRpGKKIiKhFgiDgVFElfvt7xOlCadvBKcTDHrf39cFtMT4I9Ww8BVelUmPD0XysP5qHo1es62RvJcP4aG9MjfNFQk83znUio8QQRUREOoIg4MylSvx2vHHEKaukus39A1ztcFuMN27v66MbcRIEAX9ml2H90TwkpRfp7nUnlQC3hHlgapwvxkR6cVFMMnoMUUREhLOXKvHr8SL8drwQ5y+3HZy85TaYGN0YnGL85LpFLi8q6rAxNR8/pOQj+4rwFexhj3v6+2NqP194Otl0aTuIuhNDFBGRmTpXXNU44pReiLOX2l7Hyd3BChOjvXFbXx/EB7joTr+pNVrsPH0J3/2Ziz/OXtYtomlvJcNtMT64Z4Af4gJcuJo4mSSGKCIiM5J1uSk4FeH0xco293W2s8T4Pl64PcYHg4LdILti3tLlShW+P5KLdYdzUXjFelADg1xxd38/TIj2hr01v2LItPE3nIjIxOWV1eDnvwrx6/EinCpStrmvvZUMY6O8cHusD24OdYel7J9lBgRBQEpOOb4+mIMtGUVo0DQOO7naW+Ge/v64p78f13Qis8IQRURkgmrrNdh6oggbjubjwPnSNve1lEkwLNwTU/r5YFTvHrCxbD7hu6ZejR+PFeKbQznNQli/AGc8MDgQE6K9YW3BSeJkfhiiiIhMhCAISM2twA8pefj1ryJUqtRt7j+wpysmx/pgQh9vuNhbXfP8+ctVWHMoBz+k5KPy7wUxbSylmNzXFzMHB6KPr7xL2kFkLBiiiIiM3CXlP1fFZV3nyroIL0dMjvXFpFgf+DrbXvN840TxYnxzMAf7zpXotge62WFmQiDujveH3M6y09tAZIxMJkRVVVXhlVdewfr161FWVoaIiAjMmzcP06dPF7s0IqJOp1JrsPNUMTYczWt2VVxLfOQ2mBTriyn9fBDh5dTiPiVVKnx/JA9rD+XoJopLJMDICE/cnxCIoWEeXBCT6ComE6KmTp2KI0eOYMmSJQgPD8e6deswY8YMaLVa3HvvvWKXR0TUKTIKFPghJR8/phWgoqah1f3ktpaYGOONyX19MCDItdUAlJZXgS/3Z+O39H8mirvYWWLagADcNygA/q52XdIOIlMgEQShjX+/GIekpCRMnDhRF5yajBkzBidOnEBubi5ksvZNelQqlZDL5VAoFHByavlfbERE3am8uh6bjxVgQ0p+m1fXWVtIMSqyB6bE+mJYuEerN/Bt0GiRlF6ELw9cwLHcCt32WP9/JopfPbmcyNCJ8f1tEiNRmzdvhoODA+6+++5m2x9++GHce++9OHz4MIYMGSJSdUREN06jFbDvXAnWH8nD9pOXUK/Rtrn/kqnRmBjjDUeb1ucrlVapsO5wLtYczsElpQoAYCWT4ra+3nhoSBBi/Jw7swlEJs8kQlRGRgZ69+4NC4vmzYmJidE931qIUqlUUKlUup+VyrbXUCEi6kp5ZTXYcDQPP6TkN1vEsi3/N7YXpg8MaPX5jAIFvjxwAT//VYh6dWMY83C0xv2DAnHvoAB4OFp3Su1E5sYkQlRpaSmCg4Ov2e7q6qp7vjWLFy/GokWLuqw2IqLrqWvQ4PcTF/H9kbzrrul0pWB3e3z/+OAWQ5Bao8W2k5fw5f4L+PNCmW57Xz85Hr6pJyZEe7d6uo+I2sckQhSANu/L1NZz8+fPx9y5c3U/K5VK+Pv7d2ptREQtyShQ4PsjefgprQDKurbXdLrS7X19sGBCb3jJr72Zb7VKjW//zMXq/RdQUFELALCQSjA+2hsP3xSEfv7OvI8dUScxiRDl5ubW4mhTWVnjv76aRqRaYm1tDWtrDmUTUfeoqKnHj8cKsP5oPk5e5xYsV3t5QgTuHRQIhxbuSVdRU4+vDuRg9YFs3VV7rvZWuHdgAO5PCGwxcBFRx5hEiIqOjsa3334LtVrdbF5Ueno6AKBPnz5ilUZEBK1WwP7zJVh/NB+/n7iom5fUHu4O1lh4eyTG9/GCheza02/Fyjp8vi8baw7loLpeAwAIcrPDE8NCMKWfL6+yI+pCJhGi7rjjDqxcuRIbN27EtGnTdNu/+uor+Pj4YNCgQSJWR0TmqqCiFuuPNE4Sbzq11l7De3ngiWEhGNjTtcXTb4UVtVi2+xzWH83XhbIIL0fMHh6KCdHekHFhTKIuZxIhavz48Rg9ejSefPJJKJVKhIaG4ttvv8XWrVuxZs2adq8RRUTUURqtgN1nirHucC52nSlucyXxlswYGIBHbu6JUE+HFp8vrqzDsl3nse5wrm7Zg/hAF8weHoLhvTw534moG5lEiAKATZs2YcGCBXjttdd0t3359ttvedsXIuoWFxV1+P5IHr4/kttsaQI/F1sUKeqgaSNNyW0t8fBNQZiZEAg3h5bnaJZX1+PTPefx1YELqGtoDE+Derri+dHhGNTKaBURdS2TWLG8M3HFciJqL41WwJ7My1h3OBfJp4t1QcnZzhJBbvZIy6to8/geTtZ4dmQY7ozza3XukrKuAav2ZuOLfdmoUjVewdcvwBkvjOmFISFuDE9Ef+OK5URERiCjQIHbPt53zfb4QBdU1jXg7KUqpNVUtHp8Xz85nh4RhpERnq3e066mXo0vD1zAZ39kQVHbeLVdpLcTXhgbztN2RAaCIYqIqB00WgGz16Zi64mL1zx3Z5wftmQUISWnvM3XGBzshmdGhGJwGyNIDRotvj+Shw93ZKKkqvFuCmGeDpg7Ohxjo7xaDV1E1P0YooiI2nA8vwKT/re/1efd7K2wMTW/zdcYEeGJ2cNDER/o0uo+giBga8ZF/Pf3M8gqqQYABLjaYe7ocNze14dX2xEZIIYoIqKrqNQavLI5AxtSWg5H46K8cEecLx7/JgWl1fWtvs7EaG88NTwEUT7yNt/vUFYpFm85jb/+nkPlZm+F50aFYfqAAN6ahciAMUQREV1h/7kS3LfqcIvPvTKxN3p5OeLdbWfx+Dcprb7G1DhfPHVraKvLFDQ5fVGJ/2w9g+TTxQAAOysZHrslGI8NDW5xVXIiMizspUREAOrVWizffR4f7Djb6j5v/naqzde4b1AAnhgWAn9Xuzb3K6ioxQfbz2Jjaj4EofHedjMGBuCZkaHwdOTtWYiMBUMUEZm9tLwK/N+Gv5BZXHXdfa0spM1u22JrKcO9gwLw2C3B170/nbKuAZ8kn8PqAxd0rzEx2hsvjO2Fnu72HWsEEXU7higiMmsZBQrctfwA1NdZWtxKJsXd/f3w9IhQvLwpHam5FZiZEIiHbwpqdYHMJhqtgO+P5OG9bWd0c6gG9XTF/Am9Eevv3FlNIaJuxhBFRGZNpda0GaA2PTUEcQHNr6pb/fDAdr/+wfOleP3XkzhVpAQAhHjYY8HE3lzricgEMEQRkVk7X1zd4vabQt2waFIUQj0d9Xrd3NIavJ10SreulJONBeaMCsfMwYGwlPGKOyJTwBBFRGZt3Z+5zX5e9UB/3NrLAxZ6Bp0qlRr/Sz6HL/Zlo16jhUwqwX2DAjBnVDhc7a06o2QiMhAMUURk1pZOj8W+cyW4va8PnGws9X4dQRCQlH4Ri345geLKxpXGbwlzx6u3RSK8h36jWURk2BiiiMisBbrZI9CtY1fG5ZXV4NWfMrD7zGUAQJCbHV6ZGImRvTnviciUMUQREempQaPF/E3p+OHvlc2tZFI8NTwET94aAmsLmcjVEVFXY4giItJDam45pi47oPvZz8UWX80aiBCPtlcpJyLTwRBFRHSDyqvrmwUoANj9wq16T0YnIuPEHk9EdAMOnCvB+KV7dT8PDHJF1tsTGKCIzBBHooiI2qFercX728/isz3nIQhAsIc9PpreD3185WKXRkQiYYgiIrqOnNJqPPvtMfyVrwAAzBgYgFdv6w07K/4JJTJn/AtARNSGH48V4JUfM1ClUkNua4l37ozGuD7eYpdFRAaAIYqIqAVVKjVe+ykDm1ILADTOffpgeix8nW1FroyIDAVDFBHRVU4VKfHU2lRkl1RDKgGeHRmGp4eHcvI4ETXDEEVEdIVNqfl4eXM66hq08JHbYOmMfhgQ5Cp2WURkgBiiiIgAqNQavP7LSaw93HhD4mHhHvhwWixceNNgImoFQxQRmb388hrMXpuKv/IVkEiAOSPD8cyIUEilvO8dEbWOIYqIzNofZy/jue+OoaKmAc52lvhwWixu7eUpdllEZAQYoojILAmCgGW7z+PdbWcgCECMnxzL7ouDn4ud2KURkZFgiCIis1PXoMG8jcfxY1ohAODeQQFYeHskrC1kIldGRMaEIYqIzEqxsg6PfZOCv/IqIJNKsGhSFO5PCBS7LCIyQgxRRGQ20vMVeOzro7iorIPc1hLL74vDkFB3scsiIiPFEEVEZuG340X494Y01DVoEerpgFUP9EeQu73YZRGREWOIIiKT1jSB/L+/nwEA3NrLAx/N6AcnG0uRKyMiY8cQRUQmS6MVsPDnDKw51LiA5iM398TLE3pDxvWfiKgTMEQRkUmqa9DgmW+PYfvJS5BIgNdui8TDN/UUuywiMiEMUURkcirrGvDoV0dxOLsMVhZSLJ0Wi/HR3mKXRUQmhiGKiExKWXU9Hlr9J47nK+BgbYFVD/ZHQrCb2GURkQmSil1AZ0lOTsasWbMQEREBe3t7+Pr6YvLkyUhJSRG7NCLqJhcVdZj22UEcz1fA1d4K3z6WwABFRF3GZELU8uXLceHCBTz33HNISkrC0qVLUVxcjISEBCQnJ4tdHhF1sfzyGtz92QFkFlfBy8kG6x9PQLSfXOyyiMiESQRBEMQuojMUFxfD07P5TUOrqqoQGhqKPn36YMeOHe16HaVSCblcDoVCAScnp64olYg6WZGiFtM+O4TcshoEudnhm0cGwd+V98AjMidifH+bzJyoqwMUADg4OCAyMhJ5eXkiVERE3aFYWYf7Vh5GblkNAt3s8N2/BsNLbiN2WURkBkwmRLVEoVAgNTUVI0aMaHUflUoFlUql+1mpVHZHaUTUCUqrVLhv1WFklVTD19kW6x5LYIAiom5jMnOiWjJ79mxUV1djwYIFre6zePFiyOVy3cPf378bKyQifVXWNWDm53/q5kB9+1gCfJ1txS6LiMyIyYaoV199FWvXrsUHH3yA+Pj4VvebP38+FAqF7sFTf0SGr0GjxVNrU3GySAl3B2use2wQAtw4B4qIupdJns5btGgR3nzzTbz11lt4+umn29zX2toa1tbW3VQZEXWUIAiYtzEdezNLYGclw+qHBiDYw0HssojIDJncSNSiRYuQmJiIxMREvPzyy2KXQ0Sd7IMdmdiYmg+pBPjk3jguY0BEojGpEPXGG28gMTERr7zyChYuXCh2OUTUyX5KK8BHOzMBAG9OicbwiGuvyiUi6i4mczrvvffew2uvvYZx48Zh4sSJOHToULPnExISRKqMiDrDiUIFXtp4HADwxLAQ3DsoQOSKiMjcmUyI+uWXXwAAW7duxdatW6953kTWFCUyS+XV9Xj8mxTUNWgxNNwD/ze2l9glERGZTojavXu32CUQURfQaAU8+90x5JfXIsDVDh9Nj4VMKhG7LCIi05oTRUSm57+/n8HezBLYWsrw2cx4ONtZiV0SEREAhigiMmDbTlzEp3+cBwC8c1cMenvzfpZEZDgYoojIIF1U1OHFvyeSP3JzT0zq6yNyRUREzTFEEZHB0WoF/HtDGipqGtDH1wkvjYsQuyQiomswRBGRwVm5Nwv7z5XC1lKGpdP7wcqCf6qIyPDwLxMRGZT0fAXe3XYGALDw9kiE8JYuRGSgGKKIyGDU1Kvx3HfH0KARMC7KC9MG+ItdEhFRqxiiiMhgvPHrSWSVVMPLyQZL7oyGRML1oIjIcDFEEZFB2JpxEd/+mQeJBHh/Wl+uB0VEBo8hiohEd1FRh3mbGpczeHxoCIaEuItcERHR9ZnMbV+Iukq9WovLVSpcUtZBUdOASpUa1So1GjRaaLQCtAJgKZPA2kIKawsZHKwt4OZgBXcHa7jaW8HOSsbTUm3QagXMXd+4nEG0rxxzR4eLXRIRUbswRBFdoV6txakiJdLyKnAstxxpeRW4UFrTode0sZTCzd4a7o7WCHC1Q5CbHYLc7BHkbo8gNzu42luZdcj6dM95HDjfuJzBh9NjuZwBERkNhigyW4IgIL+8FsfyKpCWW4G0vHJkFCpRr9Zes6+lTAJPRxu42lvBwdoC9tYWsLaQQiIBJBIJNFotVA1a1Kk1qKxTo7SqHiVVKqjUWtQ1aFFQUYuCilr8lVdxzWs72lgg2MMBkd6OiPR2QqSPHBFejrC3Nv3umXz6Ev77O5czICLjZPp/pYn+VlnXgOP5imajTCVV9dfs52xniVh/Z/Tzd0FsgDMivZ3gZm8FqfTGRosEQUBNvQalVfUorVbhklKF3LJqZJfUIKe0GhdKqlGoqENlnRp/5VU0C1gSCdDTzR69fZzQx0eOfgHOiPGTw87KdLps5qVKPPttGgQBmDEwgMsZEJHRMZ2/yEQtyChQYM73aThXXAWJBBCE5s9bSCWI9HFCP39nxAY0BqdAN7tOOb0mkUhg//eoVYCbXYv71DVokFtWg7OXKnGqSImThUqcLFLiklKFrJJqZJVU47fjRQAAmVSC3t6O6OfvgrhAZ8QFuCDAtXNq7W55ZTV44Is/UaVSY2BPVyyaFGWU7SAi8yYRhKu/VsybUqmEXC6HQqGAkxPvGG/svjmUg1d/zND97O5ghUHBbujn74x+Ac6I8pHDxlImYoUtK6lS4VSREicKlTieX4HUnApcVNZds5+bvRX6BTijX4AL+gU4o6+fs8GfBsy8VImHVh9BQUUtQjzssf7xwXBzsBa7LCIycmJ8fzNEXYUhyrQ0aLSYvTYV205eAgD4u9rixbERmBDtDdkNnp4TW2FFLVJzy3EstwKpueU4UaBEvab5/C2ZVIIIL0fEB7ogPtAFcQEu8HOxNZhRnu0nL+Hf69OgrFMj2MMe3z2WAE8nG7HLIiITwBBlABiiTI8gCPj9xCUs+uUEihSNozkhHvZ4ZkQYbovxhoXMOK8GU6k1yChQ4tjfwSolp7zF0SoPR2vEBTjrQlUf3+4ffSurrsd/fz+Db//MBQDEB7pg1QP94WLPBTWJqHMwRBkAhijTVa1SY9XebHy+LwvKOjUAoKe7PWYPD8XkWB9YGmmYulLTaFVKTjlScytwokABtbZ5F7eUSRDlI0dfPzmifOXo4yNHWA+HLmn/5UoVvv0zF5/vy4aitgEA8NgtPfF/YyO4lAERdSqGKAPAEGX6Kusa8PXBHKzcm4WKmsYvdn9XW8y+NRR3xvuZRJhqUtegQXqBojFU5ZQjNbe8xSsSrWRSBHvYI6yHI8I9HRDq6QB/Vzv4udhCbmvZ7tOBGq2AnNJq7DtXgj1nL2PP2RLdKccIL0csmhSFQcFundpGIiKAIcogMESZjyqVGmsO5WDlniyUVjcGiwBXO8wdHY7b+/oY3Zyp9hAEAXlltTiWV470fAUyChU4UaBEpUrd6jH2VjL4ONvCxc4KTrYWcLKxhL21BTSCgHq1Fg0aLZS1Dcgpq0F+We0187T6BTjjoSFBmBhtvKdOicjwMUQZAIYo81NTr8a6w7n49I/zulGaPr5O+OCeWIT1cBS5uq6n1TYuOppZXInM4iqcvVSJ85erUVBei5Iq1Q2/npWFFPEBLrgl3B23hnsi0of9iIi6HkOUAWCIMl819Wqs3n8Bn/5xHpV1alhbSLH64QFmfTPc2noNCipqUaSohbJWDWVdA5S1DahWqSGTSmFpIYGVTNq4FparHQLd7OAttzXJUTwiMmwMUQaAIYouKeswd30a9p8rhZ+LLfa+ONxglgggIqKWifH9zQkKRFfp4WSDVyZGAgDyy2uhauFeekRERAxRRFcpVtbh6XWpAICbQt0MckVzIiISH0MU0RVKqlSYsfIQzl+uho/cBu/cGSN2SUREZKAYooj+pqhtwAOf/6kLUN/9azD8XFq+cTARERFDFBEaF4l8ck0KThYp4e5ghTWPDkKAGwMUERG1jiGKCMBXBy7gwPlS2FnJ8PWsQQj2cBC7JCIiMnAMUWT2tFoBX+zPBgDMn9Cbi0MSEVG7MESR2csqqUZ+eS1sLKW4K85P7HKIiMhIMESR2Wtab9bGUgZbKy5nQERE7cMQRWav6aa4Gg0X7yciovZjiCKzZ/H3fd4atFyZnIiI2o8hisyehawxRGm0HIkiIqL2M9kQtWrVKkgkEjg48FJ1apvt37d1adAIOH+5SuRqiIjIWJhkiCooKMALL7wAHx8fsUshI+BsZ4Vbe3kAAOau/wsNGp7WIyKi6zPJEPXEE09g6NChGD16tNilkJF4645oONlY4K+8Crz7+xmxyyEiIiNgciFqzZo1+OOPP7Bs2TKxSyEj4utsi//c1Xiz4c/2ZOGHlHyRKyIiIkNnUiGquLgYc+bMwZIlS+Dn175FE1UqFZRKZbMHmadxfbzxxLAQAMC8jcfxx9nLIldERESGzKRC1FNPPYVevXrhySefbPcxixcvhlwu1z38/f27sEIydC+O7YUpsT5Q/31D4pScMrFLIiIiA2UyIWrjxo345ZdfsHLlSkgkknYfN3/+fCgUCt0jLy+vC6skQyeVSvCfu/ri5lB31NRr8MDnf+LoBQYpIiK6lkmEqKqqKsyePRvPPPMMfHx8UFFRgYqKCtTX1wMAKioqUF1d3eKx1tbWcHJyavYg82ZlIcXKB/pjSIgbqus1ePCLP3GEQYqIiK4iEZpuHGbELly4gJ49e7a5z+TJk/Hjjz9e97WUSiXkcjkUCgUDlZmrrdfg0a+PYP+5UthYSrHsvjiMiOghdllERNQCMb6/TSJE1dXV4dChQ9dsX7JkCf744w9s2bIF7u7u6NOnz3VfiyGKrlTXoMFTa1ORfLoYMqkE79wZg7vi23fRAhERdR+GqE720EMP4YcffkBVVftXoWaIoqs1aLSYtzEdG1Mblz2YPz4Cj/99FR8RERkGMb6/TWJOFFFXspRJ8e7dMXh8aDAAYPGW00j8+QTUXNmciMismfRIlD44EkVtWbHnPN5OOg0AGBrugf/d2w9ONpYiV0VERByJIjJw/xoaguX3xcHGUoo9Zy/jzmUHkFtaI3ZZREQkAoYoohs0PtobGx4fgh5O1sgsrsKUZftxOKtU7LKIiKibMUQR6SHaT46fZt+MaF85yqrrce+qw1i5Jws8O05EZD4Yooj05CW3wfrHB2NyrA80WgFvJZ3Ck2tSUVnXIHZpRETUDRiiiDrA1kqGD6fF4o3JUbCUSbD1xEVM+t9+nLlYKXZpRETUxRiiiDpIIpFg5uAgrH98MHzkNsguqcaUT/Zj09/rShERkWliiCLqJP0CXPDrs7fgljB31DZoMHf9X3j++zQoeXqPiMgkdThEffXVV51RB5FJcLW3wpcPD8Tzo8IhlQCbjxVgwtK9OMobGBMRmZwOh6glS5YgNjYW27Zt64x6iIyeTCrBc6PCsOGJwfB3tUV+eS3u+ewg3t9+lqucExGZkA6HqIyMDDz11FOYNWsWRo8ejbS0tE4oi8j4xQe6IunZWzA1zhdaAfhoZybu/uwgF+ckIjIRHQ5RMpkM//rXv3Du3DmMHDkSo0aNwsyZM5GXl9cZ9REZNUcbS7x/Tyw+mtEPjjYWOJZbgfFL9+CbQznQarmmFBGRMeu0ieU2NjaYN28eMjMz4e3tjZiYGLz44ouoqKjorLcgMlqT+vpgy3O3YGBPV1TXa/DqjxmYvvIQskuqxS6NiIj01GU3IM7IyMADDzyA3NxclJSUdMVbdAnegJi6klYr4OuDF/Cf38+gpl4Dawsp/j0mHLNu6gkLGS+WJSLSlxjf3x0OUUVFRTh16tQ1j9LSUgQHByMqKgo//PBDZ9Xb5RiiqDvkldVg/qZ07DvX+A+MaF853pzSB339ncUtjIjISBlliLK0tERYWBh69+6NqKgoREVFITIyEr169YKVlVVn1dltGKKouwiCgPVH8/Dmb6dQWaeGRALMGBiAF8f2grOd8fUdIiIxGWWIqq+vN8qw1BqGKOpuxZV1WJx0GpuPFQBoXGtq3rgI3BXvB6lUInJ1RETGweBDVJ8+fRAfH4/4+HjExcWhX79+sLe378r6uh1DFInlUFYpXv0xA5nFVQCA+EAXvD45ClE+cpErIyIyfAYfoqRSKSSSxn8ZC4IAqVSK8PBwxMXF6cJVv3794Ojo2GUFdzWGKBJTg0aL1fuz8eGOTNTUayCVAA8MDsLcMeFwsrEUuzwiIoNl8CEqLCwMly5dwvTp0zF48GAcP34cqampSEtLQ2VlJSQSCSQSCUJCQhAfH49169Z1Ze1dgiGKDEGRohZv/nYKvx0vAgB4OFrj5QkRmBLrq/uHDBER/cPgQ1RDQwPee+89vP322/D398eHH36I0aNHAwAyMzORkpKC1NRUpKSkIC0tDaWlpV1WeFdhiCJDsjfzMhb+dAJZf68nFRfgjIW3R/EqPiKiqxh8iGpSVFSEl156CWvXrsWkSZPw/vvvo2fPnl1RX7djiCJDo1JrsGpvNj7ZdQ419RoAwF3xfnhxbC94OtmIXB0RkWEQ4/tbr9X9vL298fXXX2P//v0oKChAZGQkFixYgOpqrr5M1NmsLWSYPTwUu164FVPjfAEAP6TkY/i7u7F893mo1BqRKyQiMk+dsmL56tWr8fLLL0Mmk2HZsmWYNGlSZ9QmCo5EkaE7lluOxF9O4q+8CgBAoJsdFkzojdGRPThfiojMltGMRDWpqKjA3r17UVtbi+HDh6OwsBBffvllJ5VGRC3pF+CCzU8Owfv39IWnozVySmvwr29S8MAXf+LspUqxyyMiMhs3NBK1Zs0apKen6x6FhYUQBAGurq7o27cvYmNjMX78eIwaNaora+5SHIkiY1KtUmPZ7nNYuScb9RotZFIJpsT64pkRoQhyN6013IiI2mLwE8ulUinkcjlGjhypC02xsbHw9/fvyhq7FUMUGaPc0hq8lXQSv5+4BACQSoAp/XzxzIgw9GSYIiIzIMb3t8WNHqBQKLBv3z7U1NRApVJBo9FAEAQEBAR0RX1E1A4Bbnb4bGZ/pOVV4KOdmUg+XYxNqQX48VgBJsf64ukRoQjxcBC7TCIik3JDI1FbtmxBamqqbi2o3NzcxheRSODi4qJbuTwuLg5xcXEICQnpssK7CkeiyBQcz6/A0h2Z2Hm6GEDjyNSkvj54ekQYQj0ZpojI9Bj86byrlZWV6QJVU7jKysqCIAiQSCTQaIzv0muGKDIl6fkKLN2ZiR2nGk/zSSTA7TE+eHZkKEI9jff2TEREVzO6ENUSpVKJlJQUHDt2DHPnzu3Ml+4WDFFkijIKGsPU9pP/hKnbYnzw7IhQhPVgmCIi42cSIcrYMUSRKTtRqMBHOzN1E9AlEmBitDeeHRmGcIYpIjJiDFEGgCGKzMHJQiU+2pmJrScuAmgMUxP6NIapXl4MU0RkfBiiDABDFJmTU0WNYWpLxkXdtgnRXnh2ZBgivPj7T0TGgyHKADBEkTk6fVGJj3eew2/pRbpt46Iaw1SkD/sBERk+higDwBBF5uzMxUp8lJyJpPQiNP1lGBvVA8+ODEOUj1zc4oiI2sAQZQAYooiAs5cq8XHyOfx6vFAXpsZENoapPr4MU0RkeIzuBsSGZt++fZgwYQJcXFxga2uLsLAwvPHGG2KXRWR0wns44uMZ/bBtzlBM6usDiQTYdvISbvt4Hx796igyChRil0hEJDqTCVHr1q3DsGHDIJfL8fXXXyMpKQkvvfQSONBGpL+wHo74aEY/bH9+KCbH+kAqAXacagpTR5CezzBFRObLJE7nFRQUoFevXnjggQewbNmyDr0WT+cRte785Sr8L/kcfkorgPbvvxwjIjzx3Mgw9PV3FrU2IjJvPJ2np1WrVqG6uhovvfSS2KUQmbQQDwd8MC0WO+YOw9R+vpBKgOTTxZj8yX48vPpPpOVViF0iEVG3MYkQtWfPHri6uuL06dOIjY2FhYUFPD098cQTT0CpVIpdHpHJCfZwwPvTYrHz37fizjg/yKQS7DpzGVM+2Y+HVv+JY7nlYpdIRNTlTOJ0XkREBHJycmBpaYn58+dj8ODBOHLkCBYuXIi4uDjs3bsXEomkxWNVKhVUKpXuZ6VSCX9/f57OI7oBF0qq8b9d57D5WAE0f5/nGxrugedGhiE+0EXk6ojIHHCJAz2Fh4cjMzMTixcvxrx583Tbly5dijlz5mD79u0YNWpUi8cmJiZi0aJF12xniCK6cTml1fhf8jlsuiJM3RLmjjmjwhAf6CpydURkyjgnSk9ubm4AgLFjxzbbPn78eABAampqq8fOnz8fCoVC98jLy+u6QolMXKCbPf57d1/s+vetmNbfHxZSCfZmluDO5Qdx/6rDOHqhTOwSiYg6jUmEqJiYmBa3Nw2ySaWtN9Pa2hpOTk7NHkTUMQFudnjnrhjseuFWTB/QGKb2nSvBXZ82hqmUHIYpIjJ+JhGi7rzzTgDAli1bmm1PSkoCACQkJHR7TUQE+LvaYcmdjWFqxsAAXZi6c/lBzPz8MFI5AZ2IjJhJzIkCgEmTJmHbtm145ZVXkJCQgKNHj2LRokUYNWoUfvnll3a/DteJIuo6eWU1+GTXOfyQkg/133Ombu3lgedHhXOdKSLqEE4s74Da2losWrQI69atQ1FREXx8fHDfffdh4cKFsLa2bvfrMEQRdb3c0hr8b1cmNqb+MwF9RIQnnh8Vjmg/3puPiG4cQ5QBYIgi6j45pdX4OPkcNqXm61ZAH9XbE3NGhfNGx0R0QxiiDABDFFH3yy6pxsc7M/HjFbeTGRPZA3NGhSPSh/2QiK6PIcoAMEQRief85Sp8vDMTP/1ViKa/TOOivDBndBgivNgfiah1DFEGgCGKSHzniiuxdOc5/Hr8nzA1IdoLz4wIQ29v9ksiuhZDlAFgiCIyHGcvVWLpzkz8drxIt21sVA88MyKMc6aIqBmGKAPAEEVkeM5crMTHyZn4Lb1INzI1qrcnnhkRxqURiAgAQ5RBYIgiMlzniivxv+Rz+PmvQt0E9GHhHniWNzomMnsMUQaAIYrI8GVdrsInu87jx7R/1pm6OdQdz4wIxaBgN5GrIyIxMEQZAIYoIuORW1qDZbubr4A+qKcrnhsZhsEhbpBIJCJXSETdhSHKADBEERmf/PIafPrHeaw/ko96jRYA0D/QBc+MDMPQMHeGKSIzwBBlABiiiIxXkaIWn/2RhXV/5qJe3Rim+vo747mRoRjey5NhisiEMUQZAIYoIuNXrKzDZ3uysPZwDuoaGsNUH18nPDMiDKN794BUyjBFZGoYogwAQxSR6SipUmHl3ix8czAHNfUaAECElyOeHRmGcVFeDFNEJoQhygAwRBGZnrLqenyxLxtfHriAKpUaABDm6YDZw0NxW4w3LGRSkSskoo5iiDIADFFEpktR04Av9mfji/3ZqKxrDFOBbnZ4YlgIpsb5wtpCJnKFRKQvhigDwBBFZPqUdQ345mAOPt+XjbLqegCAt9wG/xoajOkDAmBrxTBFZGwYogwAQxSR+aipV2Pd4Vys3JuFS0oVAMDdwQqP3ByM+xMC4GhjKXKFRNReDFEGgCGKyPyo1Br8kJKP5bvPI7+8FgDgZGOBh27qiYeHBMHF3krkConoehiiDABDFJH5atBo8XNaIZbtPofzl6sBAHZWMsxMCMQjt/SEp6ONyBUSUWsYogwAQxQRabQCfj9xER8nn8OpIiUAwNpCimkD/PH4sBD4OtuKXCERXY0hygAwRBFRE0EQsOtMMf6XfA6puRUAAAupBFPjfPHkraHo6W4vboFEpMMQZQAYoojoaoIg4GBWKf6XfA4HzpcCAKQSYGKMD54bGYpQT0eRKyQihigDwBBFRG1JySnHJ7vOIfl0MYDGMDWlny/mjAxHgJudyNURmS+GKAPAEEVE7XGiUIGPdmbi9xOXADSe5ps+0B/PjAhDDydOQCfqbgxRBoAhiohuxF95FXh32xnszSwB0DgB/cEhQXhiWAhcuTQCUbdhiDIADFFEpI9DWaV49/czOJpTDgBwsLbArJt74rFbenLRTqJuwBBlABiiiEhfgiBg99nLeG/bGWQUNC6N4GpvhWdHhOLeQYGwsuCNjom6CkOUAWCIIqKOEgQBWzMu4t1tZ3SLdga62eHFsRGYEO0FiUQicoVEpochygAwRBFRZ1FrtPj+aB4+2J6JkqrGe/PF+jtjwcTeGBDkKnJ1RKaFIcoAMEQRUWerVqmxcm8WVuzJQk29BgAwOrIHXhoXgVBPB5GrIzINDFEGgCGKiLpKsbIOH+7MxPdH8qDRCpBJJZg+wB/PjQrjffmIOoghygAwRBFRVztXXIklW85gx6nGNabsrGT419BgPHZLMOytLUSujsg4MUQZAIYoIuouh7NK8faW0/grrwIA4OFojbmjw3FPf3/IpJx8TnQjGKIMAEMUEXUnQRCQlH4R//n9NHJKawAAEV6OeGViJG4Ocxe5OiLjwRBlABiiiEgM9WotvjmUg492ZkJR2wAAGBnhifkTenPyOVE7MEQZAIYoIhJTRU09PtyRiTWHcqDWCrCQSnB/QiCeGxkGF95GhqhVDFEGgCGKiAzB+ctVWJx0CjtOFQMAnGws8OzIMDwwOIgrnxO1gCHKADBEEZEh2X+uBG/8ehKnL1YCAMJ7OODtO6LRn4t1EjUjxve3yfxz5tixY5gyZQp8fHxgZ2eHiIgIvP7666ipqRG7NCIivd0U6o7fnr0FS6ZGw83eCmcvVeGuTw/i5c3pUNQ0iF0ekVkziRB18uRJDBkyBBcuXMCHH36IX3/9FdOnT8frr7+OGTNmiF0eEVGHyKQSTB8YgB1zh2Faf38AwLrDuRj9wR/YfaZY5OqIzJdJrOq2bt061NXVYePGjQgJCQEAjBgxAkVFRVixYgXKy8vh4uIicpVERB3jYm+Fd+6KwdQ4X8zfnI6sy9V4aPUR3J8QgJcn9IadlUn8SScyGiYxEmVpaQkAkMvlzbY7OztDKpXCyopXtBCR6RgU7IakZ2/BQ0OCAABrDuVi4kf7cKJQIW5hRGbGJELUgw8+CGdnZzz55JPIyspCZWUlfv31V3z22WeYPXs27O3tWz1WpVJBqVQ2exARGTobSxkSJ0VhzSOD4C23QXZJNaYuO4D1R/LELo3IbJhEiAoKCsLBgweRkZGBkJAQODk54fbbb8eDDz6IpUuXtnns4sWLIZfLdQ9/f/9uqpqIqONuDnNH0rO3YHgvD6jUWry48Tj+b8NfqGvQiF0akckziSUOLly4gNGjR6NHjx54/vnn4eHhgcOHD+PNN9/EXXfdhc8//7zVY1UqFVQqle5npVIJf39/LnFAREZFqxWw/I/zeG/bGWgFIC7AGZ/N7A8PR2uxSyPqFlwnSk/Tp0/Hrl27kJWV1ezU3erVqzFr1izs3r0bw4YNa9drcZ0oIjJm+zJL8NTaFCjr1PB1tsXnD/VHhBf/lpHp4zpRekpLS0NkZOQ1c58GDBgAAMjIyBCjLCKibndzmDs2z74JPd3tUVBRi7uWH8T+cyVil0VkkkwiRPn4+ODEiROoqqpqtv3gwYMAAD8/PzHKIiISRYiHAzY/NQSDerqiSqXGQ6v/xM9/FYpdFpHJMYkQNWfOHJSUlGD06NFYv349kpOT8fbbb2Pu3LmIjIzE+PHjxS6RiKhbOdtZ4etHBmJitDcaNAKe/fYY1h7OEbssIpNiEnOiAGDXrl1YsmQJjh8/DoVCAX9/f9x+++2YP38+3Nzc2v06nBNFRKZEoxXw+i8n8NXBxgD15pQ+uD8hUOSqiDofJ5YbAIYoIjI1giDg7aRTWLk3GwDw1h19cN8gBikyLZxYTkREnU4ikeDlCb3xr6HBAIBXfszAr8c5R4qooxiiiIjMgEQiwfzxEZiZEAhBAJ7/Pg37MnnVHlFHMEQREZkJiUSCxElRmBjTONn88W+OIvNSpdhlERkthigiIjMik0rw/j19kRDsiup6Df71TQoUtQ1il0VklBiiiIjMjLWFDJ/cGwdfZ1tkl1Tj+e/TwGuMiG4cQxQRkRlyc7DGZzPjYW0hRfLpYmxMLRC7JCKjwxBFRGSm+vjKMXd0OADg7aRTqKipF7kiIuPCEEVEZMZm3dwTYZ4OKKuux4c7MsUuh8ioMEQREZkxS5kUr90eCQD4/kgeFDWcZE7UXgxRRERm7uZQd0R4OaK2QYPvjuSKXQ6R0WCIIiIycxKJBA8NCQIAbOIEc6J2Y4giIiKM6+MFqQQ4c6kSBRW1YpdDZBQYooiICM52VogLcAEA7Mu8LHI1RMaBIYqIiAAA8YGNIepEoVLkSoiMA0MUEREBAII97AEAF0prRK6EyDgwRBEREQAg0K0xROWWVotcCZFxYIgiIiIAgLfcBgCQU8aRKKL2YIgiIiIAgLOtFQBAEIC6Bo3I1RAZPoYoIiICANRrtLr/LZGIWAiRkWCIIiIiAEB+eeNpPG+5DawtZCJXQ2T4GKKIiAgAkFfeuMimv4udyJUQGQeGKCIiAvDPSJSfi63IlRAZB4YoIiICAOSVNY5E+blyJIqoPRiiiIgIAEeiiG4UQxQREQEA8jkniuiGMEQRERG0WgEFf4cojkQRtQ9DFBERobhShXqNFjKpRLdyORG1jSGKiIiarRFlIeNXA1F7sKcQERHy/g5RnA9F1H4MUUREhPwyzociulEMUURE9M9IFNeIImo3higiItItb8CRKKL2Y4giIiKORBHpgSGKiMjMqTVaFFXUAeBIFNGNYIgiIjJzF5V1UGsFWMok6OHINaKI2oshiojIzDXdeNjX2RZSqUTkaoiMh8GHqMrKSrz44osYM2YMPDw8IJFIkJiY2OK+qampGDVqFBwcHODs7IypU6ciKyurewsmIjIy+ZwPRaQXgw9RpaWlWLFiBVQqFaZMmdLqfqdPn8att96K+vp6rF+/Hl988QXOnj2LW265BZcvX+6+gomIjEwer8wj0ouF2AVcT2BgIMrLyyGRSFBSUoJVq1a1uN9rr70Ga2tr/Prrr3BycgIAxMfHIywsDO+++y7eeeed7iybiMhoNI1E+XG1cqIbYvAjURKJBBJJ2+fo1Wo1fv31V9x55526AAU0BrDhw4dj8+bNXV0mEZHR4mrlRPox+BDVHufPn0dtbS1iYmKueS4mJgbnzp1DXV1di8eqVCoolcpmDyIic8I5UUT6MYkQVVpaCgBwdXW95jlXV1cIgoDy8vIWj128eDHkcrnu4e/v36W1EhEZknq1FkXKxn9k8ubDRDfGJEJUk7ZO+7X23Pz586FQKHSPvLy8riqPiMjgFClqIQiAjaUU7g5WYpdDZFQMfmJ5e7i5uQH4Z0TqSmVlZZBIJHB2dm7xWGtra1hbW3dleUREBitPNx/K7rrzT4moOZMYiQoJCYGtrS3S09OveS49PR2hoaGwseEqvEREV/vnyjxOKie6USYRoiwsLHD77bdj06ZNqKys1G3Pzc3Frl27MHXqVBGrIyIyXLobD3M+FNENM4rTeVu2bEF1dbUuIJ08eRI//PADAGDChAmws7PDokWLMGDAANx2222YN28e6urq8Nprr8Hd3R3//ve/xSyfiMhg5XOhTSK9GUWIevLJJ5GTk6P7ecOGDdiwYQMAIDs7G0FBQYiIiMDu3bvx0ksv4a677oKFhQVGjBiBd999Fx4eHmKVTkRk0PLKuLwBkb6MIkRduHChXfvFx8djx44dXVsMEZEJ4UgUkf5MYk4UERHduLoGDYorVQA4J4pIHwxRRERmqqCicRTK3koGZztLkashMj4MUUREZuqPM5cBNM6H4hpRRDeOIYqIyAxptAJW7MkCAEyN8xW5GiLjxBBFRGSGDp4vxUVlHeS2lnhwSJDY5RAZJYYoIiIztCk1HwBwe19vWFvIRK6GyDgxRBERmZkqlRpbMi4CAKbG+YlcDZHxYogiIjIzWzMuorZBg57u9ujn7yx2OURGiyGKiMjMNJ3KuzPOl1flEXUAQxQRkRkpqKjFwaxSAMCUfrwqj6gjGKKIiMzIj8cKIAhAQrAr/LhKOVGHMEQREZkJQRCw8e9TeYeyyjDz88MoUtSKXBWR8WKIIiIyE3/lK5B1uVr3897MEpwvrm7jCCJqC0MUEZGZKKxoPuokkQBDQtxEqobI+DFEERGZieG9PNGrh6Pu5z3/NxxSKa/OI9KXhdgFEBFR97C1kuH354eipl6NE4VK+LtyYjlRR3AkiojIzNhZWWBAkKvYZRAZPYYoIiIiIj0wRBERERHpgSGKiIiISA8MUURERER6YIgiIiIi0gNDFBEREZEeGKKIiIiI9MAQRURERKQHhigiIiIiPTBEEREREemBIYqIiIhIDwxRRERERHpgiCIiIiLSA0MUERERkR4YooiIiIj0wBBFREREpAeGKCIiIiI9MEQRERER6YEhioiIiEgPBh+iKisr8eKLL2LMmDHw8PCARCJBYmJis300Gg3ef/99jBs3Dn5+frCzs0Pv3r0xb948VFRUiFI3ERERmTaDD1GlpaVYsWIFVCoVpkyZ0uI+tbW1SExMRGBgID788EMkJSXhsccew4oVK3DTTTehtra2e4smIiIik2chdgHXExgYiPLyckgkEpSUlGDVqlXX7GNra4vs7Gy4ubnptt16660ICAjA3XffjY0bN+L+++/vzrKJiIjIxBl8iJJIJNfdRyaTNQtQTQYOHAgAyMvL6/S6iIiIyLwZfIjqiOTkZABAVFRUq/uoVCqoVCrdz0qlssvrIiIiIuNn8HOi9FVQUIB58+ahf//+uO2221rdb/HixZDL5bqHv79/N1ZJRERExsokQ1RZWRkmTJgAQRDw/fffQyptvZnz58+HQqHQPXjqj4iIiNrD5E7nlZeXY/To0SgoKEBycjKCg4Pb3N/a2hrW1tbdVB0RERGZCpMKUeXl5Rg1ahSys7Oxc+dOxMTEiF0SERERmSiTCVFNASorKwvbt29Hv379xC6JiIiITJhRhKgtW7aguroalZWVAICTJ0/ihx9+AABMmDABEokEY8eOxbFjx/Dhhx9CrVbj0KFDuuM9PDwQEhIiSu1ERERkmiSCIAhiF3E9QUFByMnJafG57OxsAEDPnj1bPf7BBx/El19+2a73UiqVkMvlUCgUcHJyuuFaiYiIqPuJ8f1tFCNRFy5cuO4+RpAFiYiIyISY5BIHRERERF2NIYqIiIhIDwxRRERERHpgiCIiIiLSg1FMLO9OTRPUeSNiIiIi49H0vd2dF5oxRF2ltLQUAHgjYiIiIiNUWloKuVzeLe/FEHUVV1dXAEBubm63fQiGQKlUwt/fH3l5eWa1PhbbzXabA7ab7TYHCoUCAQEBuu/x7sAQdRWptHGamFwuN6tfviZOTk5stxlhu80L221ezLXdTd/j3fJe3fZORERERCaEIYqIiIhIDwxRV7G2tsbChQthbW0tdindiu1mu80B2812mwO2u/vabRQ3ICYiIiIyNByJIiIiItIDQxQRERGRHhiiiIiIiPRgciGqqqoKc+bMgY+PD2xsbBAbG4vvvvuuXccWFxfjoYcegru7O+zs7DB48GDs3LmzxX137NiBwYMHw87ODu7u7njooYdQXFzcmU25Ifq2e9OmTZgxYwZCQ0Nha2uLoKAg3HfffcjMzLxm31tvvRUSieSax7hx47qiSe2ib7u//PLLFtsikUhw8eLFa/Y3lc+7tc+wpbYb4uddWVmJF198EWPGjIGHhwckEgkSExPbfbyx9vGOtNuY+3hH2m3Mfbwj7TbWPp6cnIxZs2YhIiIC9vb28PX1xeTJk5GSktKu48Xq2ya32ObUqVNx5MgRLFmyBOHh4Vi3bh1mzJgBrVaLe++9t9XjVCoVRo4ciYqKCixduhSenp745JNPMG7cOOzYsQPDhg3T7fvHH39g/PjxmDhxIn766ScUFxfjpZdewsiRI3H06FFRrojQt93vvPMOvLy8sGDBAgQHByMvLw9vv/024uLicOjQIURFRTXbPzg4GGvXrm22zdnZuSua1C76trvJ6tWrERER0Wybm5tbs59N6fNetmzZNfeFrKmpwbhx4xAfHw8vL69mzxna511aWooVK1agb9++mDJlClatWtXuY425j3ek3cbcxzvS7ibG2Mc70m5j7ePLly9HaWkpnnvuOURGRuLy5ct47733kJCQgN9//x0jRoxo9VhR+7ZgQn777TcBgLBu3bpm20ePHi34+PgIarW61WM/+eQTAYBw4MAB3baGhgYhMjJSGDhwYLN9BwwYIERGRgoNDQ26bfv37xcACMuWLeuk1rRfR9p96dKla7YVFBQIlpaWwiOPPNJs+7Bhw4SoqKjOKboTdKTdq1evFgAIR44cue77mNLn3ZIvv/xSACCsWrWq2XZD+7wFQRC0Wq2g1WoFQRCEy5cvCwCEhQsXtutYY+7jHWm3MffxjrTbmPt4R9rdEmPo4y39nlZWVgo9evQQRo4c2eaxYvZtkzqdt3nzZjg4OODuu+9utv3hhx9GYWEhDh8+3OaxvXr1wuDBg3XbLCwscP/99+PPP/9EQUEBAKCgoABHjhzBzJkzYWHxz0DekCFDEB4ejs2bN3dyq66vI+329PS8ZpuPjw/8/PyQl5fX6bV2po60u71M7fNuyeeffw4HBwdMmzatM8vsEk2nG/RhzH28I+025j7ekXa3l6l93i0xhj7e0u+pg4MDIiMjr/t7KmbfNqkQlZGRgd69ezf7PwcAYmJidM+3dWzTfi0de+LEiWav0dq+bb1HV+lIu1uSlZWFnJyca4b5AeD8+fNwdXWFhYUFQkJCsGDBAtTW1upffAd0Rrtvu+02yGQyuLq6YurUqdccY+qfd2ZmJvbu3Yvp06fDwcHhmucN6fPuKGPu453NWPp4ZzDGPt6ZjLmPKxQKpKamtvh7eiUx+7ZJzYkqLS1FcHDwNdub7uhcWlra5rEt3fn56mOb/tvavm29R1fpSLuvplar8cgjj8DBwQHPP/98s+duvvlmTJs2DREREaitrcWWLVvwn//8B/v27cOuXbu69aaPQMfa3TRHJCEhAU5OTkhPT8eSJUuQkJCA/fv3o2/fvs1ew1Q/788//xwA8Mgjj1zznKF93h1lzH28MxlTH+8IY+7jncmY+/js2bNRXV2NBQsWtLmfmH3bpEIUgDaHQK83PHojx7a2b1cPPbemI+1uIggCHnnkEezduxcbN26Ev79/s+fffPPNZj9PmDABQUFBeOGFF/DTTz/hjjvuuPHCO0jfdo8bN67ZFShDhw7FxIkTER0djddeew0//fRTu17LmD9vtVqNr776ClFRUUhISLjmeUP8vDvKmPt4ZzDGPq4vY+/jncGY+/irr76KtWvX4uOPP0Z8fPx19xerbxvPPyvawc3NrcUkWVZWBqDl9HmjxzZd1dHavm29R1fpSLubCIKARx99FGvWrMGXX36JyZMnt+u977//fgDAoUOHbqDiztEZ7b5SUFAQbr755mZtMdXPGwCSkpJw8eJFPProo+1+bzE/744y5j7eGYyxj3c2Y+njncVY+/iiRYvw5ptv4q233sLTTz993f3F7NsmFaKio6Nx6tQpqNXqZtvT09MBAH369Gnz2Kb92jq26b+t7dvWe3SVjrQb+OeP6+rVq7Fq1SpdJ7oRYgz7drTdLREEoVlbTPHzbvL555/DysoKM2fOvOEaxB7m14cx9/GOMtY+3hWMoY93FmPs44sWLUJiYiISExPx8ssvt+sYUfv2DV/PZ8CSkpIEAMJ3333XbPu4ceOue+n3smXLBADCoUOHdNsaGhqEqKgoYdCgQc32HThwoNCnT59mr3fw4EEBgLB8+fJOak37daTdWq1WeOSRRwSJRCKsWLHiht/7nXfeEQAIP/744w0f21EdaXdLsrKyBAcHB2HKlCnNtpvS592kqKhIsLCwEO65554bem8xP++r3eil38bcx690o+025j5+pc641N9Y+viV9G23Mfbx119/XQAgvPLKKzd0nJh926RClCA0rpXj4uIirFixQkhOThYee+wxAYCwZs0a3T6zZs0SZDKZcOHCBd22uro6ISoqSvD39xfWrl0rbN++XbjjjjsECwsLYffu3c3eY9euXYKFhYVwxx13CNu3bxfWrl0r+Pv7C3369BHq6uq6ra1X0rfdTz/9tABAmDVrlnDw4MFmj9TUVN1+e/bsEcaOHSt8+umnwrZt24Sff/5ZePLJJwWZTCaMGDFC0Gg03dreJvq2e+TIkcKiRYuEzZs3Czt37hQ+/PBDwcfHR3B0dBTS09ObvYcpfd5NlixZIgAQtm3b1uLrG+rnLQiNIXLDhg3CF198IQAQ7r77bmHDhg3Chg0bhOrqakEQTLOP69tuY+/j+rbb2Pu4vu1uYmx9/N133xUACOPGjbvm9/TgwYO6/Qytb5tciKqsrBSeffZZwcvLS7CyshJiYmKEb7/9ttk+Dz74oABAyM7Obrb94sWLwgMPPCC4uroKNjY2QkJCgrB9+/YW32fbtm1CQkKCYGNjI7i6ugoPPPBAi4uFdRd92x0YGCgAaPERGBio2y8zM1OYMGGC4OvrK1hbWws2NjZCdHS08NZbb4n2R0YQ9G/3nDlzhMjISMHR0VGwsLAQfHx8hPvvv184c+ZMi+9jKp93k/DwcCEoKEi3oN/VDPXzFoS2f2eb2mqKfVzfdht7H9e33cbexzvyey4IxtfHhw0b1mp7rzxpZmh9WyIIgnAjp/+IiIiIyMQmlhMRERF1F4YoIiIiIj0wRBERERHpgSGKiIiISA8MUURERER6YIgiIiIi0gNDFBEREZEeGKKIiIiI9MAQRURERKQHhigiMlhJSUmQSCTNHk5OToiPj8f69eu77H1nzZoFS0tL1NfXt7rPhAkTYGdnh9zc3C6rg4gMG0MUERms1NRUAMBPP/2EgwcP4sCBA1i5ciWqq6sxY8YMHD9+vEveNyYmBmq1GmfOnGnx+d9//x1btmzBvHnzEBAQ0CU1EJHhsxC7ACKi1qSmpkIul2PSpEm6bYMHD4Zarcb999+PY8eOISYmptPft+k1MzIyEB0d3ew5tVqNuXPnIigoCC+++GKnvzcRGQ+GKCIyWCkpKYiNjb1me35+PgCgd+/eXfK+V4aoq3366ac4efIkNm7cCBsbmy55fyIyDjydR0QGqbS0FLm5uejbty/UajXUajWKi4vxzTff4K233sKjjz6KgQMHdsl7u7u7w8vL65oQVV5ejsTERIwcORJTp07tkvcmIuMhEQRBELsIIqKrbd++HWPGjLlmu4WFBRITE7FgwYIuff+xY8fi/PnzOHfunG7bc889h2XLliEtLQ1RUVFd+v5EZPg4EkVEBiklJQUAsGnTJhw5cgRHjhzB1q1bMXHiRLz22mvYtGlTi8ft3r37miv6WnukpaW1+v7R0dHIzs5GTU0NAODMmTNYvnw5Zs+ezQBFRAA4J4qIDFRqaipsbW0xadIkyGQy3fZhw4bByckJK1eubPGUWq9evbBy5cp2vUdbV9bFxMRAq9Xi5MmT6N+/P+bOnQtnZ2ckJibecFuIyDQxRBGRQUpNTUV0dHSzAAUAlpaWkMlkqK2tbfE4b29vPProox1+/ysnl5eWliIpKQkrV66Es7Nzh1+biEwDT+cRkcFRKBTIyspq8cq8n376CXV1dRg6dGiX1tC7d29YWFggLS0Nc+fORXx8PGbNmtWl70lExoUjUURkcFJTUyEIAuzt7XHo0CEAjVfGHThwAB988AFiYmLwwgsvdGkN1tbWCA8Px4oVK1BXV4f9+/dDKuW/O4noHwxRRGRwmlYq/+CDD/DBBx8AAOzt7REeHo6XX34Zc+bMgZ2dXZfXERMTg5MnT2LmzJkYPHhwl78fERkXLnFAREREpAeOTRMRERHpgSGKiIiISA8MUURERER6YIgiIiIi0gNDFBEREZEeGKKIiIiI9MAQRURERKQHhigiIiIiPTBEEREREemBIYqIiIhIDwxRRERERHr4fyZNTnBkcSr+AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Evaluate absolute magnitudes\n",
    "\n",
    "M_B = -2.5*np.log10(F_obs_B)\n",
    "M_V = -2.5*np.log10(F_obs_V)\n",
    "\n",
    "# Plot the CMD\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "plt.plot(M_B-M_V, M_V)\n",
    "\n",
    "plt.xlim(0, 2)\n",
    "plt.ylim(12.5,-2.5)\n",
    "\n",
    "plt.xlabel('$B-V$')\n",
    "plt.ylabel('$M_V$');"
   ]
  },
  {
   "cell_type": "raw",
   "id": "fb835442",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/restructuredtext",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    ".. rubric:: Footnotes\n",
    "\n",
    ".. [#isochrone] Generated using the `MIST <mist_>`__ web interpolator for an age :math:`\\log(t/\\yr) = 9.6` (non-rotating, solar-metallicity)."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}