{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 作业11"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "假设在红移为1处去观测一个现在红移为2的星系，其红移为多少？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def get_z_obs(z_emit, z_observer):\n",
    "    return (1 + z_emit) / (1 + z_observer) - 1\n",
    "get_z_obs(2, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "假设宇宙的膨胀是一个线性过程，请根据今天的哈勃常数估算宇宙的年龄。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "估计的宇宙年龄为：13.9685 Gyr\n"
     ]
    }
   ],
   "source": [
    "from astropy import units as u\n",
    "H0 = 70 * u.km / (u.s * u.Mpc)  # 假设宇宙线性膨胀，那么H(z)就假设为H0\n",
    "# 计算哈勃时间（宇宙年龄），即哈勃常数的倒数\n",
    "T_Hubble = (1 / H0).to(u.Gyr)\n",
    "print(f\"估计的宇宙年龄为：{T_Hubble:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 作业12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一团1000个太阳质量的中性气体云（全部为H原子）, 温度为30K，当云的密度大于多少时，气体云将发生塌缩？塌缩（自由下落）时标为多少？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "分子云坍缩临界条件：$\\frac{GMm_H}{R}=\\frac{3KT}{2}$\n",
    "\n",
    "密度$n$与半径$R$的关系：$n\\dfrac{4\\pi R^3}{3}m_H = 1000M_{\\odot}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R: 3.576e+19 cm\n",
      "大于rho_crit: 6.202 1 / cm3时坍缩\n",
      "坍缩时标tau: 0.026 Gyr\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from astropy import units as u\n",
    "from astropy.constants import k_B, G, M_sun\n",
    "from astropy.units import Quantity\n",
    "\n",
    "M_J = 1000 * M_sun\n",
    "T = 30 * u.K\n",
    "m_H = 1.674e-27 * u.kg  # 氢原子质量\n",
    "\n",
    "R = (G*M_J*m_H)/(3/2*k_B*T)\n",
    "print(f\"R: {R.to(u.cm):.3e}\")\n",
    "n = M_J/((4/3)*np.pi*R**3*m_H)\n",
    "print(f\"大于rho_crit: {n.to(1/u.cm**3):.3f}时坍缩\")\n",
    "\n",
    "tau = np.sqrt(2*R**3/(G*M_J))\n",
    "print(f\"坍缩时标tau: {tau.to(u.Gyr):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这团气体中能形成的最大质量的恒星的质量是多少？\n",
    "\n",
    "• 假设分子云中的恒星形成遵循Salpter初始质量函数（最小质量的恒星为0.08太阳质量）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$M_\\text{total} = \\frac{k_2}{0.35}(M_\\text{min}^{-0.35} - M_\\text{max}^{-0.35})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这个应该不对"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10000/10000 [00:03<00:00, 3075.55it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximal mass: 192.16 M_sun\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAGwCAYAAAB4h2vpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABBHElEQVR4nO3dfVhVdb7//9cWELSQUoQNhYCOlUY1Ch6ECbEbUTTNtKKaSPLmGqYbQ3JStL46Tke0IYfL8W4qtTyVcTpY40mmwFMwNmKjiOY4jGNXKKbsYaAESwPE9fvDH3va7Q0CwdojPB/Xta72/qz3+qzPWi7Yr9ZarG0xDMMQAAAAulwvdw8AAACgpyB4AQAAmITgBQAAYBKCFwAAgEkIXgAAACYheAEAAJiE4AUAAGAST3cPoDu7cOGCTp06JV9fX1ksFncPBwAAtIFhGDpz5oyCg4PVq1fnnqMieHWhU6dOKSQkxN3DAAAAHXDixAlde+21ndonwasL+fr6Srr4D9evXz83jwYAALRFXV2dQkJC7J/jnYng1YWaLy/269eP4AUAwGWmK24T4uZ6AAAAkxC8AAAATELwAgAAMAnBCwAAwCQELwAAAJMQvAAAAExC8AIAADAJwQsAAMAkBC8AAACTELwAAABMQvACAAAwCcELAADAJAQvAAAAkxC8AAAATELwAgAAMImnuweAjgtbuMPdQ2i3YysmuXsIAAC4DWe8AAAATOL24LVu3TqFh4fLx8dHkZGR2rVrV6v1RUVFioyMlI+PjwYPHqwNGzY4zH/55ZcVFxenq6++WldffbXuvPNO/fnPf273eg3D0NKlSxUcHKw+ffpo7NixOnz48A/fYAAA0GO5NXjl5OQoLS1NixcvVmlpqeLi4pSYmKiKigqX9eXl5Zo4caLi4uJUWlqqRYsWae7cucrNzbXXFBYW6sEHH9RHH32k4uJiDRo0SAkJCTp58mS71vvCCy9o1apVWrNmjfbu3Sur1apx48bpzJkzXbdDAABAt2YxDMNw18qjo6M1cuRIrV+/3t42bNgwTZ06VZmZmU71CxYs0Pbt21VWVmZvS01N1cGDB1VcXOxyHU1NTbr66qu1Zs0aPfLII21ar2EYCg4OVlpamhYsWCBJqq+vV2BgoFauXKmf/exnLtdVX1+v+vp6+/u6ujqFhISotrZW/fr1a8eeaRvu8QIAoPPV1dXJz8+vSz6/3XbGq6GhQSUlJUpISHBoT0hI0O7du10uU1xc7FQ/fvx47du3T42NjS6XOXv2rBobG9W/f/82r7e8vFw2m82hxtvbW/Hx8S2OTZIyMzPl5+dnn0JCQlqsBQAAPY/bgld1dbWampoUGBjo0B4YGCibzeZyGZvN5rL+/Pnzqq6udrnMwoULdc011+jOO+9s83qb/9uesUlSRkaGamtr7dOJEydarAUAAD2P2x8nYbFYHN4bhuHUdql6V+3Sxfu0tm7dqsLCQvn4+LR7ve0dm7e3t7y9vVucDwAAeja3nfHy9/eXh4eH0xmkqqoqpzNNzaxWq8t6T09PDRgwwKE9KytLy5cvV35+vm6++eZ2rddqtUpSu8YGAABwKW4LXr1791ZkZKQKCgoc2gsKChQbG+tymZiYGKf6/Px8RUVFycvLy97261//Wr/61a/0/vvvKyoqqt3rDQ8Pl9VqdahpaGhQUVFRi2MDAAC4FLdeakxPT1dycrKioqIUExOjl156SRUVFUpNTZV08Z6pkydPasuWLZIu/gXjmjVrlJ6erjlz5qi4uFgbN27U1q1b7X2+8MILeu655/Tmm28qLCzMftbqyiuv1JVXXtmm9VosFqWlpWn58uUaOnSohg4dquXLl6tv37566KGHzNxFAACgG3Fr8EpKSlJNTY2WLVumyspKRUREKC8vT6GhoZKkyspKh2drhYeHKy8vT/PmzdPatWsVHBys1atXa/r06faadevWqaGhQffee6/DupYsWaKlS5e2ab2S9Mwzz+jcuXN67LHH9NVXXyk6Olr5+fny9fXtwj0CAAC6M7c+x6u768rngEg8xwsAgK7QLZ/jBQAA0NMQvAAAAExC8AIAADAJwQsAAMAkBC8AAACTELwAAABMQvACAAAwCcELAADAJAQvAAAAkxC8AAAATELwAgAAMAnBCwAAwCQELwAAAJMQvAAAAExC8AIAADAJwQsAAMAkBC8AAACTELwAAABMQvACAAAwCcELAADAJAQvAAAAkxC8AAAATELwAgAAMAnBCwAAwCQELwAAAJMQvAAAAExC8AIAADAJwQsAAMAkBC8AAACTELwAAABMQvACAAAwCcELAADAJG4PXuvWrVN4eLh8fHwUGRmpXbt2tVpfVFSkyMhI+fj4aPDgwdqwYYPD/MOHD2v69OkKCwuTxWJRdna2Ux/N874/Pf744/aalJQUp/mjR4/ulG0GAAA9k1uDV05OjtLS0rR48WKVlpYqLi5OiYmJqqiocFlfXl6uiRMnKi4uTqWlpVq0aJHmzp2r3Nxce83Zs2c1ePBgrVixQlar1WU/e/fuVWVlpX0qKCiQJN13330OdRMmTHCoy8vL66QtBwAAPZGnO1e+atUqzZo1S7Nnz5YkZWdn64MPPtD69euVmZnpVL9hwwYNGjTIfhZr2LBh2rdvn7KysjR9+nRJ0qhRozRq1ChJ0sKFC12ud+DAgQ7vV6xYoSFDhig+Pt6h3dvbu8XwBgAA0F5uO+PV0NCgkpISJSQkOLQnJCRo9+7dLpcpLi52qh8/frz27dunxsbGDo/j9ddf18yZM2WxWBzmFRYWKiAgQNddd53mzJmjqqqqVvuqr69XXV2dwwQAANDMbcGrurpaTU1NCgwMdGgPDAyUzWZzuYzNZnNZf/78eVVXV3doHO+++65Onz6tlJQUh/bExES98cYb+vDDD/Xiiy9q7969uv3221VfX99iX5mZmfLz87NPISEhHRoTAADontx6qVGS01kmwzCc2i5V76q9rTZu3KjExEQFBwc7tCclJdlfR0REKCoqSqGhodqxY4emTZvmsq+MjAylp6fb39fV1RG+AACAnduCl7+/vzw8PJzOblVVVTmd1WpmtVpd1nt6emrAgAHtHsPx48e1c+dObdu27ZK1QUFBCg0N1dGjR1us8fb2lre3d7vHAQAAega3XWrs3bu3IiMj7X9R2KygoECxsbEul4mJiXGqz8/PV1RUlLy8vNo9hs2bNysgIECTJk26ZG1NTY1OnDihoKCgdq8HAABAcvPjJNLT0/XKK69o06ZNKisr07x581RRUaHU1FRJFy/dPfLII/b61NRUHT9+XOnp6SorK9OmTZu0ceNGzZ8/317T0NCgAwcO6MCBA2poaNDJkyd14MABffbZZw7rvnDhgjZv3qwZM2bI09PxxN/XX3+t+fPnq7i4WMeOHVNhYaEmT54sf39/3XPPPV24RwAAQHfm1nu8kpKSVFNTo2XLlqmyslIRERHKy8tTaGioJKmystLhmV7h4eHKy8vTvHnztHbtWgUHB2v16tX2R0lI0qlTpzRixAj7+6ysLGVlZSk+Pl6FhYX29p07d6qiokIzZ850GpeHh4cOHTqkLVu26PTp0woKCtJtt92mnJwc+fr6dsGeAAAAPYHFaL47HZ2urq5Ofn5+qq2tVb9+/Tq9/7CFOzq9z652bMWlL+sCAOBOXfn57favDAIAAOgpCF4AAAAmIXgBAACYhOAFAABgEoIXAACASQheAAAAJiF4AQAAmITgBQAAYBKCFwAAgEkIXgAAACYheAEAAJiE4AUAAGASghcAAIBJCF4AAAAmIXgBAACYhOAFAABgEoIXAACASQheAAAAJiF4AQAAmITgBQAAYBKCFwAAgEkIXgAAACYheAEAAJiE4AUAAGASghcAAIBJCF4AAAAmIXgBAACYhOAFAABgEoIXAACASQheAAAAJiF4AQAAmITgBQAAYBK3B69169YpPDxcPj4+ioyM1K5du1qtLyoqUmRkpHx8fDR48GBt2LDBYf7hw4c1ffp0hYWFyWKxKDs726mPpUuXymKxOExWq9WhxjAMLV26VMHBwerTp4/Gjh2rw4cP/+DtBQAAPZdbg1dOTo7S0tK0ePFilZaWKi4uTomJiaqoqHBZX15erokTJyouLk6lpaVatGiR5s6dq9zcXHvN2bNnNXjwYK1YscIpTH3XjTfeqMrKSvt06NAhh/kvvPCCVq1apTVr1mjv3r2yWq0aN26czpw50zkbDwAAehy3Bq9Vq1Zp1qxZmj17toYNG6bs7GyFhIRo/fr1Lus3bNigQYMGKTs7W8OGDdPs2bM1c+ZMZWVl2WtGjRqlX//613rggQfk7e3d4ro9PT1ltVrt08CBA+3zDMNQdna2Fi9erGnTpikiIkKvvfaazp49qzfffLPzdgAAAOhR3Ba8GhoaVFJSooSEBIf2hIQE7d692+UyxcXFTvXjx4/Xvn371NjY2K71Hz16VMHBwQoPD9cDDzygzz//3D6vvLxcNpvNYV3e3t6Kj49vcWySVF9fr7q6OocJAACgmduCV3V1tZqamhQYGOjQHhgYKJvN5nIZm83msv78+fOqrq5u87qjo6O1ZcsWffDBB3r55Zdls9kUGxurmpoa+3qa+27r2CQpMzNTfn5+9ikkJKTNYwIAAN2f22+ut1gsDu8Nw3Bqu1S9q/bWJCYmavr06brpppt05513aseOHZKk11577QeNLSMjQ7W1tfbpxIkTbR4TAADo/jzdtWJ/f395eHg4nUGqqqpyOtPUzGq1uqz39PTUgAEDOjyWK664QjfddJOOHj1qX4908cxXUFBQm8YmXbwc2dp9ZQAAoGdz2xmv3r17KzIyUgUFBQ7tBQUFio2NdblMTEyMU31+fr6ioqLk5eXV4bHU19errKzMHrLCw8NltVod1tXQ0KCioqIWxwYAAHApbjvjJUnp6elKTk5WVFSUYmJi9NJLL6miokKpqamSLl66O3nypLZs2SJJSk1N1Zo1a5Senq45c+aouLhYGzdu1NatW+19NjQ06K9//av99cmTJ3XgwAFdeeWV+tGPfiRJmj9/viZPnqxBgwapqqpKzz//vOrq6jRjxgxJFy8xpqWlafny5Ro6dKiGDh2q5cuXq2/fvnrooYfM3EUAAKAbcWvwSkpKUk1NjZYtW6bKykpFREQoLy9PoaGhkqTKykqHZ3qFh4crLy9P8+bN09q1axUcHKzVq1dr+vTp9ppTp05pxIgR9vdZWVnKyspSfHy8CgsLJUlffPGFHnzwQVVXV2vgwIEaPXq09uzZY1+vJD3zzDM6d+6cHnvsMX311VeKjo5Wfn6+fH19u3ivAACA7spiNN+djk5XV1cnPz8/1dbWql+/fp3ef9jCHZ3eZ1c7tmKSu4cAAECruvLz2+1/1QgAANBTELwAAABMQvACAAAwCcELAADAJAQvAAAAkxC8AAAATELwAgAAMAnBCwAAwCQELwAAAJMQvAAAAExC8AIAADAJwQsAAMAkBC8AAACTELwAAABMQvACAAAwCcELAADAJAQvAAAAkxC8AAAATOLp7gEAAAD3C1u4w91DaLdjKya5ewjtxhkvAAAAkxC8AAAATELwAgAAMAnBCwAAwCQELwAAAJMQvAAAAExC8AIAADAJwQsAAMAkBC8AAACTELwAAABMQvACAAAwCcELAADAJAQvAAAAkxC8AAAATOL24LVu3TqFh4fLx8dHkZGR2rVrV6v1RUVFioyMlI+PjwYPHqwNGzY4zD98+LCmT5+usLAwWSwWZWdnO/WRmZmpUaNGydfXVwEBAZo6daqOHDniUJOSkiKLxeIwjR49+gdvLwAA6LncGrxycnKUlpamxYsXq7S0VHFxcUpMTFRFRYXL+vLyck2cOFFxcXEqLS3VokWLNHfuXOXm5tprzp49q8GDB2vFihWyWq0u+ykqKtLjjz+uPXv2qKCgQOfPn1dCQoK++eYbh7oJEyaosrLSPuXl5XXexgMAgB7HsyMLvfrqq7r//vvVt2/fH7TyVatWadasWZo9e7YkKTs7Wx988IHWr1+vzMxMp/oNGzZo0KBB9rNYw4YN0759+5SVlaXp06dLkkaNGqVRo0ZJkhYuXOhyve+//77D+82bNysgIEAlJSUaM2aMvd3b27vF8OZKfX296uvr7e/r6uravCwAAOj+OnTGKyMjQ1arVbNmzdLu3bs7tOKGhgaVlJQoISHBoT0hIaHFPouLi53qx48fr3379qmxsbFD45Ck2tpaSVL//v0d2gsLCxUQEKDrrrtOc+bMUVVVVav9ZGZmys/Pzz6FhIR0eEwAAKD76VDw+uKLL/T666/rq6++0m233aYbbrhBK1eulM1ma3Mf1dXVampqUmBgoEN7YGBgi/3YbDaX9efPn1d1dXX7N0SSYRhKT0/XrbfeqoiICHt7YmKi3njjDX344Yd68cUXtXfvXt1+++0OZ7S+LyMjQ7W1tfbpxIkTHRoTAADonjp0qdHDw0NTpkzRlClTVFVVpddff12vvvqqnnvuOU2YMEGzZs3S5MmT1avXpXOdxWJxeG8YhlPbpepdtbfVE088oU8//VQff/yxQ3tSUpL9dUREhKKiohQaGqodO3Zo2rRpLvvy9vaWt7d3h8YBAAC6vx98c31AQIB+8pOfKCYmRr169dKhQ4eUkpKiIUOGqLCwsMXl/P395eHh4XR2q6qqyumsVjOr1eqy3tPTUwMGDGj32J988klt375dH330ka699tpWa4OCghQaGqqjR4+2ez0AAADSDwhe//jHP5SVlaUbb7xRY8eOVV1dnd577z2Vl5fr1KlTmjZtmmbMmNHi8r1791ZkZKQKCgoc2gsKChQbG+tymZiYGKf6/Px8RUVFycvLq81jNwxDTzzxhLZt26YPP/xQ4eHhl1ympqZGJ06cUFBQUJvXAwAA8F0dCl6TJ09WSEiIXn31Vc2ZM0cnT57U1q1bdeedd0qS+vTpo6effvqS9zilp6frlVde0aZNm1RWVqZ58+apoqJCqampki7eM/XII4/Y61NTU3X8+HGlp6errKxMmzZt0saNGzV//nx7TUNDgw4cOKADBw6ooaFBJ0+e1IEDB/TZZ5/Zax5//HG9/vrrevPNN+Xr6yubzSabzaZz585Jkr7++mvNnz9fxcXFOnbsmAoLCzV58mT5+/vrnnvu6cguAwAA6Ng9XgEBASoqKlJMTEyLNUFBQSovL2+1n6SkJNXU1GjZsmWqrKxURESE8vLyFBoaKkmqrKx0eKZXeHi48vLyNG/ePK1du1bBwcFavXq1/VESknTq1CmNGDHC/j4rK0tZWVmKj4+3X/pcv369JGns2LEO49m8ebNSUlLk4eGhQ4cOacuWLTp9+rSCgoJ02223KScnR76+vm3aRwAAAN9nMZrvTm+HLVu2KCkpyelG8oaGBr311lsOZ6l6srq6Ovn5+am2tlb9+vXr9P7DFu7o9D672rEVk9w9BACAC3ym/EtXfn536FLjo48+an/21XedOXNGjz766A8eFAAAQHfUoeDV0iMfvvjiC/n5+f3gQQEAAHRH7brHa8SIEfYvjL7jjjvk6fmvxZuamlReXq4JEyZ0+iABAAC6g3YFr6lTp0qSDhw4oPHjx+vKK6+0z+vdu7fCwsIcbnQHAADAv7QreC1ZskSSFBYWpqSkJPn4+HTJoAAAALqjDj1OorUHowIAAMC1Ngev/v376+9//7v8/f119dVXt/rdiF9++WWnDA4AAKA7aXPw+s1vfmN/eOhvfvObDn8pNQAAQE/V5uD13cuLKSkpXTEWAACAbq3Nwauurq7NnXbFU9oBAAAud20OXlddddUlLy82P1i1qanpBw8MAACgu2lz8Proo4+6chwAAADdXpuDV3x8fFeOAwAAoNtrc/D69NNPFRERoV69eunTTz9ttfbmm2/+wQMDAADobtocvH784x/LZrMpICBAP/7xj2WxWGQYhlMd93gBAAC41ubgVV5eroEDB9pfAwAAoH3aHLxCQ0NdvgYAAEDbdOi7GiXpyJEj+u1vf6uysjJZLBbdcMMNevLJJ3X99dd35vgAAAC6jV4dWeh//ud/FBERoZKSEt1yyy26+eabtX//fkVEROjtt9/u7DECAAB0Cx064/XMM88oIyNDy5Ytc2hfsmSJFixYoPvuu69TBgcAANCddOiMl81m0yOPPOLU/vDDD8tms/3gQQEAAHRHHQpeY8eO1a5du5zaP/74Y8XFxf3gQQEAAHRHbb7UuH37dvvrKVOmaMGCBSopKdHo0aMlSXv27NHbb7+tX/7yl50/SgAAgG7AYrh6CqoLvXq17eQYD1D9l7q6Ovn5+am2tlb9+vXr9P7DFu7o9D672rEVk9w9BACAC3ym/EtXfn63+YzXhQsXOnXFAAAAPU2H7vECAABA+3X4AarffPONioqKVFFRoYaGBod5c+fO/cEDAwAA6G46FLxKS0s1ceJEnT17Vt9884369++v6upq9e3bVwEBAQQvAAAAFzp0qXHevHmaPHmyvvzyS/Xp00d79uzR8ePHFRkZqaysrM4eIwAAQLfQoeB14MABPf300/Lw8JCHh4fq6+sVEhKiF154QYsWLersMQIAAHQLHQpeXl5eslgskqTAwEBVVFRIkvz8/OyvAQAA4KhD93iNGDFC+/bt03XXXafbbrtN/+///T9VV1frv/7rv3TTTTd19hgBAAC6hQ6d8Vq+fLmCgoIkSb/61a80YMAA/fznP1dVVZVeeumldvW1bt06hYeHy8fHR5GRkS6/iui7ioqKFBkZKR8fHw0ePFgbNmxwmH/48GFNnz5dYWFhslgsys7O7tB6DcPQ0qVLFRwcrD59+mjs2LE6fPhwu7YNAADguzoUvKKionTbbbdJkgYOHKi8vDzV1dVp//79uuWWW9rcT05OjtLS0rR48WKVlpYqLi5OiYmJLV6uLC8v18SJExUXF6fS0lItWrRIc+fOVW5urr3m7NmzGjx4sFasWCGr1drh9b7wwgtatWqV1qxZo71798pqtWrcuHE6c+ZMm7cPAADgu9r8lUGuVFVV6ciRI7JYLLr++us1cODAdi0fHR2tkSNHav369fa2YcOGaerUqcrMzHSqX7BggbZv366ysjJ7W2pqqg4ePKji4mKn+rCwMKWlpSktLa1d6zUMQ8HBwUpLS9OCBQskSfX19QoMDNTKlSv1s5/9rE3bx1cGOeMrgwDg3xOfKf/SlZ/fHTrjVVdXp+TkZF1zzTWKj4/XmDFjFBwcrIcffli1tbVt6qOhoUElJSVKSEhwaE9ISNDu3btdLlNcXOxUP378eO3bt0+NjY2dtt7y8nLZbDaHGm9vb8XHx7c4NuliOKurq3OYAAAAmnUoeM2ePVuffPKJ3nvvPZ0+fVq1tbV67733tG/fPs2ZM6dNfVRXV6upqUmBgYEO7YGBgbLZbC6XsdlsLuvPnz+v6urqTltv83/bMzZJyszMlJ+fn30KCQlp05gAAEDP0KHgtWPHDm3atEnjx49Xv3795Ovrq/Hjx+vll1/Wjh3tO1XZ/FiKZoZhOLVdqt5Ve2est71jy8jIUG1trX06ceJEu8YEAAC6tw49TmLAgAHy8/Nzavfz89PVV1/dpj78/f3l4eHhdAapqqrK6UxTM6vV6rLe09NTAwYM6LT1Nt+Ub7PZ7H+9eamxSRcvR3p7e7dpHAAAoOfp0BmvZ599Vunp6aqsrLS32Ww2/eIXv9Bzzz3Xpj569+6tyMhIFRQUOLQXFBQoNjbW5TIxMTFO9fn5+YqKipKXl1enrTc8PFxWq9WhpqGhQUVFRS2ODQAA4FLafMZrxIgRDpfZjh49qtDQUA0aNEiSVFFRIW9vb/3zn/9s81/9paenKzk5WVFRUYqJidFLL72kiooKpaamSrp46e7kyZPasmWLpIt/wbhmzRqlp6drzpw5Ki4u1saNG7V161Z7nw0NDfrrX/9qf33y5EkdOHBAV155pX70ox+1ab0Wi0VpaWlavny5hg4dqqFDh2r58uXq27evHnroobbuMgAAAAdtDl5Tp07t9JUnJSWppqZGy5YtU2VlpSIiIpSXl6fQ0FBJUmVlpcOztcLDw5WXl6d58+Zp7dq1Cg4O1urVqzV9+nR7zalTpzRixAj7+6ysLGVlZSk+Pl6FhYVtWq8kPfPMMzp37pwee+wxffXVV4qOjlZ+fr58fX07fT8AAICe4Qc9xwut4zlezniOFwD8e+Iz5V+68vO7QzfXNyspKVFZWZksFouGDx/ucKYJAAAAjjoUvKqqqvTAAw+osLBQV111lQzDUG1trW677Ta99dZb7X6CPQAAQE/Qob9qfPLJJ1VXV6fDhw/ryy+/1FdffaW//OUvqqur09y5czt7jAAAAN1Ch854vf/++9q5c6eGDRtmbxs+fLjWrl3r9FU8AAAAuKhDZ7wuXLjg8rlZXl5eunDhwg8eFAAAQHfUoeB1++2366mnntKpU6fsbSdPntS8efN0xx13dNrgAAAAupMOBa81a9bozJkzCgsL05AhQ/SjH/1I4eHhOnPmjH7729929hgBAAC6hQ7d4xUSEqL9+/eroKBAf/vb32QYhoYPH64777yzs8cHAADQbbQ7eJ0/f14+Pj46cOCAxo0bp3HjxnXFuAAAALqddl9q9PT0VGhoqJqamrpiPAAAAN1Wh+7xevbZZ5WRkaEvv/yys8cDAADQbXXoHq/Vq1frs88+U3BwsEJDQ3XFFVc4zN+/f3+nDA4AAKA76VDwmjp1qiwWi/h+bQAAgLZrV/A6e/asfvGLX+jdd99VY2Oj7rjjDv32t7+Vv79/V40PAACg22jXPV5LlizRq6++qkmTJunBBx/Uzp079fOf/7yrxgYAANCttOuM17Zt27Rx40Y98MADkqSf/vSn+slPfqKmpiZ5eHh0yQABAAC6i3ad8Tpx4oTi4uLs7//jP/5Dnp6eDl8dBAAAANfaFbyamprUu3dvhzZPT0+dP3++UwcFAADQHbXrUqNhGEpJSZG3t7e97dtvv1VqaqrDIyW2bdvWeSMEAADoJtoVvGbMmOHU9vDDD3faYAAAALqzdgWvzZs3d9U4AAAAur0OfWUQAAAA2o/gBQAAYBKCFwAAgEkIXgAAACYheAEAAJiE4AUAAGASghcAAIBJCF4AAAAmIXgBAACYhOAFAABgEoIXAACASdwevNatW6fw8HD5+PgoMjJSu3btarW+qKhIkZGR8vHx0eDBg7VhwwanmtzcXA0fPlze3t4aPny43nnnHYf5YWFhslgsTtPjjz9ur0lJSXGaP3r06M7ZaAAA0CO5NXjl5OQoLS1NixcvVmlpqeLi4pSYmKiKigqX9eXl5Zo4caLi4uJUWlqqRYsWae7cucrNzbXXFBcXKykpScnJyTp48KCSk5N1//3365NPPrHX7N27V5WVlfapoKBAknTfffc5rG/ChAkOdXl5eV2wFwAAQE9hMQzDcNfKo6OjNXLkSK1fv97eNmzYME2dOlWZmZlO9QsWLND27dtVVlZmb0tNTdXBgwdVXFwsSUpKSlJdXZ3+8Ic/2GsmTJigq6++Wlu3bnU5jrS0NL333ns6evSoLBaLpItnvE6fPq133323w9tXV1cnPz8/1dbWql+/fh3upyVhC3d0ep9d7diKSe4eAgDABT5T/qUrP7/ddsaroaFBJSUlSkhIcGhPSEjQ7t27XS5TXFzsVD9+/Hjt27dPjY2Nrda01GdDQ4Nef/11zZw50x66mhUWFiogIEDXXXed5syZo6qqqla3qb6+XnV1dQ4TAABAM7cFr+rqajU1NSkwMNChPTAwUDabzeUyNpvNZf358+dVXV3dak1Lfb777rs6ffq0UlJSHNoTExP1xhtv6MMPP9SLL76ovXv36vbbb1d9fX2L25SZmSk/Pz/7FBIS0mItAADoeTzdPYDvn2UyDMOp7VL1329vT58bN25UYmKigoODHdqTkpLsryMiIhQVFaXQ0FDt2LFD06ZNc9lXRkaG0tPT7e/r6uoIXwAAwM5twcvf318eHh5OZ6Kqqqqczlg1s1qtLus9PT01YMCAVmtc9Xn8+HHt3LlT27Ztu+R4g4KCFBoaqqNHj7ZY4+3tLW9v70v2BQAAeia3XWrs3bu3IiMj7X9R2KygoECxsbEul4mJiXGqz8/PV1RUlLy8vFqtcdXn5s2bFRAQoEmTLn1zXk1NjU6cOKGgoKBL1gIAALji1sdJpKen65VXXtGmTZtUVlamefPmqaKiQqmpqZIuXrp75JFH7PWpqak6fvy40tPTVVZWpk2bNmnjxo2aP3++veapp55Sfn6+Vq5cqb/97W9auXKldu7cqbS0NId1X7hwQZs3b9aMGTPk6el44u/rr7/W/PnzVVxcrGPHjqmwsFCTJ0+Wv7+/7rnnnq7bIQAAoFtz6z1eSUlJqqmp0bJly1RZWamIiAjl5eUpNDRUklRZWenwTK/w8HDl5eVp3rx5Wrt2rYKDg7V69WpNnz7dXhMbG6u33npLzz77rJ577jkNGTJEOTk5io6Odlj3zp07VVFRoZkzZzqNy8PDQ4cOHdKWLVt0+vRpBQUF6bbbblNOTo58fX27aG8AAIDuzq3P8erueI6XM57jBQD/nvhM+Zdu+RwvAACAnobgBQAAYBKCFwAAgEkIXgAAACYheAEAAJjE7V8ZhJ6Fv5pBSzg2zHE57ufL0eV4bMAcnPECAAAwCcELAADAJAQvAAAAkxC8AAAATELwAgAAMAnBCwAAwCQELwAAAJMQvAAAAExC8AIAADAJwQsAAMAkBC8AAACTELwAAABMQvACAAAwCcELAADAJAQvAAAAkxC8AAAATELwAgAAMAnBCwAAwCQELwAAAJMQvAAAAExC8AIAADAJwQsAAMAkBC8AAACTELwAAABMQvACAAAwCcELAADAJAQvAAAAk7g9eK1bt07h4eHy8fFRZGSkdu3a1Wp9UVGRIiMj5ePjo8GDB2vDhg1ONbm5uRo+fLi8vb01fPhwvfPOOw7zly5dKovF4jBZrVaHGsMwtHTpUgUHB6tPnz4aO3asDh8+/MM3GAAA9FhuDV45OTlKS0vT4sWLVVpaqri4OCUmJqqiosJlfXl5uSZOnKi4uDiVlpZq0aJFmjt3rnJzc+01xcXFSkpKUnJysg4ePKjk5GTdf//9+uSTTxz6uvHGG1VZWWmfDh065DD/hRde0KpVq7RmzRrt3btXVqtV48aN05kzZzp/RwAAgB7BrcFr1apVmjVrlmbPnq1hw4YpOztbISEhWr9+vcv6DRs2aNCgQcrOztawYcM0e/ZszZw5U1lZWfaa7OxsjRs3ThkZGbrhhhuUkZGhO+64Q9nZ2Q59eXp6ymq12qeBAwfa5xmGoezsbC1evFjTpk1TRESEXnvtNZ09e1Zvvvlmi9tTX1+vuro6hwkAAKCZ24JXQ0ODSkpKlJCQ4NCekJCg3bt3u1ymuLjYqX78+PHat2+fGhsbW635fp9Hjx5VcHCwwsPD9cADD+jzzz+3zysvL5fNZnPox9vbW/Hx8S2OTZIyMzPl5+dnn0JCQlrZAwAAoKdxW/Cqrq5WU1OTAgMDHdoDAwNls9lcLmOz2VzWnz9/XtXV1a3WfLfP6OhobdmyRR988IFefvll2Ww2xcbGqqamxt5H83JtHZskZWRkqLa21j6dOHGitV0AAAB6GE93D8BisTi8NwzDqe1S9d9vv1SfiYmJ9tc33XSTYmJiNGTIEL322mtKT0/v8Ni8vb3l7e3d4nwAANCzue2Ml7+/vzw8PJzOIFVVVTmdaWpmtVpd1nt6emrAgAGt1rTUpyRdccUVuummm3T06FF7H5La3Q8AAEBr3Ba8evfurcjISBUUFDi0FxQUKDY21uUyMTExTvX5+fmKioqSl5dXqzUt9SldvCm+rKxMQUFBkqTw8HBZrVaHfhoaGlRUVNRqPwAAAK1x66XG9PR0JScnKyoqSjExMXrppZdUUVGh1NRUSRfvmTp58qS2bNkiSUpNTdWaNWuUnp6uOXPmqLi4WBs3btTWrVvtfT711FMaM2aMVq5cqbvvvlu///3vtXPnTn388cf2mvnz52vy5MkaNGiQqqqq9Pzzz6uurk4zZsyQdPESY1pampYvX66hQ4dq6NChWr58ufr27auHHnrIxD0EAAC6E7cGr6SkJNXU1GjZsmWqrKxURESE8vLyFBoaKkmqrKx0eKZXeHi48vLyNG/ePK1du1bBwcFavXq1pk+fbq+JjY3VW2+9pWeffVbPPfechgwZopycHEVHR9trvvjiCz344IOqrq7WwIEDNXr0aO3Zs8e+Xkl65plndO7cOT322GP66quvFB0drfz8fPn6+pqwZwAAQHdkMZrvTkenq6urk5+fn2pra9WvX79O7z9s4Y5O7xPOjq2Y5O4h9AiX4/F8OR4bl+N+vhxxbJijq/ZzV35+u/0rgwAAAHoKghcAAIBJCF4AAAAmIXgBAACYhOAFAABgEoIXAACASQheAAAAJiF4AQAAmITgBQAAYBKCFwAAgEkIXgAAACYheAEAAJiE4AUAAGASghcAAIBJCF4AAAAmIXgBAACYhOAFAABgEoIXAACASQheAAAAJiF4AQAAmITgBQAAYBKCFwAAgEkIXgAAACYheAEAAJiE4AUAAGASghcAAIBJCF4AAAAmIXgBAACYhOAFAABgEoIXAACASQheAAAAJiF4AQAAmMTtwWvdunUKDw+Xj4+PIiMjtWvXrlbri4qKFBkZKR8fHw0ePFgbNmxwqsnNzdXw4cPl7e2t4cOH65133nGYn5mZqVGjRsnX11cBAQGaOnWqjhw54lCTkpIii8XiMI0ePfqHbzAAAOix3Bq8cnJylJaWpsWLF6u0tFRxcXFKTExURUWFy/ry8nJNnDhRcXFxKi0t1aJFizR37lzl5ubaa4qLi5WUlKTk5GQdPHhQycnJuv/++/XJJ5/Ya4qKivT4449rz549Kigo0Pnz55WQkKBvvvnGYX0TJkxQZWWlfcrLy+uaHQEAAHoET3eufNWqVZo1a5Zmz54tScrOztYHH3yg9evXKzMz06l+w4YNGjRokLKzsyVJw4YN0759+5SVlaXp06fb+xg3bpwyMjIkSRkZGSoqKlJ2dra2bt0qSXr//fcd+t28ebMCAgJUUlKiMWPG2Nu9vb1ltVrbvD319fWqr6+3v6+rq2vzsgAAoPtz2xmvhoYGlZSUKCEhwaE9ISFBu3fvdrlMcXGxU/348eO1b98+NTY2tlrTUp+SVFtbK0nq37+/Q3thYaECAgJ03XXXac6cOaqqqmp1mzIzM+Xn52efQkJCWq0HAAA9i9uCV3V1tZqamhQYGOjQHhgYKJvN5nIZm83msv78+fOqrq5utaalPg3DUHp6um699VZFRETY2xMTE/XGG2/oww8/1Isvvqi9e/fq9ttvdzij9X0ZGRmqra21TydOnGh5BwAAgB7HrZcaJclisTi8NwzDqe1S9d9vb0+fTzzxhD799FN9/PHHDu1JSUn21xEREYqKilJoaKh27NihadOmuezL29tb3t7eLY4dAAD0bG4LXv7+/vLw8HA6E1VVVeV0xqqZ1Wp1We/p6akBAwa0WuOqzyeffFLbt2/XH//4R1177bWtjjcoKEihoaE6evToJbcNAADAFbddauzdu7ciIyNVUFDg0F5QUKDY2FiXy8TExDjV5+fnKyoqSl5eXq3WfLdPwzD0xBNPaNu2bfrwww8VHh5+yfHW1NToxIkTCgoKatP2AQAAfJ9bHyeRnp6uV155RZs2bVJZWZnmzZuniooKpaamSrp4z9Qjjzxir09NTdXx48eVnp6usrIybdq0SRs3btT8+fPtNU899ZTy8/O1cuVK/e1vf9PKlSu1c+dOpaWl2Wsef/xxvf7663rzzTfl6+srm80mm82mc+fOSZK+/vprzZ8/X8XFxTp27JgKCws1efJk+fv765577jFn5wAAgG7Hrfd4JSUlqaamRsuWLVNlZaUiIiKUl5en0NBQSVJlZaXDM73Cw8OVl5enefPmae3atQoODtbq1avtj5KQpNjYWL311lt69tln9dxzz2nIkCHKyclRdHS0vWb9+vWSpLFjxzqMZ/PmzUpJSZGHh4cOHTqkLVu26PTp0woKCtJtt92mnJwc+fr6duEeAQAA3Znbb65/7LHH9Nhjj7mc9+qrrzq1xcfHa//+/a32ee+99+ree+9tcX7zDfkt6dOnjz744INWawAAANrL7V8ZBAAA0FMQvAAAAExC8AIAADAJwQsAAMAkBC8AAACTELwAAABMQvACAAAwCcELAADAJAQvAAAAkxC8AAAATELwAgAAMAnBCwAAwCQELwAAAJMQvAAAAExC8AIAADAJwQsAAMAkBC8AAACTELwAAABMQvACAAAwCcELAADAJAQvAAAAkxC8AAAATELwAgAAMAnBCwAAwCQELwAAAJMQvAAAAExC8AIAADAJwQsAAMAkBC8AAACTELwAAABMQvACAAAwCcELAADAJG4PXuvWrVN4eLh8fHwUGRmpXbt2tVpfVFSkyMhI+fj4aPDgwdqwYYNTTW5uroYPHy5vb28NHz5c77zzTrvXaxiGli5dquDgYPXp00djx47V4cOHf9jGAgCAHs2twSsnJ0dpaWlavHixSktLFRcXp8TERFVUVLisLy8v18SJExUXF6fS0lItWrRIc+fOVW5urr2muLhYSUlJSk5O1sGDB5WcnKz7779fn3zySbvW+8ILL2jVqlVas2aN9u7dK6vVqnHjxunMmTNdt0MAAEC3ZjEMw3DXyqOjozVy5EitX7/e3jZs2DBNnTpVmZmZTvULFizQ9u3bVVZWZm9LTU3VwYMHVVxcLElKSkpSXV2d/vCHP9hrJkyYoKuvvlpbt25t03oNw1BwcLDS0tK0YMECSVJ9fb0CAwO1cuVK/exnP2vT9tXV1cnPz0+1tbXq169fO/ZM24Qt3NHpfcLZsRWT3D2EHuFyPJ4vx2PjctzPlyOODXN01X7uys9vz07trR0aGhpUUlKihQsXOrQnJCRo9+7dLpcpLi5WQkKCQ9v48eO1ceNGNTY2ysvLS8XFxZo3b55TTXZ2dpvXW15eLpvN5rAub29vxcfHa/fu3S0Gr/r6etXX19vf19bWSrr4D9gVLtSf7ZJ+4air/v3g6HI8ni/HY+Ny3M+XI44Nc3TVfm7utyvOTbkteFVXV6upqUmBgYEO7YGBgbLZbC6XsdlsLuvPnz+v6upqBQUFtVjT3Gdb1tv8X1c1x48fb3GbMjMz9ctf/tKpPSQkpMVl8O/PL9vdI8C/K44NtIRjwxxdvZ/PnDkjPz+/Tu3TbcGrmcVicXhvGIZT26Xqv9/elj47q+a7MjIylJ6ebn9/4cIFffnllxowYECry31fXV2dQkJCdOLEiS65RHm5YX84Y584Y584Y584Y584Yn84a94nf/3rXxUcHNzp/bstePn7+8vDw8Pp7FZVVZXTmaZmVqvVZb2np6cGDBjQak1zn21Zr9VqlXTxzFdQUFCbxiZdvBzp7e3t0HbVVVe1WH8p/fr14wfhO9gfztgnztgnztgnztgnjtgfzq655hr16tX5f4Potr9q7N27tyIjI1VQUODQXlBQoNjYWJfLxMTEONXn5+crKipKXl5erdY099mW9YaHh8tqtTrUNDQ0qKioqMWxAQAAXIpbLzWmp6crOTlZUVFRiomJ0UsvvaSKigqlpqZKunjp7uTJk9qyZYuki3/BuGbNGqWnp2vOnDkqLi7Wxo0b7X+tKElPPfWUxowZo5UrV+ruu+/W73//e+3cuVMff/xxm9drsViUlpam5cuXa+jQoRo6dKiWL1+uvn376qGHHjJxDwEAgG7FcLO1a9caoaGhRu/evY2RI0caRUVF9nkzZsww4uPjHeoLCwuNESNGGL179zbCwsKM9evXO/X59ttvG9dff73h5eVl3HDDDUZubm671msYhnHhwgVjyZIlhtVqNby9vY0xY8YYhw4d6pyNvoRvv/3WWLJkifHtt9+asr5/d+wPZ+wTZ+wTZ+wTZ+wTR+wPZ129T9z6HC8AAICexO1fGQQAANBTELwAAABMQvACAAAwCcELAADAJAQvN8nMzNSoUaPk6+urgIAATZ06VUeOHHGoSUlJkcVicZhGjx7tphF3vaVLlzptb/PDbKWL3xywdOlSBQcHq0+fPho7dqwOHz7sxhF3rbCwMKf9YbFY9Pjjj0vqGcfHH//4R02ePFnBwcGyWCx69913Hea35Zior6/Xk08+KX9/f11xxRWaMmWKvvjiCxO3onO1tk8aGxu1YMEC3XTTTbriiisUHBysRx55RKdOnXLoY+zYsU7HzgMPPGDylnSeSx0nbflZ6UnHiSSXv1ssFot+/etf22u603HSls9cs36fELzcpKioSI8//rj27NmjgoICnT9/XgkJCfrmm28c6iZMmKDKykr7lJeX56YRm+PGG2902N5Dhw7Z573wwgtatWqV1qxZo71798pqtWrcuHE6c+aMG0fcdfbu3euwL5of6HvffffZa7r78fHNN9/olltu0Zo1a1zOb8sxkZaWpnfeeUdvvfWWPv74Y3399de666671NTUZNZmdKrW9snZs2e1f/9+Pffcc9q/f7+2bdumv//975oyZYpT7Zw5cxyOnd/97ndmDL9LXOo4kS79s9KTjhNJDvuisrJSmzZtksVi0fTp0x3qustx0pbPXNN+n3TJQyrQblVVVYYkp+eY3X333e4blMmWLFli3HLLLS7nXbhwwbBarcaKFSvsbd9++63h5+dnbNiwwaQRutdTTz1lDBkyxLhw4YJhGD3v+JBkvPPOO/b3bTkmTp8+bXh5eRlvvfWWvebkyZNGr169jPfff9+0sXeV7+8TV/785z8bkozjx4/b2+Lj442nnnqqawfnJq72yaV+VjhODOPuu+82br/9doe27nycfP8z18zfJ5zx+jdRW1srSerfv79De2FhoQICAnTddddpzpw5qqqqcsfwTHP06FEFBwcrPDxcDzzwgD7//HNJUnl5uWw2mxISEuy13t7eio+P1+7du901XNM0NDTo9ddf18yZMx2+cL2nHR/f1ZZjoqSkRI2NjQ41wcHBioiI6BHHjXTxd4vFYnH63tg33nhD/v7+uvHGGzV//vxue+a4WWs/Kz39OPnHP/6hHTt2aNasWU7zuutx8v3PXDN/n7j1K4NwkWEYSk9P16233qqIiAh7e2Jiou677z6FhoaqvLxczz33nG6//XaVlJQ4fRl3dxAdHa0tW7bouuuu0z/+8Q89//zzio2N1eHDh+1fav79LykPDAzU8ePH3TFcU7377rs6ffq0UlJS7G097fj4vrYcEzabTb1799bVV1/tVNO8fHf27bffauHChXrooYccvgD5pz/9qf07af/yl78oIyNDBw8edPoO2+7iUj8rPf04ee211+Tr66tp06Y5tHfX48TVZ66Zv08IXv8GnnjiCX366acO3ycpSUlJSfbXERERioqKUmhoqHbs2OH0A9IdJCYm2l/fdNNNiomJ0ZAhQ/Taa6/Zb4T97tke6eIP0PfbuqONGzcqMTFRwcHB9raedny0pCPHRE84bhobG/XAAw/owoULWrduncO8OXPm2F9HRERo6NChioqK0v79+zVy5Eizh9rlOvqz0hOOE0natGmTfvrTn8rHx8ehvbseJy195krm/D7hUqObPfnkk9q+fbs++ugjXXvtta3WBgUFKTQ0VEePHjVpdO51xRVX6KabbtLRo0ftf934/f+rqKqqcvo/lO7m+PHj2rlzp2bPnt1qXU87PtpyTFitVjU0NOirr75qsaY7amxs1P3336/y8nIVFBQ4nO1yZeTIkfLy8uoxx873f1Z66nEiSbt27dKRI0cu+ftF6h7HSUufuWb+PiF4uYlhGHriiSe0bds2ffjhhwoPD7/kMjU1NTpx4oSCgoJMGKH71dfXq6ysTEFBQfbT3d89xd3Q0KCioiLFxsa6cZRdb/PmzQoICNCkSZNaretpx0dbjonIyEh5eXk51FRWVuovf/lLtz1umkPX0aNHtXPnTg0YMOCSyxw+fFiNjY095tj5/s9KTzxOmm3cuFGRkZG65ZZbLll7OR8nl/rMNfX3yQ/6swB02M9//nPDz8/PKCwsNCorK+3T2bNnDcMwjDNnzhhPP/20sXv3bqO8vNz46KOPjJiYGOOaa64x6urq3Dz6rvH0008bhYWFxueff27s2bPHuOuuuwxfX1/j2LFjhmEYxooVKww/Pz9j27ZtxqFDh4wHH3zQCAoK6rb7wzAMo6mpyRg0aJCxYMECh/aecnycOXPGKC0tNUpLSw1JxqpVq4zS0lL7X+i15ZhITU01rr32WmPnzp3G/v37jdtvv9245ZZbjPPnz7trs36Q1vZJY2OjMWXKFOPaa681Dhw44PC7pb6+3jAMw/jss8+MX/7yl8bevXuN8vJyY8eOHcYNN9xgjBgxolvuk7b+rPSk46RZbW2t0bdvX2P9+vVOy3e34+RSn7mGYd7vE4KXm0hyOW3evNkwDMM4e/askZCQYAwcONDw8vIyBg0aZMyYMcOoqKhw78C7UFJSkhEUFGR4eXkZwcHBxrRp04zDhw/b51+4cMFYsmSJYbVaDW9vb2PMmDHGoUOH3DjirvfBBx8YkowjR444tPeU4+Ojjz5y+XMyY8YMwzDadkycO3fOeOKJJ4z+/fsbffr0Me66667Lej+1tk/Ky8tb/N3y0UcfGYZhGBUVFcaYMWOM/v37G7179zaGDBlizJ0716ipqXHvhv0Are2Ttv6s9KTjpNnvfvc7o0+fPsbp06edlu9ux8mlPnMNw7zfJ5b/f0AAAADoYtzjBQAAYBKCFwAAgEkIXgAAACYheAEAAJiE4AUAAGASghcAAIBJCF4AAAAmIXgBAACYhOAFAABgEoIXgG4hJSVFFotFqampTvMee+wxWSwWpaSkOC2zcOFCSZLFYpHFYtGePXscaurr6zVgwABZLBYVFhZ21fAB9BAELwDdRkhIiN566y2dO3fO3vbtt99q69atGjRokEPthQsXtGPHDt19990Oy2/evNmh7p133tGVV17ZtQMH0GMQvAB0GyNHjtSgQYO0bds2e9u2bdsUEhKiESNGONT+6U9/Uq9evRQdHW1vmzFjhlNw27Rpk2bMmNHmMTQ0NOiJJ55QUFCQfHx8FBYWpszMTEnSsWPHZLFYdODAAXv96dOnHc6mFRYWymKx6P/+7/8UFRWlvn37KjY2VkeOHGnPrgDwb4rgBaBbefTRRx3OWm3atEkzZ850qtu+fbsmT56sXr3+9WswMjJS4eHhys3NlSSdOHFCf/zjH5WcnNzm9a9evVrbt2/Xf//3f+vIkSN6/fXXFRYW1u7tWLx4sV588UXt27dPnp6eLrcBwOWH4AWgW0lOTtbHH3+sY8eO6fjx4/rTn/6khx9+2Klu+/btDpcZmz366KPatGmTJGnz5s2aOHGiBg4c2Ob1V1RUaOjQobr11lsVGhqqW2+9VQ8++GC7t+M///M/FR8fr+HDh2vhwoXavXu3vv3223b3A+DfC8ELQLfi7++vSZMm6bXXXtPmzZs1adIk+fv7O9SUlZXpiy++0J133um0/MMPP6zi4mJ9/vnnevXVV9t9piklJUUHDhzQ9ddfr7lz5yo/P79D23HzzTfbXwcFBUmSqqqqOtQXgH8fBC8A3c7MmTP16quv6rXXXmvxMuO4cePUp08fp3kDBgzQXXfdpVmzZunbb79VYmJiu9Y9cuRIlZeX61e/+pXOnTun+++/X/fee68k2S9rGoZhr29sbHTZj5eXl/21xWKRdPEPAgBc3gheALqdCRMmqKGhQQ0NDRo/frzT/N///veaMmVKi8vPnDlThYWFeuSRR+Th4dHu9ffr109JSUl6+eWXlZOTo9zcXH355Zf2S5aVlZX22u/eaA+g+/N09wAAoLN5eHiorKzM/vq7qqqqtHfvXr377rstLj9hwgT985//VL9+/dq97t/85jcKCgrSj3/8Y/Xq1Utvv/22rFarrrrqKvXq1UujR4/WihUrFBYWpurqaj377LPtXgeAyxdnvAB0S/369XMZnP73f/9X0dHRCggIaHFZi8Uif39/9e7du93rvfLKK7Vy5UpFRUVp1KhROnbsmPLy8uyXGTdt2qTGxkZFRUXpqaee0vPPP9/udQC4fFmM795sAADd3JQpU3TrrbfqmWeecfdQAPRAnPEC0KN09PEOANAZOOMFAO2wfPlyLV++3OW8uLg4/eEPfzB5RAAuJwQvAGiHL7/8Ul9++aXLeX369NE111xj8ogAXE4IXgAAACbhHi8AAACTELwAAABMQvACAAAwCcELAADAJAQvAAAAkxC8AAAATELwAgAAMMn/BwCTuOIeDbXmAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import emcee\n",
    "import matplotlib.pyplot as plt\n",
    "alpha = 2.35\n",
    "M_min = 0.08\n",
    "M_max_limit = 200\n",
    "M_cloud = 1000\n",
    "def log_prob(M):\n",
    "    if M > M_min:\n",
    "        return -2.35 * np.log(M)\n",
    "    return -np.inf\n",
    "\n",
    "nwalkers = 16\n",
    "ndim = 1\n",
    "nsteps = 10000\n",
    "initial_positions = M_min + (M_max_limit - M_min) * np.random.rand(nwalkers, ndim)\n",
    "sampler = emcee.EnsembleSampler(nwalkers, ndim, log_prob)\n",
    "sampler.run_mcmc(initial_positions, nsteps, progress=True)\n",
    "samples = sampler.get_chain(flat=True)\n",
    "cumulative_mass = np.cumsum(samples)\n",
    "valid_samples = samples[cumulative_mass <= M_cloud]\n",
    "max_mass = np.max(valid_samples)\n",
    "print(f\"Maximal mass: {max_mass:.2f} M_sun\")\n",
    "plt.hist(valid_samples, bins=10, density=True)\n",
    "plt.xlabel(\"M/M_sun\")\n",
    "plt.ylabel(\"Probability\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Mtot(M>M*)/Mtot(M<M*)v= M*/(1000-M*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "太阳中的氢大概有10%在主序阶段被燃烧，请估算太阳处于主序阶段的时间？太阳表面的温度大概是~5500K，请由此估算地球表面的温度。地球的反射率~0.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "太阳处于主序阶段的时间为：0.725 Gyr\n"
     ]
    }
   ],
   "source": [
    "from astropy.constants import M_sun, L_sun, c, au, R_earth, sigma_sb\n",
    "M_H = 0.7346 * M_sun\n",
    "M_H_burn = 0.1 * M_H\n",
    "# 每千克氢核聚变释放的能量\n",
    "energy_per_kg = 6e13 * u.J / u.kg\n",
    "E_total = M_H_burn * energy_per_kg\n",
    "t_main_seq = E_total / L_sun\n",
    "print(f\"太阳处于主序阶段的时间为：{t_main_seq.to(u.Gyr):.3f}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Etot 和energy_per_kg可以自己算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "太阳辐射到地球的能量为$(1-0.3)*\\frac{L_{\\odot}}{4\\pi d^2} \\pi R^2$，地球到太阳的距离为1au\n",
    "\n",
    "和地球辐射能量为$4\\pi R^2 \\sigma T_{\\text{earth}}^4$相等\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "地球的温度约为: 254.59 K\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from astropy import units as u\n",
    "from astropy.constants import L_sun, au, sigma_sb\n",
    "A = 0.3\n",
    "d = au\n",
    "T_earth = (((1 - A) * L_sun) / (16 * np.pi * sigma_sb * d**2))**0.25\n",
    "print(f\"地球的温度约为: {T_earth:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 可以直接从太阳的温度和视面积计算"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "stylegan3",
   "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}