jremington:
Great circle calculations won't work as well as the Haversine function, due to the inaccuracy of the single precision trigonometric functions.
I think for very small distances it should be no problem to calculate the course bearing using Pythagoras' theorem instead of using great circle trigonometric calculations.
jremington:
However, due to the spherical approximation for the ellipsoidal Earth, all these approximations suffer from about 3 meters/km, or 0.3% error on average for distances greater than a few km. Considering that it is difficult to get GPS coordinates for any location with positional accuracy better than 2-3 m, it may not be worth the effort to take this much further.
OK, to keep things simple I left the distance calculations based on earth to be a perfect sphere and not a WGS84 ellipsoid. Of course, this will create a small percentage of error.
But the main thing is, that the calculations are now highly precise using an 8-bit Atmega, without the need of a highly accurate 8-byte double data type.
I just did some comparison of the tinyGPSplus distance calculation against my distance calculation. What I found is that:
tinyGPSplus method with 'float' variable input:
- calculated distance is often staying the same with very small coordinate differences
- calculated distance sometimes is creating a big jump in the distance
- calculation error for distance is always smaller than 1,00 meter (when fed with accurate float input)
optimized haversine method with 'long' variable input:
- even the smallest difference in coordinates lead to a small calculated distance
- calculated distance is proportional to difference in coordinates
- no significant calculation error compared against 64-bit 'double' math calculations
Here is some code for comparison:
void nmeaDegrees2Decimal(char* nmeaDegrees, long &value)
{ // converts NMEA encoded degrees into long value
// the last digit of the returned long "value" represents 0.0000001 degrees
// i.e. NMEA string "5133.8123" ==> 51 degrees and 33.8123 minutes
// = 51.5635383 degrees ==> returned value == 515635383
value= atol(nmeaDegrees); // retrieves 5133 (= 51 degrees, 33 minutes)
long arcminutes=value%100; // 33
value=(value/100)*10000000; // 510000000
nmeaDegrees= strchr(nmeaDegrees,'.'); // find decimal point in string
// handle fractions of minutes (up to 5 digits after decimal point)
for (byte b=0;b<5;b++) // five times multiply by 10 and add digit
{
if (isdigit(nmeaDegrees[1]))
{
arcminutes*=10;
arcminutes+=nmeaDegrees[1]-'0';
nmeaDegrees++;
}
else arcminutes*=10;
}
arcminutes= (5*arcminutes+1)/3; // 6000000 arc minute units = 1 degree = 10000000 decimal units, do conversion with rounding
value+= arcminutes; // add fraction units for final value
}
double calcDistanceHaversine(long lat1, long lon1, long lat2, long lon2){
// lat/lon params must be 'long' values in 10E-7 degree units,
// i.e. 12.0000000° = 120000000
// 'double' is actually 'float' with 8-bit Atmega controllers
double degreeUnits= 1E-7;
double dlon,dlat,a,c;
int earthradius = 6371; /*km*/
dlon = radians((lon2-lon1)*degreeUnits);
dlat = radians((lat2-lat1)*degreeUnits);
a = sin(dlat/2)*sin(dlat/2) + sin(dlon/2)*sin(dlon/2) * cos(radians(lat1*degreeUnits)) * cos(radians(lat2*degreeUnits));
c = 2 * atan2(sqrt(a), sqrt(1-a));
return earthradius*c;
}
double calcDistanceTinyGPSPlus(double lat1, double long1, double lat2, double long2)
{
// returns distance in meters between two positions, both specified
// as signed decimal-degrees latitude and longitude. Uses great-circle
// distance computation for hypothetical sphere of radius 6372795 meters.
// Because Earth is no exact sphere, rounding errors may be up to 0.5%.
// Courtesy of Maarten Lamers
double delta = radians(long1-long2);
double sdlong = sin(delta);
double cdlong = cos(delta);
lat1 = radians(lat1);
lat2 = radians(lat2);
double slat1 = sin(lat1);
double clat1 = cos(lat1);
double slat2 = sin(lat2);
double clat2 = cos(lat2);
delta = (clat1 * slat2) - (slat1 * clat2 * cdlong);
delta = sq(delta);
delta += sq(clat2 * sdlong);
delta = sqrt(delta);
double denom = (slat1 * slat2) + (clat1 * clat2 * cdlong);
delta = atan2(delta, denom);
return delta * 6372795;
}
void setup()
{
Serial.begin(9600);
Serial.println();
long lat1, lon1, lat2, lon2;
nmeaDegrees2Decimal("5100.00000", lat1); // 51 degrees
nmeaDegrees2Decimal("1000.00000", lon1); // 10 degrees;
lat2=lat1;
lon2=lon1;
Serial.println("i\tlat1\t\tlon1\t\tlat2\t\tlon2\t\ttinyGPS\tjursGPS\tdifference");
for (int i=0;i<=1000;i++)
{
Serial.print(i);Serial.print('\t');
Serial.print(lat1);Serial.print('\t');
Serial.print(lon1);Serial.print('\t');
Serial.print(lat2);Serial.print('\t');
Serial.print(lon2);Serial.print('\t');
double f=1E-7;
double distance1=calcDistanceTinyGPSPlus(lat1*f, lon1*f, lat2*f, lon2*f);
Serial.print(distance1);Serial.print('\t');
double distance2=1000*calcDistanceHaversine(lat1, lon1, lat2, lon2);
Serial.print(distance2);Serial.print('\t');
Serial.println(abs(distance2-distance1));
if (lat2%3==0) lat2+=2; else lat2++;
lon2+=1;
}
}
void loop()
{
}
TinyGPS is great for getting GPS data but the calculations are useless for distances < 100m )