Calculating distance between latitude longitudes.

[latecomer]

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.

Any suggestions/advices are welcome.

Thanks in advance.

 

[xan-der]

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

// 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:

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

   w *= EARTH_AVERAGE_RADIUS;

   return w;
}

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

Xander
[/code]