HajimeKawahara/exojax

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison the AM geometric opacity with PieMieScatt "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We import `pdb` with `MgSiO3` to use the refractive index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/kawahara/exojax/documents/tutorials/.database/particulates/virga/virga.zip  exists. Remove it if you wanna re-download and unzip.\n",
      "Refractive index file found:  /home/kawahara/exojax/documents/tutorials/.database/particulates/virga/MgSiO3.refrind\n",
      "Miegrid file exists: /home/kawahara/exojax/documents/tutorials/.database/particulates/virga/miegrid_lognorm_MgSiO3.mg.npz\n"
     ]
    }
   ],
   "source": [
    "from exojax.spec.pardb import PdbCloud\n",
    "miedir = \"/home/kawahara/exojax/documents/tutorials/.database/particulates/virga\"\n",
    "#miedir = \"/home/exoplanet01/exojax/documents/tutorials/.database/particulates/virga\"\n",
    "pdb = PdbCloud(\"MgSiO3\", path=miedir)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pdb` has the information of refractive index in `pdb.refraction_index` and `pdb.refraction_index_wavelength_nm`.\n",
    "These values were taken from `VIRGA`.\n",
    "\n",
    "Notice the refactive index has the form of $m = n + ik$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHLCAYAAAA0kLlRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB4YklEQVR4nO3dd3hTZfsH8G+SNt17L2ihUFZZZRUEylBERYYDEQQBceFABBV/iuLiVWT4qq+IylRERHCAskH23nvTFjrp3mlyfn+cJjR0kLRJTpp8P9fVq8nJc865KWl79xn3IxMEQQARERGRjZBLHQARERGRKTG5ISIiIpvC5IaIiIhsCpMbIiIisilMboiIiMimMLkhIiIim8LkhoiIiGwKkxsiIiKyKUxuiIiIyKYwuSEimyeTyfD+++9LHQYRWQiTGyIyi8WLF0Mmk0Emk2HXrl1VXhcEAREREZDJZHjooYfqdI9du3Zh4MCBCAsLg7OzMxo1aoRBgwZh+fLlBp1/+vRpjBo1CmFhYXByckJoaChGjhyJ06dPV9v2scceQ5MmTeDq6gp/f3/06tULf/31V51iJyLzcZA6ACKybc7Ozli+fDnuueceveP//vsvkpOT4eTkVKfr/vrrrxg+fDjat2+PV199FT4+Prh69Sp27NiB7777Dk8++aSubXFxMRwc9H/crV69GiNGjICvry/Gjx+PqKgoXLt2DT/88ANWrVqFFStWYOjQobr2169fR35+PsaMGYPQ0FAUFRXht99+w8MPP4xvv/0Wzz77bJ3+HURkejJunElE5rB48WKMHTsWw4YNw44dO5CSkqKXYDz77LM4cuQIMjMz0aZNG6xdu9ao67du3RoymQxHjhyBUqnUey09PR2BgYE1nnv58mW0bdsWjRo1wo4dOxAQEKB7LTMzEz179kRSUhJOnDiBJk2a1HgdtVqNuLg4lJSU4Ny5c0bFT0Tmw2EpIjKrESNG4NatW9i0aZPuWFlZGVatWqXXu6J169YtPPXUU/D09IS3tzfGjBmD48ePQyaTYfHixbp2ly9fRufOnaskNgCqJDZ3zrmZNWsWioqKsGDBAr3EBgD8/f3x7bfforCwEJ999lmt/zaFQoGIiAjk5OTU2o6ILIvJDRGZVWRkJOLj4/Hzzz/rjv3zzz/Izc3FE088oddWo9Fg0KBB+PnnnzFmzBh8/PHHSElJwZgxY6pct3HjxtiyZQuSk5ONjumvv/5CZGQkevbsWe3rvXr1QmRkJNatW1fltcLCQmRmZuLy5cuYO3cu/vnnH/Tr18/oGIjIfJjcEJHZPfnkk/j9999RXFwMAPjpp5/Qu3dvhIaG6rX7/fffsXfvXsyePRtffvklJk6ciPXr18PPz6/KNd98800kJSWhadOm6Nu3L6ZPn45du3ZBo9HUGktubi5u3ryJdu3a1dqubdu2SE5ORn5+vt7x119/HQEBAYiOjsaUKVMwdOhQfPXVV4Z8GYjIQpjcEJHZPf744yguLsbatWuRn5+PtWvXVjsktX79ejg6OmLChAm6Y3K5HBMnTqzSdty4cVi/fj0SEhKwa9cufPjhh+jZsyeaNWuGPXv21BiLNlnx8PCoNWbt63l5eXrHJ02ahE2bNmHJkiUYOHAg1Go1ysrKar0WEVkWkxsiMruAgAD0798fy5cvx+rVq6FWq/Hoo49WaXf9+nWEhITA1dVV73h0dHS11x0wYAA2bNiAnJwc7NixAxMnTsT169fx0EMPIT09vdpztEnLnT0yd6opCWrRogX69++P0aNHY+3atSgoKMCgQYPAtRlE1oPJDRFZxJNPPol//vkH8+fPx8CBA+Ht7W2ya7u6uqJnz5746quv8M477yA7Oxv//PNPtW29vLwQEhKCEydO1HrNEydOICwsDJ6enrW2e/TRR3Hw4EFcuHChzvETkWkxuSEiixg6dCjkcjn27dtX7ZAUIE4STklJQVFRkd7xS5cuGXyfTp06AQBSUlJqbPPQQw/h6tWr1RYXBICdO3fi2rVrBhUX1M4jys3NNThGIjIvJjdEZBHu7u745ptv8P7772PQoEHVthkwYABUKhW+++473TGNRoOvv/66StstW7ZUe42///4bABATE1NjLFOnToWLiwuee+453Lp1S++1rKwsPP/883B1dcXUqVN1x6sb5lKpVFi6dClcXFzQqlWrGu9HRJbFCsVEZDHVLemubMiQIejSpQtef/11XLp0CS1atMCff/6JrKwsAGK9Gq3BgwcjKioKgwYNQtOmTVFYWIjNmzfjr7/+QufOnWtMoACgWbNmWLJkCUaOHInY2NgqFYozMzPx888/o2nTprpznnvuOeTl5aFXr14ICwtDamoqfvrpJ5w7dw6zZ8+Gu7t7Pb86RGQqTG6IyGooFAqsW7cOr776KpYsWQK5XI6hQ4fivffeQ48ePeDs7Kxr+/333+OPP/7AypUrcfPmTQiCgCZNmuD//u//8Oabb1bZbuFOjz32GFq0aIGZM2fqEho/Pz/06dMHb7/9Ntq0aaPXfvjw4fjhhx/wzTff4NatW/Dw8EBcXBw+/fRTPPzww2b5ehBR3XD7BSKyer///juGDh2KXbt2oUePHlKHQ0RWjskNEVmV4uJiuLi46J6r1Wrcd999OHToEFJTU/VeIyKqDoeliMiqvPzyyyguLkZ8fDxKS0uxevVq7NmzB5988gkTGyIyCHtuiMiqLF++HLNnz8alS5dQUlKC6OhovPDCC3jppZekDo2IGggmN0RERGRTWOeGiIiIbAqTGyIiIrIpdjehWKPR4ObNm/Dw8NArCEZERETWSxAE5OfnIzQ0FHJ57X0zdpfc3Lx5ExEREVKHQURERHWQlJSE8PDwWtvYXXLj4eEBQPzi3G23XyIiIrIOeXl5iIiI0P0er43dJTfaoShPT08mN0RERA2MIVNKOKGYiIiIbAqTGyIiIrIpdjcsZSi1Wg2VSiV1GGRjlErlXWf5ExFR/TC5uYMgCEhNTUVOTo7UoZANksvliIqKglKplDoUIiKbxeTmDtrEJjAwEK6urqyFQyajrbGUkpKCRo0a8b1FRGQmTG4qUavVusTGz89P6nDIBgUEBODmzZsoLy+Ho6Oj1OEQEdkkDv5Xop1j4+rqKnEkZKu0w1FqtVriSIiIbBeTm2pwuIDMhe8tIiLzY3JDRERENoXJDREREdkUJjdUq2vXrkEmk+HYsWNSh0JERGQQJjdk1Z5++mkMGTJE6jCIiKgB4VJwG1VWVtagC8Wp1WpOviWiBkOtEaDWCFA6iH0GuUUq7Lt6C3nFKuSVlFd8ViGvuBzFqnIMaB2Mwe3DAADpeSWY/sdp3bUECHrX7tsiEMM7NwIA5BSV4Y1VJyq11dezmT9Gx0cCAIrKyvHKz0dvt72jcZcoXzzXuykAoFytwYSlh2q8bocIH7zav5nu+fjFB1GuEapt2zrUE2/e3wJSYnJjoKKy8hpfk8tkcHZUmLStq9K4/5qEhAS0adMGDg4O+PHHHxEbG4tt27bh1KlTmDp1Knbu3Ak3Nzfcd999mDt3Lvz9/QEA69evx0cffYRTp05BoVAgPj4eX3zxBZo2bWrwvSMjIzF+/HicOXMGf/75J7y9vfH2229j4sSJujZz5szBokWLcOXKFfj6+mLQoEH47LPP4O7uDgBYvHgxJk2ahKVLl+Ktt97ChQsXMGrUKCxZsgTA7VVG27ZtQ0JCglFfGyKyL+l5Jdh3NQuCIEAjCNBoIH4WBGgEoGMjH8QEewAA0vJK8PfJFGgE3G4viO0FAejWxBdxjX0BAFczC/HfLReRX5GkiMmKmLwUlJbjjftj8GJCNAAgKbsIzy07XGOMkX5uuseFZWqsP51aY9sQLxfd49JyDTaeSauxrZ/b7T9qyzUCNp9Nr7Gtm9Pt3zMCgG3nM2ps63DHtjE7LmZApb4zrRFpNNUftyQmNwZqNX1Dja/1iQnAorFddM/jPtyMYlX1dUy6Rvnil+fidc/v+XQbsgrLqrS79p8HjY5xyZIleOGFF7B7924AQE5ODvr27YtnnnkGc+fORXFxMd588008/vjj2Lp1KwCgsLAQkydPRtu2bVFQUIDp06dj6NChOHbsmFF7IM2aNQtvv/02ZsyYgQ0bNuDVV19F8+bNce+99wIQtx3473//i6ioKFy5cgUvvvgi3njjDfzvf//TXaOoqAiffvopvv/+e/j5+SEkJATFxcXIy8vDokWLAAC+vr5Gf12IyHapNQKOJWUjyNMZ4T5ijbKzqfl6PRZ3mv5QK11yk5hVhBl/namx7dQBMbrkprC0HGuO3qixbV7x7T9W/dyV6NjIG54ujvB0doSHs4PusatSgdhwL11bXzclPhzSRu9alfutW1TECgAezg74ZGisfttKjZv4306anB0U+PSRO9pWunKE7+2abgqZDJ892rbGGEK9XfRe+/SRtno9QZVjCPBwgtSY3NiQZs2a4bPPPtM9/+ijj9ChQwd88sknumMLFy5EREQELly4gObNm+ORRx7Ru8bChQsREBCAM2fOoE0b/W+22vTo0QNvvfUWAKB58+bYvXs35s6dq0tuJk2apGsbGRmJjz76CM8//7xecqNSqfC///0P7dq10x1zcXFBaWkpgoODDY6FiGxbbpEK/17MwLZz6dh+Ph3ZRSq82q8ZXru3OQDAx9UR3Zv6QS6TQSYTe8zlFZ9lMpneL3VfNyUGtQut9Lp++1Yhnrq2Yd4umDawhS5J8XRxqPjsCE9nB3g43646HuLlgtUv9jDo3+Pl4oinujU2qK2r0gFPdm1kUFulg1w3nHU3crkMj3eKMKgtAAzrGG5wWykwuTHQmQ8G1Pia/I65IYff7W9w211v9qlfYJXExcXpPT9+/Di2bdumG/qp7PLly2jevDkuXryI6dOnY//+/cjMzIRGowEAJCYmGpXcxMfHV3k+b9483fPNmzdj5syZOHfuHPLy8lBeXo6SkhIUFRXpKkIrlUq0bav/lwMRESAO4S/bex1bzqXj8PVsqCsNfXg6O+jN+2gb7o3lE7oZdN2mAe74ckQHg9r6uCl1c1TIukma3MycOROrV6/GuXPn4OLigu7du+PTTz9FTExMjecsXrwYY8eO1Tvm5OSEkpISs8ZqzBwYc7W9Gzc3N73nBQUFGDRoED799NMqbUNCQgAAgwYNQuPGjfHdd98hNDQUGo0Gbdq0QVlZ1aGyurp27RoeeughvPDCC/j444/h6+uLXbt2Yfz48SgrK9MlNy4uLpxETEQAgBKVGsnZRYgOFIdkHBVyfLXtEvJLxKGfZoHu6NsyEP1aBKFjI284KLj4l26TNLn5999/MXHiRHTu3Bnl5eV4++23cd999+HMmTNVflFX5unpifPnz+ue8xdi9Tp27IjffvsNkZGRcHCo+l9969YtnD9/Ht999x169uwJANi1a1ed7rVv374qz1u2bAkAOHz4MDQaDWbPnq2bx7Ny5UqDrqtUKrkPE5Gd0GgE/H0qBX8cu4ldFzPh567Ezjf6QCaTwVEhx4sJ0XBVKtC3RaDe0BLRnSRNbtavX6/3fPHixQgMDMThw4fRq1evGs+TyWScg2GAiRMn4rvvvsOIESPwxhtvwNfXF5cuXcKKFSvw/fffw8fHB35+fliwYAFCQkKQmJiomzdjrN27d+Ozzz7DkCFDsGnTJvz6669Yt24dACA6OhoqlQpffvklBg0ahN27d2P+/PkGXTcyMhIbNmzA+fPn4efnBy8vL+6mTWRjBEHApjNpmLPpAs6l5uuOq9QaZBaU6SaovpDAISEyjFX14+Xm5gK4+4qYgoICNG7cGBERERg8eDBOnz5dY9vS0lLk5eXpfdiL0NBQ7N69G2q1Gvfddx9iY2MxadIkeHt7Qy6XQy6XY8WKFTh8+DDatGmD1157DbNmzarTvV5//XUcOnQIHTp0wEcffYQ5c+ZgwABxnlK7du0wZ84cfPrpp2jTpg1++uknzJw506DrTpgwATExMejUqRMCAgJ0K8GIyDYcTczGkP/twbPLDuNcaj48nB3wUp9orH35Huyb1s8qVt5QwyMThDvL+khDo9Hg4YcfRk5OTq1DI3v37sXFixfRtm1b5Obm4vPPP8eOHTtw+vRphIdXnb39/vvvY8aMGVWO5+bmwtPTU+9YSUkJrl69iqioKDg7O9f/H2UnIiMjMWnSJL0VUVQ9vseI9B24moXHv90LV6UCY3tE4tmeTeHlyt5ZqiovLw9eXl7V/v6+k9Wslpo4cSJOnTp11zkf8fHxeitzunfvjpYtW+Lbb7/Fhx9+WKX9tGnTMHnyZN3zvLw8REQYvtyNiIhM59SNXJxPzccjceIfo12ifPHhkDa4v3Uwe2nIZKwiuXnppZewdu1a7Nixo9rel9o4OjqiQ4cOuHTpUrWvOzk5wcmJ3zBERFK6mJaPOZsu4J9TqXBxVKBnc38Eeoi9l4bWeCEylKTJjSAIePnll7FmzRps374dUVFRRl9DrVbj5MmTeOCBB8wQIRni2rVrUodARFZKEAS8/+dpLN13HYIgVrId0DoIFSW1iMxC0uRm4sSJWL58Of744w94eHggNVXcW8PLywsuLmKp59GjRyMsLEw3AfWDDz5At27dEB0djZycHMyaNQvXr1/HM888I9m/g4iIqrf6yA0s2XsdAHB/62C8dm9z3dYHROYiaXLzzTffAECVjRAXLVqEp59+GoBYKbfyHkfZ2dmYMGECUlNT4ePjg7i4OOzZswetWrWyVNhERGSA7MIyfPz3WQDiHk0T+0RLHBHZC8mHpe5m+/btes/nzp2LuXPnmikiIiIyla3n0pFVWIbmQe54tlcTqcMhO2IVE4qJiMj2PBIXjghfVygd5HDk9ghkQUxuiIjIbLpE1V6UlcgcmEoTEZFJrTuRgqSsIqnDIDvG5MZGJCQkWE2F4MjISMybN0/qMIhIAtdvFWLyymO4d+6/uJiWf/cTiMyAw1I2YvXq1VazoeTBgwdr3dXdlslkMqxZswZDhgyROhQiixMEAe/8fgql5Rr0iPZDdKC71CGRnWJyYyPuttmoJQUEBJj9HmVlZVAqlWa/j6GsLR4iKfx1IgU7L2ZC6SDHh4PbQCaTSR0S2SkOS92NIABlhdJ8GLGn6Z3DUpGRkfjoo48wevRouLu7o3Hjxvjzzz+RkZGBwYMHw93dHW3btsWhQ4d059y6dQsjRoxAWFgYXF1dERsbi59//lnvPvn5+Rg5ciTc3NwQEhKCuXPnVnvvysNSMpkM33//PYYOHQpXV1c0a9YMf/75p+51tVqN8ePHIyoqCi4uLoiJicEXX3yhd9+nn34aQ4YMwccff4zQ0FDExMTggw8+QJs2bap8Ldq3b49333232q/T9u3bIZPJsG7dOrRt2xbOzs7o1q0bTp06ZdTXISEhAS+99BImTZoEf39/DBgwAJGRkQCAoUOHQiaT6Z4T2YPcYhU++OsMAGBiQjSaBLDXhqTDnpu7URUBn4RKc++3bwLKug/vzJ07F5988gneffddzJ07F0899RS6d++OcePGYdasWXjzzTcxevRonD59GjKZDCUlJYiLi8Obb74JT09PrFu3Dk899RSaNm2KLl26AAAmT56M3bt3488//0RQUBCmT5+OI0eOoH379rXGMmPGDHz22WeYNWsWvvzyS4wcORLXr1+Hr68vNBoNwsPD8euvv8LPzw979uzBs88+i5CQEDz++OO6a2zZsgWenp7YtGkTALGS9YwZM3Dw4EF07twZAHD06FGcOHECq1evrjWeqVOn4osvvkBwcDDefvttDBo0CBcuXICjo6NBXwcAWLJkCV544QXs3r0bgNh7FhgYiEWLFuH++++HQqEw+v+MqKGateEcMgtK0STADc8nsKYNSYvJjQ174IEH8NxzzwEApk+fjm+++QadO3fGY489BgB48803ER8fj7S0NAQHByMsLAxTpkzRnf/yyy9jw4YNWLlyJbp06YL8/HwsWbIEy5cvR79+/QCI1aRDQ++e/D399NMYMWIEAOCTTz7Bf//7Xxw4cAD3338/HB0dMWPGDF3bqKgo7N27FytXrtRLbtzc3PD999/rDf8MGDAAixYt0iU3ixYtQu/evdGkSe0/XN977z3ce++9AMQkJTw8HGvWrMHjjz9+16+DVrNmzfDZZ59Vuba3tzeCg4Pv+jUhshXHk3Lw0/5EAMBHQ9rAyYGJPUmLyc3dOLqKPShS3bse2rZtq3scFBQEAIiNja1yLD09HcHBwVCr1fjkk0+wcuVK3LhxA2VlZSgtLYWrqxjHlStXoFKp9H7Be3l5ISYmxqhY3Nzc4OnpifT0dN2xr7/+GgsXLkRiYiKKi4tRVlZWpTcoNja2yryWCRMmYNy4cZgzZw7kcjmWL19uUAXr+Ph43WNfX1/ExMTg7FmxTPzdvg5acXFxd70PkT2ICfbAy32bISO/BN2b+ksdDhGTm7uSyeo1NCSlyquntBP7qjumqdied9asWfjiiy8wb948xMbGws3NDZMmTUJZWZlJY9HeW3vfFStWYMqUKZg9ezbi4+Ph4eGBWbNmYf/+/XrnVLcCa9CgQXBycsKaNWugVCqhUqnw6KOP1itWQ78O9roijOhOzo4KTL63uUFb6hBZApMb0tm9ezcGDx6MUaNGARCTngsXLug2JW3SpAkcHR1x8OBBNGrUCACQm5uLCxcuoFevXvW6b/fu3fHiiy/qjl2+fNmgcx0cHDBmzBgsWrQISqUSTzzxhG5H+drs27dP92/Izs7GhQsX0LJlS108tX0dauPo6Ai1Wm1Q7EQNXW6RCm5OCjhUbK3A1VFkLZjckE6zZs2watUq7NmzBz4+PpgzZw7S0tJ0v9Q9PDwwZswYTJ06VTd59r333oNcLq/XD7VmzZph6dKl2LBhA6KiorBs2TIcPHgQUVFRBp3/zDPP6CUmhvjggw/g5+eHoKAg/N///R/8/f11tWnu9nWoTWRkJLZs2YIePXrAyckJPj4+BsVD1NAIgoDJK48hNa8Enz/WDi1DPKUOiUiHS8FJ55133kHHjh0xYMAAJCQkIDg4uEoxujlz5iA+Ph4PPfQQ+vfvjx49eqBly5Zwdnau832fe+45DBs2DMOHD0fXrl1x69YtvV6cu2nWrBm6d++OFi1aoGvXrgad85///Aevvvoq4uLikJqair/++ks3n8eQr0NNZs+ejU2bNiEiIgIdOnQw+N9A1NBsOJ2KLefScSEtH44K9tiQdZEJdjZImpeXBy8vL+Tm5sLTU/8vjZKSEly9ehVRUVH1+mVtTwoLCxEWFobZs2dj/PjxksQgCAKaNWuGF198EZMnT6617fbt29GnTx9kZ2fD29vbMgFWwvcY2YKC0nL0n/0vUvNK8FKfaEwZcPdFBUT1Vdvv7ztxWIqMcvToUZw7dw5dunRBbm4uPvjgAwDA4MGDJYknIyMDK1asQGpqKsaOHStJDET2ZvbG80jNK0FjP1e81Dda6nCIqmByQ0b7/PPPcf78eSiVSsTFxWHnzp3w95dm+WdgYCD8/f2xYMECzm8hsoBTN3KxZM81AMCHg9vA2ZE1bcj6MLkho3To0AGHDx+WOgwdY0dVExISuFyVqB7mbroAjQAMaheKXs3Nv48cUV1wQjERERmkuEyN9PxSyGXA5HubSx0OUY3Yc1MN/mVP5sL3FjVkLkoF/nypB86m5CPKn0UsyXqx56YSbRXdoqIiiSMhW6WtcsxNNamhkslkaBXKmjZk3dhzU4lCoYC3t7duzyNXV1dW3CST0Wg0yMjIgKurKxwc+K1HDcvG06mIb+oHD2fHuzcmkhh/wt5Bu5tz5U0diUxFLpejUaNGTJqpQbmSUYDnfzwMD2dHbJuSAF835d1PIpIQk5s7yGQyhISEIDAwECqVSupwyMYolUrI5RwNpoblf9svQyMAnRr7MLGhBoHJTQ0UCgXnRRCR3UvKKsKaozcAgAX7qMHgn5BERFSjb/69DLVGQM9m/ujQiIUyqWFgckNERNVKyS3GqkPJAICX+zaTOBoiwzG5ISKian377xWUqTXoEuWLLlG+UodDZDAmN0REVK0SlRpyGfAKe22ogZEJdlYy1Zgt04mI7F1SVhHCfVxYvoAkZ8zvb66WIiKiGkX4ukodApHROCxFRER61p9KxZWMAqnDIKozJjdERKSTW6zC1F+Po/+cf3EsKUfqcIjqhMkNERHpLN1zDfml5YgOdEfbMC+pwyGqEyY3REQEACgoLccPu68CACb2iYZczknE1DAxuSEiIgDAj/uuI6dIhSb+bniobajU4RDVGZMbIiJCcZka3++8AgB4sU80FOy1oQaMyQ0REeHnA4nILChDuI8LBrdnrw01bExuiIgIDgoZvFwc8WJCNBwV/NVADRuL+BEREUbHR2JYx3AomdiQDWByQ0REAAB3J/5KINvAFJ2IyI4duJqFHRcyoNbY1TaDZOOY3BAR2bF5my9g9MIDupVSRLaAyQ0RkZ26kVOMvVduAQAebBsicTREpsPkhojICJfS85GUVSR1GCax5kgyBAGIb+KHcB/u/k22g7PHiIgMlFlQiv5zdgAALn08EA4NeGWRIAj47cgNAMAjceESR0NkWg33O5OIyMKOJuboHp9Py5cuEBM4kpiNq5mFcFUqMLBNsNThEJkUkxsiIgOdupGre1w50WmIVh0We20GtgmBG5eAk41hckNEZKBnekahS5QvALHno6ESBAEXKnqeHokLkzgaItNjuk5EZCAPZ0e8kNAUB65mNeieG5lMhlXPx+PUjTy0DvWUOhwik2NyQ0RkhA4R3mgX4Y2Ojbyh1ggNdvdsmUyG2HAvqcMgMgsmN0REBth35RY2nUlDr+YB+GNiD6nDqbOC0nIoZDK4KBVSh0JkNpxzQ0RkgH8vZOCHXVex/lSq1KHUy4/7rqPzx5vx7b+XpQ6FyGyY3BARGeBEcg4AoG3FUE5RWTnO3MyTMCLjCYKA3w4no6C0HJ4ujlKHQ2Q2TG6IiO5CoxFwIllcBt423AtXMwsR+/5GPDZ/T4PacPJEci4uphfAyUHO7RbIpkma3MycOROdO3eGh4cHAgMDMWTIEJw/f/6u5/36669o0aIFnJ2dERsbi7///tsC0RKRvbqeVYT8knI4OcjRPMgDjXxd4eKoQGGZGhfTG04xv9+OJAMABrQOhqcze27Idkma3Pz777+YOHEi9u3bh02bNkGlUuG+++5DYWFhjefs2bMHI0aMwPjx43H06FEMGTIEQ4YMwalTpywYORHZE+2QVKtQTzgq5FDIZWgXIQ5PHbmeI11gRigtV+PP4zcBAI9yuwWycZImN+vXr8fTTz+N1q1bo127dli8eDESExNx+PDhGs/54osvcP/992Pq1Klo2bIlPvzwQ3Ts2BFfffWVBSMnInuiHZJqF+6tO9YhwgdAwynmt/VsOnKKVAj2dEaPaH+pwyEyK6uac5ObK/4A8fX1rbHN3r170b9/f71jAwYMwN69e6ttX1pairy8PL0PIiJjXL8l9ibHht2uC9OxsTcA4GgDSW60Q1JDOoQ12No8RIaymjo3Go0GkyZNQo8ePdCmTZsa26WmpiIoKEjvWFBQEFJTq1+eOXPmTMyYMcOksRKRffludCek5pXo7cHUvqLn5nJGIXKKyuDtqpQqPINMf6g1Wod64eH2oVKHQmR2VtNzM3HiRJw6dQorVqww6XWnTZuG3Nxc3UdSUpJJr09Etk8mkyHEy0VvEq6vmxJR/m4AgKNJORJFZrhGfq547d7maBrgLnUoRGZnFT03L730EtauXYsdO3YgPLz2iW7BwcFIS0vTO5aWlobg4OBq2zs5OcHJyclksRIRaY27Jwrlag2aB3lIHQoRVSJpz40gCHjppZewZs0abN26FVFRUXc9Jz4+Hlu2bNE7tmnTJsTHx5srTCKyY7M3nsczSw5i18XMKq891a0xxvaIQpi3iwSRGebMzTw8s+QQNp9Ju3tjIhshaXIzceJE/Pjjj1i+fDk8PDyQmpqK1NRUFBcX69qMHj0a06ZN0z1/9dVXsX79esyePRvnzp3D+++/j0OHDuGll16S4p9ARDZu35Vb2Hw2HbcKS6UOpU5WHU7G5rNpWH00WepQiCxG0uTmm2++QW5uLhISEhASEqL7+OWXX3RtEhMTkZKSonvevXt3LF++HAsWLEC7du2watUq/P7777VOQiYiqquMfDGpCfGqvnfm1I1c/HHsBnKLVZYMyyAqtQZ/HLsBAHikI2vbkP2QdM6NINy9bPn27durHHvsscfw2GOPmSEiIiJ9twrKAAB+7tWvhpq4/Aiu3yrCzxO6Ib6pnyVDu6vt5zNwq7AM/u5O6NU8QOpwiCzGalZLERFZmxKVGvml5QAAf7fqFyZoV0xdzay5srpUfjtcUdumfSgcFfxxT/aD73YiohrcKhR7bRwVMni6VN/RfTu5KbBYXIbILizDlnPiJOJHuN0C2RkmN0RENcismG/j5+YEmaz6qr5NKurGXMmwrp6bLefSoVILaBniiZYhnlKHQ2RRVlHnhojIGhWUlsNVqYC/R83Vh5tY6bCUq1KBzpE+6NbEuuYBEVkCkxsiohr0iPbHmQ/uR1m5psY2TQLE5CYxqwgqtcZq5rY8EBuCB2JDpA6DSBLW8V1IRGTFlA41/6gM8nCGi6MC5RoBydnFNbYjIsthzw0RUT3I5TJ8NKQNfNwcEeBhHVu9lJarodYIcFXyRzzZJ/bcEBHVYM6mCxi3+CC2n0+vtd0jceHo2yII7k7WkUzsvpSJ1u9twPjFB6UOhUgSTG6IiGpw5Ho2tp5LR1bFkvCG4mxKPgQBcHe2jmSLyNL4ziciqkFmgbgU3N+99uGm7MIy7L6ciVKVxipqypxNyQMAtAjmEnCyT0xuiIhqkHmXrRe0rmcV4aXlRxHo4WRVyU3LEA+JIyGSBoeliIiqodYIyKrYCTzgLj03UX7icvD0/FIUVGzXIJUSlVpXc4fF+8heMbkhIqpGTlEZNBV7+/q41d5z4+XqCL+KNtckLuZ3IS0fGgHwdVMi0EpWbxFZGpMbIqJqaIekvF0dDSrMpy3mdzlD2j2mzqXkAwBaBHvUuGUEka1jckNEVI2CUhXclIq7TibWspbdwSN8XfFoXDj6tQySNA4iKXFCMRFRNeIa++L0XbZeqCzKX9xAU+rkJr6pH+Kbcj8psm/suSEiqkVtWy9Upu25sbbdwYnsEXtuiIhMoFOkDxY8FYfoQHfJYigoLUdydhGaBrhbzQaeRFLgu5+IqBpfb7uEsYsOYNOZNIPa+7s74b7WwWgSIF1ys//KLdw/bycGf7VbshiIrAGTGyKiahxNzMa28xnIyC+VOhSDnUsVV0o1C5IuwSKyBhyWIiKqRkbFUnD/u1QnruzQtSzsv5qFuMY+6NbE8pN6z+gqE7N4H9k39twQEVXjVsW+Un4GLgUHgHUnUzBrw3lsPVf7LuLmcpbJDREAJjdERFUIgqDbNPNuWy9U1kTCFVPFZWpddeSWwdxTiuwbkxsiojsUlalRohLr2/h7GD4spZ1MfCXT8lWKtdsu+LkpEcBtF8jOMbkhIrqDttfGxVEBV6XhUxO1tW4SbxWhXG1Y8T9TqTwkxW0XyN5xQjER0R3yS8rh7uQAHzdHo84L9nSGs6McJSoNkrOLEVmR7FhC+0bemHJfc4R6u1jsnkTWiskNEdEd2oR54dSMAQZvvaAll8vQxN8dZ1LycCYlz6LJTYtgT7QI5kRiIoDDUkRENTJ064XKtEvAD1zNMnU4RGQg9twQEZnQmO6NMaxjGFqHWq4XJaeoDPuuZKF1qCcifF0tdl8ia8WeGyKiOyzYcRlPLzqAv0+mGH1uYz83tAnzsuik3sPXs/H8j4cxYekhi92TyJoxuSEiusPx5FxsP5+B1NwSqUMxiHalVAvWtyECwOSGiKiKzIr9pPzrWC/mRk4xJv9yDCMW7DNlWDU6W7GnFCsTE4k454aI6A63Civ2lXIzvIBfZa6OCvx+7AY0ApCcXYRwH/POg+G2C0T62HNDRHQHbRG/uvbc+Lgp0amxLwBgy1nz7jNVeduFFiEcliICmNwQEelRqTXIKVIBAPyN2FfqTv1aBgIANp9NM0lcNdFuu+DvrkSgh7NZ70XUUBid3Lz//vvQaKoWtsrNzcWIESNMEhQRkVSyKoakFHIZvF2Mq1BcWb+WQQCA/VeyUFBabpLYqnN7MjGHpIi0jE5ufvjhB9xzzz24cuWK7tj27dsRGxuLy5cvmzQ4IiJLyytWwcPJAb5uSsjldV/O3TTADZF+rihTa7DzQoYJI9SXEBOIL0d0wPh7osx2D6KGxujk5sSJEwgPD0f79u3x3XffYerUqbjvvvvw1FNPYc+ePeaIscEoUakx8vt9yC9RSR0KEdVRsyAPnJwxALvf7Fuv68hkMl3vzWYzzrsJ9nLGoHah6NMi0Gz3IGpojF4t5ePjg5UrV+Ltt9/Gc889BwcHB/zzzz/o16+fOeJrUP69kIHM/DJ4ON/uyt52Lh3Rge6sGkrUwNRl64U79WsZiN2XMtE8yN0EERGRoWSCIAjGnvTll1/irbfewpAhQ3D48GEoFAosX74c7dq1M0eMJpWXlwcvLy/k5ubC09O0Y9SZBaW4kV2MdhHeAMRVDB0/3IRilRoxQR7o1zIQ/VoGon2EDxT16O4mIgKAvBIVVh1KRqS/K/q2CJI6HCKzMub3t9E9N/fffz8OHTqEJUuW4NFHH0VxcTEmT56Mbt26YcaMGXjjjTfqHHhD5+/upLe6Ij2/BLHhXjh8PRvn0/JxPi0f/9t+Gb5uSiQ0D8CjceHoHu0vYcREdKfFu69i2/kMDOsYhsHtw6QOp1aX0wvwwdozCPZ0xr63mdwQaRnd76pWq3HixAk8+uijAAAXFxd88803WLVqFebOnWvyABuyxn5uWPlcPA6/0x9fPNEeg9qFwsPZAVmFZVh99AZO38zTtc0uLMOey5koLlNLGDERnbqZh38vZOBGTrHJrllUVo79V26Z7HpaiVlFAIBGfhz2JqrM6J6bTZs2VXv8wQcfxMmTJ+sdkC3ydlVicHvxr0CVWoOD17Kw40Kmrg4GINbCmLrqBBzkMrQJ80LnSB90ivRFu3BvBHk6WXQTPiJ7llssLgjwdqlbdeLqrtfl480oLdfg0Dv961U7506JtyqSG87pI9JTpxlzO3fuxKhRoxAfH48bN24AAJYtW4Zz586ZNDhb5KiQo3tTf7w1sAWaBNyeZFiiUiPI0wnlGgHHknLw3c6reG7ZYXSbuQXtP9iEk8m5urbJ2UW4mJaPojLz1c4gsld5FcmNp4tpdqfxcnFEdKD4vb7zommXhGt7bhozuSHSY/R372+//YannnoKI0eOxNGjR1FaKpYpz83NxSeffIK///7b5EHag6fiIzGqW2MkZxfj0PUsHLqWjUPXsnExPR+5xSqEeN+uPLpkzzV8t/MqAMDH1RFhPi4I9XJBmI8Lwrxd8GhcOLxdTfNXJ5G9ySsR/2jwdK57Ab879W4egNM387D9fAaGdgg32XU5LEVUPaOTm48++gjz58/H6NGjsWLFCt3xHj164KOPPjJpcPZGJpMhwtcVEb6uuh+AJSo1rmYWVunK9nB2QH5JObKLVMguUuHUjdvzdwa1C9U9nr3xPP4+mYIwH1eEeTvDy0UJD2cHeLo4wtPZAfe1CoaLUgEAKCgth1wGuDgqOAxGdut2z43pkpuEmED8b/tl7LiQAbVGMNlqSV1yw54bIj1GJzfnz59Hr169qhz38vJCTk6OKWKiSpwdFVV2+v2/B1vh/x5shbwSFW5kF+NmTjFu5BSLj3NLEFApEbqcUYDLGYW4nFFY7fWPvHuvLrn5zz9n8eO+RDjIZfB0cYSXiyP83ZXwc3OCv4cSk++NgW/FLskpucUoVWng7+EENyWTIbIduuTG2TTDUgDQoZE3PJwckF2kwskbuWhfUS6iPkpUaqTmlQBgckN0J6O/e4ODg3Hp0iVERkbqHd+1axeaNGliqrjIAJ7OjvAMcayS/FT2fw+2wsiujXEjuxgpuSXIK1Ehr1iF/JJy5JWo4FHpB3hBRXd8uUZAVmEZsgrLcDXzdlI0+d4Y3eNvtl/G0r3XAQDOjnLdMnh/dycEeCgx5b4Y+FUkWSeSc5CUVQylgxyuSgV83ZTwc1fC11UJBwX3biXrodEI0Bb+8jJhz42jQo57mvnjn1Op2H4+3STJjYNchjUv9kBiVpHujw4iEhmd3EyYMAGvvvoqFi5cCJlMhps3b2Lv3r2YMmUK3n33XXPESPUQ5i3OwzHE3OHt8fHQWF3ik11YhsyCMmQWlCKzoFRvE0FBEIevilVqlKg0SM4uRnL27aWzbwxooXv884Ek/Hwgscr9ZDLA28URf750j66C87oTKTh0PQsezuKwmbuTg/jYxQE+rko0C3KHk4Oirl8OolrJ5TKcmjEA5WqNyQtt9m4egH9OpeLfCxmY1L95va/noJCjfYS3SRIlIltjdHLz1ltvQaPRoF+/figqKkKvXr3g5OSEKVOm4OWXXzZHjGQhMpkMbk4OcHNyQLCXc61tPxzSBh8OaYPC0nJd8qNLhPLL9P7qjfRzRZcoX5SVa1BUVq7rFdIIQHaRCl6ut9vuvpyJ5furJkJa26ckINLfDQCwcNdV/H0ypaInyKliCE0JX3cn+Lkp0bGRj27IbfelTBy5no2UvBLkFaugVMjhqJBD6SB+frlvNHz41y9VMEePYt8WgXj3oVZIiAkw+bWJSF+dtl8AgLKyMly6dAkFBQVo1aoV3N0bxt4p5tx+gQyn1gjIKRJ7hpoHuevm7Gw8nYojiTkoKBWHzgpKynU9SbcKy7Dl9d66VSzTVp/AzweSarzHzjf66HqE3l5zstakafdbfXU9XHsv30JOURmaBXkg3McFzo7sKSLrs/5UKlJyi3FPtD+aBXlIHQ6R2Zl1+wUtpVKJVq1a1fV0snMKuQx+7k66eTla97UOxn2tgw26xtgeUejVLACZhWXIKijDrcJS3KroPcoqLNNb7dI1yhflag1CvFzg5eKIco0GKrWA0nINVGqN3uTR+f9exr8XbtcjCfBwQqi3CxzlMrgoFVg2vqvutbfXnMTRxBzIAET4uqBNqBfahHkhJtgDbk4Oej1Yp2/mIiO/FAWlYtJWWKZGiUoNXzclmgW6o0Mj7jkmtZPJufh843m0CPHAtIEtpQ6nVqsOJ2Pz2TR8OKQNkxuiOxiU3AwbNszgC65evbrOwRAZo3mQB5ob+ENdWyHaEC1CPJCeX4qkrCIUlJYjI78UGfliPScPJ/1vmcRbRTibIi7DP5OShw2n03Sv+bopceTde3XP315zCseTcqq9p5ODHGc+uF/3fPn+RJxNyYNMBsgrerVclQoEeTojyNMZ97cxLAEk49zIKcK/FzJQWGqeApklKjX+PHYTB69l4dNH2kJej2Q2MUuc7M8CfkRVGZTceHl56R4LgoA1a9bAy8sLnTp1AgAcPnwYOTk5RiVBRNZq2sCWmDZQfK/nFquQlFWMm7nFEAQBSgf9uRhv3B+DZ4uaQC0IuJxegFM3cnHqZh6uZBTAUaH/i6tpgBtU5Rp4ODvAw9kBrkoHODvKkZpXCrkMer02y/Zd1yVNd/JzU+olN5+tP4ecYpVucmm4jwtclaZbxmxP8oorCviZcKVUZXKZDDP+Oo3CMjXGdI9EmzCvu59UDUEQWOOGqBYG/QRctGiR7vGbb76Jxx9/HPPnz4dCIc5FUKvVePHFFzmHhWyKTCaDt6sS3q5KxIZX/0uobbi37nGfmNt7hQmCALVGfzrbnMfbG3zvh9uF4t5WQYAgQCMAAgQUlJQjPb8ULnfMAfrj2E3cyCnWm1Pk4eSAAE8ntArxxFdPdtQd33UxEw4KGQI9nBDk6Qw3JyZBleWaocZNZUoHObpH+2PTmTRsP59e5+QmI78UJSoN5DIg1MDVkET2xOjv4IULF2LXrl26xAYAFAoFJk+ejO7du2PWrFkGX2vHjh2YNWsWDh8+jJSUFKxZswZDhgypsf327dvRp0+fKsdTUlIQHMxuerIeMpkMDoq6Dzm8kNDUoHaCIOCdB1viWFIOjibl4PSNXBSWqZFfWo78jPIqw2hvrzmp+4sfANydHBDo4YRATye0CPbE+w+31r2WlFUEP3elXfUC5ZWYvjrxnRJiAiqSmwy81LdZna6h/T8M9Xap0ptIRHVIbsrLy3Hu3DnExMToHT937hw0Go1R1yosLES7du0wbtw4o4a0zp8/r9dLFBgYWEtrItslk8kwMDYEA2NDAIjJTkGp2MOTnleKOwtHNwlwg4NchnTtxOaKjyuZhShR6X//PvXDfly7VYQQL2e0DvVEXGNfdIr0QWyYl82uINNWJzZlAb879W4uLgU/kpiNwtLyOvWeXedu4ES1Mvq7auzYsRg/fjwuX76MLl26AAD279+P//znPxg7dqxR1xo4cCAGDhxobAgIDAyEt7e3QW1LS0t1m3sC4lIyIlslk8ng4ewID2dHNA2oWp5h8dguuseFFUlQWl4J0vNL4VypB0CtEVBQqgYApOSWICW3BJvPpgMAlAo5HogNxrwnOpj5X2N55tg0807hPq4I9HBCen4pzqXmI66xj9HX4HwbotoZndx8/vnnCA4OxuzZs5GSkgIACAkJwdSpU/H666+bPMDqtG/fHqWlpWjTpg3ef/999OjRo8a2M2fOxIwZMywSF1FD4ubkgCgnB0RVFEWsTCGX4dA7/ZFdWIYrmQU4mpgj7lR/PRuZBaV6QyFqjYDXVx5D96b+6NcysMry/oakRCUmdJ4u5h2KaxXqifTzGTiTklen5GZcjyj0jgmAO+dMEVWrzkX8gNu9IKaYSCyTye465+b8+fPYvn07OnXqhNLSUnz//fdYtmwZ9u/fj44dO1Z7TnU9NxERESziR1QHgiDohkS0laIPXcvCo/P3AgDkMqBTpC8GtA7GgNZBCPdpeD0L5WoNBIj7QZnLZ+vP4Zt/L2NiQjSmDIi5+wlEZFQRv3olN6ZkSHJTnd69e6NRo0ZYtmyZQe1ZoZjItG7mFGPV4WRsOJ2K0zf1h30fjA3BG/fHoLFf1d4he5ZTVAZHhZyr1YiMYMzvb6P/NElLS8NTTz2F0NBQODg4QKFQ6H1YWpcuXXDp0iWL35eIRKHeLnilXzOse6Undr7RB+8+1Apdo3whkwHrTqYgKav47hexM96uyjonNsVlany49gyW7r1WpdwAEYmM/u56+umnkZiYiHfffRchISG6PYGkcuzYMYSEhEgaAxGJInxdMf6eKIy/JwrnUvPw98lU3NPMX/f66Zu5iAnyMMvGlKbw7NJDcHZU4P2HW8PXSjdSTcwqwg+7rsLT2QGj4yOlDofIKhmd3OzatQs7d+5E+/bt633zgoICvV6Xq1ev4tixY/D19UWjRo0wbdo03LhxA0uXLgUAzJs3D1FRUWjdujVKSkrw/fffY+vWrdi4cWO9YyEi02oR7IkWwbe7jtPzS/DEt/sQ5uOCDwa3QZcoXwmjq0ql1mDjGXH7jBmV6v2Yyw+7rmL9qRSMv6eJUdtpaFdKcaiPqGZG//kUEREBU03TOXToEDp06IAOHcQlpZMnT0aHDh0wffp0AGJxvsTE21VXy8rK8PrrryM2Nha9e/fG8ePHsXnzZvTr188k8RCR+VxKL4BCIcO51Hw8/u1eTFt9EsVlaqnD0tHWuAEADzNVKK7sSkYBDl7LxrEa9huryfVb4p5SXAZOVDOjv4PnzZuHt956C99++y0iIyPrdfOEhIRaE6XFixfrPX/jjTfwxhtv1OueRCSN7k39se31BHy24TxWHEzEzwcScfh6Fr56sqPBG6Cak7bGjbuTg0WGzVqFir1aNe0hVpMkbY0bPyY3RDUx+jt4+PDh2L59O5o2bQoPDw/4+vrqfRAR1cTHTYmZw2KxbFxX+Ls74UJaAR7+ahd+OZh495PNLM/M+0rdqVWImNycMTK5YQE/orurU88NEVF93NPMH/+82hOTVx7DzouZOHUjD8M7SxuTJfaVqiwm2AMymbgJZnp+CQI9nA067zqTG6K7Mjq5GTNmjDniICI7E+DhhCVju+CXQ0kY2iFM6nCQV1yx9YKFkhtXpVgd+kpGIc6m5BuU3Gg0ApIrltYzuSGqmUHJTV5enq5gzt32ZmJhPCIylFwuw4gujXTP1RoBX2+7hDHdI826eWV1CkvNv6/UnVqFeFYkN3m6DTVrI5MBu9/qi8SsQoR6u1ggQqKGyaDkxsfHBykpKboNK6urbSMIAmQyGdRq61n9QEQNywd/ncaSvdex5Vw6lo3vYtFE4/HOERjaMQwqtebujU2kVagnDl7LgqELUGUyGQI8nBDg0XD37yKyBIOSm61bt+omC2/bts2sARGR/XqiSyP8cfwmjiflYMzCA1g6rgs8LJjgOCrkZt1T6k7P9WqKFxOiLXY/InthNXtLWQr3liKybqdv5uLJ7/Yjt1iFuMY+WDKuC3e/rvDF5osoKFVheOcIRAdKv3yeyJLMurcUEZE5tQ71wk/PdIWnswMOX8/GuEUHUaIy/3D3/H8v45Wfj2LPpUyz36s6GgP2ifr1cBK+23kVaXmlFoiIqOFickNEVqdNmBd+fKYrPJwdcOBaFv7zzzmz33PflVv48/hN3Mix7Eafn284j84fb8aKg0m1trtVUIrkbDG22HAvS4RG1GAxuSEiq9Q23BtfjugAT2cHxDX2Mfv9costW+dGS6XWICO/FGdScmttd+KG+HqTADeLTrQmaog4kE1EVishJhA73+xrkWXh2grFll6Cfnsbhvxa251IEpObduHe5g6JqMGrU89NeXk5Nm/ejG+//Rb5+eI35M2bN1FQUGDS4IiIKicbtwpKoTZgbkpdaPeWsnSviHYbhrMpebXOuzmenAMAaMshKaK7Mjq5uX79OmJjYzF48GBMnDgRGRkZAIBPP/0UU6ZMMXmAREQAsPNiBgbM24n/bbtkluvr9pZysWyHdpS/G5QOchSVqXVbK9xJEAScqEhu2kV4Wy44ogbK6OTm1VdfRadOnZCdnQ0Xl9sVMocOHYotW7aYNDgiIq2M/FJkFpTiy22XkJxdfRJQVyUqNUrLxeJ9lp5z46CQo0WwuKz7zM3qK8DnFKkgCICDXKbr6SGimhmd3OzcuRPvvPMOlEql3vHIyEjcuHHDZIEREVU2tEMYujXxRVm5BrM2nDfptbWbZsplgLvS8lMRb+8QXv2kYh83JQ690x973uoLZ0eFJUMjapCMTm40Gk21WywkJyfDw4NFpYjIPGQyGd55sBVkMuCPYzdxNDHbZNcO9HDGxY8H4vA790Iur7q9jLl1bOyDrlG+CPepeTNMmUyGQE/Ddg4nsndGJzf33Xcf5s2bp3suk8lQUFCA9957Dw888IApYyMi0tMmzAuPdAwHAMz827S1bxwVcvi4Ke/e0Awe7xSBX56L19tElIjqzujkZvbs2di9ezdatWqFkpISPPnkk7ohqU8//dQcMRIR6Uy5LwZKhRwHrmXpJtnaMo1GQL/Z2zFh6SFkFZZJHQ5Rg2D04HJ4eDiOHz+OFStW4MSJEygoKMD48eMxcuRIvQnGRETmEOzljAfbhmDN0RvYeTETbU1Q92XflVtYvj8R7SO8Me6eqPoHWUcFpeUoKi3XG366klmAyxmFuJFTDE9nliYjMoTR3yklJSVwdnbGqFGjzBEPEdFdvdKvGZ7t1QQtTbRy6GJ6Af48fhOl5WrJkpsvNl/E3M0X0DzIHWtf7gmlg9ixfryieF9smBccLLhjOVFDZvR3SmBgIMaMGYNNmzZBo9GYIyYiolpF+buZLLEBKtW4kXBbg9HxjeHnpsSFtALM//ey7vgJXfE+b2kCI2qAjE5ulixZgqKiIgwePBhhYWGYNGkSDh06ZI7YiIjuKrOgFCp1/f7Q0i4Ft3SNm8p83JSYPqgVAOCrrZdwKV2s/n4sWey5YWViIsMZndwMHToUv/76K9LS0vDJJ5/gzJkz6NatG5o3b44PPvjAHDESEVXr43Vn0H3mVmw4nVqv6+QVi1svWHpfqTs93C4UfWICUKbWYNrqkygtV+NsRWG/9qxMTGSwOg/genh4YOzYsdi4cSNOnDgBNzc3zJgxw5SxERHVysVRgTK1Bkv2XKvXdXKKxFVIUic3MpkMHw2NhatSgYPXsvHeH6dRptbA29URjXxrroFDRPrqnNyUlJRg5cqVGDJkCDp27IisrCxMnTrVlLEREdVqZLfGcJDLcPBaNk7frL66ryGyK5Ibb1dpkxsACPN2wRsDYgAAKw4mIczbBZ0jfSGTWb64IFFDZfRqqQ0bNmD58uX4/fff4eDggEcffRQbN25Er169zBEfEVGNgjydcX+bYKw9kYJfDyWj9cN1m5eSUyTOufFxlaaI352eio/ExjNpSIgJwLgeUVBIUDWZqCEzOrkZOnQoHnroISxduhQPPPAAHB2l/0uHiOyXNrk5Uo/tGNa+fA9yi1Vwc7KOOjIKuQw/PdOVvTVEdWT0d3JaWhr3kCIiq9GuYon02ZQ8lJar4eRg/MaSDgo5/NydTBxZ/TCxIao7g5KbvLw8eHqKNSUEQUBeXl6NbbXtiIgsIdzHBb5uSmQVluFsSj5XFRGRYcmNj48PUlJSEBgYCG9v72r/ohAEATKZrNodw4mIzEUmk+H53k3g5KBAqJfxu2Zn5Jfig7VnEOjhhHcfamWGCInI0gxKbrZu3QpfX18AwLZt28waEBGRsZ7t1bTO56blleCv4zcRwOSGyGYYlNz07t1b9zgqKgoRERFVem8EQUBSUpJpoyMiMrPbK6W4OILIVhhd5yYqKgoZGRlVjmdlZSEqSrrddInIvl1Iy8fKQ0koKC036rzbNW6sYxk4EdWf0aultHNr7lRQUABnZ+PHu4mITGHc4oNIzi5GuLcLukf7G3yetjoxe26IbIfByc3kyZMBiJP33n33Xbi63i4FrlarsX//frRv397kARIRGaJduDeSs4txPDnXqOQm28oK+BFR/Rmc3Bw9ehSA2HNz8uRJKJW3fxAolUq0a9cOU6ZMMX2EREQGaBfhhXUnU3A8Kceo8zgsRWR7DE5utKukxo4diy+++IL1bIjIqrStKOZ3PDnHqPNyOaGYyOYYPedm3rx5KC+vOmEvKysLDg4OTHqISBKxYV6Qy4CU3BKk55Ug0NOwOYCzHmuHdx5qBQcFKwIT2QqjV0s98cQTWLFiRZXjK1euxBNPPGGSoIiIjOXm5IDoQHcAwPFkw3cIV8hl8HVTwtOZPTdEtsLo5Gb//v3o06dPleMJCQnYv3+/SYIiIqoL7T5TJ4wcmiIi22L0sFRpaWm1w1IqlQrFxcUmCYqIqC6eim+MB9qGoH1FkmOIN1Ydh5ODAq/d2xy+bpxUTGQLjO656dKlCxYsWFDl+Pz58xEXF2eSoIiI6qJtuDf6xATCx8AkRa0R8OvhZCzbdx0aQTBzdERkKUb33Hz00Ufo378/jh8/jn79+gEAtmzZgoMHD2Ljxo0mD5CIyFxyi1XQ5jTeLpxzQ2QrjO656dGjB/bu3YuIiAisXLkSf/31F6Kjo3HixAn07NnTHDESERls2/l0zNl4Hqdu3H1SsbbGjYezAxwURv84JCIrZXTPDQC0b98eP/30k6ljISKqt18PJeHvk6nwdHFEmzCvWtve3nqBc22IbEmdkhutkpISlJWV6R1jnRsiklLTAHE5+OWMwru2zS5kAT8iW2R0P2xRURFeeuklBAYGws3NDT4+PnofRERSup3cFNy1LbdeILJNRic3U6dOxdatW/HNN9/AyckJ33//PWbMmIHQ0FAsXbrUHDESERlMm9xcMSC5yS1mzw2RLTJ6WOqvv/7C0qVLkZCQgLFjx6Jnz56Ijo5G48aN8dNPP2HkyJHmiJOIyCBNAtwAAJkFZcgtUsGrlsRlbI8oPNIxXH8ZuFoFZF8D/KIBGbdkIGqIjO65ycrKQpMmTQCI82uysrIAAPfccw927Nhh2uiIiIzk5uSA4Ip9pS5n1t57o5DL4OOmhJ+7k3jg6g7gm+7AV52AU7+ZO1QiMhOjk5smTZrg6tWrAIAWLVpg5cqVAMQeHW9vb5MGR0RUF00Dxd6bS+l3H5oCAJTkAr9NAJYMAjIviMeSDpgpOiIyN6OHpcaOHYvjx4+jd+/eeOuttzBo0CB89dVXUKlUmDNnjjliJCIyynuDWkOpkCPcx6XWdnM3XUBWYRleE5bC9+RKADIgpC2QchzIumKZYInI5IxObl577TXd4/79++PcuXM4fPgwoqOj0bZtW5MGR0RUF82DPAxq98vBJKTmlWBqkzPigQdmAQEtgCUPAVmXzRihCZ1cBRSkAV2fB+QKqaMhsgpGJTcqlQr3338/5s+fj2bNmgEAGjdujMaNG5slOCIicylRqZGaVwIAcCtJFQ/6Nwf8moqPs6+Lk4sVVrySqvAWsPpZQFAD13YDj3wPKF2ljopIckbNuXF0dMSJEydMdvMdO3Zg0KBBCA0NhUwmw++//37Xc7Zv346OHTvCyckJ0dHRWLx4scniISLbUFquxpdbLuK1X46hrFxTbZukrCIAgIeTAvL8m+JBzzDAIwRwcBEThpxES4VcN+f/FuMEgPPrxB6nggxpYyKyAkZPKB41ahR++OEHk9y8sLAQ7dq1w9dff21Q+6tXr+LBBx9Enz59cOzYMUyaNAnPPPMMNmzYYJJ4iMg2KBVyfPPvZaw5egOJFUnMna7fEo+38AFkqopqxp4h4vJvX3FFqNXPuzn7l/i55cOAiw9w4zDwQ3/gVgMZUiMyE6Pn3JSXl2PhwoXYvHkz4uLi4Obmpve6MZOKBw4ciIEDBxrcfv78+YiKisLs2bMBAC1btsSuXbswd+5cDBgwwODrEJFtk8lkaBrgjpM3cnE5owDRge5V2lyvSHraehYAOQCcvQFlxc8zvyZA+mkxSWh2r8XiNkpJHnBlm/i4z/+J821+fESs0fN9f+DJX4CILpKGSCQVo5ObU6dOoWPHjgCACxcu6L0mM3PBq71796J///56xwYMGIBJkybVeE5paSlKS0t1z/Py8swVHhFZkSYBbjh5IxdXathjKvGWeLy5S8Vycc+w2y/6Vsy7seZJxRc3AuoywK8ZEBAj9jg9sxlYPhy4eURc1j7sO6DVw1JHSmRxBiU3J06cQJs2bSCXy7Ft2zZzx1Sj1NRUBAUF6R0LCgpCXl4eiouL4eJSddnnzJkzMWPGDEuFSERW4m57TKXliX/0RCqzxQNelZObBjAspRuSGnS7krJ7IPD0WmDVeODCP8DK0cDDXwIdn5IuTiIJGDTnpkOHDsjMzAQgFvG7deuWWYMypWnTpiE3N1f3kZSUJHVIRGQBd0tuvhnVEUffvRftPbXzbUJvv6hdMWWtc1dUxcDFTeLjloP0X1O6AU/8BHQaD0AA/p4CpJ+zeIhEUjIoufH29tZVJb527Ro0mupXH5hbcHAw0tLS9I6lpaXB09Oz2l4bAHBycoKnp6feBxHZPm2V4svpBRAq7x1VQSYTt15wKqpYBl7dsFROorgc3Npc3gqoCgHPcCC0Q9XX5QrgwdlAdH+gvARYPQEoL7N8nEQSMWhY6pFHHkHv3r0REhICmUyGTp06QaGovljUlSvm68aNj4/H33//rXds06ZNiI+PN9s9iahhivRzg0wGFJapkVOkgo+bsvqGedpl4JV6bjyCAUdXQFUkJjjanhxrUd2Q1J1kMmDw18D/4oHUE8C//wH6TbdcjEQSMii5WbBgAYYNG4ZLly7hlVdewYQJE+DhYVgF0NoUFBTg0qVLuudXr17FsWPH4Ovri0aNGmHatGm4ceMGli5dCgB4/vnn8dVXX+GNN97AuHHjsHXrVqxcuRLr1q2rdyxEZFucHRXY+noCwrxdoHTQ76Q+mpiNL7deQlxjH0ysLrnRLgdPOyUOTVlTcqNWifVtgKpDUnfyCAYGfQGsfArYNRdodh/QqJv5YySSmMGrpe6//34AwOHDh/Hqq6+aJLk5dOgQ+vTpo3s+efJkAMCYMWOwePFipKSkIDHxdhGtqKgorFu3Dq+99hq++OILhIeH4/vvv+cycCKqVpS/W7XHz6XmY+u5dKg1Aibm3RAPeobrN9ImN9a2YuraTnGjT1d/wxKVVg8D7Z4Eji8H1jwHPL8LcKr/z28ia2b0UvBFixYBAC5duoTLly+jV69ecHFxgSAIRi8FT0hIqHYsXKu66sMJCQk4evSoUfchIqpMW8CvubcAJFaUh/AM0W9krSumtENSLR40fC+pgf8Bru0Sa+CsnwYM/sps4RFZA6MrFGdlZaFfv35o3rw5HnjgAaSkpAAAxo8fj9dff93kARIR1dXpm7mY8utxvP/nab3jiVniCqkY13zxgJNX1d4Ma1wxpVEDZ9eKj1saUb/G2QsY+g0AGXB0GXCOQ/lk24xObiZNmgRHR0ckJibC1fX2Bm3Dhw/H+vXrTRocEVF9KBVyrDqcjGX7riO9YpNM4HbPTZRS22sTWvVkayzkd/MYUJgOOHkCUb2MOzfyHqD7y+LjP18BCtJNHh6RtTA6udm4cSM+/fRThIfrj083a9YM169fN1lgRET11SzIA3GNfaDWCPj1cDIAQBAEJFYkN2GKLLFhtclNxbBUTqL1LKPOrKgKH9oecKhh9Vdt+r4DBLUBijKBP18GapkWQNSQGZ3cFBYW6vXYaGVlZcHJyckkQRERmcqILo0AACsOJkKjEZBdpEJ+aTkAwE8tFifVq06s5REMOLoBgsZ6dgfPqfgD0rtx3c53cAKGLQAUSuDCeuDIEtPFRmRFjE5uevbsqVuaDYiFsDQaDT777DO9lU9ERNbgwdgQeDg7ICmrGLsvZyI9vwRuSgWCPJ3gWKBdBl5NcqO3O7iVDE1lVyQ3PnVMbgAgqPXtejfr37auOUVEJmJ0cvPZZ59hwYIFGDhwIMrKyvDGG2+gTZs22LFjBz799FNzxEhEVGcuSgWGdRCTl58PJKJFsCdOzRiATZN73y7g5xFS/cm+UeJna0kAdD03kfW7TreJQGRPscrxmucAdXm9QyOyJkYnN23atMGFCxdwzz33YPDgwSgsLMSwYcNw9OhRNG1qRYWuiIgqjOgqDk1dzSxCVmEZZDIZPJ0dgcIMsYF7UPUnaldMWctycFP03ACAXA4M+Z84MTn5oFjgj8iGGFXnRqVS4f7778f8+fPxf//3f+aKiYjIpFoEe6JTYx+E+bjAx9Xx9gtFFROK3fyrP9GaVkyVlwHagoN1nXNTmXcj4IHPgTXPilszRPcFwuLqf10iK2BUz42joyNOnDhhrliIiMzmqyc7YmCbYGgqLxAqqphQ7Opb/UnaOTfWMCyVmwRAABxcAPdA01yz7eNAqyGAphxYPhxIP2ua6xJJzOhhqVGjRuGHH34wRyxERGYT7OWM+9uEQCGvqKReViRujAmIWxlURzsslZsk/XJw3XybRjVvlmksmQwYNA8IjhWH6BY/BKSdvutpRNbO6O0XysvLsXDhQmzevBlxcXFwc9Pfu2XOnDkmC46IyGyKbomf5Y4177XkHiQuB1cVismFfzPLxXcnU823uZOLDzD6T2DZECDluJjgjPlTTHiIGiijk5tTp06hY8eOAIALFy7ovWbs3lJERJLRJjdu/jX3hOh2Bz8pTiqWMrmpb42b2rj6AqP/AJYNBW4eBZYMEp+HtDP9vYgswOjkZtu2beaIg4jIsnTzbfxqb+cVJiY3+Snmj6k25uq50XLxAZ76HfjxEeDGodsJTmgH89yPyIyMnnNDRGQTtCul7pbcaFdSaZeNS8WcPTdaLt7AU6uB8C5ASS6wZDBw47D57kdkJkxuiMg+FRrYc+NWsTKpQOLkxtw9N1rOXmKCE9ENKM0Flg4Bkg6a955EJsbkhojsU+U5N7VxCxA/S9lzU1pwexjNnD03Wk4ewKhVQKPuQGmeOBcncb/570tkIkxuiMg+aZObu/XcaGvKSJncaDfudPYSh44sQZvgRPYEyvKBH4cB1/da5t5E9cTkhojsk6ETiq1hzo0l5ttUR+kGPLkSiOoFlBWIk42v7bZsDER1wOSGiOyTwROKtXNu0s0bT20sNd+mOkpXYMQvQJM+Yr2fnx4Fru60fBxERmByQ0T2yeAJxRVzboqzpNs9W6qeGy2lKzDiZ6BpP7Gq80+PsQeHrBqTGyKyT4ZOKHb1BWQVPyq1Q1mWpp1z4xMpzf0BwNEFeGI50Ow+oLwY+ONFQFUiXTxEtWByQ0T2R6MRe2KAu/fcyBW320g17yZb4p4bLUdn4NFFgEcIkH0N2PultPEQ1YDJDRHZn5IcQNCIj++W3ADSzrsRhNvDUlLMubmTkztw74fi451zgNxkaeMhqgaTGyKyP9r5Nk5egMLx7u11K6YkGJYqzhZrzQDijuDWIPZRscifqgjYNF3qaIiqYHJDRPZHN9/GgF4boFKtGwl6brS9Nu5B4rwXayCTAQ98BkAGnPqNk4vJ6jC5ISL7Y2gBPy0pqxRby3ybO4W0A+KeFh//84Z0K8mIqsHkhojsj66A311WSmlpkxsp9peypvk2d+r7rlg1Oe0UcGSx1NEQ6TC5ISL7w54b03DzA/q8Iz7e+tHtwohEEmNyQ0T2p1Cb3Pga1t4a5txYY88NAHQaBwS2Eic+b/tY6miIADC5ISJ7ZGgBPy0pV0tZc88NACgcgIGfio8PLQRST0obDxGY3BCRPTJ000wtt0o7gwuCeWKqjkZTqTqxlSY3gLixZqshYu2gf9607NeIqBpMbojI/ujm3Bg5oVhdBpTkmiem6hSkAepSQKYAPMMtd9+6uO9DwMEFuL4bOL1a6mjIzjG5ISL7U2jkhGJHZ8DJs+JcC04q1s638QoTh3+smXcj4J7XxMcb3wXKCqWNh+wakxsisj/GFvEDKs27sWByY+3zbe7U4xXAqxGQdwPYNVfqaMiOMbkhIvuiKgZUFb0KhvbcANLsL2XtK6Xu5OgCDKhYMbX7v0DWVWnjIbvF5IaI7Iu210bueHuoyRCS9txEWu6e9dVyEBDVW5wrtPEdqaMhO8XkhojsS+UCfjKZ4efpat1YcDl4Q+u5AcSv6cDPxEnQ59YCl7ZIHRHZISY3RGRfCo1cBq6lq1JsyWGpimXg1rIbuKECWwBdnxMfr38LUKukjYfsDpMbIrIv2i0CjJlMDFh+CwaNGsi7KT72svJl4NXp/aa41D7zAnBggdTRkJ1hckNE9sXYAn5alt48syAd0KjE4R33YMvc05RcvIH+74mPt//HshOxye4xuSEi+2JsAT8t90pVii0h74b42SPE+mvc1KT9KCC0A1CaB2yeIXU0ZEeY3BCRfan3nBsLJTe5SeLnhjgkpSWXi5OLAeDYj8C13dLGQ3aDyQ0R2RdjN83U0iY3pXmAqsS0MVUnt6LnxivM/Pcyp4guQMcx4uM/XxbrDBGZGZMbIrIv2gnFrr7GnefsBSiU4mNL9N7kJoufG3LPjda9H4jzhrIuA/9+KnU0ZAeY3BCRfdFNKDay50Yms+zQVF5FcmPtG2YawsUbeGiO+Hj3f4GU45KGQ7aPyQ0R2ZciIzfNrMySVYptqecGAFo8CLQaAghq4I+XAHW51BGRDWNyQ0T2Q6OpVOfGyJ4b4Pb+UhZJbmxkzk1lD8wCnL2B1BPA3i+ljoZsGJMbIrIfJTlizwEAuBg55waoVOvGzDVbVCW3KyF7RZj3XpbkHgjcP1N8vP0/QOYlaeMhm8Xkhojsh3ZIyskTcFAaf767ds6NmfeX0ta4cXQFXHzMey9LazcCaNoXKC8B/npF7E0jMjEmN0RkP+oz3waw3P5S2uTGM8y4zT0bApkMeGge4OgGXN8NHFksdURkg5jcEJH9qGsBPy1LzbmxtcnEd/JpDPR7V3y8cfrt+UVEJsLkhojsR10L+GlpzzP3/lK2OJn4Tl2eBcI7A2X5wLrJgCBIHRHZECY3RGQ/6jssZan9pXRbL9jQZOI7yRXAw18Cckfgwnrg1G9SR0Q2hMkNEdkPU825KcoENGrTxFSdynNubFlgS6DXVPHxP28ChbekjYdsBpMbIrIf9U1uXP0ByABBc/ta5mDrc24qu+c1ILCVmDBumCZ1NGQjrCK5+frrrxEZGQlnZ2d07doVBw4cqLHt4sWLIZPJ9D6cnZ0tGC0RNVjaCcV1nXOjcKg078ZMK6YEwb6SGwcl8PBXgEwOnPgFuLhJ6ojIBkie3Pzyyy+YPHky3nvvPRw5cgTt2rXDgAEDkJ5e8w8OT09PpKSk6D6uX79uwYiJqMGqb88NcHvFVEFa/eOpTkkuUFYgPrb1YSmt8Dig6wvi478mAaX5koZDDZ/kyc2cOXMwYcIEjB07Fq1atcL8+fPh6uqKhQsX1niOTCZDcHCw7iMoKMiCERNRg1XXTTMr004qNlfPjbbXxsUXULqa5x7WqO//Ad6NxQ1Dt3wgdTTUwEma3JSVleHw4cPo37+/7phcLkf//v2xd+/eGs8rKChA48aNERERgcGDB+P06dM1ti0tLUVeXp7eBxHZKe2+Uq512HpBy73ijylzFfLTTia2hyGpypRuwKAvxMcHFgDHlksbDzVokiY3mZmZUKvVVXpegoKCkJqaWu05MTExWLhwIf744w/8+OOP0Gg06N69O5KTk6ttP3PmTHh5eek+IiJseGklEdVMVXJ7uKc+w1LuZt5fyh6WgdekaR8g/iXx8R8TgdO/SxoONVySD0sZKz4+HqNHj0b79u3Ru3dvrF69GgEBAfj222+rbT9t2jTk5ubqPpKSkiwcMRFZBe18G7kD4OxV9+toe27MNefGHgr41ea+j4AOT4kr0n57hhOMqU4kTW78/f2hUCiQlqb/QyItLQ3BwcEGXcPR0REdOnTApUvV7y7r5OQET09PvQ8iskOVJxPXZ78mXXJj5jk39jYspSWTicNTrYcBGhXwyyjg2i6po6IGRtLkRqlUIi4uDlu2bNEd02g02LJlC+Lj4w26hlqtxsmTJxESEmKuMInIFphiMjFwu5CfuZMbe1kpVR25Ahi2AGh+v7h7+PLhQPJhqaOiBkTyYanJkyfju+++w5IlS3D27Fm88MILKCwsxNixYwEAo0ePxrRptws7ffDBB9i4cSOuXLmCI0eOYNSoUbh+/TqeeeYZqf4JRNQQmGIyMWD+Yak8bc+NHc65qUzhCDy2BIjqJc6V+nEYkHpK6qiogXCQOoDhw4cjIyMD06dPR2pqKtq3b4/169frJhknJiZCLr+dg2VnZ2PChAlITU2Fj48P4uLisGfPHrRq1UqqfwIRNQT1LeCnpU1uirMAtUr8JWwqGjWQd1N8bK9zbipzdAae+BlYNhRIPgAsGwKMXQ/4R0sdGVk5mSDY11aseXl58PLyQm5uLuffENmTrR8DOz4DOj8DPDi77tfRaIAP/QFBDUw+C3iGmi7GvBRgTgtApgDeSRcrIhNQnAMseQhIPQl4hgPj/gG8G0kdFVmYMb+/JR+WIiKyCFPNuZHLKxXyM/HQlHa+jUcIE5vKXLyBUWsA/+bisN3SwUB+9eVCiAAmN0RkL0yx9YKWLrnJqP+1Ksuz85VStXEPAJ76XeyxyboCLB1yex4V0R2Y3BCRfSjUJjf1nFAMmG9SsW4ZOOfbVMsrDBj9p9izlXFWnGRcwqrzVBWTGyKyD9rtEtxNsBeduTbPzLXTrReM4RsFjP5D7IG7eVRcJl5WJHVUZGWY3BCRfdAmIqZIbrTDUoUmHpay560XjBEQAzy1BnDyAhL3AL+MBMpLpY6KrAiTGyKyfaoSoCRXfKxNTOrD3MNS9lzAz1Ah7YCRvwKOrsDlrcCqcYC6XOqoyEowuSEi26cdklIo67evlJa5Ns+01x3B66pRV+CJ5eL/67m14mabGo3UUZEVYHJDRLavoNJ8m/rsK6Vljv2lVCW3h7mY3BiuaR+xkrFMAZxYAfz9OmBf5duoGkxuiMj26ebbmGBICjBPcqPttXF0BVx8THdde9DiAXEvKsiAQwuBTdOZ4Ng5JjdEZPsKTLhSCri9eWZprtjjYgqV59uYonfJ3sQ+Ku4mDgB7/ivOwTF1HSJqMJjcEJHt0yU3Juq5cfYCFE7i40IT9d5wvk39xY0BBn4mDlGdXg183Rk4/gt7cewQkxsisn2mXAYOiD0rph6aYgE/0+j6HDBhCxAUCxRnA2ueBX56DMhJkjoysiAmN0Rk+7TJjXY4yRR0K6ZMtBycNW5MJ7QD8Ow2oO+74kqqS5uA/3UDDnzH1VR2gskNEdk+U8+5qXwtU/XcZF8TP7PGjWkoHIFeU4DndwMR3YCyAuDvKcDiB4DMi1JHR2bG5IaIbJ8pt17Q0m2eaYLkpiQXSNwnPg7vXP/r0W0BzYGx/wAPfA4o3YHEvcA3PYCdswG1SuroyEyY3BCRbRME008oBkxbpfjiJkBdBvg3BwJb1P96pE8uB7pMAF7cC0T3B9SlwJYPgO/6ADePSR0dmQGTGyKybWUFgKpiY0VTJjfa+TumWC119k/xc8tB9b8W1cy7ETByFTD0W7GWUOpJ4Lu+wOb3AVWx1NGRCTG5ISLbpu21UboDSjfTXddUc25UxWLPDQC0eKh+16K7k8mAdk8AEw8CrYcBghrYNVccqrq2W+royESY3BCRbTN1dWItUw1LXd4m9ix5hourfMgy3AOAxxaJe1N5hABZl8XJxmsnAyV5UkdH9cTkhohsm6lr3GjploLXswru2b/Ezy0HsTKxFFo8CLy4D+g4Rnx+6Adx2fiFDdLGRfXC5IaIbJs2+TB1z41bxfVUhUBpQd2uoVYB5/8WH3O+jXRcvIGH/wuM+QvwiRKrRS9/HPjtGaAwU+roqA6Y3BCRbTNXz42TuziPp/I9jHVtF1CSA7j6A426mSw0qqOoXsALe4DuLwMyOXDyV+DrLsDJVdzCoYFhckNEtk1XndjEPTfA7RVTdZ1UrB2SavEgIFeYJiaqH6UrcN9HwDObgcDWQNEt4LfxwPLhQOopqaMjAzG5ISLbZo4aN1ra3qC6LAfXaIBza8XHLR82XUxkGmFxwLPbgT7viFs4XNwAzO8BLBsGXNnOnhwrx+SGiGybuYalgPpVKU4+KMbm5CkOh5D1cVACvacCz+8CWg8Vh6oubwGWDgYW9BaHq9TlUkdJ1WByQ0S2rdBME4orX7Muc260hfua3y/+EiXrFRADPLYYePkI0OVZwMEFSDkuDld92QHYNx8oK5Q6SqqEyQ0R2S6NxjybZmrVtZCfIFQakuIqqQbDNwp4YBbw2mkg4W3A1Q/ISQTWvwnMaQVs/ch0G6lSvTC5ISLbVZIDaCo2R9RO/jWlug5LpZ0SdwF3cAGi+5k8LDIzNz8g4U0xyXlwDuDbRHyv7ZgFzG0D/PUqkHlJ6ijtGpMbIrJd2uEiFx/zDP241XFYSrtKKrqfabeEIMtydAE6jwdeOgQ8vgwI6yRuynl4MfBVJ2DFSCDpgNRR2iUmN0Rku8w5mbjydQvSDFs9IwhA8iHgxC/icw5J2Qa5Amj1sLh8fOw/QPOBACqGHn+4F/hhAHBunThMShbhIHUARERmk5cifjZXcuMdAcgUYkXbX0YBD38JuPpWbZebDBxfIX7cuigeU3oAzQeYJy6ShkwGNO4ufmScB/b8FzixEkjaB6zYB/g1EwsEth0OODpLHa1NY88NEdmu3CTxs3eEea7vHgg88BkgdxT/Sv+mO3DlX/G10gLg2M/AkofFeRhbPxQTGwcX8Zfb2HXicBnZpoAYYPDXwKSTwD2vAU5e4v//X68Ac1sBq58Fjv9S/73JqFoyQbCvSkR5eXnw8vJCbm4uPD09pQ6HiMzpj5eAo8uAPv8H9H7DfPdJOQ6sGl/RKyMDmvYFEveJ+05pNb4HaD8CaDUYcPIwXyxknUrzgSNLgb3/A/KS9V8LbivOv2raD4joytIANTDm9zeTGyKyXUseBq7+CwyZLyYW5lRWCKyfBhxZcvuYTxTQ/kmg7eOAT6R5708Ng1olJr6XtwCXtgCpJ/RfV7oDkT0rkp2+gF9TaeK0QkxuasHkhsiO/LcDkHUFeHodEHmPZe557m/g+m5xsnBEV3EeBlFNCtKBy9vEZOfy1ttFJ7V8IsUeneh+YiVrO+71Y3JTCyY3RHZCowE+DgLUZeK8B+9GUkdEVDuNBkg7KfboXN4q9vBo6zQBgNxBTJib9hWTneB2gNx+ps4yuakFkxsiO5GXAsxpIa5meicdUHBxKDUwpQXAtZ0Vyc4WsReyMld/oGkfsWenaV/Aw0yrAq2EMb+/+d1ORLZJu1LKM5SJDTVMTu5AzEDxAwCyrlbM1dkKXN0BFGUCJ38VPwAgKBaI7gtE9gIiOgPOXtLFLjF+xxORbcpJFD9zOIpshW8U4PsM0PkZcWJy0oHbE5NTjotDWmkngd1fAJABQa2BRt2ARvHiZ69wqf8FFsPkhohskza58TJTjRsiKSkcgcge4ke/6UBhJnBle8Vcnb3iEFbaKfHj4PfiOZ7hFclORcIT2FKsrmyDmNwQkW0ydwE/Imvi5g/EPip+AEB+mlgZOXG/mOykHBfr65xaJX4AYmHBiM63k52wOHG/LBvA5IaIbFOONrnhsBTZIY8gsWBkq8Hi87JCcV+zxH1i0pN0ACjNBS5tFj8AsdJ2aHtxRZZ2KMvNX7J/Qn0wuSEi28RhKaLblG5Ak97iBwCoy4H002Kyk7hP7N3JTwGSD4ofe78S2/k1Axppk514wLdJg6jdxOSGiGyPIFQalmLPDVEVCgcgpJ340fU58XsmJ/F2opO0H0g/I24pcusicPRH8Ty3ALFHJ6JiKCukrTj/x8owuSEi21N0C1AViY/taIUIUZ3JZIBPY/Gj3XDxWFGW2IuTuFdMem4cESson/1L/ADEjWDDO92eqBzeBXCWvoYckxsisj3aISn3YMDBSdpYiBoqV1+g+QDxAwDKS4Gbx24nO0n7gOJssdDgtZ1iG5lcXILeJAG47yOpImdyQ0Q2iCuliEzPwali/k1X8blGIw5ZaZOdxL1A9jUg9SSglHYPLCY3RGR7WMCPyPzkciAgRvyIe1o8lp8qJjoOzpKGxuSGiGxPygnxs28TaeMgsjcewUDrIVJHAfvZTpSI7INGI5akB8RxfyKyO0xuiMi23DwqrpZy8hSLkRGR3WFyQ0S25eJG8XOTBKusv0FE5sfkhohsy8UN4udm90kbBxFJhskNEdmOM3+Iw1JyB6DZvVJHQ0QSYXJDRLahOBv4e6r4uMckcdUGEdklq0huvv76a0RGRsLZ2Rldu3bFgQMHam3/66+/okWLFnB2dkZsbCz+/vtvC0VKRFbp1mVg6WCgIE3c6K/XVKkjIiIJSV7n5pdffsHkyZMxf/58dO3aFfPmzcOAAQNw/vx5BAYGVmm/Z88ejBgxAjNnzsRDDz2E5cuXY8iQIThy5AjatGkjwb+AiGqk0QAaFaBWVXwur/S8vNLxO59XaqcuA1TFQHmJ+FlVDJQXiyuiCjLEhCb9rHjMxRcY9i3gKG0BMSKSlkwQBEHKALp27YrOnTvjq6/E7dU1Gg0iIiLw8ssv46233qrSfvjw4SgsLMTatWt1x7p164b27dtj/vz5VdqXlpaitLRU9zwvLw8RERHIzc2Fp6cJN/dKPwv8Ncl017MYSf/7607at20d3RFzlX+DYNhr9T63HtcVNBUfavGzptJzvcSkIjkRNLCYxj2AYd8BXmGWuycRWUxeXh68vLwM+v0tac9NWVkZDh8+jGnTpumOyeVy9O/fH3v37q32nL1792Ly5Ml6xwYMGIDff/+92vYzZ87EjBkzTBZzjUoLxE3EiKh2cgdA7igu05Y7VHx2BBQ1HXcEHF3Ecu6OrmKvjIML4OIDuAcC7kGAZygQ3FYsB09Edk/S5CYzMxNqtRpBQUF6x4OCgnDu3Llqz0lNTa22fWpqarXtp02bppcMaXtuTM6vKTD8R9Nf1yJkUgdQN7KGGPcdMVf5N8gMe+2u5xpzX2OuKxc/5ApApqh4rD3mUE2SUk3S0iD/34ioIZF8zo25OTk5wcnJyfw3cvUFWg4y/32IiIioVpL24fr7+0OhUCAtLU3veFpaGoKDq1/GGRwcbFR7IiIisi+SJjdKpRJxcXHYsmWL7phGo8GWLVsQHx9f7Tnx8fF67QFg06ZNNbYnIiIi+yL5sNTkyZMxZswYdOrUCV26dMG8efNQWFiIsWPHAgBGjx6NsLAwzJw5EwDw6quvonfv3pg9ezYefPBBrFixAocOHcKCBQuk/GcQERGRlZA8uRk+fDgyMjIwffp0pKamon379li/fr1u0nBiYiLklVZAdO/eHcuXL8c777yDt99+G82aNcPvv//OGjdEREQEwArq3FiaMevkiYiIyDoY8/ubRSGIiIjIpjC5ISIiIpvC5IaIiIhsCpMbIiIisilMboiIiMimMLkhIiIim8LkhoiIiGwKkxsiIiKyKZJXKLY0bc3CvLw8iSMhIiIiQ2l/bxtSe9jukpv8/HwAQEREhMSREBERkbHy8/Ph5eVVaxu7235Bo9Hg5s2b8PDwgEwmQ+fOnXHw4EGL3NvU96rP9Yw919D2hrS7W5uaXs/Ly0NERASSkpIa3NYZ9vg+q8t5fJ/VD99npj2H77PqSfU+EwQB+fn5CA0N1dtzsjp213Mjl8sRHh6ue65QKCz2xjL1vepzPWPPNbS9Ie3u1uZur3t6eja4Hwb2+D6ry3l8n9UP32emPYfvs+pJ+T67W4+Nlt1PKJ44cWKDvVd9rmfsuYa2N6Td3dpY8v/EUuzxfVaX8/g+qx++z0x7Dt9n1WsI7zO7G5aihos7upMl8H1GlsD3mXnZfc8NNRxOTk5477334OTkJHUoZMP4PiNL4PvMvNhzQ0RERDaFPTdERERkU5jcEBERkU1hckNEREQ2hckNERER2RQmN0RERGRTmNxQg5eTk4NOnTqhffv2aNOmDb777jupQyIbVlRUhMaNG2PKlClSh0I2KjIyEm3btkX79u3Rp08fqcNpkOxu+wWyPR4eHtixYwdcXV1RWFiINm3aYNiwYfDz85M6NLJBH3/8Mbp16yZ1GGTj9uzZA3d3d6nDaLDYc0MNnkKhgKurKwCgtLQUgiCA5ZvIHC5evIhz585h4MCBUodCRLVgckOS27FjBwYNGoTQ0FDIZDL8/vvvVdp8/fXXiIyMhLOzM7p27YoDBw7ovZ6Tk4N27dohPDwcU6dOhb+/v4Wip4bCFO+zKVOmYObMmRaKmBoiU7zPZDIZevfujc6dO+Onn36yUOS2hckNSa6wsBDt2rXD119/Xe3rv/zyCyZPnoz33nsPR44cQbt27TBgwACkp6fr2nh7e+P48eO4evUqli9fjrS0NEuFTw1Efd9nf/zxB5o3b47mzZtbMmxqYEzx82zXrl04fPgw/vzzT3zyySc4ceKEpcK3HQKRFQEgrFmzRu9Yly5dhIkTJ+qeq9VqITQ0VJg5c2a113jhhReEX3/91ZxhUgNXl/fZW2+9JYSHhwuNGzcW/Pz8BE9PT2HGjBmWDJsaGFP8PJsyZYqwaNEiM0Zpm9hzQ1atrKwMhw8fRv/+/XXH5HI5+vfvj7179wIA0tLSkJ+fDwDIzc3Fjh07EBMTI0m81DAZ8j6bOXMmkpKScO3aNXz++eeYMGECpk+fLlXI1AAZ8j4rLCzU/TwrKCjA1q1b0bp1a0nibci4WoqsWmZmJtRqNYKCgvSOBwUF4dy5cwCA69ev49lnn9VNJH755ZcRGxsrRbjUQBnyPiOqL0PeZ2lpaRg6dCgAQK1WY8KECejcubPFY23omNxQg9elSxccO3ZM6jDIjjz99NNSh0A2qkmTJjh+/LjUYTR4HJYiq+bv7w+FQlFlgnBaWhqCg4MliopsDd9nZAl8n1kOkxuyakqlEnFxcdiyZYvumEajwZYtWxAfHy9hZGRL+D4jS+D7zHI4LEWSKygowKVLl3TPr169imPHjsHX1xeNGjXC5MmTMWbMGHTq1AldunTBvHnzUFhYiLFjx0oYNTU0fJ+RJfB9ZiWkXq5FtG3bNgFAlY8xY8bo2nz55ZdCo0aNBKVSKXTp0kXYt2+fdAFTg8T3GVkC32fWQSYIrFNPREREtoNzboiIiMimMLkhIiIim8LkhoiIiGwKkxsiIiKyKUxuiIiIyKYwuSEiIiKbwuSGiIiIbAqTGyIiIrIpTG6IiIjIpjC5ISLJyGQy/P7771KHAQB4//330b59+zqd+9RTT+GTTz4xbUDVeOutt/Dyyy+b/T5EDR2TGyKyO6ZMqo4fP46///4br7zyikmuV5spU6ZgyZIluHLlitnvRdSQMbkhIqqHL7/8Eo899hjc3d3Nfi9/f38MGDAA33zzjdnvRdSQMbkhsgNr166Ft7c31Go1AODYsWOQyWR46623dG2eeeYZjBo1CgBw69YtjBgxAmFhYXB1dUVsbCx+/vlnXdsFCxYgNDQUGo1G7z6DBw/GuHHjdM//+OMPdOzYEc7OzmjSpAlmzJiB8vLyGuNMSkrC448/Dm9vb/j6+mLw4MG4du2a7vWnn34aQ4YMweeff46QkBD4+flh4sSJUKlUujYpKSl48MEH4eLigqioKCxfvhyRkZGYN28eACAyMhIAMHToUMhkMt1zrWXLliEyMhJeXl544oknkJ+fX2O8arUaq1atwqBBg/SOR0ZG4pNPPsG4cePg4eGBRo0aYcGCBbrXr127BplMhpUrV6Jnz55wcXFB586dceHCBRw8eBCdOnWCu7s7Bg4ciIyMDL1rDxo0CCtWrKgxJiICIPW25ERkfjk5OYJcLhcOHjwoCIIgzJs3T/D39xe6du2qaxMdHS189913giAIQnJysjBr1izh6NGjwuXLl4X//ve/gkKhEPbv3y8IgiBkZWUJSqVS2Lx5s+78W7du6R3bsWOH4OnpKSxevFi4fPmysHHjRiEyMlJ4//33decAENasWSMIgiCUlZUJLVu2FMaNGyecOHFCOHPmjPDkk08KMTExQmlpqSAIgjBmzBjB09NTeP7554WzZ88Kf/31l+Dq6iosWLBAd83+/fsL7du3F/bt2yccPnxY6N27t+Di4iLMnTtXEARBSE9PFwAIixYtElJSUoT09HRBEAThvffeE9zd3YVhw4YJJ0+eFHbs2CEEBwcLb7/9do1f1yNHjggAhNTUVL3jjRs3Fnx9fYWvv/5auHjxojBz5kxBLpcL586dEwRBEK5evSoAEFq0aCGsX79eOHPmjNCtWzchLi5OSEhIEHbt2iUcOXJEiI6OFp5//nm9a589e1YAIFy9erWW/3Ei+8bkhshOdOzYUZg1a5YgCIIwZMgQ4eOPPxaUSqWQn58vJCcnCwCECxcu1Hj+gw8+KLz++uu654MHDxbGjRune/7tt98KoaGhglqtFgRBEPr16yd88sknetdYtmyZEBISonteOblZtmyZEBMTI2g0Gt3rpaWlgouLi7BhwwZBEMTkpnHjxkJ5ebmuzWOPPSYMHz5cEITbv/i1SZwgCMLFixcFALrk5s77ar333nuCq6urkJeXpzs2depUvQTwTmvWrBEUCoVezIIgJjejRo3SPddoNEJgYKDwzTffCIJwO7n5/vvvdW1+/vlnAYCwZcsW3bGZM2cKMTExetfOzc0VAAjbt2+vMS4ie8dhKSI70bt3b2zfvh2CIGDnzp0YNmwYWrZsiV27duHff/9FaGgomjVrBkAcbvnwww8RGxsLX19fuLu7Y8OGDUhMTNRdb+TIkfjtt99QWloKAPjpp5/wxBNPQC4Xf6wcP34cH3zwAdzd3XUfEyZMQEpKCoqKiqrEd/z4cVy6dAkeHh669r6+vigpKcHly5d17Vq3bg2FQqF7HhISgvT0dADA+fPn4eDggI4dO+pej46Oho+Pj0Ffo8jISHh4eFR77eoUFxfDyckJMpmsymtt27bVPZbJZAgODq5yrcptgoKCAACxsbF6x+48x8XFBQCq/RoSkchB6gCIyDISEhKwcOFCHD9+HI6OjmjRogUSEhKwfft2ZGdno3fv3rq2s2bNwhdffIF58+YhNjYWbm5umDRpEsrKynRtBg0aBEEQsG7dOnTu3Bk7d+7E3Llzda8XFBRgxowZGDZsWJVYnJ2dqxwrKChAXFwcfvrppyqvBQQE6B47OjrqvSaTyarM/akrY6/t7++PoqIilJWVQalUGn2tym20CdKdx+48JysrC4D+14SI9DG5IbITPXv2RH5+PubOnatLZBISEvCf//wH2dnZeP3113Vtd+/ejcGDB+smGGs0Gly4cAGtWrXStXF2dsawYcPw008/4dKlS4iJidHrMenYsSPOnz+P6Ohog+Lr2LEjfvnlFwQGBsLT07NO/8aYmBiUl5fj6NGjiIuLAwBcunQJ2dnZeu0cHR11k6vrQ1sX58yZM3WukWOsU6dOwdHREa1bt7bI/YgaIg5LEdkJHx8ftG3bFj/99BMSEhIAAL169cKRI0dw4cIFvZ6bZs2aYdOmTdizZw/Onj2L5557DmlpaVWuOXLkSKxbtw4LFy7EyJEj9V6bPn06li5dihkzZuD06dM4e/YsVqxYgXfeeafa+EaOHAl/f38MHjwYO3fuxNWrV7F9+3a88sorSE5ONujf2KJFC/Tv3x/PPvssDhw4gKNHj+LZZ5+Fi4uL3tBRZGQktmzZgtTU1CqJjzECAgLQsWNH7Nq1q87XMNbOnTt1K6yIqHpMbojsSO/evaFWq3XJja+vL1q1aoXg4GDExMTo2r3zzjvo2LEjBgwYgISEBAQHB2PIkCFVrte3b1/4+vri/PnzePLJJ/VeGzBgANauXYuNGzeic+fO6NatG+bOnYvGjRtXG5urqyt27NiBRo0a6eYDjR8/HiUlJUb15CxduhRBQUHo1asXhg4digkTJsDDw0NvKGz27NnYtGkTIiIi0KFDB4OvXZ1nnnmm2qE0c1mxYgUmTJhgsfsRNUQyQRAEqYMgIjKX5ORkREREYPPmzejXr5/Jr19cXIyYmBj88ssviI+PN/n1K/vnn3/w+uuv48SJE3Bw4KwCoprwu4OIbMrWrVtRUFCA2NhYpKSk4I033kBkZCR69epllvu5uLhg6dKlyMzMNMv1KyssLMSiRYuY2BDdBXtuiMimbNiwAa+//jquXLkCDw8PdO/eHfPmzatxOIyIbA+TGyIiIrIpnFBMRERENoXJDREREdkUJjdERERkU5jcEBERkU1hckNEREQ2hckNERER2RQmN0RERGRTmNwQERGRTfl/SckcpFeqF4QAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = pdb.refraction_index\n",
    "mwav = pdb.refraction_index_wavelength_nm\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(mwav,m.real, label=\"real part\", ls=\"dashed\")\n",
    "plt.plot(mwav,m.imag, label=\"imaginary part\")\n",
    "plt.legend()\n",
    "plt.xscale(\"log\")\n",
    "plt.xlabel(\"wavelength (nm)\")\n",
    "plt.ylabel(\"refractive index\")\n",
    "plt.title(\"MgSiO3\")\n",
    "plt.savefig(\"mgsio3_refractive_index.png\")\n",
    "#plt.savefig(\"mgsio3_refractive_index.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute the opacity of the condensate for the inicident light of 2 $\\mu\\mathrm{m}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Incident light:  268.0 nm =  0.268 um\n"
     ]
    }
   ],
   "source": [
    "imie = 195\n",
    "print(\"Incident light: \" ,mwav[imie],\"nm = \", mwav[imie]*1.e-3,\"um\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "rg_um = 0.05  # 0.1um = 100nm\n",
    "sigmag = 2.0\n",
    "cm2um = 1.0e4\n",
    "cm2nm = 1.0e7\n",
    "\n",
    "rg = rg_um / cm2um  # in cgs\n",
    "rg_nm = rg * cm2nm\n",
    "\n",
    "N0 = 1.0 # cm-3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Qext': 0.6246146481677305,\n",
       " 'Qsca': 0.6222127380166314,\n",
       " 'Qabs': 0.0024019101510991403,\n",
       " 'g': 0.31392139184771023,\n",
       " 'Qpr': 0.4292887594241749,\n",
       " 'Qback': 0.3552936962793397,\n",
       " 'Qratio': 0.5710164298658941}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from PyMieScatt import MieQ\n",
    "MieQ(m[imie], mwav[imie], 2.0 * rg_nm, asDict=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "268.7035\n"
     ]
    }
   ],
   "source": [
    "from exojax.special.lognormal import cubeweighted_mean\n",
    "rmean = cubeweighted_mean(rg, sigmag)\n",
    "print(rmean*cm2nm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`PyMieScatt.Mie_Lognormal` can be used to compute the extinction coefficient [Mm-1] etc.\n",
    "Note that the integration range lower - upper [nm] is very important. Be sure if the range is sufficient to cover the main body of the lognormal distribution "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from PyMieScatt import Mie_Lognormal\n",
    "from exojax.spec.mie import mie_lognormal_pymiescatt\n",
    "from exojax.spec.mie import auto_rgrid\n",
    "\n",
    "rgrid = auto_rgrid(rg_nm, sigmag)\n",
    "coeff = mie_lognormal_pymiescatt(\n",
    "    m[imie], mwav[imie], sigmag, rg_nm, N0, rgrid\n",
    ")  # geoMean is a diameter in PyMieScatt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.056542985376912956,\n",
       " 0.0562450133928908,\n",
       " 0.00029797198402215647,\n",
       " 0.6040582321868927,\n",
       " 0.022567722017475235,\n",
       " 0.07088795617178653,\n",
       " 0.029948661060947998)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coeff #Bext, Bsca, Babs, bigG, Bpr, Bback, Bratio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do not forget the unit of `Bext`, `Bsca`, and `Babs` is $\\mathrm{Mm}^{-1}$ i.e. the inverse of mega meter. To convert the values in cgs ($\\mathrm{cm^{-1}}$), just multiply $10^{-8}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_ext = coeff[0]*1.e-8 #Mm-1 to cm-1 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computes the optical depth for L = 10 km and $n = 10^7 \\mathrm{cm^{-3}}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from exojax.atm.amclouds import geometric_radius\n",
    "rgeo = geometric_radius(rg, sigmag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming the large size limit ($Q_e = 2$), we estimate the extinction coefficient from the geometric radius. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(4.106162e-10, dtype=float32, weak_type=True)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import jax.numpy as jnp\n",
    "Qe = 2 # large size limit\n",
    "Qe*jnp.pi*rgeo**2 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is close to the extinction coefficient computed using `PyMieScatt`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.654298537691296e-10"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta_ext "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}