subroutine aptrotv (ax, ay, az, bx, by, bz, tol, rotm, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTROTV
c
c     call aptrotv (ax, ay, az, bx, by, bz, tol, rotm, nerr)
c
c     Version:  aptrotv  Updated    1991 July 25 14:30.
c               aptrotv  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", for rotating
c               the vector a = (ax, ay, az) to be parallel to the vector
c               b = (bx, by, bz), around an axis perpendicular to both.
c               Any components of "rotm" within tol of -1.0, 0.0, or 1.0,
c               will be truncated to those values.
c               Flag nerr indicates any input error.
c
c     History:  1990 March 14.  Changed tol to 0.0 in call to unit vector
c               module aptvunb.  Allows small magnitudes.
c               Protected costh from absolute values exceeding 1.0.
c               1991 July 25.  Added coding to find an arbitrary rotation
c               axis perpendicular to vectors "a" and "b" when they are
c               parallel or antiparallel.  Previously, the resulting operator
c               "rotm" was an inversion operator, not a rotation operator.
c
c     Input:    ax, ay, az, bx, by, bz, tol.
c
c     Output:   rotm, nerr.
c
c     Calls: aptvxun, aptvdot, aptvunb 
c               
c
c
c     Glossary:
c
c     ax,ay,az  Input    The x, y, z components of vector "a".
c
c     bx,by,bz  Input    The x, y, z components of vector "b".
c
c     nerr      Output   Indicates an input error, if not 0.
c                          1 if the magnitude of vector "a" or "b" is too small.
c
c     rotm      Output   Rotation operator (a unitary 3 x 3 matrix).
c                          Must be sized rotm(3,3).
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---- Cosine of rotation angle.
      common /laptrotv/ costh
c---- 1.0 - cosine of rotation angle.
      common /laptrotv/ costhm
c---- Index of matrix row.
      common /laptrotv/ i
c---- Index of matrix column.
      common /laptrotv/ j
c---- Sine of rotation angle.
      common /laptrotv/ sinth
c---- Component x of vector "ua".
      common /laptrotv/ uax
c---- Component y of vector "ua".
      common /laptrotv/ uay
c---- Component z of vector "ua".
      common /laptrotv/ uaz
c---- Component x of vector "ub".
      common /laptrotv/ ubx
c---- Component y of vector "ub".
      common /laptrotv/ uby
c---- Component z of vector "ub".
      common /laptrotv/ ubz
c---- Component x of unit axis vector.
      common /laptrotv/ ux
c---- Component y of unit axis vector.
      common /laptrotv/ uy
c---- Component z of unit axis vector.
      common /laptrotv/ uz
c---- Magnitude of a vector.
      common /laptrotv/ vlen
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptrotv.  tol=',1pe13.5,'  Rotating vector' /
cbug     &  '  ax,ay,az=',1p3e22.14 /
cbug     &  '  around mutually perpendicular axis onto vector' /
cbug     &  '  bx,by,bz=',1p3e22.14)
cbug      write ( 3, 9901) tol, ax, ay, az, bx, by, bz
cbugc***DEBUG ends.

c.... Initialize.

      nerr = 0

c.... Find the unit vector parallel to "a".

      call aptvunb (ax, ay, az, 1, 0., uax, uay, uaz, vlen, nerr)

c---- Vector magnitude too small.
      if (vlen .le. tol) then
        nerr = 1
        go to 210
      endif

c.... Find the unit vector parallel to "b".

      call aptvunb (bx, by, bz, 1, 0., ubx, uby, ubz, vlen, nerr)

c---- Vector magnitude too small.
      if (vlen .le. tol) then
        nerr = 1
        go to 210
      endif

c.... Find the unit vector parallel to the axis of rotation.

      call aptvxun (uax, uay, uaz, ubx, uby, ubz, 1, tol,
     &              ux, uy, uz, vlen, nerr)

c.... Find an arbitrary rotation axis if "a" and "b" are parallel
c....   or antiparallel.

c---- Vectors "a" and "b" aligned.
      if (vlen .eq. 0.0) then
        if (uax .ne. 0.0) then
          vlen = sqrt (uax**2 + uay**2)
          ux   = -uay / vlen
          uy   =  uax / vlen
        elseif (uay .ne. 0.0) then
          vlen = sqrt (uay**2 + uaz**2)
          uy   = -uaz / vlen
          uz   =  uay / vlen
        elseif (uaz .ne. 0.0) then
          vlen = sqrt (uaz**2 + uax**2)
          uz   = -uax / vlen
          ux   =  uaz / vlen
        endif
c---- Tested rotation axis vector.
      endif

c.... Find the needed functions of the angle of rotation.
c....   costh = ua .dot. ub (scalar product).

      call aptvdot (uax, uay, uaz, ubx, uby, ubz, 1, tol, costh, nerr)

      costh = amax1 (-1.0, amin1 (costh, 1.0))

      costhm = 1.0 - costh
      sinth  = sqrt (costhm * (1.0 + costh))

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

      rotm(1,1) = costhm * ux * ux + costh
      rotm(1,2) = costhm * ux * uy - sinth * uz
      rotm(1,3) = costhm * ux * uz + sinth * uy
      rotm(2,1) = costhm * uy * ux + sinth * uz
      rotm(2,2) = costhm * uy * uy + costh
      rotm(2,3) = costhm * uy * uz - sinth * ux
      rotm(3,1) = costhm * uz * ux - sinth * uy
      rotm(3,2) = costhm * uz * uy + sinth * ux
      rotm(3,3) = costhm * uz * uz + costh

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 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptrotv results:  rotm=',3(/ 11x,1p3e22.14))
cbug      write ( 3, 9902) ((rotm(i,j), j = 1, 3), i = 1, 3)
cbugc***DEBUG ends.
      return

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

UCRL-WEB-209832