HajimeKawahara/exojax

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Terminal velocity of Cloud particles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we compute a terminal velocity as a function of the condensate (cloud) radius."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import exojax.atm.viscosity as vc\n",
    "import jax.numpy as jnp\n",
    "\n",
    "g=980.0 #gravity cm/s2\n",
    "drho=1.0 #condensate density g/cm3\n",
    "rho=1.29*1.e-3 #atmosphere density g/cm3\n",
    "vfactor,Tr=vc.calc_vfactor(atm=\"Air\") #we use Air\n",
    "eta=vc.eta_Rosner(300.0,vfactor) #dynamic viscosity at 300K\n",
    "r=jnp.logspace(-5,0,70) # radius array (cm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The terminal velocity of the cloud particle can be computed using atm.vterm.vf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from exojax.atm.vterm import vf\n",
    "vterminal=vf(r,g,eta,drho,rho)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAs3UlEQVR4nO3deXhV1bnH8e/LEGaSIPMQBkEQBAXCIA611rmiVayKWkUQRKudbgd7vR1t69D29mpVFAUjDigqVlAU66wFDKOMooAMYZA5zGR67x/nYNNIwklyTvYZfp/nydOclX32fpen5M3aa+31mrsjIiJSnlpBByAiIvFNiUJERCqkRCEiIhVSohARkQopUYiISIWUKEREpEJ1gg4gFpo3b+6dOnUKOgwRkYQyf/787e7eomx7UiaKTp06MW/evKDDEBFJKGa27mjtSXXrycyGmtn4/Pz8oEMREUkaSZUo3H26u49JT08POhQRkaSRVIlCRESiL6kShW49iYhEX1IlCt16EhGJvqRKFCIiEn1KFCIiSWBz/kFeW7yZfYeLon7upHyOQkQkmRUUlbB88x4WrNvF/PW7WLhuF5vyDwHw1KiBnNHta8/MVUtSJQozGwoM7dq1a9ChiIhEzba9h1mwfhcL1u1iwfpdLM7L53BRCQBt0+vTr2MmN2Vl0q9jJj3bNI369S0ZK9xlZ2e7nswWkURUVFzCp1v2fpUY5q/fxYadBwGoW9s4qV06/bIyQ18dM2iT3iBq1zaz+e6eXbY9qUYUIiKJZse+wyxcv5v54cSwOC+fg4XFALRqWo9+WZlcP7gT/Tpm0KttOvXr1q7xGJUoRERqSHGJs/LIaCGcGNbuOABAnVpGr7ZNuWpAB/p1zKRfVgbtMhpgZgFHrUQhIhIz+QcKWbAhNNk8f/0uPtmQ/9WqpOaN0+iblclVA7Lo3zGT3u3SaZBW86OFSMR9ojCzE4EfAs2Bt919XMAhiYh8TUmJs3rbPuaHJ5wXrN/Nqq37AKhl0KN1Uy7r245+HTPol5VJVrOGcTFaiEQgicLMJgIXA1vd/aRS7RcA9wO1gcfd/R53XwGMNbNawCRAiUJEArf3UCGLNuxmwbrQ/MKi9bvYcyg0WshoWJd+WZl855S29OuYycntM2hUL+7/Li9XUJHnAA8S+sUPgJnVBh4CzgXygLlmNs3dl5vZJcAtwFMBxCoiKc7dWbN9/1fLUxes281nW/fiDmbQvVUTvt2nLf2yMujXMZMuzRslzGghEoEkCnf/wMw6lWkeCKxy9zUAZvYccCmw3N2nAdPM7DXg2aOd08zGAGMAsrKyYhW6iKSA/YeL+GTD7q9uIS1Yv4vdBwoBaFK/Dv2yMrmodxv6dczglA4ZNKlfN+CIYyuexkLtgA2lXucBg8zsLOByoB4wo7w3u/t4YDyEnqOIWZQiklTcnfU7D3w1Upi/bhefbtlDSfi3SNeWjTmvZyv6dww9u3B8i8bUqpU8o4VIxFOiOCp3fw94L5Jj9WS2iBzLwYJiFuft/mqksHD9LrbvKwCgUVpt+mZlcts3u9K3YyZ9O2SQ0TAt4IiDF0+JYiPQodTr9uE2EZEqcXc27j4YSgrh+YXlm/ZQFB4udG7eiDNPaPHVaOGEVk2onWKjhUjEU6KYC3Qzs86EEsTVwDWVOYG7TwemZ2dnj45BfCIS5w4VFrNsU35oieq60Ihh697DADSoW5uTO6Qz5swu9MvKpG9WBsc1rhdwxIkhqOWxk4GzgOZmlgf8xt0nmNltwExCy2MnuvuySp5Xt55EUsjm/INfJYT563axbFM+hcWh0UKHZg0Ycvxx4aecM+nRugl1aquyQlUEteppeDntM6hgwlpEUtueQ4W8tngzH63a/h9ba9erU4s+7dMZeVpn+nUMjRZaNqkfcLTJQ7vHikhcKylxZq/ZwZR5G3hj6RYOF5XQJr3+V/MK/TtmcmKbpqTV0WihurR7rIgklN0HCpicu4Gn56xj4+6DNK1fh+9mt+fK7A70bpeeVA+0xbukShSaoxBJfJ99uZcn/rWWlxfmcaiwhFO7HMcvLuzBeT1bBbLFtiRZotCqJ5HEVFLivP/5NiZ+9AUffr6denVqcVnfdow4rRM9Wke/YptUTlIlCo0oRBLLocJiXl64kQkffcGqrfto1bQePzu/O8MHZtGskR50ixeazBaRGrdzfwGTZq9l0ux17NxfQK+2TbnpjM58u3dbTUoHSJPZIhK4dTv28/iHX/DC/A0cKizhWz1aMvrMLgzq3EyT03FMiUJEYm5JXj6PvL+aGUs3U7dWLb7Tty2jz+hCt1ZNgg5NIpBUiUJzFCLxw92ZvXoH495fzYefb6dJvTrcfObxjDytEy2b6mG4RKI5ChGJqpIS560VX/LQe6v5ZMNuWjSpx6jTO3PNoCyaJnndhkSnOQoRiami4hJeW7KZh99dzcov95LVrCF/vOwkhvVrr+cfEpwShYhUS2FxCS8v3MjD765i7Y4DdGvZmP+76hQu7tNGm/AliaRKFJqjEKk5BUUlvDg/j4ffW0XeroOc1K4pj1zXn/N6tkq5CnDJTnMUIlIph4uKmTIvj3HvrmJT/iFO6ZDBD7/VjbO6t9AS1wSnOQoRqZaCohKmzNvAw+EE0b9jJvcM68MZ3ZorQSQ5JQoRqVBhcegW04PvrGLj7oP0y8rg3iv6cHpXJYhUoUQhIkdVXOL8Y+FG7n/7c9bvPMApHTL40+W9OVMjiJQT94nCzL4DfBtoCkxw9zeDjUgkuZWUODOWbuZv//yM1dv206ttUyaOyOab3VsqQaSooGpmTwQuBra6+0ml2i8A7idUM/txd7/H3f8B/MPMMoG/AEoUIjHg7ry3chv3zVzJis176NayMY9c14/zerbWKqYUF9SIIgd4EJh0pMHMagMPAecCecBcM5vm7svDh/xP+OciEmVz1+7kvjc+Ze7aXWQ1a8jfrjqZS05uR20lCCGgROHuH5hZpzLNA4FV7r4GwMyeAy41sxXAPcDr7r6gZiMVSW6fbtnDfW+s5J1Pt9KiST3u+s5JXJXdQVt9y3+IpzmKdsCGUq/zgEHA7cA5QLqZdXX3R472ZjMbA4wByMrKinGoIolt4+6D/O+bnzF1YR5N6tXh5xd058YhnWmQpq025OviKVEclbs/ADwQwXHjzWwzMDQtLa1/7CMTSTy7DxTw0LureHL2OgDGnNGFW846noyGqiYn5YunRLER6FDqdftwm4hU06HCYibNXsuD76xi7+EirujXnh+fewJtMxoEHZokgHhKFHOBbmbWmVCCuBq4pjIncPfpwPTs7OzRMYhPJOGUlDjTF2/izzNXkrfrIGd1b8EdF/agR+umQYcmCSSo5bGTgbOA5maWB/zG3SeY2W3ATELLYye6+7JKnlebAoqE5X6xkz+8tpzFefn0bNOUZ27qw2ldmwcdliSgCjcFNLP2hP6yPwNoCxwElgKvEVqFVFITQVaWNgWUVLZux37uef1TXl+6hTbp9fnped25rG87PQshx1TpTQHN7AlCK5FeBe4FtgL1gROAC4A7zewOd/8gNiFXnkYUksr2HCrkwXdWkfOvtdSpbfzXuSdw0xldtJJJqq3cEYWZneTuS8t9o1kakOXuq2IVXFVpRCGppLjEeX7uBv765kp2Hijgu/3b81/ndaeV6lJLJVV6RHG0JBHeRqODuy929wIgrpKERhSSamav3sHvX13Ois17GNipGU8O7clJ7dKDDkuSzDELF5nZe8AlhJLKfEK3oGa5+49jHl0VaUQhyS5v1wH+NGMFM5ZsoV1GA/77ohO5qHdrbdon1VKdwkXp7r7HzG4CJrn7b8xscfRDrD6NKCTZHSos5pH3VzPuvdWYwU/OPYExZ3ahfl3NQ0jsRJIo6phZG+BK4M4Yx1Mteo5CkpW788bSLfzhtRVs3H2Qi/u04b8vOlEPzEmNiCRR/J7Qsw0fuftcM+sCfB7bsETkiNXb9vHbacv48PPt9GjdhOfGDGZwl+OCDktSSEXLY4cDb7r7C8ALR9rDu7sOq4HYRFLa/sNF/P2dVUz4aA3169bmt0N7ct3gjtSprZ1dpWZVNKLIAl4ws7rA28DrQK4fa/Y7QJqjkGQxZ80Ofvz8IjbnH+KK/u35xQU9aNGkXtBhSYqKZNVTE0LbfF9AqGbECuANYKa7fxnzCKtAq54kkb2yaCM/e2Ex7Zs14M9X9KF/x2ZBhyQposqrntx9L/By+Asz6wlcSKg63flRjlMkZbk7495fzX1vrGRQ52aM/1426Q3rBh2WSGSbAppZH6BTqeO/cHclCZEoKSou4bfTl/H0nPUMPbktf/luH+rV0ZJXiQ/HTBRmNhHoAywDjmwC6MDUGMYlkjIOFhRz++QFvLViK2O/cTw/P7+7NvCTuBLJiGKwu/eMeSRRoMlsSTQ79xcw6sm5LNqwm7su7cX3Tu0UdEgiXxPJOrvZ4XmJuOfu0919THq69rqR+Ldh5wGuGDeLZZv2MO7a/koSErciGVFMIpQstgCHAQPc3fvENDKRJLZsUz4jnpjL4cJinrlpEAM6aWWTxK9IEsUE4HvAEv49RyEiVTRr1XbGPDWfpvXr8OwtQ+jWqknQIYlUKJJEsc3dp8U8knKEtwy5k9DmhFcEFYdINEz7ZBP/NWURXZo3JmfkANqka68miX+RzFEsNLNnzWy4mV1+5Ks6FzWziWa21cyWlmm/wMxWmtkqM7sDQluGuPuo6lxPJB48/uEafjB5IX2zMpky9lQlCUkYkYwoGhCamzivVFt1l8fmAA8Smv8AwMxqAw8B5wJ5wFwzm+buy6txHZHAlZQ4f5qxgsc/+oKLerfmf688RduCS0KJ5MnsG6N9UXf/wMw6lWkeCKwKbzqImT0HXAooUUjCKigq4WcvfsIrizZxw6kd+fXQXtTWMxKSYI5568nMnjSzjFKvM8MP4UVbO2BDqdd5QDszO87MHgH6mtkvK4hzjJnNM7N527Zti0F4IpWz91AhN+bk8sqiTfz8gu789hIlCUlMkdx66uPuu4+8cPddZtY3diH9J3ffAYyN4LjxZrYZGJqWltY/9pGJlG/rnkOMeGIun325l79+92SG9W8fdEgiVRZJoqhlZpnuvgvAzJpF+L7K2gh0KPW6fbhNJKGs2baP6yfmsnN/AY/fkM1Z3VsGHZJItUSy6umvhB64u8vM7gJmAffFIJa5QDcz62xmacDVQKWW5erJbAnawvW7GDZuFgcLinluzGAlCUkKx0wU7j4JuBz4Mvx1ubs/VZ2LmtlkYDbQ3czyzGyUuxcBtxEqu7oCmOLuyyp53qFmNj4/P7864YlUydsrvmT4Y3No2qAuL90yhD7tM4IOSSQqyi1cZGaN3X1fhW+O4JggqHCR1LTnctfz3y8v4aR26UwcMYDmjVWNThJPeYWLKhpRvGJmfzWzM82sUakTdTGzUWY2k1DVu7ihEYXUNHfn/rc+546pSzi9Wwsmjx6sJCFJp9xE4e7fIlQr+2ZgmZntMbMdwNNAa+AGd3+xZsKMjOYopCYVFZdw5z+W8re3PuPyfu2YcEM2jerFYp2HSLAq/H+1u88AZtRQLNWmehRSU0LFhhby1oovufWs4/nZ+d0x0zMSkpwiWfWUMDSikJqw+0AB1034mLc//ZLfXdKLn1/QQ0lCkprGySKVkLfrADdMzGXDroM8dE0/LurdJuiQRGIuqRKFbj1JLK3YvIcRT+RyoKCYp0YOZFCX44IOSaRGRLLX01/NrFdNBFNduvUksTJ79Q6ufGQ2hvHi2CFKEpJSIpmjWAGMN7OPzWysmem3sKSUVxdv4oaJubROr8/UW4fQvbUq0klqieTJ7Mfd/TTgeqATsDhcyOibsQ6usvQchUTbE//6gtsnL+TkDum8MPZU2mao2JCknohWPYWLCvUIf20HPgF+Eq4ZETd060mipaTEuef1T/nd9OWc17MVT40aREbDtKDDEgnEMSezzexvwMXAO8Cf3D03/KN7zWxlLIMTCUJBUQm/eGkxLy/cyLWDsvj9pSepjoSktEhWPS0G/sfd9x/lZwOjHI9IoPYdLuKWp+fz4efb+el5J/D9b3bVMxKS8iK59XRd2SRhZm8DuLsmAyRpbNt7mOHj5zBr9Q7uG9aH287upiQhQgUjCjOrDzQEmptZJnDkX0xTQmVL446eo5CqWrt9P9dPzGXb3sM8dn1/zu7RKuiQROJGRSOKm4H5hCawF4S/nw+8AjwY+9AqT5PZUhWLNuxm2LhZ7DtcxLOjBylJiJRR7ojC3e8H7jez29397zUYk0iNeXflVm59egHNm6Tx5I0D6dKicdAhicSdim49ne3u7wAbzezysj9396kxjUwkxl6Yt4E7pi6hR+smPHHjAFo2qR90SCJxqaJVT98gtCR26FF+5oAShSQkd+ehd1fxlzc/4/SuzRl3XT+a1K8bdFgicauiW0+/Cf/vjTUXzteFq+s9DBQA77n7M0HGI4mtuMT57bRlPDVnHZee0pY/X3EyaXWSard9kaiLZFPAP5lZRqnXmWb2h+pc1MwmmtlWM1tapv0CM1tpZqvM7I5w8+XAi+4+GrikOteV1HaosJjvP7OAp+as4+Yzu/C3K09RkhCJQCT/Si50991HXrj7LuCial43hzL1tsPbhDwEXAj0BIabWU+gPbAhfFhxNa8rKSr/QCHXT8hl5vIt/PrinvzyohOppaetRSISSaKobWZfVYs3swZAtarHu/sHwM4yzQOBVe6+xt0LgOeAS4E8QsmiwnjNbIyZzTOzedu2batOeJJkNu0+yHcfncWiDbv5+/C+jDy9c9AhiSSUSLbweAZ428yeCL++EXgyBrG0498jBwgliEHAA8CDZvZtYHp5b3b38Wa2GRialpbWPwbxSQJauWUvN0zMZf/hInJGDmDI8c2DDkkk4RwzUbj7vWb2CXBOuOkud58Z27D+4/r7CSWnSI6dDkzPzs4eHduoJBF8vGYHoyfNo37d2kwZeyontmkadEgiCSnSUqgLgbqElsUujFEsG4EOpV63D7dFTFt4yBGvL9nMD59fRIfMBjw5ciDtMxsGHZJIwopk1dOVQC5wBXAl8LGZXRGDWOYC3cyss5mlAVcD02JwHUlyk2av5dZnF3BS26a8OHaIkoRINUUymX0nMMDdb3D36wlNOv+qOhc1s8nAbKC7meWZ2Sh3LwJuA2YSKr86xd2XVea82usptbk7f575Kb9+ZRnf6tGKZ24aTGYjFRsSqa5Ibj3VcvetpV7vIMLKeOVx9+HltM8AZlT1vLr1lLoKi0v45dQlvDg/j+EDs7jr0l7Uqa1nJESiIZJ/SW+Y2UwzG2FmI4DXqMYv81jSiCI1HSgoYvSkebw4P48fn3MCf7rsJCUJkSiKZNXTz8xsGHBauGm8u78c27CqRiOK1LNj32FG5sxlycZ87r68N8MHZgUdkkjSMXcPOoaoy87O9nnz5gUdhsTYuh37uWFiLlv2HOLB4f04p6fqSIhUh5nNd/fssu0VbTO+l9By2K/9CHB316J0CcySvHxuzMmlqMR55qbB9O+YGXRIIkmrot1jm9RkINGgW0+p4YPPtjH26flkNkzjuZED6dpSxYZEYimiGT8zO93Mbgx/39zM4nKzHE1mJ7+pC/IYmTOXrGYNmXrrECUJkRpwzMlsM/sNkA10B54A0oCn+ffktkjMuTuPfrCGe17/lFO7HMej1/enqYoNidSISJ6juAzoCywAcPdNZhaXt6V06yk5lZQ4v391OTmz1jL05Lb85bt9qFendtBhiaSMSG49FXhoaZTDVxXn4pJuPSWfQ4XF3D55ITmz1nLT6Z25/6pTlCREalgkI4opZvYokGFmo4GRwGOxDUsE8g8WMmbSPD7+Yid3XnQio8/sEnRIIikpkgfu/mJm5wJ7CM1T/Nrd/xnzyCSlbck/xIgnclm9bR/3X30Kl57SLuiQRFJWJJPZPwGeV3KQmrJq616un5DLnkNF5Nw4kNO6qtiQSJAimaNoArxpZh+a2W1mFrePv5rZUDMbn5+fH3QoUkXz1u5k2LjZFBQ7z40ZrCQhEgeOmSjc/Xfu3gv4PtAGeN/M3op5ZFWgyezENnPZFq59/GOaNUrj5VuHcFI7fY4i8SDSCncAW4EthLYZbxmbcCRVPT1nHb9+ZSl92mcwccQAmqmOhEjciGSO4lZCle1aAC8Ao919eawDk9Tg7vztn5/xwDurOLtHSx68pi8N0yrz94uIxFok/yI7AD9y90UxjkVSTFFxCXe+vJTn523gyuz2/Omy3qojIRKHIlke+8uaCKQ8ZtaFUDnWdHePRa1uCcCBgiJue3Yh73y6lR+c3ZUfn3sCZhZ0WCJyFDH9883MJprZVjNbWqb9AjNbaWarzOyOis7h7mvcfVQs45SatXN/Adc89jHvrdzKH75zEj85r7uShEgci/XN4BzgQWDSkQYzqw08BJwL5AFzzWwaUBu4u8z7R5ap1y0JbsPOA9wwMZeNuw8y7rr+nN+rddAhicgxxDRRuPsHZtapTPNAYJW7rwEws+eAS939buDiql7LzMYAYwCyslQOMx4t25TPiCfmcriwmKdvGsSATs2CDklEIlDurScz22tme47ytdfM9lTjmu2ADaVe54XbyovjODN7BOhrZuXOl7j7eHfPdvfsFi1aVCM8iYV/rdrOVY/OoU4t46VbhihJiCSQuK9w5+47gLGRHKttxuPTK4s28tMXPqFL88bkjBxAm/QGQYckIpUQ8a0nM2sJ1D/y2t3XV/GaGwktuT2ifbhNktBjH6zhjzNWMKhzM8Zfn016AxUbEkk0x1z1ZGaXmNnnwBfA+8Ba4PVqXHMu0M3MOptZGnA1MK0a5/uKtvCIHyUlzh9eXc4fZ6zg273b8OTIgUoSIgkqkuWxdwGDgc/cvTPwLWBOJCc3s8nAbKC7meWZ2Sh3LwJuA2YCK4Ap7r6sStF//XraFDAOHC4q5kfPL+Lxj75gxJBO/H14X+rXVbEhkUQVya2nQnffYWa1zKyWu79rZv8XycndfXg57TOAGZWIUxLE3kOF3PzUfGat3sEdF/bg5jO76BkJkQQXyYhit5k1Bj4AnjGz+4H9sQ2ranTrKVhb9xziykfnkPvFTv73ypMZ+43jlSREkkAkI4pLgUPAj4FrgXTg97EMShLP6m37uH5CLrsOFDBhxAC+cYKWKIski0j2eio9engyhrFUm5bHBmP+ul2MenIudWoZz40ZTJ/2GUGHJCJRFMmqp8vN7HMzy4/SA3cxo1tPNe/tFV9y7eNzSG9Ql5duGaIkIZKEIpmjuA+4xN3T3b2puzdx96axDqwqtOqpZj2Xu57Rk+ZxQqsmvHTLEDoe1yjokEQkBiJJFF+6+4qYRxIFGlHUDHfn/rc+546pSzijWwsmjx5M88b1gg5LRGIkksnseWb2PPAP4PCRRnefGqugJH4VFZfwq1eWMTl3PcP6teeeYb2pq2JDIkktkkTRFDgAnFeqzQElihRzsKCY2ycv5K0VX/L9bx7PT1VHQiQlRLLq6caaCCQatOopdnbtL2DUk3NZuGE3v7ukFzcM6RR0SCJSQ8pNFGb2c3e/z8z+TmgE8R/c/QcxjawK3H06MD07O3t00LEkk7xdoWJDG3Ye5KFr+nFR7zZBhyQiNaiiEcWRCex5NRGIxKflm/Yw4olcDhYWM2nUQAZ3OS7okESkhlVUj2J6+H/j+iE7iZ1Zq7dz86T5NKpXhxfGnkqP1nG5KlpEYuyYcxRmlg3cCXQsfby794lhXFWiOYroeXXxJn7y/Cd0PK4hT44cSNsMFRsSSVXm/rXph/88wGwl8DNgCVBypN3d18U2tKrLzs72efN0x6yqJn70BXe9tpzsjpk8dn02GQ3Tgg5JRGqAmc139+yy7ZEsj93m7lEpLCTxraTEufeNT3n0gzWc36sV91+tOhIiElmi+I2ZPQ68jR64S1oFRSX8/MVP+MeiTXxvcEd+e0kvatfSMxIiElmiuBHoAdTl37ee9MBdEtl3uIhbnp7Ph59v56fnncD3v9lVD9KJyFciSRQD3L17zCOpgJl9B/g2oafEJ7j7m0HGk0y27j3EyJy5rNi8l/uu6MOV2R2CDklE4kwkm/TMMrOeVb2AmU00s61mtrRM+wVmttLMVpnZHRWdw93/4e6jgbHAVVWNRf7TF9v3M2zcLFZv3c/j12crSYjIUUUyohgMLDKzLwjNURjglVgemwM8CEw60mBmtYGHgHOBPGCumU0DagN3l3n/SHffGv7+f8Lvk2patGE3I3PmAjB5zGBO6ZARbEAiErcqTBQWulF9M1DlpbDu/oGZdSrTPBBY5e5rwtd5DrjU3e8GLi4njnuA1919QVVjkZB3V27l1qcX0LxJGpNGDqJzc9WREJHyVZgo3N3N7CF37x3l67YDNpR6nQcMquD424FzgHQz6+ruj5Q9wMzGAGMAsrKyohhqcnlh3gbumLqEE9s04YkRA2nRRHUkRKRikdx6WmBmA9x9bsyjKYe7PwA8cIxjxpvZZmBoWlpa/5qJLHG4Ow+9u4q/vPkZZ3Rrzrjr+tO4XiQfv4ikukh+UwwCrjOztcB+Kj9HcTQbgdIzp+3DbRIDxSXOb6ct46k567isbzvuHdaHtDoqNiQikYkkUZwfg+vOBbqZWWdCCeJq4JrqnlTbjH/docJifvTcIt5YtoWbv9GFX5zfg1p6kE5EKuGYf1aG93TqAJwd/v5AJO87wswmA7OB7maWZ2aj3L0IuA2YSWg78ynuvqwqHShzraFmNj4/P7+6p0oK+QcK+d6Ej5m5fAu/vrgnv7zwRCUJEam0SDYF/A2QDXR39xPMrC3wgrufVhMBVoU2BYRNuw9yw8Rc1u04wF+vPJmhJ7cNOiQRiXPV2RTwMqAvsADA3TeZWZMoxxcV2mY8ZOWWvdwwMZf9h4vIGTmAIcc3DzokEUlgkdxCKvDQsMMBzCxuF927+3R3H5Oenh50KIGZs2YHVzwyixJ3pow9VUlCRKotkkQxxcweBTLMbDTwFvB4bMOqmlSfo5ixZDPXT8ylRZN6TL11CCe2UUU6Eam+Y85RAJjZucB5hJbGznT3f8Y6sOpIxTmKJ2et5bfTl9G3QwYTbhhAZiMVGxKRyqnyHIWZ3evuvwD+eZS2uJKKcxTuzp9nruTh91Zzzomt+PvwvjRIU7EhEYmeSG49nXuUtgujHUg0pNocRWFxCT99YTEPv7ea4QOzeOS6fkoSIhJ15Y4ozOwW4Fagi5ktLvWjJsC/Yh2YVGz/4SJufWYB73+2jR+fcwI/+JaKDYlIbFR06+lZ4HVC236Xrhex1913xjQqqdD2fYcZmTOXpRvzuefy3lw9UJsgikjslJso3D0fyAeG11w41ZMKcxTrduzn+om5fLnnEOO/l805PVsFHZKIJLmk2hku2ecoluTlM2zcLPYcLOTZ0YOVJESkRmif6QTxwWfbGPv0fDIbpjFp1ECOb9E46JBEJEUoUSSAqQvy+PmLi+nWqglP3jiAlk3rBx2SiKQQJYo45u488v4a7n3jU07rehyPXNefJvXrBh2WiKSYpEoUyTSZXVLi/P7V5eTMWsslJ7flL989WcWGRCQQSfWbJ1kmsw8VFnPb5AXkzFrLTad35v+uOkVJQkQCk1QjimSQf7CQMZPm8fEXO7nzohMZfWaXoEMSkRSnRBFHtuQf4oaJuazZvo/7rz6FS09pF3RIIiLxnyjM7ETgh0Bz4G13HxdwSDHx+ZehYkN7DhWRc+NATuuqOhIiEh9ieuPbzCaa2VYzW1qm/QIzW2lmq8zsjvLeD+DuK9x9LHAlELflV6tj7tqdXPHIbApLnOdvHqwkISJxJdYzpDnABaUbzKw28BChHWh7AsPNrKeZ9TazV8t8tQy/5xLgNWBGjOOtcTOXbeG6xz/muEZpTL1lCL3aJvZEvIgkn5jeenL3D8ysU5nmgcAqd18DYGbPAZe6+93AxeWcZxowzcxeI7RZYVJ4es46fv3KUvq0z2DiiAE0U7EhEYlDQcxRtAM2lHqdBwwq72AzOwu4HKhHBSMKMxsDjAHIyorv3VTdnb/98zMeeGcVZ/doyYPX9KVhWtxPF4lIior7307u/h7wXgTHjTezzcDQtLS0/rGOq6qKikv475eXMGVeHldmt+dPl/WmTm09IyEi8SuI31AbgQ6lXrcPt1VbvD9wd6CgiDFPzWfKvDx+cHZX7h3WR0lCROJeECOKuUA3M+tMKEFcDVwTjRPH8xYeO/cXMDJnLovzdvPHy07i2kEdgw5JRCQisV4eOxmYDXQ3szwzG+XuRcBtwExgBTDF3ZfFMo6gbdh5gCvGzWLF5j2Mu66/koSIJBRz96BjiLrs7GyfN29e0GEAsHRjPjfmzKWgqIQJN2ST3alZ0CGJiByVmc139+yy7Ul1g9zMhprZ+Pz8/KBDAeCjz7dz9fg5pNWuxUu3nKokISIJKakSRTxNZr+yaCM35uTSPrMBL90yhK4tmwQdkohIlSRVooiXEcVjH6zhh88tol9WJs/ffCqt01WRTkQSV1IliqBHFCUlzl2vLuePM1bw7d5teHLkQNIbqCKdiCS2uH/gLlEcLirmpy8sZvonmxgxpBO/vrgntWpZ0GGJiFRbUiWKoJ6j2HuokJufms+s1Tu448Ie3HxmF8yUJEQkOejWUzVt3XOIKx+dQ+4XO/nfK09m7DeOV5IQkaSSVCOKmrZ62z6un5DLrgMFTBwxgDNPaBF0SCIiUZdUiaImbz3NX7eLUU/OpU4t4/kxp9K7ffBLckVEYkG3nqrgn8u/5NrH55DRoC4v3TJESUJEklpSjShqwuTc9dz58hJ6t0tnwogBNG9cL+iQRERiSokiQu7O/W9/zv+99TlndW/BQ9f0o1E9/ecTkeSn33QRKCou4VevLGVy7gau6N+euy/vTV3VkRCRFJFUiSIWk9kHC4q5ffJC3lrxJbd9syv/dd4JWv4qIiklqf4sjvZk9q79BVw34WPe/vRLfn9pL356fnclCRFJOUk1ooimvF0HuGFiLht2HeTha/pxYe82QYckIhIIJYqjWL5pDyOeyOVQYTFPjxrEwM6qIyEiqUuJooxZq7dz86T5NKpXhxfGDqF7a9WREJHUlhBzFGbWyMzmmdnFsbzO9E82MWLiXFqn12fqrUoSIiIQ40RhZhPNbKuZLS3TfoGZrTSzVWZ2RwSn+gUwJTZR/tvSTfmc0iGDF8cOoW1Gg1hfTkQkIcT61lMO8CAw6UiDmdUGHgLOBfKAuWY2DagN3F3m/SOBk4HlQMzLxP3i/B4UlpRQr07tWF9KRCRhxDRRuPsHZtapTPNAYJW7rwEws+eAS939buBrt5bM7CygEdATOGhmM9y95CjHjQHGAGRlZVUp3lq1jHq1lCREREoLYjK7HbCh1Os8YFB5B7v7nQBmNgLYfrQkET5uvJltBoampaX1j164IiKpLSEmswHcPcfdXz3GMYHWzBYRSUZBJIqNQIdSr9uH26rNzIaa2fj8/PxonE5ERAgmUcwFuplZZzNLA64GpgUQh4iIRCDWy2MnA7OB7maWZ2aj3L0IuA2YCawAprj7smhcT7eeRESiL9arnoaX0z4DmBHt69VkKVQRkVSRMJPZkdCIQkQk+szdg44h6sxsG7AOSAdKz2yXfl3e982B7VEIo+y1q3NseT+vqH9Ha0uFPkf6mSdKnyNpS9Q+R/oZH61NfT56n6vb347u3uJrre6etF/A+PJeV/D9vFhcuzrHlvfzivqXqn2uxGeeEH2OpC1R+xzpZ6w+R97naPW37FdS3Xo6iukVvC7v+1hduzrHlvfzivp3tLZU6HOkn3m0xLrPkbQlap8j/YyP1qY+x77PX0nKW0/VYWbz3D076DhqkvqcGtTn5Ber/ib7iKIqxgcdQADU59SgPie/mPRXIwoREamQRhQiIlIhJQoREamQEoWIiFRIiaISzKyLmU0wsxeDjiWWwjXKnzSzx8zs2qDjqQmp8tkeYWbfCX++z5vZeUHHUxPM7EQze8TMXjSzW4KOp6aE/z3PM7OvFYaLVMokimjU73b3Ne4+KraRxkYl+3858KK7jwYuqfFgo6QyfU7kz/aISvb3H+HPdyxwVRDxRkMl+7zC3ccCVwKnBRFvNFThd9kvgCnVuWbKJApC9bsvKN1Qqn73hYRKrQ43s55m1tvMXi3z1bLmQ46qHCLsP6EaIUeqEBbXYIzRlkPkfU4GOVS+v/8T/nmiyqESfTazS4DXiMGmpDUoh8h/l50LLAe2VueCQZRCDYRHoX53IqtM/wmVp20PLCKB/5ioZJ+X13B4UVeZ/prZCuAe4HV3X1CzkUZPZT9jd58GTDOz14BnazTYKKlknxsDjQglj4NmNsPLKSddkYT9JRAlR6vf3a68g83sODN7BOhrZr+MdXA1oLz+TwWGmdk4Yrw1QACO2uck/GyPKO8zvh04B7jCzMYGEVgMlfcZn2VmD5jZoyT2iOJojtpnd7/T3X9EKCk+VpUkASk0oogGd99B6J5uUnP3/cCNQcdRk1Llsz3C3R8AHgg6jprk7u8B7wUcRiDcPac670/1EUXM6ncniFTsf6r1OdX6C+ozRLnPqZ4oUr1+dyr2P9X6nGr9BfU56n1OmURhNVy/O96kYv9Trc+p1l9Qn2uqz9oUUEREKpQyIwoREakaJQoREamQEoWIiFRIiUJERCqkRCEiIhVSohARkQopUYiISIWUKEQqyUKFb7pU4vhsM4vpvkpm9hczOzuW15DUpQfuRMphZkbo30hJqbZewB/c/bIYXre2u1eqDoiZdSS0O2hKVKuTmqURhUgpZtYpXCVsErCU/9xoDeBa4JVSx+8zsz+b2TIze8vMBprZe2a2Jlwkh/D21q+Gv29sZk+Y2RIzW2xmw0qd569m9glwqpn9xMyWhr9+VCq2FRYqYbrMzN40swYA7r4OOM7MWsf4P5GkICUKka/rBjzs7r3Cv4BLOw2YX+p1I+Add+8F7AX+AJwLXAb8/ijn/hWQ7+693b0P8E6p83zs7icDBwlt8z4IGAyMNrO+pWJ7KHy93cCwUudeQAKX+JT4pXoUIl+3zt3nlPOzNsC2Uq8LgDfC3y8BDrt7oZktATod5f3nENrZEwB33xX+thh4Kfz96cDL4bogmNlU4AxCu4F+4e6LwsfNL3ONrUDbY/RNpNI0ohD5uv0V/OwgUL/U60L/90RfCXAYIDyvUZk/xA5FOC9xuNT3xWWuUT8cn0hUKVGIVM4KoGs13v9P4PtHXphZ5lGO+RD4jpk1NLNGhG5jfRjBuU8gNK8iElVKFCKV8xpwVjXe/wcgMzxJ/QnwzbIHuPsCIAfIBT4GHnf3hRWd1MzqEkpg86oRm8hRaXmsSCWEVxm9C5xW2SWssWRmlwH93P1XQcciyUcjCpFKcPeDwG+AdkHHUkYd4K9BByHJSSMKERGpkEYUIiJSISUKERGpkBKFiIhUSIlCREQqpEQhIiIV+n/ij/IkSGHokAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(r*1e4,vterminal)\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "\n",
    "plt.xlabel(\"r (micron)\")\n",
    "plt.ylabel(\"terminal velocity (cm/s)\")\n",
    "plt.savefig(\"vterm.pdf\", bbox_inches=\"tight\", pad_inches=0.0)\n",
    "plt.savefig(\"vterm.png\", bbox_inches=\"tight\", pad_inches=0.0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above figure, you see three regimes of the r dependence on the terminal velocity.\n",
    "For small r, the terminal velocity obeys the Stokes flow.\n",
    "In the mid - large r (i.e. medium and large Reynolds number Nre), exojax uses a phenomenological modeling as explained below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start from importing the the drag coefficients as Nre from Table 10.1 p381 in Pruppacher and Klett\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv(\"~/exojax/data/clouds/drag_force.txt\",comment=\"#\",delimiter=\",\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's fit logarithms of Davies number $N_D = C_d N_{re}^2$ and Reynolds number $N_{re}$ by a polynomial equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nre=data[\"Nre\"].values\n",
    "logNre=np.log(Nre) #Reynolds number\n",
    "Cd=(data[\"Cd_rigid\"].values)\n",
    "logNd=np.log(Nre**2*Cd)\n",
    "\n",
    "Cdinf=0.45\n",
    "Nreinf=np.logspace(3,5,30)\n",
    "logNreinf=np.log(Nreinf)\n",
    "logNdinf=np.log(Nreinf**2*Cdinf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.00883374,  0.84514511, -2.49105354])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coeff=np.polyfit(logNd,logNre,2)\n",
    "coeff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the coefficient we use in exojax in the mid Nre regime.\n",
    "\n",
    "i.e.\n",
    "\n",
    "$\\log{N_{re}} = 0.0088 \\log^2{N_{D}} + 0.85 \\log{N_{D}} + 2.49$\n",
    "\n",
    "Davies number can be computed using the following function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Davies number= 400.34301797889896\n"
     ]
    }
   ],
   "source": [
    "from exojax.atm.vterm import Ndavies\n",
    "\n",
    "g=980.0 #gravity cm/s2\n",
    "drho=1.0 #condensate density g/cm3\n",
    "rho=1.29*1.e-3 #atmosphere density g/cm3\n",
    "vfactor,Tr=vc.calc_vfactor(atm=\"Air\") #we use Air\n",
    "eta=vc.eta_Rosner(300.0,vfactor) #dynamic viscosity at 300K\n",
    "r=0.01 #cm\n",
    "print(\"Davies number=\",Ndavies(r,g,eta,drho,rho))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would obtain a boundary  between the mid Nre regime and the Stokes flow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#boundary between the Stokes flow and the mid Nre regime\n",
    "#-0.00883374*xarr**2+(0.84514511-1)*xarr-2.49105354 +log(24) = 0\n",
    "a=-0.0088 #coeff[0]\n",
    "b=0.85-1 #coeff[1]-1\n",
    "c=-2.49+np.log(24.) #coeff[2]+np.log(24.)\n",
    "logNdc=(-b-np.sqrt(b*b-4*a*c))/(2*a)\n",
    "Ndc=np.exp(logNdc)   #boundary for Davies number\n",
    "Nrec=np.exp(logNdc-np.log(24.)) #boundary for Reynolds number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.7583482270854875, 42.87754348901474, 1.7865643120422807)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logNdc, Ndc, Nrec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, for large Nre, we assume Cd=0.45 following Akerman and Marley 2001. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#boundary between the mid and large Nre regime\n",
    "#-0.00883374*xarr**2+(0.84514511-0.5)*xarr-2.49105354 +0.5*log(0.45) = 0\n",
    "a=-0.0088 #coeff[0]\n",
    "b=0.85-0.5 #coeff[1]-0.5\n",
    "c=-2.49+0.5*np.log(0.45) #coeff[2]+0.5*np.log(0.45)\n",
    "logNde=(-b+np.sqrt(b*b-4*a*c))/(2*a)\n",
    "Nde=np.exp(logNde)\n",
    "Nree=np.exp(0.5*logNde-0.5*np.log(0.45))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(11.692270778931425, 119643.38181447262, 515.629888398587)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logNde, Nde, Nree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure shows Davies number - Reynolds number relation we assume in exojax."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb8AAAENCAYAAACfEWsMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABWg0lEQVR4nO3deXxU9b3/8dc3e0hCEkLCkgQSthCyExAUZBEEFyqIWrWguNXaWuvW1lu3tre119a61/68alt6hVpbtUrFDSpYFlkS1hDWLEA2spB9n5nP748JUwJJSEjCTMjn+XjkQeacM+d8zkDy5pzzXYyIoJRSSvUnbs4uQCmllLrQNPyUUkr1Oxp+Siml+h0NP6WUUv2Ohp9SSql+x8PZBfSUwYMHS1RUlLPLUEop5ULS09NLRST0zOUXTfhFRUWRlpbm7DKUUkq5EGPM0baW621PpZRS/Y6Gn1JKqX5Hw08ppVS/o+GnlFKq39HwU0op1e+4bPgZYx42xuwzxmQYY94xxvg4uyallFIXB5cMP2NMOPADYJKIxAPuwC3OrUoppdTFwiXDr4UH4GuM8QAGAAVOrkcppc5itdooPlrFycJaZ5eiusAlO7mLSL4x5rfAMaAe+EJEvjhzO2PMvcC9ACNGjLiwRSql+rXGegvFOVWUHK/G0mQlJNyfQcP8nF3WRUNEMMb02v5dMvyMMcHAQiAaqAD+boxZKiIrTt9ORN4A3gCYNGmSzsqrlOpVIkJNeSMncqooL6pFBIKHDmDoqED8g72dXV6fZ7VYKMs7RnFuNiERkQwbE9Nrx3LJ8APmAjkiUgJgjPkAuAxY0eG7lFKqF9hsQnlhLUXZldRWNOLu6cbQ6EDCogLwHuDp7PL6NBGh5mQZxblZlB47itXSjI9/AO4evfu5umr4HQOmGmMGYL/tOQfQgTuVUhdUc5OVkqPVFB+toqnego+/JyPjBzM40h93D1duMuH6mhsaKD6aQ3FuFvVVlbi5exASEUlY1GgGhob16i1PcNHwE5Gtxpj3gB2ABdhJy+1NpZTqbfU1TZzIrqI0rwab1UZgqC9RCYMJDPPt9V/KFzOx2SgvKqQ45wjlhQWI2PAfNJjRqVMIiRyJh+eFu4p2yfADEJGfAj91dh1Kqf5BRKgqbaAou5LK4jrc3A0h4f4MiQ5kwEAvZ5fXp9VVVVKSm0PJ0WyaGurx9PZh2NgYwqJHM2BgoFNqctnwU0qpC8FmtVGWX8uJnErqqprw8HYnPCaYsJED8fR2d3Z5fZbV0kzp8WMU52ZRXVqCMW4EDRtOdNQogocNx83NuZ+thp9Sql9qbrRSfLSKE7lVWBqtDBjoRXRSKCHhfri56/O88yEiVJeVUJyTTVneMayWZnwDBjIyIYXQqGi8fHydXaKDhp9Sql+pr26iKKeKsrxqbFYhMMzeVWHgYB99nneemurrKDmaa2+8Ul2Fu4cnIZEjCYuKJiAk1CU/Vw0/pdRFr83neREBDI0eiG+APs87HzablfLCAopzs6loabwSMDiU0TFTGRw5ote7KnSXhp9S6qJlswknC2ooyq6irrJRn+f1gLqqSopzsig5mkNzYwNePr4MHxdLWPQofAMGOru8TtPwU0pddCzN9v55J3Lt/fN8A/R5XndYmpspO36UEzlZ1JwsxRg3goeHExY1muChwzBufe8z1fBTSl00GuuaKcqpouRYNTaLjYHaP++8iQhVJcUU52ZRlnccm9WC78BAopImEjoiGk+fvj3LnIafUqrPq61opDCrkvIi+8wKg4b7MXRUIH6BOt5mVzXW1VGSm03x0Wwaaqpx9/AkdGQUYVGj8R8UctH8J0LDTynVJ4kIlcX1FGZVUl1Wj7unG0OiBzI0OhAvX/3V1hU2q5XywnxO5GRRUVQICIFhQ4mIjSckYgTuHhff53nxnZFS6qJ2qlN6UXYl9dVNePl6EDkhhLARAbh79r1nT85UW1FOcW42JUdzsDQ14uU7gIjYOEKjRuHrH+Ds8nqVhp9Sqk+wNFspPlrNiewqmhstDAj0ZlRKGIOG++HmdnHcirsQmpsaKTtmb7xSW3ES4+bGoOGRhEWPIihsaJ9svHI+NPyUUi6tsd7CiexKilsasQSG+jJ0dKh2Su8C+y3iExTnZnEy7zg2mxW/oGCikycxeMRIPL37duOV86Hhp5RySXVVTRRmVXCyQBuxnK+G2hp745XcHBrranD39CIsejRh0aPxDx7k7PKcSsNPKeUy7GNDNlCY1TISi4cbQ6IGMmRUIN7aiKVTbFYrZfnHKc7JorL4BACBYUMYmZDEoPBI3Ny1cz9o+CmlXIDYhJNFtRRl2WdK9/B2J2L8IMJGBuDhpb+sz0VE7I1XcrIoOZaLtbkJ7wH+RE6IJzRqFD5+/s4u0eVo+CmlnMZqtVGWV0NhViWNtc14+3kSlTiYwRH+OhJLJzQ3NlB67CjFuVnUVpTj5ubOoJbZ0APDhugz0Q5o+CmlLjhLk73lZlFOJZZGK35B3kSmDiF46ACMttzskNhsVBQXUZyTzcmC44jNhn9wCKNSJhMyYiSeXvpMtDM0/JRSF0xTvYWinEqKj7a03AwbwLDRgQSEaMvNc6mvqW5pvJJNU30dHl7eDB09jrCoUfgFBTu7vD5Hw08p1evqa5ooPFJJWX4NAIOG+TF0tLbcPBerxUJZ3jGKc7OpKjkBGIKGDiM6OZXgYeHaeKUbNPyUUr2mpryRwqwKyovqcHOD0BEBDB0ViI+fa8/15kwiQk15GcU5WZQeO4rV0oyPfwAj4pMIHTkK7wEDnF3iRUHDTynVo05NHFt4pIKqUvuYm8PGBDI0OlDn0OtAc0MDxUdz7LOhV1Xi5u5BSEvjlYGhYXpbuIdp+CmleoTYhPKiOgqzKqitaMTTx4PI2EGEjRyoY262Q2w2yosKKc7NorwgHxEb/oMGMzp1CiGRI/Hw1Cvk3qLhp5TqFptNWrorVNBQo90VOqO+uorinGxKjmbT1FCPp7cPw8aMIyx6NAMCg5xdXr+g4aeUOi9Wi42SY9UUZVfSVG8faHr0xDAGDfPT7gptsFqaKT1+jOLcLKpLSzDGjaBhw4mOGkXwsOG4uekt4QtJw08p1SX2PnpVFGVXYWmyEhDiQ1TiYAJDdbb0M9mHayuhODeHsuOnGq8MZERCMmEjo/Hy1cYrzqLhp5TqlKYGCydyqig+WoW12UbQkAEMGxNEwKD+NyPAuTTV11FyNNfeeKW6CncPT0IiRhAWPYqAkFD9T4IL0PBTSnWosa6ZwqxKSo9XY7PZ++gNG6N99M5ks1kpLyygODebisICRGwEDA5ldMxUBkeOwN1DG6+4Eg0/pVSbzuyYPjjCn2Gjg/Dx11/ip6ttquWrjf/AVlzFIPdAvHx8GT4ultCoaAYMDHR2eaodGn5KqVZqKxspPFLJycJa3NwMYSMHMnS0Til0prrmOnYV7yKjLIPG/KOMCI5i/MRZBA8d1m9mQ+/L9F+zUgqA6pP2jukVJ+q0Y3oH6prr2Fm8k31l+7DarIwNHsvEG25g0IAQZ5emusBlw88YEwS8BcQDAtwlIl87tSilLjKnJo8tOGwfjcXDy53wmGCGRA3UefTOUNtcaw+90n3YsDEucCyJ5QEE1AbgM1KDr69x2fADXgY+E5EbjTFegLYJVqqHiAiVxfUUHKmg5mQDnt4eRE4IIWxkAO4eesvudLXNtew4sYPMssyW0BvDpGGXEOgdSOWeVTSaYnwmTHB2maqLXDL8jDGBwAzgDgARaQKanFmTUhcDEfsQZAWHK6irbMTL14ORCYMJjdTRWM50ZujFBMeQWDMIs24XAbfaW7oGzJ+P8fJycqXqfLhk+AHRQAnwJ2NMEpAOPCgitadvZIy5F7gXYMSIERe8SKX6CrEJJwtrKThcQX11E95+nkQnhRIS4Y+bjsbSypmhF+sxgqThEwkOHoaltJS6oYVIczP4+ODmrd09+iojIs6u4SzGmEnAFmCaiGw1xrwMVInIU+29Z9KkSZKWlnbBalSqL7DZhLL8GgqP2Mfd9A3wYvjYIB2CrA1nht744PGkBMVhWfk+PvHx+F9+ubNLVOfBGJMuIpPOXO6qV355QJ6IbG15/R7wX06sR6k+xWYTSo9XU3ikksa6ZgYEejMmdQjBwwbo6CJnOLMhS6x1CHHWIYSOmAZAwxVz8Awf7uQqVU9zyfATkSJjzHFjTIyIHATmAJnOrkspV2ez2ig5br/Sa6q34BfkzYi4IQQN0dA7U11zHTuKdzhCLyY4htQhqbhvz6Dx4AEkaTLGywufmHHOLlX1ApcMvxYPACtbWnpmA3c6uR6lXJbVap9hofBIJc0NFvwH6WDT7TnVT++LrDSOnazh8sBRXFXYROjceDy9A7FNnoTfJZO1IctFzmXDT0R2AWfdp1VK/cepaYUKj1TS3GghIMSXUcmhDBzso6F3hrrmOnaV7CKjNIPjZdV8vLGRyoZxpIk3qRENDG5oANBGLP2Ey4afUqp9VouN4twqCrMrsTRaGRjqy+ixoQwM8XV2aS6n3lJvH4asNAOLzcLY4LEEfrSXov1NrBkxgDoDG2MvJ0lbjPcrGn5K9SHWZhsncqsoyq7E0mQlMGwAw8fqtEJtabA0sLtkN3tK9mBtamR8lT+Jl97AIN9B7JgRxBulGbgjeHq4M3WUjtDS32j4KdUH2EOv0jGBbGDYAMLHBeEfrKF3pkZrI3tK9rC7ZDdN1iZGB40mqSIQ931pBMTUQzhMnDmJ30SNZkt2GVNHhZA6MtjZZasLTMNPKRdmabZSnPufWdODhgxg+Nhg/IP1udSZmqxN7C3dy67iXTQ11TM+D2JHT2Fo1GTEasUSNgLP4f/pspA6MlhDrx/T8FPKBVmarZzIsd/etDbbCBo6gPCxwfgFaeidqdnWTEZpBjuLd9LQVM/IoCgmj5qM56F/4Vlub8Ri3N1bBZ9SGn5KuRANvc6z2CzsK9vHzhM7qbPUMabIMK7AjZG3z8d4eCA33KDdFVS7NPyUcgFnhl7wUD+GjwvCL1BD70xWm5X9J/eTfiKdhqpyhgRFMD9qPiEhFhrMfqS5GePhocGnOqThp5QTWZttFOVUauh1gk1sHDx5kLQTaVQ3VTPcBDF9ayPDLo9mgP8w8AevyEhnl6n6CA0/pZzgzNabenuzfTaxcaTiCNuLtlNTXszQBm9mTlpAZEAkDbIbr+hoZ5eo+iANP6UuIKvF5ri9qaHXMREhpzKHbUXbONlwkhCfEK44MZjgkgZCZgzHGINvcrKzy1R9lIafUheAY0SWrErtsnAOIsKx6mNsLdxK+ckCwg9XMHH2IsZGJGEbXg3GYDz0V5fqHv0XpFQvslptFOdWU5hVgaVRO6efS35NPtsKt1FYU0CA90BmDL+cwbu3M7DBD2MM7gMHOrtEdZHQ8FOqF9isNoqPVVN42D7g9MBQX8LHBeswZO0oqi1iW9E28qrzGL63iBn+0cQu/Bbubu7Y7krATVtuqh6m4adUDzpzPr2Bg30ZnaoDTrentL6UbYXbOFZ8GG//gVw2/DKi62pwtwpuxg1Ag0/1Cg0/pXrAqZnTCw7bQ89/kE/L1EIaem2paKhg+4ntHCk/QsCJai7dXcW422/CN2w4hDm7OtUfaPgp1Q1iE0rzayg4VEFjXTP+wT5EJw1m4GCdRLYt1U3VpBWlcTh/N55WSBk9laRxE7AO2IOXf6Czy1P9iIafUudBRDhZUEv+oXIaapoZEOjNuEuGEhimodeWuuY6dhTvIKM0A0RI3VbOyIg4wqZPtW9w+XTnFqj6HQ0/pbpARKg4UUf+wXLqqprwDfBizKQhBA8doKHXhkZro30i2eNp+OQUMW7qLCYPuwTvwZXaclM5lYafUp0gIlQW15N/qJzaikZ8/L0YPTGMQcP9NPTa0GxrZm/JXnYW76TR2sj4ah/GFgUx1DcBT68AiAhwdomqn9PwU+ocqkrryTtYTs3JBrwHeBKdHMrgcH+Mm4bemaw2Kx/s28qawxtJKC5iXMwEUqd+kxDvQdhSq3APCnJ2iUoBGn5KtaumvIG8A+VUldbj5etBVOJgBkcG4Kahdxab2DhcfpgPMtfz9vYDWBoC8d1nmB4WyWDfwQAafMqlaPgpdYbaykbyD5ZTcaIOD293RsSFEDYyADd3N2eX5nJEhNyqXLYWbqXx8GECN+fT7D6FZkson4bbGBc8mhRnF6lUGzT8lGpRX9NE/sEKThbU4O7pRsT4QQyJHoi7h4ZeW/Jr8tmat5kTNUUM9BvEpdGzqC3O443DPlQbg5unJ1NHhTi7TKXapOGn+r3GumbyD1VQmleDm7th+Nggho4OxMPT3dmluaSSuhK2FG4hvyyHqC8PMmPyLMZPXIy7mzvEwx+OlrMlu4ypo0JIHRns7HKVapOGn+q3mhosFByuoORYNcbA0OiBDBsThKe3hl5bKhoq2FawhaPHMnALCWZq1Ayip6fgGzXKHnwtUkcGa+gpl6fhp/odS5OVwqxKTuRUISKERgYwfGwQXr7649CW2uZa0orSyDyZScju40w64cb4796Dr3+QDkWm+iz9aVf9hn0i2UoKsyqxWoSQcD/CxwXj4+fp7NJcUoOlgZ0ndnA4cyMN/l7ERU4k5eqFeJbX4O2nQ5Gpvk3DT130bFYbxUerKThin1MveKgf4THBDBioswW0pdnWTEZpBjtO7MBSU03S9hNEXzafIREz7BuEOrc+pXqChp+6aIlNKM2rIf9QuWN6oYjJwTqRbDtsYuPAyQPsPLAeKShiSOolTB2zkIGRjXgOHeLs8pTqURp+6qIjIpQX1pF3sJyGmib8grx1poUOiAg5lTlsKdxCRWMFo/KqGFcSSFT4lbh5e0OEsytUque5dPgZY9yBNCBfRBY4ux7l+ipL6sk7cJLaikYddLoT8mvy+Xv651Rt3s7AxBium3YdI8cPx4A9+JS6SLl0+AEPAvsBHf5ddaimvJG8gyepKrEPRRadFMrgCB1/sz2l9aVsLdzK18cO8M+tRczbbWVbnjeXxwczKkhvC6uLn8uGnzEmArgWeAZ4xMnlKBdVX9NE3oFyygtr8fDSocjOpbqpmm1F2yhI20hAcTWMnEF1xVj+MUIQdw+2ZJdpHz3VL7hs+AEvAT8G2p37xBhzL3AvwIgRIy5MVcolNNVbyD9UTslx+6gs4eOCGToqEHdPDb221Fvq2VmYzt6yDDCG5MDRjDa+RI6dynub02nGhqeHmw5HpvoNlww/Y8wCoFhE0o0xs9rbTkTeAN4AmDRpklyY6pQzWZqsFBypoDi3ChEYEj2Q4ToqS7tOzau3O2czIeszGD/9ciZech3+nv4YYxgMrLxnqg5HpvqdboefMWYk4CYiOT1QzynTgOuMMdcAPsBAY8wKEVnag8dQfYjVauNEdhWFWRVYLcLgcH/CY4LwHqAd1NtyqttC2tHN1Lg1MXLwSJJiBhM2chpeXq1vpuhwZKo/6lb4GWOeB6KB0caYicCNIvJud4sSkZ8AP2k5xizghxp8/ZPNJpQeryb/UAXNDRaChg4gImaQdlBvh4hwtOooWwq3wJadRJRZGHfPg0QEjoAxzq5OKdfR3Su/ySIywxizTkSsxpi7gW6Hn1IiwsnCWvIPltNQ04z/IB/GpIYRMEhbIranqLaIrYfXUWAtY+CAYKakLmBY0wAG+A13dmlKuZzuhl+TMcYHOPW8rcf/Oy4i64H1Pb1f5bqqSus5vv8/ffXGTh5C0BDtq9eeysZKthRu4WjePqLW7Gf67KuYkLKo1UwLSqnWuht+vwW+AEKMMbf0QD2qH6utbCTvQDmVxXX2vnrJoQwO17567alrriM9fyuHs7djHRxM8ujpjPeZjl9MrAafUufQ3fDLAO4EbgXigNu6XZHqdxrrmsk/WEFpfg0enm5ETghhSJT21WtPs62ZPSV72Fm8k8AtB0io8GTC9+7Fzy8Ihjm7OqX6hu6G31sichXwy54oRvUvzU1WCg9XcCK3CmNg2OhAho0OxMNLr1raYhMbB8sOsHvn51QMdGNk2DguufYKAmzeePkFObs8pfqU7obfHmPMeBE50CPVqH7BarVxIqeKwiMt3RYi/AmPCcZbJ5NtU1ruST4/lIGn/xGCbBXEbM7ispkLGBF9tbNLU6rP6u5vm9nA3caYNcDXQJqIbOp+WepidOYUQ0FDBhAxXrstdORfhw7z43f/j8H1+RwfNIxnr7mRqXfeiFd4uLNLU6pP61b4ichkY4wH9ud9E4FvAhp+qhURobLY3oKzvto+xdCo5FAGDvZ1dmkuq7qpmq2FW1mxbwsxJceILHLnaFQqx08E4D1Bh/JTqrs6FX7GmDnYr+oqz1wnIhZgd8uXUq3UVjRyfP9Jqkrr8fbzZHRqGIOG+Wm3hXY0WhtJL9hG7pa1NAwN4oqoKfzWfzx7RwDevjr2plI9xIice0hMY0wt9mHGsrDPr7e95c90Eanr1Qo7adKkSZKWlubsMlSLhtpm8g6Uc7KgBg9vd8LHBhM6MgA37bbQJqvNSkZxBpXFlfjYvPFqsOLh44u7tw9NFhuNFiveHu54eWgLWKXa4uPjQ0REBJ6erYc8NMaki8ikM7fv7G3PgdhvbU4CJgPfAv4HcDfGHAS2i8id3apcXRSam6wUHLYPPG2MYfjYIIaNDtLZFtohImRXZpO+93OGE8XIMXEMCx2Gp3HHuGurV6U6Q0QoKysjLy+P6OjoTr2nU+EnIlZgT8vXHwGMMZ7AIuBXwO3Y+/upfsrW0oKzoKUFZ2ikvQWnl4+24GxPUW0RXxd8TWFtIRGl1QwMH0RkaDhuHvqZKdUVxhhCQkIoKSnp9Hu6/FNmjAnH3rDlZiAJ+BL4RVf3oy4OIkJZfi15B07SVG8hMGwAkbHagrMjlY2VbD22iZObv0IihjEz+UrGx43l0OEjGnxKnaeutiPobIOXIcBNwC3Yb32uA/4X+IeIVHStRHWxqCqt51jmSeoqGxkQqC04z6XR2kh6UTp7SvfgZoPJ1QGMCphKYEicfQNtBKTUBdPZBzF5wHex3/IcJiJXi8ifNPj6p/rqJg5tK+LA14VYmqyMSgkj7vLhGnztsNqs7CnZw3v/epWCj/7O2MDRfCtuKcnfeYzAS6Y6uzyHsrIykpOTSU5OZujQoYSHhzteNzU1tdr2jjvu4L333jtrH+vXr2fBggVdOu5VV11FUFDQWe/LyclhypQpjBkzhptvvvmsGk7VPHv2bPz9/fn+97/f7jHuueceMjMzO13T8uXLCQ0NJTk5mQkTJvDmm292/oQuEH9/f2eX4JCbm0t8fPw5l7/55pukpqZSXl7e7r+h0y1fvpyCggLH65deeom6up5pY9nZ8CsBYoFngZXGmP82xnzDGKMjCfYjTQ0WcvaUsverfKpPNhAZO4jE2REMjvDXrgttONWY5a8H/8rG/I2EeAxkyoBYZoVMxd/LH+PpWhPxhoSEsGvXLnbt2sV9993Hww8/7Hjt5dV7t7F/9KMf8fbbb5+1/LHHHuPhhx/myJEjBAcH84c//OGsbXx8fPjFL37Bb3/72w6P8dZbbzFhwoQu1XXzzTeza9cu1q9fz+OPP86JEydarbdYLF3anytxRu1vv/02r776Kp9//jnBwZ2bPNnp4Sciw4EI4NvADmAKsBzIM8bkGWM+7JFqlEuyWm3kHypnz7o8So9XMyRqIIlXRDJsTJAOPt2O4rpiPjzwPml/ew2/7BNcO+pa5s26h4jb78E9IODcO+ik9KPlvLbuCOlHy3tsn6d78803mTx5MklJSdxwww2tfvGsXbuWSZMmMW7cOD7++OOz3ltbW8tdd93FJZdcQkpKCh999FGbx5gzZw4BZ3wmIsKXX37JjTfeCMCyZcv48MMPz3qvn58f06dPx8en43keZ82axamuUP7+/jzxxBMkJSUxderUs0LtTGFhYYwePZqjR49yxx13cN999zFlyhR+/OMf87Of/axV8MbHx5Obm0tubi7jx49nyZIlxMbGcuONNzo+u6ioKH784x+TkJDAJZdcwpEjRwD45z//yZQpU0hJSWHu3LmOumpqarjzzjtJSEggMTGR999/33G8ts6jpKSEG264gcmTJzN58mQ2bbKPO/Kzn/2M2267jWnTpnHbba3nIKipqWHOnDlMnDiRhIQEx99Vbm4usbGxfPvb3yYuLo558+ZRX18PQHp6OklJSSQlJfHaa691+Bn+7W9/49lnn+WLL75g8ODBZ61PT09n5syZpKamMn/+fAoLC3nvvfdIS0tjyZIlJCcn8/LLL1NQUMDs2bOZPXt2h8frjE7/5hKRAhH5SESeFJH5IhICjAN+hL3/n7oINdQ2s+fLPPIPlhMY6kvCzAhGxofgqYNPt6m6qZp/5a7lvUPvUWGpIs5nFHOHXM7IgSNxc3PDuPXcfxbSj5az5K0tPP/FQZa8taVXAnDx4sVs376d3bt3Exsb2+rqKzc3l23btrF69Wruu+8+GhoaWr33mWee4YorrmDbtm2sW7eOH/3oR9TW1nbquGVlZQQFBeHR0gAoIiKC/Pz8Hjmn2tpapk6dyu7du5kxY8Y5b2lmZ2eTnZ3NmDFjAMjLy2Pz5s288MILHb7v4MGDfO9732P//v0MHDiQ3//+9451gYGB7N27l+9///s89NBDAEyfPp0tW7awc+dObrnlFn7zm98A8Itf/MKx/Z49e7jiiis6PI8HH3yQhx9+mO3bt/P+++9zzz33OI6bmZnJ2rVreeedd1rV6uPjwz/+8Q927NjBunXrePTRRznVB/zw4cPcf//97Nu3j6CgIEf43nnnnbz66qvs3t3x+CZHjx7l+9//Pl988QVDhw49a31zczMPPPAA7733Hunp6dx111088cQT3HjjjUyaNImVK1eya9cuHnzwQYYPH866detYt25dh8fsjO4Ob5aFPfjeOde2qm/y9vUgMNSX0BEBOot6B5qsTews3smBPesJ3nuMlBtvJjVyKp7xnr12S3hLdhlNFhs2gWaLjS3ZZaSO7NztpM7KyMjgySefpKKigpqaGubPn+9Y981vfhM3NzfGjh3LqFGjOHCg9fj2X3zxBatWrXJcGTU0NHDs2DFiY2N7tMau8vLycjxfTE1NZc2aNW1u9+6777Jx40a8vb353//9XwYNGgTATTfdhHsn+mBGRkYybdo0AJYuXcorr7zCD3/4QwBuvfVWx58PP/wwYA/Vm2++mcLCQpqamhz91dauXctf//pXx35P3TJs7zzWrl3b6vlmVVUVNTU1AFx33XX4+p79bF5EePzxx/n3v/+Nm5sb+fn5jivJ6OhokpOTHcfJzc2loqKCiooKZsyYAcBtt93Gp59+2ubnEBoayqBBg/jb3/7mONfTHTx4kIyMDK688koArFYrw4b1/hM1bVetOmTcDKOSQ51dhsuyiY0DJw+wrWArddZ6YgaPYsLwYYQFJeLu3rvdPaaOCsHLw41miw1PD7deGfrsjjvu4MMPPyQpKYnly5ezfv16x7ozQ/3M1yLC+++/T0xMTJePGxISQkVFBRaLBQ8PD/Ly8gjvocG8PT3/8x8Sd3f3dp9/3Xzzzfzud787a7mfn5/jew8PD2w2m+P16Ve/HX0+bX3/wAMP8Mgjj3Ddddexfv16fvazn53XedhsNrZs2dLmreDTaz/dypUrKSkpIT09HU9PT6Kiohzn4u3t7djO3d3dcduzswYMGMAnn3zC5ZdfTlhYGEuWLGm1XkSIi4vj66+/7tJ+u6tL92CMMSPa+RrSWwUq5aqOVx/n7/vfZd/f32LowVIWj13MnNSbGHrrbbgHBfX68VNHBrPynqk8Mi+GlfdM7fGrPoDq6mqGDRtGc3MzK1eubLXu73//OzabjaysLLKzs88Kufnz5/Pqq686bp/t3Lmz08c1xjB79mxHa8A///nPLFy4sJtn0/OioqLYsWMHADt27CAnJ8ex7tixY45f6H/5y1+YPn26Y927777r+PPSSy8FoLKy0hHwf/7znx3bXnnlla2eqZWXd3x7e968ebz66quO17t27TrneVRWVhIWFoanpyfr1q3j6NGjHW4fFBREUFAQGzduBDjr38aZwsLC+Oyzz3j88cf5/PPPW62LiYmhpKTE8Vk1Nzezb98+AAICAqiurnZse+br7ujqA4hcIKeNrwJjTJ0x5o/GmIE9UplSLqq8oZxPjnzMP7P+SRMWEsMncXnUbIb62Z9nXMiWr6kjg7l/9pheCT6wP2+aMmUK06ZNY/z48a3WjRgxgksuuYSrr76a119//awrjaeeeorm5mYSExOJi4vjqaeeavMYl19+OTfddBP/+te/iIiIcPxy/PWvf80LL7zAmDFjKCsr4+677wZg1apVPP300473R0VF8cgjj7B8+XIiIiK61KWhu2644QZOnjxJXFwcv/vd7xg3bpxjXUxMDK+99hqxsbGUl5fz3e9+17GuvLycxMREXn75ZV588UXA3iDlpptuIjU1tVWjkCeffJLy8nLi4+NJSko65/OuV155hbS0NBITE5kwYQKvv/76Oc9jyZIlpKWlkZCQwP/93/+d9Xfdlj/96U/cf//9JCcn05kxoqOjo1m1ahV33XUX27Ztcyz38vLivffe47HHHiMpKYnk5GQ2b94M4GhglJycTH19Pffeey9XXXVVjzR46dTA1o6NjbkLuA37iC5HgZHAE8DfsAfjL4BdInJvtyvrIh3YWvW2eks9aUVpHNm/mbD0owy78VYSR12Kh1vPPD3Yv3+/05+HqZ6Rm5vLggULyMjIOGtdVFQUaWlpbbZ6VN3T1s9Qdwe2PuXHwHQRKW15nWWM2QtsEJHxxpjDwL/Pp2ilXJXVZiWjLIPtBdtoFgtxkUmMr4xi0KC4Hgs+pdSF1dWf3KHAmU8761uWIyLZxpjAnihMKWcTEXKrctmcvwnPTTsZ5R1C4k33Mth3MJz7rpDqx6Kiotq86gP7VaFyvq6G3wZguTHmh8BxYAT2UV82ABhjEoCiHq1QKScorS9l8/GN5NUVEOQdxKRxVxDmMxg/H51MVqmLQVfD7x7gL9gbuZx6WLge+/x+YG9A8+0eqUwpJ6hrrmNb0TaystIYtjmL6Qu/SVzMDNxjtVO/UheTLoWfiJwA5rRMaxQO5ItI/mnrO+7qr5SLstgsvJexmbVHNjE8xId5UZMYXzmawMFxuLtp8Cl1sTnfp/UWwNbyp1J9loiQU5nDu/vWsvuDDQwq9+BvEdew8O4kBi/sne4DSinn61L4GWOCgbeBa1oWiTHmE2CZiJzs6eKU6k2l9aVsPPZvCuoKOVZmo7Q5kSbcsTZ698pQYUop19HVTu4vtvw5HvDEPs2RAB2P8KqUC6lrrmPdsXV8sHMFXqu+ZEZzNHclfou8wHh2D4nB3cuzV4YKU0q5jq6G3zxgiYgcEhGriBwClgHzz/E+pZzOYrOws3gnf9m3ggPlB4iLTGVK3HzGDUtgclRIrw8V1lc888wzxMXFkZiYSHJyMlu3bqWioqLVrATtaW9S0+545ZVXiI2NdYwJ2duTuD766KMkJSXxwAMPUF9fz8yZM7Fare1u39TUxIwZM1xifr/ufjbHjx9n9uzZTJgwgbi4OF5++eVW661WKykpKV2esNgVdfWZn+E/rTxPsbUs7zHGmEjg/4AhLcd7Q0Re7vhdSrVNRMipyuHrgq+xZR4i5ngT8Xc/SnDAYIj8z3apI4P7degBfP3113z88cfs2LEDb29vSktLaWpqcoTf9773vQte0+9//3vWrl1LRERErx8rKyuLTZs2Oabpee2111i8eHGHszh4eXkxZ84c3n333bMGbe5rPDw8eP7555k4cSLV1dWkpqZy5ZVXOiYCfvnll4mNjaWqqsrJlXZfV6/81gBvG2NGGWPcjDGjsE9q+0UP12UBHhWRCcBU4H5jTNemYVYKKKsv459Zq/gs61PcjBsz4haQFHcFgZ69e/XQVxUWFjJ48GDHSP6DBw9m+PDh/Nd//RdZWVkkJyfzox/9CIAXXniB+Ph44uPjeemll87aV3Z2NikpKWzfvh2AFStWcMkll5CcnMx3vvMdrFYrtbW1XHvttSQlJREfH+8Y8PmU++67j+zsbK6++mrHGJina6uG5557jldeeQWAhx9+2DH/3ZdfftlhOB08eJBZs2Zx9OhRUlJSqK2tZeXKla0G1J49e7Zj6qAnn3ySBx54AIBFixadc3DnjrS33+5o7+/nF7/4BTExMUyfPp1bb7211WS8w4YNY+LEiYB9EOnY2FjHPIp5eXmsXr261fyAfVlXr/wewt7P7wj2KzKDPfh6tG+fiBQChS3fVxtj9mPvWnHhRqxVfVq9pZ7tRdvJLNpD+KYjXB53KROSb7R3W+gLo7McXgs1Hc8w3mX+Q2Ds3A43mTdvHv/93//NuHHjmDt3LjfffDMzZ87k2WefJSMjwzFDQHp6On/605/YunUrIsKUKVOYOXOmY665gwcPcsstt7B8+XKSkpLYv38/7777Lps2bcLT05Pvfe97rFy5Ej8/P4YPH87q1asB++wCp3v99df57LPPWLdu3VljYbZXw+WXX87zzz/PD37wA9LS0mhsbKS5uZkNGzY45p+75ppreOuttxg+fLhjfzExMSxbtoyoqCjuuecempqayM7OJioqyrHNz3/+c55++mmKi4vZuXMnq1atAuwzuJ8K+dNdfvnlbc5C8Nvf/pa5c//zd9Hefs9Xe5+NxWLh/fffZ/fu3TQ3NzNx4kRSU1Pb3Edubi47d+5kypQpADz00EP85je/6bFZFZytq/38TgJXGWOGYb9hdLwlqHqNMSYKSAG2trHuXuBesI8wr5TVZmVf2T62F2yjiWbihiYSFzcGv4go7a/XCf7+/qSnp7NhwwbWrVvHzTffzLPPPsusWbNabbdx40auv/56x/xwixcvZsOGDVx33XWUlJSwcOFCPvjgA8ftsn/961+kp6czefJkAOrr6wkLC+Nb3/oWjz76KI899hgLFizg8ssv73St7dXw3e9+l/T0dKqqqvD29mbixImkpaWxYcMGxxXhJ5980uY+9+7d67jSKy0tJeiMqalmzJiBiPDCCy+wfv16x+1Qd3d3vLy8qK6uJiAgwLH9hg0bOnUu7e33lLlz51JUdPbgWc8880ybUz2199nYbDYWLlyIj48PPj4+fOMb32iznpqaGm644QZeeuklBg4cyMcff0xYWBipqamt5nTsy84ZfsaYx8+xHgAR+VUP1XT6vv2B94GHROSsm8wi8gbwBthndejp46u+5XjVcTYWbKQhN4exmRXELnuA0JBI6P1HRT3vHFdovcnd3Z1Zs2Yxa9YsEhIS+POf/3xW+HUkMDCQESNGsHHjRkf4iQjLli3jf/7nf87afseOHXzyySc8+eSTzJkzp9V0RefD09OT6Oholi9fzmWXXUZiYiLr1q3jyJEj55w1Y9++fY4GO76+vq0mpwV7OBYWFhISEtIq5AAaGxvPmtaps1d+He0X7LOzXyjNzc3ccMMNLFmyhMWLFwOwadMmVq1axSeffEJDQwNVVVUsXbqUFStWXLC6epyIdPgFrOvE15fn2k9Xv7B3pfgceKQz26emporqnyoaKmR11mp5bcfvZEXmCsk6ukvK//GhWMrLnV1al2RmZjq7BDlw4IAcOnTI8fqJJ56Q+++/X0pLS2XEiBGO5enp6ZKQkCC1tbVSU1MjcXFxsmPHDsnJyZG4uDipqamRadOmycqVK0VEZN++fTJmzBg5ceKEiIiUlZVJbm6u5OfnS319vYiI/POf/5SFCxeeVdPIkSOlpKTE8drPz6/DGkREfvrTn0pkZKSsWbNGioqKJDIyUhYtWtThuVdVVcn48eNbLYuIiHDUV1BQIAkJCZKZmSlz586VTz/91LFdaWmpxMTEdPzhtqOj/XbVuT6bbdu2SUpKitTX10t1dbWMHTtWnnvuOcf7bTab3HbbbfLggw+2e4x169bJtddee9419qa2foaANGkjM8555Sci3Z81sIuM/XLyD8B+EdE+hKpNTdYm0k6ksad4N6E7jjIteAxxC2+2TzM0IsnZ5fVJNTU1PPDAA1RUVODh4cGYMWN44403CAkJYdq0acTHx3P11Vfz3HPPcccdd3DJJZcAcM8995CSkuKYscDPz4+PP/6YK6+8En9/f6677jp++ctfMm/ePGw2G56enrz22mtUVlbyox/9CDc3Nzw9Pfl//+//dbrWiRMntlkD2K+4nnnmGS699FL8/Pzw8fFpdUu1rWd+GRkZZ3XTmDdvHhs3buSyyy5j8eLFPP/888TGxvLUU0/x2GOPcdVVVwGwbt06rr322i5/3nV1dR3u93x19Nlcd911JCYmMmTIEBISEggM/M9EPJs2beLtt98mISGB5ORkAH71q19xzTXXnHWMvq5Lk9leKMaY6dhnitiLvSsFwOMi0vaNenQy2/5ERDhYfpAt+V9TZ60nZlAMybkGH09fBkydekFnUu9JOpmt69mxYwcvvvgib7/9dofbLV68mGeffbbVTO6uqqamBn9/f+rq6pgxYwZvvPGGo4VnX9ebk9leECKykR7uO6guDkW1RWzM30h50VFG7SpmwsJlDBsxwT65llI9bOLEicyePRur1dpuX7+mpiYWLVrUJ4IP4N577yUzM5OGhgaWLVt20QRfV7lk+Cl1ptrmWrYUbuHgyYMM8BjA5WOvZEjhfvzcde5k1bvuuuuuDtd7eXlx++23X6Bquu8vf/mLs0twCRp+yqVZbVb2lO4hrSgN36wCLqkLIuHmW/H28Ea+mdBnb3EqpZxLw0+5rKNVR9mYv5HKxkpGDhzJ5IgJeBeW4WmzB54Gn1LqfGn4KZdT0VDBpoJNHC/NInJfKdMmXUXUqMuQaHvjLA09pVR3afgpl9FkbSL9RDq7S3bj4ebB1MjpRB45wACbfZQKDT2lVE/R8FNOJyIcKj/ElsItNJ84QUKBIfnG+/DzCUBuSca4dXX8daWU6piGn3KqkroSNuZvpLC2kLABYVw6dDY+R/fhXdcMPmjwKaV6hYafcop6Sz3bCreRWbqPkGOVzApLIDbpWowxyNhkjKens0tUSl3ENPzUBWUTG5llmWwt3EqTrYmE0ATG55zAs8w4nulp8CmlepuGn7pgCmoK2Ji/kbKqE4w52kjKnG8xOHg4ttBGjJeXs8tTqkd9+OGHrF69mqqqKu6++27mzZvn7JLUafSBiup1tc21rD26lg+PfEiDpYG5g6aSfGIA/iU1ALh5e2tLThfz6KOPkpSUxLe//W1mzpyJ1Wptd9umpiZmzJiBxWK5YPV99tlnxMTEMGbMGJ599tkub9fR+1988UXi4uKIj4/n1ltvdUxr1N7y9ixatIg333yT119//awZ6nv6XNurLSoqyjFI9aRJZw1v2W3Hjx9n9uzZTJgwgbi4OF5++eV2t7VaraSkpLBgwYJWy19++WXi4+OJi4trNeN8r2trqoe++KVTGrkei9UiO07skDd2vyFvbnhBtn31rjRZmuzrqmucXJ3rcYUpjUREjhw5IlOmTBERkd/97nfy0ksvnfM9P/vZz2TFihW9XZqIiFgsFhk1apRkZWVJY2OjJCYmyr59+zq9XUfvz8vLk6ioKKmrqxMRkZtuukn+9Kc/tbu8Mx555BFJT09vc926detk2bJl3TrXjmo7czqonlZQUOA4t6qqKhk7dmybfxciIs8//7zceuutraZD2rt3r8TFxUltba00NzfLnDlz5PDhw+ddT1emNNIrP9Urjlcd592D7/J1wdeE+4fzjbpxjDpUhbvV3lHd3d/PyRWqthw8eJBZs2Zx9OhRUlJSeOutt1rNFD579mzWrFkDwJNPPskDDzwA2K9yVq5ceUFq3LZtG2PGjGHUqFF4eXlxyy238NFHH3V6u3O932KxUF9fj8Vioa6uzjHtUXvL2/tMRITHHnuMq6+++rwHj+7subZXW0faq7srhg0b5ji3gIAAYmNjyc/PP2u7vLw8Vq9ezT333NNq+f79+5kyZQoDBgzAw8ODmTNn8sEHH3S5jvOhz/xUj6pqqmJz/mayK7MJPWnl6jGziB6RgG14A2Kx4qbP9jplY/5GSutLe3Sfg30HMz18eofbxMTEsGzZMqKiorj99tsZMWIEUVFRjvU///nPefrppykuLmbnzp2sWrUKgPj4eLZv337W/jo7k3lX5OfnExkZ6XgdERHB1q1bO71dR+8PDw/nhz/8ISNGjMDX15d58+Y5ntW1t7y9z+TVV19l7dq1VFZWcuTIEe67775eOdeOajbGMG/ePIwxfOc73+Hee+91vK+9us9Xbm4uO3fuZMqUKWete+ihh/jNb35z1r+F+Ph4nnjiCcrKyvD19eWTTz7plduzbdHwUz3CYrOwq3gXO4p3ADBlUAqRW3bjSwmMADcfHydXqDpr7969LFy4kNLSUoKCglqtmzFjBiLCCy+8wPr16x3T/Li7u+Pl5UV1dTUBAQGO7Tds2NClY8+dO5eioqKzlj/zzDOtrkB7S3l5OR999BE5OTkEBQVx0003sWLFCq699to2ly9durTdz+QHP/gBP/jBD9o8zpQpU2hsbKSmpoaTJ086Jo799a9/zfz583uk5qVLl7Jx40bCw8MpLi7myiuvZPz48cyYMQNo/+/ylK78XdTU1HDDDTfw0ksvMXDgwFbrPv74Y8LCwkhNTWX9+vWt1sXGxvLYY48xb948/Pz8SE5ObnfqqJ6m4ae6Lbcyl00Fm6isr2Bc/UCmTF5EgFcAzdePwiMkxNnl9UnnukLrTfv27SM+Pp6GhoazGnXs3buXwsJCQkJCWoUcQGNjIz5n/Cenq1d+a9euPWd94eHhHD9+3PE6Ly+P8PDwTm/X0fvXrl1LdHQ0oaGhgH2S2s2bN+Pt7d3m8qVLl3b4mbTn1NXb+vXrWb58OcuXLz/vc22v5qVLlzq2DQsL4/rrr2fbtm2O8DtX3Z35uwBobm7mhhtuYMmSJSxevPis9Zs2bWLVqlV88sknNDQ0UFVVxdKlS1mxYgUAd999N3fffTcAjz/+OBEREZ06bre19SCwL35pg5cLr6KhQj7O+lhe2/marMxcKTnb/iXFr7wqTfn5zi6tT3KFBi9VVVUyfvx4x+uIiAipr68XEXvjhoSEBMnMzJS5c+fKp59+6tiutLRUYmJiLkiNzc3NEh0dLdnZ2Y5GIBkZGZ3erqP3b9myRSZMmCC1tbVis9nk9ttvl1deeaXd5R19Jp1xrgYvnTnX9mqrqamRqqoqERGpqamRSy+91FFfd+s+xWazyW233SYPPvhgp8/39AYvIiInTpwQEZGjR49KTEyMlJeXn1ctIl1r8OL00OqpLw2/C6fJ2iRbC7bK67telzfT/5/sOPSVWKwWsVks0nDkiNhsNmeX2Ce5Qvht3rxZbrzxRsfru+66S9asWSO1tbUydepU+eKLL0RE5KuvvpKpU6c6tvv73/8ujzzyyAWrc/Xq1TJ27FgZNWqU/PKXv2y17uqrr5b8lv+AtbddR+9/+umnJSYmRuLi4mTp0qXS0NDQ5vKKiooOP5POOFf4dVTr6efZVs1ZWVmSmJgoiYmJMmHCBMd7z/V32RUbNmwQQBISEiQpKUmSkpJk9erVZ9V3+vmeGX7Tp0+X2NhYSUxMlLVr155XHad0JfyMfV3fN2nSJElLS3N2GRc1ESG3KpdN+ZuoaqpiTNAYEreV4mWB4Ftv1b563bR//35iY2OdXUYrO3bs4MUXX+Ttt9/ucLvFixfz7LPPMm7cuAtUmVJna+tnyBiTLiJntaLRZ36qUyobK9mYv5GjVUcJkQF8Y9QCIgNH0OSRj3F30+C7SE2cOJHZs2djtVrbbYjQ1NTEokWLNPhUn6LhpzrUbGlg594V7LTV4ebpy2V+CQxfm0GARyUkgFfE2Q0N1MXlrrvu6nC9l5cXt99++wWqRqmeoeGn2tdQRcPev7Erfx2jhl/GZYnfYoDHAOrKvPAaMcLZ1Sml1HnT8FNtK8+FzI8IsDZzQ8B1kFGHb6I7xtPgN/XsTqxKKdWXaPip1kTg6GbI3QADQiD5egLrod49A/S5nlLqIqHhp/6juR72/xPKsmDIBBh3NXh44eEHAbNmObs6pZTqMRp+yq6qAPZ9CE01MHYehE/UKz2l1EVLw6+/E4GCHXDkX+DlBylLYeC5R4RXSqm+TKc06s8sTbB/FRz6AoKjIPVODT7VJcuXL6egoMDx+p577iEzMxOAX/3qV622dXd3Jzk5mfj4eL7xjW9QUVHRpWPNmjWLUwNZnLlvpbpKw6+/qi2FHX+G4v0QPQMSbgKvAc6uSvUxZ4bfW2+9xYQJE4CzA8rX15ddu3aRkZHBoEGDeO211877uO2Fn4hgs9nOe7+q/3DZ8DPGXGWMOWiMOWKM+S9n13NRObEP0pdDcx0k3gxR0/T5ngLsc7KNHz+eJUuWEBsby4033khdXR3p6enMnDmT1NRU5s+fT2FhIe+99x5paWksWbKE5ORk6uvrHVdn//Vf/0V9fT3JycksWbLkrONceumljklPt23bxqWXXkpKSgqXXXYZBw8eBKC+vp5bbrmF2NhYrr/+eurr6wHO2ndubi4xMTHcfvvtxMfHc/z4cb773e8yadIk4uLi+OlPf+o4blRUFD/96U+ZOHEiCQkJHDhwAICvvvqK5ORkkpOTSUlJaXMmCnWRaWvAT2d/Ae5AFjAK8AJ2AxM6eo8ObN0JlmaRg5+LfPkrkfT/E6mvdHZF6jRnDspb/v4HUt+yzGax2F8fOGB/3dRkf33woIiIWBsapPz9D6ThyBH767o6++usbPvrmppO1ZCTkyOAbNy4UURE7rzzTvnNb34jl156qRQXF4uIyF//+le58847RURk5syZsn37dsf7T3/t5+fXat+nXlssFrnxxhsdMwlUVlZKc3OziIisWbNGFi9eLCIizz//vOM4u3fvFnd39zb3nZOTI8YY+frrrx3LysrKHMeaOXOm7N69W0RERo4cKa+88oqIiLz22mty9913i4jIggULHOdcXV3tqEf1LV0Z2NpVG7xcAhwRkWwAY8xfgYVAplOr6svqKyDzQ6gqhMjJMGo2uF2YSSNV3xIZGcm0adMAWLp0Kb/61a/IyMjgyiuvBMBqtTJs2LAu7/fU1Vp+fj6xsbGO/VVWVrJs2TIOHz6MMYbm5mYA/v3vfzsmg01MTCQxMbHdfY8cOZKpU6c6Xv/tb3/jjTfewGKxUFhYSGZmpuP9p+acS01N5YMPPgBg2rRpPPLII4456S7YnHLKaVw1/MKB46e9zgN0WJHzVZZlb9giNohfDKExzq5IdULQ4usd3xt399avPT1bvXbz9m792te39Ws/v04f98xBygMCAoiLi+Prr7/uUv1nOvXMr66ujvnz5/Paa6/xgx/8gKeeeorZs2fzj3/8g9zcXGadR59Sv9POLycnh9/+9rds376d4OBg7rjjjlaT8np7ewP2BjgWiwWw30q99tpr+eSTT5g2bRqff/4548eP79b5Ktfmss/8OsMYc68xJs0Yk1ZSUuLsclyPzQbZX8Gev4H3QHtrTg0+dQ7Hjh1zBN1f/vIXpk6dSklJiWNZc3Mz+/btA+zB2N7zMU9PT8dV3OkGDBjAK6+8wvPPP4/FYqGystIx4/jpM5rPmDGDv/zlLwBkZGSwZ8+ec+4boKqqCj8/PwIDAzlx4gSffvrpOc85KyuLhIQEHnvsMSZPnux4FqguXq4afvlA5GmvI1qWtSIib4jIJBGZFBoaesGK6xOaamHPu/ahyoYlwsTbYcAgZ1el+oCYmBhee+01YmNjKS8v54EHHuC9997jscceIykpieTkZDZv3gzAHXfcwX333edo8HK6e++9l8TExDYbvKSkpJCYmMg777zDj3/8Y37yk5+QkpLiuBID+O53v0tNTQ2xsbE8/fTTpKamdmrfSUlJpKSkMH78eL71rW85buF25KWXXiI+Pp7ExEQ8PT25+uqrO/15qb7JJSezNcZ4AIeAOdhDbzvwLRHZ1957dDLb01Qctz/fa26AsVfC8GRnV6Q6wRUms83NzWXBggVkZGQ4tQ6lzkefn8xWRCzGmO8Dn2Nv+fnHjoJPtRCBvO2QtQ58AmHiNyFgiLOrUkopl+OS4QcgIp8Anzi7jj6juQEOroaSQxA6DmKuBU8fZ1el+pioqCi96lP9gsuGn+qC6hOw7x/QUAmjr4DIS7TTulJKdUDDr68r3AOHPrdf5SV/C4Iiz/0epZTq5zT8+iprMxxeA4W7IXgkxF4H3v7OrkoppfoEDb++qO6k/TZnTTGMvAyiLgc3V+21opRSrkfDr68pOQgHPgbjZp+JYfAYZ1eklFJ9joZfX2GzQvZ6OL4NAoZC3PXgG+TsqpRSqk/S8OsLGqsh8yN75/XwiTB6DrjrX51SrubDDz9k9erVVFVVcffddzNv3jxnl6TaoQ+KXF15LqT9EaqLYMJ1MG6+Bp/qEz777DNiYmIYM2YMzz77bJvbvPjii8TFxREfH8+tt97aagDqC1lHR9vcddddhIWFER8ff85jLVq0iDfffJPXX3+dd999t81t6uvrmTlzJlar1bHsww8/xBhz1piixhiWLl3qeG2xWAgNDWXBggXnrL2tupuampgxY0arYeT6Kw0/VyViH5dz91/BwxdS74Ahcc6uSqlOsVqt3H///Xz66adkZmbyzjvvkJnZekay/Px8XnnlFdLS0sjIyMBqtfLXv/71gtdxrm3uuOMOPvvssy4d95e//CX3339/m+v++Mc/snjxYtzd/zOl2DvvvMP06dN55513Wm3r5+dHRkaGY9zUNWvWOAYBP1ftbdXt5eXFnDlz2g3m/kTDzxU118Pe9+wzMoSOtwef32BnV6X6gYyMDC677DLH6x07djBnzpwu72fbtm2MGTOGUaNG4eXlxS233MJHH3101nYWi4X6+nosFgt1dXUMHz7csW727NmsWbMGgCeffJIHHnigV+o41zYzZsxg0KCzB4Vvqz4R4bHHHuPqq69m4sSJbda0cuVKFi5c6HhdU1PDxo0b+cMf/tBm+F9zzTWsXr0asIfkrbfe2qna26t70aJFrFy5st3PrL/Q+2eupqrQ3o2hqQbGzrM/49PRWvqdo/vKqKts6tF9Dgj0YmRcSIfbTJgwgezsbKxWK+7u7jzyyCO88MILrba5/PLL25zG6Le//S1z584F7Fd1kZH/GXAhIiKCrVu3tto+PDycH/7wh4wYMQJfX1/mzZvX6hnZz3/+c55++mmKi4vZuXMnq1at6vI5d6aOzmzTlrbqe/XVV1m7di2VlZUcOXKE++67r9V7mpqayM7OJioqyrHso48+4qqrrmLcuHGEhISQnp7eagaLW265hf/+7/9mwYIF7Nmzh7vuuosNGzacd+3x8fFs3779nOd3sdPwcxUiULADjvwLvPwgZSkMHH7u9ynVg9zc3IiLi2Pfvn0cPnyYkSNHnnUFc+oXb3eVl5fz0UcfkZOTQ1BQEDfddBMrVqxwPOOaMWMGIsILL7zA+vXrW90mnDt3LkVFRWft85lnnml1VdWb2qrvBz/4gWP2+baUlpYSFBTUatk777zDgw8+CNiD7p133mkVfomJieTm5vLOO+9wzTXXdLtud3d3vLy8qK6uJiAgoNv766s0/FyBpQkOfQonMiFkNIxfAF4DnF2VcqJzXaH1pqlTp7Jp0yZ+//vft/msqzNXfuHh4Rw/ftyxLi8vr9WzKoC1a9cSHR3Nqbk4Fy9ezObNmx3ht3fvXgoLCwkJCTnrl/TatWs7dS6dqaMz27Slo/ra4+vr26pRz8mTJ/nyyy/Zu3cvxhisVivGGJ577jnMaXd8rrvuOn74wx+yfv16ysrKul17Y2MjPj79fOB7EbkovlJTU6VPqikV2fqGyLr/EcnZKGKzObsi5SSZmZnOLkFERFatWiWDBg2Sp5566rz30dzcLNHR0ZKdnS2NjY2SmJgoGRkZrbbZsmWLTJgwQWpra8Vms8ntt98ur7zyioiIFBQUSEJCgmRmZsrcuXPl008/7bU6OrNNTk6OxMXFOV53p76IiAipr68XEZH//d//lXvvvbfV+hkzZshXX30lIiJ+fn4iInL8+HF5+eWXRURk3bp1cu2113aq9jPrFhEpLS2VmJiYTtfbl7T1MwSkSRuZ4fTQ6qmvPhl+RftEvnpOZONLImXZzq5GOZmrhN+hQ4dk2LBhUlNT0639rF69WsaOHSujRo2SX/7yl47lV199teTn54uIyNNPPy0xMTESFxcnS5culYaGBqmtrZWpU6fKF198ISIiX331lUydOrVX62hvGxGRW265RYYOHSoeHh4SHh4uv/vd77pV31133SVr1qwREZFZs2adFZwvv/yy3HfffSLyn/A73enh11HtZ9b91ltviYjI3//+d3nkkUc6XW9f0pXwc8mZ3M9Hn5rJ3WqBrC8hPx0CI2DCQvAZ6OyqlJO5wkzuAN///veZPHkyy5Ytc3YpF6UdO3bw4osv8vbbbzvl+IsXL+bZZ59l3LhxTjl+b+rKTO7a1eFCq6+AXSvswRc52T4NkQafcgFZWVmMHz+e+vp6Db5eNHHiRGbPnt2qk/uF0tTUxKJFiy7K4OsqbfByIZVlwf5VIDaIXwyhMc6uSCmH0aNHnzXCiOodd911l1OO6+Xlxe233+6UY7saDb8LwWaDoxvtI7b4DYa4xTDg7M6nSimlLgwNv97WVAuZq+xjdA5LtHdcd/d0dlVKKdWvafj1porjkPkhNDfA+GtgWJKzK1JKKYWGX+8QgbztkLUOfAJh4jchYIizq1JKKdVCw6+nWRrhwGr7jOuDx9pHa/Hs5yMpKKWUi9Hw60nVJ+yDUjdUwugrIPISHZRaKaVckIZfTyncA4c+t1/lJd8KQSOcXZFSSql2aPh1l7UZDq+Bwt0QPBJirwNvf2dXpZRSqgMaft1Rd9J+m7OmGEZeClEzwE0HzVFKKVenv6nPV8khSF8OjVWQ+E0YNUuDT100Hn30UZKSkvj2t7/NzJkzOxyKq6mpiRkzZmCxWM77eJ999hkxMTGMGTOGZ599tt3toqKiSEhIIDk5mUmT/jNcY319/TnrPB/+/t2/i9PZc7NaraSkpLBgwYJWy9s65574zPs7/W3dVTarfcLZjPfBNxhS77TPwafURSIrK4tNmzaxe/dukpOTWbx4cauJZM/k5eXFnDlzePfdd8/reFarlfvvv59PP/2UzMxM3nnnHTIzM9vdft26dezatYvTB7L/4x//eM4627J+/XruuOOO86q7M7pybi+//HK7A5ufec7d/cyV3vbsmsZq2PchVOZB+EQYPQfc9SNUPS9nVzq1FeU9uk+/oGCik1M73ObgwYPMnTsXi8VCSkoKAP/4xz8c62fPns3jjz/OlVdeyZNPPkllZSWvvvoqixYt4ic/+QlLlizpcl3btm1jzJgxjBo1CrDPZv7RRx8xYcKETu9j5cqV/OUvfzlnnd3xwgsv8Mc//hGAe+65h4ceegiAX/ziF6xYsYLQ0FAiIyNJTU3lhz/8YZfOLS8vj9WrV/PEE0/wwgsvdKqe7nzmygXDzxjzHPANoAnIAu4UkQqnFgX24ckyP7I3cJlwHQyJc3ZFSvW4mJgYli1bRlRUFLfffjsjRowgKirKsf7nP/85Tz/9NMXFxezcuZNVq1YBEB8fz/bt28/aX2dmfc/PzycyMtKxLiIigq1bt7ZZnzGGefPmYYzhO9/5Dvfeey9NTU1kZ2d3qs7zlZ6ezp/+9Ce2bt2KiDBlyhRmzpyJxWLh/fffZ/fu3TQ3NzNx4kRSU//zH4zOnttDDz3Eb37zmzY/q7bOGdr/zFXnuFz4AWuAn4iIxRjza+AnwGNOq0YEjn0NOf8G30GQvMQ+OLVSvehcV2i9ae/evSxcuJDS0lKCgoJarZsxYwYiwgsvvMD69esdtxnd3d3x8vKiurqagIAAx/YbNmzo0do2btxIeHg4xcXFXHnllYwfP54xY8Z0us5TpkyZQmNjIzU1NZw8eZLk5GQAfv3rXzN//vw2j3v99dfj5+cH2OfE27BhAzabjYULF+Lj44OPjw/f+MY3unxOH3/8MWFhYaSmprJ+/fpOnfOMGTPa/cxV57hc+InIF6e93ALc6KxaaK6H/R9D2REIi4WYa8DDy2nlKHUh7Nu3j/j4eBoaGmhoaGi1bu/evRQWFhISEnLWL9zGxkZ8fFqPZtSZK7/w8HCOHz/uWJeXl0d4eHibtZ1aHhYWxvXXX8+2bdtISEjoUp2A4+pr/fr1LF++nOXLl7d5vO7qzLlt2rSJVatW8cknn9DQ0EBVVRVLly5lxYoVjn1A63OeMWMG0PZnrjqprendXeUL+CewtIP19wJpQNqIESPOc+L7dlQWiGx+TWT9r0WOp4nYbD27f6XOkJmZ6ewSpKqqSsaPH+94HRERIfX19SIiUlBQIAkJCZKZmSlz586VTz/91LFdaWmpxMTEnNcxm5ubJTo6WrKzs6WxsVESExMlIyPjrO1qamqkqqrK8f2ll17qqKGzdZ5p3bp1smzZsnbX+/n5iYhIenq6JCQkSG1trdTU1EhcXJzs2LFDtm3bJikpKVJfXy/V1dUyduxYee6557p8bqfXc+2113bqnLvzmV+s2voZAtKkjfxwSmtPY8xaY0xGG18LT9vmCcACrGxvPyLyhohMEpFJoaGhPVOciH2W9Z1vAwIpSyEiVYcpU/1CRkYG8fHxjtfz5s1j48aN1NXVsXjxYp5//nliY2N56qmn+PnPf+7Ybt26dVx77bXndUwPDw9+97vfMX/+fGJjY/nmN79JXJz9mfo111xDQUEBACdOnGD69OkkJSVxySWXcO2113LVVVd1qc7zNXHiRO644w4uueQSpkyZwj333ENKSgqTJ0/muuuuIzExkauvvpqEhAQCAwO7fG7t6eicu/OZK1zzyg+4A/gaGNDZ96SmpnbnPwx2zY0i+z4S+fJXIrv+KtJY2/19KtVJrnDld6b09HRZunTpObe7/vrr5eDBgxegorZ1ts7eUF1dLSIitbW1kpqaKunp6RfkuM7+zF1RV678XO6ZnzHmKuDHwEwRqbtgB26ogj3vQl0ZRF8OI6fp1Z7q9yZOnMjs2bOxWq3t9qFrampi0aJFjBs37gJX9x+dqbO33HvvvWRmZtLQ0MCyZcuYOHFirx/TFT7zvs7Yg9F1GGOOAN5AWcuiLSJy37neN2nSJDm902uX2az2ocrCU2FQ9PnvR6nztH///nY7OSulzq2tnyFjTLqITDpzW5e78hORMU45sJs7JDivYalSSqkLR4c3U0op1e9o+CnlQlztMYRSfUVXf3Y0/JRyET4+PpSVlWkAKtVFIkJZWVmXOvy73DM/pfqriIgI8vLyKCkpcXYpSvU5Pj4+REREdHp7DT+lXISnpyfR0drSWKkLQW97KqWU6nc0/JRSSvU7Gn5KKaX6HZcb4eV8GWNKgKPOrqOTBgOlzi6ih+i5uJ6L5TxAz8UV9bXzGCkiZ818cNGEX19ijElra7idvkjPxfVcLOcBei6u6GI5D73tqZRSqt/R8FNKKdXvaPg5xxvOLqAH6bm4novlPEDPxRVdFOehz/yUUkr1O3rlp5RSqt/R8FNKKdXvaPhdQMaYSGPMOmNMpjFmnzHmQWfX1B3GGHdjzE5jzMfOrqU7jDFBxpj3jDEHjDH7jTGXOrum82WMebjl31aGMeYdY0znh7l3MmPMH40xxcaYjNOWDTLGrDHGHG75M9iZNXZGO+fxXMu/rz3GmH8YY4KcWGKntXUup6171BgjxpjBzqituzT8LiwL8KiITACmAvcbYyY4uabueBDY7+wiesDLwGciMh5Ioo+ekzEmHPgBMElE4gF34BbnVtUly4Grzlj2X8C/RGQs8K+W165uOWefxxogXkQSgUPATy50UedpOWefC8aYSGAecOxCF9RTNPwuIBEpFJEdLd9XY/8lG+7cqs6PMSYCuBZ4y9m1dIcxJhCYAfwBQESaRKTCqUV1jwfga4zxAAYABU6up9NE5N/AyTMWLwT+3PL9n4FFF7Km89HWeYjIFyJiaXm5Bej83DtO1M7fCcCLwI+BPttiUsPPSYwxUUAKsNXJpZyvl7D/47c5uY7uigZKgD+13MJ9yxjj5+yizoeI5AO/xf6/8UKgUkS+cG5V3TZERApbvi8ChjizmB5yF/Cps4s4X8aYhUC+iOx2di3doeHnBMYYf+B94CERqXJ2PV1ljFkAFItIurNr6QEewETg/4lIClBL37i1dpaW52ELsQf6cMDPGLPUuVX1HLH3y+qzVxoAxpgnsD/+WOnsWs6HMWYA8DjwtLNr6S4NvwvMGOOJPfhWisgHzq7nPE0DrjPG5AJ/Ba4wxqxwbknnLQ/IE5FTV+DvYQ/DvmgukCMiJSLSDHwAXObkmrrrhDFmGEDLn8VOrue8GWPuABYAS6TvdrAejf0/V7tbfv4jgB3GmKFOreo8aPhdQMYYg/3Z0n4RecHZ9ZwvEfmJiESISBT2BhVfikifvMIQkSLguDEmpmXRHCDTiSV1xzFgqjFmQMu/tTn00cY7p1kFLGv5fhnwkRNrOW/GmKuwPya4TkTqnF3P+RKRvSISJiJRLT//ecDElp+jPkXD78KaBtyG/UppV8vXNc4uSvEAsNIYswdIBn7l3HLOT8vV63vADmAv9p/vPjMUlTHmHeBrIMYYk2eMuRt4FrjSGHMY+5Xts86ssTPaOY/fAQHAmpaf+9edWmQntXMuFwUd3kwppVS/o1d+Siml+h0NP6WUUv2Ohp9SSql+R8NPKaVUv6Php5RSqt/R8FNKKdXvaPgppZTqdzT8lHIhxpj1xpgnnV2HUhc7DT+l+gFjTIExxnpqnMyWZR7GmDpjzJXOrE0pZ9DwU+oi1zLJ7TDsk6jedNqqOMAXSHNGXUo5k4afUi7KGBNijPk/Y0xRy9efjTGDTls/1BjzT2NMpTHmkDHmbmOMtMwVebrJQBnwPHDzacsnAVkiUt7rJ6OUi9HwU8p1rQSCgdiWr8HA22esbwIigenYB01vy2TsV3cfAJOMMZFnLFeq39HwU8oFGWOGA/OBR0SkvOXq7BHgGmPMMGNMBHAF8CMRqRKRYuAX7exuMpAmIieBdcA3W5ZPAra3c/zBxpi+Pgu8Uu3S8FPKNZ26Oss5bVnWaevCW74/dtr6o+3s6/SQ+xtwszHGC0ig/Su/JGBnVwpWqi/R8FPKNR1v+TPqtGWjTluX3/L9iNPWn/49AMaYMdhvnZ4KuX9gD7ZFgAeQftq21xljvjbGrMM+SbGGn7poafgp5YJEpAD4AnjeGBNkjAnG3mDlUxEpFJE8YD3wrDEmwBgTCrTVP3AyUCgi+S37LQf+hX3C3oMiUgPQ0kjmYeyzv88FrkTDT13ENPyUcl1LgWrgIHAAqABuP239t4ABQB6wCfh7y/LG07Zpq1HLu8DoM5YvBf4gInUiYgXqgcM9chZKuSCdyV2pi4QxZj7wEeArXfzBNsb8D3BYRP5ojFkCfE9EpvVGnUq5Ag0/pfooY0wyYAP2AtHYr+gyRWTZeewrFngf+7NEAQ6JyPd7rlqlXIuGn1J9lDFmNvAm9tFbKoFPgUdFpMKZdSnVF2j4KaWU6ne0wYtSSql+R8NPKaVUv6Php5RSqt/R8FNKKdXvaPgppZTqdzT8lFJK9Tsafkoppfqd/w9mXMqmkK1hegAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7,4))\n",
    "plt.plot(logNd,logNre,\".\",label=\"Table 10.1 in Pruppacher and Klett\")\n",
    "\n",
    "xarr=np.linspace(1,logNdc,100)\n",
    "plt.plot(xarr,xarr - np.log(24.),alpha=0.5,label=\"Stokes flow: $f(x)=x-\\log{24}$\")\n",
    "xarr=np.linspace(logNdc,logNde,100)\n",
    "plt.plot(xarr,-0.0088*xarr**2+0.85*xarr-2.49,alpha=0.5,label=\"$f(x)=-0.0088 x^2+0.85 x-2.49$\")\n",
    "plt.plot(xarr,-2.7905+0.9209*xarr-0.0135*xarr**2,label=\"petitRadtrans\",ls=\"dotted\",alpha=0.5)\n",
    "plt.plot(xarr,0.8*xarr-0.01*xarr**2,label=\"$y=0.8x-0.01x^2$ (AM01)\",alpha=0.5)\n",
    "\n",
    "xarr=np.linspace(logNde,15,100)\n",
    "plt.plot(xarr,0.5*(xarr-np.log(0.45)) ,alpha=0.5,label=\"$f(x)=0.5(x+\\\\log{0.45})$  \")\n",
    "plt.xlabel(\"$\\\\log{N_d}$\",fontsize=13)\n",
    "plt.ylabel(\"$\\\\log{N_{re}}$\",fontsize=13)\n",
    "plt.legend(loc=\"lower right\")\n",
    "plt.savefig(\"davies_reynolds.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that there is a typo (?) in Akerman and Marley (2001), tagged by \"AM01\". \n",
    "\n",
    "Using this relation, we can compute the Reynolds number, then we can also compute the terminal velocity using\n",
    "\n",
    "$v_f(r) = \\frac{2}{9 \\eta}  g r^2 (\\rho_c - \\rho) \\left( \\frac{C_d N_{re}}{24} \\right)^{-1}$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's how exojax compute the terminal velocity in exojax.atm.vterm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}