GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
cunk1.f
Go to the documentation of this file.
1  SUBROUTINE cunk1(Z, FNU, KODE, MR, N, Y, NZ, TOL, ELIM, ALIM)
2 C***BEGIN PROLOGUE CUNK1
3 C***REFER TO CBESK
4 C
5 C CUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
6 C RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
7 C UNIFORM ASYMPTOTIC EXPANSION.
8 C MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
9 C NZ=-1 MEANS AN OVERFLOW WILL OCCUR
10 C
11 C***ROUTINES CALLED CS1S2,CUCHK,CUNIK,R1MACH
12 C***END PROLOGUE CUNK1
13  COMPLEX CFN, CK, CONE, CRSC, CS, CSCL, CSGN, CSPN, CSR, CSS,
14  * CWRK, CY, CZERO, C1, C2, PHI, RZ, SUM, S1, S2, Y, Z,
15  * ZETA1, ZETA2, ZR, PHID, ZETA1D, ZETA2D, SUMD
16  REAL ALIM, ANG, APHI, ASC, ASCLE, BRY, CPN, C2I, C2M, C2R, ELIM,
17  * FMR, FN, FNF, FNU, PI, RS1, SGN, SPN, TOL, X, R1MACH
18  INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG,
19  * KK, KODE, MR, N, NW, NZ, J, IPARD, INITD, IC
20  dimension bry(3), init(2), y(n), sum(2), phi(2), zeta1(2),
21  * zeta2(2), cy(2), cwrk(16,3), css(3), csr(3)
22  DATA czero, cone / (0.0e0,0.0e0) , (1.0e0,0.0e0) /
23  DATA pi / 3.14159265358979324e0 /
24 C
25  kdflg = 1
26  nz = 0
27 C-----------------------------------------------------------------------
28 C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
29 C THE UNDERFLOW LIMIT
30 C-----------------------------------------------------------------------
31  cscl = cmplx(1.0e0/tol,0.0e0)
32  crsc = cmplx(tol,0.0e0)
33  css(1) = cscl
34  css(2) = cone
35  css(3) = crsc
36  csr(1) = crsc
37  csr(2) = cone
38  csr(3) = cscl
39  bry(1) = 1.0e+3*r1mach(1)/tol
40  bry(2) = 1.0e0/bry(1)
41  bry(3) = r1mach(2)
42  x = REAL(Z)
43  zr = z
44  IF (x.LT.0.0e0) zr = -z
45  j=2
46  DO 70 i=1,n
47 C-----------------------------------------------------------------------
48 C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
49 C-----------------------------------------------------------------------
50  j = 3 - j
51  fn = fnu + float(i-1)
52  init(j) = 0
53  CALL cunik(zr, fn, 2, 0, tol, init(j), phi(j), zeta1(j),
54  * zeta2(j), sum(j), cwrk(1,j))
55  IF (kode.EQ.1) GO TO 20
56  cfn = cmplx(fn,0.0e0)
57  s1 = zeta1(j) - cfn*(cfn/(zr+zeta2(j)))
58  GO TO 30
59  20 CONTINUE
60  s1 = zeta1(j) - zeta2(j)
61  30 CONTINUE
62 C-----------------------------------------------------------------------
63 C TEST FOR UNDERFLOW AND OVERFLOW
64 C-----------------------------------------------------------------------
65  rs1 = REAL(S1)
66  IF (abs(rs1).GT.elim) GO TO 60
67  IF (kdflg.EQ.1) kflag = 2
68  IF (abs(rs1).LT.alim) GO TO 40
69 C-----------------------------------------------------------------------
70 C REFINE TEST AND SCALE
71 C-----------------------------------------------------------------------
72  aphi = cabs(phi(j))
73  rs1 = rs1 + alog(aphi)
74  IF (abs(rs1).GT.elim) GO TO 60
75  IF (kdflg.EQ.1) kflag = 1
76  IF (rs1.LT.0.0e0) GO TO 40
77  IF (kdflg.EQ.1) kflag = 3
78  40 CONTINUE
79 C-----------------------------------------------------------------------
80 C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
81 C EXPONENT EXTREMES
82 C-----------------------------------------------------------------------
83  s2 = phi(j)*sum(j)
84  c2r = REAL(S1)
85  c2i = aimag(s1)
86  c2m = exp(c2r)*REAL(CSS(KFLAG))
87  s1 = cmplx(c2m,0.0e0)*cmplx(cos(c2i),sin(c2i))
88  s2 = s2*s1
89  IF (kflag.NE.1) GO TO 50
90  CALL cuchk(s2, nw, bry(1), tol)
91  IF (nw.NE.0) GO TO 60
92  50 CONTINUE
93  cy(kdflg) = s2
94  y(i) = s2*csr(kflag)
95  IF (kdflg.EQ.2) GO TO 75
96  kdflg = 2
97  GO TO 70
98  60 CONTINUE
99  IF (rs1.GT.0.0e0) GO TO 290
100 C-----------------------------------------------------------------------
101 C FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
102 C-----------------------------------------------------------------------
103  IF (x.LT.0.0e0) GO TO 290
104  kdflg = 1
105  y(i) = czero
106  nz=nz+1
107  IF (i.EQ.1) GO TO 70
108  IF (y(i-1).EQ.czero) GO TO 70
109  y(i-1) = czero
110  nz=nz+1
111  70 CONTINUE
112  i=n
113  75 CONTINUE
114  rz = cmplx(2.0e0,0.0e0)/zr
115  ck = cmplx(fn,0.0e0)*rz
116  ib = i+1
117  IF (n.LT.ib) GO TO 160
118 C-----------------------------------------------------------------------
119 C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW, SET SEQUENCE TO ZERO
120 C ON UNDERFLOW
121 C-----------------------------------------------------------------------
122  fn = fnu+float(n-1)
123  ipard = 1
124  IF (mr.NE.0) ipard = 0
125  initd = 0
126  CALL cunik(zr,fn,2,ipard,tol,initd,phid,zeta1d,zeta2d,sumd,
127  *cwrk(1,3))
128  IF (kode.EQ.1) GO TO 80
129  cfn=cmplx(fn,0.0e0)
130  s1=zeta1d-cfn*(cfn/(zr+zeta2d))
131  GO TO 90
132  80 CONTINUE
133  s1=zeta1d-zeta2d
134  90 CONTINUE
135  rs1=REAL(S1)
136  IF (abs(rs1).GT.elim) GO TO 95
137  IF (abs(rs1).LT.alim) GO TO 100
138 C-----------------------------------------------------------------------
139 C REFINE ESTIMATE AND TEST
140 C-----------------------------------------------------------------------
141  aphi=cabs(phid)
142  rs1=rs1+alog(aphi)
143  IF (abs(rs1).LT.elim) GO TO 100
144  95 CONTINUE
145  IF (rs1.GT.0.0e0) GO TO 290
146 C-----------------------------------------------------------------------
147 C FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
148 C-----------------------------------------------------------------------
149  IF (x.LT.0.0e0) GO TO 290
150  nz=n
151  DO 96 i=1,n
152  y(i) = czero
153  96 CONTINUE
154  RETURN
155  100 CONTINUE
156 C-----------------------------------------------------------------------
157 C RECUR FORWARD FOR REMAINDER OF THE SEQUENCE
158 C-----------------------------------------------------------------------
159  s1 = cy(1)
160  s2 = cy(2)
161  c1 = csr(kflag)
162  ascle = bry(kflag)
163  DO 120 i=ib,n
164  c2 = s2
165  s2 = ck*s2 + s1
166  s1 = c2
167  ck = ck + rz
168  c2 = s2*c1
169  y(i) = c2
170  IF (kflag.GE.3) GO TO 120
171  c2r = REAL(C2)
172  c2i = aimag(c2)
173  c2r = abs(c2r)
174  c2i = abs(c2i)
175  c2m = amax1(c2r,c2i)
176  IF (c2m.LE.ascle) GO TO 120
177  kflag = kflag + 1
178  ascle = bry(kflag)
179  s1 = s1*c1
180  s2 = c2
181  s1 = s1*css(kflag)
182  s2 = s2*css(kflag)
183  c1 = csr(kflag)
184  120 CONTINUE
185  160 CONTINUE
186  IF (mr.EQ.0) RETURN
187 C-----------------------------------------------------------------------
188 C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0E0
189 C-----------------------------------------------------------------------
190  nz = 0
191  fmr = float(mr)
192  sgn = -sign(pi,fmr)
193 C-----------------------------------------------------------------------
194 C CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
195 C-----------------------------------------------------------------------
196  csgn = cmplx(0.0e0,sgn)
197  inu = int(fnu)
198  fnf = fnu - float(inu)
199  ifn = inu + n - 1
200  ang = fnf*sgn
201  cpn = cos(ang)
202  spn = sin(ang)
203  cspn = cmplx(cpn,spn)
204  IF (mod(ifn,2).EQ.1) cspn = -cspn
205  asc = bry(1)
206  kk = n
207  iuf = 0
208  kdflg = 1
209  ib = ib-1
210  ic = ib-1
211  DO 260 k=1,n
212  fn = fnu + float(kk-1)
213 C-----------------------------------------------------------------------
214 C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
215 C FUNCTION ABOVE
216 C-----------------------------------------------------------------------
217  m=3
218  IF (n.GT.2) GO TO 175
219  170 CONTINUE
220  initd = init(j)
221  phid = phi(j)
222  zeta1d = zeta1(j)
223  zeta2d = zeta2(j)
224  sumd = sum(j)
225  m = j
226  j = 3 - j
227  GO TO 180
228  175 CONTINUE
229  IF ((kk.EQ.n).AND.(ib.LT.n)) GO TO 180
230  IF ((kk.EQ.ib).OR.(kk.EQ.ic)) GO TO 170
231  initd = 0
232  180 CONTINUE
233  CALL cunik(zr, fn, 1, 0, tol, initd, phid, zeta1d,
234  * zeta2d, sumd, cwrk(1,m))
235  IF (kode.EQ.1) GO TO 190
236  cfn = cmplx(fn,0.0e0)
237  s1 = -zeta1d + cfn*(cfn/(zr+zeta2d))
238  GO TO 200
239  190 CONTINUE
240  s1 = -zeta1d + zeta2d
241  200 CONTINUE
242 C-----------------------------------------------------------------------
243 C TEST FOR UNDERFLOW AND OVERFLOW
244 C-----------------------------------------------------------------------
245  rs1 = REAL(S1)
246  IF (abs(rs1).GT.elim) GO TO 250
247  IF (kdflg.EQ.1) iflag = 2
248  IF (abs(rs1).LT.alim) GO TO 210
249 C-----------------------------------------------------------------------
250 C REFINE TEST AND SCALE
251 C-----------------------------------------------------------------------
252  aphi = cabs(phid)
253  rs1 = rs1 + alog(aphi)
254  IF (abs(rs1).GT.elim) GO TO 250
255  IF (kdflg.EQ.1) iflag = 1
256  IF (rs1.LT.0.0e0) GO TO 210
257  IF (kdflg.EQ.1) iflag = 3
258  210 CONTINUE
259  s2 = csgn*phid*sumd
260  c2r = REAL(S1)
261  c2i = aimag(s1)
262  c2m = exp(c2r)*REAL(CSS(IFLAG))
263  s1 = cmplx(c2m,0.0e0)*cmplx(cos(c2i),sin(c2i))
264  s2 = s2*s1
265  IF (iflag.NE.1) GO TO 220
266  CALL cuchk(s2, nw, bry(1), tol)
267  IF (nw.NE.0) s2 = cmplx(0.0e0,0.0e0)
268  220 CONTINUE
269  cy(kdflg) = s2
270  c2 = s2
271  s2 = s2*csr(iflag)
272 C-----------------------------------------------------------------------
273 C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
274 C-----------------------------------------------------------------------
275  s1 = y(kk)
276  IF (kode.EQ.1) GO TO 240
277  CALL cs1s2(zr, s1, s2, nw, asc, alim, iuf)
278  nz = nz + nw
279  240 CONTINUE
280  y(kk) = s1*cspn + s2
281  kk = kk - 1
282  cspn = -cspn
283  IF (c2.NE.czero) GO TO 245
284  kdflg = 1
285  GO TO 260
286  245 CONTINUE
287  IF (kdflg.EQ.2) GO TO 265
288  kdflg = 2
289  GO TO 260
290  250 CONTINUE
291  IF (rs1.GT.0.0e0) GO TO 290
292  s2 = czero
293  GO TO 220
294  260 CONTINUE
295  k = n
296  265 CONTINUE
297  il = n - k
298  IF (il.EQ.0) RETURN
299 C-----------------------------------------------------------------------
300 C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
301 C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
302 C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
303 C-----------------------------------------------------------------------
304  s1 = cy(1)
305  s2 = cy(2)
306  cs = csr(iflag)
307  ascle = bry(iflag)
308  fn = float(inu+il)
309  DO 280 i=1,il
310  c2 = s2
311  s2 = s1 + cmplx(fn+fnf,0.0e0)*rz*s2
312  s1 = c2
313  fn = fn - 1.0e0
314  c2 = s2*cs
315  ck = c2
316  c1 = y(kk)
317  IF (kode.EQ.1) GO TO 270
318  CALL cs1s2(zr, c1, c2, nw, asc, alim, iuf)
319  nz = nz + nw
320  270 CONTINUE
321  y(kk) = c1*cspn + c2
322  kk = kk - 1
323  cspn = -cspn
324  IF (iflag.GE.3) GO TO 280
325  c2r = REAL(CK)
326  c2i = aimag(ck)
327  c2r = abs(c2r)
328  c2i = abs(c2i)
329  c2m = amax1(c2r,c2i)
330  IF (c2m.LE.ascle) GO TO 280
331  iflag = iflag + 1
332  ascle = bry(iflag)
333  s1 = s1*cs
334  s2 = ck
335  s1 = s1*css(iflag)
336  s2 = s2*css(iflag)
337  cs = csr(iflag)
338  280 CONTINUE
339  RETURN
340  290 CONTINUE
341  nz = -1
342  RETURN
343  END
mod
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:860
abs
static T abs(T x)
Definition: pr-output.cc:1696
cmplx
double complex cmplx
Definition: Faddeeva.cc:217
dimension
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)