program nms_prb
!
!*******************************************************************************
!
!! NMS_PRB calls the NMS tests.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'NMS_PRB'
  write ( *, '(a)' ) '  Run the NMS tests.'

  call test001
  call test002
  call test003
  call test004
  call test005
  call test006
  call test007
  call test008
  call test009
  call test010

  call test011
  call test012
  call test013
  call test014
  call test015
  call test016
  call test017
  call test018
  call test019
  call test020

  call test021
! call test022
  call test023
  call test024
  call test025
  call test026
  call test027
  call test028
  call test029
  call test030

  call test031
  call test032
  call test033
  call test034
  call test035
  call test036

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'NMS_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test001
!
!*******************************************************************************
!
!! TEST001 tests ALNGAM.
!
  implicit none
!
  real alngam
  real fx
  real fx2
  integer n
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST001:'
  write ( *, '(a)' ) '  ALNGAM evaluates the log of the Gamma function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X       Exact F       ALNGAM(X)'
  write ( *, '(a)' ) ' '

  n = 0

  do

    call gamma_values ( n, x, fx )

    if ( n == 0 ) then
      exit
    end if

    fx = log ( fx )
    fx2 = alngam ( x )

    write ( *, '(f8.4,2g14.6)' ) x, fx, fx2

  end do

  return
end
subroutine test002
!
!*******************************************************************************
!
!! TEST002 tests BESI0.
!
  implicit none
!
  real besi0
  real fx
  real fx2
  integer n
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST002:'
  write ( *, '(a)' ) '  BESI0 evaluates the Bessel I0 function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X	   Exact F	 BESI0(X)'
  write ( *, '(a)' ) ' '
  n = 0

  do

    call besi0_values ( n, x, fx )

    if ( n == 0 ) then
      exit
    end if

    if ( x < 0.0E+00 ) then
      cycle
    end if

    fx2 = besi0 ( x )

    write ( *, '(f8.4,2g14.6)' ) x, fx, fx2

  end do

  return
end
subroutine test003
!
!*******************************************************************************
!
!! TEST003 tests BESJ.
!
  implicit none
!
  real, parameter :: alpha = 0.0E+00
  real fx
  real fx2(1)
  integer n
  integer nz
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST003:'
  write ( *, '(a)' ) '  BESJ evaluates the Bessel function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X       Exact F       BESJ(0)(X)'
  write ( *, '(a)' ) ' '
  n = 0

  do

    call besj0_values ( n, x, fx )

    if ( n == 0 ) then
      exit
    end if

    if ( x < 0.0E+00 ) then
      cycle
    end if

    call besj ( x, alpha, 1, fx2(1), nz )

    write ( *, '(f8.4,2g14.6)' ) x, fx, fx2(1)

  end do

  return
end
subroutine test004
!
!*******************************************************************************
!
!! TEST004 tests BESJ.
!
  implicit none
!
  real, parameter :: alpha = 1.0E+00
  real fx
  real fx2(1)
  integer n
  integer nz
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST004:'
  write ( *, '(a)' ) '  BESJ evaluates the Bessel function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X       Exact F       BESJ(1)(X)'
  write ( *, '(a)' ) ' '
  n = 0

  do

    call besj1_values ( n, x, fx )

    if ( n == 0 ) then
      exit
    end if

    if ( x < 0.0E+00 ) then
      cycle
    end if

    call besj ( x, alpha, 1, fx2(1), nz )

    write ( *, '(f8.4,2g14.6)' ) x, fx, fx2(1)

  end do

  return
end
subroutine test005
!
!*******************************************************************************
!
!! TEST005 tests BESJ.
!
  implicit none
!
  real alpha
  real fx
  real fx2(1)
  integer n
  integer nu
  integer nz
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST005:'
  write ( *, '(a)' ) '  BESJ evaluates the Bessel function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  NU     X       Exact F       BESJ(NU)(X)'
  write ( *, '(a)' ) ' '

  n = 0

  do

    call besjn_values ( n, nu, x, fx )

    if ( n == 0 ) then
      exit
    end if

    if ( x < 0.0E+00 ) then
      cycle
    end if

    alpha = real ( nu )
    call besj ( x, alpha, 1, fx2(1), nz )

    write ( *, '(i4,f8.4,2g14.6)' ) nu, x, fx, fx2(1)

  end do

  return
end
subroutine test006
!
!*******************************************************************************
!
!! TEST006 tests CFFTB.
!! TEST006 tests CFFTI.
!! TEST006 tests CFFTF.
!
!  Find the autocorrelation to El Nino data using FFT methods.
!
  implicit none
!
  integer, parameter :: n = 168
!
  real acov(0:n-1)
  complex cel(0:2*n-1)
  complex corr(0:2*n-1)
  real el(0:2*n-1)
  real el_sum
  logical ex
  integer i
  real wsave(4*(2*n)+15)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST006'
  write ( *, '(a)' ) '  For Fourier transforms of complex data,'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  CFFTI initializes,'
  write ( *, '(a)' ) '  CFFTF forward transforms data,'
  write ( *, '(a)' ) '  CFFTB backward transforms coefficient.'
!
!  Check to see if the required data file exists.
!
  inquire ( file = 'elnino.dat', exist = ex )

  if ( .not. ex ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST006 - Fatal error!'
    write ( *, '(a)' ) '  Cannot find the data file: elnino.dat '
    return
  end if

  open ( unit = 8, file = 'elnino.dat', status = 'old' )
!
!  Read the data, find the mean of the data.
!
  do i = 0, n-1
    read ( 8, * ) el(i)
  end do

  close ( unit = 8 )

  el_sum = sum ( el(0:n-1) )
!
!  Subtract the mean, and append N zeroes.
!
  el(0:n-1) = el(0:n-1) - el_sum / real ( n )

  el(n:2*n-1) = 0.0E+00
!
!  Make a complex copy of EL.
!
  cel(0:2*n-1) = cmplx ( el(0:2*n-1), 0.0E+00 )
!
!    compute fft of data of length 2n.
!    compute square of magnitude of transform components and place
!       in complex array as real parts.
!    compute inverse transform, throwing away second half and
!       imaginary parts (which are zero), and multiply by length of
!       sequence, 2n.
!
  call cffti ( 2*n, wsave )

  call cfftf ( 2*n, cel, wsave )
!
!  CFFTF returns unscaled transforms.
!  The actual transforms are output divided by 2*N.
!
  corr(0:2*n-1) = abs ( cel(0:2*n-1) / real ( 2 * n ) ) **2
!
!  Since we compute transform times its conjugate, we must divide by
!  (2n) for each, i.e., (2n)**2.
!
  call cfftb ( 2*n, corr, wsave )

  acov(0:n-1) = real ( corr(0:n-1) ) * real ( 2 * n )
!
!  Autocovariance is the inverse transform times the sequence length, 2*N.
!
!  Normally, all the scaling would be done once
!  by dividing by 2*N.  We've broken it up for exposition.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Autocorrelation  by the complex FFT method.'
  write ( *, '(a)' ) ' '
  write ( *, '(5e14.6)' ) acov(0:19) / acov(0)

  return
end
subroutine test007
!
!*******************************************************************************
!
!! TEST007 tests CFFTB_2D.
!! TEST007 tests CFFTF_2D.
!
!  Plot the image and transform of an 8 by 8 unit source
!  in a 64 by 64 array.
!
  implicit none
!
  integer, parameter :: n = 64
  integer, parameter :: lda = n
!
  real a(n,n)
  real dat
  real ermax
  real err
  integer i
  complex image(lda,n)
  complex image2(lda,n)
  integer j
  real wsave(4*n+15)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST007'
  write ( *, '(a)' ) '  For two dimensional complex data:'
  write ( *, '(a)' ) '  CFFTF_2D computes the forward FFT transform;'
  write ( *, '(a)' ) '  CFFTB_2D computes the backward FFT transform.'
  write ( *, '(a)' ) ' '
!
!  Initialize WSAVE.
!
  call cffti ( n, wsave )
!
!  Set up the data.  
!
!  IMAGE is the original data.  
!
!  IMAGE2 is IMAGE scaled by (-1)**(I+J), to place the Fourier coefficients 
!  in the correct place for viewing.
!
  do i = 1, n
    do j = 1, n

      if ( i >=(n/2-4) .and.  i<=(n/2+4) .and. &
           j>=(n/2-4) .and. j<=(n/2+4) ) then
        a(i,j) = 1.0E+00
      else
        a(i,j) = 0.0E+00
      end if

      image(i,j) = cmplx ( a(i,j), 0.0E+00 )
      image2(i,j) = image(i,j) * (-1)**(i-1+j-1)

    end do
  end do
!
!  Compute the forward Fourier transform of IMAGE and IMAGE2.
!
  call cfftf_2d ( lda, n, image, wsave )

  call cfftf_2d ( lda, n, image2, wsave )
!
!  Compute the magnitude of the components of transforms.
!  The actual transforms are unscaled and need to be divided by N*N
!  to be correct.
!
  a(1:n,1:n) = abs ( image(1:n,1:n) )
!
!  Compute the inverse Fourier transform of IMAGE.
!
  call cfftb_2d ( lda, n, image, wsave )
!
!  The transforms need to be divided by N*N to be correct.
!
  image(1:n,1:n) = image(1:n,1:n) / real ( n**2 )
!
!  See if the result agrees with the original data.
!
  ermax = 0.0E+00

  do i = 1, n
    do j = 1, n

      if ( i >=(n/2-4) .and. i <=(n/2+4) .and. &
           j >=(n/2-4) .and. j <=(n/2+4)) then
        dat = 1.0E+00
      else
        dat = 0.0E+00
      end if

      err = abs ( dat - abs ( image(i,j) ) )

      ermax = max ( ermax, err )

    end do
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Maximum error in CFFT2D calculation:'
  write ( *, '(a)' ) ' '
  write ( *, '(g14.6)' ) ermax

  return
end
subroutine test008
!
!*******************************************************************************
!
!! TEST008 tests CFFTF.
!! TEST008 tests CFFTI.
!
!  Using complex discrete Fourier transform, find the approximate
!  Fourier coefficients to Runge's function on [-1,1] with N=16 and
!  N=17.
!
  implicit none
!
  complex coeff(0:16)
  real del
  real f
  integer j
  integer n
  real pi
  real, external :: runge
  complex sqtm1
  real wsave(150)
  real x0
  real xj
!
  sqtm1 = cmplx ( 0.0E+00, -1.0E+00 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST008'
  write ( *, '(a)' ) '  CFFTI initializes the complex FFT routines.'
  write ( *, '(a)' ) '  CFFTF does a forward Fourier transform on'
  write ( *, '(a)' ) '    complex data.'
  write ( *, '(a)' ) ' '

  x0 = -1.0E+00

  do n = 16, 17

    call cffti ( n, wsave )
!
!  Function assumed to be periodic on [-1,1], an interval of length 2.
!
    del = 2.0E+00 / real ( n )
    f = 2.0E+00 * pi() / ( real ( n ) * del )
!
!  First sample point at -1, last at 1-del
!
    do j = 0, n-1
      xj = (-1.0E+00) + j * del
      coeff(j) = cmplx ( runge(xj), 0.0E+00 )
    end do

    call cfftf ( n, coeff, wsave )
!
!  Returned coefficients must be divided by N for correct
!  normaliziation.  Also, note repetition after N/2 in original
!  coefficients.  Scaling because X0 not at origin destroys this
!  to some extent.
!
    write ( *, '(a)' ) ' '
    write ( *, '(a,i4)' ) '  Results for n=' , n
    write ( *, '(a)' ) ' '
    write ( *, '(a,2g14.6)' ) ' czero=', coeff(0)/n*2
    write ( *, '(a)' ) ' j          output from cfftf,             scaled coeffs'

    do j = 1, n-1
       write ( *, '(i5,2e15.6,5x,2e15.6)' ) &
         j, coeff(j), exp(-sqtm1*j*f*x0) * coeff(j)/n *2
    end do

  end do

  return
end
subroutine test009
!
!*******************************************************************************
!
!! TEST009 tests CHKDER.
!
  implicit none
!
  integer, parameter :: n = 5
  integer, parameter :: m = n
  integer, parameter :: ldfjac = n
!
  real err(m)
  real fjac(ldfjac,n)
  real fvec(m)
  real fvecp(m)
  integer i
  integer ido
  integer iflag
  integer j
  integer mode
  real x(n)
  real xp(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST009'
  write ( *, '(a)' ) '  CHKDER compares a user supplied jacobian'
  write ( *, '(a)' ) '  and a finite difference approximation to it'
  write ( *, '(a)' ) '  and judges whether the jacobian is correct.'

  do ido = 1, 2

    if ( ido == 1 ) then

       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) '  On the first test, use a correct jacobian.'

    else if ( ido == 2 ) then

       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) '  Repeat the test, but use a "bad" jacobian'
       write ( *, '(a)' ) '  and see if the routine notices!'
       write ( *, '(a)' ) ' '

     end if
!
!  Set the point at which the test is to be made:
!
    x(1:n) = 0.5E+00

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Evaluation point X:'
    write ( *, '(a)' ) ' '
    do i = 1, n
      write ( *, '(g14.6)' ) x(i)
    end do
 
    mode = 1
    call chkder ( m, n, x, fvec, fjac, ldfjac, xp, fvecp, mode, err )

    iflag = 1

    call f33 ( n, x, fvec, fjac, ldfjac, iflag )
    call f33 ( n, xp, fvecp, fjac, ldfjac, iflag )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Sampled function values F(X) and F(XP)'
    write ( *, '(a)' ) ' '
    do i = 1, m
      write ( *, '(i3,2g14.6)' ) i, fvec(i), fvecp(i)
    end do

    iflag = 2
    call f33 ( n, x, fvec, fjac, ldfjac, iflag )
!
!  Here's where we put a mistake into the jacobian, on purpose.
!
    if ( ido == 2 ) then
      fjac(1,1) = 1.01E+00 * fjac(1,1)
      fjac(2,3) = - fjac(2,3)  
    end if

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Computed jacobian'
    write ( *, '(a)' ) ' '
    do i = 1, m
      write ( *, '(5g14.6)' ) fjac(i,1:n)
    end do

    mode = 2
    call chkder ( m, n, x, fvec, fjac, ldfjac, xp, fvecp, mode, err )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  CHKDER error estimates:'
    write ( *, '(a)' ) '     > 0.5, gradient component is probably correct.'
    write ( *, '(a)' ) '     < 0.5, gradient component is probably incorrect.'
    write ( *, '(a)' ) ' '
    do i = 1, m
      write ( *, '(i6,g14.6)' ) i, err(i)
    end do

  end do

  return
end
subroutine f33 ( n, x, fvec, fjac, ldfjac, iflag )
!
!*******************************************************************************
!
!! F33 is a function/jacobian routine.
!
!
!  Modified:
!
!    17 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of variables.
!
!    Input, real X(N), the variable values.
!
!    Output, real FVEC(N), the function values at X, if IFLAG = 1.
!
!    Output, real FJAC(LDFJAC,N), the N by N jacobian at X, if IFLAG = 2.
!
!    Input, integer LDFJAC, the leading dimension of FJAC, which must
!    be at least N.
!
!    Input, integer IFLAG:
!    1, please compute F(I) (X).
!    2, please compute FJAC(I,J) (X).
!
  implicit none
!
  integer ldfjac
  integer n
!
  real fjac(ldfjac,n)
  real fvec(n)
  integer i
  integer iflag
  integer j
  real prod
  real x(n)
!
!  If IFLAG is 1, we are supposed to evaluate F(X).
!
  if ( iflag == 1 ) then

    do i = 1, n-1
      fvec(i) = x(i) - real ( n + 1 ) + sum ( x(1:n) )
    end do

    fvec(n) = product ( x(1:n) ) - 1.0E+00
!
!  If IFLAG is 2, we are supposed to evaluate FJAC(I,J) = d F(I)/d X(J)
!
  else if ( iflag == 2 ) then

    fjac(1:n-1,1:n) = 1.0E+00

    do i = 1, n-1
      fjac(i,i) = 2.0E+00
    end do

    prod = product ( x(1:n) )

    fjac(n,1:n) = prod / x(1:n)

  end if

  return
end
subroutine test010
!
!*******************************************************************************
!
!! TEST010 tests ERF.
!
  implicit none
!
  real erf
  real fx
  real fx2
  integer n
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST010:'
  write ( *, '(a)' ) '  ERF evaluates the Error function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X       Exact F       ERF(X)'
  write ( *, '(a)' ) ' '

  n = 0

  do

    call erf_values ( n, x, fx )

    if ( n == 0 ) then
      exit
    end if

    fx2 = erf ( x )

    write ( *, '(f8.4,2g14.6)' ) x, fx, fx2

  end do

  return
end
subroutine test011
!
!*******************************************************************************
!
!! TEST011 tests ERFC.
!
  implicit none
!
  real erfc
  real fx
  real fx2
  integer n
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST011:'
  write ( *, '(a)' ) '  ERFC evaluates the Complementary Error function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X       Exact F       ERFC(X)'
  write ( *, '(a)' ) ' '

  n = 0

  do

    call erf_values ( n, x, fx )

    if ( n == 0 ) then
      exit
    end if

    fx = 1.0E+00 - fx
    fx2 = erfc ( x )

    write ( *, '(f8.4,2g14.6)' ) x, fx, fx2

  end do

  return
end
subroutine test012
!
!*******************************************************************************
!
!! TEST012 tests EZFFTB.
!! TEST012 tests EZFFTI.
!! TEST012 tests EZFFTB.
!
!  Find the autocorrelation to El Nino data using real FFT methods.
!
  implicit none
!
  integer, parameter :: n = 168
!
  real a(2*n)
  real acovr(0:2*n-1)
  real azero
  real b(2*n)
  real el(0:2*n-1)
  real el_sum
  logical ex
  integer i
  real wsave(4*(2*n)+15)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST012'
  write ( *, '(a)' ) '  For the FFT of a real data sequence:'
  write ( *, '(a)' ) '  EZFFTI initializes,'
  write ( *, '(a)' ) '  EZFFTF does forward transforms,'
  write ( *, '(a)' ) '  EZFFTB does backward transforms.'
!
!  Check to see if the required data file exists.
!
  inquire ( file = 'elnino.dat', exist = ex )

  if ( .not. ex ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST012 - Fatal error!'
    write ( *, '(a)' ) '  Cannot find the data file: elnino.dat '
    return
  end if

  open ( unit = 8, file = 'elnino.dat', status = 'old' )
!
!  Read the data, find the mean of the data.
!
  do i = 0, n-1
    read ( 8, * ) el(i)
  end do

  close ( unit = 8 )

  el_sum = sum ( el(0:n-1) )
!
!  Subtract the mean, and append N zeroes.
!
  el(0:n-1) = el(0:n-1) - el_sum / real ( n )

  el(n:2*n-1) = 0.0E+00
!
!  fft approach (real).
!
!  Compute FFT of data of length 2n.
!  ezfftf produces correctly scaled a's and b's so no extra scaling
!  needed to get transform.
!
  call ezffti ( 2*n, wsave )

  call ezfftf ( 2*n, el, azero, a, b, wsave )
!
!  Compute array of square of each frequency component and place
!  in cosine array (a's) to be back transformed. set b's to 0.
!  There are n a's, and n b's.
!  Note that care must be taken to compute magnitude correctly,
!    0.5*(a(i)**2+b(i)**2) for i < n, twice that for i=n.
!
  azero = azero**2

  a(1:n-1) = ( a(1:n-1)**2 + b(1:n-1)**2 ) / 2.0E+00
  a(n) = a(n)**2 + b(n)**2

  b(1:n) = 0.0E+00
!
!  Compute the back transform, throwing away its second half.
!
  call ezfftb ( 2*n, acovr, azero, a, b, wsave )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Autocorrelation (real fft) output reduced.'
  write ( *, '(a)' ) ' '
  write ( *, '(5e14.6)' ) acovr(0:19) / acovr(0)

  return
end
subroutine test013
!
!*******************************************************************************
!
!! TEST013 tests EZFFTB.
!! TEST013 tests EZFFTF.
!! TEST013 tests EZFFTI.
!
!  Using the real discrete fourier transform, find the approximate
!  fourier coefficients to runge's function on [-1,1] with n=16 and
!  n=17.
!
  implicit none
!
  integer, parameter :: mcoef = 17
!
  real a(mcoef/2)
  real azero
  real b(mcoef/2)
  real c(mcoef/2)
  logical, parameter :: debug = .true.
  real del
  real dfta(mcoef/2)
  real dftb(mcoef/2)
  real error
  real f
  integer j
  integer k
  integer m
  integer n
  real pi
  real r(mcoef)
  real, external :: runge
  real s(mcoef/2)
  real tn
  real wsave(3*mcoef+15)
  real x
  real x0
  real xj
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST013'
  write ( *, '(a)' ) '  The "EZ" FFT package:'
  write ( *, '(a)' ) '  EZFFTI initializes,'
  write ( *, '(a)' ) '  EZFFTF does forward transforms,'
  write ( *, '(a)' ) '  EZFFTB does backward transforms.'

  x0 = -1.0E+00

  do n = mcoef-1, mcoef

    call ezffti ( n, wsave )
!
!  Function assumed to be periodic on [-1,1], of length 2.
!
    del = 2.0E+00 / n
    f = 2.0E+00 * pi() / ( real ( n ) * del )
!
!  The first sample point is at -1, The last at 1-del
!
    do j = 1, n

      xj = (-1.0E+00) + (j-1) * del
      r(j) = runge ( xj )
!
!  Compute sines and cosines to adjust output of ezfftf to give
!  approximate fourier coefficients.
!
      if ( j <= n/2 ) then
        c(j) = cos ( j * f * x0 )
        s(j) = sin ( j * f * x0 )
      end if

    end do

    call ezfftf ( n, r, azero, a, b, wsave )
!
!  As a convenience this loop can go to N/2.  If N is even, the last B is zero.
!
    do j = 1, n/2
      dfta(j) = a(j) * c(j) - b(j) * s(j)
      dftb(j) = a(j) * s(j) + b(j) * c(j)
    end do

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'EZFFTF results'
    write ( *, '(a,i4)' ) '  N = ' , n
    write ( *, '(a,g14.6)' )  ' AZERO = ', azero
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  j             dfta(j)              dftb(j) '
    write ( *, '(a)' ) ' '

    do j = 1, n/2
      write ( *, '(i6,2g14.6)' ) j, dfta(j), dftb(j)
    end do
!
!  Evaluate interpolant at points on [-1,1]
!
    if ( debug ) then

      m = 21

      write ( *, '(a)' ) ' '
      write ( *, '(a,i4)' ) ' Trigonometric polynomial order n= ', n
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  X    Interpolant    Runge     Error'
      write ( *, '(a)' ) ' '

      do k = 1, m

        x = - 1.0E+00 + 2.0E+00 * ( k - 1.0E+00 ) / real ( m - 1 )

        tn = azero
        do j = 1, n/2
          tn = tn + dfta(j) * cos(j*f*x) + dftb(j) * sin(j*f*x)
        end do

        error = tn - runge ( x )
        write ( *, '(4g14.6)' ) x, tn, runge(x), error

      end do

    end if

  end do

  return
end
subroutine test014
!
!*******************************************************************************
!
!! TEST014 tests FMIN.
!
  implicit none
!
  real a
  real b
  real fmin
  real fp05
  real, external :: fx05
  real tol
  real xstar
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST014'
  write ( *, '(a)' ) '  FMIN, function minimizer.'
  write ( *, '(a)' ) '  Find a minimizer of F(X) = X^3 - 2 * X - 5.'
  write ( *, '(a)' ) ' '

  a = 0.1E+00
  b = 0.9E+00
  tol = 1.0E-06

  xstar = fmin ( a, b, fx05, tol )

  write ( *, '(a)' ) '  Results:'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  X* =      ', xstar
  write ( *, '(a,g14.6)' ) '  F(X*) =   ', fx05 ( xstar)
  write ( *, '(a,g14.6)' ) '  F''(X*) = ', fp05 ( xstar )

  return
end
function fx05 ( x )
!
!*******************************************************************************
!
!! FX05 is a function to be minimized.
!
  implicit none
!
  real fx05
  real x
!
  fx05 = x * ( x * x - 2.0E+00 ) - 5.0E+00

  return
end
function fp05 ( x )
!
!*******************************************************************************
!
!! FP05 is the derivative of a function to be minimized.
!
  implicit none
!
  real fp05
  real x
!
  fp05 = 3.0E+00 * x * x - 2.0E+00

  return
end
subroutine test015
!
!*******************************************************************************
!
!! TEST015 tests FZERO.
!
  implicit none
!
  real ae
  real b
  real c
  real, external :: fx06
  integer iflag
  real re
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST015'
  write ( *, '(a)' ) '  FZERO, single nonlinear equation solver.'
  write ( *, '(a)' ) '  F(X) = X^3 - 2 * X - 5'
  write ( *, '(a)' ) ' '

  b = 2.0E+00
  c = 3.0E+00
  ae = 1.0E-06
  re = 1.0E-06

  write ( *, '(a)' ) '  Initial interval: '
  write ( *, '(2g14.6)' ) b, c
  write ( *, '(a,g14.6)' ) '  Absolute error tolerance=', ae
  write ( *, '(a,g14.6)' ) '  Relative error tolerance=', re

  call fzero ( fx06, b, c, c, re, ae, iflag )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  FZERO results'
  write ( *, '(a)' ) ' '

  if ( iflag /= 1 ) then
    write ( *, '(a,i4)' ) '  FZERO returned error code =', iflag
  end if

  write ( *, '(a,g14.6)' ) '  Estimate of zero = ', b
  write ( *, '(a,g14.6)' ) '  Function value=    ', fx06(b)

  return
end
function fx06 ( x )
!
!*******************************************************************************
!
!! FX06 is a function whose zero is desired.
!
  implicit none
!
  real fx06
  real x
!
  fx06 = x * ( x * x - 2.0E+00 ) - 5.0E+00

  return
end
subroutine test016
!
!*******************************************************************************
!
!! TEST016 tests GAMMA.
!
  implicit none
!
  real fx
  real fx2
  real gamma
  integer n
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST016:'
  write ( *, '(a)' ) '  GAMMA evaluates the Gamma function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X       Exact F       GAMMA(X)'
  write ( *, '(a)' ) ' '

  n = 0

  do

    call gamma_values ( n, x, fx )

    if ( n == 0 ) then
      exit
    end if

    fx2 = gamma ( x )

    write ( *, '(f8.4,2g14.6)' ) x, fx, fx2

  end do

  return
end
subroutine test017
!
!*******************************************************************************
!
!! TEST017 tests PCHEV.
!! TEST017 tests PCHEZ.
!! TEST017 tests PCHQA.
!
  implicit none
!
  integer, parameter :: n = 21
  integer, parameter :: nwk = 42
  integer, parameter :: ne = 101
!
  real a
  real b
  real d(n)
  real error
  real errord
  real f(n)
  real fd(ne)
  real fe(ne)
  integer i
  integer ierr
  real pchqa
  real q
  real, external :: runge
  real, external :: rungep
  logical spline
  real wk(nwk)
  real x(n)
  real xe(ne)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST017'
  write ( *, '(a)' ) '  PCHEZ carries out piecewise cubic '
  write ( *, '(a)' ) '    Hermite interpolation.'
  write ( *, '(a)' ) '  PCHEV evaluates the interpolant.'
  write ( *, '(a)' ) '  PCHQA integrates the interpolant.'
  write ( *, '(a)' ) ' '
!
!  Compute Runge's function at N points in [-1,1].
!
  do i = 1, n
    x(i) = -1.0E+00 + real ( i - 1 ) / 10.0E+00
  end do

  f(1:n) = runge ( x(1:n) )

  spline = .false.
!
!  Compute cubic hermite interpolant because spline is .false.
!
  call pchez ( n, x, f, d, spline, wk, nwk, ierr )

  if ( ierr < 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST017 - Error!'
    write ( *, '(a,i4)' ) '  PCHEZ returned error code IERR = ',ierr
    return
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  PCHEZ has set up the interpolant.'
  end if
!
!  Evaluate interpolant and derivative at NE points from -1 to 0.
!
  do i = 1, ne
    xe(i) = -1.0E+00 + real ( i - 1 ) / real ( ne - 1 )
  end do

  call pchev ( n, x, f, d, ne, xe, fe, fd, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST017 - Error!'
    write ( *, '(a,i4)' ) '  PCHEV returned error code IERR = ',ierr
    return
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  PCHEV has evaluated the interpolant.'
  end if
!
!  Examine the error in the function and derivative.
!
  do i = 1, ne
    error = fe(i) - runge ( xe(i) )
    errord = fd(i) - rungep ( xe(i) )
    write ( *, '(2f8.4,2g12.4)' ) xe(i), fe(i), error, errord
  end do
!
!  Compute the integral over the interval [0,1].
!
  a = 0.0E+00
  b = 1.0E+00
  q = pchqa ( n, x, f, d, a, b, ierr )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  PCHQA estimates the integral from A to B.'
  write ( *, '(a,g14.6)' ) '  A = ', a
  write ( *, '(a,g14.6)' ) '  B = ', b
  write ( *, '(a,g14.6)' ) '  Integral estiamte = ', q
  write ( *, '(a,i4)' ) '  Return code IERR = ', ierr

  return
end
subroutine test018
!
!*******************************************************************************
!
!! TEST018 tests Q1DA.
!
  implicit none
!
  real a
  real b
  real e
  real eps
  real, external :: f10
  integer iflag
  integer kf
  real r
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST018'
  write ( *, '(a)' ) '  Q1DA, quadrature routine.'
  write ( *, '(a)' ) ' '

  a = 0.0E+00
  b = 1.0E+00

  eps = 0.001E+00

  call q1da ( f10, a, b, eps, r, e, kf, iflag )

  write ( *, '(a)' ) 'Q1DA results: a, b, eps, r, e, kf, iflag'
  write ( *, '(3f7.4,2e16.8,2i4)' ) a,b,eps,r,e,kf,iflag

  return
end
function f10 ( x )
!
!*******************************************************************************
!
!! F10 is a function to be integrated.
!
  implicit none
!
  real f10
  real x
!
  f10 = sin ( 2.0E+00 * x ) - sqrt ( x )

  return
end
subroutine test019
!
!*******************************************************************************
!
!! TEST019 tests Q1DAX.
!
  implicit none
!
  integer, parameter :: nmax = 50
!
  real a
  real b
  real e
  real eps
  real, external :: f105
  real fmax
  real fmin
  integer i
  integer iflag
  integer int_num
  integer kf
  real r
  logical rst
  real w(nmax,6)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST019'
  write ( *, '(a)' ) '  Q1DAX estimates the integral of a function over a'
  write ( *, '(a)' ) '  a finite interval, allowing more flexibility than Q1DA.'
  write ( *, '(a)' ) ' '

  a = 0.0E+00
  b = 1.0E+00
!
!  Set up an initial partition of 2 intervals, with an internal
!  partition point at 0.3.
!
  w(1,1) = a
  w(2,1) = 0.3E+00
  w(3,1) = b

  int_num = 2

  do i = 1, 2

    if ( i == 1 ) then
      eps = 0.001E+00
      rst = .false.
    else
      eps = 0.0001E+00
      rst = .true.
    end if

    call q1dax ( f105, a, b, eps, r, e, int_num, rst, w, nmax, fmin, fmax, &
      kf, iflag )

    if ( iflag >= 3 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a,i4)' ) '  Q1DAX error flag = ', iflag
      exit
    end if

    write ( *, '(a)' ) ' '
    write ( *, '(a,g14.6)' ) '  Error tolerance      ', eps
    write ( *, '(a,g14.6)' ) '  Integral estimate    ', r
    write ( *, '(a,g14.6)' ) '  Error estimate       ', e

  end do

  return
end
function f105 ( x )
!
!*******************************************************************************
!
!! F105 is a function to be integrated.
!
  implicit none
!
  real f105
  real x
!
  if ( x < 0.3E+00 ) then
    f105 = x**( 0.2E+00 ) * log ( x )
  else
    f105 = sin ( x )
  end if

  return
end

subroutine test020
!
!*******************************************************************************
!
!! TEST020 tests QAGI.
!
!  compute integral of exp(-x)*cos(x*x)**2 on [0,infinity)
!  Correct result is 0.70260...
!
  implicit none
!
  integer, parameter :: limit = 100
  integer, parameter :: lenw = limit * 4
!
  real abserr
  real bound
  real epsabs
  real epsrel
  real, external :: f11
  integer ier
  integer inf
  integer iwork(limit)
  integer last
  integer neval
  real result
  real work(lenw)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST020'
  write ( *, '(a)' ) '  QAGI estimates an integral on a semi-infinite interval.'

  bound = 0.0E+00
  inf = 1
  epsabs = 0.0E+00
  epsrel = 1.0E-05

  call qagi ( f11, bound, inf, epsabs, epsrel, result, abserr, neval, &
    ier, limit, lenw, last, iwork, work )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Estimated integral =   ', result
  write ( *, '(a,g14.6)' ) '  Estimated error =      ', abserr
  write ( *, '(a,i4)' ) '  Function evaluations = ', neval
  write ( *, '(a,i4)' ) '  Return code IER =      ', ier

  return
end
function f11 ( x )
!
!*******************************************************************************
!
!! F11 is a function to be integrated.
!
  implicit none
!
  real f11
  real x
!
  f11 = exp ( -x ) * cos ( x**2 )**2

  return
end
subroutine test021
!
!*******************************************************************************
!
!! TEST021 tests QK15.
!
!  Compute erf(1), i.e. integral of 2/sqrt(pi) * exp(-x*x) from 0
!  to 1.0E+00
!
  implicit none
!
  real a
  real abserr
  real b
  real, external :: f12
  real pi
  real resabs
  real resasc
  real result
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST021'
  write ( *, '(a)' ) '  QK15 estimates an integral using '
  write ( *, '(a)' ) '  Gauss-Kronrod integration.'
  write ( *, '(a)' ) ' '

  a = 0.0E+00
  b = 1.0E+00

  call qk15 ( f12, a, b, result, abserr, resabs, resasc )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  QK15 estimate of ERF(1) '
  write ( *, '(a)' ) '  2.0/sqrt(pi())*result,      abserr'
  write ( *, '(2g14.6)' ) 2.0/sqrt(pi())*result, abserr

  return
end
function f12 ( x )
!
!*******************************************************************************
!
!! F12 is a function to be integrated.
!
  implicit none
!
  real f12
  real x
!
  f12 = exp ( - x**2 )

  return
end
subroutine test022
!
!*******************************************************************************
!
!! TEST022 tests Q1DA.
!
!  double integral by two one dimensional subroutines
!
  implicit none
!
  real a_x
  real b_x
  real eps_x
  real err
  real, external :: g13
  integer iflag
  integer kf
  real result
  real x_fixed
!
  common /comm13/ x_fixed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST022'
  write ( *, '(a)' ) '  Demonstration of two-dimensional quadrature'
  write ( *, '(a)' ) '  using one-dimensional techniques.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Q1DA is called to integrate'
  write ( *, '(a)' ) '  G(X) = Integral ( F(X,Y) dY ).'
 
  a_x = 0.0E+00
  b_x = 1.0E+00
  eps_x = 1.0E-04

  call q1da ( g13, a_x, b_x, eps_x, result, err, kf, iflag )

  write ( *, '(a,2g14.6)' ) eps_x, result

  return
end
function g13 ( x )
!
!*******************************************************************************
!
!! G13 returns G13(X) = Integral ( F(X,Y) dY )
!
!
!  The routine is used as part of an effort to integrate F(X,Y)
!  over a two dimensional region.
!
  implicit none
!
  real a_y
  real b_y
  real eps
  real err
  real, external :: f13
  real g13
  integer iflag
  integer kf
  real result
  real x
  real x_fixed
!
  common /comm13/ x_fixed
!
  x_fixed = x

  a_y = 0.0E+00
  b_y = 2.0E+00
  eps = 1.0E-04

  call q1da ( f13, a_y, b_y, eps, result, err, kf, iflag )
 
  g13 = result

  return
end
function f13 ( y )
!
!*******************************************************************************
!
!! F13 is a function of two variables to be integrated.
!
!
!  Discussion:
!
!    The function is being integrated over a Y range, with a fixed
!    value of X.  The fixed value of X is passed through a common block.
!
  implicit none
!
  real f13
  real x_fixed
  real y
!
  common /comm13/ x_fixed
!
  f13 = exp ( - x_fixed**2 * y**2 )

  return
end
subroutine test023
!
!*******************************************************************************
!
!! TEST023 tests RNOR.
!
  implicit none
!
  integer, parameter :: nbins = 32
!
  real, parameter :: a = -3.0E+00
  real, parameter :: b = 3.0E+00
  integer h(nbins)
  integer i
  integer inbin
  integer iseed
  integer j
  integer, parameter :: nr = 10000
  real r
  real rnor
  real rseed
  real rstart
  real width
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST023'
  write ( *, '(a)' ) '  RNOR, normal random number generator.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Number of normal values to compute is ', nr
  write ( *, '(a,i6)' ) '  Number of bins is ', nbins

  iseed = 305
  rseed = rstart ( iseed )
  width = ( b - a ) / real  ( nbins - 2 )

  h(1:nbins) = 0

  do i = 1, nr
    r = rnor()
    j = inbin ( r, nbins, a, b, width )
    h(j) = h(j) + 1
  end do

  write ( *, '(a)' ) 'Histogram for rnor: number in bin 1,...,32'
  write ( *, '(a)' ) '   (-infinity,-3],(-3,-2.8],...,(2.8,3],(3,infinity)'
  write ( *, '(a)' ) ' (values are slightly computer dependent)'
  write ( *, '(9i8)' ) h(1:nbins)

  return
end
subroutine test024
!
!*******************************************************************************
!
!! TEST024 tests RNOR.
!
  implicit none
!
  integer i
  integer iseed
  real r
  real rnor
  real rseed
  real rstart
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST024'
  write ( *, '(a)' ) '  RNOR generates random normal numbers.'
  write ( *, '(a)' ) ' '
!
!  set initial seed
!
  iseed = 305
  rseed = rstart ( iseed )
!
!  RSTART returns a floating echo of ISEED.
!
  write ( *, '(a,i20)' ) ' ISEED = ', iseed
  write ( *, '(a,g14.6)' ) ' RSEED = ', rseed
  write ( *, '(a)' ) ' '

  do i = 1, 3
    r = rnor()
    write ( *, '(g14.6)' )r
  end do

  return
end
subroutine test025
!
!*******************************************************************************
!
!! TEST025 tests SDRIV1.
!
!
!  An example of the use of the ODE solver SDRIV1.
!
!  Here we solve the simple system
!
!    Y1' = Y2
!    Y2' = -Y1
!
!  with initial conditions
!
!    Y1(0) = 0
!    Y2(0) = 1
!
!  with exact solution
!
!    Y1(T) = SIN(T)
!    Y2(T) = COS(T)
!
  implicit none
!
  integer, parameter :: n = 2

  integer, parameter :: lenw = n * n + 11 * n + 225
!
  real eps
  integer i
  integer mstate
  integer nstep
  real pi
  real t
  real tout
  real work(lenw)
  real y(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST025'
  write ( *, '(a)' ) '  SDRIV1 is a simple interface to the ODE solver.'
  write ( *, '(a)' ) ' '
!
!  Set the error tolerance.
!
  eps = 0.00001E+00
!
!  Set the initial time.
!
  t = 0.0E+00
  tout = t
!
!  Set the initial conditions
!
  y(1) = 0.0E+00
  y(2) = 1.0E+00
!
!  Set the number of steps we will take in the DO loop.
!
  nstep = 12
!
!  Tell SDRIV1 that this is the first call for this problem.
!
  mstate = 1
!
!  Print a header for the results.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Results'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '   t        y(1)      y(2)'
  write ( *, '(a)' ) '            sin(t)    cos(t)'
  write ( *, '(a)' ) '            error     error'
!
!  Call SDRIV1 NSTEP+1 times.
!
  do i = 0, nstep

    tout = real ( 2 * i ) * pi() / real ( nstep ) 

    call sdriv1 ( n, t, y, tout, mstate, eps, work, lenw )

    write ( *, '(a)' ) ' '
    write ( *, '(3f11.5)' ) t, y(1), y(2)
    write ( *, '(11x,2f11.5)' ) sin(t), cos(t)
    write ( *, '(11x,2f11.5)' ) y(1)-sin(t), y(2)-cos(t)
!
!  Cancel the computation if we get any output code but 1 or 2.
!
    if ( mstate /= 1 .and. mstate /= 2 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST16 - Fatal error!'
      write ( *, '(a,i4)' ) 'SDRIV1 returned MSTATE=',mstate
      write ( *, '(a)' ) 'The computation is being cancelled.'
      return
    end if

  end do

  return
end
subroutine f ( n, t, y, ydot )
!
!*******************************************************************************
!
!! F evaluates the right hand sides of the ODE's.
!
  implicit none
!
  integer n
!
  real t
  real y(n)
  real ydot(n)
!
  ydot(1) = y(2)
  ydot(2) = -y(1)

  return
end
subroutine test026
!
!*******************************************************************************
!
!! TEST026 tests SDRIV2.
!
  implicit none
!
  integer, parameter :: n = 2
  integer, parameter :: nroot = 1
!
  integer, parameter :: lw = n*n+10*n+2*nroot+204
  integer, parameter :: liw = 23
!
  real eps
  real ewt
  external fsub
  real, external :: gfun
  real, parameter :: h = 10.0E+00
  integer iw(liw)
  real, parameter :: mass = 0.125E+00
  integer, parameter :: mint = 2
  integer ms
  real t
  real tout
  real w(lw)
  real y(n+1)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST026'
  write ( *, '(a)' ) '  SDRIV2 is an ODE solver.'
  write ( *, '(a)' ) ' '

  eps = 1.0E-05
!
!  Set initial point
!
  t = 0.0E+00
  tout = t
!
!  Set for pure relative error
!
  ewt = 0.0E+00
!
!  Set the initial conditions.
!
  y(1) = h
  y(2) = 0.0E+00
!
!  Set the parameter value.
!
  y(3) = mass
!
!  Set MS to 1, signaling the beginning of the run.
!
  ms = 1

  write ( *, '(a)' ) 'sdriv2 results'
  write ( *, '(a)' ) '   t,         y(1),      y(2),     ms '

  do

    call sdriv2 ( n, t, y, fsub, tout, ms, nroot, eps, ewt, mint, w, &
      lw, iw, liw, gfun )

    tout = tout + 0.1E+00

    if ( ms == 5 )then
      write ( *, '(3f11.5,i4,a,f11.5)' ) t, y(1), y(2), ms, ' <-- y=0 at t= ', t
      exit
    else
      write ( *, '(3f11.5,i4)' ) t,y(1),y(2),ms
!
!  Stop if any output code but 1 or 2.
!
      if ( ms > 2 ) then
        exit
      end if

    end if

  end do

  return
end
subroutine fsub ( n, t, y, ydot )
!
!*******************************************************************************
!
!! FSUB returns the right hand side of an ODE.
!
  implicit none
!
  integer n
!
  real, parameter :: g = 32.0E+00
  real t
  real y(n)
  real ydot(n)
!
  ydot(1) = y(2)
  ydot(2) = - g - y(2) / y(3)

  return
end
function gfun ( n, t, y, iroot )
!
!*******************************************************************************
!
!! GFUN ??
!
  implicit none
!
  integer n
!
  real gfun
  integer iroot
  real t
  real y(n)
!
  gfun = y(1)

  return
end
subroutine test027
!
!*******************************************************************************
!
!! TEST027 tests SINT.
!! TEST027 tests SINTI.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  real c(n)
  integer i
  real wsave((5*n+30)/2)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST027'
  write ( *, '(a)' ) '  For sine analysis of real data,'
  write ( *, '(a)' ) '  SINTI initializes the FFT routines.'
  write ( *, '(a)' ) '  SINT does a forward or backward FFT.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  NOTE THAT SOMETHING IS WRONG,'
  write ( *, '(a)' ) '  EITHER WITH THE PROBLEM, OR THE CODE.'
  write ( *, '(a)' ) '  THE TRANSFORM IS NOT PROPERLY INVERTED.'
!
!  Set the data values.
!
  call rvec_random ( alo, ahi, n, c )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The first 10 data values:'
  write ( *, '(a)' ) ' '
  do i = 1, 10
    write ( *, '(g14.6)' ) c(i)
  end do
!
!  Initialize the WSAVE array.
!
  call sinti ( n, wsave )
!
!  Compute the coefficients.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Compute the sine coefficients from data.'

  call sint ( n, c, wsave )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The first 10 sine coefficients:'
  write ( *, '(a)' ) ' '
  do i = 1, 10
    write ( *, '(g14.6)' ) c(i)
  end do
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Retrieve data from coeficients.'

  call sint ( n, c, wsave )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The first 10 data values:'
  write ( *, '(a)' ) ' '
  do i = 1, 10
    write ( *, '(g14.6)' ) c(i) / real ( 2 * ( n + 1 ) )
  end do

  return
end
subroutine test028
!
!*******************************************************************************
!
!! TEST028 tests SNSQE.
!
  implicit none
!
  integer, parameter :: n = 2
  integer, parameter :: lw = 19
!
  external f18
  real fvec(n)
  integer iflag
  integer info
  integer iopt
  external j18
  integer nprint
  real tol
  real w(lw)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST028'
  write ( *, '(a)' ) '  SNSQE, nonlinear equation system solver.'
!
!  Set the parameters for SNSQE.
!
  tol = 1.0E-05

  x(1:2) = (/ 2.0E+00, 3.0E+00 /)

  call f18 ( n, x, fvec, iflag )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Initial solution estimate X0:'
  write ( *, '(a)' ) ' '
  write ( *, '(4x,2g14.6)' ) x(1:2)
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Function value F(X0):'
  write ( *, '(a)' ) ' '
  write ( *, '(4x,2g14.6)' ) fvec(1:2)

  iopt = 2
  nprint = 0
!
!  Solve the nonlinear equations.
!
  call snsqe ( f18, j18, iopt, n, x, fvec, tol, nprint, info, w, lw )
!
!  Print the results.
!
  if ( info /= 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SNSQE INFO flag = ', info
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  SNSQE solution estimate X:'
  write ( *, '(a)' ) ' '
  write ( *, '(4x,2g14.6)' ) x(1:2)
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Function value F(X):'
  write ( *, '(a)' ) ' '
  write ( *, '(4x,2g14.6)' ) fvec(1:2)

  return
end
subroutine f18 ( n, x, fvec, iflag )
!
!*******************************************************************************
!
!! F18 evaluates a set of nonlinear equations whose zero is sought.
!
  implicit none
!
  integer n
!
  real fvec(n)
  integer iflag
  real x(n)
!
  fvec(1) = x(1) * x(2) - x(2)**3 - 1.0E+00
  fvec(2) = x(1)**2 * x(2) + x(2) - 5.0E+00

  return
end
subroutine j18 ( n, x, fvec, fjac, ldfjac, iflag )
!
!*******************************************************************************
!
!! J18 is a dummy routine.
!
  implicit none
!
  integer ldfjac
  integer n
!
  real fjac(ldfjac,n)
  real fvec(n)
  integer iflag
  real x(n)
!
  return
end
subroutine test029
!
!*******************************************************************************
!
!! TEST029 tests SQRLS.
!
  implicit none
!
  integer, parameter :: mm = 5
  integer, parameter :: nn = 3
!
  real a(mm,nn)
  real b(mm)
  integer i
  integer ind
  integer itask
  integer j
  integer jpvt(nn)
  integer kr
  integer m
  integer n
  real qraux(nn)
  real tol
  real work(nn)
  real x(nn)
!
!  Set up least-squares problem
!  quadratic model, equally-spaced points
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST029'
  write ( *, '(a)' ) '  SQRLS solves linear systems in the least squares sense.'
  write ( *, '(a)' ) ' '
!
  m = 5
  n = 3
  do i = 1, m
    a(i,1) = 1.0E+00
    do j = 2, n
      a(i,j) = a(i,j-1) * i
    end do
  end do

  b(1:5) = (/ 1.0E+00, 2.3E+00, 4.6E+00, 3.1E+00, 1.2E+00 /)

  tol = 1.0E-06

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Coefficient matrix'
  write ( *, '(a)' ) ' '
  do i = 1, m
    write ( *, '(3f12.6)' ) a(i,1:n)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Right-hand side'
  write ( *, '(a)' ) ' '
  write ( *, '(5f12.6)' ) b(1:m)
!
!  Solve least-squares problem
!
  itask = 1
  call sqrls ( a, mm, m, n, tol, kr, b, x, b, jpvt, qraux, work, itask, ind )
!
!  Print results
!
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  Error code =', ind
  write ( *, '(a,i4)' ) '  Estimated matrix rank =', kr
  write ( *, '(a)' ) '  Parameters'
  write ( *, '(3f12.6)' ) x(1:n)
  write ( *, '(a)' ) '  Residuals'
  write ( *, '(5f12.6)' ) b(1:m)

  return
end
subroutine test030
!
!*******************************************************************************
!
!! TEST030 tests SSVDC.
!
  implicit none
!
  integer, parameter :: ldx = 8
  integer, parameter :: n = 8
  integer, parameter :: p = 3
  integer, parameter :: ldu = n
  integer, parameter :: ldv = p
  integer, parameter :: job = 11
!
  real c(p)
  real e(p)
  integer i
  integer info
  integer j
  integer jb
  real, dimension ( n ) :: pop = (/ 75.994575E+00, 91.972266E+00, &
    105.710620E+00, 122.775046E+00, 131.669275E+00, 150.697361E+00, & 
    179.323175E+00, 203.235298E+00 /)
  real pop80
  real r
  real relerr
  real ri
  real rsq
  real s(p)
  real sum2
  real tol
  real u(ldu,ldu)
  real v(ldv,p)
  real w(n)
  real x(ldx,p)
  real y(n)
  real year
!
!  C contains coefficients of polynomial
!
!    c(1)*1+c(2)*t+c(3)*t*t
!
!  where t=year (1900 etc.).
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST030'
  write ( *, '(a)' ) '  SSVDC computes the singular value decomposition.'

  do i = 1, 8
    y(i) = 1900.0E+00 + real ( ( i - 1 ) * 10 )
  end do

  x(1:8,1) = 1.0E+00
  x(1:8,2) = y(1:8)
  x(1:8,3) = y(1:8)**2

  call ssvdc ( x, ldx, n, p, s, e, u, ldu, v, ldv, w, job, info )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Computed singular values: '
  write ( *, '(a)' ) ' '
  write ( *, '(5g12.4)' ) s(1:p)

  c(1:p) = 0.0E+00
!
!  RELERR reflects number of accurate digits in data
!  e.g. 6 digits ==> relerr=1.0E-06, and so on.
!  Making relerr larger increases residuals.
!
  relerr = 1.0E-06
  tol = relerr * s(1)
!
!  Multiply U' * pop, and solve for coefficients c(i)
!
  do j = 1, p

    if ( s(j) > tol ) then

      sum2 = dot_product ( pop(1:n), u(1:n,j) ) / s(j)

      c(1:p) = c(1:p) + sum2 * v(1:p,j)

    end if

  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Computed polynomial coefficients:'
  write ( *, '(a)' ) ' '
  write ( *, '(5g12.4)' ) c(1:p)
!
!  Evaluate the model using Horner's rule, and residuals at
!  year =1900,...,1980
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '         Model       True'
  write ( *, '(a)' ) 'Year  Population  Population  Error'
  write ( *, '(a)') ' '

  r = 0.0E+00

  do i = 1, 9

    year = 1900.0E+00 + real ( i - 1 ) * 10.0E+00

    pop80 = 0.0E+00
    do j = p, 1, -1
      pop80 = year * pop80 + c(j)
    end do

    if ( i < 9 ) then
      r = r + ( pop(i) - pop80 )**2
      write ( *, '(i4,3f10.2)' ) int(year), pop80, pop(i), pop(i) - pop80
    else
      write ( *, '(i4,f10.3)' ) int(year), pop80
    end if

  end do

  r = sqrt ( r )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  RMS error is ', r

  return
end
subroutine f21 ( n, x, f )
!
!*******************************************************************************
!
!! F21 is a function to be minimized.
!
!
!  Parameters:
!
!    Input, integer N, the size of the X vector.
!
!    Input, real X(N), the value of the X vector.
!
!    Output, real F, the value of the function F(X).
!
  implicit none
!
  integer n
!
  real f
  integer i
  real t1
  real t2
  real x(n)
!
  t1 = 0.0E+00
  do i = 2, n
    t1 = t1 + ( x(i) - x(i-1)**2 )**2
  end do

  t2 = sum ( ( 1.0E+00 - x(1:n-1) )**2 )

  f = 1.0E+00 + 100.0E+00 * t1 + t2

  return
end
subroutine test031
!
!*******************************************************************************
!
!! TEST031 tests SGEFS.
!
  implicit none
!
  integer, parameter :: lda = 10
!
  real a(lda,lda)
  real b(lda)
  integer i
  integer ind
  integer itask
  integer iwork(lda)
  integer j
  integer n
  real rcond
  real work(lda)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST031'
  write ( *, '(a)' ) '  SGEFS solves a system of linear equations.'
  write ( *, '(a)' ) ' '
!
!  Set the number of equations.
!
  n = 3

  itask = 1
!
!  Set the coefficient matrix A.
!
  a(1,1) = 10.0E+00
  a(2,1) = -3.0E+00
  a(3,1) =  5.0E+00
  a(1,2) = -7.0E+00
  a(2,2) =  2.0E+00
  a(3,2) = -1.0E+00
  a(1,3) =  0.0E+00
  a(2,3) =  6.0E+00
  a(3,3) =  5.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Coefficient matrix A:'
  write ( *, '(a)' ) ' '

  do i = 1, n
    write ( *, '(3f12.6)' ) a(i,1:n)
  end do
!
!  Set the right hand side vector B.
!
  b(1:3) = (/ 7.0E+00, 4.0E+00, 6.0E+00 /)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Right-hand side B:'
  write ( *, '(a)' ) ' '
  write ( *, '(3f12.6)') b(1:n)
!
!  Solve the linear system A*x=b.
!
  call sgefs ( a, lda, n, b, itask, ind, work, iwork, rcond )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'SGEFS results:'
  write ( *, '(a)' ) ' '

  if ( ind == -10 ) then
    write ( *, '(a,i4)' ) '  Error code =',ind
  else if ( ind < 0 ) then
    write ( *, '(a,i4)' ) '  Error code =',ind
    return
  else
    write ( *, '(a,i4)' ) 'Estimated number of accurate digits =', ind
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Solution:'
  write ( *, '(a)' ) ' '
  write ( *, '(3f12.6)' ) b(1:n)

  return
end
subroutine test032
!
!*******************************************************************************
!
!! TEST032 tests UNCMIN.
!
  implicit none
!
  integer, parameter :: n = 2
  integer, parameter :: lwork = n*(n+10)
!
  real f
  external f8
  integer i
  integer ierror
  real work(lwork)
  real x(n)
  real x0(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST032'
  write ( *, '(a)' ) '  UNCMIN, unconstrained minimization code.'
  write ( *, '(a)' ) ' '
!
!  Specify an initial estimate of the solution.
!
  x0(1) = 1.0E+00
  x0(2) = 1.0E+00
!
!  Minimize function
!
  call uncmin ( n, x0, f8, x, f, ierror, work, lwork )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  UNCMIN return code =', ierror
  write ( *, '(a,g14.6)' ) '  f(x*) =', f
  write ( *, '(a)' ) '  x* ='
  write ( *, '(4g14.6)' ) x(1:n)

  return
end
subroutine f8 ( n, x, f )
!
!*******************************************************************************
!
!! F8 is a function to be minimized.
!
  implicit none
!
  integer, parameter :: m = 4
  integer n
!
  real, parameter, dimension ( m ) :: b = &
    (/ 20.0E+00, 9.0E+00, 3.0E+00, 1.0E+00 /)
  real f
  integer j
  real, parameter, dimension ( m ) :: t = &
    (/ 0.0E+00, 1.0E+00, 2.0E+00, 3.0E+00 /)
  real x(n)
!
  f = 0.0E+00
  do j = 1, m
    f = f + ( b(j) - x(1) * exp ( x(2) * t(j) ) )**2
  end do

  return
end
subroutine test033
!
!*******************************************************************************
!
!! TEST033 tests UNCMIN.
!
  implicit none
!
  integer, parameter :: n = 10
  integer, parameter :: lwork = n*(n+10)
!
  real f
  integer i
  integer ierror
  real work(lwork)
  real x(n)
  real x0(n)
!
  external f21
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST033'
  write ( *, '(a)' ) '  UNCMIN carries out unconstrained minimization'
  write ( *, '(a)' ) '  of a scalar function of several variables.'

  do i = 1, n
    x0(i) = real ( i ) / real ( n + 1 )
  end do
!
!  minimize function
!
  call uncmin ( n, x0, f21, x, f, ierror, work, lwork )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  Return code =', ierror
  write ( *, '(a,g14.6)' ) '  f(x*) =', f
  write ( *, '(a)' ) '  x* ='
  write ( *, '(5f12.6)' ) x(1:n)

  return
end
subroutine test034
!
!*******************************************************************************
!
!! TEST034 tests UNI.
!
  implicit none
!
  integer i
  integer iseed
  real u
  real uni
  real ustart
  real useed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST034'
  write ( *, '(a)' ) '  UNI is a uniform random number generator.'
  write ( *, '(a)' ) ' '
!
!  Set the initial seed.
!
  iseed = 305
  useed = ustart ( iseed )
!
!  USTART returns floating echo of iseed.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a,i20)' ) '  The seed value ISEED is ', iseed
  write ( *, '(a,g14.6)' ) '  The starting value is ', useed

  do i = 1, 1000
    u = uni()
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  The 1000-th random number generated is ', u

  return
end
subroutine test035
!
!*******************************************************************************
!
!! TEST035 computes an autocorrelation using a direct method.
!
  implicit none
!
  integer, parameter :: n = 168
!
  real acov(0:n-1)
  real el(0:2*n-1)
  real el_sum
  logical ex
  integer i
  integer j
  integer m
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST035'
  write ( *, '(a)' ) '  Compute the autocorrelation of El Nino data'
  write ( *, '(a)' ) '  using a direct method.'
!
!  Check to see if the required data file exists.
!
  inquire ( file = 'elnino.dat', exist = ex )

  if ( .not. ex ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST004 - Fatal error!'
    write ( *, '(a)' ) '  Cannot find the data file: elnino.dat '
    return
  end if

  open ( unit = 8, file = 'elnino.dat', status = 'old' )
!
!  Read the data, find the mean of the data.
!
  do i = 0, n-1
    read ( 8, * ) el(i)
  end do

  close ( unit = 8 )

  el_sum = sum ( el(0:n-1) )
!
!  Subtract the mean, and append N zeroes.
!
  el(0:n-1) = el(0:n-1) - el_sum / real ( n )

  el(n:2*n-1) = 0.0E+00
!
!  Direct calculation.
!  Only sum as far as there is data.
!  Simple, but slow.
!
  do j = 0, n-1
    acov(j) = 0.0E+00
    do m = 0, n-1-j
      acov(j) = acov(j) + el(m) * el(m+j)
    end do
  end do
!
!  Write the scaled autocorrelation.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Autocorrelation by the direct method.'
  write ( *, '(a)' ) ' '
  write ( *, '(5e14.6)' ) acov(0:19) / acov(0)

  return
end
subroutine test036
!
!*******************************************************************************
!
!! TEST036 is the reactor shielding problem.
!
  implicit none
!
  logical absorb
  real azm
  real d
  real dist2c
  real e
  real ea
  real er
  real et
  integer i
  integer iexit
  real mu
  integer na
  integer npart
  integer nr
  integer nt
  integer ntot
  real sa
  real sr
  real st
  real thick
  real x
  real y
  real z
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST036'
  write ( *, '(a)' ) '  The reactor shielding problem'
  write ( *, '(a)' ) '  Monte Carlo simulation using the UNI code'
  write ( *, '(a)' ) '  for uniform random numbers.'
  write ( *, '(a)' ) ' '
!
!  Set total number of particles.
!
  ntot = 10
!
!  Set the slab thickness.
!
  thick = 2.0E+00
!
!  Initialize the particle counts and energy tallies.
!
  call init ( na, nt, nr, npart, ea, et, er, sa, sr, st )
!
!  Main loop over ntot particles
!
    1 continue

  npart = npart + 1
!
!  Finished, compute and print averages, standard deviations
!
  if ( npart == ntot ) then
    call output ( na, ea, sa, nr, er, sr, nt, et, st, ntot )
    return
  end if
!
!  Source generates a new particle with position, direction, energy
!
  call source ( e, mu, azm, x, y, z )

    2 continue
!
!  Compute-distance-to-collision
!
  d = dist2c ( e )
!
!  Update position to decide if particle exits or collides
!
  call update ( mu, azm, d, x, y, z )

  i = iexit ( thick, x, y, z )
!
!  Returns -1, 0 , 1 for 'out on left', 'in', 'out on right'
!
  if ( i < 0 ) then
!
!  Exit on left (reflected), tally and goto source
!
    nr = nr + 1
    er = er + e
    sr = sr + e * e
    go to 1

  else if ( i > 0 ) then
!
!  Exit on right (transmitted), tally and goto source
!
    nt = nt + 1
    et = et + e
    st = st + e * e
    go to 1

  else if ( i == 0 .and. absorb() ) then
!
!  Collision and absorbed, tally and goto source
!
    na = na + 1
    ea = ea + e
    sa = sa + e * e
    go to 1

  else
!
!  Collision and scattered, find scattering angle and energy
!
    call scatt ( mu, azm, e )
!
!  Go to compute-distance-to-collision
!
    go to 2

  end if

  return
end
function absorb ( )
!
!*******************************************************************************
!
!! ABSORB returns TRUE if the particle is absorbed.
!
!
!  Discussion:
!
!    Absorption occurs with probability PA.
!
!  Parameters:
!
!    Output, logical ABSORB, is TRUE if the particle is absorbed.
!
  implicit none
!
  logical absorb
  real, parameter :: pa = 0.1E+00
  real uni
!
  if ( uni() <= pa ) then
    absorb = .true.
  else
    absorb = .false.
  end if

  return
end
function cross ( e )
!
!*******************************************************************************
!
!! CROSS returns cross section (fictional) for energy in range [emin,emax]
!
  implicit none
!
  real cross
  real e
  real s
  real y
!
  s = abs ( sin ( 100.0E+00 * ( exp ( e ) - 1.0E+00 ) ) &
    + sin ( 18.81E+00 * ( exp ( e ) - 1.0E+00 ) ) )

  y = max ( 0.02E+00, s )
  cross = 10.0E+00 * exp ( -0.1E+00 / y )

  return
end
function dist2c ( e )
!
!*******************************************************************************
!
!! DIST2C returns the distance to collision.
!
!
!  Discussion:
!
!    An exponential distribution is used with parameter `cross section'.
!
  implicit none
!
  real cross
  real dist2c
  real e
  real uni
!
  dist2c = - log ( uni() ) / cross ( e )

  return
end
function energy ( )
!
!*******************************************************************************
!
!! ENERGY returns energy, with distribution const/sqrt(energy) over [emin,emax]
!  use inverse function approach to compute this
!
!  energy min, max in mev
!
  implicit none
!
  real, parameter :: emin = 1.0E-03
  real, parameter :: emax = 2.5E+00
!
  real c
  real energy
  real uni
!
  c = 1.0E+00 / ( 2.0E+00 * ( sqrt ( emax ) - sqrt ( emin ) ) )
  energy = ( uni() / ( 2.0E+00 * c ) + sqrt ( emin ) )**2

  return
end
function iexit ( thick, x, y, z )
!
!*******************************************************************************
!
!! IEXIT determines if the particle where the particle is.
!
!
!  Parameters:
!
!    Input, real THICK, the thickness of the wall.
!
!    Input, real X, Y, Z, the coordinates of the particle.
!
!    Output, integer IEXIT:
!    -1, the particle is to the left of the wall;
!     0, the particle is in the wall;
!    +1, the particle is to the right of the wall.
!
  implicit none
!
  integer iexit
  real thick
  real x
  real y
  real z
!
  if ( x < 0.0E+00 )then
    iexit = -1
  else if ( x <= thick ) then
    iexit = 0
  else
    iexit = 1
  end if

  return
end
subroutine init ( na, nt, nr, npart, ea, et, er, sa, sr, st )
!
!*******************************************************************************
!
!! INIT initializes data for the reactor shielding problem.
!
  implicit none
!
  real ea
  real er
  real et
  integer na
  integer npart
  integer nr
  integer nt
  real sa
  real sr
  real st
!
  na = 0
  nt = 0
  nr = 0
  npart = 0
  ea = 0.0E+00
  et = 0.0E+00
  er = 0.0E+00
  sa = 0.0E+00
  sr = 0.0E+00
  st = 0.0E+00

  return
end
subroutine output ( na, ea, sa, nr, er, sr, nt, et, st, ntot )
!
!*******************************************************************************
!
!! OUTPUT prints the results of the reactor shielding problem.
!
  implicit none
!
  real ea
  real er
  real et
  integer na
  integer nr
  integer nt
  integer ntot
  real sa
  real sr
  real st
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Reactor Shielding Problem Tallies:'

  if ( na > 0 ) then
    ea = ea / na
    sa = sqrt ( sa / na - ea**2 )
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '% absorbed, energy, std dev: '
  write ( *, '(3g14.6)' ) real(na) / ntot*100, ea, sa

  if ( nr > 0 ) then
    er = er / nr
    sr = sqrt ( sr / nr - er * er )
  end if

  write ( *, '(a)' ) '% reflected, energy, sd: '
  write ( *, '(3g14.6)' )  real(nr)/ntot*100, er, sr

  if ( nt > 0 ) then
    et = et / nt
    st = sqrt ( st / nt - et * et )
  end if

  write ( *, '(a)' ) '% transmitted, energy, sd: '
  write ( *, '(3g14.6)' )  real(nt)/ntot*100, et, st

  write ( *, '(a)' ) ' '

  return
end
subroutine scatt ( mu, azm, e )
!
!*******************************************************************************
!
!! SCATT returns scattering angle and energy.
!
  implicit none
!
  real azm
  real e
  real mu
  real pi
  real uni
!
!  Isotropic scattering, i.e., uniform on sphere
!
  mu = -1.0E+00 + 2.0E+00 * uni()
  azm = uni() * 2.0E+00 * pi()
!
!  Find scattering energy, uniform in [.3*e,e]
!
  e = ( uni() * 0.7E+00 + 0.3E+00 ) * e

  return
end
subroutine source ( e, mu, azm, x, y, z )
!
!*******************************************************************************
!
!! SOURCE ??
!
!  emitter returns mu uniform on [0,1] as the cosine of an angle,
!  and azm uniform on [0,2*pi] as asimuthal angle
!
  implicit none
!
  real azm
  real e
  real energy
  real mu
  real pi
  real uni
  real x
  real y
  real z
!
  mu = uni()
  azm = uni() * 2.0E+00 * pi()
  x = 0.0E+00
  y = 0.0E+00
  z = 0.0E+00
!
!  Returns energy from source energy distribution
!
  e = energy()

  return
end
subroutine update ( mu, azm, d, x, y, z )
!
!*******************************************************************************
!
!! UPDATE ??
!
  real azm
  real d
  real mu
  real r
  real x
  real y
  real z
!
  x = x + d * mu
  r = d * sqrt ( 1.0E+00 - mu**2 )
  y = y + r * cos ( azm )
  z = z + r * sin ( azm )

  return
end