text:
module GasDispersion

export Scenario, plume, puff, scenario_builder

using Interpolations

# helpful utilities
include("utils/scenario.jl")
include("utils/scenario_builder.jl")
include("utils/utils.jl")

# plume models
include("models/gaussian_plume.jl")
include("models/britter_mcquaid_plume.jl")

# puff models
include("models/gaussian_puff.jl")
include("models/britter_mcquaid_puff.jl")

"""
    plume(scenario::Scenario; model=:gaussian, kwargs...)

Runs the plume dispersion model on the given scenario and returns a function
giving the concentration of the form
    c(x, y, z[, t])

If `model` is unspecified, defaults to gaussian, `kwargs` are passed to the
plume model. All model parameters are assumed to be in SI base units (i.e.
distances in m, velocities in m/s, mass in kg, etc.)
"""
function plume(scenario::Scenario; model=:gaussian, kwargs...)
    if model==:gaussian
        return gaussian_plume_factory(scenario; kwargs...)
    elseif model==:brittermcquaid
        return britter_plume_factory(scenario; kwargs...)
    else
        error_string = "plume dispersion model $model is not currently implemented"
        error(error_string)
    end
end

"""
    puff(scenario::Scenario; model=:gaussian, kwargs...)

Runs the puff dispersion model on the given scenario and returns a function
giving the concentration of the form
    c(x, y, z, t)

If `model` is unspecified, defaults to gaussian, `kwargs` are passed to the
puff model. All model parameters are assumed to be in SI base units (i.e.
distances in m, velocities in m/s, mass in kg, etc.)
"""
function puff(scenario::Scenario; model=:gaussian, kwargs...)
    if model==:gaussian
        return gaussian_puff_factory(scenario; kwargs...)
    # elseif model=="britter-mcquaid"
    #     return britter_plume_factory(scenario; kwargs...)
    else
        error_string = "puff dispersion model $model is not currently implemented"
        error(error_string)
    end
end

end


text:
function energy_balance(ps_m::CanonicalModel,
                        time_range::UnitRange{Int64},
                        initial_conditions::Array{Tuple{String,Float64},1},
                        p_eff_data::Array{Tuple{String,Float64},1},
                        cons_name::Symbol,
                        var_names::Tuple{Symbol,Symbol,Symbol})

    ps_m.constraints[cons_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, [i[1] for i in initial_conditions], time_range)
    # println(initial_conditions)
    for (ix,i) in enumerate(initial_conditions)

        ps_m.constraints[cons_name][i[1], 1] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_names[3]][i[1], 1] == i[2] + (ps_m.variables[var_names[1]][i[1], 1])/p_eff_data[ix][2] - (ps_m.variables[var_names[2]][i[1], 1])*p_eff_data[ix][2])

    end

    for t in time_range[2:end], (ix,i) in enumerate(initial_conditions)

        ps_m.constraints[cons_name][i[1], t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_names[3]][i[1], t] == ps_m.variables[var_names[3]][i[1], t-1] + (ps_m.variables[var_names[1]][i[1], t])/p_eff_data[ix][2]  - (ps_m.variables[var_names[2]][i[1], t])*p_eff_data[ix][2])

    end

    return

end

ext: jl