R + plotly: solid of revolution

This doesn't answer your question, but it will give a result you can interact with in a web page: don't use plot_ly, use rgl. For example,

library(rgl)

# Your initial values...

r <- function(x) x^2
int <- c(1, 3)
nx <- 20
ntheta <- 36

# Set up x and colours for each x

x <- seq(int[1], int[2], length.out = nx)
cols <- colorRampPalette(c("blue", "yellow"), space = "Lab")(nx)

clear3d()
shade3d(turn3d(x, r(x), n = ntheta,  smooth = TRUE, 
        material = list(color = rep(cols, each = 4*ntheta))))
aspect3d(1,1,1)
decorate3d()
rglwidget()

You could do better on the colours with some fiddling: you probably want to create a function that uses x or r(x) to set the colour instead of just repeating the colours the way I did.

Here's the result:

interesting question, I've struggled to use the surface density to improve on your solution. There is a hack you could do with layering multiple lines, that comes out nice for this e.g. Only changes made to the original eg is to use lots more x points: nx to 1000, and change add_markers to add_lines. Might not be scalable, but works fine for this size of data :)

library(dplyr)
library(plotly)

# radius depends on x
r <- function(x) x^2

# interval of interest
int <- c(1, 3)

# number of points along the x-axis
nx <- 1000

# number of points along the rotation
ntheta <- 36

# set x points and get corresponding radii
coords <- data_frame(x = seq(int[1], int[2], length.out = nx), r = r(x))

# for each x: rotate r to get y and z coordinates
# edit: ensure 0 and pi are both amongst the angles used
coords %<>%
  rowwise() %>%
  do(data_frame(x = .$x, r = .$r,
                theta = seq(0, pi, length.out = ntheta / 2 + 1) %>%
                  c(pi + .[-c(1, length(.))]))) %>%

  ungroup %>%
  mutate(y = r * cos(theta), z = r * sin(theta))

# plot points to make sure the coordinates define the desired shape
coords %>%
  plot_ly(x = ~x, y = ~y, z = ~z, color = ~x) %>%
  add_lines()

Best, Jonny

I have had another crack at it and have a closer solution, using the "surface" type. What helped was looking at the results of your first surface plot with nx = 5 and ntheta = 18. The reason it's jaggardy is because of the way its linking up the columns in zs (across the x points). It's having to link from part way up the larger ring around it, and this causes the density to spike up to meet this point.

I can't get rid of this jaggardy behaviour 100%. I've made these changes:

1. add some small points to theta around the edges: where the two densities are joined. This reduces the size of the jaggardy part as there are some more points close to the boundary
2. calculation to mod zs to zs2: ensure that each ring has an equal dimension to the ring outside, by adding the 0's in.
3. increased nx to 40 and reduced ntheta to 18 - more x's makes step smaller. reduce ntheta for run time, as I've added on more points

the steps come in how it tries to join up the x rings. In theory if you have more x rings it should remove this jaggardiness, but that's time consuming to run.

I don't think this answers the Q 100%, and I'm unsure if this library is the best for this job. Get in touch if any Q's.

library(dplyr)
library(plotly)

# radius depends on x
r <- function(x) x^2

# interval of interest
int <- c(1, 3)

# number of points along the x-axis
nx <- 40

# number of points along the rotation
ntheta <- 18

# set x points and get corresponding radii
coords <- data_frame(x = seq(int[1], int[2], length.out = nx), r = r(x))

# theta: add small increments at the extremities for the density plot
theta <- seq(0, pi, length.out = ntheta / 2 + 1)
theta <- c(theta, pi + theta)
theta <- theta[theta != 2*pi]
inc <- 0.00001
theta <- c(theta, inc, pi + inc, pi - inc, 2*pi - inc)
theta <- sort(theta)

coords %<>%
  rowwise() %>%
  do(data_frame(x = .$x, r = .$r, theta = theta)) %>%
  ungroup %>%
  mutate(y = r * cos(theta), z = r * sin(theta))

# get all x & y values used (sort to connect halves on the side)
xs <-
  unique(coords$x) %>%
  sort
ys <-
  unique(coords$y) %>%
  sort

# for each possible x/y pair: get z^2 value
coords <-
  expand.grid(x = xs, y = ys) %>%
  as_data_frame %>%
  mutate(r = r(x), z2 = r^2 - y^2)

# format z coordinates above x/y plane as matrix where columns
# represent x and rows y
zs <- matrix(sqrt(coords$z2), ncol = length(xs), byrow = TRUE)
zs2 <- zs

L <- ncol(zs)
for(i in (L-1):1){
  w <- which(!is.na(zs[, (i+1)]) & is.na(zs[, i]))
  zs2[w, i] <- 0
}

# format x coordiantes as matrix as above (for color gradient)
gradient <-
  rep(xs, length(ys)) %>%
  matrix(ncol = length(xs), byrow = TRUE)

# plot upper half of shape as surface
p <- plot_ly(x = xs, y = ys, z = zs2, surfacecolor = gradient,
             type = "surface", colorbar = list(title = 'x'))

# plot lower have of shape as second surface
p %>%
  add_surface(z = -zs2, showscale = FALSE)