GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
cunk2.f
Go to the documentation of this file.
1 SUBROUTINE cunk2(Z, FNU, KODE, MR, N, Y, NZ, TOL, ELIM, ALIM)
2C***BEGIN PROLOGUE CUNK2
3C***REFER TO CBESK
4C
5C CUNK2 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
6C RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
7C UNIFORM ASYMPTOTIC EXPANSIONS FOR H(KIND,FNU,ZN) AND J(FNU,ZN)
8C WHERE ZN IS IN THE RIGHT HALF PLANE, KIND=(3-MR)/2, MR=+1 OR
9C -1. HERE ZN=ZR*I OR -ZR*I WHERE ZR=Z IF Z IS IN THE RIGHT
10C HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC-
11C ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
12C NZ=-1 MEANS AN OVERFLOW WILL OCCUR
13C
14C***ROUTINES CALLED CAIRY,CS1S2,CUCHK,CUNHJ,R1MACH
15C***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)/
35C
36 kdflg = 1
37 nz = 0
38C-----------------------------------------------------------------------
39C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
40C THE UNDERFLOW LIMIT
41C-----------------------------------------------------------------------
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
73C-----------------------------------------------------------------------
74C K(FNU,Z) IS COMPUTED FROM H(2,FNU,-I*Z) WHERE Z IS IN THE FIRST
75C QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
76C CONJUGATION SINCE THE K FUNCTION IS REAL ON THE POSITIVE REAL AXIS
77C-----------------------------------------------------------------------
78 j = 2
79 DO 70 i=1,n
80C-----------------------------------------------------------------------
81C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
82C-----------------------------------------------------------------------
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
94C-----------------------------------------------------------------------
95C TEST FOR UNDERFLOW AND OVERFLOW
96C-----------------------------------------------------------------------
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
101C-----------------------------------------------------------------------
102C REFINE TEST AND SCALE
103C-----------------------------------------------------------------------
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
112C-----------------------------------------------------------------------
113C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
114C EXPONENT EXTREMES
115C-----------------------------------------------------------------------
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
138C-----------------------------------------------------------------------
139C FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
140C-----------------------------------------------------------------------
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
157C-----------------------------------------------------------------------
158C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW, SET SEQUENCE TO ZERO
159C ON UNDERFLOW
160C-----------------------------------------------------------------------
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
175C-----------------------------------------------------------------------
176C REFINE ESTIMATE AND TEST
177C-----------------------------------------------------------------------
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
184C-----------------------------------------------------------------------
185C FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
186C-----------------------------------------------------------------------
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
194C-----------------------------------------------------------------------
195C SCALED FORWARD RECURRENCE FOR REMAINDER OF THE SEQUENCE
196C-----------------------------------------------------------------------
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
225C-----------------------------------------------------------------------
226C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0E0
227C-----------------------------------------------------------------------
228 nz = 0
229 fmr = float(mr)
230 sgn = -sign(pi,fmr)
231C-----------------------------------------------------------------------
232C CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
233C-----------------------------------------------------------------------
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
242C-----------------------------------------------------------------------
243C CS=COEFF OF THE J FUNCTION TO GET THE I FUNCTION. I(FNU,Z) IS
244C COMPUTED FROM EXP(I*FNU*HPI)*J(FNU,-I*Z) WHERE Z IS IN THE FIRST
245C QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
246C CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS
247C-----------------------------------------------------------------------
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
259C-----------------------------------------------------------------------
260C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
261C FUNCTION ABOVE
262C-----------------------------------------------------------------------
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
287C-----------------------------------------------------------------------
288C TEST FOR UNDERFLOW AND OVERFLOW
289C-----------------------------------------------------------------------
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
294C-----------------------------------------------------------------------
295C REFINE TEST AND SCALE
296C-----------------------------------------------------------------------
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)
321C-----------------------------------------------------------------------
322C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
323C-----------------------------------------------------------------------
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
349C-----------------------------------------------------------------------
350C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
351C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
352C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
353C-----------------------------------------------------------------------
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
cmplx
double complex cmplx
Definition Faddeeva.cc:227
cairy
subroutine cairy(z, id, kode, ai, nz, ierr)
Definition cairy.f:2
cs1s2
subroutine cs1s2(zr, s1, s2, nz, ascle, alim, iuf)
Definition cs1s2.f:2
cuchk
subroutine cuchk(y, nz, ascle, tol)
Definition cuchk.f:2
cunhj
subroutine cunhj(z, fnu, ipmtr, tol, phi, arg, zeta1, zeta2, asum, bsum)
Definition cunhj.f:3
cunk2
subroutine cunk2(z, fnu, kode, mr, n, y, nz, tol, elim, alim)
Definition cunk2.f:2
real
ColumnVector real(const ComplexColumnVector &a)
Definition dColVector.cc:137
mod
T mod(T x, T y)
Definition lo-mappers.h:294