CalcIBesselFunc Subroutine

public subroutine CalcIBesselFunc(BessFuncArg, BessFuncOrd, IBessFunc, ErrorCode)

Arguments

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

Help
Source Code

CalcIBesselFunc
Source Code

SUBROUTINE CalcIBesselFunc(BessFuncArg,BessFuncOrd,IBessFunc,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 modified Bessel Function from order 0 to BessFuncOrd
          ! BessFuncArg    ARGUMENT OF BESSEL FUNCTION
          ! BessFuncOrd    ORDER OF BESSEL FUNCTION, GREATER THAN OR EQUAL TO ZERO
          ! IBessFunc   RESULTANT VALUE OF I BESSEL FUNCTION
          ! ErrorCode  RESULTANT ERROR CODE:
          !       ErrorCode = 0   NO ERROR
          !       ErrorCode = 1   BessFuncOrd .LT. 0
          !       ErrorCode = 2   BessFuncArg .LT. 0
          !       ErrorCode = 3   IBessFunc .LT. 10**(-30),     IBessFunc IS SET TO 0
          !       ErrorCode = 4   BessFuncArg .GT. BessFuncOrd & BessFuncArg .GT. 90,  IBessFunc IS SET TO 10**38

          ! REFERENCES:
          ! First found in MODSIM.

          ! USE STATEMENTS:
          ! na

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

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

          ! SUBROUTINE PARAMETER DEFINITIONS:
    REAL(r64), PARAMETER :: ErrorTol = 1.0d-06

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
    INTEGER :: LoopCount

    REAL(r64) :: FI
    REAL(r64) :: FK
    REAL(r64) :: TERM

    ErrorCode=0
    IBessFunc=1.d0
    IF (BessFuncArg .EQ. 0.0d0 .AND. BessFuncOrd .EQ. 0) RETURN

    IF (BessFuncOrd .LT. 0) THEN
      ErrorCode=1
      RETURN
    ELSE IF (BessFuncArg .LT. 0.0d0) THEN
      ErrorCode = 2
      RETURN
    ELSE IF (BessFuncArg .GT. 12.d0 .AND. BessFuncArg .GT. BessFuncOrd) THEN
      IF (BessFuncArg .GT. 90.d0) THEN
        ErrorCode = 4
        IBessFunc = 10.d0**30
        RETURN
      ENDIF
      TERM = 1.d0
      IBessFunc = 1.d0
      DO LoopCount = 1, 30  !Start of 1st LoopCount Loop
        IF (ABS(TERM) .LE. ABS(ErrorTol * IBessFunc)) THEN
          IBessFunc = IBessFunc * EXP(BessFuncArg) / SQRT(2.d0 * PI * BessFuncArg)
          RETURN
        ENDIF
        TERM = TERM * 0.125d0 / BessFuncArg * ((2 * LoopCount - 1)**2 - &
               4 * BessFuncOrd * BessFuncOrd) / REAL(LoopCount,r64)
        IBessFunc = IBessFunc + TERM
      END DO  ! End of 1st LoopCount loop
    ENDIF

    TERM = 1.d0
    IF (BessFuncOrd .GT. 0) THEN
      DO LoopCount = 1, BessFuncOrd !Start of 2nd LoopCount Loop
        FI = LoopCount
        IF (ABS(TERM) .LT. 1.d-30 * FI / (BessFuncArg * 2.0d0)) THEN
          ErrorCode = 3
          IBessFunc = 0.0d0
          RETURN
        ENDIF
        TERM = TERM * BessFuncArg / (2.d0 * FI)
      END DO !End of 2nd LoopCount loop
    ENDIF

    IBessFunc = TERM
    DO LoopCount = 1,1000 !Start of 3rd LoopCount Loop
      IF (ABS(TERM) .LE. ABS(IBessFunc * ErrorTol)) RETURN
      FK = LoopCount * (BessFuncOrd + LoopCount)
      TERM = TERM * (BessFuncArg**2 / (4.d0 * FK))
      IBessFunc = IBessFunc + TERM
    END DO !End of  3rd LoopCount loop
    RETURN

END SUBROUTINE CalcIBesselFunc
All Procedures
