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"))

# Source 3 - 4 ends


# side-view cross-section plots Using base R graphics to plot the matrix plots
# ('side-view cross-section plots')

# 'side view from the 'left'' add r = 15 to r = 40 for h_grid, W_grid
matplot(h_grid, W_grid, type = "l", xlab = "h (m)", ylab = "W (N)")

# 'side view from the 'right'' add h = 25 to h = 50 for r_cylinder_grid, W_grid
matplot(t(r_cylinder_grid), t(W_grid), type = "l", xlab = "r_{cylinder} (mm)", ylab = "W (N)")



Solution Using GNU Octave


rho = 1000;
g = 9.81;
r_cylinder = 15:5:40;
h = 25:5:50;

[r_cylinder_grid, h_grid] = meshgrid(r_cylinder, h);

W_grid = rho*g*(pi*(r_cylinder_grid/1000).^2 .* h_grid);

mesh(r_cylinder_grid, h_grid, W_grid);
xlabel('r_{cylinder}, mm');
ylabel('h, m');
zlabel('W, N');

print("example1.png")




rho = 1000;
g = 9.81;
r_cylinder = 15:5:40;
h = 25:5:50;

[r_cylinder_grid, h_grid] = meshgrid(r_cylinder, h);

W_grid = rho*g*(pi*(r_cylinder_grid/1000).^2 .* h_grid);

meshz(r_cylinder_grid, h_grid, W_grid);
xlabel('r_{cylinder}, mm');
ylabel('h, m');
zlabel('W, N');

print("example2.png")




rho = 1000;
g = 9.81;
r_cylinder = 15:5:40;
h = 25:5:50;

[r_cylinder_grid, h_grid] = meshgrid(r_cylinder, h);

W_grid = rho*g*(pi*(r_cylinder_grid/1000).^2 .* h_grid);

surf(r_cylinder_grid, h_grid, W_grid);
xlabel('r_{cylinder}, mm');
ylabel('h, m');
zlabel('W, N');

print("example3.png")




rho = 1000;
g = 9.81;
r_cylinder = 15:5:40;
h = 25:5:50;

[r_cylinder_grid, h_grid] = meshgrid(r_cylinder, h);

W_grid = rho*g*(pi*(r_cylinder_grid/1000).^2 .* h_grid);

[C, hd1] = contour(r_cylinder_grid, h_grid, W_grid);
clabel(C, hd1);
xlabel('r_{cylinder}, mm');
ylabel('h, m');

% save the object C in the "MATLAB v4 binary data format"
% save("-v4", "C.mat", "C")

print("example4.png")




rho = 1000;
g = 9.81;
r_cylinder = 15:5:40;
h = 25:5:50;

[r_cylinder_grid, h_grid] = meshgrid(r_cylinder, h);

W_grid = rho*g*(pi*(r_cylinder_grid/1000).^2 .* h_grid);

% add r = 15 to r = 40 for h_grid, W_grid
plot(h_grid, W_grid);
xlabel('h, m');
ylabel('W, N');

% add h = 25 to h = 50 for r_cylinder_grid, W_grid
plot(r_cylinder_grid', W_grid')
xlabel('r_{cylinder}, mm');
ylabel('W, N');

print("example5.png")



Sources used in the R code

Source 1
Statistical tools for high-throughput data analysis (STHDA): Impressive package for 3D and 4D graph - R software and data visualization. See http://www.sthda.com/english/wiki/impressive-package-for-3d-and-4d-graph-r-software-and-data-visualization.

Source 2
Fifty ways to draw a volcano using package plot3D by Karline Soetaert. See https://cloud.r-project.org/web/packages/plot3D/vignettes/volcano.pdf.

Source 3
ggplot2 - R: How to Label Specific Contours using direct.label - Stack Overflow answered by a different ben on Jul 21 2014. See https://stackoverflow.com/questions/10675387/r-how-to-label-specific-contours-using-direct-label.

Source 4
R adding legend and directlabels to ggplot2 contour plot - Stack Overflow answered and edited by Sandy Muspratt on Jul 2 2016. See https://stackoverflow.com/questions/38154679/r-adding-legend-and-directlabels-to-ggplot2-contour-plot.

Source 5
r - Superscript and subscript axis labels in ggplot2 - Stack Overflow answered and edited by akrun on Dec 12 2014. See https://stackoverflow.com//questions/27445452/superscript-and-subscript-axis-labels-in-ggplot2.

Source 6
r - RMarkdown Octave Plot - Stack Overflow answered by danlooo on Mar 17, 2022. See https://stackoverflow.com/questions/71512367/rmarkdown-octave-plot



Works Cited

Jay B. Brockman, Introduction to Engineering: Modeling and Problem Solving, Hoboken, New Jersey: John Wiley & Sons, Inc., 2009, pages 470-473.



Donations accepted with Liberapay

If you would like to contribute to the continued development of Irucka Embry’s R packages and/or Irucka Embry’s R Examples, please feel free to donate via the link below:

Please feel free to review Irucka Embry (iaembry)’s profile on Liberapay.