CCS C Software and Maintenance Offers
FAQFAQ   FAQForum Help   FAQOfficial CCS Support   SearchSearch  RegisterRegister 

ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

CCS does not monitor this forum on a regular basis.

Please do not post bug reports on this forum. Send them to CCS Technical Support

Problems with math calculations

 
Post new topic   Reply to topic    CCS Forum Index -> General CCS C Discussion
View previous topic :: View next topic  
Author Message
vafonkin



Joined: 20 Mar 2012
Posts: 16

View user's profile Send private message

Problems with math calculations
PostPosted: Mon Apr 16, 2012 5:12 am     Reply with quote

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

View user's profile Send private message

PostPosted: Mon Apr 16, 2012 6:01 am     Reply with quote

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

View user's profile Send private message

PostPosted: Mon Apr 16, 2012 6:09 am     Reply with quote

Always post your compiler version number!
Ttelmah



Joined: 11 Mar 2010
Posts: 19539

View user's profile Send private message

PostPosted: Mon Apr 16, 2012 7:29 am     Reply with quote

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

View user's profile Send private message

PostPosted: Tue Apr 17, 2012 2:56 am     Reply with quote

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

View user's profile Send private message

PostPosted: Tue Apr 17, 2012 6:22 am     Reply with quote

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

View user's profile Send private message

PostPosted: Tue Apr 17, 2012 2:03 pm     Reply with quote

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!
Display posts from previous:   
Post new topic   Reply to topic    CCS Forum Index -> General CCS C Discussion All times are GMT - 6 Hours
Page 1 of 1

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