IntBesselK Function

  1. korc_experimental_pdf.f90
  2. korc_experimental_pdf
  3. IntBesselK

private function IntBesselK(a, b)

@brief Extended trapezoidal rule for integrating the modified Bessel function of second kind. 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

Return Value real(kind=rp)


Calls

proc~~intbesselk~2~~CallsGraph proc~intbesselk~2 IntBesselK proc~besselk~2 besselk proc~intbesselk~2->proc~besselk~2

Called by

proc~~intbesselk~2~~CalledByGraph proc~intbesselk~2 IntBesselK proc~p_integral~2 P_integral proc~p_integral~2->proc~intbesselk~2 proc~pr~2 PR proc~pr~2->proc~p_integral~2 proc~frexpr fRExPR proc~frexpr->proc~pr~2 proc~fre_hxpr fRE_HxPR proc~fre_hxpr->proc~pr~2

Contents

Source Code

IntBesselK

Source Code

  FUNCTION IntBesselK(a,b)
    REAL(rp), INTENT(IN) :: a
    REAL(rp), INTENT(IN) :: b
    REAL(rp) :: IntBesselK
    REAL(rp) :: Iold
    REAL(rp) :: Inew
    REAL(rp) :: rerr
    REAL(rp) :: sum_f
    REAL(rp) :: v,h,z
    INTEGER :: ii,jj,npoints
    LOGICAL :: flag

    v = 5.0_rp/3.0_rp
    h = b - a
    sum_f = 0.5*(besselk(v,a) + besselk(v,b))

    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 + besselk(v,z)
       end do

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