Matlab实现月球天平动与日下点计算遇问题求助

你的代码存在三个关键问题，导致计算结果异常：

1. 旋转矩阵与向量的乘法顺序错误 #

planetEphemeris返回的sun_pos是3x1列向量，旋转矩阵作用于列向量时需要左乘（rotm * sun_vec），你之前使用的右乘（sun_vec * rotm）仅适用于行向量，会完全扭曲坐标变换结果。

2. 未显式指定欧拉角旋转顺序 #

moonLibration返回的[phi, theta, psi]对应**Z-Y-X（3-2-1）**旋转顺序，虽然eul2rotm默认也是该顺序，但显式指定可以避免版本差异或逻辑混淆，确保旋转矩阵生成正确。

3. 经度范围未做常规调整（可选） #

cart2sph返回的方位角（经度）范围是[-π, π]，而月球经度通常使用[0, 360°]表示，需要做简单的范围转换。

修正后的代码 #

mission_time = juliandate(2022, 1, 1);
% 获取ICRF坐标系下太阳相对月球的位置（3x1列向量）
sun_pos = planetEphemeris(mission_time, 'Moon', 'Sun');
% 获取月球天平动角：[经度天平动, 纬度天平动, 物理经度天平动]
moon_rot = moonLibration(mission_time);

% 生成月球固定系到ICRF的旋转矩阵，显式指定Z-Y-X旋转顺序
rotm = eul2rotm(moon_rot, 'ZYX');
% 转置得到ICRF到月球固定系的旋转矩阵（旋转矩阵的逆等于转置）
rotm = rotm';

% 归一化太阳方向向量（ICRF坐标系）
sun_vec = sun_pos / norm(sun_pos);
% 将向量转换到月球固定系（列向量左乘旋转矩阵）
sun_vec_moon = rotm * sun_vec;

% 转换为球坐标：方位角=经度，仰角=纬度
[ss_long_rad, ss_lat_rad, ~] = cart2sph(sun_vec_moon(1), sun_vec_moon(2), sun_vec_moon(3));

% 转换为角度并调整经度范围到0-360°
ss_long = rad2deg(ss_long_rad);
ss_lat = rad2deg(ss_lat_rad);
if ss_long < 0
    ss_long = ss_long + 360;
end

fprintf("Subsolar Lat: %2.2f°	Subsolar Long: %2.2f°\n", ss_lat, ss_long)

验证说明 #

修正后，2022年1月1日的日下点纬度会落在±1.5°左右（符合月球黄赤交角的范围），经度结果也会符合预期的天文观测值。

内容的提问来源于stack exchange，提问作者Big Dumb Rat