ieva.f90

FDS为火灾动力学模拟软件源代码,该软件为开源项目,代码语言主要为FORTRAN,可在WINDOWS和LINUX下编译运行,详细说明可参考http://fire.nist.gov/fds/官方网址

    !  Parameters:
    !
    !    Input, real ( kind = 8 ) A, B, the parameters of the function.
    !    A and B should be nonnegative.
    !
    !    Input, real ( kind = 8 ) X, Y, ?
    !
    !    Input, integer N, ?
    !
    !    Input, real ( kind = 8 ) EPS, the tolerance.
    !
    !    Output, real ( kind = 8 ) BETA_UP, the value of IX(A,B) - IX(A+N,B).
    !
    Implicit None

    Real ( kind = 8 ) a
    Real ( kind = 8 ) ap1
    Real ( kind = 8 ) apb
    Real ( kind = 8 ) b
    !    Real ( kind = 8 ) beta_rcomp1
    Real ( kind = 8 ) beta_up
    Real ( kind = 8 ) d
    Real ( kind = 8 ) eps
    !    Real ( kind = 8 ) exparg
    Integer i
    Integer k
    Real ( kind = 8 ) l
    Integer mu
    Integer n
    Real ( kind = 8 ) r
    Real ( kind = 8 ) t
    Real ( kind = 8 ) w
    Real ( kind = 8 ) x
    Real ( kind = 8 ) y
    !
    !  Obtain the scaling factor EXP(-MU) AND
    !  EXP(MU)*(X**A*Y**B/BETA(A,B))/A
    !
    apb = a + b
    ap1 = a + 1.0D+00
    mu = 0
    d = 1.0D+00

    If ( n /= 1 ) Then

       If ( 1.0D+00 <= a ) Then

          If ( 1.1D+00 * ap1 <= apb ) Then
             mu = Abs ( exparg ( 1 ) )
             k = exparg ( 0 )
             If ( k < mu ) Then
                mu = k
             End If
             t = mu
             d = Exp ( - t )
          End If

       End If

    End If

    beta_up = beta_rcomp1 ( mu, a, b, x, y ) / a

    If ( n == 1 .Or. beta_up == 0.0D+00 ) Then
       Return
    End If

    w = d
    !
    !  Let K be the index of the maximum term.
    !
    k = 0

    If ( 1.0D+00 < b ) Then

       If ( y <= 0.0001D+00 ) Then

          k = n - 1

       Else

          r = ( b - 1.0D+00 ) * x / y - a

          If ( 1.0D+00 <= r ) Then
             k = n - 1
             t = n - 1
             If ( r < t ) Then
                k = r
             End If
          End If

       End If
       !
       !  Add the increasing terms of the series.
       !
       Do i = 1, k
          l = i - 1
          d = ( ( apb + l ) / ( ap1 + l ) ) * x * d
          w = w + d
       End Do

    End If
    !
    !  Add the remaining terms of the series.
    !
    Do i = k+1, n-1
       l = i - 1
       d = ( ( apb + l ) / ( ap1 + l ) ) * x * d
       w = w + d
       If ( d <= eps * w ) Then
          beta_up = beta_up * w
          Return
       End If
    End Do

    beta_up = beta_up * w

    Return
  End Function beta_up

  Subroutine binomial_cdf_values ( n_data, a, b, x, fx )

    !*****************************************************************************80
    !
    !! BINOMIAL_CDF_VALUES returns some values of the binomial CDF.
    !
    !  Discussion:
    !
    !    CDF(X)(A,B) is the probability of at most X successes in A trials,
    !    given that the probability of success on a single trial is B.
    !
    !  Modified:
    !
    !    27 May 2001
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Milton Abramowitz, Irene Stegun,
    !    Handbook of Mathematical Functions,
    !    US Department of Commerce, 1964.
    !
    !    Daniel Zwillinger,
    !    CRC Standard Mathematical Tables and Formulae,
    !    30th Edition, CRC Press, 1996, pages 651-652.
    !
    !  Parameters:
    !
    !    Input/output, integer 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, integer A, real ( kind = 8 ) B, integer X, the arguments 
    !    of the function.
    !
    !    Output, real ( kind = 8 ) FX, the value of the function.
    !
    Implicit None

    Integer, Parameter :: n_max = 17

    Integer a
    Integer, Save, Dimension ( n_max ) :: a_vec = (/ &
         2,  2,  2,  2, &
         2,  4,  4,  4, &
         4, 10, 10, 10, &
         10, 10, 10, 10, &
         10 /)
    Real ( kind = 8 ) b
    Real ( kind = 8 ), Save, Dimension ( n_max ) :: b_vec = (/ &
         0.05D+00, 0.05D+00, 0.05D+00, 0.50D+00, &
         0.50D+00, 0.25D+00, 0.25D+00, 0.25D+00, &
         0.25D+00, 0.05D+00, 0.10D+00, 0.15D+00, &
         0.20D+00, 0.25D+00, 0.30D+00, 0.40D+00, &
         0.50D+00 /)
    Real ( kind = 8 ) fx
    Real ( kind = 8 ), Save, Dimension ( n_max ) :: fx_vec = (/ &
         0.9025D+00, 0.9975D+00, 1.0000D+00, 0.2500D+00, &
         0.7500D+00, 0.3164D+00, 0.7383D+00, 0.9492D+00, &
         0.9961D+00, 0.9999D+00, 0.9984D+00, 0.9901D+00, &
         0.9672D+00, 0.9219D+00, 0.8497D+00, 0.6331D+00, &
         0.3770D+00 /)
    Integer n_data
    Integer x
    Integer, Save, Dimension ( n_max ) :: x_vec = (/ &
         0, 1, 2, 0, &
         1, 0, 1, 2, &
         3, 4, 4, 4, &
         4, 4, 4, 4, &
         4 /)

    If ( n_data < 0 ) Then
       n_data = 0
    End If

    n_data = n_data + 1

    If ( n_max < n_data ) Then
       n_data = 0
       a = 0
       b = 0.0D+00
       x = 0
       fx = 0.0D+00
    Else
       a = a_vec(n_data)
       b = b_vec(n_data)
       x = x_vec(n_data)
       fx = fx_vec(n_data)
    End If

    Return
  End Subroutine binomial_cdf_values

  Subroutine cdfbet ( which, p, q, x, y, a, b, status, bound )

    !*****************************************************************************80
    !
    !! CDFBET evaluates the CDF of the Beta Distribution.
    !
    !  Discussion:
    !
    !    This routine calculates any one parameter of the beta distribution 
    !    given the others.
    !
    !    The value P of the cumulative distribution function is calculated 
    !    directly by code associated with the reference.
    !
    !    Computation of the other parameters involves a seach for a value that
    !    produces the desired value of P.  The search relies on the
    !    monotonicity of P with respect to the other parameters.
    !
    !    The beta density is proportional to t^(A-1) * (1-t)^(B-1).
    !
    !  Modified:
    !
    !    09 June 2004
    !
    !  Reference:
    !
    !    Armido DiDinato, Alfred Morris,
    !    Algorithm 708: 
    !    Significant Digit Computation of the Incomplete Beta Function Ratios,
    !    ACM Transactions on Mathematical Software, 
    !    Volume 18, 1993, pages 360-373.
    !
    !  Parameters:
    !
    !    Input, integer WHICH, indicates which of the next four argument
    !    values is to be calculated from the others.
    !    1: Calculate P and Q from X, Y, A and B;
    !    2: Calculate X and Y from P, Q, A and B;
    !    3: Calculate A from P, Q, X, Y and B;
    !    4: Calculate B from P, Q, X, Y and A.
    !
    !    Input/output, real ( kind = 8 ) P, the integral from 0 to X of the
    !    chi-square distribution.  Input range: [0, 1].
    !
    !    Input/output, real ( kind = 8 ) Q, equals 1-P.  Input range: [0, 1].
    !
    !    Input/output, real ( kind = 8 ) X, the upper limit of integration 
    !    of the beta density.  If it is an input value, it should lie in
    !    the range [0,1].  If it is an output value, it will be searched for
    !    in the range [0,1].
    !
    !    Input/output, real ( kind = 8 ) Y, equal to 1-X.  If it is an input
    !    value, it should lie in the range [0,1].  If it is an output value,
    !    it will be searched for in the range [0,1].
    !
    !    Input/output, real ( kind = 8 ) A, the first parameter of the beta
    !    density.  If it is an input value, it should lie in the range
    !    (0, +infinity).  If it is an output value, it will be searched
    !    for in the range [1D-300,1D300].
    !
    !    Input/output, real ( kind = 8 ) B, the second parameter of the beta
    !    density.  If it is an input value, it should lie in the range
    !    (0, +infinity).  If it is an output value, it will be searched
    !    for in the range [1D-300,1D300].
    !
    !    Output, integer STATUS, reports the status of the computation.
    !     0, if the calculation completed correctly;
    !    -I, if the input parameter number I is out of range;
    !    +1, if the answer appears to be lower than lowest search bound;
    !    +2, if the answer appears to be higher than greatest search bound;
    !    +3, if P + Q /= 1;
    !    +4, if X + Y /= 1.
    !
    !    Output, real ( kind = 8 ) BOUND, is only defined if STATUS is nonzero.
    !    If STATUS is negative, then this is the value exceeded by parameter I.
    !    if STATUS is 1 or 2, this is the search bound that was exceeded.
    !
    Implicit None

    Real ( kind = 8 ) a
    Real ( kind = 8 ), Parameter :: atol = 1.0D-50
    Real ( kind = 8 ) b
    Real ( kind = 8 ) bound
    Real ( kind = 8 ) ccum
    Real ( kind = 8 ) cum
    Real ( kind = 8 ) fx
    Real ( kind = 8 ), Parameter :: inf = 1.0D+300
    Real ( kind = 8 ) p
    Real ( kind = 8 ) q
    Logical qhi
    Logical qleft
    Integer status
    Real ( kind = 8 ), Parameter :: tol = 1.0D-08
    Integer which
    Real ( kind = 8 ) x
    Real ( kind = 8 ) xhi
    Real ( kind = 8 ) xlo
    Real ( kind = 8 ) y

    If ( which < 1 ) Then
       bound = 1.0D+00
       status = -1
       Write ( *, '(a)' ) ' '
       Write ( *, '(a)' ) 'CDFBET - Fatal error!'
       Write ( *, '(a)' ) '  Input parameter WHICH is out of range.'
       Return
    End If

    If ( 4 < which ) Then
       bound = 4.0D+00
       status = -1
       Write ( *, '(a)' ) ' '
       Write ( *, '(a)' ) 'CDFBET - Fatal error!'
       Write ( *, '(a)' ) '  Input parameter WHICH is out of range.'
       Return
    End If
    !
    !  Unless P is to be computed, make sure it is legal.
    !
    If ( which /= 1 ) Then
       If ( p < 0.0D+00 ) Then
          bound = 0.0D+00
          status = -2
          Write ( *, '(a)' ) ' '
          Write ( *, '(a)' ) 'CDFBET - Fatal error!'
          Write ( *, '(a)' ) '  Input parameter P is out of range.'
          Return
       Else If ( 1.0D+00 < p ) Then
          bound = 1.0D+00
          status = -2
          Write ( *, '(a)' ) ' '
          Write ( *, '(a)' ) 'CDFBET - Fatal error!'
          Write ( *, '(a)' ) '  Input parameter P is out of range.'
          Return
       End If
    End If
    !
    !  Unless Q is to be computed, make sure it is legal.
    !
    If ( which /= 1 ) Then
       If ( q < 0.0D+00 ) Then
          bound = 0.0D+00
          status = -3
          Write ( *, '(a)' ) ' '
          Write ( *, '(a)' ) 'CDFBET - Fatal error!'
          Write ( *, '(a)' ) '  Input parameter Q is out of range.'
          Return
       Else If ( 1.0D+00 < q ) Then
          bound = 1.0D+00
          status = -3
          Write ( *, '(a)' ) ' '
          Write ( *, '(a)' ) 'CDFBET - Fatal error!'
          Write ( *, '(a)' ) '  Input parameter Q is out of range.'
          Return
       End If
    End If
    !
    !  Unless X is to be computed, make sure it is legal.
    !
    If ( which /= 2 ) Then
       If ( x < 0.0D+00 ) Then
          bound = 0.0D+00
          status = -4
          Write ( *, '(a)' ) ' '
          Write ( *, '(a)' ) 'CDFBET - Fatal error!'
          Write ( *, '(a)' ) '  Input parameter X is out of range.'
          Return
       Else If ( 1.0D+00 < x ) Then
          bound = 1.0D+00
          status = -4
          Write ( *, '(a)' ) ' '
          Write ( *, '(a)' ) 'CDFBET - Fatal error!'
          Write ( *, '(a)' ) '  Input parameter X is out of range.'
          Return
       End If
    End If
    !
    !  Unless Y is to be computed, make sure it is legal.
    !
    If ( which /= 2 ) Then
       If ( y < 0.0D+00 ) Then
          bound = 0.0D+00
          status = -5
          Write ( *, '(a)' ) ' '
          Write ( *, '(a)' ) 'CDFBET - Fatal error!'
          Write ( *, '(a)' ) '  Input parameter Y is out of range.'
          Return
       Else If ( 1.0D+00 < y ) Then
          bound = 1.0D+00
          status = -5
          Write ( *, '(a)' ) ' '
          Write ( *, '(a)' ) 'CDFBET - Fatal error!'
          Write ( *, '(a)' ) '  Input parameter Y is out of range.'
          Return
       End If
    End If
    !
    !  Unless A is to be computed, make sure it is legal.
    !
    If ( which /= 3 ) Then
       If ( a <= 0.0D+00 ) Then
          bound = 0.0D+00
          status = -6
          Write ( *, '(a)' ) ' '
          Write ( *, '(a)' ) 'CDFBET - Fatal error!'
          Write ( *, '(a)' ) '  Input parameter A is out of range.'
          Return
       End If
    End If
    !
    !  Unless B is to be computed, make sure it is legal.
    !
    If ( which /= 4 ) Then
       If ( b <= 0.0D+00 ) Then
          bound = 0.0D+00
          status = -7
          Write ( *, '(a)' ) ' '
          Write ( *, '(a)' ) 'CDFBET - Fatal error!'
          Write ( *, '(a)' ) '  Input parameter B is out of range.'
          Return
       End If
    End If
    !
    !  Check that P + Q = 1.
    !
    If ( which /= 1 ) Then
       If ( 3.0D+00 * Epsilon ( p ) < Abs ( ( p + q ) - 1.0D+00 ) ) Then
          status = 3
          Write ( *, '(a)' ) ' '
          Write ( *, '(a)' ) 'CDFBET - Fatal error!'
          Write ( *, '(a)' ) '  P + Q /= 1.'
          Return
       End If
    End If
    !
    !  Check that X + Y = 1.
    !
    If ( which /= 2 ) Then
       If ( 3.0D+00 * Epsilon ( x ) < Abs ( ( x + y ) - 1.0D+00 ) ) Then
          status = 4
          Write ( *, '(a)' ) ' '
          Write ( *, '(a)' ) 'CDFBET - Fatal error!'
          Write ( *, '(a)' ) '  X + Y /= 1.'
          Return
       End If
    End If
    !
    !  Compute P and Q.
    !
    If ( which == 1 ) Then

       Call cumbet ( x, y, a, b, p, q )
       status = 0
       !
       !  Compute X and Y.
       !
    Else If ( which == 2 ) Then

       Call dstzr ( 0.0D+00, 1.0D+00, atol, tol )

       If ( p <= q ) Then

          status = 0
          fx = 0.0D+00
          Call dzror ( status, x, fx, xlo, xhi, qleft, qhi )
          y = 1.0D+00 - x

          Do

             If ( status /= 1 ) Then
                Exit
             End If

             Call cumbet ( x, y, a, b, cum, ccum )
             fx = cum - p
             Call dzror ( status, x, fx, xlo, xhi, qleft, qhi )
             y = 1.0D+00 - x

          End Do

       Else

          status = 0
          fx = 0.0D+00
          Call dzror ( status, y, fx, xlo, xhi, qleft, qhi )
          x = 1.0D+00 - y

          Do

             If ( status /= 1 ) Then
                Exit
             End If

             Call cumbet ( x, y, a, b, cum, ccum )
             fx = ccum - q
             Call dzror ( status, y, fx, xlo, xhi, qleft, qhi )
             x = 1.0D+00 - y

          End Do

       End If

       If ( status == -1 ) Then
          If ( qleft ) Then
             status = 1
             bound = 0.0D+00
             Write ( *, '(