# 给定方位角和距离问题，计算坐标

I've been attempting to implement the function described here in C++. It seems to work pretty well, but breaks down at large vales of `distance`, or when close to a pole. The angle at which the abrupt change in the first test case happens is around 243 degrees, which doesn't seem significant mathematically.

I'm using hamstermap to map out the function output and jdoodle to quickly compile/run the code - I can't zoom out enough to get good pictures, but to recreate the issue, here are a couple of test cases that demonstrate it (enter values in the function call in main() on line 35):

1: Set `distance` to 500,000, `startingLocation.lat` to -88 and `startingLocation.lon` to 173. Note that the distance value presents no problem at latitude values closer to the equator.

2: Set `distance` to 10,000,000, `startingLocation.lat` to -15 and `startingLocation.lon` to 173

``````// Value for pi
#define PI 3.14159265358979323846
#define REARTH 6371008.7714
// Threshold for equality in floating point operations
#define EPSILON 0.000001

#include <iostream>
#include <iomanip>
#include <string>
#include <cmath>
using namespace std;

struct LLlocation
{
double lat;
double lon;
};

void CalculatePosition(struct LLlocation *Location_0, struct LLlocation *Location_1, double, double, double=REARTH);

struct LLlocation startingLocation, finalLocation;

int main()
{
startingLocation.lat = -88.000000000;
startingLocation.lon = 173.000000000;

cout << fixed << setprecision(9);

for (int i = 0; i <= 360; i++)
{
CalculatePosition(&startingLocation, &finalLocation, 500000, i);
//cout << "Leaving at " << i << " degrees\n";
cout << finalLocation.lat << "," << finalLocation.lon << "\n";
}
}

void CalculatePosition(struct LLlocation *Location_0, struct LLlocation *Location_1, double distance, double azimuth, double rearth)
{
/*
INPUTS
Location_0.lat = starting latitude in decimal degrees
Location_0.lon = starting longitude in decimal degrees
distance = distance in meters
azimuth = angle from North in decimal degrees
*/

// Convert angles to radians and normalize distance
double rdistance = distance/rearth;

// Use law of haversines to calculate new lat and lon
double rlat2 = asin(sin(rlat1)*cos(rdistance)+cos(rlat1)*sin(rdistance)*cos(rbearing));

double rlon2;
if(cos(rlat2) == 0 || abs(cos(rlat2) < EPSILON))
{
//endpoint is a pole
rlon2 = rlon1;
}
else
{
rlon2 = (fmod((rlon1-asin(sin(rbearing)*sin(rdistance)/cos(rlat2))+PI),(2*PI)))-PI;
}