%A=ones(L,L)

%L=input('Enter lattice size L: ');
L=10;
T=input('Enter temperature T: '); %in units of J, i.e. kb*T

for i=1:L
    for j=1:L
        A(i,j)=-1+2*round(rand(1)); %+1=A, -1=B
    end
end

%get energy contacts
%go through A and find AA contacts
contact=0;
for i=1:L
    for j=1:L
        if A(i,j)==1
            %check all neighbors for contacts, periodic bc
            contacts
            contactAA=contact/2; %divide by 2 to correct for double counting
        end
    end
end

%go through A to find BB contacts
contact=0;
for i=1:L
    for j=1:L
        if A(i,j)==-1
            contacts
            contactBB=contact/2;
        end
    end
end

%go through A to find AB contacts
contact=0;
for i=1:L
    for j=1:L
        if A(i,j)==-1
            contactsAB
            contactAB=contact/2;
        end
    end
end

%get numbers you need
N=L*L;
z=4; %for 2D lattice
EAA=1;
EBB=1;
EAB=(EAA+EBB)/2;

phiA=(L^2/2+sum(sum(A))/2)/N;
phiB=1-phiA;

E=N*z*(1/2*EAA*phiA^2+1/2*EBB*phiB^2+EAB*phiA*phiB);
Omega=factorial(N)/(factorial(round(N*phiA))*factorial(round(N*phiB)));
kb=1; %set kb=1, arbitrary units

%results
Q=Omega*exp(-E/(kb*T)) 
AN=(-kb*T*(phiA*log(phiA)+phiB*log(phiB))+z*((EAA+EBB)/2-EAB)*phiA*phiB)/N
GNRT=phiA*log(phiA)+phiB*log(phiB)
phiA
phiB