This routine will take an 8 bit integer that corresponds to the numerator of a fraction whose denominator is 256 and find its arctangent. So the input ranges from 0 to 255 which corresponds to 0 to 255/256 = 0.996 . The output for an arctangent routine that returns a floating point number would be from 0 (atan(0)) to 0.783 (atan(255/256)) radians; or if you prefer, 0 to 44.89 degrees. However, this routine scales the output so that pi/4 radians (or 45 degrees) corresponds to 256. So for the input range of 0 to 255 you get an output of 0 to 255 ( atan(255/256) * 256 / (pi/4) is about 255). It's probably a little more interesting to see an intermediate data point or two:
Intger Float x | atan(x) | x | atan(x) ------+----------+-------+--------- 0x4a | 0x5b | .289 | .281 0x62 | 0x77 | .383 | .366 0x6f | 0x84 | .434 | .409 0xa6 | 0xbb | .648 | .575 0xdb | 0xe6 | .855 | .707
The only thing that's left is combining the fractional division and the swapping of the x and y values if y is greater than x (and then subtracting the result from pi/2 or actually 512 in this case).
;----------------------------------------------------------
;
;arctan (as adapted from the similar arcsin function)
;
; The purpose of this routine is to take the arctan of an
;8-bit number that ranges from 0 < x < 255/256. In other
;words, the input, x, is an unsigned fraction whose implicit
;divisor is 256.
; The output is in a conveniently awkward format of binary
;radians (brads?). The output corresponds to the range of zero
;to pi/4 for the normal arctan function. Specifically, this
;algorithm computes:
;
; arctan(x) = real_arctan(x/256) * 256 / (pi/4)
; for 0 <= x <= 255
;
; where, real_arctan returns the real arctan of its argument
;in radians.
;
; The algorithm is a table look-up algorithm plus first order
;linear interpolation. The psuedo code is:
;
;unsigned char arctan(unsigned char x)
;{
; unsigned char i;
;
; i = x >> 4;
; return(arctan[i] + ((arctan[i+1] - arctan[i]) * (x & 0xf))/16);
;}
;
;
arctan
SWAPF x,W
ANDLW 0xf
ADDLW 1
MOVWF temp ;Temporarily store the index
CALL arc_tan_table ;Get a2=atan( (x>>4) + 1)
MOVWF result ;Store temporarily in result
DECF temp,W ;Get the saved index
CALL arc_tan_table ;Get a1=atan( (x>>4) )
SUBWF result,W ;W=a2-a1, This is always positive.
SUBWF result,F ;a1 = a1 - (a1-W) = W
CLRF temp ;Clear the product
CLRC
BTFSC x,0
ADDWF temp,F
RRF temp,F
CLRC
BTFSC x,1
ADDWF temp,F
RRF temp,F
CLRC
BTFSC x,2
ADDWF temp,F
RRF temp,F
CLRC
BTFSC x,3
ADDWF temp,F
RRF temp,W
ADDWF result,F
RETURN
arc_tan_table
ADDWF PCL,F
RETLW 0
RETLW 20 ;atan(1/16) = 3.576deg * 256/45
RETLW 41
RETLW 60
RETLW 80
RETLW 99
RETLW 117
RETLW 134
RETLW 151
RETLW 167
RETLW 182
RETLW 196
RETLW 210
RETLW 222
RETLW 234
RETLW 245
RETLW 0 ;atan(32/32) = 45deg * 256/45
The other part of the problem is implementing the division
FRAC_DIV:
;-------------------
;Fractional division
;
; Given x,y this routine finds:
; a = 256 * y / x
;
movlw 8 ;number of bits in the result
movwf cnt
clrf a ; the result
movf x,w
L1:
clrc
rlf y,f ;if msb of y is set we know x<y
rlf a,f ;and that the lsb of 'a' should be set
subwf y,f ;But we still need to subtract the
;divisor from the dividend just in
;case y is less than 256.
skpnc ;If y>x, but y<256
bsf a,0 ; we still need to set a:0
btfss a,0 ;If y<x then we shouldn't have
addwf y,f ;done the subtraction
decfsz cnt,f
goto L1
return
It's easy enough to combine these two routines to obtain a 4-quadrant arctan(y/x) routine. However, you do need to keed in mind that the arctangent routine posted above is only valid over 1/8th of the unit circle. To obtain the other 7/8th's you'll need to apply the appropriate trig identities.
// pic routines:
extern int arctan(int x);
extern int frac_div(int y, int x);
// Untested c code that implements a 4-quadrant arctan function
int arctan(int x, int y)
{
int f, swapped;
int reference_angle;
swapped = 0;
if(x < 0) {
if(y < 0)
reference_angle = 256 * 2; // pi
else
reference_angle = 256; // pi/2
} else {
if( y < 0)
reference_angle = 256 * 3; // 3*pi/2
else
reference_angle = 0;
}
if (x<=y) {
f = y;
y = x;
x = f;
swapped = 1;
}
f = frac_div(y,x);
f = arctan(f);
if (swapped)
f = 256 - f;
return f + reference_angle;
}