Нахождение координат Земли (широта, долгота), расстояние (в метрах) и азимут (угол)

Мне нужно иметь дело с координатами Земли по-разному. В C / C ++ нет функции, которая делает это сразу.
Ниже приведены вопросы:

1. Python: получить широту / долготу с учетом текущей точки, расстояния и
   подшипник
2. C: Заданная точка (широта, долгота), расстояние и азимут, как
   чтобы получить новую широту и
   долгота

С 1-го и сайт скриптов подвижного типа, Я обнаружил, что ниже приведены формулы:

Найти подшипник (угол) между 2 координатами

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

Найти расстояние (в метрах) между 2 координатами

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

Найти координату из координаты + расстояние + угол

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

Ниже псевдокод должен проверить выше 3 уравнения взаимно:

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

Теперь на самом деле ответ приходит почти рядом, но не правильно. то есть cX не равно c2!
Всегда есть разница 0.0005 разница в значении долготы.
например

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

Я не очень разбираюсь в математике Havesine Forumla. Но я знаю, что сайт fcc.gov, ответ всегда приходит правильно.

Что я делаю неправильно?

Код только для справки

Хотя синтаксис в C ++, но все математические функции взяты из C и легко переносимы в C (следовательно, помечены для обоих)

#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";
}