N = 10^2;

m = 20; T = 10;
results = zeros(N,1);
X = ones(m,m);

for steps = 1:N    
    %Make bonds
    horizontal_bonds = zeros(m,m);
    vertical_bonds = zeros(m,m);
    for i = 1:m
        for j = 1:m
            %Horizontal bond
            U = exp( (1 + X(i,j) * X(i, mod(j,m) + 1)) / T) * rand;
            horizontal_bonds(i,j) = (U > 1);
            %Vertical bond
            U = exp( (1 + X(i,j) * X(mod(i,m) + 1, j)) / T) * rand;
            vertical_bonds(i,j) = (U > 1);
        end
    end
    
    %Find clusters
    S = 1:m^2;
    while isempty(S) == false
    current_state = min(S);
    C = [current_state];
    states_to_check = [];
    
    i = state_i(current_state,m);
    j = state_j(current_state,m);
    if horizontal_bonds(i,j) == 1
        add_state = ij_state(i, mod(j,m) +1,m);
        states_to_check = [states_to_check add_state];
    end
    if vertical_bonds(i,j) == 1
        add_state = ij_state(mod(i,m) + 1, j,m);
        states_to_check = [states_to_check add_state];
    end
    
    i_up = mod(i-2,m) + 1;
    j_left = mod(j-2,m) + 1;
    if horizontal_bonds(i,j_left) == 1
        add_state = ij_state(i, j_left,m);
        states_to_check = [states_to_check add_state];
    end
    if vertical_bonds(i_up,j) == 1
        add_state = ij_state(i_up, j,m);
        states_to_check = [states_to_check add_state];
    end
    
    while isempty(states_to_check) == false        
        current_state = min(states_to_check);
        C = [C current_state];
        states_to_check = setdiff(states_to_check, current_state);
        
        i = state_i(current_state,m);
        j = state_j(current_state,m);
        if horizontal_bonds(i,j) == 1
            add_state = ij_state(i, mod(j,m) +1,m);
            if (ismember(add_state,C) == false) && (ismember(add_state,states_to_check) == false)
                states_to_check = [states_to_check add_state];
            end
        end
        if vertical_bonds(i,j) == 1
            add_state = ij_state(mod(i,m) + 1, j,m);
            if (ismember(add_state,C) == false) && (ismember(add_state,states_to_check) == false)
                states_to_check = [states_to_check add_state];
            end
        end
        
        i_up = mod(i-2,m) + 1;
        j_left = mod(j-2,m) + 1;
        if horizontal_bonds(i,j_left) == 1
            add_state = ij_state(i, j_left,m);
            if (ismember(add_state,C) == false) && (ismember(add_state,states_to_check) == false)
                states_to_check = [states_to_check add_state];
            end
        end
        if vertical_bonds(i_up,j) == 1
            add_state = ij_state(i_up, j,m);
            if (ismember(add_state,C) == false) && (ismember(add_state,states_to_check) == false)
                states_to_check = [states_to_check add_state];
            end
        end        
    end    
    S = setdiff(S,S(ismember(S,C)));
     %Flip states
    if rand < 1/2
        for c_step = 1:length(C)
            i = state_i(C(c_step),m);
            j = state_j(C(c_step),m);
            X(i,j) = -1*X(i,j);
        end
    end
    end
    
   results(steps) = sum(sum(X))/m^2;
end

running_mean = cumsum(results')./(1:length(results));
HeatMap(X,'colormap', 'autumn')
plot(running_mean)