{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "76e4c6f0-7178-434f-943a-0c4737beb903",
   "metadata": {},
   "source": [
    "# Garnet Geotherms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c10667c8-945a-4958-9020-8e59687e6bd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import Thermobar as pt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "377d03fb-0a18-4150-a7fa-ad4c3de85f1c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/sinan/.local/lib/python3.10/site-packages/pandas/core/indexing.py:2115: FutureWarning: In a future version, the Index constructor will not infer numeric dtypes when passed object-dtype sequences (matching Series behavior)\n",
      "  new_ix = Index(new_ix)\n",
      "/home/sinan/.local/lib/python3.10/site-packages/pandas/core/indexing.py:2115: FutureWarning: In a future version, the Index constructor will not infer numeric dtypes when passed object-dtype sequences (matching Series behavior)\n",
      "  new_ix = Index(new_ix)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>P_kbar_calc</th>\n",
       "      <th>T_K_calc</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>62.186165</td>\n",
       "      <td>1665.508418</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>59.185112</td>\n",
       "      <td>1638.988377</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>60.174362</td>\n",
       "      <td>1663.121565</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>61.142703</td>\n",
       "      <td>1644.163358</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>62.367269</td>\n",
       "      <td>1680.010277</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>60.658036</td>\n",
       "      <td>1653.928696</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>60.041928</td>\n",
       "      <td>1648.787123</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>62.232799</td>\n",
       "      <td>1654.135943</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>65.247289</td>\n",
       "      <td>1667.306600</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>57.612930</td>\n",
       "      <td>1653.004617</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>61.387758</td>\n",
       "      <td>1627.492544</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>49.661312</td>\n",
       "      <td>1435.644159</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>57.120652</td>\n",
       "      <td>1646.498971</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    P_kbar_calc     T_K_calc\n",
       "0     62.186165  1665.508418\n",
       "1     59.185112  1638.988377\n",
       "2     60.174362  1663.121565\n",
       "3     61.142703  1644.163358\n",
       "4     62.367269  1680.010277\n",
       "5     60.658036  1653.928696\n",
       "6     60.041928  1648.787123\n",
       "7     62.232799  1654.135943\n",
       "8     65.247289  1667.306600\n",
       "9     57.612930  1653.004617\n",
       "10    61.387758  1627.492544\n",
       "11    49.661312  1435.644159\n",
       "12    57.120652  1646.498971"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Here let's recalculate T-P from the study of Tappe et al. (2021).\n",
    "#First creating the composition arrays.\n",
    "out=pt.import_excel('Tappe2021.xlsx', sheet_name=\"Sheet1\")\n",
    "myOpxs = out['Opxs']\n",
    "myCpxs = out['Cpxs']\n",
    "myGts = out['Gts']\n",
    "#Iterative solution of Brey and Köhler (1990) for the given dataset.\n",
    "calc=pt.calculate_pyroxenes_garnet_press_temp(opx_comps=myOpxs, cpx_comps=myCpxs,\n",
    "                                              gt_comps = myGts, equationP=\"P_Brey1990\",\n",
    "                                              equationT=\"T_Brey1990\",T_K_guess = 1300)\n",
    "calc\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4aa22b59",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUAAAAFCCAYAAACAWAShAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5W0lEQVR4nO3deZyO9f7H8dfH2PcsKUdMWSplVx2kDJWiRL8TSdEMDQmVrWQnSxNGlN2g6HBKixyhU2MrJctkL1Q6SKflVJYw+Pz+uO9xBjO3a2bu+77uue/P8/G4Hu71ut5fw2eu77V8v6KqGGNMJMrjdgBjjHGLFUBjTMSyAmiMiVhWAI0xEcsKoDEmYlkBNMZErLxuB0ivTJkyGh0d7XaMixw7dowiRYq4HcMvrC2hK5zaE2pt2bRp08+qWvbC10OqAEZHR7Nx40a3Y1xk1apVNGnSxO0YfmFtCV3h1J5Qa4uI7M/odesCG2MilhVAY0zEsgJojIlYIXUM0JjcJDU1lQMHDnDixAm/rK9EiRLs2rXLL+tym1ttKViwIBUqVCBfvnyOPm8F0JhsOnDgAMWKFSM6OhoRyfH6jhw5QrFixfyQzH1utEVV+eWXXzhw4ABXX321o+9YF9iYbDpx4gSlS5f2S/EzOScilC5dOkt75FYAjckBK36hJas/DyuAxuRiUVFR1K5dm1q1alG3bl0+/fTTbK1n4sSJHD9+PMP3li5dSp06dahVqxbVq1dn+vTpPtc1bNgwJk2a5PMzKSkpLFu27NzzJUuWMHbs2KwHzyE7BmhMLlaoUCFSUlIAWLFiBQMGDGD16tVZXs/EiRN55JFHKFy48Hmvp6amEh8fz4YNG6hQoQInT57ku+++y3HulJQUNm7cSIsWLQBo1aoVrVq1yvF6syrX7gGePHmSdevWuR3DGMcWbFtA9MRo8gzPQ/TEaBZsW+DX9f/xxx9cdtll556/9NJL3HTTTdSsWZOhQ4cCnlvUWrZsSa1atbjxxhtZtGgRkyZN4tChQ8TExBATE3PeOo8cOcLp06cpXbo0AAUKFODaa68F4LvvvqNp06bUrFmTZs2a8f3331+UqUmTJufu7vr555+Jjo7m1KlTDBkyhEWLFlG7dm0WLVrE3Llz6dGjh8/1PvbYY/Tq1YuGDRtyzTXX8NZbb+X47yzXFsABAwZwxx13hM1lAya8Ldi2gPj349n/+34UZf/v+4l/Pz7HRfDPP/+kdu3aXHfddXTp0oXBgwcDsHLlSvbs2cOGDRtISUlh06ZNrFmzhuXLl1O+fHm+/PJLtm/fzt13302vXr0oX748ycnJJCcnn7f+UqVK0apVKypVqkT79u1ZsGABZ8+eBaBnz5506tSJrVu30qFDB3r16uUoc/78+RkxYgTt2rUjJSWFdu3anfe+r/X+8MMPrFu3jqVLl/Lcc8/l5K8OyMUF8Nlnn6Vo0aJ07NiR06dPux3HGJ8GfjSQ46nnH2M7nnqcgR8NzNF607rAu3fvZvny5XTs2BFVZeXKlaxcuZI6depQt25ddu/ezZ49e6hRowYffvghzz77LGvXrqVEiRKX3MasWbP46KOPuPnmmxk3bhxxcXEArF+/nocffhiARx991G89Ml/rbd26NXny5KF69er8+OOPOd5Wri2A5cqVY9q0aWzcuNGVg6fGZMX3v1/cPfT1enY0aNCAn3/+mZ9++glVZcCAAaSkpJCSksLevXvp3Lkz1apVY/PmzdSoUYNBgwYxYsQIR+uuUaMGzzzzDB9++CGLFy92nClv3rzn9hj9ccF4gQIFzj32x4RuubYAAvztb3+jffv2DB8+/NyBYGNCUcUSFbP0enbs3r2bM2fOULp0aZo3b05SUhJHjx4F4ODBg/znP//h0KFDFC5cmEceeYR+/fqxefNmAIoVK8aRI0cuWufRo0dZtWrVuecpKSlUqlQJgIYNG7Jw4UIAFixYQOPGjS/6fnR0NJs2bQI475hdZttzul5/yfVngV955RVWrVpFp06d2LBhw3m/IYwJFaOajSL+/fjzusGF8xVmVLNROVpv2jFA8OwRzZs3j6ioKO666y527dpFgwYNAChatCjz589n79699OvXjzx58pAvXz6mTp0KQHx8PHffffe5Y4FpVJWEhAS6du1KoUKFKFKkCHPnzgVg8uTJxMbG8tJLL1G2bFnmzJlzUb6+ffvStm1bZsyYQcuWLc+9HhMTw9ixY6lduzYDBgw47ztO1us3qhoyS7169TQ7li5dqoA+//zz2fr+pSQnJwdkvW6wtvjPzp07s/T5+Vvna6XESirDRCslVtL5W+ef9/4ff/zhz3iucrMtGf1cgI2aQc3J9XuAAC1btqRz586MHTuWVq1accstt7gdyZiLdKjRgQ41Orgdw6STq48BpjdhwgQqVKhAx44dM72i3Rhj0gubAli8eHHmzJnD119/zcCBObu0wBgTGcKmAAI0bdqUHj168PLLL2frdiBjskr9cCmG8Z+s/jzCqgACjB07lsqVKxMbG5vpaXZj/KFgwYL88ssvVgRDhHrHAyxYsKDj7wTsJIiIXAssSvfSNcAQVZ0YqG0C507TN27cmH79+jFt2rRAbs5EsAoVKnDgwAF++uknv6zvxIkTWfrPG8rcakvaiNBOBawAqupXQG0AEYkCDgLvBGp76TVq1Ig+ffowbtw42rRpQ/PmzYOxWRNh8uXL53jkYSdWrVpFnTp1/LY+N+WWtgSrC9wM2KeqGc7NGQgjR47k+uuvp3Pnzvz222/B2qwxJhcJVgF8CPh7kLYFeHaFX3vtNQ4fPszTTz8dzE0bY3IJCfQBXBHJDxwCblDVi4ZvEJF4IB6gXLly9dLuAfSXpKQkXn/9dV544QUaNWqUrXUcPXqUokWL+jWXW6wtoSuc2hNqbYmJidmkqvUveiOj20P8uQD3AyudfDa7t8L5cvLkSa1du7aWK1dOf/7552ytw+1brvzJ2hK6wqk9odYWMrkVLhhd4PYEufubXv78+Zk3bx6//vorTz75pFsxjDEhKKAFUESKAHcCbwdyO5dSs2ZNhg0bxqJFi/jHP/7hZhRjTAi5ZAEUkctFpI2IPCkicSJys4g4KpyqekxVS6vq7zmPmjP9+/fn5ptvpnv37hw+fNjtOMaYEJBpIRORGBFZAfwTuAe4EqgODAK2ichwESkenJg5lzdvXubOncvRo0fp1q2bXb1vjPG5B9gCeFxVb1LVeFUdpKp9VbUVUAvYgqd7m2tcf/31jBo1ivfee4/58+e7HccY47JMC6Cq9lPVDCcsUNXTqvquqjqfHCBEPP300zRq1IiePXty8OBBt+MYY1zk5BhgSRHpJSITRGRS2hKMcIEQFRXF3LlzSU1NpUuXLtYVNiaCOTmZsQyIBrYBm9ItuVaVKlV48cUXWb58ObNnz3Y7jjHGJU4GQyioqr0DniTIunfvzttvv03v3r258847z810ZYyJHE72AF8XkcdF5EoRKZW2BDxZgOXJk4ekpCRUlbi4uHNzlxpjIoeTAngKeAlYz/+6vxsDGSpYoqOjGT9+PB9//LGNG2hMBHJSAPsAVVQ1WlWv9i7XBDpYsDz++OPcdddd9OvXj3379rkdxxgTRE4K4F4gbKdZExFmzZpF3rx5iY2Nta6wMRHESQE8BqSIyPRwuAwmI1dddRUvv/wya9euZdKksGqaMcYHJ2eB3/UuYa1Tp04sXryYAQMG0KJFC6pVq+Z2JGNMgPksgN65PB5T1Zgg5XGNiDB9+nRuvPFGHnvsMdauXUtUVJTbsYwxAeSzC6yqZ4CzIlIiSHlcVb58eSZPnsz69euZMGGC23GMMQHmpAt8FM/oLx/iOR4IgKr2ClgqFz388MMsXryYQYMG0bJlS6pXr+52JGNMgDgpgG/j8oCmwSQiTJs2jRtuuIHHHnuMTz/91O1IxpgAuWQBVNV5IlIIqKieuX7D3uWXX86UKVNo27YtCQkJNGzY0O1IxpgAcDIazH1ACrDc+7y2iCwJcC7XPfjgg7Rt25Zhw4bxzTffuB3HGBMATq4DHAbcDPwGoKopQNjcCeLLq6++ymWXXcbYsWNJTU11O44xxs+cFMDUDOb0iIjbJcqUKcO0adPYs2cPY8aMcTuOMcbPnBTAHSLyMBAlIlVFZDIQMWcG2rRpQ7NmzRg5ciQpKSluxzHG+JGTAtgTuAE4iWd+3z+ApwOYKeT06tWLMmXK0KlTJ06dOuV2HGOMn1yyAKrqcVUd6J0cqb738YlghAsVxYsXZ8aMGWzdupUXXnjB7TjGGD9xcha4mojMEJGVIvJx2hKMcKHkvvvuo2PHjowePZpNm3L1jADGGC8nF0K/CUwDZgFnAhsntE2cOJF//etfdOrUiU2bNlGgQAG3IxljcsDJMcDTqjpVVTeo6qa0JeDJQtBll13GzJkz2bFjB8OGDXM7jjEmhzItgOnm/nhfRLqH25wg2dWiRQvi4uJISEhgw4YNbscxxuSArz3AtLk/OgH98Fz6kqU5QbxzCr8lIrtFZJeINMhp4FAwYcIEypcvT6dOnThxIqLOBxkTVjItgOnm/rg+3VwgV6vq1YDTIVJeBpar6nVALWBXziO7r0SJEsyaNYvdu3czdOhQt+MYY7LJyTHAjC56vuSF0N4xBG8DZgOo6ilV/S1L6UJY8+bNefzxxxk3bhzr1693O44xJht8HQO8QkTqAYVEpI6I1PUuTYDCDtZ9NfATMEdEtojILBEp4pfUIWLcuHFUqFCB2NhY/vzzT7fjGGOySFQ14zdEOgGPAfU5/5jfEWCuqvocI1BE6gOfAY1U9XMReRn4Q1UHX/C5eCAeoFy5cvUWLlyYzaYEztGjRylatGiG723atIm+ffvStm1bnnjiiSAnyzpfbcltwqktEF7tCbW2xMTEbFLV+he9oao+F+D/LvWZTL53BfBduueNgX/6+k69evU0FCUnJ/t8v2vXrioi+sknnwQnUA5cqi25STi1RTW82hNqbQE2agY1x8kxwI9EZIKIbPQu453MEaKqh4F/i8i13peaATudVOvc5qWXXqJixYrWFTYml3FSAGfj6fa29S5/AHMcrr8nsEBEtgK1gdHZyBjyihUrRlJSEl9//TWDBg1yO44xxiEnt8JVVtX/S/d8uIikOFm5egZPvbjfHYaaNm3KE088QWJiIg888ACNGjVyO5Ix5hKc7AH+KSK3pj0RkUaA9fMykJCQQKVKlXjsscc4fvy423GMMZfgpAA+AbwqIt+JyHfAK0C3gKbKpYoWLcqcOXPYu3cvAwcOdDuOMeYSnIwHmKKqtYCaQE1VraOqXwY+Wu7UpEkTnnzySV5++WXWrl3rdhxjjA9OxgMsJyKzgUWq+oeIVBeRzkHIlmuNHTuW6Oho4uLiOHbs2KW/YIxxhZMu8FxgBVDe+/xrImxI/KwqWrQoSUlJ7N27l+eff97tOMaYTDgpgGVU9R94Z4JT1dNE+MCoTjRp0oQePXowadIk1qxZ43YcY0wGnBTAYyJSGlAAEfkrcOE0mSYDY8eO5ZprriE2Nta6wsaEICcFsDewBKgsIp8Ar+G5wNlcQpEiRUhKSuKbb76xrrAxIcjJWeDNwO1AQ6ArcIOqbg10sHBx++2307NnT+sKGxOCfA2HVVxEqsK5437VgbpAexEpF6R8YWHMmDFUrlzZusLGhBhfe4DjgPT3c43Bc1vbbcDwQIYKN+m7wgMGDHA7jjHGy1cBvAmYl+75EVXtpapdgBsDGyv83HbbbfTq1YvJkyezevVqt+MYY/BdAPN6x9FK82i6xyUDEye8jR49msqVK9sF0saECF8F8KyIXJH2RFW3A4jIX/BeE2iyxrrCxoQWXwXwJTxzAt8mIsW8y+3Au973TDbcdttt9OzZ07rCxoQAX9NizgcGAy8A3wHfAiOAIar6elDShakxY8ZwzTXXWFfYGJf5vA5QVZer6m2qWlpVy6jq7ar6QbDChSvrChsTGpzcCWICIO0C6cmTJ9sF0sa4xAqgi6wrbIy7rAC6qEiRIsyePZt9+/bZCNLGuCBLBVBElgYqSKRKP2yWjSBtTHBldQ/wLwFJEeHGjBlzbgRpm0zJmODJagHcEpAUES79CNLWFTYmeLJUAFU1LlBBIl36yZTWrVvndhxjIoKdBAkhY8eOpVKlSsTGxlpX2JggsAIYQtJ3hQcPHux2HGPCnpNpMWsEI4jxiImJ4YknniAxMZFPP/3U7TjGhDUne4BTRGSDiHQXkRIBT2R48cUXqVixIrGxsfz5559uxzEmbDmZE6Qx0AG4CtgkIm+IyJ1OVi4i34nINhFJEZGNOcwaMYoVK8bs2bP5+uuvrStsTAA5OgaoqnuAQcCzeCZImiQiu0XkAQdfj1HV2qpaPwc5I06zZs3o2rUriYmJrF+/3u04xoQlJ8cAa4pIIrALaArcp6rXex8nBjhfREtISKBChQrWFTYmQJzsAU4GNgO1VPVJ7zSZqOohPHuFviiwUkQ2iUh8zqJGnuLFizNr1iy++uorhg4d6nYcY8KOnD/txwVvikQBr6vqw9lauchfVPWgiFwOfAj0VNU1F3wmHogHKFeuXL2FCxdmZ1MBdfToUYoWLera9seNG8cHH3zA5MmTqV69eo7W5XZb/Cmc2gLh1Z5Qa0tMTMymDA/DqarPBVgL5L/U5xysZxjQ19dn6tWrp6EoOTnZ1e3//vvvetVVV+n111+vf/75Z47W5XZb/Cmc2qIaXu0JtbYAGzWDmuOkC/wt8ImIDBaR3mnLpb4kIkVEpFjaY+AuYLvDgm3SKV68ODNnzmTXrl2MGDHC7TjGhA0nBXAfsNT72WLexcm+bTlgnYh8CWwA/qmqy7MbNNI1b96cuLg4XnzxRb744gu34xgTFvI6+MxOVX0z/Qsi8uClvqSq3wC1shvMXGz8+PGsWLGC2NhYNm3aRIECBdyOZEyu5mQPMKNZe2wmHxeULFmS6dOns2PHDl544QW34xiT62W6Bygi9wAtgL+IyKR0bxUHTgc6mMlYy5Yt6dixI2PGjKFNmzbUrVvX7UjG5Fq+9gAPARuBE8CmdMsSoHngo5nMJCYmUrZsWWJjYzl16pTbcYzJtXxNjP6lqs4DqgD/AD5T1Xmq+raq/jdoCc1FSpUqxfTp09m6dSujR492O44xuZaTY4B3AynAcgARqS0iSwIZylxaq1atePjhhxk1ahRffvml23GMyZWcFMBhwM3AbwCqmgJcHbBExrFJkyZRqlQpYmNjSU1NdTuOMbmOkwKYqqq/X/Ba5vfPmaApXbo0U6dOZcuWLSQkJLgdx5hcx0kB3CEiDwNRIlJVRCYDNlRxiHjggQdo27Ytw4cPZ/t2u9HGmKxwUgB7AjcAJ4G/A38ATwcwk8miV155hRIlShAbG8vp03aFkjFOORkR+riqDlTVm1S1vvfxiWCEM86ULVuWV199lY0bNzJ+/Hi34xiTazgZELW+iLwtIptFZGvaEoxwxrkHH3yQBx54gKFDh7J792634xiTKzjpAi8A5gL/B9yXbjEhRESYMmUKRYoUITY2ljNnzrgdyZiQ56QA/qSqS1T1W1Xdn7YEPJnJsnLlyjFp0iQ+++wzXn75ZbfjGBPynBTAoSIyS0Tai8gDaUvAk5lsefjhh7nvvvsYOHAge/bscTuOMSHNSQGMBWrjuSMkrft7bwAzmRwQEaZNm0bBggWJi4vj7NmzbkcyJmQ5KYBpZ387qWqsd4kLeDKTbeXLlycxMZF169bx6quvuh3HmJDlpAB+KiI5m4nHBF2nTp245557eO655/jmm2/cjmNMSHJSAP8KpIjIV95LYLbZZTChT0SYPn06UVFRdOnSxbrCxmTA6WgwVfFMapR2/M8ug8kFrrrqKsaPH09ycjIzZsxwO44xIcfJnSD7gZL87wRISbsMJvfo0qULd9xxB/369WP/fvuxGZOekztBnsJzMfTl3mW+iPQMdDDjHyLCzJkzUVXi4+PT5mg2xuCsC9wZuEVVh6jqEDzHBB8PbCzjT9HR0SQkJLBy5Uo++OADt+MYEzKcFEAB0t9Xdcb7mslFunXrxu23386UKVM4ePCg23GMCQlOCuAc4HMRGSYiw4DPgNkBTWX8Lk+ePMyaNYvTp0/TrVs36wobg7OTIBPw3A3yq3eJVdWJAc5lAqBKlSp06dKFpUuXsmDBArfjGOM6J3uAAN8Cq4B1gIiITUabS7Vp04YGDRrQq1cvDh8+7HYcY1zl5CzwSGArMAkY713GBTiXCZCoqCiSkpI4fvw43bt3t66wiWhO9gDbApVVtYmqxniXpk43ICJRIrJFRJZmP6bxp+uuu47hw4fzzjvv8Oabb7odxxjXOCmA2/FcCJ1dTwG7cvB9EwB9+vThpptu4sknn+Snn35yO44xrnBSAMcAW0RkhYgsSVucrFxEKgAtgVk5CWn8L2/evCQlJfH777/Ts6dd124iU14Hn5kHvAhsA7J6R/1EoD9QLIvfM0Fw4403MnjwYIYMGcJDDz1E69at3Y5kTFA5KYDHVXVSVlcsIvcC/1HVTSLSxMfn4oF48AzpvmrVqqxuKuCOHj0akrmy48K2NGjQgCpVqtC5c2fy5MlD8eLF3QuXReH0c4Hwak+uaYuq+lyACXi6wQ2AummLg++NAQ4A3wGHgePAfF/fqVevnoai5ORktyP4TUZt2bx5s+bNm1c7deoU9Dw5EU4/F9Xwak+otQXYqBnUHCfHAOvguf93NFm4DEZVB6hqBVWNBh4CPlbVR7JQm02Q1KlTh+eee4558+axbNkyt+MYEzQ+C6CIRAFL9H+Xv2T5MhiTOwwaNIjq1avTtWtXfv/9d7fjGBMUPgugqp4B2ud0I6q6SlVtIqUQVqBAAebMmcOhQ4fo37+/23GMCQonXeBPROQVEWksInXTloAnM0F3880306dPH2bMmMFHH33kdhxjAs7JWeDa3j9HpHtNAesGh6Hhw4fz7rvv0qVLF7Zt20bRokXdjmRMwDgZDebC4392DDCMFSpUiKSkJPbv38/zzz/vdhxjAsrJYAjlRGS2iHzgfV5dRDoHPppxy6233kqPHj2YPHkya9eudTuOMQHj5BjgXGAFUN77/Gvg6QDlMSFizJgxXH311XTu3Jnjx4+7HceYgHBSAMuo6j/w3ganqqc5f4h8E4aKFCnCzJkz2bNnD0OGDHE7jjEB4aQAHhOR0nhOfCAifwXsQrEI0KxZM+Lj40lMTOTzzz93O44xfuekAPYGlgCVReQT4DXAhg+JEAkJCZQvX564uDhOnjzpdhxj/MrJWeDNwO1AQ6ArcIOqbg10MBMaSpQowfTp09m5cyejRo1yO44xfuV0TpCbgVp4BkJoLyIdAxfJhJoWLVrQsWNHxowZQ0pKittxjPEbJ5fBvI5n8INbgZu8S/0A5zIhJjExkTJlyhAXF0dqaqrbcYzxCyd7gPWBRqraXVV7epdegQ5mQkupUqWYOnUqW7ZsISEhwe04xviF0zlBrgh0EBP6WrduTbt27RgxYgQ7duxwO44xOZZpARSR971zf5QBdmZnThATfiZPnkzx4sWJi4vjzBm7HNTkbr4GQ7C5f81FypYty+TJk2nfvj0TJ06kT58+bkcyJtt8FcCBqnpX0JKYXKNdu3YsXLiQQYMG0apVK6pWrep2JGOyxdcxwDJBS2FyFRFhypQpFCxYkM6dO3P2bFYnCzQmNPgqgCVF5IHMlqAlNCGpfPnyJCYmsnbtWqZMmeJ2HGOyxVcXuARwLyAZvKfA2wFJZHKNTp06sXDhQp577jnuvfdeoqOj3Y5kTJb4KoD7VTUuaElMriMizJgxgxtuuIHHH3+clStXIpLR70tjQpOvLrD9SzaXVLFiRRISEvjXv/5FUlKS23GMyRJfBfCS9/uK/bo3QNeuXbn99tvp3bs3Bw4ccDuOMY75KoCTRaSniFRM/6KI5BeRpiIyD+gU2HgmN8iTJw+zZs0iNTWVJ554AlV1O5IxjvgqgHfjGfn57yJySER2isg3wB48cwVPVNW5QchocoEqVaowatQoli5dyhtvvOF2HGMcybQAquoJVZ2iqo2ASkAzoK6qVlLVx1V1S9BSmlyhV69eNGjQgF69evHjjz+6HceYS3I0HqCqpqrqD6r6W4DzmFwsKiqKpKQkjh07xpNPPul2HGMuyemAqMY4ct111zFs2DAWL17MW2+95XYcY3yyAmj8rm/fvtSrV48nn3ySX375xe04xmTK13BYK3OyYhEpKCIbRORLEdkhIsNzsj6Te+TNm5ekpCR+/fVXnnrqKbfjGJMpX3uAZXO47pNAU1WtBdQG7vZOqWkiQM2aNRk4cCALFixg6dKlbscxJkM+7wX2NeiBqvq8F1g9F4Md9T7N513sArEI8vzzz7N48WK6du3Kjh07KFmypNuRjDmPrz3AtMEQ7stgudfJykUkSkRSgP8AH6qqza4dQfLnz09SUhKHDx+mX79+bscx5iKS2VX7IrJZVev6ZSMiJYF3gJ6quv2C9+KBeIBy5crVW7hwoT826VdHjx6laNGibsfwCzfaMn36dBYuXMhLL71E/fr+m1AwnH4uEF7tCbW2xMTEbFLVi//xqWqGC7Als/eyswBDgL6+PlOvXj0NRcnJyW5H8Bs32nL8+HGtVq2aVqpUSY8cOeK39YbTz0U1vNoTam0BNmoGNcdXF/jRnFRcESnr3fNDRAoBdwK7c7JOkzsVKlSIpKQkvv/+ewYMGOB2HGPO8XUr3PbM3nPoSiBZRLYCX+A5BminAyNUo0aN6NmzJ6+88gpr1651O44xQAAvhFbVrapaR1VrquqNqjoiUNsyucPo0aO5+uqriYuL4/jx427HMcbnhdAfef98MXhxTDgrUqQIs2bNYu/evQwdOtTtOMb43AO8UkQaAq1EpI6I1E2/BCugCS9NmzYlPj6eCRMm8PnndlWUcZevC6GHAIOBCsCEC95ToGmgQpnwlpCQwLJly4iLi2Pz5s0UKFDA7UgmQvk6CfKWqt4DJKhqzAWLFT+TbSVKlGD69Ons3LmTkSNHuh3HRLBLngRR1ZEi0kpExnkXR3eBGONLixYt6NixI2PHjmXLFhtb17jjkgVQRMYATwE7vctTIjI60MFM+EtMTKRs2bLExcWRmprqdhwTgZxcBtMSuFNVk1Q1Cc9cIbYXaHKsVKlSTJ06lZSUFBISEtyOYyKQ0+sAS6Z7XCIAOUyEat26Ne3atWP48OHs2LHD7TgmwjgpgGOALSIy1zsV5iZgVGBjmUgyefJkSpQoQWxsLKdPn3Y7jokgTk6C/B34K/A2sBhooKqLAh3MRI6yZcsyefJkvvjiCxITE92OYyKI01nhflDVJd7lcKBDmcjTrl07WrduzeDBg/nqq6/cjmMihE2KZEKCiDBlyhQKFy5M586dOXPmjNuRTASwAmhCxpVXXsnEiRP55JNPePXVV92OYyKAk+sAX3fymjH+8Oijj3LPPfcwYMAA9u3b53YcE+ac7AHekP6JiEQB9QITx0Q6EWHGjBnkzZuXLl26cPbsWbcjmTDmazisASJyBKgpIn94lyN4Jjh6L2gJTcSpUKEC48ePZ9WqVUyfPt3tOCaM+RoMYYyqFgNeUtXi3qWYqpZWVRvX3ARU586dueOOO+jfvz/79+93O44JU772AK/zPnzzwrEAbTxAE2giwsyZM1FVHn/88bSJtYzxK1/HAHt7/xyfwTIuwLmMITo6moSEBD788ENmz57tdhwThnwNiPqh98/OqvpNMMIYc6Fu3brx5ptv0qdPH5o3b85VV13ldiQTRnztAaYd53srGEGMyUiePHmYNWsWp0+fpmvXrtYVNn7lqwD+IiIrgatFZMmFS7ACGlO5cmVGjx7NBx98wGuvveZ2HBNGfHWBWwJ1gdfxHPczxjU9e/bkrbfe4umnn+bOO++kfPnybkcyYcDXZTCnVPUzoKGqrlbV1cBaYIv3sTFBkydPHpKSkjhx4gTdunWzrrDxCyd3grwsIsVFpAiwHdgpIv0CnMuYi1StWpVRo0bx/vvv88Ybb7gdx4QBJwWwuqr+AbQGPgCuBh4NZChjMvPUU0/RoEEDevbsya+//up2HJPLOSmA+UQkH54CuERVU/HMC2xM0EVFRTFnzhyOHz9OYmKidYVNjjgpgNOB74AiwBoRqQT8cakvichVIpIsIjtFZIeIPJWzqMZ4XHvttYwcOZJ169axcOFCt+OYXMzJkPiTVPUvqtpCPfYDMQ7WfRroo6rV8Qyp/6SIVM9hXmMA6N27N9dffz09evTgxx9/dDuOyaV83Qv8iPfP3hcuQK9Lrdg7jP5m7+MjwC7gL37KbSJcVFQUzz77LMeOHeOJJ56wrrDJFl97gEW8fxbLZHFMRKKBOsDnWY9oTMYqVarEiBEjeOedd6wrbLJFAv2bU0SKAquBUar6dgbvxwPxAOXKlasXiv+Qjx49StGiRd2O4Rfh1pZChQrRo0cPDh06xJw5cyhVqpTbsbIt3H42odSWmJiYTapa/6I3VNXngueylwl4psVckrZc6nve7+YDVgC9nXy+Xr16GoqSk5PdjuA34diWnTt3aoECBbRNmzZ69uxZd0PlQDj+bEIFsFEzqDlOzgK/i+cs8GTOHxLLJxERYDawS1UnONiOMdly/fXXW1fYZIuTAnhCPWeCk9V7S5w6uxWuEZ4LppuKSIp3aZGzuMZkrE+fPtxyyy306NGDw4dt6mrjjNNb4YaKSIOsjAitqutUVVS1pqrW9i7L/JDZmItERUUxd+5cjh07ZvcKG8ecFMAawOPAWGxEaBPCrrvuOkaNGsV7771n9wobR5wUwAeBa1T1dlWN8S5NAx3MmOx4+umnadiw4bkzw8b44qQAbgdKBjiHMX6Rdq/wyZMniY+Pt66w8clJASwJ7BaRFTYitMkNqlWrxpgxY/jnP//J3Llz3Y5jQpivEaHTDA14CmP8rGfPnrzzzjs8/fTT3HHHHTaZksmQk8EQVme0BCOcMdmVNoL0mTNn6Ny5s3WFTYacdIGNyZWuueYaxo0bx4cffsi0adPcjmNCkBVAE9a6du3KXXfdRd++fdm3b5/bcUyIcVwARSSfiNQRkcsDGcgYfxIRZs+eTb58+ejUqRNnzpxxO5IJIb7GA5wmIjd4H5cAvgReA7aISPsg5TMmxypUqMDkyZP55JNPmDDBbks3/+NrD7Cxqu7wPo4FvlbVGkA9oH/AkxnjR4888ght2rRh0KBBbN++3e04JpsOHjzI3r17/bY+XwXwVLrHd+IZFQZVtTvNTa4jIkyfPp2SJUvyyCOPcOrUqUt/ybhu//79vP7663Tp0oWqVatSoUIFhgwZ4rf1+7oO8DcRuRc4iGdkl84AIpIXKOS3BMYESdmyZZk5cyb3338/w4cPZ9SoUW5HMumcPXuW3bt3s27dOtasWcPatWv5/vvvAShZsiS33XYb3bt354477vDbNn0VwK7AJOAK4Ol0e37NgH/6LYExQdSqVSvi4uIYO3YsLVu2pGHDhm5HilgnTpxg48aNfPLJJ+eWtLmey5UrR+PGjenbty+33XYbNWrUIE8e/1+0kmkBVNWvgbszeH0FnlGejcmVEhMT+fjjj+nYsSMpKSkhNXR7uFJV/v3vf/PZZ5+xfv161q9fz+bNm0lNTQU8ty+2bt2aW2+9lUaNGlG1alU8YyoHVqYFUEQm+fqiql5yZjhjQlHx4sV57bXXuP322+nduzczZsxwO1LYOXbsGMnJyWzYsIHPP/+czz///NzoPAULFuSmm26id+/eNGzYkIYNG1KmTBlXcvrqAnfDMxLMP4BDQODLsTFB0rhxY/r378+LL77IvffeS6tWrdyOlGsdP36clJQUNm3axBdffMHGjRvZvXv3udsPq1SpQkxMDLfccgsNGjSgVq1a5MuXz+XUHr4K4JV4xgJsh2eS80XAW6r6WxByGRNwI0aMYMWKFXTp0oVt27ZRrlw5tyPlGj/99BPz589n/vz5pKSkcPbsWQCuuOIK6tevz1//+lfatWtH/fr1KV26tMtpM5fpUUVV/UVVp6lqDJ7rAEsCO0Xk0WCFMyaQ8ufPz4IFCzhy5AixsbE2YIIDGzdupH379pQvX57evXuTN29eBg4cyLvvvsu///1vDh06xPvvv0/Hjh1p3rx5SBc/cDAclnf+j/Z4rgX8ANgU6FDGBEv16tUZN24cPXr04NVXX6VHjx5uRwo5//3vf3nnnXeYN28ea9asoXjx4vTs2ZPOnTtzww03uB0vR3ydBBkBtAR2AQuBAap6OljBjAmW7t27s2zZMvr27UuTJk248cYb3Y7kul9//ZUlS5awePFiVqxYQWpqKpUrV2b8+PF06dKF4sWLux3RL3ztAQ4CvgVqeZfR3tPSAqiq1gx8PGMCT0RISkqiZs2aPPzww3z++ecUKhR51/p/++23vP322yxdupS1a9dy5swZKlasSM+ePWnfvj316tULyqUpweSrAF4dtBTGuKxcuXLMmzePe+65h379+vHKK6+4HSkoVJXk5GQmTpzI0qVLUVVq1KhB//79eeCBB8Ky6KXn60Lo/Rm9LiJ58BwTzPB9Y3Kru+++m969ezNhwgTuvPNO7r//frcjBczBgwd57bXXmDNnDnv27KFs2bIMGjSIuLg4oqOj3Y4XNL6OARYHngT+AiwBPgR6AH3wDI21IBgBjQmm0aNHs2rVKuLi4khJSQmruUSOHDnCG2+8wcKFC1m9ejWqym233cbAgQNp164dBQsWdDti0Pm6ue514FpgG9AFSAb+BrRW1fD91WgiWoECBVi4cCGnTp3ioYceOnerVm63bNkyqlevTrdu3fjhhx8YOnQoe/bsYfXq1XTq1Ckiix/4LoDXqOpjqjodT5e3OtBcVVOCkswYl1StWpUZM2bw6aefMnjwYLfj5Mi+ffvo0KEDLVu2pHjx4qxdu5Zdu3YxdOhQqlSp4nY81/kqgOd+9anqGeCAqp5wumIRSRKR/4iIjT5pcp327dvTtWtXXnzxRZYtW+Z2nCzbt28fjz76KNWqVWPx4sUMGTKEzZs3c+utt4b1SY2s8lUAa4nIH97lCFAz7bGI/OFg3XPJYDQZY3KLxMREatWqxaOPPsp3333ndhzHlixZQp06dXj77bd55pln+Pbbbxk+fDgFChRwO1rI8XUrXJSqFvcuxVQ1b7rHl7wKUlXXAL/6Na0xQVSoUCHeeustzpw5w//93/9x4oTjDpArzp49y/Dhw7n//vupVq0au3btYty4cVx55ZVuRwtZNi2mMT5UqVKF119/nc2bN4f0bXInT56kQ4cODBs2jI4dO7J27VoqVqzodqyQJ4G8AVxEooGlqprpvUUiEg/EA5QrV67ewoULA5Ynu44ePRo2g2ZaW7Jn9uzZzJ8/n2eeeSZgQ2dltz1Hjhxh0KBBbN26lfj4eB566CHXj/OF2r+zmJiYTapa/6I3VDVgCxANbHf6+Xr16mkoSk5OdjuC31hbsuf06dPaokULzZs3b8C2m531bt++XatVq6b58+fXN954w/+hsinU/p0BGzWDmmNdYGMciIqK4o033qBKlSr87W9/49tvv3U7EosWLeKWW27h999/51//+hft29t03VkVsAIoIn8H1gPXisgBEekcqG0ZEwwlSpRgyZIlnD17lvvuu4/ffvvNlRyqypAhQ3jooYeoVasWmzdvpnHjxq5kye0CVgBVtb2qXqmq+VS1gqrODtS2jAmWqlWrsnjxYr7++mtat24d9DPDqampxMbGMnLkSDp37kxycjLly5cPaoZwYl1gY7IoJiaGefPmsXr1ajp27HhuOPhAS01N5f7772fevHmMGDGCmTNnkj9//qBsO1xdckRoY8zF2rdvz6FDh+jbty8lS5Zk2rRpAZm3No2q0r17dz744AOmT59OfHx8wLYVSawAGpNNvXv35tdff2X06NGICFOnTg1YERw/fjyzZs1i0KBBVvz8yAqgMdkkIrzwwguoKmPGjEFVmTp1KlFRUX7dzvz58+nfvz9t27Zl+PDhfl13pLMCaEwOiAijRo1CRBg9ejQHDx7k73//u1/mzFBVRo0axeDBg2nSpAlz584NaDc7EtnfpjE5lFYEp02bxooVK2jUqBH79u3L0TqPHTtGbGwsgwcP5pFHHmH58uUROU9JoFkBNMZPunbtyvLlyzlw4AA1a9ZkwoQJnD6d9YkU169fT+3atZk3bx5Dhw7ltddes5FcAsQKoDF+dMcdd7Bt2zaaNWtGnz59qFu3LsOGDWPVqlWcPHny3OcWbFtA9MRo8gzPQ/TEaBZsW8DBgwfp1q0bt956K6mpqSQnJzNs2DDX7+sNZ3YM0Bg/q1ChAu+99x5vvvkmCQkJjBw5kuHDh1O8eHHuu+8+SjcpzYwfZnDirOci6v2/76fTW504884ZCnxdgCeeeILRo0eHzdy7ocwKoDEBICK0bduWtm3b8ttvv7FmzRree+893n33XX4t+yuUPP/zZ/KcoXCrwuzrsY8rrrjClcyRyLrAxgRYyZIladWqFbNnz+bw4cNIyYy7tH/m+9OKX5BZATQmiPLly0fFEhkPVHp5gcuDnMZYATQmyEY1G0XhfIXPe61wvsJ0ubqLS4kilx0DNCbIOtToAMDAjwby/e/fU7FERUY1G8VffvmLy8kijxVAY1zQoUaHc4UwzapVq9wJE8GsC2yMiVhWAI0xEcsKoDEmYlkBNMZELCuAxpiIZQXQGBOxrAAaYyKWeCZNDw0i8hOw3+0cGSgD/Ox2CD+xtoSucGpPqLWlkqqWvfDFkCqAoUpENqpqfbdz+IO1JXSFU3tyS1usC2yMiVhWAI0xEcsKoDMz3A7gR9aW0BVO7ckVbbFjgMaYiGV7gMaYiGUFMB0RKSgiG0TkSxHZISLDva+LiIwSka9FZJeI9HI7qxM+2tNMRDaLSIqIrBORKm5ndUpEokRki4gs9T6/WkQ+F5G9IrJIRPK7ndGpDNqyQES+EpHtIpIkIvnczujUhW1J9/okETnqVq5LsQJ4vpNAU1WtBdQG7haRvwKPAVcB16nq9cBC1xJmTWbtmQp0UNXawBvAINcSZt1TwK50z18EElW1CvBfoLMrqbLnwrYsAK4DagCFgNw0RPSFbUFE6gOXuRPHGSuA6ahH2m+rfN5FgSeAEap61vu5/7gUMUt8tEeBtDkXSwCHXIiXZSJSAWgJzPI+F6Ap8Jb3I/OA1q6Ey6IL2wKgqsu8PzMFNgAV3MqXFRm1RUSigJeA/m7lcsIK4AW8u/IpwH+AD1X1c6Ay0E5ENorIByJS1dWQWZBJe7oAy0TkAPAoMNbFiFkxEc9/qLPe56WB31T1tPf5ASC3jCs/kfPbco636/sosDzImbJrIhe3pQewRFV/cCWRQ1YAL6CqZ7xdwwrAzSJyI1AAOOG9sn0mkORixCzJpD3PAC1UtQIwB5jgYkRHRORe4D+qusntLDnloC1TgDWqujaIsbIlo7aISHngQWCya8EcsjlBMqGqv4lIMnA3nj2Lt71vvYOnaOQq6dpzD1DLuycIsIjcsafRCGglIi2Agni68C8DJUUkr3cvsAJw0MWMTl3UFhGZr6qPiMhQoCzQ1dWEzmX0c9mB5/jzXs9RCgqLyF7vcdqQYnuA6YhIWREp6X1cCLgT2A28C8R4P3Y78LUb+bIqk/bsAkqISDXvx9JeC2mqOkBVK6hqNPAQ8LGqdgCSgb95P9YJeM+liI5l0pZHRKQL0Bxon3a8OdRl0pbLVPUKVY32vn48FIsf2B7gha4E5nkP4OYB/qGqS0VkHbBARJ4BjpJ7zs5l1p7HgcUichbPmdM4N0Pm0LPAQhF5AdgCzHY5T05MwzMa0nrvntPbqjrC3Ujhze4EMcZELOsCG2MilhVAY0zEsgJojIlYVgCNMRHLCqAxJmJZATTGRCwrgMaYiGUXQhuTi4hIazwjrxQHZqvqSncT5W62B5hLiMhA76CmW70Dmd6Sxe/38g7muiCAGY9e8PwxEXklm+sqKSLdfbx/xvv3kLZEB2LgTREpJCKrvXfTOPl8HhH5Ld3z2iKiabceikhRETkoIsW8z6eJSCOneVT1XVV9HOgGtPOuI7+IrBER26HJIvsLywVEpAFwL1BXVU+KSBnA0cjH3jHzBOgO3KGqBwKX1K9K4sk8JZP3//SOcnOO9/Yxf4vDc0vaGScfVtWzF+Toime8xbTxFx8G3lfVI97nfwWezEauQcCr3m2eEpGP8BTEgP2CC0e2B5g7XAn8rKonAVT1Z1U95N3r2Z72IRHpKyLDvK9/JSKvAdvx3B97DfCB935mRORdEdnk3auMT7eOjt69zC9F5HXva4+IZ2j9FBGZ7nRvKL3M1pFZDjxjFFb2fv6lbGyvt3iGlt8uIk97X+sn3ukMRCRRRD72Pm7qY8+4A+kGWBCRN0XkFfFMJbBfRG4VkdfFM11C2n3Ix0SkiIgUBRrjGUGomPe9x/Hc84uIXA98nVZcnaxbPF4EPlDVzelyvuvNarJCVW0J8QUoCqTgGYVmCnC79/VoYHu6z/UFhnlfPwv8Nd173wFl0j0v5f2zEJ4iWRq4wbuNMmmfAa4H3gfyeV+bAnTMJOcZb8605XvgFV/ryChHRm27xLbe8b521PtnPWAbUMT7d7cDqINnb+tN72fW4hl1OR8wFOiawTbyA4cveG030Nv7+HngKzy/oPICh/GMHbkLuAKIBwYDI4H7vRnWp1tXbyAui+vuBWzCU0S7pftuFPCT2/9Wc9tiXeBcQFWPikg9PHsTMcAiEXkOWOXja/tV9TMf7/cSkTbex1cBVYGb8BSIn73b/VVEHsZTUL7wdu0K4RldOiPndUtF5DGgPtDMxzoyyvGLj9wZbusCt+Ipise8Od7G83c3FagnIsXxjFe32ZuvMZ7CcqEywG/p2lMQT9d8ovclxXMi4gfv+2eAU8AfePb44oA2eEZ3LoanOzw93fqbA7FZWbeqTgImXRhUVc+IyCkRKab/616bS7ACmEuop5u0ClglItvwjH33L84/jFEw3eNjma1LRJoAdwANVPW4iKy64LvnfRyYp6oDsps9s3VkMUeOqWqqiHyLZ5KrT4GteH6hVCHjMRH/vCDPDcBm/d9YfbXwFNW0eTEOqaqKyB945io5oKo/iMgRPHtyd+MZjRsRKQyUVNVDWVn3JZpYADhxyb8Ic44dA8wFRORaOX8ektp4xo37EbhcREqLSAE8J0qcKAH811t0rsPTNQT4GHhQREp7t1sK+Aj4m4hcnvaaiFTKYhMyW0dmOQCO8L/jZlm1FmgtIoVFpAievbC16d7rC6zxPu4GbMmouKjqf4Eo794ZeGZr+zLdR2riKaLgKVhpj//A072dlq4tXYB3VfVP72sxeAZzTeN03Rny/sx+VtVUX58z57MCmDsUxTOw6U4R2QpUB4Z5/7GPwHMs60M8x5CcWA7kFZFdeE42fAagqjuAUcBqEfkSmKCqO/GccVzp3faHePZmHPOxjgxzeL/zC/CJ9yRGlk6CqOfkwFw8fy+fA7NUdYv37bXeba9X1R/x7DH5mntjJZ4uNXiKVAqc67IW8hZJOL9g/YFnr/cj7/MjQDX+VxDBMzVB+qkInK47MzHAPy/xGXMBGxDVGB9EpC7wjKo+6uf1bgZu8dcem/c453OqmiumawgVtgdojA/evcnk7Fz6c4n11vVj8cuPp3ttxS+LbA/QGBOxbA/QGBOxrAAaYyKWFUBjTMSyAmiMiVhWAI0xEcsKoDEmYlkBNMZELCuAxpiI9f9f0QSWDNOV3gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "T_ext = calc['T_K_calc']\n",
    "P_ext = calc['P_kbar_calc'] / 10.0\n",
    "\n",
    "#Solving for SHF values between 35 (SHF_start) and 45 (SHF_end) for the function with increments of 0.1 (SHF increment)\n",
    "#to minimise the RMS misfit function find the best fitting geotherm. \n",
    "#Adiabat=True adds adiabat to the end of the conductive geotherm.\n",
    "#Kinked=True adds kinked geotherm at the temperature BDL_T, parallel to the D-G transition.\n",
    "#max_depth=maximum depth until the geotherm is calculated.\n",
    "#plot_solution=True returns a graph of minimisation process.\n",
    "shf_solution, T_solution, depth_solution, p_solution, misfit_solution = pt.invert_generalised_mantle_geotherm(P_sample = P_ext, T_sample = T_ext, std_P = 0.3, std_T = 50,\n",
    " SHF_start = 35, SHF_end=45, SHF_increment=0.1, max_depth=300, kinked=False, BDL_T = 170, adiabat = True,\n",
    " plot_solution = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b0186f43",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARYAAAJVCAYAAAAbadoQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABLp0lEQVR4nO3dd3gUVdvH8e/ZJJsAoTcjNdIRI01RRCkqIkpTfARFQNSAoUgJUkRQEaSDUhQQEVDAgo8VC6+CPKhUpfcqoQiElpCePe8fu4npJGF2Z8v9ua692J0zM/tjs7lzpp1RWmuEEMJIFrMDCCG8jxQWIYThpLAIIQwnhUUIYTgpLEIIw0lhEUIYzu0Ki1KqnVLqgFLqsFJqpNl5hBAFp9zpPBallB9wEHgQiAK2AN211ntNDSaEKBB367HcCRzWWh/VWicBK4FOJmcSQhSQuxWWSsDJDK+jHNOEEB7E3+wABaWUCgfCHS+bFC9RgvIVKgDKxFSgtEYrczNk5E55JEvO3CkLFCzP0cOHLmity+fW7m6F5RRQJcPryo5p6bTWC4AFAKXLlNGXL12iwb33ETF5Gv4BAa5LmkXgscMkhtY07f2zcqc8kiVn7pQFCpane50aJ/Jqd7dNoS1ALaVUqFLKCnQDvs5t5jJly9E98mX++O5bpvfvR2J8vMuCCiFy51aFRWudAgwAfgT2AZ9qrffktUzHF/ry/PgJ7Fj/K5Oef5a4mBhXRBVC5MGtCguA1nq11rq21rqG1npCfpa5/z/dGDjjbQ7t2M74nk9z9WK0s2MKIfLgdoWlsO5u/wiR8+Zz+ugRXn+qGxdOnzY7khA+y2sKC0DD+1oyatESLl84z2tP/Yczx46ZHUkIn+RVhQWgbtOmvLpsOcmJibz21H84tjfPXTRCCCfwusICUL1efV5b/inWoCDGP/MU+7duNTuSED7FKwsLQEhoKK8t/4TS5Svw1nO92P7rOrMjCeEzvLawAJQNuZlxy1dSqUYNpkX05Y/V35odSQif4NWFBaBEmbKMWfIxtRo1ZvbQwfz8yQqzIwnh9by+sAAULV6cUe8vpmHLVrw/dgxfL5hvdiQhvJpPFBYAa1AQQ+e8S/NHO7Bi+hRWTJuCO41FI4Q3cbeLEJ3KPyCA/lNnULR4Cb5eOJ+42FieHfsaFovP1FchXMKnCguAxWKhz7jXKRoczNcL55OcmEj4mxOx+PmZHU0Ir+FzhQVAKUW3YcMJCApi1ey3SUlK4sXJU/Hz98mPQwjD+exvklKKrgMGYbUGsmL6FJKTkxg4bSb+VqvZ0YTweD5bWNJ0DO9LQKCVpRPfZGZSEi+9PQdrYKDZsYTwaLLXEni417P0eW08f679hWkvhsuAUULcICksDg92f4p+b01m9x+/Mzn8OeJjY82OJITHksKSQcvHutJ/6nQObNvKW8/3Ji5WRqMTojCksGRxz6MdeWnWbI7u2sXk5/tIz0WIQpDCkoM72z7EwBlvc3jnDiaHP0fCtWtmRxLCo0hhyUWzh9oxcPosDm3/i8l9nychLs7sSEJ4DCksebjr4fYMmDaDA9u2MrXv83K0SIh8ksJyHXe3f5T+U6azb+sWpvZ7QYqLEPkghSUf7unQkRcnTWXvpo1M79+PpMREsyMJ4daksOTTvZ0603fiJHb9toGZAyNITpLiIkRupLAUQMvHuvL8+Als/3Udb780iJSkJLMjCeGWpLAU0P3/6cazY19n2y//x+xhg0lJTjY7khBuRwpLIbR9ugc9R49h808/Mm9EJKkpKWZHEsKt+PzVzYX1cK9nSU1J4eMpk7BY/Bjcr6/ZkYRwG1JYbsCjz71AamoqK6dPxRJ7lb5z3pPBooRANoVuWKfwfnQf9jL/W7uWuS8Pk80iIZAeiyE6hvfF/3I0yxYtAq3pP3WG9FyET5Nvv0G6PPEfbGXL8/GUSdhsmgHTZuAfEGB2LCFMIYXFQI8+9wLKYuGjSRNBawZMnynFRfgkKSwGe+TZ51BKseytCWhtY+CMt6W4CJ8jO2+doH3vPunnubwzZJCcRCd8jhQWJ3m417P0euVVtqz5iXkjIrGlppodSQiXkU0hJ2rXszdJCYmsmD6FwKAivPDmRLmdq/AJUlicrGN4XxIT4vli7mysRYLoPWYcSimzYwnhVFJYXKDrwJdIjI/juw8WERhUhO6RL0txEV5NCosLKKV4+uVRJMYn8M37CwgsUpTHBww0O5YQTiOFxUWUUjw79jWSExP4fPYsAosE8ehzL5gdSwinkMLiQhaLhfA33yIxPp6Pp0zCGlSEtk/3MDuWEIaTwuJiFj8/+k+dQXJiEh+Of43gUiVp/kgHs2MJYSg59mkC/4AABs18mzpNmjJvxHB2/f6b2ZGEMJQUFpNYg4KIfHcBN4fewoz+L3J09y6zIwlhGCksJipWogQj319McOlSTA5/jrMnjpsdSQhDSGExWZmKFRn1/ofYUm289dyzXD5/3uxIQtwwKSxu4OZbbmHEwkVcuXCeSS/0IS42xuxIQtwQKSxuombY7Qx+Zy5Rhw4yPaKf3BBNeDQpLG6k4X0t6TtxMns3beS9USPQWpsdSYhCkfNY3My9nTpz8Z+zrJw+lZuqVuOJl4aYHUmIApPC4oY6vtCXf06c4It5c6hQtSotuzxudiQhCkQKixtSStHntTc4fyqKha++QvmbK1G/2V1mxxIi32Qfi5vyDwhg8DtzualqNWYMjOD00aNmRxIi36SwuLFiJUrw8vz38fPzY3L4c1y9eNHsSELkixQWN1ehShUi313ApX/O8vbggTIwt/AIUlg8QK2GjXh+/AT2btrIR5Mnmh1HiOuSnbce4r7Oj3Fi/35WL15EtTr1aP3Ef8yOJESupMfiQZ6KfJnb7mnBotfHcmDbVrPjCJErKSwexM/fn0Ez36H8zTczY2AE0WdOmx1JiBxJYfEwwSVLEvnuApITEpk5aIBcUyTckhQWD1SpRk36TZrCkZ07WPaW7MwV7kcKi4e6s+1DPNrnedYs/4gNX39ldhwhMvHoo0KnT0XxxjNPZZp2V7v2tH26B4nx8UwOfy7bMi27PEbLx7py9eJFZr00IFv7g92f4u72jxJ95jRzX47M1v7Is8/RpM39nD56lPfHjUmfbomPx1akCF1e7M9tze/h+L69LJ34Zrbluw0ZRu3GTTj45zZWzpyerb3n6DFUr1efXb//xn/fnZut/fnX3+TmW25h2y8/c3jXTooUL868EZF8v2QxgUWL0n/KNMqG3MyGX9fx/dix2ZYf/PYcSpQpw69ffM6v//0iW/uIBYsILFKEnz7+iI0/rM7WPnbZcgC+XbSQP9etzdRmDQxk5PuLAfhi7mx2b/wj02dTvFQphsyeB8CK6VM5tP2vTMuXqXgTA6bNAGDJhPGc2L8vU3tI9eq8MN7eQ1v46mjOHD+eqb1a3Xr0euVVAOZEDuXiP2cztddq2Ijej3UBYObACGIuX87U3uCuu3msv/1+T5Oef5akxMybmY1btU6/ZUvW7x0U/LuX9rmkKex3L82Nfvee79WLm0Nr5uu7dz0eXVgAkm22TK+vpaQQnZhIUmJitjaAWEd7bFJSju0xyfb2S7m2JxOdmMjlLO1WNMk2G1eTk4hOTORKLstfcSx/JTk55/Yk+/JXk3Ne/nJSEoGJicQkJ5OiNTfdUoO/9+zm5KGDVKl/K5eSkiAxEZtN57j8paQkkhMTiU1JybH9YmIiVouFa7m0Rzt+2XJqVzZbentcamp6e9pnk5j6b3t8avblk2yp6e0JGZZPk5Bh+YRUWw7t/y6fZMu+fHxqCimOjIk5LB+XYflkW/b2tO9WWntWBf3upX0uaQr73Utzo9+9VMdnk5/v3nVprT32Ubt2be0u1q5da9p7//rrr9rPz08/88wzbpEnK8mSM3fKonXB8gBbdR6/m7KPxQvcd999vPzyyyxbtoydO3eaHUcIKSzeYvjw4ZQsWZJx48aZHUUIKSzeonTp0gwdOpQvv/ySbdu2mR1H+DgpLF5k8ODBlClThrE5HA0SwpWksHiREiVKMHz4cFavXs2ePXvMjiN8mBQWLzNgwADKly/P4sWLzY4ifJgUFi8THBzMyJEj2bZtG+vXrzc7jvBRUli80IsvvkjZsmV59dVX5d5EwhRSWLxQkSJFePrpp1m/fj0//3z906+FMJoUFi/1yCOPUKVKFem1CFNIYfFSVquVMWPGsHHjRr7//nuz4wgfI4XFiz377LOEhoYyduxY6bUIl5LC4sUCAgIYO3Ys27Zt46uvZMwW4TpSWLxcjx49qF27NmPHjsWWw6XwQjiDFBYv5+/vz7hx49i1axeff/652XGEj5DC4gOefPJJ6tevz2uvvUZqaqrZcYQPMKWwKKWOK6V2KaW2K6W2OqaVUUqtUUodcvxb2oxs3sjPz4/XX3+dffv2sWLFCrPjCB9gZo+ltda6oda6qeP1SOBnrXUt4GfHa2GQxx57jNtvv53XX3+dlJQUs+MIL+dOm0KdgCWO50uAzuZF8T4Wi4U33niDw4cPs3TpUrPjCC9nVmHRwE9KqW1KqXDHtIpa6zOO52eBiuZE814dOnTgjjvu4I033iApKcnsOMKLKTNOnFJKVdJan1JKVQDWAAOBr7XWpTLMc0lrnW0/i6MQhQOUL1++yaeffuqi1HmLjY0lODjY7BjpcsuzefNmRowYwZAhQ+jYsaOpWcwgWXJXkDytW7felmE3RnZ5jbTtigfwGhAJHABCHNNCgAPXW1ZG6c9dbnlsNptu3ry5rlSpko6Pjzc1ixkkS+48epR+pVQxpVTxtOdAW2A38DXQyzFbL0BOFXUCpRTjx4/n1KlTLFiwwOw4wkuZsY+lIrBBKbUD2Ax8p7X+AZgEPKiUOgQ84HgtnKBNmza0atWKiRMnEhcXZ3Yc4YVcXli01ke11rc7HrdqrSc4pkdrre/XWtfSWj+gtb7o6my+ZPz48fzzzz/MmzfP7CjCC7nT4WbhQi1atKBt27ZMnjyZmJgYs+MILyOFxYe98cYbXLhwgdmzZ5sdRXgZKSw+rFmzZjzyyCNMmzaNK1eumB1HeBEpLD7ujTfe4NKlS8ycOdPsKMKLSGHxcY0bN6Zz587MmjWLy5cvmx1HeAkpLIKxY8dy5coV3nnnHbOjCC8hhUXQqFEjOnbsyMyZM2VfizCEFBYB2Hstly9fliNEwhBSWAQATZo04dFHH2XGjBlcvXrV7DjCw0lhEenGjRvHpUuXmDNnjtlRhIeTwiLSNW3alPbt2zN9+nQ5G1fcECksIpNx48Zx8eJF5s6da3YU4cGksIhM7rzzTtq1a8e0adOIjY01O47wUFJYRDbjxo0jOjparnwWhSaFRWRz11130bZtW6ZNm8a1a9fMjiM8kBQWkaNx48Zx/vx53n33XbOjCA8khUXkqHnz5jzwwANMnTpVRpkTBSaFReRq3LhxnDt3jvfee8/sKMLDSGERuWrRogVt2rRh8uTJ0msRBSKFReQprdcyf/58s6MIDyKFReTpvvvuo3Xr1kyZMoX4+Hiz4wgPIYVFXNe4ceM4e/as3IdI5JsUFnFdLVu2pGXLlkyePJmEhASz4wgPIIVF5Mu4ceM4c+YMCxcuNDuK8ABSWES+tGrVinvvvZdJkyZJr0VclxQWkS9KKcaNG8fp06dZtGiR2XGEm5PCIvKtTZs23HPPPUyaNInExESz4wg3JoVF5FtaryUqKooPPvjA7DjCjUlhEQXywAMPcPfdd/PWW29Jr0XkSgqLKJC0XsvJkyf58MMPzY4j3JQUFlFgbdu2pVmzZkycOJGkpCSz4wg3JIVFFFhar+Xvv/9myZIlZscRbkgKiyiUdu3acccddzBx4kSSk5PNjiPcjBQWUShpvZbjx4+zdOlSs+MINyOFRRRa+/btadq0KRMmTJBei8hECosoNKUUY8eO5dixY3z00UdmxxFuRAqLuCGPPvoojRs3ZsKECaSkpJgdR7gJKSzihiilGDNmDEeOHGHlypVmxxFuQgqLuGGdOnWiQYMGTJgwgdTUVLPjCDcghUXcMIvFwpgxY9i/fz+rVq0yO45wA1JYhCG6du1KnTp1ePPNN7HZbGbHESaTwiIM4efnxyuvvMKuXbv45ptvzI4jTCaFRRime/fu3HLLLYwfPx6ttdlxhImksAjD+Pv7M2rUKLZt28YPP/xgdhxhIikswlA9e/akatWq0mvxcVJYhKGsVisjRozgjz/+YO3atWbHESaRwiIM16dPH0JCQhg/frzZUYRJpLAIwwUFBTF8+HDWrVvHhg0bzI4jTCCFRThF3759KV++vPRafJQUFuEURYsWZdiwYfz000/s27fP7DjCxaSwCKeJiIigdOnSMqSCD5LCIpymePHiDB48mN9//53t27ebHUe4kBQW4VSDBg2iWLFivPnmm2ZHES4khUU4ValSpejSpQurVq1iz549ZscRLiKFRThd165dKVasGBMnTjQ7inARKSzC6UqWLElERAQrV67k0KFDZscRLiCFRbjE0KFDsVqt0mvxEVJYhEvcdNNNhIeHs2zZMo4dO2Z2HOFkUliEy7z88sv4+fkxefJks6MIJ5PCIlymUqVK9OnTh8WLFxMVFWV2HOFEUliES40YMQKbzcaUKVPMjiKcSAqLcKnq1avzzDPPsHDhQs6ePWt2HOEkUliEy40ePZqkpCSmT59udhThJFJYhMvVrFmT7t278+6773LhwgWz4wgnkMIiTPHKK68QFxfHzJkzzY4inEAKizBFvXr16Nq1K7Nnz+by5ctmxxEGk8IiTDN69GhiYmKYO3eu2VGEwaSwCNM0bNiQ9u3bM2vWLK5du2Z2HGEgKSzCVK+88goXLlxg4cKFZkcRBpLCIkzVvHlzWrZsybRp00hMTDQ7jjCIFBZhutGjR3Pq1CmWLVtmdhRhECkswnQPPvggTZs2ZdKkSaSkpJgdRxhACoswnVKK0aNHc+TIET777DOz4wgDSGERbqFTp07Ur1+fiRMnYrPZzI4jbpAUFuEWLBYLo0aNYvfu3Xz77bdmxxE3SAqLcBvdunUjNDSUCRMmoLU2O464AVJYhNvw9/dnxIgRbN68mV9++cXsOOIGSGERbqVXr16EhITIoNseTgqLcCtBQUEMGzaMX375hY0bN5odRxSSFBbhdvr27UuZMmWk1+LBpLAItxMcHMxLL73EN998w86dO82OIwrBaYVFKfWBUuqcUmp3hmlllFJrlFKHHP+WdkxXSql3lFKHlVI7lVKNnZVLeIaBAwcSHBzMpEmTzI4iCsGZPZYPgXZZpo0EftZa1wJ+drwGeBio5XiEA+86MZfwAKVLlyYiIoJPPvmEw4cPmx1HFJDTCovWej1wMcvkTsASx/MlQOcM05dqu41AKaVUiLOyCc8wZMgQAgIC5AZnHsjV+1gqaq3POJ6fBSo6nlcCTmaYL8oxTfiwm266ieeff54lS5Zw8uTJ6y8g3IZy5hmOSqnqwLda6waO15e11qUytF/SWpdWSn0LTNJab3BM/xkYobXemsM6w7FvLlG+fPkmn376qdPyF0RsbCzBwcFmx0jnTnluJMvZs2fp0aMHnTt3ZsCAAaZmMZo7ZYGC5WnduvU2rXXTXGfQWjvtAVQHdmd4fQAIcTwPAQ44ns8Huuc0X16P2rVra3exdu1asyNk4k55bjRLr169dJEiRfS5c+dMz2Ikd8qidcHyAFt1Hr+brt4U+hro5XjeC/gqw/SejqNDdwFX9L+bTMLHjRw5koSEBGbNmmV2FJFPzjzcvAL4A6ijlIpSSj0HTAIeVEodAh5wvAZYDRwFDgMLgQhn5RKep27dujz++OPMmTOHK1eumB1H5IMzjwp111qHaK0DtNaVtdaLtNbRWuv7tda1tNYPaK0vOubVWuv+WusaWuvbdA77VoRvGz16NFevXmXevHlmRxH5IGfeCo/QqFEjHn74YWbOnElcXJzZccR1SGERHmP06NGcP3+eDz74wOwo4jqksAiP0aJFC+655x6mTZtGcnKy2XFEHqSwCI8yYsQITpw4gbucvyRyJoVFeJRHHnmEW2+9lUmTJsnwlW5MCovwKBaLhREjRrB7925Wr15tdhyRCykswuN069aNqlWrypAKbkwKi/A4AQEBREZGsmHDBjZs2GB2HJEDKSzCIz333HOUK1dOhlRwU1JYhEcqWrQogwYN4ttvv2XXrl1mxxFZSGERHqt///4UK1aMKVOmmB1FZCGFRXisMmXK0LdvX1asWMGJEyfMjiMykMIiPNqQIUOwWCxMnz7d7CgiAykswqNVrlyZHj168P7773P+/Hmz4wgHKSzC4w0fPpyEhARmz55tdhThIIVFeLx69erRuXNn5syZQ0xMjNlxBFJYhJcYMWIEly5dYuHChWZHEUhhEV6iWbNmtG7dmunTp5OYmGh2HJ8nhUV4jZEjR3L69Gk+/vhjs6P4PCkswms8+OCDNGrUiClTpmCz2cyO49OksAivoZRi5MiRHDhwgK+++ur6CwinkcIivMrjjz9OjRo1ZCAok0lhEV7Fz8+Pl19+mc2bN7Nu3Tqz4/gsKSzC6/Ts2ZOKFSvKQFAmksIivE5QUBBDhgzhp59+4s8//zQ7jk+SwiK8Ur9+/ShRooT0WkwihUV4pZIlSxIREcGqVas4cuSI2XF8jhQW4bUGDRqEv78/M2bMMDuKz5HCIrxWSEgIPXr0YPHixVy4cMHsOD5FCovwapGRkcTHxzN37lyzo/gUKSzCq9WrV48OHTowe/Zs4uLizI7jM6SwCK83fPhwoqOj+fDDD82O4jOksAiv16JFC+666y6mT59Oamqq2XF8ghQW4fWUUgwfPpyjR4/yv//9z+w4PkEKi/AJnTp1ombNmnzyySdycaILSGERPsHPz49hw4axf/9+1q9fb3YcryeFRfiMXr16UapUKblzogtIYRE+o0iRInTp0oXVq1ezZ88es+N4NSkswqd06tSJokWLMm3aNLOjeDUpLMKnlCxZkj59+vDxxx9z6tQps+N4LSkswucMHTqU1NRU3n77bbOjeC0pLMLnhIaG8sQTTzB//nyuXLlidhyvJIVF+KThw4dz9epVFixYYHYUrySFRfikJk2a0KZNG95++22SkpLMjuN1pLAInzV8+HBOnTrFihUrzI7idaSwCJ/10EMPcdtttzF16lQ5zd9gUliEz1JKERkZyZ49e/j+++/NjuNVpLAIn9atWzcqV67M1KlTzY7iVaSwCJ9mtVoZPHgw69atY8uWLWbH8RpSWITPe+GFFyhRooT0WgwkhUX4vBIlStC3b19WrVrF8ePHzY7jFaSwCIH9HkQWi0VO8zeIFBYhgMqVK/Pkk0/y/vvvc/nyZbPjeDwpLEI4DB06lNjYWN5//32zo3g8KSxCODRu3JhWrVrx9ttvk5ycbHYcjyaFRYgMhg0bRlRUFJ999pnZUTyaFBYhMmjfvj116tRh+vTpcpr/DZDCIkQGFouFIUOG8Oeff8po/jdACosQWfTs2ZNy5coxffp0s6N4rFwLi1Lq6nUeMUqpg64MK4QrFClShIiICL755hsOHpSveGHk1WM5orUukcejOHDNVUGFcKWIiAgCAwOZOXOm2VE8Ul6F5fF8LJ+feYTwOBUrVqRHjx58+OGHXLhwwew4Hsc/twat9dGMr5VSJTLOr7W+mHUed5CcnExUVBQJCQkufd+SJUuyb98+l75nXtwpT1qWoKAgKleuTEBAgNmR8mXIkCEsWrSI9957jzFjxpgdx6PkWljSKKX6Aq8DCUDa8TcN3OLEXIUWFRVF8eLFqV69Okopl71vTEwMxYsXd9n7XY875YmJiSE4OJjo6GiioqIIDQ01O1K+3HrrrbRr1445c+YQGRlJUFCQ2ZE8Rn6OCkUCDbTW1bXWoY6HWxYVgISEBMqWLevSoiKuTylF2bJlXd6TvFHDhg3jn3/+Yfny5WZH8Sj5KSxHgDhnBzGSFBX35Ik/l/vvv5+wsDBmzJghJ8wVQH4Kyyjgd6XUfKXUO2kPZwfzZEopevTokf46JSWF8uXL8+ijj6ZP+/LLLwkLC6NevXrcdtttfPnll+ltvXv35vPPP8+0zuDg4PTne/bsoU2bNtSpU4datWoxfvx4+dI7iVKKoUOHsmfPHn766Sez43iM/BSW+cAvwEZgW4aHyEWxYsXYvXs38fHxAKxZs4ZKlSqlt+/YsYPIyEi++uor9u3bx9dff01kZCQ7d+687rrj4+Pp2LEjI0eO5MCBA+zYsYPff/+defPmOe3/4+u6d+9OSEiInDBXAPkpLAFa66Fa68Va6yVpD6cn83Dt27fnu+++A2DFihV07949vW3atGmMHj06fSdmaGgoo0aNytfQiMuXL+eee+6hbdu2ABQtWpQ5c+YwadIkJ/wvBNjHxR04cCBr1qxh165dZsfxCPkpLN8rpcKVUiFKqTJpD6cn83DdunVj5cqVJCQksHPnTpo1a5betmfPHpo0aZJp/qZNm7Jnz57018OHD6dhw4bpj7yWrVGjBrGxsVy9etU5/xlB3759KVq0KDNmzDA7ike47uFmIO1P7agM09z2cHNWrVq1yjbtP//5DxEREcTFxdG+ffts7b1796Z3795cuHCBrl27Zmpbt25dvt43LCyM48ePs2LFihzf43qmTp2a6b0z7mMRrlemTBmeffZZFixYwMSJEwkJCTE7klvLT4+lRobDzKFa61CgvrODeYOOHTsSGRmZaTMIoH79+mzblnk31bZt27j11luvu86clj169CjBwcGUKFHixkOLXA0ePJiUlBTmzp1rdhT3p7XO8wF8kOV1MeDn6y3nikft2rV1Vnv37s02zRWuXr2a/rxYsWJaa61Pnjyp3377ba211mvXrtWPPPKI1lrrv/76S9esWVMfO3ZMa631sWPHdI0aNfRff/2ltda6V69e+rPPPsu0/rR1xsXF6dDQUL1mzZr014888oh+5513cs1jtoxZzPr5pFm7du0NLd+lSxddpkwZHRsba3oWoxUkD7BV5/G7mZ8eyyml1DwApVRpYA3wkZPqnFepXLkygwYNyja9YcOGTJ48mQ4dOlC3bl06dOjAlClTMu1LyU2RIkX46quvePPNN6lTpw633XYbd9xxBwMGDHDC/0BkNXToUC5evMjSpUvNjuLe8qo6aQ9gCvAesAV4PD/LuOLhrj0Wd+BOebypx2Kz2XTTpk113bp1dWpqqqlZjOaSHotS6rG0B7AJuAv4C9COaUL4HKUUgwcPZv/+/XLCXB7y2hTqkOHxKPaiEpDhtRA+6YknniAkJIRZs2aZHcVt5TVswrOuDCKEp7BarfTv358xY8awd+9e6teXg6RZ5bUpFH69hfMzjxDeKDw8nKCgIN55Ry6by0leJ8iNVErlNXSWAl4CFhgbSQj3V758eXr06MHSpUuZMGECZcuWNTuSW8lrH8uvZN7PkvXxKPZDz0L4pJdeeon4+HgWLlxodhS3I/tYnCA4OJjY2Ngc2wYPHsxnn33GyZMnsVjsdf3DDz9k+PDhVKpUieTkZOrVq8fSpUspWrSoK2OLAmrQoAEPPPAAc+bMYdiwYR4z5KYrOO2+QkqpD5RS55RSuzNMe00pdUoptd3xaJ+hbZRS6rBS6oBS6iFn5TKTzWbjv//9L1WqVOHXX3/N1Pbkk0+yfft29uzZg9Vq5ZNPPjEppSiIwYMHc+rUKVatWmV2FLfizBuWfQi0y2H6TK11Q8djNYBSqj7QDbjVscw8pZSfE7OZYt26ddx66628+OKLrFixIsd5UlJSuHbtGqVLl3ZxOlEYDz/8MLVq1ZJDz1k4rbBordcDF/M5eydgpdY6UWt9DDgM3OmsbGZJG5elS5cufPfddyQnJ6e3ffLJJzRs2JBKlSpx8eJFOnToYGJSkV8Wi4WXXnqJTZs2sXHjRrPjuI38jNIfiP3+QdXJfPuPNwr5ngOUUj2BrcAwrfUloBL2EerSRDmm3bDVnTplmxbaqRP1+vQhJS6On7JceQxQq1s3anXvTkJ0NL/06ZOprf1XXxUqR1JSEqtXr2bGjBkUL16cZs2a8eOPP6YPV/nkk08yZ84ctNb079+fqVOnMnLkyEK9l3CtXr168corrzBr1ixWrlxpdhy3kJ/xWL4CrmAfjjLxBt/vXWA89vFcxgPTgT55LpGF49yZcLAf8ss6PkrJkiWJiYlJf52ampptHQkJCcTExJASH59ne2JsbLb2jOvOKDU1NVNb1vm+//57Ll++TIMGDQCIi4vD39+fli1bkpCQQFJSUvoy999/P/Pnz6d///65fQzXlTWPmTJmSUhIyPeYNs4QGxvrlPdv164dn332GY899hgVKlQwNUthGZonrwuJ7Ncasft68+SxbPXcls/Yhn0QqVEZ2n4E7r7e+t31IsS0IQ4y6t69u16+fHn669jYWF2+fHl97do1vXjxYt2/f//0ttGjR+sBAwYYlsds3nQRYm6OHz+uLRaLHjFihOlZCsvVwyb8rpS6zYgippTKOOxWFyDtiNHXQDelVKBSKhSoBWw24j3NEBcXR+XKldMfEydO5IcffuCRRx5Jn6dYsWK0aNGCb775Bvh3H0tYWBh//fUXr776qlnxRSFUq1aNxx57jAULFnDtmtzSPNdNIaXULuybLP7As0qpo9g3hRSgtdZhea1YKbUCaAWUU0pFAeOAVkqpho71Hgf6Yl/ZHqXUp8BeIAXor7XOvo3iIWw2W7Zpo0ePzjbtiy++SH/eu3dvZ0YSLjB48GA+//xzli1bRr9+/cyOY6q89rHc0BXMWuvse0VhUR7zTwAm3Mh7CmGm5s2b07RpU95++23Cw8PTT4D0Rbn+z7XWJ7TWJ4A3055nnOa6iEJ4Bhmr5V/5KamZRnh2nLjWJJd5hfBpMlaLXV7DJoxSSsUAYUqpq0qpGMfrc9gPQQshskgbq+XHH39k7969ZscxTV6bQm9prYsDU7XWJbTWxR2PslrrUbktJ4Svk7Fa8rcpNNox9u0MpdR0pVRnZ4cSwpNlHKslOjra7DimyE9hmQv0A3ZhP++kn1LK6+7YdPDgQUPWc/LkSUJDQ7l40X6Z1KVLlwgNDeX48eOGrB/g+PHj6Wfwbt26Nf0WI+vWreP3339Pn6937958/vnn112fn58fDRs25NZbb+X2229n+vTpOR4yz5ph+fLlN/C/8B42m42dO3fy3XffsXPnTmw2W/pYLQsW+OY4aPk5pb8NUM9xth1KqSXAnrwX8TwHDhygdu3aN7yeKlWq8OKLLzJy5EgWLFjAyJEjCQ8Pp3r16jceMgdNmzaladOmgL2wBAcH07x58wKto0iRImzfvh2Ac+fO8dRTT3H16lVef/31XJdJKyxPPfVUobN7A5vNxuzp04ja8gfVigSwLj6ZynfczcBhkbRp04Z3332X4cOH4++fn18175GfHsthoGqG11Uc07xC2l+bLVu2pP+1uVFDhgxh48aNzJo1iw0bNhAZGQnAkSNHaNeuHU2aNOHee+9l//79gL1nMWjQIJo3b84tt9yS3svQWjN8+HAaNGjAbbfdluMYLevWrePRRx/l+PHjvPfee8ycOZOGDRum91zWr1+fbb15qVChAgsWLEi/IPL48ePce++9NG7cmMaNG6evd+TIkfzvf/+jYcOGzJw5M9f5vN3u3buJ2vIHvRvX44Fba9O7cT1Obv6d3bt3M2jQIE6ePMlXhbxw1ZPlp4wWB/YppTZjP2P2TmCrUuprAK11Ryfmc6qMf22KxF5m2f5d6X9tbuTkpoCAAKZOnUq7du346aef0kcWCw8P57333qNWrVps2rSJiIgIfvnlFwDOnDnDhg0b2L9/Px07dqRr16588cUXbN++nR07dnDhwgXuuOMO7rvvvhzfs3r16vTr14/g4GAiIyOJiYlhxYoVOa73em655RZSU1M5d+4cFSpUYM2aNQQFBXHo0CG6d+/O1q1bmTRpEtOmTePbb78F7Jcx5DSftzl48CAHDhxIf71lyxaKxF4m6uTJ9GlFr11h1apVNGrUiAoVKjBu3Dhuu+02Q3rEniI/hWWs01OYJONfm6iTJ6lcpQqLHX9twsLyvGLhur7//ntCQkLYvXs3Dz74ILGxsfz+++888cQT6fMkJv57sXjnzp2xWCzUr1+ff/75B4ANGzbQvXt3/Pz8qFixIi1btmTLli0FypbTegsiOTmZAQMGsH37dvz8/HLdF5Xf+Txd7dq1MxWIatWq2f8gVamCn8VCqs1G3LkYej7+OGFhYRw6dIiXX36ZhIQEE1O73nULi9b6V6VUNaCW1vr/lFJFAH+ttXtck18IaX91Mv61OXf+PPDvX5sTJ05Qp06dQv2V2b59O2vWrGHjxo20aNGCbt26UaxYMUqVKpW+LyOrwMDA9OeO3VmGKMx6jx49ip+fHxUqVOD111+nYsWK7NixA5vNRlBQUI7LzJw5M1/zeZsGDRpQ+Y67Wbz5d6oXtXI8LokqdzZP37n+3HPPMW7cOGbPnu1Tg25ft7+vlHoB+ByY75hUGfjSiZmcrnbt2nTo0IGuXbsSH1yKylWqUKF8eSpXqUJcsZI8/vjjdOjQoVBFRWvNiy++yKxZs6hatSrDhw8nMjKSEiVKEBoaymeffZY+344dO/Jc17333ssnn3xCamoq58+fZ/369dx5Z+4D6xUvXvyGx2A5f/48/fr1Y8CAASiluHLlCiEhIVgsFpYtW5Y+Pk3W98ptPm9nsVgYOCySnq++wS1P9qbnq29k2pQuU6YMTz/9NB9//HH6kUJfkJ8dCf2Be4CrAFrrQ0D+RrJxc+l/bbbt5bdjUSzetjfTX5vCWLhwIVWrVuXBBx8EICIign379vHrr7/y8ccfs2jRIm6//XZuvfXW6+7U69KlC2FhYdx+++20adOGKVOmcNNNN+U6f4cOHfjvf/+baedtfsTHx6cfbn7ggQdo27Yt48aNS8+/ZMkSbr/9dvbv30+xYsUACAsLw8/Pj9tvv52ZM2fmOp8vsFgshIWF0b59e8LCwrLtnxs4cCDx8fEsWpTrNbjeJ6/BWhxd502Of/9y/OsP7Lzecq54GDHQU2pqqt6xY4ceO3as3rFjh05NTS3Q8mncaWAlrd0rjy8M9HQ9LVu21NWqVdMpKSmmZ8mNqwd6+lUpNRooopR6EPgM+MY5Zc710v7aNG3aNMe/NkIYYeDAgZw4cSJ9YC9vl5/fopHAeexn3vYFVgNjnBnKDHXq1DE7gvBinTp1okqVKsyePdvsKC5x3cKitbZh31kbobXuqrVe6OgKeRVfOsdAuJ6/vz8vvvgiv/zyC3v2eN2J69nkNWyCcty58AJwADiglDqvlHL781q8sO55BV//ubzwwgsEBgYyZ84cs6M4XV49liHYjwbdobUuo7UuAzQD7lFKDXFJukIICgoiOjra57/E7kZrTXR0tM+c35KTcuXK8dRTT7F06VIuX75sdhynyusEuWeAB7XWF9ImaK2PKqV6AD8BM50drjAqV65MVFQU5x0nvLlKQkKCW/3SuFOetCxBQUFUrlzZ7DimGjhwIIsXL+aDDz6gcePGZsdxmrwKS0DGopJGa31eKRXgxEw3JCAggNDQUJe/77p162jUqJHL3zc37pTHnbKYrVGjRtxzzz3MnTvXq4dUyGtTKKmQbUKIPAwaNIijR4+yebPH3jrruvIqLLc7xrrN+ogBDLmBmRC+qEuXLlSqVCnTfaW8TV5j3vpp+1i3WR/FtdZuuykkhLsLCAigX79+bN26NX1MHm8jp5kKYYLw8HACAgK89tCzFBYhTFChQgVatWrFkiVLuHr1qtlxDCeFRQiTPPbYY8TGxvLhhx+aHcVwUliEMEndunW56667mDNnjiFjLbsTKSxCmGjgwIEcOnTI6+71LIVFCBN17dqVm266yeuuepbCIoSJrFYrL7zwAt9//z3Hjh0zO45hpLAIYbLw8HAsFgvvvfee2VEMI4VFCJNVrlyZTp06sWjRIq+5TYgUFiHcQP/+/YmOjubTTz81O4ohpLAI4QZat25N3bp1mTt3rtlRDCGFRQg3oJQiIiKCzZs3e8WtaaWwCOEmevbsSbFixZg3b57ZUW6YFBYh3ETJkiXp0aMHK1as8Pi7JkphEcKNREREkJCQwOLFi82OckOksAjhRsLCwmjRogXvvvuuR18/JIVFCDcTERHBkSNHPPr6ISksQriZxx9/nAoVKnj0TlwpLEK4mbTrh7799luOHz9udpxCkcIihBvq27cvSinmz59vdpRCkcIihBuqUqUKHTt25P333/fI64eksAjhpvr378+FCxf4/PPPzY5SYFJYhHBTbdq0oXbt2rleP2Sz2di5cyffffcdO3fudKvD03ndYlUIYSKLxUJERASDBw/mzz//zHSvZ5vNxuzp04ja8gfVigSwLj6ZynfczcBhkVgs5vcXzE8ghMhVr169KFq0aLZDz7t37yZqyx/0alSPxiHlaFOuGHt/+ZGdO3ealDQz6bEI4cZKlSrF008/zdKlS7n//vsJDg4GYMuWLQTFXOKPDf8j5eplivv7wdlLjB4+nL4DBqT3WurUqUPt2rVdnlsKixBuLiIigoULF3L27FmGDBkCQLVq1Zi35XdCk67RoMrN2LRmQ0wy5UkhNDSUsLAwUzNLYRHCzTVs2JDmzZvzzjvvUKNGDZRS2Gw2olIVZ4+c5sSVGE7EJVGsVGmq+mtWrVrFiRMn0pc3o9cihUUID9CvXz969uxJ8eLFad26NWDvtcx6eSg33VyOxqVLckv5ciz5ax89H3/c9B6L7LwVwgN07dqV0qVLZxrJPywsjEYPtuPItSROXrzCkr/2UeXO5jRo0MDEpHbSYxHCAxQpUoTevXsze/Zs/vnnHypWrIjFYmHgsEj7EaKoKO6rXJkGDRrI4WYhRP717duXlJQUPvjgg/RpFouFsLAw2rdvT1hYmFsUFZDCIoTHqFOnDq1bt2bBggVudZZtTqSwCOFB+vXrx/Hjx/nhhx/c9nR+kH0sQniUzp07U758eV4d8TJt6tVyy9P5QXosQngUq9VKhw4dSPr7CE/UqcoDt9amd+N6nNz8O7t37zY7XjrpsQjh5g4ePMiBAwfSXxctWpRqRaxs2byZmjVr2qddu+IWJ8alkcIihJurXbt2pgJRrVo1Xlr9JadOn6Z1mzZorYk7F+MWJ8alkcIihIdp0KABte+7n69XfYLfut+IDy7lNifGpZF9LEJ4GIvFwuz5CzgXVJxVJ85R7b77CalW3a2OEEmPRQgP5O/vT7NGDTn+x/84snwRa6IvU7RkKSpXr06VO5ubfoRIeixCeKDdu3dTPdCP+8sHE2qL5+XGtSmtk2l+czm3OEIkPRYhPEDakaFTp05x6tQpDh06RMXofwguVoyka7EcOnmSoIQ4Nmz7k1QNb731FjVr1qRSpUpUqlTJ5UeIpLAI4QGyHhnauXMnS98cS8WL/lw7E4VfSjIHrsRRqfRNlLjpJkaNGmXqESLZFBLCAzVo0AD/yqGsPnqKUwkpfBkVTUhwUc6cPYN/1VtMP0IkPRYhPETGE+VsNhvbd+4kWNsgIIAQq41qZYpTJqAo1mLBfPfdd+nLyQhyQohc1a5dm5o1a7J7927Wr1+P9VoM5cqUpustIRw9cphLSYlcLVGRZ554wvQT5aSwCOHm0noqNpuN77/6kvijBymVHI/l3FkOYuGDvfGUs6Xw14WLBJWsyrFjx9JP7ZdR+oUQOUrbcbtz506Kx15mYLuWXLl8mf3btrDp4jXq16tL7NWr/P33/+h5X0s6depkdmQpLEK4q6wXH27ZsoUisZeJOnkSrSHVGoSOO8fG/YcJCA4mxlqE5cuX061bN9MuPkwjhUUIN5XTxYfL9u+icpUq+FksVK5Sha22TbR49DHuueceyi9bxowZM9i1axc1a9Y09cxbKSxCeIgGDRpQ+Y67Wbz5d6oXtXI8Lol6rdty9913M++dt7l2/AiPVCzBrJHD+PG7b4kY9JJp4+BKYRHCQ2Qdlb/FzTez9scfWPDKywSeOYHF38oR/Ckfcwnrtv/x3qiT1GvT1pTrhqSwCOFBDh8+nH7EZ/369Wz//hvuLx+MpXQxAv39OHrKRpOSQVQu6o812J//++4r5gQGERoa6tIjRHLmrRAeIuvO3AMHDlAqOZ6E+ASi4+KJjY+nbnAgV1PhXEwsCfHxlE6O5+DBg+nzf/PNN+mvnUl6LEJ4iBx35r55muaN6rFv9y7OnfybvbGJNClfkkvxcTSvU4eDf5/jufBwl58wJz0WITxU2s7cD//cy1n/IHYULUdqtdqcr1CVA9eSmfzDOtNGlpMeixAeKuvO3JaVK1O/fn327t3LsGHD+O2331gZ3teUo0JOe0elVBWl1Fql1F6l1B6l1EuO6WWUUmuUUocc/5Z2TFdKqXeUUoeVUjuVUo2dlU0Ib5H1Fqv+/v6EhYUxfvx44uPjWblypTm5nLjuFGCY1ro+cBfQXylVHxgJ/Ky1rgX87HgN8DBQy/EIB951YjYhvFqzZs0ICwtj/vz5aK1d/v5OKyxa6zNa6z8dz2OAfUAloBOwxDHbEqCz43knYKm22wiUUkqFOCufEN5MKUW/fv3466+/2LJli8vf3yUbX0qp6kAjYBNQUWt9xtF0FqjoeF4JOJlhsSjHNCFEITz99NMUK1aM+fPnu/y9nb7zVikVDKwCBmutryql0tu01lopVaB+mlIqHPumEuXLl2fdunUGpi282NhYt8kC7pVHsuTMFVlatWrF8uXL6dKlC8HBwa7Lo7V22gMIAH4EhmaYdgAIcTwPAQ44ns8Huuc0X26P2rVra3exdu1asyNk4k55JEvOXJFl69atGtCzZ882NA+wVefxu+m0Houyd00WAfu01jMyNH0N9AImOf79KsP0AUqplUAz4Ir+d5MpR8mnT7M6y9gToZ06Ua9PH1Li4vipe/dsy9Tq1o1a3buTEB3NL336ZGuv27s3t3TpQuypU6yPiMjW3iAigqoPPcSVw4f5bdiw9OmXL19mdalS3D50KJVatiR61y42jRmTbfkmr7xCxTvv5J/Nm9k2YUK29mZvvknZ227j1K+/smPGjGzt90yfTsmaNfn7xx/ZPW9etvb75s0juFIlYn/7jdUzZ2Zrb/PBBwSVLcuhFSs4lMMRg7YrVuBftCj7PviAY199la29vWParrlzOfnTT5na/IKCeOiTTwDYPn06p9evz/TZBJYuzf0ffgjA1vHjObd1a6bli918My3fte+z3/jKK1zMcguLEjVq0MLxmWwYOpSrR45kai/ToAF3OT7TX198kWunT2dqr9C0Kdx7LwA/9+5N4qVLmdpvvu8+Gjp+pj8++SSpCQmZ2qu0bctt/fsDZPveQcG/e2mfS5rCfvfS5Pbde6NSJaInTuTsHXdwU7NmuX73/Dp3hlat8vXdux5nbgrdAzwD7FJKbXdMG429oHyqlHoOOAH8x9G2GmgPHAbigGev+w5aE3cqc+05++tvEG8jNSkpWxvA6Z/Xk3LxGsnXruXYfuqntSSejibxyuUc209+9xPXDp8k/sL5TO02qz9xp87w91ffc3X3Ia6dOZ3j8idWfcPFLTuJ+ftEju3HPvkv59Zv4vKRwzm2H/n4M4qUK8/FA/tybD+8ZAWBJUuRfDWG1BzaDy76iIBixTj317Ycl98//0P8rFbObv4jx/Z9c98H4Nxvm7K1WwIC0tvPb9qa3p722SRdupLefuHPHdmWT4mJS2+/uGM3cWczt9sSktPbL+/ZT0L0hczhUv/Nd+XAYZKuXsnUfMGyA+rWY9/c94k5epyUuLhM7ec3bU1f/trfUdiSkzO1n/ttE/sIBMjxsynody/tc0lT2O9emty+e+WKFuPYuX/4vznv0WTrrly/e4GXLrNv7vv5+u5dj9ImHIoySu0aNfQPE6aYHQOAw4F+1ExMNTtGOnfKI1ly5qos1xLiaR4RzoNN72RaxEBD8tTo3nWb1rppbu1ySr8QXq5YUBGa1qnHwZN/u+w95ZR+IbyAzWbj4Mm/ORMdTUjZstSuUjXTqfxFAq0kZdm0cyYpLEJ4OJvNxrLV33DhxGGqFQ3gz7hkNlWryTPtO6QXl8AAK4lSWIQQ+XXw5N9cOHGYHnWr4WexkGqz8dH+Qxw8+Td1q1UHwBoQIIVFCJG7Y2dOczTDofSdRw4TnBzH2ejo9GnByfH8sHkTpy7Yj5ydu3SJuMSEbOtyFiksQniY0JCbCQ25Of11pXLl+O7bU9xUtmx6jyX2fAyP3Nksvceyad8eNu3b47KMUliE8HC1q1RlU7WafLT/ENWKWjkRl0S5arWoXaVq+jyBAVYSk5LRWpPxshpnkcIihIezWCw8076D/ajQxYs0KFMm21GhwIAAbNpGSmoqAf7O/7WXwiKEF7BYLNStVj190yerQGsAAEkpyS4pLHKCnBA+IDDACkBikmuODElhEcIHBAbYeyyuOuQshUUIH5DeY0lOcsn7SWERwgek7WORHosQwjDpm0JJ0mMRQhjE6u84KiQ9FiGEUWRTSAhhuLSdt9JjEUIYRg43CyEMZ5XDzUIIo8k+FiGE4f7dFJIeixDCIHKtkBDCcLLzVghhOH8/P5RSsikkhDCOUopAFw6oLYVFCB8RGOC6ewtJYRHCR0iPRQhhuEBrgOxjEUIYKzDASpIcbhZCGEk2hYQQhrPfZlU2hYQQBnLljeGlsAjhI+w7b6WwCCEMZPUPkPNYhBDGkp23QgjDBVqtsvNWCGGswIAAGTZBCGEs2RQSQhjOfhGibAoJIQwUGBBAUkoKNpvN6e8lhUUIH5E2oHZSivM3h6SwCOEjrOk3LUtx+ntJYRHCR7hypH4pLEL4CKsLB9SWwiKEj0jvsbjgXBYpLEL4iEAX3mZVCosQPsKV9xaSwiKEj/j3/s3SYxFCGMSVt1mVwiKEj0jbFJIT5IQQhvl3560UFiGEQdL3scimkBDCKFZ/2XkrhDCYHG4WQhgu/epm6bEIIYySvikk+1iEEEaxWCxY/f1lU0gIYSxrgGtG6pfCIoQPCQwIkIGehBDGst9mVXosQggDuerG8FJYhPAh9nsLSY9FCGEgq4tuWiaFRQgf4qrbrEphEcKHyKaQEMJw9tusSo9FCGEg+3ksUliEEAaSnbdCCMNJYRFCGE72sQghDGcNkKubhRAGC5Srm4UQRrMG+JNqs5GSmurU95HCIoQPSbsFiLP3s0hhEcKHuGpAbSksQvgQa9rdEKWwCCGMIj0WIYTh/r3NqnOPDElhEcKHWAP8AdkUEkIYyFU3hpfCIoQP8fh9LEqpKkqptUqpvUqpPUqplxzTX1NKnVJKbXc82mdYZpRS6rBS6oBS6iFnZRPCV7nqqJC/E9edAgzTWv+plCoObFNKrXG0zdRaT8s4s1KqPtANuBW4Gfg/pVRtrbVzTxEUwod4fI9Fa31Ga/2n43kMsA+olMcinYCVWutErfUx4DBwp7PyCeGLvOqokFKqOtAI2OSYNEAptVMp9YFSqrRjWiXgZIbFosi7EAkhCujfo0LOvRuiMzeFAFBKBQOrgMFa66tKqXeB8YB2/Dsd6FOA9YUD4QDly5fncKCf8aELIUHhNlnAvfJIlpyZkSW6WBAAUTol23sbmcephUUpFYC9qHystf4CQGv9T4b2hcC3jpengCoZFq/smJaJ1noBsACgdo0aumaie+yCORzoh7tkAffKI1lyZkaW8tpeOErEJ2Z7byPzOPOokAIWAfu01jMyTA/JMFsXYLfj+ddAN6VUoFIqFKgFbHZWPiF8kTccFboHeAbYpZTa7pg2GuiulGqIfVPoONAXQGu9Ryn1KbAX+xGl/nJESAhjuerMW6cVFq31BkDl0LQ6j2UmABOclUkIX+dn8cPfz887jgoJIdyH/d5Czj0qJIVFCB/jinFvpbAI4WNcMVK/FBYhfIzVBfcWksIihI8JdMHdEKWwCOFj7LdZlX0sQggDSY9FCGE4++FmKSxCCANZpccihDCabAoJIQwXKIebhRBGk00hIYThAuVwsxDCaFZ/OSokhDBYoFUKixDCYFb/AJJSUtBaO+09pLAI4WMCXTA8pRQWIXxM+ri3KVJYhBAGSb8bYpIUFiGEQaTHIoQw3L/3b3beuLdSWITwMf/eW8h5J8lJYRHCx/zbY5FNISGEQaz+crhZCGGwtE2hRNl5K4QwSvoJcnK4WQhhFDncLIQwnOy8FUIYzirXCgkhjCY9FiGE4eRwsxDCcIFW2XkrhDBYWo9Frm4WQhjGYrEQ4OcvPRYhhLGsAf6y81YIYSxn37RMCosQPkh6LEIIwzn73kJSWITwQc6+MbwUFiF8kDUgQI4KCSGMZQ0IkPNYhBDGkh6LEMJwgQGy81YIYTCr7LwVQhgtUDaFhBBGs/rLzlshhMECrVbpsQghjGX1l1P6hRAGs8pRISGE0dJO6ddaO2X9UliE8EFW/wC01qSkpjpl/VJYhPBBaePeOms/ixQWIXyQs0fql8IihA8KdPJtVqWwCOGDrE6+aZkUFiF8UGCAFZBNISGEgawB/oD0WIQQBpIeixDCcNJjEUIYTnosQgjDWf2lxyKEMJhVzmMRQhgtbVPIWYM9SWERwgel7byVHosQwjCy81YIYTg53CyEMJz0WIQQhvP388OiLNJjEUIYK9DqvHsLSWERwkc5895CUliE8FHOvBuiFBYhfJT9/s1JTlm3FBYhfFRgQABJySlOWbcUFiF8lP2mZdJjEUIYyH6bVemxCCEMFBjgvBvDS2ERwkdZA/zlzFshhLECA6xSWIQQxrIG+Msp/UIIY0mPRQhhOPsJclJYhBAGsp8g52GFRSkVpJTarJTaoZTao5R63TE9VCm1SSl1WCn1iVLK6pge6Hh92NFe3VnZhBCOixA9rbAAiUAbrfXtQEOgnVLqLmAyMFNrXRO4BDznmP854JJj+kzHfEIIJ7F64rAJ2i7W8TLA8dBAG+Bzx/QlQGfH806O1zja71dKKWflE8LXBfoHkJKaSqot1fB1O3Ufi1LKTym1HTgHrAGOAJe11mnnEUcBlRzPKwEnARztV4CyzswnhC9Lv7eQE07r9zd8jRlorVOBhkqpUsB/gbo3uk6lVDgQ7niZWKN71903uk6DlAMumB0iA3fKI1ly5hZZGvR+Ou1pQfJUy6vRqYUljdb6slJqLXA3UEop5e/olVQGTjlmOwVUAaKUUv5ASSA6h3UtABYAKKW2aq2buuL/cD3ulAXcK49kyZk7ZQFj8zjzqFB5R08FpVQR4EFgH7AW6OqYrRfwleP5147XONp/0VprZ+UTQjiPM3ssIcASpZQf9gL2qdb6W6XUXmClUupN4C9gkWP+RcAypdRh4CLQzYnZhBBO5LTCorXeCTTKYfpR4M4cpicATxTwbRYULp1TuFMWcK88kiVn7pQFDMyjZGtDCGE0OaVfCGE4jy0sSql2SqkDjksARrrg/aoopdYqpfY6LlF4yTH9NaXUKaXUdsejfYZlRjnyHVBKPWRwnuNKqV2O99zqmFZGKbVGKXXI8W9px3SllHrHkWWnUqqxgTnqZPi/b1dKXVVKDXbl56KU+kApdU4ptTvDtAJ/FkqpXo75DymleuX0XoXMMlUptd/xfv/NcFCjulIqPsNn9F6GZZo4fr6HHXkLfLJoLlkK/HMp1O+a1trjHoAf9pPtbgGswA6gvpPfMwRo7HheHDgI1AdeAyJzmL++I1cgEOrI62dgnuNAuSzTpgAjHc9HApMdz9sD3wMKuAvY5MSfy1ns5zi47HMB7gMaA7sL+1kAZYCjjn9LO56XNihLW8Df8XxyhizVM86XZT2bHfmUI+/DBmUp0M+lsL9rntpjuRM4rLU+qrVOAlZivyTAabTWZ7TWfzqex2A/dF4pj0U6ASu11ola62PAYXLYaW2wjJdFZL1cYqm224j9XKIQJ7z//cARrfWJ62Q09HPRWq/HfiQx6/sU5LN4CFijtb6otb6E/UzxdkZk0Vr/pP8923wj9vO3cuXIU0JrvVHbf+uXZsh/Q1nykNvPpVC/a55aWNJP/3fIeGmA0yn7ldeNgE2OSQMc3dwP0rrcLsiogZ+UUtuU/WxkgIpa6zOO52eBii7KkqYbsCLDazM+lzQF/SxclasP9h5ImlCl1F9KqV+VUvdmyBjlxCwF+bkU6nPx1MJiGqVUMLAKGKy1vgq8C9TAfgX3GWC6i6K00Fo3Bh4G+iul7svY6PhL57JDfso+/EVH4DPHJLM+l2xc/VnkRin1CpACfOyYdAaoqrVuBAwFliulSjg5hkt+Lp5aWNJO/0+T8dIAp1FKBWAvKh9rrb8A0Fr/o7VO1VrbgIX82613akat9SnHv+ewX4d1J/BP2iaO499zrsji8DDwp9b6H0cuUz6XDAr6WTg1l1KqN/Ao8LSj0OHY7Ih2PN+GfV9Gbcf7ZtxcMixLIX4uhfpcPLWwbAFqKfugUVbsXfCvnfmGjr3yi4B9WusZGaZn3FfRBUjbA/810E3ZB7AKBWph3yFnRJZiSqniac+x7xzcTebLIrJeLtHTcUTkLuBKhs0Eo3Qnw2aQGZ9LFgX9LH4E2iqlSjs2D9o6pt0wpVQ74GWgo9Y6LsP08sp+ZjpKqVuwfxZHHXmuKqXucnzvembIf6NZCvpzKdzvWkH3NLvLA/ve/YPYq/wrLni/Fti70zuB7Y5He2AZsMsx/WsgJMMyrzjyHaAQe/XzyHIL9r3zO4A9af9/7MNM/AwcAv4PKOOYroC5jiy7gKYGfzbFsF8wWjLDNJd9LtgL2hkgGfs+gOcK81lg3/9x2PF41sAsh7Hvp0j73rznmPdxx89vO/An0CHDeppi/6U/AszBcTKrAVkK/HMpzO+anHkrhDCcp24KCSHcmBQWIYThpLAIIQwnhUUIYTgpLEIIw0lhEUIYTgqL8ClKqYeVUhOUUvLddyL5cL2EUqpshjE2zmYZc8Nqdr6cKKVKKaUinLj+9PFOMky+D/vZpHdnmfcJZb+173ZlH29nnFKqiON1klKqnLNyeiMpLF5Cax2ttW6otW4IvIf9NrYNHY8ks3I5Tp3P7XtWCihwYbnOOrM64vhM0qQCT2Mf9iJtfb2AEcDjjnnvAC5qreMdr08XNKOvk8LiI5RSPZRSmx1/gecr+10qqyv7yGYfKqUOKqU+Vko9oJT6TdlHUbvTsWzafB8rpfYppT5XShW9znoPKKWWYj8tvYpS6ktlH+Jhj/p3mIdJQA3HslMdy2Uc7SxSKfVahgxZ15ntva/3OWitx2itn9BaX3SstwQwA/iP1jrKMU+c1nq2QR+9T5LC4gOUUvWAJ4F7HH+B0/5qA9TEful8XcfjKezXRUUCozOspg4wT2tdD7gKRFxnvbUc89+q7QM/9dFaN8F+DcwgpVRZ7CO7HXH0qobn47+Svk6gaB7vXRCdsY8id7QQy4pcuOROiMJ09wNNgC32i2Upgn0YgfXAMa31LgCl1B7gZ621Vkrtwj50YpqTWuvfHM8/AgYBCXms94S2j9CWZpBSqovjeRXsReJsAf8fGdeZ2/+poBpgvwhQGEgKi29QwBKt9ahME+0j4SVmmGTL8NpG5u9H1qtV9XXWey3D61bAA8DdWus4pdQ6ICiHnClk7kVnnedahuc5vnchXMNelISBZFPIN/wMdFVKVYD0EezzvKl3DqoqpdKOpDwFbCjAeksClxxFpS72QaIBYrAPTJ7mH6CC4whXIPaBkZz5fwL7MJFPKKUqOtYTqJR6oRDrERlIYfEBWuu9wBjsY+TuxD5QdEEH0z6AfQjMfdhHsX+3AOv9AfB3LDsJ+4DSaPvoab8ppXYrpaZqrZOBN7APMLQG2O/k/xNa683YR67/0bGe7UCFgq5HZCbjsYjrcmzafKu1bmB2loIwKrdS6jj2AaEuGJHLF0iPRXizVKBklhPk8i3tBDkgAPs+J5FP0mMRQhhOeixCCMNJYRFCGE4KixDCcFJYhBCGk8IihDCcFBYhhOGksAghDCeFRQhhuP8HHok2jbhzcr0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Making a plot of the solution.\n",
    "pt.mantle_geotherm_plot(T = T_solution, P = p_solution, Depth = depth_solution, \n",
    "                        plot_style = 'Depth', Temp_unit = 'Celsius', \n",
    "                        T_Sample = T_ext, P_Sample = P_ext, T_std = 50, \n",
    "                        P_std = 0.3, plot_type = 'show', max_depth = 250, \n",
    "                        moho = 38, lab = 202, leg = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1d360b70-c577-432e-b3b8-51232d6bf8de",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function mantle_geotherm_plot in module Thermobar.plotting:\n",
      "\n",
      "mantle_geotherm_plot(T, P, Depth, plot_style, Temp_unit, T_Sample, P_Sample, T_std, P_std, max_depth, plot_type, **kwargs)\n",
      "    A function to plot calculate geotherm alongside the thermobarometric\n",
      "    calculations.\n",
      "    \n",
      "    ###Parameters###\n",
      "    T: Temperature array of the geotherm.\n",
      "    \n",
      "    P: Pressure array of the geotherm.\n",
      "    \n",
      "    Depth: Depth array of the geotherm in meters.\n",
      "    \n",
      "    plot_style: String parameter for the y-axis of the geotherm plot 'Pressure' or 'Depth'.\n",
      "    \n",
      "    Temp_unit: String parameter for the temperature unit, 'Celsius' or 'Kelvin'.\n",
      "    \n",
      "    T_Sample: Array of temperature of the thermobarometric solutions.\n",
      "    \n",
      "    P_Sample: Array of pressure of the thermobarometric solutions in GPa.\n",
      "    \n",
      "    T_std: Standart deviation of thermobarometric temperature estimation. Could be array or a single value.\n",
      "    \n",
      "    P_std: Standart deviation of thermobarometric pressure estimation. Could be array or a single value.\n",
      "    \n",
      "    max_depth: Maximum depth to show the plot.\n",
      "    \n",
      "    leg: Boolean parameter to set existence of a legend.\n",
      "    \n",
      "    plot_type: 'show' or 'save' the figure.\n",
      "    \n",
      "    moho: moho depth in km.\n",
      "    \n",
      "    lab: lab depth in km.\n",
      "    \n",
      "    Depth_Sample: Array of depths of the thermobarometric solutions.\n",
      "    \n",
      "    filename_save: string parameter for filename to save the figures.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(pt.mantle_geotherm_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "901d6f02",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}