GNU Octave  4.0.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
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
cairy
subroutine cairy(Z, ID, KODE, AI, NZ, IERR)
Definition: cairy.f:1
workspace_element::dimension
std::string dimension(void) const
Definition: workspace-element.h:74
mod
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:959
base_det::exp
int exp(void) const
Definition: DET.h:64
cuchk
subroutine cuchk(Y, NZ, ASCLE, TOL)
Definition: cuchk.f:1
cunk2
subroutine cunk2(Z, FNU, KODE, MR, N, Y, NZ, TOL, ELIM, ALIM)
Definition: cunk2.f:1
cmplx
double complex cmplx
Definition: Faddeeva.cc:226
octave_value::sin
octave_value sin(void) const
Definition: ov.h:1198
cs1s2
subroutine cs1s2(ZR, S1, S2, NZ, ASCLE, ALIM, IUF)
Definition: cs1s2.f:1
cunhj
subroutine cunhj(Z, FNU, IPMTR, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM)
Definition: cunhj.f:1
octave_value::cos
octave_value cos(void) const
Definition: ov.h:1170
abs
T abs(T x)
Definition: pr-output.cc:3062
Generated on Thu Jun 4 2015 23:30:23 for GNU Octave by   doxygen 1.8.8