lxsd / fvn

Blame view

fvn_fnlib/besri.f 2.98 KB
edit raw normal view history
d55dcfb5a   wdaniau   1) Added single p... 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
 
        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