square.m

这里面有很多有用的东东

function square(n)
%SQUARE
%Draws a fractal square.
%Calling format: square(n), where n is desired # iterations. 

%Copyright Gareth Williams, Stetson University 
[e-mail], http://www.stetson.edu/~gwilliam
%Accompanies "Linear Algebra with Applications" by Gareth Williams


disp(' ')
disp('   - To Stop At Any Time -') 
disp('    command period (Mac)')
disp('       Ctrl C (IBM)')
disp(' ')
disp('- press Enter to get picture - ')
disp(' ')
pause

clf;

%matrices of affine transformations e.g. AX + E
A=[.5 0;0 .5];E=[.5;.5];
B=[.5 0;0 .5];F=[0;.5];
C=[.5 0;0 .5];G=[.5;0];
D=[.5 0;0 .5];H=[0;0];
%accumulated probabilities
P=[.25 .5 .75 1];
%Ranges of x & y axes
V=[-.1 1.1 -.1 1.1];
axis(V);
%Plot the points generated by the affine transforms
hold on
X=[0;0];Y=[0;0];
figure(gcf)
for i=1:n,
 pk=rand;
 if pk<=P(1)
   Y(1,1) = A(1,1)*X(1,1)+A(1,2)*X(2,1)+E(1,1);
   Y(2,1) = A(2,1)*X(1,1)+A(2,2)*X(2,1)+E(2,1);
  elseif pk<=P(2)
   Y(1,1) = B(1,1)*X(1,1)+B(1,2)*X(2,1)+F(1,1);
   Y(2,1) = B(2,1)*X(1,1)+B(2,2)*X(2,1)+F(2,1); 
  elseif pk<=P(3) 
   Y(1,1) = C(1,1)*X(1,1)+C(1,2)*X(2,1)+G(1,1);
   Y(2,1) = C(2,1)*X(1,1)+C(2,2)*X(2,1)+G(2,1);
  elseif pk<=P(4)
   Y(1,1) = D(1,1)*X(1,1)+D(1,2)*X(2,1)+H(1,1);
   Y(2,1) = D(2,1)*X(1,1)+D(2,2)*X(2,1)+H(2,1);
 end
 plot(Y(1,1),Y(2,1),'.')
 X=Y;
end

hold off