text:
module TestShow
using Test
using Dates
using CategoricalArrays

@testset "show() for CategoricalPool{String} and its values" begin
    pool = CategoricalPool(["c", "b", "a"])
    opool = CategoricalPool(["c", "b", "a"], true)

    nv1 = CategoricalValue(1, pool)
    nv2 = CategoricalValue(2, pool)
    nv3 = CategoricalValue(3, pool)

    ov1 = CategoricalValue(1, opool)
    ov2 = CategoricalValue(2, opool)
    ov3 = CategoricalValue(3, opool)

    @test sprint(show, pool) == "$CategoricalPool{String,UInt32}([\"c\", \"b\", \"a\"])"
    @test sprint(show, opool) == "$CategoricalPool{String,UInt32}([\"c\", \"b\", \"a\"]) with ordered levels"

    @test sprint(show, nv1) == "$CategoricalValue{String,UInt32} \"c\""
    @test sprint(show, nv2) == "$CategoricalValue{String,UInt32} \"b\""
    @test sprint(show, nv3) == "$CategoricalValue{String,UInt32} \"a\""

    @test sprint(show, ov1) == "$CategoricalValue{String,UInt32} \"c\" (1/3)"
    @test sprint(show, ov2) == "$CategoricalValue{String,UInt32} \"b\" (2/3)"
    @test sprint(show, ov3) == "$CategoricalValue{String,UInt32} \"a\" (3/3)"

    @test sprint(show, nv1, context=:typeinfo=>typeof(nv1)) == "\"c\""
    @test sprint(show, nv2, context=:typeinfo=>typeof(nv2)) == "\"b\""
    @test sprint(show, nv3, context=:typeinfo=>typeof(nv3)) == "\"a\""

    @test sprint(show, ov1, context=:typeinfo=>typeof(ov1)) == "\"c\""
    @test sprint(show, ov2, context=:typeinfo=>typeof(ov2)) == "\"b\""
    @test sprint(show, ov3, context=:typeinfo=>typeof(ov3)) == "\"a\""

    @test sprint(print, nv1) == sprint(print, ov1) == "c"
    @test sprint(print, nv2) == sprint(print, ov2) == "b"
    @test sprint(print, nv3) == sprint(print, ov3) == "a"

    @test string(nv1) == string(ov1) == "c"
    @test string(nv2) == string(ov2) == "b"
    @test string(nv3) == string(ov3) == "a"

    @test String(nv1) == String(ov1) == "c"
    @test String(nv2) == String(ov2) == "b"
    @test String(nv3) == String(ov3) == "a"

    @test repr(nv1) == repr(ov1) == "\"c\""
    @test repr(nv2) == repr(ov2) == "\"b\""
    @test repr(nv3) == repr(ov3) == "\"a\""

    b = IOBuffer()
    @test write(b, nv1) == 1
    @test String(take!(b)) == "c"
    @test write(b, nv2) == 1
    @test String(take!(b)) == "b"
    @test write(b, nv3) == 1
    @test String(take!(b)) == "a"
    @test write(b, ov1) == 1
    @test String(take!(b)) == "c"
    @test write(b, ov2) == 1
    @test String(take!(b)) == "b"
    @test write(b, ov3) == 1
    @test String(take!(b)) == "a"
end

@testset "show() for CategoricalPool{Date} and its values" begin
    levs = [Date(1999, 12), Date(1991, 8), Date(1993, 10)]
    pool = CategoricalPool(levs)
    opool = CategoricalPool(levs, true)

    nv1 = CategoricalValue(1, pool)
    nv2 = CategoricalValue(2, pool)
    nv3 = CategoricalValue(3, pool)

    ov1 = CategoricalValue(1, opool)
    ov2 = CategoricalValue(2, opool)
    ov3 = CategoricalValue(3, opool)

    if VERSION >= v"1.5.0-DEV"
        @test sprint(show, pool) == "$CategoricalPool{$Date,UInt32}($levs)"
        @test sprint(show, opool) == "$CategoricalPool{$Date,UInt32}($levs) with ordered levels"

        @test sprint(show, nv1) == "$CategoricalValue{$Date,UInt32} $Date(\"1999-12-01\")"
        @test sprint(show, nv2) == "$CategoricalValue{$Date,UInt32} $Date(\"1991-08-01\")"
        @test sprint(show, nv3) == "$CategoricalValue{$Date,UInt32} $Date(\"1993-10-01\")"

        @test sprint(show, ov1) == "$CategoricalValue{$Date,UInt32} $Date(\"1999-12-01\") (1/3)"
        @test sprint(show, ov2) == "$CategoricalValue{$Date,UInt32} $Date(\"1991-08-01\") (2/3)"
        @test sprint(show, ov3) == "$CategoricalValue{$Date,UInt32} $Date(\"1993-10-01\") (3/3)"

        @test sprint(show, nv1, context=:typeinfo=>typeof(nv1)) == "$Date(\"1999-12-01\")"
        @test sprint(show, nv2, context=:typeinfo=>typeof(nv2)) == "$Date(\"1991-08-01\")"
        @test sprint(show, nv3, context=:typeinfo=>typeof(nv3)) == "$Date(\"1993-10-01\")"

        @test sprint(show, ov1, context=:typeinfo=>typeof(ov1)) == "$Date(\"1999-12-01\")"
        @test sprint(show, ov2, context=:typeinfo=>typeof(ov2)) == "$Date(\"1991-08-01\")"
        @test sprint(show, ov3, context=:typeinfo=>typeof(ov3)) == "$Date(\"1993-10-01\")"
    else
        @test sprint(show, pool) == "$CategoricalPool{$Date,UInt32}([1999-12-01, 1991-08-01, 1993-10-01])"
        @test sprint(show, opool) == "$CategoricalPool{$Date,UInt32}([1999-12-01, 1991-08-01, 1993-10-01]) with ordered levels"

        @test sprint(show, nv1) == "$CategoricalValue{$Date,UInt32} 1999-12-01"
        @test sprint(show, nv2) == "$CategoricalValue{$Date,UInt32} 1991-08-01"
        @test sprint(show, nv3) == "$CategoricalValue{$Date,UInt32} 1993-10-01"

        @test sprint(show, ov1) == "$CategoricalValue{$Date,UInt32} 1999-12-01 (1/3)"
        @test sprint(show, ov2) == "$CategoricalValue{$Date,UInt32} 1991-08-01 (2/3)"
        @test sprint(show, ov3) == "$CategoricalValue{$Date,UInt32} 1993-10-01 (3/3)"

        @test sprint(show, nv1, context=:typeinfo=>typeof(nv1)) == "1999-12-01"
        @test sprint(show, nv2, context=:typeinfo=>typeof(nv2)) == "1991-08-01"
        @test sprint(show, nv3, context=:typeinfo=>typeof(nv3)) == "1993-10-01"

        @test sprint(show, ov1, context=:typeinfo=>typeof(ov1)) == "1999-12-01"
        @test sprint(show, ov2, context=:typeinfo=>typeof(ov2)) == "1991-08-01"
        @test sprint(show, ov3, context=:typeinfo=>typeof(ov3)) == "1993-10-01"
    end

    @test sprint(print, nv1) == sprint(print, ov1) == "1999-12-01"
    @test sprint(print, nv2) == sprint(print, ov2) == "1991-08-01"
    @test sprint(print, nv3) == sprint(print, ov3) == "1993-10-01"

    @test string(nv1) == string(ov1) == "1999-12-01"
    @test string(nv2) == string(ov2) == "1991-08-01"
    @test string(nv3) == string(ov3) == "1993-10-01"

    if VERSION >= v"1.5.0-DEV"
        @test repr(nv1) == repr(ov1) == "$Date(\"1999-12-01\")"
        @test repr(nv2) == repr(ov2) == "$Date(\"1991-08-01\")"
        @test repr(nv3) == repr(ov3) == "$Date(\"1993-10-01\")"
    else
        @test repr(nv1) == repr(ov1) == "1999-12-01"
        @test repr(nv2) == repr(ov2) == "1991-08-01"
        @test repr(nv3) == repr(ov3) == "1993-10-01"
    end
end

using JSON
@testset "JSON.lower" for pool in (CategoricalPool(["a"]),
                                   CategoricalPool([1]),
                                   CategoricalPool([1.0]))
    v = CategoricalValue(1, pool)
    @test JSON.lower(v) == JSON.lower(get(v))
    @test typeof(JSON.lower(v)) == typeof(JSON.lower(get(v)))
end

end


text:
# ---
# title: 45. Jump Game II
# id: problem45
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array, Greedy
# link: <https://leetcode.com/problems/jump-game-ii/description/>
# hidden: true
# ---
#
# Given an array of non-negative integers `nums`, you are initially positioned
# at the first index of the array.
#
# Each element in the array represents your maximum jump length at that
# position.
#
# Your goal is to reach the last index in the minimum number of jumps.
#
# You can assume that you can always reach the last index.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [2,3,1,1,4]
#     Output: 2
#     Explanation: The minimum number of jumps to reach the last index is 2. Jump 1 step from index 0 to 1, then 3 steps to the last index.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [2,3,0,1,4]
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 3 * 104`
#   * `0 <= nums[i] <= 105`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
#=
Provides a function to solve for asset prices under optimistic beliefs in the
Harrison -- Kreps model.

Authors: Shunsuke Hori
=#
"""
Function to Solve Optimistic Beliefs
"""

function price_optimisticbeliefs(transitions, dividend_payoff;
                                beta=.75,max_iter=50000, tol=1e-16)

    # We will guess an initial price vector of [0, 0]
    p_new = [0,0]
    p_old = [10.0,10.0]

    # We know this is a contraction mapping, so we can iterate to conv
    for i in 1:max_iter
        p_old = p_new
        temp=[maximum((q* p_old) + (q* dividend_payoff))
                               for q in transitions]
        p_new = beta * temp

        # If we succed in converging, break out of for loop
        if maximum(sqrt((p_new - p_old).^2)) < 1e-12
            break
        end
    end

    temp=[minimum((q* p_old) + (q* dividend_payoff))
                               for q in transitions]
    ptwiddle = beta * temp

    phat_a = [p_new[1], ptwiddle[2]]
    phat_b = [ptwiddle[1], p_new[2]]

    return p_new, phat_a, phat_b
end


text:
include("./util.jl")
using Test: @test

box_specs = util.load_input(2)

struct Box
    l::Int
    w::Int
    h::Int
end

function Box(s::AbstractString)::Box
    l, w, h = tuple(map(part -> parse(Int, part), split(s, 'x'))...)
    Box(l, w, h)
end

surface_area(b::Box)::Int = 2 * b.l * b.w + 2 * b.w * b.h + 2 * b.h * b.l

paper_area(b::Box)::Int = surface_area(b) + min(b.l * b.w, b.w * b.h, b.h * b.l)

@test paper_area(Box("2x3x4")) == 58
@test paper_area(Box("1x1x10")) == 43

perimeters(b::Box)::Tuple{Int,Int,Int} = (2 * (b.l + b.w), 2 * (b.w + b.h), 2 * (b.h + b.l))

ribbon_length(b::Box)::Int = min(perimeters(b)...) + b.l * b.w * b.h

boxes = map(Box, split(box_specs, '\n'))

total_sqft = sum(map(paper_area, boxes))
println("Part 1: $total_sqft")

total_ribbon_ft = sum(map(ribbon_length, boxes))
println("Part 2: $total_ribbon_ft")


text:
struct TauLeapingConstantCache <: StochasticDiffEqConstantCache end

@cache struct TauLeapingCache{uType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  newrate::rateType
  EEstcache::rateType
end

alg_cache(alg::TauLeaping,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) = TauLeapingConstantCache()

function alg_cache(alg::TauLeaping,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  tmp = zero(u)
  newrate = zero(jump_rate_prototype)
  EEstcache = zero(jump_rate_prototype)
  TauLeapingCache(u,uprev,tmp,newrate,EEstcache)
end

alg_cache(alg::CaoTauLeaping,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) = TauLeapingConstantCache()

function alg_cache(alg::CaoTauLeaping,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  tmp = zero(u)
  TauLeapingCache(u,uprev,tmp,nothing,nothing)
end


text:
## update inputs vars using calculated gradients

function getfield_nested(s, path)
    v = s
    for name in path
        v = getfield(v, name)
    end
    return v
end

function setfield_nested!(s, path, val)
    container = s
    for name in path[1:end-1]
        container = getfield(container, name)
    end
    setfield!(container, path[end], val)
end

minus_updater(x::AbstractArray, gx) = (x .- gx)
minus_updater(x::Real, gx) = (x - gx)

function update!(x::AbstractArray{T,N}, gx::AbstractArray{T,N}, fn::Function=minus_updater) where {T,N}
    x .= fn(x, gx)
end

function update!(x::Real, gx::Real, fn::Function=minus_updater)
    error("Can't update value of $(typeof(x)) in-place!")
end

function update!(x, gx, fn::Function=minus_updater)
    @assert isstruct(x) "Expected mutable struct as 1st argument"
    @assert gx isa Dict "Gradient for structs should be a dict of (field path -> grad value)"
    for (path, gv) in gx
        v = getfield_nested(x, path)
        setfield_nested!(x, path, fn(v, gv))
    end
end


text:
using FileIO
using CSVFiles
using DataFrames

function load_datatable(file_path::AbstractString; skip::Int = 0)
    return load_DataFrame(file_path;
        skip_lines = skip, delimiter = space::Delimiter, header = true)
end

function load_DataFrame(
    file_path::AbstractString;
    skip_lines::Int = 0,
    delimiter::Delimiter = space::Delimiter,
    header::Bool = true, column_names::Dict{Int, AbstractString} = Dict{Int, AbstractString}())
    @assert(isfile(file_path), ArgumentError("The file '$file_path' does not exist or is not a file."))
    file_object = query(file_path)
    @assert(!unknown(file_object), ArgumentError("The file '$file_path' has unknown extension and/or format."))
    file_format = typeof(file_object).parameters[1]
    @info "Detected the format $file_format for the file '$file_path'. Loading to DataFrame..."

    data = if isempty(column_names)
        load(file_object;
            skiplines_begin = skip_lines,
            spacedelim=(delimiter == space::Delimiter),
            header_exists=header,
            type_detect_rows=22) |> DataFrame
    else
        load(file_object;
            skiplines_begin = skip_lines,
            spacedelim=(delimiter == space::Delimiter),
            header_exists=header, colnames=column_names,
            type_detect_rows=22) |> DataFrame
    end
    if data == nothing
        @warn "... No data has been loaded."
        return nothing
    end
    @info "... Data loaded to a new DataFrame with size $(size(data, 1))×$(size(data, 2))."
    return data
end


text:
function LeastSquares(name::String, wn_dict::Dict{Any,Any}, wn::PyCall.PyObject, link_results::Dict{Any,Any}, node_results::Dict{Any,Any}, num_timesteps::Int64)
    energy = metric.pump_energy(link_results["flowrate"], node_results["head"], wn)[name][:values][1:num_timesteps]
    flow = link_results["flowrate"][name][:values][1:num_timesteps]
    flows = Array{Float64}(0)
    energies = Array{Float64}(0)
    for i= 1:length(energy)
        if energy[i] != 0 && flow[i] !=0
            flows = vcat(flows,flow[i])
            energies = vcat(energies, energy[i])
        end
    end
    #if pump never turns on, delete control that inhibits pump, calulation energy, then replace control
    if length(energies) == 0
        # warn("$name does not turn on. Ajusting controls to allow for turn on of pump.")
        controls_dict = wn_dict["controls"]
        control_info = Array{PyCall.PyObject}(0)
        control_names = Array{String}(0)
        for control in keys(controls_dict)
            push!(control_info, controls_dict[control])
            push!(control_names, control)
            wn[:remove_control](control)
        end
        wns = sim.EpanetSimulator(wn) #Simulation either Demand Driven(DD) or Presure Dependent Demand(PDD), default DD
        results = wns[:run_sim]()
        node_results = results[:node] #dictionary of head, pressure, demand, quality
        link_results = results[:link]
        energy = metric.pump_energy(link_results["flowrate"], node_results["head"], wn)[name][:values][1:num_timesteps]
        flow = link_results["flowrate"][name][:values][1:num_timesteps]
        flows = Array{Float64}(0)
        energies = Array{Float64}(0)
        for i= 1:length(energy)
            if energy[i] != 0 && flow[i] !=0
                flows = vcat(flows,flow[i])
                energies = vcat(energies, energy[i])
            end
        end
        for (ix, control) in enumerate(control_names)
            wn[:add_control](control, control_info[ix])
        end
    end

    avg_energy = mean(energies)
    a, b = linear_fit(flows, energies)
    error = Array{Float64}(0)
    for i = 1:length(energies)
        push!(error, energies[i] - (avg_energy +b *flows[i]))
    end
    avg_error = mean(error)
    intercept = avg_energy + avg_error
    return intercept, b
    # total_error = energies - (intercept + b * flows)

end


text:
# Julia NRLMSISE00 implementation is Dominik Brodowski's C implmenation adapted
# for Julia.
#
# Ref: https://www.brodo.de/space/nrlmsise/

##############
# Model Data #
##############

include("nrlmsise00_data.jl")

###################
# Model Structure #
###################

"""
Internal data structure used for setting model parameters
"""
mutable struct NRLMSISE_Flags
    switches::Array{<:Int, 1}
    sw::Array{<:Real, 1}
    swc::Array{<:Real, 1}

    function NRLMSISE_Flags(switches::Array{<:Int, 1}=zeros(Int, 24),
            sw::Array{<:Real, 1}=zeros(Float32, 24),
            swc::Array{<:Real, 1}=zeros(Float32, 24),)
        new(switches, sw, swc)
    end
end

"""
Internal structure for storing model input

ap_array contains the following magnetic values:
- 1 : daily AP
- 2 : 3 hr AP index for current time
- 3 : 3 hr AP index for 3 hrs before current time
- 4 : 3 hr AP index for 6 hrs before current time
- 5 : 3 hr AP index for 9 hrs before current time
- 6 : Average of eight 3 hr AP indices from 12 to 33 hrs prior to current time
- 7 : Average of eight 3 hr AP indices from 36 to 57 hrs prior to current time
"""
mutable struct NRLMSISE_Input
    year::Int       # Year, currently ignored
    doy::Int        # Day of year
    sec::Float64    # Seconds in day (UT)
    alt::Float64    # Altitude in [km]
    g_lat::Float64  # Geodetic latitude
    g_lon::Float64  # Geodetic longitude
    lst::Float64    # Local apparent solar time (hours)
    f107A::Float64  # 81 day average of F10.7cm flux (centered on day)
    f107::Float64   # Daily F10.7cm flux for previous day
    ap::Float64     # Magnetic index (daily)
    ap_array::Array{<:Real, 1} # Magnetic index array

    function NRLMSISE_Input(year::Int=2000, doy::Int=1, sec::Float64=0.0,
                alt::Float64=0.0, g_lat::Float64=0.0, g_lon::Float64=0.0,
                lst::Float64=0.0, f107A::Float64=0.0, f107::Float64=0.0,
                ap::Float64=0.0, ap_array::Array{<:Real, 1}=zeros(Float64, 7))
        new(year, doy, sec, alt, g_lat, g_lon, lst, f107A, f107, ap, ap_array)
    end
end

"""
# Output variables
d = zeros(Float64, 9)
t = zeros(Float64, 2)
"""
mutable struct NRLMSISE_Output
    d::Array{Float64, 1}
    t::Array{Float64, 1}

    function NRLMSISE_Output(d::Array{<:Real, 1}=zeros(Float64, 9), t::Array{<:Real, 1}=zeros(Float64, 2))
        new(d, t)
    end
end

#####################
# Utility Functions #
#####################

"""
Select internal model flags
"""
function tselec!(flags::NRLMSISE_Flags)
    for i in 1:24
        if i != 10
            flags.switches[i] == 1 ? flags.sw[i] = 1 : flags.sw[i] = 0
            flags.switches[i] >  0 ? flags.swc[i] = 1 : flags.swc[i] = 0

        else
            flags.sw[i]  = flags.switches[i]
            flags.swc[i] = flags.switches[i]
        end
    end
end

function glatf(lat::Real)
    dgtr = 1.74533E-2

    c2  = cos(2.0*dgtr*lat)

    gv   = 980.616 * (1.0 - 0.0026373 * c2)
    reff = 2.0 * (gv) / (3.085462E-6 + 2.27E-9 * c2) * 1.0E-5

    return gv, reff
end

"""
Chemistry/Dissociation correction for MSIS models

Arguments:
- `alt::Real` Altitude
- `r::Real` Target ratio
- `h1::Real` Transition scale length
- `zh::Real` Altitude of 1/2 R

Returns:
- `e::Real` Correction coefficient
"""
function ccor(alt::Real, r::Real, h1::Real, zh::Real)
    e = (alt - zh)/h1

    if e > 70
        return exp(0)
    elseif e < -70
        return exp(r)
    end

    ex = exp(e)
    e  = r/(1.0+ex)
    return exp(e)
end

"""
Chemistry/Dissociation correction for MSIS models

Arguments:
- `alt::Real` Altitude
- `r::Real` Target ratio
- `h1::Real` Transition scale length
- `zh::Real` Altitude of 1/2 R
- `h2::Real` Transition scale length 2

Returns:
- `e::Real` Correction coefficient
"""
function ccor2(alt::Real, r::Real, h1::Real, zh::Real, h2::Real)
    e1 = (alt - zh)/h1
    e2 = (alt - zh)/h2

    if (e1 > 70) || (e2 > 70)
        return exp(0)
    elseif (e1 < -70) && (e2 < -70) # TODO: Should this be && not ||?
        return exp(r)
    end

    ex1 = exp(e1)
    ex2 = exp(e2)

    ccor2v = r/(1.0+0.5*(ex1+ex2))
    return exp(ccor2v)
end

"""
Compute scale height
"""
function scaleh(alt::Real, xm::Real, temp::Real; gsurf::Real, re::Real)
    rgas = 831.4
    g = gsurf/(1.0 + alt/re)^2
    g = rgas*temp/(g*xm)

    return g
end

"""
Turbopause correction for MSIS models

Arguments:
- `dd::Real` diffusive density
- `dm::Real` full mixed density density
- `zhm::Real` transition scale length
- `xmm::Real` full mixed molecular weight
- `xm::Real` species molecular weight
- `dnet::Real` combined density
"""
function dnet(dd::Real, dm::Real, zhm::Real, xmm::Real, xm::Real)
    a = zhm/(xmm-xm)

    if !((dm>0) && (dd>0))
        @error "dnet log error $dm $dd $xm"
        if (dd==0) && (dm==0)
            dd = 1
        end

        if dm == 0
            return dd
        end

        if dd == 0
            return dm
        end
    end

    ylog = a * log(dm/dd)

    if ylog < -10
        return dd
    end

    if ylog > 10
        return dm
    end

    a = dd*(1.0 + exp(ylog))^(1.0/a)

    return a
end

"""
Integrate cubic spline from xa[1] to x

Arguments:
-`xa::Array{Real, 1}` Array of tabulated function inputs in ascending order by x
-`ya::Array{Real, 1}` Array of tabulated function outputs in ascending order by x
-`y2a::Array{Real, 1}` Array of second derivatives
-`x::Real` Array output point

Returns:
-`y::Real` Array output value
"""
function splini(xa::Array{<:Real, 1}, ya::Array{<:Real, 1}, y2a::Array{<:Real, 1}, n::Int, x::Real)
    yi  = 0
    klo = 1
    khi = 2

    while (x > xa[klo]) && (khi <= n)
        xx = x
        if (khi < n)
            if x < xa[khi]
                xx = x
            else
                xx = xa[khi]
            end
        end
        h    = xa[khi] - xa[klo]
        a    = (xa[khi] - xx)/h
        b    = (xx - xa[klo])/h
        a2   = a^2
        b2   = b^2
        yi  += ((1.0 - a2) * ya[klo] / 2.0 + b2 * ya[khi] / 2.0 + ((-(1.0+a2*a2)/4.0 + a2/2.0) * y2a[klo] + (b2*b2/4.0 - b2/2.0) * y2a[khi]) * h * h / 6.0) * h
        klo += 1
        khi += 1
    end

    return yi
end

"""
Interpolate cubic spline from to x

Arguments:
-`xa::Array{Real, 1}` Array of tabulated function inputs in ascending order by x
-`ya::Array{Real, 1}` Array of tabulated function outputs in ascending order by x
-`y2a::Array{Real, 1}` Array of second derivates
-`x::Real` Array output point

Returns:
-`y::Real` Array output value
"""
function splint(xa::Array{<:Real, 1}, ya::Array{<:Real, 1}, y2a::Array{<:Real, 1}, n::Int, x::Real)
    klo = 1
    khi = n
    k   = 0
    while (khi - klo) > 1
        k = floor(Int, (khi + klo)/2)
        if xa[k] > x
            khi = k
        else
            klo = k
        end
    end

    h = xa[khi] - xa[klo]
    if h == 0
        @error "Bad array input to splint. xa[$khi]: $(xa[khi]), xa[$klo]: $(xa[klo])"
        @error "Array Input: $xa, $ya"
    end

    a = (xa[khi] - x)/h
    b = (x - xa[klo])/h

    y = a * ya[klo] + b * ya[khi] + ((a*a*a - a) * y2a[klo] + (b*b*b - b) * y2a[khi]) * h * h/6.0

    return y
end

"""
Interpolate cubic spline from to x

Arguments:
-`x::Array{Real, 1}` Array of tabulated function inputs in ascending order by x
-`y::Array{Real, 1}` Array of tabulated function outputs in ascending order by x
-`yp1::Real` Specified derivative at x[1]   (yp1 >= 1e30 signals second derivative is zero)
-`ypn::Real` Specified derivative at x[end] (ypn >= 1e30 signals second derivative is zero)
-`x::Real` Array output point

Returns:
-`y::Array{Float64, 1}` Array second derivatives
"""
function spline(x::Array{<:Real, 1}, y::Array{<:Real, 1}, n::Int, yp1::Real, ypn::Real)
    u  = zeros(Float64, n) # Need the + 1 here because C implementation already allocates extra
    y2 = zeros(Float64, n)

    if yp1 > 0.99e30
        y2[1] = 0
        u[1]  = 0
    else
        y2[1] = -0.5
        u[1]  = (3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1)
    end

    # println("u: $u")
    # println("y2: $y2")

    for i in 2:n-1
        sig = (x[i]-x[i-1])/(x[i+1] - x[i-1])
        p     = sig*y2[i-1] + 2.0
        y2[i] = (sig - 1.0)/p
        u[i]  = (6.0 * ((y[i+1] - y[i])/(x[i+1] - x[i]) -(y[i] - y[i-1]) / (x[i] - x[i-1]))/(x[i+1] - x[i-1]) - sig * u[i-1])/p
    end

    # println("u: $u")
    # println("y2: $y2")

	if ypn>0.99E30
		qn = 0
		un = 0
    else
		qn = 0.5
		un = (3.0 / (x[n] - x[n-1])) * (ypn - (y[n] - y[n-1])/(x[n] - x[n-1]))
    end

    # println("un: $un")

    y2[n] = (un - qn * u[n-1]) / (qn * y2[n-1] + 1.0)

    # println("y2: $y2")

    for k in (n-1):-1:1
        y2[k] = y2[k] * y2[k+1] + u[k]
    end

    # println("y2: $y2")

    return y2
end

@inline function zeta(zz::Real, zl::Real; re::Real)
    return ((zz-zl)*(re+zl)/(re+zz))
end

"""
Calculate temperature and density profiles for lower atmosphere

Arguments:
- `alt::Real`
- `dm28m`
- `xnm`
- `mn3`
- `zn3`
- `mn2`
- `zn2`

Returns:
- `n2` N2 number density
- `tz`
- `meso_tn3`
- `meso_tgn3`
- `meso_tn2`
- `meso_tgn2`
"""
function densm!(alt::Real, d0::Real, xm::Real, tz::Real,
                mn3::Int, zn3::Array{<:Real, 1}, tn3::Array{<:Real, 1}, tgn3::Array{<:Real, 1},
                mn2::Int, zn2::Array{<:Real, 1}, tn2::Array{<:Real, 1}, tgn2::Array{<:Real, 1};
                gsurf::Real, re::Real)

    # @printf "alt: %.10f xm: %f tz: %.10f\n" alt xm tz

    xs = zeros(Float64, 10)
    ys = zeros(Float64, 10)

    rgas = 831.4

    densm_tmp = d0

    if alt > zn2[1]
        if xm == 0.0
            return tz, tz
        else
            # println("densm return 1")
            return tz, d0
        end
    end

    # Straosphere/Mesosphere Temperature
    if alt > zn2[mn2]
        z = alt
    else
        z = zn2[mn2]
    end
    mn = mn2
    z1 = zn2[1]
    z2 = zn2[mn]
    t1 = tn2[1]
    t2 = tn2[mn]
    zg = zeta(z, z1, re=re)
    zgdif = zeta(z2, z1, re=re)

    # Setup spline nodes
    for k in 1:mn
        xs[k] = zeta(zn2[k], z1, re=re)/zgdif
        ys[k] = 1.0/tn2[k]
    end
    yd1 = -tgn2[1] / (t1*t1) * zgdif
    yd2 = -tgn2[2] / (t2*t2) * zgdif * ((re+z2)/(re+z1))^2

    # Calculate spline coefficients
    y2out = spline(xs, ys, mn, yd1, yd2)
    x     = zg/zgdif
    y     = splint(xs, ys, y2out, mn, x)

    # Temperature at altitude
    tz = 1.0 / y

    # @printf "tz: %.10f densm_tmp: %.10f\n" tz densm_tmp

    if xm != 0
        # Calculate stratosphere/mesosphere density
        glb  = gsurf / (1.0 + z1/re)^2
        gamm = xm * glb * zgdif / rgas

        # Integrate temperature profiles
        yi = splini(xs, ys,y2out, mn, x)
        expl = gamm*yi
        if expl > 50
            expl = 50
        end

        # Density at altitude
        densm_tmp = densm_tmp * (t1 / tz) * exp(-expl)
    end

    if alt > zn3[1]
        if xm == 0.0
            return tz, tz
        else
            # println("densm return 2")
            # @printf "tz: %.10f densm_tmp: %.10f\n" tz densm_tmp
            return tz, densm_tmp
        end
    end

    # Troposphere / stratosphere temperatures
    z  = alt
    mn = mn3
    z1 = zn3[1]
    z2 = zn3[mn]
    t1 = tn3[1]
    t2 = tn3[mn]
    zg = zeta(z, z1, re=re)
    zgdif = zeta(z2, z1, re=re)

    # @printf "mn: %d z1: %f z2: %f t1: %f t2: %f zg: %f zgdif: %f\n" mn z1 z2 t1 t2 zg zgdif

    # Setup spline nodes
    for k in 1:mn
        xs[k] = zeta(zn3[k], z1, re=re)/zgdif
        ys[k] = 1.0/tn3[k]
    end
    yd1 = -tgn3[1] / (t1*t1) * zgdif
    yd2 = -tgn3[2] / (t2*t2) * zgdif * ((re+z2)/(re+z1))^2

    # Calculate spline coefficients
    # println("mn: $mn")
    # println("xs: $xs")
    # println("ys: $ys")
    # println("yd1: $yd1")
    # println("yd2: $yd2")
    y2out = spline(xs, ys, mn, yd1, yd2)
    # println("y2out: $y2out")
    x = zg/zgdif
    y = splint(xs, ys, y2out, mn, x)
    # println("x: $x")
    # println("y: $y")

    # @printf "tz: %.10f densm_tmp: %.10f\n" tz densm_tmp

    # Temperature at altitude
    tz = 1.0 / y

    # @printf "tz: %.10f densm_tmp: %.10f\n" tz densm_tmp

    if xm != 0
        # Calculate stratosphere/mesosphere density
        glb  = gsurf / (1.0 + z1/re)^2
        gamm = xm * glb * zgdif / rgas

        # Integrate temperature profiles
        yi = splini(xs, ys, y2out, mn, x)
        expl = gamm*yi
        if expl > 50
            expl = 50
        end

        # Density at altitude
        densm_tmp = densm_tmp * (t1 / tz) * exp(-expl)
    end

    if xm == 0.0
        return tz, tz
    else
        # println("densm return 3")
        # @printf "tz: %.10f densm_tmp: %.10f\n" tz densm_tmp
        return tz, densm_tmp
    end
end

function densu!(alt::Real, dlb::Real, tinf::Real, tlb::Real, xm::Real,
               alpha::Real, tz::Real, zlb::Real, s2::Real,
               mn1::Int, zn1::Array{<:Real, 1}, tn1::Array{<:Real, 1}, tgn1::Array{<:Real, 1};
               gsurf::Real, re::Real)

    x          = 0
    rgas       = 831.4
    densu_temp = 1.0
    z1         = 0
    t1         = 0
    zgdif      = 0

    xs    = zeros(Float64, 5)
    ys    = zeros(Float64, 5)
    y2out = zeros(Float64, 5)

    # Join altitudes of Bates and spline
    za = zn1[1]
    if alt > za
        z = alt
    else
        z = za
    end

    # println("z: $z, zlb: $zlb, tinf: $tinf, tlb: $tlb")

    # Geopotential altitude difference from ZLB
    zg2 = zeta(z, zlb, re=re)

    # Bates temperatures
    tt = tinf - (tinf - tlb) * exp(-s2*zg2)
    ta = tt
    tz = tt
    densu_temp = tz

    # println("alt: $alt, za: $za")
    if alt < za
        # Calculate temperature below za
        # Get temperaturre gradient at za from Bates profile
        dta = (tinf - ta) * s2 * ((re+zlb)/(re+za))^2

        tgn1[1] = dta
        tn1[1]  = ta

        if alt > zn1[end]
            z = alt
        else
            z = zn1[end]
        end

        mn = mn1
        z1 = zn1[1]
        z2 = zn1[end]
        t1 = tn1[1]
        t2 = tn1[end]

        # Geopotential difference from z1
        zg    = zeta(z, z1, re=re)
        zgdif = zeta(z2, z1, re=re)

        # Setup spline nodes
        for k in 1:mn
            xs[k] = zeta(zn1[k], z1, re=re)/zgdif
            ys[k] = 1.0/tn1[k]
        end

        # End node derivatives
        yd1 = -tgn1[1] / (t1*t1) * zgdif
        yd2 = -tgn1[2] / (t2*t2) * zgdif * ((re+z2)/(re+z1))^2

        # Calculate spline coefficients
        y2out[5] = 0.0
        # println("xs: $xs")
        # println("ys: $ys")
        # println("mn: $mn")
        # println("yd1: $yd1")
        # println("yd2: $yd2")
        # println("before spline call")
        y2out = spline(xs, ys, mn, yd1, yd2)
        x     = zg/zgdif
        y     = splint(xs, ys, y2out, mn, x)

        # Temperature at altitude
        tz = 1.0 / y
        densu_temp = tz
        # println("tz: $tz, densu_temp: $densu_temp")
    end

    # println("break 1")
    # println("tz: $tz, densu_temp: $densu_temp")

    if xm == 0
        # println("return 1")
        # println("tz: $tz, densu_temp: $densu_temp")
        return tz, densu_temp
    end

    # Calculate density above za
    glb = gsurf / (1.0+zlb/re)^2
    gamma = xm *glb / (s2 * rgas * tinf)
    expl = exp(-s2 * gamma * zg2)

    if expl > 50
        expl = 50
    end

    if tt <= 0
        expl = 50
    end

    # println("break 2")
    # println("tz: $tz, densu_temp: $densu_temp")

    # @printf "dlb: %.10f tlb: %.10f tt: %.10f alpha: %.10f gamma: %.10f expl: %.10f\n" dlb tlb tt alpha gamma expl

    # density at altitude
    densa      = dlb * (tlb/tt)^(1.0+alpha+gamma)*expl
    densu_temp = densa

    # println("break 3")
    # println("tz: $tz, densu_temp: $densu_temp")

    if alt >= za
        # println("return 2")
        # println("tz: $tz, densu_temp: $densu_temp")
        return tz, densu_temp
    end

    # println("branch here")
    # println("tz: $tz, densu_temp: $densu_temp")

    # Calculate density below za
    glb  = gsurf / (1.0 + z1/re)^2.0
    gamm = xm * glb * zgdif / rgas

    # Integrate spline temperatures
    yi = splini(xs, ys, y2out, mn, x)
    # println("xs: $xs")
    # println("ys: $ys")
    # println("y2out: $y2out")
    # println("mn: $mn")
    # println("x: $x")

    expl = gamm * yi
    # @printf "%.10f\n" yi
    # println("glb: $glb, gamm: $gamm, yi: $yi, expl: $expl")

    if expl > 50.0
        expl = 50.0
    end

    if tz <= 0
        expl = 50.0
    end

    # density at altitude
    # println("densu_temp: $densu_temp, t1: $t1, tz: $tz, alpha: $alpha, expl: $expl")
    densu_temp = densu_temp * (t1/tz)^(1.0+alpha) * exp(-expl)
    # println("return 3")
    # println("tz: $tz, densu_temp: $densu_temp")
    return tz, densu_temp
end

# Equation A24a
@inline function g0(a::Real, p::Array{<:Real, 1})
    return (a - 4.0 + (p[26] - 1.0) * (a - 4.0 + (exp(-sqrt(p[25]*p[25]) * (a - 4.0)) - 1.0) / sqrt(p[25]*p[25])))
end

# Equation A24c
@inline function sumex(ex::Real)
    return (1.0 + (1.0 - ex^19.0) / (1.0 - ex) * ex^0.5)
end

# Equation A24a
@inline function sg0(ex::Real, p::Array{<:Real, 1}, ap::Array{<:Real, 1})
    return (g0(ap[2],p) + (g0(ap[3],p)*ex + g0(ap[4],p)*ex*ex +
            g0(ap[5],p)*ex^3.0	+ (g0(ap[6],p)*ex^4.0 +
            g0(ap[7],p)*ex^12.0)*(1.0- ex^8.0)/(1.0-ex)))/sumex(ex)
end

##################
# NRLMSISE Model #
##################

"""
Calculate G(L) function
"""
function globe7(p::Array{<:Real, 1}, input::NRLMSISE_Input, flags::NRLMSISE_Flags)
    # Working variables
    t = zeros(Float64, 15)

    # LPOLY Variables
    dfa = 0.0
    plg = Array{Float64, 1}[zeros(Float64, 9) for i in 1:4]
    ctloc,  stloc  = 0.0, 0.0
    c2tloc, s2tloc = 0.0, 0.0
    s3tloc, c3tloc = 0.0, 0.0
    apdf = 0.0
    apt  = zeros(Float64, 4)

    # Upper therrmosphere parameters
    sr   = 7.2722E-5
    dgtr = 1.74533E-2
    dr   = 1.72142E-2
    hr   = 0.2618

    tloc = input.lst
    for j in 1:14
        t[j] = 0
    end

    # Calculate legendre polynomials
    c  = sin(input.g_lat * dgtr)
    s  = cos(input.g_lat * dgtr)
    c2 = c*c
    c4 = c2*c2
    s2 = s*s

    plg[1][2] = c
	plg[1][3] = 0.5*(3.0*c2 -1.0)
	plg[1][4] = 0.5*(5.0*c*c2-3.0*c)
	plg[1][5] = (35.0*c4 - 30.0*c2 + 3.0)/8.0
	plg[1][6] = (63.0*c2*c2*c - 70.0*c2*c + 15.0*c)/8.0
	plg[1][7] = (11.0*c*plg[1][6] - 5.0*plg[1][5])/6.0
    # plg[1][8] = (13.0*c*plg[1][7] - 6.0*plg[1][6])/7.0
	plg[2][2] = s
	plg[2][3] = 3.0*c*s
	plg[2][4] = 1.5*(5.0*c2-1.0)*s
	plg[2][5] = 2.5*(7.0*c2*c-3.0*c)*s
	plg[2][6] = 1.875*(21.0*c4 - 14.0*c2 +1.0)*s
	plg[2][7] = (11.0*c*plg[2][6]-6.0*plg[2][5])/5.0
    # plg[2][8] = (13.0*c*plg[2][7]-7.0*plg[2][6])/6.0
    # plg[2][9] = (15.0*c*plg[2][8]-8.0*plg[2][7])/7.0
	plg[3][3] = 3.0*s2
	plg[3][4] = 15.0*s2*c
	plg[3][5] = 7.5*(7.0*c2 -1.0)*s2
	plg[3][6] = 3.0*c*plg[3][5]-2.0*plg[3][4]
	plg[3][7] = (11.0*c*plg[3][6]-7.0*plg[3][5])/4.0
	plg[3][8] = (13.0*c*plg[3][7]-8.0*plg[3][6])/5.0
	plg[4][4] = 15.0*s2*s
	plg[4][5] = 105.0*s2*s*c
	plg[4][6] = (9.0*c*plg[4][5]-7.0*plg[4][4])/2.0
	plg[4][7] = (11.0*c*plg[4][6]-8.0*plg[4][5])/3.0

    if !(((flags.sw[8] == 0) && (flags.sw[9] == 0)) && (flags.sw[15] == 0))
        stloc  = sin(hr*tloc)
        ctloc  = cos(hr*tloc)
        s2tloc = sin(2.0*hr*tloc)
        c2tloc = cos(2.0*hr*tloc)
        s3tloc = sin(3.0*hr*tloc)
        c3tloc = cos(3.0*hr*tloc)
    end

    cd32 = cos(dr*(input.doy - p[32]))
    cd18 = cos(2.0*dr*(input.doy - p[18]))
    cd14 = cos(dr*(input.doy - p[14]))
    cd39 = cos(2.0*dr*(input.doy - p[39]))

    # F10.7 Effect
    df   = input.f107 - input.f107A
    dfa = input.f107A - 150.0
    t[1] = p[20]*df*(1.0+p[60]*dfa) + p[21]*df*df + p[22]*dfa + p[30]*dfa^2.0
    f1   = 1.0 + (p[48]*dfa +p[20]*df+p[21]*df*df)*flags.swc[2]
    f2   = 1.0 + (p[50]*dfa+p[20]*df+p[21]*df*df)*flags.swc[2]

    # Time independent
    t[2] = (p[2]*plg[1][3]+ p[3]*plg[1][5]+p[23]*plg[1][7]) +
           (p[15]*plg[1][3])*dfa*flags.swc[2] +p[27]*plg[1][2]

    # Symmetrical annual
    t[3] = p[19]*cd32

    # Symmetrical semiannual
    t[4] =  (p[16]+p[17]*plg[1][3])*cd18

    # Asymmetrical annual
    t[5] =  f1*(p[10]*plg[1][2]+p[11]*plg[1][4])*cd14

    # Asymmetrical semiannual
    t[6] = p[38]*plg[1][2]*cd39

    # Diurnal
    if flags.sw[7] != 0
        t71  = (p[12]*plg[2][3])*cd14*flags.swc[6]
        t72  = (p[13]*plg[2][3])*cd14*flags.swc[6]
        t[7] = (f2*((p[4]*plg[2][2] + p[5]*plg[2][4] + p[28]*plg[2][6] + t71) *
               ctloc + (p[7]*plg[2][2] + p[8]*plg[2][4] + p[29]*plg[2][6] + t72)*stloc))
    end

    # Semidiurnal
    if flags.sw[9] != 0
        t81  = (p[24]*plg[3][4]+p[36]*plg[3][6])*cd14*flags.swc[6]
		t82  = (p[34]*plg[3][4]+p[37]*plg[3][6])*cd14*flags.swc[6]
        t[8] = (f2*((p[6]*plg[3][3]+ p[42]*plg[3][5] + t81)*c2tloc +
               (p[9]*plg[3][3] + p[43]*plg[3][5] + t82)*s2tloc))
    end

    # Terdiurnal
    if flags.sw[15] != 0
        t[14] = (f2*((p[40]*plg[4][4]+(p[94]*plg[4][5]+p[47]*plg[4][7])*cd14*flags.swc[6])*s3tloc
                 + (p[41]*plg[4][4]+(p[95]*plg[4][5]+p[49]*plg[4][7])*cd14*flags.swc[6])*c3tloc))
    end

    # Magnetic activity based on daily AP
    if flags.sw[10] == -1
        if p[52] != 0
            exp1 = exp(-10800.0*sqrt(p[52]*p[52])/(1.0+p[139]*(45.0-sqrt(input.g_lat*input.g_lat))))
        end

        if exp1 > 0.99999
            exp1 = 0.99999
        end

        if p[25] < 1.0e-4
            p[25] = 1.0e-4
        end

        apt[1] = sg0(exp1, p, input.ap_array)

        if flags.sw[10] != 0
            t[8] = (apt[1]*(p[51]+p[97]*plg[1][3]+p[55]*plg[1][5] +
                    (p[126]*plg[1][2]+p[127]*plg[1][4]+p[128]*plg[1][6])*cd14*flags.swc[6] +
                    (p[129]*plg[2][2]+p[130]*plg[2][4]+p[131]*plg[2][6])*flags.swc[8] *
                    cos(hr*(tloc-p[132]))))
        end
    else
        apd = input.ap - 4.0
        p44 = p[44]
        p45 = p[45]
        if p44 < 0
            p44 = 1.0e-5
        end
        apdf = apd + (p45-1.0)*(apd + (exp(-p44 * apd) - 1.0)/p44)
        if flags.sw[10] != 0
            t[9] = (apdf*(p[33]+p[46]*plg[1][3]+p[35]*plg[1][5] +
                    (p[101]*plg[1][2]+p[102]*plg[1][4]+p[103]*plg[1][6])*cd14*flags.swc[6] +
                    (p[122]*plg[2][2]+p[123]*plg[2][4]+p[124]*plg[2][6])*flags.swc[8]*
                    cos(hr*(tloc-p[125]))))
        end
    end

    if flags.sw[11] != 0 && (input.g_lon > - 1000.0) != 0
        # Longitudinal
        if flags.sw[12] != 0
			t[11] = ((1.0 + p[81]*dfa*flags.swc[2])*
                    ((p[65]*plg[2][3]+p[66]*plg[2][5]+p[67]*plg[2][7]
                    + p[104]*plg[2][2]+p[105]*plg[2][4]+p[106]*plg[2][6]
                    + flags.swc[6]*(p[110]*plg[2][2]+p[111]*plg[2][4]+p[112]*plg[2][6])*cd14)*
                        cos(dgtr*input.g_lon)
                    + (p[91]*plg[2][3]+p[92]*plg[2][5]+p[93]*plg[2][7]
                    + p[107]*plg[2][2]+p[108]*plg[2][4]+p[109]*plg[2][6]
                    + flags.swc[6]*(p[113]*plg[2][2]+p[114]*plg[2][4]+p[115]*plg[2][6])*cd14)*
                    sin(dgtr*input.g_lon)))
        end

        # UT and mixed UT, longitude
        if flags.sw[13] != 0
            t[12] = ((1.0+p[96]*plg[1][2])*(1.0+p[82]*dfa*flags.swc[2])*
				    (1.0+p[120]*plg[1][2]*flags.swc[6]*cd14)*
				    ((p[69]*plg[1][2]+p[70]*plg[1][4]+p[71]*plg[1][6])*
                    cos(sr*(input.sec-p[72]))))

			t[12] += (flags.swc[12]*(p[77]*plg[3][4]+p[78]*plg[3][6]+p[79]*plg[3][8])*
				    cos(sr*(input.sec-p[80])+2.0*dgtr*input.g_lon)*(1.0+p[138]*dfa*flags.swc[2]))
        end

        # UT longitude magnetic activity
        if flags.sw[14] != 0
            if flags.sw[10] == -1
                if p[52] != 0.0
                    t[13] = (apt[1]*flags.swc[12]*(1.0+p[133]*plg[1][2])*
                            ((p[53]*plg[2][3]+p[99]*plg[2][5]+p[68]*plg[2][7])*
                            cos(dgtr*(input.g_lon-p[98])))
                            +apt[1]*flags.swc[12]*flags.swc[6]*
                            (p[134]*plg[2][2]+p[135]*plg[2][4]+p[136]*plg[2][6])*
                            cd14*cos(dgtr*(input.g_lon-p[137]))
                            +apt[1]*flags.swc[13]*
                            (p[56]*plg[1][2]+p[57]*plg[1][4]+p[58]*plg[1][6])*
                            cos(sr*(input.sec-p[59])))
                end
            else
                t[13] = (apdf*flags.swc[12]*(1.0+p[121]*plg[1][2])*
					((p[61]*plg[2][3]+p[62]*plg[2][5]+p[63]*plg[2][7])*
					cos(dgtr*(input.g_lon-p[64])))
					+ apdf*flags.swc[12]*flags.swc[6] *
					(p[116]*plg[2][2]+p[117]*plg[2][4]+p[118]*plg[2][6])*
					cd14*cos(dgtr*(input.g_lon-p[119]))
					+ apdf*flags.swc[13]*
					(p[84]*plg[1][2]+p[85]*plg[1][4]+p[86]*plg[1][6])*
                    cos(sr*(input.sec-p[76])))
            end
        end

    end

    # Params not used: 82, 89, 99, 139-149
    tinf = p[31]
    for i in 1:14
        tinf = tinf + abs(flags.sw[i+1])*t[i]
    end

    # println("t: $t")
    # println("tinf: $tinf dfa: $dfa ctloc: $ctloc stloc: $stloc c2tloc: $c2tloc s2tloc: $s2tloc c3tloc: $c3tloc s3tloc: $s3tloc apdf: $apdf")

    return tinf, dfa, plg, ctloc, stloc, c2tloc, s2tloc, c3tloc, s3tloc, apdf, apt
end

"""
Calculate G(L) function for lower atmosphere
"""
function glob7s(p::Array{<:Real, 1}, input::NRLMSISE_Input, flags::NRLMSISE_Flags;
                dfa::Real, plg::Array{<:Array{<:Real,1},1}, ctloc::Real, stloc::Real, c2tloc::Real, s2tloc::Real, s3tloc::Real, c3tloc::Real)
    # Working variables
    pset = 2.0
    t    = zeros(Float64, 14)
    dr   = 1.72142E-2
    dgtr = 1.74533E-2

    # Confirm parameter set
    if p[100] == 0
        p[100] = pset
    end

    if p[100] != pset
        @error "Wrong parameter set for glob7s"
        return -1
    end

    for j in 1:14
        t[j] = 0.0
    end

	cd32 = cos(dr*(input.doy - p[32]))
	cd18 = cos(2.0*dr*(input.doy - p[18]))
	cd14 = cos(dr*(input.doy - p[14]))
	cd39 = cos(2.0*dr*(input.doy - p[39]))

    # F10.7
    t[1] = p[22]*dfa

    # Time independent
    t[2] = p[2]*plg[1][3] + p[3]*plg[1][5] + p[23]*plg[1][7] + p[27]*plg[1][2] + p[15]*plg[1][4] + p[60]*plg[1][6]

    # Symmetrical annual
    t[3] = (p[19]+p[48]*plg[1][3]+p[30]*plg[1][5])*cd32

    # Symmetrical semi annual
    t[4]=(p[16]+p[17]*plg[1][3]+p[31]*plg[1][5])*cd18

    # Asymmetrical annual
    t[5] = (p[10]*plg[1][2]+p[11]*plg[1][4]+p[21]*plg[1][6])*cd14

    # Asymmetrical semiannual
    t[6] = (p[38]*plg[1][2])*cd39

    # Diurnal
    if flags.sw[8] != 0
        t71  = p[12]*plg[2][3]*cd14*flags.swc[6]
		t72  = p[13]*plg[2][3]*cd14*flags.swc[6]
		t[7] = ((p[4]*plg[2][2] + p[5]*plg[2][4] + t71) * ctloc + (p[7]*plg[2][2] + p[8]*plg[2][4] + t72) * stloc)
    end

    # Semidiurnal
    if flags.sw[9] != 0
        t81  = (p[24]*plg[3][4]+p[36]*plg[3][6])*cd14*flags.swc[6]
		t82  = (p[34]*plg[3][4]+p[37]*plg[3][6])*cd14*flags.swc[6]
        t[8] = ((p[6]*plg[3][3] + p[42]*plg[3][5] + t81) * c2tloc + (p[9]*plg[3][3] + p[43]*plg[3][5] + t82) * s2tloc)
    end

    # Terdiurnal
    if flags.sw[15] != 0
        t[14] = p[40] * plg[4][4] * s3tloc + p[41] * plg[4][4] * c3tloc
    end

    # Magnetic activity
    if !((flags.sw[11]==0) || (flags.sw[12]==0) || (input.g_lon <= -1000.0))
        t[11] = (1.0 + plg[1][2]*(p[81]*flags.swc[6]*cos(dr*(input.doy-p[82]))
		        + p[86]*flags.swc[7]*cos(2.0*dr*(input.doy-p[87])))
			    + p[84]*flags.swc[4]*cos(dr*(input.doy-p[85]))
			    + p[88]*flags.swc[5]*cos(2.0*dr*(input.doy-p[89])))*((p[65]*plg[2][3]+p[66]*plg[2][5]+p[67]*plg[2][7]
			    + p[75]*plg[2][2]+p[76]*plg[2][4]+p[77]*plg[2][6])*cos(dgtr*input.g_lon)
			    + (p[91]*plg[2][3]+p[92]*plg[2][5]+p[93]*plg[2][7]
			    + p[78]*plg[2][2]+p[79]*plg[2][4]+p[80]*plg[2][6])*sin(dgtr*input.g_lon))
    end

    tt = 0
    for i in 1:14
        tt += abs(flags.sw[i+1])*t[i]
    end

    return tt
end
!
function gtd7!(input::NRLMSISE_Input, flags::NRLMSISE_Flags, output::NRLMSISE_Output)
    # Working output varriable
    soutput = NRLMSISE_Output()
    tz = 0.0
    dm28m = 0.0

    # MESO7
    meso_tn1  = zeros(Float64, 5)
    meso_tn2  = zeros(Float64, 4)
    meso_tn3  = zeros(Float64, 5)
    meso_tgn1 = zeros(Float64, 2)
    meso_tgn2 = zeros(Float64, 2)
    meso_tgn3 = zeros(Float64, 2)

    # Working variables
    mn3  = 5
    mn2  = 4
    zn3  = [32.5, 20.0, 15.0, 10.0, 0.0]
    zn2  = [72.5, 55.0, 45.0, 32.5]
    zmix = 62.5

    # Set Input
    tselec!(flags)

    # Lattitude variation of gravity (none for flags.sw[3] = 0)
    xlat = input.g_lat
    if flags.sw[3] == 0
        xlat = 45.0
    end

    gsurf, re = glatf(xlat)

    xmm = pdm[3][5]

    # Thermosphere / mesosphere (above zn2[1])
    altt = 0.0
    if input.alt > zn2[1]
        altt = input.alt
    else
        altt = zn2[1]
    end

    tmp = input.alt
    input.alt = altt

    dm28, meso_tn1, meso_tgn1, dfa, plg, ctloc, stloc, c2tloc, s2tloc, s3tloc, c3tloc = gts7!(input, flags, soutput, gsurf=gsurf, re=re)

    # println(soutput.d)
    # println(soutput.t)

    altt = input.alt
    input.alt = tmp

    if flags.sw[1] != 0 # Metric adjustment != 0
        dm28m = dm28*1.0e6
    else
        dm28m = dm28
    end

    output.t[1] = soutput.t[1]
    output.t[2] = soutput.t[2]

    if input.alt >= zn2[1]
        for i in 1:9
            output.d[i] = soutput.d[i]
        end
        # for i in 1:8
        #     @printf "output.d[%d]: %.10f\n" i output.d[i]
        # end
        return
    end

    # Low Mesosphere/Upper stratosphere between zn3[1] and zn2[1]
    meso_tgn2[1] = meso_tgn1[2]
    meso_tn2[1]  = meso_tn1[5]
    meso_tn2[2]  = pma[1][1]*pavgm[1]/(1.0-flags.sw[21]*glob7s(pma[1], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))
    meso_tn2[3]  = pma[2][1]*pavgm[2]/(1.0-flags.sw[21]*glob7s(pma[2], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))
    meso_tn2[4]  = pma[3][1]*pavgm[3]/(1.0-flags.sw[21]*flags.sw[23]*glob7s(pma[3], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))
    meso_tgn2[2] = pavgm[9]*pma[10][1]*(1.0+flags.sw[21]*flags.sw[23]*glob7s(pma[10], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))*meso_tn2[4]*meso_tn2[4]/((pma[3][1]*pavgm[3])^2.0)
    meso_tn3[1]  = meso_tn2[4]

    # Lower stratrosphere and troposphere below (zn3[1])
    if input.alt < zn3[1]
        meso_tgn3[1] = meso_tgn2[2]
        meso_tn3[2]  = pma[4][1]*pavgm[4]/(1.0-flags.sw[23]*glob7s(pma[4], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))
        meso_tn3[3]  = pma[5][1]*pavgm[5]/(1.0-flags.sw[23]*glob7s(pma[5], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))
        meso_tn3[4]  = pma[6][1]*pavgm[6]/(1.0-flags.sw[23]*glob7s(pma[6], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))
        meso_tn3[5]  = pma[7][1]*pavgm[7]/(1.0-flags.sw[23]*glob7s(pma[7], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))
        meso_tgn3[2] = pma[8][1]*pavgm[8]*(1.0+flags.sw[23]*glob7s(pma[8], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))*meso_tn3[5]*meso_tn3[5]/((pma[7][1]*pavgm[7])^2.0)
    end

    # println("meso_tn2: $meso_tn2")
    # println("meso_tn3: $meso_tn3")
    # println("meso_tgn2: $meso_tgn2")
    # println("meso_tgn3: $meso_tgn3")

    # Linear transition to full mixing below zn2[0]
    dmc = 0
    if input.alt > zmix
        dmc = 1.0 - (zn2[1]-input.alt)/(zn2[1] - zmix)
    end
    dz28 = soutput.d[3]

    # N2 Density
    dmr = soutput.d[3] / dm28m - 1.0
    tz, output.d[3] = densm!(input.alt, dm28m, xmm, tz,
                        mn3, zn3, meso_tn3, meso_tgn3,
                        mn2, zn2, meso_tn2, meso_tgn2,
                        gsurf=gsurf, re=re)
    output.d[3] = output.d[3] * (1.0 + dmr*dmc)

    # HE density
    dmr         = soutput.d[1] / (dz28 * pdm[1][2]) - 1.0
    output.d[1] = output.d[3] * pdm[1][2] * (1.0 + dmr*dmc)

    # O density
    output.d[2] = 0
    output.d[9] = 0

    # O2 density
    dmr = soutput.d[4] / (dz28 * pdm[4][2]) - 1.0
    output.d[4] = output.d[3] * pdm[4][2] * (1.0 + dmr*dmc)

    # AR density
    dmr = soutput.d[5] / (dz28 * pdm[5][2]) - 1.0
    output.d[5] = output.d[3] * pdm[5][2] * (1.0 + dmr*dmc)

    # H density
    output.d[7] = 0

    # Atomic N density
    output.d[8] = 0

    # Total mass density
    output.d[6] = 1.66e-24 * (  4.0 * output.d[1]
                             + 16.0 * output.d[2]
                             + 28.0 * output.d[3]
                             + 32.0 * output.d[4]
                             + 40.0 * output.d[5]
                             +  1.0 * output.d[7]
                             + 14.0 * output.d[8])

    # Correct units
    if flags.sw[1] != 0
        output.d[6] = output.d[6]/1000
    end

    # for i in 1:8
    #     @printf "output.d[%d]: %.10f\n" i output.d[i]
    # end

    # Temperature at altitude
    tz, dd = densm!(input.alt, 1.0, 0, tz,
                    mn3, zn3, meso_tn3, meso_tgn3,
                    mn2, zn2, meso_tn2, meso_tgn2,
                    gsurf=gsurf, re=re)
    output.t[2] = tz
end

function gtd7d!(input::NRLMSISE_Input, flags::NRLMSISE_Flags, output::NRLMSISE_Output)
    gtd7!(input, flags, output)

    output.d[6] = 1.66e-24 * (  4.0 * output.d[1]
                             + 16.0 * output.d[2]
                             + 28.0 * output.d[3]
                             + 32.0 * output.d[4]
                             + 40.0 * output.d[5]
                             +  1.0 * output.d[7]
                             + 14.0 * output.d[8])

    if flags.sw[1] != 0
        output.d[6] = output.d[6]/1000
    end
end

"""
Thermospheric portion of NRLMSISE-00.
See GTD7 for more extensive comments
For alt > 72.5 km
"""
function gts7!(input::NRLMSISE_Input, flags::NRLMSISE_Flags, output::NRLMSISE_Output; gsurf::Real, re::Real)
    dm28  = 0.0
    tz    = 0.0
    mn1   = 5
    zn1   = [120.0, 110.0, 100.0, 90.0, 72.5]
    dgtr  = 1.74533E-2
    dr    = 1.72142E-2
    alpha = [-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0]
    altl  = [200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0]

    # MESO7
    meso_tn1  = zeros(Float64, 5)
    meso_tgn1 = zeros(Float64, 2)

    za     = pdl[2][16]
    zn1[1] = za

    for j in 1:9
        output.d[j] = 0.0
    end

    # Tinf variations not important below za or zn[1]
    if input.alt > zn1[1]
        tinf_g7, dfa, plg, ctloc, stloc, c2tloc, s2tloc, c3tloc, s3tloc, apdf, apt = globe7(pt,input,flags)
        # println("globe7 temp: $tinf_g7")
        tinf = ptm[1]*pt[1] * (1.0+flags.sw[17]*tinf_g7)
    else
        tinf = ptm[1]*pt[1]
        # println("branch 2")
    end
    output.t[1] = tinf

    # println("tinf: $tinf")

    # Gradient variations not important below zn1[6]
    if input.alt > zn1[5]
        tinf_g7, dfa, plg, ctloc, stloc, c2tloc, s2tloc, c3tloc, s3tloc, apdf, apt = globe7(ps, input, flags)
        g0 = ptm[4]*ps[1] * (1.0 + flags.sw[20]*tinf_g7)
    else
        g0 = ptm[4]*ps[1]
    end

    tinf_g7, dfa, plg, ctloc, stloc, c2tloc, s2tloc, c3tloc, s3tloc, apdf, apt = globe7(pd[4], input, flags)
    tlb = ptm[2] * (1.0 + flags.sw[18]*tinf_g7)*pd[4][1]
    s   = g0 / (tinf - tlb)

    # Lower thermosphere temp variations now significant for density above 300km
    if input.alt < 300.0
        meso_tn1[2]  = ptm[7]*ptl[1][1]/(1.0-flags.sw[19]*glob7s(ptl[1], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))
        meso_tn1[3]  = ptm[3]*ptl[2][1]/(1.0-flags.sw[19]*glob7s(ptl[2], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))
        meso_tn1[4]  = ptm[8]*ptl[3][1]/(1.0-flags.sw[19]*glob7s(ptl[3], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))
        meso_tn1[5]  = ptm[5]*ptl[4][1]/(1.0-flags.sw[19]*flags.sw[21]*glob7s(ptl[4], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))
        meso_tgn1[2] = ptm[9]*pma[9][1]*(1.0+flags.sw[19]*flags.sw[21]*glob7s(pma[9], input, flags, dfa=dfa, plg=plg, ctloc=ctloc, stloc=stloc, c2tloc=c2tloc, s2tloc=s2tloc, s3tloc=s3tloc, c3tloc=c3tloc))*meso_tn1[5]*meso_tn1[5]/((ptm[5]*ptl[4][1])^2.0)
    else
        meso_tn1[2]  = ptm[7]*ptl[1][1]
        meso_tn1[3]  = ptm[3]*ptl[2][1]
        meso_tn1[4]  = ptm[8]*ptl[3][1]
        meso_tn1[5]  = ptm[5]*ptl[4][1]
        meso_tgn1[2] = ptm[9]*pma[9][1]*meso_tn1[5]*meso_tn1[5]/((ptm[5]*ptl[4][1])^2.0)
    end

    # N2 variations factor at Zlb
    tinf_g7, dfa, plg, ctloc, stloc, c2tloc, s2tloc, c3tloc, s3tloc, apdf, apt = globe7(pd[3], input, flags)
    g28 = flags.sw[22]*tinf_g7

    # Variation of turbopause height
    zhf         = pdl[2][25]*(1.0+flags.sw[6]*pdl[1][25]*sin(dgtr*input.g_lat)*cos(dr*(input.doy-pt[14])))
	output.t[1] = tinf
	xmm         = pdm[3][5]
    z           = input.alt

    # N2 Density
    # Diffusive density at Zlb
    db28 = pdm[3][1]*exp(g28)*pd[3][1]

    # Diffusive density at Alt
    output.t[2], output.d[3] = densu!(z, db28, tinf, tlb, 28.0, alpha[3], output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
    dd = output.d[3]

    # println("output.d[3]: $(output.d[3])")
    # return 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0

    # Turbopause
    zh28  = pdm[3][3]*zhf
    zhm28 = pdm[3][4]*pdl[2][6]
    xmd   = 28.0 - xmm

    # Mixed density at Zlb
    tz, b28 = densu!(zh28, db28, tinf, tlb, xmd, (alpha[3]-1.0), tz,ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
    if flags.sw[16] != 0 && (z <= altl[3]) != 0
        # Mixed density at alt
        tz, dm28 = densu!(z, b28, tinf, tlb, xmm, alpha[3], tz, ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
        # Net density at alt
        output.d[3] = dnet(output.d[3], dm28, zhm28, xmm, 28.0)
    end

    # println("output.d[3]: $(output.d[3])")

    # HE density
    # Density variation factor at Zlb
    tinf_g7, dfa, plg, ctloc, stloc, c2tloc, s2tloc, c3tloc, s3tloc, apdf, apt = globe7(pd[1], input, flags)
    g4 = flags.sw[22]*tinf_g7

    # Diffusive denity at Zlb
    db04 = pdm[1][1]*exp(g4)*pd[1][1]

    # Diffusive density at alt
    tz, output.d[1] = densu!(z, db04, tinf, tlb, 4.0, alpha[1], output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
    dd = output.d[1]
    if flags.sw[16] != 0 && (z < altl[1]) != 0
        # Turbo pause
        zh04 = pdm[1][3]

        # Mixed density at Zlb
        tz, b04 = densu!(zh04, db04, tinf, tlb, 4.0-xmm, alpha[1]-1.0, output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)

        # Mixed density at alt
        tz, dm04  = densu!(z, b04, tinf, tlb, xmm, 0.0, output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
        zhm04 = zhm28

        # Net density at alt
        output.d[1] = dnet(output.d[1], dm04, zhm04, xmm, 4.0)

        # Correction to specified mixing ratio at ground
        rl = log(b28*pdm[1][2]/b04)
        zc04 = pdm[1][5]*pdl[2][1]
        hc04 = pdm[1][6]*pdl[2][2]

        # Net density corrected at alt
        output.d[1] = output.d[1]*ccor(z, rl, hc04, zc04)
    end

    # O density

    #  Density variation factor at Zlb
    tinf_g7, dfa, plg, ctloc, stloc, c2tloc, s2tloc, c3tloc, s3tloc, apdf, apt = globe7(pd[2], input, flags)
	g16 = flags.sw[22]*tinf_g7

    #  Diffusive density at Zlb
	db16 = pdm[2][1]*exp(g16)*pd[2][1]

    # Diffusive density at Alt
	tz, output.d[2] = densu!(z, db16, tinf, tlb, 16.0, alpha[2],output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
	dd=output.d[2]
	if flags.sw[16] != 0 && (z <= altl[2])
		# Turbopause
        zh16=pdm[2][3]

		# Mixed density at Zlb
        tz, b16 = densu!(zh16, db16, tinf, tlb, 16.0-xmm, (alpha[2]-1.0), output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)

		# Mixed density at Alt
		tz, dm16  = densu!(z, b16, tinf, tlb, xmm, 0.0, output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
        zhm16 = zhm28

		# Net density at Alt
		output.d[2]=dnet(output.d[2], dm16, zhm16, xmm, 16.0)
        rl    = pdm[2][2]*pdl[2][17]*(1.0+flags.sw[2]*pdl[1][24]*(input.f107A - 150.0))
        hc16  = pdm[2][6]*pdl[2][4]
        zc16  = pdm[2][5]*pdl[2][3]
        hc216 = pdm[2][6]*pdl[2][5]
		output.d[2] = output.d[2]*ccor2(z, rl, hc16, zc16, hc216)

		# Chemistry correction
		hcc16 = pdm[2][8]*pdl[2][14]
		zcc16 = pdm[2][7]*pdl[2][13]
		rc16  = pdm[2][4]*pdl[2][15]

		# Net density corrected at Alt
		output.d[2] = output.d[2]*ccor(z, rc16, hcc16, zcc16)
    end

    # O2 Density

    # Density variation factor at Zlb
    tinf_g7, dfa, plg, ctloc, stloc, c2tloc, s2tloc, c3tloc, s3tloc, apdf, apt = globe7(pd[5], input, flags)
    # @printf "globe7: %.18f\n" tinf_g7
	g32 = flags.sw[22]*tinf_g7

    # Diffusive density at Zlb
    db32 = pdm[4][1]*exp(g32)*pd[5][1]

    # Diffusive density at Alt
	tz, output.d[4] = densu!(z, db32, tinf, tlb, 32.0, alpha[4], output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
    dd=output.d[4]

    if flags.sw[16] != 0
		if z <= altl[4]
			# Turbopause
            zh32 = pdm[4][3]

			# Mixed density at Zlb
            tz, b32 = densu!(zh32, db32, tinf, tlb, 32.0-xmm, alpha[4]-1.0, output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)

			# Mixed density at Alt
			tz, dm32  = densu!(z, b32, tinf, tlb, xmm, 0.0, output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
            zhm32 = zhm28

			# Net density at Alt
            output.d[4] = dnet(output.d[4], dm32, zhm32, xmm, 32.0)

			# Correction to specified mixing ratio at ground
			rl = log(b28*pdm[4][2]/b32)
			hc32 = pdm[4][6]*pdl[2][8]
			zc32 = pdm[4][5]*pdl[2][7]
			output.d[4] = output.d[4]*ccor(z, rl, hc32, zc32)
        end

        # Correction for general departure from diffusive equilibrium above Zlb
		hcc32  = pdm[4][8]*pdl[2][23]
		hcc232 = pdm[4][8]*pdl[1][23]
		zcc32  = pdm[4][7]*pdl[2][22]
		rc32   = pdm[4][4]*pdl[2][24]*(1.0+flags.sw[2]*pdl[1][24]*(input.f107A - 150.0))

        # Net density corrected at Alt
		output.d[4]=output.d[4]*ccor2(z, rc32, hcc32, zcc32, hcc232)
    end

    # AR density
    # Density variation factor at Zlb
    tinf_g7, dfa, plg, ctloc, stloc, c2tloc, s2tloc, c3tloc, s3tloc, apdf, apt = globe7(pd[6], input, flags)
	g40 = flags.sw[22]*tinf_g7

    # Diffusive density at Zlb
	db40 = pdm[5][1]*exp(g40)*pd[6][1]

    # Diffusive density at Alt
	tz, output.d[5] = densu!(z, db40, tinf, tlb, 40.0, alpha[5], output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
    dd = output.d[5]

    # println("AR calc begin")
    # println("output.d[5]: $(output.d[5])")

    if flags.sw[16] != 0 && (z <= altl[5]) != 0
		# Turbopause
		zh40 = pdm[5][3]

        # Mixed density at Zlb
		tz, b40 = densu!(zh40, db40, tinf, tlb, 40.0-xmm, alpha[5]-1.0, output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)

        # Mixed density at Alt
		tz, dm40  = densu!(z, b40, tinf, tlb, xmm, 0.0, output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
		zhm40 = zhm28

        # Net density at Alt
        output.d[5] = dnet(output.d[5], dm40, zhm40, xmm,40.0)
        # println("output.d[5]: $(output.d[5])")

        # Correction to specified mixing ratio at ground
		rl   = log(b28*pdm[5][2]/b40)
		hc40 = pdm[5][6]*pdl[2][10]
		zc40 = pdm[5][5]*pdl[2][9]

        # Net density corrected at Alt
        output.d[5] = output.d[5]*ccor(z, rl, hc40, zc40)
        # println("output.d[5]: $(output.d[5])")
    end

    # println("AR calc end")

    # Hydrogen Density
    # Density variation factor at Zlb
    tinf_g7, dfa, plg, ctloc, stloc, c2tloc, s2tloc, c3tloc, s3tloc, apdf, apt = globe7(pd[7], input, flags)
    g1 = flags.sw[22]*tinf_g7

    # Diffusive density at Zlb
    db01 = pdm[6][1]*exp(g1)*pd[7][1]

    # Diffusive density at Alt
	tz, output.d[7] = densu!(z, db01, tinf, tlb, 1.0, alpha[7], output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
    dd = output.d[7]

	if flags.sw[16] != 0 && ( z <= altl[7])
		# Turbopause
		zh01 = pdm[6][3]

        # Mixed density at Zlb
		tz, b01 = densu!(zh01, db01, tinf, tlb, 1.0-xmm, alpha[7]-1.0, output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)

        # Mixed density at Alt
		tz, dm01  = densu!(z, b01, tinf, tlb, xmm, 0.0, output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
		zhm01 = zhm28

        # Net density at Alt
		output.d[7] = dnet(output.d[7], dm01, zhm01, xmm, 1.0)

        # Correction to specified mixing ratio at ground
		rl   = log(b28*pdm[6][2]*sqrt(pdl[2][18]*pdl[2][18])/b01)
		hc01 = pdm[6][6]*pdl[2][12]
		zc01 = pdm[6][5]*pdl[2][11]
		output.d[7] = output.d[7]*ccor(z, rl, hc01, zc01)

        # Chemistry correction
		hcc01 = pdm[6][8]*pdl[2][20]
		zcc01 = pdm[6][7]*pdl[2][19]
		rc01  = pdm[6][4]*pdl[2][21]

        # Net density corrected at Alt
		output.d[7] = output.d[7]*ccor(z, rc01, hcc01, zcc01)
    end

    # Atomic Nitrogen density
    # Density variation factor at Zlb
    tinf_g7, dfa, plg, ctloc, stloc, c2tloc, s2tloc, c3tloc, s3tloc, apdf, apt = globe7(pd[8], input, flags)
	g14 = flags.sw[22]*tinf_g7

    # Diffusive density at Zlb
	db14 = pdm[7][1]*exp(g14)*pd[8][1]

    # Diffusive density at Alt
	tz, output.d[8] = densu!(z, db14, tinf, tlb, 14.0, alpha[8], output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
	dd = output.d[8]

    if flags.sw[16] != 0 && (z <= altl[8]) != 0
		# Turbopause
		zh14 = pdm[7][3]

        # Mixed density at Zlb
		tz, b14 = densu!(zh14, db14, tinf, tlb, 14.0-xmm, alpha[8]-1.0, output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)

        # Mixed density at Alt
		tz, dm14  = densu!(z, b14, tinf, tlb, xmm, 0.0, output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
		zhm14 = zhm28

        # Net density at Alt
        output.d[8] = dnet(output.d[8], dm14, zhm14, xmm, 14.0)

        # Correction to specified mixing ratio at ground
        rl   = log(b28*pdm[7][2]*sqrt(pdl[1][3]*pdl[1][3])/b14)
		hc14 = pdm[7][6]*pdl[1][2]
        zc14 = pdm[7][5]*pdl[1][1]
        output.d[8] = output.d[8]*ccor(z, rl, hc14, zc14)

        # Chemistry correction
		hcc14 = pdm[7][8]*pdl[1][5]
		zcc14 = pdm[7][7]*pdl[1][4]
		rc14  = pdm[7][4]*pdl[1][6]

        # Net density corrected at Alt
        output.d[8] = output.d[8]*ccor(z, rc14, hcc14, zcc14)

    end

    # Anomalous Oxygen density
    tinf_g7, dfa, plg, ctloc, stloc, c2tloc, s2tloc, c3tloc, s3tloc, apdf, apt = globe7(pd[9], input, flags)
    g16h  = flags.sw[22]*tinf_g7
    db16h = pdm[8][1]*exp(g16h)*pd[9][1]
    tho   = pdm[8][10]*pdl[1][7]
    output.t[2], dd  = densu!(z, db16h, tho, tho, 16.0, alpha[9], output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)
    zsht  = pdm[8][6]
    zmho  = pdm[8][5]
    zsho  = scaleh(zmho, 16.0, tho, gsurf=gsurf, re=re)
    output.d[9] = dd*exp(-zsht/zsho*(exp(-(z-zmho)/zsht)-1.0))

    # Total mass density
    output.d[6] = 1.66e-24 * (  4.0 * output.d[1]
                             + 16.0 * output.d[2]
                             + 28.0 * output.d[3]
                             + 32.0 * output.d[4]
                             + 40.0 * output.d[5]
                             +  1.0 * output.d[7]
                             + 14.0 * output.d[8])

    # Temperature
    z = abs(input.alt)
    output.t[2], ddum = densu!(z, 1.0, tinf, tlb, 0.0, 0.0, output.t[2], ptm[6], s, mn1, zn1, meso_tn1, meso_tgn1, gsurf=gsurf, re=re)

    if flags.sw[1] != 0
        for i in 1:9
            output.d[i] = output.d[i]*1.0e6
            # @printf "output.d[%d]: %.10f\n" i output.d[i]
        end
        output.d[6] = output.d[6]/1000
        # @printf "output.d[%d]: %.10f\n" 6 output.d[6]
    end

    return dm28, meso_tn1, meso_tgn1, dfa, plg, ctloc, stloc, c2tloc, s2tloc, s3tloc, c3tloc
end

export density_nrlmsise00
"""
Computes the local atmospheric density using the NRLMSISE00 atmosphere model.

Arguments:
- `epc::Epoch`: Epoch of computation. Used to lookup space weather data
- `x::Array{<:Real, 1}`: Satellite state in geodetic coordinates [lon, lat, alt]
- `use_degrees:Bool`: If `true` interprets geodetic inputs as being in degrees

Returns:
- `rho:Float64`: Local atmospheric density [kg/m^3]

Notes:
1. Uses the gtd7d subroutine of the NRLMSISE00 to compute local atmospheric density. This subroutine includes the contribution of anomalous Oxygen in local density which is important for satellites above 500 km altitude.

References:
1. _Picone, JM, et al._ NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues _Journal of Geophysical Research: Space Physics_
"""
function density_nrlmsise00(epc::Epoch, x::Array{<:Real, 1}; use_degrees::Bool=false)
    # Create NRLMSISE00 model input
    input  = NRLMSISE_Input()
    flags  = NRLMSISE_Flags()
    output = NRLMSISE_Output()

    # Convert geodetic position
    lon = x[1]
    lat = x[2]
    alt = x[3]

    if !use_degrees
        lon *= 180.0/pi
        lat *= 180.0/pi
    end

    # Look-up space weather data
    mjd_ut1 = mjd(epc, tsys="UT1")
    doy     = floor(Int, day_of_year(epc, tsys="UT1"))
    seconds = (mjd_ut1 - floor(mjd_ut1))*86400.0

    ap_array    = zeros(Float64, 7)
    ap_array[1] = ApDailyIndex(mjd_ut1)
    ap_array[2] = ApIndex(epc)
    ap_array[3] = ApIndex(epc - 3*3600)
    ap_array[4] = ApIndex(epc - 6*3600)
    ap_array[5] = ApIndex(epc - 9*3600)

    # Ap average for 8 3-hour segments before current time
    n = 0
    for i in 12:3:33
        n += 1
        ap_array[6] += ApIndex(epc - i*3600)
    end
    ap_array[6] = ap_array[6]/n

    # Ap average for 8 3-hour segments before current time
    n = 0
    for i in 36:3:57
        n += 1
        ap_array[7] += ApIndex(epc - i*3600)
    end
    ap_array[7] = ap_array[7]/n

    # Set model flags to default
    flags.switches = ones(Int, 24)

    # Set input values
    input.year     = 0         # Unused in model
    input.doy      = doy       # Day of Year
    input.sec      = seconds   # UT
    input.alt      = alt / 1000.0      # [km]
    input.g_lat    = lat       # [deg]
    input.g_lon    = lon       # [deg]
    input.lst      = (input.sec/3600.0 + input.g_lon/15.0) # Set LST to be self consistent
    input.f107A    = f107ObservedAvg(mjd_ut1)  # 81 day average of F10.7 flux, centered on day
    input.f107     = f107Observed(mjd_ut1)     # Daily F10.7 flux for previous day
    input.ap       = ap_array[1] # Daily magnetic index
    input.ap_array = ap_array

    # Run NRLMSISE00 Atmospheric model - including anomalous oxygen
    gtd7d!(input, flags, output)

    # Return local density
    rho = output.d[6]
    return rho
end


text:
# Autogenerated wrapper script for CALCEPH_jll for i686-w64-mingw32
export libcalceph

JLLWrappers.@generate_wrapper_header("CALCEPH")
JLLWrappers.@declare_library_product(libcalceph, "libcalceph-1.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcalceph,
        "bin\\libcalceph-1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
const data_repo_path = "https://github.com/max-andr/MIPVerify_data/raw/master"

function prep_data_file(relative_dir::String, filename::String)::String
    absolute_dir = joinpath(dependencies_path, relative_dir)
    if !ispath(absolute_dir)
        mkpath(absolute_dir)
    end

    relative_file_path = joinpath(relative_dir, filename)
    absolute_file_path = joinpath(dependencies_path, relative_file_path)
    if !isfile(absolute_file_path)
        url = joinpath(data_repo_path, relative_file_path)
        download(url, absolute_file_path)
    end

    return absolute_file_path
end


text:
motions = Dict{Char, Any}(
    'i' => (buf) -> Motion(buf),
    'I' => (buf) -> line_begin(buf),
    'a' => (buf) -> begin
        return if !eof(buf)
            return Motion(position(buf), position(buf) + 1)
        else
            return Motion(buf)
        end
    end,
    'A' => (buf) -> begin
        motion = line_end(buf)
        return if !eof(buf)
            return Motion(motion.start, motion.stop + 1)
        else
            return motion
        end
    end,
    'h' => (buf) -> Motion(position(buf), position(buf) - 1), # , exclusive
    'l' => (buf) -> Motion(position(buf), position(buf) + 1),# exclusive
    'j' => down,
    'k' => up,
    'w' => word_next, # exclusive),
    'W' => word_big_next, # exclusive),
    'e' => word_end, # inclusive,
    'E' => word_big_end, # )
    'b' => word_back, # exclusive)
    'B' => word_big_back, # exclusive)
    '^' => line_begin, # exclusive)
    '$' => line_end, # inclusive)
)

"""
    Generate a Motion object for the given `name`
"""
function gen_motion(buf, name :: Char) :: Motion
    # motions = Motion[]
    fn_name = get_safe_name(name)
    fn = if name in keys(motions)
        motions[name]
    else
        log("$name has no mapped function")
        (buf) -> Motion(buf)
    end
    # call the command's function to generate the motion object
    motion = fn(buf)
    return motion
end
"""
Generate motion for the given `name` which is a TextObject
"""
function gen_motion(buf, name :: String) :: Motion

end

function double_quote(mode::NormalMode, s::LE.MIState) :: Action
    @log vim.mode = SelectRegister()
end
special_keys = Dict(
    '`' => "backtic",
    '~' => "tilde",
    '!' => "bang",
    '@' => "at",
    '#' => "hash",
    '$' => "dollar",
    '%' => "percent",
    '^' => "caret",
    '&' => "ampersand",
    '*' => "asterisk",
    '(' => "open_paren",
    ')' => "close_paren",
    '-' => "dash",
    '_' => "underscore",
    '=' => "equals",
    '+' => "plus",
    '\\' => "backslash",
    '|' => "bar",
    '[' => "open_bracket",
    ']' => "close_bracket",
    '{' => "open_curly_brace",
    '}' => "close_curly_brace",
    "'"[1] => "single_quote",
    '"' => "double_quote",
    ';' => "semicolon",
    ':' => "colon",
    ',' => "comma",
    '<' => "open_angle_bracket",
    '>' => "close_angle_bracket",
    '.' => "dot",
    '/' => "slash",
    '?' => "question_mark"
)

all_keys = Char[collect(keys(special_keys));
                collect('a':'z');
                collect('A':'Z');
                collect('0':'9')]

"""
    Get the function-safe name for the character c
"""
function get_safe_name(c :: Char) :: Symbol
    get(special_keys, c, string(c)) |> Symbol
end

"""
    Get the function-safe name for the string s, which must be
a 1 character string
"""
function get_safe_name(s :: AbstractString) :: Symbol
    if length(s) != 1
        error("length of given name is $(length(s)). Length must be 1")
    end
    return get_safe_name(s[1])
end


text:
function solve()
    n, m = [parse(Int, x) for x in split(readline())]
    h = [parse(Int, x) for x in split(readline())]

    good = ones(Bool, n)
    for _ in 1:m
        a, b = [parse(Int, x) for x in split(readline())]
        if h[a] < h[b]
            good[a] = false
        elseif h[b] < h[a]
            good[b] = false
        else
            good[a] = good[b] = false
        end
    end
    sum(good)
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# Test the mut_evolve() function from src/Evolve.jl
# Worked 5/29/21

#include("../src/CGP.jl")
using Main.CGP
MyInt = Main.CGP.MyInt

# This is a test function but is included here for convenience
function test_mut_evolve(p::Parameters, funcs::Vector{Func}, max_steps::Int64, goallist::GoalList; use_robustness::Bool=false,
      avgfitness::Bool=false )
  #println("ngoals: ",ngoals)
  c = random_chromosome(p,funcs)
  println("gl: ",gl)
  fit_limit = Float64(c.params.numoutputs)
  (c,step,worse,same,better,output,matched_goals,matched_goals_list) =
      mut_evolve( c, gl, funcs, max_steps,  use_robustness=use_robustness, avgfitness=avgfitness, fit_limit=fit_limit, print_improvements=false )
  if step < max_steps
    @assert output == output_values(c)
  end
  (c,step,worse,same,better,output,gl,matched_goals,matched_goals_list)
end

avgfitness = true
numinputs = 4
numoutputs = 1
numinteriors = 10
use_robust=false
funcs = default_funcs(numinputs)
p = Parameters( numinputs=numinputs, numoutputs=numoutputs, numinteriors=numinteriors, numlevelsback=numinteriors )
context = construct_context(numinputs)
c = random_chromosome( p, funcs )
print_circuit(c)
maxsteps = 200000
ngoals = 1   # number goals in goallist
gl = randgoallist(ngoals,c.params.numinputs,c.params.numoutputs)
#(steps,worse,same,c,output,goallist,matched_goals,matched_goals,matched_goals_list) =
#      test_mut_evolve(p,funcs,maxsteps,ngoals, avgfitness=avgfitness)
(ec,steps,worse,same,better,output,gl,matched_goals,matched_goals_list ) =
      test_mut_evolve(p,funcs,maxsteps,gl,use_robustness=use_robust, avgfitness=avgfitness )
print_circuit(ec)
#(nc,steps) = neutral_evolution(c,gl[1],maxsteps)
#print_circuit(nc)
#println("result: ",(c,steps,worse,same,better,output,gl,matched_goals_list))
#(c,steps,worse,same,better,output,gl,matched_goals_list)


text:
module ZygoteRules

export @adjoint, @adjoint!

"""
    ZygoteRules.literal_getproperty(x, ::Val{f})

In Zygote, differentiation of property access is defined by defining adjoint of
`ZygoteRules.literal_getproperty` rather than of `Base.getproperty`.
"""
literal_getproperty(x, ::Val{f}) where f = getproperty(x, f)

include("adjoint.jl")

end


text:
include("MaterialType.jl")
include("DislocationType.jl")
include("IntegrationType.jl")
include("FEMType.jl")

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
# MadNLP.jl
# Created by Sungho Shin (sungho.shin@wisc.edu)

module Mc68

import ..MadNLPHSL: @kwdef, libhsl

@kwdef mutable struct Control
    f_array_in::Cint = 0
    f_array_out::Cint = 0
    min_l_workspace::Clong = 0
    lp::Cint = 0
    wp::Cint = 0
    mp::Cint = 0
    nemin::Cint = 0
    print_level::Cint = 0
    row_full_thresh::Cint = 0
    row_search::Cint = 0
end

@kwdef mutable struct Info
    flag::Cint = 0
    iostat::Cint = 0
    stat::Cint = 0
    out_range::Cint = 0
    duplicate::Cint = 0
    n_compressions::Cint = 0
    n_zero_eigs::Cint = 0
    l_workspace::Clong = 0
    zb01_info::Cint = 0
    n_dense_rows::Cint = 0
end

function get_mc68_default_control()
    control = Control(0,0,0,0,0,0,0,0,0,0)
    mc68_default_control_i(control)
    return control
end

mc68_default_control_i(control::Control) = ccall((:mc68_default_control_i,libhsl),
                                                 Nothing,
                                                 (Ref{Control},),
                                                 control)

mc68_order_i(ord::Int32,n::Int32,ptr::Array{Int32,1},row::Array{Int32,1},
             perm::Array{Int32,1},control::Control,info::Info) = ccall(
                 (:mc68_order_i,libhsl),
                 Nothing,
                 (Cint,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ref{Control},Ref{Info}),
                 ord,n,ptr,row,perm,control,info)

end #module


text:
# function from baggepinnen/ControlSystemIdentification
function wcfft(y, u; n = length(y) ÷ 10, noverlap = n ÷ 2, window = hamming)
    win, norm2 = DSP.Periodograms.compute_window(window, n)
    uw = arraysplit(u, n, noverlap, nextfastfft(n), win)
    yw = arraysplit(y, n, noverlap, nextfastfft(n), win)
    Syy = zeros(length(uw[1]) ÷ 2 + 1)
    Suu = zeros(length(uw[1]) ÷ 2 + 1)
    Syu = zeros(ComplexF64, length(uw[1]) ÷ 2 + 1)
    for i in eachindex(uw)
        xy = rfft(yw[i])
        xu = rfft(uw[i])
        Syu .+= xy .* conj.(xu)
        Syy .+= abs2.(xy)
        Suu .+= abs2.(xu)
    end
    return Syy, Suu, Syu, norm2
end

"""
    tfbrs(RR::Vector{<:Real}, SBP::Vector{<:Real}; ...)

Transfer function based BRS measure as defined by Robbe et al.

# Keyword Arguments

* n: length of hamming window for spectral estimation, defaults to length(RR) ÷ 10
* minCoh: minimal valid coherence, defaults to 0.5
* LF: The frequency range defined as low frequency, defaults to (0.04, 0.15)
"""
function tfbrs(RR::Vector{<:Real}, SBP::Vector{<:Real}; n = length(RR) ÷ 10, minCoh = 0.5, LF::Tuple{Float64, Float64} = (0.04, 0.15))
    0 < minCoh < 1 || throw(DomainError("Coherence is defined between 0 and 1"))
    # detrend
    SBPm = SBP .- mean(SBP)
    RRm = RR .- mean(RR)

    # Spectra estimate
    Syy, Suu, Syu, norma2 = wcfft(RRm, SBPm, n = n)

    # TF
    tf = Syu ./ Suu
    # tfest  = arx(iddata(RRm, SBPm), 30, 30)
    # mag, _ = bode(tfest,  range(0, π, length = length(Syy)))
    # tf = vec(mag)
    freqs = range(0, 0.5, length = length(tf))

    # Coherence κ
    κ = abs2.(Syu) ./ (Suu .* Syy)

    # valid indices
    indsLF = findall(x -> LF[1] <= x <= LF[2], freqs)
    inds = indsLF[findall(x -> κ[x] >= minCoh, indsLF)]
    tfBRSv = mean(abs.(tf[inds]))

    return tfBRS(value = tfBRSv, minCoh = minCoh, κ = κ, tfMag = abs.(tf), inds = inds)
end

"""
Struct that stores all information regarding the tfBRS estimation. The final result is stored in 'value'. It can be plotted for visual inspection.
"""
@with_kw mutable struct tfBRS
    value::Float64 = 0
    minCoh::Float64 = 0
    LF::Tuple{Float64, Float64} = (0.04, 0.15)
    κ::Vector{<:Real} = Float64[]
    tfMag::Vector{<:Real} = Float64[]
    inds::Vector{Int} = Int[]
end

@recipe function f(res::tfBRS)

    layout := (2,1)
    freqs = range(0, 0.5, length = length(res.tfMag))
    @series begin
        subplot := 1
        yguide := "ms/mmHg"
        label := ""
        color := :grey

        freqs, res.tfMag
    end
    @series begin
        subplot := 1
        seriestype := :vline
        yguide := "ms/mmHg"
        label := ""
        color := :red

        [res.LF...]
    end
    @series begin
        subplot := 1
        title := "tfBRS = $(round(res.value, digits = 2)) ms/mmHg"
        seriestype := :scatter
        marker := :d
        yguide := "ms/mmHg"
        label := ""
        color := :red

        freqs[res.inds], res.tfMag[res.inds]
    end

    @series begin
        subplot := 2
        yguide := "Coherence"
        label := ""
        color := :grey

        freqs, res.κ
    end
    @series begin
        subplot := 2
        seriestype := :hline
        yguide := "Coherence"
        label := ""
        color := :grey

        [res.minCoh]
    end
    @series begin
        subplot := 2
        seriestype := :vline
        yguide := "Coherence"
        xguide := "f [c/b]"
        label := ""
        color := :red

        [res.LF...]
    end

end


text:
module TestAdHocRF

using Test
using Transducers
using Transducers: wheninit, whencomplete, whencombine, complete, whenstart
using MicroCollections: EmptyVector

@testset "setters" begin
    rf = nothing |> wheninit(1) |> whenstart(2) |> whencomplete(3) |> whencombine(4)
    @test rf.next === nothing
    @test rf.oninit === 1
    @test rf.start === 2
    @test rf.complete === 3
    @test rf.combine === 4
    @test wheninit(10, rf).oninit === 10
    @test wheninit(10, rf).start === 2
end

@testset "collector!!" begin
    collector!! = push!! |> whencombine(append!!) |> wheninit(EmptyVector)
    @test collector!!.next === push!!
    @test collector!!.combine === append!!
    @test collector!!.oninit === EmptyVector
    @test complete(collector!!, Val(1)) === Val(1)
    @test foldl(collector!!, Filter(isodd), 1:5) == 1:2:5
    @test foldxl(collector!!, Filter(isodd), 1:5) == 1:2:5
    @test foldxt(collector!!, Filter(isodd), 1:5; basesize = 1) == 1:2:5
end

getoninit(rf) = rf.oninit

@testset "inference" begin
    rf = wheninit(EmptyVector, append!!)
    @test @inferred(getoninit(rf)) === EmptyVector
end

end  # module


text:
"""
    HMCDA(n_iters::Int, n_adapts::Int, delta::Float64, lambda::Float64)

Hamiltonian Monte Carlo sampler with Dual Averaging algorithm.

Usage:

```julia
HMCDA(1000, 200, 0.65, 0.3)
```

Arguments:

- `n_iters::Int` : Number of samples to pull.
- `n_adapts::Int` : Numbers of samples to use for adaptation.
- `delta::Float64` : Target acceptance rate. 65% is often recommended.
- `lambda::Float64` : Target leapfrop length.

Example:

```julia
# Define a simple Normal model with unknown mean and variance.
@model gdemo(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0, sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end

sample(gdemo([1.5, 2]), HMCDA(1000, 200, 0.65, 0.3))
```

For more information, please view the following paper ([arXiv link](https://arxiv.org/abs/1111.4246)):

Hoffman, Matthew D., and Andrew Gelman. "The No-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo." Journal of Machine Learning Research 15, no. 1 (2014): 1593-1623.
"""
mutable struct HMCDA{AD, T} <: AdaptiveHamiltonian{AD}
    n_iters   ::  Int       # number of samples
    n_adapts  ::  Int       # number of samples with adaption for epsilon
    delta     ::  Float64   # target accept rate
    lambda    ::  Float64   # target leapfrog length
    space     ::  Set{T}    # sampling space, emtpy means all
    gid       ::  Int       # group ID
end
HMCDA(args...) = HMCDA{ADBackend()}(args...)
function HMCDA{AD}(n_adapts::Int, delta::Float64, lambda::Float64, space...) where AD
    _space = isa(space, Symbol) ? Set([space]) : Set(space)
    return HMCDA{AD, eltype(_space)}(1, n_adapts, delta, lambda, _space, 0)
end
function HMCDA{AD}(n_iters::Int, delta::Float64, lambda::Float64) where AD
    n_adapts_default = Int(round(n_iters / 2))
    n_adapts = n_adapts_default > 1000 ? 1000 : n_adapts_default
    return HMCDA{AD, Any}(n_iters, n_adapts, delta, lambda, Set(), 0)
end
function HMCDA{AD}(n_iters::Int, n_adapts::Int, delta::Float64, lambda::Float64) where AD
    return HMCDA{AD, Any}(n_iters, n_adapts, delta, lambda, Set(), 0)
end
function HMCDA{AD}(n_iters::Int, n_adapts::Int, delta::Float64, lambda::Float64, space...) where AD
    _space = isa(space, Symbol) ? Set([space]) : Set(space)
    return HMCDA{AD, eltype(_space)}(n_iters, n_adapts, delta, lambda, _space, 0)
end
function HMCDA{AD1}(alg::HMCDA{AD2, T}, new_gid::Int) where {AD1, AD2, T}
    return HMCDA{AD1, T}(alg.n_iters, alg.n_adapts, alg.delta, alg.lambda, alg.space, new_gid)
end
function HMCDA{AD, T}(alg::HMCDA, new_gid::Int) where {AD, T}
    return HMCDA{AD, T}(alg.n_iters, alg.n_adapts, alg.delta, alg.lambda, alg.space, new_gid)
end

function hmc_step(θ, lj, lj_func, grad_func, H_func, ϵ, alg::HMCDA, momentum_sampler::Function;
                  rev_func=nothing, log_func=nothing)
    θ_new, lj_new, is_accept, τ_valid, α = _hmc_step(
                θ, lj, lj_func, grad_func, H_func, ϵ, alg.lambda, momentum_sampler; rev_func=rev_func, log_func=log_func)
    return θ_new, lj_new, is_accept, α
end


text:
using RingBuffers
if VERSION >= v"0.5.0-"
    using Base.Test
else
    using BaseTestNext
end

@testset "RingBuffer Tests" begin
    include("lockfree.jl")
    @testset "Construction" begin
        r = RingBuffer(Float64, 8, 1)
        @test eltype(r) == Float64
        @test RingBuffers.read_space(r) == 0
        @test RingBuffers.write_space(r) == 8
    end
    @testset "Basic read/write" begin
        r = RingBuffer(Int, 8, 2)
        data = reshape(1:12, 6, 2)
        write(r, data[1:4, :])
        @test read(r, 4) == data[1:4, :]
        write(r, data)
        @test read(r, 6) == data
    end
    @testset "Overwriting overflow" begin
        r = RingBuffer(Int, 8, 2; overflow=OVERWRITE)
        data = reshape(1:12, 6, 2)
        @test write(r, data[1:4, :]) == 4
        @test write(r, data) == 6
        @test read(r, 4) == [3 9; 4 10; 1 7; 2 8]
    end
    @testset "Single-write overwrite overflow" begin
        r = RingBuffer(Int, 8, 2; overflow=OVERWRITE)
        data = reshape(1:20, 10, 2)
        @test write(r, data) == 10
        @test read(r, 4) == [3 13; 4 14; 5 15; 6 16]
    end
    @testset "truncating overflow" begin
        r = RingBuffer(Int, 8, 2; overflow=TRUNCATE)
        data = reshape(1:12, 6, 2)
        @test write(r, data[1:4, :]) == 4
        @test write(r, data) == 4
        @test write(r, data) == 0
        @test read(r, 5) == [1 7; 2 8; 3 9; 4 10; 1 7]
    end
    @testset "truncating underflow" begin
        r = RingBuffer(Int, 8, 2; underflow=TRUNCATE)
        data = reshape(1:6, 3, 2)
        @test write(r, data) == 3
        @test read(r, 5) == [1 4; 2 5; 3 6]
    end

    @testset "padded underflow" begin
        r = RingBuffer(Int, 8, 2; underflow=PAD)
        data = reshape(1:6, 3, 2)
        @test write(r, data) == 3
        @test read(r, 5) == [1 4; 2 5; 3 6; 0 0; 0 0]
    end

    @testset "queuing blocking writes" begin
        r = RingBuffer(Int, 4, 2)
        data = reshape(1:12, 6, 2)
        res1 = @async write(r, data)
        res2 = @async write(r, data)
        res3 = @async write(r, data)
        yield()
        @test read(r, 5) == [1 7; 2 8; 3 9; 4 10; 5 11]
        @test read(r, 13) == [6 12; 1 7; 2 8; 3 9; 4 10;
                              5 11; 6 12; 1 7; 2 8; 3 9;
                              4 10; 5 11; 6 12]
        @test wait(res1) == 6
        @test wait(res2) == 6
        @test wait(res3) == 6
    end
    @testset "queueing blocking reads" begin
        r = RingBuffer(Int, 8, 2)
        data = reshape(1:12, 6, 2)
        res1 = @async read(r, 2)
        res2 = @async read(r, 2)
        res3 = @async read(r, 2)
        yield()
        write(r, data)
        @test wait(res1) == [1 7; 2 8]
        @test wait(res2) == [3 9; 4 10]
        @test wait(res3) == [5 11; 6 12]
    end
end


text:
# This adds mutators to the Vimes.jl library.
# It uses this as a guide:
# http://pitest.org/quickstart/mutators/
using InteractiveUtils
using SourceWalk
using Vimes

# Checks whether a library replacement works on a given expression.
function show_replacements(ex, fs)
    match_any = false
    Vimes.pathwalk(ex) do p, x
      for f in fs
        if Vimes.matches(f, x)
          println("starts as $(x)")
          println("matches to $(f(x))")
          match_any = true
        end
      end
      return x
    end
    match_any
end

"""
    check_replacement(example_file, mutation)

Check that a patch works on a file and runs when patched.
The `example_file` is the path to a file to mutate.
The mutation is a single mutation function. The example
file must contain a function called `eg_func` that has
a single argument. It will be called with `eg_func(3)`
and should return one value unmutated and a different
value when mutated.
"""
function check_replacement(example_file, mutation)
    include(example_file)
    original_result = Base.invokelatest(eg_func, 3)

    pf = Vimes.parsefile(example_file)
    fidx = Vimes.index(pf, [mutation])
    length(fidx) > 0 || return (false, original_result, nothing)
    mut_fn = "$(tempname()).jl"
    cp(example_file, mut_fn, force = true)
    path, patch = first(fidx)
    Vimes.apply!(mut_fn, path, patch[1])

    mutated_result = original_result
    try
        include(mut_fn)
        mutated_result = Base.invokelatest(eg_func, 3)
    catch e
        rethrow()
    finally
        rm(mut_fn)
    end
    (true, original_result, mutated_result)
end

# Turns > into >=, >= into >. The same for less-than.
function conditionals_boundary(x)
    if isexpr(x, :call)
        if x.args[1] == :>
            Expr(x.head, :>=, x.args[2:end]...)
        elseif x.args[1] == :>=
            Expr(x.head, :>, x.args[2:end]...)
        elseif x.args[1] == :<
            Expr(x.head, :<=, x.args[2:end]...)
        elseif x.args[1] == :<=
            Expr(x.head, :<, x.args[2:end]...)
        else
            nothing
        end
    else
        nothing
    end
end

# This takes the unary - to a positive.
function invert_negatives(x)
    isexpr(x) && x.args[1] == :- && length(x.args) == 2 || return
    x.args[2]
end

# Turns + -> -, - -> +, etc, for integer and float.
function math_mutator(x)
    isexpr(x) && length(x.args) == 3 || return
    a, b = x.args[2:3]
    if x.args[1] == :-
        :(isa($a, Real) && isa($b, Real) ? $a + $b : $x)
    elseif x.args[1] == :+
        :(isa($a, Real) && isa($b, Real) ? $a - $b : $x)
    elseif x.args[1] == :*
        :(isa($a, Real) && isa($b, Real) ? $a / $b : $x)
    elseif x.args[1] == :/
        :(isa($a, Real) && isa($b, Real) ? $a * $b : $x)
    elseif x.args[1] == :%
        :(isa($a, Real) && isa($b, Real) ? $a * $b : $x)
    elseif x.args[1] == :*
        :(isa($a, Real) && isa($b, Real) ? $a % $b : $x)
    elseif x.args[1] == :^
        :(isa($a, InteRealger) && isa($b, Real) ? $a * $b : $x)
    elseif x.args[1] == :&
        :(isa($a, Integer) && isa($b, Integer) ? $a | $b : $x)
    elseif x.args[1] == :|
        :(isa($a, Integer) && isa($b, Integer) ? $a & $b : $x)
    elseif x.args[1] == :<<
        :(isa($a, Integer) && isa($b, Integer) ? $a >> $b : $x)
    elseif x.args[1] == :>>
        :(isa($a, Integer) && isa($b, Integer) ? $a << $b : $x)
    elseif x.args[1] == :>>>
        :(isa($a, Integer) && isa($b, Integer) ? $a << $b : $x)
    else
        nothing
    end
end

function incode_math_mutator(a, b)
    isa(a, Real) && isa(b, Real) ? a - b : a + b
end

empty_return(x::AbstractString) = ""
empty_return(x::Real) = zero(x)
empty_return(x::Set{T}) where {T} = Set{T}()
empty_return(x::Array{T,1}) where {T} = Array{T,1}()
empty_return(x::Dict{K,V}) where {K,V} = Dict{K,V}()
empty_return(x::Char) = '\0'
empty_return(x) = x

function empty_returns(x)
    isexpr(x, :return) || return
    Expr(:return, Expr(:call, :empty_return, x.args[1]))
end

"""
Swap a returned bool for its negation. This won't do anything unless it's a bool.
It's a good example of the difficulty of doing mutation at the
syntax level when you don't know the types.
"""
function flip_returns(x)
    isexpr(x, :return) || return
    Expr(:return,
        Expr(:if,
            Expr(:call,
                :isa,
                x.args[1],
                :Bool
            ),
            Expr(:call,
                :!,
                x.args[1]
            ),
            x.args[1]
        )
    )
end

function return_nothing(x)
    isexpr(x, :return) && x.args[1] !== nothing || return
    Expr(:return, nothing)
end

function increment_integer(x)
    x isa Integer && !(x isa Bool) && return x + 1
    nothing
end

function decrement_integer(x)
    x isa Integer && !(x isa Bool) && return x - 1
    nothing
end

function scale_float(x)
    x isa AbstractFloat && return 0.9 * x
    nothing
end

function flip_and_or(x)
    isexpr(x, :&&) && return Expr(:||, x.args...)
    isexpr(x, :||) && return Expr(:&&, x.args...)
end


text:
module SetfieldImpl

export @set!!, prefermutation

using ..BangBang: setproperty!!, setindex!!
using Setfield: ComposedLens, DynamicIndexLens, IndexLens, Lens, PropertyLens, Setfield, set

struct Lens!!{T<:Lens} <: Lens
    lens::T
end

Setfield.get(obj, lens::Lens!!) = get(obj, lens.lens)

# Default to immutable:
Setfield.set(obj, lens::Lens!!, value) = set(obj, lens.lens, value)

Setfield.set(obj, lens::Lens!!{<:ComposedLens}, value) =
    set(obj, Lens!!(lens.lens.outer) ∘ Lens!!(lens.lens.inner), value)

Setfield.set(obj, ::Lens!!{<:PropertyLens{fieldname}}, value) where {fieldname} =
    setproperty!!(obj, fieldname, value)

indicesfor(lens::IndexLens, _) = lens.indices
indicesfor(lens::DynamicIndexLens, obj) = lens.f(obj)

if isdefined(Setfield, :ConstIndexLens)
    using Setfield: ConstIndexLens
    const SupportedIndexLens = Union{ConstIndexLens,DynamicIndexLens,IndexLens}
    indicesfor(::ConstIndexLens{I}, _) where {I} = I
else
    const SupportedIndexLens = Union{DynamicIndexLens,IndexLens}
end

Setfield.set(obj, lens::Lens!!{<:SupportedIndexLens}, value) =
    setindex!!(obj, value, indicesfor(lens.lens, obj)...)

"""
    prefermutation(lens::Lens) :: Lens

See also [`@set!!`](@ref).
"""
prefermutation
prefermutation(lens::Lens) = Lens!!(lens)

"""
    @set!! assignment

Like `Setfield.@set`, but prefer mutation if possible.

# Examples
```jldoctest
julia> using BangBang

julia> mutable struct Mutable
           a
           b
       end

julia> x = orig = Mutable((x=Mutable(1, 2), y=3), 4);

julia> @set!! x.a.x.a = 10;

julia> @assert x.a.x.a == orig.a.x.a == 10
```
"""
:(@set!!)
macro set!!(ex)
    Setfield.setmacro(prefermutation, ex, overwrite = true)
end

end  # module


text:
"""
Two dimensional square lattice.
"""
struct SquareLattice <: AbstractLattice
    L::Int
    sites::Int
    neighs::Matrix{Int} # row = up, right, down, left; col = siteidx
    neighs_cartesian::Array{Int, 3} # row (1) = up, right, down, left; cols (2,3) = cartesian siteidx
    lattice::Matrix{Int}

    # for generic checkerboard decomposition
    # OPT: implement Assaad's two-group square lattice version
    n_bonds::Int
    bonds::Matrix{Int} # src, trg, type
end

# constructors
"""
    SquareLattice(L::Int)

Create a square lattice with linear dimension `L`.
"""
function SquareLattice(L::Int)
    sites = L^2
    lattice = convert(Array, reshape(1:L^2, (L, L)))
    neighs, neighs_cartesian = build_neighbortable(SquareLattice, lattice, L)

    # for generic checkerboard decomposition
    n_bonds = 2*sites
    bonds = zeros(n_bonds, 3)
    bondid = 1
    for src in lattice
        nup = neighs[1, src]
        bonds[bondid,:] .= [src,nup,0]
        bondid += 1

        nright = neighs[2, src]
        bonds[bondid,:] .= [src,nright,0]
        bondid += 1
    end

    return SquareLattice(L, sites, neighs, neighs_cartesian, lattice, n_bonds, bonds)
end

function build_neighbortable(::Type{SquareLattice}, lattice, L)
    up = circshift(lattice,(-1,0))
    right = circshift(lattice,(0,-1))
    down = circshift(lattice,(1,0))
    left = circshift(lattice,(0,1))
    neighs = vcat(up[:]',right[:]',down[:]',left[:]')

    neighs_cartesian = Array{Int, 3}(undef, 4, L, L)
    neighs_cartesian[1,:,:] = up
    neighs_cartesian[2,:,:] = right
    neighs_cartesian[3,:,:] = down
    neighs_cartesian[4,:,:] = left
    return neighs, neighs_cartesian
end

# Implement AbstractLattice interface: mandatory
@inline Base.length(s::SquareLattice) = s.sites

# Implement AbstractLattice interface: optional
@inline neighbors_lookup_table(s::SquareLattice) = copy(s.neighs)

# HasNeighborsTable and HasBondsTable traits
has_neighbors_table(::SquareLattice) = HasNeighborsTable()
has_bonds_table(::SquareLattice) = HasBondsTable()


text:
using CompScienceMeshes
using BEAST
using LinearAlgebra

# G = meshsphere(1.0, 0.25)
G = meshsphere(1.0, 0.30)

c = 1.0
S = BEAST.HH3DSingleLayerTDBIO(c)
D = BEAST.HH3DDoubleLayerTDBIO(speed_of_light=c)
Id = BEAST.Identity()

# width, delay, scaling = 24.0, 36.0, 1.0
width, delay, scaling = 16.0, 24.0, 1.0
gaussian = creategaussian(width, delay, scaling)
fgaussian = fouriertransform(gaussian)
e = BEAST.planewave(point(0,0,1), c, gaussian)
de = BEAST.planewave(point(0,0,1), c, derive(gaussian))
h = BEAST.gradient(e)

# X = lagrangecxd0(G)
X = lagrangecxd0(G)
Y = duallagrangec0d1(G)
# Y = lagrangecxd0(G)

# X = duallagrangecxd0(G, boundary(G))
# Y = lagrangec0d1(G)

# Δt, Nt = 0.9032321, 301
# Δt, Nt = 1.2032321, 301
Δt, Nt = 0.16, 300
# T = timebasisc0d1(Δt, Nt)
P = timebasiscxd0(Δt, Nt)
H = timebasisc0d1(Δt, Nt)
δ = timebasisdelta(Δt, Nt)

# assemble the right hand side

bd  = assemble(n⋅h,     X ⊗ P)
Z1d = assemble(Id ⊗ Id, X ⊗ P, X ⊗ P, Val{:bandedstorage})
Z0d = assemble(D,       X ⊗ P, X ⊗ P, Val{:bandedstorage})
Zd = Z0d + (-0.5)*Z1d
u = marchonintime(inv(Zd[:,:,1]), Zd, bd, Nt)

bs = assemble(e, X ⊗ δ)
Zs = assemble(S, X ⊗ δ, X ⊗ P, Val{:bandedstorage})
v = marchonintime(inv(Zs[:,:,1]), Zs, -bs, Nt)

U, Δω, ω0 = fouriertransform(u, Δt, 0.0, 2)
V, Δω, ω0 = fouriertransform(v, Δt, 0.0, 2)

ω = collect(ω0 .+ (0:Nt-1)*Δω)
_, i1 = findmin(abs.(ω.-1.0)); ω1 = ω[i1]
U_td = U[:,i1] / fgaussian(ω1)
V_td = V[:,i1] / fgaussian(ω1)

fd_S = Helmholtz3D.singlelayer(wavenumber=1.0)
fd_e = Helmholtz3D.planewave(wavenumber=1.0, direction=ẑ)
fd_eh = assemble(strace(fd_e,G), X)
fd_Sh = assemble(fd_S, X, X)
U_fd = fd_Sh \ fd_eh

fcru_td, _ = facecurrents(U_td, X)
fcrv_td, _ = facecurrents(V_td, X)
fcru_fd, _ = facecurrents(U_fd, X)

using Plots
plot()
plot!(u[1,:])
plot!(v[1,:])

# import PlotlyJS
# ptch = patch(G, norm.(fcru_fd))
# PlotlyJS.plot(ptch)

Sstatic = Helmholtz3D.singlelayer(gamma=0.0)
Zstatic = assemble(Sstatic, X, X)


text:
# Distribution of the one-sided Kolmogorov-Smirnov test statistic:
#    D^+_n = \sup_x (\hat{F}_n(x) -F(x))
immutable KSOneSided <: ContinuousUnivariateDistribution
    n::Int
end

insupport(d::KSOneSided, x::Real) = 0.0 <= x <= 1.0

# formula of Birnbaum and Tingey (1951)
function ccdf(d::KSOneSided,x::Real)
    if x >= 1.0
        return 0.0
    elseif x <= 0.0
        return 1.0
    end
    n = d.n
    s = 0.0
    for j = 0:ifloor(n-n*x)
        p = x+j/n
        s += pdf(Binomial(n,p),j) / p
    end
    s*x
end
cdf(d::KSOneSided,x::Real) = 1.0 - ccdf(d,x)


text:
module Pretalx
using HTTP
using JSON2

export get_all_submissions

function get_all_raw(endpoint; pretalx_base, pretalx_token)
    endpoint = joinpath(pretalx_base, endpoint)
    results = Vector{NamedTuple}()
    while true
        # GOLDPLATE: This would be better done with asyncs
        resp = HTTP.get(endpoint, ["Authorization"=>pretalx_token])
        raw = JSON2.read(String(resp.body))

        append!(results, raw.results)
        raw.next === nothing && return results
        endpoint = raw.next
    end
end

get_all_submissions_raw(;kwargs...) = get_all_raw("submissions"; kwargs...)
get_all_speakers_raw(;kwargs...) = get_all_raw("speakers"; kwargs...)
get_all_reviews_raw(;kwargs...) = get_all_raw("reviews"; kwargs...)

function get_all_submissions(
        allowed_state = state -> state == "confirmed" || state=="accepted" ;
        include_raw=false, ## set this to add the raw data as well
        kwargs...
)
    submissions = get_all_submissions_raw(; kwargs...)

    viable_submissions = Vector{NamedTuple}()
    for sub in submissions
        if !allowed_state(sub.state)
            continue
        end

        submission_type = sub.submission_type isa String ? sub.submission_type : sub.submission_type.en
        speakers = getproperty.(sub.speakers, :name)
        shortsub = (
            state = sub.state,
            code=sub.code,
            title=sub.title,
            submission_type=submission_type,
            abstract = sub.abstract,
            speakers = speakers,
        )
        if include_raw
            shortsub =(shortsub..., raw=sub)
        end
        push!(viable_submissions, shortsub)
    end
    return viable_submissions
end

end # module