A_CS=@(theta,r) ( (r^2/2)*(theta-sin(theta)) );
X_CS=@(theta,r) ( (4*r*sin(theta/2)^3)/(3*(theta-sin(theta))) );
xInt=@(r1,r2,d) ( (d^2-r2^2+r1^2)/(2*d) );
yInt=@(r1,r2,d) ( sqrt( (4*d^2*r1^2-(d^2-r2^2+r1^2)^2)/(4*d^2) ) );
theta1=@(x,y) atan(y/x);
theta2=@(x,y,d) atan( y / (d-x) );
Nom=@(r1,r2,d) ( d*pi*r2^2 ...
- A_CS(theta1(xInt(r1,r2,d),yInt(r1,r2,d)),r1)*X_CS(theta1(xInt(r1,r2,d),yInt(r1,r2,d)),r1) ...
- A_CS(theta2(xInt(r1,r2,d),yInt(r1,r2,d),d),r2)*(d-X_CS(theta2(xInt(r1,r2,d),yInt(r1,r2,d),d),r2)) );
Denom= @(r1,r2,d) ( pi*r2^2 ...
- A_CS(theta1(xInt(r1,r2,d),yInt(r1,r2,d)),r1) ...
- A_CS(theta2(xInt(r1,r2,d),yInt(r1,r2,d),d),r2) );
xHatch=@(r1,r2,d) (Nom(r1,r2,d)/Denom(r1,r2,d));
r1=50/2;
r2=52/2;
D3=40;
[D1,fval] = fsolve(@(d1) (D3-xHatch(r1,r2,d1)), D3-1 );
D2=D3-D1;
ezplot(@(x,y) (x.^2+y.^2-r1.^2),[-r1, D1+r2,-max(r1,r2),max(r1,r2)])
hold on
ezplot(@(x,y) ((x-D1).^2+y.^2-r2.^2),[-r1, D1+r2,-max(r1,r2),max(r1,r2)])
title('')
axis equal
line([xInt(r1,r2,D1) xInt(r1,r2,D1)],ylim,'Color','r')
Best Answer