前三篇介绍了坐标系和矩阵的数学知识，从本篇开始，我们试图运用这些知识来解决实际问题。
如上图，模拟了一个以球心为原点的固定坐标系，该坐标系有一个名称地心地固坐标系(ECEF)，对应我们之前介绍的坐标系
，而平面场景在我们生活中更为直观，上北下南，左东右西，对应上图中绿色的切平面，简称NEU坐标系，对应之前介绍的坐标系
。于是，给定一点
，我们需要计算一个矩阵
，实现两个坐标系的转换。
这里对应两个环节，(1)球心坐标系的单位换算， 从经纬度
到米单位的笛卡尔坐标
；(2)从ECEF到NEU，从全球坐标系
到本地坐标系
。
整体来看，默认初始时
，从F到M大概需要三个过程：(1)沿
逆时针旋转
，如上图，橙色对应的是
，红色对应的是
，方向均向内；(2)沿着新坐标系中的红轴逆时针旋转
；（3）沿新坐标系的
方向平移到绿色坐标系的原点。
前两个旋转矩阵对应的是：
这样，只要知道平移
对应的三个分量，可以轻松的得到最终的矩阵。然后这里就有一个问题，地球是椭球而不是圆球，如下图所示，维度是
，但准确的切平面对应的法线是
，前者是geocentric normal，后者是geodetic normal，通常情况下，我们都是指的
。
椭球下经纬度和笛卡尔坐标转换的问题暂时先不做处理，假设已知某点对应球心的位置
，椭球在X，Y，Z方向的半径分别为a, b，c。则，up方向也就是该点对应椭球切面的法线(geodetic normal)：
又因为地球对应的椭球体中，a==b。North轴指向北极，而East轴和Up轴，North轴正交，可得East轴必须垂直于
和
：
而
，因此，我们可以获取ENU坐标系三个轴的向量
，这样，对应的转换公式为：
这样，我们在ENU本地坐标系上的一点
，对应球心坐标系上的点
，满足：
如上，我们实现了ECEF和ENU之间的转化，下面，我们讲一下经纬度到ECEF之间的转换，该问题可以抽象为已知经纬度 高度
，这里的
对应ECEF坐标系下的
，求其在笛卡尔坐标系下对应的位置
，这块并不涉及到坐标系的概念，不感兴趣的可以略过。
不失一般性，我们认为地球的椭球为：
因为
是geodetic normal，所以，在点
处的法线为
，同样，由椭圆的属性可知，该点的法线（未归一化）同样可以是
，则可得如下关系：
将
三个未知变量放到等式左侧：
同时，因为
在椭球体上，带入椭球体的公式，可得：
此时，我们可以得出
三个值，也就是
点的未知，最终就可以得到绿色点的位置：
接下来，自然就是逆运算了，已知
，如何求
。
如上图，一旦我们知道
：
已知
，可得：
同样的思路，将
三个未知变量放到等式左侧：
这里，我们定义一个函数S：
这是求解一元四次方程，我觉得应该可以直接公式计算解析解吧？不过标准方式是采用的牛顿迭代法。初始值就是以geocentric和
连线与椭球体的焦点，此时
：
这样，我们可以计算函数S相对于
的导数，用牛顿迭代法不断逼近，找到满足自己要求的
值：
我个人对这个迭代对应的几何意义理解大概如下，迭代调整α，不断的缩小alpha*n_s的距离：
Iteration 1：初始值，得到
Iteration 2：以
的geodetic normal方向进行scale，得到某个未知，连接r，得到
Iteration 3：同理，不断逼近，得到满足容限的近似解
本文到此结束，主要介绍了球心坐标系转换的相关内容，因为涉及到椭球而变得有些复杂，但在坐标系转换这方面并不复杂。下一篇介绍运动学中常用的Denavit-Hartenberg Algorithm。
参考资料：Cesium
3D Engine Design for Virtual Globes第二章
