function matrix = melt_hex(i,j,A,c)

[ni,nj] = getNeighbors(i,j);
ne = 1;
n = 2;
nw = 3;
sw = 4;
s = 5;
se = 6;

%Get center positions of neighbors in these directions
NE = c(ni(ne),nj(ne));
N = c(ni(n),nj(n));
NW = c(ni(nw),nj(nw));
SW = c(ni(sw),nj(sw));
S = c(ni(s),nj(s));
SE = c(ni(se),nj(se));

%The profuse conditions below are to avoid, as much 
%as possible, creating inconsitencies in symmetry
%during the melting process. If there used to be
%a neighboring cell in any of the six directions,
%we do not remove the points connecting this cell
%to the neighboring one, so that, if the neighboring
%cell is still frozen later on, it will not appear as
%half a cell, for instance. 

if (mod(j,2) == 0)
    A(i,j) = 0;
    %bottom row
    if (SW ~= 0 && S ~= 0)
        A(i-1,j-1) = 0;
    end
    if (S ~= 0 && SE ~= 0)
        A(i,j-1) = 0;
    end
    %middle row
    if (NW ~= 0 && SW ~= 0)
        A(i-1,j) = 0;
    end
    if (NE ~= 0 && SE ~= 0)
        A(i+1,j) = 0;
    end
    %top row
    if (N ~= 0 && NW ~= 0)
        A(i-1,j+1) = 0;
    end
    if (NE ~= 0 && NW ~= 0)
        A(i,j+1) = 0;
    end
else
    A(i,j) = 0;
    
    %bottom row
    if (SW ~= 0 && S ~= 0)
        A(i,j-1) = 0;
    end
    if (S ~= 0 && SE ~= 0)
        A(i+1,j-1) = 0;
    end
    
    %middle row
    if (NW ~= 0 && SW ~= 0)
        A(i-1,j) = 0;
    end
    if (NE ~= 0 && SE ~= 0)
        A(i+1,j) = 0;
    end

    %top row
    if (N ~= 0 && NW ~= 0)
        A(i,j+1) = 0;
    end
    if (NE ~= 0 && NW ~= 0)
        A(i+1,j+1) = 0;
    end
end

matrix = A;