1 SUBROUTINE cbiry(Z, ID, KODE, BI, IERR)
120 COMPLEX BI, CONE, CSQ, CY, S1, S2, TRM1, TRM2, Z, ZTA, Z3
121 REAL AA, AD, AK, ALIM, ATRM, AZ, AZ3, BB, BK, CK, COEF, C1, C2,
122 * dig, dk, d1, d2, elim, fid, fmr, fnu, fnul, pi, rl, r1m5, sfac,
123 * tol, tth, zi, zr, z3i, z3r, r1mach
124 INTEGER ID, IERR, K, KODE, K1, K2, NZ, I1MACH
126 DATA tth, c1, c2, coef, pi /6.66666666666666667e-01,
127 * 6.14926627446000736e-01,4.48288357353826359e-01,
128 * 5.77350269189625765e-01,3.14159265358979324e+00/
129 DATA cone / (1.0e0,0.0e0) /
133 IF (id.LT.0 .OR. id.GT.1) ierr=1
134 IF (kode.LT.1 .OR. kode.GT.2) ierr=1
135 IF (ierr.NE.0)
RETURN
137 tol = amax1(r1mach(4),1.0e-18)
139 IF (az.GT.1.0e0) go to 60
145 IF (az.LT.tol) go to 110
147 IF (aa.LT.tol/az) go to 40
154 bk = 3.0e0 - fid - fid
156 dk = 3.0e0 + fid + fid
160 ak = 24.0e0 + 9.0e0*fid
161 bk = 30.0e0 - 9.0e0*fid
165 trm1 = trm1*
cmplx(z3r/d1,z3i/d1)
167 trm2 = trm2*
cmplx(z3r/d2,z3i/d2)
173 IF (atrm.LT.tol*ad) go to 40
178 IF (id.EQ.1) go to 50
179 bi = s1*
cmplx(c1,0.0e0) + z*s2*
cmplx(c2,0.0e0)
180 IF (kode.EQ.1)
RETURN
181 zta = z*csqrt(z)*
cmplx(tth,0.0e0)
187 bi = s2*
cmplx(c2,0.0e0)
188 IF (az.GT.tol) bi = bi + z*z*s1*
cmplx(c1/(1.0e0+fid),0.0e0)
189 IF (kode.EQ.1)
RETURN
190 zta = z*csqrt(z)*
cmplx(tth,0.0e0)
199 fnu = (1.0e0+fid)/3.0e0
214 k = min0(iabs(k1),iabs(k2))
215 elim = 2.303e0*(float(k)*r1m5-3.0e0)
218 dig = amin1(aa,18.0e0)
220 alim = elim + amax1(-aa,-41.45e0)
221 rl = 1.2e0*dig + 3.0e0
222 fnul = 10.0e0 + 6.0e0*(dig-3.0e0)
227 bb=float(i1mach(9))*0.5e0
230 IF (az.GT.aa) go to 190
234 zta=z*csq*
cmplx(tth,0.0e0)
242 IF (zr.GE.0.0e0) go to 70
247 IF (zi.EQ.0.0e0 .AND. zr.LE.0.0e0) zta =
cmplx(0.0e0,ak)
249 IF (kode.EQ.2) go to 80
254 IF (bb.LT.alim) go to 80
255 bb = bb + 0.25e0*alog(az)
257 IF (bb.GT.elim) go to 170
260 IF (aa.GE.0.0e0 .AND. zr.GT.0.0e0) go to 90
262 IF (zi.LT.0.0e0) fmr = -pi
269 CALL
cbinu(zta, fnu, kode, 1, cy, nz, rl, fnul, tol, elim, alim)
270 IF (nz.LT.0) go to 180
272 z3 =
cmplx(sfac,0.0e0)
274 fnu = (2.0e0-fid)/3.0e0
275 CALL
cbinu(zta, fnu, kode, 2, cy, nz, rl, fnul, tol, elim, alim)
281 s2 = cy(1)*
cmplx(fnu+fnu,0.0e0)/zta + cy(2)
284 IF (id.EQ.1) go to 100
286 bi = s1*
cmplx(1.0e0/sfac,0.0e0)
290 bi = s1*
cmplx(1.0e0/sfac,0.0e0)
293 aa = c1*(1.0e0-fid) + fid*c2
301 IF(nz.EQ.(-1)) go to 170
std::string dimension(void) const
subroutine cbiry(Z, ID, KODE, BI, IERR)
octave_value sin(void) const
octave_value cos(void) const
octave_value sqrt(void) const
subroutine cbinu(Z, FNU, KODE, N, CY, NZ, RL, FNUL, TOL, ELIM, ALIM)