{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3a8c36bc",
   "metadata": {},
   "source": [
    "作业16\n",
    "如果某质量为m的天体，从无穷远处到被黑洞吸积到最小稳定轨道（3Rs），释放的能量是多少?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "757b34c6",
   "metadata": {},
   "source": [
    "$R_s = \\frac{2GM}{c^2}$\n",
    "\n",
    "无穷远处的能量：$E=mc^2$\n",
    "\n",
    "在圆轨道$3R_s$处的能量为：？\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f19d9596",
   "metadata": {},
   "source": [
    "如果一个星系中的AGN光度超过恒星光度，其黑洞质量与恒星质量的比值最低是多少？"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8af31d05",
   "metadata": {},
   "source": [
    "$L_E \\sim 30000 \\times \\frac{M_\\text{BH}}{M_{\\odot}}L_{\\odot}$\n",
    "\n",
    "$L_* = \\frac{M_*}3\\times L_\\odot$\n",
    "\n",
    "假设$L_{AGN}\\approx L_{Edd} > L_*$\n",
    "\n",
    "$30000\\times \\frac{M_\\text{BH}}{M_{\\odot}}L_\\odot > \\dfrac{M_*}{3}L_\\odot$\n",
    "\n",
    "得到$\\frac{M_\\text{BH}}{M_*}>\\dfrac{1}{90000}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30abbd6b",
   "metadata": {},
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.cosmology import Planck18 as cosmo\n",
    "from astropy.io import fits\n",
    "from scipy.ndimage import zoom\n",
    "\n",
    "# 假设我们有一个低红移星系的 FITS 图像\n",
    "low_redshift_image_file = '149.74358_2.125144_sci.fits'\n",
    "hdul = fits.open(low_redshift_image_file)\n",
    "image_data = hdul[1].data\n",
    "header = hdul[0].header\n",
    "# 显示原始图像\n",
    "plt.imshow(image_data, cmap='gray', origin='lower')\n",
    "plt.colorbar(label='Flux')\n",
    "plt.title('raw iamge')\n",
    "plt.show()\n",
    "\n",
    "def k_correction(z_old, z_new, wave_length_rest):\n",
    "    \"\"\"\n",
    "    简单的 K-修正计算函数（假设一个粗略的模型）。\n",
    "    \n",
    "    参数:\n",
    "    - z_old: 原始红移\n",
    "    - z_new: 目标红移\n",
    "    - wave_length_rest: 本征波长（单位：微米）\n",
    "    \n",
    "    返回:\n",
    "    - k_corr: K-修正因子\n",
    "    \"\"\"\n",
    "    # 简化的 K-修正：可以使用一个线性模型或经验公式\n",
    "    delta_z = z_new - z_old\n",
    "    k_corr = 2.5 * np.log10(1 + delta_z)  # 一个非常粗略的估计\n",
    "\n",
    "    return k_corr\n",
    "\n",
    "def apply_redshift_and_k_correction(image_data, z_old, z_new, wave_length_rest):\n",
    "    \"\"\"\n",
    "    将星系从低红移 z_old 移动到更高红移 z_new，并应用红移效应和 K-修正。\n",
    "    \n",
    "    参数:\n",
    "    - image_data: 原始星系图像数据 (2D numpy array)\n",
    "    - z_old: 原始红移\n",
    "    - z_new: 目标红移\n",
    "    - wave_length_rest: 本征波长（单位：微米）\n",
    "    \n",
    "    返回:\n",
    "    - new_image_data: 应用红移效应和 K-修正后的星系图像数据\n",
    "    \"\"\"\n",
    "    # 计算表面亮度衰减因子\n",
    "    surface_brightness_factor = (1 + z_old)**4 / (1 + z_new)**4\n",
    "\n",
    "    # 计算 K-修正\n",
    "    k_corr = k_correction(z_old, z_new, wave_length_rest)\n",
    "\n",
    "    # 调整图像亮度，考虑表面亮度衰减和 K-修正\n",
    "    new_image_data = image_data * surface_brightness_factor * 10**(-0.4 * k_corr)\n",
    "\n",
    "    # 缩放空间分辨率\n",
    "    scale_factor = (1 + z_old) / (1 + z_new)\n",
    "    from scipy.ndimage import zoom\n",
    "    new_image_data = zoom(new_image_data, scale_factor)\n",
    "\n",
    "    return new_image_data\n",
    "\n",
    "# 设置原始和目标红移\n",
    "z_old = 0.05  # 低红移\n",
    "z_new = 1.0   # 更高红移\n",
    "wave_length_rest = 0.5  # 假设的本征波长（单位：微米）\n",
    "\n",
    "# 应用红移效应和 K-修正\n",
    "new_image_data = apply_redshift_and_k_correction(image_data, z_old, z_new, wave_length_rest)\n",
    "\n",
    "# 显示应用红移效应和 K-修正后的星系图像\n",
    "plt.imshow(new_image_data, cmap='gray', origin='lower')\n",
    "plt.colorbar(label='Flux')\n",
    "plt.title(f'new image')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77ed93b9",
   "metadata": {},
   "source": [
    "task_18. 对于作业5或者作业15中的星系，采取至少一种方法来（定量）评估该星系所处的环境"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4f7c13d",
   "metadata": {},
   "source": [
    "要评估一个星系所处的环境，给定星系的坐标（赤经 RA、赤纬 DEC）和红移 ( z )，我们通常需要计算该星系在宇宙中的位置，并分析其周围其他星系的分布密度。这里我们可以通过以下步骤来评估该星系所处的环境：\n",
    "\n",
    "1. 通过红移估算星系的距离\n",
    "通过红移 ( z )，我们可以计算星系的距离。最简单的距离计算方法是基于“哈勃定律”或者使用一个宇宙学模型。哈勃定律公式为： [ D = \\frac{c \\times z}{H_0} ] 其中：\n",
    "\n",
    "( D ) 是星系的距离，单位为 Mpc（百万秒差距）。\n",
    "( c ) 是光速，约为 300,000 km/s。\n",
    "( z ) 是红移。\n",
    "( H_0 ) 是哈勃常数，常用值为 ( H_0 = 70 , \\text{km/s/Mpc} )。\n",
    "2. 将星系坐标从天球坐标转换为空间坐标\n",
    "使用星系的赤经和赤纬坐标，我们可以将它们转换为三维空间坐标。结合估算的距离（通过红移），我们可以得到星系的三维空间位置。\n",
    "\n",
    "3. 计算局部密度\n",
    "为了评估星系的环境，可以通过计算它周围一定半径内的邻近星系的数量（或质量密度）来进行局部密度估算。这些星系的分布可以帮助我们了解目标星系所在的密集或稀疏区域。\n",
    "\n",
    "代码实现\n",
    "假设你有目标星系的坐标（赤经 RA，赤纬 DEC）和红移 ( z )，我们可以按照以下步骤来实现这一过程：\n",
    "\n",
    "1. 计算星系距离\n",
    "首先，计算星系的距离。我们可以使用简单的哈勃定律来进行距离估算。\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "# 计算距离的函数\n",
    "def calculate_distance(z, H0=70, c=300000):\n",
    "    # 计算哈勃定律距离（单位：Mpc）\n",
    "    return (c * z) / H0\n",
    "\n",
    "# 目标星系的红移\n",
    "z = 0.03  # 假设红移为 0.03\n",
    "\n",
    "# 计算距离\n",
    "distance = calculate_distance(z)\n",
    "print(f\"目标星系的距离为 {distance:.2f} Mpc\")\n",
    "2. 将天球坐标转换为三维空间坐标\n",
    "假设星系的坐标是赤经 RA 和赤纬 DEC，我们需要将其转换为三维空间坐标。\n",
    "\n",
    "# 将天球坐标 (RA, DEC) 转换为笛卡尔坐标 (x, y, z)\n",
    "def sky_to_cartesian(ra, dec, z):\n",
    "    # 计算星系的距离\n",
    "    distance = calculate_distance(z)\n",
    "\n",
    "    # 将角度转换为弧度\n",
    "    ra_rad = np.radians(ra)\n",
    "    dec_rad = np.radians(dec)\n",
    "\n",
    "    # 计算笛卡尔坐标\n",
    "    x = distance * np.cos(dec_rad) * np.cos(ra_rad)\n",
    "    y = distance * np.cos(dec_rad) * np.sin(ra_rad)\n",
    "    z = distance * np.sin(dec_rad)\n",
    "\n",
    "    return np.array([x, y, z])\n",
    "\n",
    "# 目标星系的坐标\n",
    "ra = 150  # 假设赤经 RA 为 150 度\n",
    "dec = 2.5  # 假设赤纬 DEC 为 2.5 度\n",
    "\n",
    "# 转换为三维坐标\n",
    "position = sky_to_cartesian(ra, dec, z)\n",
    "print(f\"目标星系的三维空间坐标为 {position}\")\n",
    "3. 计算局部密度\n",
    "为了评估星系的环境，我们需要知道周围的星系分布。通常可以通过以下方式来评估局部密度：\n",
    "\n",
    "使用附近的邻近星系数目：在给定半径范围内计算邻近星系的数目。这个方法通常使用一个球形半径（例如 1 Mpc）来计算星系的局部密度。\n",
    "我们可以通过给定半径范围，计算目标星系附近的其他星系数量。\n",
    "\n",
    "from scipy.spatial import KDTree\n",
    "\n",
    "# 假设我们有多个星系的三维空间坐标\n",
    "# 这里使用假数据作为示例\n",
    "ra_list = [150, 151, 152, 160, 170]  # 假设的星系赤经\n",
    "dec_list = [2.5, 3.0, 2.8, 1.5, 2.6]  # 假设的星系赤纬\n",
    "z_list = [0.03, 0.04, 0.05, 0.02, 0.03]  # 假设的星系红移\n",
    "\n",
    "# 计算所有星系的三维坐标\n",
    "positions = np.array([sky_to_cartesian(ra, dec, z) for ra, dec, z in zip(ra_list, dec_list, z_list)])\n",
    "\n",
    "# 使用 KDTree 计算星系的最近邻\n",
    "tree = KDTree(positions)\n",
    "\n",
    "# 目标星系的索引\n",
    "target_index = 0\n",
    "\n",
    "# 查询目标星系周围 1 Mpc 半径内的 10 个最近邻星系\n",
    "distances, indices = tree.query(positions[target_index], k=10)\n",
    "\n",
    "# 计算在 1 Mpc 内的星系数量\n",
    "radius = 1.0  # 半径为 1 Mpc\n",
    "neighbors_within_radius = sum(distances < radius)\n",
    "print(f\"目标星系周围 {radius} Mpc 内的星系数量：{neighbors_within_radius}\")\n",
    "评估环境：\n",
    "如果目标星系周围的星系数量很高，表明该星系位于一个密集的区域，可能是星系团的一部分。\n",
    "如果星系数量较少，说明目标星系可能位于一个较为稀疏的区域。\n",
    "小结：\n",
    "通过上述步骤，我们可以：\n",
    "\n",
    "计算星系的距离，从而将其位置映射到三维空间中。\n",
    "评估星系周围的密度，通过查询邻近星系的分布来定量评估星系所处的环境。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72bfcd1e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}