科学计算：分形几何_随机科赫雪花

%********************************************************************/

%* 科学计算：分形几何_随机科赫雪花                                                         */

%* 挑灯看剑[e-mail] 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

运行结果测试如下：