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:

Python: Get lat/long given current point, distance and bearing
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";
}

And you think that will compile as C? It **is** not C. There is no "a little bit pregnant"! — too honest for this site, Aug 19 '15 at 13:35
@Olaf, The code is for additional reference, not the center point of the question. It won't compile in C as it is, but it's quite easily portable. The other question I referred is actually from C. To get wider audience, both C and C++ are tagged. — iammilind, Aug 19 '15 at 13:36
`double` ... floating point arithmetic is potentially imprecise. — Daniel Jour, Aug 19 '15 at 14:02
Try converting the numbers into fixed point representation to capture the full given precision and calculate it. But even before, calculate it manually. — Eugene Sh., Aug 19 '15 at 14:03
@EugeneSh., I din't get your comment. At which particular equation are you referring? Can you give an example of converting the number into fixed point representation? There can be inaccuracies due to `double`, but `0.0005` is quite large in that sense. The result is same for `double` and `long double`. — iammilind, Aug 19 '15 at 14:08
At all of them. If you want to check whether it is an issue of the double precision and not the calculation itself, you should first try more precise methods. The large *resulting* inaccuracy is accumulated through the steps. https://en.wikipedia.org/wiki/Fixed-point_arithmetic — Eugene Sh., Aug 19 '15 at 14:09
To test whether it's floating point rounding errors, you could try to change the rounding mode: http://stackoverflow.com/a/6867722/1116364 This won't solve your issue, but may confirm that the imprecision is the actual cause. — Daniel Jour, Aug 19 '15 at 15:56

score 2 · Accepted Answer · edited Aug 20 '15 at 06:39

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.

score 0 · Answer 2 · edited May 23 '17 at 10:29

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.)

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

Code only for reference

2 Answers2