{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20b98fc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#作业11\n",
    "#假设在红移为1处去观测一个现在红移为2的星系，其红移为多少？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "84090167",
   "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": "code",
   "execution_count": 4,
   "id": "77d2ab68",
   "metadata": {},
   "outputs": [],
   "source": [
    "#假设宇宙的膨胀是一个线性过程，请根据今天的\n",
    "#哈勃常数估算宇宙的年龄。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0e70f4bf",
   "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": "code",
   "execution_count": 6,
   "id": "c0c56d3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#作业12\n",
    "#一团1000个太阳质量的中性气体云（全部为H原子）, 温度为30K，当云的密度大于多少时，气体云将发生塌缩？塌缩（自由下落）时标为多少？\n",
    "\n",
    "#分子云坍缩临界条件：𝐺𝑀𝑚𝐻𝑅=3𝐾𝑇2\n",
    "\n",
    "#密度𝑛\n",
    "#与半径𝑅\n",
    "#的关系：𝑛4𝜋𝑅33𝑚𝐻=1000𝑀⊙"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d13c034a",
   "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": "code",
   "execution_count": 9,
   "id": "1d11a917",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHHCAYAAABDUnkqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAV8dJREFUeJzt3XdYFFfbBvB7d4GlSRGkaGg2sPeC3UiCgho1dhMpxsSIlZjEjiVqLLHG8sZXweS1K3ajiSgalRhLMCYq0Qjip4JiAekK5/uDsLrSdnFxZLl/1zXXNXvmzM6zOyIP55w5RyaEECAiIiLSE3KpAyAiIiLSJSY3REREpFeY3BAREZFeYXJDREREeoXJDREREekVJjdERESkV5jcEBERkV5hckNERER6hckNERER6RUmN0Q60KlTJ3Tq1EnqMIjeGP7+/nB1dZU6DKqgmNxQhXTp0iX07dsXLi4uMDY2RrVq1fDOO+9gxYoVUodWpNOnT2PGjBl4/Pjxa7vmjBkzIJPJkJSUpCrz9/eHTCaDhYUFMjIyCpxz7do1yGQyyGQyLFq0SFUeGRmpKn95GzhwoE7j1vT+zp07F7t379bptYsTFxdX5HfQunXr1xaHrty5cwczZsxAdHS01KEQqTGQOgCi1+306dPo3LkznJ2dMXz4cDg4OODWrVv49ddfsWzZMowePVrqEAt1+vRpzJw5E/7+/rCyspI0FgMDA6Snp2Pfvn3o37+/2rGNGzfC2NgYmZmZhZ47ZswYtGjRQq1Ml3/ha3N/586di759+6JXr146u74mBg0aBB8fH7WyKlWqvNYYdOHOnTuYOXMmXF1d0bhxY7Vja9euRW5urjSBUYXH5IYqnDlz5sDS0hJnz54tkCTcu3dPmqAklJ6eDlNTU63OUSqVaNu2LTZv3lwgudm0aRN8fX2xc+fOQs9t3749+vbtW+p4SyL1/U1LS4OZmVmxdZo2bYoPPvhA59fOzMyEkZER5HLpG+UNDQ2lDoEqMOl/Aohes3/++Qf16tUrtPXDzs5O7XVoaCjefvtt2NnZQalUom7duli9enWJ18jvgtm6dSsmT54MBwcHmJmZoWfPnrh161aB+mfOnEHXrl1haWkJU1NTdOzYEadOnVIdnzFjBj7//HMAgJubm6orIy4uTlXnf//7H5o1awYTExNUrlwZAwcOLHCtTp06oX79+jh//jw6dOgAU1NTTJ48ucTPU5jBgwfjxx9/VOsmO3v2LK5du4bBgweX6j11QdP7K5PJkJaWhg0bNqi+T39/fwDAzZs3MXLkSLi7u8PExAQ2Njbo16+f2vcNAGFhYZDJZDh+/DhGjhwJOzs7vPXWW6/8GW7cuIF+/fqhcuXKMDU1RevWrXHgwAG1Ovn/xrZs2YKpU6eiWrVqMDU1RUpKCvz9/WFubo74+Hh0794d5ubmqFatGlauXAkgr9vu7bffhpmZGVxcXLBp0ya193748CEmTJiABg0awNzcHBYWFujWrRsuXryodv38FriAgADVdxgWFgag8DE3aWlp+Oyzz+Dk5ASlUgl3d3csWrQIQgi1ejKZDKNGjcLu3btRv359KJVK1KtXD4cOHXrl75YqBrbcUIXj4uKCqKgo/Pnnn6hfv36xdVevXo169eqhZ8+eMDAwwL59+zBy5Ejk5uYiKCioxGvNmTMHMpkMX375Je7du4elS5fCy8sL0dHRMDExAQAcPXoU3bp1Q7NmzRASEgK5XK5Kqn755Re0bNkSffr0wd9//43NmzdjyZIlsLW1BfC8K2POnDmYNm0a+vfvj48++gj379/HihUr0KFDB/z+++9qv+gfPHiAbt26YeDAgfjggw9gb29fqu+xT58+GDFiBMLDwxEYGAggr9XGw8MDTZs2LfK8J0+eqI3hAYDKlSvrrLVB0/v7ww8/4KOPPkLLli3x8ccfAwBq1KgBIC9JO336NAYOHIi33noLcXFxWL16NTp16oTLly8XaOkaOXIkqlSpgunTpyMtLa3EGNPT0wt8B5aWljA0NERiYiLatGmD9PR0jBkzBjY2NtiwYQN69uyJHTt2oHfv3mrnzZ49G0ZGRpgwYQKysrJgZGQEAMjJyUG3bt3QoUMHLFiwABs3bsSoUaNgZmaGKVOmYMiQIejTpw/WrFmDoUOHwtPTE25ubgDykqvdu3ejX79+cHNzQ2JiIv7zn/+gY8eOuHz5MqpWrYo6depg1qxZmD59Oj7++GO0b98eANCmTZtCP7MQAj179sSxY8cwbNgwNG7cGIcPH8bnn3+O27dvY8mSJWr1T548ifDwcIwcORKVKlXC8uXL8f777yM+Ph42NjYlfsdUwQmiCuann34SCoVCKBQK4enpKb744gtx+PBhkZ2dXaBuenp6gTJvb29RvXp1tbKOHTuKjh07ql4fO3ZMABDVqlUTKSkpqvJt27YJAGLZsmVCCCFyc3NFrVq1hLe3t8jNzVW7rpubm3jnnXdUZQsXLhQARGxsrNq14+LihEKhEHPmzFErv3TpkjAwMFAr79ixowAg1qxZU8w39FxISIgAIO7fv68q8/PzE2ZmZkIIIfr27Su6dOkihBAiJydHODg4iJkzZ4rY2FgBQCxcuLDAd1LY9vJnehXa3F8zMzPh5+dXoLyw+x4VFSUAiO+//15VFhoaKgCIdu3aiWfPnpUYW/73Uth27NgxIYQQ48aNEwDEL7/8ojrvyZMnws3NTbi6uoqcnBwhxPPvs3r16gXi9fPzEwDE3LlzVWWPHj0SJiYmQiaTiS1btqjKr169KgCIkJAQVVlmZqbqOi/GrlQqxaxZs1RlZ8+eFQBEaGhogc/q5+cnXFxcVK93794tAIivvvpKrV7fvn2FTCYT169fV5UBEEZGRmplFy9eFADEihUrClyL6GXslqIK55133kFUVBR69uyJixcvYsGCBfD29ka1atWwd+9etbr5rSsAkJycjKSkJHTs2BE3btxAcnJyidcaOnQoKlWqpHrdt29fODo64uDBgwCA6OhoVTfOgwcPkJSUhKSkJKSlpaFLly44ceJEiYMyw8PDkZubi/79+6vOT0pKgoODA2rVqoVjx46p1VcqlQgICCgxdk0MHjwYkZGRSEhIwNGjR5GQkFBil9T06dPx888/q20ODg46iQfQ7v4W5cX7/vTpUzx48AA1a9aElZUVLly4UKD+8OHDoVAoNI7x448/LvAdNGrUCABw8OBBtGzZEu3atVPVNzc3x8cff4y4uDhcvnxZ7b38/PzU4n3RRx99pNq3srKCu7s7zMzM1MZJubu7w8rKCjdu3FCVKZVKVUtaTk4OHjx4AHNzc7i7uxf6+TVx8OBBKBQKjBkzRq38s88+gxACP/74o1q5l5eXqiUNABo2bAgLCwu1OImKwm4pqpBatGiB8PBwZGdn4+LFi9i1axeWLFmCvn37Ijo6GnXr1gUAnDp1CiEhIYiKikJ6erraeyQnJ8PS0rLY69SqVUvttUwmQ82aNVVjN65duwYg7xdUUZKTk2FtbV3k8WvXrkEIUeBa+V4e2FmtWjVV18Wr8vHxQaVKlbB161ZER0ejRYsWap+vMA0aNICXl5fG18jIyCiQSJaUDGl6f4u75rx58xAaGorbt2+rjQkpLKnN787RVK1atYr8Dm7evIlWrVoVKK9Tp47q+IvdbUVd29jYuMATWJaWlnjrrbcgk8kKlD969Ej1Ojc3F8uWLcOqVasQGxuLnJwc1bHSdgndvHkTVatWVUv2X/5cL3J2di7wHtbW1mpxEhWFyQ1VaEZGRmjRogVatGiB2rVrIyAgANu3b0dISAj++ecfdOnSBR4eHli8eDGcnJxgZGSEgwcPYsmSJTp5zDX/PRYuXFjgUdp85ubmJb6HTCbDjz/+WGjrwcvnF/VXfmkolUr06dMHGzZswI0bNzBjxgydvXe+rVu3FmhpEi8NQC1Kcfe3OKNHj0ZoaCjGjRsHT09PWFpaqubjKey+6/I71VZR1y6qJamo8he/07lz52LatGkIDAzE7NmzVWOixo0b99oe79YkTqKiMLkh+lfz5s0BAHfv3gUA7Nu3D1lZWdi7d6/aX5Evd/MUJ79lJp8QAtevX0fDhg0BPB/AamFhUWJrxst/beerUaMGhBBwc3ND7dq1NY5NVwYPHoz169dDLpfrfDI+APD29sbPP//8yu/z8v0Fiv5Od+zYAT8/P3zzzTeqsszMzNcygaKLiwtiYmIKlF+9elV1vKzt2LEDnTt3xrp169TKHz9+rBrMDhT9/RXGxcUFR44cwZMnT9Rab17n56KKg2NuqMI5duxYoX/95Y+DcXd3B/D8L8eXuyRCQ0M1vtb333+PJ0+eqF7v2LEDd+/eRbdu3QAAzZo1Q40aNbBo0SKkpqYWOP/+/fuq/fy5U17+BdunTx8oFArMnDmzwOcSQuDBgwcax1sanTt3xuzZs/Htt9/qdOxMPkdHR3h5ealtxdH0/gJ532lhCYtCoSjwHitWrFDrnikrPj4++O233xAVFaUqS0tLw3fffQdXV9cSu9R0obDPv337dty+fVutrKh/k4Xx8fFBTk4Ovv32W7XyJUuWQCaTqX4miHSBLTdU4YwePRrp6eno3bs3PDw8kJ2djdOnT2Pr1q1wdXVVdYG8++67MDIyQo8ePfDJJ58gNTUVa9euhZ2dndpf/8WpXLky2rVrh4CAACQmJmLp0qWoWbMmhg8fDgCQy+X473//i27duqFevXoICAhAtWrVcPv2bRw7dgwWFhbYt28fgLxECACmTJmCgQMHwtDQED169ECNGjXw1VdfYdKkSYiLi0OvXr1QqVIlxMbGYteuXfj4448xYcKEMvgmofoMU6dOLbP315am9xfI+06PHDmCxYsXo2rVqnBzc0OrVq3QvXt3/PDDD7C0tETdunURFRWFI0eOvJZHkCdOnIjNmzejW7duGDNmDCpXrowNGzYgNjYWO3fufC0T9HXv3h2zZs1CQEAA2rRpg0uXLmHjxo2oXr26Wr0aNWrAysoKa9asQaVKlWBmZoZWrVoVOg6oR48e6Ny5M6ZMmYK4uDg0atQIP/30E/bs2YNx48apDR4memUSPKFFJKkff/xRBAYGCg8PD2Fubi6MjIxEzZo1xejRo0ViYqJa3b1794qGDRsKY2Nj4erqKubPny/Wr19f4PHloh4F37x5s5g0aZKws7MTJiYmwtfXV9y8ebNATL///rvo06ePsLGxEUqlUri4uIj+/fuLiIgItXqzZ88W1apVE3K5vEAMO3fuFO3atRNmZmbCzMxMeHh4iKCgIBETE6MWZ7169TT+rkp6FLwoxT0Kvn37do2vXxra3N+rV6+KDh06CBMTEwFA9Vj4o0ePREBAgLC1tRXm5ubC29tbXL16Vbi4uKg9Op7/KPjZs2c1iq2w76Uw//zzj+jbt6+wsrISxsbGomXLlmL//v1qdYr7Pou6R0XdfxcXF+Hr66t6nZmZKT777DPh6OgoTExMRNu2bUVUVFSBf+dCCLFnzx5Rt25dYWBgoPZY+MuPgguR90j7+PHjRdWqVYWhoaGoVauWWLhwodo0CELkPQoeFBRUaJyFPbpP9DKZEBydRaRrkZGR6Ny5M7Zv316mSw0QEVFBHHNDREREeoXJDREREekVJjdERESkVzjmhoiIiPQKW26IiIhIrzC5ISIiIr1S4Sbxy83NxZ07d1CpUiWtpg4nIiIi6Qgh8OTJE1StWrXEySwrXHJz584dODk5SR0GERERlcKtW7fw1ltvFVunwiU3+Qu23bp1CxYWFhJHQ0REpEfS0oCqVfP279wB/l1/TBdSUlLg5OSktvBqUSpccpPfFWVhYcHkhoiISJf+XXAYAGBhodPkJp8mQ0o4oJiIiIj0CpMbIiIi0isVrluKiIiIyoiBAeDn93xfqjAkuzIRERHpF6USCAuTOgp2SxEREZF+YcsNERER6YYQQHp63r6pKSDRZLlsuSEiIiLdSE8HzM3ztvwkRwJMboiIiEivMLkhIiIivcLkhoiIiPQKkxsiIiLSK0xuiIiISK8wuSEiIiK9wnluiIiISDcUCqBv3+f7EmFyQ0RERLphbAxs3y51FExudG1Y2Fm11+v8W0gUCRERUcXEMTdERESkV5jcEBERkW6kpeWtJyWT5e1LhMkNERER6RUmN0RERKRXmNwQERGRXpE0uTlx4gR69OiBqlWrQiaTYffu3SWeExkZiaZNm0KpVKJmzZoICwsr8ziJiIio/JA0uUlLS0OjRo2wcuVKjerHxsbC19cXnTt3RnR0NMaNG4ePPvoIhw8fLuNIiYiIqLyQdJ6bbt26oVu3bhrXX7NmDdzc3PDNN98AAOrUqYOTJ09iyZIl8Pb2LqswiYiIqBwpV2NuoqKi4OXlpVbm7e2NqKgoiSIiIiIiFYUC8PHJ27j8gmYSEhJgb2+vVmZvb4+UlBRkZGTAxMSkwDlZWVnIyspSvU5JSSnzOImIiCokY2PgwAGpoyhfLTelMW/ePFhaWqo2JycnqUMiIiKiMlSukhsHBwckJiaqlSUmJsLCwqLQVhsAmDRpEpKTk1XbrVu3XkeoREREJJFy1S3l6emJgwcPqpX9/PPP8PT0LPIcpVIJpVJZ1qERERFRWhpgZ5e3f+8eYGYmSRiSttykpqYiOjoa0dHRAPIe9Y6OjkZ8fDyAvFaXoUOHquqPGDECN27cwBdffIGrV69i1apV2LZtG8aPHy9F+ERERPSy9PS8TUKSJjfnzp1DkyZN0KRJEwBAcHAwmjRpgunTpwMA7t69q0p0AMDNzQ0HDhzAzz//jEaNGuGbb77Bf//7Xz4GTkRERCqSdkt16tQJQogijxc2+3CnTp3w+++/l2FUREREVJ6VqwHFRERERCVhckNERER6hckNERER6ZVy9Sg4ERERvcHkcqBjx+f7EmFyQ0RERLphYgJERkodBbuliIiISL8wuSEiIiK9wuSGiIiIdCMtDahSJW9LS5MsDI65ISIiIt1JSpI6ArbcEBERkX5hckNERER6hckNERER6RUmN0RERKRXmNwQERGRXuHTUkRERKQbcjnQvPnzfYkwuSEiIiLdMDEBzp6VOgp2SxEREZF+YXJDREREeoXJDREREelGejrg6pq3padLFgbH3BAREZFuCAHcvPl8XyJsuSEiIiK9wuSGiIiI9AqTGyIiItIrTG6IiIhIrzC5ISIiIr3Cp6WIiIhIN2QyoG7d5/sSYXJDREREumFqCvz1l9RRsFuKiIiI9AuTGyIiItIrTG6IiIhIN9LTgXr18jYuv0BERETlnhDA5cvP9yXClhsiIiLSK0xuiIiISK8wuSEiIiK9wuSGiIiI9AqTGyIiItIrfFqKiIiIdEMmA1xcnu9LhMkNERER6YapKRAXJ3UU7JYiIiIi/cLkhoiIiPQKkxsiIiLSjYwMoEWLvC0jQ7IwOOaGiIiIdCM3Fzh37vm+RNhyQ0RERHqFyQ0RERHpFSY3REREpFeY3BAREZFeYXJDREREeoVPSxEREZHu2NpKHQGTGyIiItIRMzPg/n2po2C3FBEREekXJjdERESkV5jcEBERkW5kZACdOuVtXH6BiIiIyr3cXOD48ef7EmHLDREREekVJjdERESkV5jcEBERkV5hckNERER6hckNERER6RXJk5uVK1fC1dUVxsbGaNWqFX777bdi6y9duhTu7u4wMTGBk5MTxo8fj8zMzNcULRERERXL1DRvk5Ckyc3WrVsRHByMkJAQXLhwAY0aNYK3tzfu3btXaP1NmzZh4sSJCAkJwZUrV7Bu3Tps3boVkydPfs2RExERUQFmZkBaWt5mZiZZGJImN4sXL8bw4cMREBCAunXrYs2aNTA1NcX69esLrX/69Gm0bdsWgwcPhqurK959910MGjSoxNYeIiIiqjgkS26ys7Nx/vx5eHl5PQ9GLoeXlxeioqIKPadNmzY4f/68Kpm5ceMGDh48CB8fnyKvk5WVhZSUFLWNiIiI9JdkMxQnJSUhJycH9vb2auX29va4evVqoecMHjwYSUlJaNeuHYQQePbsGUaMGFFst9S8efMwc+ZMncZOREREhcjMBN5/P29/507A2FiSMCQfUKyNyMhIzJ07F6tWrcKFCxcQHh6OAwcOYPbs2UWeM2nSJCQnJ6u2W7duvcaIiYiIKpCcHODgwbwtJ0eyMCRrubG1tYVCoUBiYqJaeWJiIhwcHAo9Z9q0afjwww/x0UcfAQAaNGiAtLQ0fPzxx5gyZQrk8oK5mlKphFKp1P0HICIiojeSZC03RkZGaNasGSIiIlRlubm5iIiIgKenZ6HnpKenF0hgFAoFAEAIUXbBEhERUbkh6argwcHB8PPzQ/PmzdGyZUssXboUaWlpCAgIAAAMHToU1apVw7x58wAAPXr0wOLFi9GkSRO0atUK169fx7Rp09CjRw9VkkNEREQVm6TJzYABA3D//n1Mnz4dCQkJaNy4MQ4dOqQaZBwfH6/WUjN16lTIZDJMnToVt2/fRpUqVdCjRw/MmTNHqo9AREREbxiZqGD9OSkpKbC0tERycjIsLCx0/v7Dws6qvV7n30Ln1yAiInojpaUB5uZ5+6mpOp3IT5vf3+XqaSkiIiKikkjaLUVERER6xMwMeAM6hNhyQ0RERHqFyQ0RERHpFSY3REREpBuZmUC/fnlbZqZkYTC5ISIiIt3IyQF27MjbJFx+gckNERER6RUmN0RERKRXmNwQERGRXmFyQ0RERHqFyQ0RERHpFSY3REREpFe4/AIRERHphqlp3oKZ+fsSYXJDREREuiGT6XQl8NJitxQRERHpFSY3REREpBtZWYC/f96WlSVZGExuiIiISDeePQM2bMjbnj2TLAytk5vQ0FCkp6eXRSxEREREr0zr5GbixIlwcHDAsGHDcPr06bKIiYiIiKjUtE5ubt++jQ0bNiApKQmdOnWCh4cH5s+fj4SEhLKIj4iIiEgrWic3BgYG6N27N/bs2YNbt25h+PDh2LhxI5ydndGzZ0/s2bMHubm5ZRErERERUYleaUCxvb092rVrB09PT8jlcly6dAl+fn6oUaMGIiMjdRQiERERkeZKldwkJiZi0aJFqFevHjp16oSUlBTs378fsbGxuH37Nvr37w8/Pz9dx0pERERUIq1nKO7RowcOHz6M2rVrY/jw4Rg6dCgqV66sOm5mZobPPvsMCxcu1GmgRERE9IYzNQXu3Xu+LxGtkxs7OzscP34cnp6eRdapUqUKYmNjXykwIiIiKmdkMqBKFamj0L5bqmPHjmjatGmB8uzsbHz//fcAAJlMBhcXl1ePjoiIiEhLWic3AQEBSE5OLlD+5MkTBAQE6CQoIiIiKoeysoCgoLytPC2/IISATCYrUP5///d/sLS01ElQREREVA49ewasWpW3Sbj8gsZjbpo0aQKZTAaZTIYuXbrAwOD5qTk5OYiNjUXXrl3LJEgiIiIiTWmc3PTq1QsAEB0dDW9vb5ibm6uOGRkZwdXVFe+//77OAyQiIiLShsbJTUhICADA1dUVAwYMgLGxcZkFRURERFRaWj8Kzsn5iIiI6E2mUXJTuXJl/P3337C1tYW1tXWhA4rzPXz4UGfBEREREWlLo+RmyZIlqFSpkmq/uOSGiIiISEoaJTcvdkX5+/uXVSxERERUnpmYAPkrFJiYSBaG1vPcXLhwAZcuXVK93rNnD3r16oXJkycjOztbp8ERERFROSKXA66ueZu8VGtz6yYMbU/45JNP8PfffwMAbty4gQEDBsDU1BTbt2/HF198ofMAiYiIiLShdXLz999/o3HjxgCA7du3o2PHjti0aRPCwsKwc+dOXcdHRERE5UV2NvD553mbhL05pVp+ITc3FwBw5MgR+Pj4AACcnJyQlJSk2+iIiIio/Hj6FFi0KG97+lSyMLRObpo3b46vvvoKP/zwA44fPw5fX18AQGxsLOzt7XUeIBEREZE2tE5uli5digsXLmDUqFGYMmUKatasCQDYsWMH2rRpo/MAiYiIiLSh9QzFDRs2VHtaKt/ChQuhUCh0EhQRERFRaWmd3OTLzs7GvXv3VONv8jk7O79yUERERESlpXVy8/fff2PYsGE4ffq0WrkQAjKZDDk5OToLjoiIiEhbWic3AQEBMDAwwP79++Ho6MilGIiIiOiNonVyEx0djfPnz8PDw6Ms4iEiIqLyysQE+PPP5/sS0Tq5qVu3LuezISIiooLkcqBePamj0P5R8Pnz5+OLL75AZGQkHjx4gJSUFLWNiIiISEpat9x4eXkBALp06aJWzgHFREREFVx2NjB3bt7+5MmAkZEkYWid3Bw7dqws4iAiIqLy7ulTYObMvP3PPy8/yU3Hjh3LIg4iIiIindB6zA0A/PLLL/jggw/Qpk0b3L59GwDwww8/4OTJkzoNjoiIiEhbWic3O3fuhLe3N0xMTHDhwgVkZWUBAJKTkzE3v5+NiIiISCJaJzdfffUV1qxZg7Vr18LQ0FBV3rZtW1y4cEGnwRERERFpS+vkJiYmBh06dChQbmlpicePH+siJiIiIqJS0zq5cXBwwPXr1wuUnzx5EtWrV9dJUERERESlpfXTUsOHD8fYsWOxfv16yGQy3LlzB1FRUZgwYQKmTZtWFjESERFReWBsDPz22/N9iWjdcjNx4kQMHjwYXbp0QWpqKjp06ICPPvoIn3zyCUaPHq11ACtXroSrqyuMjY3RqlUr/Jb/pRTh8ePHCAoKgqOjI5RKJWrXro2DBw9qfV0iIiLSMYUCaNEib1MoJAtD65YbmUyGKVOm4PPPP8f169eRmpqKunXrwtzcXOuLb926FcHBwVizZg1atWqFpUuXwtvbGzExMbCzsytQPzs7G++88w7s7OywY8cOVKtWDTdv3oSVlZXW1yYiIiL9pHVyA+QttZCSkgJ7e3vUrVu31BdfvHgxhg8fjoCAAADAmjVrcODAAaxfvx4TJ04sUH/9+vV4+PAhTp8+rXpSy9XVtdTXJyIiIh3KzgaWLcvbHztWshmKteqWSkhIwNChQ2FtbQ17e3vY2dnB2toagYGBSExM1OrC2dnZOH/+vGqtKgCQy+Xw8vJCVFRUoefs3bsXnp6eCAoKgr29PerXr4+5c+cWu55VVlYWF/ckIiJ6HZ4+Bb74Im97+lSyMDRuuUlJSUGbNm2QmpqKgIAAeHh4QAiBy5cvY/PmzTh58iQuXLigcfdUUlIScnJyYG9vr1Zub2+Pq1evFnrOjRs3cPToUQwZMgQHDx7E9evXMXLkSDx9+hQhISGFnjNv3jzMzF/ngoiIiPSexsnNsmXLoFAo8Ndff6FKlSpqx6ZOnYq2bdti+fLlmDx5ss6DzJebmws7Ozt89913UCgUaNasGW7fvo2FCxcWmdxMmjQJwcHBqtcpKSlwcnIqsxiJiIhIWhp3Sx04cACTJ08ukNgAgJ2dHSZNmoR9+/ZpfGFbW1soFIoC3VmJiYlwcHAo9BxHR0fUrl0bihdGYNepUwcJCQnIzs4u9BylUgkLCwu1jYiIiPSXxsnN33//jTZt2hR5vE2bNoiJidH4wkZGRmjWrBkiIiJUZbm5uYiIiICnp2eh57Rt2xbXr19Hbm6uWlyOjo4wkmjQEhEREb1ZNE5uUlJSin3k2srKSuvBusHBwVi7di02bNiAK1eu4NNPP0VaWprq6amhQ4di0qRJqvqffvopHj58iLFjx+Lvv//GgQMHMHfuXAQFBWl1XSIiItJfGo+5EUJALi86F5LJZBBCaHXxAQMG4P79+5g+fToSEhLQuHFjHDp0SDXIOD4+Xu2aTk5OOHz4MMaPH4+GDRuiWrVqGDt2LL788kutrktERET6SyY0zEjkcjksLS0hk8kKPZ4/901xj2W/CVJSUmBpaYnk5OQyGX8zLOys2ut1/i10fg0iIqI3Uk4O8Msvefvt2+t0lmJtfn9r3HITGhr6yoERERGRHlMogE6dpI5C8+TGz8+vLOMgIiIi0olSLb9AREREVMDTp8B33+Xtf/wx8O9SSa8bkxsiIiLSjexsYNSovH1/f8mSG63WliIiIiJ60zG5ISIiIr2idXJz7NixsoiDiIiISCe0Tm66du2KGjVq4KuvvsKtW7fKIiYiIiKiUtM6ubl9+zZGjRqFHTt2oHr16vD29sa2bduKXLiSiIiI6HXSOrmxtbXF+PHjER0djTNnzqB27doYOXIkqlatijFjxuDixYtlEScRERGRRl5pQHHTpk0xadIkjBo1CqmpqVi/fj2aNWuG9u3b46+//tJVjERERFQeKJXA/v15m1IpWRilSm6ePn2KHTt2wMfHBy4uLjh8+DC+/fZbJCYm4vr163BxcUG/fv10HSsRERG9yQwMAF/fvM1Auqn0tL7y6NGjsXnzZggh8OGHH2LBggWoX7++6riZmRkWLVqEqlWr6jRQIiIiIk1ondxcvnwZK1asQJ8+faAsosnJ1taWj4wTERFVNE+fAhs35u0PGVJ+ZigOCQlBv379CiQ2z549w4kTJwAABgYG6Nixo24iJCIiovIhOxsICMjbJHyKWuvkpnPnznj48GGB8uTkZHTu3FknQRERERGVltbJjRACMpmsQPmDBw9gZmamk6CIiIiISkvjMTd9+vQBAMhkMvj7+6t1S+Xk5OCPP/5AmzZtdB8hERERkRY0Tm4sLS0B5LXcVKpUCSYmJqpjRkZGaN26NYYPH677CImIiIi0oHFyExoaCgBwdXXFhAkT2AVFREREbyStHwUPCQkpiziIiIiIdEKj5KZp06aIiIiAtbU1mjRpUuiA4nwXLlzQWXBERERUjiiVwLZtz/clolFy895776kGEPfq1ass4yEiIqLyysAAeAOWX9IouXmxK4rdUkRERPQmk25VKyIiItIvz54Bu3bl7ffuLdnimRpd1drauthxNi8qbPZiIiIiqgCysoD+/fP2U1Pf7ORm6dKlZRwGERERkW5olNz4+fmVdRxEREREOqFRcpOSkgILCwvVfnHy6xERERFJQeMxN3fv3oWdnR2srKwKHX+Tv6BmTk6OzoMkIiIi0pRGyc3Ro0dRuXJlAMCxY8fKNCAiIiKiV6FRctOxY8dC94mIiIjeNKV6RuvRo0dYt24drly5AgCoW7cuAgICVK07REREVAEZGQH/LrQNIyPJwpBre8KJEyfg6uqK5cuX49GjR3j06BGWL18ONzc3nDhxoixiJCIiovLA0BDw98/bDA0lC0PrlpugoCAMGDAAq1evhkKhAADk5ORg5MiRCAoKwqVLl3QeJBEREZGmtG65uX79Oj777DNVYgMACoUCwcHBuH79uk6DIyIionLk2TPgwIG87dkzycLQuuWmadOmuHLlCtzd3dXKr1y5gkaNGuksMCIiIipnsrKA7t3z9t/05Rf++OMP1f6YMWMwduxYXL9+Ha1btwYA/Prrr1i5ciW+/vrrsomSiIiISEMyIYQoqZJcLodMJkNJVcvDJH4pKSmwtLREcnJymcymPCzsrNrrdf4tdH4NIiKiN1JaGmBunrefmgqYmensrbX5/a1Ry01sbKxOAiMiIiIqaxolNy4uLmUdBxEREZFOlHqkz+XLlxEfH4/s7Gy18p49e75yUERERESlpXVyc+PGDfTu3RuXLl1SG4eTv5jmmz7mhoiIiPSb1vPcjB07Fm5ubrh37x5MTU3x119/4cSJE2jevDkiIyPLIEQiIiIqF4yMgG+/zdskXH5B65abqKgoHD16FLa2tpDL5ZDL5WjXrh3mzZuHMWPG4Pfffy+LOImIiOhNZ2gIBAVJHYX2LTc5OTmoVKkSAMDW1hZ37twBkDfoOCYmRrfREREREWlJ65ab+vXr4+LFi3Bzc0OrVq2wYMECGBkZ4bvvvkP16tXLIkYiIiIqD3JygF9+ydtv3x54Yamm10nr5Gbq1KlIS0sDAMyaNQvdu3dH+/btYWNjg61bt+o8QCIiIionMjOBzp3z9nU8iZ82tE5uvL29Vfs1a9bE1atX8fDhQ1hbW6uemCIiIiKSyiutaHXr1i0AgJOTk06CISIiInpVWg8ofvbsGaZNmwZLS0u4urrC1dUVlpaWmDp1Kp4+fVoWMRIRERFpTOuWm9GjRyM8PBwLFiyAp6cngLzHw2fMmIEHDx5g9erVOg+SiIiISFNaJzebNm3Cli1b0K1bN1VZw4YN4eTkhEGDBjG5ISIiIklp3S2lVCrh6upaoNzNzQ1GEs5GSERERASUIrkZNWoUZs+ejaysLFVZVlYW5syZg1GjRuk0OCIiIipHDA2BBQvyNkNDycLQqFuqT58+aq+PHDmCt956C40aNQIAXLx4EdnZ2ejSpYvuIyQiIqLywcgI+PxzqaPQLLmxtLRUe/3++++rveaj4ERERPSm0Ci5CQ0NLdMgVq5ciYULFyIhIQGNGjXCihUr0LJlyxLP27JlCwYNGoT33nsPu3fvLtMYiYiIqAQ5OcCFC3n7TZtKtvyC1mNu8t2/fx8nT57EyZMncf/+/VIHsHXrVgQHByMkJAQXLlxAo0aN4O3tjXv37hV7XlxcHCZMmID27duX+tpERESkQ5mZQMuWeVtmpmRhaJ3cpKWlITAwEI6OjujQoQM6dOiAqlWrYtiwYUhPT9c6gMWLF2P48OEICAhA3bp1sWbNGpiammL9+vVFnpOTk4MhQ4Zg5syZXKyTiIiI1Gid3AQHB+P48ePYt28fHj9+jMePH2PPnj04fvw4PvvsM63eKzs7G+fPn4eXl9fzgORyeHl5ISoqqsjzZs2aBTs7OwwbNqzEa2RlZSElJUVtIyIiIv2l9SR+O3fuxI4dO9CpUydVmY+PD0xMTNC/f3+tJvFLSkpCTk4O7O3t1crt7e1x9erVQs85efIk1q1bh+joaI2uMW/ePMycOVPjmIiIiKh807rlJj09vUAyAgB2dnal6pbSxpMnT/Dhhx9i7dq1sLW11eicSZMmITk5WbXlL/ZJRERE+knrlhtPT0+EhITg+++/h7GxMQAgIyMDM2fOVK01pSlbW1soFAokJiaqlScmJsLBwaFA/X/++QdxcXHo0aOHqiw3NzfvgxgYICYmBjVq1FA7R6lUQqlUahUXERERlV9aJzdLly5F165dC0ziZ2xsjMOHD2v1XkZGRmjWrBkiIiLQq1cvAHnJSkRERKGzHXt4eODSpUtqZVOnTsWTJ0+wbNkyzrdDRERE2ic3DRo0wLVr17Bx40bVuJhBgwZhyJAhMDEx0TqA4OBg+Pn5oXnz5mjZsiWWLl2KtLQ0BAQEAACGDh2KatWqYd68eTA2Nkb9+vXVzreysgKAAuVERET0mhkaAiEhz/clolVy8/TpU3h4eGD//v0YPny4TgIYMGAA7t+/j+nTpyMhIQGNGzfGoUOHVON64uPjIZeXejoeIiIiel2MjIAZM6SOQrvkxtDQEJllMCnPqFGjilx0MzIysthzw8LCdB4PERERlV9aN4kEBQVh/vz5ePbsWVnEQ0REROVVbi7w1195278P/EhB6zE3Z8+eRUREBH766Sc0aNAAZmZmasfDw8N1FhwRERGVIxkZQP4Y2NRU4KUc4XXROrmxsrIqsCo4ERER0ZtC6+SmrFcIJyIiInoVGo+5yc3Nxfz589G2bVu0aNECEydOREZGRlnGRkRERKQ1jZObOXPmYPLkyTA3N0e1atWwbNkyBAUFlWVsRERERFrTOLn5/vvvsWrVKhw+fBi7d+/Gvn37sHHjRtXyB0RERERvAo2Tm/j4ePj4+Khee3l5QSaT4c6dO2USGBEREVFpaDyg+NmzZ6qFMvMZGhri6dOnOg+KiIiIyiFDQ2DChOf7EtE4uRFCwN/fX22F7czMTIwYMUJtrhvOc0NERFRBGRkBCxdKHYXmyY2fn1+Bsg8++ECnwRARERG9Ko2TG85vQ0RERMXKzQXi4/P2nZ0BiRa+1noSPyIiIqJCZWQAbm55+xIuvyBNSkVERERURpjcEBERkV5hckNERER6hckNERER6RUmN0RERKRXmNwQERGRXuGj4ERERKQbBgbAyJHP96UKQ7IrExERkX5RKoGVK6WOgt1SREREpF/YckNERES6IQSQlJS3b2sLyGSShMHkhoiIiHQjPR2ws8vb5/ILRERERLrB5IaIiIj0CpMbIiIi0itMboiIiEivMLkhIiIivcLkhoiIiPQKHwUnIiIi3TAwAPz8nu9LFYZkVyYiIiL9olQCYWFSR8FuKSIiItIvbLkhIiIi3RAib5ZiADA1lWz5BbbcEBERkW6kpwPm5nlbfpIjASY3REREpFeY3BAREZFeYXJDREREeoXJDREREekVJjdERESkV5jcEBERkV7hPDdERESkGwoF0Lfv832JMLkhIiIi3TA2BrZvlzoKdksRERGRfmFyQ0RERHqFyQ0RERHpRlpa3npSMlnevkSY3BAREZFeYXJDREREeoXJDREREekVJjdERESkV5jcEBERkV5hckNERER6hTMUExERkW4oFICPz/N9iTC5ISIiIt0wNgYOHJA6CnZLERERkX5hckNERER6hckNERER6UZaGmBmlrdV9OUXVq5cCVdXVxgbG6NVq1b47bffiqy7du1atG/fHtbW1rC2toaXl1ex9YmIiOg1Sk/P2yQkeXKzdetWBAcHIyQkBBcuXECjRo3g7e2Ne/fuFVo/MjISgwYNwrFjxxAVFQUnJye8++67uH379muOnIiIiN5EMiGEkDKAVq1aoUWLFvj2228BALm5uXBycsLo0aMxceLEEs/PycmBtbU1vv32WwwdOrTE+ikpKbC0tERycjIsLCxeOf6XDQs7q/Z6nX8LnV+DiIjojZSWBpib5+2npuZ1T+mINr+/JW25yc7Oxvnz5+Hl5aUqk8vl8PLyQlRUlEbvkZ6ejqdPn6Jy5cqFHs/KykJKSoraRkRERPpL0uQmKSkJOTk5sLe3Vyu3t7dHQkKCRu/x5ZdfomrVqmoJ0ovmzZsHS0tL1ebk5PTKcRMREdGbS/IxN6/i66+/xpYtW7Br1y4YGxsXWmfSpElITk5Wbbdu3XrNURIREdHrJOkMxba2tlAoFEhMTFQrT0xMhIODQ7HnLlq0CF9//TWOHDmChg0bFllPqVRCqVTqJF4iIiIqhlwOdOz4fF+qMCS7MgAjIyM0a9YMERERqrLc3FxERETA09OzyPMWLFiA2bNn49ChQ2jevPnrCJWIiIhKYmICREbmbSYmkoUh+dpSwcHB8PPzQ/PmzdGyZUssXboUaWlpCAgIAAAMHToU1apVw7x58wAA8+fPx/Tp07Fp0ya4urqqxuaYm5vDPH+ENhEREVVYkic3AwYMwP379zF9+nQkJCSgcePGOHTokGqQcXx8POQvNG2tXr0a2dnZ6Nu3r9r7hISEYMaMGa8zdCIiInoDST7PzevGeW6IiIjKSFoa4Oqatx8XJ9k8N5K33BAREZEeSUqSOoLy/Sg4ERER0cuY3BAREZFeYXJDREREeoXJDREREekVJjdERESkV/i0FBEREemGXA7krxwg4fILTG6ICpGTk4OnT59KHQZRhWNkZKQ2cSuVMyYmwNmzJdcrY0xuiF4ghEBCQgIeP34sdShEFZJcLoebmxuMjIykDoXKMSY3RC/IT2zs7OxgamoKmUwmdUhEFUZubi7u3LmDu3fvwtnZmT9/VGpMboj+lZOTo0psbGxspA6HqEKqUqUK7ty5g2fPnsHQ0FDqcEhb6elA3bp5+5cvA6amkoTB5IboX/ljbEwl+mEkIqi6o3JycpjclEdCADdvPt+XCEdtEb2ETeFE0uHPH+kCkxsiIiLSK0xuiKhM+fv7o1evXjp/37CwMFhZWalez5gxA40bN9b5dQq7FhG92TjmhkgDw8Je37wN6/xbaH2Ov78/Hj9+jN27d6teb9iwAZ988gnWrFmjVjcoKAirVq2Cn58fwsLC1Oq/7Nq1a6hZs2aB8sjISHTu3BlAXjdCpUqVUL16dbzzzjsYP348HB0dVXWXLVsGoWHf+8ufozgDBgyAj4+PRu+rDVdXV4wbNw7jxo0r82sRUdlgyw2RnnJycsKWLVuQkZGhKsvMzMSmTZvg7OxcoH7Xrl1x9+5dtc3Nza3Ya8TExODOnTs4e/YsvvzySxw5cgT169fHpUuXVHUsLS113urx9OlTmJiYwM7OTqfvW5TXeS0ienVMboj0VNOmTeHk5ITw8HBVWXh4OJydndGkSZMC9ZVKJRwcHNQ2hUJR7DXs7Ozg4OCA2rVrY+DAgTh16hSqVKmCTz/9VFXn5W6pHTt2oEGDBjAxMYGNjQ28vLyQlpaGGTNmYMOGDdizZw9kMhlkMhkiIyMRFxcHmUyGrVu3omPHjjA2NsbGjRuL7Cr6z3/+AycnJ5iamqJ///5ITk5WHevUqZNaiwwA9OrVC/7+/qrjN2/exPjx41UxAIV3S61evRo1atSAkZER3N3d8cMPP6gdl8lk+O9//4vevXvD1NQUtWrVwt69e4v9PonKPZks71HwunXz9iXC5IZIjwUGBiI0NFT1ev369QgICCiz65mYmGDEiBE4deoU7t27V+D43bt3MWjQIAQGBuLKlSuIjIxEnz59IITAhAkT0L9/f7UWpDZt2qjOnThxIsaOHYsrV67A29u70Otfv34d27Ztw759+3Do0CH8/vvvGDlypMbxh4eH46233sKsWbNUMRRm165dGDt2LD777DP8+eef+OSTTxAQEIBjx46p1Zs5cyb69++PP/74Az4+PhgyZAgePnyocTxE5Y6pKfDXX3mbhNNqMLkh0mMffPABTp48iZs3b+LmzZs4deoUPvjgg0Lr7t+/H+bm5qqtX79+pbqmh4cHACAuLq7Asbt37+LZs2fo06cPXF1d0aBBA4wcOVJ1TRMTE7UWpBen4B83bhz69OkDNzc3tTE9L8rMzMT333+Pxo0bo0OHDlixYgW2bNmChIQEjWKvXLkyFAoFKlWqpIqhMIsWLYK/vz9GjhyJ2rVrIzg4GH369MGiRYvU6vn7+2PQoEGoWbMm5s6di9TUVPz2228axUJEpccBxUR6rEqVKvD19UVYWBiEEPD19YWtrW2hdTt37ozVq1erXpuZmZXqmvmDhwubr6RRo0bo0qULGjRoAG9vb7z77rvo27cvrK2tS3zf5vkrDRfD2dkZ1apVU7329PREbm4uYmJiikxUSuPKlSv4+OOP1cratm2LZcuWqZU1bNhQtW9mZgYLC4tCW7SISLeY3BDpucDAQIwaNQoAsHLlyiLrmZmZFfpklLauXLkCIO+po5cpFAr8/PPPOH36NH766SesWLECU6ZMwZkzZ0ocvFzaZOtFcrm8wJNbZbn6+8sz7MpkMuTm5pbZ9Ygkl54OtPj3ic+zZyXrmmK3FJGe69q1K7Kzs/H06dMix6roSkZGBr777jt06NABVapUKbSOTCZD27ZtMXPmTPz+++8wMjLCrl27AORNvZ+Tk1Pq68fHx+POnTuq17/++ivkcjnc3d0B5LVkvTiOJicnB3/++afae2gSQ506dXDq1Cm1slOnTqFu/po6RBWVEHlrSl2+LOnyC2y5IdJzCoVC1ZpS0tNP2rp37x4yMzPx5MkTnD9/HgsWLEBSUpLaE1ovOnPmDCIiIvDuu+/Czs4OZ86cwf3791GnTh0Aea09hw8fRkxMDGxsbGBpaalVPMbGxvDz88OiRYuQkpKCMWPGoH///qouqbfffhvBwcE4cOAAatSogcWLF+Px48dq7+Hq6ooTJ05g4MCBUCqVhXbjff755+jfvz+aNGkCLy8v7Nu3D+Hh4Thy5IhW8RJR2WByQ1QBWFhYlMn7uru7QyaTwdzcHNWrV8e7776L4ODgIse3WFhY4MSJE1i6dClSUlLg4uKCb775Bt26dQMADB8+HJGRkWjevDlSU1Nx7NixQru3ilKzZk306dMHPj4+ePjwIbp3745Vq1apjgcGBuLixYsYOnQoDAwMMH78eNVkhPlmzZqFTz75BDVq1EBWVlahExD26tULy5Ytw6JFizB27Fi4ubkhNDQUnTp10jhWIio7MqHp1KF6IiUlBZaWlkhOTi6T//Bfnsm2NLPNkjQyMzMRGxsLNzc3GBsbSx0OUYXEn8NyLi0NMDfP209NBXQwVi6fNr+/OeaGiIiI9AqTGyIiItIrHHNDREREuiGTAS4uz/clwuSGiIiIdMPUFChkdvLXjd1SREREpFeY3BAREZFeYXJDREREupGRkbf8QosWefsS4ZgbIiIi0o3cXODcuef7EmHLDREREekVJjdERESkV5jcEOkBf39/yGQyjBgxosCxoKAgyGQy+Pv7l2kMYWFhkMlkqkUwX7R9+3bIZDKt1ol605w9exZdunSBlZUVrK2t4e3tjYsXL6rVOXz4MFq3bo1KlSqhSpUqeP/99xFXzGOxkZGRkMlkhW5nzz5fymXbtm1o3LgxTE1N4eLigoULF5YYb/77/Prrr2rlWVlZsLGxgUwmQ2RkpFbfAVF5weSGSE84OTlhy5YtyHhhEF9mZiY2bdoEZ2fn1xKDmZkZ7t27h6ioKLXydevWvbYYykJqaiq6du0KZ2dnnDlzBidPnkSlSpXg7e2Np0+fAgBiY2Px3nvv4e2330Z0dDQOHz6MpKQk9OnTp8j3bdOmDe7evau2ffTRR3Bzc0Pz5s0BAD/++COGDBmCESNG4M8//8SqVauwZMkSfPvttyXG7eTkhNDQULWyXbt2wTx/7R8iPcXkhkhPNG3aFE5OTggPD1eVhYeHw9nZGU2aNFGre+jQIbRr1w5WVlawsbFB9+7d8c8//6iOf//99zA3N8e1a9dUZSNHjoSHhwfS09OLjMHAwACDBw/G+vXrVWX/93//h8jISAwePFit7j///IP33nsP9vb2MDc3R4sWLXDkyBG1OqtWrUKtWrVgbGwMe3t79O3bV3Vsx44daNCgAUxMTGBjYwMvLy+kpaVp+G1p5+rVq3j48CFmzZoFd3d31KtXDyEhIUhMTMTNmzcBAOfPn0dOTg6++uor1KhRA02bNsWECRMQHR2tSoBeZmRkBAcHB9VmY2ODPXv2ICAgALJ/Z3f94Ycf0KtXL4wYMQLVq1eHr68vJk2ahPnz5xe6YvmL/Pz8CiS869evh5+fn8afPTs7G6NGjYKjoyOMjY3h4uKCefPmAQDi4uIgk8kQHR2tqv/48WO1VqH81qmIiAg0b94cpqamaNOmDWJiYjSOgUhbTG6INJGWVvSWmal53ZcfjSyszisIDAxU+0t9/fr1CAgIKOTjpCE4OBjnzp1DREQE5HI5evfujdx/n24YOnQofHx8MGTIEDx79gwHDhzAf//7X2zcuBGmpqYlxrBt2zZVEhQWFoauXbvC3t5erV5qaip8fHwQERGB33//HV27dkWPHj0QHx8PADh37hzGjBmDWbNmISYmBocOHUKHDh0AAHfv3sWgQYMQGBiIK1euIDIyEn369Cn2l725uXmxW2Fdevnc3d1hY2ODdevWITs7GxkZGVi3bh3q1Kmj6mpr1qwZ5HI5QkNDkZOTg+TkZPzwww/w8vKCoaFhsd9Zvr179+LBgwdq9ywrK6vA6tgmJib4v//7P1ViVZRmzZrB1dUVO3fuBADEx8fjxIkT+PDDDzWKBwCWL1+OvXv3Ytu2bYiJicHGjRtL1b04ZcoUfPPNNzh37hwMDAwQGBio9XtQOWFrm7dJSVQwycnJAoBITk4uk/cPDP1NbaPyIyMjQ1y+fFlkZGQUPAgUvfn4qNc1NS26bseO6nVtbQvWKQU/Pz/x3nvviXv37gmlUini4uJEXFycMDY2Fvfv3xfvvfee8PPzK/L8+/fvCwDi0qVLqrKHDx+Kt956S3z66afC3t5ezJkzp9gYQkNDhaWlpRBCiMaNG4sNGzaI3NxcUaNGDbFnzx6xZMkS4eLiUux71KtXT6xYsUIIIcTOnTuFhYWFSElJKVDv/PnzAoCIi4sr9v1edO3atWK3xMTEYs+/dOmSqFGjhpDL5UIulwt3d/cC14+MjBR2dnZCoVAIAMLT01M8evRI4xi7desmunXrplb2n//8R5iamoojR46InJwcERMTIzw8PAQAcfr06SLfC4DYtWuXWLp0qejcubMQQoiZM2eK3r17i0ePHgkA4tixYyXGNHr0aPH222+L3NzcAsdiY2MFAPH777+ryl5+72PHjgkA4siRI6o6Bw4cEAAK/Vkr9ueQKjRtfn+z5YZIj1SpUgW+vr4ICwtDaGgofH19YVvIX1DXrl3DoEGDUL16dVhYWKj+Es9vNQEAa2trrFu3DqtXr0aNGjUwceJEjePIb0E6fvw40tLS4OPjU6BOamoqJkyYgDp16sDKygrm5ua4cuWKKoZ33nkHLi4uqF69Oj788ENs3LhR1RrUqFEjdOnSBQ0aNEC/fv2wdu1aPHr0qNiYatasWexmZ2dX5LkZGRkYNmwY2rZti19//RWnTp1C/fr14evrq+rySUhIwPDhw+Hn54ezZ8/i+PHjMDIyQt++fUvsPgLyuu8OHz6MYcOGqZUPHz4co0aNQvfu3WFkZITWrVtj4MCBAAC5vOT/wj/44ANERUXhxo0bCAsL07rFxN/fH9HR0XB3d8eYMWPw008/aXV+voYNG6r2HR0dAQD37t0r1XsRlYTJDZEmUlOL3v5t8le5d6/ouj/+qF43Lq5gnVcUGBiIsLAwbNiwochfZD169MDDhw+xdu1anDlzBmfOnAGQN77iRSdOnIBCocDdu3e1Gs8yZMgQ/Prrr5gxYwY+/PBDGBgUnC90woQJ2LVrF+bOnYtffvkF0dHRaNCggSqGSpUq4cKFC9i8eTMcHR0xffp0NGrUCI8fP4ZCocDPP/+MH3/8EXXr1sWKFSvg7u6O2NjYImN6lW6pTZs2IS4uDqGhoWjRogVat26NTZs2ITY2Fnv27AEArFy5EpaWlliwYAGaNGmCDh064H//+x8iIiJU329xQkNDYWNjg549e6qVy2QyzJ8/H6mpqbh58yYSEhLQsmVLAED16tVLfN/8MVXDhg1DZmYmunXrVuI5L2ratCliY2Mxe/ZsZGRkoH///qqxT/nJ1YvJW1Hji17smssfT5Qr4SRvpN84Q3EZGxZ2tkDZOv8WEkRCr8TMTPq6GuratSuys7Mhk8ng7e1d4PiDBw8QExODtWvXon379gCAkydPFqh3+vRpzJ8/H/v27cOXX36JUaNGYcOGDRrFULlyZfTs2RPbtm3DmjVrCq1z6tQp+Pv7o3fv3gDyWnJefmzawMAAXl5e8PLyQkhICKysrHD06FH06dMHMpkMbdu2Rdu2bTF9+nS4uLhg165dCA4OLvR6Lw56LYyFhUWRx9LT0yGXy1W/lAGoXuf/gs6v8yKFQgGg5F/iQgiEhoZi6NChRY7PUSgUqFatGgBg8+bN8PT0RJUqVYp933yBgYHw8fHBl19+qYpJGxYWFhgwYAAGDBiAvn37omvXrnj48KHq+nfv3lUNWi/peyY9l5EB5CfQP/4ImJhIEgaTGyI9o1AocOXKFdX+y6ytrWFjY4PvvvsOjo6OiI+PL9Dl9OTJE3z44YcYM2YMunXrhrfeegstWrRAjx491J5YKk5YWBhWrVoFGxubQo/XqlUL4eHh6NGjB2QyGaZNm6aWBOzfvx83btxAhw4dYG1tjYMHDyI3Nxfu7u44c+YMIiIi8O6778LOzg5nzpzB/fv3C51jJ1/NmjU1irsw77zzDj7//HMEBQVh9OjRyM3Nxddffw0DAwN07twZAODr64slS5Zg1qxZGDRoEJ48eYLJkyfDxcVF9Yv/t99+w9ChQxEREaFKVADg6NGjiI2NxUcffVTg2klJSdixYwc6deqEzMxMhIaGYvv27Th+/LjG8Xft2hX3798vNoEryuLFi+Ho6IgmTZpALpdj+/btcHBwgJWVFeRyOVq3bo2vv/4abm5uuHfvHqZOnar1NUiP5OYC+f82ufwCEemShYVFkb/I5HI5tmzZgvPnz6N+/foYP358gUnhxo4dCzMzM8ydOxcA0KBBA8ydOxeffPIJbt++rVEM+Y9oF2Xx4sWwtrZGmzZt0KNHD3h7e6Np06aq41ZWVggPD8fbb7+NOnXqYM2aNdi8eTPq1asHCwsLnDhxAj4+PqhduzamTp2Kb775RusuF015eHhg3759+OOPP+Dp6Yn27dvjzp07OHTokGr8yNtvv41NmzZh9+7daNKkCbp27QqlUolDhw7B5N+/XtPT0xETE1Og62bdunVo06YNPDw8Cr3+hg0b0Lx5c7Rt2xZ//fUXIiMjVV1TmpDJZLC1tYWRkZHWn71SpUpYsGABmjdvjhYtWiAuLg4HDx5UtVKtX78ez549Q7NmzTBu3Dh89dVXWl+DSNdkQpORbnokJSUFlpaWSE5OLtVfMSUprBvqZeyWejNlZmYiNjYWbm5uBR69JaLXgz+H5VxaGpA/SWRqqk673rX5/c2WGyIiItIrTG6IiCqwuXPnFvkEWVl18xGVNQ4oJiKqwEaMGIH+/fsXesxEoiddiF4VkxsiogqscuXKqFy5stRhkD4pYYmW14HJDREREemGmdkrr5GnCxxzQ/SSCvYAIdEbhT9/pAtMboj+lT8zbP76RUT0+uUvv1GamZSJ8rFbiuhfCoUCVlZWqsX8TE1N1abbJ6KylZubi/v378PU1LTQ9cioHMjMBN5/P29/505AormK+K+H6AUODg4AuFoxkVTkcjmcnZ35h0V5lZMDHDz4fF8iTG4k8PIsxpyx+M0hk8ng6OgIOzu7Ilc3JqKyY2RkVGABUiJtvRHJzcqVK7Fw4UIkJCSgUaNGWLFiRbHrpmzfvh3Tpk1DXFwcatWqhfnz58PHx+c1Rkz6TqFQsM+fiKickjw93rp1K4KDgxESEoILFy6gUaNG8Pb2LrJb4PTp0xg0aBCGDRuG33//Hb169UKvXr3w559/vubIiYiI6E0k+cKZrVq1QosWLfDtt98CyBtQ5uTkhNGjR2PixIkF6g8YMABpaWnYv3+/qqx169Zo3Lgx1qxZU+L13oSFMzXBrioiIip33pCFMyXtlsrOzsb58+cxadIkVZlcLoeXlxeioqIKPScqKgrBwcFqZd7e3ti9e3dZhvracVwOERFR6Uia3CQlJSEnJwf29vZq5fb29rh69Wqh5yQkJBRaPyEhodD6WVlZyMrKUr1OTk4GkJcBloXsjNQyed8PVx8rk/fV1MohzbQ+J2jj+TJ5XyIiekO9ODtxSopOn5jK/72tSYfTGzGguCzNmzcPM2fOLFDu5OQkQTTl1/9Glq/3JSIiiVWtWiZv++TJE1haWhZbR9LkxtbWFgqFAomJiWrliYmJqvlGXubg4KBV/UmTJql1Y+Xm5uLhw4ewsbHR6TwKKSkpcHJywq1bt8pkLA+VDu/Lm4n35c3De/Jm4n15TgiBJ0+eoKoGSZOkyY2RkRGaNWuGiIgI9OrVC0Be8hEREYFRo0YVeo6npyciIiIwbtw4VdnPP/8MT0/PQusrlUoolUq1MisrK12EXygLC4sK/w/wTcT78mbifXnz8J68mXhf8pTUYpNP8m6p4OBg+Pn5oXnz5mjZsiWWLl2KtLQ0BAQEAACGDh2KatWqYd68eQCAsWPHomPHjvjmm2/g6+uLLVu24Ny5c/juu++k/BhERET0hpA8uRkwYADu37+P6dOnIyEhAY0bN8ahQ4dUg4bj4+PVZqts06YNNm3ahKlTp2Ly5MmoVasWdu/ejfr160v1EYiIiOgNInlyAwCjRo0qshsqMjKyQFm/fv3Qr1+/Mo5KO0qlEiEhIQW6wEhavC9vJt6XNw/vyZuJ96V0JJ/Ej4iIiEiXJF9+gYiIiEiXmNwQERGRXmFyQ0RERHqFyQ0RERHpFSY3OrJy5Uq4urrC2NgYrVq1wm+//SZ1SBXGvHnz0KJFC1SqVAl2dnbo1asXYmJi1OpkZmYiKCgINjY2MDc3x/vvv19gpmsqW19//TVkMpnaBJy8L9K4ffs2PvjgA9jY2MDExAQNGjTAuXPnVMeFEJg+fTocHR1hYmICLy8vXLt2TcKI9VtOTg6mTZsGNzc3mJiYoEaNGpg9e7baGkq8J1oS9Mq2bNkijIyMxPr168Vff/0lhg8fLqysrERiYqLUoVUI3t7eIjQ0VPz5558iOjpa+Pj4CGdnZ5GamqqqM2LECOHk5CQiIiLEuXPnROvWrUWbNm0kjLpi+e2334Srq6to2LChGDt2rKqc9+X1e/jwoXBxcRH+/v7izJkz4saNG+Lw4cPi+vXrqjpff/21sLS0FLt37xYXL14UPXv2FG5ubiIjI0PCyPXXnDlzhI2Njdi/f7+IjY0V27dvF+bm5mLZsmWqOrwn2mFyowMtW7YUQUFBqtc5OTmiatWqYt68eRJGVXHdu3dPABDHjx8XQgjx+PFjYWhoKLZv366qc+XKFQFAREVFSRVmhfHkyRNRq1Yt8fPPP4uOHTuqkhveF2l8+eWXol27dkUez83NFQ4ODmLhwoWqssePHwulUik2b978OkKscHx9fUVgYKBaWZ8+fcSQIUOEELwnpcFuqVeUnZ2N8+fPw8vLS1Uml8vh5eWFqKgoCSOruJKTkwEAlStXBgCcP38eT58+VbtHHh4ecHZ25j16DYKCguDr66v2/QO8L1LZu3cvmjdvjn79+sHOzg5NmjTB2rVrVcdjY2ORkJCgdl8sLS3RqlUr3pcy0qZNG0RERODvv/8GAFy8eBEnT55Et27dAPCelMYbMUNxeZaUlIScnBzVchH57O3tcfXqVYmiqrhyc3Mxbtw4tG3bVrUkR0JCAoyMjAosmGpvb4+EhAQJoqw4tmzZggsXLuDs2bMFjvG+SOPGjRtYvXo1goODMXnyZJw9exZjxoyBkZER/Pz8VN99Yf+n8b6UjYkTJyIlJQUeHh5QKBTIycnBnDlzMGTIEADgPSkFJjekV4KCgvDnn3/i5MmTUodS4d26dQtjx47Fzz//DGNjY6nDoX/l5uaiefPmmDt3LgCgSZMm+PPPP7FmzRr4+flJHF3FtG3bNmzcuBGbNm1CvXr1EB0djXHjxqFq1aq8J6XEbqlXZGtrC4VCUeAJj8TERDg4OEgUVcU0atQo7N+/H8eOHcNbb72lKndwcEB2djYeP36sVp/3qGydP38e9+7dQ9OmTWFgYAADAwMcP34cy5cvh4GBAezt7XlfJODo6Ii6deuqldWpUwfx8fEAoPru+X/a6/P5559j4sSJGDhwIBo0aIAPP/wQ48ePx7x58wDwnpQGk5tXZGRkhGbNmiEiIkJVlpubi4iICHh6ekoYWcUhhMCoUaOwa9cuHD16FG5ubmrHmzVrBkNDQ7V7FBMTg/j4eN6jMtSlSxdcunQJ0dHRqq158+YYMmSIap/35fVr27ZtgakS/v77b7i4uAAA3Nzc4ODgoHZfUlJScObMGd6XMpKeng65XP3XsUKhQG5uLgDek1KRekSzPtiyZYtQKpUiLCxMXL58WXz88cfCyspKJCQkSB1ahfDpp58KS0tLERkZKe7evava0tPTVXVGjBghnJ2dxdGjR8W5c+eEp6en8PT0lDDqiunFp6WE4H2Rwm+//SYMDAzEnDlzxLVr18TGjRuFqamp+N///qeq8/XXXwsrKyuxZ88e8ccff4j33nuPjx2XIT8/P1GtWjXVo+Dh4eHC1tZWfPHFF6o6vCfaYXKjIytWrBDOzs7CyMhItGzZUvz6669Sh1RhACh0Cw0NVdXJyMgQI0eOFNbW1sLU1FT07t1b3L17V7qgK6iXkxveF2ns27dP1K9fXyiVSuHh4SG+++47teO5ubli2rRpwt7eXiiVStGlSxcRExMjUbT6LyUlRYwdO1Y4OzsLY2NjUb16dTFlyhSRlZWlqsN7oh2ZEC9MgUhERERUznHMDREREekVJjdERESkV5jcEBERkV5hckNERER6hckNERER6RUmN0RERKRXmNwQERGRXmFyQ0TlVocOHbBp0yapwyiVQ4cOoXHjxqop9olId5jcEFGh/P39IZPJMGLEiALHgoKCIJPJ4O/v//oD+9fevXuRmJiIgQMHqspcXV0hk8mwZcuWAvXr1asHmUyGsLCw1xhl0bp27QpDQ0Ns3LhR6lCI9A6TGyIqkpOTE7Zs2YKMjAxVWWZmJjZt2gRnZ2cJIwOWL1+OgICAAgsOOjk5ITQ0VK3s119/RUJCAszMzF5niCXy9/fH8uXLpQ6DSO8wuSGiIjVt2hROTk4IDw9XlYWHh8PZ2RlNmjRRq3vo0CG0a9cOVlZWsLGxQffu3fHPP/+ojmdnZ2PUqFFwdHSEsbExXFxcMG/ePAB5K7vPmDEDzs7OUCqVqFq1KsaMGVNkXPfv38fRo0fRo0ePAseGDBmC48eP49atW6qy9evXY8iQITAwMNDoc5cUj0wmw+7du9XOsbKyUrUKxcXFQSaTITw8HJ07d4apqSkaNWqEqKgotXN69OiBc+fOqX1PRPTqmNwQUbECAwPVWkLWr1+PgICAAvXS0tIQHByMc+fOISIiAnK5HL1791aNKVm+fDn27t2Lbdu2ISYmBhs3boSrqysAYOfOnViyZAn+85//4Nq1a9i9ezcaNGhQZEwnT56Eqakp6tSpU+CYvb09vL29sWHDBgBAeno6tm7disDAQI0/s7bxFGXKlCmYMGECoqOjUbt2bQwaNAjPnj1THXd2doa9vT1++eUXrd+biIqm2Z8xRFRhffDBB5g0aRJu3rwJADh16hS2bNmCyMhItXrvv/++2uv169ejSpUquHz5MurXr4/4+HjUqlUL7dq1g0wmg4uLi6pufHw8HBwc4OXlBUNDQzg7O6Nly5ZFxnTz5k3Y29sX6JLKFxgYiM8++wxTpkzBjh07UKNGDTRu3Fjjz6xtPEWZMGECfH19AQAzZ85EvXr1cP36dXh4eKjqVK1aVfXdEpFusOWGiIpVpUoV+Pr6IiwsDKGhofD19YWtrW2BeteuXcOgQYNQvXp1WFhYqFpl4uPjAeSNL4mOjoa7uzvGjBmDn376SXVuv379kJGRgerVq2P48OHYtWuXWgvHyzIyMmBsbFzkcV9fX6SmpuLEiRNYv369Vq02pYmnKA0bNlTtOzo6AgDu3bunVsfExATp6elavzcRFY3JDRGVKDAwEGFhYdiwYUORiUKPHj3w8OFDrF27FmfOnMGZM2cA5I21AfLG78TGxmL27NnIyMhA//790bdvXwB5g4BjYmKwatUqmJiYYOTIkejQoQOePn1a6LVsbW3x6NGjIuM1MDDAhx9+iJCQEJw5cwZDhgzR6vOWFI9MJoMQQu2cwmI1NDRU7ctkMgAo8Oj3w4cPUaVKFa3iI6LiMbkhohJ17doV2dnZePr0Kby9vQscf/DgAWJiYjB16lR06dIFderUKTT5sLCwwIABA7B27Vps3boVO3fuxMOHDwHktWD06NEDy5cvR2RkJKKionDp0qVC42nSpAkSEhKKTXACAwNx/PhxvPfee7C2ttb6MxcXT5UqVXD37l1V3WvXrpWq9SUzMxP//PNPgcHZRPRqOOaGiEqkUChw5coV1f7LrK2tYWNjg++++w6Ojo6Ij4/HxIkT1eosXrwYjo6OaNKkCeRyObZv3w4HBwfVU0Y5OTlo1aoVTE1N8b///Q8mJiZq43Je1KRJE9ja2uLUqVPo3r17oXXq1KmDpKQkmJqaav15S4rn7bffxrfffgtPT0/k5OTgyy+/VGul0dSvv/4KpVIJT09Prc8loqKx5YaINGJhYQELC4tCj8nlcmzZsgXnz59H/fr1MX78eCxcuFCtTqVKlbBgwQI0b94cLVq0QFxcHA4ePAi5XA4rKyusXbsWbdu2RcOGDXHkyBHs27cPNjY2hV5PoVAgICCgxAnwbGxsYGJiovVnLSmeb775Bk5OTmjfvj0GDx6MCRMmlCqJ2rx5M4YMGVKqc4moaDLxcscxEVE5kJCQgHr16uHChQtFtvC8yZKSkuDu7o5z587Bzc1N6nCI9ApbboioXHJwcMC6detUT2OVN3FxcVi1ahUTG6IywJYbIqpwNm7ciE8++aTQYy4uLvjrr79ec0REpEtMboiownny5AkSExMLPWZoaFguu7mI6DkmN0RERKRXOOaGiIiI9AqTGyIiItIrTG6IiIhIrzC5ISIiIr3C5IaIiIj0CpMbIiIi0itMboiIiEivMLkhIiIivfL/zi2qpE3Q1I8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "88.78850237041276"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#task12.2 这团气体中能形成的最大质量的恒星的质量是多少?\n",
    "#• 假设分子云中的恒星形成遵循Salpter初始质量函数(最小质量的恒星为0.08太阳质量) • 这道题不是很好解析求解(需要估算一个合理的ΔM)，或者采用MC方法\n",
    "#• http://cluster.shao.ac.cn/~shen/Lecture/IMF.pdf\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import integrate\n",
    "\n",
    "# 定义Salpeter IMF函数\n",
    "def salpeter_IMF(M, alpha=-2.35):\n",
    "    \"\"\"\n",
    "    Salpeter initial mass function, proportional to M^-2.35.\n",
    "    \"\"\"\n",
    "    return M ** alpha\n",
    "\n",
    "# 归一化因子\n",
    "def normalization_factor(M_min, M_max, alpha=-2.35):\n",
    "    \"\"\"\n",
    "    计算Salpeter IMF的归一化因子\n",
    "    \"\"\"\n",
    "    return 1 / integrate.quad(salpeter_IMF, M_min, M_max)[0]\n",
    "\n",
    "# 蒙特卡罗模拟生成恒星质量\n",
    "def monte_carlo_sampling(num_samples, M_min, M_max, alpha=-2.35):\n",
    "    \"\"\"\n",
    "    使用蒙特卡罗方法从Salpeter IMF中采样恒星质量\n",
    "    \"\"\"\n",
    "    # 归一化常数\n",
    "    norm = normalization_factor(M_min, M_max, alpha)\n",
    "    \n",
    "    # 随机生成样本\n",
    "    samples = []\n",
    "    while len(samples) < num_samples:\n",
    "        M_random = np.random.uniform(M_min, M_max)\n",
    "        probability = salpeter_IMF(M_random, alpha) * norm\n",
    "        if np.random.rand() < probability:\n",
    "            samples.append(M_random)\n",
    "    \n",
    "    return np.array(samples)\n",
    "\n",
    "# 设置参数\n",
    "M_min = 0.08  # 最小质量 (太阳质量)\n",
    "M_max = 100   # 假设最大质量为100太阳质量\n",
    "num_samples = 10000  # 采样数量\n",
    "\n",
    "# 运行蒙特卡罗模拟\n",
    "samples = monte_carlo_sampling(num_samples, M_min, M_max)\n",
    "\n",
    "# 计算最大质量\n",
    "max_mass = samples.max()\n",
    "\n",
    "# 可视化结果\n",
    "plt.hist(samples, bins=100, density=True, alpha=0.7, label=\"IMF Distribution\")\n",
    "plt.axvline(max_mass, color='r', linestyle='--', label=f\"Max Mass = {max_mass:.2f} M_sun\")\n",
    "plt.xlabel(\"Mass (M_sun)\")\n",
    "plt.ylabel(\"Probability Density\")\n",
    "plt.title(\"Salpeter IMF - Star Formation\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "# 输出最大质量\n",
    "max_mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9493e03a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "太阳处于主序阶段的时间为：0.725 Gyr\n"
     ]
    }
   ],
   "source": [
    "#太阳中的氢大概有10%在主序阶段被燃烧，请估算太阳处于主序阶段的时间？太阳表面的温度大概是~5500K，请由此估算地球表面的温度。地球的反射率~0.3\n",
    "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}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d06f15a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#地球表面的温度。\n",
    "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}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}