mrtommyb/Why mass and radius are not consistent with log(g).ipynb

## Why mass and radius are not consistent with log(g).ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "df4bd5f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "073eabae",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "true_mass = 0.2, 0.1 # mean and standard deviation of mass\n",
    "true_radius = 0.2, 0.1 # mean and standard deviation of radius [though we don't use the mean]\n",
    "lower_bound = 0.08 # the smallest a star can be\n",
    "upper_bound = 1 # random high number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "6a1297d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate a truncated normal - truncated to keep mass physical\n",
    "stellar_mass = stats.truncnorm.rvs((lower_bound - true_mass[0]) / true_mass[1],\n",
    "                             (upper_bound - true_mass[0]) / true_mass[1],\n",
    "                             loc=true_mass[0], scale=true_mass[1], size=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "f34687cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make stellar radius covariant with stellar mass, keeping the\n",
    "\n",
    "stellar_radius = stats.truncnorm.rvs((lower_bound - stellar_mass) / true_radius[1],\n",
    "                             (upper_bound - stellar_mass) / true_radius[1],\n",
    "                             loc=stellar_mass, scale=true_radius[1], size=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "7e88c465",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x163f47d60>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABFl0lEQVR4nO3dcXAV130v8K+QhFCphbCxVcDkDqPYMrFiA/IYAcUmphZDmNa4f5iYwdht3Cm0DsUMnQeDX2U57chNcezx1FLDODgTExNeaydOC6+gqU2MTWmfNaKOY4I9JJYoiBAckDCjChD7/lD2eu9qd+85u2f3nN37/cxobF3u3T179mrPb8/5nbNllmVZICIiItJknO4CEBERUWljMEJERERaMRghIiIirRiMEBERkVYMRoiIiEgrBiNERESkFYMRIiIi0orBCBEREWlVobsAIq5evYpTp07hmmuuQVlZme7iEBERkQDLsnDhwgVMmzYN48b593+kIhg5deoUZsyYobsYREREFMKJEydw4403+v57KoKRa665BsDowdTU1GguDREREYkYHBzEjBkz8u24n1QEI/bQTE1NDYMRIiKilCmWYsEEViIiItKKwQgRERFpxWCEiIiItGIwQkRERFoxGCEiIiKtGIwQERGRVgxGiIiISCsGI0RERKQVgxEiIiLSisEIERERacVghIiIiLRiMEJERERaMRghIqKSsfNwLxY+/QZ2Hu7VXRRyYDBCREQlo/PAcZw8P4TOA8d1F4UcGIwQUcniXXLpWbe4HtNrq7Fucb3uopBDmWVZlu5CFDM4OIhJkyZhYGAANTU1uotDRBmx8Ok3cPL8EKbXVuOdzffoLg5R5oi23+wZIaKSxbtkIjOwZ4SIiIhiwZ4RIiIiSgUGI0RERKQVgxEiIiLSisEIERERacVghIioCK5HQhQvBiNEREVw1U6ieDEYISIqguuREMWL64wQERFRLLjOCBEREaUCgxEiIiLSisEIERERacVghIiIiLRiMEJERERaMRghIiIirRiMEBERkVYMRoiIiEgrBiNERESkFYMRIiIi0ipUMNLR0YGZM2diwoQJaGpqwsGDBwPfPzw8jK1btyKXy6Gqqgr19fXYsWNHqAITERFRtlTIfmD37t3YsGEDOjo6sHDhQnzrW9/CsmXL8MEHH+Bzn/uc52ceeOAB/PKXv8S3v/1tfP7zn8eZM2dw5cqVyIUnIiKi9JN+UN68efMwd+5cdHZ25l+bNWsWVqxYgfb29jHv/9d//Vd85Stfwc9//nNce+21oQrJB+URERGlTywPyrt06RK6u7vR0tJS8HpLSwsOHTrk+Zkf/ehHuOOOO/CNb3wD06dPx80334xNmzZhaGjIdz/Dw8MYHBws+CEiIqJskhqmOXv2LEZGRlBXV1fwel1dHU6fPu35mZ///Od4++23MWHCBPzgBz/A2bNn8Wd/9mf49a9/7Zs30t7ejra2NpmiERERUUqFSmAtKysr+N2yrDGv2a5evYqysjJ873vfw5133okvf/nL+OY3v4nvfOc7vr0jW7ZswcDAQP7nxIkTYYpJREREKSAVjEyZMgXl5eVjekHOnDkzprfENnXqVEyfPh2TJk3KvzZr1ixYloX//u//9vxMVVUVampqCn6IiCi8nYd7sfDpN7DzcK/uohCNIRWMjB8/Hk1NTejq6ip4vaurCwsWLPD8zMKFC3Hq1Cl8+umn+dc+/PBDjBs3DjfeeGOIIhMRkazOA8dx8vwQOg8c110UojGkh2k2btyIF198ETt27MDRo0fx+OOPo6+vD2vXrgUwOsSyZs2a/PtXrVqF6667Dn/0R3+EDz74AG+99Rb+8i//En/8x3+M6upqdUdCRES+1i2ux/TaaqxbXK+7KERjSK8zsnLlSnzyySd46qmn0N/fj8bGRuzduxe5XA4A0N/fj76+vvz7f/u3fxtdXV342te+hjvuuAPXXXcdHnjgAfz1X/+1uqMgIqJAq5tzWN2c010MIk/S64zowHVGiIgoLXYe7kXngeNYt7i+5APAWNYZISIiomDMz5HHYISIiEgh5ufI4zANERERxYLDNERERJQKDEaIUowLWRFRFjAYIUoxJsoRURYwGCFKMSbKEVEWMIGViIiIYsEEViIiIkoFBiNERESkFYMRIiIi0orBCBEREWnFYISIiIi0YjBCREREWjEYISIiIq0YjBAREZFWDEaIiIhIKwYjREREpBWDESIiItKKwQgRERFpxWCEiIiItGIwQkRERFoxGCEiIiKtGIwQERGRVgxGiIiISCsGI0RERKQVgxEiA+w83IuFT7+BnYd7dReFiChxDEaIDNB54DhOnh9C54HjuotCRJQ4BiNEBli3uB7Ta6uxbnG97qIQESWuzLIsS3chihkcHMSkSZMwMDCAmpoa3cUhIioJOw/3ovPAcaxbXI/VzTndxaEUEm2/2TNCRESeOHxISWEwQkREnjh8SEnhMA0RERHFgsM0RERElAoMRoiIiEgrBiNEREQRceHCaBiMEBERRcSZR9EwGCEiIoqIM4+i4WwaIiLK40JnpBJn0xARkTQON5AODEaIiCiPww2kA4dpiIiIKBYcpiEiIqJUYDBCRKnHNR6I0o3BCBGlHpMuidKNwQgRpR6TLkkl9rQljwmsREREDguffgMnzw9hem013tl8j+7ipBoTWImIiEJgT1vy2DNCREREsYi1Z6SjowMzZ87EhAkT0NTUhIMHD/q+98CBAygrKxvz87Of/SzMromIiChjpIOR3bt3Y8OGDdi6dSt6enqwaNEiLFu2DH19fYGfO3bsGPr7+/M/N910U+hCE5EeTOwjojhIByPf/OY38dWvfhWPPvooZs2aheeeew4zZsxAZ2dn4OduuOEG/M7v/E7+p7y8PHShiUqdrqCAU2iJKA5SwcilS5fQ3d2NlpaWgtdbWlpw6NChwM/OmTMHU6dOxZIlS/Dmm28Gvnd4eBiDg4MFP0T0GR1Bwc7Dvbg4fAW11ZVM7CMipaSCkbNnz2JkZAR1dXUFr9fV1eH06dOen5k6dSq2b9+OV199Fa+99hoaGhqwZMkSvPXWW777aW9vx6RJk/I/M2bMkCkmUebpyPbvPHAc54cuY2JVBR8tT0RKVYT5UFlZWcHvlmWNec3W0NCAhoaG/O/z58/HiRMnsG3bNtx1112en9myZQs2btyY/31wcJABCZHD6uZc4gHBusX16DxwnL0iRKScVM/IlClTUF5ePqYX5MyZM2N6S4I0Nzfjo48+8v33qqoq1NTUFPwQkV6rm3N4Z/M97BUhT2lNbk5rubNGKhgZP348mpqa0NXVVfB6V1cXFixYILydnp4eTJ06VWbXRFQCSqFhyOoxpjW5Oa3lzhrp2TQbN27Eiy++iB07duDo0aN4/PHH0dfXh7Vr1wIYHWJZs2ZN/v3PPfccfvjDH+Kjjz7CT3/6U2zZsgWvvvoqHnvsMXVHQUSZUAoNQ1aPMa2rlqa13FkjnTOycuVKfPLJJ3jqqafQ39+PxsZG7N27F7ncaNdtf39/wZojly5dwqZNm3Dy5ElUV1fj1ltvxZ49e/DlL39Z3VEQUSaUQl5KVo9RRx6TCmktd9ZwOXgiIiKKBR+UR0SkQVZzQojixGCEiDyxUQ1HV04IzxelGYMRIvKU1UTLuOlKiOT5ojRjMEJEnjjLIBxd67HwfGVHKfZyMYGViIi02nm4Nz/DiDNbgIVPv4GT54cwvbYa72y+R3dxImECKxEZpxTv+Kg4DjEVKsVeLgYjRJQYZ6NTaoGJqcdrQrlKsfENUoqPXmAwQpRxJjQ2NmejU2p3w6Yeb9zlEvn+lWLjS4UYjBBlnEmNoLPRKbW7YVOPN+5ymfT9I3MxgZUo45gcOBbrJDms69Im2n4zGCGikrLzcC9aX38fIxYyMVuByGScTUNE5KHzwHGMWEB5GYwbMiEqVQxGUsSkREQikwX9rdg5Em33NXLYgMgQDEZShIlgZCrTAuWgvxXO3CAyD4ORFDE1G5/SxR04qAgkTAuUZf9WTAumiEoNE1iJSox7qWkVS0+nfcZElpbfJjIJE1iJyJO710BFj1vahz7Y60ikF3tGyFhpv9smInPweqIHe0Yo9UzLQyCiUWnMseH1xGwMRshY7DonWSY1kjrLEve+09iw83piNgYjZKy05yGQv7gaS5MaSZ1liXvfaWzYeT0xG4MRohglcXdsUm+AqLgaS5MaSZ1liXvfbNhJNSawEsUoiSmjaZyWymRCihO/X+ZgAiuRg67egyTujk3qDRBV7M46jb09ZA6ThutIDIMRKgm6Lk5xd2dn9Q6QjQlFkcYAvdQxGKGSkNWLU1Yb7bScL/bgmIk5LenDnBFKpaz2CMhiPeiVxnwdoiQxZ4QyLas9ArKyeAeYpt6GtPTgEJmOwQilEhuB7EpToJnFYJBIBwYjlEpsBLIr64Fmmnp+dGNdlQ4GI0SUqGINjB1oAshkQ5Smnh/dWFelg8EIUQJ4h/cZ0QYmqw1R1nt+VGJdlQ4GI0QJyGrDGoZoA5PVhkhmiLHUg1gOx5YOTu0lSkAcU3DjntbLacP6ceowpR2n9hIZJI674Si9LSL7YG+OfnbvUFNuckn3kFD2MRghiolMF7vzvaJBQJRhDJF9ZHWYJE3sILa79xwDQ8o0BiNEMZHpWXC+VzQIiDKeLrIPjtebg4EhZR1zRohiIpNzwfwMIsoi0fabwQgRlQwGffqw7ksTE1iJyGg6pq0yKVcf1j0FYTBCRFroaJyYe6EP656CMBghKqLUF56Ky7rF9aitrsTF4SuJ1W2pJuWa8B0u1bpPAxO+HwxGiIpg93I8VjfnMLGqAueHLo+pW9UXR1XbM+GiHQa/wxTEhO8HgxGiIti9HB+/3hGvi2OUQMC9vbDbMuGiHQa/wxTEhO8HgxGiIti97E1FL4Ff74jXxTFKIODeXthtmXDRDoPfYQpiwveDU3uJyFfQdExVz00RnfKpcmpoUtNMdx7uxbZ9xwAAm5Y2ZC4Y4HRdKobrjBBRZEEBBxui4uz6A5DJh93xQX5UDNcZISpxKoZRgma8mNC162Zagqldf7XVlakb2hGR1mErMk+onpGOjg783d/9Hfr7+3Hrrbfiueeew6JFi4p+7p133sHdd9+NxsZGHDlyRHh/7BkhkqfqrjUNd792L83F4Ss4P3TZ6LISlZLYekZ2796NDRs2YOvWrejp6cGiRYuwbNky9PX1BX5uYGAAa9aswZIlS2R3SUQhqLprTcPdr52QCsD4spIepvWaUSHpnpF58+Zh7ty56OzszL82a9YsrFixAu3t7b6f+8pXvoKbbroJ5eXl+OEPf8ieEQMxB6D0ZOWcp/E40ljmNEtDD18WxdIzcunSJXR3d6OlpaXg9ZaWFhw6dMj3cy+99BKOHz+O1tZWof0MDw9jcHCw4Ifil9Y1FCi8rJxzE/NXislK3adFGnr4SplUMHL27FmMjIygrq6u4PW6ujqcPn3a8zMfffQRNm/ejO9973uoqKgQ2k97ezsmTZqU/5kxY4ZMMSmkYn+s7ObMHl6g9WnKTUZ52eh/KX5pDFhLSajZNGVlZQW/W5Y15jUAGBkZwapVq9DW1oabb75ZePtbtmzBwMBA/ufEiRNhikmSiv2x8k4ue8JcoBmURrfzcC/2vHcKIxbQ3XtOd3GItBPrqviNKVOmoLy8fEwvyJkzZ8b0lgDAhQsX8O6776KnpwePPfYYAODq1auwLAsVFRXYv38/7rln7NhdVVUVqqqqZIpGCVi3uD4/xk2lyxmU+gUxznwI+zPMjfhM54HjGLGA8jLw74kIksHI+PHj0dTUhK6uLtx///3517u6unDfffeNeX9NTQ1+8pOfFLzW0dGBN954A//0T/+EmTNnhiw26bC6OcfGhISCUncvWrHgpdQ465B1QiQZjADAxo0b8dBDD+GOO+7A/PnzsX37dvT19WHt2rUARodYTp48ie9+97sYN24cGhsbCz5/ww03YMKECWNeJ4oTZy6oIxKUugMW9qgVYmBPVEg6GFm5ciU++eQTPPXUU+jv70djYyP27t2LXG70D6u/v7/omiOUDllqwEWGFkgdd2NbKnUe199Mlv4Wibzw2TTkK0vz8nkxN08Wz0lcfzNZ+luk0sJn01BkSU37TGJ2Rtqn9WVxBksWZ2fF9TfDKdiUdQxGyFdcDbi7YU1Do6Q7GEhDHclyNrC669d0aQ+miYphMEKJczesabjr0x0MuOsoC423s4HVXb+qZOU4iJLGYIQS525Y03DXpztgctdR1hq9uOs3qeBN9/eEKK2YwEqUQllM/oyTOwHU5PozuWxEspjASpRhaehNMom7x8LkniWTy0YUFwYj5CkLOQlENnfwpno4ReXfS5qHenjdoLAYjJAn3p3JSctFOKicaTkGFcL2LPnVUdi/F6/tyZbNpPPG6waFxWCEPKX57kwl0Qt9Wi7CQeVMyzHEReRc+9VR2L8XFXUus424A5eo1w2TAitKFoMR8sRHy48SvdCnJXgLKqepx5DU90rkXPvVUdieFhV1LrONuAPOqLlMpR4QlzLOpiFlsrhktT2zoSk3Gd2952KZ4cDZE8GS+l6Vwnkw/RhNLx/JE22/GYyQMlm+kIg0iGGP3952eRnQdl9jYnXnVV4Tz6FomUwse1awbiksTu2lxGV5uqlIV3jYLuZ1i+tRXgaMWEi0e9qrvCZ2k4t+r0TLnsXhxLiZ+L2gbGEwkhK8gOol0iCGHf9f3ZxD232NiedreJXX7xi8vn+mfSdFy86GVZ6p+USUHRymSYks5mNQenh9/9LyndS5+mpahzfSWm4yD4dpMoZ3JtlmWi+Dm0wvShgqj9+9LZ3PQkprL0xay03pxZ4RyrSwd3hJ3xmmpZchLiqP36S6TGsPQ1rLTeZhzwgRPrvD27bvmNSdd9J3hnH2fJne6wKoPX6TehHTmtSd1nJTejEYoUyzGyYAUsFF0g2a18VfVRAhG1jJ7ldFOd3HH2WbaWtI0xAsEsWNwQhlmt0wbVraIBVcmNCghe3VcZMNrGSDlzh6kVRu0/TGnvkZRAxGqESYEFzICtur4yZ77LLBi+rnkew83IuLw1dQW12ppGfK9MY+7l4404MxIoAJrESxUJkAmPVkQnfCqeoE1KzXXzEmJfRS6WECK5FGOh8pX2x7pnH3DKxbXI/a6kpcHL6ipNym94qZ/iRdoiQwGCEj7Tzci9lt+zG7bb/0RdqEBlhmJVMn1UMKpg9RAGODhdXNOUysqsD5octGl1sV+xy1vv5+LN9Z04MxIoDBCBmq88BxnB+6HKpBcjfAUYKTsJ/1awCKBQeq72LTelec1nKHoevZREQmYTBCRrK76sMkMbobsii9A6p7Foo1skF3sWECI5Pvik3owRIRdzl1PZuIyCRMYKXMshMXm3KT0d17TiqBMcpnVfBKunQmIq5bXG9cUqZsomhQYqVJSZcmlYUobZjASiXP7tWwg4nOA8elV2Dt7j0X2LMQ112zV4+Ms1dFRy6I6nyXoF4ik4ZpTCoLUVaxZ8RHqU8HzALnObQbStG7W9HzH9ddc7H965g6XOxY+TdDRG6i7TeDER/smjVXmEYvroYyCw2w6Hc9C8dKRMniME1E7JpNjuxQR5ghCpFnv4R5JksWGmfR77rJybAy0pI4S1RKGIz4yMqFV6Uk8yOCqAoU3fs14ZksqomcM5UPqYtSjqSk4bwRlRoGIyRM9UXcbqCacpO1PMTOa+VPVc9kWb+rB/Vb9mD9rp5IZYxq275j+YftidL54Lsoi92J4rNgiMzDnBESpjpXI2xejq7hEZn91m/ZgxELKC8DjrcvT6iEY81u24/zQ5dRW12JI60tQp9RVb/O7QCQSpIFkNp8LeabEX2GOSMUStBdXZgeiaA74rB3qLq62WX2u/y2aSgvG/1v0pzncNPSBkyvrcampQ3Cn1fV8+SsL9FtRlnszk1XDwXzzYjksWeECqThiakm9YzIliWJKbum3JkHHUsS59CUeiAqZZzaS6FkZYYIkHyDJ7IqarEGUkUDavI5tMt2cfgKzg9djjVQUBE8ElE0HKahUFY356RXKw2ioqvcaxsir8kkTYpO8XX/m+yqqMW68J3/rvohfSaw6whA7EMZXvWgcyYNE1uJ/DEYicjvApPmC4/KC7aKbXltQ+Q10bF7mSm+7n9zNngi+ysWKDj/PYtTUO062rS0QThgkgkORfevMggSLYN9Pltffz+V1wWiODEYicivwUhzQ6Lygq1iW17bEHlNJmlSdIpv0L+p7pHIYiKk6iRo2b+zOHqNRMuwbnE9ysuAEQupvC4QxYk5IxH5jUEnMTbN8W9/rJvs0J0IG6V8Ud5LlAVMYC0BaZotkPRFOE1145Slxsr0IIKI4scE1hKQpm78pIet0lQ3Tmke3nNTObwSlzTndpWiuM8Xvw/6MBhJMZNnTbglHRzorJtiF7SgmUBeS+On9QIZNvcmSaYERSQm7vPF74M+DEYoEWkKnKIqdkELmgnU3XvOqOmospyBU9A5N+X7UCwoSmsgmFVxB7GmBMmliMGIZrzYheOsN9PqUGYtEZHP6JyOKitNgRNQPCjidFyzxB3EmhIklyImsGqW1kRL3Zz1BkB5HWY9wVL2eydaH0nWW1Iz1lpffx8jVnof3EekExNYU4LdguE46y2OOlR9h5+G3pugMorWR5J3lkmco9XNObTd18i/UaKYhQpGOjo6MHPmTEyYMAFNTU04ePCg73vffvttLFy4ENdddx2qq6txyy234Nlnnw1d4KyJ++JtWiOoStz1pjrAMW34QnapdBOD5qTOEbvuieInHYzs3r0bGzZswNatW9HT04NFixZh2bJl6Ovr83z/xIkT8dhjj+Gtt97C0aNH8cQTT+CJJ57A9u3bIxeeijOtEYxDHMeoYzVV3YFjkqvLqhC1TEHPGRL9DBGpIZ0zMm/ePMydOxednZ3512bNmoUVK1agvb1daBt/+Id/iIkTJ+Lll18Wen+Wc0bilvXcB0DsGNNQD6WWP6T7nPjVd1C5Su0cEUUVS87IpUuX0N3djZaWloLXW1pacOjQIaFt9PT04NChQ7j77rtldk0hmXhHq5rIMaahh2jd4nrUVlfi4vCVkrjzdp4THT0Ofj0haRmuClNn7NkhU0kFI2fPnsXIyAjq6uoKXq+rq8Pp06cDP3vjjTeiqqoKd9xxB/78z/8cjz76qO97h4eHMTg4WPBD8VB1cTL9Ipd0I7LzcC9mt+3H7Lb9wk+bXd2cw8SqCpwfuqytgU6S85zIBIuq6sUviE3LcFWYADsNQTmVplAJrGVlZQW/W5Y15jW3gwcP4t1338U//MM/4LnnnsOuXbt839ve3o5Jkyblf2bMmBGmmCRA1cVJ90Vu/a4e1G/Zg/W7ejz/PelGpPPAcZwfupwPLNy27TuGk+eHsG3fsYLXwzbQYYk07HEFRc5zItMrFHe9mBRwBAkTYJvUs0PkJBWMTJkyBeXl5WN6Qc6cOTOmt8Rt5syZ+OIXv4g/+ZM/weOPP44nn3zS971btmzBwMBA/ufEiRMyxSQJqi5Oui9ye947hRFr9L+iRBvZMI2x3bjWVldK1Ym7gY67TkUa9jgaf3edunuFguj+rpkiTNCUlkCLSk+FzJvHjx+PpqYmdHV14f7778+/3tXVhfvuu094O5ZlYXh42Pffq6qqUFVVJVM0Cml1c07JhUl0O3ElLS6/bRr2vHcKy2+bJrwPZyOr4n1Oxepj09KGfBnDbiMKu46acpMBwLMcIu8Jy6tO7d6gYvuJs16ISA+pYAQANm7ciIceegh33HEH5s+fj+3bt6Ovrw9r164FMNqrcfLkSXz3u98FALzwwgv43Oc+h1tuuQXA6Loj27Ztw9e+9jWFh0Fp4dUIqQhQnn9wDp5/cE7BipnFgoem3GScHhjKN7Z+RBtJGbobVPs8APCdFSLynrC86lSmTnTPxEkz1h2ZSDoYWblyJT755BM89dRT6O/vR2NjI/bu3YtcbvRL3d/fX7DmyNWrV7Flyxb84he/QEVFBerr6/H000/jT//0T9UdBaWGVyMUpufBT+eB4xixgPKy4nfy3b3nMGKN/jeIXSZ7+CALF3CRACuOIMwWNRhT+Z0pNaw7MhGfTaOZzrsUU+6QVJZDZlsy7+X6EqN2Hu7NJ91uWtqg7XsT53fXlL+LuGT9+Mgsou03gxHNdDZybGDFrd/Vk89Jef7BOcKfc1/4424I4n6gnf2dAYIfHBd2+yY0lPy7IFKHD8pLCZ0zA6LsO+trYLiJDum4uWeiRJ2ZUmz9EtHthy2H6Cwhke17fYd0TxEHOFuHSAcGI4rJNtI6p9pF2bdfoyFz/GkKaMI2UE25ySgvQ8GMlCgNXbH1S0TX6yhWDr9zs7o5hyOtLTjS2hL4vRHZfuvr74/5DoV9PozK7xKnvxIlj8M0is1u24/zQ5dRW12JI60txT+QUn7d6TJd3PZ7a6srMbGqQqpr3oTufBGqu/xFcjZk9qniPIZh/52UAfj6isZQQ0XOsnFohchMHKYhaTJ3l2GW0naz3wsg0rLWpvWwOMsjWh+ix+DsmQDg+RmZc2DX47Z9x6SfYKvCpOpK6WDSq2xer5n2vSAifyUfjKi+YG1a2oDptdXYtLRB6X6TuLCqGK+X6eK232vXWdhlrVXlGaiqY/fUSZH60PWcEb+AMO6hCtG/Ey9eZfN6TXf+CYMhInElH4yovmDF1fgkcWHVlbgXdVlrVeVWVcdJPTMkzFNnvZZhDxsQRpFEXobuRNSwSbxEpajkc0Z05R7I7jctORJpKacX08ruXI69u/ec1HfFL6+kFPJ0TCFSX8x1oazjOiMJMPHi7FempMoaR8KmqUFb3Puy67K8DBixgtf18Pus+zN2mS8OX8H5octKtqmaCX9XSZXBhGMlihMTWBOge0zai1+Zkiqr6kRCE4ezktjXzsO9uDh8BbXVlVh+2zRlQziq8nTiZMLfVVJl4DRiolEMRiLQPSbtxa9M9utNucmxjlGrSCRcv6sH9Vv2YP2uHuk6TvKcqN6XM2iz1xOZWFWBO2deK70tr/Pg3L7Jj59P4hwWC5BN/NsmyjIO05QYHWPUsl3R9Vv25B92d7x9eeTtpaEr3Pm0YecsIedsoSjrhri3ryM/waTzYEKuhkn1QRQXDtOklMrseq9txXk370f2jnr5bdNQXjb6Xy8mD93IcPeEOJ82HHa2UNAwnejTjKMcR5iy6WBCz4dJ9UGkG4MRw6i8QHltS3VXexwX1OcfnIPj7ct9H0hn8tCNW1BD7aw7u4xt941djVTmnBUbpvPaflSi3wETAgCbCbkaJtUHkW4cpjGMyq7bJJ4QK/I4+VLujnZPpXVO0/0//+8E3js5gNumT8KPvva7Wsqn4tyU8vl1Y10QFeLU3hIjGhioJDrurmt83oSGwT4vA0OX4fxDm15bjdMDQ4G5MSLbjjqN24TciSxhfRIVYs5Iidm271jg01zj4H4irZ+g7ug4V6Dctu9Y/rkrSe7XaXVzDhOrKgoCETtvo1huTDEqpnHrnrmStRVIOfRivqx957KCwUgEJn6py6A+QdFPd+85jFij/w3iNz7v9xj5JCSZPGg3UH9w+7SCvI2g3BiR71ax/BD79aBtJZE7EVTXoufBtL81v/KYkItCwZg4bCYGIxGY9KW2F7GSfRx7FFEXOItzhgcQ/DC2YnewKhs/u4F6/sE5wg2VyHfLr+Fzv677expU16I9CbqPwc208pA49l6ZiTkjEZiQkxCnMMcnM2Zucv0lOfbvVQ9h6kb3owDiZNry7FmoU6IkMIE1YVm8OIVpkNNaD+5yJ3kcqgIf93bYsMpjAiqRWkxgTVgc3ba6x8nDdGfGMWaeRD24z5+q44iS+yG7Pfd2RL+TMt9d2XOh+zssi134RHowGFEkjotYXOPSog2ESIO883AvZrftx+y2/coaHHf5khifj6sRipL7Ibs993ZEj0nFCq+q3q8bE1CJ9GAwAjV3b3FcxHQ2kEG8Huimckqxu3xx1IP7nMfVCImW3Q7qZv3v/xsY2MnUhegxFXuoXtj9h3k/EZUm5owgfePEUR8UFzWXwFlf6xbXK19sLWwOQ5oX+7LLY9NdLtPqh4jSiTkjEnTcvUXpjYnaVS56x+y3H2d9rW7O4UhrC460tkRaTly2l8Kr/kxZ7CvMuV23uB611ZWorhyH2upK7T0JUadtUzx4Diir2DOiSZQ7z6g9I3HtJ+xnw9SF12eK7VPldFmRstnPoyn22fW7erDnvVNYfts034cDmkBVb4mzTgFwJo8E9lhR2rBnxHBR7sxlew6CVkANusuKkkcRdy+F12eKldevTKJP1pUtGwChz+557xRGrNH/ipZLB1W9Sc46TVuCq27MwaGsYs9IRoncQcV5l2Xi2hUiOTDuelDZO+T3+3UTx+OnpwbG9Iyk5S44Sk8dwJ4RoizjomclTqSBUB0wmBiAiEiq3O7goliwkdRwXFRJBU2qji+t31OiNOIwTYkTGWJRPZ1V5ZThODn3E0dAJrooWbEu92LnR8d6LF6SGjpQdXwcGiIyD4OREufVeIq+5ha1UYq7kbCPYdu+Y/n92P9vT0+OSmZRsqjBYNT1WFQFfyqC2iS+Xyq3Y1o+j5c0lJHIxmBEE1MuFF6Np+hrblEbJZkFwsLUnX0MAGK7k09igTb796bcZNRWV+Li8JVQicoqgj9V3+Ogstj7AFBwfGH3rSJ4SkPvShrKSGRjMKKJyIUiiYDFq/EUfS1ImLJHXf+k2P7tY9i0tKEgt6G2uhKbljYIlzNIHCu5uo/X/r279xwmVlXg/NBltL7+/pi69uoJclIR/Klq8ILK4rcPnY1tGma1pKGMRDYmsGoikquQhtkUMjNUkkxA1D2bSKassttwz0Zpff19jFhjV22VXe/ET1yzjfyIrhjMRFQi83E2TQZEvdgmMRtjdtt+nB+6jNrqShxpbcm/7lzI686Z16LzwHFcHL6C80OXAxv/JAOWJBrSuAOenYd7fZfjT+vskzQE4TZTAyJTy0Wlh7NpMuLi8BVs23dMOpkU0PuE1e7ecxixRv8rk6uhqgw6ZhMByTzkz72/80OXMbGqYsxxqDq+pJ9km6bhBVPzMkwtF5EfBiMGsxsa9xNxRS80qp+w6hUEbVrakM/D8NuWO1cjqFFLU0PkxV3+uBvytNeXl6SDHyfZXCdT69/UchH54TCNIby6Vdfv6sE//9cpTKgch63Lv1AwiyDLi1uJSHs3tGweRBy5JzRWkt9xng8qBRymSRmv3o7u3nOwAPzP5asFQzW67hzDzr5QPSto5+FetL7+fr6+TJkmbfMrj/N12RkiKrrd4+i6N63uo0qyR4FDKUSfYTBiCL/ptOVlgAWMGaoJElcDEXbqreppzJ0HjmPEAsrwWU6N7EVdZn+y9el1vO4Ayj7fTbnJBdv2awxVNJJxNLRZa1BlA/0of2scSiH6DIMRQ3hdBFc359B2XyNqqytRW105puHyI9pAqAxanNvyWvbcuThXlDI7tweMBmkAim4/yv5kG1yvRsYOoMrLkO+Wf2fzPejuPVewbb/GUEVvmOzy8iLS3qBG/RuIEozpzI0hMg2DkRiobORXN+dwpLUFR1pbxjRcfkQbCJV3tc5trW7OYd3i+vwQyurmXH5xLr8VNi8OX0FtdaVQo2Zvz8Jo475paUPg9qPuT7bB9Wpk7G0sv21avl78tq3q+6OiR6eYtDeoUf8G0h6MEZmCwUgM4uq6Fr3wuRsIv0bJb3uyjZhX4y4zxTVoeqofe3tt9zXmgx+/7buPR3Z/Mg2uX93J9ISo+v6o6NHJuqjHnPZgjMgUDEZiENfzScJm3vs9DM7vQmo3Ytv2HRMeFnI37jJTXMPUl9dD55y9MV7Ho2rtjyjLo4sMWanqLVHRo5NmInWWtWMmSisGIzEIc4ErduFMMlHQbsQAhB4WkqkD0Z6cYuw6cj+nRfXaH37nQmT4p9iQlV/5/JJi/erJpGmjumbcZC25lijLGIwYQvaOWuYC77cwmR+7MbQ/JzssFFXYRsSefTRioeCzqoId53686kV0+CdMz4xfUqxfPfkFZmFFqTNdQUEpDjsRpRUXPTOE7MPfABizABmg9k48zDN17CGou26+Ht295wo+G/V5MaLlsd/XlJuM7t5z+f+KHEeY+rOPe/jKVVRVjCt4No09ldjrAXphRFkMLOjYTOrBISL1+KA8Q8TVSAMw6iKuc3VWe9+Ad8Prfpif+5wUO0eyT7+132/30ojUSdj6Czr2OL57MgGWCJNW9SUi9WJdgbWjowMzZ87EhAkT0NTUhIMHD/q+97XXXsO9996L66+/HjU1NZg/fz727dsXZreppLKL2jncYFrinc4ucXsIS2ZqsF137sXI/LYfJodm+W3ThFeslZlq7N6X37Gr/I74zQaKikMpRASE6BnZvXs3HnroIXR0dGDhwoX41re+hRdffBEffPABPve5z415/4YNGzBt2jR86UtfQm1tLV566SVs27YN//Ef/4E5c+YI7VN3z0iUO8y4H1MPqOshkR2OMKVXppig8jp7MexpwsW249c7ELZeRHoHwm5b9blK27knIr1iG6aZN28e5s6di87Ozvxrs2bNwooVK9De3i60jVtvvRUrV67EX/3VXwm9X3cwYlpXsrMBvWZCJc4PXVaS/yB6nEHvc+ZvOHMYkiI7BBOmcfU7fvfrorkSQPGAMsx30CtvhMEEESUplmGaS5cuobu7Gy0tLQWvt7S04NChQ0LbuHr1Ki5cuIBrr73W9z3Dw8MYHBws+NFJ1eJgoopt1zlrBJBbCj1o2Ei0y1xkATN76urOw72Y3bYfs9v2Bz44ThXZ5+KEWdCsKTdZ6PkxIrNd7BVri5UhzHCGewl6535F15ApJmsPyiMiPaSCkbNnz2JkZAR1dXUFr9fV1eH06dNC23jmmWdw8eJFPPDAA77vaW9vx6RJk/I/M2bMkCmmcsUWB5MdP4+6poj9zBp7uq7MUuhBjZpow1xsATNnDoM7OJE5ThHuunQfn0gjLtqg2uXt7j0n9PyYoH0nsSCZe5Va52uAWP5LMUFrrjBIISJRFWE+VFZWVvC7ZVljXvOya9cuPPnkk3j99ddxww03+L5vy5Yt2LhxY/73wcFB7QGJF7uxlU2+c98Vh9muncTq3KZo8mac3fNe27eHbbx6EsLUn5O7LsMcX7HzEba8QWWJ+zz47cN+zT1MFJZfnYjWqYk4lEWkgSVheHjYKi8vt1577bWC19evX2/dddddgZ/9/ve/b1VXV1v/8i//IrNLy7Isa2BgwAJgDQwMSH9WlZf//WNrQfu/WS//+8eJbUvlPuOmo35E3rug/d+s3P/6F2tB+78p2V+c28iKl//9Y+v2J/dZtz+5L5X1IfKdISIxou231DDN+PHj0dTUhK6uroLXu7q6sGDBAt/P7dq1C4888gheeeUVLF++PEzMpJ27OzpKN7Rol3ualrN2ljXJx7IXq0uR4RAVU2DTdK7iJroSbVQqhoK8tsHpxkTJk15nZOPGjXjxxRexY8cOHD16FI8//jj6+vqwdu1aAKNDLGvWrMm/f9euXVizZg2eeeYZNDc34/Tp0zh9+jQGBgbUHUUCZJIT49qnyZxlLVY3fo2IM0G0troSv7447Jn4KkPlWhtBjV+UcxW2UU0iLyOJB/SFpeJv0Gsbpq3hQ1QKQq3A2tHRgW984xvo7+9HY2Mjnn32Wdx1110AgEceeQQff/wxDhw4AABYvHgxfvzjH4/ZxsMPP4zvfOc7QvvTPbXXi+4prCYTXdE0aGosgMBVVXXwK3exKbzFvieyK7wWK49Kpk1rd1KR28H8EKJ4cTn4BDgv1M5EPtUXtSQumGGeB6N6IbjPnrUyAqAMgIWqinJjgj2/Rc+CGuxiS9U7t3tx+IqyNWNUYWNNRFHEuhw8jZIZmogiar6KyPtly79t37H8ehWy7G5wAAXlWt2cw8SqCgxdvoqhyyO4dmIVjrS2FDz8zW/NEtv6XT2o37IH63f1SJdLtNzuJdGLTeEttlS97FOS3Z8LEwyKfn84ZEFESWAwEoHzQq16nNzZYETNVxF5v478FK9yBTXeQWuW2Pa8dwoj1uh/bSJBjAx3XQU12KubczjS2lIQVPnx2k4SC8MREekWap0RGkv1uhHOBsPdQMmudxFm3ZIgdsNYW12JTUsbhD4jWq6gcqxbXO+7ZoldrvEV4/A/l69i+W3T8q/bQYz9/1GHoZJYI8QWx3odKtZ3ISJSiTkjhtI5Vh82AVW3oARTmWRjr2e66MKcDSJKM+aMxCRqt7nI54vNzoh7OmexNUOiDunEdQx+5ZIZKgG8n+ni5C6/1/Go+p4ASDRnw6Rl3E0qCxHFi8EI5C56UcfbRT4v+oC1qPxyKYol5kZNnLSTX53bVNHw2OX6z1/8WjiJNSjYcj7TxUnkYXxB5ymOhGJV4tyv7DlmbgtR6WAwArmLXtReAZHPq3zAWhC/hNC4EnPteh6+MoLyMqApN3nMv/n1xsg0ZF5JrF7s4RjZYEvkYXzFnmxsYkJx3PuVDS7cZWFPCVF2MWcEpTEu73WMSSzc5n4gm996Gs732Y2W899lFgZbv6sHe947heW3TcPzD87xLZu9zfIy+PaCxEHV9y1t39uo5TU1V4mI/HHRs4TE2SAU23bQ4mHuxbm8Gvi4GzO/RNAwxxXHwmAyxx9UJl3BQBoaZ5V1pLu+iUgeE1gTkuRiZ37/3vr6+wVd1/bre947lf+8V/e73+ejcuaGeCWC2sM/9pCMm3tRtPW7evLHcNfN148Z4gkSVIfu4ZigYQDZvJAkpOHZRSrriAuwEWUXg5GIijUIorNnvN5TbNvrFtejvAwYsTBm4bDptdVYftu0/Oe9LuR+n4/KboCA0bVIrplQ6fsekUReZ1DV3XsOIxbQ3XtOqCwyDXZQmWTzQpKQhsZZdx0RUTpwmCZmIl3pUbrbo3Zdx9H1XSz/Q3S/XsNNALQMi6ka7gmzbx1MKw8RpRNzRjTwSgiVaXSLXfhNbCDC5rWE2ZZOfg9FBKBkgTTT8j9MKw8RpRNzRmLkN6ziNVVWpCtdtLtdd46Cl2JlkhlKSPr4ZKaK+q29UmyBNFGmDWeYVh4iyjY+myYEv+eFFHt2SlQmPFPEPXRiJ5KqKFPSx+d3HgHvXpqLw1ewbd8x3HXz9fny2tuJ2puT5PNuRJhWHiLKNg7ThJDElFjThyvsxNc0dOMXmwJtvx6U6zK7bX/+YXtRjtnkc0tEpBqHaVJM1XBFHCtWes3UUVUuv+Xpo/KrT/cUYucy9X7DFGUI3wvkt+IrEVGpYzASQty5DaoeROdsXKNMMXb+GzD64LbnH5yTzwVR9cwRv+XpoypWn86pyH5ToTctbcD02mp8fUX4lVqj5peYthy6aeUhovRiMCJp5+FeXBy+gtrqythyG1Q9iG74ytX8AmEy63rY73E2NkGfD/qcF7/gYN3ietRWV0rVrUiD6O4B8VvPZdPSBt969zonso1xsQfwFaMyCFYRSJiYUE1E6cRgRJJ99z6xqsK4MX+7cbCDEMDKLxAW5gF9zsam2MP7aqsrcXH4SkHg4reyq1+wtbo5hyOtLTjS2iJctzINYrHhmmLTqt2Nt+i+nb1KQSvPFhPHQwujBBKccUNEqjAYkaTrAuxsDP1yK+yyVVWMw4gFVFWUB67AWozzWP0+bwcfAPLDK14ru7obc1Vd/DLnI8q582q8RbfnngocNghQueKqiu9xGlaAJaJ0KOnZNKbMbBAph3MRKgAFOQ7umR1hjyvMQld+T9N1l8G9bZWLakU5jyoXnPObnRN1BVmvJx/HtWIuAwsiUokrsApQ1SAm8Wh05z7+8xe/xj//1ylMqByHrcu/EHujFDQFFhBrHIO2EbX8fqujimxXZVAUV8DlFYiqnFLN1VaJKC6c2itApqta9omuqsvh7BJ/68NfwcLoMIzKO1m/p9g6Z+UAYxcLC0oOdW7bmS8Rx5BDU26y9NRZlcNu7m2p2rZzO2G2GTahmIgoKSXdMyLKXh/Cb5GvpLu57QW4aqsrcaS1Jbb9iA7BON8b9EA8v4fmqS5veRmEZ6yEOXdpG9ZgzwcR6cKeEYWKrQ+RdCKfvebFpqUNAIIXEYuSJOo35dV5vPY+mnKT83fXzv3avSrb9h2L/Q48zNTZML1aJkxpDftcHYqGa6sQxYM9I78R9yPi4+R35yvSW6EyX8M5XGO/9uuLwxi6fBXVleNw9OvLjKvLtPaMsLdDD9Y7kRz2jEgKutvVOYVx/a4e1G/Zg/W7enzfE7SImNfrsouUBfHah/O1qopyAKP5LTsP9+J///D9fE+JCcKcWxOmtKahtyOLvQjF6j2Lx0yUBPaM/EaSd7sys1Hqt+zJDxEdb1+ufP9e025V8soZAYDa6kpsWtqgtYch6DyY0PuRdqXYi1CKx0wUhD0jkpK825VZBGv5bdNQXjb6Xz+yd2PunI84l7d37su53LsdiIjmXsRxxxl0Hpy/8243nDT03qhWisdMpAKDEQ1kpmo+/+AcHG9fjucfnOO7vSgJlSLL26tojO2ehk1LG/IzgGSCoGJLzAft16/sQefB+bu97237jikJSuwyrd/VE3uQozOQMmE4K2mleMxEKnCYJgPCDinYs12A0Rk6Iqu/hu16jrogWLHp1aL7leFeQfXi8BWcH7qsbBEzW5xd+hw2ICKdOEyjQFx3laq3G/ZuTPShfyq6nqMuCLa6OYe2+xqly6HieTTdvefwzuZ78lOqVSxiNvogQ//p4qpw2ICI0oA9IwHiuqs05W41y0maKo4tzvrJct2HxTohyh72jCgQ112lc7uqe0lEtud8pH3ax7f9jlfFwmTuRF+vJyWr2DaNMmExOSLSg8EI/Bu0uBoM53ZVXoDtvIpi29N50VcdfPkdi+pActu+Yzg/dBnnhy6zsYwJh5SISheDEehtnFU++KzYsvVR9qmKzMwUkcDFni58cfhKwfviCiTL4F23nP4bHXuLiEoXgxHE3zgHdfH7PSk3qFEr1htQ7NksXs+WCdOIynzW/QwbAMI9OEGBy+rmHCZWVSjpsQg6Hjt59esrGj3rjUMMREThMRiB/B2ZbANuz1qxG0z3552/izRqdtDRlJuc/1zY5L+oa5SIflZkZoq7XuzjBIIDF1XBpMwjAdzv5RADEVF4DEYcRJMUZRtw58qjzkW0vFb7FGnU7Iaxu/ec8EquQWUL24h6fdYvUHO/1ysA/Js9H+Dk+SH8zZ4PCt5TbEqtqu59mbpQFXwk1TOVBlk7HiISx6m9v+FcVAvwXojKvQhW2CmI7l6MYr+LbOc/f/Fr7HnvFJbfNi1wtVa/MqgSZdryzM17YH8Z7SXwRY5FRNzTRsMed7HPBZXblCniqmTteIiIU3ul2cmfZYDvEuXuoYawjVqxO/lt+44FPtnWa2pud+85jFhAd+85oTKonsVj9yjZOSFhegx+//Zp+cXARizgn//rlLI75bhyOty5MLLHXayHJajcWRsaytrxEJG4Ct0FMIU9fBJ05+x8jyrOHhm7wRkYuhz4GWcDZZdVtmxNuck4PTCEptzkaAeAz3JiAOQDtTCef3AO7px5Lf5mz1H8z+URTKgsH3OcYcnWj2hPivNJxGGOe3Vzrujqt37lLvbZtMna8RCROPaM/IZI3kGY3IRi4+Du6bidB47Dwujvm5Y2eG4rSu+DTbYnJYg7J6aYoDrpPHAcQ5dHMK22GluXz/I9zrBPKgYg9DnRnpS47+Y53ZWISgFzRmImmxMQNkcgzIPndC29HVRW0XLFlaMhWw4iIvIn2n4zGImZykYtaFumP4tFdj/F3hPlScUMMoiIksFgpMSFaXRNms1gUlmyhgEZESWFs2lSJOz6CsVyL2Rnj+iazeB1HDJrmMRZjiziarFEZJpQwUhHRwdmzpyJCRMmoKmpCQcPHvR9b39/P1atWoWGhgaMGzcOGzZsCFtWI8TRYIVtHFRP+9SVLOl1HF5libsRLZVGmlNoicg00sHI7t27sWHDBmzduhU9PT1YtGgRli1bhr6+Ps/3Dw8P4/rrr8fWrVtx++23Ry6wbqrX54hrjYo0zcIIOg5n8FesEY0aKCbVSOvugUnTd4OISoN0zsi8efMwd+5cdHZ25l+bNWsWVqxYgfb29sDPLl68GLNnz8Zzzz0nVUiTckZUjrfHlRex83BvfsG0TUsbUt3oeNWR3zlIMs8kyvfAr5wm5nKYWCYiSo9YckYuXbqE7u5utLS0FLze0tKCQ4cOhStpyqi8q4zrTnzbvmOBD+ZLE686KvbU4qhrnYiI0kPmV04Th4lMLBMRZY9UMHL27FmMjIygrq6u4PW6ujqcPn1aWaGGh4cxODhY8GOysA1b3N3lZfhsIbU0NChe9ehVR36NuUx9Rq2TKIGkXzlNzOUwsUxElD2hEljLysoKfrcsa8xrUbS3t2PSpEn5nxkzZijbtmr2cu4mNfb2U26/vqIRq5tzqWlQRAME2dVUvUStkzgCSRNzOUwsExFlj1QwMmXKFJSXl4/pBTlz5syY3pIotmzZgoGBgfzPiRMnlG1bNfdy7jJEelTc7xH5jLsBSUuDYgcITbnJQscsErz4fdarTtI8nEVElGZSwcj48ePR1NSErq6ugte7urqwYMECZYWqqqpCTU1NwY+p7Aa07b5G6cZepDF1vyfq8IKOBlc0oLIDhO7eczh5fgitr7+fT6AMmyMiU1/2e+39EhFRMqSHaTZu3IgXX3wRO3bswNGjR/H444+jr68Pa9euBTDaq7FmzZqCzxw5cgRHjhzBp59+il/96lc4cuQIPvjgAzVHoFmUXgeRxtT9nqjDC1GCmbCBjHOfIsNa6xbXo7wM+ScZR8kRkakv9351Y08NEZWKCtkPrFy5Ep988gmeeuop9Pf3o7GxEXv37kUuN9og9Pf3j1lzZM6cOfn/7+7uxiuvvIJcLoePP/44WukTFMcUR5FHprvfE+Ux6zsP9+Li8BXhp+u6OYMKmTLYSbT2f0d+M5m8KTfZ8/32tt31bQcIMvuWqS/3fnULW99ERGnDZ9MIimttiCTXcYi6Doeqh/G1vv4+RixIPTnXbpiTfFaN7jU2dO+fiCgqPptGMRVrQ+w83IvZbfsxu21/vus96POqu+nDDvHY5QAQORF2dXMObfc1CpVj275jOHl+CNv2HdMyI0jFlOgo5zCuxOOsDf9k7XiIShGDEUEq1oboPHC8YDGyYp9351pEveCGbdxUL4Ef5m5ftuwq6ktFABR3UmyY40zLujOisnY8RKWIwYgP0Yu8TCO5bnE9aqsrC3I2gj7vbAx1XnCj9qg461DmOOz1UjYtbZAus4r6UtEzEXdSrMlPZ06qxyIt6+gQkT/mjPhI8jknItKSP1AszyPMcST1mbjEWRaTjtPNtL8hIkqeaPvNYMSHaRd508pjc5fLboAA4Lbpk/DJxUv5O1bVD5aj+KhKVjbxO0tEyWEwkjGmNsjucu083Isnfvg+gNFVaY+3L/d8nww2askz9ftGROnC2TQZY+q4uLtcq5tz+IPbp6G8DFh+2zTf98lIy3L2WWLq942Isok9I0RERBQL9owQERFRKjAYISIiIq0YjBAREZFWDEaIiIhIKwYjREREpBWDESIiItKKwQgRERFpxWCEiIiItGIwQkRERFoxGCEiIiKtGIwQERGRVgxGiIiISCsGI0RERKRVhe4CiLAfLDw4OKi5JERERCTKbrftdtxPKoKRCxcuAABmzJihuSREREQk68KFC5g0aZLvv5dZxcIVA1y9ehWnTp3CNddcg7KyMt3FCTQ4OIgZM2bgxIkTqKmp0V2cksfzYR6eE/PwnJgnK+fEsixcuHAB06ZNw7hx/pkhqegZGTduHG688UbdxZBSU1OT6i9Q1vB8mIfnxDw8J+bJwjkJ6hGxMYGViIiItGIwQkRERFoxGFGsqqoKra2tqKqq0l0UAs+HiXhOzMNzYp5SOyepSGAlIiKi7GLPCBEREWnFYISIiIi0YjBCREREWjEYISIiIq0YjEjq6OjAzJkzMWHCBDQ1NeHgwYO+7+3v78eqVavQ0NCAcePGYcOGDckVtITInJPXXnsN9957L66//nrU1NRg/vz52LdvX4KlLQ0y5+Ttt9/GwoULcd1116G6uhq33HILnn322QRLWxpkzonTO++8g4qKCsyePTveApYgmXNy4MABlJWVjfn52c9+lmCJ48NgRMLu3buxYcMGbN26FT09PVi0aBGWLVuGvr4+z/cPDw/j+uuvx9atW3H77bcnXNrSIHtO3nrrLdx7773Yu3cvuru78aUvfQm///u/j56enoRLnl2y52TixIl47LHH8NZbb+Ho0aN44okn8MQTT2D79u0Jlzy7ZM+JbWBgAGvWrMGSJUsSKmnpCHtOjh07hv7+/vzPTTfdlFCJY2aRsDvvvNNau3ZtwWu33HKLtXnz5qKfvfvuu62/+Iu/iKlkpSvKObF94QtfsNra2lQXrWSpOCf333+/tXr1atVFK1lhz8nKlSutJ554wmptbbVuv/32GEtYemTPyZtvvmkBsM6dO5dA6ZLHnhFBly5dQnd3N1paWgpeb2lpwaFDhzSVqrSpOCdXr17FhQsXcO2118ZRxJKj4pz09PTg0KFDuPvuu+MoYskJe05eeuklHD9+HK2trXEXseRE+TuZM2cOpk6diiVLluDNN9+Ms5iJSsWD8kxw9uxZjIyMoK6uruD1uro6nD59WlOpSpuKc/LMM8/g4sWLeOCBB+IoYsmJck5uvPFG/OpXv8KVK1fw5JNP4tFHH42zqCUjzDn56KOPsHnzZhw8eBAVFWwmVAtzTqZOnYrt27ejqakJw8PDePnll7FkyRIcOHAAd911VxLFjhW/ZZLKysoKfrcsa8xrlKyw52TXrl148skn8frrr+OGG26Iq3glKcw5OXjwID799FMcPnwYmzdvxuc//3k8+OCDcRazpIiek5GREaxatQptbW24+eabkypeSZL5O2loaEBDQ0P+9/nz5+PEiRPYtm0bg5FSMmXKFJSXl4+JWs+cOTMmuqVkRDknu3fvxle/+lX84z/+I37v934vzmKWlCjnZObMmQCAL37xi/jlL3+JJ598ksGIArLn5MKFC3j33XfR09ODxx57DMDocKZlWaioqMD+/ftxzz33JFL2rFLVnjQ3N2Pnzp2qi6cFc0YEjR8/Hk1NTejq6ip4vaurCwsWLNBUqtIW9pzs2rULjzzyCF555RUsX7487mKWFFV/J5ZlYXh4WHXxSpLsOampqcFPfvITHDlyJP+zdu1aNDQ04MiRI5g3b15SRc8sVX8nPT09mDp1quri6aExeTZ1vv/971uVlZXWt7/9beuDDz6wNmzYYE2cONH6+OOPLcuyrM2bN1sPPfRQwWd6enqsnp4eq6mpyVq1apXV09Nj/fSnP9VR/EySPSevvPKKVVFRYb3wwgtWf39//uf8+fO6DiFzZM/J3//931s/+tGPrA8//ND68MMPrR07dlg1NTXW1q1bdR1C5oS5djlxNo16sufk2WeftX7wgx9YH374ofX+++9bmzdvtgBYr776qq5DUIrBiKQXXnjByuVy1vjx4625c+daP/7xj/P/9vDDD1t33313wfsBjPnJ5XLJFjrjZM7J3Xff7XlOHn744eQLnmEy5+T555+3br31Vuu3fuu3rJqaGmvOnDlWR0eHNTIyoqHk2SV77XJiMBIPmXPyt3/7t1Z9fb01YcIEa/Lkydbv/u7vWnv27NFQ6niUWZZlaeqUISIiImLOCBEREenFYISIiIi0YjBCREREWjEYISIiIq0YjBAREZFWDEaIiIhIKwYjREREpBWDESIiItKKwQgRERFpxWCEiIiItGIwQkRERFoxGCEiIiKt/j9gaoJKEgTFfAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(stellar_mass, stellar_radius, s=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "576a7f51",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_logg(stellar_mass, stellar_radius):\n",
    "    G = 6.67430e-8  # Gravitational constant in cm^3 g^-1 s^-2\n",
    "\n",
    "    # Convert input values to SI units\n",
    "    mass_g = stellar_mass * 1.989e30  * 1.e3 # Convert solar masses to g\n",
    "    radius_cm = stellar_radius * 6.955e8  * 1.e2 # Convert solar radii to cm\n",
    "\n",
    "    # Calculate surface gravity (log(g))\n",
    "    surface_gravity = np.log10((G * mass_g) / (radius_cm ** 2))\n",
    "    return surface_gravity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "7bb40265",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.023365205848966\n"
     ]
    }
   ],
   "source": [
    "# log(g) from the sampled mass and radius\n",
    "print(np.median(get_logg(stellar_mass, stellar_radius)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "697c1229",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.137416242462084\n"
     ]
    }
   ],
   "source": [
    "# the true log(g)\n",
    "print(get_logg(true_mass[0], true_radius[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "f3525408",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.046208739553704\n"
     ]
    }
   ],
   "source": [
    "# log(g) if you compute using the median of the samples of mass and radius\n",
    "print(get_logg(np.median(stellar_mass), np.median(stellar_radius)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61ca6eb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the log(g) we get from stellar mass and radius != to the true value and != to the value\n",
    "# from the distribution medians"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}