File: BESI.FT of Tape: Various/ETH/eth11-1

(Source file text) 

---

C
C     ..................................................................
C
C        SUBROUTINE BESI
C
C        PURPOSE
C           COMPUTE THE I BESSEL FUNCTION FOR A GIVEN ARGUMENT AND ORDER
C
C        USAGE
C           CALL BESI(X,N,BI,IER)
C
C        DESCRIPTION OF PARAMETERS
C           X  -THE ARGUMENT OF THE I BESSEL FUNCTION DESIRED
C           N  -THE ORDER OF THE I BESSEL FUNCTION DESIRED
C           BI -THE RESULTANT I BESSEL FUNCTION
C           IER-RESULTANT ERROR CODE WHERE
C              IER=0 NO ERROR
C              IER=1 N IS NEGATIVE
C              IER=2 X IS NEGATIVE
C              IER=3 UNDERFLOW, BI .LT. 1.E-38, BI SET TO 0.0
C              IER=4 OVERFLOW, X .GT. 60 WHERE X .GT. N
C
C        REMARKS
C           N AND X MUST BE .GE. ZERO
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           COMPUTES I BESSEL FUNCTION USING SERIES OR ASYMPTOTIC
C           APPROXIMATION DEPENDING ON RANGE OF ARGUMENTS.
C
C     ..................................................................
C
      SUBROUTINE BESI(X,N, BI,IER)
C
C     CHECK FOR ERRORS IN N AND X AND EXIT IF ANY ARE PRESENT
C
      IER=0
      BI=1.0
      IF(N)150,15,10
   10 IF(X)160,20,20
   15 IF(X)160,17,20
   17 RETURN
C
C     DEFINE TOLERANCE
C
   20 TOL=1.E-6
C
C     IF ARGUMENT GT 12 AND GT N, USE ASYMPTOTIC FORM
C
      IF(X-12.)40,40,30
   30 IF(X-FLOAT(N))40,40,110
C
C     COMPUTE FIRST TERM OF SERIES AND SET INITIAL VALUE OF THE SUM
C
   40 XX=X/2.
   50 TERM=1.0
      IF(N) 70,70,55
   55 DO 60 I=1,N
      FI=I
      IF(ABS(TERM)-1.E-36)56,60,60
   56 IER=3
      BI=0.0
      RETURN
   60 TERM=TERM*XX/FI
   70 BI=TERM
      XX=XX*XX
C
C     COMPUTE TERMS, STOPPING WHEN ABS(TERM) LE ABS(SUM OF TERMS)
C     TIMES TOLERANCE
C
      DO 90 K=1,1000
      IF(ABS(TERM)-ABS(BI*TOL))100,100,80
   80 FK=K*(N+K)
      TERM=TERM*(XX/FK)
   90 BI=BI+TERM
C
C     RETURN BI AS ANSWER
C
  100 RETURN
C
C     X GT 12 AND X GT N, SO USE ASYMPTOTIC APPROXIMATION
C
  110 FN=4*N*N
C  THE LIMIT OF 60. MAY BE ABLE TO BE REVISED UPWARD
      IF(X-60.0)115,111,111
  111 IER=4
      RETURN
  115 XX=1./(8.*X)
      TERM=1.
      BI=1.
      DO 130 K=1,30
      IF(ABS(TERM)-ABS(TOL*BI))140,140,120
  120 FK=(2*K-1)**2
      TERM=TERM*XX*(FK-FN)/FLOAT(K)
  130 BI=BI+TERM
C
C     SIGNIFICANCE LOST AFTER 30 TERMS, TRY SERIES
C
      GO TO 40
  140 PI=3.141592653
      BI=BI*EXP(X)/SQRT(2.*PI*X)
      GO TO 100
  150 IER=1
      GO TO 100
  160 IER=2
      GO TO 100
      END

---