Matlab 版 (精华区)

发信人: zjliu (秋天的萝卜), 信区: Matlab
标  题: Re: 求外切圆
发信站: BBS 哈工大紫丁香站 (Wed Nov 17 22:05:48 2004)

今天又改了一下程序，现在可以迭代到n=2，但是对于n=3还是不行，
excircle(2,2,0,2,0,6,sqrt(6^2+2^2)-2,-2,0,2);
是可以的，这样
excircle(3,2,0,2,0,6,sqrt(6^2+2^2)-2,-2,0,2);
就不行了，以下是程序：

function excircle(n,x1,y1,r1,x2,y2,r2,x3,y3,r3);

% Examples:
%      excircle(2,2,0,2,0,6,sqrt(6^2+2^2)-2,-2,0,2);
% See also inline
% close all

W=[x1,x2,x3;y1,y2,y3;r1,r2,r3];
t=0;
for p=1:n;
    Q=[];
    for g=1:(size(W,2)/3);
        K=W(:,((g-1)*3+1):((g-1)*3+3));
        C=reshape(K,1,9);
        x1=C(1);
        y1=C(2);
        r1=C(3);
        x2=C(4);
        y2=C(5);
        r2=C(6);
        x3=C(7);
        y3=C(8);
        r3=C(9);
        [x,y,r]=ivanhit2p(x1,y1,r1,x2,y2,r2,x3,y3,r3);
        t=t+1
        Cq=[x,y,r]';
        Q=[Q,K(:,1),K(:,2),Cq,K(:,1),K(:,3),Cq,K(:,3),K(:,2),Cq];
    end
    clear W
    W=Q;
end     
figure;
hold on
t=0:.1:2.1*pi;
for z=1:size(W,2);
plot(W(1,z)+W(3,z)*cos(t),W(2,z)+W(3,z)*sin(t),'black');
end
axis equal

function [x,y,r]=ivanhit2p(x1,y1,r1,x2,y2,r2,x3,y3,r3);
% 三个圆的外切圆算法
% 已知两两外切的三个圆(圆心和半径已知),求这三个圆的外切圆
% (即求圆心坐标和半径)
% Examples:
%       [x,y,r]=ivanhitp(2,0,2,0,6,sqrt(6^2+2^2)-2,-2,0,2);
%       [x,y,r]=ivanhit2p(2,0,2,-2,0,2,0,1.2891,0.3794);
% See also inline
% close all

T1=[x1,x2,x3];
T2=[y1,y2,y3];
T3=[r1,r2,r3];
[T2,v]=sort(T2);
T1=T1(v);
T3=T3(v);
x1=T1(1);
x2=T1(3);
x3=T1(2);
y1=T2(1);
y2=T2(3);
y3=T2(2);
r1=T3(1);
r2=T3(3);
r3=T3(2);

FF=inline('sqrt((x-x1)^2+(y-y1)^2)-r1','x','y','x1','y1','r1');
x0=mean([x1,x3]);
y0=mean([y1,y3]);
xp=x2-x0;
yp=y2-y0;
A=atan((y3-y1)/(x3-x1));
K=[cos(A),sin(A);-sin(A),cos(A)]*[xp;yp];
x22=K(1);
y22=K(2);
c=sqrt((x1-x3)^2+(y1-y3)^2)/2;
a=(r1-r3)/2*sign(x3-x1);
b=sqrt(c^2-a^2);
xf=inline('a*sqrt(y^2/b^2+1)','y','a','b');
ymin=0;    
ymax=y22;
y=mean([ymin,ymax]);
x=xf(y,a,b);
x11=c*sign(x1-x3);
y11=0;
tf1=FF(x,y,x11,y11,r3);
tf2=FF(x,y,x22,y22,r2);
while abs(tf1-tf2)>10^(-2);
    if tf1>tf2;
        ymax=y;
    else
        ymin=y;
    end
    y=mean([ymin,ymax]);
    x=xf(y,a,b);
    tf1=FF(x,y,x11,y11,r1);
    tf2=FF(x,y,x22,y22,r2);
end
K=[cos(-A),sin(-A);-sin(-A),cos(-A)]*[x;y];
x=K(1);
y=K(2);
x=x+x0;
y=y+y0;
r=tf1;
% 
% t=0:.1:2.1*pi;
% % figure;
% plot(x1+r1*cos(t),y1+r1*sin(t));
% hold on
% plot(x2+r2*cos(t),y2+r2*sin(t));
% plot(x3+r3*cos(t),y3+r3*sin(t));
% plot(x+r*cos(t),y+r*sin(t));
% axis equal
【 在 ivanhit (ivan) 的大作中提到: 】
: 不愧是大斑竹，您先忙，谢谢！！！



--
╔═══════════════════╗
║★★★★★友谊第一  比赛第二★★★★★║
╚═══════════════════╝


※ 来源:·哈工大紫丁香 http://bbs.hit.edu.cn·[FROM: 202.118.229.*]