{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "作业3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy.units as u\n",
    "# 这里的G是以Mpc为单位的引力常数\n",
    "G = 4.301e-9 * (u.km**2 * u.Mpc / u.Msun / u.s**2)  # km^2 Mpc Msun^-1 s^-2\n",
    "def compute(mass, r_eff, N):\n",
    "    r_eff = r_eff.to(u.Mpc)\n",
    "    vel = np.sqrt(G * mass / r_eff).to(u.km/u.s)\n",
    "    sigma = np.sqrt(G * mass / (5 * r_eff)).to(u.km/u.s)\n",
    "    t_cross = (r_eff / vel).to(u.year)\n",
    "    t_relax = (t_cross * N / (6 * np.log(N / 2))).to(u.year)\n",
    "    return sigma, t_cross, t_relax\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Globular cluster: sigma=0.293 km / s, t_cross=1.490946e+09 yr, t_relax=1.893645e+13 yr\n"
     ]
    }
   ],
   "source": [
    "# https://en.wikipedia.org/wiki/Globular_cluster\n",
    "gc_r_eff = 1 * u.kpc  # 典型有效半径\n",
    "gc_mass = 1e5 * u.Msun  # 星团质量\n",
    "gc_n = 1e6  # 粗略估计恒星数量，使用 parsec\n",
    "gc_sigma, gc_t_cross, gc_t_relax = compute(gc_mass, gc_r_eff, gc_n)\n",
    "print(f\"Globular cluster: sigma={gc_sigma:.3f}, t_cross={gc_t_cross:2e}, t_relax={gc_t_relax:2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 星团的半径10pc量级，恒星数目 10^5量级"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elliptical galaxy: sigma=29.329 km / s, t_cross=1.490946e+08 yr, t_relax=1.008679e+17 yr\n"
     ]
    }
   ],
   "source": [
    "# https://en.wikipedia.org/wiki/Elliptical_galaxy\n",
    "eg_r_eff = 10 * u.kpc  # 典型有效半径\n",
    "eg_mass = 1e10 * u.Msun  # 星团质量\n",
    "eg_n = 1e11  # 粗略估计恒星数量，使用 parsec\n",
    "eg_sigma, eg_t_cross, eg_t_relax = compute(eg_mass, eg_r_eff, eg_n)\n",
    "print(f\"Elliptical galaxy: sigma={eg_sigma:.3f}, t_cross={eg_t_cross:2e}, t_relax={eg_t_relax:2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 椭圆星系的速度弥散度在100km/s 左右，10^10太阳质量是小质量，半径没有10kpc大"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Galaxy cluster: sigma=2.933 km / s, t_cross=1.490946e+11 yr, t_relax=7.877837e+22 yr\n"
     ]
    }
   ],
   "source": [
    "# https://en.wikipedia.org/wiki/Galaxy_cluster\n",
    "gc_r_eff = 1 * u.Mpc  # 典型有效半径\n",
    "gc_mass = 1e10 * u.Msun  # 星团质量\n",
    "gc_n = 1e14  # 粗略估计恒星数量，使用 parsec\n",
    "gc_sigma, gc_t_cross, gc_t_relax = compute(gc_mass, gc_r_eff, gc_n)\n",
    "print(f\"Galaxy cluster: sigma={gc_sigma:.3f}, t_cross={gc_t_cross:2e}, t_relax={gc_t_relax:2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 星系团的粒子（星系）数1000个左右，速度弥散度在1000km/s量级"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "作业4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5304.542055981422 km2 kg kpc / (solMass m3)\n",
      "0.08231756814038453\n"
     ]
    }
   ],
   "source": [
    "from astropy.constants import G\n",
    "V_c = 220 * u.km/u.s\n",
    "R_d = 3 * u.kpc # 指数盘\n",
    "M = 5.8e11 * u.Msun # https://zh.wikipedia.org/zh-cn/%E9%93%B6%E6%B2%B3%E7%B3%BB\n",
    "E = - M * V_c**2 / 2 # 维利定理计算总能量\n",
    "Sigma_0 = 500 * u.Msun / u.pc**2 # 假设表面密度\n",
    "M_d = 2 * np.pi * Sigma_0 * R_d**2\n",
    "m_d = M_d / M # (恒星+气体)盘质量占总质量的比例\n",
    "j_d = m_d.copy()\n",
    "J_d = 2 * M_d * R_d * V_c\n",
    "J = J_d / j_d\n",
    "\n",
    "galaxy_lambda = J * np.abs(E) ** 0.5 * G**(-1) * M**(-5/2) \n",
    "print(galaxy_lambda)\n",
    "print(galaxy_lambda.to(u.dimensionless_unscaled))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## 面密度怎么来的？"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "stylegan3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}