text:
"""
`module  JAC.StrongField
    ... a submodel of JAC that contains all methods to set-up and perform strong-field computations.
"""
module StrongField

    using    GSL, HypergeometricFunctions, Plots, Printf, SpecialFunctions, ..AngularMomentum, ..Basics, ..Continuum, ..Defaults, ..InteractionStrength, ..Radial, ..ManyElectron,
             ..Nuclear, ..Pulse, ..TableStrings

    export   aaaa

    """
    `abstract type StrongField.AbstractSFAObservable`
        ... defines an abstract and a number of types for the observables that can be computed with SFA amplitudes.

        + struct SfaNoObservable            ... an empty instance of an observable that does not help compute anything.
        + struct SfaEnergyDistribution      ... to compute the energy distribution of the photoelectrons.
        + struct SfaMomentumDistribution    ... to compute the momentum distribution of the photoelectrons.
        + struct SfaAzimuthalAngularDistribution     ... to compute the angular distribution (azimuthal angle phi) of photoelectrons with fixed energy and fixed polar angle.
        + struct SfaPolarAngularDistribution     ... to compute the angular distribution (polar angle theta) of photoelectrons with fixed energy and fixed azimuthal angle.
    """
    abstract type  AbstractSFAObservable                                 end
    struct   SfaNoObservable     <: StrongField.AbstractSFAObservable    end

    """
    `struct  StrongField.SfaEnergyDistribution  <: StrongField.AbstractSFAObservable`
        ... to compute in SFA the energy distribution of photoelectrons at given energies and angles.

        + theta         ::Float64             ... polar angle of the energy distribution
        + phi           ::Float64             ... azimuthal angle of the energy distribution
        + energies      ::Array{Float64,1}    ... specifies the photoelectron energy in the current units.
    """
    struct   SfaEnergyDistribution  <: StrongField.AbstractSFAObservable
        theta           ::Float64
        phi             ::Float64
        energies        ::Array{Float64,1}
    end

    """
    `StrongField.SfaEnergyDistribution(theta::Float64, phi::Float64, NoEnergies::Int64, maxEnergy::Float64)`
        ... defines an energy distribution for given (theta,phi) and for NoEnergies between 0 < energy <= maxEnergy;
            a dist::SfaEnergyDistribution is returned.
    """
    function SfaEnergyDistribution(theta::Float64, phi::Float64, NoEnergies::Int64, maxEnergy::Float64)
        energies = Float64[];      for  i=1:NoEnergies   push!(energies, i*maxEnergy / NoEnergies)   end
        SfaEnergyDistribution(theta, phi, energies)
    end

    function Base.string(obs::SfaEnergyDistribution)
        sa = "Compute photo-electron energy distribution at angles (theta,phi) = ($(obs.theta), $(obs.phi)) for the energies [a.u.]:"
        return( sa )
    end

    function Base.show(io::IO, obs::SfaEnergyDistribution)
        sa = string(obs);       print(io, sa, "\n");     print(io, "   ", obs.energies)
    end

    """
    `struct  StrongField.SfaMomentumDistribution  <: StrongField.AbstractSFAObservable`
        ... to compute in SFA the momentum distribution of photoelectrons for given polar angle theta and at given energies and azimuthal angles.

        + theta         ::Float64             ... polar angle of the energy distribution
        + phi           ::Array{Float64,1}    ... specifies the azimuthal angles
        + energies      ::Array{Float64,1}    ... specifies the photoelectron energy in the current units.
    """
    struct   SfaMomentumDistribution  <: StrongField.AbstractSFAObservable
        theta           ::Float64
        phis            ::Array{Float64,1}
        energies        ::Array{Float64,1}
    end

    """
    `StrongField.SfaMomentumDistribution(theta::Float64, NoPhi::Int64, NoEnergies::Int64, maxEnergy::Float64)`
        ... defines a momentum distribution for given theta and for NoEnergies between 0 < energy <= maxEnergy and NoPhi between 0 <= phi < 2pi;
            a dist::SfaMomentumDistribution is returned.
    """
    function SfaMomentumDistribution(theta::Float64, NoPhi::Int64, NoEnergies::Int64, maxEnergy::Float64)
        energies = Float64[];      for  i=1:NoEnergies   push!(energies, i*maxEnergy / NoEnergies)   end
        phis = Float64[];          for  i=1:NoPhi        push!(phis, i*2pi / NoPhi)                  end
        SfaMomentumDistribution(theta, phis, energies)
    end

    function Base.string(obs::SfaMomentumDistribution)
        sa = "Compute photo-electron momentum distribution at polar angle theta = $(obs.theta) for (energy [a.u.], phi):"
        return( sa )
    end

    function Base.show(io::IO, obs::SfaMomentumDistribution)
        sa = string(obs);       print(io, sa, "\n");     print(io, "   ", obs.energies);    print(io, "   ", obs.phi)
    end

    """
    `struct  StrongField.SfaAzimuthalAngularDistribution  <: StrongField.AbstractSFAObservable`
        ... to compute in SFA the angular distribution of photoelectrons for given polar angle theta and at given azimuthal angles for a fixed energy

        + theta         ::Float64             ... polar angle of the energy distribution
        + phi           ::Array{Float64,1}    ... specifies the azimuthal angles
        + energy        ::Float64             ... specifies the photoelectron energy in the current units.
    """
    struct   SfaAzimuthalAngularDistribution  <: StrongField.AbstractSFAObservable
        theta           ::Float64
        phis            ::Array{Float64,1}
        energy          ::Float64
    end

    """
    `StrongField.SfaAzimuthalAngularDistribution(theta::Float64, NoPhi::Int64, energy::Float64)`
        ... defines a angular distribution for given theta and for NoPhi between 0 <= phi < 2pi and a fixed energy
            a dist::SfaAzimuthalAngularDistribution is returned.
    """
    function SfaAzimuthalAngularDistribution(theta::Float64, NoPhi::Int64, energy::Float64)
        phis = Float64[];          for  i=1:NoPhi        push!(phis, i*2pi / NoPhi)                  end
        SfaAzimuthalAngularDistribution(theta, phis, energy)
    end

    function Base.string(obs::SfaAzimuthalAngularDistribution)
        sa = "Compute photoelectron angular distribution at polar angle theta = $(obs.theta) and energy $(obs.energy) for phi:"
        return( sa )
    end

    function Base.show(io::IO, obs::SfaAzimuthalAngularDistribution)
        sa = string(obs);       print(io, sa, "\n");    print(io, "   ", obs.phi)
    end

    """
    `struct  StrongField.SfaPolarAngularDistribution  <: StrongField.AbstractSFAObservable`
        ... to compute in SFA the angular distribution of photoelectrons for given azimuthal angle phi and at given polar angles for a fixed energy

        + phi         ::Float64             ... phi angle of the energy distribution
        + theta       ::Array{Float64,1}    ... specifies the polar angles
        + energy      ::Float64             ... specifies the photoelectron energy in the current units.
    """
    struct   SfaPolarAngularDistribution  <: StrongField.AbstractSFAObservable
        phi           ::Float64
        thetas        ::Array{Float64,1}
        energy        ::Float64
    end

    """
    `StrongField.SfaPolarAngularDistribution(phi::Float64, NoTheta::Int64, energy::Float64)`
        ... defines a angular distribution for given phi and for NoPhi between 0 <= phi < 2pi and a fixed energy
            a dist::SfaPolarAngularDistribution is returned.
    """
    function SfaPolarAngularDistribution(phi::Float64, NoTheta::Int64, energy::Float64)
        thetas = Float64[];          for  i=1:NoTheta        push!(thetas, i*2pi / NoTheta)                  end
        SfaPolarAngularDistribution(phi, thetas, energy)
    end

    function Base.string(obs::SfaPolarAngularDistribution)
        sa = "Compute photoelectron angular distribution at azimuthal angle phi = $(obs.phi) and energy $(obs.energy) for theta:"
        return( sa )
    end

    function Base.show(io::IO, obs::SfaPolarAngularDistribution)
        sa = string(obs);       print(io, sa, "\n");    print(io, "   ", obs.theta)
    end

    """
    `abstract type StrongField.AbstractVolkovState`
        ... defines an abstract and a number of types for dealing with the Volkov states in the computation of the SFA amplitudes.

        + struct FreeVolkov       ... to apply the free-Volkov states.
        + struct CoulombVolkov    ... apply the Coulomb-Volkov states.
        + struct DistortedVolkov  ... apply the Distorted-Volkov states.
    """
    abstract type  AbstractVolkovState                         end
    struct         FreeVolkov        <:  AbstractVolkovState   end

    """
    `struct StrongField.CoulombVolkov     <:  StrongField.AbstractVolkovState`

        + Z ::Float64   ... Charge that generates the Coulomb potential
    """
    struct  CoulombVolkov   <: StrongField.AbstractVolkovState
        Z   ::Float64
    end

    struct         DistortedVolkov   <:  AbstractVolkovState   end

    function Base.string(volkov::FreeVolkov)         return( "Free-Volkov states" )         end
    function Base.string(volkov::CoulombVolkov)      return( "Coulomb-Volkov states in potential of charge Z = $(volkov.Z)" )      end
    function Base.string(volkov::DistortedVolkov)    return( "Distorted-Volkov states" )    end

    """
    `struct  StrongField.Settings`
        ... defines a type for the details and parameters of computing SFA amplitudes.

        + multipoles            ::Array{EmMultipoles}  ... Multipoles of the radiation field that are to be included.
        + gauges                ::Array{UseGauge}      ... Specifies the gauges to be included into the computations.
        + printAmplitudes       ::Bool                 ... True: the amplitudes are printed during computation
        + coupledBasis          ::Bool                 ... True: compute in |j l mj> angular momentum basis; False: compute in |l ml> angular momentum basis
        + hydrogenic            ::Bool                 ... True: Use hydrogenic wave functions for the initial state
        + hydrogenic1s          ::Bool                 ... True: Use hydrogenic wave functions for the initial state
        + mAverage              ::Bool                 ... True: Average over initialState projections mj or ml and sum over continuum spin projection msp
    """
    struct Settings
        multipoles              ::Array{EmMultipole,1}
        gauges                  ::Array{UseGauge}
        printAmplitudes         ::Bool
        coupledBasis            ::Bool
        hydrogenic              ::Bool
        hydrogenic1s            ::Bool
        mAverage                ::Bool
    end

    """
    `StrongField.Settings()`
        ... constructor for the default values of SFA computations.
    """
    function Settings()
        Settings([E1], [UseCoulomb], true, true, false, false, true)
    end

    # `Base.show(io::IO, settings::StrongField.Settings)`  ... prepares a proper printout of the variable settings::StrongField.Settings.
    function Base.show(io::IO, settings::StrongField.Settings)
        println(io, "multipoles:                 $(settings.multipoles)  ")
        println(io, "use-gauges:                 $(settings.gauges)  ")
        println(io, "printBefore:                $(settings.printBefore)  ")
        println(io, "printAmplitudes:            $(settings.printAmplitudes)  ")
        println(io, "coupledBasis:               $(settings.coupledBasis)  ")
        println(io, "hydrogenic:                 $(settings.hydrogenic)  ")
        println(io, "mAverage:                   $(settings.mAverage)  ")
    end

    """
    `struct  StrongField.Computation`
        ... defines a type for the computation of strong-field amplitudes and observables (properties).

        + observable         ::AbstractSFAObservable        ... The obserable to be calculated in SFA.
        + nuclearModel       ::Nuclear.Model                ... Model, charge and parameters of the nucleus.
        + grid               ::Radial.Grid                  ... The radial grid to be used for the computation.
        + initialLevel       ::Level                        ... Initial level of the atom
        + finalLevel         ::Level                        ... Final level of the atom
        + beam               ::Pulse.AbstractBeam           ... Beam character and properties of the incident light pulse.
        + envelope           ::Pulse.AbstractEnvelope       ... Envelope of the incident light pulse.
        + polarization       ::Basics.AbstractPolarization  ... Envelope of the incident light pulse.
        + volkov             ::AbstractVolkovState          ... Specify the treatment/approach of the Volkov states.
        + settings           ::StrongField.Settings         ... Settings for the controlling the SFA amplitudes.
    """
    struct  Computation
        observable           ::AbstractSFAObservable
        nuclearModel         ::Nuclear.Model
        grid                 ::Radial.Grid
        initialLevel         ::Level
        finalLevel           ::Level
        beam                 ::Pulse.AbstractBeam
        envelope             ::Pulse.AbstractEnvelope
        polarization         ::Basics.AbstractPolarization
        volkov               ::AbstractVolkovState
        settings             ::StrongField.Settings
    end

    """
    `StrongField.Computation()`  ... constructor for an `empty` instance of StrongField.Computation().
    """
    function Computation()
        Computation( SfaNoObservable(), Nuclear.Model(1.0), Radial.Grid(true), Orbital("1s_1/2",0.5), Orbital("1s_1/2",0.5),
                     Pulse.PlaneWaveBeam(), Pulse.InfiniteEnvelope(), Basics.RightCircular(), FreeVolkov(), Settings()  )
    end

    # `Base.string(comp::StrongField.Computation)`  ... provides a String notation for the variable comp::StrongField.Computation.
    function Base.string(comp::StrongField.Computation)
        sa = "Strong-field computation:  " * string(comp.observable) * " for Z = $(comp.nuclearModel.Z)  with " * string(comp.volkov) * "\n"
        sa = sa * " initial level with (J, P, energy) = ($(comp.initialLevel.J), $(comp.initialLevel.parity), $(comp.initialLevel.energy)) \n"
        sa = sa * " final level with   (J, P, energy) = ($(comp.finalLevel.J), $(comp.finalLevel.parity), $(comp.finalLevel.energy))   \n\n"
        sa = sa * "The incident laser pulse is described by the: \n "
        sa = sa * string(comp.beam) * "\n " * string(comp.polarization) * "\n " * string(comp.envelope)
        return( sa )
    end

    # `Base.show(io::IO, comp::StrongField.Computation)`  ... prepares a proper printout of the variable comp::StrongField.Computation.
    function Base.show(io::IO, comp::StrongField.Computation)
        sa = Base.string(comp);             print(io, sa, "\n\n")
        println(io, "Settings:              \n$(comp.settings)    ")
        println(io, "nuclearModel:          $(comp.nuclearModel)  ")
        println(io, "grid:                  $(comp.grid)  ")
    end

    """
    `struct  StrongField.SphericalAmplitude`
        ... to keep the amplitude at a given  energy-angular point (energy, theta, phi) in momentum space.

        + energy        ::Float64      ... Kinetic energy of the (outgoing) photoelectron.
        + theta         ::Float64      ... Polar angle of the (outgoing) photoelectron.
        + phi           ::Float64      ... Azimuthal angle of the (outgoing) photoelectron.
        + value         ::ComplexF64   ... (Total) Amplitude.
    """
    struct   SphericalAmplitude
        energy          ::Float64
        theta           ::Float64
        phi             ::Float64
        value           ::ComplexF64
    end

    function Base.string(amp::StrongField.SphericalAmplitude)
        sa = "Total SFA amplitude at (energy, theta, phi) = ($(amp.energy),$(amp.theta),$(amp.phi)) is $(amp.value)."
        return( sa )
    end

    function Base.show(io::IO, amp::StrongField.SphericalAmplitude)
        sa = string(amp);       print(io, sa)
    end

    """
    `struct  StrongField.OutcomeEnergyDistribution`
        ... to comprise the energy distribution of photoelectrons at given energies and angles, for instance,
            for later graphical representation.

        + theta         ::Float64             ... polar angle of the energy distribution
        + phi           ::Float64             ... azimuthal angle of the energy distribution
        + energies      ::Array{Float64,1}    ... selected energies of the distribution.
        + probabilities ::Array{Float64,1}    ... calculated probabilities of the energy distribution.
    """
    struct   OutcomeEnergyDistribution
        theta           ::Float64
        phi             ::Float64
        energies        ::Array{Float64,1}
        probabilities   ::Array{Float64,1}
    end

    # `Base.show(io::IO, outcome::StrongField.OutcomeEnergyDistribution)`  ... prepares a proper printout of the variable outcome::StrongField.OutcomeEnergyDistribution.
    function Base.show(io::IO, outcome::StrongField.OutcomeEnergyDistribution)
        println(io, "theta:             $(outcome.theta)  ")
        println(io, "phi:               $(outcome.phi)  ")
        println(io, "energies:          $(outcome.energies)  ")
        println(io, "probabilities:     $(outcome.probabilities)  ")
    end

    """
    `struct  StrongField.OutcomeMomentumDistribution`
        ... to comprise the momentum distribution of photoelectrons at given energies and angles, for instance,
            for later graphical representation.

        + theta         ::Float64             ... polar angle of the momentum distribution
        + momenta       ::Array{Float64}      ... pairs [px py] of momenta of the momentum distribution
        + probabilities ::Array{Float64,1}    ... calculated probabilities of the momentum distribution.
    """
    struct   OutcomeMomentumDistribution
        theta           ::Float64
        momenta         ::Array{Float64}
        probabilities   ::Array{Float64,1}
    end

    # `Base.show(io::IO, outcome::StrongField.OutcomeMomentumDistribution)`  ... prepares a proper printout of the variable outcome::StrongField.OutcomeMomentumDistribution.
    function Base.show(io::IO, outcome::StrongField.OutcomeMomentumDistribution)
        println(io, "theta:             $(outcome.theta)  ")
        println(io, "momenta:           $(outcome.momenta)  ")
        println(io, "probabilities:     $(outcome.probabilities)  ")
    end

    """
    `struct  StrongField.OutcomeAzimuthalAngularDistribution`
        ... to comprise the angular distribution of photoelectrons at given energy and angles,
            for later graphical representation.

        + theta         ::Float64             ... polar angle of the angular distribution
        + energy        ::Float64             ... photoelectron energy
        + phis          ::Array{Float64,1}      ... azimuthal angles phi
        + probabilities ::Array{Float64,1}    ... calculated probabilities of the angular distribution.
    """
    struct   OutcomeAzimuthalAngularDistribution
        theta           ::Float64
        energy          ::Float64
        phis            ::Array{Float64,1}
        probabilities   ::Array{Float64,1}
    end

    # `Base.show(io::IO, outcome::StrongField.OutcomeAzimuthalAngularDistribution)`  ... prepares a proper printout of the variable outcome::StrongField.OutcomePolarAngularDistribution.
    function Base.show(io::IO, outcome::StrongField.OutcomeAzimuthalAngularDistribution)
        println(io, "theta:             $(outcome.theta)  ")
        println(io, "energy:            $(outcome.energy)  ")
        println(io, "phis:              $(outcome.phis)  ")
        println(io, "probabilities:     $(outcome.probabilities)  ")
    end

    """
    `struct  StrongField.OutcomePolarAngularDistribution`
        ... to comprise the angular distribution of photoelectrons at given energy and angles,
            for later graphical representation.

        + phi           ::Float64             ... azimuthal angle of the angular distribution
        + energy        ::Float64             ... photoelectron energy
        + thetas        ::Array{Float64,1}      ... polar angles phi
        + probabilities ::Array{Float64,1}    ... calculated probabilities of the angular distribution.
    """
    struct   OutcomePolarAngularDistribution
        phi             ::Float64
        energy          ::Float64
        thetas          ::Array{Float64,1}
        probabilities   ::Array{Float64,1}
    end

    # `Base.show(io::IO, outcome::StrongField.OutcomePolarAngularDistribution)`  ... prepares a proper printout of the variable outcome::StrongField.OutcomePolarAngularDistribution.
    function Base.show(io::IO, outcome::StrongField.OutcomePolarAngularDistribution)
        println(io, "phi:             $(outcome.phi)  ")
        println(io, "energy:            $(outcome.energy)  ")
        println(io, "thetas:              $(outcome.thetas)  ")
        println(io, "probabilities:     $(outcome.probabilities)  ")
    end

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

    """
    `StrongField.VolkovP(epsilonp::Float64, lp::Int, rGrid::Array{Float64,1})`
        ... returns the (plane-wave) Volkov radial wave function
                at the radial grid points rGrid
    """
    function  VolkovP(epsilonp::Float64, lp::Int, rGrid::Array{Float64,1})
        P = Float64[]

        for j = 1:length(rGrid)
            r = rGrid[j]
            p = r * GSL.sf_bessel_jl( lp, sqrt(2*epsilonp)*r )
            push!( P, p )
        end

        return( P )
    end

    """
    `StrongField.CoulombVolkovP(epsilonp::Float64, lp::Int, ZContinuum::Float64, rGrid::Array{Float64,1})`
        ... returns the Coulomb-Volkov radial wave function
                at the radial grid points rGrid
    """
    function  CoulombVolkovP(epsilonp::Float64, lp::Int, ZContinuum::Float64, rGrid::Array{Float64,1})
        sqrtTwoEps = sqrt(2*epsilonp)

        P = ComplexF64[]

        for j = 1:length(rGrid)
            r = rGrid[j]
            rho = sqrtTwoEps*r
            eta = -ZContinuum/sqrtTwoEps
            GammaValue = gamma( lp + 1 + eta * im )
            sigma = angle( GammaValue )
            Fl = 2.0^lp * exp(-0.5*pi*eta) * abs( GammaValue ) / factorial(2*lp+1) * rho^(lp+1) * exp( rho * im ) * HypergeometricFunctions.drummond1F1(lp+1+eta*im,2*lp+2,-2*rho*im)

            p = Fl * exp( sigma * im ) / sqrtTwoEps

            push!( P, p )
        end

        return( P )
    end

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

    """
    `StrongField.computeEnvelopeIntegrals(envelope::Pulse.AbstractEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                          thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)`
        ... computes the pulse-envelope integrals for the given laser pulse; this pulse is completely specified by its beam character
            and parameters, the pulse-envelope parameters and the polarization of the light. A tuple of the two integrals
            (fVolkov, fVolkovSquared)::Tuple{Complex{Float64},Complex{Float64}} is returned.
    """
    function computeEnvelopeIntegrals(envelope::Pulse.AbstractEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                      thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)
        fVolkovPlus    = Pulse.envelopeVolkovIntegral(true,  envelope, beam, polarization, thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)
        fVolkovMinus   = Pulse.envelopeVolkovIntegral(false, envelope, beam, polarization, thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)
        fVolkovSquared = Pulse.envelopeQuadVolkovIntegral(envelope, beam, polarization,    thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)

        #println("> Envelope integrals for a $(string(envelope)) with A0=$(beam.A0), omega=$(beam.omega) [a.u.], cep=$(beam.cep)" *
        #        " and for $(string(polarization)) light are:" )
        #println("    F^(Volkov) [+; omega; f^(env); A]   = $fVolkovPlus" )
        #println("    F^(Volkov) [-; omega; f^(env); A]   = $fVolkovMinus" )
        #println("    F^(quad, Volkov) [f^(env); A]       = $fVolkovSquared" )

        return( (fVolkovPlus, fVolkovMinus, fVolkovSquared) )
    end

    """
    `StrongField.computeOutcome(observable::StrongField.SfaEnergyDistribution, amplitudes::Array{SphericalAmplitude,3})`
        ... computes the requested photoelectron energy distribution and returns the results in the variable
            outcome::StrongField.OutcomeEnergyDistribution
    """
    function  computeOutcome(obs::StrongField.SfaEnergyDistribution, amplitudes::Array{SphericalAmplitude,3})
        probabilities = Float64[]
        dimensions = size(amplitudes)
        for  jAmp  = 1:dimensions[3]
            pp = sqrt(2*amplitudes[1,1,jAmp].energy)
            amp = 0.
            for jmi = 1:dimensions[1]
                for jmf = 1:dimensions[2]
                    amp += amplitudes[jmi,jmf,jAmp].value * conj(amplitudes[jmi,jmf,jAmp].value) #sum of probabilities
                end
            end
            push!(probabilities, pp * amp )
        end
        outcome = OutcomeEnergyDistribution(obs.theta, obs.phi, obs.energies, probabilities)
        return( outcome )
    end

    """
    `StrongField.computeOutcome(observable::StrongField.SfaMomentumDistribution, amplitudes::Array{SphericalAmplitude,3})`
        ... computes the requested photoelectron momentum distribution and returns the results in the variable
            outcome::StrongField.OutcomeMomentumDistribution
    """
    function  computeOutcome(obs::StrongField.SfaMomentumDistribution, amplitudes::Array{SphericalAmplitude,3})
        probabilities   = Float64[]
        dimensions      = size(amplitudes)
        momenta         = Array{Float64}(undef,0,2)

        for  jAmp  = 1:dimensions[3]
            pp = sqrt(2*amplitudes[1,1,jAmp].energy)
            px = pp*sin(amplitudes[1,1,jAmp].theta)*cos(amplitudes[1,1,jAmp].phi)
            py = pp*sin(amplitudes[1,1,jAmp].theta)*sin(amplitudes[1,1,jAmp].phi)
            momenta = [momenta ; [px py]]
            amp = 0.
            for jmi = 1:dimensions[1]
                for jmf = 1:dimensions[2]
                    amp += amplitudes[jmi,jmf,jAmp].value * conj(amplitudes[jmi,jmf,jAmp].value) #sum of probabilities
                end
            end
            push!(probabilities, pp * amp )
        end
        outcome = OutcomeMomentumDistribution(obs.theta, momenta, probabilities)
        return( outcome )
    end

    """
    `StrongField.computeOutcome(observable::StrongField.SfaAzimuthalAngularDistribution, amplitudes::Array{SphericalAmplitude,3})`
        ... computes the requested photoelectron angular distribution and returns the results in the variable
            outcome::StrongField.OutcomeAzimuthalAngularDistribution
    """
    function  computeOutcome(obs::StrongField.SfaAzimuthalAngularDistribution, amplitudes::Array{SphericalAmplitude,3})
        probabilities = Float64[]
        dimensions = size(amplitudes)
        for  jAmp  = 1:dimensions[3]
            pp = sqrt(2*amplitudes[1,1,jAmp].energy)
            amp = 0.
            for jmi = 1:dimensions[1]
                for jmf = 1:dimensions[2]
                    amp += amplitudes[jmi,jmf,jAmp].value * conj(amplitudes[jmi,jmf,jAmp].value) #sum of probabilities
                end
            end
            push!(probabilities, pp * amp )
        end
        outcome = OutcomeAzimuthalAngularDistribution(obs.theta, obs.energy, obs.phis, probabilities)
        return( outcome )
    end

    """
    `StrongField.computeOutcome(observable::StrongField.SfaPolarAngularDistribution, amplitudes::Array{SphericalAmplitude,3})`
        ... computes the requested photoelectron angular distribution and returns the results in the variable
            outcome::StrongField.OutcomePolarAngularDistribution
    """
    function  computeOutcome(obs::StrongField.SfaPolarAngularDistribution, amplitudes::Array{SphericalAmplitude,3})
        probabilities = Float64[]
        dimensions = size(amplitudes)
        for  jAmp  = 1:dimensions[3]
            pp = sqrt(2*amplitudes[1,1,jAmp].energy)
            amp = 0.
            for jmi = 1:dimensions[1]
                for jmf = 1:dimensions[2]
                    amp += amplitudes[jmi,jmf,jAmp].value * conj(amplitudes[jmi,jmf,jAmp].value) #sum of probabilities
                end
            end
            push!(probabilities, pp * amp )
        end
        outcome = OutcomePolarAngularDistribution(obs.phi, obs.energy, obs.thetas, probabilities)
        return( outcome )
    end

    #---------------------------------------------------------------------------------------------------
    #--------------------------------------FORMULATION in j-l-mj-BASIS:---------------------------------
    #---------------------------------------------------------------------------------------------------

    """
    `StrongField.pReducedME(Pepsplp::Array{ComplexF64,1}, lp::Int, jp::Float64, n::Int, l::Int, j::Float64, initialOrbital::Orbital, grid::Radial.Grid)`
        ... computes the reduced matrix elements of the momentum operator <epsilonp lp jp||p||n l j> in the one-particle picture
            where the radial wave function Pepsplp of |epsilonp lp jp> needs to be provided as an argument on the grid rGrid
    """
    function  pReducedME(Pepsplp::Array{ComplexF64,1}, lp::Int, jp::Float64, l::Int, j::Float64, initialOrbital::Orbital, grid::Radial.Grid)

        #---- Compute the factor in <epsilonp lp jp||p||n l j> = fac * <epsilonp lp||p||n l>-----
        fac = 0.
        for    kms = -1:0 #ms = 1/2 + kms
            ms = 1/2 + kms
            for    ml = -l:l
                mj = ml + ms
                for    mlp = -lp:lp
                    mjp = mlp + ms
                    fac = fac + AngularMomentum.ClebschGordan_old(  AngularJ64(Rational(lp)),   AngularM64(Rational(mlp)) ,
                                                                    AngularJ64(1//2),           AngularM64(Rational(ms)),
                                                                    AngularJ64(Rational(jp)),   AngularM64(Rational(mjp))       )       *
                                AngularMomentum.ClebschGordan_old(  AngularJ64(Rational(l)),    AngularM64(Rational(ml)),
                                                                    AngularJ64(1//2),           AngularM64(Rational(ms)),
                                                                    AngularJ64(Rational(j)),    AngularM64(Rational(mj))       )
                end
            end
        end
        #----------------------------------------------------------------------------------------

        #--------------------------------Compute <epsilonp lp||p||n l>---------------------------
        rvalues = grid.r
        rweights = grid.wr
        orderGL = min(size(Pepsplp, 1), size(initialOrbital.P, 1))

        integral = 0. * im

        #Sum over grid and compute Gauss-Legendre sum
        for    k = 1:orderGL
            r = rvalues[k]
            integrand = conj( Pepsplp[k] )/r * ( r*initialOrbital.Pprime[k] - ((lp-l)*(lp+l+1))/2 * initialOrbital.P[k] )

            #Gauss-Legendre sum
            integral = integral + rweights[k] * integrand
        end

        #Note that GSL.sf_coupling_3j takes takes the input (2*j1,2*j2,2*j3,2*m1,2*m2,2*m3)
        integral = integral * (-im)^(lp+1) * (-1)^lp * GSL.sf_coupling_3j( 2*lp, 2*1, 2*l, 0, 0, 0 )
        #----------------------------------------------------------------------------------------

        return( fac * integral )
    end

    """
    `StrongField.scalarProdBoundCont(Pepsplp::Array{ComplexF64,1}, l::Int, initialOrbital::Orbital, grid::Radial.Grid)`
        ... computes the scalar product of the bound and continuum states in the one-particle picture
            where the radial wave function Pepsplp of |epsilonp lp jp> needs to be provided as an argument on the grid rGrid
    """
    function  scalarProdBoundCont(Pepsplp::Array{ComplexF64,1}, l::Int, initialOrbital::Orbital, grid::Radial.Grid)
        rvalues = grid.r
        rweights = grid.wr

        orderGL = min(size(Pepsplp, 1), size(initialOrbital.P, 1))

        integral = 0. * im

        #Sum over grid and compute Gauss-Legendre sum
        for    k = 1:orderGL
            r = rvalues[k]
            integrand = conj( Pepsplp[k] ) * initialOrbital.P[k]

            #Gauss-Legendre sum
            integral = integral + rweights[k] * integrand
        end

        return((-im)^l * integral)
    end

    """
    `StrongField.computeSphericalAmplitudes(comp::StrongField.Computation)`
        ... computes all necessary spherical amplitudes for the given initial and final levels, the Volkov states
            and polarization and by using the pulse envelope integrals - using a coupled angular-momentum basis for the active electron,
            A sfaAmplitudes::Array{SphericalAmplitude,3} is returned.
    """
    function  computeSphericalAmplitudes(comp::StrongField.Computation)

        #Extract the initial orbital of the active electron from the many-electron comp.initialLevel and set quantum numbers
        initialOrbitals = comp.initialLevel.basis.orbitals

        #Find highest lying orbital (smallest ionization potential)
        defaultSubshell = [sh for (sh,or) in initialOrbitals][1] #This is not nice; must be a better way to simply get a default element from a Dict
        initialOrbital = initialOrbitals[defaultSubshell]
        minIonPotential = abs(initialOrbital.energy)
        for (subshell,orbital) in initialOrbitals
            if   abs(orbital.energy) < minIonPotential
                initialOrbital = orbital
                minIonPotential = abs(orbital.energy)
            end
        end

        initialEneV = convertUnits("energy: from atomic to eV", initialOrbital.energy)
        println("")
        println("Selected orbital $(initialOrbital.subshell) as active electron for the strong-field ionization process. Binding energy: $(initialEneV) eV.")
        println("")

        ls = LevelSymmetry(initialOrbital.subshell)
        n = initialOrbital.subshell.n;      l = (ls.J.num+1)/2;    j = ls.J.num/2;     initialEn = initialOrbital.energy

        if  (sign((-1)^l) == -1 && ls.parity == plus::Parity) || (sign((-1)^l) == 1 && ls.parity == minus::Parity)
            l = l - 1
        end
        l = floor(Int,l)

        if comp.settings.hydrogenic && comp.settings.hydrogenic1s
            n = 1
            l = 0
            j = 1/2
        end

        lpMin = abs(l-1)
        lpMax = l+1

        #Determine which spherical SFA amplitudes need to be computed before the actual computation starts
        mjNum = 1; mspNum = 1    #If comp.Average == false
        mjAverageFactor = 1.     #If comp.Average == false
        if( comp.settings.mAverage )
            mjNum = Int(2*j + 1)
            mspNum = 2
            mjAverageFactor = 1.0/sqrt( mjNum )
        end

        sfaAmplitudes = StrongField.determineSphericalAmplitudes(comp.observable,mjNum,mspNum)

        #-----Generate continuum wave functions for all energies and values of lp that are needed below-----

        println("")
        println("Computing continuum wave functions...")

        energies = Array{Float64}(undef,length(sfaAmplitudes[1,1,:]))
        for  jAmp = 1:length(sfaAmplitudes[1,1,:])  #Get all different energy values
            energies[jAmp] = sfaAmplitudes[1,1,jAmp].energy
        end
        energies = unique(energies)
        lpArray = unique( [lpMin l lpMax] )

        Pepsplp = Array{ComplexF64}(undef,size(energies)[1],size(lpArray)[1],size(comp.grid.r)[1])

        for jlp = 1:size(lpArray)[1]

            lp = lpArray[jlp]

            for jE = 1:size(energies)[1]

                epsilonp = energies[jE]

                if  typeof(comp.volkov) == FreeVolkov

                    Pepsplp[jE,jlp,:] = VolkovP( epsilonp, lp, comp.grid.r )

                elseif  typeof(comp.volkov) == CoulombVolkov

                    Pepsplp[jE,jlp,:] = CoulombVolkovP( epsilonp, lp, comp.volkov.Z, comp.grid.r )

                elseif  typeof(comp.volkov) == DistortedVolkov

                    eta = (-1)^(lp+1)
                    kappa = eta*(lp+1/2)-1/2

                    nrContinuum = Continuum.gridConsistency(epsilonp, comp.grid)
                    contSettings = Continuum.Settings(false, comp.grid.NoPoints)

                    newiLevel = Basics.generateLevelWithSymmetryReducedBasis(comp.initialLevel, comp.initialLevel.basis.subshells)
                    newfLevel = Basics.generateLevelWithSymmetryReducedBasis(comp.finalLevel, newiLevel.basis.subshells)
                    newiLevel = Basics.generateLevelWithExtraSubshell(Subshell(101, kappa), newiLevel)
                    cOrbital, phase  = Continuum.generateOrbitalForLevel(epsilonp, Subshell(101, kappa), newfLevel, comp.nuclearModel, comp.grid, contSettings)

                    jr = size(cOrbital.P)[1]
                    while jr < size(comp.grid.r)[1]
                        push!(cOrbital.P,0)
                        jr += 1
                    end
                    Pepsplp[jE,jlp,:] = cOrbital.P * exp(im*phase)

                end

            end #end for energy

        end #end for lp

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

        println("")
        println("Computing SFA amplitudes...")
        println("")

        # Compute the requested amplitudes
        for jmj = 1:mjNum
            mj = -j + (jmj-1)
            for jmsp = 1:mspNum
                msp = -1/2 + (jmsp-1)
                for  jAmp = 1:length(sfaAmplitudes[jmj,jmsp,:]) #amp in  sfaAmplitudes[jmj,jmsp]
                    amp = sfaAmplitudes[jmj,jmsp,jAmp]
                    thetap        = amp.theta;   phip = amp.phi;
                    energyp = amp.energy; jE = findall(x->x==energyp, vec(energies))[1]
                    envIntegrals  = StrongField.computeEnvelopeIntegrals(comp.envelope, comp.beam, comp.polarization, thetap, phip, energyp, initialEn)
                    fVolkovPlus, fVolkovMinus, fVolkovSquared = envIntegrals
                    #

                    wminus = 0.0im;    wplus = 0.0im
                    # Collect contributions from all l_p, j_p, q terms
                    for  jlp = 1:size(lpArray)[1]
                        lp = lpArray[jlp]

                        kjpMin = floor(Int, ( 2 * abs(lp-1/2) - 1 ) / 2) # jp = (2*kjp + 1)/2, so that the sum below only runs over the correct values. THERE IS A BETTER WAY
                        kjpMax = floor(Int, ( 2 * abs(lp+1/2) - 1 ) / 2)

                        #Compute the sum over total angular momentum jp
                        for kjp = kjpMin:kjpMax
                            jp = 1/2 * (2 * kjp + 1)

                            KWignerEckart = 1.0/sqrt(2*jp+1) #Factor in the Wigner-Eckart theorem (CONVENTION that needs to fit the reduced ME!)

                            if comp.settings.hydrogenic
                                reducedME = pReducedMEHydrogenic(Pepsplp[jE,jlp,:], lp, jp, initialEn, n, l, j, comp.grid)
                            else
                                reducedME = pReducedME(Pepsplp[jE,jlp,:], lp, jp, l, j, initialOrbital, comp.grid)
                            end

                            #Compute the sum over spherical basis components
                            #Note: sphericalYlm needs an integer in the second argument. mj-msp-q is always a whole number, so we can convert to Int using floor without problems
                            for  q in [-1, 0, 1]
                                wplus = wplus +  KWignerEckart * (-1)^q * Basics.determinePolarizationVector(q, comp.polarization, star=false)   *
                                                    AngularMomentum.sphericalYlm(lp, floor(Int,mj-msp-q), amp.theta, amp.phi)                           *
                                                    AngularMomentum.ClebschGordan_old(  AngularJ64(Rational(lp)),   AngularM64(Rational(mj-msp-q)),
                                                                                AngularJ64(1//2),           AngularM64(Rational(msp)),
                                                                                AngularJ64(Rational(jp)),   AngularM64(Rational(mj-q)) )        *
                                                    AngularMomentum.ClebschGordan_old(  AngularJ64(Rational(j)),    AngularM64(Rational(mj)),
                                                                                AngularJ64(1),              AngularM64(Rational(-q)),
                                                                                AngularJ64(Rational(jp)),   AngularM64(Rational(mj-q)) )         *
                                                    reducedME

                                wminus = wminus +  KWignerEckart * Basics.determinePolarizationVector(q, comp.polarization, star=true)           *
                                                    AngularMomentum.sphericalYlm(lp, floor(Int,mj-msp+q), amp.theta, amp.phi)               *
                                                    AngularMomentum.ClebschGordan_old(  AngularJ64(Rational(lp)),   AngularM64(Rational(mj-msp+q)),
                                                                                AngularJ64(1//2),           AngularM64(Rational(msp)),
                                                                                AngularJ64(Rational(jp)),   AngularM64(Rational(mj+q)) )        *
                                                    AngularMomentum.ClebschGordan_old(  AngularJ64(Rational(j)),    AngularM64(Rational(mj)),
                                                                                AngularJ64(1),              AngularM64(Rational(q)),
                                                                                AngularJ64(Rational(jp)),   AngularM64(Rational(mj+q)) )         *
                                                    reducedME
                            end

                        end # for jp

                    end # for lp

                    #Compute the scalar product in the term proportional to F_2 in the spherical amplitude
                    if comp.settings.hydrogenic
                        jlp = findall(x->x==l, vec(lpArray))[1]
                        scalarProd = scalarProdBoundContHydrogenic(Pepsplp[jE,jlp,:], initialEn, n, l, comp.grid)
                    else
                        jlp = findall(x->x==l, vec(lpArray))[1]
                        scalarProd = scalarProdBoundCont(Pepsplp[jE,jlp,:], l, initialOrbital, comp.grid)
                    end

                    #Compute the total amplitude
                    wa = (fVolkovPlus * wminus + fVolkovMinus * wplus) + fVolkovSquared * AngularMomentum.sphericalYlm(l, floor(Int,mj-msp), amp.theta, amp.phi) *
                                                                            AngularMomentum.ClebschGordan_old(    AngularJ64(Rational(l)),    AngularM64(Rational(mj-msp)),
                                                                                AngularJ64(1//2),           AngularM64( Rational(msp) ),
                                                                                AngularJ64(Rational(j)),    AngularM64(Rational(mj))       )  *
                                                                            scalarProd

                    wa = - im/sqrt(2*pi) * wa
                    wa = mjAverageFactor * wa

                    sfaAmplitudes[jmj,jmsp,jAmp] = SphericalAmplitude(amp.energy, amp.theta, amp.phi, wa)

                    if( comp.settings.printAmplitudes )
                        println(">> $(SphericalAmplitude(amp.energy, amp.theta, amp.phi, wa))")
                    end
                end #for amp
            end #for jmsp
        end #for jmj

        return( sfaAmplitudes )
    end

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

    """
    `StrongField.determineSphericalAmplitudes(observable::StrongField.SfaEnergyDistribution, miNum::Int64, mfNum::Int64)`
        ... determines which direct (and other) SFA amplitudes need to be computed; these amplitudes are not yet computed here but can be arranged
            so that only a minimum number of Volkov states and/or reduced many-electron matrix elements need to be computed.
            A list of amplitudes::Array{StrongField.SphericalAmplitude,1} is returned.
    """
    function  determineSphericalAmplitudes(observable::StrongField.SfaEnergyDistribution, miNum::Int64, mfNum::Int64 )
        ENum = length(observable.energies)
        amplitudes = Array{StrongField.SphericalAmplitude}(undef, miNum, mfNum, ENum)

        for jmi = 1:miNum
            for jmf = 1:mfNum
                for jE = 1:ENum
                    energy = observable.energies[jE]
                    amplitudes[jmi,jmf,jE] = SphericalAmplitude(observable.energies[jE], observable.theta, observable.phi, 0.)
                end
            end
        end

        println("> A total of $(length(amplitudes)) spherical amplitudes need to be calculated.")
        return( amplitudes )
    end

    """
    `StrongField.determineSphericalAmplitudes(observable::StrongField.SfaMomentumDistribution, miNum::Int64, mfNum::Int64)`
        ... determines which direct (and other) SFA amplitudes need to be computed; these amplitudes are not yet computed here but can be arranged
            so that only a minimum number of Volkov states and/or reduced many-electron matrix elements need to be computed.
            A list of amplitudes::Array{StrongField.SphericalAmplitude,1} is returned.
    """
    function  determineSphericalAmplitudes(observable::StrongField.SfaMomentumDistribution, miNum::Int64, mfNum::Int64)
        ENum = length(observable.energies)
        phiNum = length(observable.phis)

        amplitudes = Array{StrongField.SphericalAmplitude}(undef, miNum, mfNum, ENum*phiNum)

        for jmi = 1:miNum
            for jmf = 1:mfNum
                ctr = 1
                for  jE = 1:ENum
                    for  jphi = 1:phiNum
                        amplitudes[jmi,jmf,ctr] = SphericalAmplitude(observable.energies[jE], observable.theta, observable.phis[jphi], 0.)
                        ctr = ctr + 1
                    end
                end
            end
        end

        println("> A total of $(length(amplitudes)) spherical amplitudes need to be calculated.")
        return( amplitudes )
    end

    """
    `StrongField.determineSphericalAmplitudes(observable::StrongField.SfaAzimuthalAngularDistribution, miNum::Int64, mfNum::Int64)`
        ... determines which direct (and other) SFA amplitudes need to be computed; these amplitudes are not yet computed here but can be arranged
            so that only a minimum number of Volkov states and/or reduced many-electron matrix elements need to be computed.
            A list of amplitudes::Array{StrongField.SphericalAmplitude,1} is returned.
    """
    function  determineSphericalAmplitudes(observable::StrongField.SfaAzimuthalAngularDistribution, miNum::Int64, mfNum::Int64)
        phiNum = length(observable.phis)
        amplitudes = Array{StrongField.SphericalAmplitude}(undef, miNum, mfNum, phiNum)

        for jmi = 1:miNum
            for jmf = 1:mfNum
                for  jphi = 1:phiNum
                    amplitudes[jmi,jmf,jphi] = SphericalAmplitude(observable.energy, observable.theta, observable.phis[jphi], 0.)
                end
            end
        end

        println("> A total of $(length(amplitudes)) spherical amplitudes need to be calculated.")
        return( amplitudes )
    end

    """
    `StrongField.determineSphericalAmplitudes(observable::StrongField.SfaPolarAngularDistribution, miNum::Int64, mfNum::Int64)`
        ... determines which direct (and other) SFA amplitudes need to be computed; these amplitudes are not yet computed here but can be arranged
            so that only a minimum number of Volkov states and/or reduced many-electron matrix elements need to be computed.
            A list of amplitudes::Array{StrongField.SphericalAmplitude,1} is returned.
    """
    function  determineSphericalAmplitudes(observable::StrongField.SfaPolarAngularDistribution, miNum::Int64, mfNum::Int64)
        thetaNum = length(observable.thetas)
        amplitudes = Array{StrongField.SphericalAmplitude}(undef, miNum, mfNum, thetaNum)

        for jmi = 1:miNum
            for jmf = 1:mfNum
                for  jtheta = 1:thetaNum
                    amplitudes[jmi,jmf,thetaNum] = SphericalAmplitude(observable.energy, observable.thetas[jtheta], observable.phi, 0.)
                end
            end
        end

        println("> A total of $(length(amplitudes)) spherical amplitudes need to be calculated.")
        return( amplitudes )
    end

    """
    `StrongField.perform(comp::StrongField.Computation; output::Bool=false)`
        ... to perform a computation of (one) selected observable that is related to SFA, such as a photoelectron energy
            distribution, momentum distribution or the computation of sidebands.
    """
    function perform(comp::StrongField.Computation; output::Bool=false)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        nModel = comp.nuclearModel

        if       typeof(comp.observable) == StrongField.SfaEnergyDistribution
            if  comp.settings.coupledBasis
                sfaAmplitudes = StrongField.computeSphericalAmplitudes(comp)
            else
                sfaAmplitudes = StrongField.computeSphericalAmplitudesUncoupled(comp)
            end
            sfaOutcome    = StrongField.computeOutcome(comp.observable, sfaAmplitudes)
            if output    results = Base.merge( results, Dict("computation" => comp, "energy distribution" => sfaOutcome) )  end
        elseif   typeof(comp.observable) == StrongField.SfaMomentumDistribution
            if  comp.settings.coupledBasis
                sfaAmplitudes = StrongField.computeSphericalAmplitudes(comp)
            else
                sfaAmplitudes = StrongField.computeSphericalAmplitudesUncoupled(comp)
            end
            sfaOutcome    = StrongField.computeOutcome(comp.observable, sfaAmplitudes)
            if output    results = Base.merge( results, Dict("computation" => comp, "momentum distribution" => sfaOutcome) )  end
        elseif   typeof(comp.observable) == StrongField.SfaAzimuthalAngularDistribution
            if  comp.settings.coupledBasis
                sfaAmplitudes = StrongField.computeSphericalAmplitudes(comp)
            else
                sfaAmplitudes = StrongField.computeSphericalAmplitudesUncoupled(comp)
            end
            sfaOutcome    = StrongField.computeOutcome(comp.observable, sfaAmplitudes)
            if output    results = Base.merge( results, Dict("computation" => comp, "angular distribution" => sfaOutcome) )  end
        elseif   typeof(comp.observable) == StrongField.SfaPolarAngularDistribution
            if  comp.settings.coupledBasis
                sfaAmplitudes = StrongField.computeSphericalAmplitudes(comp)
            else
                sfaAmplitudes = StrongField.computeSphericalAmplitudesUncoupled(comp)
            end
            sfaOutcome    = StrongField.computeOutcome(comp.observable, sfaAmplitudes)
            if output    results = Base.merge( results, Dict("computation" => comp, "angular distribution" => sfaOutcome) )  end
        else     error("Undefined observable for strong-field computations.")
        end

        if  output   return( results )   else    return( nothing )   end
    end

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

    include("module-StrongField-inc-hydrogenic.jl") #StrongField routines for hydrogenic initial states
    include("module-StrongField-inc-uncoupled.jl")  #StrongField routines for computations in the basis of "uncoupled" angular momentum l,ml instead of j,l,mj
    ## include("module-StrongField-inc-postProcessing.jl")       #StrongField.plot(...) routine to plot results of StrongField computations

end # module


text:
@symbol_func function cross_sectional_area(cur_reactor::AbstractReactor)
  cur_volume = cur_reactor.pi

  cur_volume *= a(cur_reactor) ^ 2

  cur_volume *= cur_reactor.kappa_95

  cur_volume *= g(cur_reactor)

  cur_volume
end

ext: jl


text:
const bullet_lib_path = joinpath(ENV["BULLET3"], "src");
const bullet_objf_path = joinpath(ENV["BULLET3"], "bullet-build/src");
const bullet_include_paths = [bullet_lib_path, bullet_objf_path, dirname(@__FILE__())]  # TODO(tzahroof): add the header and cpp(? maybe .o? Maybe look into this) files

for p in bullet_include_paths
  addHeaderDir(p, kind=C_System)
end

# add *.so files from /usr/local/lib
include_files = ["libBulletCollision.so","libLinearMath.so","libBulletDynamics.so"];
so_hdr = "/usr/lib/x86_64-linux-gnu"
for fn in include_files
  Libdl.dlopen(joinpath(so_hdr,fn), Libdl.RTLD_GLOBAL);
end

cxx"""
#include <iostream>
#include <sstream>
#include "btBulletCollisionCommon.h"
#include "BulletCollision/CollisionShapes/btConvex2dShape.h"
#include "BulletCollision/NarrowPhaseCollision/btGjkPairDetector.h"
#include "BulletCollision/NarrowPhaseCollision/btGjkEpaPenetrationDepthSolver.h"
#include "BulletCollision/NarrowPhaseCollision/btPointCollector.h"
#include "BulletDynamics/Featherstone/btMultiBody.h"
#include "distanceComputation.h"
#include "LinearMath/btConvexHull.h"

#ifndef BT_BINARY_COLLISION_CALLBACKS
#define BT_BINARY_COLLISION_CALLBACKS

struct BinaryCollisionCallback : public btCollisionWorld::ContactResultCallback
{
  bool is_collision;

  BinaryCollisionCallback() : btCollisionWorld::ContactResultCallback(), is_collision(false) {}

  virtual btScalar addSingleResult(btManifoldPoint& cp,
    const btCollisionObjectWrapper* colObj0Wrap, int partId0, int index0,
    const btCollisionObjectWrapper* colObj1Wrap, int partId1, int index1) {
    is_collision = true;
    return 0;
  }
};

struct BinarySweptConvexCollisionCallback : public btCollisionWorld::ConvexResultCallback
{
  bool is_collision;

  BinarySweptConvexCollisionCallback() : btCollisionWorld::ConvexResultCallback(), is_collision(false) {}

  virtual btScalar addSingleResult(btCollisionWorld::LocalConvexResult& convexResult, bool normalInWorldSpace) {
    is_collision = true;
    return 0;
  }
};
#endif // BT_BINARY_COLLISION_CALLBACKS
"""


text:
mutable struct DistSparseMatrix{T} <: ElementalMatrix{T}
    obj::Ptr{Cvoid}
end

for (elty, ext) in ((:ElInt, :i),
                    (:Float32, :s),
                    (:Float64, :d),
                    (:ComplexF32, :c),
                    (:ComplexF64, :z))
    @eval begin

        # destructor to be used in finalizer. Don't call explicitly
        function destroy(A::DistSparseMatrix{$elty})
            ElError(ccall(($(string("ElDistSparseMatrixDestroy_", ext)), libEl), Cuint,
                (Ptr{Cvoid},), A.obj))
            return nothing
        end

        function DistSparseMatrix(::Type{$elty}, comm::ElComm = MPI.CommWorld[])
            obj = Ref{Ptr{Cvoid}}(C_NULL)
            ElError(ccall(($(string("ElDistSparseMatrixCreate_", ext)), libEl), Cuint,
                (Ref{Ptr{Cvoid}}, ElComm),
                obj, comm))
            A = DistSparseMatrix{$elty}(obj[])
            finalizer(destroy, A)
            return A
        end

        function comm(A::DistSparseMatrix{$elty})
            cm = Ref{ElComm}()
            ElError(ccall(($(string("ElDistSparseMatrixComm_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ref{ElComm}),
                A.obj, cm))
            return cm[]
        end

        function globalRow(A::DistSparseMatrix{$elty}, iLoc::Integer)
            i = Ref{ElInt}(0)
            ElError(ccall(($(string("ElDistSparseMatrixGlobalRow_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, Ref{ElInt}),
                A.obj, iLoc-1, i))
            return i[]+1
        end

        function height(A::DistSparseMatrix{$elty})
            i = Ref{ElInt}(0)
            ElError(ccall(($(string("ElDistSparseMatrixHeight_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ref{ElInt}),
                A.obj, i))
            return i[]
        end

        function localHeight(A::DistSparseMatrix{$elty})
            i = Ref{ElInt}(0)
            ElError(ccall(($(string("ElDistSparseMatrixLocalHeight_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ref{ElInt}),
                A.obj, i))
            return i[]
        end

        function numLocalEntries(A::DistSparseMatrix{$elty})
            n = Ref{ElInt}(0)
            ElError(ccall(($(string("ElDistSparseMatrixNumLocalEntries_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ref{ElInt}),
                A.obj, n))
            return n[]
        end

        function processQueues(A::DistSparseMatrix{$elty})
            ElError(ccall(($(string("ElDistSparseMatrixProcessQueues_", ext)), libEl), Cuint,
                (Ptr{Cvoid},),
                A.obj))
            return nothing
        end

        function queueLocalUpdate(A::DistSparseMatrix{$elty}, localRow::Integer, col::Integer, value::$elty)
            ElError(ccall(($(string("ElDistSparseMatrixQueueLocalUpdate_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, ElInt, $elty),
                A.obj, localRow-1, col-1, value))
            return nothing
        end

        function queueUpdate(A::DistSparseMatrix{$elty}, row::Integer, col::Integer, value::$elty, passive::Bool = true)
            ElError(ccall(($(string("ElDistSparseMatrixQueueUpdate_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, ElInt, $elty, Bool),
                A.obj, row-1, col-1, value, passive))
            return nothing
        end

        function reserve(A::DistSparseMatrix{$elty}, numLocalEntries::Integer, numRemoteEntries::Integer = 0)
            ElError(ccall(($(string("ElDistSparseMatrixReserve_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, ElInt),
                A.obj, numLocalEntries, numRemoteEntries))
            return nothing
        end

        function resize!(A::DistSparseMatrix{$elty}, height::Integer, width::Integer = 1) # to mimic vector behavior
            ElError(ccall(($(string("ElDistSparseMatrixResize_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, ElInt),
                A.obj, height, width))
            return A
        end

        function width(A::DistSparseMatrix{$elty})
            i = Ref{ElInt}(0)
            ElError(ccall(($(string("ElDistSparseMatrixWidth_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ref{ElInt}),
                A.obj, i))
            return i[]
        end
    end
end

# The other constructors don't have a version with dimensions. Should they, or should this one go?
function DistSparseMatrix(::Type{T}, m::Integer, n::Integer, comm::ElComm = MPI.CommWorld[]) where {T}
    A = DistSparseMatrix(T, comm)
    resize!(A, m, n)
    return A
end

# Julia convenience

function Base.show(io::IO, ::MIME"text/plain", A::DistSparseMatrix)
    print(io, "$(size(A, 1))x$(size(A, 2)) $(typeof(A))")
end


text:
puzzleinput = readlines("inputs/day13.txt")

earliesttimestamp = parse(Int, puzzleinput[1])

timetable = map(x -> tryparse(Int64, x), split(puzzleinput[2], ","))

# Firt Half

dt(id::Nothing, ets) = typemax(Int64)
dt(id::Int64, ets) = ((ets ÷ id + 1) * id) - ets
# Calculates the time to next bus id from an earliest time stamp ets

waittimes = [[t dt(t, earliesttimestamp)] for t in timetable]
waittimes = reduce(vcat, waittimes)
# Matrix with in first column bus ids and
# in second colums wait times from the earliest time stamp

waittimes = waittimes[sortperm(waittimes[:, 2]), :]
# Sort by wait times

waittimes[1, 1] * waittimes[1, 2]

# Second Half

#= =======================================================================
This solution cames from an original idea of a friend of mines,
Matteo Bernardini:

https://github.com/matteberna/advent-of-code-2020/blob/master/Solutions/Day13.cs

Thx Whiskee!
======================================================================== =#

function findets(tt)
    ids = filter((!) ∘ isnothing, tt)

    departs = [i - 1 for i = 1:length(tt) if !isnothing(tt[i])]
    # Minutes after 0 at which each bus in the time table has to depart

    dt = 1
    t = ids[1]

    for (i, id) in enumerate(ids)
        while true
            if (t + departs[i]) % id == 0
                dt *= id
                # Given the fact that ids are all prime numbers:
                # if the depart of bus ids[i] is ok a this time stamp t
                # then the times stamp for the subsequent in
                # the time table (eg. ids[i+1]) has to be necessarily
                # t + n * ids[i] * ids[i+1] + departs[i],
                # if we want it to be ok also for the previous ids

                break
            else
                t += dt
            end
        end
    end

    return t
end

@time findets(timetable)

# With input:
# 17,x,x,x,x,x,x,41,x,x,x,37,x,x,x,x,x,367,x,x,x,x,x,x,x,x,x,x,x,x,x,x,x,x,x,x,19,x,x,x,23,x,x,x,x,x,29,x,613,x,x,x,x,x,x,x,x,x,x,x,x,13
# in Wolfram Alpha you could type this to obtain the solution
#  {41x_2-17x_1=7;37x_3-17x_1=11;367x_4-17 x_1=17;19x_5- 17x_1=36;23x_6-17x_1=40;29x_7-17 x_1=46;613 x_8-17x_1=48;13 x_9-17x_1=61} integer solution


text:
#==============================================================================#
# AWSS3/test/runtests.jl
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

using AWSS3
using Base.Test
using Retry

AWSCore.set_debug_level(1)

function test_without_catch(f)

    @protected try
        f()
    catch e
        @ignore if isa(e, Base.Test.Error)
            e = e.err
            rethrow(e)
        end
    end
end

#-------------------------------------------------------------------------------
# Load credentials...
#-------------------------------------------------------------------------------

aws = AWSCore.default_aws_config()
aws[:region] = "ap-southeast-2"

#-------------------------------------------------------------------------------
# S3 tests
#-------------------------------------------------------------------------------

# Delete old test files...

for b in s3_list_buckets()

    if ismatch(r"^ocaws.jl.test", b)
        @protected try
            println("Cleaning up old test bucket: " * b)
            @sync for v in s3_list_versions(aws, b)
                @async s3_delete(aws, b, v["Key"]; version = v["VersionId"])
            end
            s3_delete_bucket(aws, b)
        catch e
            @ignore if isa(e, AWSCore.AWSException) &&
                       e.code == "NoSuchBucket" end
        end
    end
end

import AWSS3.HTTP
HTTP.ConnectionPool.showpool(STDOUT)

# Temporary bucket name...

bucket_name = "ocaws.jl.test." * lowercase(Dates.format(now(Dates.UTC),
                                                        "yyyymmddTHHMMSSZ"))

# Test exception code for deleting non existent bucket...

@protected try

    s3_delete_bucket(aws, bucket_name)

catch e
     @ignore if isa(e, AWSCore.AWSException) &&
                e.code == "NoSuchBucket" end
end

# Create bucket...

s3_create_bucket(aws, bucket_name)

@repeat 4 try

    # Turn on object versioning for this bucket...

    s3_enable_versioning(aws, bucket_name)

    # Check that the new bucket is returned in the list of buckets...

    @test bucket_name in s3_list_buckets(aws)

    test_without_catch() do
        # Check that our test keys do not exist yet...
        @test !s3_exists(aws, bucket_name, "key 1")
        @test !s3_exists(aws, bucket_name, "key2")
        @test !s3_exists(aws, bucket_name, "key3")
    end

catch e
    @delay_retry if isa(e, AWSCore.AWSException) &&
                    e.code == "NoSuchBucket" end
end

# Bucket tagging...

@test isempty(s3_get_tags(aws, bucket_name))
tags = Dict("A" => "1", "B" => "2", "C" => "3")
s3_put_tags(aws, bucket_name, tags)
@test s3_get_tags(aws, bucket_name) == tags
s3_delete_tags(aws, bucket_name)
@test isempty(s3_get_tags(aws, bucket_name))

# Create test objects...

s3_put(aws, bucket_name, "key 1", "data1.v1")
s3_put(bucket_name, "key2", "data2.v1", tags = Dict("Key" => "Value"))
s3_put(aws, bucket_name, "key3", "data3.v1")
s3_put(aws, bucket_name, "key3", "data3.v2")
s3_put(aws, bucket_name, "key3", "data3.v3"; metadata = Dict("foo" => "bar"))
s3_put_tags(aws, bucket_name, "key3", Dict("Left" => "Right"))

@test isempty(s3_get_tags(aws, bucket_name, "key 1"))
@test s3_get_tags(aws, bucket_name, "key2")["Key"] == "Value"
@test s3_get_tags(aws, bucket_name, "key3")["Left"] == "Right"
s3_delete_tags(aws, bucket_name, "key2")
@test isempty(s3_get_tags(aws, bucket_name, "key2"))

# Check that test objects have expected content...

@test s3_get(aws, bucket_name, "key 1") == b"data1.v1"
@test s3_get(aws, bucket_name, "key2") == b"data2.v1"
@test s3_get(bucket_name, "key3") == b"data3.v3"
@test s3_get_meta(bucket_name, "key3")["x-amz-meta-foo"] == "bar"

@testset "test coroutine reading" begin
    @sync begin
        for i in 1:2
            @async begin
                @test s3_get(bucket_name, "key3") == b"data3.v3"
                println("success ID: $i")
            end
        end
    end
end

# Check raw return of XML object...
xml = "<?xml version='1.0'?><Doc><Text>Hello</Text></Doc>"
s3_put(aws, bucket_name, "file.xml", xml, "text/xml")
@test String(s3_get(aws, bucket_name, "file.xml", raw=true)) == xml
@test s3_get(aws, bucket_name, "file.xml")["Text"] == "Hello"

# Check object copy function...

s3_copy(bucket_name, "key 1";
        to_bucket = bucket_name, to_path = "key 1.copy")

@test s3_get(aws, bucket_name, "key 1.copy") == b"data1.v1"

url = s3_sign_url(aws, bucket_name, "key 1")
curl_output = ""
@repeat 3 try
    curl_output = readstring(`curl -s -o - $url`)
catch e
    @delay_retry if true end
end
@test curl_output == "data1.v1"

fn = "/tmp/jl_qws_test_key1"
if isfile(fn)
    rm(fn)
end
@repeat 3 try
    s3_get_file(aws, bucket_name, "key 1", fn)
catch e
    sleep(1)
    @retry if true end
end
@test readstring(fn) == "data1.v1"
rm(fn)

# Check exists and list objects functions...

for key in ["key 1", "key2", "key3", "key 1.copy"]
    @test s3_exists(bucket_name, key)
    @test key in [o["Key"] for o in s3_list_objects(aws, bucket_name)]
end

# Check delete...

s3_delete(aws, bucket_name, "key 1.copy")

@test !("key 1.copy" in [o["Key"] for o in s3_list_objects(aws, bucket_name)])

# Check metadata...

meta = s3_get_meta(aws, bucket_name, "key 1")
@test meta["ETag"] == "\"68bc8898af64159b72f349b391a7ae35\""

# Check versioned object content...

versions = s3_list_versions(aws, bucket_name, "key3")
@test length(versions) == 3
@test (s3_get(aws, bucket_name, "key3"; version = versions[3]["VersionId"])
      == b"data3.v1")
@test (s3_get(aws, bucket_name, "key3"; version = versions[2]["VersionId"])
      == b"data3.v2")
@test (s3_get(aws, bucket_name, "key3"; version = versions[1]["VersionId"])
      == b"data3.v3")

@testset "default Content-Type" begin
# https://github.com/samoconnor/AWSS3.jl/issues/24

    ctype(key) = s3_get_meta(bucket_name, key)["Content-Type"]

    for k in [
        "file.foo",
        "file",
        "file_html",
        "file/html",
        "foobar.html/file.htm"]

        s3_put(aws, bucket_name, k, "x")
        @test ctype(k) == "application/octet-stream"
    end

    for (k, t) in [
        ("foo/bar/file.html",  "text/html"),
        ("x.y.z.js",           "application/javascript"),
        ("downalods/foo.pdf",  "application/pdf"),
        ("data/foo.csv",       "text/csv"),
        ("this.is.a.file.txt", "text/plain"),
        ("my.log",             "text/plain"),
        ("big.dat",            "application/octet-stream"),
        ("some.tar.gz",        "application/octet-stream"),
        ("data.bz2",           "application/octet-stream")]

        s3_put(aws, bucket_name, k, "x")
        @test ctype(k) == t
    end
end

# Check pruning of old versions...

s3_purge_versions(aws, bucket_name, "key3")
versions = s3_list_versions(aws, bucket_name, "key3")
@test length(versions) == 1
@test s3_get(aws, bucket_name, "key3") == b"data3.v3"

HTTP.ConnectionPool.showpool(STDOUT)

# Create objects...

max = 1000
sz = 10000
objs = [rand(UInt8(65):UInt8(75), sz) for i in 1:max]

asyncmap(x->AWSS3.s3(aws, "PUT", bucket_name;
                          path = "obj$(x[1])", content = x[2]),
         enumerate(objs);
         ntasks=30)
HTTP.ConnectionPool.showpool(STDOUT)

asyncmap(x->begin
    o = AWSS3.s3(aws, "GET", bucket_name; path = "obj$(x[1])")

    @test o == x[2]
end,
enumerate(objs);
ntasks=30)

HTTP.ConnectionPool.showpool(STDOUT)

#==============================================================================#
# End of file.
#==============================================================================#


text:
@testset "Domains" begin
  @testset "PointSet" begin
    # testing different ways to build a PointSet results in the same PointSet
    ps1 = PointSet([1. 0.; 0. 1.])            # as matrix
    ps2 = PointSet([(1.0, 0.0), (0.0, 1.0)])  # as vector of tuples

    @test ps1.coords == ps2.coords

    # building a PointSet of tuples with different lengths should fail
    @test_throws MethodError PointSet([(1.0, 0.0), (0.0, 1.0 , 4.0)])

    # more tests
    ps = PointSet([1. 0.; 0. 1.])

    @test nelms(ps) == 2
    @test coordinates(ps, 1) == [1., 0.]
    @test coordinates(ps, 2) == [0., 1.]

    @test sprint(show, ps) == "2 PointSet{Float64,2}"
    @test sprint(show, MIME"text/plain"(), ps) == "2 PointSet{Float64,2}\n 1.0  0.0\n 0.0  1.0"

    if visualtests
      Random.seed!(2019)
      @plottest plot(PointSet(rand(1,10))) joinpath(datadir,"pset1D.png") !istravis
      @plottest plot(PointSet(rand(2,10))) joinpath(datadir,"pset2D.png") !istravis
      @plottest plot(PointSet(rand(3,10))) joinpath(datadir,"pset3D.png") !istravis
      @plottest plot(PointSet(rand(1,10)),1:10) joinpath(datadir,"pset1D-data.png") !istravis
      @plottest plot(PointSet(rand(2,10)),1:10) joinpath(datadir,"pset2D-data.png") !istravis
      @plottest plot(PointSet(rand(3,10)),1:10) joinpath(datadir,"pset3D-data.png") !istravis
    end
  end

  @testset "RegularGrid" begin
    grid = RegularGrid{Float32}(200,100)
    @test ncoords(grid) == 2
    @test coordtype(grid) == Float32
    @test size(grid) == (200,100)
    @test nelms(grid) == 200*100
    @test coordinates(grid, 1) == [0.,0.]
    @test origin(grid) == [0f0, 0f0]
    @test spacing(grid) == [1f0, 1f0]

    grid = RegularGrid((200,100,50), (0.,0.,0.), (1.,1.,1.))
    @test ncoords(grid) == 3
    @test coordtype(grid) == Float64
    @test size(grid) == (200,100,50)
    @test nelms(grid) == 200*100*50
    @test coordinates(grid, 1) == [0.,0.,0.]
    @test origin(grid) == [0.,0.,0.]
    @test spacing(grid) == [1.,1.,1.]

    grid = RegularGrid((-1.,-1.), (1.,1.), dims=(200,100))
    @test ncoords(grid) == 2
    @test coordtype(grid) == Float64
    @test size(grid) == (200,100)
    @test nelms(grid) == 200*100
    @test coordinates(grid, 1) == [-1.,-1.]
    @test coordinates(grid, 200*100) == [1.,1.]
    @test origin(grid) == [-1.,-1.]

    grid = RegularGrid(100,200)
    @test sprint(show, grid) == "100×200 RegularGrid{Float64,2}"
    @test sprint(show, MIME"text/plain"(), grid) == "100×200 RegularGrid{Float64,2}\n  origin:  (0.0, 0.0)\n  spacing: (1.0, 1.0)"

    if visualtests
      @plottest plot(RegularGrid(10)) joinpath(datadir,"grid1D.png") !istravis
      @plottest plot(RegularGrid(10,20)) joinpath(datadir,"grid2D.png") !istravis
      @plottest plot(RegularGrid(10,20,30)) joinpath(datadir,"grid3D.png") !istravis
      @plottest plot(RegularGrid(10),[1,2,3,4,5,5,4,3,2,1]) joinpath(datadir,"grid1D-data.png") !istravis
      @plottest plot(RegularGrid(10,10),1:100) joinpath(datadir,"grid2D-data.png") !istravis
      # @plottest plot(RegularGrid(10,10,10),collect(1:1.0:1000)) joinpath(datadir,"grid3D-data.png") !istravis
    end
  end

  @testset "StructuredGrid" begin
    nx, ny, nz = 20, 10, 10

    X = [1.,5.,6.,10.]
    g1 = StructuredGrid(X)
    @test nelms(g1) == 4
    @test coordinates(g1) == [1. 5. 6. 10.]

    @test sprint(show, g1) == "4 StructuredGrid{Float64,1}"
    @test sprint(show, MIME"text/plain"(), g1) == "4 StructuredGrid{Float64,1}\n 1.0  5.0  6.0  10.0"

    X  = [x for x in range(0,10,length=nx), j in 1:ny]
    Y  = sin.(X) .+ [0.5j for i in 1:nx, j in 1:ny]
    g2 = StructuredGrid(X, Y)
    @test nelms(g2) == nx*ny
    @test size(g2) == (nx, ny)
    @test size(coordinates(g2)) == (2, nx*ny)

    @test sprint(show, g2) == "20×10 StructuredGrid{Float64,2}"

    X = [x for x in range(0,10,length=nx), j in 1:ny, k in 1:nz]
    Y = sin.(X) .+ [0.5j for i in 1:nx, j in 1:ny, k in 1:nz]
    Z = [1.0(k-1) for i in 1:nx, j in 1:ny, k in 1:nz]
    g3 = StructuredGrid(X, Y, Z)
    @test nelms(g3) == nx*ny*nz
    @test size(g3) == (nx, ny, nz)
    @test size(coordinates(g3)) == (3, nx*ny*nz)

    @test sprint(show, g3) == "20×10×10 StructuredGrid{Float64,3}"

    if visualtests
      @plottest plot(g1) joinpath(datadir,"sgrid1D.png") !istravis
      @plottest plot(g2) joinpath(datadir,"sgrid2D.png") !istravis
      @plottest plot(g3,camera=(30,60)) joinpath(datadir,"sgrid3D.png") !istravis
      @plottest plot(g1,[1.,2.,2.,1]) joinpath(datadir,"sgrid1D-data.png") !istravis
      @plottest plot(g2,1:nx*ny) joinpath(datadir,"sgrid2D-data.png") !istravis
      @plottest plot(g3,1:nx*ny*nz) joinpath(datadir,"sgrid3D-data.png") !istravis
    end
  end
end


text:
#todo:
#   write decorators? for original function (or anything which can make a
#       corresponding call in Julia)

using BitSAD

#done
#args are m, u, and v respectively in C version
function mod2sparse_mulvec(H, received_codeword, synd)
    M = mod2sparse_rows(H)
    N = mod2sparse_cols(H)

    for 1:M
        synd[i] = 0
    end

    for j in 1:N
        if received_codeword[j] == 1
            for e in mod2sparse_cols[j]
                synd[e.row] ^= 1
            end
        end
    end
end

function mod2sparse_first_in_row(m, i)
    return m.rows[i].right
end
function mod2sparse_first_in_col(m, j)
    return m.cols[j].down
end
function mod2sparse_last_in_row(m, i)
    return m.rows[i].left
end
function mod2sparse_last_in_col(m, j)
    return m.cols[j].up
end

function mod2sparse_next_in_row(e)
    return e.right
end
function mod2sparse_next_in_col(e)
    return e.down
end
function mod2sparse_prev_in_row(e)
    return e.left
end
function mod2sparse_prev_in_col(e)
    return e.up
end

function bp_decode_cy()
    bp_decode_log_prob_ratios()

function bp_decode_log_prob_ratios(dec)
#variables that are typically held in the bp_decoder class:
# ch[] (the channel probs)
# H (??)
# max_iter
# synd[] (??)
# log_prob_ratios[]
# bp_decoding[]
# bp_decoding_synd
# N
# M

# for each of the above: am I reading or read/writing? Where is it initialized,
#   and what to? Determine if it needs to be defined before bp_decode_log_prob_ratios
#   (ie as part of the decoder struct) or if it can be defined in this function.
# update: above does not matter, any values which are written to need to be returned
#   (i.e. pass by reference) and any which are read are defined outside. So all in struct.

    for j in 1:dec.N #fix syntax
        e = mod2sparse_first_in_col(dec.H, j)
        while !(mod2sparse_at_end(e))
            e.bit_to_check = log((1 - dec.channel_probs[j])/(dec.channel_probs[j]))
            e = mod2sparse_next_in_col(e)

    dec.converge = 0
    for iteration in 1:(dec.max_iter+1)
        #iter = iteration # probably redundant/unused
        #if bp_method==2 #this is the setting for our trial run
        #product sum check_to_bit messages
        for i in 1:dec.M
            e = mod2sparse_first_in_row(dec.H, i)
            temp = 1.0
            while !(mod2sparse_at_end(e))
                e.check_to_bit = temp
                temp *= tanh(e.bit_to_check/2)
                e = mod2sparse_next_in_row(e)

            e = mod2sparse_last_in_row(dec.H, i)
            temp = 1.0
            while !(mod2sparse_at_end(e))
                e.check_to_bit *= temp
                e.check_to_bit = (((-1)^dec.synd[i]) * log((1 + e.check_to_bit) / (1 - e.check_to_bit)))
                temp *= tanh(e.bit_to_check/2)
                e = mod2sparse_prev_in_row(e)

        # bit-to-check messages
        for j in 1:dec.N
            e = mod2sparse_first_in_col(dec.H, j)
            temp = log((1-dec.channel_probs[j]) / (dec.channel_probs[j]))

            while !(mod2sparse_at_end(e)):
                e.bit_to_check = temp
                temp += e.check_to_bit
                e = mod2sparse_next_in_col(e)

            dec.log_prob_ratios[j] = temp
            if temp <= 0
                dec.bp_decoding[j]=1
            else
                dec.bp_decoding[j]=0

            e = mod2sparse_last_in_col(dec.H, j)
            temp = 0.0
            while !(mod2sparse_at_end)
                e.bit_to_check += temp
                temp += e.check_to_bit
                e = mod2sparse_prev_in_col(e)

        #all of these params need to be defined still, as do N and M for the message loops.
        mod2sparse_mulvec(dec.H, dec.bp_decoding, dec.bp_decoding_synd)

        equal = 1
        for check in 1:dec.M
            if dec.synd[check] != dec.bp_decoding_synd[check]
                equal = 0
                break
        if equal == 1
            dec.converge = 1 #what is this flag used for?
            return 1

        return 0

function mod2sparse_at_end(e)
    return (e.row < 0)

mutable struct mod2entry()
    bit_to_check
    check_to_bit
    row
    col
    sgn

    left
    right
    up
    down
end

mutable struct mod2sparse()
    n_rows::int
    n_cols::int
    rows
    cols
end

mutable struct decoder()
    channel_probs
    H
    max_iter
    synd
    log_prob_ratios #these are defined in bp_decode_log_prob_ratios
    bp_decoding
    bp_decoding_synd
    N
    M
    converge
end

# complicated stuff:
# 2-step process:
#   1. Run the modified decoder with sample parameters. Don't care about result,
#       but as it runs it will populate a julia file with equivalent julia code
#       (including calls to mod2sparse stuff and bp_decode_log_prob_ratios)
#   2. Call generatehw on the populated julia file, write the results to a verilog
#       file.

function programTrace()
    # make original call
    #

function main()
#add the following lines to L
    f_verilog, f_circuit = generatehw(programTrace)
    io = open("hwfile.vl", "w")
    write(io, f_verilog)
    close(io)

    jl_file = open("jl_code.jl", 'a')
    jl_file.writelines(L) #where L is an array of all the lines of julia to be written
    jl_file.close()

max_stars_repo_name: alexandrupaler/ldpc

max_issues_repo_name: alexandrupaler/ldpc

max_forks_repo_name: alexandrupaler/ldpc


text:
#!/usr/bin/env julia

# Use the EnsembleSampler
# This means that each likelihood evaluation is *serial*, while the walkers are parallelized.

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    # "--opt1"
    # help = "an option with an argument"
    "--p", "-p"
    help = "number of processes to add"
    arg_type = Int
    default = 0
    "--run_index", "-r"
    help = "Output run index"
    arg_type = Int
    "--chain"
    help = "Write the chain to ~/web/ directory?"
    action = :store_true
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

# Since we've made this a #!/usr/bin/env julia script, we can no longer specify the extra
# processors via the julia prompt, so we need to ghost this behavior back in.
if parsed_args["p"] > 0
    addprocs(parsed_args["p"])
end

import YAML
config = YAML.load(open(parsed_args["config"]))

outfmt(run_index::Int) = config["out_base"] * @sprintf("run%02d/", run_index)

# This code is necessary for multiple simultaneous runs on odyssey
# so that different runs do not write into the same output directory
if parsed_args["run_index"] == nothing
    run_index = 0
    outdir = outfmt(run_index)
    while ispath(outdir)
        println(outdir, " exists")
        run_index += 1
        outdir = outfmt(run_index)
    end
else
    run_index = parsed_args["run_index"]
    outdir = outfmt(run_index)
    println("Deleting old $outdir")
    run(`rm -rf $outdir`)
end

# make the output directory
println("Creating ", outdir)
mkdir(outdir)

using DiskJockey.constants
using DiskJockey.visibilities
using DiskJockey.image
using DiskJockey.gridding
using DiskJockey.model
using Base.Test

# load data and figure out how many channels
dvarr = DataVis(config["data_file"])
nchan = length(dvarr)

if haskey(config, "exclude")
    exclude = config["exclude"]
    # which channels of the dset to fit
    keylist = filter(x->(!in(x, exclude)), Int[i for i=1:nchan])
else
    keylist = Int[i for i=1:nchan]
end

# go through any previously created directories and remove them
function cleardirs!(keylist::Vector{Int})
    println("Removing old directories")
    for key in keylist
        keydir = basedir * "jud$key"
        run(`rm -rf $keydir`)
    end
    println("Removed directories")
end

# make the values of run_index and config available on all processes
global run_id=run_index
global cfg=config
global kl=keylist

# Now, redo this to only load the dvarr for the keys that we need, and conjugate
dvarr = DataVis(cfg["data_file"], kl)
visibilities.conj!(dvarr)
nchan = length(dvarr)

const global species = cfg["species"]
basefmt(id::Int) = cfg["base_dir"] * @sprintf("run%02d/", id)
const global basedir = basefmt(run_id)

const global npix = cfg["npix"] # number of pixels, can alternatively specify x and y separately

# Keep track of the current home directory
const global homedir = pwd() * "/"

# make the internal Judith directory, if it doesn't exist
if !ispath(basedir)
    println("Creating ", basedir)
    mkdir(basedir)
end

# Clear all directories
cleardirs!(keylist)

# Delete the old log file (if it exists)
const logfile = outdir * "log.log"
if isfile(logfile)
    rm(logfile)
end

using Logging
# change the default logger
Logging.configure(filename=logfile, level=DEBUG)

debug("Created logfile.")

# Create the model grid
grd = cfg["grid"]
global const grid = Grid(grd["nr"], grd["ntheta"], grd["r_in"], grd["r_out"], true)

# Regenerate all of the static files (e.g., amr_grid.inp)
# so that they may be later copied
debug("Writing grid")
write_grid(basedir, grid)
debug("Wrote grid")

# Here, if we actually are going to be fixing distance, evaluate one image to get the interpolation
# closures

if cfg["fix_d"]
    # @everywhere pixsize = cfg["pixsize"] # [cm]
    # @everywhere pix_AU = pixsize/AU # [AU]

    # Simply calculate pix_AU as 1.1 * (2 * r_out) / npix
    # This is assuming that RADMC always calculates the image the same
    @everywhere pix_AU = (1.1 * 2 * grd["r_out"]) / cfg["npix"] # [AU/pixel]

    @everywhere dl = sin(pix_AU/cfg["parameters"]["dpc"][1] * arcsec)

    @everywhere uu = fftshift(fftfreq(npix, dl)) * 1e-3 # [kλ]
    @everywhere vv = fftshift(fftfreq(npix, dl)) * 1e-3 # [kλ]

    # For each channel, also calculate the interpolation closures
    @everywhere int_arr = Array(Function, nchan)
    @everywhere for (i, dset) in enumerate(dvarr)
        int_arr[i] = plan_interpolate(dset, uu, vv)
    end
end

# This function is fed to the EnsembleSampler
# That means, using the currently available global processes, like the data visibilities,
# and it must create it's own temporary directory to write the necessary files for
# RADMC to run.
function fprob(p::Vector{Float64})

    # println("Starting in ", pwd())

    # Each walker needs to create it's own temporary directory
    # where all RADMC files will reside and be driven from
    # It only needs to last for the duration of this function, so let's use a tempdir
    keydir = mktempdir() * "/"

    # Copy all relevant configuration scripts to this subdirectory
    # these are mainly setup files that will be static throughout the run
    # they were written by DJInitialize.jl and write_grid()
    for fname in ["radmc3d.inp", "wavelength_micron.inp", "lines.inp", "molecule_" * molnames[species] * ".inp"]
        ff = homedir * fname
        run(`cp $ff $keydir`)
    end
    # run(`cp radmc3d.inp $keydir`)
    # run(`cp wavelength_micron.inp $keydir`)
    # run(`cp lines.inp $keydir`)
    #
    ag = basedir * "amr_grid.inp"
    run(`cp $ag $keydir`)

    # mf = "molecule_" * molnames[species] * ".inp"
    # run(`cp $mf $keydir`)

    # change the subprocess to reside in this directory for the remainder of the run
    # where it will drive its own independent RADMC3D process for a subset of channels
    cd(keydir)
    # println("Changed to ", pwd())

    # Fix the following arguments: gamma, dpc
    gamma = 1.0 # surface temperature gradient exponent

    if cfg["fix_d"]
        dpc = cfg["parameters"]["dpc"][1] # [pc] distance
        M_star, r_c, T_10, q, logM_gas, ksi, incl, PA, vel, mu_RA, mu_DEC = p
    else
        M_star, r_c, T_10, q, logM_gas, ksi, dpc, incl, PA, vel, mu_RA, mu_DEC = p
    end

    # Enforce hard priors on physical parameters
    # Short circuit evaluation if we know the RADMC won't be valid.
    if ksi <= 0. || T_10 <= 0. || r_c <= 0.0 || M_star <= 0.0 || T_10 > 1500. || q < 0. || q > 1.0
        return -Inf
    end

    if incl < 0. || incl > 180.
        return -Inf
    end

    M_gas = 10^logM_gas

    # If we are going to fit with some parameters dropped out, here's the place to do it
    pars = Parameters(M_star, r_c, T_10, q, gamma, M_gas, ksi, dpc, incl, PA, vel, mu_RA, mu_DEC)

    # Compute parameter file using model.jl, write to disk in current directory
    write_model(pars, keydir, grid, species)

    # copied from f
    # dvarr, int_arr = data

    # We are using the following conventions: inclination ranges from
    # 0 to 180 degrees. 0 means face on, angular momentum vector pointing
    # at observer; 90 means edge on; and 180 means face on, angular momentum
    # vector pointing away from observer.
    # These are the same as the RADMC-3D conventions.
    incl = pars.incl # [deg]

    # We also adopt the RADMC-3D convention for position angle, which defines position angle
    # by the angular momentum vector.
    # A positive PA angle means the disk angular momentum vector will be
    # rotated counter clockwise (from North towards East).
    PA = pars.PA # [deg]

    vel = pars.vel # [km/s]

    # Doppler shift the dataset wavelengths to rest-frame wavelength
    beta = vel/c_kms # relativistic Doppler formula
    lams = Array(Float64, nchan)
    for i=1:nchan
        lams[i] =  dvarr[i].lam * sqrt((1. - beta) / (1. + beta)) # [microns]
    end

    write_lambda(lams, keydir) # write into current directory

    # println("Now in $keydir")
    # run(`ls`)

    # println("p ", p)
    # Run RADMC-3D, redirect output to /dev/null
    tic()
    run(`radmc3d image incl $incl posang $PA npix $npix loadlambda`)
    println("RADMC3D call")
    toc()

    # Read the RADMC-3D images from disk (we should already be in sub-directory)
    im = imread()

    if cfg["fix_d"]
        # After the fact, we should be able to check that the pixel size of the image is the
        # same as the one we originally calculated from the outer disk radius.
        @test_approx_eq_eps im.pixsize_x/AU pix_AU 1e-4
    end

    # Convert raw images to the appropriate distance
    skim = imToSky(im, pars.dpc)

    # Apply the gridding correction function before doing the FFT
    # No shift needed, since we will shift the resampled visibilities
    corrfun!(skim)

    tic()
    lnprobs = Array(Float64, nchan)
    # Do the Fourier domain stuff per channel
    for i=1:nchan
        dv = dvarr[i]
        # FFT the appropriate image channel
        vis_fft = transform(skim, i)

        if cfg["fix_d"]
            # Interpolate the `vis_fft` to the same locations as the DataSet
            mvis = int_arr[i](dv, vis_fft)
        else
            mvis = ModelVis(dv, vis_fft)
        end

        # Apply the phase shift here
        phase_shift!(mvis, pars.mu_RA, pars.mu_DEC)

        # Calculate chi^2 between these two
        lnprobs[i] = lnprob(dv, mvis)
    end

    println("Interpolation time")
    toc()
    # remove the temporary directory in which we currently reside
    # println("removing $keydir")
    run(`rm -rf $keydir`)

    # Sum them all together and feed back to the master process
    lnp = sum(lnprobs)

    return lnp

end

debug("Initializing MCMC")
using Distributions
using PDMats

pp = config["parameters"]
# The parameters we'll be using
if config["fix_d"]
    params = ["M_star", "r_c", "T_10", "q", "logM_gas", "ksi", "incl", "PA", "vel", "mu_RA", "mu_DEC"]
else
    params = ["M_star", "r_c", "T_10", "q", "logM_gas", "ksi", "dpc", "incl", "PA", "vel", "mu_RA", "mu_DEC"]
end

nparam = length(params)
starting_param = Array(Float64, nparam)
jumps = Array(Float64, nparam)

for i=1:nparam
    starting_param[i], jumps[i] = pp[params[i]]
end

# given these made up parameters, call fprob
println(fprob(starting_param))

# println("Now for the profiling.")
#
# clear_malloc_data()
#
# fprob(starting_param)
#
# @profile fprob(starting_param)
# Profile.print()

tic()
fprob(starting_param)
println("Full f Timing")
toc()

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
abstract type Consumer end

mutable struct Residential <: Consumer
    periods::Float64
    econsumption::Float64 #kWh-month
    growth::Float64
    peak_power::Union{Nothing,Array{Float64}}
    load_curve::Array{Float64}
    tariff::Tariff
    decision_timeframe::Int64
    rate_return::Float64
end

mutable struct CommIndus <: Consumer
    periods::Float64
    econsumption::Float64 #kWh-month
    growth::Float64
    peak_power::Union{Nothing,Array{Float64}}
    load_curve::Array{Float64}
    tariff::Tariff
    decision_timeframe::Int64
    rate_return::Float64
end

mutable struct TMT <: Consumer
    periods::Float64
    econsumption::Float64 #kWh-month
    growth::Float64
    peak_power::Union{Nothing,Array{Float64}}
    load_curve::Array{Float64}
    tariff::Tariff
    decision_timeframe::Int64
    rate_return::Float64
end

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
"""
getIndRand(cc)

Get the row of the indicator for column `cc`, else return a random row.

# Examples

```
getIndRand(cc)
```
"""
function getIndRand(cc)
    if all(cc .== 0)
        println("zeros encountered")
        return rand(1:length(cc))
    else
        return findfirst(x->x==1,cc)
    end
end

# """
# invEpInd(epochinds,Z)

# Invert the partitioned vector of indicator matrices `Z` and return an array similar to `getA()`.

# # Examples

# ```
# invEpInd(epochinds,Z)
# ```
# """
# function invEpInd(epochinds,Z)
#     inds = reduce(vcat,[reduce(vcat,epi) for epi in epochinds])
#     a = reduce(hcat,convert.(Array,Z))[:,sortperm(inds)]
#     map(t->getIndRand(t), eachcol(a))
# end

"""
    initShared(corpus,partition,k,l,alphalocal,betalocal)

Initialize:
1) shared data vectors `T`,`D`
2) shared assignment vector `Z`
3) shared indicator matrices `Z`,`J`,`W`
4) shared sufficient statistics `nj`
5) global sufficient statistics `nwglobal`,`nglobal`
6) shared increment matrices `nwinc`,`ninc`
7) vectors of 'local indices': `indt`,`inda`,`indd`,`indnj`,`indnw` corresponding to each worker id in `wkids`

# Examples

```
using Distributions
corp = map(x->Word.(rand(1:100,x),Categorical([0.5,0.3,0.2])),[10,50,100])
dist = [[1,2],[3]]
initShared(corpus,partition,3,100,fill(0.1,3),fill(0.1,100))
```
"""
function initShared(corpus,partition,k,l,alphalocal,betalocal;statefloat=Float32,stateint=Int32,sstatfloat=Float32,sstatint=Int32)
    Z = SharedArray(ones(Bool,k,sum(map(x->length(x),corpus)))) # size Z is k by # corpus tokens
    J = SharedArray(ones(Bool,sum(map(x->length(x),corpus)),length(partition))) # size J is # corpus tokens by # docs
    W = SharedArray(ones(Bool,sum(map(x->length(x),corpus)),l)) # size W is # corpus tokens by # terms
    T = SharedVector(ones(stateint,sum(map(x->length(x),corpus))))
    A = SharedVector(ones(stateint,sum(map(x->length(x),corpus))))
    D = SharedVector(ones(stateint,sum(map(x->length(x),corpus))))
    ends = cumsum(map(x->length(x),corpus)) # ends of the segments for each doc
    starts = cumsum(map(x->length(x),corpus)) .- map(x->length(x),corpus) .+ 1 # the starts of the segments for each doc
    nwstarts = map(x->x-l,cumsum(fill(l,length(partition)))).+1
    nwends = cumsum(fill(l,length(partition)))
    nj = SharedArray(zeros(sstatint,k,length(partition)))
    nw = SharedArray(zeros(sstatint,k,l,length(partition)))
    nwinc = SharedArray(zeros(sstatint,k,l*length(partition)))
    indnw = [(1:k,nwstarts[1]:nwends[1])]
    indt = [starts[1]:ends[1]]
    for i in 2:length(partition)
        push!(indnw,(1:k,nwstarts[i]:nwends[i]))
        push!(indt,starts[i]:ends[i])
    end
    indt=SharedArray(indt)
    indnw=SharedArray(indnw)

    @sync @distributed for i in 1:length(partition)
        ba = getBitarrays(corpus,partition[i],ttype=stateint,atype=stateint,dtype=stateint)
        T[starts[i]:ends[i]] .= vec(ba[:T])
        A[starts[i]:ends[i]] .= vec(ba[:A])
        D[starts[i]:ends[i]] .= vec(ba[:D])
        Z[1:k,starts[i]:ends[i]] = getZ(ba[:A],k)
        J[starts[i]:ends[i],1:length(partition)] = getJ(corpus[partition[i]],partition[i],length(partition))
        W[starts[i]:ends[i],1:l] = getW(corpus[partition[i]],l)
        nw[:,:,i] .= Z[1:k,starts[i]:ends[i]] * W[starts[i]:ends[i],1:l]
        nj[1:k,i] .= vec(Z[1:k,starts[i]:ends[i]]*J[starts[i]:ends[i],1:length(partition)][:,partition[i]])
    end
    nw = SharedArray{sstatint,2}(sum(nw,dims=3)[:,:,1])

    println("initialization complete for:")
    println("Z{",typeof(Z),"}: ",Base.summarysize(Z)," bytes")
    println("J{",typeof(J),"}: ",Base.summarysize(J)," bytes")
    println("W{",typeof(W),"}: ",Base.summarysize(W)," bytes")
    println("T{",typeof(T),"}: ",Base.summarysize(T)," bytes")
    println("A{",typeof(A),"}: ",Base.summarysize(A)," bytes")
    println("D{",typeof(D),"}: ",Base.summarysize(D)," bytes")
    println("nwinc{",typeof(nwinc),"}: ",Base.summarysize(nwinc)," bytes")
    # copy the global counts to separate sampling vectors
    if length(partition) > 1
        for p in 2:length(partition)
            nw = hcat(nw,nw[:,1:l])
        end
        println("allocated nw{",typeof(nw),"}: ",sizeof(nw)," bytes")
    end

    return (T=T,A=A,D=D,Z=Z,J=J,W=W,
            nj=nj,nwinc=nwinc,nw=nw,
            alpha=SharedArray{sstatfloat,1}(alphalocal),beta=SharedArray{sstatfloat,2}(betalocal),
            indt=indt,indnw=indnw)
end

# """
#   getwids(arr)

# Get the worker ids of a DArray or SharedArray.

# # Examples
# ```
# getwids(arr)
# ```
# """
# function getwids(arr)
#     if typeof(arr) <: DArray
#         ids = [(i,@fetchfrom i DistributedArrays.localindices(arr)) for i in workers()]
#     elseif typeof(arr) <: SharedArray
#         ids = [(i,@fetchfrom i SharedArrays.localindices(arr)) for i in workers()]
#     end
#     ids = filter(!isnothing,[length(i[2][1]) > 0 ? i[1] : nothing for i in ids])
# end

"""
  splitApprox(corpus,n)

Re-order `corpus` into `n` chunks each containing approximately the same number of documents.
Return the shuffled view `c_view`, a vector `partition` of unit-ranges for each chunk, and vector `partmap`, where `c_view[i] = corpus[partmap[i]]`.

# Examples
```
splitApprox(corpus,n)
```
"""
function splitApprox(corpus,n)
    parts = [findall(x->x==i,mod.([1:length(corpus);],n) .+ 1) for i in 1:n]
    partmap = reduce(vcat,parts)
    ind1 = [0,cumsum([length(i) for i in parts])[1:end-1]...] .+ 1
    ind2 = cumsum([length(i) for i in parts])
    partition = UnitRange.(ind1,ind2)
    c_perm = corpus[partmap]
    c_perm,partition,partmap
end

"""
  epochIndex(corpus,n,l)

Make random, approximately equal-length partition of docs and terms.

# Examples
```
epochIndex(corpus,n,l)
```
"""
function epochIndex(corpus,n,l)
    # partition the documents
    docindex = [findall(x->x==i,mod.([1:length(corpus);],n) .+ 1) for i in 1:n]
    # partition the terms
    termindex = [findall(x->x==i,mod.([1:l;],n) .+ 1) for i in 1:n]
    return (docindex=docindex,termindex=termindex)
end

"""
  filterTokens(n,D,T,docindex,termindex)

Generate global token epoch indices by filtering concatenated tokens.

# Examples
```
filterTokens(n,D,T,docindex,termindex)
```
"""
function filterTokens(n,D,T,docindex,termindex)
    epochs = []
    for i in 1:n
        termindexShft = circshift(termindex,i-1)
        chunks = []
        for j in 1:n
            docs = reduce(vcat,[findall(d->d==di,vec(D)) for di in docindex[j]])
            terms = reduce(vcat,[findall(t->t==ti,vec(T)) for ti in termindexShft[j]])
            inds = intersect(docs,terms)
            push!(chunks,inds)
        end
        push!(epochs,chunks)
    end
    epochs
end


text:
function GeographicZoneDetector(arg0::BodyShape, arg1::SphericalPolygonsSet, arg2::jdouble)
    return GeographicZoneDetector((BodyShape, SphericalPolygonsSet, jdouble), arg0, arg1, arg2)
end

function GeographicZoneDetector(arg0::jdouble, arg1::jdouble, arg2::BodyShape, arg3::SphericalPolygonsSet, arg4::jdouble)
    return GeographicZoneDetector((jdouble, jdouble, BodyShape, SphericalPolygonsSet, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function g(obj::GeographicZoneDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_body(obj::GeographicZoneDetector)
    return jcall(obj, "getBody", BodyShape, ())
end

function get_margin(obj::GeographicZoneDetector)
    return jcall(obj, "getMargin", jdouble, ())
end

function get_zone(obj::GeographicZoneDetector)
    return jcall(obj, "getZone", SphericalPolygonsSet, ())
end

function with_margin(obj::GeographicZoneDetector, arg0::jdouble)
    return jcall(obj, "withMargin", GeographicZoneDetector, (jdouble,), arg0)
end