      subroutine randg ( l, x, r, y )

c*********************************************************************72
c
c  subroutine randg
c
c  purpose:
c
c    compute random numbers from any general distribution.
c
c  usage:
c
c    call randg ( l, x, r, y )
c
c  description of parameters:
c
c    input
c
c      l  a nonzero odd random integer.
c
c      x  vector of length 257 containing ordinate points
c         separated by equal probability on desired distribution.
c         (can be calculated in randgi.)
c
c      r  vector of length 256 containing ratios of derivative*dx
c         to area/dx for each ordinate point in x.
c         (can be calculated in randgi.)
c
c    output
c
c      l  has been updated so that it may be used as a seed for
c      a subsequent call.
c
c      y  random number.
c
c  remarks:
c
c    quadratic approximation of cdf (cumulative density function)
c    which implies linear approximation of pdf (probability density
c    function).
c
c  subroutines and function subprograms required:
c
c    none directly.  randgi may be used for initialization.
c
c  method:
c
c    table lookup plus uniform and triangular distribution
c    variables are used.
c
      implicit none

      real ak2
      integer k1
      integer l
      integer l1
      integer l2
      real r(256)
      real x(257)
      real y
c
c  generate two uniform random numbers on interval ( 1, 2**31 ).
c  any good generator may be substituted.
c
      l1 = iabs ( 65539 * l )
      l = iabs ( 65539 * l1 )
      l2 = l
c
c  calculate two uniform random numbers
c  k1 integer on interval ( 1, 256 ),
c  ak2 real on interval ( 0.0, 1.0 ).
c
      k1 = l1 / 8388608 + 1
      ak2 = float ( mod ( l1, 8388608 ) ) * 1.192093e-07
      if ( ak2 - abs ( r(k1) ) ) 8, 8, 30
8     if ( r(k1) ) 20, 10, 10
c
c  calculate triangular random skewed left.
c
10    y = x(k1) + ( x(k1+1) - x(k1) ) * amax0 ( l1, l2 ) * 4.656613e-10
      return
c
c  calculate triangular random skewed right.
c
20    y = x(k1) + ( x(k1+1) - x(k1) ) * amin0 ( l1, l2 ) * 4.656613e-10
      return
c
c  calculate uniform random.
c
30    y = x(k1) + ( x(k1+1) - x(k1) ) * float ( l2 ) * 4.656613e-10

      return
      end
      subroutine randgi ( n, x, y, p, q, r, k, ier )

c*********************************************************************72
c
c  subroutine randgi
c
c  purpose:
c
c    compute initializing vectors for randg.
c
c  usage:
c
c    call randgi ( n, x, y, p, q, r, k, ier )
c
c  description of parameters:
c
c    input
c
c      n   number of (x,y) points of approximation to pdf
c          (probability density function) or cdf
c          (cumulative density function), n.ge.3.
c
c      x   vector of length n containing ordinate of pdf or cdf.
c
c      y   vector of length n containing abscissa of pdf or cdf.
c
c    output
c
c      p   work vector of length n.
c
c      q   vector or length 257 containing ordinate points
c          separated by equal probability on desired distribution.
c      r   vector of length 256 containing ratios of derivatives
c          to area/dx for each ordinate point in q.
c
c      k   k should be set to 1 if y represents a cdf.  otherwise,
c          y will be interpreted as a pdf.
c
c      ier error indicator
c          1, error in scaling.  i.e., total area of pdf not equal to
c             1.  assuming estimation errors, a fudge factor is used
c             to scale a result.
c          2, density not positive.  i.e., some y(i) .lt. 0.  abort.
c          3, not in sort.  i.e., some x(i) .le. x(i-1).  abort.
c          4, search error.  should never occur because of fudge
c             factor used.  this means some p(i) not large enough
c             for search of proper q.  investigation is needed.
c          5, n < 3.  abort.
c
c  remarks:
c
c    none
c
c  subroutines and function subprograms needed:
c
c    none
c
c  method:
c
c    linear approximation of pdf to find cdf (cumulative density
c    function) and icdf (inverse cdf).  quadratic interpolation on
c    icdf to find q and r.
c
      implicit none

      real f
      integer i
      integer ier
      integer j
      integer j1
      integer k
      integer n
      real p(300)
      real q(257)
      real r(256)
      real t1
      real t2
      real v
      real x(300)
      real xt1
      real xt2
      real xt3
      real xv1
      real xv2
      real xv3
      real y(300)
c
c  calculate cumulative probabilities.
c
      ier = 0

      if ( n < 3 ) then
        ier = 5
        return
      end if

      if ( y(1) ) 5, 10, 10
c
c  error 2.
c
5     ier = 2
      return
10    p(1) = 0.0e+00

      do 15 i = 2, n

        if ( y(i) ) 5, 11, 11
11      if ( x(i) - x(i-1) ) 6, 6, 12
c
c  error 3.
c
6       ier = 3
        return
12      if ( k .eq. 1 ) go to 13
c
c  y is a pdf.
c
        p(i) = ( y(i) + y(i-1) ) * ( x(i) - x(i-1) ) * 0.5e+00 + p(i-1)
        go to 15
c
c  y is a cdf.
c
13      p(i) = y(i)

15    continue

      if ( p(n) - 0.996094e+00 ) 7, 7, 16
16    if ( p(n) - 1.003906e+00 ) 3, 7, 7
c
c  error 1.
c
7     ier = 1
3     f = 1.0e+00 / p(n)
      do i = 2, n
        p(i) = p(i) * f
      end do
c
c  calculate x points for equal-distant cumulative probabilities.
c
      v = 0.0e+00
      q(1) = x(1)
      t1 = y(1)
      j1 = 2

      do 150 i = 2, 257
        if ( i - 257 ) 102, 103, 103
102     v = v + 3.90625e-03
c
c  locate best point for interpolation.
c
        do j = j1, n
          if ( p(j) - v ) 101, 104, 105
101       continue
        end do
c
c  error 4.
c
        ier = 4
103     j = n
104     q(i) = x(j)
        t2 = y(j)
        go to 125

105     if ( j .le. 3 ) go to 113
        j1 = j - 2
        go to 120
113     j1 = 1
c
c  quadratic interpolation of p inverse for q.
c
120     xt2 = p(j1+2) - p(j1)
        xt3 = p(j1+2) - p(j1+1)
        xt1 = p(j1+1) - p(j1)
        xv1 = v - p(j1)
        xv2 = v - p(j1+1)
        xv3 = v - p(j1+2)
        q(i) = ( xv3 * xv2 * x(j1)   ) / ( xt1 * xt2 )
     &       - ( xv3 * xv1 * x(j1+1) ) / ( xt1 * xt3 )
     &       + ( xv2 * xv1 * x(j1+2) ) / ( xt2 * xt3 )
c
c  linear interpolation of y for t2 and r.
c
        t2 = ( y(j) - y(j-1) ) * ( q(i) - x(j-1) )
     &    / ( x(j) - x(j-1) ) + y(j-1)
125     r(i-1) = ( t2 - t1 ) / ( t2 + t1 )
        t1 = t2
        j1 = j
150   continue

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 January 2007
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    None
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ',
     &  'May      ', 'June     ', 'July     ', 'August   ',
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

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

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

      return
      end
