!    u,v =   variables used to scale a and b for computing r
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) aa
  real ( kind = 8 ) b
  real ( kind = 8 ) bb
  real ( kind = 8 ) c
  real ( kind = 8 ) r
  real ( kind = 8 ) s
  real ( kind = 8 ) u
  real ( kind = 8 ) v

  aa = a
  bb = b
!
!  Abs(b) < abs(a)
!  Note that R has the sign of A, 0 < C, and S has sign(a)*sign(b).
!
  if ( abs ( bb ) < abs ( aa ) ) then
    u = aa + aa
    v = bb / u
    r = sqrt ( 0.25D+00 + v * v ) * u
    c = aa / r
    s = v * ( c + c )
    b = s
    a = r
    return
  end if
!
!  abs(a) <= abs(b)
!
  if ( bb == 0.0D+00 ) then
    c = 1.0D+00
    s = 0.0D+00
    return
  end if

  u = bb + bb
  v = aa / u
!
!  Store R in A.
!
  a = sqrt ( 0.25D+00 + v * v ) * u
  s = bb / a
  c = v * ( s + s )
!
!  Note that R has the sign of b, 0 < S, and c has sign(a)*sign(b).
!
  b = 1.0D+00
  if ( c /= 0.0D+00 ) then
    b = 1.0D+00 / c
  end if

  return
end
subroutine rotate ( n, c, s, x, y )

!***********************************************************************
!
!! ROTATE applies a Givens rotation to two vectors.
!
!  Discussion:
!
!    This routine applies the Givens rotation
!
!      ( c  s)
!      (-s  c)
! 
!    to the 2 by n matrix
!
!      (x(1) ... x(n))
!      (y(1) ... y(n))
!
!  Author:
!
!    Robert Renka,
!    University of North Texas,
!    (817) 565-2767.
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 661: QSHEP3D, Quadratic Shepard method for trivariate
!    interpolation of scattered data,
!    ACM Transactions on Mathematical Software,
!    Volume 14, 1988, pages 151-152.
!
!  Parameters:
!
!    Input, integer N, the number of columns to be rotated.
!
!    Input, real ( kind = 8 ) C, S, the elements of the Givens rotation.
!    These may be determined by subroutine GIVENS.
!
!    Input/output, real ( kind = 8 ) X(N), Y(N), the vectors to be rotated.
!
  implicit none

  integer n

  real ( kind = 8 ) c
  integer i
  real ( kind = 8 ) s
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) xi
  real ( kind = 8 ) y(n)
  real ( kind = 8 ) yi

  if ( c == 1.0D+00 .and. s == 0.0D+00 ) then
    return
  end if

  do i = 1, n
    xi = x(i)
    yi = y(i)
    x(i) =  c * xi + s * yi
    y(i) = -s * xi + c * yi
  end do

  return
end
subroutine setup3 ( xk, yk, zk, fk, xi, yi, zi, fi, s1, s2, r, row )

!***********************************************************************
!
!! SETUP3 sets up the weighted least-squares fit of the data.
!
!  Discussion:
!
!    This routine sets up the I-th row of an augmented regression matrix 
!    for a weighted least-squares fit of a quadratic function Q(X,Y,Z) 
!    to a set of data values F, where Q(XK,YK,ZK) = FK.  
!
!    The first 6 columns (quadratic terms) are scaled by 1/S2, and columns 
!    7, 8, and 9 (linear terms) are scaled by 1/S1.  The weight is 
!    (R-D)/(R*D) if D < R, and 0 if R <= D, where D is the distance 
!    between nodes I and K.
!
!  Author:
!
!    Robert Renka,
!    University of North Texas,
!    (817) 565-2767.
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 661: QSHEP3D, Quadratic Shepard method for trivariate
!    interpolation of scattered data,
!    ACM Transactions on Mathematical Software,
!    Volume 14, 1988, pages 151-152.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) XK, YK, ZK, FK = coordinates and data value
!    at node K (interpolated by q).
!
!    Input, real ( kind = 8 ) XI, YI, ZI, FI = coordinates and data value
!    at node I.
!
!    Input, real ( kind = 8 ) S1, S2 = reciprocals of the scale factors.
!
!    Input, real ( kind = 8 ) R = radius of influence about node K defining the
!    weight.
!
!    Output, real ( kind = 8 ) ROW(10), a row of the augmented
!    regression matrix.
!
!  Local parameters:
!
!    d =   distance between nodes k and i
!    w =   weight associated with the row
!    w1 =   w/s1
!    w2 =   w/s2
!
  implicit none

  real ( kind = 8 ) d
  real ( kind = 8 ) dx
  real ( kind = 8 ) dxsq
  real ( kind = 8 ) dy
  real ( kind = 8 ) dysq
  real ( kind = 8 ) dz
  real ( kind = 8 ) dzsq
  real ( kind = 8 ) fi
  real ( kind = 8 ) fk
  integer i
  real ( kind = 8 ) r
  real ( kind = 8 ) row(10)
  real ( kind = 8 ) s1
  real ( kind = 8 ) s2
  real ( kind = 8 ) w
  real ( kind = 8 ) w1
  real ( kind = 8 ) w2
  real ( kind = 8 ) xi
  real ( kind = 8 ) xk
  real ( kind = 8 ) yi
  real ( kind = 8 ) yk
  real ( kind = 8 ) zi
  real ( kind = 8 ) zk

  dx = xi - xk
  dy = yi - yk
  dz = zi - zk

  d = sqrt ( dx**2 + dy**2 + dz**2 )

  if ( d <= 0.0D+00 .or. r <= d ) then
    row(1:10) = 0.0D+00
    return
  end if

  w = ( r - d ) / r / d
  w1 = w / s1
  w2 = w / s2

  row(1) = dx * dx * w2
  row(2) = dx * dy * w2
  row(3) = dy * dy * w2
  row(4) = dx * dz * w2
  row(5) = dy * dz * w2
  row(6) = dz * dz * w2

  row(7) = dx * w1
  row(8) = dy * w1
  row(9) = dz * w1

  row(10) = ( fi - fk ) * w

  return
end
subroutine store3 ( n, x, y, z, nr, lcell, lnext, xyzmin, xyzdel, ier )

!***********************************************************************
!
!! STORE3 sets up a data structure for N scattered nodes in 3D.
!
!  Discussion:
!
!    Given a set of N arbitrarily distributed nodes in three-space, 
!    this subroutine creates a data structure for a cell-based method of 
!    solving closest-point problems.  The smallest box containing the nodes 
!    is partitioned into an NR by NR by NR uniform grid of cells, and 
!    nodes are associated with cells.  In particular, the data structure
!    stores the indices of the nodes contained in each cell.  For a 
!    uniform random distribution of nodes, the nearest node to an 
!    arbitrary point can be determined in constant expected time.
!
!  Author:
!
!    Robert Renka,
!    University of North Texas,
!    (817) 565-2767.
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 661: QSHEP3D, Quadratic Shepard method for trivariate
!    interpolation of scattered data,
!    ACM Transactions on Mathematical Software,
!    Volume 14, 1988, pages 151-152.
!
!  Parameters:
!
!    Input, integer N, the number of nodes.  2 <= N.
!
!    Input, real ( kind = 8 ) X(N), Y(N), Z(N), the coordinates of the nodes.
!
!    Input, integer NR, the number of rows, columns, and planes in the
!    grid.  The cell density (average number of nodes per cell) is 
!    D = N/(NR**3).  A recommended value, based on empirical evidence, 
!    is D = 3, so NR = (N/3)**(1/3).  1 <= NR.
!
!    Output, integer LCELL(NR,NR,NR), a cell array such that LCELL(I,J,K)
!    contains the index (for X, Y, and Z) of the first node (node with 
!    smallest index) in cell (I,J,K), or LCELL(I,J,K) = 0 if no nodes are
!    contained in the cell.  The orner of cell (I,J,K) which is farthest
!    from the box corner defined by XYZMIN has coordinates
!    (XMIN+I*DX,YMIN+J*DY,ZMIN+K*DZ), where (XMIN,YMIN,ZMIN) are the 
!    elements of XYZMIN.  LCELL is not defined if IER is returned nonzero.
!
!    Output, integer LNEXT(N),  next-node indices such that
!    lnext(l) contains the index of the next node
!    in the cell which contains node l, or
!    lnext(l) = l if l is the last node in the
!    cell for l = 1,...,n.  (the nodes contained
!    in a cell are ordered by their indices.)
!    if, for example, cell (i,j,k) contains nodes
!    2, 3, and 5 (and no others), then
!    lcell(i,j,k) = 2, lnext(2) = 3, lnext(3) =
!    5, and lnext(5) = 5.  lnext is not defined
!    if ier /= 0.
!
!    Output, real ( kind = 8 ) XYZMIN(3), the minimum
!    nodal coordinates xmin, ymin, and zmin (in
!    that order) unless ier = 1.  the opposite
!    corner of the box defined by the nodes is
!    (xmin+nr*dx,ymin+nr*dy,zmin+nr*dz).
!
!    Output, real ( kind = 8 ) XYZDEL(3), the dimensions
!    of the cells unless ier = 1.  xyzdel(1) =
!    (xmax-xmin)/nr, xyzdel(2) = (ymax-ymin)/nr,
!    and xyzdel(3) = (zmax-zmin)/nr, where xmin,
!    xmax, ymin, ymax, zmin, and zmax are the
!    extrema of x, y, and z.
!
!    Output, integer IER,  = error indicator --
!    0, if no errors were encountered.
!    1, if n < 2 or nr < 1.
!    2, if a component of xyzdel is not positive.
!
  implicit none

  integer n
  integer nr

  real ( kind = 8 ) delx
  real ( kind = 8 ) dely
  real ( kind = 8 ) delz
  integer i
  integer ier
  integer j
  integer k
  integer l
  integer lb
  integer lcell(nr,nr,nr)
  integer ll
  integer lnext(n)
  integer nn
  integer nnr
  integer np1
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) xmn
  real ( kind = 8 ) xmx
  real ( kind = 8 ) xyzdel(3)
  real ( kind = 8 ) xyzmin(3)
  real ( kind = 8 ) y(n)
  real ( kind = 8 ) ymx
  real ( kind = 8 ) ymn
  real ( kind = 8 ) z(n)
  real ( kind = 8 ) zmx
  real ( kind = 8 ) zmn

  ier = 0
  nn = n
  nnr = nr

  if ( nn < 2 .or. nnr < 1 ) then
    ier = 1
    return
  end if
!
!  Compute the dimensions of the box containing the nodes.
!
  xmn = minval ( x(1:nn) )
  xmx = maxval ( x(1:nn) )
  ymn = minval ( y(1:nn) )
  ymx = maxval ( y(1:nn) )
  zmn = minval ( z(1:nn) )
  zmx = maxval ( z(1:nn) )

  xyzmin(1) = xmn
  xyzmin(2) = ymn
  xyzmin(3) = zmn
!
!  Compute cell dimensions and test for zero area.
!
  delx = ( xmx - xmn ) / real ( nnr, kind = 8 )
  dely = ( ymx - ymn ) / real ( nnr, kind = 8 )
  delz = ( zmx - zmn ) / real ( nnr, kind = 8 )

  xyzdel(1) = delx
  xyzdel(2) = dely
  xyzdel(3) = delz

  if ( delx == 0.0D+00 .or. dely == 0.0D+00 .or. delz == 0.0D+00 ) then
    ier = 2
    return
  end if
!
!  Initialize LCELL.
!
  lcell(1:nnr,1:nnr,1:nnr) = 0
!
!  Loop on nodes, storing indices in LCELL and LNEXT.
!
  np1 = nn + 1

  do ll = 1, nn

    lb = np1 - ll

    i = int ( ( x(lb) - xmn ) / delx ) + 1
    if ( nnr < i ) then
      i = nnr
    end if

    j = int ( ( y(lb) - ymn ) / dely ) + 1
    if ( nnr < j ) then
      j = nnr
    end if

    k = int ( ( z(lb) - zmn ) / delz ) + 1
    if ( nnr < k ) then
      k = nnr
    end if

    l = lcell(i,j,k)

    if ( l == 0 ) then
      lnext(lb) = lb
    else
      lnext(lb) = l
    end if

    lcell(i,j,k) = lb

  end do

  return
end
subroutine timestamp ( )

!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone

  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end