%% Matlab script for testing the coordinates of touching points of a tangent on a circle. 
%% Circle are randomly selected, tangents are computed and displayed on the figure
%% Written by Philippe Lucidarme
%% http://www.lucidarme.me

close all;
clear all;
clc;

% Circle A
xA=randi(20)-10;
yA=randi(20)-10;
RA=randi(10);

% Circle B
xB=randi(20)-10;
yB=randi(20)-10;
RB=randi(10);


% Draw circles
a=[0:pi/16:2*pi];
Xa=xA+RA*cos(a);
Ya=yA+RA*sin(a);
plot (Xa,Ya,'k');
hold on;
axis square equal;
grid on;
plot (xA,yA,'xk');


Xb=xB+RB*cos(a);
Yb=yB+RB*sin(a);
plot (Xb,Yb,'k');
plot (xB,yB,'xk');


% Compute distance between circle centers
D=sqrt ( (xB-xA)^2 + (yB-yA)^2 );

%% First case : process external tangents

if (D^2 - (RA-RB)^2>=0)
    
    % Compute the length of the tangents
    L=sqrt(D^2 - (RA-RB)^2);

    % Compute the parameters
    R1=sqrt(L^2+RB^2);    
    Sigma1= (1/4) * sqrt ( ( D+RA+R1 )*( D+RA-R1 )*( D-RA+R1 )*( -D+RA+R1 )  );
    R2=sqrt(L^2+RA^2);
    Sigma2= (1/4) * sqrt ( ( D+RB+R2 )*( D+RB-R2 )*( D-RB+R2 )*( -D+RB+R2 )  );
    
    
    % Compute the first tangent
    x11= (xA+xB)/2 + (xB-xA)*(RA^2-R1^2)/(2*D^2) + 2*(yA-yB)*Sigma1/(D^2);
    y11= (yA+yB)/2 + (yB-yA)*(RA^2-R1^2)/(2*D^2) - 2*(xA-xB)*Sigma1/(D^2);
    x21= (xB+xA)/2 + (xA-xB)*(RB^2-R2^2)/(2*D^2) - 2*(yB-yA)*Sigma2/(D^2);
    y21= (yB+yA)/2 + (yA-yB)*(RB^2-R2^2)/(2*D^2) + 2*(xB-xA)*Sigma2/(D^2);   
    
    % Display tangent
    plot (x11,y11,'ob');
    plot (x21,y21,'ob');
    line ([x11 x21],[y11 y21],'Color','b');
       
    
    % Compute second tangent
    x12= (xA+xB)/2 + (xB-xA)*(RA^2-R1^2)/(2*D^2) - 2*(yA-yB)*Sigma1/(D^2);
    y12= (yA+yB)/2 + (yB-yA)*(RA^2-R1^2)/(2*D^2) + 2*(xA-xB)*Sigma1/(D^2);
    x22= (xB+xA)/2 + (xA-xB)*(RB^2-R2^2)/(2*D^2) + 2*(yB-yA)*Sigma2/(D^2);
    y22= (yB+yA)/2 + (yA-yB)*(RB^2-R2^2)/(2*D^2) - 2*(xB-xA)*Sigma2/(D^2);   
    
    % Display tangent
    plot (x12,y12,'om');
    plot (x22,y22,'om');
    line ([x12 x22],[y12 y22],'Color','m');
else
    disp ('No external tangents');
end






%% Second case : process internal tangents

if (D^2 - (RA+RB)^2>=0)
    
    % Compute the length of the tangents
    L=sqrt(D^2 - (RA+RB)^2);

    % Compute the parameters
    R1=sqrt(L^2+RB^2);    
    Sigma1= (1/4) * sqrt ( ( D+RA+R1 )*( D+RA-R1 )*( D-RA+R1 )*( -D+RA+R1 )  );
    R2=sqrt(L^2+RA^2);
    Sigma2= (1/4) * sqrt ( ( D+RB+R2 )*( D+RB-R2 )*( D-RB+R2 )*( -D+RB+R2 )  );
    
    
    % Compute the first tangent
    x11= (xA+xB)/2 + (xB-xA)*(RA^2-R1^2)/(2*D^2) + 2*(yA-yB)*Sigma1/(D^2);
    y11= (yA+yB)/2 + (yB-yA)*(RA^2-R1^2)/(2*D^2) - 2*(xA-xB)*Sigma1/(D^2);
    x21= (xB+xA)/2 + (xA-xB)*(RB^2-R2^2)/(2*D^2) + 2*(yB-yA)*Sigma2/(D^2);
    y21= (yB+yA)/2 + (yA-yB)*(RB^2-R2^2)/(2*D^2) - 2*(xB-xA)*Sigma2/(D^2);   
    
    % Display tangent
    plot (x11,y11,'og');
    plot (x21,y21,'og');
    line ([x11 x21],[y11 y21],'Color','g');
       
    
    % Compute second tangent
    x12= (xA+xB)/2 + (xB-xA)*(RA^2-R1^2)/(2*D^2) - 2*(yA-yB)*Sigma1/(D^2);
    y12= (yA+yB)/2 + (yB-yA)*(RA^2-R1^2)/(2*D^2) + 2*(xA-xB)*Sigma1/(D^2);
    x22= (xB+xA)/2 + (xA-xB)*(RB^2-R2^2)/(2*D^2) - 2*(yB-yA)*Sigma2/(D^2);
    y22= (yB+yA)/2 + (yA-yB)*(RB^2-R2^2)/(2*D^2) + 2*(xB-xA)*Sigma2/(D^2);   
    
    % Display tangent
    plot (x12,y12,'or');
    plot (x22,y22,'or');
    line ([x12 x22],[y12 y22],'Color','r');
else
    disp ('No internal tangents');
end