{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b4afbead",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Globular Cluster': {'Velocity Dispersion (km/s)': np.float64(6.556655229839465),\n",
       "  'Crossing Time (Myr)': np.float64(1.4922853830119807),\n",
       "  'Relaxation Time (Gyr)': np.float64(1.62022826816746)},\n",
       " 'Elliptical Galaxy': {'Velocity Dispersion (km/s)': np.float64(293.22253597901096),\n",
       "  'Crossing Time (Myr)': np.float64(16.684257791220492),\n",
       "  'Relaxation Time (Gyr)': np.float64(8233955.787929997)},\n",
       " 'Galaxy Cluster': {'Velocity Dispersion (km/s)': np.float64(463.62553749217267),\n",
       "  'Crossing Time (Myr)': np.float64(4220.820455173343),\n",
       "  'Relaxation Time (Gyr)': np.float64(1636677707844.781)}}"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Constants\n",
    "G = 6.67e-11  # Gravitational constant (m^3 kg^-1 s^-2)\n",
    "Msun = 1.989e30  # Solar mass (kg)\n",
    "pc = 3.086e16  # Parsec (m)\n",
    "kpc = 1e3 * pc  # Kiloparsec (m)\n",
    "Mpc = 1e6 * pc  # Megaparsec (m)\n",
    "\n",
    "# Typical values\n",
    "systems = {\n",
    "    \"Globular Cluster\": {\"M\": 1e5 * Msun, \"R\": 10 * pc},\n",
    "    \"Elliptical Galaxy\": {\"M\": 1e11 * Msun, \"R\": 5 * kpc},\n",
    "    \"Galaxy Cluster\": {\"M\": 1e14 * Msun, \"R\": 2 * Mpc},\n",
    "}\n",
    "\n",
    "# Calculations\n",
    "def velocity_dispersion(M, R):\n",
    "    return np.sqrt(G * M / R)  # m/s\n",
    "\n",
    "def crossing_time(R, sigma_v):\n",
    "    return R / sigma_v  # seconds\n",
    "\n",
    "def relaxation_time(N, t_cross):\n",
    "    return (N / (8 * np.log(N))) * t_cross  # seconds\n",
    "\n",
    "# Stellar mass (assume 1 solar mass per star)\n",
    "m_star = Msun\n",
    "\n",
    "# Calculate for each system\n",
    "results = {}\n",
    "for system, params in systems.items():\n",
    "    M = params[\"M\"]\n",
    "    R = params[\"R\"]\n",
    "    \n",
    "    # Number of stars\n",
    "    N = M / m_star\n",
    "    \n",
    "    # Velocity dispersion\n",
    "    sigma_v = velocity_dispersion(M, R)\n",
    "    \n",
    "    # Crossing time\n",
    "    t_cross = crossing_time(R, sigma_v)\n",
    "    \n",
    "    # Relaxation time\n",
    "    t_relax = relaxation_time(N, t_cross)\n",
    "    \n",
    "    results[system] = {\n",
    "        \"Velocity Dispersion (km/s)\": sigma_v / 1e3,  # Convert to km/s\n",
    "        \"Crossing Time (Myr)\": t_cross / (3.154e13),  # Convert to million years\n",
    "        \"Relaxation Time (Gyr)\": t_relax / (3.154e16),  # Convert to billion years\n",
    "    }\n",
    "\n",
    "results\n",
    "\n",
    "#这些数值表明，从球状星团到星系团，随着系统质量和尺度的增大，速度弥散度增加，\n",
    "#而穿越时标和弛豫时标也显著增加，尤其是星系团的弛豫时标极为漫长。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "195e5d7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "### SSY： 星系团的粒子是星系不是恒星"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ab0971bf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1972.9083155243254"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#计算银河系的角动量，自旋参数\n",
    "# 设定常数\n",
    "G = 4.3009e-6  # 万有引力常数，单位：kpc (km/s)^2 / M_sun\n",
    "Vc = 200  # 圆盘速度，单位：km/s\n",
    "Rd = 3  # 指数盘半径，单位：kpc\n",
    "Rvir = 200  # 晕半径，单位：kpc\n",
    "M = 1e12  # 银河系总质量，单位：M_sun\n",
    "md = 0.1  # 假设盘质量比例 (盘质量占总质量的比例)\n",
    "jd = 0.1  # 假设盘角动量比例 (盘角动量占总角动量的比例)\n",
    "\n",
    "# 计算盘的角动量 J_d\n",
    "Md = md * M  # 盘质量 M_d = m_d * M\n",
    "Jd = Md * Rd * Vc  # 盘的角动量 J_d = M_d * R_d * V_c\n",
    "\n",
    "# 计算总角动量 J\n",
    "J = Jd / jd  # 假设 m_d ~ j_d，因此 J = J_d / j_d\n",
    "\n",
    "# 计算银河系的总能量 E\n",
    "E = -G * M**2 / Rvir  # 系统的总引力势能 E = - G * M^2 / R_vir\n",
    "\n",
    "# 计算自旋参数 λ\n",
    "lambda_param = J * (Rvir**0.5) / (G * M**(3/2))\n",
    "\n",
    "# 输出自旋参数\n",
    "lambda_param"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ccc3ed4",
   "metadata": {},
   "outputs": [],
   "source": [
    "### SSY:想想lambda的物理含义,肯定是个小于1的数，哪儿错了？"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "base"
  },
  "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}