subroutine ptrap ( a, b, n, fn, gn, vint )

c*********************************************************************72
c
cc ptrap() carries out product-type trapezoidal quadrature.
c
c  this subroutine uses the product type trapezoidal rule
c  compounded n times to approximate the integral from a to b
c  of the function fn(x) * gn(x).  fn and gn are function
c  subprograms which must be supplied by the user.  the
c  result is stored in vint.
c
      double precision a, ag, am(2,2), b, f(2), fn, g(2),
     &  gn, h, vint, x, dble
      data am(1,1), am(2,2) /2 * 2.d0 /, am(1,2), am(2,1)
     &  / 2 * 1.d0/
      h = ( b - a ) / dble ( float ( n ) )
      vint = 0.d0
      x = a
      f(2) = fn ( a )
      g(2) = gn ( a )

      do i = 1, n
        f(1) = f(2)
        g(1) = g(2)
        x = x + h
        f(2) = fn ( x )
        g(2) = gn ( x )
        do j = 1, 2
          ag = 0.d0
          do k = 1, 2
            ag = ag + am(j,k) * g(k)
          end do
          vint = vint + f(j) * ag
        end do
      end do

      vint = h * vint / 6.d0

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
c
c  Discussion:
c
c    This FORTRAN77 version is made available for cases where the
c    FORTRAN90 version cannot be used.
c
c  Modified:
c
c    16 September 2005
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    None
c
      implicit none

      character ( len = 8 ) date
      character ( len = 10 ) time

      call date_and_time ( date, time )

      write ( *, '(a8,2x,a10)' ) date, time

      return
      end