function squares(n);
% 謝爾賓斯基“地毯”;
% n 為遞歸次數(shù)，亦即分形的層數(shù);

clc;close all;clear;
if nargin==0;       % 設(shè)定遞歸次數(shù)n
    n=3;
end
rand('state',3);    % 將當(dāng)前狀態(tài)重設(shè)到狀態(tài)2，
                    % 使每次程序運(yùn)行時(shí)生成的各個(gè)分形正方形的顏色不變
C=rand(n+5,3);      % 生成一個(gè)(n+4)*3的隨機(jī)矩陣C，對(duì)各個(gè)分形填色
figure;
axis square equal;  % 保持縱、橫坐標(biāo)的刻度單位一致
axis([-1 9.5 -1.2 1.2]); % 設(shè)定標(biāo)度范圍
grid;hold on;
a=-pi/4;            % 設(shè)定正方形的偏轉(zhuǎn)角度。a=-pi/4時(shí)，正方形有一條邊平行于x軸
p=-2;                % 設(shè)定正方形中心位置。p=j時(shí)，中心位置為(j,0)
r=sqrtm(2);                % 設(shè)定正方形的半徑，即中心到頂角的距離
[p,r,n,a]=squsqu(p,r,n,a,C);

function [p,r,n,a]=squsqu(p,r,n,a,C);
% 繪制正方形，該函數(shù)通過(guò)遞歸實(shí)現(xiàn)對(duì)正方形的分形;
% p 是正方形的中心位置;
% r 是正方形的半徑;
% n 是遞歸次數(shù);
% a 是正方形的偏轉(zhuǎn)角度;
% C 是顏色矩陣;

s=2.5;              % 設(shè)定正方形的中心偏移量,從而可得到分形的生成過(guò)程
z=p+r*exp(i*([0:4]*pi*1/2+a))+s;        % z 是主正方形的中心位置
z1=p+r*exp(i*([0:4]*pi*1/2))/2+s;     % zr 是主正方形內(nèi)各個(gè)分形的中心位置
z2=p+sqrtm(2)*r*exp(i*([0:4]*pi*1/2+a))/2+s;
pf=fill(real(z),imag(z),C(n+2,:));      % 生成一個(gè)已填色的正方形
set(pf,'EdgeColor',C(n+2,:));           % 設(shè)定正方形邊緣的顏色
if n>0;
%   [p,r,n,a]=squsqu(p,r/2,n-1,a+pi/3,C);       % 繪制主正方形中心的分形，     
%   n=n+1;r=r*2;a=a-pi/3;
%   正方形中這一部分是挖空的，故這兩行程序可省去
    [z1(1),r,n,a]=squsqu(z1(1),r/3,n-1,a,C);    % 繪制主正方形右下角的分形
    n=n+1;r=r*3;
    [z1(2),r,n,a]=squsqu(z1(2),r/3,n-1,a,C);    % 繪制主正方形上部的分形
    n=n+1;r=r*3;
    [z1(3),r,n,a]=squsqu(z1(3),r/3,n-1,a,C);    % 繪制主正方形左下角的分形
    n=n+1;r=r*3;
    [z1(4),r,n,a]=squsqu(z1(4),r/3,n-1,a,C);    % 繪制主正方形左下角的分形
    n=n+1;r=r*3;
    [z2(1),r,n,a]=squsqu(z2(1),r/3,n-1,a,C);    % 繪制主正方形右下角的分形
    n=n+1;r=r*3;
    [z2(2),r,n,a]=squsqu(z2(2),r/3,n-1,a,C);    % 繪制主正方形上部的分形
    n=n+1;r=r*3;
    [z2(3),r,n,a]=squsqu(z2(3),r/3,n-1,a,C);    % 繪制主正方形左下角的分形
    n=n+1;r=r*3;
    [z2(4),r,n,a]=squsqu(z2(4),r/3,n-1,a,C);    % 繪制主正方形左下角的分形
    n=n+1;r=r*3;
end