{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f6c70e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#task8.对于作业5中的星系，依据fits图像画出其一维的表面亮度轮廓，并计算其星等\n",
    "#• 注意流量单位\n",
    "#• 一维轮廓：\n",
    "        #✓椭圆测光：椭圆是等面亮度轮廓，半径沿着长轴方向\n",
    "        #✓圆测光: 圆孔径内的平均表面亮度\n",
    "#• 星等：\n",
    "     #✓积分（求和）到某个等亮度半径处\n",
    "     #✓根据面亮度轮廓模型，积分到无穷远处\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import fits\n",
    "from photutils.isophote import Ellipse, EllipseGeometry, build_ellipse_model\n",
    "\n",
    "# 加载FITS文件\n",
    "file_path = '149.74358_2.125144_sci.fits'  # 替换为你的FITS文件路径\n",
    "\n",
    "\n",
    "with fits.open(file_path) as hdul:\n",
    "    hdul.info()\n",
    "    image_data = hdul[1].data\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import fits\n",
    "from photutils.detection import find_peaks\n",
    "from photutils.isophote import EllipseGeometry, Ellipse\n",
    "\n",
    "# 1. 加载FITS文件并提取图像数据\n",
    "file_path = '149.74358_2.125144_sci.fits'  # 替换为你的FITS文件路径\n",
    "with fits.open(file_path) as hdul:\n",
    "    image_data = hdul[1].data\n",
    "\n",
    "# 检查 image_data 是否为空\n",
    "if image_data is None:\n",
    "    raise ValueError(\"图像数据未正确加载，请检查FITS文件。\")\n",
    "\n",
    "# 2. 提取亮点区域\n",
    "# 设定一个阈值来检测亮点\n",
    "threshold = np.percentile(image_data, 99)  # 选择99%亮度作为阈值（可调整）\n",
    "\n",
    "# 使用 photutils 的 find_peaks 函数来查找亮点\n",
    "peaks = find_peaks(image_data, threshold=threshold)\n",
    "\n",
    "# 3. 绘制原始图像和亮点区域\n",
    "plt.figure(figsize=(8, 8))\n",
    "plt.imshow(image_data, origin='lower', cmap='gray', vmin=np.min(image_data), vmax=np.max(image_data))\n",
    "plt.colorbar()\n",
    "plt.title(\"原始图像及亮点区域\")\n",
    "\n",
    "# 标出亮点区域\n",
    "plt.scatter(peaks['x_peak'], peaks['y_peak'], color='red', label='Peak locations', marker='x')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "# 4. 选择第一个亮点区域进行椭圆拟合\n",
    "# 假设第一个亮点是我们的目标\n",
    "x_center, y_center = peaks['x_peak'][0], peaks['y_peak'][0]\n",
    "\n",
    "# 定义椭圆拟合的初始参数\n",
    "sma = 10  # 半长轴（以像素为单位）\n",
    "eps = 0.1  # 离心率\n",
    "pa = 0.0  # 位置角（以弧度为单位）\n",
    "\n",
    "# 创建椭圆几何体\n",
    "geometry = EllipseGeometry(x0=x_center, y0=y_center, sma=sma, eps=eps, pa=pa)\n",
    "\n",
    "# 使用 Ellipse 进行拟合\n",
    "ellipse = Ellipse(image_data, geometry)\n",
    "isolated_fits = ellipse.fit_isophote(sma)\n",
    "\n",
    "# 5. 绘制拟合的椭圆轮廓\n",
    "plt.figure(figsize=(8, 8))\n",
    "plt.imshow(image_data, origin='lower', cmap='gray', vmin=np.min(image_data), vmax=np.max(image_data))\n",
    "plt.colorbar()\n",
    "plt.title(\"椭圆轮廓拟合\")\n",
    "\n",
    "# 获取拟合椭圆的坐标，使用 isolated_fits.sample.x 和 isolated_fits.sample.y\n",
    "sample = isolated_fits.sample\n",
    "\n",
    "# 从椭圆样本中提取 x 和 y 坐标\n",
    "x_coords = sample.x\n",
    "y_coords = sample.y\n",
    "\n",
    "# 绘制拟合的椭圆轮廓\n",
    "plt.plot(x_coords, y_coords, color='red')\n",
    "\n",
    "# 标出亮点位置\n",
    "plt.scatter(peaks['x_peak'], peaks['y_peak'], color='blue', label='Peak locations', marker='x')\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f31e267c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distance in pc: 2992264.64 +/- 110389.83\n"
     ]
    }
   ],
   "source": [
    "#作业9\n",
    "import numpy as np\n",
    "m = 25\n",
    "M = -3\n",
    "sigma_m = 0.1\n",
    "Rv = 3.1\n",
    "E_BV = 0.2\n",
    "sigma_E_BV = 0.05\n",
    "\n",
    "A_v = Rv * E_BV\n",
    "m_corrected = m - A_v\n",
    "\n",
    "distance_parsecs = 10 ** ((m_corrected - M + 5) / 5)\n",
    "\n",
    "sigma_Av = Rv * sigma_E_BV\n",
    "sigma_m_corrected = np.sqrt(sigma_m**2 + sigma_Av**2)\n",
    "sigma_distance_parsecs = (distance_parsecs / 5) * sigma_m_corrected\n",
    "\n",
    "print(f\"Distance in pc: {distance_parsecs:.2f} +/- {sigma_distance_parsecs:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a7e2b54",
   "metadata": {},
   "outputs": [],
   "source": [
    "#作业10\n",
    "from astropy.cosmology import FlatLambdaCDM\n",
    "cosmo = FlatLambdaCDM(H0=10, Om0=0.28)\n",
    "\n",
    "halpha = 6563\n",
    "halpha_z = 7000\n",
    "z = halpha_z / halpha - 1\n",
    "\n",
    "D_L = cosmo.luminosity_distance(z).value\n",
    "D_L_pc = D_L * 1e6  # 转换为pc\n",
    "flux = 0.6 * (6700 - 5500) * 1e-10  # 1e-10用于转换为erg/s/cm^2\n",
    "AB_mag=20\n",
    "M_r = AB_mag - 5 * (np.log10(D_L_pc) - 1) - 2.5 * np.log10(flux)  # 将D_L转换为pc\n",
    "M_r\n"
   ]
  }
 ],
 "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
}