# Calculating distance between latitude longitudes.

Status
Not open for further replies.

#### latecomer

##### Member level 4
Hello all,
I need to calculate the distance between two sets of Longitudes and Latitudes. I need this for a GPS based system which I am designing which will show the distance travelled by a vehicle.

Can anyone help me out in this ? I have seen on the net that Haversine formula is good but I am having problem porting it to my PIC16F877A controller using HT-PIC compiler.

Assuming that you want the geodetic distance, this routine works well and with good speed:

Code:
// Lat/Long inputs in radians, range output in meters
/* This version of the distance between two points on the earth
is based on the formula given in Appendix A of Sodano, E.M.,
"General non-iterative solution of the inverse and direct geodetic
problems", Bulletin Geodesique, vol 75 (March 1965), pp 69-84.
*/
// Constants used by Sodano's method
const double A    =		6378137.0;	// semi-major axis of ellipsoid
const double RECF =		298.257223563;	// reciprocal of flattening (1/F)
const double F    =		(1.0 / RECF);	// flattening f
const double ES  = ((2.0 - F) * F);		// First eccentricity squared
const double B   = sqrt(A * A * (1.0 - ES));		// Semi-major axis
const double F2 = F * F;
const double F2_2 = 0.5 * F2;
const double F_1 = 1.0 - F;

static void
ComputeRangeSodano(double pt1Lat, double pt1Lon,
double pt2Lat, double pt2Lon,
double *range)
{
double L = pt2Lon - pt1Lon;

const double TESTV = 1.0E-11;
if(fabs(pt1Lat - pt2Lat) < TESTV && fabs(L) < TESTV)
{
*range = 0.0;
return; // Points coincident
}

double beta1 = atan(F_1 * tan(pt1Lat));
double beta2 = atan(F_1 * tan(pt2Lat));

double a = sin(beta1) * sin(beta2);
double b = cos(beta1) * cos(beta2);

double cos_phi = a + b * cos(L);

double temp1, temp2;
double sin_phi, phi;

temp1 = sin(L) * cos(beta2);
temp2 = sin(beta2) * cos(beta1) - sin(beta1) * cos(beta2) * cos(L);
sin_phi = sqrt(temp1 * temp1 + temp2 * temp2);
phi = atan2(sin_phi, cos_phi);

double c = (b * sin(L)) / sin_phi;

double m = 1.0 - c * c;

double csc_phi = 1.0 / sin_phi;
double cot_phi = cos_phi / sin_phi;

double cos_phi2 = cos_phi * cos_phi;
double sin_cos_phi = sin_phi * cos_phi;

double phi2 = phi * phi;

double Sdb =
(1.0 + F + F2) * phi
+ a * ((F + F2) * sin_phi - F2_2 * phi2 * csc_phi)
+ m * 0.5 * (-(F + F2) * (phi + sin_cos_phi) + F2 * phi2 * cot_phi)
- a * a * F2_2 * (sin_phi * cos_phi)
+ m * m * F2_2 * (0.125 * (phi + sin_cos_phi) - phi2 * cot_phi - 0.25 * sin_cos_phi * cos_phi2)
+ a * m * F2_2 * (phi2 * csc_phi + sin_phi * cos_phi2);

double S = Sdb * B;

*range = S;
}

A geocentric calculation can be performed much more quickly and might even be ok for points that are close together. It looks like this:

Code:
// Geocentric distance constants
const double EARTH_AVERAGE_RADIUS = 6367443.8; // Average radius of the earth in Meters

// Distance on a perfect sphere using simple trig
static double
geocentric_ground_distance1(double Lat1, double Lon1, double Lat2, double Lon2)
{
double w;

w = acos(cos(Lat1)*cos(Lat2)*cos((Lon2-Lon1)) + sin(Lat1)*sin(Lat2));

return w;
}

You would need to convert from geodetic to geocentric before calling.

Xander
[/code]

Status
Not open for further replies.