{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "32bbcf29",
   "metadata": {},
   "source": [
    "# Unit Conversions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8388d068",
   "metadata": {},
   "source": [
    "Pylanetary implements various solar-system-specific unit conversions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b0055d49",
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports\n",
    "from pylanetary.utils import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e58e99c7",
   "metadata": {},
   "source": [
    "## conversions for radio data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "06662d80",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result is 1.74e-01 Jy/beam\n",
      "which translates to 4.09e+10 Jy/sr.\n",
      "The above is using the Rayleigh-Jeans approximation, and could have been done in one step: 4.08e+10 Jy/sr\n",
      "Using the full Planck function we get: 3.97e+10 Jy/sr\n",
      "The angular area of the beam is 0.18 arcsec^2 assuming Gaussian shape\n"
     ]
    }
   ],
   "source": [
    "Tb = 100 #K\n",
    "freq = 115.271 #GHz\n",
    "beamx = 0.4 #arcsec\n",
    "beamy = beamx\n",
    "\n",
    "jybm = tb_to_jybm(Tb, freq, beamx, beamy)\n",
    "print(f'Result is {jybm:.2e} Jy/beam')\n",
    "\n",
    "jysr = jybm_to_jysr(jybm, beamx, beamy)\n",
    "print(f'which translates to {jysr:.2e} Jy/sr.')\n",
    "\n",
    "jysr_rj = rayleigh_jeans(Tb, freq)\n",
    "print(f'The above is using the Rayleigh-Jeans approximation, and could have been done in one step: {jysr_rj:.2e} Jy/sr')\n",
    "\n",
    "jysr = planck(Tb, freq)\n",
    "print(f'Using the full Planck function we get: {jysr:.2e} Jy/sr')\n",
    "\n",
    "beamA = beam_area(beamx, beamy)\n",
    "print(f'The angular area of the beam is {beamA:.2f} arcsec^2 assuming Gaussian shape')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef08acf5",
   "metadata": {},
   "source": [
    "## I/F calculation\n",
    "\n",
    "definition from Hammel et al 1989, DOI:10.1016/0019-1035(89)90149-8\n",
    "\n",
    "Let's calculate the I/F of Proteus in H-band from a measurement in erg s-1 cm-2 um-1.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "960b4d4a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Proteus solid angle is 6.55e-15 sr = 0.000278 arcsec^2\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG0CAYAAADO5AZFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA++0lEQVR4nO3de3hU1fn28XtmEpJwCIJICBANSgWRoyAxogVsJFXkYFulgEIBaVEoQgSVqiCoTbXloP2hFISqLQii1VJBKqaAICjlEF+xAiIiiCSAHAJBA5nZ7x+wNxkJkAnJ7FnD93NduS6zs2dmbUkmd579rLU8lmVZAgAAiBJetwcAAABQkQg3AAAgqhBuAABAVCHcAACAqEK4AQAAUYVwAwAAogrhBgAARJUYtwcQboFAQN98841q1Kghj8fj9nAAAEAZWJalw4cPq379+vJ6z16bueDCzTfffKOUlBS3hwEAAMph586datiw4VnPueDCTY0aNSSd+J+TmJjo8mgAAEBZFBQUKCUlxfk9fjYXXLixb0UlJiYSbgAAMExZWkpoKAYAAFGFcAMAAKIK4QYAAEQVwg0AAIgqhBsAABBVCDcAACCqEG4AAEBUIdwAAICoQrgBAABRhXADAACiiqvh5v3331e3bt1Uv359eTwevfXWW+d8zLJly3TNNdcoLi5OjRs31ksvvVTp4wQAAOZwNdwUFhaqVatWmjp1apnO//LLL9W1a1d17txZubm5GjFihO655x79+9//ruSRAgAAU7i6ceYtt9yiW265pcznT5s2TY0aNdLEiRMlSVdddZVWrlypyZMnKzMzs7KGCVSa4/6A8gu+d3sYwBnFx/pUp3qc28MAQmLUruCrV69WRkZG0LHMzEyNGDHijI8pKipSUVGR83lBQUFlDQ8IiWVZ6vbnldqUd9jtoQBn9czPW+rOa1PcHgZQZkaFm7y8PCUlJQUdS0pKUkFBgb777jslJCSc9pjs7GyNHz8+XEMEyqyoOOAEmyoxXnlcHg/wQ8UBS/6ApY+/Pki4gVGMCjflMWbMGGVlZTmfFxQUKCWFH1K4zx+wnP/+f+O6KD7W5+JogNP9OedzTVyyRQHLOvfJQAQxKtzUq1dP+fn5Qcfy8/OVmJhYatVGkuLi4hQXx/1iRJ7iEuHG56Vug8jj8534viz2E25gFqPWuUlPT1dOTk7QsSVLlig9Pd2lEQHlFygZbjyEG0Qe+/vST+UGhnE13Bw5ckS5ubnKzc2VdGKqd25urnbs2CHpxC2lfv36OecPGTJE27Zt04MPPqhNmzbp+eef12uvvaaRI0e6MXzgvJSs3Hip3CAC2RXFkrdQARO4Gm7Wrl2rNm3aqE2bNpKkrKwstWnTRmPHjpUk7d692wk6ktSoUSMtXLhQS5YsUatWrTRx4kS9+OKLTAOHkew+hhiCDSIU4QamcrXnplOnTrLOUu4sbfXhTp06acOGDZU4KiA87MoNVRtEKsINTGVUzw0QTeyeGyo3iFSEG5iKcAO4xK7c0EyMSOU0FBNuYBjCDeAS+xeGPd0WiDRO5YbZUjAM4QZwiZ/KDSIct6VgKsIN4BIn3NBzgwhFuIGpCDeAS+yp4IQbRCrCDUxFuAFcUkzlBhEuhnADQxFuAJdwWwqRzsv2CzAU4QZwCeEGkS7GR+UGZiLcAC5hthQinZd1bmAowg3gEio3iHQx3hO/Igg3MA3hBnCJn9lSiHAnsw3hBsYh3AAu8QcCkthbCpGLyg1MRbgBXOI/kW3YFRwRy2dXbpgtBcMQbgCX2JUbGooRqeyG4mI/4QZmIdwALrErN/TcIFLZt6UCVG5gGMIN4JJiu3JDuEGEshuKi+m5gWEIN4BL2FsKkc6p3BBuYBjCDeASu4+BcINI5aNyA0MRbgCX2JUbpoIjUvmo3MBQhBvAJfZfw15mSyFC+dg4E4Yi3AAusf8atjcnBCKN7+T3JrelYBrCDeASKjeIdHblhttSMA3hBnCJvaQ9PTeIVHazO5UbmIZwA7jEbihm+wVEqpIz+ajewCSEG8AlxVRuEOFKhhuqNzAJ4QZwif2XMOvcIFIFVW6YMQWDEG4AlxQTbhDhYqjcwFCEG8AlTuWG2VKIUCVn8vkJNzAI4QZwyanKDT+GiEwlKzeEG5iEd1XAJX5n40yXBwKcgZdwA0Pxtgq4xO9nKjgin90TRriBSQg3gEv8bJwJAzjhhtlSMAjhBnCJn4ZiGMDZPNNPuIE5CDeAS/w0FMMAMVRuYCDeVQGXnAo3Lg8EOAuv03MTcHkkQNnxtgq4hMoNTOBUbsg2MAjvqoBLqNzABF5mS8FAvK0CLjm1zg0/hohcMYQbGIh3VcAlzgrFTJZCBLO3YKChGCYh3AAucfaW4r4UIliMj4ZimId3VcAlrHMDEzjr3JBtYBDCDeASO9ywQjEimb1CcTGVGxiEcAO4xO5hYG8pRDI73JBtYBLCDeASKjcwAZUbmIhwA7jEDjdUbhDJnMoNs6VgEMIN4JJiGophAHsqeDEbZ8IghBvAJc5UcCo3iGAxVG5gIMIN4JJiwg0M4HV6bgg3MAfhBnCJ/ZcwDcWIZGy/ABMRbgCX2D0MNBQjkvkINzAQ4QZwCZUbmIBwAxMRbgCX2D0MXmZLIYLZs/loKIZJCDeAS+zZUjFsC44I5qOhGAYi3AAuoXIDE5zafoFwA3MQbgCXsP0CTEDlBiYi3AAu8bPODQxAQzFM5Hq4mTp1qlJTUxUfH6+0tDStWbPmrOdPmTJFTZo0UUJCglJSUjRy5Eh9//33YRotUHHsXcEJN4hkhBuYyNVwM2/ePGVlZWncuHFav369WrVqpczMTO3Zs6fU8+fMmaOHH35Y48aN02effaaZM2dq3rx5+t3vfhfmkQPnj8oNTGDPlvIzWwoGcTXcTJo0SYMHD9aAAQPUrFkzTZs2TVWrVtWsWbNKPX/VqlXq0KGD+vTpo9TUVHXp0kW9e/c+a7WnqKhIBQUFQR9AJCDcwAT2bD4/G2fCIK6Fm2PHjmndunXKyMg4NRivVxkZGVq9enWpj7n++uu1bt06J8xs27ZNixYt0q233nrG18nOzlbNmjWdj5SUlIq9EKCcAuwKDgN4qdzAQDFuvfC+ffvk9/uVlJQUdDwpKUmbNm0q9TF9+vTRvn37dMMNN8iyLBUXF2vIkCFnvS01ZswYZWVlOZ8XFBQQcBAR2DgTJmBvKZjI9YbiUCxbtky///3v9fzzz2v9+vX6xz/+oYULF+qJJ54442Pi4uKUmJgY9AFEAhqKYQIv4QYGcq1yU6dOHfl8PuXn5wcdz8/PV7169Up9zGOPPaa7775b99xzjySpRYsWKiws1K9//Ws98sgj8nqNymq4wNFzAxM4DcWEGxjEtTRQpUoVtW3bVjk5Oc6xQCCgnJwcpaenl/qYo0ePnhZgfD6fJMnifjAMYlkW4QZG8PkINzCPa5UbScrKylL//v3Vrl07tW/fXlOmTFFhYaEGDBggSerXr58aNGig7OxsSVK3bt00adIktWnTRmlpadq6dasee+wxdevWzQk5gAlK/p6goRiRzP7+ZIVimMTVcNOrVy/t3btXY8eOVV5enlq3bq3Fixc7TcY7duwIqtQ8+uij8ng8evTRR7Vr1y5dcskl6tatm5566im3LgEol5J/BfvYOBMRzG4oZldwmMTVcCNJw4YN07Bhw0r92rJly4I+j4mJ0bhx4zRu3LgwjAyoPEHhhsoNIpiXvaVgIDpwAReUXDOEnhtEshh2BYeBCDeAC0qu9kq4QSRjKjhMRLgBXBBUueG2FCIYi/jBRIQbwAXFgYAkyeM59ZcxEInYfgEmItwALjiZbZy/ioFIFUNDMQxEuAFcYFduvNySQoTz0VAMAxFuABdQuYEpfCfXGqNyA5MQbgAXOJUbwg0inO/kbwkqNzAJ4QZwgb3aK5UbRDoqNzAR4QZwgf/kbSnWuEGkcyo3zJaCQQg3gAtoKIYp7O/RYj/hBuYg3AAuoKEYpog5eVuKdW5gEsIN4AIaimEK+7YUKxTDJIQbwAU0FMMUdkMx4QYmIdwALrD7F6jcINJRuYGJCDeAC/xUbmAIKjcwEeEGcIH9i4LZUoh09q71TAWHSQg3gAvscBPjI9wgsvnYOBMGItwALrDDjY/KDSIcG2fCRIQbwAVOuKHnBhGOyg1MRLgBXEC4gSns71EaimESwg3gAnu2FOEGkS6GcAMDEW4AF1C5gSnsGX1svwCTEG4AF5wKN/wIIrLZM/qo3MAkvLMCLih2Zku5PBDgHJzKDeEGBiHcAC4IULmBIei5gYl4ZwVccKqh2OWBAOfAbCmYiLdWwAU0FMMUXsINDES4AVxAQzFM4dyWYrYUDMI7K+ACPw3FMETJhmKLgANDEG4AF1C5gSliStw65c4UTME7K+ACZyo4P4GIcN4S4aY4EHBxJEDZ8dYKuICp4DBFUOWGbAND8M4KuIDKDUxRckYfTcUwBW+tgAsCJ39JxFC5QYQLCjd+wg3MwDsr4AK7cmPPRAEilc9D5QbmIdwALrB7bmKYC44I5/V6ZOcbGophCsIN4AIqNzCJXb0h28AUhBvABfY6NzFsvwAD2H03VG5gCsIN4AI73HgJNzCAHW7INjAF4QZwgd+icgNzULmBaQg3gAvsKbXsCg4TOJUbZkvBEIQbwAV25YaGYpjAbiguZnMpGIJwA7iAhmKYxK7c+Ak3MAThBnABDcUwCeEGpiHcAC6goRgmIdzANIQbwAV2QzGVG5iAcAPTEG4AF1C5gUkINzAN4QZwgf1LwsdsKRjA/j4l3MAUhBvABU64oXIDAziVG9a5gSEIN4ALCDcwCbelYJqY8jzo4MGDWrNmjfbs2aPAD5bj7tevX4UMDIhmhBuYJIZwA8OEHG7+9a9/qW/fvjpy5IgSExPlKdEz4PF4CDdAGRBuYBIv4QaGCfm21AMPPKCBAwfqyJEjOnjwoA4cOOB87N+/vzLGCEQdu3eBcAMTULmBaUION7t27dLw4cNVtWrVyhgPcEEoZrYUDGLvgUZDMUwRcrjJzMzU2rVrK2MswAUjYIcbH+EGkS/GR+UGZgk53HTt2lWjR4/W448/rjfeeEMLFiwI+gjV1KlTlZqaqvj4eKWlpWnNmjVnPf/gwYMaOnSokpOTFRcXpyuvvFKLFi0K+XUBN1G5gUm8rHMDw4TcUDx48GBJ0oQJE077msfjkd/vL/NzzZs3T1lZWZo2bZrS0tI0ZcoUZWZmavPmzapbt+5p5x87dkw333yz6tatq9dff10NGjTQV199pYsuuijUywBcFWBXcBjE/j4tJtzAECGHmx9O/T4fkyZN0uDBgzVgwABJ0rRp07Rw4ULNmjVLDz/88Gnnz5o1S/v379eqVasUGxsrSUpNTT3raxQVFamoqMj5vKCgoMLGD5RX8cmfI/aWggnsxvcA4QaGcG0Rv2PHjmndunXKyMg4NRivVxkZGVq9enWpj1mwYIHS09M1dOhQJSUlqXnz5vr9739/1mpRdna2atas6XykpKRU+LUAobJ/RzBbCiawb0tRuYEpyhVuli9frm7duqlx48Zq3LixunfvrhUrVoT0HPv27ZPf71dSUlLQ8aSkJOXl5ZX6mG3btun111+X3+/XokWL9Nhjj2nixIl68sknz/g6Y8aM0aFDh5yPnTt3hjROoDLYlRvCDUxgNxQHmC0FQ4Qcbv7+978rIyNDVatW1fDhwzV8+HAlJCToJz/5iebMmVMZY3QEAgHVrVtX06dPV9u2bdWrVy898sgjmjZt2hkfExcXp8TExKAPwG323V0aimECp3LjJ9zADCH33Dz11FN65plnNHLkSOfY8OHDNWnSJD3xxBPq06dPmZ6nTp068vl8ys/PDzqen5+vevXqlfqY5ORkxcbGyufzOceuuuoq5eXl6dixY6pSpUqolwO4ghWKYRK7oZjKDUwRcuVm27Zt6tat22nHu3fvri+//LLMz1OlShW1bdtWOTk5zrFAIKCcnBylp6eX+pgOHTpo69atQU3NW7ZsUXJyMsEGRikm3MAgXmZLwTAhh5uUlJSgQGJ77733Qm7WzcrK0owZM/Tyyy/rs88+07333qvCwkJn9lS/fv00ZswY5/x7771X+/fv1/33368tW7Zo4cKF+v3vf6+hQ4eGehmAq+y/gJkKDhOw/QJME/JtqQceeEDDhw9Xbm6urr/+eknSBx98oJdeeknPPvtsSM/Vq1cv7d27V2PHjlVeXp5at26txYsXO03GO3bskNd7Kn+lpKTo3//+t0aOHKmWLVuqQYMGuv/++/XQQw+FehmAq4r9TAWHOZgKDtOEHG7uvfde1atXTxMnTtRrr70m6UTfy7x589SjR4+QBzBs2DANGzas1K8tW7bstGPp6en68MMPQ34dIJLYvyOo3MAEPm5LwTAhhxtJuv3223X77bdX9FiAC4aziB+zpWAAe1YfDcUwhWuL+AEXMrsnPoaNM2EA38n2ACo3MEWZKje1a9fWli1bVKdOHdWqVUues/y1uX///gobHBCtnEX8qNzAAL6TfwbTcwNTlCncTJ48WTVq1HD++2zhBsDZWZbF9gswCpUbmKZM4aZ///7Of//qV7+qrLEAF4SS02kJNzCBXblhKjhMEXLPzfr16/XJJ584n//zn/9Uz5499bvf/U7Hjh2r0MEB0chvEW5gFrtyQ7iBKUION7/5zW+0ZcsWSSdWK+7Vq5eqVq2q+fPn68EHH6zwAQLRhsoNTGP3hvmZLQVDhBxutmzZotatW0uS5s+fr44dO2rOnDl66aWX9MYbb1T0+ICoUzLcMBUcJnBuS7FxJgwRcrixLMvZ2+m9997TrbfeKunE6sH79u2r2NEBUahkuGERP5jAuS1F5QaGCDnctGvXTk8++aT+9re/afny5eratask6csvv3S2TQBwZtyWgmloKIZpQg43U6ZM0fr16zVs2DA98sgjaty4sSTp9ddfd/aaAnBm9i8Ir0csqwAj0FAM04S8/ULLli2DZkvZ/vjHP8rn81XIoIBoZpf2qdrAFPZC2oQbmCLkys3OnTv19ddfO5+vWbNGI0aM0CuvvKLY2NgKHRwQjYr9hBuYxeejcgOzhBxu+vTpo6VLl0qS8vLydPPNN2vNmjV65JFHNGHChAofIBBt7M0H2XoBprC/V1mhGKYIOdxs3LhR7du3lyS99tprat68uVatWqXZs2frpZdequjxAVHH/uuXyg1MYc/qY1dwmCLkcHP8+HHFxcVJOjEVvHv37pKkpk2bavfu3RU7OiAKEW5gGu/J71VuS8EUIYebq6++WtOmTdOKFSu0ZMkS/fSnP5UkffPNN7r44osrfIBAtDnVUBzyjx/gihjCDQwT8rvr008/rb/85S/q1KmTevfurVatWkmSFixY4NyuAnBmpxqKXR4IUEZUbmCakKeCd+rUSfv27VNBQYFq1arlHP/1r3+tqlWrVujggGhk9y3EULmBIajcwDQhhxtJ8vl8QcFGklJTUytiPEDUs2eckG1gCi8bZ8IwZQo311xzjXJyclSrVi21adPmrKuqrl+/vsIGB0SjQIDKDcxiV26YCg5TlCnc9OjRw5kh1bNnz8ocDxD1iktsvwCYwJ7ZFyDcwBBlCjfjxo0r9b8BhI7KDUzjo3IDw5Sr58Z25MgRBQKBoGOJiYnnNSAg2p3quaF0AzNQuYFpQv7T8csvv1TXrl1VrVo11axZU7Vq1VKtWrV00UUXndZkDOB0p9a5cXkgQBl5ne0XAuc4E4gMIVdu7rrrLlmWpVmzZikpKemszcUATuf3s4gfzBLjs7dfcHkgQBmFHG4+/vhjrVu3Tk2aNKmM8QBRz6nc8HcBDEHlBqYJ+U/Ha6+9Vjt37qyMsQAXBD8NxTCMs3Em2QaGCLly8+KLL2rIkCHatWuXmjdvrtjY2KCvt2zZssIGB0QjP4v4wTCnZkuRbmCGkMPN3r179cUXX2jAgAHOMY/HI8uy5PF45Pf7K3SAQLShcgPT+JztF1weCFBGIYebgQMHqk2bNnr11VdpKAbKwc9UcBjGmQrO9gswRMjh5quvvtKCBQvUuHHjyhgPEPX8zsaZhBuYwbktRekGhgi5Ln7TTTfp448/royxABcEp3JD1ROG8HmYCg6zhFy56datm0aOHKlPPvlELVq0OK2huHv37hU2OCAaneq5IdzADDQUwzQhh5shQ4ZIkiZMmHDa12goBs7NDjc+wg0M4WMqOAwTcrj54V5SAEJDuIFpYqjcwDAVMhf14MGDFfE0wAWBcAPTeL2nem4sZkzBACGHm6efflrz5s1zPr/jjjtUu3ZtNWjQgEZjoAxObZxJuIEZSvaH+ekqhgFCDjfTpk1TSkqKJGnJkiV67733tHjxYt1yyy0aPXp0hQ8QiDZO5YbZUjBEyTWZ/FRuYICQe27y8vKccPP222/rzjvvVJcuXZSamqq0tLQKHyAQbZxww86ZMASVG5gm5MpNrVq1nI0zFy9erIyMDEkn7sMyUwo4t2IqNzBMyTWZCDcwQciVm5/97Gfq06ePfvSjH+nbb7/VLbfcIknasGEDqxYDZRCgoRiG8VG5gWFCDjeTJ09Wamqqdu7cqWeeeUbVq1eXJO3evVv33XdfhQ8QiDbFhBsYxkflBoYJOdzExsZq1KhRpx0fOXJkhQwIiHYBZkvBMF6vRx6PZFmEG5gh5HAjSZ9//rmWLl2qPXv2nLao39ixYytkYEC0KvYTbmCeGK9Hx/0Ws6VghJDDzYwZM3TvvfeqTp06qlevnjwlypUej4dwA5yDU7mhoRgGOdFUbDnhHIhkIYebJ598Uk899ZQeeuihyhgPEPXsJeyp3MAkMV6PinQqnAORLOSp4AcOHNAdd9xRGWMBLgj+k3dyCTcwib2QHz03MEHI4eaOO+7Qu+++WxljAS4ITAWHiWIINzBIyLelGjdurMcee0wffvihWrRoodjY2KCvDx8+vMIGB0QjpoLDRPb3Kw3FMEHI4Wb69OmqXr26li9fruXLlwd9zePxEG6Ac7B7FmIINzCIHW5oKIYJQg43X375ZWWMA7hg2JUbL7OlYBB7dh8NxTBByD03AM6P3XMTw8aZMIi90WsxPTcwQLkW8fv666+1YMEC7dixQ8eOHQv62qRJkypkYEC0sqeCU7mBSZzKDeEGBgg53OTk5Kh79+66/PLLtWnTJjVv3lzbt2+XZVm65pprKmOMQFSxp4LTcwOTOD03hBsYIOTbUmPGjNGoUaP0ySefKD4+Xm+88YZ27typjh07sv4NUAZ+u3JDuIFB7HBD5QYmCDncfPbZZ+rXr58kKSYmRt99952qV6+uCRMm6Omnny7XIKZOnarU1FTFx8crLS1Na9asKdPj5s6dK4/Ho549e5brdQE32JNN2H4BJrFvo1K5gQlCDjfVqlVz+mySk5P1xRdfOF/bt29fyAOYN2+esrKyNG7cOK1fv16tWrVSZmam9uzZc9bHbd++XaNGjdKNN94Y8msCbrIrNzQUwyT29yvr3MAEIYeb6667TitXrpQk3XrrrXrggQf01FNPaeDAgbruuutCHsCkSZM0ePBgDRgwQM2aNdO0adNUtWpVzZo164yP8fv96tu3r8aPH6/LL7885NcE3ORnKjgMZFca/axzAwOEHG4mTZqktLQ0SdL48eP1k5/8RPPmzVNqaqpmzpwZ0nMdO3ZM69atU0ZGxqkBeb3KyMjQ6tWrz/i4CRMmqG7duho0aNA5X6OoqEgFBQVBH4Cb7HBDQzFMwgrFMElIs6X8fr++/vprtWzZUtKJW1TTpk0r94vv27dPfr9fSUlJQceTkpK0adOmUh+zcuVKzZw5U7m5uWV6jezsbI0fP77cYwQqmlO5IdzAID72loJBQqrc+Hw+denSRQcOHKis8ZzV4cOHdffdd2vGjBmqU6dOmR4zZswYHTp0yPnYuXNnJY8SODsqNzAR4QYmCXmdm+bNm2vbtm1q1KjReb94nTp15PP5lJ+fH3Q8Pz9f9erVO+38L774Qtu3b1e3bt2cYwG7OTMmRps3b9YVV1wR9Ji4uDjFxcWd91iBimKX9ancwCTOVHBuS8EAIffcPPnkkxo1apTefvtt7d69+7z6WapUqaK2bdsqJyfHORYIBJSTk6P09PTTzm/atKk++eQT5ebmOh/du3dX586dlZubq5SUlFAvBwg7e+NBKjcwic974tcFG2fCBGWu3EyYMEEPPPCAbr31VklS9+7d5Skx28OyLHk8Hvn9/pAGkJWVpf79+6tdu3Zq3769pkyZosLCQg0YMECS1K9fPzVo0EDZ2dmKj49X8+bNgx5/0UUXSdJpx4FIZf/lyzo3MIm9cgENxTBBmcPN+PHjNWTIEC1durRCB9CrVy/t3btXY8eOVV5enlq3bq3Fixc7TcY7duyQ18v+noge9iJoPio3MIhduaHnBiYoc7ixTqb1jh07Vvgghg0bpmHDhpX6tWXLlp31sS+99FKFjweoTAHCDQzkO/k3JuEGJgipJOKhjA6cN7usT7iBSWKo3MAgIc2WuvLKK88ZcPbv339eAwKinb3CK+EGJvEyFRwGCSncjB8/XjVr1qyssQAXBCo3MFEM4QYGCSnc/PKXv1TdunUrayzABcFPzw0MZO+FxmwpmKDMPTf02wAVgxWKYSIqNzBJmcONRVoHKkQxu4LDQPTcwCRlvi1lb3MA4PwwFRwmsqeCFxNuYABWxwPCjEX8YCJ7KniAcAMDEG6AMAswWwoGsm+jUrmBCQg3QJhRuYGJYnzsCg5zEG6AMAoELNm/G9g4EyZxKjfsCg4DEG6AMCq5RkgMG8LCIPZUcCo3MAHvrkAYlZxGS7aBSZgKDpPw9gqEUclfDFRuYBK7ckNDMUzAuysQRiVvS5FtYBK7AZ6p4DABb69AGPn9VG5gJh+VGxiEd1cgjIIqN0yWgkHs2X00FMMEhBsgjEpuvcBmtDAJlRuYhHADhJGzgB/BBoah5wYmIdwAYeRndWIY6lTlhk2UEfkIN0AYEW5gKp+zzo3LAwHKgHADhJHdUEy2gWnsW6l+KjcwAOEGCCO7chPj40cPZnEqN7TcwAC8wwJhZIcbLw3FMMyp21JUbhD5CDdAGDmVG+5LwTA+9paCQQg3QBjRUAxTEW5gEsINEEbFhBsYinADkxBugDCyl64n3MA0zmwpsg0MQLgBwqjYT7iBmXw+GophDsINEEZO5YbZUjDMqXVuXB4IUAaEGyCM6LmBqWKYCg6DEG6AMAoQbmAoLw3FMAjhBggjKjcwVQzhBgYh3ABhxDo3MJVTubEIN4h8hBsgjAg3MJVTuWEuOAxAuAHCyM9sKRjK66FyA3MQboAwCji7ghNuYJYYHz03MAfhBgijYnYFh6FOrXNDuEHkI9wAYcRUcJjKbiguJtzAAIQbIIyYCg5T2Q3FAcINDEC4AcKIhmKYyr6VSuUGJiDcAGHkP7kxj4+GYhjGbigOMFsKBiDcAGFkLxFC5Qam8VG5gUEIN0AY2ZsOxtBzA8PYfWKWJVlUbxDhCDdAGJ28K+XMPAFMUbIJnungiHSEGyCMqNzAVCXDDbemEOkIN0AYUbmBqUqGG5qKEekIN0AYUbmBqajcwCSEGyCM7HVu2H4Bpik5w4+F/BDpCDdAGNl/8VK5gWmo3MAkhBsgjNhbCqbyeDyyv22p3CDSEW6AMGJvKZgsxnviVwaVG0Q6wg0QRlRuYLKT2YZ1bhDxCDdAGNFQDJPZTcWEG0Q6wg0QRn4aimEwu+LoZ50bRDjCDRBGdrhhET+YyAk3VG4Q4SIi3EydOlWpqamKj49XWlqa1qxZc8ZzZ8yYoRtvvFG1atVSrVq1lJGRcdbzgUjCVHCYzHey6YZwg0jneriZN2+esrKyNG7cOK1fv16tWrVSZmam9uzZU+r5y5YtU+/evbV06VKtXr1aKSkp6tKli3bt2hXmkQOho6EYJvPRUAxDuB5uJk2apMGDB2vAgAFq1qyZpk2bpqpVq2rWrFmlnj979mzdd999at26tZo2baoXX3xRgUBAOTk5YR45EDqmgsNkMVRuYAhXw82xY8e0bt06ZWRkOMe8Xq8yMjK0evXqMj3H0aNHdfz4cdWuXbvUrxcVFamgoCDoA3CLveEg4QYmsqeCs84NIp2r4Wbfvn3y+/1KSkoKOp6UlKS8vLwyPcdDDz2k+vXrBwWkkrKzs1WzZk3nIyUl5bzHDZRXsZ9wA3PZlRt2BUekc/221Pn4wx/+oLlz5+rNN99UfHx8qeeMGTNGhw4dcj527twZ5lECpziVG9a5gYHsTM5tKUS6GDdfvE6dOvL5fMrPzw86np+fr3r16p31sX/605/0hz/8Qe+9955atmx5xvPi4uIUFxdXIeMFzhc9NzAZPTcwhauVmypVqqht27ZBzcB2c3B6evoZH/fMM8/oiSee0OLFi9WuXbtwDBWoEH7CDQzmZZ0bGMLVyo0kZWVlqX///mrXrp3at2+vKVOmqLCwUAMGDJAk9evXTw0aNFB2drYk6emnn9bYsWM1Z84cpaamOr051atXV/Xq1V27DqAsCDcwWQzhBoZwPdz06tVLe/fu1dixY5WXl6fWrVtr8eLFTpPxjh075PWeKjC98MILOnbsmH7xi18EPc+4ceP0+OOPh3PoQMgINzAZlRuYwvVwI0nDhg3TsGHDSv3asmXLgj7fvn175Q8IqCTsLQWT2d+3TAVHpDN6thRgGnYFh8nsWX5MBUekI9wAYeRUbnyEG5jHR+UGhiDcAGHk7ApO5QYGssNNgHCDCEe4AcKIhmKYzEvlBoYg3ABhRLiByWKo3MAQhBsgjPxsvwCD2bdTqdwg0hFugDCioRgmcxbxY7YUIhzhBggjGophMvt2qt8fcHkkwNkRboAwOrWIHz96MI8TbijcIMLxDguEkVO54ScPBmIqOEzBWywQRlRuYDIW8YMpeIcFwsiZLcVPHgzE9gswBW+xQBj5/Xa44UcP5vGdnOVXTNMNIhzvsEAYsc4NTGZ/3zIVHJGOcAOEkd2r4GOdGxjImS0VYCo4IhvhBggje5YJlRuY6FS4cXkgwDkQboAwKmZvKRgshsoNDEG4AcKk5NoghBuYyEvlBoYg3ABhUnJtEG5LwUROQzGVG0Q4wg0QJiXXBqGhGCbysXEmDEG4AcKEyg1Md6qhmHCDyEa4AcLET88NDEe4gSkIN0CY0FAM07G3FExBuAHCpOQvBLINTBTDruAwBOEGCJOAdWqNGw89NzCQ10PlBmYg3ABhwgJ+MF2Mj13BYQbCDRAmbL0A03k9NBTDDIQbIEzsyk0MlRsYKobZUjAE4QYIE/sXgpdwA0N5CTcwBOEGCBM/lRsYLoap4DAE4QYIEyo3MJ3dDE9DMSId4QYIEyo3MJ2ziJ+fcIPIRrgBwsTebNDLbCkYyp7pR+UGkY5wA4SJPxCQdGqtEMA0bL8AUxBugDDxn8g2rHMDY/nYfgGGINwAYVJ8snJDQzFM5aVyA0MQboAwOZltaCiGsVjED6Yg3ABh4lRuuC0FQ/nYfgGGINwAYWLPMKGhGKaye278zJZChCPcAGFiNxRTuYGpfNyWgiEIN0CYOFPB6bmBoQg3MAXhBggTp3JDuIGhmAoOUxBugDAppnIDw7GIH0xBuAHCxG4o9hFuYCg2zoQpCDdAmNibDRJuYKoYKjcwBOEGCBOncsNsKRjKyzo3MAThBggT+69dKjcwVYz3xK8Mwg0iHeEGCJMA4QaGO5ltCDeIeIQbIEyo3MB0VG5gCsINECZ+wg0M51RumC2FCEe4AcLECTc0FMNQ9veuZbGQHyIb4QYIEz/r3MBw9m0pieoNIhvhBggTP+vcwHAlsg19N4hohBsgTKjcwHRBlRvCDSIY4QYIExqKYbqSlRtWKUYkI9wAYUK4gelKVm5oKEYki4hwM3XqVKWmpio+Pl5paWlas2bNWc+fP3++mjZtqvj4eLVo0UKLFi0K00iB8mO2FExXMpdTuUEkcz3czJs3T1lZWRo3bpzWr1+vVq1aKTMzU3v27Cn1/FWrVql3794aNGiQNmzYoJ49e6pnz57auHFjmEcOhMYJNz7CDczk8XjYGRxG8FiWu9+haWlpuvbaa/V///d/kqRAIKCUlBT99re/1cMPP3za+b169VJhYaHefvtt59h1112n1q1ba9q0aed8vYKCAtWsWVOHDh1SYmJihV1HUbFfew8XVdjzIfo8+97nmr/ua93X6Qo9+NOmbg8HKJcrH3lHx/wBvXHv9UpKjHN7OIhQVWK8qlsjvkKfM5Tf3zEV+sohOnbsmNatW6cxY8Y4x7xerzIyMrR69epSH7N69WplZWUFHcvMzNRbb71V6vlFRUUqKjoVOgoKCs5/4KX49JsC/ez5VZXy3Igu9NzAZD6vR/JLP3+B9zuc2TWXXqR/3NfBtdd3Ndzs27dPfr9fSUlJQceTkpK0adOmUh+Tl5dX6vl5eXmlnp+dna3x48dXzIDPwiMpLsb1u3yIcDXiY/XjKy9xexhAufVoXV9vbtjl9jAQ4WJ97v4+dDXchMOYMWOCKj0FBQVKSUmp8Ndpc2ktbX7ylgp/XgCIJH/4eUv94ect3R4GcFauhps6derI5/MpPz8/6Hh+fr7q1atX6mPq1asX0vlxcXGKi+O+MAAAFwpX60ZVqlRR27ZtlZOT4xwLBALKyclRenp6qY9JT08POl+SlixZcsbzAQDAhcX121JZWVnq37+/2rVrp/bt22vKlCkqLCzUgAEDJEn9+vVTgwYNlJ2dLUm6//771bFjR02cOFFdu3bV3LlztXbtWk2fPt3NywAAABHC9XDTq1cv7d27V2PHjlVeXp5at26txYsXO03DO3bskLfEqpjXX3+95syZo0cffVS/+93v9KMf/UhvvfWWmjdv7tYlAACACOL6OjfhVlnr3AAAgMoTyu9v5i4DAICoQrgBAABRhXADAACiCuEGAABEFcINAACIKoQbAAAQVQg3AAAgqhBuAABAVCHcAACAqOL69gvhZi/IXFBQ4PJIAABAWdm/t8uyscIFF24OHz4sSUpJSXF5JAAAIFSHDx9WzZo1z3rOBbe3VCAQ0DfffKMaNWrI4/FU6HMXFBQoJSVFO3fujMp9q6L9+qTov0auz3zRfo1cn/kq6xoty9Lhw4dVv379oA21S3PBVW68Xq8aNmxYqa+RmJgYtd+0UvRfnxT918j1mS/ar5HrM19lXOO5KjY2GooBAEBUIdwAAICoQripQHFxcRo3bpzi4uLcHkqliPbrk6L/Grk+80X7NXJ95ouEa7zgGooBAEB0o3IDAACiCuEGAABEFcINAACIKoQbAAAQVQg3ZfT++++rW7duql+/vjwej956662znv+Pf/xDN998sy655BIlJiYqPT1d//73v8Mz2HIK9RpXrlypDh066OKLL1ZCQoKaNm2qyZMnh2ew5RDq9ZX0wQcfKCYmRq1bt6608Z2vUK9v2bJl8ng8p33k5eWFZ8DlUJ5/w6KiIj3yyCO67LLLFBcXp9TUVM2aNavyB1sOoV7fr371q1L/Da+++urwDLgcyvNvOHv2bLVq1UpVq1ZVcnKyBg4cqG+//bbyB1sO5bm+qVOn6qqrrlJCQoKaNGmiV155pfIHWk7Z2dm69tprVaNGDdWtW1c9e/bU5s2bz/m4+fPnq2nTpoqPj1eLFi20aNGiSh0n4aaMCgsL1apVK02dOrVM57///vu6+eabtWjRIq1bt06dO3dWt27dtGHDhkoeafmFeo3VqlXTsGHD9P777+uzzz7To48+qkcffVTTp0+v5JGWT6jXZzt48KD69eunn/zkJ5U0sopR3uvbvHmzdu/e7XzUrVu3kkZ4/spzjXfeeadycnI0c+ZMbd68Wa+++qqaNGlSiaMsv1Cv79lnnw36t9u5c6dq166tO+64o5JHWn6hXuMHH3ygfv36adCgQfr00081f/58rVmzRoMHD67kkZZPqNf3wgsvaMyYMXr88cf16aefavz48Ro6dKj+9a9/VfJIy2f58uUaOnSoPvzwQy1ZskTHjx9Xly5dVFhYeMbHrFq1Sr1799agQYO0YcMG9ezZUz179tTGjRsrb6AWQibJevPNN0N+XLNmzazx48dX/IAqQXmv8fbbb7fuuuuuih9QBQvl+nr16mU9+uij1rhx46xWrVpV6rgqSlmub+nSpZYk68CBA2EZU0UryzW+8847Vs2aNa1vv/02PIOqQOX5GXzzzTctj8djbd++vXIGVcHKco1//OMfrcsvvzzo2HPPPWc1aNCgEkdWMcpyfenp6daoUaOCjmVlZVkdOnSoxJFVnD179liSrOXLl5/xnDvvvNPq2rVr0LG0tDTrN7/5TaWNi8pNmAQCAR0+fFi1a9d2eyiVZsOGDVq1apU6duzo9lAqzF//+ldt27ZN48aNc3solaZ169ZKTk7WzTffrA8++MDt4VSoBQsWqF27dnrmmWfUoEEDXXnllRo1apS+++47t4dWKWbOnKmMjAxddtllbg+lwqSnp2vnzp1atGiRLMtSfn6+Xn/9dd16661uD61CFBUVKT4+PuhYQkKC1qxZo+PHj7s0qrI7dOiQJJ31d9vq1auVkZERdCwzM1OrV6+utHERbsLkT3/6k44cOaI777zT7aFUuIYNGyouLk7t2rXT0KFDdc8997g9pArx+eef6+GHH9bf//53xcRE3x6zycnJmjZtmt544w298cYbSklJUadOnbR+/Xq3h1Zhtm3bppUrV2rjxo168803NWXKFL3++uu677773B5ahfvmm2/0zjvvRM3Pn61Dhw6aPXu2evXqpSpVqqhevXqqWbNmyLdfI1VmZqZefPFFrVu3TpZlae3atXrxxRd1/Phx7du3z+3hnVUgENCIESPUoUMHNW/e/Izn5eXlKSkpKehYUlJSpfb3Rd87dgSaM2eOxo8fr3/+858R3c9QXitWrNCRI0f04Ycf6uGHH1bjxo3Vu3dvt4d1Xvx+v/r06aPx48fryiuvdHs4laJJkyZBvSfXX3+9vvjiC02ePFl/+9vfXBxZxQkEAvJ4PJo9e7azm/CkSZP0i1/8Qs8//7wSEhJcHmHFefnll3XRRRepZ8+ebg+lQv3vf//T/fffr7FjxyozM1O7d+/W6NGjNWTIEM2cOdPt4Z23xx57THl5ebruuutkWZaSkpLUv39/PfPMM/J6I7v+MHToUG3cuFErV650eyinIdxUsrlz5+qee+7R/PnzTyvLRYtGjRpJklq0aKH8/Hw9/vjjxoebw4cPa+3atdqwYYOGDRsm6cQvSsuyFBMTo3fffVc33XSTy6OseO3bt4/IN6rySk5OVoMGDZxgI0lXXXWVLMvS119/rR/96Ecujq7iWJalWbNm6e6771aVKlXcHk6Fys7OVocOHTR69GhJUsuWLVWtWjXdeOONevLJJ5WcnOzyCM9PQkKCZs2apb/85S/Kz89XcnKypk+frho1auiSSy5xe3hnNGzYML399tt6//331bBhw7OeW69ePeXn5wcdy8/PV7169SptfJEdCw336quvasCAAXr11VfVtWtXt4cTFoFAQEVFRW4P47wlJibqk08+UW5urvMxZMgQNWnSRLm5uUpLS3N7iJUiNzfX+F8WJXXo0EHffPONjhw54hzbsmWLvF7vOd+QTbJ8+XJt3bpVgwYNcnsoFe7o0aOnVTB8Pp+kE6EuWsTGxqphw4by+XyaO3eubrvttois3FiWpWHDhunNN9/Uf/7zH+eP27NJT09XTk5O0LElS5YoPT29soZJ5aasjhw5oq1btzqff/nll8rNzVXt2rV16aWXasyYMdq1a5ezPsGcOXPUv39/Pfvss0pLS3PuLSYkJAT9FRlJQr3GqVOn6tJLL1XTpk0lnZj+/qc//UnDhw93ZfznEsr1eb3e0+4h161bV/Hx8We9t+ymUP/9pkyZokaNGunqq6/W999/rxdffFH/+c9/9O6777p1CecU6jX26dNHTzzxhAYMGKDx48dr3759Gj16tAYOHBiRt6RCvT7bzJkzlZaWFrHfmyWFeo3dunXT4MGD9cILLzi3pUaMGKH27durfv36bl3GGYV6fVu2bNGaNWuUlpamAwcOaNKkSdq4caNefvllty7hrIYOHao5c+bon//8p2rUqOH8bqtZs6bzM9WvXz81aNBA2dnZkqT7779fHTt21MSJE9W1a1fNnTtXa9eurdxlQyptHlaUsafN/vCjf//+lmVZVv/+/a2OHTs653fs2PGs50eiUK/xueees66++mqratWqVmJiotWmTRvr+eeft/x+vzsXcA6hXt8PRfpU8FCv7+mnn7auuOIKKz4+3qpdu7bVqVMn6z//+Y87gy+j8vwbfvbZZ1ZGRoaVkJBgNWzY0MrKyrKOHj0a/sGXQXmu7+DBg1ZCQoI1ffr08A+4HMpzjc8995zVrFkzKyEhwUpOTrb69u1rff311+EffBmEen3/+9//rNatW1sJCQlWYmKi1aNHD2vTpk3uDL4MSrs2SdZf//pX55yOHTue9rvutddes6688kqrSpUq1tVXX20tXLiwUsfpOTlYAACAqBB5N/QAAADOA+EGAABEFcINAACIKoQbAAAQVQg3AAAgqhBuAABAVCHcAACAqEK4AQAAUYVwAwAAogrhBgAARBXCDQBXderUSSNGjHB7GI7yjufbb79V3bp1tX379gof0w/98pe/1MSJEyv9dQBTEW6AC8C0adNUo0YNFRcXO8eOHDmi2NhYderUKejcZcuWyePx6IsvvgjzKMOrokPVU089pR49eig1NbXCnvNMHn30UT311FM6dOhQpb8WYCLCDXAB6Ny5s44cOaK1a9c6x1asWKF69erpo48+0vfff+8cX7p0qS699FJdccUVbgzVSEePHtXMmTM1aNCgsLxe8+bNdcUVV+jvf/97WF4PMA3hBrgANGnSRMnJyVq2bJlzbNmyZerRo4caNWqkDz/8MOh4586dJUmLFy/WDTfcoIsuukgXX3yxbrvttqCKzvTp01W/fn0FAoGg1+vRo4cGDhwoSQoEAsrOzlajRo2UkJCgVq1a6fXXXz/jWMtyfqdOnTR8+HA9+OCDql27turVq6fHH3886JzDhw+rb9++qlatmpKTkzV58mSnWvOrX/1Ky5cv17PPPiuPxyOPxxN0OykQCJz1uX9o0aJFiouL03XXXRd0fOXKlYqNjQ0Kj9u3b5fH49FXX32lTp066be//a1GjBihWrVqKSkpSTNmzFBhYaEGDBigGjVqqHHjxnrnnXdOe81u3bpp7ty5Zx0XcKEi3AAXiM6dO2vp0qXO50uXLlWnTp3UsWNH5/h3332njz76yAk3hYWFysrK0tq1a5WTkyOv16vbb7/dCTN33HGHvv3226Dn3b9/vxYvXqy+fftKkrKzs/XKK69o2rRp+vTTTzVy5EjdddddWr58eanjLOv5L7/8sqpVq6aPPvpIzzzzjCZMmKAlS5Y4X8/KytIHH3ygBQsWaMmSJVqxYoXWr18vSXr22WeVnp6uwYMHa/fu3dq9e7dSUlLK/Nw/tGLFCrVt2/a047m5ubrqqqsUHx/vHNuwYYNq1aqlyy67zHmtOnXqaM2aNfrtb3+re++9V3fccYeuv/56rV+/Xl26dNHdd9+to0ePBj13+/bttWbNGhUVFZ1xXMAFywJwQZgxY4ZVrVo16/jx41ZBQYEVExNj7dmzx5ozZ4714x//2LIsy8rJybEkWV999VWpz7F3715LkvXJJ584x3r06GENHDjQ+fwvf/mLVb9+fcvv91vff/+9VbVqVWvVqlVBzzNo0CCrd+/elmVZVseOHa3777/fsiyrTOfbj7nhhhuCzrn22muthx56yLIsyyooKLBiY2Ot+fPnO18/ePCgVbVqVee1Sr5uSed67tL88P+B7Z577rH69esXdGzs2LFWp06dSn2t4uJiq1q1atbdd9/tHNu9e7clyVq9enXQ83z88ceWJGv79u1nHBdwoaJyA1wgOnXqpMLCQv33v//VihUrdOWVV+qSSy5Rx44dnb6bZcuW6fLLL9ell14qSfr888/Vu3dvXX755UpMTHSaZXfs2OE8b9++ffXGG284FYTZs2frl7/8pbxer7Zu3aqjR4/q5ptvVvXq1Z2PV155pdSG5VDOb9myZdDnycnJ2rNnjyRp27ZtOn78uNq3b+98vWbNmmrSpEmZ/l+d7blL89133wVVZ2y5ublq3bp10LENGzYEHSv5Wj6fTxdffLFatGjhHEtKSpKk014/ISFBkk6r6ACQYtweAIDwaNy4sRo2bKilS5fqwIED6tixoySpfv36SklJ0apVq7R06VLddNNNzmO6deumyy67TDNmzHB6a5o3b65jx44FnWNZlhYuXKhrr71WK1as0OTJkyWdmJElSQsXLlSDBg2CxhMXF3faGEM5PzY2Nuhzj8dzWu9PeYX63HXq1NGBAweCjvn9fm3cuFFt2rQJOr5+/Xr9/Oc/P+trlTzm8Xgk6bTX379/vyTpkksuOdflABccwg1wAencubOWLVumAwcOaPTo0c7xH//4x3rnnXe0Zs0a3XvvvZJOrNuyefNmzZgxQzfeeKOkEw2yPxQfH6+f/exnmj17trZu3aomTZrommuukSQ1a9ZMcXFx2rFjhxOmzibU88/k8ssvV2xsrP773/86VahDhw5py5Yt+vGPfyxJqlKlivx+f7lfo6Q2bdqcNnNp8+bN+v7771W/fn3n2OrVq7Vr167TqjnlsXHjRjVs2FB16tQ57+cCog3hBriAdO7cWUOHDtXx48eDwkPHjh01bNgwHTt2zGkmrlWrli6++GJNnz5dycnJ2rFjhx5++OFSn7dv37667bbb9Omnn+quu+5yjteoUUOjRo3SyJEjFQgEdMMNN+jQoUP64IMPlJiYqP79+wc9T6jnn0mNGjXUv39/jR49WrVr11bdunU1btw4eb1epxKSmpqqjz76SNu3b1f16tVVu3Zteb3lu1OfmZmpMWPG6MCBA6pVq5akE7ekJOnPf/6zhg8frq1bt2r48OGSFFT5Kq8VK1aoS5cu5/08QDSi5wa4gHTu3FnfffedGjdu7PRySCfCzeHDh50p45Lk9Xo1d+5crVu3Ts2bN9fIkSP1xz/+sdTnvemmm1S7dm1t3rxZffr0CfraE088occee0zZ2dm66qqr9NOf/lQLFy5Uo0aNSn2uUM8/k0mTJik9PV233XabMjIy1KFDh6CZS6NGjZLP51OzZs10ySWXBPURhapFixa65ppr9NprrznHcnNzlZmZqW3btqlFixZ65JFHNH78eCUmJuq5554r92tJ0vfff6+33npLgwcPPq/nAaKVx7Isy+1BAEBlKywsVIMGDTRx4sRKWWxv4cKFGj16tDZu3Civ16vMzExde+21evLJJyv8tV544QW9+eabevfddyv8uYFowG0pAFFpw4YN2rRpk9q3b69Dhw5pwoQJkk4sMFgZunbtqs8//1y7du1SSkqKPv74Y2chw4oWGxurP//5z5Xy3EA0oHIDICpt2LBB99xzjzZv3qwqVaqobdu2mjRpUtA068qSl5en5ORkffrpp2rWrFmlvx6AYIQbAAAQVWgoBgAAUYVwAwAAogrhBgAARBXCDQAAiCqEGwAAEFUINwAAIKoQbgAAQFQh3AAAgKhCuAEAAFGFcAMAAKLK/wclMjUZwo5HhAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Proteus I/F at 1.60 um = 0.11\n"
     ]
    }
   ],
   "source": [
    "observed_flux = 7e-14 #erg s-1 cm-2 um-1\n",
    "neptune_earth_dist = 4.6e9 #km\n",
    "neptune_sun_dist = 4.3e9 #km\n",
    "\n",
    "r = 210 #km; radius of Proteus\n",
    "dS = np.pi*r**2\n",
    "\n",
    "# find the solid angle of Proteus\n",
    "Omega = solid_angle(dS, neptune_earth_dist)\n",
    "print(f'Proteus solid angle is {Omega:.2e} sr = {Omega*4.25e10:.6f} arcsec^2')\n",
    "\n",
    "# example transmission curve for roughly H-band\n",
    "wls = np.linspace(1.2, 2.0, 100) #microns\n",
    "trans = np.zeros(wls.shape)\n",
    "trans[wls < 1.75] = 1.0\n",
    "trans[wls < 1.45] = 0.0\n",
    "bp = np.array([wls, trans])\n",
    "\n",
    "plt.plot(wls, trans)\n",
    "plt.xlabel(r'Wavelength ($\\mu$m)')\n",
    "plt.ylabel('Transmission')\n",
    "plt.show()\n",
    "\n",
    "# find I/F in that filter\n",
    "wl_eff, IoverF = I_over_F(observed_flux, bp, neptune_sun_dist/1.5e8, Omega)\n",
    "print(f'Proteus I/F at {wl_eff:.2f} um = {IoverF:.2f}')"
   ]
  }
 ],
 "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}