dofile("../definitions.Quanty")

-- we define an Hamiltonian
Hamiltonian = OppU + 0.000001 * (2*OppSz+OppLz)

-- number of states calculated
Npsi=15
StartRestrictions = {Nf, Nb,  {"111111",2,2}}
psiList = Eigensystem(Hamiltonian, StartRestrictions, Npsi)

oppList={Hamiltonian, OppSsqr, OppLsqr, OppJsqr, OppSz, OppLz}

print(" <E>     <S^2>   <L^2>   <J^2>   <S_z>   <L_z>");
for i=1,#psiList do
  for j=1,#oppList do
    io.write(string.format("%7.4f ",psiList[i]*oppList[j]*psiList[i]))
  end
  io.write("\n")
end

Hamiltonian = Hamiltonian + 0.1 * Oppldots

oppList={Hamiltonian, OppSsqr, OppLsqr, OppJsqr, OppSz, OppLz}


-- Recalculate eigen states with spin orbit
Eigensystem(Hamiltonian,psiList)

print(" <E>     <S^2>   <L^2>   <J^2>   <S_z>   <L_z>");
for i=1,#psiList do
  for j=1,#oppList do
    io.write(string.format("%7.4f ",psiList[i]*oppList[j]*psiList[i]))
  end
  io.write("\n")
end