%Anthony Goodrow
%ChE 240 Thermo. 

clear A; clear m; clear t;

L=input('Enter lattice size L: ');
T=input('Enter dimensionless T: ');
B=input('Enter magnetic field B: ');
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;
trial=1;

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)+B);
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)+B);
    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)+B);
        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)+B);
            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)+B);
                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)+B);
                    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)+B);
                        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)+B);
                            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)+B);
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

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

end

%Determine mean magnetization per spin, <m>=<M>/N at the end of each complete time step
m(k,1)=sum(sum(A))/(L*L);
t(k,1)=t(k-1,1)+1;
k=k+1;

end

figure(1)
plot(t(1:length(t)),m(1:length(m),1))
title('Mean Magnetization Per Spin, <m>=<M>/N, vs. Time')
xlabel('Time')
ylabel('Mean Magnetization Per Spin, <m>')

r=50; %average over the last r values
average(trial,1)=sum(m(length(m)-r:length(m),1))/(r+1)