text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to the EOP data.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/transformations/eop.jl
# ==================================

@testset "EOP data (IAU1980)" begin
    # Compare the fetched data to the one available locally
    # ==========================================================================

    # Fetch EOP from the internet.
    eop = get_iers_eop()

    # Parse the EOP from the local file.
    eop_local = read_iers_eop("./eop_IAU1980.txt")

    # The data in `eop_IAU1980.txt` is final. Hence, it must never change.
    JD = date_to_jd(2004, 4, 20, 19, 19, 19)
    @test eop.x(JD)           == eop_local.x(JD)
    @test eop.y(JD)           == eop_local.y(JD)
    @test eop.UT1_UTC(JD)     == eop_local.UT1_UTC(JD)
    @test eop.LOD(JD)         == eop_local.LOD(JD)
    @test eop.dPsi(JD)        == eop_local.dPsi(JD)
    @test eop.dEps(JD)        == eop_local.dEps(JD)
    @test eop.x_err(JD)       == eop_local.x_err(JD)
    @test eop.y_err(JD)       == eop_local.y_err(JD)
    @test eop.UT1_UTC_err(JD) == eop_local.UT1_UTC_err(JD)
    @test eop.LOD_err(JD)     == eop_local.LOD_err(JD)
    @test eop.dPsi_err(JD)    == eop_local.dPsi_err(JD)
    @test eop.dEps_err(JD)    == eop_local.dEps_err(JD)
    @test eop.x(JD)           == eop_local.x(JD)
    @test eop.y(JD)           == eop_local.y(JD)
    @test eop.UT1_UTC(JD)     == eop_local.UT1_UTC(JD)
    @test eop.LOD(JD)         == eop_local.LOD(JD)
    @test eop.dPsi(JD)        == eop_local.dPsi(JD)
    @test eop.dEps(JD)        == eop_local.dEps(JD)

    # Show
    # ==========================================================================

    buf = IOBuffer()
    io = IOContext(buf)
    show(io, MIME"text/plain"(), eop_local)
    result = String(take!(buf))

    expected = """
        EOPData_IAU1980:
             Data │ Timespan
         ─────────┼──────────────────────────────────────────────
                x │ 2004-04-01T00:00:00 -- 2004-04-30T00:00:00
                y │ 2004-04-01T00:00:00 -- 2004-04-30T00:00:00
          UT1-UTC │ 2004-04-01T00:00:00 -- 2004-04-30T00:00:00
              LOD │ 2004-04-01T00:00:00 -- 2004-04-30T00:00:00
               dψ │ 2004-04-01T00:00:00 -- 2004-04-30T00:00:00
               dϵ │ 2004-04-01T00:00:00 -- 2004-04-30T00:00:00"""

    @test result == expected

    buf = IOBuffer()
    io = IOContext(buf, :color => true)
    show(io, MIME"text/plain"(), eop_local)
    result = String(take!(buf))

    expected = """
        EOPData_IAU1980:
        \e[1m     Data \e[90m│ \e[0m\e[1mTimespan\e[0m
        \e[90m ─────────┼──────────────────────────────────────────────\e[0m
        \e[1m        x \e[90m│ \e[0m2004-04-01T00:00:00 -- 2004-04-30T00:00:00
        \e[1m        y \e[90m│ \e[0m2004-04-01T00:00:00 -- 2004-04-30T00:00:00
        \e[1m  UT1-UTC \e[90m│ \e[0m2004-04-01T00:00:00 -- 2004-04-30T00:00:00
        \e[1m      LOD \e[90m│ \e[0m2004-04-01T00:00:00 -- 2004-04-30T00:00:00
        \e[1m       dψ \e[90m│ \e[0m2004-04-01T00:00:00 -- 2004-04-30T00:00:00
        \e[1m       dϵ \e[90m│ \e[0m2004-04-01T00:00:00 -- 2004-04-30T00:00:00"""

    @test result == expected
end

@testset "EOP data (IAU2000A)" begin
    # Compare the fetched data to the one available locally
    # ==========================================================================

    # Fetch EOP from the internet.
    eop = get_iers_eop(Val(:IAU2000A))

    # Parse the EOP from the local file.
    eop_local = read_iers_eop("./eop_IAU2000A.txt", Val(:IAU2000A))

    # The data in `eop_IAU1980.txt` is final. Hence, it must never change.
    JD = date_to_jd(2004, 4, 20, 19, 19, 19)
    @test eop.x(JD)           == eop_local.x(JD)
    @test eop.y(JD)           == eop_local.y(JD)
    @test eop.UT1_UTC(JD)     == eop_local.UT1_UTC(JD)
    @test eop.LOD(JD)         == eop_local.LOD(JD)
    @test eop.dX(JD)          == eop_local.dX(JD)
    @test eop.dY(JD)          == eop_local.dY(JD)
    @test eop.x_err(JD)       == eop_local.x_err(JD)
    @test eop.y_err(JD)       == eop_local.y_err(JD)
    @test eop.UT1_UTC_err(JD) == eop_local.UT1_UTC_err(JD)
    @test eop.LOD_err(JD)     == eop_local.LOD_err(JD)
    @test eop.dX_err(JD)      == eop_local.dX_err(JD)
    @test eop.dY_err(JD)      == eop_local.dY_err(JD)
    @test eop.x(JD)           == eop_local.x(JD)
    @test eop.y(JD)           == eop_local.y(JD)
    @test eop.UT1_UTC(JD)     == eop_local.UT1_UTC(JD)
    @test eop.LOD(JD)         == eop_local.LOD(JD)
    @test eop.dX(JD)          == eop_local.dX(JD)
    @test eop.dY(JD)          == eop_local.dY(JD)

    # Show
    # ==========================================================================

    buf = IOBuffer()
    io = IOContext(buf)
    show(io, MIME"text/plain"(), eop_local)
    result = String(take!(buf))

    expected = """
        EOPData_IAU2000A:
             Data │ Timespan
         ─────────┼──────────────────────────────────────────────
                x │ 2004-04-01T00:00:00 -- 2004-04-30T00:00:00
                y │ 2004-04-01T00:00:00 -- 2004-04-30T00:00:00
          UT1-UTC │ 2004-04-01T00:00:00 -- 2004-04-30T00:00:00
              LOD │ 2004-04-01T00:00:00 -- 2004-04-30T00:00:00
               dX │ 2004-04-01T00:00:00 -- 2004-04-30T00:00:00
               dY │ 2004-04-01T00:00:00 -- 2004-04-30T00:00:00"""

    @test result == expected

    buf = IOBuffer()
    io = IOContext(buf, :color => true)
    show(io, MIME"text/plain"(), eop_local)
    result = String(take!(buf))

    expected = """
        EOPData_IAU2000A:
        \e[1m     Data \e[90m│ \e[0m\e[1mTimespan\e[0m
        \e[90m ─────────┼──────────────────────────────────────────────\e[0m
        \e[1m        x \e[90m│ \e[0m2004-04-01T00:00:00 -- 2004-04-30T00:00:00
        \e[1m        y \e[90m│ \e[0m2004-04-01T00:00:00 -- 2004-04-30T00:00:00
        \e[1m  UT1-UTC \e[90m│ \e[0m2004-04-01T00:00:00 -- 2004-04-30T00:00:00
        \e[1m      LOD \e[90m│ \e[0m2004-04-01T00:00:00 -- 2004-04-30T00:00:00
        \e[1m       dX \e[90m│ \e[0m2004-04-01T00:00:00 -- 2004-04-30T00:00:00
        \e[1m       dY \e[90m│ \e[0m2004-04-01T00:00:00 -- 2004-04-30T00:00:00"""

    @test result == expected
end


text:
################################################################################
#
# Convenience function used during testing.
#
################################################################################

function setup_ckks(poly_modulus_degree, primes)
    # Set the encryption scheme to CKKS.
    params = EncryptionParameters(SchemeType.ckks)

    # Set the degree of the cyclotomic polynomial to be a power of 2.
    set_poly_modulus_degree!(params, poly_modulus_degree)

    # Set the ciphertext coefficient modulus (a large integer) as the product of
    # large distinct prime numbers.
    set_coeff_modulus!(params, coeff_modulus_create(poly_modulus_degree, primes))

    # Now that all the parameters are set, we are ready to construct a SEALContext
    # object. This is a heavy class that checks the validity and properties of the
    # parameters we just set.
    context = SEALContext(params)

    # Generate the public and private keys by using the KeyGenerator, which
    # automatically generates the public and private keys for us.
    keygen = KeyGenerator(context)

    pub_key = PublicKey()
    create_public_key!(pub_key, keygen)

    pri_key = secret_key(keygen)

    relin_keys = RelinKeys()
    create_relin_keys!(relin_keys, keygen)

    galois_keys = GaloisKeys()
    create_galois_keys!(galois_keys, keygen)

    # Create an encryptor, evaluator and decryptor.
    encryptor = Encryptor(context, pub_key)
    evaluator = Evaluator(context)
    decryptor = Decryptor(context, pri_key)

    # Create a CKKSEncoder for encoding values and vectors.
    encoder = CKKSEncoder(context)

    (
        context = context,
        public_key = pub_key,
        private_key = pri_key,
        relinearization_keys = relin_keys,
        galois_keys = galois_keys,
        encryptor = encryptor,
        evaluator = evaluator,
        decryptor = decryptor,
        encoder = encoder
    )
end

function encrypt_array(arr, initial_scale, SEALParams)
    output = Array{Ciphertext}(undef, length(arr))

    for (idx, elem) in enumerate(arr)
        p = Plaintext()
        encode!(p, elem, initial_scale, SEALParams.encoder)

        c = Ciphertext()
        encrypt!(c, p, SEALParams.encryptor)

        output[idx] = c
    end

    output
end

function decrypt_array(arr, SEALParams)
    output = []

    for elem in arr
        p = Plaintext()
        decrypt!(p, elem, SEALParams.decryptor)

        m = similar(zeros(slot_count(SEALParams.encoder)))
        decode!(m, p, SEALParams.encoder)

        push!(output, m[1])
    end

    output
end

function encrypt_vector(arr, initial_scale, SEALParams)
    p = Plaintext()
    encode!(p, arr, initial_scale, SEALParams.encoder)

    c = Ciphertext()
    encrypt!(c, p, SEALParams.encryptor)

    c
end

function decrypt_vector(arr, SEALParams)
    p = Plaintext()
    decrypt!(p, arr, SEALParams.decryptor)

    m = similar(zeros(slot_count(SEALParams.encoder)))
    decode!(m, p, SEALParams.encoder)

    m
end

ext: jl


text:
NamedDims.dimnames(::Type{<:StremingContainer{T,N,SN,<:NamedDimsArray{L,T,N}}}) where {L,T,N,SN} = L

function ImageAxes.isstreamedaxis(name::Symbol, S::StreamingContainer{T,N,saxnames}) where {name,T,N,saxnames}
    return in(name, saxnames)
end


text:
function update(me::EarthMoonSystem, time::Float64, h::Float64)
    me.time = time

    update(me.moon, time)
    update(me.command_module, time, h)

    return me
end

max_stars_repo_name: rodrigolece/grupo-fisica-compu

max_issues_repo_name: rodrigolece/grupo-fisica-compu

max_forks_repo_name: rodrigolece/grupo-fisica-compu


text:
# include("../../datastructs/src/datastructs_module.jl")
include("../src/geom_module.jl")

module dbench
using ..datastructs
using Random
using BenchmarkTools
using ..geom

function create_random_points_in_circle(N, R = 10.0)

    r_ = R .* Random.rand(N)
    phi_ = (2*π) .* Random.rand(N)
    x_ = r_ .* cos.(phi_)
    y_ = r_ .* sin.(phi_)
    point_array = map(geom.MutablePoint, x_, y_)

    return point_array
end

function test_random_in_circle(N)
    R = 10.0
    point_array = create_random_points_in_circle(N, R)

    b = @benchmark geom.build_delaunay_mesh!($point_array)

    println("=== test_random_in_circle ===")
    println("N = $(N)")
    show(
        stdout,
        MIME("text/plain"),
        b )
    println("\n\n")
end

function testraw_random_in_circle(N)
    R = 10.0
    point_array = create_random_points_in_circle(N, R)
    b = geom.build_delaunay_mesh!(point_array)
    return nothing
end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using DiffEqDiffTools
using Test, LinearAlgebra

@time begin
include("finitedifftests.jl")
end


text:
function load_pointcloud(filename::AbstractString; kwargs...)
    @info "Loading pointcloud from $filename"
    (_, ext) = splitext(filename)
    if lowercase(ext) == ".h5" || lowercase(ext) == ".hdf5"
        pc = load_h5(filename; kwargs...)
    elseif lowercase(ext) == ".las" || lowercase(ext) == ".laz"
        pc = load_las(filename; kwargs...)
    else
        error("No loader for file extension $ext")
    end

    @info "Successfully loaded pointcloud from $filename"
    return pc
end

"""
    ATTRIBUTE_PATH

Dictionary mapping Julia column name to HDF5 dataset location
"""
const ATTRIBUTE_PATH = Dict{Symbol, String}()

# LAS file attributes
ATTRIBUTE_PATH[:position] = "LAS/Position" # 3D
ATTRIBUTE_PATH[:intensity] = "LAS/Intensity"
ATTRIBUTE_PATH[:gpstime] = "LAS/GpsTime"
ATTRIBUTE_PATH[:returnnumber] = "LAS/ReturnNumber"
ATTRIBUTE_PATH[:numberofreturns] = "LAS/NumberOfReturns"
ATTRIBUTE_PATH[:pointsourceid] = "LAS/PointSourceId"
ATTRIBUTE_PATH[:classification] = "LAS/Classification"
ATTRIBUTE_PATH[:color] = "LAS/Color" # RGB
ATTRIBUTE_PATH[:scananglerank] = "LAS/ScanAngleRank"
ATTRIBUTE_PATH[:scandirectionflag] = "LAS/ScanDirectionFlag"

# Some Julia-specific attributes
ATTRIBUTE_PATH[:featurevector] = "FeatureVector"
ATTRIBUTE_PATH[:featurename] = "FeatureName"
ATTRIBUTE_PATH[:flightstripid] = "FlightStripId"
ATTRIBUTE_PATH[:swathnum] = "Marine/SwathNum"

# Other attributes supported by PointTool, etc (see hpc repo: libs/PointBuffer/include/points/PointAttribute.h)
ATTRIBUTE_PATH[:kp] = "Marine/Distance"

ATTRIBUTE_PATH[:filter] = "Filter"
ATTRIBUTE_PATH[:stringid] = "StringId"

ATTRIBUTE_PATH[:floatposition] = "Float/Position" # 3D

ATTRIBUTE_PATH[:dtmheight] = "DtmHeight"
ATTRIBUTE_PATH[:clusterid] = "ClusterId"
ATTRIBUTE_PATH[:agl] = "AGL"
ATTRIBUTE_PATH[:distance] = "Distance"
ATTRIBUTE_PATH[:normal] = "Normal" # 3D
ATTRIBUTE_PATH[:tile] = "Tile" # 3D

function load_h5(filename::AbstractString; spacing = nothing, format = nothing)
    if format isa AbstractString
        return load_h5_repo2(filename, spacing = spacing, format = format)
    end

    pointcloud = h5open(filename, "r") do f_h # f_h = file handle
        attrs = NamedTuple()
        for (jlname, h5name) in ATTRIBUTE_PATH
            if exists(f_h, h5name)
                attrs = merge(attrs, NamedTuple{(jlname,)}((_read(f_h, jlname, h5name),)))
            end
        end

        if !haskey(attrs, :position)
            @warn "Did not find \"LAS/Position\" in HDF5 file"
        end

        if spacing === nothing
            return Table(attrs)
        else
            return Table(merge(attrs, (position = accelerate(attrs[:position], GridIndex; spacing = spacing),)))
        end
    end
    return pointcloud
end

function load_h5_repo2(filename::AbstractString; spacing = nothing, format::AbstractString = "")
    if format == "" || !isascii(format)
        error("Need to specify point format - see `man ExtractPoints` for help")
    end

    return h5open(filename, "r") do fh
        matrix = read(fh["Points"])
        n_points = isempty(matrix) ? 0 : size(matrix, 2)

        X_i = findall(isequal('X'), format)
        Y_i = findall(isequal('Y'), format)
        Z_i = findall(isequal('Z'), format)
        G_i = findall(isequal('G'), format)
        U_i = findall(isequal('U'), format)
        R_i = findall(isequal('R'), format)
        P_i = findall(isequal('P'), format)
        S_i = findall(isequal('S'), format)
        r_i = findall(isequal('r'), format)
        g_i = findall(isequal('g'), format)
        b_i = findall(isequal('b'), format)
        I_i = findall(isequal('I'), format)

        if isempty(X_i) || isempty(Y_i) || isempty(Z_i)
            error("The XYZ coordinates are mandatory")
        end
        if isempty(r_i) != isempty(g_i) || isempty(r_i) != isempty(b_i)
            error("Either all of 'rgb' or none of 'rgb' must be included in the format string")
        end

        n_attrs = 3
        nt = NamedTuple()

        if n_points == 0
            if spacing === nothing
                nt = (position = SVector{3, Float64}[],)
            else
                nt = (position = accelerate(SVector{3, Float64}[], GridIndex; spacing = spacing),)
            end

            if !isempty(G_i)
                nt = merge(nt, (agl = Float64[],))
                n_attrs += 1
            end
            if !isempty(U_i)
                nt = merge(nt, (clusterid = Int64[],))
                n_attrs += 1
            end
            if !isempty(R_i)
                nt = merge(nt, (returnnumber = Int64[],))
                n_attrs += 1
            end
            if !isempty(P_i)
                nt = merge(nt, (numberofreturns = Int64[],))
                n_attrs += 1
            end
            if !isempty(S_i)
                nt = merge(nt, (flightstripid = Int64[],))
                n_attrs += 1
            end
            if !isempty(I_i)
                nt = merge(nt, (intensity = Int64[],))
                n_attrs += 1
            end
            if !isempty(r_i)
                nt = merge(nt, (color = RGB{N0f8}[],))
                n_attrs += 3
            end
        else
            Xi = X_i[]
            Yi = Y_i[]
            Zi = Z_i[]

            pos = [@inbounds(SVector(matrix[Xi, i], matrix[Yi, i], matrix[Zi, i])) for i in 1:n_points]

            if spacing === nothing
                nt = (position = pos,)
            else
                nt = (position = accelerate(pos, GridIndex; spacing = spacing),)
            end

            if !isempty(G_i)
                Gi = G_i[]
                nt = merge(nt, (agl = [@inbounds(matrix(Zi, i) - matrix(Gi, i)) for i in 1:n_points],))
                n_attrs += 1
            end

            if !isempty(U_i)
                Ui = U_i[]
                nt = merge(nt, (clusterid = Int64[@inbounds(matrix[Ui, i]) for i in 1:n_points],))
                n_attrs += 1
            end

            if !isempty(R_i)
                Ri = R_i[]
                nt = merge(nt, (returnnumber = Int64[@inbounds(matrix[Ri, i]) for i in 1:n_points],))
                n_attrs += 1
            end

            if !isempty(P_i)
                Pi = P_i[]
                nt = merge(nt, (numberofreturns = Int64[@inbounds(matrix[Pi, i]) for i in 1:n_points],))
                n_attrs += 1
            end

            if !isempty(S_i)
                Si = S_i[]
                nt = merge(nt, (flightstripid = Int64[@inbounds(matrix[Si, i]) for i in 1:n_points],))
                n_attrs += 1
            end

            if !isempty(r_i)
                ri = r_i[]
                gi = g_i[]
                bi = b_i[]

                nt = merge(nt, (color = RGB{N0f8}[@inbounds(RGB(reinterpret(N0f8, convert(UInt8, matrix[ri, i])),
                                                                reinterpret(N0f8, convert(UInt8, matrix[gi, i])),
                                                                reinterpret(N0f8, convert(UInt8, matrix[bi, i])))) for i in 1:n_points],))

               n_attrs += 3
            end

            if !isempty(I_i)
                Ii = I_i[]
                nt = merge(nt, (intensity = Int64[@inbounds(matrix[Ii, i]) for i in 1:n_points],))
                n_attrs += 1
            end
        end

        if n_attrs != length(format)
            @warn "Some attribues in $format were ignored"
        end

        return Table(nt)
    end
end

function _read(f_h, jlname::Symbol, h5name::String)
    hdf5_to_data(read(f_h[h5name]), Val(jlname))
end

hdf5_to_data(h5data::AbstractVector, ::Val) = h5data

function hdf5_to_data(h5data::AbstractMatrix, ::Val)
    [SVector(h5data[:, i]...) for i in 1:size(h5data, 2)]
end

function hdf5_to_data(h5data::AbstractMatrix, ::Val{:position})
    @assert size(h5data, 1) == 3
    @inbounds return [SVector(h5data[1, i], h5data[2, i], h5data[3, i]) for i in 1:size(h5data, 2)]
end

function hdf5_to_data(h5data::AbstractMatrix, ::Val{:floatposition})
    @assert size(h5data, 1) == 3
    @inbounds return [SVector(h5data[1, i], h5data[2, i], h5data[3, i]) for i in 1:size(h5data, 2)]
end

function hdf5_to_data(h5data::AbstractMatrix, ::Val{:normal})
    @assert size(h5data, 1) == 3
    @inbounds return [SVector(h5data[1, i], h5data[2, i], h5data[3, i]) for i in 1:size(h5data, 2)]
end

function hdf5_to_data(h5data::AbstractMatrix, ::Val{:tile})
    @assert size(h5data, 1) == 3
    @inbounds return [SVector(h5data[1, i], h5data[2, i], h5data[3, i]) for i in 1:size(h5data, 2)]
end

function hdf5_to_data(h5data::AbstractMatrix{UInt8}, ::Val{:color})
    @assert size(h5data, 1) == 3
    @inbounds return [RGB(reinterpret(N0f8, h5data[1, i]), reinterpret(N0f8, h5data[2, i]), reinterpret(N0f8, h5data[3, i])) for i in 1:size(h5data, 2)]
end

function hdf5_to_data(h5data::AbstractMatrix{UInt16}, ::Val{:color})
    @assert size(h5data, 1) == 3
    @inbounds return [RGB(reinterpret(N0f16, h5data[1, i]), reinterpret(N0f16, h5data[2, i]), reinterpret(N0f16, h5data[3, i])) for i in 1:size(h5data, 2)]
end

function load_las(filename::AbstractString; spacing = nothing, getHeader = false)
    (header, points) = load(filename)

    pointcloud = make_table(points,
                            SVector(header.x_offset, header.y_offset, header.z_offset),
                            SVector(header.x_scale, header.y_scale, header.z_scale))

    if spacing === nothing
        if getHeader
            return pointcloud, header
        else
            return pointcloud
        end
    else
        if getHeader
            return Table(merge(columns(pointcloud), (position = accelerate(pointcloud.position, GridIndex; spacing = spacing),))), header
        else
            return Table(merge(columns(pointcloud), (position = accelerate(pointcloud.position, GridIndex; spacing = spacing),)))
        end
    end
end

function make_table(points::AbstractVector{LasPoint0}, offset, scale)
    position = map(points) do p
        @inbounds SVector(p.x*scale[1] + offset[1], p.y*scale[2] + offset[2], p.z*scale[3] + offset[3])
    end

    intensity = map(p -> p.intensity, points)
    returnnumber = map(p -> p.flag_byte & 0b00000111, points)
    numberofreturns = map(p -> (p.flag_byte & 0b00111000) >> 3, points)
    classification = map(p -> p.raw_classification, points)
    pointsourceid = map(p -> p.pt_src_id, points)

    return Table(position = position,
                 intensity = intensity,
                 returnnumber = returnnumber,
                 numberofreturns = numberofreturns,
                 classification = classification,
                 pointsourceid = pointsourceid)
end

function make_table(points::AbstractVector{LasPoint1}, offset, scale)
    position = map(points) do p
        @inbounds SVector(p.x*scale[1] + offset[1], p.y*scale[2] + offset[2], p.z*scale[3] + offset[3])
    end

    intensity = map(p -> p.intensity, points)
    returnnumber = map(p -> p.flag_byte & 0b00000111, points)
    numberofreturns = map(p -> (p.flag_byte & 0b00111000) >> 3, points)
    classification = map(p -> p.raw_classification, points)
    pointsourceid = map(p -> p.pt_src_id, points)
    gpstime = map(p -> p.gps_time, points)

    return Table(position = position,
                 intensity = intensity,
                 returnnumber = returnnumber,
                 numberofreturns = numberofreturns,
                 classification = classification,
                 pointsourceid = pointsourceid,
                 gpstime = gpstime)
end

function make_table(points::AbstractVector{LasPoint2}, offset, scale)
    position = map(points) do p
        @inbounds SVector(p.x*scale[1] + offset[1], p.y*scale[2] + offset[2], p.z*scale[3] + offset[3])
    end

    intensity = map(p -> p.intensity, points)
    returnnumber = map(p -> p.flag_byte & 0b00000111, points)
    numberofreturns = map(p -> (p.flag_byte & 0b00111000) >> 3, points)
    classification = map(p -> p.raw_classification, points)
    pointsourceid = map(p -> p.pt_src_id, points)

    color = map(points) do p
        @inbounds RGB(p.red, p.green, p.blue)
    end

    return Table(position = position,
                 intensity = intensity,
                 returnnumber = returnnumber,
                 numberofreturns = numberofreturns,
                 classification = classification,
                 pointsourceid = pointsourceid,
                 color = color)
end

function make_table(points::AbstractVector{LasPoint3}, offset, scale)
    position = map(points) do p
        @inbounds SVector(p.x*scale[1] + offset[1], p.y*scale[2] + offset[2], p.z*scale[3] + offset[3])
    end

    intensity = map(p -> p.intensity, points)
    returnnumber = map(p -> p.flag_byte & 0b00000111, points)
    numberofreturns = map(p -> (p.flag_byte & 0b00111000) >> 3, points)
    classification = map(p -> p.raw_classification, points)
    pointsourceid = map(p -> p.pt_src_id, points)
    gpstime = map(p -> p.gps_time, points)
    userdata = map(p -> p.user_data, points)

    color = map(points) do p
        @inbounds RGB(p.red, p.green, p.blue)
    end

    return Table(position = position,
                 intensity = intensity,
                 returnnumber = returnnumber,
                 numberofreturns = numberofreturns,
                 classification = classification,
                 pointsourceid = pointsourceid,
                 gpstime = gpstime,
                 userdata = userdata,
                 color = color)
end

function make_table(points::AbstractVector{LasPoint7}, offset, scale)
    position = map(points) do p
        @inbounds SVector(p.x*scale[1] + offset[1], p.y*scale[2] + offset[2], p.z*scale[3] + offset[3])
    end

    intensity = map(p -> p.intensity, points)
    returnnumber = map(p -> p.flag_byte_1 & 0b00001111, points)
    numberofreturns = map(p -> (p.flag_byte_1 & 0b11110000) >> 4, points)
    classification = map(p -> p.classification, points)
    pointsourceid = map(p -> p.pt_src_id, points)
    gpstime = map(p -> p.gps_time, points)
    userdata = map(p -> p.user_data, points)

    color = map(points) do p
        @inbounds RGB(p.red, p.green, p.blue)
    end

    return Table(position = position,
                 intensity = intensity,
                 returnnumber = returnnumber,
                 numberofreturns = numberofreturns,
                 classification = classification,
                 pointsourceid = pointsourceid,
                 gpstime = gpstime,
                 userdata = userdata,
                 color = color)
end

# Saving

function save_pointcloud(filename::String, pc::AbstractVector{<:NamedTuple}; kwargs...)
    @info "Saving pointcloud to $filename"
    (_, ext) = splitext(filename)
    if lowercase(ext) == ".h5" || lowercase(ext) == ".hdf5"
        pc = save_h5(filename, pc; kwargs...)
    elseif lowercase(ext) == ".las" || lowercase(ext) == ".laz"
        pc = save_las(filename, pc; kwargs...)
    else
        error("No loader for file extension $ext")
    end

    @info "Successfully saved pointcloud to $filename"
    return pc
end

function save_h5(filename::AbstractString, pc::AbstractVector{<:NamedTuple}; format = nothing)
    if format isa AbstractString
        return save_h5_repo2(filename, pc; format = format)
    end

    h5open(filename, "w") do f_h # f_h = file handle
        cols = columns(pc)
        for jlname in propertynames(cols)
            if !haskey(ATTRIBUTE_PATH, jlname)
                @info "Not saving column: $jlname"
                continue
            end
            h5name = ATTRIBUTE_PATH[jlname]
            f_h[h5name] = data_to_h5(getproperty(cols, jlname), Val(jlname))
        end
    end
end

function save_h5_repo2(filename::AbstractString, pc::AbstractVector{<:NamedTuple}; format::AbstractString = "")
    if format == "" || !isascii(format)
        error("Need to specify point format - see `man ExtractPoints` for help")
    end

    h5open(filename, "w") do f_h
        points = Matrix{Float64}(undef, (length(format), length(pc)))
        n_attrs = 0

        for i in 1:length(format)
            if format[i] == 'X'
                points[i,:] .= (p -> p[1]).(pc.position)
                n_attrs += 1
            elseif format[i] == 'Y'
                points[i,:] .= (p -> p[2]).(pc.position)
                n_attrs += 1
            elseif format[i] == 'Z'
                points[i,:] .= (p -> p[3]).(pc.position)
                n_attrs += 1
            elseif format[i] == 'G'
                points[i,:] .= ((agl,p)->p[3]-agl).(pc.agl, pc.position)
                n_attrs += 1
            elseif format[i] == 'U'
                points[i,:] .= pc.clusterid
                n_attrs += 1
            elseif format[i] == 'R'
                points[i,:] .= pc.returnnumber
                n_attrs += 1
            elseif format[i] == 'P'
                points[i,:] .= pc.numberofreturns
                n_attrs += 1
            elseif format[i] == 'S'
                points[i,:] .= pc.flightstripid
                n_attrs += 1
            elseif format[i] == 'r'
                points[i, :] .= round.(255.0 .* red.(pc.color))   # Save 8 bit colors
                n_attrs += 1
            elseif format[i] == 'g'
                points[i, :] .= round.(255.0 .* green.(pc.color)) # Save 8 bit colors
                n_attrs += 1
            elseif format[i] == 'b'
                points[i, :] .= round.(255.0 .* blue.(pc.color))  # Save 8 bit colors
                n_attrs += 1
            elseif format[i] == 'I'
                points[i,:] .= pc.intensity
                n_attrs += 1
            end
        end

        if n_attrs != length(format)
            @warn "Some attribues in $format were ignored"
        end

        f_h["Points"] = points
    end
end

data_to_h5(v::AbstractVector{T}, ::Val) where {T} = convert(Vector{T}, v)

function data_to_h5(v::AbstractVector{<:StaticVector{n}}, ::Val) where {n}
    return [@inbounds(v[i][j]) for j in 1:n, i in firstindex(v):lastindex(v)]
end

function data_to_h5(v::AbstractVector{RGB{N0f8}}, ::Val) where {n}
    return [reinterpret(UInt8, getfield(@inbounds(v[i]), channel)) for channel in [:r, :g, :b], i in firstindex(v):lastindex(v)]
end

function data_to_h5(v::AbstractVector{RGB{N0f16}}, ::Val) where {n}
    return [reinterpret(UInt16, getfield(@inbounds(v[i]), channel)) for channel in [:r, :g, :b], i in firstindex(v):lastindex(v)]
end

function save_las(filename::AbstractString, pc::AbstractVector{<:NamedTuple}; x_scale = 0.001, y_scale = 0.001, z_scale = 0.001, global_encoding = UInt16(0), variable_length_records=[])
    # Get spatial bounds
    box = boundingbox(getproperty(:position).(pc))
    x_min = box.xmin
    y_min = box.ymin
    z_min = box.zmin
    x_max = box.xmax
    y_max = box.ymax
    z_max = box.zmax
    x_offset = x_scale * div((x_min + x_max) / 2, x_scale)
    y_offset = y_scale * div((y_min + y_max) / 2, y_scale)
    z_offset = z_scale * div((z_min + z_max) / 2, z_scale)

    records_count = Int32(length(pc))
    point_return_count = UInt32[0,0,0,0,0]
    if isempty(pc)
        data_format_id = 0
        data_record_length = 20
        data = LasPoint0[]
    else
        if haskey(first(pc), :gpstime)
            if haskey(first(pc), :color)
                if haskey(first(pc), :userdata) # Though user data is also present in point type 3 let's ue it as a differentiater for now
                    data_format_id = 7
                    data_record_length = 36 # in bytes
                    data = Vector{LasPoint7}(undef, records_count)
                    @inbounds for i in 1:records_count
                        p = pc[i]
                        data[i] = laspoint7(p, x_offset, y_offset, z_offset, x_scale, y_scale, z_scale)
                        if haskey(p, :returnnumber)
                            point_return_count[min(p.returnnumber, length(point_return_count))] += 1
                        else
                            point_return_count[1] += 1
                        end
                    end
                else
                    data_format_id = 3
                    data_record_length = 34
                    data = Vector{LasPoint3}(undef, records_count)
                    @inbounds for i in 1:records_count
                        p = pc[i]
                        data[i] = laspoint3(p, x_offset, y_offset, z_offset, x_scale, y_scale, z_scale)
                        if haskey(p, :returnnumber)
                            point_return_count[p.returnnumber] += 1
                        else
                            point_return_count[1] += 1
                        end
                    end
                end
            else
                data_format_id = 1
                data_record_length = 28
                data = Vector{LasPoint1}(undef, records_count)
                @inbounds for i in 1:records_count
                    p = pc[i]
                    data[i] = laspoint1(p, x_offset, y_offset, z_offset, x_scale, y_scale, z_scale)
                    if haskey(p, :returnnumber)
                        point_return_count[p.returnnumber] += 1
                    else
                        point_return_count[1] += 1
                    end
                end
            end
        else
            if haskey(first(pc), :color)
                data_format_id = 2
                data_record_length = 26
                data = Vector{LasPoint2}(undef, records_count)
                @inbounds for i in 1:records_count
                    p = pc[i]
                    data[i] = laspoint2(p, x_offset, y_offset, z_offset, x_scale, y_scale, z_scale)
                    if haskey(p, :returnnumber)
                        point_return_count[p.returnnumber] += 1
                    else
                        point_return_count[1] += 1
                    end
                end
            else
                data_format_id = 0
                data_record_length = 20
                data = Vector{LasPoint0}(undef, records_count)
                @inbounds for i in 1:records_count
                    p = pc[i]
                    data[i] = laspoint0(p, x_offset, y_offset, z_offset, x_scale, y_scale, z_scale)
                    if haskey(p, :returnnumber)
                        point_return_count[p.returnnumber] += 1
                    else
                        point_return_count[1] += 1
                    end
                end
            end
        end
    end
    version = data_format_id>3 ? 4 : 2
    header_size = UInt16(data_format_id>3 ? 375 : 227)
    t = now()
    header = LasIO.LasHeader(UInt16(0), # file_source_id
                             global_encoding, # global_encoding
                             UInt32(0), # guid_1
                             UInt16(0), # guid_2
                             UInt16(0), # guid_3
                             "", # guid_4
                             UInt8(1), # LAS version major
                             UInt8(version), # LAS version minor
                             "OTHER", # System
                             "ROAMES Julia writer", # Software
                             UInt16(dayofyear(t)), # creation_year
                             UInt16(year(t)), # creation_dayofyear
                             header_size, # header_size
                             UInt32(header_size+(!isempty(variable_length_records) ? sum(sizeof, variable_length_records) : 0)), # data_offset
                             UInt32(length(variable_length_records)), # n_vlr
                             UInt8(data_format_id), # data_format_id
                             UInt16(data_record_length), # data_format_id
                             UInt32(records_count),
                             point_return_count,
                             x_scale,
                             y_scale,
                             z_scale,
                             x_offset,
                             y_offset,
                             z_offset,
                             x_max,
                             x_min,
                             y_max,
                             y_min,
                             z_max,
                             z_min,
                             UInt64(0),
                             UInt64(0),
                             UInt32(0),
                             UInt64(0),
                             zeros(UInt64,15),
                             variable_length_records,
                             []) #user_defined_bytes

    save(filename, header, data)
end

function laspoint7(p::NamedTuple, x_offset, y_offset, z_offset, x_scale, y_scale, z_scale)
    position = p.position
    @inbounds x = round(Int32, (position[1] - x_offset) / x_scale)
    @inbounds y = round(Int32, (position[2] - y_offset) / y_scale)
    @inbounds z = round(Int32, (position[3] - z_offset) / z_scale)
    intensity = haskey(p, :intensity) ? convert(UInt16, p.intensity) : 0x0000
    flagbyte1 = (haskey(p, :numberofreturns) ? convert(UInt8, p.numberofreturns) << 4 : 0x10) | (haskey(p, :returnnumber) ? convert(UInt8, p.returnnumber) : 0x01)
    classification = haskey(p, :classification) ? convert(UInt8, p.classification) : 0x00
    scan_angle = 0x00
    user_data = haskey(p, :userdata) ? p.userdata : 0x00
    pt_src_id = haskey(p, :pointsourceid) ? convert(UInt16, p.pointsourceid) : 0x0000
    gps_time = haskey(p, :gpstime) ? convert(Float64, p.gpstime) : 0.0
    red = haskey(p, :color) ? convert(N0f16, p.color.r) : 0N0f16
    green = haskey(p, :color) ? convert(N0f16, p.color.g) : 0N0f16
    blue = haskey(p, :color) ? convert(N0f16, p.color.b) : 0N0f16

    return LasPoint7(x, y, z, intensity, flagbyte1, 0x00, classification, user_data, scan_angle, pt_src_id, gps_time, red, green, blue)
end

function laspoint3(p::NamedTuple, x_offset, y_offset, z_offset, x_scale, y_scale, z_scale)
    position = p.position
    @inbounds x = round(Int32, (position[1] - x_offset) / x_scale)
    @inbounds y = round(Int32, (position[2] - y_offset) / y_scale)
    @inbounds z = round(Int32, (position[3] - z_offset) / z_scale)
    intensity = haskey(p, :intensity) ? convert(UInt16, p.intensity) : 0x0000
    flagbyte = (haskey(p, :numberofreturns) ? convert(UInt8, p.numberofreturns) << 3 : 0x08) | (haskey(p, :returnnumber) ? convert(UInt8, p.returnnumber) : 0x01)
    raw_classification = haskey(p, :classification) ? convert(UInt8, p.classification) : 0x00
    scan_angle = 0x00
    user_data = haskey(p, :userdata) ? p.userdata : 0x00
    pt_src_id = haskey(p, :pointsourceid) ? convert(UInt16, p.pointsourceid) : 0x0000
    gps_time = haskey(p, :gpstime) ? convert(Float64, p.gpstime) : 0.0
    red = haskey(p, :color) ? convert(N0f16, p.color.r) : 0N0f16
    green = haskey(p, :color) ? convert(N0f16, p.color.g) : 0N0f16
    blue = haskey(p, :color) ? convert(N0f16, p.color.b) : 0N0f16

    return LasPoint3(x, y, z, intensity, flagbyte, raw_classification, scan_angle, user_data, pt_src_id, gps_time, red, green, blue)
end

function laspoint2(p::NamedTuple, x_offset, y_offset, z_offset, x_scale, y_scale, z_scale)
    position = p.position
    @inbounds x = round(Int32, (position[1] - x_offset) / x_scale)
    @inbounds y = round(Int32, (position[2] - y_offset) / y_scale)
    @inbounds z = round(Int32, (position[3] - z_offset) / z_scale)
    intensity = haskey(p, :intensity) ? convert(UInt16, p.intensity) : 0x0000
    flagbyte = (haskey(p, :numberofreturns) ? convert(UInt8, p.numberofreturns) << 3 : 0x08) | (haskey(p, :returnnumber) ? convert(UInt8, p.returnnumber) : 0x01)
    raw_classification = haskey(p, :classification) ? convert(UInt8, p.classification) : 0x00
    scan_angle = 0x00
    user_data = 0x00
    pt_src_id = haskey(p, :pointsourceid) ? convert(UInt16, p.pointsourceid) : 0x0000
    red = haskey(p, :color) ? convert(N0f16, p.color.r) : 0N0f16
    green = haskey(p, :color) ? convert(N0f16, p.color.g) : 0N0f16
    blue = haskey(p, :color) ? convert(N0f16, p.color.b) : 0N0f16

    return LasPoint2(x, y, z, intensity, flagbyte, raw_classification, scan_angle, user_data, pt_src_id, red, green, blue)
end

function laspoint1(p::NamedTuple, x_offset, y_offset, z_offset, x_scale, y_scale, z_scale)
    position = p.position
    @inbounds x = round(Int32, (position[1] - x_offset) / x_scale)
    @inbounds y = round(Int32, (position[2] - y_offset) / y_scale)
    @inbounds z = round(Int32, (position[3] - z_offset) / z_scale)
    intensity = haskey(p, :intensity) ? convert(UInt16, p.intensity) : 0x0000
    flagbyte = (haskey(p, :numberofreturns) ? convert(UInt8, p.numberofreturns) << 3 : 0x08) | (haskey(p, :returnnumber) ? convert(UInt8, p.returnnumber) : 0x01)
    raw_classification = haskey(p, :classification) ? convert(UInt8, p.classification) : 0x00
    scan_angle = 0x00
    user_data = 0x00
    pt_src_id = haskey(p, :pointsourceid) ? convert(UInt16, p.pointsourceid) : 0x0000
    gps_time = haskey(p, :gpstime) ? convert(Float64, p.gpstime) : 0.0

    return LasPoint1(x, y, z, intensity, flagbyte, raw_classification, scan_angle, user_data, pt_src_id, gps_time)
end

function laspoint0(p::NamedTuple, x_offset, y_offset, z_offset, x_scale, y_scale, z_scale)
    position = p.position
    @inbounds x = round(Int32, (position[1] - x_offset) / x_scale)
    @inbounds y = round(Int32, (position[2] - y_offset) / y_scale)
    @inbounds z = round(Int32, (position[3] - z_offset) / z_scale)
    intensity = haskey(p, :intensity) ? convert(UInt16, p.intensity) : 0x0000
    flagbyte = (haskey(p, :numberofreturns) ? convert(UInt8, p.numberofreturns) << 3 : 0x08) | (haskey(p, :returnnumber) ? convert(UInt8, p.returnnumber) : 0x01)
    raw_classification = haskey(p, :classification) ? convert(UInt8, p.classification) : 0x00
    scan_angle = 0x00
    user_data = 0x00
    pt_src_id = haskey(p, :pointsourceid) ? convert(UInt16, p.pointsourceid) : 0x0000

    return LasPoint0(x, y, z, intensity, flagbyte, raw_classification, scan_angle, user_data, pt_src_id)
end


text:
import Base: length, eltype
using Base: HasLength, HasEltype
using Phylo.API

# For a single tree
import Base.IteratorSize
function IteratorSize(::Type{<: AbstractTree{OneTree}})
    return HasLength()
end

import Base.IteratorEltype
function IteratorEltype(::Type{<: AbstractTree{OneTree}})
    return HasEltype()
end
length(::AbstractTree{OneTree}) = 1
eltype(ni::T) where T <: AbstractTree = T

abstract type AbstractTreeIterator{T <: AbstractTree} end

function IteratorSize(::Type{AbstractTreeIterator})
    return HasLength()
end

function IteratorEltype(::Type{AbstractTreeIterator})
    return HasEltype()
end

abstract type AbstractNodeIterator{T <: AbstractTree} <: AbstractTreeIterator{T} end

function length(ni::It) where It <: AbstractNodeIterator
    return ni.filterfn === nothing ? _nnodes(ni.tree) :
        count(val -> ni.filterfn(ni.tree, _getnode(ni.tree, val)), ni)
end

abstract type AbstractBranchIterator{T <: AbstractTree} <: AbstractTreeIterator{T} end

function length(bi::It) where It <: AbstractBranchIterator
    return bi.filterfn === nothing ? _nbranches(bi.tree) :
        count(val -> bi.filterfn(bi.tree, _getbranch(bi.tree, val)), bi)
end

struct NodeIterator{T <: AbstractTree} <: AbstractNodeIterator{T}
    tree::T
    filterfn::Union{Function, Nothing}
end
"""
    nodeiter(tree::AbstractTree)

Returns an iterator over the nodes of any tree.
"""
nodeiter(tree::T) where T <: AbstractTree = NodeIterator{T}(tree, nothing)

"""
    nodefilter(filterfn::Function, tree::AbstractTree)

Returns an iterator over the nodes of any tree, where the
`AbstractNode` is filtered by the function `filterfn`.
"""
nodefilter(filterfn::Function, tree::T) where T <: AbstractTree =
    NodeIterator{T}(tree, filterfn)

eltype(ni::NodeIterator{T}) where T <: AbstractTree = nodetype(T)

struct NodeNameIterator{T <: AbstractTree} <: AbstractNodeIterator{T}
    tree::T
    filterfn::Union{Function, Nothing}
end

"""
    nodenameiter(tree::AbstractTree)

Returns an iterator over the names of the nodes of any tree.
"""
nodenameiter(tree::T) where T <: AbstractTree =
    NodeNameIterator{T}(tree, nothing)

"""
    nodenamefilter(filterfn::Function, tree::AbstractTree)

Returns an iterator over the nodenames of any tree, where the
`AbstractNode` itself is filtered by the function `filterfn`.
"""
nodenamefilter(filterfn::Function, tree::T) where T <: AbstractTree =
    NodeNameIterator{T}(tree, filterfn)

eltype(ni::NodeNameIterator{T}) where T <: AbstractTree = nodenametype(T)

struct BranchIterator{T <: AbstractTree} <: AbstractBranchIterator{T}
    tree::T
    filterfn::Union{Function, Nothing}
end

"""
    branchiter(tree::AbstractTree)

Returns an iterator over the branches of any tree.
"""
branchiter(tree::T) where T <: AbstractTree =
    BranchIterator{T}(tree, nothing)

"""
    branchfilter(filterfn::Function, tree::AbstractTree)

Returns an iterator over the branches of any tree, where the
`AbstractBranch` is filtered by the function `filterfn`.
"""
branchfilter(filterfn::Function, tree::T) where T <: AbstractTree =
    BranchIterator{T}(tree, filterfn)

eltype(bi::BranchIterator{T}) where T <: AbstractTree = branchtype(T)

struct BranchNameIterator{T <: AbstractTree} <: AbstractBranchIterator{T}
    tree::T
    filterfn::Union{Function, Nothing}
end

"""
    branchnameiter(tree::AbstractTree)

Returns an iterator over the names of branches of any tree.
"""
branchnameiter(tree::T) where T <: AbstractTree =
    BranchNameIterator{T}(tree, nothing)

"""
    branchnamefilter(filterfn::Function, tree::AbstractTree)

Returns an iterator over the names of the branches of any tree, where
the `AbstractBranch` is filtered by the function `filterfn`.
"""
branchnamefilter(filterfn::Function, tree::T) where T <: AbstractTree =
    BranchNameIterator{T}(tree, filterfn)

eltype(bi::BranchNameIterator{T}) where T <: AbstractTree = branchnametype(T)

import Base: iterate
function iterate(tree::AbstractTree, state = nothing)
    if state === nothing
        return first(gettrees(tree)), 1
    elseif ntrees(tree) > state
        return collect(gettrees(tree))[state], state + 1
    else
        return nothing
    end
end

function iterate(ni::NodeIterator, state = nothing)
    nodes = _getnodes(ni.tree)
    if state === nothing
        result = iterate(nodes)
    else
        result = iterate(nodes, state)
    end

    result === nothing && return nothing

    if ni.filterfn === nothing
        return _getnode(ni.tree, result[1]), result[2]
    end

    val, state = result
    node = _getnode(ni.tree, val)
    while !ni.filterfn(ni.tree, node)
        result = iterate(nodes, state)
        result === nothing && return nothing
        val, state = result
        node = _getnode(ni.tree, val)
    end

    return node, state
end

function iterate(ni::NodeNameIterator, state = nothing)
    nodes = _getnodes(ni.tree)
    if state === nothing
        result = iterate(nodes)
    else
        result = iterate(nodes, state)
    end

    result === nothing && return nothing

    if ni.filterfn === nothing
        return _getnodename(ni.tree, result[1]), result[2]
    end

    val, state = result
    node = _getnode(ni.tree, val)
    while !ni.filterfn(ni.tree, node)
        result = iterate(nodes, state)
        result === nothing && return nothing
        val, state = result
        node = _getnode(ni.tree, val)
    end

    name = _getnodename(ni.tree, val)
    return name, state
end

function iterate(bi::BranchIterator, state = nothing)
    branches = _getbranches(bi.tree)
    if state === nothing
        result = iterate(branches)
    else
        result = iterate(branches, state)
    end

    result === nothing && return nothing

    if bi.filterfn === nothing
        return _getbranch(bi.tree, result[1]), result[2]
    end

    val, state = result
    branch = _getbranch(bi.tree, val)
    while !bi.filterfn(bi.tree, branch)
        result = iterate(branches, state)
        result === nothing && return nothing
        val, state = result
        branch = _getbranch(bi.tree, val)
    end

    return branch, state
end

function iterate(bi::BranchNameIterator, state = nothing)
    branches = _getbranches(bi.tree)
    if state === nothing
        result = iterate(branches)
    else
        result = iterate(branches, state)
    end

    result === nothing && return nothing

    if bi.filterfn === nothing
        return _getbranchname(bi.tree, result[1]), result[2]
    end

    val, state = result
    branch = _getbranch(bi.tree, val)
    while !bi.filterfn(bi.tree, branch)
        result = iterate(branches, state)
        result === nothing && return nothing
        val, state = result
        branch = _getbranch(bi.tree, val)
    end
    name = _getbranchname(bi.tree, val)

    return name, state
end


text:
#src # Adding a new equation: scalar conservation laws

# If you want to use Trixi for your own research, you might be interested in
# a new physics model that's not already included in Trixi.jl. In this tutorial,
# we will implement the cubic conservation law
# ```math
# \partial_t u(t,x) + \partial_x u(t,x)^3 = 0
# ```
# in a periodic domain in one space dimension. In Trixi.jl, such a mathematical model
# is encoded as a subtype of [`Trixi.AbstractEquations`](@ref).

# ## Basic setup
using Trixi

struct CubicEquation <: Trixi.AbstractEquations{1 #= number of spatial dimensions =#,
                                                1 #= number of primary variables, i.e. scalar =#};
end

# We create `CubicEquation` as an empty `struct` since we do not use any parameters
# for this equation. Other models could bundle arbitrary parameters, e.g., the
# ideal gas constant for the compressible Euler equations.

# Next, we define the physical flux `f(u) = u^3` using the calling structure
# used in Trixi.jl.

Trixi.flux(u, orientation, equation::CubicEquation) = u.^3
Trixi.varnames(_, ::CubicEquation) = ("scalar",)

# In Trixi.jl, the conserved variables `u` are usually passed as `SVector`s of variables
# at a single physical location. Hence, we must use `u.^3` instead of the scalar
# operation `u^3`.

# That's already enough to run a simple simulation with a standard DGSEM discretization
# using the non-dissipative central flux at interfaces.

using OrdinaryDiffEq

## Create a simulation setup
equation = CubicEquation()

initial_condition_sine(x, t, equation::CubicEquation) = SVector(sinpi(x[1]))

mesh = TreeMesh(-1.0, 1.0, # min/max coordinates
                initial_refinement_level=4,
                n_cells_max=10^4)

solver = DGSEM(3 #= polynomial degree =#, flux_central)

semi = SemidiscretizationHyperbolic(mesh, equation, initial_condition_sine, solver)

# We wrap the return value of the `initial_condition_sine` inside an `SVector` since that's the approach
# used in Trixi.jl also for systems of equations. We need to index the spatial coordinate `x[1]`,
# since it is an `SVector` with one component. In multiple space dimensions, all spatial coordinates
# are passed together.

# Next, we create an `ODEProblem` from the SciML/DifferentialEquations ecosystem.
# We can solve this ODE numerically using any time integration method,
# e.g., `SSPRK43` from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl).
# Before, we set up a [callback](@ref callbacks-id) to summarize the simulation setup.

## Create ODE problem with given time span
tspan = (0.0, 0.09)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
callbacks = CallbackSet(summary_callback)

## OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = solve(ode, SSPRK43(),
            save_everystep=false, callback=callbacks);

# That's it, you ran your first simulation using your new equation with Trixi! Now, we can plot
# the solution at the final time using Plots.jl.
using Plots
plot(sol)

# You can already see that discontinuities will develop and oscillations start to
# occur around steep parts of the wave. That's expected from our central discretization.
# To avoid these issues, we need to use dissipative numerical fluxes (approximate
# Riemann solvers) at interfaces.

# ## Advanced setup

# Thus, we add a Godunov's flux for our cubic equation. That is easy for this equation
# since the wave speed `f'(u) = 3u^2` is always non-negative.

@inline Trixi.flux_godunov(u_ll, u_rr, orientation, equation::CubicEquation) = flux(u_ll, orientation, equation)

# Let's run the example again but with a dissipative numerical flux at interfaces.
# `remake` will recreate the semidiscretization we used before and only change
# selected parameters, in this case the `solver`.

## A new setup with dissipation
semi = remake(semi, solver=DGSEM(3, flux_godunov))
ode = semidiscretize(semi, tspan)
sol = solve(ode, SSPRK43(), save_everystep=false)
plot!(sol)

# You can see that there are fewer oscillations, in particular around steep edges.
# Now let's increase the final time (and also the spatial resolution).

## A larger final time: Nonclassical shocks develop (you can even increase the refinement to 12)
semi = remake(semi, mesh=TreeMesh(-1.0, 1.0, initial_refinement_level=8, n_cells_max=10^5))
ode = semidiscretize(semi, (0.0, 0.5) #= tspan =#)
sol = solve(ode, SSPRK43(), save_everystep=false)
plot(sol)

# You can observe that nonclassical shocks develop and are stable under grid refinement,
# e.g. for `initial_refinement_level=12`. In this case, these nonclassical shocks
# can be avoided by using an entropy-dissipative semidiscretization. Thus, we need
# to define an entropy-conservative numerical flux

@inline function Trixi.flux_ec(u_ll, u_rr, orientation, equation::CubicEquation)
  return SVector(0.25 * (u_ll[1]^3 + u_ll[1]^2 * u_rr[1] + u_ll[1] * u_rr[1]^2 + u_rr[1]^3))
end

# and use a [`VolumeIntegralFluxDifferencing`](@ref) instead of the standard
# [`VolumeIntegralWeakForm`](@ref) in the DGSEM.

## Let's use a provably entropy-dissipative semidiscretization
semi = remake(semi, solver=DGSEM(3, flux_godunov, VolumeIntegralFluxDifferencing(flux_ec)))
ode = semidiscretize(semi, (0.0, 0.5))
sol = solve(ode, SSPRK43(), save_everystep=false);
plot(sol)

# Possible next steps could be
# - to define `Trixi.max_abs_speeds(u, equations::CubicEquation) = 3 * u[1]^2`
#   to use CFL-based time step control via a [`StepsizeCallback`](@ref)
# - to define quantities of interest like `Trixi.entropy(u, equations::CubicEquation) = u[1]^2`
#   and integrate them in a simulation using the [`AnalysisCallback`](@ref)
# - to experiment with shock-capturing volume integrals [`VolumeIntegralShockCapturingHG`](@ref)
#   and adaptive mesh refinement [`AMRCallback`](@ref)

# ## Summary of the code

# To sum up, here is the complete code that we used (without the callbacks since these create a
# lot of unnecessary output in the doctests of this tutorial).
# In addition, we create the `struct` inside the new module `CubicConservationLaw`. That
# ensures that we can re-create `struct`s defined therein without having to restart Julia.

## Define new physics
module CubicConservationLaw

using Trixi

struct CubicEquation <: Trixi.AbstractEquations{1 #= number of spatial dimensions =#,
                                                1 #= number of primary variables, i.e. scalar =#}
end

@inline Trixi.flux(u, orientation, equation::CubicEquation) = u.^3
Trixi.varnames(_, ::CubicEquation) = ("scalar",)

@inline Trixi.flux_godunov(u_ll, u_rr, orientation, equation::CubicEquation) = flux(u_ll, orientation, equation)
@inline function Trixi.flux_ec(u_ll, u_rr, orientation, equation::CubicEquation)
  return SVector(0.25 * (u_ll[1]^3 + u_ll[1]^2 * u_rr[1] + u_ll[1] * u_rr[1]^2 + u_rr[1]^3))
end

end # module

## Create a simulation setup
import .CubicConservationLaw
using Trixi
using OrdinaryDiffEq
using Plots

equation = CubicConservationLaw.CubicEquation()

initial_condition_sine(x, t, equation::CubicConservationLaw.CubicEquation) = SVector(sinpi(x[1]))

mesh = TreeMesh(-1.0, 1.0, # min/max coordinates
                initial_refinement_level=4,
                n_cells_max=10^4)

solver = DGSEM(3 #= polynomial degree =#, flux_central)

semi = SemidiscretizationHyperbolic(mesh, equation, initial_condition_sine, solver)

## Create ODE problem with given time span
tspan = (0.0, 0.1)
ode = semidiscretize(semi, tspan)

## OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = solve(ode, SSPRK43(), save_everystep=false)
plot(sol)

## A new setup with dissipation
semi = remake(semi, solver=DGSEM(3, flux_godunov))
ode = semidiscretize(semi, tspan)
sol = solve(ode, SSPRK43(), save_everystep=false)
plot!(sol)

## A larger final time: Nonclassical shocks develop (you can even increase the refinement to 12)
semi = remake(semi, mesh=TreeMesh(-1.0, 1.0, initial_refinement_level=8, n_cells_max=10^5))
ode = semidiscretize(semi, (0.0, 0.5))
sol = solve(ode, SSPRK43(), save_everystep=false)
plot(sol)

## Let's use a provably entropy-dissipative semidiscretization
semi = remake(semi, solver=DGSEM(3, flux_godunov, VolumeIntegralFluxDifferencing(flux_ec)))
ode = semidiscretize(semi, (0.0, 0.5))
sol = solve(ode, SSPRK43(), save_everystep=false)
plot(sol)


text:
using DSGE, JLD2, FileIO, Distributions, Nullables, SparseArrays
import Test: @test, @testset, @test
import DataStructures: OrderedDict

### Model
m = OneAssetHANK()
# Test ouput of steadystate! has not changed
#=
out = load("reference/new_steadystate_output.jld2")
@testset "steadystate!(m) output" begin
    steadystate!(m)
    N = length(m.parameters)
    if (m.steady_state[N+1].key == :V_ss)
        N += 1
    end
    for param in out["parameters"]
        @test m.parameters[m.keys[param.key]].value == param.value
    end
    for ss_param in out["steady_state"]
        @test m.steady_state[m.keys[ss_param.key] - N].value == ss_param.value
    end
    @test m.keys == out["keys"]
end
=#

# Test ouput of eqcond has not changed
out = load("reference/eqcond_output.jld2")
Γ0, Γ1, Ψ, Π, C = eqcond(m)
@testset "eqcond(m) output" begin
    @test size(Γ0) == size(out["output"][1])
    @test size(Γ1) == size(out["output"][2])
    @test typeof(Γ0) == typeof(out["output"][1])
    @test typeof(Γ1) == typeof(out["output"][2])
    @test Γ0 == out["output"][1]
    @test Γ1 ≈ out["output"][2] # Note: this does not hold with perfect equality
    @test Ψ  == out["output"][3]
    @test Π  == out["output"][4]
    @test C  == out["output"][5]
end

# Test ouput of solve has not changed
#out = load("reference/solve_output_interns.jld2")
out = load("reference/one_asset_after_upgrade.jld2")
TTT, RRR, CCC, inverse_basis = solve(m)
@testset "solve(m) output" begin
    @test TTT           = out["TTT"]
    @test RRR           = out["RRR"]
    @test CCC           = out["CCC"]
    @test inverse_basis = out["inverse_basis"]
end


text:
# Masks

const SF_FORMAT_ENDMASK   = 0x30000000
const SF_FORMAT_TYPEMASK  = 0x0FFF0000
const SF_FORMAT_SUBMASK   = 0x0000FFFF

# Endian-ness options

const SF_ENDIAN_FILE      = 0x00000000 # Default file endian-ness.
const SF_ENDIAN_LITTLE    = 0x10000000 # Force little endian-ness.
const SF_ENDIAN_BIG       = 0x20000000 # Force big endian-ness.
const SF_ENDIAN_CPU       = 0x30000000 # Force CPU endian-ness.

# Major Formats

const SF_FORMAT_WAV       = 0x00010000 # Microsoft WAV format (little endian).
const SF_FORMAT_AIFF      = 0x00020000 # Apple/SGI AIFF format (big endian).
const SF_FORMAT_AU        = 0x00030000 # Sun/NeXT AU format (big endian).
const SF_FORMAT_RAW       = 0x00040000 # RAW PCM data.
const SF_FORMAT_PAF       = 0x00050000 # Ensoniq PARIS file format.
const SF_FORMAT_SVX       = 0x00060000 # Amiga IFF / SVX8 / SV16 format.
const SF_FORMAT_NIST      = 0x00070000 # Sphere NIST format.
const SF_FORMAT_VOC       = 0x00080000 # VOC files.
const SF_FORMAT_IRCAM     = 0x000A0000 # Berkeley/IRCAM/CARL
const SF_FORMAT_W64       = 0x000B0000 # Sonic Foundry's 64 bit RIFF/WAV
const SF_FORMAT_MAT4      = 0x000C0000 # Matlab (tm) V4.2 / GNU Octave 2.0
const SF_FORMAT_MAT5      = 0x000D0000 # Matlab (tm) V5.0 / GNU Octave 2.1
const SF_FORMAT_PVF       = 0x000E0000 # Portable Voice Format
const SF_FORMAT_XI        = 0x000F0000 # Fasttracker 2 Extended Instrument
const SF_FORMAT_HTK       = 0x00100000 # HMM Tool Kit format
const SF_FORMAT_SDS       = 0x00110000 # Midi Sample Dump Standard
const SF_FORMAT_AVR       = 0x00120000 # Audio Visual Research
const SF_FORMAT_WAVEX     = 0x00130000 # MS WAVE with WAVEFORMATEX
const SF_FORMAT_SD2       = 0x00160000 # Sound Designer 2
const SF_FORMAT_FLAC      = 0x00170000 # FLAC lossless file format
const SF_FORMAT_CAF       = 0x00180000 # Core Audio File format
const SF_FORMAT_WVE       = 0x00190000 # Psion WVE format
const SF_FORMAT_OGG       = 0x00200000 # Xiph OGG container
const SF_FORMAT_MPC2K     = 0x00210000 # Akai MPC 2000 sampler
const SF_FORMAT_RF64      = 0x00220000 # RF64 WAV file

# SubFormats

const SF_FORMAT_PCM_S8    = 0x00000001 # Signed 8 bit data
const SF_FORMAT_PCM_16    = 0x00000002 # Signed 16 bit data
const SF_FORMAT_PCM_24    = 0x00000003 # Signed 24 bit data
const SF_FORMAT_PCM_32    = 0x00000004 # Signed 32 bit data
const SF_FORMAT_PCM_U8    = 0x00000005 # Unsigned 8 bit data (WAV and RAW only)
const SF_FORMAT_FLOAT     = 0x00000006 # 32 bit float data
const SF_FORMAT_DOUBLE    = 0x00000007 # 64 bit float data
const SF_FORMAT_ULAW      = 0x00000010 # U-Law encoded.
const SF_FORMAT_ALAW      = 0x00000011 # A-Law encoded.
const SF_FORMAT_IMA_ADPCM = 0x00000012 # IMA ADPCM.
const SF_FORMAT_MS_ADPCM  = 0x00000013 # Microsoft ADPCM.
const SF_FORMAT_GSM610    = 0x00000020 # GSM 6.10 encoding.
const SF_FORMAT_VOX_ADPCM = 0x00000021 # Oki Dialogic ADPCM encoding.
const SF_FORMAT_G721_32   = 0x00000030 # 32kbs G721 ADPCM encoding.
const SF_FORMAT_G723_24   = 0x00000031 # 24kbs G723 ADPCM encoding.
const SF_FORMAT_G723_40   = 0x00000032 # 40kbs G723 ADPCM encoding.
const SF_FORMAT_DWVW_12   = 0x00000040 # 12 bit Delta Width Variable Word encoding.
const SF_FORMAT_DWVW_16   = 0x00000041 # 16 bit Delta Width Variable Word encoding.
const SF_FORMAT_DWVW_24   = 0x00000042 # 24 bit Delta Width Variable Word encoding.
const SF_FORMAT_DWVW_N    = 0x00000043 # N bit Delta Width Variable Word encoding.
const SF_FORMAT_DPCM_8    = 0x00000050 # 8 bit differential PCM (XI only)
const SF_FORMAT_DPCM_16   = 0x00000051 # 16 bit differential PCM (XI only)
const SF_FORMAT_VORBIS    = 0x00000060 # Xiph Vorbis encoding.

# Library flags

# const SF_SEEK_SET = 0
# const SF_SEEK_CUR = 1
# const SF_SEEK_END = 2

const SFM_READ = Int32(0x10)
const SFM_WRITE = Int32(0x20)

const SF_SEEK_SET = Int32(0)
const SF_SEEK_CUR = Int32(1)
const SF_SEEK_END = Int32(2)

formatcode(::Type{format"WAV"}) = SF_FORMAT_WAV
formatcode(::Type{format"FLAC"}) = SF_FORMAT_FLAC
formatcode(::Type{format"OGG"}) = SF_FORMAT_OGG

subformatcode(::Type{PCM16Sample}) = SF_FORMAT_PCM_16
subformatcode(::Type{PCM32Sample}) = SF_FORMAT_PCM_32
subformatcode(::Type{Float32}) = SF_FORMAT_FLOAT
subformatcode(::Type{Float64}) = SF_FORMAT_DOUBLE

"""Take a LibSndFile format code and return a suitable sample type"""
function fmt_to_type(fmt)
    mapping = Dict{UInt32, Type}(
        SF_FORMAT_PCM_S8 => PCM16Sample,
        SF_FORMAT_PCM_U8 => PCM16Sample,
        SF_FORMAT_PCM_16 => PCM16Sample,
        SF_FORMAT_PCM_24 => PCM32Sample,
        SF_FORMAT_PCM_32 => PCM32Sample,
        SF_FORMAT_FLOAT => Float32,
        SF_FORMAT_DOUBLE => Float64,
        SF_FORMAT_VORBIS => Float32,
    )

    masked = fmt & SF_FORMAT_SUBMASK
    masked in keys(mapping) || error("Format code $masked not recognized by LibSndFile.jl")

    mapping[masked]
end

const sf_count_t = Int64

mutable struct SF_INFO
    frames::sf_count_t
    samplerate::Int32
    channels::Int32
    format::Int32
    sections::Int32
    seekable::Int32
end

SF_INFO() = SF_INFO(0, 0, 0, 0, 0, 0)

function sf_open(fname::String, mode, sfinfo)
    filePtr = ccall((:sf_open, libsndfile), Ptr{Cvoid},
                    (Cstring, Int32, Ref{SF_INFO}),
                    fname, mode, sfinfo)

    if filePtr == C_NULL
        error("LibSndFile.jl error while opening $fname: ", sf_strerror(C_NULL))
    end

    filePtr
end

# internals to get the virtual IO interface working
include("virtualio.jl")

function sf_open(io::T, mode, sfinfo) where T <: IO
    virtio = SF_VIRTUAL_IO(T)
    filePtr = ccall((:sf_open_virtual, libsndfile), Ptr{Cvoid},
                    (Ref{SF_VIRTUAL_IO}, Int32, Ref{SF_INFO}, Ptr{T}),
                    virtio, mode, sfinfo, pointer_from_objref(io))
    if filePtr == C_NULL
        error("LibSndFile.jl error while opening stream: ", sf_strerror(C_NULL))
    end

    filePtr
end

function sf_close(filePtr)
    err = ccall((:sf_close, libsndfile), Int32, (Ptr{Cvoid},), filePtr)
    if err != 0
        error("LibSndFile.jl error: Failed to close file: ", sf_strerror(filePtr))
    end
end

"""
Wrappers for the family of sf_readf_* functions, which read the given number
of frames into the given array. Returns the number of frames read.
"""
function sf_readf end

sf_readf(filePtr, dest::Array{T}, nframes) where T <: Union{Int16, PCM16Sample} =
    ccall((:sf_readf_short, libsndfile), Int64,
        (Ptr{Cvoid}, Ptr{T}, Int64),
        filePtr, dest, nframes)

sf_readf(filePtr, dest::Array{T}, nframes) where T <: Union{Int32, PCM32Sample} =
    ccall((:sf_readf_int, libsndfile), Int64,
        (Ptr{Cvoid}, Ptr{T}, Int64),
        filePtr, dest, nframes)

sf_readf(filePtr, dest::Array{Float32}, nframes) =
    ccall((:sf_readf_float, libsndfile), Int64,
        (Ptr{Cvoid}, Ptr{Float32}, Int64),
        filePtr, dest, nframes)

sf_readf(filePtr, dest::Array{Float64}, nframes) =
    ccall((:sf_readf_double, libsndfile), Int64,
        (Ptr{Cvoid}, Ptr{Float64}, Int64),
        filePtr, dest, nframes)

"""
Wrappers for the family of sf_writef_* functions, which write the given number
of frames in the source array to the file. Returns the number of frames written.
"""
function sf_writef end

sf_writef(filePtr, src::Array{T}, nframes) where T <: Union{Int16, PCM16Sample} =
    ccall((:sf_writef_short, libsndfile), Int64,
                (Ptr{Cvoid}, Ptr{T}, Int64),
                filePtr, src, nframes)

sf_writef(filePtr, src::Array{T}, nframes) where T <: Union{Int32, PCM32Sample} =
    ccall((:sf_writef_int, libsndfile), Int64,
                (Ptr{Cvoid}, Ptr{T}, Int64),
                filePtr, src, nframes)

sf_writef(filePtr, src::Array{Float32}, nframes) =
    ccall((:sf_writef_float, libsndfile), Int64,
                (Ptr{Cvoid}, Ptr{Float32}, Int64),
                filePtr, src, nframes)

sf_writef(filePtr, src::Array{Float64}, nframes) =
    ccall((:sf_writef_double, libsndfile), Int64,
                (Ptr{Cvoid}, Ptr{Float64}, Int64),
                filePtr, src, nframes)

function sf_strerror(filePtr)
    errmsg = ccall((:sf_strerror, libsndfile), Ptr{UInt8}, (Ptr{Cvoid},), filePtr)
    unsafe_string(errmsg)
end

sf_seek(filePtr, frames::sf_count_t, whence::Integer) =
    ccall((:sf_seek, libsndfile), Int64,
                (Ptr{Cvoid}, Int64, Int32),
                filePtr, frames, whence)

function version()
  SFC_GET_LIB_VERSION	= 0x1000
  buf = zeros(Cchar,256)
  v = Cstring(pointer(buf))
  ccall((:sf_command, libsndfile), Int64,
        (Ptr{Cvoid}, UInt, Cstring, UInt),
        C_NULL, SFC_GET_LIB_VERSION, v, sizeof(buf))
  unsafe_string(v)
end


text:
using ValueOrientedRiskManagementInsurance

using Distributions
using DataFrames
using Base.Test

include("S2NonLife.jl")

println("Start S2NonLife test")

@test_approx_eq(round(scr_prem_res, 2), 210.95)
@test_approx_eq(round(scr_lapse, 2), 0.0)
@test_approx_eq(round(scr_cat_fire, 2), 60.00)
@test_approx_eq(round(scr_cat_liab, 2), 239.65)
@test_approx_eq(round(scr_cat_man_made, 2), 247.05)
@test_approx_eq(round(scr_cat_other, 2), 0.0)
@test_approx_eq(round(scr_cat_nprop, 2), 0.0)
@test_approx_eq(round(scr_cat, 2), 247.05)
@test_approx_eq(round(scr_nl, 2), 362.76)

println("End S2NonLife test")


text:
function default_fallback_function(
    data::I18nData,
    default_languages::AbstractDict
)
    default_languages = Dict(
        string(key) => parse_locale_name.(value)
        for (key, value) in default_languages
    )
    if !haskey(default_languages, "_")
        default_languages["_"] = String[]
    end
    cached = Dict{String, Vector{String}}()
    (lang) -> get!(cached, lang) do
        t = split(lang, '_')
        base_lang = t[1]
        fallbacks = if length(t) > 1 && haskey(data, base_lang)
            String[base_lang]
        else
            String[]
        end
        append!(fallbacks, String[
            key for key in keys(data)
            if key != lang && contains(key, "_") && startswith(key, base_lang)
        ])
        append!(fallbacks, GlobalI18nConfig.fallback)
        append!(fallbacks, get(default_languages, lang, default_languages["_"]))
        fallbacks
    end
end

get_fallback_function(data::I18nData, fallback::AbstractString) = get_fallback_function(data, [fallback])
get_fallback_function(data::I18nData, fallback::AbstractVector{<:AbstractString}) =
    default_fallback_function(data, Dict("_" => fallback))
get_fallback_function(_::I18nData, fallback::AbstractDict) =
    default_fallback_function(data, fallback)
get_fallback_function(_::I18nData, fallback::Function) = fallback


text:
module mbas1
using FinEtools
using FinEtoolsVibInFluids
using LinearAlgebra
using Test
function test()
	n = fill(0.0, 3)
	loc = fill(0.0, 3, 1)
	J =  [0.233725  0.852953
	0.914648  0.182704
	0.309126  0.522979 ]
	# cross(J[:,1], J[:,2])
	# 	FinEtoolsVibInFluids.WetVibBEMModule.getnormal!(n, loc, J)
	true
end
end
using .mbas1
mbas1.test()

max_stars_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_issues_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_forks_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids


text:
module IPhreeqc

if Sys.iswindows()
    const Lib_IPhreeqc_path=joinpath(@__DIR__, "..", "lib", "libiphreeqc.dll")
elseif Sys.islinux()
    const Lib_IPhreeqc_path=joinpath(@__DIR__, "..", "lib", "libiphreeqc.so")
else
    error("No binary files for your OS. Only Windows and Linux at the moment!")
end

VRESULT = Cint
IPQ_RESULT = Cint
VAR = Union{Clong, Cint, Cdouble, Cstring}

function VarAllocString(pSrc)
    ccall((:VarAllocString, Lib_IPhreeqc_path), Cstring, (Cstring,), pSrc)
end

function VarClear(pvar)
    ccall((:VarClear, Lib_IPhreeqc_path), VRESULT, (Ptr{VAR},), pvar)
end

function VarCopy(pvarDest, pvarSrc)
    ccall((:VarCopy, Lib_IPhreeqc_path), VRESULT, (Ptr{VAR}, Ptr{VAR}), pvarDest, pvarSrc)
end

function VarFreeString(pSrc)
    ccall((:VarFreeString, Lib_IPhreeqc_path), Nothing, (Cstring,), pSrc)
end

function VarInit(pvar)
    ccall((:VarInit, Lib_IPhreeqc_path), Nothing, (Ptr{VAR},), pvar)
end

function AccumulateLine(id::Int, line::AbstractString)
    # ccall((:AccumulateLine, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cstring), id, line)
    IPQ_RESULT = ccall((:AccumulateLine, Lib_IPhreeqc_path), Cint, (Cint, Cstring), id, line)
end

function AddError(id::Int, error_msg)
    ccall((:AddError, Lib_IPhreeqc_path), Cint, (Cint, Cstring), id, error_msg)
end

function AddWarning(id::Int, warn_msg)
    ccall((:AddWarning, Lib_IPhreeqc_path), Cint, (Cint, Cstring), id, warn_msg)
end

function ClearAccumulatedLines(id::Int)
    ccall((:ClearAccumulatedLines, Lib_IPhreeqc_path), IPQ_RESULT, (Cint,), id)
end

function CreateIPhreeqc()
    id = ccall((:CreateIPhreeqc, Lib_IPhreeqc_path), Cint, ())
    id = convert(Int, id)
    return id
end

function DestroyIPhreeqc(id::Int)
    ccall((:DestroyIPhreeqc, Lib_IPhreeqc_path), IPQ_RESULT, (Cint,), id)
end

function GetComponent(id::Int, n::Int)
    ccall((:GetComponent, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
end

function GetComponentCount(id::Int)
    n_components = ccall((:GetComponentCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
    return convert(Int, n_components)
end

function GetCurrentSelectedOutputUserNumber(id::Int)
    ccall((:GetCurrentSelectedOutputUserNumber, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetDumpFileName(id::Int)
    ccall((:GetDumpFileName, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetDumpFileOn(id::Int)
    ccall((:GetDumpFileOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetDumpString(id::Int)
    ccall((:GetDumpString, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetDumpStringLine(id::Int, n::Int)
    ccall((:GetDumpStringLine, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
end

function GetDumpStringLineCount(id::Int)
    ccall((:GetDumpStringLineCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetDumpStringOn(id::Int)
    ccall((:GetDumpStringOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetErrorFileName(id::Int)
    ccall((:GetErrorFileName, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetErrorFileOn(id::Int)
    ccall((:GetErrorFileOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetErrorString(id::Int)
    ccall((:GetErrorString, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetErrorStringLine(id::Int, n::Int)
    ccall((:GetErrorStringLine, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
end

function GetErrorStringLineCount(id::Int)
    ccall((:GetErrorStringLineCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetErrorStringOn(id::Int)
    ccall((:GetErrorStringOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetLogFileName(id::Int)
    ccall((:GetLogFileName, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetLogFileOn(id::Int)
    ccall((:GetLogFileOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetLogString(id::Int)
    ccall((:GetLogString, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetLogStringLine(id::Int, n::Int)
    ccall((:GetLogStringLine, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
end

function GetLogStringLineCount(id::Int)
    ccall((:GetLogStringLineCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetLogStringOn(id::Int)
    ccall((:GetLogStringOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetNthSelectedOutputUserNumber(id::Int, n::Int)
    ccall((:GetNthSelectedOutputUserNumber, Lib_IPhreeqc_path), Cint, (Cint, Cint), id, n)
end

function GetOutputFileName(id::Int)
    ccall((:GetOutputFileName, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetOutputFileOn(id::Int)
    ccall((:GetOutputFileOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetOutputString(id::Int)
    res = ccall((:GetOutputString, Lib_IPhreeqc_path), Cstring, (Cint,), id)
    return unsafe_string(res)
end

function GetOutputStringLine(id::Int, n::Int)
    res = ccall((:GetOutputStringLine, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
    return unsafe_string(res)
end

function GetOutputStringLineCount(id::Int)
    ccall((:GetOutputStringLineCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetOutputStringOn(id::Int)
    ccall((:GetOutputStringOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputColumnCount(id::Int)
    ccall((:GetSelectedOutputColumnCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputCount(id::Int)
    ccall((:GetSelectedOutputCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputFileName(id::Int)
    ccall((:GetSelectedOutputFileName, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetSelectedOutputFileOn(id::Int)
    ccall((:GetSelectedOutputFileOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputRowCount(id::Int)
    ccall((:GetSelectedOutputRowCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputString(id::Int)
    ccall((:GetSelectedOutputString, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetSelectedOutputStringLine(id::Int, n::Int)
    ccall((:GetSelectedOutputStringLine, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
end

function GetSelectedOutputStringLineCount(id::Int)
    ccall((:GetSelectedOutputStringLineCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputStringOn(id::Int)
    ccall((:GetSelectedOutputStringOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputValue(id::Int, row::Int, col::Int, pVAR)
    ccall((:GetSelectedOutputValue, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint, Cint, Ptr{VAR}), id, row, col, pVAR)
end

function GetSelectedOutputValue2(id::Int, row::Int, col::Int, vtype, dvalue, svalue, svalue_length::UInt32)
    ccall((:GetSelectedOutputValue2, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint, Cint, Ptr{Cint}, Ptr{Cdouble}, Cstring, UInt32), id, row, col, vtype, dvalue, svalue, svalue_length)
end

function GetVersionString()
    ccall((:GetVersionString, Lib_IPhreeqc_path), Cstring, ())
end

function GetWarningString(id::Int)
    ccall((:GetWarningString, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetWarningStringLine(id::Int, n::Int)
    ccall((:GetWarningStringLine, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
end

function GetWarningStringLineCount(id::Int)
    ccall((:GetWarningStringLineCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function LoadDatabase(id::Int, filename)
    ccall((:LoadDatabase, Lib_IPhreeqc_path), Cint, (Cint, Cstring), id, filename)
end

function LoadDatabaseString(id::Int, input)
    ccall((:LoadDatabaseString, Lib_IPhreeqc_path), Cint, (Cint, Cstring), id, input)
end

function OutputAccumulatedLines(id::Int)
    ccall((:OutputAccumulatedLines, Lib_IPhreeqc_path), Nothing, (Cint,), id)
end

function OutputErrorString(id::Int)
    ccall((:OutputErrorString, Lib_IPhreeqc_path), Nothing, (Cint,), id)
end

function OutputWarningString(id::Int)
    ccall((:OutputWarningString, Lib_IPhreeqc_path), Nothing, (Cint,), id)
end

function RunAccumulated(id::Int)
    ccall((:RunAccumulated, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function RunFile(id::Int, filename)
    ccall((:RunFile, Lib_IPhreeqc_path), Cint, (Cint, Cstring), id, filename)
end

function RunString(id::Int, input)
    ccall((:RunString, Lib_IPhreeqc_path), Cint, (Cint, Cstring), id, input)
end

function SetBasicCallback(id::Int, fcn, cookie1)
    ccall((:SetBasicCallback, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Ptr{Nothing}, Ptr{Nothing}), id, fcn, cookie1)
end

function SetBasicFortranCallback(id::Int, fcn)
    ccall((:SetBasicFortranCallback, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Ptr{Nothing}), id, fcn)
end

function SetCurrentSelectedOutputUserNumber(id::Int, n::Int)
    ccall((:SetCurrentSelectedOutputUserNumber, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, n)
end

function SetDumpFileName(id::Int, filename)
    ccall((:SetDumpFileName, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cstring), id, filename)
end

function SetDumpFileOn(id::Int, dump_on::Int)
    ccall((:SetDumpFileOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, dump_on)
end

function SetDumpStringOn(id::Int, dump_string_on::Int)
    ccall((:SetDumpStringOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, dump_string_on)
end

function SetErrorFileName(id::Int, filename)
    ccall((:SetErrorFileName, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cstring), id, filename)
end

function SetErrorFileOn(id::Int, error_on::Int)
    ccall((:SetErrorFileOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, error_on)
end

function SetErrorStringOn(id::Int, error_string_on::Int)
    ccall((:SetErrorStringOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, error_string_on)
end

function SetLogFileName(id::Int, filename)
    ccall((:SetLogFileName, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cstring), id, filename)
end

function SetLogFileOn(id::Int, log_on::Int)
    ccall((:SetLogFileOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, log_on)
end

function SetLogStringOn(id::Int, log_string_on::Int)
    ccall((:SetLogStringOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, log_string_on)
end

function SetOutputFileName(id::Int, filename)
    ccall((:SetOutputFileName, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cstring), id, filename)
end

function SetOutputFileOn(id::Int, output_on::Int)
    ccall((:SetOutputFileOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, output_on)
end

function SetOutputStringOn(id::Int, output_string_on::Int)
    ccall((:SetOutputStringOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, output_string_on)
end

function SetSelectedOutputFileName(id::Int, filename)
    ccall((:SetSelectedOutputFileName, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cstring), id, filename)
end

function SetSelectedOutputFileOn(id::Int, sel_on::Int)
    ccall((:SetSelectedOutputFileOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, sel_on)
end

function SetSelectedOutputStringOn(id::Int, sel_string_on::Int)
    ccall((:SetSelectedOutputStringOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, sel_string_on)
end

function VarAllocString(pSrc)
    ccall((:VarAllocString, Lib_IPhreeqc_path), Cstring, (Cstring,), pSrc)
end

function VarClear(pvar)
    ccall((:VarClear, Lib_IPhreeqc_path), VRESULT, (Ptr{VAR},), pvar)
end

function VarCopy(pvarDest, pvarSrc)
    ccall((:VarCopy, Lib_IPhreeqc_path), VRESULT, (Ptr{VAR}, Ptr{VAR}), pvarDest, pvarSrc)
end

function VarFreeString(pSrc)
    ccall((:VarFreeString, Lib_IPhreeqc_path), Nothing, (Cstring,), pSrc)
end

function VarInit(pvar)
    ccall((:VarInit, Lib_IPhreeqc_path), Nothing, (Ptr{VAR},), pvar)
end

function AccumulateLine(id::Int, line)
    ccall((:AccumulateLine, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cstring), id, line)
end

function AddError(id::Int, error_msg)
    ccall((:AddError, Lib_IPhreeqc_path), Cint, (Cint, Cstring), id, error_msg)
end

function AddWarning(id::Int, warn_msg)
    ccall((:AddWarning, Lib_IPhreeqc_path), Cint, (Cint, Cstring), id, warn_msg)
end

function ClearAccumulatedLines(id::Int)
    ccall((:ClearAccumulatedLines, Lib_IPhreeqc_path), IPQ_RESULT, (Cint,), id)
end

function CreateIPhreeqc()
    ccall((:CreateIPhreeqc, Lib_IPhreeqc_path), Cint, ())
end

function DestroyIPhreeqc(id::Int)
    ccall((:DestroyIPhreeqc, Lib_IPhreeqc_path), IPQ_RESULT, (Cint,), id)
end

function GetComponent(id::Int, n::Int)
    ccall((:GetComponent, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
end

function GetComponentCount(id::Int)
    ccall((:GetComponentCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetCurrentSelectedOutputUserNumber(id::Int)
    ccall((:GetCurrentSelectedOutputUserNumber, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetDumpFileName(id::Int)
    ccall((:GetDumpFileName, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetDumpFileOn(id::Int)
    ccall((:GetDumpFileOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetDumpString(id::Int)
    ccall((:GetDumpString, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetDumpStringLine(id::Int, n::Int)
    ccall((:GetDumpStringLine, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
end

function GetDumpStringLineCount(id::Int)
    ccall((:GetDumpStringLineCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetDumpStringOn(id::Int)
    ccall((:GetDumpStringOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetErrorFileName(id::Int)
    ccall((:GetErrorFileName, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetErrorFileOn(id::Int)
    ccall((:GetErrorFileOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetErrorString(id::Int)
    ccall((:GetErrorString, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetErrorStringLine(id::Int, n::Int)
    ccall((:GetErrorStringLine, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
end

function GetErrorStringLineCount(id::Int)
    ccall((:GetErrorStringLineCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetErrorStringOn(id::Int)
    ccall((:GetErrorStringOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetLogFileName(id::Int)
    ccall((:GetLogFileName, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetLogFileOn(id::Int)
    ccall((:GetLogFileOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetLogString(id::Int)
    ccall((:GetLogString, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetLogStringLine(id::Int, n::Int)
    ccall((:GetLogStringLine, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
end

function GetLogStringLineCount(id::Int)
    ccall((:GetLogStringLineCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetLogStringOn(id::Int)
    ccall((:GetLogStringOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetNthSelectedOutputUserNumber(id::Int, n::Int)
    ccall((:GetNthSelectedOutputUserNumber, Lib_IPhreeqc_path), Cint, (Cint, Cint), id, n)
end

function GetOutputFileName(id::Int)
    ccall((:GetOutputFileName, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetOutputFileOn(id::Int)
    ccall((:GetOutputFileOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetOutputString(id::Int)
    ccall((:GetOutputString, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetOutputStringLine(id::Int, n::Int)
    ccall((:GetOutputStringLine, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
end

function GetOutputStringLineCount(id::Int)
    ccall((:GetOutputStringLineCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetOutputStringOn(id::Int)
    ccall((:GetOutputStringOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputColumnCount(id::Int)
    ccall((:GetSelectedOutputColumnCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputCount(id::Int)
    ccall((:GetSelectedOutputCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputFileName(id::Int)
    ccall((:GetSelectedOutputFileName, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetSelectedOutputFileOn(id::Int)
    ccall((:GetSelectedOutputFileOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputRowCount(id::Int)
    ccall((:GetSelectedOutputRowCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputString(id::Int)
    ccall((:GetSelectedOutputString, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetSelectedOutputStringLine(id::Int, n::Int)
    ccall((:GetSelectedOutputStringLine, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
end

function GetSelectedOutputStringLineCount(id::Int)
    ccall((:GetSelectedOutputStringLineCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputStringOn(id::Int)
    ccall((:GetSelectedOutputStringOn, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function GetSelectedOutputValue(id::Int, row::Int, col::Int, pVAR)
    ccall((:GetSelectedOutputValue, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint, Cint, Ptr{VAR}), id, row, col, pVAR)
end

function GetSelectedOutputValue2(id::Int, row::Int, col::Int, vtype, dvalue, svalue, svalue_length::UInt32)
    ccall((:GetSelectedOutputValue2, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint, Cint, Ptr{Cint}, Ptr{Cdouble}, Cstring, UInt32), id, row, col, vtype, dvalue, svalue, svalue_length)
end

function GetVersionString()
    ccall((:GetVersionString, Lib_IPhreeqc_path), Cstring, ())
end

function GetWarningString(id::Int)
    ccall((:GetWarningString, Lib_IPhreeqc_path), Cstring, (Cint,), id)
end

function GetWarningStringLine(id::Int, n::Int)
    ccall((:GetWarningStringLine, Lib_IPhreeqc_path), Cstring, (Cint, Cint), id, n)
end

function GetWarningStringLineCount(id::Int)
    ccall((:GetWarningStringLineCount, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function LoadDatabase(id::Int, filename)
    ccall((:LoadDatabase, Lib_IPhreeqc_path), Cint, (Cint, Cstring), id, filename)
end

function LoadDatabaseString(id::Int, input)
    ccall((:LoadDatabaseString, Lib_IPhreeqc_path), Cint, (Cint, Cstring), id, input)
end

function OutputAccumulatedLines(id::Int)
    ccall((:OutputAccumulatedLines, Lib_IPhreeqc_path), Nothing, (Cint,), id)
end

function OutputErrorString(id::Int)
    ccall((:OutputErrorString, Lib_IPhreeqc_path), Nothing, (Cint,), id)
end

function OutputWarningString(id::Int)
    ccall((:OutputWarningString, Lib_IPhreeqc_path), Nothing, (Cint,), id)
end

function RunAccumulated(id::Int)
    ccall((:RunAccumulated, Lib_IPhreeqc_path), Cint, (Cint,), id)
end

function RunFile(id::Int, filename)
    ccall((:RunFile, Lib_IPhreeqc_path), Cint, (Cint, Cstring), id, filename)
end

function RunString(id::Int, input)
    ccall((:RunString, Lib_IPhreeqc_path), Cint, (Cint, Cstring), id, input)
end

function SetBasicCallback(id::Int, fcn, cookie1)
    ccall((:SetBasicCallback, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Ptr{Nothing}, Ptr{Nothing}), id, fcn, cookie1)
end

function SetBasicFortranCallback(id::Int, fcn)
    ccall((:SetBasicFortranCallback, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Ptr{Nothing}), id, fcn)
end

function SetCurrentSelectedOutputUserNumber(id::Int, n::Int)
    ccall((:SetCurrentSelectedOutputUserNumber, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, n)
end

function SetDumpFileName(id::Int, filename)
    ccall((:SetDumpFileName, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cstring), id, filename)
end

function SetDumpFileOn(id::Int, dump_on::Int)
    ccall((:SetDumpFileOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, dump_on)
end

function SetDumpStringOn(id::Int, dump_string_on::Int)
    ccall((:SetDumpStringOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, dump_string_on)
end

function SetErrorFileName(id::Int, filename)
    ccall((:SetErrorFileName, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cstring), id, filename)
end

function SetErrorFileOn(id::Int, error_on::Int)
    ccall((:SetErrorFileOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, error_on)
end

function SetErrorStringOn(id::Int, error_string_on::Int)
    ccall((:SetErrorStringOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, error_string_on)
end

function SetLogFileName(id::Int, filename)
    ccall((:SetLogFileName, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cstring), id, filename)
end

function SetLogFileOn(id::Int, log_on::Int)
    ccall((:SetLogFileOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, log_on)
end

function SetLogStringOn(id::Int, log_string_on::Int)
    ccall((:SetLogStringOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, log_string_on)
end

function SetOutputFileName(id::Int, filename)
    ccall((:SetOutputFileName, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cstring), id, filename)
end

function SetOutputFileOn(id::Int, output_on::Int)
    ccall((:SetOutputFileOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, output_on)
end

function SetOutputStringOn(id::Int, output_string_on::Int)
    ccall((:SetOutputStringOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, output_string_on)
end

function SetSelectedOutputFileName(id::Int, filename)
    ccall((:SetSelectedOutputFileName, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cstring), id, filename)
end

function SetSelectedOutputFileOn(id::Int, sel_on::Int)
    ccall((:SetSelectedOutputFileOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, sel_on)
end

function SetSelectedOutputStringOn(id::Int, sel_string_on::Int)
    ccall((:SetSelectedOutputStringOn, Lib_IPhreeqc_path), IPQ_RESULT, (Cint, Cint), id, sel_string_on)
end
# Julia wrapper for header: Var.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function VarAllocString(pSrc)
    ccall((:VarAllocString, Lib_IPhreeqc_path), Cstring, (Cstring,), pSrc)
end

function VarClear(pvar)
    ccall((:VarClear, Lib_IPhreeqc_path), VRESULT, (Ptr{VAR},), pvar)
end

function VarCopy(pvarDest, pvarSrc)
    ccall((:VarCopy, Lib_IPhreeqc_path), VRESULT, (Ptr{VAR}, Ptr{VAR}), pvarDest, pvarSrc)
end

function VarFreeString(pSrc)
    ccall((:VarFreeString, Lib_IPhreeqc_path), Nothing, (Cstring,), pSrc)
end

function VarInit(pvar)
    ccall((:VarInit, Lib_IPhreeqc_path), Nothing, (Ptr{VAR},), pvar)
end

end # module


text:
__precompile__(false)
module Dolark

    import YAML
    import Dolang
    import Dolo
    import OrderedCollections: OrderedDict
    using StaticArrays

    # extracts annotated yaml structure from model file
    # the model file is supposed to contain two documents
    # - a valid dolo model (a.k.a. agent)
    # - a description of aggregate model
    function get_hmodel_source(txt::AbstractString)
        cons = YAML.Constructor()
        YAML.add_multi_constructor!((c,s,m)->m, cons, "tag:yaml.org")
        YAML.add_multi_constructor!((c,s,m)->m, cons, "!")
        data = YAML.load_all(txt, cons)
        agent, model = data
        return (;agent, model)
    end

    struct HModel{ID, AgentType}
        filename
        source
        symbols
        calibration
        exogenous
        equations
        factories
        agent::AgentType # dolo model
    end

    import Base: show
    Base.show(io::IO, model::HModel) = print(io, "HModel")

    function get_symbols(yaml_node)
        symbols_src = yaml_node[:symbols]
        exogenous = [Symbol(e.value) for e in symbols_src[:exogenous]]
        aggregate = [Symbol(e.value) for e in symbols_src[:aggregate]]
        parameters = [Symbol(e.value) for e in symbols_src[:parameters]]
        return (;exogenous, aggregate, parameters)
    end

    function get_h_equations(yaml_node)
        projection = yaml_node[:projection].value
        equilibrium = yaml_node[:equilibrium].value
        return (;projection, equilibrium)
    end

    ### the next two functions are almost duplicates

    function gen_projection_ff(hmodel)
        projs_block = Dolang.parse_assignment_block(hmodel.equations.projection)

        eqs = OrderedDict{Symbol, Any}()
        for ch in projs_block.children
            lhs, rhs = ch.children
            name = Dolang.convert(Expr, lhs)
            eqs[name] = Dolang.convert(Expr, rhs)
        end

        arguments = OrderedDict{Symbol, Vector{Symbol}}()
        arguments[:y] = hmodel.symbols.aggregate
        arguments[:z] = hmodel.symbols.exogenous
        arguments[:p] = hmodel.symbols.parameters

        preamble = OrderedDict()

        return Dolang.FunctionFactory(eqs, arguments, preamble, :projection)
    end

    function gen_projection_ff(agent, symbols, equations)
        projs_block = Dolang.parse_assignment_block(equations.projection)

        eqs = OrderedDict{Symbol, Any}()
        for ch in projs_block.children
            lhs, rhs = ch.children
            name = Dolang.convert(Expr, lhs)
            eqs[name] = Dolang.convert(Expr, rhs)
        end

        arguments = OrderedDict{Symbol, Vector{Symbol}}()
        arguments[:y] = [Dolang.stringify(e,0) for e in symbols.aggregate]
        arguments[:z] = [Dolang.stringify(e,0) for e in symbols.exogenous]
        arguments[:p] = [Dolang.stringify(e) for e in symbols.parameters]

        preamble = OrderedDict()

        return Dolang.FunctionFactory(eqs, arguments, preamble, :projection)
    end

    import Dolang: Tree

    function gen_equilibrium_ff(agent, symbols, equations)

        eq_block = Dolang.parse_equation_block(equations.equilibrium)

        eqs = OrderedDict{Symbol, Any}()
        for (i,ch) in enumerate(eq_block.children)
            eq = Tree("sub", [ch.children[2], ch.children[1]])
            name = Symbol("out_", i, "_")
            eqs[name] = Dolang.convert(Expr, eq)
        end

        arguments = OrderedDict{Symbol, Vector{Symbol}}()
        arguments[:s] = [Dolang.stringify(e,0) for e in agent.symbols[:states]]
        arguments[:x] = [Dolang.stringify(e,0) for e in agent.symbols[:controls]]
        arguments[:y] = [Dolang.stringify(e,0) for e in symbols.aggregate]
        arguments[:z] = [Dolang.stringify(e,0) for e in symbols.exogenous]
        arguments[:Y] = [Dolang.stringify(e,1) for e in symbols.aggregate]
        arguments[:Z] = [Dolang.stringify(e,1) for e in symbols.exogenous]
        arguments[:p] = [Dolang.stringify(e) for e in symbols.parameters]
        preamble = OrderedDict()

        return Dolang.FunctionFactory(eqs, arguments, preamble, :equilibrium)

    end

    module ModelSpace

    end

    # import Dolark.ModelSpace

    function HModel(filename::AbstractString)

        txt = open(f->read(f, String), filename)

        source = get_hmodel_source(txt)
        symbols = get_symbols(source.model)
        calibration = Dolo.get_calibration(source.model)
        exogenous = Dolo.get_exogenous(source.model, symbols.exogenous, calibration.flat)
        equations = get_h_equations(source.model)

        agent = Dolo.Model(source.agent)
        AgentType = typeof(agent)
        ID = gensym()

        projection = gen_projection_ff(agent, symbols, equations)
        equilibrium = gen_equilibrium_ff(agent, symbols, equations)

        factories = (;projection, equilibrium)

        hmodel = HModel{ID, AgentType}(filename, source, symbols, calibration, exogenous, equations, factories, agent)

        tm = typeof(hmodel)

        funp = Dolang.gen_generated_gufun(factories.projection; funname=:projection, dispatch=tm)
        fune = Dolang.gen_generated_gufun(factories.equilibrium; funname=:equilibrium, dispatch=tm)

        Core.eval(Dolark, funp)
        Core.eval(Dolark, fune)

        return hmodel

    end

    struct DModel
        hmodel
        F
        G
    end

    function discretize(hmodel::HModel, sol_agent)

        F = Dolo.Euler(hmodel.agent)
        G = Dolo.distG(hmodel.agent, sol_agent)

        DModel(hmodel, F, G)

    end

    struct Unknown
        μ
        p
        x::Dolo.MSM
        y
    end

    using Dolang: _get_oorders

    include("temp.jl")

end # module