汇总三式，可得：
根据前面的推导过程，可以用个人擅长的编程语言进行编码。
其最关键的还是计算大地纬度B时的迭代过程，其余的计算都只是套公式。数值计算中的很多算法都是采用迭代趋近的方法来趋近一个最佳解。
在作者的个人测试中，最后的运行结果如下：
Web墨卡托投影是横轴墨卡托投影的特化版，要完全搞清楚Web墨卡托投影就必须得先搞清楚横轴墨卡托投影，不过横轴墨卡托投影实在太复杂了，但是我们可以定性地去理解。
它的计算过程大概可以这样理解：
- 在X方向上，为了保证投影到平面后经线和纬线仍然垂直，那么每条纬线都会按照赤道周长展开，也就是由于原点位于平面中心，那么可以算得X轴的取值范围：[-20037508.3427892,20037508.3427892]。经度与投影后X长是简单的线性关系；
- 在Y方向上，则需要借助于墨卡托投影公式。为了保证投影的结果是正方形，那么就把Y轴的取值范围也取值成[-20037508.3427892,20037508.3427892]之间。这样做没什么道理，纯粹是为了希望投影的结果是正方形，便于切片。最后，通过墨卡托投影公式进行反算，得到的经纬度范围就是[-85.05112877980659，85.05112877980659]。也就是这种投影方式，大于这个范围是失效的。
参考Cesium的具体实现如下：
#include <iostream>
//#include <eigen3/Eigen/Eigen>
//#include <osgEarth/GeoData>
using namespace std;
const double epsilon = 0.000000000000001;
const double pi = 3.14159265358979323846;
const double d2r = pi / 180;
const double r2d = 180 / pi;
const double a = 6378137.0; //椭球长半轴
const double f_inverse = 298.257223563; //扁率倒数
const double b = a - a / f_inverse;
//const double b = 6356752.314245; //椭球短半轴
const double e = sqrt(a * a - b * b) / a;
//墨卡托范围[-PI, PI]->大地纬度范围[-PI/2, PI/2]
static double mercatorAngleToGeodeticLatitude(double mercatorAngle)
{
return pi / 2.0 - (2.0 * atan(exp(-mercatorAngle)));
//return 2.0 * atan(exp(mercatorAngle)) - pi / 2.0;
}
//Web墨卡托投影所支持的最大纬度（北和南）
static double maximumLatitude = mercatorAngleToGeodeticLatitude(pi);
//大地纬度范围[-PI/2, PI/2]->墨卡托范围[-PI, PI]
static double geodeticLatitudeToMercatorAngle(double latitude)
{
// Clamp the latitude coordinate to the valid Mercator bounds.
if (latitude > maximumLatitude)
{
latitude = maximumLatitude;
}
else if (latitude < -maximumLatitude)
{
latitude = -maximumLatitude;
}
double sinLatitude = sin(latitude);
return 0.5 * log((1.0 sinLatitude) / (1.0 - sinLatitude));
}
void Blh2Wmc(double &x, double &y, double &z)
{
x = x * d2r * a;
y = geodeticLatitudeToMercatorAngle(y * d2r) * a;
}
void Wmc2Blh(double &x, double &y, double &z)
{
//var oneOverEarthSemimajorAxis = this._oneOverSemimajorAxis;
x = x / a * r2d;
y = mercatorAngleToGeodeticLatitude(y / a) * r2d;
}
int main()
{
double x = 113.6;
double y = 38.8;
double z = 100;
printf("%.10lf\n", maximumLatitude * r2d);
printf("原大地经纬度坐标：%.10lf\t%.10lf\t%.10lf\n", x, y, z);
Blh2Wmc(x, y, z);
printf("Web墨卡托坐标：%.10lf\t%.10lf\t%.10lf\n", x, y, z);
Wmc2Blh(x, y, z);
printf("转回大地经纬度坐标：%.10lf\t%.10lf\t%.10lf\n", x, y, z);
}
最终运行的结果：
令选取的站心点为P，其大地经纬度坐标为
