Problem with the Haversine Formula in C++ Help

[DeepSeaBiscuit]

Hello, I have the haversine function within my C++ program to calculate the distance between to points, the two tester points I have chosen are Swansea and Cardiff both located in Wales.
Swansea lat = 51.622559, long = -3.934534
Cardiff lat = 51.475661, long = -3.174688

The problem I am having is that the actual line distance for these two points is 55.02 km but my program keeps giving out 28.30 km and I am not sure why?

My function that I am using to calculate the distance is #include <cmath>
using namespace std;

double calcDistance(double lat2, double lat1, double long2, double long1)
{
double toRadians = 3.1415 / 180;
double dLat = (lat1-lat2)*toRadians;
double dLong = (long1-long2)*toRadians;
double a = pow(sin(dLat / 2), 2) + cos(lat1)*cos(lat2)*pow(sin(dLong/ 2),2);
double c = 2 * atan2(sqrt(a), sqrt(1-a));
double R = 6371;
double distance = R * c;
return distance;
}

Any help?

 

[Mark44]

Hello, I have the haversine function within my C++ program to calculate the distance between to points, the two tester points I have chosen are Swansea and Cardiff both located in Wales.
Swansea lat = 51.622559, long = -3.934534
Cardiff lat = 51.475661, long = -3.174688

The problem I am having is that the actual line distance for these two points is 55.02 km but my program keeps giving out 28.30 km and I am not sure why?

My function that I am using to calculate the distance is#include <cmath>
using namespace std;

double calcDistance(double lat2, double lat1, double long2, double long1)
{
double toRadians = 3.1415 / 180;
double dLat = (lat1-lat2)*toRadians;
double dLong = (long1-long2)*toRadians;
double a = pow(sin(dLat / 2), 2) + cos(lat1)*cos(lat2)*pow(sin(dLong/ 2),2);
double c = 2 * atan2(sqrt(a), sqrt(1-a));
double R = 6371;
double distance = R * c;
return distance;
}

Any help?

What you calculate for a looks OK as far as it goes. The wikipedia page (http://en.wikipedia.org/wiki/Haversine_formula) differs from what you have for the code after your calculation of a.

It gives this for the great circle distance: d = 2r * arcsin($\sqrt{a}$).

BTW, your value for $\pi$ is not very precise. You could do much better using 4.0 * atan(1.0).

 

[DeepSeaBiscuit]

I have altered the function according to what you suggested however this still hasn't solved the issue some how? As a calculated distance I get 5926.76 km which isn't too far off the radius of the earth! I have been stuck on this problem for a day now I really can't figure out the problem? Could there be something else I'm missing?

double toRadians = 3.14159265359 / 180;
double dLat = (lat1-lat2)*toRadians;
double dLong = (long1-long2)*toRadians;
double a = pow(sin(dLat / 2), 2) + cos(lat1)*cos(lat2)*pow(sin(dLong/ 2),2);

double R = 6371;
double distance = 2*R*asin(sqrt(a));

return distance;

 

[DeepSeaBiscuit]

I may have figured the problem out, I converted the lat and long individually from deg to radians and the program works perfectly! Thank you for the help!

 

[jim mcnamara]

Your value for R will cause you grief - the Earth is a spheroid - like a ball with middle age spread. The WGS-84 ellipsoid has:
Equatorial radius (6,378.1370 km)
Polar radius (6,356.7523 km)

Which, going from the equator to the poles, with all other considerations ignored is 3m/km. Beyond a few hundred km this error can become unacceptably large. Which is the exact idea Mark44 is trying to convey about the PI constant you have. In C/C++ the default precision for a double is 15 digits. So you may want to tweak your code to actually work to something closer to that level of accuracy. The radius for your example should be correct to 8 significant digits. Ditto PI. Not four. Nine would be better for all possible inputs.

Google for Meridional radius or Gauss radius of curvature maybe. Usually a decent app using haversines will calculate distances over small changes in Lat with a recalculated radius for each Lat delta.