Volume of a part of a sphere cutted using 3 planes using R-Language (integral)

Question

I want to calculate the volume (V) of a part of sphere, which is the result of intersetion of the sphere with three palnes (x=0, y=0 and z=1.5). I am using R-Language and this is my code. I tried 2 different methods using cartesian and polar coordinates. Both of them deliver negative answers.

##  Compute the Volume between 3 planes x=0, y=0 and z=1.5 and a sphere
library("pracma", lib.loc="~/R/win-library/3.1")
f <- function(x, y)  (sqrt(4 -x^2 - y^2) - 1.5 ) # here the function(x,y)  is subtracted with -1.5 which represents the plane z=1.5 
xmin <- 0
xmax <- 2
ymin <- 0 
ymax <- function(x) (sqrt(4 - x^2))
I <- integral2(f, xmin, xmax, ymin, ymax)
I$Q   # Integral over sector is not so exact

# Exact Volume from AutoCAD V=0.3600



##  Volume of the sphere: use polar coordinates
f0 <- function(x, y) (sqrt(4 - x^2 - y^2)-1.5)  # for x^2 + y^2 <= 4 the f(x,y) means z changes between zmin=1 and zmax= sqrt(4-x^2-y^2)
fp <- function(t, r) r * f0(r*cos(t), r*sin(t))
quad2d(fp, 0, pi/2, 0, 2, n = 101)  # -0.523597

The correct answer is V= 0.3600 . Can anyone give a me hint, please?

Cheers

Answer 1

Your X-Y region of integration covers areas where f(x,y)-1.5 is negative, as well as positive. The intersection of your sphere with the line z=1.5 is a circle of radius sqrt(7/4) (using Pythagoras), so adjusting your limits appropriately, you get:

library(pracma)
f <- function(x, y)  (sqrt(4 -x^2 - y^2) - 1.5 ) # here the function(x,y)  is subtracted with -1.5 which represents the plane z=1.5 
xmin <- 0
xmax <- sqrt(7/4)
ymin <- 0 
ymax <- function(x) (sqrt(7/4 - x^2))
I <- integral2(f, xmin, xmax, ymin, ymax)
I$Q   # Integral over sector is not so exact
# [1] 0.3599741

Pretty close to what you're expecting.

Answer 2

I have solved it for the sphere with r=2 and the volume of intersection between sphere shell and 3 planes x=1,y=1 and z=1.

##  Compute the Volume between 3 planes x=1.0, y=1.0 and z=1.0 and a sphere
library(pracma)
f <- function(x, y)  (sqrt(4 -x^2 - y^2) - 1 ) # here the function(x,y) is  subtracted with -1.5 which represents the plane z=1.5 
xmin <- 1
xmax <- sqrt(2)
ymin <- 1 
ymax <- function(x) (sqrt(3 - x^2))
I <- integral2(f, xmin, xmax, ymin, ymax)
I$Q   # 
# [1] 0.01520549
# Exact Volume from AutoCAD: 0.0152