1 /* @(#)atan.c 4.1 12/25/82 */ 2 3 4 /* 5 floating-point arctangent 6 7 atan returns the value of the arctangent of its 8 argument in the range [-pi/2,pi/2]. 9 10 atan2 returns the arctangent of arg1/arg2 11 in the range [-pi,pi]. 12 13 there are no error returns. 14 15 coefficients are #5077 from Hart & Cheney. (19.56D) 16 */ 17 18 19 double static sq2p1 =2.414213562373095048802e0; 20 static double sq2m1 = .414213562373095048802e0; 21 static double pio2 =1.570796326794896619231e0; 22 static double pio4 = .785398163397448309615e0; 23 static double p4 = .161536412982230228262e2; 24 static double p3 = .26842548195503973794141e3; 25 static double p2 = .11530293515404850115428136e4; 26 static double p1 = .178040631643319697105464587e4; 27 static double p0 = .89678597403663861959987488e3; 28 static double q4 = .5895697050844462222791e2; 29 static double q3 = .536265374031215315104235e3; 30 static double q2 = .16667838148816337184521798e4; 31 static double q1 = .207933497444540981287275926e4; 32 static double q0 = .89678597403663861962481162e3; 33 34 35 /* 36 atan makes its argument positive and 37 calls the inner routine satan. 38 */ 39 40 double 41 atan(arg) 42 double arg; 43 { 44 double satan(); 45 46 if(arg>0) 47 return(satan(arg)); 48 else 49 return(-satan(-arg)); 50 } 51 52 53 /* 54 atan2 discovers what quadrant the angle 55 is in and calls atan. 56 */ 57 58 double 59 atan2(arg1,arg2) 60 double arg1,arg2; 61 { 62 double satan(); 63 64 if((arg1+arg2)==arg1) 65 if(arg1 >= 0.) return(pio2); 66 else return(-pio2); 67 else if(arg2 <0.) 68 if(arg1 >= 0.) 69 return(pio2+pio2 - satan(-arg1/arg2)); 70 else 71 return(-pio2-pio2 + satan(arg1/arg2)); 72 else if(arg1>0) 73 return(satan(arg1/arg2)); 74 else 75 return(-satan(-arg1/arg2)); 76 } 77 78 /* 79 satan reduces its argument (known to be positive) 80 to the range [0,0.414...] and calls xatan. 81 */ 82 83 static double 84 satan(arg) 85 double arg; 86 { 87 double xatan(); 88 89 if(arg < sq2m1) 90 return(xatan(arg)); 91 else if(arg > sq2p1) 92 return(pio2 - xatan(1.0/arg)); 93 else 94 return(pio4 + xatan((arg-1.0)/(arg+1.0))); 95 } 96 97 /* 98 xatan evaluates a series valid in the 99 range [-0.414...,+0.414...]. 100 */ 101 102 static double 103 xatan(arg) 104 double arg; 105 { 106 double argsq; 107 double value; 108 109 argsq = arg*arg; 110 value = ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0); 111 value = value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0); 112 return(value*arg); 113 } 114