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
|
*DECK DBKISR
SUBROUTINE DBKISR (X, N, SUM, IERR)
C***BEGIN PROLOGUE DBKISR
C***SUBSIDIARY
C***PURPOSE Subsidiary to DBSKIN
C***LIBRARY SLATEC
C***TYPE DOUBLE PRECISION (BKISR-S, DBKISR-D)
C***AUTHOR Amos, D. E., (SNLA)
C***DESCRIPTION
C
C DBKISR computes repeated integrals of the K0 Bessel function
C by the series for N=0,1, and 2.
C
C***SEE ALSO DBSKIN
C***ROUTINES CALLED D1MACH, DPSIXN
C***REVISION HISTORY (YYMMDD)
C 820601 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890911 Removed unnecessary intrinsics. (WRB)
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900328 Added TYPE section. (WRB)
C 910722 Updated AUTHOR section. (ALS)
C***END PROLOGUE DBKISR
INTEGER I, IERR, K, KK, KKN, K1, N, NP
DOUBLE PRECISION AK, ATOL, BK, C, FK, FN, HX, HXS, POL, PR, SUM,
* TKP, TOL, TRM, X, XLN
DOUBLE PRECISION DPSIXN, D1MACH
DIMENSION C(2)
SAVE C
C
DATA C(1), C(2) /1.57079632679489662D+00,1.0D0/
C***FIRST EXECUTABLE STATEMENT DBKISR
IERR=0
TOL = MAX(D1MACH(4),1.0D-18)
IF (X.LT.TOL) GO TO 50
PR = 1.0D0
POL = 0.0D0
IF (N.EQ.0) GO TO 20
DO 10 I=1,N
POL = -POL*X + C(I)
PR = PR*X/I
10 CONTINUE
20 CONTINUE
HX = X*0.5D0
HXS = HX*HX
XLN = LOG(HX)
NP = N + 1
TKP = 3.0D0
FK = 2.0D0
FN = N
BK = 4.0D0
AK = 2.0D0/((FN+1.0D0)*(FN+2.0D0))
SUM = AK*(DPSIXN(N+3)-DPSIXN(3)+DPSIXN(2)-XLN)
ATOL = SUM*TOL*0.75D0
DO 30 K=2,20
AK = AK*(HXS/BK)*((TKP+1.0D0)/(TKP+FN+1.0D0))*(TKP/(TKP+FN))
K1 = K + 1
KK = K1 + K
KKN = KK + N
TRM = (DPSIXN(K1)+DPSIXN(KKN)-DPSIXN(KK)-XLN)*AK
SUM = SUM + TRM
IF (ABS(TRM).LE.ATOL) GO TO 40
TKP = TKP + 2.0D0
BK = BK + TKP
FK = FK + 1.0D0
30 CONTINUE
GO TO 80
40 CONTINUE
SUM = (SUM*HXS+DPSIXN(NP)-XLN)*PR
IF (N.EQ.1) SUM = -SUM
SUM = POL + SUM
RETURN
C-----------------------------------------------------------------------
C SMALL X CASE, X.LT.WORD TOLERANCE
C-----------------------------------------------------------------------
50 CONTINUE
IF (N.GT.0) GO TO 60
HX = X*0.5D0
SUM = DPSIXN(1) - LOG(HX)
RETURN
60 CONTINUE
SUM = C(N)
RETURN
80 CONTINUE
IERR=2
RETURN
END
|
