Haversine 的方位角计算结果不准确

Question

我正在尝试在我正在编写的一个小 GPS 程序中实现 Haversine 公式。距离计算似乎是准确的。但是，我相信方位是以弧度计算的，我不知道如何正确地将结果转换为指南针方向（0 表示北，90 表示东等）。

任何帮助都将不胜感激，所有这些关于余符号和反正切的讨论都让我头疼！我是程序员，不是数学家！

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void FindDistance (double latHome, double lonHome, double latDest, double lonDest)
{
    double pi=3.141592653589793;
    int R=6371; //Radius of the Earth in kilometers
    //Keep the parameters passed to the function immutable
    double latHomeTmp=(pi/180)*(latHome);
    double latDestTmp=(pi/180)*(latDest);
    double differenceLon= (pi/180)*(lonDest- lonHome);
    double differenceLat=(pi/180)*(latDest- latHome);
    double a= sin (differenceLat/2.)*sin (differenceLat/2.)+cos   (latHomeTmp)*cos (latDestTmp)*sin (differenceLon/2.)*sin (differenceLon/2.);
    double c=2*atan2 (sqrt (a), sqrt (1-a));
    double Distance=R*c;
    printf ("Distance is %f\n", Distance);

    double RadBearing=atan2 (sin (differenceLon)*cos (latDestTmp), cos   (latHomeTmp)*sin (latDestTmp)-sin (latHomeTmp)*cos (latDestTmp)*cos     (differenceLon));
    double DegBearing=RadBearing*57.2958;
    if (DegBearing<0) DegBearing=360+DegBearing;
    printf ("Bearing is %f\n", DegBearing);
} //Function FindDistance

int main (void) {
    puts ("LA to NY");
    FindDistance (34.052235, -118.243683, 40.748817, -73.985428);
    puts ("NY to LA");
    FindDistance (40.748817, -73.985428, 34.052235, -118.243683);
} //Function main

gcc -o gps -lm gps.c

从 LA 到 NY 的方位角为 65，从 NY 到 LA 的方位角为 273。

如果我们将轴承加在一起，我们会得到 338，这是不可能的 - 它不应该等于 360 吗？
还是我完全出去吃午饭了？

无论如何，如您所见，我总是同时计算距离和方位。如果您还可以建议一种清理代码的方法，使其不会执行不必​​要的计算，那将非常出色！我在一个小型微处理器上运行它，我喜欢计算每个周期！

Answer 1

没问题。

考虑相同纬度但经度不同的 2 个位置。在做空大圈路线时，一个不在另一个的正东（90）。也不是第一次到期（西 270）。轴承不一定是互补的。

FindDistance(34.0, -118.0, 34.0, -73.0);
FindDistance(34.0, -73.0, 34.0, -118.0);

Distance is 4113.598081
Bearing is 76.958824
Distance is 4113.598081
Bearing is 283.041176

@user3386109 增加了更多好的信息。

根据@M Oehm 的site 建议，您的代码基本正确。

---

根据 OP 请求，一些模组可能会略微提高速度。

void FindDistance(double latHome, double lonHome, double latDest,
    double lonDest) {
  // A few extra digits sometimes is worth it - rare is adding more digits a problem
  // double pi=3.141592653589793;
  static const double pi_d180 = 3.1415926535897932384626433832795 / 180;
  static const double d180_pi = 180 / 3.1415926535897932384626433832795;

  // int R=6371; //Radius of the Earth in kilometers
  // (Compiler may do this all ready)
  static const double R = 6371.0; // better to make FP to avoid the need to convert 

  //Keep the parameters passed to the function immutable
  double latHomeTmp = pi_d180 * (latHome);
  double latDestTmp = pi_d180 * (latDest);
  double differenceLon = pi_d180 * (lonDest - lonHome);
  double differenceLat = pi_d180 * (latDest - latHome);

  double a = sin(differenceLat / 2.) * sin(differenceLat / 2.)
      + cos(latHomeTmp) * cos(latDestTmp) * sin(differenceLon / 2.)
          * sin(differenceLon / 2.);

  double c = 2 * atan2(sqrt(a), sqrt(1 - a));
  double Distance = R * c;
  printf("Distance is %f\n", Distance);

  double RadBearing = atan2(sin(differenceLon) * cos(latDestTmp),
      cos(latHomeTmp) * sin(latDestTmp)
          - sin(latHomeTmp) * cos(latDestTmp) * cos(differenceLon));

  // double DegBearing = RadBearing * 57.2958;
  double DegBearing = RadBearing * d180_pi;

  // Why is this even needed?
  if (DegBearing < 0) DegBearing = 360 + DegBearing;

  printf("Bearing is %f\n", DegBearing);
} //Function FindDistance