IntGamma Function

  1. korc_experimental_pdf.f90
  2. korc_experimental_pdf
  3. IntGamma

private function IntGamma(a, b, k, t)

@brief Extended trapezoidal rule for integrating the Gamma PDF. See Sec. 4.2 of Numerical Recipies in Fortran 77.

Arguments

Type IntentOptional AttributesName
real(kind=rp), intent(in) :: a
real(kind=rp), intent(in) :: b
real(kind=rp), intent(in) :: k
real(kind=rp), intent(in) :: t

Return Value real(kind=rp)


Calls

proc~~intgamma~~CallsGraph proc~intgamma IntGamma proc~fgamma~2 fGamma proc~intgamma->proc~fgamma~2

Called by

proc~~intgamma~~CalledByGraph proc~intgamma IntGamma proc~initialize_params~2 initialize_params proc~initialize_params~2->proc~intgamma proc~get_experimentalg_distribution get_experimentalG_distribution proc~get_experimentalg_distribution->proc~initialize_params~2

Contents

Source Code

IntGamma

Source Code

  FUNCTION IntGamma(a,b,k,t)
    REAL(rp), INTENT(IN) :: a
    REAL(rp), INTENT(IN) :: b
    REAL(rp), INTENT(IN) :: k ! shape factor
    REAL(rp), INTENT(IN) :: t ! scale factor
    REAL(rp) :: IntGamma
    REAL(rp) :: Iold
    REAL(rp) :: Inew
    REAL(rp) :: rerr
    REAL(rp) :: sum_f
    REAL(rp) :: h,z
    INTEGER :: ii,jj,npoints
    LOGICAL :: flag

    h = b - a
    sum_f = 0.5*(fGamma(a,k,t) + fGamma(b,k,t))

    Iold = 0.0_rp
    Inew = sum_f*h

    ii = 1_idef
    flag = .TRUE.
    do while (flag)
       Iold = Inew

       ii = ii + 1_idef
       npoints = 2_idef**(ii-2_idef)
       h = 0.5_rp*(b-a)/REAL(npoints,rp)
       sum_f = 0.0_rp
       do jj=1_idef,npoints
          z = a + h + 2.0_rp*(REAL(jj,rp) - 1.0_rp)*h
          sum_f = sum_f + fGamma(z,k,t)
       end do

       Inew = 0.5_rp*Iold + sum_f*h
       rerr = ABS((Inew - Iold)/Iold)
       flag = .NOT.(rerr.LT.Tol)
    end do
    IntGamma = Inew
  END FUNCTION IntGamma