GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
cunk2.f
Go to the documentation of this file.
1  SUBROUTINE cunk2(Z, FNU, KODE, MR, N, Y, NZ, TOL, ELIM, ALIM)
2 C***BEGIN PROLOGUE CUNK2
3 C***REFER TO CBESK
4 C
5 C CUNK2 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 EXPANSIONS FOR H(KIND,FNU,ZN) AND J(FNU,ZN)
8 C WHERE ZN IS IN THE RIGHT HALF PLANE, KIND=(3-MR)/2, MR=+1 OR
9 C -1. HERE ZN=ZR*I OR -ZR*I WHERE ZR=Z IF Z IS IN THE RIGHT
10 C HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC-
11 C ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
12 C NZ=-1 MEANS AN OVERFLOW WILL OCCUR
13 C
14 C***ROUTINES CALLED CAIRY,CS1S2,CUCHK,CUNHJ,R1MACH
15 C***END PROLOGUE CUNK2
16  COMPLEX AI, ARG, ASUM, BSUM, CFN, CI, CIP,
17  * CK, CONE, CRSC, CR1, CR2, CS, CSCL, CSGN, CSPN, CSR, CSS, CY,
18  * CZERO, C1, C2, DAI, PHI, RZ, S1, S2, Y, Z, ZB, ZETA1,
19  * ZETA2, ZN, ZR, PHID, ARGD, ZETA1D, ZETA2D, ASUMD, BSUMD
20  REAL AARG, AIC, ALIM, ANG, APHI, ASC, ASCLE, BRY, CAR, CPN, C2I,
21  * C2M, C2R, ELIM, FMR, FN, FNF, FNU, HPI, PI, RS1, SAR, SGN, SPN,
22  * TOL, X, YY, R1MACH
23  INTEGER I, IB, IFLAG, IFN, IL, IN, INU, IUF, K, KDFLG, KFLAG, KK,
24  * KODE, MR, N, NAI, NDAI, NW, NZ, IDUM, J, IPARD, IC
25  dimension bry(3), y(n), asum(2), bsum(2), phi(2), arg(2),
26  * zeta1(2), zeta2(2), cy(2), cip(4), css(3), csr(3)
27  DATA czero, cone, ci, cr1, cr2 /
28  1 (0.0e0,0.0e0),(1.0e0,0.0e0),(0.0e0,1.0e0),
29  1(1.0e0,1.73205080756887729e0),(-0.5e0,-8.66025403784438647e-01)/
30  DATA hpi, pi, aic /
31  1 1.57079632679489662e+00, 3.14159265358979324e+00,
32  1 1.26551212348464539e+00/
33  DATA cip(1),cip(2),cip(3),cip(4)/
34  1 (1.0e0,0.0e0), (0.0e0,-1.0e0), (-1.0e0,0.0e0), (0.0e0,1.0e0)/
35 C
36  kdflg = 1
37  nz = 0
38 C-----------------------------------------------------------------------
39 C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
40 C THE UNDERFLOW LIMIT
41 C-----------------------------------------------------------------------
42  cscl = cmplx(1.0e0/tol,0.0e0)
43  crsc = cmplx(tol,0.0e0)
44  css(1) = cscl
45  css(2) = cone
46  css(3) = crsc
47  csr(1) = crsc
48  csr(2) = cone
49  csr(3) = cscl
50  bry(1) = 1.0e+3*r1mach(1)/tol
51  bry(2) = 1.0e0/bry(1)
52  bry(3) = r1mach(2)
53  x = REAL(Z)
54  zr = z
55  IF (x.LT.0.0e0) zr = -z
56  yy = aimag(zr)
57  zn = -zr*ci
58  zb = zr
59  inu = int(fnu)
60  fnf = fnu - float(inu)
61  ang = -hpi*fnf
62  car = cos(ang)
63  sar = sin(ang)
64  cpn = -hpi*car
65  spn = -hpi*sar
66  c2 = cmplx(-spn,cpn)
67  kk = mod(inu,4) + 1
68  cs = cr1*c2*cip(kk)
69  IF (yy.GT.0.0e0) GO TO 10
70  zn = conjg(-zn)
71  zb = conjg(zb)
72  10 CONTINUE
73 C-----------------------------------------------------------------------
74 C K(FNU,Z) IS COMPUTED FROM H(2,FNU,-I*Z) WHERE Z IS IN THE FIRST
75 C QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
76 C CONJUGATION SINCE THE K FUNCTION IS REAL ON THE POSITIVE REAL AXIS
77 C-----------------------------------------------------------------------
78  j = 2
79  DO 70 i=1,n
80 C-----------------------------------------------------------------------
81 C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
82 C-----------------------------------------------------------------------
83  j = 3 - j
84  fn = fnu + float(i-1)
85  CALL cunhj(zn, fn, 0, tol, phi(j), arg(j), zeta1(j), zeta2(j),
86  * asum(j), bsum(j))
87  IF (kode.EQ.1) GO TO 20
88  cfn = cmplx(fn,0.0e0)
89  s1 = zeta1(j) - cfn*(cfn/(zb+zeta2(j)))
90  GO TO 30
91  20 CONTINUE
92  s1 = zeta1(j) - zeta2(j)
93  30 CONTINUE
94 C-----------------------------------------------------------------------
95 C TEST FOR UNDERFLOW AND OVERFLOW
96 C-----------------------------------------------------------------------
97  rs1 = REAL(S1)
98  IF (abs(rs1).GT.elim) GO TO 60
99  IF (kdflg.EQ.1) kflag = 2
100  IF (abs(rs1).LT.alim) GO TO 40
101 C-----------------------------------------------------------------------
102 C REFINE TEST AND SCALE
103 C-----------------------------------------------------------------------
104  aphi = cabs(phi(j))
105  aarg = cabs(arg(j))
106  rs1 = rs1 + alog(aphi) - 0.25e0*alog(aarg) - aic
107  IF (abs(rs1).GT.elim) GO TO 60
108  IF (kdflg.EQ.1) kflag = 1
109  IF (rs1.LT.0.0e0) GO TO 40
110  IF (kdflg.EQ.1) kflag = 3
111  40 CONTINUE
112 C-----------------------------------------------------------------------
113 C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
114 C EXPONENT EXTREMES
115 C-----------------------------------------------------------------------
116  c2 = arg(j)*cr2
117  CALL cairy(c2, 0, 2, ai, nai, idum)
118  CALL cairy(c2, 1, 2, dai, ndai, idum)
119  s2 = cs*phi(j)*(ai*asum(j)+cr2*dai*bsum(j))
120  c2r = REAL(S1)
121  c2i = aimag(s1)
122  c2m = exp(c2r)*REAL(CSS(KFLAG))
123  s1 = cmplx(c2m,0.0e0)*cmplx(cos(c2i),sin(c2i))
124  s2 = s2*s1
125  IF (kflag.NE.1) GO TO 50
126  CALL cuchk(s2, nw, bry(1), tol)
127  IF (nw.NE.0) GO TO 60
128  50 CONTINUE
129  IF (yy.LE.0.0e0) s2 = conjg(s2)
130  cy(kdflg) = s2
131  y(i) = s2*csr(kflag)
132  cs = -ci*cs
133  IF (kdflg.EQ.2) GO TO 75
134  kdflg = 2
135  GO TO 70
136  60 CONTINUE
137  IF (rs1.GT.0.0e0) GO TO 300
138 C-----------------------------------------------------------------------
139 C FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
140 C-----------------------------------------------------------------------
141  IF (x.LT.0.0e0) GO TO 300
142  kdflg = 1
143  y(i) = czero
144  cs = -ci*cs
145  nz=nz+1
146  IF (i.EQ.1) GO TO 70
147  IF (y(i-1).EQ.czero) GO TO 70
148  y(i-1) = czero
149  nz=nz+1
150  70 CONTINUE
151  i=n
152  75 CONTINUE
153  rz = cmplx(2.0e0,0.0e0)/zr
154  ck = cmplx(fn,0.0e0)*rz
155  ib = i + 1
156  IF (n.LT.ib) GO TO 170
157 C-----------------------------------------------------------------------
158 C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW, SET SEQUENCE TO ZERO
159 C ON UNDERFLOW
160 C-----------------------------------------------------------------------
161  fn = fnu+float(n-1)
162  ipard = 1
163  IF (mr.NE.0) ipard = 0
164  CALL cunhj(zn,fn,ipard,tol,phid,argd,zeta1d,zeta2d,asumd,bsumd)
165  IF (kode.EQ.1) GO TO 80
166  cfn=cmplx(fn,0.0e0)
167  s1=zeta1d-cfn*(cfn/(zb+zeta2d))
168  GO TO 90
169  80 CONTINUE
170  s1=zeta1d-zeta2d
171  90 CONTINUE
172  rs1=REAL(S1)
173  IF (abs(rs1).GT.elim) GO TO 95
174  IF (abs(rs1).LT.alim) GO TO 100
175 C-----------------------------------------------------------------------
176 C REFINE ESTIMATE AND TEST
177 C-----------------------------------------------------------------------
178  aphi=cabs(phid)
179  aarg = cabs(argd)
180  rs1=rs1+alog(aphi)-0.25e0*alog(aarg)-aic
181  IF (abs(rs1).LT.elim) GO TO 100
182  95 CONTINUE
183  IF (rs1.GT.0.0e0) GO TO 300
184 C-----------------------------------------------------------------------
185 C FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
186 C-----------------------------------------------------------------------
187  IF (x.LT.0.0e0) GO TO 300
188  nz=n
189  DO 96 i=1,n
190  y(i) = czero
191  96 CONTINUE
192  RETURN
193  100 CONTINUE
194 C-----------------------------------------------------------------------
195 C SCALED FORWARD RECURRENCE FOR REMAINDER OF THE SEQUENCE
196 C-----------------------------------------------------------------------
197  s1 = cy(1)
198  s2 = cy(2)
199  c1 = csr(kflag)
200  ascle = bry(kflag)
201  DO 120 i=ib,n
202  c2 = s2
203  s2 = ck*s2 + s1
204  s1 = c2
205  ck = ck + rz
206  c2 = s2*c1
207  y(i) = c2
208  IF (kflag.GE.3) GO TO 120
209  c2r = REAL(C2)
210  c2i = aimag(c2)
211  c2r = abs(c2r)
212  c2i = abs(c2i)
213  c2m = amax1(c2r,c2i)
214  IF (c2m.LE.ascle) GO TO 120
215  kflag = kflag + 1
216  ascle = bry(kflag)
217  s1 = s1*c1
218  s2 = c2
219  s1 = s1*css(kflag)
220  s2 = s2*css(kflag)
221  c1 = csr(kflag)
222  120 CONTINUE
223  170 CONTINUE
224  IF (mr.EQ.0) RETURN
225 C-----------------------------------------------------------------------
226 C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0E0
227 C-----------------------------------------------------------------------
228  nz = 0
229  fmr = float(mr)
230  sgn = -sign(pi,fmr)
231 C-----------------------------------------------------------------------
232 C CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
233 C-----------------------------------------------------------------------
234  csgn = cmplx(0.0e0,sgn)
235  IF (yy.LE.0.0e0) csgn = conjg(csgn)
236  ifn = inu + n - 1
237  ang = fnf*sgn
238  cpn = cos(ang)
239  spn = sin(ang)
240  cspn = cmplx(cpn,spn)
241  IF (mod(ifn,2).EQ.1) cspn = -cspn
242 C-----------------------------------------------------------------------
243 C CS=COEFF OF THE J FUNCTION TO GET THE I FUNCTION. I(FNU,Z) IS
244 C COMPUTED FROM EXP(I*FNU*HPI)*J(FNU,-I*Z) WHERE Z IS IN THE FIRST
245 C QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
246 C CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS
247 C-----------------------------------------------------------------------
248  cs = cmplx(car,-sar)*csgn
249  in = mod(ifn,4) + 1
250  c2 = cip(in)
251  cs = cs*conjg(c2)
252  asc = bry(1)
253  kk = n
254  kdflg = 1
255  ib = ib-1
256  ic = ib-1
257  iuf = 0
258  DO 270 k=1,n
259 C-----------------------------------------------------------------------
260 C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
261 C FUNCTION ABOVE
262 C-----------------------------------------------------------------------
263  fn = fnu+float(kk-1)
264  IF (n.GT.2) GO TO 180
265  175 CONTINUE
266  phid = phi(j)
267  argd = arg(j)
268  zeta1d = zeta1(j)
269  zeta2d = zeta2(j)
270  asumd = asum(j)
271  bsumd = bsum(j)
272  j = 3 - j
273  GO TO 190
274  180 CONTINUE
275  IF ((kk.EQ.n).AND.(ib.LT.n)) GO TO 190
276  IF ((kk.EQ.ib).OR.(kk.EQ.ic)) GO TO 175
277  CALL cunhj(zn, fn, 0, tol, phid, argd, zeta1d, zeta2d,
278  * asumd, bsumd)
279  190 CONTINUE
280  IF (kode.EQ.1) GO TO 200
281  cfn = cmplx(fn,0.0e0)
282  s1 = -zeta1d + cfn*(cfn/(zb+zeta2d))
283  GO TO 210
284  200 CONTINUE
285  s1 = -zeta1d + zeta2d
286  210 CONTINUE
287 C-----------------------------------------------------------------------
288 C TEST FOR UNDERFLOW AND OVERFLOW
289 C-----------------------------------------------------------------------
290  rs1 = REAL(S1)
291  IF (abs(rs1).GT.elim) GO TO 260
292  IF (kdflg.EQ.1) iflag = 2
293  IF (abs(rs1).LT.alim) GO TO 220
294 C-----------------------------------------------------------------------
295 C REFINE TEST AND SCALE
296 C-----------------------------------------------------------------------
297  aphi = cabs(phid)
298  aarg = cabs(argd)
299  rs1 = rs1 + alog(aphi) - 0.25e0*alog(aarg) - aic
300  IF (abs(rs1).GT.elim) GO TO 260
301  IF (kdflg.EQ.1) iflag = 1
302  IF (rs1.LT.0.0e0) GO TO 220
303  IF (kdflg.EQ.1) iflag = 3
304  220 CONTINUE
305  CALL cairy(argd, 0, 2, ai, nai, idum)
306  CALL cairy(argd, 1, 2, dai, ndai, idum)
307  s2 = cs*phid*(ai*asumd+dai*bsumd)
308  c2r = REAL(S1)
309  c2i = aimag(s1)
310  c2m = exp(c2r)*REAL(CSS(IFLAG))
311  s1 = cmplx(c2m,0.0e0)*cmplx(cos(c2i),sin(c2i))
312  s2 = s2*s1
313  IF (iflag.NE.1) GO TO 230
314  CALL cuchk(s2, nw, bry(1), tol)
315  IF (nw.NE.0) s2 = cmplx(0.0e0,0.0e0)
316  230 CONTINUE
317  IF (yy.LE.0.0e0) s2 = conjg(s2)
318  cy(kdflg) = s2
319  c2 = s2
320  s2 = s2*csr(iflag)
321 C-----------------------------------------------------------------------
322 C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
323 C-----------------------------------------------------------------------
324  s1 = y(kk)
325  IF (kode.EQ.1) GO TO 250
326  CALL cs1s2(zr, s1, s2, nw, asc, alim, iuf)
327  nz = nz + nw
328  250 CONTINUE
329  y(kk) = s1*cspn + s2
330  kk = kk - 1
331  cspn = -cspn
332  cs = -cs*ci
333  IF (c2.NE.czero) GO TO 255
334  kdflg = 1
335  GO TO 270
336  255 CONTINUE
337  IF (kdflg.EQ.2) GO TO 275
338  kdflg = 2
339  GO TO 270
340  260 CONTINUE
341  IF (rs1.GT.0.0e0) GO TO 300
342  s2 = czero
343  GO TO 230
344  270 CONTINUE
345  k = n
346  275 CONTINUE
347  il = n-k
348  IF (il.EQ.0) RETURN
349 C-----------------------------------------------------------------------
350 C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
351 C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
352 C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
353 C-----------------------------------------------------------------------
354  s1 = cy(1)
355  s2 = cy(2)
356  cs = csr(iflag)
357  ascle = bry(iflag)
358  fn = float(inu+il)
359  DO 290 i=1,il
360  c2 = s2
361  s2 = s1 + cmplx(fn+fnf,0.0e0)*rz*s2
362  s1 = c2
363  fn = fn - 1.0e0
364  c2 = s2*cs
365  ck = c2
366  c1 = y(kk)
367  IF (kode.EQ.1) GO TO 280
368  CALL cs1s2(zr, c1, c2, nw, asc, alim, iuf)
369  nz = nz + nw
370  280 CONTINUE
371  y(kk) = c1*cspn + c2
372  kk = kk - 1
373  cspn = -cspn
374  IF (iflag.GE.3) GO TO 290
375  c2r = REAL(CK)
376  c2i = aimag(ck)
377  c2r = abs(c2r)
378  c2i = abs(c2i)
379  c2m = amax1(c2r,c2i)
380  IF (c2m.LE.ascle) GO TO 290
381  iflag = iflag + 1
382  ascle = bry(iflag)
383  s1 = s1*cs
384  s2 = ck
385  s1 = s1*css(iflag)
386  s2 = s2*css(iflag)
387  cs = csr(iflag)
388  290 CONTINUE
389  RETURN
390  300 CONTINUE
391  nz = -1
392  RETURN
393  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)