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