View previous topic :: View next topic |
Author |
Message |
vafonkin
Joined: 20 Mar 2012 Posts: 16
|
Intersection of two paths given start points and bearings? |
Posted: Tue Mar 20, 2012 8:08 pm |
|
|
I spent about 2 weeks trying different approaches to get calculated distance between 2 GPS coordinates, intersection of two paths given start points and bearings, and calculate bearing to the giving point from current location and so far was able to calculate distance using "Flat Earth" approach by converting GPS coordinates to Cartesian and calculate distance in Cartesian.
What I was NOT able to complete is getting intersection and bearing.
I tried cordic library - this sounds like a valid approach - but even with that not able to get results.
My project does not require spherical earth - since all points a located within maximum 4km from each other and intersection point can be anywhere withing 4km by 4 km area.
Please advise if anyone was able successfully calculate intersection and bearing for Flat Earth model with any PIC controller, or this just a dead end road.
Thank you in advance. |
|
|
RF_Developer
Joined: 07 Feb 2011 Posts: 839
|
|
Posted: Wed Mar 21, 2012 2:15 am |
|
|
Flat, spherical or even spheroidal earth should be no problem. Its straightforward and well documented trigonometry, which the PIC is quite capable of, albeit sluggishly slowly. The CCS C float precision should be enough for what you want to do.
So, where is your problem? Is this a CCS C specific problem - in which case post what you've got - or is it that you simply don't know how to do the maths - in which case use Google.
RF Developer |
|
|
Mike Walne
Joined: 19 Feb 2004 Posts: 1785 Location: Boston Spa UK
|
|
Posted: Wed Mar 21, 2012 2:23 am |
|
|
What is the required accuracy?
You may be able to get away with mainly integer maths.
Like RF developer says, you've not explained what your real problem is.
Mike |
|
|
vafonkin
Joined: 20 Mar 2012 Posts: 16
|
|
Posted: Wed Mar 21, 2012 6:15 am |
|
|
10 meters for required accuracy would be fine. I cannot compile/get valid results from following foxLocation function. It implements Haversine and PIC MCU appears cannot handle it. I tried to split code to smaller chuncks - does not help:
Code: |
#define R 6371
#define TO_RAD (3.1415926 / 180)
#define TO_DEG (180 / 3.1415926)
float foxLat,foxLon;
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
//:: Function to calculate cross bearing :
//:: Takes first bearing location lat and lon and true bearing :
//:: second bearing location lat and lon and true bearing :
//:: :
//:: Assignes values to cross bearing location foxLat and foxLon :
//:: Example: :
//:: foxLocation(42.307524,-71.201024,260,42.304524,-71.203024,350); :
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/*
void foxLocation (float lat1,float lon1,int16 brng1,float lat2,float lon2,int16 brng2) {
float dLat=0;
float dLon=0;
int16 brng12=0;
int16 brng21=0;
int16 dist12=0;
int16 brngA=0;
int16 brngB=0;
float alpha1=0;
float alpha2=0;
float alpha3=0;
int16 dist13=0;
float dLon13=0;
lat1 = TO_RAD*lat1;
lon1 = TO_RAD*lon1;
lat2 = TO_RAD*lat2;
lon2 = TO_RAD*lon2;
brng1 = TO_RAD*(brng1/10);
brng2 = TO_RAD*(brng2/10);
dist12=(lat1,lon1,lat2,lon2);
// initial/final bearings between points
brngA = acos((sin(lat2) - sin(lat1)*cos(dist12)) / (sin(dist12)*cos(lat1)));
brngB = acos((sin(lat1) - sin(lat2)*cos(dist12)) / (sin(dist12)*cos(lat2)));
if (sin(lon2-lon1) > 0) {
brng12 = brngA;
brng21 = 2*3.1415926 - brngB;
} else {
brng12 = 2*3.1415926 - brngA;
brng21 = brngB;
}
alpha1 = fmod((brng1 - brng12 + 3.1415926),(2*3.1415926)) - 3.1415926; // angle 2-1-3
alpha2 = fmod((brng21 - brng2 + 3.1415926),(2*3.1415926)) - 3.1415926; // angle 1-2-3
alpha3 = acos(-cos(alpha1)*cos(alpha2) + sin(alpha1)*sin(alpha2)*cos(dist12));
dist13 = atan2(sin(dist12)*sin(alpha1)*sin(alpha2),cos(alpha2)+cos(alpha1)*cos(alpha3));
foxLat = asin(sin(lat1)*cos(dist13) + cos(lat1)*sin(dist13)*cos(brng1) );
dLon13 = atan2(sin(brng1)*sin(dist13)*cos(lat1),cos(dist13)-sin(lat1)*sin(foxLat));
foxLon = lon1+dLon13;
foxLon = fmod((foxLon+3.1415926),(2*3.1415926)) - 3.1415926; // normalise lon3 to -180..180º
foxLon = TO_DEG*foxLon;
foxLat = TO_DEG*foxLat;
}
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
//:: Calculate distance between 2 points, assume Earthis flat :
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
int16 dist(float distlat1, float distlon1, float distlat2, float distlon2){
signed int32 x1,y1,x2,y2;
float cc1,cc2;
int16 d;
// convert to Cartesian coordinates
cc1=cos(TO_RAD*distlat1);
x1=1852*cc1;
x1=x1*60;
x1=x1*distlon1;
y1=distlat1*60*1852;
cc2=cos(TO_RAD*distlat2);
x2=1852*cc2;
x2=x2*60;
x2=x2*distlon2;
y2=distlat2*60*1852;
// calculate distance
d=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
return(d);
}
|
|
|
|
RF_Developer
Joined: 07 Feb 2011 Posts: 839
|
|
Posted: Wed Mar 21, 2012 7:26 am |
|
|
It compiled OK for me on PCWH 4.130. Lots of code of course, and it would run like a dog, but it compiled OK.
I did a lot of this sort of thing with a radar project I was involved with in the 90's. There are many approaches, but I found that once I'd sorted out the spherical geometry (which you can get from the net) it was actually simpler to do all the math in spherical form. I had object methods to give me range and direction from point to point, and to project the lat long on to maps in various projections.
This points to the other relatively simple way to do this: project the points on to a vritual mercator map and then, once you've sussed the scale at a suitable point (and assumed it to be valid for the entire region of interest), simply do the crossing point math in straightforward planar geometry. That should be accurate enough for you, though doing it in spherical geometry would give you a better result (I was working over several tens of km). This approach is the virtual equivalent of doing it by hand on a paper navigational chart.
Anyway, it compiles for me. I cant say whether the results are right, and I know from experience that there are always all sorts of special cases to consider to make this sort of thing work over all possible combinations of input points.
RF Developer
Anyway, |
|
|
vafonkin
Joined: 20 Mar 2012 Posts: 16
|
|
Posted: Wed Mar 21, 2012 7:40 am |
|
|
Which PIC MCU processor do you use? |
|
|
RF_Developer
Joined: 07 Feb 2011 Posts: 839
|
|
Posted: Wed Mar 21, 2012 8:34 am |
|
|
No, which do YOU use? For the record, I compile the code for PIC18F4520 using PCWH 4.130. 34% of ROM was used, after all its full-on floating point, what do you expect?
RF Developer |
|
|
vafonkin
Joined: 20 Mar 2012 Posts: 16
|
|
Posted: Wed Mar 21, 2012 8:49 am |
|
|
16F886 ... obviously not enough of power.
Will move to new processor and test it out there.
Thanks a lot!
73!
KB1RLI |
|
|
dorinm
Joined: 07 Jan 2006 Posts: 38
|
|
Posted: Wed Mar 21, 2012 12:57 pm |
|
|
try converting your coordinates in UTM, then it's straightforward and easy even for a 12f to do the mat;
UTM coordinates are in cm, so your 10km would be expressed in 1.000.000cm domain, always positive, so you would use unsigned32 or signed, either one won't overflow and produce correct result. |
|
|
vafonkin
Joined: 20 Mar 2012 Posts: 16
|
|
Posted: Wed Mar 21, 2012 2:28 pm |
|
|
dorinm wrote: | try converting your coordinates in UTM, then it's straightforward and easy even for a 12f to do the mat;
UTM coordinates are in cm, so your 10km would be expressed in 1.000.000cm domain, always positive, so you would use unsigned32 or signed, either one won't overflow and produce correct result. |
Thank you for advice. I tried it already - using cordic.h for functions and signed int32 for calculations - not enough RAM. It simply will not fit into 368 bytes available. In my app I am reading data from compass over I2C, reading and converting GPS NMEA sentence over rs232, show results on LCD over another rs232... storing compass bearing, do some calculations...but everything works until I start cross bearing or bearing calculations where a lot of trig functions used. Those are too heavy for 16F family. I placed an order for PIC24FJ64GA002 and PIC18LF2520 - both has plenty of memory and support 8x8 or 17x17 hardware multiplication. I would have to redesign a PCB, but if it works I will be happy. |
|
|
vafonkin
Joined: 20 Mar 2012 Posts: 16
|
|
Posted: Tue Mar 27, 2012 6:01 am |
|
|
It appears triginimetric functions provide incorrect results starting with 4th decimal point. It cannot be used for navigation formula due to lack of precision.
.
I need to calculate intersection point given 2 gps locations and bearings, and I need to get bearing from current location to given point.
Area for calculation is small 4 by 4 km, precision should be not greater then 10 meters.
I am sure someone already done it and I will appreciate any help
If free help is not available I am willing to pay for solution. |
|
|
RF_Developer
Joined: 07 Feb 2011 Posts: 839
|
|
Posted: Tue Mar 27, 2012 6:13 am |
|
|
4th decimal place huh? In other words 6 or 7 significant places. Yes, that's about right for single length float precision, and CCS C doesn't have a double type. I don't know if other PIC Cs have a double type. If so, you'll have to use them. Either that or implement your own double length floating point... remembering that famously even Intel got that wrong on the first pentiums.
RF Developer. |
|
|
vafonkin
Joined: 20 Mar 2012 Posts: 16
|
|
Posted: Tue Mar 27, 2012 6:39 am |
|
|
I believe there are some approaches around:
1. Cordic
2. Lookup tables
3. Some other high math functions
It is probably possible to find intersection by solving simple line formula:
ax + by + c = 0
y = kx + d, где k=-(a/b), d=-(c/b)
k=tan(a)
d=y-tan(a)*x
b=1
k1x + d1 = y
k2x + d2 = y
I done it, however tan function badly degrades precision and introduces big errors i.e. instead of 0.013756 it returns 0.013744 |
|
|
vafonkin
Joined: 20 Mar 2012 Posts: 16
|
tan() function results on different MCUs |
Posted: Tue Mar 27, 2012 1:06 pm |
|
|
I want to ask for a little help from community - I do not have a lot of processors available but want to get a test results for the following function for different PICs:
tan()
atan()
atan2()
For 89 decimal degree. You need to convert it to Rad to get result.
Can someone help me?
The purpose is to compare calculation precision.
Thanks a lot in advance! |
|
|
Ttelmah
Joined: 11 Mar 2010 Posts: 19540
|
|
Posted: Tue Mar 27, 2012 2:40 pm |
|
|
The only PIC's offering a higher maths precision, are the DsPIC's, which do offer a double type. The library for all the others is basically the same.
Best Wishes |
|
|
|