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