function bivprb ( h, k, r ) 

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

      double precision bivprb 
      double precision h
      integer ifail
      double precision k
      double precision q
      double precision r
      double precision ri
      double precision rr
      double precision rroot2
      parameter ( rroot2 = 0.70710678118654752440D+00 )
      double precision rtwopi
      parameter ( rtwopi = 0.15915494309189533577D+00 )
      double precision znorm2

      if ( r .eq. 0.0D+00 ) then

        bivprb = znorm2 ( h ) * znorm2 ( k )

      else 

        rr = 1.0D+00 - r * r 

        if ( 0.0D+00 .lt. rr ) then

          ri = 1.0D+00 / sqrt ( rr )

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

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

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

        else if ( 1.0D+00 .le. r ) then 

          if ( k .le. 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 )

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

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

      zabs = dabs ( z )
c
c  37 < |Z|.
c
      if ( 37.0D+00 .lt. zabs ) then

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

        expntl = dexp ( - 0.5D+00 * zabs * zabs )
        pdf = expntl / root2pi
c
c  |Z| < CUTOFF = 10 / sqrt(2).
c
        if ( zabs .lt. 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 )
c
c  CUTOF <= |Z|.
c
        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 .lt. 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 )

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

      integer n_max
      parameter ( n_max = 22 )

      double precision a
      double precision a_vec(n_max) 
      double precision h
      double precision h_vec(n_max) 
      integer n_data
      double precision t
      double precision t_vec(n_max) 

      save a_vec
      save h_vec
      save t_vec

      data 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 /
      data 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 /
      data 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 .lt. 0 ) then
        n_data = 0
      end if

      n_data = n_data + 1

      if ( n_max .lt. 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 ) 

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

      double precision ah
      double precision ahh
      double precision h
      integer ifail
      double precision q
      double precision rroot2
      parameter ( rroot2 = 0.70710678118654752440D+00 )
      double precision t
      double precision tfun
      double precision x
      double precision znorm1
      double precision znorm2

      if ( 1.0D+00 .lt. 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 ) 

c*********************************************************************72
c
cc T computes Owen's T function for arbitrary H and A.
c
c  Discussion:
c
c    As printed, this routine had an error in the second case.
c
c  Modified:
c
c    16 January 2008
c
c  Author:
c
c    Original version by Mike Patefield, David Tandy.
c    Modifications by John Burkardt.
c
c  Reference:
c
c    Mike Patefield, David Tandy,
c    Fast and Accurate Calculation of Owen's T Function,
c    Journal of Statistical Software,
c    Volume 5, Number 5, 2000, pages 1-25.
c
c  Parameters:
c
c    Input, double precision H, A, the arguments.
c
c    Output, double precision T, the value of Owen's T function.
c
      implicit none

      double precision a
      double precision absa
      double precision absh
      double precision ah
      double precision cut
      parameter ( cut = 0.67D+00 )
      double precision h
      integer ifail
      double precision normah
      double precision normh
      double precision rroot2
      parameter ( rroot2 = 0.70710678118654752440D+00 )
      double precision t
      double precision tfun
      double precision x
      double precision znorm1
      double precision znorm2

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

      if ( absa .le. 1.0D+00 ) then 

        t = tfun ( absh, absa, ah )
c
c  In the printed paper, the formula for T that follows was incorrect.
c
      else if ( absh .le. 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 .lt. 0.0D+00 ) then
        t = - t 
      end if

      return 
      end 
      function tfun ( h, a, ah ) 

c*********************************************************************72
c
cc TFUN computes Owen's T function for a restricted range of parameters.
c
c  Discussion:
c
c    This routine computes Owen's T-function of H and A.
c
c    Originally named "TF", this routine was renamed "TFUN" to avoid a
c    conflict with a MATLAB built in function.
c
c  Modified:
c
c    16 January 2008
c
c  Author:
c
c    Original version by Mike Patefield, David Tandy.
c    Modifications by John Burkardt.
c
c  Reference:
c
c    Mike Patefield, David Tandy,
c    Fast and Accurate Calculation of Owen's T Function,
c    Journal of Statistical Software,
c    Volume 5, Number 5, 2000, pages 1-25.
c
c  Parameters:
c
c    Input, double precision H, the H argument of the function.
c    0 <= H.
c
c    Input, double precision A, the A argument of the function.
c    0 <= A <= 1.
c
c    Input, double precision AH, the value of A*H.
c
c    Output, double precision TF, the value of Owen's T function.
c
      implicit none

      double precision a
      double precision ah
      double precision ai
      double precision aj
      double precision arange(7)
      double precision as
      double precision c2(21)
      double precision dhs
      double precision dj
      double precision gj
      double precision h
      double precision hrange(14)
      double precision hs
      integer i
      integer iaint
      integer icode
      integer ifail
      integer ihint
      integer ii
      integer j
      integer jj
      integer m
      integer maxii
      integer meth(18)
      double precision normh
      integer ord(18)
      double precision pts(13)
      double precision r
      double precision rroot2
      parameter ( rroot2 = 0.70710678118654752440D+00 )
      double precision rrtpi
      parameter ( rrtpi = 0.39894228040143267794D+00 )
      double precision rtwopi
      parameter ( rtwopi = 0.15915494309189533577D+00 )
      integer select(15,8)
      double precision tfun
      double precision vi
      double precision wts(13)
      double precision x
      double precision y
      double precision yi
      double precision z
      double precision zi
      double precision znorm1
      double precision znorm2

      data arange /
     &  0.025D+00, 0.09D+00, 0.15D+00, 0.36D+00, 0.5D+00, 
     &  0.9D+00, 0.99999D+00 / 

      data 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 / 

      data 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 / 

      data meth / 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6 / 

      data ord /  2, 3, 4, 5, 7,10,12,18,10,20,30,20, 4, 7, 8,20,13, 0 / 

      data 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 / 

      data select / 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 / 

      data 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 / 
c 
c  Determine appropriate method from t1...t6 
c
      ihint = 15

      do i = 1, 14 
        if ( h .le. hrange(i) ) then
          ihint = i
          go to 20 
        end if
      end do

   20 continue

      iaint = 8

      do i = 1, 7 
        if ( a .le. arange(i) ) then
          iaint = i
          go to 40 
        end if
      end do 

   40 continue

      icode = select(ihint,iaint) 
      m = ord(icode) 
c 
c  t1(h, a, m) ; m = 2, 3, 4, 5, 7, 10, 12 or 18 
c  jj = 2j - 1 ; gj = exp(-h*h/2) * (-h*h/2)**j / j! 
c  aj = a**(2j-1) / (2*pi) 
c 
      if ( meth(icode) .eq. 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
 
  110   continue

        tfun = tfun + dj * aj / dble ( jj ) 

        if ( m .le. j ) then
          return
        end if 

        j = j + 1 
        jj = jj + 2 
        aj = aj * as 
        dj = gj - dj 
        gj = gj * hs / dble ( j ) 
        go to 110 
c 
c  t2(h, a, m) ; m = 10, 20 or 30 
c  z = (-1)**(i-1) * zi ; ii = 2i - 1 
c  vi = (-1)**(i-1) * a**(2i-1) * exp[-(a*h)**2/2] / sqrt(2*pi) 
c 
      else if ( meth(icode) .eq. 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 

  210   continue

        tfun = tfun + z 

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

        z = y * ( vi - dble ( ii ) * z ) 
        vi = as * vi 
        ii = ii + 2 
        go to 210 
c 
c  t3(h, a, m) ; m = 20 
c  ii = 2i - 1 
c  vi = a**(2i-1) * exp[-(a*h)**2/2] / sqrt(2*pi) 
c 
      else if ( meth(icode) .eq. 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 

  310   continue

        tfun = tfun + zi * c2(i) 

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

        zi = y  * ( dble ( ii ) * zi - vi ) 
        vi = as * vi 
        i = i + 1 
        ii = ii + 2 
        go to 310 
c 
c  t4(h, a, m) ; m = 4, 7, 8 or 20;  ii = 2i + 1 
c  ai = a * exp[-h*h*(1+a*a)/2] * (-a*a)**i / (2*pi) 
c 
      else if ( meth(icode) .eq. 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

  410   continue

        tfun = tfun + ai * yi 

        if ( maxii .le. ii ) then
          return
        end if
 
        ii = ii + 2 
        yi = ( 1.0D+00 - hs * yi ) / dble ( ii ) 
        ai = ai * as 
        go to 410 
c 
c  t5(h, a, m) ; m = 13 
c  2m - point gaussian quadrature 
c 
      else if ( meth(icode) .eq. 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 
c 
c  t6(h, a);  approximation for a near 1, (a.le.1) 
c 
      else if ( meth(icode) .eq. 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 .ne. 0.0D+00 ) then
          tfun = tfun - rtwopi * r * exp ( - 0.5D+00 * y * h * h / r ) 
        end if

      end if

      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    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
      function znorm1 ( z )

c*********************************************************************72
c
cc ZNORM1 evaluates the normal CDF from -oo to Z.
c
c  Modified:
c
c    16 January 2008
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, double precision Z, the upper limit.
c
c    Output, double precision ZNORM1, the probability that a standard
c    normal variable will lie between -oo and Z.
c
      implicit none

      double precision p
      double precision pdf
      double precision q
      double precision z
      double precision znorm1

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

      return
      end
      function znorm2 ( z )

c*********************************************************************72
c
cc ZNORM2 evaluates the normal CDF from Z to +oo.
c
c  Modified:
c
c    16 January 2008
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, double precision Z, the lower limit.
c
c    Output, double precision ZNORM2, the probability that a standard
c    normal variable will lie between Z and +oo.
c
      implicit none

      double precision p
      double precision pdf
      double precision q
      double precision z
      double precision znorm2

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

      return
      end