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

i need deal earth coordinates in various ways. there no function in c/c++ straight away.
referred below questions:

1. python: lat/long given current point, distance , bearing
2. c: given point of (latitude,longitude), distance , bearing, how new latitude , longitude

from 1st 1 , movable type scripts website, found below 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 degrees , -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 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 answer comes almost near not correct. i.e. cx not equal c2!
there 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 knowledge of mathematics' havesine forumla. know the website fcc.gov, answer comes correct.

what doing wrong?

code reference

though syntax in c++, math functions c , portable in c (hence tagged 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"; } 

in coordinatetocoordinate use sin(coordinate.first) in degrees. use sin(degreetoradian(coordinate.first)).

or 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 problem. live demo.