问题描述
编写程序,实现如下功能:给定平面上两个圆,根据圆心位置和半径,判别它们是否是公切线。如果有的话,绘制两个圆的所有公切线。 要求:
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
实验结果
![](https://i-blog.csdnimg.cn/blog_migrate/a3cb44252817fc3d85cf973d66819a6f.png)