clear;

%Nx=input('Kastenlänge in x:');
Nx=21;
%Ny=input('Kastenlänge in y:');
Ny=21;
%iTs=input('Türposition in x:');
iTs=11;
%jTs=input('Türposition in y:');
jTs=1;
S=Sfeld(Nx,Ny,iTs,jTs);
D=zeros(Nx,Ny);
Bosons=zeros(Nx,Ny);
xi=zeros(Nx,Ny);
for i=2:Nx-1
    for j=2:Ny-1
        xi(i,j)=1;
    end;
end;
for k=1:size(iTs)
    xi(iTs(k),jTs(k))=1;
end;
%M=input('Anzahl der Teilchen:');
M=110;
%Tmax=input('maximale Zeitschritte:');
Tmax=11000;
%for v=0:10
%ks=input('ks:');
ks=1;
%kd=input('kd:');
kd=2.5;
%delta=input('delta:');
delta=0.3;
%alpha=input('alpha:');
alpha=0.3;
%Dmax=input('Dmax:');
Dmax=20;
%Dpic=input('Dfeld-Bilder zu zeitpunkten:');
Dpic=[5;50;100;150;200;250;300;350;400;450];
dpicture=zeros(Nx,Ny,1);

scrsz = get(0,'ScreenSize');
figure('Position',[scrsz(3)/6 scrsz(4)/6 2*scrsz(3)/3 2*scrsz(4)/3])
clear scrsz;
%for u=1:5
n=zeros(Nx,Ny);
k=0;
while k<M
    x=ceil(rand*(Nx-2))+1;
    y=ceil(rand*(Ny-2))+1;
    if n(x,y)~=1 && xi(x,y)~=0
        n(x,y)=1;
        k=k+1;
    end;
end;
%disp(sum(sum(n)));
t=0;
while t<Tmax
    t=t+1;
    D=dfeld(Nx,Ny,D,Bosons,xi,delta,alpha,Dmax);
    Bosons=zeros(Nx,Ny);
    sprung=zeros(Nx,Ny);
    for i=1:Nx
        for j=1:Ny
            if n(i,j)==1
                ppre(1)=exp(ks*S(i+1,j))*exp(kd*D(i+1,j))*(1-sign(n(i+1,j)))*xi(i+1,j);
                ppre(2)=exp(ks*S(i,j+1))*exp(kd*D(i,j+1))*(1-sign(n(i,j+1)))*xi(i,j+1);
                ppre(3)=exp(ks*S(i-1,j))*exp(kd*D(i-1,j))*(1-sign(n(i-1,j)))*xi(i-1,j);
                ppre(4)=exp(ks*S(i,j-1))*exp(kd*D(i,j-1))*(1-sign(n(i,j-1)))*xi(i,j-1);
                ppre(5)=exp(ks*S(i,j))*exp(kd*D(i,j))*xi(i,j);
                N=1/sum(ppre);
                p=[0 0 0 0 0];
                for k=1:5
                    p(k)=N*ppre(k);
                end;
                zuf=rand;
                if zuf<p(1)
                    sprung(i,j)=1;
                elseif zuf<p(2)+p(1)
                    sprung(i,j)=2;
                elseif zuf<p(3)+p(2)+p(1)
                    sprung(i,j)=3;
                elseif zuf<p(4)+p(3)+p(2)+p(1)
                    sprung(i,j)=4;
                end;
            end;
        end;
    end;
    for i=2:Nx-1
        for j=2:Ny-1
            k=0;
            b=[0,0,0,0];
            if sprung(i-1,j)==1
                k=k+1;
                b(1)=k;
            end;
            if sprung(i,j-1)==2
                k=k+1;
                b(2)=k;
            end;
            if sprung(i+1,j)==3
                k=k+1;
                b(3)=k;
            end;
            if sprung(i,j+1)==4
                k=k+1;
                b(4)=k;
            end;
            l=1;
            if k>1
                zuf=rand;
                if zuf<2/k && zuf>=1/k
                    l=2;
                elseif zuf<3/k && zuf>=2/k
                    l=3;
                elseif zuf>=3/k
                    l=4;
                end;
                for a=1:4
                    if b(a)~=l && b(a)~=0;
                        if a==1
                            sprung(i-1,j)=0;
                        end;
                        if a==2
                            sprung(i,j-1)=0;
                        end;
                        if a==3
                            sprung(i+1,j)=0;
                        end;
                        if a==4
                            sprung(i,j+1)=0;
                        end;
                    end;
                end;
            end;
        end;
    end;
    for i=1:Nx
        for j=1:Ny
            if sprung(i,j)~=0
                n(i,j)=0;
                Bosons(i,j)=1;
                if sprung(i,j)==1
                    n(i+1,j)=1;
                elseif sprung(i,j)==2
                    n(i,j+1)=1;
                elseif sprung(i,j)==3
                    n(i-1,j)=1;
                else
                    n(i,j-1)=1;
                end;
            end;
        end;
    end;
    clear sprung;
    plotmatrix=zeros(Nx+1,Ny+1);
    for i=1:Nx
       for j=1:Ny
          plotmatrix(i,j) = xi(i,j)+n(i,j);
          plotmatrix(Nx+1,j) = 1;
       end
       plotmatrix(i,Ny+1)=1;
    end
    pcolor(plotmatrix)    
    view(2)
    colormap gray
    axis xy
    axis equal
    axis off
    drawnow
    clear plotmatrix;
    for k=1:size(iTs)
        if n(iTs(k),jTs(k))==1
            n(iTs(k),jTs(k))=0;
            if max(max(n))==0
                disp(t);
                disp('evakuiert')
                for a=1:size(dpicture,3)
                    pause;
                    plotmatrix=zeros(Nx,Ny);
                    if a>=2
                    for i=1:Nx
                        for j=1:Ny
                            plotmatrix(i,j) = dpicture(i,j,a)+xi(i,j)-1;
                        end
                    end
                    h=max(max(max(dpicture)));
                    else
                    for i=1:Nx
                        for j=1:Ny
                            plotmatrix(i,j) = S(i,j);
                        end
                    end
                    h=max(max(S));
                    end;                    
                    surf(plotmatrix)
                    colormap gray
                    axis xy
                    axis([1,Nx+1,1,Ny+1,-1,h])
                    drawnow
                    clear plotmatrix;
                end;
                %clear;
                return;
            end;
        end;
    end;  
%    if max(max(n))==0
%               z(u,v+1)=t;
%               disp(t);
%               t=Tmax;
%               continue;
%    end;
    if find(Dpic==t)>0
        dpicture(:,:,size(dpicture,3)+1)=D;
    end;
    pause(0.02);
end;
disp(t);
disp(sum(sum(n)));
disp('noch im Raum');
for a=2:size(dpicture,3)
                    plotmatrix=zeros(Nx+1,Ny+1);
                    h=max(max(max(dpicture)));
                    for i=1:Nx
                        for j=1:Ny
                            plotmatrix(i,j) = dpicture(i,j,a)+xi(i,j)-1;
                        end
                    end
                    surf(plotmatrix)
                    colormap gray
                    axis xy
                    axis([1,Nx+1,1,Ny+1,-1,h])
                    drawnow
                    clear plotmatrix;
    pause; 
end;
%end;
%end;
%clear;