ieva.f90

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

    End If
    !
    !  Procedure for A0 <= 1 or B0 <= 1
    !
    If ( 0.5D+00 < x ) Then
       ind = 1
       a0 = b
       b0 = a
       x0 = y
       y0 = x
    End If

    If ( b0 < Min ( eps, eps * a0 ) ) Then
       go to 90
    End If

    If ( a0 < Min ( eps, eps * b0 ) .And. b0 * x0 <= 1.0D+00 ) Then
       go to 100
    End If

    If ( 1.0D+00 < Max ( a0, b0 ) ) Then
       go to 20
    End If

    If ( Min ( 0.2D+00, b0 ) <= a0 ) Then
       go to 110
    End If

    If ( x0**a0 <= 0.9D+00 ) Then
       go to 110
    End If

    If ( 0.3D+00 <= x0 ) Then
       go to 120
    End If

    n = 20
    go to 140

20  Continue

    If ( b0 <= 1.0D+00 ) Then
       go to 110
    End If

    If ( 0.3D+00 <= x0 ) Then
       go to 120
    End If

    If ( 0.1D+00 <= x0 ) Then
       go to 30
    End If

    If ( (x0 * b0 )**a0 <= 0.7D+00 ) Then
       go to 110
    End If

30  Continue

    If ( 15.0D+00 < b0 ) Then
       go to 150
    End If

    n = 20
    go to 140
    !
    !  PROCEDURE for 1 < A0 and 1 < B0.
    !
40  Continue

    If ( a <= b ) Then
       lambda = a - ( a + b ) * x
    Else
       lambda = ( a + b ) * y - b
    End If

    If ( lambda < 0.0D+00 ) Then
       ind = 1
       a0 = b
       b0 = a
       x0 = y
       y0 = x
       lambda = Abs ( lambda )
    End If

70  Continue

    If ( b0 < 40.0D+00 .And. b0 * x0 <= 0.7D+00 ) Then
       go to 110
    End If

    If ( b0 < 40.0D+00 ) Then
       go to 160
    End If

    If ( b0 < a0 ) Then
       go to 80
    End If

    If ( a0 <= 100.0D+00 ) Then
       go to 130
    End If

    If ( 0.03D+00 * a0 < lambda ) Then
       go to 130
    End If

    go to 200

80  Continue

    If ( b0 <= 100.0D+00 ) Then
       go to 130
    End If

    If ( 0.03D+00 * b0 < lambda ) Then
       go to 130
    End If

    go to 200
    !
    !  Evaluation of the appropriate algorithm.
    !
90  Continue

    w = fpser ( a0, b0, x0, eps )
    w1 = 0.5D+00 + ( 0.5D+00 - w )
    go to 250

100 Continue

    w1 = apser ( a0, b0, x0, eps )
    w = 0.5D+00 + ( 0.5D+00 - w1 )
    go to 250

110 Continue

    w = beta_pser ( a0, b0, x0, eps )
    w1 = 0.5D+00 + ( 0.5D+00 - w )
    go to 250

120 Continue

    w1 = beta_pser ( b0, a0, y0, eps )
    w = 0.5D+00 + ( 0.5D+00 - w1 )
    go to 250

130 Continue

    w = beta_frac ( a0, b0, x0, y0, lambda, 15.0D+00 * eps )
    w1 = 0.5D+00 + ( 0.5D+00 - w )
    go to 250

140 Continue

    w1 = beta_up ( b0, a0, y0, x0, n, eps )
    b0 = b0 + n

150 Continue

    Call beta_grat ( b0, a0, y0, x0, w1, 15.0D+00 * eps, ierr1 )
    w = 0.5D+00 + ( 0.5D+00 - w1 )
    go to 250

160 Continue

    n = b0
    b0 = b0 - n

    If ( b0 == 0.0D+00 ) Then
       n = n - 1
       b0 = 1.0D+00
    End If

170 Continue

    w = beta_up ( b0, a0, y0, x0, n, eps )

    If ( x0 <= 0.7D+00 ) Then
       w = w + beta_pser ( a0, b0, x0, eps )
       w1 = 0.5D+00 + ( 0.5D+00 - w )
       go to 250
    End If

    If ( a0 <= 15.0D+00 ) Then
       n = 20
       w = w + beta_up ( a0, b0, x0, y0, n, eps )
       a0 = a0 + n
    End If

190 Continue

    Call beta_grat ( a0, b0, x0, y0, w, 15.0D+00 * eps, ierr1 )
    w1 = 0.5D+00 + ( 0.5D+00 - w )
    go to 250

200 Continue

    w = beta_asym ( a0, b0, lambda, 100.0D+00 * eps )
    w1 = 0.5D+00 + ( 0.5D+00 - w )
    go to 250
    !
    !  Termination of the procedure.
    !
250 Continue

    If ( ind /= 0 ) Then
       t = w
       w = w1
       w1 = t
    End If

    Return
    !
    !  Procedure for A and B < 0.001 * EPS
    !
260 Continue

    w = b / ( a + b )
    w1 = a / ( a + b )

    Return
  End Subroutine beta_inc

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

    !*****************************************************************************80
    !
    !! BETA_INC_VALUES returns some values of the incomplete Beta function.
    !
    !  Discussion:
    !
    !    The incomplete Beta function may be written
    !
    !      BETA_INC(A,B,X) = Integral (0 to X) T**(A-1) * (1-T)**(B-1) dT
    !                      / Integral (0 to 1) T**(A-1) * (1-T)**(B-1) dT
    !
    !    Thus,
    !
    !      BETA_INC(A,B,0.0) = 0.0
    !      BETA_INC(A,B,1.0) = 1.0
    !
    !    Note that in Mathematica, the expressions:
    !
    !      BETA[A,B]   = Integral (0 to 1) T**(A-1) * (1-T)**(B-1) dT
    !      BETA[X,A,B] = Integral (0 to X) T**(A-1) * (1-T)**(B-1) dT
    !
    !    and thus, to evaluate the incomplete Beta function requires:
    !
    !      BETA_INC(A,B,X) = BETA[X,A,B] / BETA[A,B]
    !
    !  Modified:
    !
    !    17 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Milton Abramowitz, Irene Stegun,
    !    Handbook of Mathematical Functions,
    !    US Department of Commerce, 1964.
    !
    !    Karl Pearson,
    !    Tables of the Incomplete Beta Function,
    !    Cambridge University Press, 1968.
    !
    !  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, real ( kind = 8 ) A, B, X, the arguments of the function.
    !
    !    Output, real ( kind = 8 ) FX, the value of the function.
    !
    Implicit None

    Integer, Parameter :: n_max = 30

    Real ( kind = 8 ) a
    Real ( kind = 8 ), Save, Dimension ( n_max ) :: a_vec = (/ &
         0.5D+00,  0.5D+00,  0.5D+00,  1.0D+00, &
         1.0D+00,  1.0D+00,  1.0D+00,  1.0D+00, &
         2.0D+00,  2.0D+00,  2.0D+00,  2.0D+00, &
         2.0D+00,  2.0D+00,  2.0D+00,  2.0D+00, &
         2.0D+00,  5.5D+00, 10.0D+00, 10.0D+00, &
         10.0D+00, 10.0D+00, 20.0D+00, 20.0D+00, &
         20.0D+00, 20.0D+00, 20.0D+00, 30.0D+00, &
         30.0D+00, 40.0D+00 /)
    Real ( kind = 8 ) b
    Real ( kind = 8 ), Save, Dimension ( n_max ) :: b_vec = (/ &
         0.5D+00,  0.5D+00,  0.5D+00,  0.5D+00, &
         0.5D+00,  0.5D+00,  0.5D+00,  1.0D+00, &
         2.0D+00,  2.0D+00,  2.0D+00,  2.0D+00, &
         2.0D+00,  2.0D+00,  2.0D+00,  2.0D+00, &
         2.0D+00,  5.0D+00,  0.5D+00,  5.0D+00, &
         5.0D+00, 10.0D+00,  5.0D+00, 10.0D+00, &
         10.0D+00, 20.0D+00, 20.0D+00, 10.0D+00, &
         10.0D+00, 20.0D+00 /)
    Real ( kind = 8 ) fx
    Real ( kind = 8 ), Save, Dimension ( n_max ) :: fx_vec = (/ &
         0.0637686D+00, 0.2048328D+00, 1.0000000D+00, 0.0D+00,       &
         0.0050126D+00, 0.0513167D+00, 0.2928932D+00, 0.5000000D+00, &
         0.028D+00,     0.104D+00,     0.216D+00,     0.352D+00,     &
         0.500D+00,     0.648D+00,     0.784D+00,     0.896D+00,     &
         0.972D+00,     0.4361909D+00, 0.1516409D+00, 0.0897827D+00, &
         1.0000000D+00, 0.5000000D+00, 0.4598773D+00, 0.2146816D+00, &
         0.9507365D+00, 0.5000000D+00, 0.8979414D+00, 0.2241297D+00, &
         0.7586405D+00, 0.7001783D+00 /)
    Integer n_data
    Real ( kind = 8 ) x
    Real ( kind = 8 ), Save, Dimension ( n_max ) :: x_vec = (/ &
         0.01D+00, 0.10D+00, 1.00D+00, 0.0D+00,  &
         0.01D+00, 0.10D+00, 0.50D+00, 0.50D+00, &
         0.1D+00,  0.2D+00,  0.3D+00,  0.4D+00,  &
         0.5D+00,  0.6D+00,  0.7D+00,  0.8D+00,  &
         0.9D+00,  0.50D+00, 0.90D+00, 0.50D+00, &
         1.00D+00, 0.50D+00, 0.80D+00, 0.60D+00, &
         0.80D+00, 0.50D+00, 0.60D+00, 0.70D+00, &
         0.80D+00, 0.70D+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
       a = 0.0D+00
       b = 0.0D+00
       x = 0.0D+00
       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 beta_inc_values

  Function beta_log ( a0, b0 )

    !*****************************************************************************80
    !
    !! BETA_LOG evaluates the logarithm of the beta function.
    !
    !  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, real ( kind = 8 ) A0, B0, the parameters of the function.
    !    A0 and B0 should be nonnegative.
    !
    !    Output, real ( kind = 8 ) BETA_LOG, the value of the logarithm
    !    of the Beta function.
    !
    Implicit None

    Real ( kind = 8 ) a
    Real ( kind = 8 ) a0
    !    Real ( kind = 8 ) algdiv
    !    Real ( kind = 8 ) alnrel
    Real ( kind = 8 ) b
    Real ( kind = 8 ) b0
    !    Real ( kind = 8 ) bcorr
    Real ( kind = 8 ) beta_log
    Real ( kind = 8 ) c
    Real ( kind = 8 ), Parameter :: e = 0.918938533204673D+00
    !    Real ( kind = 8 ) gamma_log
    !    Real ( kind = 8 ) gsumln
    Real ( kind = 8 ) h
    Integer i
    Integer n
    Real ( kind = 8 ) u
    Real ( kind = 8 ) v
    Real ( kind = 8 ) w
    Real ( kind = 8 ) z

    a = Min ( a0, b0 )
    b = Max ( a0, b0 )

    If ( 8.0D+00 <= a ) Then
       go to 100
    End If

    If ( 1.0D+00 <= a ) Then
       go to 20
    End If
    !
    !  Procedure when A < 1
    !
    If ( b < 8.0D+00 ) Then
       beta_log = gamma_log ( a ) + ( gamma_log ( b ) - gamma_log ( a + b ) )
    Else
       beta_log = gamma_log ( a ) + algdiv ( a, b )
    End If

    Return
    !
    !  Procedure when 1 <= A < 8
    !
20  Continue

    If ( 2.0D+00 < a ) Then
       go to 40
    End If

    If ( b <= 2.0D+00 ) Then
       beta_log = gamma_log ( a ) + gamma_log ( b ) - gsumln ( a, b )
       Return
    End If

    w = 0.0D+00
    If ( b < 8.0D+00 ) Then
       go to 60
    End If

    beta_log = gamma_log ( a ) + algdiv ( a, b )
    Return
    !
    !  Reduction of A when B <= 1000.
    !
40  Continue

    If ( 1000.0D+00 < b ) Then
       go to 80
    End If

    n = a - 1.0D+00
    w = 1.0D+00
    Do i = 1, n
       a = a - 1.0D+00
       h = a / b
       w = w * ( h / ( 1.0D+00 + h ) )
    End Do
    w = Log ( w )

    If ( 8.0D+00 <= b ) Then
       beta_log = w + gamma_log ( a ) + algdiv ( a, b )
       Return
    End If
    !
    !  Reduction of B when B < 8.
    !
60  Continue

    n = b - 1.0D+00
    z = 1.0D+00
    Do i = 1, n
       b = b - 1.0D+00
       z = z * ( b / ( a + b ))
    End Do

    beta_log = w + Log ( z ) + ( gamma_log ( a ) + ( gamma_log ( b ) &
         - gsumln ( a, b ) ) )
    Return
    !
    !  Reduction of A when 1000 < B.
    !
80  Continue

    n = a - 1.0D+00
    w = 1.0D+00
    Do i = 1, n
       a = a - 1.0D+00
       w = w * ( a / ( 1.0D+00 + a / b ))
    End Do

    beta_log = ( Log ( w ) - n * Log ( b ) ) + ( gamma_log ( a ) + algdiv ( a, b ) )

    Return
    !
    !  Procedure when 8 <= A.
    !
100 Continue

    w = bcorr ( a, b )
    h = a / b
    c = h / ( 1.0D+00 + h )
    u = - ( a - 0.5D+00 ) * Log ( c )
    v = b * alnrel ( h )

    If ( v < u ) Then
       beta_log = ((( -0.5D+00 * Log ( b ) + e ) + w ) - v ) - u
    Else
       beta_log = ((( -0.5D+00 * Log ( b ) + e ) + w ) - u ) - v
    End If

    Return
  End Function beta_log

  Function beta_pser ( a, b, x, eps )

    !*****************************************************************************80
    !
    !! BETA_PSER uses a power series expansion to evaluate IX(A,B)(X).
    !
    !  Discussion:
    !
    !    BETA_PSER is used when B <= 1 or B*X <= 0.7.
    !
    !  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, real ( kind = 8 ) A, B, the parameters.