-- A function to calculate the Auger-Meitner selfenergy
-- Input:
--      phiX: list of X-Ray excited states before Auger-Meitner scattering
--      psi:  list of states after Auger-Meitner scattering
--      EAuger: list of energies of the psi
--      EDf, EDf: large and small part of radial wavefunctions of EDs; this is a list of InterpolatingFunction
--      Econt: list of energies corresponding to the above EDs
--      valence_orbitals: labels for active orbitals
--      kappa: relativistic quantum number of emittetd Auger-Meitner electron
--      wMin, wMax, dw: the selfenergy is evaluated at equidistant points between wMin and wMax in steps of dw
-- Output:
--      selfenergy: table of compact matrices representing the selfenergy evaluated at points specified by  wMin, wMax, dw
function Calculate_Selfenergy(phiX, psi,EAuger, EDg, EDf, Econt, valence_orbitals, kappa, wMin, wMax, dw,gammaSE)
    TimeStart("Calc Self-Energy")
    -- first we need to add indices for unbound 1-particle states:
    local nContIndices = 2*(abs(kappa))
    local m = {}
    for i = 1, nContIndices do
        m[i] = NF + i -1
    end
    NF = NF + nContIndices
    
    -- second we need to add an Auger-Meitner electron to the states psi
    -- this is done with the Tensor product
    contdet = "1"
    for i = 2, nContIndices do contdet = contdet .."0" end
    psiCont = NewWavefunction(nContIndices,0, {{contdet,1}})
    for i = 1, #psi do
        psi[i] = Simon.TensorProduct(psi[i],psiCont)
    end
    OppLadder = {}
    for i = 2, #m do
        OppLadder[i-1] = NewOperator("Number",NF,m[i],m[1])
    end
    psiAll = {psi}
    for i = 1, #OppLadder do
        psiAll[i+1] = {}
        for j = 1, #psi do
            psiAll[i+1][j] = OppLadder[i] * psi[j]
            EAuger[#EAuger+1] = EAuger[j]
        end
    end
    psi = Simon.Flatten(psiAll)
    
    -- we also need to modify the states phiX without Auger-Electrons such that they are in the same Hilbert space
    zeroes = "0"
    for i = 2, nContIndices do zeroes = zeroes .."0" end
    psiCont = NewWavefunction(nContIndices,0, {{zeroes,1}})
    for i = 1, #phiX do
        phiX[i] = Simon.TensorProduct(phiX[i],psiCont)
    end

    -- third we calculate the Coulomb-Operators without the Slaterintegrals:
    Uk = Martin.CreateAugerMeitnerU(NF,valence_orbitals,m,kappa)
    -- and then we calculate the Slater integrals for the Auger-Meitner scattering:
    FkNf = Martin.CreateContinuumSlaterIntegrals(G_vale,F_vale,valence_orbitals.kappas,EDg,EDf,kappa)
    -- note that calculating Uk and Fk separately saves a lot of computations now and later on (see PhD thesis of M Brass for further explanation)
    -- for these later purposes we need to flatten this list:
    Fk = {}
    for i = 1, #FkNf do
        Fk[i] = Simon.Flatten(FkNf[i])
    end
    -- we transform the Hartree units to eV:
    for i = 1, #Fk do
        for j = 1, #Fk[i] do
            Fk[i][j] = Fk[i][j] * EnergyUnits.Ha.value
        end
    end
    
    -- finally we calculate the selfenergy:
    selfenergy = Martin.CNSE(psi,phiX,Uk,wMin,wMax,dw,gammaSE,Fk,Econt,EAuger)
    TimeEnd("Calc Self-Energy")
    return selfenergy
end
                             

                        































