%********************************************************************/
%* 科学计算:分形几何_随机科赫雪花 */
%* 挑灯看剑-shuchangs@126.com 2010-5 */
%* 云歌国际(Cloud Singers International) www.cocoral.com */
%********************************************************************/
function RandKoch
k=3; %Koch曲线的迭代次数
p=[0 0;10 0]; %存放结点坐标,每行一个点,初始值为两节点的坐标
n=1; %存放线段的数量,初始值为1
AU=[cos(pi/3) sin(pi/3);-sin(pi/3) cos(pi/3)]; %计算新的节点,向外(上)生长
AD=[cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3)]; %计算新的节点,向里(下)生长
c=1;%图像序号计数器
for s=1:k %计算所有的节点坐标
j=0;
for i=1:n %每条边计算一次,每条边会新增3个点
q1=p(i,:); %目前线段的起点坐标
disp('目前线段的起点坐标');disp(q1);
q2=p(i+1,:); %目前线段的终点坐标
disp('目前线段的终点坐标');disp(q2);
d=(q2-q1)/3;
j=j+1;r(j,:)=q1; %原起点存入r
disp('起点坐标(x,y)');disp(r(j,:));
%***************************************
j=j+1;r(j,:)=q1+d; %新1点存入r
disp('新1点坐标(x,y)');disp(r(j,:));
%***************************************
%---------------------------------------
%选择变换矩阵
if coinTest()==1
A=AU;
elseif coinTest()==-1
A=AD;
end
disp('向量d阵[1*2]为');disp(d);
disp('向量d的变换阵:A阵[2*2]');disp(A);
disp('d*A阵[1*2]');disp(d*A);
%显示向量d的变换过程图
%dM=[0 0;d];
%dT=[0 0;d*A];
%subplot(5,2,c),plot(dM(:,1),dM(:,2)),title('向量d'),
%subplot(5,2,c+1),plot(dT(:,1),dT(:,2)),title('变换后的向量d');
%c=c+2;
%pause(2);%延迟2秒
j=j+1;r(j,:)=q1+d+d*A;%新2点存入r
disp('新2点坐标(x,y)');disp(r(j,:));
%---------------------------------------
%***************************************
j=j+1;r(j,:)=q1+2*d; %新3点存入r
disp('新3点坐标(x,y)');disp(r(j,:));
%***************************************
disp('暂存坐标点(x,y)');disp(r);
disp('*********************************')
end
n=4*n;
clear p
p=[r;q2];
end
plot(p(:,1),p(:,2));
disp('科赫曲线坐标+--x--+--y--+');disp(p);
axis equal
end
%定义选择1或-1函数
function rtnV=coinTest() %抛币试验函数,返回值为1(正)或-1(反)
x=rand;
if x<0.5
rtnV=-1;
elseif x>0.5
rtnV=1;
else
rtnV=coinTest();
end
end
运行结果测试如下: