Besides the formula approach, you can do an arbitrarily-fine sum using arrays to get a close-enough approximation. You will have to fine-tune the grid and complete two functions: how to determine whether a point on a grid is within the cuboid, and how to approximate the # of points on the grid that would be present in the whole sphere. Both should not be too difficult.
sphere <- list(x = 0.5, y = 0.4, z = 0.1, radius = 0.8)
# The equations of the six planes
cuboid <- as.character(outer(c("x","y","z"), c(0,1), paste, sep = ' = '))
# "x = 0" "y = 0" "z = 0" "x = 1" "y = 1" "z = 1"
in_sphere <- function(x, y, z)
(x - sphere$x)^2 + (y - sphere$y)^2 + (z - sphere$z)^2 <= sphere$radius^2
in_cuboid <- function(x, y, z) {
# TODO: This seems tricky, fill it in.
# It's easy to do for a box, but you may have to do some linear algebra to
# to do it for a cuboid. For example, an easy computational way to tell if
# a point lies "inside" a cuboid is to iterate over each side, and take
# the line between the side and the point that gives the closest distance of the
# point to the side. Then take the intersection of this line with the remaining
# sides and you should get only one other assuming the cuboid is convex. Then,
# if the distance between the two points on the two sides is less than the distances
# between the tested point and either side, the point lies outside the cuboid.
# Otherwise, it is inside the cuboid.
}
# You need to pre-determine bounds on your grid to search.
# This may be an interesting sub-problem. Note it is hard in general because
# two sides of a cuboid can intersect arbitrarily far in 3-space with
# arbitrarily small modifications of their coefficients.
dims <- c(10, 10, 10) # The number of points in each dimensions to check
box <- c(2, 2, 2) # The actual x-y-z dimensions of the box.
grid <- array(dim = dims, logical(Reduce(`*`, dims)))
iterate <- function(expr, over_points = TRUE) { # Iterate an expression over the grid.
force(over_points)
eval.parent(substitute({
.seqs <- lapply(seq_along(dims), function(i) seq(0, to = box[i], length.out = dims[i]))
for(.x in seq_len(dims[1]))
for(.y in seq_len(dims[2]))
for(.z in seq_len(dims[3])) {
if (over_points) { x <- .seqs[[1]][.x]; y <- .seqs[[2]][.y]; z <- .seqs[[3]][.z] }
else { x <- .x; y <- .y; z <- .z }
expr
}
}))
}
# Mark which points lie in the intersection
iterate(grid[.x, .y, .z] <- in_sphere(x, y, z) && in_cuboid(x, y, z))
# Count the points lying on each surface component.
cuboid_points <- 0; sphere_points <- 0
threshold <- 2*box / min(dims) # If the distance from a point is less than this
# we count it as lying "on" a surface component.
# This is not perfect but should work for a large
# enough grid -- basic calculus.
# Replace the cuboid side equations with the epsilon equations.
cuboid_epsilons <- sapply(cuboid, gsub, pattern = "^(.*)=(.*)$", replacement = paste0("abs(\\1 - \\2) <= ", threshold))
# > cuboid_epsilons
# x = 0 y = 0 z = 0 x = 1 y = 1 z = 1
# "abs(x - 0) <= 0.1" "abs(y - 0) <= 0.1" "abs(z - 0) <= 0.1" "abs(x - 1) <= 0.1" "abs(y - 1) <= 0.1" "abs(z - 1) <= 0.1"
# Whether or not a point lies close to any of the sides.
cuboid_epsilon_condition <- parse(text = paste(cuboid_epsilons, collapse = ' || '))
# expression(abs(x - 0) <= 0.1 || abs(y - 0) <= 0.1 || abs(z - 0) <= 0.1 || abs(x - 1) <= 0.1 || abs(y - 1) <= 0.1 || abs(z - 1) <= 0.1)
# Whether or not a point lies close to the surface of the sphere.
sphere_epsilon_condition <- parse(text =
paste0('abs(', paste(sapply(c('x','y','z'), function(i) paste0('(', i, ' - ', sphere[[i]], ')^2')), collapse = ' + '),
' - ', sphere$radius, '^2) < ', threshold))
# expression(abs((x - 0.5)^2 + (y - 0.4)^2 + (z - 0.1)^2 - 0.8) < 0.1^2)
iterate({
if (grid[.x, .y, .z]) { # In the intersection of cuboid and sphere
if (eval(cuboid_epsilon_condition)) # Lies close to any of the sides
cuboid_points <- cuboid_points + 1
else if (eval(sphere_epsilon_condition)) # Lies close to sphere's surface
sphere_points <- sphere_points + 1
}
})
number_of_points_in_a_grid_of_that_size_on_the_total_surface_of_the_sphere <-
# TODO: Fill this in. Shouldn't be too hard.
cat("The ratio of the intersection of the surface of the sphere with the cuboid to the",
"total area of the sphere is approximately: ",
sphere_points / number_of_points_in_a_grid_of_that_size_on_the_total_surface_of_the_sphere)
