double precision / TinyGPS++

PaulS is right

however you can split the whole part (use int) and decimal part (use unsigned long in steps of 1E-8),

especially if you only add this could give you a " 4.8 "format

but as PaulS says, you must rewrite the math.

phil_rudd69:
I need to submit double values with 8 decimal places like "12.12345678" via Serial Monitor to my Arduino Mega 2560. For this purpose I read the serial buffer and put the seperate entrys (1,2,.,1,.2,3,4...) into a char-array.

Arduino 'float' precision is 6-7 significant digits.
Arduino 'double' is exactly the same as 'float', there is no difference with 8-bit Atmegas.

If you need to handle more significant digits, you can easily use 'long' for example.
The 'long' variable type offers 9-10 significant digits.

So if you receive "12.12345678", ignore the decimal point and read "1212345678" into a long.

Then it would be easy to compare that value against a destination value of let's say 1212345675.

So if you receive "12.12345678", ignore the decimal point and read "1212345678" into a long.

In your scheme, how do you plan to keep track of the decimal point? It cannot be ignored.

Why we cannot have doubles in arduino? The lib is about 1.5x larger and calcs are only 2x slower. Just an idea for a future release :slight_smile:

jremington:
In your scheme, how do you plan to keep track of the decimal point? It cannot be ignored.

I'd suggest that the lowest digit in the 'long' value represents 1E-7 degrees = 0.0000001°

So latitude -90.0000000 ... 90.0000000 is actually a long from -900000000 ... 900000000

And longitude from -179.9999999 ... 180.0000000 is actually a long from -1799999999 ... 1800000000

So all valid latitude/longitude values will fit with the long range from −2147483648 up to 2147483647

If Serial sends less than 7 digits after the decimal point, fill up with 0s.
If Serial sends more than 7 digits after the decimal point, ignore (or round the last valid digit).

Calculating small distances is easy then:

  1. Either use Pythagoras’ theorem (using float calculation then) after finding out how much distance 1 unit of latitude and longitude is in kilometers/meters at the given location. For small distances there is no big difference to the "correct" spheric distance value when calculating a flat triangle instead.

  2. Or use the haversine formula, which is an algorithm that gives very accurate results even if the platform has limited floating point numbers calculating capabilities.

With Haversine the resolution with 6 digits lat/long is something like ~10m. For example lat/long distance

15.0000 50.0000
15.0000 50.0001
gives you 10m with double

15.0000 50.0000
15.0001 50.0001
gives you 13m with single (15m with double).

// Haversine formula for distance[m] and initial bearing[degree]
// of 2 points on Earth based on LAT[degree] and LON[degree]
// Pito 6/2015

#define R 6371000.0
#define PI 3.1415926535897932384626433832795

// lat and long in decimal degree, distance in meters
double distance (double lat1, double lon1, double lat2, double lon2) {
double phi1 = lat1 * PI / 180.0;
double phi2 = lat2 * PI / 180.0;
double dphi = (lat2 - lat1) * PI / 180.0;
double dlambda = (lon2 - lon1)* PI / 180.0;
double a1 = sin(dphi/2.0)*sin(dphi/2.0);
double a2 = cos(phi1)*cos(phi2);
double a3 = sin(dlambda/2.0)*sin(dlambda/2.0);
double a = a1 + (a2 * a3);
double c = 2.0 * atan2(sqrt(a), sqrt(1.0-a));
double d = R * c;
return d;
}

double bearing (double lat1, double lon1, double lat2, double lon2) {
double phi1 = lat1 * PI / 180.0;
double phi2 = lat2 * PI / 180.0;
double dlambda = (lon2 - lon1) * PI / 180.0;
double y = sin(dlambda)*cos(phi2);
double x = cos(phi1)*sin(phi2)-sin(phi1)*cos(phi2)*cos(dlambda);
double brn = (atan2(y,x)) * 180.0 / PI;
return brn;
}

// lat1  lon1    lat2  lon2 (degree)
// 50.0 15.0     50.0001 15.0001   (13.2185m, 32.73222deg)

I made something crazy (and perhaps relevant) here.
http://forum.arduino.cc/index.php?topic=335534.0

It maybe isn't the best, but it appears to work.
I could improve it, maybe, but not today.

phil_rudd69:
I ask because I also use the TinyGPS++ library within my project and I think this library can handle double values with more then 6 decimal places correctly. But that could be wrong as well :slight_smile:
Unfortunately I don't understand what is going on within the functions of TinyGPS++.

I've just looked up the source code of TinyGPSplus and found, that the library tries to keep the received data in a very precise data structure, which is declared like this:

struct RawDegrees
{
   uint16_t deg;
   uint32_t billionths;
   bool negative;
public:
   RawDegrees() : deg(0), billionths(0), negative(false)
   {}
};

Despite of that the "billionths" part of the RawDegrees is a fake with 8-bit Atmega controllers, this is not the worst thing and could lead to accurate calculations and output.

Unfortunately, not a single function of that library is making any use of the precision provided in the raw data. As soon as you retrieve the data from the library, this will happen (latitude example):

double TinyGPSLocation::lat()
{
   updated = false;
   double ret = rawLatData.deg + rawLatData.billionths / 1000000000.0;
   return rawLatData.negative ? -ret : ret;
}

The (precise) rawLatData structure is converted into a (not precise) 'double', which is indeed only a 'float' with Atmega 8-bit controllers.

So the internal storage precision of locations in a special data structure of TinyGPSplus will be absolutely WORTHLESS with any calculations and data retrievel you might do using that library. So if you want to do any "precision data handling" with tinyGPSplus, you would have to access the "rawLat()" and "rawLng()" high precision data of that library directly and use them with your own high precision calculations.

So if I wanted to do high precision calculations of small GPS distances, I'd completely avoid tinyGPSplus and use some own parsing routines to parse the received NMEA data directly into 'long' variables with 9-10 digits precision like explained in reply #7. Then use an optimized haversine formula to calculate a small and precise distance.

Even with 8-bit Atmegas and only 'long' and 'float' calculations it should be possible to find an algorithm to calculate small GPS distances accurate up to less than 10 millimeters (0,01m).

Please let me know if you need some assistance in writing code.

Please let me know if you need some assistance in writing code.

Yes, it would be excellent if you would post a library where the trig and common navigational functions are recoded to take full advantage of the expected precision and accuracy of your scheme.

jremington:
Yes, it would be excellent if you would post a library where the trig and common navigational functions are recoded to take full advantage of the expected precision and accuracy of your scheme.

What "common navigational functions" are you thinking about?

The thread starter thought about "distance to" and "bearing to" functions.

What else would be needed to calculate?

As a feasibility check for my suggestin I provide some test code with two functions:

  • nmeaDegrees2Decimal()
  • calcDistanceHalversine()

The nmeaDegrees2Decimal() function takes a string with NMEA encoded degrees and converts to a decimal "long" value as described in reply #7.

The calcDistanceHalversine() function takes such "long" coordinates and calculates the distance in kilometers.

The example shows a result of 0.000007 km, which is 0.007m or 7 millimeters distance when calculated on an 8-bit Atmega328.

When calculating the distance on a PC with Javascript on page Calculate distance and bearing between two Latitude/Longitude points using haversine formula in JavaScript between
Point-1: 50.5000000 / 10.2057612
Point-2: 50.5000000 / 10.2057613
then I cannot find a big difference. The results with Javascript on a 64-bit Windows PC seem to be pretty the near to the results on a 8-bit Atmega328.

Or what do you think?

This is the test code for NMEA coordinate conversion and for distance calculation between two locations:

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 calcDistanceHalversine(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;
}


void setup()
{
  Serial.begin(9600);
  Serial.println();
  long lat1, lon1, lat2, lon2;
  nmeaDegrees2Decimal("5030.00000", lat1 ); // 50 degrees 30.0 arc minutes
  Serial.print("lat1= ");Serial.println(lat1);
  
  nmeaDegrees2Decimal("1012.34567", lon1 ); // 10 degrees 12.34567 arc minutes
  Serial.print("lon1= ");Serial.println(lon1);
  
  lat2=lat1;
  Serial.print("lat2= ");Serial.println(lat2);
  
  nmeaDegrees2Decimal("1012.34568", lon2 ); // 10 degrees 12.34568 arc minutes
  Serial.print("lon2= ");Serial.println(lon2);
  Serial.println();
  Serial.print("Distance [km] = ");
  double distance=calcDistanceHalversine(lat1, lon1, lat2, lon2);
  Serial.println(distance,6);
  Serial.print("Distance [m] = ");
  Serial.println(distance*1000,6);
}

void loop()
{
}

Certainly good enough for amateur robotics!

Great circle calculations won't work as well as the Haversine function, due to the inaccuracy of the single precision trigonometric functions.

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.

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()
{
}

BTW: Does anybody know where the earth radius of 6372795 meters = 6372.795 km was taken from, which is used in calculations of tinyGPSplus?

In my calculations I'm using an earth radius of 6371 km, which is taken from Wikipedia and can be calculated by different ways to be the mean radius of earth (including mean radius of the WGS84 ellipsoid).

I am going to guess they use some standard spheroid instead of the actual value. In 'Big Red' It lists the value as 6378388 meters, and says it is based on the "International Hayford Spheroid".

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.

See reply #2 for an example. In fact, it works much, much better.

Try this one for exact GPS calculations :frowning: TinyGPS is great for getting GPS data but the calculations are useless for distances < 100m )

and as 64bit float calculation are VERY CPU expensive you may actually need to do the calculations is background thread as long calculation will prevent you reading from devices on time..

e.g. you need also :

and if you are still brave enough take a look at real-world application example :

vtomanov:
and as 64bit float calculation are VERY CPU expensive you may actually need to do the calculations is background thread as long calculation will prevent you reading from devices on time..

That is because used Math64 lib is not very efficient. It use also Taylor series known converges very slowly, tan is calculated by sin/cos, sqrt with exp function, etc. You have used your own replacements, however that is still too much multiplications and divisions...

The famous HP-35 from 1972 used a custom made chip to process BCD numbers with execution time of just 250us for one clock per instruction, but still was able to process sqrt and trig functions (sin, cos, tan and natural log including its inversions) with 10 digits precision in max ~500ms. That was really impressive for the time, knowing all was accomplished with addition and shifting with in only 768 "words" of code, with the chip had only 5 registers!

Imagine what can be done with 16MHz clock and modern (just) 8-bit architecture CPU having multiplication instruction, at least 32K ROM, 2K RAM and 1K EEPROM...

I'm working for a while on fast double precision (to the last bit) lib for 8-bits based MCUs as I have a need as well for max precision for my own not trivial at all GPS navigator project (including maps etc). Really there is no need for 32-bit MCU platform.

FWIW, I'm not sure 64-bit math is required for GPS... NeoGPS achieves 10 significant digits, about 1cm (interesting discussion here). This is better than almost any GPS device can produce, and much better than the common TinyGPS, TinyGPS+ and Adafruit_GPS libraries. Furthermore, the NeoGPS distance/bearing calculations are structured to retain this accuracy, so 64-bit math isn't required for those operations, either.

Remember that everything is based on an approximation to the Earth's surface. Not that 64-bit math isn't interesting... :wink:

Cheers,
/dev

It is true, but Haversine and Vincenty formula cannot be compared in accuracy.

As well, from public GPS transponder now send signal which allow precision under 2.5m, but that was more than 11 in the past. I'm not sure what is current policy for GLONASS and Galileo, but it is possible as well that USA military allow for GPS maximum precision need for autonomous vehicles...

Secondly, used float trig functions have mostly precision up to 6 digits, but be aware that every calculation degrade precision of result moderately in some complex trig. formulas or any plain polynomial expression as some Taylor series.