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

11 次查看（过去 30 天）
Sadeep Thilakarathna
评论： John D'Errico 2021-1-19
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.
  4 个评论
显示 2更早的评论
Prudhvi Peddagoni
Prudhvi Peddagoni 2021-1-18
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 2021-1-19
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

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 3-D Volumetric Image Processing 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by