2024MathModeling/Problem3.py at master · hengxt/2024MathModeling · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
## 问题3  在建立脉冲星光子到达航天器（卫星等）与太阳系质心之间的精确转换时延模型时，需要考虑脉冲星自行的影响以及几个关键的时延因素：几何传播时延、Shapiro时延、引力红移时延和狭义相对论的动钟变慢效应。在考虑脉冲星自行以及上述时延因素下，请建立脉冲星光子到达航天器与太阳系质心的精确转换时延模型。若光子到达探测器的时刻对应MJD为58119.1651507519，根据问题1中卫星的位置、速度以及附件4提供的脉冲星位置参考历元、自行参数信息，计算脉冲星光子到达航天器与太阳系质心间的时延（附件3提供了太阳系天体位置信息的DE系列历表）。

"""
#### 问题3  在建立脉冲星光子到达航天器（卫星等）与太阳系质心之间的精确转换时延模型时，需要考虑脉冲星自行的影响以及几个关键的时延因素：几何传播时延、Shapiro时延、引力红移时延和狭义相对论的动钟变慢效应。在考虑脉冲星自行以及上述时延因素下，请建立脉冲星光子到达航天器与太阳系质心的精确转换时延模型。若光子到达探测器的时刻对应MJD为58119.1651507519，根据问题1中卫星的位置、速度以及附件4提供的脉冲星位置参考历元、自行参数信息，计算脉冲星光子到达航天器与太阳系质心间的时延（附件3提供了太阳系天体位置信息的DE系列历表）。

#### 解决步骤总结：
1. 载入太阳系天体位置数据:
 - 使用jplephem库中的SPK类打开并加载太阳系天体位置的DE系列历表文件。
2. 获取观测时刻的地球位置:
- 根据给定的儒略日（JD），调用compute_earth_position函数计算地球相对于太阳系质心的位置向量。
3. 更新卫星位置:
- 将卫星在地心惯性参考系中的位置向量与地球位置向量相加，得到卫星在太阳系质心坐标系中的位置向量。
4. 更新脉冲星位置:
- 考虑脉冲星的自行参数，使用update_pulsar_position函数更新脉冲星的赤经和赤纬。
5. 计算指向脉冲星的单位矢量:
- 利用更新后的脉冲星位置信息，调用compute_unit_vector函数计算指向脉冲星方向的单位矢量。
6. 计算几何传播时延（Roemer时延）:
- 通过compute_roemer_delay函数基于卫星位置和指向脉冲星的方向来计算光子传播的时间延迟。
7. 计算Shapiro时延:
- 首先获取太阳的位置向量，然后利用compute_shapiro_delay函数计算由于太阳引力导致的光传播路径弯曲所产生的额外时间延迟。
8. 计算引力红移时延:
- 使用compute_gravitational_redshift函数根据地球的位置计算由太阳引力场引起的频率变化所对应的时延。
9. 计算动钟变慢效应:
- 依据卫星的速度，调用compute_relativistic_effect函数估算狭义相对论效应下的时间膨胀。
10. 汇总所有时延:
- 将上述各部分时延相加以获得总的时延，并输出各个分项及时延总和的结果。

"""
from jplephem.spk import SPK
import numpy as np

# 常量
c = 299792458  # 光速，单位：米/秒
G = 6.67430e-11  # 万有引力常数，单位：m^3 kg^(-1) s^(-2)
M_sun = 1.989e30  # 太阳质量，单位：kg

# 载入太阳系天体位置历表 (DE系列)
kernel = SPK.open("附件/附件3-de200.bsp")


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


def compute_satellite_position(r_sat_gcrs, r_earth):
    """计算卫星在太阳系质心坐标系中的位置"""
    r_sat_bcrs = r_sat_gcrs + r_earth
    return r_sat_bcrs


def compute_unit_vector(ra, dec):
    """计算指向脉冲星的单位矢量"""
    ra_rad = np.deg2rad(ra)
    dec_rad = np.deg2rad(dec)
    n_hat = np.array([np.cos(dec_rad) * np.cos(ra_rad),
                      np.cos(dec_rad) * np.sin(ra_rad),
                      np.sin(dec_rad)])
    return n_hat


def compute_roemer_delay(r_sat_bcrs, n_hat, c):
    """计算几何传播时延（Roemer时延）"""
    delta_t_roemer = np.dot(r_sat_bcrs, n_hat) / c
    return delta_t_roemer


def compute_shapiro_delay(r_sat_bcrs, r_sun, G, M_sun, c):
    """计算Shapiro时延"""
    r_psr = np.linalg.norm(r_sat_bcrs - r_sun)
    delta_t_shapiro = (2 * G * M_sun / c ** 3) * np.log((r_psr + np.linalg.norm(r_sat_bcrs)) / r_psr)
    return delta_t_shapiro


def compute_gravitational_redshift(r_earth, G, M_sun, c):
    """计算引力红移时延"""
    U_sun = G * M_sun / np.linalg.norm(r_earth)
    delta_t_grav_redshift = U_sun / c ** 2
    return delta_t_grav_redshift


def compute_relativistic_effect(v_sat, c):
    """计算动钟变慢效应"""
    delta_t_relativistic = 0.5 * (v_sat / c) ** 2
    return delta_t_relativistic


def update_pulsar_position(ra, dec, dt_ra, dt_dec, jd, pch_pm):
    """更新脉冲星位置"""
    dt_ra_rad = np.deg2rad(dt_ra / 1000 / 3600)  # 转换为弧度/年
    dt_dec_rad = np.deg2rad(dt_dec / 1000 / 3600)  # 转换为弧度/年
    delta_t = (jd - 2400000.5) - pch_pm  # 单位：儒略日
    ra_updated = ra + dt_ra_rad * (delta_t / 365.25)  # 考虑自行影响
    dec_updated = dec + dt_dec_rad * (delta_t / 365.25)  # 考虑自行影响
    return ra_updated, dec_updated


# 观测时刻的儒略日
jd = 58119.1651507519 + 2400000.5  # 对应MJD为58119.1651507519

# 卫星在地心惯性参考系中的位置和速度 (通过轨道根数计算)
r_sat_gcrs = np.array([1274913.41509, -1848851.96499, 6507262.65528])  # 单位：米
v_sat_gcrs = np.array([-6210.79517, 3746.12731, 2276.33088])  # 单位：米/秒

# 脉冲星位置和自行参数
ra_pulsar = 83.63308333  # 脉冲星赤经，单位：度
dec_pulsar = 22.0145  # 脉冲星赤纬，单位：度
dt_ra = -14.7  # 赤经自行，单位：毫角秒/年
dt_dec = 2.0  # 赤纬自行，单位：毫角秒/年
pch_pm = 57715.000000295  # par文件的参考历元

# 计算地球位置
r_earth = compute_earth_position(jd)

# 更新脉冲星位置
ra_pulsar, dec_pulsar = update_pulsar_position(ra_pulsar, dec_pulsar, dt_ra, dt_dec, jd, pch_pm)

# 计算卫星在太阳系质心坐标系中的位置
r_sat_bcrs = compute_satellite_position(r_sat_gcrs, r_earth)

# 计算指向脉冲星的单位矢量
n_hat = compute_unit_vector(ra_pulsar, dec_pulsar)

# 计算几何传播时延（Roemer时延）
delta_t_roemer = compute_roemer_delay(r_sat_bcrs, n_hat, c)

# 计算太阳位置
r_sun = kernel[0, 10].compute(jd) * 1000  # 太阳相对于太阳系质心的位置，单位：米

# 计算Shapiro时延
delta_t_shapiro = compute_shapiro_delay(r_sat_bcrs, r_sun, G, M_sun, c)

# 计算引力红移时延
delta_t_grav_redshift = compute_gravitational_redshift(r_earth, G, M_sun, c)

# 计算动钟变慢效应
v_sat = np.linalg.norm(v_sat_gcrs)
delta_t_relativistic = compute_relativistic_effect(v_sat, c)

# 总时延
total_delay = delta_t_roemer + delta_t_shapiro + delta_t_grav_redshift + delta_t_relativistic

# 输出结果
print("几何传播时延: ", delta_t_roemer)
print("Shapiro时延: ", delta_t_shapiro)
print("引力红移时延: ", delta_t_grav_redshift)
print("动钟变慢效应: ", delta_t_relativistic)
print("脉冲星自行: ", ra_pulsar, dec_pulsar)
print("=================================")
print(f"总时延：{total_delay} 秒")