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 support@ccsinfo.com
Compute distance using Lat & long data

 
Post new topic   Reply to topic    CCS Forum Index -> General CCS C Discussion
View previous topic :: View next topic  
Author Message
minoan_war
Guest
Compute distance using Lat & long data
PostPosted: Wed Aug 09, 2006 2:33 pm     Reply with quote
HI All!
I work in a small program to get the distance between 2 points using latitud and longitud data,
this is the program:

#define EARTH 6371
#define PI 3.141592654
float d_lat1=19.508016;
float d_lat2=19.674716;
float d_lon1=99.234144;
float d_lon2=99.234000;
float distance;

main()
{
clrscr(); // turbo c
// degrees to radians!
d_lat1=d_lat1*PI/180;
d_lat2=d_lat2*PI/180;
d_lon1=d_lon1*PI/180;
d_lon2=d_lon2*PI/180;

distance=acos( (sin(d_lat1)*sin(d_lat2)) + (cos(d_lat1)*cos(d_lat2)*cos(d_lon2-d_lon1)))*EARTH ;
printf("\n\nDistance is %f",distance);
getch();
}

Where EARTH is the Curvature of the EARTH

using Turbo c the result is OK! 18.5 km
using Hi tech compiler the result is 17.8 Km
using CCS compiler the result is 15.2 Km

I can live with the result of CCS compiler, but... when the distance between points is more small the error math is bigger

Anyone has idea Why this discrepances?

Thanks
Douglas Kennedy



Joined: 07 Sep 2003
Posts: 755
Location: Florida

View user's profile Send private message AIM Address

PostPosted: Wed Aug 09, 2006 3:39 pm     Reply with quote
This is probably a precision issue with a 23 bit mantissa in floating point binary with CCS. Take a look at the CORDIC algorithm it can get you as much precision as you want. The CORDIC got the Lunar lander down and it wasn't out 2 or 3 kilometers.
Below is a 13bit example using degrees
you can get the sine cosine arctan and sqrt(x^2+y^2) with nothing more than addition. Watch out for crossing quadrants.
Code:

#include <18F452.h>

#DEVICE *=16 ICD=TRUE




#fuses HS,NOWDT,NOPROTECT,NOPUT,NOLVP,NOBROWNOUT,NOWRT


#use delay(clock=40 000 000)
//#use rs232(baud=57600, bits=8, parity = N, xmit=PIN_C6, rcv=PIN_C7)
#use rs232(DEBUGGER)

////////////        cordic   ///////////////////////////////////
////////////////////////////////////////////////////////////////
///               theory                                   /////
///                                                        /////
///    coords (x,y) are rotated thru angle a to (x',y')    /////
///                                                        /////
///           x'=x.cos(a)-y.sin(a)                         /////
///           y'=y.cos(a)+x.sin(a)                         /////
///                                                        /////
///           x'/cos(a)=x-y.tan(a)                         /////
///           y'/cos(a)=y+x.tan(a)                         /////
///                                                        /////
///           now if we choose tan(a) to be                /////
///            1/2,1/4,1/8......1/2*n                      /////
///                                                        /////
///           we can rotate thru sucessive values          /////
///           ex 30.00 degs = 45-26.57+14.03-7.13+3.58     ////
///                           +17.9,.90,.45                ////
///          where a minus angle is a counter rotation     ////
///          after a rotation of 30 degs                   ////
///          x' is x.cos(30) if we start with y=0          ////
///          y' is x.sin(30)                               ////
///                                                        ////
///          angle is in 1/100th units 4500 is 45 degrees  ////
///          lengths are in aggregate constant units       ////
///                                                        ////
///  aggregate const cos(atan(1)*cos(atan(1/2)*cos(atan(1/4).....)))))
///  or .607252935                                         ////
///  so 10000 is 6073

long int const atan[13]={4500,2656,1404,712,358,179,90,45,
                 22,12,6,3,1};/// atan(1/2n)
////
////
//// to clac sin(A) cos(A) set y=0 a=angle and drive a to zero
//// via iterations
//// x value is the cos(a) providing x=6073
//// y value is the sign(a) providing x=6073
////
//// to calc atan(a) set a=0 and drive y to zero
//// a value is atan(y/x)
//// x value is sqrt(x^2+y^2) providing x=6073

signed long int cos(signed long a)
{
int i;
signed long  x,y,dx,dy,da;;
x=6073;
y=0;
for (i=0;i<13>=0) {dx=X >>i;} else { dx=-X >> i;dx=-dx;}
if (Y>=0) {dy=Y >>i;} else { dy=-Y >> i;dy=-dy;}

  da=atan[i];
  if (a>=0)
          {
          a=a-da;
          x=x-dy;
          y=y+dx;
          }
 else     {
          a=a+da;
          x=x+dy;
          y=y-dx;
          }
 }
return (x);
}
signed long int sin(signed long a)
{
int i;
signed long  x,y,dx,dy,da;;
x=6073;
y=0;
for (i=0;i<13>=0) {dx=X >>i;} else { dx=-X >> i;dx=-dx;}
if (Y>=0) {dy=Y >>i;} else { dy=-Y >> i;dy=-dy;}

  da=atan[i];
  if (a>=0)
          {
          a=a-da;
          x=x-dy;
          y=y+dx;
          }
 else     {
          a=a+da;
          x=x+dy;
          y=y-dx;
          }
 }
return (y);
}


signed long int arctan(signed long x,signed long y)
{
int i;
signed long  a,dx,dy,da;;

a=0;
for (i=0;i<13>=0) {dx=X >>i;} else { dx=-X >> i;dx=-dx;}
if (Y>=0) {dy=Y >>i;} else { dy=-Y >> i;dy=-dy;}

  da=atan[i];
  if (y<0)
          {
          a=a-da;
          x=x-dy;
          y=y+dx;
          }
 else     {
          a=a+da;
          x=x+dy;
          y=y-dx;
          }
 }
return (a);
}

signed long int sqrt_x2_y2(signed long x,signed long y)
{
int i;
signed long  a,dx,dy,da;;

a=0;
for (i=0;i<13>=0) {dx=X >>i;} else { dx=-X >> i;dx=-dx;}
if (Y>=0) {dy=Y >>i;} else { dy=-Y >> i;dy=-dy;}

  da=atan[i];
  if (y<0)
          {
          a=a-da;
          x=x-dy;
          y=y+dx;
          }
 else     {
          a=a+da;
          x=x+dy;
          y=y-dx;
          }
 }
return (x);
}


main()
{
signed int16 result;

result= cos(6000);

printf("cosine of 60 is %ld \n\r",result);

result= sin(6000);

printf("sine of 60 is %ld \n\r",result);

result= arctan(10000,6000);

printf("arctan of 0.60 is %ld \n\r",result);

result=sqrt_x2_y2(10000,6000);

printf("sqrt of (6^2+10^2)/0.6073 is %ld \n\r",result);
while(true);

}


Last edited by Douglas Kennedy on Wed Aug 09, 2006 3:46 pm; edited 1 time in total
Ttelmah
Guest
PostPosted: Wed Aug 09, 2006 3:42 pm     Reply with quote
The big problem is rounding errors. For instance, your 'sin' value, is being held as 3.1415927, and this is one of the least 'rounded' of the values involved... Turbo C, will probably be using a 'double', internally for the sin calculations.
Your 'Earth' diameter is slightly below the normally accepted average value (6372.795), but the effect of this is small.
Generally, you might be better off using the haversine formula, which shows less rounding errors for small angles. You would need to add testing for the denominator being zero in the final atan function for this.

Best Wishes
minoan_war
Guest
Haversine implementation
PostPosted: Wed Aug 09, 2006 3:59 pm     Reply with quote
Thanks Douglas, Ttelmah!
I have done the Haversine implementation, ands look a little better

this is!

tmp1=d_lat2-d_lat1;
tmp2=d_lon2-d_lon1;

tmp3=(sin(tmp1/2) * sin(tmp1/2)) + cos(d_lat1) * cos(d_lat2) * (sin(tmp2/2) * sin(tmp2/2));
tmp4=2 * atan2(sqrt(tmp3),sqrt(1-tmp3));
tmp5=EARTH * tmp4;
printf("\nHaversine %f Km",tmp5);

prints "Haversine 18.8 Km"

it have also a error of .4 km, I will test it, extensively
I think that for a $3USD chip is ok!
ryan.reeve



Joined: 23 Jul 2006
Posts: 20

View user's profile Send private message

whats wrong in haversine formula ?
PostPosted: Thu Sep 28, 2006 4:30 am     Reply with quote
Hi
I am using the Haversine formula, but value is too much erronious.Not even near the expected value.

float lat1,lat2,dlat,long1,long2,dlong,tmp1,tmp2,d;
lat1=31.000100;
lat2=31.000000;
long1=74.260000;
long2=74.260000;
dlat=lat1-lat2;
dlong=long1-long2;

tmp1=(sin(dlat/2)*sin(dlat/2))+cos(lat1)*cos(lat2)*(sin(dlong/2*sin(dlong/2)));
tmp2=2*atan2(sqrt(tmp1),sqrt(1-tmp1));
d=EARTH*tmp2;
Ttelmah
Guest
PostPosted: Sun Oct 01, 2006 5:04 am     Reply with quote
Remember that in C, sin, cos etc., work in _radians_, not degrees....

Best Wishes
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