-- Create Eigendifferentials (EDs) with energies from Emin up to Emin + Ned * dE 
-- units for input and output are Hartree atomic units
-- Input:
--      Emin = minimal Energy for EDs
--      dE = energy difference between two EDs
--      Ned = Number of EDs
--      relativistic quantum number kappa
--      grid for radial wavefunctions
-- Output:
--      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 in Hartree atomic units

function CreateEigenDifferentials(Emin,dE,Ned,kappa,grid)
    Econt = {}
    for i = 1, Ned do
        Econt[i] = (i + 1/2)*dE
    end
    EDg , EDf = Martin.CreateEigenDifferentials(Emin,dE,Ned,kappa,grid,1)
    return EDg, EDf, Econt
end

-- Create Eigendifferentials (EDs) with energies from Emin up to Emin + Ned * dE which are orthogonal to bound states G, F
-- Input:
--      same as above
--      G, F large and small parts of bound state radial wavefunctions; this is a list of InterpolatingFunction
function CreateOrthogonalEigenDifferentials(Emin,dE,Ned,kappa,indices,grid,G,F)
    local basisfuns = {}
    local derivfuns = {}
    for i = 1, # G do
        Gval = G[i](grid) 
        Fval = F[i](grid)
        basisfuns[i] = Join(Gval,Fval)
        Gval = G[i](grid,1) 
        Fval = F[i](grid,1)
        derivfuns[i] = Join(Gval,Fval)
    end
    local positions_bound_states_same_kappa = {}
    for i = 1, #indices do 
        if kappa == indices.kappas[i] then
            positions_bound_states_same_kappa[#positions_bound_states_same_kappa+1] = i
        end
    end
    Econt, BaseContinuum = Martin.ContinuumBase(basisfuns,derivfuns,grid,{kappa,positions_bound_states_same_kappa},Emin,dE,Ned) 
    EDg, EDf = InterpolateEDList(BaseContinuum,grid)
    return EDg, EDf, Econt
end

-- just some utility for the above function, to interpolate the EDs on their grid:
function InterpolateED(ed, x)
    local gval = {}
    local fval = {}
    for i = 1, #x do
        gval[i] = ed[i]
        fval[i] = ed[i+#x]
    end
    local gint = InterpolatingFunction.Spline(x,gval)
    local fint = InterpolatingFunction.Spline(x,fval)
    return gint, fint
end

function InterpolateEDList(eds, x)
    local glist = {}
    local flist = {}
    for i = 1, #eds do
        glist[i], flist[i] = InterpolateED(eds[i],x)
    end
    return glist, flist
end


