double precision / TinyGPS++

[phil_rudd69]

Hi,

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. I generate the double value from these entries with the -'0' and the pow(10,n) function by multipicating the seperate entries with the correct power n and then adding them up to the correct number:

Example: 1*10^1 + 2*10^0 + 1*10^(-1) + 2*10^-2) + ... = 12.12346578

The problem is that I loose precision after the 6th decimal place. I know this calculation is a bit unhandy but I already tried atof()-commands and they all failed as well.

I tried the following to get to the problem: I added the following numbers as double-values

   0.000003
+ 0.0000003
+ 0.00000003
= 0.00000333

which was calculated correctly. But as I added numbers like

   1.000003
+ 0.0000003
+ 0.00000003
= 1.00000331675

I got crap like this. Obviously because of the additional value left of the decimal point.

Any ideas on how to fix this problem?

Note: I already googled the problem and read, that MEGA boards could only handle 6 decimal places or less.

Nonetheless I wanted to ask if there are any workarounds existing?
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 
Unfortunately I don't understand what is going on within the functions of TinyGPS++. 

My Project:
I want my robot to drive from a point A to point B via GPS. For that I submit GPS-destination-coordinates like "12.12345678" via Serial Monitor to my Arduino Mega 2560. I use the TinyGPS++ library to get information like "distance to destination B" and "current orientation of robot". I want my robot to find the destination at least within a 5 meter-radius. But this precision needs more decimal places then 6.

[PaulS]

but I already tried atof()-commands and they all failed as well.

No. The failure is on your part. A double on the Arduino is the same size, and precision, as a float - 4 bytes. A 4 byte float has a precision of 6 or 7 decimal places. Expecting better than that is a mistake.

For that I submit GPS-destination-coordinates like "12.12345678" via Serial Monitor to my Arduino Mega 2560.

That value is a combination of degrees, minutes, and seconds. You could parse that string (as a string, not as a double) to get degrees, minutes, and seconds. You could parse the GPS string (as a string, not as a double) to get degrees, minutes, and seconds. You can not expect seconds to 5 decimal places, though. You could then write your own functions for getting distance based on degrees, minutes, and seconds.

[jremington]

If you store GPS coordinates as degrees*1000000 in a long integer, you can preserve the accuracy needed for global positioning. However, calculations using trigonometric functions are limited to single precision floating point on the Arduino.

It is possible to get around this in various ways. For example, you can calculate the distance and bearing between two (lat, long) points using the equirectangular approximation, which treats the Earth locally as a flat surface. This is good for points up to about 10 km distant, with precision and accuracy of better than 1 meter (better than GPS accuracy). See the code below for an excellent web site that describes the various approximations.

Example code:

Code: [Select]


// find the bearing and distance from point 1 to 2, using the equirectangular approximation
// lat and lon are degrees*1.0e6
// bearing is returned in degrees, distance in meters

// reference: http://www.movable-type.co.uk/scripts/latlong.html

double course_to(long lat1, long lon1, long lat2, long lon2, double* distance) {

double dlam,dphi,radius;

double deg2rad = 3.141592654/180.0;
double rad2deg = 180.0/3.141592654;
radius = 6371000.0; //approximate Earth radius in meters.

dphi = deg2rad*(lat1+lat2)*0.5e-6; //average latitude in radians
double cphi=cos(dphi);

dphi = deg2rad*(lat2-lat1)*1.0e-6; //differences in degrees
dlam = deg2rad*(lon2-lon1)*1.0e-6;

dlam *= cphi;  //correct for latitude

double bearing=rad2deg*atan2(dlam,dphi);
if(bearing<0) bearing=bearing+360.0;

*distance = radius * sqrt(dphi*dphi + dlam*dlam);

return bearing;
}


[robtillaart]

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.

[jurs]

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.

[jremington]

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.

[pito]

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

[jurs]

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.

[pito]

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

Code: [Select]


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

[odometer]

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.

[jurs]

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

Code: [Select]


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

Code: [Select]


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.

[jremington]

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.

[jurs]

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 http://www.movable-type.co.uk/scripts/latlong.html 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:

Code: [Select]


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

[jremington]

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.

[jurs]

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.

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:

Code: [Select]


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