-- import Utility:
dofile("Utility_Functions.Quanty")
dofile("Utility_Operators.Quanty")

-- define the labels of atomic orbitals considered in the calculation:
orbital_labels = {"1s1/2","2s1/2","3s1/2","3p1/2","3p3/2","2p1/2","2p3/2","3d3/2","3d5/2"}

-- transform labels into indices for later usage in built-in functions:
indices, NF = CreateAtomicIndicesList(orbital_labels,{{"Kappas",true}}) 
-- ind contains lists of indices corresponding to the above labels, this is needed as input for other functions
-- ind.kappas contains the relativistic quantum numbers for the above orbitals: kappa = -l-1 if j = l+1/2 and kappa = l if j = l-1/2
-- NF = Number of fermionic 1-particle states
-- also have a look at CreateAtomicIndicesDict(orbital_labels)
NB = 0 -- no Bosons in calculation

-- nuclear charge of Atom:
Z = 22
Valence = 4 -- Ti is 4+
-- import 1-particle radial wave-functions from FPLO stored in +fval.*
--dir = "/scratch/brass/Resource/Ti4"
dir = "."
G, F, grid = ReadFPLOBasisFunctions(orbital_labels,dir.."/+fval.001.1", dir.."/+fval.001.2")    
-- G and F are radial wavefunctions of large and small part of a Dirac spinor that solves the relativistic Kohn-Sham equations
-- they are stored as list of InterpolatingFunction ordered according to orbitalLabels

-- define core and valence orbitals:
core_orbitals = Take(indices,1,2) 
core_orbitals.kappas = Take(indices.kappas,1,2) 
valence_orbitals = Take(indices,3,9) 
valence_orbitals.kappas = Take(indices.kappas,3,9) 
G_core = Take(G,1,2)
G_vale = Take(G,3,9) 
F_core = Take(F,1,2)
F_vale = Take(F,3,9) 

-- now we calculate all bound states before and after Auger-Meitner scattering
dofile("SrTiO_Calculate_Bound_States.Quanty")
-- this file created:
--      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
--      poles: list of poles of XAS Green's function
--	residues: list of weights corresponding to above poles

-- now we need to add the Auger-Meitner electrons to the above states:
-- the first argument supplied to Quanty defines kappa of unbound states
kappa_unbound = tonumber(args[1])
-- note that for a full calculation one has to loop over all possible values for kappa
-- this can be done via a batch script to paralleize this calculation on many computation nodes

-- we now create a single-particle basis for the unbound states using functions from:
print("------------------------------ Start Eigendifferentials ------------------------------")
dofile("Eigendifferentials.Quanty")
TimeStart("Eigendifferentials")
-- we define the energy discretization for the Eigendifferentials:
eV = 1/EnergyUnits.Ha.value
Emin = 0 -- smalest unbound energy
dE   = 1 * eV  -- energy spacing of 1eV in Hartree atomic units
dE = 10 *eV -- TODO: testing only
Ned  = 2000 -- number of Eigendifferentials
Ned = 50 -- TODO: testing only
gammaSE = 1 --3 -- small imaginary part in denominator of self-energy; 3*dE in eV is always a good value to get smooth spectra


if IsElement(kappa_unbound,indices.kappas) then
    EDg, EDf, Econt = CreateOrthogonalEigenDifferentials(Emin,dE,Ned,kappa_unbound,indices,grid,G,F)
else
    EDg, EDf, Econt = CreateEigenDifferentials(Emin,dE,Ned,kappa_unbound,grid)
end
-- we transform the Hartree units to eV:
for i = 1, #Econt do
    Econt[i] = Econt[i] * EnergyUnits.Ha.value
end
TimeEnd("Eigendifferentials")
print("------------------------------   End Eigendifferentials ------------------------------")

wMin = 460-- -10
wMax = 490 --20
dw = 1E-2--2.5E-2

print("------------------------------ Start selfenergy ------------------------------")
dofile("SrTiO_Calculate_Auger_Meitner_Selfenergy.Quanty")
selfenergy = Calculate_Selfenergy(phiX, psi,EAuger, EDg, EDf, Econt, valence_orbitals, kappa_unbound, wMin, wMax, dw,gammaSE)
print("------------------------------   End selfenergy ------------------------------")

-- now save output:
--Martin.WriteData("CompactMatrixType_Table","SrTiO_SE_V4_Dq26" ,selfenergy)

-- or plot the specrtum:
TimeStart("Plot")
Nw = floor((wMax-wMin)/dw) + 1
spec = {}
for i = 1, Nw do -- #cnse
    w = wMin + (i-1)*dw
    M = Matrix.ToTable(selfenergy[i]) 
    M0= Matrix.Identity(#M) -- for comparrison: spectrum without self-energy
    for j = 1, #M do
        M[j][j] = M[j][j] - w + poles[j] 
        M0[j][j]= poles[j]- w -I/115 -- small imaginary part for comparrison
    end
    R = Matrix.Inverse(M)
    R0 = Matrix.Inverse(M0)
    spec[i] = {w,0,0}
    for m = 1, #R do
        for j = 1, #R do
            spec[i][2] = spec[i][2] + Complex.Im( R[m][j] * Conjugate(residues[m][1]) * residues[j][1])
            spec[i][3] = spec[i][3] + Complex.Im(R0[m][j] * Conjugate(residues[m][1]) * residues[j][1])
        end
    end
end

TimeEnd("Plot")

file = io.open("SrTiO_XAS_CF_SE_V4_Dq26.dat","w")
io.output(file)

for i = 1, #spec do
   io.write(spec[i][1],"\t",spec[i][2],"\t",spec[i][3],"\n")
end
io.close(file)


TimePrint()































