ECI坐标转换为ECEF坐标matlab代码

     ECI坐标是指地球固定坐标系，ECEF坐标是指地心地固坐标系。将ECI坐标转换为ECEF坐标需要考虑地球自转的影响。

    具体的转换步骤如下：

1. 获取当前时间的UT1(世界时1)时间和地球自转角。

2. 根据UT1时间计算格林尼治平恒星时(GMST)。

3. 将ECI坐标系下的三维坐标转换为四元数。

4. 根据GMST计算出相应的转换矩阵。

5. 利用转换矩阵将EC坐标系下的坐标转换为ECEF坐标系下的坐标。

6. 最后根据地球半径和所得到的ECEF坐标计算出相应的地理坐标.

7. 需要注意的是，该转换过程需要考虑一些细节问题，如时间系统的选择、转换矩阵的构建等。

      详细理论方法可见：Orbital Mechanics with Numerit, http://www.cdeagle.com/omnum/pdf/csystems.pdf

       提供详细的matlab代码供大家使用：

function [r_ECEF v_ECEF a_ECEF] = ECItoECEF(JD,r_ECI,v_ECI,a_ECI)
%Enforce JD to be [N x 1]
JD = JD(:);
%Calculate the Greenwich Apparent Sideral Time (THETA)
%See http://www.cdeagle.com/omnum/pdf/csystems.pdf equation 27
THETA = JD2GAST(JD);
%Average inertial rotation rate of the earth radians per second
omega_e = 7.29211585275553e-005;
%Assemble the transformation matricies to go from ECI to ECEF

r_ECEF = squeeze(MultiDimMatrixMultiply(T3D(THETA),r_ECI));
v_ECEF = squeeze(MultiDimMatrixMultiply(T3D(THETA),v_ECI) + ...
    MultiDimMatrixMultiply(Tdot3D(THETA,omega_e),r_ECI));
a_ECEF = squeeze(MultiDimMatrixMultiply(T3D(THETA),a_ECI) + ...
    2.*MultiDimMatrixMultiply(Tdot3D(THETA,omega_e),v_ECI) + ...
    MultiDimMatrixMultiply(Tddot3D(THETA,omega_e),r_ECI));
%---------------------------------------------------------------------%'

function C = MultiDimMatrixMultiply(A,B)
C = sum(bsxfun(@times,A,repmat(permute(B',[3 2 1]),[size(A,2) 1 1])),2);
function T = T3D(THETA)
T = zeros([3 3 length(THETA)]);
T(1,1,:) = cosd(THETA);
T(1,2,:) = sind(THETA);
T(2,1,:) = -T(1,2,:);
T(2,2,:) = T(1,1,:);
T(3,3,:) = ones(size(THETA));
function Tdot = Tdot3D(THETA,omega_e)
Tdot = zeros([3 3 length(THETA)]);
Tdot(1,1,:) = -omega_e.*sind(THETA);
Tdot(1,2,:) = omega_e.*cosd(THETA);
Tdot(2,1,:) = -Tdot(1,2,:);
Tdot(2,2,:) = Tdot(1,1,:);
function Tddot = Tddot3D(THETA,omega_e)
Tddot = zeros([3 3 length(THETA)]);
Tddot(1,1,:) = -(omega_e.^2).*cosd(THETA);
Tddot(1,2,:) = -(omega_e.^2).*sind(THETA);
Tddot(2,1,:) = -Tddot(1,2,:);
Tddot(2,2,:) =  Tddot(1,1,:);
%------------------------------------------------------------------%
function GAST = JD2GAST(JD)
%THETAm is the mean siderial time in degrees
THETAm = JD2GMST(JD);
%Compute the number of centuries since J2000
T = (JD - 2451545.0)./36525;
%Mean obliquity of the ecliptic (EPSILONm)
% see http://www.cdeagle.com/ccnum/pdf/demogast.pdf equation 3
% also see Vallado, Fundamentals of Astrodynamics and Applications, second edition.
%pg. 214 EQ 3-53
EPSILONm = 23.439291-0.0130111.*T - 1.64E-07.*(T.^2) + 5.04E-07.*(T.^3);
%Nutations in obliquity and longitude (degrees)
% see http://www.cdeagle.com/ccnum/pdf/demogast.pdf equation 4
L = 280.4665 + 36000.7698.*T;
dL = 218.3165 + 481267.8813.*T;
OMEGA = 125.04452 - 1934.136261.*T;
%Calculate nutations using the following two equations:
% see http://www.cdeagle.com/ccnum/pdf/demogast.pdf equation 5
dPSI = -17.20.*sind(OMEGA) - 1.32.*sind(2.*L) - .23.*sind(2.*dL) ...
    + .21.*sind(2.*OMEGA);
dEPSILON = 9.20.*cosd(OMEGA) + .57.*cosd(2.*L) + .10.*cosd(2.*dL) - ...
    .09.*cosd(2.*OMEGA);
%Convert the units from arc-seconds to degrees
dPSI = dPSI.*(1/3600);
dEPSILON = dEPSILON.*(1/3600);
%(GAST) Greenwhich apparent sidereal time expression in degrees
% see http://www.cdeagle.com/ccnum/pdf/demogast.pdf equation 1
GAST = mod(THETAm + dPSI.*cosd(EPSILONm+dEPSILON),360);
%------------------------------------------------------------------%
function GMST = JD2GMST(JD)
%Find the Julian Date of the previous midnight, JD0
JD0 = NaN(size(JD));
JDmin = floor(JD)-.5;
JDmax = floor(JD)+.5;
JD0(JD > JDmin) = JDmin(JD > JDmin);
JD0(JD > JDmax) = JDmax(JD > JDmax);
H = (JD-JD0).*24;       %Time in hours past previous midnight
D = JD - 2451545.0;     %Compute the number of days since J2000
D0 = JD0 - 2451545.0;   %Compute the number of days since J2000
T = D./36525;           %Compute the number of centuries since J2000
%Calculate GMST in hours (0h to 24h) ... then convert to degrees
GMST = mod(6.697374558 + 0.06570982441908.*D0  + 1.00273790935.*H + ...
    0.000026.*(T.^2),24).*15;
%----------------------------End Function----------------------------------

函数具体输入输出参量：
JD                         [1 x N]                           儒略日向量
r_ECI                    [3 x N]                           ECI坐标系中的位置向量
v_ECI                   [3 x N]                           ECI坐标系中的速度向量
a_ECI                   [3 x N]                           ECI坐标系中的加速度向量
输出：
r_ECEF                 [3 x N]                          ECEF坐标系中的位置向量
v_ECEF                [3 x N]                          ECEF坐标系中的速度向量
a_ECEF                [3 x N]                          ECEF坐标系中的加速度向量