text:
let

    # Complex and Real Matrix mults, BlockDiagonal, UDT
    for type in (Float64, ComplexF64)
        M = rand(type, 8, 8)
        E = exp(M)
        T = MonteCarlo.taylor_exp(M)
        @test all(abs.(T .- E) .< 10_000eps.(abs.(E)))

        @testset "avx multiplications ($type)" begin
            A = rand(type, 8, 8)
            B = rand(type, 8, 8)
            C = rand(type, 8, 8)
            atol = 100eps(Float64)
            type == ComplexF64 && (atol *= sqrt(2))

            vmul!(C, A, B)
            @test A * B ≈ C atol = atol

            vmul!(C, A, B')
            @test A * B' ≈ C atol = atol

            vmul!(C, A', B)
            @test A' * B ≈ C atol = atol

            vmul!(C, A', B')
            @test A' * B' ≈ C atol = atol

            D = Diagonal(rand(8))
            vmul!(C, A, D)
            @test A * D ≈ C atol = atol

            copyto!(C, A)
            rvmul!(C, D)
            @test A * D ≈ C atol = atol

            copyto!(C, A)
            lvmul!(D, C)
            @test D * A ≈ C atol = atol

            copyto!(C, A)
            rvadd!(C, D)
            @test A + D ≈ C atol = atol

            copyto!(C, A)
            rvadd!(C, B)
            @test A + B ≈ C atol = atol

            vsub!(C, A, I)
            @test A - I ≈ C atol = atol

            if type == Float64
                v = rand(8) .+ 0.5
                w = copy(v)

                vmin!(v, w)
                @test v ≈ min.(1.0, w) atol = atol

                vmininv!(v, w)
                @test v ≈ 1.0 ./ min.(1.0, w) atol = atol

                vmax!(v, w)
                @test v ≈ max.(1.0, w) atol = atol

                vmaxinv!(v, w)
                @test v ≈ 1.0 ./ max.(1.0, w) atol = atol

                v = copy(w)
                vinv!(w)
                @test w ≈ 1.0 ./ v atol = atol
            end
        end

        @testset "UDT transformations + rdivp! ($type)" begin
            U = Matrix{Float64}(undef, 8, 8)
            D = Vector{Float64}(undef, 8)
            T = rand(8, 8)
            X = copy(T)
            MonteCarlo.udt_AVX!(U, D, T)
            @test U * Diagonal(D) * T ≈ X

            U = Matrix{type}(undef, 8, 8)
            T = rand(type, 8, 8)
            X = copy(T)
            pivot = Vector{Int64}(undef, 8)
            tempv = Vector{type}(undef, 8)
            udt_AVX_pivot!(U, D, T)
            @test U * Diagonal(D) * T ≈ X

            copyto!(T, X)
            pivot = Vector{Int64}(undef, 8)
            tempv = Vector{type}(undef, 8)
            udt_AVX_pivot!(U, D, T, pivot, tempv, Val(false))
            # pivoting matrix
            P = zeros(length(pivot), length(pivot))
            for (i, j) in enumerate(pivot)
                P[i, j] = 1.0
            end
            @test U * Diagonal(D) * UpperTriangular(T) * P ≈ X

            u = copy(U)
            t = copy(T)
            tmp = similar(T)
            rdivp!(u, t, tmp, pivot)
            @test u ≈ U * P' / UpperTriangular(T)
        end

        @testset "BlockDiagonal ($type)" begin
            # setindex!
            B = BlockDiagonal(zeros(4, 4), zeros(4, 4))
            for i in 1:8, j in 1:8
                if div(i, 5) == div(j, 5) # cause 5 / 5 = 1 is where we need to jump
                    B[i, j] = 8i+j
                    @test B[i, j] == 8i+j
                else
                    @test_throws BoundsError B[i, j] = 8i+j
                end
            end

            # typing
            b1 = rand(type, 4, 4)
            b2 = rand(type, 4, 4)
            atol = 100eps(Float64)

            B = BlockDiagonal(b1, b2)
            @test B isa BlockDiagonal{type, 2, Matrix{type}}

            # getindex
            for i in 1:4, j in 1:4
                @test B[i, j] == b1[i, j]
                @test B[4+i, 4+j] == b2[i, j]
                @test B[4+i, j] == 0.0
                @test B[i, 4+j] == 0.0
            end

            # Matrix()
            B1 = copy(B)
            M1 = Matrix(B1)
            @test M1 == B1

            B2 = BlockDiagonal(rand(type, 4, 4), rand(type, 4, 4))
            M2 = Matrix(B2)
            B3 = BlockDiagonal(rand(type, 4, 4), rand(type, 4, 4))
            M3 = Matrix(B3)

            # taylor exp
            E = exp(Matrix(B))
            T = MonteCarlo.taylor_exp(B)
            @test all(abs.(T .- E) .< 10_000eps.(abs.(E)))

            # Test (avx) multiplications
            vmul!(B1, B2, B3)
            vmul!(M1, M2, M3)
            @test M1 ≈ B1 atol = atol

            vmul!(B1, B2, adjoint(B3))
            vmul!(M1, M2, adjoint(M3))
            @test M1 ≈ B1

            vmul!(B1, adjoint(B2), B3)
            vmul!(M1, adjoint(M2), M3)
            @test M1 ≈ B1 atol = atol

            D = Diagonal(rand(8))
            vmul!(B1, B2, D)
            vmul!(M1, M2, D)
            @test M1 ≈ B1 atol = atol

            rvmul!(B1, D)
            rvmul!(M1, D)
            @test M1 ≈ B1 atol = atol

            lvmul!(D, B1)
            lvmul!(D, M1)
            @test M1 ≈ B1 atol = atol

            rvadd!(B1, D)
            rvadd!(M1, D)
            @test M1 ≈ B1 atol = atol

            rvadd!(B1, B2)
            rvadd!(M1, M2)
            @test M1 ≈ B1 atol = atol

            vsub!(B1, B2, I)
            @test M2 - I ≈ B1 atol = atol

            x = rand()
            @test B2 * x ≈ M2 * x atol = atol

            # @test log(B2) ≈ log(M2) atol = atol

            # @test det(B2) ≈ det(M2) atol = atol

            transpose!(B1, B2)
            @test B1 ≈ transpose(M2) atol = atol

            vmul!(B1, B2', B3')
            @test B1 ≈ M2' * M3' atol = atol

            # Test UDT and rdivp!
            M2 = Matrix(B2)
            D = rand(8)
            pivot = Vector{Int64}(undef, 8)
            tempv = Vector{type}(undef, 8)
            udt_AVX_pivot!(B1, D, B2, pivot, tempv, Val(false))

            P = zeros(length(pivot), length(pivot))
            # pivot is per block, temporary make it global (in the sense of matrix indices)
            n = size(B1.blocks[1], 1); N = length(B.blocks)
            for i in 1:N; pivot[(i-1)*n+1 : i*n] .+= (i-1)*n; end
            for (i, j) in enumerate(pivot)
                P[i, j] = 1.0
            end
            for i in 1:N; pivot[(i-1)*n+1 : i*n] .-= (i-1)*n; end

            @test Matrix(B1) * Diagonal(D) * UpperTriangular(Matrix(B2)) * P ≈ M2

            M1 = Matrix(B1)
            M2 = Matrix(B2)
            rdivp!(B1, B2, B3, pivot)
            @test B1 ≈ M1 * P' / UpperTriangular(M2)
        end
    end

    @testset "Complex StructArray" begin
        M1 = rand(ComplexF64, 8, 8)
        C1 = StructArray(M1)
        M2 = rand(ComplexF64, 8, 8)
        C2 = StructArray(M2)
        M3 = rand(ComplexF64, 8, 8)
        C3 = StructArray(M3)

        # taylor exp
        E = exp(Matrix(M1))
        T = MonteCarlo.taylor_exp(C1)
        @test all(abs.(T .- E) .< 10_000eps.(abs.(E)))

        @test C1 isa CMat64
        @test M1 == C1

        # Test avx multiplications (adjoint)
        vmul!(C1, C2, C3)
        vmul!(M1, M2, M3)
        @test M1 ≈ C1

        vmul!(C1, C2, adjoint(C3))
        vmul!(M1, M2, adjoint(M3))
        @test M1 ≈ C1

        vmul!(C1, adjoint(C2), C3)
        vmul!(M1, adjoint(M2), M3)
        @test M1 ≈ C1

        vmul!(C1, adjoint(C2), adjoint(C3))
        vmul!(M1, adjoint(M2), adjoint(M3))
        @test M1 ≈ C1

        # Diagonal
        D = Diagonal(rand(8))
        vmul!(C1, C2, D)
        vmul!(M1, M2, D)
        @test M1 == C1

        DC = Diagonal(rand(ComplexF64, 8))
        DCSA = Diagonal(StructArray(DC.diag))
        vmul!(C1, C2, DCSA)
        vmul!(M1, M2, DC)
        @test M1 ≈ C1

        rvmul!(C1, D)
        rvmul!(M1, D)
        @test M1 ≈ C1

        @test_throws ErrorException rvmul!(C1, DCSA)

        lvmul!(D, C1)
        lvmul!(D, M1)
        @test M1 ≈ C1

        @test_throws ErrorException lvmul!(DCSA, C1)

        copyto!(M1, M2)
        copyto!(C1, C2)
        rvadd!(C1, D)
        rvadd!(M1, D)
        @test C1 ≈ M1

        copyto!(M1, M2)
        copyto!(C1, C2)
        rvadd!(C1, C3)
        rvadd!(M1, M3)
        @test C1 ≈ M1

        vsub!(C1, C2, I)
        vsub!(M1, M2, I)
        @test C1 ≈ M1

        # Test UDT backend

        # indmax
        i1, sqn1 = MonteCarlo.indmaxcolumn(C2)
        i2 = LinearAlgebra.indmaxcolumn(M2)
        sqn2 = norm(M2[:, i2])^2
        @test i1 == i2
        @test sqn1 ≈ sqn2

        # reflector
        x1 = MonteCarlo.reflector!(C2, norm(C2[:, 1])^2, 1, size(C2, 1), C1)
        x2 = LinearAlgebra.reflector!(view(M2, :, 1))
        @test x1 ≈ x2
        @test C2 ≈ M2

        # reflectorApply
        C2 = StructArray(M2)
        @views LinearAlgebra.reflectorApply!(M2[:, 1], norm(M2[:, 1]), M2[:, 2:end])
        MonteCarlo.reflectorApply!(norm(C2[:, 1]), C2, 1, size(C2, 1))
        @test M2 ≈ C2

        # Test UDT's, rdivp!

        M2 = rand(ComplexF64, 8, 8)
        C2 = StructArray(M2)
        D = rand(Float64, 8)
        pivot = Vector{Int64}(undef, 8)
        tempv = Vector{ComplexF64}(undef, 8)
        udt_AVX_pivot!(C1, D, C2, pivot, tempv, Val(false))
        P = zeros(length(pivot), length(pivot))
        for (i, j) in enumerate(pivot)
            P[i, j] = 1.0
        end
        @test Matrix(C1) * Diagonal(D) * UpperTriangular(Matrix(C2)) * P ≈ M2

        # rdivp!
        M1 = Matrix(C1)
        M2 = Matrix(C2)
        rdivp!(C1, C2, C3, pivot)
        @test C1 ≈ M1 * P' / UpperTriangular(M2)

        C2 = StructArray(M2)
        udt_AVX_pivot!(C1, D, C2, pivot, tempv, Val(true))
        @test Matrix(C1) * Diagonal(D) * Matrix(C2) ≈ M2

    end
end


text:
# ---
# title: 589. N-ary Tree Preorder Traversal
# id: problem589
# author: zhwang
# date: 2022-03-10
# difficulty: Easy
# categories: Tree
# link: <https://leetcode.com/problems/n-ary-tree-preorder-traversal/description/>
# hidden: true
# ---
#
# Given an n-ary tree, return the _preorder_ traversal of its nodes' values.
#
# _Nary-Tree input serialization  is represented in their level order traversal,
# each group of children is separated by the null value (See examples)._
#
#
#
# **Follow up:**
#
# Recursive solution is trivial, could you do it iteratively?
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2018/10/12/narytreeexample.png)
#
#
#
#     Input: root = [1,null,3,2,4,null,5,6]
#     Output: [1,3,5,6,2,4]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2019/11/08/sample_4_964.png)
#
#
#
#     Input: root = [1,null,2,3,4,5,null,null,6,7,null,8,null,9,10,null,null,11,null,12,null,13,null,null,14]
#     Output: [1,2,3,6,7,11,14,4,8,12,5,9,13,10]
#
#
#
#
# **Constraints:**
#
#   * The height of the n-ary tree is less than or equal to `1000`
#   * The total number of nodes is between `[0, 10^4]`
#
#
## @lc code=start
using LeetCode
preorder(::Nothing) = Int[]
function preorder(root::TreeNode{Int})::Vector{Int}
    stack, ans = [root], Int[]
    while !isempty(stack)
        node = pop!(stack)
        push!(ans, node.val)
        append!(stack, reverse(node.children))
    end
    ans
end

## @lc code=end

## @lc test=start
@testset "589.n-ary-tree-preorder-traversal.jl" begin
    @test
end
## @lc test=end

## @lc info=start
# link: [solution to 589]()
## @lc info=end

ext: jl


text:
"""
Li, D., Wang, Z., Cao, C. and Liu, Y., 2018.
Information entropy based sample reduction for support vector data description.
Applied Soft Computing, 71, pp.1153-1160.

Information Entropy based Sample Reduction for Support Vector Data Description (IESRSVDD)Parameter guideline from the paper, given a data set of size N:
eps (on benchmark data eps ∈ [0.5, 0.6, 0.7])

Simplification of equations:

Dis_ij = K(x_i * x_i) + K(x_j * x_j) - 2 K(x_i * x_j) (Eq. 10)
       = 2 (1 - K(x_i * x_j)) (with RBF Kernel identity == 1)

p_ij = dis_ij / sum(dis_ik for k in 1:n) (Eq. 11)
     = 2 (1 - K(x_i * x_j)) / sum(2 (1 - K(x_i * x_k)) for k in 1:n)
     = (1 - K(x_i * x_j)) / sum(1 - K(x_i * x_k) for k in 1:n)

p_ii = 1

H_i = - sum(p_ik log_2(p_ik) for k in 1:n) (Eq. 12)
    = - sum(p_ik log_2(p_ik) for k in 1:n if k != i) (p_ii == 1 and log_2(1) = 0)
"""
mutable struct IESRSVDD <: Sampler
    kde_cache::Union{KDECache, Nothing}
    gamma::Union{Symbol, Real}
    eps::Real
    function IESRSVDD(kde_cache::Union{KDECache, Nothing}, gamma::Union{Symbol, Real}, eps::Real=0.5)
        (eps >= 0 && eps <= 1) || throw(ArgumentError("Invalid threshold $(eps)."))
        return new(kde_cache, gamma, eps)
    end
end

function IESRSVDD(kde_cache::KDECache, eps::Real=0.5)
    return IESRSVDD(kde_cache, :precomputed, eps)
end

function IESRSVDD(gamma::Union{Symbol, Real}, eps::Real=0.5)
    return IESRSVDD(nothing, gamma, eps)
end

function sample(sampler::IESRSVDD, data::DataSet, labels::Vector{Symbol})::BitArray{1}
    if sampler.kde_cache === nothing || sampler.kde_cache.data != data
        sampler.kde_cache = KDECache(data, sampler.gamma)
    end
    dis_ij = 1 .- sampler.kde_cache.K
    p_ij = dis_ij ./ sum(dis_ij, dims=1)
    p_ij[diagind(p_ij)] .= 1
    h_i = [-sum(p_ij[i, k] * log(2, p_ij[i, k]) for k in 1:size(p_ij, 2)) for i in 1:size(p_ij, 2)]
    h_i_min, h_i_max = extrema(h_i)
    θ = h_i_max - sampler.eps * (h_i_max - h_i_min)
    sample_mask = h_i .> θ
    return sample_mask
end


text:
function add_dual_cone_constraint(dual_model::MOI.ModelLike, primal_model::MOI.ModelLike, vi::Vector{VI},
                                  ci::CI{F, S}) where {F <: MOI.AbstractScalarFunction, S <: MOI.AbstractScalarSet}
    # In this case vi should have only one entry
    return MOI.add_constraint(dual_model, SVF(vi[1]), MOI.dual_set(get_set(primal_model, ci)))
end

function add_dual_cone_constraint(dual_model::MOI.ModelLike, primal_model::MOI.ModelLike, vi::Vector{VI},
                                  ci::CI{F, MOI.EqualTo{T}}) where {T, F <: MOI.AbstractScalarFunction}
    return
end

function add_dual_cone_constraint(dual_model::MOI.ModelLike, primal_model::MOI.ModelLike, vis::Vector{VI},
                                  ci::CI{F, S}) where {F <: MOI.AbstractVectorFunction, S <: MOI.AbstractVectorSet}
    return MOI.add_constraint(dual_model, VVF(vis), MOI.dual_set(get_set(primal_model, ci)))
end

function add_dual_cone_constraint(dual_model::MOI.ModelLike, primal_model::MOI.ModelLike, vis::Vector{VI},
                                  ci::CI{F, MOI.Zeros}) where {F <: MOI.AbstractVectorFunction}
    return
end


text:
function conform_tree_and_data(data::DataFrame, newick::AbstractString)
    taxa_data = data[!, TAXON_NAME]
    tree = readTopology(newick)
    taxa_tree = tipLabels(tree)

    # remove taxa that are not present in both tree and data from the tree
    not_in_data = setdiff(taxa_tree, taxa_data)
    for taxon in not_in_data
        @warn "removing $taxon from tree"
        deleteleaf!(tree, taxon)
    end

    in_both = intersect(taxa_data, taxa_tree)
    in_both_inds = [findfirst(isequal(x), taxa_data) for x in in_both]
    trimmed_df = data[in_both_inds, :]

    return (data = trimmed_df, newick = writeTopology(tree))
end

function merge_traitdata(data::DataFrame...; cut_taxa::Bool = false)
    n = length(data)
    taxon_names = [check_valid(d) for d in data]

    taxon_name = taxon_names[1]
    for i = 2:n
        rename!(data[i], taxon_names[i] => taxon_name)
    end

    args = (on = taxon_name, makeunique = true)

    big_df = cut_taxa ? innerjoin(data...; args...) : outerjoin(data...; args...)
    return big_df
end

function conform_traitdata(data::DataFrame...;kw_args...)
    n = length(data)
    dims = size.(data, 2)

    joint_df = merge_traitdata(data...; kw_args...)
    joint_names = names(joint_df)
    taxon_name = check_valid(joint_df)
    taxa = joint_df[!, taxon_name]

    new_data = Vector{DataFrame}(undef, n)
    ind = 1
    for i = 1:n
        df = DataFrame()
        df[!, taxon_name] = taxa
        for j = 1:(dims[i] - 1) # each of the dims includes a single taxon column
            ind += 1
            df[!, joint_names[ind]] = joint_df[!, ind]
        end

        new_data[i] = df
    end

    return new_data
end


text:
using CBLS
using Test

@testset "CBLS.jl" begin
    include("MOI_wrapper.jl")
    include("JuMP.jl")
end


text:
using LinearAlgebra
using Symbolics

# Problem One

A = Matrix{Int8}([0 5 0; 8 3 7; 9 -2 9])
B = Matrix{Int8}([4 6 -2; 7 2 3; 1 3 -4])
C = Vector{Int8}([-1, 2, 5])

println("Result of A*B: ")
display(A*B)
println("")

println("Result of A*C: ")
display(A*C)
println("")

println("Identity Matrix of Matrix A: ")
identityA = Matrix{Int8}(I,3,3)
display(identityA)
println("")

println("Result of I*A: ")
display(identityA*A)
println("")

println("Result of A*I: ")
display(A*identityA)
println("")

println("Matrix A")
display(A)
println("")

#Problem Two

A = Matrix{Int8}([5  7  2  0  3  5;
     3  8 -3 -5  0  8;
     1  4  0  7 15  9;
     0 10  5 12  3 -1;
     2 -5  9  2 18 -10])

B = Matrix{Int8}([ 2 10  0;
      8  7  5;
     -5  2 -4;
      4  8 13;
      3 12  0;
      1  5  7])

println("Result C = A*B: ")
display(A*B)
println("")

println("Transpose of Matrix A: ")
display(transpose(A))
println("")

# Problem Three

A = Matrix{Int8}([1 5 0; 8 3 7; 6 -2 9])

println("Determinant of Matrix A: ")
display(det(A))
println("")

println("Determinant of the Transpose of Matrix A: ")
display(det(transpose(A)))
println("")

println("Inverse of Matrix A: ")
display(inv(A))
println("")

# Problem Four
@variables k1
A = [0 0 0 0;
     0 k1 0 -k1;
     0 0 0 0;
     0 -k1 0 k1;]
println("Local Matrix A: ")
display(A)
println("")

G = [0 0 0 0 0 0 0 0 0 0;
     0 0 0 0 0 0 0 0 0 0;
     0 0 0 0 0 0 0 0 0 0;
     0 0 0 k1 0 0 0 -k1 0 0;
     0 0 0 0 0 0 0 0 0 0;
     0 0 0 0 0 0 0 0 0 0;
     0 0 0 0 0 0 0 0 0 0;
     0 0 0 -k1 0 0 0 k1 0 0;
     0 0 0 0 0 0 0 0 0 0;
     0 0 0 0 0 0 0 0 0 0;]
println("Mapped global matrix, G: ")
display(G)


text:
module MLDataPattern
using StatsBase
using LearnBase
using MLLabelUtils

using LearnBase: ObsDimension
import StatsBase: nobs
import LearnBase: getobs, getobs!, gettarget, gettargets, targets, datasubset, default_obsdim, DataView,
    AbstractObsView, AbstractBatchView, DataIterator, AbstractDataIterator, ObsIterator, BatchIterator

using Base.Cartesian
using Random
using SparseArrays

export

    ObsDim,

    nobs,
    getobs,
    getobs!,

    randobs,

    DataSubset,
    datasubset,

    shuffleobs,
    splitobs,

    ObsView,
    BatchView,
    obsview,
    batchview,
    batchsize,

    SlidingWindow,
    slidingwindow,

    targets,
    eachtarget,

    stratifiedobs,

    oversample,
    undersample,
    upsample,
    downsample,

    FoldsView,
    kfolds,
    leaveout,

    RandomObs,
    BalancedObs,
    RandomBatches,
    BufferGetObs,
    eachobs,
    eachbatch

obsdim_string(::ObsDim.First) = ":first"
obsdim_string(::ObsDim.Last) = ":last"
obsdim_string(::ObsDim.Constant{D}) where {D} = string(D)
obsdim_string(::ObsDim.Undefined) = "\"NA\""
obsdim_string(obsdim::Tuple) = string("(", join(map(obsdim_string, obsdim), ", "), ")")

include("container.jl")
include("datasubset.jl")
include("randobs.jl")
include("shuffleobs.jl")
include("splitobs.jl")
include("dataview.jl")
include("slidingwindow.jl")
include("targets.jl")
include("stratifiedobs.jl")
include("resample.jl")
include("folds.jl")
include("dataiterator.jl")

end # module


text:
#!/usr/bin/env julia

using Luxor

function ellipse2c(focus1::Point, focus2::Point, k; stepvalue=pi/100)
    # k is sum of distances to focii of any points on the ellipse)
    a = k/2  # Major axis
    cpoint = midpoint(focus1, focus2)
    dc = norm(focus1, cpoint) # Distance of the foci to the center
    b = sqrt(abs(a^2 - dc^2)) # Minor axis
    phi = slope(O, focus2) # Angle between the major axis and the x-axis
    drawing = false
    for t in 0:stepvalue:2pi
        xt = cpoint.x + a * cos(t) * cos(phi) - b * sin(t) * sin(phi)
        yt = cpoint.y + a * cos(t) * sin(phi) + b * sin(t) * cos(phi)
        if drawing
            line(Point(xt, yt))
        else
            move(Point(xt, yt))
            drawing = true
        end
    end
    closepath()
    stroke()
end

"""
Draw a logo with ellipse and colored circles
"""
function draw_logo(fname)
    Drawing(500, 500, fname)
    setline(3)
    origin()
    cpoints = ngon(O, 200, 3, vertices=true)

    ellipse2c(O, cpoints[1], 280)
    ellipse2c(O, cpoints[2], 280)
    ellipse2c(O, cpoints[3], 280)

    darker_purple = (0.584, 0.345, 0.698)
    lighter_purple  = (0.667, 0.475, 0.757)
    darker_green  = (0.22, 0.596, 0.149)
    lighter_green  = (0.376, 0.678, 0.318)
    darker_red  = (0.796, 0.235, 0.2)
    lighter_red  = (0.835, 0.388, 0.361)

    ballradius = 40
    sethue(darker_red...)
    circle(cpoints[1], 0.75ballradius, :fill)
    sethue(lighter_red...)
    circle(cpoints[1], 0.65ballradius, :fill)

    sethue(darker_green...)
    circle(cpoints[2], 0.75ballradius, :fill)
    sethue(lighter_green...)
    circle(cpoints[2], 0.65ballradius, :fill)

    sethue(darker_purple...)
    circle(cpoints[3], 0.75ballradius, :fill)
    sethue(lighter_purple...)
    circle(cpoints[3], 0.65ballradius, :fill)

    sethue("gold")
    star(O, 35, 12, 0.7, pi/3, :fill)
    sethue("yellow")
    star(O, 35, 15, 0.3, 0, :fill)
    finish()
    preview()
end

draw_logo("/tmp/logo.png")


text:
# ---------------------------------------------------------------------------- #
#
#   quadratureTet.jl
#
#   Get Legendre quadrature points for Tetrahedron
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

function quadratureTet( P::PG1 )

  nq = 1
  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)
  zc     = Array{Float64}(nq,1)

  weight[1] = 1
  xc[1]     = 1./4.
  yc[1]     = 1./4.
  zc[1]     = 1./4.

  pts = hcat( xc, yc, zc )

  weight = weight / 6

  return pts, weight

end

function quadratureTet( P::PG2 )

  nq  = 4
  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)
  zc     = Array{Float64}(nq,1)

  w = 1./4.
  a = (5 - sqrt(5))/20.

  weight[1] = w
  xc[1]     = 1 - 3*a
  yc[1]     = a
  zc[1]     = a

  weight[2] = w
  xc[2]     = a
  yc[2]     = 1 - 3*a
  zc[2]     = a

  weight[3] = w
  xc[3]     = a
  yc[3]     = a
  zc[3]     = 1 - 3*a

  weight[4] = 1./4.
  xc[4] = a
  yc[4] = a
  zc[4] = a

  pts = hcat( xc, yc, zc )

  weight = weight / 6

  return pts, weight

end

function quadratureTet( P::PG3 )

  nq  = 8
  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)
  zc     = Array{Float64}(nq,1)

  xc[ 1] =   3.2805469671142667e-01;  yc[ 1] =   3.2805469671142667e-01;  zc[ 1] =   1.5835909865719922e-02
  xc[ 2] =   3.2805469671142667e-01;  yc[ 2] =   1.5835909865719922e-02;  zc[ 2] =   3.2805469671142667e-01
  xc[ 3] =   1.5835909865719922e-02;  yc[ 3] =   3.2805469671142667e-01;  zc[ 3] =   3.2805469671142667e-01
  xc[ 4] =   3.2805469671142667e-01;  yc[ 4] =   3.2805469671142667e-01;  zc[ 4] =   3.2805469671142667e-01
  xc[ 5] =   1.0695227393293068e-01;  yc[ 5] =   1.0695227393293068e-01;  zc[ 5] =   6.7914317820120795e-01
  xc[ 6] =   1.0695227393293068e-01;  yc[ 6] =   6.7914317820120795e-01;  zc[ 6] =   1.0695227393293068e-01
  xc[ 7] =   6.7914317820120795e-01;  yc[ 7] =   1.0695227393293068e-01;  zc[ 7] =   1.0695227393293068e-01
  xc[ 8] =   1.0695227393293068e-01;  yc[ 8] =   1.0695227393293068e-01;  zc[ 8] =   1.0695227393293068e-01

  weight[ 1] =   2.3087994418643690e-02# * 6.0
  weight[ 2] =   2.3087994418643690e-02# * 6.0
  weight[ 3] =   2.3087994418643690e-02# * 6.0
  weight[ 4] =   2.3087994418643690e-02# * 6.0
  weight[ 5] =   1.8578672248022978e-02# * 6.0
  weight[ 6] =   1.8578672248022978e-02# * 6.0
  weight[ 7] =   1.8578672248022978e-02# * 6.0
  weight[ 8] =   1.8578672248022978e-02# * 6.0

  pts = hcat( xc, yc, zc )

  return pts, weight

end

function quadratureTet( P::PG4 )

  nq = 14
  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)
  zc     = Array{Float64}(nq,1)

  xc[ 1] =   9.2735250310891221e-02;  yc[ 1] =   9.2735250310891221e-02;  zc[ 1] =   7.2179424906732637e-01
  xc[ 2] =   9.2735250310891221e-02;  yc[ 2] =   7.2179424906732637e-01;  zc[ 2] =   9.2735250310891221e-02
  xc[ 3] =   7.2179424906732637e-01;  yc[ 3] =   9.2735250310891221e-02;  zc[ 3] =   9.2735250310891221e-02
  xc[ 4] =   9.2735250310891221e-02;  yc[ 4] =   9.2735250310891221e-02;  zc[ 4] =   9.2735250310891221e-02
  xc[ 5] =   3.1088591926330061e-01;  yc[ 5] =   3.1088591926330061e-01;  zc[ 5] =   6.7342242210098213e-02
  xc[ 6] =   3.1088591926330061e-01;  yc[ 6] =   6.7342242210098213e-02;  zc[ 6] =   3.1088591926330061e-01
  xc[ 7] =   6.7342242210098213e-02;  yc[ 7] =   3.1088591926330061e-01;  zc[ 7] =   3.1088591926330061e-01
  xc[ 8] =   3.1088591926330061e-01;  yc[ 8] =   3.1088591926330061e-01;  zc[ 8] =   3.1088591926330061e-01
  xc[ 9] =   4.5503704125649649e-02;  yc[ 9] =   4.5449629587435036e-01;  zc[ 9] =   4.5449629587435036e-01
  xc[10] =   4.5449629587435036e-01;  yc[10] =   4.5503704125649649e-02;  zc[10] =   4.5449629587435036e-01
  xc[11] =   4.5449629587435036e-01;  yc[11] =   4.5449629587435036e-01;  zc[11] =   4.5503704125649649e-02
  xc[12] =   4.5503704125649649e-02;  yc[12] =   4.5449629587435036e-01;  zc[12] =   4.5503704125649649e-02
  xc[13] =   4.5503704125649649e-02;  yc[13] =   4.5503704125649649e-02;  zc[13] =   4.5449629587435036e-01
  xc[14] =   4.5449629587435036e-01;  yc[14] =   4.5503704125649649e-02;  zc[14] =   4.5503704125649649e-02

  weight[ 1] =   1.2248840519393659e-02# * 6.0
  weight[ 2] =   1.2248840519393659e-02# * 6.0
  weight[ 3] =   1.2248840519393659e-02# * 6.0
  weight[ 4] =   1.2248840519393659e-02# * 6.0
  weight[ 5] =   1.8781320953002643e-02# * 6.0
  weight[ 6] =   1.8781320953002643e-02# * 6.0
  weight[ 7] =   1.8781320953002643e-02# * 6.0
  weight[ 8] =   1.8781320953002643e-02# * 6.0
  weight[ 9] =   7.0910034628469112e-03# * 6.0
  weight[10] =   7.0910034628469112e-03# * 6.0
  weight[11] =   7.0910034628469112e-03# * 6.0
  weight[12] =   7.0910034628469112e-03# * 6.0
  weight[13] =   7.0910034628469112e-03# * 6.0
  weight[14] =   7.0910034628469112e-03# * 6.0

  pts = hcat( xc, yc, zc )

  return pts, weight

end

function quadratureTet( P::PG5 )

  nq = 14
  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)
  zc     = Array{Float64}(nq,1)

  xc[ 1] =   3.1088591926330061e-01;  yc[ 1] =   3.1088591926330061e-01;  zc[ 1] =   6.7342242210098213e-02
  xc[ 2] =   3.1088591926330061e-01;  yc[ 2] =   6.7342242210098213e-02;  zc[ 2] =   3.1088591926330061e-01
  xc[ 3] =   6.7342242210098213e-02;  yc[ 3] =   3.1088591926330061e-01;  zc[ 3] =   3.1088591926330061e-01
  xc[ 4] =   3.1088591926330061e-01;  yc[ 4] =   3.1088591926330061e-01;  zc[ 4] =   3.1088591926330061e-01
  xc[ 5] =   9.2735250310891221e-02;  yc[ 5] =   9.2735250310891221e-02;  zc[ 5] =   7.2179424906732637e-01
  xc[ 6] =   9.2735250310891221e-02;  yc[ 6] =   7.2179424906732637e-01;  zc[ 6] =   9.2735250310891221e-02
  xc[ 7] =   7.2179424906732637e-01;  yc[ 7] =   9.2735250310891221e-02;  zc[ 7] =   9.2735250310891221e-02
  xc[ 8] =   9.2735250310891221e-02;  yc[ 8] =   9.2735250310891221e-02;  zc[ 8] =   9.2735250310891221e-02
  xc[ 9] =   4.5503704125649649e-02;  yc[ 9] =   4.5449629587435036e-01;  zc[ 9] =   4.5449629587435036e-01
  xc[10] =   4.5449629587435036e-01;  yc[10] =   4.5503704125649649e-02;  zc[10] =   4.5449629587435036e-01
  xc[11] =   4.5449629587435036e-01;  yc[11] =   4.5449629587435036e-01;  zc[11] =   4.5503704125649649e-02
  xc[12] =   4.5503704125649649e-02;  yc[12] =   4.5449629587435036e-01;  zc[12] =   4.5503704125649649e-02
  xc[13] =   4.5503704125649649e-02;  yc[13] =   4.5503704125649649e-02;  zc[13] =   4.5449629587435036e-01
  xc[14] =   4.5449629587435036e-01;  yc[14] =   4.5503704125649649e-02;  zc[14] =   4.5503704125649649e-02

  weight[ 1] =   1.8781320953002643e-02# * 6.0
  weight[ 2] =   1.8781320953002643e-02# * 6.0
  weight[ 3] =   1.8781320953002643e-02# * 6.0
  weight[ 4] =   1.8781320953002643e-02# * 6.0
  weight[ 5] =   1.2248840519393659e-02# * 6.0
  weight[ 6] =   1.2248840519393659e-02# * 6.0
  weight[ 7] =   1.2248840519393659e-02# * 6.0
  weight[ 8] =   1.2248840519393659e-02# * 6.0
  weight[ 9] =   7.0910034628469112e-03# * 6.0
  weight[10] =   7.0910034628469112e-03# * 6.0
  weight[11] =   7.0910034628469112e-03# * 6.0
  weight[12] =   7.0910034628469112e-03# * 6.0
  weight[13] =   7.0910034628469112e-03# * 6.0
  weight[14] =   7.0910034628469112e-03# * 6.0

  pts = hcat( xc, yc, zc )

  return pts, weight

end

function quadratureTet( P::PG6 )

  nq = 24
  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)
  zc     = Array{Float64}(nq,1)

  xc[ 1] =   2.1460287125915203e-01;  yc[ 1] =   2.1460287125915203e-01;  zc[ 1] =   3.5619138622254387e-01
  xc[ 2] =   2.1460287125915203e-01;  yc[ 2] =   3.5619138622254387e-01;  zc[ 2] =   2.1460287125915203e-01
  xc[ 3] =   3.5619138622254387e-01;  yc[ 3] =   2.1460287125915203e-01;  zc[ 3] =   2.1460287125915203e-01
  xc[ 4] =   2.1460287125915203e-01;  yc[ 4] =   2.1460287125915203e-01;  zc[ 4] =   2.1460287125915203e-01
  xc[ 5] =   4.0673958534611351e-02;  yc[ 5] =   4.0673958534611351e-02;  zc[ 5] =   8.7797812439616596e-01
  xc[ 6] =   4.0673958534611351e-02;  yc[ 6] =   8.7797812439616596e-01;  zc[ 6] =   4.0673958534611351e-02
  xc[ 7] =   8.7797812439616596e-01;  yc[ 7] =   4.0673958534611351e-02;  zc[ 7] =   4.0673958534611351e-02
  xc[ 8] =   4.0673958534611351e-02;  yc[ 8] =   4.0673958534611351e-02;  zc[ 8] =   4.0673958534611351e-02
  xc[ 9] =   3.2233789014227548e-01;  yc[ 9] =   3.2233789014227548e-01;  zc[ 9] =   3.2986329573173601e-02
  xc[10] =   3.2233789014227548e-01;  yc[10] =   3.2986329573173601e-02;  zc[10] =   3.2233789014227548e-01
  xc[11] =   3.2986329573173601e-02;  yc[11] =   3.2233789014227548e-01;  zc[11] =   3.2233789014227548e-01
  xc[12] =   3.2233789014227548e-01;  yc[12] =   3.2233789014227548e-01;  zc[12] =   3.2233789014227548e-01
  xc[13] =   6.3661001875017525e-02;  yc[13] =   6.0300566479164919e-01;  zc[13] =   2.6967233145831571e-01
  xc[14] =   6.3661001875017525e-02;  yc[14] =   2.6967233145831571e-01;  zc[14] =   6.0300566479164919e-01
  xc[15] =   6.0300566479164919e-01;  yc[15] =   6.3661001875017525e-02;  zc[15] =   2.6967233145831571e-01
  xc[16] =   6.0300566479164919e-01;  yc[16] =   2.6967233145831571e-01;  zc[16] =   6.3661001875017525e-02
  xc[17] =   2.6967233145831571e-01;  yc[17] =   6.3661001875017525e-02;  zc[17] =   6.0300566479164919e-01
  xc[18] =   2.6967233145831571e-01;  yc[18] =   6.0300566479164919e-01;  zc[18] =   6.3661001875017525e-02
  xc[19] =   6.3661001875017525e-02;  yc[19] =   6.3661001875017525e-02;  zc[19] =   2.6967233145831571e-01
  xc[20] =   6.3661001875017525e-02;  yc[20] =   2.6967233145831571e-01;  zc[20] =   6.3661001875017525e-02
  xc[21] =   2.6967233145831571e-01;  yc[21] =   6.3661001875017525e-02;  zc[21] =   6.3661001875017525e-02
  xc[22] =   6.3661001875017525e-02;  yc[22] =   6.3661001875017525e-02;  zc[22] =   6.0300566479164919e-01
  xc[23] =   6.3661001875017525e-02;  yc[23] =   6.0300566479164919e-01;  zc[23] =   6.3661001875017525e-02
  xc[24] =   6.0300566479164919e-01;  yc[24] =   6.3661001875017525e-02;  zc[24] =   6.3661001875017525e-02

  weight[ 1] =   6.6537917096945827e-03# * 6.0
  weight[ 2] =   6.6537917096945827e-03# * 6.0
  weight[ 3] =   6.6537917096945827e-03# * 6.0
  weight[ 4] =   6.6537917096945827e-03# * 6.0
  weight[ 5] =   1.6795351758867739e-03# * 6.0
  weight[ 6] =   1.6795351758867739e-03# * 6.0
  weight[ 7] =   1.6795351758867739e-03# * 6.0
  weight[ 8] =   1.6795351758867739e-03# * 6.0
  weight[ 9] =   9.2261969239424545e-03# * 6.0
  weight[10] =   9.2261969239424545e-03# * 6.0
  weight[11] =   9.2261969239424545e-03# * 6.0
  weight[12] =   9.2261969239424545e-03# * 6.0
  weight[13] =   8.0357142857142867e-03# * 6.0
  weight[14] =   8.0357142857142867e-03# * 6.0
  weight[15] =   8.0357142857142867e-03# * 6.0
  weight[16] =   8.0357142857142867e-03# * 6.0
  weight[17] =   8.0357142857142867e-03# * 6.0
  weight[18] =   8.0357142857142867e-03# * 6.0
  weight[19] =   8.0357142857142867e-03# * 6.0
  weight[20] =   8.0357142857142867e-03# * 6.0
  weight[21] =   8.0357142857142867e-03# * 6.0
  weight[22] =   8.0357142857142867e-03# * 6.0
  weight[23] =   8.0357142857142867e-03# * 6.0
  weight[24] =   8.0357142857142867e-03# * 6.0

  pts = hcat( xc, yc, zc )

  return pts, weight

end

function quadratureTet( P::PG7 )

  nq = 35
  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)
  zc     = Array{Float64}(nq,1)

  xc[ 1] =   2.5000000000000000e-01;  yc[ 1] =   2.5000000000000000e-01;  zc[ 1] =   2.5000000000000000e-01
  xc[ 2] =   3.1570114977820279e-01;  yc[ 2] =   3.1570114977820279e-01;  zc[ 2] =   5.2896550665391562e-02
  xc[ 3] =   3.1570114977820279e-01;  yc[ 3] =   5.2896550665391562e-02;  zc[ 3] =   3.1570114977820279e-01
  xc[ 4] =   5.2896550665391562e-02;  yc[ 4] =   3.1570114977820279e-01;  zc[ 4] =   3.1570114977820279e-01
  xc[ 5] =   3.1570114977820279e-01;  yc[ 5] =   3.1570114977820279e-01;  zc[ 5] =   3.1570114977820279e-01
  xc[ 6] =   5.0489822598396371e-02;  yc[ 6] =   4.4951017740160365e-01;  zc[ 6] =   4.4951017740160365e-01
  xc[ 7] =   4.4951017740160365e-01;  yc[ 7] =   5.0489822598396371e-02;  zc[ 7] =   4.4951017740160365e-01
  xc[ 8] =   4.4951017740160365e-01;  yc[ 8] =   4.4951017740160365e-01;  zc[ 8] =   5.0489822598396371e-02
  xc[ 9] =   5.0489822598396371e-02;  yc[ 9] =   4.4951017740160365e-01;  zc[ 9] =   5.0489822598396371e-02
  xc[10] =   5.0489822598396371e-02;  yc[10] =   5.0489822598396371e-02;  zc[10] =   4.4951017740160365e-01
  xc[11] =   4.4951017740160365e-01;  yc[11] =   5.0489822598396371e-02;  zc[11] =   5.0489822598396371e-02
  xc[12] =   1.8883383102600104e-01;  yc[12] =   5.7517163758700007e-01;  zc[12] =   4.7160700360997843e-02
  xc[13] =   1.8883383102600104e-01;  yc[13] =   4.7160700360997843e-02;  zc[13] =   5.7517163758700007e-01
  xc[14] =   5.7517163758700007e-01;  yc[14] =   1.8883383102600104e-01;  zc[14] =   4.7160700360997843e-02
  xc[15] =   5.7517163758700007e-01;  yc[15] =   4.7160700360997843e-02;  zc[15] =   1.8883383102600104e-01
  xc[16] =   4.7160700360997843e-02;  yc[16] =   1.8883383102600104e-01;  zc[16] =   5.7517163758700007e-01
  xc[17] =   4.7160700360997843e-02;  yc[17] =   5.7517163758700007e-01;  zc[17] =   1.8883383102600104e-01
  xc[18] =   1.8883383102600104e-01;  yc[18] =   1.8883383102600104e-01;  zc[18] =   4.7160700360997843e-02
  xc[19] =   1.8883383102600104e-01;  yc[19] =   4.7160700360997843e-02;  zc[19] =   1.8883383102600104e-01
  xc[20] =   4.7160700360997843e-02;  yc[20] =   1.8883383102600104e-01;  zc[20] =   1.8883383102600104e-01
  xc[21] =   1.8883383102600104e-01;  yc[21] =   1.8883383102600104e-01;  zc[21] =   5.7517163758700007e-01
  xc[22] =   1.8883383102600104e-01;  yc[22] =   5.7517163758700007e-01;  zc[22] =   1.8883383102600104e-01
  xc[23] =   5.7517163758700007e-01;  yc[23] =   1.8883383102600104e-01;  zc[23] =   1.8883383102600104e-01
  xc[24] =   2.1265472541483248e-02;  yc[24] =   8.1083024109854851e-01;  zc[24] =   1.4663881381848498e-01
  xc[25] =   2.1265472541483248e-02;  yc[25] =   1.4663881381848498e-01;  zc[25] =   8.1083024109854851e-01
  xc[26] =   8.1083024109854851e-01;  yc[26] =   2.1265472541483248e-02;  zc[26] =   1.4663881381848498e-01
  xc[27] =   8.1083024109854851e-01;  yc[27] =   1.4663881381848498e-01;  zc[27] =   2.1265472541483248e-02
  xc[28] =   1.4663881381848498e-01;  yc[28] =   2.1265472541483248e-02;  zc[28] =   8.1083024109854851e-01
  xc[29] =   1.4663881381848498e-01;  yc[29] =   8.1083024109854851e-01;  zc[29] =   2.1265472541483248e-02
  xc[30] =   2.1265472541483248e-02;  yc[30] =   2.1265472541483248e-02;  zc[30] =   1.4663881381848498e-01
  xc[31] =   2.1265472541483248e-02;  yc[31] =   1.4663881381848498e-01;  zc[31] =   2.1265472541483248e-02
  xc[32] =   1.4663881381848498e-01;  yc[32] =   2.1265472541483248e-02;  zc[32] =   2.1265472541483248e-02
  xc[33] =   2.1265472541483248e-02;  yc[33] =   2.1265472541483248e-02;  zc[33] =   8.1083024109854851e-01
  xc[34] =   2.1265472541483248e-02;  yc[34] =   8.1083024109854851e-01;  zc[34] =   2.1265472541483248e-02
  xc[35] =   8.1083024109854851e-01;  yc[35] =   2.1265472541483248e-02;  zc[35] =   2.1265472541483248e-02

  weight[ 1] =   1.5914214910688475e-02# * 6.0
  weight[ 2] =   7.0549302016611713e-03# * 6.0
  weight[ 3] =   7.0549302016611713e-03# * 6.0
  weight[ 4] =   7.0549302016611713e-03# * 6.0
  weight[ 5] =   7.0549302016611713e-03# * 6.0
  weight[ 6] =   5.3161546388095964e-03# * 6.0
  weight[ 7] =   5.3161546388095964e-03# * 6.0
  weight[ 8] =   5.3161546388095964e-03# * 6.0
  weight[ 9] =   5.3161546388095964e-03# * 6.0
  weight[10] =   5.3161546388095964e-03# * 6.0
  weight[11] =   5.3161546388095964e-03# * 6.0
  weight[12] =   6.2011884547224366e-03# * 6.0
  weight[13] =   6.2011884547224366e-03# * 6.0
  weight[14] =   6.2011884547224366e-03# * 6.0
  weight[15] =   6.2011884547224366e-03# * 6.0
  weight[16] =   6.2011884547224366e-03# * 6.0
  weight[17] =   6.2011884547224366e-03# * 6.0
  weight[18] =   6.2011884547224366e-03# * 6.0
  weight[19] =   6.2011884547224366e-03# * 6.0
  weight[20] =   6.2011884547224366e-03# * 6.0
  weight[21] =   6.2011884547224366e-03# * 6.0
  weight[22] =   6.2011884547224366e-03# * 6.0
  weight[23] =   6.2011884547224366e-03# * 6.0
  weight[24] =   1.3517951383172236e-03# * 6.0
  weight[25] =   1.3517951383172236e-03# * 6.0
  weight[26] =   1.3517951383172236e-03# * 6.0
  weight[27] =   1.3517951383172236e-03# * 6.0
  weight[28] =   1.3517951383172236e-03# * 6.0
  weight[29] =   1.3517951383172236e-03# * 6.0
  weight[30] =   1.3517951383172236e-03# * 6.0
  weight[31] =   1.3517951383172236e-03# * 6.0
  weight[32] =   1.3517951383172236e-03# * 6.0
  weight[33] =   1.3517951383172236e-03# * 6.0
  weight[34] =   1.3517951383172236e-03# * 6.0
  weight[35] =   1.3517951383172236e-03# * 6.0

  pts = hcat( xc, yc, zc )

  return pts, weight

end

function quadratureTet( P::PG9 )

  nq = 46
  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)
  zc     = Array{Float64}(nq,1)

  xc[ 1] =   3.9675423070389904e-02;  yc[ 1] =   3.9675423070389904e-02;  zc[ 1] =   8.8097373078883034e-01
  xc[ 2] =   3.9675423070389904e-02;  yc[ 2] =   8.8097373078883034e-01;  zc[ 2] =   3.9675423070389904e-02
  xc[ 3] =   8.8097373078883034e-01;  yc[ 3] =   3.9675423070389904e-02;  zc[ 3] =   3.9675423070389904e-02
  xc[ 4] =   3.9675423070389904e-02;  yc[ 4] =   3.9675423070389904e-02;  zc[ 4] =   3.9675423070389904e-02
  xc[ 5] =   3.1448780069809629e-01;  yc[ 5] =   3.1448780069809629e-01;  zc[ 5] =   5.6536597905711083e-02
  xc[ 6] =   3.1448780069809629e-01;  yc[ 6] =   5.6536597905711083e-02;  zc[ 6] =   3.1448780069809629e-01
  xc[ 7] =   5.6536597905711083e-02;  yc[ 7] =   3.1448780069809629e-01;  zc[ 7] =   3.1448780069809629e-01
  xc[ 8] =   3.1448780069809629e-01;  yc[ 8] =   3.1448780069809629e-01;  zc[ 8] =   3.1448780069809629e-01
  xc[ 9] =   1.0198669306270330e-01;  yc[ 9] =   1.0198669306270330e-01;  zc[ 9] =   6.9403992081189014e-01
  xc[10] =   1.0198669306270330e-01;  yc[10] =   6.9403992081189014e-01;  zc[10] =   1.0198669306270330e-01
  xc[11] =   6.9403992081189014e-01;  yc[11] =   1.0198669306270330e-01;  zc[11] =   1.0198669306270330e-01
  xc[12] =   1.0198669306270330e-01;  yc[12] =   1.0198669306270330e-01;  zc[12] =   1.0198669306270330e-01
  xc[13] =   1.8420369694919150e-01;  yc[13] =   1.8420369694919150e-01;  zc[13] =   4.4738890915242546e-01
  xc[14] =   1.8420369694919150e-01;  yc[14] =   4.4738890915242546e-01;  zc[14] =   1.8420369694919150e-01
  xc[15] =   4.4738890915242546e-01;  yc[15] =   1.8420369694919150e-01;  zc[15] =   1.8420369694919150e-01
  xc[16] =   1.8420369694919150e-01;  yc[16] =   1.8420369694919150e-01;  zc[16] =   1.8420369694919150e-01
  xc[17] =   6.3436287754539897e-02;  yc[17] =   4.3656371224546009e-01;  zc[17] =   4.3656371224546009e-01
  xc[18] =   4.3656371224546009e-01;  yc[18] =   6.3436287754539897e-02;  zc[18] =   4.3656371224546009e-01
  xc[19] =   4.3656371224546009e-01;  yc[19] =   4.3656371224546009e-01;  zc[19] =   6.3436287754539897e-02
  xc[20] =   6.3436287754539897e-02;  yc[20] =   4.3656371224546009e-01;  zc[20] =   6.3436287754539897e-02
  xc[21] =   6.3436287754539897e-02;  yc[21] =   6.3436287754539897e-02;  zc[21] =   4.3656371224546009e-01
  xc[22] =   4.3656371224546009e-01;  yc[22] =   6.3436287754539897e-02;  zc[22] =   6.3436287754539897e-02
  xc[23] =   2.1690162067728006e-02;  yc[23] =   7.1993192203946599e-01;  zc[23] =   2.3668775382507801e-01
  xc[24] =   2.1690162067728006e-02;  yc[24] =   2.3668775382507801e-01;  zc[24] =   7.1993192203946599e-01
  xc[25] =   7.1993192203946599e-01;  yc[25] =   2.1690162067728006e-02;  zc[25] =   2.3668775382507801e-01
  xc[26] =   7.1993192203946599e-01;  yc[26] =   2.3668775382507801e-01;  zc[26] =   2.1690162067728006e-02
  xc[27] =   2.3668775382507801e-01;  yc[27] =   2.1690162067728006e-02;  zc[27] =   7.1993192203946599e-01
  xc[28] =   2.3668775382507801e-01;  yc[28] =   7.1993192203946599e-01;  zc[28] =   2.1690162067728006e-02
  xc[29] =   2.1690162067728006e-02;  yc[29] =   2.1690162067728006e-02;  zc[29] =   2.3668775382507801e-01
  xc[30] =   2.1690162067728006e-02;  yc[30] =   2.3668775382507801e-01;  zc[30] =   2.1690162067728006e-02
  xc[31] =   2.3668775382507801e-01;  yc[31] =   2.1690162067728006e-02;  zc[31] =   2.1690162067728006e-02
  xc[32] =   2.1690162067728006e-02;  yc[32] =   2.1690162067728006e-02;  zc[32] =   7.1993192203946599e-01
  xc[33] =   2.1690162067728006e-02;  yc[33] =   7.1993192203946599e-01;  zc[33] =   2.1690162067728006e-02
  xc[34] =   7.1993192203946599e-01;  yc[34] =   2.1690162067728006e-02;  zc[34] =   2.1690162067728006e-02
  xc[35] =   2.0448008063679571e-01;  yc[35] =   5.8057719012880926e-01;  zc[35] =   1.0462648597599378e-02
  xc[36] =   2.0448008063679571e-01;  yc[36] =   1.0462648597599378e-02;  zc[36] =   5.8057719012880926e-01
  xc[37] =   5.8057719012880926e-01;  yc[37] =   2.0448008063679571e-01;  zc[37] =   1.0462648597599378e-02
  xc[38] =   5.8057719012880926e-01;  yc[38] =   1.0462648597599378e-02;  zc[38] =   2.0448008063679571e-01
  xc[39] =   1.0462648597599378e-02;  yc[39] =   2.0448008063679571e-01;  zc[39] =   5.8057719012880926e-01
  xc[40] =   1.0462648597599378e-02;  yc[40] =   5.8057719012880926e-01;  zc[40] =   2.0448008063679571e-01
  xc[41] =   2.0448008063679571e-01;  yc[41] =   2.0448008063679571e-01;  zc[41] =   1.0462648597599378e-02
  xc[42] =   2.0448008063679571e-01;  yc[42] =   1.0462648597599378e-02;  zc[42] =   2.0448008063679571e-01
  xc[43] =   1.0462648597599378e-02;  yc[43] =   2.0448008063679571e-01;  zc[43] =   2.0448008063679571e-01
  xc[44] =   2.0448008063679571e-01;  yc[44] =   2.0448008063679571e-01;  zc[44] =   5.8057719012880926e-01
  xc[45] =   2.0448008063679571e-01;  yc[45] =   5.8057719012880926e-01;  zc[45] =   2.0448008063679571e-01
  xc[46] =   5.8057719012880926e-01;  yc[46] =   2.0448008063679571e-01;  zc[46] =   2.0448008063679571e-01

  weight[ 1] =   1.0661912966503870e-03# * 6.0
  weight[ 2] =   1.0661912966503870e-03# * 6.0
  weight[ 3] =   1.0661912966503870e-03# * 6.0
  weight[ 4] =   1.0661912966503870e-03# * 6.0
  weight[ 5] =   6.6984080034943614e-03# * 6.0
  weight[ 6] =   6.6984080034943614e-03# * 6.0
  weight[ 7] =   6.6984080034943614e-03# * 6.0
  weight[ 8] =   6.6984080034943614e-03# * 6.0
  weight[ 9] =   4.0513292507950533e-03# * 6.0
  weight[10] =   4.0513292507950533e-03# * 6.0
  weight[11] =   4.0513292507950533e-03# * 6.0
  weight[12] =   4.0513292507950533e-03# * 6.0
  weight[13] =   9.1431487356162400e-03# * 6.0
  weight[14] =   9.1431487356162400e-03# * 6.0
  weight[15] =   9.1431487356162400e-03# * 6.0
  weight[16] =   9.1431487356162400e-03# * 6.0
  weight[17] =   5.9532687056831970e-03# * 6.0
  weight[18] =   5.9532687056831970e-03# * 6.0
  weight[19] =   5.9532687056831970e-03# * 6.0
  weight[20] =   5.9532687056831970e-03# * 6.0
  weight[21] =   5.9532687056831970e-03# * 6.0
  weight[22] =   5.9532687056831970e-03# * 6.0
  weight[23] =   1.1971984496420898e-03# * 6.0
  weight[24] =   1.1971984496420898e-03# * 6.0
  weight[25] =   1.1971984496420898e-03# * 6.0
  weight[26] =   1.1971984496420898e-03# * 6.0
  weight[27] =   1.1971984496420898e-03# * 6.0
  weight[28] =   1.1971984496420898e-03# * 6.0
  weight[29] =   1.1971984496420898e-03# * 6.0
  weight[30] =   1.1971984496420898e-03# * 6.0
  weight[31] =   1.1971984496420898e-03# * 6.0
  weight[32] =   1.1971984496420898e-03# * 6.0
  weight[33] =   1.1971984496420898e-03# * 6.0
  weight[34] =   1.1971984496420898e-03# * 6.0
  weight[35] =   2.7286969908865196e-03# * 6.0
  weight[36] =   2.7286969908865196e-03# * 6.0
  weight[37] =   2.7286969908865196e-03# * 6.0
  weight[38] =   2.7286969908865196e-03# * 6.0
  weight[39] =   2.7286969908865196e-03# * 6.0
  weight[40] =   2.7286969908865196e-03# * 6.0
  weight[41] =   2.7286969908865196e-03# * 6.0
  weight[42] =   2.7286969908865196e-03# * 6.0
  weight[43] =   2.7286969908865196e-03# * 6.0
  weight[44] =   2.7286969908865196e-03# * 6.0
  weight[45] =   2.7286969908865196e-03# * 6.0
  weight[46] =   2.7286969908865196e-03# * 6.0

  pts = hcat( xc, yc, zc )

  return pts, weight

end


text:
using AstroUtils
using Test, SafeTestsets

@time begin
@time @safetestset "cartToKep tests..." begin include("cartToKepTests.jl") end
end


text:
#=
This file is part of Pseudospectra.jl., whose LICENSE file applies.

Julia translation
copyright (c) 2017 Ralph Smith

Portions derived from EigTool
Copyright (c) 2002-2014, The Chancellor, Masters and Scholars
of the University of Oxford, and the EigTool Developers. All rights reserved.

SPDX-License-Identifier: BSD-3-Clause
License-Filename: LICENSES/BSD-3-Clause_Eigtool
=#

"""
    modeplot(ps_data, pkey [,z])

Extract and plot an eigenmode or pseudo-eigenmode for the matrix
encapsulated in the Pseudospectra object `ps_data`.
Use the value `z` if provided or prompt for one.
If `pkey` is 1, find the pseudoeigenmode for `z`;
otherwise find the eigenmode for the eigenvalue closest to `z`.
"""
function modeplot(ps_data::PSAStruct, pkey, z=NaN; gs=defaultgs(),
                  zgetter=replzdlg, zmarker=addmark, verbosity = 0)
    # based on switch_pseudomode()
    # FIXME: a lot of logic is missing here
    ps_dict = ps_data.ps_dict
    pseudo = (pkey == 1) # in lieu of Bool arg FSO GUI
    # get a point
    if isnan(z)
        z = zgetter(gs, what=(pseudo ? "a z-value" : "approx. eigenvalue"))
    end
    if isnan(z)
        warn("modeplot: invalid z")
        return
    end
    z = complex(z)
    if !pseudo
        eigA = get(ps_dict,:proj_ews,ps_dict[:ews])
        isempty(eigA) && throw(ArgumentError("ps_data holds no eigenvalues"))
        dists = [abs(eig - z) for eig in eigA]
        idx = argmin(dists)
        z = eigA[idx]
    end
    if zmarker != nothing
        zmarker(gs,z,(pseudo ? :pseudo : :eigen))
    end
    # FIXME: check for direct method, etc.
    if !isempty(get(ps_dict,:schur_mtx,zeros(0)))
        A = ps_dict[:schur_mtx]
        U = ps_dict[:schur_unitary_mtx]
    elseif !isempty(get(ps_dict,:proj_matrix,zeros(0))) # && etc.
            A = ps_dict[:proj_matrix]
            U = ps_dict[:proj_unitary_mtx]
    # elseif sparse_direct (no Schur)
    else
        A = ps_data.input_matrix
        U = ps_data.input_unitary_mtx
    end
    # in Arnoldi case extract upper square mtx
    mm,nn = size(A)
    if mm == nn
        plotmode(gs,z,A,U,pseudo,false,verbosity) # avoid pointless copy
    else
        plotmode(gs,z,A[1:nn,1:nn],U,pseudo,false,verbosity)
    end
    # TODO: store marker info in gs
    nothing
end

"""
    oneeigcond(T,ew,verbosity;dlg=replqdlg) -> cond,evr,evl

Obtain condition and eigenvector(s) for eigenvalue `ew` associated with
Schur factor `T`.
"""
function oneeigcond(T0,ew0,verbosity; dlg=replqdlg, max_its=3)
    quest_str = ("The eigenvector could not be determined by inverse iteration."
                 * "\nDo you want to do a (possibly slow) dense SVD instead?")
    n,m = size(T0)
    @assert n==m
    tol = 1e-14
    # perturb eigenvalue to avoid exact singularity
    ew = ifelse(n<200,ew0,ew0*(1+2*eps(abs(ew0))))
    T = copy(T0)
    T = T - ew * I
    # for small matrices, use SVD
    if n < 200
        U,S,V = svd(Matrix(T))
        evr = V[:,end]
        evl = conj(U[:,end])
    else
        # inverse iteration
        if issparse(T)
            F = lufact(T)
        end
        wb_size = 200
        v0 = normalize!(randn(n) .+ 0im)
        infs_found = false # indicator for inv-iter success
        evr = copy(v0)
        local resid
        for i=1:max_its
            oldev = copy(evr)
            try
                if issparse(T)
                    evr = F \ oldev
                else
                    evr = T \ oldev
                end
            catch JE
                # FIXME: check exception type: LAPACKException or ??
                infs_found = true
            end
            any(isinf.(evr)) && (infs_found = true)
            if infs_found
                break
            end
            λ = (oldev' * evr)[1]
            resid = norm(evr - λ * oldev) / abs(λ)
            normalize!(evr)
            (verbosity > 2) && println("evr resid: $resid")
            if resid < tol
                break
            end
        end
        if !(infs_found || (resid > tol))
            # no problem w/ right ev, so do left one
            evl = v0
            for i=1:max_its
                oldev = copy(evl)
                try
                    if issparse(T)
                        ldiv!(evl,transpose(F),oldev)
                    else
                        evl = transpose(T) \ oldev # At_ldiv_B!(evl,T,oldev)
                    end
                catch JE
                    infs_found = true
                end
                any(isinf.(evl)) && (infs_found = true)
                if infs_found
                    break
                end
                λ = (oldev' * evl)[1]
                resid = norm(evl - λ * oldev) / abs(λ)
                normalize!(evl)
                (verbosity > 2) && println("evl resid: $resid")
                if resid < tol
                    break
                end
            end
        end
        if infs_found || (resid > tol)
            (verbosity > 1) && println("Infs found: $infs_found resid: $resid")
            if dlg("Compute using full SVD?",quest_str) == 1 # 1:yes
                U,S,V = svd(Matrix(T))
                evr = V[:,end]
                evl = conj(U[:,end])
            else
                evr = NaN
                evl = NaN
            end
        end
    end
    z0 = transpose(evl)*evr
    the_cond = abs(1/z0[1])
    return the_cond, evr, evl
end

"""
     psmode_inv_lanczos(A,q,z,tol,num_its) -> Z,q

use inverse-Lanczos iteration to compute a pseudomode

 A         the matrix to use to find min σ(zI-A)
 q         starting vector for the iteration
 z         point in the complex plane to use
 tol       tolerance to use
 num_its   Maximum number of iterations

 Z         the smallest singular values
 q         the smallest singular vector
"""
function psmode_inv_lanczos(A,q0::Vector,z::Complex,tol,num_its;
                            fallbacksvd=false)
    m,n = size(A)
    (m == n) || throw(ArgumentError("Matrix must be square"))
    diaga = diag(A)
    cdiaga = conj(diaga)
    q = q0 .+ 0im
    if n < 200
        # just use plain SVD
        A1 = copy(A) .+ 0im
        A1[1:m+1:end] = diaga .- z
        U,S,V = svd(Matrix(A1))
        Z = minimum(S)
        @assert S[end] == Z
        q = V[:,end]
    else
        H = zeros(num_its+1,num_its+1)
        if issparse(A)
            F = lufact(A-z*I)
            # use UMFPACK factor method to avoid this:
            # Fp = lufact(A'-conj(z)*I)
        else
            T1 = A - z * I
            T2 = T1'
        end
        qold = zeros(n)
        Q = copy(q0)
        β = 0.0
        sigold = 0.0
        local σ,α
        l = 0
        for it=1:num_its
            l += 1
            if issparse(A)
                # v = F \ (Fp \ q) - β * qold
                w = zeros(q+0im)
                v = similar(w)
                ldiv!(w,adjoint(F),q)
                ldiv!(v,F,q)
                v -= β * qold
            else
                v = T1 \ (T2 \ q) - β * qold
            end
            α = real((q' * v)[1])
            v = v - α * q
            β = norm(v)
            qold = q
            q = v / β
            Q = hcat(Q,q)
            H[l+1,l] = β
            H[l,l+1] = β
            H[l,l] = α
            # calculate eigenvalues of H, but if error is too big,
            # set σ to a huge value
            try
                E = eigen(H[1:l,1:l])
                σ = maximum(E.values)
            catch
                σ = 1e308
            end
            resid = abs(sigold/σ - 1)
            if (resid < tol) || (β == 0)
                break
            end
            sigold = σ
        end
        H = H[1:l,1:l]
        Z = 1 / sqrt(σ)
        # get eigenvector from inverse Lanczos basis
        # if matrix is too non-normal, revert
        try
            E = eigen(H)
            q = Q[:,1:end-1] * E.vectors[:,end]
        catch
            l = num_its
        end
        if l >= num_its
            if fallbacksvd
                A1 = A - z * I
                U,S,V = svd(Matrix(A1))
                Z = minimum(S)
                q = V[:,end]
            else
                Z = NaN
                q = NaN
            end
        end
    end # if n < 200
    return Z,q
end


text:
if(Pkg.installed("PyPlot") == nothing)
    Pkg.add("PyPlot")
end
using PyPlot
close("all")
if !isfile("gravity/gravity_bc2016.oifits")
    if !isdir("gravity")
        mkdir("gravity")
    end
    download("http://www.opticalinterferometry.com/s/Object1-v3.oifits","gravity/gravity_bc2016.oifits")
end

function plotfunction(PDATA::PAINTER.PAINTER_Data,OIDATA::PAINTER.PAINTER_Input)
    nx = OIDATA.nx
    nw = OIDATA.nw
    wvl = OIDATA.wvl
    FOV = OIDATA.FOV
    x = PDATA.x
    w = PDATA.w
    if sum(w)==0
      w = 1.
    end
    indpix = linspace(-(FOV / 2), (FOV / 2), nx)
    pos = round(Int, [1, round(Int,nx / 4), round(Int,nx / 2), round(Int,nx * 3 / 4), nx])

    count_y = 0
    count_x = 0
    SubRow  = 5
    SubColumn = 5

    X3D = x .* max(0, w) .* max(x,0)
    VMIN = minimum(X3D)
    VMAX = maximum(X3D)

    figure("3D Image Reconstruction ")
    if nw == 250
      set = collect(1:10:250)
    elseif nw == 25
      set = collect(1:25)
    elseif nw == 50
      set = collect(1:2:50)
    end
    n2 = 0
    for n in set
        n2+=1
        subplot(SubColumn, SubRow, n2)
        imshow( sqrt( rotr90(X3D[:, :, n]))  , origin = "lower")
        titlestring = @sprintf("%2.4f µm", wvl[n2] * 1e6)
        title(titlestring)
        xticks([])
        yticks([])

        if( n2 == (length(set) + 1 - SubRow + count_x) )
            xticks(collect(pos - 1), round(Int,indpix[pos] * 1000) )
            xlabel("FOV (mas)")
            count_x += 1
        end
        if( n2 == (1 + count_y * SubRow))
            yticks(collect(pos - 1), round(Int,indpix[pos] * 1000) )
            ylabel("FOV (mas)")
            count_y += 1
        end
    end

    figure("Projection")

    clf()

    subplot(222); imshow(sqrt(squeeze(mean(X3D,1),1)),origin = "lower",aspect="auto");
    yticks([])
    xlabel("channels")
    title("Spectrum 2")

    subplot(221); imshow(sqrt(rotr90(squeeze(mean(X3D,3),3))),origin = "lower",aspect="auto");
    yticks(collect(pos - 1), round(Int,indpix[pos] * 1000) )
    xticks([])
    ylabel("FOV (mas)")
    title("Gray")

    subplot(223); imshow( sqrt(rotr90(squeeze(mean(X3D,2) ,2))),aspect="auto");
    xticks(collect(pos - 1), round(Int,indpix[pos] * 1000) )
    xlabel("FOV (mas)")
    ylabel("channels")
    title("spectrum 1")

end
    # To change size of the simulation
    # nx pixels
    # files contains 250 wavelength, the set can be reduced to collect(1:10:250)
    # the plot function is done for 25 wavelengths
    # choose collect(1:10:250) or collect(1:250) or collect(1:25)

    println("")
    println("-----------------------------------------")
    println("| Gravity - DEMO - Beauty Constest 2016 |")
    println("-----------------------------------------")
    println("")

    nx = 64
    indwvl = 1:5:250 # will plot 1/10

    Folder = "gravity"
    savepath = "gravity_bc2016.jld"

    dptype = "phase"

    CountPlot = 10
    nbitermax = 400

    aff = true
    admm = true

    PlotFct = plotfunction

    FOV = 0.05

    rho_y = 10.
    rho_spat = 4.
    rho_spec = .5
    rho_ps =  rho_spat

    epsilon = 1e-6

    eps1 = 1e-3
    eps2 = 1e-3

    mask3D = PAINTER.mask(nx,round(Int, nx/2 - 2 ))
    xinit3D = PAINTER.mask(nx,round(Int, nx/2  ))

# initialize algorithm and run admm
    OIDATA, PDATA = PAINTER.painter(nbitermax = nbitermax, nx = nx,
      rho_y = rho_y, rho_spat = rho_spat,rho_spec = rho_spec, rho_ps = rho_ps,
      eps1 = eps1, eps2 = eps2, FOV = FOV, indwvl = indwvl,
      PlotFct = PlotFct, aff = aff, dptype = dptype, admm = admm,
      xinit3D = xinit3D, mask3D = mask3D, Folder = Folder)

# save data struture in .jld files
    println("save results of gravity BC2016 data")
    PAINTER.paintersave(savepath,PDATA,OIDATA)

# load data struture in .jld files
    println("load results of gravity BC2016 data")
    PDATA, OIDATA = PAINTER.painterload(savepath)
    nothing


text:
module CatIndices

export BidirectionalVector, PinIndices, deletetail!, deletehead!

using OffsetArrays

include("bidirectional.jl")
include("cat.jl")

end # module


text:
### Utilities

export compute_critical_dt, compute_abcsize

function extend(m::Model, abc_geom::ABC_Geom; opt_rep::Bool = true)

	return extend(m.m, abc_geom; opt_rep = opt_rep)

end

function extend(m::Array{Float32, 2}, abc_geom::ABC_Geom; opt_rep::Bool = true) # Array size: nz, nb

	# Initialize
	m_ext = Array{Float32, 2}(undef, size(m, 1)+sum(abc_geom.size), size(m, 2))

	# Setting values on regular domain
	m_ext[abc_geom.size[1]+1:end-abc_geom.size[2], :] .= m

	# Setting values on extended domain
	if opt_rep
		m_ext[1:abc_geom.size[1],         :] .= m[1:1,     :]
		m_ext[end-abc_geom.size[2]+1:end, :] .= m[end:end, :]
	else
		m_ext[1:abc_geom.size[1],         :] .= 0f0
		m_ext[end-abc_geom.size[2]+1:end, :] .= 0f0
	end

	return m_ext

end

function taper(abc_geom::ABC_Geom)

	# Initialize
	tap_fcn = Array{Float32, 1}(undef, abc_geom.nz+sum(abc_geom.size))

	# Values on regular domain
	tap_fcn[abc_geom.size[1]+1:end-abc_geom.size[2]] .= 0f0

	# Values on extended domain
	tap_fcn[1:abc_geom.size[1]] = range(-1f0, stop = 0f0, length = abc_geom.size[1]).^2f0
	tap_fcn[end-abc_geom.size[2]+1:end] = range(0f0, stop = 1f0, length = abc_geom.size[2]).^2f0

	return tap_fcn*abc_geom.fact

end

function compute_critical_dt(m::Model, dat_pars::DataParams)

	return 1.73f0/2f0*m.geom.dz*sqrt(min(m.m...)-max(dat_pars.p...)^2f0)

end

function compute_abcsize(m::Model, dat_pars::DataParams, freq::Float32)

	return Int64(ceil(1f0/(sqrt(min(m.m...)-max(dat_pars.p...)^2f0)*freq*m.geom.dz)))

end


text:
const DEFAULT_JOB_NAME = ""

const DEFAULT_DT_NULL = DateTime(0)
const DEFAULT_DT_NEXT_FIRE = DEFAULT_DT_NULL

"""
    JobConfig([misfire_grace_period])

Job default configuration

# Optional arguments
- misfire_grace_period (defaults to Dates.Second(1)): grace period for which a task can still be fired.
"""
mutable struct JobConfig
    misfire_grace_period

    JobConfig(; misfire_grace_period=Dates.Second(1)) = new(misfire_grace_period)
end

"""
    Job(id, action, trigger, name, priority, dt_created, dt_updated, dt_next_fire, n_triggered, config)

A job is an internal structure which store what action
should be executed when triggered.

It also store several properties such as priority level, number of time a job is triggered, when will next trigger should occur...
"""
mutable struct Job
    id::String
    action::Action
    trigger::AbstractTrigger
    name::String
    priority::Int
    dt_created::DateTime
    dt_updated::DateTime
    #dt_last_fire::DateTime
    dt_next_fire::DateTime
    n_triggered::BigInt
    config::JobConfig
end

"""
    hasnextfire(job)

Returns `true` if a task should be fired in the future.
"""
function hasnextfire(job::Job)
    job.dt_next_fire != DEFAULT_DT_NEXT_FIRE
end


text:
module LongwaveModePropagator

import Base: @_inline_meta, @propagate_inbounds, @_propagate_inbounds_meta
import Base: ==
using Dates, LinearAlgebra

using StaticArrays
using StaticArrays: promote_tuple_eltype, convert_ntuple
using OrdinaryDiffEq, DiffEqCallbacks
using Parameters, ProgressLogging
using JSON3, StructTypes
using Interpolations
import FunctionWrappers: FunctionWrapper

using PolynomialRoots: roots!

using RootsAndPoles, Romberg, ModifiedHankelFunctionsOfOrderOneThird
export GRPFParams

# LongwaveModePropagator.jl
export propagate

# EigenAngle.jl
export EigenAngle
export attenuation, phasevelocity, referencetoground, setea

# Geophysics.jl
export BField, Species, Fields, Ground, GROUND
export waitprofile, electroncollisionfrequency, ioncollisionfrequency

# IO.jl
export ExponentialInput, TableInput, BatchInput, BasicOutput, BatchOutput

# Samplers.jl
export Receiver, Sampler, GroundSampler

# Emitters.jl
export Transmitter, Dipole, VerticalDipole, HorizontalDipole, Frequency
export inclination, azimuth

# modefinder.jl
export findmodes, PhysicalModeEquation

# Waveguides.jl
export HomogeneousWaveguide, SegmentedWaveguide

"""
    IntegrationParams{T}

Parameters passed to `OrdinaryDiffEq.jl` during the integration of the ionosphere reflection
coefficient matrix in `modefinder.jl`.

# Fields

- `tolerance::Float64 = 1e-5`: integration `atol` and `rtol`.
- `solver::T = Vern7()`: a `DifferentialEquations.jl` solver.
- `dt::Float64 = 1.0`: height step in meters (many methods use a variable step size).
- `force_dtmin::Bool = false`: if true, continue integration when solver reaches minimum
    step size.
- `maxiters::Int = 100_000`: maximum number of iterations before stopping.
"""
@with_kw struct IntegrationParams{T}
    tolerance::Float64 = 1e-5
    solver::T = Vern7()
    dt::Float64 = 1.0
    force_dtmin::Bool = false
    maxiters::Int = 100_000
end
export IntegrationParams

const DEFAULT_GRPFPARAMS = GRPFParams(100000, 1e-5, true)

"""
    LMPParams{T,T2,H <: AbstractRange{Float64}}

Parameters for the `LongwaveModePropagator` module with defaults:

- `topheight::Float64 = 110e3`: starting height for integration of the ionosphere reflection
    coefficient.
- `earthradius::Float64 = 6369e3`: Earth radius in meters.
- `earthcurvature::Bool = true`: toggle inclusion of Earth curvature in calculations. This
    is not supported by all functions.
- `curvatureheight::Float64 = 50e3`: reference height for Earth curvature in meters. At this
    height, the index of refraction is 1, and is therefore the reference height for
    eigenangles.
- `approxsusceptibility::Bool = false`: use a cubic interpolating spline representation of
    [`susceptibility`](@ref) during the integration of [`dRdz`](@ref).
- `susceptibilitysplinestep::Float64 = 10.0`: altitude step in meters used to build the
    spline representation of [`susceptibility`](@ref) if `approxsusceptibility == true`.
- `grpfparams::GRPFParams = GRPFParams(100000, 1e-5, true)`: parameters for the `GRPF`
    complex root-finding algorithm.
- `integrationparams::IntegrationParams{T} =
    IntegrationParams(solver=Vern7(), tolerance=1e-5)`:
    parameters passed to `DifferentialEquations.jl` for integration of the ionosphere
    reflection coefficient.
- `wavefieldheights::H = range(topheight, 0, length=513)`: heights in meters at which
    wavefields will be integrated.
- `wavefieldintegrationparams::IntegrationParams{T2} =
    IntegrationParams(solver=Tsit5(), tolerance=1e-6)`:
    parameters passed to `DifferentialEquations.jl` for integration of the wavefields
    used in mode conversion. The solver cannot be lazy.

The struct is created using `Parameters.jl` `@with_kw` and supports that package's
instantiation capabilities, e.g.:

```julia
p = LMPParams()
p2 = LMPParams(earth_radius=6370e3)
p3 = LMPParams(p2; grpf_params=GRPFParams(100000, 1e-6, true))
```

See also: [`IntegrationParams`](@ref)
"""
@with_kw struct LMPParams{T,T2,H<:AbstractRange{Float64}}
    topheight::Float64 = 110e3
    earthradius::Float64 = 6369e3  # m
    earthcurvature::Bool = true
    curvatureheight::Float64 = 50e3  # m
    approxsusceptibility::Bool = false
    susceptibilitysplinestep::Float64 = 10.0
    grpfparams::GRPFParams = DEFAULT_GRPFPARAMS
    integrationparams::IntegrationParams{T} = IntegrationParams()
    wavefieldheights::H = range(topheight, 0; length=513)
    wavefieldintegrationparams::IntegrationParams{T2} =
        IntegrationParams(solver=Tsit5(), tolerance=1e-6)
end
export LMPParams

"""
    ModeEquation

Functions can dispatch on subtypes of `ModeEquation`, although currently only
`PhysicalModeEquation` is supported.

Future work might include the `ModifiedModeEquation` of [Morfitt1976].

# References

[Morfitt1976]: D. G. Morfitt and C. H. Shellman, “‘MODESRCH’, an improved computer program
    for obtaining ELF/VLF/LF mode constants in an Earth-ionosphere waveguide,” Naval
    Electronics Laboratory Center, San Diego, CA, NELC/IR-77T, Oct. 1976. [Online].
    Available: http://www.dtic.mil/docs/citations/ADA032573.
"""
abstract type ModeEquation end  # concrete subtypes defined in modefinder files

"Indicates derivative with respect to eigenangle ``θ``."
struct Dθ end

"""
Bottom of reflection coefficient integration. `LongwaveModePropagator.jl` assumes
`BOTTOMHEIGHT = 0.0`.
"""
const BOTTOMHEIGHT = 0.0  # WARNING: if this isn't zero, many assumptions break

#
include("utils.jl")
include("Antennas.jl")
include("Emitters.jl")  # must be before EigenAngles.jl
include("EigenAngles.jl")
include("Geophysics.jl")
include("Samplers.jl")
include("TMatrix.jl")
include("Waveguides.jl")
include("Wavefields.jl")

include("bookerquartic.jl")
include("magnetoionic.jl")
include("modeconversion.jl")
include("modefinder.jl")
include("modesum.jl")

include("IO.jl")

"""
    propagate(waveguide::HomogeneousWaveguide, tx::Emitter, rx::AbstractSampler;
              modes::Union{Nothing,Vector{EigenAngle}}=nothing, mesh=nothing,
              params=LMPParams())

Compute electric field `E`, `amplitude`, and `phase` at `rx`.

Precomputed waveguide `modes` can optionally be provided as a `Vector{EigenAngle}`. By
default modes are found with [`findmodes`](@ref).

If `mesh = nothing`, use [`defaultmesh`](@ref) to generate `mesh` for the
mode finding algorithm. This is ignored if `modes` is not `nothing`.
"""
function propagate(waveguide::HomogeneousWaveguide, tx::Emitter, rx::AbstractSampler;
    modes::Union{Nothing,Vector{EigenAngle}}=nothing, mesh=nothing, unwrap=true,
    params=LMPParams())

    if isnothing(modes)
        if isnothing(mesh)
            mesh = defaultmesh(tx.frequency)
        else
            if minimum(imag(mesh)) < deg2rad(-31)
                @warn "imaginary component less than -0.5410 rad (-31°) may cause wave"*
                    "fields calculated with modified Hankel functions to overflow."
            end
        end

        modeequation = PhysicalModeEquation(tx.frequency, waveguide)
        modes = findmodes(modeequation, mesh; params=params)
    end

    E = Efield(modes, waveguide, tx, rx; params=params)

    amplitude, phase = amplitudephase(E)

    if unwrap
        unwrap!(phase)
    end

    return E, amplitude, phase
end

"""
    propagate(waveguide::SegmentedWaveguide, tx::Emitter, rx::AbstractSampler;
              mesh=nothing, params=LMPParams())

Compute electric field `E`, `amplitude`, and `phase` at `rx` through a `SegmentedWaveguide`.

If `mesh = nothing`, use [`defaultmesh`](@ref) to generate `mesh` for the
mode finding algorithm.
"""
function propagate(waveguide::SegmentedWaveguide, tx::Emitter, rx::AbstractSampler;
    mesh=nothing, unwrap=true, params=LMPParams())

    if isnothing(mesh)
        mesh = defaultmesh(tx.frequency)
    else
        if minimum(imag(mesh)) < deg2rad(-31)
            @warn "imaginary component less than -0.5410 rad (-31°) may cause wave fields"*
                "calculated with modified Hankel functions to overflow."
        end
    end

    J = length(waveguide)  # number of waveguide segments

    heighttype = typeof(params.wavefieldheights)
    wavefields_vec = Vector{Wavefields{heighttype}}(undef, J)
    adjwavefields_vec = Vector{Wavefields{heighttype}}(undef, J)

    # Calculate wavefields and adjoint wavefields for each segment of waveguide
    for j in 1:J
        wvg = waveguide[j]

        modeequation = PhysicalModeEquation(tx.frequency, wvg)

        modes = findmodes(modeequation, mesh; params=params)

        # adjoint wavefields are wavefields through adjoint waveguide, but for same modes
        # as wavefield
        adjwvg = adjoint(wvg)

        wavefields = Wavefields(params.wavefieldheights, modes)
        adjwavefields = Wavefields(params.wavefieldheights, modes)

        calculate_wavefields!(wavefields, adjwavefields, tx.frequency, wvg, adjwvg;
                              params=params)

        wavefields_vec[j] = wavefields
        adjwavefields_vec[j] = adjwavefields
    end

    E = Efield(waveguide, wavefields_vec, adjwavefields_vec, tx, rx; params=params)

    # Efield for SegmentedWaveguides doesn't have a specialized form for AbstractSamplers
    # of Number type, but for consistency we will return scalar E.
    # Although now this function is technically type unstable, it has a practically
    # unmeasurable affect on the total runtime.
    if length(E) == 1
        E = only(E)
    end

    amplitude, phase = amplitudephase(E)

    if unwrap
        unwrap!(phase)
    end

    return E, amplitude, phase
end

"""
    propagate(file::AbstractString, outfile=missing; incrementalwrite=false, append=false,
              mesh=nothing)

Run the model scenario described by `file` and save the results as `outfile`.

If `outfile = missing`, the output file name will be `\$(file)_output.json`.
"""
function propagate(file::AbstractString, outfile=missing; incrementalwrite=false,
    append=false, mesh=nothing, unwrap=true, params=LMPParams())

    ispath(file) || error("$file is not a valid file name")

    if ismissing(outfile)
        basepath = dirname(file)
        filename, fileextension = splitext(basename(file))

        outfile = joinpath(basepath, filename)*"_output.json"
    end

    s = parse(file)
    if incrementalwrite
        s isa BatchInput || throw(ArgumentError("incrementalwrite only supported for"*
                                                "BatchInput files"))
        output = buildrunsave(outfile, s; append=append, mesh=mesh, unwrap=unwrap, params=params)
    else
        output = buildrun(s; mesh=mesh, unwrap=unwrap, params=params)

        json_str = JSON3.write(output)

        open(outfile, "w") do f
            write(f, json_str)
        end
    end

    return output
end

end


text:
@userplot Ternary_Axes

@recipe function f(ta::Ternary_Axes)
    # empty plot without axes
    framestyle := :none
    aspect_ratio := true
    dist_from_graph = get(plotattributes, :arrow_length, 0.04)
    arrow_length = get(plotattributes, :arrow_length, 0.4)
    arrow_pos = (1 - arrow_length) / 2, 1 - (1 - arrow_length) / 2
    tick_length = get(plotattributes, :tick_length, 0.015)

    # calculate where arrows and labels should be
    x_arrow_src, x_arrow_dst, x_guide_pos = calculate_guide_pos("x", arrow_pos, dist_from_graph)
    y_arrow_src, y_arrow_dst, y_guide_pos = calculate_guide_pos("y", arrow_pos, dist_from_graph)
    z_arrow_src, z_arrow_dst, z_guide_pos = calculate_guide_pos("z", arrow_pos, dist_from_graph)

    x_arrow_xs, x_arrow_ys = coords_to_axes(x_arrow_src, x_arrow_dst)
    y_arrow_xs, y_arrow_ys = coords_to_axes(y_arrow_src, y_arrow_dst)
    z_arrow_xs, z_arrow_ys = coords_to_axes(z_arrow_src, z_arrow_dst)
    arrow_axes_xs_yes = [(x_arrow_xs, x_arrow_ys), (y_arrow_xs, y_arrow_ys), (z_arrow_xs, z_arrow_ys)]

    if !get(plotattributes, :axis_arrows, true)
        x_guide_pos = (x_arrow_src .+ x_arrow_dst) ./ 2
        y_guide_pos = (y_arrow_src .+ y_arrow_dst) ./ 2
        z_guide_pos = (z_arrow_src .+ z_arrow_dst) ./ 2
    else
        # draw axis arrows
        for xs_ys in arrow_axes_xs_yes
            @series begin
                seriestype := :path
                primary := false
                seriescolor := :black
                arrow := true
                xs_ys
            end
        end
    end

    # labels
    xguide = get(plotattributes, :xguide, "")
    @series begin
        seriestype := :path
        primary := false
        series_annotations := [(xguide)] #[text(xguide),]
        tupvec(x_guide_pos)
    end

    yguide = get(plotattributes, :yguide, "")
    @series begin
        seriestype := :path
        primary := false
        series_annotations := [(yguide, -60.0)]
        tupvec(y_guide_pos)
    end

    zguide = get(plotattributes, :zguide, "")
    @series begin
        seriestype := :path
        primary := false
        series_annotations := [(zguide, 60.0)] #[text(zguide, 60.0),]
        tupvec(z_guide_pos)
    end

    # grids
    if get(plotattributes, :gridmajor, true)
        if get(plotattributes, :xgrid, true)
            ticks = get(plotattributes, :xticks, 0.2:0.2:0.8)
            src = tern2cart.(1 .- ticks, ticks, 0)
            dst = tern2cart.(1 .- ticks, 0, ticks)
            sa = get(plotattributes, :xgridalpha, 0.5)
            lw = get(plotattributes, :xgridlinewidth, 0.5)
            ls = get(plotattributes, :xgridstyle, :solid)
            for i in eachindex(ticks)
                @series begin
                    primary := false
                    seriestype := :path
                    seriescolor := :black
                    seriesalpha := sa
                    linewidth := lw
                    linestyle := ls
                    arrow := false
                    [src[i], dst[i]]
                end
            end
        end

        if get(plotattributes, :ygrid, true)
            ticks = get(plotattributes, :yticks, 0.2:0.2:0.8)
            src = tern2cart.(1 .- ticks, ticks, 0)
            dst = tern2cart.(0, ticks, 1 .- ticks)
            sa = get(plotattributes, :ygridalpha, 0.5)
            lw = get(plotattributes, :ygridlinewidth, 0.5)
            ls = get(plotattributes, :ygridstyle, :solid)
            for i in eachindex(ticks)
                @series begin
                    primary := false
                    seriestype := :path
                    seriescolor := :black
                    seriesalpha := sa
                    linewidth := lw
                    linestyle := ls
                    [src[i], dst[i]]
                end
            end
        end

        if get(plotattributes, :zgrid, true)
            ticks = get(plotattributes, :zticks, 0.2:0.2:0.8)
            src = tern2cart.(1 .- ticks, 0, ticks)
            dst = tern2cart.(0, 1 .- ticks, ticks)
            sa = get(plotattributes, :zgridalpha, 0.5)
            lw = get(plotattributes, :zgridlinewidth, 0.5)
            ls = get(plotattributes, :zgridstyle, :solid)
            for i in eachindex(ticks)
                @series begin
                    primary := false
                    seriestype := :path
                    seriescolor := :black
                    seriesalpha := sa
                    linewidth := lw
                    linestyle := ls
                    [src[i], dst[i]]
                end
            end
        end
    end

    tick_direction = get(plotattributes, :tick_direction, :in)
    # axis ticks
    if tick_direction != :none
        if tick_direction == :in
            tickmult = -1
        else
            tickmult = 1
        end
        xticks = get(plotattributes, :xticks, 0.2:0.2:0.8)
        yticks = get(plotattributes, :yticks, 0.2:0.2:0.8)
        zticks = get(plotattributes, :zticks, 0.2:0.2:0.8)

        function text_structure(axis)
            if axis == "x"
                return 0.4, -0.8
            elseif axis == "y"
                return -1, 0.2
            elseif axis == "z"
                return 0.5, 1
            else
                throw(ArgumentError("Invalid axis (\"x, y, z\""))
            end
        end

        function axis_angle(axis)
            if axis == "x"
                return 0.0
            elseif axis == "y"
                return -60.0
            elseif axis == "z"
                return 60.0
            else
                throw(ArgumentError("Invalid axis (\"x, y, z\""))
            end
        end

        for (axis, axis_ticks) in zip(("x", "y", "z"), (xticks, yticks, zticks))
            angle = axis_angle(axis)
            ts = text_structure(axis)
            # Better to calculate the offset once by calculating the tick_structure once
            # then just adding the offset to each start_pos
            calculated_ticks = calculate_ticks(axis, axis_ticks, tick_length, tickmult)
            vec_tick_textpos = [start_pos .+ (dist_from_graph .* ts) for start_pos in calculated_ticks[1]]

            vec_tick_coords = coords_to_axes.(calculated_ticks...)
            if axis == "x"
                axis_ticks = reverse(axis_ticks)
            end
            for tick_coords in vec_tick_coords
                @series begin
                    primary := false
                    seriestype := :path
                    seriescolor := :black
                    tick_coords
                end
            end
            for (tick, tick_textpos) in zip(collect(axis_ticks), vec_tick_textpos)
                @series begin
                    seriestype := :path
                    primary := false
                    series_annotations := [(tick, angle, 10)] #[text(tick, angle, 10),]
                    [tick_textpos]
                end
            end
        end
    end

    size --> (580, 550)
    seriestype := :path
    primary := false
    seriescolor := :black
    return [1 0 1/2; 0 1/2 1], [0 0 √3/2; 0 √3/2 0]
end


text:
pybytes_fromdata(x::Ptr, n::Integer) = pynew(errcheck(C.PyBytes_FromStringAndSize(x, n)))
pybytes_fromdata(x) = pybytes_fromdata(pointer(x), sizeof(x))

"""
    pybytes(x)

Convert `x` to a Python `bytes`.
"""
pybytes(x) = pynew(errcheck(@autopy x C.PyObject_Bytes(getptr(x_))))
pybytes(x::Vector{UInt8}) = pybytes_fromdata(x)
pybytes(x::Base.CodeUnits{UInt8, String}) = pybytes_fromdata(x)
pybytes(x::Base.CodeUnits{UInt8, SubString{String}}) = pybytes_fromdata(x)
pybytes(::Type{T}, x) where {Vector{UInt8} <: T <: Vector} = (b=pybytes(x); ans=pybytes_asvector(b); pydel!(b); ans)
pybytes(::Type{T}, x) where {Base.CodeUnits{UInt8,String} <: T <: Base.CodeUnits} = (b=pybytes(x); ans=Base.CodeUnits(pybytes_asUTF8string(b)); pydel!(b); ans)
export pybytes

pyisbytes(x) = pytypecheckfast(x, C.Py_TPFLAGS_BYTES_SUBCLASS)

function pybytes_asdata(x::Py)
    ptr = Ref(Ptr{Cchar}(0))
    len = Ref(C.Py_ssize_t(0))
    errcheck(C.PyBytes_AsStringAndSize(getptr(x), ptr, len))
    ptr[], len[]
end

function pybytes_asvector(x::Py)
    ptr, len = pybytes_asdata(x)
    unsafe_wrap(Array, Ptr{UInt8}(ptr), len)
end

function pybytes_asUTF8string(x::Py)
    ptr, len = pybytes_asdata(x)
    unsafe_string(Ptr{UInt8}(ptr), len)
end

pyconvert_rule_bytes(::Type{Vector{UInt8}}, x::Py) = pyconvert_return(copy(pybytes_asvector(x)))
pyconvert_rule_bytes(::Type{Base.CodeUnits{UInt8,String}}, x::Py) = pyconvert_return(codeunits(pybytes_asUTF8string(x)))

pyconvert_rule_fast(::Type{Vector{UInt8}}, x::Py) = pyisbytes(x) ? pyconvert_rule_bytes(Vector{UInt8}, x) : pyconvert_unconverted()
pyconvert_rule_fast(::Type{Base.CodeUnits{UInt8,String}}, x::Py) = pyisbytes(x) ? pyconvert_rule_bytyes(Base.CodeUnits{UInt8,String}, x) : pyconvert_unconverted()


text:
module Spectroscopy

const VecI  = AbstractVector # Input  Vector
const VecO  = AbstractVector # Output Vector
const VecB  = AbstractVector # Buffer Vector
const VecIO = AbstractVector # In/Out Vector
const MatI  = AbstractMatrix # Input  Matrix
const MatO  = AbstractMatrix # Output Matrix
const MatB  = AbstractMatrix # Buffer Matrix
const MatIO = AbstractMatrix # In/Out Matrix

using DelimitedFiles

abstract type TimeResolvedData             end
abstract type TimeAxis <: TimeResolvedData end

macro def(genre::Symbol, ex::Union{Expr,Symbol}, vars::Symbol...)
    o = genre == :prop ? :(::) : genre == :vars ? :(=) : error("@def(genre = $genre, ...) is invalid.")
    n = length(vars)
    e = Vector{Expr}(undef, n)
    @inbounds for i in 1:n
        e[i] = Expr(o, vars[i], ex)
    end
    return Expr(:escape, Expr(:block, e...))
end

macro get(obj::Symbol, vars::Symbol...)
    n = length(vars)
    e = Vector{Expr}(undef, n)
    @inbounds for i in 1:n
        vari = vars[i]
        e[i] = :($vari = $obj.$vari)
    end
    return Expr(:escape, Expr(:block, e...))
end

#### ceiling `x` to a power-of-2 integer
function clp2(x::Int)
    x == 0 && return 1
    x == 1 && return 2
    x = x - 1
    x = x | (x >> 1)
    x = x | (x >> 2)
    x = x | (x >> 4)
    x = x | (x >> 8)
    x = x | (x >> 16)
    return x + 1
end

function sqr2(x::Real, y::Real)
    isnan(x) && return x
    isnan(y) && return y
    # general case
    xabs = abs(x)
    yabs = abs(y)
    w = max(xabs, yabs)
    z = min(xabs, yabs)
    iszero(z) && return abs2(w)
    return abs2(w) * (1.0 + abs2(z / w))
end

function apy2(x::Real, y::Real)
    isnan(x) && return x
    isnan(y) && return y
    # general case
    xabs = abs(x)
    yabs = abs(y)
    w = max(xabs, yabs)
    z = min(xabs, yabs)
    iszero(z) && return w
    return w * sqrt(1.0 + abs2(z / w))
end

# (xn + i * yn) / (xd + i * yd)
function cdiv(xn::Real, yn::Real, xd::Real, yd::Real)
    m = sqr2(xd, yd)
    return (xn * xd + yn * yd) / m, (yn * xd - xn * yd) / m
end

function polyinterp(x::Real, xv::VecI, yv::VecI, bv::VecB, n::Int)
    one2n = eachindex(1:n)
    @inbounds for i in one2n
        x == xv[i] && return yv[i]
    end
    Δx = Inf
    yp = 0.0
    @inbounds for i in one2n
        δx = abs(x - xv[i])
        δx < Δx && (Δx = δx; yp = yv[i])
    end
    @simd for i in one2n
        @inbounds bv[i] = yv[i] - yp
    end
    @inbounds for k in 1:n-1, i in 1:n-k
        bv[i] += (bv[i] - bv[i+1]) * (x - xv[i]) / (xv[i] - xv[i+k])
    end
    @inbounds return yp + bv[1]
end

function polyinterp!(ydes::VecO{Ty}, xdes::VecO{Tx}, ysrc::VecI, xsrc::VecI; order::Int=3) where {Tx<:Real,Ty<:Real}
    lx = 1                     # left index
    rx = sz = order + 1        # right index
    ia = sz >> 1 + 1           # anchor index
    mx = size(xsrc, 1)         # upper bound of window
    xb = Vector{Tx}(undef, sz) # xsrc buffer
    yb = Vector{Ty}(undef, sz) # ysrc buffer
    tb = Vector{Ty}(undef, sz) # interpolation buffer
    @inbounds for i in 1:sz
        xb[i] = xsrc[i]
        yb[i] = ysrc[i]
    end
    @inbounds xanchor = xb[ia]
    @inbounds for i in eachindex(xdes)
        xi = xdes[i]
        #### moving interpolation window
        while xi > xanchor && rx < mx
            lx += 1
            rx += 1
            ix  = 1
            jx  = 2
            while ix < sz
                @inbounds xb[ix] = xb[jx]
                @inbounds yb[ix] = yb[jx]
                ix += 1
                jx += 1
            end
            @inbounds xb[ix]  = xsrc[rx]
            @inbounds yb[ix]  = ysrc[rx]
            @inbounds xanchor = xb[ia]
        end
        ydes[i] = polyinterp(xi, xb, yb, tb, sz)
    end
    return nothing
end

export polyinterp!

include("./TimeDelay.jl")
include("./DifferentialData.jl")
include("./TimeDomainData.jl")

end # module


text:
import GeometricMultigrid: δ,inside

# set up examples
function static(;n=32,N=2,T::Type=Float32,m=rand(T,N))
    dims = ntuple(i->n+2,N)
    L,x = zeros(T,dims...,N), zeros(T,dims)
    m ./= √sum(abs2,m)
    for I ∈ inside(x)
        x[I] = sum( @. m*I.I)
        for i ∈ 1:N
            I.I[i] !=2 && (L[I,i] = 1)
        end
    end
    A,x = Poisson(L),FieldVector(x)
    A,A*x
end

function dipole(;n=32,N=2,T::Type=Float32,x0=rand(T,N),m=rand(T,N),s=rand(),Δt=1)
    dims = ntuple(i->n+2,N)
    L,b = zeros(T,dims...,N), zeros(T,dims)
    x0 .*= n/3; x0 .+= (n/3); s *= n/4; m ./= √sum(abs2,m)
    f(x) = sum(@. m*(x-x0))*exp(-s^2*sum(abs2,x.-x0))
    for I ∈ inside(b)
        b[I] = f(I.I)
        for i ∈ 1:N
            I.I[i] !=2 && (L[I,i] = Δt)
        end
    end
    b .-= sum(b)/n^N
    Poisson(L),FieldVector(b)
end

@fastmath kern₀(d) = 0.5+0.5d+0.5sin(π*d)/π
μ₀(d) = kern₀(clamp(d,-1,1))
function sphere(;n=32,N=2,T::Type=Float32,x0=rand(T,N),m=rand(T,N),s=rand(),Δt=1)
    dims = ntuple(i->n+2,N)
    L,b = zeros(T,dims...,N), zeros(T,dims)
    x0 .*= n/3; x0 .+= (n/3); s *= n/4; m ./= √sum(abs2,m)
    d(x) = √sum(abs2,x.-x0)-s
    for I ∈ inside(b), i ∈ 1:N
        μp = μ₀(d(I.I .+ 0.5 .* δ(i,N).I))
        μm = μ₀(d(I.I .- 0.5 .* δ(i,N).I))
        I.I[i] !=2 && (L[I,i] = Δt*μm)
        b[I] += Δt*m[i]*(μp-μm)
    end
    Poisson(L),FieldVector(b)
end

using DataStructures
function create_synthetic(;len=100,kw...)
    OrderedDict(
        "2D-static"=>[static(N=2;kw...) for i ∈ 1:len],
        "3D-static"=>[static(N=3;kw...) for i ∈ 1:len],
        "2D-dipole"=>[dipole(N=2;kw...) for i ∈ 1:len],
        "3D-dipole"=>[dipole(N=3;kw...) for i ∈ 1:len],
        "2D-sphere"=>[sphere(N=2;kw...) for i ∈ 1:len],
        "3D-sphere"=>[sphere(N=3;kw...) for i ∈ 1:len]
    )
end


text:
import Base: +

#
#       NORTH
# WEST    +    EAST
#       SOUTH
#

const NORTH = :north
const EAST = :east
const SOUTH = :south
const WEST = :west

function rotate(position, direction)
    # Rotate:   Left  -1  <- ->  +1  Right
    compass = (:north, :east, :south, :west)

    idx = findfirst(x -> x == position, compass)

    if direction < 0 && idx == 1
        return compass[end]
    elseif direction > 0 && idx == length(compass)
        return compass[1]
    else
        return compass[idx + direction]
    end
end

struct Point{T<:Integer}
    x::T
    y::T
end

+(a::Point, b::Point) = Point(a.x + b.x, a.y + b.y)

mutable struct Robot
    position::Point
    direction::Symbol

    Robot(position::Tuple{Integer, Integer}, direction::Symbol) = new(Point(position...), direction)
end

position(r::Robot) = r.position
heading(r::Robot) = r.direction

function turn_right!(r::Robot)
    r.direction = rotate(r.direction, +1)
    return r
end

function turn_left!(r::Robot)
    r.direction = rotate(r.direction, -1)
    return r
end

function advance!(r::Robot)
    movements = Dict(
        :north => Point(0, 1),
        :east => Point(1, 0),
        :south => Point(0, -1),
        :west => Point(-1, 0)
    )
    r.position += movements[r.direction]
    return r
end

function move!(r::Robot, path::String)
    for m in path
        if m == 'L'
            turn_left!(r)
        elseif m == 'R'
            turn_right!(r)
        elseif m == 'A'
            advance!(r)
        else
            error("Wrong input")
        end
    end
    return r
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
using Test

@time @testset "Periodic Operators" begin include("periodic_operators_test.jl") end
@time @testset "Non-Periodic Operators" begin include("SBP_operators_test.jl") end
@time @testset "Dissipation Operators" begin include("dissipation_operators_test.jl") end
@time @testset "Variable Coefficient Operators" begin include("var_coef_operators_test.jl") end
@time @testset "Banded Matrices" begin include("banded_matrices_test.jl") end
@time @testset "Fourier Operators" begin include("fourier_operators_test.jl") end
@time @testset "Legendre Operators" begin include("legendre_operators_test.jl") end
@time @testset "Sum of Operators" begin include("sum_of_operators_test.jl") end
@time @testset "Conservation Laws" begin
    include("conservation_laws/burgers_test.jl")
    include("conservation_laws/cubic_test.jl")
    include("conservation_laws/variable_linear_advection_test.jl")
end
@time @testset "Second Order Equations" begin
    include("second_order_eqs/wave_eq_test.jl")
end


text:
export iauTcgtt
"""
Time scale transformation:  Geocentric Coordinate Time, TCG, to
Terrestrial Time, TT.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  canonical.

Given:
   tcg1,tcg2  double    TCG as a 2-part Julian Date

Returned:
   tt1,tt2    double    TT as a 2-part Julian Date

Returned (function value):
               int       status:  0 = OK

Note:

   tcg1+tcg2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tcg1 is the Julian
   Day Number and tcg22 is the fraction of a day.  The returned
   tt1,tt2 follow suit.

References:

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),.
   IERS Technical Note No. 32, BKG (2004)

   IAU 2000 Resolution B1.9

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauTcgtt(double tcg1, double tcg2, double *tt1, double *tt2)

function iauTcgtt(tcg1::Real, tcg2::Real)
   ref_tt1 = Ref{Float64}(0.0)
   ref_tt2 = Ref{Float64}(0.0)

   status = ccall((:iauTcgtt, libsofa_c), Cint,
                   (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
                   convert(Float64, tcg1), convert(Float64, tcg2), ref_tt1, ref_tt2)

   return status, ref_tt1[], ref_tt2[]
end


text:
module ContinuumAbsorption

using ..Korg: ionization_energies, @species_str, _data_dir # not sure that this is the best idea
using ..Korg: Interval, closed_interval, contained, contained_slice, λ_to_ν_bound
include("../constants.jl") # I'm not thrilled to duplicate this, but I think it's probably alright

# define helper functions
include("bounds_checking.jl")
include("hydrogenic_bf_ff.jl")
include("stancil_tables.jl")

export H_I_bf, H_I_ff, Hminus_bf, Hminus_ff, H2plus_bf_and_ff
include("absorption_H.jl")

export He_II_bf, He_II_ff, Heminus_ff
include("absorption_He.jl")

include("scattering.jl")

# the following are only imported for computing experimental metal bf continuum opacities
using ..Korg: partition_funcs, Species, Formula, ismolecule, get_roman_numeral
export absorption_coef_bf_TOPBase
include("absorption_metal.jl")

export total_continuum_absorption

"""
    total_continuum_absorption(νs, T, nₑ, number_densities, partition_funcs; error_oobounds)

The total continuum linear absoprtion coefficient, α, at many frequencies, ν.

# Arguments

- `νs` are frequencies in Hz
- `T` is temperature in K
- `nₑ` is the electron number density in cm^-3
- `number_densities` is a `Dict` mapping each `Species` to its number density
- `partition_funcs` is a `Dict` mapping each `Species` to its partition function (e.g.
  `Korg.partition_funcs`)
- `error_oobounds::Bool` specifies the behavior of most continnum absorption sources when passed
   frequencies or temperature values that are out of bounds for their implementation. When `false`
   (the default), those absorption sources are ignored at those values. Otherwise, an error is
   thrown.

!!! note
    For efficiency reasons, `νs` must be sorted. While this function technically supports any
    sorted `AbstractVector`, it is most effient when passed an  `AbstractRange`.
"""
function total_continuum_absorption(νs::AbstractVector{F}, T::F, nₑ::F, number_densities::Dict,
                                    partition_funcs::Dict; error_oobounds = false) where F <: Real
    α = zeros(promote_type(F, valtype(number_densities)), length(νs))

    kwargs = Dict(:out_α => α, :error_oobounds => error_oobounds)

    # Hydrogen continuum absorption
    nH_I = number_densities[species"H_I"]
    nH_I_div_U = nH_I / partition_funcs[species"H_I"](T)
    H_I_bf(νs, T, nH_I_div_U; kwargs...)
    H_I_ff(νs, T, number_densities[species"H_II"], nₑ; kwargs...)
    Hminus_bf(νs, T, nH_I_div_U, nₑ; kwargs...)
    Hminus_ff(νs, T, nH_I_div_U, nₑ; kwargs...)
    H2plus_bf_and_ff(νs, T, nH_I, number_densities[species"H_II"]; kwargs...)

    # He continuum opacities
    He_II_bf(νs, T, number_densities[species"H_II"] / partition_funcs[species"H_II"](T); kwargs...)
    He_II_ff(νs, T, number_densities[species"He_III"], nₑ; kwargs...)
    Heminus_ff(νs, T, number_densities[species"He_I"] / partition_funcs[species"He_I"](T), nₑ;
               kwargs...)

    # scattering
    α .+= electron_scattering(nₑ)
    α .+= rayleigh(νs, nH_I, number_densities[species"He_I"], number_densities[species"H2"])

    α
end

end