{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 第二章 作业1\n",
    "\n",
    "1.单位平方度内，有多少r<20.5 mag 的星系？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "单位平方度内 r < 20.5 的星系数量为: 3596.56\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import quad\n",
    "from astropy.cosmology import FlatLambdaCDM, z_at_value\n",
    "from astropy.constants import c\n",
    "from astropy import units as u\n",
    "\n",
    "cosmo = FlatLambdaCDM(H0=70, Om0=0.3)\n",
    "\n",
    "phi_star = 1.46e-2  # h^3 Mpc^-3\n",
    "M_star = -20.83\n",
    "alpha = -1.20\n",
    "r_lim = 20.5\n",
    "\n",
    "def schechter(M, M_star, phi_star, alpha):\n",
    "    factor = 0.4 * np.log(10) * phi_star\n",
    "    term1 = 10**(0.4 * (M_star - M) * (alpha + 1))\n",
    "    term2 = np.exp(-10**(0.4 * (M_star - M)))\n",
    "    return factor * term1 * term2\n",
    "\n",
    "def absolute_magnitude_limit(z, r_lim):\n",
    "    Dl = cosmo.luminosity_distance(z).to(u.pc).value\n",
    "    return r_lim - 5 * np.log10(Dl / 10)\n",
    "\n",
    "# The differential comoving volume element dV_c/dz/dSolidAngle.\n",
    "def dV_dOmega_dz(z):\n",
    "    return cosmo.differential_comoving_volume(z).value\n",
    "\n",
    "def number_density_per_sq_degree():\n",
    "    z_max = 1\n",
    "    def integrand(z):\n",
    "        M_lim = absolute_magnitude_limit(z, r_lim)\n",
    "        M_min = -30  # 低光度极限\n",
    "        # 对光度函数积分\n",
    "        phi_integral, _ = quad(schechter, M_min, M_lim, args=(M_star, phi_star, alpha))\n",
    "        # 乘以体积元素\n",
    "        return phi_integral * dV_dOmega_dz(z)\n",
    "    result, _ = quad(integrand, 0, z_max)\n",
    "    return result * (np.pi / 180)**2\n",
    "\n",
    "N_per_sq_degree = number_density_per_sq_degree()\n",
    "print(f\"单位平方度内 r < 20.5 的星系数量为: {N_per_sq_degree:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2.SDSS给出的本地星系的Schechter光度函数的参数为Mr*=-20.83，α=-1.0，请问单位体积范围内Mr=-16mag的星系数目是Mr=-22mag星系数目的多少倍？SDSS是星等极限样本（r<17.77）,请问在SDSS的观测样本中， Mr=-22 mag的星系数目是Mr=-16mag星系数目的多少倍？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "单位体积内 M_r = -22 星系数目比是 M_r = -16 星系数目比的 0.0536 倍\n"
     ]
    }
   ],
   "source": [
    "cosmo = FlatLambdaCDM(H0=70, Om0=0.3)\n",
    "phi_star = 1.46e-2\n",
    "M_star = -20.83\n",
    "alpha = -1.0\n",
    "M_r_1 = -22\n",
    "M_r_2 = -16\n",
    "\n",
    "# 计算单位体积中的星系数量比\\phi(M1) / \\phi(M2)\n",
    "def volume_density_ratio(M1, M2):\n",
    "    phi1 = schechter(M1, M_star, phi_star, alpha)\n",
    "    phi2 = schechter(M2, M_star, phi_star, alpha)\n",
    "    return phi1 / phi2\n",
    "\n",
    "R_volume = volume_density_ratio(M_r_1, M_r_2)\n",
    "print(f\"单位体积内 M_r = -22 星系数目比是 M_r = -16 星系数目比的 {R_volume:.4f} 倍\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SDSS 观测样本中 M_r = -22 mag 的星系数量是 M_r = -16 mag 星系数量的 133.3253 倍\n"
     ]
    }
   ],
   "source": [
    "def max_redshift(M_r, r_limit):\n",
    "    \"\"\"\n",
    "    r = M_r + 5 \\log_{10}(D_L / 10 pc)\n",
    "    D_L = 10^{(r - M_r)/5 + 1}\n",
    "    \"\"\"\n",
    "    D_L = (10 ** ((r_limit - M_r) / 5 + 1)) * u.pc  # 光度距离 (pc)\n",
    "    z_max = z_at_value(cosmo.luminosity_distance, D_L.to(u.Mpc), zmax=5)  # 转换为 Mpc\n",
    "    return z_max\n",
    "\n",
    "def observed_density_ratio(M_r_1, M_r_2, r_limit, M_star, phi_star, alpha):\n",
    "    z_max_1 = max_redshift(M_r_1, r_limit)\n",
    "    z_max_2 = max_redshift(M_r_2, r_limit)\n",
    "    def integrand(z, M_r):\n",
    "        phi = schechter(M_r, M_star, phi_star, alpha)\n",
    "        volume_element = cosmo.differential_comoving_volume(z).value  # 单位: Mpc^3/sr/dz\n",
    "        return phi * volume_element\n",
    "    \n",
    "    N1, _ = quad(integrand, 0, z_max_1, args=(M_r_1,))\n",
    "    N2, _ = quad(integrand, 0, z_max_2, args=(M_r_2,))\n",
    "    return N1 / N2\n",
    "\n",
    "phi_star = 1.46e-2\n",
    "M_star = -20.83\n",
    "alpha = -1.0\n",
    "r_limit = 17.77\n",
    "M_r_1 = -22\n",
    "M_r_2 = -16\n",
    "R_observed = observed_density_ratio(M_r_1=M_r_1, M_r_2=M_r_2, r_limit=r_limit, \n",
    "                                   M_star=M_star, phi_star=phi_star, alpha=alpha)\n",
    "\n",
    "print(f\"SDSS 观测样本中 M_r = -22 mag 的星系数量是 M_r = -16 mag 星系数量的 {R_observed:.4f} 倍\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# good~"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. 根据以上光度函数，画出一个Mr=-22 mag星系的红移分布概率（不考虑星系的K改正和演化效应）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1cAAAIOCAYAAABUNPd7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB7EklEQVR4nO3dd3iUVfrG8XtmUgkphPQCCb2FhBrpxSh2sSK6itjW3wIrorjgqqxljX1RQbGXXVHWggUVgdAhtASQ3iEhIRVIJW1mfn8EolmKCSS8k+T7ua65Ft458+YZd5S5Oec8x2S32+0CAAAAAFwQs9EFAAAAAEBjQLgCAAAAgDpAuAIAAACAOkC4AgAAAIA6QLgCAAAAgDpAuAIAAACAOkC4AgAAAIA6QLgCAAAAgDrgZHQBjshmsyk9PV2enp4ymUxGlwMAAADAIHa7XQUFBQoJCZHZfO65KcLVGaSnpys8PNzoMgAAAAA4iNTUVIWFhZ1zDOHqDDw9PSVV/gP08vIyuBoAAAAARsnPz1d4eHhVRjgXwtUZnFoK6OXlRbgCAAAAUKPtQjS0AAAAAIA6QLgCAAAAgDpAuAIAAACAOsCeKwAAAKAW7Ha7KioqZLVajS4FdcBiscjJyalOjmAiXAEAAAA1VFZWpiNHjqi4uNjoUlCHmjVrpuDgYLm4uFzQfQhXAAAAQA3YbDYdOHBAFotFISEhcnFxqZPZDhjHbrerrKxM2dnZOnDggNq3b/+HBwWfC+EKAAAAqIGysjLZbDaFh4erWbNmRpeDOuLu7i5nZ2cdOnRIZWVlcnNzO+970dACAAAAqIULmdmAY6qr/0/5ZAAAAABAHSBcAQAAAEAdIFwBAAAAQB0gXAEAAABNxJAhQ2QymfT8889Xu2632xUbGyuTyaRnnnnGoOpqJj4+Xn369JGnp6cCAgI0cuRI7dq1q9Zj6gPhCgAAAGgC7Ha7Nm7cqNatW2vLli3Vnvvkk0+Unp4uSerZs6cR5dXYsmXLNG7cOK1Zs0YLFy5UeXm5Lr/8chUVFdVqTH2gFTsAAADQBOzZs0cFBQV6+OGH9fXXX1ddLygo0NSpU3XffffpueeeU69evQys8o/Nnz+/2u8//vhjBQQEKCkpSYMHD67xmPpAuAIAAADOk91u14lyqyE/293ZUqtDjJOSktSsWTONHj1azz//vMrKyuTi4qJnn31WvXv3lr+/v4KCghQcHPyH93r++edPW1r4v7Zv365WrVrVuL7zlZeXJ0ny9fW9oDF1gXAFAAAAnKcT5VZ1eeoXQ3729mdGqJlLzb/OJycnq3v37urYsaPc3Ny0c+dOubu76+2331ZycrL++c9/1nhJ4IMPPqhbb731nGNCQkJqXNv5stlsmjhxogYMGKBu3bqd95i6Yvieq5kzZyoiIkJubm6KjY3VunXrzjp227ZtuummmxQRESGTyaTp06ef894vvPCCTCaTJk6cWLdFA1BJuVVpx09od2aBklOOacWebCXsyFTCjkwt3pmpJTuztHx3tpJTjmlPZoGO5J1QYWmF7Ha70aUDANAkJScnq2fPnjKZTOrevbu2bNmihx9+WP/3f/+n9u3bKykp6YxLAq3W02fmfH191a5du3M+nJzOHfymTJkik8l0zsfOnTvPeY9x48Zp69at+uKLLy5oTF0xdOZqzpw5mjRpkmbNmqXY2FhNnz5dI0aM0K5duxQQEHDa+OLiYrVp00a33HKLHn744XPee/369XrnnXfUvXv3+iofaPTyTpRrb1aBdmcWandmgQ7kFCkjr0QZ+SU6Xlx+Xvds5mJRkJebAr3cFOTtprAW7mrj76E2fs3Vxt9Dnm7OdfwuAACoP+7OFm1/ZoRhP7s2kpOTdfvtt0uSYmJiNH36dKWmpurzzz9XSUmJdu7cWTVzdd111yksLEzr16/Xn//8Z913333V7lUXywIfeeQR3X333ee8R5s2bc763Pjx4zVv3jwtX75cYWFh5z2mLhkarl577TXdf//9Gjt2rCRp1qxZ+vHHH/Xhhx9qypQpp43v06eP+vTpI0lnfP6UwsJC3XHHHXrvvff03HPP1U/xQCOUklusNQdylXTwmJJSjmlvVuE5xztbTGru6iQPVyc1d3WSq1PlZPipuamyCpuKyipUWFKhgpIKVdjsKi6zan9OkfbnnLlbT4CnqzoFeykq1EtRod7qFuqtUB/3Wq0pBwDgYjGZTLVammeU/fv36/jx41XhqUePHnrrrbf04YcfytPTU2vXrlVFRUXVzNWWLVt06aWX6q233jrj/epiWaC/v7/8/f1r/V7sdrsmTJiguXPnaunSpYqMjDyvMfXBsE9CWVmZkpKSNHXq1KprZrNZcXFxSkxMvKB7jxs3TldffbXi4uJqFK5KS0tVWlpa9fv8/PwL+vlAQ1FhtWndgaNavDNLi3dlaX/26YEn2NtN7QM91T6gudr6N1eIj5uCvd0V5OUmL3enGoeeUxt+M/NLlZFXosz8Eh3JK1HK0WLtzy7U/pwiZReUKqugVFkF2Vq+O7vqtS2aOatX6xbqG+mrvpEt1S3ES04Ww1c1AwDQYCQlJcnFxaVqz9GYMWM0cuRItWzZUlLlrJa/v7/Cw8NVUFAgq9Wqhx566Kz38/X1rffmEGczbtw4zZ49W9999508PT2VkZEhSfL29pa7u3uNx9QHw8JVTk6OrFarAgMDq10PDAz8w7WV5/LFF18oOTlZ69evr/Fr4uPj9fTTT5/3zwQaErvdri1pefp2Y7q+35yunMLf/mLByWxSj1Y+6h3hq16tWqhn6xby9XCpk5976m/2Iv2cFOnnccYx+SXl2pdVqG3p+dqalqctaXnalVGgY8XlWrQjS4t2ZEmqXFrYq3UL9W/rpyEd/NU52JOZLQAAziE5OVndunWTs3Pl8ntnZ2f5+flVe75Hjx6SKvsc9O/f35A6a+Ltt9+WJA0dOrTa9Y8++qhqmWFNxtQHx5/DrIXU1FQ99NBDWrhwodzc3Gr8uqlTp2rSpElVv8/Pz1d4eHh9lAgYprisQt9uTNcnqw9qV2ZB1XVfDxcN7xSg4Z0CNLC9n7wM3PPk5easHq1aqEerFlXXSsqt2nEkXxsOHtPaA0e1/uBR5Z0o14o9OVqxJ0cvzt+pQC9XDengryEdKt+Dtzv7tgAA+L34+HjFx8ef9fn33nuv6tdbtmxRVFTUxSjrvNSkOZZRDbQMC1d+fn6yWCzKzMysdj0zM1NBQUHndc+kpCRlZWVVayFptVq1fPlyzZgxQ6WlpbJYTt/45+rqKldX1/P6mYCjyyoo0QcrDuiL9anKO1HZhMLVyazLuwbphh4hGtTeX84OvMTOzdlSFbjuH9xGNptduzILtGZ/rlbsydHqfTnKzC/Vfzcc1n83HJbFbNIlbXx1RbdgjegaqADPmv9FCwAAqAxXcXFxRpfRIBkWrlxcXNSrVy8lJCRo5MiRkip70CckJGj8+PHndc9LL71UW7ZsqXZt7Nix6tSpk/72t7+dMVgBjVVWQYneWbZf/1lzSKUVNklSK99mGtM/Qjf3Cmuwsztms0mdg73UOdhLYwdEqqTcqvUHj2rprmwt3ZWlfdlFWrU3V6v25uqp77aqd+sWGtE1SFd0C1JYi2ZGlw8AgMN74403jC6hwTJ0WeCkSZM0ZswY9e7dW3379tX06dNVVFRU1T3wrrvuUmhoaNUUZllZmbZv317167S0NG3atEnNmzdXu3bt5OnpedrBYB4eHmrZsmW9HxgGOIqi0gq9tXSvPlh5QCXllaGqZysf/WVoOw3rFCCLuXHtTXJztmhQe38Nau+vJ6/pooM5RZq/LUM/b83Q5tTjWn/wmNYfPKbnftyh3q1b6Iaeobo6Klg+zepmLxkAAMAphoarUaNGKTs7W0899ZQyMjIUExOj+fPnVzW5SElJkdn823Kl9PT0qo12kvTKK6/olVde0ZAhQ7R06dKLXT7gUGw2u+ZuTNOL83cqq6CySUWPVj56OK6DBrX3azINHyL8PPTgkLZ6cEhbpR8/oV9OBq31B49qw6Fj2nDomJ7+fruGdfLXDT1CNaxTgFydmNUGAAAXzmQ3areXA8vPz5e3t7fy8vLk5eVldDnAH9qbVai/ff2rkg4dk1S5/O/vV3fW5V0Cm0yo+iMZeSX6fnOa5m5M144jvx234OXmpGujQzS6byt1C/U2sEIAgKMrKSnRgQMHFBkZWavmaXB85/r/tjbZoFF1CwSamgqrTe+u2K/pi/aorMImDxeLxg9vr3sGRjAb8z+CvN30wOC2emBwW+3MyNfcjWn6bmO6MvJL9NnaFH22NkVRod66rW+4rosOkaeBXRMBAEDDxMzVGTBzhYbgYE6RHvpiozYfzpMkDengr/gboxTiU38H4zU2Vptdifty9cX6FC3Ylqkya+UeNXdni66NDtZtfVupR7gPs38AAEm/zW5ERETU60G0uPhOnDihgwcPMnMFNEU/bE7X1G+2qLC0Ql5uTnrymi66uVcYIaCWLGaTBrb308D2fjpaVKZvkg/r83Up2pddVNXavXOwl8b2j9B1MSFyc2Y2EACaslMH8BYXFxOuGpni4mJJv/1/fL6YuToDZq7gqErKrXpm3nbNXpsiSeob4avXR8co2Jv/wNcVu92u9QeP6Yt1Kfpxy5GqNvYtmjnrtr6t9KdLWiuU2UEAaLKOHDmi48ePKyAgQM2aNeMvNhs4u92u4uJiZWVlycfHR8HBwaeNqU02IFydAeEKjiiroEQPfJqkTanHZTJJ44e100OXtpeTAx8A3NAdLy7TnPWp+jTxkNKOn5AkmU3SiK5BGtM/QrGRvvyhCgBNjN1uV0ZGho4fP250KahDPj4+CgoKOuOf64SrC0S4gqPZmpanBz7doPS8Enm7O+vN0T00uIO/0WU1GVabXYt2ZOrjVQeVuD+36nqnIE/dP6iNro0OkYsTIRcAmhKr1ary8nKjy0AdcHZ2lsVy9qX/hKsLRLiCI1m4PVN//XyjTpRb1cbfQx+M6aNIPw+jy2qydmUU6JPEg/om+XDVIc1BXm66Z2CERvdtRZdBAAAaGcLVBSJcwVF8uSFVf/v6V9ns0qD2fppxe095u/Pl3RHkFZfrs3WH9NGqg8o+eWizp6uTbo9tpbEDIhXkzfknAAA0BoSrC0S4giN4f8V+PffjDknSrb3D9PwNUeyvckClFVZ9tzFd767Yr71ZhZIkZ4tJ10WH6oHBbdQxyNPgCgEAwIUgXF0gwhWMZLfb9a+Fu/XG4r2SpPsHRerxqzrTOMHB2Wx2LdmVpXeW79e6A0errl/aKUDjhrdTz1YtDKwOAACcL8LVBSJcwUjTF+3W9EV7JEmTR3TUX4a2JVg1MBtTjund5fs1f1uGTv0XdmA7P40f3o4OgwAANDCEqwtEuIJRZi7Zq5d/2SVJeuLqzrpvUBuDK8KF2J9dqLeX7tPcjWmqsFX+p7Z36xYaP7ydhnTwJ2QBANAAEK4uEOEKRvj9HqvHruiovwxtZ3BFqCupR4v1zvJ9+u/6wyqzVnYYjAr11vjh7XRZ50CZzYQsAAAcFeHqAhGucLF9nXRYj3y5WZL0cFwHPRTX3uCKUB8y80v03vL9+mxtik6UWyVJHQM9NW54O10dFSwLIQsAAIdDuLpAhCtcTCv2ZGvsR+tVYbPrgcFtNPXKTiwXa+RyC0v14aoD+nT1IRWUVkiS2gc018S4DrqyWxAzWQAAOBDC1QUiXOFi2Zaep1HvrFFhaYWujQ7R66Ni+GLdhOSdKNcnqw/q/RX7lV9SGbI6BXlqYlwHjegaSMgGAMABEK4uEOEKF8ORvBO6fsYqZRWUql+blvr4nj5ydbIYXRYMkHeiXB+uPKAPVx6omsnqGuKlh+M66NLOAYQsAAAMRLi6QIQr1LeScqtufSdRvx7OU4fA5vrywf7ydnc2uiwY7HhxmT44GbKKyir3ZHUP89bDcR00tCPdBQEAMALh6gIRrlCf7Ha7HvnvZn2zMU0+zZz1w/iBCvdtZnRZcCBHi8r03or9+njVwarGFzHhPpp0WQcNau9HyAIA4CIiXF0gwhXq0wcrD+jZedtlMZv06T19NaCdn9ElwUHlFJbq3eX79WniQZWUV7Zw7xPRQpNHdFLfSF+DqwMAoGkgXF0gwhXqS+K+XP3pg7Wy2ux66pouumdgpNEloQHIKijRrKX79dnaQyqtqAxZwzr669ERHdU1xNvg6gAAaNwIVxeIcIX6kFNYqqteX6GsglLd2DNUr94SzfIu1EpGXoleT9ij/25IldVW+Z/ua6NDNOmyDor08zC4OgAAGifC1QUiXKGu2Wx2jf14vZbtzla7gOb6fvwANXNxMrosNFAHcor02sLd+mFzuiTJYjbp1t7heujS9grydjO4OgAAGpfaZAPzRaoJaNLeW7Ffy3Zny9XJrJm39yRY4YJE+nnozdE99ONfB2pYR39ZbXZ9vi5FQ15eoud/2qFjRWVGlwgAQJPEzNUZMHOFurQx5ZhumZWoCptdz98QpdtjWxldEhqZ9QeP6qX5O7X+4DFJkqerk+4f3Eb3DIxUc1eCPAAAF4JlgReIcIW6UlxWoateX6GDucW6pnuw3hzdg31WqBd2u11Ld2XrpV92aceRfElSSw8XjRvWTndc0ooDqgEAOE8sCwQcxEvzd+lgbrGCvNz0zxuiCFaoNyaTScM6BejHCQP1xugeimjZTLlFZXpm3nYNf2WZ5m48LJuNv0sDAKA+Ea6AerJ6X44+Xn1QkvTizd3l7e5sbEFoEsxmk66LDtHCSUMUf2OUAr1clXb8hB6es1nXvLlSy3dnG10iAACNFuEKqAcFJeWa/OWvkqTRfVtpSAd/gytCU+NsMWt031Za+ugwTR7RUZ6uTtp+JF93fbhOf3p/rbam5RldIgAAjQ7hCqgHz/+0U2nHTyishbv+fnVno8tBE+buYtG4Ye207LFhumdApJwtJq3cm6Nr3lyph77YqNSjxUaXCABAo0G4AurYugNH9fm6FEnSyzdH060NDsHXw0VPXdtFix8ZqutjQiRJ321K16WvLtMzP2ynfTsAAHWAcAXUobIKmx6fu0WSdFufcPVr29LgioDqwn2b6fXbemjehIEa2M5PZVabPlx1QINfWqKZS/bqRJnV6BIBAGiwCFdAHXpvxX7tzSpUSw8XTbmyk9HlAGfVLdRb/7kvVp/e01ddgr1UUFqhl3/ZpWGvLNWc9Smy0lkQAIBaI1wBdeRQbpHeSNgjSXrims7yaeZicEXAHxvcwV/zJgzUv0ZFK9THXRn5Jfrb11t0xfTlWrQ9UxyFCABAzRGugDpgt9v15HfbVFph04B2LTUyJtTokoAaM5tNuqFHmBIeGaInru4sb3dn7ckq1H2fbtCod9ZoY8oxo0sEAKBBIFwBdWDB9kwt350tF4tZz17fjcOC0SC5OVt036A2Wv7YMD04pK1cncxad/Cobnhrtf7yWZIO5BQZXSIAAA6NcAVcoNIKq57/aYck6f7BkWrj39zgioAL4+3urClXdtKSR4fqll5hMpmkn7Zk6LLXlmnad1uVW1hqdIkAADgkwhVwgT5adVCHcosV4OmqvwxtZ3Q5QJ0J8XHXy7dE6+eHBmlYR39V2Oz6JPGQhry8lM6CAACcAeEKuADZBaWasXivJOmxKzrJgzOt0Ah1CvLSR2P7avZ9seoW6qXCk50Fh76yRP9dn0pnQQAATiJcARfglV92qbC0QtFh3rqxB00s0Lj1b+en78cN1Ou3xSishbsy80v12Ne/6srXl2vJziw6CwIAmjzCFXCetqXn6b9JqZKkp67tIrOZJhZo/Mxmk66PCa3WWXB3ZqHGfrxet7+3Vr8ePm50iQAAGIZwBZynl+bvkt0uXdM9WL1a+xpdDnBRuTqd7Cw4eZj+PLiNXJzMStyfq+tmrNJfP9+o1KPFRpcIAMBFR7gCzsOa/blatjtbTmaTJo/oaHQ5gGG8mzlr6lWdtfiRIbqxR6hMJun7zeka/upSPTtvu44VlRldIgAAFw3hCqglu92ul+bvlCTd1jdcrVt6GFwRYLywFs302qgY/TB+oAa281O51a4PVh7Q4JeXaNayfSopp7MgAKDxI1wBtbRoR5aSU47Lzdmsvw5vb3Q5gEPpFuqt/9wXq0/v6atOQZ4qKKnQCz/v1PBXlurrpMOy0VkQANCIEa6AWrDa7Hrll12SpHsGRCrAy83gigDHNLiDv3786yC9eku0QrzdlJ5Xoke+3Kyr31ypFXuyjS4PAIB6QbgCauG7TWnalVkgLzcn/XlwW6PLARyaxWzSTb3CtPjRoZpyZSd5ujlpx5F83fnBOt35wVptS88zukQAAOoU4QqooQqrTW8k7JEkPTi0rbybORtcEdAwuDlb9OCQtlo+eZjuHRgpZ4tJK/bk6Jo3V2rSfzcp7fgJo0sEAKBOEK6AGvrh13QdzC1Wi2bOGtMvwuhygAanhYeLnrymixImDdV10SGy26VvktM07JWliv95h/JOlBtdIgAAF8TwcDVz5kxFRETIzc1NsbGxWrdu3VnHbtu2TTfddJMiIiJkMpk0ffr008bEx8erT58+8vT0VEBAgEaOHKldu3bV4ztAU2C12TVj8V5J0n2D2sjD1cngioCGq1XLZnpjdA99N26AYiN9VVZh0zvL9mvIy0v0/or9Kq2gsyAAoGEyNFzNmTNHkyZN0rRp05ScnKzo6GiNGDFCWVlZZxxfXFysNm3a6IUXXlBQUNAZxyxbtkzjxo3TmjVrtHDhQpWXl+vyyy9XUVFRfb4VNHI/bz2ifdlF8nZ31l39WhtdDtAoRIf76IsHLtGHd/dW+4DmOl5crud+3KFLX12m7zal0VkQANDgmOx2u2F/esXGxqpPnz6aMWOGJMlmsyk8PFwTJkzQlClTzvnaiIgITZw4URMnTjznuOzsbAUEBGjZsmUaPHhwjerKz8+Xt7e38vLy5OXlVaPXoPGy2ey68vUV2pVZoIfjOuihONqvA3WtwmrT18mH9eqC3coqKJUkRYV6a+pVndS/rZ/B1QEAmrLaZAPDZq7KysqUlJSkuLi434oxmxUXF6fExMQ6+zl5eZXdqHx9fc86prS0VPn5+dUewCkLtmdqV2aBPF2ddPeACKPLARolJ4tZo/q00tLJQ/Xo5R3U3NVJW9LydPt7a3XPx+u1K6PA6BIBAPhDhoWrnJwcWa1WBQYGVrseGBiojIyMOvkZNptNEydO1IABA9StW7ezjouPj5e3t3fVIzw8vE5+Pho+u92uNxdXdggc0z9C3u50CATqUzMXJ40f3l5LJw/VmH6t5WQ2afHOLF35+nL97atflZFXYnSJAACcleENLerTuHHjtHXrVn3xxRfnHDd16lTl5eVVPVJTUy9ShXB0q/bmalt6vtydLbpnYKTR5QBNhl9zVz19fTctnDREV0UFyWaX5mxI1dBXluiVX3apoITOggAAx2NYuPLz85PFYlFmZma165mZmWdtVlEb48eP17x587RkyRKFhYWdc6yrq6u8vLyqPQBJemf5PknSqD7h8vVwMbgaoOmJ9PPQW3f00tf/11+9W7dQSblNM5bs1ZCXl+qT1QdVVmEzukQAAKoYFq5cXFzUq1cvJSQkVF2z2WxKSEhQv379zvu+drtd48eP19y5c7V48WJFRjLbgPOz40i+VuzJkdkk3cusFWCoXq1b6MsH++ndO3upjb+HjhaVadr323T5v5bppy1HZGBvJgAAqhh6WM+kSZM0ZswY9e7dW3379tX06dNVVFSksWPHSpLuuusuhYaGKj4+XlJlE4zt27dX/TotLU2bNm1S8+bN1a5dO0mVSwFnz56t7777Tp6enlX7t7y9veXu7m7Au0RD9d6K/ZKkK6OCFe7bzOBqAJhMJl3eNUjDOwVozoZU/WvhHh3MLdZfPktWTLiPHr+qs/pGnr15EQAA9c3QVuySNGPGDL388svKyMhQTEyM3njjDcXGxkqShg4dqoiICH388ceSpIMHD55xJmrIkCFaunSppMo/fM/ko48+0t13312jmmjFjiN5JzToxSWqsNn1/fgB6h7mY3RJAP5HYWmF3lu+X+8u368T5ZUHD1/WJVB/u6KT2gU0N7g6AEBjUZtsYHi4ckSEK8T/tEPvLN+v2Ehfzfnz+S9TBVD/svJLND1hj+asT5XVZpfFbNKoPuGaGNdeAZ5uRpcHAGjgGsQ5V4CjKigp1+y1KZKkPw9pY3A1AP5IgJebnr8hSr9MHKTLugTKarNr9toUDX15qaYv2q2i0gqjSwQANBGEK+B/fJV0WAWlFWoX0FxDOwQYXQ6AGmoX4Kn37uqt//65n2LCfVRcZtX0RXs05OWl+mztIVVY6SwIAKhfhCvgd2w2u/6deEiSdHf/CJnNZ97DB8Bx9Y301dy/9Ndbd/RU65bNlFNYqr/P3arLpy/Xgm0ZdBYEANQbwhXwOyv35mh/TpE8XZ10Q49Qo8sBcJ5MJpOuigrWwoeH6OnrusrXw0X7s4v0wL+TdOs7iUpOOWZ0iQCARohwBfzOJ6sPSpJu7h0mD1dDTyoAUAdcnMwa0z9CyyYP1fhh7eTmbNb6g8d041ur9ZfPknQgp8joEgEAjQjhCjgpJbdYi3dlSZLuvKS1wdUAqEuebs56dERHLX10mG7tHSazSfppS4Yue22Zpn23VTmFpUaXCABoBAhXwEn/WXtIdrs0uIO/2vhzRg7QGAV5u+mlm6P100ODNKyjvypsdn2SeEhDX16qGYv36ESZ1egSAQANGOEKkHSizKo561MlSWP6MWsFNHadgrz00di+mn1frLqFeqmwtEKvLNitoa8s0Zz1KbLaaHoBAKg9whUg6fvNaco7Ua5wX3cN7Uj7daCp6N/OT9+PG6jXb4tRWAt3ZeaX6m9fb9GVry/Xkp1ZdBYEANQK4QqQ9J81lYcG/ym2tSy0XweaFLPZpOtjQpXwyBA9cXVnebs7a3dmocZ+vF63v7dWvx4+bnSJAIAGgnCFJm9rWp62pOXJ2WLSLb3DjS4HgEFcnSy6b1AbLZ88TH8e0kYuTmYl7s/VdTNW6a+fb1Tq0WKjSwQAODjCFZq8U3utLu8aJF8PF4OrAWA072bOmnplZy15dKhu7Bkqk0n6fnO6hr+6VM/O265jRWVGlwgAcFCEKzRpJ8qs+nZTmiRpdJ9WBlcDwJGE+rjrtVtjNG/CQA1q76dyq10frDygwS8v0axl+1RSTmdBAEB1hCs0aT9tOaKCkgqFtXBX/7YtjS4HgAPqGuKtf98bq0/v6avOwV4qKKnQCz/v1PBXlurrpMN0FgQAVCFcoUk7tSRwVO9wmWlkAeAcBnfw17wJA/XqLdEK8XZTel6JHvlys658fbl+2ZZBZ0EAAOEKTdferEKtO3hUZpNoZAGgRixmk27qFabFjw7VlCs7ycvNSbszC/Xnfydp5FurtWpvjtElAgAMRLhCk/XfDZWzVsM6BijI283gagA0JG7OFj04pK1W/G24xg9rJ3dnizanHtcd76/VHe+v0caUY0aXCAAwAOEKTVJZhU1fJx2WJN3Wl0YWAM6Pt7uzHh3RUcsfG6a7+0fIxWLWqr25uuGt1br/0w3alVFgdIkAgIuIcIUmaemuLOUWlcnf01XDOvobXQ6ABs7f01X/uK6rFj86RLf0CpPZJC3cnqkrXl+uh+ds0qHcIqNLBABcBIQrNEnfJFe2Xx8ZEyInC/8aAKgbYS2a6eVborXg4SG6OipYdrs0d2OaLn11mf4+d4sy80uMLhEAUI/4Vokm53hxmRbvzJIk3dgzzOBqADRG7QKaa+YdPTVvwkAN6eCvCptdn61N0eCXlij+px0cRAwAjRThCk3OvF+PqMxqU+dgL3UO9jK6HACNWLdQb31yT1/NeeAS9W7dQqUVNr2zfL8Gv7REbyTsUWFphdElAgDqEOEKTc43yZWNLG7sEWpwJQCaitg2LfXlg/300d191CXYSwWlFXpt4W4NeWmJPlh5QCXlVqNLBADUAcIVmpSDOUVKTjkus0m6PibE6HIANCEmk0nDOgVo3oSBmnF7D7Xx81BuUZmenbddw15ZqtlrU1RutRldJgDgAhCu0KR8s7GykcWg9v4K8OJsKwAXn9ls0jXdQ7Tg4cF68aYohXi76UheiR6fu0XDX12qLzekqoKQBQANEuEKTYbdbtfcjSeXBPZkSSAAYzlZzBrVp5UWPzpUT13TRX7NXZV69IQmf/WrLvvXcn27MU1Wm93oMgEAtUC4QpOx4dAxpR49IQ8Xiy7vEmR0OQAgSXJztuiegZFa8dgwPX5VJ/l6uOhATpEmztmkK6Yv14+/HpGNkAUADQLhCk3G3JNLAq+MCpa7i8XgagCgOncXix4Y3FbLHxumySM6ytvdWXuyCjVudrKuemOFftmWIbudkAUAjoxwhSah3GrTz1uOSJJGxrAkEIDjau7qpHHD2mnF34ZpYlx7ebo6aWdGgf787yRdN2OVluzMImQBgIMiXKFJWL0vV8eKy9XSw0WXtPE1uhwA+ENebs6aGNdBK/42TOOGtVUzF4u2pOVp7MfrdePbq7ViTzYhCwAcDOEKTcIPm9MlSVdGBcnJwsceQMPh08xFk0d00orHhunPg9vIzdmsjSnHdecH6zTqnTVK3JdrdIkAgJP4lolGr7TCql+2ZUiSru3O2VYAGqaWzV019arOWv7YMN0zIFIuTmatO3hUo99bo9vfW6MNB48aXSIANHmEKzR6K3bnqKCkQoFeruoTwZJAAA1bgKebnrq2i5ZPHqY7L2ktZ4tJq/fl6uZZibrzg7VKOkTIAgCjEK7Q6P3wa+WSwKuigmU2mwyuBgDqRpC3m54d2U1LHh2q0X3D5WQ2acWeHN30dqL+9P5arWcmCwAuOsIVGrUTZVYt2p4pSbo2miWBABqfsBbNFH9jdy1+ZKhG9a4MWSv35uiWWYm6/b01WrufPVkAcLEQrtCoLdmVpaIyq0J93NUj3MfocgCg3rRq2Uwv3ty92kzW6n25GvXuGt32biKNLwDgIiBcoVGbd3JJ4DXRwTKZWBIIoPEL962cyVo6eajuiG0lZ4tJa/ZXNr649Z1ErdqbQwt3AKgnhCs0WkWlFUrYkSWJLoEAmp6wFs30zxuitOxk4wsXi1nrDhzVHe+v1S2zEjknCwDqAeEKjdbSXdkqrbCpdctm6hriZXQ5AGCIEB93PTuym5Y9NlRj+rWWi5NZGw4d050frNNNb6/Wst2ELACoK4QrNFqnzra6omsQSwIBNHnB3u56+vpuWvHYMI0dECFXJ7OSU45rzIfrdMNbq7VkZxYhCwAuEOEKjVJphVWLd1YuCRzRLcjgagDAcQR6uWnatV214rFhundgpNyczdqUelxjP16vkTNXadH2TEIWAJwnwhUapdV7c1VYWqEAT1fFhPkYXQ4AOJwALzc9eU0XLX9smO4fVBmyNh/O032fbtCVr6/QD5vTZbURsgCgNghXaJROLQkc0TWIg4MB4BwCPN3096u7aOXfhuvPQ9rIw8WinRkFmvD5RsW9tkz/XZ+qsgqb0WUCQINAuEKjY7XZteDkwcFXsCQQAGrEr7mrpl7ZWaumDNfEuPbydnfWgZwiPfb1rxr2ylJ9mnhQJeVWo8sEAIdGuEKjs/7gUR0tKpNPM2f1jfQ1uhwAaFB8mrloYlwHrZoyXFOv7CS/5q5KO35CT323TQNfXKJ3lu1TYWmF0WUCgEMiXKHRObUk8NJOgXK28BEHgPPR3NVJfx7SViv/NkzPXN9VoT7uyiksVfzPOzXghcWavmi3jheXGV0mADgUvnmiUbHb7VqwjSWBAFBX3JwtuqtfhJZOHqqXbu6uNn4eyjtRrumL9mjAC4sV//MOZReUGl0mADgEwhUala1p+Uo7fkLNXCwa1N7P6HIAoNFwtph1a+9wLZw0RG+O7qFOQZ4qKrPqnWX7NfDFxZr23ValHT9hdJkAYCjDw9XMmTMVEREhNzc3xcbGat26dWcdu23bNt10002KiIiQyWTS9OnTL/ieaFzmbzsiSRrWMUBuzhaDqwGAxsdiNuna6BD9/NAgvX9Xb8WE+6i0wqZPEg9pyEtL9NhXm7U3q9DoMgHAEIaGqzlz5mjSpEmaNm2akpOTFR0drREjRigrK+uM44uLi9WmTRu98MILCgo685Kv2t4TjcvCk10CL+8aaHAlANC4mUwmxXUJ1Ny/9Ndn98WqX5uWqrDZ9d8Nh3XZv5bpgU83KDnlmNFlAsBFZbIbeAx7bGys+vTpoxkzZkiSbDabwsPDNWHCBE2ZMuWcr42IiNDEiRM1ceLEOrvnKfn5+fL29lZeXp68vLxq/8ZgiNSjxRr00hJZzCYlP3GZvJs5G10SADQpSYeOadayfVV/0SVJfSN99X9D2mpoR3+ZTJw7CKDhqU02MGzmqqysTElJSYqLi/utGLNZcXFxSkxMdJh7ouFI2FH5h3nv1i0IVgBggF6tW+i9u3pr0aTBuqVXmJwtJq07cFRjP16vK19foW83pqnCyoHEABovw8JVTk6OrFarAgOrL98KDAxURkbGRb1naWmp8vPzqz3Q8CTsrFz6GdeZJYEAYKR2AZ56+ZZoLX9smO4fFCkPF4t2ZhRo4pxNGvLyUn286oBOlHEgMYDGx/CGFo4gPj5e3t7eVY/w8HCjS0ItFZZWaO3+o5Kk4Z0DDK4GACBJwd7u+vvVXbR6yqV69PIOaunhorTjJ/SPH7ar/wsJmr5ot44VcVYWgMbDsHDl5+cni8WizMzMatczMzPP2qyivu45depU5eXlVT1SU1PP6+fDOCt2Z6vMalNEy2Zq4+dhdDkAgN/xbuas8cPba9WU4Xp2ZDe18m2mY8WVZ2X1f2Gxnv5hG23cATQKhoUrFxcX9erVSwkJCVXXbDabEhIS1K9fv4t6T1dXV3l5eVV7oGE5tSTw0s6BbJgGAAfl5mzRnZe01uJHKs/K6hLspRPlVn206qCGvLREk+Zs0q6MAqPLBIDz5mTkD580aZLGjBmj3r17q2/fvpo+fbqKioo0duxYSdJdd92l0NBQxcfHS6psWLF9+/aqX6elpWnTpk1q3ry52rVrV6N7ovGx2uxacipcdWJJIAA4OieLWddGh+ia7sFasSdHs5bt0+p9ufpmY5q+2ZimIR38df+gNhrQriV/YQagQTE0XI0aNUrZ2dl66qmnlJGRoZiYGM2fP7+qIUVKSorM5t8m19LT09WjR4+q37/yyit65ZVXNGTIEC1durRG90Tjs/nwceUWlcnT1Ul9In2NLgcAUEMmk0mDO/hrcAd/bU49rneW79P8rRlatjtby3Znq3Owl+4fFKlruofIxYlt4gAcn6HnXDkqzrlqWF7+ZadmLtmnq7sHa+btPY0uBwBwAQ7lFunDlQf03w2HdaK8sqNgoJer7u4fqdtjW8nbnaM2AFxcDeKcK6CuJOw41YKdJYEA0NC1bumhp6/vpsSpwzV5REcFeLoqM79UL87fqX7xCXr6h21KPVpsdJkAcEbMXJ0BM1cNx+FjxRr44hKZTdKGJy6Tr4eL0SUBAOpQaYVV329K1/srDmhXZmWzC7NJurJbsO4bFKkerVoYXCGAxq422cDQPVfAhVqyK1uS1LNVC4IVADRCrk4W3dI7XDf3CtOKPTl6b8V+rdiTox+3HNGPW46od+sWum9QG13WJVAWM80vABiLcIUGbfnuynA1jC6BANCo/b75xY4j+Xp/xQF9vzlNGw4d04ZDSYpo2Uz3DozUzb3C5e5iMbpcAE0UywLPgGWBDUNZhU09nlmgojKrfhg/UFFh3kaXBAC4iDLzS/TJ6oP6z5pDyi+pkCT5NHPW7X1b6c5+rRXs7W5whQAag9pkA8LVGRCuGobEfbka/d4atfRw0fq/x8nMchAAaJKKSiv05YZUfbDqgFKPnpAkOZlNujIqWGMHRKgn+7IAXAD2XKFJWHZySeDgDv4EKwBowjxcnXT3gEjd2S9CC7dn6MNVB7XuwFH9sDldP2xOV0y4j8YOiNBVUcFyttAoGUD94b8waLBO7bca0sHf4EoAAI7AYjbpim7B+u+f+2nehIG6uVeYXCxmbUo9roe+2KSBLy7WzCV7dbSozOhSATRSLAs8A5YFOr6s/BL1fT5BkrThiTj5NXc1uCIAgCPKLijV7LUp+veaQ8opLJUkuTqZdUOPUI0dEKmOQZ4GVwjA0XGIMBq95XtyJElRod4EKwDAWfl7uuqhuPZaNWWYXrs1WlGh3iqtsOmL9akaMX25bn9vjRZtz5TNxt81A7hw7LlCg7SMJYEAgFpwdbLoxp5huqFHqJIOHdOHqw5o/tYMrd6Xq9X7ctW6ZTPd3T9CN/cKk6ebs9HlAmigCFdocKw2u1buORmuOhKuAAA1ZzKZ1DvCV70jfHX4WLH+nXhIn69L0aHcYj39w3a9umC3bu4Vpjv7tVZb/+ZGlwuggWHP1Rmw58qxbUo9rpEzV8nTzUkbn7xMTnR+AgBcgOKyCn2TnKaPVh3QvuyiquuD2vvprn4RGt4pQBa60gJNFq3Y0agt21U5azWgrR/BCgBwwZq5OOlPl7TW7X1baeXeHH2aeEgJOzO1Yk+OVuzJUaiPu/50SWuN6hMuXw8Xo8sF4MAIV2hwlu3OksSSQABA3TKbTRrcwV+DO/gr9Wix/rP2kOasT1Xa8RN6cf5O/WvRbl3bPURj+rdW9zAfo8sF4IBYFngGLAt0XHnF5erx7ALZ7NKqKcMV6uNudEkAgEaspNyqHzan69PEQ9qSlld1PTrcR2P6tdZVUcFyc7YYWCGA+sayQDRaiftzZbNLbfw9CFYAgHrn5mzRLb3DdXOvMG1KPa5PEw/px1+PaHPqcU1KPa7nftyh2/qE645LWvPnEoDan3P10Ucfqbi4uD5qAf7Q6n2V51sNbOdncCUAgKbEZDKpR6sW+teoGK2eOlyTR3RUsLebjhaV6a2l+zToxcV64NMNWrU3RywKApquWi8LDAwM1IkTJ3TLLbfo3nvvVf/+/eurNsOwLNBxDX91qfZnF2nWn3rpim5BRpcDAGjCKqw2LdqRpU8TD2r1vtyq6238PXRHbGvd1DNUPs1ogAE0dLXJBrWeuUpLS9Mnn3yinJwcDR06VJ06ddKLL76ojIyM8y4YqIkjeSe0P7tIZpPUr01Lo8sBADRxThazrugWpNn3X6KFDw/WXf1ay8PFov3ZRXp23nbFPp+gR/67WUmHjjGbBTQRF9TQIjMzU//5z3/0ySefaOfOnbriiit077336tprr5XZ3HBbZDNz5Zi+SjqsR7/crOgwb303fqDR5QAAcJrC0gp9tylN/1mToh1H8quudwry1B2XtNbImBB5ujkbWCGA2qrXmavfCwwM1MCBA9WvXz+ZzWZt2bJFY8aMUdu2bbV06dILuTVwmtV7K/dbDWC/FQDAQTV3ddIdsa31018Hau5f+uvmXmFydTJrZ0aBnvx2q2KfT9DUb7Zo6+86DwJoPM4rXGVmZuqVV15R165dNXToUOXn52vevHk6cOCA0tLSdOutt2rMmDF1XSuaMLvdrpWEKwBAA3GqAcYrt0Rr3eNxmnZtF7ULaK7iMqs+X5eia95cqetnrNR/16equKzC6HIB1JFaLwu89tpr9csvv6hDhw667777dNddd8nX17famKysLAUFBclms9VpsRcLywIdz96sAsW9tlyuTmZtnnY5Z4oAABocu92udQeO6rO1Kfp56xGVWyu/gnm6OemmnmG6PbaVOgR6GlwlgP9Vr+dcBQQEaNmyZerXr99Zx/j7++vAgQO1vTVwViv3VM5a9Y5oQbACADRIJpNJsW1aKrZNS+UUdtFXSYf1+boUHcot1serD+rj1QfVJ6KF7ohtrSu6BfHnHdAA1XpZ4JAhQ9SzZ8/TrpeVlenTTz+VVPkfj9atW194dcBJq062uGVJIACgMfBr7qoHh7TVkkeG6t/39tUVXYNkMZu0/uAxTZyzSbHPJ+gf32/Tzoz8P74ZAIdR62WBFotFR44cUUBAQLXrubm5CggIkNVqrdMCjcCyQMdSYbWpxzMLVVBaoe/GDVB0uI/RJQEAUOcy80v03/Wp+mJ9qtKOn6i6Hh3uo9v6hOva6BA1d631oiMAF6helwXa7XaZTKbTrh8+fFje3t61vR3wh7ak5amgtEJebk7qFspnDADQOAV6uWnCpe31l2HttHJvjuasT9HC7ZnanHpcm1OP69l523Vt9xCN6huuHuE+Z/w+BsBYNQ5XPXr0kMlkkslk0qWXXionp99earVadeDAAV1xxRX1UiSatlUnuwT2a9tSFjN/kAAAGjeL2aQhHfw1pIO/cgpLNTc5TZ+vT9H+7CLN2ZCqORtS1SGwuUb1aaUbeoTK18PF6JIBnFTjcDVy5EhJ0qZNmzRixAg1b9686jkXFxdFRETopptuqvMCgVV7K/dbDWS/FQCgifFr7qr7B7fRfYMiteHQMX2xLlU/bknX7sxCPTtvu178eacu7xqo2/q0Uv+2LWXmLyEBQ9V6z9Unn3yiUaNGyc3Nrb5qMhx7rhxHSblV3Z9eoLIKmxIeGaK2/s3/+EUAADRi+SXl+n5TuuasT9WW3x1GHNbCXaN6h+uW3uEK8m6839OAi6022aDW4aopIFw5jsR9uRr93hoFeLpq7eOXsr4cAIDf2ZqWpznrU/XtpjQVlFQeRmw2SUM7BuiWXmEa3jlArk60dAcuRJ03tPD19dXu3bvl5+enFi1anPML7tGjR2tXLXAOa/ZXLgmMbdOSYAUAwP/oFuqtbqHeevyqzvp56xF9sT5V6w4c1eKdWVq8M0s+zZw1MiZUN/cKU9cQL/4sBepZjcLVv/71L3l6elb9mn8xcbGsPVAZri5p42twJQAAOC53F4tu7BmmG3uGaV92ob5KOqxvkg8rM7+06oDiTkGeurlXmEb2CJVfc1ejSwYaJZYFngHLAh0D+60AADh/VptdK/Zk66ukw1qwPVNlFTZJkpPZpGGdKpcNDusUIGeL2eBKAcdWr+dcJScny9nZWVFRUZKk7777Th999JG6dOmif/zjH3JxoR0o6sbm1OMqq7DJr7mr2vh5GF0OAAANisVs0tCOARraMUB5xeX6/td0fbUhVZsP52nh9kwt3J6plh4uGtmjctlg52D+Qhm4ULX+q4o///nP2r17tyRp//79GjVqlJo1a6Yvv/xSjz32WJ0XiKZrzf7K/XuXtPFlKSoAABfAu5mz7ryktb4bP1ALHh6sBwa3kV9zV+UWlemDlQd05esrdM2bK/TxqgM6VlRmdLlAg1XrZYHe3t5KTk5W27Zt9eKLL2rx4sX65ZdftGrVKt12221KTU2tr1ovGpYFOobR765R4v5cPTeym/50SWujywEAoFGpsNq0bHflssFFOzJVbq38SuhsMSmuc6Bu7hWmwR38WTaIJq9elwXa7XbZbJVrdhctWqRrrrlGkhQeHq6cnJzzKBc4XWmFVckpxyTRzAIAgPrgZDHr0s6BurRzoI4Wlen7TWn6MumwtqXn6+etGfp5a4Zaerjo2ugQ3dAjVN3DvFlJAvyBWoer3r1767nnnlNcXJyWLVumt99+W5J04MABBQYG1nmBaJo2p+aptMImv+YuNLIAAKCe+Xq46O4Bkbp7QKS2p+frq6TD+n5zmnIKy6q6Dbbx99CNPUI1skeowlo0M7pkwCHVOlxNnz5dd9xxh7799lv9/e9/V7t27SRJX331lfr371/nBaJpWnvqfKtIzrcCAOBi6hLipadCuujxqzppxd4czU1O04LtGdqfXaRXFuzWKwt2q2+kr27sEaoro4Ll7e5sdMmAw6izVuwlJSWyWCxydm74/4Kx58p4d7y/Rqv25urZ67vqzn4RRpcDAECTVlBSrvlbMzR3Y5oS9+fq1LdHFyezLuscqBt6hGpIR/ZnoXGq1z1Xp5SVlSkrK6tq/9UprVq1Ot9bApKksgqbkg6d2m/V0uBqAACAp5uzbukdrlt6hyv9+Al9tyldczce1u7MQv245Yh+3HJEvh4uurZ7sG7oGaZo9mehiar1zNXu3bt17733avXq1dWu2+12mUwmWa3WOi3QCMxcGWvDwaO6eVaiWnq4aMMTcfzHGQAAB2S327UtPV9zN6bpu03pyiksrXqujZ+Hbji5Pyvcl/1ZaNjqdeZq7NixcnJy0rx58xQcHMwXX9S5Naf2W3G+FQAADstkMqlbqLe6hXpr6pWdtHJvjuZuTNMv2zK0P6dIry7crVcX7lbPVj66PiZUV0UFy9/T1eiygXpV63C1adMmJSUlqVOnTvVRD6C1B04dHsySQAAAGgIni1lDOwZoaMcAFZZWnNyfdViJ+3KVnHJcySnH9fQP2zSgnZ+uiw7RiG5B8nJr+Pv0gf9V63DVpUsXzrNCvSm32rThYOV+q9hIwhUAAA1Nc1cn3dwrTDf3ClNWfonm/XpE329O16bU41qxJ0cr9uTo799u1fCOAbo+JkTDOgXIzdlidNlAnaj1nqvFixfriSee0PPPP6+oqKjTugM2hj1K7LkyzubU47p+5ir5NHNW8hOXyWxmWSAAAI3Bodwifb8pXd9tTtferMKq681dnTSia5CuiwnRgLYt5UTHQTiY2mSDWocrs7nyA/+/e2FoaIG68P6K/Xruxx2K6xyg98f0MbocAABQx+x2u3YcKdD3m9P1w+Z0pR0/UfVcSw8XXd09WNfHhKhnqxbsvYZDqNeGFkuWLDnvws5k5syZevnll5WRkaHo6Gi9+eab6tu371nHf/nll3ryySd18OBBtW/fXi+++KKuuuqqqucLCws1ZcoUffvtt8rNzVVkZKT++te/6sEHH6zTulE/Ti0J7B3ha3AlAACgPphMJnUJ8VKXEC89NqKjklOO6btN6fpxyxHlFpXp08RD+jTxkEJ93HVdTIiuiw5RpyBPghYahDo7RPh8zJkzR3fddZdmzZql2NhYTZ8+XV9++aV27dqlgICA08avXr1agwcPVnx8vK655hrNnj1bL774opKTk9WtWzdJ0gMPPKDFixfr/fffV0REhBYsWKC//OUv+uabb3TdddfVqC5mroxht9vV55+LlFNYpq//r596tSZgAQDQVJRbbVq1N0ffb07XL1szVFT222qodgHNdXVUsK7pHqz2gZ4GVommqF6XBUrSihUr9M4772j//v368ssvFRoaqn//+9+KjIzUwIEDa3yf2NhY9enTRzNmzJAk2Ww2hYeHa8KECZoyZcpp40eNGqWioiLNmzev6toll1yimJgYzZo1S5LUrVs3jRo1Sk8++WTVmF69eunKK6/Uc889V6O6CFfGOJBTpGGvLJWLk1lb/nG5XJ3Y3AoAQFNUUm7V4p1Z+m5TmpbszFaZ1Vb1XIfA5ro6KkRXdw9Wu4DmBlaJpqI22aDWOwa//vprjRgxQu7u7kpOTlZpaeWBcXl5eXr++edrfJ+ysjIlJSUpLi7ut2LMZsXFxSkxMfGMr0lMTKw2XpJGjBhRbXz//v31/fffKy0tTXa7XUuWLNHu3bt1+eWXn7WW0tJS5efnV3vg4lt/sgV7TJgPwQoAgCbMzdmiq6KC9c6dvbXhyTj9a1S04joHyNli0u7MQv1r0W7FvbZMV0xfrjcT9mh/duEf3xS4CGodrp577jnNmjVL7733XrVOgQMGDFBycnKN75OTkyOr1arAwMBq1wMDA5WRkXHG12RkZPzh+DfffFNdunRRWFiYXFxcdMUVV2jmzJkaPHjwWWuJj4+Xt7d31SM8PLzG7wN1Z/3BynDVO6KFwZUAAABH4eXmrBt6hOn9MX204YnL9Oot0RrW0V9OZpN2ZhTo1YW7NfzVZbry9RWauWSvDuYUGV0ymrBaN7TYtWvXGYOKt7e3jh8/Xhc1XZA333xTa9as0ffff6/WrVtr+fLlGjdunEJCQk6b9Tpl6tSpmjRpUtXv8/PzCVgG2HCosplFH5pZAACAM/B2d9ZNvcJ0U68w5RWX65ftGfrx1yNatTdHO47ka8eRfL38yy51DfHS1d2DdU1UiFq1bGZ02WhCah2ugoKCtHfvXkVERFS7vnLlSrVp06bG9/Hz85PFYlFmZma165mZmQoKCjrrzz7X+BMnTujxxx/X3LlzdfXVV0uSunfvrk2bNumVV145a7hydXWVq6trjWtH3csuKNWBnCKZTFLPVsxcAQCAc/Nu5qxbe4fr1t7hOlZUpgXbMzTv1yNavS9X29LztS09Xy/N36XuYd66KipYV0cFK9yXoIX6Vetlgffff78eeughrV27ViaTSenp6frss8/06KOP6v/+7/9qfB8XFxf16tVLCQkJVddsNpsSEhLUr1+/M76mX79+1cZL0sKFC6vGl5eXq7y8vOosrlMsFotsNpvguJIOVS4J7BjoKe9mzn8wGgAA4DctPFw0qk8r/fveWK3/e5zib4zSwHZ+MpukXw/n6YWfd2rQS0t0zZuVSwf3sUcL9aTWM1dTpkyRzWbTpZdequLiYg0ePFiurq569NFHNWHChFrda9KkSRozZox69+6tvn37avr06SoqKtLYsWMlSXfddZdCQ0MVHx8vSXrooYc0ZMgQvfrqq7r66qv1xRdfaMOGDXr33XclSV5eXhoyZIgmT54sd3d3tW7dWsuWLdOnn36q1157rbZvFRfR+qrzrZi1AgAA58/Xw0Wj+7bS6L6tlFNYql+2ZWje5iNaeyBXW9PytTWtculgh8DmuqJbsK7oGqTOwZyjhbpx3udclZWVae/evSosLFSXLl3UvPn5tcKcMWNG1SHCMTExeuONNxQbGytJGjp0qCIiIvTxxx9Xjf/yyy/1xBNPVB0i/NJLL1U7RDgjI0NTp07VggULdPToUbVu3VoPPPCAHn744Rr/S0Mr9ovv+hkrtflwnl6/LUbXx4QaXQ4AAGhkcgpLtXB7pn7emqHVe3NUYfvtK3BEy2Ya0S1IV3YLVnSYN0EL1dT7OVd2u125ubkymUxq2bLleRfqqAhXF1dxWYWi/rFAVptdq6YMV6iPu9ElAQCARiyvuFyLdmRq/rYMLdudrbKK37aPhHi7VQWtXq1byGImaDV1tckGtVoWmJGRoccee0zff/+9CgoKJFUuxbvhhhsUHx9/Wpt0oCY2pRyX1WZXiLcbwQoAANQ772a/dR0sKq3Qkl1Z+nlrhpbszFJ6Xok+WnVQH606KL/mrrq8a6Cu7BakS9q0lLOl1u0K0MTUOFzl5+erf//+Kiws1NixY9WpUyfZ7XZt375dn3/+uVauXKnk5OTzXh6Ipuu3/Va0YAcAABeXh6uTrukeomu6h6ik3KoVe3L089YjWrQ9UzmFpZq9NkWz16bI291ZcZ0rg9bA9n5yc7YYXTocUI3D1euvvy6LxaJt27bJ39+/2nNPPPGEBgwYoDfeeEOPP/54nReJxm3DyU6BfWhmAQAADOTmbNFlXQJ1WZdAlVXYlLg/V/O3ZmjBtgzlFpXp6+TD+jr5sJq5WDS4vb8u7xqo4Z0C5NPMxejS4SBqvOfqkksu0Z///OeqTn7/68MPP9R7772nxMTEOi3QCOy5ungqrDZFP71ARWVW/fzQIHUO5p83AABwLFabXesPHtX8rRn6ZVuGjuSVVD1nMZvUN8JXl3etDGVhLThLq7Gpl4YWvr6+SkxMVMeOHc/4/M6dO9W/f38dPXq09hU7GMLVxbM1LU/XvLlSnm5O2vTU5WwaBQAADs1ut2trWr4Wbs/Qgu2Z2plRUO35LsFeuqxLoC7vGqguwV50HmwE6qWhRX5+vnx8fM76vI+Pj/Lz82tcJCBJySmV+616tKIbDwAAcHwmk0lRYd6KCvPWpMs7KiW3WAu2Z2jh9kytP3hU24/ka/uRfL2esEehPu6VQatLoPpE+tIQowmocbiy2+0ym8/+gTCZTDrPI7PQhCUfqgxXPVv5GFsIAADAeWjVspnuG9RG9w1qo6NFZVq8M0sLtmVo+Z5spR0/oY9XH9THqw/K291ZwzsF6PIugRrcwV8errVq2o0GolbhqkOHDmed2iRY4XwkpxyXJPVsRTMLAADQsPl6uOjmXmG6uVeYTpRZtXJvjhZsy1DCziwdLSrT3I1pmrsxTS5OZg1s56e4zoG6tHOAAr3cjC4ddaTG4eqjjz6qzzrQBOUUlirlaLFMJimGmSsAANCIuLv81nnQarMr6dCxqn1ah3KLtXhnlhbvzJLmSlGh3hreKUCXdg5QtxBvmdkq0WDVuKFFU0JDi4tjwbYMPfDvJHUIbK4FDw8xuhwAAIB6Z7fbtTuzUAu3Z2jRjixtPnxcv/82HuDpquGdAjS8U4AGtvdTMxeWDxqtXhpaAHWNJYEAAKCpMZlM6hjkqY5Bnho/vL2yC0q1ZFeWFu/I0oo92coqKNUX61P1xfpUuTiZ1b9tS13aufI8rVAfd6PLxx8gXMEwpzoFEq4AAEBT5e/pqlt7h+vW3uEqrbBq7f6jWrwzS4t2ZOrwsRNauitbS3dl60lJnYI8Fdc5UMM7BygmzIflgw6IZYFnwLLA+ldutSnqH7+opNymRZMGq12Ap9ElAQAAOAy73a49WYVK2JGlhB2ZSk45JtvvvrX7NXfR0I4BiuscoIHt/dWc7oP1hmWBcHg7jxSopNwmLzcntfFrbnQ5AAAADsVkMqlDoKc6BHrq/4a21dGiMi3bnaVFO7K0fFe2cgrL9FXSYX2VdFjOFpNiI1tqaEd/De0YoLb+HhxebJBah6slS5Zo2LBh9VELmpBTSwJjWrVgShsAAOAP+Hq46IYeYbqhR5jKrTatP3BUCTsrZ7UO5hZr5d4crdybo+d+3KGwFu4a2tFfwzoGqF/bljTFuIhqvSzQ1dVVYWFhGjt2rMaMGaPw8PD6qs0wLAusfw99sVHfbUrXxLj2mhjXwehyAAAAGiS73a79OUVasjNLy3Zna+3+oyqz2qqed7GYFdvGV0M7BmhoR3+18WNWq7Zqkw1qHa5ycnL073//W5988om2bdum4cOH695779XIkSPl4uJyQYU7CsJV/Rv00mKlHj2hT+/pq8Ed/I0uBwAAoFEoLqvQ6r25Wro7S0t3ZevwsRPVng/3ddewk0GrXxs/ubtYDKq04ajXcPV7ycnJ+uijj/T5559Lkm6//Xbde++9io6OPt9bOgTCVf3KLihVn38ukskkbZ52ubzcnI0uCQAAoNGx2+3al11Y1XFw3YH/mdVyMuuSNi01tIO/hnb0VySzWmd00cKVJKWnp+vdd9/VCy+8ICcnJ5WUlKhfv36aNWuWunbteiG3Ngzhqn5xeDAAAMDFV1RaocR9uVqyq3JWK+149VmtVr7NNOxkU4zYNr7s1Tqp3rsFlpeX67vvvtOHH36ohQsXqnfv3poxY4ZGjx6t7OxsPfHEE7rlllu0ffv283oDaNw4PBgAAODi83B1UlyXQMV1Caya1VqyM1tLd2dp3YGjSjlarE8SD+mTxENysZjVO6KFBrX31+AOfuoc5EUTshqo9czVhAkT9Pnnn8tut+vOO+/Ufffdp27dulUbk5GRoZCQENlstrPcxbExc1W/bn0nUesOHNVLN3XXrX0aX0MUAACAhqawtEKr9+Zo6e5sLTvDrJZfcxcNbOenwR38NbC9nwI83Qyq9OKr15mr7du3680339SNN94oV1fXM47x8/PTkiVLantrNAHlVpt+PXxcktSztY+htQAAAKBSc1cnXd41SJd3DZLdbteBnCKt2JOjFXuytXpfrnIKy/TtpnR9uyldktQpyFODO/hrcHt/9Y5oITdnGmNI5zFztXz5cvXv319OTtVzWUVFhVavXq3BgwfXaYFGYOaq/mw5nKdrZ6yUl5uTNj11OdPLAAAADq6swqbklGNasSdby3fnaGt6nn6fIFydzIpt01KD21fObLUPaN6oGmPUa0MLi8WiI0eOKCAgoNr13NxcBQQEyGq11r5iB0O4qj+frD6oad9v05AO/vrknr5GlwMAAIBayi0s1ap9uVq+O1sr9mQrM7+02vOBXq4n92r5a2A7P/l6NOzjmup1WaDdbj9jEs3NzZWHh0dtb4cmJjnlmCSaWQAAADRULZu76rroEF0XHSK73a49WYVavjtby/fkaO3+XGXml+qrpMP6KumwTCapW4i3BrTz08B2fo1+CWGNw9WNN94oSTKZTLr77rur7beyWq369ddf1b9//7qvEI3KptTjkqQerXwMrQMAAAAXzmQyqUOgpzoEeuq+QW1UUm7V+oNHtWJPjpbvztbOjAJtScvTlrQ8zVq2Ty5OZvWJaKEB7fw0oK2fuoV6y9KItonUOFx5e3tLqpy58vT0lLu7e9VzLi4uuuSSS3T//ffXfYVoNI4VlelQbrEkKTrMx9hiAAAAUOfcnC0a1N5fg9r76/GrOisrv0Sr9uVo5Z5crdqbo4z8Eq3am6tVe3Ml7ZKXm5P6t/XTgPaVM1sRLZs16P1aNQ5XH330kSQpIiJCjz76KEsAUWubT3YJbOPnIe9mzsYWAwAAgHoX4OWmG3qE6YYeYSfP1irSqr05Wrk3R2v25Sq/pELzt2Vo/rYMSVKoj7sGtGupAe381L+tn/w9z9yd3FHVuqFFU0BDi/oxfdFuTV+0Rzf0CNW/RsUYXQ4AAAAMVGG16de0PK0+GbaSDh1TufW3aOLp6qRN0y43fNlgnTe06NmzpxISEtSiRQv16NHjnFN1ycnJtasWTcbmk/utYsJ9DK0DAAAAxnOymNWzVQv1bNVC44e3V3FZhdYfPFY5s7UnR6Et3A0PVrVVo3B1/fXXVzWwGDlyZH3Wg0bKbrdr8+E8SVI04QoAAAD/o5mLk4Z08NeQDv6SKme2GhqWBZ4BywLrXurRYg16aYmcLSZtfXqEXJ0abwtOAAAANB61yQbmi1QTmrhTLdi7BHsRrAAAANAo1WhZYIsWLWrcEvHo0aMXVBAap1P7rVgSCAAAgMaqRuFq+vTp9VwGGrtTbdg53woAAACNVY3C1ZgxY+q7DjRiFVabtqTRzAIAAACNW43CVX5+ftXmrfz8/HOOpQEE/teuzAKVlNvk6eakNn4cPg0AAIDGqcZ7ro4cOaKAgAD5+Piccf+V3W6XyWSS1Wqt8yLRsG1OPTlrFeYjcwM7qwAAAACoqRqFq8WLF8vX11eStGTJknotCI3Pb80svI0tBAAAAKhHNQpXQ4YMOeOvgZqgmQUAAACaghqFq/917NgxffDBB9qxY4ckqUuXLho7dmzV7BZwSlFphXZnFkiSYmhmAQAAgEas1ocIL1++XBEREXrjjTd07NgxHTt2TG+88YYiIyO1fPny+qgRDdjWtDzZ7FKwt5sCvNyMLgcAAACoN7WeuRo3bpxGjRqlt99+WxaLRZJktVr1l7/8RePGjdOWLVvqvEg0XCwJBAAAQFNR65mrvXv36pFHHqkKVpJksVg0adIk7d27t06LQ8O36WQzi5hWPobWAQAAANS3Woernj17Vu21+r0dO3YoOjq6TopC4/H7NuwAAABAY1ajZYG//vpr1a//+te/6qGHHtLevXt1ySWXSJLWrFmjmTNn6oUXXqifKtEgZRWUKO34CZlMUlQYbdgBAADQuJnsdrv9jwaZzWaZTCb90dDGcohwfn6+vL29lZeXJy8vL6PLabAWbc/UfZ9uUIfA5lrwMC38AQAA0PDUJhvUaObqwIEDdVIYmhaaWQAAAKApqVG4at26dX3XgUboVDOLaM63AgAAQBNQ64YWp2zfvl3z58/X999/X+1RWzNnzlRERITc3NwUGxurdevWnXP8l19+qU6dOsnNzU1RUVH66aefThuzY8cOXXfddfL29paHh4f69OmjlJSUWteG82e32/XrYZpZAAAAoOmo9TlX+/fv1w033KAtW7ZU24dlMpkkqVZ7rubMmaNJkyZp1qxZio2N1fTp0zVixAjt2rVLAQEBp41fvXq1Ro8erfj4eF1zzTWaPXu2Ro4cqeTkZHXr1k2StG/fPg0cOFD33nuvnn76aXl5eWnbtm1yc+MA24sp5Wix8k6Uy8XJrI5BnkaXAwAAANS7GjW0+L1rr71WFotF77//viIjI7Vu3Trl5ubqkUce0SuvvKJBgwbV+F6xsbHq06ePZsyYIUmy2WwKDw/XhAkTNGXKlNPGjxo1SkVFRZo3b17VtUsuuUQxMTGaNWuWJOm2226Ts7Oz/v3vf9fmbVVDQ4sL98PmdE34fKOiw3303bgBRpcDAAAAnJfaZINaLwtMTEzUM888Iz8/P5nNZpnNZg0cOFDx8fH661//WuP7lJWVKSkpSXFxcb8VYzYrLi5OiYmJZ/3Zvx8vSSNGjKgab7PZ9OOPP6pDhw4aMWKEAgICFBsbq2+//factZSWlio/P7/aAxdma1rlksCoUMIpAAAAmoZahyur1SpPz8plXn5+fkpPT5dU2fRi165dNb5PTk6OrFarAgMDq10PDAxURkbGGV+TkZFxzvFZWVkqLCzUCy+8oCuuuEILFizQDTfcoBtvvFHLli07ay3x8fHy9vaueoSHh9f4feDMTu236h7qY2whAAAAwEVS6z1X3bp10+bNmxUZGanY2Fi99NJLcnFx0bvvvqs2bdrUR401ZrPZJEnXX3+9Hn74YUlSTEyMVq9erVmzZmnIkDOftTR16lRNmjSp6vf5+fkErAtgs9l/m7ni8GAAAAA0EbUOV0888YSKiookSc8884yuueYaDRo0SC1bttScOXNqfB8/Pz9ZLBZlZmZWu56ZmamgoKAzviYoKOic4/38/OTk5KQuXbpUG9O5c2etXLnyrLW4urrK1dW1xrXj3A7mFqmgtEKuTma1D2hudDkAAADARVHrZYEjRozQjTfeKElq166ddu7cqZycHGVlZWn48OE1vo+Li4t69eqlhISEqms2m00JCQnq16/fGV/Tr1+/auMlaeHChVXjXVxc1KdPn9OWJ+7evZuzui6iLSdnrbqGeMnJct7d/gEAAIAGpdYzV7+XmpoqSee9hG7SpEkaM2aMevfurb59+2r69OkqKirS2LFjJUl33XWXQkNDFR8fL0l66KGHNGTIEL366qu6+uqr9cUXX2jDhg169913q+45efJkjRo1SoMHD9awYcM0f/58/fDDD1q6dOmFvFXUQtV+K863AgAAQBNS62mFiooKPfnkk/L29lZERIQiIiLk7e2tJ554QuXl5bW616hRo/TKK6/oqaeeUkxMjDZt2qT58+dXNa1ISUnRkSNHqsb3799fs2fP1rvvvqvo6Gh99dVX+vbbb6vOuJKkG264QbNmzdJLL72kqKgovf/++/r66681cODA2r5VnKdTM1fdQtlvBQAAgKaj1udc/d///Z+++eYbPfPMM1XL8RITE/WPf/xDI0eO1Ntvv10vhV5MnHN1/qw2u7r/4xcVlVm14OHB6hDIAcIAAABouGqTDWq9LHD27Nn64osvdOWVV1Zd6969u8LDwzV69OhGEa5w/g7kFKqozCp3Z4va+tPMAgAAAE1HrZcFurq6KiIi4rTrkZGRcnFxqYua0ICd2m/VLdRLFrPJ4GoAAACAi6fW4Wr8+PF69tlnVVpaWnWttLRU//znPzV+/Pg6LQ4Nz6lwFcXhwQAAAGhiarQs8FTr9VMWLVqksLAwRUdHS5I2b96ssrIyXXrppXVfIRqUU4cHd+fwYAAAADQxNQpX3t7VvyjfdNNN1X5/vq3Y0bhUWG3alp4viU6BAAAAaHpqFK4++uij+q4DjcC+7CKdKLfKw8WiNn4eRpcDAAAAXFTnfYhwdna2du3aJUnq2LGj/P3966woNEy/Hj4uqXLWykwzCwAAADQxtW5oUVRUpHvuuUfBwcEaPHiwBg8erJCQEN17770qLi6ujxrRQGxhvxUAAACasFqHq0mTJmnZsmX64YcfdPz4cR0/flzfffedli1bpkceeaQ+akQDcSpcRYX5GFsIAAAAYIBaLwv8+uuv9dVXX2no0KFV16666iq5u7vr1ltv5RDhJqrcatP2k80somhmAQAAgCao1jNXxcXFCgwMPO16QEAAywKbsD2ZhSqtsMnTzUmtfZsZXQ4AAABw0dU6XPXr10/Tpk1TSUlJ1bUTJ07o6aefVr9+/eq0ODQcW9KOS6qctaKZBQAAAJqiWi8LnD59uq644orTDhF2c3PTL7/8UucFomH49fCp/VYsCQQAAEDTVOtwFRUVpT179uizzz7Tzp07JUmjR4/WHXfcIXd39zovEA3D1lOdAkN9jC0EAAAAMEitwlV5ebk6deqkefPm6f7776+vmtDAlFXYtONIgSTasAMAAKDpqtWeK2dn52p7rQBJ2p1ZoDKrTd7uzgprwewlAAAAmqZaN7QYN26cXnzxRVVUVNRHPWiATu236h7mLZOJZhYAAABommq952r9+vVKSEjQggULFBUVJQ8Pj2rPf/PNN3VWHBqG33cKBAAAAJqqWocrHx8f3XTTTfVRCxqoLWm/zVwBAAAATVWtw9VHH31UH3WggSopt2pXRmUzi6gwH2OLAQAAAAxU4z1XNptNL774ogYMGKA+ffpoypQpOnHiRH3WhgZgV0aByq12tfRwUYi3m9HlAAAAAIapcbj65z//qccff1zNmzdXaGioXn/9dY0bN64+a0MD8OvJJYHdQmlmAQAAgKatxuHq008/1VtvvaVffvlF3377rX744Qd99tlnstls9VkfHNyWw8clsd8KAAAAqHG4SklJ0VVXXVX1+7i4OJlMJqWnp9dLYWgYtqTlS6qcuQIAAACashqHq4qKCrm5Vd9T4+zsrPLy8jovCg1DSblVezJPNrMgXAEAAKCJq3G3QLvdrrvvvluurq5V10pKSvTggw9WO+uKc66ajl0ZBaqwVTazCKaZBQAAAJq4GoerMWPGnHbtT3/6U50Wg4ZlC80sAAAAgCo1Dlecb4X/tbUqXHkZXAkAAABgvBrvuQL+16mZK/ZbAQAAAIQrnKfSCqt2n2xmQadAAAAAgHCF87Q7o1DlVrtaNHNWqI+70eUAAAAAhiNc4bzQzAIAAACojnCF8/L7cAUAAACAcIXztJVmFgAAAEA1hCvUWlmFTbsyKptZEK4AAACASoQr1NruzAKVWW3ydndWWAuaWQAAAAAS4Qrn4feHB9PMAgAAAKhEuEKt0cwCAAAAOB3hCrW2NT1fktQthHAFAAAAnEK4Qq2UW23acaQyXNHMAgAAAPgN4Qq1siezUGUVNnm6Oal1y2ZGlwMAAAA4DMIVaqWqmUWIN80sAAAAgN8hXKFWtqb/1ikQAAAAwG8IV6gVOgUCAAAAZ0a4Qo1V0MwCAAAAOCvCFWpsb3ahSsptau7qpIiWHkaXAwAAADgUwhVqbGta5axVlxAvmc00swAAAAB+j3CFGjvVKZAlgQAAAMDpHCJczZw5UxEREXJzc1NsbKzWrVt3zvFffvmlOnXqJDc3N0VFRemnn34669gHH3xQJpNJ06dPr+Oqm54thCsAAADgrAwPV3PmzNGkSZM0bdo0JScnKzo6WiNGjFBWVtYZx69evVqjR4/Wvffeq40bN2rkyJEaOXKktm7detrYuXPnas2aNQoJCanvt9HoWW12bU+vXBZIp0AAAADgdIaHq9dee03333+/xo4dqy5dumjWrFlq1qyZPvzwwzOOf/3113XFFVdo8uTJ6ty5s5599ln17NlTM2bMqDYuLS1NEyZM0GeffSZnZ+eL8VYatf3ZhTpRblUzF4si/WhmAQAAAPwvQ8NVWVmZkpKSFBcXV3XNbDYrLi5OiYmJZ3xNYmJitfGSNGLEiGrjbTab7rzzTk2ePFldu3b9wzpKS0uVn59f7YHqTi0J7BriJQvNLAAAAIDTGBqucnJyZLVaFRgYWO16YGCgMjIyzviajIyMPxz/4osvysnJSX/9619rVEd8fLy8vb2rHuHh4bV8J40fhwcDAAAA52b4ssC6lpSUpNdff10ff/yxTKaazbBMnTpVeXl5VY/U1NR6rrLhoVMgAAAAcG6Ghis/Pz9ZLBZlZmZWu56ZmamgoKAzviYoKOic41esWKGsrCy1atVKTk5OcnJy0qFDh/TII48oIiLijPd0dXWVl5dXtQd+Y7PZtY1mFgAAAMA5GRquXFxc1KtXLyUkJFRds9lsSkhIUL9+/c74mn79+lUbL0kLFy6sGn/nnXfq119/1aZNm6oeISEhmjx5sn755Zf6ezON2P6cIhWXWeXubFFb/+ZGlwMAAAA4JCejC5g0aZLGjBmj3r17q2/fvpo+fbqKioo0duxYSdJdd92l0NBQxcfHS5IeeughDRkyRK+++qquvvpqffHFF9qwYYPeffddSVLLli3VsmXLaj/D2dlZQUFB6tix48V9c43EqSWBXWhmAQAAAJyV4eFq1KhRys7O1lNPPaWMjAzFxMRo/vz5VU0rUlJSZDb/NsHWv39/zZ49W0888YQef/xxtW/fXt9++626detm1Fto9Dg8GAAAAPhjJrvdbje6CEeTn58vb29v5eXlsf9K0q3vJGrdgaN6+ebuuqU3nRQBAADQdNQmGzS6boGoWzabXdtPNrOICmPmCgAAADgbwhXO6WBukQpLK+TqZFY7mlkAAAAAZ0W4wjmd2m/VOdhLThY+LgAAAMDZ8G0Z58ThwQAAAEDNEK5wTlvTTu63IlwBAAAA50S4wlnZ7XZtTa+cuepGuAIAAADOiXCFszqUW6yCkgq5OJnVPpBmFgAAAMC5EK5wVlXNLII85UwzCwAAAOCc+MaMs2JJIAAAAFBzhCucFZ0CAQAAgJojXOGM7HZ7VadAZq4AAACAP0a4whmlHj2hvBPlcrGY1SHQ0+hyAAAAAIdHuMIZndpv1THIUy5OfEwAAACAP8K3ZpzRqU6B3UK9DK4EAAAAaBgIVzijXw8flyR1D/MxtA4AAACgoSBc4TR2u12/HqZTIAAAAFAbhCuc5lBusQpKKuTiZFbHIJpZAAAAADVBuMJpNp9cEtgl2EvOFj4iAAAAQE3wzRmn2XJySWD3MJYEAgAAADVFuMJpfk1jvxUAAABQW4QrVGO12bX1ZLiKDvcxthgAAACgASFcoZr92YUqLrPK3dmitv7NjS4HAAAAaDAIV6jmVAv2bqFesphNBlcDAAAANByEK1SzJe1UMwsfYwsBAAAAGhjCFao51YadToEAAABA7RCuUKXcatP29HxJdAoEAAAAaotwhSp7MgtVWmGTp5uTIlp6GF0OAAAA0KAQrlDl15NLAqNCvWWmmQUAAABQK4QrVKk6PJj9VgAAAECtEa5QZcvJNuzRdAoEAAAAao1wBUlSaYVVOzNoZgEAAACcL8IVJEk7jxSo3GpXi2bOCmvhbnQ5AAAAQINDuIKk3/ZbdQ/zkclEMwsAAACgtghXkCT9mnpcEocHAwAAAOeLcAVJ0pZTnQLZbwUAAACcF8IVdKLMqt2ZBZKk6HAfY4sBAAAAGijCFbQ1PU82uxTg6apALzejywEAAAAaJMIVtCnluCQphlkrAAAA4LwRrqBNJ5tZxLTyMbQOAAAAoCEjXOG3cMXMFQAAAHDeCFdNXFZBidKOn5DJVHnGFQAAAIDzQ7hq4k7tt+oQ4Knmrk7GFgMAAAA0YISrJo4lgQAAAEDdIFw1cTSzAAAAAOoG4aoJs9rs+vVwniRmrgAAAIALRbhqwvZlF6qwtELNXCzqEOhpdDkAAABAg0a4asJONbOICvWWxWwythgAAACggSNcNWEb2W8FAAAA1BnCVRN2qplFD/ZbAQAAABfMIcLVzJkzFRERITc3N8XGxmrdunXnHP/ll1+qU6dOcnNzU1RUlH766aeq58rLy/W3v/1NUVFR8vDwUEhIiO666y6lp6fX99toUIrLKrQrI1+SFBPewuBqAAAAgIbP8HA1Z84cTZo0SdOmTVNycrKio6M1YsQIZWVlnXH86tWrNXr0aN17773auHGjRo4cqZEjR2rr1q2SpOLiYiUnJ+vJJ59UcnKyvvnmG+3atUvXXXfdxXxbDm/L4TzZ7FKQl5uCvN2MLgcAAABo8Ex2u91uZAGxsbHq06ePZsyYIUmy2WwKDw/XhAkTNGXKlNPGjxo1SkVFRZo3b17VtUsuuUQxMTGaNWvWGX/G+vXr1bdvXx06dEitWrX6w5ry8/Pl7e2tvLw8eXl5nec7c2yzlu3TCz/v1BVdgzTrzl5GlwMAAAA4pNpkA0NnrsrKypSUlKS4uLiqa2azWXFxcUpMTDzjaxITE6uNl6QRI0acdbwk5eXlyWQyycfH54zPl5aWKj8/v9qjsTvVKZBmFgAAAEDdMDRc5eTkyGq1KjAwsNr1wMBAZWRknPE1GRkZtRpfUlKiv/3tbxo9evRZk2Z8fLy8vb2rHuHh4efxbhoOu92upJRjkmhmAQAAANQVw/dc1afy8nLdeuutstvtevvtt886burUqcrLy6t6pKamXsQqL77Dx04ou6BUzhaToglXAAAAQJ1wMvKH+/n5yWKxKDMzs9r1zMxMBQUFnfE1QUFBNRp/KlgdOnRIixcvPuf6SFdXV7m6up7nu2h4Nhw6KknqGuItN2eLwdUAAAAAjYOhM1cuLi7q1auXEhISqq7ZbDYlJCSoX79+Z3xNv379qo2XpIULF1YbfypY7dmzR4sWLVLLli3r5w00UBsOVi4J7NWaFuwAAABAXTF05kqSJk2apDFjxqh3797q27evpk+frqKiIo0dO1aSdNdddyk0NFTx8fGSpIceekhDhgzRq6++qquvvlpffPGFNmzYoHfffVdSZbC6+eablZycrHnz5slqtVbtx/L19ZWLi4sxb9SBJB2qDFe9CVcAAABAnTE8XI0aNUrZ2dl66qmnlJGRoZiYGM2fP7+qaUVKSorM5t8m2Pr376/Zs2friSee0OOPP6727dvr22+/Vbdu3SRJaWlp+v777yVJMTEx1X7WkiVLNHTo0IvyvhxVfkm5dmUWSJJ6RRCuAAAAgLpi+DlXjqgxn3O1bHe2xny4Tq18m2n5Y8OMLgcAAABwaA3mnCtcfEkHK5tZsN8KAAAAqFuEqybm1PlWhCsAAACgbhGumpAKq00bU45Lknqz3woAAACoU4SrJmRnRoGKy6zydHNShwBPo8sBAAAAGhXCVROy4eR+q56tWshsNhlcDQAAANC4EK6akKSTSwLZbwUAAADUPcJVE2G327XuQK4k9lsBAAAA9YFw1UQcyi1WZn6pXCxm9WxFuAIAAADqGuGqiVizv3LWKibcR27OFoOrAQAAABofwlUTsfZAZTOL2Da+BlcCAAAANE6EqybAbrdXzVxd0qalwdUAAAAAjRPhqglIPXpCR/JK5Gwxsd8KAAAAqCeEqybg1KxVdJiP3F3YbwUAAADUB8JVE7DmZAt29lsBAAAA9Ydw1QSs3V/ZzIL9VgAAAED9IVw1cqlHi5V2/ISczCb1as1+KwAAAKC+EK4aucST+62iwrzVzMXJ4GoAAACAxotw1cit3JMjSRrYzs/gSgAAAIDGjXDViNlsdq3aS7gCAAAALgbCVSO2IyNfuUVlauZiUQ/OtwIAAADqFeGqEVtxcklgvzYt5eLE/9UAAABAfeIbdyNWtd+qPUsCAQAAgPpGuGqkSsqtWnew8nyrQYQrAAAAoN4RrhqpdQeOqqzCpmBvN7X1b250OQAAAECjR7hqpFbsyZZU2SXQZDIZXA0AAADQ+BGuGqnFO7MkSYM7+BtcCQAAANA0EK4aoYM5RdqXXSQns4lwBQAAAFwkhKtGKOHkrFWfCF95uzsbXA0AAADQNBCuGqHFOzMlSZd2DjC4EgAAAKDpIFw1Mvkl5Vq7v7IFe1znQIOrAQAAAJoOwlUjs3x3tipsdrXx91CEn4fR5QAAAABNBuGqkUnYUbnfilkrAAAA4OIiXDUiJeVWLdpeud/q8i6EKwAAAOBiIlw1Iiv25KigtELB3m7q2aqF0eUAAAAATQrhqhH58dd0SdJVUcEym00GVwMAAAA0LYSrRqKk3KqFJ5cEXt092OBqAAAAgKaHcNVILNudraIyq0J93NUj3MfocgAAAIAmh3DVSHyTfFiSdFVUkEwmlgQCAAAAFxvhqhHIKSytasF+c69wg6sBAAAAmibCVSPw7cY0Vdjsig7zVscgT6PLAQAAAJokwlUDZ7fb9d8NqZKkW3ozawUAAAAYhXDVwG04dEy7Mwvl6mTWtdEhRpcDAAAANFmEqwbugxUHJEk39AiVt7uzwdUAAAAATRfhqgFLPVqsBdszJEn3DIw0uBoAAACgaSNcNWAfrDwgm10a1N5PHQJpZAEAAAAYiXDVQKUdP6HZa1MkSQ8OaWtwNQAAAAAIVw3UG4v2qMxqU782LdW/bUujywEAAACaPMJVA7Qp9bj+m1TZfv3RER1kMpkMrggAAAAA4aqBKauwaeo3W2S3SyNjQtSrta/RJQEAAACQg4SrmTNnKiIiQm5uboqNjdW6devOOf7LL79Up06d5ObmpqioKP3000/Vnrfb7XrqqacUHBwsd3d3xcXFac+ePfX5Fi6a537crh1H8uXt7qwnrulidDkAAAAATjI8XM2ZM0eTJk3StGnTlJycrOjoaI0YMUJZWVlnHL969WqNHj1a9957rzZu3KiRI0dq5MiR2rp1a9WYl156SW+88YZmzZqltWvXysPDQyNGjFBJScnFelt1zm63a/qi3fo08ZAk6bVbo+XX3NXgqgAAAACcYrLb7XYjC4iNjVWfPn00Y8YMSZLNZlN4eLgmTJigKVOmnDZ+1KhRKioq0rx586quXXLJJYqJidGsWbNkt9sVEhKiRx55RI8++qgkKS8vT4GBgfr444912223/WFN+fn58vb2Vl5enry8vOronZ6fTanHtTerUHM3HtaqvbmSpMev6qQHBtMhEAAAAKhvtckGhs5clZWVKSkpSXFxcVXXzGaz4uLilJiYeMbXJCYmVhsvSSNGjKgaf+DAAWVkZFQb4+3trdjY2LPes7S0VPn5+dUejmLyl5v16JebtWpvrlyczHp2ZDeCFQAAAOCAnIz84Tk5ObJarQoMDKx2PTAwUDt37jzjazIyMs44PiMjo+r5U9fONuZ/xcfH6+mnnz6v91Dfeke0kL+nq7qFeuv2vq0U4edhdEkAAAAAzsDQcOUopk6dqkmTJlX9Pj8/X+Hh4QZW9Jv4G7sbXQIAAACAGjB0WaCfn58sFosyMzOrXc/MzFRQUNAZXxMUFHTO8af+tzb3dHV1lZeXV7UHAAAAANSGoeHKxcVFvXr1UkJCQtU1m82mhIQE9evX74yv6devX7XxkrRw4cKq8ZGRkQoKCqo2Jj8/X2vXrj3rPQEAAADgQhm+LHDSpEkaM2aMevfurb59+2r69OkqKirS2LFjJUl33XWXQkNDFR8fL0l66KGHNGTIEL366qu6+uqr9cUXX2jDhg169913JUkmk0kTJ07Uc889p/bt2ysyMlJPPvmkQkJCNHLkSKPeJgAAAIBGzvBwNWrUKGVnZ+upp55SRkaGYmJiNH/+/KqGFCkpKTKbf5tg69+/v2bPnq0nnnhCjz/+uNq3b69vv/1W3bp1qxrz2GOPqaioSA888ICOHz+ugQMHav78+XJzc7vo7w8AAABA02D4OVeOyJHOuQIAAABgnAZzzhUAAAAANBaEKwAAAACoA4QrAAAAAKgDhCsAAAAAqAOEKwAAAACoA4QrAAAAAKgDhCsAAAAAqAOEKwAAAACoA4QrAAAAAKgDhCsAAAAAqAOEKwAAAACoA4QrAAAAAKgDhCsAAAAAqANORhfgiOx2uyQpPz/f4EoAAAAAGOlUJjiVEc6FcHUGBQUFkqTw8HCDKwEAAADgCAoKCuTt7X3OMSZ7TSJYE2Oz2ZSeni5PT0+ZTCZDa8nPz1d4eLhSU1Pl5eVlaC1oGPjMoLb4zKC2+MygtvjMoLYc6TNjt9tVUFCgkJAQmc3n3lXFzNUZmM1mhYWFGV1GNV5eXoZ/sNCw8JlBbfGZQW3xmUFt8ZlBbTnKZ+aPZqxOoaEFAAAAANQBwhUAAAAA1AHClYNzdXXVtGnT5OrqanQpaCD4zKC2+MygtvjMoLb4zKC2GupnhoYWAAAAAFAHmLkCAAAAgDpAuAIAAACAOkC4AgAAAIA6QLgCAAAAgDpAuHJwM2fOVEREhNzc3BQbG6t169YZXRIcVHx8vPr06SNPT08FBARo5MiR2rVrl9FloYF44YUXZDKZNHHiRKNLgYNLS0vTn/70J7Vs2VLu7u6KiorShg0bjC4LDshqterJJ59UZGSk3N3d1bZtWz377LOilxp+b/ny5br22msVEhIik8mkb7/9ttrzdrtdTz31lIKDg+Xu7q64uDjt2bPHmGJrgHDlwObMmaNJkyZp2rRpSk5OVnR0tEaMGKGsrCyjS4MDWrZsmcaNG6c1a9Zo4cKFKi8v1+WXX66ioiKjS4ODW79+vd555x11797d6FLg4I4dO6YBAwbI2dlZP//8s7Zv365XX31VLVq0MLo0OKAXX3xRb7/9tmbMmKEdO3boxRdf1EsvvaQ333zT6NLgQIqKihQdHa2ZM2ee8fmXXnpJb7zxhmbNmqW1a9fKw8NDI0aMUElJyUWutGZoxe7AYmNj1adPH82YMUOSZLPZFB4ergkTJmjKlCkGVwdHl52drYCAAC1btkyDBw82uhw4qMLCQvXs2VNvvfWWnnvuOcXExGj69OlGlwUHNWXKFK1atUorVqwwuhQ0ANdcc40CAwP1wQcfVF276aab5O7urv/85z8GVgZHZTKZNHfuXI0cOVJS5axVSEiIHnnkET366KOSpLy8PAUGBurjjz/WbbfdZmC1Z8bMlYMqKytTUlKS4uLiqq6ZzWbFxcUpMTHRwMrQUOTl5UmSfH19Da4EjmzcuHG6+uqrq/23Bjib77//Xr1799Ytt9yigIAA9ejRQ++9957RZcFB9e/fXwkJCdq9e7ckafPmzVq5cqWuvPJKgytDQ3HgwAFlZGRU+zPK29tbsbGxDvt92MnoAnBmOTk5slqtCgwMrHY9MDBQO3fuNKgqNBQ2m00TJ07UgAED1K1bN6PLgYP64osvlJycrPXr1xtdChqI/fv36+2339akSZP0+OOPa/369frrX/8qFxcXjRkzxujy4GCmTJmi/Px8derUSRaLRVarVf/85z91xx13GF0aGoiMjAxJOuP34VPPORrCFdAIjRs3Tlu3btXKlSuNLgUOKjU1VQ899JAWLlwoNzc3o8tBA2Gz2dS7d289//zzkqQePXpo69atmjVrFuEKp/nvf/+rzz77TLNnz1bXrl21adMmTZw4USEhIXxe0GixLNBB+fn5yWKxKDMzs9r1zMxMBQUFGVQVGoLx48dr3rx5WrJkicLCwowuBw4qKSlJWVlZ6tmzp5ycnOTk5KRly5bpjTfekJOTk6xWq9ElwgEFBwerS5cu1a517txZKSkpBlUERzZ58mRNmTJFt912m6KionTnnXfq4YcfVnx8vNGloYE49Z23IX0fJlw5KBcXF/Xq1UsJCQlV12w2mxISEtSvXz8DK4OjstvtGj9+vObOnavFixcrMjLS6JLgwC699FJt2bJFmzZtqnr07t1bd9xxhzZt2iSLxWJ0iXBAAwYMOO2Ih927d6t169YGVQRHVlxcLLO5+ldNi8Uim81mUEVoaCIjIxUUFFTt+3B+fr7Wrl3rsN+HWRbowCZNmqQxY8aod+/e6tu3r6ZPn66ioiKNHTvW6NLggMaNG6fZs2fru+++k6enZ9VaZG9vb7m7uxtcHRyNp6fnafvxPDw81LJlS/bp4awefvhh9e/fX88//7xuvfVWrVu3Tu+++67effddo0uDA7r22mv1z3/+U61atVLXrl21ceNGvfbaa7rnnnuMLg0OpLCwUHv37q36/YEDB7Rp0yb5+vqqVatWmjhxop577jm1b99ekZGRevLJJxUSElLVUdDR0Irdwc2YMUMvv/yyMjIyFBMTozfeeEOxsbFGlwUHZDKZznj9o48+0t13331xi0GDNHToUFqx4w/NmzdPU6dO1Z49exQZGalJkybp/vvvN7osOKCCggI9+eSTmjt3rrKyshQSEqLRo0frqaeekouLi9HlwUEsXbpUw4YNO+36mDFj9PHHH8tut2vatGl69913dfz4cQ0cOFBvvfWWOnToYEC1f4xwBQAAAAB1gD1XAAAAAFAHCFcAAAAAUAcIVwAAAABQBwhXAAAAAFAHCFcAAAAAUAcIVwAAAABQBwhXAAAAAFAHCFcAAAAAUAcIVwCAJmfo0KGaOHFinY39/Ri73a4HHnhAvr6+MplM2rRp0wXVCgBoOAhXAACHdPfdd8tkMslkMsnZ2VmRkZF67LHHVFJSYnRpp/nmm2/07LPPSpLmz5+vjz/+WPPmzdORI0fUrVu3WoU5AEDD5WR0AQAAnM0VV1yhjz76SOXl5UpKStKYMWNkMpn04osvGl1aNb6+vlW/3rdvn4KDg9W/f38DKwIAGIGZKwCAw3J1dVVQUJDCw8M1cuRIxcXFaeHChVXP22w2xcfHKzIyUu7u7oqOjtZXX31V7R5FRUW666671Lx5cwUHB+vVV1+t9vxXX32lqKgoubu7q2XLloqLi1NRUVG1MTabTY899ph8fX0VFBSkf/zjH9WePzUzdffdd2vChAlKSUmRyWRSRESE7r77bi1btkyvv/561UzcwYMHz/h+n3/++aoxv39Mnz79vP8ZAgAuHsIVAKBB2Lp1q1avXi0XF5eqa/Hx8fr00081a9Ysbdu2TQ8//LD+9Kc/admyZVVjJk+erGXLlum7777TggULtHTpUiUnJ0uSjhw5otGjR+uee+7Rjh07tHTpUt14442y2+3VfvYnn3wiDw8PrV27Vi+99JKeeeaZaiHvlNdff13PPPOMwsLCdOTIEa1fv16vv/66+vXrp/vvv19HjhzRkSNHFB4efsb3OGHChKoxR44c0f3336/WrVvr5ptvrot/hACAesayQACAw5o3b56aN2+uiooKlZaWymw2a8aMGZKk0tJSPf/881q0aJH69esnSWrTpo1Wrlypd955R0OGDFFhYaE++OAD/ec//9Gll14qqTIohYWFSaoMVxUVFbrxxhvVunVrSVJUVNRpdXTv3l3Tpk2TJLVv314zZsxQQkKCLrvssmrjvL295enpKYvFoqCgoKrrLi4uatasWbVrZ+Lp6SlPT09J0pNPPlkVBk/VCwBwbIQrAIDDGjZsmN5++20VFRXpX//6l5ycnHTTTTdJkvbu3avi4uLTAk5ZWZl69OghqXL/U1lZmWJjY6ue9/X1VceOHSVJ0dHRuvTSSxUVFaURI0bo8ssv180336wWLVpUu2f37t2r/T44OFhZWVl1/n5Peeqpp/Tvf/9bS5cuVURERL39HABA3SJcAQAcloeHh9q1aydJ+vDDDxUdHa0PPvhA9957rwoLCyVJP/74o0JDQ6u9ztXVtUb3t1gsWrhwoVavXq0FCxbozTff1N///netXbtWkZGRVeOcnZ2rvc5kMslms13IWzuradOm6dNPPyVYAUADxJ4rAECDYDab9fjjj+uJJ57QiRMn1KVLF7m6uiolJUXt2rWr9ji1p6lt27ZydnbW2rVrq+5z7Ngx7d69u+r3JpNJAwYM0NNPP62NGzfKxcVFc+fOrdPaXVxcZLVa/3DctGnT9MknnxCsAKCBYuYKANBg3HLLLZo8ebJmzpypRx99VI8++qgefvhh2Ww2DRw4UHl5eVq1apW8vLw0ZswYNW/eXPfee68mT56sli1bKiAgQH//+99lNlf+3eLatWuVkJCgyy+/XAEBAVq7dq2ys7PVuXPnOq07IiJCa9eu1cGDB9W8eXP5+vpW1XDKc889p7ffflvff/+93NzclJGRIUlq0aJFjWfiAADGIlwBABoMJycnjR8/Xi+99JL+7//+T88++6z8/f0VHx+v/fv3y8fHRz179tTjjz9e9ZqXX35ZhYWFuvbaa+Xp6alHHnlEeXl5kiQvLy8tX75c06dPV35+vlq3bq1XX31VV155ZZ3W/eijj2rMmDHq0qWLTpw4oQMHDlSbmbLb7Xr55ZeVn59f1ZzjlHXr1qlPnz51Wg8AoH6Y7P/bbxYAAAAAUGvsuQIAAACAOkC4AgAAAIA6QLgCAAAAgDpAuAIAAACAOkC4AgAAAIA6QLgCAAAAgDpAuAIAAACAOkC4AgAAAIA6QLgCAAAAgDpAuAIAAACAOkC4AgAAAIA6QLgCAAAAgDrw/wRvvh5Nk3W/AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "M_target = -22\n",
    "M_star = -20.83\n",
    "z_min = 0.0\n",
    "z_max = 10.0\n",
    "num_z = 10000\n",
    "z = np.linspace(z_min, z_max, num_z)\n",
    "# 体积元\n",
    "dV_dz = cosmo.differential_comoving_volume(z).value\n",
    "n_M = schechter(M_target, M_star, phi_star, alpha)\n",
    "\n",
    "h = 0.7\n",
    "n_M = n_M / h**3\n",
    "sky_area = 4 * np.pi\n",
    "N_z = n_M * dV_dz * sky_area  # 星系数目 / z\n",
    "prob_z = N_z / np.trapz(N_z, z)\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(z, prob_z, label=r'$M_r = -22$')\n",
    "plt.xlabel('Redshift $z$')\n",
    "plt.ylabel('Probability Density')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 第二章 作业2\n",
    "\n",
    "某旋涡星系和椭圆星系，其观测到的I波段的视星等都是18等；半光度半径内的表面亮度都是20mag/arcsec^2；其中旋涡星系（图像椭率为0.5 ）观测得到的HI的速度展宽是200km/s，椭圆星系观测到的中心速度弥散度是200km/s，请问这两个星系的距离分别是多少？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "旋涡星系距离： 347.51666546986405 Mpc\n"
     ]
    }
   ],
   "source": [
    "m_spiral = 18\n",
    "sfb_spiral = 20\n",
    "sigma_HI_spiral = 200\n",
    "\n",
    "# 椭率e=0.5, 倾角sin i = 1-e = 0.5\n",
    "W_spiral = sigma_HI_spiral / 0.5\n",
    "# TF关系\n",
    "# 根据课件里的近似取a=2.4/-20.4, b=-19\n",
    "M_spiral = (2.4/(-20.4)) * np.log(W_spiral) - 19\n",
    "\n",
    "d_spiral = 10**((m_spiral - M_spiral)/5 + 1)\n",
    "print(\"旋涡星系距离：\", (d_spiral*u.pc).to(u.Mpc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "椭圆星系距离： 1584.8931924611109 Mpc\n"
     ]
    }
   ],
   "source": [
    "# 椭圆星系绝对星等根据表面亮度和绝对星等的基本面得到大致约为-23等\n",
    "m_ellipse = 18\n",
    "M_ellipse = -23\n",
    "d_ellipse = 10**((m_ellipse - M_ellipse)/5 + 1)\n",
    "print(\"椭圆星系距离：\", (d_ellipse*u.pc).to(u.Mpc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. 旋涡星系的 Tully-Fisher关系以及椭圆星系的Fundamental Plane的内秉弥散在0.2个星等左右，请问以上距离计算的误差是多少？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "主要计算公式就是$10^{(5\\times x)}$进行误差传递，x为星等差，单位为mag，基本面给出弥散为0.2mag\n",
    "\n",
    "对这个公式求导得到$5 \\times\\ln(10)\\times 10^{(5x)}\\times \\Delta x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "椭圆星系距离误差： 145.97405755794787 Mpc\n",
      "旋涡星系距离误差： 32.00746773911631 Mpc\n"
     ]
    }
   ],
   "source": [
    "def error_prop(x):\n",
    "    return np.log(10) * 10**x * (0.2/5) \n",
    "\n",
    "x_ellipse = (m_ellipse - M_ellipse) / 5 + 1\n",
    "delta_x = 0.2 / 5  # 假设 Δ(m - M) = 0.2\n",
    "print(\"椭圆星系距离误差：\", (error_prop(x_ellipse) * u.pc).to(u.Mpc))\n",
    "\n",
    "x_spiral = (m_spiral - M_spiral) / 5 + 1\n",
    "print(\"旋涡星系距离误差：\", (error_prop(x_spiral) * u.pc).to(u.Mpc))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ai4galaxy",
   "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.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}