Make your own free website on Tripod.com

 


 

/******************************************************************************

sine.c - an integer sine function using a 16 bit integer table with

         interpolation.

(C language program, Rev. June 23, 1986, from 1962 FORTRAN program by JS)

cyberanimation@lycos.com, cyberanimation.tripod.com, © 2003 Cybernetic Cinema

/*

         subprograms include:  SINE, ANGINC, ROTATE (all use the TABL.C module)

/*

Sine (a),  a = given angle +180 to -180 degrees, as: +32767 to -32767, the

result is the sine value, range: +32767 to -32767.  (Range selected to

maximize 16 bit integer angle & sine resolution and speed.)  Table size = 512

elements per 90 degrees (2048 per circle) or .17578125 deg per element.  The

sine increment is Pi 3.14159265 / 1024 = .00306796; the most significant 9

bits (6 thru 14) selects the table element, bit 15 determines the sign,

bit 14 on indicates +-90 thru 180 deg & a reverse of table direction.  The

least significant 6 bits (0 thru 5) used for interpolation between table

sine values.  The cosine is the sine of the (angle + 90 degrees or 16384).

/*

Example:    sinval = sine (ang);

/*

Sin (a),  Same as "sine" function described but without interpolation, used

        when extra angular resolution (greater than .1757 deg) not needed.

/*

Anginc (a,i,s,c);   a = angle +-32767,  i = angle increment (input),

                s = sine result,    c = cosine result (+-32767).

                Example:         anginc (ang, incr, sax, cax);

/*

Rotate ();  Given a point in global memory (x, y, z) and sine/cosine

              values for 3 axis, then rotate point with result in xr,yr,zr

            and take perspective given focal length into xr,yr.

            Example:      rotate ();

/*

sintab - an integer table of sine values, 512 elements extending from 0 to 90

       degrees.  0 to +1. is scaled to the integer range of 0 to 32767.

                                    Rev. July 7, 1986, JS

/*

         The angular increment = 3.14159 / 1024 = .00306796 radians.

         The 0 to 1 result is multiplied by 32767 then rounded to the

         nearest integer and written to a disk file by a Basic (BASICA)

         program approximately as follows:

                A=0   D=0

            OPEN "SINDAT.TXT" FOR OUTPUT AS #1

       400  B = SIN(A)

            C = B * 32767.

            A = A + .00306796

            LPRINT USING " #####.######";A,B,C,D

            PRINT#1,USING "#####.,";C

            D = D + 1

            IF (A > 1.6) STOP

            GO TO 400

            END

/*

******************************************************************************

*/

 unsigned short sintab[] =  {0,101,201,302,402,503,603,704,804,905,1005,

 1106,1206,1307,1407,1507,1608,1708,1809,1909,2009,2110,2210,2310,2410,2511,

 2611,2711,2811,2911,3012,3112,3212,3312,3412,3512,3612,3712,3811,3911,4011,

 4111,4210,4310,4410,4509,4609,4708,4808,4907,5007,5106,5205,5305,5404,5503,

 5602,5701,5800,5899,5998,6096,6195,6294,6393,6491,6590,6688,6786,6885,6983,

 7081,7179,7277,7375,7473,7571,7669,7767,7864,7962,8059,8157,8254,8351,8448,

 8545,8642,8739,8836,8933,9030,9126,9223,9319,9416,9512,9608,9704,9800,9896,

 9992,10087,10183,10278,10374,10469,10564,10659,10754,10849,10944,11039,

11133,11228,11322,11417,11511,11605,11699,11793,11886,11980,12074,12167,

12260,12353,12446,12539,12632,12725,12817,12910,13002,13094,13187,13279,

13370,13462,13554,13645,13736,13828,13919,14010,14101,14191,14282,14372,

14462,14553,14643,14732,14822,14912,15001,15090,15180,15269,15358,15446,

15535,15623,15712,15800,15888,15976,16063,16151,16238,16325,16413,16499,

16586,16673,16759,16846,16932,17018,17104,17189,17275,17360,17445,17530,

17615,17700,17784,17869,17953,18037,18121,18204,18288,18371,18454,18537,

18620,18703,18785,18868,18950,19032,19113,19195,19276,19357,19438,19519,

19600,19680,19761,19841,19921,20000,20080,20159,20238,20317,20396,20475,

20553,20631,20709,20787,20865,20942,21019,21096,21173,21250,21326,21403,

21479,21554,21630,21705,21781,21856,21930,22005,22079,22154,22228,22301,

22375,22448,22521,22594,22667,22739,22812,22884,22956,23027,23099,23170,

23241,23312,23382,23452,23522,23592,23662,23731,23801,23870,23938,24007,

24075,24143,24211,24279,24346,24413,24480,24547,24613,24680,24746,24811,

24877,24942,25007,25072,25137,25201,25265,25329,25393,25456,25520,25582,

25645,25708,25770,25832,25894,25955,26016,26077,26138,26198,26259,26319,

26378,26438,26497,26556,26615,26674,26732,26790,26848,26905,26962,27019,

27076,27133,27189,27245,27301,27356,27411,27466,27521,27575,27629,27683,

27737,27790,27844,27896,27949,28001,28053,28105,28157,28208,28259,28310,

28360,28411,28460,28510,28560,28609,28658,28706,28755,28803,28850,28898,

28945,28992,29039,29085,29131,29177,29223,29268,29313,29358,29403,29447,

29491,29535,29578,29621,29664,29706,29749,29791,29832,29874,29915,29956,

29997,30037,30077,30117,30156,30195,30234,30273,30311,30349,30387,30424,

30462,30498,30535,30571,30607,30643,30679,30714,30749,30783,30818,30852,

30885,30919,30952,30985,31017,31050,31082,31113,31145,31176,31207,31237,

31267,31297,31327,31356,31385,31414,31442,31470,31498,31526,31553,31580,

31607,31633,31659,31685,31710,31736,31760,31785,31809,31833,31857,31880,

31903,31926,31949,31971,31993,32014,32036,32057,32077,32098,32118,32137,

32157,32176,32195,32213,32232,32250,32267,32285,32302,32318,32335,32351,

32367,32382,32397,32412,32427,32441,32455,32469,32482,32495,32508,32521,

32533,32545,32556,32567,32578,32589,32599,32609,32619,32628,32637,32646,

32655,32663,32671,32678,32685,32692,32699,32705,32711,32717,32722,32728,

32732,32737,32741,32745,32748,32752,32755,32757,32759,32761,32763,32765,

32766,32766,32767,32767,32767,32766,32766,32765,32763,32761,32759,32757,

32755,32752,32748,32745,32741,32737,32732,32728,32722,32717,32711 };

/*

/*    int  xcomp, ycomp, nbrend;/* Output from COMPNT, new compass point */

/*    float  fdis;              /* Input to COMPNT, orig compass distance*/

/*    float  dist (int, int, int, int); /* Define a floating point funct */

/*    extern int kul, sintab[]; /* Sine values 0 to 90 degre, 0 to 32767 */

/*    extern int xlolimit, ylolimit, xhilimit, yhilimit, xc, yc, frmchg;

/*    int xt, yt, zt;                     /* Test variables **************/

#include <stdio.h>

#include <math.h>

/* All variables have to be initialized only for use on the Recog Module: */

      short  xc1=0, yc1=0, xp1=0, yp1=0;

      long   int  sax=0, cax=0, say=0, cay=0, saz=0, caz=0, xr=0, yr=0,

               zr=0, focal=0, sv=0, sr=0, xoff=0, yoff=0, zoff=0;

      int    xs=0, ys=0, zs=0, xper=0, yper=0;

/*

/*    main () {        /* TEST Integer SINE Function, print a table.

                 Sine error for angle +32768 or 8000 hex or -32768

                 Test range: 32752 to 32768, 16 unit interpolation range. */

/*    int j, k;

      k = 32750;

      for (j=0; j<=120; ++j) {

      sine (k);

      ++k; }

      return;  }

*/

    sine (angl) short angl; { /* Sine pgm., table look-up & interpolation*/

            short  aindx, sign, s1, s2, dif, sien;

                      /*      printi (" SINE angl=%",angl);

                      /* If angle is - set sign=-1 and angle to positive */

      if (angl < 0) { sign = -1;  angl = -angl; }  else  { sign = 0; }

           /* If angle 0 to 90 deg (0 to 16383) use table foreward as is */

      aindx = angl >> 5;     /* Right shift, get 10 signif bits for index*/

            /* printi(" aindx=%",aindx); printi(" sign=%",sign); */

      if ((angl&0x4000) == 0) goto S20; /* If 0 to +-90 deg use table fwd*/

           /* Angle is 90 to 180 deg (+ or -) index from the table bottom*/

      aindx = 1024 - aindx;   /* Reverse index direction, bottom to top  */

            /* printi (" 2nd aindx=%",aindx); */

      s1 = sintab[aindx];     /* Obtain sine value from "tabl" by look-up*/

      s2 = sintab[aindx-1];   /* Get next sine value from reversed table */

      goto S30;         /* Compute sine differences for interpltion.

            If the table index AINDX is negative, this is because the

            angle (which varies between +32k and -32k for 360 deg) added

            up to +32768 which is equivilant to -32768 as 16 bit integer,

            however, when an attempt is made to take the absolute value

            of -32768 an impossibility occures, so the absolute value

            comes out negative: -32768 or 8000 hex, then when shifted 5

            bits right the result is an array index of -1024 instead of a

            positive index of 1024.  Since this is a unique case that can

            occure and the SIEN result should be 0, set it and exit.     */

S20:  if (aindx < 0) {sien=0; goto S50;} /* If AINDX=-, SIEN=0, and exit */

      s1 = sintab[aindx];     /* Obtain sine value from "tabl" by look-up*/

      s2 = sintab[aindx+1];   /* Get next sine value from foreward table */

S30:  dif = s2 - s1;          /* Compute difference between sine values  */

      sien = s1 + ((dif * (angl & 0x1F)) >> 5);  /* Linear interpolation */

      if (sign < 0) sien = -sien; /* If angle is negative set sine to "-"*/

S50:  /*    printi (" aindx=%",aindx);

      /*    printi (" s1=%",s1);      printi (" s2=%",s2);

      /*    printi (" dif=%",dif);    printi (" sien=%%",sien); */

      /*printf("\nSine angl=%d sign=%d aindx=%d s1=%d s2=%d dif=%d sien=%d",

                  angl,sign,aindx,s1,s2,dif,sien); */

      return (sien); }        /* Return functional sine value to caller. */

/*

*****************************************************************************

anginc - given a 16 bit angle (-32k to 32k = -180 to +180 degrees), increment

         the angle, and return integer Sine and Cosine values (-32k to +32k).

*****************************************************************************

*/

      anginc (angle, anginc, sinang, cosang) {/* Incrmt angle get sin-cos*/

      angle = angle + anginc;         /* Increment the given angle       */

      sinang = sine (angle);          /* Obtain sine value of the angle  */

      cosang = sine (angle + 16384);  /* Cosine =sine of angle+90 degrees*/

      return (angle, sinang, cosang);}/* Return a new: angle, sin, & cos */

/*

/*    rotest () { /*---------------------------------------TEST Rotation */

/*          int j, angtst=0, angtst1=0, angtst2=0, xv, yv;

/*          float ang;

            xs = 95;  ys = 5;  zs = 8;  focal = 100;

            sax = cax = say = cay = saz = caz = 0;

            ang = 1.570796;

      for (j=0; j<=4096; ++j) {

      saz = sine(angtst);

      caz = sine(angtst+16384);

      say = sine(angtst1);

      cay = sine(angtst1+16384);

      sax = sine(angtst2);

      cax = sine(angtst2+16384);

/*          sv = saz * 95;

/*          sv = caz * 95;

/*          crtwr (xv+160, yv+100, kul);

            printi(" CIRCLE xv=%",xv); printi(" yv=%",yv);

      xv = (cos(ang) * 90.) + .5;

      yv = (sin(ang) * 90.) + .5;

      crtwr (xv+120, yv+100, kul+1);

      printi(" FLOAT xv=%",xv); printi(" yv=%%",yv);

      ang = ang - .0981747;                      /*   3.14159265 / 32   */

/*          rotate();

            crtwr (xt+160, yt+100, kul); 

/*          crtwr (xper+180, yper+100, kul+1);  */

/*    angtst  = angtst + 1024;

      angtst1 = angtst1 + 400;

      angtst2 = angtst2 + 300; }

      return; }   */

/*

***************************************************************************

rotate - Three axis rotation for an x, y, z point.

***************************************************************************

*/

      rotate () {

/**/

           int foc, i, isaz;  /* Printing test variables *****************/

           long int sv;

      xr = xs << 4;           /* Mult by 16 to reduce fractional errors  */

      yr = ys << 4;           /* Convert 16 bit integer to 32 bit integer*/

      zr = zs << 4;

      sv = xr;                /* Use XR in its unmodified to compute YR. */

      sr = (xr * caz) - (yr * saz);       /* "Z" axis rotation, given sin*/

      i = sr & 0x4000;

      if (i != 0) {xr=(sr>>15)+1;} else {xr=sr>>15;}          /* Round X */

      sr = (sv * saz) + (yr * caz);       /* xr & yr = point rotated on Z*/

      i = sr & 0x4000;

      if (i != 0) {yr=(sr>>15)+1;} else {yr=sr>>15;}          /* Round Y */

ROTY:                         /* isaz=xr;  printi("Rotated:Z xr=%",isaz);

                                         isaz=yr;  printi(" yr=%",isaz); */

      sv = xr;

      sr = (xr * cay) - (zr * say);       /* "Y" axis rotation, given sin*/

      i = sr & 0x4000;

      if (i != 0) {xr=(sr>>15)+1;} else {xr=sr>>15;}          /* Round X */

      sr = (sv * say) + (zr * cay);       /* xr & zr = point rotated on Y*/

      i = sr & 0x4000;

      if (i != 0) {zr=(sr>>15)+1;} else {zr=sr>>15;}

ROTX:       /* isaz=xr; printi("RZ xr=%",isaz);

            isaz=zr;    printi(" zr=%",isaz); */

      sv = yr;

      sr = (yr * cax) - (zr * sax);       /* "X" axis rotation, given sin*/

                                  /* isaz=sr>>10;  printi(" sr=%",isaz); */

      i = sr & 0x4000;                        /* printi(" sr&4000=%",i); */

      if (i != 0) {yr=(sr>>15)+1;} else {yr=sr>>15;}           /* Round Y*/

            /* isaz=yr; printi(" yr=%",isaz); */

      sr = (sv * sax) + (zr * cax);       /* yr & zr = point rotated on X*/

            /* isaz=sr>>10;   printi(" 2nd sr=%",isaz); */

      i = sr & 0x4000;  /* printi(" sr&4000=%",i); */

      if (i != 0) {zr=(sr>>15)+1;} else {zr=sr>>15;}           /* Round Z*/

            /* isaz=yr; printi("RX yr=%",isaz);

            isaz=zr;    printi(" zr=%",isaz); */

             /* Compute perspective on rotated point given "focal" length*/

ROTP: /* xt=xr>>4;  yt=yr>>4;  zt=zr>>4;  /* Printing test ***************/

        /* printi(" ROT xt=%",xt); printi(" yt=%",yt); printi(" zt=%",zt);

                                      foc=focal; printi(" focal=%",foc); */

      xr = xr + (xoff << 4);         /* Add X offset * 16 to rotated "X" */

      yr = yr + (yoff << 4);

      zr = zr + (zoff << 4);

               /* Prevent: divide by zero CPU error, & showing -Z points */

      if (zr <= 0) {xper = -32768; goto RTN;} /* Flag XPER outside window*/

      xper = ((xr * focal)/zr);      /* "X" perspective on rotated point */

      yper = ((yr * focal)/zr);      /* "Y" perspective on rotated point */

RTN:          /* printi (" PERS xper=%",xper);  printi(" yper=%%",yper); */

      return; }               /* Return x, y from rotate and perspective.*/

/**/


E-mail: info@cyberneticcinema.com, Web site: www.cyberneticcinema.com, © 2003 Cybernetic Cinema, updated Tuesday, May 13, 2003 3:43 AM