Problem2.py

import numpy as np
from jplephem.spk import SPK

# 常量
c = 299792458  # 光速，单位：米/秒


# 切换坐标系，以太阳为中心， Barycentric Coordinate System, BCS
def pos_to_BCS(pos_GCRS, TT):
    kernel = SPK.open("de421.bsp")
    jd = TT + 2400000.5  # 转换为儒略日
    # 计算地球（id=3）相对于太阳系质心（id=0）的位置
    r_earth = kernel[0, 3].compute(jd)  # 返回地心在BCRS中的位置，单位：km
    r_earth = r_earth * 1000  # 转换为米
    # 计算卫星在太阳系质心坐标系中的位置
    r_sat_BCS = pos_GCRS + r_earth  # 卫星在BCRS中的位置，单位：米
    # 卫星位置转换
    return r_sat_BCS


# 1、计算Roemer时延
def roemer_func(pos_BCS, k):
    return np.dot(pos_BCS, k) / c


# 2、获取卫星在太阳系质心坐标系(BCS)中的位置
# 输入卫星在GCRS中的位置和速度 (单位：米) Geocentric Celestial Reference System，GCRS
p_sat_GCS = np.array([1274913.41509, -1848851.96499, 6507262.65528])  # 位置 米
v_sat_GCS = np.array([-6210.79517, 3746.12731, 2276.33088])  # 速度 米/S
print("卫星在GCRS中的位置：", p_sat_GCS)
p_sat_BCS = pos_to_BCS(p_sat_GCS, 57062.0)
print("卫星在太阳系质心坐标系中的位置：", p_sat_BCS)

# 3、计算指向脉冲星的单位矢量
## 根据下面的参数计算出单位矢量
# 给定的赤经和赤纬
ra_deg = 83.63307631324167  # deg
dec_deg = 22.014493269173464  # deg
# 转换为弧度
ra_rad = np.radians(ra_deg)
dec_rad = np.radians(dec_deg)
# 计算单位矢量 k
k = np.array([
    np.cos(dec_rad) * np.cos(ra_rad),
    np.cos(dec_rad) * np.sin(ra_rad),
    np.sin(dec_rad)
])
# 归一化确保 k 是单位矢量
k /= np.linalg.norm(k)

# 4、计算传播时间差
delay_sat = roemer_func(p_sat_BCS, k)
print("传播时间差: {} s".format(delay_sat))

"""

# 载入太阳系天体位置历表 (DE系列)
kernel = SPK.open("de421.bsp")

# 输入光子到达卫星的时间对应的MJD（约化儒略日）
mjd = 57062.0
jd = mjd + 2400000.5  # 转换为儒略日

# 计算地球（id=3）相对于太阳系质心（id=0）的位置
r_earth = kernel[0, 3].compute(jd)  # 返回地心在BCRS中的位置，单位：km
r_earth = r_earth * 1000  # 转换为米

# 计算卫星在太阳系质心坐标系中的位置
r_sat_BCRS = p_sat_GCS + r_earth  # 卫星在BCRS中的位置，单位：米

# 假设脉冲星信号为平行光，定义脉冲星方向单位向量
# 这里以假设的方向向量为例（需要根据实际情况进行修改）
n_hat = np.array([0.0, 0.0, 1.0])  # 假设脉冲星沿Z轴方向辐射

# 计算几何传播时延（Roemer时延）
delta_t_roemer = np.dot(r_sat_BCRS, n_hat) / c

# 输出结果
print(f"几何传播时延 (Roemer时延)：{delta_t_roemer:.10f} 秒")

"""