ieva.f90

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

    Real ( kind = 8 ) s
    Real ( kind = 8 ) sum1
    Real ( kind = 8 ) t
    Real ( kind = 8 ) t0
    Real ( kind = 8 ) t1
    Real ( kind = 8 ) u
    Real ( kind = 8 ) w
    Real ( kind = 8 ) w0
    Real ( kind = 8 ) z
    Real ( kind = 8 ) z0
    Real ( kind = 8 ) z2
    Real ( kind = 8 ) zn
    Real ( kind = 8 ) znm1

    beta_asym = 0.0D+00

    If ( a < b ) Then
       h = a / b
       r0 = 1.0D+00 / ( 1.0D+00 + h )
       r1 = ( b - a ) / b
       w0 = 1.0D+00 / Sqrt ( a * ( 1.0D+00 + h ))
    Else
       h = b / a
       r0 = 1.0D+00 / ( 1.0D+00 + h )
       r1 = ( b - a ) / a
       w0 = 1.0D+00 / Sqrt ( b * ( 1.0D+00 + h ))
    End If

    f = a * rlog1 ( - lambda / a ) + b * rlog1 ( lambda / b )
    t = Exp ( - f )
    If ( t == 0.0D+00 ) Then
       Return
    End If

    z0 = Sqrt ( f )
    z = 0.5D+00 * ( z0 / e1 )
    z2 = f + f

    a0(1) = ( 2.0D+00 / 3.0D+00 ) * r1
    c(1) = -0.5D+00 * a0(1)
    d(1) = -c(1)
    j0 = ( 0.5D+00 / e0 ) * error_fc ( 1, z0 )
    j1 = e1
    sum1 = j0 + d(1) * w0 * j1

    s = 1.0D+00
    h2 = h * h
    hn = 1.0D+00
    w = w0
    znm1 = z
    zn = z2

    Do n = 2, num, 2

       hn = h2 * hn
       a0(n) = 2.0D+00 * r0 * ( 1.0D+00 + h * hn ) &
            / ( n + 2.0D+00 )
       np1 = n + 1
       s = s + hn
       a0(np1) = 2.0D+00 * r1 * s / ( n + 3.0D+00 )

       Do i = n, np1

          r = -0.5D+00 * ( i + 1.0D+00 )
          b0(1) = r * a0(1)
          Do m = 2, i
             bsum = 0.0D+00
             mm1 = m - 1
             Do j = 1, mm1
                mmj = m - j
                bsum = bsum + ( j * r - mmj ) * a0(j) * b0(mmj)
             End Do
             b0(m) = r * a0(m) + bsum / m
          End Do

          c(i) = b0(i) / ( i + 1.0D+00 )

          dsum = 0.0
          Do j = 1, i-1
             dsum = dsum + d(i-j) * c(j)
          End Do
          d(i) = - ( dsum + c(i) )

       End Do

       j0 = e1 * znm1 + ( n - 1.0D+00 ) * j0
       j1 = e1 * zn + n * j1
       znm1 = z2 * znm1
       zn = z2 * zn
       w = w0 * w
       t0 = d(n) * w * j0
       w = w0 * w
       t1 = d(np1) * w * j1
       sum1 = sum1 + ( t0 + t1 )

       If ( ( Abs ( t0 ) + Abs ( t1 )) <= eps * sum1 ) Then
          u = Exp ( - bcorr ( a, b ) )
          beta_asym = e0 * t * u * sum1
          Return
       End If

    End Do

    u = Exp ( - bcorr ( a, b ) )
    beta_asym = e0 * t * u * sum1

    Return
  End Function beta_asym

  Function beta_frac ( a, b, x, y, lambda, eps )

    !*****************************************************************************80
    !
    !! BETA_FRAC evaluates a continued fraction expansion for IX(A,B).
    !
    !  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 of the function.
    !    A and B should be nonnegative.  It is assumed that both A and
    !    B are greater than 1.
    !
    !    Input, real ( kind = 8 ) X, Y.  X is the argument of the
    !    function, and should satisy 0 <= X <= 1.  Y should equal 1 - X.
    !
    !    Input, real ( kind = 8 ) LAMBDA, the value of ( A + B ) * Y - B.
    !
    !    Input, real ( kind = 8 ) EPS, a tolerance.
    !
    !    Output, real ( kind = 8 ) BETA_FRAC, the value of the continued
    !    fraction approximation for IX(A,B).
    !
    Implicit None

    Real ( kind = 8 ) a
    Real ( kind = 8 ) alpha
    Real ( kind = 8 ) an
    Real ( kind = 8 ) anp1
    Real ( kind = 8 ) b
    Real ( kind = 8 ) beta
    Real ( kind = 8 ) beta_frac
    !    Real ( kind = 8 ) beta_rcomp
    Real ( kind = 8 ) bn
    Real ( kind = 8 ) bnp1
    Real ( kind = 8 ) c
    Real ( kind = 8 ) c0
    Real ( kind = 8 ) c1
    Real ( kind = 8 ) e
    Real ( kind = 8 ) eps
    Real ( kind = 8 ) lambda
    Real ( kind = 8 ) n
    Real ( kind = 8 ) p
    Real ( kind = 8 ) r
    Real ( kind = 8 ) r0
    Real ( kind = 8 ) s
    Real ( kind = 8 ) t
    Real ( kind = 8 ) w
    Real ( kind = 8 ) x
    Real ( kind = 8 ) y
    Real ( kind = 8 ) yp1

    beta_frac = beta_rcomp ( a, b, x, y )

    If ( beta_frac == 0.0D+00 ) Then
       Return
    End If

    c = 1.0D+00 + lambda
    c0 = b / a
    c1 = 1.0D+00 + 1.0D+00 / a
    yp1 = y + 1.0D+00

    n = 0.0D+00
    p = 1.0D+00
    s = a + 1.0D+00
    an = 0.0D+00
    bn = 1.0D+00
    anp1 = 1.0D+00
    bnp1 = c / c1
    r = c1 / c
    !
    !  Continued fraction calculation.
    !
    Do

       n = n + 1.0D+00
       t = n / a
       w = n * ( b - n ) * x
       e = a / s
       alpha = ( p * ( p + c0 ) * e * e ) * ( w * x )
       e = ( 1.0D+00 + t ) / ( c1 + t + t )
       beta = n + w / s + e * ( c + n * yp1 )
       p = 1.0D+00 + t
       s = s +  2.0D+00 
       !
       !  Update AN, BN, ANP1, and BNP1.
       !
       t = alpha * an + beta * anp1
       an = anp1
       anp1 = t
       t = alpha * bn + beta * bnp1
       bn = bnp1
       bnp1 = t

       r0 = r
       r = anp1 / bnp1

       If ( Abs ( r - r0 ) <= eps * r ) Then
          beta_frac = beta_frac * r
          Exit
       End If
       !
       !  Rescale AN, BN, ANP1, and BNP1.
       !
       an = an / bnp1
       bn = bn / bnp1
       anp1 = r
       bnp1 = 1.0D+00

    End Do

    Return
  End Function beta_frac

  Subroutine beta_grat ( a, b, x, y, w, eps, ierr )

    !*****************************************************************************80
    !
    !! BETA_GRAT evaluates an asymptotic expansion for IX(A,B).
    !
    !  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 of the function.
    !    A and B should be nonnegative.  It is assumed that 15 <= A 
    !    and B <= 1, and that B is less than A.
    !
    !    Input, real ( kind = 8 ) X, Y.  X is the argument of the
    !    function, and should satisy 0 <= X <= 1.  Y should equal 1 - X.
    !
    !    Input/output, real ( kind = 8 ) W, a quantity to which the
    !    result of the computation is to be added on output.
    !
    !    Input, real ( kind = 8 ) EPS, a tolerance.
    !
    !    Output, integer IERR, an error flag, which is 0 if no error
    !    was detected.
    !
    Implicit None

    Real ( kind = 8 ) a
    !    Real ( kind = 8 ) algdiv
    !    Real ( kind = 8 ) alnrel
    Real ( kind = 8 ) b
    Real ( kind = 8 ) bm1
    Real ( kind = 8 ) bp2n
    Real ( kind = 8 ) c(30)
    Real ( kind = 8 ) cn
    Real ( kind = 8 ) coef
    Real ( kind = 8 ) d(30)
    Real ( kind = 8 ) dj
    Real ( kind = 8 ) eps
    !    Real ( kind = 8 ) gam1
    Integer i
    Integer ierr
    Real ( kind = 8 ) j
    Real ( kind = 8 ) l
    Real ( kind = 8 ) lnx
    Integer n
    Real ( kind = 8 ) n2
    Real ( kind = 8 ) nu
    Real ( kind = 8 ) p
    Real ( kind = 8 ) q
    Real ( kind = 8 ) r
    Real ( kind = 8 ) s
    Real ( kind = 8 ) sum1
    Real ( kind = 8 ) t
    Real ( kind = 8 ) t2
    Real ( kind = 8 ) u
    Real ( kind = 8 ) v
    Real ( kind = 8 ) w
    Real ( kind = 8 ) x
    Real ( kind = 8 ) y
    Real ( kind = 8 ) z

    bm1 = ( b - 0.5D+00 ) - 0.5D+00
    nu = a + 0.5D+00 * bm1

    If ( y <= 0.375D+00 ) Then
       lnx = alnrel ( - y )
    Else
       lnx = Log ( x )
    End If

    z = -nu * lnx

    If ( b * z == 0.0D+00 ) Then
       ierr = 1
       Return
    End If
    !
    !  Computation of the expansion.
    !
    !  Set R = EXP(-Z)*Z**B/GAMMA(B)
    !
    r = b * ( 1.0D+00 + gam1 ( b ) ) * Exp ( b * Log ( z ))
    r = r * Exp ( a * lnx ) * Exp ( 0.5D+00 * bm1 * lnx )
    u = algdiv ( b, a ) + b * Log ( nu )
    u = r * Exp ( - u )

    If ( u == 0.0D+00 ) Then
       ierr = 1
       Return
    End If

    Call gamma_rat1 ( b, z, r, p, q, eps )

    v = 0.25D+00 * ( 1.0D+00 / nu )**2
    t2 = 0.25D+00 * lnx * lnx
    l = w / u
    j = q / r
    sum1 = j
    t = 1.0D+00
    cn = 1.0D+00
    n2 = 0.0D+00

    Do n = 1, 30

       bp2n = b + n2
       j = ( bp2n * ( bp2n + 1.0D+00 ) * j &
            + ( z + bp2n + 1.0D+00 ) * t ) * v
       n2 = n2 +  2.0D+00 
       t = t * t2
       cn = cn / ( n2 * ( n2 + 1.0D+00 ))
       c(n) = cn
       s = 0.0D+00

       coef = b - n
       Do i = 1, n-1
          s = s + coef * c(i) * d(n-i)
          coef = coef + b
       End Do

       d(n) = bm1 * cn + s / n
       dj = d(n) * j
       sum1 = sum1 + dj

       If ( sum1 <= 0.0D+00 ) Then
          ierr = 1
          Return
       End If

       If ( Abs ( dj ) <= eps * ( sum1 + l ) ) Then
          ierr = 0
          w = w + u * sum1
          Return
       End If

    End Do

    ierr = 0
    w = w + u * sum1

    Return
  End Subroutine beta_grat

  Subroutine beta_inc ( a, b, x, y, w, w1, ierr )

    !*****************************************************************************80
    !
    !! BETA_INC evaluates the incomplete beta function IX(A,B).
    !
    !  Author:
    !
    !    Alfred Morris,
    !    Naval Surface Weapons Center,
    !    Dahlgren, Virginia.
    !
    !  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 of the function.
    !    A and B should be nonnegative.
    !
    !    Input, real ( kind = 8 ) X, Y.  X is the argument of the
    !    function, and should satisy 0 <= X <= 1.  Y should equal 1 - X.
    !
    !    Output, real ( kind = 8 ) W, W1, the values of IX(A,B) and
    !    1-IX(A,B).
    !
    !    Output, integer IERR, the error flag.
    !    0, no error was detected.
    !    1, A or B is negative;
    !    2, A = B = 0;
    !    3, X < 0 or 1 < X;
    !    4, Y < 0 or 1 < Y;
    !    5, X + Y /= 1;
    !    6, X = A = 0;
    !    7, Y = B = 0.
    !
    Implicit None

    Real ( kind = 8 ) a
    Real ( kind = 8 ) a0
    !    Real ( kind = 8 ) apser
    Real ( kind = 8 ) b
    Real ( kind = 8 ) b0
    !    Real ( kind = 8 ) beta_asym
    !    Real ( kind = 8 ) beta_frac
    !    Real ( kind = 8 ) beta_pser
    !    Real ( kind = 8 ) beta_up
    Real ( kind = 8 ) eps
    !    Real ( kind = 8 ) fpser
    Integer ierr
    Integer ierr1
    Integer ind
    Real ( kind = 8 ) lambda
    Integer n
    Real ( kind = 8 ) t
    Real ( kind = 8 ) w
    Real ( kind = 8 ) w1
    Real ( kind = 8 ) x
    Real ( kind = 8 ) x0
    Real ( kind = 8 ) y
    Real ( kind = 8 ) y0
    Real ( kind = 8 ) z

    eps = Epsilon ( eps )
    w = 0.0D+00
    w1 = 0.0D+00

    If ( a < 0.0D+00 .Or. b < 0.0D+00 ) Then
       ierr = 1
       Return
    End If

    If ( a == 0.0D+00 .And. b == 0.0D+00 ) Then
       ierr = 2
       Return
    End If

    If ( x < 0.0D+00 .Or. 1.0D+00 < x ) Then
       ierr = 3
       Return
    End If

    If ( y < 0.0D+00 .Or. 1.0D+00 < y ) Then
       ierr = 4
       Return
    End If

    z = ( ( x + y ) - 0.5D+00 ) - 0.5D+00

    If ( 3.0D+00 * eps < Abs ( z ) ) Then
       ierr = 5
       Return
    End If

    ierr = 0

    If ( x == 0.0D+00 ) Then
       w = 0.0D+00
       w1 = 1.0D+00
       If ( a == 0.0D+00 ) Then
          ierr = 6
       End If
       Return
    End If

    If ( y == 0.0D+00 ) Then
       If ( b == 0.0D+00 ) Then
          ierr = 7
          Return
       End If
       w = 1.0D+00
       w1 = 0.0D+00
       Return
    End If

    If ( a == 0.0D+00 ) Then
       w = 1.0D+00
       w1 = 0.0D+00
       Return
    End If

    If ( b == 0.0D+00 ) Then
       w = 0.0D+00
       w1 = 1.0D+00
       Return
    End If

    eps = Max ( eps, 1.0D-15 )

    If ( Max ( a, b ) < 0.001D+00 * eps ) Then
       go to 260
    End If

    ind = 0
    a0 = a
    b0 = b
    x0 = x
    y0 = y

    If ( 1.0D+00 < Min ( a0, b0 ) ) Then
       go to 40