How to calculate the volume of a segment of spherical cap intersected by a plane ?
12 ビュー (過去 30 日間)
古いコメントを表示
Sadeep Thilakarathna
2021 年 1 月 14 日
コメント済み: 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.
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 件のコメント
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
採用された回答
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?
2 件のコメント
John D'Errico
2021 年 1 月 19 日
Note that it is not that difficult to compute the area inside a circular cap analytically. And that would speed things up.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Surface and Mesh Plots についてさらに検索
製品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!