How to calculate the volume of a segment of spherical cap intersected by a plane ?
12 views (last 30 days)
Show older comments
Sadeep Thilakarathna
on 14 Jan 2021
Commented: John D'Errico
on 19 Jan 2021
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/487628/image.png)
Thanks in advance.
4 Comments
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
Accepted Answer
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?
2 Comments
John D'Errico
on 19 Jan 2021
Note that it is not that difficult to compute the area inside a circular cap analytically. And that would speed things up.
More Answers (0)
See Also
Categories
Find more on Surface and Mesh Plots in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!