Calculating coordinates given a bearing and distance issues

Question

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

What might be causing this issue?

// Value for pi
#define PI 3.14159265358979323846
// Mean radius of earth in m, from https://en.wikipedia.org/wiki/Earth_radius#Average_distance_from_center_to_surface
#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);
double deg2rad(double);
double rad2deg(double);

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 rlat1 = deg2rad(Location_0->lat);
    double rlon1 = deg2rad(Location_0->lon);
    double rbearing = deg2rad(azimuth);
    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;
    }

    Location_1->lat = rad2deg(rlat2);
    Location_1->lon = rad2deg(rlon2);

}

double deg2rad(double d){
    return d*PI/180.0;
}

double rad2deg(double r){
    return r*180.0/PI;
}