编写程序，实现如下功能：给定平面上两个圆，根据圆心位置和半径，判别它们是否是公切线。如果有的话，绘制两个圆的所有公切线。 要求：

1、用鼠标交互绘制两个圆，程序获取圆心坐标和半径；

2、通过对两个圆心连线的几何变换绘制所有的公切线。

利用matlab功能检测鼠标按下，鼠标左键描点，当已有三个或以上的点时鼠标右键以前三点画圆，当有两个圆时鼠标中键做出公切线。通过已有函数disp_point获得鼠标坐标，利用三点计算圆心公式求出圆心坐标与半径后画圆，已知两圆圆心及半径根据点到直线距离公式列出二元二次方程，求解得到公切线方程并绘制。

function draw_lines

list11=[0 0];

list22=[0 0];

figure('WindowButtonDownFcn',@wbdcb)

ah = axes('SortMethod','childorder');

axis ([1 100 1 100])

function wbdcb(src,callbackevent)

if strcmp(get(src,'SelectionType'),'normal')

[x,y,str] = disp_point(ah);

text(x,y,str,'VerticalAlignment','bottom');

listk=[x y];

list11=[list11 listk];

drawnow

hold on;

elseif strcmp(get(src,'SelectionType'),'alt')

[x0,y0,r]=calk();

drawnow

hold on;

elseif strcmp(get(src,'SelectionType'),'extend')

qiexian(list22);

end

end

function [x,y,str] = disp_point(ah)

cp = get(ah,'CurrentPoint');

x = cp(1,1);y = cp(1,2);

str = ['(',num2str(x,'%0.3g'),', ',num2str(y,'%0.3g'),')'];

end

function [x0,y0,r]=calk()

A=[list11(1,3) list11(1,4)];%坐标

B=[list11(1,5) list11(1,6)];

C=[list11(1,7) list11(1,8)];

judge1=(B(1)-C(1))*(B(2)-A(2))-(B(1)-A(1))*(B(2)-C(2));

judge2=B(2)-A(2);

I1=0.5*(C(2)-A(2))*(B(2)-C(2))*(B(2)-A(2));

I2=0.5*(A(1)^2-B(1)^2)*(B(2)-C(2));

I3=0.5*(B(1)^2-C(1)^2)*(B(2)-A(2));

if(judge1==0)

fprintf('输入的三点构成一条直线，不能画一个圆。');

else

x0=(I1+I2+I3)/judge1;

if(judge2==0)

y0=-(B(1)-C(1))*(x0-0.5*(B(1)+C(1)))/(B(2)-C(2))+0.5*(B(2)+C(2));

else

y0=-(B(1)-A(1))*(x0-0.5*(B(1)+A(1)))/(B(2)-A(2))+0.5*(B(2)+A(2));

end

r=sqrt((x0-A(1))^2+(y0-A(2))^2);

theta=0:0.01:2*pi;

x=x0+r*cos(theta);

y=y0+r*sin(theta);

plot([x0,A(1)],[y0,A(2)],'g');

plot([x0,B(1)],[y0,B(2)],'g');

plot([x0,C(1)],[y0,C(2)],'g');

plot(x,y,'r','linewidth',2);

list22=[list22 x0 y0 r];

end

list11=[0 0];

end

function qiexian(listkk)

x1=listkk(1,3);

y1=listkk(1,4);

r1=listkk(1,5);

x2=listkk(1,6);

y2=listkk(1,7);

r2=listkk(1,8);

d = sqrt((x1-x2)^2+(y1-y2)^2);

d1 = abs(r1-r2);

d2 = r1+r2;

if d < d1

error('大圆包含小圆，没有公切线')

end

if isequal([x1 y1 r1],[x2 y2 r2])

error('两圆重合，有无穷多条公切线')

end

syms k b

[k,b]=vpasolve([((k*x1-y1+b)^2)/(k*k+1)==r1^2,((k*x2-y2+b)^2)/(k*k+1)==r2^2],[k,b]);

w=-50:200;

q=k(1)*w+b(1);

plot(w,q);

q=k(2)*w+b(2);

plot(w,q);

q=k(3)*w+b(3);

plot(w,q);

q=k(4)*w+b(4);

plot(w,q);

end

end