a = 1;
Nx = 150;
Ny = 150/2;
Nyy = Ny*2;

A = zeros(Nx,Nyy);
centers = A; %when calling set_hex(i,j,A), must also set centers(i,j) = 1

%Get state sets
state1 = perms([0,0,0,0,0,1]);
state2 = perms([0,0,0,0,1,1]);
state3 = perms([0,0,0,1,0,1]);
state4 = perms([0,0,0,1,1,1]);
state5 = perms([0,0,1,0,0,1]);
state6 = perms([0,0,1,0,1,1]);
state7 = perms([0,0,1,1,1,1]);
state8 = perms([0,1,0,1,0,1]);
state9 = perms([0,1,0,1,1,1]);
state10 = perms([0,1,1,0,1,1]);
state11 = perms([0,1,1,1,1,1]);
state12 = perms([1,1,1,1,1,1]);

%snowflake1
% p_freeze = [0.2, 0, 0.6, 0, 0.3, 0.15, 0, 0.2, 0, 0.2, 0, 0.8];
% p_melt = [0, 0.7, 0.5, 0.7, 0.7, 0.5, 0.3, 0.5, 0.3, 0.2, 1, 0];

%snowflake2 & better_snowflake4
% p_freeze = [0.2, 0.4, 0.6, 0.2, 0.3, 0.15, 0.4, 0.5, 0.34, 0.2, 0.3, 0.3];
% p_melt = [0.1, 0.1, 0.2, 0.3, 0.4, 0.5, 0.7, 0.5, 0.7, 0.6, 0.9, 1];

%snowflake3 & better_snowflake3
% p_freeze = [0.2, 0.4, 0.6, 0.2, 0.3, 0.15, 0.4, 0.5, 0.34, 0.2, 0.3, 0.3];
% p_melt = [1, 0.9, 0.87, 0.78, 0.7, 0.8, 0.88, 0.78, 0.5, 0.3, 0.2, 0.1];

%snowflake4 & better_snowflake1
% p_freeze = [0.2, 0.5, 0.53, 0.43, 0.33, 0.25, 0.3, 0.45, 0.34, 0.2, 0.3, 0.3];
% p_melt = [1, 0.9, 0.87, 0.78, 0.7, 0.8, 0.88, 0.78, 0.5, 0.3, 0.2, 0.1];

%better_snowflake2
% p_freeze = [0.1, 0.3, 0.43, 0.45, 0.43, 0.55, 0.5, 0.75, 0.84, 0.8, 0.7, 0.63];
% p_melt = [1, 0.9, 0.87, 0.78, 0.7, 0.8, 0.88, 0.78, 0.5, 0.3, 0.2, 0.1];

%better_snowflake5: harder to freeze
% p_freeze = [0.2, 0.3, 0.43, 0.51, 0.53, 0.55, 0.6, 0.75, 0.84, 0.8, 0.8, 0.83];
% p_melt = [1, 0.9, 0.87, 0.78, 0.7, 0.8, 0.88, 0.78, 0.5, 0.3, 0.2, 0.1];

%better_snowflake6: even harder to freeze
% p_freeze = [0.3, 0.35, 0.5, 0.51, 0.53, 0.6, 0.65, 0.75, 0.84, 0.8, 0.9, 0.93];
% p_melt = [0.1, 0.9, 0.87, 0.78, 0.7, 0.8, 0.68, 0.58, 0.5, 0.3, 0.2, 0.1];

%better_snowflake7: random
p_freeze = [0.3, 0.6, 0.3, 0.4, 0.8, 0.65, 0.35, 0.75, 0.74, 0.2, 0.78, 0.93];
p_melt = [0.14, 0.39, 0.67, 0.78, 0.54, 0.8, 0.69, 0.23, 0.5, 0.56, 0.8, 0.3];

A = set_hex(Nx/2,Nyy/2,A);
centers(Nx/2,Nyy/2) = 1;
newA = A;
oldCenters = centers;
numIter = 36;

for i = 1:numIter
%     i
    rand_set = [rand,rand,rand,rand,rand,rand,rand,rand,rand,rand,rand,rand];%generate 12 random numbers per iteration
    for x = 5:Nx-5
        for y = 5:Nyy-5
            [nextdoor_x,nextdoor_y] = getNeighbors(x,y);
            
            %get set state of this ice cell
            state_formation = arrayfun(@(i,j) oldCenters(i,j), nextdoor_x, nextdoor_y);
            
            if (all(state_formation == 0)) %can't grow something from nothing
                continue
            end
            
            state_num = 0;
            if (ismember(state_formation,state1,'rows'))
                state_num = 1;
            elseif (ismember(state_formation,state2,'rows'))
                state_num = 2;
            elseif (ismember(state_formation,state3,'rows'))
                state_num = 3;
            elseif (ismember(state_formation,state4,'rows'))
                state_num = 4;
            elseif (ismember(state_formation,state5,'rows'))
                state_num = 5;
            elseif (ismember(state_formation,state6,'rows'))
                state_num = 6;
            elseif (ismember(state_formation,state7,'rows'))
                state_num = 7;
            elseif (ismember(state_formation,state8,'rows'))
                state_num = 8;
            elseif (ismember(state_formation,state9,'rows'))
                state_num = 9;
            elseif (ismember(state_formation,state10,'rows'))
                state_num = 10;
            elseif (ismember(state_formation,state11,'rows'))
                state_num = 11;
            elseif (ismember(state_formation,state12,'rows'))
                state_num = 12;
            end
            %end get set state
            
            if (oldCenters(x,y) == 1) %previously frozen
                if (rand_set(state_num) > p_melt(state_num))
                    centers(x,y) = 0;
                    newA = melt_hex(x,y,newA,oldCenters);
                end
            else %previously no ice here
                if (rand_set(state_num) > p_freeze(state_num))
                    centers(x,y) = 1;
                    newA = set_hex(x,y,newA);
                end
            end
        end
    end
    plot_hexagonal_lattice(newA);
    whitebg(1,'k');
    set(gcf,'color','black');
    set(gcf, 'InvertHardCopy', 'off');
    axis off
    drawnow;
    saveas(gcf,['better_snowflake7_',sprintf('%05d',i),'.png']);
    oldCenters = centers;
end
