%
%    Eenvoudige Monte Carlo simulatie van het groei/oplos proces 
%                    van het Kosselkristal
%
%                     Willem van Enckevort
%
% --------------------------------------------------------------------
%
%
% Invoermogelijkheden: - Lengte en breedte 2D array (bijv. 40 x 20)
%                      - Starthoogte (bijv. 10)
%                      - Aantal treden (bijv. 0, 1, 2, 3)
%                      - Bindingssterkte phi/kT (van 0.2 tot 1.5)
%                      - Drijvende kracht Dmu/kT (van -1 via 0 (evenwicht)
%                        tot 3)
%                      - Aantal malen resultaat weergeven (NVieuws)
%                        (Druk willekeurige toets in na iedere view)
%                      - Aantal cycli per tussenresultaat (NCycl)
%
% --------------------------------------------------------------------
%
% Uitvoer: - Bovenaanzicht kristaloppervlak bij ieder tussenresultaat
%          - Interface dichtheids profiel
%          - Groei- of oplossnelheid
%
% --------------------------------------------------------------------
%
%                         Start programma
%
% Maak alles schoon
clear
%
% INVOER DATA:
%
Len = 40;      % Lengte Array
Width = 20;    % Breedte Array 
Height = 10;    % Starthoogte Array
%
NSteps = 0;    % Aantal treden 
%
PhikT = 1.0;   % Phi/kT
DMukT = 1.0;   % DMu/kT
%
% Aantal views
NVieuws = 4;
% Aantal cycli per view
NCycl = 50000;
%
% --------------------------------------------------------------------
%
% Voorbereiding: initialisatie parameters 
%
% Definieer de kansen
Pp = exp(-4*PhikT)*exp(DMukT);
for i = 1:5
    Pm(i)=exp(-2*(i-1)*PhikT);
end
% Normaliseer de kansen (Maximale kans: Pm(1)=1, d.w.z 0 buren)
PpN = Pp/(Pp+Pm(1));
%
% Startoppervlak zonder treden
Surface(1:Width,1:Len) = Height;
% Toevoeging treden
for i = 1 : Len
    for j = 1 : Width
        Surface(j,i) = Surface(j,i) + floor(NSteps*i/Len);
    end
end
%
% Geef het originele startoppervlak weer
figure(1)
colormap(Copper)  
imagesc(Surface)
colorbar
%
pause
%
% Initialiseer tellers
%
Nadd = 0;
Nrem = 0;
%
% -----------------------------------------------------------------
%
% Start de simulatie
%
for Vieuws = 1 : NVieuws
    for Cycle = 1:NCycl
        % Kies willekeurige positie
        X = randi(Len);  
        Y = randi(Width);
        % Groei of ets of niets
        Toss1 = rand(1);
        % Groei
        if Toss1 <= PpN
            Surface(Y,X) = Surface(Y,X) + 1;
            Nadd = Nadd + 1;
        else
        % Ets of niets
        %
        % Periodieke randvoorwaarden: Tel de buren
        %
        NBuur = 0;
        %
        % Linker buur
        if X == 1
            XL = Len;
            DZ = -NSteps;
        else
            XL = X-1;
            DZ = 0;
        end
        if Surface(Y,X) <= Surface(Y,XL)+DZ
            NBuur = NBuur + 1;
        end
        %
        % Rechter buur
        if X == Len
            XR = 1;
            DZ = NSteps;
        else
            XR = X+1;
            DZ = 0;
        end
        if Surface(Y,X) <= Surface(Y,XR)+DZ
            NBuur = NBuur + 1;
        end
        %
        % Bovenbuur
        if Y == Width
            Yup = 1;
        else
            Yup = Y+1;
        end
        if Surface(Y,X) <= Surface(Yup,X)
            NBuur = NBuur + 1;
        end
        %
        % Onderbuur
        if Y == 1
            Ylo = Width;
        else
            Ylo = Y-1;
        end
        if Surface(Y,X) <= Surface(Ylo,X)
            NBuur = NBuur + 1;
        end
        %
        % Etsen of niets?
            Toss2 = rand(1);
            if Toss2 <= Pm(NBuur+1)
                Surface(Y,X) = Surface(Y,X) - 1;
                Nrem = Nrem + 1;
            end
        end    
    end
%
% Geef het (tussen)resultaat weer
colormap(Copper)    
imagesc(Surface)
colorbar
pause
end
% 
% --------------------------------------------------------------------
% Geef de resultaten weer
%
% Print de groeisnelheid (sticking fraction)
    Stickingfraction = (Nadd-Nrem)/Nadd
%
% Bereken interface profile bij NSteps = 0
if NSteps == 0
    MinHeight = min(min(Surface));
    MinHeight = MinHeight - 3;
    MaxHeight = MinHeight + 15;
    HeightProfile (1:15)= 0;
    for H = 1:15
        Height(H)= H;
    end    
    for X = 1:Len
        for Y= 1:Width
            for H = 1:15
                if Surface(Y,X) == H + MinHeight
                    for i = 1:H
                        HeightProfile(i) = HeightProfile(i) + 1;
                    end    
                end
            end
        end
    end
    HeightProfile = HeightProfile/(Width*Len);
    %
    % Maak Plotje
    figure(2)
    plot(Height,HeightProfile,'bx-')
    TITLE('Interface density profile');
    XLABEL('Height');
    YLABEL('Density');
end    
%             
% Einde programma