Listing 1 CORDIC.C

/*  CORDIC.C : Demonstrates the integer-based CORDIC
 *  system for calculating sines and cosines. The
 *  vertices of a regular hexagon are calculated using
 *  CORDIC trig and the hexagon itself is rotated.
 *  Make in Borland C++'s internal environment.
 *  (c) 1991 Michael Bertrand.
 */

#include <stdio.h>     /* printf */
#include <math.h>      /* sqrt, atan */
#include <conio.h>     /* getch */
#include <graphics.h>  /* BGI */

typedef struct
{
  int x;
  int y;
} POINT;

typedef unsigned int WORD;

void CalcHexPtsCORDIC(POINT center, POINT vertex);
void DrawHexagon(POINT center, POINT vertex);
void DrawCross(POINT pt, int colr);
void SinCosSetup(void);
void SinCos(WORD theta, int *sin, int *cos);

#define ESC 0x1B

/* Quadrants for angles. */
#define QUAD1 1
#define QUAD2 2
#define QUAD3 3
#define QUAD4 4

/* NBITS is number of bits for CORDIC units. */
#define NBITS 14

/* NUM_PTS is number of vertices in polygon. */
#define NUM_PTS 6

int  ArcTan[NBITS];       /* angles for rotation */
int  xInit;               /* initial x projection */
WORD CordicBase;          /* base for CORDIC units */
WORD HalfBase;            /* CordicBase / 2 */
WORD Quad2Boundary;       /* 2 * CordicBase */
WORD Quad3Boundary;       /* 3 * CordicBase */
POINT HexPts[NUM_PTS+1];  /* calculated poly points */

void main(void)
{
  int   driver;   /* for initgraph() */
  int   mode;     /* for initgraph() */
  WORD  theta;    /* CORDIC angle */
  int   sine;     /* sine of CORDIC angle */
  int   cosine;   /* cosine of CORDIC angle */
  POINT center;   /* center of hexagon */
  POINT vertex;   /* hexagon's original base vertex */
  POINT vertex1;  /* hexagon's changing base vertex */
  POINT del;      /* vertex - center (radial spoke) */
  int   radius;   /* radius of circumscribing circle */

  driver = VGA;     /* for initgraph() */
  mode   = VGAHI;   /* mode 0x12 : 640x480 16 color */

  if (registerbgidriver(EGAVGA_driver) < 0)
    {
    printf("couldn't find VGA driver"); return;
    }
  initgraph(&driver, &mode, NULL);

  printf("Press ENTER to rotate, ESC to quit.");

  center.x = 320; center.y = 240;
  vertex.x = 470; vertex.y = 240;
  radius = vertex.x - center.x;
  /* Calculate the radial spoke : vertex - center. */
  del.x = vertex.x - center.x;
  del.y = vertex.y - center.y;

  setwritemode(XOR_PUT);

  /* Draw circumscribing circle. */
  setcolor(RED);
  circle(center.x, center.y, radius);

  /* Draw small cross at center. */
  DrawCross(center, YELLOW);
  setcolor(WHITE);

  /* Setup CORDIC system, initialize theta = 0. */
  SinCosSetup();
  theta = 0;

  /* Draw initial hexagon. */
  DrawHexagon(center, vertex1 = vertex);

  /* Rotate hexagon. vertex is fixed; vertex1 rotates
   * clockwise around vertex in increments of 650
   * CORDIC units (3.57 deg). CORDIC sines/cosines are
   * used to find vertex1. DrawHexagon() also uses
   * CORDIC sines/cosines to calculate the remaining
   * vertices for each hexagon so they can be drawn.
   */

  while (getch() ! = ESC)
    {
    /* Erase last hexagon. */
    DrawHexagon(center, vertex1);
    /* Inc theta by 650 CORDIC units (3.57 deg). */
    theta += 650;
    SinCos(theta, &sine, &cosine);
    /* Calc new vertex by rotating around center. */
    vertex1.x = (int)
      (((long) del.x * cosine - (long) del.y * sine +
      HalfBase) >> NBIIS) + center.x;
    vertexl.y = (int)
      (((long) del.x * sine + (long) del.y * cosine +
      HalfBase) >> NBITS) + center.y;
    /* Draw new hexagon. */
    DrawHexagon(center, vertex1);
    }

  closegraph();
}

void CalcHexPtsCORDIC(POINT center, POINT vertex)
/*
  USE:  Calc array of hex vertices using CORDIC calcs.
  IN:   center = center of hexagon.
       vertex = One of the hexagon vertices
  NOTE: Loads global array HexPoints[] with other 5
   vertices of regular hexagon. Uses CORDIC routines
   for trig and long integer calculations.
*/
{
  int   sine;     /* sine of central angle */
  int   cosine;   /* cosine of central angle */
  int   corner;   /* index for vertices of polygon */
  POINT del;      /* vertex - center (radial spoke) */

  /* 60 deg = 10923 CORDIC units. */
  SinCos(10923, &sine, &cosine);

  /* Set initial and final point to incoming vertex. */
  HexPts[0].x = HexPts[NUM PTS].x = vertex.x;
  HexPts[0].y = HexPts[NUM_PTS].y = vertex.y;

  /* Go clockwise around circle to calc hex points. */
  for (corner = 1; corner < NUM_PTS; corner++)
    {
    /* Calculate the radial spoke : vertex - center. */
    del.x = vertex.x - center.x;
    del.y = vertex.y - center.y;
    /* calc new vertex by rotating around center. */
    vertex.x = (int)
      (((long) del.x * cosine - (long) del.y * sine +
      HalfBase) >> NBITS) + center.x;
    vertex.y = (int)
      (((long) del.x * sine + (long) del.y * cosine +
      HalfBase) >> NBITS) + center.y;
    /* Store new vertex in array. */
    HexPts[corner].x = vertex.x;
    HexPts[corner].y = vertex.y;
    }
}

void DrawHexagon(POINT center, POINT vertex)
/*
  USE:  Draw Hexagon given center and one vertex.
  IN:   center = center of hexagon.
       vertex = one of the hexagon vertices
  NOTE: Call CalcHexPtsCORDIC() to load global array
   HexPts[] with hexagon vertices.
*/
{
  CalcHexPtsCORDIC(center, vertex);
  drawpoly(NUM_PTS+1, (int far *)HexPts);
}

void DrawCross(POINT pt, int colr)
/*
  USE: Draw cross on screen at pt with given color.
*/
{
  int oldColor;

  setwritemode(COPY_PUT);
  oldColor = getcolor();
  setcolor(colr);
  moveto(pt.x - 2, pt.y); lineto(pt.x + 2, pt.y);
  moveto(pt.x, pt.y - 2); lineto(pt.x, pt.y + 2);
  setcolor(oldColor);
  setwritemode(XOR_PUT);
}

void SinCosSetup(void)
/*
  USE : Load globals used by SinCos().
  OUT : Loads globals used in SinCos() :
   CordicBase    = base for CORDIC units
   HalfBase      = Cordicbase / 2
   Quad2Boundary = 2 * CordicBase
   Quad3Boundary = 3 * CordicBase
   ArcTan[]      = the arctans of 1/(2^i)
   xInit         = initial value for x projection
  NOTE: Must be called once only to initialize before
   calling SinCos(). xInit is sufficiently less than
   CordicBase to exactly compensate for the expansion
   in each rotation.
*/
{
  int i;        /* to index ArcTan[] */
  double f;     /* to calc initial x projection */
  long powr;    /* powers of 2 up to 2^(2*(NBITS-1)) */

  CordicBase = 1 << NBITS;
  HalfBase = CordicBase >> 1;
  Quad2Boundary = CordicBase << 1;
  Quad3Boundary = CordicBase + Quad2Boundary;

  /* ArcTan's are diminishingly small angles. */
  powr = 1;
  for (i = 0; i < NBITS; i++)
    {
    ArcTan[i] = (int)
      (atan(1.0/powr)/(M_PI/2)*CordicBase + 0.5);
    powr << = 1;
    }
  /* xInit is initial value of x projection to
   * compensate for expansions. f = 1/sqrt(2/1 * 5/4 * ...
   * Normalize as an NBITS binary fraction (multiply by
   * CordicBase) and store in xInit. Get f = 0.607253
   * and xInit = 9949 = 0x26DD for NBITS = 14.
   */
  f = 1.0;
  powr = 1;
  for (i = 0; i < NBITS; i++)
   {
   f = (f * (powr + 1)) / powr;
   powr <<= 2;
   }
  f = 1.0/sqrt(f);
  xInit = (int) (CordicBase * f + 0.5);
}

void SinCos(WORD theta, int *sin, int *cos)
/*
  USE : Calc sin and cos with integer CORDIC routine.
  IN  : theta = incoming angle (in CORDIC angle units)
  OUT : sin = ptr to sin (in CORDIC fixed point units)
       cos = ptr to cos (in CORDIC fixed point units)
  NOTE: The incoming angle theta is in CORDIC angle
   units, which subdivide the circle into 64K parts,
   with 0 deg = 0, 90 deg = 16384 (CordicBase), 180 deg
   = 32768, 270 deg = 49152, etc. The calculated sine
   and cosine are in CORDIC fixed point units : an int
   considered as a fraction of 16384 (CordicBase).
*/
{
  int quadrant;  /* quadrant of incoming angle */
  int z;         /* incoming angle moved to 1st quad */
  int i;         /* to index rotations : one per bit */
  int x, y;      /* projections onto axes */
  int xl, yl;    /* projections of rotated vector */

  /* Determine quadrant of incoming angle, translate to
   * 1st quadrant. Note use of previously calculated
   * values CordicBase, etc. for speed.
   */

  if (theta < CordicBase)
    {
    quadrant = QUAD1;
    z = (int) theta;
    }
  else if (theta < Quad2Boundary)
    {
    quadrant = QUAD2;
    z = (int) (Quad2Boundary - theta);
    }
  else if (theta < Quad3Boundary)
    {
    quadrant = QUAD3;
    z = (int) (theta - Quad2Boundary);
    }
  else
    {
    quadrant = QUAD4;
    z = - ((int) theta);
    }

  /* Initialize projections. */
  x = xInit;
  y = 0;

  /* Negate z, so same rotations taking angle z to 0
   * will take (x, y) = (xInit, 0) to (*cos, *sin).
   */
  z = -z;

  /* Rotate NBITS times. */
  for (i = 0; i < NBITS; i++)
    {
    if (z < 0)   
      {
      /* Counter-clockwise rotation. */
      z += ArcTan[i];
      y1 = y + (x >> i);
      x1 = x - (y >> i);
      }
    else
      {
      /* Clockwise rotation. */
      z -= ArcTan[i];
      y1 = y - (x >> i);
      x1 = x + (y >> i);
      }

    /* Put new projections into (x,y) for next go. */
    x = x1;
    y = y1;
    }  /* for i */

  /* Attach signs depending on quadrant. */
  *cos = (quadrant==QUAD1 || quadrant==QUAD4) ? x : -x;
  *sin = (quadrant==QUAD1 || quadrant==QUAD2) ? y : -y;
}
/* End of File */