Finding Earth Coordinates (latitude,longitude), Distance (meters) and Bearing (angle)

Question

I need to deal with Earth coordinates in various ways. There is no function in C/C++ which does this straight away.
Referred below questions:

1. Python: Get lat/long given current point, distance and bearing
2. C: Given point of (latitude,longitude), distance and bearing, How to get the new latitude and longitude

From 1st one and movable type scripts website, I found that below are the formulas:

Find Bearing(angle) between 2 coordinates

x = cos(lat1Rad)*sin(lat2Rad) - sin(lat1Rad)*cos(lat2Rad)*cos(lon2Rad-lon1Rad);
y = sin(lon2Rad-lon1Rad) * cos(lat2Rad);
bearing = atan2(y, x);  // In radians; 
// Convert to degrees and for -ve add 360

Find Distance(meters) between 2 coordinates

PI = 3.14159265358979323846, earthDiameterMeters = 2*6371*1000;
x = sin((lat2Rad-lat1Rad) / 2);
y = sin((lon2Rad-lon1Rad) / 2);
meters = earthDiameterMeters * asin(sqrt(x*x + y*y*cos(lat1Rad)*cos(lat2Rad)));

Find Coordinate from Coordinate+Distance+Angle

meters *= 2 / earthDiameterMeters;
lat2Rad = asin(sin(lat1Rad)*cos(meters) + cos(lat1Rad)*sin(meters)*cos(bearing));
lon2Rad = lon1Rad + atan2(sin(bearing)*sin(meters)*cos(lat1Rad), 
                          cos(meters) - sin(lat1Rad)*sin(lat2Rad));

Below pseudo code should verify above 3 equations mutually:

struct Coordinate { double lat, lon; } c1, c2;  
auto degree = FindBearing(c1, c2);
auto meters = FindDistance(c1, c2);
auto cX = FindCoordiante(c1, degree, meters);

Now actually the answer comes almost near but not the correct. i.e. cX is Not equal to c2!
There is always a difference of 0.0005 difference in longitude value. e.g.

c1 = (12.968460,77.641308)  
c2 = (12.967862,77.653130)  
angle = 92.97         ^^^
distance = 1282.74  
cX = (12.967862,77.653613)
                      ^^^

I don't have much knowledge of mathematics' Havesine Forumla. But what I know is that from the website fcc.gov, the answer always comes correct.

What am I doing wrong?

Code only for reference

Though the syntax is in C++, but all the math functions are from C and is easily portable in C as well (hence tagged for both)

#include<iostream>
#include<iomanip>
#include<cmath>


// Source: http://www.movable-type.co.uk/scripts/latlong.html

static const double PI = 3.14159265358979323846, earthDiameterMeters = 6371.0 * 2 * 1000;

double degreeToRadian (const double degree) { return (degree * PI / 180); };
double radianToDegree (const double radian) { return (radian * 180 / PI); };

double CoordinatesToAngle (double latitude1,
                           const double longitude1,
                           double latitude2,
                           const double longitude2)
{
  const auto longitudeDifference = degreeToRadian(longitude2 - longitude1);
  latitude1 = degreeToRadian(latitude1);
  latitude2 = degreeToRadian(latitude2);

  using namespace std;
  const auto x = (cos(latitude1) * sin(latitude2)) -
                 (sin(latitude1) * cos(latitude2) * cos(longitudeDifference));
  const auto y = sin(longitudeDifference) * cos(latitude2);

  const auto degree = radianToDegree(atan2(y, x));
  return (degree >= 0)? degree : (degree + 360);
}

double CoordinatesToMeters (double latitude1,
                            double longitude1,
                            double latitude2,
                            double longitude2)
{
  latitude1 = degreeToRadian(latitude1);
  longitude1 = degreeToRadian(longitude1);
  latitude2 = degreeToRadian(latitude2);
  longitude2 = degreeToRadian(longitude2);

  using namespace std;
  auto x = sin((latitude2 - latitude1) / 2), y = sin((longitude2 - longitude1) / 2);
#if 1
  return earthDiameterMeters * asin(sqrt((x * x) + (cos(latitude1) * cos(latitude2) * y * y)));
#else
  auto value = (x * x) + (cos(latitude1) * cos(latitude2) * y * y);
  return earthDiameterMeters * atan2(sqrt(value), sqrt(1 - value));
#endif
}

std::pair<double,double> CoordinateToCoordinate (double latitude,
                                                 double longitude,
                                                 double angle,
                                                 double meters)
{
  latitude = degreeToRadian(latitude);
  longitude = degreeToRadian(longitude);
  angle = degreeToRadian(angle);
  meters *= 2 / earthDiameterMeters;

  using namespace std;
  pair<double,double> coordinate;

  coordinate.first = radianToDegree(asin((sin(latitude) * cos(meters))
                             + (cos(latitude) * sin(meters) * cos(angle))));
  coordinate.second = radianToDegree(longitude
                    + atan2((sin(angle) * sin(meters) * cos(latitude)),
                    cos(meters) - (sin(latitude) * sin(coordinate.first))));

  return coordinate;
}

int main ()
{
  using namespace std;
  const auto latitude1 = 12.968460, longitude1 = 77.641308,
             latitude2 = 12.967862, longitude2 = 77.653130;

  cout << std::setprecision(10);
  cout << "(" << latitude1 << "," << longitude1 << ") --- "
          "(" << latitude2 << "," << longitude2 << ")\n";

  auto angle = CoordinatesToAngle(latitude1, longitude1, latitude2, longitude2);
  cout << "Angle =  " << angle << endl;

  auto meters = CoordinatesToMeters(latitude1, longitude1, latitude2, longitude2);
  cout << "Meters = " << meters << endl;

  auto coordinate = CoordinateToCoordinate(latitude1, longitude1, angle, meters);
  cout << "Destination = (" << coordinate.first << "," << coordinate.second << ")\n";
}

Answer 1

In CoordinateToCoordinate you use sin(coordinate.first) which is already in degrees. Use sin(degreeToRadian(coordinate.first)).

Or to be more cleaner:

... CoordinateToCoordinate (...)
{
  ...
  coordinate.first = asin((sin(latitude) * cos(meters))
                        + (cos(latitude) * sin(meters) * cos(angle)));
  coordinate.second = longitude + atan2((sin(angle) * sin(meters) * cos(latitude)), 
         cos(meters) - (sin(latitude) * sin(coordinate.first)));

  coordinate.first = radianToDegree(coordinate.first);
  coordinate.second = radianToDegree(coordinate.second);

  return coordinate;
}

This fixes the problem. Live Demo.

Answer 2

Partial answer

With an angle of 92.97° then converted to radians, a call to sin/cos/tan will effectively change the angle to 2.97°. This step alone loses 6 bits of precision as the period reduction occurs after degrees to radians conversion and in the trig function call.

Precision of trigonometric functions with large angles in degrees can be enhanced. Use the fortunate aspect the there are exactly 360.0 degrees in a circle. Perform a "modulo 45°", maybe with remquo(angle, 45, &octant), and then converting to radians before calling the trig function with a radian argument.

Example sind()

---

Your answers of 77.641308 and 77.653130 differs by about 1 part in 6600 (~13 bits of precision). This answer may not fully explain that, yet should help. (If some usage of float occurs somewhere, that should be made double.)