I don't know what is your D3. So go to the bottom of the code and set D3 to the one that you have.
% Area of a circular segment
A_CS=@(theta,r) ( (r^2/2)*(theta-sin(theta)) );
% Center of mass of a circular segment
X_CS=@(theta,r) ( (4*r*sin(theta/2)^3)/(3*(theta-sin(theta))) );
% x/y of circle-circle intersect
% It is assumed the coordinate system is at the center of the circle on the right
% r1 is the radius of the circle on the right,
% r2 is the radius of the circle on the left.
% d is the distance between the two centers.
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));
%%SET THESE NUMBERS ACCORDINGLY
r1=50/2; %[50mm]
r2=52/2; %[52mm]
D3=40; %[80mm]
%solving for D1
[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')
Just note that I don't check if the r1 and r2 and D3 are valid numbers. I am assuming that they are valid. For example if D3 > (r1+r2) there is no solution at all and this program generates an error or gives a wrong answer. So make sure that the numbers are set properly. Also if the left circle is too much inside the the right circle this code gives you wrong answer. Make sure that D1<D3 at the end. and xInt(r1,r2,D1)>0;
