text:
    export compute

    """
    `Basics.compute("angular coefficients: e-e, Ratip2013", csfa::CsfR, csfb::CsfR)`
        ... to compute the angular coefficients in the decomposition of a (reduced) many-electron matrix element with a general
            rank-0 electron-electron interaction operator ``⟨csf_a ||V(e-e)|| csfb⟩ = \\sum_t  T(a_t, b_t, c_t, d_t) * R^k (a_t, b_t, c_t, d_t)``
            by a call to the Fortran procedure `anco_calculate_csf_pair` of the RATIP program; a Tuple{Array{AngularTcoeff,1},Array{AngularVcoeff,1}}`
            is returned.
    """
    function Basics.compute(sa::String, csfa::CsfR, csfb::CsfR)
        if sa == "angular coefficients: e-e, Ratip2013"
            if csfa == csfb
                JAC.AngularCoefficientsRatip2013.load_csl(csfa)
                t_coeffs, v_coeffs = JAC.AngularCoefficientsRatip2013.angular_coefficients_pair(1,1)
            else
                JAC.AngularCoefficientsRatip2013.load_csl(csfa, csfb)
                t_coeffs, v_coeffs = JAC.AngularCoefficientsRatip2013.angular_coefficients_pair(1,2)
            end
        else
            error("Unsupported keystring = ", sa)
        end
        return t_coeffs, v_coeffs
    end

    """
    `Basics.compute("angular coefficients: 1-p, Ratip2013", rank::Int64, csfa::CsfR, csfb::CsfR)`
        ... to compute the the angular coefficients in the decomposition of a (reduced) many-electron matrix element with a general
            single-particle operator of the given rank ``⟨csf_a ||O^rank|| csfb⟩ = \\sum_t  T(a_t, b_t) * R (a_t, b_t)``  by a call
            to the Fortran procedure `anco_calculate_csf_pair_1p` of the RATIP program; an `Array{AngularTcoeff,1}` is returned.
    """
    function Basics.compute(sa::String, rank, csfa::CsfR, csfb::CsfR)
        if      sa == "angular coefficients: 1-p, Ratip2013"
            if csfa == csfb
                JAC.AngularCoefficientsRatip2013.load_csl(csfa)
                t_coeffs = JAC.AngularCoefficientsRatip2013.angular_coefficients_pair_1p(rank,1,1)
            else
                JAC.AngularCoefficientsRatip2013.load_csl(csfa, csfb)
                t_coeffs = JAC.AngularCoefficientsRatip2013.angular_coefficients_pair_1p(rank,1,2)
            end
        else    error("Unsupported keystring = ", sa)
        end
        return t_coeffs
    end

    """
    `Basics.compute("angular coefficients: 1-p, Grasp92", parity, rank::Integer, csfa::CsfR, csfb::CsfR)`
        ... to compute the the angular coefficients in the decomposition of a (reduced) many-electron matrix element with a general
            single-particle operator of the given parity and rank by a call to the Fortran procedure `mct_generate_coefficients`
            of the RATIP program; an `Array{AngularTcoeff,1}` is returned.
    """
    function Basics.compute(sa::String, parity, rank::Integer, csfa::CsfR, csfb::CsfR)
        if      sa == "angular coefficients: 1-p, Grasp92"
            if csfa == csfb
                subshells  = JAC.AngularCoefficientsRatip2013.load_csl(csfa)
                mct_coeffs = JAC.AngularCoefficientsRatip2013.mct_generate_coefficients(1, 1, 1, 1, Int32(parity), rank)
            else
                # Add 'zeros' to the fields if the length of occupation does not agree
                if       ( nz = length(csfa.occupation) - length(csfb.occupation) ) <  0   csfa = Basics.addZerosToCsfR( -nz, csfa)
                elseif   ( nz = length(csfa.occupation) - length(csfb.occupation) ) >  0   csfb = Basics.addZerosToCsfR(  nz, csfb)    end
                subshells  = JAC.AngularCoefficientsRatip2013.load_csl(csfa, csfb)
                mct_coeffs = JAC.AngularCoefficientsRatip2013.mct_generate_coefficients(1, 1, 2, 2, Int32(parity), rank)
            end
        else    error("Unsupported keystring = ", sa)
        end
        return map(t -> JAC.AngularCoefficientsRatip2013.AngularTcoeff(t, subshells), mct_coeffs) # Convert Fmctcoefficient to AngularTcoeff
    end

    """
    `Basics.compute("matrix: CI, J^P symmetry", JP::LevelSymmetry, basis::Basis, nuclearModel::Nuclear.Model, grid::Radial.Grid,
                                               settings::AsfSettings; printout::Bool=true)`
        ... to compute the CI matrix for a given J^P symmetry block of basis and by making use of the nuclear model and the grid;
            a matrix::Array{Float64,2} is returned.
    """
    function Basics.compute(sa::String, JP::LevelSymmetry, basis::Basis, nuclearModel::Nuclear.Model, grid::Radial.Grid,
                            settings::AsfSettings; printout::Bool=true)
        !(sa == "matrix: CI, J^P symmetry")   &&   error("Not supported keystring")

        # Determine the dimension of the CI matrix and the indices of the CSF with J^P symmetry in the basis
        idx_csf = Int64[]
        for  idx = 1:length(basis.csfs)
            if  basis.csfs[idx].J == JP.J   &&   basis.csfs[idx].parity == JP.parity    push!(idx_csf, idx)    end
        end
        n = length(idx_csf)
        if printout    print("Compute CI matrix of dimension $n x $n for the symmetry $(string(JP.J))^$(string(JP.parity)) ...")    end

        # Generate an effective nuclear charge Z(r) on the given grid
        potential = Nuclear.nuclearPotential(nuclearModel, grid)
        if  settings.qedModel in [QedPetersburg(), QedSydney()]
            meanPot = potential
            ## meanPot = compute("radial potential: Dirac-Fock-Slater", grid, basis)
            ## meanPot = Basics.add(potential, meanPot)
        end

        matrix = zeros(Float64, n, n)
        keep = true
        InteractionStrength.XL_Coulomb_reset_storage(keep, printout=false)
        InteractionStrength.XL_Breit_reset_storage(keep)
        print(">> performCI() ... ")
        @time   for  r = 1:n
            for  s = 1:n
                wa = compute("angular coefficients: e-e, Ratip2013", basis.csfs[idx_csf[r]], basis.csfs[idx_csf[s]])
                me = 0.
                for  coeff in wa[1]
                    jj = Basics.subshell_2j(basis.orbitals[coeff.a].subshell)
                    me = me + coeff.T * sqrt( jj + 1) * RadialIntegrals.GrantIab(basis.orbitals[coeff.a], basis.orbitals[coeff.b], grid, potential)
                    if  settings.qedModel != NoneQed()
                        me = me + InteractionStrengthQED.qedLocal(basis.orbitals[coeff.a], basis.orbitals[coeff.b], nuclearModel,
                                                                  settings.qedModel, meanPot, grid)
                    end
                end

                for  coeff in wa[2]
                    if  settings.eeInteractionCI in [CoulombInteraction(), CoulombBreit()]
                        me = me + coeff.V * InteractionStrength.XL_Coulomb(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid, keep=keep)
                    elseif  false
                        xl1 = InteractionStrength.XL_Coulomb(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid, keep=false)
                        xl2 = InteractionStrength.XL_Coulomb_WO(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid)
                        if abs(xl1 - xl2) > 1.0e-12  println("XL_Coulomb differ: $xl1   $xl2")      end
                    elseif  false
                        me = me + coeff.V * InteractionStrength.XL_Coulomb_WO(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid)
                    end

                    if  settings.eeInteractionCI in [BreitInteraction(), CoulombBreit()]
                        me = me + coeff.V * InteractionStrength.XL_Breit(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid, keep=keep)
                    elseif  false
                        xl1 = InteractionStrength.XL_Breit(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid, keep=false)
                        xl2 = InteractionStrength.XL_Breit_WO(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid)
                        if abs(xl1 - xl2) > 1.0e-12  println("XL_Breit differ: $xl1   $xl2")      end
                    elseif  false
                        me = me + coeff.V * InteractionStrength.XL_Breit_WO(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid)
                    end
                end
                matrix[r,s] = me
            end
        end
        if printout    println("   ... done.")    end

        return( matrix )
    end

    """
    `Basics.compute("matrix: CI for plasma, J^P symmetry", JP::LevelSymmetry, basis::Basis, nuclearModel::Nuclear.Model, grid::Radial.Grid,
                                                          settings::AsfSettings, plasmaSettings::PlasmaShift.Settings; printout::Bool=true)`
        ... to compute the CI matrix for a given J^P symmetry block of basis and by making use of the nuclear model and the grid;
            a matrix::Array{Float64,2} is returned.
    """
    function Basics.compute(sa::String, JP::LevelSymmetry, basis::Basis, nuclearModel::Nuclear.Model, grid::Radial.Grid,
                                settings::AsfSettings, plasmaSettings::PlasmaShift.Settings; printout::Bool=true)
        !(sa == "matrix: CI for plasma, J^P symmetry")   &&   error("Not supported keystring")

        # Determine the dimension of the CI matrix and the indices of the CSF with J^P symmetry in the basis
        idx_csf = Int64[]
        for  idx = 1:length(basis.csfs)
            if  basis.csfs[idx].J == JP.J   &&   basis.csfs[idx].parity == JP.parity    push!(idx_csf, idx)    end
        end
        n = length(idx_csf)

        if  settings.eeInteractionCI in [BreitInteraction(), CoulombBreit()]
                                             error("No Breit interaction supported for plasma computations; use breitCI=false  in the asfSettings.")    end
        if  settings.qedModel != NoneQed()   error("No QED estimates supported for plasma computations; use qedModel=NoneQed()  in the asfSettings.")   end

        # Now distinguis the CI matrix for different plasma models
        if  plasmaSettings.plasmaModel == PlasmaShift.DebyeHueckel()
            if printout    print("Compute DebyeHueckel-CI matrix of dimension $n x $n for the symmetry $(string(JP.J))^$(string(JP.parity)) ...")    end

            # Generate an effective nuclear charge Z(r) for a screened Debye-Hueckel potential on the given grid
            potential = Nuclear.nuclearPotentialDH(nuclearModel, grid, plasmaSettings.lambdaDebye)

            matrix = zeros(Float64, n, n)
            for  r = 1:n
                for  s = 1:n
                    wa = compute("angular coefficients: e-e, Ratip2013", basis.csfs[idx_csf[r]], basis.csfs[idx_csf[s]])
                    me = 0.
                    for  coeff in wa[1]
                        jj = Basics.subshell_2j(basis.orbitals[coeff.a].subshell)
                        me = me + coeff.T * sqrt( jj + 1) * RadialIntegrals.GrantIab(basis.orbitals[coeff.a], basis.orbitals[coeff.b], grid, potential)
                    end

                    for  coeff in wa[2]
                        if  settings.eeInteractionCI in [CoulombInteraction(), CoulombBreit()]
                            me = me + coeff.V * InteractionStrength.XL_Coulomb_DH(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid,
                                                                                plasmaSettings.lambdaDebye)   end
                    end
                    matrix[r,s] = me
                end
            end
        else
            error("Unsupported plasma model = $(plasmaSettings.plasmaModel)")
        end

        if printout    println("   ... done.")    end

        return( matrix )
    end

    """
    `Basics.compute("radial orbital: NR, Bunge (1993)", subshell::Subshell, Z::Int64)`
        ... to compute a radial orbital::Orbital for the given subshell and nuclear charge by using the Roothan-Hartree-Fock data by
            Bunge et al., Atomic Data and Nuclear Data Tables 53 (1993) 113, as obtained for a non-relativistic RHF computation of the
            neutral atom. These functions are used a large component, and the small component is obtained from the kinetic balance
            condition. Radial orbitals can be obtained for the ground-state configuration for all elements with 2 <= Z <= 54.
    `Basics.compute("radial orbital: NR, McLean (1981)", subshell::Subshell, Z::Int64)`
        ... to compute a radial orbital::Orbital for the given subshell and nuclear charge by using the Roothan-Hartree-Fock data by
            McLean and McLean, Atomic Data and Nuclear Data Tables 26 (1981) 197., as obtained for a non-relativistic RHF computation of
            the neutral atom. These functions are used a large component, and the small component is obtained from the kinetic balance
            condition. Radial orbitals can be obtained for the ground-state configuration for all elements with 55 <= Z <= 92.
    """
    function Basics.compute(sa::String, subshell::Subshell, Z::Int64)
        if      sa == "radial orbital: NR, Bunge (1993)"          wa = Radial.OrbitalBunge1993(subshell,Z)
        elseif  sa == "radial orbital: NR, McLean (1981)"         wa = Radial.OrbitalMcLean1981(subshell,Z)
        elseif  sa == "radial orbital: hydrogenic"                wa = Radial.OrbitalHydrogenic(subshell,Z)
        elseif  sa == "radial orbital: Thomas-Fermi"              wa = Radial.OrbitalThomasFermi(subshell,Z)
        else    error("Unsupported keystring = $sa ")
        end

        return( wa )
    end

    """
    `Basics.compute("radial potential: core-Hartree", grid::Radial.Grid, level::Level)`
        ... to compute a (radial) core-Hartree potential from the given list of orbitals; cf. Basics.computePotentialCoreHartree.
            A potential::RadialPotential is returned.

    `Basics.compute("radial potential: Kohn-Sham", grid::Radial.Grid, level::Level)`
        ... to compute a (radial) Kohn-Sham potential from the given list of orbitals; cf. Basics.computePotentialKohnSham.
            A potential::RadialPotential is returned.

    `Basics.compute("radial potential: Dirac-Fock-Slater", grid::Radial.Grid, level::Level)`
        ... to compute a (radial) Dirac-Fock-Slater potential from the given list of orbitals; this potential is rather simple but
            includes some undesired self-interaction and exhibits an asymptotic behaviour. Cf. Basics.computePotentialDFS.
            A potential::RadialPotential is returned.

    `Basics.compute("radial potential: Dirac-Fock-Slater", grid::Radial.Grid, basis::Basis)`
        ... to compute a (radial) Dirac-Fock-Slater potential but for mean occupation of the given basis. Cf. Basics.computePotentialDFS.
            A potential::RadialPotential is returned.

    `Basics.compute("radial potential: extended Hartree", grid::Radial.Grid, level::Level)`
        ... to compute a (radial) extended Hartree potential from the given list of orbitals; it is a local potential that is based on
            direct and exchange coefficients as obtained from the configuration-averaged energy. Cf. Basics.computePotentialExtendedHartree.
            A potential::RadialPotential is returned.
    """
    function Basics.compute(sa::String, grid::Radial.Grid, level::Level)
        if      sa == "radial potential: core-Hartree"               wa = Basics.computePotentialCoreHartree(grid, level)
        elseif  sa == "radial potential: Hartree"                    wa = Basics.computePotentialHartree(grid, level)
        elseif  sa == "radial potential: Hartree-Slater"             wa = Basics.computePotentialHartreeSlater(grid, level)
        elseif  sa == "radial potential: Kohn-Sham"                  wa = Basics.computePotentialKohnSham(grid, level)
        elseif  sa == "radial potential: Dirac-Fock-Slater"          wa = Basics.computePotentialDFS(grid, level)
        elseif  sa == "radial potential: extended-Hartree"           wa = Basics.computePotentialExtendedHartree(grid, level)
        else    error("Unsupported keystring = $sa ")
        end

        return( wa )
    end

    function Basics.compute(sa::String, grid::Radial.Grid, basis::Basis)
        if      sa == "radial potential: Dirac-Fock-Slater"          wa = Basics.computePotentialDFS(grid, basis)
        else    error("Unsupported keystring = $sa ")
        end

        return( wa )
    end

    """
    `Basics.computeDiracEnergy(sh::Subshell, Z::Float64)`
        ... computes the Dirac energy for the hydrogenic subshell sh and for a point-like nucleus with nuclear charge Z;
            a energy::Float64 in atomic units and without the rest energy of the electron is returned. That is the binding
            energy of a 1s_1/2 electron for Z=1 is -0.50000665.
    """
    function Basics.computeDiracEnergy(sh::Subshell, Z::Float64)
        if  Z <= 0.1    error("Requires nuclear charge Z >= 0.1")    end
        # Compute the energy from the Dirac formula
        jPlusHalf = (Basics.subshell_2j(sh) + 1) / 2;   nr = sh.n - jPlusHalf;    alpha = Defaults.getDefaults("alpha")
        wa = sqrt(jPlusHalf^2 - Z^2 * alpha^2 )
        wa = sqrt(1.0 + Z^2 * alpha^2 / (nr + wa)^2)
        wa = Defaults.getDefaults("speed of light: c")^2 * (1/wa - 1.0)
        return( wa )
    end

    """
    `Basics.computeMeanSubshellOccupation(sh::Subshell, levels::Array{Level,1})`
        ... computes the mean subshell occupation for the subshell sh and for the given levels; a q::Float64 is returned.
    """
    function Basics.computeMeanSubshellOccupation(sh::Subshell, levels::Array{Level,1})
        q = 0.
        if  length(levels) < 1   error("stop a")    end
        for  level in levels
            subshells = level.basis.subshells;    nsh = 0
            for  i = 1:length(level.basis.subshells)    if  sh == subshells[i]   nsh = i;   break   end    end
            if   nsh == 0   error("Subshell not found in CSF basis.")   end
            #
            for  i = 1:length(level.basis.csfs)   q = q + abs(level.mc[i])^2 * level.basis.csfs[i].occupation[nsh]    end
        end
        return( q/length(levels) )
    end

    """
    `Basics.computeMeanSubshellOccupation(sh::Subshell, basis::Basis)`
        ... computes the mean subshell occupation for the subshell sh and for the given CSF in the basis; a q::Float64 is returned.
    """
    function Basics.computeMeanSubshellOccupation(sh::Subshell, basis::Basis)
        q = 0.
        for  csf in basis.csfs
            subshells = basis.subshells;    nsh = 0
            for  i = 1:length(basis.subshells)    if  sh == subshells[i]   nsh = i;   break   end    end
            if   nsh == 0   error("Subshell not found in CSF basis.")   end
            #
            for  i = 1:length(basis.csfs)   q = q + basis.csfs[i].occupation[nsh]    end
        end
        return( q/length(basis.csfs) )
    end

    """
    `Basics.computeMultipletForGreenApproach(approach::AtomicState.SingleCSFwithoutCI, basis::Basis, nModel::Nuclear.Model, grid::Radial.Grid,
                                             asfSettings::AsfSettings, greenSettings::GreenSettings; printout::Bool=false)`
        ... computes the (Green channel) multiplet from the given basis with the SingleCSFwithoutCI approach.
    """
    function Basics.computeMultipletForGreenApproach(approach::AtomicState.SingleCSFwithoutCI, basis::Basis, nModel::Nuclear.Model,
                                                     grid::Radial.Grid, asfSettings::AsfSettings, greenSettings::GreenSettings; printout::Bool=false)
        print("Compute a Green function multiplet in $approach approach ... ")
        # In the SingleCSFwithoutCI, only the diagonal ME are included into the Hamiltonian matrix
        # (1) Check that all CSF have same (level) symmetry; issue an error if not
        sym = LevelSymmetry(0, Basics.plus)
        for  (r, csf)  in enumerate(basis.csfs)
            if     r == 1   sym  = LevelSymmetry( csf.J, csf.parity )
            elseif          sym != LevelSymmetry( csf.J, csf.parity )  error("stop a")
            end
        end

        # (2) Compute Hamiltonian matrix with only diagonal ME with given settings; issue a message which flags are considered.
        if  asfSettings.qedModel != NoneQed()   println("   ++ No QED terms included for this Green function approach.")               end
        if  asfSettings.jjLS.makeIt             println("   ++ No jj-LS transformation included for this Green function approach.")    end
        potential = Nuclear.nuclearPotential(nModel, grid)
        ncsf      = length(basis.csfs);    matrix = zeros(Float64, ncsf, ncsf)
        for  (r, csf)  in enumerate(basis.csfs)
            wa = compute("angular coefficients: e-e, Ratip2013", basis.csfs[r], basis.csfs[r])
            me = 0.
            for  coeff in wa[1]
                jj = Basics.subshell_2j(basis.orbitals[coeff.a].subshell)
                me = me + coeff.T * sqrt( jj + 1) * RadialIntegrals.GrantIab(basis.orbitals[coeff.a], basis.orbitals[coeff.b], grid, potential)
            end

            for  coeff in wa[2]
                if  asfSettings.coulombCI
                    me = me + coeff.V * InteractionStrength.XL_Coulomb(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid)   end
                if  asfSettings.breitCI
                    me = me + coeff.V * InteractionStrength.XL_Breit(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                               basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid)     end
            end
            matrix[r,r] = me
        end

        # (3) Diagonalize matrix with Julia;   assign a multiplet
        eigen  = Basics.diagonalize("matrix: LinearAlgebra", matrix)
        levels = Level[]
        for  ev = 1:length(eigen.values)
            level = Level( sym.J, AngularM64(sym.J.num//sym.J.den), sym.parity, 0, eigen.values[ev], 0., true, basis, eigen.vectors[ev] )
            push!( levels, level)
        end
        println("done with $(length(levels)) levels")

        multiplet = Multiplet("SingleCSFwithoutCI multiplet for $sym", levels)
        return( multiplet )
    end

    """
    `Basics.computeMultipletForGreenApproach(approach::AtomicState.CoreSpaceCI, basis::Basis, nModel::Nuclear.Model, grid::Radial.Grid,
                                             asfSettings::AsfSettings, greenSettings::GreenSettings; printout::Bool=false)`
        ... computes the (Green channel) multiplet from the given basis with the CoreSpaceCI approach in which the electron-electron
            interaction is taken into account only between the bound-state orbitals.
    """
    function Basics.computeMultipletForGreenApproach(approach::AtomicState.CoreSpaceCI, basis::Basis, nModel::Nuclear.Model,
                                                     grid::Radial.Grid, asfSettings::AsfSettings, greenSettings::GreenSettings; printout::Bool=false)
        print("Compute a Green function multiplet in $approach approach ... ")
        # In the SingleCSFwithoutCI, only the diagonal ME are included into the Hamiltonian matrix
        # (1) Check that all CSF have same (level) symmetry; issue an error if not
        sym = LevelSymmetry(0, Basics.plus)
        for  (r, csf)  in enumerate(basis.csfs)
            if     r == 1   sym  = LevelSymmetry( csf.J, csf.parity )
            elseif          sym != LevelSymmetry( csf.J, csf.parity )  error("stop a")
            end
        end

        # (2) Compute Hamiltonian matrix with only diagonal ME with given settings; issue a message which flags are considered.
        if  asfSettings.qedModel != NoneQed()   println("   ++ No QED terms included for this Green function approach.")               end
        if  asfSettings.jjLS.makeIt             println("   ++ No jj-LS transformation included for this Green function approach.")    end
        potential = Nuclear.nuclearPotential(nModel, grid)
        ncsf      = length(basis.csfs);    matrix = zeros(Float64, ncsf, ncsf)
        for  (r, rcsf)  in enumerate(basis.csfs)
            for  (s, scsf)  in enumerate(basis.csfs)
                wa = compute("angular coefficients: e-e, Ratip2013", basis.csfs[r], basis.csfs[s])
                me = 0.
                for  coeff in wa[1]
                    jj = Basics.subshell_2j(basis.orbitals[coeff.a].subshell)
                    if  basis.orbitals[coeff.a].isBound   &&   basis.orbitals[coeff.b].isBound
                        me = me + coeff.T * sqrt( jj + 1) * RadialIntegrals.GrantIab(basis.orbitals[coeff.a], basis.orbitals[coeff.b], grid, potential)
                    end
                end

                for  coeff in wa[2]
                    if  basis.orbitals[coeff.a].isBound   &&   basis.orbitals[coeff.b].isBound   &&
                        basis.orbitals[coeff.c].isBound   &&   basis.orbitals[coeff.d].isBound
                        if  asfSettings.coulombCI
                            me = me + coeff.V * InteractionStrength.XL_Coulomb(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid)   end
                        if  asfSettings.breitCI
                            me = me + coeff.V * InteractionStrength.XL_Breit(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid)     end
                    end
                end
                matrix[r,s] = me
            end
        end

        # (3) Diagonalize matrix with Julia;   assign a multiplet
        eigen  = Basics.diagonalize("matrix: LinearAlgebra", matrix)
        levels = Level[]
        for  ev = 1:length(eigen.values)
            level = Level( sym.J, AngularM64(sym.J.num//sym.J.den), sym.parity, 0, eigen.values[ev], 0., true, basis, eigen.vectors[ev] )
            push!( levels, level)
        end
        println("done with $(length(levels)) levels")

        multiplet = Multiplet("CoreSpaceCI multiplet for $sym", levels)
        return( multiplet )
    end

    """
    `Basics.computeMultipletForGreenApproach(approach::AtomicState.DampedSpaceCI, basis::Basis, nModel::Nuclear.Model, grid::Radial.Grid,
                                             asfSettings::AsfSettings, greenSettings::GreenSettings; printout::Bool=false)`
        ... computes the (Green channel) multiplet from the given basis with the DampedSpaceCI approach in which the electron-electron
            interaction strength are all damped by the factor exp(- dampingTau * r).
    """
    function Basics.computeMultipletForGreenApproach(approach::AtomicState.DampedSpaceCI, basis::Basis, nModel::Nuclear.Model,
                                                     grid::Radial.Grid, asfSettings::AsfSettings, greenSettings::GreenSettings;
                                                     printout::Bool=false)
        print("Compute a Green function multiplet in $approach approach ... ")
        # In the SingleCSFwithoutCI, only the diagonal ME are included into the Hamiltonian matrix
        # (1) Check that all CSF have same (level) symmetry; issue an error if not
        sym = LevelSymmetry(0, Basics.plus)
        for  (r, csf)  in enumerate(basis.csfs)
            if     r == 1   sym  = LevelSymmetry( csf.J, csf.parity )
            elseif          sym != LevelSymmetry( csf.J, csf.parity )  error("stop a")
            end
        end

        # (2) Compute Hamiltonian matrix with only diagonal ME with given settings; issue a message which flags are considered.
        if  asfSettings.qedModel != NoneQed()   println("   ++ No QED terms included for this Green function approach.")               end
        if  asfSettings.jjLS.makeIt             println("   ++ No jj-LS transformation included for this Green function approach.")    end
        potential = Nuclear.nuclearPotential(nModel, grid)
        ncsf      = length(basis.csfs);    matrix = zeros(Float64, ncsf, ncsf);    tau = greenSettings.dampingTau
        for  (r, rcsf)  in enumerate(basis.csfs)
            for  (s, scsf)  in enumerate(basis.csfs)
                wa = compute("angular coefficients: e-e, Ratip2013", basis.csfs[r], basis.csfs[s])
                me = 0.
                for  coeff in wa[1]
                    jj = Basics.subshell_2j(basis.orbitals[coeff.a].subshell)
                    me = me + coeff.T * sqrt( jj + 1) * RadialIntegrals.GrantIabDamped(tau, basis.orbitals[coeff.a], basis.orbitals[coeff.b], grid, potential)
                end

                for  coeff in wa[2]
                    if  asfSettings.eeInteraction in [CoulombInteraction(), CoulombBreit()]
                        me = me + coeff.V * InteractionStrength.XL_CoulombDamped(tau, coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid)   end
                    if  asfSettings.eeInteraction in [BreitInteraction(), CoulombBreit()]
                        me = me + coeff.V * InteractionStrength.XL_BreitDamped(tau, coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid)     end
                end
                matrix[r,s] = me
            end
        end

        # (3) Diagonalize matrix with Julia;   assign a multiplet
        eigen  = Basics.diagonalize("matrix: LinearAlgebra", matrix)
        levels = Level[]
        for  ev = 1:length(eigen.values)
            level = Level( sym.J, AngularM64(sym.J.num//sym.J.den), sym.parity, 0, eigen.values[ev], 0., true, basis, eigen.vectors[ev] )
            push!( levels, level)
        end
        println("done with $(length(levels)) levels")

        multiplet = Multiplet("DampedSpaceCI multiplet for $sym", levels)
        return( multiplet )
    end

    """
    `Basics.computePotentialCoreHartree(grid::Radial.Grid, level::Level)`
        ... to compute a (radial) core-Hartree potential for the given level; a potential::RadialPotential is returned.
            The core-Hartree potential is defined by

                V_CH(r) = int_0^infty dr'  rho_c(r') / r_>

            with r_> = max(r,r')  and  rho_c(r) =  sum_a (Pa^2(r) + Qa^2(r)) the charge density of the core-electrons.
            An Radial.Potential with -r * V_CH(r)  is returned to be consistent with an effective charge Z(r).
    """
    function Basics.computePotentialCoreHartree(grid::Radial.Grid, level::Level)
        basis = level.basis;    npoints = grid.NoPoints
        println("coreSubshells = $(basis.coreSubshells);    subshells = $(basis.subshells)")
        rhoc = zeros( npoints );    wx = zeros( npoints );    wb = zeros( npoints )
        # Compute the charge density of the core orbitals for the given level
        for  sh in basis.coreSubshells
        ## for  sh in basis.subshells
            orb  = basis.orbitals[sh]
            occ  = Basics.computeMeanSubshellOccupation(sh, [level])    ## (Basics.subshell_2j(sh) + 1.0)  for closed-shells
            nrho = length(orb.P)
            for    i = 1:nrho   rhoc[i] = rhoc[i] + occ * (orb.P[i]^2 + orb.Q[i]^2)    end
        end
        # Define the integrant and take the integrals
        for i = 1:npoints
            for j = 1:npoints    rg = max( grid.r[i],  grid.r[j]);    wx[j] = rhoc[j]/rg   end
            wb[i] = RadialIntegrals.V0(wx, npoints, grid::Radial.Grid)
        end
        # Define the potential with regard to Z(r)
        for i = 1:npoints    wb[i] = - wb[i] * grid.r[i]   end
        #
        wc = Radial.Potential("CoreHartree", wb, grid)
        return( wc )
    end

    """
    `Basics.computePotentialHartree(grid::Radial.Grid, level::Level)`
        ... to compute a (radial) Hartree potential for the given level; a potential::RadialPotential is returned.
            Here, the Hartree potential is simply defined by

                V_H(r) = - Sum_a  q_a^bar  Y0_aa(r)

            An Radial.Potential with V_H(r)  is returned that is consistent with an effective charge Z(r).
    """
    function Basics.computePotentialHartree(grid::Radial.Grid, level::Level)
        basis = level.basis;    npoints = grid.NoPoints;    wx = zeros( npoints )
        # Sum_a ...
        for  a in basis.subshells
            occa = Basics.computeMeanSubshellOccupation(a, [level])
            orba = basis.orbitals[a];   nrho = length(orba.P);      rhoaa  = zeros(nrho)
            for  i = 1:nrho    rhoaa[i] = orba.P[i]^2 + orba.Q[i]^2    end
            for  i = 1:nrho    wx[i]    = wx[i] - occa * RadialIntegrals.Yk_ab(0, grid.r[i], rhoaa, nrho, grid)   end
            Yk = RadialIntegrals.Yk_ab(0, grid.r[nrho], rhoaa, nrho, grid)
            for  i = nrho+1:npoints       wx[i] = wx[i] - occa * Yk    end
        end

        wc = Radial.Potential("Hartree", wx, grid)
        return( wc )
    end

    """
    `Basics.computePotentialHartreeSlater(grid::Radial.Grid, level::Level)`
        ... to compute a (radial) Hartree-Slater potential for the given level; a potential::RadialPotential is returned.
            Here, the Hartree-Slater potential is defined by

                                                        3  [     3          ]^1/3      r
                V_HS(r) = = - Sum_a  q_a^bar  Y0_aa(r)  +  -  [ ---------  rho ]       *  -
                                                        2  [ 4pi^2 r^2      ]          2

            with rho(r) = Sum_a  q_a^bar * (Pa^2(r) + Qa^2(r))   ... charge density of all electrons.
            An Radial.Potential with V_HS(r)  is returned that is consistent with an effective charge Z(r).
    """
    function Basics.computePotentialHartreeSlater(grid::Radial.Grid, level::Level)
        basis = level.basis;    npoints = grid.NoPoints;    wx = zeros( npoints );    rho = zeros( npoints )
        # Sum_a ...
        for  a in basis.subshells
            occa = Basics.computeMeanSubshellOccupation(a, [level])
            orba = basis.orbitals[a];   nrho = length(orba.P);      rhoaa  = zeros(nrho)
            for  i = 1:nrho    rhoaa[i] = orba.P[i]^2 + orba.Q[i]^2    end
            for  i = 1:nrho    wx[i]    = wx[i]  - occa * RadialIntegrals.Yk_ab(0, grid.r[i], rhoaa, nrho, grid)   end
            for  i = 1:nrho    rho[i]   = rho[i] + occa * rhoaa[i]     end
            Yk = RadialIntegrals.Yk_ab(0, grid.r[nrho], rhoaa, nrho, grid)
            for  i = nrho+1:npoints       wx[i] = wx[i] - occa * Yk    end
        end
        ## for  i = 2:npoints     rho[i]   = rho[i] / (4pi * grid.r[i])   end;       rho[1] = 0.
        for  i = 2:npoints     wx[i]    = wx[i] + (3/2) * (3/(4pi^2 * grid.r[i]^2) * rho[i])^(1/3) * grid.r[i] / 2   end
        #
        wc = Radial.Potential("Hartree-Slater", wx, grid)
        return( wc )
    end

    """
    `Basics.computePotentialKohnSham(grid::Radial.Grid, level::Level)`
        ... to compute a (radial) Kohn-Sham potential for the given level; a potential::RadialPotential is returned.
            The Kohn-Sham potential is defined by

                                                                2     [   81                   ]^(1/3)
                V_KS(r) =  int_0^infty dr'  rho_t(r') / r_>  --  ---  [--------   r * rho_t(r) ]
                                                                3 r   [ 32 pi^2                ]

            with r_> = max(r,r')  and  rho_t(r) = sum_a (Pa^2(r) + Qa^2(r))   ... charge density of all electrons.
            An Radial.Potential with -r * V_KS(r)  is returned to be consistent with an effective charge Z(r).
    """
    function Basics.computePotentialKohnSham(grid::Radial.Grid, level::Level)
        basis = level.basis;       npoints = grid.NoPoints
        rhot = zeros( npoints );   wb = zeros( npoints );   wx = zeros( npoints )
        # Compute the charge density of the core orbitals for the given level
        for  sh in basis.subshells
            orb  = basis.orbitals[sh]
            occ  = Basics.computeMeanSubshellOccupation(sh, [level])
            nrho = length(orb.P)
            for    i = 1:nrho   rhot[i] = rhot[i] + occ * (orb.P[i]^2 + orb.Q[i]^2)    end
        end
        # Define the integrant and take the integrals (without alpha)
        for i = 1:npoints
            for j = 1:npoints    rg = max( grid.r[i],  grid.r[j]);    wx[j] = rhot[j]/rg   end
            wb[i] = JAC.RadialIntegrals.V0(wx, npoints, grid::Radial.Grid)
            wb[i] = wb[i] - 2 / (3grid.r[i]) * (81 /(32 * pi^2) * grid.r[i] * rhot[i])^(1/3)
        end
        # Define the potential with regard to Z(r)
        for i = 2:npoints    wb[i] = - wb[i] * grid.r[i]   end;    wb[1] = 0.
        #
        wc = JAC.Radial.Potential("Kohn-Sham", wb, grid)
        return( wc )
    end

    """
    `Basics.computePotentialDFS(grid::Radial.Grid, level::Level)`
        ... to compute a (radial) Dirac-Fock-Slater potential for the given level; a potential::RadialPotential is returned.
            The Dirac-Fock-Slater potential is defined by

                                                                [   3                   ]^(1/3)
                V_DFS(r) = int_0^infty dr'  rho_t(r') / r_>  -- [-----------   rho_t(r) ]
                                                                [ 4 pi^2 r^2            ]

            with r_> = max(r,r')  and  rho_t(r) = sum_a (Pa^2(r) + Qa^2(r))   ... charge density of all electrons.
            An Radial.Potential with -r * V_DFS(r)  is returned to be consistent with an effective charge Z(r).
    """
    function Basics.computePotentialDFS(grid::Radial.Grid, level::Level)
        basis = level.basis;    npoints = grid.NoPoints
        rhot = zeros( npoints );    wb = zeros( npoints );    wx = zeros( npoints )
        # Compute the charge density of the core orbitals for the given level
        for  sh in basis.subshells
            orb  = basis.orbitals[sh]
            occ  = Basics.computeMeanSubshellOccupation(sh, [level])
            nrho = length(orb.P)
            for    i = 1:nrho   rhot[i] = rhot[i] + occ * (orb.P[i]^2 + orb.Q[i]^2)    end
        end
        # Define the integrant and take the integrals
        for i = 1:npoints
            for j = 1:npoints    rg = max( grid.r[i],  grid.r[j]);    wx[j] = rhot[j]/rg   end
            wb[i] = JAC.RadialIntegrals.V0(wx, npoints, grid::Radial.Grid)
            wb[i] = wb[i] - (3 /(4*pi^2 * grid.r[i]^2) * rhot[i])^(1/3)
        end
        # Define the potential with regard to Z(r)
        for i = 2:npoints    wb[i] = - wb[i] * grid.r[i]   end;    wb[1] = 0.
        #
        wc = JAC.Radial.Potential("DFS", wb, grid)
        return( wc )
    end

    """
    `Basics.computePotentialDFS(grid::Radial.Grid, basis::Basis)`
        ... to compute the same but for the mean occupation of the orbitals in the given basis.
    """
    function Basics.computePotentialDFS(grid::Radial.Grid, basis::Basis)
        npoints = grid.NoPoints
        rhot = zeros( npoints );    wb = zeros( npoints );    wx = zeros( npoints )
        # Compute the charge density of the core orbitals for the given level
        for  sh in basis.subshells
            orb  = basis.orbitals[sh]
            occ  = Basics.computeMeanSubshellOccupation(sh, basis)
            nrho = length(orb.P)
            for    i = 1:nrho   rhot[i] = rhot[i] + occ * (orb.P[i]^2 + orb.Q[i]^2)    end
        end
        # Define the integrant and take the integrals
        for i = 1:npoints
            for j = 1:npoints    rg = max( grid.r[i],  grid.r[j]);    wx[j] = rhot[j]/rg   end
            wb[i] = JAC.RadialIntegrals.V0(wx, npoints, grid::Radial.Grid)
            wb[i] = wb[i] - (3 /(4*pi^2 * grid.r[i]^2) * rhot[i])^(1/3)
        end
        # Define the potential with regard to Z(r)
        for i = 2:npoints    wb[i] = - wb[i] * grid.r[i]   end;    wb[1] = 0.
        #
        wc = JAC.Radial.Potential("DFS for CSF basis", wb, grid)
        return( wc )
    end

    """
    `Basics.computePotentialExtendedHartree(grid::Radial.Grid, level::Level)`
        ... to compute a (radial) extended-Hartree potential for the given level; a potential::RadialPotential is returned.
            The extended-Hartreepotential is defined by Gu (2008, Eq.(8)) and is based on the Y^k_ab (r) function and various
            direct and exchange coefficients as the arise from average energy of the configuration.
    """
    function Basics.computePotentialExtendedHartree(grid::Radial.Grid, level::Level)
        basis = level.basis;    npoints = grid.NoPoints
        rhoab = zeros( npoints );    wx = zeros( npoints );    wg = zeros( npoints )
        # First term Sum_ab ...
        for  a in basis.subshells
            occa = Basics.computeMeanSubshellOccupation(a, [level])
            orba = basis.orbitals[a]
            nrho = length(orba.P);      rhoa  = zeros(nrho)
            for  i = 1:nrho    rhoa[i] = orba.P[i]^2 + orba.Q[i]^2    end
            for  b in basis.subshells
                occb = Basics.computeMeanSubshellOccupation(b, [level])
                orbb = basis.orbitals[b]
                nrho = length(orbb.P);   rhobb = zeros(nrho)
                if  a != b  wc = occa*occb   else   wc = occa * (occb -1.0)   end
                if  wc < 0.  error("Improper potential for subshells with mean occupation < 1; stop a")    end
                for  i = 1:nrho    rhobb[i] = orbb.P[i]^2 + orbb.Q[i]^2       end
                for  i = 1:length(orba.P)    wx[i] = wx[i] + wc * JAC.RadialIntegrals.Yk_ab(0, grid.r[i], rhobb, nrho, grid) * rhoa[i]   end
            end
        end
        # Second term Sum_a ...
        for  a in basis.subshells
            occa = Basics.computeMeanSubshellOccupation(a, [level])
            orba = basis.orbitals[a]
            nrho = length(orba.P);      rhoa  = zeros(nrho)
            for  i = 1:nrho    rhoa[i] = orba.P[i]^2 + orba.Q[i]^2   end
            wc = occa * (occa - 1.0)
            if  wc < 0.  error("Improper potential for subshells with mean occupation < 1; stop b")    end
            for  k = 1:10
                fk = JAC.AngularMomentum.fk_ab(a,a);    if  fk == 0   cycle    end
                for  i = 1:length(orba.P)    wx[i] = wx[i] + wc * fk * JAC.RadialIntegrals.Yk_ab(k, grid.r[i], rhoa, nrho, grid) * rhoa[i]   end
            end
        end
        # Third term Sum_a != b ...
        for  a in basis.subshells,   b in basis.subshells
            if  a == b   continue   end
            occa = Basics.computeMeanSubshellOccupation(a, [level])
            occb = Basics.computeMeanSubshellOccupation(b, [level])
            orba = basis.orbitals[a]
            orbb = basis.orbitals[b]
            nrho = min(length(orba.P), length(orbb.P));      rhoab  = zeros(nrho)
            for  i = 1:nrho    rhoab[i] = orba.P[i]*orbb.P[i] + orba.Q[i]*orbb.Q[i]   end
            wc = occa * occb
            for  k = 0:10
                gk = JAC.AngularMomentum.gk_ab(a,b);    if  gk == 0   continue    end
                for  i = 1:nrho    wx[i] = wx[i] + wc * gk * JAC.RadialIntegrals.Yk_ab(k, grid.r[i], rhoab, nrho, grid) * rhoab[i]   end
            end
        end
        # Weight factor
        for  a in basis.subshells
            occa = Basics.computeMeanSubshellOccupation(a, [level])
            orba = basis.orbitals[a]
            nrho = length(orba.P);      rhoa  = zeros(nrho)
            for  i = 1:nrho    rhoa[i] = orba.P[i]^2 + orba.Q[i]^2   end
            wc = occa
            for  i = 1:length(orba.P)    wg[i] = wg[i] + wc * rhoa[i] * grid.r[i]   end
        end
        #
        for  i = 1:npoints   if  wg[i] > 0.    wx[i] = wx[i]/wg[i]   else   wx[i] = 0.   end    end
        return( wx )
    end

    """
    `Basics.computeScfCoefficients(field::Basics.ALField, basis::Basis, subsh::Subshell)`
        ... to compute and sort out all angular coefficients that contribute to the self-consistent field of the subshell sh.
            For an average-level field, the energy functional is given by the sum of all diagonal matrix elements of the given
            basis; the function also adds all coefficients that give rise to the same XL_Coulomb interaction strengths.
            A tuple of two lists Tuple{Array{AngularTcoeff,1},Array{AngularVcoeff,1}} is returned.
    """
    function Basics.computeScfCoefficients(field::Basics.ALField, basis::Basis, subsh::Subshell)
        #  Select only angular coefficients for the Coulomb interaction
        function isCoulomb(L::Int64, sha::Subshell, shb::Subshell, shc::Subshell, shd::Subshell)
            la = Basics.subshell_l(sha);    ja2 = Basics.subshell_2j(sha)
            lb = Basics.subshell_l(shb);    jb2 = Basics.subshell_2j(shb)
            lc = Basics.subshell_l(shc);    jc2 = Basics.subshell_2j(shc)
            ld = Basics.subshell_l(shd);    jd2 = Basics.subshell_2j(shd)

            if  AngularMomentum.triangularDelta(ja2+1,jc2+1,L+L+1) * AngularMomentum.triangularDelta(jb2+1,jd2+1,L+L+1) == 0   ||
                rem(la+lc+L,2) == 1   ||   rem(lb+ld+L,2) == 1
                return( false )
            end
            return( true )
        end
        #
        Tcoeffs = JAC.AngularCoefficientsRatip2013.AngularTcoeff[];   allTcoeffs = JAC.AngularCoefficientsRatip2013.AngularTcoeff[]
        Vcoeffs = JAC.AngularCoefficientsRatip2013.AngularVcoeff[];   allVcoeffs = JAC.AngularCoefficientsRatip2013.AngularVcoeff[]
        # Compute all angular coefficients for the diagonal matrix elements
        for  csf in basis.csfs
            wa = Basics.compute("angular coefficients: e-e, Ratip2013", csf, csf)
            append!(allTcoeffs, wa[1])
            append!(allVcoeffs, wa[2])
        end
        ##x @show allTcoeffs
        ##x @show allVcoeffs
        # Determine a list reduced angular that just depend on subshell subsh
        tcoeffs = Tuple{Int64,Subshell,Subshell}[]
        vcoeffs = Tuple{Int64,Subshell,Subshell,Subshell,Subshell}[]
        for  coeff in allTcoeffs
            if   subsh == coeff.a
                wa = (coeff.nu, coeff.a, coeff.b)
                if      wa in tcoeffs
                else    push!(tcoeffs, wa)
                end
            end
        end
        for  coeff in allVcoeffs
            if   subsh == coeff.a   &&   isCoulomb(coeff.nu, coeff.a, coeff.b, coeff.c, coeff.d)
                wa = (coeff.nu, coeff.a, coeff.b, coeff.c, coeff.d)
                if      wa in vcoeffs
                else    push!(vcoeffs, wa)
                end
            end
            if   subsh == coeff.b   &&   isCoulomb(coeff.nu, coeff.b, coeff.a, coeff.d, coeff.c)
                wa = (coeff.nu, coeff.b, coeff.a, coeff.d, coeff.c)
                if      wa in vcoeffs
                else    push!(vcoeffs, wa)
                end
            end
        end
        # Now collect all contributions that belong to the selected reduced coefficients
        for  redCoeff in tcoeffs
            redT = 0.
            for  coeff  in  allTcoeffs
                wa = (coeff.nu, coeff.a, coeff.b)
                if  wa == redCoeff      redT = redT + coeff.T   end
            end
            push!( Tcoeffs, JAC.AngularCoefficientsRatip2013.AngularTcoeff( redCoeff[1], redCoeff[2], redCoeff[3], redT) )
        end
        for  redCoeff in vcoeffs
            redV = 0.
            for coeff  in  allVcoeffs
                ## V = coeff.V
                if      coeff.a == coeff.d   &&  coeff.a != coeff.b     V = coeff.V / 2.
                elseif  coeff.b == coeff.c   &&  coeff.a != coeff.b     V = coeff.V / 2.
                else                                                                V = coeff.V
                end
                wa = (coeff.nu, coeff.a, coeff.b, coeff.c, coeff.d);
                if  wa == redCoeff      redV = redV + V    end
                wa = (coeff.nu, coeff.b, coeff.a, coeff.d, coeff.c)
                if  wa == redCoeff      redV = redV + V    end
            end
            push!( Vcoeffs, JAC.AngularCoefficientsRatip2013.AngularVcoeff( redCoeff[1], redCoeff[2], redCoeff[3], redCoeff[4], redCoeff[5],  redV) )
        end
        ##x @show " "
        ##x @show Tcoeffs
        ##x @show Vcoeffs

        return( Tcoeffs, Vcoeffs )
    end


text:
"""
    FiniteDiffBackend <: AbstractDiffBackend

A type to specify / use differentiation backend based on FiniteDiff package.

# Constructor
    FiniteDiffBackend(method::Val{Symbol} = Val{:central})
"""
struct FiniteDiffBackend{TM<:Val} <: AbstractDiffBackend
    method::TM
end

FiniteDiffBackend() = FiniteDiffBackend(Val(:central))

function _derivative(f, p, ::FiniteDiffBackend{Method}) where {Method}
    return FiniteDiff.finite_difference_derivative(f, p, Method)
end

function _gradient(f, p, ::FiniteDiffBackend{Method}) where {Method}
    return FiniteDiff.finite_difference_gradient(f, p, Method)
end

function _gradient!(f, X, p, ::FiniteDiffBackend{Method}) where {Method}
    return FiniteDiff.finite_difference_gradient!(X, f, p, Method)
end

function _jacobian(f, p, ::FiniteDiffBackend{Method}) where {Method}
    return FiniteDiff.finite_difference_jacobian(f, p, Method)
end

function _jacobian!(f, X, p, ::FiniteDiffBackend{Method}) where {Method}
    return FiniteDiff.finite_difference_jacobian!(X, f, p, Method)
end

if default_differential_backend() === NoneDiffBackend()
    set_default_differential_backend!(FiniteDiffBackend())
end

ext: jl


text:
# Export
season = ""
param = "gas_26,45"
session_time = Dates.now()
session = Dates.format(session_time, "yyyy-mm-dd__HH-MM-SS")

# Create session folder for output
mkdir("output/" * session)

# ##########
# # CM info
# output_CM_info = DataFrame(parameter=String[],value=Float64[])
# push!(output_CM_info, ["Total system costs" JuMP.objective_value(CM_mod)])
# push!(output_CM_info, ["P_up" sum(CM_ΔP_up_opt)])
# push!(output_CM_info, ["P_dn" sum(CM_ΔP_dn_opt)])
# push!(output_CM_info, ["P_R_up" sum(CM_ΔP_R_up_opt)])
# push!(output_CM_info, ["P_R_dn" sum(CM_ΔP_R_dn_opt)])
#
# electrondisplay(output_CM_info)
# CSV.write("output/CM_info_" * season * param * ".csv", output_CM_info)
#
# # CM PtG info
# output_CM_PtG_info = DataFrame(parameter=String[],value=Float64[])
# push!(output_CM_PtG_info, ["Total system costs" JuMP.objective_value(CM_PtG_mod)])
# push!(output_CM_PtG_info, ["P_up" sum(CM_PtG_ΔP_up_opt)])
# push!(output_CM_PtG_info, ["P_dn" sum(CM_PtG_ΔP_dn_opt)])
# push!(output_CM_PtG_info, ["P_R_up" sum(CM_PtG_ΔP_R_up_opt)])
# push!(output_CM_PtG_info, ["P_R_dn" sum(CM_PtG_ΔP_R_dn_opt)])
# push!(output_CM_PtG_info, ["P_syn" sum(CM_PtG_P_syn_opt)])
# push!(output_CM_PtG_info, ["D_PtG" sum(CM_PtG_D_PtG_opt)])
# push!(output_CM_PtG_info, ["L_syn" sum(CM_PtG_L_syn_opt)])
#
# electrondisplay(output_CM_PtG_info)
# CSV.write("output/CM_PtG_info_" * season * param *".csv", output_CM_PtG_info)

function has(x, y)
	x == y
end

function getKeyVector(df::Dict, value)
	key_value = [k for (k,v) in df if has(v, value)]
	return key_value
end

###################
# ED ##############
###################

# Dispatchable power plants
output_ED_P_opt = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[])
for t in time
	for g in pp.unit
		val = ED_P_opt[t, g]
		fuel =  pp.fuel[g]
		push!(output_ED_P_opt, [t g fuel val])
	end
end
CSV.write("output/" * session * "/" * "ED_P.csv", output_ED_P_opt)

# RES generation
output_ED_P_R_opt = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[])
for t in time
	for r in res.unit
		val = ED_P_R_opt[t, r]
		tech =  res.fuel[r]
		push!(output_ED_P_R_opt, [t r tech val])
	end
end
CSV.write("output/" * session * "/" * "ED_P_R.csv", output_ED_P_R_opt)

# Maximum available RES infeed from input data
output_MaxRES = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[],capacity=Float64[])
for t in time
	for r in res.unit
		output = res.generation[r][t]
		capacity = res.pmax[r]
		tech =  res.fuel[r]
		push!(output_MaxRES, [t r tech output capacity])
	end
end
CSV.write("output/" * session * "/" * "MaxRES.csv", output_MaxRES)

# Load
output_load = DataFrame(time=Int64[],load=Float64[])
for t in time
	load = nodes.systemload[t]
	push!(output_load, [t load])
end
CSV.write("output/" * session * "/" * "load.csv", output_load)

# Prices
output_ED_price = DataFrame(time=Int64[],price=Float64[])
for t in time
	price =  ED_price[t]
	push!(output_ED_price, [t price])
end
CSV.write("output/" * session * "/" * "ED_price.csv", output_ED_price)

# Must-run
output_pmin = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[])
for t in time
	for g in pp.unit
		val = pp.pmin[g]
		fuel =  pp.fuel[g]
		push!(output_pmin, [t g fuel val])
	end
end
CSV.write("output/" * session * "/" * "pmin.csv", output_pmin)

###################
# CM ##############
###################

# Dispatchable power plants
output_CM_ΔP_up_opt = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[])
for t in time
	for g in pp.unit
		val = CM_ΔP_up_opt[t, g]
		fuel =  pp.fuel[g]
		push!(output_CM_ΔP_up_opt, [t g fuel val])
	end
end
CSV.write("output/" * session * "/" * "CM_P_up.csv", output_CM_ΔP_up_opt)

output_CM_ΔP_dn_opt = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[])
for t in time
	for g in pp.unit
		val = -1*CM_ΔP_dn_opt[t, g]
		fuel =  pp.fuel[g]
		push!(output_CM_ΔP_dn_opt, [t g fuel val])
	end
end
CSV.write("output/" * session * "/" * "CM_P_dn.csv", output_CM_ΔP_dn_opt)

# RES generation
output_CM_ΔP_R_up_opt = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[])
for t in time
	for r in res.unit
		val = CM_ΔP_R_up_opt[t, r]
		tech =  res.fuel[r]
		push!(output_CM_ΔP_R_up_opt, [t r tech val])
	end
end
CSV.write("output/" * session * "/" * "CM_P_R_up.csv", output_CM_ΔP_R_up_opt)

output_CM_ΔP_R_dn_opt = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[])
for t in time
	for r in res.unit
		val = -1*CM_ΔP_R_dn_opt[t, r]
		tech =  res.fuel[r]
		push!(output_CM_ΔP_R_dn_opt, [t r tech val])
	end
end
CSV.write("output/" * session * "/" * "CM_P_R_dn.csv", output_CM_ΔP_R_dn_opt)

###################
# CM PtG ##########
###################

# Dispatchable power plants
output_CM_PtG_ΔP_up_opt = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[])
for t in time
	for g in pp.unit
		val = CM_PtG_ΔP_up_opt[t, g]
		fuel =  pp.fuel[g]
		push!(output_CM_PtG_ΔP_up_opt, [t g fuel val])
	end
end
CSV.write("output/" * session * "/" * "CM_PtG_P_up.csv", output_CM_PtG_ΔP_up_opt)

output_CM_PtG_ΔP_dn_opt = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[])
for t in time
	for g in pp.unit
		val = -1*CM_PtG_ΔP_dn_opt[t, g]
		fuel =  pp.fuel[g]
		push!(output_CM_PtG_ΔP_dn_opt, [t g fuel val])
	end
end
CSV.write("output/" * session * "/" * "CM_PtG_P_dn.csv", output_CM_PtG_ΔP_dn_opt)

# RES generation
output_CM_PtG_ΔP_R_up_opt = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[])
for t in time
	for r in res.unit
		val = CM_PtG_ΔP_R_up_opt[t, r]
		tech =  res.fuel[r]
		push!(output_CM_PtG_ΔP_R_up_opt, [t r tech val])
	end
end
CSV.write("output/" * session * "/" * "CM_PtG_P_R_up.csv", output_CM_PtG_ΔP_R_up_opt)

output_CM_PtG_ΔP_R_dn_opt = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[])
for t in time
	for r in res.unit
		val = -1*CM_PtG_ΔP_R_dn_opt[t, r]
		tech =  res.fuel[r]
		push!(output_CM_PtG_ΔP_R_dn_opt, [t r tech val])
	end
end
CSV.write("output/" * session * "/" * "CM_PtG_P_R_dn.csv", output_CM_PtG_ΔP_R_dn_opt)

# Electricity demand from PtG units
output_CM_PtG_D_PtG_opt = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[])
for t in time
	for r in res.unit
		val = CM_PtG_D_PtG_opt[t, r]
		fuel = "PtG"
		push!(output_CM_PtG_D_PtG_opt, [t r fuel val])
	end
end
CSV.write("output/" * session * "/" * "CM_PtG_D_PtG.csv", output_CM_PtG_D_PtG_opt)

# Electricity generation from synthetic methane
output_CM_PtG_P_syn_opt = DataFrame(time=Int64[],unit=Int64[],fuel=String[],output=Float64[])
for t in time
	for e in getKeyVector(pp.fuel, "NaturalGas")
		val = CM_PtG_P_syn_opt[t, e]
		fuel = "SNG"
		push!(output_CM_PtG_P_syn_opt, [t e fuel val])
	end
end
CSV.write("output/" * session * "/" * "CM_PtG_P_syn.csv", output_CM_PtG_P_syn_opt)

# Synthetic methane storage level
output_CM_PtG_L_syn_opt = DataFrame(time=Int64[],storage=Float64[])
for t in time
	storage = CM_PtG_L_syn_opt[t]
	push!(output_CM_PtG_L_syn_opt, [t storage])
end
CSV.write("output/" * session * "/" * "CM_PtG_L_syn.csv", output_CM_PtG_L_syn_opt)

# System costs
output_system_cost = DataFrame(model=String[],system_cost=Float64[])
push!(output_system_cost, ["DCOPF" JuMP.objective_value(DCOPF_mod)])
push!(output_system_cost, ["Economic Dispatch" JuMP.objective_value(ED_mod)])
push!(output_system_cost, ["Congestion Management" JuMP.objective_value(CM_mod)])
push!(output_system_cost, ["Congestion Management (PtG)" JuMP.objective_value(CM_PtG_mod)])
CSV.write("output/" * session * "/" * "system_cost.csv", output_system_cost)

# ### Electron display
# electrondisplay(output_system_cost)
#
# # ED
# electrondisplay(output_ED_P_opt)
# electrondisplay(output_ED_P_R_opt)
# electrondisplay(output_MaxRES)
# electrondisplay(output_load)
# electrondisplay(output_ED_price)
#
# # CM
# electrondisplay(output_CM_ΔP_up_opt)
# electrondisplay(output_CM_ΔP_dn_opt)
# electrondisplay(output_CM_ΔP_R_up_opt)
# electrondisplay(output_CM_ΔP_R_dn_opt)
#
# # CM + PtG
# electrondisplay(output_CM_PtG_ΔP_up_opt)
# electrondisplay(output_CM_PtG_ΔP_dn_opt)
# electrondisplay(output_CM_PtG_ΔP_R_up_opt)
# electrondisplay(output_CM_PtG_ΔP_R_dn_opt)
# electrondisplay(output_CM_PtG_D_PtG_opt)
# electrondisplay(output_CM_PtG_P_syn_opt)
# electrondisplay(output_CM_PtG_L_syn_opt)


text:
let
    global af_lib
    if !isdefined(:af_lib)
        lib = is_unix()  ? begin "libaf" end : is_windows() ? "af" : ""
        try
            Libdl.dlopen(lib)
        catch
            error("ArrayFire library not found")
        end
        @eval const af_lib = $lib
    end
end


text:
using BinaryBuilder

name = "kubectl"
version = v"1.20.0"

# Links to downloads can found at: https://kubernetes.io/docs/tasks/tools/
sources = [
    FileSource("https://dl.k8s.io/release/v$(version)/bin/linux/amd64/kubectl", "a5895007f331f08d2e082eb12458764949559f30bcc5beae26c38f3e2724262c"; filename="x86_64-linux-gnu-kubectl"),
    FileSource("https://dl.k8s.io/release/v$(version)/bin/linux/arm64/kubectl", "25e4465870c99167e6c466623ed8f05a1d20fbcb48cab6688109389b52d87623"; filename="aarch64-linux-gnu-kubectl"),
    FileSource("https://dl.k8s.io/release/v$(version)/bin/darwin/amd64/kubectl", "82046a4abb056005edec097a42cc3bb55d1edd562d6f6f38c07318603fcd9fca"; filename="x86_64-apple-darwin14-kubectl"),
    FileSource("https://dl.k8s.io/release/v$(version)/bin/windows/amd64/kubectl.exe", "ee7be8e93349fb0fd1db7f5cdb5985f5698cef69b7b7be012fc0e6bed06b254d"; filename="x86_64-w64-mingw32-kubectl"),
    FileSource("https://raw.githubusercontent.com/kubernetes/kubernetes/v$(version)/LICENSE", "cfc7749b96f63bd31c3c42b5c471bf756814053e847c10f3eb003417bc523d30")
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/
mkdir -p ${bindir}

install_license LICENSE

mv "${target}-kubectl" "${bindir}/kubectl${exeext}"
chmod 755 "${bindir}/kubectl${exeext}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "windows"),
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("kubectl", :kubectl),
]

dependencies = Dependency[]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
export subtract_PBC, subtract_PBC!

"""
    subtract_PBC( x1::Real, x0::Real; axis_size )

Subtraction defined on a _ring_ with periodic boundary conditions of size `axis_size`.

# Examples
```jldocstest
julia> subtract_PBC( 2., 1.; axis_size = 10 )
1.0

julia> subtract_PBC( 7., 1.; axis_size = 10 )
-4.0

julia> subtract_PBC( 1., 8.; axis_size = 10 )
3.0
```
"""
function subtract_PBC( x1::Real, x0::Real; axis_size )
    value = x1 - x0
    while value > axis_size/2
        value -= axis_size
    end
    while value < -axis_size/2
        value += axis_size
    end
    return value
end

"""
    subtract_PBC!( A::Vector2D, B::Vector2D; Lx, Ly = Lx )

In-place `Vector2D` subtraction with periodic boundary conditions.

# Examples
```jldocstest
julia> A = Vector2D(1., 2.)
Vector2D{Float64}([1.0, 2.0])

julia> B = Vector2D(7., 9.)
Vector2D{Float64}([7.0, 9.0])

julia> subtract_PBC!(A, B; Lx = 10)

julia> A
Vector2D{Float64}([4.0, 3.0])
```
"""
function subtract_PBC!( A::Vector2D, B::Vector2D; Lx, Ly = Lx )
    broadcast!( (a, b, ax) -> subtract_PBC(a, b; axis_size=ax), A.vec, A.vec, B.vec, (Lx, Ly) )
    return nothing
end