besri.f 2.98 KB
SUBROUTINE BESRI(X, NB, B)
C
C THIS ROUTINE CALCULATES MODIFIED BESSEL FUNCTIONS I OF REAL
C ARGUMENT AND INTEGER ORDER.
C
C      EXPLANATION OF VARIABLES IN THE CALLING SEQUENCE -
C
C X  -  REAL ARGUMENT FOR WHICH THE I*S
C       ARE TO BE CALCULATED.
C
C NB -  INTEGER TYPE.  1 + HIGHEST ORDER TO BE CALCULATED.
C       IT MUST BE POSITIVE.
C
C B  -  REAL VECTOR OF LENGTH NB, NEED NOT BE
C       INITIALIZED BY USER.  IF THE ROUTINE TERMINATES
C       NORMALLY IT RETURNS I-SUB-ZERO
C       THROUGH I-SUB-NB-MINUS-ONE OF X IN THIS
C       VECTOR.
C
C                  ACCURACY OF THE COMPUTED VALUES -
C
C       IN CASE OF AN ERROR, NOT ALL I*S
C       ARE CALCULATED TO THE DESIRED ACCURACY.
C
C       THE SUBPROGRAM CALLED BY BESRI, B1SLR,
C       RETURNS IN THE VARIABLE, NCALC, THE NUMBER CALCULATED CORRECTLY.
C
C       LET NTEN BE THE LARGEST INTEGER K SUCH THAT 10**K IS MACHINE-
C       REPRESENTABLE IN REAL.
C       THEN NB.GT.NCALC.GT.0 WILL OCCUR IF NB.GT.MAGX AND ABS(I-
C       SUB-NB-OF-X/I-SUB-MAGX+NP-OF-X).LT.10.**(NTEN/2), I.E. NB
C       IS MUCH GREATER THAN MAGX.  IN THIS CASE, B(N) IS CALCU-
C       LATED TO THE DESIRED ACCURACY FOR N.LE.NCALC, BUT FOR
C       NCALC.LT.N.LE.NB, PRECISION IS LOST.  IF N.GT.NCALC AND
C       ABS(B(NCALC)/B(N)).EQ.10**-K, THEN THE LAST K SIGNIFICANT
C       FIGURES OF B(N) ARE ERRONEOUS.  IF THE USER WISHES TO
C       CALCULATE B(N) TO HIGHER ACCURACY, HE SHOULD USE AN
C       ASYMPTOTIC FORMULA FOR LARGE ORDER.
C
      REAL X,B(NB),R1MACH
C
C CHECK INPUT VALUES
C
C       AN UPPER LIMIT OF 10000 IS SET ON THE MAGNITUDE OF X.
C       BEAR IN MIND THAT IF ABS(X)=N, THEN AT LEAST N ITERATIONS
C       OF THE BACKWARD RECURSION WILL BE EXECUTED.
C
      MAGX = ABS(X)
      MEXP = ALOG(R1MACH(2))
C
C/6S
C     IF  (MAGX .GT. 10000 .OR. MAGX .GT. MEXP) CALL SETERR(
C    1    33H BESRI - X IS TOO BIG (MAGNITUDE),33,1,2)
C/7S
      IF  (MAGX .GT. 10000 .OR. MAGX .GT. MEXP) CALL SETERR(
     1    ' BESRI - X IS TOO BIG (MAGNITUDE)',33,1,2)
C/
C
C/6S
C     IF  (NB .LT. 1) CALL SETERR(
C    1    28H BESRI - NB SHOULD = ORDER+1,28,2,2)
C/7S
      IF  (NB .LT. 1) CALL SETERR(
     1    ' BESRI - NB SHOULD = ORDER+1',28,2,2)
C/
C
C BESRJ CALLS ON THE SUBPROGRAM,B1SLR,
C WHICH IS SOOKNES  ORIGINAL BESLRI.
C
C THE ADDITIONAL INPUT ARGUMENTS REQUIRED FOR IT ARE -
C
C IZE   INTEGER TYPE.  ZERO IF J*S ARE TO BE CALCULATED, 1
C       IF I*S ARE TO BE CALCULATED.(THIRD ARGUMENT BELOW)
C
C NCALC INTEGER TYPE, NEED NOT BE INITIALIZED BY USER.
C       BEFORE USING THE RESULTS, IT SHOULD BE CHECKED THAT
C       NCALC=NB, I.E. ALL ORDERS HAVE BEEN CALCULATED TO
C       THE DESIRED ACCURACY.
C
      CALL B1SLR (X, NB, 1, B, NCALC)
C
C TEST IF ALL GOT COMPUTED OK
C (SINCE SOME VALUES MAY BE OK, THIS IS A RECOVERABLE ERROR.)
C
      IF (NB .EQ. NCALC) RETURN
C
      NCALC = NCALC+10
C/6S
C     CALL SETERR(
C    1    38H BESRI - ONLY THIS MANY ANSWERS ARE OK,38,NCALC,1)
C/7S
      CALL SETERR(
     1    ' BESRI - ONLY THIS MANY ANSWERS ARE OK',38,NCALC,1)
C/
C
      RETURN
      END