在 C 中使用逆文森蒂公式的错误结果

Question

我已经编写了一个C脚本来实现逆Vincenty公式，以基于https://en.wikipedia.org/wiki/Vincenty%27s_formulae

然而，我的结果与这个在线计算器给出的结果不同https://www.cqsrg.org/tools/GCDistance/和谷歌地图。我的结果始终是在线计算器结果的1.18倍左右。

我的功能如下，任何关于我可能出错的提示都将不胜感激！

double get_distance(double lat1, double lon1, double lat2, double lon2)
{
double rad_eq = 6378137.0;             //Radius at equator
double flattening = 1 / 298.257223563; //flattenig of earth
double rad_pol = (1 - flattening) * rad_eq;     //Radius at poles

double U1,U2,L,lambda,old_lambda,sigma,sin_sig,cos_sig,alpha,cos2sigmam,A,B,C,u_sq,delta_s,dis;
//Convert to radians
lat1=M_PI*lat1/180.0;
lat2=M_PI*lat2/180.0;
lon1=M_PI*lon1/180.0;
lon2=M_PI*lon2/180.0;
//Calculate U1 and U2
U1=atan((1-flattening)*tan(lat1));
U2=atan((1-flattening)*tan(lat2));

L=lon2-lon1;

lambda=L;

double tolerance=pow(10.,-12.);//iteration tollerance should give 0.6mm
double diff=1.;

while (abs(diff)>tolerance)
{

sin_sig=sqrt(pow(cos(U2)*sin(lambda),2.)+pow(cos(U1)*sin(U2)-(sin(U1)*cos(U2)*cos(lambda)),2.));


cos_sig=sin(U1)*cos(U2)+cos(U1)*cos(U2)*cos(lambda);

sigma=atan(sin_sig/cos_sig);
alpha=asin((cos(U1)*cos(U2)*sin(lambda))/(sin_sig));
cos2sigmam=cos(sigma)-(2*sin(U1)*sin(U2))/((pow(cos(alpha),2.)));
C=(flattening/16)*pow(cos(alpha),2.)*(4+(flattening*(4-(3*pow(cos(alpha),2.)))));
old_lambda=lambda;
lambda=L+(1-C)*flattening*sin(alpha)*(sigma+C*sin_sig*(cos2sigmam+C*cos_sig*(-1+2*pow(cos2sigmam,2.))));
diff=abs(old_lambda-lambda);
}

u_sq=pow(cos(alpha),2.)*((pow(rad_eq,2.)-pow(rad_pol,2.))/(pow(rad_pol,2.)));
A=1+(u_sq/16384)*(4096+(u_sq*(-768+(u_sq*(320-(175*u_sq))))));
B=(u_sq/1024)*(256+(u_sq*(-128+(u_sq*(74-(47*u_sq))))));
delta_s=B*sin_sig*(cos2sigmam+(B/4)*(cos_sig*(-1+(2*pow(cos2sigmam,2.)))-(B/6)*cos2sigmam*(-3+(4*pow(sin_sig,2.)))*(-3+(4*pow(cos2sigmam,2.)))));
dis=rad_pol*A*(sigma-delta_s);

//Returns distance in metres
return dis;
}

Answer 1

这个公式不是对称的：

cos_sig = sin(U1)*cos(U2)
+ cos(U1)*cos(U2) * cos(lambda);

结果证明是错误的，缺少一个sin。

另一种格式样式(包括一些空白(也会有所帮助。

除了abs的fabs和cos的一个sin外；有两个CCD_ 6-调用，并且CCD_。

我插入了一个printf来查看该值的进展情况。

有些括号可以省略。这些公式真的很难实现。在这个嵌套数学运算的丛林中，一些更多的辅助变量可能会很有用。

do {
sin_sig = sqrt(pow(        cos(U2) * sin(lambda), 2)
+ pow(cos(U1)*sin(U2)
-    (sin(U1)*cos(U2) * cos(lambda))
, 2)
);
cos_sig = sin(U1) * sin(U2)
+ cos(U1) * cos(U2) * cos(lambda);
sigma = atan2(sin_sig, cos_sig);
alpha = asin(cos(U1) * cos(U2) * sin(lambda)
/ sin_sig
);
double cos2alpha = cos(alpha)*cos(alpha);      // helper var.
cos2sigmam = cos(sigma) - 2*sin(U1)*sin(U2) / cos2alpha;
C = (flat/16) * cos2alpha * (4 + flat * (4 - 3*cos2alpha));
old_lambda = lambda;
lambda = L + (1-C) * flat * sin(alpha)
*(sigma + C*sin_sig
*(cos2sigmam + C*cos_sig
*(2 * pow(cos2sigmam, 2) - 1)
)
);
diff = fabs(old_lambda - lambda);
printf("%.12fn", diff);
} while (diff > tolerance);

对于80，80，0，0，输出为(以公里为单位(：

0.000885870048
0.000000221352
0.000000000055
0.000000000000
9809.479224

其对应于具有WGS-84的毫米。