GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
zunik.f
Go to the documentation of this file.
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
d1mach
double precision function d1mach(i)
Definition: d1mach.f:23
xzlog
subroutine xzlog(AR, AI, BR, BI, IERR)
Definition: xzlog.f:2
xzsqrt
subroutine xzsqrt(AR, AI, BR, BI)
Definition: xzsqrt.f:2
zdiv
subroutine zdiv(AR, AI, BR, BI, CR, CI)
Definition: zdiv.f:2
zunik
subroutine zunik(ZRR, ZRI, FNU, IKFLG, IPMTR, TOL, INIT, PHIR, PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
Definition: zunik.f:3