[R] rgl: cylinder3d() with elliptical cross-section

[Ben Bolker]

Michael Friendly <friendly <at> yorku.ca> writes:

> 
> For a paper dealing with generalized ellipsoids, I want to illustrate in 
> 3D an ellipsoid that is unbounded
> in one dimension, having the shape of an infinite cylinder along, say, 
> z, but whose cross-section in (x,y)
> is an ellipse, say, given by the 2x2 matrix cov(x,y).
> 
> I've looked at rgl:::cylinder3d, but don't see any way to make it 
> accomplish this.  Does anyone have
> any ideas?
> 

  I haven't quite got this in the form you want, but maybe it's
a start (I started from the ellipsoid3d() function in demo("shapes3d") ...

cyl3d <- function(rx=1,ry=1,n=30,ctr=c(0,0,0),
                        zmax=1,
                        qmesh=FALSE,
                        trans = par3d("userMatrix"),...) {
   if (missing(trans) && !rgl.cur()) trans <- diag(4)
   degvec <- seq(0,2*pi,length=n)
   zvec <- seq(0,zmax,length=n)
   ecoord2 <- function(p) {
     c(rx*cos(p[1]),ry*sin(p[1]),p[2])
   }
   v <- apply(expand.grid(degvec,zvec),1,ecoord2)
   if (qmesh) v <- rbind(v,rep(1,ncol(v))) ## homogeneous
   e <- expand.grid(1:(n-1),1:n)
   i1 <- apply(e,1,function(z)z[1]+n*(z[2]-1))
   i2 <- i1+1
   i3 <- (i1+n-1) %% n^2 + 1
   i4 <- (i2+n-1) %% n^2 + 1
   i <- rbind(i1,i2,i4,i3)
   if (!qmesh)
     quads3d(v[1,i],v[2,i],v[3,i],...)
   else return(rotate3d(qmesh3d(v,i,material=...),matrix=trans))
 }

cyl3d(rx=1,ry=2,zmax=5,col="blue")
cyl3d(rx=3,ry=2,zmax=5,col="red",alpha=0.3,ctr=c(2,1,0),
      trans=rotationMatrix(pi/6,0,0,1))