1 SUBROUTINE ZUNHJ(ZR, ZI, FNU, IPMTR, TOL, PHIR, PHII, ARGR, ARGI, 2 * ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) 3C***BEGIN PROLOGUE ZUNHJ 4C***REFER TO ZBESI,ZBESK 5C 6C REFERENCES 7C HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ AND I.A. 8C STEGUN, AMS55, NATIONAL BUREAU OF STANDARDS, 1965, CHAPTER 9. 9C 10C ASYMPTOTICS AND SPECIAL FUNCTIONS BY F.W.J. OLVER, ACADEMIC 11C PRESS, N.Y., 1974, PAGE 420 12C 13C ABSTRACT 14C ZUNHJ COMPUTES PARAMETERS FOR BESSEL FUNCTIONS C(FNU,Z) = 15C J(FNU,Z), Y(FNU,Z) OR H(I,FNU,Z) I=1,2 FOR LARGE ORDERS FNU 16C BY MEANS OF THE UNIFORM ASYMPTOTIC EXPANSION 17C 18C C(FNU,Z)=C1*PHI*( ASUM*AIRY(ARG) + C2*BSUM*DAIRY(ARG) ) 19C 20C FOR PROPER CHOICES OF C1, C2, AIRY AND DAIRY WHERE AIRY IS 21C AN AIRY FUNCTION AND DAIRY IS ITS DERIVATIVE. 22C 23C (2/3)*FNU*ZETA**1.5 = ZETA1-ZETA2, 24C 25C ZETA1=0.5*FNU*CLOG((1+W)/(1-W)), ZETA2=FNU*W FOR SCALING 26C PURPOSES IN AIRY FUNCTIONS FROM CAIRY OR CBIRY. 27C 28C MCONJ=SIGN OF AIMAG(Z), BUT IS AMBIGUOUS WHEN Z IS REAL AND 29C MUST BE SPECIFIED. IPMTR=0 RETURNS ALL PARAMETERS. IPMTR= 30C 1 COMPUTES ALL EXCEPT ASUM AND BSUM. 31C 32C***ROUTINES CALLED ZABS,ZDIV,ZLOG,ZSQRTAMOS,D1MACH 33C***END PROLOGUE ZUNHJ 34C COMPLEX ARG,ASUM,BSUM,CFNU,CONE,CR,CZERO,DR,P,PHI,PRZTH,PTFN, 35C *RFN13,RTZTA,RZTH,SUMA,SUMB,TFN,T2,UP,W,W2,Z,ZA,ZB,ZC,ZETA,ZETA1, 36C *ZETA2,ZTH 37 EXTERNAL ZABS 38 DOUBLE PRECISION ALFA, ANG, AP, AR, ARGI, ARGR, ASUMI, ASUMR, 39 * ATOL, AW2, AZTH, BETA, BR, BSUMI, BSUMR, BTOL, C, CONEI, CONER, 40 * CRI, CRR, DRI, DRR, EX1, EX2, FNU, FN13, FN23, GAMA, GPI, HPI, 41 * PHII, PHIR, PI, PP, PR, PRZTHI, PRZTHR, PTFNI, PTFNR, RAW, RAW2, 42 * RAZTH, RFNU, RFNU2, RFN13, RTZTI, RTZTR, RZTHI, RZTHR, STI, STR, 43 * SUMAI, SUMAR, SUMBI, SUMBR, TEST, TFNI, TFNR, THPI, TOL, TZAI, 44 * TZAR, T2I, T2R, UPI, UPR, WI, WR, W2I, W2R, ZAI, ZAR, ZBI, ZBR, 45 * ZCI, ZCR, ZEROI, ZEROR, ZETAI, ZETAR, ZETA1I, ZETA1R, ZETA2I, 46 * ZETA2R, ZI, ZR, ZTHI, ZTHR, ZABS, AC, D1MACH 47 INTEGER IAS, IBS, IPMTR, IS, J, JR, JU, K, KMAX, KP1, KS, L, LR, 48 * LRP1, L1, L2, M, IDUM 49 DIMENSION AR(14), BR(14), C(105), ALFA(180), BETA(210), GAMA(30), 50 * AP(30), PR(30), PI(30), UPR(14), UPI(14), CRR(14), CRI(14), 51 * DRR(14), DRI(14) 52 DATA AR(1), AR(2), AR(3), AR(4), AR(5), AR(6), AR(7), AR(8), 53 1 AR(9), AR(10), AR(11), AR(12), AR(13), AR(14)/ 54 2 1.00000000000000000D+00, 1.04166666666666667D-01, 55 3 8.35503472222222222D-02, 1.28226574556327160D-01, 56 4 2.91849026464140464D-01, 8.81627267443757652D-01, 57 5 3.32140828186276754D+00, 1.49957629868625547D+01, 58 6 7.89230130115865181D+01, 4.74451538868264323D+02, 59 7 3.20749009089066193D+03, 2.40865496408740049D+04, 60 8 1.98923119169509794D+05, 1.79190200777534383D+06/ 61 DATA BR(1), BR(2), BR(3), BR(4), BR(5), BR(6), BR(7), BR(8), 62 1 BR(9), BR(10), BR(11), BR(12), BR(13), BR(14)/ 63 2 1.00000000000000000D+00, -1.45833333333333333D-01, 64 3 -9.87413194444444444D-02, -1.43312053915895062D-01, 65 4 -3.17227202678413548D-01, -9.42429147957120249D-01, 66 5 -3.51120304082635426D+00, -1.57272636203680451D+01, 67 6 -8.22814390971859444D+01, -4.92355370523670524D+02, 68 7 -3.31621856854797251D+03, -2.48276742452085896D+04, 69 8 -2.04526587315129788D+05, -1.83844491706820990D+06/ 70 DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10), 71 1 C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18), 72 2 C(19), C(20), C(21), C(22), C(23), C(24)/ 73 3 1.00000000000000000D+00, -2.08333333333333333D-01, 74 4 1.25000000000000000D-01, 3.34201388888888889D-01, 75 5 -4.01041666666666667D-01, 7.03125000000000000D-02, 76 6 -1.02581259645061728D+00, 1.84646267361111111D+00, 77 7 -8.91210937500000000D-01, 7.32421875000000000D-02, 78 8 4.66958442342624743D+00, -1.12070026162229938D+01, 79 9 8.78912353515625000D+00, -2.36408691406250000D+00, 80 A 1.12152099609375000D-01, -2.82120725582002449D+01, 81 B 8.46362176746007346D+01, -9.18182415432400174D+01, 82 C 4.25349987453884549D+01, -7.36879435947963170D+00, 83 D 2.27108001708984375D-01, 2.12570130039217123D+02, 84 E -7.65252468141181642D+02, 1.05999045252799988D+03/ 85 DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32), 86 1 C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40), 87 2 C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/ 88 3 -6.99579627376132541D+02, 2.18190511744211590D+02, 89 4 -2.64914304869515555D+01, 5.72501420974731445D-01, 90 5 -1.91945766231840700D+03, 8.06172218173730938D+03, 91 6 -1.35865500064341374D+04, 1.16553933368645332D+04, 92 7 -5.30564697861340311D+03, 1.20090291321635246D+03, 93 8 -1.08090919788394656D+02, 1.72772750258445740D+00, 94 9 2.02042913309661486D+04, -9.69805983886375135D+04, 95 A 1.92547001232531532D+05, -2.03400177280415534D+05, 96 B 1.22200464983017460D+05, -4.11926549688975513D+04, 97 C 7.10951430248936372D+03, -4.93915304773088012D+02, 98 D 6.07404200127348304D+00, -2.42919187900551333D+05, 99 E 1.31176361466297720D+06, -2.99801591853810675D+06/ 100 DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56), 101 1 C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64), 102 2 C(65), C(66), C(67), C(68), C(69), C(70), C(71), C(72)/ 103 3 3.76327129765640400D+06, -2.81356322658653411D+06, 104 4 1.26836527332162478D+06, -3.31645172484563578D+05, 105 5 4.52187689813627263D+04, -2.49983048181120962D+03, 106 6 2.43805296995560639D+01, 3.28446985307203782D+06, 107 7 -1.97068191184322269D+07, 5.09526024926646422D+07, 108 8 -7.41051482115326577D+07, 6.63445122747290267D+07, 109 9 -3.75671766607633513D+07, 1.32887671664218183D+07, 110 A -2.78561812808645469D+06, 3.08186404612662398D+05, 111 B -1.38860897537170405D+04, 1.10017140269246738D+02, 112 C -4.93292536645099620D+07, 3.25573074185765749D+08, 113 D -9.39462359681578403D+08, 1.55359689957058006D+09, 114 E -1.62108055210833708D+09, 1.10684281682301447D+09/ 115 DATA C(73), C(74), C(75), C(76), C(77), C(78), C(79), C(80), 116 1 C(81), C(82), C(83), C(84), C(85), C(86), C(87), C(88), 117 2 C(89), C(90), C(91), C(92), C(93), C(94), C(95), C(96)/ 118 3 -4.95889784275030309D+08, 1.42062907797533095D+08, 119 4 -2.44740627257387285D+07, 2.24376817792244943D+06, 120 5 -8.40054336030240853D+04, 5.51335896122020586D+02, 121 6 8.14789096118312115D+08, -5.86648149205184723D+09, 122 7 1.86882075092958249D+10, -3.46320433881587779D+10, 123 8 4.12801855797539740D+10, -3.30265997498007231D+10, 124 9 1.79542137311556001D+10, -6.56329379261928433D+09, 125 A 1.55927986487925751D+09, -2.25105661889415278D+08, 126 B 1.73951075539781645D+07, -5.49842327572288687D+05, 127 C 3.03809051092238427D+03, -1.46792612476956167D+10, 128 D 1.14498237732025810D+11, -3.99096175224466498D+11, 129 E 8.19218669548577329D+11, -1.09837515608122331D+12/ 130 DATA C(97), C(98), C(99), C(100), C(101), C(102), C(103), C(104), 131 1 C(105)/ 132 2 1.00815810686538209D+12, -6.45364869245376503D+11, 133 3 2.87900649906150589D+11, -8.78670721780232657D+10, 134 4 1.76347306068349694D+10, -2.16716498322379509D+09, 135 5 1.43157876718888981D+08, -3.87183344257261262D+06, 136 6 1.82577554742931747D+04/ 137 DATA ALFA(1), ALFA(2), ALFA(3), ALFA(4), ALFA(5), ALFA(6), 138 1 ALFA(7), ALFA(8), ALFA(9), ALFA(10), ALFA(11), ALFA(12), 139 2 ALFA(13), ALFA(14), ALFA(15), ALFA(16), ALFA(17), ALFA(18), 140 3 ALFA(19), ALFA(20), ALFA(21), ALFA(22)/ 141 4 -4.44444444444444444D-03, -9.22077922077922078D-04, 142 5 -8.84892884892884893D-05, 1.65927687832449737D-04, 143 6 2.46691372741792910D-04, 2.65995589346254780D-04, 144 7 2.61824297061500945D-04, 2.48730437344655609D-04, 145 8 2.32721040083232098D-04, 2.16362485712365082D-04, 146 9 2.00738858762752355D-04, 1.86267636637545172D-04, 147 A 1.73060775917876493D-04, 1.61091705929015752D-04, 148 B 1.50274774160908134D-04, 1.40503497391269794D-04, 149 C 1.31668816545922806D-04, 1.23667445598253261D-04, 150 D 1.16405271474737902D-04, 1.09798298372713369D-04, 151 E 1.03772410422992823D-04, 9.82626078369363448D-05/ 152 DATA ALFA(23), ALFA(24), ALFA(25), ALFA(26), ALFA(27), ALFA(28), 153 1 ALFA(29), ALFA(30), ALFA(31), ALFA(32), ALFA(33), ALFA(34), 154 2 ALFA(35), ALFA(36), ALFA(37), ALFA(38), ALFA(39), ALFA(40), 155 3 ALFA(41), ALFA(42), ALFA(43), ALFA(44)/ 156 4 9.32120517249503256D-05, 8.85710852478711718D-05, 157 5 8.42963105715700223D-05, 8.03497548407791151D-05, 158 6 7.66981345359207388D-05, 7.33122157481777809D-05, 159 7 7.01662625163141333D-05, 6.72375633790160292D-05, 160 8 6.93735541354588974D-04, 2.32241745182921654D-04, 161 9 -1.41986273556691197D-05, -1.16444931672048640D-04, 162 A -1.50803558053048762D-04, -1.55121924918096223D-04, 163 B -1.46809756646465549D-04, -1.33815503867491367D-04, 164 C -1.19744975684254051D-04, -1.06184319207974020D-04, 165 D -9.37699549891194492D-05, -8.26923045588193274D-05, 166 E -7.29374348155221211D-05, -6.44042357721016283D-05/ 167 DATA ALFA(45), ALFA(46), ALFA(47), ALFA(48), ALFA(49), ALFA(50), 168 1 ALFA(51), ALFA(52), ALFA(53), ALFA(54), ALFA(55), ALFA(56), 169 2 ALFA(57), ALFA(58), ALFA(59), ALFA(60), ALFA(61), ALFA(62), 170 3 ALFA(63), ALFA(64), ALFA(65), ALFA(66)/ 171 4 -5.69611566009369048D-05, -5.04731044303561628D-05, 172 5 -4.48134868008882786D-05, -3.98688727717598864D-05, 173 6 -3.55400532972042498D-05, -3.17414256609022480D-05, 174 7 -2.83996793904174811D-05, -2.54522720634870566D-05, 175 8 -2.28459297164724555D-05, -2.05352753106480604D-05, 176 9 -1.84816217627666085D-05, -1.66519330021393806D-05, 177 A -1.50179412980119482D-05, -1.35554031379040526D-05, 178 B -1.22434746473858131D-05, -1.10641884811308169D-05, 179 C -3.54211971457743841D-04, -1.56161263945159416D-04, 180 D 3.04465503594936410D-05, 1.30198655773242693D-04, 181 E 1.67471106699712269D-04, 1.70222587683592569D-04/ 182 DATA ALFA(67), ALFA(68), ALFA(69), ALFA(70), ALFA(71), ALFA(72), 183 1 ALFA(73), ALFA(74), ALFA(75), ALFA(76), ALFA(77), ALFA(78), 184 2 ALFA(79), ALFA(80), ALFA(81), ALFA(82), ALFA(83), ALFA(84), 185 3 ALFA(85), ALFA(86), ALFA(87), ALFA(88)/ 186 4 1.56501427608594704D-04, 1.36339170977445120D-04, 187 5 1.14886692029825128D-04, 9.45869093034688111D-05, 188 6 7.64498419250898258D-05, 6.07570334965197354D-05, 189 7 4.74394299290508799D-05, 3.62757512005344297D-05, 190 8 2.69939714979224901D-05, 1.93210938247939253D-05, 191 9 1.30056674793963203D-05, 7.82620866744496661D-06, 192 A 3.59257485819351583D-06, 1.44040049814251817D-07, 193 B -2.65396769697939116D-06, -4.91346867098485910D-06, 194 C -6.72739296091248287D-06, -8.17269379678657923D-06, 195 D -9.31304715093561232D-06, -1.02011418798016441D-05, 196 E -1.08805962510592880D-05, -1.13875481509603555D-05/ 197 DATA ALFA(89), ALFA(90), ALFA(91), ALFA(92), ALFA(93), ALFA(94), 198 1 ALFA(95), ALFA(96), ALFA(97), ALFA(98), ALFA(99), ALFA(100), 199 2 ALFA(101), ALFA(102), ALFA(103), ALFA(104), ALFA(105), 200 3 ALFA(106), ALFA(107), ALFA(108), ALFA(109), ALFA(110)/ 201 4 -1.17519675674556414D-05, -1.19987364870944141D-05, 202 5 3.78194199201772914D-04, 2.02471952761816167D-04, 203 6 -6.37938506318862408D-05, -2.38598230603005903D-04, 204 7 -3.10916256027361568D-04, -3.13680115247576316D-04, 205 8 -2.78950273791323387D-04, -2.28564082619141374D-04, 206 9 -1.75245280340846749D-04, -1.25544063060690348D-04, 207 A -8.22982872820208365D-05, -4.62860730588116458D-05, 208 B -1.72334302366962267D-05, 5.60690482304602267D-06, 209 C 2.31395443148286800D-05, 3.62642745856793957D-05, 210 D 4.58006124490188752D-05, 5.24595294959114050D-05, 211 E 5.68396208545815266D-05, 5.94349820393104052D-05/ 212 DATA ALFA(111), ALFA(112), ALFA(113), ALFA(114), ALFA(115), 213 1 ALFA(116), ALFA(117), ALFA(118), ALFA(119), ALFA(120), 214 2 ALFA(121), ALFA(122), ALFA(123), ALFA(124), ALFA(125), 215 3 ALFA(126), ALFA(127), ALFA(128), ALFA(129), ALFA(130)/ 216 4 6.06478527578421742D-05, 6.08023907788436497D-05, 217 5 6.01577894539460388D-05, 5.89199657344698500D-05, 218 6 5.72515823777593053D-05, 5.52804375585852577D-05, 219 7 5.31063773802880170D-05, 5.08069302012325706D-05, 220 8 4.84418647620094842D-05, 4.60568581607475370D-05, 221 9 -6.91141397288294174D-04, -4.29976633058871912D-04, 222 A 1.83067735980039018D-04, 6.60088147542014144D-04, 223 B 8.75964969951185931D-04, 8.77335235958235514D-04, 224 C 7.49369585378990637D-04, 5.63832329756980918D-04, 225 D 3.68059319971443156D-04, 1.88464535514455599D-04/ 226 DATA ALFA(131), ALFA(132), ALFA(133), ALFA(134), ALFA(135), 227 1 ALFA(136), ALFA(137), ALFA(138), ALFA(139), ALFA(140), 228 2 ALFA(141), ALFA(142), ALFA(143), ALFA(144), ALFA(145), 229 3 ALFA(146), ALFA(147), ALFA(148), ALFA(149), ALFA(150)/ 230 4 3.70663057664904149D-05, -8.28520220232137023D-05, 231 5 -1.72751952869172998D-04, -2.36314873605872983D-04, 232 6 -2.77966150694906658D-04, -3.02079514155456919D-04, 233 7 -3.12594712643820127D-04, -3.12872558758067163D-04, 234 8 -3.05678038466324377D-04, -2.93226470614557331D-04, 235 9 -2.77255655582934777D-04, -2.59103928467031709D-04, 236 A -2.39784014396480342D-04, -2.20048260045422848D-04, 237 B -2.00443911094971498D-04, -1.81358692210970687D-04, 238 C -1.63057674478657464D-04, -1.45712672175205844D-04, 239 D -1.29425421983924587D-04, -1.14245691942445952D-04/ 240 DATA ALFA(151), ALFA(152), ALFA(153), ALFA(154), ALFA(155), 241 1 ALFA(156), ALFA(157), ALFA(158), ALFA(159), ALFA(160), 242 2 ALFA(161), ALFA(162), ALFA(163), ALFA(164), ALFA(165), 243 3 ALFA(166), ALFA(167), ALFA(168), ALFA(169), ALFA(170)/ 244 4 1.92821964248775885D-03, 1.35592576302022234D-03, 245 5 -7.17858090421302995D-04, -2.58084802575270346D-03, 246 6 -3.49271130826168475D-03, -3.46986299340960628D-03, 247 7 -2.82285233351310182D-03, -1.88103076404891354D-03, 248 8 -8.89531718383947600D-04, 3.87912102631035228D-06, 249 9 7.28688540119691412D-04, 1.26566373053457758D-03, 250 A 1.62518158372674427D-03, 1.83203153216373172D-03, 251 B 1.91588388990527909D-03, 1.90588846755546138D-03, 252 C 1.82798982421825727D-03, 1.70389506421121530D-03, 253 D 1.55097127171097686D-03, 1.38261421852276159D-03/ 254 DATA ALFA(171), ALFA(172), ALFA(173), ALFA(174), ALFA(175), 255 1 ALFA(176), ALFA(177), ALFA(178), ALFA(179), ALFA(180)/ 256 2 1.20881424230064774D-03, 1.03676532638344962D-03, 257 3 8.71437918068619115D-04, 7.16080155297701002D-04, 258 4 5.72637002558129372D-04, 4.42089819465802277D-04, 259 5 3.24724948503090564D-04, 2.20342042730246599D-04, 260 6 1.28412898401353882D-04, 4.82005924552095464D-05/ 261 DATA BETA(1), BETA(2), BETA(3), BETA(4), BETA(5), BETA(6), 262 1 BETA(7), BETA(8), BETA(9), BETA(10), BETA(11), BETA(12), 263 2 BETA(13), BETA(14), BETA(15), BETA(16), BETA(17), BETA(18), 264 3 BETA(19), BETA(20), BETA(21), BETA(22)/ 265 4 1.79988721413553309D-02, 5.59964911064388073D-03, 266 5 2.88501402231132779D-03, 1.80096606761053941D-03, 267 6 1.24753110589199202D-03, 9.22878876572938311D-04, 268 7 7.14430421727287357D-04, 5.71787281789704872D-04, 269 8 4.69431007606481533D-04, 3.93232835462916638D-04, 270 9 3.34818889318297664D-04, 2.88952148495751517D-04, 271 A 2.52211615549573284D-04, 2.22280580798883327D-04, 272 B 1.97541838033062524D-04, 1.76836855019718004D-04, 273 C 1.59316899661821081D-04, 1.44347930197333986D-04, 274 D 1.31448068119965379D-04, 1.20245444949302884D-04, 275 E 1.10449144504599392D-04, 1.01828770740567258D-04/ 276 DATA BETA(23), BETA(24), BETA(25), BETA(26), BETA(27), BETA(28), 277 1 BETA(29), BETA(30), BETA(31), BETA(32), BETA(33), BETA(34), 278 2 BETA(35), BETA(36), BETA(37), BETA(38), BETA(39), BETA(40), 279 3 BETA(41), BETA(42), BETA(43), BETA(44)/ 280 4 9.41998224204237509D-05, 8.74130545753834437D-05, 281 5 8.13466262162801467D-05, 7.59002269646219339D-05, 282 6 7.09906300634153481D-05, 6.65482874842468183D-05, 283 7 6.25146958969275078D-05, 5.88403394426251749D-05, 284 8 -1.49282953213429172D-03, -8.78204709546389328D-04, 285 9 -5.02916549572034614D-04, -2.94822138512746025D-04, 286 A -1.75463996970782828D-04, -1.04008550460816434D-04, 287 B -5.96141953046457895D-05, -3.12038929076098340D-05, 288 C -1.26089735980230047D-05, -2.42892608575730389D-07, 289 D 8.05996165414273571D-06, 1.36507009262147391D-05, 290 E 1.73964125472926261D-05, 1.98672978842133780D-05/ 291 DATA BETA(45), BETA(46), BETA(47), BETA(48), BETA(49), BETA(50), 292 1 BETA(51), BETA(52), BETA(53), BETA(54), BETA(55), BETA(56), 293 2 BETA(57), BETA(58), BETA(59), BETA(60), BETA(61), BETA(62), 294 3 BETA(63), BETA(64), BETA(65), BETA(66)/ 295 4 2.14463263790822639D-05, 2.23954659232456514D-05, 296 5 2.28967783814712629D-05, 2.30785389811177817D-05, 297 6 2.30321976080909144D-05, 2.28236073720348722D-05, 298 7 2.25005881105292418D-05, 2.20981015361991429D-05, 299 8 2.16418427448103905D-05, 2.11507649256220843D-05, 300 9 2.06388749782170737D-05, 2.01165241997081666D-05, 301 A 1.95913450141179244D-05, 1.90689367910436740D-05, 302 B 1.85533719641636667D-05, 1.80475722259674218D-05, 303 C 5.52213076721292790D-04, 4.47932581552384646D-04, 304 D 2.79520653992020589D-04, 1.52468156198446602D-04, 305 E 6.93271105657043598D-05, 1.76258683069991397D-05/ 306 DATA BETA(67), BETA(68), BETA(69), BETA(70), BETA(71), BETA(72), 307 1 BETA(73), BETA(74), BETA(75), BETA(76), BETA(77), BETA(78), 308 2 BETA(79), BETA(80), BETA(81), BETA(82), BETA(83), BETA(84), 309 3 BETA(85), BETA(86), BETA(87), BETA(88)/ 310 4 -1.35744996343269136D-05, -3.17972413350427135D-05, 311 5 -4.18861861696693365D-05, -4.69004889379141029D-05, 312 6 -4.87665447413787352D-05, -4.87010031186735069D-05, 313 7 -4.74755620890086638D-05, -4.55813058138628452D-05, 314 8 -4.33309644511266036D-05, -4.09230193157750364D-05, 315 9 -3.84822638603221274D-05, -3.60857167535410501D-05, 316 A -3.37793306123367417D-05, -3.15888560772109621D-05, 317 B -2.95269561750807315D-05, -2.75978914828335759D-05, 318 C -2.58006174666883713D-05, -2.41308356761280200D-05, 319 D -2.25823509518346033D-05, -2.11479656768912971D-05, 320 E -1.98200638885294927D-05, -1.85909870801065077D-05/ 321 DATA BETA(89), BETA(90), BETA(91), BETA(92), BETA(93), BETA(94), 322 1 BETA(95), BETA(96), BETA(97), BETA(98), BETA(99), BETA(100), 323 2 BETA(101), BETA(102), BETA(103), BETA(104), BETA(105), 324 3 BETA(106), BETA(107), BETA(108), BETA(109), BETA(110)/ 325 4 -1.74532699844210224D-05, -1.63997823854497997D-05, 326 5 -4.74617796559959808D-04, -4.77864567147321487D-04, 327 6 -3.20390228067037603D-04, -1.61105016119962282D-04, 328 7 -4.25778101285435204D-05, 3.44571294294967503D-05, 329 8 7.97092684075674924D-05, 1.03138236708272200D-04, 330 9 1.12466775262204158D-04, 1.13103642108481389D-04, 331 A 1.08651634848774268D-04, 1.01437951597661973D-04, 332 B 9.29298396593363896D-05, 8.40293133016089978D-05, 333 C 7.52727991349134062D-05, 6.69632521975730872D-05, 334 D 5.92564547323194704D-05, 5.22169308826975567D-05, 335 E 4.58539485165360646D-05, 4.01445513891486808D-05/ 336 DATA BETA(111), BETA(112), BETA(113), BETA(114), BETA(115), 337 1 BETA(116), BETA(117), BETA(118), BETA(119), BETA(120), 338 2 BETA(121), BETA(122), BETA(123), BETA(124), BETA(125), 339 3 BETA(126), BETA(127), BETA(128), BETA(129), BETA(130)/ 340 4 3.50481730031328081D-05, 3.05157995034346659D-05, 341 5 2.64956119950516039D-05, 2.29363633690998152D-05, 342 6 1.97893056664021636D-05, 1.70091984636412623D-05, 343 7 1.45547428261524004D-05, 1.23886640995878413D-05, 344 8 1.04775876076583236D-05, 8.79179954978479373D-06, 345 9 7.36465810572578444D-04, 8.72790805146193976D-04, 346 A 6.22614862573135066D-04, 2.85998154194304147D-04, 347 B 3.84737672879366102D-06, -1.87906003636971558D-04, 348 C -2.97603646594554535D-04, -3.45998126832656348D-04, 349 D -3.53382470916037712D-04, -3.35715635775048757D-04/ 350 DATA BETA(131), BETA(132), BETA(133), BETA(134), BETA(135), 351 1 BETA(136), BETA(137), BETA(138), BETA(139), BETA(140), 352 2 BETA(141), BETA(142), BETA(143), BETA(144), BETA(145), 353 3 BETA(146), BETA(147), BETA(148), BETA(149), BETA(150)/ 354 4 -3.04321124789039809D-04, -2.66722723047612821D-04, 355 5 -2.27654214122819527D-04, -1.89922611854562356D-04, 356 6 -1.55058918599093870D-04, -1.23778240761873630D-04, 357 7 -9.62926147717644187D-05, -7.25178327714425337D-05, 358 8 -5.22070028895633801D-05, -3.50347750511900522D-05, 359 9 -2.06489761035551757D-05, -8.70106096849767054D-06, 360 A 1.13698686675100290D-06, 9.16426474122778849D-06, 361 B 1.56477785428872620D-05, 2.08223629482466847D-05, 362 C 2.48923381004595156D-05, 2.80340509574146325D-05, 363 D 3.03987774629861915D-05, 3.21156731406700616D-05/ 364 DATA BETA(151), BETA(152), BETA(153), BETA(154), BETA(155), 365 1 BETA(156), BETA(157), BETA(158), BETA(159), BETA(160), 366 2 BETA(161), BETA(162), BETA(163), BETA(164), BETA(165), 367 3 BETA(166), BETA(167), BETA(168), BETA(169), BETA(170)/ 368 4 -1.80182191963885708D-03, -2.43402962938042533D-03, 369 5 -1.83422663549856802D-03, -7.62204596354009765D-04, 370 6 2.39079475256927218D-04, 9.49266117176881141D-04, 371 7 1.34467449701540359D-03, 1.48457495259449178D-03, 372 8 1.44732339830617591D-03, 1.30268261285657186D-03, 373 9 1.10351597375642682D-03, 8.86047440419791759D-04, 374 A 6.73073208165665473D-04, 4.77603872856582378D-04, 375 B 3.05991926358789362D-04, 1.60315694594721630D-04, 376 C 4.00749555270613286D-05, -5.66607461635251611D-05, 377 D -1.32506186772982638D-04, -1.90296187989614057D-04/ 378 DATA BETA(171), BETA(172), BETA(173), BETA(174), BETA(175), 379 1 BETA(176), BETA(177), BETA(178), BETA(179), BETA(180), 380 2 BETA(181), BETA(182), BETA(183), BETA(184), BETA(185), 381 3 BETA(186), BETA(187), BETA(188), BETA(189), BETA(190)/ 382 4 -2.32811450376937408D-04, -2.62628811464668841D-04, 383 5 -2.82050469867598672D-04, -2.93081563192861167D-04, 384 6 -2.97435962176316616D-04, -2.96557334239348078D-04, 385 7 -2.91647363312090861D-04, -2.83696203837734166D-04, 386 8 -2.73512317095673346D-04, -2.61750155806768580D-04, 387 9 6.38585891212050914D-03, 9.62374215806377941D-03, 388 A 7.61878061207001043D-03, 2.83219055545628054D-03, 389 B -2.09841352012720090D-03, -5.73826764216626498D-03, 390 C -7.70804244495414620D-03, -8.21011692264844401D-03, 391 D -7.65824520346905413D-03, -6.47209729391045177D-03/ 392 DATA BETA(191), BETA(192), BETA(193), BETA(194), BETA(195), 393 1 BETA(196), BETA(197), BETA(198), BETA(199), BETA(200), 394 2 BETA(201), BETA(202), BETA(203), BETA(204), BETA(205), 395 3 BETA(206), BETA(207), BETA(208), BETA(209), BETA(210)/ 396 4 -4.99132412004966473D-03, -3.45612289713133280D-03, 397 5 -2.01785580014170775D-03, -7.59430686781961401D-04, 398 6 2.84173631523859138D-04, 1.10891667586337403D-03, 399 7 1.72901493872728771D-03, 2.16812590802684701D-03, 400 8 2.45357710494539735D-03, 2.61281821058334862D-03, 401 9 2.67141039656276912D-03, 2.65203073395980430D-03, 402 A 2.57411652877287315D-03, 2.45389126236094427D-03, 403 B 2.30460058071795494D-03, 2.13684837686712662D-03, 404 C 1.95896528478870911D-03, 1.77737008679454412D-03, 405 D 1.59690280765839059D-03, 1.42111975664438546D-03/ 406 DATA GAMA(1), GAMA(2), GAMA(3), GAMA(4), GAMA(5), GAMA(6), 407 1 GAMA(7), GAMA(8), GAMA(9), GAMA(10), GAMA(11), GAMA(12), 408 2 GAMA(13), GAMA(14), GAMA(15), GAMA(16), GAMA(17), GAMA(18), 409 3 GAMA(19), GAMA(20), GAMA(21), GAMA(22)/ 410 4 6.29960524947436582D-01, 2.51984209978974633D-01, 411 5 1.54790300415655846D-01, 1.10713062416159013D-01, 412 6 8.57309395527394825D-02, 6.97161316958684292D-02, 413 7 5.86085671893713576D-02, 5.04698873536310685D-02, 414 8 4.42600580689154809D-02, 3.93720661543509966D-02, 415 9 3.54283195924455368D-02, 3.21818857502098231D-02, 416 A 2.94646240791157679D-02, 2.71581677112934479D-02, 417 B 2.51768272973861779D-02, 2.34570755306078891D-02, 418 C 2.19508390134907203D-02, 2.06210828235646240D-02, 419 D 1.94388240897880846D-02, 1.83810633800683158D-02, 420 E 1.74293213231963172D-02, 1.65685837786612353D-02/ 421 DATA GAMA(23), GAMA(24), GAMA(25), GAMA(26), GAMA(27), GAMA(28), 422 1 GAMA(29), GAMA(30)/ 423 2 1.57865285987918445D-02, 1.50729501494095594D-02, 424 3 1.44193250839954639D-02, 1.38184805735341786D-02, 425 4 1.32643378994276568D-02, 1.27517121970498651D-02, 426 5 1.22761545318762767D-02, 1.18338262398482403D-02/ 427 DATA EX1, EX2, HPI, GPI, THPI / 428 1 3.33333333333333333D-01, 6.66666666666666667D-01, 429 2 1.57079632679489662D+00, 3.14159265358979324D+00, 430 3 4.71238898038468986D+00/ 431 DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / 432C 433 RFNU = 1.0D0/FNU 434C----------------------------------------------------------------------- 435C OVERFLOW TEST (Z/FNU TOO SMALL) 436C----------------------------------------------------------------------- 437 TEST = D1MACH(1)*1.0D+3 438 AC = FNU*TEST 439 IF (DABS(ZR).GT.AC .OR. DABS(ZI).GT.AC) GO TO 15 440 ZETA1R = 2.0D0*DABS(DLOG(TEST))+FNU 441 ZETA1I = 0.0D0 442 ZETA2R = FNU 443 ZETA2I = 0.0D0 444 PHIR = 1.0D0 445 PHII = 0.0D0 446 ARGR = 1.0D0 447 ARGI = 0.0D0 448 RETURN 449 15 CONTINUE 450 ZBR = ZR*RFNU 451 ZBI = ZI*RFNU 452 RFNU2 = RFNU*RFNU 453C----------------------------------------------------------------------- 454C COMPUTE IN THE FOURTH QUADRANT 455C----------------------------------------------------------------------- 456 FN13 = FNU**EX1 457 FN23 = FN13*FN13 458 RFN13 = 1.0D0/FN13 459 W2R = CONER - ZBR*ZBR + ZBI*ZBI 460 W2I = CONEI - ZBR*ZBI - ZBR*ZBI 461 AW2 = ZABS(W2R,W2I) 462 IF (AW2.GT.0.25D0) GO TO 130 463C----------------------------------------------------------------------- 464C POWER SERIES FOR CABS(W2).LE.0.25D0 465C----------------------------------------------------------------------- 466 K = 1 467 PR(1) = CONER 468 PI(1) = CONEI 469 SUMAR = GAMA(1) 470 SUMAI = ZEROI 471 AP(1) = 1.0D0 472 IF (AW2.LT.TOL) GO TO 20 473 DO 10 K=2,30 474 PR(K) = PR(K-1)*W2R - PI(K-1)*W2I 475 PI(K) = PR(K-1)*W2I + PI(K-1)*W2R 476 SUMAR = SUMAR + PR(K)*GAMA(K) 477 SUMAI = SUMAI + PI(K)*GAMA(K) 478 AP(K) = AP(K-1)*AW2 479 IF (AP(K).LT.TOL) GO TO 20 480 10 CONTINUE 481 K = 30 482 20 CONTINUE 483 KMAX = K 484 ZETAR = W2R*SUMAR - W2I*SUMAI 485 ZETAI = W2R*SUMAI + W2I*SUMAR 486 ARGR = ZETAR*FN23 487 ARGI = ZETAI*FN23 488 CALL ZSQRTAMOS(SUMAR, SUMAI, ZAR, ZAI) 489 CALL ZSQRTAMOS(W2R, W2I, STR, STI) 490 ZETA2R = STR*FNU 491 ZETA2I = STI*FNU 492 STR = CONER + EX2*(ZETAR*ZAR-ZETAI*ZAI) 493 STI = CONEI + EX2*(ZETAR*ZAI+ZETAI*ZAR) 494 ZETA1R = STR*ZETA2R - STI*ZETA2I 495 ZETA1I = STR*ZETA2I + STI*ZETA2R 496 ZAR = ZAR + ZAR 497 ZAI = ZAI + ZAI 498 CALL ZSQRTAMOS(ZAR, ZAI, STR, STI) 499 PHIR = STR*RFN13 500 PHII = STI*RFN13 501 IF (IPMTR.EQ.1) GO TO 120 502C----------------------------------------------------------------------- 503C SUM SERIES FOR ASUM AND BSUM 504C----------------------------------------------------------------------- 505 SUMBR = ZEROR 506 SUMBI = ZEROI 507 DO 30 K=1,KMAX 508 SUMBR = SUMBR + PR(K)*BETA(K) 509 SUMBI = SUMBI + PI(K)*BETA(K) 510 30 CONTINUE 511 ASUMR = ZEROR 512 ASUMI = ZEROI 513 BSUMR = SUMBR 514 BSUMI = SUMBI 515 L1 = 0 516 L2 = 30 517 BTOL = TOL*(DABS(BSUMR)+DABS(BSUMI)) 518 ATOL = TOL 519 PP = 1.0D0 520 IAS = 0 521 IBS = 0 522 IF (RFNU2.LT.TOL) GO TO 110 523 DO 100 IS=2,7 524 ATOL = ATOL/RFNU2 525 PP = PP*RFNU2 526 IF (IAS.EQ.1) GO TO 60 527 SUMAR = ZEROR 528 SUMAI = ZEROI 529 DO 40 K=1,KMAX 530 M = L1 + K 531 SUMAR = SUMAR + PR(K)*ALFA(M) 532 SUMAI = SUMAI + PI(K)*ALFA(M) 533 IF (AP(K).LT.ATOL) GO TO 50 534 40 CONTINUE 535 50 CONTINUE 536 ASUMR = ASUMR + SUMAR*PP 537 ASUMI = ASUMI + SUMAI*PP 538 IF (PP.LT.TOL) IAS = 1 539 60 CONTINUE 540 IF (IBS.EQ.1) GO TO 90 541 SUMBR = ZEROR 542 SUMBI = ZEROI 543 DO 70 K=1,KMAX 544 M = L2 + K 545 SUMBR = SUMBR + PR(K)*BETA(M) 546 SUMBI = SUMBI + PI(K)*BETA(M) 547 IF (AP(K).LT.ATOL) GO TO 80 548 70 CONTINUE 549 80 CONTINUE 550 BSUMR = BSUMR + SUMBR*PP 551 BSUMI = BSUMI + SUMBI*PP 552 IF (PP.LT.BTOL) IBS = 1 553 90 CONTINUE 554 IF (IAS.EQ.1 .AND. IBS.EQ.1) GO TO 110 555 L1 = L1 + 30 556 L2 = L2 + 30 557 100 CONTINUE 558 110 CONTINUE 559 ASUMR = ASUMR + CONER 560 PP = RFNU*RFN13 561 BSUMR = BSUMR*PP 562 BSUMI = BSUMI*PP 563 120 CONTINUE 564 RETURN 565C----------------------------------------------------------------------- 566C CABS(W2).GT.0.25D0 567C----------------------------------------------------------------------- 568 130 CONTINUE 569 CALL ZSQRTAMOS(W2R, W2I, WR, WI) 570 IF (WR.LT.0.0D0) WR = 0.0D0 571 IF (WI.LT.0.0D0) WI = 0.0D0 572 STR = CONER + WR 573 STI = WI 574 CALL ZDIV(STR, STI, ZBR, ZBI, ZAR, ZAI) 575 CALL ZLOG(ZAR, ZAI, ZCR, ZCI, IDUM) 576 IF (ZCI.LT.0.0D0) ZCI = 0.0D0 577 IF (ZCI.GT.HPI) ZCI = HPI 578 IF (ZCR.LT.0.0D0) ZCR = 0.0D0 579 ZTHR = (ZCR-WR)*1.5D0 580 ZTHI = (ZCI-WI)*1.5D0 581 ZETA1R = ZCR*FNU 582 ZETA1I = ZCI*FNU 583 ZETA2R = WR*FNU 584 ZETA2I = WI*FNU 585 AZTH = ZABS(ZTHR,ZTHI) 586 ANG = THPI 587 IF (ZTHR.GE.0.0D0 .AND. ZTHI.LT.0.0D0) GO TO 140 588 ANG = HPI 589 IF (ZTHR.EQ.0.0D0) GO TO 140 590 ANG = DATAN(ZTHI/ZTHR) 591 IF (ZTHR.LT.0.0D0) ANG = ANG + GPI 592 140 CONTINUE 593 PP = AZTH**EX2 594 ANG = ANG*EX2 595 ZETAR = PP*DCOS(ANG) 596 ZETAI = PP*DSIN(ANG) 597 IF (ZETAI.LT.0.0D0) ZETAI = 0.0D0 598 ARGR = ZETAR*FN23 599 ARGI = ZETAI*FN23 600 CALL ZDIV(ZTHR, ZTHI, ZETAR, ZETAI, RTZTR, RTZTI) 601 CALL ZDIV(RTZTR, RTZTI, WR, WI, ZAR, ZAI) 602 TZAR = ZAR + ZAR 603 TZAI = ZAI + ZAI 604 CALL ZSQRTAMOS(TZAR, TZAI, STR, STI) 605 PHIR = STR*RFN13 606 PHII = STI*RFN13 607 IF (IPMTR.EQ.1) GO TO 120 608 RAW = 1.0D0/DSQRT(AW2) 609 STR = WR*RAW 610 STI = -WI*RAW 611 TFNR = STR*RFNU*RAW 612 TFNI = STI*RFNU*RAW 613 RAZTH = 1.0D0/AZTH 614 STR = ZTHR*RAZTH 615 STI = -ZTHI*RAZTH 616 RZTHR = STR*RAZTH*RFNU 617 RZTHI = STI*RAZTH*RFNU 618 ZCR = RZTHR*AR(2) 619 ZCI = RZTHI*AR(2) 620 RAW2 = 1.0D0/AW2 621 STR = W2R*RAW2 622 STI = -W2I*RAW2 623 T2R = STR*RAW2 624 T2I = STI*RAW2 625 STR = T2R*C(2) + C(3) 626 STI = T2I*C(2) 627 UPR(2) = STR*TFNR - STI*TFNI 628 UPI(2) = STR*TFNI + STI*TFNR 629 BSUMR = UPR(2) + ZCR 630 BSUMI = UPI(2) + ZCI 631 ASUMR = ZEROR 632 ASUMI = ZEROI 633 IF (RFNU.LT.TOL) GO TO 220 634 PRZTHR = RZTHR 635 PRZTHI = RZTHI 636 PTFNR = TFNR 637 PTFNI = TFNI 638 UPR(1) = CONER 639 UPI(1) = CONEI 640 PP = 1.0D0 641 BTOL = TOL*(DABS(BSUMR)+DABS(BSUMI)) 642 KS = 0 643 KP1 = 2 644 L = 3 645 IAS = 0 646 IBS = 0 647 DO 210 LR=2,12,2 648 LRP1 = LR + 1 649C----------------------------------------------------------------------- 650C COMPUTE TWO ADDITIONAL CR, DR, AND UP FOR TWO MORE TERMS IN 651C NEXT SUMA AND SUMB 652C----------------------------------------------------------------------- 653 DO 160 K=LR,LRP1 654 KS = KS + 1 655 KP1 = KP1 + 1 656 L = L + 1 657 ZAR = C(L) 658 ZAI = ZEROI 659 DO 150 J=2,KP1 660 L = L + 1 661 STR = ZAR*T2R - T2I*ZAI + C(L) 662 ZAI = ZAR*T2I + ZAI*T2R 663 ZAR = STR 664 150 CONTINUE 665 STR = PTFNR*TFNR - PTFNI*TFNI 666 PTFNI = PTFNR*TFNI + PTFNI*TFNR 667 PTFNR = STR 668 UPR(KP1) = PTFNR*ZAR - PTFNI*ZAI 669 UPI(KP1) = PTFNI*ZAR + PTFNR*ZAI 670 CRR(KS) = PRZTHR*BR(KS+1) 671 CRI(KS) = PRZTHI*BR(KS+1) 672 STR = PRZTHR*RZTHR - PRZTHI*RZTHI 673 PRZTHI = PRZTHR*RZTHI + PRZTHI*RZTHR 674 PRZTHR = STR 675 DRR(KS) = PRZTHR*AR(KS+2) 676 DRI(KS) = PRZTHI*AR(KS+2) 677 160 CONTINUE 678 PP = PP*RFNU2 679 IF (IAS.EQ.1) GO TO 180 680 SUMAR = UPR(LRP1) 681 SUMAI = UPI(LRP1) 682 JU = LRP1 683 DO 170 JR=1,LR 684 JU = JU - 1 685 SUMAR = SUMAR + CRR(JR)*UPR(JU) - CRI(JR)*UPI(JU) 686 SUMAI = SUMAI + CRR(JR)*UPI(JU) + CRI(JR)*UPR(JU) 687 170 CONTINUE 688 ASUMR = ASUMR + SUMAR 689 ASUMI = ASUMI + SUMAI 690 TEST = DABS(SUMAR) + DABS(SUMAI) 691 IF (PP.LT.TOL .AND. TEST.LT.TOL) IAS = 1 692 180 CONTINUE 693 IF (IBS.EQ.1) GO TO 200 694 SUMBR = UPR(LR+2) + UPR(LRP1)*ZCR - UPI(LRP1)*ZCI 695 SUMBI = UPI(LR+2) + UPR(LRP1)*ZCI + UPI(LRP1)*ZCR 696 JU = LRP1 697 DO 190 JR=1,LR 698 JU = JU - 1 699 SUMBR = SUMBR + DRR(JR)*UPR(JU) - DRI(JR)*UPI(JU) 700 SUMBI = SUMBI + DRR(JR)*UPI(JU) + DRI(JR)*UPR(JU) 701 190 CONTINUE 702 BSUMR = BSUMR + SUMBR 703 BSUMI = BSUMI + SUMBI 704 TEST = DABS(SUMBR) + DABS(SUMBI) 705 IF (PP.LT.BTOL .AND. TEST.LT.BTOL) IBS = 1 706 200 CONTINUE 707 IF (IAS.EQ.1 .AND. IBS.EQ.1) GO TO 220 708 210 CONTINUE 709 220 CONTINUE 710 ASUMR = ASUMR + CONER 711 STR = -BSUMR*RFN13 712 STI = -BSUMI*RFN13 713 CALL ZDIV(STR, STI, RTZTR, RTZTI, BSUMR, BSUMI) 714 GO TO 120 715 END 716