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
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
Generated on Mon Dec 30 2013 03:04:45 for GNU Octave by   doxygen 1.8.1.2