GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
zunhj.f
Go to the documentation of this file.
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 XZABS,ZDIV,XZLOG,XZSQRT,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, xzabs, 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 = xzabs(w2r,w2i)
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 xzsqrt(sumar, sumai, zar, zai)
488  CALL xzsqrt(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 xzsqrt(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 xzsqrt(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 xzlog(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 = xzabs(zthr,zthi)
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 xzsqrt(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
d1mach
double precision function d1mach(i)
Definition: d1mach.f:23
xzabs
double precision function xzabs(ZR, ZI)
Definition: xzabs.f:2
xzlog
subroutine xzlog(AR, AI, BR, BI, IERR)
Definition: xzlog.f:2
xzsqrt
subroutine xzsqrt(AR, AI, BR, BI)
Definition: xzsqrt.f:2
zdiv
subroutine zdiv(AR, AI, BR, BI, CR, CI)
Definition: zdiv.f:2
zunhj
subroutine zunhj(ZR, ZI, FNU, IPMTR, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
Definition: zunhj.f:3