C# 方位和距离计算（使用 Vincenty 的公式）结果为 NaN？

Question

按照本文计算经纬度和方位角和距离的点，当运行结果为 NaN 时？

以下文章：

给出了使用Vincenty公式计算椭球地球模型目的地点的方法：

将起点纬度'lat1'（在-90 到90 范围内）转换为弧度。 lat1 = lat1 * PI/180 将起点经度'lon1'（在-180 到180 范围内）转换为弧度。 lon1 = lon1 * PI/180 将方位 'brg'（在 0 到 360 范围内）转换为弧度。 brg = brg * PI/180 椭球体地球模型

注：下面的变量'flat'表示各种椭球体模型中使用的地球极地扁平化。对于常用的WGS-84，设flat = 298.257223563.

给定以米为单位的距离 s，以米为单位的半长轴 'a'，以米为单位的半短轴 'b' 和极平面 'flat'。使用 Vincenty 公式计算目标点。使用缩短的变量名称。

f = 1/flat
sb=sin(brg)
cb=cos(brg)
tu1=(1-f)*tan(lat1)
cu1=1/sqrt((1+tu1*tu1))
su1=tu1*cu1
s2=atan2(tu1, cb)
sa = cu1*sb
csa=1-sa*sa
us=csa*(a*a - b*b)/(b*b)
A=1+us/16384*(4096+us*(-768+us*(320-175*us)))
B = us/1024*(256+us*(-128+us*(74-47*us)))
s1=s/(b*A)
s1p = 2*PI

当条件为真时循环执行以下内容。

while (abs(s1-s1p) > 1e-12)
cs1m=cos(2*s2+s1)
ss1=sin(s1)
cs1=cos(s1)
ds1=B*ss1*(cs1m+B/4*(cs1*(-1+2*cs1m*cs1m)- B/6*cs1m*(-3+4*ss1*ss1)*(-3+4*cs1m*cs1m)))
s1p=s1
s1=s/(b*A)+ds1

循环结束后继续计算

t=su1*ss1-cu1*cs1*cb
lat2=atan2(su1*cs1+cu1*ss1*cb, (1-f)*sqrt(sa*sa + t*t))
l2=atan2(ss1*sb, cu1*cs1-su1*ss1*cb)
c=f/16*csa*(4+f*(4-3*csa))
l=l2-(1-c)*f*sa* (s1+c*ss1*(cs1m+c*cs1*(-1+2*cs1m*cs1m)))
d=atan2(sa, -t)
finalBrg=d+2*PI
backBrg=d+PI
lon2 = lon1+l;
Convert lat2, lon2, finalBrg and backBrg to degrees
lat2 = lat2 * 180/PI
lon2 = lon2 * 180/PI
finalBrg = finalBrg * 180/PI
backBrg = backBrg * 180/PI

如果 lon2 在 -180 到 180 的范围之外，加或减 360 以使其回到该范围内。如果 finalBrg 或 backBrg 超出 0 到 360 的范围，则加或减 360 以将它们带回到该范围内。注意：上面的变量'a'、'b'和'flat'有如下关系：

b = a - (a/flat)
flat = a / (a - b)

所以将这个公式移植到 C# 我最终得到：

        double Latitude = 50.390202;
        double Longitude = -3.9204310000000078;
        double Bearing = 225;

        double lat1 = Latitude * (Math.PI / 180.0);
        double lon1 = Longitude * (Math.PI / 180.0);
        double brg = Bearing * (Math.PI / 180.0);
        double s = 1000;

        double a = 6378137.0;
        double b = 6356752.314245;

        double f = 1 / 298.257223563;
        double sb = Math.Sin(brg);
        double cb = Math.Cos(brg);
        double tu1 = (1-f) * Math.Tan(lat1);
        double cu1 = 1 / Math.Sqrt((1+tu1*tu1));
        double su1 = tu1 * cu1;
        double s2 = Math.Atan2(tu1, cb);
        double sa = cu1 * sb;
        double csa = 1 - sa * sa;
        double us = csa * (a * a - b * b) / (b * b);
        double A = 1 + us / 16384 * (4096 + us * (- 768 + us * (320 - 175 * us)));
        double B = us / 1024 * (256 + us * (-128 + us * (74 - 47 * us)));
        double s1 = s / (b * A);
        double s1p = 2 * Math.PI;

        while (Math.Abs(s1-s1p) > 1e-12)
        {
            cs1m = Math.Cos(2 * s2 + s1);
            ss1 = Math.Sin(s1);
            cs1 = Math.Cos(s1);
            double ds1 = B * ss1 * (cs1m + B / 4 * (cs1 * (-1 + 2 * cs1m * cs1m) - B / 6 * cs1m * (-3 + 4 * ss1 * ss1) * (-3 + 4 * cs1m * cs1m)));
            s1p = s1;
            s1 = s / (b * A) + ds1;
        }

        double t = su1 * ss1 - cu1 * cs1 * cb;
        double lat2 = Math.Atan2(su1 * cs1 + cu1 * ss1 * cb, (1 - f) * Math.Sqrt(sa * sa + t * t));
        double l2 = Math.Atan2(ss1 * sb, cu1 * cs1 - su1 * ss1 * cb);
        double c = f / 16 * csa * (4 + f * (4 - 3 * csa));
        double l = l2 - (1 - c) * f * sa * (s1 + c * ss1 * (cs1m + c * cs1 * (-1 + 2 * cs1m * cs1m)));
        double d = Math.Atan2(sa, -t);
        double finalBrg = d + 2 * Math.PI;
        double backBrg = d + Math.PI;
        double lon2 = lon1 + l;

        lat2 = lat2 * 180 / Math.PI;
        lon2 = lon2 * 180/ Math.PI;
        finalBrg = finalBrg * 180/ Math.PI;
        backBrg = backBrg * 180 / Math.PI;

        if (lon2 < -180)
        {
            lon2 = lon2 + 360;
        }
        if (lon2 > 180)
        {
            lon2 = lon2 - 360;
        }

        if (finalBrg < 0)
        {
            finalBrg = finalBrg + 360;
        }
        if (finalBrg > 360)
        {
            finalBrg = finalBrg - 360;
        }

        Console.WriteLine("{0} {1}", lat2, lon2);

但是当运行导致 lat2 和 lon2 都为 NaN 时？

预期结果应为：

纬度：50°23′02″N 50.38384479 经度：3°55′49″W -3.93037298 最终轴承：224°59′32″ 224.99234101 后轴承：44°59′32″ 44.99234101

c# calc 是否做了一些不同的事情，我在那里犯了一个我看不到的错误:(

谢谢

Answer 1

根据维基百科 (https://en.wikipedia.org/wiki/Vincenty%27s_formulae) b = (1-f)a

所以你的第 10 行代码应该是：

   double b = a * (1 - 1 / flat);

注意 - 我没有检查任何其他问题，但它确实给出了非 NaN 的值。

C# 方位和距离计算（使用 Vincenty 的公式）结果为 NaN？

C# Bearing and distance calculation (using Vincenty's formula) results in NaN?

c#

wgs84

formula