|
3217
|
1 SUBROUTINE ZUNIK(ZRR, ZRI, FNU, IKFLG, IPMTR, TOL, INIT, PHIR, |
|
|
2 * PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI) |
|
|
3 C***BEGIN PROLOGUE ZUNIK |
|
|
4 C***REFER TO ZBESI,ZBESK |
|
|
5 C |
|
|
6 C ZUNIK COMPUTES PARAMETERS FOR THE UNIFORM ASYMPTOTIC |
|
|
7 C EXPANSIONS OF THE I AND K FUNCTIONS ON IKFLG= 1 OR 2 |
|
|
8 C RESPECTIVELY BY |
|
|
9 C |
|
|
10 C W(FNU,ZR) = PHI*EXP(ZETA)*SUM |
|
|
11 C |
|
|
12 C WHERE ZETA=-ZETA1 + ZETA2 OR |
|
|
13 C ZETA1 - ZETA2 |
|
|
14 C |
|
|
15 C THE FIRST CALL MUST HAVE INIT=0. SUBSEQUENT CALLS WITH THE |
|
|
16 C SAME ZR AND FNU WILL RETURN THE I OR K FUNCTION ON IKFLG= |
|
|
17 C 1 OR 2 WITH NO CHANGE IN INIT. CWRK IS A COMPLEX WORK |
|
|
18 C ARRAY. IPMTR=0 COMPUTES ALL PARAMETERS. IPMTR=1 COMPUTES PHI, |
|
|
19 C ZETA1,ZETA2. |
|
|
20 C |
|
|
21 C***ROUTINES CALLED ZDIV,XZLOG,XZSQRT,D1MACH |
|
|
22 C***END PROLOGUE ZUNIK |
|
|
23 C COMPLEX CFN,CON,CONE,CRFN,CWRK,CZERO,PHI,S,SR,SUM,T,T2,ZETA1, |
|
|
24 C *ZETA2,ZN,ZR |
|
|
25 DOUBLE PRECISION AC, C, CON, CONEI, CONER, CRFNI, CRFNR, CWRKI, |
|
|
26 * CWRKR, FNU, PHII, PHIR, RFN, SI, SR, SRI, SRR, STI, STR, SUMI, |
|
|
27 * SUMR, TEST, TI, TOL, TR, T2I, T2R, ZEROI, ZEROR, ZETA1I, ZETA1R, |
|
|
28 * ZETA2I, ZETA2R, ZNI, ZNR, ZRI, ZRR, D1MACH |
|
|
29 INTEGER I, IDUM, IKFLG, INIT, IPMTR, J, K, L |
|
|
30 DIMENSION C(120), CWRKR(16), CWRKI(16), CON(2) |
|
|
31 DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / |
|
|
32 DATA CON(1), CON(2) / |
|
|
33 1 3.98942280401432678D-01, 1.25331413731550025D+00 / |
|
|
34 DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10), |
|
|
35 1 C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18), |
|
|
36 2 C(19), C(20), C(21), C(22), C(23), C(24)/ |
|
|
37 3 1.00000000000000000D+00, -2.08333333333333333D-01, |
|
|
38 4 1.25000000000000000D-01, 3.34201388888888889D-01, |
|
|
39 5 -4.01041666666666667D-01, 7.03125000000000000D-02, |
|
|
40 6 -1.02581259645061728D+00, 1.84646267361111111D+00, |
|
|
41 7 -8.91210937500000000D-01, 7.32421875000000000D-02, |
|
|
42 8 4.66958442342624743D+00, -1.12070026162229938D+01, |
|
|
43 9 8.78912353515625000D+00, -2.36408691406250000D+00, |
|
|
44 A 1.12152099609375000D-01, -2.82120725582002449D+01, |
|
|
45 B 8.46362176746007346D+01, -9.18182415432400174D+01, |
|
|
46 C 4.25349987453884549D+01, -7.36879435947963170D+00, |
|
|
47 D 2.27108001708984375D-01, 2.12570130039217123D+02, |
|
|
48 E -7.65252468141181642D+02, 1.05999045252799988D+03/ |
|
|
49 DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32), |
|
|
50 1 C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40), |
|
|
51 2 C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/ |
|
|
52 3 -6.99579627376132541D+02, 2.18190511744211590D+02, |
|
|
53 4 -2.64914304869515555D+01, 5.72501420974731445D-01, |
|
|
54 5 -1.91945766231840700D+03, 8.06172218173730938D+03, |
|
|
55 6 -1.35865500064341374D+04, 1.16553933368645332D+04, |
|
|
56 7 -5.30564697861340311D+03, 1.20090291321635246D+03, |
|
|
57 8 -1.08090919788394656D+02, 1.72772750258445740D+00, |
|
|
58 9 2.02042913309661486D+04, -9.69805983886375135D+04, |
|
|
59 A 1.92547001232531532D+05, -2.03400177280415534D+05, |
|
|
60 B 1.22200464983017460D+05, -4.11926549688975513D+04, |
|
|
61 C 7.10951430248936372D+03, -4.93915304773088012D+02, |
|
|
62 D 6.07404200127348304D+00, -2.42919187900551333D+05, |
|
|
63 E 1.31176361466297720D+06, -2.99801591853810675D+06/ |
|
|
64 DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56), |
|
|
65 1 C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64), |
|
|
66 2 C(65), C(66), C(67), C(68), C(69), C(70), C(71), C(72)/ |
|
|
67 3 3.76327129765640400D+06, -2.81356322658653411D+06, |
|
|
68 4 1.26836527332162478D+06, -3.31645172484563578D+05, |
|
|
69 5 4.52187689813627263D+04, -2.49983048181120962D+03, |
|
|
70 6 2.43805296995560639D+01, 3.28446985307203782D+06, |
|
|
71 7 -1.97068191184322269D+07, 5.09526024926646422D+07, |
|
|
72 8 -7.41051482115326577D+07, 6.63445122747290267D+07, |
|
|
73 9 -3.75671766607633513D+07, 1.32887671664218183D+07, |
|
|
74 A -2.78561812808645469D+06, 3.08186404612662398D+05, |
|
|
75 B -1.38860897537170405D+04, 1.10017140269246738D+02, |
|
|
76 C -4.93292536645099620D+07, 3.25573074185765749D+08, |
|
|
77 D -9.39462359681578403D+08, 1.55359689957058006D+09, |
|
|
78 E -1.62108055210833708D+09, 1.10684281682301447D+09/ |
|
|
79 DATA C(73), C(74), C(75), C(76), C(77), C(78), C(79), C(80), |
|
|
80 1 C(81), C(82), C(83), C(84), C(85), C(86), C(87), C(88), |
|
|
81 2 C(89), C(90), C(91), C(92), C(93), C(94), C(95), C(96)/ |
|
|
82 3 -4.95889784275030309D+08, 1.42062907797533095D+08, |
|
|
83 4 -2.44740627257387285D+07, 2.24376817792244943D+06, |
|
|
84 5 -8.40054336030240853D+04, 5.51335896122020586D+02, |
|
|
85 6 8.14789096118312115D+08, -5.86648149205184723D+09, |
|
|
86 7 1.86882075092958249D+10, -3.46320433881587779D+10, |
|
|
87 8 4.12801855797539740D+10, -3.30265997498007231D+10, |
|
|
88 9 1.79542137311556001D+10, -6.56329379261928433D+09, |
|
|
89 A 1.55927986487925751D+09, -2.25105661889415278D+08, |
|
|
90 B 1.73951075539781645D+07, -5.49842327572288687D+05, |
|
|
91 C 3.03809051092238427D+03, -1.46792612476956167D+10, |
|
|
92 D 1.14498237732025810D+11, -3.99096175224466498D+11, |
|
|
93 E 8.19218669548577329D+11, -1.09837515608122331D+12/ |
|
|
94 DATA C(97), C(98), C(99), C(100), C(101), C(102), C(103), C(104), |
|
|
95 1 C(105), C(106), C(107), C(108), C(109), C(110), C(111), |
|
|
96 2 C(112), C(113), C(114), C(115), C(116), C(117), C(118)/ |
|
|
97 3 1.00815810686538209D+12, -6.45364869245376503D+11, |
|
|
98 4 2.87900649906150589D+11, -8.78670721780232657D+10, |
|
|
99 5 1.76347306068349694D+10, -2.16716498322379509D+09, |
|
|
100 6 1.43157876718888981D+08, -3.87183344257261262D+06, |
|
|
101 7 1.82577554742931747D+04, 2.86464035717679043D+11, |
|
|
102 8 -2.40629790002850396D+12, 9.10934118523989896D+12, |
|
|
103 9 -2.05168994109344374D+13, 3.05651255199353206D+13, |
|
|
104 A -3.16670885847851584D+13, 2.33483640445818409D+13, |
|
|
105 B -1.23204913055982872D+13, 4.61272578084913197D+12, |
|
|
106 C -1.19655288019618160D+12, 2.05914503232410016D+11, |
|
|
107 D -2.18229277575292237D+10, 1.24700929351271032D+09/ |
|
|
108 DATA C(119), C(120)/ |
|
|
109 1 -2.91883881222208134D+07, 1.18838426256783253D+05/ |
|
|
110 C |
|
|
111 IF (INIT.NE.0) GO TO 40 |
|
|
112 C----------------------------------------------------------------------- |
|
|
113 C INITIALIZE ALL VARIABLES |
|
|
114 C----------------------------------------------------------------------- |
|
|
115 RFN = 1.0D0/FNU |
|
|
116 C----------------------------------------------------------------------- |
|
|
117 C OVERFLOW TEST (ZR/FNU TOO SMALL) |
|
|
118 C----------------------------------------------------------------------- |
|
|
119 TEST = D1MACH(1)*1.0D+3 |
|
|
120 AC = FNU*TEST |
|
|
121 IF (DABS(ZRR).GT.AC .OR. DABS(ZRI).GT.AC) GO TO 15 |
|
|
122 ZETA1R = 2.0D0*DABS(DLOG(TEST))+FNU |
|
|
123 ZETA1I = 0.0D0 |
|
|
124 ZETA2R = FNU |
|
|
125 ZETA2I = 0.0D0 |
|
|
126 PHIR = 1.0D0 |
|
|
127 PHII = 0.0D0 |
|
|
128 RETURN |
|
|
129 15 CONTINUE |
|
|
130 TR = ZRR*RFN |
|
|
131 TI = ZRI*RFN |
|
|
132 SR = CONER + (TR*TR-TI*TI) |
|
|
133 SI = CONEI + (TR*TI+TI*TR) |
|
|
134 CALL XZSQRT(SR, SI, SRR, SRI) |
|
|
135 STR = CONER + SRR |
|
|
136 STI = CONEI + SRI |
|
|
137 CALL ZDIV(STR, STI, TR, TI, ZNR, ZNI) |
|
|
138 CALL XZLOG(ZNR, ZNI, STR, STI, IDUM) |
|
|
139 ZETA1R = FNU*STR |
|
|
140 ZETA1I = FNU*STI |
|
|
141 ZETA2R = FNU*SRR |
|
|
142 ZETA2I = FNU*SRI |
|
|
143 CALL ZDIV(CONER, CONEI, SRR, SRI, TR, TI) |
|
|
144 SRR = TR*RFN |
|
|
145 SRI = TI*RFN |
|
|
146 CALL XZSQRT(SRR, SRI, CWRKR(16), CWRKI(16)) |
|
|
147 PHIR = CWRKR(16)*CON(IKFLG) |
|
|
148 PHII = CWRKI(16)*CON(IKFLG) |
|
|
149 IF (IPMTR.NE.0) RETURN |
|
|
150 CALL ZDIV(CONER, CONEI, SR, SI, T2R, T2I) |
|
|
151 CWRKR(1) = CONER |
|
|
152 CWRKI(1) = CONEI |
|
|
153 CRFNR = CONER |
|
|
154 CRFNI = CONEI |
|
|
155 AC = 1.0D0 |
|
|
156 L = 1 |
|
|
157 DO 20 K=2,15 |
|
|
158 SR = ZEROR |
|
|
159 SI = ZEROI |
|
|
160 DO 10 J=1,K |
|
|
161 L = L + 1 |
|
|
162 STR = SR*T2R - SI*T2I + C(L) |
|
|
163 SI = SR*T2I + SI*T2R |
|
|
164 SR = STR |
|
|
165 10 CONTINUE |
|
|
166 STR = CRFNR*SRR - CRFNI*SRI |
|
|
167 CRFNI = CRFNR*SRI + CRFNI*SRR |
|
|
168 CRFNR = STR |
|
|
169 CWRKR(K) = CRFNR*SR - CRFNI*SI |
|
|
170 CWRKI(K) = CRFNR*SI + CRFNI*SR |
|
|
171 AC = AC*RFN |
|
|
172 TEST = DABS(CWRKR(K)) + DABS(CWRKI(K)) |
|
|
173 IF (AC.LT.TOL .AND. TEST.LT.TOL) GO TO 30 |
|
|
174 20 CONTINUE |
|
|
175 K = 15 |
|
|
176 30 CONTINUE |
|
|
177 INIT = K |
|
|
178 40 CONTINUE |
|
|
179 IF (IKFLG.EQ.2) GO TO 60 |
|
|
180 C----------------------------------------------------------------------- |
|
|
181 C COMPUTE SUM FOR THE I FUNCTION |
|
|
182 C----------------------------------------------------------------------- |
|
|
183 SR = ZEROR |
|
|
184 SI = ZEROI |
|
|
185 DO 50 I=1,INIT |
|
|
186 SR = SR + CWRKR(I) |
|
|
187 SI = SI + CWRKI(I) |
|
|
188 50 CONTINUE |
|
|
189 SUMR = SR |
|
|
190 SUMI = SI |
|
|
191 PHIR = CWRKR(16)*CON(1) |
|
|
192 PHII = CWRKI(16)*CON(1) |
|
|
193 RETURN |
|
|
194 60 CONTINUE |
|
|
195 C----------------------------------------------------------------------- |
|
|
196 C COMPUTE SUM FOR THE K FUNCTION |
|
|
197 C----------------------------------------------------------------------- |
|
|
198 SR = ZEROR |
|
|
199 SI = ZEROI |
|
|
200 TR = CONER |
|
|
201 DO 70 I=1,INIT |
|
|
202 SR = SR + TR*CWRKR(I) |
|
|
203 SI = SI + TR*CWRKI(I) |
|
|
204 TR = -TR |
|
|
205 70 CONTINUE |
|
|
206 SUMR = SR |
|
|
207 SUMI = SI |
|
|
208 PHIR = CWRKR(16)*CON(2) |
|
|
209 PHII = CWRKI(16)*CON(2) |
|
|
210 RETURN |
|
|
211 END |
