ieva.f90

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

    !
    !    Input, real ( kind = 8 ) X, the point where the function
    !    is to be evaluated.
    !
    !    Input, real ( kind = 8 ) EPS, the tolerance.
    !
    !    Output, real ( kind = 8 ) BETA_PSER, the approximate value of IX(A,B)(X).
    !
    Implicit None

    Real ( kind = 8 ) a
    Real ( kind = 8 ) a0
    !    Real ( kind = 8 ) algdiv
    Real ( kind = 8 ) apb
    Real ( kind = 8 ) b
    Real ( kind = 8 ) b0
    !    Real ( kind = 8 ) beta_log
    Real ( kind = 8 ) beta_pser
    Real ( kind = 8 ) c
    Real ( kind = 8 ) eps
    !    Real ( kind = 8 ) gam1
    !    Real ( kind = 8 ) gamma_ln1
    Integer i
    Integer m
    Real ( kind = 8 ) n
    Real ( kind = 8 ) sum1
    Real ( kind = 8 ) t
    Real ( kind = 8 ) tol
    Real ( kind = 8 ) u
    Real ( kind = 8 ) w
    Real ( kind = 8 ) x
    Real ( kind = 8 ) z

    beta_pser = 0.0D+00

    If ( x == 0.0D+00 ) Then
       Return
    End If
    !
    !  Compute the factor X**A/(A*BETA(A,B))
    !
    a0 = Min ( a, b )

    If ( 1.0D+00 <= a0 ) Then

       z = a * Log ( x ) - beta_log ( a, b )
       beta_pser = Exp ( z ) / a

    Else

       b0 = Max ( a, b )

       If ( b0 <= 1.0D+00 ) Then

          beta_pser = x**a
          If ( beta_pser == 0.0D+00 ) Then
             Return
          End If

          apb = a + b

          If ( apb <= 1.0D+00 ) Then
             z = 1.0D+00 + gam1 ( apb )
          Else
             u = a + b - 1.0D+00
             z = ( 1.0D+00 + gam1 ( u ) ) / apb
          End If

          c = ( 1.0D+00 + gam1 ( a ) ) &
               * ( 1.0D+00 + gam1 ( b ) ) / z
          beta_pser = beta_pser * c * ( b / apb )

       Else If ( b0 < 8.0D+00 ) Then

          u = gamma_ln1 ( a0 )
          m = b0 - 1.0D+00

          c = 1.0D+00
          Do i = 1, m
             b0 = b0 - 1.0D+00
             c = c * ( b0 / ( a0 + b0 ))
          End Do

          u = Log ( c ) + u
          z = a * Log ( x ) - u
          b0 = b0 - 1.0D+00
          apb = a0 + b0

          If ( apb <= 1.0D+00 ) Then
             t = 1.0D+00 + gam1 ( apb )
          Else
             u = a0 + b0 - 1.0D+00
             t = ( 1.0D+00 + gam1 ( u ) ) / apb
          End If

          beta_pser = Exp ( z ) * ( a0 / a ) &
               * ( 1.0D+00 + gam1 ( b0 )) / t

       Else If ( 8.0D+00 <= b0 ) Then

          u = gamma_ln1 ( a0 ) + algdiv ( a0, b0 )
          z = a * Log ( x ) - u
          beta_pser = ( a0 / a ) * Exp ( z )

       End If

    End If

    If ( beta_pser == 0.0D+00 .Or. a <= 0.1D+00 * eps ) Then
       Return
    End If
    !
    !  Compute the series.
    !
    sum1 = 0.0D+00
    n = 0.0D+00
    c = 1.0D+00
    tol = eps / a

    Do

       n = n + 1.0D+00
       c = c * ( 0.5D+00 + ( 0.5D+00 - b / n ) ) * x
       w = c / ( a + n )
       sum1 = sum1 + w

       If ( Abs ( w ) <= tol ) Then
          Exit
       End If

    End Do

    beta_pser = beta_pser * ( 1.0D+00 + a * sum1 )

    Return
  End Function beta_pser

  Function beta_rcomp ( a, b, x, y )

    !*****************************************************************************80
    !
    !! BETA_RCOMP evaluates X**A * Y**B / Beta(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 Beta function.
    !    A and B should be nonnegative.
    !
    !    Input, real ( kind = 8 ) X, Y, define the numerator of the fraction.
    !
    !    Output, real ( kind = 8 ) BETA_RCOMP, the value of X**A * Y**B / Beta(A,B).
    !
    Implicit None

    Real ( kind = 8 ) a
    Real ( kind = 8 ) a0
    !    Real ( kind = 8 ) algdiv
    !    Real ( kind = 8 ) alnrel
    Real ( kind = 8 ) apb
    Real ( kind = 8 ) b
    Real ( kind = 8 ) b0
    !    Real ( kind = 8 ) bcorr
    !    Real ( kind = 8 ) beta_log
    Real ( kind = 8 ) beta_rcomp
    Real ( kind = 8 ) c
    Real ( kind = 8 ), Parameter :: const = 0.398942280401433D+00
    Real ( kind = 8 ) e
    !    Real ( kind = 8 ) gam1
    !    Real ( kind = 8 ) gamma_ln1
    Real ( kind = 8 ) h
    Integer i
    Real ( kind = 8 ) lambda
    Real ( kind = 8 ) lnx
    Real ( kind = 8 ) lny
    Integer n
    !    Real ( kind = 8 ) rlog1
    Real ( kind = 8 ) t
    Real ( kind = 8 ) u
    Real ( kind = 8 ) v
    Real ( kind = 8 ) x
    Real ( kind = 8 ) x0
    Real ( kind = 8 ) y
    Real ( kind = 8 ) y0
    Real ( kind = 8 ) z

    beta_rcomp = 0.0D+00
    If ( x == 0.0D+00 .Or. y == 0.0D+00 ) Then
       Return
    End If

    a0 = Min ( a, b )

    If ( a0 < 8.0D+00 ) Then

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

       z = a * lnx + b * lny

       If ( 1.0D+00 <= a0 ) Then
          z = z - beta_log ( a, b )
          beta_rcomp = Exp ( z )
          Return
       End If
       !
       !  Procedure for A < 1 or B < 1
       !
       b0 = Max ( a, b )

       If ( b0 <= 1.0D+00 ) Then

          beta_rcomp = Exp ( z )
          If ( beta_rcomp == 0.0D+00 ) Then
             Return
          End If

          apb = a + b

          If ( apb <= 1.0D+00 ) Then
             z = 1.0D+00 + gam1 ( apb )
          Else
             u = a + b - 1.0D+00
             z = ( 1.0D+00 + gam1 ( u ) ) / apb
          End If

          c = ( 1.0D+00 + gam1 ( a ) ) &
               * ( 1.0D+00 + gam1 ( b ) ) / z
          beta_rcomp = beta_rcomp * ( a0 * c ) &
               / ( 1.0D+00 + a0 / b0 )

       Else If ( b0 < 8.0D+00 ) Then

          u = gamma_ln1 ( a0 )
          n = b0 - 1.0D+00

          c = 1.0D+00
          Do i = 1, n
             b0 = b0 - 1.0D+00
             c = c * ( b0 / ( a0 + b0 ))
          End Do
          u = Log ( c ) + u

          z = z - u
          b0 = b0 - 1.0D+00
          apb = a0 + b0

          If ( apb <= 1.0D+00 ) Then
             t = 1.0D+00 + gam1 ( apb )
          Else
             u = a0 + b0 - 1.0D+00
             t = ( 1.0D+00 + gam1 ( u ) ) / apb
          End If

          beta_rcomp = a0 * Exp ( z ) * ( 1.0D+00 + gam1 ( b0 ) ) / t

       Else If ( 8.0D+00 <= b0 ) Then

          u = gamma_ln1 ( a0 ) + algdiv ( a0, b0 )
          beta_rcomp = a0 * Exp ( z - u ) 

       End If

    Else

       If ( a <= b ) Then
          h = a / b
          x0 = h / ( 1.0D+00 + h )
          y0 = 1.0D+00 / (  1.0D+00 + h )
          lambda = a - ( a + b ) * x
       Else
          h = b / a
          x0 = 1.0D+00 / ( 1.0D+00 + h )
          y0 = h / ( 1.0D+00 + h )
          lambda = ( a + b ) * y - b
       End If

       e = -lambda / a

       If ( Abs ( e ) <= 0.6D+00 ) Then
          u = rlog1 ( e )
       Else
          u = e - Log ( x / x0 )
       End If

       e = lambda / b

       If ( Abs ( e ) <= 0.6D+00 ) Then
          v = rlog1 ( e )
       Else
          v = e - Log ( y / y0 )
       End If

       z = Exp ( - ( a * u + b * v ) )
       beta_rcomp = const * Sqrt ( b * x0 ) * z * Exp ( - bcorr ( a, b ))

    End If

    Return
  End Function beta_rcomp

  Function beta_rcomp1 ( mu, a, b, x, y )

    !*****************************************************************************80
    !
    !! BETA_RCOMP1 evaluates exp(MU) * X**A * Y**B / Beta(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, integer MU, ?
    !
    !    Input, real ( kind = 8 ) A, B, the parameters of the Beta function.
    !    A and B should be nonnegative.
    !
    !    Input, real ( kind = 8 ) X, Y, ?
    !
    !    Output, real ( kind = 8 ) BETA_RCOMP1, the value of
    !    exp(MU) * X**A * Y**B / Beta(A,B).
    !
    Implicit None

    Real ( kind = 8 ) a
    Real ( kind = 8 ) a0
    !    Real ( kind = 8 ) algdiv
    !    Real ( kind = 8 ) alnrel
    Real ( kind = 8 ) apb
    Real ( kind = 8 ) b
    Real ( kind = 8 ) b0
    !    Real ( kind = 8 ) bcorr
    !    Real ( kind = 8 ) beta_log
    Real ( kind = 8 ) beta_rcomp1
    Real ( kind = 8 ) c
    Real ( kind = 8 ), Parameter :: const = 0.398942280401433D+00
    Real ( kind = 8 ) e
    !    Real ( kind = 8 ) esum
    !    Real ( kind = 8 ) gam1
    !    Real ( kind = 8 ) gamma_ln1
    Real ( kind = 8 ) h
    Integer i
    Real ( kind = 8 ) lambda
    Real ( kind = 8 ) lnx
    Real ( kind = 8 ) lny
    Integer mu
    Integer n
    !    Real ( kind = 8 ) rlog1
    Real ( kind = 8 ) t
    Real ( kind = 8 ) u
    Real ( kind = 8 ) v
    Real ( kind = 8 ) x
    Real ( kind = 8 ) x0
    Real ( kind = 8 ) y
    Real ( kind = 8 ) y0
    Real ( kind = 8 ) z

    a0 = Min ( a, b )

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

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

    z = a * lnx + b * lny

    If ( 1.0D+00 <= a0 ) Then
       z = z - beta_log ( a, b )
       beta_rcomp1 = esum ( mu, z )
       Return
    End If
    !
    !  Procedure for A < 1 OR B < 1
    !
40  Continue

    b0 = Max ( a, b )

    If ( 8.0D+00 <= b0 ) Then
       go to 120
    End If

    If ( 1.0D+00 < b0 ) Then
       go to 70
    End If
    !
    !  Algorithm for B0 <= 1
    !
    beta_rcomp1 = esum ( mu, z )
    If ( beta_rcomp1 == 0.0D+00 ) Then
       Return
    End If

    apb = a + b

    If ( apb <= 1.0D+00 ) Then
       z = 1.0D+00 + gam1 ( apb )
    Else
       u = Real ( a, kind = 8 ) + Real ( b, kind = 8 ) - 1.0D+00
       z = ( 1.0D+00 + gam1 ( u )) / apb
    End If

    c = ( 1.0D+00 + gam1 ( a ) ) &
         * ( 1.0D+00 + gam1 ( b ) ) / z
    beta_rcomp1 = beta_rcomp1 * ( a0 * c ) / ( 1.0D+00 + a0 / b0 )
    Return
    !
    !  Algorithm for 1 < B0 < 8
    !
70  Continue

    u = gamma_ln1 ( a0 )
    n = b0 - 1.0D+00

    c = 1.0D+00
    Do i = 1, n
       b0 = b0 - 1.0D+00
       c = c * ( b0 / ( a0 + b0 ) )
    End Do
    u = Log ( c ) + u

    z = z - u
    b0 = b0 - 1.0D+00
    apb = a0 + b0

    If ( apb <= 1.0D+00 ) Then
       t = 1.0D+00 + gam1 ( apb )
    Else
       u = a0 + b0 - 1.0D+00
       t = ( 1.0D+00 + gam1 ( u ) ) / apb
    End If

    beta_rcomp1 = a0 * esum ( mu, z ) &
         * ( 1.0D+00 + gam1 ( b0 ) ) / t
    Return
    !
    !  Algorithm for 8 <= B0.
    !
120 Continue

    u = gamma_ln1 ( a0 ) + algdiv ( a0, b0 )
    beta_rcomp1 = a0 * esum ( mu, z-u )
    Return
    !
    !  Procedure for 8 <= A and 8 <= B.
    !
130 Continue

    If ( a <= b ) Then
       h = a / b
       x0 = h / ( 1.0D+00 + h )
       y0 = 1.0D+00 / ( 1.0D+00 + h )
       lambda = a - ( a + b ) * x
    Else
       h = b / a
       x0 = 1.0D+00 / ( 1.0D+00 + h )
       y0 = h / ( 1.0D+00 + h )
       lambda = ( a + b ) * y - b
    End If

    e = -lambda / a

    If ( Abs ( e ) <= 0.6D+00 ) Then
       u = rlog1 ( e )
    Else
       u = e - Log ( x / x0 )
    End If

    e = lambda / b

    If ( Abs ( e ) <= 0.6D+00 ) Then
       v = rlog1 ( e )
    Else
       v = e - Log ( y / y0 )
    End If

    z = esum ( mu, - ( a * u + b * v ))
    beta_rcomp1 = const * Sqrt ( b * x0 ) * z * Exp ( - bcorr ( a, b ) )

    Return
  End Function beta_rcomp1

  Function beta_up ( a, b, x, y, n, eps )

    !*****************************************************************************80
    !
    !! BETA_UP evaluates IX(A,B) - IX(A+N,B) where N is a positive integer.
    !
    !  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.
    !