text:
module TestExamples2DShallowWater

using Test
using Trixi

include("test_trixi.jl")

EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem")

@testset "Shallow Water" begin
  @trixi_testset "elixir_shallowwater_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"),
      l2   = [0.9911807999456177, 0.7340358820001013, 0.7446688375175905, 0.5875351036989047],
      linf = [2.0118414065026426, 2.9946841579032526, 2.6555844309739545, 3.0],
      tspan = (0.0, 0.25))
  end

  @trixi_testset "elixir_shallowwater_well_balanced.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"),
      l2   = [0.9130579602987144, 1.0602847041965408e-14, 1.082225645390032e-14, 0.9130579602987147],
      linf = [2.113062037615659, 4.6613606802974e-14, 5.4225772771633196e-14, 2.1130620376156584],
      tspan = (0.0, 0.25))
  end

  @trixi_testset "elixir_shallowwater_source_terms.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
      l2   = [0.001868474306068482, 0.01731687445878443, 0.017649083171490863, 6.274146767717023e-5],
      linf = [0.016962486402209986, 0.08768628853889782, 0.09038488750767648, 0.0001819675955490041],
      tspan = (0.0, 0.025))
  end

  @trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
      l2   = [0.0018957692481057034, 0.016943229710439864, 0.01755623297390675, 6.274146767717414e-5],
      linf = [0.015156105797771602, 0.07964811135780492, 0.0839787097210376, 0.0001819675955490041],
      tspan = (0.0, 0.025), surface_flux=(flux_hll, flux_nonconservative_fjordholm_etal))
  end
end

end # module


text:
"""
    parse_file(file; import_all)

Parses a Matpower .m `file` or PTI (PSS(R)E-v33) .raw `file` into a
PowerModels data structure. All fields from PTI files will be imported if
`import_all` is true (Default: false).
"""

function parse_file(file::String; import_all=false, validate=true)
    pm_data = open(file) do io
        pm_data = parse_file(io; import_all=import_all, validate=validate, filetype=split(lowercase(file), '.')[end])
    end
    return pm_data
end

"Parses the iostream from a file"
function parse_file(io::IO; import_all=false, validate=true, filetype="json")
    if filetype == "m"
        pm_data = parse_matpower(io, validate=validate)
    elseif filetype == "raw"
        @info("The PSS(R)E parser currently supports buses, loads, shunts, generators, branches, transformers, and dc lines")
        pm_data = parse_psse(io; import_all=import_all, validate=validate)
    elseif filetype == "json"
        pm_data = parse_json(io; validate=validate)
    else
        @info("Unrecognized filetype")
    end

    # TODO:  not sure if this relevant for all three file types, or only .m, JJS 3/7/19
    move_genfuel_and_gentype!(pm_data)

    return pm_data
end

"""
Runs various data quality checks on a PowerModels data dictionary.
Applies modifications in some cases.  Reports modified component ids.
"""
function correct_network_data!(data::Dict{String,<:Any})
    mod_bus = Dict{Symbol,Set{Int}}()
    mod_gen = Dict{Symbol,Set{Int}}()
    mod_branch = Dict{Symbol,Set{Int}}()
    mod_dcline = Dict{Symbol,Set{Int}}()

    check_conductors(data)
    check_connectivity(data)
    check_status(data)
    check_reference_bus(data)

    make_per_unit!(data)

    mod_branch[:xfer_fix] = correct_transformer_parameters!(data)
    # mod_branch[:vad_bounds] = correct_voltage_angle_differences!(data)
    mod_branch[:mva_zero] = correct_thermal_limits!(data)
    mod_branch[:orientation] = correct_branch_directions!(data)
    check_branch_loops(data)

    mod_dcline[:losses] = correct_dcline_limits!(data)

    mod_bus[:type] = correct_bus_types!(data)
    check_voltage_setpoints(data)

    check_storage_parameters(data)
    check_switch_parameters(data)

    gen, dcline = correct_cost_functions!(data)
    mod_gen[:cost_pwl] = gen
    mod_dcline[:cost_pwl] = dcline

    simplify_cost_terms!(data)

    return Dict(
        "bus" => mod_bus,
        "gen" => mod_gen,
        "branch" => mod_branch,
        "dcline" => mod_dcline
    )
end

ext: jl


text:
include("help.jl")
function count_hor_borders(r::Robot)::Int
    moves!(r, Sud)
    moves!(r, West)

    count = 0
    side = Ost

    while !isborder(r, Nord)
        move!(r, side)
        is_hor_border = false

        while !isborder(r, side)
            if isborder(r, Nord) && !is_hor_border
                count += 1
                is_hor_border = true
            elseif !isborder(r, Nord)
                is_hor_border = false
            end

            move!(r, side)
        end

        side = inverse(side)
        move!(r, Nord)
    end

    return count
end

max_stars_repo_name: MCPavlo/mirea-progs

max_issues_repo_name: MCPavlo/mirea-progs

max_forks_repo_name: MCPavlo/mirea-progs


text:
struct Cache{A3,A4,A5,A6,A7}
    lat::A3
    JC::A4
    J::A5
    X::A6
    dXdxIF::A6
    dXdxIC::A6
    dXdx::A7
    dXdxI::A7
end
function Cache(OrdPoly, OrdPolyZ, nz, NumFaces)
    # 3
    lat    = zeros(OrdPoly+1,OrdPoly+1,NumFaces)
    # 4
    JC     = zeros(OrdPoly+1,OrdPoly+1,NumFaces,nz)
    # 5
    J      = zeros(OrdPoly+1,OrdPoly+1,OrdPolyZ+1,NumFaces,nz)
    # 6
    X      = zeros(OrdPoly+1,OrdPoly+1,OrdPolyZ+1,3,NumFaces,nz)
    dXdxIF = zeros(OrdPoly+1,OrdPoly+1,NumFaces,nz+1,3,3)
    dXdxIC = zeros(OrdPoly+1,OrdPoly+1,NumFaces,nz,3,3)
    # 7
    dXdx   = zeros(OrdPoly+1,OrdPoly+1,OrdPolyZ+1,NumFaces,nz,3,3)
    dXdxI  = zeros(OrdPoly+1,OrdPoly+1,OrdPolyZ+1,NumFaces,nz,3,3)
    # dXdx   = nothing
    # dXdxI  = nothing

    A3 = typeof(lat)
    A4 = typeof(JC)
    A5 = typeof(J)
    A6 = typeof(dXdxIF)
    A7 = typeof(dXdx)
    return Cache{A3,A4,A5,A6,A7}(
        lat,
        JC,
        J,
        X,
        dXdxIF,
        dXdxIC,
        dXdx,
        dXdxI,
    )
end

Base.@kwdef mutable struct PhysParametersStruct{C}
    cache::C
    Lx = nothing
    Ly = nothing
    H = nothing
    OrdPoly = nothing
    hS = nothing
    Grid = nothing
    ProfRho = nothing
    ProfRhoBGrd = nothing
    ProfpBGrd = nothing
    ProfTheta = nothing
    ProfVel = nothing
    Damping = nothing
    Coriolis = nothing
    Buoyancy = nothing
    Source = nothing
    xc = nothing
    d = nothing
    a = nothing
    T0 = nothing
    DeltaT = nothing
    uMax = nothing
    vMax = nothing
    NBr = nothing
    Equation = nothing
    TopoS = nothing
    NumV = nothing
    RhoPos = nothing
    uPos = nothing
    vPos = nothing
    wPos = nothing
    ThPos = nothing
    Thermo = nothing
    HyperVisc = nothing
    HyperDCurl = nothing
    HyperDGrad = nothing
    HyperDDiv = nothing
    Upwind = nothing
    Flat = nothing
    vtkFileName = nothing
    vtk = nothing
    cNames = nothing
    TBGrd = nothing
    T = nothing
    RhoBGrdF = nothing
    pBGrd = nothing
    EndTime = nothing
    RK = nothing
    ROS = nothing
    Grav = nothing
    cS = nothing
    Cpd = nothing
    Cvd = nothing
    Rd = nothing
    p0 = nothing
    Cpv = nothing
    Gamma = nothing
    kappa = nothing
    Omega = nothing
    RadEarth = nothing
    JF = nothing
    dXdxIC = nothing
    dXdxIC11 = nothing
    dXdxIC12 = nothing
    dXdxIC21 = nothing
    dXdxIC22 = nothing
    dXdxIC31 = nothing
    dXdxIC32 = nothing
    dXdxIC33 = nothing
    dXdxIF = nothing
    dXdxIF11 = nothing
    dXdxIF12 = nothing
    dXdxIF21 = nothing
    dXdxIF22 = nothing
    dXdxIF31 = nothing
    dXdxIF32 = nothing
    dXdxIF33 = nothing
    latN = nothing
    nPanel = nothing
    ModelType = nothing
    Deep = nothing
    HeightLimit = nothing
    T0E = nothing
    T0P = nothing
    B = nothing
    K = nothing
    LapseRate = nothing
    U0 = nothing
    PertR = nothing
    Up = nothing
    PertExpR = nothing
    PertLon = nothing
    PertLat = nothing
    PertZ = nothing
    StrideDamp = nothing
    Relax = nothing
    CoriolisType = nothing
    Th0 = nothing
    ExpDist = nothing
    lat0 = nothing
    lon0 = nothing
    day = nothing
    k_a = nothing
    k_f = nothing
    k_s = nothing
    DeltaT_y = nothing
    DeltaTh_z = nothing
    T_equator = nothing
    T_min = nothing
    sigma_b = nothing
    z_D = nothing
    level = nothing
    fig = nothing
    SliceXY = nothing
    RadPrint = nothing
    RefProfile = nothing
    Pres = nothing
    KE = nothing
    xC0 = nothing
    zC0 = nothing
    rC0 = nothing
    DeltaTh = nothing
    hC = nothing
    x0C = nothing
    aC = nothing
    alphaG = nothing
    betaG = nothing
    hH = nothing
    H0G = nothing
    uM = nothing
    lat0G = nothing
    lat1G = nothing
    eN = nothing
    CacheC1 = nothing
    CacheC2 = nothing
    CacheC3 = nothing
    CacheC4 = nothing
    CacheC5 = nothing
    CacheC6 = nothing
    CacheC7 = nothing
    CacheF1 = nothing
    CacheF2 = nothing
    CacheF3 = nothing
    CacheF4 = nothing
    CacheF5 = nothing
    CacheF6 = nothing
    CacheF7 = nothing
    Cache1 = nothing
    Cache2 = nothing
    Cache3 = nothing
    Cache4 = nothing
    FCG = nothing
    fV = nothing
    k = nothing
    Vn = nothing
    RhoCG = nothing
    v1CG = nothing
    v2CG = nothing
    wCG = nothing
    wCCG = nothing
    ThCG = nothing
    J = nothing
end

function PhysParameters(cache)
# Physical parameters
Grav=9.81;
cS=360;
Grav=9.81e0;
Cpd=1004.0e0;
Cvd=717.0e0;
Rd=Cpd-Cvd;
p0=1.0e5;
Cpv=1885.0e0;
Gamma=Cpd/Cvd;
kappa=Rd/Cpd;

# Sphere
Omega=2*pi/(24*3600);
RadEarth=6.37122e+6;

return PhysParametersStruct{typeof(cache)}(;
cache,
Grav,
cS,
Cpd,
Cvd,
Rd,
p0,
Cpv,
Gamma,
kappa,
Omega,
RadEarth)
end


text:
struct SnailfishNumber
    triplets::Vector{Tuple{Int, Int, Int}}  # value, depth, weight
end

SnailfishNumber(str::String) = parse(SnailfishNumber, str)

function Base.parse(::Type{SnailfishNumber}, str::String)
    elements = eval(Meta.parse(str))
    triplets = tripletize!(Tuple{Int, Int, Int}[], elements, 0, 1)
    SnailfishNumber(triplets)
end

function tripletize!(triplets, elements, depth, weight)
    left = (elements[1], depth + 1, weight * 3)
    right = (elements[2], depth + 1, weight * 2)
    first(left) isa Int ? push!(triplets, left) : tripletize!(triplets, left...)
    first(right) isa Int ? push!(triplets, right) : tripletize!(triplets, right...)
    triplets
end

function Base.:+(x::SnailfishNumber, y::SnailfishNumber)
    triplets = vcat(
        [(value, depth + 1, weight * 3) for (value, depth, weight) in x.triplets],
        [(value, depth + 1, weight * 2) for (value, depth, weight) in y.triplets],
    )
    reduce!(SnailfishNumber(triplets))
end

function explode!(x::SnailfishNumber)
    for i in 1:length(x.triplets) - 1
        lvalue, ldepth, lweight = x.triplets[i]
        rvalue, rdepth, rweight = x.triplets[i + 1]
        is_exploding_pair = ldepth == rdepth == 5 && lweight > rweight
        if is_exploding_pair
            i > 1 && splice!(x.triplets, i - 1, [x.triplets[i - 1] .+ (lvalue, 0, 0)])
            i < length(x.triplets) - 1 && splice!(x.triplets, i + 2, [x.triplets[i + 2] .+ (rvalue, 0, 0)])
            deleteat!(x.triplets, i + 1)
            deleteat!(x.triplets, i)
            insert!(x.triplets, i, (0, ldepth - 1, div(lweight, 3)))
            return true
        end
    end
    false
end

function split!(x::SnailfishNumber)
    for i in 1:length(x.triplets)
        value, depth, weight = x.triplets[i]
        if value > 9
            deleteat!(x.triplets, i)
            insert!(x.triplets, i, (ceil(Int, value / 2), depth + 1, weight * 2))
            insert!(x.triplets, i, (floor(Int, value / 2), depth + 1, weight * 3))
            return true
        end
    end
    false
end

function reduce!(x::SnailfishNumber)
    while true
        explode!(x) && continue
        split!(x) && continue
        break
    end
    x
end

magnitude(x::SnailfishNumber) = sum(first.(x.triplets) .* last.(x.triplets))

numbers = SnailfishNumber.(readlines("day18/numbers.txt"))

p1 = magnitude(sum(numbers))

@show p1

p2 = maximum(magnitude.(x + y for x in numbers for y in numbers if x != y))

@show p2


text:
module Minpack
# created Summer 2015
# Andrew Ning

export hybrd

const minpack = joinpath(dirname(@__FILE__), "../deps/src/libminpack")

function residual_wrapper(nw::Ptr{Cint}, xw::Ptr{Cdouble}, fw::Ptr{Cdouble}, iflag::Ptr{Cint})
    """
    Callback function from hybrd.  Unpacks the Fortran pointers, calls a
    residual function written in Julia, then stores results in another Fortran
    pointer.
    """

    # initialize
    nvec = unsafe_load(nw)  # length of arrays
    xvec = zeros(nvec)  # initialize copy

    # copy x values into Julia
    for i = 1:nvec
        xvec[i] = unsafe_load(xw, i)
    end

    # call residual function
    # note that res and res_args are global b.c. hybrd does not provide a
    # hook to pass arbitrary data through.  same approach used in scipy.optimize
    fvec = res(xvec, res_args...)

    # copy f values into C pointer
    for i = 1:nvec
        unsafe_store!(fw, fvec[i], i)
    end

    return
end

function hybrd(residual, x0, args, tol=1e-8)
    """
    Inputs
    ------
    residual : function
        function handle to a function of the form
        out = residual(x, arg1, arg2, ...)
        where out is an array that should equal all zeros
        x is the evaluation point and
        arg1, ... are other necessary parameters

    x0 : array
        an initial estimate of the solution vector

    args : tuple of other arguments that residual takes (arg1, arg2, ...)

    tol : float
        nonnegative input variable. termination occurs
        when the algorithm estimates that the relative error
        between x and the solution is at most tol.

    Outputs
    -------
    x : array of length(x0)
        the final estimate of the solution vector.

    f : array of length(x0)
        contains the functions evaluated at the output x.
    """

    # initialize
    info_in = Cint[0]  # need to set as C pointer to get data back
    x = copy(x0)  # copy so we don't modify input
    n = length(x)
    f = zeros(x)
    lwa = round(Int, (n*(3*n+13))/2)  # cast to int
    wa = zeros(lwa)

    # unfortunately using global vars is necessary because minpack does not
    # provide hooks to pass arbitrary pointers as many modern C/fortran
    # functions do. Note that scipy.optimize uses global variables for this same
    # purpose.

    global res = residual
    global res_args = args

    # define the callback function
    const res_func = cfunction(residual_wrapper, Void, (Ptr{Cint}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cint}))

    # call hybrd.  must pass by reference for Fortran
    # compilation command I used (OS X with gfortran):
    # gfortran -shared -O2 *.f -o libhybrd.dylib -fPIC
    ccall( (:hybrd1_, minpack), Void, (Ref{Void}, Ref{Cint}, Ref{Cdouble},
        Ref{Cdouble}, Ref{Cdouble}, Ref{Cint}, Ref{Cdouble}, Ref{Cint}),
        res_func, n, x, f, tol, info_in, wa, lwa)
    info = info_in[1]  # extract termination info

    if info == 0
        error("improper input parameters.")
    elseif info == 2
        warn("number of calls to fcn has reached or exceeded 200*(n+1).")
    elseif info == 3
        warn("tol is too small. no further improvement in the approximate solution x is possible.")
    elseif info == 4
        warn("iteration is not making good progress.")
    end

    return x, f
end

end # module


text:
struct IntegralEngine{B <: AbstractBlackScholesProcess} <: PricingEngine
  process::B
end

struct Integrand{P <: StrikedTypePayoff} <: Function
  payoff::P
  s0::Float64
  drift::Float64
  variance::Float64
end

function (integrand::Integrand)(x::Float64)
  temp = integrand.s0 * exp(x)
  result = integrand.payoff(temp)
  return result * exp(-(x - integrand.drift) * (x - integrand.drift) / (2.0 * integrand.variance))
end

function _calculate!(pe::IntegralEngine, opt::EuropeanOption)
  ex = opt.exercise
  payoff = opt.payoff
  variance = black_variance(pe.process.blackVolatility, ex.dates[end], payoff.strike)

  dividendDiscount = discount(pe.process.dividendYield, ex.dates[end])
  riskFreeDiscount = discount(pe.process.riskFreeRate, ex.dates[end])

  drift_ = log(dividendDiscount / riskFreeDiscount) - 0.5 * variance

  f = Integrand(payoff, state_variable(pe.process).value, drift_, variance)

  integrator = SegmentIntegral(5000)
  infinity = 10.0 * sqrt(variance)

  opt.results.value = discount(pe.process.riskFreeRate, ex.dates[end]) / sqrt(2.0 * pi * variance) *
                      integrator(f, drift_ - infinity, drift_ + infinity)

  return opt
end


text:
global plots_inited = false

function _init_plots()
    if ! plots_inited
        try
            @eval using Plots
            plots_inited = true
            return true
        catch
            error("Unable to load 'Plots'. If the module is not installed, try `Pkg.add(\"Plots\")`")
            return false
        end
    end
    true
end

@mkapprule Plot nodefault => true

@sjdoc Plot """
    Plot(args..., color => red, xlabel => "x label")

plot using `args` and keywords. Args may be lists, numbers, strings, etc.

    Plot( :( sin ) , Table(Pi * RandomReal(), [100]), color => red)

plot using a Julia function.

    Plot( Compile( Exp(-x) * x^(1/2) ) , Range(0.0,5.0,.01) )

plot a compiled Symata expression. Plotting Symata expressions directly is not yet supported.

Other examples.

```
Plot([1,2,3])
Plot([1,2,3], [3,1,2])
```

Only a part of the functionality of `Plot.jl` is implemented.
"""

# TODO: use wrap_symata as in NIntegrate
# Unfortunatley, we need a lot of logic to mimic dispatch on types that plot relies on.

# ugh. freesyms only works on Mxpr Fix this.
function wrap_expr(expr::Mxpr)
    symerror("Plotting expressions unsupported.")
    println("Checking ", expr)
    if isa(expr, Mxpr)
        syms = freesyms(expr)
#        println("Got syms ", syms)
        length(syms) != 1 && symerror("Expected only one free symbol. Got ", syms)
        expr = wrap_symata(expr,syms[1])
#        println(" Nowdsf ", expr)
    end
    expr
end

@doap function Plot(args...)
    (! _init_plots()) && return
    (nonkeywords, keyword_dict) = separate_rules(mx)
    listargs = Any[]
    for arg in nonkeywords
        if isa(arg, Mxpr{:List})
            a = margs(arg)
            na = [ x for x in a ]  # force Array of specific type. Eg. Any[cos, sin] fails.
            push!(listargs, na)
        elseif isa(arg,Mxpr)
            push!(listargs, wrap_expr(arg))
        else
            push!(listargs, arg)
        end
    end
    for (k,v) in keyword_dict
        if is_Mxpr(v)
            keyword_dict[k] = [ x for x in margs(v)]
        end
    end
    # println("list args ", listargs)
    # println("keyword args ", keyword_dict)
    plot(listargs... ; keyword_dict...)
end


text:
function bc(x,i)
    (i > 0 && i <= length(x)) || error("bc failed")
end

@inline function mygetindex(x,i)
    #@boundscheck bc(x,i)
    @boundscheck error("Always an error")
    @inbounds return x[i]
end

foo(x,i) = @inbounds return mygetindex(x,i)
bar(x,i) = return mygetindex(x,i)

#=function getindex1(n, i)
    @boundscheck error("YOU SHALL NEVER CHECK MY BOUNDS!!")
    @inbounds return n[i]
end

#foo(x) = @inbounds return getindex1(x,1)
function foo(x)
    @inbounds out = getindex1(x,1)
    out
end

julia> getindex1((1,2,3),1)
ERROR: YOU SHALL NEVER CHECK MY BOUNDS!!
 in getindex1(::Tuple{Int64,Int64,Int64}, ::Int64) at /home/aferris/.julia/v0.5/StaticArrays/perf/boundscheck.jl:2

julia> f1((1,2,3))
ERROR: StackOverflowError:
# I have also got, by modifying the above in trivial ways:
#   * ERROR: ReadOnlyMemoryError()
#   * signal (11): Segmentation fault (... quits)

julia> @code_warntype foo((1,2,3))
Variables:
  #self#::#foo
  x::Tuple{Int64,Int64,Int64}
  out::Union{}

Body:
  begin
      $(Expr(:inbounds, true))
      $(Expr(:inbounds, false))
      # meta: location /home/aferris/.julia/v0.5/StaticArrays/perf/boundscheck.jl getindex1 2
      $(Expr(:boundscheck, true))
      (Base.throw)($(Expr(:new, :(Core.ErrorException), "YOU SHALL NEVER CHECK MY BOUNDS!!")))::Union{}
      $(Expr(:boundscheck, :pop))
      $(Expr(:inbounds, true))
      $(Expr(:inbounds, :pop))
      # meta: pop location
      $(Expr(:inbounds, :pop))
      out::Union{} = nothing
      $(Expr(:inbounds, :pop))
  end::Union{}

=#
#=immutable NeverInBounds
    x::Tuple{Int,Int,Int}
end

function getindex1(n::NeverInBounds, i)
    @boundscheck error("YOU SHALL NEVER CHECK MY BOUNDS!!")
    @inbounds return n.x[1]
end

n = NeverInBounds((1,2,3))

f0(x) = getindex1(x, 1)
#function f1(x)
#    @inbounds return getindex1(x, 1)
#end
f1(x) = @inbounds return getindex1(x,1)

try
    f0(n)
    error("There should have been an error")
catch
    println("f0 throws an error")
end

#println("f1 returns $(f1(n))")
=#


text:
using Distributions: Categorical

abstract type AbstractCPEnv{TS <: AbstractTrajectoryState} <: AbstractEnv end

struct CPEnv{TS} <: AbstractCPEnv{TS}
    reward::Float32
    terminal::Bool
    state::TS
    actions_index::Vector{Int}
    legal_actions::Vector{Int}
    legal_actions_mask::Vector{Bool}
end

RLBase.action_space(env::CPEnv) = env.actions_index
RLBase.legal_action_space(env::CPEnv) = env.legal_actions
RLBase.legal_action_space_mask(env::CPEnv) = env.legal_actions_mask
RLBase.reward(env::CPEnv) = env.reward
RLBase.is_terminated(env::CPEnv) = env.terminal
RLBase.state(env::CPEnv) = env.state
RLBase.ActionStyle(::CPEnv) = FULL_ACTION_SET
RLBase.state_space(::SeaPearl.CPEnv, ::Observation{Any}, ::DefaultPlayer) = nothing

struct unmaskedCPEnv{TS} <: AbstractCPEnv{TS}
    reward::Float32
    terminal::Bool
    state::TS
    actions_index::Array{Int64, 1}
end

RLBase.action_space(env::unmaskedCPEnv) = env.actions_index
RLBase.reward(env::unmaskedCPEnv) = env.reward
RLBase.is_terminated(env::unmaskedCPEnv) = env.terminal
RLBase.state(env::unmaskedCPEnv) = env.state
RLBase.ActionStyle(::unmaskedCPEnv) = MINIMAL_ACTION_SET
RLBase.state_space(::unmaskedCPEnv, ::Observation{Any}, ::DefaultPlayer) = nothing

function (learner::DQNLearner)(env::AbstractCPEnv{TS}) where {TS <: NonTabularTrajectoryState}
    env |>
    state |>
    x ->
        send_to_device(device(learner), x) |>
        learner.approximator |>
        vec |>
        send_to_host
end

function (learner::A2CLearner)(env::AbstractCPEnv{ST}) where {ST <: NonTabularTrajectoryState}
    env |>
    state |>
    x ->
        send_to_device(device(learner), x) |>
        learner.approximator.actor |>
        vec |>
        send_to_host
end

function Array(state::T) where {T<:NonTabularTrajectoryState}
    state
end


text:
using ANISPROU
using LinearAlgebra

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

function launchall()

    ##===========================================================
    ## do the fitting of data with Beta functions
    inpdir = "inputdata/"
    protein = "IM7"
    betamix,dobs = invertITCdata(inpdir,protein)

    ## plot straight lines defined by parameters as a function of concentration
    plotparamlines(betamix)

    ##===========================================================
    ## plot single components for fixed protein concentration
    plotbetacomp1D(betamix,0.148)

    ##===========================================================
    # plot the 3D surface
    plotsurface3D(dobs,betamix,yscal=70,ymin=0.0)

    ##===========================================================
    ## plot fit to single experiments
    outdir="figs"
    plotsingleexperiments(outdir,dobs,betamix)

    ##===========================================================
    ## compute area for all Beta components at
    ##   requested protein concentration
    protcon = 0.08
    area,errarea = area_enthalpy(betamix,protcon)
    println("\nArea at protein concentration $protcon for each component: $area, error $errarea\n")

    ##===========================================================
    ## Area for a set of protein concentrations
    N = 15
    protcons = collect(LinRange(0.0,0.14,N))
    outdir = "output"
    areas,erras,linfitres,resstdev = areasvsprotcon(betamix,protcons,outdir,protein)
    ## plot areas only
    plotareavsprotcon(protein,protcons,areas,linfitres,resstdev,outdir)

    ## plot parameters lines and areas
    plotparamlines(betamix,protcons,areas)

    ##===========================================================
    ## compute volume for all Beta components within
    ##   requested bounds of protein concentration
    minprotcon = betamix.betpar.ymin
    maxprotcon = betamix.betpar.ymax
    volume,errvol = volume_enthalpy(betamix,minprotcon,maxprotcon)
    println("\nVolume within bounds for each component: $volume, error $errvol\n")

    ##===========================================================
    ## construct binding isotherm for Beta mix using stationary and inflection points
    freeSDS,Nbound = bindingisotherm_betamix(betamix,dobs)

    return betamix,dobs #,freeSDS,Nbound
end

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

function invertITCdata(inpdir::String,protein::String)

    ##===========================================
    ## read data
    ## data = readallexperiments("inputdata/",["IM7","FN3","sFN3"])
    data = readallexperiments(inpdir,[protein],discninitrows=0)

    #protein = "FN3" # "IM7" "FN3" "sFN3"
    sdscon = data[protein]["sdscon"]
    procon = data[protein]["procon"]
    enthalpy = data[protein]["enout"]
    idxdata = data[protein]["idxdata"]

    ## x axis is SDS concentration, y axis is protein concentration
    ## xy = [sdscon procon]
    dobs = ITCObsData(protein=protein,enthalpy=enthalpy,idxdata=idxdata,
                      sdsprotcon=[sdscon procon])
    # plot only the observed data
    #plotobsdata(dobs)

    ##=========================
    ## forward test
    a = 0.99*minimum(dobs.sdsprotcon[:,1]) # lower bound for beta domain
    b = 1.05*maximum(dobs.sdsprotcon[:,1]) # upper bound for beta domain

    minprotcon = minimum(dobs.sdsprotcon[:,2])
    maxprotcon = 1.05*maximum(dobs.sdsprotcon[:,2]) # 1.1*

    @show a,b,minprotcon,maxprotcon

    ## define the parameters of Beta 2D functions
    modefuny = "linear"
    konfuny  = "linear" #"constant"
    ampfuny  = "linear"
    betpar = ScaledBeta2DParams(modefuny=modefuny,konfuny=konfuny,
                                ampfuny=ampfuny,a=a,b=b,
                                ymin=minprotcon,ymax=maxprotcon)

    ##=======================================
    ## Only testing initial guess?
    ## Used to visualise the calculated data from the initial guess
    ##   compared with the actual observed data, in order to manually
    ##   find a good initial model
    onlytestinitialguess = false #true

    ##===========================================
    ## Starting model

    if protein=="IM7"
        ## starting model IM7
        ## Each row is one component (one 2D Beta function)
        ##  Just add (remove) rows to increase (decrease) the number of components
        ## Elements are: 2 for mode, 2 for the confidence parameter and
        ##   2 for the amplitude parameter

        if betpar.konfuny=="constant"
            # constant kon
            comp1 = [0.6,  1.5,  30.0, -2.5,  -5.0 ]
            comp2 = [1.7,  4.8,  50.0, -1.6,  -4.0 ]
            comp3 = [4.5,  10.0, 35.0, 0.12, 0.16 ]
            comp4 = [6.2,  15.3, 50.0, -1.6, -2.0 ]

        else
            comp1 = [0.6,  1.5,  35.0, 30.0, -2.5,  -5.0 ]
            comp2 = [1.7,  4.8,  60.0, 40.0, -1.6,  -4.0 ]
            comp3 = [4.5,  10.0, 40.0, 30.0, 0.12, 0.16 ]
            comp4 = [6.2,  15.3, 60.0, 40.0, -1.6, -2.0 ]
        end

     end

    ## mstart is a 2D array where each column represents one component
    mstart = [comp1 comp2 comp3 comp4]
    @show size(mstart)

    ## plotparamlines(BetaMix2D(betpar,mstart,dobs.protein))

    ############################################
    if onlytestinitialguess==true
        # Visualize initial guess
        plotinitialguess(betpar,dobs,mstart)

        ## stop here when testing initial models
        return nothing,nothing
    end

    ###############################################
    ## Set constraints
    ## Elements are: 2 for mode, 2 for the confidence parameter and
    ##   2 for the amplitude parameter
    # ## lower constraints
    # lowc = [betpar.a,  betpar.a, 2.1, 2.1, -20.0, -20.0]
    # ## upper constraints
    # upc  = [betpar.b, betpar.b, 500.0, 500.0, 10.0, 10.0]
    # lowconstr,upconstr = setconstraints(betpar,mstart,lowc,upc)

    if betpar.konfuny=="constant"
        ## constant kon
        # lower constraints [confidence must be >2.0]
        lcs1 = [betpar.a, betpar.a, 2.1, -20.0, -20.0]
        lcs2 = [betpar.a, betpar.a, 2.1, -20.0, -20.0]
        lcs3 = [betpar.a, betpar.a, 2.1,   0.0,   0.0]
        lcs4 = [betpar.a, betpar.a, 2.1, -20.0, -20.0]
        lowconstr = [lcs1 lcs2 lcs3 lcs4]
        # upper constraints
        ucs1 = [betpar.b, betpar.b, 500.0,  0.0,  0.0]
        ucs2 = [betpar.b, betpar.b, 500.0,  0.0,  0.0]
        ucs3 = [betpar.b, betpar.b, 500.0, 10.0, 10.0]
        ucs4 = [betpar.b, betpar.b, 500.0,  0.0,  0.0]
        upconstr = [ucs1 ucs2 ucs3 ucs4]

    else
        # lower constraints [confidence must be >2.0]
        lcs1 = [betpar.a, betpar.a, 2.1, 2.1, -20.0, -20.0]
        lcs2 = [betpar.a, betpar.a, 2.1, 2.1, -20.0, -20.0]
        lcs3 = [betpar.a, betpar.a, 2.1, 2.1,   0.0,   0.0]
        lcs4 = [betpar.a, betpar.a, 2.1, 2.1, -20.0, -20.0]
        lowconstr = [lcs1 lcs2 lcs3 lcs4]
        # upper constraints
        ucs1 = [betpar.b, betpar.b, 500.0, 500.0,  0.0,  0.0]
        ucs2 = [betpar.b, betpar.b, 500.0, 500.0,  0.0,  0.0]
        ucs3 = [betpar.b, betpar.b, 500.0, 500.0, 10.0, 10.0]
        ucs4 = [betpar.b, betpar.b, 500.0, 500.0,  0.0,  0.0]
        upconstr = [ucs1 ucs2 ucs3 ucs4]

    end

    ## check mstart...
    # @show mstart.<=lowconstr
    # @show mstart.>=upconstr
    @assert all(mstart.>=lowconstr)
    @assert all(mstart.<=upconstr)

    ###############################################
    ## run the Newton optimization
    nobs = length(dobs.enthalpy)
    stdobs = 0.2 .* ones(nobs) # standard deviation of the error on measured data
    stdobs[1:2] .= 0.8 # specify higher uncertainty for first two points
    stdobs[end-5:end] .= 0.5
    #@show stdobs
    invCd = inv(diagm(stdobs.^2)) # in this case a diagonal precision matrix

    ## IPNewton from Optim.jl, box constraints
    outdir = "output"
    betamix = solveinvprob(betpar,dobs,invCd,mstart,lowconstr,upconstr,
                           outdir,applynonlinconstr=true,constrarea=6.0) ## 2.0

    ##================================
    ## plot results
    outdir = "figs"
    plotresults(betamix,dobs,mstart,outdir)

    return betamix,dobs
end

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

function bindingisotherm_betamix(betamix::BetaMix2D,dobs::ITCObsData)

    ##===========================================
    ## Find stationary and inflection points
    ny = 4
    protcon = collect(LinRange(betamix.betpar.ymin,betamix.betpar.ymax,ny))
    statpts,inflpts = findcurvefeatures(betamix,protcon)

    # plot found points/features
    outdir = "figs"
    plotfoundfeatures(betamix,protcon,statpts,inflpts,outdir)

    ##===========================================
    ## Selection of local minima and maxima
    selectstatpts = Vector{Array{<:Real,2}}(undef,0)

    push!(selectstatpts, [ statpts[1][2] protcon[1];
	                   statpts[2][2] protcon[2];
                           statpts[3][2] protcon[3];
                           statpts[4][2] protcon[4] ] )

    push!(selectstatpts, [ statpts[2][3] protcon[2];
                           statpts[3][3] protcon[3];
                           statpts[4][3] protcon[4] ] )

    push!(selectstatpts, [ statpts[1][3] protcon[1];
                           statpts[2][5] protcon[2];
                           statpts[3][5] protcon[3];
                           statpts[4][5] protcon[4] ] )

    push!(selectstatpts, [ statpts[1][4] protcon[1];
                           statpts[2][6] protcon[2];
                           statpts[3][6] protcon[3];
                           statpts[4][6] protcon[4] ] )

    # ===========================================
    # Selection of inflection points
    selectinflpts = Vector{Array{<:Real,2}}(undef,0)

    push!(selectinflpts, [ inflpts[1][1] protcon[1];
                           inflpts[2][2] protcon[2];
                           inflpts[3][2] protcon[3];
                           inflpts[4][2] protcon[4] ] )

    push!(selectinflpts, [ inflpts[1][2] protcon[1];
                           inflpts[2][3] protcon[2];
                           inflpts[3][3] protcon[3];
                           inflpts[4][3] protcon[4] ] )

    push!(selectinflpts, [ inflpts[1][4] protcon[1];
                           inflpts[2][5] protcon[2];
                           inflpts[3][5] protcon[3];
                           inflpts[4][5] protcon[4] ] )

    ##===========================================
    ## Linear regression to get the parameters of straight lines

    ## find the straight line by least squares
    ## lstlines = [selectstatpts..., selectinflpts...]
    outdir = "output"
    freeSDS,Nbound,resstdev = calcfreeSDSNbound(protcon,statpts,inflpts,
                                                selectstatpts,selectinflpts,
                                                outdir,dobs.protein)

    ##===========================================
    ## Plot stuff
    outdir = "figs"
    plotbindisotherm(betamix,protcon,dobs,statpts,inflpts,freeSDS,Nbound,outdir,resstdev=resstdev)

    return freeSDS,Nbound
end

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


text:
abstract Node

type ScalarNode <: Node
    tag::AbstractString
    value::AbstractString
    start_mark::(@compat Union{Mark, Void})
    end_mark::(@compat Union{Mark, Void})
    style::(@compat Union{Char, Void})
end

type SequenceNode <: Node
    tag::AbstractString
    value::Vector
    start_mark::(@compat Union{Mark, Void})
    end_mark::(@compat Union{Mark, Void})
    flow_style::Bool
end

type MappingNode <: Node
    tag::AbstractString
    value::Vector
    start_mark::(@compat Union{Mark, Void})
    end_mark::(@compat Union{Mark, Void})
    flow_style::Bool
end


text:
deploydocs(
    deps = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
)


text:
"""
Struct to hold parameter and values for the energy model close to the one in
Monteith and Unsworth (2013)

# Arguments

- `aₛₕ = 2`: number of faces of the object that exchange sensible heat fluxes
- `aₛᵥ = 1`: number of faces of the object that exchange latent heat fluxes (hypostomatous => 1)
- `ε = 0.955`: emissivity of the object
- `maxiter = 10`: maximal number of iterations allowed to close the energy balance
- `ΔT = 0.01` (°C): maximum difference in object temperature between two iterations to
consider convergence

# Examples

```julia
energy_model = Monteith() # a leaf in an illuminated chamber
```
"""
struct Monteith{T,S} <: AbstractEnergyModel
    aₛₕ::S
    aₛᵥ::S
    ε::T
    maxiter::S
    ΔT::T
end

function Monteith(;aₛₕ = 2, aₛᵥ = 1, ε = 0.955, maxiter = 10, ΔT = 0.01)
    param_int = promote(aₛₕ, aₛᵥ, maxiter)
    param_float = promote(ε, ΔT)
    Monteith(param_int[1], param_int[2], param_float[1], param_int[3], param_float[2])
end

function inputs(::Monteith)
    (:Rₛ, :skyFraction, :d)
end

function outputs(::Monteith)
    (:Tₗ, :Rn, :Rₗₗ, :H, :λE, :Cₛ, :Cᵢ, :A, :Gₛ, :Gbₕ, :Dₗ, :Gbc, :iter)
end

Base.eltype(x::Monteith) = typeof(x).parameters[1]

"""
    net_radiation!(leaf::LeafModels{I,<:Monteith,A,Gs,S},meteo::Atmosphere,constants = Constants())
    net_radiation(leaf::LeafModels{I,<:Monteith,A,Gs,S},meteo::Atmosphere,constants = Constants())

Leaf energy balance according to Monteith and Unsworth (2013), and corrigendum from
Schymanski et al. (2017). The computation is close to the one from the MAESPA model (Duursma
et al., 2012, Vezy et al., 2018) here. The leaf temperature is computed iteratively to close
the energy balance using the mass flux (~ Rn - λE).

# Arguments

- `leaf::LeafModels{.,<:Monteith,.,.,.}`: A [`LeafModels`](@ref) struct holding the parameters for
the model with initialisations for:
    - `Rₛ` (W m-2): net shortwave radiation (PAR + NIR). Often computed from a light interception model
    - `skyFraction` (0-2): view factor between the object and the sky for both faces (see details).
    - `d` (m): characteristic dimension, *e.g.* leaf width (see eq. 10.9 from Monteith and Unsworth, 2013).
- `meteo`: meteorology structure, see [`Atmosphere`](@ref)
- `constants = Constants()`: physical constants. See [`Constants`](@ref) for more details

# Details

The skyFraction in the variables is equal to 2 if all the leaf is viewing is sky (e.g. in a
controlled chamber), 1 if the leaf is *e.g.* up on the canopy where the upper side of the
leaf sees the sky, and the side bellow sees soil + other leaves that are all considered at
the same temperature than the leaf, or less than 1 if it is partly shaded.

# Notes

If you want the algorithm to print a message whenever it does not reach convergence, use the
debugging mode by executing this in the REPL: `ENV["JULIA_DEBUG"] = PlantBiophysics`.

More information [here](https://docs.julialang.org/en/v1/stdlib/Logging/#Environment-variables).

# Examples

```julia
meteo = Atmosphere(T = 22.0, Wind = 0.8333, P = 101.325, Rh = 0.4490995)

# Using a constant value for Gs:
leaf = LeafModels(energy = Monteith(),
            photosynthesis = Fvcb(),
            stomatal_conductance = ConstantGs(0.0, 0.0011),
            Rₛ = 13.747, skyFraction = 1.0, d = 0.03)
PlantBiophysics.net_radiation!(leaf,meteo)
leaf.status.Rn
julia> 12.902547446281233

# Using the model from Medlyn et al. (2011) for Gs:
leaf = LeafModels(energy = Monteith(),
            photosynthesis = Fvcb(),
            stomatal_conductance = Medlyn(0.03, 12.0),
            Rₛ = 13.747, skyFraction = 1.0, PPFD = 1500.0, d = 0.03)

PlantBiophysics.net_radiation!(leaf,meteo)
leaf.status.Rn
leaf.status.Rₗₗ
leaf.status.A
leaf.status.Gₛ
leaf.status.Cₛ
leaf.status.Cᵢ
leaf.status.Gbc
```

# References

Duursma, R. A., et B. E. Medlyn. 2012. « MAESPA: a model to study interactions between water
limitation, environmental drivers and vegetation function at tree and stand levels, with an
example application to [CO2] × drought interactions ». Geoscientific Model Development 5 (4):
919‑40. https://doi.org/10.5194/gmd-5-919-2012.

Monteith, John L., et Mike H. Unsworth. 2013. « Chapter 13 - Steady-State Heat Balance: (i)
Water Surfaces, Soil, and Vegetation ». In Principles of Environmental Physics (Fourth Edition),
edited by John L. Monteith et Mike H. Unsworth, 217‑47. Boston: Academic Press.

Schymanski, Stanislaus J., et Dani Or. 2017. « Leaf-Scale Experiments Reveal an Important
Omission in the Penman–Monteith Equation ». Hydrology and Earth System Sciences 21 (2): 685‑706.
https://doi.org/10.5194/hess-21-685-2017.

Vezy, Rémi, Mathias Christina, Olivier Roupsard, Yann Nouvellon, Remko Duursma, Belinda Medlyn,
Maxime Soma, et al. 2018. « Measuring and modelling energy partitioning in canopies of varying
complexity using MAESPA model ». Agricultural and Forest Meteorology 253‑254 (printemps): 203‑17.
https://doi.org/10.1016/j.agrformet.2018.02.005.
"""
function net_radiation!(leaf::LeafModels{I,<:Monteith,A,Gs,S}, meteo::Atmosphere, constants = Constants()) where {I,A,Gs,S}

    # Initialisations
    leaf.status.Tₗ = meteo.T - 0.2
    Tₗ_new = zero(meteo.T)
    leaf.status.Cₛ = meteo.Cₐ
    leaf.status.Dₗ = meteo.VPD
    γˢ = Rbₕ = Δ = zero(meteo.T)
    leaf.status.Rn = leaf.status.Rₛ
    iter = 0
    # ?NB: We use iter = 0 and not 1 because to get the right number of iterations at the end
    # of the for loop, becaue we use iter += 1 at the end (so it increments once again)

    # Iterative resolution of the energy balance
    for i in 1:leaf.energy.maxiter

        # Update A, Gₛ, Cᵢ from leaf.status:
        assimilation!(leaf, meteo, constants)

        # Stomatal resistance to water vapor
        Rsᵥ = 1.0 / (gsc_to_gsw(mol_to_ms(leaf.status.Gₛ, meteo.T, meteo.P, constants.R, constants.K₀),
                                constants.Gsc_to_Gsw))

        # Re-computing the net radiation according to simulated leaf temperature:
        leaf.status.Rₗₗ = net_longwave_radiation(leaf.status.Tₗ,meteo.T,leaf.energy.ε,meteo.ε,
                                                leaf.status.skyFraction,constants.K₀,constants.σ)
        #= ? NB: we use the sky fraction here (0-2) instead of the view factor (0-1) because:
            - we consider both sides of the leaf at the same time (1 -> leaf sees sky on one face)
            - we consider all objects in the scene have the same temperature as the leaf
            of interest except the atmosphere. So the leaf exchange thermal energy only with
            the atmosphere. =#
        # leaf.status.Rₗₗ = (grey_body(meteo.T,1.0) - grey_body(leaf.status.Tₗ, 1.0))*leaf.status.skyFraction

        leaf.status.Rn = leaf.status.Rₛ + leaf.status.Rₗₗ

        # Leaf boundary conductance for heat (m s-1), one sided:
        leaf.status.Gbₕ = gbₕ_free(meteo.T, leaf.status.Tₗ, leaf.status.d, constants.Dₕ₀) +
                             gbₕ_forced(meteo.Wind, leaf.status.d)
        # NB, in MAESPA we use Rni so we add the radiation conductance also (not here)

        # Leaf boundary resistance for heat (s m-1):
        Rbₕ = 1 / leaf.status.Gbₕ

        # Leaf boundary resistance for water vapor (s m-1):
        Rbᵥ = 1 / gbh_to_gbw(leaf.status.Gbₕ)

        # Leaf boundary resistance for CO₂ (mol[CO₂] m-2 s-1):
        leaf.status.Gbc = ms_to_mol(leaf.status.Gbₕ, meteo.T, meteo.P, constants.R, constants.K₀) /
            constants.Gbc_to_Gbₕ

        # Update Cₛ using boundary layer conductance to CO₂ and assimilation:
        leaf.status.Cₛ = min(meteo.Cₐ, meteo.Cₐ - leaf.status.A / (leaf.status.Gbc * leaf.energy.aₛᵥ))

        # Apparent value of psychrometer constant (kPa K−1)
        γˢ = γ_star(meteo.γ, leaf.energy.aₛₕ, leaf.energy.aₛᵥ, Rbᵥ, Rsᵥ, Rbₕ)

        leaf.status.λE = latent_heat(leaf.status.Rn, meteo.VPD, γˢ, Rbₕ, meteo.Δ, meteo.ρ,
                                        leaf.energy.aₛₕ, constants.Cₚ)

        # If potential evaporation is needed, here is how to compute it:
        # γˢₑ = γ_star(meteo.γ, energy.aₛₕ, 1, Rbᵥ, 1.0e-9, Rbₕ) # Rsᵥ is inf. low
        # Ev = latent_heat(leaf.status.Rn, meteo.VPD, γˢₑ, Rbₕ, meteo.Δ, meteo.ρ, energy.aₛₕ, constants.Cₚ)

        Tₗ_new = meteo.T + (leaf.status.Rn - leaf.status.λE) /
                (meteo.ρ * constants.Cₚ * (leaf.energy.aₛₕ / Rbₕ))

        if abs(Tₗ_new - leaf.status.Tₗ) <= leaf.energy.ΔT break end

        leaf.status.Tₗ = Tₗ_new

        # Vapour pressure difference between the surface and the saturation vapour pressure:
        leaf.status.Dₗ = e_sat(leaf.status.Tₗ) - e_sat(meteo.T) *  meteo.Rh

        iter += 1
    end

    leaf.status.H = sensible_heat(leaf.status.Rn, meteo.VPD, γˢ, Rbₕ, meteo.Δ, meteo.ρ,
                                    leaf.energy.aₛₕ, constants.Cₚ)

    leaf.status.iter = iter

    @debug begin
        if iter == leaf.energy.maxiter
            "`net_radiation!` algorithm did not converge. Please check the value."
        end
    end

    # Transpiration (mol[H₂O] m-2 s-1):
    # ET = leaf.status.λE / meteo.λ * constants.Mₕ₂ₒ
    # ET / constants.Mₕ₂ₒ to get mm s-1 <=> kg m-2 s-1 <=> l m-2 s-1

    nothing
end

"""
    latent_heat(Rn, VPD, γˢ, Rbₕ, Δ, ρ, aₛₕ, Cₚ)
    latent_heat(Rn, VPD, γˢ, Rbₕ, Δ, ρ, aₛₕ)

λE -the latent heat flux (W m-2)- using the Monteith and Unsworth (2013) definition corrected by
Schymanski et al. (2017), eq.22.

- `Rn` (W m-2): net radiation. Carefull: not the isothermal net radiation
- `VPD` (kPa): air vapor pressure deficit
- `γˢ` (kPa K−1): apparent value of psychrometer constant (see [`γ_star`](@ref))
- `Rbₕ` (s m-1): resistance for heat transfer by convection, i.e. resistance to sensible heat
- `Δ` (KPa K-1): rate of change of saturation vapor pressure with temperature (see [`e_sat_slope`](@ref))
- `ρ` (kg m-3): air density of moist air.
- `aₛₕ` (1,2): number of sides that exchange energy for heat (2 for leaves)
- `Cₚ` (J K-1 kg-1): specific heat of air for constant pressure

# References

Monteith, J. and Unsworth, M., 2013. Principles of environmental physics: plants, animals, and the atmosphere. Academic Press. See eq. 13.33.

Schymanski et al. (2017), Leaf-scale experiments reveal an important omission in the Penman–Monteith equation,
Hydrology and Earth System Sciences. DOI: https://doi.org/10.5194/hess-21-685-2017. See equ. 22.

# Examples

```julia
Tₐ = 20.0 ; P = 100.0 ;
ρ = air_density(Tₐ, P) # in kg m-3
Δ = e_sat_slope(Tₐ)

latent_heat(300.0, 2.0, 0.1461683, 50.0, Δ, ρ, 2.0)
```
"""
function latent_heat(Rn, VPD, γˢ, Rbₕ, Δ, ρ, aₛₕ, Cₚ)
  (Δ * Rn + ρ * Cₚ * VPD * (aₛₕ / Rbₕ)) / (Δ + γˢ)
end

function latent_heat(Rn, VPD, γˢ, Rbₕ, Δ, ρ, aₛₕ)
    latent_heat(Rn, VPD, γˢ, Rbₕ, Δ, ρ, aₛₕ, Constants().Cₚ)
end

"""
    sensible_heat(Rn, VPD, γˢ, Rbₕ, Δ, ρ, aₛₕ, Cₚ)
    sensible_heat(Rn, VPD, γˢ, Rbₕ, Δ, ρ, aₛₕ)

H -the sensible heat flux (W m-2)- using the Monteith and Unsworth (2013) definition corrected by
Schymanski et al. (2017), eq.22.

- `Rn` (W m-2): net radiation. Carefull: not the isothermal net radiation
- `VPD` (kPa): air vapor pressure deficit
- `γˢ` (kPa K−1): apparent value of psychrometer constant (see [`γ_star`](@ref))
- `Rbₕ` (s m-1): resistance for heat transfer by convection, i.e. resistance to sensible heat
- `Δ` (KPa K-1): rate of change of saturation vapor pressure with temperature (see [`e_sat_slope`](@ref))
- `ρ` (kg m-3): air density of moist air.
- `aₛₕ` (1,2): number of sides that exchange energy for heat (2 for leaves)
- `Cₚ` (J K-1 kg-1): specific heat of air for constant pressure

# References

Monteith, J. and Unsworth, M., 2013. Principles of environmental physics: plants, animals, and the atmosphere. Academic Press. See eq. 13.33.

Schymanski et al. (2017), Leaf-scale experiments reveal an important omission in the Penman–Monteith equation,
Hydrology and Earth System Sciences. DOI: https://doi.org/10.5194/hess-21-685-2017. See equ. 22.

# Examples

```julia
Tₐ = 20.0 ; P = 100.0 ;
ρ = air_density(Tₐ, P) # in kg m-3
Δ = e_sat_slope(Tₐ)

sensible_heat(300.0, 2.0, 0.1461683, 50.0, Δ, ρ, 2.0)
```
"""
function sensible_heat(Rn, VPD, γˢ, Rbₕ, Δ, ρ, aₛₕ, Cₚ)
  (γˢ * Rn - ρ * Cₚ * VPD * (aₛₕ / Rbₕ)) / (Δ + γˢ)
end

function sensible_heat(Rn, VPD, γˢ, Rbₕ, Δ, ρ, aₛₕ)
  sensible_heat(Rn, VPD, γˢ, Rbₕ, Δ, ρ, aₛₕ, Constants().Cₚ)
end


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

include("delays.jl")

using Actors, Test, .Delays

const fname = "test.x"

d = Dict(
    :test  => (4, 5, 6),
    :test1 => ("a", "b", "c")
    )

cp = checkpointing(1, fname)
checkpoint(cp, :test, 1,2,3)
@test restore(cp, :test) == (1,2,3)
checkpoint(cp, :test, 4,5,6)
@test restore(cp, :test) == (4,5,6)
checkpoint(cp, :test1, "a","b","c")

@test get_checkpoints(cp) == d
save_checkpoints(cp)
@test @delayed isfile(fname)
exit!(cp)
@test @delayed info(cp) == :done
cp = checkpointing(1, fname)
load_checkpoints(cp, fname)
@test get_checkpoints(cp) == d

rm(fname)


text:
# ellipsoid_im.jl

using MIRT: ellipsoid_im, shepp_logan_3d_parameters
using MIRT: image_geom
using MIRTjim: jim
import MIRT: _ellipsoid_im_check_fov

using Plots
using Test: @test, @test_throws

"""
    ellipsoid_im_show()
"""
function ellipsoid_im_show(ig)
	x1 = ellipsoid_im(ig, :spheroid)
	p1 = jim(x1, title="spheroid")

	x2 = ellipsoid_im(ig, :zhu; hu_scale=1000)
	p2 = jim(x2, title="zhu", clim=(900,1100))

	plot(p1, p2)
end

ig = image_geom(nx=512, nz=64*2, dz=0.625, fov=500)
ig = ig.down(8)

#test different ellipses
#@inferred # todo
ellipsoid_im(ig)
ellipsoid_im(ig, :zhu)
ellipsoid_im(ig, :kak, oversample=2)
ellipsoid_im(ig, :spheroid, return_params=true)
ellipsoid_im(ig; checkfov=true)
ellipsoid_im(ig; showmem=true)

@test_throws String ellipsoid_im(ig, :bad)
@test_throws String ellipsoid_im(ig, :spheroid, how=:bad)
@test_throws String shepp_logan_3d_parameters(0, 0, 0, :bad)

tmp = [1, 1, 1, 0, 0, 0]
@test !_ellipsoid_im_check_fov(2, 9, 9, [3 0 0 1 1 1 0 0 1], tmp...)
@test !_ellipsoid_im_check_fov(9, 2, 9, [0 3 0 1 1 1 0 0 1], tmp...)
@test !_ellipsoid_im_check_fov(9, 9, 2, [0 0 3 1 1 1 0 0 1], tmp...)

#ell1 = ellipsoid_im(ig, :zhu; how=:fast) # fast doesn't work
#ell2 = ellipsoid_im(ig, :zhu; how=:lowmem) # lowmem calls fast - doesn't work
# todo: test slow vs fast

ellipsoid_im_show(ig)


text:
module TensorProductQuadraturesTests

using Gridap.ReferenceFEs

using Test

degrees = (2,2)
quad = Quadrature(QUAD,tensor_product,degrees)
test_quadrature(quad)
@test sum(get_weights(quad)) ≈ 1

D = 2
degree = 3
quad = Quadrature(QUAD,tensor_product,degree)
test_quadrature(quad)
@test sum(get_weights(quad)) ≈ 1

@test isa(get_name(quad),String)

end # module


text:
export marginalrule

@marginalrule MvNormalMeanPrecision(:out_μ_Λ) (m_out::MvNormalMeanPrecision, m_μ::PointMass, m_Λ::PointMass) = begin
    return (out = prod(ProdAnalytical(), MvNormalMeanPrecision(mean(m_μ), mean(m_Λ)), m_out), μ = m_μ, Λ = m_Λ)
end

@marginalrule MvNormalMeanPrecision(:out_μ_Λ) (m_out::PointMass, m_μ::MvNormalMeanPrecision, m_Λ::PointMass) = begin
    return (out = m_out, μ = prod(ProdAnalytical(), m_μ, MvNormalMeanPrecision(mean(m_out), mean(m_Λ))), Λ = m_Λ)
end

@marginalrule MvNormalMeanPrecision(:out_μ) (m_out::MvNormalMeanPrecision, m_μ::MvNormalMeanPrecision, q_Λ::Any) = begin
    W_y  = invcov(m_out)
    xi_y = W_y * mean(m_out)

    W_m  = invcov(m_μ)
    xi_m = W_m * mean(m_μ)

    W_bar = mean(q_Λ)

    Λ  = [ W_y + W_bar -W_bar; -W_bar W_m + W_bar ]
    μ  = cholinv(Λ) * [ xi_y; xi_m ]

    return MvNormalMeanPrecision(μ, Λ)
end


text:
# The declaremissings imputor is sufficiently different in its input and behaviour
# that we don't bother using the ImputorTester here.
@testset "DeclareMissings" begin
    # List a couple known missing data values people might use.
    values = (NaN, 0.0, Nothing, "", 9999, -99, 0, DateTime(0))
    imp = DeclareMissings(; values=values)

    @testset "Vector" begin
        @testset "disallowmissing" begin
            a = collect(1.0:1.0:20.0)
            a[[2, 3, 7]] .= [NaN, 0.0, NaN]

            result = apply(a, imp)
            @test eltype(result) == Union{Float64, Missing}
            @test all(ismissing, result[[2, 3, 7]])

            # In-place operation don't work when the source array doesn't allow missings.
            b = copy(a)
            @test_throws MethodError apply!(b, imp)
        end

        @testset "allowmissing" begin
            a = allowmissing(collect(1.0:1.0:20.0))
            a[[2, 3, 7]] .= [NaN, 0.0, NaN]

            result = apply(a, imp)
            @test eltype(result) == Union{Float64, Missing}
            @test all(ismissing, result[[2, 3, 7]])

            # In-place operation work when the source array allows missings.
            b = copy(a)
            result2 = apply!(b, imp)
            @test eltype(result2) == Union{Float64, Missing}
            @test all(ismissing, result2[[2, 3, 7]])

            c = copy(a)
            result3 = Impute.declaremissings!(c; values=values)
            @test eltype(result3) == Union{Float64, Missing}
            @test all(ismissing, result3[[2, 3, 7]])
        end

        @testset "All missing" begin
            # Test having only missing data
            c = fill(missing, 10)
            @test isequal(apply(c, imp), c)
        end
    end

    @testset "Matrix" begin
        @testset "disallowmissing" begin
            a = collect(1.0:1.0:20.0)
            a[[2, 3, 7]] .= [NaN, 0.0, NaN]
            m = collect(reshape(a, 5, 4))

            result = apply(m, imp)
            @test eltype(result) == Union{Float64, Missing}
            @test all(ismissing, result[[2, 3, 7]])

            # In-place operation don't work when the source array doesn't allow missings.
            n = copy(m)
            @test_throws MethodError apply!(n, imp)
        end

        @testset "allowmissing" begin
            a = allowmissing(collect(1.0:1.0:20.0))
            a[[2, 3, 7]] .= [NaN, 0.0, NaN]
            m = collect(reshape(a, 5, 4))

            result = apply(m, imp)
            @test eltype(result) == Union{Float64, Missing}
            @test all(ismissing, result[[2, 3, 7]])

            # In-place operation don't work when the source array doesn't allow missings.
            n = copy(m)
            result2 = apply!(n, imp)
            @test eltype(result2) == Union{Float64, Missing}
            @test all(ismissing, result2[[2, 3, 7]])
        end

        @testset "All missing" begin
            # Test having only missing data
            c = fill(missing, 5, 4)
            @test isequal(apply(c, imp), c)
        end
    end
    @testset "Tables" begin
        @testset "DataFrame" begin
            table = DataFrame(
                :time => [DateTime(0), DateTime(2020, 02, 02), DateTime(2121, 12, 12)],
                :loc => [12, -5, 9999],
                :age => [4, -99, 18],
                :val => [1.5, NaN, 3.0],
                :desc => ["foo", "bar", ""],
            )
            mtable = DataFrame(
                :time => allowmissing([DateTime(0), DateTime(2020, 02, 02), DateTime(2121, 12, 12)]),
                :loc => allowmissing([12, -5, 9999]),
                :age => allowmissing([4, -99, 18]),
                :val => allowmissing([1.5, NaN, 3.0]),
                :desc => allowmissing(["foo", "bar", ""]),
            )
            expected = DataFrame(
                :time => [missing, DateTime(2020, 02, 02), DateTime(2121, 12, 12)],
                :loc => [12, -5, missing],
                :age => [4, missing, 18],
                :val => [1.5, missing, 3.0],
                :desc => ["foo", "bar", missing],
            )

            result = apply(table, imp)
            @test isequal(result, expected)
        end

        @testset "Column Table" begin
            table = (
                time = [DateTime(0), DateTime(2020, 02, 02), DateTime(2121, 12, 12)],
                loc = [12, -5, 9999],
                age = [4, -99, 18],
                val = [1.5, NaN, 3.0],
                desc = ["foo", "bar", ""],
            )
            mtable = (
                time = allowmissing([DateTime(0), DateTime(2020, 02, 02), DateTime(2121, 12, 12)]),
                loc = allowmissing([12, -5, 9999]),
                age = allowmissing([4, -99, 18]),
                val = allowmissing([1.5, NaN, 3.0]),
                desc = allowmissing(["foo", "bar", ""]),
            )
            expected = (
                time = [missing, DateTime(2020, 02, 02), DateTime(2121, 12, 12)],
                loc = [12, -5, missing],
                age = [4, missing, 18],
                val = [1.5, missing, 3.0],
                desc = ["foo", "bar", missing],
            )

            result = apply(table, imp)
            @test isequal(result, expected)
        end

        @testset "Row Table" begin
            table = Tables.rowtable((
                time = [DateTime(0), DateTime(2020, 02, 02), DateTime(2121, 12, 12)],
                loc = [12, -5, 9999],
                age = [4, -99, 18],
                val = [1.5, NaN, 3.0],
                desc = ["foo", "bar", ""],
            ))
            mtable = Tables.rowtable((
                time = allowmissing([DateTime(0), DateTime(2020, 02, 02), DateTime(2121, 12, 12)]),
                loc = allowmissing([12, -5, 9999]),
                age = allowmissing([4, -99, 18]),
                val = allowmissing([1.5, NaN, 3.0]),
                desc = allowmissing(["foo", "bar", ""]),
            ))
            expected = Tables.rowtable((
                time = [missing, DateTime(2020, 02, 02), DateTime(2121, 12, 12)],
                loc = [12, -5, missing],
                age = [4, missing, 18],
                val = [1.5, missing, 3.0],
                desc = ["foo", "bar", missing],
            ))

            result = apply(table, imp)
            @test isequal(result, expected)
        end
    end
end


text:
using DocPlotTest
using Documenter

DocMeta.setdocmeta!(DocPlotTest, :DocTestSetup, :(using DocPlotTest); recursive=true)

makedocs(;
    modules=[DocPlotTest],
    authors="Thomas Poulsen",
    repo="https://github.com/tp2750/DocPlotTest.jl/blob/{commit}{path}#{line}",
    sitename="DocPlotTest.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://tp2750.github.io/DocPlotTest.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Documenter plot" => "01_documenter-plot.md",
    ],
)

deploydocs(;
    repo="github.com/tp2750/DocPlotTest.jl",
    devbranch="main",
)


text:
#   This file is part of UnitSystems.jl
#   It is licensed under the MIT license
#   UnitSystems Copyright (C) 2020 Michael Reed
#       _           _                         _
#      | |         | |                       | |
#   ___| |__   __ _| | ___ __ __ ___   ____ _| | __ _
#  / __| '_ \ / _` | |/ / '__/ _` \ \ / / _` | |/ _` |
# | (__| | | | (_| |   <| | | (_| |\ V / (_| | | (_| |
#  \___|_| |_|\__,_|_|\_\_|  \__,_| \_/ \__,_|_|\__,_|
#
#   https://github.com/chakravala
#   https://crucialflow.com

@doc """
    μₑᵤ, μₚᵤ, μₚₑ, αinv, αG, ΩΛ

Physical measured dimensionless `Coupling` values with uncertainty are the electron to proton mass ratio `μₑᵤ`, proton to atomic mass ratio `μₚᵤ`, proton to electron mass ratio `μₚₑ`, inverted fine structure constant `αinv`, and the gravitaional coupling constant `αG`.

```Julia
julia> μₑᵤ # electronunit(Universe)
$μₑᵤ

julia> μₚᵤ # protonunit(Universe)
$μₚᵤ

julia> μₚₑ # protonelectron(Universe)
$μₚₑ

julia> αinv # 1/finestructure(Universe)
$αinv

julia> αG # coupling(Universe)
$αG

julia> ΩΛ # darkenergydensity(Universe)
$ΩΛ
```
""" Universe, μₑᵤ, μₚᵤ, μₚₑ, αinv, αG, meu, mpu, mpe, ainv, aG, ΩΛ, electronunit, protonunit, protonelectron, finestructure, coupling, darkenergydensity

@doc """
$(unitext(:turn,"2π/angle(U)"))

Complete rotation `angle` of revolution from a full circle.
```Julia
julia> turn(MetricEngineering)
$(turn(MetricEngineering))
```
""" turn

@doc """
$(unitext(:spat,"4π/solidangle(U)"))

Complete spherical `solidangle` around point from a full sphere.
```Julia
julia> spat(MetricEngineering)
$(spat(MetricEngineering))
```
""" spat

@doc """
$(unitext(:luminousefficacy,"Kcd*power(U)\nluminousefficacy(U::UnitSystem{𝟏}) = 𝟏"))

Luminous efficacy of monochromatic radiation `Kcd` of frequency 540 THz (lm⋅W⁻¹).
```Julia
julia> luminousefficacy(Metric) # lm⋅W⁻¹
$(luminousefficacy(Metric))

julia> luminousefficacy(CODATA) # lm⋅W⁻¹
$(luminousefficacy(CODATA))

julia> luminousefficacy(Conventional) # lm⋅W⁻¹
$(luminousefficacy(Conventional))

julia> luminousefficacy(International) # lm⋅W⁻¹
$(luminousefficacy(International))

julia> luminousefficacy(British) # lm⋅s³⋅slug⋅ft⁻²
$(luminousefficacy(British))
```
""" luminousefficacy, Kcd

@doc """
$(unitext(:molarmass,"avogadro(U)*electronmass(U)/electronunit(U)"))

Molar mass constant `Mᵤ` is the ratio of the `molarmass` and `relativemass` of a chemical.
```Julia
julia> molarmass(CGS) # g⋅mol⁻¹
$(molarmass(CGS))

julia> molarmass(Metric) # kg⋅mol⁻¹
$(molarmass(Metric))

julia> molarmass(SI2019) # kg⋅mol⁻¹
$(molarmass(SI2019))

julia> molarmass(International) # kg⋅mol⁻¹
$(molarmass(International))
```
""" molarmass, Mᵤ, Mu

@doc """
$(unitext(:avogadro, "molargas(x)/boltzmann(x) # Mᵤ/atomicmass(x)"))

Avogadro `NA` is `molarmass(x)/atomicmass(x)` number of atoms in a 12 g sample of C₁₂.
```Julia
julia> avogadro(SI2019) # mol⁻¹
$(avogadro(SI2019))

julia> avogadro(Metric) # mol⁻¹
$(avogadro(Metric))

julia> avogadro(CODATA) # mol⁻¹
$(avogadro(CODATA))

julia> avogadro(Conventional) # mol⁻¹
$(avogadro(Conventional))

julia> avogadro(English) # lb-mol⁻¹
$(avogadro(English))

julia> avogadro(British) # slug-mol⁻¹
$(avogadro(British))
```
""" avogadro, NA

@doc """
$(unitext(:planckreduced,"planck(x)/turn(x)"))

Reduced Planck constant `ħ` is a Planck per radian (J⋅s⋅rad⁻¹ or ft⋅lb⋅s⋅rad⁻¹).

```Julia
julia> planckreduced(SI2019) # J⋅s⋅rad⁻¹
$(planckreduced(SI2019))

julia> planckreduced(SI2019)*lightspeed(SI2019) # J⋅m⋅rad⁻¹
$(planckreduced(SI2019)*lightspeed(SI2019))

julia> planckreduced(CODATA) # J⋅s⋅rad⁻¹
$(planckreduced(CODATA))

julia> planckreduced(Conventional) # J⋅s⋅rad⁻¹
$(planckreduced(Conventional))

julia> planckreduced(SI2019)/elementarycharge(SI2019) # eV⋅s⋅rad⁻¹
$(planckreduced(SI2019)/elementarycharge(SI2019))

julia> planckreduced(SI2019)*lightspeed(SI2019)/elementarycharge(SI2019) # eV⋅m⋅rad⁻¹
$(planckreduced(SI2019)*lightspeed(SI2019)/elementarycharge(SI2019))

julia> planckreduced(British) # ft⋅lb⋅s⋅rad⁻¹
$(planckreduced(British))
```
""" planckreduced, ħ

@doc """
$(unitext(:planck,"turn(x)*planckreduced(x)"))

Planck constant `𝘩` is energy per electromagnetic frequency (J⋅s or ft⋅lb⋅s).

```Julia
julia> planck(SI2019) # J⋅s
$(planck(SI2019))

julia> planck(SI2019)*lightspeed(SI2019) # J⋅m
$(planck(SI2019)*lightspeed(SI2019))

julia> planck(CODATA) # J⋅s
$(planck(CODATA))

julia> planck(Conventional) # J⋅s
$(planck(Conventional))

julia> planck(SI2019)/elementarycharge(SI2019) # eV⋅s
$(planck(SI2019)/elementarycharge(SI2019))

julia> planck(SI2019)*lightspeed(SI2019)/elementarycharge(SI2019) # eV⋅m
$(planck(SI2019)*lightspeed(SI2019)/elementarycharge(SI2019))

julia> planck(British) # ft⋅lb⋅s
$(planck(British))
```
""" planck, 𝘩, hh

@doc """
$(unitext(:boltzmann,"molargas(x)/avogadro(x)"))

Boltzmann constant `kB` is the entropy amount of a unit number microstate permutation.

```Julia
julia> boltzmann(SI2019) # J⋅K⁻¹
$(boltzmann(SI2019))

julia> boltzmann(Metric) # J⋅K⁻¹
$(boltzmann(Metric))

julia> boltzmann(SI2019)/elementarycharge(SI2019) # eV⋅K⁻¹
$(boltzmann(SI2019)/elementarycharge(SI2019))

julia> boltzmann(SI2019)/planck(SI2019) # Hz⋅K⁻¹
$(boltzmann(SI2019)/planck(SI2019))

julia> boltzmann(CGS) # erg⋅K⁻¹
$(boltzmann(CGS))

julia> boltzmann(SI2019)/calorie(SI2019) # calᵢₜ⋅K⁻¹
$(boltzmann(SI2019)/calorie(SI2019))

julia> boltzmann(SI2019)*°R/calorie(SI2019) # calᵢₜ⋅°R⁻¹
$(boltzmann(SI2019)*°R/calorie(SI2019))

julia> boltzmann(Brtish) # ft⋅lb⋅°R⁻¹
$(boltzmann(British))

julia> boltzmann(SI2019)/planck(SI2019)/lightspeed(SI2019) # m⁻¹⋅K⁻¹
$(boltzmann(SI2019)/planck(SI2019)/lightspeed(SI2019))

julia> avogadro(SI2019)*boltzmann(SI2019)/calorie(SI2019) # calᵢₜ⋅mol⁻¹⋅K⁻¹
$(avogadro(SI2019)*boltzmann(SI2019)/calorie(SI2019))

julia> dB(boltzmann(SI2019)) # dB(W⋅K⁻¹⋅Hz⁻¹)
$(dB(boltzmann(SI2019)))
```
""" boltzmann, kB

@doc """
$(unitext(:lightspeed,"𝟏/sqrt(vacuumpermeability(U)*vacuumpermittivity(U))/lorentz(U)"))

Speed of light in a vacuum `𝘤` for massless particles (m⋅s⁻¹ or ft⋅s⁻¹).

```Julia
julia> lightspeed(Metric) # m⋅s⁻¹
$(lightspeed(Metric))

julia> lightspeed(English) # ft⋅s⁻¹
$(lightspeed(English))

julia> lightspeed(IAU) # au⋅D⁻¹
$(lightspeed(IAU))
```
""" lightspeed, 𝘤, cc

@doc """
$(unitext(:vacuumpermeability,"𝟏/vacuumpermittivity(U)/(lightspeed(U)*lorentz(U))^2"))

Magnetic permeability in a classical vacuum defined as `μ₀` in SI units (H⋅m⁻¹, kg⋅m²⋅C⁻²).

```Julia
julia> vacuumpermeability(Metric) # H⋅m⁻¹
$(vacuumpermeability(Metric))

julia> vacuumpermeability(Conventional) # H⋅m⁻¹
$(vacuumpermeability(Conventional))

julia> vacuumpermeability(CODATA) # H⋅m⁻¹
$(vacuumpermeability(CODATA))

julia> vacuumpermeability(SI2019) # H⋅m⁻¹
$(vacuumpermeability(SI2019))

julia> vacuumpermeability(International) # H⋅m⁻¹
$(vacuumpermeability(International))

julia> vacuumpermeability(EMU) # abH⋅cm⁻¹
$(vacuumpermeability(EMU))

julia> vacuumpermeability(ESU) # statH⋅cm⁻¹
$(vacuumpermeability(ESU))
```
""" vacuumpermeability, μ₀, m0

@doc """
$(unitext(:lorentz,"spat(U)*biotsavart(U)/vacuumpermeability(U)/rationalization(U)"))

Electromagnetic proportionality constant `αL` for the Lorentz's law force (dimensionless).

```Julia
julia> lorentz(Metric)
$(lorentz(Metric))

julia> lorentz(LorentzHeaviside)
$(lorentz(LorentzHeaviside))

julia> lorentz(Gauss)
$(lorentz(Gauss))
```
""" lorentz, αL, aL, C

@doc """
$(unitext(:rationalization,"spat(U)*biotsavart(U)/vacuumpermeability(U)/lorentz(U)"))

Constant of magnetization and polarization density or `spat(U)*coulomb(U)*permittivity(U)`.

```Julia
julia> rationalization(Metric)
$(rationalization(Metric))

julia> rationalization(Gauss)
$(rationalization(Gauss))
```
""" rationalization, Λ

@doc """
$(unitext(:electronmass,"protonmass(U)/protonelectron(U) # αinv^2*R∞*2𝘩/𝘤"))

Electron rest mass `mₑ` of subatomic particle with `-𝘦` elementary charge  (kg or slugs).
```Julia
julia> electronmass(Metric) # kg
$(electronmass(Metric))

julia> electronmass(CODATA) # kg
$(electronmass(CODATA))

julia> electronmass(Conventional) # kg
$(electronmass(Conventional))

julia> electronmass(International) # kg
$(electronmass(International))

julia> electronmass(Metric)/atomicmass(Metric) # Da
$μₑᵤ

julia> electronmass(Metric)*lightspeed(Metric)^2 # J
$(electronmass(Metric)*lightspeed(Metric)^2)

julia> electronmass(SI2019)*lightspeed(SI2019)^2/elementarycharge(SI2019) # eV⋅𝘤⁻²
$(electronmass(SI2019)*lightspeed(SI2019)^2/elementarycharge(SI2019))

julia> electronmass(English) # lb
$(electronmass(English))
```
""" electronmass, mₑ, me

@doc """
$(unitext(:atomicmass,"Mᵤ/avogadro(U) # $(molarmass(SI2019)) ≈ 0.001-3.5e-13"))

Atomic mass unit `mᵤ` of 1/12 of the C₁₂ carbon-12 atom's mass  (kg or slugs).
```Julia
julia> atomicmass(Metric) # kg
$(atomicmass(Metric))

julia> atomicmass(Metric)/electronmass(Metric) # mₑ
$(atomicmass(Metric)/electronmass(Metric))

julia> atomicmass(Metric)*lightspeed(Metric)^2 # J
$(atomicmass(Metric)*lightspeed(Metric)^2)

julia> atomicmass(SI2019)*lightspeed(SI2019)^2/elementarycharge(SI2019) # eV⋅𝘤⁻²
$(atomicmass(SI2019)*lightspeed(SI2019)^2/elementarycharge(SI2019))

julia> atomicmass(British) # lb
$(atomicmass(British))
```
""" atomicmass, mᵤ, mu

@doc """
$(unitext(:protonmass,"protonunit(U)*atomicmass(U)"))

Proton mass `mₚ` of subatomic particle with `+𝘦` elementary charge  (kg or mass).
```Julia
julia> protonmass(Metric) # kg
$(protonmass(Metric))

julia> protonmass(SI2019)*lightspeed(SI2019)^2/elementarycharge(SI2019) # eV⋅𝘤⁻²
$(protonmass(SI2019)*lightspeed(SI2019)^2/elementarycharge(SI2019))

julia> protonmass(Metric)/atomicmass(Metric) # mᵤ
$(protonmass(Metric)/atomicmass(Metric))

julia> protonmass(Metric)/electronmass(Metric) # mₑ
$(protonmass(Metric)/electronmass(Metric))
```
""" protonmass, mₚ, mp

@doc """
$(unitext(:planckmass,"electronmass(U)/sqrt(coupling(U))"))

Planck mass factor `mP` from the gravitational coupling constant `αG` (kg or slugs).
```Julia
juila> planckmass(Metric)*lightspeed(Metric)^2/elementarycharge(Metric) # eV⋅𝘤⁻²
$(planckmass(Metric)*lightspeed(Metric)^2/elementarycharge(Metric))

juila> planckmass(Metric) # kg
$(planckmass(Metric))

juila> planckmass(Metric)/atomicmass(Metric) # mᵤ
$(planckmass(Metric)/atomicmass(Metric))

juila> planckmass(Metric)*lightspeed(Metric)^2/elementarycharge(Metric)/sqrt(𝟐^2*τ) # eV⋅𝘤⁻²
$(planckmass(Metric)*lightspeed(Metric)^2/elementarycharge(Metric)/sqrt(𝟐^2*τ))

juila> planckmass(Metric)/sqrt(𝟐^2*τ) # kg
$(planckmass(Metric)/sqrt(𝟐^2*τ))
```
""" planckmass, mP

@doc """
$(unitext(:gravitation,"lightspeed(U)*planckreduced(U)/planckmass(U)^2"))

Universal gravitational constant `G` of Newton's law (m³⋅kg⁻¹⋅s⁻² or ft³⋅slug⁻¹⋅s⁻²).
```Julia
juila> gravitation(Metric) # m³⋅kg⁻¹⋅s⁻²
$(gravitation(Metric))

julia> gravitation(English) # ft³⋅lbm⁻¹⋅s⁻²
$(gravitation(English))

julia> gravitation(PlanckGauss)
$(gravitation(PlanckGauss))
```
""" gravitation, G, GG

@doc """
$(unitext(:einstein,"𝟐*spat(U)*gravitation(U)/lightspeed(U)^4"))

Einstein's gravitational constant from the Einstein field equations (s⋅²⋅m⁻¹⋅kg⁻¹).
```Julia
julia> einstein(Metric) # s²⋅m⁻¹⋅kg⁻¹
$(einstein(Metric))

julia> einstein(IAU) # day²⋅au⁻¹⋅M☉⁻¹
$(einstein(IAU))
```
""" einstein, κ

#=@doc """
$(unitext(:einstein2,"𝟐*spat(U)*gravitation(U)/lightspeed(U)^2"))

Einstein's gravitational constant from the Einstein field equations (m⋅kg⁻¹).
```Julia
julia> einstein2(Metric) # m⋅kg⁻¹
$(einstein2(Metric))

julia> einstein2(IAU) # au⋅M☉⁻¹
$(einstein2(IAU))
```
""" einstein2=#

@doc """
$(unitext(:gravity,"# mass*acceleration/force"))

Gravitational force reference used in technical engineering units (kg⋅m⋅N⁻¹⋅s⁻²).
```Julia
julia> gravity(Metric)
$(gravity(Metric))

julia> gravity(MetricEngineering) # m⋅kg⋅N⁻¹⋅s⁻²
$(gravity(MetricEngineering))

julia> gravity(English) # ft⋅lbm⋅lbf⁻¹⋅s⁻²
$(gravity(English))
```
""" gravity

@doc """
$(unitext(:molargas,"boltzmann(x)*avogadro(x)"))

Universal gas constant `Rᵤ` is factored into specific `gasconstant(x)*molarmass(x)` values.

```Julia
julia> molargas(SI2019) # J⋅K⁻¹⋅mol⁻¹
$(molargas(SI2019))

julia> molargas(English)/𝟐^4/𝟑^2 # psi⋅ft³⋅°R⁻¹⋅lb-mol⁻¹
$(molargas(English)/𝟐^4/𝟑^2)

julia> molargas(English)/atmosphere(English) # atm⋅ft³⋅R⁻¹⋅lb-mol⁻¹
$(molargas(English)/atmosphere(English))

julia> molargas(English)/thermalunit(English) # BTU⋅°R⁻¹⋅lb-mol⁻¹
$(molargas(English)/thermalunit(English))

julia> molargas(Metric)/atmosphere(Metric) # atm⋅m³⋅K⁻¹⋅mol⁻¹
$(molargas(Metric)/atmosphere(Metric))

julia> molargas(Metric)/torr(Metric) # m³⋅torr⋅K⁻¹⋅mol⁻¹
$(molargas(Metric)/torr(Metric))

julia> molargas(English)/torr(English) # ft³⋅torr⋅°R⁻¹⋅lb-mol⁻¹
$(molargas(English)/torr(English))

julia> molargas(CGS) # erg⋅K⁻¹⋅mol⁻¹
$(molargas(CGS))

julia> molargas(English) # ft⋅lb⋅°R⁻¹⋅lb-mol⁻¹
$(molargas(English))

julia> molargas(British) # ft⋅lb⋅°R⁻¹⋅slug-mol⁻¹
$(molargas(British))

julia> molargas(SI1976) # J⋅K⁻¹⋅mol⁻¹ (US1976 Standard Atmosphere)
$(molargas(SI1976))
```
""" molargas, Rᵤ, Ru

@doc """
$(unitext(:stefan,"π^4/𝟐*spat(U)*boltzmann(U)^4/(15planck(U)^3*lightspeed(U)^2)"))

Stefan-Boltzmann proportionality `σ` of black body radiation (W⋅m⁻²⋅K⁻⁴ or ?⋅ft⁻²⋅°R⁻⁴).

```Julia
julia> stefan(SI2019) # W⋅m⁻²⋅K⁻⁴
$(stefan(SI2019))

julia> stefan(Metric) # W⋅m⁻²⋅K⁻⁴
$(stefan(Metric))

julia> stefan(Conventional) # W⋅m⁻²⋅K⁻⁴
$(stefan(Conventional))

julia> stefan(CODATA) # W⋅m⁻²⋅K⁻⁴
$(stefan(CODATA))

julia> stefan(Metric)*day(Metric)/(calorie(Metric)*100^2) # cal⋅cm⁻²⋅day⁻¹⋅K⁻⁴
$(stefan(Metric)*day(Metric)/calorie(Metric))

julia> stefan(English) # lb⋅s⁻¹⋅ft⁻³⋅°R⁻⁴
$(stefan(English))
```
""" stefan, σ, SB

@doc """
$(unitext(:radiationdensity,"𝟐^2*stefan(U)/lightspeed(U)"))

Raditation density constant of black body radiation (J⋅m⁻³⋅K⁻⁴ or lb⋅ft⁻²⋅°R⁻⁴).

```Julia
julia> radiationdensity(Metric) # J⋅m⁻³⋅K⁻⁴
$(radiationdensity(Metric))

julia> radiationdensity(SI2019) # J⋅m⁻³⋅K⁻⁴
$(radiationdensity(SI2019))

julia> radiationdensity(Conventional) # J⋅m⁻³⋅K⁻⁴
$(radiationdensity(Conventional))

julia> radiationdensity(CODATA) # J⋅m⁻³⋅K⁻⁴
$(radiationdensity(CODATA))

julia> radiationdensity(International) # J⋅m⁻³⋅K⁻⁴
$(radiationdensity(International))
```
""" radiationdensity

@doc """
$(unitext(:vacuumpermittivity,"𝟏/vacuumpermeability(U)/(lightspeed(U)*lorentz(U))^2"))

Dielectric permittivity constant `ε₀` of a classical vacuum (C²⋅N⁻¹⋅m⁻²).

```Julia
julia> vacuumpermittivity(Metric) # F⋅m⁻¹
$(vacuumpermittivity(Metric))

julia> vacuumpermittivity(Conventional) # F⋅m⁻¹
$(vacuumpermittivity(Conventional))

julia> vacuumpermittivity(CODATA) # F⋅m⁻¹
$(vacuumpermittivity(CODATA))

julia> vacuumpermittivity(SI2019) # F⋅m⁻¹
$(vacuumpermittivity(SI2019))

julia> vacuumpermittivity(International) # F⋅m⁻¹
$(vacuumpermittivity(International))

julia> vacuumpermittivity(EMU) # abF⋅cm⁻¹
$(vacuumpermittivity(EMU))

julia> vacuumpermittivity(ESU) # statF⋅cm⁻¹
$(vacuumpermittivity(ESU))

julia> vacuumpermittivity(SI2019)/elementarycharge(SI2019) # 𝘦²⋅eV⁻¹⋅m⁻¹
$(vacuumpermittivity(SI2019)/elementarycharge(SI2019))
```
""" vacuumpermittivity, ε₀, ϵ₀, e0

@doc """
$(unitext(:electrostatic,"rationalization(U)/spat(U)/vacuumpermittivity(U)"))

Electrostatic proportionality constant `kₑ` for the Coulomb's law force (N⋅m²⋅C⁻²).

```Julia
julia> electrostatic(Metric) # N⋅m²⋅C⁻²
$(electrostatic(Metric))

julia> electrostatic(CODATA) # N·m²⋅C⁻²
$(electrostatic(CODATA))

julia> electrostatic(SI2019) # N·m²⋅C⁻²
$(electrostatic(SI2019))

julia> electrostatic(Conventional) # N·m²⋅C⁻²
$(electrostatic(Conventional))

julia> electrostatic(International) # N·m²⋅C⁻²
$(electrostatic(International))

julia> electrostatic(EMU) # dyn⋅cm²⋅abC⁻²
$(electrostatic(EMU))

julia> electrostatic(ESU) # dyn⋅cm²⋅statC⁻²
$(electrostatic(ESU))

julia> electrostatic(HLU) # dyn⋅cm²⋅hlC⁻²
$(electrostatic(HLU))
```
""" electrostatic, kₑ, ke

@doc """
$(unitext(:biotsavart,"vacuumpermeability(U)*lorentz(U)*rationalization(U)/spat(U)"))

Magnetostatic proportionality constant `αB` for the Biot-Savart's law (H/m).

```Julia
julia> biotsavart(Metric) # H⋅m⁻¹
$(biotsavart(Metric))

julia> biotsavart(CODATA) # H⋅m⁻¹
$(biotsavart(CODATA))

julia> biotsavart(SI2019) # H⋅m⁻¹
$(biotsavart(SI2019))

julia> biotsavart(Conventional) # H⋅m⁻¹
$(biotsavart(Conventional))

julia> biotsavart(International) # H⋅m⁻¹
$(biotsavart(International))

julia> biotsavart(InternationalMean) # H⋅m⁻¹
$(biotsavart(InternationalMean))

julia> biotsavart(EMU) # abH⋅cm⁻¹
$(biotsavart(EMU))

julia> biotsavart(ESU) # statH⋅cm⁻¹
$(biotsavart(ESU))

julia> biotsavart(Gauss) # abH⋅cm⁻¹
$(biotsavart(Gauss))

julia> biotsavart(HLU) # hlH⋅cm⁻¹
$(biotsavart(HLU))
```
""" biotsavart, αB, aB

@doc """
$(unitext(:magnetostatic,"lorentz(U)*biotsavart(U) # electrostatic(U)/lightspeed(U)^2"))

Magnetic proportionality constant `kₘ` for the Ampere's law force (N·s²⋅C⁻²).

```Julia
julia> magnetostatic(Metric) # H⋅m⁻¹
$(magnetostatic(Metric))

julia> magnetostatic(CODATA) # H⋅m⁻¹
$(magnetostatic(CODATA))

julia> magnetostatic(SI2019) # H⋅m⁻¹
$(magnetostatic(SI2019))

julia> magnetostatic(Conventional) # H⋅m⁻¹
$(magnetostatic(Conventional))

julia> magnetostatic(International) # H⋅m⁻¹
$(magnetostatic(International))

julia> magnetostatic(EMU) # abH⋅m⁻¹
$(magnetostatic(EMU))

julia> magnetostatic(ESU) # statH⋅m⁻¹
$(magnetostatic(ESU))

julia> magnetostatic(HLU) # hlH⋅m⁻¹
$(magnetostatic(HLU))
```
""" magnetostatic, kₘ, km

@doc """
$(unitext(:vacuumimpedance,"vacuumpermeability(U)*lightspeed(U)*rationalization(U)*lorentz(U)^2"))

Vacuum impedance of free space `Z₀` is magnitude ratio of electric to magnetic field (Ω).
```Julia
julia> vacuumimpedance(Metric) # Ω
$(vacuumimpedance(Metric))

julia> vacuumimpedance(Conventional) # Ω
$(vacuumimpedance(Conventional))

julia> vacuumimpedance(CODATA) # Ω
$(vacuumimpedance(CODATA))

julia> vacuumimpedance(SI2019) # Ω
$(vacuumimpedance(SI2019))

julia> vacuumimpedance(International) # Ω
$(vacuumimpedance(International))

julia> vacuumimpedance(InternationalMean) # Ω
$(vacuumimpedance(InternationalMean))

julia> 120π # 3e8*μ₀ # Ω
$(120π)

julia> vacuumimpedance(EMU) # abΩ
$(vacuumimpedance(EMU))

julia> vacuumimpedance(ESU) # statΩ
$(vacuumimpedance(ESU))

julia> vacuumimpedance(HLU) # hlΩ
$(vacuumimpedance(HLU))
```
""" vacuumimpedance, Z₀, Z0

@doc """
$(unitext(:elementarycharge,"√(𝟐*planck(U)*finestructure(U)/vacuumimpedance(U))"))

Quantized elementary charge `𝘦` of a proton or electron `2/(klitzing(U)*josephson(U))` (C).
```Julia
julia> elementarycharge(SI2019) # C
$(elementarycharge(SI2019))

julia> elementarycharge(Metric) # C
$(elementarycharge(Metric))

julia> elementarycharge(CODATA) # C
$(elementarycharge(CODATA))

julia> elementarycharge(Conventional) # C
$(elementarycharge(Conventional))

julia> elementarycharge(International) # C
$(elementarycharge(International))

julia> elementarycharge(EMU) # abC
$(elementarycharge(EMU))

julia> elementarycharge(ESU) # statC
$(elementarycharge(ESU))

julia> elementarycharge(Planck) # sqrt(4π/αinv)
$(elementarycharge(Planck))
```
""" elementarycharge, 𝘦, ee

@doc """
$(unitext(:faraday,"elementarycharge(U)*avogadro(U)"))

Electric charge per mole of electrons `𝔉` based on elementary charge (C⋅mol⁻¹).
```Julia
julia> faraday(SI2019) # C⋅mol⁻¹
$(faraday(SI2019))

julia> faraday(Metric) # C⋅mol⁻¹
$(faraday(Metric))

julia> faraday(CODATA) # C⋅mol⁻¹
$(faraday(CODATA))

julia> faraday(Conventional) # C⋅mol⁻¹
$(faraday(Conventional))

julia> faraday(International) # C⋅mol⁻¹
$(faraday(International))

julia> faraday(InternationalMean) # C⋅mol⁻¹
$(faraday(InternationalMean))

julia> faraday(EMU) # abC⋅mol⁻¹
$(faraday(EMU))

julia> faraday(ESU) # statC⋅mol⁻¹
$(faraday(ESU))

julia> faraday(Metric)/kilocalorie(Metric) # kcal⋅(V-g-e)⁻¹
$(faraday(Metric)/kilocalorie(Metric))

julia> faraday(Metric)/3600 # A⋅h⋅mol⁻¹
$(faraday(Metric)/HOUR)
```
""" faraday, 𝔉, FF

@doc """
$(unitext(:josephson,"𝟐*elementarycharge(U)*lorentz(U)/planck(U) # 𝟏/magneticfluxquantum(U)"))

Josephson constant `KJ` relating potential difference to irradiation frequency (Hz⋅V⁻¹).
```Julia
julia> josephson(SI2019) # Hz⋅V⁻¹
$(josephson(SI2019))

julia> josephson(Metric) # Hz⋅V⁻¹
$(josephson(Metric))

julia> josephson(Conventional) # Hz⋅V⁻¹
$(josephson(Conventional))

julia> josephson(CODATA) # Hz⋅V⁻¹
$(josephson(CODATA))

julia> josephson(International) # Hz⋅V⁻¹
$(josephson(International))

julia> josephson(EMU) # Hz⋅abV⁻¹
$(josephson(EMU))

julia> josephson(ESU) # Hz⋅statV⁻¹
$(josephson(ESU))
```
""" josephson, KJ

@doc """
$(unitext(:magneticfluxquantum,"planck(U)/𝟐/elementarycharge(U)/lorentz(U)"))

Magnetic flux quantum `Φ₀` is `𝟏/josephson(U)` (Wb).
```Julia
julia> magneticfluxquantum(SI2019) # Wb
$(magneticfluxquantum(SI2019))

julia> magneticfluxquantum(Metric) # Wb
$(magneticfluxquantum(Metric))

julia> magneticfluxquantum(Conventional) # Wb
$(magneticfluxquantum(Conventional))

julia> magneticfluxquantum(International) # Wb
$(magneticfluxquantum(International))

julia> magneticfluxquantum(InternationalMean) # Wb
$(magneticfluxquantum(InternationalMean))

julia> magneticfluxquantum(EMU) # Mx
$(magneticfluxquantum(EMU))

julia> magneticfluxquantum(ESU) # statWb
$(magneticfluxquantum(ESU))
```
""" magneticfluxquantum, Φ₀

@doc """
$(unitext(:klitzing,"planck(U)/elementarycharge(U)^2"))

Quantized Hall resistance `RK` (Ω).
```Julia
julia> klitzing(SI2019) # Ω
$(klitzing(SI2019))

julia> klitzing(Metric) # Ω
$(klitzing(Metric))

julia> klitzing(Conventional) # Ω
$(klitzing(Conventional))

julia> klitzing(International) # Ω
$(klitzing(International))

julia> klitzing(CODATA) # Ω
$(klitzing(CODATA))

julia> klitzing(EMU) # abΩ
$(klitzing(EMU))

julia> klitzing(ESU) # statΩ
$(klitzing(ESU))
```
""" klitzing, RK

@doc """
$(unitext(:conductancequantum,"𝟐*elementarycharge(U)^2/planck(U) # 2/klitzing(U)"))

Conductance quantum `G₀` is a quantized unit of electrical conductance (S).
```Julia
julia> conductancequantum(SI2019) # S
$(conductancequantum(SI2019))

julia> conductancequantum(Metric) # S
$(conductancequantum(Metric))

julia> conductancequantum(Conventional) # S
$(conductancequantum(Conventional))

julia> conductancequantum(CODATA) # S
$(conductancequantum(CODATA))

julia> conductancequantum(International) # S
$(conductancequantum(International))

julia> conductancequantum(InternationalMean) # S
$(conductancequantum(InternationalMean))

julia> conductancequantum(EMU) # abS
$(conductancequantum(EMU))

julia> conductancequantum(ESU) # statS
$(conductancequantum(ESU))
```
""" conductancequantum, G₀, G0

@doc """
$(unitext(:hartree,"electronmass(U)/gravity(U)*(lightspeed(U)*finestructure(U))^2"))

Hartree electric potential energy `Eₕ` of the hydrogen atom at ground state is `2R∞*𝘩*𝘤` (J).
```Julia
julia> hartree(SI2019)/elementarycharge(SI2019) # eV
$(hartree(SI2019)/elementarycharge(SI2019))

julia> hartree(Metric) # J
$(hartree(Metric))

julia> hartree(CGS) # erg
$(hartree(CGS))

julia> hartree(Metric)*avogadro(Metric)/kilo # kJ⋅mol⁻¹
$(hartree(Metric)*avogadro(Metric)/kilo)

julia> hartree(Metric)*avogadro(Metric)/kilocalorie(Metric) # kcal⋅mol⁻¹
$(hartree(Metric)*avogadro(Metric)/kilocalorie(Metric))

julia> 𝟐*centi*rydberg(Metric) # Eₕ/𝘩/𝘤/100 cm⁻¹
$(hartree(Metric)/planck(Metric)/lightspeed(Metric)/(𝟐*𝟓)^2)

julia> hartree(Metric)/planck(Metric) # Hz
$(hartree(Metric)/planck(Metric))

julia> hartree(Metric)/boltzmann(Metric) # K
$(hartree(Metric)/boltzmann(Metric))
```
In a Gaussian unit system where `4π*ε₀ == 1` the Hartree energy is `𝘦^2/a₀`.
""" hartree, Eₕ, Eh

@doc """
$(unitext(:rydberg,"hartree(U)/2planck(U)/lightspeed(U) # Eₕ/2𝘩/𝘤"))

Rydberg constant `R∞` is lowest energy photon capable of ionizing atom at ground state (m⁻¹).
```Julia
julia> rydberg(Metric) # m⁻¹
$(rydberg(Metric))
```
The Rydberg constant for hydrogen `RH` is `R∞*mₚ/(mₑ+mₚ)` (m⁻¹).
```Julia
julia> rydberg(Metric)*protonmass(Metric)/(electronmass(Metric)+protonmass(Metric)) # m⁻¹
$(rydberg(Metric)*protonmass(Metric)/(electronmass(Metric)+protonmass(Metric)))
```
Rydberg unit of photon energy `Ry` is `𝘩*𝘤*R∞` or `Eₕ/2` (J).
```Julia
julia> hartree(Metric)/2 # J
$(hartree(Metric)/𝟐)

julia> hartree(SI2019)/𝟐/elementarycharge(SI2019) # eV
$(hartree(SI2019)/𝟐/elementarycharge(SI2019))
```
Rydberg photon frequency `𝘤*R∞` or `Eₕ/2𝘩` (Hz).
```Julia
julia> lightspeed(Metric)*rydberg(Metric) # Hz
$(lightspeed(Metric)*rydberg(Metric))
```
Rydberg wavelength `1/R∞` (m).
```Julia
julia> 𝟏/rydberg(Metric) # m
$(𝟏/rydberg(Metric))

julia> 𝟏/rydberg(Metric)/τ # m⋅rad⁻¹
$(𝟏/rydberg(Metric)/τ)
```
Precision measurements of the Rydberg constants are within a relative standard uncertainty of under 2 parts in 10¹², and is chosen to constrain values of other physical constants.
""" rydberg, R∞, RH, Ry

@doc """
$(unitext(:bohr,"planckreduced(U)*gravity(U)/electronmass(U)/lightspeed(U)/finestructure(U)"))

Bohr radius of the hydrogen atom in its ground state `a₀` (m).
```Julia
julia> bohr(Metric) # m
$(bohr(Metric))

julia> bohr(IPS) # in
$(bohr(IPS))

julia> bohr(Hartree) # a₀
$(bohr(Hartree))
```
""" bohr, a₀, a0
#julia> bohr(Metric)/length(PlanckGauss) # ℓP
#$(bohr(Metric)/length(PlanckGauss))

#=@doc """
$(unitext(:bohrreduced,"bohr(U)*(1+1/protonelectron(U))"))

Reduced Bohr radius including the effect of reduced mass in hydrogen atom (m).
```Julia
julia> bohrreduced(Metric) # m
$(bohrreduced(Metric))

julia> bohrreduced(Metric)/bohr(Metric) # a₀
$(bohrreduced(Metric)/bohr(Metric))
```
""" bohrreduced=#

@doc """
$(unitext(:electronradius,"finestructure(U)*planckreduced(U)*gravity(U)/electronmass(U)/lightspeed(U)"))

Classical electron radius or Lorentz radius or Thomson scattering length (m).
```Julia
julia> electronradius(Metric) # m
$(electronradius(Metric))

julia> electronradius(CODATA) # m
$(electronradius(CODATA))

julia> electronradius(Conventional) # m
$(electronradius(Conventional))
```
""" electronradius, rₑ, re

@doc """
$(unitext(:magneton,"elementarycharge(U)*planckreduced(U)*lorentz(U)/2electronmass(U)"))

Bohr magneton `μB` natural unit for expressing magnetic moment of electron (J⋅T⁻¹).
```Julia
julia> magneton(SI2019) # J⋅T⁻¹
$(magneton(SI2019))

julia> magneton(Metric) # J⋅T⁻¹
$(magneton(Metric))

julia> magneton(CODATA) # J⋅T⁻¹
$(magneton(CODATA))

julia> magneton(Conventional) # J⋅T⁻¹
$(magneton(Conventional))

julia> magneton(International) # J⋅T⁻¹
$(magneton(International))

julia> magneton(ESU) # statA⋅cm²
$(magneton(ESU))

julia> magneton(SI2019)/elementarycharge(SI2019) # eV⋅T⁻¹
$(magneton(SI2019)/elementarycharge(SI2019))

julia> magneton(Hartree) # 𝘤⋅ħ⋅mₑ⁻¹
$(magneton(Hartree))
```
""" magneton, μB

@doc """
$(unitext(:wienwavelength,"planck(U)*lightspeed(U)/boltzmann(U)/(𝟓+W₀(-𝟓*exp(-𝟓)))"))

Wien wavelength displacement law constant based on Lambert `W₀` evaluation (m⋅K or ft⋅°R).
```Julia
julia> wienwavelength(Metric) # m⋅K
$(wienwavelength(Metric))

julia> wienwavelength(SI2019) # m⋅K
$(wienwavelength(SI2019))

julia> wienwavelength(Conventional) # m⋅K
$(wienwavelength(Conventional))

julia> wienwavelength(CODATA) # m⋅K
$(wienwavelength(CODATA))

julia> wienwavelength(English) # ft⋅°R
$(wienwavelength(English))
```
""" wienwavelength

@doc """
$(unitext(:wienfrequency,"(𝟑+W₀(-𝟑*exp(-𝟑)))*boltzmann(U)/planck(U)"))

Wien frequency radiation law constant based on Lambert `W₀` evaluation (Hz⋅K⁻¹).
```Julia
julia> wienfrequency(Metric) # Hz⋅K⁻¹
$(wienfrequency(Metric))

julia> wienfrequency(SI2019) # Hz⋅K⁻¹
$(wienfrequency(SI2019))

julia> wienfrequency(Conventional) # Hz⋅K⁻¹
$(wienfrequency(Conventional))

julia> wienfrequency(CODATA) # Hz⋅K⁻¹
$(wienfrequency(CODATA))

julia> wienfrequency(English) # Hz⋅°R⁻¹
$(wienfrequency(English))
```
""" wienfrequency

@doc """
$(unitext(:loschmidt,"atmosphere(U)/boltzmann(U)/temperature(T₀,SI2019,U)"))

Number of molecules (number density) of an ideal gas in a unit volume (m⁻³ or ft⁻³).
```Julia
julia> loschmidt(SI2019) # m⁻³
$(loschmidt(SI2019))

julia> loschmidt(Metric,atm,T₀) # m⁻³
$(loschmidt(Metric,atm,T₀))

julia> loschmidt(Conventional,atm,T₀) # m⁻³
$(loschmidt(Conventional,atm,T₀))

julia> loschmidt(CODATA,atm,T₀) # m⁻³
$(loschmidt(CODATA,atm,T₀))

julia> loschmidt(SI1976,atm,T₀) # m⁻³
$(loschmidt(SI1976,atm,T₀))

julia> loschmidt(English) # ft⁻³
$(loschmidt(English))

julia> loschmidt(IAU) # au⁻³
$(loschmidt(IAU))
```
""" loschmidt

@doc """
$(unitext(:amagat,"loschmidt(U)/avogadro(U)"))

Number of moles of an ideal gas in a unit volume (mol⋅m⁻³ or lb-mol⋅ft⁻³).
```Julia
julia> amagat(Metric) # mol⋅m⁻³
$(amagat(Metric))

julia> amagat(SI2019) # mol⋅m⁻³
$(amagat(SI2019))

julia> amagat(English) # slug-mol⋅ft⁻³
$(amagat(English))
```
""" amagat

@doc """
$(unitext(:hyperfine,"frequency($ΔνCs,U)"))

Unperturbed groundstate hyperfine transition frequency `ΔνCs` of caesium-133 atom (Hz).
```Julia
julia> hyperfine(Metric) # Hz
$(hyperfine(Metric))
```
""" hyperfine, ΔνCs

@doc """
$(unitext(:hubble,"time(U,Hubble)"))

Hubble universe expansion frequency parameter.
```Julia
julia> hubble(Metric)
$(hubble(Metric))

julia> hubble(Hubble)
$(hubble(Hubble))

julia> hubble(Cosmological)
$(hubble(Cosmological))

julia> 𝟏/hubble(Metric)/year(Metric)
$(𝟏/hubble(Metric)/year(Metric))
```
""" hubble

@doc """
$(unitext(:cosmological,"𝟑*darkenergydensity(U)*(hubble(U)/lightspeed(U))^2"))

Cosmological constant from Einstein's controversial theory expanded on by Hubble.
```Julia
julia> cosmological(Metric)
$(cosmological(Metric))

julia> cosmological(Hubble)
$(cosmological(Hubble))

julia> cosmological(Cosmological)
$(cosmological(Cosmological))
```
""" cosmological

@doc """
$(unitext(:solarmass,"mass($(GM☉/G),U)"))

Solar `mass` estimated from gravitational constant estimates (kg or slug).
```Julia
julia> solarmass(Metric) # kg
$(solarmass(Metric))

julia> solarmass(British) # slug
$(solarmass(British))

julia> solarmass(English) # lb
$(solarmass(English))

julia> solarmass(IAUE) # ME
$(solarmass(IAUE))

julia> solarmass(IAUJ) # MJ
$(solarmass(IAUJ))
```
""" solarmass, mₛ

@doc """
$(unitext(:earthmass,"mass($(GME/G),U)"))

Earth `mass` estimated from gravitational constant estimates (kg or slug).
```Julia
julia> earthmass(Metric) # kg
$(earthmass(Metric))

julia> earthmass(British) # slug
$(earthmass(British))

julia> earthmass(English) # lb
$(earthmass(English))

julia> earthmass(IAU) # M☉
$(earthmass(IAU))

julia> earthmass(IAUJ) # MJ
$(earthmass(IAUJ))
```
""" earthmass

@doc """
$(unitext(:jupitermass,"mass($(GMJ/G),U)"))

Jupiter `mass` estimated from gravitational constant estimates (kg or slug).
```Julia
julia> jupitermass(Metric) # kg
$(jupitermass(Metric))

julia> jupitermass(British) # slug
$(jupitermass(British))

julia> jupitermass(English) # lb
$(jupitermass(English))

julia> jupitermass(IAU) # M☉
$(jupitermass(IAU))

julia> jupitermass(IAUE) # ME
$(jupitermass(IAUE))
```
""" jupitermass

@doc """
$(unitext(:lunarmass,"earthmass(U)/μE☾"))

Lunar `mass` estimated from `μE☾` Earth-Moon mass ratio (kg or slug).
```Julia
julia> lunarmass(Metric) # kg
$(lunarmass(Metric))

julia> lunarmass(British) # slug
$(lunarmass(British))

julia> lunarmass(English) # lb
$(lunarmass(English))

julia> lunarmass(IAU) # M☉
$(lunarmass(IAU))

julia> lunarmass(IAUE) # ME
$(lunarmass(IAUE))

julia> lunarmass(IAUJ) # MJ
$(lunarmass(IAUJ))
```
""" lunarmass

@doc """
$(unitext(:gaussgravitation,"sqrt(gravitation(U)*solarmass(U)/astronomicalunit(U)^3)"))

Gaussian  gravitational constant `k` of Newton's laws (Hz or rad⋅D⁻¹).
```Julia
julia> gaussgravitation(Metric)
$(gaussgravitation(Metric))

juila> gaussgravitation(MPH)
$(gaussgravitation(MPH))

julia> gaussgravitation(IAU)
$(gaussgravitation(IAU))
```
""" gaussgravitation, k, kG

@doc """
$(unitext(:gaussianyear,"turn(U)/gaussgravitation(U)"))

Orbit `time` defined by `gaussgravitation` constant `kG` for neglible `mass` satellite.
```Julia
julia> gaussianyear(Metric) # s
$(gaussianyear(Metric))

julia> gaussianyear(MPH) # h
$(gaussianyear(MPH))

julia> gaussianyear(IAU) # D
$(gaussianyear(IAU))
```
""" gaussianyear

@doc """
$(unitext(:siderealyear,"gaussianyear(U)/√(𝟏+earthmass(IAU)+lunarmass(IAU))"))

Orbit `time` defined by `gaussgravitation` constant `kG` and the Earth-Moon system `mass`.
```Julia
julia> siderealyear(Metric) # s
$(siderealyear(Metric))

julia> siderealyear(MPH) # h
$(siderealyear(MPH))

julia> siderealyear(IAU) # D
$(siderealyear(IAU))
```
""" siderealyear

@doc """
$(unitext(:jovianyear,"τ*day(U)*√(jupiterdistance(U)^3/solarmass(U)/gravitation(U))/√(𝟏+jupitermass(IAU))"))

Orbit `time` defined by `jupiterdistance` and the Sun-Jupiter system `mass`.
```Julia
julia> jovianyear(Metric) # s
$(jovianyear(Metric))

julia> jovianyear(MPH) # h
$(jovianyear(MPH))

julia> jovianyear(IAU) # D
$(jovianyear(IAU))
```
""" jovianyear

@doc """
$(unitext(:gaussianmonth,"τ*sqrt(lunardistance(U)^3/earthmass(U)/gravitation(U))"))

Orbit `time` defined by `lunardistance` and `earthmass` for neglible `mass` satellite.
```Julia
julia> gaussianmonth(Metric) # s
$(gaussianmonth(Metric))

julia> gaussianmonth(MPH) # h
$(gaussianmonth(MPH))

julia> gaussianmonth(IAU) # D
$(gaussianmonth(IAU))
```
""" gaussianmonth

@doc """
$(unitext(:siderealmonth,"gaussianmonth(U)/√(𝟏+lunarmass(IAU))"))

Orbit `time` defined by standard `lunardistance` and the Earth-Moon system `mass`.
```Julia
julia> siderealmonth(Metric) # s
$(siderealmonth(Metric))

julia> siderealmonth(MPH) # h
$(siderealmonth(MPH))

julia> siderealmonth(IAU) # D
$(siderealmonth(IAU))
```
""" siderealmonth

@doc """
$(unitext(:synodicmonth,"𝟏/(𝟏/siderealmonth(U)-𝟏/siderealyear(U))"))

Orbit `time` defined by `siderealmonth` and `siderealyear` of the Sun-Earth-Moon system.
```Julia
julia> synodicmonth(Metric) # s
$(synodicmonth(Metric))

julia> synodicmonth(MPH) # h
$(synodicmonth(MPH))

julia> synodicmonth(IAU) # D
$(synodicmonth(IAU))
```
""" synodicmonth

@doc """
$(unitext(:earthradius,"sqrt(earthmass(U)*gravitation(U)/gforce(U))"))

Approximate `length` of standard Earth two-body radius consistent with units (m or ft).
```Julia
julia> earthradius(KKH) # km
$(earthradius(KKH))

julia> earthradius(Nautical) # nm
$(earthradius(Nautical))

julia> earthradius(IAU) # au
$(earthradius(IAU))
```
""" earthradius

@doc """
$(unitext(:greatcircle,"τ*earthradius(U)"))

Approximate `length` of standard Earth two-body circle consistent with units (m or ft).
```Julia
julia> greatcircle(KKH) # km
$(greatcircle(KKH))

julia> greatcircle(Nautical) # nm
$(greatcircle(Nautical))

julia> greatcircle(IAU) # au
$(greatcircle(IAU))
```
""" greatcircle

@doc """
    sackurtetrode(U::UnitSystem,P=atm,T=𝟏,m=mᵤ) = log(kB*T/P*sqrt(m*kB*T/τ/ħ^2)^3)+5/2
    dimensionless : [𝟙], [𝟙], [𝟙], [𝟙], [𝟙]

Ideal gas entropy density for pressure `P`, temperature `T`, atomicmass `m` (dimensionless).
```Julia
julia> sackurtetrode(Metric)
$(sackurtetrode(Metric))

julia> sackurtetrode(SI2019)
$(sackurtetrode(SI2019))

julia> sackurtetrode(Conventional)
$(sackurtetrode(Conventional))

julia> sackurtetrode(CODATA)
$(sackurtetrode(CODATA))

julia> sackurtetrode(SI2019,𝟏𝟎^5)
$(sackurtetrode(SI2019,𝟏𝟎^5))
```
""" sackurtetrode

@doc """
    mechanicalheat(U::UnitSystem) = molargas(U)/molargas(Metric)*calorie(Metric)
    energy : [FL], [FL], [ML²T⁻²], [ML²T⁻²], [ML²T⁻²]

Heat to raise 1 `mass` unit of water by 1 `temperature` unit, or $(normal(molargas(Metric)/calorie(Metric))) `mechanicalheat` per `molaramount` per `temperature` units (J or ft⋅lb).
```Julia
julia> mechanicalheat(Metric) # J
$(mechanicalheat(Metric))

julia> mechanicalheat(English) # ft⋅lb
$(mechanicalheat(English))

julia> mechanicalheat(British) # ft⋅lb
$(mechanicalheat(British))
```
""" mechanicalheat

include("derivedocs.jl")


text:
module models

import ergm.stats

export ExponentialFamily, get_params, update_params, log_likelihood, set_state, update_state, test_state

mutable struct ExponentialFamily
    stats
    params

    function ExponentialFamily(stats, params)
        new(stats, params)
    end
end

function get_params(model :: ExponentialFamily)
    model.params
end

function update_params(model :: ExponentialFamily, params)
    model.params = params
end

function log_likelihood(model :: ExponentialFamily)
    sum(model.params .* stats.get_stats(model.stats))
end

function set_state(model :: ExponentialFamily, state)
    stats.set_state(model.stats, state)
end

function update_state(model :: ExponentialFamily, update)
    stats.update_state(model.stats, update)
end

function test_state(model :: ExponentialFamily, update)
    test_stats = stats.test_state(model.stats, update)
    sum(model.params .* test_stats)
end

function Base.copy(model :: ExponentialFamily)
    ExponentialFamily(copy(model.stats), copy(model.params))
end

end


text:
function zad37(a,b)

    return ceil((sin(pi/2 + a)/a)* sqrt(abs(a+1))* exp(b)),
            round(factorial(a)/(cos(pi-a)+a ^ (b+1)))
end


text:
#=
Power-law slope plot for the paper.
=#
using DrWatson
@quickactivate "VelocityMTD"
using MusicManipulations, PyPlot, Statistics
include(srcdir("style.jl"))
include(srcdir("drums_processing.jl"))
include(srcdir("vel_normalization.jl"))
include(srcdir("PSD.jl"))

# %% ############# Tapping #######################
files = readdir(datadir("exp_raw", "tapping"))
# constants
middles = fill(NaN, length(files))
βs = zeros(length(files))
δs = fill(NaN, length(files))

close("all")

for (i, file) in enumerate(files)
    midi = readMIDIFile(datadir("exp_raw", "tapping", file))
    notes = getnotes(midi)
    @show length(notes)
    middle = middles[i]
    minsub = 1

    figure(figsize = (10,10))
    ax = gca()
    slope1, slope2 = psdfit!(ax, notes, nothing; minsub = minsub, tpq = midi.tpq)
    βs[i] = slope1
    ax.set_title("$(file[1:end-4])_sub=$(minsub)")
    tight_layout()
    savefig(plotsdir("psd", "tapping", file*".png"))
end

bs = middles

tapping = @dict βs δs bs

# %% ############# PLAYALONGS #######################
allnotes = aligned_drums_notes()

# %%
files = readdir(datadir("exp_raw", "playalongs"))
# middles: must be optimized
middles = fill(2.0, length(allnotes))
middles[2] = 4
middles[3] = 4
middles[9] = 3

βs = zeros(length(files))
δs = zeros(length(files))

close("all")

for (i, notes) in enumerate(allnotes)
    notes = combine(notes)
    @show length(notes)

    middle = middles[i]
    file = files[i][1:end-4]
    minsub = 12
    figure(figsize = (10,10))
    ax = gca()
    slope1, slope2 = psdfit!(ax, notes, middle; minsub = minsub)
    # β is large time scales!
    βs[i] = slope2
    δs[i] = slope1
    ax.set_title("playalong_$(file)_middle=$(middle)")
    tight_layout()
    savefig(plotsdir("psd", "playalongs", file*".png"))
end

bs = middles

playalongs = @dict βs δs bs

# %% ############## UWE #######################
files = readdir(datadir("exp_raw", "uwe"))
# constants
middles = fill(2.0, length(files))
middles[2] = 3
middles[3] = 1.5
middles[4] = 3

βs = zeros(length(files))
δs = zeros(length(files))

close("all")

for (i, file) in enumerate(files)

    midi = readMIDIFile(joinpath(datadir("exp_raw", "uwe"), file))
    piano = getnotes(midi, 2)
    d = estimate_delay(piano, 0:1//3:1)
    notes = translate(piano, -round(Int, d))
    @show length(notes)

    middle = middles[i]
    minsub = 6

    figure(figsize = (10,10))
    ax = gca()
    # notice the different window factor and overlap used for the Uwe dataset
    # (because it has smaller amount of points)
    slope1, slope2 = psdfit!(ax, notes, middle; minsub = minsub,
    window_factor = 2, overlap_factor = 4)
    βs[i] = slope2
    δs[i] = slope1
    ax.set_title("$(file[1:end-4])_sub=$(minsub)")
    tight_layout()
    savefig(plotsdir("psd", "uwe", file*".png"))
end

bs = middles

uwe = @dict βs δs bs

# %% ############## PG Music #######################
files = readdir(datadir("exp_raw", "pgmusic"))
# constants
middles = fill(2.0, length(files))
middles[1] = 1.5
middles[15] = 3
middles[17] = 4
βs = zeros(length(files))
δs = zeros(length(files))

close("all")

for (i, file) in enumerate(files)
    midi = readMIDIFile(datadir("exp_raw", "pgmusic", file))
    notes = getnotes(midi)
    @show length(notes)
    middle = middles[i]
    minsub = 12

    figure(figsize = (10,10))
    ax = gca()
    slope1, slope2 = psdfit!(ax, notes, middle; minsub = minsub, tpq = midi.tpq)
    βs[i] = slope2
    δs[i] = slope1
    ax.set_title("$(file[1:end-4])_sub=$(minsub)")
    tight_layout()
    savefig(plotsdir("psd", "pgmusic", file*".png"))
end

bs = middles

pgmusic = @dict βs δs bs

# %% Save data:
@tagsave(datadir("psd", "alldatasets.bson"), @strdict uwe pgmusic playalongs tapping)

max_stars_repo_name: Datseris/VelocityMTD

max_issues_repo_name: Datseris/VelocityMTD

max_forks_repo_name: Datseris/VelocityMTD


text:
input = readlines("d13.in")

function validids(ids)
    return map(x -> parse(Int, x), filter(x -> tryparse(Int, x) != nothing, ids))
end

function validpositions(ids)
    return findall(!=(nothing), tryparse.(Int, ids)) .- 1
end

function findearliest(timestamp, ids)
    delay, i = findmin(mod.(-timestamp,ids))
    return delay, ids[i]
end

function crt(n, a)
    Π = prod(n)
    return mod(sum(ai * invmod(Π ÷ ni, ni) * Π ÷ ni for (ni, ai) in zip(n, -a)), Π)
end

timestamp = parse(Int, input[1])
ids = split(input[2],",")

println("P1: ", prod(findearliest(timestamp, validids(ids))))
println("P2: ", crt(validids(ids), validpositions(ids)))

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode