module Particles
export AbstractParticles

    abstract type AbstractParticles end
    using Juriba.Constants
    # const J = Juriba.Constants
    # using GSL
    # const GSL_CONST = GSL.Const

    # const cee  =           GSL_CONST.CGSM_SPEED_OF_LIGHT
    # const emgm=            GSL_CONST.CGSM_MASS_ELECTRON
    # const pmgm  =          GSL_CONST.CGSM_MASS_PROTON
    # const kboltz    =      GSL_CONST.CGSM_BOLTZMANN
    # const kboltz_kev2erg = 1.6022e-9
    # const gr_to_kev   =    5.6095883571872e+29
    # const me_kev    =      511.
    # const emerg   =        (GSL_CONST.CGSM_MASS_ELECTRON*GSL_CONST.CGSM_SPEED_OF_LIGHT^2.0)
    # const charg   =        4.8e-10
    # const sigtom =         GSL_CONST.CGSM_THOMSON_CROSS_SECTION
    # const herg   =         GSL_CONST.CGSM_PLANCKS_CONSTANT_H
    # const hkev   =         (GSL_CONST.CGSM_PLANCKS_CONSTANT_H*6.2415e8)
    # const mjy   =          1.e-26
    # const re0    =         2.81794e-13
    # const gconst  =        GSL_CONST.CGSM_GRAVITATIONAL_CONSTANT
    # const sbconst  =       GSL_CONST.CGSM_STEFAN_BOLTZMANN_CONSTANT
    # const aconst   =       7.56e-15

    struct pl_params 
        s::Float64 
        n::Float64
    end

    struct bkn_params 
        s1::Float64
        s2::Float64
        brk::Float64
        max::Float64 
        m::Float64
    end

    struct th_params 
        t::Float64 
        n::Float64
        m::Float64 
    end

    struct k_params 
        t::Float64
        k::Float64 
    end

    struct injection_mixed_params
        s::Float64 
        t::Float64 
        nth::Float64
        npl::Float64
        m::Float64 
        min::Float64
        max::Float64
        cutoff::Float64
    end

    struct injection_kappa_params
        t::Float64 
        k::Float64
        n::Float64 
        m::Float64
    end

    struct injection_pl_params
        s::Float64
        n::Float64 
        m::Float64 
        max::Float64
    end

    struct injection_bkn_params
        s1::Float64
        s2::Float64
        brk::Float64
        max::Float64 
        m::Float64
        n::Float64
    end

    function count_particles(p::AbstractParticles)
        temp = 0.0
        for i ∈ 1:p.size-1
            temp = temp + 0.5*(p.p[i+1] - p.p[i])*(p.ndens[i+1] + p.ndens[i])
        end
        return temp
    end

    function count_particles_energy(p::AbstractParticles)
        temp = 0.0
        for i ∈ 1:p.size-1
            temp = temp + 0.5*(p.p[i+1] - p.p[i])*(p.gdens[i+1] + p.gdens[i])
        end
        return temp
    end

    function av_p(p::AbstractParticles)
        temp = 0.0
        for i ∈ 1:p.size-1
            temp = temp + 0.5*(p.p[i+1] - p.p[i])*
                            (p.p[i+1]*p.ndens[i+1] + p.p[i]*p.ndens[i])
        end
        return temp/count_particles(p)
    end

    function av_gamma(p::AbstractParticles)
        return sqrt((av_p(p)/(p.mass_gr*cee))^2 + 1)
    end

    function av_psq(p::AbstractParticles)
        temp = 0.0
        for i ∈ 1:p.size - 1
            temp = temp + 0.5*(p.p[i+1] - p.p[i])*p.p[i+1]^2*
                        p.ndens[i+1]*p[i]^2*p.ndens[i]
        end

        return temp/count_particles(p)
    end

    function av_gammasq(p::AbstractParticles)

        return sqrt(av_psq(p)/((p.mass_gr*cee)^2) + 1)

    end

    function initialize_gdens(p::AbstractParticles)
        for i ∈ 1:p.size
            p.gdens[i] = p.ndens[i]*p.gamma[i]*p.mass_gr*cee/sqrt(p.gamma[i]^2 - 1)
        end
    end

    function initialize_pdens(p::AbstractParticles)
        for i ∈ 1:p.size
            ndens[i] = p.gdens[i]*p.p[i]/((p.mass_gr*cee)^2*
                    sqrt(p.p[i]/(p.mass_gr*cee)^2 + 1))
        end
    end

    function gdens_differentiate(p::AbstractParticles)
        temp = Vector{Float64}(undef, p.size)

        for i in eachindex(temp)
            temp[i] = p.gdens[i]/p.gamma[i]
        end

        for i ∈ 1:p.size - 1
            p.gdens_diff[i] = (temp[i+1] - temp[i])/
                            (p.mass_gr*cee^2*p.gamma[i+1] - p.gamma[i])
        end

        p.gdens_diff[end] = p.gdens_diff[end-1]
    end

    function set_mass(p::AbstractParticles, m)
        p.mass_gr = m
        p.mass_kev = m*gr_to_kev
    end

end