Continue to Site

Welcome to EDAboard.com

Welcome to our site! EDAboard.com is an international Electronics Discussion Forum focused on EDA software, circuits, schematics, books, theory, papers, asic, pld, 8051, DSP, Network, RF, Analog Design, PCB, Service Manuals... and a whole lot more! To participate you need to register. Registration is free. Click here to register now.

Calculating distance between latitude longitudes.

Status
Not open for further replies.

latecomer

Member level 4
Member level 4
Joined
Jun 6, 2001
Messages
77
Helped
3
Reputation
6
Reaction score
1
Trophy points
1,288
Location
USA
Activity points
642
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.
 

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

   w *= EARTH_AVERAGE_RADIUS;

   return w;
}

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

Xander
[/code]
 

Status
Not open for further replies.

Similar threads

Part and Inventory Search

Welcome to EDABoard.com

Sponsor

Back
Top