%Anthony Goodrow
%ChE 240 Thermo.

L=input('Enter lattice size L: ');
T=input('Enter dimensionless T: ');
tmax=input('Enter tmax: ');

A=zeros(L);

%Populate with +1 or -1
for i=1:L
    for j=1:L
        rnd=rand(1);
        if rnd<0.5
            A(i,j)=1;
        else
            A(i,j)=-1;
        end
    end
end

%Initial <m>
m(1,1)=sum(sum(A))/(L*L);
t(1,1)=0;
k=2;

for run=1:tmax
    
for time=1:L*L %one time step
%Pick random site, flip spin "artificially", and check probability for acceptance
%Use periodic boundary conditions

%Random location between 1 and L for x and y directions
i=ceil(L*rand(1));
j=ceil(L*rand(1));

%Flip
A(i,j)=-A(i,j);

%Sum spins, check for edges first, use flipped value to find energy change

if j~=1 & j~=L & i~=1 & i~=L %not on edge or corner
    S=2*A(i,j)*(A(i,j-1)+A(i,j+1)+A(i-1,j)+A(i+1,j));
else
    if j==1 & i~= 1 & i~=L %on left hand wall, but not in corner
        S=2*A(i,j)*(A(i,L)+A(i,j+1)+A(i-1,j)+A(i+1,j));
    else
        if j==L & i~=1 & i~=L %on right hand wall, but not in corner
            S=2*A(i,j)*(A(i,j-1)+A(i,1)+A(i-1,j)+A(i+1,j));
        else
            if i==1 & j~=1 & j~=L % on top wall, but not in corner
                S=2*A(i,j)*(A(i,j-1)+A(i,j+1)+A(L,j)+A(i+1,j));
            else
                if i==L & j~=1 & j~=L %on bottom wall, but not in corner
                    S=2*A(i,j)*(A(i,j-1)+A(i,j+1)+A(i-1,j)+A(1,j));
                else
                    if i==1 & j==1 %top left corner
                        S=2*A(i,j)*(A(i,L)+A(i,j+1)+A(L,j)+A(i+1,j));
                    else
                        if i==L & j==1 %bottom left corner
                            S=2*A(i,j)*(A(i,L)+A(i,j+1)+A(i-1,j)+A(1,j));
                        else
                            if i==1 & j==L %top right corner
                                S=2*A(i,j)*(A(i,j-1)+A(i,1)+A(L,j)+A(i+1,j));
                            else
                                if i==L & j==L %bottom right corner
                                    S=2*A(i,j)*(A(i,j-1)+A(i,1)+A(i-1,j)+A(1,j));
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

move=rand(1);
P=exp(S/T);
if P<move
    %don't make move, change back to original value
    A(i,j)=-A(i,j);
    %else accept move
end
    
end

    if mod(run,10)==0
       figure(1)
%     hold on
     imagesc(A)
        
    end

end