GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
cmplx
double complex cmplx
Definition: Faddeeva.cc:217
cunhj
subroutine cunhj(Z, FNU, IPMTR, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM)
Definition: cunhj.f:3
real
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137
abs
static T abs(T x)
Definition: pr-output.cc:1678
r1mach
real function r1mach(i)
Definition: r1mach.f:23