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!

You're out to lunch. Assuming bearings that are in opposite directions, subtracting the smaller bearing from the larger bearing should yield 180. Adding the bearings results in a number between 180 and 540. For example, 45 and 225 are opposites, the sum is 270, but the difference is 180. Also, when moving along a great circle, the bearing is not constant. So starting from LA, you begin with a north-eastern direction, and end going almost due east. Starting in NY, you begin almost due west, but end in a south-western path. — user3386109, Mar 30 '16 at 18:57
On the curved earth surface, I do not think the bearings need to be 180 degree complements. — chux - Reinstate Monica, Mar 30 '16 at 19:01
In addition to opposite directions differing by 180°, the bearing as angle to the meridians changes as you move along the great circles, so that the incoming bearing isn't the opposite of the outgoing bearing, [see here](http://www.movable-type.co.uk/scripts/latlong.html). The notion of opposite angles comes from mentally imagining the shortest distance as a straight line on a map, but the shortest path isn't a straight line in most projections. — M Oehm, Mar 30 '16 at 19:07

score 3 · Answer 1 · edited May 23 '17 at 10:28

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

Thanks friends for all the input. I will count my implementation as correct then, at least, correct enough for now. It sounds like a couple of you have some great thoughts on how the code could be optimized but, alas, I do not know how to incorporate the changes in my program. Never was terribly good at trig. Happy Thursday everyone! — programer8472, Mar 31 '16 at 20:52
@programer8472 Some minor improvements posted, not likely to make things much faster. 2 more ideas: `use float` and `float` functions `sinf()`, etc. Post _working_ code on [Code Review](http://codereview.stackexchange.com/questions/ask) and ask about _performance_ considerations. Only leave in `printf()` if that is part of the true code. — chux - Reinstate Monica, Mar 31 '16 at 21:20

Inaccurate Result from Haversine's Bearing Calculation

1 Answers1