subroutine aptrots (n1, th1, n2, th2, n3, th3, ku, tol,
     &                    rotm, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTROTS
c
c     call aptrots (n1, th1, n2, th2, n3, th3, ku, tol,
c    &              rotm, nerr)
c
c     Version:  aptrots  Updated    1989 November 13 15:20.
c               aptrots  Originated 1989 November 2 14:10.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  To find the rotation matrix operator rotm, to do a sequential
c               rotation of angle th1 around axis n1, angle th2 around axis n2,
c               and angle th3 around axis n3.  All angles are measured
c               counterclockwise, with the axis pointed at the observer.
c               The axes may be x (1), y (2), or z (3).
c               Angles may be in degrees (ku = 0) or radians (ku = 1).
c               Flag nerr indicates any input error.
c
c     Input:    n1, th1, n2, th2, n3, th3, ku, tol.
c
c     Output:   rotm, nerr.
c
c     Glossary:
c
c     ku        Input    Indicates theta units are degrees (0) or radians (1).
c
c     n1        Input    Indicates first axis is x (1), y (2), or z (3).
c                          May not be 0, but th1 may be 0.
c
c     n2        Input    Indicates second axis is x (1), y (2), or z (3).
c                          May not be 0, but th2 may be 0.
c
c     n3        Input    Indicates third axis is x (1), y (2), or z (3).
c                          May not be 0, but th3 may be 0.
c
c     nerr      Output   Indicates an input error, if not 0.
c                          1 if n1, n2, or n3 not in range 1-3, or not unique.
c                          2 if ku is not 0 or 1.
c
c     rotm      Output   Rotation operator (a unitary 3 x 3 matrix).
c                          Must be sized rotm(3,3).
c
c     th1       Input    Angle of rotation around axis n1, counterclockwise
c                          when rotation axis is pointed at observer.  Units
c                          are degrees (ku = 0) or radians (ku = 1).
c
c     th2       Input    Angle of rotation around axis n2.  See th1.
c
c     th3       Input    Angle of rotation around axis n3.  See th1.
c
c     tol       Input    Numerical tolerance limit.  Used to test and adjust
c                          unit vector components and point coordinates.
c                          On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Dimensioned arguments.

c---- Rotation matrix operator.
      dimension rotm    (3,3)

c.... Local variables.

c---- Rotation angle around axis.
      common /laptrots/ ang     (3)
c---- Cosine of rotation angle ang.
      common /laptrots/ costh   (3)
c---- Index of matrix row.
      common /laptrots/ i
c---- Index of matrix column.
      common /laptrots/ j
c---- Parity of axis sequence.
      common /laptrots/ npar
c---- Sum of axis numbers.
      common /laptrots/ nsum
c---- Parity of axis sequence.
      common /laptrots/ par

c---- Mathematical constant, pi.
      common /laptrots/ pi

c---- Sine of rotation angle ang.
      common /laptrots/ sinth(3)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptrots.  Rotate sequentially around axis by angle:' /
cbug     &   5x,3(i2,f13.6))
cbug      write ( 3, 9901) n1, th1, n2, th2, n3, th3
cbugc***DEBUG ends.

c.... Initialize.

c---- Mathematical constant, pi.
      pi = 3.14159265358979323

      nerr = 0

c.... Find parity of sequence of axes.

      nsum = n1 + n2 + n3

c---- Bad axis numbers.
      if (nsum .ne. 6) then
        nerr = 1
        go to 210
      endif

      npar = (n2 - n1) * (n3 - n2) * (n3 - n1)

c---- 1,2,3 or 2,3,1 or 3,1,2.
      if (npar .eq. 2) then
        par = 1.0
c---- 3,2,1 or 2,1,3 or 1,3,2.
      elseif (npar .eq. -2) then
        par = -1.0
      else
        nerr = 1
        go to 210
      endif

c.... Get rotation angles in radians, and the sines and cosines.

      if (ku .eq. 0) then
        ang(n1) = th1 * pi / 180.0
        ang(n2) = th2 * pi / 180.0
        ang(n3) = th3 * pi / 180.0
      elseif (ku .eq. 1) then
        ang(n1) = th1
        ang(n2) = th2
        ang(n3) = th3
      else
        nerr = 2
        go to 210
      endif

      costh(n1) = cos (ang(n1))
      costh(n2) = cos (ang(n2))
      costh(n3) = cos (ang(n3))
      sinth(n1) = sin (ang(n1))
      sinth(n2) = sin (ang(n2))
      sinth(n3) = sin (ang(n3))

c.... Form the components of the rotation matrix.

      rotm(n1,n1) =        costh(n2) * costh(n3)
      rotm(n1,n2) =        sinth(n1) * sinth(n2) * costh(n3) -
     &               par * costh(n1) * sinth(n3)
      rotm(n1,n3) =        sinth(n1) * sinth(n3) +
     &               par * costh(n1) * sinth(n2) * costh(n3)
      rotm(n2,n1) =  par * costh(n2) * sinth(n3)
      rotm(n2,n2) =        costh(n1) * costh(n3) +
     &               par * sinth(n1) * sinth(n2) * sinth(n3)
      rotm(n2,n3) =        costh(n1) * sinth(n2) * sinth(n3) -
     &               par * sinth(n1) * costh(n3)
      rotm(n3,n1) = -par * sinth(n2)
      rotm(n3,n2) =  par * sinth(n1) * costh(n2)
      rotm(n3,n3) =        costh(n1) * costh(n2)

c---- Adjust components near 0 or 1.
      if (tol .gt. 0.0) then

        do 120 i = 1, 3
          do 110 j = 1, 3
            if (abs (rotm(i,j)) .le. tol) then
              rotm(i,j) = 0.0
            elseif ((abs (abs (rotm(i,j)) - 1.0)) .le. tol) then
              rotm(i,j) = sign (1.0, rotm(i,j))
            endif
  110     continue
  120   continue

c---- Tested tol.
      endif

  210 return

c.... End of subroutine aptrots.      (+1 line.)
      end

UCRL-WEB-209832