Replicate the R code

Install the packages

Note: If you wish to replicate the R code below, then you will need to copy and paste the following command in R first to make sure you have all the packages and their dependencies installed:

install.packages("install.load", "import", "plot3D", "ggplot2", "rmatio", "reshape2", "directlabels", "gsubfn")
# install the install.load and import packages

## ContourFunctions # contour plot

## use for contour plots -- https://cran.r-project.org/web/packages/ContourFunctions/vignettes/Introduction_to_the_cf_R_package.html



Load the packages

Note: If you wish to replicate the R code below, then you will need to copy and paste the following commands in R first to make sure you have all the packages loaded:

install.load::load_package("plot3D", "ggplot2", "rmatio", "reshape2", "directlabels",
    "gsubfn")
# load needed packages using the load_package function from the install.load
# package (it is assumed that you have already installed these packages)


import::from(pracma, meshgrid)
# import meshgrid from the pracma package



Introduction

The following example, "Force on the Piston of a Pump versus Well Depth and Cylinder Radius", is solved using R (incomplete graphical solution) and GNU Octave. (Brockman 470-473)

The following R code is incomplete as I am still searching for the best way to obtain the same plots in R as with GNU Octave/MATLAB. I posted this question - Q: r - mesh, meshz, surf, and contour plots similar to GNU Octave/MATLAB - online at Stack Overflow, but I have not received an answer yet. Please feel free to post an answer or leave comments. Thank you.



Definition of variables

"W = load force"

"r_cylinder = radius of the cylinder, which varies from 15 - 40 mm"

"h = depth of well, which varies from 25 - 50 m"

"rho, density of water, = 1000 kg/m^3"

"W = W(rcylinder, h) = rho * g * (pi * r_cylinder^2 * h)"

Source: Brockman 470 - 473



Partial Solution Using R

# R plotmath grDevices - subscripts

# GNU Octave xlabel ylabel subscripts _ underscore

rho <- 1000

g <- 9.81

r_cylinder <- seq(15, 40, by = 5)

h <- seq(25, 50, by = 5)


# mesh plot 'for piston force example data'

list[r_cylinder_grid, h_grid] <- meshgrid(r_cylinder, h)

W_grid <- rho * g * (pi * (r_cylinder_grid/1000)^2 * h_grid)

persp3D(r_cylinder_grid, h_grid, W_grid, xlab = "r_{cylinder}, mm", ylab = "h, m",
    zlab = "W, N", theta = 100, phi = -30, expand = 0.5, ticktype = "detailed", ltheta = 120,
    bty = "g", facets = FALSE, col = NULL)  # Source 1 - 2

# plot_ly(x = r_cylinder, y = h_grid, z = W_grid, type = 'mesh3d')




# meshz plot 'adds 'curtain' around plot to highlight z-values at edges'

# Not available



# surface plot

surf3D(r_cylinder_grid, h_grid, W_grid, xlab = "r_{cylinder}, mm", ylab = "h, m",
    zlab = "W, N", theta = 100, phi = -30, expand = 0.5, ticktype = "detailed", ltheta = 120,
    bty = "g", col = gg.col(50))  # Source 1

# contour plot

# values obtained from GNU Octave
breaks <- read.mat("C.mat")

breaks <- round(breaks$C, digits = 0)

x <- melt(r_cylinder_grid)

x <- x$value

y <- melt(h_grid)

y <- y$value

z <- melt(W_grid)

z <- z$value


df <- data.frame(x, y, z)

# Source 3 - 4 begins
p <- ggplot(df, aes(x, y, z = z))

p <- p + stat_contour(aes(colour = ..level..), breaks = breaks) + scale_colour_gradient(low = "green",
    high = "purple") + theme_bw() + labs(x = bquote(~r[cylinder] ~ ", mm"), y = "h, m",
    title = "Contour plot")
direct.label(p, list("far.from.others.borders", "calc.boxes", "enlarge.box", box.color = NA,
    fill = "transparent", "draw.rects"))