074-RTKlib中ECEF坐标转换为大地坐标问题

RTKlib中ECEF坐标转换为大地坐标的方法跟平时使用的不同，特此记录。
常用的大地坐标转换为ECEF坐标的公式有：

$$\{\begin{array}{l}x=(R_n+h)cosLcos\lambda \\ y=(R_n+h)cosLsin\lambda \\ z=[R_n(1-e^2)+h]sinL\end{array}\text{\hspace{1em}(1)}$$

各符号具体含义不再赘述。

整理（1）式中第三式可得：

$$(R_n+h)sinL=z+R_n{e}^{2}sinL\text{\hspace{1em}(2)}$$

为直观表现，参考下图：

为与图中字母对应，改写公式（2）为：

$$(N+h)sin\varphi =z+N{e}^{2}sin\varphi \text{\hspace{1em}(3)}$$

RTKlib程序：

extern void ecef2pos(const double *r, double *pos)
{
    double e2=FE_WGS84*(2.0-FE_WGS84),r2=dot(r,r,2),z,zk,v=RE_WGS84,sinp; 
    for (z=r[2],zk=0.0;fabs(z-zk)>=1E-4;) {
        zk=z;
        sinp = z / sqrt(r2 + z*z);
        v=RE_WGS84/sqrt(1.0-e2*sinp*sinp); 
        z=r[2]+v*e2*sinp;
    }
    pos[0]=r2>1E-12?atan(z/sqrt(r2)):(r[2]>0.0?PI/2.0:-PI/2.0);
    pos[1]=r2>1E-12?atan2(r[1],r[0]):0.0;
    pos[2]=sqrt(r2+z*z)-v;
}

由于$sin\varphi =AC/(N+h)$，所以对于sinp = z / sqrt(r2 + z*z);，得到的为小于$\varphi$的角度的正弦值，相应的，程序中的卯酉圈曲率半径v，即公式中的$N$也是小于实际值。
循环中的最后一行z=r[2]+v*e2*sinp;为z向图中$AC$的逼近。注意：z越逼近$AC$，得到的纬度越准确。最终逼近$AC$，从而得到较为准确的$\varphi$，便可进行下面的解算。