-
git-svn-id: https://lxsd.femto-st.fr/svn/fvn@58 b657c933-2333-4658-acf2-d3c7c2708721
dbesri.f
3.04 KB
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 DBESRI(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 - DOUBLE PRECISION 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 - DOUBLE PRECISION 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 DBESRI, DB1SLR,
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 DOUBLE PRECISION.
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
DOUBLE PRECISION X,B(NB),D1MACH, DLOG
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 = DABS(X)
MEXP = DLOG(D1MACH(2))
C
C/6S
C IF (MAGX .GT. 10000 .OR. MAGX .GT. MEXP) CALL SETERR(
C 1 33HDBESRI - X IS TOO BIG (MAGNITUDE),33,1,2)
C/7S
IF (MAGX .GT. 10000 .OR. MAGX .GT. MEXP) CALL SETERR(
1 'DBESRI - X IS TOO BIG (MAGNITUDE)',33,1,2)
C/
C
C/6S
C IF (NB .LT. 1) CALL SETERR(
C 1 28HDBESRI - NB SHOULD = ORDER+1,28,2,2)
C/7S
IF (NB .LT. 1) CALL SETERR(
1 'DBESRI - NB SHOULD = ORDER+1',28,2,2)
C/
C
C DBESRJ CALLS ON THE SUBPROGRAM,DB1SLR,
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 DB1SLR (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 38HDBESRI - ONLY THIS MANY ANSWERS ARE OK,38,NCALC,1)
C/7S
CALL SETERR(
1 'DBESRI - ONLY THIS MANY ANSWERS ARE OK',38,NCALC,1)
C/
C
RETURN
END