How to calculate the volume of a segment of spherical cap intersected by a plane ?

Hi,

I want to calculate the volume of a spherical cap, intersected with a plane perpendicular to the spherical cap.

https://au.mathworks.com/matlabcentral/fileexchange/18641-spheresegmentvolume

The above code can calculate the volume of the spherical cap. Any idea how to calculate a segment of volume of that cap intersected by a plane. Please refer to the following figure.

Thanks in advance.

[James Tursa]

Arbitrary plane? Plane passing through center of sphere?

[Sadeep Thilakarathna]

Arbitrary plane. But, it is vertical. Does not need to be a slanted one.

[Prudhvi Peddagoni]

Hi,

When there is only one plane and we need to find volume of hemisphere, we integrate the area of circle from x to R,

syms x R r
expr = pi*r^2;                      % area of a circle with r as radius
actualRadius = sqrt(R^2-x^2);       % radius of circle if a plane intersects 
                                    % a sphere of radius R at a distance x from the center
expr = subs(expr,r,actualRadius);   % substituting the radius
result = int(expr,x,x,R);         % integrate the expression 

But if you want a section of the hemisphere, you need to integrate the part of a circle instead of the whole circle.

% Area of a section of hemisphere cut by planes x=x1, y=y1. The radius of sphere is R1
function result = volumehem(x1,y1,R1)
syms x y r R
R=R1;
expr = 2*sqrt(r^2-x^2);         % calculating the area of a section created 
                                % by a line which is at a distance of y from the center
p = int(expr,x,y, r);
p = subs(p,r,sqrt(R^2-x^2));
result = int(p,x,x,R);
result = subs(result,x,x1);
result = subs(result,y,y1);
end

But the issue is, the result becomes very complex. so you might have to use vpa function to evaluate.

Hope this helps

[Sadeep Thilakarathna]

Thanks for the answer. However, when evaluated the result becomes a complex number. I found the algorithm listed in the following paper useful.

https://www.sciencedirect.com/science/article/pii/S0021999116000577

Answer 1

[John D'Errico]

I did suggest this scheme in my response to your question. And it seems pretty easy to do, so I am not sure where the problem lies.

For example, consider the sphere of radius 4, centered arbitrarily at the origin. If the sphere is centered at any other point, then a translation is trivial.

Now I wish to find the volume of the partial cap, where x >= 1, and y >= 2. I will do this by reducing the problem to a numerical integration, of caps on a sequence of circles. Essentially, integrate from x == 1 to 4, As x increases, then we can visualize a plane that cuts the sphere, resulting in a circle. That circle will have radius that decreases as a function of x, but each of those circles will also be cut by a line, at fixed values of y. spheresegmentvolume can compute the area of the circular caps, and then we will just integrate that result.

R = 4;
Y0 = 2;
X0 = 1;
function caparea = cap(x,X0,Y0,R)
  r = @(x) sqrt(R^2 - x.^2);
  caparea = zeros(size(x));
  for i = 1:numel(x)
     caparea(i) = spheresegmentvolume([Y0,r(x(i))],2,r(x(i)));
  end
end

integral(@(x) cap(x,X0,Y0,R),X0,4)
ans =
          11.4642080722576

Is this correct? A simple test is to use a Monte Carlo. A million points should give me a few correct digits in the volume.

n = 1e6;
el = asin(2*rand(n,1)-1);
az = 2*pi*rand(n,1);
rad = R*nthroot(rand(n,1),3);
[x,y,z] = sph2cart(az,el,rad);

k = (x >= X0) & (y >= Y0);

approxvol = 4/3*pi*R^3 * sum(k)/n
approxvol =
          11.4420323027512

So it seems to have worked well enough. WTP?

[Sadeep Thilakarathna]

Thanks. Works perfectly.

[John D'Errico]

Note that it is not that difficult to compute the area inside a circular cap analytically. And that would speed things up.