{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "2dece90f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-20T11:24:30.907510Z",
     "start_time": "2024-11-20T11:24:30.769561Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import quad\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ba534a6",
   "metadata": {},
   "source": [
    "作业1\n",
    "\n",
    "单位平方度内，有多少r<20.5 mag 的星系？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "595fe226",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-20T10:54:54.025536Z",
     "start_time": "2024-11-20T10:54:54.021588Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "单位平方度内有: 22 个r小于20.5mag的星系\n"
     ]
    }
   ],
   "source": [
    "deltam = 20.5-15\n",
    "N = deltam*10**(0.6)\n",
    "print('单位平方度内有: {:.0f} 个r小于20.5mag的星系'.format(N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdc642b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "##是不是太少了点？"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cde25567",
   "metadata": {},
   "source": [
    "SDSS给出的本地星系的Schechter光度函数的参数为Mr*=-20.83，α=-1.0，请问单位体积范围内Mr=-16mag的星系数目是Mr=-22mag星系数目的多少倍？\n",
    "SDSS是星等极限样本（r<17.77）,请问在SDSS的观测样本中， Mr=-22mag的星系数目是Mr=-16mag星系数目的多少倍？\n",
    "• 根据以上光度函数，画出一个Mr=-22mag星系的红移分布概率（不考虑星系的K改正和演化效应）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "5543acca",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-20T11:18:47.530565Z",
     "start_time": "2024-11-20T11:18:47.523616Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "单位体积范围内Mr=-16mag的星系数目是Mr=-22mag星系数目的: 18.652倍\n"
     ]
    }
   ],
   "source": [
    "M_ = -20.83\n",
    "Mr1 = -16\n",
    "Mr2 = -22\n",
    "factor1 = np.exp(-10**(0.4*(M_-Mr1)))\n",
    "factor2 = np.exp(-10**(0.4*(M_-Mr2)))\n",
    "\n",
    "factor = factor1/factor2\n",
    "\n",
    "print('单位体积范围内Mr=-16mag的星系数目是Mr=-22mag星系数目的: {:.3f}倍'.format(factor))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "97be4cad",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-20T11:18:51.158305Z",
     "start_time": "2024-11-20T11:18:51.147754Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r<17.77下,单位体积范围内Mr=-16mag的星系数目是Mr=-22mag星系数目的: 0.074倍\n"
     ]
    }
   ],
   "source": [
    "M_lim = 17.77\n",
    "d_lim1 = 10**(-(Mr1-M_lim-5)/5)/1e6\n",
    "d_lim2 = 10**(-(Mr2-M_lim-5)/5)/1e6\n",
    "\n",
    "V_lim1 = 4*np.pi/3*d_lim1*d_lim1\n",
    "V_lim2 = 4*np.pi/3*d_lim2*d_lim2\n",
    "\n",
    "factor_lim = factor*V_lim1/V_lim2\n",
    "print('r<17.77下,单位体积范围内Mr=-16mag的星系数目是Mr=-22mag星系数目的: {:.3f}倍'.format(factor_lim))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43959d71",
   "metadata": {},
   "source": [
    "#4/3pid**3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "ddbd334d",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-20T11:29:08.111041Z",
     "start_time": "2024-11-20T11:29:07.875059Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_53643/905177659.py:9: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  fai = np.log(10)/2.5*factor2/Volume\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$\\\\phi/\\\\phi_{*}$')"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnIAAAHECAYAAAC9VcdgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de3hU1b3/8c9MwiQIJBEpkwS5RAvlIiYIZIhtwecQDILRqD0CIlAOB6o/oCqgghWo/Vm5iIogx1RPW/S0KvDYH5Ko1DTgpRIDJiAEEKmGW2VCIWYC1FzIrN8fnEwdCBgkF9bk/XqeedKs/d17r5Utnc+z9qw9DmOMEQAAAKzjbO4OAAAA4LshyAEAAFiKIAcAAGApghwAAIClCHIAAACWIsgBAABYiiAHAABgKYIcAACApQhyAAAAliLIAQAAWCpkgtyKFSvUrVs3RUZGyuPxaPPmzeetX7NmjXr27KnIyEj17dtXb731VtB2Y4zmzZunuLg4tW7dWqmpqdq7d29Qza9//Wtdf/31uuyyyxQTE1PneQ4cOKCRI0fqsssuU8eOHfXggw/q1KlTFzdYAAAAhUiQW7VqlWbMmKH58+ersLBQiYmJSktL05EjR+qs37Rpk8aMGaNJkyZp69atysjIUEZGhoqKigI1ixcv1rJly5SZman8/Hy1adNGaWlpqqioCNRUVVXp3//933XvvffWeZ6amhqNHDlSVVVV2rRpk1566SWtXLlS8+bNa9g/AAAAaJEcxhjT3J24WB6PRwMHDtRzzz0nSfL7/ercubOmT5+u2bNnn1U/atQonTx5UtnZ2YG2QYMGKSkpSZmZmTLGKD4+XjNnztSsWbMkST6fT263WytXrtTo0aODjrdy5Urdf//9KisrC2p/++23dfPNN+vLL7+U2+2WJGVmZurhhx/WP/7xD7lcrgb9OwAAgJYlvLk7cLGqqqpUUFCgOXPmBNqcTqdSU1OVl5dX5z55eXmaMWNGUFtaWprWrl0rSSouLpbX61Vqampge3R0tDwej/Ly8s4KcueSl5envn37BkJc7Xnuvfde7dy5U/369atzv8rKSlVWVgZ+9/v9Ki0t1RVXXCGHw1GvcwMAgOZljNHx48cVHx8vp7NxboJaH+SOHj2qmpqaoLAkSW63W59++mmd+3i93jrrvV5vYHtt27lq6uNc5/nmOeqyYMECPfbYY/U+DwAAuHQdPHhQV155ZaMc2/ogF4rmzJkTNGPo8/nUpUsXHTx4UFFRUc3YMwAAUF/l5eXq3Lmz2rVr12jnsD7IdejQQWFhYSopKQlqLykpUWxsbJ37xMbGnre+9mdJSYni4uKCapKSkurdt9jY2LNWz9ae91x9k6SIiAhFRESc1R4VFUWQAwDAMo35sSjrV626XC71799fubm5gTa/36/c3FylpKTUuU9KSkpQvSTl5OQE6hMSEhQbGxtUU15ervz8/HMe81zn2bFjR9Dq2ZycHEVFRal37971Pg4AAEBdrJ+Rk6QZM2ZowoQJGjBggJKTk7V06VKdPHlSEydOlCSNHz9enTp10oIFCyRJ9913n4YMGaKnnnpKI0eO1GuvvaaPP/5YL7zwgqTTyfn+++/X448/ru7duyshIUFz585VfHy8MjIyAuc9cOCASktLdeDAAdXU1Gjbtm2SpO9///tq27atbrzxRvXu3Vvjxo3T4sWL5fV69eijj2rq1Kl1zrgBAABcEBMili9fbrp06WJcLpdJTk42H330UWDbkCFDzIQJE4LqV69ebXr06GFcLpfp06ePefPNN4O2+/1+M3fuXON2u01ERIQZOnSo2bNnT1DNhAkTjKSzXhs3bgzU7Nu3z9x0002mdevWpkOHDmbmzJmmurr6gsbm8/mMJOPz+S5oPwAA0Hya4v07JJ4jF+rKy8sVHR0tn8/HZ+QAALBEU7x/W/8ZOQAAgJaKIAcAAGApghwAAIClCHIAAACWIsgBAABYiiAHAABgKYIcAACApQhyAAAAliLIAQAAWIogBwAAYCmCHAAAgKUIcgAAAJYiyAEAAFiKIAcAAGApghwAAIClCHIAAACWIsgBAABYiiAHAABgKYIcAACApQhyAAAAliLIAQAAWIogBwAAYCmCHAAAgKUIcgAAAJYiyAEAAFiKIAcAAGApghwAAIClCHIAAACWIsgBAABYiiAHAABgKYIcAACApQhyAAAAliLIAQAAWIogBwAAYCmCHAAAgKUIcgAAAJYiyAEAAFiKIAcAAGApghwAAIClCHIAAACWIsgBAABYiiAHAABgKYIcAACApQhyAAAAliLIAQAAWIogBwAAYCmCHAAAgKUIcgAAAJYiyAEAAFiKIAcAAGApghwAAIClCHIAAACWIsgBAABYiiAHAABgKYIcAACApQhyAAAAliLIAQAAWIogBwAAYCmCHAAAgKUIcgAAAJYiyAEAAFiKIAcAAGApghwAAIClQibIrVixQt26dVNkZKQ8Ho82b9583vo1a9aoZ8+eioyMVN++ffXWW28FbTfGaN68eYqLi1Pr1q2VmpqqvXv3BtWUlpZq7NixioqKUkxMjCZNmqQTJ04E1fz5z3/WoEGD1K5dO33ve9/THXfcoX379jXImAEAQMsWEkFu1apVmjFjhubPn6/CwkIlJiYqLS1NR44cqbN+06ZNGjNmjCZNmqStW7cqIyNDGRkZKioqCtQsXrxYy5YtU2ZmpvLz89WmTRulpaWpoqIiUDN27Fjt3LlTOTk5ys7O1vvvv68pU6YEthcXF+vWW2/Vv/3bv2nbtm3685//rKNHj+r2229vvD8GAABoOUwISE5ONlOnTg38XlNTY+Lj482CBQvqrL/zzjvNyJEjg9o8Ho/52c9+Zowxxu/3m9jYWPPkk08GtpeVlZmIiAjz6quvGmOM2bVrl5FktmzZEqh5++23jcPhMH//+9+NMcasWbPGhIeHm5qamkDNunXrjMPhMFVVVfUen8/nM5KMz+er9z4AAKB5NcX7t/UzclVVVSooKFBqamqgzel0KjU1VXl5eXXuk5eXF1QvSWlpaYH64uJieb3eoJro6Gh5PJ5ATV5enmJiYjRgwIBATWpqqpxOp/Lz8yVJ/fv3l9Pp1O9//3vV1NTI5/Ppf/7nf5SamqpWrVqdc0yVlZUqLy8PegEAAJzJ+iB39OhR1dTUyO12B7W73W55vd469/F6veetr/35bTUdO3YM2h4eHq727dsHahISEvTOO+/okUceUUREhGJiYnTo0CGtXr36vGNasGCBoqOjA6/OnTuftx4AALRM1ge5S5nX69XkyZM1YcIEbdmyRe+9955cLpd+8pOfyBhzzv3mzJkjn88XeB08eLAJew0AAGwR3twduFgdOnRQWFiYSkpKgtpLSkoUGxtb5z6xsbHnra/9WVJSori4uKCapKSkQM2ZiylOnTql0tLSwP4rVqxQdHS0Fi9eHKj5wx/+oM6dOys/P1+DBg2qs38RERGKiIj41rEDAICWzfoZOZfLpf79+ys3NzfQ5vf7lZubq5SUlDr3SUlJCaqXpJycnEB9QkKCYmNjg2rKy8uVn58fqElJSVFZWZkKCgoCNRs2bJDf75fH45Ek/fOf/5TTGfwnDgsLC/QRAADgojTaMoom9Nprr5mIiAizcuVKs2vXLjNlyhQTExNjvF6vMcaYcePGmdmzZwfqP/zwQxMeHm6WLFlidu/ebebPn29atWplduzYEahZuHChiYmJMW+88YbZvn27ufXWW01CQoL5+uuvAzXDhw83/fr1M/n5+eavf/2r6d69uxkzZkxge25urnE4HOaxxx4zn332mSkoKDBpaWmma9eu5p///Ge9x8eqVQAA7NMU798hEeSMMWb58uWmS5cuxuVymeTkZPPRRx8Ftg0ZMsRMmDAhqH716tWmR48exuVymT59+pg333wzaLvf7zdz5841brfbREREmKFDh5o9e/YE1Rw7dsyMGTPGtG3b1kRFRZmJEyea48ePB9W8+uqrpl+/fqZNmzbme9/7nrnlllvM7t27L2hsBDkAAOzTFO/fDmPO86l7XBLKy8sVHR0tn8+nqKio5u4OAACoh6Z4/7b+M3IAAAAtFUEOAADAUgQ5AAAASxHkAAAALEWQAwAAsBRBDgAAwFIEOQAAAEsR5AAAACxFkAMAALAUQQ4AAMBSBDkAAABLEeQAAAAsRZADAACwFEEOAADAUgQ5AAAASxHkAAAALEWQAwAAsBRBDgAAwFIEOQAAAEsR5AAAACxFkAMAALAUQQ4AAMBSBDkAAABLEeQAAAAsRZADAACwFEEOAADAUgQ5AAAASxHkAAAALEWQAwAAsBRBDgAAwFIEOQAAAEsR5AAAACxFkAMAALAUQQ4AAMBSBDkAAABLEeQAAAAsRZADAACwFEEOAADAUgQ5AAAASxHkAAAALEWQAwAAsBRBDgAAwFIEOQAAAEsR5AAAACxFkAMAALAUQQ4AAMBSBDkAAABLEeQAAAAsRZADAACwFEEOAADAUgQ5AAAASxHkAAAALEWQAwAAsBRBDgAAwFIEOQAAAEsR5AAAACxFkAMAALAUQQ4AAMBSBDkAAABLEeQAAAAsRZADAACwFEEOAADAUgQ5AAAAS4VMkFuxYoW6deumyMhIeTwebd68+bz1a9asUc+ePRUZGam+ffvqrbfeCtpujNG8efMUFxen1q1bKzU1VXv37g2qKS0t1dixYxUVFaWYmBhNmjRJJ06cOOs4S5YsUY8ePRQREaFOnTrp17/+dcMMGgAAtGghEeRWrVqlGTNmaP78+SosLFRiYqLS0tJ05MiROus3bdqkMWPGaNKkSdq6dasyMjKUkZGhoqKiQM3ixYu1bNkyZWZmKj8/X23atFFaWpoqKioCNWPHjtXOnTuVk5Oj7Oxsvf/++5oyZUrQue677z7993//t5YsWaJPP/1U69atU3JycuP8IQAAQIviMMaY5u7ExfJ4PBo4cKCee+45SZLf71fnzp01ffp0zZ49+6z6UaNG6eTJk8rOzg60DRo0SElJScrMzJQxRvHx8Zo5c6ZmzZolSfL5fHK73Vq5cqVGjx6t3bt3q3fv3tqyZYsGDBggSVq/fr1GjBihQ4cOKT4+Xrt379a1116roqIi/eAHP/jO4ysvL1d0dLR8Pp+ioqK+83EAAEDTaYr3b+tn5KqqqlRQUKDU1NRAm9PpVGpqqvLy8urcJy8vL6hektLS0gL1xcXF8nq9QTXR0dHyeDyBmry8PMXExARCnCSlpqbK6XQqPz9fkpSVlaWrrrpK2dnZSkhIULdu3fSf//mfKi0tPe+YKisrVV5eHvQCAAA4k/VB7ujRo6qpqZHb7Q5qd7vd8nq9de7j9XrPW1/789tqOnbsGLQ9PDxc7du3D9R88cUX2r9/v9asWaOXX35ZK1euVEFBgX7yk5+cd0wLFixQdHR04NW5c+fz1gMAgJbJ+iB3KfP7/aqsrNTLL7+sH//4x7rhhhv029/+Vhs3btSePXvOud+cOXPk8/kCr4MHDzZhrwEAgC2sD3IdOnRQWFiYSkpKgtpLSkoUGxtb5z6xsbHnra/9+W01Zy6mOHXqlEpLSwM1cXFxCg8PV48ePQI1vXr1kiQdOHDgnGOKiIhQVFRU0KtWjd8o7/NjemPb35X3+THV+K3/iCMAAPiOrA9yLpdL/fv3V25ubqDN7/crNzdXKSkpde6TkpISVC9JOTk5gfqEhATFxsYG1ZSXlys/Pz9Qk5KSorKyMhUUFARqNmzYIL/fL4/HI0n64Q9/qFOnTunzzz8P1Hz22WeSpK5du17wWHN2efWjRRs05sWPdN9r2zTmxY/0o0UbtL7o8AUfCwAA2C8kVq2uWrVKEyZM0G9+8xslJydr6dKlWr16tT799FO53W6NHz9enTp10oIFCySdfvzIkCFDtHDhQo0cOVKvvfaannjiCRUWFuqaa66RJC1atEgLFy7USy+9pISEBM2dO1fbt2/Xrl27FBkZKUm66aabVFJSoszMTFVXV2vixIkaMGCAXnnlFUmnA+XAgQPVtm1bLV26VH6/X1OnTlVUVJTeeeedeo+vdtVLl/tXyxFxWdA2x//+fP7u6zT8mriL/EsCAICG0hSrVsMb5ahNbNSoUfrHP/6hefPmyev1KikpSevXrw8sVjhw4ICczn9NPl5//fV65ZVX9Oijj+qRRx5R9+7dtXbt2kCIk6SHHnpIJ0+e1JQpU1RWVqYf/ehHWr9+fSDESdIf//hHTZs2TUOHDpXT6dQdd9yhZcuWBbY7nU5lZWVp+vTpGjx4sNq0aaObbrpJTz311Hcap9G/gtuZbY9l7dKw3rEKc55ZAQAAQlVIzMiFutpE3/n+1XKeMSP3Ta9OHqSUq69owp4BAIBz4TlyuCBHjld8exEAAAgZBLkQ0rFd5LcXAQCAkBESn5FrKc716TeHpNjoSCUntG/K7gAAgGbGjJxlzgxztb/PT+/NQgcAAFqYJglyZWVl5/0d9fP0qETFRgffPo2NjuTRIwAAtFCNfmv14MGDuu6667R792516NBBR48eVa9evbR161ZdeeWVjX36kDKsd6wykrtrc3GpjhyvUMd2p2+nMhMHAEDL1CSPH7nnnnvkcrm0bNkyTZ8+XdXV1crMzGzs04aMpli+DAAAGlbIPBB44cKF6tu3r4YPH661a9eqqKioKU4LAAAQ0ho1yP3qV78K/O+EhATddtttGjRokJ599tlA+7x58xqzCwAAACGrUYNccXFx4H+fPHlS1dXVOnHiRKDd4eCzXQAAAN9Vk3xGrqKiQtdcc41mzZqlp556SkVFRYqIiGjs04YMPiMHAIB9QuYrup544gklJyfrnnvu0cCBA/XrX/+6KU4LAAAQ0hp9Ru7gwYO65pprtH37dnXt2lXFxcVKSkrSzp07efxIPTEjBwCAfUJiRq5z587Kz89X165dJZ1e9JCfn0+IAwAAuEhNcmu1Z8+e5/0dAAAAF65BV63u2rVLzz33nDZt2qSvvvpK4eHh6tq1q/r166ehQ4cqLS1NYWFhDXlKAACAFqvBPiP39NNP66GHHlLt4Wp/fvMRI263Wz//+c81c+ZMtWrVqiFO2yLwGTkAAOzTFO/fDRLk3n77bY0cOTLwe/v27dWrVy85nU4dOnRI+/btCwp2vXr1UlZWlhISEi721C0CQQ4AAPtYs9jh6aefliS1bdtWL7/8skpKSvTBBx/ovffe0+eff66jR49q1apVGjZsmIwx2rVrl66//nrt27evIU4PAADQIl3QjNxVV12lxMREXXvttUpMTFRiYqKuvvpqXXHFFSorK9OyZcs0derU8x4jJydH48ePV0lJiRITE1VQUCCns0nWXFiLGTkAAOzTFO/fF7TYYd++fdq/f7/WrVsXaGvTpo1OnDghh8OhsrIy7dixQ7179z7nooZhw4bp3Xff1aBBg7R9+3a99NJLmjhx4sWNAgAAoAW6oBm5yZMna/v27SoqKtLXX3/9r4Oc8Z2pLpdLffr0UVJSkpKSktSvXz8lJiaqbdu2gZonn3xSDz/8sG644QZt2LChAYYSupiRAwDAPpfcjNyLL74o6fSK1L179+qTTz7RJ598ooULF8oYE1jQUFlZqcLCQm3dujWwr8Ph0FVXXRUId7WrVrdv395QYwEAAGhRGmTValxcnI4cOaJVq1YpMjJS27Zt07Zt27R169agFatS8OydMUYOh0MjRozQtddeq759+6pv377q2bMnz5v7BmbkAACwjzWPH7ntttu0bt06PfDAA1qyZEnQtuPHjwdCXe3PXbt2qbq6+qxbsrVatWqlXr16Bc3otWQEOQAA7GNNkMvKytKtt96qVq1aKTs7W8OGDTtvfXV1tV5//XXdddddcjqd8ng8Kioq0vHjx//VMYdDNTU1F9u1kECQAwDAPtY8Ry49PV0jRoxQdXW10tPT9eSTT6q6uvqc9a1atdKOHTskSR06dNCHH34on8+n4uJiZWVl6YknntDo0aMbomsAAAAhq8G+osvn82nw4MHasWOHHA6HOnXqpLvvvls33nijEhMTdfnll0uSvvjiC2VmZuqZZ56R3+/X7bffrjVr1jREF0IWM3IAANjHmlurtY4fP65x48YFnjP3zc/ARUREyBijqqoqSacXOrhcLuXn5ysxMbGhuhCSCHIAANjHmlurtdq1a6e1a9cqKytLgwcPDjySxBijiooKVVZWBn6PiYnR66+/TogDAAD4jhp0Ru5Mhw4d0oYNG1RYWKh9+/apqqpKl19+uVJSUjRu3DhFR0c31qlDCjNyAADY55K8tfrJJ58wi9bECHIAANjnkry1et1116lbt26aOnWq1q9fH/jMGwAAAJrWBc/IOZ2ns1/tQobLLrtMqampSk9P18iRI+V2uxu+ly0cM3IAANjnkry1+uWXXyo7O1tZWVnKzc1VRUXF6QM5HHI4HOrfv7/S09OVnp7OLdgGQpADAMA+l2SQ+6avv/5af/nLX5Sdna3s7GwdPnz49EH/d7buyiuv1MiRI5Wenq6hQ4fK5XI1TK9bGIIcAAD2ueSD3JkKCgqUlZWl7OxsFRYWnj4Bt2AvGkEOAAD7WBfkvumbt2A3bNigr7/++vQJuQV7wQhyAADYx+og900VFRX6y1/+EpitO/MWbKdOnZSenq577rlHffv2bezuWIcgBwCAfZri/Tu8UY56hsjISN188826+eabJZ2+BVs7W1dYWKhDhw4pMzNTbrebIAcAAFBPTRLkztS/f3/1799f8+fP15dffqmsrCy9+eabuuyyy5qjOwAAAFZqkluruDjcWgUAwD7W3VrdtWuXnnvuOW3atElfffWVwsPD1bVrV/Xr109Dhw5VWlqawsLCGvKUAAAALVaDzcg9/fTTeuihh1R7uNqftQsaJMntduvnP/+5Zs6cqVatWjXEaVsEZuQAALCPNatW3377bY0cOTLwe/v27dWrVy85nU4dOnRI+/btCwp2vXr1UlZWlhISEi721C0CQQ4AAPs0xfu3syEO8vTTT0uS2rZtq5dfflklJSX64IMP9N577+nzzz/X0aNHtWrVKg0bNkzGGO3atUvXX3+99u3b1xCnBwAAaJEuaEbuqquuUmJioq699lolJiYqMTFRV199ta644gqVlZVp2bJlmjp16nmPkZOTo/Hjx6ukpESJiYkqKCiQ09kgeTJkMSMHAIB9LrnFDvv27dP+/fu1bt26QFubNm104sQJORwOlZWVaceOHerdu/c5FzUMGzZM7777rgYNGqTt27frpZde0sSJEy9uFAAAAC3QBc3ITZ48Wdu3b1dRUVHgK7ek4AUNkuRyudSnTx8lJSUpKSlJ/fr1U2Jiotq2bRuoefLJJ/Xwww/rhhtu0IYNGxpgKKGLGTkAAOxzyc3Ivfjii5JOr0jdu3evPvnkE33yySdauHChjDGBBQ2VlZUqLCzU1q1bA/s6HA5dddVVgXBXu2p1+/btDTUWAACAFqVBVq3GxcXpyJEjWrVqlSIjI7Vt2zZt27ZNW7duDVqxKgXP3hlj5HA4NGLECF177bXq27ev+vbtq549e/K8uW9gRg4AAPtY8/iR2267TevWrdMDDzygJUuWBG07fvx4INTV/ty1a5eqq6vPuiVbq1WrVurVq1fQjF5LRpADAMA+1gS5rKws3XrrrWrVqpWys7M1bNiw89ZXV1fr9ddf11133SWn0ymPx6OioiIdP378Xx1zOFRTU3OxXQsJBDkAAOxjzXPk0tPTNWLECFVXVys9PV1PPvmkqqurz1nfqlUr7dixQ5LUoUMHffjhh/L5fCouLlZWVpaeeOIJjR49uiG6BgAAELIa7Cu6fD6fBg8erB07dsjhcKhTp066++67deONNyoxMVGXX365JOmLL75QZmamnnnmGfn9ft1+++1as2ZNQ3QhZDEjBwCAfay5tVrr+PHjGjduXOA5c9/8DFxERISMMaqqqpJ0eqGDy+VSfn6+EhMTG6oLIYkgBwCAfay5tVqrXbt2Wrt2rbKysjR48ODAI0mMMaqoqFBlZWXg95iYGL3++uuEOAAAgO+oQWfkznTo0CFt2LBBhYWF2rdvn6qqqnT55ZcrJSVF48aNU3R0dGOdOqQwIwcAgH2su7WKxkGQAwDAPtbdWgUAAEDTIcgBAABYiiAHAABgKYIcAACApUIqyK1YsULdunVTZGSkPB6PNm/efN76NWvWqGfPnoqMjFTfvn311ltvBW03xmjevHmKi4tT69atlZqaqr179wbVlJaWauzYsYqKilJMTIwmTZqkEydO1Hm+v/3tb2rXrp1iYmIubqAAAAAKoSC3atUqzZgxQ/Pnz1dhYaESExOVlpamI0eO1Fm/adMmjRkzRpMmTdLWrVuVkZGhjIwMFRUVBWoWL16sZcuWKTMzU/n5+WrTpo3S0tJUUVERqBk7dqx27typnJwcZWdn6/3339eUKVPOOl91dbXGjBmjH//4xw0/eAAA0CKFzONHPB6PBg4cqOeee06S5Pf71blzZ02fPl2zZ88+q37UqFE6efKksrOzA22DBg1SUlKSMjMzZYxRfHy8Zs6cqVmzZkk6/TVkbrdbK1eu1OjRo7V792717t1bW7Zs0YABAyRJ69ev14gRI3To0CHFx8cHjv3www/ryy+/1NChQ3X//ferrKys3mPj8SMAANiHx4/UU1VVlQoKCpSamhpoczqdSk1NVV5eXp375OXlBdVLUlpaWqC+uLhYXq83qCY6OloejydQk5eXp5iYmECIk6TU1FQ5nU7l5+cH2jZs2KA1a9ZoxYoV9RpPZWWlysvLg14AAABnCokgd/ToUdXU1Mjtdge1u91ueb3eOvfxer3nra/9+W01HTt2DNoeHh6u9u3bB2qOHTumn/70p1q5cmW90/iCBQsUHR0deHXu3Lle+wEAgJYlJILcpWzy5Mm66667NHjw4HrvM2fOHPl8vsDr4MGDjdhDAABgq5AIch06dFBYWJhKSkqC2ktKShQbG1vnPrGxseetr/35bTVnLqY4deqUSktLAzUbNmzQkiVLFB4ervDwcE2aNEk+n0/h4eH63e9+V2ffIiIiFBUVFfQCAAA4U0gEOZfLpf79+ys3NzfQ5utSBtUAABrZSURBVPf7lZubq5SUlDr3SUlJCaqXpJycnEB9QkKCYmNjg2rKy8uVn58fqElJSVFZWZkKCgoCNRs2bJDf75fH45F0+nN027ZtC7x+9atfqV27dtq2bZtuu+22hvkDAACAFim8uTvQUGbMmKEJEyZowIABSk5O1tKlS3Xy5ElNnDhRkjR+/Hh16tRJCxYskCTdd999GjJkiJ566imNHDlSr732mj7++GO98MILkiSHw6H7779fjz/+uLp3766EhATNnTtX8fHxysjIkCT16tVLw4cP1+TJk5WZmanq6mpNmzZNo0ePDqxY7dWrV1A/P/74YzmdTl1zzTVN9acBAAAhKmSC3KhRo/SPf/xD8+bNk9frVVJSktavXx9YrHDgwAE5nf+agLz++uv1yiuv6NFHH9Ujjzyi7t27a+3atUEB66GHHtLJkyc1ZcoUlZWV6Uc/+pHWr1+vyMjIQM0f//hHTZs2TUOHDpXT6dQdd9yhZcuWNd3AAQBAixUyz5ELZTxHDgAA+/AcOQAAAJwTQQ4AAMBSBDkAAABLEeQAAAAsRZADAACwFEEOAADAUgQ5AAAASxHkAAAALEWQAwAAsBRBDgAAwFIEOQAAAEsR5AAAACxFkAMAALAUQQ4AAMBSBDkAAABLEeQAAAAsRZADAACwFEEOAADAUgQ5AAAASxHkAAAALEWQAwAAsBRBDgAAwFIEOQAAAEsR5AAAACxFkAMAALAUQQ4AAMBSBDkAAABLEeQAAAAsRZADAACwFEEOAADAUgQ5AAAASxHkAAAALEWQAwAAsBRBDgAAwFIEOQAAAEsR5AAAACxFkAMAALAUQQ4AAMBSBDkAAABLEeQAAAAsRZADAACwFEEOAADAUgQ5AAAASxHkAAAALEWQAwAAsBRBDgAAwFIEOQAAAEsR5AAAACxFkAMAALAUQQ4AAMBSBDkAAABLEeQAAAAsRZADAACwFEEOAADAUgQ5AAAASxHkAAAALEWQAwAAsBRBDgAAwFIEOQAAAEsR5AAAACxFkAMAALAUQQ4AAMBSBDkAAABLhVSQW7Fihbp166bIyEh5PB5t3rz5vPVr1qxRz549FRkZqb59++qtt94K2m6M0bx58xQXF6fWrVsrNTVVe/fuDaopLS3V2LFjFRUVpZiYGE2aNEknTpwIbH/33Xd16623Ki4uTm3atFFSUpL++Mc/NtygAQBAixUyQW7VqlWaMWOG5s+fr8LCQiUmJiotLU1Hjhyps37Tpk0aM2aMJk2apK1btyojI0MZGRkqKioK1CxevFjLli1TZmam8vPz1aZNG6WlpamioiJQM3bsWO3cuVM5OTnKzs7W+++/rylTpgSd59prr9Xrr7+u7du3a+LEiRo/fryys7Mb748BAABaBIcxxjR3JxqCx+PRwIED9dxzz0mS/H6/OnfurOnTp2v27Nln1Y8aNUonT54MClSDBg1SUlKSMjMzZYxRfHy8Zs6cqVmzZkmSfD6f3G63Vq5cqdGjR2v37t3q3bu3tmzZogEDBkiS1q9frxEjRujQoUOKj4+vs68jR46U2+3W7373u3qNrby8XNHR0fL5fIqKirqgvwsAAGgeTfH+HRIzclVVVSooKFBqamqgzel0KjU1VXl5eXXuk5eXF1QvSWlpaYH64uJieb3eoJro6Gh5PJ5ATV5enmJiYgIhTpJSU1PldDqVn59/zv76fD61b9/+nNsrKytVXl4e9AIAADhTSAS5o0ePqqamRm63O6jd7XbL6/XWuY/X6z1vfe3Pb6vp2LFj0Pbw8HC1b9/+nOddvXq1tmzZookTJ55zPAsWLFB0dHTg1blz53PWAgCAliskgpwtNm7cqIkTJ+rFF19Unz59zlk3Z84c+Xy+wOvgwYNN2EsAAGCLkAhyHTp0UFhYmEpKSoLaS0pKFBsbW+c+sbGx562v/fltNWcupjh16pRKS0vPOu97772n9PR0PfPMMxo/fvx5xxMREaGoqKigFwAAwJlCIsi5XC71799fubm5gTa/36/c3FylpKTUuU9KSkpQvSTl5OQE6hMSEhQbGxtUU15ervz8/EBNSkqKysrKVFBQEKjZsGGD/H6/PB5PoO3dd9/VyJEjtWjRoqAVrQAAABcjvLk70FBmzJihCRMmaMCAAUpOTtbSpUt18uTJwGfRxo8fr06dOmnBggWSpPvuu09DhgzRU089pZEjR+q1117Txx9/rBdeeEGS5HA4dP/99+vxxx9X9+7dlZCQoLlz5yo+Pl4ZGRmSpF69emn48OGaPHmyMjMzVV1drWnTpmn06NGBFasbN27UzTffrPvuu0933HFH4LNzLpfrvAseAAAAvpUJIcuXLzddunQxLpfLJCcnm48++iiwbciQIWbChAlB9atXrzY9evQwLpfL9OnTx7z55ptB2/1+v5k7d65xu90mIiLCDB061OzZsyeo5tixY2bMmDGmbdu2JioqykycONEcP348sH3ChAlG0lmvIUOG1HtcPp/PSDI+n6/+fwwAANCsmuL9O2SeIxfKeI4cAAD24TlyAAAAOCeCHAAAgKUIcgAAAJYiyAEAAFiKIAcAAGApghwAAIClCHIAAACWIsgBAABYKmS+oqslq/EbbS4u1ZHjFerYLlLJCe0V5nQ0d7cAAEAjI8hZbn3RYT2WtUuHfRWBtrjoSM1P763h18Q1Y88AAEBj49aqxdYXHda9fygMCnGS5PVV6N4/FGp90eFm6hkAAGgKBDlL1fiNHsvapbq+KLe27bGsXarx81W6AACEKoKcpTYXl541E/dNRtJhX4U2F5c2XacAAECTIshZ6sjxc4e471IHAADsQ5CzVMd2kQ1aBwAA7EOQs1RyQnvFRUfqXA8Zcej06tXkhPZN2S0AANCECHKWCnM6ND+9tySdFeZqf5+f3pvnyQEAEMIIchYbfk2cnr/7OsVGB98+jY2O1PN3X8dz5AAACHE8ENhyw6+J07DesXyzAwAALRBBLgSEOR1KufqK5u4GAABoYtxaBQAAsBRBDgAAwFIEOQAAAEsR5AAAACxFkAMAALAUQQ4AAMBSBDkAAABLEeQAAAAsRZADAACwFEEOAADAUgQ5AAAASxHkAAAALEWQAwAAsBRBDgAAwFIEOQAAAEsR5AAAACxFkAMAALAUQQ4AAMBS4c3dATSOGr/R5uJSHTleoY7tIpWc0F5hTkdzdwsAADQgglwIWl90WI9l7dJhX0WgLS46UvPTe2v4NXHN2DMAANCQuLUaYtYXHda9fygMCnGS5PVV6N4/FGp90eFm6hkAAGhoBLkQUuM3eixrl0wd22rbHsvapRp/XRUAAMA2BLkQsrm49KyZuG8ykg77KrS5uLTpOgUAABoNQS6EHDl+7hD3XeoAAMCljSAXQjq2i2zQOgAAcGkjyIWQ5IT2iouO1LkeMuLQ6dWryQntm7JbAACgkRDkQkiY06H56b0l6awwV/v7/PTePE8OAIAQQZALMcOvidPzd1+n2Ojg26ex0ZF6/u7reI4cAAAhhAcCh6Dh18RpWO/Ys77ZQZLyPj/Gtz0AABAiCHIhKszpUMrVVwR+59seAAAIPdxabQH4tgcAAEITQS7E8W0PAACELoJciOPbHgAACF18Ri7E1fdbHN7+39urLIAAAMAeBLkQV99vcXg5b79eztvPAggAACzCrdUQ923f9nAmFkAAAGAPglyIO9+3PdTF/O/rkf+3Q/9v69+V9/kxFkIAAHCJchhjeJe+xJWXlys6Olo+n09RUVHf6Rh1PUeuvuKiIzV3ZC9d3iaChwkDAFBPDfH+/W0IchZoqP8QavxGm4tL9XbRYb2ct/+i+kS4AwDg/JoiyLHYoQX55rc9XGyQO+yr0P95ZWtQG+EOAICmRZBrgWoXQHh9FXU+KPi7upBwJ6nO74KtTxvBEACA07i1aoHGmJqt/douSQ0a5uoj5rJWkqSyf1ZfcFtjBMNLua1/18tVsP+rJu8HYRkALh6fkYOkxvsP4WIWQFxKLiYYXuptTof0zUXDTXHOlhaWL6W25grutHG9aWucto3b92tYvwSCXH2tWLFCTz75pLxerxITE7V8+XIlJyefs37NmjWaO3eu9u3bp+7du2vRokUaMWJEYLsxRvPnz9eLL76osrIy/fCHP9Tzzz+v7t27B2pKS0s1ffp0ZWVlyel06o477tCzzz6rtm3bBmq2b9+uqVOnasuWLfre976n6dOn66GHHqr3uBoz0dcugPD6vtb/fXO3vjpZ1eQzdLj0XUrhNpTbmiO408b1pq3x2kq/8ung0jsJcvWxatUqjR8/XpmZmfJ4PFq6dKnWrFmjPXv2qGPHjmfVb9q0SYMHD9aCBQt0880365VXXtGiRYtUWFioa665RpK0aNEiLViwQC+99JISEhI0d+5c7dixQ7t27VJk5OlvTLjpppt0+PBh/eY3v1F1dbUmTpyogQMH6pVXXpF0OoT16NFDqampmjNnjnbs2KH/+I//0NKlSzVlypR6ja0ppmal5r3dCgBAqPFX/pMgV18ej0cDBw7Uc889J0ny+/3q3Lmzpk+frtmzZ59VP2rUKJ08eVLZ2dmBtkGDBikpKUmZmZkyxig+Pl4zZ87UrFmzJEk+n09ut1srV67U6NGjtXv3bvXu3VtbtmzRgAEDJEnr16/XiBEjdOjQIcXHx+v555/XL37xC3m9XrlcLknS7NmztXbtWn366af1GltTBTkpdG63AgDQ3JoiyIXEqtWqqioVFBRozpw5gTan06nU1FTl5eXVuU9eXp5mzJgR1JaWlqa1a9dKkoqLi+X1epWamhrYHh0dLY/Ho7y8PI0ePVp5eXmKiYkJhDhJSk1NldPpVH5+vm677Tbl5eVp8ODBgRBXe55Fixbpq6++0uWXX35W3yorK1VZWRn43efzSTod6Brb9V3a6K17B6hg31f6x4kKfa9tpMr+WaVFf/5UJeWV334AAAAg6XSQk05/VKuxhESQO3r0qGpqauR2u4Pa3W73OWe9vF5vnfVerzewvbbtfDVn3rYNDw9X+/btg2oSEhLOOkbttrqC3IIFC/TYY4+d1d65c+c6xwIAAC5dx44dU3R0dKMcOySCXKiZM2dO0GxhWVmZunbtqgMHDjTafwion/LycnXu3FkHDx5s9NvcOD+uxaWDa3Hp4FpcWnw+n7p06aL27ds32jlCIsh16NBBYWFhKikpCWovKSlRbGxsnfvExsaet772Z0lJieLi4oJqkpKSAjVHjhwJOsapU6dUWloadJy6zvPNc5wpIiJCERERZ7VHR0fzD/MSERUVxbW4RHAtLh1ci0sH1+LS4nQ6G+/YjXbkJuRyudS/f3/l5uYG2vx+v3Jzc5WSklLnPikpKUH1kpSTkxOoT0hIUGxsbFBNeXm58vPzAzUpKSkqKytTQUFBoGbDhg3y+/3yeDyBmvfff1/V1dVB5/nBD35Q521VAACA+gqJICdJM2bM0IsvvqiXXnpJu3fv1r333quTJ09q4sSJkqTx48cHLYa47777tH79ej311FP69NNP9ctf/lIff/yxpk2bJklyOBy6//779fjjj2vdunXasWOHxo8fr/j4eGVkZEiSevXqpeHDh2vy5MnavHmzPvzwQ02bNk2jR49WfHy8JOmuu+6Sy+XSpEmTtHPnTq1atUrPPvvsWQstAAAALpgJIcuXLzddunQxLpfLJCcnm48++iiwbciQIWbChAlB9atXrzY9evQwLpfL9OnTx7z55ptB2/1+v5k7d65xu90mIiLCDB061OzZsyeo5tixY2bMmDGmbdu2JioqykycONEcP348qOaTTz4xP/rRj0xERITp1KmTWbhw4QWNq6KiwsyfP99UVFRc0H5oeFyLSwfX4tLBtbh0cC0uLU1xPULmOXIAAAAtTcjcWgUAAGhpCHIAAACWIsgBAABYiiAHAABgKYJcM1ixYoW6deumyMhIeTwebd68+bz1a9asUc+ePRUZGam+ffvqrbfeCtpujNG8efMUFxen1q1bKzU1VXv37m3MIYSMhrwW1dXVevjhh9W3b1+1adNG8fHxGj9+vL788svGHkZIaOh/F990zz33yOFwaOnSpQ3d7ZDVGNdj9+7duuWWWxQdHa02bdpo4MCBOnDgQGMNIWQ09LU4ceKEpk2bpiuvvFKtW7dW7969lZmZ2ZhDCBkXci127typO+64Q926dTvv//9c6PU9S6Oth0WdXnvtNeNyuczvfvc7s3PnTjN58mQTExNjSkpK6qz/8MMPTVhYmFm8eLHZtWuXefTRR02rVq3Mjh07AjULFy400dHRZu3ateaTTz4xt9xyi0lISDBff/11Uw3LSg19LcrKykxqaqpZtWqV+fTTT01eXp5JTk42/fv3b8phWakx/l3U+tOf/mQSExNNfHy8eeaZZxp7KCGhMa7H3/72N9O+fXvz4IMPmsLCQvO3v/3NvPHGG+c8Jk5rjGsxefJkc/XVV5uNGzea4uJi85vf/MaEhYWZN954o6mGZaULvRabN282s2bNMq+++qqJjY2t8/9/LvSYdSHINbHk5GQzderUwO81NTUmPj7eLFiwoM76O++804wcOTKozePxmJ/97GfGmNPPuouNjTVPPvlkYHtZWZmJiIgwr776aiOMIHQ09LWoy+bNm40ks3///obpdIhqrGtx6NAh06lTJ1NUVGS6du1KkKunxrgeo0aNMnfffXfjdDiENca16NOnj/nVr34VVHPdddeZX/ziFw3Y89Bzodfim871/z8Xc8xa3FptQlVVVSooKFBqamqgzel0KjU1VXl5eXXuk5eXF1QvSWlpaYH64uJieb3eoJro6Gh5PJ5zHhONcy3q4vP55HA4FBMT0zAdD0GNdS38fr/GjRunBx98UH369Gmczoegxrgefr9fb775pnr06KG0tDR17NhRHo9Ha9eubbyBhIDG+rdx/fXXa926dfr73/8uY4w2btyozz77TDfeeGPjDCQEfJdr0VTHJMg1oaNHj6qmpkZutzuo3e12y+v11rmP1+s9b33tzws5JhrnWpypoqJCDz/8sMaMGcOXV59HY12LRYsWKTw8XD//+c8bvtMhrDGux5EjR3TixAktXLhQw4cP1zvvvKPbbrtNt99+u957773GGUgIaKx/G8uXL1fv3r115ZVXyuVyafjw4VqxYoUGDx7c8IMIEd/lWjTVMcO/09kBnFd1dbXuvPNOGWP0/PPPN3d3WpyCggI9++yzKiwslMPhaO7utHh+v1+SdOutt+qBBx6QJCUlJWnTpk3KzMzUkCFDmrN7Lc7y5cv10Ucfad26deratavef/99TZ06VfHx8WfN5uHSx4xcE+rQoYPCwsJUUlIS1F5SUqLY2Ng694mNjT1vfe3PCzkmGuda1KoNcfv371dOTg6zcd+iMa7FBx98oCNHjqhLly4KDw9XeHi49u/fr5kzZ6pbt26NMo5Q0RjXo0OHDgoPD1fv3r2Danr16sWq1fNojGvx9ddf65FHHtHTTz+t9PR0XXvttZo2bZpGjRqlJUuWNM5AQsB3uRZNdUyCXBNyuVzq37+/cnNzA21+v1+5ublKSUmpc5+UlJSgeknKyckJ1CckJCg2Njaopry8XPn5+ec8JhrnWkj/CnF79+7VX/7yF11xxRWNM4AQ0hjXYty4cdq+fbu2bdsWeMXHx+vBBx/Un//858YbTAhojOvhcrk0cOBA7dmzJ6jms88+U9euXRt4BKGjMa5FdXW1qqur5XQGv/2HhYUFZk5xtu9yLZrsmPVeFoEG8dprr5mIiAizcuVKs2vXLjNlyhQTExNjvF6vMcaYcePGmdmzZwfqP/zwQxMeHm6WLFlidu/ebebPn1/n40diYmLMG2+8YbZv325uvfVWHj9SDw19Laqqqswtt9xirrzySrNt2zZz+PDhwKuysrJZxmiLxvh3cSZWrdZfY1yPP/3pT6ZVq1bmhRdeMHv37jXLly83YWFh5oMPPmjy8dmkMa7FkCFDTJ8+fczGjRvNF198YX7/+9+byMhI81//9V9NPj6bXOi1qKysNFu3bjVbt241cXFxZtasWWbr1q1m79699T5mfRDkmsHy5ctNly5djMvlMsnJyeajjz4KbBsyZIiZMGFCUP3q1atNjx49jMvlMn369DFvvvlm0Ha/32/mzp1r3G63iYiIMEOHDjV79uxpiqFYryGvRXFxsZFU52vjxo1NNCJ7NfS/izMR5C5MY1yP3/72t+b73/++iYyMNImJiWbt2rWNPYyQ0NDX4vDhw+anP/2piY+PN5GRkeYHP/iBeeqpp4zf72+K4VjtQq7Fud4ThgwZUu9j1ofDGGO+05wgAAAAmhWfkQMAALAUQQ4AAMBSBDkAAABLEeQAAAAsRZADAACwFEEOAADAUgQ5AAAASxHkAAAALEWQAwAAsBRBDgAAwFIEOQAAAEsR5AAAACxFkAMAALAUQQ4AmsDKlSvlcDjq/frlL3/Z3F0GYAGCHAAAgKXCm7sDANASZGRkaMCAAeetefDBB7V+/XpJUteuXZuiWwAsR5ADgCYQExOjmJiYc25fsWJFIMSNHTtWEydObKquAbCYwxhjmrsTANCS5ebmavjw4Tp16pSSk5P13nvvKTIysrm7BcACBDkAaEZ79+6Vx+PRV199pU6dOmnLli2Ki4tr7m4BsASLHQCgmZSVlSk9PV1fffWVWrdurTfeeIMQB+CCEOQAoBnU1NRo1KhR2rNnjxwOh1auXKn+/fs3d7cAWIYgBwDN4IEHHtA777wjSZo7d67uvPPOZu4RABvxGTkAaGIvvPCCfvazn0mSfvKTn2j16tVyOBzN3CsANiLIAUATevfdd3XjjTequrpa/fr101//+ldddtllzd0tAJYiyAFAE/n888/l8Xh07Ngxud1ubdmyRZ07d27ubgGwGJ+RA4Am4PP5lJ6ermPHjikiIkJr164lxAG4aHyzAwA0genTp2v37t2SpPvvv19t27ZVUVHROes7duyojh07NlX3AFiKW6sA0ARuuOEGvffee/Wunz9/vn75y182XocAhARurQIAAFiKGTkAAABLMSMHAABgKYIcAACApQhyAAAAliLIAQAAWIogBwAAYCmCHAAAgKUIcgAAAJYiyAEAAFiKIAcAAGApghwAAIClCHIAAACWIsgBAABYiiAHAABgKYIcAACApQhyAAAAlvr/wnxNoE2DAN4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "z = np.arange(0,10,0.001)\n",
    "C = 3*1e5 #km/s\n",
    "H = 67    #Mpc\n",
    "V = z*C   #km/s\n",
    "Distance = V/H   #Mpc\n",
    "\n",
    "Volume = 4*np.pi/3*Distance*Distance\n",
    "\n",
    "fai = np.log(10)/2.5*factor2/Volume\n",
    "\n",
    "plt.scatter(z,fai)\n",
    "plt.xlim([0,0.1])\n",
    "plt.ylim([0,0.001])\n",
    "plt.xlabel('z',size=20)\n",
    "plt.ylabel(r'$\\phi/\\phi_{*}$',size=20)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fccec7f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e56fef3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "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.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}