Plot rotational surface

Problem

Given the transect profile, plot the surface, that is formed by rotating the transect profile around the z-axis.

Cylindrical coordinates

Cylindrical coordinates are described here link.

The three coordinates (\(\rho\), \(\phi\), \(z\)) of a point P are defined as:

• The axial distance or radial distance \(\rho\) is the Euclidean distance from the z-axis to the point P.
• The azimuth \(\phi\) is the angle between the reference direction on the chosen plane and the line from the origin to the projection of P on the plane.
• The axial coordinate or height \(z\) is the signed distance from the chosen plane to the point P.

Transformation of Cartesian coordinates to cylindrical coordinates

cyl2cart <- function (rho, phi, z){
  x <- rho * cos(phi)
  y <- rho * sin(phi)
  z <- z
  return(list(x = x, y = y, z = z))
}

We will use package plot3d

library(plot3D)

Example

Define the transect (\(x\), \(y\))

x <- 0:3
y <- c(5,2,2.5,2)
plot(x,y, type="b",pch=16)

Number of points

n <- length(x)*10
xs <- seq(min(x),max(x),length=n)
ys <- approx(x,y,xout= xs)$y
plot(xs,ys)
points(x,y,col="red",pch=16)

We need to generate a polar mesh and then attribute the measured profile to \(z\).

nc <- 60
mc <- length(xs)
phi <- seq(0, 2* pi, length=nc)
rho <- xs
M <- mesh(rho,phi)
names(M) <- c("rho", "phi")
z <- matrix(rep(ys,nc),n,nc)

Convert to Cartesian coordinates

cart <- cyl2cart(M$rho, M$phi, z)

Plot the surface

par(mar = c(0, 0, 0, 0))
surf3D(cart$x, cart$y, cart$z, border = "black",
       colkey = FALSE, bty = "f",
       phi = 20, theta = 30,scale=FALSE)

par(mar = c(0, 0, 0, 0))
surf3D(cart$x, cart$y, cart$z, border = NA,
       colkey = FALSE, bty = "f",
       phi = 20, theta = 30,scale=FALSE)

par(mar = c(0, 0, 0, 0))
surf3D(cart$x, cart$y, cart$z, border = NA,
       colkey = FALSE, bty = "f",
       phi = 30, theta = 30,scale=FALSE)

All together now …

#' Plot rotated surface with a given profile.
#'
#' @param x numeric vector Profile x coordinates.
#' @param y numeric vector Profile y coordinates.
#' @param nphi numeric Number of angles for rotation.
#' @param ps numeric Start angle.
#' @param pe numeric End angle.
#' @param nout numeric Factor of approximated points for each node in profile.
#' @param ... Additional parameters, see functions \code{\link{persp3D}} and \code{\link{persp}}.
#' @return 
#' @export
#' @seealso \code{\link{persp3D}}, \code{\link{persp}}
#' @note Requires package \code{plot3D}. 
#' @references
#' @keywords function, graphics
#' @title
#' @author Andrej Blejec \email{andrej.blejec@nib.si}
#' @examples
#' # x coordinates
#' x <- 0:3
#' # profile (y coordiantes)
#' y <- c(5,2,2.5,2)
#' plot( x, y, type="b",pch=16)
#' # Plot rotated surface
#' rotated3D(x, y, nout=10,  border = NA,
#'        colkey = FALSE, bty = "f",
#'        phi = 20, theta = 30,scale=FALSE)
#' # Half hull
#' rotated3D(x, y, ps = 0, pe = pi, nout=10,  border = NA,
#'        colkey = FALSE, bty = "f",
#'        phi = 20, theta = 30,scale=FALSE)
#'

rotated3D <- function(x,y, nphi = 60, ps=0, pe = 2*pi, nout = 1,  ...) {
require(plot3D)  
# Function to transform cylindrical coordinates to Cartesian coordinates
cyl2cart <- function (rho, phi, z){
  x <- rho * cos(phi)
  y <- rho * sin(phi)
  z <- z
  return(list(x = x, y = y, z = z))
}

# Approximate data
n <- length(x)*nout
xs <- seq(min(x),max(x),length=n)
ys <- approx(x,y,xout= xs)$y
# Generate cylindrical coordinates
mc <- length(xs)
phi <- seq(ps, pe, length=nphi)
rho <- xs
M <- plot3D::mesh(rho,phi)
names(M) <- c("rho", "phi")
z <- matrix(rep(ys,nphi),n,nphi)


# Convert to Cartesian coordinates

cart <- cyl2cart(M$rho, M$phi, z)
# Plot
surf3D(cart$x, cart$y, cart$z, ...)
}

Define profile and plot surface

# x coordinates
x <- 0:3
# profile (y coordiantes)
y <- c(5,2,2.5,2)
plot( x, y, type="b",pch=16)

Call function rotated3D

rotated3D(x, y,  nout=10,  border = NA,
       colkey = FALSE, bty = "f",
       phi = 20, theta = 30, scale=FALSE)

rotated3D(x, y, ps=0, pe=pi, nout=10,  border = NA,
       colkey = FALSE, bty = "f",
       phi = 20, theta = 30,scale=FALSE)