Nothing
#### At first, simple interfaces to the main subroutines
#### in ../src/zbsubs.f
#### ~~~~~~~~~~~~~~~
## 10: zbesh(zr, zi, fnu, kode, m, n, cyr, cyi, nz, ierr)
## 360: zbesi(zr, zi, fnu, kode, n, cyr, cyi, nz, ierr)
## 631: zbesj(zr, zi, fnu, kode, n, cyr, cyi, nz, ierr)
## 899: zbesk(zr, zi, fnu, kode, n, cyr, cyi, nz, ierr)
## 1182: zbesy(zr, zi, fnu, kode, n, cyr, cyi, nz, cwrkr, cwrki, ierr)
##
## 1467: zairy(zr, zi, id, kode, air, aii, nz, ierr)
## 1862: zbiry(zr, zi, id, kode, bir, bii, ierr)
## 360: zbesi(zr, zi, fnu, kode, n, cyr, cyi, nz, ierr)
BesselI <- function(z, nu, expon.scaled = FALSE, nSeq = 1, verbose = 0)
{
nz <- length(z)
if(nz == 0) return(z)
isNum <- is.numeric(z)
if(isNum) {
zr <- as.double(z)
zi <- numeric(nz)
} else if(is.complex(z)) {
zr <- Re(z)
zi <- Im(z)
} else stop("'z' must be complex or numeric")
nu <- as.double(nu)
stopifnot(length(nu) == 1, length(verbose) == 1,
nSeq >= 1, nSeq == as.integer(nSeq))
if(nu < 0) {
## I(-nu,z) = I(nu,z) + (2/pi)*sin(pi*nu)*K(nu,z) [ = A.& S. 9.6.2, p.375 ]
if(nu == round(nu)) ## <==> sin(pi*nu) == 0
return(BesselI(z, -nu, expon.scaled, nSeq=nSeq, verbose=verbose))
## else
nu. <- -nu + seq_len(nSeq) - 1
kf <- rep(2/pi*sin(pi*nu.), each=nz)
if (expon.scaled) kf <- kf * rep(exp(-z-abs(zr)), nSeq)
return( BesselI(z, -nu, expon.scaled, nSeq=nSeq, verbose=verbose) +
kf*BesselK(z, -nu, expon.scaled, nSeq=nSeq, verbose=verbose))
}
## else nu >= 0 :
r <- if(isNum) numeric(nz * nSeq) else complex(nz * nSeq)
if(nSeq > 1) r <- matrix(r, nz, nSeq)
for(i in seq_len(nz)) {
ri <- .C(zbesi,
zr[i], zi[i],
fnu = nu,
kode= as.integer(1L + as.logical(expon.scaled)),
# 1 or 2, exactly as desired
n = as.integer(nSeq),
cyr = double(nSeq),
cyi = double(nSeq),
nz = integer(1),
ierr = as.integer(verbose))
if(ri$ierr) {
f.x <- sprintf("'zbesi(%g %s %gi, nu=%g)'", zr[i],
c("-","+")[1+(zi[i] >= 0)], abs(zi[i]), nu)
if(ri$ierr == 3)
warning(gettextf(
"%s large arguments -> precision loss (of at least half machine accuracy)", f.x),
domain = NA)
else if(ri$ierr == 2) {
if(verbose)
message(gettextf("%s -> overflow ; returning Inf\n", f.x))
ri$cyr <- ri$cyi <- Inf
}
else if(ri$ierr == 4) {
warning(gettextf("%s -> ierr=4: |z| or nu too large\n", f.x),
domain = NA)
## FIXME: In some cases, the answer should just be 'Inf' without any warning
ri$cyr[] <- NaN
ri$cyi[] <- if(isNum) 0 else NaN
}
else stop(gettextf("%s unexpected error 'ierr = %d'", f.x, ri$ierr),
domain = NA)
}
rz <- if(isNum && all(!is.na(ri$cyi) & ri$cyi == 0.)) ri$cyr
else complex(real = ri$cyr,
imaginary = ri$cyi)
if(nSeq > 1) r[i,] <- rz else r[i] <- rz
}
r
} ## I()
## 631: zbesj(zr, zi, fnu, kode, n, cyr, cyi, nz, ierr)
BesselJ <- function(z, nu, expon.scaled = FALSE, nSeq = 1, verbose = 0)
{
nz <- length(z)
if(nz == 0) return(z)
isNum <- is.numeric(z)
if(isNum) {
zr <- as.double(z)
zi <- numeric(nz)
} else if(is.complex(z)) {
zr <- Re(z)
zi <- Im(z)
} else stop("'z' must be complex or numeric")
nu <- as.double(nu)
stopifnot(length(nu) == 1, length(verbose) == 1,
nSeq >= 1, nSeq == as.integer(nSeq))
if(nu < 0) {
## J(-fnu,z) = J(fnu,z)*cos(pi*fnu) - Y(fnu,z)*sin(pi*fnu)
if(expon.scaled)
stop("'expon.scaled=TRUE' not yet implemented for nu < 0")
nu. <- rep(-nu + seq_len(nSeq) - 1, each=nz)
return(BesselJ(z, -nu, nSeq=nSeq, verbose=verbose)*cospi(nu.) -
if(nu == round(nu)) 0 else
BesselY(z, -nu, nSeq=nSeq, verbose=verbose)*sinpi(nu.))
}
## else nu >= 0 :
r <- if(isNum) numeric(nz * nSeq) else complex(nz * nSeq)
if(nSeq > 1) r <- matrix(r, nz, nSeq)
for(i in seq_len(nz)) {
ri <- .C(zbesj,
zr[i], zi[i],
fnu = nu,
kode= as.integer(1L + as.logical(expon.scaled)),
# 1 or 2, exactly as desired
n = as.integer(nSeq),
cyr = double(nSeq),
cyi = double(nSeq),
nz = integer(1),
ierr = as.integer(verbose))
if(ri$ierr) {
f.x <- sprintf("'zbesj(%g %s %gi, nu=%g)'", zr[i],
c("-","+")[1+(zi[i] >= 0)], abs(zi[i]), nu)
if(ri$ierr == 3)
warning(sprintf(
"%s large arguments -> precision loss (of at least half machine accuracy)",
f.x))
else if(ri$ierr == 2) {
if(verbose)
message(sprintf("%s -> overflow ; returning Inf\n", f.x))
ri$cyr <- ri$cyi <- Inf
}
else if(ri$ierr == 4) {
warning(gettextf("%s -> ierr=4: |z| or nu too large\n", f.x),
domain = NA)
## FIXME: In some cases, the answer should just be Inf or 0 (w/o warning!)
ri$cyr[] <- NaN
ri$cyi[] <- if(isNum) 0 else NaN
}
else stop(gettextf("%s unexpected error 'ierr = %d'", f.x, ri$ierr),
domain = NA)
}
rz <- if(isNum) ri$cyr
else complex(real = ri$cyr,
imaginary = ri$cyi)
if(nSeq > 1) r[i,] <- rz else r[i] <- rz
}
r
} ## J()
## 899: zbesk(zr, zi, fnu, kode, n, cyr, cyi, nz, ierr)
BesselK <- function(z, nu, expon.scaled = FALSE, nSeq = 1, verbose = 0)
{
nz <- length(z)
if(nz == 0) return(z)
isNum <- is.numeric(z)
if(isNum) {
zr <- as.double(z)
zi <- numeric(nz)
} else if(is.complex(z)) {
zr <- Re(z)
zi <- Im(z)
} else stop("'z' must be complex or numeric")
nu <- as.double(nu)
stopifnot(length(nu) == 1, length(verbose) == 1,
nSeq >= 1, nSeq == as.integer(nSeq))
if(nu < 0) {
## K(-nu,z) = K(nu,z)
return(BesselK(z, -nu, expon.scaled, nSeq=nSeq, verbose=verbose))
}
## else nu >= 0 :
r <- if(isNum) numeric(nz * nSeq) else complex(nz * nSeq)
if(nSeq > 1) r <- matrix(r, nz, nSeq)
for(i in seq_len(nz)) {
ri <- .C(zbesk,
zr[i], zi[i],
fnu = nu,
kode= as.integer(1L + as.logical(expon.scaled)),
# 1 or 2, exactly as desired
n = as.integer(nSeq),
cyr = double(nSeq),
cyi = double(nSeq),
nz = integer(1),
ierr = as.integer(verbose))
if(ri$ierr) {
f.x <- sprintf("'zbesk(%g %s %gi, nu=%g)'", zr[i],
c("-","+")[1+(zi[i] >= 0)], abs(zi[i]), nu)
if(ri$ierr == 3)
warning(sprintf(
"%s large arguments -> precision loss (of at least half machine accuracy)",
f.x))
else if(ri$ierr == 2) {
if(verbose)
message(sprintf("%s -> overflow ; returning Inf\n", f.x))
ri$cyr <- ri$cyi <- Inf
}
else if(ri$ierr == 4) {
warning(gettextf("%s -> ierr=4: |z| or nu too large\n", f.x),
domain = NA)
## FIXME: In some cases, the answer should just be Inf or 0 (w/o warning!)
ri$cyr[] <- NaN
ri$cyi[] <- if(isNum) 0 else NaN
}
else stop(gettextf("%s unexpected error 'ierr = %d'", f.x, ri$ierr),
domain = NA)
}
rz <- if(isNum && all(!is.na(ri$cyi) & ri$cyi == 0.)) ri$cyr
else complex(real = ri$cyr,
imaginary = ri$cyi)
if(nSeq > 1) r[i,] <- rz else r[i] <- rz
}
r
} ## K()
## 1182: zbesy(zr, zi, fnu, kode, n, cyr, cyi, nz, cwrkr, cwrki, ierr)
BesselY <- function(z, nu, expon.scaled = FALSE, nSeq = 1, verbose = 0)
{
nz <- length(z)
if(nz == 0) return(z)
isNum <- is.numeric(z)
if(isNum) {
zr <- as.double(z)
zi <- numeric(nz)
} else if(is.complex(z)) {
zr <- Re(z)
zi <- Im(z)
} else stop("'z' must be complex or numeric")
nu <- as.double(nu)
stopifnot(length(nu) == 1, length(verbose) == 1,
nSeq >= 1, nSeq == as.integer(nSeq))
if(nu < 0) {
## Y(-fnu,z) = Y(fnu,z)*cos(pi*fnu) + J(fnu,z)*sin(pi*fnu)
if(expon.scaled)
stop("'expon.scaled=TRUE' not yet implemented for nu < 0")
nu. <- rep(-nu + seq_len(nSeq) - 1, each=nz)
return(BesselY(z, -nu, nSeq=nSeq, verbose=verbose)*cospi(nu.) +
if(nu == round(nu)) 0 else
BesselJ(z, -nu, nSeq=nSeq, verbose=verbose)*sinpi(nu.))
}
## else nu >= 0 :
isNum <- isNum && all(zr >= 0)
r <- if(isNum) numeric(nz * nSeq) else complex(nz * nSeq)
if(nSeq > 1) r <- matrix(r, nz, nSeq)
for(i in seq_len(nz)) {
if(zr[i] == 0 && zi[i] == 0) {
rz <- if(isNum) -Inf else
## A limit z -> 0 depends on the *direction*; it is
## ``a version of Inf"', i.e. the only *complex* Inf, if you think
## of the complex sphere. --> we use the same as 1/(0+0i):
1/(0+0i)
} else {
## 1182: zbesy(zr, zi, fnu, kode, n, cyr, cyi, nz, cwrkr, cwrki, ierr)
ri <- .C(zbesy,
zr[i], zi[i],
fnu = nu,
kode= as.integer(1L + as.logical(expon.scaled)),
# 1 or 2, exactly as desired
n = as.integer(nSeq),
cyr = double(nSeq),
cyi = double(nSeq),
nz = integer(1),
cwrkr= double(nSeq),
cwrki= double(nSeq),
ierr = as.integer(verbose))
if(ri$ierr) {
f.x <- sprintf("'zbesy(%g %s %gi, nu=%g)'", zr[i],
c("-","+")[1+(zi[i] >= 0)], abs(zi[i]), nu)
if(ri$ierr == 3)
warning(sprintf(
"%s large arguments -> precision loss (of at least half machine accuracy)",
f.x))
else if(ri$ierr == 2) {
if(verbose)
message(sprintf("%s -> overflow ; returning Inf\n", f.x))
ri$cyr <- ri$cyi <- Inf
}
else if(ri$ierr == 4) {
warning(gettextf("%s -> ierr=4: |z| or nu too large\n", f.x),
domain = NA)
## FIXME: In some cases, the answer should just be Inf or 0 (w/o warning!)
ri$cyr[] <- NaN
ri$cyi[] <- if(isNum) 0 else NaN
}
else stop(gettextf("%s unexpected error 'ierr = %d'", f.x, ri$ierr),
domain = NA)
}
rz <- if(isNum) ri$cyr
else complex(real = ri$cyr,
imaginary = ri$cyi)
}
if(nSeq > 1) r[i,] <- rz else r[i] <- rz
}
r
} ## Y()
##---------------- Hankel function H() ------------------
## 10: zbesh(zr, zi, fnu, kode, m, n, cyr, cyi, nz, ierr)
BesselH <- function(m, z, nu, expon.scaled = FALSE, nSeq = 1, verbose = 0)
{
## c***keywords H-bessel functions,bessel functions of complex argument,
## c bessel functions of third kind,hankel functions
## c***author amos, donald e., sandia national laboratories
## c***purpose to compute the H-bessel functions of a complex argument
## c***description
## c
## c ***a double precision routine***
## c on kode=1, zbesh computes an n member sequence of complex
## c hankel (bessel) functions cy(j) = H(m,fnu+j-1,z) for kinds m=1 or 2,
## c real, nonnegative orders fnu+j-1, j=1,...,n, and complex
## c z != cmplx(0.,0.) in the cut plane -pi < arg(z) <= pi.
## c on kode=2, zbesh returns the scaled Hankel functions
## c
## c cy(i)= exp(-mm*z*i) * H(m,fnu+j-1,z) mm=3-2*m, i**2=-1.
## c
## c which removes the exponential behavior in both the upper and
## c lower half planes. definitions and notation are found in the
## c nbs handbook of mathematical functions (ref. 1).
m <- as.integer(m)
stopifnot(length(m) == 1, m == 1 || m == 2)
nz <- length(z)
if(nz == 0) return(z)
isNum <- is.numeric(z)
if(isNum) {
zr <- as.double(z)
zi <- numeric(nz)
} else if(is.complex(z)) {
zr <- Re(z)
zi <- Im(z)
} else stop("'z' must be complex or numeric")
nu <- as.double(nu)
stopifnot(length(nu) == 1, length(verbose) == 1,
nSeq >= 1, nSeq == as.integer(nSeq))
if(nu < 0) {
## H(1,-fnu,z) = H(1,fnu,z)*cexp( pi*fnu*i)
## H(2,-fnu,z) = H(2,fnu,z)*cexp(-pi*fnu*i) ; i^2=-1
if(expon.scaled)
stop("'expon.scaled=TRUE' not yet implemented for nu < 0")
## exp(i * pi*nu) = cos(pi*nu) + i* sin(pi*nu) -- should use cospi() / sinpi() ?
pnu <- rep(c(1i,-1i)[m] * pi*(-nu + seq_len(nSeq) - 1), each=nz)
return(BesselH(m, z, -nu, nSeq=nSeq, verbose=verbose)* exp(pnu))
}
## else nu >= 0 :
r <- if(isNum) numeric(nz * nSeq) else complex(nz * nSeq)
if(nSeq > 1) r <- matrix(r, nz, nSeq)
for(i in seq_len(nz)) {
## zbesh(zr, zi, fnu, kode, m, n, cyr, cyi, nz, ierr)
ri <- .C(zbesh,
zr[i], zi[i],
fnu = nu,
kode= as.integer(1L + as.logical(expon.scaled)),
# 1 or 2, exactly as desired
m = m,
n = as.integer(nSeq),
cyr = double(nSeq),
cyi = double(nSeq),
nz = integer(1),
ierr = as.integer(verbose))
if(ri$ierr) {
f.x <- sprintf("'zbesh(%g %s %gi, nu=%g)'", zr[i],
c("-","+")[1+(zi[i] >= 0)], abs(zi[i]), nu)
if(ri$ierr == 3)
warning(sprintf(
"%s large arguments -> precision loss (of at least half machine accuracy)",
f.x))
else if(ri$ierr == 2) {
if(verbose)
message(sprintf("%s -> overflow ; returning Inf\n", f.x))
ri$cyr <- ri$cyi <- Inf
}
else if(ri$ierr == 4) {
warning(gettextf("%s -> ierr=4: |z| or nu too large\n", f.x),
domain = NA)
## FIXME: In some cases, the answer should just be Inf or 0 (w/o warning!)
ri$cyr[] <- NaN
ri$cyi[] <- if(isNum) 0 else NaN
}
else stop(gettextf("%s unexpected error 'ierr = %d'", f.x, ri$ierr),
domain = NA)
}
rz <- if(isNum && all(!is.na(ri$cyi) & ri$cyi == 0.)) ri$cyr
else complex(real = ri$cyr,
imaginary = ri$cyi)
if(nSeq > 1) r[i,] <- rz else r[i] <- rz
}
r
} ## H()
##---------------- Airy Functions Ai() Bi() ------------------
## 1471: zairy(zr, zi, id, kode, air, aii, nz, ierr)
## 1867: zbiry(zr, zi, id, kode, bir, bii, ierr)
AiryA <- function(z, deriv = 0, expon.scaled = FALSE, verbose = 0)
{
## Purpose: to compute airy functions ai(z) and dai(z) for complex z
## airy function : "Bessel functions of order one third"
deriv <- as.integer(deriv)
stopifnot(length(deriv) == 1, deriv == 0 || deriv == 1, length(verbose) == 1)
nz <- length(z)
if(nz == 0) return(z)
isNum <- is.numeric(z)
if(isNum) {
zr <- as.double(z)
zi <- numeric(nz)
} else if(is.complex(z)) {
zr <- Re(z)
zi <- Im(z)
} else stop("'z' must be complex or numeric")
r <- if(isNum) numeric(nz) else complex(nz)
for(i in seq_len(nz)) {
## zairy(zr, zi, id, kode, air, aii, nz, ierr)
ri <- .C(zairy,
zr[i], zi[i],
id = deriv,
kode= as.integer(1L + as.logical(expon.scaled)),
# 1 or 2, exactly as desired
air = double(1),
aii = double(1),
nz = integer(1),
ierr = as.integer(verbose))
if(ri$ierr) {
f.x <- sprintf("'zairy(%g %s %gi, deriv=%d)'", zr[i],
c("-","+")[1+(zi[i] >= 0)], abs(zi[i]), deriv)
if(ri$ierr == 3)
warning(sprintf(
"%s large arguments -> precision loss (of at least half machine accuracy)",
f.x))
else if(ri$ierr == 2) {
if(verbose)
message(sprintf("%s -> overflow ; returning Inf\n", f.x))
ri$air <- ri$aii <- Inf
}
else if(ri$ierr == 4) {
warning(gettextf("%s -> ierr=4: |z| too large\n", f.x), domain = NA)
## FIXME: In some cases, the answer should just be Inf or 0 or .... (w/o warning!)
ri$air <- NaN
ri$aii <- if(isNum) 0 else NaN
}
else stop(gettextf("%s unexpected error 'ierr = %d'", f.x, ri$ierr),
domain = NA)
}
r[i] <- if(isNum) ri$air else complex(real = ri$air, imaginary = ri$aii)
}
r
} ## AiryA()
AiryB <- function(z, deriv = 0, expon.scaled = FALSE, verbose = 0)
{
## Purpose: to compute airy functions bi(z) and dbi(z) for complex z
## airy function : "Bessel functions of order one third"
deriv <- as.integer(deriv)
stopifnot(length(deriv) == 1, deriv == 0 || deriv == 1, length(verbose) == 1)
nz <- length(z)
if(nz == 0) return(z)
isNum <- is.numeric(z)
if(isNum) {
zr <- as.double(z)
zi <- numeric(nz)
} else if(is.complex(z)) {
zr <- Re(z)
zi <- Im(z)
} else stop("'z' must be complex or numeric")
r <- if(isNum) numeric(nz) else complex(nz)
for(i in seq_len(nz)) {
## zairy(zr, zi, id, kode, air, aii, nz, ierr)
ri <- .C(zbiry,
zr[i], zi[i],
id = deriv,
kode= as.integer(1L + as.logical(expon.scaled)),
# 1 or 2, exactly as desired
bir = double(1),
bii = double(1),
ierr = as.integer(verbose))
if(ri$ierr) {
f.x <- sprintf("'zairy(%g %s %gi, deriv=%d)'", zr[i],
c("-","+")[1+(zi[i] >= 0)], abs(zi[i]), deriv)
if(ri$ierr == 3)
warning(sprintf(
"%s large arguments -> precision loss (of at least half machine accuracy)",
f.x))
else if(ri$ierr == 2) {
if(verbose)
message(sprintf("%s -> overflow ; returning Inf\n", f.x))
ri$bir <- ri$bii <- Inf
}
else if(ri$ierr == 4) {
warning(gettextf("%s -> ierr=4: |z| too large\n", f.x), domain = NA)
## FIXME: In some cases, the answer should just be Inf or 0 or .... (w/o warning!)
ri$bir <- NaN
ri$bii <- if(isNum) 0 else NaN
}
else stop(gettextf("%s unexpected error 'ierr = %d'", f.x, ri$ierr),
domain = NA)
}
r[i] <- if(isNum) ri$bir else complex(real = ri$bir, imaginary = ri$bii)
}
r
} ## AiryB()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.