GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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
Generated on Mon Dec 30 2013 03:04:45 for GNU Octave by   doxygen 1.8.1.2