functionW=J2000_2_WGS84_matrix(JD_time)%计算J2000坐标系转换WGS84坐标系的旋转矩阵%目前暂未考虑地球极移矩阵%输入参数：%JD_time：儒略日期%输出参数：%W：WGS84坐标系与J2000.0坐...

function W = J2000_2_WGS84_matrix(JD_time)

% 计算J2000坐标系转换WGS84坐标系的旋转矩阵

% 目前暂未考虑地球极移矩阵

% 输入参数：

% JD_time： 儒略日期

% 输出参数：

% W： WGS84坐标系与J2000.0坐标系的转换矩阵，

%% -------------------- 计算极移矩阵 A --------------------------

%-- load: 读取地球的地极坐标(xp,yp)。

[non,Ry,non] = rotate_matrix( -deg2radian(0,0,-0.15) );

[Rx,non,non] = rotate_matrix( -deg2radian(0,0,0.3) );

%-- A = Ry(-xp)*Rx(-yp);

A = Ry*Rx;

%% ------------------- 计算自转矩阵 B----------------------------

T0 = ((JD_time - 2400000.5) - 51544.5)/36525.0; %-- 从J2000起到t时刻的儒略世纪数

Tu = (JD_time - 2451545.0)/36525.0;

%--

e_M = deg2radian(23,26,21.448) - deg2radian(0,0,46.8150)*T0 - deg2radian(0,0,0.00059)*T0^2 + deg2radian(0,0,0.001813)*T0^3; %弧度角

DE_Blocks = DE405_load(JD_time); %读取DE405数据包

DE_angle =Planet_Position_J2000_Earth(JD_time,DE_Blocks.Blocks,14); %获取对应时间的黄经章动和交角章动

delta_fai = DE_angle(1);

delta_e = DE_angle(2);

%--

% seta_G = ( 100.46061837 + 36000.770053608*Tu + 0.000387933*Tu^2 - Tu^3/38710000 )/180*pi ...

% + delta_fai *cos(e_M + delta_e) + 0.178401471734655; %-- 计算格林尼治恒星时[弧度];

seta_G = ( deg2radian(0,0,67310.54841) + deg2radian(876600,0,8640184.812866)*Tu + ...

deg2radian(0,0,0.093104)*Tu^2 - deg2radian(0,0,0.62e-5)*Tu^3)*15 + ...

delta_fai *cos(e_M + delta_e) + 2.61036020374891 ; %-- 计算格林尼治恒星时[弧度];

[non,non,B] = rotate_matrix(seta_G); %-调用旋转矩阵

%% -------------------- 计算章动矩阵 C ----------------------------

[Rx_e,non,non] = rotate_matrix(-e_M-delta_e);

[non,non,Rz_fai] = rotate_matrix(-delta_fai);

[Rx_eM,non,non] = rotate_matrix(e_M);

C = Rx_e * Rz_fai * Rx_eM;

%% -------------------- 计算岁差矩阵 D --------------------------

ebucailong_p = deg2radian(0,0,2306.2181)*T0 + deg2radian(0,0,0.30188)*T0^2 + deg2radian(0,0,0.017998)*T0^3;

seta_p = deg2radian(0,0,2004.3109)*T0 - deg2radian(0,0,0.42665)*T0^2 + deg2radian(0,0,0.041833)*T0^3;

Z_p = deg2radian(0,0,2306.2181)*T0 + deg2radian(0,0,1.09468)*T0^2 + deg2radian(0,0,0.018203)*T0^3;

%-

[non,non,Rz_Z] = rotate_matrix(-Z_p);

[non,Ry_seta,non] = rotate_matrix(seta_p);

[non,non,Rz_ep] = rotate_matrix(-ebucailong_p);

D = Rz_Z*Ry_seta*Rz_ep;

%% ------------------- 计算旋转矩阵 W -----------------------

W = A*B*C*D;

展开