IntBesselK Function

  1. korc_simple_equilibrium_pdf.f90
  2. korc_simple_equilibrium_pdf
  3. IntBesselK

private function IntBesselK(a, b)

Arguments

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

Return Value real(kind=rp)


Calls

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

Called by

proc~~intbesselk~~CalledByGraph proc~intbesselk IntBesselK proc~p_integral P_integral proc~p_integral->proc~intbesselk proc~pr PR proc~pr->proc~p_integral

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