package stamp.math; import stamp.core.*; /** * This class provides highly accurate sin,cos,tan,arcsin,arccos,arctan and atan2 functions * without using floating point calculations. * Angles are specified in 0.1 degree units, range 0 to 3599. * * @author Peter Verkaik (peterverkaik@boselectro.nl) * @version 1.1 December 21, 2005 */ public class Trig { static int[] sintable = new int[]{ 0, 175, 349, 523, 698, 872, 1045, 1219, 1392, 1564, 1736, 1908, 2079, 2250, 2419, 2588, 2756, 2924, 3090, 3256, 3420, 3584, 3746, 3907, 4067, 4226, 4384, 4540, 4695, 4848, 5000, 5150, 5299, 5446, 5592, 5736, 5878, 6018, 6157, 6293, 6428, 6561, 6691, 6820, 6947, 7071, 7193, 7314, 7431, 7547, 7660, 7771, 7880, 7986, 8090, 8192, 8290, 8387, 8480, 8572, 8660, 8746, 8829, 8910, 8988, 9063, 9135, 9205, 9272, 9336, 9397, 9455, 9511, 9563, 9613, 9659, 9703, 9744, 9781, 9816, 9848, 9877, 9903, 9925, 9945, 9962, 9976, 9986, 9994, 9998,10000}; /* Calculate sin(x) sin(x) is calculated from a table with 91 entries, each entry represents the sine value for a whole degree in the first quadrant. Sine values for fractional degrees are calculated by lineair interpolation between the lower and higher whole degree. */ /** * Calculate sin(x) * * @param x Angle in 0.1 degree units (-3600 to +3600) * @return sin(x) in 0.0001 units (-10000 to +10000) */ public static int sin(int x) { x = (x + 3600)%3600; int q = x%900; switch (x/900) { case 1: return (q == 0) ? 10000 : sin(900-q); case 2: return -sin(q); case 3: return (q == 0) ? -10000 : -sin(900-q); } int index = x/10; int j = sintable[index]; return j + ((sintable[index+1]-j)*(x%10)/10); } /** * Calculate cos(x) * * @param x Angle in 0.1 degree units (-3600 to +3600) * @return cos(x) in 0.0001 units (-10000 to +10000) */ public static int cos(int x) { x = (x + 3600)%3600; int q = x%900; switch (x/900) { case 1: return -sin(q); case 2: return -sin(900-q); case 3: return sin(q); } return sin(900-q); } /** * Calculate tan(x) * * @param x Angle in 0.1 degree units (-3600 to +3600) * @return tan(x) in 0.01 units (-29409 to +29409 for angles -89.8 to +89.8) * angles 89.9,90.0 and -90.1 return 32767, angles -89.9,-90.0 and 90.1 return -32767 */ public static int tan(int x) { x = (x + 3600)%3600; int q = x%900; switch (x/900) { case 1: return (q == 0) ? 32767 : -tan(900-q); case 2: return tan(q); case 3: return (q == 0) ? -32767 : -tan(900-q); } if (q == 899) return 32767; int a = sin(900-q); //cos(q) int b = sin(q); //tan(x)=100*sin(x)/cos(x)=(I + F/65536)*sin(x) return ((100/a)*b) + UnsignedIntMath.umulf(b,UnsignedIntMath.ufrac(100,a)); } /* Calculate arcsin(x) angle = arcsin(x/10000) -10000 <= x <= 10000, 0 <= angle <= 3599 sin(angle) = x/10000 sin(angle)*sin(angle) + cos(angle)*cos(angle) = 1 cos(angle) = sqrt(10000*10000-x*x) / 10000 = 100*sqrt(100*100-(x/100)*(x/100)) / 10000 = 100*sqrt((100+x/100)*(100-x/100)) / 10000 = 100*sqrt(((10000+x)/100)*((10000-x)/100)) / 10000 = 100*sqrt((I+F/65536)*(J+H/65536)) / 10000 = 100*sqrt(I*J + J*F/65536 + I*H/65536) / 10000 tan(angle) = sin(angle)/cos(angle) = x / y angle = atan2(x,y) */ private static int arc(int x) { int I = (10000+x)/100; int F = UnsignedIntMath.ufrac(10000+x,100); int J = (10000-x)/100; int H = UnsignedIntMath.ufrac(10000-x,100); int r = (I*J) + UnsignedIntMath.umulf(I,H) + UnsignedIntMath.umulf(J,F); r = UnsignedIntMath.usqrt(r); return 100*(r>>>8) + UnsignedIntMath.umulf(100,r<<8); } /** * Calculate arcsin(x) * * @param x Sine value in 0.0001 units (-10000 to +10000) * @return Angle in 0.1 degree units (-1799 to +1800) */ public static int arcsin(int x) { return atan2(x,arc(x)); } /* Calculate arccos(x) angle = arccos(x/10000) -10000 <= x <= 10000, 0 <= angle <= 3599 cos(angle) = x/10000 sin(angle)*sin(angle) + cos(angle)*cos(angle) = 1 sin(angle) = sqrt(10000*10000-x*x) / 10000 = 100*sqrt(100*100-(x/100)*(x/100)) / 10000 = 100*sqrt((100+x/100)*(100-x/100)) / 10000 = 100*sqrt(((10000+x)/100)*((10000-x)/100)) / 10000 = 100*sqrt((I+F/65536)*(J+H/65536)) / 10000 = 100*sqrt(I*J + J*F/65536 + I*H/65536) / 10000 tan(angle) = sin(angle)/cos(angle) = y / x angle = atan2(y,x) */ /** * Calculate arccos(x) * * @param x Cosine value in 0.0001 units (-10000 to +10000) * @return Angle in 0.1 degree units (-1799 to +1800) */ public static int arccos(int x) { return atan2(arc(x),x); } /* Calculate arctan(x) angle = arctan(x/100) -32768 <= x <= 32767, 0 <= angle <= 3599 angle = atan2(x,100) */ /** * Calculate arctan(x) * * @param x Tangens value in 0.01 units (-32768 to +32767) * @return Angle in 0.1 degree units (-900 to +900) */ public static int arctan(int x) { return atan2(x,100); } /* arctan approximation (y/x) 1800 atn = ------------------------ * ---- 0.1 deg units, for |y/x|<1 1 + 0.280872*(y/x)*(y/x) pi (y/x) 1800 = ----------------------------- * ---- 0.1 deg units, for |y/x|<1 1 + (18407/65536)*(y/x)*(y/x) pi |y/x|<1, so y/x = F/65536 (F/65536) 1800 atn = ------------------------------------- * ---- 0.1 degree units, for 0<=F<= 65535 1 + (18407/65536)*(F/65536)*(F/65536) pi F 1800 = --------------------------------- * ---- 0.1 degree units, for 0<=F<= 65535 65536 + 18407*(F/65536)*(F/65536) pi F/2 1800 = -------------------------------- * ---- 0.1 degree units, for 0<=F<= 65535 32768 + 9204*(F/65536)*(F/65536) pi 1/pi = 20861/65536 for |y/x|>1 use atn(y/x) = 900 - atn(x/y) */ /** * Calculate arctan(y/x) * * @param y Vertical signed value * @param x Horizontal signed value * @return Angle in 0.1 degree units (-1799 to +1800) */ public static int atan2(int y, int x) { int q = (x<0) ? 1 : 0; //calculate quadrant q ^= (y<0) ? 3 : 0; y = UnsignedIntMath.abs(y); x = UnsignedIntMath.abs(x); if (y == x) { q = q*900 + 450; //this returns 450 in case y=x=0 return (q>1800) ? q-3600 : q; } if (x == 0) return (q==3) ? -900 : 900; boolean inverse = (y > x); int f = (inverse) ? UnsignedIntMath.ufrac(x,y) : UnsignedIntMath.ufrac(y,x); int a = UnsignedIntMath.umulf(9204,f); // 9204*(F/65536) int n = f>>>1; int d = (short)32768 + UnsignedIntMath.umulf(a,f); int r = UnsignedIntMath.umulf(1800,UnsignedIntMath.ufrac(n,d)); // 1800*(n/d) r = UnsignedIntMath.umulf(r,20861); // 1800*(n/d)*(1/pi) if (inverse) r = 900-r; switch (q) { case 1: r = 1800-r; break; case 2: r = 1800+r; break; case 3: r = 3600-r; } return (r>1800) ? r-3600 : r; } }