/******************************************************************************
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.*/