-- This function takes a ligand field Hamiltonian and adds NShels ligand orbitals
-- this turns a Ligand field Hamiltonian into an Anderson Impurity Hamiltonian
-- Note that the number of many particle determinants explodes
-- You can tame this by setting restrictions. In order to do this we use bit masks
-- You can define a set of bitmasks for the ligand field model. The function extends
-- them for the additional bath sites.

function ExpandLigandBasis(NShels, NF, IndexDn_L, IndexUp_L, tMat, epsilon0, BitMasks)
  -- We add fermions to the one particle basis
  local NFReturn = NF + NShells * (#IndexDn_L + #IndexUp_L)

  -- We extend the bit masks to also include the new bath sites
  local BathBitMask=""
  for i=1,NF do
    BathBitMask = BathBitMask .. "0"
  end
  BathBitMask = BathBitMask .. " "
  for j=1,#BitMasks do
    BitMasks[j] = BitMasks[j] .. " "
  end
  for i=NF+1,NFReturn do
    BathBitMask = BathBitMask .. "1"
    for j=1,#BitMasks do
      BitMasks[j] = BitMasks[j] .. "0"
    end
  end
  BitMasks[#BitMasks+1] = BathBitMask

  -- We create the one particle Hamiltonian for the extended bath.
  -- The dimension of IndexDn_L and IndexUp_L need to be the same
  -- tMat and epsilon need to be square matrices of dimension IndexDn_L
  HBathMat = Matrix.Zero(NShells+1)
  HBathMat[1][1] = Matrix.Zero(NF)
  for i=1,NShells do
    HBathMat[i+1][i+1] = Matrix.AddSpin(epsilon0)
  end
  for i=1,NShells-1 do
    HBathMat[i+1][i+2] = Matrix.AddSpin(tMat)
    HBathMat[i+2][i+1] = Matrix.AddSpin(tMat)
  end
  HBathMat = Matrix.Flatten(HBathMat)
  -- Lastly we need to couple the extended bath to the previous system
  for i=1,#IndexDn_L do
    for j=1,#IndexDn_L do
      HBathMat[1+IndexDn_L[i]][NF+2*j-1] = tMat[i][j]
      HBathMat[NF+2*j-1][1+IndexDn_L[i]] = tMat[i][j]
      HBathMat[1+IndexUp_L[i]][NF+2*j] = tMat[i][j]
      HBathMat[NF+2*j][1+IndexUp_L[i]] = tMat[i][j]
    end
  end
  -- And we convert the matrix to an operator  
  HBath = Matrix.ToOperator(HBathMat)
  return NFReturn, HBath, BitMasks
end
