Inaccurate Result from Haversine's Bearing Calculation

Question

I am trying to implement the Haversine Formula in a little GPS program I'm writing. The distance calculations appear to be spot-on. However, I believe the bearing is being computed in radians, and I don't know how to properly convert the result to compass directions (0 for North, 90 for East, etc).

Any help would be greatly appreciated, all this talk of cosigns and arctangents is giving me a major headache! I'm a coder, not a mathematician!

#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

It returns a bearing of 65 from LA to NY, and a bearing of 273 from NY to LA.

If we add the bearings together, we get 338 which can't be right - shouldn't it equal 360?
Or am I completely out to lunch?

Anyway, as you can see, I always compute both distance and bearing at the same time. If you could also suggest a way to clean up the code so it doesn't perform unnecessary calculations, that would be so very outstanding! I'm running this on a small microprocessor where I like to make every cycle count!

Answer 1

Not a problem.

Consider 2 locations at the same latitude, yet differ in longitude. One is not due east (90) of the other when going in the shorting great circle route. Nor is the first due (west 270). The bearings are not necessarily complementary.

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 adds more good information.

Per the site suggest by @M Oehm your code is about correct.

---

Per OP request, some mods that may have slight speed improvement.

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