坐标系统

一. 地心惯性坐标系 (ECI)

惯性坐标系是在空间固定不动或做匀速直线运动的坐标系，可用于描述天体和卫星的运行位置和状态 (地球上的点在该坐标系下一直处于变化中)。根据坐标系原点位置的不同，可分为日心惯性坐标系(原点在太阳质心)和地心惯性坐标系(原点在地球质心)(ECI)。

二. 地心地固坐标系 (ECEI)

地固坐标系是固定在地球上与地球一起旋转的坐标系，可用于描述接收机的地理位置。根据坐标系原点位置的不同，可分为参心坐标系(原点与参考椭球中心重合)和地心坐标系(原点与地球质心重合)，后者又称为地心地固坐标系(ECEI)。

地心/参心坐标系均包括大地坐标系(BLH)和空间直角坐标系(XYZ)两种形式。

以下公式待检核正确性

直角坐标系转为大地坐标系：

$$\begin{cases} B = arctan(\frac{Z(N+H)}{\sqrt{X^2+Y^2}(N(1-e^2)+H)})\\ L = arctan(\frac{Y}{X})\\ H = Z/sinB-N(1-e^2)\\ N = a/\sqrt{1-e^2sin^2B} \end{cases}$$

function [B, L, H] = XYZtoBLH(X, Y, Z)
    a = 6378137;
    e2 = 0.0066943799013;
    
    % L　in [0°, 360°]
    L = atand( abs(Y/X) );
    if Y > 0
        if X > 0
        else
            L = 180 - L;
        end
    else
        if X > 0
            L = 360 - L;
        else
            L = 180 + L;
        end
    end
    % L　in [-180°, 180°]
    if L > 180
       L = L - 360;
    end
    
    B0 = atand( Z / sqrt(X^2 + Y^2) );
    while 1
        N = a / sqrt( 1 - e2 * sind(B0)*sind(B0) );
        H = Z / sind(B0) - N * (1 - e2);
        B = atand( Z * (N + H) / ( sqrt(X^2 + Y^2) * (N * (1-e2) + H) ) );
        if abs(B-B0) < 1e-4, break; end
        B0 = B;
    end
end

大地坐标系转为直角坐标系：

$$\begin{cases} X = (N+H)cosBcosL\\ Y = (N+H)cosBsinL\\ Z = [N(1-e^2)+H]sinB\\ N = a/\sqrt{1-e^2sin^2B} \end{cases}$$

% BLH->XYZ
function [X, Y, Z] = BLHtoXYZ(B, L, H)
    a = 6378137;
    f = 1 / 298.257223563;
    e2 = 2 * f - f * f;
    N = a / sqrt(1 - e2 * sind(B) * sind(B));
    
    X = (N + H) * cosd(B) * cosd(L);
    Y = (N + H) * cosd(B) * sind(L);
    Z = (N * (1 - e2) + H) * sind(B);
end

三. 站心坐标系

站心坐标系是与测站固联的坐标系，可用于描述卫星相对于测站的高度角和方位角，如下图所示。

设GPS卫星在直角坐标系和站心坐标系下的坐标分别为($X^S$, $Y^S$, $Z^S$)和($N$, $E$, $U$)，接收机在直角坐标系和大地坐标系下的坐标分别为($X_R$, $Y_R$, $Z_R$)和($B$, $L$, $H$)，则卫星在站心坐标系下的坐标为：

$$\left[\begin{array}{c}N\\E\\U\end{array}\right] = \left[\begin{array}{ccc}-sinBcosL&-sinBsinL&cosB\\-sinL&cosL&0\\cosBcosL&cosBsinL&sinB\end{array}\right] \left[\begin{array}{c}X^S-X_R\\Y^S-Y_R\\Z^S-Z_R\end{array}\right]$$

计算高度角:

$$Elev = arctan(\frac{U}{\sqrt{N^2+E^2}})$$

计算方位角:

$$Azim = arctan(\frac{E}{N})$$

% 计算高度角和方位角
function [E, A] = get_elev(site_B, site_L, site_H, sate_B, sate_L, sate_H)
    % BLH->XYZ
    [site_x, site_y, site_z] = BLHtoXYZ(site_B, site_L, site_H);
    [sate_x, sate_y, sate_z] = BLHtoXYZ(sate_B, sate_L, sate_H);

    % 转化矩阵
    T = [-sind(site_B)*cosd(site_L), -sind(site_B)*sind(site_L), cosd(site_B)
              -sind(site_L)        ,       cosd(site_L)        ,      0
          cosd(site_B)*cosd(site_L),  cosd(site_B)*sind(site_L), sind(site_B)];

    % △XYZ->△NEU
    deta_xyz = ([sate_x, sate_y, sate_z] - [site_x, site_y, site_z])';
    deta_NEU = T * deta_xyz;

    % 高度角
    E = atand( deta_NEU(3) / sqrt( deta_NEU(1)^2 + deta_NEU(2)^2) );
    
    % 方位角
    A = atand( abs( deta_NEU(1) / deta_NEU(2) ) );
    if deta_NEU(2) > 0
        if deta_NEU(1) > 0
            A = 90 - A;
        else
            A = 90 + A;
        end
    else
        if deta_NEU(1) < 0
            A = 270 - A;
        else
            A = 270 + A;
        end
    end
end