##
## e l l i p . R Elliptic Integrals
##
ellipke <- function(m, tol = .Machine$double.eps) {
stopifnot(is.numeric(m))
m <- c(m)
if (any(m < 0) || any(m > 1))
stop("Some elements of argument 'm' are out of range.")
a0 <- 1
b0 <- sqrt(1-m)
s0 <- m
i1 <- 0
mm <- 1
while (mm > tol) {
a1 <- (a0+b0)/2
b1 <- sqrt(a0*b0)
c1 <- (a0-b0)/2
i1 <- i1 + 1
w1 <- 2^i1 * c1^2
mm <- max(w1)
s0 <- s0 + w1
a0 <- a1
b0 <- b1
}
k <- pi / (2*a1)
e <- k * (1-s0/2)
im <- finds(m == 1)
if (!isempty(im)) {
e[im] <- ones(length(im), 1)
k[im] <- Inf
}
return(list(k = k, e = e))
}
## Jacobi elliptic functions
ellipj <- function(u, m, tol = .Machine$double.eps) {
stopifnot(is.numeric(u), is.numeric(m) || is.complex(m))
u <- c(u); m <- c(m)
if (length(u) == 1) {
u <- rep(u, length(m))
} else if (length(m) == 1) {
m <- rep(m, length(u))
} else {
if (length(u) != length(m))
stop("Arguments 'u' and 'm' must be of the same length.")
}
if (any(m < 0) || any(m > 1))
stop("Some elements of argument 'm' are out of range.")
mmax <- length(u); chunk <- 10
cn <- sn <- dn <- numeric(mmax)
a <- b <- cc <- matrix(0, nrow = chunk, ncol = mmax)
a[1, ] <- 1
b[1, ] <- sqrt(1-m)
cc[1, ] <- sqrt(m)
n <- numeric(mmax)
i <- 1
while (any(abs(cc[i, ]) > tol)) {
i <- i + 1
if (i > nrow(a)) {
a <- rbind(a, matrix(0, chunk, mmax))
b <- rbind(b, matrix(0, chunk, mmax))
cc <- rbind(cc, matrix(0, chunk, mmax))
}
a[i, ] <- 0.5 * (a[i-1, ] + b[i-1, ])
b[i, ] <- sqrt(a[i-1, ] * b[i-1, ])
cc[i, ] <- 0.5 * (a[i-1, ] - b[i-1, ])
inds <- which(abs(cc[i, ]) <= tol & abs(cc[i-1, ]) > tol)
if (!isempty(inds)) {
mi <- 1; ni <- length(inds) # [mi,ni] <- size(inds)
n[inds] <- rep(i-1, ni) # repmat((i-1), mi, ni)
}
}
phin <- matrix(0, nrow = i, ncol = mmax)
phin[i, ] <- 2^n * a[i, ] * u
while (i > 1) {
i <- i - 1
inds <- which(n >= i)
phin[i, ] <- phin[i+1, ]
if (!isempty(inds)) {
phin[i, inds] <- 0.5 * (asin(cc[i+1, inds] *
sin(phin[i+1, inds] %% (2*pi)) / a[i+1, inds]) + phin[i+1,inds])
}
}
# the general case
sn <- sin(phin[1, ] %% (2*pi))
cn <- cos(phin[1, ] %% (2*pi))
dn <- sqrt(1 - m * sn^2)
# some special cases
m1 <- which(m == 1) # special case m = 1
sn[m1] <- tanh(u[m1])
cn[m1] <- sech(u[m1])
dn[m1] <- sech(u[m1])
dn[m == 0] <- 1 # special case m = 0
return(list(sn = sn, cn = cn, dn = dn))
}