text:
"""
Collection of basic rendering functions for Cairo
"""

Cairo.move_to(ctx::CairoContext, P::VecE2) = move_to(ctx, P.x, P.y)
Cairo.line_to(ctx::CairoContext, P::VecE2) = line_to(ctx, P.x, P.y)

function grayscale_transform(color::Colorant)
    g = 1.0 - gray(convert(Gray, color)) # convert to grayscale and invert
    if isa(color, TransparentColor)
        a = alpha(color)
        GrayA(g, a)
    else
        Gray(g)
    end
end

function Cairo.set_source_rgba(ctx::CairoContext, color::RGB)

    # g = convert(Float64, gray(grayscale_transform(color)))
    # r = g = b = g

    r = convert(Float64, red(color))
    g = convert(Float64, green(color))
    b = convert(Float64, blue(color))

    set_source_rgba(ctx, r, g, b, 1.0)
end
function Cairo.set_source_rgba(ctx::CairoContext, color::TransparentColor)

    # grayscale = grayscale_transform(color)
    # g = convert(Float64, gray(grayscale))
    # r = g = b = g
    # a = convert(Float64, alpha(grayscale))

    r = convert(Float64, red(color))
    g = convert(Float64, green(color))
    b = convert(Float64, blue(color))
    a = convert(Float64, alpha(color))

    set_source_rgba(ctx, r, g, b, a)
end
function Cairo.set_source_rgba(ctx::CairoContext, color₀::Colorant, color₁::Colorant, t::Real)

    r₀ = convert(Float64, red(color₀))
    g₀ = convert(Float64, green(color₀))
    b₀ = convert(Float64, blue(color₀))
    a₀ = convert(Float64, alpha(color₀))

    r₁ = convert(Float64, red(color₁))
    g₁ = convert(Float64, green(color₁))
    b₁ = convert(Float64, blue(color₁))
    a₁ = convert(Float64, alpha(color₁))

    r = r₀ + (r₁ - r₀)*t
    g = g₀ + (g₁ - g₀)*t
    b = b₀ + (b₁ - b₀)*t
    a = a₀ + (a₁ - a₀)*t
    set_source_rgba(ctx, r, g, b, a)
end

function render_paint(
    ctx          :: CairoContext,
    color        :: Colorant,
    )

    r = convert(Float64, red(color))
    g = convert(Float64, green(color))
    b = convert(Float64, blue(color))
    a = convert(Float64, alpha(color))

    save(ctx)
    set_source_rgba(ctx, color)
    paint_with_alpha(ctx, alpha(color))
    restore(ctx)
end

function render_text(
    ctx          :: CairoContext,
    text         :: AbstractString,
    x            :: Real,
    y            :: Real,
    fontsize     :: Real,
    color        :: Colorant,
    align_center :: Bool        = false,
    fontfamily   :: AbstractString      = "monospace" # ∈ "serif", "sans-serif", "cursive", "fantasy", "monospace"
    )

    save(ctx)
    select_font_face(ctx, fontfamily, Cairo.FONT_SLANT_NORMAL, Cairo.FONT_WEIGHT_NORMAL)
    set_font_size(ctx, fontsize)
    set_source_rgba(ctx, color)

    if align_center
        extents = text_extents(ctx, text)
        x -= (extents[3]/2 + extents[1]);
        y -= (extents[4]/2 + extents[2]);
    end

    move_to(ctx, x, y)
    show_text(ctx, text)
    restore(ctx)
end

function render_circle(
    ctx          :: CairoContext,
    x            :: Real,
    y            :: Real,
    radius       :: Real,
    color_fill   :: Colorant,
    color_stroke :: Colorant = color_fill,
    line_width   :: Real = 1.0
    )

    save(ctx)
    translate(ctx, x, y)

    arc(ctx, 0.0, 0.0, radius, 0, 2pi)
    set_source_rgba(ctx, color_fill)
    fill(ctx)

    arc(ctx, 0.0, 0.0, radius, 0, 2pi)
    set_line_width(ctx, line_width)
    set_source_rgba(ctx, color_stroke)
    stroke(ctx)

    restore(ctx)
end
function render_circle(
    ctx          :: CairoContext,
    P            :: VecE2,
    radius       :: Real,
    color_fill   :: Colorant,
    color_stroke :: Colorant = color_fill,
    line_width   :: Real = 1.0
    )

    render_circle(ctx, P.x, P.y, radius, color_fill, color_stroke, line_width)
end
function render_arc(
    ctx          :: CairoContext,
    x            :: Real,
    y            :: Real,
    radius       :: Real,
    start_ang    :: Real,
    end_ang      :: Real,
    color_fill   :: Colorant,
    color_stroke :: Colorant = color_fill,
    line_width   :: Real = 1.0
    )

    line_width = user_to_device_distance!(ctx, [line_width,0])[1]

    save(ctx)
    translate(ctx, x, y)

    arc(ctx, 0.0, 0.0, radius, start_ang, end_ang)
    set_source_rgba(ctx, color_fill)
    fill(ctx)

    arc(ctx, 0.0, 0.0, radius, start_ang, end_ang)
    set_line_width(ctx, line_width)
    set_source_rgba(ctx, color_stroke)
    stroke(ctx)

    restore(ctx)
end
function render_rect(ctx::CairoContext,
                    x::Real,
                    y::Real,
                    width::Real,
                    height::Real,
                    color_fill::Colorant,
                    ffill         :: Bool        = false,
                    fstroke       :: Bool        = true,
                    color_stroke  :: Colorant    = color_fill,
                    line_width    :: Real        = 2.0
                    )
    Cairo.save(ctx)
    rectangle(ctx, x, y, width, height)
    if ffill
        set_source_rgba(ctx, color_fill)
        fstroke ? fill_preserve(ctx) : fill(ctx)
    end

    if fstroke
        set_source_rgba(ctx, color_stroke)
        set_line_width(ctx, line_width)
        stroke(ctx)
    end

    restore(ctx)
end
function render_round_rect(
    ctx           :: CairoContext,
    x             :: Real,
    y             :: Real,
    width         :: Real,
    height        :: Real,
    aspect        :: Real,
    corner_radius :: Real,
    color_fill    :: Colorant,
    ffill         :: Bool        = false,
    fstroke       :: Bool        = true,
    color_stroke  :: Colorant    = color_fill,
    line_width    :: Real        = 2.0
    )
    # (x,y) are the center of the rectangle

    save(ctx)
    radius = corner_radius / aspect

    d2r = pi/180

    new_sub_path(ctx)
    arc(ctx, x + width/2 - radius, y - height/2 + radius, radius, -90 * d2r,   0 * d2r)
    arc(ctx, x + width/2 - radius, y + height/2 - radius, radius,   0 * d2r,  90 * d2r)
    arc(ctx, x - width/2 + radius, y + height/2 - radius, radius,  90 * d2r, 180 * d2r)
    arc(ctx, x - width/2 + radius, y - height/2 + radius, radius, 180 * d2r, 270 * d2r)
    close_path(ctx)

    if ffill
        set_source_rgba(ctx, color_fill)
        fstroke ? fill_preserve(ctx) : fill(ctx)
    end

    if fstroke
        set_source_rgba(ctx, color_stroke)
        # TODO: bug - line_width seems to be changing based on vehicle heading angle
        # from Cairo.user_to_device_distance! - this may be the bug
        line_width = abs(user_to_device_distance!(ctx, [line_width,0.0])[1])
        set_line_width(ctx, line_width)
        stroke(ctx)
    end

    restore(ctx)
end
function render_car(
    ctx           :: CairoContext,
    x             :: Real, # center of vehicle [m]
    y             :: Real, # center of vehicle [m]
    yaw           :: Real, # counterclockwise [rad]
    color_fill    :: Colorant,
    color_stroke  :: Colorant = color_fill;

    color_arrow   :: Colorant = RGB(1.0,1.0,1.0),
    car_length    :: Float64 = 4.6, # [m]
    car_width     :: Float64 = 2.0, # [m]
    corner_radius :: Float64 = 0.3, # [m]
    corner_aspect :: Float64 = 1.0, # [m]
    arrow_width   :: Float64 = 0.3, # [% car width]
    arrow_length  :: Float64 = 0.6, # [% car length]
    arrow_headlen :: Float64 = 0.5, # [% arrow len that head takes up]
    line_width    :: Float64 = 0.1, # [m]
    )

    # renders a car (rounded rectangle w/ arrow) at the given location

    save(ctx)

    translate(ctx, x, y)
    rotate(ctx, yaw)

    render_round_rect(ctx, 0, 0, car_length, car_width, corner_aspect, corner_radius, color_fill, true, true, color_stroke, line_width)

    # render the arrow
    wid = car_width*arrow_width
    len = car_length*arrow_length
    hed = len*arrow_headlen

    new_sub_path(ctx)
    move_to(ctx, 0+len/2,     0      )
    line_to(ctx, 0+len/2-hed, 0+hed/2)
    line_to(ctx, 0+len/2-hed, 0+wid/2)
    line_to(ctx, 0-len/2,     0+wid/2)
    line_to(ctx, 0-len/2,     0-wid/2)
    line_to(ctx, 0+len/2-hed, 0-wid/2)
    line_to(ctx, 0+len/2-hed, 0-hed/2)
    close_path(ctx)

    set_source_rgba(ctx, color_arrow)
    fill(ctx)

    restore(ctx)
end
function render_vehicle(
    ctx           :: CairoContext,
    x             :: Real, # x-pos of the center of the vehicle
    y             :: Real, # y-pos of the center of the vehicle
    yaw           :: Real, # heading angle [rad]
    length        :: Real, # vehicle length
    width         :: Real, # vehicle width
    color_fill    :: Colorant,
    color_stroke  :: Colorant = color_fill,

    color_arrow   :: Colorant = RGB(1.0,1.0,1.0),
    corner_radius :: Float64 = 0.5,
    corner_aspect :: Float64 = 1.0,
    arrow_width   :: Float64 = 0.3, # [% car width]
    arrow_length  :: Float64 = 0.6, # [% car length]
    arrow_headlen :: Float64 = 0.5, # [% arrow len that head takes up]
    line_width    :: Float64 = 0.3,
    )

    # renders a car (rounded rectangle w/ arrow) at the given location
    # (x,y) are in meters and yaw is the radians, counter-clockwise from pos x axis

    save(ctx)

    # translate(ctx, x - 0.5length*cos(yaw), y - 0.5length*sin(yaw))
    translate(ctx, x, y)
    rotate(ctx, yaw)

    render_round_rect(ctx, 0, 0, length, width, corner_aspect, corner_radius, color_fill, true, true, color_stroke, line_width)

    # render the arrow
    wid = width*arrow_width
    len = length*arrow_length
    hed = min(len*arrow_headlen, width)

    new_sub_path(ctx)
    move_to(ctx, 0+len/2,     0      )
    line_to(ctx, 0+len/2-hed, 0+hed/2)
    line_to(ctx, 0+len/2-hed, 0+wid/2)
    line_to(ctx, 0-len/2,     0+wid/2)
    line_to(ctx, 0-len/2,     0-wid/2)
    line_to(ctx, 0+len/2-hed, 0-wid/2)
    line_to(ctx, 0+len/2-hed, 0-hed/2)
    close_path(ctx)

    set_source_rgba(ctx, color_arrow)
    fill(ctx)

    restore(ctx)
end

# aggregate
function render_point_trail(
    ctx           :: CairoContext,
    pts           :: AbstractVector{T}, # 2×n
    color         :: Colorant,
    circle_radius :: Real = 0.25 ) where {T<:Real}

    save(ctx)

    for i = 1 : size(pts,2)
        arc(ctx, pts[1,i], pts[2,i], circle_radius, 0, 2pi)
        set_line_width(ctx, 1.0)
        set_source_rgba(ctx, color)
        stroke(ctx)
    end

    restore(ctx)
end
function render_line(
    ctx        :: CairoContext,
    pts        :: AbstractArray{T}, # 2×n
    color      :: Colorant,
    line_width :: Real = 1.0,
    line_cap   :: Integer=Cairo.CAIRO_LINE_CAP_ROUND, # CAIRO_LINE_CAP_BUTT, CAIRO_LINE_CAP_ROUND, CAIRO_LINE_CAP_SQUARE
    ) where {T<:Real}

    line_width = user_to_device_distance!(ctx, [line_width,0])[1]

    save(ctx)
    set_source_rgba(ctx,color)
    set_line_width(ctx,line_width)
    set_line_cap(ctx, line_cap)

    move_to(ctx, pts[1,1], pts[2,1])
    for i = 2 : size(pts,2)
        line_to(ctx, pts[1,i], pts[2,i])
    end
    stroke(ctx)
    restore(ctx)
end
function render_line(
    ctx        :: CairoContext,
    pts        :: AbstractVector{VecE2{T}}, # 2×n
    color      :: Colorant,
    line_width :: Real = 1.0,
    line_cap   :: Integer=Cairo.CAIRO_LINE_CAP_ROUND, # CAIRO_LINE_CAP_BUTT, CAIRO_LINE_CAP_ROUND, CAIRO_LINE_CAP_SQUARE
    ) where {T<:Real}

    line_width = user_to_device_distance!(ctx, [line_width,0])[1]

    save(ctx)
    set_source_rgba(ctx,color)
    set_line_width(ctx,line_width)
    set_line_cap(ctx, line_cap)

    move_to(ctx, pts[1])
    for i in 2 : length(pts)
        line_to(ctx, pts[i])
    end
    stroke(ctx)
    restore(ctx)
end
function render_closed_line(
    ctx        :: CairoContext,
    pts        :: AbstractArray{T}, # 2×n
    color      :: Colorant,
    line_width :: Real = 1.0,
    fill_color :: Colorant = RGBA(0.0,0.0,0.0,0.0),
    ) where {T<:Real}

    line_width = user_to_device_distance!(ctx, [line_width,0])[1]

    save(ctx)

    set_source_rgba(ctx, color)
    set_line_width(ctx,line_width)

    move_to(ctx, pts[1,1], pts[2,1])
    for i in 2 : size(pts,2)
        line_to(ctx, pts[1,i], pts[2,i])
    end
    close_path(ctx)

    if alpha(fill_color) > 0.0
        stroke_preserve(ctx)
        set_source_rgba(ctx, fill_color)
        fill(ctx)
    else
        stroke(ctx)
    end

    restore(ctx)
end
function render_closed_line(
    ctx        :: CairoContext,
    pts        :: AbstractVector{VecE2{T}},
    color      :: Colorant,
    line_width :: Real = 1.0,
    fill_color :: Colorant = RGBA(0.0,0.0,0.0,0.0),
    ) where {T<:Real}

    line_width = user_to_device_distance!(ctx, [line_width,0])[1]

    save(ctx)

    set_source_rgba(ctx, color)
    set_line_width(ctx,line_width)

    move_to(ctx, pts[1])
    for i in 2 : length(pts)
        line_to(ctx, pts[i])
    end
    close_path(ctx)

    if alpha(fill_color) > 0.0
        stroke_preserve(ctx)
        set_source_rgba(ctx, fill_color)
        fill(ctx)
    else
        stroke(ctx)
    end

    restore(ctx)
end
function render_fill_region(
    ctx        :: CairoContext,
    pts        :: AbstractArray{T}, # 2×n
    color      :: Colorant,
    ) where {T<:Real}

    save(ctx)
    set_source_rgba(ctx,color)
    move_to(ctx, pts[1,1], pts[2,1])
    for i = 2 : size(pts,2)
        line_to(ctx, pts[1,i], pts[2,i])
    end
    close_path(ctx)
    fill(ctx)
    restore(ctx)

    save(ctx)
    set_source_rgba(ctx,color)
    set_line_width(ctx,1.0)
    move_to(ctx, pts[1,1], pts[2,1])
    for i = 2 : size(pts,2)
        line_to(ctx, pts[1,i], pts[2,i])
    end
    close_path(ctx)
    stroke(ctx)
    restore(ctx)
end
function render_line_segment(
    ctx        :: CairoContext,
    x1         :: Float64,
    y1         :: Float64,
    x2         :: Float64,
    y2         :: Float64,
    color      :: Colorant,
    line_width :: Real = 1.0
    )

    line_width = user_to_device_distance!(ctx, [line_width,0])[1]

    save(ctx)
    set_source_rgba(ctx,color)
    set_line_width(ctx,line_width)
    set_line_cap(ctx, Cairo.CAIRO_LINE_CAP_ROUND)

    move_to(ctx, x1, y1)
    line_to(ctx, x2, y2)
    stroke(ctx)
    restore(ctx)
end
function render_line_segment(
    ctx        :: CairoContext,
    A          :: VecE2,
    B          :: VecE2,
    color      :: Colorant,
    line_width :: Real = 1.0
    )

    line_width = user_to_device_distance!(ctx, [line_width,0])[1]

    save(ctx)
    set_source_rgba(ctx,color)
    set_line_width(ctx,line_width)
    set_line_cap(ctx, Cairo.CAIRO_LINE_CAP_ROUND)

    move_to(ctx, A)
    line_to(ctx, B)
    stroke(ctx)
    restore(ctx)
end
function render_dashed_line(
    ctx          :: CairoContext,
    pts          :: AbstractArray{T}, # 2×n
    color        :: Colorant,
    line_width_in   :: Real = 1.0,
    dash_length_in  :: Real = 1.0,
    dash_spacing_in :: Real = 1.0,
    dash_offset_in  :: Real = 0.0,
    line_cap :: Integer = Cairo.CAIRO_LINE_CAP_SQUARE ,
    ) where {T<:Real}

    line_width   = user_to_device_distance!(ctx, [line_width_in,  0])[1]
    dash_length  = user_to_device_distance!(ctx, [dash_length_in, 0])[1]
    dash_spacing = user_to_device_distance!(ctx, [dash_spacing_in,0])[1]
    dash_offset  = user_to_device_distance!(ctx, [dash_offset_in, 0])[1]
    dashes = [dash_length, dash_spacing]

    save(ctx)
    set_source_rgba(ctx,color)
    set_line_width(ctx,line_width)
    set_line_cap(ctx, convert(Int32, line_cap))
    set_dash(ctx, dashes, dash_offset)

    move_to(ctx, pts[1,1], pts[2,1])
    for i = 2 : size(pts,2)
        line_to(ctx, pts[1,i], pts[2,i])
    end
    stroke(ctx)
    restore(ctx)
end
function render_dashed_arc(
    ctx          :: CairoContext,
    x            :: Real,
    y            :: Real,
    radius       :: Real,
    start_ang    :: Real,
    end_ang      :: Real,
    color_fill   :: Colorant,
    color_stroke :: Colorant   = color_fill,
    line_width_in   :: Real = 1.0,
    dash_length_in  :: Real = 1.0,
    dash_spacing_in :: Real = 1.0,
    dash_offset_in  :: Real = 0.0
    )

    line_width   = user_to_device_distance!(ctx, [line_width_in,  0])[1]
    dash_length  = user_to_device_distance!(ctx, [dash_length_in, 0])[1]
    dash_spacing = user_to_device_distance!(ctx, [dash_spacing_in,0])[1]
    dash_offset  = user_to_device_distance!(ctx, [dash_offset_in, 0])[1]
    dashes = [dash_length, dash_spacing]

    save(ctx)
    translate(ctx, x, y)

    arc(ctx, 0.0, 0.0, radius, start_ang, end_ang)
    set_source_rgba(ctx, color_fill)
    fill(ctx)

    arc(ctx, 0.0, 0.0, radius, start_ang, end_ang)
    set_source_rgba(ctx, color_stroke)
    set_line_width(ctx, line_width)
    set_line_cap(ctx, Cairo.CAIRO_LINE_CAP_ROUND)
    set_dash(ctx, dashes, dash_offset)
    stroke(ctx)

    restore(ctx)
end
function render_arrow(
    ctx               :: CairoContext,
    pts               :: AbstractArray{T}, # 2×n
    color             :: Colorant,
    line_width        :: Real,
    arrowhead_len     :: Real;
    ARROW_WIDTH_RATIO :: AbstractFloat = 0.8,
    ARROW_ALPHA       :: AbstractFloat = 0.1pi,
    ARROW_BETA        :: AbstractFloat = 0.8pi
    ) where {T<:Real}

    @assert(size(pts,2) > 1)

    line_width = user_to_device_distance!(ctx, [line_width,0])[1]

    save(ctx)
    set_source_rgba(ctx, color)
    set_line_width(ctx, line_width)
    set_line_cap(ctx, Cairo.CAIRO_LINE_CAP_ROUND)

    move_to(ctx, pts[1,1], pts[2,1])
    for i = 2 : size(pts,2)
        line_to(ctx, pts[1,i], pts[2,i])
    end
    stroke(ctx)

    # orientation of the arrowhead
    theta = atan(pts[2,end]-pts[2,end-1], pts[1,end]-pts[1,end-1])

    arrowhead_width = arrowhead_len * ARROW_WIDTH_RATIO
    whatev =    pi - ARROW_BETA
    delta = 0.5pi - whatev
    epsilon =  pi - ARROW_BETA - ARROW_ALPHA
    len_prime = 0.5arrowhead_width*sin(delta + epsilon)/sin(ARROW_ALPHA) - arrowhead_len

    R = [ cos(theta) -sin(theta);
          sin(theta)  cos(theta) ]

    # render the arrowhead
    O = pts[:,end] # origin, center of the arrowhead
    A = O + R * [ arrowhead_len/2 0]'
    B = O + R * [-arrowhead_len/2 0]'
    C = O + R * [-arrowhead_len/2-len_prime  arrowhead_width/2]'
    D = O + R * [-arrowhead_len/2-len_prime -arrowhead_width/2]'

    new_sub_path(ctx)
    move_to(ctx, A[1], A[2])
    line_to(ctx, C[1], C[2])
    line_to(ctx, B[1], B[2])
    line_to(ctx, D[1], D[2])
    close_path(ctx)
    fill(ctx)

    restore(ctx)
end

function render_colormesh(
    ctx::CairoContext,
    C::AbstractArray{T}, # n×m matrix of 0->1 values
    X::AbstractVector{S}, # n+1 vector of x bin boundaries
    Y::AbstractVector{U} # m+1 vector of y bin boundaries
    ) where {T<:Real, S<:Real, U<:Real}

    n,m = size(C)
    @assert(length(X) == n+1)
    @assert(length(Y) == n+1)

    save(ctx)

    x₂ = X[1]
    for i = 1 : n
        x₁, x₂ = x₂, X[i+1]
        y₂ = Y[1]
        for j = 1 : m
            y₁, y₂ = y₂, Y[j+1]

            c = C[i,j]

            new_sub_path(ctx)
            move_to(ctx, x₁, y₁)
            line_to(ctx, x₂, y₁)
            line_to(ctx, x₂, y₂)
            line_to(ctx, x₁, y₂)
            close_path(ctx)

            set_source_rgba(ctx, c, c, c, 1.0)
            fill(ctx)
        end
    end

    restore(ctx)
end
function render_colormesh(
    ctx::CairoContext,
    C::AbstractArray{T}, # n×m matrix of 0->1 values
    X::AbstractVector{S}, # n+1 vector of x bin boundaries
    Y::AbstractVector{U}, # m+1 vector of y bin boundaries
    color₀::Colorant, # color for c = 0
    color₁::Colorant, # color for c = 1
    ) where {T<:Real, S<:Real, U<:Real}

    n,m = size(C)
    @assert(length(X) == n+1)
    @assert(length(Y) == m+1)

    save(ctx)

    x₂ = X[1]
    for i = 1 : n
        x₁, x₂ = x₂, X[i+1]
        y₂ = Y[1]
        for j = 1 : m
            y₁, y₂ = y₂, Y[j+1]

            set_source_rgba(ctx, color₀, color₁, C[i,j])

            new_sub_path(ctx)
            move_to(ctx, x₁, y₁)
            line_to(ctx, x₂, y₁)
            line_to(ctx, x₂, y₂)
            line_to(ctx, x₁, y₂)
            close_path(ctx)

            fill(ctx)
        end
    end

    restore(ctx)
end


text:
# ========== Single Step Autocorrelation Wavelet Transform ==========
## 1D ##
"""
    acdwt_step(v, d, h, g)

Performs one level of the autocorrelation discrete wavelet transform (acdwt) on the vector
`v`, which is the j-th level scaling coefficients (Note the 0th level scaling coefficients
is the raw signal). The vectors `h` and `g` are the detail and scaling filters.

# Arguments
- `v::AbstractArray{T} where T<:Number`: Array of coefficients from a node at level `d`.
- `d::Integer`: Depth level of `v`.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
- `w₁::AbstractVector{T} where T<:Number` or `w₁::AbstractMatrix{T} where T<:Number`: Vector
  allocation for output from low pass filter (1D case); or matrix allocation for output from
  low + low pass filter (2D case).
- `w₂::AbstractVector{T} where T<:Number` or `w₂::AbstractMatrix{T} where T<:Number`: Vector
  allocation for output from high pass filter (1D case); or matrix allocation for output
  from low + high pass filter (2D case).
- `w₃::AbstractVector{T} where T<:Number` or `w₃::AbstractMatrix{T} where T<:Number`: Matrix
  allocation for output from high + low pass filter (2D case).
- `w₄::AbstractVector{T} where T<:Number` or `w₄::AbstractMatrix{T} where T<:Number`: Matrix
  allocation for output from high + high pass filter (2D case).

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.make_acreverseqmfpair(wt)

# One step of ACDWT
ACWT.acdwt_step(v, 0, h, g)
```

**See also:** [`acdwt`](@ref), [`iacdwt`](@ref)
"""
function acdwt_step(v::AbstractVector{T}, d::Integer, h::Array{S,1}, g::Array{S,1}) where
                   {T<:Number, S<:Number}
    n = length(v)
    w₁ = zeros(T,n)
    w₂ = zeros(T,n)

    acdwt_step!(w₁, w₂, v, d, h, g)
    return w₁, w₂
end

"""
    acdwt_step!(w₁, w₂, v, j, h, g)

Same with `acdwt_step` but without array allocation.

# Arguments
- `w₁::AbstractVector{T} where T<:Number` or `w₁::AbstractMatrix{T} where T<:Number`: Vector
  allocation for output from low pass filter (1D case); or matrix allocation for output from
  low + low pass filter (2D case).
- `w₂::AbstractVector{T} where T<:Number` or `w₂::AbstractMatrix{T} where T<:Number`: Vector
  allocation for output from high pass filter (1D case); or matrix allocation for output
  from low + high pass filter (2D case).
- `w₃::AbstractVector{T} where T<:Number` or `w₃::AbstractMatrix{T} where T<:Number`: Matrix
  allocation for output from high + low pass filter (2D case).
- `w₄::AbstractVector{T} where T<:Number` or `w₄::AbstractMatrix{T} where T<:Number`: Matrix
  allocation for output from high + high pass filter (2D case).
- `v::AbstractArray{T} where T<:Number`: Array of coefficients from a node at level `d`.
- `d::Integer`: Depth level of `v`.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
- `w₁::AbstractVector{T} where T<:Number` or `w₁::AbstractMatrix{T} where T<:Number`: Output
  from low pass filter (1D case); or output from low + low pass filter (2D case).
- `w₂::AbstractVector{T} where T<:Number` or `w₂::AbstractMatrix{T} where T<:Number`: Output
  from high pass filter (1D case); or output from low + high pass filter (2D case).
- `w₃::AbstractVector{T} where T<:Number` or `w₃::AbstractMatrix{T} where T<:Number`: Output
  from high + low pass filter (2D case).
- `w₄::AbstractVector{T} where T<:Number` or `w₄::AbstractMatrix{T} where T<:Number`: Output
  from high + high pass filter (2D case).

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
w₁ = similar(v)
w₂ = similar(v)
wt = wavelet(WT.haar)
g, h = WT.make_acreverseqmfpair(wt)

# One step of ACDWT
ACWT.acdwt_step!(w₁, w₂, v, 0, h, g)
```

**See also:** [`acdwt_step!`](@ref), [`acdwt`](@ref), [`iacdwt`](@ref)
"""
function acdwt_step!(w₁::AbstractVector{T},
                     w₂::AbstractVector{T},
                     v::AbstractVector{T},
                     d::Integer,
                     h::Array{T,1},
                     g::Array{T,1}) where {T<:Number, S<:Number}
    # Sanity check
    @assert length(w₁) == length(w₂) == length(v)
    @assert length(h) == length(g)

    # Setup
    N = length(v)
    L = length(h)

    # One step of autocorrelation transform
    for i in 1:N
        t = i+(1<<d) |> t -> t>N ? mod1(t,N) : t
        i = mod1(i + (L÷2+1) * 2^d,N) # Need to shift by half the filter size because of periodicity assumption
        @inbounds w₁[i] = g[1] * v[t]
        @inbounds w₂[i] = h[1] * v[t]
        for n in 2:L
            t = t+(1<<d) |> t -> t>N ? mod1(t,N) : t
            @inbounds w₁[i] += g[n] * v[t]
            @inbounds w₂[i] += h[n] * v[t]
        end
    end
    return w₁, w₂
end

"""
    iacdwt_step(w₁, w₂)
    iacdwt_step(w₁, w₂, w₃, w₄)

Perform one level of the inverse autocorrelation discrete wavelet transform (IACDWT) on the
vectors `w₁` and `w₂`, which are the `j+1`-th level scaling coefficients (Note that the 0th
level scaling coefficients is the raw signal).

# Arguments
- `w₁::AbstractVector{T} where T<:Number` or `w₁::AbstractMatrix{T} where T<:Number`:
  Coefficients of left child node (1D case); or coefficients from top left child node (2D
  case).
- `w₂::AbstractVector{T} where T<:Number` or `w₂::AbstractMatrix{T} where T<:Number`:
  Coefficients of right child node (1D case); or coefficients from top right child node (2D
  case).
- `w₃::AbstractVector{T} where T<:Number` or `w₃::AbstractMatrix{T} where T<:Number`:
  Coefficients from bottom left child node (2D case).
- `w₄::AbstractVector{T} where T<:Number` or `w₄::AbstractMatrix{T} where T<:Number`:
  Coefficients from bottom right child node (2D case).

# Returns
- `v::Array{T}`: Reconstructed coefficients.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.make_acreverseqmfpair(wt)

# One step of ACDWT
w₁, w₂ = ACWT.acdwt_step(v, 0, h, g)

# One step of IACDWT
v̂ = ACWT.iacdwt_step(w₁, w₂)
```

**See also:** [`iacdwt_step!`](@ref), [`acdwt_step`](@ref), [`iacdwt`](@ref)
"""
function iacdwt_step(w₁::AbstractVector{T}, w₂::AbstractVector{T}) where T<:Number
    v = similar(w₁)
    iacdwt_step!(v, w₁, w₂)
    return v
end

"""
    iacdwt_step!(v, w₁, w₂)

Same as `iacdwt_step` but without array allocation.

# Arguments
- `v::AbstractArray{T} where T<:Number`: Array allocation for reconstructed coefficients.
- `w₁::AbstractVector{T} where T<:Number` or `w₁::AbstractMatrix{T} where T<:Number`:
  Coefficients of left child node (1D case); or coefficients from top left child node (2D
  case).
- `w₂::AbstractVector{T} where T<:Number` or `w₂::AbstractMatrix{T} where T<:Number`:
  Coefficients of right child node (1D case); or coefficients from top right child node (2D
  case).
- `w₃::AbstractVector{T} where T<:Number` or `w₃::AbstractMatrix{T} where T<:Number`:
  Coefficients from bottom left child node (2D case).
- `w₄::AbstractVector{T} where T<:Number` or `w₄::AbstractMatrix{T} where T<:Number`:
  Coefficients from bottom right child node (2D case).

# Returns
- `v::Array{T}`: Reconstructed coefficients.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.make_acreverseqmfpair(wt)

# One step of ACDWT
w₁, w₂ = ACWT.acdwt_step(v, 0, h, g)

# One step of IACDWT
v̂ = similar(v)
ACWT.iacdwt_step!(v̂, w₁, w₂)
```

**See also:** [`iacdwt_step`](@ref), [`acdwt_step`](@ref), [`iacdwt`](@ref)
"""
function iacdwt_step!(v::AbstractVector{T},
                      w₁::AbstractVector{T},
                      w₂::AbstractVector{T}) where T<:Number
    @assert length(v) == length(w₁) == length(w₂)
    for i in eachindex(v)
        @inbounds v[i] = (w₁[i]+w₂[i]) / √2
    end
end

## 2D ##
# Forward transform with allocation
function acdwt_step(v::AbstractMatrix{T}, d::Integer, h::Vector{S}, g::Vector{S}) where
                   {T<:Number, S<:Number}
    n, m = size(v)
    w₁ = Matrix{T}(undef, (n,m))
    w₂ = Matrix{T}(undef, (n,m))
    w₃ = Matrix{T}(undef, (n,m))
    w₄ = Matrix{T}(undef, (n,m))
    temp = Array{T,3}(undef, (n,m,2))
    acdwt_step!(w₁, w₂, w₃, w₄, v, d, h, g, temp)
    return w₁, w₂, w₃, w₄
end
# Forward transform without allocation
function acdwt_step!(w₁::AbstractMatrix{T}, w₂::AbstractMatrix{T},
                     w₃::AbstractMatrix{T}, w₄::AbstractMatrix{T},
                     v::AbstractMatrix{T},
                     d::Integer,
                     h::Vector{S}, g::Vector{S},
                     temp::AbstractArray{T,3}) where {T<:Number, S<:Number}
    # Sanity check
    @assert size(v) == size(w₁) == size(w₂) == size(w₃) == size(w₄)
    @assert ndims(temp) == 3
    @assert size(temp,1) == size(v,1)
    @assert size(temp,2) == size(v,2)

    # Setup
    n, m = size(w₁)

    # --- Transform ---
    # Compute acdwt for all columns
    for j in 1:m
        @inbounds temp₁ⱼ = @view temp[:,j,1]
        @inbounds temp₂ⱼ = @view temp[:,j,2]
        @inbounds vⱼ = @view v[:,j]
        @inbounds acdwt_step!(temp₁ⱼ, temp₂ⱼ, vⱼ, d, h, g)
    end
    # Compute acdwt for all rows
    for i in 1:n
        @inbounds temp₁ᵢ = @view temp[i,:,1]
        @inbounds w₁ᵢ = @view w₁[i,:]
        @inbounds w₂ᵢ = @view w₂[i,:]
        @inbounds acdwt_step!(w₁ᵢ, w₂ᵢ, temp₁ᵢ, d, h, g)

        @inbounds temp₂ᵢ = @view temp[i,:,2]
        @inbounds w₃ᵢ = @view w₃[i,:]
        @inbounds w₄ᵢ = @view w₄[i,:]
        @inbounds acdwt_step!(w₃ᵢ, w₄ᵢ, temp₂ᵢ, d, h, g)
    end
    return w₁, w₂, w₃, w₄
end

# Inverse transform with allocation
function iacdwt_step(w₁::AbstractMatrix{T}, w₂::AbstractMatrix{T},
                     w₃::AbstractMatrix{T}, w₄::AbstractMatrix{T}) where T<:Number
    n, m = size(w₁)
    v = Matrix{T}(undef, (n,m))
    temp = Array{T,3}(undef, (n,m,2))
    iacdwt_step!(v, w₁, w₂, w₃, w₄, temp)
    return v
end
# Inverse transform without allocation
function iacdwt_step!(v::AbstractMatrix{T},
                      w₁::AbstractMatrix{T}, w₂::AbstractMatrix{T},
                      w₃::AbstractMatrix{T}, w₄::AbstractMatrix{T},
                      temp::AbstractArray{T,3}) where T<:Number
    # Sanity check
    @assert size(v) == size(w₁) == size(w₂) == size(w₃) == size(w₄)
    @assert ndims(temp) == 3
    @assert size(temp,1) == size(v,1)
    @assert size(temp,2) == size(v,2)

    # Setup
    n, m = size(w₁)

    # --- Transform ---
    # Compute iacdwt for all rows
    for i in 1:n
        @inbounds temp₁ᵢ = @view temp[i,:,1]
        @inbounds w₁ᵢ = @view w₁[i,:]
        @inbounds w₂ᵢ = @view w₂[i,:]
        @inbounds iacdwt_step!(temp₁ᵢ, w₁ᵢ, w₂ᵢ)

        @inbounds temp₂ᵢ = @view temp[i,:,2]
        @inbounds w₃ᵢ = @view w₃[i,:]
        @inbounds w₄ᵢ = @view w₄[i,:]
        @inbounds iacdwt_step!(temp₂ᵢ, w₃ᵢ, w₄ᵢ)
    end
    # Compute iacdwt for all columns
    for j in 1:m
        @inbounds temp₁ⱼ = @view temp[:,j,1]
        @inbounds temp₂ⱼ = @view temp[:,j,2]
        @inbounds vⱼ = @view v[:,j]
        @inbounds iacdwt_step!(vⱼ, temp₁ⱼ, temp₂ⱼ)
    end
    return v
end

ext: jl


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))
res = ŷ .- y
stem(res)

savefig(joinpath(@OUTPUT, "ISL-lab-3-res.svg")) # hide


text:
extension(f::AbstractString) = splitext(f)[2]

######################
## WORD BREAKING

struct SolidString
    c::Array{Char, 1}
    function SolidString(s::AbstractString, l::Integer)
        l > length(s) && error("Offset larger than string length.")

        c = Char[]
        i = firstindex(s)
        count = 0
        while !(i > ncodeunits(s)) && count < l
            push!(c, s[i])
            (k, i) = iterate(s, i)
            count +=1
        end
        new(c)
    end
    SolidString(s::AbstractString) = SolidString(s, length(s))
end

import Base: length, getindex
length(s::SolidString) = length(s.c)
endof(s::SolidString) = length(s)
getindex(s::SolidString, i::Integer) = s.c[i]
getindex(s::SolidString, i::UnitRange) = string(s.c[i]...)

# maybe not the most efficient way of doing this.
global const _word_bounardy = [
    ' ', '\n','\t','(',')','[',']',',','\'',
    '*','+','-','/','\\','%','{','}','#',':',
    '&','|','?','!','"','$','=','>','<'
]
global const _word_bounardy_dot = [_word_bounardy; '.']#include dot in function of the context

function is_word_boundary(s::Char, include_dot::Bool)
    w = include_dot ? _word_bounardy_dot : _word_bounardy
    for c in w
        s == c && return true
    end
    false
end

function extend_word_backward(it::Integer, txt::SolidString, include_dot::Bool)
    it <= 1 && return 1

    while !is_word_boundary(txt[it], include_dot)
        it == 1 && return it
        it = it-1
    end
    return it+1 #I stopped at the boundary
end
function extend_word_forward(it::Integer, txt::SolidString, include_dot::Bool)
    it >= length(txt) && return length(txt)

    while !is_word_boundary(txt[it], include_dot)
        it == length(txt) && return it
        it = it+1
    end
    return it-1 #I stopped at the boundary
end

function extend_word(txt, pos, include_dot)
    stxt = SolidString(txt)#this is a bit of a mess
    i = extend_word_backward(pos, stxt, include_dot)
    j = extend_word_forward(pos, stxt, include_dot)

    if j < length(txt) && stxt[j+1] == '!' #allow for a single ! at the end of words
       j = j + 1
    end
    stxt[i:j], i, j
end

function select_word(it::GtkTextIter, buffer::GtkTextBuffer, include_dot::Bool)#include_dot means we include "." in word boundary def

    (txt, line_start, line_end) = get_line_text(buffer, it)

    pos = offset(it) - offset(line_start) +1#not sure about the +1 but it feels better
    if pos <= 0
#        warn("negative position $pos ($(offset(it)) - $(offset(line_start)) )")
        return ("", GtkTextIter(buffer, offset(it)),
        GtkTextIter(buffer, offset(it)))
    end

    word, i, j = extend_word(txt, pos, include_dot)

    its = GtkTextIter(buffer, i + offset(line_start) )
    ite = GtkTextIter(buffer, j + offset(line_start) + 1)

    return (word, its, ite)
end
select_word(it::GtkTextIter, buffer::GtkTextBuffer) = select_word(it, buffer, true)

function select_word_backward(it::GtkTextIter, buffer::GtkTextBuffer, include_dot::Bool)

    (txt, line_start, line_end) = get_line_text(buffer, it)
    pos = offset(it) - offset(line_start) #position of cursor in txt

    if pos <= 0 || length(txt) == 0
#        warn("negative positon $pos ($(offset(it)) - $(offset(line_start)))")
        return ("", GtkTextIter(buffer, offset(it)),
        GtkTextIter(buffer, offset(it)))
    end

    txt = SolidString(txt, pos)
    (i, j) = select_word_backward(pos, txt, include_dot)

    its = GtkTextIter(buffer, i + offset(line_start) )
    ite = GtkTextIter(buffer, offset(it))

    return (txt[i:j], its, it)
end
function select_word_backward(pos::Integer, txt::SolidString, include_dot::Bool)

    j = pos
    #allow for autocomplete on functions and arrays
    pos = txt[pos] == '(' ? pos-1 : pos
    pos = txt[pos] == '!' ? pos-1 : pos
    pos = txt[pos] == '[' ? pos-1 : pos

    i = extend_word_backward(pos, txt, include_dot)

    #allow for \alpha and such
    i = (i > 1 && txt[i-1] == '\\') ? i-1 : i

    return (i, j)
end
select_word_backward(pos::Integer, txt::AbstractString, include_dot::Bool) =
select_word_backward(pos, SolidString(txt, pos), include_dot)

######################
## Utility functions

function select_tuple(it::GtkTextIter, buffer::GtkTextBuffer)

    (txt, line_start, line_end) = get_line_text(buffer, it)
    pos = offset(it) - offset(line_start) #position of cursor in txt

    if pos <= 1 || length(txt) < 2 || pos > length(txt)
        return (false, nothing, nothing)
    end
    txt = txt[1:pos]

    i = findlast(isequal('('), txt)
    isnothing(i) && return (false, nothing, nothing)

    its = GtkTextIter(buffer, i + offset(line_start))
    return (true, txt[i:pos], its)
end

function text_iter_line_start(it::GtkTextIter)
    b = Gtk.buffer(it)
    (txt, line_start, line_end) = get_line_text(b, it)
    i = lstrip_idx(txt)
    i > length(txt) && return it

    return GtkTextIter(b, offset(line_start) + i)
end

function lstrip_idx(s::AbstractString)
    j = firstindex(s)
    for (i, c) in pairs(s)
        !isspace(c) && return j
        j = i
    end
    firstindex(s)
end

get_buffer(view::GtkTextView) = get_gtk_property(view, :buffer, GtkTextBuffer)

function select_on_ctrl_shift(direction, buffer::GtkSourceBuffer)

    (found, its, ite) = selection_bounds(buffer)

    if direction == :start
        ite, its = its, ite
    end

    its = found ? nonmutable(buffer, its) : get_text_iter_at_cursor(buffer)

    direction == :start && move_cursor_to_sentence_start(buffer)
    direction == :end && move_cursor_to_sentence_end(buffer)

    ite = get_text_iter_at_cursor(buffer)
    select_range(buffer, ite, its)#invert here so the cursor end up on the far right
end


text:
"""
    fk5hip()

FK5 to Hipparcos rotation and spin.

### Returned ###

- `r5h`: R-matrix: FK5 rotation wrt Hipparcos (Note 2)
- `s5h`: R-vector: FK5 spin wrt Hipparcos (Note 3)

### Notes ###

1. This function models the FK5 to Hipparcos transformation as a
   pure rotation and spin;  zonal errors in the FK5 catalogue are
   not taken into account.

2. The r-matrix r5h operates in the sense:

         P_Hipparcos = r5h x P_FK5

   where P_FK5 is a p-vector in the FK5 frame, and P_Hipparcos is
   the equivalent Hipparcos p-vector.

3. The r-vector s5h represents the time derivative of the FK5 to
   Hipparcos rotation.  The units are radians per year (Julian,
   TDB).

### Called ###

- `eraRv2m`: r-vector to r-matrix

### Reference ###

- F.Mignard & M.Froeschle, Astron. Astrophys. 354, 732-739 (2000).

"""
function fk5hip()
    r5h = zeros((3, 3))
    s5h = zeros(3)
    ccall((:eraFk5hip, liberfa), Cvoid,
          (Ptr{Cdouble}, Ptr{Cdouble}),
          r5h, s5h)
    r5h, s5h
end

"""
    fk5hz(r5, d5, date1, date2)

Transform an FK5 (J2000.0) star position into the system of the
Hipparcos catalogue, assuming zero Hipparcos proper motion.

### Given ###

- `r5`: FK5 RA (radians), equinox J2000.0, at date
- `d5`: FK5 Dec (radians), equinox J2000.0, at date
- `date1`, `date2`: TDB date (Notes 1,2)

### Returned ###

- `rh`: Hipparcos RA (radians)
- `dh`: Hipparcos Dec (radians)

### Notes ###

1. This function converts a star position from the FK5 system to
   the Hipparcos system, in such a way that the Hipparcos proper
   motion is zero.  Because such a star has, in general, a non-zero
   proper motion in the FK5 system, the function requires the date
   at which the position in the FK5 system was determined.

2. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

3. The FK5 to Hipparcos transformation is modeled as a pure
   rotation and spin;  zonal errors in the FK5 catalogue are not
   taken into account.

4. The position returned by this function is in the Hipparcos
   reference system but at date date1+date2.

5. See also eraFk52h, eraH2fk5, eraHfk5z.

### Called ###

- `eraS2c`: spherical coordinates to unit vector
- `eraFk5hip`: FK5 to Hipparcos rotation and spin
- `eraSxp`: multiply p-vector by scalar
- `eraRv2m`: r-vector to r-matrix
- `eraTrxp`: product of transpose of r-matrix and p-vector
- `eraPxp`: vector product of two p-vectors
- `eraC2s`: p-vector to spherical
- `eraAnp`: normalize angle into range 0 to 2pi

### Reference ###

- F.Mignard & M.Froeschle, 2000, Astron.Astrophys. 354, 732-739.

"""
function fk5hz(r5, d5, date1, date2)
    rh = Ref(0.0)
    dh = Ref(0.0)
    ccall((:eraFk5hz, liberfa), Cvoid,
          (Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
          r5, d5, date1, date2, rh, dh)
    rh[], dh[]
end

"""
    fw2xy(gamb, phib, psi, eps)

CIP X,Y given Fukushima-Williams bias-precession-nutation angles.

### Given ###

- `gamb`: F-W angle gamma_bar (radians)
- `phib`: F-W angle phi_bar (radians)
- `psi`: F-W angle psi (radians)
- `eps`: F-W angle epsilon (radians)

### Returned ###

- `x`, `y`: CIP unit vector X,Y

### Notes ###

1. Naming the following points:

         e = J2000.0 ecliptic pole,
         p = GCRS pole
         E = ecliptic pole of date,
   and   P = CIP,

   the four Fukushima-Williams angles are as follows:

      gamb = gamma = epE
      phib = phi = pE
      psi = psi = pEP
      eps = epsilon = EP

2. The matrix representing the combined effects of frame bias,
   precession and nutation is:

      NxPxB = R_1(-epsA).R_3(-psi).R_1(phib).R_3(gamb)

   The returned values x,y are elements [2][0] and [2][1] of the
   matrix.  Near J2000.0, they are essentially angles in radians.

### Called ###

- `eraFw2m`: F-W angles to r-matrix
- `eraBpn2xy`: extract CIP X,Y coordinates from NPB matrix

### Reference ###

- Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351

"""
function fw2xy(gamb, phib, psi, eps)
    x = Ref(0.0)
    y = Ref(0.0)
    ccall((:eraFw2xy, liberfa), Cvoid,
          (Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
          gamb, phib, psi, eps, x, y)
    x[], y[]
end

"""
    fad03(t)

Fundamental argument, IERS Conventions (2003):

mean elongation of the Moon from the Sun.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned ###

- `D`, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003) and
   is from Simon et al. (1994).

### References ###

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

- Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

"""
fad03

"""
    fae03(t)

Fundamental argument, IERS Conventions (2003): Mean longitude of Earth.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned ###

- Mean longitude of Earth, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003) and
   comes from Souchay et al. (1999) after Simon et al. (1994).

### References ###

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

- Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

- Souchay, J., Loysel, B., Kinoshita, H., Folgueira, M. 1999,
    Astron.Astrophys.Supp.Ser. 135, 111

"""
fae03

"""
    faf03(t)

Fundamental argument, IERS Conventions (2003): Mean longitude of the Moon minus
mean longitude of the ascending node.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned ###

- `F`, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003) and
   is from Simon et al. (1994).

### References ###

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

- Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

"""
faf03

"""
    faju03(t)

Fundamental argument, IERS Conventions (2003): Mean longitude of Jupiter.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned ###

- Mean longitude of Jupiter, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003) and
   comes from Souchay et al. (1999) after Simon et al. (1994).

### References ###

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

- Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

- Souchay, J., Loysel, B., Kinoshita, H., Folgueira, M. 1999,
    Astron.Astrophys.Supp.Ser. 135, 111

"""
faju03

"""
    fal03(t)

Fundamental argument, IERS Conventions (2003):

mean anomaly of the Moon.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned ###

- `l`, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003) and
   is from Simon et al. (1994).

### References ###

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

- Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

"""
fal03

"""
    falp03(t)

Fundamental argument, IERS Conventions (2003):

mean anomaly of the Sun.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned ###

- `l'`, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003) and
   is from Simon et al. (1994).

### References ###

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

- Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

"""
falp03

"""
    fama03(t)

Fundamental argument, IERS Conventions (2003): Mean longitude of Mars.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned ###

- Mean longitude of Mars, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003) and
   comes from Souchay et al. (1999) after Simon et al. (1994).

### References ###

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

- Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

- Souchay, J., Loysel, B., Kinoshita, H., Folgueira, M. 1999,
    Astron.Astrophys.Supp.Ser. 135, 111

"""
fama03

"""
    fame03(t)

Fundamental argument, IERS Conventions (2003): Mean longitude of Mercury.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned ###

- Mean longitude of Mercury, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003) and
   comes from Souchay et al. (1999) after Simon et al. (1994).

### References ###

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

- Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

- Souchay, J., Loysel, B., Kinoshita, H., Folgueira, M. 1999,
    Astron.Astrophys.Supp.Ser. 135, 111

"""
fame03

"""
    fane03(t)

Fundamental argument, IERS Conventions (2003): Mean longitude of Neptune.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned ###

- Mean longitude of Neptune, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003) and
   is adapted from Simon et al. (1994).

### References ###

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

- Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

"""
fane03

"""
    faom03(t)

Fundamental argument, IERS Conventions (2003): Mean longitude of the Moon's ascending node.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned ###

- `Omega`, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003) and
   is from Simon et al. (1994).

### References ###

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

- Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

"""
faom03

"""
    fapa03(t)

Fundamental argument, IERS Conventions (2003):

general accumulated precession in longitude.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned ###

- General precession in longitude, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003).  It
   is taken from Kinoshita & Souchay (1990) and comes originally
   from Lieske et al. (1977).

### References ###

- Kinoshita, H. and Souchay J. 1990, Celest.Mech. and Dyn.Astron.
    48, 187

- Lieske, J.H., Lederle, T., Fricke, W. & Morando, B. 1977,
    Astron.Astrophys. 58, 1-16

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

"""
fapa03

"""
    fasa03(t)

Fundamental argument, IERS Conventions (2003): Mean longitude of Saturn.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned ###

- Mean longitude of Saturn, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003) and
   comes from Souchay et al. (1999) after Simon et al. (1994).

### References ###

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

- Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

- Souchay, J., Loysel, B., Kinoshita, H., Folgueira, M. 1999,
    Astron.Astrophys.Supp.Ser. 135, 111

"""
fasa03

"""
    faur03(t)

Fundamental argument, IERS Conventions (2003): Mean longitude of Uranus.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned  ###

- Mean longitude of Uranus, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003) and
   is adapted from Simon et al. (1994).

### References ###

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

- Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

"""
faur03

"""
    fave03(t)

Fundamental argument, IERS Conventions (2003): Mean longitude of Venus.

### Given ###

- `t`: TDB, Julian centuries since J2000.0 (Note 1)

### Returned ###

- Mean longitude of Venus, radians (Note 2)

### Notes ###

1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

2. The expression used is as adopted in IERS Conventions (2003) and
   comes from Souchay et al. (1999) after Simon et al. (1994).

### References ###

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

- Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

- Souchay, J., Loysel, B., Kinoshita, H., Folgueira, M. 1999,
    Astron.Astrophys.Supp.Ser. 135, 111

"""
fave03

for name in ("fad03",
             "fae03",
             "faf03",
             "faju03",
             "fal03",
             "falp03",
             "fama03",
             "fame03",
             "fane03",
             "faom03",
             "fapa03",
             "fasa03",
             "faur03",
             "fave03")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval ($f)(d) = ccall(($fc, liberfa), Cdouble, (Cdouble,), d)
end

"""
    fk52h(ra, dec, dra, ddec, px, rv)

Transform FK5 (J2000.0) star data into the Hipparcos system.

### Given (all FK5, equinox J2000.0, epoch J2000.0) ###

- `r5`: RA (radians)
- `d5`: Dec (radians)
- `dr5`: Proper motion in RA (dRA/dt, rad/Jyear)
- `dd5`: Proper motion in Dec (dDec/dt, rad/Jyear)
- `px5`: Parallax (arcsec)
- `rv5`: Radial velocity (km/s, positive = receding)

### Returned (all Hipparcos, epoch J2000.0) ###

- `rh`: RA (radians)
- `dh`: Dec (radians)
- `drh`: proper motion in RA (dRA/dt, rad/Jyear)
- `ddh`: proper motion in Dec (dDec/dt, rad/Jyear)
- `pxh`: parallax (arcsec)
- `rvh`: radial velocity (km/s, positive = receding)

### Notes ###

1. This function transforms FK5 star positions and proper motions
   into the system of the Hipparcos catalog.

2. The proper motions in RA are dRA/dt rather than
   cos(Dec)*dRA/dt, and are per year rather than per century.

3. The FK5 to Hipparcos transformation is modeled as a pure
   rotation and spin;  zonal errors in the FK5 catalog are not
   taken into account.

4. See also eraH2fk5, eraFk5hz, eraHfk5z.

### Called ###

- `eraStarpv`: star catalog data to space motion pv-vector
- `eraFk5hip`: FK5 to Hipparcos rotation and spin
- `eraRxp`: product of r-matrix and p-vector
- `eraPxp`: vector product of two p-vectors
- `eraPpp`: p-vector plus p-vector
- `eraPvstar`: space motion pv-vector to star catalog data

### Reference ###

- F.Mignard & M.Froeschle, Astron. Astrophys. 354, 732-739 (2000).

"""
function fk52h(ra, dec, dra, ddec, px, rv)
    r = Ref(0.0)
    d = Ref(0.0)
    dr = Ref(0.0)
    dd = Ref(0.0)
    p = Ref(0.0)
    v = Ref(0.0)
    ccall((:eraFk52h, liberfa), Cvoid,
            (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
            ra, dec, dra, ddec, px, rv, r, d, dr, dd, p, v)
    r[], d[], dr[], dd[], p[], v[]
end

"""
    fw2m(x, y, s, t)

Form rotation matrix given the Fukushima-Williams angles.

### Given ###

- `gamb`: F-W angle gamma_bar (radians)
- `phib`: F-W angle phi_bar (radians)
- `psi`: F-W angle psi (radians)
- `eps`: F-W angle epsilon (radians)

### Returned ###

- `r`: Rotation matrix

### Notes ###

1. Naming the following points:

         e = J2000.0 ecliptic pole,
         p = GCRS pole,
         E = ecliptic pole of date,
   and   P = CIP,

   the four Fukushima-Williams angles are as follows:

      gamb = gamma = epE
      phib = phi = pE
      psi = psi = pEP
      eps = epsilon = EP

2. The matrix representing the combined effects of frame bias,
   precession and nutation is:

      NxPxB = R_1(-eps).R_3(-psi).R_1(phib).R_3(gamb)

3. Three different matrices can be constructed, depending on the
   supplied angles:

   o  To obtain the nutation x precession x frame bias matrix,
      generate the four precession angles, generate the nutation
      components and add them to the psi_bar and epsilon_A angles,
      and call the present function.

   o  To obtain the precession x frame bias matrix, generate the
      four precession angles and call the present function.

   o  To obtain the frame bias matrix, generate the four precession
      angles for date J2000.0 and call the present function.

   The nutation-only and precession-only matrices can if necessary
   be obtained by combining these three appropriately.

### Called ###

- `eraIr`: initialize r-matrix to identity
- `eraRz`: rotate around Z-axis
- `eraRx`: rotate around X-axis

### Reference ###

- Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351

"""
function fw2m(x, y, s, t)
    r = zeros((3, 3))
    ccall((:eraFw2m, liberfa), Cvoid,
            (Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
            x, y, s, t, r)
    r
end


text:
function calc_B_0(cur_reactor::AbstractReactor)
  cur_constraint = cur_reactor.constraint

  cur_symbol = Symbol("B_0_from_$(cur_constraint)")

  cur_function = getfield(Fussy, cur_symbol)

  cur_B = cur_function(cur_reactor)

  cur_B
end


text:
# Blender
export freecad

#=

To dribble Blender, do:

freecad_exe_path = raw"C:\Program Files\Blender Foundation\Blender 2.91\freecad.exe"
using OutputCollectors
oc = OutputCollector(`$freecad_exe_path $["--python", joinpath(@__DIR__, "KhepriServer.py")]`, verbose=true)
=#

#=
=#

# FRC is a subtype of Python
parse_signature(::Val{:FRCAD}, sig::T) where {T} = parse_signature(Val(:PY), sig)
encode(::Val{:FRCAD}, t::Val{T}, c::IO, v) where {T} = encode(Val(:PY), t, c, v)
decode(::Val{:FRCAD}, t::Val{T}, c::IO) where {T} = decode(Val(:PY), t, c)
encode(ns::Val{:FRCAD}, t::Tuple{T1,T2,T3}, c::IO, v) where {T1,T2,T3} =
  begin
    encode(ns, T1(), c, v[1])
    encode(ns, T2(), c, v[2])
    encode(ns, T3(), c, v[3])
  end
decode(ns::Val{:FRCAD}, t::Tuple{T1,T2,T3}, c::IO) where {T1,T2,T3} =
  (decode(ns, T1(), c),
   decode(ns, T2(), c),
   decode(ns, T3(), c))

@encode_decode_as(:FRCAD, Val{:Id}, Val{:size})
@encode_decode_as(:FRCAD, Val{:MatId}, Val{:size})
@encode_decode_as(:FRCAD, Val{:Length}, Val{:float})

encode(::Val{:FRCAD}, t::Union{Val{:Point3d},Val{:Vector3d}}, c::IO, p) =
  encode(Val(:PY), Val(:float3), c, raw_point(p))
decode(::Val{:FRCAD}, t::Val{:Point3d}, c::IO) =
  xyz(decode(Val(:PY), Val(:float3), c)..., world_cs)
decode(::Val{:FRCAD}, t::Val{:Vector3d}, c::IO) =
  vxyz(decode(Val(:PY), Val(:float3), c)..., world_cs)

encode(ns::Val{:FRCAD}, t::Val{:Frame3d}, c::IO, v) = begin
  encode(ns, Val(:Point3d), c, v)
  t = v.cs.transform
  encode(Val(:PY), Val(:float3), c, (t[1,1], t[2,1], t[3,1]))
  encode(Val(:PY), Val(:float3), c, (t[1,2], t[2,2], t[3,2]))
  encode(Val(:PY), Val(:float3), c, (t[1,3], t[2,3], t[3,3]))
end

decode(ns::Val{:FRCAD}, t::Val{:Frame3d}, c::IO) =
  u0(cs_from_o_vx_vy_vz(
      decode(ns, Val(:Point3d), c),
      decode(ns, Val(:Vector3d), c),
      decode(ns, Val(:Vector3d), c),
      decode(ns, Val(:Vector3d), c)))

freecad_api = @remote_functions :FRCAD """
def find_or_create_collection(name:str, active:bool, color:RGBA)->str:
def get_current_collection()->str:
def set_current_collection(name:str)->None:
def delete_all_shapes_in_collection(name:str)->None:
def delete_all_shapes()->None:
def delete_shape(name:Id)->None:
def select_shape(name:Id)->None:
def deselect_shape(name:Id)->None:
def deselect_all_shapes()->None:
def get_material(name:str)->MatId:
def get_freecadkit_material(ref:str)->MatId:
def new_material(name:str, diffuse_color:RGBA, metallic:float, specular:float, roughness:float, clearcoat:float, clearcoat_roughness:float, ior:float, transmission:float, transmission_roughness:float, emission:RGBA, emission_strength:float)->MatId:
def new_metal_material(name:str, color:RGBA, roughness:float, ior:float)->MatId:
def new_glass_material(name:str, color:RGBA, roughness:float, ior:float)->MatId:
def new_mirror_material(name:str, color:RGBA)->MatId:
def line(ps:List[Point3d], closed:bool, mat:MatId)->Id:
def spline(ps:List[Point3d], closed:bool, mat:MatId)->Id:
def nurbs(order:int, ps:List[Point3d], closed:bool, mat:MatId)->Id:
def objmesh(verts:List[Point3d], edges:List[Tuple[int,int]], faces:List[List[int]], smooth:bool, mat:MatId)->Id:
def trig(p1:Point3d, p2:Point3d, p3:Point3d, mat:MatId)->Id:
def quad(p1:Point3d, p2:Point3d, p3:Point3d, p4:Point3d, mat:MatId)->Id:
def quad_strip(ps:List[Point3d], qs:List[Point3d], smooth:bool, mat:MatId)->Id:
def quad_strip_closed(ps:List[Point3d], qs:List[Point3d], smooth:bool, mat:MatId)->Id:
def ngon(ps:List[Point3d], pivot:Point3d, smooth:bool, mat:MatId)->Id:
def polygon(ps:List[Point3d], mat:MatId)->Id:
def polygon_with_holes(pss:List[List[Point3d]], mat:MatId)->Id:
def quad_surface(ps:List[Point3d], nu:int, nv:int, closed_u:bool, closed_v:bool, smooth:bool, mat:MatId)->Id:
def circle(c:Point3d, v:Vector3d, r:Length, mat:MatId)->Id:
def cuboid(verts:List[Point3d], mat:MatId)->Id:
def pyramid_frustum(bs:List[Point3d], ts:List[Point3d], smooth:bool, bmat:MatId, tmat:MatId, smat:MatId)->Id:
def sphere(center:Point3d, radius:Length, mat:MatId)->Id:
def cone_frustum(b:Point3d, br:Length, t:Point3d, tr:Length, bmat:MatId, tmat:MatId, smat:MatId)->Id:
def box(p:Point3d, vx:Vector3d, vy:Vector3d, dx:Length, dy:Length, dz:Length, mat:MatId)->Id:
def text(txt:str, p:Point3d, vx:Vector3d, vy:Vector3d, size:float)->Id:
def area_light(p:Point3d, v:Vector3d, size:float, color:RGBA, strength:float)->Id:
def sun_light(p:Point3d, v:Vector3d)->Id:
def light(p:Point3d, type:str)->Id:
def camera_from_view()->None:
def set_view(camera:Point3d, target:Point3d, lens:float)->None:
def get_view()->Tuple[Point3d, Point3d, float]:
def set_render_size(width:int, height:int)->None:
def set_render_path(filepath:str)->None:
def default_renderer()->None:
def cycles_renderer(samples:int, denoising:bool, motion_blur:bool, transparent:bool, exposure:float)->None:
def freestylesvg_renderer(thickness:float, crease_angle:float)->None:
def clay_renderer(samples:int, denoising:bool, motion_blur:bool, transparent:bool)->None:
def set_sun(latitude:float, longitude:float, elevation:float, year:int, month:int, day:int, time:float, UTC_zone:float, use_daylight_savings:bool)->None:
def set_sky(turbidity:float)->None:
def set_max_repeated(n:Int)->Int:
def freecad_cmd(expr:str)->None:
def wall(ps:List[Point3d], height:Length, align:str, mat:MatId)->Id:
"""

abstract type FRCADKey end
const FRCADId = Union{Int32,String} # Although shapes and materials are ints, layers are strings
const FRCADIds = Vector{FRCADId}
const FRCADRef = NativeRef{FRCADKey, FRCADId}
const FRCADRefs = Vector{FRCADRef}
const FRCADEmptyRef = EmptyRef{FRCADKey, FRCADId}
const FRCADUniversalRef = UniversalRef{FRCADKey, FRCADId}
const FRCADUnionRef = UnionRef{FRCADKey, FRCADId}
const FRCADSubtractionRef = SubtractionRef{FRCADKey, FRCADId}
const FRCAD = SocketBackend{FRCADKey, FRCADId}

const KhepriServerPath = Parameter(abspath(@__DIR__, "FreeCADServer.py"))
export headless_freecad
const headless_freecad = Parameter(false)
const starting_freecad = Parameter(false)

start_freecad() =
  starting_freecad() ?
    sleep(1) : # Just wait a little longer
    let freecad_cmd = Sys.iswindows() ?
		  let versions = filter(p->!isnothing(match(r"FreeCAD", p)), readdir("C:/Program Files/", join=true))
			isempty(versions) ?
			  error("Could not find FreeCAD!") :
			  joinpath(versions[1], "bin/FreeCAD.exe")
		  end :
    	  "freecad"
	  starting_freecad(true)
      run(detach(headless_freecad() ?
            `$(freecad_cmd) -P $(@__DIR__) --console $(KhepriServerPath())` :
      	    `$(freecad_cmd) -P $(@__DIR__) $(KhepriServerPath())`),
    	  wait=false)
    end

#
#=
To support reload in FreeCAD, redefine:

start_freecad() =
  starting_freecad() ?
    sleep(1) : # Just wait a little longer
    let start = abspath(@__DIR__, "StartKhepri.py"),
		    freecad_cmd = Sys.iswindows() ?
		      let versions = filter(p->!isnothing(match(r"FreeCAD", p)), readdir("C:/Program Files/", join=true))
			      isempty(versions) ?
			        error("Could not find FreeCAD!") :
			        joinpath(versions[1], "bin/FreeCAD.exe")
		      end :
    	    "freecad"
	  starting_freecad(true)
      run(detach(headless_freecad() ?
            `$(freecad_cmd) -P $(@__DIR__) --console $(start)` :
      	    `$(freecad_cmd) -P $(@__DIR__) $(start)`),
    	  wait=false)
    end
=#

KhepriBase.retry_connecting(b::FRCAD) =
  (@info("Starting $(b.name)."); start_freecad(); sleep(2))

KhepriBase.after_connecting(b::FRCAD) =
  begin
	starting_freecad(false)
	# #set_material(freecad, material_basic, )
	# set_material(freecad, material_metal, "asset_base_id:f1774cb0-b679-46b4-879e-e7223e2b4b5f asset_type:material")
	# #set_material(freecad, material_glass, "asset_base_id:ee2c0812-17f5-40d4-992c-68c5a66261d7 asset_type:material")
	# set_material(freecad, material_glass, "asset_base_id:ffa3c281-6184-49d8-b05e-8c6e9fe93e68 asset_type:material")
	# set_material(freecad, material_wood, "asset_base_id:d5097824-d5a1-4b45-ab5b-7b16bdc5a627 asset_type:material")
	# #set_material(freecad, material_concrete, "asset_base_id:0662b3bf-a762-435d-9407-e723afd5eafc asset_type:material")
	# set_material(freecad, material_concrete, "asset_base_id:df1161da-050c-4638-b376-38ced992ec18 asset_type:material")
	# set_material(freecad, material_plaster, "asset_base_id:c674137d-cfae-45f1-824f-e85dc214a3af asset_type:material")
	#
	# #set_material(freecad, material_grass, "asset_base_id:97b171b4-2085-4c25-8793-2bfe65650266 asset_type:material")
	# #set_material(freecad, material_grass, "asset_base_id:7b05be22-6bed-4584-a063-d0e616ddea6a asset_type:material")
	# set_material(freecad, material_grass, "asset_base_id:b4be2338-d838-433b-9f0d-2aa9b97a0a8a asset_type:material")
	# set_material(freecad, material_clay, b -> b_plastic_material(b, "Clay", rgb(0.9, 0.9, 0.9),	1.0))
	# We will use the same view as Rhino
	b_set_view(b, xyz(43.11,-74.67,49.78), xyz(-0.19,0.33,-0.22), 50, 22)
  end

const freecad = FRCAD("FreeCAD", freecad_port, freecad_api)

KhepriBase.has_boolean_ops(::Type{FRCAD}) = HasBooleanOps{true}()

KhepriBase.backend(::FRCADRef) = freecad
KhepriBase.void_ref(b::FRCAD) = FRCADRef(-1 % Int32)

# Primitives

KhepriBase.b_line(b::FRCAD, ps, mat) =
  @remote(b, line(ps, false, mat))

KhepriBase.b_polygon(b::FRCAD, ps, mat) =
	@remote(b, line(ps, true, mat))

KhepriBase.b_nurbs_curve(b::FRCAD, order, ps, knots, weights, closed, mat) =
  @remote(b, nurbs(order, ps, closed, mat))

KhepriBase.b_trig(b::FRCAD, p1, p2, p3, mat) =
  @remote(b, trig(p1, p2, p3, mat))

KhepriBase.b_quad(b::FRCAD, p1, p2, p3, p4, mat) =
	@remote(b, quad(p1, p2, p3, p4, mat))

KhepriBase.b_ngon(b::FRCAD, ps, pivot, smooth, mat) =
	@remote(b, ngon(ps, pivot, smooth, mat))

KhepriBase.b_quad_strip(b::FRCAD, ps, qs, smooth, mat) =
  @remote(b, quad_strip(ps, qs, smooth, mat))

KhepriBase.b_quad_strip_closed(b::FRCAD, ps, qs, smooth, mat) =
  @remote(b, quad_strip_closed(ps, qs, smooth, mat))

KhepriBase.b_surface_polygon(b::FRCAD, ps, mat) =
  @remote(b, polygon(ps, mat))

KhepriBase.b_surface_polygon_with_holes(b::FRCAD, ps, qss, mat) =
  @remote(b, polygon_with_holes([ps, qss...], mat))

KhepriBase.b_surface_circle(b::FRCAD, c, r, mat) =
  @remote(b, circle(c, vz(1, c.cs), r, mat))

KhepriBase.b_surface_grid(b::FRCAD, ptss, closed_u, closed_v, smooth_u, smooth_v, mat) =
  let (nu, nv) = size(ptss)
	smooth_u && smooth_v ?
	  @remote(b, quad_surface(vcat(ptss...), nu, nv, closed_u, closed_v, true, mat)) :
	  smooth_u ?
	  	(closed_u ?
          vcat([b_quad_strip_closed(b, ptss[:,i], ptss[:,i+1], true, mat) for i in 1:nv-1],
	           closed_v ? [b_quad_strip_closed(b, ptss[:,end], ptss[:,1], true, mat)] : []) :
	      vcat([b_quad_strip(b, ptss[:,i], ptss[:,i+1], true, mat) for i in 1:nv-1],
	           closed_v ? [b_quad_strip(b, ptss[:,end], ptss[:,1], true, mat)] : [])) :
 	    (closed_v ?
           vcat([b_quad_strip_closed(b, ptss[i,:], ptss[i+1,:], smooth_v, mat) for i in 1:nu-1],
  	         	closed_u ? [b_quad_strip_closed(b, ptss[end,:], ptss[1,:], smooth_v, mat)] : []) :
  	       vcat([b_quad_strip(b, ptss[i,:], ptss[i+1,:], smooth_v, mat) for i in 1:nu-1],
  	          	closed_u ? [b_quad_strip(b, ptss[end,:], ptss[1,:], smooth_v, mat)] : []))
  end

KhepriBase.b_generic_pyramid_frustum(b::FRCAD, bs, ts, smooth, bmat, tmat, smat) =
  @remote(b, pyramid_frustum(bs, ts, smooth, bmat, tmat, smat))

KhepriBase.b_cone(b::FRCAD, cb, r, h, bmat, smat) =
  @remote(b, cone_frustum(cb, r, add_z(cb, h), 0, bmat, bmat, smat))

KhepriBase.b_cone_frustum(b::FRCAD, cb, rb, h, rt, bmat, tmat, smat) =
  @remote(b, cone_frustum(cb, rb, add_z(cb, h), rt, bmat, tmat, smat))

KhepriBase.b_cylinder(b::FRCAD, cb, r, h, bmat, tmat, smat) =
  @remote(b, cone_frustum(cb, r, add_z(cb, h), r, bmat, tmat, smat))

KhepriBase.b_box(b::FRCAD, c, dx, dy, dz, mat) =
  @remote(b, box(c, vx(1, c.cs), vy(1, c.cs), dx, dy, dz, mat))

KhepriBase.b_cuboid(b::FRCAD, pb0, pb1, pb2, pb3, pt0, pt1, pt2, pt3, mat) =
  @remote(b, cuboid([pb0, pb1, pb2, pb3, pt0, pt1, pt2, pt3], mat))

KhepriBase.b_sphere(b::FRCAD, c, r, mat) =
  @remote(b, sphere(c, r, mat))

# BIM

KhepriBase.b_wall(b::FRCAD, w_path, w_height, l_thickness, r_thickness, family) =
  path_length(w_path) < path_tolerance() ?
  	void_ref(b) :
    @remote(b, wall(path_vertices(w_path), w_height, "Center", -1))

# Materials

KhepriBase.b_get_material(b::FRCAD, spec::AbstractString) =
  startswith(spec, "asset_base_id") ?
    @remote(b, get_freecadkit_material(spec)) :
    @remote(b, get_material(spec))

KhepriBase.b_get_material(b::FRCAD, spec::Function) = spec(b)

#=
Important source of materials:

1. Activate Blender's freecadkit addon:
https://www.freecadkit.com/get-freecadkit/

2. Browse BlenderKit's material database:
https://www.freecadkit.com/asset-gallery?query=category_subtree:material

3. Select material and copy reference, e.g.:
asset_base_id:ced25dc0-d461-42f7-aa03-85cb88f671a1 asset_type:material

4. Install material and retrive its id with:
b_get_material(freecad, "asset_base_id:ced25dc0-d461-42f7-aa03-85cb88f671a1 asset_type:material")
=#

KhepriBase.b_new_material(b::FRCAD, name,
						  base_color,
						  metallic, specular, roughness,
	                 	  clearcoat, clearcoat_roughness,
						  ior,
						  transmission, transmission_roughness,
	                 	  emission_color,
						  emission_strength) =
  @remote(b, new_material(name,
  						  convert(RGBA, base_color),
						  metallic, specular, roughness,
  						  clearcoat, clearcoat_roughness,
  				  		  ior,
  				  		  transmission, transmission_roughness,
						  convert(RGBA, emission_color), emission_strength))

KhepriBase.b_plastic_material(b::FRCAD, name, color, roughness) =
  @remote(b, new_material(name, convert(RGBA, color), 0.0, 1.0, roughness, 0.0, 0.0, 1.4, 0.0, 0.0, RGBA(0.0, 0.0, 0.0, 1.0), 0.0))

KhepriBase.b_metal_material(b::FRCAD, name, color, roughness, ior) =
  @remote(b, new_metal_material(name, convert(RGBA, color), roughness, ior))

KhepriBase.b_glass_material(b::FRCAD, name, color, roughness, ior) =
  @remote(b, new_glass_material(name, convert(RGBA, color), roughness, ior))

KhepriBase.b_mirror_material(b::FRCAD, name, color) =
  @remote(b, new_mirror_material(name, convert(RGBA, color)))

#KhepriBase.b_translucent_material(b::FRCAD, name, diffuse, specular, roughness, reflect, transmit, bump_map)
#KhepriBase.b_substrate_material(name, diffuse, specular, roughness, bump_map)

#=

Default families

=#
export freecad_family_materials
freecad_family_materials(m1, m2=m1, m3=m2, m4=m3) = (materials=(m1, m2, m3, m4), )

KhepriBase.b_layer(b::FRCAD, name, active, color) =
  @remote(b, find_or_create_collection(name, active, color))
KhepriBase.b_current_layer(b::FRCAD) =
  @remote(b, get_current_collection())
KhepriBase.b_current_layer(b::FRCAD, layer) =
  @remote(b, set_current_collection(layer))
KhepriBase.b_all_shapes_in_layer(b::FRCAD, layer) =
  @remote(b, all_shapes_in_collection(layer))
KhepriBase.b_delete_all_shapes_in_layer(b::FRCAD, layer) =
  @remote(b, delete_all_shapes_in_collection(layer))

realize(b::FRCAD, s::EmptyShape) =
  FRCADEmptyRef()
realize(b::FRCAD, s::UniversalShape) =
  FRCADUniversalRef()

KhepriBase.b_set_view(b::FRCAD, camera::Loc, target::Loc, lens::Real, aperture::Real) =
  @remote(b, set_view(camera, target, lens))

KhepriBase.b_get_view(b::FRCAD) =
  @remote(b, get_view())

KhepriBase.b_zoom_extents(b::FRCAD) = @remote(b, ZoomExtents())

KhepriBase.b_set_view_top(b::FRCAD) = @remote(b, ViewTop())

KhepriBase.b_set_time_place(b::FRCAD, date, latitude, longitude, elevation, meridian) =
  @remote(b, set_sun(latitude, longitude, elevation, year(date), month(date), day(date), hour(date)+minute(date)/60, meridian, false))

KhepriBase.b_set_sky(b::FRCAD, turbidity, sun) =
  @remote(b, set_sky(turbidity)) #Add withsun

KhepriBase.b_delete_ref(b::FRCAD, r::FRCADId) =
  @remote(b, delete_shape(r))

KhepriBase.b_delete_all_refs(b::FRCAD) =
  @remote(b, delete_all_shapes())

#=
backend_set_length_unit(b::FRCAD, unit::String) = @remote(b, SetLengthUnit(unit))

# Dimensions

const FRCADDimensionStyles = Dict(:architectural => "_ARCHTICK", :mechanical => "")

backend_dimension(b::FRCAD, p0::Loc, p1::Loc, p::Loc, scale::Real, style::Symbol) =
    @remote(b, CreateAlignedDimension(p0, p1, p,
        scale,
        FRCADDimensionStyles[style]))

backend_dimension(b::FRCAD, p0::Loc, p1::Loc, sep::Real, scale::Real, style::Symbol) =
    let v = p1 - p0
        angle = pol_phi(v)
        dimension(p0, p1, add_pol(p0, sep, angle + pi/2), scale, style, b)
    end

# Blocks

realize(b::FRCAD, s::Block) =
    @remote(b, CreateBlockFromShapes(s.name, collect_ref(s.shapes)))

realize(b::FRCAD, s::BlockInstance) =
    @remote(b, CreateBlockInstance(
        collect_ref(s.block)[1],
        center_scaled_cs(s.loc, s.scale, s.scale, s.scale)))
=#
#=

# Manual process
@time for i in 1:1000 for r in 1:10 circle(x(i*10), r) end end

# Create block...
Khepri.create_block("Foo", [circle(radius=r) for r in 1:10])

# ...and instantiate it
@time for i in 1:1000 Khepri.instantiate_block("Foo", x(i*10), 0) end

=#
#=
# Lights
KhepriBase.b_pointlight(b::FRCAD, loc::Loc, color::RGB, range::Real, intensity::Real) =
  # HACK: Fix this
  @remote(b, SpotLight(loc, intensity, range, loc+vz(-1)))

KhepriBase.b_spotlight(b::FRCAD, loc::Loc, dir::Vec, hotspot::Real, falloff::Real) =
    @remote(b, SpotLight(loc, hotspot, falloff, loc + dir))

KhepriBase.b_ieslight(b::FRCAD, file::String, loc::Loc, dir::Vec, alpha::Real, beta::Real, gamma::Real) =
    @remote(b, IESLight(file, loc, loc + dir, vxyz(alpha, beta, gamma)))

# User Selection

backend_shape_from_ref(b::FRCAD, r) =
  let c = connection(b),
      code = @remote(b, ShapeCode(r)),
      ref = DynRefs(b=>FRCADRef(r))
    if code == 1 # Point
        point(@remote(b, PointPosition(r)),
              ref=ref)
    elseif code == 2
        circle(loc_from_o_vz(@remote(b, CircleCenter(r)), @remote(b, CircleNormal(r))),
               @remote(b, CircleRadius(r)),
               ref=ref)
    elseif 3 <= code <= 6
        line(@remote(b, LineVertices(r)),
             ref=ref)
    elseif code == 7
        let tans = @remote(b, SplineTangents(r))
            if length(tans[1]) < 1e-20 && length(tans[2]) < 1e-20
                closed_spline(@remote(b, SplineInterpPoints(r))[1:end-1],
                              ref=ref)
            else
                spline(@remote(b, SplineInterpPoints(r)), tans[1], tans[2],
                       ref=ref)
            end
        end
    elseif code == 9
        let start_angle = mod(@remote(b, ArcStartAngle(r)), 2pi),
            end_angle = mod(@remote(b, ArcEndAngle(r)), 2pi)
            arc(loc_from_o_vz(@remote(b, ArcCenter(r)), @remote(b, ArcNormal(r))),
                @remote(b, ArcRadius(r)), start_angle, mod(end_angle - start_angle, 2pi),
                ref=ref)
        #=    if end_angle > start_angle
                arc(maybe_loc_from_o_vz(@remote(b, ArcCenter(r)), @remote(b, ArcNormal(r))),
                    @remote(b, ArcRadius(r)), start_angle, end_angle - start_angle,
                    ref=ref)
            else
                arc(maybe_loc_from_o_vz(@remote(b, ArcCenter(r)), @remote(b, ArcNormal(r))),
                    @remote(b, ArcRadius(r)), end_angle, start_angle - end_angle,
                    ref=ref)
            end=#
        end
    elseif code == 10
        let str = @remote(b, TextString(r)),
            height = @remote(b, TextHeight(r)),
            loc = @remote(b, TextPosition(r))
            text(str, loc, height, ref=ref)
        end
    elseif code == 11
        let str = @remote(b, MTextString(r)),
            height = @remote(b, MTextHeight(r)),
            loc = @remote(b, MTextPosition(r))
            text(str, loc, height, ref=ref)
        end
    elseif code == 16
        let pts = @remote(b, MeshVertices(r)),
            (type, n, m, n_closed, m_closed) = @remote(b, PolygonMeshData(r))
            surface_grid(reshape(pts, (n, m)), n_closed == 1, m_closed == 1, ref=ref)
        end
    elseif 12 <= code <= 14
        surface(Shapes1D[], ref=ref)
    elseif 103 <= code <= 106
        polygon(@remote(b, LineVertices(r)),
                ref=ref)
    elseif code == 107
        closed_spline(@remote(b, SplineInterpPoints(r))[1:end-1],
                      ref=ref)
    else
        #unknown(ref=ref)
        unknown(r, ref=ref) # To force copy
        #error("Unknown shape with code $(code)")
    end
  end
#
=#
#=
In case we need to realize an Unknown shape, we just copy it
=#
#=
realize(b::FRCAD, s::Unknown) =
    @remote(b, Copy(s.baseref))

backend_select_position(b::FRCAD, prompt::String) =
  begin
    @info "$(prompt) on the $(b) backend."
    let ans = @remote(b, GetPosition(prompt))
      length(ans) > 0 ? ans[1] : nothing
    end
  end

backend_select_positions(b::FRCAD, prompt::String) =
  let sel() =
    let p = select_position(prompt, b)
      if p == nothing
        []
      else
        [p, sel()...]
      end
    end
    sel()
  end

# HACK: The next operations should receive a set of shapes to avoid re-creating already existing shapes

backend_select_point(b::FRCAD, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetPoint)
backend_select_points(b::FRCAD, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetPoints)

backend_select_curve(b::FRCAD, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetCurve)
backend_select_curves(b::FRCAD, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetCurves)

backend_select_surface(b::FRCAD, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetSurface)
backend_select_surfaces(b::FRCAD, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetSurfaces)

backend_select_solid(b::FRCAD, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetSolid)
backend_select_solids(b::FRCAD, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetSolids)

backend_select_shape(b::FRCAD, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetShape)
backend_select_shapes(b::FRCAD, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetShapes)

backend_captured_shape(b::FRCAD, handle) =
  backend_shape_from_ref(b, @remote(b, GetShapeFromHandle(handle)))
backend_captured_shapes(b::FRCAD, handles) =
  map(handles) do handle
      backend_shape_from_ref(b, @remote(b, GetShapeFromHandle(handle)))
  end

backend_generate_captured_shape(b::FRCAD, s::Shape) =
    println("captured_shape(freecad, $(@remote(b, GetHandleFromShape(ref(s).value))))")
backend_generate_captured_shapes(b::FRCAD, ss::Shapes) =
  begin
    print("captured_shapes(freecad, [")
    for s in ss
      print(@remote(b, GetHandleFromShape(ref(s).value)))
      print(", ")
    end
    println("])")
  end

# Register for notification

backend_register_shape_for_changes(b::FRCAD, s::Shape) =
    let conn = connection(b)
        @remote(b, RegisterForChanges(ref(s).value))
        @remote(b, DetectCancel())
        s
    end

backend_unregister_shape_for_changes(b::FRCAD, s::Shape) =
    let conn = connection(b)
        @remote(b, UnregisterForChanges(ref(s).value))
        @remote(b, UndetectCancel())
        s
    end

backend_waiting_for_changes(b::FRCAD, s::Shape) =
    ! @remote(b, WasCanceled())

backend_changed_shape(b::FRCAD, ss::Shapes) =
    let conn = connection(b)
        changed = []
        while length(changed) == 0 && ! @remote(b, WasCanceled())
            changed =  @remote(b, ChangedShape())
            sleep(0.1)
        end
        if length(changed) > 0
            backend_shape_from_ref(b, changed[1])
        else
            nothing
        end
    end

# HACK: This should be filtered on the plugin, not here.
backend_all_shapes(b::FRCAD) =
  Shape[backend_shape_from_ref(b, r)
        for r in filter(r -> @remote(b, ShapeCode(r)) != 0, @remote(b, GetAllShapes()))]

backend_all_shapes_in_layer(b::FRCAD, layer) =
  Shape[backend_shape_from_ref(b, r) for r in @remote(b, GetAllShapesInLayer(layer))]

=#

KhepriBase.b_highlight_ref(b::FRCAD, r::FRCADId) =
  @remote(b, select_shape(r))

KhepriBase.b_unhighlight_ref(b::FRCAD, r::FRCADId) =
  @remote(b, deselect_shape(r))

KhepriBase.b_unhighlight_all_refs(b::FRCAD) =
  @remote(b, deselect_all_shapes())
#=
backend_pre_selected_shapes_from_set(ss::Shapes) =
  length(ss) == 0 ? [] : pre_selected_shapes_from_set(ss, backend(ss[1]))

# HACK: This must be implemented for all backends
backend_pre_selected_shapes_from_set(ss::Shapes, b::Backend) = []

backend_pre_selected_shapes_from_set(b::FRCAD, ss::Shapes) =
  let refs = map(id -> @remote(b, GetHandleFromShape(id)), @remote(b, GetPreSelectedShapes()))
    filter(s -> @remote(b, GetHandleFromShape(ref(s).value)) in refs, ss)
  end
backend_disable_update(b::FRCAD) =
  @remote(b, DisableUpdate())

backend_enable_update(b::FRCAD) =
  @remote(b, EnableUpdate())

# Render

#render exposure: [-3, +3] -> [-6, 21]
convert_render_exposure(b::FRCAD, v::Real) = -4.05*v + 8.8
#render quality: [-1, +1] -> [+1, +50]
convert_render_quality(b::FRCAD, v::Real) = round(Int, 25.5 + 24.5*v)

=#

KhepriBase.b_render_view(b::FRCAD, path::String) =
  let (camera, target, lens) = @remote(b, get_view())
    @remote(b, set_camera_view(camera, target, lens))
    @remote(b, set_render_size(render_width(), render_height()))
    @remote(b, set_render_path(path))
  	@remote(b, cycles_renderer(1200, true, false, false, render_exposure()))
  end

export render_svg
render_svg(b::FRCAD, path) =
  let (camera, target, lens) = @remote(b, get_view())
    @remote(b, set_camera_view(camera, target, lens))
    @remote(b, set_render_size(render_width(), render_height()))
    @remote(b, set_render_path(path))
    @remote(b, freestylesvg_renderer(1.0, deg2rad(135), 0.001, 0.0))
  end

#=
with_clay_model(f, level::Real=0) =
  begin
  	with(f, default_material, material_clay)
	area_light
=#


text:
@everywhere function putmsg(pid)
    mypid = myid()
    msg = "Hi from $mypid"
    rr = RemoteChannel(pid)
    put!(rr, msg)
    println(myid(), " sent ", msg, " to ", pid)
    return rr
end

@everywhere function getmsg(rr)
    msg  = fetch(rr)
    println(myid(), " got: ", msg)
end

rr = remotecall_fetch(putmsg, 2, 3)
remotecall_wait(getmsg, 3, rr)


text:
@testset "queue.jl" begin

    @testset "fifo" begin
        queue = Queue{Int64}()
        @test enqueue!(queue, 2) == [2]
        @test enqueue!(queue, 3) == [2, 3]
        @test add!(queue, 4) == [2, 3, 4]
        @test first(queue) == 2
        @test last(queue) == 4
        @test length(queue) == 3
        @test dequeue!(queue) == 2
        @test dequeue!(queue) == 3
        @test dequeue!(queue) == 4
        @test isempty(queue) == true
    end

    @testset "lifo" begin
        lifo = LifoQueue{Int64}()
        @test enqueue!(lifo, 2) == [2]
        @test enqueue!(lifo, 3) == [2, 3]
        @test first(lifo) == 2
        @test last(lifo) == 3
        @test length(lifo) == 2
        @test dequeue!(lifo) == 3
        @test dequeue!(lifo) == 2
        @test isempty(lifo) == true
    end

end


text:
"""
Bridges `CP.VectorDomain` to MILP by adding one binary variable per
possible combination.
"""
struct VectorDomain2MILPBridge{T} <: MOIBC.AbstractBridge
    vars::Vector{MOI.VariableIndex}
    vars_bin::Vector{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    con_choose_one::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}
    cons_values::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
end

function MOIBC.bridge_constraint(
    ::Type{VectorDomain2MILPBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.VectorDomain{T},
) where {T}
    return MOIBC.bridge_constraint(
        VectorDomain2MILPBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{VectorDomain2MILPBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.VectorDomain{T},
) where {T}
    vars, vars_bin = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:length(s.values)])

    con_choose_one = MOI.add_constraint(
        model,
        sum(one(T) .* vars),
        MOI.EqualTo(one(T))
    )

    f_scalars = MOIU.scalarize(f)
    values = collect(s.values)

    cons_values = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[
        MOI.add_constraint(
            model,
            f_scalars[i] - sum(one(T) * vars[j] * values[j][i] for j in 1:length(s.values)),
            MOI.EqualTo(zero(T))
        )
        for i in 1:s.dimension
    ]

    return VectorDomain2MILPBridge(vars, vars_bin, con_choose_one, cons_values)
end

function MOI.supports_constraint(
    ::Type{VectorDomain2MILPBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.VectorDomain{T}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{VectorDomain2MILPBridge{T}}) where {T}
    return [(MOI.ZeroOne,)]
end

function MOIB.added_constraint_types(::Type{VectorDomain2MILPBridge{T}}) where {T}
    return [
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
    ]
end

function MOI.get(b::VectorDomain2MILPBridge, ::MOI.NumberOfVariables)
    return length(b.vars)
end

function MOI.get(
    b::VectorDomain2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return length(b.vars_bin)
end

function MOI.get(
    b::VectorDomain2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return 1 + length(b.cons_values)
end

function MOI.get(
    b::VectorDomain2MILPBridge{T},
    ::MOI.ListOfVariableIndices,
) where {T}
    return copy(b.vars)
end

function MOI.get(
    b::VectorDomain2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return copy(b.vars_bin)
end

function MOI.get(
    b::VectorDomain2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return [b.con_choose_one, b.cons_values...]
end


text:
n = 10000000
randnstuple = [(x, y) for (x, y) in zip(rand(n), rand(n))]

num_inside = length(filter(x -> x[1]^2 + x[2]^2 <= 1, randnstuple))
println(num_inside / n * 4)


text:
include("delta.jl")

function epsilon(ii::Int64, jj::Int64, ksoc::Float64, Omega::Float64)

    energy0 = 0 # Int64
    njj = ceil(Int64,jj/2) - 1
    nii = ceil(Int64,ii/2) - 1

    if ii == jj

       energy0 = njj + 1/2 # Float64

    elseif iseven(ii+jj)

       energy0 = 1im*ksoc/sqrt(2)*(-1)^(ii)*(sqrt(njj+1)*delta(ii-jj-2) - sqrt(njj)*delta(ii-jj+2)) # Float64
       # with ii odd, (-1)^ii=-1 and ceil(jj/2)=(jj+1)/2 since jj is odd
       # with ii even, (-1)^ii=1 and ceil(jj/2)=jj/2 since jj is even

    elseif abs(ii-jj) == 1 && nii == njj # && isodd(ii+jj)

       energy0 = Omega/2.0 # Float64

    end

    return energy0

end

max_stars_repo_name: leios/SOCED

max_issues_repo_name: leios/SOCED

max_forks_repo_name: leios/SOCED


text:
using ArgParse # To have access to the command-line arguments
using SphericalHarmonics # To compute Legendre associate functions
using StaticArrays # To have access to fast arrays
using LinearAlgebra # To compute eigenvalues
using Plots # To have acces to plotting tools

max_stars_repo_name: KerwannTEP/LiRA

max_issues_repo_name: KerwannTEP/LiRA

max_forks_repo_name: KerwannTEP/LiRA


text:
mutable struct BasisSet # list of CGBFs
  bfs::Array{CGBF,1}
end

basisset() = BasisSet(CGBF[])

function push!(basis::BasisSet,cbf::CGBF)
  Base.push!(basis.bfs,cbf)
end

function build_basis(mol::Molecule,name="sto3g")
    data = basis_set_data[name]
    basis_set = basisset()
    for atom in mol.atomlist
        for btuple in data[atom.atno]
            sym,primlist = btuple
            for (I,J,K) in sym2power[sym]
                cbf = init_CGBF( (atom.x,atom.y,atom.z), (I,J,K) )
                push!(basis_set,cbf)
                # add primitives
                for (expn,coef) in primlist
                  push!(cbf,expn,coef)
                end
            end
        end
    end
    return basis_set
end

const sym2power = Dict{Any,Any}(
  'S' => [(0,0,0)],
  'P' => [(1,0,0),(0,1,0),(0,0,1)],
  'D' => [(2,0,0),(0,2,0),(0,0,2),(1,1,0),(1,0,1),(0,1,1)]
  )