function [data100] = potts_100

%%%%%%%%%%%%%%%%%%% parameters %%%%%%%%%%%%%%%%%%%%%
q = 3;  
J = 1;  %%coupling
L = 100; %%size

%These values specify data taking for each temperature

N = 10;  %%number of bins of data for each temperature
Ntherm = 2000; %%number of lattice sweeps before taking data
Nfreq = 20; %%number of lattice sweeps before each 'measurement'
Nbin = 5; %%number of measurements in a bin



%matrix of data to be output
data100=[];


%%%%%%%%% sweep through temperatures around critical Temp %%%%%%%%%
f = 0;
for T = 0.95:0.002:1.05

%%initial state 
grid = randint(L,L,q)+1;


%%thermalize spin matrix
    
    grid = lattice_pass(grid,Ntherm,T);
    
%%grand sums for various quantities set to zero
Esum = 0;
E2sum = 0;
E4sum = 0;
Msum = 0;
M2sum = 0;
M4sum = 0;
Esigmabinsquaresum = 0;
Msigmabinsquaresum = 0;
EMsum = 0;
EM2sum = 0;
UMsum = 0;
UEsum = 0;
LgMsum = 0;
LgM2sum = 0;
Heatsum = 0;
Sussum = 0;

%%sweep through all of the bins
for t=1:N
    
    %%sums for each bin set to zero
    Ebin = 0;
    E2bin = 0;
    E4bin = 0;
    Mbin = 0;
    M2bin = 0;
    M4bin = 0;
    EMbin = 0;
    EM2bin = 0;
    
    %%collect measurements for each bin
    for g=1:Nbin
        
        
        %%Nfreq whole lattice passes between each measurement
        grid = lattice_pass(grid,Nfreq,T);
       
        %%make measurements of energy and magnetization of system, divide
        %%by L^2 to make values per spin, and increment bin sums
        Esites = ( (grid==circshift(grid, [ 0 1])) + (grid==circshift(grid, [ 0 -1])) + (grid==circshift(grid, [ 1 0])) + (grid==circshift(grid, [-1 0])) );
        E = sum(sum(J.*Esites))/(2*L^2);
       
        for j=1:q
            B(j)=length(find(grid==j));
        end
            M = (max(B)-(L^2-max(B))/(q-1))/L^2;
                    
        Ebin = Ebin + E;
        E2bin = E2bin + E^2;
        E4bin = E4bin + E^4;
        Mbin = Mbin + M;
        M2bin = M2bin + M^2;
        M4bin = M4bin + M^4;
        EMbin = EMbin + E*M;
        EM2bin = EM2bin + E*M^2;
          
    end
   
    %%normalize averages by number of measurements per
    %%bin, calculate bin error
    
    Ebin = Ebin/Nbin;
    E2bin = E2bin/Nbin;
    E4bin = E4bin/Nbin;
    Esigmabinsquare = (E2bin - Ebin^2)/Nbin;
    
    Mbin = Mbin/Nbin;
    M2bin = M2bin/Nbin;
    M4bin = M4bin/Nbin;
    Msigmabinsquare = (M2bin - Mbin^2)/Nbin;  
    
    EMbin = EMbin/Nbin;
    EM2bin = EM2bin/Nbin;
    
    %%increment grand sums
    Esigmabinsquaresum = Esigmabinsquaresum + Esigmabinsquare;
    Esum = Esum + Ebin;
    E2sum = E2sum + E2bin;
    E4sum = E4sum + E4bin;
    
    Msigmabinsquaresum = Msigmabinsquaresum + Msigmabinsquare;
    Msum = Msum + Mbin;
    M2sum = M2sum + M2bin;
    M4sum = M4sum + M4bin;
    
    EMsum = EMsum + EMbin;
    EM2sum = EM2sum + EM2bin;
    
    UMsum = UMsum + (1 - M4bin/(3*M2bin^2));
    UEsum = UEsum + (1 - E4bin/(3*E2bin^2));
    LgMsum = LgMsum + (EMbin/Mbin - Ebin);
    LgM2sum = LgM2sum + (EM2bin/M2bin - Ebin);
    
    Heatsum = Heatsum + (E2bin - Ebin^2);
    Sussum = Sussum + (M2bin - Mbin^2);
    
end

%%normalize grand averages by number of bins, calculate average error of
%%all bins and grand error
Esum = Esum/N;
E2sum = E2sum/N;
E4sum = E4sum/N;
Esigmasquare = (E2sum - Esum^2)/N;

Eerrorquad = sqrt(Esigmabinsquaresum/N);
Eerror = sqrt(Esigmasquare);


Msum = Msum/N;
M2sum = M2sum/N;
M4sum = M4sum/N;
Msigmasquare = (M2sum - Msum^2)/N;

Merrorquad = sqrt(Msigmabinsquaresum/N);
Merror = sqrt(Msigmasquare);

EMsum = EMsum/N;
EM2sum = EM2sum/N;
UMsum = UMsum/N;
UEsum = UEsum/N;
LgMsum = LgMsum/N;
LgM2sum = LgM2sum/N;

UMgrand = (1 - M4sum/(3*M2sum^2));
UEgrand = (1 - E4sum/(3*E2sum^2));
LgMgrand = (EMsum/Msum - Esum);
LgM2grand = (EM2sum/M2sum - Esum);

Heatsum = Heatsum/N;
Sussum = Sussum/N;

Heatgrand = E2sum - Esum^2;
Susgrand = M2sum - Msum^2;


disp('Esum E2sum E4sum Eerror Eerrorquad Msum M2sum M4sum Merror Merrorquad UMsum UMgrand UEsum UEgrand LgMsum LgMgrand LgM2sum LgM2grand Heatsum Heatgrand Sussum Susgrand');
disp([Esum E2sum E4sum Eerror Eerrorquad Msum M2sum M4sum Merror Merrorquad UMsum UMgrand UEsum UEgrand LgMsum LgMgrand LgM2sum LgM2grand Heatsum Heatgrand Sussum Susgrand]);

f=f+1;

%%save data for current temperature, increment temp, and repeat
data100(f,:)=[Esum E2sum E4sum Eerror Eerrorquad Msum M2sum M4sum Merror Merrorquad UMsum UMgrand UEsum UEgrand LgMsum LgMgrand LgM2sum LgM2grand Heatsum Heatgrand Sussum Susgrand];
end

save data100 data100




%%%%% subroutine to do Npass 'passes through lattice' %%%%%%%%%%%%%%%
function [grid] = lattice_pass(grid,Npass,T)
 
for tt=1:Npass  
    
%%Choose random possibility for each spin in lattice to flip to 
%%(excluding flipping to itself), calculate local change in energy 
%%for each spin if it does flip holding spins around it constant, 
%%calulate bolztman weight by change in energy, and perform metropolis
%%decision
rr = randint(L,L,q-1)+1;
grid2 = (rr<grid).*rr + ((rr==grid) + (rr>grid)).*(rr+1);
E1 = (grid==circshift(grid, [ 0 1])) + (grid==circshift(grid, [ 0 -1])) + (grid==circshift(grid, [ 1 0])) + (grid==circshift(grid, [-1 0]));
E2 = (grid2==circshift(grid, [ 0 1])) + (grid2==circshift(grid, [ 0 -1])) + (grid2==circshift(grid, [ 1 0])) + (grid2==circshift(grid, [-1 0]));
DeltaE = -J * (E2 - E1);

%%metropolis decision
p_trans = exp(-DeltaE/(T));

%%execute sping flips according to metropolis decision
R = rand(L);
transitions = (R < p_trans ) .* grid2 + (R > p_trans ) .* grid;

%%hold sublattice 1 constant,ie, set it back to original values
C = zeros(L,1);
C(1:2:length(C))=1;
toe = toeplitz(C,C);
pos = find(toe==1);
transitions(pos)=grid(pos);

%%spin lattice after one sublattice has flipped according to metropolis
%%decision
grid = transitions;



%%repeat process holding sublattice 2 fixed 
rr = randint(L,L,q-1)+1;
grid2 = (rr<grid).*rr + ((rr==grid) + (rr>grid)).*(rr+1);
E1 = (grid==circshift(grid, [ 0 1])) + (grid==circshift(grid, [ 0 -1])) + (grid==circshift(grid, [ 1 0])) + (grid==circshift(grid, [-1 0]));
E2 = (grid2==circshift(grid, [ 0 1])) + (grid2==circshift(grid, [ 0 -1])) + (grid2==circshift(grid, [ 1 0])) + (grid2==circshift(grid, [-1 0]));
DeltaE = -J * (E2 - E1);

p_trans = exp(-DeltaE/(T));

R = rand(L);
transitions = (R < p_trans ) .* grid2 + (R > p_trans ) .* grid;

C = zeros(L,1);
C(2:2:length(C)-1)=1;
toe = toeplitz(C,C);
pos = find(toe==1);
transitions(pos)=grid(pos);

%%spin lattice after both sublattices have flipped according to metropolis
%%decision
grid = transitions;
end
end
end