function dandelin(r1,r2,alpha)
% DANDELIN(R1,R2,ALPHA) reprŽsente 2 sph�res de Dandelin dans un c™ne.
% EXEMPLE: dandelin(1.4,3.4,70)

% Jean-Bernard Roux, coll�ge de Saussure, 9.02

% Initialisation
m=tan(pi*alpha/180);    % pente du c™tŽ du c™ne
h1=r1*sqrt(1+m^2);      % hauteur de la sph�re no 1
h2=r2*sqrt(1+m^2);      % hauteur de la sph�re no 2
C1=[0 0 h1];              % centres des sph�res
C2=[0 0 h2];

% M est le point d'intersction entre l'axe du c™ne et la tangente
% aux deux sph�res
MC1=(r1*(h2-h1))/(r1+r2);   % distance du centre C1 au point M
MC2=MC1*(r2/r1);            % distance du centre C2 au point M
M=[0 MC1+h1];

% intersctions sph�res avec le plan
T1=C1+[0 -r1*sqrt(1-((r1^2)/MC1^2)) (r1^2)/MC1];
T2=C2+[0 r2*sqrt(1-((r2^2)/MC2^2)) -(r2^2)/MC2];

% dessin du c™ne
ht=1.2*(h2+r2);
theta=linspace(pi/4,305*pi/180,40);
r=linspace(0,ht,40);
x=r'*cos(theta)./m;
y=r'*sin(theta)./m;
z=r'*ones(size(theta));
surfl(x,y,z)
hold on

% dessin des sph�res
theta=linspace(0,2*pi,40);
phi=linspace(0,pi,40);
r=r1-0.01;
x1=r.*cos(theta')*sin(phi);
y1=r.*sin(theta')*sin(phi);
z1=r.*ones(size(theta'))*cos(phi)+h1;
r=r2-0.01;
x2=r.*cos(theta')*sin(phi);
y2=r.*sin(theta')*sin(phi);
z2=r.*ones(size(theta'))*cos(phi)+h2;
surfl(x1,y1,z1)
surfl(x2,y2,z2)

% dessin du plan tangent
[X,Y]=meshgrid(-1.2*r2:0.5:1.2*r2,-1.2*r2:0.5:1.2*r2);
a=sqrt(m^2*((r2-r1)^2)-4*r1*r2)/(r1+r2);
b=(2*r1*r2*sqrt(1+m^2))/(r1+r2);
Z=a*Y+0*X+b;
surfl(X,Y,Z)

% dessin de l'ellipse et des cercles
theta=linspace(0,2*pi,60);
r=b./(1-a*sin(theta)./m);   % intersection du c™ne et du plan
x3=r.*cos(theta)./m;
y3=r.*sin(theta)./m;
z3=r;
plot3(x3,y3,z3,'linewidth',5)
t=linspace(0,2*pi,50);  % cercle autour de la sph�re 1
x4=1.07*r1*cos(t);      % 1.07 pour agrandir le cercle      
y4=1.07*r1*sin(t);
z4=h1*ones(size(t));
plot3(x4,y4,z4,'linewidth',5)
t=linspace(0,2*pi,50);  % cercle autour de la sph�re 2
x5=1.07*r2*cos(t);
y5=1.07*r2*sin(t);
z5=h2*ones(size(t));
plot3(x5,y5,z5,'linewidth',5)


% dessin des foyers de l'ellipse
c=(T1+T2)./2;   % point milieu entre les points tangents
k=0.2;          % facteur de rŽduction de l'ellipse
xr=k*(x3-c(1)); % on dŽplace l'ellipse ˆ l'origine
yr=k*(y3-c(2)); % on rŽduit l'ellipse et on
zr=k*(z3-c(3)); % la remonte
fill3(T1(1)+xr,T1(2)+yr,T1(3)+zr+0.01,'m')
fill3(T2(1)+xr,T2(2)+yr,T2(3)+zr+0.01,'m')
view(84,10)
axis equal
axis off
hold off
