text:
"""
`module  JAC.Cascade
    ... a submodel of JAC that contains all methods to set-up and process (simple) atomic cascade computations.
"""
module Cascade

    using Dates, JLD2, Printf, ..AngularMomentum, ..AutoIonization, ..Basics, ..Bsplines, ..Continuum, ..Defaults,
                               ..DecayYield, ..Dielectronic, ..ElectronCapture, ..Radial, ..ManyElectron, ..Nuclear,
                               ..PhotoEmission, ..PhotoExcitation, ..PhotoIonization, ..Semiempirical, ..TableStrings

    """
    `abstract type Cascade.AbstractCascadeScheme`
        ... defines an abstract type to distinguish different excitation, ionization and decay schemes of an atomic cascade; see also:

        + struct StepwiseDecayScheme
            ... to represent a standard decay scheme, starting from the levels of one or several initial multiplets.
        + struct PhotonIonizationScheme
            ... to represent a (prior) photoionization part of a cascade for a list of given photon energies.
        + struct PhotonExcitationScheme
            ... to represent a (prior) photo-excitation part of a cascade for a specified set of excitations.
        + struct DielectronicCaptureScheme
            ... to represent the (dielectronic) capture of an electron for a maximum excitation energy and a given list of
                subshells; the excitation energy refers to the levels of the reference configurations, and cascade blocks are built
                only with the given subshells.
        + struct HollowIonScheme
            ... to represent the capture of one or several electrons into a list of subshells; various distributions of the
                electron among these shells are supported. For the subsequent decay, the list of decay shells need to be specified
                as well.
        + struct ElectronExcitationScheme
            ... to represent a (prior) electron-impact excitation part of a cascade for given electron energies (not yet).
    """
    abstract type  AbstractCascadeScheme       end

    """
    `struct  Cascade.StepwiseDecayScheme  <:  Cascade.AbstractCascadeScheme`
        ... a struct to represent (and generate) a mean-field orbital basis.

        + processes             ::Array{Basics.AbstractProcess,1}
            ... List of the atomic processes that are supported and should be included into the cascade.
        + maxElectronLoss       ::Int64
            ... (Maximum) Number of electrons in which the initial- and final-state configurations can differ from each other;
                this also determines the maximal steps of any particular decay path.
        + chargeStateShifts     ::Dict{Int64,Float64}
            ... (N => en) total energy shifts of all levels with N electrons; these shifts [in a.u.] help open/close decay
                channels by simply shifting the total energies of all levels.
        + NoShakeDisplacements  ::Int64
            ... Maximum number of electron displacements due to shake-up  or shake-down processes in any individual step of cascade.
        + shakeFromShells       ::Array{Shell,1}        ... List of shells from which shake transitions may occur.
        + shakeToShells         ::Array{Shell,1}        ... List of shells into which shake transitions may occur.
    """
    struct   StepwiseDecayScheme  <:  Cascade.AbstractCascadeScheme
        processes               ::Array{Basics.AbstractProcess,1}
        maxElectronLoss         ::Int64
        chargeStateShifts       ::Dict{Int64,Float64}
        NoShakeDisplacements    ::Int64
        shakeFromShells         ::Array{Shell,1}
        shakeToShells           ::Array{Shell,1}
    end

    """
    `Cascade.StepwiseDecayScheme()`  ... constructor for an 'default' instance of a Cascade.StepwiseDecayScheme.
    """
    function StepwiseDecayScheme()
        StepwiseDecayScheme([Radiative()], 0, Dict{Int64,Float64}(), 0, Shell[], Shell[] )
    end

    # `Base.string(scheme::StepwiseDecayScheme)`  ... provides a String notation for the variable scheme::StepwiseDecayScheme.
    function Base.string(scheme::StepwiseDecayScheme)
        sa = "Stepwise decay (scheme) of an atomic cascade with:"
        return( sa )
    end

    # `Base.show(io::IO, scheme::StepwiseDecayScheme)`  ... prepares a proper printout of the scheme::StepwiseDecayScheme.
    function Base.show(io::IO, scheme::StepwiseDecayScheme)
        sa = Base.string(scheme);                print(io, sa, "\n")
        println(io, "processes:                  $(scheme.processes)  ")
        println(io, "maxElectronLoss:            $(scheme.maxElectronLoss)  ")
        println(io, "chargeStateShifts:          $(scheme.chargeStateShifts)  ")
        println(io, "NoShakeDisplacements:       $(scheme.NoShakeDisplacements)  ")
        println(io, "shakeFromShells:            $(scheme.shakeFromShells)  ")
        println(io, "shakeToShells:              $(scheme.shakeToShells)  ")
    end

    """
    `struct  Cascade.PhotonIonizationScheme  <:  Cascade.AbstractCascadeScheme`
        ... a struct to represent (and generate) a mean-field orbital basis.

        + processes             ::Array{Basics.AbstractProcess,1}
            ... List of the atomic processes that are supported and should be included into the cascade.
        + maxPhotoElectrons     ::Int64
            ... (Maximum) Number of photo-electrons in which the initial-and (photo-ionized) final-state configurations can differ
                from each other; this affects the number of ionized multiplets in the cascade.
        + photonEnergies        ::Array{Float64,1}
            ... List of photon energies for which this photo-ionization scheme is to be calculated.
    """
    struct   PhotonIonizationScheme  <:  Cascade.AbstractCascadeScheme
        processes               ::Array{Basics.AbstractProcess,1}
        maxPhotoElectrons       ::Int64
        photonEnergies          ::Array{Float64,1}
    end

    """
    `Cascade.PhotonIonizationScheme()`  ... constructor for an 'default' instance of a Cascade.PhotonIonizationScheme.
    """
    function PhotonIonizationScheme()
        PhotonIonizationScheme([PhotoIonization], 0, Float64[] )
    end

    # `Base.string(scheme::PhotonIonizationScheme)`  ... provides a String notation for the variable scheme::PhotonIonizationScheme.
    function Base.string(scheme::PhotonIonizationScheme)
        sa = "Photo-ionization (scheme) as prior part of an atomic cascade:"
        return( sa )
    end

    # `Base.show(io::IO, scheme::PhotonIonizationScheme)`  ... prepares a proper printout of the scheme::PhotonIonizationScheme.
    function Base.show(io::IO, scheme::PhotonIonizationScheme)
        sa = Base.string(scheme);                print(io, sa, "\n")
        println(io, "processes:                  $(scheme.processes)  ")
        println(io, "maxPhotoElectrons:          $(scheme.maxPhotoElectrons)  ")
        println(io, "photonEnergies:             $(scheme.photonEnergies)  ")
    end

    """
    `struct  Cascade.PhotonExcitationScheme  <:  Cascade.AbstractCascadeScheme`
        ... a struct to define and describe a photo-excitation calculation for an atom in some initial state/configuration
            and for a given set of shell-excitations

        + processes             ::Array{Basics.AbstractProcess,1}
            ... List of the atomic processes that are supported and should be included into the cascade.
        + multipoles            ::Array{EmMultipole}
            ... Multipoles of the radiation field that are to be included into the excitation processes.
        + minPhotonEnergy       ::Float64
            ... Minimum photon energy [in a.u.] that restrict the number of excited configurations to be taken into accout.
        + maxPhotonEnergy       ::Float64
            ... Maximum photon energy [in a.u.] that restrict the number of excited configurations to be taken into accout.
        + NoExcitations         ::Int64
            ... (Maximum) Number of electron replacements with regard to the initial configurations/multiplets.
        + excitationFromShells  ::Array{Shell,1}
            ... List of shells from which photo-excitations are to be considered.
        + excitationToShells    ::Array{Shell,1}
            ... List of shells into which photo-excitations are to be considered, including possibly already occupied shells.
    """
    struct   PhotonExcitationScheme  <:  Cascade.AbstractCascadeScheme
        processes               ::Array{Basics.AbstractProcess,1}
        multipoles              ::Array{EmMultipole}
        minPhotonEnergy         ::Float64
        maxPhotonEnergy         ::Float64
        NoExcitations           ::Int64
        excitationFromShells    ::Array{Shell,1}
        excitationToShells      ::Array{Shell,1}
    end

    """
    `Cascade.PhotonExcitationScheme()`  ... constructor for an 'default' instance of a Cascade.PhotonExcitationScheme.
    """
    function PhotonExcitationScheme()
        PhotonExcitationScheme([PhotoExc()], [E1], 1.0, 100., 0, Shell[], Shell[] )
    end

    # `Base.string(scheme::PhotonExcitationScheme)`  ... provides a String notation for the variable scheme::PhotoAbsorptionScheme.
    function Base.string(scheme::PhotonExcitationScheme)
        sa = "Photon-excitation (scheme) due to processes:"
        return( sa )
    end

    # `Base.show(io::IO, scheme::PhotonExcitationScheme)`  ... prepares a proper printout of the scheme::PhotonExcitationScheme.
    function Base.show(io::IO, scheme::PhotonExcitationScheme)
        sa = Base.string(scheme);                print(io, sa, "\n")
        println(io, "processes:                  $(scheme.processes)  ")
        println(io, "multipoles:                 $(scheme.multipoles)  ")
        println(io, "minPhotonEnergy:            $(scheme.minPhotonEnergy)  ")
        println(io, "maxPhotonEnergy:            $(scheme.maxPhotonEnergy)  ")
        println(io, "NoExcitations :             $(scheme.NoExcitations )  ")
        println(io, "excitationFromShells:       $(scheme.excitationFromShells)  ")
        println(io, "excitationToShells:         $(scheme.excitationToShells)  ")
    end

    """
    `struct  Cascade.DielectronicCaptureScheme  <:  Cascade.AbstractCascadeScheme`
        ... a struct to define and describe the (dielectronic) capture of electrons for an atom in some initial state/configuration;
            for such a scheme, the doubly-excited configurations due to the electron capture are generated automatically due to
            given maximal numbers of the (into-) shells (nl) as well as the maximum displacement with regard to the initial configuration.
            An additional maximum excitation energy need to be provided due to the maximum temperatures for which DR plasma rate coefficients
            are to be determined, cf. Basics.convert().

        + multipoles            ::Array{EmMultipole}
            ... Multipoles of the radiation field that are to be included into the radiative stabilization processes.
        + maxExcitationEnergy   ::Float64
            ... Maximum excitation energy [in a.u.] with regard to the reference configurations/levels that restrict the number
                of excited configurations to be taken into accout. This maximum excitation energy has to be derived from the maximum
                temperature for which DR coefficients need to be derived and is typically set to 5x T_e,max.
        + electronEnergyShift   ::Float64
            ... Energy shift for all resonance energies; this is realized by shifting the initial level energies by the negative amount.
                The shift is taken in the user-defined units.
        + minPhotonEnergy       ::Float64
            ... Minimum (mean) photon energy [in a.u.] in the radiative stabilization of the doubly-excited configurations;
                If cascade blocks are separated be less than this energy, the radiative stabilization is neglected.
        + NoExcitations         ::Int64
            ... (Maximum) Number of electron replacements in the doubly-excited configuration with regard to the initial
                configurations/multiplets, apart from one additional electron due to the electron capture itself.
        + excitationFromShells  ::Array{Shell,1}
            ... List of shells from which excitations are to be considered.
        + excitationToShells  ::Array{Shell,1}
            ... List of shells to which (core-shell) excitations are to be considered.
        + intoShells            ::Array{Shell,1}
            ... List of shells into which electrons are initially placed (captured).
        + decayShells           ::Array{Shell,1}
            ... List of shells into which electrons the electrons can decay (apart from the core shells).
    """
    struct   DielectronicCaptureScheme  <:  Cascade.AbstractCascadeScheme
        multipoles              ::Array{EmMultipole}
        maxExcitationEnergy     ::Float64
        electronEnergyShift     ::Float64
        minPhotonEnergy         ::Float64
        NoExcitations           ::Int64
        excitationFromShells    ::Array{Shell,1}
        excitationToShells      ::Array{Shell,1}
        intoShells              ::Array{Shell,1}
        decayShells             ::Array{Shell,1}
    end

    """
    `Cascade.DielectronicCaptureScheme()`  ... constructor for an 'default' instance of a Cascade.DielectronicCaptureScheme.
    """
    function DielectronicCaptureScheme()
        DielectronicCaptureScheme([E1], 1.0, 0., 0., 1, Shell[], Shell[], Shell[], Shell[] )
    end

    # `Base.string(scheme::DielectronicCaptureScheme)`  ... provides a String notation for the variable scheme::DielectronicCaptureScheme.
    function Base.string(scheme::DielectronicCaptureScheme)
        sa = "Dielectronic capture & stabilization (scheme):"
        return( sa )
    end

    # `Base.show(io::IO, scheme::DielectronicCaptureScheme)`  ... prepares a proper printout of the scheme::DielectronicCaptureScheme.
    function Base.show(io::IO, scheme::DielectronicCaptureScheme)
        sa = Base.string(scheme);                print(io, sa, "\n")
        println(io, "multipoles:                 $(scheme.multipoles)  ")
        println(io, "maxExcitationEnergy:        $(scheme.maxExcitationEnergy)  ")
        println(io, "electronEnergyShift:        $(scheme.electronEnergyShift)  ")
        println(io, "minPhotonEnergy:            $(scheme.minPhotonEnergy)  ")
        println(io, "NoExcitations:              $(scheme.NoExcitations)  ")
        println(io, "excitationFromShells:       $(scheme.excitationFromShells)  ")
        println(io, "excitationToShells:         $(scheme.excitationToShells)  ")
        println(io, "intoShells:                 $(scheme.intoShells)  ")
        println(io, "decayShells:                $(scheme.decayShells)  ")
    end

    """
    `struct  Cascade.HollowIonScheme  <:  Cascade.AbstractCascadeScheme`
        ... a struct to define and describe the formation and decay of a hollow ion, e.g. an electronic core configuration
            into which one or several additional electrons are captured into (high) nl shells. Both the shell (lists) for
            the initial capture (intoShells) and the subsequent decay (decayShells) need to be specified explicitly to
            readily control the size of the computations.

        + processes             ::Array{Basics.AbstractProcess,1}
            ... List of the atomic processes that are supported and should be included into the decay scheme.
        + multipoles            ::Array{EmMultipole,1}
            ... Multipoles of the radiation field that are to be included into the radiative stabilization processes.
        + NoCapturedElectrons   ::Int64
            ... Number of captured electrons, e.g. placed in the intoShells.
        + intoShells            ::Array{Shell,1}
            ... List of shells into which electrons are initially placed (captured).
        + decayShells           ::Array{Shell,1}
            ... List of shells into which electrons the electrons can decay (apart from the core shells).
    """
    struct   HollowIonScheme  <:  Cascade.AbstractCascadeScheme
        processes               ::Array{Basics.AbstractProcess,1}
        multipoles              ::Array{EmMultipole}
        NoCapturedElectrons     ::Int64
        intoShells              ::Array{Shell,1}
        decayShells             ::Array{Shell,1}
    end

    """
    `Cascade.HollowIonScheme()`  ... constructor for an 'default' instance of a Cascade.HollowIonScheme.
    """
    function HollowIonScheme()
        HollowIonScheme([Auger(), Radiative()], [E1], 0, Shell[], Shell[] )
    end

    # `Base.string(scheme::HollowIonScheme)`  ... provides a String notation for the variable scheme::HollowIonScheme.
    function Base.string(scheme::HollowIonScheme)
        sa = "Hollow ion (scheme) due to processes:"
        return( sa )
    end

    # `Base.show(io::IO, scheme::HollowIonScheme)`  ... prepares a proper printout of the scheme::HollowIonScheme.
    function Base.show(io::IO, scheme::HollowIonScheme)
        sa = Base.string(scheme);                print(io, sa, "\n")
        println(io, "processes:                  $(scheme.processes)  ")
        println(io, "multipoles:                 $(scheme.multipoles)  ")
        println(io, "NoCapturedElectrons:        $(scheme.NoCapturedElectrons)  ")
        println(io, "intoShells:                 $(scheme.intoShells)  ")
        println(io, "decayShells:                $(scheme.decayShells)  ")
    end

    """
    `struct  Cascade.ElectronExcitationScheme  <:  Cascade.AbstractCascadeScheme`
        ... a struct to represent (and generate) a mean-field orbital basis.

        + processes             ::Array{Basics.AbstractProcess,1}
            ... List of the atomic processes that are supported and should be included into the cascade.
        + electronEnergies      ::Array{Float64,1}
            ... List of electron energies for which this electron-impact excitation scheme is to be calculated.
    """
    struct   ElectronExcitationScheme  <:  Cascade.AbstractCascadeScheme
        processes               ::Array{Basics.AbstractProcess,1}
        electronEnergies        ::Array{Float64,1}
    end

    """
    `Cascade.ElectronExcitationScheme()`  ... constructor for an 'default' instance of a Cascade.ElectronExcitationScheme.
    """
    function ElectronExcitationScheme()
        ElectronExcitationScheme([Eimex], Float64[] )
    end

    # `Base.string(scheme::ElectronExcitationScheme)`  ... provides a String notation for the variable scheme::ElectronExcitationScheme.
    function Base.string(scheme::ElectronExcitationScheme)
        sa = "Electron-impact excitation (scheme) as prior part of an atomic cascade:"
        return( sa )
    end

    # `Base.show(io::IO, scheme::ElectronExcitationScheme)`  ... prepares a proper printout of the scheme::ElectronExcitationScheme.
    function Base.show(io::IO, scheme::ElectronExcitationScheme)
        sa = Base.string(scheme);                 print(io, sa, "\n")
        println(io, "processes:                   $(scheme.processes)  ")
        println(io, "electronEnergies:            $(scheme.electronEnergies)  ")
    end

    """
    `abstract type Cascade.AbstractCascadeApproach`
        ... defines an abstract and a number of singleton types for the computational approach/model that is applied in order to
            generate and evaluate all many-electron amplitudes of a given cascade.

      + struct AverageSCA
        ... all levels in the cascade are described in single-configuration and single-CSF approximation; this (rather crude) approach
            neglects all configuration-interactions and also applies just a single set of one-electron orbitals (from the least-ionized charge
            state) for all considered charge states.

      + struct SCA
        ... all levels in the cascade are described in single-configuration approximation but with 'mixtures' within the configuration;
            an individual mean-field is generated for each charge state and all continuum orbitals are generated for the correct transition
            energy in the field of the remaining ion. Moreover, all the fine-structure transitions are calculated individually.
    """
    abstract type  AbstractCascadeApproach                   end
    struct         AverageSCA  <:  AbstractCascadeApproach   end
    struct         SCA         <:  AbstractCascadeApproach   end
    struct         UserMCA     <:  AbstractCascadeApproach   end

    """
    `struct  Cascade.Block`
        ... defines a type for an individual block of configurations that are treatet together within the cascade. Such an block is given
            by a list of configurations that may occur as initial- and/or final-state configurations in some step of the canscade and that
            give rise to a common multiplet in order to allow for configuration interactions but to avoid 'double counting' of individual
            levels in the cascade.

        + NoElectrons     ::Int64                     ... Number of electrons in this block.
        + confs           ::Array{Configuration,1}    ... List of one or several configurations that define the multiplet.
        + hasMultiplet    ::Bool
            ... true if the (level representation in the) multiplet has already been computed and false otherwise.
        + multiplet       ::Multiplet                 ... Multiplet of the this block.
    """
    struct  Block
        NoElectrons       ::Int64
        confs             ::Array{Configuration,1}
        hasMultiplet      ::Bool
        multiplet         ::Multiplet
    end

    """
    `Cascade.Block()`  ... constructor for an 'empty' instance of a Cascade.Block.
    """
    function Block()
        Block( )
    end

    # `Base.show(io::IO, block::Cascade.Block)`  ... prepares a proper printout of the variable block::Cascade.Block.
    function Base.show(io::IO, block::Cascade.Block)
        println(io, "NoElectrons:        $(block.NoElectrons)  ")
        println(io, "confs :             $(block.confs )  ")
        println(io, "hasMultiplet:       $(block.hasMultiplet)  ")
        println(io, "multiplet:          $(block.multiplet)  ")
    end

    """
    `struct  Cascade.Step`
        ... defines a type for an individual step of an excitation and/or decay cascade. Such a step is determined by the two lists
            of initial- and final-state configuration as well as by the atomic process, such as Auger, PhotoEmission, or others, which related
            the initial- and final-state levels to each other. Since the (lists of) initial- and final-state configurations treated (each)
            by a single multiplet (for parities and total angular momenta), a cascade step supports full configuration interaction within
            the multiplet but also help  avoid 'double counting' of individual levels. Indeed, each electron configuration may occur only in
            one cascade block. In contrast, each list of initial- and final-state multiplets (cascade blocks) can occur in quite different
            steps due to the considered processes and parallel decay pathes in a cascade.

        + process          ::JBasics.AbstractProcess   ... Atomic process that 'acts' in this step of the cascade.
        + settings         ::Union{PhotoEmission.Settings, AutoIonization.Settings, PhotoIonization.Settings, PhotoExcitation.Settings,
                                   DielectronicCapture.Settings}
                                                       ... Settings for this step of the cascade.
        + initialConfigs   ::Array{Configuration,1}    ... List of one or several configurations that define the initial-state multiplet.
        + finalConfigs     ::Array{Configuration,1}    ... List of one or several configurations that define the final-state multiplet.
        + initialMultiplet ::Multiplet                 ... Multiplet of the initial-state levels of this step of the cascade.
        + finalMultiplet   ::Multiplet                 ... Multiplet of the final-state levels of this step of the cascade.
    """
    struct  Step
        process            ::Basics.AbstractProcess
        settings           ::Union{PhotoEmission.Settings, AutoIonization.Settings, PhotoIonization.Settings, PhotoExcitation.Settings}
        initialConfigs     ::Array{Configuration,1}
        finalConfigs       ::Array{Configuration,1}
        initialMultiplet   ::Multiplet
        finalMultiplet     ::Multiplet
    end

    """
    `Cascade.Step()`  ... constructor for an 'empty' instance of a Cascade.Step.
    """
    function Step()
        Step( Basics.NoProcess, PhotoEmission.Settings, Configuration[], Configuration[], Multiplet(), Multiplet())
    end

    # `Base.show(io::IO, step::Cascade.Step)`  ... prepares a proper printout of the variable step::Cascade.Step.
    function Base.show(io::IO, step::Cascade.Step)
        println(io, "process:                $(step.process)  ")
        println(io, "settings:               $(step.settings)  ")
        println(io, "initialConfigs:         $(step.initialConfigs)  ")
        println(io, "finalConfigs:           $(step.finalConfigs)  ")
        println(io, "initialMultiplet :      $(step.initialMultiplet )  ")
        println(io, "finalMultiplet:         $(step.finalMultiplet)  ")
    end

    """
    `struct  Cascade.Computation`
        ... defines a type for a cascade computation, i.e. for the computation of a whole photon excitation, photon ionization and/or
            decay cascade. The -- input and control -- data from this computation can be modified, adapted and refined to the practical needs,
            and before the actual computations are carried out explictly. Initially, this struct just contains the physical meta-data about the
            cascade, that is to be calculated, but a new instance of the same Cascade.Computation gets later enlarged in course of the
            computation in order to keep also wave functions, level multiplets, etc.

        + name               ::String                          ... A name for the cascade
        + nuclearModel       ::Nuclear.Model                   ... Model, charge and parameters of the nucleus.
        + grid               ::Radial.Grid                     ... The radial grid to be used for the computation.
        + asfSettings        ::AsfSettings                     ... Provides the settings for the SCF process.
        + scheme             ::Cascade.AbstractCascadeScheme   ... Scheme of the atomic cascade (photoionization, decay, ...)
        + approach           ::Cascade.AbstractCascadeApproach
            ... Computational approach/model that is applied to generate and evaluate the cascade; possible approaches are:
                {AverageSCA(), SCA(), ...}
        + initialConfs       ::Array{Configuration,1}
            ... List of one or several configurations that contain the level(s) from which the cascade starts.
        + initialMultiplets  ::Array{Multiplet,1}
            ... List of one or several (initial) multiplets; either initialConfs 'xor' initialMultiplets can  be specified
                for a given cascade computation.
    """
    struct  Computation
        name                 ::String
        nuclearModel         ::Nuclear.Model
        grid                 ::Radial.Grid
        asfSettings          ::AsfSettings
        scheme               ::Cascade.AbstractCascadeScheme
        approach             ::Cascade.AbstractCascadeApproach
        initialConfigs       ::Array{Configuration,1}
        initialMultiplets    ::Array{Multiplet,1}
    end

    """
    `Cascade.Computation()`  ... constructor for an 'default' instance of a Cascade.Computation.
    """
    function Computation()
        Computation("Default cascade computation",  Nuclear.Model(10.), Radial.Grid(), AsfSettings(), Cascade.StepwiseDecayScheme(),
                    Cascade.AverageSCA(), Configuration[], Multiplet[] )
    end

    """
    `Cascade.Computation(comp::Cascade.Computation;`

                name=..,               nuclearModel=..,             grid=..,              asfSettings=..,
                scheme=..,             approach=..,                 initialConfigs=..,    initialMultiplets=..)

        ... constructor for re-defining the computation::Cascade.Computation.
    """
    function Computation(comp::Cascade.Computation;
        name::Union{Nothing,String}=nothing,                                  nuclearModel::Union{Nothing,Nuclear.Model}=nothing,
        grid::Union{Nothing,Radial.Grid}=nothing,                             asfSettings::Union{Nothing,AsfSettings}=nothing,
        scheme::Union{Nothing,Cascade.AbstractCascadeScheme}=nothing,         approach::Union{Nothing,Cascade.AbstractCascadeApproach}=nothing,
        initialConfigs::Union{Nothing,Array{Configuration,1}}=nothing,        initialMultiplets::Union{Nothing,Array{Multiplet,1}}=nothing)

        if  initialConfigs != nothing   &&   initialMultiplets != nothing
            error("Only initialConfigs=..  'xor'  initialMultiplets=..  can be specified for a Cascade.Computation().")
        end

        if  name                 == nothing   namex                 = comp.name                    else  namex = name                                   end
        if  nuclearModel         == nothing   nuclearModelx         = comp.nuclearModel            else  nuclearModelx = nuclearModel                   end
        if  grid                 == nothing   gridx                 = comp.grid                    else  gridx = grid                                   end
        if  asfSettings          == nothing   asfSettingsx          = comp.asfSettings             else  asfSettingsx = asfSettings                     end
        if  scheme               == nothing   schemex               = comp.scheme                  else  schemex = scheme                               end
        if  approach             == nothing   approachx             = comp.approach                else  approachx = approach                           end
        if  initialConfigs       == nothing   initialConfigsx       = comp.initialConfigs          else  initialConfigsx = initialConfigs               end
        if  initialMultiplets    == nothing   initialMultipletsx    = comp.initialMultiplets       else  initialMultipletsx = initialMultiplets         end

    	Computation(namex, nuclearModelx, gridx, asfSettingsx, schemex, approachx, initialConfigsx, initialMultipletsx)
    end

    # `Base.string(comp::Cascade.Computation)`  ... provides a String notation for the variable comp::Cascade.Computation.
    function Base.string(comp::Cascade.Computation)
        if        typeof(comp.scheme) == Cascade.StepwiseDecayScheme       sb = "stepwise decay scheme"
        elseif    typeof(comp.scheme) == Cascade.PhotonIonizationScheme    sb = "(prior) photo-ionization scheme"
        elseif    typeof(comp.scheme) == Cascade.PhotonExcitationScheme    sb = "photo-excitation scheme"
        elseif    typeof(comp.scheme) == Cascade.DielectronicCaptureScheme     sb = "(di-) electronic capture scheme"
        elseif    typeof(comp.scheme) == Cascade.HollowIonScheme           sb = "hollow ion scheme"
        else      error("unknown typeof(comp.scheme)")
        end

        sa = "Cascade computation   $(comp.name)  for a $sb,  in $(comp.approach) approach as well as "
        sa = sa * "for Z = $(comp.nuclearModel.Z) and initial configurations: \n "
        if  comp.initialConfigs    != Configuration[]  for  config  in  comp.initialConfigs     sa = sa * string(config) * ",  "     end     end
        if  comp.initialMultiplets != Multiplet[]
            for  mp      in  comp.initialMultiplets  sa = sa * mp.name * " with " * string(length(mp.levels)) * " levels,  "         end     end
        return( sa )
    end

    # `Base.show(io::IO, comp::Cascade.Computation)`  ... prepares a proper printout comp::Cascade.Computation.
    function Base.show(io::IO, comp::Cascade.Computation)
        sa = Base.string(comp)
        sa = sa * "\n ... in addition, the following parameters/settings are defined: ";       print(io, sa, "\n")
        println(io, "> cascade scheme:           $(comp.scheme)  ")
        println(io, "> nuclearModel:             $(comp.nuclearModel)  ")
        println(io, "> grid:                     $(comp.grid)  ")
        println(io, "> asfSettings:            \n$(comp.asfSettings) ")
    end

    #######################################################################################################################################
    #######################################################################################################################################
    #######################################################################################################################################

    """
    abstract type  Cascade.AbstractData`
        ... defines an abstract type to distinguish different output data of a cascade computation; see also:

        + struct DecayData       ... to comprise the amplitudes/rates from a stepswise decay cascade.
        + struct PhotoIonData    ... to comprise the amplitudes/rates from a photo-ionization part of a cascade.
        + struct ExcitationData  ... to comprise the amplitudes/rates from a photo-excitation part of a cascade.
        + struct CaptureData     ... to comprise the amplitudes/rates from a electron-capture part of a cascade.
    """
    abstract type  AbstractData      end

    """
    `struct  Cascade.DecayData  <:  Cascade.AbstractData`  ... defines a type for an atomic cascade, i.e. lists of radiative and Auger lines.

        + linesR         ::Array{PhotoEmission.Line,1}          ... List of radiative lines.
        + linesA         ::Array{AutoIonization.Line,1}         ... List of Auger lines.
    """
    struct  DecayData  <:  Cascade.AbstractData
        linesR           ::Array{PhotoEmission.Line,1}
        linesA           ::Array{AutoIonization.Line,1}
    end

    """
    `Cascade.DecayData()`  ... (simple) constructor for cascade  decay data.
    """
    function DecayData()
        DecayData(Array{PhotoEmission.Line,1}[], Array{AutoIonization.Line,1}[] )
    end

    # `Base.show(io::IO, data::Cascade.DecayData)`  ... prepares a proper printout of the variable data::Cascade.DecayData.
    function Base.show(io::IO, data::Cascade.DecayData)
        println(io, "linesR:                  $(data.linesR)  ")
        println(io, "linesA:                  $(data.linesA)  ")
    end

    """
    `struct  Cascade.PhotoIonData  <:  Cascade.AbstractData`  ... defines a type for an atomic cascade, i.e. lists of photoionization lines.

        + photonEnergy   ::Float64                              ... Photon energy for this (part of the) cascade.
        + linesP         ::Array{PhotoIonization.Line,1}        ... List of photoionization lines.
    """
    struct  PhotoIonData  <:  Cascade.AbstractData
        photonEnergy     ::Float64
        linesP           ::Array{PhotoIonization.Line,1}
    end

    """
    `Cascade.PhotoIonData()`  ... (simple) constructor for cascade photo-ionization data.
    """
    function PhotoIonData()
        PhotoIonData(0., Array{PhotoIonization.Line,1}[] )
    end

    # `Base.show(io::IO, data::Cascade.PhotoIonData)`  ... prepares a proper printout of the variable data::Cascade.PhotoIonData.
    function Base.show(io::IO, data::Cascade.PhotoIonData)
        println(io, "photonEnergy:            $(data.photonEnergy)  ")
        println(io, "linesP:                  $(data.linesP)  ")
    end

    """
    `struct  Cascade.ExcitationData  <:  Cascade.AbstractData`  ... defines a type for an atomic cascade, i.e. lists of excitation lines.

        + linesE         ::Array{PhotoExcitation.Line,1}        ... List of photoionization lines.
    """
    struct  ExcitationData  <:  Cascade.AbstractData
        linesE           ::Array{PhotoExcitation.Line,1}
    end

    """
    `Cascade.ExcitationData()`  ... (simple) constructor for cascade excitation data.
    """
    function ExcitationData()
        ExcitationData(Array{PhotoExcitation.Line,1}[] )
    end

    # `Base.show(io::IO, data::Cascade.ExcitationData)`  ... prepares a proper printout of the variable data::Cascade.ExcitationData.
    function Base.show(io::IO, data::Cascade.ExcitationData)
        println(io, "linesE:                  $(data.linesE)  ")
    end

    #######################################################################################################################################
    #######################################################################################################################################
    #######################################################################################################################################

    """
    `abstract type  Cascade.AbstractSimulationProperty`
        ... defines an abstract and various singleton types for the different properties that can be obtained from the simulation of
            cascade data.

        + struct DecayPathes                ... determine the major 'decay pathes' of the cascade.
        + struct DrRateCoefficients         ... simulate the DR (plasma) rate coefficients for given plasma temperatures.
        + struct ElectronCoincidence        ... simulate electron-coincidence spectra (not yet).
        + struct FinalLevelDistribution     ... simulate the 'final-level distribution' as it is found after all cascade
                                                processes are completed.
        + struct IonDistribution            ... simulate the 'ion distribution' as it is found after all cascade processes are completed.
        + struct MeanLineWidths             ... simulate the mean line widths of a line near to a given energy (not yet).
        + struct MeanRelaxationTime         ... simulate the mean relaxation times in which 70%, 80%, of the occupied levels
                                                decay to the ground configuration.
        + struct ElectronIntensities        ... simulate the electron-line intensities as function of electron energy.
        + struct PhotoAbsorption            ... simulate the photoabsorption cross sections for a given set of photo-excitation
                                                and ionization processes.
        + struct PhotonIntensities          ... simulate the photon-line intensities as function of electron energy.
        + struct TimeBinnedPhotonIntensity  ... simulate the photon-line intensities as function of electron energy
                                                in a given time interval (not yet).
    """
    abstract type  AbstractSimulationProperty                              end

    """
    `struct  Cascade.DrRateCoefficients   <:  Cascade.AbstractSimulationProperty`
        ... defines a type for simulating the DR plasma rate coefficients as function of the (free) electron energy and
            plasma temperature.

        + initialLevelNo      ::Int64       ... Level No of initial level for which rate coefficients are to be computed.
        + electronEnergyShift ::Float64
            ... (total) energy shifts that apply to all resonances when alpha^(DR) is computed.
        + temperatures        ::Array{Float64,1}
            ... temperatures [K] for which the DR plasma rate coefficieints to be calculated.
    """
    struct  DrRateCoefficients   <:  Cascade.AbstractSimulationProperty
        initialLevelNo        ::Int64
        electronEnergyShift   ::Float64
        temperatures          ::Array{Float64,1}
    end

    """
    `Cascade.DrRateCoefficients()`  ... (simple) constructor for cascade DrRateCoefficients.
    """
    function DrRateCoefficients()
        DrRateCoefficients(1, 0.,  Float64[])
    end

    # `Base.show(io::IO, dist::Cascade.DrRateCoefficients)`  ... prepares a proper printout of the variable data::Cascade.DrRateCoefficients.
    function Base.show(io::IO, dist::Cascade.DrRateCoefficients)
        println(io, "initialLevelNo:           $(dist.initialLevelNo)  ")
        println(io, "electronEnergyShift:      $(dist.electronEnergyShift)  ")
        println(io, "temperatures:             $(dist.temperatures)  ")
    end

    """
    `struct  Cascade.FinalLevelDistribution   <:  Cascade.AbstractSimulationProperty`
        ... defines a type for simulating the 'final-level distribution' as it is found after all cascade processes are completed.

        + initialOccupations  ::Array{Tuple{Int64,Float64},1}
            ... List of one or several (tupels of) levels in the overall cascade tree together with their relative population.
    """
    struct  FinalLevelDistribution   <:  Cascade.AbstractSimulationProperty
        initialOccupations    ::Array{Tuple{Int64,Float64},1}
    end

    """
    `struct  Cascade.ElectronIntensities   <:  Cascade.AbstractSimulationProperty`
        ... defines a type for simulating the electron-line intensities as function of electron energy.

        + minElectronEnergy   ::Float64     ... Minimum electron energy for the simulation of electron spectra.
        + maxElectronEnergy   ::Float64     ... Maximum electron energy for the simulation of electron spectra.
        + initialOccupations  ::Array{Tuple{Int64,Float64},1}
            ... List of one or several (tupels of) levels in the overall cascade tree together with their relative population.
    """
    struct  ElectronIntensities   <:  Cascade.AbstractSimulationProperty
        minElectronEnergy     ::Float64
        maxElectronEnergy     ::Float64
        initialOccupations    ::Array{Tuple{Int64,Float64},1}
    end

    """
    `Cascade.ElectronIntensities()`  ... (simple) constructor for cascade ElectronIntensities.
    """
    function ElectronIntensities()
        ElectronIntensities(0., 1.0e6,  [(1, 1.0)])
    end

    # `Base.show(io::IO, data::Cascade.ElectronIntensities)`  ... prepares a proper printout of the variable data::Cascade.ElectronIntensities.
    function Base.show(io::IO, dist::Cascade.ElectronIntensities)
        println(io, "minElectronEnergy:        $(dist.minElectronEnergy)  ")
        println(io, "maxElectronEnergy:        $(dist.maxElectronEnergy)  ")
        println(io, "initialOccupations:       $(dist.initialOccupations)  ")
    end

    """
    `Cascade.FinalLevelDistribution()`  ... (simple) constructor for cascade FinalLevelDistribution.
    """
    function FinalLevelDistribution()
        FinalLevelDistribution([(1, 1.0)])
    end

    # `Base.show(io::IO, data::Cascade.FinalLevelDistribution)`
    #       ... prepares a proper printout of the variable data::Cascade.FinalLevelDistribution.
    function Base.show(io::IO, dist::Cascade.FinalLevelDistribution)
        println(io, "initialOccupations:       $(dist.initialOccupations)  ")
    end

    """
    `struct  Cascade.IonDistribution   <:  Cascade.AbstractSimulationProperty`
        ... defines a type for simulating the 'ion distribution' as it is found after all cascade processes are completed.

        + initialOccupations  ::Array{Tuple{Int64,Float64},1}
            ... List of one or several (tupels of) levels in the overall cascade tree together with their relative population.
        + leadingConfigs      ::Array{Configuration,1}
            ... List of leading configurations whose levels are equally populated, either initially or ....
    """
    struct  IonDistribution   <:  Cascade.AbstractSimulationProperty
        initialOccupations    ::Array{Tuple{Int64,Float64},1}
        leadingConfigs        ::Array{Configuration,1}
    end

    """
    `Cascade.IonDistribution()`  ... (simple) constructor for cascade IonDistribution data.
    """
    function IonDistribution()
        IonDistribution([(1, 1.0)], Configuration[])
    end

    # `Base.show(io::IO, data::Cascade.IonDistribution)`  ... prepares a proper printout of the variable data::Cascade.IonDistribution.
    function Base.show(io::IO, dist::Cascade.IonDistribution)
        println(io, "initialOccupations:       $(dist.initialOccupations)  ")
        println(io, "leadingConfigs:           $(dist.leadingConfigs)  ")
    end

    """
    `struct  Cascade.MeanRelaxationTime   <:  Cascade.AbstractSimulationProperty`
        ... defines a type for simulating the mean relaxation times in which 70%, 80%, of the occupied levels decay to the
            ground configuration.

        + timeStep            ::Float64     ... Time-step for following the decay of levels [a.u.]
        + initialOccupations  ::Array{Tuple{Int64,Float64},1}
            ... List of one or several (tupels of) levels in the overall cascade tree together with their relative population.
        + leadingConfigs      ::Array{Configuration,1}
            ... List of leading configurations whose levels are equally populated, either initially or ....
        + groundConfigs       ::Array{Configuration,1}
            ... List of ground configurations into which the decay is considered.
    """
    struct  MeanRelaxationTime   <:  Cascade.AbstractSimulationProperty
        timeStep              ::Float64
        initialOccupations    ::Array{Tuple{Int64,Float64},1}
        leadingConfigs        ::Array{Configuration,1}
        groundConfigs         ::Array{Configuration,1}
    end

    """
    `Cascade.MeanRelaxationTime()`  ... (simple) constructor for cascade MeanRelaxationTime.
    """
    function MeanRelaxationTime()
        MeanRelaxationTime(10., [(1, 1.0)], Configuration[], Configuration[])
    end

    # `Base.show(io::IO, dist::Cascade.MeanRelaxationTime)`  ... prepares a proper printout of the variable data::Cascade.MeanRelaxationTime.
    function Base.show(io::IO, dist::Cascade.MeanRelaxationTime)
        println(io, "timeStep:                 $(dist.timeStep)  ")
        println(io, "initialOccupations:       $(dist.initialOccupations)  ")
        println(io, "leadingConfigs:           $(dist.leadingConfigs)  ")
        println(io, "groundConfigs:            $(dist.groundConfigs)  ")
    end

    """
    `struct  Cascade.PhotonIntensities   <:  Cascade.AbstractSimulationProperty`
        ... defines a type for simulating the photon-line intensities as function of photon energy.

        + minPhotonEnergy     ::Float64     ... Minimum photon energy for the simulation of photon spectra.
        + maxPhotonEnergy     ::Float64     ... Maximum photon energy for the simulation of photon spectra.
        + initialOccupations  ::Array{Tuple{Int64,Float64},1}
            ... List of one or several (tupels of) levels in the overall cascade tree together with their relative population.
        + leadingConfigs      ::Array{Configuration,1}
            ... List of leading configurations whose levels are equally populated, either initially or ....
    """
    struct  PhotonIntensities   <:  Cascade.AbstractSimulationProperty
        minPhotonEnergy       ::Float64
        maxPhotonEnergy       ::Float64
        initialOccupations    ::Array{Tuple{Int64,Float64},1}
        leadingConfigs        ::Array{Configuration,1}
    end

    """
    `Cascade.PhotonIntensities()`  ... (simple) constructor for cascade PhotonIntensities.
    """
    function PhotonIntensities()
        PhotonIntensities(0., 1.0e6,  [(1, 1.0)], Configuration[])
    end

    # `Base.show(io::IO, dist::Cascade.PhotonIntensities)`  ... prepares a proper printout of the variable data::Cascade.PhotonIntensities.
    function Base.show(io::IO, dist::Cascade.PhotonIntensities)
        println(io, "minPhotonEnergy:          $(dist.minPhotonEnergy)  ")
        println(io, "maxPhotonEnergy:          $(dist.maxPhotonEnergy)  ")
        println(io, "initialOccupations:       $(dist.initialOccupations)  ")
        println(io, "leadingConfigs:           $(dist.leadingConfigs)  ")
    end

    struct   PhotoAbsorption              <:  AbstractSimulationProperty   end
    struct   DecayPathes                  <:  AbstractSimulationProperty   end
    struct   ElectronCoincidence          <:  AbstractSimulationProperty   end
    struct   TimeBinnedPhotonIntensity    <:  AbstractSimulationProperty   end
    struct   MeanLineWidth                <:  AbstractSimulationProperty   end

    """
    abstract type  Cascade.AbstractSimulationMethod`
        ... defines a abstract and a list of singleton data types for the properties that can be 'simulated' from a given
            set of a set Cascade.AbstractData.

        + struct ProbPropagation     ... to propagate the (occupation) probabilites of the levels until no further changes occur.
        + struct MonteCarlo          ... to simulate the cascade decay by a Monte-Carlo approach of possible pathes (not yet considered).
        + struct RateEquations       ... to solve the cascade by a set of rate equations (not yet considered).
    """
    abstract type  AbstractSimulationMethod                  end
    struct   ProbPropagation  <:  AbstractSimulationMethod   end
    struct   MonteCarlo       <:  AbstractSimulationMethod   end
    struct   RateEquations    <:  AbstractSimulationMethod   end

    """
    `struct  Cascade.SimulationSettings`  ... defines settings for performing the simulation of some cascade (data).

        + printTree           ::Bool        ... Print the cascade tree in a short form
        + printLongTree       ::Bool        ... Print the cascade tree in a long form
        + initialPhotonEnergy ::Float64     ... Photon energy for which photoionization data are considered.
    """
    struct  SimulationSettings
        printTree             ::Bool
        printLongTree         ::Bool
        initialPhotonEnergy   ::Float64
    end

    """
    `Cascade.SimulationSettings()`  ... constructor for an 'empty' instance of a Cascade.Block.
    """
    function SimulationSettings()
        SimulationSettings(false, false, 0.)
    end

    # `Base.show(io::IO, settings::SimulationSettings)`  ... prepares a proper printout of the variable settings::SimulationSettings.
    function Base.show(io::IO, settings::SimulationSettings)
        println(io, "printTree:                $(settings.printTree)  ")
        println(io, "printLongTree:            $(settings.printLongTree)  ")
        println(io, "initialPhotonEnergy:      $(settings.initialPhotonEnergy)  ")
    end

    """
    `struct  Cascade.Simulation`  ... defines a simulation on some given cascade (data).

        + name            ::String                              ... Name of the simulation
        + property        ::Cascade.AbstractSimulationProperty
            ... Property that is to be considered in this simulation of the cascade (data).
        + method          ::Cascade.AbstractSimulationMethod
            ... Method that is used in the cascade simulation; cf. Cascade.SimulationMethod.
        + settings        ::Cascade.SimulationSettings          ... Settings for performing these simulations.
        + computationData ::Array{Dict{String,Any},1}           ... Date on which the simulations are performed
    """
    struct  Simulation
        name              ::String
        property          ::Cascade.AbstractSimulationProperty
        method            ::Cascade.AbstractSimulationMethod
        settings          ::Cascade.SimulationSettings
        computationData   ::Array{Dict{String,Any},1}
    end

    """
    `Cascade.Simulation()`  ... constructor for an 'default' instance of a Cascade.Simulation.
    """
    function Simulation()
        Simulation("Default cascade simulation", Cascade.PhotoAbsorption(), Cascade.ProbPropagation(),
                   Cascade.SimulationSettings(), Array{Dict{String,Any},1}[] )
    end

    """
    `Cascade.Simulation(sim::Cascade.Simulation;`

                name=..,               property=..,             method=..,              settings=..,     computationData=.. )

        ... constructor for re-defining the computation::Cascade.Simulation.
    """
    function Simulation(sim::Cascade.Simulation;
        name::Union{Nothing,String}=nothing,                                  property::Union{Nothing,Cascade.AbstractSimulationProperty}=nothing,
        method::Union{Nothing,Cascade.AbstractSimulationMethod}=nothing,      settings::Union{Nothing,Cascade.SimulationSettings}=nothing,
        computationData::Union{Nothing,Array{Dict{String,Any},1}}=nothing )

        if  name            == nothing   namex            = sim.name              else  namex            = name                end
        if  property        == nothing   propertyx        = sim.property          else  propertyx        = property            end
        if  method          == nothing   methodx          = sim.method            else  methodx          = method              end
        if  settings        == nothing   settingsx        = sim.settings          else  settingsx        = settings            end
        if  computationData == nothing   computationDatax = sim.computationData   else  computationDatax = computationData     end

    	Simulation(namex, propertyx, methodx, settingsx, computationDatax)
    end

    # `Base.show(io::IO, simulation::Cascade.Simulation)`  ... prepares a proper printout of the variable simulation::Cascade.Simulation.
    function Base.show(io::IO, simulation::Cascade.Simulation)
        println(io, "name:              $(simulation.name)  ")
        println(io, "property:          $(simulation.property)  ")
        println(io, "method:            $(simulation.method)  ")
        println(io, "computationData:   ... based on $(length(simulation.computationData)) cascade data sets ")
        println(io, "> settings:      \n$(simulation.settings)  ")
    end

    """
    `struct  Cascade.LineIndex`  ... defines a line index with regard to the various lineLists of data::Cascade.LineIndex.

        + lineSet      ::Cascade.AbstractData    ... refers to the data set for which this index is defined.
        + process      ::Basics.AbstractProcess  ... refers to the particular lineList of cascade (data).
        + index        ::Int64                   ... index of the corresponding line.
    """
    struct  LineIndex
        lineSet        ::Cascade.AbstractData
        process        ::Basics.AbstractProcess
        index          ::Int64
    end

    # `Base.show(io::IO, index::Cascade.LineIndex)`  ... prepares a proper printout of the variable index::Cascade.LineIndex.
    function Base.show(io::IO, index::Cascade.LineIndex)
        println(io, "lineSet (typeof):      $(typeof(index.lineSet))  ")
        println(io, "process:               $(index.process)  ")
        println(io, "index:                 $(index.index)  ")
    end

    """
    `mutable struct  Cascade.Level`  ... defines a level specification for dealing with cascade transitions.

        + energy       ::Float64                     ... energy of the level.
        + J            ::AngularJ64                  ... total angular momentum of the level
        + parity       ::Basics.Parity               ... total parity of the level
        + NoElectrons  ::Int64                       ... total number of electrons of the ion to which this level belongs.
        + relativeOcc  ::Float64                     ... relative occupation
        + parents      ::Array{Cascade.LineIndex,1}  ... list of parent lines that (may) populate the level.
        + daugthers    ::Array{Cascade.LineIndex,1}  ... list of daugther lines that (may) de-populate the level.
    """
    mutable struct  Level
        energy         ::Float64
        J              ::AngularJ64
        parity         ::Basics.Parity
        NoElectrons    ::Int64
        relativeOcc    ::Float64
        parents        ::Array{Cascade.LineIndex,1}
        daugthers      ::Array{Cascade.LineIndex,1}
    end

    # `Base.show(io::IO, level::Cascade.Level)`  ... prepares a proper printout of the variable level::Cascade.Level.
    function Base.show(io::IO, level::Cascade.Level)
        println(io, "energy:        $(level.energy)  ")
        println(io, "J:             $(level.J)  ")
        println(io, "parity:        $(level.parity)  ")
        println(io, "NoElectrons:   $(level.NoElectrons)  ")
        println(io, "relativeOcc:   $(level.relativeOcc)  ")
        println(io, "parents:       $(level.parents)  ")
        println(io, "daugthers:     $(level.daugthers)  ")
    end

    """
    `Base.:(==)(leva::Cascade.Level, levb::Cascade.Level)`  ... returns true if both levels are equal and false otherwise.
    """
    function  Base.:(==)(leva::Cascade.Level, levb::Cascade.Level)
        if  leva.energy == levb.energy   &&   leva.J == levb.J  && leva.parity == levb.parity  && leva.NoElectrons == levb.NoElectrons
                return( true )
        else    return( false )
        end
    end

    """
    `struct  Cascade.AbsorptionCrossSection`
        ... defines the absorption cross section for a particular (incident) photon energy in terms of its discrete and
            (direct photoionization) contributions. Of course, this absorption cross section depends on the relative population
            of the initial levels.

        + photonEnergy ::Float64              ... incident photon energy/photon-energy dependence of the absorption spectrum.
        + excitationCS ::Basics.EmProperty    ... contribution due to discrete excitation processes.
        + ionizationCS ::Basics.EmProperty    ... contribution due to contineous ionization processes.
    """
    struct  AbsorptionCrossSection
        photonEnergy   ::Float64
        excitationCS   ::Basics.EmProperty
        ionizationCS   ::Basics.EmProperty
    end

    # `Base.show(io::IO, cs::Cascade.AbsorptionCrossSection)`  ... prepares a proper printout of the variable cs::Cascade.AbsorptionCrossSection.
    function Base.show(io::IO, cs::Cascade.AbsorptionCrossSection)
        println(io, "photonEnergy:      $(cs.photonEnergy)  ")
        println(io, "excitationCS:      $(cs.excitationCS)  ")
        println(io, "ionizationCS:      $(cs.ionizationCS)  ")
    end

    include("module-Cascade-inc-computations.jl")
    include("module-Cascade-inc-dielectronic-capture.jl")
    include("module-Cascade-inc-hollow-ion.jl")
    include("module-Cascade-inc-photoabsorption.jl")
    include("module-Cascade-inc-photoexcitation.jl")
    include("module-Cascade-inc-photoionization.jl")
    include("module-Cascade-inc-stepwise-decay.jl")
    include("module-Cascade-inc-simulations.jl")

end # module


text:
@cache mutable struct DImplicitEulerCache{uType,duType,uNoUnitsType,N} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  du::duType
  atmp::uNoUnitsType
  nlsolver::N
end

mutable struct DImplicitEulerConstantCache{duType,N} <: OrdinaryDiffEqConstantCache
  du::duType
  nlsolver::N
end

function alg_cache(alg::DImplicitEuler,du,u,res_prototype,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  γ, c = 1, 1
  α = 1
  nlsolver = build_nlsolver(alg,u,uprev,p,t,dt,f,res_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,α,Val(false))

  DImplicitEulerConstantCache(du,nlsolver)
end

function alg_cache(alg::DImplicitEuler,du,u,res_prototype,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})

  γ, c = 1, 1
  α = 1
  nlsolver = build_nlsolver(alg,u,uprev,p,t,dt,f,res_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,α,Val(true))

  atmp = similar(u,uEltypeNoUnits)

  DImplicitEulerCache(u,uprev,uprev2,du,atmp,nlsolver)
end

@cache mutable struct DABDF2ConstantCache{duType,N,dtType,rate_prototype} <: OrdinaryDiffEqConstantCache
  du::duType
  nlsolver::N
  eulercache::DImplicitEulerConstantCache
  dtₙ₋₁::dtType
  fsalfirstprev::rate_prototype
end

function alg_cache(alg::DABDF2,du,u,res_prototype,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  γ, c = 1//1, 1
  α = 1//1
  nlsolver = build_nlsolver(alg,u,uprev,p,t,dt,f,res_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,α,Val(false))
  eulercache = DImplicitEulerConstantCache(du,nlsolver)

  dtₙ₋₁ = one(dt)
  fsalfirstprev = rate_prototype

  DABDF2ConstantCache(du,nlsolver, eulercache, dtₙ₋₁, fsalfirstprev)
end

@cache mutable struct DABDF2Cache{uType,duType,rateType,uNoUnitsType,N,dtType} <: OrdinaryDiffEqMutableCache
  uₙ::uType
  uₙ₋₁::uType
  uₙ₋₂::uType
  du::duType
  fsalfirst::rateType
  fsalfirstprev::rateType
  atmp::uNoUnitsType
  nlsolver::N
  eulercache::DImplicitEulerCache
  dtₙ₋₁::dtType
end

function alg_cache(alg::DABDF2,du,u,res_prototype,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  γ, c = 1//1, 1
  α = 1//1
  nlsolver = build_nlsolver(alg,u,uprev,p,t,dt,f,res_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,α,Val(true))
  fsalfirst = zero(rate_prototype)

  fsalfirstprev = zero(rate_prototype)
  atmp = similar(u,uEltypeNoUnits)

  eulercache = DImplicitEulerCache(u,uprev,uprev2,du,atmp,nlsolver)

  dtₙ₋₁ = one(dt)

  DABDF2Cache(u,uprev,uprev2,du,fsalfirst,fsalfirstprev,atmp,
              nlsolver,eulercache,dtₙ₋₁)
end

ext: jl


text:
using SecureSessions
using Base.Test
using HttpServer
using URIParser
using Requests

##################################################
### Test username and password permissibility
username1 = "John Smith"
username2 = "<script>do_something_terrible();</script>"
@test  username_is_permissible(username1)
@test !username_is_permissible(username2)

password1 = "v@l1d_pasSw0rd"
password2 = "bad_password"
@test  password_is_permissible(password1)
@test !password_is_permissible(password2)

##################################################
### Test secure cookies

# Override default globals
encrypted_sessions_only = false    # Avoids requiring https for this test

# Define app
function app(req::Request)
    res = Response()
    uri = URI(req.resource)
    if uri.path == "/set_secure_cookie"
	data = bytestring(req.data)
	create_secure_session_cookie(data, res, "sessionid")
    elseif uri.path == "/read_cookie"
        res.data = get_session_cookie_data(req, "sessionid")
    else
	res.status = 404
    end
    res
end

# Define and run server
server = Server((req, res) -> app(req))
@async run(server, port = 8000)
sleep(1.0)

# Post data "John Smith", which will be encrypted and included in the "sessionid" cookie.
username = "John Smith"
res1     = Requests.post("http://localhost:8000/set_secure_cookie"; data = username)

# Extract the cookie from the response and post it back to the server to be decrypted.
# Compare the returned data to the data originally encrypted.
cookie = res1.cookies["sessionid"]
res2   = Requests.post("http://localhost:8000/read_cookie"; cookies = [cookie])
@test bytestring(res2.data) == username

##################################################
### Test password hashing
sp = StoredPassword("pwd_alice")
@test  password_is_valid("pwd_alice", sp)
@test !password_is_valid("pwd_bob",   sp)

# EOF


text:
import Patchwork.div
import Base: size

export inset,
       empty,
       container,
       offset,
       width,
       height,
       minwidth,
       minheight,
       maxwidth,
       maxheight,
       size,
       topleft,
       midtop,
       topright,
       midleft,
       middle,
       midright,
       bottomleft,
       midbottom,
       bottomright,
       vertical,
       horizontal,
       depth,
       left,
       right,
       top,
       bottom,
       inward,
       outward,
       floating,
       flow,
       hbox,
       vbox,
       vskip,
       hskip,
       grow,
       shrink,
       flex,
       wrap,
       axisstart,
       axisend,
       center,
       spacebetween,
       spacearound,
       stretch,
       packitems,
       packacross,
       packlines,
       space,
       pad,
       padcontent,
       hidden,
       visible,
       scroll,
       auto,
       clip

# 0. Width and height

@api width => (Width <: Tile) begin
    doc("Set the width of a tile")
    typedarg(prefix::AbstractString="",
        doc=md"""either `""`, `"min"` or `"max"`. See `minwidth` and `maxwidth`.""")
    arg(width::Length, doc="The width")
    curry(tile::Tile, doc="the tile to set the height of")
end
render(t::Width, state) = begin
    prefixed = t.prefix == "" ? "width" : t.prefix * "Width"
    render(t.tile, state) & style(@d(prefixed => t.width))
end

@api height => (Height <: Tile) begin
    doc("Set the height of a tile")
    typedarg(prefix::AbstractString="",
        doc=md"""either `""`, `"min"` or `"max"`. See `minheight` and `maxheight`.""")
    arg(height::Length, doc="the height")
    curry(tile::Tile, doc="the tile to set the height of")
end

render(t::Height, state) =
    render(t.tile, state) &
        style(@d((t.prefix == "" ? "height" : t.prefix * "Height") => t.height))

minwidth(w, x...) = width("min", w, x...)
minheight(h, x...) = height("min", h, x...)

@apidoc minwidth => (Width <: Tile) begin
    doc("Set the minimum width of a tile")
    arg(width::Length, doc="the width")
    curry(tile::Tile, doc="the tile to set the width of")
end

@apidoc minheight => (Height <: Tile) begin
    doc("Set the minimum height of a tile")
    arg(height::Length, doc="the height")
    curry(tile::Tile, doc="the tile to set the height of")
end

maxwidth(w, x...) = width("max", w, x...)
maxheight(h, x...) = height("max", h, x...)

@apidoc maxwidth => (Width <: Tile) begin
    doc("Set the maximum width of a tile")
    arg(width::Length, doc="the width")
    curry(tile::Tile, doc="the tile to set the width of")
end

@apidoc maxheight => (Height <: Tile) begin
    doc("Set the maximum height of a tile")
    arg(height::Length, doc="the height")
    curry(tile::Tile, doc="the tile to set the height of")
end

size(w::Length, h::Length, t) =
    width(w, height(h, t))
size(w::Length, h::Length) =
    t -> size(w, h, t)

@apidoc size => (Height <: Tile) begin
    doc("Set the width and the height of a tile")
    arg(width::Length, doc="the width")
    arg(height::Length, doc="the height")
    curry(tile::Tile, doc="the tile to set the size of")
end
container(w, h) =
    empty |> size(w, h)

@apidoc container => (Height <: Tile) begin
    doc("Make an empty tile of a given size")
    arg(width::Length, doc="the width")
    arg(height::Length, doc="the height")
end

# 1. Placing a Tile inside another
abstract Position
abstract Corner <: Position

@terms Corner begin
    topleft     => TopLeft
    midtop      => MidTop
    topright    => TopRight
    midleft     => MidLeft
    middle      => Middle
    midright    => MidRight
    bottomleft  => BottomLeft
    midbottom   => MidBottom
    bottomright => BottomRight
end

name{C <: Corner}(::C) = string(C)

@api offset => (Relative{T <: Corner} <: Position) begin
    doc("Create an offset relative to a corner")
    arg(c::T, doc="The corner to offset from")
    arg(x::Length, doc="The horizontal offset")
    arg(y::Length, doc="The vertical offset")
    # z::Length
end

@api inset => (Inset <: Tile) begin
    doc("Position a tile inside a container tile at a corner or an offset")
    typedarg(position::Position=topleft, doc="The corner or the offset")
    arg(container::Tile, doc="The container")
    curry(contained::Tile, doc="The tile to be placed inside")
end

render_position(p::TopLeft, x, y) =
    @d(:top => y, :left => x)
render_position(p::MidTop, x, y) =
    @d(:left =>  50cent, :top => y,
     :transform => "translate(-50%)",
     :marginLeft => x)
render_position(p::TopRight, x, y) =
    @d(:top => x, :right => y)
render_position(p::MidLeft, x, y) =
    @d(:top => 50cent, :left => x,
     :marginTop => y,
     :transform => "translate(0, -50%)")
render_position(p::Middle, x, y) =
    @d(:top => 50cent, :left=>50cent,
     :marginLeft => x, :marginTop => y,
     :transform => "translate(-50%, -50%)")
render_position(p::MidRight, x, y) =
    @d(:top => 50cent,
    :transform => "translate(0, -50%)",
    :marginTop => y, :right => x)
render_position(p::BottomLeft, x, y) =
    @d(:bottom => y, :left => x)
render_position(p::MidBottom, x, y) =
    @d(:left => 50cent, :bottom => y,
     :marginLeft => x,
     :transform => "translate(-50%)")
render_position(p::BottomRight, x, y) =
    @d(:bottom => y, :right => x)

render_position(c::Corner) = style(render_position(c, 0, 0))
render_position{C <: Corner}(p::Relative{C}) =
    style(render_position(C(), p.x, p.y))

render(tile::Inset, state) = begin
    outer = render(tile.container, state)
    inner = render(tile.contained, state)

    outer &= style(@d(:position => :relative))
    inner &= style(@d(:position => :absolute))

    outer << (inner & render_position(tile.position))
end

# 2. Axes, Directions and Flow

abstract Axis

abstract FixedAxis <: Axis
@terms FixedAxis begin
    horizontal => Horizontal
    vertical => Vertical
    depth => Depth
end

abstract  Side{Perpendicular <: Axis}
@terms Side{Horizontal} begin
    left => Left
    right => Right
end

@terms Side{Vertical} begin
    top => TopSide
    bottom => Bottom
end

name(s::Left) = "Left"
name(s::Right) = "Right"
name(s::TopSide) = "Top"
name(s::Bottom) = "Bottom"

@terms Side{Depth} begin
    inward => Inward
    outward => Outward
end

# TODO: render inward, outward flow

abstract FlowRelativeAxis <: Axis
@terms FlowRelativeAxis begin
    mainaxis => MainAxis
    crossaxis => CrossAxis
end

@terms Side{MainAxis} begin
    mainstart => MainStart
    mainend => MainEnd
end

@terms Side{CrossAxis} begin
    crossstart => CrossStart
    crossend => CrossEnd
end

abstract FlexContainer <: Tile

render(f::FlexContainer, state) =
    addclasses(render(f.tile, state), classes(f))

@api flow => (Flow{A <: FixedAxis} <: FlexContainer) begin
    doc("Flow a list of tiles along the horizontal or the vertical axis")
    typedarg(axis::A, doc="The axis to layout along")
    curry(tiles::TileList, doc="The list of tiles to layout")
    kwarg(reverse::Bool=false, doc="Should the layout be in the reverse order?")
end
Flow{T}(x::T, y, z) = Flow{T}(x, y, z) # Julia issue 10641

flow(axis::FixedAxis, tiles...; reverse=false) =
    flow(axis, [t for t in tiles]; reverse=reverse)

Base.reverse(f::Flow) =
    flow(f.axis, f.tiles, !f.reverse)

classes(f::Flow{Horizontal}) =
    f.reverse ? "flow flow-reverse horizontal" : "flow horizontal"
classes(f::Flow{Vertical}) =
    f.reverse ? "flow flow-reverse vertical" : "flow vertical"

render(f::Flow, state) =
    addclasses(render(f.tiles, :div, state), classes(f))

hbox(arg::(@compat Union{AbstractArray, Tuple})) = flow(horizontal, arg)
vbox(arg::(@compat Union{AbstractArray, Tuple})) = flow(vertical, arg)

hbox(args...) = hbox(args)
vbox(args...) = vbox(args)

@apidoc hbox => (Flow <: Tile) begin
    doc(md"Arrange tiles horizontally. `hbox(args...)`
is equivalent to `flow(horizontal, args...)`")
    arg(tiles::TileList)
end

@apidoc vbox => (Flow <: Tile) begin
    doc("Arrange tiles vertically. `vbox(args...)`
is equivalent to `flow(vertical, args...)`")
    arg(tiles::TileList)
end

vskip(y) = size(0px, y, empty)
hskip(x) = size(x, 0px, empty)

@api ordering => (FlowOrder <: FlexContainer) begin
    arg(ordering::AbstractArray)
    curry(flow::FlexContainer)
end

# TODO: render ordering

@api wrap => (Wrap <: FlexContainer) begin
    doc(md"Wrap a `flow` of tiles")
    arg(tile::FlexContainer, doc=md"Either a `vbox` or a `hbox`")
    kwarg(reverse::Bool=false, doc="Should the wrapping be reversed in direction?")
end
wrapreverse(t) = wrap(t, reverse=true)

classes(f::Wrap) =
    f.reverse ? "flex-wrap-reverse" : "flex-wrap"

# 3. Flexing and alignment

@api floating => (FloatingTile{T <: Side{Horizontal}} <: Tile) begin
    typedarg(side::T)
    curry(tile::Tile)
end

render(f::FloatingTile, state) =
    render(f.tile, state) & style(@d(:float => lowercase(name(f.side))))

@api grow => (Grow <: Tile) begin
    doc(md"""Expand a tile along the main axis to fit extra space in the parent
    `hbox` or `vbox`.""")

    arg(
        factor::Float64,
        doc="""The relative rate at which this tile will expand compared to other
        tiles that can grow.""",
    )

    curry(
        tile::Tile,
        doc=md"The tile to stretch. This tile should go inside a `vbox` or `hbox`",
    )
end
grow(t::Tile) = grow(1.0, t)
grow(t::AbstractVector) = map(grow, t)

render(t::Grow, state) =
    render(t.tile, state) & style(@d(:flexGrow => t.factor))

@api shrink => (Shrink <: Tile) begin
    doc(md"Shrink a tile along the main axis to accomodate space in the parent `hbox` or `vbox`.")
    arg(factor::Float64,
        doc="The relative rate at which this tile will shrink compared to other tiles that can shrink.")
    curry(tile::Tile,
        doc=md"The tile to stretch. This tile should go inside a `vbox` or `hbox`")
end
shrink(t::Tile) = shrink(1.0, t)
shrink(t::AbstractVector) = map(shrink, t)

render(t::Shrink, state) =
    render(t.tile, state) & style(@d(:flexShrink => t.factor))

@api flexbasis => (FlexBasis <: Tile) begin
    arg(basis::(@compat Union{Length, Symbol}))
    curry(tile::Tile)
end

render(t::FlexBasis, state) =
    render(t.tile, state) & style(@d(:flexBasis => t.basis))

# Flex ignores the width and distributes forcefully
flex(factor::Real, t) =
    flexbasis(0mm, grow(factor, t))
flex(t) = flex(1.0, t)
flex() = flex(empty)
flex{T <: Real}(factor::AbstractVector{T}, t) =
    map(flex, factor, t)
flex(t::AbstractVector) =
    map(flex, t)

flex(factor::Real) = t -> flex(factor, t)
flex{T <: Real}(factor::AbstractVector{T}) = t -> flex(factor, t)

@apidoc flex => (FlexBasis <: Tile) begin
    doc("Ignore the width (in a hbox) or height (in a vbox) and stretch or shrink a tile to fill / distribute remaining space")
    arg(tile::Tile, doc="The tile to flex")
end

# Flow alignment
abstract Packing

@terms Packing begin
    axisstart => AxisStart
    axisend => AxisEnd
    center => AxisCenter
    stretch => Stretch
    baseline => Baseline
    spacebetween => SpaceBetween
    spacearound => SpaceAround
end

@api packlines => (PackedLines{T <: Packing} <: FlexContainer) begin
    doc("Pack wrapped lines of tiles across the cross axis")
    typedarg(packing::T, doc="The kind of packing to use")
    curry(tile::FlexContainer, doc="The flex container")
end

classes(t::PackedLines{AxisStart}) = "pack-lines-start"
classes(t::PackedLines{AxisEnd}) = "pack-lines-end"
classes(t::PackedLines{AxisCenter}) = "pack-lines-center"
classes(t::PackedLines{Stretch}) = "pack-lines-stretch"
classes(t::PackedLines{SpaceBetween}) = "pack-lines-space-between"
classes(t::PackedLines{SpaceAround}) = "pack-lines-space-around"

@api packitems => (PackedItems{T <: Packing} <: FlexContainer) begin
    doc("Pack items in a flex container along the main axis")
    typedarg(packing::T, doc="The kind of packing to use")
    curry(tile::FlexContainer, doc="The flex container")
end

classes(t::PackedItems{AxisStart}) = "pack-start"
classes(t::PackedItems{AxisEnd}) = "pack-end"
classes(t::PackedItems{AxisCenter}) = "pack-center"
classes(t::PackedItems{SpaceBetween}) = "pack-space-between"
classes(t::PackedItems{SpaceAround}) = "pack-space-around"

@api packacross => (PackedAcross{T <: Packing} <: FlexContainer) begin
    doc("Stretch or provide spacing around items in the cross axis")
    typedarg(packing::T, doc="The kind of packing to use")
    curry(tile::FlexContainer, doc="The flex container")
end

classes(t::PackedAcross{AxisStart}) = "pack-across-start"
classes(t::PackedAcross{AxisEnd}) = "pack-across-end"
classes(t::PackedAcross{AxisCenter}) = "pack-across-center"
classes(t::PackedAcross{Stretch}) = "pack-across-stretch"
classes(t::PackedAcross{Baseline}) = "pack-across-baseline"

# 4. Padding

immutable Container <: Tile
    tile::Tile
end

render(cont::Container, state) = Elem(:div, render(cont.tile, state))

@api padcontent => (PadContent <: Tile) begin
    typedarg(sides::AbstractVector=allsides)
    arg(length::Length)
    curry(tile::Tile)
end

render(t::PadContent, state) =
    render(t.tile, state) &
        style(mapparts(allsides, t.sides, "padding", "", t.length))

pad(len::Length, tile) =
    padcontent(len, Container(tile))

pad(len::Length) =
    t -> pad(len, t)

pad(sides::AbstractVector, len::Length, tile) =
    padcontent(sides, len, Container(tile))

pad(sides::AbstractVector, len::Length) =
    tile -> padcontent(sides, len, Container(tile))

@apidoc pad => (PadContent <: Tile) begin
    doc("Wrap a tile in a container with the specified padding")
    typedarg(sides::AbstractVector=allsides, doc="Sides to pad")
    arg(length::Length, doc="Amount of padding")
    curry(tile::Tile, doc="The tile to pad")
end

# Clipping

abstract Overflow

@terms Overflow begin
    hidden => Hidden
    visible => Visible
    scroll => Scroll
    auto => AutoClip
end

@api clip => (Clip <: Tile) begin
    doc(md"Clip a tile to the dimensions set using `width` and `height`.")
    typedarg(
        overflow::Overflow=auto,
        doc=md"""The method for clipping. Valid values are `hidden`, `visible`,
                 `scroll` and `auto`."""
        )
    curry(tile::Tile, doc="The tile to clip.")
end

name(::Hidden) = "hidden"
name(::Visible) = "visible"
name(::Scroll) = "scroll"
name(::AutoClip) = "auto"

render(t::Clip, state) =
    render(Container(t.tile), state) &
        @d(:style => @d(:overflow => name(t.overflow)), :className => "scrollbar")


text:
export CoefficientHomotopy

"""
    CoefficientHomotopy(F::Union{AbstractSystem,System}; start_coefficients, target_coefficients)
    CoefficientHomotopy(F::Union{AbstractSystem,System}, start_coefficients, target_coefficients)

Construct the homotopy ``H(x, t) = ∑_{a ∈ Aᵢ} (c_a t + (1-t)d_a) x^a`` where ``c_a`` are
the start coefficients and ``d_a`` the target coefficients.
"""
struct CoefficientHomotopy{T<:AbstractSystem} <: AbstractHomotopy
    F::T
    start_coeffs::Vector{ComplexF64}
    target_coeffs::Vector{ComplexF64}
    #cache
    x::Vector{ComplexF64}
    t_cache::Base.RefValue{ComplexF64}
    t_taylor_cache::Base.RefValue{ComplexF64}
    coeffs::Vector{ComplexF64}
    dt_coeffs::Vector{ComplexF64}
    taylor_coeffs::TaylorVector{2,ComplexF64}
end

function CoefficientHomotopy(
    F;
    start_coefficients::AbstractVector,
    target_coefficients::AbstractVector,
)
    CoefficientHomotopy(F, start_coefficients, target_coefficients)
end
function CoefficientHomotopy(
    F::ModelKit.System,
    start_coeffs::AbstractVector,
    target_coeffs::AbstractVector,
)
    @assert length(start_coeffs) == length(target_coeffs) == length(F.parameters)
    CoefficientHomotopy(ModelKitSystem(F), start_coeffs, target_coeffs)
end
function CoefficientHomotopy(
    F::AbstractSystem,
    start_coeffs::AbstractVector,
    target_coeffs::AbstractVector,
)
    CoefficientHomotopy(
        F,
        Vector{ComplexF64}(start_coeffs),
        Vector{ComplexF64}(target_coeffs),
    )
end

function CoefficientHomotopy(
    F::AbstractSystem,
    start_coeffs::Vector{ComplexF64},
    target_coeffs::Vector{ComplexF64},
)
    m = length(target_coeffs)
    @assert length(start_coeffs) == m
    CoefficientHomotopy(
        F,
        start_coeffs,
        target_coeffs,
        zeros(ComplexF64, size(F, 2)),
        Ref(complex(NaN)),
        Ref(complex(NaN)),
        zeros(ComplexF64, m),
        start_coeffs - target_coeffs,
        TaylorVector{2}(ComplexF64, length(target_coeffs)),
    )
end

Base.size(H::CoefficientHomotopy) = size(H.F)

start_parameters!(H::CoefficientHomotopy, p::AbstractVector) = (H.start_coeffs .= p; H)
target_parameters!(H::CoefficientHomotopy, p::AbstractVector) = (H.target_coeffs .= p; H)

function coeffs!(H::CoefficientHomotopy, t)
    t == H.t_cache[] && return H.coeffs

    if isreal(t)
        s = real(t)
        s1 = 1.0 - s
        @inbounds for i = 1:length(H.start_coeffs)
            H.coeffs[i] = s * H.start_coeffs[i] + s1 * H.target_coeffs[i]
        end
    else
        t1 = 1.0 - t
        @inbounds for i = 1:length(H.start_coeffs)
            H.coeffs[i] = t * H.start_coeffs[i] + t1 * H.target_coeffs[i]
        end
    end
    H.t_cache[] = t

    H.coeffs
end

function evaluate!(u, H::CoefficientHomotopy, x, t)
    evaluate!(u, H.F, x, coeffs!(H, t))
end

function evaluate_and_jacobian!(u, U, H::CoefficientHomotopy, x, t)
    evaluate_and_jacobian!(u, U, H.F, x, coeffs!(H, t))
end

function taylor!(u, v::Val{1}, H::CoefficientHomotopy, x, t)
    evaluate!(u, H.F, x, H.dt_coeffs)
end

function taylor_coeffs!(H::CoefficientHomotopy, t::Union{ComplexF64,Float64})
    t == H.t_taylor_cache[] && return H.taylor_coeffs

    coeffs!(H, t)
    for i = 1:length(H.taylor_coeffs)
        H.taylor_coeffs[i, 1] = H.coeffs[i]
        H.taylor_coeffs[i, 2] = H.dt_coeffs[i]
    end
    H.t_taylor_cache[] = t

    H.taylor_coeffs
end

function taylor!(u, v::Val, H::CoefficientHomotopy, tx::TaylorVector, t)
    taylor!(u, v, H.F, tx, taylor_coeffs!(H, t))
end


text:
# -*- coding: utf-8 -*-
# ## Grafting models together
# Much of scientific thought and scholarship involves reading papers and textbooks with different models and synthesizing a novel model from components found in existing, well studied models.
#
# We refer to this as model grafting, where a component is taking from one model and grafted onto another model. SemanticModels supports automating much of the low level code details from this task to free scientists to think about what features to combine instead of the mechanical aspects of changing the code by hand.
#
# This notebook is an example based on the SEIR model and the ScalingModel examples in the epirecipes cookbook.

using DifferentialEquations
using SemanticModels.Parsers
using SemanticModels.ModelTools
using SemanticModels.ModelTools.ExpODEModels

# ## Loading the original model
# We use parsefile to load the model into an expression. The original model is an SEIR model which has 4 states suceptible, exposed, infected, and recovered. It has parameters $\beta, \gamma, \mu, \sigma$.

expr1 = parsefile("../examples/epicookbook/src/SEIRmodel.jl")
model1 = model(ExpODEModel, expr1)

module1 = eval(model1.expr)

# ## Running our baseline model
#
# The code that defines the baseline model creates a module for that model to run in. This ensures that the code will not have unintended sideeffects when run in a julia process with other models. The entrypoint to this module is called `main`, which is a function that has no arguments and does the setup and execution of the model.

module1.main()

# seir_ode is the name of the function we want to modify
# an ODEProblem is defined by the right hand side of the equation.
# $du/dt = f(u, t)$

# The ScalingModel provides a population growth component that we want to graft onto the SEIR model to create an SEIR model with population dynamics. We load that model from its source file. You can inspect this file to see the definition of $dS/dt = r * (1 - S / K) * S - \beta * S * I$ which includes a population growth rate parameter $r$.

expr2 = parsefile("../examples/epicookbook/src/ScalingModel.jl")

# Once the ASTs are processed into a structured representation we can manipulate with regular julia code, we are able to write manipulations of the models that operate on a higher level than textual changes to the code.

model2 = model(ExpODEModel, expr2)
fluxes(x::ExpODEModel) = x.variables[1].flux

# Find the expression we want to graft
# vital dynamics S rate expression

fluxvar = fluxes(model2)[1].args[2].args[1]
popgrowth = replacevar(findassign(model2.funcs[1], fluxvar)[1], :K, :N).args[2].args[2]
ex = model1.variables[1].flux[1]
ex.args[2] = :($(popgrowth)+$(ex.args[2]))

# define N as the sum of the entries of Y ie. S+E+I

@assert argslist(:(function foo(x, y); return x+y; end)) == [:foo, :x, :y]

pushfirst!(bodyblock(model1.funcs[1]), :(N = sum(Y)))

# we need to add a new paramter to the function we are going to define
# this signature doesn't match the old signature so we are going to
# do some surgery on the main function to add that parameter
# this parameter instead could be added to the vector of parameters.

pusharg!(model1.funcs[1], :r)
# gensym gives us a unique name for the new function
g_func = gensym(argslist(model1.funcs[1])[1])
argslist(model1.funcs[1])[1] = g_func

# ## Model Augmentations often require new parameters
#
# When we add the population growth term to the SEIR model, we introduce a new parameter $r$
# that needs to be supplied to the model. One problem with approaches that require scientists
# to modify source code is the fact that adding the new features necessitates changes to the
# APIs provided by the original author. SemanticModels.ModelTools provides a higher level API
# for making these changes that assist in propagating the necessary changes to the API.
#
# For example, in this code we need to add an argument to the entrypoint function `main` and
# provide an anonymous function that conforms to the API that `DifferentialEquations` expects
# from its inputs.

mainx = findfunc(model1.expr, :main)[end]
pusharg!(mainx, :λ)

# An `ODEProblem` expects the user to provide a function $f(du, u, p, t)$ which takes the current fluxes, current system state, parameters, and current time as its arguments and updates the value of `du`. Since our new function `g_func` does not satisfy this interface, we need to introduce a wrapper function that does.
#
# Here is an instance where having a smart compiler helps julia. In many dynamic languages where this kind of metaprogramming would be easy, the runtime is not smart enough to inline these anonymous functions, which means that there is additional runtime performance overhead to metaporgramming like this. Julia's compiler (and LLVM) can inline these functions which drastically reduces that overhead.

setarg!(model1.calls[end], :seir_ode, :((du,u,p,t)->$g_func(du,u,p,t,λ)))
@show model1.expr
NewModule = eval(model1.expr)

# ## Modeling configuration
#
# The following code sets up our modeling configuration with initial conditions and parameters. It represents the entry point to solving the model.

findfunc(model1.expr, :main)[end]

# ## Solving the new model
#
# Once we have changed the function `seir_ode` and adapted the API of `main` to suit we can do a parameter sweep over our new parameter by solving the problem with different values of $\lambda$.

newsol = NewModule.main(1)

# ## Parameter Estimation
# Adding a capability to a model usually introduces additional parameters that must be chosen. Analyzing a model requires developing a procedure for estimating those parameters or characterizing the effect of that parameter on the behavior of the modeled system.

# Here we  sweep over population growth rates to show what happens when the population growth rate changes.

scalegrowth(λ=1.0) = NewModule.main(λ)

println("S\tI\tR")
for λ in [1.0,1.1,1.2,1.3,1.4,1.5]
    S,I,R = scalegrowth(λ)(365)
    println("$S\t$I\t$R")
end

# ## It Works!
#
# This simulation allows an epidemiologist to examine the effects of population growth on an SEIR disease outbreak. A brief analysis of this simulation shows that as you increase the population growth rate, you increase the final population of infected people. More sophisticated analysis could be employed to show something more interesting about this model.
#
# We have shown how you can use SemanticModels.jl to combine features of various ODE systems and solve them with a state of the art solver to increase the capabilities of a code that implements a scientific model. We call this combination process grafting and believe that it supports a frequent use case of scientific programming.


text:
test("Add path to offset") do
    path = Vector{IntPoint}()

    push!(path, IntPoint(0, 0))
    push!(path, IntPoint(0, 1))

    c = ClipperOffset()

    add_path!(c, path, JoinTypeRound, EndTypeClosedPolygon)
end

test("Add paths to offset") do
    path1 = Vector{IntPoint}()

    push!(path1, IntPoint(0, 0))
    push!(path1, IntPoint(0, 1))

    path2 = Vector{IntPoint}()

    push!(path2, IntPoint(5, 0))
    push!(path2, IntPoint(2, 1))

    paths = Vector{IntPoint}[path1, path2]

    c = ClipperOffset()

    add_paths!(c, paths, JoinTypeRound, EndTypeClosedPolygon)
end

test("Clear") do
    path = Vector{IntPoint}()

    push!(path, IntPoint(0, 0))
    push!(path, IntPoint(0, 1))

    c = ClipperOffset()

    add_path!(c, path, JoinTypeRound, EndTypeClosedPolygon)

    Clipper.clear!(c)
end

test("Offset") do
  path = Vector{IntPoint}()

  push!(path, IntPoint(0, 0))
  push!(path, IntPoint(0, 1))

  c = ClipperOffset()

  add_path!(c, path, JoinTypeRound, EndTypeOpenRound)

  poly = execute(c, 1.0)

  @test string(poly) == "Array{Clipper.IntPoint,1}[Clipper.IntPoint[[1,2],[-1,2],[-1,-1],[1,-1]]]"
end


text:
###############################################################################
#
#   nmod.jl : Nemo nmod (integers modulo small n)
#
###############################################################################

###############################################################################
#
#   Type and parent object methods
#
###############################################################################

parent_type(::Type{nmod}) = NmodRing

elem_type(::Type{NmodRing}) = nmod

base_ring(a::NmodRing) = Union{}

base_ring(a::nmod) = Union{}

parent(a::nmod) = a.parent

function check_parent(a::nmod, b::nmod)
   a.parent != b.parent && error("Operations on distinct residue rings not supported")
end

###############################################################################
#
#   Basic manipulation
#
###############################################################################

function Base.hash(a::nmod, h::UInt)
   b = 0x1812aa3492cbf4d9%UInt
   return xor(xor(hash(a.data), h), b)
end

function zero(R::NmodRing)
   return nmod(UInt(0), R)
end

function one(R::NmodRing)
   if R.n == 1
      return nmod(UInt(0), R)
   else
      return nmod(UInt(1), R)
   end
end

iszero(a::nmod) = a.data == 0

isone(a::nmod) = a.parent.n == 1 ? a.data == 0 : a.data == 1

isunit(a::nmod) = a.parent.n == 1 ? a.data == 0 : gcd(a.data, a.parent.n) == 1

modulus(R::NmodRing) = R.n

function deepcopy_internal(a::nmod, dict::ObjectIdDict)
   R = parent(a)
   return nmod(deepcopy(a.data), R)
end

###############################################################################
#
#   Canonicalisation
#
###############################################################################

function canonical_unit(x::nmod)
  #the simple return x does not work
  # - if x == 0, this is not a unit
  # - if R is not a field....
  if iszero(x)
    return parent(x)(0)
  end
  g = gcd(modulus(x), data(x))
  u = divexact(data(x), g)
  a, b = ppio(modulus(x), u)
  if isone(a)
    r = u
  elseif isone(b)
    r = b
  else
    r = crt(fmpz(1), fmpz(a), fmpz(u), fmpz(b))
  end
  return parent(x)(r)
end

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, R::NmodRing)
   print(io, "Integers modulo ", signed(widen(R.n)))
end

function show(io::IO, a::nmod)
   print(io, signed(widen(a.data)))
end

needs_parentheses(x::nmod) = false

isnegative(x::nmod) = false

show_minus_one(::Type{nmod}) = true

###############################################################################
#
#   Unary operations
#
###############################################################################

function -(x::nmod)
   if x.data == 0
      return deepcopy(x)
   else
      R = parent(x)
      return nmod(R.n - x.data, R)
   end
end

###############################################################################
#
#   Binary operations
#
###############################################################################

function +(x::nmod, y::nmod)
   check_parent(x, y)
   R = parent(x)
   n = modulus(R)
   d = x.data + y.data - n
   if d > x.data
      return nmod(d + n, R)
   else
      return nmod(d, R)
   end
end

function -(x::nmod, y::nmod)
   check_parent(x, y)
   R = parent(x)
   n = modulus(R)
   d = x.data - y.data
   if d > x.data
      return nmod(d + n, R)
   else
      return nmod(d, R)
   end
end

function *(x::nmod, y::nmod)
   check_parent(x, y)
   R = parent(x)
   d = ccall((:n_mulmod2_preinv, :libflint), UInt, (UInt, UInt, UInt, UInt),
             x.data, y.data, R.n, R.ninv)
   return nmod(d, R)
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function *(x::Integer, y::nmod)
   R = parent(y)
   return R(widen(x)*signed(widen(y.data)))
end

*(x::nmod, y::Integer) = y*x

function *(x::Int, y::nmod)
   R = parent(y)
   if x < 0
      d = ccall((:n_mulmod2_preinv, :libflint), UInt, (UInt, UInt, UInt, UInt),
             UInt(-x), y.data, R.n, R.ninv)
      return -nmod(d, R)
   else
      d = ccall((:n_mulmod2_preinv, :libflint), UInt, (UInt, UInt, UInt, UInt),
             UInt(x), y.data, R.n, R.ninv)
      return nmod(d, R)
   end
end

*(x::nmod, y::Int) = y*x

function *(x::UInt, y::nmod)
   R = parent(y)
   d = ccall((:n_mulmod2_preinv, :libflint), UInt, (UInt, UInt, UInt, UInt),
             UInt(x), y.data, R.n, R.ninv)
   return nmod(d, R)
end

*(x::nmod, y::UInt) = y*x

+(x::nmod, y::Integer) = x + parent(x)(y)

+(x::Integer, y::nmod) = y + x

-(x::nmod, y::Integer) = x - parent(x)(y)

-(x::Integer, y::nmod) = parent(y)(x) - y

###############################################################################
#
#   Powering
#
###############################################################################

function ^(x::nmod, y::Int)
   R = parent(x)
   if y < 0
      x = inv(x)
      y = -y
   end
   d = ccall((:n_powmod2_preinv, :libflint), UInt, (UInt, Int, UInt, UInt),
             UInt(x.data), y, R.n, R.ninv)
   return nmod(d, R)
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(x::nmod, y::nmod)
   check_parent(x, y)
   return x.data == y.data
end

###############################################################################
#
#   Ad hoc comparison
#
###############################################################################

==(x::nmod, y::Integer) = x == parent(x)(y)

==(x::Integer, y::nmod) = parent(y)(x) == y

==(x::nmod, y::fmpz) = x == parent(x)(y)

==(x::fmpz, y::nmod) = parent(y)(x) == y

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(x::nmod)
   R = parent(x)
   (x == 0 && R.n != 1) && throw(DivideError())
   if R.n == 1
      return deepcopy(x)
   end
   s = [UInt(0)]
   g = ccall((:n_gcdinv, :libflint), UInt, (Ref{UInt}, UInt, UInt),
         pointer(s), x.data, R.n)
   g != 1 && error("Impossible inverse in ", R)
   return nmod(s[], R)
end

###############################################################################
#
#   Coprime bases
#
###############################################################################

# Bernstein, "Factoring into coprimes in essentially linear time"
# ppio(a,b) = (c,n) where v_p(c) = v_p(a) if v_p(b) != 0, 0 otherwise
# c*n = a or c = gcd(a, b^infty), n = div(a, c).
# This is used in various Euclidean domains for Chinese remaindering.

doc"""
   ppio(a::T, b::T)

> Split $a$ into $c*d$ where $c = gcd(a, b^\infty)$.
"""
function ppio(a::E, b::E) where E <: Integer
   c = gcd(a, b)
   n = div(a, c)
   g = gcd(c, n)
   while !isone(g)
      c *= g
      n = div(n, g)
      g = gcd(c, n)
   end
   return c, n
end

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::nmod, y::nmod)
   check_parent(x, y)
   fl, q = divides(x, y)
   if !fl
     error("Impossible inverse in ", R)
   end
   return q
end

function divides(a::nmod, b::nmod)
   check_parent(a, b)
   if iszero(a)
      return true, a
   end
   A = data(a)
   B = data(b)
   R = parent(a)
   m = modulus(R)
   gb = gcd(B, m)
   q, r = divrem(A, gb)
   if r != 0
      return false, b
   end
   ub = div(B, gb)
   _, x = ppio(m, ub)
   r = R(q)*inv(R(ub))
   return true, r
end

###############################################################################
#
#   GCD
#
###############################################################################

function gcd(x::nmod, y::nmod)
   check_parent(x, y)
   R = parent(x)
   d = gcd(gcd(x.data, R.n), y.data)
   if d == R.n
      return nmod(0, R)
   else
      return nmod(d, R)
   end
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(z::nmod)
   R = parent(z)
   return nmod(UInt(0), R)
end

function mul!(z::nmod, x::nmod, y::nmod)
   return x*y
end

function addeq!(z::nmod, x::nmod)
   return z + x
end

function add!(z::nmod, x::nmod, y::nmod)
   return x + y
end

###############################################################################
#
#   Random functions
#
###############################################################################

function rand(R::NmodRing)
   n = rand(UInt(0):R.n - 1)
   return nmod(n, R)
end

function rand(R::NmodRing, b::UnitRange{Int64})
   n = rand(b)
   return R(n)
end

###############################################################################
#
#   Promotions
#
###############################################################################

promote_rule(::Type{nmod}, ::Type{T}) where T <: Integer = nmod

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (R::NmodRing)()
   return nmod(UInt(0), R)
end

function (R::NmodRing)(a::Integer)
   n = R.n
   d = a%signed(widen(n))
   if d < 0
      d += n
   end
   return nmod(UInt(d), R)
end

function (R::NmodRing)(a::Int)
   n = R.n
   ninv = R.ninv
   if reinterpret(Int, n) > 0 && a < 0
      a %= Int(n)
   end
   d = reinterpret(UInt, a)
   if a < 0
      d += n
   end
   if d >= n
      d = ccall((:n_mod2_preinv, :libflint), UInt, (UInt, UInt, UInt),
             d, n, ninv)
   end
   return nmod(d, R)
end

function (R::NmodRing)(a::UInt)
   n = R.n
   ninv = R.ninv
   a = ccall((:n_mod2_preinv, :libflint), UInt, (UInt, UInt, UInt),
             a, n, ninv)
   return nmod(a, R)
end

function (R::NmodRing)(a::fmpz)
   d = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ref{fmpz}, UInt),
             a, R.n)
   return nmod(d, R)
end

function (R::NmodRing)(a::nmod)
   return a
end

###############################################################################
#
#   nmod constructor
#
###############################################################################

function ResidueRing(R::FlintIntegerRing, n::Int; cached::Bool=true)
   n <= 0 && throw(DomainError())
   return NmodRing(UInt(n), cached)
end

function ResidueRing(R::FlintIntegerRing, n::UInt; cached::Bool=true)
   return NmodRing(n, cached)
end