subroutine aptrefl (pu, pv, pw, au, av, bu, bv, np, tol,
     &                    qu, qv, qw, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTREFL
c
c     call aptrefl (pu, pv, pw, au, av, bu, bv, np, tol,
c    &              qu, qv, qw, nerr)
c
c     Version:  aptrefl  Updated    1990 November 29 15:30.
c               aptrefl  Originated 1990 January 10 13:50.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  To find, for each of np sets of input data, the vector
c               q = (qu, qv, qw) resulting from the reflection of vector
c               p = (pu, pv, pw) from the surface perpendicular to the uv plane
c               and through the points a = (au, av, 0), and b = (bu, bv, 0),
c               when incident on the surface on the line "ab".
c               When point "a" coincides with point "b", based on tol,
c               no reflection takes place, and vector "q" equals vector "p".
c               Coordinates u, v and w may be any orthogonal coordinates.
c               Flag nerr indicates any input error.
c
c     Method:   The vector "sn" normal to the reflecting surface is the cross
c               product of the vector normal to the uv plane, (0, 0, 1), and
c               the vector parallel to the line "ab".  From simple geometry,
c               "q" = "p" - 2.0 * ("p" dot "u") * "u", where "u" is the unit
c               vector parallel to "sn".
c
c     Input:    pu, pv, pw, au, av, bu, bv, np, tol.
c
c     Output:   qu, qv, qw, nerr.
c
c     Glossary:
c
c     au, av    Input    The u and v coordinates of point "a" in the uv plane.
c                          Must differ from "b", based on tol.  Size np.
c
c     bu, bv    Input    The u and v coordinates of point "b" in the uv plane.
c                          Must differ from "a", based on tol.  Size np.
c
c     nerr      Output   Indicates an input error, if not 0.
c                          1 if np is not positive.
c
c     np        Input    Size of arrays.
c
c     pu,pv,pw  Input    The u, v, w components of incident vector "p".
c                          Size np.
c
c     qu,qv,qw  Output   The u, v, w components of reflected vector "q".
c                          Components of "q" less than the estimated error in
c                          their calculation, based on tol, will be truncated
c                          to zero.
c                          Note that qw = pw.  Size np.
c
c     tol       Input    Numerical tolerance limit.
c                          On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Dimensioned arguments.

c---- Coordinate u of point "a".
      dimension au      (1)
c---- Coordinate v of point "a".
      dimension av      (1)
c---- Coordinate u of point "b".
      dimension bu      (1)
c---- Coordinate v of point "b".
      dimension bv      (1)
c---- Component u of incident vector "p".
      dimension pu      (1)
c---- Component v of incident vector "p".
      dimension pv      (1)
c---- Component w of incident vector "p".
      dimension pw      (1)
c---- Component u of reflected vector "q".
      dimension qu      (1)
c---- Component v of reflected vector "q".
      dimension qv      (1)
c---- Component w of reflected vector "q".
      dimension qw      (1)

c.... Local variables.

c---- Ratio of magnitudes of "p", "q".
      common /laptrefl/ fact
c---- Factor 2.0 * pdots / vlen2.
      common /laptrefl/ fan

c---- A very small number.
      common /laptrefl/ fuz

c---- First index of subset of data.
      common /laptrefl/ n1
c---- Last index of subset of data.
      common /laptrefl/ n2
c---- Index in external array.
      common /laptrefl/ nn
c---- Size of current subset of data.
      common /laptrefl/ ns
c---- Dot product of "p" and "sn".
      common /laptrefl/ pdots
c---- Estimated error in pdots.
      common /laptrefl/ pdserr
c---- Estimated error in qu.
      common /laptrefl/ querr
c---- Estimated error in qv.
      common /laptrefl/ qverr
c---- Component u of normal vector.
      common /laptrefl/ snu
c---- Estimated error in snu.
      common /laptrefl/ snuerr
c---- Component v of normal vector.
      common /laptrefl/ snv
c---- Estimated error in snv.
      common /laptrefl/ snverr
c---- Square of magnitude of vector "sn".
      common /laptrefl/ vlen2
c---- Square of magnitude of vector "p".
      common /laptrefl/ vlenp2
c---- Square of magnitude of vector "q".
      common /laptrefl/ vlenq2
cbugc***DEBUG begins.
cbugc---- Index in arrays.
cbug      common /laptrefl/ n
cbug 9901 format (/ 'aptrefl finding reflection of the vector',
cbug     &  ' (tol=',1pe13.5,'):' /
cbug     &  (i3,' pu,pv,pw=',1p3e22.14 /
cbug     &  '    from the line (in the uv plane) through:' /
cbug     &  '    au,av=   ',1p2e22.14 /
cbug     &  '    bu,bv=   ',1p2e22.14))
cbug      write ( 3, 9901) tol, (n, pu(n), pv(n), pw(n),
cbug     &  au(n), av(n), bu(n), bv(n), n = 1, np)
cbugc***DEBUG ends.

c.... Initialize.

c---- A very small number.
      fuz = 1.e-99

      nerr = 0

c.... Test for input errors.

      if (np .le. 0) then
        nerr = 1
        go to 210
      endif

c.... Set up the indices of the first subset of data.

      n1 = 1
      n2 = min (np, 64)

  110 ns = n2 - n1 + 1

c.... Find the reflection of incident vector "p" in the plane with normal
c....   vector "sn", which is perpendicular to line "ab" and the uv plane.

c---- No truncation tests.
      if (tol .le. 0.0) then

c---- Loop over subset of data.
        do 120 nn = n1, n2

c....     Find the vector "sn" normal to the reflecting surface.

          snu    = av(nn) - bv(nn)
          snv    = bu(nn) - au(nn)
          vlen2  = snu**2 + snv**2

          pdots  = pu(nn) * snu + pv(nn) * snv

c....     Find the reflected vector "q".

          qu(nn) = pu(nn) - 2.0 * pdots * snu / (vlen2 + fuz)
          qv(nn) = pv(nn) - 2.0 * pdots * snv / (vlen2 + fuz)
          qw(nn) = pw(nn)

c---- End of loop over subset of data.
  120   continue

c---- Test for numerical truncation.
      else

c---- Loop over subset of data.
        do 130 nn = n1, n2

c....     Find the vector "sn" normal to the reflecting surface.

          snu    = av(nn) - bv(nn)
          snuerr = tol * (abs (av(nn)) + abs (bv(nn)))
          if (abs (snu) .lt. snuerr) then
            snu = 0.0
          endif

          snv    = bu(nn) - au(nn)
          snverr = tol * (abs (av(nn)) + abs (bv(nn)))
          if (abs (snv) .lt. snverr) then
            snv = 0.0
          endif

          vlen2  = snu**2 + snv**2

          pdots  = pu(nn) * snu + pv(nn) * snv
          pdserr = tol * (abs (pu(nn) * snu) + abs (pv(nn) * snv))
          if (abs (pdots) .lt. pdserr) then
            pdots = 0.0
          endif

          fan    = 2.0 * pdots / (vlen2 + fuz)

c....     Find the reflected vector "q".

          qu(nn) = pu(nn) - fan * snu
          querr  = tol * (abs (pu(nn)) + abs (fan * snu))
          if (abs (qu(nn)) .lt. querr) then
            qu(nn) = 0.0
          endif

          qv(nn) = pv(nn) - fan * snv
          qverr  = tol * (abs (pv(nn)) + abs (fan * snv))
          if (abs (qv(nn)) .lt. qverr) then
            qv(nn) = 0.0
          endif

          qw(nn) = pw(nn)

c....     Make the magnitude of vector "q" the same as vector "p".

          vlenp2 = pu(nn)**2 + pv(nn)**2 + pw(nn)**2
          vlenq2 = qu(nn)**2 + qv(nn)**2 + qw(nn)**2
          fact   = sqrt ((vlenp2 + fuz) / (vlenq2 + fuz))
          qu(nn) = fact * qu(nn)
          qv(nn) = fact * qv(nn)
          qw(nn) = fact * qw(nn)

c---- End of loop over subset of data.
  130   continue

c---- Tested tol.
      endif

c.... See if all data subsets are done.

c---- Do another subset of data.
      if (n2 .lt. np) then
        n1 = n2 + 1
        n2 = min (np, n1 + 63)
        go to 110
      endif
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptrefl results:' /
cbug     &  (i3,' qu,qv,qw=',1p3e22.14))
cbug      write ( 3, 9902) (n, qu(n), qv(n), qw(n), n = 1, np)
cbugc***DEBUG ends.

  210 return

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

UCRL-WEB-209832