Directory of image this file is from
This file as a plain text file

C
C     ..................................................................
C
C        SUBROUTINE INUE
C
C        PURPOSE
C           COMPUTE THE MODIFIED BESSEL FUNCTIONS I FOR ORDERS 1 TO N
C
C        USAGE
C           CALL INUE(X,N,ZI,RI)
C
C        DESCRIPTION OF PARAMETERS
C           X     -GIVEN ARGUMENT OF THE BESSEL FUNCTIONS I
C           N     -GIVEN MAXIMUM ORDER OF BESSEL FUNCTIONS I
C           ZI    -GIVEN VALUE OF BESSEL FUNCTION I OF ORDER ZERO
C                  FOR ARGUMENT X
C           RI    -RESULTANT VECTOR OF DIMENSION N, CONTAINING THE
C                  VALUES OF THE FUNCTIONS I FOR ORDERS 1 TO N
C
C        REMARKS
C           THE VALUE OF ZI MAY BE CALCULATED USING SUBROUTINE I0.
C           USING A DIFFERENT VALUE HAS THE EFFECT THAT ALL VALUES OF
C           BESSEL FUNCTIONS I ARE MULTIPLIED BY THE  FACTOR ZI/I(0,X)
C           WHERE I(0,X) IS THE VALUE OF I FOR ORDER 0 AND ARGUMENT X.
C           THIS MAY BE USED DISADVANTAGEOUSLY IF ONLY THE RATIOS OF I
C           FOR DIFFERENT ORDERS ARE REQUIRED.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           THE VALUES ARE OBTAINED USING BACKWARD RECURRENCE RELATION
C           TECHNIQUE. THE RATIO I(N+1,X)/I(N,X) IS OBTAINED FROM A
C           CONTINUED FRACTION.
C           FOR REFERENCE SEE
C           G. BLANCH,'NUMERICAL EVALUATION OF CONTINUED FRACTIONS',
C           SIAM REVIEW, VOL.6,NO.4,1964,PP.383-421.
C
C     ..................................................................
C
      SUBROUTINE INUE(X,N,ZI,RI)
      DIMENSION RI(1)
      IF(N)10,10,1
    1 FN=N+N
      Q1=X/FN
      IF(ABS(X)-5.E-4)6,6,2
    2 A0=1.
      A1=0.
      B0=0.
      B1=1.
      FI=FN
    3 FI=FI+2.
      AN=FI/ABS(X)
      A=AN*A1+A0
      B=AN*B1+B0
      A0=A1
      B0=B1
      A1=A
      B1=B
      Q0=Q1
      Q1=A/B
      IF(ABS((Q1-Q0)/Q1)-1.E-6)4,4,3
    4 IF(X)5,6,6
    5 Q1=-Q1
    6 K=N
    7 Q1=X/(FN+X*Q1)
      RI(K)=Q1
      FN=FN-2.
      K=K-1
      IF(K)8,8,7
    8 FI=ZI
      DO 9 I=1,N
      FI=FI*RI(I)
    9 RI(I)=FI
   10 RETURN
      END