text:
"""
`module  JAC.Basics`
    ... a submodel of JAC that contains many basic types/struct that are specific to the JAC module; this module also defines
        a number of basic functions/methods that are later extended. We here provide proper docstrings for all abstract and
        conrete types (struct) in order to allow the user to easy understand and access the individual fields of each type/struct
        definition.
"""
module Basics

    using Printf

    export  add, analyze, AngularJ, AngularJ64, AngularM, AngularM64, Auger,
            compute, diagonalize,
            estimate, EmMultipole, E1, M1, E2, M2, E3, M3, E4, M4, EmProperty, ExpStokes, EmStokes,
            generate,
            ##x HalfInt, HalfInteger,
            interpolate, integrate,
            LevelSymmetry, LevelKey,
            modify,
            oplus,
            Parity, perform, provide,
            RadialMesh, Radiative,
            SolidAngle, Shell, Subshell, subshell_2j,
            tabulate, tools, TensorComp,
            UseCoulomb, UseBabushkin

    ##x include("inc-halfintegers.jl")

    abstract type  AngularJ  end

    """
    `struct  AngularJ64 <: AngularJ`  ... defines a type for angular momenta j.

        + num  ::Int64              ... numerator
        + den  ::Int64              ... denominator, must be 1 or 2
    """
    struct  AngularJ64 <: AngularJ
        num      ::Int64
        den      ::Int64
    end

    """
    `Basics.AngularJ64(m::Integer)`  ... constructor for a given integer (numerator).
    """
    function AngularJ64(j::Integer)
        j < 0   &&   error("j must be positive.")
        AngularJ64(j, 1)
    end

    """
    `Basics.AngularJ64(rational::Rational{Int64})`  ... constructor for a given  Rational{Int64}.
    """
    function AngularJ64(rational::Rational{Int64})
        !(rational.den in [1,2])   &&   error("Denominator must be 1 or 2.")
        rational.num < 0           &&   error("j must be positive.")
        AngularJ64(rational.num, rational.den)
    end

    # `Base.show(io::IO, j::AngularJ64)`  ... prepares a proper printout of the variable j::AngularJ64.
    function Base.show(io::IO, j::AngularJ64)
        if      j.den == 1    print(io, j.num)
        elseif  j.den == 2    print(io, j.num, "/2")
        else    error("stop a")
        end
    end

    """
    `Basics.oplus(ja::AngularJ64, jb::AngularJ64)`
        ... adds the angular momenta ja `oplus` jb and returns a list::Array{AngularJ64,1} of j-valuescin the interval
            |ja - jb| <= j <= ja + jb.
    """
    function  oplus(ja::AngularJ64, jb::AngularJ64)
        if  ja.den == 1   ja2 = 2ja.num   else   ja2 = ja.num   end
        if  jb.den == 1   jb2 = 2jb.num   else   jb2 = jb.num   end
        jList = AngularJ64[]
        for  j = abs(ja2 - jb2):2:ja2+jb2    push!(jList, AngularJ64(j//2) )    end
        return( jList )
    end

    ##x oplus(ja, jb) = oplus(HalfInt(ja), HalfInt(jb))
    ##x oplus(ja::Union{HalfInt,Integer}, jb::Union{HalfInt,Integer}) = abs(ja-jb):ja+jb

    ##x const ⊕ = oplus

    """
    `Basics.projections(ja::AngularJ64)`
        ... returns all allowed projections of a given angular momenta ja as a list::Array{AngularM64,1} of m-values, i.e. -ja, -ja+1, ..., j.
    """
    function  projections(ja::AngularJ64)
        if  ja.den == 1   ja2 = 2ja.num   else   ja2 = ja.num   end
        mList = AngularM64[]
        for  m = -ja2:2:ja2    push!(mList, AngularM64(m//2) )    end
        return( mList )
    end

    ##x projections(j) = projections(HalfInt(j))
    ##x projections(j::Union{HalfInt,Integer}) = j≥0 ? (-j:j) : throw(DomainError(j, "angular momentum j must be non-negative."))

    abstract type  AngularM  end

    """
    `struct  Basics.AngularM64 <: AngularM`  ... defines a type for angular momentum projections m.

        + num  ::Int64              ... numerator
        + den  ::Int64              ... denominator, must be 1 or 2
    """
    struct  AngularM64 <: AngularM
        num      ::Int64
        den      ::Int64
    end

    """
    `Basics.AngularM64(m::Integer)`  ... constructor for a given integer (numerator).
    """
    function AngularM64(m::Integer)
        AngularM64(m, 1)
    end

    """
    `Basics.AngularM64(m::Integer, j::AngularJ64)`
        ... constructor for a given integer (numerator) that is consistent with a given j-value.
    """
    function AngularM64(m::Integer, j::AngularJ64)
        !(j.den == 1)      &&  error("m must be integer for j = $(j).")
        j.num < abs(m)   &&  error("abs(m) must be <= j = $(j).")
        AngularM64(m, 1)
    end

    """
    `Basics.AngularM64(rational::Rational{Int64})`  ... constructor for a given  Rational{Int64}.
    """
    function AngularM64(rational::Rational{Int64})
        !(rational.den in [1,2])  &&  error("Denominator must be 1 or 2.")
        AngularM64(rational.num, rational.den)
    end

    """
    `Basics.AngularM64(rational::Rational{Int64}, j::AngularJ64)`
        ... constructor for a given  Rational{Int64} that is consistent with a given m-value.
    """
    function AngularM64(rational::Rational{Int64}, j::AngularJ64)
        !(rational.den in [1,2])      &&   error("Denominator must be 1 or 2.")
        !(j.den == rational.den)      &&   error("j,m must be both integer or half-integer.")
        j.num < abs(rational.num)   &&   error("abs(m) must be <= j = $(j).")
        AngularM64(rational.num, rational.den)
    end

    """
    `Basics.AngularM64(j::AngularJ64)`  ... constructor for a given j::AngularJ64 to define m = j.
    """
    function AngularM64(j::AngularJ64)
        AngularM64(j.num, j.den)
    end

    # `Base.show(io::IO, m::AngularM64)`  ... prepares a proper printout of the variable m::AngularM64.
    function Base.show(io::IO, m::AngularM64)
        if      m.den == 1    print(io, m.num)
        elseif  m.den == 2    print(io, m.num, "/2")
        else    error("stop a")
        end
    end

    """
    `Basics.add(ma::AngularM64, mb::AngularM64)`
        ... adds the projections of the angular momenta ma + mb and returns a mc::AngularM64.
    """
    function  add(ma::AngularM64, mb::AngularM64)
        if  ma.den == 1   ma2 = 2ma.num   else   ma2 = ma.num   end
        if  mb.den == 1   mb2 = 2mb.num   else   mb2 = mb.num   end
        return( AngularM64( (ma2+mb2)//2 ) )
    end

    # Conversion between HalfInt and AngularJ64, AngularM64

    twice(x::Union{AngularJ64,AngularM64}) = ifelse(isone(x.den), twice(x.num), x.num)
    twice(x) = x + x

    ##x AngularJ64(x::HalfInt) = isinteger(x) ? AngularJ64(Integer(x)) : AngularJ64(twice(x), 2)
    ##x AngularM64(x::HalfInt) = isinteger(x) ? AngularM64(Integer(x)) : AngularM64(twice(x), 2)

    """
    `struct  Basics.Eigen`
        ... defines a simple struct to communicate eigenvalues and eigenvectors if different diagonalization procedures are used.

        + values   ::Array{Float64,1}            ... List of eigenvalues.
        + vectors  ::Array{Vector{Float64},1}    ... List of eigenvectors.
    """
    struct  Eigen
        values     ::Array{Float64,1}
        vectors    ::Array{Vector{Float64},1}
    end

    """
    `Basics.isEqual(a::Eigen, b::Eigen, accuracy::Float64)`
        ... determines whether a == b with the requested accuracy; a value::Bool is returned. The test is made element-wise on
            abs(a.number - b.number) > accuracy.
    """
    function isEqual(a::Eigen, b::Eigen, accuracy::Float64)
        # Check the dimensions
        if  accuracy <= 0.                          error("Requires accuracy > 0.")                                    end
        if  size(a.values) != size(b.values)   ||   size(a.vectors) != size(b.vectors)                return( false )  end
        # Check the eigenvalues and vectors
        for  i = 1:length(a.values)
            if  abs(a.values[i] - b.values[i]) > accuracy                                             return( false )  end
            for  j = 1:length(a.vectors[i])   if  abs(a.vectors[i][j] - b.vectors[i][j]) > accuracy   return( false )  end    end
        end
        return( true )
    end

    """
    `struct  Basics.SolidAngle`  ... defines a type for a solid angle Omega = (theta, phi).

        + theta  ::Float64       ... polar angle theta with regard to the z-axis; 0 <= theta <= pi.
        + phi    ::Float64       ... azimuthal angle phi (with regard to the x-z plane in mathematical positive direction); 0 <= phi <= 2 pi.
    """
    struct  SolidAngle
        theta    ::Float64
        phi      ::Float64
    end

    # `Base.show(io::IO, Omega::SolidAngle)`  ... prepares a proper printout of the variable Omega::SolidAngle.
    function Base.show(io::IO, Omega::SolidAngle)
        sa = "Omega == (" * string(Omega.theta) * ", " * string(Omega.phi) * ")"
        print(io, sa)
    end

    """
    `@enum   Parity`  ... defines a enumeration for the allowed values of parity.

        + plus, minus               ... with obvious meaning
    """
    @enum   Parity    plus   minus

    """
    `Basics.Parity(sa::String)`  ... constructor for a given String.
    """
    function Parity(sa::String)
        if       sa == "+"    wa = plus
        elseif   sa == "-"    wa = minus
        else     error("stop a")
        end

        return( wa )
    end

    """
    `Basics.invertParity(p::Parity)`  ... inverts the given parity plus <--> minus.
    """
    function invertParity(p::Parity)
        if       p == plus    return( minus )
        elseif   p == minus   return( plus )
        else     error("stop a")
        end
    end

    #============================================
    # `Base.show(io::IO, p::Parity)`  ... prepares a proper printout of the variable p::Parity.
    function Base.show(io::IO, p::Parity)
        print(io, string(p) )
    end  ==================================#

    # `Base.string(p::Parity)`  ... provides a proper printout of the variable p::Parity.
    function Base.string(p::Parity)
        if      p == plus     return( "+" )
        elseif  p == minus    return( "-" )
        else    error("stop a")
        end
    end

    """
    `struct  Basics.LevelSymmetry`  ... defines a struct for defining the overall J^P symmetry of a level.

        + J          ::AngularJ64  ... total angular momentum of a level
        + parity     ::Parity      ... total parity of the level
    """
    struct  LevelSymmetry
        J            ::AngularJ64
        parity       ::Parity
    end

    """
    `Basics.LevelSymmetry(rational::Rational{Int64}, parity::Parity)`  ... constructor for a given (Rational,Parity).
    """
    function  LevelSymmetry(rational::Rational{Int64}, parity::Parity)
        !(rational.den in [1,2])      &&   error("Denominator must be 1 or 2.")
        LevelSymmetry( AngularJ64(rational.num, rational.den), parity )
    end

    """
    `Basics.LevelSymmetry(i::Int64, parity::Parity)`  ... constructor for a given (Int64,Parity).
    """
    function  LevelSymmetry(i::Int64, parity::Parity)
        LevelSymmetry( AngularJ64(i), parity )
    end

    """
    `Basics.LevelSymmetry(rational::Rational{Int64}, sa::String)`  ... constructor for a given (Rational,String).
    """
    function  LevelSymmetry(rational::Rational{Int64},sa::String)
        !(rational.den in [1,2])      &&   error("Denominator must be 1 or 2.")
        LevelSymmetry( AngularJ64(rational.num, rational.den), Parity(sa) )
    end

    """
    `Basics.LevelSymmetry(i::Int64, sa::String)`  ... constructor for a given (Int64,String).
    """
    function  LevelSymmetry(i::Int64,sa::String)
        LevelSymmetry( AngularJ64(i), Parity(sa) )
    end

    # `Base.show(io::IO, sym::LevelSymmetry)`  ... prepares a proper printout of the variable sym::LevelSymmetry.
    function Base.show(io::IO, sym::LevelSymmetry)
        print(io, string(sym) )
    end

    # `Base.string(sym::LevelSymmetry)`  ... provides a proper printout of the variable sym::LevelSymmetry.
    function Base.string(sym::LevelSymmetry)
        if      sym.parity == plus     return( "$(sym.J) +" )
        elseif  sym.parity == minus    return( "$(sym.J) -" )
        else    error("stop a")
        end
    end

    """
    `struct  Basics.LevelKey`  ... defines a struct for identifying a level by its symmetry, energy, etc.

        + sym          ::LevelSymmetry    ... Level symmetry.
        + index        ::Int64            ... Index of this level in its original multiplet, or 0.
        + energy       ::Float64          ... Total energy.
        + relativeOcc  ::Float64          ... Relative occupation of this level, if part of some multiplet/cascade.
    """
    struct  LevelKey
        sym            ::LevelSymmetry
        index          ::Int64
        energy         ::Float64
        relativeOcc    ::Float64
    end

    """
    `Basics.LevelKey()`  ... constructor for an empty instance of a LevelKey.
    """
    function  LevelKey()
        LevelKey( LevelSymmetry(0,Basics.plus), 0, 0., 0.)
    end

    # `Base.show(io::IO, key::LevelKey)`  ... prepares a proper printout of the variable key::LevelKey.
    function Base.show(io::IO, key::LevelKey)
        print(io, string(key) )
    end

    # `Base.string(key::LevelKey)`  ... provides a proper printout of the variable key::LevelKey.
    function Base.string(key::LevelKey)
        return( "LevelKey[$(key.sym) ($(key.index)), en=$(key.energy), occ=$(key.relativeOcc)]")
    end

    """
    `struct  Basics.LineKey`
        ... defines a struct for identifying a line by the keys of the initial and final level, its transition energy,
            decay strength, etc.

        + iLevelKey    ::LevelKey         ... Key of the initial level.
        + fLevelKey    ::LevelKey         ... Key of the final level.
        + energy       ::Float64          ... Total energy.
        + strength     ::Float64          ... Strength of the given line.
    """
    struct  LineKey
        iLevelKey      ::LevelKey
        fLevelKey      ::LevelKey
        energy         ::Float64
        strength       ::Float64
    end

    """
    `Basics.LineKey()`  ... constructor for an empty instance of a LineKey.
    """
    function  LineKey()
        LineKey( LevelKey(), LevelKey(), 0., 0.)
    end

    # `Base.show(io::IO, key::LineKey)`  ... prepares a proper printout of the variable key::LineKey.
    function Base.show(io::IO, key::LineKey)
        print(io, string(key) )
    end

    # `Base.string(key::LineKey)`  ... provides a proper printout of the variable key::LevelKey.
    function Base.string(key::LineKey)
        return( "LineKey[i=$(key.iLevelKey) --> f=$(key.fLevelKey), en=$(key.energy), strength=$(key.strength)]")
    end

    """
    `struct  Basics.Shell`  ... defines a enumeration for the allowed values of a non-relativistic shell.

        + n     ::Int64  ... principal quantum number
        + l     ::Int64  ... orbital angular quantum number
    """
    struct  Shell
        n       ::Int64
        l       ::Int64
    end

    """
    `Basics.Shell(sa::Union{String,SubString{String}})`  ... constructor for a given String, such as 1s, 2s, 2p, ... .
    """
    function Shell(sa::Union{String,SubString{String}})
        wa = strip( Base.convert(String, sa) )
        n  = parse(Int64, wa[1:end-1]);    l = shellNotation( string(wa[end]) )
        Shell(n, l)
    end

    # `Base.show(io::IO, sh::Shell)`  ... prepares a proper printout of the variable sh::Shell.
    function Base.show(io::IO, sh::Shell)
        print(io, string(sh) )
    end

    # `Base.string(sh::Shell)`  ... provides a proper printout of the variable sh::Shell.
    function Base.string(sh::Shell)
        return( string(sh.n) * shellNotation(sh.l) )
    end

    """
    `Basics.shellNotation(l::Int64)`  ... returns the corresponding spectroscopy letter for a given orbital angular quantum number l.
    """
    function shellNotation(l::Int64)
        !(0 <= l <= 25)      &&   error("Orbital QN 0 <= l <= 25; l = $l")
        wa = [ "s", "p", "d", "f", "g", "h", "i", "j", "k", "l", "m", "n", "o", "q", "r", "t", "u", "v", "w", "x", "y", "z"]
        return( wa[l+1] )
    end

    """
    `Basics.shellNotation(sa::String)`  ... returns for a given spectroscopy letter the orbital angular quantum number l.
    """
    function shellNotation(sa::String)
        wa = [ "s", "p", "d", "f", "g", "h", "i", "j", "k", "l", "m", "n", "o", "q", "r", "t", "u", "v", "w", "x", "y", "z"]
        wb = findall(in([sa]), wa)
        !(length(wb)==1)      &&   error("Unrecognized orbital letter = $sa")
        return( wb[1]-1)
    end

    """
    `Basics.shellSplitIntoSubshells(sh::Shell)'  ... to convert a non-relativistic shell into a list of Subshell[..]
    """
    function shellSplitIntoSubshells(sh::Shell)
        subshells = Subshell[]
        if     sh.l == 0   push!(subshells, Subshell(sh.n, -1) )
        else   j2 = 2sh.l - 1;    kappa = Int64(  (j2+1)/2 );   push!(subshells, Subshell(sh.n, kappa) )
            j2 = 2sh.l + 1;    kappa = Int64( -(j2+1)/2 );   push!(subshells, Subshell(sh.n, kappa) )
        end
        return( subshells )
    end

    """
    `Basics.shellSplitOccupation(sh::Shell, occ::Int64)'
        ... to split the occupation of a shell into proper subshell occupations; an  Array{Dict{Basics.Subshell,Int64},1} is returned.
    """
    function shellSplitOccupation(sh::Shell, occ::Int64)

        subshells = shellSplitIntoSubshells(sh)
        wa        = Dict{Subshell,Int64}[]

        for  k1 in 0:occ
            for  k2 in 0:occ
                if  k1 + k2 == occ
                    if       length( subshells ) == 1
                        !(0 <= occ <= Basics.subshell_2j( subshells[1] ) + 1)   &&    error("Wrong pair of subshell = $(subshells[1]) and occupation = $occ")
                        wb = Dict( subshells[1] => occ )
                        push!(wa, wb)
                        return( wa )

                        elseif  length( subshells ) == 2
                            j2_1 = Basics.subshell_2j(subshells[1])
                        j2_2 = Basics.subshell_2j(subshells[2])
                        if  k1 <= j2_1 + 1   &&    k2 <= j2_2 + 1
                        wb = Dict( subshells[1] => k1, subshells[2] => k2 )
                        # Delete contributions with zero occupation
                        ## if  k1 == 0    delete!(wb, subshells[1])   end
                        ## if  k2 == 0    delete!(wb, subshells[2])   end
                        push!(wa, wb)
                        end
                    else    error("stop a")
                    end
                end
            end
        end

        return( wa )
    end

    """
    `struct  Basics.Subshell`  ... defines a type for the allowed values of a relativistic subhell.

        + n        ::Int64  ... principal quantum number
        + kappa    ::Int64  ... relativistic angular quantum number
    """
    struct  Subshell
        n          ::Int64
        kappa      ::Int64
    end

    """
    `Basics.Subshell(sa::String)`  ... constructor for a given String, such as 1s_1/2, 2s_1/2, 2p_3/2, ... .
    """
    function Subshell(sa::String)
        wa = strip( sa );
        wb = findnext("_", wa, 1);    if  wb == nothing   error("No underscore in sa = $sa")           else   wb = wb[1]   end
        wc = findnext("/", wa, wb);   if  wc == nothing   error("No / in wa[] = $(wa[wb[1]+1:end])")   else   wc = wc[1]   end
        wd = string( wa[wb+1:wc-1] )
        n  = parse(Int64, wa[1:wb-2]);                  l  = shellNotation( string(wa[ wb-1 ]) )
        j2 = parse(Int64, wd );                         !(string(wa[end]) == "2")    &&  error("Unrecognized subshell string sa = $sa")
        if  l - j2/2 > 0     kappa = Int64( (j2+1)/2 )    else    kappa = Int64( -(j2+1)/2 )   end
        Subshell(n, kappa)
    end

    # `Base.show(io::IO, sh::Subshell)`  ... prepares a proper printout of the variable sh::Subshell.
    function Base.show(io::IO, sh::Subshell)
        print(io, string(sh) )
    end

    # `Base.string(sh::Subshell)`  ... provides a proper printout of the variable sh::Subshell.
    function Base.string(sh::Subshell)
        sa = string(sh.n) * shellNotation( Basics.subshell_l(sh) ) * "_" * string( Basics.subshell_j(sh) )
        return( sa )
    end

    """
    `Basics.subshellGrasp(sa::Stringl)`  ... returns a sh::Subshell if sa denotes a (relativistic) Subshell in Grasp-notation.
    """
    function subshellGrasp(sa::String)
            if   sa == "1s"   return( Subshell("1s_1/2")  )
        elseif   sa == "2s"   return( Subshell("2s_1/2")  )   elseif   sa == "2p-"   return( Subshell("2p_1/2") )
        elseif   sa == "2p"   return( Subshell("2p_3/2")  )
        elseif   sa == "3s"   return( Subshell("3s_1/2")  )   elseif   sa == "3p-"   return( Subshell("3p_1/2") )
        elseif   sa == "3p"   return( Subshell("3p_3/2")  )   elseif   sa == "3d-"   return( Subshell("3d_3/2") )
        elseif   sa == "3d"   return( Subshell("3d_5/2")  )
        elseif   sa == "4s"   return( Subshell("4s_1/2")  )   elseif   sa == "4p-"   return( Subshell("4p_1/2") )
        elseif   sa == "4p"   return( Subshell("4p_3/2")  )   elseif   sa == "4d-"   return( Subshell("4d_3/2") )
        elseif   sa == "4d"   return( Subshell("4d_5/2")  )   elseif   sa == "4f-"   return( Subshell("4f_5/2") )
        elseif   sa == "4f"   return( Subshell("4f_7/2")  )
        elseif   sa == "5s"   return( Subshell("5s_1/2")  )   elseif   sa == "5p-"   return( Subshell("5p_1/2") )
        elseif   sa == "5p"   return( Subshell("5p_3/2")  )   elseif   sa == "5d-"   return( Subshell("5d_3/2") )
        elseif   sa == "5d"   return( Subshell("5d_5/2")  )   elseif   sa == "5f-"   return( Subshell("5f_5/2") )
        elseif   sa == "5f"   return( Subshell("5f_7/2")  )   elseif   sa == "5g-"   return( Subshell("5g_7/2") )
        elseif   sa == "5g"   return( Subshell("5g_9/2")  )
        elseif   sa == "6s"   return( Subshell("6s_1/2")  )   elseif   sa == "6p-"   return( Subshell("6p_1/2") )
        elseif   sa == "6p"   return( Subshell("6p_3/2")  )   elseif   sa == "6d-"   return( Subshell("6d_3/2") )
        elseif   sa == "6d"   return( Subshell("6d_5/2")  )   elseif   sa == "6f-"   return( Subshell("6f_5/2") )
        elseif   sa == "6f"   return( Subshell("6f_7/2")  )   elseif   sa == "6g-"   return( Subshell("6g_7/2") )
        elseif   sa == "6g"   return( Subshell("6g_9/2")  )
        elseif   sa == "7s"   return( Subshell("7s_1/2")  )   elseif   sa == "7p-"   return( Subshell("7p_1/2") )
        elseif   sa == "7p"   return( Subshell("7p_3/2")  )   elseif   sa == "7d-"   return( Subshell("7d_3/2") )
        else     error("Grasp notation ... not yet implemented; sa = $sa")
        end
    end

    """
    `Basics.subshell_l(sh::Subshell)`  ... returns the orbital angular quantum number l for a given sh::Subshell; an Int64 is returned
    """
    function subshell_l(sh::Subshell)
        if   sh.kappa < 0   l = abs(sh.kappa) -1   else   l = sh.kappa   end
        return( l )
    end

    """
    `Basics.subshell_j(sh::Subshell)`  ... returns the total angular quantum number j for a given sh::Subshell; an AngularJ64 is returned.
    """
    function subshell_j(sh::Subshell)
        return( AngularJ64(abs(sh.kappa) - 1//2) )
    end

    """
    `Basics.subshell_2j(sh::Subshell)`  ... returns the total angular quantum number j for a given sh::Subshell; an AngularJ64 is returned.
    """
    function subshell_2j(sh::Subshell)
        return( 2* abs(sh.kappa) - 1 )
    end

    """
    `Basics.subshellsFromClosedShellConfiguration("[Ne]")`
        ... to provide a list of (relativistic) subshells for the given closed-shell configuration.
    """
    function subshellsFromClosedShellConfiguration(sa::String)
        if       sa == "[He]"    wa = [ Subshell("1s_1/2")]
        elseif   sa == "[Ne]"    wa = [ Subshell("1s_1/2"), Subshell("2s_1/2"), Subshell("2p_1/2"), Subshell("2p_3/2")]
        elseif   sa == "[Mg]"    wa = [ Subshell("1s_1/2"), Subshell("2s_1/2"), Subshell("2p_1/2"), Subshell("2p_3/2"), Subshell("3s_1/2")]
        elseif   sa == "[Ar]"    wa = [ Subshell("1s_1/2"), Subshell("2s_1/2"), Subshell("2p_1/2"), Subshell("2p_3/2"),
                                                            Subshell("3s_1/2"), Subshell("3p_1/2"), Subshell("3p_3/2")]
        elseif   sa == "[Kr]"    wa = [ Subshell("1s_1/2"), Subshell("2s_1/2"), Subshell("2p_1/2"), Subshell("2p_3/2"),
                                                            Subshell("3s_1/2"), Subshell("3p_1/2"), Subshell("3p_3/2"),
                                                            Subshell("3d_3/2"), Subshell("3d_5/2"),
                                        Subshell("4s_1/2"), Subshell("4p_1/2"), Subshell("4p_3/2")]
        elseif   sa == "[Xe]"    wa = [ Subshell("1s_1/2"), Subshell("2s_1/2"), Subshell("2p_1/2"), Subshell("2p_3/2"),
                                                            Subshell("3s_1/2"), Subshell("3p_1/2"), Subshell("3p_3/2"),
                                                            Subshell("3d_3/2"), Subshell("3d_5/2"),
                                        Subshell("4s_1/2"), Subshell("4p_1/2"), Subshell("4p_3/2"),
                                                            Subshell("4d_3/2"), Subshell("4d_5/2"),
                                        Subshell("5s_1/2"), Subshell("5p_1/2"), Subshell("5p_3/2")]
        else
            error("Unsupported keystring = $sa ")
        end

        return( wa )
    end

    """
    `struct   Basics.SubshellStateR`
        ... defines a struct for the relativistic antisymmetric subshell states in the seniority scheme; this struct can be accessed
            only internally and, therefore, the only the standard constructor is supported.

        + subshell     ::Subshell        ... to refer to a particular subshell.
        + occ          ::Int64           ... occupation of the subshell
        + nu           ::Int64           ... seniority of the subshell state, often called nu
        + Jsub2        ::Int64           ... 2*Jsub
    """
    struct  SubshellStateR
        subshell       ::Subshell
        occ            ::Int64
        nu             ::Int64
        Jsub2          ::Int64
    end

    # `Base.show(io::IO, state::SubshellStateR)`  ... prepares a proper printout of the variable state::SubshellStateR.
    function Base.show(io::IO, state::SubshellStateR)
        print(io, string(state) )
    end

    # `Base.string(state::SubshellStateR)`  ... provides a proper printout of the variable state::SubshellStateR.
    function Base.string(state::SubshellStateR)
        sa = "SubshellState: [" * string(state.subshell) * "^$(state.occ), nu=$(state.nu), 2*J_sub=$(state.Jsub2)]"
    end

    """
    `Basics.subshellStateString(subshell::String, occ::Int64, seniority::Int64, Jsub::AngularJ64, X::AngularJ64)`
        ... to provide a string of a given subshell state in the form '[2p_1/2^occ]_(seniority, J_sub), X=Xo' ... .
    """
    function subshellStateString(subshell::String, occ::Int64, seniority::Int64, Jsub::AngularJ64, X::AngularJ64)
        sa = "[" * subshell * "^$occ]_($seniority, " * string(Jsub) * ") X=" * string(X)
        return( sa )
    end

    """
    `@enum   Basics.AtomicProcess`
        ... defines a enumeration for the atomic processes that are supported in the JAC program in the computation of transition
            arrays or whole excitation/decay cascades.

        + Auger         ... Auger transitions, i.e. single autoionization or the emission of a single free electron into the continuum.
        + AugerInPlasma ... Auger transitions but calculated for a specified plasma model.
        + Compton       ... Rayleigh-Compton scattering cross sections.
        + Coulex        ... Coulomb-excitation of target or projeticle electrons by fast, heavy ions.
        + Coulion       ... Coulomb-ionization of target or projeticle electrons by fast, heavy ions.
        + Dierec        ... di-electronic recombination, i.e. the dielectronic capture of a free electron and the subsequent emission of a photon.
        + MultiPhotonDE ... multi-photon excitation and decay rates, including 2-photon, etc. processes.
        + MultiPI       ... multi-photon (single-electron) ionization.
        + MultiPDI      ... multi-photon (single-electron) double ionization.
        + Photo         ... Photoionization processes, i.e. the emission of a single free electron into the continuum due to an external light field.
        + PhotoExc      ... Photoexcitation rates.
        + PhotoExcFluor ... photoexcitation fluorescence rates and cross sections.
        + PhotoExcAuto  ... photoexcitation autoionization cross sections and collision strengths.
        + PhotoInPlasma ... Photoionization processes but calculated for a specified plasma model.
        + PhotoIonFluor ... photoionization fluorescence rates and cross sections.
        + PhotoIonAuto  ... photoionization autoionization cross sections and collision strengths.
        + Radiative    ... Radiative (multipole) transitions between bound-state levels of the same charge state.
        + Rec           ... radiative electron capture, i.e. the capture of a free electron with the simultaneous emission of a photon.
        + Eimex         ... electron-impact excitation cross sections and collision strengths.
        + RAuger        ... Radiative Auger rates.
    """
    @enum   AtomicProcess  NoProcess  Auger  AugerInPlasma  Compton  Coulex  Dierec  Eimex  ImpactExcAuto  InternalConv  MultiPhotonDE  MultiPI  MultiPDI=
                    20  Photo  PhotoExc  PhotoExcAuto  PhotoExcFluor  PhotoInPlasma  PhotoIonAuto  PhotoIonFluor   Radiative  RAuger=
                    40  Rec  PairA1P  Coulion
    export  AtomicProcess, NoProcess, Auger, AugerInPlasma, Compton, Coulex, Dierec, Eimex, ImpactExcAuto, InternalConv, MultiPhotonDE, MultiPI, MultiPDI,
                        Photo, PhotoExc, PhotoExcAuto, PhotoExcFluor, PhotoInPlasma, PhotoIonAuto, PhotoIonFluor,  Radiative, RAuger,
                        Rec, PairA1P, Coulion

    """
    `Basics.AtomicProcess(sa::String)`  ... constructor for a given String.
    """
    function AtomicProcess(sa::String)
        if       sa in [ "none", "NoProcess"]                       wa = NoProcess
        elseif   sa in [ "Auger"]                                   wa = Auger
        elseif   sa in [ "Auger in plasma"]                         wa = AugerInPlasma
        elseif   sa in [ "Compton", "Rayleigh"]                     wa = Compton
        elseif   sa in [ "Coulex", "Coulomb excitation"]            wa = Coulex
        elseif   sa in [ "Coulion", "Coulomb ionization"]           wa = Coulion
        elseif   sa in [ "dierec", "Dierec"]                        wa = Dierec
        elseif   sa in [ "eimex", "Eimex"]                          wa = Eimex
        elseif   sa in [ "internal conversion"]                     wa = InternalConv
        elseif   sa in [ "multi-DE"]                                wa = MultiPhotonDE
        elseif   sa in [ "multi-PI"]                                wa = MultiPI
        elseif   sa in [ "multi-PDI"]                               wa = MultiPDI
        elseif   sa in [ "photo", "Photo"]                          wa = Photo
        elseif   sa in [ "photo in plasma"]                         wa = PhotoInPlasma
        elseif   sa in [ "PhotoExc", "photoexcitation"]             wa = PhotoExc
        elseif   sa in [ "photo-EA"]                                wa = PhotoExcAuto
        elseif   sa in [ "photo-IA"]                                wa = PhotoIonAuto
        elseif   sa in [ "photo-IF"]                                wa = PhotoIonFluor
        elseif   sa in [ "photo-EF"]                                wa = PhotoExcFluor
        elseif   sa in [ "radiative", "Radiative"]                  wa = Radiative
        elseif   sa in [ "radiativeAuger", "RAuger"]                wa = RAuger
        elseif   sa in [ "rec", "Rec"]                              wa = Rec
        elseif   sa in [ "pair-annihilation-1-photon", "PA1P"]      wa = PairA1P
        elseif   sa in [ "impact-EA"]                               wa = ImpactExcAuto
        elseif   sa in [ "rec", "Rec"]                              wa = Rec
        else     error("stop a")
        end

        AtomicProcess(wa)
    end

    #===========================
    # `Base.show(io::IO, process::AtomicProcess)`  ... prepares a proper printout of the variable process::AtomicProcess.
    function Base.show(io::IO, process::AtomicProcess)
        print(io, string(process) )
    end   =======================#

    # `Base.string(process::AtomicProcess)`  ... provides a proper printout of the variable process::AtomicProcess.
    function Base.string(process::AtomicProcess)
        if      process == NoProcess       return( "no process" )
        elseif  process == Auger           return( "Auger" )
        elseif  process == AugerInPlasma   return( "Auger in plasma" )
        elseif  process == Compton         return( "Rayleigh-Compton" )
        elseif  process == Dierec          return( "Dielectronic recombination" )
        elseif  process == Eimex           return( "Eimex" )
        elseif  process == MultiPhotonDE   return( "multi-photon excitation & decay" )
        elseif  process == Photo           return( "Photo" )
        elseif  process == PhotoExcFluor   return( "Photo-Excitation-Fluoresence" )
        elseif  process == PhotoExcAuto    return( "Photo-Excitation-Autoionization" )
        elseif  process == PhotoInPlasma   return( "Photo in plasma" )
        elseif  process == PhotoIonFluor   return( "Photo-Ionization-Fluoresence" )
        elseif  process == PhotoIonAuto    return( "Photo-Ionization-Autoionization" )
        elseif  process == PhotoExc        return( "Photo-Excitation" )
        elseif  process == Radiative       return( "Radiative" )
        elseif  process == RAuger          return( "Radiative Auger" )
        elseif  process == Rec             return( "Rec" )
        else    error("stop a")
        end
    end

    """
    `struct  Basics.EmMultipole`  ... defines a struct for the em multipoles, E1, M1, E2, ...

        + L           ::Int64    ... multipolarity
        + electric    ::Bool     ... True if the multipole is an electric one, and false otherwise.
    """
    struct  EmMultipole
        L             ::Int64
        electric      ::Bool
    end

    const  E1 = EmMultipole(1, true);     const  M1 = EmMultipole(1, false)
    const  E2 = EmMultipole(2, true);     const  M2 = EmMultipole(2, false)
    const  E3 = EmMultipole(3, true);     const  M3 = EmMultipole(3, false)
    const  E4 = EmMultipole(4, true);     const  M4 = EmMultipole(4, false)

    """
    `Basics.EmMultipole(sa::String)`  ... constructor for a given String.
    """
    function EmMultipole(sa::String)
        if       sa == "E1"    L = 1;   electric = true
        elseif   sa == "M1"    L = 1;   electric = false
        elseif   sa == "E2"    L = 2;   electric = true
        elseif   sa == "M2"    L = 2;   electric = false
        elseif   sa == "E3"    L = 3;   electric = true
        elseif   sa == "M3"    L = 3;   electric = false
        elseif   sa == "E4"    L = 4;   electric = true
        elseif   sa == "M4"    L = 4;   electric = false
        else     error("stop a")
        end

        EmMultipole(L, electric)
    end

    # `Base.show(io::IO, mp::EmMultipole)`  ... prepares a proper printout of the variable mp::EmMultipole.
    function Base.show(io::IO, mp::EmMultipole)
        if  mp.electric   sa = "E"*string(mp.L)    else    sa = "M"*string(mp.L)   end
        print(io, string(mp) )
    end

    # `Base.string(mp::EmMultipole)`  ... provides a proper printout of the variable mp::EmMultipole.
    function Base.string(mp::EmMultipole)
        mp.L < 1    &&    error("Improper multipolarity L = $(mp.L)")
        if  mp.electric   sa = "E"*string(mp.L)    else    sa = "M"*string(mp.L)   end
        return( sa )
    end

    """
    `Basics.multipole_p(mp::EmMultipole)`  ... returns the kind integer p=0 (magnetic) or p=1 (electric) of the given  mp::EmMultipole.
    """
    function multipole_p(mp::EmMultipole)
        if  mp.electric   return(1)    else    return(0)   end
    end

    """
    `@enum   Basics.EmGauge`
        ... defines a enumeration for the (allowed) gauges to deal with the radiation field; this differs from the UseGauge which can
            be selected for explicit computations.

            + NoGauge      ... No gauge define (yet).
            + Coulomb      ... Coulomb gauge (= velocity gauge in the non-relativistic limit; Grant, 1974).
            + Babushkin    ... Babushkin gauge (= length gauge in the non-relativistic limit; Grant, 1974).
            + Magnetic     ... Magnetic gauge (transverse velocity gauge).
            + Velocity     ... Transverse velocity gauge (Johnson, 2007).
            + Length       ... Length gauge (Johnson, 2007).
    """
    @enum   EmGauge    NoGauge    Coulomb   Babushkin   Magnetic   Velocity   Length
    export  EmGauge,   NoGauge,   Coulomb,  Babushkin,  Magnetic,  Velocity,  Length

    """
    `Basics.EmGauge(sa::String)`  ... constructor for a given String.
    """
    function EmGauge(sa::String)
        if       sa in ["none", "NoGauge"]         wa = NoGauge
        elseif   sa == "Coulomb"                   wa = Coulomb
        elseif   sa == "Babushkin"                 wa = Babushkin
        elseif   sa in ["Magnetic", "magnetic"]    wa = Magnetic
        else     error("stop a")
        end

        return( wa )
    end

    #===========================
    # `Base.show(io::IO, gauge::EmGauge)`  ... prepares a proper printout of the variable gauge::EmGauge.
    function Base.show(io::IO, gauge::EmGauge)
        print(io, string(gauge) )
    end  ====================#

    # `Base.string(gauge::EmGauge)`  ... provides a proper printout of the variable gauge::EmGauge.
    function Base.string(gauge::EmGauge)
        if      gauge == NoGauge     return( "no gauge" )
        elseif  gauge == Coulomb     return( "Coulomb" )
        elseif  gauge == Babushkin   return( "Babushkin" )
        elseif  gauge == Magnetic    return( "Magnetic" )
        else    error("stop a")
        end
    end

    """
    `@enum   Basics.UseGauge`  ... defines a enumeration for the (allowed) gauges to be selected in explicit computations
    """
    @enum   UseGauge    UseCoulomb   UseBabushkin
    export  UseGauge,   UseCoulomb,  UseBabushkin

    """
    `Basics.UseGauge(sa::String)`  ... constructor for a given String.
    """
    function UseGauge(sa::String)
        if       sa == "Coulomb"                   wa = UseCoulomb
        elseif   sa == "Babushkin"                 wa = UseBabushkin
        else     error("stop a")
        end

        UseGauge(wa)
    end

    #===========================
    # `Base.show(io::IO, gauge::UseGauge)`  ... prepares a proper printout of the variable gauge::UseGauge.
    function Base.show(io::IO, gauge::UseGauge)
        print(io, string(gauge) )
    end    ==========================#

    # `Base.string(gauge::UseGauge)`  ... provides a proper printout of the variable gauge::EmGauge.
    function Base.string(gauge::UseGauge)
        if      gauge == UseCoulomb     return( "Coulomb" )
        elseif  gauge == UseBabushkin   return( "Babushkin" )
        else    error("stop a")
        end
    end

    """
    `struct  EmProperty`  ... defines a type to maintain two gauge forms of a computed result that depends on the radiation field.

        + Coulomb         ::Float64    ... Value for the Coulomb gauge of the radiation field.
        + Babushkin       ::Float64    ... Value for the Coulomb gauge of the radiation field.
    """
    struct  EmProperty
        Coulomb           ::Float64
        Babushkin         ::Float64
    end

    """
    `Basics.EmProperty(wa::Float64)`  ... constructor for an `constant` instance of EmProperty.
    """
    function EmProperty(wa::Float64)
        EmProperty(wa, wa)
    end

    Base.:+(a::EmProperty, b::EmProperty) = EmProperty(a.Coulomb + b.Coulomb, a.Babushkin + b.Babushkin)
    Base.:+(a::EmProperty, b) = EmProperty(a.Coulomb + b, a.Babushkin + b)
    Base.:+(a, b::EmProperty) = b + a
    Base.:/(a::EmProperty, b::EmProperty) = EmProperty(a.Coulomb / b.Coulomb, a.Babushkin / b.Babushkin)
    Base.:/(a, b::EmProperty) = EmProperty(a / b.Coulomb, a / b.Babushkin)

    # `Base.show(io::IO, property::EmProperty)`  ... prepares a proper printout of the variable property::EmProperty.
    function Base.show(io::IO, property::EmProperty)
        sa = Base.string(property);                print(io, sa)
    end

    # `Base.string(property::EmProperty)`  ... provides a String notation for the variable property::EmProperty.
    function Base.string(property::EmProperty)
        sa = "$(property.Coulomb) [Coulomb],  $(property.Babushkin) [Babushkin]"
        return( sa )
    end

    """
    `struct  Basics.ExpStokes`
        ... defines a type to maintain the (experimentally) given Stokes parameter for the polarization of incoming
            photons or electrons.

        + P1              ::Float64    ... Stokes P1 parameter.
        + P2              ::Float64    ... Stokes P2 parameter.
        + P3              ::Float64    ... Stokes P3 parameter.
    """
    struct  ExpStokes
        P1                ::Float64
        P2                ::Float64
        P3                ::Float64
    end

    """
    `Basics.ExpStokes()`  ... constructor for an (empty or unpolarized) Stokes vector.
    """
    function ExpStokes()
        ExpStokes(0., 0., 0.)
    end

    # `Base.show(io::IO, ps::ExpStokes)`  ... prepares a proper printout of the variable ps::ExpStokes.
    function Base.show(io::IO, ps::ExpStokes)
        sa = Base.string(ps);                print(io, sa)
    end

    # `Base.string(ps::ExpStokes)`  ... provides a String notation for the variable ps::ExpStokes.
    function Base.string(ps::ExpStokes)
        sa = "exp. Stokes parameters P1 = $(ps.P1),  P2 = $(ps.P2),  P3 = $(ps.P3)"
        return( sa )
    end

    """
    `struct  Basics.EmStokes`
        ... defines a type to maintain the (computed) given Stokes parameter for the polarization of emitted photons or electrons.

        + P1              ::EmProperty    ... Stokes P1 parameter in Coulomb and Babushkin gauge
        + P2              ::EmProperty    ... Stokes P2 parameter.
        + P3              ::EmProperty    ... Stokes P3 parameter.
    """
    struct  EmStokes
        P1                ::EmProperty
        P2                ::EmProperty
        P3                ::EmProperty
    end

    # `Base.show(io::IO, ps::EmStokes)`  ... prepares a proper printout of the variable ps::EmStokes.
    function Base.show(io::IO, ps::EmStokes)
        sa = Base.string(ps);                print(io, sa)
    end

    # `Base.string(ps::EmStokes)`  ... provides a String notation for the variable ps::EmStokes.
    function Base.string(ps::EmStokes)
        sa = "Stokes parameter P1 = $(ps.P1),  P2 = $(ps.P2),  P3 = $(ps.P3)  "
        return( sa )
    end

    """
    `struct  Basics.TensorComp`  ... defines a type for a component of the statistical tensor as associated with an atomic or ionic level.

        + k               ::Int64          ... rank of the tensor component
        + q               ::Int64          ... q-component
        + comp            ::ComplexF64     ... value of the tensor component
    """
    struct  TensorComp
        k                 ::Int64
        q                 ::Int64
        comp              ::ComplexF64
    end

    # `Base.show(io::IO, st::TensorComp)`  ... prepares a proper printout of the variable st::TensorComp.
    function Base.show(io::IO, st::TensorComp)
        sa = Base.string(st);                print(io, sa)
    end

    # `Base.string(st::TensorComp)`  ... provides a String notation for the variable st::TensorComp.
    function Base.string(st::TensorComp)
        sa = "Statistical tensor component rho_($(st.k),$(st.q)) = $(st.comp)  "
        return( sa )
    end

    """
    `Basics.tensorComp(k::Int64, q::Int64, tcomps::Array{TensorComp,1};  withZeros::Bool = false)`
        ... returns the value of the statistical tensor component rho_kq if it is contained in tcomps.
    """
    function tensorComp(k::Int64, q::Int64, tcomps::Array{TensorComp,1};  withZeros::Bool = false)
        for  tc in  tcomps
            if   k == tc.k    &&    q == tc.q    return( tc.comp )    end
        end

        if   withZeros    return( ComplexF64(0.) )    else    error("Statistical tensor component not found for k = $k and q = $q ")    end
    end

    """
    `@enum   Basics.AtomicLevelProperty`  ... defines a enumeration for the atomic level properties that are supported in the JAC program.

        + NoProperty      ... No level property defined.
        + Einstein        ... Einstein A, B coefficients and oscillator strength; although not a 'level property', `Einstein` is treated within
                            a single basis and does not support to include relaxation effects, etc. It can be used for a quick overwiew to
                            transition probabilities or in cascade computations, however,
        + HFS             ... Hyperfine A and B parameters.
        + Isotope         ... Isotope shift M and F parameters.
        + LandeJ          ... Lande g_J factors.
        + LandeF          ... Lande g_F factors.
        + Polarizibility  ... static and dynamic polarizibilities.
        + Plasma          ... CI computations including interactions from various plasma models.
        + Zeeman          ... Zeeman splitting of fine-structure levels.
    """
    @enum   AtomicLevelProperty   NoProperty   AlphaX   EinsteinX   FormF   Green  HFS   Isotope   LandeJ   LandeF  Plasma=
                            20   Polarity   Yields  Zeeman
    export  AtomicLevelProperty,  NoProperty,  AlphaX,  EinsteinX,  FormF,  Green, HFS,  Isotope,  LandeJ,  LandeF, Plasma,
                                Polarity,  Yields,  Zeeman

    """
    `Basics.AtomicLevelProperty(sa::Union{String,SubString{String}})`  ... constructor for a given String.
    """
    function AtomicLevelProperty(sa::Union{String,SubString{String}})
        sa = strip( Base.convert(String, sa) )
        if       sa in ["none", "NoProperty"]                wa = NoProperty
        elseif   sa in ["alpha"]                             wa = AlphaX
        elseif   sa in ["Einstein"]                          wa = EinsteinX
        elseif   sa in ["form factor"]                       wa = FormF
        elseif   sa in ["Einstein"]                          wa = EinsteinX
        elseif   sa in ["Green", "Greens function"]          wa = Green
        ##x elseif   sa in ["IJF_expansion", "IJF"]              wa = IJF_expansion
        elseif   sa == ["isotope shift", "IS"]               wa = Isotope
        elseif   sa == ["LandeJ", "g_J"]                     wa = LandeJ
        elseif   sa == ["LandeF", "g_F"]                     wa = LandeF
        elseif   sa == ["Plasma", "plasma"]                  wa = Plasma
        elseif   sa == ["polarity", "polarizibility"]        wa = Polarity
        elseif   sa == ["yields", "fluorescence yield"]      wa = Yields
        elseif   sa == ["Zeeman"]                            wa = Zeeman
        else     error("stop a")
        end

        AtomicLevelProperty(wa)
    end

    #===========================
    # `Base.show(io::IO, property::AtomicLevelProperty)`  ... prepares a proper printout of the variable property::AtomicLevelProperty.
    function Base.show(io::IO, property::AtomicLevelProperty)
        print(io, string(property) )
    end   ======================#

    # `Base.string(property::AtomicLevelProperty)`  ... provides a proper printout of the variable property::AtomicLevelProperty.
    function Base.string(property::AtomicLevelProperty)
        if      property == NoProperty    return( "no property" )
        elseif  property == AlphaX        return( "alpha variation" )
        elseif  property == EinsteinX     return( "Einstein A, B coefficients" )
        elseif  property == FormF         return( "form & scattering factors" )
        elseif  property == Green         return( "approximate Greens function" )
        elseif  property == HFS           return( "HFS A, B parameters" )
        elseif  property == Isotope       return( "isotope shift M and F parameters" )
        elseif  property == LandeJ        return( "Lande g_J factors" )
        elseif  property == LandeF        return( "Lande g_F factors" )
        elseif  property == Plasma        return( "CI computations with plasma-type interactions" )
        elseif  property == Polarity      return( "multipole polarizibility" )
        elseif  property == Yields        return( "fluorescence & Auger yields" )
        elseif  property == Zeeman        return( "Zeeman splitting in an external magnetic field" )
        else    error("stop a")
        end
    end

    """
    `@enum   Basics.RadialMesh`  ... defines a enumeration for dealing with the radial mesh points and integration scheme.

        + MeshGrasp       ... use the exponential or exponential-linearized grid from Grasp92.
        + MeshGL          ... use a Gauss-Legendre grid with grid points due to the orderGL.
    """
    @enum   RadialMesh    MeshGrasp    MeshGL
    export  RadialMesh,   MeshGrasp,   MeshGL

    #==
    """
    `@enum   Basics.QedModel`  ... defines a enumeration for dealing with the radiative (QED) contribution to the Hamiltonian matrix.

        + QedSydney       ... to use the local QED potentials by Flambaum & Ginges PRA 72, 052115 (2005).
        + QedPetersburg   ... to use the local QED potential similar as suggested by Shabaev et al., PRA 88, 012513 (2013).
    """
    @enum   QedModel    QedSydney    QedPetersburg
    export  QedModel,   QedSydney,   QedPetersburg   ==#

    """
    `@enum   Basics.ContinuumNormalization`  ... defines a enumeration for dealing with the normalization of continuum orbitals.

        + PureSine       ... normalize with regard to an (asymtotic) pure sine funtion, sin(kr).
        + CoulombSine    ... normalize with regard to an (asymtotic) Coulombic-sine funtion, sin(kr + ...).
        + OngRussek      ... normalize by following Ong & Russek (1973).
    """
    @enum   ContinuumNormalization    PureSine    CoulombSine    OngRussek
    export  ContinuumNormalization,   PureSine,   CoulombSine,   OngRussek

    """
    `@enum   Basics.ContinuumSolutions`  ... defines a enumeration for solving the continuum orbitals in a given potential.

        + ContBessel              ... generate a pure Bessel function for the large component together with kinetic balance.
        + ContSine                ... generate a pure Sine function for the large component together with kinetic balance.
        + NonrelativisticCoulomb  ... generate a non-relativistic Coulomb function for the large component together with kinetic balance.
        + AsymptoticCoulomb       ... generate a pure (asymptotic) Coulombic function for both components.
        + BsplineGalerkin         ... generate a continuum orbital with the Galerkin method.
    """
    @enum   ContinuumSolutions    ContBessel    ContSine     AsymptoticCoulomb    NonrelativisticCoulomb    BsplineGalerkin
    export  ContinuumSolutions,   ContBessel,   ContSine,    AsymptoticCoulomb,   NonrelativisticCoulomb,   BsplineGalerkin

    """
    `@enum   Basics.Warnings`      ... defines a enumeration for dealing with warnings that are made during a run or REPL session.
                                Cf. Defaults.warn().

        + AddWarning        ... add a Warning to a warningList.
        + PrintWarnings     ... print all warnings into a jac-warn.report file.
        + ResetWarnings     ... reset (empty) the warningList, usually at the beginning of a new run.
    """
    @enum   Warnings    AddWarning    PrintWarnings    ResetWarnings
    export  Warnings,   AddWarning,   PrintWarnings,   ResetWarnings

    """
    `@enum   Basics.Guint`      ... defines a enumeration for dealing with graphical user interfaces (GUI).

        + Gui            ... use a graphical (interactive) user interface.
        + GuiSettings    ... use a graphical (interactive) user interface to define some settings.
        + Cascade        ... use a graphical (interactive) user interface to define some settings.
    """
    @enum   Guint    Gui    GuiSettings   GuiCascade
    export  Guint,   Gui,   GuiSettings,  GuiCascade

    """
    `abstract type Basics.AbstractExcitationScheme`
        ... defines an abstract and a number of singleton types to distinguish between different schemes for
            generating configuration lists as they frequently occur in Green function and cascade computations.

      + struct NoExcitationScheme
        ... dummy scheme for (unsupported) initialization of this abstract tpye.

      + struct DeExciteSingleElectron
        ... generates all excitations and de-excitations of a single electron from a given list of bound
            electron configurations. The number of electrons of the generated configurations is the same as
            for the given bound configurations.

      + struct DeExciteTwoElectrons
        ... generates all excitations and de-excitations of one or two electrons from a given list of bound
            electron configurations. The number of electrons of the generated configurations is the same as
            for the given bound configurations.

      + struct AddSingleElectron
        ... generates configurations by just adding a single electrons to a given list of bound
            electron configurations. The number of electrons of the generated configurations is N+1.
    """
    abstract type  AbstractExcitationScheme                               end
    struct         NoExcitationScheme      <:  AbstractExcitationScheme   end
    struct         DeExciteSingleElectron  <:  AbstractExcitationScheme   end
    struct         DeExciteTwoElectrons    <:  AbstractExcitationScheme   end
    struct         AddSingleElectron       <:  AbstractExcitationScheme   end

    # Functions/methods that are later added to the module Basics
    function add                                            end
    function addZerosToCsfR                                 end
    function analyze                                        end
    function analyzeConvergence                             end
    function compute                                        end
    function computeDiracEnergy                             end
    function computeMeanSubshellOccupation                  end
    function computePotentialCoreHartree                    end
    function computePotentialHartree                        end
    function computePotentialHartreeSlater                  end
    function computePotentialKohnSham                       end
    function computePotentialDFS                            end
    function computePotentialExtendedHartree                end
    function determineEnergySharings                        end
    function determineHoleShells                            end
    function determineParity                                end
    function determineSelectedLines                         end
    function determineSelectedPathways                      end
    function diagonalize                                    end
    function display                                        end
    function excludeDoubles                                 end
    function extractNoOpenShells                            end
    function extractNonrelativisticShellList                end
    function extractNonrelativisticConfigurations           end
    function extractNonrelativisticConfigurationFromCsfR    end
    function extractOpenShellQNfromCsfNR                    end
    function extractOpenShellQNfromCsfR                     end
    function extractShellOccupationFromCsfR                 end
    function generate()                                     end
    function generateConfigurationsForExcitationScheme      end
    function generateLevelWithExtraElectron                 end
    function generateLevelWithExtraSubshell                 end
    function generateLevelWithSymmetryReducedBasis          end
    function generateShellList                              end
    function integrate                                      end
    function integrateOnGridNewtonCotes                     end
    function integrateOnGridSimpsonRule                     end
    function integrateOnGridTrapezRule                      end
    function interpolateOnGridGrasp92                       end
    function interpolateOnGridTrapezRule                    end
    function isSimilar                                      end
    function isStandardSubshellList                         end
    function merge                                          end
    function modifyLevelEnergies                            end
    function perform                                        end
    function plot                                           end
    function read                                           end
    function readCslFileGrasp92                             end
    function readOrbitalFileGrasp92                         end
    function readMixFileRelci                               end
    function readMixingFileGrasp92                          end
    function recast                                         end
    function sortByEnergy                                   end
    function tabulate                                       end
    function tabulateKappaSymmetryEnergiesDirac             end
    function tools                                          end
    function yesno                                          end

end # module


text:
using Distributions
using PyPlot
using Random
using Statistics

function plot_p_sigmoid(; lb, ub, alpha, filename=nothing)
    chisq_dist = Chisq(3.0)
    xs = [x for x in lb:0.05:ub]

    figure(figsize=(9.0, 6.0))
    pdfs = pdf.(chisq_dist, xs)
    plt.plot(xs, 5.0 * pdfs, label="cost", lw=1.0)

    alpha_quantile = quantile(chisq_dist, 1.0 - alpha)
    ys = [y for y in 0.0:0.05:1.2]
    plt.plot(fill(alpha_quantile, length(ys)), ys, label="α-quantile", linestyle="--")

    sigmoid_lw = 0.7
    sigmoid_1_0 = (x) -> 1 / (1 + exp.(-(x - alpha_quantile)))
    plt.plot(xs, sigmoid_1_0.(xs), label="sigmoid1", lw=sigmoid_lw)
    sigmoid_0_5 = (x) -> 1 / (1 + exp.(-0.5 * (x - alpha_quantile)))
    plt.plot(xs, sigmoid_0_5.(xs), label="sigmoid2", lw=sigmoid_lw)
    sigmoid_1_5 = (x) -> 1 / (1 + exp.(-2.0 * (x - alpha_quantile)))
    plt.plot(xs, sigmoid_1_5.(xs), label="sigmoid3", lw=sigmoid_lw)

    plt.legend()

    if !isnothing(filename)
        plt.savefig(filename, dpi=500)
    end
end

plot_p_sigmoid(lb=0.0, ub=15.0, alpha=0.05, filename="/Users/kykim/Desktop/sigmoid.png")

ext: jl


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

const ERR_BUF = 10

_terminate!(A::_ACT, reason) = isnothing(A.term) || Base.invokelatest(A.term.f, reason)
_sticky(A::_ACT) = A.mode in (:system, :sticky, :supervisor)
function trysend(lk::Link, msg...)
    try
        send(lk, msg...)
    catch
        nothing
    end
end

function saveerror(t::Task)
    err = try
        task_local_storage("_ERR")
    catch exc
        exc isa KeyError ?
            task_local_storage("_ERR", Task[]) :
            rethrow()
    end
    length(err) ≥ ERR_BUF && popfirst!(err)
    push!(err, t)
end
saveerror(::Nothing) = nothing
function errored()
    try
        task_local_storage("_ERR")
    catch exc
        exc isa KeyError ? nothing : rethrow()
    end
end

function onerror(A::_ACT, exc)
    for c in A.conn
        c isa Monitor ?
            trysend(c.lk, Down(self(), exc, current_task())) :
            c isa Super ?
                trysend(c.lk, Exit(exc, self(), current_task(), A)) :
                trysend(c.lk, Exit(exc, self(), current_task(), nothing))
    end
end

#
# Actors protocol on Down, Exit
#

# Down
function onmessage(A::_ACT, msg::Down)
    ix = 0
    for (i, c) in enumerate(A.conn)
        if c.lk == msg.from && c isa Monitored
            !isnothing(c.action) ?
                c.action(msg.reason) :
                log_warn(msg, "monitored")
            ix = i
        end
    end
    ix != 0 ?
        deleteat!(A.conn, ix) :
        log_warn(msg)
end

# Exit
function onmessage(A::_ACT, msg::Exit)
    for c in A.conn
        if c.lk != msg.from
            c isa Monitor ?
                trysend(c.lk, Down(self(), msg.reason, current_task())) :
                send(c.lk, Exit(msg.reason, self(), msg.task, A))
        end
    end
    _terminate!(A, msg.reason)
end
onmessage(A::_ACT, ::Val{:sticky}, msg::Exit) = onmessage(A, Val(:system), msg)
function onmessage(A::_ACT, ::Val{:system}, msg::Exit)
    if msg.reason != :normal
        saveerror(msg.task)
        log_warn(msg, "connected")
    end
    ix = findfirst(c->c.lk==msg.from, A.conn)
    isnothing(ix) ? # Exit not from a connection
        A.mode = Symbol(string(A.mode)*"∇") :
        deleteat!(A.conn, ix)
end
function onmessage(A::_ACT, ::Val{:supervisor}, msg::Exit)
    !isnormal(msg.reason) && saveerror(msg.task)
    ix = findfirst(c->c.lk==msg.from, A.conn)
    if !isnothing(ix) && A.conn[ix] isa Child
        log_warn(msg, "supervised")
        A.bhv(msg)
    elseif !isnothing(ix) && A.conn[ix] isa Peer
        msg.reason != :normal && log_warn(msg, "connected")
    else
        if msg.reason ∉ (:done, :timed_out)
            for c in A.conn
                if c.lk != msg.from
                    c isa Monitor ? trysend(c.lk, Down(self(), msg.reason)) :
                    c isa Child   ? shutdown_child(c) :
                        send(c.lk, Exit(msg.reason, self(), msg.task, A))
                end
            end
            log_warn(msg, "supervisor")
            _terminate!(A, msg.reason)
            A.mode = Symbol(string(A.mode)*"∇")
        end
    end
end


text:
using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials, ForwardDiff, Plots, StaticArrays
import ForwardDiff: derivative, jacobian, Dual
import SemiclassicalOrthogonalPolynomials: Weighted
import ClassicalOrthogonalPolynomials: associated, affine

Base.floatmin(::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T,V,N}(floatmin(V))
Base.big(d::Dual{T,V,N}) where {T,V,N} = Dual{T}(big(d.value), ForwardDiff.Partials(map(big,d.partials.values)))

####
#
# We first do a single interval.
# Equilibrium measures for a symmetric potential
#  with one interval of support
# a measure w(x) supported on
# [-b,b]
# such that
#       1. H*w == V'
#       2. sum(w) == 1
#       3.  w is bounded
#
#  rescaling x == b*t these becomes find
# a measure w̃(t) supported on
# [-1,1]
# such that
#       1. H*w̃ == V'(b*t)
#       2. sum(w̃) == 1/b
#       3.  w is bounded
#
# Note (1) and (2) can always be satisfied
# thus the constraint comes from (3).
# The following gives the evaluation of the
# unweighted-component of the measure evaluated
# at (a,b)
#####

V = x -> x^2
function equilibriumcoefficients(T, b)
    U = associated(T)
    W = Weighted(T)
    t = axes(W,1)
    H = @. inv(t - t')
    H̃ = U \H*W
    [1/(b*sum(W[:,1])); 2H̃[:,2:end] \ ( U \ derivative.(V, b*t))]
end
function equilibriumendpointvalue(b::Number)
    T = ChebyshevT{typeof(b)}()
    dot(T[end,:], equilibriumcoefficients(T,b))
end

function equilibrium(b::Number)
    T = ChebyshevT{typeof(b)}()
    U = ChebyshevU{typeof(b)}()
    # convert to Weighted(U) to make value at ±b accurate
    Weighted(U)[affine(-b..b,axes(T,1)),:] * ((Weighted(T) \ Weighted(U))[3:end,:] \ equilibriumcoefficients(T,b)[3:end])
end

μ = equilibrium(sqrt(2))

T = Chebyshev()
b = sqrt(2)
μ = Weighted(T) * equilibriumcoefficients(T, b)
x = axes(μ,1)

plot(μ)

xx = 0.7; 2b*(log.(abs.(x .- x'))*μ)[xx] - V(b*xx)

b = 1.0 # initial guess
for _ = 1:10
    b -= derivative(equilibriumendpointvalue,b) \ equilibriumendpointvalue(b)
end
b

plot(equilibrium(b))

#####
# Equilibrium measures for a symmetric potential
# with two intervals of support consists of finding
# a measure w(x) supported on
# [-b,-a] ∪ [a,b]
# such that
#       1. H*w == V'
#       2. sum(w) == 1
#       3.  w is bounded
#
#  rescaling x == b*t these becomes find
# a measure w̃(t) supported on
# [-1,-a/b] ∪ [a/b,1]
# such that
#       1. H*w̃ == V'(b*t)
#       2. sum(w̃) == 1/b
#       3.  w is bounded
#
# Note (1) and (2) can always be satisfied
# thus the two constraints come from (3).
# The following gives the evaluation of the
# unweighted-component of the measure evaluated
# at (a,b)
#####
V = x -> x^4 - 10x^2
function equilibriumcoefficients(P,a,b)
    W = Weighted(P)
    Q = associated(P)
    t = axes(W,1)
    x = axes(Q,1)
    H = @. inv(x - t')
    H̃ = Q \ H*W
    [1/(b*sum(W[:,1])); 2H̃[:,2:end] \( Q \ derivative.(V, b*x))]
end
function equilibriumendpointvalues(ab::SVector{2})
    a,b = ab
    # orthogonal polynomials w.r.t.
    # abs(x) / (sqrt(1-x^2) * sqrt(x^2 - ρ^2))
    P = TwoBandJacobi(a/b, -one(a)/2, -one(a)/2, one(a)/2)
    Vector(P[[a/b,1],:] * equilibriumcoefficients(P,a,b))
end

function equilibrium(ab)
    a,b = ab
    P = TwoBandJacobi(a/b, -1/2, -1/2, 1/2)
    Weighted(P) * equilibriumcoefficients(P,a,b)
end

ab = SVector(2.,3.)
ab -= jacobian(equilibriumendpointvalues,ab) \ equilibriumendpointvalues(ab)
a,b = ab
xx = range(-4,4;length=1000)
μ = equilibrium(ab)
μx = x -> a < abs(x) < b ? μ[x/b] : 0.0
plot!(xx, μx.(xx))

plot(equilibrium(ab))


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "MFEM_seq"
version = v"4.4.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/mfem/mfem.git",
              "a1f6902ed72552f3e680d1489f1aa6ade2e0d3b2"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/mfem*

# Backport https://github.com/mfem/mfem/pull/2995 to be able to compile with Musl
atomic_patch -p1 ../patches/0001-use_overloaded_function_definitions_instead_of_unreliable_macros.diff

mkdir build; cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_SHARED_LIBS=1 \
    ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

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

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"8")


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# ToDo: Add `density_logvalgrad!` to support HMC, etc.

doc"""
    AbstractDensity

The following functions must be implemented for subtypes:

* `BAT.nparams`
* `BAT.unsafe_density_logval`

In some cases, it may be desirable to override the default implementations
of the functions

* `BAT.exec_capabilities`
* `BAT.unsafe_density_logval!`
"""
abstract type AbstractDensity end
export AbstractDensity

Base.rand(rng::AbstractRNG, density::AbstractDensity, T::Type{<:AbstractFloat}) =
    rand!(rng, density, Vector{T}(nparams(density)))

Base.rand(rng::AbstractRNG, density::AbstractDensity, T::Type{<:AbstractFloat}, n::Integer) =
    rand!(rng, density, Matrix{T}(nparams(density), n))

Base.rand!(rng::AbstractRNG, density::AbstractDensity, x::StridedVecOrMat{<:Real}) =
    rand!(rng, sampler(density), x)

doc"""
    param_bounds(density::AbstractDensity)::AbstractParamBounds

Get the parameter bounds of `density`. See `density_logval!` for the
implications and handling of the bounds.

Use

   new_density = density[bounds::ParamVolumeBounds]

to create a new density function with additional bounds.
"""
function param_bounds(density::AbstractDensity)
    NoParamBounds(nparams(density))
end
export param_bounds

# doc"""
#     getindex(density::AbstractDensity, bounds::ParamVolumeBounds)
#
# Limit `density` to `bounds`. See `param_bounds` and `density_logval!`.
# """
# Base.getindex(density::AbstractDensity, bounds::ParamVolumeBounds) =
#     density * ConstDensity(bounds, false)

Base.convert(::Type{AbstractDensity}, bounds::ParamVolumeBounds) =
    ConstDensity(bounds, one)

doc"""
    density_logval(
        density::AbstractDensity,
        params::AbstractVector{<:Real},
        exec_context::ExecContext = ExecContext()
    )

Version of `density_logval` for a single parameter vector.

Do not implement `density_logval` directly for subtypes of
`AbstractDensity`, implement `BAT.unsafe_density_logval` instead.

See `ExecContext` for thread-safety requirements.
"""
function density_logval(
    density::AbstractDensity,
    params::AbstractVector{<:Real},
    exec_context::ExecContext = ExecContext()
)
    !(size(params, 1) == nparams(density)) && throw(ArgumentError("Invalid length of parameter vector"))
    r = unsafe_density_logval(density, params, exec_context)
    isnan(r) && throw(ErrorException("Return value of unsafe_density_logval must not be NaN"))
    r
end
export density_logval

# Assume that density_logval isn't always thread-safe, but usually remote-safe:
exec_capabilities(::typeof(density_logval), density::AbstractDensity, params::AbstractVector{<:Real}) =
    exec_capabilities(unsafe_density_logval, density, params)

doc"""
    BAT.unsafe_density_logval(
        density::AbstractDensity,
        params::AbstractVector{<:Real},
        exec_context::ExecContext = ExecContext()
    )

Unsafe variant of `density_logval`, implementations may rely on

* `size(params, 1) == nparams(density)`

The caller *must* ensure that these conditions are met!
"""
function unsafe_density_logval end

# Assume that density_logval isn't always thread-safe, but usually remote-safe:
exec_capabilities(::typeof(unsafe_density_logval), density::AbstractDensity, args...) =
    ExecCapabilities(1, false, 1, true)

doc"""
    density_logval!(
        r::AbstractArray{<:Real},
        density::AbstractDensity,
        params::AbstractMatrix{<:Real},
        exec_context::ExecContext = ExecContext()
    )

Compute log of values of a density function for multiple parameter value
vectors.

Input:

* `density`: density function
* `params`: parameter values (column vectors)
* `exec_context`: Execution context

Output is stored in

* `r`: Array of log-result values, length must match, shape is ignored

Array size requirements:

    size(params, 2) == length(r)

Note: `density_logval!` must not be called with out-of-bounds parameter
vectors (see `param_bounds`). The result of `density_logval!` for parameter
vectors that are out of bounds is implicitly `-Inf`, but for performance
reasons the output is left undefined: `density_logval!` may fail or store
arbitrary values in `r`.

Do not implement `density_logval!` directly for subtypes of
`AbstractDensity`, implement `BAT.unsafe_density_logval!` instead.

See `ExecContext` for thread-safety requirements.
"""
function density_logval!(
    r::AbstractVector{<:Real},
    density::AbstractDensity,
    params::AbstractMatrix{<:Real},
    exec_context::ExecContext = ExecContext()
)

    !(size(params, 1) == nparams(density)) && throw(ArgumentError("Invalid length of parameter vector"))
    !(indices(params, 2) == indices(r, 1)) && throw(ArgumentError("Number of parameter vectors doesn't match length of result vector"))
    unsafe_density_logval!(r, density, params, exec_context)
    any(isnan, r) && throw(ErrorException("unsafe_density_logval! must not set any return value to NaN"))
    r
end
export density_logval!

exec_capabilities(::typeof(density_logval!), r::AbstractArray{<:Real}, density::AbstractDensity, params::AbstractMatrix{<:Real}) =
    exec_capabilities(unsafe_density_logval!, r, density, params)

doc"""
    BAT.unsafe_density_logval!(
        r::AbstractArray{<:Real},
        density::AbstractDensity,
        params::AbstractMatrix{<:Real},
        exec_context::ExecContext
    )

Unsafe variant of `density_logval!`, implementations may rely on

* `size(params, 1) == nparams(density)`
* `size(params, 2) == length(r)`

The caller *must* ensure that these conditions are met!
"""
function unsafe_density_logval!(
    r::AbstractArray{<:Real},
    density::AbstractDensity,
    params::AbstractMatrix{<:Real},
    exec_context::ExecContext
)
    # TODO: Support for parallel execution
    single_ec = exec_context # Simplistic, will have to change for parallel execution
    for i in eachindex(r, indices(params, 2))
        p = view(params, :, i) # TODO: Avoid memory allocation
        r[i] = density_logval(density, p, single_ec)
    end
    r
end

# ToDo: Derive from exec_capabilities(density_logval, density, ...)
exec_capabilities(::typeof(unsafe_density_logval!), r::AbstractArray{<:Real}, density::AbstractDensity, args...) =
    ExecCapabilities(1, false, 1, true) # Change when default implementation of density_logval! for AbstractDensity becomes multithreaded.

doc"""
    GenericDensity{F} <: AbstractDensity

Constructors:

    GenericDensity(log_f, nparams::Int)

Turns the logarithmic density function `log_f` into a
BAT-compatible `AbstractDensity`. `log_f` must support

    `log_f(params::AbstractVector{<:Real})::Real`

with `length(params) == nparams`.

It must be safe to execute `log_f` in parallel on multiple threads and
processes.
"""
struct GenericDensity{F} <: AbstractDensity
    log_f::F
    nparams::Int
end

export GenericDensity

Base.parent(density::GenericDensity) = density.log_f

nparams(density::GenericDensity) = density.nparams

function unsafe_density_logval(
    density::GenericDensity,
    params::AbstractVector{<:Real},
    exec_context::ExecContext = ExecContext()
)
    size(params, 1) != nparams(density) && throw(ArgumentError("Invalid number of parameters"))
    density.log_f(params)
end

exec_capabilities(::typeof(unsafe_density_logval), density::GenericDensity, params::AbstractVector{<:Real}) =
    ExecCapabilities(1, true, 1, true)

#=

mutable struct TransformedDensity{...}{
    SO<:AbstractDensity,
    SN<:AbstractDensity
} <: AbstractDensity
   before::SO
   after::SN
   # ... transformation, Jacobi matrix of transformation, etc.
end

export TransformedDensity

...

=#


text:
import TriangleMesh
using Printf
using LoopVectorization

Base.@kwdef mutable struct FCFV_Mesh
    type   ::Union{String, Missing}          = missing
    nel    ::Union{Int64,  Missing}          = missing
    nf     ::Union{Int64,  Missing}          = missing
    nv     ::Union{Int64,  Missing}          = missing
    nn_el  ::Union{Int64,  Missing}          = missing
    nf_el  ::Union{Int64,  Missing}          = missing
    xv     ::Union{Vector{Float64}, Missing} = missing # node x coordinate
    yv     ::Union{Vector{Float64}, Missing} = missing # node y coordinate
    xf     ::Union{Vector{Float64}, Missing} = missing # face x coordinate
    yf     ::Union{Vector{Float64}, Missing} = missing # face y coordinate
    bc     ::Union{Vector{Int64},   Missing} = missing # face y coordinate
    xc     ::Union{Vector{Float64}, Missing} = missing # cent x coordinate
    yc     ::Union{Vector{Float64}, Missing} = missing # cent y coordinate
    e2v    ::Union{Matrix{Int64},   Missing} = missing # cell 2 node numbering
    e2f    ::Union{Matrix{Int64},   Missing} = missing # cell 2 face numbering
    vole   ::Union{Vector{Float64}, Missing} = missing # volume of element
    n_x    ::Union{Matrix{Float64}, Missing} = missing # normal 2 face x
    n_y    ::Union{Matrix{Float64}, Missing} = missing # normal 2 face y
    dA     ::Union{Matrix{Float64}, Missing} = missing # face length
end

function MakeTriangleMesh( nx, ny, xmin, xmax, ymin, ymax )

dx   = (xmax-xmin)/nx
dy   = (ymax-ymin)/ny
area = dx*dy

# Four corners of the domain
px   = [xmin xmax xmax xmin]
py   = [ymin ymin ymax ymax]
sx   = [ 1 2 3 4 ]
sy   = [ 2 3 4 1 ]
st   = [ 1 1 1 1 ]          # segment markers

p    = vcat(px, py)         # points
s    = vcat(sx, sy)         # segments

# inum = 0;
# dx
# px = Float64[]
# py = Float64[]
# sx = Int64[]
# sy = Int64[]
# st = Int64[]
# for i=1:nx+1
#     inum+=1
#     px   = push!(px, (i-1)*dx)
#     py   = push!(py,    ymin)
#     sx   = push!(sx, inum   )
#     sy   = push!(sy, inum+1 )
#     st   = push!(st, 1      )
# end

# for i=1:nx+1
#     inum+=1
#     px   = push!(px, (i-1)*dx)
#     py   = push!(py,    ymax)
#     sx   = push!(sx, inum   )
#     sy   = push!(sy, inum+1 )
#     st   = push!(st, 1      )
# end

# for i=1:ny+1
#     inum+=1
#     py   = push!(py, (i-1)*dy)
#     px   = push!(px,    xmin)
#     sx   = push!(sx, inum   )
#     sy   = push!(sy, inum+1 )
#     st   = push!(st, 1      )
# end

# for i=1:ny+1
#     inum+=1
#     py   = push!(py, (i-1)*dy)
#     px   = push!(px,    xmax)
#     sx   = push!(sx, inum   )
#     sy   = push!(sy, inum+1 )
#     st   = push!(st, 1      )
# end

# sy[end] = 1;

# p    = hcat(px, py)         # points
# s    = hcat(sx, sy)         # segments

# p = p'
# s = s'

st = st[:]

# Triangulation
holes    = Array{Float64}(undef,2,0)
domain   = TriangleMesh.Polygon_pslg(size(p,2), p, 0, Array{Int64}(undef,2,0), 0, Array{Float64}(undef,2,0),  size(s,2), s, st, 0, holes)
astring = @sprintf("%0.10lf", area)
switches = "Dpenq33o2IAa$(astring)"

println(switches)
trimesh  = TriangleMesh.create_mesh(domain, switches)
nvert_el = 3 # vertices per element

mesh        = FCFV_Mesh()
mesh.type   = "triangle"
mesh.nel    = trimesh.n_cell
e2v         = trimesh.cell[1:3,:]
mesh.nv     = maximum(e2v)
e2f         = trimesh.cell[4:6,:] .- mesh.nv
mesh.nf     = maximum(e2f)
mesh.xv     = trimesh.point[1,1:mesh.nv]
mesh.yv     = trimesh.point[2,1:mesh.nv]
mesh.xf     = trimesh.point[1,mesh.nv+1:end]
mesh.yf     = trimesh.point[2,mesh.nv+1:end]
mesh.bc     = trimesh.point_marker[mesh.nv+1:end]

nel  = trimesh.n_cell
vole = zeros(nel)
xc   = zeros(nel)
yc   = zeros(nel)

@avx for iel=1:nel
    # Compute volumes of triangles - use vertices coordinates
    x1 = mesh.xv[e2v[1,iel]]
    y1 = mesh.yv[e2v[1,iel]]
    x2 = mesh.xv[e2v[2,iel]]
    y2 = mesh.yv[e2v[2,iel]]
    x3 = mesh.xv[e2v[3,iel]]
    y3 = mesh.yv[e2v[3,iel]]
    a         = sqrt((x1-x2)^2 + (y1-y2)^2)
    b         = sqrt((x2-x3)^2 + (y2-y3)^2)
    c         = sqrt((x1-x3)^2 + (y1-y3)^2)
    s         = 1/2*(a+b+c)
    vole[iel] = sqrt(s*(s-a)*(s-b)*(s-c))
    xc[iel]   = 1.0/3.0*(x1+x2+x3)
    yc[iel]   = 1.0/3.0*(y1+y2+y3)
end

mesh.e2v    = e2v'
mesh.e2f    = e2f'
mesh.nn_el  = 3
mesh.nf_el  = 3
mesh.xc     = xc
mesh.yc     = yc
mesh.vole   = vole

nodeA = [2 3 1]
nodeB = [3 1 2]
nodeC = [1 2 3]

# Compute normal to faces
mesh.n_x = zeros(Float64,mesh.nel,mesh.nf_el)
mesh.n_y = zeros(Float64,mesh.nel,mesh.nf_el)
mesh.dA  = zeros(Float64,mesh.nel,mesh.nf_el)

 # Assemble FCFV elements
 for iel=1:mesh.nel

    # println("element: ",  iel)

    @avx for ifac=1:mesh.nf_el

        nodei  = mesh.e2f[iel,ifac]

        # Vertices
        vert1  = mesh.e2v[iel,nodeA[ifac]]
        vert2  = mesh.e2v[iel,nodeB[ifac]]
        vert3  = mesh.e2v[iel,nodeC[ifac]]
        bc     = mesh.bc[nodei]
        dx     = (mesh.xv[vert1] - mesh.xv[vert2] );
        dy     = (mesh.yv[vert1] - mesh.yv[vert2] );
        dAi    = sqrt(dx^2 + dy^2);

        # println(bc)
        # @printf("face node, x = %2.2e y = %2.2e\n", mesh.xf[nodei], mesh.yf[nodei])
        # @printf("vert1    , x = %2.2e y = %2.2e\n", mesh.xv[vert1], mesh.yv[vert1])
        # @printf("vert2    , x = %2.2e y = %2.2e\n", mesh.xv[vert2], mesh.yv[vert2])

        # Face normal
        n_x  = -dy/dAi
        n_y  =  dx/dAi

        # Third vector
        v_x  = mesh.xf[nodei] - mesh.xc[iel]
        v_y  = mesh.yf[nodei] - mesh.yc[iel]
        # v_x  = mesh.xv[vert1] - mesh.xv[vert3]
        # v_y  = mesh.yv[vert1] - mesh.yv[vert3]

        # Check wether the normal points outwards
        dot                 = n_x*v_x + n_y*v_y
        mesh.n_x[iel,ifac]  = ((dot>=0.0)*n_x - (dot<0.0)*n_x)
        mesh.n_y[iel,ifac]  = ((dot>=0.0)*n_y - (dot<0.0)*n_y)
        mesh.dA[iel,ifac]   = dAi
    end
end

return mesh

end

#--------------------------------------------------------------------#

function MakeQuadMesh( nx, ny, xmin, xmax, ymin, ymax)
    # Generates a 2D rectangular mesh of nx*ny cells
    ncell      = nx*ny
    # Twice resolution mesh
    nx2 = 2*nx + 1
    ny2 = 2*ny + 1
    # 1D axis
    dx,dy = (xmax-xmin)/nx, (xmax-xmin)/ny
    xc  = collect(LinRange(xmin+dx/2.0,xmax-dx/2.0,nx))
    yc  = collect(LinRange(ymin+dy/2.0,ymax-dy/2.0,ny))
    x2  = LinRange(xmin,xmax,nx2)
    y2  = LinRange(ymin,ymax,ny2)
    # 2D mesh: fake ngrid for P-T space ;)
    x2d   = repeat(x2, 1, length(y2))'
    y2d   = repeat(y2, 1, length(x2))
    nodes = zeros(Int64, size(x2d))
    # Number active dofs (midfaces) - Face node numbering
    inum = 0
    for i=1:nx2
        for j=1:ny2
            if (mod(j,2)==1 && mod(i,2)==0) || (mod(j,2)==0 && mod(i,2)==1)
                inum += 1
                nodes[i,j] = inum
            end
        end
    end
    nface = inum;
    # Number vertices
    inum = 0
    for i=1:nx2
        for j=1:ny2
            if (mod(j,2)==1 && mod(i,2)==1) || (mod(j,2)==1 && mod(i,2)==1)
                inum += 1
                nodes[i,j] = inum
            end
        end
    end
    # Face nodes
    xf = zeros(nface)
    yf = zeros(nface)
    tf = zeros(Int64,nface) # node type - 1 boundary
    for i=1:nx2
        for j=1:ny2
            if (mod(j,2)==1 && mod(i,2)==0) || (mod(j,2)==0 && mod(i,2)==1)
                xf[nodes[i,j]] = x2d[i,j]
                yf[nodes[i,j]] = y2d[i,j]
                if (i==1 || i==nx2 || j==1 || j==ny2)
                    tf[nodes[i,j]] = 1
                end
            end
        end
    end
    # Vertices nodes
    xn = zeros( (nx+1)*(ny+1) )
    yn = zeros( (nx+1)*(ny+1) )
    tn = zeros(Int64, (nx+1)*(ny+1) ) # node type - 1 boundary
    for i=1:nx+1
        for j=1:ny+1
            ii = i + (j-1)*(ny+1)
                xn[ii] = xmin + (i-1)*dx
                yn[ii] = ymin + (j-1)*dy
                if (i==1 || i==nx+1 || j==1 || j==ny+1)
                    tn[ii] = 1
                end
        end
    end
    # Cell 2 face numbering - for matrix connectivity
    face = zeros(Int64, 4, ncell)
    @avx for i=1:nx
        for j=1:ny
            k  = j + (i-1)*ny
            jc = i+1 + (i-1)*1
            ic = j+1 + (j-1)*1
            face[1,k] =  nodes[ic,jc-1] # South
            face[2,k] =  nodes[ic-1,jc] # West
            face[3,k] =  nodes[ic+1,jc] # East
            face[4,k] =  nodes[ic,jc+1] # North
            # println(face[:,k])
        end
    end
    # Cell 2 vertices - used for visualisation of quads
    vert = zeros(Int64, 4, ncell)
    @avx for i=1:nx
        for j=1:ny
            k  = j + (i-1)*ny
            jc = i+1 + (i-1)*1
            ic = j+1 + (j-1)*1
            vert[1,k] =  nodes[ic-1,jc+1] # South
            vert[2,k] =  nodes[ic-1,jc-1] # East
            vert[3,k] =  nodes[ic+1,jc-1] # North
            vert[4,k] =  nodes[ic+1,jc+1] # West
        end
    end
    # Centroids
    xc = zeros(ncell)
    yc = zeros(ncell)
    w  = 0.25
    @avx for iel=1:ncell
        tempx = 0
        tempy = 0
        for j=1:4
            tempx += w * xf[face[j,iel]]
            tempy += w * yf[face[j,iel]]
        end
        xc[iel] = tempx
        yc[iel] = tempy
    end

    # Fill structure
    mesh        = FCFV_Mesh()
    mesh.type   = "quad"
    mesh.nel    = ncell
    mesh.nf     = nface
    mesh.nf_el  = 4
    mesh.nf_el  = 4
    mesh.nv     = (nx+1)*(ny+1)
    mesh.xv     = xn
    mesh.yv     = yn
    mesh.xf     = xf
    mesh.yf     = yf
    mesh.e2v    = vert' # cell nodes - not dofs! Dofs are in 'cell'
    mesh.e2f    = face'
    mesh.xc     = xc
    mesh.yc     = yc
    mesh.vole   = dx*dy*ones(ncell)
    mesh.bc     = tf

    nodeA = [2 1 3 4]
    nodeB = [3 2 4 1]

    # Compute normal to faces
    mesh.n_x = zeros(Float64,mesh.nel,mesh.nf_el)
    mesh.n_y = zeros(Float64,mesh.nel,mesh.nf_el)
    mesh.dA  = zeros(Float64,mesh.nel,mesh.nf_el)

     # Assemble FCFV elements
    @avx for iel=1:mesh.nel

        # println("element: ",  iel)

        for ifac=1:mesh.nf_el

            nodei  = mesh.e2f[iel,ifac]

            # Vertices
            vert1  = mesh.e2v[iel,nodeA[ifac]]
            vert2  = mesh.e2v[iel,nodeB[ifac]]
            bc     = mesh.bc[nodei]
            dx     = abs(mesh.xv[vert1] - mesh.xv[vert2] );
            dy     = abs(mesh.yv[vert1] - mesh.yv[vert2] );
            dAi    = sqrt(dx^2 + dy^2);

            # if iel==1
            #     println(bc)
            #     @printf("face node, x = %2.2e y = %2.2e\n", mesh.xf[nodei], mesh.yf[nodei])
            #     @printf("vert1    , x = %2.2e y = %2.2e\n", mesh.xv[vert1], mesh.yv[vert1])
            #     @printf("vert2    , x = %2.2e y = %2.2e\n", mesh.xv[vert2], mesh.yv[vert2])
            # end

            # Face normal
            n_x  =  dy/dAi
            n_y  = -dx/dAi

            # Third vector
            v_x  = mesh.xf[nodei] - mesh.xc[iel]
            v_y  = mesh.yf[nodei] - mesh.yc[iel]

            # Check wether the normal points outwards
            dot                 = n_x*v_x + n_y*v_y
            mesh.n_x[iel,ifac]  = (dot>=0.0)*n_x - (dot<0.0)*n_x
            mesh.n_y[iel,ifac]  = (dot>=0.0)*n_y - (dot<0.0)*n_y
            mesh.dA[iel,ifac]   = dAi
        end
    end

    return mesh
end


text:
using Units
using Base.Test

l = Quantity(Micro, Meter, 380)
@test prefix(l) == Micro
@test base(l) == Meter
pshow(l)
println()
lshow(l)
println()
fshow(l)
println()
lcm = convert(Unit(Centi,Meter), l)
@test lcm.value == l.value/10^4

temp = parse_quantity("32F")
tempC = convert(Unit(SINone,Celsius), temp)
println(temp, " = ", tempC)

q1 = parse_quantity("3.2 μm")
q2 = parse_quantity("5.7 m")
@test q2/q1 == (5.7/3.2)*1e6
@test 5*q1 == Quantity(Micro, Meter, 5*3.2)
@test q2/7 == Quantity(SINone, Meter, 5.7/7)

q1 = parse_quantity("7 s")
q2 = parse_quantity("4 d")
@test q1/q2 == 7/4/86400

q1 = parse_quantity("1 g")
q2 = parse_quantity("1 Da")
@test q2/q1 == 1.66053892173e-24

q1 = parse_quantity("3.2 μm")
q2 = parse_quantity("4 d")
@test_throws q1/q2


text:
function read_output(file_path::String)
	file_stream::IOStream = open(file_path, "r")
	file_string::String = readstring(file_stream)
	close(file_stream)

	D_empirical::Float64 = read_key(file_string, "D_empirical", Float64)
	t::Array{Float64, 1} = read_key(file_string, "t", Array{Float64, 1})
	msd::Array{Float64, 1} = read_key(file_string, "msd", Array{Float64, 1})
	msd_x::Array{Float64, 1} = read_key(file_string, "msd_x", Array{Float64, 1})
	msd_y::Array{Float64, 1} = read_key(file_string, "msd_y", Array{Float64, 1})
	msd_z::Array{Float64, 1} = read_key(file_string, "msd_z", Array{Float64, 1})
	t_exec::Float64 = read_key(file_string, "t_exec", Float64)

	return (
		D_empirical,
		t,
		msd,
		msd_x,
		msd_y,
		msd_z,
		t_exec)
end


text:
using Images, MXNet, JLD

nnet = mx.load_checkpoint("weights/mobilenet-v2/MobiletNet-FE", 0, mx.FeedForward; context = mx.gpu());

PROJECT_IMG_FOLDER = joinpath("data", "book-covers");
DATASETS = ["train"];
BATCH_SIZE = 100;

if ~isdir(PROJECT_IMG_FOLDER) mkdir(PROJECT_IMG_FOLDER) end;

for dataset_name in DATASETS

    categories = Int[]

    current_folder = joinpath(PROJECT_IMG_FOLDER, dataset_name)
    file_name = joinpath(PROJECT_IMG_FOLDER, "bookcover30-labels-$dataset_name.txt")
    lines = readlines(file_name);

    results = zeros(Float16, (1280, size(lines, 1)))

    for i=1:BATCH_SIZE:size(lines, 1)

        println("BATCH: $i")

        mx_data = mx.zeros((224, 224, 3, BATCH_SIZE));

        for j=1:BATCH_SIZE

            if j % 50 == 0
                println("\tIMG: $j")
            end

            img = imresize(load(joinpath(current_folder, split(lines[i + j - 1])[1])), (224, 224));
            img = permutedims(Float16.(channelview(RGB.(img))), (3, 2, 1));

            img[:, :, :] *= 256.0;
            img[:, :, :] -= 128.;
            img[:, :, :] /= 128.;

            img = reshape(img, 224, 224, 3, 1);
            mx_data[j:j] = img;
        end

        data_provider = mx.ArrayDataProvider(:data => mx_data);
        results[:, i:i+(BATCH_SIZE-1)] = mx.predict(nnet, data_provider)
    end

    JLD.save(joinpath(PROJECT_IMG_FOLDER, "dump-$dataset_name.jld"), "data", results)
end


text:
using NCDatasets
using Formatting
using ArgParse
using Statistics
using JSON

include("./lib/map_transform.jl")
include("constants.jl")
include("CESMReader.jl")

using .CESMReader
using .MapTransform

function mreplace(x)
    return replace(x, missing=>NaN)
end

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file-prefix"
            help = "Data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--data-file-timestamp-form"
            help = "Data filename timestamp form. Either `YEAR` or `YEAR_MONTH`."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--domain-file"
            help = "Ocn domain file."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["domain-file"], "r") do ds
    global area = mreplace(ds["area"][:])
    global mask = mreplace(ds["mask"][:])
    global lat  = mreplace(ds["yc"][:])

    area .*= 4π * Re^2 / sum(area)
end

lat_bnd = collect(Float64, -90:1:90) # 181 elements

r = MapTransform.Relation(
    lat = lat,
    area = area,
    mask = mask,
    lat_bnd = lat_bnd,
)

r_aqua = MapTransform.Relation(
    lat = lat,
    area = area,
    mask = mask * 0 .+ 1,
    lat_bnd = lat_bnd,
)

_1flux = area * 0 .+ 1.0
sum_valid_area = MapTransform.∫∂a(r, _1flux)[end]

println("Sum of valid area: ", sum_valid_area, "; ratio: ", sum_valid_area / sum(area))

ocn_frac = MapTransform.f∂a(r, _1flux) ./ MapTransform.f∂a(r_aqua, _1flux)

let
    global TFLUX_DIV_implied, OHT, TSAS_clim, OHT_TSAS_clim, TSAS_clim_mean

    if parsed["data-file-timestamp-form"] == "YEAR"
        filename_format = format("{:s}{{:04d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR
    elseif parsed["data-file-timestamp-form"] == "YEAR_MONTH"
        filename_format = format("{:s}{{:04d}}-{{:02d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR_MONTH
    end

    fh = FileHandler(filename_format=filename_format, form=form)

    beg_t = (parsed["beg-year"] - 1) * 12 + 1
    end_t = (parsed["end-year"] - 1) * 12 + 12

    Nt = end_t - beg_t + 1

    TFLUX_DIV_implied = zeros(Float64, length(r.lat_bnd)-1, Nt)
    OHT               = zeros(Float64, length(r.lat_bnd)  , Nt)
    TSAS_clim         = zeros(Float64, length(r.lat_bnd)-1, Nt)
    OHT_TSAS_clim     = zeros(Float64, length(r.lat_bnd)  , Nt)
    TSAS_clim_mean    = zeros(Float64, Nt)

    TFLUX_DIV, seaice_nudge_energy, _TSAS_clim = getData(fh, ["TFLUX_DIV_implied", "seaice_nudge_energy", "TSAS_clim"], (parsed["beg-year"], parsed["end-year"]), (:, :))

    adjust_TFLUX_DIV = TFLUX_DIV * ρc + seaice_nudge_energy
    for t = 1:Nt
        _data = view(adjust_TFLUX_DIV, :, :, t)
        TFLUX_DIV_implied[:, t] = MapTransform.transform(r, _data)
        OHT[:, t] = MapTransform.∫∂a(r, _data)
    end

    adjust_TSAS_clim = _TSAS_clim * ρc
    for t = 1:Nt
        _data = view(adjust_TSAS_clim, :, :, t)
        TSAS_clim[:, t] = MapTransform.transform(r, _data)
        TSAS_clim_mean[t] = MapTransform.mean(r, _data)
        OHT_TSAS_clim[:, t] = - MapTransform.∫∂a(r, _data .- TSAS_clim_mean[t])
    end

end

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "time", Inf)
    defDim(ds, "lat_bnd", length(r.lat_bnd))
    defDim(ds, "lat",     length(r.lat_bnd)-1)

    ds.attrib["TSAS_clim_mean_MEAN"] = mean(TSAS_clim_mean)

    for (varname, vardata, vardim, attrib) in [
        ("ocn_frac",          ocn_frac,          ("lat_bnd", ), Dict()),
        ("TFLUX_DIV_implied", TFLUX_DIV_implied, ("lat", "time", ), Dict()),
        ("OHT",               OHT,               ("lat_bnd", "time", ), Dict()),
        ("TSAS_clim",         TSAS_clim,         ("lat", "time", ), Dict()),
        ("TSAS_clim_mean",    TSAS_clim_mean,    ("time", ), Dict()),
        ("OHT_TSAS_clim",     OHT_TSAS_clim,     ("lat_bnd", "time", ), Dict()),
        ("OHT_TOTAL",         OHT + OHT_TSAS_clim,     ("lat_bnd", "time", ), Dict()),

        ("TFLUX_DIV_implied_MEAN", mean(TFLUX_DIV_implied, dims=2)[:, 1],    ("lat",), Dict()),
        ("OHT_MEAN",               mean(OHT, dims=2)[:, 1],                  ("lat_bnd",), Dict()),
        ("TSAS_clim_MEAN",         mean(TSAS_clim, dims=2)[:, 1],            ("lat",), Dict()),
        ("OHT_TSAS_clim_MEAN",     mean(OHT_TSAS_clim, dims=2)[:, 1],        ("lat_bnd",), Dict()),
        ("OHT_TOTAL_MEAN",         mean(OHT + OHT_TSAS_clim, dims=2)[:, 1],  ("lat_bnd",), Dict()),
    ]

        println("Doing var: ", varname)

        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""PriorityLevelConfiguration represents the configuration of a priority level.

    IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration(;
        apiVersion=nothing,
        kind=nothing,
        metadata=nothing,
        spec=nothing,
        status=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta : `metadata` is the standard object's metadata. More info: https://git.k8s.io/community/contributors/devel/api-conventions.md#metadata
    - spec::IoK8sApiFlowcontrolV1alpha1PriorityLevelConfigurationSpec : `spec` is the specification of the desired behavior of a \"request-priority\". More info: https://git.k8s.io/community/contributors/devel/api-conventions.md#spec-and-status
    - status::IoK8sApiFlowcontrolV1alpha1PriorityLevelConfigurationStatus : `status` is the current status of a \"request-priority\". More info: https://git.k8s.io/community/contributors/devel/api-conventions.md#spec-and-status
"""
mutable struct IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1PriorityLevelConfigurationSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1PriorityLevelConfigurationStatus } # spec name: status

    function IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration

const _property_map_IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiFlowcontrolV1alpha1PriorityLevelConfigurationSpec", Symbol("status")=>"IoK8sApiFlowcontrolV1alpha1PriorityLevelConfigurationStatus")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration)
    true
end

function validate_property(::Type{ IoK8sApiFlowcontrolV1alpha1PriorityLevelConfiguration }, name::Symbol, val)
end