ieva.f90

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

!*********************************************************
!*****   FDS5-Evac: Modules and misc routines        *****
!*********************************************************
!
! VTT Technical Research Centre of Finland            2007
!
! Author:  Timo Korhonen
! Date:    1.3.2007
! Version: 5.00
!
! This file (ieva.f90) contains:
! * dcdflib.f90 (Netlib cumulative density function library)
! * stat.f90 (random number generation etc, by T.K./Pfs3.0)
! * miscellaneous Evac routines, like initializations
!
!*********************************************************
Module DCDFLIB
  Implicit None
  CHARACTER(255), PARAMETER :: ievaid='$Id: ieva.f90 709 2007-09-28 17:47:43Z mcgratta $'
  CHARACTER(255), PARAMETER :: ievarev='$Revision: 709 $'
  CHARACTER(255), PARAMETER :: ievadate='$Date: 2007-09-28 13:47:43 -0400 (Fri, 28 Sep 2007) $'

  Private
  Public cdfbet,cdfgam,cdfnor,GET_REV_ieva

  !
  !*********************************************************
  !*****             DCDFLIB                           *****
  !*********************************************************
  !
  ! The Fortran90 version form:
  ! http://people.scs.fsu.edu/~burkardt/f_src/dcdflib/dcdflib.html
  !
  ! Original F77 version from Netlib: http://www.netlib.org/random:
  !
  !                                    DCDFLIB
  ! 
  !             Library of Fortran Routines for Cumulative Distribution
  !                  Functions, Inverses, and Other Parameters
  ! 
  !                                 (February, 1994)
  !                     Summary Documentation of Each Routine
  ! 
  !                             Compiled and Written by:
  ! 
  !                                  Barry W. Brown
  !                                   James Lovato
  !                                   Kathy Russell
  ! 
  !     Department of Biomathematics, Box 237
  !     The University of Texas, M.D. Anderson Cancer Center
  !     1515 Holcombe Boulevard
  !     Houston, TX      77030
  ! 
  !  This work was supported by grant CA-16672 from the National Cancer Institute.
  ! 
  ! 
  !                           SUMMARY OF DCDFLIB
  ! 
  ! This  library  contains routines  to compute  cumulative  distribution
  ! functions, inverses, and    parameters  of the  distribution  for  the
  ! following set of statistical distributions:
  ! 
  !     (1) Beta
  !     (2) Binomial
  !     (3) Chi-square
  !     (4) Noncentral Chi-square
  !     (5) F
  !     (6) Noncentral F
  !     (7) Gamma
  !     (8) Negative Binomial
  !     (9) Normal
  !     (10) Poisson
  !     (11) Student's t
  ! 
  ! Given values of all but one parameter of a distribution, the other is
  ! computed.  These calculations are done with  FORTRAN Double Precision 
  ! variables.
  !*********************************************************
  ! 
  !
  !*********************************************************

!!$  Real ( kind = 8 ) :: algdiv, alnrel, apser, bcorr, beta, beta_asym, &
!!$       beta_frac, beta_log, beta_pser, beta_rcomp, beta_rcomp1, &
!!$       beta_up, dbetrm, dexpm1, dinvnr, dlanor, dstrem, dt1, &
!!$       error_f, error_fc, esum, eval_pol, exparg, fpser, &
!!$       gam1, gamma, gamma_ln1, gamma_log, gsumln, ipmpar, &
!!$       psi, rcomp, rexp, rlog, rlog1, stvaln

Contains

  Function algdiv ( a, b )

    !*****************************************************************************80
    !
    !! ALGDIV computes ln ( Gamma ( B ) / Gamma ( A + B ) ) when 8 <= B.
    !
    !  Discussion:
    !
    !    In this algorithm, DEL(X) is the function defined by
    !
    !      ln ( Gamma(X) ) = ( X - 0.5 ) * ln ( X ) - X + 0.5 * ln ( 2 * PI ) 
    !                      + DEL(X).
    !
    !  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, define the arguments.
    !
    !    Output, real ( kind = 8 ) ALGDIV, the value of ln(Gamma(B)/Gamma(A+B)).
    !
    Implicit None

    Real ( kind = 8 ) a
    Real ( kind = 8 ) algdiv
    !    Real ( kind = 8 ) alnrel
    Real ( kind = 8 ) b
    Real ( kind = 8 ) c
    Real ( kind = 8 ), Parameter :: c0 =  0.833333333333333D-01
    Real ( kind = 8 ), Parameter :: c1 = -0.277777777760991D-02
    Real ( kind = 8 ), Parameter :: c2 =  0.793650666825390D-03
    Real ( kind = 8 ), Parameter :: c3 = -0.595202931351870D-03
    Real ( kind = 8 ), Parameter :: c4 =  0.837308034031215D-03
    Real ( kind = 8 ), Parameter :: c5 = -0.165322962780713D-02
    Real ( kind = 8 ) d
    Real ( kind = 8 ) h
    Real ( kind = 8 ) s11
    Real ( kind = 8 ) s3
    Real ( kind = 8 ) s5
    Real ( kind = 8 ) s7
    Real ( kind = 8 ) s9
    Real ( kind = 8 ) t
    Real ( kind = 8 ) u
    Real ( kind = 8 ) v
    Real ( kind = 8 ) w
    Real ( kind = 8 ) x
    Real ( kind = 8 ) x2

    If ( b < a ) Then
       h = b / a
       c = 1.0D+00 / ( 1.0D+00 + h )
       x = h / ( 1.0D+00 + h )
       d = a + ( b - 0.5D+00 )
    Else
       h = a / b
       c = h / ( 1.0D+00 + h )
       x = 1.0D+00 / ( 1.0D+00 + h )
       d = b + ( a - 0.5D+00 )
    End If
    !
    !  Set SN = (1 - X**N)/(1 - X).
    !
    x2 = x * x
    s3 = 1.0D+00 + ( x + x2 )
    s5 = 1.0D+00 + ( x + x2 * s3 )
    s7 = 1.0D+00 + ( x + x2 * s5 )
    s9 = 1.0D+00 + ( x + x2 * s7 )
    s11 = 1.0D+00 + ( x + x2 * s9 )
    !
    !  Set W = DEL(B) - DEL(A + B).
    !
    t = ( 1.0D+00 / b )**2
    w = (((( &
         c5 * s11  * t &
         + c4 * s9 ) * t &
         + c3 * s7 ) * t &
         + c2 * s5 ) * t &
         + c1 * s3 ) * t &
         + c0

    w = w * ( c / b )
    !
    !  Combine the results.
    !
    u = d * alnrel ( a / b )
    v = a * ( Log ( b ) - 1.0D+00 )

    If ( v < u ) Then
       algdiv = ( w - v ) - u
    Else
       algdiv = ( w - u ) - v
    End If

    Return
  End Function algdiv

  Function alnrel ( a )

    !*****************************************************************************80
    !
    !! ALNREL evaluates the function ln ( 1 + A ).
    !
    !  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, the argument.
    !
    !    Output, real ( kind = 8 ) ALNREL, the value of ln ( 1 + A ).
    !
    Implicit None

    Real ( kind = 8 ) a
    Real ( kind = 8 ) alnrel
    Real ( kind = 8 ), Parameter :: p1 = -0.129418923021993D+01
    Real ( kind = 8 ), Parameter :: p2 =  0.405303492862024D+00
    Real ( kind = 8 ), Parameter :: p3 = -0.178874546012214D-01
    Real ( kind = 8 ), Parameter :: q1 = -0.162752256355323D+01
    Real ( kind = 8 ), Parameter :: q2 =  0.747811014037616D+00
    Real ( kind = 8 ), Parameter :: q3 = -0.845104217945565D-01
    Real ( kind = 8 ) t
    Real ( kind = 8 ) t2
    Real ( kind = 8 ) w
    Real ( kind = 8 ) x

    If ( Abs ( a ) <= 0.375D+00 ) Then

       t = a / ( a +  2.0D+00  )
       t2 = t * t

       w = ((( p3 * t2 + p2 ) * t2 + p1 ) * t2 + 1.0D+00 ) &
            / ((( q3 * t2 + q2 ) * t2 + q1 ) * t2 + 1.0D+00 )

       alnrel =  2.0D+00  * t * w

    Else

       x = 1.0D+00 + Real ( a, kind = 8 )
       alnrel = Log ( x )

    End If

    Return
  End Function alnrel
  
  Function apser ( a, b, x, eps )

    !*****************************************************************************80
    !
    !! APSER computes the incomplete beta ratio I(SUB(1-X))(B,A).
    !
    !  Discussion:
    !
    !    APSER is used only for cases where
    !
    !      A <= min ( EPS, EPS * B ), 
    !      B * X <= 1, and 
    !      X <= 0.5.
    !
    !  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, X, the parameters of the
    !    incomplete beta ratio.
    !
    !    Input, real ( kind = 8 ) EPS, a tolerance.
    !
    !    Output, real ( kind = 8 ) APSER, the computed value of the
    !    incomplete beta ratio.
    !
    Implicit None

    Real ( kind = 8 ) a
    Real ( kind = 8 ) aj
    Real ( kind = 8 ) apser
    Real ( kind = 8 ) b
    Real ( kind = 8 ) bx
    Real ( kind = 8 ) c
    Real ( kind = 8 ) eps
    Real ( kind = 8 ), Parameter :: g = 0.577215664901533D+00
    Real ( kind = 8 ) j
    !    Real ( kind = 8 ) psi
    Real ( kind = 8 ) s
    Real ( kind = 8 ) t
    Real ( kind = 8 ) tol
    Real ( kind = 8 ) x

    bx = b * x
    t = x - bx

    If ( b * eps <= 0.02D+00 ) Then
       c = Log ( x ) + psi ( b ) + g + t
    Else
       c = Log ( bx ) + g + t
    End If

    tol = 5.0D+00 * eps * Abs ( c )
    j = 1.0D+00
    s = 0.0D+00

    Do

       j = j + 1.0D+00
       t = t * ( x - bx / j )
       aj = t / j
       s = s + aj

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

    End Do

    apser = -a * ( c + s )

    Return
  End Function apser

  Function bcorr ( a0, b0 )

    !*****************************************************************************80
    !
    !! BCORR evaluates DEL(A0) + DEL(B0) - DEL(A0 + B0).
    !
    !  Discussion:
    !
    !    The function DEL(A) is a remainder term that is used in the expression:
    !
    !      ln ( Gamma ( A ) ) = ( A - 0.5 ) * ln ( A ) 
    !        - A + 0.5 * ln ( 2 * PI ) + DEL ( A ),
    !
    !    or, in other words, DEL ( A ) is defined as:
    !
    !      DEL ( A ) = ln ( Gamma ( A ) ) - ( A - 0.5 ) * ln ( A ) 
    !        + A + 0.5 * ln ( 2 * PI ).
    !
    !  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 arguments.
    !    It is assumed that 8 <= A0 and 8 <= B0.
    !
    !    Output, real ( kind = 8 ) BCORR, the value of the function.
    !
    Implicit None

    Real ( kind = 8 ) a
    Real ( kind = 8 ) a0
    Real ( kind = 8 ) b
    Real ( kind = 8 ) b0
    Real ( kind = 8 ) bcorr
    Real ( kind = 8 ) c
    Real ( kind = 8 ), Parameter :: c0 =  0.833333333333333D-01
    Real ( kind = 8 ), Parameter :: c1 = -0.277777777760991D-02
    Real ( kind = 8 ), Parameter :: c2 =  0.793650666825390D-03
    Real ( kind = 8 ), Parameter :: c3 = -0.595202931351870D-03
    Real ( kind = 8 ), Parameter :: c4 =  0.837308034031215D-03
    Real ( kind = 8 ), Parameter :: c5 = -0.165322962780713D-02
    Real ( kind = 8 ) h
    Real ( kind = 8 ) s11
    Real ( kind = 8 ) s3
    Real ( kind = 8 ) s5
    Real ( kind = 8 ) s7
    Real ( kind = 8 ) s9
    Real ( kind = 8 ) t
    Real ( kind = 8 ) w
    Real ( kind = 8 ) x
    Real ( kind = 8 ) x2

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

    h = a / b
    c = h / ( 1.0D+00 + h )
    x = 1.0D+00 / ( 1.0D+00 + h )
    x2 = x * x
    !
    !  Set SN = (1 - X**N)/(1 - X)
    !
    s3 = 1.0D+00 + ( x + x2 )
    s5 = 1.0D+00 + ( x + x2 * s3 )
    s7 = 1.0D+00 + ( x + x2 * s5 )
    s9 = 1.0D+00 + ( x + x2 * s7 )
    s11 = 1.0D+00 + ( x + x2 * s9 )
    !
    !  Set W = DEL(B) - DEL(A + B)
    !
    t = ( 1.0D+00 / b )**2

    w = (((( &
         c5 * s11  * t &
         + c4 * s9 ) * t &
         + c3 * s7 ) * t &
         + c2 * s5 ) * t &
         + c1 * s3 ) * t &
         + c0

    w = w * ( c / b )
    !
    !  Compute  DEL(A) + W.
    !
    t = ( 1.0D+00 / a )**2

    bcorr = ((((( &
         c5   * t &
         + c4 ) * t &
         + c3 ) * t &
         + c2 ) * t &
         + c1 ) * t &
         + c0 ) / a + w

    Return
  End Function bcorr

  Function beta ( a, b )

    !*****************************************************************************80
    !
    !! BETA evaluates the beta function.
    !
    !  Modified:
    !
    !    03 December 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, real ( kind = 8 ) A, B, the arguments of the beta function.
    !
    !    Output, real ( kind = 8 ) BETA, the value of the beta function.
    !
    Implicit None

    Real ( kind = 8 ) a
    Real ( kind = 8 ) b
    Real ( kind = 8 ) beta
    !    Real ( kind = 8 ) beta_log

    beta = Exp ( beta_log ( a, b ) )

    Return
  End Function beta

  Function beta_asym ( a, b, lambda, eps )

    !*****************************************************************************80
    !
    !! BETA_ASYM computes an asymptotic expansion for IX(A,B), for large A and 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 or equal to 15.
    !
    !    Input, real ( kind = 8 ) LAMBDA, the value of ( A + B ) * Y - B.
    !    It is assumed that 0 <= LAMBDA.
    !
    !    Input, real ( kind = 8 ) EPS, the tolerance.
    !
    Implicit None

    Integer, Parameter :: num = 20

    Real ( kind = 8 ) a
    Real ( kind = 8 ) a0(num+1)
    Real ( kind = 8 ) b
    Real ( kind = 8 ) b0(num+1)
    !    Real ( kind = 8 ) bcorr
    Real ( kind = 8 ) beta_asym
    Real ( kind = 8 ) bsum
    Real ( kind = 8 ) c(num+1)
    Real ( kind = 8 ) d(num+1)
    Real ( kind = 8 ) dsum
    Real ( kind = 8 ), Parameter :: e0 = 1.12837916709551D+00
    Real ( kind = 8 ), Parameter :: e1 = 0.353553390593274D+00
    Real ( kind = 8 ) eps
    !    Real ( kind = 8 ) error_fc
    Real ( kind = 8 ) f
    Real ( kind = 8 ) h
    Real ( kind = 8 ) h2
    Real ( kind = 8 ) hn
    Integer i
    Integer j
    Real ( kind = 8 ) j0
    Real ( kind = 8 ) j1
    Real ( kind = 8 ) lambda
    Integer m
    Integer mm1
    Integer mmj
    Integer n
    Integer np1
    Real ( kind = 8 ) r
    Real ( kind = 8 ) r0
    Real ( kind = 8 ) r1
    !    Real ( kind = 8 ) rlog1