function twodcell2(steps,sav,rule,top)
%for life twodcell2(200,0,1,2)
%for periodic demons twodcell2(200,0,2,8)
%for chemical waves use twodcell2(1000,0,3,30) % with g=7
%for time tunnel twodcell2(1000,0,4,30) 
%for fractal 2-D twodcell2(100,0,5,4)
%program will save sequence as jpgs if sav is 1
mysiz=100;
myh=mysiz;
myw=mysiz+20;
switch rule
    case 4
        mymat=zeros(myh,myw);
        mymat(floor(myh/2),floor(myw/2))=1;
        mymat(floor(myh/2),floor(myw/2)+1)=2;
        mymat(floor(myh/2),floor(myw/2)-1)=2;
        mymat(floor(myh/2)+1,floor(myw/2))=2;
        mymat(floor(myh/2)-1,floor(myw/2))=2;
        mymat(floor(myh/2)+1,floor(myw/2)+1)=1;
        mymat(floor(myh/2)-1,floor(myw/2)+1)=1;
        mymat(floor(myh/2)+1,floor(myw/2)-1)=1;
        mymat(floor(myh/2)-1,floor(myw/2)-1)=1;
    case 5
        mymat=zeros(myh,myw);
        mymat(floor(myh/2),floor(myw/2))=1;
    otherwise
        mymat=floor(top*rand(myh,myw));
end
mymat2=mymat;
show(mymat,0,sav,top);
myn=zeros(3);
for s=1:steps
    for i=1:myh
        for j=1:myw
            a=i-1;
            b=i+1;
            if i==1
                a=myh;
            elseif i==myh
                b=1;
            end
            l=j-1;
            r=j+1;
            if j==1
                l=myw;
            elseif j==myw
                r=1;
            end
            myn(1,1)=mymat(a,l);
            myn(1,2)=mymat(a,j);
            myn(1,3)=mymat(a,r);
            myn(2,1)=mymat(i,l);
            myn(2,2)=mymat(i,j);
            myn(2,3)=mymat(i,r);
            myn(3,1)=mymat(b,l);
            myn(3,2)=mymat(b,j);
            myn(3,3)=mymat(b,r);
            switch rule
                case 1
                    mymat2(i,j)=life(myn,top);
                case 2
                    mymat2(i,j)=demon(myn,top);
                case 3
                    mymat2(i,j)=infect(myn,top);
                case 4
                    mymat2(i,j)=time(myn,top);
                case 5

                    mymat2(i,j)=bbc(myn,top);            
            end
        end
    end
    mymat=mymat2;
    show(mymat,s,sav,top);
end

function show(mymat,s,sav,top)
[mylen width]=size(mymat);
myimat=zeros(mylen,width,3);
for i=1:mylen
    for j=1:width
        k=8*pi*mymat(i,j)/top;
        myimat(i,j,1) = abs(sin(k/10));
        myimat(i,j,2) = abs(sin(k/10+pi/4));
        myimat(i,j,3) = abs(cos(k/10));    
    end
end
%image([1 mylen],[1 width],myimat);
image(myimat);
getframe;
%F(i)=getframe;
if sav
    imwrite(myimat,strcat('im',num2str(s),'.jpg'),'jpg'); 
end

function out=infect(in,top)
A=0; B=0; S=0; g=7;
k1=2; k2=3;
top=top-1;
for i=1:3
    for j=1:3
        S=S+in(i,j);
        if in(i,j)>0
            A=A+1;
        end
        if in(i,j)==top
            B=B+1;
        end
    end
end
if in(2,2)>0 
    A=A-1;
end
if in(2,2)==top
    B=B-1;
end
if in(2,2)==0
    out=floor(A/k1)+floor(B/k2);
else
    out=floor(S/A)+g;
end
if in(2,2)==top
    out=0;
end
if out>top
    out=top;
end

function out=demon(in,top)
if in(1,2)==0
    in(1,2)=top;
end
if in(3,2)==0
    in(3,2)=top;
end
if in(2,1)==0
    in(2,1)=top;
end
if in(2,3)==0
    in(2,3)=top;
end
    
if in(1,2)==in(2,2)+1 || in(3,2)==in(2,2)+1 || in(2,1)==in(2,2)+1 || in(2,3)==in(2,2)+1
    out=in(2,2)+1;
else
    out=in(2,2);
end
if out==top
    out=0;
end

function out=life(in,top)
s=0;
for i=1:3
    for j=1:3
        s=s+in(i,j);
    end
end
s=s-in(2,2);
out=0;
if in(2,2)==1
    if s==2 || s==3
        out=1;
    end
elseif s==3
    out=1;
end

function out=time(in,top)
%in=uint8(in);
%s=mod(in(1,2),2)+mod(in(2,2),2)+mod(in(3,2),2)+mod(in(2,1),2)+mod(in(2,3),2);
%s=bitget(uint8(in(1,2)),1)+bitget(uint8(in(2,2)),1)+bitget(uint8(in(3,2)),1)+bitget(uint8(in(2,1)),1)+bitget(uint8(in(2,3)),1);
s=bitget(in(1,2),1)+bitget(in(2,2),1)+bitget(in(3,2),1)+bitget(in(2,1),1)+bitget(in(2,3),1);
if s>0 && s<5
    out=xor(in(2,2)>1,1);
else
    out=xor(in(2,2)>1,0);
end
%out=out+2*mod(in(2,2),2);
out=out+2*bitget(in(2,2),1); 

function out=bbc(in,top)
s=(in(1,2)~=0)+(in(2,2)~=0)+(in(3,2)~=0)+(in(2,1)~=0)+(in(2,3)~=0);
if in(2,2)~=0
    out=in(2,2)-1;
else
    if bitand(uint8(s),1)
        out=3;
    else
        out=0;
    end
end