php cookbook有一个计算两个经纬点的例子，代码就照着抄了一遍，多次检查，多次运行的结果都是略小于预期值的一半，楞是费了一番，才发现是经纬度参数的顺序反了，在英文里，是纬度在前，经度在后，不似耳熟的东经多少多少，北纬多少多少。既然到了这了，也学着梳理一下计算经纬点的数学步骤。对这个，园里softfair的博文很给力，图文对照，https://www.cnblogs.com/softfair/p/distance_of_two_latitude_and_longitude_points.html，是以此文仅是加深一下学习的印象。
（以下论述纯自我理解，必有错有误）地球的自转轴与公转平面约23.4的夹角，沿赤道线自转，受离心力和月球引力的影响，地球的赤道面略宽，两极略扁的不规则椭圆球体，表面地形的高低起伏，有山高海平，故而这里的经纬点计算仅是假想地球的规则的球体，忽略地形的近似结果，误差究竟有多少呢，有的说几米，有的说0.5%（若参照的是1km，0.5%就是5米，还对得上；若是参照实际距离的0.5%，比如超百公里，就偏差过大），海拔相差百米内，距离千公里内的两点，就想当然预计是10米内了，日常测距中可忽略不计。
一，HALF-VERSINE，就是sin角的一半,三角函数是初中就接触了的，也都荒废了，有个别的公式比较直观，一看之下就能想起来，比如这个sinθ2+cosθ2=1,还有这个sin(θ/2)^2=(1-cosθ)/2,后者转换一下就是sin(θ/2)=sqrt((1-cosθ)/2),这就是需要用到的HALF-VERSINE了。引用博文的图片很直观，就直接复制粘贴了。
二，平面圆通过半径和圆心角求弦
如图是设半径为1，弦长为L，等腰三角形作高垂直于弦，根据三角函数的定义，直观得到sin(θ/2)=1/2L，得L=2sin(θ/2)
三，球面上的两个点，坐标P1(lat1,lon1),p2(lat2,lon2)，通过互换经纬度坐标，很容易得到另两个对称的点P3(lat1,lon2),P4(lat2,lon1),这四个就像刀切苹果一样，连起来就是去除了外弧线部分的梯形
经纬度有现实的意义的，两个点的度差，就是规则的角度差，经纬皆是，像经度是整个圆分布，就是0180度，球体的两面就是东经和西经；纬度由赤道开始沿向两极，就是半圆，是090度，视觉上下就是北纬和南纬。如引图，就理解到：角AOC是A点与C点的纬度差dlat。角EOF是经度E点和经度F点的差dlon，EF所在的线是赤道(equator)，纬度是0。EF的距离是EF=2sin(dlon/2)。
既然角度有了，就根据三角函数的定义，得到弦AC的长度，AC=2sin(dlat/2)，弦BD和AC是等腰梯形对称的两个腰，故AC=BD（补充：OA和OC都是球体的半径哦，故三角形AOC是等腰的，由上文二的申明直观得到）。
对AD的弦长，引文详细，通过A点作垂线于赤道上的半径，由三角函数推得，这里简单的想当然一下，因赤道到两极的纬线是渐次缩小为零，直观上等价于呈余弦的曲率变化，姑A点上的半径就是赤道的半径rcos(lat1),设半径为1，就是cos(lat1)，B点就是cos(lat2)。如此弦长AD= 2sin(dlon/2)cos(lat1)，类推出CB=2sin(dlon/2)*cos(lat2)。
四，如图，现在需求得AD的弦长，以进一步获得AD的弧长。
AH^2 = AC^2 - CH^2 = AC^2 - (CB-AD)^2/4
HB= AD+CH = AD + (CB-AD)/2 = (CB+AD)/2
AB^2 = AH^2 + HB^2 = AC^2 - (CB-AD)^2/4 + (CB+AD)^2/4 = AC^2 + CB*AD
根据前面求得的AC、AD和CB的长度，代入公式得到：
AB^2 = 4(sin^2(dlat/2) + 4cos(lat1)cos(lat2)sin^2(dlon/2))
再需求得代表AB长度的角度AOB：
若规则球体，圆点到任意球面点的距离都是半径，前边已设r=1，则有
tan(<AOC) = AC/OC = (AB/2)/sqrt(1-(AB/2)^2)
∠AOC = arctan((AB/2)/sqrt(1-(AB/2)^2))
弧长AB = EARTH_RADIUS * ∠AOC = 2r * arctan((AB/2)/sqrt(1-(AB/2)^2))
到了这，问题看似解决了，只是和预期的表达式相差过大，是哪里没转换到吗？一时不了，还是动手运行一下看看：
多番调试纠错，排除了和cookbook上冲突的代码，得到的结果是预期的2倍，奇怪，哪里错了呢，再发现引文在得出最终结论时，用了这个表达：
弧长AB = EARTH_RADIUS * ∠AOC = 2r * arctan((AB/2)/sqrt(1-(AB/2)^2))
实际
c = 2 * arctan(sqrt(a)/sqrt(1-a))，
这一项说明c就是∠AOC全角的度数了，而不是半角，结论那根本不用乘二，怕是博主埋的瓜，修改至此得证引文的论述是成立的。
