1 SUBROUTINE zbiry(ZR, ZI, ID, KODE, BIR, BII, IERR)
124 DOUBLE PRECISION AA, AD, AK, ALIM, ATRM, AZ, AZ3, BB, BII, BIR,
125 * BK, CC, CK, COEF, CONEI, CONER, CSQI, CSQR, CYI, CYR, C1, C2,
126 * DIG, DK, D1, D2, EAA, ELIM, FID, FMR, FNU, FNUL, PI, RL, R1M5,
127 * SFAC, STI, STR, S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I,
128 * TRM2R, TTH, ZI, ZR, ZTAI, ZTAR, Z3I, Z3R, D1MACH, XZABS
129 INTEGER ID, IERR, K, KODE, K1, K2, NZ, I1MACH
131 DATA tth, c1, c2, coef, pi /6.66666666666666667d-01,
132 * 6.14926627446000736d-01,4.48288357353826359d-01,
133 * 5.77350269189625765d-01,3.14159265358979324d+00/
134 DATA coner, conei /1.0d0,0.0d0/
138 IF (id.LT.0 .OR. id.GT.1) ierr=1
139 IF (kode.LT.1 .OR. kode.GT.2) ierr=1
140 IF (ierr.NE.0)
RETURN 142 tol = dmax1(d1mach(4),1.0d-18)
143 fid = dble(float(id))
144 IF (az.GT.1.0e0)
GO TO 70
152 IF (az.LT.tol)
GO TO 130
154 IF (aa.LT.tol/az)
GO TO 40
162 z3r = str*zr - sti*zi
163 z3i = str*zi + sti*zr
166 bk = 3.0d0 - fid - fid
168 dk = 3.0d0 + fid + fid
172 ak = 24.0d0 + 9.0d0*fid
173 bk = 30.0d0 - 9.0d0*fid
175 str = (trm1r*z3r-trm1i*z3i)/d1
176 trm1i = (trm1r*z3i+trm1i*z3r)/d1
180 str = (trm2r*z3r-trm2i*z3i)/d2
181 trm2i = (trm2r*z3i+trm2i*z3r)/d2
189 IF (atrm.LT.tol*ad)
GO TO 40
194 IF (id.EQ.1)
GO TO 50
195 bir = c1*s1r + c2*(zr*s2r-zi*s2i)
196 bii = c1*s1i + c2*(zr*s2i+zi*s2r)
197 IF (kode.EQ.1)
RETURN 198 CALL xzsqrt(zr, zi, str, sti)
199 ztar = tth*(zr*str-zi*sti)
200 ztai = tth*(zr*sti+zi*str)
210 IF (az.LE.tol)
GO TO 60
212 str = s1r*zr - s1i*zi
213 sti = s1r*zi + s1i*zr
214 bir = bir + cc*(str*zr-sti*zi)
215 bii = bii + cc*(str*zi+sti*zr)
217 IF (kode.EQ.1)
RETURN 218 CALL xzsqrt(zr, zi, str, sti)
219 ztar = tth*(zr*str-zi*sti)
220 ztai = tth*(zr*sti+zi*str)
231 fnu = (1.0d0+fid)/3.0d0
246 k = min0(iabs(k1),iabs(k2))
247 elim = 2.303d0*(dble(float(k))*r1m5-3.0d0)
249 aa = r1m5*dble(float(k1))
250 dig = dmin1(aa,18.0d0)
252 alim = elim + dmax1(-aa,-41.45d0)
253 rl = 1.2d0*dig + 3.0d0
254 fnul = 10.0d0 + 6.0d0*(dig-3.0d0)
259 bb=dble(float(i1mach(9)))*0.5d0
262 IF (az.GT.aa)
GO TO 260
265 CALL xzsqrt(zr, zi, csqr, csqi)
266 ztar = tth*(zr*csqr-zi*csqi)
267 ztai = tth*(zr*csqi+zi*csqr)
273 IF (zr.GE.0.0d0)
GO TO 80
279 IF (zi.NE.0.0d0 .OR. zr.GT.0.0d0)
GO TO 90
284 IF (kode.EQ.2)
GO TO 100
289 IF (bb.LT.alim)
GO TO 100
290 bb = bb + 0.25d0*dlog(az)
292 IF (bb.GT.elim)
GO TO 190
295 IF (aa.GE.0.0d0 .AND. zr.GT.0.0d0)
GO TO 110
297 IF (zi.LT.0.0d0) fmr = -pi
305 CALL zbinu(ztar, ztai, fnu, kode, 1, cyr, cyi, nz, rl, fnul, tol,
307 IF (nz.LT.0)
GO TO 200
312 s1r = (str*cyr(1)-sti*cyi(1))*z3r
313 s1i = (str*cyi(1)+sti*cyr(1))*z3r
314 fnu = (2.0d0-fid)/3.0d0
315 CALL zbinu(ztar, ztai, fnu, kode, 2, cyr, cyi, nz, rl, fnul, tol,
324 CALL zdiv(cyr(1), cyi(1), ztar, ztai, str, sti)
325 s2r = (fnu+fnu)*str + cyr(2)
326 s2i = (fnu+fnu)*sti + cyi(2)
330 s1r = coef*(s1r+s2r*str-s2i*sti)
331 s1i = coef*(s1i+s2r*sti+s2i*str)
332 IF (id.EQ.1)
GO TO 120
333 str = csqr*s1r - csqi*s1i
334 s1i = csqr*s1i + csqi*s1r
340 str = zr*s1r - zi*s1i
341 s1i = zr*s1i + zi*s1r
347 aa = c1*(1.0d0-fid) + fid*c2
356 IF(nz.EQ.(-1))
GO TO 190
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)