matlab 两个球体相交的体积

主要参考了一个网上的一些方法下面是原理

设 d 为两球球心之间的距离。

1、d >= R+ri ：两球不相交，即相交部分体积=0.

2、d+ri = R ：小球在大球里面，即相交部分体积 = 小球体积 = $\frac{4}{3}*\pi* (ri)^{3}$ .

3、R-ri < d < R+ri ：两球相交，相交部分体积：

设 $cos\alpha = \frac{R^{2}+d^{2}-(ri)^{2}}{2*R*d}$ ， $h1 = R*(1-cos\alpha )$ .

$cos\beta = \frac{(ri)^{2}+d^{2}-R^{2}}{2*ri*d}$ ， $h2 = ri*(1-cos\beta )$ .

V = $\frac{\pi }{3}*(3*R-h1)*h1^{2}+\frac{\pi }{3}*(3*ri-h2)*h2^{2}$ .

证明：

相交部分体积是由2块构成的，分别属于两个球体。其中一块的体积公式为（以大球为例）：

$V=\int_{0}^{h}\pi (R^{2}-(Rcos\alpha +x)^{2})dx$

$=\int_{0}^{h}\pi (R^{2}-R^{2}cos^{2}\alpha -x^{2}-2Rcos\alpha x)dx$

$=\int_{0}^{h}\pi (R^{2}sin^{2}\alpha -x^{2}-2Rcos\alpha x)dx$

$= \pi (- \frac{1}{3}h^{3}-Rcos\alpha h^{2}+R^{2}sin^{2}\alpha h)$

$= \pi (Rh(Rsin^{2}\alpha -hcos\alpha )- \frac{1}{3}h^{3})$

$= \pi (Rh(R(1-cos^{2}\alpha ) -hcos\alpha )- \frac{1}{3}h^{3})$

$= \pi (Rh(R-Rcos^{2}\alpha -hcos\alpha )- \frac{1}{3}h^{3})$

$= \pi (Rh(R-cos\alpha (Rcos\alpha +h) )- \frac{1}{3}h^{3})$

$= \pi (Rh(R-Rcos\alpha )- \frac{1}{3}h^{3})$

$=\pi (Rh^{2}- \frac{1}{3}h^{3})$

代码：

 1 function [] = sphere()    
 2 clc
 3 clear
 4  pi=3.1415926;
 5 Circle_A.x=0;
 6 Circle_A.y=0;
 7 Circle_A.z=0;
 8 Circle_A.r=5;
 9 % % % % % % % % % % 
10 Circle_B.x=0;
11 Circle_B.y=0;
12 Circle_B.z=5;
13 Circle_B.r=5;
14 int16 T=20;
15 int16 Case=1;
16 V=0.0;
17 d=sqrt((Circle_A.x - Circle_B.x)*(Circle_A.x - Circle_B.x) + (Circle_A.y - Circle_B.y)*(Circle_A.y - Circle_B.y) + (Circle_A.z - Circle_B.z)*(Circle_A.z - Circle_B.z));
18 if(d>=Circle_A.r+Circle_B.r)
19 V=0;
20 elseif(d+Circle_A.r<=Circle_B.r)
21 V=V+(4.0 / 3)*pi*Circle_A.r*Circle_A.r*Circle_A.r;
22 elseif(d+Circle_B.r<=Circle_A.r)
23 V=V+(4.0 / 3)*pi*Circle_B.r*Circle_B.r*Circle_B.r;
24 else
25 co = double((Circle_A.r*Circle_A.r+d*d-Circle_B.r*Circle_B.r)/(2.0*d*Circle_A.r));
26 h = double(Circle_A.r*(1-co));
27 V=V+(1.0/3)*pi*(3.0*Circle_A.r-h)*h*h;
28 co=double((Circle_B.r*Circle_B.r+d*d-Circle_A.r*Circle_A.r)/(2.0*d*Circle_B.r));
29 h=double(Circle_B.r*(1-co));
30 V=V+(1.0/3)*pi*(3.0*Circle_B.r-h)*h*h;
31 end
32 V