On Sat, 3 Jul 1999, Dave VanHorn wrote:

> I'm on a quest, to find an elegant way to convert from rectangular
> coordinates to polar coordinates.
>
> What I have, is a pair of 8 bit signed values, for the X and Y coordinates.
> I have the radius, which wasn't as ugly as I thought it would be.
> Finding the angle in any sort of elegant manner is proving elusive.
> Worse, the only manner i know to calculate it is A=csc(X/radius), which
> ignores the Y value completely, and has poor resolution where X is near it's
> maxima and minima.
>

As others have commented already, it sounds as though if you want an
arctangent as opposed to a 1/sine (though I think you may've meant
arcsine). As I see it, or I guess saw it, there are two problems you need
to solve: 1) Computing the arctangent and 2) performing the division. I
say saw since that's the way I ended up doing it for Jason Wolfson back in
April. Here's what I wrote (privately) to Jason:


;--------------------------------------------------------------

Jason,

I really don't want to 'spoil your fun'. So if you're not interested in a
solution, then don't look at the code that's below :).

I had a little time this afternoon so I decided to hack my arcsine routine
to create the arctan routine. What you see below does not solve your
problem completely - but it is the major part. 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

Anyway, I thought I'd share this in case you're interested. 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).

Regards,
Scott




;----------------------------------------------------------
;
;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.


Off the top of my head and in C:


// 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;
}


Whew!

Scott