J2000toWGS84.rar

J2000坐标系转换WGS84坐标系代码。代码使用MATLAB编写。% 计算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: 读取地球的地极坐标（xpyp）。
[nonRynon] = rotate_matrix（ -deg2radian（00-0.15） ）;
[Rxnonnon] = rotate_matrix（ -deg2radian（000.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（232621.448） - deg2radian（0046.8150）*T0 - deg2radian（000.00059）*T0^2 + deg2radian（000.001813）*T0^3;   %弧度角
DE_Blocks = DE405_load（JD_time）;           %读取DE405数据包
DE_angle =Planet_Position_J2000_Earth（JD_timeDE_Blocks.Blocks14）;   %获取对应时间的黄经章动和交角章动
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（0067310.54841） + deg2radian（87660008640184.812866）*Tu + ...
            deg2radian（000.093104）*Tu^2 - deg2radian（000.62e-5）*Tu^3）*15 + ...
            delta_fai *cos（e_M + delta_e） + 2.61036020374891  ;   %--  计算格林尼治恒星时[弧度];
    
[nonnonB] = rotate_matrix（seta_G）;  %-调用旋转矩阵
%% -------------------- 计算章动矩阵 C ----------------------------
[Rx_enonnon] = rotate_matrix（-e_M-delta_e）;
[nonnonRz_fai] = rotate_matrix（-delta_fai）;
[Rx_eMnonnon] = rotate_matrix（e_M）;
C = Rx_e * Rz_fai * Rx_eM;
%% -------------------- 计算岁差矩阵 D --------------------------
ebucailong_p = deg2radian（002306.2181）*T0 + deg2radian（000.30188）*T0^2 + deg2radian（000.017998）*T0^3;
seta_p       = deg2radian（002004.3109）*T0 - deg2radian（000.42665）*T0^2 + deg2radian（000.041833）*T0^3;
Z_p          = deg2radian（002306.2181）*T0 + deg2radian（001.09468）*T0^2 + deg2radian（000.018203）*T0^3;
%-
[nonnonRz_Z] = rotate_matrix（-Z_p）;
[nonRy_setanon] = rotate_matrix（seta_p）;
[nonnonRz_ep] = rotate_matrix（-ebucailong_p）;
D = Rz_Z*Ry_seta*Rz_ep;
%% ------------------- 计算旋转矩阵 W  -----------------------
 W = A*B*C*D;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

    I.A....      2443  2016-02-18 21:41  J2000toWGS84.m

----------- ---------  ---------- -----  ----

                 2443                    1