text:
"""
    OLS{TF<:AbstractFloat} <: RegressionModel

Data from an ordinary least squares regression.
"""
struct OLS{TF<:AbstractFloat} <: RegressionModel
    X::Matrix{TF}
    invXX::Matrix{TF}
    coef::VecOrMat{TF}
    resid::VecOrMat{TF}
    score::Matrix{TF}
end

function OLS(Y::AbstractVecOrMat, X::AbstractMatrix)
    X = convert(Matrix, X)
    crossx = cholesky!(X'X)
    coef = X'Y
    ldiv!(crossx, coef)
    invXX = inv!(crossx)
    resid = X * coef
    resid .= Y .- resid
    score = getscore(X, resid)
    return OLS(X, invXX, coef, resid, score)
end

modelmatrix(m::OLS) = m.X
coef(m::OLS) = m.coef
residuals(m::OLS) = m.resid

show(io::IO, ::OLS) = print(io, "OLS regression")

"""
    AbstractEstimator

Supertype for all estimators.
"""
abstract type AbstractEstimator end

@fieldequal AbstractEstimator

show(io::IO, e::AbstractEstimator) = print(io, typeof(e).name.name)

"""
    LeastSquaresLP <: AbstractEstimator

Ordinary least squares estimator for local projecitons.

# Reference
Jordà, Òscar. 2005. "Estimation and Inference of Impulse Responses by Local Projections."
American Economic Review 95 (1): 161-182.
"""
struct LeastSquaresLP <: AbstractEstimator end

show(io::IO, ::MIME"text/plain", ::LeastSquaresLP) =
    print(io, "Ordinary least squares local projection")

"""
    AbstractEstimatorResult

Supertype for all estimator-specific additional results.
"""
abstract type AbstractEstimatorResult end

show(io::IO, er::AbstractEstimatorResult) = print(io, typeof(er).name.name)

"""
    LeastSquaresLPResult{TF<:AbstractFloat} <: AbstractEstimatorResult

Additional results from estimating least squares local projections.
See also [`LeastSquaresLP`](@ref) and [`LocalProjectionResult`](@ref).

# Field
- `ms::Vector{OLS{TF}}`: data from the OLS regressions.
"""
struct LeastSquaresLPResult{TF<:AbstractFloat} <: AbstractEstimatorResult
    ms::Vector{OLS{TF}}
end

"""
    VarName

A type union of all accepted types for indexing variables by name.
"""
const VarName = Union{Symbol,TransformedVar{Symbol}}

"""
    LocalProjectionResult

Results from local projection estimation.
"""
struct LocalProjectionResult{TE<:AbstractEstimator,
        ER<:Union{AbstractEstimatorResult,Nothing},
        VCE<:CovarianceEstimator, TF<:AbstractFloat} <: StatisticalModel
    B::Array{TF,3}
    V::Array{TF,3}
    T::Vector{Int}
    est::TE
    estres::ER
    vce::VCE
    ynames::Vector{VarName}
    xnames::Vector{VarName}
    wnames::Vector{VarName}
    lookupy::Dict{VarName,Int}
    lookupx::Dict{VarName,Int}
    lookupw::Dict{VarName,Int}
    nlag::Int
    minhorz::Int
    subset::Union{BitVector,Nothing}
    normnames::Union{Vector{VarName},Nothing}
    normtars::Union{Vector{VarName},Nothing}
    normmults::Union{Vector{TF},Nothing}
    endonames::Union{Vector{VarName},Nothing}
    ivnames::Union{Vector{VarName},Nothing}
    firststagebyhorz::Bool
    nocons::Bool
end

"""
    coef(r::LocalProjectionResult, horz::Int, xwname::VarName; kwargs...)

Return the coefficient estimate at horizon `horz` for variable with column name `xwname`.

# Keywords
- `yname::VarName=r.ynames[1]`: name of the outcome variable.
- `lag::Int=0`: lag of variable `xwname`; being 0 means the variable is contemporaneous.
"""
function coef(r::LocalProjectionResult, horz::Int, xwname::VarName;
        yname::VarName=r.ynames[1], lag::Int=0)
    if lag > 0
        k = length(r.xnames)+r.lookupw[xwname]+(lag-1)*length(r.wnames)
        return r.B[k, r.lookupy[yname], horz]
    else
        return r.B[r.lookupx[xwname], r.lookupy[yname], horz]
    end
end

"""
    vcov(r::LocalProjectionResult, horz::Int, xwname1::VarName; kwargs...)

Return the variance of the estimate at horizon `horz` for variable with column name `xwname1`.
If `xwname2` is specified, return the covariance between `xwname1` and `xwname2`.

# Keywords
- `yname1::VarName=r.ynames[1]`: name of the outcome variable for `xwname1`.
- `lag1::Int=0`: lag of variable `xwname1`; being 0 means the variable is contemporaneous.
- `xwname2::VarName=xwname1`: the second variable involved in the covariance.
- `yname2::VarName=yname1`: name of the outcome variable for `xwname1`.
- `lag2::Int=lag1`: lag of variable `xwname2`; being 0 means the variable is contemporaneous.
"""
function vcov(r::LocalProjectionResult, horz::Int, xwname1::VarName;
        yname1::VarName=r.ynames[1], lag1::Int=0, xwname2::VarName=xwname1,
        yname2::VarName=yname1, lag2::Int=lag1)
    if lag1 > 0
        k1 = length(r.xnames)+r.lookupw[xwname1]+(lag1-1)*length(r.wnames)
    else
        k1 = r.lookupx[xwname1]
    end
    if lag2 > 0
        k2 = length(r.xnames)+r.lookupw[xwname2]+(lag2-1)*length(r.wnames)
    else
        k2 = r.lookupx[xwname2]
    end
    n1 = r.lookupy[yname1]
    n2 = r.lookupy[yname2]
    return r.V[(k1-1)*length(r.ynames)+n1, (k2-1)*length(r.ynames)+n2, horz]
end

"""
    VarIndex

A type union of all accepted types for indexing variables.
"""
const VarIndex = Union{Integer, Symbol, TransformedVar{<:Union{Integer,Symbol}}}

"""
    VarIndexPair

`Pair` type with both type parameters being a subtype of [`VarIndex`](@ref).
"""
const VarIndexPair = Pair{<:VarIndex,<:VarIndex}

_checknames(names) = all(n isa VarIndex for n in names)

_toint(data, name::Symbol) = Tables.columnindex(data, name)
_toint(data, i::Integer) = Int(i)

_toname(data, name::Symbol) = name
_toname(data, i::Integer) = Tables.columnnames(data)[i]

function _getcols(data, ynames, xnames, wnames, states, addylag, nocons; TF=Float64)
    # The names must be iterable
    ynames isa VarIndex && (ynames = (ynames,))
    xnames isa VarIndex && (xnames = (xnames,))
    wnames isa VarIndex && (wnames = (wnames,))
    states isa VarIndex && (states = (states,))
    argmsg = ", must contain either integers, `Symbol`s or `TransformedVar`s"
    _checknames(ynames) || throw(ArgumentError("invalid ynames"*argmsg))
    length(xnames)==0 || _checknames(xnames) || throw(ArgumentError("invalid xnames"*argmsg))
    length(wnames)==0 || _checknames(wnames) || throw(ArgumentError("invalid wnames"*argmsg))
    states===nothing || _checknames(states) || throw(ArgumentError("invalid states"*argmsg))

    # Convert all column indices to Int
    ynames = VarIndex[_toint(data, n) for n in ynames]
    xnames = VarIndex[_toint(data, n) for n in xnames]
    wnames = VarIndex[_toint(data, n) for n in wnames]
    addylag && union!(wnames, ynames)

    ys = Any[getcolumn(data, n) for n in ynames]
    xs = Any[getcolumn(data, n) for n in xnames]
    ws = Any[getcolumn(data, n) for n in wnames]
    sts = states === nothing ? nothing : Any[getcolumn(data, n) for n in states]
    if !nocons
        push!(xs, ones(TF, length(ys[1])))
        xnames = (xnames..., :constant)
    end
    return ynames, xnames, wnames, states, ys, xs, ws, sts
end

function _makeYX(ys, xs, ws, sts, nlag::Int, horz::Int, subset::Union{BitVector,Nothing};
        firststage::Bool=false, TF=Float64)
    ny = length(ys)
    ny > 0 || throw(ArgumentError("ys cannot be empty"))
    Tfull = size(ys[1],1)
    nlag > 0 || throw(ArgumentError("nlag must be at least 1"))
    # Number of rows possibly involved in estimation
    T = Tfull - nlag - horz
    nx = length(xs)
    nw = length(ws)
    nw > 0 || throw(ArgumentError("ws cannot be empty"))
    nX = sts === nothing ? nx+nw*nlag : nx+nw*nlag*length(sts)
    T > nX || throw(ArgumentError(
        "not enough observations for nlag=$nlag and horz=$horz"))
    # Indicators for valid observations within the T rows
    esampleT = trues(T)
    # A cache for indicators
    aux = BitVector(undef, T)
    # Construct matrices
    Y = Matrix{TF}(undef, T, ny)
    for j in 1:ny
        if firststage
            v = view(vec(ys[j], subset, :x, horz, TF), nlag+1:Tfull-horz)
            subset === nothing || j > 1 || (esampleT .&= view(subset, nlag+1:Tfull-horz))
        else
            v = view(vec(ys[j], subset, :y, horz, TF), nlag+horz+1:Tfull)
            subset === nothing || j > 1 || (esampleT .&= view(subset, nlag+horz+1:Tfull))
        end
        _esample!(esampleT, aux, v)
        copyto!(view(Y,esampleT,j), view(v,esampleT))
    end
    X = Matrix{TF}(undef, T, nX)
    if nx > 0
        for j in 1:nx
            v = view(vec(xs[j], subset, :x, horz, TF), nlag+1:Tfull-horz)
            subset === nothing || j > 1 || (esampleT .&= view(subset, nlag+1:Tfull-horz))
            _esample!(esampleT, aux, v)
            copyto!(view(X,esampleT,j), view(v,esampleT))
        end
    end
    if sts === nothing
        for j in 1:nw
            for l in 1:nlag
                v = view(vec(ws[j], subset, :x, horz, TF), nlag+1-l:Tfull-horz-l)
                subset === nothing || j > 1 ||
                    (esampleT .&= view(subset, nlag+1-l:Tfull-horz-l))
                _esample!(esampleT, aux, v)
                # Variables with the same lag are put together
                copyto!(view(X,esampleT,nx+(l-1)*nw+j), view(v,esampleT))
            end
        end
    else
        ns = length(sts)
        ssts = Vector{AbstractVector}(undef, ns)
        for j in 1:ns
            v = view(vec(sts[j], subset, :x, horz, TF), nlag+1:Tfull-horz)
            _esample!(esampleT, aux, v)
            ssts[j] = v
        end
        for j in 1:nw
            for l in 1:nlag
                for s in 1:ns
                    v = view(vec(ws[j], subset, :x, horz, TF), nlag+1-l:Tfull-horz-l)
                    subset === nothing || j > 1 ||
                        (esampleT .&= view(subset, nlag+1-l:Tfull-horz-l))
                    s > 1 || _esample!(esampleT, aux, v)
                    X[esampleT, nx+(l-1)*nw*ns+(j-1)*ns+s] .=
                        view(ssts[s],esampleT) .* view(v,esampleT)
                end
            end
        end
    end
    T1 = sum(esampleT)
    if T1 < T
        T1 > size(X,2) || throw(ArgumentError(
            "not enough observations for nlag=$nlag and horz=$(horz)"))
        T = T1
        Y = Y[esampleT, :]
        X = X[esampleT, :]
    end
    return Y, X, T, esampleT
end

function _lp(ys, xs, ws, sts, nlag::Int, horz::Int, vce::Union{CovarianceEstimator,Nothing},
        subset::Union{BitVector,Nothing}; TF=Float64)
    Y, X, T, esampleT = _makeYX(ys, xs, ws, sts, nlag, horz, subset; TF=TF)
    m = OLS(Y, X)
    return coef(m), vcov(m, vce), T, m
end

function _normalize!(data, normalize, xnames, xs, ws, sts, nlag, minhorz, subset; TF=Float64)
    snames = normalize isa Pair ? (normalize[1],) : (p[1] for p in normalize)
    ix_s = Vector{Int}(undef, length(snames))
    for (i, n) in enumerate(snames)
        fr = findfirst(x->x==_toint(data, n), xnames)
        fr === nothing && throw(ArgumentError(
            "Regressor $(n) is not found among x variables"))
        ix_s[i] = fr
    end
    nnames = normalize isa Pair ? (normalize[2],) : (p[2] for p in normalize)
    yn = Any[getcolumn(data, n) for n in nnames]
    Bn, Vn, Tn, _ = _lp(yn, xs, ws, sts, nlag, minhorz, nothing, subset; TF=TF)
    normmults = TF[Bn[ix,i] for (i,ix) in enumerate(ix_s)]
    xs[ix_s] .*= normmults
    normnames = VarName[_toname(data, n) for n in snames]
    normtars = VarName[_toname(data, n) for n in nnames]
    return normnames, normtars, normmults
end

_normalize!(data, normalize::Nothing, xnames, xs, ws, sts, nlag, minhorz, subset; TF=Float64) =
    (nothing, nothing, nothing)

function _firststage(ys, xs, ws, sts, nlag::Int, horz::Int, subset::Union{BitVector,Nothing};
        TF=Float64)
    Y, X, T, esampleT = _makeYX(ys, xs, ws, sts, nlag, horz, subset; firststage=true, TF=TF)
    # Get first-stage estimates
    bf = X'Y
    ldiv!(cholesky!(X'X), bf)
    # Make the fitted values, need to have the same length as original data
    Tfull = size(ys[1],1)
    nY = length(ys)
    fitted = ((fill(convert(TF, NaN), Tfull) for _ in 1:nY)...,)
    Xb = X * bf
    for i in 1:nY
        copyto!(view(view(fitted[i], nlag+1:Tfull-horz), esampleT), view(Xb,:,i))
    end
    return fitted
end

function _iv!(data, iv, firststagebyhorz, xnames, xs, ws, sts, nlag, minhorz, subset;
        TF=Float64)
    endonames = iv[1]
    endonames isa VarIndex && (endonames = (endonames,))
    ix_iv = Vector{Int}(undef, length(endonames))
    for (i, n) in enumerate(endonames)
        fr = findfirst(x->x==_toint(data, n), xnames)
        fr === nothing && throw(ArgumentError(
            "Endogenous variable $(n) is not found among x variables"))
        ix_iv[i] = fr
    end
    ivnames = iv[2]
    ivnames isa VarIndex && (ivnames = (ivnames,))
    length(ivnames)==0 && throw(ArgumentError("invalid specification of option iv"))
    argmsg = ", must contain either integers, `Symbol`s or `TransformedVar`s"
    _checknames(ivnames) || throw(ArgumentError("invalid ivnames"*argmsg))
    # Collect columns used for first-stage regression
    yfs = xs[ix_iv]
    xfs = Any[(xs[i] for i in 1:length(xs) if !(i in ix_iv))...,
        (getcolumn(data, n) for n in ivnames)...]
    if !firststagebyhorz
        any(x->x isa Cum, endonames) &&
            @warn "firststagebyhorz=false while endogenous variables contain Cum"
        # Replace the endogenous variable with the fitted values
        xs[ix_iv] .= _firststage(yfs, xfs, ws, sts, nlag, minhorz, subset; TF=TF)
    end
    endonames = VarName[_toname(data, n) for n in endonames]
    ivnames = VarName[_toname(data, n) for n in ivnames]
    return endonames, ivnames, ix_iv, yfs, xfs
end

_iv!(data, iv::Nothing, firststagebyhorz, xnames, xs, ws, sts, nlag, minhorz, subset;
    TF=Float64) = (nothing, nothing, nothing, nothing, nothing)

function _est(::LeastSquaresLP, data, xnames, ys, xs, ws, sts,
        nlag, minhorz, nhorz, vce, subset,
        iv, ix_iv, yfs, xfs, firststagebyhorz; TF=Float64)
    ny = length(ys)
    nstate =sts === nothing ? 1 : length(sts)
    nr = length(xs) + length(ws)*nlag*nstate
    B = Array{TF,3}(undef, nr, ny, nhorz)
    V = Array{TF,3}(undef, nr*ny, nr*ny, nhorz)
    T = Vector{Int}(undef, nhorz)
    M = Vector{OLS{TF}}(undef, nhorz)
    for h in minhorz:minhorz+nhorz-1
        if iv !== nothing && firststagebyhorz
            xs[ix_iv] .= _firststage(yfs, xfs, ws, sts, nlag, h, subset; TF=TF)
        end
        i = h - minhorz + 1
        Bh, Vh, T[i], M[i] = _lp(ys, xs, ws, sts, nlag, h, vce, subset; TF=TF)
        B[:,:,i] = reshape(Bh, nr, ny, 1)
        V[:,:,i] = reshape(Vh, nr*ny, nr*ny, 1)
    end
    return B, V, T, LeastSquaresLPResult(M)
end

"""
    lp([estimator], data, ynames; kwargs...)

Estimate local projections with the specified `estimator`
for outcome variable(s) with column name(s) `ynames` in `data` table.
If the `estimator` is not specified, [`LeastSquaresLP`](@ref) is assumed.
The input `data` must be `Tables.jl`-compatible.

# Keywords
- `xnames=()`: indices of contemporaneous regressors from `data`.
- `wnames=()`: indices of lagged control variables from `data`.
- `nlag::Int=4`: number of lags to be included for the lagged control variables.
- `nhorz::Int=1`: total number of horizons to be estimated.
- `minhorz::Int=0`: the minimum horizon involved in estimation.
- `normalize::Union{VarIndexPair,Vector{VarIndexPair},Nothing}=nothing`: normalize the magnitude of the specified regressor(s) based on their initial impact on the respective targeted outcome variable.
- `iv::Union{Pair,Nothing}=nothing`: endogenous varable(s) paired with instrument(s).
- `firststagebyhorz::Bool=false`: estimate the first-stage regression separately for each horizon.
- `states=nothing`: variable(s) representing the states for a state-dependent specification.
- `vce::CovarianceEstimator=HARVCE(EWC())`: the variance-covariance estimator.
- `subset::Union{BitVector,Nothing}=nothing`: subset of `data` to be used for estimation.
- `addylag::Bool=true`: include lags of the outcome variable(s).
- `nocons::Bool=false`: do not add the constant term.
- `TF::Type=Float64`: numeric type used for estimation.
"""
function lp(estimator, data, ynames;
        xnames=(), wnames=(), nlag::Int=4, nhorz::Int=1, minhorz::Int=0,
        normalize::Union{VarIndexPair,Vector{VarIndexPair},Nothing}=nothing,
        iv::Union{Pair,Nothing}=nothing, firststagebyhorz::Bool=false, states=nothing,
        vce::CovarianceEstimator=HARVCE(EWC()),
        subset::Union{BitVector,Nothing}=nothing,
        addylag::Bool=true, nocons::Bool=false, TF::Type=Float64)

    checktable(data)
    ynames, xnames, wnames, states, ys, xs, ws, sts =
        _getcols(data, ynames, xnames, wnames, states, addylag, nocons)

    if any(x->x isa Cum, ynames)
        addylag && @warn "addylag=true while outcome variables contain Cum"
    end

    normalize !== nothing && iv !== nothing &&
        throw(ArgumentError("options normalize and iv cannot be specified at the same time"))
    normnames, normtars, normmults =
        _normalize!(data, normalize, xnames, xs, ws, sts, nlag, minhorz, subset, TF=TF)
    endonames, ivnames, ix_iv, yfs, xfs =
        _iv!(data, iv, firststagebyhorz, xnames, xs, ws, sts, nlag, minhorz, subset, TF=TF)

    B, V, T, er = _est(estimator, data, xnames, ys, xs, ws, sts,
        nlag, minhorz, nhorz, vce, subset, iv, ix_iv, yfs, xfs, firststagebyhorz, TF=TF)

    ynames = VarName[_toname(data, i) for i in ynames]
    xnames = VarName[_toname(data, i) for i in xnames]
    wnames = VarName[_toname(data, i) for i in wnames]
    return LocalProjectionResult(B, V, T, estimator, er, vce, ynames, xnames, wnames,
        Dict{VarName,Int}(n=>i for (i,n) in enumerate(ynames)),
        Dict{VarName,Int}(n=>i for (i,n) in enumerate(xnames)),
        Dict{VarName,Int}(n=>i for (i,n) in enumerate(wnames)),
        nlag, minhorz, subset, normnames, normtars, normmults,
        endonames, ivnames, firststagebyhorz, nocons)
end

lp(data, ynames; kwargs...) = lp(LeastSquaresLP(), data, ynames; kwargs...)

"""
    lp(r::LocalProjectionResult, vce::CovarianceEstimator)

Reestimate the variance-covariance matrices with `vce`.
"""
function lp(r::LocalProjectionResult{LeastSquaresLP}, vce::CovarianceEstimator)
    V = similar(r.V)
    K = size(V, 1)
    for h in 1:length(r.T)
        V[:,:,h] = reshape(vcov(r.estres.ms[h], vce), K, K, 1)
    end
    return LocalProjectionResult(r.B, V, r.T, r.est, r.estres, vce,
        r.ynames, r.xnames, r.wnames, r.lookupy, r.lookupx, r.lookupw,
        r.nlag, r.minhorz, r.subset, r.normnames, r.normtars, r.normmults,
        r.endonames, r.ivnames, r.firststagebyhorz, r.nocons)
end

show(io::IO, r::LocalProjectionResult) = print(io, typeof(r).name.name)

_vartitle(::LocalProjectionResult) = "Variable Specifications"

function _varinfo(r::LocalProjectionResult, halfwidth::Int)
    namex = "Regressor"*(length(r.xnames) > 1 ? "s" : "")
    xs = join(r.xnames, " ")
    namex = length(namex) + length(xs) > halfwidth ? (namex,) : namex
    namew = "Lagged control"*(length(r.wnames) > 1 ? "s" : "")
    ws = join(r.wnames, " ")
    namew = length(namew) + length(ws) > halfwidth ? (namew,) : namew
    info = Pair[
        "Outcome variable"*(length(r.ynames) > 1 ? "s" : "") => join(r.ynames, " "),
        "Minimum horizon" => r.minhorz,
        namex => xs,
        namew => ws]
    if r.endonames !== nothing
        push!(info, "Endogenous variable"*(length(r.endonames) > 1 ? "s" : "") =>
            join(r.endonames, " "),
            "Instrument"*(length(r.ivnames) > 1 ? "s" : "") =>
            join(r.ivnames, " "))
    end
    return info
end

_estimatortitle(::LocalProjectionResult) = nothing
_estimatorinfo(::LocalProjectionResult, ::Int) = nothing

function show(io::IO, ::MIME"text/plain", r::LocalProjectionResult;
        totalwidth::Int=78, interwidth::Int=4+mod(totalwidth,2))
    H = length(r.T)
    print(io, "$(typeof(r).name.name) with $(r.nlag) lag")
    print(io, r.nlag > 1 ? "s " : " ", "over $H horizon")
    println(io, H > 1 ? "s:" : ":")
    halfwidth = div(totalwidth-interwidth, 2)
    titles = (_vartitle(r), _estimatortitle(r))
    blocks = (_varinfo(r, halfwidth), _estimatorinfo(r, halfwidth))
    for ib in 1:length(blocks)
        if blocks[ib] !== nothing
            println(io, repeat('─', totalwidth))
            if titles[ib] !== nothing
                println(io, titles[ib])
                println(io, repeat('─', totalwidth))
            end
            # Use count to determine whether a newline is needed
            count = 0
            for e in blocks[ib]
                # Allow the flexibility of printing on the entire row
                if e[1] isa Tuple
                    isodd(count) && print(io, '\n')
                    print(io, e[1][1], ':')
                    println(io, lpad(e[2], totalwidth - length(e[1][1]) - 1))
                    count = 0
                else
                    # Allow having vacancy by specifying an empty string
                    if e[1] != ""
                        print(io, e[1], ':')
                        print(io, lpad(e[2], halfwidth - length(e[1]) - 1))
                    end
                    count += 1
                    print(io, isodd(count) ? repeat(' ', interwidth) : '\n')
                end
            end
            isodd(count) && print(io, '\n')
        end
    end
    print(io, repeat('─', totalwidth))
end


text:
export orbitdiagram

"""
    orbitdiagram(ds::DiscreteDynamicalSystem, i, p_index, pvalues; kwargs...)
Compute the orbit diagram (also called bifurcation diagram) of the given system,
saving the `i` variable(s) for parameter values `pvalues`. The `p_index` specifies
which parameter of the equations of motion is to be changed.

`i` can be `Int` or `AbstractVector{Int}`.
If `i` is `Int`, returns a vector of vectors. Else
it returns vectors of vectors of vectors.
Each entry are the points at each parameter value.

## Keyword Arguments
* `Ttr::Int = 1000` : Transient steps;
  each orbit is evolved for `Ttr` first before saving output.
* `n::Int = 100` : Amount of points to save for each initial condition.
* `dt = 1` : Stepping time. Changing this will give you the orbit diagram of
  the `dt` order map.
* `u0 = get_state(ds)` : Initial condition. Besides a vector you can also give
  a vector of vectors such that `length(u0) == length(pvalues)`. Then each parameter
  has a different initial condition.
* `ulims = (-Inf, Inf)` : only record system states within `ulims`
  (only valid if `i isa Int`).

See also [`poincaresos`](@ref) and [`produce_orbitdiagram`](@ref).
"""
function orbitdiagram(
        ds::DDS{IIP, S, D}, idxs, p_index, pvalues;
        n::Int = 100, Ttr::Int = 1000, u0 = get_state(ds), dt = 1, ulims = nothing
    ) where {IIP, S, D}

    p0 = ds.p[p_index]
    if D == 1
        idxs != 1 &&
        error("You have a 1D system and yet you gave `i=$i`. What's up with that!?")
    end

    typeof(u0) <: Vector{<:AbstractVector} && @assert length(u0)==length(pvalues)
    i = idxs isa Int ? idxs : SVector{length(idxs), Int}(idxs...)
    !isnothing(ulims) && i isa SVector && error("If `i` is a vector, you can't use `ulims`.")

    output = _initialize_od_output(ds.u0, i, n, length(pvalues))
    integ = integrator(ds)
    _fill_orbitdiagram!(output, integ, i, pvalues, p_index, n, Ttr, u0, dt, ulims)
    ds.p[p_index] = p0
    return output
end

function _initialize_od_output(u::S, i::Int, n, l) where {S}
    output = [zeros(eltype(S), n) for k in 1:l]
end
function _initialize_od_output(u::S, i::SVector, n, l) where {S}
    s = u[i]
    output = [Vector{typeof(s)}(undef, n) for k in 1:l]
end

function _fill_orbitdiagram!(output, integ, i, pvalues, p_index, n, Ttr, u0, dt, ulims)
    isavector = i isa AbstractVector
    for (j, p) in enumerate(pvalues)
        integ.p[p_index] = p
        st = (u0 isa Vector{ <: _innertype(integ.u) }) ? u0[j] : u0
        reinit!(integ, st)
        step!(integ, Ttr)
        if isavector || isnothing(ulims) # if-clause gets compiled away (I hope)
            @inbounds for k in 1:n
                step!(integ, dt)
                output[j][k] = integ.u[i]
            end
        else
            k = 1
            while k ≤ n
                step!(integ, dt)
                u = integ.u[i]
                @inbounds if ulims[1] ≤ u ≤ ulims[2]
                    output[j][k] = u
                    k += 1
                end
            end
        end
    end
    return output
end

_innertype(::Real) = Real
_innertype(::AbstractVector) = AbstractVector

ext: jl


text:
using NeighbourLists
using DelimitedFiles: readdlm
using JuLIP: r_sum
using LinearAlgebra: rmul!

export EAM

# =================== General Single-Species EAM Potential ====================

"""
`struct EAM1`

Generic Single-species EAM potential, to specify it, one needs to
specify the pair potential ϕ, the electron density ρ and the embedding
function F.

The most convenient constructor is likely via tabulated values,
more below.

# Constructors:
```
EAM(pair::PairPotential, eden::PairPotential, embed)
EAM(fpair::AbstractString, feden::AbstractString, fembed::AbstractString; kwargs...)
```

## Constructing an EAM potential from tabulated values

At the moment only the .plt format is implemented. Files can e.g. be
obtained from
* [NIST](https://www.ctcms.nist.gov/potentials/)

Use the `EAM(fpair, feden, fembed)` constructure. The keyword arguments specify
details of how the tabulated values are fitted; see
`?SplinePairPotential` for more details.

TODO: implement other file formats.
"""
struct EAM1{T1, T2, T3, T4} <: SitePotential
   ϕ::T1    # pair potential
   ρ::T2    # electron density potential
   F::T3    # embedding function
   info::T4
end

@pot EAM1

EAM1(ϕ, ρ, F) = EAM1(ϕ, ρ, F, nothing)

cutoff(V::EAM1) = max(cutoff(V.ϕ), cutoff(V.ρ))

function evaluate!(tmp, V::EAM1, R::AbstractVector{JVec{T}}) where {T}
   if  length(R) == 0
      return V.F(T(0.0))
   else
      return V.F(sum(V.ρ ∘ norm, R)) + T(0.5) * sum(V.ϕ ∘ norm, R)
   end
end

function evaluate_d!(dEs::AbstractVector{JVec{T}}, tmp, V::EAM1, Rs) where {T}
   if length(Rs) == 0; return dEs; end
   ρ̄ = sum(V.ρ ∘ norm, Rs)
   dF = @D V.F(ρ̄)
   for (i, R) in enumerate(Rs)
      r = norm(R)
      dEs[i] = ((T(0.5) * (@D V.ϕ(r)) + dF * (@D V.ρ(r))) / r) * R
   end
   return dEs
end

alloc_temp_dd(V::EAM1, N::Integer) =
      ( ∇ρ = zeros(JVecF, N),
         r = zeros(Float64, N) )

evaluate_dd!(hEs, tmp, V::EAM1, R) = _hess_!(hEs, tmp, V, R, identity)

# ff preconditioner specification for EAM potentials
#   (just replace id with abs and hess with precon in the hessian code)
precon!(hEs, tmp, V::EAM1, R, stab=0.0) = _hess_!(hEs, tmp, V, R, abs, stab)

function _hess_!(hEs, tmp, V::EAM1, R::AbstractVector{JVec{T}}, fabs, stab=T(0)
                 ) where {T}
   for i = 1:length(R)
      r = norm(R[i])
      tmp.r[i] = r
      tmp.∇ρ[i] = grad(V.ρ, r, R[i])
   end
   # precompute some stuff
   ρ̄ = sum(V.ρ, tmp.r)
   dF = @D V.F(ρ̄)
   ddF = @DD V.F(ρ̄)
   # assemble
   for i = 1:length(R)
      for j = 1:length(R)
         hEs[i,j] = (1-stab) * fabs(ddF) * tmp.∇ρ[i] * tmp.∇ρ[j]'
      end
      r = tmp.r[i]
      S = R[i] / r
      dϕ = @D V.ϕ(r)
      dρ = @D V.ρ(r)
      ddϕ = @DD V.ϕ(r)
      ddρ = @DD V.ρ(r)
      a = fabs(0.5 * ddϕ + dF * ddρ)
      b = fabs((0.5 * dϕ + dF *  dρ) / r)
      hEs[i,i] += ( (1-stab) * ( (a-b) * S * S' + b * I )
                   + stab  * ( (a+b) * I ) )
   end
   return hEs
end

# implementation of EAM models using data files

function EAM(fpair::AbstractString, feden::AbstractString,
             fembed::AbstractString; kwargs...)
   pair = SplinePairPotential(fpair; kwargs...)
   eden = SplinePairPotential(feden; kwargs...)
   embed = SplinePairPotential(feden; fixcutoff = false, kwargs...)
   return EAM1(pair, eden, embed)
end

#
# Load EAM file from .fs file format
#
function EAM(fname::AbstractString; kwargs...)

   if fname[end-3:end] == ".eam"
      error(".eam is not yet implemented, please file an issue")
   elseif fname[end-6:end] == ".eam.fs"
      error(".eam.fs is not yet implemented, please file an issue")
   elseif fname[end-2:end] == ".fs"
      return eam_from_fs(fname; kwargs...)
   end

   error("unknwon EAM file format, please file an issue")
end

# ================= Finnis-Sinclair Potential =======================

mutable struct FSEmbed end
@pot FSEmbed
evaluate(V::FSEmbed, ρ̄) = - sqrt(ρ̄)
evaluate_d(V::FSEmbed, ρ̄::T) where {T} = - T(0.5) / sqrt(ρ̄)
evaluate_dd(V::FSEmbed, ρ̄::T) where {T} = T(0.25) * ρ̄^(T(-3/2))

"""
`FinnisSinclair`: constructs an EAM potential with embedding function
-√ρ̄.
"""
FinnisSinclair(pair::PairPotential, eden::PairPotential) =
   EAM1(pair, eden, FSEmbed())

function FinnisSinclair(fpair::AbstractString, feden::AbstractString; kwargs...)
   pair = SplinePairPotential(fpair; kwargs...)
   eden = SplinePairPotential(feden; kwargs...)
   return FinnisSinclair(pair, eden)
end

# ================= Various File Loaders =======================

"""
`eam_from_fs(fname; kwargs...) -> EAM`

Read a `.fs` file specifying and EAM / Finnis-Sinclair potential.
"""
function eam_from_fs(fname; kwargs...)
   F, ρfun, ϕfun, ρ, r, info = read_fs(fname)
   return EAM1( SplinePairPotential(r, ϕfun; kwargs...) * (@analytic r -> 1/r),
               SplinePairPotential(r, ρfun; kwargs...),
               SplinePairPotential(ρ, F; fixcutoff= false, kwargs...),
               info )
end

"""
`read_fs(fname)` -> F, ρfun, ϕfun, ρ, r, info

Read a `.fs` file specifying and EAM / Finnis-Sinclair potential.
The description of the file format is taken from
   http://lammps.sandia.gov/doc/pair_eam.html
see also
   https://sites.google.com/a/ncsu.edu/cjobrien/tutorials-and-guides/eam
"""
function read_fs(fname)
   f = open(fname)
   # ignore the first three lines
   for n = 1:3
      readline(f)
   end

   # line 4: Nelements Element1 Element2 ... ElementN
   L4 = readline(f) |> chomp |> IOBuffer |> readdlm
   if L4[1] != 1
      error("""`read_fs`: the file `$fname` is for alloys, but only
               the single species potential is implemented so far.""")
   end
   info = Dict()
   info[:species] = L4[2]

   # line 5: Nrho, drho, Nr, dr, cutoff
   L5 = readline(f) |> IOBuffer |> readdlm
   Nrho, drho, Nr, dr, cutoff = Int(L5[1]), L5[2], Int(L5[3]), L5[4], L5[5]
   info[:cutoff] = cutoff

   # line 6: atomic number, mass, lattice constant, lattice type (e.g. FCC)
   L6 = readline(f) |> IOBuffer |> readdlm
   info[:number], info[:mass], info[:a0], info[:lattice] = tuple(L6...)

   # all the data
   data = readdlm(f)
   @assert length(data) == Nrho+2*Nr
   ρ = range(0, stop=(Nrho-1)*drho, length=Nrho)
   r = range(cutoff - (Nr-1)*dr, stop=cutoff, length=Nr)

   # embedding function
   F = data[1:Nrho]
   # density function
   ρfun  = data[Nrho+1:Nrho+Nr]
   # interatomic potential
   ϕfun  = data[Nrho+Nr+1:Nrho+2*Nr]

   return F, ρfun, ϕfun, ρ, r, info
end


text:
using Weave
using Test

function weavestring(source; doctype = "pandoc", informat="markdown", mod=Main)
    p1 = Weave.parse_doc(source, informat)
    doc = Weave.WeaveDoc("dummy1.jmd", p1, Dict())
    return Weave.run(doc, doctype=doctype, mod=mod)
end

smod = """
```julia
module TestMod
    x = 3

    function printx()
        print("x")
    end

    export printx
end
```

```julia
using .TestMod
printx()
```
"""

wdoc = weavestring(smod)
@test wdoc.chunks[1].output == "Main.TestMod\n"
@test wdoc.chunks[2].output == "x"

sdoc = weavestring(smod, mod=:sandbox)
@test occursin(r"Main.WeaveSandBox[0-9]*.TestMod\n", sdoc.chunks[1].output)
@test sdoc.chunks[2].output == "x"

include_weave(Main, joinpath(@__DIR__, "documents/include_test.jmd"))
@test Test1.x == 10
@test Test1.y == 20
@test Test1.testing("weave") == "weave"
@test Test2.d["a"] == "α"
@test split(Test2.doc, r"\r\n|\n")[1] == "---"


text:
function ikegami(m_prefs::Vector{Vector{Int}}, f_prefs::Vector{Vector{Int}})
    m = length(m_prefs)
    n = length(f_prefs)

    m_prefs_2d = Matrix{Int64}(n+1, m)
    n_prefs_2d = Matrix{Int64}(m+1, n)

    for (t,i) in enumerate(m_prefs)
        if length(i) != n
            new_element = Array(Int, n+1-length(i))
            new_element[1] = 0
            new_element[2:end] = [j for j in range(1,n) if !(j in i)]
            m_prefs[t] = append!(i, new_element)
        else
            m_prefs[t] = push!(i, 0)
        end
        m_prefs_2d[:, t] = m_prefs[t]
    end

    for (t,i) in enumerate(n_prefs)
        if length(i) != m
            new_element = Array(Int, m+1-length(i))
            new_element[1] = 0
            new_element[2:end] = [j for j in range(1,m) if !(j in i)]
            n_prefs[t] = append!(i, new_element)
        else
            n_prefs[t] = push!(i, 0)
        end
        n_prefs_2d[:, t] = n_prefs[t]
    end

    return ikegamida(m_prefs_2d, n_prefs_2d)
end

function ikegamida(m_prefs, f_prefs)

    m_num, f_num = size(m_prefs)[2], size(f_prefs)[2]
    m_pool, f_pool = collect(1:m_num), collect(1:f_num)
    match_or_unmatch = Array(Bool, m_num)
    for i in 1:m_num
        match_or_unmatch[i] = true
    end
    m_matched, f_matched = ones(Int64, m_num) + f_num, ones(Int64, f_num) + m_num

    m_rank = Array(Int64, (f_num + 1, m_num))
    f_rank = Array(Int64, (m_num + 1, f_num))

    sorting = Array(Int64, (f_num + 1, 2))
    sorting[:, 1] = collect(1:(f_num + 1))

    for i in 1:m_num
        sorting[:, 2] = m_prefs[:, i]
        m_rank[:,i] = sortrows(sorting, by = x->(x[2]))[:,1]
        sorting = sortrows(sorting, by = x->(x[1]))
    end

    sorting2 = Array(Int64, (m_num + 1, 2))
    sorting2[:, 1] = collect(1:(m_num + 1))

    for i in 1:f_num
        sorting2[:, 2] = f_prefs[:, i]
        f_rank[:,i] = sortrows(sorting2, by = x->(x[2]))[:,1]
        sorting2 = sortrows(sorting2, by = x->(x[1]))
    end

    while sum(match_or_unmatch) > 0
        for i in 1:m_num
            if match_or_unmatch[i] == true

                for j in m_prefs[:, i]

                    if j == 0
                        m_matched[i] = 0
                        match_or_unmatch[i] = false
                        break

                    else
                        no_marriage = f_rank[1, j]

                        if f_rank[i+1, j] < no_marriage

                            if f_matched[j] == m_num + 1
                                f_matched[j] = i
                                m_matched[i] = j
                                match_or_unmatch[i] = false
                                break

                            else

                                if f_rank[f_matched[j]+1, j] > f_rank[i+1, j]
                                    match_or_unmatch[f_matched[j]] = true
                                    f_matched[j] = i
                                    m_matched[i] = j
                                    match_or_unmatch[i] = false
                                    break
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    for (i,t) in enumerate(f_matched)
        if t == m_num + 1
            f_matched[i] = 0
        end
    end

     for (i,t) in enumerate(m_matched)
        if t == f_num + 1
            m_matched[i] = 0
        end
    end

    return m_matched, f_matched
end


text:
module Poly

export LoopKernel, Domain, Instruction

"""
Represents an instruction
Args:
    iname: symbol representing instruction
    body: expression in the instruction
    dependencies: ids of instructions that this instruction depends on
    cond: expression for condition of this instruction, if instruction is in an if/elseif/else block

Example Instruction:
>>> instructions = [Instruction(:mult, :(out[i, j] += A[i, k] * B[k, j], :()), Set()),
                    Instruction(:double, :(out[i, j] *= 2), Set([:mult])), :()]
"""
struct Instruction
    iname::Symbol
    body::Expr
    dependencies::Set{Symbol}
    cond::Expr
end

# Base.show(io::IO, d::Instruction) = print(io, d.iname)

"""
Represents an iteration domain
Args:
    iname: symbol representing the simple set
    lowerbound: lowerbound of the simple set (inclusive)
    upperbound: upperbound of the simple set (inclusive)
    step: step of loop iteration
    dependencies: dependencies of this domain
    instructions: instructions in the body of this domain (assigned after initialization)

Example Domain:
>>> domains = [Domain(:i, 1, :(size(out, 1)), 1, Set(), []),
               Domain(:j, 1, :(size(out, 2)), 1, Set(), []),
               Domain(:k, 1, :(size(A, 2)), 1, Set(), [])]
"""
struct Domain
    iname::Symbol
    lowerbound::Union{Number, Symbol, Expr}
    upperbound::Union{Number, Symbol, Expr}
    step::Union{Number, Symbol, Expr}
    dependencies::Set{Symbol}
    instructions::Vector{Union{Instruction, Domain}}
end

# Base.show(io::IO, d::Domain) = print(io, d.iname)

"""
Represents a kernel
Args:
    instructions: list of instructions in kernel. Includes assignments, variables, function calls, etc
    domains: list of ISL SimpleSet domains
    args: arguments to kernel
    argtypes: types of arguments to kernel
    consts: constants in kernel

Example LoopKernel:
>>> kern = LoopKernel(instructions, domains, [:out, :A, :B], [Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}], [])
"""
struct LoopKernel
    instructions::Vector{Instruction}
    domains::Vector{Domain}
    args::Vector{Symbol}
    argtypes::Vector{Type}
    consts::Vector{Symbol}
end

"""
Represents the AST of instruction and loop ordering for a kernel
Args:
    self: either Instruction or Domain object (or nothing)
    children: list of dependent Instructions or Domains
    parents: list of Instructions or Domains dependent on
"""
struct AST
    self::Union{Domain, Instruction, Nothing}
    children::Vector{AST}
    parents::Vector{AST}
end

function get_id(ast::AST)
    if ast.self == nothing
        return nothing
    end
    return ast.self.iname
end

using AbstractTrees
AbstractTrees.children(d::AST) = d.children
AbstractTrees.parentlinks(::AST) = AbstractTrees.StoredParents()
AbstractTrees.printnode(io::IO, d::AST) = print(io, get_id(d))
Base.show(io::IO, d::AST) = print_tree(io, d)

include("compilation.jl")
include("macros.jl")
include("polyhedral.jl")
include("striding_analysis.jl")
include("tiling.jl")

end


text:
module Euclid
import ProjectiveGeometricAlgebra3d as PGA
using .PGA
using .PGA: sandwich
using .PGA: MultiVector3, MultiVector1, MultiVector2, MultiVectorEven, ∨
using .PGA: normalize, scalartype, normalize_ideal
using .PGA: random
using ArgCheck

export Point, Plane, Direction, Line
export meet, orthogonal_through, join, distance, project
export Motor, translator, rotator
export random

function pga end

struct Point{T}
    data::MultiVector3{T}
    function Point(data::MultiVector3)
        data1 = normalize(data)
        T = scalartype(typeof(data1))
        return new{T}(data1)
    end
end
Point(pt::Point) = pt
function Point(itr)
    @argcheck length(itr) == 3
    x,y,z = itr
    Point(x,y,z)
end
function Point(x,y,z)
    Point(MultiVector3(e023=-x, e013=y, e012=-z, e123=oneunit(x)))
end

struct Direction{T}
    data::MultiVector3{T}
    function Direction(data::MultiVector3)
        data1 = normalize_ideal(data)
        T = scalartype(typeof(data1))
        return new{T}(data1)
    end
end
Direction(o::Direction) = o
function Direction(itr)
    @argcheck length(itr) == 3
    x,y,z = itr
    Direction(x,y,z)
end
function Direction(x,y,z)
    data::MultiVector3 = MultiVector3(e023=-x, e013=y, e012=-z)
    Direction(data)
end

struct Line{T}
    data::MultiVector2{T}
    function Line(data::MultiVector2)
        data1 = normalize(data)
        T = scalartype(typeof(data1))
        return new{T}(data1)
    end
end
Line(o::Line) = o
function Line(;point=nothing, direction=nothing)
    has_point = point !== nothing
    has_direction = direction !== nothing
    if has_point && has_direction
        Line_point_direction(point, direction)
    else
        throw(ArgumentError("Cannot construct Line from arguments"))
    end
end
function Line_point_direction(point, direction)
    data::MultiVector2 = pga(Point(point)) ∨ pga(Direction(direction))
    Line(data)
end

struct Plane{T}
    data::MultiVector1{T}
    function Plane(data::MultiVector1)
        data1 = normalize(data)
        T = scalartype(typeof(data1))
        return new{T}(data1)
    end
end
Plane(o::Plane) = o
function Plane(;point, directions)
    Plane_point_directions(point, directions)
end
function Plane_point_directions(point, directions)
    @argcheck length(directions) == 2
    dir1, dir2 = directions
    pt::MultiVector3 = pga(Point(point))
    d1::MultiVector3 = pga(Direction(dir1))
    d2::MultiVector3 = pga(Direction(dir2))
    data::MultiVector1 = pt ∨ d1 ∨ d2
    Plane(data)
end

################################################################################
##### Motor
################################################################################
struct Motor{T}
    data::MultiVectorEven{T}
end
Motor(m::Motor) = m

function rotator(;axis, angle)
    l = Line(axis)
    s,c = sincos(1//2*angle)
    data = c + s*pga(l)
    return Motor(data)
end

function translator end

(m::Motor)(o::Point) = Point(sandwich(pga(o), pga(m)))
(m::Motor)(o::Direction) = Direction(sandwich(pga(o), pga(m)))
(m::Motor)(o::Line)  = Line(sandwich(pga(o), pga(m)))
(m::Motor)(o::Plane) = Plane(sandwich(pga(o), pga(m)))
(m::Motor)(o::Motor) = Motor(sandwich(pga(o), pga(m)))

Base.:(∘)(m2::Motor, m1::Motor) = Motor(pga(m1) * pga(m2))

pgatype(::Type{Point{T}}     ) where {T} = MultiVector3{T}
pgatype(::Type{Direction{T}} ) where {T} = MultiVector3{T}
pgatype(::Type{Line{T}}      ) where {T} = MultiVector2{T}
pgatype(::Type{Plane{T}}     ) where {T} = MultiVector1{T}
pgatype(::Type{Motor{T}}     ) where {T} = MultiVectorEven{T}

for Geo in [:Point, :Direction, :Line, :Plane, :Motor]
    @eval function pga(x::$Geo)
        x.data
    end
    @eval function PGA.random(rng::PGA.AbstractRNG, G::Type{<:$Geo})
        $Geo(PGA.random(rng, pgatype(G)))::G
    end
    @eval function Base.isapprox(o1::$Geo, o2::$Geo; kw...)
        isapprox(pga(o1), pga(o2); kw...)
    end
end

meet(o1::Plane, o2::Plane) = Line( pga(o1) ∧ pga(o2))
meet(o1::Line, o2::Plane)  = Point(pga(o1) ∧ pga(o2))
meet(o1::Plane, o2::Line)  = Point(pga(o1) ∧ pga(o2))

# TODO join
join(o::Union{Point, Line, Plane}) = o
join(o1::Point, o2::Point) = Line(pga(o1) ∨ pga(o2))
join(o1::Point, o2::Line) = Plane(pga(o1) ∨ pga(o2))
join(o1::Line, o2::Point) = Plane(pga(o1) ∨ pga(o2))
join(o1::Point, o2::Point, o3::Point) = Plane(pga(o1) ∨ pga(o2) ∨ pga(o3))

orthogonal_through(orth::Plane, through::Line) = Plane(inner(orth, through))
orthogonal_through(orth::Plane, through::Point) = Line(inner(orth, through))
orthogonal_through(orth::Line, through::Point) = Plane(inner(orth, through))

project(onto::Union{Point, Line, Plane}, obj::Point) = Point(PGA.project(pga(onto), pga(obj)))
project(onto::Plane, obj::Line) = Line(PGA.project(pga(onto), pga(obj)))

################################################################################
##### metric
################################################################################
function cos_angle(o1::Plane, o2::Plane)
    inner(pga(o1), pga(o2)).e
end
function cos_angle(o1::Line, o2::Plane)
    inner(pga(o1), pga(o2)).e
end
function cos_angle(o1::Plane, o2::Line)
    inner(pga(o1), pga(o2)).e
end

function distance(o1::Point, o2::Point)
    norm(pga(o1) ∨ pga(o2))
end

end #module


text:
Base.@kwdef struct TableMetaData
    name::Union{Nothing,String} = nothing
    id::Union{Nothing,String} = nothing
    provider::Union{Nothing,String} = nothing
    reference::Union{Nothing,String} = nothing
    content_type::Union{Nothing,String} = nothing
    description::Union{Nothing,String} = nothing
    comments::Union{Nothing,String} = nothing
    source_path::Union{Nothing,String} = nothing
end

mort_max_dur = 121
mort_max_issue_age = 120


text:
function evaluate(obj::ContinuedFraction, arg0::jdouble)
    return jcall(obj, "evaluate", jdouble, (jdouble,), arg0)
end

function evaluate(obj::ContinuedFraction, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "evaluate", jdouble, (jdouble, jdouble), arg0, arg1)
end

function evaluate(obj::ContinuedFraction, arg0::jdouble, arg1::jdouble, arg2::jint)
    return jcall(obj, "evaluate", jdouble, (jdouble, jdouble, jint), arg0, arg1, arg2)
end

function evaluate(obj::ContinuedFraction, arg0::jdouble, arg1::jint)
    return jcall(obj, "evaluate", jdouble, (jdouble, jint), arg0, arg1)
end


text:
using Blink
using Base.Test

# open window and wait for it to initialize
w = Window(Blink.@d(:show => false)); sleep(5.0)

# make sure the window is really active
@test_approx_eq @js(w, Math.log(10)) log(10)


text:
"""
    SurfaceFluxes

## Interface
  - [`surface_conditions`](@ref) computes
    - Monin-Obukhov length
    - Potential temperature flux (if not given) using Monin-Obukhov theory
    - transport fluxes using Monin-Obukhov theory
    - friction velocity/temperature scale/tracer scales
    - exchange coefficients

## References
 - [Nishizawa2018](@cite)
 - [Byun1990](@cite)

"""
module SurfaceFluxes

import RootSolvers
const RS = RootSolvers

import KernelAbstractions
const KA = KernelAbstractions

using DocStringExtensions
const DSE = DocStringExtensions

import Thermodynamics
const TD = Thermodynamics

import CLIMAParameters
const CP = CLIMAParameters
const CPP = CP.Planet
const APS = CP.AbstractEarthParameterSet
const CPSGS = CP.SubgridScale

import StaticArrays
const SA = StaticArrays

include("UniversalFunctions.jl")
import .UniversalFunctions
const UF = UniversalFunctions

abstract type SolverScheme end
struct FVScheme <: SolverScheme end
struct FDScheme <: SolverScheme end

# Allow users to skip error on non-convergence
# by importing:
# ```julia
# import SurfaceFluxes
# SurfaceFluxes.error_on_non_convergence() = false
# ```
# Error on convergence must be the default
# behavior because this can result in printing
# very large logs resulting in CI to seemingly hang.
error_on_non_convergence() = true

"""
    SurfaceFluxConditions

Surface flux conditions, returned from `surface_conditions`.

# Fields

$(DSE.FIELDS)
"""
struct SurfaceFluxConditions{FT <: AbstractFloat}
    L_MO::FT
    shf::FT
    lhf::FT
    buoy_flux::FT
    ρτxz::FT
    ρτyz::FT
    ustar::FT
    Cd::FT
    Ch::FT
end

function Base.show(io::IO, sfc::SurfaceFluxConditions)
    println(io, "----------------------- SurfaceFluxConditions")
    println(io, "L_MO                   = ", sfc.L_MO)
    println(io, "Sensible Heat Flux     = ", sfc.shf)
    println(io, "Latent Heat Flux       = ", sfc.lhf)
    println(io, "Buoyancy Flux          = ", sfc.buoy_flux)
    println(io, "Friction velocity u⋆   = ", sfc.ustar)
    println(io, "C_drag                 = ", sfc.Cd)
    println(io, "C_heat                 = ", sfc.Ch)
    println(io, "-----------------------")
end

"""
    SurfaceValues

Input container for state variables at the ground level.

# Fields

$(DSE.FIELDS)
"""
struct SurfaceValues{FT <: AbstractFloat, A, TS <: TD.ThermodynamicState}
    z::FT
    u::A
    ts::TS
end

"""
    InteriorValues

Input container for state variables at the first interior node.

# Fields

$(DSE.FIELDS)
"""
struct InteriorValues{FT <: AbstractFloat, A, TS <: TD.ThermodynamicState}
    z::FT
    u::A
    ts::TS
end

abstract type AbstractSurfaceConditions{FT <: AbstractFloat, VI <: InteriorValues, VS <: SurfaceValues} end

"""
    Fluxes

Input container for state variables, latent and sensible heat fluxes roughness lengths,
initial obukhov length and gustiness.

# Fields

$(DSE.FIELDS)
"""
Base.@kwdef struct Fluxes{FT, VI, VS} <: AbstractSurfaceConditions{FT, VI, VS}
    state_in::VI
    state_sfc::VS
    shf::FT
    lhf::FT
    z0m::FT
    z0b::FT
    L_MO_init::FT = FT(100)
    gustiness::FT = FT(1)
end

function Fluxes{FT}(; state_in, state_sfc, kwargs...) where {FT}
    types = (FT, typeof(state_in), typeof(state_sfc))
    return Fluxes{types...}(; state_in, state_sfc, kwargs...)
end

"""
    FluxesAndFrictionVelocity

Input container, given surface state variables, latent and sensible heat fluxes,
and the friction velocity, roughness lengths,
initial obukhov length and gustiness.

# Fields

$(DSE.FIELDS)
"""
Base.@kwdef struct FluxesAndFrictionVelocity{FT, VI, VS} <: AbstractSurfaceConditions{FT, VI, VS}
    state_in::VI
    state_sfc::VS
    shf::FT
    lhf::FT
    ustar::FT
    z0m::FT
    z0b::FT
    gustiness::FT = FT(1)
end

function FluxesAndFrictionVelocity{FT}(; state_in, state_sfc, kwargs...) where {FT}
    types = (FT, typeof(state_in), typeof(state_sfc))
    return FluxesAndFrictionVelocity{types...}(; state_in, state_sfc, kwargs...)
end

"""
    Coefficients

Input container, given surface state variables, and exchange coefficients,roughness lengths,
initial obukhov length and gustiness.

# Fields

$(DSE.FIELDS)
"""
Base.@kwdef struct Coefficients{FT, VI, VS} <: AbstractSurfaceConditions{FT, VI, VS}
    state_in::VI
    state_sfc::VS
    Cd::FT
    Ch::FT
    z0m::FT
    z0b::FT
    gustiness::FT = FT(1)
end

function Coefficients{FT}(; state_in, state_sfc, kwargs...) where {FT}
    types = (FT, typeof(state_in), typeof(state_sfc))
    return Coefficients{types...}(; state_in, state_sfc, kwargs...)
end

"""
    ValuesOnly

Input container, given only surface state variables, roughness lengths,
initial obukhov length and gustiness.

# Fields

$(DSE.FIELDS)
"""
Base.@kwdef struct ValuesOnly{FT, VI, VS} <: AbstractSurfaceConditions{FT, VI, VS}
    state_in::VI
    state_sfc::VS
    z0m::FT
    z0b::FT
    L_MO_init::FT = FT(100)
    gustiness::FT = FT(1)
end

function ValuesOnly{FT}(; state_in, state_sfc, kwargs...) where {FT}
    types = (FT, typeof(state_in), typeof(state_sfc))
    return ValuesOnly{types...}(; state_in, state_sfc, kwargs...)
end

ts_in(sc::AbstractSurfaceConditions) = sc.state_in.ts
ts_sfc(sc::AbstractSurfaceConditions) = sc.state_sfc.ts

z_in(sc::AbstractSurfaceConditions) = sc.state_in.z
z_sfc(sc::AbstractSurfaceConditions) = sc.state_sfc.z
Δz(sc::AbstractSurfaceConditions) = z_in(sc) - z_sfc(sc)

z0(sc::AbstractSurfaceConditions, ::UF.HeatTransport) = sc.z0b
z0(sc::AbstractSurfaceConditions, ::UF.MomentumTransport) = sc.z0m

Δu1(sc::AbstractSurfaceConditions) = sc.state_in.u[1] - sc.state_sfc.u[1]
Δu2(sc::AbstractSurfaceConditions) = sc.state_in.u[2] - sc.state_sfc.u[2]

qt_in(sc::AbstractSurfaceConditions) = TD.total_specific_humidity(ts_in(sc))
qt_sfc(sc::AbstractSurfaceConditions) = TD.total_specific_humidity(ts_sfc(sc))
Δqt(sc::AbstractSurfaceConditions) = qt_in(sc) - qt_sfc(sc)

u_in(sc::AbstractSurfaceConditions) = sc.state_in.u
u_sfc(sc::AbstractSurfaceConditions) = sc.state_sfc.u

function windspeed(sc::AbstractSurfaceConditions)
    return max(hypot(Δu1(sc), Δu2(sc)), sc.gustiness)
end

"""
    surface_conditions(
        param_set::CLIMAParameters.AbstractParameterSet,
        sc::SurfaceFluxes.AbstractSurfaceConditions{FT},
        uft::UniversalFunctions.AbstractUniversalFunctionType,
        scheme::SurfaceFluxes.SolverScheme = FVScheme();
        tol::RootSolvers.AbstractTolerance =
            RootSolvers.SolutionTolerance(max(FT(1e-2), sqrt(eps(FT)))),
        maxiter::Int = 100,
    ) where {FT}

The main user facing function of the module.
It computes the surface conditions
based on the Monin-Obukhov similarity functions. Requires
information about thermodynamic parameters (`param_set`)
the surface state `sc`, the universal function type and
the discretisation `scheme`.

Result struct of type SurfaceFluxConditions{FT} contains:
  - L_MO:   Monin-Obukhov lengthscale
  - shf:    Sensible Heat Flux
  - lhf:    Latent Heat Flux
  - ρτxz:   Momentum Flux (Eastward component)
  - ρτyz:   Momentum Flux (Northward component)
  - ustar:  Friction velocity
  - Cd:     Momentum Exchange Coefficient
  - Ch:     Thermal Exchange Coefficient
"""
function surface_conditions(
    param_set::APS,
    sc::AbstractSurfaceConditions{FT},
    uft::UF.AUFT,
    scheme::SolverScheme = FVScheme();
    tol::RS.AbstractTolerance = RS.SolutionTolerance(max(FT(1e-2), sqrt(eps(FT)))),
    maxiter::Int = 100,
) where {FT}
    # Compute output variables for SurfaceFluxConditions
    L_MO = obukhov_length(param_set, sc, uft, scheme; tol, maxiter)
    ustar = compute_ustar(param_set, L_MO, sc, uft, scheme)
    Cd = momentum_exchange_coefficient(param_set, L_MO, sc, uft, scheme)
    Ch = heat_exchange_coefficient(param_set, L_MO, sc, uft, scheme)
    shf = sensible_heat_flux(param_set, Ch, sc, scheme)
    lhf = latent_heat_flux(param_set, Ch, sc, scheme)
    buoy_flux = compute_buoyancy_flux(param_set, shf, lhf, ts_in(sc), ts_sfc(sc), scheme)
    ρτxz, ρτyz = momentum_fluxes(param_set, Cd, sc, scheme)
    return SurfaceFluxConditions{FT}(L_MO, shf, lhf, buoy_flux, ρτxz, ρτyz, ustar, Cd, Ch)
end

"""
    obukhov_length(sfc::SurfaceFluxConditions)

    obukhov_length( # internal method
        param_set::AbstractParameterSet,
        sc::AbstractSurfaceConditions,
        uft,
        scheme;
        tol::RootSolvers.AbstractTolerance =
             RootSolvers.SolutionTolerance(max(1e-2, sqrt(eps))),
        maxiter::Int = 100
    )

Compute and return the Monin-Obukhov lengthscale (LMO).

The internal method for computing LMO depends on the
particular surface condition `sc`, of which there are
several options:

 - `FluxesAndFrictionVelocity`
 - `Coefficients`

## `AbstractSurfaceConditions` (fallback)

The Monin-Obukhov length is computed by solving a non-linear
equation given a tolerance `tol` and maximum iterations `maxiter`.

## `FluxesAndFrictionVelocity`

Surface fluxes and friction velocity are known.
Iterations are not needed to determine LMO.

## `Coefficients`

Exchange coefficients are known.
Iterations are not needed to determine LMO.
"""
function obukhov_length end

obukhov_length(sfc::SurfaceFluxConditions) = sfc.L_MO

function obukhov_length(
    param_set,
    sc::AbstractSurfaceConditions{FT},
    uft::UF.AUFT,
    scheme;
    tol::RS.AbstractTolerance = RS.SolutionTolerance(max(FT(1e-2), sqrt(eps(FT)))),
    maxiter::Int = 100,
) where {FT}

    function root_l_mo(x_lmo)
        root = x_lmo - local_lmo(param_set, x_lmo, sc, uft, scheme)
        return root
    end
    sol = RS.find_zero(root_l_mo, RS.NewtonsMethodAD(sc.L_MO_init), RS.CompactSolution(), tol, maxiter)
    if sol.converged
        L_MO = sol.root
    else
        if error_on_non_convergence()
            KA.@print("-----------------------------------------\n")
            KA.@print("maxiter reached in SurfaceFluxes.jl:\n")
            KA.@print(", T_in = ", TD.air_temperature(ts_in(sc)))
            KA.@print(", T_sfc = ", TD.air_temperature(ts_sfc(sc)))
            KA.@print(", q_in = ", TD.total_specific_humidity(ts_in(sc)))
            KA.@print(", q_sfc = ", TD.total_specific_humidity(ts_sfc(sc)))
            KA.@print(", u_in = ", u_in(sc))
            KA.@print(", u_sfc = ", u_sfc(sc))
            error("Unconverged Surface Fluxes.")
        else
            # KA.@print("Warning: Unconverged Surface Fluxes. Returning neutral condition solution \n")
            L_MO = sol.root
        end
    end
    return L_MO
end

function obukhov_length(param_set, sc::FluxesAndFrictionVelocity{FT}, uft::UF.AUFT, scheme; kwargs...) where {FT}
    return -sc.ustar^3 / FT(CPSGS.von_karman_const(param_set)) / compute_buoyancy_flux(param_set, sc, scheme)
end

function obukhov_length(param_set, sc::Coefficients{FT}, uft::UF.AUFT, scheme; kwargs...) where {FT}
    lhf = latent_heat_flux(param_set, sc.Ch, sc, scheme)
    shf = sensible_heat_flux(param_set, sc.Ch, sc, scheme)
    ustar = sqrt(sc.Cd) * windspeed(sc)
    buoyancy_flux = compute_buoyancy_flux(param_set, shf, lhf, ts_in(sc), ts_sfc(sc), scheme)
    return -ustar^3 / FT(CPSGS.von_karman_const(param_set)) / buoyancy_flux
end

"""
    local_lmo(param_set, x_lmo, sc, uft, scheme)

Helper function for the iterative solver with the Monin-Obukhov length equation with buoyancy flux.
"""
function local_lmo(param_set, x_lmo, sc::Fluxes, uft::UF.AUFT, scheme)
    κ = CPSGS.von_karman_const(param_set)
    u_scale = compute_physical_scale_coeff(param_set, sc, x_lmo, UF.MomentumTransport(), uft, scheme)
    return -(windspeed(sc) * u_scale)^3 / κ / compute_buoyancy_flux(param_set, sc, scheme)
end

"""
    local_lmo(param_set, x_lmo, sc, uft, scheme)

Helper function for the iterative solver with the Monin-Obukhov
length equation with buoyancy star.
"""
function local_lmo(param_set, x_lmo, sc::ValuesOnly, uft::UF.AUFT, scheme)
    κ = CPSGS.von_karman_const(param_set)
    u_scale = compute_physical_scale_coeff(param_set, sc, x_lmo, UF.MomentumTransport(), uft, scheme)
    return (windspeed(sc) * u_scale)^2 / κ / compute_bstar(param_set, x_lmo, sc, uft, scheme)
end

"""
    compute_buoyancy_flux(param_set, shf, lhf, ts_in, ts_sfc, scheme)

Returns the buoyancy flux when the surface fluxes are known.
"""
function compute_buoyancy_flux(param_set, shf::FT, lhf::FT, ts_in, ts_sfc, scheme) where {FT}

    grav::FT = CPP.grav(param_set)
    ε_vd::FT = CPP.molmass_ratio(param_set)

    cp_m = TD.cp_m(ts_in)
    L_v = TD.latent_heat_vapor(ts_in)
    ρ_sfc = TD.air_density(ts_sfc)
    T_in = TD.air_temperature(ts_in)
    return grav / ρ_sfc * (shf / cp_m / T_in + (ε_vd - 1) * lhf / L_v)
end

function compute_buoyancy_flux(param_set, sc::Union{FluxesAndFrictionVelocity, Fluxes}, scheme)
    return compute_buoyancy_flux(param_set, sc.shf, sc.lhf, ts_in(sc), ts_sfc(sc), scheme)
end

"""
    compute_bstar(param_set, L_MO, sc, uft, scheme)

Returns buoyancy star based on known air densities.
"""
function compute_bstar(param_set, L_MO, sc::AbstractSurfaceConditions{FT}, uft::UF.AUFT, scheme) where {FT}
    uf = UF.universal_func(uft, L_MO, param_set)
    grav::FT = CPP.grav(param_set)

    ρ_in = TD.air_density(ts_in(sc))
    ρ_sfc = TD.air_density(ts_sfc(sc))

    ρ_star = compute_physical_scale_coeff(param_set, sc, L_MO, UF.HeatTransport(), uft, scheme) * (ρ_in - ρ_sfc)

    return -grav / ρ_in * ρ_star
end

"""
    compute_bstar(param_set, L_MO, sc, uft, scheme)

Returns buoyancy star based on known friction velocity  and fluxes.
"""
compute_bstar(param_set, L_MO, sc::FluxesAndFrictionVelocity, uft, scheme) =
    -compute_buoyancy_flux(param_set, sc, scheme) / sc.ustar

"""
    compute_bstar(param_set, L_MO, sc, uft, scheme)

Return buoyancy star based on known fluxes.
"""
compute_bstar(param_set, L_MO, sc::Fluxes, uft, scheme) =
    -compute_buoyancy_flux(param_set, sc, scheme) / compute_ustar(param_set, L_MO, sc, uft, scheme)

#=
    compute_ustar(
        param_set::AbstractParameterSet,
        L_MO,
        sc::AbstractSufaceCondition,
        uft,
        scheme
    )

Return the friction velocity. This method is dispatched
by the surface condition:

## `sc::FluxesAndFrictionVelocity`

Friction velocity is known.

## `sc::Fluxes`

Compute given the Monin-Obukhov lengthscale.

## `sc::Coefficients`

Compute given the exchange coefficients.

## `sc::ValuesOnly`
Compute given the Monin-Obukhov lengthscale.
=#
function compute_ustar end

compute_ustar(param_set, L_MO, sc::FluxesAndFrictionVelocity, uft, scheme) = sc.ustar

compute_ustar(param_set, L_MO, sc::Fluxes, uft, scheme) =
    windspeed(sc) * compute_physical_scale_coeff(param_set, sc, L_MO, UF.MomentumTransport(), uft, scheme)

compute_ustar(param_set, L_MO, sc::Coefficients, uft, scheme) = sqrt(sc.Cd) * (windspeed(sc))

compute_ustar(param_set, L_MO, sc::ValuesOnly, uft, scheme) =
    windspeed(sc) * compute_physical_scale_coeff(param_set, sc, L_MO, UF.MomentumTransport(), uft, scheme)

"""
    momentum_exchange_coefficient(param_set, L_MO, sc, uft, scheme)

Compute and return Cd, the momentum exchange coefficient, given the
Monin-Obukhov lengthscale.
"""
function momentum_exchange_coefficient(
    param_set,
    L_MO,
    sc::Union{Fluxes, ValuesOnly, FluxesAndFrictionVelocity},
    uft::UF.AUFT,
    scheme,
)
    κ = CPSGS.von_karman_const(param_set)
    transport = UF.MomentumTransport()
    FT = eltype(L_MO)
    ρ_in = TD.air_density(ts_in(sc))
    ρ_sfc = TD.air_density(ts_sfc(sc))
    Δρ = ρ_in - ρ_sfc
    ustar = compute_ustar(param_set, L_MO, sc, uft, scheme)
    if Δρ ≈ FT(0)
        Cd = (κ / log(Δz(sc) / z0(sc, transport)))^2
    else
        Cd = ustar^2 / windspeed(sc)^2
    end
    return Cd
end

"""
    momentum_exchange_coefficient(param_set, L_MO, sc, uft, scheme)

Return Cd, the momentum exchange coefficient.
"""
function momentum_exchange_coefficient(param_set, L_MO, sc::Coefficients, uft, scheme)
    return sc.Cd
end

"""
    heat_exchange_coefficient(param_set, L_MO, sc, uft, scheme)

Compute and return Ch, the heat exchange coefficient given the
Monin-Obukhov lengthscale.
"""
function heat_exchange_coefficient(
    param_set,
    L_MO,
    sc::Union{Fluxes, ValuesOnly, FluxesAndFrictionVelocity},
    uft,
    scheme,
)
    FT = eltype(L_MO)
    Δθ = TD.dry_pottemp(ts_in(sc)) - TD.dry_pottemp(ts_sfc(sc))
    ustar = compute_ustar(param_set, L_MO, sc, uft, scheme)
    θstar = Δθ * compute_physical_scale_coeff(param_set, sc, L_MO, UF.HeatTransport(), uft, scheme)
    if Δθ ≈ FT(0) || windspeed(sc) ≈ FT(0)
        Ch = FT(0)
    else
        Ch = ustar * θstar / windspeed(sc) / Δθ
    end
    return Ch
end

"""
    heat_exchange_coefficient(param_set, L_MO, sc, uft, scheme)

Return Ch, the heat exchange coefficient given the
Monin-Obukhov lengthscale.
"""
function heat_exchange_coefficient(param_set, L_MO, sc::Coefficients, uft, scheme)
    return sc.Ch
end

"""
    momentum_fluxes(param_set, Cd, sc, scheme)

Compute and return the momentum fluxes
## Arguments
  - param_set: Abstract Parameter Set containing physical, thermodynamic parameters.
  - Cd: Momentum exchange coefficient
  - sc: Container for surface conditions based on known combination
        of the state vector, and {fluxes, friction velocity, exchange coefficients} for a given experiment
  - scheme: Discretization scheme (currently supports FD and FV)
"""
function momentum_fluxes(param_set, Cd, sc::AbstractSurfaceConditions, scheme)
    ρ_sfc = TD.air_density(ts_sfc(sc))
    ρτxz = -ρ_sfc * Cd * Δu1(sc) * windspeed(sc)
    ρτyz = -ρ_sfc * Cd * Δu2(sc) * windspeed(sc)
    return (ρτxz, ρτyz)
end

"""
    sensible_heat_flux(param_set, Ch, sc, scheme)

Compute and return the sensible heat fluxes
## Arguments
  - param_set: Abstract Parameter Set containing physical, thermodynamic parameters.
  - Ch: Thermal exchange coefficient
  - sc: Container for surface conditions based on known combination
        of the state vector, and {fluxes, friction velocity, exchange coefficients} for a given experiment
  - scheme: Discretization scheme (currently supports FD and FV)
"""
function sensible_heat_flux(param_set, Ch::FT, sc::Union{ValuesOnly, Coefficients}, scheme) where {FT}
    grav::FT = CPP.grav(param_set)
    cp_d::FT = CPP.cp_d(param_set)
    R_d::FT = CPP.R_d(param_set)
    T_0::FT = CPP.T_0(param_set)
    cp_m = TD.cp_m(ts_in(sc))
    ρ_sfc = TD.air_density(ts_sfc(sc))
    T_in = TD.air_temperature(ts_in(sc))
    T_sfc = TD.air_temperature(ts_sfc(sc))
    ΔT = T_in - T_sfc
    hd_sfc = cp_d * (T_sfc - T_0) + R_d * T_0
    ΔΦ = grav * Δz(sc)
    E = evaporation(sc, param_set, Ch)
    return -ρ_sfc * Ch * windspeed(sc) * (cp_m * ΔT + ΔΦ) - (hd_sfc) * E
end

function evaporation(sc::Union{ValuesOnly, Coefficients}, param_set, Ch)
    ρ_sfc = TD.air_density(ts_sfc(sc))
    return -ρ_sfc * Ch * windspeed(sc) * Δqt(sc)
end

"""
    sensible_heat_flux(param_set, Ch, sc, scheme)

In cases where surface fluxes are known,
return the known sensible heat flux.
"""
function sensible_heat_flux(param_set, Ch, sc::Union{Fluxes, FluxesAndFrictionVelocity}, scheme)
    return sc.shf
end

"""
    latent_heat_flux(param_set, Ch, sc, scheme)

In cases where surface fluxes are known,
return the known latent heat flux.
"""
function latent_heat_flux(param_set, L_MO, sc::Union{Fluxes, FluxesAndFrictionVelocity}, scheme)
    return sc.lhf
end

"""
    latent_heat_flux(param_set, Ch, sc, scheme)

Compute and return the latent heat flux
## Arguments
  - param_set: Abstract Parameter Set containing physical, thermodynamic parameters.
  - Ch: Thermal exchange coefficient
  - sc: Container for surface conditions based on known combination
        of the state vector, and {fluxes, friction velocity, exchange coefficients} for a given experiment
  - scheme: Discretization scheme (currently supports FD and FV)
"""
function latent_heat_flux(param_set, Ch::FT, sc::Union{ValuesOnly, Coefficients}, scheme) where {FT}
    grav::FT = CPP.grav(param_set)
    ρ_sfc = TD.air_density(ts_sfc(sc))
    cp_v::FT = CPP.cp_v(param_set)
    Lv_0::FT = CPP.LH_v0(param_set)
    T_0::FT = CPP.T_0(param_set)
    T_sfc = TD.air_temperature(ts_sfc(sc))
    hv_sfc = cp_v * (T_sfc - T_0) + Lv_0
    Φ_sfc = grav * z_sfc(sc)
    E = evaporation(sc, param_set, Ch)
    lhf = (hv_sfc + Φ_sfc) * E
    return lhf
end

"""
    compute_physical_scale_coeff(param_set, sc, L_MO, transport, uft, ::FVScheme)

Computes the coefficient for the physical scale of a variable based on Nishizawa(2018)
for the FV scheme.

## Arguments
  - param_set: Abstract Parameter Set containing physical, thermodynamic parameters.
  - sc: Container for surface conditions based on known combination
        of the state vector, and {fluxes, friction velocity, exchange coefficients} for a given experiment
  - L_MO: Monin-Obukhov length
  - transport: Transport type, (e.g. Momentum or Heat, used to determine physical scale coefficients)
  - uft: A Universal Function type, (returned by, e.g. Businger())
  - scheme: Discretization scheme (currently supports FD and FV)
"""
function compute_physical_scale_coeff(
    param_set::APS,
    sc::AbstractSurfaceConditions,
    L_MO::FT,
    transport,
    uft,
    ::FVScheme,
) where {FT}
    von_karman_const::FT = CPSGS.von_karman_const(param_set)
    uf = UF.universal_func(uft, L_MO, param_set)
    _π_group = FT(UF.π_group(uf, transport))
    _π_group⁻¹ = (1 / _π_group)
    R_z0 = 1 - z0(sc, transport) / Δz(sc)
    denom1 = log(Δz(sc) / z0(sc, transport))
    denom2 = -UF.Psi(uf, Δz(sc) / uf.L, transport)
    denom3 = z0(sc, transport) / Δz(sc) * UF.Psi(uf, z0(sc, transport) / uf.L, transport)
    denom4 = R_z0 * (UF.psi(uf, z0(sc, transport) / uf.L, transport) - 1)
    Σterms = denom1 + denom2 + denom3 + denom4
    return _π_group⁻¹ * von_karman_const / Σterms
end

"""
    compute_physical_scale_coeff(param_set, sc, L_MO, transport, uft, ::FDScheme)

Computes the coefficient for the physical scale of a variable based on Byun(1990)
for the Finite Differneces scheme.

## Arguments
  - param_set: Abstract Parameter Set containing physical, thermodynamic parameters.
  - sc: Container for surface conditions based on known combination
        of the state vector, and {fluxes, friction velocity, exchange coefficients} for a given experiment
  - L_MO: Monin-Obukhov length
  - transport: Transport type, (e.g. Momentum or Heat, used to determine physical scale coefficients)
  - uft: A Universal Function type, (returned by, e.g., Businger())
  - scheme: Discretization scheme (currently supports FD and FV)
"""
function compute_physical_scale_coeff(
    param_set,
    sc::AbstractSurfaceConditions,
    L_MO::FT,
    transport,
    uft,
    ::FDScheme,
) where {FT}
    von_karman_const::FT = CPSGS.von_karman_const(param_set)
    uf = UF.universal_func(uft, L_MO, param_set)
    _π_group = FT(UF.π_group(uf, transport))
    _π_group⁻¹ = (1 / _π_group)
    denom1 = log(Δz(sc) / z0(sc, transport))
    denom2 = -UF.psi(uf, Δz(sc) / uf.L, transport)
    denom3 = UF.psi(uf, z0(sc, transport) / uf.L, transport)
    Σterms = denom1 + denom2 + denom3
    return _π_group⁻¹ * von_karman_const / Σterms
end

"""
    recover_profile(param_set, sc, L_MO, Z, X_in, X_sfc, transport, uft, scheme)

Recover profiles of variable X given values of Z coordinates. Follows Nishizawa equation (21,22)
## Arguments
  - param_set: Abstract Parameter Set containing physical, thermodynamic parameters.
  - sc: Container for surface conditions based on known combination
        of the state vector, and {fluxes, friction velocity, exchange coefficients} for a given experiment
  - L_MO: Monin-Obukhov length
  - Z: Z coordinate(s) (within surface layer) for which variable values are required
  - X_in,X_sfc: For variable X, values at interior and surface nodes
  - transport: Transport type, (e.g. Momentum or Heat, used to determine physical scale coefficients)
  - uft: A Universal Function type, (returned by, e.g., Businger())
  - scheme: Discretization scheme (currently supports FD and FV)

# TODO: add tests
"""
function recover_profile(
    param_set::APS,
    sc::AbstractSurfaceConditions,
    L_MO,
    Z,
    X_in,
    X_sfc,
    transport,
    uft::UF.AUFT,
    scheme::Union{FVScheme, FDScheme},
) where {FT}
    @assert isless.(Z, sc.vals_in.z)
    uf = UF.universal_func(uft, L_MO, param_set)
    von_karman_const::FT = CPSGS.von_karman_const(param_set)
    _π_group = FT(UF.π_group(uf, transport))
    _π_group⁻¹ = (1 / _π_group)
    num1 = log(Z / z0(sc, transport))
    num2 = -UF.psi(uf, Z / L_MO, transport)
    num3 = UF.psi(uf, z0(sc, transport) / L_MO, transport)
    Σnum = num1 + num2 + num3
    ΔX = X_in - X_sfc
    return Σnum * compute_physical_scale_coeff(param_set, sc, L_MO, transport, uft, scheme) * _π_group⁻¹ * ΔX + X_sfc
end

end # SurfaceFluxes module


text:
@testset "indexing" begin
    L32 = @inferred(range(SInt32(1), stop=SInt32(4), length=SInt(4)))
    L64 = @inferred(range(SInt64(1), stop=SInt64(4), length=SInt(4)))
    @test @inferred(L32[1]) == 1.0 && @inferred(L64[1]) == 1.0
    @test L32[2] == 2 && L64[2] == 2
    @test L32[3] == 3 && L64[3] == 3
    @test L32[4] == 4 && L64[4] == 4

    @test @inferred(range(SVal(1.0), stop=SVal(2.0), length=SVal(2)))[1] == 1.0
    @test @inferred(range(SVal(1.0f0), stop=SVal(2.0f0), length=SVal(2)))[1] == 1.0f0
    @test @inferred(range(SVal(Float16(1.0)), stop=SVal(Float16(2.0)), length=SVal(2)))[1] == Float16(1.0)

    let r = SInt(5):SInt(-1):SInt(1)
        @test r[1]==5
        @test r[2]==4
        @test r[3]==3
        @test r[4]==2
        @test r[5]==1
    end
    # FIXME @test @inferred(range(SVal(0.1), step=SVal(0.1), stop=SVal(0.3))[2]) == 0.2
    @test @inferred(range(SVal(0.1f0), step=SVal(0.1f0), stop=SVal(0.3f0))[2]) == 0.2f0

    @test @inferred(SVal(1):SVal(5))[1:4] == 1:SVal(1):4
    #FIXME @test @inferred(SVal(1.0):SVal(5))[1:4] == 1.0:4
    @test (SVal(2):SVal(6))[1:4] == 2:5
    @test (SOne:SVal(6))[2:5] == 2:5
    @test (SOne:SVal(6))[2:2:5] == 2:2:4
    @test (SOne:SVal(2):SVal(13))[2:6] == 3:2:11
    @test (SOne:SVal(2):SVal(13))[2:3:7] == 3:6:13

    #@test isempty(srange(1:4)[5:4])
    #@test_throws BoundsError range(1:10)[8:-1:-2]

    let r = typemax(SInt)-SVal(5):typemax(SInt)-SOne
        @test_throws BoundsError r[7]
    end
    @testset "indexing range with empty range (#4309)" begin
        @test (SVal(3):SVal(6))[5:4] == 7:6
        @test_throws BoundsError (SVal(3):SVal(6))[5:5]
        @test_throws BoundsError (SVal(3):SVal(6))[5]
        @test (SZero:SVal(2):SVal(10))[7:6] == 12:2:10
        @test_throws BoundsError (SZero:SVal(2):SVal(10))[7:7]
    end
end

@testset "intersect" begin
    @test intersect(SOne:SVal(5), SVal(2):SVal(3)) == 2:3
    @test intersect(SVal(-3):SVal(5), SVal(2):SVal(8)) == 2:5
    @test intersect(SVal(-8):SVal(-3), SVal(-8):SVal(-3)) == -8:-3
    @test intersect(SOne:SVal(5), SVal(5):SVal(13)) == 5:5
    @test isempty(intersect(SVal(-8):SVal(-3), SVal(-2):SVal(2)))
    @test isempty(intersect(SVal(-3):SVal(7), SVal(2):SVal(1)))
    @test intersect(SOne:SVal(11), SVal(-2):SVal(3):SVal(15)) == 1:3:10
    @test intersect(SOne:SVal(11), SVal(-2):SVal(2):SVal(15)) == 2:2:10
    @test intersect(SOne:SVal(11), SVal(-2):SVal(1):SVal(15)) == 1:11
    @test intersect(SOne:SVal(11), SVal(15):SVal(-1):SVal(-2)) == 1:11
    @test intersect(SOne:SVal(11), SVal(15):SVal(-4):SVal(-2)) == 3:4:11
    @test intersect(SVal(-20):SVal(-5), SVal(-10):SVal(3):SVal(-2)) == -10:3:-7
    @test isempty(intersect(SVal(-5):SVal(5), SVal(-6):SVal(13):SVal(20)))
    @test isempty(intersect(SOne:SVal(11), SVal(15):SVal(4):SVal(-2)))
    @test isempty(intersect(SVal(11):SOne, SVal(15):SVal(-4):SVal(-2)))
    #@test intersect(-5:5, 1+0*(1:3)) == 1:1
    #@test isempty(intersect(-5:5, 6+0*(1:3)))
    @test intersect(SVal(-15):SVal(4):SVal(7), SVal(-10):SVal(-2)) == SVal(-7):SVal(4):SVal(-3)
    @test intersect(SVal(13):SVal(-2):SVal(1), SVal(-2):SVal(8)) == (SVal(7):SVal(-2):SVal(1))
    @test isempty(intersect(SVal(13):SVal(2):SOne, SVal(-2):SVal(8)))
    @test isempty(intersect(SVal(13):SVal(-2):SOne, SVal(8):SVal(-2)))
    #@test intersect(5+0*(1:4), 2:8) == 5+0*(1:4)
    #@test isempty(intersect(5+0*(1:4), -7:3))
    @test intersect(SZero:SVal(3):SVal(24), SZero:SVal(4):SVal(24)) == 0:12:24
    @test intersect(SZero:SVal(4):SVal(24), SZero:SVal(3):SVal(24)) == 0:12:24
    @test intersect(SZero:SVal(3):SVal(24), SVal(24):SVal(-4):SZero) == 0:12:24
    @test intersect(SVal(24):SVal(-3):SZero, SZero:SVal(4):SVal(24)) == 24:-12:0
    @test intersect(SVal(24):SVal(-3):SZero, SVal(24):SVal(-4):SZero) == 24:-12:0
    @test intersect(SVal(1):SVal(3):SVal(24), SZero:SVal(4):SVal(24)) == 4:12:16
    @test intersect(SZero:SVal(6):SVal(24), SZero:SVal(4):SVal(24)) == 0:12:24
    @test isempty(intersect(SOne:SVal(6):SVal(2400), SZero:SVal(4):SVal(2400)))
    @test intersect(SVal(-51):SVal(5):SVal(100), SVal(-33):SVal(7):SVal(125)) == -26:35:79
    @test intersect(SVal(-51):SVal(5):SVal(100), SVal(-32):SVal(7):SVal(125)) == -11:35:94
    #@test intersect(0:6:24, 6+0*(0:4:24)) == 6:6:6
    #@test intersect(12+0*(0:6:24), 0:4:24) == AbstractRange(12, 0, 5)
    #@test isempty(intersect(6+0*(0:6:24), 0:4:24))
    @test intersect(SVal(-10):SVal(3):SVal(24), SVal(-10):SVal(3):SVal(24)) == -10:3:23
    @test isempty(intersect(SVal(-11):SVal(3):SVal(24), SVal(-10):SVal(3):SVal(24)))

    #= TODO: this overflows when finding length, which happens in base too
    @test intersect(srange(typemin(Int):2:typemax(Int)), srange(1:10)) == srange(2:2:10)
    @test intersect(srange(1:10), srange(typemin(Int):2:typemax(Int))) == srange(2:2:10)

    @test intersect(reverse(srange(typemin(Int):2:typemax(Int))),srange(typemin(Int):2:typemax(Int))) == reverse(srange(typemin(Int):2:typemax(Int)))
    @test intersect(srange(typemin(Int):2:typemax(Int)), reverse(srange(typemin(Int):2:typemax(Int)))) == srange(typemin(Int):2:typemax(Int))
    =#

    @test intersect(SOne:SInt(2), SVal(3)) == 3:2
    @test intersect(SOne:SVal(2), SOne:SVal(5), SVal(3):SVal(7), SVal(4):SVal(6)) == 4:3

    @test intersect(SOne:SVal(3), SVal(2)) == intersect(SVal(2), SOne:SVal(3)) == 2:2
    @test intersect(SVal(1.0):SVal(3.0), SVal(2)) == intersect(SVal(2), SVal(1.0):SVal(3.0)) == [2.0]
end

@testset "findall(::Base.Fix2{typeof(in)}, ::Array)" begin
    @test findall(in(SVal(3):SVal(20)), [5.2, 3.3]) == findall(in(Vector(SInt(3):SInt(20))), [5.2, 3.3])

    let span = SVal(5):SVal(20),
        r = SVal(-7):SVal(3):SVal(42)
        @test findall(in(span), r) == 5:10
        r = SVal(15):SVal(-2):SVal(-38)
        @test findall(in(span), r) == 1:6
    end
end

@testset "findfirst" begin
    @test findfirst(isequal(7), SOne:SVal(2):SVal(10)) == 4
    @test findfirst(==(7), SOne:SVal(2):SVal(10)) == 4
    @test findfirst(==(10), SOne:SVal(2):SVal(10)) == nothing
    @test findfirst(==(11), SOne:SVal(2):SVal(10)) == nothing
end

@testset "reverse" begin
    @test reverse(reverse(SOne:SVal(10))) == 1:10
    @test reverse(reverse(typemin(SInt):typemax(SInt))) == typemin(Int):typemax(Int)
    @test reverse(reverse(typemin(SInt):SVal(2):typemax(SInt))) == typemin(Int):2:typemax(Int)
end

@testset "sort/sort!/partialsort" begin
    @test sort(SOne:SVal(2)) == 1:2
    @test sort!(SOne:SVal(2)) == 1:2
    @test sort(SOne:SVal(10), rev=true) == 10:-1:1
    @test sort(SVal(-3):SVal(3), by=abs) == [0,-1,1,-2,2,-3,3]
    @test partialsort(SOne:SVal(10), 4) == 4
end

@testset "length" begin
    @test length(SVal(.1):SVal(.1):SVal(.3)) == 3
    @test length(SVal(1.1):SVal(1.1):SVal(3.3)) == 3
    @test length(SVal(1.1):SVal(1.3):SVal(3)) == 2
    @test length(SVal(1):SVal(1):SVal(1.8)) == 1
    @test length(SVal(1):SVal(.2):SVal(2)) == 6
    @test length(SVal(1.):SVal(.2):SVal(2.)) == 6
    @test length(SVal(2):SVal(-.2):SVal(1)) == 6
    @test length(SVal(2.):SVal(-.2):SVal(1.)) == 6
    @test length(SVal(2):SVal(.2):(1)) == 0
    @test length(SVal(2.):SVal(.2):SVal(1.)) == 0

    @test length(SOne:SZero) == 0
    @test length(SVal(0.0):SVal(-0.5)) == 0
    @test length(SOne:SVal(2):SZero) == 0
# FIXME    @test length(srange(Char(0):Char(0x001fffff))) == 2097152
# FIXME    @test length(srange(typemax(SUInt64)//one(SUInt64):1:typemax(SUInt64)//one(SUInt64))) == 1
end

@testset "in" begin
    @test 0 in SUInt(0):SInt(100):typemax(SUInt)
    @test last(SVal(UInt(0)):SVal(100):SVal(typemax(UInt))) in SVal(UInt(0)):SVal(100):SVal(typemax(UInt))
    @test -9223372036854775790 in SVal(-9223372036854775790):SVal(100):SVal(9223372036854775710)
    @test -9223372036854775690 in SVal(-9223372036854775790):SVal(100):SVal(9223372036854775710)
    @test -90 in SVal(-9223372036854775790):SVal(100):SVal(9223372036854775710)
    @test 10 in SVal(-9223372036854775790):SVal(100):SVal(9223372036854775710)
    @test 110 in SVal(-9223372036854775790):SVal(100):SVal(9223372036854775710)
    @test 9223372036854775610 in SVal(-9223372036854775790):SVal(100):SVal(9223372036854775710)
    @test 9223372036854775710 in SVal(-9223372036854775790):SVal(100):SVal(9223372036854775710)

    @test !(SVal(3.5) in SOne:SVal(5))
    @test (3 in SOne:SVal(5))
    @test (3 in SVal(5):SVal(-1):SOne)
    #@test (3 in 3+0*(1:5))
    #@test !(4 in 3+0*(1:5))

    let r = SVal(0.0):SVal(0.01):SVal(1.0)
        @test (r[30] in r)
    end
    let r = (-SInt(4)*SInt64(maxintfloat(Int === Int32 ? Float32 : Float64))):SInt(5)
        @test (3 in r)
        @test (3.0 in r)
    end

    @test !(1 in SOne:SZero)
    @test !(1.0 in SVal(1.0):SVal(0.0))
end

# TODO: reimplement all below with SRational and SComplex
@testset "in() works across types, including non-numeric types (#21728)" begin
    @test 1//1 in SOne:SVal(3)
    @test 1//1 in SVal(1.0):SVal(3.0)
    @test !(5//1 in SOne:SVal(3))
    @test !(5//1 in SVal(1.0):SVal(3.0))

    @test Complex(1, 0) in SOne:SVal(3)
    @test Complex(1, 0) in SVal(1.0):SVal(3.0)
    @test Complex(1.0, 0.0) in SOne:SVal(3)
    @test Complex(1.0, 0.0) in SVal(1.0):SVal(3.0)
    @test !(Complex(1, 1) in SOne:SVal(3))
    @test !(Complex(1, 1) in SVal(1.0):SVal(3.0))
    @test !(Complex(1.0, 1.0) in SOne:SVal(3))
    @test !(Complex(1.0, 1.0) in SVal(1.0):SVal(3.0))
    @test !(π in SOne:SVal(3))
    @test !(π in SVal(1.0):SVal(3.0))
    @test !("a" in SOne:SVal(3))
    @test !("a" in SVal(1.0):SVal(3.0))

# FIXME   @test !(1 in srange(Date(2017, 01, 01):Dates.Day(1):Date(2017, 01, 05)))
# FIXME   @test !(Complex(1, 0) in srange(Date(2017, 01, 01):Dates.Day(1):Date(2017, 01, 05)))
# FIXME   @test !(π in srange(Date(2017, 01, 01):Dates.Day(1):Date(2017, 01, 05)))
# FIXME   @test !("a" in srange(Date(2017, 01, 01):Dates.Day(1):Date(2017, 01, 05)))
end


text:
module Trajectories

export AbstractPairedArray, Trajectory, trajectory, issynchron
export interpolate, Linear, Left, Right
export piecewise

include("definitions.jl")
include("interpolate.jl")
include("piecewise.jl")
include("tables.jl")
include("recipe.jl")

end # module