text:
module M3M6AppliedComplexAnalysis
using Weave

weave("src/Lecture1.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture2.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture3.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture4.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture5.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture6.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture7.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture8.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture9.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture10.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture11.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture12.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture14.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture15.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture16.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture17.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture18.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture19.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture20.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture21.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture22.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture23.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture24.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture25.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture26.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture27.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Lecture28.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")

weave("src/Solutions1.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Solutions2.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Solutions3.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
weave("src/Solutions4.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")

end # module


text:
# %% Delay time impact
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))
using GLMakie, DynamicalSystems, InteractiveDynamics

figure = Figure(resolution = (1000, 800))
display(figure)
ax = figure[1, :] = Axis3(figure)
sll = labelslider!(figure, "τ =", 1:100)
figure[2, :] = sll.layout
τ = sll.slider.value

ds = Systems.lorenz()
x = trajectory(ds, 100; Ttr = 100)[:, 1]

R = embed(x, 3, τ[])
Robs = Observable(R.data)

js = (206, 946)

Pobs = Observable([Point3f0(Robs[][j]) for j in js])

lines!(ax, Robs; color = COLORS[1])
scatter!(ax, Pobs; color = :red, markersize = 1000)
ax.azimuth = 5.555530633326979
ax.xlabel = "x(t)"
ax.ylabel = "x(t+τ)"
ax.zlabel = "x(t+2τ)"

on(τ) do t
    Robs[] = embed(x, 3, t).data
    Pobs[] = [Point3f0(Robs[][j]) for j in js]
end

# record_interaction(figure, string(@__FILE__)[1:end-2]*"mp4"; total_time = 15)

ext: jl

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# use the DX reverse mode with:
# * an infinite far plane
# * a float z-buffer
# * glClearDepth(0.0)
# * glDepthFunc(GL_GREATER)
# * either:
#     * glClipControl(GL_LOWER_LEFT, GL_ZERO_TO_ONE) (DX style mapping of the Z values)
#   or:
#     * glDepthRangedNV(-1.0, 1.0)
#
# this way near will get mapped to 1.0 and infinite to 0.0
@enum ProjectionMode begin
    IsometricProjection = 0
    PerspectiveProjectionOpenGL    # maps to -1..1
    PerspectiveProjectionDXReverse # maps to  1..0, ideal for Vulkan too
end

@enum StereoMode begin
    Mono = 0
    ParallelShift
    OffAxis
    ToeIn
end

mutable struct GenericCamera <: CameraBase
    mPosition::Vec3f0
    mRotationMatrix::Mat3f0
    mLookAtDistance::Float32
    mAspectRatio::Float32
    mStereoMode::StereoMode
    mNearClippingPlane::Float32
    mFarClippingPlane::Float32
    mProjectionMode::ProjectionMode
    mHorizontalFieldOfView::Float32

    function GenericCamera()
        this = new(zero(Vec3f0), Mat3f0(1I), 500.0, 4.0 / 3.0, Mono, 0.1, Inf, PerspectiveProjectionDXReverse, 75.0)
        setRotationMatrix(this, Mat3f0(1I))
        return this
    end
end

@class GenericCamera

function setRotationMatrix(this::GenericCamera, _matrix::Mat3f0)
    this.mRotationMatrix = _matrix
    @assert isOrthonormalMatrix(this.mRotationMatrix) "Rotation matrix must be orthonormal!"
end

function getRotationMatrix4(this::GenericCamera)::Mat4f0
    return make_mat4(this.mRotationMatrix)
end

function getTranslationMatrix4(this::GenericCamera)::Mat4f0
    trans = make_mat4(Mat3f0(1I), -this.mPosition.x, -this.mPosition.y, -this.mPosition.z, 1f0)
    return trans
end

function setTarget(this::GenericCamera, _target::Vec3f0)
    setTarget(this, _target, getUpDirection(this))
end

function getUpDirection(this::GenericCamera)::Vec3f0
    up = Vec3f0(this.mRotationMatrix[1, 2], this.mRotationMatrix[2, 2], this.mRotationMatrix[3, 2])
    @assert (distance(col_multiple(getInverseRotationMatrix3(this), Vec3f0(0.0, 1.0, 0.0)), up) < .01)
    return up
end

function getInverseRotationMatrix3(this::GenericCamera)::Mat3f0
    return transpose(this.mRotationMatrix)
end

function setTarget(this::GenericCamera, _target::Vec3f0, _up::Vec3f0)
    forwardVector::Vec3f0 = _target - this.mPosition
    this.mLookAtDistance = norm(forwardVector)
    if (this.mLookAtDistance < .0001) # in case target == mPosition
        this.mLookAtDistance = .0001
        forwardVector = Vec3f0(this.mLookAtDistance, 0, 0)
    end

    forwardVector = forwardVector / this.mLookAtDistance # normalize
    rightVector::Vec3f0 = normalize(cross(forwardVector, _up))
    upVector::Vec3f0 = cross(rightVector, forwardVector)

    rotMatrix = Mat3f0([rightVector[1] upVector[1] -forwardVector[1];
                    rightVector[2] upVector[2] -forwardVector[2];
                    rightVector[3] upVector[3] -forwardVector[3]])

    setRotationMatrix(this, rotMatrix)
end

function getTarget(this::GenericCamera)::Vec3f0
    return this.mPosition + getForwardDirection(this) * getLookAtDistance(this)
end

function getForwardDirection(this::GenericCamera)::Vec3f0
    forward = Vec3f0(-this.mRotationMatrix[1, 3], -this.mRotationMatrix[2, 3], -this.mRotationMatrix[3, 3])
    inv = getInverseRotationMatrix3(this)
    @assert (distance(col_multiple(inv, Vec3f0(0, 0, -1)), forward) < .01)
    return forward
end

function getLookAtDistance(this::GenericCamera)::Float32
    return this.mLookAtDistance
end

function setAspectRatio(this::GenericCamera, _aspectRatio::Float32)
    this.mAspectRatio = _aspectRatio
end

function setPosition(this::GenericCamera, _position::Vec3f0)
    this.mPosition = _position
end

function getPosition(this::GenericCamera)::Vec3f0
    return this.mPosition
end

function rotateAroundTarget_GlobalAxes(this::GenericCamera, _x::Float32, _y::Float32, _z::Float32)
    # move camera so, that the target is the center, then rotate around the
    # global coordinate system
    rotateAroundTaget_helper(this, _x, _y, _z, Mat3f0(1I))
end

function rotateAroundTaget_helper(this::GenericCamera, _x::Float32, _y::Float32, _z::Float32, _rotationAxes::Mat3f0)
    T = make_vec4(getTarget(this), 1.0)
    P = make_vec4(getPosition(this), 1.0)

    tempPos = P - T
    newRotation = rotate(RotMatrix(Mat3f0(1I)), _x, getRow(_rotationAxes, 1))
    newRotation = rotate(newRotation, _y, getRow(_rotationAxes, 2))
    newRotation = rotate(newRotation, _z, getRow(_rotationAxes, 3))
    newRotation = make_mat4(Mat3f0(newRotation))
    tempPos = col_multiple(newRotation, tempPos)

    P = tempPos + T # new position
    N = make_vec4(getUpDirection(this), 1.0)
    N = col_multiple(newRotation, N)

    setLookAtMatrix(this, make_vec3(P), make_vec3(T), make_vec3(N))
end

function setLookAtMatrix(this::GenericCamera, _position::Vec3f0, _target::Vec3f0, _up::Vec3f0)
    setPosition(this, _position)
    setTarget(this, _target, _up)
end

function getViewMatrix(this::GenericCamera)::Mat4f0
    if (this.mStereoMode == Mono)
        return getMonoViewMatrix(this)
    else
        # all kinds of stereo
        eyeIsLeftEye = (getEye(this) == EyeLeft)
        return getStereoViewMatrix(this, eyeIsLeftEye, this.mStereoMode)
    end
end

function getMonoViewMatrix(this::GenericCamera)::Mat4f0
    m = this.mRotationMatrix
    m41 = -(m[1, 1] * this.mPosition.x + m[2, 1] * this.mPosition.y +
                m[3, 1] * this.mPosition.z)
    m42 = -(m[1, 2] * this.mPosition.x + m[2, 2] * this.mPosition.y +
                m[3, 2] * this.mPosition.z)
    m43 = -(m[1, 3] * this.mPosition.x + m[2, 3] * this.mPosition.y +
                m[3, 3] * this.mPosition.z)
    m = make_mat4(m, m41, m42, m43, 1f0)
    @assert (isApproxEqual(col_multiple(getRotationMatrix4(this), getTranslationMatrix4(this)), m))
    return m
end

function getStereoViewMatrix(this::GenericCamera, _leftEye::Bool, _stereoMode::StereoMode)::Mat4f0
    # The view matrix is independent of the projection mode (isometric or
    # perspective)
    # so only the stereo mode has to be checked.
    @assert _stereoMode != Mono "mono is not a stereo mode!"

    cameraPositionShiftValue::Float32 = this.mInterpupillaryDistance * 0.5 # shift to the right
    if (_leftEye)
        cameraPositionShiftValue *= -1.0 # if left eye shift to the left
    end

    if ((_stereoMode == ParallelShift) || (_stereoMode == OffAxis))
        #
        # parallel shift and off-axis have the same view matrices:
        # just shift the camera position to the left/right by half the
        # eye-distance
        #

        # ACGL::Utils::debug(this) << "WARNING: getStereoViewMatrix(this) is not tested
        # yet" << std::endl // remove after
        # testing

        inverseRotation = getInverseRotationMatrix3(this)
        eyePosition::Vec3f0 = this.mPosition + col_multiple(inverseRotation, Vec3f0(cameraPositionShiftValue, 0.0, 0.0))

        m = make_mat4(this.mRotationMatrix)
        m[4, 1] = -(m[1, 1] * eyePosition.x + m[2, 1] * eyePosition.y +
                    m[3, 1] * eyePosition.z)
        m[4, 2] = -(m[1, 2] * eyePosition.x + m[2, 2] * eyePosition.y +
                    m[3, 2] * eyePosition.z)
        m[4, 3] = -(m[1, 3] * eyePosition.x + m[2, 3] * eyePosition.y +
                    m[3, 3] * eyePosition.z)
        return m
    end

    # else it has to be toe-in:
    @assert (_stereoMode == ToeIn) "unsupported stereo mode!"
    #
    # Toe-in: shift the camera position to the left/right by half the
    # eye-distance and
    #         rotate a bit inwards so that the two cameras focus the same point
    #         at the look-at distance (focal point)

    @assert false "getStereoViewMatrix(this) for TOE_IN is not implemented yet!"
    return Mat4f0(1I)
end

function getProjectionMatrix(this::GenericCamera)::Mat4f0
    if (this.mStereoMode == Mono)
        return getMonoProjectionMatrix(this)
    else
        # all kinds of stereo
        eyeIsLeftEye = getEye(this) == EyeLeft
        return getStereoProjectionMatrix(this, eyeIsLeftEye, this.mStereoMode)
    end
end

function getProjectionMode(this::GenericCamera)::ProjectionMode
    return this.mProjectionMode
end

function getVerticalFieldOfView(this::GenericCamera)::Float32
    rad2deg(atan(tan(deg2rad(0.5 * this.mHorizontalFieldOfView)) / this.mAspectRatio) * 2.0)
end

function getAspectRatio(this::GenericCamera)::Float32
    return this.mAspectRatio
end

function getMonoProjectionMatrix(this::GenericCamera)::Mat4f0
    projectionMatrix = Mat4f0(1I) # identity matrix

    if (getProjectionMode(this) == IsometricProjection)
        # we don't set the left/right/top/bottom values explicitly, so we want
        # that
        # all object at our focal distance appear the same in perspective and
        # isometric view
        right::Float32 = tan(deg2rad(getHorizontalFieldOfView(this) * 0.5)) * this.mLookAtDistance
        left::Float32 = -right
        top::Float32 = tan(deg2rad(getVerticalFieldOfView(this) * 0.5f)) * this.mLookAtDistance
        bottom::Float32 = -top

        # we do the same here as a glOrtho call would do, but with flipped y
        projectionMatrix[1, 1] = 2.0 / (right - left)
        projectionMatrix[2, 2] = -2.0 / (top - bottom)
        projectionMatrix[3, 3] = -2.0 / (this.mFarClippingPlane - mNearClippingPlane)
        projectionMatrix[1, 4] = -(right + left) / (right - left)
        projectionMatrix[2, 4] = -(top + bottom) / (top - bottom)
        projectionMatrix[3, 4] = -(this.mFarClippingPlane + mNearClippingPlane) /
                                  (this.mFarClippingPlane - mNearClippingPlane)
        projectionMatrix[4, 4] = 1.0
    elseif (this.mProjectionMode == PerspectiveProjectionOpenGL)
        if (std::isinf(this.mFarClippingPlane))
            e = 1.01 / tan(deg2rad(getVerticalFieldOfView(this) * 0.5))
            a = getAspectRatio(this)

            # infinite Perspective matrix reversed mapping to 1..-1
            projectionMatrix = Mat4f0([e / a 0.0 0.0  0.0;
                                     0.0   -e  0.0  0.0;
                                     0.0   0.0 -1.0 -1.0;
                                     0.0   0.0 -2.0 * this.mNearClippingPlane 0.0])
        else
            projectionMatrix = glm::perspective(
                deg2rad(getHorizontalFieldOfView(this)),
                getAspectRatio(this), this.mNearClippingPlane,
                this.mFarClippingPlane)
            projectionMatrix[1, 2] *= -1.0
            projectionMatrix[2, 2] *= -1.0
            projectionMatrix[3, 2] *= -1.0
            projectionMatrix[4, 2] *= -1.0
        end
    elseif (this.mProjectionMode == PerspectiveProjectionDXReverse)
        if (isinf(this.mFarClippingPlane))
            e = 1.0 / tan(deg2rad(getVerticalFieldOfView(this) * 0.5))
            a = getAspectRatio(this)

            # infinite Perspective matrix reversed mapping to 1..0, but flip y
            projectionMatrix = Mat4f0([e / a 0.0 0.0  0.0;  #
                                     0.0   -e  0.0  0.0;  #
                                     0.0   0.0 0.0 -1.0; #
                                     0.0   0.0 this.mNearClippingPlane 0.0])
        else
            @assert false "unsupported projection mode"
        end
    else
        @assert false "unsupported projection mode"
    end

    return projectionMatrix
end

function getStereoProjectionMatrix(this::GenericCamera, _leftEye::Bool, _stereoMode::StereoMode)::Mat4f0
    @assert (_stereoMode != Mono) "mono is not a stereo mode!"

    if (getProjectionMode(this) == IsometricProjection)
        # very unusual, prepare for headaches!
        return getMonoProjectionMatrix(this)
    end

    if ((_stereoMode == ParallelShift) || (_stereoMode == ToeIn))
        # the view matrix changes but the projection matrix stays the same
        return getMonoProjectionMatrix(this)
    end

    # so off-axis it is!
    @assert (_stereoMode == OffAxis) "unknown projection mode!"
    #
    # In this mode the camera positions (view matrix) is shifted to the
    # left/right still looking
    # straight ahead. The projection is also looking ahead but the projection
    # center is
    # off (hence off-axis).
    # There is one plane in front of the cameras where the view-frusta match.
    # This should be the distance to the physical screen from the users
    # position.

    @assert false "getStereoViewMatrix(this) is not implemented for OFF_AXIS yet!"
    return Mat4f0(1I)
end


text:
"""
Sw(args...)

Calculate the water saturation

`` { Sw }^{ n }=\\frac { a\\quad { Rw } }{ { \\Phi  }^{ m }\\quad Rt } ``

## Input Example

```julia
SwArchie(a,m,n,Rw,phi,rt)
```
"""
function SwFunction(a,m,n,Rw,phi,rt;Rsh=4.0,Vsh=0,α=0.3,Method="Archie")

    if Method == "Archie"
        Sw=((a.*Rw)./(rt.*phi.^m)).^(1/n)
    elseif Method == :Smdx #https://www.spec2000.net/14-sws.htm

            C=((1 .-Vsh).*a.*Rw)./(phi.^m)
            D=C.*Vsh./(2 .* Rsh)
            E=C./rt
            Sw=((D.^2 .+ E).^0.5 .- D).^(2 ./n)
    elseif Method == "Indo"
       A=(1 ./rt).^0.5
       B=(Vsh.^(1 .-(Vsh./2)))/(Rsh.^0.5)
       C=((phi.^m)./(a*Rw)).^0.5
       Sw=(A./(B+C)).^(2 ./n)
   elseif Method == "Fertl"
       A=phi.^(-m./2)
       B=(a.*Rw)./rt
       C=((α.*Vsh)./2).^2
       D=(α.*Vsh)./2
       Sw=A.*(((B.+C).^0.5)-D)
   end

    return Sw
end


text:
struct HalfEdge
    src::Int
    dst::Int
end
src(he::HalfEdge) = he.src
dst(he::HalfEdge) = he.dst

# NOTE:
# isolated vertices can only appear without any other connection.
# thus we can just use a vertice -> face map for isolated vertices.
mutable struct PlanarMultigraph
    v2he::Dict{Int, Int}    # v_id -> he_id
    half_edges::Dict{Int, HalfEdge} # he_id -> he

    f2he::Dict{Int, Int}    # f_id -> he_id
    he2f::Dict{Int, Int}    # he_id -> f_id

    next::Dict{Int, Int}    # he_id -> he_id
    twin::Dict{Int, Int}    # he_id -> he_id

    vs_isolated::Dict{Int, Int} # v_id -> f_id

    v_max::Int
    he_max::Int
    f_max::Int
end

Base.copy(g::PlanarMultigraph) = PlanarMultigraph(
    copy(g.v2he), copy(g.half_edges), copy(g.f2he),
    copy(g.he2f), copy(g.next), copy(g.twin),
    copy(g.vs_isolated), g.v_max, g.he_max, g.f_max
)

vertices(g::PlanarMultigraph) = vcat(collect(keys(g.v2he)), collect(keys(g.vs_isolated)))
isolated_vertices(g::PlanarMultigraph) = collect(keys(g.vs_isolated))
is_isolated(g::PlanarMultigraph, v::Integer) = haskey(g.vs_isolated, v)

faces(g::PlanarMultigraph) = sort!(collect(keys(g.f2he)))
half_edges(g::PlanarMultigraph) = sort!(collect(keys(g.half_edges)))

src(g::PlanarMultigraph, he_id::Integer) = src(g.half_edges[he_id])
dst(g::PlanarMultigraph, he_id::Integer) = dst(g.half_edges[he_id])
half_edge(g::PlanarMultigraph, he_id::Integer) = g.half_edges[he_id]
face(g::PlanarMultigraph, he_id::Integer) = g.he2f[he_id]

α(g::PlanarMultigraph, he::Integer) = g.next[he]
next(g::PlanarMultigraph, he::Integer) = α(g, he)
function prev(g::PlanarMultigraph, he::Integer)
    current_he = next(g, he)
    next_he = next(g, current_he)
    while next_he != he
        current_he = next_he
        next_he = next(g, next_he)
    end
    return current_he
end

# NOTE: rename ϕ and σ
ϕ(g::PlanarMultigraph, he::Integer) = g.twin[he]
twin(g::PlanarMultigraph, he::Integer) = ϕ(g, he)

σ(g::PlanarMultigraph, he::Integer) = twin(g, prev(g, he))
σ_inv(g::PlanarMultigraph, he::Integer) = next(g, twin(g, he))

nv(g::PlanarMultigraph) = length(g.v2he) + length(g.vs_isolated)
nf(g::PlanarMultigraph) = length(g.f2he)
nhe(g::PlanarMultigraph) = length(g.half_edges)
ne(g::PlanarMultigraph) = nhe(g) ÷ 2

function out_half_edge(g::PlanarMultigraph, v::Integer)
    is_isolated(g, v) && return 0
    return g.v2he[v]
end
surrounding_half_edge(g::PlanarMultigraph, f::Integer) = g.f2he[f]

function trace_orbit(f::Function, a::T; rev::Bool = false) where T
    next = f(a)
    perm = T[a]
    while next != a
        if rev
            pushfirst!(perm, next)
        else
            push!(perm, next)
        end
        next = f(next)
    end
    return perm
end
function trace_face(g::PlanarMultigraph, f::Integer; safe_trace = false)
    !safe_trace && return trace_orbit(h -> g.next[h], surrounding_half_edge(g, f))
    hes_f = Int[]
    for (he, f_he) in g.he2f
        f_he == f && push!(hes_f, he)
    end
    return hes_f
end
function trace_vertex(g::PlanarMultigraph, v::Integer)
    is_isolated(g, v) && return Int[]
    return trace_orbit(h -> σ_inv(g, h), out_half_edge(g, v); rev = true)
end
neighbors(g::PlanarMultigraph, v::Integer) = [dst(g, he) for he in trace_vertex(g, v)]

is_boundary(g::PlanarMultigraph, he_id::Integer) = (face(g, he_id) == 0)

function rem_vertex!(g::PlanarMultigraph, v::Integer; update::Bool = true)
    for he in trace_vertex(g, v)
        rem_edge!(g, he; update = update)
    end
    delete!(g.v2he, v)
    delete!(g.vs_isolated, v)
    return g
end

function rem_edge!(g::PlanarMultigraph, he_id::Integer; update::Bool = true)
    # make sure the face of he_id is an inner face
    if is_boundary(g, he_id)
        he_id = twin(g, he_id)
    end

    # handle self-loop
    if next(g, twin(g, he_id)) == twin(g, he_id)
        if next(g, he_id) == he_id  # isolated self-loop
            v_loop = src(g, he_id)
            twin_id = twin(g, he_id)
            he_in = (surrounding_half_edge(g, face(g, he_id)) == he_id) ? he_id : twin_id
            he_out = twin(g, he_in)
            f_in = face(g, he_in)
            f_out = face(g, he_out)
            for (v_iso, f_iso) in g.vs_isolated
                f_iso == f_in && (g.vs_isolated[v_iso] = f_out)
            end
            hes_f_he = trace_face(g, f_in; safe_trace = true)
            rem_face!(g, f_in)
            for he in hes_f_he
                g.he2f[he] = f_out
            end
            for he_rm in (he_in, he_out)
                delete!(g.he2f, he_rm)
                delete!(g.half_edges, he_rm)
                delete!(g.next, he_rm)
                delete!(g.twin, he_rm)
            end
            if update
                delete!(g.v2he, v_loop)
                g.vs_isolated[v_loop] = f_out
            end
            return g
        end
        he_id = twin(g, he_id)
        is_boundary(g, he_id) && (he_id = twin(g, he_id))
    end

    twin_id = twin(g, he_id)

    if update
        he_next = next(g, he_id)
        he_prev = prev(g, he_id)
        twin_next = next(g, twin_id)
        twin_prev = prev(g, twin_id)

        face_he = face(g, he_id)
        face_twin = face(g, twin_id)

        # remove a inner face
        if face_he != face_twin
            for (v_iso, f_iso) in g.vs_isolated
                f_iso == face_he && (g.vs_isolated[v_iso] = face_twin)
            end
            hes_f_he = trace_face(g, face_he; safe_trace = true)
            rem_face!(g, face_he)
            for he in hes_f_he
                g.he2f[he] = face_twin
            end
        end

        # update f2he
        if surrounding_half_edge(g, face_twin) in (he_id, twin_id)
            new_he = nothing
            for nhe in (he_next, he_prev, twin_next, twin_prev)
                if !(nhe in (he_id, twin_id))
                    g.f2he[face_twin] = nhe
                    new_he = nhe
                    break
                end
            end
            new_he === nothing && error("surrounding half edge not founded")
        end

        # update v2he
        (out_half_edge(g, src(g, he_id)) == he_id) && (g.v2he[src(g, he_id)] = twin(g, he_prev))
        if out_half_edge(g, src(g, twin_id)) in (he_id, twin_id)
            if twin_id == twin_next
                g.v2he[src(g, twin_id)] = he_next
            else
                g.v2he[src(g, twin_id)] = twin(g, twin_prev)
            end
        end

        if he_next == he_id # he_id is the inner half edge of a self-loop
            g.next[twin_prev] = twin_next
        elseif twin_next == twin_id
            g.next[he_prev] = he_next
        else
            g.next[he_prev] = twin_next
            g.next[twin_prev] = he_next
        end

        if he_next == twin_id
            g.vs_isolated[src(g, he_next)] = face(g, he_next)
            delete!(g.v2he, src(g, he_next))
        end
        if twin_next == he_id
            g.vs_isolated[src(g, twin_next)] = face(g, twin_next)
            delete!(g.v2he, src(g, twin_next))
        end
    end
    delete!(g.next, he_id)
    delete!(g.next, twin_id)
    delete!(g.half_edges, he_id)
    delete!(g.half_edges, twin_id)
    delete!(g.twin, he_id)
    delete!(g.twin, twin_id)
    delete!(g.he2f, he_id)
    delete!(g.he2f, twin_id)

    return g
end

function rem_face!(g::PlanarMultigraph, f::Integer)
    f == 0 && error("Face 0 is for the boundary. It can not be removed.")
    half_edges_f = trace_face(g, f)
    for he in half_edges_f
        delete!(g.he2f, he)
    end
    delete!(g.f2he, f)
    return g
end

function merge_graph!(A::PlanarMultigraph, B::PlanarMultigraph)
    for (v, he_id) in B.v2he
        A.v2he[v + A.v_max] = he_id + A.he_max
    end
    for (he_id, he) in B.half_edges
        A.half_edges[he_id + A.he_max] = HalfEdge(src(he) + A.v_max, dst(he) + A.v_max)
    end
    for (f, he_id) in B.f2he
        if f != 0
            A.f2he[f + A.f_max] = he_id + A.he_max
        end
    end
    for (he_id, f) in B.he2f
        if f != 0
            A.he2f[he_id + A.he_max] = f + A.f_max
        else
            A.he2f[he_id + A.he_max] = 0
        end
    end
    for (curr, next) in B.next
        A.next[curr + A.he_max] = next + A.he_max
    end
    for (curr, twin) in B.twin
        A.twin[curr + A.he_max] = twin + A.he_max
    end
    for (v, f) in B.vs_isolated
        A.vs_isolated[v + A.v_max] = (f == 0) ? 0 : (f + A.f_max)
    end

    A.v_max += B.v_max
    A.he_max += B.he_max
    A.f_max += B.f_max

    return A
end

function check_faces(g::PlanarMultigraph)
    for f in faces(g)
        hes_f = trace_face(g, f)
        for he in hes_f
            face(g, he) == f || return false
        end
    end
    return true
end
function check_vertices(g::PlanarMultigraph)
    for v in vertices(g)
        hes_v = trace_vertex(g, v)
        for he in hes_v
            src(g, he) == v || return false
        end
    end
    return true
end
function check_combinatorial_maps(g::PlanarMultigraph)
    for he in half_edges(g)
        (he == α(g, ϕ(g, σ(g, he)))) || return false
    end
    return true
end

function update_face!(g::PlanarMultigraph, he_id)
    face_id = face(g, he_id)
    fs_rm = Int[]
    g.f2he[face_id] = he_id
    curr_he = next(g, he_id)
    while curr_he != he_id
        if g.he2f[curr_he] != face_id
            push!(fs_rm, g.he2f[curr_he])
            g.he2f[curr_he] = face_id
        end
        curr_he = next(g, curr_he)
    end
    for (v, f) in g.vs_isolated
        (f in fs_rm) && (g.vs_isolated[v] = face_id)
    end
    return g
end

function add_vertex!(g::PlanarMultigraph, f::Integer)
    haskey(g.f2he, f) || return 0
    g.v_max += 1
    v = g.v_max
    g.vs_isolated[v] = f
    return v
end

function add_edge_isolated_1!(g::PlanarMultigraph, v1::Integer, v2::Integer, f::Integer)
    f == g.vs_isolated[v1] || return (0, 0)
    hes = trace_vertex(g, v2)
    he2_in = 0
    he2_out = 0
    for he in hes
        if face(g, twin(g, he)) == f
            he2_in = twin(g, he)
            he2_out = next(g, he2_in)
            break
        end
    end
    he2_in * he2_out != 0 || return (0, 0)

    g.he_max += 2
    new_he1 = g.he_max - 1
    new_he2 = g.he_max
    g.v2he[v1] = new_he2
    g.twin[new_he1] = new_he2
    g.twin[new_he2] = new_he1
    g.next[he2_in] = new_he1
    g.next[new_he1] = new_he2
    g.next[new_he2] = he2_out
    g.he2f[new_he1] = f
    g.he2f[new_he2] = f
    g.half_edges[new_he1] = HalfEdge(v2, v1)
    g.half_edges[new_he2] = HalfEdge(v1, v2)
    delete!(g.vs_isolated, v1)

    return (new_he1, new_he2)
end
function add_edge_isolated_2!(g::PlanarMultigraph, v1::Integer, v2::Integer, f::Integer)
    f == g.vs_isolated[v1] == g.vs_isolated[v2] || return (0, 0)
    g.he_max += 2
    new_he1 = g.he_max - 1
    new_he2 = g.he_max
    g.twin[new_he1] = new_he2
    g.twin[new_he2] = new_he1
    g.next[new_he1] = new_he2
    g.next[new_he2] = new_he1
    g.v2he[v1] = new_he1
    g.v2he[v2] = new_he2
    g.he2f[new_he1] = f
    g.he2f[new_he2] = f
    g.half_edges[new_he1] = HalfEdge(v1, v2)
    g.half_edges[new_he2] = HalfEdge(v2, v1)
    delete!(g.vs_isolated, v1)
    delete!(g.vs_isolated, v2)

    return (new_he1, new_he2)
end

function add_edge!(g::PlanarMultigraph, v1::Integer, v2::Integer, f::Integer)
    if is_isolated(g, v1)
        if is_isolated(g, v2)
            return add_edge_isolated_2!(g, v1, v2, f)
        else
            return add_edge_isolated_1!(g, v1, v2, f)
        end
    elseif is_isolated(g, v2)
        return add_edge_isolated_1!(g, v2, v1, f)
    end
    hes_f = trace_face(g, f)
    he1_in, he1_out, he2_in, he2_out = (0,0,0,0)
    for he in hes_f
        dst(g, he) == v1 && (he1_in = he; he1_out = next(g, he))
        dst(g, he) == v2 && (he2_in = he; he2_out = next(g, he))
    end
    he1_in * he1_out * he2_in * he2_out != 0 || return (0, 0)
    new_he1 = g.he_max + 1
    new_he2 = g.he_max + 2
    g.he_max += 2
    g.twin[new_he1] = new_he2
    g.twin[new_he2] = new_he1
    g.half_edges[new_he1] = HalfEdge(v1, v2)
    g.half_edges[new_he2] = HalfEdge(v2, v1)
    g.next[he1_in] = new_he1
    g.next[new_he1] = he2_out
    g.next[he2_in] = new_he2
    g.next[new_he2] = he1_out
    g.he2f[new_he1] = f
    g.f2he[f] = new_he1
    g.f_max += 1
    g.he2f[new_he2] = g.f_max
    g.f2he[g.f_max] = new_he2
    update_face!(g, new_he2)
    return (new_he1, new_he2)
end

function contract_edge!(g::PlanarMultigraph, he_id::Integer)
    twin_id = twin(g, he_id)
    he_prev = prev(g, he_id)
    he_next = next(g, he_id)
    twin_prev = prev(g, twin_id)
    twin_next = next(g, twin_id)

    v1 = src(g, he_id)
    v2 = dst(g, he_id)
    if v1 == v2
        rem_edge!(g, he_id; update = true)
        return (v1, v2)
    end
    if length(trace_vertex(g, v1)) == 1
        if length(trace_vertex(g, v2)) > 1
            (v11, v22) = contract_edge!(g, twin_id)
            return (v22, v11)
        else
            # 2 isolated vertices
            error("TODO: fix this")
        end
    end

    # update out half edge of v1
    out_half_edge(g, v1) == he_id && (g.v2he[v1] = twin_next)
    for he in trace_vertex(g, v2)
        v0 = dst(g, he)
        g.half_edges[he] = HalfEdge(v1, v0)
        g.half_edges[twin(g, he)] = HalfEdge(v0, v1)
    end

    if he_next == twin_id
        g.next[he_prev] = twin_next
        g.f2he[face(g, he_id)] = he_prev
    else
        g.next[he_prev] = he_next
        g.next[twin_prev] = twin_next
        g.f2he[face(g, he_id)] = he_prev
        g.f2he[face(g, twin_id)] = twin_prev
    end
    delete!(g.next, he_id)
    delete!(g.next, twin_id)
    delete!(g.half_edges, he_id)
    delete!(g.half_edges, twin_id)
    delete!(g.twin, he_id)
    delete!(g.twin, twin_id)
    delete!(g.he2f, he_id)
    delete!(g.he2f, twin_id)
    delete!(g.v2he, v2)
    # v2 is removed
    return (v1, v2)
end

"""
    split_edge!(g::PlanarMultigraph, he_id)

Split the edge corresponding to `he` into 2 edges.
This is used for creating planar simple graphs from planar multigraphs.
"""
function split_edge!(g::PlanarMultigraph, he_id::Integer)
    he1 = he_id
    he2 = twin(g, he_id)
    next1 = next(g, he1)
    next2 = next(g, he2)
    f1 = face(g, he1)
    f2 = face(g, he2)
    s = src(g, he_id)
    d = dst(g, he_id)
    g.v_max += 1
    v = g.v_max
    g.he_max += 2
    nhe1 = g.he_max - 1
    nhe2 = g.he_max

    g.next[he1] = nhe1
    g.next[nhe1] = next1
    g.next[he2] = nhe2
    g.next[nhe2] = next2
    g.twin[nhe1] = he2
    g.twin[he2] = nhe1
    g.twin[nhe2] = he1
    g.twin[he1] = nhe2

    g.half_edges[he1] = HalfEdge(s, v)
    g.half_edges[nhe2] = HalfEdge(v, s)
    g.half_edges[nhe1] = HalfEdge(v, d)
    g.half_edges[he2] = HalfEdge(d, v)

    g.he2f[nhe1] = f1
    g.he2f[nhe2] = f2

    g.v2he[v] = nhe1

    return v, nhe1, nhe2
end

function simple_connected_planar_graph(g::PlanarMultigraph)
    g = copy(g)
    vs_iso = copy(g.vs_isolated)
    for (v, f) in vs_iso
        v_f = src(g, surrounding_half_edge(g, f))
        add_edge!(g, v, v_f, f)
    end
    if n_conn_comp(g) != 1
        @error "there are more than ONE connected components"
    end

    hes_splitted = Dict{Int, Int}() # old he -> new v
    vs = vertices(g)
    while !isempty(vs)
        v = pop!(vs)
        hes_v = trace_vertex(g, v)
        nbs = Set{Int}()
        for he in hes_v
            d = dst(g, he)
            if d == v   # self-loop
                he_twin = twin(g, he)
                new_v1, he1, _ = split_edge!(g, he)
                push!(vs, v)
                hes_splitted[he] = new_v1
                new_v2, _, _ = split_edge!(g, he1)
                hes_splitted[he_twin] = new_v2
                break
            elseif !(d in nbs)
                push!(nbs, d)
            else
                he_twin = twin(g, he)
                new_v, _, _ = split_edge!(g, he)
                push!(vs, v)
                hes_splitted[he] = new_v
                hes_splitted[he_twin] = new_v
                break
            end
        end
    end
    return g, hes_splitted
end

"""
    normalize(g)

Return a relabeled planar graph.
"""
function normalize(g::PlanarMultigraph)
    f_max = nf(g) - 1
    he_max = nhe(g)
    v_max = nv(g)

    f_map = Dict(zip(faces(g), 0:f_max))
    he_map = Dict(zip(half_edges(g), 1:he_max))
    v_map = Dict(zip(sort!(vertices(g)), 1:v_max))

    v2he = Dict{Int, Int}(v_map[v] => he_map[he] for (v, he) in g.v2he)
    halfedges = Dict{Int, HalfEdge}(
        he_map[he_id] => HalfEdge(v_map[he.src], v_map[he.dst]) for (he_id, he) in g.half_edges
    )

    f2he = Dict{Int, Int}(f_map[f] => he_map[he] for (f, he) in g.f2he)
    he2f = Dict{Int, Int}(he_map[he] => f_map[f] for (he, f) in g.he2f)

    next = Dict{Int, Int}(he_map[cur] => he_map[nxt] for (cur, nxt) in g.next)
    twin = Dict{Int, Int}(he_map[cur] => he_map[twn] for (cur, twn) in g.twin)

    vs_isolated = Dict{Int, Int}(v_map[v] => f_map[f] for (v, f) in g.vs_isolated)

    g_new = PlanarMultigraph(v2he, halfedges, f2he, he2f, next, twin, vs_isolated,
        v_max, he_max, f_max)
    return g_new, v_map, he_map, f_map
end

"""
    n_conn_comp(g::PlanarMultigraph)

Return the number of connected components.
"""
n_conn_comp(g::PlanarMultigraph) = nv(g) - ne(g) + nf(g) - 1

has_vertex(g::PlanarMultigraph, v) = haskey(g.v2he, v) || haskey(g.vs_isolated, v)
has_half_edge(g::PlanarMultigraph, he) = haskey(g.half_edges, he)
has_face(g::PlanarMultigraph, f) = haskey(g.f2he, f)


text:
module PandorasBoxLaserSensor

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/laserSensor" LaserSensor(x::Integer, y::Integer, flag::String="", color::String="White", blockLight::Bool=true, mode::String="Continuous")

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))
const modes = String["Continuous", "HitOnce"]

const placements = Ahorn.PlacementDict(
    "Laser Sensor (Pandora's Box)" => Ahorn.EntityPlacement(
        LaserSensor,
        "point"
    )
)

Ahorn.editingOptions(entity::LaserSensor) = Dict{String, Any}(
    "color" => colors,
    "mode" => modes
)

function getColor(color)
    if haskey(Ahorn.XNAColors.colors, color)
        return Ahorn.XNAColors.colors[color]

    else
        try
            return ((Ahorn.argb32ToRGBATuple(parse(Int, replace(color, "#" => ""), base=16))[1:3] ./ 255)..., 1.0)

        catch

        end
    end

    return (1.0, 1.0, 1.0, 1.0)
end

orbTexture = "objects/pandorasBox/laser/sensor/orb"
metalRingTexture = "objects/pandorasBox/laser/sensor/metal_ring"
lightRingTexture = "objects/pandorasBox/laser/sensor/light_ring"

function Ahorn.selection(entity::LaserSensor)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(metalRingTexture, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::LaserSensor, room::Maple.Room)
    colorRaw = entity.color
    color = getColor(colorRaw)

    Ahorn.drawSprite(ctx, metalRingTexture, 0, 0)
    Ahorn.drawSprite(ctx, lightRingTexture, 0, 0, tint=color)
    Ahorn.drawSprite(ctx, orbTexture, 0, 0)
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
"""
App level functionality -- loading and managing app-wide components like configs, models, initializers, etc.
"""
module App

using Genie.Configuration

if is_dev()
  @eval using Revise
end

using Genie, YAML, Macros, Logger, Inflector, Util
SEARCHLIGHT_ON && eval(:(using SearchLight, Validation))

IS_IN_APP && (isdefined(:config) || const config = Genie.config)

"""
    load_libs(dir = Genie.LIB_PATH) :: Void

Recursively adds subfolders of lib to LOAD_PATH.
"""
function load_libs(dir = Genie.LIB_PATH) :: Void
  lib_dirs = [d::String for d::String in Util.walk_dir(dir, only_dirs = true)]
  ! isempty(lib_dirs) && push!(LOAD_PATH, lib_dirs...)

  # Util.reload_modules([dir, lib_dirs...], current_module())

  nothing
end

"""
    load_resources(dir = Genie.RESOURCES_PATH) :: Void

Recursively adds subfolders of resources to LOAD_PATH.
"""
function load_resources(dir = Genie.RESOURCES_PATH) :: Void
  res_dirs = Util.walk_dir(dir, only_dirs = true)
  ! isempty(res_dirs) && push!(LOAD_PATH, res_dirs...)

  unique(LOAD_PATH)

  nothing
end

"""
    load_acl(dir::String) :: Dict{Any,Any}

Loads the ACL file associated with the invoked `controller` and returns the rules.
"""
function load_acl(dir::String) :: Dict{Any,Any}
  file_path = joinpath(dir, Genie.GENIE_AUTHORIZATOR_FILE_NAME)
  isfile(file_path) ? YAML.load(open(file_path)) : Dict{Any,Any}()
end

"""
    load_configurations() :: Void

Loads (includes) the framework's configuration files.
"""
function load_configurations() :: Void
  loggers_path = abspath("$(Genie.CONFIG_PATH)/loggers.jl")
  if isfile(loggers_path)
    include(loggers_path)
    is_dev() && Revise.track(loggers_path)
  end

  secrets_path = abspath("$(Genie.CONFIG_PATH)/secrets.jl")
  if isfile(secrets_path)
    include(secrets_path)
    is_dev() && Revise.track(secrets_path)
  end

  nothing
end

"""
    load_initializers() :: Void

Loads (includes) the framework's initializers.
"""
function load_initializers() :: Void
  dir = abspath(joinpath(Genie.CONFIG_PATH, "initializers"))

  if isdir(dir)
    f = readdir(dir)
    for i in f
      fi = joinpath(dir, i)
      include(fi)
      is_dev() && Revise.track(fi)
    end
  end

  nothing
end

"""
    secret_token() :: String

Wrapper around /config/secrets.jl SECRET_TOKEN `const`.
"""
function secret_token() :: String
  isdefined(App, :SECRET_TOKEN) ||
    throw("SECRET_TOKEN not configured - please make sure that you have a valid secrets.jl file. You can generate a new secrets.jl file with a random SECRET_TOKEN using Genie.REPL.write_secrets_file() or use the included /app/config/secrets.jl.example file as a model.")

  SECRET_TOKEN
end

if IS_IN_APP
  load_configurations()
  load_initializers()
  load_libs()
  load_resources()
end

end


text:
module AstroBase

using Reexport

@reexport using AstroTime

struct AstroException <: Exception
    msg::String
end

Base.showerror(io::IO, err::AstroException) = print(io, err.msg)

include(joinpath("constants", "Constants.jl"))
include(joinpath("interfaces", "Interfaces.jl"))
include(joinpath("util", "Util.jl"))
include(joinpath("bodies", "Bodies.jl"))
include(joinpath("earth_attitude", "EarthAttitude.jl"))
include(joinpath("two_body", "TwoBody.jl"))
include(joinpath("ephemerides", "Ephemerides.jl"))
include(joinpath("n_body", "NBody.jl"))
include(joinpath("frames", "Frames.jl"))
include(joinpath("coords", "Coords.jl"))
include(joinpath("astrometry", "Astrometry.jl"))

@reexport using .Interfaces
@reexport using .Constants
@reexport using .Util
@reexport using .Bodies
@reexport using .Ephemerides
@reexport using .EarthAttitude
@reexport using .TwoBody
@reexport using .Frames
@reexport using .Coords
@reexport using .Astrometry

end # module


text:
struct ViewerDisplay <: Display
end

function pngwrite(x, filename)
    b64 = open(filename, "w")
    show(IOContext(b64, :limit=>true), MIME"image/png"(), x)
    close(b64)
end

function svgwrite(x, filename)
    content = limitstringmime(MIME("image/svg+xml"), x)
    js = joinpath(".", "library", "plotly-latest.min.js")
    script = "<script src=\"" * js * "\"></script>"
    content = replace(content, r"<script.*plotly.*</script>", script)
    write(filename, content)
end

plotwrite(::MIME"image/png", x, filename) = pngwrite(x, filename)
plotwrite(::MIME"image/svg+xml", x, filename) = svgwrite(x, filename)

tempplot(::MIME"image/png") = String(RCall.reval("tempfile(fileext = '.png')"))
tempplot(::MIME"image/svg+xml") = String(RCall.reval("tempfile(fileext = '.html')"))

function display(d::ViewerDisplay, m::MIME, p)
    fname = tempplot(m)
    plotwrite(m, p, fname)
    RCall.reval("viewer = getOption('viewer', utils::browseURL)")
    RCall.rcall(:viewer, fname)
end

function display(d::ViewerDisplay, p)
    try
        display(d, MIME("image/svg+xml"), p)
        return
        display(d, MIME("image/png"), p)
        return
    catch e;
    end
    throw(MethodError(display, [d, p]))
end

const viewer_display = ViewerDisplay()

max_stars_repo_name: pthariensflame/JuliaCall

max_issues_repo_name: tinagui/JuliaCall

max_forks_repo_name: tinagui/JuliaCall


text:
using OffsetArrays
using Polynomials
using NCDatasets
using Interpolations
using DocStringExtensions
using Plots
using Parameters
# using DataInterpolations

include("types.jl")
include("ils_tables_io.jl")
include("prepare_ils.jl")
include("instrument_convolutions.jl")

FT = Float32

using NCDatasets
using Polynomials
file = "/net/fluo/data2/groupMembers/cfranken/oco2_L1bScGL_18689a_180105_B10003r_200316045235.h5"
ils_Δ, ils_in, dispersion = read_ils_table(file, "src/InstrumentOperators/json/ils_oco2.json");
band = 1
footprint = 1
# Footprint and band index:
extended_dims  = [footprint,band]
dispPoly = Polynomial(view(dispersion, :, extended_dims...))

ν = FT.(dispPoly.(0:1015))
# resolution (μm)
res = FT(0.001/1e3)
# Max range of ILS from center (µm)
kernel_range = FT(0.5/1e3)

# Kernel grid
grid_x = -kernel_range:res:kernel_range
ind_out = collect(1:1016);

ils_pixel   = prepare_ils_table(grid_x, ils_in, ils_Δ,extended_dims)
oco2_kernel = VariableKernelInstrument(ils_pixel, ν, ind_out)

ν_mod = collect(FT, 0.758:Float32(res / 1e3):0.771)
spectrum = zeros(FT, length(ν_mod))
spectrum[length(spectrum) ÷ 2]=1;
y_conv =  conv_spectra(oco2_kernel, ν_mod, spectrum);


text:
# Forward methods from abstracted domains to abstract interpreter #

satisfy(domain::AbstractedDomain, state::State, term::Term) =
    satisfy(domain.interpreter, domain, state, term)
satisfy(domain::AbstractedDomain, state::State, terms::AbstractVector{<:Term}) =
    satisfy(domain.interpreter, domain, state, terms)

satisfiers(domain::AbstractedDomain, state::State, term::Term) =
    satisfiers(domain.interpreter, domain, state, term)
satisfiers(domain::AbstractedDomain, state::State, terms::AbstractVector{<:Term}) =
    satisfiers(domain.interpreter, domain, state, terms)

evaluate(domain::AbstractedDomain, state::State, term::Term) =
    evaluate(domain.interpreter, domain, state, term)

initstate(domain::AbstractedDomain, problem::Problem) =
    initstate(domain.interpreter, domain, problem)
initstate(domain::AbstractedDomain, objtypes::AbstractDict) =
    initstate(domain.interpreter, domain, objtypes)
initstate(domain::AbstractedDomain, objtypes::AbstractDict, fluents) =
    initstate(domain.interpreter, domain, objtypes, fluents)

transition(domain::AbstractedDomain, state::State, action::Term; options...) =
    transition(domain.interpreter, domain, state, action; options...)
transition!(domain::AbstractedDomain, state::State, action::Term; options...) =
    transition!(domain.interpreter, domain, state, action; options...)

available(domain::AbstractedDomain, state::State, action::Action, args) =
    available(domain.interpreter, domain, state, action, args)
available(domain::AbstractedDomain, state::State) =
    available(domain.interpreter, domain, state)

execute(domain::AbstractedDomain, state::State, action::Action, args; options...) =
    execute(domain.interpreter, domain, state, action, args; options...)
execute!(domain::AbstractedDomain, state::State, action::Action, args; options...) =
    execute!(domain.interpreter, domain, state, action, args; options...)

relevant(domain::AbstractedDomain, state::State, action::Action, args) =
    relevant(domain.interpreter, domain, state, action, args)
relevant(domain::AbstractedDomain, state::State) =
    relevant(domain.interpreter, domain, state)

regress(domain::AbstractedDomain, state::State, action::Action, args; options...) =
    regress(domain.interpreter, domain, state, action, args; options...)
regress!(domain::AbstractedDomain, state::State, action::Action, args; options...) =
    regress!(domain.interpreter, domain, state, action, args; options...)

update!(domain::AbstractedDomain, state::State, diff::Diff) =
    update!(domain.interpreter, domain, state, diff)
update(domain::AbstractedDomain, state::State, diff::Diff) =
    update(domain.interpreter, domain, state, diff)

widen!(domain::AbstractedDomain, state::State, diff::Diff) =
    widen!(domain.interpreter, domain, state, diff)
widen(domain::AbstractedDomain, state::State, diff::Diff) =
    widen(domain.interpreter, domain, state, diff)


text:
using Documenter
using MultiStateSystems

makedocs(
    modules     = [MultiStateSystems],
    format      = Documenter.HTML(mathengine = Documenter.MathJax()),
    sitename    = "MultiStateSystems.jl",
    authors     = "Tom Van Acker",
    pages       = [ "Home"              => "index.md",
                    "Getting Started"   => "quickguide.md",
                    "DSL Manual"        =>
                        [ "State-Transition Diagram"=> "std.md",
                          "Distributions"           => "distribution.md",
                          "Generating Functions"    => "ugf.md",
                          "Network"                 => "network.md"],
                    "Models Manual"     =>
                        [ "Dependence"              => "dependence.md"
                          "Stochastic Processes"    => "processes.md"
                          "Generating Operators"    => "ugo.md"]
                  ]
)

deploydocs(
     repo = "github.com/timmyfaraday/MultiStateSystems.jl.git"
)


text:
include("libavformat_h.jl")

include("avformat.jl")
include("avio.jl")


text:
"""
$(SIGNATURES)

Compute the minimum sample size needed for the desired power, type I error rate and effect size

Arguments
---------
* `Test` -- The test for calculating the power
* `std` -- Standard deviation
* `power` -- Power
* `alpha` -- Desired type I error rate
* `side` -- Side of the test, "two" is two-sided, "one" is one-sided
"""
function sample_size{T <: TrialTest}(
    test::T;
    power::Union{Real, Void} = nothing,
    std::Union{Real, Tuple{Real, Real}, Void} = nothing,
    alpha::Real = 0.05,
    side::String = "two",
)
    # Transcode side
    side = lowercase(side)
    # Check arguments
    check_args_sample_size(test, power = power, std = std, alpha = alpha, side = side)
    # Calculate sample size
    gap(n::Real) = analytic_power(test, n = n, std = std, alpha = alpha, side = side) - power
    # Round to the minimum integer that is larger than the root
    # print(fzero(gap, 2, MAXSAMPLESIZE))
    return round(fzero(gap, 2, MAXSAMPLESIZE), RoundUp)

end


text:
# tests the Pumi interface to FFD

push!(LOAD_PATH, "../src/")

using FFD
using PumiInterface
using PdePumiInterface
using SummationByParts
using ODLCommonTools
using Base.Test

function getTestMesh(dim::Integer, ::Type{Tmsh}; order::Integer=1) where Tmsh
# dim: 2D or 3D
# Tmsh: datatype for mesh arrays
# order: linear or quadratic coordinate field
  degree = 1
  shape_type = 2

  if dim == 2
    sbp = getTriSBPOmega(degree=degree)
  sbpface = TriFace{Float64}(degree, sbp.cub, sbp.vtx)
  elseif dim == 3
    sbp = getTetSBPOmega(degree=degree)
    sbpface = TetFace{Float64}(degree, sbp.cub, sbp.vtx)
  end

  opts = PdePumiInterface.get_defaults()
  if order == 1
    opts["smb_name"] = "meshes/square5.smb"
  else
    opts["smb_name"] = "meshes/square_quad5.smb"
  end
  opts["dmg_name"] = ".null"
  opts["order"] = 1
  opts["coloring_distance"] = 2
  opts["numBC"] = 2
  opts["BC1"] = [0]
  opts["BC2"] = [1,2,3]
  opts["BC1_name"] = "fakeBC"
  opts["BC2_name"] = "fake2BC"

  if dim == 3
    opts["dimensions"] = 3
    opts["BC2"] = [1,2,3,4,5]
    if order == 1
      opts["smb_name"] = "meshes/tet4.smb"
    else
      opts["smb_name"] = "meshes/tet4_quad.smb"
    end
    face_verts = SummationByParts.SymCubatures.getfacevertexindices(sbp.cub)
    face_verts = SummationByParts.SymCubatures.getfacevertexindices(sbp.cub)
    edge_verts = [1 2 1 1 2 3;  # TODO: SBP should provide this
                  2 3 3 4 4 4]
    topo2 = ElementTopology2()   #TODO: make this the correct one for SBP
    topo = ElementTopology{3}(face_verts, edge_verts, topo2=topo2)

#    topo = ElementTopology{3}(face_verts)
  end

  if dim == 2
    mesh = PumiMeshDG2(Tmsh, sbp, opts, sbpface, dofpernode=1,
                      shape_type=shape_type)
  else
    mesh = PumiMeshDG3(Tmsh, sbp, opts, sbpface, topo, dofpernode=1,
                      shape_type=shape_type)
  end

  return mesh, sbp, opts
end

function test_surface(map, mesh)

  @testset "----- Testing FFD Surface Evaluation -----" begin
    vertices_orig = zeros(Complex128, mesh.dim, map.numFacePts)
    evalSurface(map, vertices_orig)

    cp_xyz = getControlPoints(map)
    # Rigid body rotation
    theta = -20*pi/180  # Rotate wall coordinates by 10 degrees
    phi = -10*pi/180  # rotation about x axis (3D only)
    if mesh.dim == 2
      rotMat = [cos(theta) -sin(theta)
                sin(theta) cos(theta)] # Rotation matrix

      for i=1:map.nctl[2]
        for j=1:map.nctl[1]
          cp_xyz[:, j, i] = rotMat*cp_xyz[:, j, i]
        end
      end
    else
      rotMat1 = [cos(theta) -sin(theta) 0
                 sin(theta) cos(theta)  0
                 0          0           1] # Rotation matrix
      rotMat2 = [1 0         0        ;
                 0 cos(phi) -sin(phi) ;
                 0 sin(phi)  cos(phi) ]

      rotMat = rotMat1*rotMat2
      # Rotate the control points
      for k = 1:map.nctl[3]
        for j = 1:map.nctl[2]
          for i = 1:map.nctl[1]
            cp_xyz[:,i,j,k] = rotMat*cp_xyz[:,i,j,k]
          end
        end
      end

    end

    # Rigid body translation
    xfac = 0.2
    yfac = 0.3
    if mesh.dim == 3
      zfac = 0.4
    else
      zfac = 0.0
    end

    if mesh.dim == 2
      cp_xyz[1, :, :] += xfac
      cp_xyz[2, :, :] += yfac
    else
      cp_xyz[1,:,:,:] += xfac
      cp_xyz[2,:,:,:] += yfac
      cp_xyz[3,:,:,:] += zfac
    end

    # set new control point values
    setControlPoints(map, cp_xyz)

    vertices = zeros(Complex128, mesh.dim, map.numFacePts)
    evalSurface(map, vertices)

    for i=1:map.numFacePts
      if mesh.dim == 2
        verts_exact = rotMat*vertices_orig[:, i] + [xfac; yfac]
      else
        verts_exact = rotMat*vertices_orig[:, i] + [xfac; yfac; zfac]
      end

      @test isapprox( norm(verts_exact - vertices[:, i]), 0.0) atol=1e-13
    end

#=
    rotMat2 = rotMat[1:2, 1:2]
    for bc=1:length(vertices)
      verts_bc = vertices[bc]
      verts_orig_bc = vertices_orig[bc]

      for i=1:size(verts_bc, 3)
        for j=1:size(verts_bc, 2)
          if mesh.dim == 2
            verts_exact = rotMat[1:2, 1:2]*verts_orig_bc[:, j, i] + [xfac; yfac]
          else
            verts_exact = rotMat*verts_orig_bc[:, j, i] + [xfac; yfac; zfac]
          end

          @test isapprox( norm(verts_exact - verts_bc[:, j, i]), 0.0) atol=1e-13
        end
      end
    end
=#

    # test getControlPoints
    cp_xyz2 = getControlPoints(map)
    fill!(cp_xyz2, 0.0)
    getControlPoints(map, cp_xyz2)

    if mesh.dim == 2
      @test isapprox( maximum(abs.(cp_xyz2 - map._cp_xyz[1:2, :, :, 1])), 0.0) atol=1e-13
    else
      @test isapprox( maximum(abs.(cp_xyz2 - map._cp_xyz)), 0.0) atol=1e-13
    end

  end  # end facts block
  return nothing
end

function test_jac2(map::PumiMapping{Tffd}, mesh) where Tffd

  @testset "----- Testing Jacobian-vector product -----" begin

    Xs_dot = zeros(Tffd, mesh.dim, map.numFacePts)
    Xs_dot2 = zeros(Xs_dot)

    vertices_orig = zeros(Xs_dot)
    evalSurface(map, vertices_orig)

    if Tffd <: Complex
      # test that map.cs produces the same coordinates as map
      map_cs = map.map_cs
      copy!(map_cs._cp_xyz, map._cp_xyz)
      vertices_cs = zeros(Xs_dot)
      evalSurface(map, vertices_cs)

      @test isapprox( maximum(abs.(vertices_cs - vertices_orig)), 0.0) atol=1e-13

      for i=1:length(map.cp_xyz)
        map._cp_xyz[i] += 1
      end

      copy!(map_cs._cp_xyz, map._cp_xyz)
      evalSurface(map, vertices_orig)
      evalSurface(map_cs, vertices_cs)
      @test isapprox( maximum(abs.(vertices_cs - vertices_orig)), 0.0) atol=1e-13
    end

    for i=1:10
      fill!(Xs_dot, 0.0)
      fill!(Xs_dot2, 0.0)

      Xcp = getControlPoints(map)
      _Xcp_dot = getControlPoints(map)
      rand!(_Xcp_dot)
      Xcp_dot = real(_Xcp_dot)

      # finite difference
      h = 1e-7
      for j=1:length(Xcp_dot)
        Xcp[j] += Xcp_dot[j]*h
#        map.cp_xyz[j] += Xcp_dot[j]*h
      end

      setControlPoints(map, Xcp)

      vertices_new = zeros(vertices_orig)
      evalSurface(map, vertices_new)

      for j=1:length(Xcp_dot)
        Xcp[j] -= Xcp_dot[j]*h
#        map.cp_xyz[j] -= Xcp_dot[j]*h
      end

      setControlPoints(map, Xcp)

      for j=1:length(Xs_dot)
        Xs_dot[j] = (vertices_new[j] - vertices_orig[j])/h
      end

      # complex step
      evaldXdControlPointProduct(map, Xcp_dot, Xs_dot2)

      for j=1:length(Xs_dot)
        @test isapprox( Xs_dot[j], Xs_dot2[j]) atol=1e-7
      end
    end

  end  # end facts block

  return nothing
end

function test_jac(map, mesh)

  @testset "----- Testing Transposed Jacobian-vector product -----" begin

    # construct entire jacobian
    cp_xyz = getControlPoints(map)
    jac = zeros(mesh.dim*map.numFacePts, length(cp_xyz))
    jac2 = zeros(jac)

    h = 1e-20
    pert = Complex128(0, h)
    Xs = zeros(Complex128, mesh.dim, map.numFacePts)
    for i=1:length(cp_xyz)
      cp_xyz[i] += pert
      setControlPoints(map, cp_xyz)
      evalSurface(map, Xs)

      for j=1:length(Xs)
        jac[j, i] = imag(Xs[j])/h
      end

      cp_xyz[i] -= pert
      setControlPoints(map, cp_xyz)
    end

      # compute reverse mode jacobian
      Xs_dot = zeros(Complex128, mesh.dim, map.numFacePts)
      B_dot = zeros(cp_xyz)
      for i=1:length(Xs_dot)
        Xs_dot[i] = 1
        fill!(B_dot, 0.0)
        evaldXdControlPointTransposeProduct(map, Xs_dot, B_dot)

        for j=1:length(B_dot)
          jac2[i, j] = real(B_dot[j])
        end
        Xs_dot[i] = 0
      end

      @test isapprox( norm(jac - jac2), 0.0) atol=1e-13

    println("\nTesting product")
    # test the product using the quadratic form dXs^T dXs/dXcp dXcp
    dXs = real(rand(size(Xs)))
    dXcp = real(rand(size(cp_xyz)))

    # forward product
    for i=1:length(cp_xyz)
      cp_xyz[i] += pert*dXcp[i]
    end
    setControlPoints(map, cp_xyz)
    fill!(Xs, 0.0)
    evalSurface(map, Xs)
    val1 = sum(dXs .* imag(Xs)/h)  # contract with dXs

    for i=1:length(cp_xyz)
      cp_xyz[i] -= pert*dXcp[i]
    end
    setControlPoints(map, cp_xyz)

    # transposed product
    fill!(B_dot, 0.0)
    evaldXdControlPointTransposeProduct(map, dXs, B_dot)
    val2 = sum(dXcp .* B_dot)

    @test isapprox( val1, val2) atol=1e-13

  end  # end facts block

  return nothing
end

function runtests()

  mesh, sbp, opts = getTestMesh(2, Complex128)

  # Free Form deformation parameters
  ndim = 2
  order = [3,3,2]  # Order of B-splines in the 3 directions
  nControlPts = [6,6,2]
  offset = [0.25, 0.25, 0.25]
  full_geom = false
  bc_nums = [1] # = opts["BC1"]
#  @assert length(geom_faces) == 1

  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_surface(map, mesh)

  # test quadratic
  mesh, sbp, opts = getTestMesh(2, Complex128, order=2)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_surface(map, mesh)

  # test jacobian2
  mesh, sbp, opts = getTestMesh(2, Complex128)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_jac2(map, mesh)

  # test jacbian2 quadratic
  mesh, sbp, opts = getTestMesh(2, Complex128, order=2)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_jac2(map, mesh)

  # test with real values map

  mesh, sbp, opts = getTestMesh(2, Float64)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_jac2(map, mesh)

  # test jacobian1
  # make new objects in case test_surface modified the old ones
  mesh, sbp, opts = getTestMesh(2, Complex128)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_jac(map, mesh)

  # test quadratic jacobian1
  mesh, sbp, opts = getTestMesh(2, Complex128, order=2)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_jac(map, mesh)

  # test with one boundary condition with several geometric entities
  bc_nums = [2]

  mesh, sbp, opts = getTestMesh(2, Complex128)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_surface(map, mesh)

  mesh, sbp, opts = getTestMesh(2, Complex128)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)
  test_jac(map, mesh)

  # test multi-geometry entity on quadratic
  mesh, sbp, opts = getTestMesh(2, Complex128, order=2)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_surface(map, mesh)

  mesh, sbp, opts = getTestMesh(2, Complex128, order=2)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)
  test_jac(map, mesh)

  # test multiple BCs
  bc_nums = [1, 2]
  mesh, sbp, opts = getTestMesh(2, Complex128)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_surface(map, mesh)

  mesh, sbp, opts = getTestMesh(2, Complex128)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)
  test_jac2(map, mesh)

  mesh, sbp, opts = getTestMesh(2, Complex128)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)
  test_jac(map, mesh)

  # test multi-BC on quadratic
  bc_nums = [1, 2]
  mesh, sbp, opts = getTestMesh(2, Complex128, order=2)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_surface(map, mesh)

  mesh, sbp, opts = getTestMesh(2, Complex128, order=2)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)
  test_jac2(map, mesh)

  mesh, sbp, opts = getTestMesh(2, Complex128, order=2)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)
  test_jac(map, mesh)

  println("\nTesting 3D")
  # now test 3D
  mesh, sbp, opts = getTestMesh(3, Complex128)

  # Free Form deformation parameters
  ndim = 3
  order = [3,3,3]  # Order of B-splines in the 3 directions
  nControlPts = [6,6,6]
  offset = [0.25, 0.25, 0.25]
  full_geom = false
  bc_nums = [1] # opts["BC1"]
#  @assert length(geom_faces) == 1

  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_surface(map, mesh)

  # test surface quadratic
  mesh, sbp, opts = getTestMesh(3, Complex128, order=2)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_surface(map, mesh)

  # test jacobian
  mesh, sbp, opts = getTestMesh(3, Complex128)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_jac(map, mesh)

  # test jacobian quadratic
  mesh, sbp, opts = getTestMesh(3, Complex128, order=2)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_jac(map, mesh)

  # test multiple BCs
  bc_nums = [1, 2]

  # surface
  mesh, sbp, opts = getTestMesh(3, Complex128)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_surface(map, mesh)

  # surface quadratic
  mesh, sbp, opts = getTestMesh(3, Complex128, order=2)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_surface(map, mesh)

  # jacobian 2
  mesh, sbp, opts = getTestMesh(3, Complex128)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_jac2(map, mesh)

  # jacobian 2 quadratic
  mesh, sbp, opts = getTestMesh(3, Complex128, order=2)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_jac2(map, mesh)

  # jacobian
  mesh, sbp, opts = getTestMesh(3, Complex128)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_jac(map, mesh)

  # jacobian quadratic
  mesh, sbp, opts = getTestMesh(3, Complex128, order=2)
  map, box = initializeFFD(mesh, sbp, order, nControlPts, offset,
                           bc_nums)

  test_jac(map, mesh)

  return nothing
end

runtests()


text:
using MajorizationExtrema
using MajorizationExtrema: entropy
using Test, LinearAlgebra

function check_simplexpt(q)
    if eltype(q) <: Rational
        @test sum(q) == one(eltype(q))
    else
        @test sum(q) ≈ one(eltype(q))
    end
    @test all(x -> x >= 0, q)
end

function check_dm(ρ)
    @test tr(ρ) ≈ one(eltype(ρ))
    @test eigmin(ρ) >= -1e-8
end

# this likely doesn't sample uniformly
function sample_L1(d)
    u = randprobvec(2d)
    return [u[i] - u[i+d] for i = 1:d]
end

# this definitely does not sample uniformly
function sample_TV(p, ϵ)
    v = sample_L1(length(p))
    q = p + ϵ * v
    @. q = clamp(q, 0, 1)
    q .= q ./ sum(q)

    # Not sure if this is necessary, but a priori it is
    TV(p, q) > ϵ && return sample_TV(p, ϵ)

    return q
end

@testset "Types" begin
    @test @inferred(randdm(3)) isa Hermitian
    ρ = randdm(3)
    @test @inferred(majmin(ρ, 0.01)) isa Hermitian
    @test @inferred(majmax(ρ, 0.01)) isa Hermitian
    @test eltype(@inferred(majmin(ρ, 1 // 10))) == Complex{Float64}

    @test typeof(@inferred(randprobvec(3))) == Vector{Float64}
    @test typeof(@inferred(randprobvec(3, 100))) == Vector{Rational{Int}}
    @test typeof(@inferred(randprobvec(3, big(100)))) == Vector{Rational{BigInt}}
    @test typeof(@inferred(randprobvec(BigFloat, 3))) == Vector{BigFloat}

    for m in (majmin, majmax)
        p = randprobvec(3)
        @test typeof(@inferred(m(p, 0.01))) == Vector{Float64}
        @test typeof(@inferred(m(p, 1 // 10))) == Vector{Float64}
        @test typeof(@inferred(m(p, big(1 // 10)))) == Vector{BigFloat}

        p = randprobvec(3, 100)
        @test typeof(@inferred(m(p, 0.01))) == Vector{Float64}
        @test typeof(@inferred(m(p, 1 // 10))) == Vector{Rational{Int}}
        @test typeof(@inferred(m(p, big(1 // 10)))) == Vector{Rational{BigInt}}

        p = randprobvec(3, big(100))
        @test typeof(@inferred(m(p, 0.01))) == Vector{BigFloat}
        @test typeof(@inferred(m(p, 1 // 10))) == Vector{Rational{BigInt}}
        @test typeof(@inferred(m(p, big(1 // 10)))) == Vector{Rational{BigInt}}
    end

end

@testset "Quantum states" begin

    X = [1.0 2.0; 0.0 0.0]
    @test_throws ArgumentError majmin(X, 0.01)
    @test_throws ArgumentError majmax(X, 0.01)

    ρ = randdm(2)
    @test_throws ArgumentError ρ ≺ X
    @test_throws ArgumentError X ≺ ρ
    for d in (2, 3, 5)
        idmat = Matrix(1.0I, d, d)

        for ϵ in (1, 2, 1 - 1 / d + 0.01)
            ρ = randdm(d)
            check_dm(ρ)
            @test majmin(ρ, ϵ) ≈ idmat / d
            @test tr(majmax(ρ, ϵ)^2) ≈ 1.0
            @test majmin(ρ, ϵ) ≺ ρ
            @test ρ ≺ majmax(ρ, ϵ)
        end

        for ϵ in (0.01, 0.05, 0.1)
            ρ = randdm(d)
            check_dm(ρ)
            @test localbound(entropy, ρ, ϵ) <= ϵ * log2(d - 1) + entropy([ϵ, 1 - ϵ])

            p = randprobvec(d)
            U = randunitary(d)
            @test U * U' ≈ idmat
            @test U' * U ≈ idmat

            ρ = U * Diagonal(p) * U'
            ρ = (ρ + ρ') / 2
            check_dm(ρ)

            @test majmin(ρ, ϵ) ≈ U * Diagonal(majmin(p, ϵ)) * U'
            @test majmax(ρ, ϵ) ≈ U * Diagonal(majmax(p, ϵ)) * U'

            @test majmin(ρ, ϵ) ≺ ρ
            @test ρ ≺ majmax(ρ, ϵ)

            if tracedist(ρ, idmat / d) > ϵ
                @test tracedist(majmin(ρ, ϵ), ρ) ≈ ϵ
            else
                @test majmin(ρ, ϵ) ≈ idmat / d
            end

            p = randprobvec(d)
            q = randprobvec(d)
            ρ = U * Diagonal(p) * U'
            σ = U * Diagonal(q) * U'
            @test tracedist(ρ, σ) ≈ tracedist(Hermitian(ρ), Hermitian(σ))
            @test tracedist(ρ, σ) ≈ TV(p, q)
        end

    end
end

@testset "Test local bound" begin
    for d in [2, 3, 5, 20]
        p = randprobvec(d)
        entropy_p = entropy(p)
        for ϵ in [0.01, 0.01, 0.1]
            LB = localbound(entropy, p, ϵ)
            p_max = majmax(p, ϵ)
            p_min = majmin(p, ϵ)
            @test (LB ≈ entropy(p) - entropy(majmax(p, ϵ))) ||
                  (LB ≈ entropy(majmin(p, ϵ)) - entropy(p))
            for _ = 1:5
                r = sample_TV(p, ϵ)
                @test p_min ≺ r
                @test r ≺ p_max
                @test abs(entropy(r) - entropy_p) <= LB
            end
        end
    end
end

@testset "Basic checks" begin

    p1 = [1 // 3, 1 // 3, 1 // 3]
    p2 = [1 // 3, 1 // 3, 1 // 3, 1 // 3]
    @test_throws DimensionMismatch p1 ≺ p2
    @test_throws DimensionMismatch float.(p1) ≺ float.(p2)
    @test_throws DimensionMismatch TV(p1, p2)
    @test_throws DimensionMismatch TV(float.(p1), float.(p2))

    p = [1 // 3, 1 // 3, 1 // 3]
    q = [1 // 4, 1 // 4, 1 // 2]

    @test p ≺ q
    @test !(q ≺ p)
    @test TV(p, q) == 1 // 6

    @test majmin(q, 1 // 20) == [11 // 40, 11 // 40, 9 // 20]
    @test majmax(q, 1 // 20) == [1 // 5, 1 // 4, 11 // 20]
    @test majmax(q, 2) == [0 // 1, 0 // 1, 1 // 1]
    @test majmin(q, 2) == [1 // 3, 1 // 3, 1 // 3]

    p = float.(p)
    q = float.(q)

    @test p ≺ q
    @test !(q ≺ p)
    @test majmin(q, 1 / 20) ≈ [11 / 40, 11 / 40, 9 / 20]
    @test majmax(q, 1 / 20) ≈ [1 / 5, 1 / 4, 11 / 20]
    @test majmax(q, 2) ≈ [0, 0, 1]
    @test majmin(q, 2) ≈ [1 / 3, 1 / 3, 1 / 3]

    @test TV(p, q) ≈ 1 / 6

end

@testset "Tests with type $T" for T in [Float64, BigFloat, Rational{Int}, Rational{BigInt}]
    for d in [2, 3, 5, 10]
        δ = zeros(T, d)
        δ[1] = one(T)

        u = ones(T, d) / T(d)

        for ϵ in [T(1 // 30), T(1 // 20), T(1 // 4)]
            if T == Float64
                r = randprobvec(d)
                tol = 1e-8
            elseif T == BigFloat
                r = randprobvec(BigFloat, d)
                tol = 1e-25
            elseif T == Rational{Int}
                r = randprobvec(d, rand(10:10000))
                tol = zero(T)
            elseif T == Rational{BigInt}
                r = randprobvec(d, big(rand(10:10000)))
                tol = zero(T)
            end

            r_max = majmax(r, ϵ)
            r_min = majmin(r, ϵ)

            @test TV(r, r_max) ≈ 0.5 * norm(r - r_max, 1)
            @test TV(r, r_min) ≈ 0.5 * norm(r - r_min, 1)

            @test TV(r, r_max) <= ϵ + tol
            @test TV(r, r_min) <= ϵ + tol

            check_simplexpt(r)
            check_simplexpt(r_max)
            check_simplexpt(r_min)

            @test eltype(r) == T
            @test eltype(r_max) == T
            @test eltype(r_min) == T

            @test r ≺ δ
            @test r ≺ r_max
            @test r_min ≺ r
            @test u ≺ r
        end
    end
end

@testset "Random tests" begin
    for _ = 1:10
        q_rational = randprobvec(3, 1000)
        check_simplexpt(q_rational)

        q_float = float.(q_rational)
        check_simplexpt(q_float)

        for ϵ in (1 // 100, 2 // 100)
            @test majmin(q_rational, ϵ) ≈ majmin(q_float, ϵ)
            @test majmax(q_rational, ϵ) ≈ majmax(q_float, ϵ)
        end

        for ϵ in (1, 2, 2 // 3)
            @test majmin(q_rational, ϵ) == [1 // 3, 1 // 3, 1 // 3]
            @test majmin(q_float, ϵ) ≈ [1 / 3, 1 / 3, 1 / 3]

            @test sort(majmax(q_rational, ϵ)) == [0 // 1, 0 // 1, 1 // 1]
            @test sort(majmax(q_float, ϵ)) ≈ [0.0, 0.0, 1.0]
        end

        for T in (Float64, BigFloat)
            p, q = @inferred(MajorizationExtrema.randincomppair(T, 4))
            @test !(p ≺ q) && !(q ≺ p)
            @test !(p ≻ q) && !(q ≻ p)
            check_simplexpt(p)
            check_simplexpt(q)
            @test eltype(p) == T
            @test eltype(q) == T

            p, q = @inferred(MajorizationExtrema.randmajpair(T, 4))
            @test (p ≺ q) && (q ≻ p)
            check_simplexpt(p)
            check_simplexpt(q)
            @test eltype(p) == T
            @test eltype(q) == T
        end

    end
end

include("probvecmults.jl")
include("subentropy.jl")


text:
using Distributed;

println(workers())

@everywhere function test()
    id = myid()
    println(id)

    try
        open("/data/cache/$(id).txt", "w") do io
            println(io, "myid() = $id")
        end
    catch err
        println(err)
    end
end

ras = [@spawnat w test() for w in workers()]

println("ras: ", ras)

@time wait.(ras)


text:
using ExamplePackage, Test

@testset "Test my_f evaluation" begin
    # 2*2+3*1 = 7
    @test my_f(2,1) == 7
    @test my_f(3,2) == 12
    @test my_f(3,1) == 9
end

@testset "Test my_f derivative" begin
    @test derivative_of_my_f(2,1) == 2
end


text:
using Base.Filesystem
include( "touch.jl")

mutable struct folder
	baseSrcroot::String
	baseTgrroot::String

	cursrcroot::String
	curtgrroot::String

	curSUBname::String #,maybe folder or file or link
	nextFullSrcName::String
	nextFullTgrName::String
	folder() =new() #new(baseSrcroot, baseTgrRoot,cursrcroot,curtgrroot,curSUBname,curfullsrc, curfulltgr)

end
const channelOut=Channel( 10 )
function getProgress( )
	global channelOut
	isready(channelOut )
	a = take!(channelOut)
	return a
end
function showProgress(args...)
	global channelOut
	s1 =""
	for a in args
		s1 =string(s1, a,"\n" )
	end
	put!(channelOut,s1)
end

function procdir(folderobj::folder)
	if ispath( folderobj.nextFullTgrName )
			;
	else
		mkdir(folderobj.nextFullTgrName)
	end
end
function procfile(folderobj::folder)
	# return
		if Filesystem.ispath(folderobj.nextFullTgrName )
			srcfileDT = Filesystem.mtime(folderobj.nextFullSrcName )
			tgrfileDT=  Filesystem.mtime(folderobj.nextFullTgrName )
			delta = abs(srcfileDT - tgrfileDT)*1000
			if  delta > 1.0
						if srcfileDT  >  tgrfileDT
							showProgress("Newer file found: ,copying....:",folderobj.nextFullTgrName )
								Filesystem.cp( folderobj.nextFullSrcName,folderobj.nextFullTgrName,force=true )
								touchfile( folderobj.nextFullSrcName,folderobj.nextFullTgrName )
					    else
								#println("Older  file found: ",folderobj.nextFullTgrName )
								#### Filesystem.cp( folderobj.nextFullSrcName,folderobj.nextFullTgrName,force=true )
								#touchfile( folderobj.nextFullSrcName,folderobj.nextFullTgrName )
						end
			else
				 #println("Same File ....")
			end

		else
			showProgress("LOST  file found,copying:",folderobj.nextFullSrcName , " \nCopying ...:" , folderobj.nextFullTgrName )
			try
				Filesystem.cp( folderobj.nextFullSrcName,folderobj.nextFullTgrName )
				touchfile( folderobj.nextFullSrcName,folderobj.nextFullTgrName )
			catch ex1
				showProgress( ex1)
			end
		end
end

function procLink(folderobj::folder)
	#println("linkname:$(folderobj.nextFullTgrName)\r\n" )
	#if   true #ispath(folderobj.nextFullTgrName)
		rm(folderobj.nextFullTgrName, force=true, recursive=true)
	#end
		#Filesystem.cp( folderobj.nextFullSrcName,folderobj.nextFullTgrName,follow_symlinks=false,force=true )

		src_LinkTgr =readlink(  folderobj.nextFullSrcName )
		tgr_LinkTgr = replace( src_LinkTgr ,folderobj.baseSrcroot =>folderobj.baseTgrroot )
		symlink( tgr_LinkTgr , folderobj.nextFullTgrName );
		#println("    Link created! Link Tgr name : $tgr_LinkTgr"   )
		#touchfile( folderobj.nextFullSrcName,folderobj.nextFullTgrName )

end

function doprocFolder( folderobj::folder )
	subfolderObj =folder()
	subfolderObj.baseSrcroot=folderobj.baseSrcroot
	subfolderObj.baseTgrroot=folderobj.baseTgrroot
	subfolderObj.cursrcroot=folderobj.nextFullSrcName
	subfolderObj.curtgrroot= folderobj.nextFullTgrName

	if ispath( subfolderObj.baseTgrroot)
					;
					#return ;
	else
		mkdir(subfolderObj.baseTgrroot)
	end

	for file in  readdir( subfolderObj.cursrcroot )#,topdown=true,follow_symlinks=false )
		#for file in files
			subfolderObj.curSUBname= file
			# println("---",file ,"  ")
			subfolderObj.nextFullSrcName = Filesystem.joinpath(  subfolderObj.cursrcroot ,file)
			subfolderObj.nextFullTgrName = Filesystem.joinpath(  subfolderObj.curtgrroot ,file)

			if islink( subfolderObj.nextFullSrcName)
				procLink(subfolderObj)
			elseif isfile(subfolderObj.nextFullSrcName)
				procfile(subfolderObj)
			elseif isdir(subfolderObj.nextFullSrcName)
				procdir(subfolderObj)
				doprocFolder(subfolderObj)
			end
		#end
	end

end
function proCopyFolder( srcfolder,tgrfolder)
	folderobj=folder()
	folderobj.baseSrcroot =abspath( srcfolder)
	folderobj.baseTgrroot =abspath(tgrfolder)
	folderobj.cursrcroot=  folderobj.baseSrcroot
	folderobj.curtgrroot=  folderobj.baseTgrroot
	folderobj.nextFullSrcName =    folderobj.cursrcroot
	folderobj.nextFullTgrName =    folderobj.curtgrroot
	doprocFolder(folderobj)
end


text:
#=
using PowerSystems
using JuMP
base_dir = dirname(dirname(pathof(PowerSystems)))
include(joinpath(base_dir, "data/data_5bus_uc.jl"))
sys5 = PSY.System(nodes5, generators5, loads5_DA, branches5, nothing, 100.0);
using PowerSimulations
const PS = PowerSimulations

simple_reserve = PSY.StaticReserve("test_reserve", vcat(sys5.generators.thermal, sys5.generators.renewable[2]), 60.0, [gen.tech for gen in sys5.generators.thermal])
#simple_reserve = PSY.StaticReserve("test_reserve", sys5.generators.thermal, 60.0, [sys5.generators.thermal[1].tech])

@test try
    Net = PSI.CopperPlatePowerModel
    m = Model();
    netinjection = PSI.instantiate_network(Net, sys5);
    PSI.construct_device!(m, netinjection, ThermalGen, PSI.ThermalDispatch, Net, sys5);
    PSI.construct_device!(m, netinjection, RenewableGen, PSI.RenewableCurtail, Net, sys5);
    PSI.construct_device!(m, netinjection, ElectricLoad, PSI.InterruptibleLoad, Net, sys5);
    PSI.construct_network!(m, [(device=Branch, formulation=PSI.PiLine)], netinjection, Net, sys5)
    PSI.construct_service!(m, simple_reserve, PSI.RampLimitedReserve, [(device = ThermalGen, formulation =PSI.ThermalDispatch),
                                                              (device = RenewableGen, formulation = PSI.RenewableCurtail)],
                                                              sys5)
    m.obj_dict
true finally end

@test try
    Net = PSI.CopperPlatePowerModel
    m = Model();
    netinjection = PSI.instantiate_network(Net, sys5);
    PSI.construct_device!(m, netinjection, ThermalGen, PSI.ThermalStandardUnitCommitment , Net, sys5);
    PSI.construct_device!(m, netinjection, RenewableGen, PSI.RenewableCurtail, Net, sys5);
    PSI.construct_device!(m, netinjection, ElectricLoad, PSI.InterruptibleLoad, Net, sys5);
    PSI.construct_network!(m, [(device=Branch, formulation=PSI.PiLine)], netinjection, Net, sys5)
    PSI.construct_service!(m, simple_reserve, PSI.RampLimitedReserve, [(device = ThermalGen, formulation =PSI.ThermalStandardUnitCommitment ),
                                                              (device = RenewableGen, formulation = PSI.RenewableCurtail)],
                                                              sys5)
    m.obj_dict
true finally end
=#


text:
__precompile__(true)

module InsarLOS

include("./projections.jl")

using Sario
import MapImages

import Glob
using SQLite
using LinearAlgebra
using HDF5

const EARTH_SMA = 6378137.0
const EARTH_E2 = 0.0066943799901499996

const SOL = 299792458

const MAP_FILENAME = "los_map.h5"

"""
    get_los_enu(lat_lons::AbstractArray{<:AbstractFloat, 2}; xyz_los_vecs=nothing, dbfile=nothing, param_dict=nothing)

Calculate the LOS vector in ENU from a lat/lon Array
Can also take in the pre-computed xyz los vectors

`lat_lons` is `2 x N` for N different points
"""
function get_los_enu(lat_lons::AbstractArray{<:AbstractFloat, 2}; xyz_los_vecs=nothing,
                     dbfile=nothing, param_dict=nothing)
    if isnothing(xyz_los_vecs)
        xyz_los_vecs = calculate_los_xyz(lat_lons, dbfile=dbfile, param_dict=param_dict)
    end
    # In projections file:
    return convert_xyz_latlon_to_enu(reshape(lat_lons, 2, :),
                                     reshape(xyz_los_vecs, 3, :))
end

function get_los_enu(lat_lons::AbstractArray{<:AbstractFloat, 1}; xyz_los_vecs=nothing,
                     dbfile=nothing, param_dict=nothing)
    return get_los_enu(reshape(lat_lons, :, 1), xyz_los_vecs=xyz_los_vecs, dbfile=dbfile, param_dict=param_dict)
end

filepath(full_path) = joinpath(splitpath(full_path)[1:end-1]...)
#
# Start of fortran converted functions:

"""Calculate the line of sight vector (from satellite to ground) in X,Y,Z coordinates at one lat/lon point"""
function calculate_los_xyz(lat::T, lon::T; dbfile::Union{String, Nothing}=nothing, param_dict=nothing) where {T<:AbstractFloat}
    if isnothing(param_dict)
        if isnothing(dbfile)
            dbfile_list = Glob.glob("*.db*")
            dbfile = dbfile_list[1]
        end
        param_dict = load_all_params(dbfile)
    end
    xyz_ground = _compute_xyz_ground(lat, lon)

    orbinfo_filename = param_dict["orbinfo"]  # The .db file doesn't save path
    dbpath = filepath(dbfile)
    orbinfo_file = joinpath(dbpath, orbinfo_filename)

    timeorbit, xx, vv = read_orbit_vector(orbinfo_file)

    for idx=1:param_dict["azimuthBursts"]
        idx = 9
        vecr = compute_burst_vec(xyz_ground, param_dict, idx, timeorbit, xx, vv)
        if !isnothing(vecr)
            return vecr
        end
    end
    return nothing
end

function calculate_los_xyz(lat::T, lon::T, dem, demrsc, param_dict, timeorbit, xx, vv) where {T<:AbstractFloat}
    xyz_ground = _compute_xyz_ground(lat, lon, dem, demrsc)
    for idx=1:param_dict["azimuthBursts"]
        idx = 9
        vecr = compute_burst_vec(xyz_ground, param_dict, idx, timeorbit, xx, vv)
        if !isnothing(vecr)
            return vecr
        end
    end
    return nothing
end

function calculate_los_xyz(lat_lon::AbstractArray{<:AbstractFloat, 1}; dbfile=nothing, param_dict=nothing)
    lat, lon = lat_lon
    return calculate_los_xyz(lat, lon, dbfile=dbfile, param_dict=param_dict)
end

function calculate_los_xyz(lat_lon_vecs::AbstractArray{<:AbstractFloat, 2}; dbfile=nothing, param_dict=nothing)
    num_vecs = size(lat_lon_vecs, 2)
    xyz_vecs = Array{eltype(lat_lon_vecs)}(undef, (3, ))
    for i=1:num_vecs
        lat, lon = lat_lon_vecs[:, i]
        xyz_vecs[:, i] .= calculate_los_xyz(lat, lon, dbfile=dbfile, param_dict=param_dict)
    end
    return xyz_vecs
end

function _compute_xyz_ground(lat, lon)
    # TODO: do i wanna get rid of this "params" file?
    # dem_file, dem_rsc_file = readlines("params")

    # TODO: what directory should this be??
    dem_rsc_file = Sario.find_rsc_file(directory="../")
    demrsc = Sario.load(dem_rsc_file)

    row, col = MapImages.latlon_to_rowcol(demrsc, lat, lon)
    # println("($lat, $lon) is at ($row, $col)")

    dem_file = replace(dem_rsc_file, ".rsc" => "")
    # Load just the 1 value from the DEM
    dem_height = Sario.load(dem_file, (row, col))

    llh = [ deg2rad(lat), deg2rad(lon), dem_height ]
    xyz_ground = llh_to_xyz(llh)
    return xyz_ground
end

function _compute_xyz_ground(lat, lon, dem, demrsc)
    row, col = MapImages.nearest_pixel(demrsc, lat, lon)

    # Load just the 1 value from the DEM
    dem_height = dem[row, col]

    llh = [ deg2rad(lat), deg2rad(lon), dem_height ]
    xyz_ground = llh_to_xyz(llh)
    return xyz_ground
end

function compute_burst_vec(xyz_ground, param_dict, idx, timeorbit, xx, vv)
    dtaz = 1.0 / param_dict["prf"]  # Nazlooks / prf
    tstart = param_dict["azimuthTimeSeconds$idx"]
    tend  = tstart + (param_dict["linesPerBurst"] - 1) * dtaz
    tmid = (tstart + tend) / 2.0

    # println("Burst $idx, Start, stop Acquisition time: $tstart,$tend")

    rngstart = param_dict["slantRangeTime"] * SOL / 2.0
    dmrg = SOL / 2.0 / param_dict["rangeSamplingRate"] # Nnum rnglooks * drho
    rngend = rngstart + (param_dict["samplesPerBurst"]-1)*dmrg
    rngmid = 0.50*(rngstart+rngend)

    # TODO: Does it matter if the lat/lon is within this specific burst?
    # If we interpolate to the middle it seems to always be the same
    #
    # latlons = bounds(tstart, tend, rngstart, rngend, timeorbit, xx, vv)
    # if lat > latlons[1] || lat < latlons[2]
    #     print("$lat not within bounds $latlons")
    #     return nothing
    # end

    xyz_mid, vel_mid = intp_orbit(timeorbit, xx, vv, tmid)

    # println("Satellite midpoint time,position,velocity: $tmid $xyz_mid $vel_mid")

    tline, range = orbitrangetime(xyz_ground, timeorbit, xx, vv, tmid, xyz_mid, vel_mid)
    if isnothing(tline) || isnothing(range)
        println("Failed on burst $idx")
        return nothing
    end

    satx, satv = intp_orbit(timeorbit, xx, vv, tline)
    # TESTING: see how much the East and North change at beginning/end of orbit
    # satx = xx[:, end]
    # satv = vv[:, end]

    # @show satx, satv
    # Note: pointing AWAY from satellite (hence -satx)
    dr = xyz_ground - satx
    vecr = dr / range

    return vecr

end

function load_table_params(dbfile, param_list, param_types=nothing; tablename="file")
    db = SQLite.DB(dbfile)

    # TODO: How do you really do string interpolation for this library :\
    param_list_str = join(param_list, "\",\"")
    querystr = """SELECT * FROM $tablename WHERE name IN ("$param_list_str") """

    table = SQLite.Query(db, querystr)
    if isnothing(param_types)
        return Dict(row.name => row.value for row in table)
    else
        return Dict(row.name => parse(param_types[row.name], row.value)
                    for row in table)
    end
end

function load_all_params(dbfile; tablename="file")
    param_dict = Dict{String, Any}()
    param_list = [
        "azimuthBursts",
        "linesPerBurst",
        "samplesPerBurst",
        "prf",
        "slantRangeTime",
        "rangeSamplingRate",
    ]
    # To parse into the correct types, since all are strings
    param_types = Dict(
        "azimuthBursts" => Int,
        "linesPerBurst" => Int,
        "samplesPerBurst" => Int,
        "prf" => Float64,
        "rangeSamplingRate" => Float64,
        "slantRangeTime" => Float64,
    )
    d = load_table_params(dbfile, param_list, param_types, tablename=tablename)
    merge!(param_dict, d)

    # For this, one of each of these per burst up to "azimuthBursts"
    for idx = 1:param_dict["azimuthBursts"]
        burst_params = [
            "azimuthTimeSeconds$idx",
        ]
        types = Dict(
            "azimuthTimeSeconds$idx" => Float64,
        )
        d = load_table_params(dbfile, burst_params, types, tablename=tablename)
        merge!(param_dict, d)
    end

    # Finally, the orbinfo file is a string
    merge!(param_dict,
           load_table_params(dbfile, ["orbinfo"], tablename=tablename))
    return param_dict
end

function read_orbit_vector(orbtiming_file)
    lines = readlines(orbtiming_file)
    # timefirst, timeend = map(line -> parse(Float64, line), lines[1:2])
    nlines, num_state_vec = map(line -> parse(Int, line), lines[3:4])

    timeorbit = zeros(num_state_vec)
    xx = Array{Float64, 2}(undef, (3, num_state_vec))
    vv = similar(xx)
    # Note: we don't use the accel
    for i = 1:num_state_vec

        floats = map(num -> parse(Float64, num), split(lines[4 + i]))
        timeorbit[i] = floats[1]
        xx[:, i] = floats[2:4]
        vv[:, i] = floats[5:7]
    end
    return timeorbit, xx, vv
end

function llh_to_xyz(llh::AbstractArray{<:AbstractFloat, 1})
    lat, lon, h = llh

    rad_earth = EARTH_SMA/sqrt(1.0 - EARTH_E2*sin(lat)^2)

    xyz = similar(llh)
    xyz[1] = (rad_earth + h)*cos(lat)*cos(lon)
    xyz[2] = (rad_earth + h)*cos(lat)*sin(lon)
    xyz[3] = (rad_earth*(1.0 - EARTH_E2) + h) * sin(lat)
    return xyz
end

function orbitrangetime(xyz, timeorbit, xx, vv, tline0, satx0, satv0, tol=1e-6)
    # starting state
    tline = tline0
    satx = satx0
    satv = satv0

    idx, max_iter = 1, 100
    tprev = tline + 1  # Need starting guess
    while abs(tline - tprev) > tol && idx < max_iter
        tprev = tline

        dr = xyz - satx
        range = norm(dr, 2)

        fn = dr' * satv
        fnprime = -norm(satv, 2)^2

        tline = tline - fn / fnprime

        satx, satv = intp_orbit(timeorbit, xx, vv, tline)
        if isnothing(satx) || isnothing(satv)
            return nothing, nothing
        end
        idx += 1
    end
    if abs(tline - tprev) > tol
        println("Warning: orbitrangetime didn't converge within $tol (residual $(abs(tline - tprev))")
    end

    dr = xyz - satx
    range = norm(dr, 2)

    return tline, range
end

function intp_orbit(timeorbit, xx, vv, time)
    if isinf(time) || isnan(time)
        return nothing, nothing
    end
    ilocation = (time - timeorbit[1]) / (timeorbit[2]-timeorbit[1])
    # @show time
    # @show ilocation
    ilocation = round(Int, clamp(ilocation, 2, length(timeorbit) - 2))

    xyz_mid, vel_mid = orbithermite(xx[:,ilocation-1:end], vv[:,ilocation-1:end],
                                    timeorbit[ilocation-1:end], time)
    # xyz_mid, vel_mid = orbithermite(xx, vv,timeorbit, time)

    return xyz_mid, vel_mid
end

"""orbithermite - hermite polynomial interpolation of orbits"""
function orbithermite(x,v,t,time)
# inputs
#  x - 3x4 matrix of positions at four times
#  v - 3x4 matrix of velocities
#  t - 4-vector of times for each of the above data points
#  time - time to interpolate orbit to
#
# outputs
#  xx - position at time time
#  vv - velocity at time time
    n = 4

    h, hdot = zeros(n), zeros(n)
    f0, f1 = zeros(n), zeros(n)
    g0, g1 = zeros(n), zeros(n)
    sum, product = 0, 0

    for i=1:n
        f1[i]=time-t[i]
        sum=0.
        for j=1:n
            if j != i
                sum += 1.0/(t[i]-t[j])
            end
        end
        f0[i]=1.0 - 2.0 * (time-t[i])*sum
    end

    for i=1:n
        product=1.0
        for k=1:n
            if k != i
                product *= (time-t[k])/(t[i]-t[k])
            end
        end
        h[i]=product

        sum=0.0
        for j=1:n
            product=1.0
            for k=1:n
               if k != i && k != j
                   product *= (time-t[k])/(t[i]-t[k])
               end
            end
            if j != i
                sum += 1.0 / (t[i]-t[j])*product
            end
        end
        hdot[i]=sum
    end

    for i=1:n
        g1[i]=h[i] + 2.0 * (time-t[i]) * hdot[i]
        sum=0.
        for j=1:n
            if i != j
                sum += 1.0 / (t[i]-t[j])
            end
        end
        g0[i]=2.0 *(f0[i]*hdot[i]-h[i]*sum)
    end

    xout = zeros(3)
    vout = zeros(3)
    for k=1:3
        sum=0.0
        for i=1:n
            sum += (x[k,i]*f0[i]+v[k,i]*f1[i])*h[i]*h[i]
        end
        xout[k]=sum

        sum=0.
        for i=1:n
            sum += (x[k,i]*g0[i]+v[k,i]*g1[i])*h[i]  #*h[i] extra in pdf
        end
        vout[k] = sum
    end

    return xout, vout
end

function _find_db_path(geo_path)
    extra_path = joinpath(geo_path, "extra_files")
    if isdir(extra_path) && length(Glob.glob(extra_path * "/*.db*")) > 0
        return extra_path
    else
        return geo_path
    end
end

"""Create a map with 3 layers (E, N, U) of the ENU vectors for every pixel within a demrsc"""
function create_los_map(;directory=".", demrsc=nothing, outfile="los_map.h5",
                        full_elevation_dem=nothing, dbpath=nothing, dbfile=nothing)
    if isnothing(dbfile)
        isnothing(dbpath) && error("need dbfile or dbpath")
        @show dbpath
        dbfile = Glob.glob("*.db*", dbpath)[1]
    end
    @show dbfile

    if isnothing(demrsc)
        demrsc = Sario.load(Sario.find_rsc_file(directory))
    end

    param_dict = load_all_params(dbfile)
    orbinfo_filename = param_dict["orbinfo"]  # The .db file doesn't save path
    dbpath = filepath(dbfile)
    orbinfo_file = joinpath(dbpath, orbinfo_filename)
    timeorbit, xorbit, vorbit = read_orbit_vector(orbinfo_file)

    if isnothing(full_elevation_dem)
        # TODO: this is brittle
        dem_file = replace(Sario.find_rsc_file(directory=".."), ".rsc" => "")
        @show dem_file
        full_elevation_dem = Sario.load(dem_file)
    end

    xx, yy = MapImages.grid(demrsc, sparse=true)
    enu_out = Array{Float64, 3}(undef, (length(yy), length(xx), 3))
    # lats = Array{Float64, 2}(undef, (length(yy), length(xx)))
    # lons = similar(lats)

    Threads.@threads for j in 1:length(xx)
        for i in 1:length(yy)
    # for (j, x) in enumerate(xx)  # Doesn't seem to work with thread?
        # for (i, y) in enumerate(yy)
            # @show i, j, y, x
            y = yy[i]
            x = xx[j]
            xyz_los_vecs = calculate_los_xyz(y, x, full_elevation_dem, demrsc, param_dict, timeorbit, xorbit, vorbit)
            # println("$y $x is at ", MapImages.latlon_to_rowcol(demrsc, y, x))
            enu_out[i, j, :] = get_los_enu([y, x], xyz_los_vecs=xyz_los_vecs)
        end
    end

    if !isnothing(outfile)
        println("Writing to $outfile dset 'stack'")
        h5open(outfile, "w") do f
            write(f, "stack", permutedims(enu_out, (2, 1, 3)))
            write(f, "lats", collect(yy))
            write(f, "lons", collect(xx))
        end
        Sario.save_dem_to_h5(outfile, demrsc)
    end

    return enu_out
end

_check_oob(arr::AbstractArray, num) = @assert (num > minimum(arr) && num < maximum(arr)) "$num out of bounds for $(extrema(arr))"

findnearest(A::AbstractArray, t) = findmin(abs.(A .- t))[2]

function read_los_map(lat_lons::AbstractArray{<:AbstractFloat, 2}, los_map_file::String=MAP_FILENAME)
    los_map = permutedims(h5read(los_map_file, "stack"), (2, 1, 3))
    gridlats = h5read(los_map_file, "lats")
    gridlons = h5read(los_map_file, "lons")
    return read_los_map(lat_lons, los_map, gridlats, gridlons)
end

# """Given an array of lat/lon points, read the E,N,U coeffs from the LOS map"""
function read_los_map(lat_lons::AbstractArray{<:AbstractFloat, 2}, los_map::AbstractArray{<:AbstractFloat, 3}, gridlats, gridlons)
    num_points = size(lat_lons, 2)
    out = Array{eltype(lat_lons), 2}(undef, (3, num_points))
    for i = 1:num_points
        lat, lon = lat_lons[:, i]
        _check_oob(gridlats, lat)
        row = findnearest(gridlats, lat)

        _check_oob(gridlons, lon)
        col = findnearest(gridlons, lon)

        out[:, i] = los_map[row, col, :]
    end
    return out
end

read_los_map(lat_lons::AbstractArray{<:AbstractFloat, 1}, los_map_file::String=MAP_FILENAME) = read_los_map(reshape(lat_lons, :, 1), los_map_file)

end # module


text:
import Base: show

function show(io::IO, pd::PDistribution)
    x,y = map(length,states(pd))
    print(io, "ProbabilityDistribution{$(x)x$(y)}")
end

function show(io::IO, vd::DBayesianNode)
    print(io, "dnode [$(vd.index), $(vd.label)] $(vd.pd)")
end

function show(io::IO, vc::CBayesianNode)
    print(io, "cnode [$(vc.index), $(vc.label)]")
end

function show(io::IO, e::BayesianEdge)
    print(io, "$(e.source) --> $(e.target)")
end

function show(io::IO, cpd::CPD)
    dist=map(bn -> string(bn), cpd.distribution)
    j_dist=join(dist, ",")

    if length(cpd.conditionals) > 0
        cond=map(bn -> string(bn), cpd.conditionals)
        j_cond=join(cond, ",")
        print(io,"CPD{P($(j_dist)|$(j_cond))}")
    else
        print(io, "CPD{P($(j_dist))}")
    end
end

function show(io::IO, bn::BayesianNetwork)
    edges = showEdgeList(bn.edges)
    nodes = showNodeList(bn.nodes)
    print(io, "edges $edges, nodes $nodes")
end

function showNodeList(nodes::Array{BayesianNode,1})
    if length(nodes) == 0
        names = "none"
    else
        names = string("[", join(map(x -> x.index, nodes), ", "), "]")
    end
    names
end

function showEdgeList(edges::Array{BayesianEdge,1})
    if length(edges) == 0
        names = "none"
    else
        names = string("[", join(map(x -> x.index, edges), ", "), "]")
    end
    names
end