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