MATLAB Answers

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

38 views (last 30 days)
Hi,
I want to calculate the volume of a spherical cap, intersected with a plane perpendicular to the spherical cap.
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.

  4 Comments

Show 1 older comment
Prudhvi Peddagoni
Prudhvi Peddagoni on 18 Jan 2021
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
Sadeep Thilakarathna on 19 Jan 2021
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

Sign in to comment.

Accepted Answer

John D'Errico
John D'Errico on 19 Jan 2021
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?

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by