Equations that involve trigonometric functions can have infinitely many solutions. For example, the solution to the equation tan(θ)=1 is θ = π/4 + kπ, where k is any integer. In order to obtain a unique solution to the equation, we define the "arc" functions: inverse trigonometric functions that return a principal value. For example, the Arctan function returns a principal value in the range (-π/2, π/2), such as Arctan(1) = π/4. In SAS software, the ATAN function computes the Arctan function.

You can extend the Arctan function in a natural way: given the coordinates (x, y) of a point in the plane, you can compute the angle formed between the x axis and the vector from the origin to (x, y). SAS provides the ATAN2 function for this computation. The ATAN2 function takes two arguments (the coordinates of an (x, y) pair) and returns the angle in the range (-π, π] that corresponds to the vector angle.

The ATAN and ATAN2 functions are supported in Base SAS, so you can call them from the DATA step or from the SAS/IML language. The following statements compute the points on the unit circle for several polar angles. The ATAN function evaluated at y/x returns the principal arctangent function. The ATAN2 function evaluated at (y, x) returns the polar angle in (-π, π]. Notice that order of the arguments for the ATAN2 function is the reverse of what you might expect!

data Arctan;
drop pi i theta;
array a[8] $6 _temporary_ 
           ("0" "pi/4" "pi/2" "3*pi/4" "pi" "5*pi/4" "3*pi/2" "7*pi/4") ;
pi = constant("pi");
do i = 0 to 7;
   Angle = a[i+1];
   theta = i*pi/4;        /* angle in [0, 2*pi) */
   x = cos(theta);        /* (x,y) on unit circle */
   y = sin(theta);
   atan = atan( y/x );    /* ATAN takes one argument */
   atan2 = atan2(y, x);   /* ATAN2 takes two arguments: Notice order! */
   output;
end;
run;
proc print; run;

Notice that the ATAN and ATAN2 functions agree for points that are in the first and fourth quadrants. Notice also that both functions return their values in radians. You can convert radians to degrees by multiplying the radian measure by 180/π.

Compute angles along polar curves

Suppose that you have data points such as are shown in the scatter plot at the left. The ATAN2 function enables you to compute an angle in (-π, π] for a given data point. But to reconstruct the curve, you need to compute angles outside that interval. The scatter plot shows points that lie along a logarithmic spiral as computed by the following program.

proc iml;
pi = constant("pi");
/* create (x,y) coordinates of polar curve for (golden) logarithmic spiral 
   r(theta) = a*exp(b*theta)  */   
a = 1;   b = -0.3063;
theta = do(0, 6*pi, pi/24)`;
x = a*cos(theta)#exp(b*theta);
y = a*sin(theta)#exp(b*theta);
title "Logarithmic Spiral";
call scatter(x,y) other="refline 0 / axis=y; refline 0 / axis=x;";

Recently I had data that looked like the scatter plot, and I needed to recover the angles that correspond to the points along the curve. If you naively use the ATAN2 function, all angles are translated into the interval (-π, π]. However, if the data are ordered according to the polar angle, as in these data, you can monitor the output of the ATAN2 function as you sequentially traverse the points along the curve. When one point is in the second quadrant but the next point is in the third quadrant, the ATAN2 function "jumps" by -2π across those two points. You can detect the discontinuity and add 2π to the value of the ATAN2 function evaluated all subsequent points along the curve. If you do this every time that the curve crosses the negative x axis, the resulting set of angles will correctly represent the polar angle along the curve.

This algorithm is implemented in the following SAS/IML module, which assumes that the polar angle increases along the curve. With slightly more complicated logic, you can derive the polar angle for points along any parametric curve in the plane.

/* assume the polar angle is an increasing function along the curve */
start IncrAngle(x, y);
   theta = atan2(y,x);            /* angle from center to point */
   dif = dif(theta);              /* difference of angles between adjacent pts */
   idx = loc(dif^=. & dif<0);     /* for which pts does angle decrease? */
   n = nrow(dif);
   pi = constant("pi");
   do i = 1 to ncol(idx);
      theta[idx[i]:n] = 2*pi + theta[idx[i]:n]; /* add 2*pi at discontinuities */
   end;
   return(theta);
finish;
 
phi = IncrAngle(x,y);         /* test function for pts along logarithmic curve */

The variable phi is a vector of increasing values. It is equal to the theta vector, which contained the original sequence of angles in [0, 6π].

Summary

For any value, the ATAN function computes the principal arctangent of that value. For a point (x, y) in the plane, the ATAN2 function enables you to compute the angle formed between the x axis and the vector from the origin to (x, y). The ATAN2 function returns an angle in the range (-π, π], and the ATAN2 function has a jump discontinuity of 2π across the negative x axis. Using these facts, you can construct a continuous function that returns the polar angle for points along a parametric curve in the plane.

In my next blog post I will use this function to analyze data that seems to lie along a parametric curve in the plane.