clear;
theta=pi*1.5;
a=[0.5 0.5 0.5];
b=[0 0 0];
c=b;
d=a;
 e=[0 1 0.5];
 f=[0 0 0.5];
 
 k1=1/2; k2=1/2;
  a1=a;
 b1=c*k1/k2;
 c1=b*k2/k1;
 d1=d;
 e1=e*k1;
 f1=f*k2;
 
 p=[0.33 0.66 1.0];
            n=5000;            
     x=0.0; y=0.0;       
            for l=1:1:n 
                pk=rand;
                 if pk <= p(1) 
                    k=1;
                 elseif pk <= p(2)
                    k=2;
                 else         k=3;
                  end
                 xnew(l) =  a1(k)*x+ b1(k)*y + e1(k);
                 ynew(l) =  c1(k)*x+ d1(k)*y + f1(k);
                
                    x=xnew(l);
                    y=ynew(l);
                    
                              
          end
          
          subplot(2,2,4); 
          plot(xnew,ynew,'.');
          xlabel('(d)');
          % axis equal;
           axis([0 4 0 4]);
          %  axis([-2 2.2 -2 2.2]);

