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
cunhj.f
Go to the documentation of this file.
1  SUBROUTINE cunhj(Z, FNU, IPMTR, TOL, PHI, ARG, ZETA1, ZETA2,
2  * asum, bsum)
3 C***BEGIN PROLOGUE CUNHJ
4 C***REFER TO CBESI,CBESK
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 CUNHJ 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 (NONE)
33 C***END PROLOGUE CUNHJ
34  COMPLEX arg, asum, bsum, cfnu, cone, cr, czero, dr, p, phi,
35  * przth, ptfn, rfn13, rtzta, rzth, suma, sumb, tfn, t2, up, w, w2,
36  * z, za, zb, zc, zeta, zeta1, zeta2, zth
37  REAL alfa, ang, ap, ar, atol, aw2, azth, beta, br, btol, c, ex1,
38  * ex2, fnu, fn13, fn23, gama, hpi, pi, pp, rfnu, rfnu2, thpi, tol,
39  * wi, wr, zci, zcr, zetai, zetar, zthi, zthr, asumr, asumi, bsumr,
40  * bsumi, test, tstr, tsti, ac
41  INTEGER ias, ibs, ipmtr, is, j, jr, ju, k, kmax, kp1, ks, l, lr,
42  * lrp1, l1, l2, m
43  dimension ar(14), br(14), c(105), alfa(180), beta(210), gama(30),
44  * ap(30), p(30), up(14), cr(14), dr(14)
45  DATA ar(1), ar(2), ar(3), ar(4), ar(5), ar(6), ar(7), ar(8),
46  1 ar(9), ar(10), ar(11), ar(12), ar(13), ar(14)/
47  2 1.00000000000000000e+00, 1.04166666666666667e-01,
48  3 8.35503472222222222e-02, 1.28226574556327160e-01,
49  4 2.91849026464140464e-01, 8.81627267443757652e-01,
50  5 3.32140828186276754e+00, 1.49957629868625547e+01,
51  6 7.89230130115865181e+01, 4.74451538868264323e+02,
52  7 3.20749009089066193e+03, 2.40865496408740049e+04,
53  8 1.98923119169509794e+05, 1.79190200777534383e+06/
54  DATA br(1), br(2), br(3), br(4), br(5), br(6), br(7), br(8),
55  1 br(9), br(10), br(11), br(12), br(13), br(14)/
56  2 1.00000000000000000e+00, -1.45833333333333333e-01,
57  3 -9.87413194444444444e-02, -1.43312053915895062e-01,
58  4 -3.17227202678413548e-01, -9.42429147957120249e-01,
59  5 -3.51120304082635426e+00, -1.57272636203680451e+01,
60  6 -8.22814390971859444e+01, -4.92355370523670524e+02,
61  7 -3.31621856854797251e+03, -2.48276742452085896e+04,
62  8 -2.04526587315129788e+05, -1.83844491706820990e+06/
63  DATA c(1), c(2), c(3), c(4), c(5), c(6), c(7), c(8), c(9), c(10),
64  1 c(11), c(12), c(13), c(14), c(15), c(16), c(17), c(18),
65  2 c(19), c(20), c(21), c(22), c(23), c(24)/
66  3 1.00000000000000000e+00, -2.08333333333333333e-01,
67  4 1.25000000000000000e-01, 3.34201388888888889e-01,
68  5 -4.01041666666666667e-01, 7.03125000000000000e-02,
69  6 -1.02581259645061728e+00, 1.84646267361111111e+00,
70  7 -8.91210937500000000e-01, 7.32421875000000000e-02,
71  8 4.66958442342624743e+00, -1.12070026162229938e+01,
72  9 8.78912353515625000e+00, -2.36408691406250000e+00,
73  a 1.12152099609375000e-01, -2.82120725582002449e+01,
74  b 8.46362176746007346e+01, -9.18182415432400174e+01,
75  c 4.25349987453884549e+01, -7.36879435947963170e+00,
76  d 2.27108001708984375e-01, 2.12570130039217123e+02,
77  e -7.65252468141181642e+02, 1.05999045252799988e+03/
78  DATA c(25), c(26), c(27), c(28), c(29), c(30), c(31), c(32),
79  1 c(33), c(34), c(35), c(36), c(37), c(38), c(39), c(40),
80  2 c(41), c(42), c(43), c(44), c(45), c(46), c(47), c(48)/
81  3 -6.99579627376132541e+02, 2.18190511744211590e+02,
82  4 -2.64914304869515555e+01, 5.72501420974731445e-01,
83  5 -1.91945766231840700e+03, 8.06172218173730938e+03,
84  6 -1.35865500064341374e+04, 1.16553933368645332e+04,
85  7 -5.30564697861340311e+03, 1.20090291321635246e+03,
86  8 -1.08090919788394656e+02, 1.72772750258445740e+00,
87  9 2.02042913309661486e+04, -9.69805983886375135e+04,
88  a 1.92547001232531532e+05, -2.03400177280415534e+05,
89  b 1.22200464983017460e+05, -4.11926549688975513e+04,
90  c 7.10951430248936372e+03, -4.93915304773088012e+02,
91  d 6.07404200127348304e+00, -2.42919187900551333e+05,
92  e 1.31176361466297720e+06, -2.99801591853810675e+06/
93  DATA c(49), c(50), c(51), c(52), c(53), c(54), c(55), c(56),
94  1 c(57), c(58), c(59), c(60), c(61), c(62), c(63), c(64),
95  2 c(65), c(66), c(67), c(68), c(69), c(70), c(71), c(72)/
96  3 3.76327129765640400e+06, -2.81356322658653411e+06,
97  4 1.26836527332162478e+06, -3.31645172484563578e+05,
98  5 4.52187689813627263e+04, -2.49983048181120962e+03,
99  6 2.43805296995560639e+01, 3.28446985307203782e+06,
100  7 -1.97068191184322269e+07, 5.09526024926646422e+07,
101  8 -7.41051482115326577e+07, 6.63445122747290267e+07,
102  9 -3.75671766607633513e+07, 1.32887671664218183e+07,
103  a -2.78561812808645469e+06, 3.08186404612662398e+05,
104  b -1.38860897537170405e+04, 1.10017140269246738e+02,
105  c -4.93292536645099620e+07, 3.25573074185765749e+08,
106  d -9.39462359681578403e+08, 1.55359689957058006e+09,
107  e -1.62108055210833708e+09, 1.10684281682301447e+09/
108  DATA c(73), c(74), c(75), c(76), c(77), c(78), c(79), c(80),
109  1 c(81), c(82), c(83), c(84), c(85), c(86), c(87), c(88),
110  2 c(89), c(90), c(91), c(92), c(93), c(94), c(95), c(96)/
111  3 -4.95889784275030309e+08, 1.42062907797533095e+08,
112  4 -2.44740627257387285e+07, 2.24376817792244943e+06,
113  5 -8.40054336030240853e+04, 5.51335896122020586e+02,
114  6 8.14789096118312115e+08, -5.86648149205184723e+09,
115  7 1.86882075092958249e+10, -3.46320433881587779e+10,
116  8 4.12801855797539740e+10, -3.30265997498007231e+10,
117  9 1.79542137311556001e+10, -6.56329379261928433e+09,
118  a 1.55927986487925751e+09, -2.25105661889415278e+08,
119  b 1.73951075539781645e+07, -5.49842327572288687e+05,
120  c 3.03809051092238427e+03, -1.46792612476956167e+10,
121  d 1.14498237732025810e+11, -3.99096175224466498e+11,
122  e 8.19218669548577329e+11, -1.09837515608122331e+12/
123  DATA c(97), c(98), c(99), c(100), c(101), c(102), c(103), c(104),
124  1 c(105)/
125  2 1.00815810686538209e+12, -6.45364869245376503e+11,
126  3 2.87900649906150589e+11, -8.78670721780232657e+10,
127  4 1.76347306068349694e+10, -2.16716498322379509e+09,
128  5 1.43157876718888981e+08, -3.87183344257261262e+06,
129  6 1.82577554742931747e+04/
130  DATA alfa(1), alfa(2), alfa(3), alfa(4), alfa(5), alfa(6),
131  1 alfa(7), alfa(8), alfa(9), alfa(10), alfa(11), alfa(12),
132  2 alfa(13), alfa(14), alfa(15), alfa(16), alfa(17), alfa(18),
133  3 alfa(19), alfa(20), alfa(21), alfa(22)/
134  4 -4.44444444444444444e-03, -9.22077922077922078e-04,
135  5 -8.84892884892884893e-05, 1.65927687832449737e-04,
136  6 2.46691372741792910e-04, 2.65995589346254780e-04,
137  7 2.61824297061500945e-04, 2.48730437344655609e-04,
138  8 2.32721040083232098e-04, 2.16362485712365082e-04,
139  9 2.00738858762752355e-04, 1.86267636637545172e-04,
140  a 1.73060775917876493e-04, 1.61091705929015752e-04,
141  b 1.50274774160908134e-04, 1.40503497391269794e-04,
142  c 1.31668816545922806e-04, 1.23667445598253261e-04,
143  d 1.16405271474737902e-04, 1.09798298372713369e-04,
144  e 1.03772410422992823e-04, 9.82626078369363448e-05/
145  DATA alfa(23), alfa(24), alfa(25), alfa(26), alfa(27), alfa(28),
146  1 alfa(29), alfa(30), alfa(31), alfa(32), alfa(33), alfa(34),
147  2 alfa(35), alfa(36), alfa(37), alfa(38), alfa(39), alfa(40),
148  3 alfa(41), alfa(42), alfa(43), alfa(44)/
149  4 9.32120517249503256e-05, 8.85710852478711718e-05,
150  5 8.42963105715700223e-05, 8.03497548407791151e-05,
151  6 7.66981345359207388e-05, 7.33122157481777809e-05,
152  7 7.01662625163141333e-05, 6.72375633790160292e-05,
153  8 6.93735541354588974e-04, 2.32241745182921654e-04,
154  9 -1.41986273556691197e-05, -1.16444931672048640e-04,
155  a -1.50803558053048762e-04, -1.55121924918096223e-04,
156  b -1.46809756646465549e-04, -1.33815503867491367e-04,
157  c -1.19744975684254051e-04, -1.06184319207974020e-04,
158  d -9.37699549891194492e-05, -8.26923045588193274e-05,
159  e -7.29374348155221211e-05, -6.44042357721016283e-05/
160  DATA alfa(45), alfa(46), alfa(47), alfa(48), alfa(49), alfa(50),
161  1 alfa(51), alfa(52), alfa(53), alfa(54), alfa(55), alfa(56),
162  2 alfa(57), alfa(58), alfa(59), alfa(60), alfa(61), alfa(62),
163  3 alfa(63), alfa(64), alfa(65), alfa(66)/
164  4 -5.69611566009369048e-05, -5.04731044303561628e-05,
165  5 -4.48134868008882786e-05, -3.98688727717598864e-05,
166  6 -3.55400532972042498e-05, -3.17414256609022480e-05,
167  7 -2.83996793904174811e-05, -2.54522720634870566e-05,
168  8 -2.28459297164724555e-05, -2.05352753106480604e-05,
169  9 -1.84816217627666085e-05, -1.66519330021393806e-05,
170  a -1.50179412980119482e-05, -1.35554031379040526e-05,
171  b -1.22434746473858131e-05, -1.10641884811308169e-05,
172  c -3.54211971457743841e-04, -1.56161263945159416e-04,
173  d 3.04465503594936410e-05, 1.30198655773242693e-04,
174  e 1.67471106699712269e-04, 1.70222587683592569e-04/
175  DATA alfa(67), alfa(68), alfa(69), alfa(70), alfa(71), alfa(72),
176  1 alfa(73), alfa(74), alfa(75), alfa(76), alfa(77), alfa(78),
177  2 alfa(79), alfa(80), alfa(81), alfa(82), alfa(83), alfa(84),
178  3 alfa(85), alfa(86), alfa(87), alfa(88)/
179  4 1.56501427608594704e-04, 1.36339170977445120e-04,
180  5 1.14886692029825128e-04, 9.45869093034688111e-05,
181  6 7.64498419250898258e-05, 6.07570334965197354e-05,
182  7 4.74394299290508799e-05, 3.62757512005344297e-05,
183  8 2.69939714979224901e-05, 1.93210938247939253e-05,
184  9 1.30056674793963203e-05, 7.82620866744496661e-06,
185  a 3.59257485819351583e-06, 1.44040049814251817e-07,
186  b -2.65396769697939116e-06, -4.91346867098485910e-06,
187  c -6.72739296091248287e-06, -8.17269379678657923e-06,
188  d -9.31304715093561232e-06, -1.02011418798016441e-05,
189  e -1.08805962510592880e-05, -1.13875481509603555e-05/
190  DATA alfa(89), alfa(90), alfa(91), alfa(92), alfa(93), alfa(94),
191  1 alfa(95), alfa(96), alfa(97), alfa(98), alfa(99), alfa(100),
192  2 alfa(101), alfa(102), alfa(103), alfa(104), alfa(105),
193  3 alfa(106), alfa(107), alfa(108), alfa(109), alfa(110)/
194  4 -1.17519675674556414e-05, -1.19987364870944141e-05,
195  5 3.78194199201772914e-04, 2.02471952761816167e-04,
196  6 -6.37938506318862408e-05, -2.38598230603005903e-04,
197  7 -3.10916256027361568e-04, -3.13680115247576316e-04,
198  8 -2.78950273791323387e-04, -2.28564082619141374e-04,
199  9 -1.75245280340846749e-04, -1.25544063060690348e-04,
200  a -8.22982872820208365e-05, -4.62860730588116458e-05,
201  b -1.72334302366962267e-05, 5.60690482304602267e-06,
202  c 2.31395443148286800e-05, 3.62642745856793957e-05,
203  d 4.58006124490188752e-05, 5.24595294959114050e-05,
204  e 5.68396208545815266e-05, 5.94349820393104052e-05/
205  DATA alfa(111), alfa(112), alfa(113), alfa(114), alfa(115),
206  1 alfa(116), alfa(117), alfa(118), alfa(119), alfa(120),
207  2 alfa(121), alfa(122), alfa(123), alfa(124), alfa(125),
208  3 alfa(126), alfa(127), alfa(128), alfa(129), alfa(130)/
209  4 6.06478527578421742e-05, 6.08023907788436497e-05,
210  5 6.01577894539460388e-05, 5.89199657344698500e-05,
211  6 5.72515823777593053e-05, 5.52804375585852577e-05,
212  7 5.31063773802880170e-05, 5.08069302012325706e-05,
213  8 4.84418647620094842e-05, 4.60568581607475370e-05,
214  9 -6.91141397288294174e-04, -4.29976633058871912e-04,
215  a 1.83067735980039018e-04, 6.60088147542014144e-04,
216  b 8.75964969951185931e-04, 8.77335235958235514e-04,
217  c 7.49369585378990637e-04, 5.63832329756980918e-04,
218  d 3.68059319971443156e-04, 1.88464535514455599e-04/
219  DATA alfa(131), alfa(132), alfa(133), alfa(134), alfa(135),
220  1 alfa(136), alfa(137), alfa(138), alfa(139), alfa(140),
221  2 alfa(141), alfa(142), alfa(143), alfa(144), alfa(145),
222  3 alfa(146), alfa(147), alfa(148), alfa(149), alfa(150)/
223  4 3.70663057664904149e-05, -8.28520220232137023e-05,
224  5 -1.72751952869172998e-04, -2.36314873605872983e-04,
225  6 -2.77966150694906658e-04, -3.02079514155456919e-04,
226  7 -3.12594712643820127e-04, -3.12872558758067163e-04,
227  8 -3.05678038466324377e-04, -2.93226470614557331e-04,
228  9 -2.77255655582934777e-04, -2.59103928467031709e-04,
229  a -2.39784014396480342e-04, -2.20048260045422848e-04,
230  b -2.00443911094971498e-04, -1.81358692210970687e-04,
231  c -1.63057674478657464e-04, -1.45712672175205844e-04,
232  d -1.29425421983924587e-04, -1.14245691942445952e-04/
233  DATA alfa(151), alfa(152), alfa(153), alfa(154), alfa(155),
234  1 alfa(156), alfa(157), alfa(158), alfa(159), alfa(160),
235  2 alfa(161), alfa(162), alfa(163), alfa(164), alfa(165),
236  3 alfa(166), alfa(167), alfa(168), alfa(169), alfa(170)/
237  4 1.92821964248775885e-03, 1.35592576302022234e-03,
238  5 -7.17858090421302995e-04, -2.58084802575270346e-03,
239  6 -3.49271130826168475e-03, -3.46986299340960628e-03,
240  7 -2.82285233351310182e-03, -1.88103076404891354e-03,
241  8 -8.89531718383947600e-04, 3.87912102631035228e-06,
242  9 7.28688540119691412e-04, 1.26566373053457758e-03,
243  a 1.62518158372674427e-03, 1.83203153216373172e-03,
244  b 1.91588388990527909e-03, 1.90588846755546138e-03,
245  c 1.82798982421825727e-03, 1.70389506421121530e-03,
246  d 1.55097127171097686e-03, 1.38261421852276159e-03/
247  DATA alfa(171), alfa(172), alfa(173), alfa(174), alfa(175),
248  1 alfa(176), alfa(177), alfa(178), alfa(179), alfa(180)/
249  2 1.20881424230064774e-03, 1.03676532638344962e-03,
250  3 8.71437918068619115e-04, 7.16080155297701002e-04,
251  4 5.72637002558129372e-04, 4.42089819465802277e-04,
252  5 3.24724948503090564e-04, 2.20342042730246599e-04,
253  6 1.28412898401353882e-04, 4.82005924552095464e-05/
254  DATA beta(1), beta(2), beta(3), beta(4), beta(5), beta(6),
255  1 beta(7), beta(8), beta(9), beta(10), beta(11), beta(12),
256  2 beta(13), beta(14), beta(15), beta(16), beta(17), beta(18),
257  3 beta(19), beta(20), beta(21), beta(22)/
258  4 1.79988721413553309e-02, 5.59964911064388073e-03,
259  5 2.88501402231132779e-03, 1.80096606761053941e-03,
260  6 1.24753110589199202e-03, 9.22878876572938311e-04,
261  7 7.14430421727287357e-04, 5.71787281789704872e-04,
262  8 4.69431007606481533e-04, 3.93232835462916638e-04,
263  9 3.34818889318297664e-04, 2.88952148495751517e-04,
264  a 2.52211615549573284e-04, 2.22280580798883327e-04,
265  b 1.97541838033062524e-04, 1.76836855019718004e-04,
266  c 1.59316899661821081e-04, 1.44347930197333986e-04,
267  d 1.31448068119965379e-04, 1.20245444949302884e-04,
268  e 1.10449144504599392e-04, 1.01828770740567258e-04/
269  DATA beta(23), beta(24), beta(25), beta(26), beta(27), beta(28),
270  1 beta(29), beta(30), beta(31), beta(32), beta(33), beta(34),
271  2 beta(35), beta(36), beta(37), beta(38), beta(39), beta(40),
272  3 beta(41), beta(42), beta(43), beta(44)/
273  4 9.41998224204237509e-05, 8.74130545753834437e-05,
274  5 8.13466262162801467e-05, 7.59002269646219339e-05,
275  6 7.09906300634153481e-05, 6.65482874842468183e-05,
276  7 6.25146958969275078e-05, 5.88403394426251749e-05,
277  8 -1.49282953213429172e-03, -8.78204709546389328e-04,
278  9 -5.02916549572034614e-04, -2.94822138512746025e-04,
279  a -1.75463996970782828e-04, -1.04008550460816434e-04,
280  b -5.96141953046457895e-05, -3.12038929076098340e-05,
281  c -1.26089735980230047e-05, -2.42892608575730389e-07,
282  d 8.05996165414273571e-06, 1.36507009262147391e-05,
283  e 1.73964125472926261e-05, 1.98672978842133780e-05/
284  DATA beta(45), beta(46), beta(47), beta(48), beta(49), beta(50),
285  1 beta(51), beta(52), beta(53), beta(54), beta(55), beta(56),
286  2 beta(57), beta(58), beta(59), beta(60), beta(61), beta(62),
287  3 beta(63), beta(64), beta(65), beta(66)/
288  4 2.14463263790822639e-05, 2.23954659232456514e-05,
289  5 2.28967783814712629e-05, 2.30785389811177817e-05,
290  6 2.30321976080909144e-05, 2.28236073720348722e-05,
291  7 2.25005881105292418e-05, 2.20981015361991429e-05,
292  8 2.16418427448103905e-05, 2.11507649256220843e-05,
293  9 2.06388749782170737e-05, 2.01165241997081666e-05,
294  a 1.95913450141179244e-05, 1.90689367910436740e-05,
295  b 1.85533719641636667e-05, 1.80475722259674218e-05,
296  c 5.52213076721292790e-04, 4.47932581552384646e-04,
297  d 2.79520653992020589e-04, 1.52468156198446602e-04,
298  e 6.93271105657043598e-05, 1.76258683069991397e-05/
299  DATA beta(67), beta(68), beta(69), beta(70), beta(71), beta(72),
300  1 beta(73), beta(74), beta(75), beta(76), beta(77), beta(78),
301  2 beta(79), beta(80), beta(81), beta(82), beta(83), beta(84),
302  3 beta(85), beta(86), beta(87), beta(88)/
303  4 -1.35744996343269136e-05, -3.17972413350427135e-05,
304  5 -4.18861861696693365e-05, -4.69004889379141029e-05,
305  6 -4.87665447413787352e-05, -4.87010031186735069e-05,
306  7 -4.74755620890086638e-05, -4.55813058138628452e-05,
307  8 -4.33309644511266036e-05, -4.09230193157750364e-05,
308  9 -3.84822638603221274e-05, -3.60857167535410501e-05,
309  a -3.37793306123367417e-05, -3.15888560772109621e-05,
310  b -2.95269561750807315e-05, -2.75978914828335759e-05,
311  c -2.58006174666883713e-05, -2.41308356761280200e-05,
312  d -2.25823509518346033e-05, -2.11479656768912971e-05,
313  e -1.98200638885294927e-05, -1.85909870801065077e-05/
314  DATA beta(89), beta(90), beta(91), beta(92), beta(93), beta(94),
315  1 beta(95), beta(96), beta(97), beta(98), beta(99), beta(100),
316  2 beta(101), beta(102), beta(103), beta(104), beta(105),
317  3 beta(106), beta(107), beta(108), beta(109), beta(110)/
318  4 -1.74532699844210224e-05, -1.63997823854497997e-05,
319  5 -4.74617796559959808e-04, -4.77864567147321487e-04,
320  6 -3.20390228067037603e-04, -1.61105016119962282e-04,
321  7 -4.25778101285435204e-05, 3.44571294294967503e-05,
322  8 7.97092684075674924e-05, 1.03138236708272200e-04,
323  9 1.12466775262204158e-04, 1.13103642108481389e-04,
324  a 1.08651634848774268e-04, 1.01437951597661973e-04,
325  b 9.29298396593363896e-05, 8.40293133016089978e-05,
326  c 7.52727991349134062e-05, 6.69632521975730872e-05,
327  d 5.92564547323194704e-05, 5.22169308826975567e-05,
328  e 4.58539485165360646e-05, 4.01445513891486808e-05/
329  DATA beta(111), beta(112), beta(113), beta(114), beta(115),
330  1 beta(116), beta(117), beta(118), beta(119), beta(120),
331  2 beta(121), beta(122), beta(123), beta(124), beta(125),
332  3 beta(126), beta(127), beta(128), beta(129), beta(130)/
333  4 3.50481730031328081e-05, 3.05157995034346659e-05,
334  5 2.64956119950516039e-05, 2.29363633690998152e-05,
335  6 1.97893056664021636e-05, 1.70091984636412623e-05,
336  7 1.45547428261524004e-05, 1.23886640995878413e-05,
337  8 1.04775876076583236e-05, 8.79179954978479373e-06,
338  9 7.36465810572578444e-04, 8.72790805146193976e-04,
339  a 6.22614862573135066e-04, 2.85998154194304147e-04,
340  b 3.84737672879366102e-06, -1.87906003636971558e-04,
341  c -2.97603646594554535e-04, -3.45998126832656348e-04,
342  d -3.53382470916037712e-04, -3.35715635775048757e-04/
343  DATA beta(131), beta(132), beta(133), beta(134), beta(135),
344  1 beta(136), beta(137), beta(138), beta(139), beta(140),
345  2 beta(141), beta(142), beta(143), beta(144), beta(145),
346  3 beta(146), beta(147), beta(148), beta(149), beta(150)/
347  4 -3.04321124789039809e-04, -2.66722723047612821e-04,
348  5 -2.27654214122819527e-04, -1.89922611854562356e-04,
349  6 -1.55058918599093870e-04, -1.23778240761873630e-04,
350  7 -9.62926147717644187e-05, -7.25178327714425337e-05,
351  8 -5.22070028895633801e-05, -3.50347750511900522e-05,
352  9 -2.06489761035551757e-05, -8.70106096849767054e-06,
353  a 1.13698686675100290e-06, 9.16426474122778849e-06,
354  b 1.56477785428872620e-05, 2.08223629482466847e-05,
355  c 2.48923381004595156e-05, 2.80340509574146325e-05,
356  d 3.03987774629861915e-05, 3.21156731406700616e-05/
357  DATA beta(151), beta(152), beta(153), beta(154), beta(155),
358  1 beta(156), beta(157), beta(158), beta(159), beta(160),
359  2 beta(161), beta(162), beta(163), beta(164), beta(165),
360  3 beta(166), beta(167), beta(168), beta(169), beta(170)/
361  4 -1.80182191963885708e-03, -2.43402962938042533e-03,
362  5 -1.83422663549856802e-03, -7.62204596354009765e-04,
363  6 2.39079475256927218e-04, 9.49266117176881141e-04,
364  7 1.34467449701540359e-03, 1.48457495259449178e-03,
365  8 1.44732339830617591e-03, 1.30268261285657186e-03,
366  9 1.10351597375642682e-03, 8.86047440419791759e-04,
367  a 6.73073208165665473e-04, 4.77603872856582378e-04,
368  b 3.05991926358789362e-04, 1.60315694594721630e-04,
369  c 4.00749555270613286e-05, -5.66607461635251611e-05,
370  d -1.32506186772982638e-04, -1.90296187989614057e-04/
371  DATA beta(171), beta(172), beta(173), beta(174), beta(175),
372  1 beta(176), beta(177), beta(178), beta(179), beta(180),
373  2 beta(181), beta(182), beta(183), beta(184), beta(185),
374  3 beta(186), beta(187), beta(188), beta(189), beta(190)/
375  4 -2.32811450376937408e-04, -2.62628811464668841e-04,
376  5 -2.82050469867598672e-04, -2.93081563192861167e-04,
377  6 -2.97435962176316616e-04, -2.96557334239348078e-04,
378  7 -2.91647363312090861e-04, -2.83696203837734166e-04,
379  8 -2.73512317095673346e-04, -2.61750155806768580e-04,
380  9 6.38585891212050914e-03, 9.62374215806377941e-03,
381  a 7.61878061207001043e-03, 2.83219055545628054e-03,
382  b -2.09841352012720090e-03, -5.73826764216626498e-03,
383  c -7.70804244495414620e-03, -8.21011692264844401e-03,
384  d -7.65824520346905413e-03, -6.47209729391045177e-03/
385  DATA beta(191), beta(192), beta(193), beta(194), beta(195),
386  1 beta(196), beta(197), beta(198), beta(199), beta(200),
387  2 beta(201), beta(202), beta(203), beta(204), beta(205),
388  3 beta(206), beta(207), beta(208), beta(209), beta(210)/
389  4 -4.99132412004966473e-03, -3.45612289713133280e-03,
390  5 -2.01785580014170775e-03, -7.59430686781961401e-04,
391  6 2.84173631523859138e-04, 1.10891667586337403e-03,
392  7 1.72901493872728771e-03, 2.16812590802684701e-03,
393  8 2.45357710494539735e-03, 2.61281821058334862e-03,
394  9 2.67141039656276912e-03, 2.65203073395980430e-03,
395  a 2.57411652877287315e-03, 2.45389126236094427e-03,
396  b 2.30460058071795494e-03, 2.13684837686712662e-03,
397  c 1.95896528478870911e-03, 1.77737008679454412e-03,
398  d 1.59690280765839059e-03, 1.42111975664438546e-03/
399  DATA gama(1), gama(2), gama(3), gama(4), gama(5), gama(6),
400  1 gama(7), gama(8), gama(9), gama(10), gama(11), gama(12),
401  2 gama(13), gama(14), gama(15), gama(16), gama(17), gama(18),
402  3 gama(19), gama(20), gama(21), gama(22)/
403  4 6.29960524947436582e-01, 2.51984209978974633e-01,
404  5 1.54790300415655846e-01, 1.10713062416159013e-01,
405  6 8.57309395527394825e-02, 6.97161316958684292e-02,
406  7 5.86085671893713576e-02, 5.04698873536310685e-02,
407  8 4.42600580689154809e-02, 3.93720661543509966e-02,
408  9 3.54283195924455368e-02, 3.21818857502098231e-02,
409  a 2.94646240791157679e-02, 2.71581677112934479e-02,
410  b 2.51768272973861779e-02, 2.34570755306078891e-02,
411  c 2.19508390134907203e-02, 2.06210828235646240e-02,
412  d 1.94388240897880846e-02, 1.83810633800683158e-02,
413  e 1.74293213231963172e-02, 1.65685837786612353e-02/
414  DATA gama(23), gama(24), gama(25), gama(26), gama(27), gama(28),
415  1 gama(29), gama(30)/
416  2 1.57865285987918445e-02, 1.50729501494095594e-02,
417  3 1.44193250839954639e-02, 1.38184805735341786e-02,
418  4 1.32643378994276568e-02, 1.27517121970498651e-02,
419  5 1.22761545318762767e-02, 1.18338262398482403e-02/
420  DATA ex1, ex2, hpi, pi, thpi /
421  1 3.33333333333333333e-01, 6.66666666666666667e-01,
422  2 1.57079632679489662e+00, 3.14159265358979324e+00,
423  3 4.71238898038468986e+00/
424  DATA czero, cone / (0.0e0,0.0e0), (1.0e0,0.0e0) /
425 C
426  rfnu = 1.0e0/fnu
427 C ZB = Z*CMPLX(RFNU,0.0E0)
428 C-----------------------------------------------------------------------
429 C OVERFLOW TEST (Z/FNU TOO SMALL)
430 C-----------------------------------------------------------------------
431  tstr = REAL(z)
432  tsti = aimag(z)
433  test = r1mach(1)*1.0e+3
434  ac = fnu*test
435  IF (abs(tstr).GT.ac .OR. abs(tsti).GT.ac) go to 15
436  ac = 2.0e0*abs(alog(test))+fnu
437  zeta1 = cmplx(ac,0.0e0)
438  zeta2 = cmplx(fnu,0.0e0)
439  phi=cone
440  arg=cone
441  RETURN
442  15 CONTINUE
443  zb = z*cmplx(rfnu,0.0e0)
444  rfnu2 = rfnu*rfnu
445 C-----------------------------------------------------------------------
446 C COMPUTE IN THE FOURTH QUADRANT
447 C-----------------------------------------------------------------------
448  fn13 = fnu**ex1
449  fn23 = fn13*fn13
450  rfn13 = cmplx(1.0e0/fn13,0.0e0)
451  w2 = cone - zb*zb
452  aw2 = cabs(w2)
453  IF (aw2.GT.0.25e0) go to 130
454 C-----------------------------------------------------------------------
455 C POWER SERIES FOR CABS(W2).LE.0.25E0
456 C-----------------------------------------------------------------------
457  k = 1
458  p(1) = cone
459  suma = cmplx(gama(1),0.0e0)
460  ap(1) = 1.0e0
461  IF (aw2.LT.tol) go to 20
462  DO 10 k=2,30
463  p(k) = p(k-1)*w2
464  suma = suma + p(k)*cmplx(gama(k),0.0e0)
465  ap(k) = ap(k-1)*aw2
466  IF (ap(k).LT.tol) go to 20
467  10 CONTINUE
468  k = 30
469  20 CONTINUE
470  kmax = k
471  zeta = w2*suma
472  arg = zeta*cmplx(fn23,0.0e0)
473  za = csqrt(suma)
474  zeta2 = csqrt(w2)*cmplx(fnu,0.0e0)
475  zeta1 = zeta2*(cone+zeta*za*cmplx(ex2,0.0e0))
476  za = za + za
477  phi = csqrt(za)*rfn13
478  IF (ipmtr.EQ.1) go to 120
479 C-----------------------------------------------------------------------
480 C SUM SERIES FOR ASUM AND BSUM
481 C-----------------------------------------------------------------------
482  sumb = czero
483  DO 30 k=1,kmax
484  sumb = sumb + p(k)*cmplx(beta(k),0.0e0)
485  30 CONTINUE
486  asum = czero
487  bsum = sumb
488  l1 = 0
489  l2 = 30
490  btol = tol*cabs(bsum)
491  atol = tol
492  pp = 1.0e0
493  ias = 0
494  ibs = 0
495  IF (rfnu2.LT.tol) go to 110
496  DO 100 is=2,7
497  atol = atol/rfnu2
498  pp = pp*rfnu2
499  IF (ias.EQ.1) go to 60
500  suma = czero
501  DO 40 k=1,kmax
502  m = l1 + k
503  suma = suma + p(k)*cmplx(alfa(m),0.0e0)
504  IF (ap(k).LT.atol) go to 50
505  40 CONTINUE
506  50 CONTINUE
507  asum = asum + suma*cmplx(pp,0.0e0)
508  IF (pp.LT.tol) ias = 1
509  60 CONTINUE
510  IF (ibs.EQ.1) go to 90
511  sumb = czero
512  DO 70 k=1,kmax
513  m = l2 + k
514  sumb = sumb + p(k)*cmplx(beta(m),0.0e0)
515  IF (ap(k).LT.atol) go to 80
516  70 CONTINUE
517  80 CONTINUE
518  bsum = bsum + sumb*cmplx(pp,0.0e0)
519  IF (pp.LT.btol) ibs = 1
520  90 CONTINUE
521  IF (ias.EQ.1 .AND. ibs.EQ.1) go to 110
522  l1 = l1 + 30
523  l2 = l2 + 30
524  100 CONTINUE
525  110 CONTINUE
526  asum = asum + cone
527  pp = rfnu*REAL(rfn13)
528  bsum = bsum*cmplx(pp,0.0e0)
529  120 CONTINUE
530  RETURN
531 C-----------------------------------------------------------------------
532 C CABS(W2).GT.0.25E0
533 C-----------------------------------------------------------------------
534  130 CONTINUE
535  w = csqrt(w2)
536  wr = REAL(w)
537  wi = aimag(w)
538  IF (wr.LT.0.0e0) wr = 0.0e0
539  IF (wi.LT.0.0e0) wi = 0.0e0
540  w = cmplx(wr,wi)
541  za = (cone+w)/zb
542  zc = clog(za)
543  zcr = REAL(zc)
544  zci = aimag(zc)
545  IF (zci.LT.0.0e0) zci = 0.0e0
546  IF (zci.GT.hpi) zci = hpi
547  IF (zcr.LT.0.0e0) zcr = 0.0e0
548  zc = cmplx(zcr,zci)
549  zth = (zc-w)*cmplx(1.5e0,0.0e0)
550  cfnu = cmplx(fnu,0.0e0)
551  zeta1 = zc*cfnu
552  zeta2 = w*cfnu
553  azth = cabs(zth)
554  zthr = REAL(zth)
555  zthi = aimag(zth)
556  ang = thpi
557  IF (zthr.GE.0.0e0 .AND. zthi.LT.0.0e0) go to 140
558  ang = hpi
559  IF (zthr.EQ.0.0e0) go to 140
560  ang = atan(zthi/zthr)
561  IF (zthr.LT.0.0e0) ang = ang + pi
562  140 CONTINUE
563  pp = azth**ex2
564  ang = ang*ex2
565  zetar = pp*cos(ang)
566  zetai = pp*sin(ang)
567  IF (zetai.LT.0.0e0) zetai = 0.0e0
568  zeta = cmplx(zetar,zetai)
569  arg = zeta*cmplx(fn23,0.0e0)
570  rtzta = zth/zeta
571  za = rtzta/w
572  phi = csqrt(za+za)*rfn13
573  IF (ipmtr.EQ.1) go to 120
574  tfn = cmplx(rfnu,0.0e0)/w
575  rzth = cmplx(rfnu,0.0e0)/zth
576  zc = rzth*cmplx(ar(2),0.0e0)
577  t2 = cone/w2
578  up(2) = (t2*cmplx(c(2),0.0e0)+cmplx(c(3),0.0e0))*tfn
579  bsum = up(2) + zc
580  asum = czero
581  IF (rfnu.LT.tol) go to 220
582  przth = rzth
583  ptfn = tfn
584  up(1) = cone
585  pp = 1.0e0
586  bsumr = REAL(bsum)
587  bsumi = aimag(bsum)
588  btol = tol*(abs(bsumr)+abs(bsumi))
589  ks = 0
590  kp1 = 2
591  l = 3
592  ias = 0
593  ibs = 0
594  DO 210 lr=2,12,2
595  lrp1 = lr + 1
596 C-----------------------------------------------------------------------
597 C COMPUTE TWO ADDITIONAL CR, DR, AND UP FOR TWO MORE TERMS IN
598 C NEXT SUMA AND SUMB
599 C-----------------------------------------------------------------------
600  DO 160 k=lr,lrp1
601  ks = ks + 1
602  kp1 = kp1 + 1
603  l = l + 1
604  za = cmplx(c(l),0.0e0)
605  DO 150 j=2,kp1
606  l = l + 1
607  za = za*t2 + cmplx(c(l),0.0e0)
608  150 CONTINUE
609  ptfn = ptfn*tfn
610  up(kp1) = ptfn*za
611  cr(ks) = przth*cmplx(br(ks+1),0.0e0)
612  przth = przth*rzth
613  dr(ks) = przth*cmplx(ar(ks+2),0.0e0)
614  160 CONTINUE
615  pp = pp*rfnu2
616  IF (ias.EQ.1) go to 180
617  suma = up(lrp1)
618  ju = lrp1
619  DO 170 jr=1,lr
620  ju = ju - 1
621  suma = suma + cr(jr)*up(ju)
622  170 CONTINUE
623  asum = asum + suma
624  asumr = REAL(asum)
625  asumi = aimag(asum)
626  test = abs(asumr) + abs(asumi)
627  IF (pp.LT.tol .AND. test.LT.tol) ias = 1
628  180 CONTINUE
629  IF (ibs.EQ.1) go to 200
630  sumb = up(lr+2) + up(lrp1)*zc
631  ju = lrp1
632  DO 190 jr=1,lr
633  ju = ju - 1
634  sumb = sumb + dr(jr)*up(ju)
635  190 CONTINUE
636  bsum = bsum + sumb
637  bsumr = REAL(bsum)
638  bsumi = aimag(bsum)
639  test = abs(bsumr) + abs(bsumi)
640  IF (pp.LT.btol .AND. test.LT.tol) ibs = 1
641  200 CONTINUE
642  IF (ias.EQ.1 .AND. ibs.EQ.1) go to 220
643  210 CONTINUE
644  220 CONTINUE
645  asum = asum + cone
646  bsum = -bsum*rfn13/rtzta
647  go to 120
648  END
Generated on Mon Dec 30 2013 03:04:45 for GNU Octave by   doxygen 1.8.1.2