|
|
View previous topic :: View next topic |
Author |
Message |
vafonkin
Joined: 20 Mar 2012 Posts: 16
|
Problems with math calculations |
Posted: Mon Apr 16, 2012 5:12 am |
|
|
Can anyone tell me why following code does not calculate properly? Result should be 711, when CCS gives 938. MPLab C30, when using long double gives correct result for the same code:
Code: |
#define TO_RAD 3.141592653589793/180.0
#define TO_DEG (180.0 / 3.141592653589793)
#define R 6371.0 // km
float dist(float distlat1, float distlon1, float distlat2, float distlon2);
void main (void){
float fDist;
while(1){
fDist=dist(42.307524, -71.201024, 42.301524, -71.198024);
printf("%f\n",fDist);
}
}
float dist(float distlat1, float distlon1, float distlat2, float distlon2){
float d;
float slat1,slat2,clat1,clat2,cl2l1;
slat1=sin(distlat1*TO_RAD);
slat2=sin(distlat2*TO_RAD);
clat1=cos(distlat1*TO_RAD);
clat2=cos(distlat2*TO_RAD);
cl2l1=cos((distlon2-distlon1)*TO_RAD);
d=acos(slat1*slat2+clat1*clat2*cl2l1)*R;
return(d);
} |
|
|
|
SherpaDoug
Joined: 07 Sep 2003 Posts: 1640 Location: Cape Cod Mass USA
|
|
Posted: Mon Apr 16, 2012 6:01 am |
|
|
A CCS float is only good to about 6 digits.
What you are trying to do is pretty common, though I have never dealt with it. I expect if you search through the code library you will find the answers you need. _________________ The search for better is endless. Instead simply find very good and get the job done. |
|
|
ckielstra
Joined: 18 Mar 2004 Posts: 3680 Location: The Netherlands
|
|
Posted: Mon Apr 16, 2012 6:09 am |
|
|
Always post your compiler version number! |
|
|
Ttelmah
Joined: 11 Mar 2010 Posts: 19496
|
|
Posted: Mon Apr 16, 2012 7:29 am |
|
|
and chip.
The float on a PIC16/18, with CCS, only supports 'single precision', but on a PIC24 for example, double is supported, but the function would have to be declared with 'double', not 'float'....
Best Wishes |
|
|
RF_Developer
Joined: 07 Feb 2011 Posts: 839
|
|
Posted: Tue Apr 17, 2012 2:56 am |
|
|
We've been over this ground before in a previous thread.
You are using a computationally naive process: one that works fine mathematically but doesn't work well under some practical computational conditions.
There are three ways to calculate this, one is the spherical rule of cosines that you are using. This is computationally tricky when the points are close together, as they are in your case, which appears to be a golf course in Boston. You must work through the numbers carefully to see what's happening, but essentially the difference of the longitudes is close to zero which makes its cosine close to but not quite one, that is messing up the result a lot.
The second mathematically equivalent process is the haversine formula. This was used in practical navigation long before the rule of cosines as its well suited to limited precision of manual calculation using tables. By using mathematical identities based on half the angles it maintains higher precision by work in the high slope region of sines and cosines. It is still a better bet for limited computation precision for processors. In fact its a much better way to calculate this than the rule of cosines if double/extended precision is not available, as on most PICs. It does involve more calculations than the rule of cosines so its slower, but at least it should work, which is a BIG advantage.
The third process is an approximation that works satisfactorily over small ranges. I suggested this in the last thread, that is to project the points on to a surface, or do so virtually and then use simple planar trigonometry to get the distance. This is what a navigator does on a chart. Over the distances you want, a few kilometers or so, it should work fine.
I looked at this webpage; http://www.movable-type.co.uk/scripts/latlong.html which has a online calculator and java code for all three processes. I put them into a C# app and tried them all with float precision. I have to say that C#'s maths routines are all double precision only and so I had to cast in and out. That does not produce the same results as doing everything in float. Even if it did, I am not certain that CCS float is the same precision as IEEE754 single length float. It probably is, but I don't know that for sure.
The results were illuminating. For your lat/longs the rule of cosines would not produce a result at all as the input to the arc cosine was greater than one, mainly as the cosine of the difference of the longitudes computed to 1.0, where it should be slightly less than one. The precision was such that it was rounded up to one. Note such rounding may well differ on the fully IEEE 754 PC to the non 754 CCS routines. This sort of subtlety is often discarded to improve the performance. The PIC cosine may well not round, and simply truncate, giving a smaller, but still wrong, result.
The haversine version gave a result of 711.549m and looked solid on the potentially trick case of the lats being the same (i.e. an east-west line).
The projected version also gave 711.549m. So over this range the errors caused by the approximation is probably smaller than the precision error and so is small enough to be ignored. As it uses no trig, just simple Pythagoras, its going to be a lot faster on the PIC.
This is all simple stuff. All I did was Google "distance on the earth" and look at the second hit. Then I did a simple test app in the PC. Easy.
As I said before you have to be careful with special cases: will it work with north-south and east-west lines? Does it matter which order the points are given? Will it work for all lats/longs? Try it at and around the poles for example. What about the trivially simple case of both points being the same? Is the distance zero? The trig version tend to have trouble with that one. What happens as the distance tends to zero? All important stuff if you are doing surveying and map-making.
This is not straightforward. The maths might be easy enough, but the computational precision implications were only taught to me at degree level. Until then I blindly and naively assumed the maths and the computations would always be the same, despite a fair bit of practical evidence to the contrary. Now I go through all such code line by line to make sure it really is doing what I expect it to do.
Here's the guts of my C# test app. Its not compilable as is and is for C# not C.
Code: |
const float TO_RAD = 3.141592653589793f/180.0f;
const float TO_DEG = 180.0f / 3.141592653589793f;
const float R = 6371.0f; // km
float dist(float distlat1, float distlon1, float distlat2, float distlon2)
{
float d;
// Rule of Cosines
float slat1, slat2, clat1, clat2, cl2l1;
slat1 = (float)Math.Sin(distlat1 * TO_RAD);
slat2 = (float)Math.Sin(distlat2 * TO_RAD);
clat1 = (float)Math.Cos(distlat1 * TO_RAD);
clat2 = (float)Math.Cos(distlat2 * TO_RAD);
cl2l1 = (float)Math.Cos((distlon2 - distlon1) * TO_RAD);
d = (float)Math.Acos(slat1 * slat2 + clat1 * clat2 * cl2l1) * R;
// Haversines
float dLat, dLon, lat1, lat2, a, c;
dLat = (distlat2 - distlat1) * TO_RAD;
dLon = (distlon2 - distlon1) * TO_RAD;
lat1 = distlat1 * TO_RAD;
lat2 = distlat2 * TO_RAD;
a = (float)(Math.Sin(dLat / 2) * Math.Sin(dLat / 2) +
Math.Sin(dLon / 2) * Math.Sin(dLon / 2) * Math.Cos(lat1) * Math.Cos(lat2));
c = 2.0f * (float)Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
d = R * c;
// Project and Pythag
float x, y;
x = (distlon2 - distlon1) * TO_RAD * (float)Math.Cos((distlat1 + distlat2) * TO_RAD / 2);
y = (distlat2 - distlat1) * TO_RAD;
d = (float)(Math.Sqrt(x * x + y * y) * R);
return d;
}
void Calculate_Distance (){
float fDist;
fDist = dist((float)Convert.ToDouble(Lat_A.Value), (float)Convert.ToDouble(Long_A.Value), (float)Convert.ToDouble(Lat_B.Value), (float)Convert.ToDouble(Long_B.Value));
//fDist=dist(42.307524f, -71.201024f, 42.301524f, -71.198024f);
Float_Result.Text = fDist.ToString();
}
|
RF Developer |
|
|
SherpaDoug
Joined: 07 Sep 2003 Posts: 1640 Location: Cape Cod Mass USA
|
|
Posted: Tue Apr 17, 2012 6:22 am |
|
|
RF-Developer, that is an excellent reply! You tell him/us what is wrong, why it is wrong, how to fix it for high precision and for adequate precision for what seems to be the goal, and how he could have found the solution himself.
That is going into my personal software archive. I have never had to solve spherical geometry problems, but someday I may. _________________ The search for better is endless. Instead simply find very good and get the job done. |
|
|
miro
Joined: 15 Jan 2011 Posts: 62
|
|
Posted: Tue Apr 17, 2012 2:03 pm |
|
|
Be careful with CCS float math calcs. The only calcs which work with a precision comparable with C30 are * / + - (both float and double). Do test trig calcs before using them! |
|
|
|
|
You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|