CalcKBesselFunc Subroutine

public subroutine CalcKBesselFunc(BessFuncArg, BessFuncOrd, KBessFunc, ErrorCode)

Arguments

Type		Name
real(kind=r64)	::	BessFuncArg
integer	::	BessFuncOrd
real(kind=r64)	::	KBessFunc
integer	::	ErrorCode
Called By

Help
Source Code

CalcKBesselFunc
Source Code

SUBROUTINE CalcKBesselFunc(BessFuncArg,BessFuncOrd,KBessFunc,ErrorCode)
          ! SUBROUTINE INFORMATION:
          !       AUTHOR   Unknown
          !       DATE WRITTEN   Unknown
          !       DATE REWRITTEN  April 1997 by Russell D. Taylor, Ph.D.
          !       MODIFIED
          !       RE-ENGINEERED

          ! PURPOSE OF THIS SUBROUTINE:
          ! To calculate the K Bessel Function for a given argument and
          ! order
          !
          !  BessFuncArg    THE ARGUMENT OF THE K BESSEL FUNCTION DESIRED
          !  BessFuncOrd    THE ORDER OF THE K BESSEL FUNCTION DESIRED
          !  KBessFunc   THE RESULTANT K BESSEL FUNCTION
          !  ErrorCode  RESULTANT ERROR CODE:
          !        ErrorCode=0  NO ERROR
          !        ErrorCode=1  BessFuncOrd IS NEGATIVE
          !        ErrorCode=2  BessFuncArg IS ZERO OR NEGATIVE
          !        ErrorCode=3  BessFuncArg .GT. 85, KBessFunc .LT. 10**-38; KBessFunc SET TO 0.
          !        ErrorCode=4  KBessFunc .GT. 10**38; KBessFunc SET TO 10**38
          !
          ! NOTE: BessFuncOrd MUST BE GREATER THAN OR EQUAL TO ZERO
          !
          ! METHOD:
          !  COMPUTES ZERO ORDER AND FIRST ORDER BESSEL FUNCTIONS USING
          !  SERIES APPROXIMATIONS AND THEN COMPUTES BessFuncOrd TH ORDER FUNCTION
          !  USING RECURRENCE RELATION.
          !  RECURRENCE RELATION AND POLYNOMIAL APPROXIMATION TECHNIQUE
          !  AS DESCRIBED BY A.J.M. HITCHCOCK, 'POLYNOMIAL APPROXIMATIONS
          !  TO BESSEL FUNCTIONS OF ORDER ZERO AND ONE AND TO RELATED
          !  FUNCTIONS,' M.T.A.C., V.11, 1957, PP. 86-88, AND G.BessFuncOrd. WATSON,
          !  'A TREATISE ON THE THEORY OF BESSEL FUNCTIONS,' CAMBRIDGE
          !  UNIVERSITY PRESS, 1958, P.62

          ! USE STATEMENTS:
          ! na

  IMPLICIT NONE    ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENT DEFINITIONS:
    REAL(r64)    :: BessFuncArg
    INTEGER :: BessFuncOrd
    REAL(r64)    :: KBessFunc
    INTEGER :: ErrorCode

          ! SUBROUTINE PARAMETER DEFINITIONS:
    REAL(r64), PARAMETER :: GJMAX = 1.0d+38

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
    INTEGER :: LoopCount
    LOGICAL :: StopLoop

    REAL(r64) :: FACT
    REAL(r64) :: G0
    REAL(r64) :: G1
    REAL(r64) :: GJ
    REAL(r64) :: HJ
    REAL(r64), DIMENSION(12) :: T
    REAL(r64) :: X2J


    KBessFunc=0.0d0
    G0=0.0d0
    GJ=0.0d0

    IF (BessFuncOrd.LT.0.0d0) THEN
      ErrorCode=1
      RETURN
    ELSE IF (BessFuncArg.LE.0.0d0) THEN
      ErrorCode=2
      RETURN
    ELSE IF(BessFuncArg.GT.85.d0) THEN
      ErrorCode=3
      KBessFunc=0.0d0
      RETURN
    ENDIF

    ErrorCode=0

  !     Use polynomial approximation if BessFuncArg > 1.

    IF (BessFuncArg.GT.1.0d0) THEN
      T(1)=1.d0/BessFuncArg
      DO LoopCount = 2,12
        T(LoopCount)=T(LoopCount-1)/BessFuncArg
      END DO !End of LoopCount Loop
      IF (BessFuncOrd.NE.1) THEN

  !     Compute K0 using polynomial approximation

        G0=EXP(-BessFuncArg)*(1.2533141d0-.1566642d0*T(1)+.08811128d0*T(2)-.09139095d0    &
           *T(3)+.1344596d0*T(4)-.2299850d0*T(5)+.3792410d0*T(6)-.5247277d0  &
           *T(7)+.5575368d0*T(8)-.4262633d0*T(9)+.2184518d0*T(10)          &
           -.06680977d0*T(11)+.009189383d0*T(12))*SQRT(1.d0/BessFuncArg)
        IF (BessFuncOrd .EQ. 0) THEN
          KBessFunc=G0
          RETURN
        ENDIF
      ENDIF

  !     Compute K1 using polynomial approximation

      G1=EXP(-BessFuncArg)*(1.2533141d0+.4699927d0*T(1)-.1468583d0*T(2)+.1280427d0*T(3)    &
         -.1736432d0*T(4)+.2847618d0*T(5)-.4594342d0*T(6)+.6283381d0*T(7)     &
         -.6632295d0*T(8)+.5050239d0*T(9)-.2581304d0*T(10)+.07880001d0*T(11)  &
         -.01082418d0*T(12))*SQRT(1.d0/BessFuncArg)
      IF (BessFuncOrd .EQ. 1) THEN
        KBessFunc=G1
        RETURN
      ENDIF
    ELSE

  !     Use series expansion if BessFuncArg <= 1.

      IF (BessFuncOrd.NE.1) THEN

  !     Compute K0 using series expansion

        G0=-(.5772157d0+LOG(BessFuncArg/2.d0))
        X2J=1.d0
        FACT=1.d0
        HJ=0.0d0
        DO LoopCount = 1,6
          X2J=X2J*BessFuncArg**2/4.0d0
          FACT=FACT*(1.0d0/REAL(LoopCount,r64))**2
          HJ=HJ+1.0d0/REAL(LoopCount,r64)
          G0=G0+X2J*FACT*(HJ-(.5772157d0+LOG(BessFuncArg/2.d0)))
        END DO !End of LoopCount Loop
        IF (BessFuncOrd.EQ.0.0d0) THEN
          KBessFunc=G0
          RETURN
        ENDIF
      ENDIF

  !     Compute K1 using series expansion

      X2J=BessFuncArg/2.0d0
      FACT=1.d0
      HJ=1.d0
      G1=1.0d0/BessFuncArg+X2J*(.5d0+(.5772157d0+LOG(BessFuncArg/2.d0))-HJ)
      DO LoopCount = 2,8
        X2J=X2J*BessFuncArg**2/4.0d0
        FACT=FACT*(1.0d0/REAL(LoopCount,r64))**2
        HJ=HJ+1.0d0/REAL(LoopCount,r64)
        G1=G1+X2J*FACT*(.5d0+((.5772157d0+LOG(BessFuncArg/2.d0))-HJ)*REAL(LoopCount,r64))
      END DO !End of LoopCount Loop
      IF (BessFuncOrd.EQ.1) THEN
        KBessFunc=G1
        RETURN
      ENDIF
    ENDIF

  !     From K0 and K1 compute KN using recurrence relation

    LoopCount = 2
    StopLoop = .FALSE.
    DO WHILE (LoopCount .LE. BessFuncOrd .AND. .NOT. StopLoop)
      GJ=2.d0*(REAL(LoopCount,r64)-1.0d0)*G1/BessFuncArg+G0
      IF (GJ-GJMAX .GT. 0.0d0) THEN
        ErrorCode=4
        GJ=GJMAX
        StopLoop = .TRUE.
      ELSE
        G0=G1
        G1=GJ
        LoopCount = LoopCount + 1
      END IF
    END DO !End of LoopCount Loop
    KBessFunc=GJ

    RETURN
END SUBROUTINE CalcKBesselFunc
All Procedures
