1ellipsePoints <- function(a,b, alpha = 0, loc = c(0,0), n = 201,
2                          keep.ab.order = FALSE)
3{
4    ## Purpose: ellipse points,radially equispaced, given geometric par.s
5    ## -------------------------------------------------------------------------
6    ## Arguments: a, b : length of half axes in (x,y) direction
7    ##            alpha: angle (in degrees) for rotation
8    ##            loc  : center of ellipse
9    ##            n    : number of points
10    ## -------------------------------------------------------------------------
11    ## Author: Martin Maechler, Date: 19 Mar 2002
12
13    stopifnot(is.numeric(a), is.numeric(b))
14    reorder <- a < b && keep.ab.order
15    B <- min(a,b)
16    A <- max(a,b)
17    ## B <= A
18    d2 <- (A-B)*(A+B) ## = A^2 - B^2
19    phi <- 2*pi*seq(0,1, len = n)
20    sp <- sin(phi)
21    cp <- cos(phi)
22    r <- a*b / sqrt(B^2 + d2 * sp^2)
23    xy <- r * if(reorder) cbind(sp, cp) else cbind(cp, sp)
24    ## xy are the ellipse points for alpha = 0 and loc = (0,0)
25    al <- alpha * pi/180
26    ca <- cos(al)
27    sa <- sin(al)
28    xy %*% rbind(c(ca, sa), c(-sa, ca)) + cbind(rep(loc[1],n),
29                                                rep(loc[2],n))
30}
31