1*DECK DBKIAS
2      SUBROUTINE DBKIAS (X, N, KTRMS, T, ANS, IND, MS, GMRN, H, IERR)
3C***BEGIN PROLOGUE  DBKIAS
4C***SUBSIDIARY
5C***PURPOSE  Subsidiary to DBSKIN
6C***LIBRARY   SLATEC
7C***TYPE      DOUBLE PRECISION (BKIAS-S, DBKIAS-D)
8C***AUTHOR  Amos, D. E., (SNLA)
9C***DESCRIPTION
10C
11C     DBKIAS computes repeated integrals of the K0 Bessel function
12C     by the asymptotic expansion
13C
14C***SEE ALSO  DBSKIN
15C***ROUTINES CALLED  D1MACH, DBDIFF, DGAMRN, DHKSEQ
16C***REVISION HISTORY  (YYMMDD)
17C   820601  DATE WRITTEN
18C   890531  Changed all specific intrinsics to generic.  (WRB)
19C   890911  Removed unnecessary intrinsics.  (WRB)
20C   891214  Prologue converted to Version 4.0 format.  (BAB)
21C   900328  Added TYPE section.  (WRB)
22C   910722  Updated AUTHOR section.  (ALS)
23C***END PROLOGUE  DBKIAS
24      INTEGER I, II, IND, J, JMI, JN, K, KK, KM, KTRMS, MM, MP, MS, N,
25     * IERR
26      DOUBLE PRECISION ANS, B, BND, DEN1, DEN2, DEN3, ER, ERR, FJ, FK,
27     * FLN, FM1, GMRN, G1, GS, H, HN, HRTPI, RAT, RG1, RXP, RZ, RZX, S,
28     * SS, SUMI, SUMJ, T, TOL, V, W, X, XP, Z
29      DOUBLE PRECISION DGAMRN, D1MACH
30      DIMENSION B(120), XP(16), S(31), H(*), V(52), W(52), T(50),
31     * BND(15)
32      SAVE B, BND, HRTPI
33C-----------------------------------------------------------------------
34C             COEFFICIENTS OF POLYNOMIAL P(J-1,X), J=1,15
35C-----------------------------------------------------------------------
36      DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10),
37     * B(11), B(12), B(13), B(14), B(15), B(16), B(17), B(18), B(19),
38     * B(20), B(21), B(22), B(23), B(24) /1.00000000000000000D+00,
39     * 1.00000000000000000D+00,-2.00000000000000000D+00,
40     * 1.00000000000000000D+00,-8.00000000000000000D+00,
41     * 6.00000000000000000D+00,1.00000000000000000D+00,
42     * -2.20000000000000000D+01,5.80000000000000000D+01,
43     * -2.40000000000000000D+01,1.00000000000000000D+00,
44     * -5.20000000000000000D+01,3.28000000000000000D+02,
45     * -4.44000000000000000D+02,1.20000000000000000D+02,
46     * 1.00000000000000000D+00,-1.14000000000000000D+02,
47     * 1.45200000000000000D+03,-4.40000000000000000D+03,
48     * 3.70800000000000000D+03,-7.20000000000000000D+02,
49     * 1.00000000000000000D+00,-2.40000000000000000D+02,
50     * 5.61000000000000000D+03/
51      DATA B(25), B(26), B(27), B(28), B(29), B(30), B(31), B(32),
52     * B(33), B(34), B(35), B(36), B(37), B(38), B(39), B(40), B(41),
53     * B(42), B(43), B(44), B(45), B(46), B(47), B(48)
54     * /-3.21200000000000000D+04,5.81400000000000000D+04,
55     * -3.39840000000000000D+04,5.04000000000000000D+03,
56     * 1.00000000000000000D+00,-4.94000000000000000D+02,
57     * 1.99500000000000000D+04,-1.95800000000000000D+05,
58     * 6.44020000000000000D+05,-7.85304000000000000D+05,
59     * 3.41136000000000000D+05,-4.03200000000000000D+04,
60     * 1.00000000000000000D+00,-1.00400000000000000D+03,
61     * 6.72600000000000000D+04,-1.06250000000000000D+06,
62     * 5.76550000000000000D+06,-1.24400640000000000D+07,
63     * 1.10262960000000000D+07,-3.73392000000000000D+06,
64     * 3.62880000000000000D+05,1.00000000000000000D+00,
65     * -2.02600000000000000D+03,2.18848000000000000D+05/
66      DATA B(49), B(50), B(51), B(52), B(53), B(54), B(55), B(56),
67     * B(57), B(58), B(59), B(60), B(61), B(62), B(63), B(64), B(65),
68     * B(66), B(67), B(68), B(69), B(70), B(71), B(72)
69     * /-5.32616000000000000D+06,4.47650000000000000D+07,
70     * -1.55357384000000000D+08,2.38904904000000000D+08,
71     * -1.62186912000000000D+08,4.43390400000000000D+07,
72     * -3.62880000000000000D+06,1.00000000000000000D+00,
73     * -4.07200000000000000D+03,6.95038000000000000D+05,
74     * -2.52439040000000000D+07,3.14369720000000000D+08,
75     * -1.64838430400000000D+09,4.00269508800000000D+09,
76     * -4.64216395200000000D+09,2.50748121600000000D+09,
77     * -5.68356480000000000D+08,3.99168000000000000D+07,
78     * 1.00000000000000000D+00,-8.16600000000000000D+03,
79     * 2.17062600000000000D+06,-1.14876376000000000D+08,
80     * 2.05148277600000000D+09,-1.55489607840000000D+10/
81      DATA B(73), B(74), B(75), B(76), B(77), B(78), B(79), B(80),
82     * B(81), B(82), B(83), B(84), B(85), B(86), B(87), B(88), B(89),
83     * B(90), B(91), B(92), B(93), B(94), B(95), B(96)
84     * /5.60413987840000000D+10,-1.01180433024000000D+11,
85     * 9.21997902240000000D+10,-4.07883018240000000D+10,
86     * 7.82771904000000000D+09,-4.79001600000000000D+08,
87     * 1.00000000000000000D+00,-1.63560000000000000D+04,
88     * 6.69969600000000000D+06,-5.07259276000000000D+08,
89     * 1.26698177760000000D+10,-1.34323420224000000D+11,
90     * 6.87720046384000000D+11,-1.81818864230400000D+12,
91     * 2.54986547342400000D+12,-1.88307966182400000D+12,
92     * 6.97929436800000000D+11,-1.15336085760000000D+11,
93     * 6.22702080000000000D+09,1.00000000000000000D+00,
94     * -3.27380000000000000D+04,2.05079880000000000D+07,
95     * -2.18982980800000000D+09,7.50160522280000000D+10/
96      DATA B(97), B(98), B(99), B(100), B(101), B(102), B(103), B(104),
97     * B(105), B(106), B(107), B(108), B(109), B(110), B(111), B(112),
98     * B(113), B(114), B(115), B(116), B(117), B(118)
99     * /-1.08467651241600000D+12,7.63483214939200000D+12,
100     * -2.82999100661120000D+13,5.74943734645920000D+13,
101     * -6.47283751398720000D+13,3.96895780558080000D+13,
102     * -1.25509040179200000D+13,1.81099255680000000D+12,
103     * -8.71782912000000000D+10,1.00000000000000000D+00,
104     * -6.55040000000000000D+04,6.24078900000000000D+07,
105     * -9.29252692000000000D+09,4.29826006340000000D+11,
106     * -8.30844432796800000D+12,7.83913848313120000D+13,
107     * -3.94365587815520000D+14,1.11174747256968000D+15,
108     * -1.79717122069056000D+15,1.66642448627145600D+15,
109     * -8.65023253219584000D+14,2.36908271543040000D+14/
110      DATA B(119), B(120) /-3.01963769856000000D+13,
111     * 1.30767436800000000D+12/
112C-----------------------------------------------------------------------
113C             BOUNDS B(M,K) , K=M-3
114C-----------------------------------------------------------------------
115      DATA BND(1), BND(2), BND(3), BND(4), BND(5), BND(6), BND(7),
116     * BND(8), BND(9), BND(10), BND(11), BND(12), BND(13), BND(14),
117     * BND(15) /1.0D0,1.0D0,1.0D0,1.0D0,3.10D0,5.18D0,11.7D0,29.8D0,
118     * 90.4D0,297.0D0,1070.0D0,4290.0D0,18100.0D0,84700.0D0,408000.0D0/
119      DATA HRTPI /8.86226925452758014D-01/
120C
121C***FIRST EXECUTABLE STATEMENT  DBKIAS
122      IERR=0
123      TOL = MAX(D1MACH(4),1.0D-18)
124      FLN = N
125      RZ = 1.0D0/(X+FLN)
126      RZX = X*RZ
127      Z = 0.5D0*(X+FLN)
128      IF (IND.GT.1) GO TO 10
129      GMRN = DGAMRN(Z)
130   10 CONTINUE
131      GS = HRTPI*GMRN
132      G1 = GS + GS
133      RG1 = 1.0D0/G1
134      GMRN = (RZ+RZ)/GMRN
135      IF (IND.GT.1) GO TO 70
136C-----------------------------------------------------------------------
137C     EVALUATE ERROR FOR M=MS
138C-----------------------------------------------------------------------
139      HN = 0.5D0*FLN
140      DEN2 = KTRMS + KTRMS + N
141      DEN3 = DEN2 - 2.0D0
142      DEN1 = X + DEN2
143      ERR = RG1*(X+X)/(DEN1-1.0D0)
144      IF (N.EQ.0) GO TO 20
145      RAT = 1.0D0/(FLN*FLN)
146   20 CONTINUE
147      IF (KTRMS.EQ.0) GO TO 30
148      FJ = KTRMS
149      RAT = 0.25D0/(HRTPI*DEN3*SQRT(FJ))
150   30 CONTINUE
151      ERR = ERR*RAT
152      FJ = -3.0D0
153      DO 50 J=1,15
154        IF (J.LE.5) ERR = ERR/DEN1
155        FM1 = MAX(1.0D0,FJ)
156        FJ = FJ + 1.0D0
157        ER = BND(J)*ERR
158        IF (KTRMS.EQ.0) GO TO 40
159        ER = ER/FM1
160        IF (ER.LT.TOL) GO TO 60
161        IF (J.GE.5) ERR = ERR/DEN3
162        GO TO 50
163   40   CONTINUE
164        ER = ER*(1.0D0+HN/FM1)
165        IF (ER.LT.TOL) GO TO 60
166        IF (J.GE.5) ERR = ERR/FLN
167   50 CONTINUE
168      GO TO 200
169   60 CONTINUE
170      MS = J
171   70 CONTINUE
172      MM = MS + MS
173      MP = MM + 1
174C-----------------------------------------------------------------------
175C     H(K)=(-Z)**(K)*(PSI(K-1,Z)-PSI(K-1,Z+0.5))/GAMMA(K) , K=1,2,...,MM
176C-----------------------------------------------------------------------
177      IF (IND.GT.1) GO TO 80
178      CALL DHKSEQ(Z, MM, H, IERR)
179      GO TO 100
180   80 CONTINUE
181      RAT = Z/(Z-0.5D0)
182      RXP = RAT
183      DO 90 I=1,MM
184        H(I) = RXP*(1.0D0-H(I))
185        RXP = RXP*RAT
186   90 CONTINUE
187  100 CONTINUE
188C-----------------------------------------------------------------------
189C     SCALED S SEQUENCE
190C-----------------------------------------------------------------------
191      S(1) = 1.0D0
192      FK = 1.0D0
193      DO 120 K=2,MP
194        SS = 0.0D0
195        KM = K - 1
196        I = KM
197        DO 110 J=1,KM
198          SS = SS + S(J)*H(I)
199          I = I - 1
200  110   CONTINUE
201        S(K) = SS/FK
202        FK = FK + 1.0D0
203  120 CONTINUE
204C-----------------------------------------------------------------------
205C     SCALED S-TILDA SEQUENCE
206C-----------------------------------------------------------------------
207      IF (KTRMS.EQ.0) GO TO 160
208      FK = 0.0D0
209      SS = 0.0D0
210      RG1 = RG1/Z
211      DO 130 K=1,KTRMS
212        V(K) = Z/(Z+FK)
213        W(K) = T(K)*V(K)
214        SS = SS + W(K)
215        FK = FK + 1.0D0
216  130 CONTINUE
217      S(1) = S(1) - SS*RG1
218      DO 150 I=2,MP
219        SS = 0.0D0
220        DO 140 K=1,KTRMS
221          W(K) = W(K)*V(K)
222          SS = SS + W(K)
223  140   CONTINUE
224        S(I) = S(I) - SS*RG1
225  150 CONTINUE
226  160 CONTINUE
227C-----------------------------------------------------------------------
228C     SUM ON J
229C-----------------------------------------------------------------------
230      SUMJ = 0.0D0
231      JN = 1
232      RXP = 1.0D0
233      XP(1) = 1.0D0
234      DO 190 J=1,MS
235        JN = JN + J - 1
236        XP(J+1) = XP(J)*RZX
237        RXP = RXP*RZ
238C-----------------------------------------------------------------------
239C     SUM ON I
240C-----------------------------------------------------------------------
241        SUMI = 0.0D0
242        II = JN
243        DO 180 I=1,J
244          JMI = J - I + 1
245          KK = J + I + 1
246          DO 170 K=1,JMI
247            V(K) = S(KK)*XP(K)
248            KK = KK + 1
249  170     CONTINUE
250          CALL DBDIFF(JMI, V)
251          SUMI = SUMI + B(II)*V(JMI)*XP(I+1)
252          II = II + 1
253  180   CONTINUE
254        SUMJ = SUMJ + SUMI*RXP
255  190 CONTINUE
256      ANS = GS*(S(1)-SUMJ)
257      RETURN
258  200 CONTINUE
259      IERR=2
260      RETURN
261      END
262