GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
cuoik.f
Go to the documentation of this file.
1  SUBROUTINE cuoik(Z, FNU, KODE, IKFLG, N, Y, NUF, TOL, ELIM, ALIM)
2 C***BEGIN PROLOGUE CUOIK
3 C***REFER TO CBESI,CBESK,CBESH
4 C
5 C CUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC
6 C EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM
7 C (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW
8 C WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING
9 C EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN
10 C THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER
11 C MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE
12 C EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)=
13 C EXP(-ELIM)/TOL
14 C
15 C IKFLG=1 MEANS THE I SEQUENCE IS TESTED
16 C =2 MEANS THE K SEQUENCE IS TESTED
17 C NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE
18 C =-1 MEANS AN OVERFLOW WOULD OCCUR
19 C IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO
20 C THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE
21 C IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO
22 C IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY
23 C ANOTHER ROUTINE
24 C
25 C***ROUTINES CALLED CUCHK,CUNHJ,CUNIK,R1MACH
26 C***END PROLOGUE CUOIK
27  COMPLEX ARG, ASUM, BSUM, CWRK, CZ, CZERO, PHI, SUM, Y, Z, ZB,
28  * ZETA1, ZETA2, ZN, ZR
29  REAL AARG, AIC, ALIM, APHI, ASCLE, AX, AY, ELIM, FNN, FNU, GNN,
30  * GNU, RCZ, TOL, X, YY
31  INTEGER I, IFORM, IKFLG, INIT, KODE, N, NN, NUF, NW
32  dimension y(n), cwrk(16)
33  DATA czero / (0.0e0,0.0e0) /
34  DATA aic / 1.265512123484645396e+00 /
35  nuf = 0
36  nn = n
37  x = real(z)
38  zr = z
39  IF (x.LT.0.0e0) zr = -z
40  zb = zr
41  yy = aimag(zr)
42  ax = abs(x)*1.7321e0
43  ay = abs(yy)
44  iform = 1
45  IF (ay.GT.ax) iform = 2
46  gnu = amax1(fnu,1.0e0)
47  IF (ikflg.EQ.1) GO TO 10
48  fnn = float(nn)
49  gnn = fnu + fnn - 1.0e0
50  gnu = amax1(gnn,fnn)
51  10 CONTINUE
52 C-----------------------------------------------------------------------
53 C ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE
54 C REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET
55 C THE SIGN OF THE IMAGINARY PART CORRECT.
56 C-----------------------------------------------------------------------
57  IF (iform.EQ.2) GO TO 20
58  init = 0
59  CALL cunik(zr, gnu, ikflg, 1, tol, init, phi, zeta1, zeta2, sum,
60  * cwrk)
61  cz = -zeta1 + zeta2
62  GO TO 40
63  20 CONTINUE
64  zn = -zr*cmplx(0.0e0,1.0e0)
65  IF (yy.GT.0.0e0) GO TO 30
66  zn = conjg(-zn)
67  30 CONTINUE
68  CALL cunhj(zn, gnu, 1, tol, phi, arg, zeta1, zeta2, asum, bsum)
69  cz = -zeta1 + zeta2
70  aarg = cabs(arg)
71  40 CONTINUE
72  IF (kode.EQ.2) cz = cz - zb
73  IF (ikflg.EQ.2) cz = -cz
74  aphi = cabs(phi)
75  rcz = real(cz)
76 C-----------------------------------------------------------------------
77 C OVERFLOW TEST
78 C-----------------------------------------------------------------------
79  IF (rcz.GT.elim) GO TO 170
80  IF (rcz.LT.alim) GO TO 50
81  rcz = rcz + alog(aphi)
82  IF (iform.EQ.2) rcz = rcz - 0.25e0*alog(aarg) - aic
83  IF (rcz.GT.elim) GO TO 170
84  GO TO 100
85  50 CONTINUE
86 C-----------------------------------------------------------------------
87 C UNDERFLOW TEST
88 C-----------------------------------------------------------------------
89  IF (rcz.LT.(-elim)) GO TO 60
90  IF (rcz.GT.(-alim)) GO TO 100
91  rcz = rcz + alog(aphi)
92  IF (iform.EQ.2) rcz = rcz - 0.25e0*alog(aarg) - aic
93  IF (rcz.GT.(-elim)) GO TO 80
94  60 CONTINUE
95  DO 70 i=1,nn
96  y(i) = czero
97  70 CONTINUE
98  nuf = nn
99  RETURN
100  80 CONTINUE
101  ascle = 1.0e+3*r1mach(1)/tol
102  cz = cz + clog(phi)
103  IF (iform.EQ.1) GO TO 90
104  cz = cz - cmplx(0.25e0,0.0e0)*clog(arg) - cmplx(aic,0.0e0)
105  90 CONTINUE
106  ax = exp(rcz)/tol
107  ay = aimag(cz)
108  cz = cmplx(ax,0.0e0)*cmplx(cos(ay),sin(ay))
109  CALL cuchk(cz, nw, ascle, tol)
110  IF (nw.EQ.1) GO TO 60
111  100 CONTINUE
112  IF (ikflg.EQ.2) RETURN
113  IF (n.EQ.1) RETURN
114 C-----------------------------------------------------------------------
115 C SET UNDERFLOWS ON I SEQUENCE
116 C-----------------------------------------------------------------------
117  110 CONTINUE
118  gnu = fnu + float(nn-1)
119  IF (iform.EQ.2) GO TO 120
120  init = 0
121  CALL cunik(zr, gnu, ikflg, 1, tol, init, phi, zeta1, zeta2, sum,
122  * cwrk)
123  cz = -zeta1 + zeta2
124  GO TO 130
125  120 CONTINUE
126  CALL cunhj(zn, gnu, 1, tol, phi, arg, zeta1, zeta2, asum, bsum)
127  cz = -zeta1 + zeta2
128  aarg = cabs(arg)
129  130 CONTINUE
130  IF (kode.EQ.2) cz = cz - zb
131  aphi = cabs(phi)
132  rcz = real(cz)
133  IF (rcz.LT.(-elim)) GO TO 140
134  IF (rcz.GT.(-alim)) RETURN
135  rcz = rcz + alog(aphi)
136  IF (iform.EQ.2) rcz = rcz - 0.25e0*alog(aarg) - aic
137  IF (rcz.GT.(-elim)) GO TO 150
138  140 CONTINUE
139  y(nn) = czero
140  nn = nn - 1
141  nuf = nuf + 1
142  IF (nn.EQ.0) RETURN
143  GO TO 110
144  150 CONTINUE
145  ascle = 1.0e+3*r1mach(1)/tol
146  cz = cz + clog(phi)
147  IF (iform.EQ.1) GO TO 160
148  cz = cz - cmplx(0.25e0,0.0e0)*clog(arg) - cmplx(aic,0.0e0)
149  160 CONTINUE
150  ax = exp(rcz)/tol
151  ay = aimag(cz)
152  cz = cmplx(ax,0.0e0)*cmplx(cos(ay),sin(ay))
153  CALL cuchk(cz, nw, ascle, tol)
154  IF (nw.EQ.1) GO TO 140
155  RETURN
156  170 CONTINUE
157  nuf = -1
158  RETURN
159  END
double complex cmplx
Definition: Faddeeva.cc:230
subroutine cuchk(Y, NZ, ASCLE, TOL)
Definition: cuchk.f:2
subroutine cunhj(Z, FNU, IPMTR, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM)
Definition: cunhj.f:3
subroutine cunik(ZR, FNU, IKFLG, IPMTR, TOL, INIT, PHI, ZETA1, ZETA2, SUM, CWRK)
Definition: cunik.f:3
subroutine cuoik(Z, FNU, KODE, IKFLG, N, Y, NUF, TOL, ELIM, ALIM)
Definition: cuoik.f:2
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137
real function r1mach(i)
Definition: r1mach.f:23