global hadamard cnot cphase sx sy sz rx ry rz N zz xx yy

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%pauli matrices

sx=[0 1; 1 0];
sy=[0 -i; i 0];
sz=[1 0; 0 -1];
si=[1 0; 0 1];

pauli={sx,sy,sz};
h=sym('h');
J=sym('J');
K=sym('K');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%two-qubit interaction terms: product space

zz=kron(sz,sz);
xx=kron(sx,sx);
yy=kron(sy,sy);
zi=kron(sz,si);
iz=kron(si,sz);
ii=kron(si,si);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%State Vectors--This will create a column vector of length 16 representing
%The 16 possible combinations of up and down

psi0=[1; 0];
psi1=[0; 1];

psi00=kron(psi0,psi0);
psi01=kron(psi0,psi1);
psi10=kron(psi1,psi0);
psi11=kron(psi1,psi1);

psi0000=kron(psi00,psi00);
psi0001=kron(psi00,psi01);
psi0010=kron(psi00,psi10);
psi0011=kron(psi00,psi11);
psi0100=kron(psi01,psi00);
psi0101=kron(psi01,psi01);
psi0110=kron(psi01,psi10);
psi0111=kron(psi01,psi11);
psi1000=kron(psi10,psi00);
psi1001=kron(psi10,psi01);
psi1010=kron(psi10,psi10);
psi1011=kron(psi10,psi11);
psi1100=kron(psi11,psi00);
psi1101=kron(psi11,psi01);
psi1110=kron(psi11,psi10);
psi1111=kron(psi11,psi11);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Basis State Vectors-This will create the state vectors symbolically
%u stands for spin up; d stands for spin down

uuuu=sym('uuuu');uuud=sym('uuud');uudu=sym('uudu');uudd=sym('uudd');
uduu=sym('uduu');udud=sym('udud');uddu=sym('uddu');uddd=sym('uddd');
duuu=sym('duuu');duud=sym('duud');dudu=sym('dudu');dudd=sym('dudd');
dduu=sym('dduu');ddud=sym('ddud');dddu=sym('dddu');dddd=sym('dddd');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Here is the basis of spin up down combinations

basis=[uuuu;uuud;uudu;uudd;uduu;udud;uddu;uddd;duuu;duud;dudu;dudd;...
        dduu;ddud;dddu;dddd];

%Test to see if it works for two spin 1/2 particles
%I believe it does

s1.s2=kron(zz,ii);
s2.s3=kron(kron(si,sz),kron(sz,si));
s3.s4=kron(ii,zz);
s4.s1=kron(kron(sz,si),kron(si,sz));

s1.s3=kron(kron(sz,si),kron(sz,si));
s2.s4=kron(kron(si,sz),kron(si,sz));

H=(J*h^2/4)*(s1.s2+s2.s3+s3.s4+s4.s1)+(K*h^2/4)*(s1.s3+s2.s4);


[evec,eval]=eig(H);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Print out eigenvectors and eigenvalues in nice fashion

for i=1:16
    a(i,1)=eval(i,i);
    a(i,2)=evec(:,i)'*basis;
end
    
disp('                     Eigenvalue                       Eigenvector');
disp(a);