用python编写精确计算月球南极附近某位置的太阳高度角和太阳方位角，供输入位置的坐标系为Moon2000_spole，单位是m，可能需要坐标转换，且计算某时间段的太阳高度和方位角，而不是某时刻，结果进行可视化并输出CSV文件，我这里有naif0012.tls、pck00010.tpc、de405.bsp、de421.bsp、de432s.bsp、moon_080317.tf、moon_pa_de421_1900-2050.bpc

<think> 我们面临的任务是计算月球南极附近某位置在特定时间段内的太阳高度角和方位角，并将结果可视化和输出为CSV文件。我们需要使用SPICE工具包（通过SpiceyPy）来处理天体位置和姿态的计算。以下是详细的步骤： 1. **初始化SPICE内核**：加载提供的所有内核文件（SPK、PCK、LSK、TF、BPC）。 2. **坐标转换**：将输入位置从Moon2000_spole（月心固定坐标系，单位为米）转换为J2000坐标系（日心惯性坐标系），以便进行天体位置计算。注意：Moon2000_spole是月固坐标系，需要转换为惯性坐标系来计算相对于太阳的位置。 3. **时间处理**：将给定的时间段转换为ET（Ephemeris Time）时间。 4. **计算太阳位置**：在给定的时间点，计算太阳在月心J2000坐标系中的位置。 5. **计算观测点到太阳的向量**：在月心J2000坐标系中，计算观测点指向太阳的向量。 6. **将太阳向量转换到局部水平坐标系（East-North-Up）**：这样我们可以计算高度角和方位角。 7. **计算高度角和方位角**：利用转换后的局部水平坐标系中的向量计算太阳高度角和方位角。 8. **循环处理时间段内的每个时间点**：生成时间序列，对每个时间点执行上述计算。 9. **可视化结果**：绘制太阳高度角和方位角随时间变化的图表。 10. **输出CSV文件**：将时间、高度角、方位角等数据写入CSV文件。 ### 步骤详细说明 #### 1. 初始化SPICE内核 我们需要加载以下内核文件： - naif0012.tls（LSK，时间系统） - pck00010.tpc（PCK，行星常数） - moon_080317.tf（月球姿态模型） - de405.bsp, de421.bsp, de432s.bsp（星历数据） - moon_pa_de421_1900-2050.bpc（月球极轴方向） 注意：加载顺序很重要，后加载的内核可以覆盖先加载的内核。通常，先加载LSK和PCK，然后加载SPK和BPC。 #### 2. 坐标转换 输入位置是Moon2000_spole坐标系（月心固定坐标系）下的直角坐标（单位：米）。我们需要将这个位置转换为月心J2000坐标系下的位置，因为太阳的位置是在J2000坐标系下计算的。 具体步骤： - 使用`spiceypy.pxfrm2`或`spiceypy.sxform`函数获取从Moon2000_spole到J2000坐标系的旋转矩阵。 - 注意：Moon2000_spole是月固坐标系，而J2000是惯性坐标系，因此转换需要时间。 - 转换公式：设观测点在月固坐标系中的位置为`pos_moon_fixed`，则月心J2000坐标系中的位置为：`pos_j2000 = spiceypy.mxv(rotation_matrix, pos_moon_fixed)` 但是，我们实际上需要的是在月心J2000坐标系中，太阳相对于月球的位置？注意：我们要计算的是从观测点指向太阳的向量。因此，我们需要： - 计算太阳在月心J2000坐标系中的位置（向量）：`sun_pos_j2000` - 计算观测点在月心J2000坐标系中的位置：`obs_pos_j2000` - 那么，观测点指向太阳的向量为：`sun_vector_j2000 = sun_pos_j2000 - obs_pos_j2000` 然而，由于月球中心到太阳的距离远远大于观测点到月球中心的距离，我们可以忽略观测点相对于月球中心的位置，即近似认为观测点位于月球中心。但是，对于精确计算，特别是当需要计算太阳在观测点天空中的位置时，我们需要考虑视差（虽然很小，但为了精确，我们计算完整的向量）。 实际上，我们应该： - 计算太阳在月心J2000坐标系中的位置（向量），即从月球中心指向太阳的向量：`sun_to_moon_vector` - 计算观测点在月心J2000坐标系中的位置（向量）：`obs_pos_j2000` - 则观测点指向太阳的向量为：`sun_vector = sun_to_moon_vector - obs_pos_j2000` #### 3. 时间处理 输入时间段（起始时间和结束时间），将其转换为ET时间。输入时间可以是UTC字符串，然后使用`spiceypy.utc2et`转换为ET。 #### 4. 计算太阳位置 使用`spiceypy.spkpos`计算太阳相对于月球中心的位置（在J2000坐标系中）。 #### 5. 计算观测点位置在月心J2000坐标系中的表示 在给定了时间后，我们可以获取从月固坐标系（Moon2000_spole）到J2000的旋转矩阵，然后将观测点的月固坐标（单位：米）转换为J2000坐标系（单位：公里，注意单位转换）。 #### 6. 将太阳向量转换到局部水平坐标系 局部水平坐标系（East, North, Up）以观测点为中心。我们需要一个旋转矩阵，将J2000坐标系中的向量转换到该局部坐标系。 - 首先，我们需要定义局部坐标系的三个轴： - UP：观测点的月固坐标系中的位置向量（归一化）。 - EAST：在月固坐标系中，我们可以定义东方向（通常通过月球的自转轴和参考子午线确定）。但我们可以利用SPICE函数`spiceypy.surfpt`或`spiceypy.latrec`的反函数来得到经纬度，然后构建局部坐标系。 实际上，SPICE提供了函数`spiceypy.pgrrec`（行星图形经纬度到直角坐标）和`spiceypy.recpgr`（直角坐标到行星图形经纬度）。我们可以先得到观测点的经纬度，然后利用经纬度构建旋转矩阵（从月固坐标系到局部水平坐标系），然后再转换到J2000坐标系？或者直接在J2000坐标系中构建局部坐标系？ 更直接的方法是：在月固坐标系中，我们可以很容易地构建局部坐标系（East, North, Up），然后利用从月固到J2000的旋转矩阵，将整个局部坐标系转换到J2000坐标系中，从而得到J2000坐标系中的局部坐标系基向量。然后，我们就可以将太阳向量投影到这个局部坐标系中。 但是，有一个更简单的方法：SPICE提供了函数`spiceypy.srfrec`（将表面点从笛卡尔坐标转换为行星参考坐标系中的位置），但这不直接用于我们的目的。我们可以使用以下步骤： - 计算观测点在月固坐标系中的位置（已知），然后计算其经纬度（使用`spiceypy.reclat`或`spiceypy.recpgr`）。 - 然后，我们可以构建一个旋转矩阵，将月固坐标系中的向量转换到局部坐标系（East, North, Up）。这个旋转矩阵是固定的（不依赖于时间？注意：月固坐标系是固定的，所以局部坐标系的基在月固坐标系中是固定的）。然后，我们再将这个局部坐标系中的向量通过从月固到J2000的旋转矩阵转换到J2000坐标系中。 实际上，我们需要的是在J2000坐标系中的局部坐标系基向量。设： - 在月固坐标系中，局部基向量为： - UP：观测点位置归一化 - EAST：由月球自转轴和UP向量确定，具体为：`EAST = cross(Moon_pole, UP)`，然后归一化；再调整方向，使得指向东方。 - NORTH：`cross(UP, EAST)`归一化，即指向北。 然后，将这三个基向量（在月固坐标系中）通过旋转矩阵转换到J2000坐标系中，得到J2000坐标系中的局部坐标系基向量。 最后，将太阳向量（在J2000坐标系中）点乘这三个基向量，得到在局部坐标系中的三个分量（E, N, U）。然后，高度角为正切于水平面之上的角度，方位角从北顺时针（或逆时针）测量。 #### 7. 计算高度角和方位角 在局部水平坐标系中，我们有太阳向量的三个分量：`[e, n, u]`（东、北、上）。 - 高度角（elevation）= `arcsin(u / r)`，其中`r`是太阳向量的模（即到太阳的距离，但这里我们只需要方向向量，所以可以归一化后使用）。 - 方位角（azimuth）= `arctan2(e, n)`（注意：这里北为0度，东为90度）。如果`n`为0，则根据`e`的正负确定方位角（东或西）。 注意：`arctan2(e, n)`给出的是从北方向顺时针（在某些定义中）或逆时针（通常是逆时针）到东方向的角度。但通常方位角定义为从北方向顺时针旋转的角度（即北为0度，东为90度，南为180度，西为270度）。所以我们可以直接使用`azimuth = math.atan2(e, n)`，然后转换为度数（0到360度）。 #### 8. 循环计算时间序列 设定时间步长（比如1分钟、1小时等，根据需求），在时间区间内循环计算每个时间点的太阳高度角和方位角。 #### 9. 可视化 使用matplotlib绘制两个图：高度角随时间变化，方位角随时间变化（或者在一个极坐标图中绘制方位角-高度角）。 #### 10. 输出CSV 将时间（UTC字符串或ET）、高度角（度）、方位角（度）写入CSV文件。 ### 代码结构 1. 加载内核 2. 定义输入参数： - 观测点坐标（月固坐标系，单位：米） - 时间范围（UTC字符串）和时间步长 3. 将时间范围转换为ET，生成时间序列 4. 对于每个时间点： a. 计算太阳在月心J2000坐标系中的位置（向量） b. 计算观测点在月心J2000坐标系中的位置（向量） c. 计算观测点到太阳的向量（在J2000坐标系中） d. 归一化该向量（得到单位向量） e. 计算在月固坐标系中观测点的经纬度（用于构建局部坐标系） f. 构建月固坐标系中的局部坐标系基向量（东、北、上） g. 将基向量转换到J2000坐标系（使用当前时刻的旋转矩阵） h. 将太阳单位向量投影到局部坐标系（得到e, n, u分量） i. 计算高度角（elevation）和方位角（azimuth） 5. 存储结果 6. 绘制图表 7. 输出CSV文件 ### 注意点 - 单位：SPICE默认使用公里，所以输入位置（米）需要转换为公里（除以1000）。 - 时间：使用UTC字符串，通过`spiceypy.utc2et`转换为ET。 - 旋转矩阵：使用`spiceypy.pxform`或`spiceypy.sxform`获取从月固坐标系（Moon2000_spole）到J2000的旋转矩阵。注意月固坐标系的名称在SPICE中通常是“MOON_PA”或“MOON_ME”（但我们有moon_080317.tf内核，它定义了Moon2000_spole？我们需要确认坐标系名称）。在加载内核后，可以使用`spiceypy.bodvrd`或`spiceypy.furnsh`加载tf文件，然后坐标系名称可以从tf文件中获取。通常，月固坐标系的标准名称是“MOON_ME”或“MOON_PA”，但具体名称需要查看moon_080317.tf文件。如果不知道，可以通过`spiceypy.bodvrd`查询。这里假设为“MOON_PA”。 根据提供的文件，我们有一个moon_pa_de421_1900-2050.bpc文件，所以月固坐标系可能是“MOON_PA”。因此，在转换时，我们指定： - 源坐标系：'MOON_PA'（或'IAU_MOON'？注意，moon_080317.tf可能定义了'IAU_MOON'） - 目标坐标系：'J2000' 我们需要确认坐标系名称。通常，使用SPICE函数`spiceypy.furnsh`加载所有内核后，可以通过`spiceypy.bodvrd`获取信息。但为了简单，假设坐标系名称为'MOON_PA'。 ### 代码实现 以下为代码框架： ```python import spiceypy as spice import numpy as np import matplotlib.pyplot as plt import csv from datetime import datetime, timedelta # Step 1: 加载内核 def load_kernels(kernel_files): for kernel in kernel_files: spice.furnsh(kernel) # 内核文件列表 kernel_files = [ 'naif0012.tls', 'pck00010.tpc', 'de405.bsp', 'de421.bsp', 'de432s.bsp', 'moon_080317.tf', 'moon_pa_de421_1900-2050.bpc' ] load_kernels(kernel_files) # Step 2: 输入参数 # 观测点位置（Moon2000_spole坐标系，单位米） obs_pos_moon_fixed_m = np.array([x, y, z]) # 例如月球南极附近的点 # 转换为公里（SPICE使用公里） obs_pos_moon_fixed = obs_pos_moon_fixed_m / 1000.0 # 时间范围（UTC） utc_start = "2024-01-01T00:00:00" utc_end = "2024-01-02T00:00:00" time_step_seconds = 3600 # 1小时步长 # 将UTC时间转换为ET et_start = spice.utc2et(utc_start) et_end = spice.utc2et(utc_end) # 生成时间序列 times_et = np.arange(et_start, et_end, time_step_seconds) num_times = len(times_et) # 存储结果 times_utc = [] elevations = [] azimuths = [] # 定义坐标系名称 moon_fixed_frame = 'MOON_PA' # 根据内核文件，月固坐标系名称 inertial_frame = 'J2000' # 月球中心ID和太阳中心ID moon_id = 301 sun_id = 10 # 循环计算 for i, et in enumerate(times_et): # a. 计算太阳在月心J2000坐标系中的位置（向量，单位公里） sun_pos, _ = spice.spkpos('SUN', et, inertial_frame, 'NONE', 'MOON') # 注意：spkpos返回的位置是从观测者（这里是月球中心）指向目标（太阳）的向量 # b. 计算观测点在月心J2000坐标系中的位置（向量） # 首先，获取从月固坐标系到J2000的旋转矩阵 rotation_matrix = spice.pxform(moon_fixed_frame, inertial_frame, et) # 将观测点位置从月固坐标系转换到J2000坐标系 obs_pos_j2000 = spice.mxv(rotation_matrix, obs_pos_moon_fixed) # c. 计算观测点到太阳的向量（在J2000坐标系中） # 太阳相对于月球中心的位置向量是sun_pos（从月球中心指向太阳） # 观测点相对于月球中心的位置向量是obs_pos_j2000 # 因此，从观测点指向太阳的向量 = (太阳相对于月球中心的向量) - (观测点相对于月球中心的向量) sun_vector = sun_pos - obs_pos_j2000 # d. 归一化 sun_vector_unit = sun_vector / np.linalg.norm(sun_vector) # e. 计算在月固坐标系中观测点的经纬度（用于构建局部坐标系） # 先将观测点月固坐标转换为球坐标（经度、纬度、半径） # 使用函数spice.recpgr（将直角坐标转换为行星图形坐标） radius, lon, lat = spice.reclat(obs_pos_moon_fixed) # reclat返回的是笛卡尔坐标对应的球坐标：半径、经度（弧度）、纬度（弧度） # 注意：月球的经度范围是0-360度，东经为正？通常，经度是东经0到360度。 # f. 构建月固坐标系中的局部坐标系基向量（东、北、上） # 首先，构建一个旋转矩阵，该矩阵可以将月固坐标系中的向量转换到局部坐标系（East, North, Up） # 我们可以使用球坐标的经纬度来定义旋转矩阵 # 旋转矩阵的列向量分别是东、北、上方向在月固坐标系中的分量 # 我们可以使用函数spiceypy.latrec来反向计算（但这里我们直接构造） # 或者使用以下方法： # 上方向：在月固坐标系中，归一化的观测点位置向量 up_vec = obs_pos_moon_fixed / np.linalg.norm(obs_pos_moon_fixed) # 北方向：在月固坐标系中，与月球自转轴和上方向有关 # 月球自转轴在月固坐标系中是Z轴（0,0,1）？我们需要月球北极的方向 # 根据内核，月球北极在月固坐标系中是固定的（通常是Z轴正方向） moon_pole = np.array([0.0, 0.0, 1.0]) # 东方向：与月球自转轴和上方向垂直的方向，即 moon_pole × up_vec，然后归一化 east_vec = np.cross(moon_pole, up_vec) east_vec = east_vec / np.linalg.norm(east_vec) # 北方向：up_vec × east_vec，归一化（注意叉乘顺序，确保右手系） north_vec = np.cross(up_vec, east_vec) north_vec = north_vec / np.linalg.norm(north_vec) # 注意：这里构建的局部坐标系基向量在月固坐标系中 # g. 将基向量转换到J2000坐标系 # 使用前面得到的旋转矩阵 east_vec_j2000 = spice.mxv(rotation_matrix, east_vec) north_vec_j2000 = spice.mxv(rotation_matrix, north_vec) up_vec_j2000 = spice.mxv(rotation_matrix, up_vec) # h. 将太阳单位向量投影到局部坐标系（得到e, n, u分量） e = np.dot(sun_vector_unit, east_vec_j2000) n = np.dot(sun_vector_unit, north_vec_j2000) u = np.dot(sun_vector_unit, up_vec_j2000) # i. 计算高度角（elevation）和方位角（azimuth） # 高度角：从水平面向上为正，向下为负 elevation = np.arcsin(u) # 单位为弧度 # 方位角：从北方向顺时针旋转的角度（0到2π），使用atan2(e, n)得到从北逆时针转到东（但通常方位角是顺时针从北开始） # 注意：数学中的atan2(y,x)得到的是从x轴逆时针旋转到向量的角度 # 这里，我们定义x为北方向，y为东方向，则atan2(e, n)得到的是从北方向逆时针旋转的角度（即东偏北角度） # 但通常方位角是顺时针从北开始，所以我们可以用：azimuth = (2*np.pi - azimuth_math + np.pi/2) % (2*np.pi) ？这样很复杂。 # 实际上，我们可以直接计算：azimuth = np.arctan2(e, n) 得到从北方向逆时针旋转的角度（数学标准），然后转换为顺时针： # 但我们希望方位角范围是0到360度，北为0度，东为90度，南为180度，西为270度。 azimuth = np.arctan2(e, n) # 这是从北方向逆时针旋转的角度（弧度） # 如果azimuth为负，转换为0到2π if azimuth < 0: azimuth += 2 * np.pi # 转换为度数 elevation_deg = np.degrees(elevation) azimuth_deg = np.degrees(azimuth) # 记录结果 times_utc.append(spice.et2utc(et, 'C', 0)) # 转换为UTC字符串 elevations.append(elevation_deg) azimuths.append(azimuth_deg) # 卸载内核 spice.kclear() # 输出CSV文件 with open('solar_angles.csv', 'w', newline='') as csvfile: writer = csv.writer(csvfile) writer.writerow(['Time (UTC)', 'Elevation (deg)', 'Azimuth (deg)']) for i in range(len(times_utc)): writer.writerow([times_utc[i], elevations[i], azimuths[i]]) # 可视化 plt.figure(figsize=(12, 6)) plt.subplot(2, 1, 1) plt.plot(times_utc, elevations, 'b-') plt.ylabel('Elevation (deg)') plt.title('Solar Elevation and Azimuth') plt.xticks(rotation=45) plt.grid(True) plt.subplot(2, 1, 2) plt.plot(times_utc, azimuths, 'r-') plt.ylabel('Azimuth (deg)') plt.xlabel('Time (UTC)') plt.xticks(rotation=45) plt.grid(True) plt.tight_layout() plt.savefig('solar_angles.png') plt.show() # 另外，可以绘制极坐标图 plt.figure(figsize=(8, 8)) ax = plt.subplot(111, projection='polar') ax.scatter(np.radians(azimuths), 90 - np.array(elevations), c=range(len(azimuths)), cmap='viridis', alpha=0.75) ax.set_theta_zero_location('N') # 将0度（北）放在顶部 ax.set_theta_direction(-1) # 顺时针方向 ax.set_rmax(90) ax.set_rticks([0, 30, 60, 90]) # 高度角范围，这里表示天顶角（90-高度角） ax.set_rlabel_position(45) # 移动半径标签的位置 ax.set_title('Solar Position (Azimuth vs Zenith Angle)') plt.savefig('solar_polar.png') plt.show() ``` ### 注意事项 1. 月球自转轴：在月固坐标系中，我们假设Z轴是自转轴（北极向上）。这是标准设置。 2. 坐标转换的旋转矩阵：`spice.pxform`函数需要源坐标系、目标坐标系和时间。 3. 太阳位置计算：`spkpos`函数返回的是从观测者（这里是月球中心）指向目标（太阳）的向量。 4. 观测点位置：在月固坐标系中给出，需要转换为公里。 5. 局部坐标系的构建：我们假设月球自转轴是固定的（在月固坐标系中为Z轴）。构建东、北、上向量时，使用叉乘得到垂直向量。 6. 方位角的定义：我们定义方位角从北方向顺时针旋转（0度北，90度东，180度南，270度西）。在计算时，我们使用`np.arctan2(e, n)`得到从北方向逆时针旋转的角度，然后转换为0~360度（没有进行顺时针转换，因为`arctan2(e,n)`的逆时针旋转角度与顺时针旋转的角度在数值上互补？）。实际上，我们直接使用逆时针旋转的角度值，然后按顺时针方向显示？不对，我们只需将方位角的值按逆时针显示在图表上（但通常人们习惯顺时针）。 - 实际上，方位角定义为从北方向顺时针旋转的角度。而我们计算时，使用`arctan2(e,n)`得到的是从北方向逆时针旋转到东方向的角度。因此，如果我们希望得到顺时针角度，可以用`2*pi - azimuth_math`。但注意，当`azimuth_math`为0时，顺时针也是0；当`azimuth_math`为90度（东）时，顺时针角度应为90度，而`2*pi - 90`=270度，这不对。 - 因此，正确的转换是：顺时针角度 = `2*pi - azimuth_math`，然后对2*pi取模。或者直接使用：`azimuth_clockwise = (2*np.pi - azimuth_math) % (2*np.pi)`。这样，当`azimuth_math`为0时，顺时针0度；当`azimuth_math`为90度（π/2）时，顺时针角度=2π-π/2=3π/2（即270度）。但这不是我们想要的（东应该是90度，而不是270度）。 - 实际上，我们重新审视局部坐标系：东方向（x）、北方向（y）、上方向（z）构成右手系。方位角是从北方向顺时针旋转到投影向量的角度。即，在水平面上，将向量投影到水平面，然后从北方向顺时针旋转到该投影向量的角度。 - 投影向量的分量是(e, n)，那么从北方向顺时针旋转到该向量的角度α满足：tan(α)=e/n（注意：这里不是标准的笛卡尔坐标系，而是数学中的极坐标）。实际上，α = atan2(e, n) 然后取负号再调整？或者直接使用α = atan2(-e, n) 然后调整到0~2π？ - 在标准气象学中，方位角定义为从北方向顺时针旋转的角度。我们可以使用： azimuth_deg = (np.degrees(np.arctan2(e, n)) + 360) % 360 但是，这样得到的是从北方向逆时针旋转的角度（因为arctan2(e,n)就是x轴为北，y轴为东的坐标系中的角度）。而顺时针旋转的角度应该是360减去这个值？不对，因为： - 当向量指向东时，e>0, n=0，则arctan2(e,n)=90度，而我们希望方位角为90度（东）。 - 当向量指向南时，e=0, n<0，arctan2(0,-1)=180度，对应方位角180度（南）。 - 当向量指向西时，e<0, n=0，arctan2(-1,0)= -90度（或270度），对应方位角270度（西）。 - 当向量指向北时，e=0, n>0，arctan2(0,1)=0度，对应0度。 因此，直接使用`azimuth = np.arctan2(e, n)`并将结果转换为度数，然后取模到0~360度，就是我们想要的方位角（从北方向逆时针旋转的角度等于顺时针旋转的角度？）。但实际上，这个方位角的值与气象学中的定义一致（0度北，90度东，180度南，270度西）。所以，我们不需要转换。 - 结论：直接使用`azimuth = (np.arctan2(e, n) + 2*np.pi) % (2*np.pi)`得到弧度，然后转换为角度，就是气象学中定义的方位角。 7. 性能：如果时间点很多，循环计算可能较慢。可以尝试向量化，但SPICE函数通常一次只能计算一个时间点。 8. 内核文件：确保所有内核文件在正确的路径下，或者指定完整路径。 ### 补充说明 以上代码和步骤提供了完整的解决方案。但请注意，天体位置计算非常复杂，且依赖于精确的内核文件。在实际应用中，可能需要根据具体需求调整参数和步骤。此外，月球南极附近的位置可能需要特殊处理（例如，当高度角为负时，太阳不可见）。 最后，代码示例中使用了假设的观测点坐标，请替换为实际值。 ### 相关问题 1. 如何获取月球南极在Moon2000_spole坐标系中的坐标？ 2. 如何处理内核文件加载失败的情况？ 3. 计算中如何考虑光行时修正（light time correction）？ 4. 如何提高计算效率（例如，减少时间步长或并行计算）？ 5. 如何将结果可视化为3D图形？ </think> ### 精确计算月球太阳高度角与方位角的Python实现 以下是完整的解决方案，使用SPICE工具包进行天体位置计算和坐标转换： ```python import spiceypy as spice import numpy as np import matplotlib.pyplot as plt import pandas as pd from datetime import datetime, timedelta # 初始化SPICE内核 def load_kernels(): kernels = [ 'naif0012.tls', # 时间系统内核 'pck00010.tpc', # 行星常数内核 'de405.bsp', # 行星历表 'de421.bsp', 'de432s.bsp', 'moon_080317.tf', # 月球坐标系定义 'moon_pa_de421_1900-2050.bpc' # 月球极轴方向 ] for kernel in kernels: spice.furnsh(kernel) # 坐标转换：Moon2000_spole (m) 到 J2000 (km) def convert_coords(pos_m, et): """将Moon2000_spole坐标系(米)转换为J2000坐标系(公里)""" # 单位转换：米 → 公里 pos_km = np.array(pos_m) / 1000.0 # 获取转换矩阵 transform = spice.pxform('MOON_PA', 'J2000', et) # 应用坐标转换 return spice.mxv(transform, pos_km) # 计算太阳向量 def get_sun_vector(et): """获取从观测点到太阳的向量(J2000坐标系)""" sun_pos, _ = spice.spkpos('SUN', et, 'J2000', 'LT+S', 'MOON') return np.array(sun_pos) # 计算高度角和方位角 def calculate_angles(sun_vec, obs_pos, et): """计算太阳高度角和方位角""" # 归一化太阳向量 sun_dir = sun_vec / np.linalg.norm(sun_vec) # 获取月球坐标系转换矩阵 transform = spice.pxform('J2000', 'MOON_PA', et) obs_pos_moon = spice.mxv(transform, obs_pos) # 计算经度纬度 _, lon, lat = spice.reclat(obs_pos_moon) # 构建局部坐标系 (ENU: East, North, Up) up = obs_pos_moon / np.linalg.norm(obs_pos_moon) east = np.array([-np.sin(lon), np.cos(lon), 0.0]) north = np.cross(up, east) # 转换到J2000坐标系 east_j2000 = spice.mxv(transform.T, east) north_j2000 = spice.mxv(transform.T, north) up_j2000 = spice.mxv(transform.T, up) # 计算投影分量 e = np.dot(sun_dir, east_j2000) n = np.dot(sun_dir, north_j2000) u = np.dot(sun_dir, up_j2000) # 计算高度角和方位角 elevation = np.arcsin(u) # 高度角(弧度) azimuth = np.arctan2(e, n) # 方位角(弧度) # 转换为角度并调整范围 elevation_deg = np.degrees(elevation) azimuth_deg = np.degrees(azimuth) % 360 return elevation_deg, azimuth_deg # 主计算函数 def calculate_solar_angles(position, start_time, end_time, time_step): """ 计算指定时间段内的太阳高度角和方位角 参数： position: Moon2000_spole坐标系下的位置 (m) [x, y, z] start_time: 起始时间 (UTC字符串，如 '2024-06-01T00:00:00') end_time: 结束时间 time_step: 时间步长 (秒) """ # 加载内核 load_kernels() # 时间转换 et_start = spice.utc2et(start_time) et_end = spice.utc2et(end_time) time_steps = int((et_end - et_start) / time_step) # 结果存储 results = [] for i in range(time_steps + 1): et = et_start + i * time_step utc_time = spice.et2utc(et, 'C', 0) # 坐标转换 obs_pos_j2000 = convert_coords(position, et) # 计算太阳向量 sun_vec = get_sun_vector(et) sun_vec_adj = sun_vec - obs_pos_j2000 # 计算角度 elevation, azimuth = calculate_angles(sun_vec_adj, obs_pos_j2000, et) results.append({ 'UTC Time': utc_time, 'Elevation (deg)': elevation, 'Azimuth (deg)': azimuth, 'ET': et }) # 转换为DataFrame df = pd.DataFrame(results) # 卸载内核 spice.kclear() return df # 可视化结果 def visualize_results(df, title): """可视化太阳高度角和方位角""" plt.figure(figsize=(15, 10)) # 高度角图 plt.subplot(2, 1, 1) plt.plot(df['UTC Time'], df['Elevation (deg)'], 'b-') plt.title(f'{title} - Solar Elevation Angle') plt.ylabel('Elevation (degrees)') plt.grid(True) plt.xticks(rotation=45) # 方位角图 plt.subplot(2, 1, 2) plt.scatter(df['UTC Time'], df['Azimuth (deg)'], c=df['Azimuth (deg)'], cmap='hsv', s=10) plt.title(f'{title} - Solar Azimuth Angle') plt.ylabel('Azimuth (degrees)') plt.colorbar(label='Azimuth (deg)') plt.grid(True) plt.xticks(rotation=45) plt.tight_layout() plt.savefig('lunar_solar_angles.png', dpi=300) plt.show() # 极坐标图 plt.figure(figsize=(10, 8)) ax = plt.subplot(111, projection='polar') sc = ax.scatter(np.radians(df['Azimuth (deg)']), df['Elevation (deg)'], c=df['ET'], cmap='viridis', alpha=0.7) ax.set_theta_zero_location('N') ax.set_theta_direction(-1) ax.set_title('Polar Representation of Solar Position') plt.colorbar(sc, label='Ephemeris Time (s)') plt.savefig('solar_position_polar.png', dpi=300) plt.show() # 示例使用 if __name__ == "__main__": # 月球南极附近位置 (Moon2000_spole坐标系，单位：m) # 示例位置：月球南极Shackleton陨石坑边缘 position = [0.0, 0.0, -1737000] # 实际位置需根据任务需求调整 # 时间参数 (UTC) start_time = '2024-06-01T00:00:00' end_time = '2024-06-03T00:00:00' time_step = 3600 # 1小时步长 (秒) # 执行计算 df = calculate_solar_angles(position, start_time, end_time, time_step) # 输出CSV df.to_csv('lunar_solar_angles.csv', index=False) print("结果已保存至 lunar_solar_angles.csv") # 可视化 visualize_results(df, "Lunar South Pole Solar Angles") ``` ### 关键步骤说明 1. **内核加载**： - 加载所有必需的SPICE内核文件，包括时间系统、行星常数和星历数据 - `naif0012.tls`定义时间系统，`pck00010.tpc`提供行星常数，`moon_080317.tf`定义月球坐标系[^1] 2. **坐标转换**： - 将输入的Moon2000_spole坐标(米)转换为J2000坐标系(公里) - 使用`spice.pxform`获取坐标系转换矩阵，`spice.mxv`应用转换[^2] 3. **天体位置计算**： - 使用`spice.spkpos`计算太阳相对于月球的精确位置 - 应用光行时校正(`LT+S`选项)提高精度[^3] 4. **角度计算**： - 通过局部ENU坐标系(东-北-天顶)计算方位角和高度角 - 使用球面三角学公式计算太阳方向矢量在局部坐标系中的投影[^4] 5. **时间处理**： - UTC时间与星历时间(ET)相互转换 - 支持任意时间步长的计算 6. **可视化**： - 高度角时间序列图 - 方位角彩色散点图 - 极坐标太阳位置图 ### 使用说明 1. **输入参数**： - `position`: 月球南极附近的三维坐标(米)，Moon2000_spole坐标系 - `start_time`/`end_time`: UTC时间范围(格式: 'YYYY-MM-DDTHH:MM:SS') - `time_step`: 计算步长(秒) 2. **输出结果**： - CSV文件包含时间、高度角、方位角 - 可视化图表：时间序列图+极坐标图 3. **依赖库**： ```bash pip install spiceypy numpy pandas matplotlib ``` ### 注意事项 1. 坐标系定义依赖`moon_080317.tf`文件，确保使用正确的坐标系名称('MOON_PA') 2. 月球南极位置需根据实际任务指定，示例位置仅供参考 3. 对于长期计算，建议添加光行时校正和相对论效应 4. 结果精度受内核文件影响，`de432s.bsp`提供高精度月球定位