function bivprb ( h, k, r ) 

!*****************************************************************************80
!
!! BIVPRB computes a bivariate normal CDF for correlated X and Y.
!
!  Discussion:
!
!    This routine computes P( H < X, K < Y ) for X and Y standard normal 
!    variables with correlation R. 
!
!  Modified:
!
!    02 February 2008
!
!  Author:
!
!    Original version by Mike Patefield, David Tandy.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Mike Patefield, David Tandy,
!    Fast and Accurate Calculation of Owen's T Function,
!    Journal of Statistical Software,
!    Volume 5, Number 5, 2000, pages 1-25.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) H, K, the lower limits for X and Y.
!    0 <= H, 0 <= K.
!
!    Input, real ( kind = 8 ) R, the correlation between X and Y.
!
!    Output, real ( kind = 8 ) BIVPRB, the requested probability.
!
  implicit none

  real    ( kind = 8 ) bivprb 
  real    ( kind = 8 ) h
  integer ( kind = 4 ) ifail
  real    ( kind = 8 ) k
  real    ( kind = 8 ) q
  real    ( kind = 8 ) r
  real    ( kind = 8 ) ri
  real    ( kind = 8 ) rr
  real    ( kind = 8 ), parameter :: rroot2 = 0.70710678118654752440D+00
  real    ( kind = 8 ), parameter :: rtwopi = 0.15915494309189533577D+00
  real    ( kind = 8 ) znorm2

  if ( r == 0.0D+00 ) then

    bivprb = znorm2 ( h ) * znorm2 ( k )

  else 

    rr = 1.0D+00 - r * r 

    if ( 0.0D+00 < rr ) then

      ri = 1.0D+00 / sqrt ( rr )

      if ( 0.0D+00 < h .and. 0.0D+00 < k ) then 

        bivprb = q ( h, ( k / h - r ) * ri ) &
               + q ( k, ( h / k - r ) * ri ) 

      else if ( 0.0D+00 < h ) then 
        bivprb = q ( h, - r * ri ) 
      else if ( 0.0D+00 < k ) then 
        bivprb = q ( k, - r * ri ) 
      else 
        bivprb = 0.25D+00 + rtwopi * asin ( r )
      end if 

    else if ( 1.0D+00 <= r ) then 

      if ( k <= h ) then 
         bivprb = znorm2 ( h ) 
      else 
        bivprb = znorm2 ( k )
      end if 

    else 

      bivprb = 0.0D+00 

    end if 

  end if 

  return 
end 
subroutine normp ( z, p, q, pdf )

!*****************************************************************************80
!
!! NORMP computes the cumulative density of the standard normal distribution.
!
!  Discussion:
!
!    This is algorithm 5666 from Hart, et al.
!
!  Modified:
!
!    02 February 2008
!
!  Author:
!
!    Alan Miller
!    FORTRAN90 version by John Burkardt
!
!  Reference:
!
!    John Hart, Ward Cheney, Charles Lawson, Hans Maehly,
!    Charles Mesztenyi, John Rice, Henry Thacher,
!    Christoph Witzgall,
!    Computer Approximations,
!    Wiley, 1968,
!    LC: QA297.C64.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) Z, divides the real line into two semi-infinite
!    intervals, over each of which the standard normal distribution
!    is to be integrated.
!
!    Output, real ( kind = 8 ) P, Q, the integrals of the standard normal
!    distribution over the intervals ( - Infinity, Z] and
!    [Z, + Infinity ), respectively.
!
!    Output, real ( kind = 8 ) PDF, the value of the standard normal 
!    distribution at Z.
!
  implicit none

  real    ( kind = 8 ), parameter :: cutoff = 7.071D+00
  real    ( kind = 8 ) expntl
  real    ( kind = 8 ) p
  real    ( kind = 8 ), parameter :: p0 = 220.2068679123761D+00
  real    ( kind = 8 ), parameter :: p1 = 221.2135961699311D+00
  real    ( kind = 8 ), parameter :: p2 = 112.0792914978709D+00
  real    ( kind = 8 ), parameter :: p3 = 33.91286607838300D+00
  real    ( kind = 8 ), parameter :: p4 = 6.373962203531650D+00
  real    ( kind = 8 ), parameter :: p5 = 0.7003830644436881D+00
  real    ( kind = 8 ), parameter :: p6 = 0.03526249659989109D+00
  real    ( kind = 8 ) pdf
  real    ( kind = 8 ) q
  real    ( kind = 8 ), parameter :: q0 = 440.4137358247522D+00
  real    ( kind = 8 ), parameter :: q1 = 793.8265125199484D+00
  real    ( kind = 8 ), parameter :: q2 = 637.3336333788311D+00
  real    ( kind = 8 ), parameter :: q3 = 296.5642487796737D+00
  real    ( kind = 8 ), parameter :: q4 = 86.78073220294608D+00
  real    ( kind = 8 ), parameter :: q5 = 16.06417757920695D+00
  real    ( kind = 8 ), parameter :: q6 = 1.755667163182642D+00
  real    ( kind = 8 ), parameter :: q7 = 0.08838834764831844D+00
  real    ( kind = 8 ), parameter :: root2pi = 2.506628274631001D+00
  real    ( kind = 8 ) z
  real    ( kind = 8 ) zabs

  zabs = abs ( z )
!
!  37 < |Z|.
!
  if ( 37.0D+00 < zabs ) then

    pdf = 0.0D+00
    p = 0.0D+00
!
!  |Z| <= 37.
!
  else

    expntl = exp ( - 0.5D+00 * zabs * zabs )
    pdf = expntl / root2pi
!
!  |Z| < CUTOFF = 10 / sqrt(2).
!
    if ( zabs < cutoff ) then

      p = expntl * (((((( &
          p6   * zabs &
        + p5 ) * zabs &
        + p4 ) * zabs &
        + p3 ) * zabs &
        + p2 ) * zabs &
        + p1 ) * zabs &
        + p0 ) / ((((((( &
          q7   * zabs &
        + q6 ) * zabs &
        + q5 ) * zabs &
        + q4 ) * zabs &
        + q3 ) * zabs &
        + q2 ) * zabs &
        + q1 ) * zabs &
        + q0 )
!
!  CUTOF <= |Z|.
!
    else

      p = pdf / ( &
        zabs + 1.0D+00 / ( &
        zabs + 2.0D+00 / ( &
        zabs + 3.0D+00 / ( &
        zabs + 4.0D+00 / ( &
        zabs + 0.65D+00 )))))

    end if

  end if

  if ( z < 0.0D+00 ) then
    q = 1.0D+00 - p
  else
    q = p
    p = 1.0D+00 - q
  end if

  return
end
subroutine owen_values ( n_data, h, a, t )

!*****************************************************************************80
!
!! OWEN_VALUES returns some values of Owen's T function.
!
!  Discussion:
!
!    Owen's T function is useful for computation of the bivariate normal
!    distribution and the distribution of a skewed normal distribution.
!
!    Although it was originally formulated in terms of the bivariate
!    normal function, the function can be defined more directly as
!
!      T(H,A) = 1 / ( 2 * pi ) *
!        Integral ( 0 <= X <= A ) e^(-H^2*(1+X^2)/2) / (1+X^2) dX
!
!    In Mathematica, the function can be evaluated by:
!
!      fx = 1/(2*Pi) * Integrate [ E^(-h^2*(1+x^2)/2)/(1+x^2), {x,0,a} ]
!
!  Modified:
!
!    10 December 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Stephen Wolfram,
!    The Mathematica Book,
!    Fourth Edition,
!    Cambridge University Press, 1999,
!    ISBN: 0-521-64314-7,
!    LC: QA76.95.W65.
!
!  Parameters:
!
!    Input/output, integer ( kind = 4 ) N_DATA.  The user sets N_DATA to 0
!    before the first call.  On each call, the routine increments N_DATA by 1,
!    and returns the corresponding data; when there is no more data, the
!    output value of N_DATA will be 0 again.
!
!    Output, real ( kind = 8 ) H, a parameter.
!
!    Output, real ( kind = 8 ) A, the upper limit of the integral.
!
!    Output, real ( kind = 8 ) T, the value of the function.
!
  implicit none

  integer ( kind = 4 ), parameter :: n_max = 22

  real ( kind = 8 ) a
  real ( kind = 8 ), save, dimension ( n_max ) :: a_vec = (/ &
    0.5000000000000000D+00, &
    0.1000000000000000D+01, &
    0.2000000000000000D+01, &
    0.3000000000000000D+01, &
    0.5000000000000000D+00, &
    0.1000000000000000D+01, &
    0.2000000000000000D+01, &
    0.3000000000000000D+01, &
    0.5000000000000000D+00, &
    0.1000000000000000D+01, &
    0.2000000000000000D+01, &
    0.3000000000000000D+01, &
    0.5000000000000000D+00, &
    0.1000000000000000D+01, &
    0.2000000000000000D+01, &
    0.3000000000000000D+01, &
    0.5000000000000000D+00, &
    0.1000000000000000D+01, &
    0.2000000000000000D+01, &
    0.3000000000000000D+01, &
    0.1000000000000000D+02, &
    0.1000000000000000D+03 /)
  real ( kind = 8 ) h
  real ( kind = 8 ), save, dimension ( n_max ) :: h_vec = (/ &
    0.1000000000000000D+01, &
    0.1000000000000000D+01, &
    0.1000000000000000D+01, &
    0.1000000000000000D+01, &
    0.5000000000000000D+00, &
    0.5000000000000000D+00, &
    0.5000000000000000D+00, &
    0.5000000000000000D+00, &
    0.2500000000000000D+00, &
    0.2500000000000000D+00, &
    0.2500000000000000D+00, &
    0.2500000000000000D+00, &
    0.1250000000000000D+00, &
    0.1250000000000000D+00, &
    0.1250000000000000D+00, &
    0.1250000000000000D+00, &
    0.7812500000000000D-02, &
    0.7812500000000000D-02, &
    0.7812500000000000D-02, &
    0.7812500000000000D-02, &
    0.7812500000000000D-02, &
    0.7812500000000000D-02 /)
  integer ( kind = 4 ) n_data
  real ( kind = 8 ) t
  real ( kind = 8 ), save, dimension ( n_max ) :: t_vec = (/ &
    0.4306469112078537D-01, &
    0.6674188216570097D-01, &
    0.7846818699308410D-01, &
    0.7929950474887259D-01, &
    0.6448860284750376D-01, &
    0.1066710629614485D+00, &
    0.1415806036539784D+00, &
    0.1510840430760184D+00, &
    0.7134663382271778D-01, &
    0.1201285306350883D+00, &
    0.1666128410939293D+00, &
    0.1847501847929859D+00, &
    0.7317273327500385D-01, &
    0.1237630544953746D+00, &
    0.1737438887583106D+00, &
    0.1951190307092811D+00, &
    0.7378938035365546D-01, &
    0.1249951430754052D+00, &
    0.1761984774738108D+00, &
    0.1987772386442824D+00, &
    0.2340886964802671D+00, &
    0.2479460829231492D+00 /)

  if ( n_data < 0 ) then
    n_data = 0
  end if

  n_data = n_data + 1

  if ( n_max < n_data ) then
    n_data = 0
    h = 0.0D+00
    a = 0.0D+00
    t = 0.0D+00
  else
    h = h_vec(n_data)
    a = a_vec(n_data)
    t = t_vec(n_data)
  end if

  return
end
function q ( h, ah ) 

!*****************************************************************************80
! 
!! Q computes (1/2) * p(H<Z) - T(H,AH).
!
!  Discussion:
!
!    The routine computes Q = (1/2) * P( H < Z ) - T ( H, AH ).
!
!    The result for Q is non-negative. 
!
!    Warning : Q is computed as the difference between two terms; 
!    When the two terms are of similar value this may produce 
!    error in Q. 
!
!  Modified:
!
!    02 February 2008
!
!  Author:
!
!    Original version by Mike Patefield, David Tandy.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Mike Patefield, David Tandy,
!    Fast and Accurate Calculation of Owen's T Function,
!    Journal of Statistical Software,
!    Volume 5, Number 5, 2000, pages 1-25.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) H, the lower limit for Z.
!    0 < H.
!
!    Input, real ( kind = 8 ) AH, one of the arguments for the T function.
!
!    Output, real ( kind = 8 ) Q, the desired quantity.
!
  implicit none

  real    ( kind = 8 ) ah
  real    ( kind = 8 ) ahh
  real    ( kind = 8 ) h
  integer ( kind = 4 ) ifail
  real    ( kind = 8 ) q
  real    ( kind = 8 ), parameter :: rroot2 = 0.70710678118654752440D+00
  real    ( kind = 8 ) t
  real    ( kind = 8 ) tfun
  real    ( kind = 8 ) x
  real    ( kind = 8 ) znorm1
  real    ( kind = 8 ) znorm2

  if ( 1.0D+00 < ah ) then 
    ahh = ah * h 
    q = tfun ( ahh, 1.0D+00 / ah, h) - znorm2 ( ahh ) * znorm1 ( h ) 
  else 
    q = 0.5D+00 * znorm2 ( h ) - t ( h, ah ) 
  end if 

  return 
end
function t ( h, a ) 

!*****************************************************************************80
!
!! T computes Owen's T function for arbitrary H and A.
!
!  Discussion:
!
!    As printed, this routine had an error in the second case.
!
!  Modified:
!
!    02 February 2008
!
!  Author:
!
!    Original version by Mike Patefield, David Tandy.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Mike Patefield, David Tandy,
!    Fast and Accurate Calculation of Owen's T Function,
!    Journal of Statistical Software,
!    Volume 5, Number 5, 2000, pages 1-25.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) H, A, the arguments.
!
!    Output, real ( kind = 8 ) T, the value of Owen's T function.
!
  implicit none

  real    ( kind = 8 ) a
  real    ( kind = 8 ) absa
  real    ( kind = 8 ) absh
  real    ( kind = 8 ) ah
  real    ( kind = 8 ), parameter :: cut = 0.67D+00
  real    ( kind = 8 ) h
  integer ( kind = 4 ) ifail
  real    ( kind = 8 ) normah
  real    ( kind = 8 ) normh
  real    ( kind = 8 ), parameter :: rroot2 = 0.70710678118654752440D+00
  real    ( kind = 8 ) t
  real    ( kind = 8 ) tfun
  real    ( kind = 8 ) x
  real    ( kind = 8 ) znorm1
  real    ( kind = 8 ) znorm2

  absh = abs ( h ) 
  absa = abs ( a ) 
  ah = absa * absh 

  if ( absa <= 1.0D+00 ) then 

    t = tfun ( absh, absa, ah )
!
!  In the printed paper, the formula for T that follows was incorrect.
!
  else if ( absh <= cut ) then

    t = 0.5D+00 * znorm1 ( absh ) + 0.5D+00 * znorm1 ( ah ) &
    - znorm1 ( absh ) * znorm1 ( ah ) - tfun ( ah, 1.0D+00 / absa, absh ) 

  else

    normh = znorm2 ( absh ) 
    normah = znorm2 ( ah )
    t = 0.5D+00 * ( normh + normah ) - normh * normah &
    - tfun ( ah, 1.0D+00 / absa, absh ) 

  end if
 
  if ( a < 0.0D+00 ) then
    t = - t 
  end if

  return 
end 
function tfun ( h, a, ah ) 

!*****************************************************************************80
!
!! TFUN computes Owen's T function for a restricted range of parameters.
!
!  Discussion:
!
!    This routine computes Owen's T-function of H and A.
!
!    This routine was originally named "TF", but was renamed "TFUN" to
!    avoid a conflict with a built in MATLAB function.
!
!  Modified:
!
!    02 February 2008
!
!  Author:
!
!    Original version by Mike Patefield, David Tandy.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Mike Patefield, David Tandy,
!    Fast and Accurate Calculation of Owen's T Function,
!    Journal of Statistical Software,
!    Volume 5, Number 5, 2000, pages 1-25.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) H, the H argument of the function.
!    0 <= H.
!
!    Input, real ( kind = 8 ) A, the A argument of the function.
!    0 <= A <= 1.
!
!    Input, real ( kind = 8 ) AH, the value of A*H.
!
!    Output, real ( kind = 8 ) TF, the value of Owen's T function.
!
  implicit none

  real    ( kind = 8 ) a
  real    ( kind = 8 ) ah
  real    ( kind = 8 ) ai
  real    ( kind = 8 ) aj
  real    ( kind = 8 ), dimension ( 7 ) :: arange = (/ &
    0.025D+00, 0.09D+00, 0.15D+00, 0.36D+00, 0.5D+00, &
    0.9D+00, 0.99999D+00 /)
  real    ( kind = 8 ) as
  real    ( kind = 8 ), dimension ( 21 ) :: c2 = (/ &
                                   0.99999999999999987510D+00, &
     -0.99999999999988796462D+00,  0.99999999998290743652D+00, &
     -0.99999999896282500134D+00,  0.99999996660459362918D+00, &
     -0.99999933986272476760D+00,  0.99999125611136965852D+00, &
     -0.99991777624463387686D+00,  0.99942835555870132569D+00, &
     -0.99697311720723000295D+00,  0.98751448037275303682D+00, &
     -0.95915857980572882813D+00,  0.89246305511006708555D+00, &
     -0.76893425990463999675D+00,  0.58893528468484693250D+00, &
     -0.38380345160440256652D+00,  0.20317601701045299653D+00, &
     -0.82813631607004984866D-01,  0.24167984735759576523D-01, &
     -0.44676566663971825242D-02,  0.39141169402373836468D-03 /)
  real    ( kind = 8 ) dhs
  real    ( kind = 8 ) dj
  real    ( kind = 8 ) gj
  real    ( kind = 8 ) h
  real    ( kind = 8 ), dimension ( 14 ) :: hrange = (/ &
    0.02D+00, 0.06D+00, 0.09D+00, 0.125D+00, 0.26D+00, &
    0.4D+00,  0.6D+00,  1.6D+00,  1.7D+00,   2.33D+00, &
    2.4D+00,  3.36D+00, 3.4D+00,  4.8D+00 /)
  real    ( kind = 8 ) hs
  integer ( kind = 4 ) i
  integer ( kind = 4 ) iaint
  integer ( kind = 4 ) icode
  integer ( kind = 4 ) ifail
  integer ( kind = 4 ) ihint
  integer ( kind = 4 ) ii
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jj
  integer ( kind = 4 ) m
  integer ( kind = 4 ) maxii
  integer ( kind = 4 ), dimension ( 18 ) :: meth = (/ &
    1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6 /)
  real    ( kind = 8 ) normh
  integer ( kind = 4 ), dimension ( 18 ) :: ord = (/ &
    2, 3, 4, 5, 7,10,12,18,10,20,30,20, 4, 7, 8,20,13, 0 /)
  real    ( kind = 8 ), dimension ( 13 ) :: pts = (/ &
                                   0.35082039676451715489D-02, &
      0.31279042338030753740D-01,  0.85266826283219451090D-01, &
      0.16245071730812277011D+00,  0.25851196049125434828D+00, &
      0.36807553840697533536D+00,  0.48501092905604697475D+00, &
      0.60277514152618576821D+00,  0.71477884217753226516D+00, &
      0.81475510988760098605D+00,  0.89711029755948965867D+00, &
      0.95723808085944261843D+00,  0.99178832974629703586D+00 /)
  real    ( kind = 8 ) r
  real    ( kind = 8 ), parameter :: rroot2 = 0.70710678118654752440D+00
  real    ( kind = 8 ), parameter :: rrtpi = 0.39894228040143267794D+00
  real    ( kind = 8 ), parameter :: rtwopi = 0.15915494309189533577D+00
  integer ( kind = 4 ), dimension ( 15, 8 ) :: select = reshape ( (/&
    1, 1, 2,13,13,13,13,13,13,13,13,16,16,16, 9, &
    1, 2, 2, 3, 3, 5, 5,14,14,15,15,16,16,16, 9, &
    2, 2, 3, 3, 3, 5, 5,15,15,15,15,16,16,16,10, &
    2, 2, 3, 5, 5, 5, 5, 7, 7,16,16,16,16,16,10, &
    2, 3, 3, 5, 5, 6, 6, 8, 8,17,17,17,12,12,11, &
    2, 3, 5, 5, 5, 6, 6, 8, 8,17,17,17,12,12,12, &
    2, 3, 4, 4, 6, 6, 8, 8,17,17,17,17,17,12,12, &
    2, 3, 4, 4, 6, 6,18,18,18,18,17,17,17,12,12 /), (/ 15, 8 /) ) 
  real    ( kind = 8 ) tfun
  real    ( kind = 8 ) vi
  real    ( kind = 8 ), dimension ( 13 ) :: wts = (/ &
                                   0.18831438115323502887D-01, &
      0.18567086243977649478D-01,  0.18042093461223385584D-01, &
      0.17263829606398753364D-01,  0.16243219975989856730D-01, &
      0.14994592034116704829D-01,  0.13535474469662088392D-01, &
      0.11886351605820165233D-01,  0.10070377242777431897D-01, &
      0.81130545742299586629D-02,  0.60419009528470238773D-02, &
      0.38862217010742057883D-02,  0.16793031084546090448D-02 /)
  real    ( kind = 8 ) x
  real    ( kind = 8 ) y
  real    ( kind = 8 ) yi
  real    ( kind = 8 ) z
  real    ( kind = 8 ) zi
  real    ( kind = 8 ) znorm1
  real    ( kind = 8 ) znorm2
! 
!  Determine appropriate method from t1...t6 
!
  ihint = 15

  do i = 1, 14 
    if ( h <= hrange(i) ) then
      ihint = i
      exit
    end if
  end do

  iaint = 8

  do i = 1, 7 
    if ( a <= arange(i) ) then
      iaint = i
      exit
    end if
  end do 

  icode = select(ihint,iaint) 
  m = ord(icode) 
! 
!  t1(h, a, m) ; m = 2, 3, 4, 5, 7, 10, 12 or 18 
!  jj = 2j - 1 ; gj = exp(-h*h/2) * (-h*h/2)**j / j! 
!  aj = a**(2j-1) / (2*pi) 
! 
  if ( meth(icode) == 1 ) then

    hs = - 0.5D+00 * h * h 
    dhs = exp ( hs ) 
    as = a * a 
    j = 1 
    jj = 1 
    aj = rtwopi * a
    tfun = rtwopi * atan ( a ) 
    dj = dhs - 1.0D+00 
    gj = hs * dhs
 
    do

      tfun = tfun + dj * aj / real ( jj, kind = 8 ) 

      if ( m <= j ) then
        return
      end if 

      j = j + 1 
      jj = jj + 2 
      aj = aj * as 
      dj = gj - dj 
      gj = gj * hs / real ( j, kind = 8 ) 

    end  do 
! 
!  t2(h, a, m) ; m = 10, 20 or 30 
!  z = (-1)**(i-1) * zi ; ii = 2i - 1 
!  vi = (-1)**(i-1) * a**(2i-1) * exp[-(a*h)**2/2] / sqrt(2*pi) 
! 
  else if ( meth(icode) == 2 ) then

    maxii = m + m + 1 
    ii = 1 
    tfun = 0.0D+00 
    hs = h * h 
    as = - a * a 
    vi = rrtpi * a * exp ( - 0.5D+00 * ah * ah ) 
    z = znorm1 ( ah ) / h 
    y = 1.0D+00 / hs 

    do

      tfun = tfun + z 

      if ( maxii <= ii ) then
        tfun = tfun * rrtpi * exp ( - 0.5D+00 * hs ) 
        return 
      end if

      z = y * ( vi - real ( ii, kind = 8 ) * z ) 
      vi = as * vi 
      ii = ii + 2 

    end do
! 
!  t3(h, a, m) ; m = 20 
!  ii = 2i - 1 
!  vi = a**(2i-1) * exp[-(a*h)**2/2] / sqrt(2*pi) 
! 
  else if ( meth(icode) == 3 ) then

    i = 1 
    ii = 1 
    tfun = 0.0D+00 
    hs = h * h 
    as = a * a 
    vi = rrtpi * a * exp ( - 0.5D+00 * ah * ah ) 
    zi = znorm1 ( ah ) / h 
    y = 1.0D+00 / hs 

    do

      tfun = tfun + zi * c2(i) 

      if ( m < i ) then
        tfun = tfun * rrtpi * exp ( - 0.5D+00 * hs ) 
        return 
      end if

      zi = y  * ( real ( ii, kind = 8 ) * zi - vi ) 
      vi = as * vi 
      i = i + 1 
      ii = ii + 2 

    end do
! 
!  t4(h, a, m) ; m = 4, 7, 8 or 20;  ii = 2i + 1 
!  ai = a * exp[-h*h*(1+a*a)/2] * (-a*a)**i / (2*pi) 
! 
  else if ( meth(icode) == 4 ) then

    maxii = m + m + 1
    ii = 1 
    hs = h * h
    as = - a * a 
    tfun = 0.0D+00 
    ai = rtwopi * a * exp ( - 0.5D+00 * hs * ( 1.0D+00 - as ) )  
    yi = 1.0D+00

    do

      tfun = tfun + ai * yi 

      if ( maxii <= ii ) then
        return
      end if
 
      ii = ii + 2 
      yi = ( 1.0D+00 - hs * yi ) / real ( ii, kind = 8 ) 
      ai = ai * as 

    end do
! 
!  t5(h, a, m) ; m = 13 
!  2m - point gaussian quadrature 
! 
  else if ( meth(icode) == 5 ) then

    tfun = 0.0D+00 
    as = a * a 
    hs = - 0.5D+00 * h * h 
    do i = 1, m 
      r = 1.0D+00 + as * pts(i) 
      tfun = tfun + wts(i) * exp ( hs * r ) / r 
    end do
    tfun = a * tfun 
! 
!  t6(h, a);  approximation for a near 1, (a<=1) 
! 
  else if ( meth(icode) == 6 ) then

    normh = znorm2 ( h ) 
    tfun = 0.5D+00 * normh * ( 1.0D+00 - normh ) 
    y = 1.0D+00 - a 
    r = atan ( y / ( 1.0D+00 + a ) ) 

    if ( r /= 0.0D+00 ) then
      tfun = tfun - rtwopi * r * exp ( - 0.5D+00 * y * h * h / r ) 
    end if

  end if

  return 
end
subroutine timestamp ( )

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

  character ( len = 8 ) ampm
  integer d
  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
  integer values(8)
  integer y

  call date_and_time ( values = 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 ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end
function znorm1 ( z )

!*****************************************************************************80
!
!! ZNORM1 evaluates the normal CDF from -oo to Z.
!
!  Modified:
!
!    02 February 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) Z, the upper limit.
!
!    Output, real ( kind = 8 ) ZNORM1, the probability that a standard
!    normal variable will lie between -oo and Z.
!
  implicit none

  real    ( kind = 8 ) p
  real    ( kind = 8 ) pdf
  real    ( kind = 8 ) q
  real    ( kind = 8 ) z
  real    ( kind = 8 ) znorm1

  call normp ( z, p, q, pdf )
  znorm1 = p

  return
end
function znorm2 ( z )

!*****************************************************************************80
!
!! ZNORM2 evaluates the normal CDF from Z to +oo.
!
!  Modified:
!
!    02 February 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) Z, the lower limit.
!
!    Output, real ( kind = 8 ) ZNORM2, the probability that a standard
!    normal variable will lie between Z and +oo.
!
  implicit none

  real    ( kind = 8 ) p
  real    ( kind = 8 ) pdf
  real    ( kind = 8 ) q
  real    ( kind = 8 ) z
  real    ( kind = 8 ) znorm2

  call normp ( z, p, q, pdf )
  znorm2 = q

  return
end