%% Intersection de deux spheres
clear all;clf

n=40;
theta=(-n:2:n)/n*pi;
phi=(-n:2:n)'/n*pi/2;
cosphi=cos(phi);cosphi(1)=0;cosphi(n+1)=0;
sintheta=sin(theta);sintheta(1)=0;sintheta(n+1)=0;

x=cosphi*cos(theta);
y=cosphi*sintheta;
z=sin(phi)*ones(1,n+1);

RSat=0.05;

figure(1);clf;whitebg('w');hold on
view(3);
C1=[-0.5 0 1.5];R1=0.8;
Sat1=surf(RSat*x+C1(1),RSat*y+C1(2),RSat*z+C1(3),10*ones(size(z)));
S1=surf(R1*x+C1(1),R1*y+C1(2),R1*z+C1(3),8*ones(size(z)));

C2=[0.5 0 1.5];R2=0.9;
Sat2=surf(RSat*x+C2(1),RSat*y+C2(2),RSat*z+C2(3),10*ones(size(z)));
S2=surf(R2*x+C2(1),R2*y+C2(2),R2*z+C2(3),8*ones(size(z)));


shading interp
alpha(0.2)
alpha(Sat1,1)
alpha(Sat2,1)
camlight;
lighting gouraud;
axis square equal


tht=0:0.01:2*pi;
xI=(R2*C1(1)+R1*C2(1))/(R1+R2);xI=-0.1;
yI=(R2*C1(2)+R1*C2(2))/(R1+R2);%yI=0;
zI=(R2*C1(3)+R1*C2(3))/(R1+R2);%zI=1.5;

D1=(xI-C1(1))^2+(yI-C1(2))^2+(zI-C1(3))^2;
D1=sqrt(D1);
D2=(xI-C2(1))^2+(yI-C2(2))^2+(zI-C2(3))^2;
D2=sqrt(D2);
R=R1^2-D1^2;R=sqrt(R);
%R=R2^2-D2^2;R=sqrt(R)


plot3(xI,yI,zI,'.r')
pI=plot3(xI+zeros(size(tht)),yI+R*cos(tht),zI+R*sin(tht),'-r');
set(pI,'linewidth',2.5)

vw=[0.9681    0.2504    0.0000   -0.6093
   -0.1973    0.7629    0.6157   -0.5906
   -0.1541    0.5961   -0.7880    8.8333
         0         0         0    1.0000];

view(vw);
axis off