arctan approx..

jondaeh

New Member
Hi!

I need to calculate the (approximate) angle between a gps unit and a waypoint(correct, it's a autopilot project). I have got all the GPS data safely in to the Picaxe 28X1, and only waiting for a calculation algorithm.

I guess it's a arctan algorithm I'm looking for. Using it like:

angle = arctan(latitude/longitude),

and a search returned some interesting code. But I find the most of it very complex and accurate, and I'm satisfied with a simple, inaccurate one which I understand, rather than just copy-paste one. I found this thread, where CPEDW described an seemingly easy formula:

http://www.picaxeforum.co.uk/showthread.php?t=2771

However, I can't understand how to use it. Consider the line:

z=abs(x/y), 0<z<1

what's abs? And what sign is < ?

If anyone got another formula/algorithm, I will be grateful to get a peek at it:)!

Thanks!

Jon
Norway
 

BeanieBots

Moderator
The 28X1 supports SIN & COS.
I'm sure you can derive ARCTAN using those. (COS/SIN??)

PS. ABS is the "absolute value of".
That is, any minus sign is ignored.
 

Jeremy Leach

Senior Member
How approximate is acceptable though? I mean, if you need the nearest 45 degrees then all you need to do is compare the LAT and LONG values.
 

Andrew Cowan

Senior Member
The 28X1 supports SIN & COS.
Amazing!

Sin/Cos = tan
Cos/Sin = cot

Not sure how you could get arctan...

Code:
SIN and COS
The sin function returns a number equivalent to the sine of the value in degrees.
The system uses a 45 step lookup table in each quadrant, giving a very fast, and
reasonably accurate, result.
The sine function only works on positive whole integers. However as all sin and
cos values repeat every 360 degrees, simply add 360 to make a negative value
positive. e.g. sin (-30) is the same as sin (330) (-30 + 360)
As the real sine value is always a value between 1 and -1, a coding system is used
to increase the accuracy when working with PICAXE whole integers. The value
returned by the sin function is actually 100 x the real sine value. Therefore in
normal mathematics sin 30 = 0.5. In PICAXE mathematics this is returned as 50
(100*0.5). This coding method provides a sine function accuracy equivalent to
two decimal places.
e.g let b1 = sin 30 (answer b1 = 50)
Negative numbers are indicated by setting bit 7 of the returned byte. This has the
effect of making negative values appear as 128 + the expected value.
e.g let b1 = sin 210 (answer b1 = 128+50 = 178)
The cos function operates in an identical manner.
A
 
Last edited:

hippy

Ex-Staff (retired)
One of the uM-FPU floating point co-processors ( FPU001, FPU003 ) may be the best option.

Worth considering is how the angle is to be used once you have it. If you are going to apply further trigonometric processing to it you may be able to go straight from latitude / longitude to whatever you need with no or simpler trigonometry.
 

skyhawk

New Member
I would recommend a table lookup. You only need to be able to calculate for values of the argument less than or equal to one. For arguments greater than one take the reciprocal and compute the complementary angle.

For example, store in a table the values of arctan for .1, .2, .3, .4, .5, .6, .7, .8, .9, and 1. Then find the pair that brackets the given argument and interpolate between the two. If lesser accuracy is acceptable, then just take the closest table entry. If you need more accuracy then add more table entries.
 

Jeremy Leach

Senior Member
I would perhaps use a lookup too (and this old post of mine did precisely that, but probably overkill here ! http://www.picaxeforum.co.uk/showthread.php?t=10362&highlight=arctan).

However what about using iterative approximations. i.e in pseudo code ...

Code:
'This loop will eventually give the result of Angle = ArcTan(TargetTan) 

Angle = InitialValue 'Set the initial value very roughly using a bit of logic
LOOP
    TanAngle = Sin(Angle)/Cos(Angle)
    TanError = TanAngle - TargetTan 
    Angle = Angle + .... (increment or decrement Angle based on TanError) 
UNTIL Abs(TanError) <=Threshold
This just gives the general idea and you'd have to flesh it out with sensible logic and write it in Picaxe Basic! But it would make use of the inbuilt Sin and Cos functions, which seem a bit of a shame to waste. You might also need to keep a check of the iterations and bomb out if too many.

However it 'might' be doable and pretty accurate with a modest amount of code. Not sure?
 
Last edited:

womai

Senior Member
Probably the simplest is a lookup table. You only need to cover the range from 0 to 90 degrees. The rest can be reduced to this case.

Or you can implement a search routine that given x calculates

tan(y) = sin(y)/cos(y)

where y is some approximation (or guess) to the final result, compare it to the required value, and based on the comparison either increment or decrement y, until x = sin(y)/cos(y) with sufficient accuracy. y is then actan(x).

Another method to calculate actan(x) is using a Taylor series (polynomial) approximation:

arctan(x) = x - x^3/3 + x^5/5 - x^7/7 + x^9/9 - x^11/11 + ....

Note that here x is in Radian and -1 < x < 1. To go from degrees to radian:

x(Radian) = x(Degrees) * pi/180 (pi = 3.1415927)

All these floating point operations will be a challenge to implement on a Picaxe, so I'd really suggest going with a lookup table.

Wolfgang

P.S.: While writing this post others have made the same suggestion...
 

Jeremy Leach

Senior Member
Couldn't resist trying some code - see snippets section for iterative routine. Seems to give pretty good results and doesn't use too much code ;) Not saying it's better than lookup though.
 

skyhawk

New Member
The OP says that the method does not need to be particularly accurate, but does not state how accurate the method needs to be. The method I previously posted will provide accuracy to one degree using a minimal amount of integer arithmetic.

I repeat that one only needs to consider the range 0 to 45 degrees corresponding to an argument of 0 to 1. Angles between 45 and 90 degrees are handled by working with the complementary angle. Remember that tan(90) is infinite.

Working with a Taylor series on a PICAXE is needlessly computationally expensive. A table lookup with linear interpolation is computationally easy. Scale the argument by 100 so that it ranges from 0 to 100. Then do a search to find which two entries in the table the argument falls between. For example, an argument of 73 (.73) would fall between 70 (.7) and .80 (.8). The corresponding table entries would be 35 and 39. So one would take 35 plus 3x(39-35)/10, which gives 36 using integer arithmetic. The exact answer is 36.13.
 
Last edited:

crossthreaded

New Member
I needed arcsin for calculating inclination using a force transducer recently and I used interpolation with 6 data points covering 45 degrees. The data points weren't uniformly spread, but chosen to give better accuracy by concentrating them at the bendy bit of the curve. I only needed enough accuracy to drive a RC servo, so one degree was plenty close enough.

The result covered 360 degrees by keeping track of the sign and using trig identities. This meant doing square-root too, for which a fast iterative technique was quite adequate.

I mocked it up in Excel first, but I haven't got the code with me right now.
 

Jeremy Leach

Senior Member
Ok, I agree Skyhawk. And I've proved it by writing out some code based on your suggestion, which is very simple with good accuracy ;) I've ammended my code in snippets section (with credit to yourself!) but also listed here:

Code:
Symbol Offset	= b0
Symbol Result	= b0
Symbol Index 	= b1
Symbol ArcTan0	= b2
Symbol ArcTan1	= b3
Symbol Tan100	= b4

EEPROM 0, (0,5,11,16,21,26,30,34,38,41,45)

Test:
	Tan100 = 82
	Gosub ArcTan
	Stop


ArcTan:
	'ON ENTRY: Tan100 holds the argument * 100
	'ON EXIT: Result holds the answer in degrees.
	'NOTE: This only works for angles 0 to 45.
	
	Index = Tan100 / 10
	Offset = Index * 10
	Offset = Tan100 - Offset
	
	Read Index,ArcTan0
	Inc Index
	Read Index,ArcTan1
	
	Result = ArcTan1 - ArcTan0 * Offset / 10 + ArcTan0
Return
Incidentally ArcSin Routine could be based on the same routine, but with different EEPROM values.
 
Last edited:

AliasMrJones

New Member
Some Code

Hi, I'm the builder of Deathpod3000, which was referenced in a previous post. I'm in the process of writing some new articles on my build blog with code samples and thought I'd share them here as well. OP didn't ask for distance, but for navigation, this is pretty useful. I use the haversine formula as other formulas don't work well for 8 bit micros and small distances. You also could use simpler trig functions and ignore the curvature of the earth for small distances. I used an Atmel AVR so this is in C, but should translate easily.
Code:
float calcDistance(float lat1, float lon1, float lat2, float lon2)
{
	float dlon, dlat, a, c;
	float dist = 0.0;
	lon1 = lon1 * -1.0; //make west = positive
	lon2 = lon2 * -1.0;
        dlon = lon2 - lon1;
        dlat = lat2 - lat1;
        a = pow(sin(dlat/2),2) + cos(lat1) * cos(lat2) * pow(sin(dlon/2),2);
        c = 2 * atan2(sqrt(a), sqrt(1-a));
        dist = 20925656.2 * c;  //radius of the earth (6378140 meters) in feet 20925656.2
	
	return(dist);
}
Lat and lon in the above function must be in radians. I wrote 2 small functions to convert between degrees for display and radians for calculations.

Code:
float dtor(float degrees)
{
	return(degrees * PI / 180);
}
float rtod(float radians)
{
	return(radians * 180.0 / PI);
}
You can read more about the haversine formula and other method of calculating distance and bearing here:

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

The code I use to calculate bearing (again lat and lon must be in radians) is:

Code:
float calcBearing(float lat1, float lon1, float lat2, float lon2)
{
	float bearing = 0.0;
	//determine angle
	bearing = atan2(sin(lon2-lon1)*cos(lat2), (cos(lat1)*sin(lat2))-(sin(lat1)*cos(lat2)*cos(lon2-lon1)));
	//convert to degrees
	bearing = rtod(bearing);
	//use mod to turn -90 = 270
	bearing = fmod((bearing + 360.0), 360);
	return bearing;
}
Be careful with the order of the parameters for atan2 as the avr gcc libs and Excel take the parameters in opposite order and this messed me up at first. This function takes parameters in radians and returns bearing in degrees. To convert negative values to correct positive 0-360 degree values, 360 is added to the bearing and then mod 360 (remainder when divided by 360).

I also created a relative bearing function so after computing bearing and reading my digital compass I could figure out which direction to turn and how many degrees.

Code:
float CalcRelBearing(float bearing, float heading)
{
	float relBearing;
	relBearing = bearing - heading;
	if(relBearing > 180.0)
	{
		relBearing = bearing - (heading + 360);
	}
	if(relBearing < -180.0)
	{
		relBearing = bearing - (heading - 360);
	}
	return(relBearing);
}
Bearing is where you want to go and heading is where you are going now. These are in degrees. The function returns a value of -180 - +180 for left 0-180 degrees or right 0-180 degrees. I then used a proportional function to send appropriate servo pulses to the steering servo to turn the car.

Lookup tables would work for rough values of some of the trig functions if they aren't available. I used that method on a previous rover where the floating point stuff took too long. I thought about using fixed point and lookup tables, but decided to try the float stuff and see how it performed. On an Atmega32 at 8mhz, using the free gcc C compiler, I was getting 5-10 times through my main navigation loop per second, which seemed plenty fast so I just stuck with the floating point math as the trig fuctions were build into the compiler library so it was easy to implement.

Deathpod3000 was one of only 2 vehicles, and the only ground vehicle, to complete the course at the recent Sparkfun Autonomous Vehicle Competition. The other vehicle was an airplane. Airplane > car in a speed race so we got second place, but also got the Engineer's Choice award. I have many details about the construction and programming of the car and pics and videos from the competition on my build blog here:

http://deathpod3000.wordpress.com

Check it out if you're interested.

Eric
 

jondaeh

New Member
Thanks for great answers!

There was some questions regarding how exact I needed the approximation. If it's inside 2-3 degrees I'm satistfied.

I get the impression that most recommend a look up table, and it seem like a good solution. But I wonder, if arctan actually is as simple as cos/sin (why didn't we learn that at school?), isn't that the simplest solution? Or will it take way too much memory from the Picaxe?


Jon.
 

womai

Senior Member
No, arctan(x) is NOT cos(x)/sin(x). That would be cot(x) (cotangens).

arctan(x) is the inverse operation to tan(x), same relationship as e.g. square root and square.

actan(x) = y, where tan(y) = x

And yes, a lookup table (in combination with linear interpolation between the data points if necessary) will be the fastest solution, and as mentioned can be applied to other functions just as well.
 
Top