text:
# some basic linear algebra stuff missing from stdlib

using LinearAlgebra.LAPACK: BlasInt, chklapackerror, @blasfunc, liblapack
using LinearAlgebra.LAPACK: checksquare

"""
    condInfest(A,F,anorm)

computes an approximation to the condition of matrix `A` in
the infinity-norm, using factorization `F` and the precomputed
infinity norm `anorm` of `A`.
"""
function condInfest(A::StridedMatrix{T},F::Factorization{T},
                    anorm=opnorm(A,Inf)) where {T}
    γ = normInfest(F) * anorm
end

"""
    norm1est!(applyA!,applyAH!,y::Vector) => γ

Estimate the 1-norm of a linear operator `A` expressed as functions which
apply `A` and `adjoint(A)` to a vector such as `y`.

cf. N.J. Higham, SIAM J. Sci. Stat. Comp. 11, 804 (1990)
"""
function norm1est!(applyA!,applyAH!,x::AbstractVector{T}) where {T}
    n = length(x)
    RT = real(T)
    x = fill(one(T)/n,n)
    y = copy(x)
    z = similar(y)
    za = Vector{RT}(undef,n)
    asign(a::Real) = a >= zero(T) ? one(T) : -one(T)
    asign(a::Complex) = a == zero(T) ? one(T) : a / abs(a)
    γ = zero(RT)
    jprev=0
    for iter=1:5
        applyA!(y)
        z = asign.(y)
        applyAH!(z)
        za .= abs2.(z)
        zam = maximum(za)
        j = findfirst(za .== zam)
        if (iter > 1) && (zam <= za[jprev])
            γ = norm(y,1)
            break
        end
        fill!(x,zero(T))
        x[j] = one(T)
        jprev = j
    end
    v,w = x,z
    v = T.((n-1:2n-2)/(n-1))
    for j=2:2:n
        v[j] = -v[j]
    end
    vnorm = norm(v,1)
    applyA!(v)
    max(γ, norm(v,1) / vnorm)
end

function norm1est(F::Factorization{T}) where {T}
    n = size(F,1)
    y = Vector{T}(undef, n)
    norm1est!(x->ldiv!(F,x),x->ldiv!(F',x),y)
end

function normInfest(F::Factorization{T}) where {T}
    n = size(F,1)
    y = Vector{T}(undef, n)
    norm1est!(x->ldiv!(F',x), x->ldiv!(F,x), y)
end

"""
    equilibrators(A) -> R,C

compute row- and column-wise scaling vectors `R,C` for a matrix `A`
such that the absolute value of the largest element in any row or
column of `Diagonal(R)*A*Diagonal(C)` is close to unity. Designed to
reduce the condition number of the working matrix.
"""
function equilibrators(A::AbstractMatrix{T}) where {T}
    abs1(x::Real) = abs(x)
    abs1(x::Complex) = abs(real(x)) + abs(imag(x))
    m,n = size(A)
    R = zeros(T,m)
    C = zeros(T,n)
    @inbounds for j=1:n
        for i=1:m
            R[i] = max(R[i],abs(A[i,j]))
        end
    end
    @inbounds for i=1:m
        if R[i] > 0
            R[i] = T(2)^floor(Int,log2(R[i]))
        end
    end
    R .= 1 ./ R
    @inbounds for i=1:m
        for j=1:n
            C[j] = max(C[j],R[i] * abs(A[i,j]))
        end
    end
    @inbounds for j=1:n
        if C[j] > 0
            C[j] = T(2)^floor(Int,log2(C[j]))
        end
    end
    C .= 1 ./ C
    R,C
end

const BlasTypes = Union{Float32,Float64,ComplexF32,ComplexF64}
# can use LAPACK.gecon for BLAS types
function condInfest(A::StridedMatrix{T},F::Factorization{T},
                    anorm=opnorm(A,Inf)) where {T<:BlasTypes}
    1/LAPACK.gecon!('I',F.factors,anorm)
end

# can use LAPACK.geequb for BLAS types
function equilibrators(A::AbstractMatrix{T}) where {T<:BlasTypes}
    Rv, Cv, rowcond, colcond, amax = geequb(A)
    Rv,Cv
end

# but first we need to wrap it...
for (geequb, elty, relty) in
    ((:dgeequb_, :Float64, :Float64),
     (:zgeequb_, :ComplexF64, :Float64),
     (:cgeequb_, :ComplexF32, :Float32),
     (:sgeequb_, :Float32, :Float32))
    @eval begin
#=
*       SUBROUTINE DGEEQUB( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
*                           INFO )
*
*       .. Scalar Arguments ..
*       INTEGER            INFO, LDA, M, N
*       DOUBLE PRECISION   AMAX, COLCND, ROWCND
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   A( LDA, * ), C( * ), R( * )
=#
        function geequb(A::AbstractMatrix{$elty})
            m,n = size(A)
            lda = max(1, stride(A,2))
            C = Vector{$relty}(undef, n)
            R = Vector{$relty}(undef, m)
            info = Ref{BlasInt}()
            rowcond = Ref{$relty}()
            colcond = Ref{$relty}()
            amax = Ref{$relty}()
            ccall((@blasfunc($geequb), liblapack), Cvoid,
                  (Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                   Ptr{$relty}, Ptr{$relty},
                   Ptr{$relty}, Ptr{$relty}, Ptr{$relty},
                   Ptr{BlasInt}),
                  m, n, A, lda, R, C, rowcond, colcond, amax, info)
            chklapackerror(info[])
            R, C, rowcond, colcond, amax
        end
    end
end


text:
export bearing, final_bearing

"""
    bearing(pos₁::Point, pos₂::Point)

Return the initial `bearing` [deg] of the great circle line between the
positions `pos₁` and `pos₂` [deg] on a sphere.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function bearing(pos₁::Point, pos₂::Point)
    Δpos = pos₂ - pos₁
    return mod(atand(sind(Δpos.λ)*cosd(pos₂.ϕ), cosd(pos₁.ϕ)*sind(pos₂.ϕ) -
    sind(pos₁.ϕ)*cosd(pos₂.ϕ)*cosd(Δpos.λ)), 360.0)
end

"""
    bearing(section::RouteSection)

Return the initial `bearing` [deg] of the route section on a sphere.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
bearing(section::RouteSection) = bearing(section.pos₁, section.pos₂)

"""
    final_bearing(pos₁::Point, pos₂::Point)

Return the `final_bearing` [deg] of the great circle line between the positions
`pos₁` and `pos₂` [deg] on a sphere.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function final_bearing(pos₁::Point, pos₂::Point)
    return mod(bearing(pos₂, pos₁) + 180.0, 360.0)
end

"""
    final_bearing(section::RouteSection)

Return the `final_bearing` [deg] of the route section on a sphere.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
final_bearing(section::RouteSection) = final_bearing(section.pos₁, section.pos₂)

ext: jl


text:
#####
##### replace the storage place of the grid information based on the architecture
#####
function replace_grid_storage(arch::Architecture, grid::RegularLatLonGrid{TX, TY, TZ}) where {TX, TY, TZ}
    dxC = grid.dxC |> array_type(arch)
    dyC = grid.dyC |> array_type(arch)
    dxF = grid.dxF |> array_type(arch)
    dyF = grid.dyF |> array_type(arch)
    Ax  = grid.Ax  |> array_type(arch)
    Ay  = grid.Ay  |> array_type(arch)
    Az  = grid.Az  |> array_type(arch)
    Vol = grid.Vol |> array_type(arch)

    return RegularLatLonGrid{TX, TY, TZ, typeof(grid.xF), typeof(dxC)}(
        grid.xC, grid.yC, grid.zC, grid.xF, grid.yF, grid.zF, grid.Δx, grid.Δy, grid.Δz,
        dxC, dyC, dxF, dyF, Ax, Ay, Az, Vol, grid.Nx, grid.Ny, grid.Nz, grid.Hx, grid.Hy, grid.Hz)
end

function replace_grid_storage(arch::Architecture, grid::VerticallyStretchedLatLonGrid{TX, TY, TZ}) where {TX, TY, TZ}
    zF  = grid.zF |> array_type(arch)
    zC  = grid.zC |> array_type(arch)
    dxC = grid.dxC |> array_type(arch)
    dyC = grid.dyC |> array_type(arch)
    dzC = grid.dzC |> array_type(arch)
    dxF = grid.dxF |> array_type(arch)
    dyF = grid.dyF |> array_type(arch)
    dzF = grid.dzF |> array_type(arch)
    Ax  = grid.Ax  |> array_type(arch)
    Ay  = grid.Ay  |> array_type(arch)
    Az  = grid.Az  |> array_type(arch)
    Vol = grid.Vol |> array_type(arch)

    return VerticallyStretchedLatLonGrid{TX, TY, TZ, typeof(grid.xF), typeof(zF), typeof(dxC), typeof(Vol)}(
        grid.xC, grid.yC, zC, grid.xF, grid.yF, zF, grid.Δx, grid.Δy, dxC, dyC, dzC,
        dxF, dyF, dzF, Ax, Ay, Az, Vol, grid.Nx, grid.Ny, grid.Nz, grid.Hx, grid.Hy, grid.Hz)
end

function replace_grid_storage(arch::Architecture, grid::RegularRectilinearGrid{TX, TY, TZ}) where {TX, TY, TZ}
    return grid
end

#####
##### adapt the grid struct to GPU
#####

Adapt.adapt_structure(to, grid::RegularLatLonGrid{TX, TY, TZ}) where {TX, TY, TZ} =
    RegularLatLonGrid{TX, TY, TZ, typeof(grid.xF), typeof(Adapt.adapt(to, grid.dxC))}(
        grid.xC, grid.yC, grid.zC,
        grid.xF, grid.yF, grid.zF,
        grid.Δx, grid.Δy, grid.Δz,
        Adapt.adapt(to, grid.dxC),
        Adapt.adapt(to, grid.dyC),
        Adapt.adapt(to, grid.dxF),
        Adapt.adapt(to, grid.dyF),
        Adapt.adapt(to, grid.Ax),
        Adapt.adapt(to, grid.Ay),
        Adapt.adapt(to, grid.Az),
        Adapt.adapt(to, grid.Vol),
        grid.Nx, grid.Ny, grid.Nz,
        grid.Hx, grid.Hy, grid.Hz)

Adapt.adapt_structure(to, grid::VerticallyStretchedLatLonGrid{TX, TY, TZ}) where {TX, TY, TZ} =
    VerticallyStretchedLatLonGrid{TX, TY, TZ, typeof(grid.xF), typeof(Adapt.adapt(to, grid.zF)),
                                  typeof(Adapt.adapt(to, grid.dxC)), typeof(Adapt.adapt(to, grid.Vol))}(
        grid.xC, grid.yC,
        Adapt.adapt(to, grid.zC),
        grid.xF, grid.yF,
        Adapt.adapt(to, grid.zF),
        grid.Δx, grid.Δy,
        Adapt.adapt(to, grid.dxC),
        Adapt.adapt(to, grid.dyC),
        Adapt.adapt(to, grid.dzC),
        Adapt.adapt(to, grid.dxF),
        Adapt.adapt(to, grid.dyF),
        Adapt.adapt(to, grid.dzF),
        Adapt.adapt(to, grid.Ax),
        Adapt.adapt(to, grid.Ay),
        Adapt.adapt(to, grid.Az),
        Adapt.adapt(to, grid.Vol),
        grid.Nx, grid.Ny, grid.Nz,
        grid.Hx, grid.Hy, grid.Hz)


text:
using Test
using ZeroRing

zeroelem::ZeroElem

@test zeroelem === zeroelem

@test zeroelem == zeroelem
@test !(zeroelem != zeroelem)
@test !(zeroelem < zeroelem)
@test hash(zeroelem) === hash(zeroelem)

@test rand(ZeroElem) === zeroelem

@test +zeroelem === zeroelem
@test -zeroelem === zeroelem

@test zeroelem + zeroelem === zeroelem
@test zeroelem - zeroelem === zeroelem

@test inv(zeroelem) === zeroelem

@test zeroelem * zeroelem === zeroelem
@test zeroelem / zeroelem === zeroelem

@test zeroelem^zeroelem === zeroelem

@test copy(zeroelem) === zeroelem

@test zeroelem' === zeroelem
@test transpose(zeroelem) === zeroelem


text:
# -------------- #
# B-Spline Basis #
# -------------- #

struct Spline <: BasisFamily end

mutable struct SplineParams{T<:AbstractVector} <: BasisParams
    breaks::T
    evennum::Int
    k::Int

    # constructor to accept spline params arguments, do some pre-processing
    function SplineParams{T}(breaks::T, evennum::Int, k::Int) where T
        # error handling
        k < 0 && error("spline order must be positive")
        length(breaks) < 2 && error("Must have at least two breakpoints")
        any(diff(breaks) .< 0) && error("Breakpoints must be non-decreasing")

        if evennum == 0  # 43
            if length(breaks) == 2  # 44
                evennum = 2
            end
        else
            if length(breaks) == 2
                breaks = linspace(breaks[1], breaks[2], evennum)
            else
                error("Breakpoint squence must contain 2 values when evennum > 0")
            end
        end
        new{T}(breaks, evennum, k)
    end

end

SplineParams(breaks::T, evennum::Int, k::Int) where {T<:AbstractVector} =
    SplineParams{T}(breaks, evennum, k)

# constructor to take a, b, n and form linspace for breaks
SplineParams(n::Int, a::Real, b::Real, k::Int=3) =
    SplineParams(linspace(a, b, n), 0, k)

## BasisParams interface
# define these methods on the type, the instance version is defined over
# BasisParams
family(::Type{T}) where {T<:SplineParams} = Spline
family_name(::Type{T}) where {T<:SplineParams} = "Spline"
Base.issparse(::Type{T}) where {T<:SplineParams} = true
function Base.eltype(::Type{SplineParams{T}}) where T
    elT = eltype(T)
    elT <: Integer ? Float64 : elT
end

# methods that only make sense for instances
Base.min(p::SplineParams) = minimum(p.breaks)
Base.max(p::SplineParams) = maximum(p.breaks)
Base.length(p::SplineParams) = length(p.breaks) + p.k - 1

function Base.show(io::IO, p::SplineParams)
    m = string("$(p.k) order spline interpoland parameters from ",
               "$(p.breaks[1]), $(p.breaks[end])")
    print(io, m)
end

"""
Construct interpolation nodes, given SplineParams

Note that `p.k - 1` additional nodes will be inserted

##### Arguments

- `p::SplineParams`: `SplineParams` instance

##### Returns

- `x::Vector`: The Vector of 1d interpolation nodes.

"""
function nodes(p::SplineParams)
    breaks, evennum, k = p.breaks, p.evennum, p.k
    a = breaks[1]  # 20
    b = breaks[end]  # 21
    n = length(breaks) + k - 1  # 22
    x = cumsum(vcat(fill(a, k), breaks, fill(b, k)))  # 23
    x = (x[1+k:n+k] - x[1:n]) / k  # 24
    x[1] = a  # 25
    x[end] = b  # 26
    x
end

# TODO: define method derivative_op(::Type{SplineSparse}, p::SplineParams, order::Int)
function derivative_op(p::SplineParams, x::AbstractArray, order=1)
    breaks, evennum, k = p.breaks, p.evennum, p.k

    any(order .> k) && error("Order of differentiation can't be greater than k")

    # 38-40
    n = length(breaks) + k - 1
    kk = max(k - 1, k - order - 1)
    augbreaks = vcat(fill(breaks[1], kk), breaks, fill(breaks[end], kk))

    D = Array{SparseMatrixCSC{basis_eltype(p, x),Int64}}(abs(order), 1)

    if order > 0  # derivative
        temp = k ./ (augbreaks[k+1:n+k-1] - augbreaks[1:n-1])
        D[1] = spdiagm((-temp, temp), 0:1, n-1, n)

        for i in 2:order
            temp = (k+1-i) ./ (augbreaks[k+1:n+k-i] - augbreaks[i:n-1])
            D[i] = spdiagm((-temp, temp), 0:1, n-i, n+1-i)*D[i-1]
        end
    else
        error("not implemented")
    end

    D, SplineParams(breaks, evennum, k-order)
end

function _chk_evalbase(p::SplineParams, x::AbstractArray, order)
    breaks, evennum, k = p.breaks, p.evennum, p.k

    # error handling
    k < 0 && error("spline order must be positive")
    !(issorted(breaks)) && error("Breakpoints must be non-decreasing")
    any(order .>= k) && error("Order of differentiation must be less than k")
    size(x, 2) > 1 && error("x must be a column vector")

    m = size(x, 1)  # 54
    minorder = minimum(order)  # 55

    # Augment the breakpoint sequence 57-59
    n = length(breaks)+k-1
    a = breaks[1]
    b = breaks[end]
    augbreaks = vcat(fill(a, k-minorder), breaks, fill(b, k-minorder))

    ind = lookup(augbreaks, x, 3)  # 69

    n, m, minorder, augbreaks, ind
end

function evalbase(p::SplineParams, x::AbstractArray=nodes(p), order::Int=0)
    evalbase(p, x, [order])[1]
end

function evalbase(
        ::Type{SplineSparse}, p::SplineParams,
        x::AbstractArray=nodes(p), order::Int=0
    )
    evalbase(SplineSparse, p, x, [order])[1]
end

"""
Evaluate spline basis matrices for a certain order derivative at x

##### Arguments

- `p::SplineParams`: A `SplineParams` summarizing spline properties
- `x(nodes(p))` : the nodes at which to evaluate the basis matrices
- `order(0)` : The order(s) of derivative for which to evaluate the basis
matrices. `order=0` corresponds to the function itself, negative numbers
correspond to integrals.

##### Returns

- `B::SparseMatrixCSC` : Matrix containing the evaluation of basis functions
at each point in `x`. Each column represents a basis function.
- `x`: Points at which the functions were evaluated

"""
function evalbase(p::SplineParams, x::AbstractArray, order::AbstractVector{Int})
    n, m, minorder, augbreaks, ind = _chk_evalbase(p, x, order)

    max_repeat = p.k-minorder + 1
    T = basis_eltype(p, x)
    bas = zeros(T, m, max_repeat)  # 73
    bas[:, 1] = one(T)  # 74
    B = Array{SparseMatrixCSC{T,Int}}(length(order))  # 75

    # 76
    if maximum(order) > 0
        D = derivative_op(p, x, maximum(order))[1]
    end

    # 77
    if minorder < 0
        I = derivative_op(p, x, minorder)[1]
    end

    # We know what the rows and columns will be, we just compute them once and
    # then extract chunks as we need them.
    # note: I benchmarked and the version with collect was (2-5)x faster
    r = repmat(collect(1:m), max_repeat)
    c = (minorder-p.k:0)' .+ ind

    for j in 1:p.k-minorder  # 78
        for jj in j:-1:1  # 79
            for ix in eachindex(ind)
                b0 = augbreaks[ind[ix]+jj-j]  # 80
                b1 = augbreaks[ind[ix]+jj]  # 81
                temp = bas[ix, jj] / (b1 - b0)  # 82
                bas[ix, jj+1] = (x[ix] - b0) * temp + bas[ix, jj+1]  # 83
                bas[ix, jj] = (b1-x[ix]) * temp  # 84
            end
        end

        # bas now contains the order `j` spline basis
        ii = findfirst(order, p.k-j)  # 87
        if ii > 0
            # Put values in appropriate columns of a sparse matrix
            N = m * (p.k-order[ii]+1)
            B[ii] = sparse(view(r, 1:N), view(c, 1:N), view(bas, 1:N), m, n-order[ii])

            # 96-100
            if order[ii] > 0
                B[ii] = B[ii] * D[order[ii]]
            elseif order[ii] < 0
                B[ii] = B[ii] * I[-order[ii]]
            end
        end
    end

    B
end

function evalbase(
        ::Type{SplineSparse}, p::SplineParams, x::AbstractArray,
        order::AbstractVector{Int}
    )
    n, m, minorder, augbreaks, ind = _chk_evalbase(p, x, order)

    max_repeat = p.k-minorder + 1
    T = basis_eltype(p, x)
    bas = zeros(T, max_repeat, m)
    bas[1, :] = one(T)
    B = Array{SplineSparse{T,Int}}(length(order))  # 75

    if maximum(order) > 0
        D = derivative_op(p, x, maximum(order))[1]
    end

    if minorder < 0
        error("not supported yet")
        I = derivative_op(p, x, minorder)[1]
    end

    for j in 1:p.k-minorder  # 78
        for ix in eachindex(ind)
            for jj in j:-1:1  # 79
                b0 = augbreaks[ind[ix]+jj-j]  # 80
                b1 = augbreaks[ind[ix]+jj]  # 81
                temp = bas[jj, ix] / (b1 - b0)  # 82
                bas[jj+1, ix] = (x[ix] - b0) * temp + bas[jj+1, ix]  # 83
                bas[jj, ix] = (b1-x[ix]) * temp  # 84
            end
        end

        # bas now contains the order `j` spline basis
        ii = findfirst(order, p.k-j)  # 87
        if ii > 0
            ord = order[ii]
            cols = (ord - p.k) + ind
            if ord > 0
                D_ord = D[ord]
                new_vals = zeros(T, m * (p.k+1))
                ix = 0
                for row in 1:m  # row, col are rows/columns of the output
                    start_bas_col = cols[row]
                    end_bas_col = start_bas_col + p.k - ord
                    for col in 1:n
                        colptr_start = D_ord.colptr[col]
                        colptr_end = D_ord.colptr[col+1] - 1
                        start_D_row = D_ord.rowval[colptr_start]
                        end_D_row = D_ord.rowval[colptr_end]

                        # stop end_D_row is smaller than start_bas_col
                        end_D_row < start_bas_col && continue

                        # also stop if end_bas_col is smaller than start_D_row
                        end_bas_col < start_D_row && continue

                        # otherwise, we have some overlap
                        ix += 1

                        for D_ptr in colptr_start:colptr_end
                            D_row = D_ord.rowval[D_ptr]
                            if D_row < start_bas_col  # haven't made it to overlap
                                continue
                            elseif D_row > end_bas_col  # gone past overlap
                                break
                            else  # in the overlap
                                bas_ix = D_row - start_bas_col + 1
                                new_vals[ix] += D_ord.nzval[D_ptr] * bas[bas_ix, row]
                            end
                        end
                    end
                end

                B[ii] = SplineSparse{T,Int,1,p.k+1}(n, new_vals, cols)
            elseif ord == 0

                B[ii] = SplineSparse{T,Int,1,p.k+1}(n-ord, vec(bas), cols)
            else
                error("Shouldn't be here. Please file an issue")
            end
        end
    end
    return B
end


text:
import .Tracker

preprocess(name, x::T, data) where {T<:Union{Tracker.TrackedReal,Tracker.TrackedArray}}= push!(data, name => Tracker.data(x))


text:
module QSimulator

using Grid,
      NumericExtensions,
      Iterators

include("base-types.jl")

include("systems.jl")

include("controls.jl")

include("interactions.jl")

include("composite-systems.jl")

#include("dissipation.jl")

include("evolution.jl")

include("read_APS_file.jl")

include("paulis.jl")

end


text:
# Test for gausslegendre().

for n=2:10  # check all special cases
    x, w = gausslegendre(n)
    @test_approx_eq dot( w,(x.^2)) 2/3
end

tol = 1e-14
n = 42
x, w = gausslegendre(n)
@test length(x) == n && length(w) == n
@test_approx_eq_eps x[37]  0.910959724904127 tol
@test_approx_eq_eps w[37] 0.030479240699603 tol

@test_approx_eq dot( w,(x.^2)) 2/3
@test_approx_eq dot( w,exp(x)) exp(1)-exp(-1)

# Test a larger n (using ASY)
n = 251
x, w = gausslegendre(n)
@test all(length(x) == n) && all(length(w) == n)
@test abs(x[37] + 0.896467746955729) < tol
@test abs(w[37] - 0.005535005742012) < tol

@test_approx_eq dot( w,(x.^2)) 2/3
@test_approx_eq dot( w,exp(x)) exp(1)-exp(-1)

x, w = gausslegendre(1013)
@test norm(x[2] - -0.999985167586110, Inf) < tol
@test norm(x[13] - -0.999218995240887, Inf) < tol
@test norm(w[2] - 1.681691163200592e-05, Inf) < tol
@test norm(w[13] - 1.224755309137936e-04, Inf) < tol

@test_approx_eq dot( w,(x.^2)) 2/3
@test_approx_eq dot( w,exp(x)) exp(1)-exp(-1)

x, w = gausslegendre(10013)
@test norm(x[2] - -0.999999848054223, Inf) < tol
@test norm(x[13] - -0.999991998242661, Inf) < tol
@test norm(w[2] - 1.722757320118474e-07, Inf) < tol
@test norm(w[13] - 1.254980540032470e-06, Inf) < tol

@test_approx_eq dot( w,(x.^2)) 2/3
@test_approx_eq dot( w,exp(x)) exp(1)-exp(-1)


text:
module ModParameters

export Param, setParameters

const BCDirichlet,BCRiemann = 1,2 # Constants for different boundary types

using TOML

# Parameters
struct Param
   nD::Int             # Dimension of the system
   LatticeType::Int    # DmQn lattice model
   Δx::Float64         # Discretized space interval
   Δt::Float64         # Discretized time interval
   α ::Float64         # Thermal diffusivity
   ω ::Float64         # Relaxation coefficient
   TotalTime::Float64  # Total simulation time
   nT::Int             # Total timesteps

   Temp::Vector{Float64} # Left boundary temperature

   nI::Int
   nJ::Int
   BCtype::Array{Int,1}         # Type of boundary conditions
   BClines::Array{Float64,2}

   DoPlot::Bool

   function Param(BClines::Array{Float64,2},
      nD::Int=2,LatticeType::Int=3,dx::Float64=1.0,dt::Float64=1.0,
      α::Float64=0.25,TotalTime::Float64=400.0,
      Temp::Vector{Float64}=[1.0, 0.0, 0.0, 0.0],
      BCtype::Vector{Int}=[1,1,1,1], nI::Int=5, nJ::Int=5,
      DoPlot::Bool=false)
      @assert(0 ≤ TotalTime, "Simulation time must be positive!")

      new(nD,LatticeType,dx,dt,α,inv(nD*dt*α/dx^2 + 0.5),TotalTime,
         Int(TotalTime/dt),Temp,nI,nJ,BCtype,BClines,DoPlot)
   end
end

"""
   setParameters()

Read parameters from PARAM.toml and construct the parameter list.
"""
function setParameters()

   paramIn = TOML.parsefile("PARAM.toml")
   nD = paramIn["Parameters"]["nD"]
   LatticeType = paramIn["Parameters"]["LatticeType"]
   dx = paramIn["Parameters"]["dx"]
   dt = paramIn["Parameters"]["dt"]
   α = paramIn["Parameters"]["α"]
   TotalTime = paramIn["Parameters"]["TotalTime"]
   Temp = paramIn["Parameters"]["Temp"]
   DoPlot = paramIn["Plots"]["DoPlot"]

   if LatticeType == "D2Q4"
      LT = 1
   elseif LatticeType == "D2Q5"
      LT = 2
   elseif LatticeType == "D2Q9"
      LT = 3
   end

   nI = paramIn["Mesh"]["nI"]
   nJ = paramIn["Mesh"]["nJ"]
   nBCtype = length(paramIn["Mesh"]["BCtype"])
   if nBCtype != size(Temp,1) # This needs to be modified!
      error("Number of boundary values does not match! Check Temp in PARAM!")
   end
   BCtype = Vector{Int}(undef,nBCtype)
   for i = 1:nBCtype
      if paramIn["Mesh"]["BCtype"][i] == "Dirichlet"
         BCtype[i] = 1
      elseif paramIn["Mesh"]["BCtype"][i] == "Riemann"
         BCtype[i] = 2
      end
   end

   BClines = hcat(paramIn["Mesh"]["BClines"]...)

   param = Param(BClines,nD,LT,dx,dt,α,TotalTime,Temp,BCtype,nI,nJ,DoPlot)

   return param
end

end


text:
# module Octo

db_keywords      = String[]
db_functionnames = String[]

"""
    @sql_keywords(args...)
"""
macro sql_keywords(args...)
    esc(sql_keywords(args))
end
function sql_keywords(s)
    for keyword in s
        push!(db_keywords, String(keyword))
    end
    :(($(s...),) = $(map(SQLKeyword, s)))
end

"""
    @sql_functions(args...)
"""
macro sql_functions(args...)
    esc(sql_functions(args))
end
function sql_functions(s)
    for funcname in s
        push!(db_functionnames, String(funcname))
    end
    :(($(s...),) = $(map(x->SQLFunctionName(x), s)))
end

# module Octo


text:
using InfrastructureSystems
using PowerSystems
using InteractiveUtils
const IS = InfrastructureSystems
const PSY = PowerSystems

IS.strip_module_name

function _check_exception(T, exceptions::Vector)
    for type_exception in exceptions
        if T <: type_exception
            return true
        end
    end
    return false
end

function _write_first_level_markdown(c::String)
    file_name = "model_library/generated_$(c).md"
            open(joinpath("docs/src", file_name), "w") do io
                print(
                    io,
                    """
                    # $(c)

                    ```@autodocs
                    Modules = [PowerSystems]
                    Pages   = ["generated/$(c).jl"]
                    Order = [:type, :function]
                    Public = true
                    ```
                    """,
                )
            end
    return file_name
end

function _write_second_level_markdown(input::DataType, subtypes::Vector{DataType}, exceptions)
    c = IS.strip_module_name(input)
    file_name = "model_library/generated_$(c).md"
            open(joinpath("docs/src", file_name), "w") do io
                print(io, "# $input\n\n")
                for T_ in subtypes
                    _check_exception(T_, exceptions) && continue
                    T = IS.strip_module_name(T_)
                    print(
                        io,
                        """
                        ## $(T)

                        ```@autodocs
                        Modules = [PowerSystems]
                        Pages   = ["/$(T).jl"]
                        Order = [:type, :function]
                        Public = true
                        ```

                        """,
                    )
                end
            end
    return file_name
end

function make_dynamics_library!(model_library;
dyn_categories =[
    PSY.DynamicGeneratorComponent,
    PSY.DynamicInverterComponent,
],
exceptions = [PSY.OuterControl,
              PSY.ActivePowerControl,
              PSY.ReactivePowerControl,],

manual_additions = Dict{String, Any}("DynamicInverterComponent" => Any["OuterControl" => "model_library/outer_control.md"])
)
    for abstract_type in dyn_categories
        @info "Making entries for subtypes of $abstract_type"
        abstract_type_string = IS.strip_module_name(abstract_type)
        addition = Dict{String, Any}()
        internal_index = Any[]
        for c_ in subtypes(abstract_type)
            c_string = IS.strip_module_name(c_)
            _check_exception(c_, exceptions) && continue
            concretes = IS.get_all_concrete_subtypes(c_)
            file_name = _write_second_level_markdown(c_,
            concretes, exceptions)
            push!(internal_index, c_string => file_name)
        end
        push!(model_library, abstract_type_string => internal_index)
        if haskey(manual_additions, abstract_type_string)
            addition = get(manual_additions, abstract_type_string, nothing)
            push!(model_library[abstract_type_string], addition...)
        end
    end
end

function make_model_library(;
    categories = [],
    exceptions = [],
    manual_additions = Dict{String, Any}()
)

    model_library = Dict{String, Any}()

    for abstract_type in categories
        @info "Making entries for subtypes of $abstract_type"
        internal_index = Any[]
        concrete = IS.get_all_concrete_subtypes(abstract_type)
        for c_ in concrete
            _check_exception(c_, exceptions) && continue
            c = IS.strip_module_name(c_)
            file_name = _write_first_level_markdown(c)
            push!(internal_index, c => file_name)
        end
        isempty(internal_index) && continue
        model_library[IS.strip_module_name(abstract_type)] = internal_index
    end

    make_dynamics_library!(model_library)

    for (k, v) in manual_additions
        if haskey(model_library, k)
            push!(model_library[k], v...)
        else
            model_library[k] = v
        end
    end

    return Any[p for p in model_library]
end


text:
@doc doc"""
    TranslationAction(
        M::Manifold,
        Rn::TranslationGroup,
        AD::ActionDirection = LeftAction(),
    )

Space of actions of the [`TranslationGroup`](@ref) $\mathrm{T}(n)$ on a Euclidean-like
manifold `M`.

The left and right actions are equivalent.
"""
struct TranslationAction{TM<:Manifold,TRn<:TranslationGroup,TAD<:ActionDirection} <:
       AbstractGroupAction{TAD}
    M::TM
    Rn::TRn
end

function TranslationAction(
    M::Manifold,
    Rn::TranslationGroup,
    ::TAD = LeftAction(),
) where {TAD<:ActionDirection}
    return TranslationAction{typeof(M),typeof(Rn),TAD}(M, Rn)
end

function show(io::IO, A::TranslationAction)
    print(io, "TranslationAction($(A.M), $(A.Rn), $(direction(A)))")
end

base_group(A::TranslationAction) = A.Rn

g_manifold(A::TranslationAction) = A.M

function switch_direction(A::TranslationAction{TM,TRN,TAD}) where {TM,TRN,TAD}
    return TranslationAction(A.M, A.Rn, switch_direction(TAD()))
end

apply(A::TranslationAction, a, x) = x + a

apply!(A::TranslationAction{M,G}, y, a, x) where {M,G} = (y .= x .+ a)
apply!(A::TranslationAction{M,G}, y, e::Identity{G}, x) where {M,G} = copyto!(y, x)

inverse_apply(A::TranslationAction, a, x) = x - a

inverse_apply!(A::TranslationAction{M,G}, y, a, x) where {M,G} = (y .= x .- a)
inverse_apply!(A::TranslationAction{M,G}, y, e::Identity{G}, x) where {M,G} = copyto!(y, x)

apply_diff(A::TranslationAction, a, x, v) = v

apply_diff!(A::TranslationAction, vout, a, x, v) = copyto!(vout, v)

inverse_apply_diff(A::TranslationAction, a, x, v) = v

inverse_apply_diff!(A::TranslationAction, vout, a, x, v) = copyto!(vout, v)


text:
module CubicEquation

export Solver

type Solver
end

function (solver::Solver)(a::Real, b::Real, c::Real, d::Real)
    if a == zero(typeof(a))
        return Base.error("Can't solve. a=0 in a*x³ + b*x² + c*x + d =0")
    end
    f = (3c/a - (b/a)^2) / 3
    g = (2(b/a)^3 - (9*b*c/a^2) + 27d/a) / 27
    h = g^2 / 4 + f^3 / 27
    if f == zero(typeof(f)) && g == zero(typeof(g)) && h == zero(typeof(h))
        x1 = x2 = x3 = -cbrt(d/a)
        return [x1, x2, x3]
    elseif  h <= zero(typeof(h))
        return three_real_root(a, b, g, h)
    elseif h > zero(typeof(h))
        return one_real_root(a, b, g, h)
    end

    return Base.error("Can't solve when a=$(a) b=$(b) c=$(c) d=$(d)")
end

function three_real_root(a::Real, b::Real, g::Float64, h::Float64)
    i = (g^2 / 4 - h)^(1/2)
    j = cbrt(i)
    K = acos(-(g/2i))
    L = -j
    M = cos(K/3)
    N = sqrt(3) * sin(K/3)
    P = -(b/3a)
    x1 = 2j * cos(K/3) - (b/3a)
    x2 = L * (M + N) + P
    x3 = L * (M - N) + P
    return sort([x1, x2, x3])
end

function one_real_root(a::Real, b::Real, g::Float64, h::Float64)
    R = -g/2 + sqrt(h)
    S = cbrt(R)
    T = -g/2 - sqrt(h)
    U = cbrt(T)
    x1 = (S + U) - (b/3a)
    x2 = -(S + U)/2 - (b/3a) - im*(S-U)*sqrt(3) / 2
    x3 = -(S + U)/2 - (b/3a) + im*(S-U)*sqrt(3) / 2
    return Union{Real,Complex}[x1, x2, x3]
end

function (solver::Solver)(a::Real, b::Real, c::Real)
    if a == zero(typeof(a))
        return Base.error("Can't solve. a=0 in a*x² + b*x + c = 0")
    end

    # Discriminant for quadratic equation
    D = b^2 - 4*a*c

    if D >= zero(typeof(D))
        return [(-b - √(D)) / 2a, (-b + √(D)) / 2a]
    elseif D < zero(typeof(D))
        return [(-b - √(-D)*im) / 2a, (-b + √(-D)*im) / 2a]
    else
        return Base.error("Can't solve when a=$(a) b=$(b) c=$(c)")
    end
end

end #end module


text:
include("exampleTrajectories.jl")
include("exampleRadial.jl")
include("exampleFieldmap.jl")
include("exampleSENSE.jl")
include("exampleIO.jl")


text:
module DA_DCSTests

using Test
using LightPropagation

τ = 10 .^(range(-10,stop=0,length=250))

# test that the sum of all pathlengths (solution given in time-domain) matches CW solution
@test isapprox(g1_DA_semiinf_TD(τ,[0.00001, 10.0], 1.0, 0.1, 10.0, N_quad = 200), g1_DA_semiinf_CW(τ, 1.0, 0.1, 10.0))

# check other optical properties and index of refraction
@test g1_DA_semiinf_TD(τ,[0.00001, 10.0], 0.5, 0.04, 14.2, n_med = 1.2, n_ext = 1.1, N_quad = 300) ≈ g1_DA_semiinf_CW(τ, 0.5, 0.04, 14.2, n_med = 1.2, n_ext = 1.1)

# test if a very narrow gate integrates close to the single point time domain solution
@test isapprox(g1_DA_semiinf_TD(τ,[1.101, 1.102], 1.0, 0.1, 10.0), g1_DA_semiinf_TD(τ,1.1015, 1.0, 0.1, 10.0), rtol = 1e-7)

#check corresponding g2 functions (these call g1 functions)
@test g2_DA_semiinf_TD(τ,[0.00001, 10.0], 1.8, 0.32, 4.1, BFi = 1.2e-7, n_ext = 1.8, n_med = 1.4, z = 0.0, λ = 820.1, β = 0.88, N_quad = 200) ≈ g2_DA_semiinf_CW(τ, 1.8, 0.32, 4.1, BFi = 1.2e-7, n_ext = 1.8, n_med = 1.4, z = 0.0, λ = 820.1, β = 0.88)
@test isapprox(g2_DA_semiinf_TD(τ,[2.101, 2.102], 0.4, 0.01, 12.0), g2_DA_semiinf_TD(τ,2.1015, 0.4, 0.01, 12.0), rtol = 1e-7)

### test layered against semi-inf

ρ = 1.0; μa = 0.1; μsp = 10.0; n_med = 1.0; n_ext = 1.0
β = 1.0; λ = 700.0; z = 0.0
BFi = 2.0e-8

si = g2_DA_semiinf_CW(τ, ρ, μa, μsp, BFi = BFi)

μa = [0.1, 0.1]; μsp = [10.0, 10.0]; n_med = [1.0, 1.0]; n_ext = 1.0
BFi = [2.0e-8, 2.0e-8]; l = [1.0, 10.0]; a = 25.0

data = Nlayer_cylinder_DCS(ρ = ρ, μa = μa, μsp = μsp, n_med = n_med, n_ext = n_ext, β = β, λ = λ, BFi = BFi, z = z, a = a, l = l, bessels = besselroots[1:2000])

@test g2_DA_Nlay_cylinder_CW(τ, data) ≈ g2_DA_Nlay_cylinder_CW(τ, ρ, μa, μsp; BFi = BFi, β = β, n_ext = n_ext, n_med = n_med, l = l, a = a, z = z, λ = λ, bessels = besselroots[1:2000])

layered = g2_DA_Nlay_cylinder_CW(τ, data)

@test si ≈ layered

# test layered TD solution for single time

τ = 10 .^(range(-10,stop=0,length=250))
ρ = 1.0; mua = 0.1; musp = 10.0; BFi = 2.0e-8; n_ext = 1.0; n_med = 1.0; z = 0.0; λ = 700.0
t = 1.0
si = g1_DA_semiinf_TD(τ, t, ρ, mua, musp; BFi = BFi, n_ext = n_ext, n_med = n_med, z = z, λ = λ)
layered = g1_DA_Nlay_cylinder_TD(τ, t, ρ, [mua, mua], [musp, musp]; BFi = [BFi, BFi], n_ext = n_ext, n_med = [n_med, n_med], z = z, λ = λ, bessels = besselroots[1:10], N_laplace = 12)

@test si ≈ layered

# test for different optical properties and SDS (limited by degeneration to semi-infinite solution have adjusted relative error)
τ = 10 .^(range(-10,stop=0,length=250))
ρ = 2.4; mua = 0.18; musp = 20.1; BFi = 1.2e-6; n_ext = 1.2; n_med = 1.1; z = 0.0; λ = 740.0
t = 2.1
si = g1_DA_semiinf_TD(τ, t, ρ, mua, musp; BFi = BFi, n_ext = n_ext, n_med = n_med, z = z, λ = λ)
layered = g1_DA_Nlay_cylinder_TD(τ, t, ρ, [mua, mua], [musp, musp]; BFi = [BFi, BFi], n_ext = n_ext, n_med = [n_med, n_med], z = z, λ = λ, bessels = besselroots[1:20], N_laplace = 12)

@test isapprox(si, layered, rtol = 1e-5)

# test layered TD solution for

# test range of time vector matches
τ = 10 .^(range(-10,stop=0,length=250))
ρ = 1.0; mua = 0.1; musp = 10.0; BFi = 2.0e-8; n_ext = 1.0; n_med = 1.0; z = 0.0; λ = 700.0
t = [1.0, 1.5]
si = g1_DA_semiinf_TD(τ, t, ρ, mua, musp; BFi = BFi, n_ext = n_ext, n_med = n_med, z = z, λ = λ)
layered = g1_DA_Nlay_cylinder_TD(τ, t, ρ, [mua, mua], [musp, musp]; BFi = [BFi, BFi], n_ext = n_ext, n_med = [n_med, n_med], z = z, λ = λ, bessels = besselroots[1:200], N_laplace = 8, N_quad = 100)
@test isapprox(si, layered, rtol = 1e-5)

# test range of time vector matches different optical properties
τ = 10 .^(range(-10,stop=0,length=250))
ρ = 2.2; mua = 0.42; musp = 18.2; BFi = 8.2e-8; n_ext = 1.2; n_med = 1.4; z = 0.0; λ = 700.0
t = [1.0, 2.6]
si = g1_DA_semiinf_TD(τ, t, ρ, mua, musp; BFi = BFi, n_ext = n_ext, n_med = n_med, z = z, λ = λ)
layered = g1_DA_Nlay_cylinder_TD(τ, t, ρ, [mua, mua], [musp, musp]; BFi = [BFi, BFi], n_ext = n_ext, n_med = [n_med, n_med], z = z, λ = λ, bessels = besselroots[1:200], N_laplace = 16, N_quad = 100)
@test isapprox(si, layered, rtol = 1e-5)

# test that the integral over whole time window range matches cw solution
τ = 10 .^(range(-10,stop=0,length=250))
ρ = 1.0; mua = 0.2; musp = 10.0; BFi = 3.0e-8; n_ext = 1.2; n_med = 1.4; z = 0.0; λ = 700.0
t = [1e-5, 10.0]
layered = g1_DA_Nlay_cylinder_TD(τ, t, ρ, [mua, mua], [musp, musp]; BFi = [BFi, BFi], n_ext = n_ext, n_med = [n_med, n_med], z = z, λ = λ, bessels = besselroots[1:700], N_laplace = 54, N_quad = 100)
si = g1_DA_semiinf_CW(τ, ρ, mua, musp; BFi = BFi, n_ext = n_ext, n_med = n_med, z = z, λ = λ)

@test isapprox(si, layered, rtol = 1e-5)

# test that g2 is working

τ = 10 .^(range(-10,stop=0,length=250))
ρ = 1.0; mua = 0.2; musp = 10.0; BFi = 3.0e-8; n_ext = 1.2; n_med = 1.4; z = 0.0; λ = 700.0
t = [1e-5, 10.0]
layered = g2_DA_Nlay_cylinder_TD(τ, t, ρ, [mua, mua], [musp, musp]; BFi = [BFi, BFi], n_ext = n_ext, n_med = [n_med, n_med], z = z, λ = λ, bessels = besselroots[1:700], N_laplace = 54, N_quad = 100)
si = g2_DA_semiinf_CW(τ, ρ, mua, musp; BFi = BFi, n_ext = n_ext, n_med = n_med, z = z, λ = λ)

@test isapprox(si, layered, rtol = 1e-5)

end


text:
@testset "util" begin
    include("util/simulate.jl")
    include("util/calibrate.jl")
end


text:
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	Scattering pattern liquid sphere
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function Pattern_LiquidSphere( rho_10::Any, k_0::Any, k_1::Any, a::Any, grid_size::Any,
					sum_size::Any, theta_inc::Any )
		#	Calculate the far-field sphere pattern in the liquid case.
		#
		#	input:	rho_10		: density quotient rho_1/rho_0
		#               k_0		: wavenumber in medium 0
		#               K_1		: wavenumber in medium 1
		#               a		: sphere radius
		#               grid_size	: angular grid points
		#               sum_size	: number of terms in the series solution
		#		theta_inc	: incidence direction
		#
		#	output:		Matrix of angles (in radians) (column 1) pattern (absolute value) (column 2)

		rho_10 = convert( AbstractFloat, rho_10 ) ;
		k_0 = convert( AbstractFloat, k_0 ) ;
		k_1 = convert( AbstractFloat, k_1 ) ;
		a = convert( AbstractFloat, a ) ;
		const Theta = linspace( 0, 2*pi , grid_size );
		const Pat = zeros( Complex128, grid_size, 1);
		const g = 1/rho_10 ; # rho/rho1 ( 'g^(-1)' with the usual 'g' in scattering theory )
		const x = k_0*a ;
		const x1 = k_1*a ;

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Main loop ( angles )
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		for j = 1 : grid_size
			acumulator = Complex(0);
			for i = 0 : sum_size
	   			Fn = Coef_LiquidSphere( i, x1 ) ;
				C = -( Sbessj(i,x)*g*Fn - x*Sbessj_p(i,x) )/
						( Sneum(i,x)*g*Fn - x*Sneum_p(i,x) );
				acumulator += (-1)^i*( 2*i + 1 )*( ( C - 1im*C*C )/( 1 + C*C ) )*
						sf_legendre_Pl( i , cos( - Theta[j] + theta_inc + pi ) )  ;
			end
			Pat[ j ] = -1/k_0*acumulator ; # Normalization
		end

		return [ Theta  map( Float64, abs( Pat ) ) ] ;
	end

	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	Prolate Spheroid coefficients (liquid)
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function Coeff_Pro( rho_10::Any, k_0::Any, k_1::Any, a::Any, b::Any, M::Any, method::Any, eta_inc::Any )

		#	Calculate the (Amn, Bmn) coefficients in the liquid case.
		#
		#	input:	rho_10	: density quotient rho_1/rho_0
		#               k_0	: wavenumber in medium 0
		#               K_1	: wavenumber in medium 1
		#               a	: spheroid major semiaxis
		#               b	: spheroid minor semiaxis
		#               M	: Max size of 'm' index
		#               method	: method of calculation. '1' for xi >> 1 and '2' for xi ~ 1
		#               eta_inc	: incidence angle
		#
		#	output:	Amn	: vector of coefficients
		#               Bmn	: vector of coefficients
		#

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Type conversion
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		convert( AbstractFloat, rho_10 ) ;
		convert( AbstractFloat, k_0 ) ;
		convert( AbstractFloat, k_1 ) ;
		convert( AbstractFloat, a ) ;
		convert( AbstractFloat, b ) ;
		convert( Int , M ) ;
		convert( Int , method ) ;
		convert( AbstractFloat, eta_inc ) ;

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Preliminaries
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		# Derived parameters

		const Size = round( Int64, (M+1)*(M+2)/2 ) ;
		const d = 2*sqrt( a^2 - b^2 ) ;
		const xi_0 = 1/sqrt( 1 - (b/a)^2 ) ; # Prolate
		const c_0 = ( d/2 )*k_0 ;
		const c_1 = ( d/2 )*k_1 ;

		# Structure declaration

		const Amn = zeros( Complex{BigFloat}, Size, 1) ;
		const Bmn = zeros( Complex{BigFloat}, Size, 1) ;
		const S1 = zeros( BigFloat, Size, 1) ;
		const S1_1 = zeros( BigFloat, Size, 1) ;
		const R1_0 = zeros( BigFloat, Size, 1) ;
		const R1p_0 = zeros( BigFloat, Size, 1) ;
		const R1_1 = zeros( BigFloat, Size, 1) ;
		const R1p_1 = zeros( BigFloat, Size, 1) ;
		const R3_0 = zeros( Complex{BigFloat}, Size, 1) ;
		const R3p_0 = zeros( Complex{BigFloat}, Size, 1) ;
		const R3_1 = zeros( Complex{BigFloat}, Size, 1) ;
		const R3p_1 = zeros( Complex{BigFloat}, Size, 1) ;
		const Alpha = zeros( BigFloat, M + 1, M + 1 ) ;

                const Norm0 = zeros(BigFloat, Size, 1) ; # 25-2-2016
		const Norm1 = zeros(BigFloat, Size, 1) ; # 25-2-2016

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Smn, Rmn calculation (Adelman-Gumerov-Duraiswami software)
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		for i = 1 : M + 1 # 'm' loop
			for j = i : M + 1 # 'n' loop
				indice = Index( M, i, j) ;
				pro_lambdamn_approx( c_0, i-1, j-1 ) ;
				run( `./pro_sphwv_S.R.sh $c_0 $(i-1) $(j-1) $eta_inc $xi_0` );
				R_0 = ReadFileToArrayBF( "Out_R.dat", ',', 0 )  ;
				S_0 = ReadFileToArrayBF( "Out_S.dat", ',', 0 )  ;
				pro_lambdamn_approx( c_1, i-1, j-1 ) ;
				run( `./pro_sphwv_S.R.sh $c_1 $(i-1) $(j-1) $eta_inc $xi_0` );
				R_1 = ReadFileToArrayBF( "Out_R.dat", ',', 0 )  ;
				S_1 = ReadFileToArrayBF( "Out_S.dat", ',', 0 )  ;
				S1[indice] = S_0[3] ;
				S1_1[indice] = S_1[3] ;
				if (method == 2)# 'Method 2'
					R1_0[indice] = R_0[6] ;
					R1p_0[indice] = R_0[7] ;
					R3_0[indice] = complex( R_0[6], R_0[12] ) ;
					R3p_0[indice] = complex( R_0[7], R_0[13] ) ;
					R1_1[indice] = R_1[6] ;
					R1p_1[indice] = R_1[7] ;
					R3_1[indice] = complex( R_1[6], R_1[12] ) ;
					R3p_1[indice] = complex( R_1[7], R_1[13] ) ;
				else # 'Method 1'
					R1_0[indice] = R_0[4] ;
					R1p_0[indice] = R_0[5] ;
					R3_0[indice]= complex( R_0[4], R_0[10] ) ;
					R3p_0[indice]= complex( R_0[5], R_0[11]) ;
					R1_1[indice]= R_1[4] ;
					R1p_1[indice]= R_1[5] ;
					R3_1[indice]= complex( R_1[4], R_1[10] ) ;
					R3p_1[indice]= complex( R_1[5], R_1[11] ) ;
				end
			end
		end

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Main algorithm ( Amn, Bmn calculation)
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		for m = 0 : M # Start with the M + 1 size matrix
			TAM = M - m + 1; # Structure's size for this 'm'
			# Structure declaration
			Q = zeros( Complex{BigFloat}, TAM, TAM );
			F = zeros( Complex{BigFloat}, 1, TAM );
			# Number of coefficients
			NRO = round(Int64,(M+1)*m -(m-1)*m/2 ) ; # 25-2-2016
			# Alpha matrix calculation
			for i = 1 : TAM # fill in the rows
				sigma = m + i - 1;
				name = @ sprintf("%08d_%03d_%03d", trunc(Int,c_0*1000), m, sigma );
				pro_lambdamn_approx( c_0, m, sigma ) ;
				run(`./pro_sphwv_dr.N.sh $c_0 $m $sigma $name`) ;
				# Norm0 = ReadFileToArrayBF( "Out_S_N.dat", '\n', 1) ;
				Norm0[ i + NRO ] = ReadFileToArrayBF( "Out_S_N.dat", '\n', 1) ;  # 25-2-2016
				DR_0 = ReadFileToArrayBF( "Out_dr.dat", '\n', 0) ;
				for j = 1 : TAM # fill in the columns
					n = m + j - 1;
					name = @ sprintf("%08d_%03d_%03d", trunc(Int,c_1*1000), m, n );
					pro_lambdamn_approx( c_1, m, n ) ;
					run(`./pro_sphwv_dr.N.sh $c_1 $m $n $name`) ;
					# Norm1 = ReadFileToArrayBF( "Out_S_N.dat", '\n', 1) ;
					Norm1[ j + NRO ] = ReadFileToArrayBF( "Out_S_N.dat", '\n', 1) ;  # 25-2-2016
					DR_1 = ReadFileToArrayBF( "Out_dr.dat", '\n', 0) ;
					size_dr = min( length(DR_1) , length(DR_0) ) ; # Size of dr matrix
					Alpha[i,j] = BigFloat(0) ;
					for k= 1 : size_dr
						r = k - 1;
						Alpha[i,j] += DR_0[k]*DR_1[k]*2*factorial( BigInt(2*(m)+r))/
								( ( 2*m + 2*r + 1 )*factorial( BigInt(r) ) );
					end
					Alpha[i,j] = Alpha[i,j]/sqrt(Norm0[ i + NRO ]*Norm1[ j + NRO ]); # 25-2-2016
				end
			end
			# Amn calculation
			S1_null_flag = 1; # Smn null for default
			for i = 1 : TAM
				if ( S1[NRO+i] != 0 )
					S1_null_flag = 0 ; # Smn is no null
					break ;
				end
			end
			if S1_null_flag == 0
				for j = 1 : TAM
					for i = 1 : TAM
						n = m + i - 1 ;
						F[j] += -1im^n*Alpha[i,j]*S1[NRO+i]*( rho_10*R1p_0[NRO+i]*R1_1[NRO+j] -
							R1p_1[NRO+j]*R1_0[NRO+i] );
					end
				end
				for i = 1 : TAM
					for j = 1 : TAM
						n = m + i - 1 ;
						Q[i,j] = 1im^n*Alpha[i,j]*S1[NRO+i]*( rho_10*R3p_0[NRO+i]*R1_1[NRO+j] -
							R1p_1[NRO+j]*R3_0[NRO+i] );
					end
				end
				if eta_inc == 0 # Q is singular
					for i = 1 : Int( floor(TAM/2) )
						Q[ 2*i, 2*i ] = 1 ;
					end
				end
 				Amn[ NRO+1 : NRO+TAM ] = \(Q',F') ;
			end
			# Bmn calculation
			for j = 1 : TAM
				if ( S1[ NRO + j ] == 0 )
					Bmn[ j + NRO ] = 0.0;
				else
					for i = 1 : TAM
						Bmn[ j + NRO ] += Alpha[i,j]*( 1/(1im) )^(i-j)*( S1[NRO+i]/S1[NRO+j] )*
								((R1_0[NRO+i] + Amn[ i+NRO ]*R3_0[NRO+i]))/R1_1[NRO+j]*
								sqrt(Norm0[ j + NRO ]/Norm1[ j + NRO ]); # 25-2-2016
					end
				end
			end
		end

		return Amn, Bmn ;
 	end

	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	Prolate spheroid far-field angular pattern
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function Pattern_LiquidPro( k_0::Any, a::Any , b::Any, M::Any,
					method::Any, eta_inc::Any, Amn::Array, deta::Any )

		#	Calculate the far field pattern in the liquid case.
		#
		#	input:	k_0	: wavenumber in medium 0
		#               a	: spheroid major semiaxis
		#               b	: spheroid minor semiaxis
		#               M	: Max size of 'm' index
		#               method	: method of calculation. '1' for xi >> 1 and '2' for xi ~ 1
		#		eta_inc	: incidence angle
		#               Amn	: vector of coefficients
		#		deta	: grid size in eta ( 2/deta = # of points in the angular pattern )
		#
		#	output:		the far-field absolute value angular pattern

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Preliminaries
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		# Parameter conversion
		convert( AbstractFloat, deta ) ;

		# Derived parameters
		const size = map( Int64, (M+1)*(M+2)/2 );
		const xi_0 = 1/sqrt( 1 - (b/a)^2 ) ;
		const d = 2*sqrt( a^2 - b^2 );
		const c_0 = (d/2)*k_0 ;

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		#	Smn calculations
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		const Smn = cell( 1, size ); # Structure for Smn(eta) indexed by 'idx'
		const Smn_inc = zeros( BigFloat, size, 1 ); # Vector for Smn(eta_inc)

		for i = 1 : M + 1 # 'm' loop
			for j = i : M + 1  # 'n' loop
 				idx = Index( M, i, j) ;
				# Calculate for all the grid and the incidence angle
				pro_lambdamn_approx( c_0, i-1, j-1 ) ;
				run(`./pro_sphwv_Sgrid.Sinc.sh $c_0 $(i-1) $(j-1) $eta_inc $deta`);
				if ( i==1 && j==1 ) # Building the grid (once)
					global ETA = reverse( readdlm( "Out_S.dat",',','\n')[:,2] ) ;
				end
				Smn[idx] = reverse( ReadFileToArrayBF2( "Out_S.dat", ',', :, 5 ) ) ;
				Smn_inc[ idx ] = ReadFileToArrayBF("Out_Sinc.dat", ',', 5 ) ;
			end
		end

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		#	Far-field pattern calculation
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		# Theta grid
		const ETAsize = length(ETA) ;
		const TETA = map( acos, ETA );
		const TETA2 = zeros( Float64, ETAsize - 1, 1 );
		for k = 1 : ETAsize - 1
			TETA2[k] = pi + TETA[k+1];
		end

		const Pat = zeros( Complex{BigFloat}, 2*ETAsize-1, 1);

		const Phi = [ 0, pi ]; # Only two values for Phi
		const Emn = cell( 1, length(Phi) ); # Emn: neumann_factor*cos( m(pi) )

		# Fill in Emn
		for k = 1 : length(Phi)
			temp = zeros( Float64, size, 1 ) ;
			for i = 1 : (M+1) # 'm' loop
				for j = i : (M+1) # 'n' loop
					idx = Index( M, i, j) ;
					if ( i==1 ) # m = 0 case
						temp[ idx ]= cos( (i-1)*Phi[k] ) ;
					else
						temp[ idx ]= 2*cos( (i-1)*Phi[k] ) ;
					end
				end
			end
			Emn[k] = d/c_0*temp ;
		end

		# Fill in the far-field pattern ( 0 - pi )
		for l = 1 : ETAsize
			for p = 1 : size
				Pat[ l ] += Emn[1][p]*Amn[p]*Smn[p][l]*Smn_inc[p];
			end
		end

		# Fill in the far-field pattern ( pi - 2*pi )
		for l = 1 : ETAsize - 1
			for p = 1 : size
				Pat[ ETAsize + l ] += Emn[2][p]*Amn[p]*Smn[p][ETAsize-l]*Smn_inc[p];
			end
		end

		# Outputting far-field absolute value
		return [ [TETA ; TETA2] map( Float64, abs( Pat ) ) ] ;
	end

	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	Oblate Spheroid coefficients (liquid)
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function Coeff_Obl( rho_10::Any, k_0::Any, k_1::Any, a::Any, b::Any, M::Any, method::Any, eta_inc::Any )

		#	Calculate the (Amn, Bmn) coefficients in the liquid case.
		#
		#	input:	rho_10	: density quotient rho_1/rho_0
		#               k_0	: wavenumber in medium 0
		#               K_1	: wavenumber in medium 1
		#               a	: spheroid major semiaxis
		#               b	: spheroid minor semiaxis
		#               M	: Max size of 'm' index
		#               method	: method of calculation. '1' for xi >> 1 and '2' for xi ~ 1
		#               eta_inc	: incidence angle
		#
		#	output:	Amn	: vector of coefficients
		#               Bmn	: vector of coefficients
		#

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Type conversion
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		convert( AbstractFloat, rho_10 ) ;
		convert( AbstractFloat, k_0 ) ;
		convert( AbstractFloat, k_1 ) ;
		convert( AbstractFloat, a ) ;
		convert( AbstractFloat, b ) ;
		convert( Int , M ) ;
		convert( Int , method ) ;
		convert( AbstractFloat, eta_inc ) ;

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Preliminaries
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		# Derived parameters

		const Size = round( Int64,(M+1)*(M+2)/2 ) ;
		const d = 2*sqrt( a^2 - b^2 ) ;
		const xi_0 = 1/sqrt( (a/b)^2 - 1 ) ; # Oblate
		const c_0 = (d/2)*k_0 ;
		const c_1 = (d/2)*k_1 ;

		# Structure declaration

		const Amn = zeros(Complex{BigFloat}, Size, 1) ;
		const Bmn = zeros(Complex{BigFloat}, Size, 1) ;
		const S1 = zeros(BigFloat, Size, 1) ;
		const S1_1 = zeros(BigFloat, Size, 1) ;
		const R1_0 = zeros(BigFloat, Size, 1) ;
		const R1p_0 = zeros(BigFloat, Size, 1) ;
		const R1_1 = zeros(BigFloat, Size, 1) ;
		const R1p_1 = zeros(BigFloat, Size, 1) ;
		const R3_0 = zeros(Complex{BigFloat}, Size, 1) ;
		const R3p_0 = zeros(Complex{BigFloat}, Size, 1) ;
		const R3_1 = zeros(Complex{BigFloat}, Size, 1) ;
		const R3p_1 = zeros(Complex{BigFloat}, Size, 1) ;
		const Alpha = zeros( BigFloat, M + 1, M + 1 ) ;

		const Norm0 = zeros(BigFloat, Size, 1) ; # 25-2-2016
		const Norm1 = zeros(BigFloat, Size, 1) ; # 25-2-2016

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Smn, Rmn calculation (Adelman-Duraiswami-Gumerov software)
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		for i = 1 : M + 1 # 'm' loop
			for j = i : M + 1 # 'n' loop
				indice = Index( M, i, j) ;
				obl_lambdamn_approx( c_0, i-1, j-1 ) ;
				run( `./obl_sphwv_S.R.sh $c_0 $(i-1) $(j-1) $eta_inc $xi_0` );
				R_0 = ReadFileToArrayBF( "Out_R.dat", ',', 0)  ;
				S_0 = ReadFileToArrayBF( "Out_S.dat", ',', 0)  ;
				obl_lambdamn_approx( c_1, i-1, j-1 ) ;
				run( `./obl_sphwv_S.R.sh $c_1 $(i-1) $(j-1) $eta_inc $xi_0` );
				R_1 = ReadFileToArrayBF( "Out_R.dat", ',', 0)  ;
				S_1 = ReadFileToArrayBF( "Out_S.dat", ',', 0)  ;
				S1[indice] = S_0[3] ;
				S1_1[indice] = S_1[3] ;
				if (method == 2) # 'Method 2'
					R1_0[indice] = R_0[5] ;
					R1p_0[indice] = R_0[6] ;
					R3_0[indice] = complex( R_0[5], R_0[11] ) ;
					R3p_0[indice] = complex( R_0[6], R_0[12] ) ;
					R1_1[indice] = R_1[5] ;
					R1p_1[indice] = R_1[6] ;
					R3_1[indice] = complex( R_1[5], R_1[11] ) ;
					R3p_1[indice] = complex( R_1[6], R_1[12] ) ;
				else # 'Method 1'
					R1_0[indice] = R_0[3] ;
					R1p_0[indice] = R_0[4] ;
					R3_0[indice]= complex( R_0[3], R_0[9] ) ;
					R3p_0[indice]= complex( R_0[4], R_0[10]) ;
					R1_1[indice]= R_1[3] ;
					R1p_1[indice]= R_1[4] ;
					R3_1[indice]= complex( R_1[3], R_1[9] ) ;
					R3p_1[indice]= complex( R_1[4], R_1[10] ) ;
				end
			end
		end

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Main algorithm ( Amn, Bmn calculation)
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		for m = 0 : M # Start with the M + 1 size matrix
			TAM = M - m + 1; # Structure's size for this 'm'
			# Structure declaration
			Q = zeros( Complex{BigFloat}, TAM, TAM );
			F = zeros( Complex{BigFloat}, 1, TAM );
			# Number of coefficients
			NRO = round(Int64,(M+1)*m -(m-1)*m/2 ) ; # 25-2-2016
			# Alpha matrix calculation
			for i = 1 : TAM # fill in the rows
				sigma = m + i - 1;
				name = @ sprintf("%08d_%03d_%03d", trunc(Int,c_0*1000), m, sigma );
				obl_lambdamn_approx( c_0, m, sigma ) ;
				run( `./obl_sphwv_dr.N.sh $c_0 $m $sigma $name` );
				# Norm0 = ReadFileToArrayBF( "Out_S_N.dat", '\n', 1) ;
				Norm0[ i + NRO ] = ReadFileToArrayBF( "Out_S_N.dat", '\n', 1) ;  # 25-2-2016
				DR_0 = ReadFileToArrayBF( "Out_dr.dat", '\n', 0) ;
				# Ciclo de 'j'
				for j = 1 : TAM # fill in the columns
					n = m + j - 1;
					name = @ sprintf("%08d_%03d_%03d", trunc(Int,c_1*1000), m, n );
					obl_lambdamn_approx( c_1, m, n ) ;
					run( `./obl_sphwv_dr.N.sh $c_1 $m $n $name` ) ;
					# Norm1 = ReadFileToArrayBF( "Out_S_N.dat", '\n', 1) ;
					Norm1[ j + NRO ] = ReadFileToArrayBF( "Out_S_N.dat", '\n', 1) ;  # 25-2-2016
					DR_1 = ReadFileToArrayBF( "Out_dr.dat", '\n', 0) ;
					size_dr = min( length(DR_1) , length(DR_0) );
					Alpha[i,j] = BigFloat(0) ;
					for k= 1 : size_dr
						r = k - 1;
						Alpha[i,j] += DR_0[k]*DR_1[k]*2*factorial( BigInt(2*(m)+r))/
								( ( 2*m + 2*r + 1 )*factorial( BigInt(r) ) );
					end
					Alpha[i,j] = Alpha[i,j]/sqrt(Norm0[ i + NRO ]*Norm1[ j + NRO ]);  # 25-2-2016
				end
			end
			# Amn calculation
			S1_null_flag = 1;  # Smn null for default
			for i = 1 : TAM
				if ( S1[NRO+i] != 0 ) # Smn is not null
					S1_null_flag = 0 ;
					break ;
				end
			end
			if S1_null_flag == 0
				for j = 1 : TAM
					for i = 1 : TAM
						n = m + i - 1 ;
						F[j] += -1im^n*Alpha[i,j]*S1[NRO+i]*( rho_10*R1p_0[NRO+i]*R1_1[NRO+j] -
							R1p_1[NRO+j]*R1_0[NRO+i] );
					end
				end
				for i = 1 : TAM
					for j = 1 : TAM
						n = m + i - 1 ;
						Q[i,j] = 1im^n*Alpha[i,j]*S1[NRO+i]*( rho_10*R3p_0[NRO+i]*R1_1[NRO+j] -
							R1p_1[NRO+j]*R3_0[NRO+i] );
					end
				end
				if eta_inc == 0 # Q is singular
					for i = 1 : Int( floor(TAM/2) )
						Q[ 2*i, 2*i ] = 1 ;
					end
				end
 				Amn[ NRO+1 : NRO+TAM ] = \(Q',F') ;
			end
			# Bmn calculation
			for j = 1 : TAM
				if ( S1[ NRO + j ] == 0 )
					Bmn[ j + NRO ] = 0.0;
				else
					for i = 1 : TAM
						Bmn[ j + NRO ] += Alpha[i,j]*( 1/(1im) )^(i-j)*( S1[NRO+i]/S1[NRO+j] )*
								((R1_0[NRO+i] + Amn[ i+NRO ]*R3_0[NRO+i]))/R1_1[NRO+j]*
								sqrt(Norm0[ j + NRO ]/Norm1[ j + NRO ]); # 25-2-2016
					end
				end
			end
		end

		return Amn, Bmn ;
 	end

	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	Oblate spheroid far-field angular pattern
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function Pattern_LiquidObl( k_0::Any, a::Any , b::Any, M::Any,
					method::Any, eta_inc::Any, Amn::Array, deta::Any )
		#	Calculate the far field pattern in the liquid case.
		#
		#	input:	k_0	: wavenumber in medium 0
		#               a	: spheroid major semiaxis
		#               b	: spheroid minor semiaxis
		#               M	: Max size of 'm' index
		#               method	: method of calculation. '1' for xi >> 1 and '2' for xi ~ 1
		#		eta_inc	: incidence angle
		#               Amn	: vector of coefficients
		#		deta	: grid size in eta ( 2/deta = # of points in the angular pattern )
		#
		#	output: 	the far-field absolute value angular pattern

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Preliminaries
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		# Parameter conversion
		convert( AbstractFloat, deta ) ;

		# Derived parameters
		const size = map( Int64, (M+1)*(M+2)/2 );
		const xi_0 = 1/sqrt( (a/b)^2 - 1 ) ;
		const d = 2*sqrt( a^2 - b^2 );
		const c_0 = (d/2)*k_0 ;

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		#	Smn calculations
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		const Smn = cell( 1, size ); # Structure for Smn(eta) indexed by 'idx'
		const Smn_inc = zeros( BigFloat, size, 1 ); # Vector for Smn(eta_inc)

		for i = 1 : M + 1 # 'm' loop
			for j = i : M + 1  # 'n' loop
 				idx = Index( M, i, j) ;
				# Calculate for all the grid and the incidence angle
				obl_lambdamn_approx( c_0, i-1, j-1 ) ;
				run( `./obl_sphwv_Sgrid.Sinc.sh $c_0 $(i-1) $(j-1) $eta_inc $deta` );
				if ( i==1 && j==1 ) # Building the grid (once)
					global ETA = reverse( readdlm( "Out_S.dat",',','\n')[:,2] ) ; # columna 2
				end
				Smn[idx] = reverse( ReadFileToArrayBF2( "Out_S.dat", ',', :, 5 ) ) ; # Levantamos Smn calculada por el método '2'. Invierto el orden
				Smn_inc[ idx ] = ReadFileToArrayBF("Out_Sinc.dat", ',', 5 ) ;
			end
		end

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		#	Far-field pattern calculation
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		# Theta grid
		const ETAsize = length(ETA) ;
		const TETA = map( acos, ETA );
		const TETA2 = zeros( Float64, ETAsize - 1, 1 );
		for k = 1 : ETAsize - 1
			TETA2[k] = pi + TETA[k+1];
		end

		const Pat = zeros( Complex{BigFloat}, 2*ETAsize-1, 1);

		const Phi = [ 0, pi ];  # Only two values for Phi
		const Emn = cell( 1, length(Phi) ); # Emn: neumann_factor*cos( m(pi) )

 		# Fill in Emn
		for k = 1 : length(Phi)
			temp = zeros( Float64, size, 1 );
			for i = 1 : (M+1) # 'm' loop
				for j = i : (M+1) # 'n' loop
					idx = Index( M, i, j) ;
					if ( i==1 ) # m = 0 case
						temp[ idx ]= cos( (i-1)*Phi[k] );
					else
						temp[ idx ]= 2*cos( (i-1)*Phi[k] );
					end
				end
			end
			Emn[k] = d/c_0*temp;
		end

		# Fill in the far-field pattern ( 0 - pi )
		for l = 1 : ETAsize
			for p = 1 : size
				Pat[ l ] += Emn[1][p]*Amn[p]*Smn[p][l]*Smn_inc[p];
			end
		end

		# Fill in the far-field pattern ( pi - 2*pi )
		for l = 1 : ETAsize - 1
			for p = 1 : size
				Pat[ ETAsize + l ] += Emn[2][p]*Amn[p]*Smn[p][ETAsize-l]*Smn_inc[p];
			end
		end

		# Outputting far-field absolute value
		return [ [TETA ; TETA2] map( Float64, abs( Pat ) ) ] ;
	end

	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	Cálculo de coeficientes Prolate Spheroid (Dirichlet-Neumann)
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function Coeff_ProDN( k_0::Any, a::Any, b::Any, M::Any, method::Any, eta_inc::Any , flag::AbstractString )
		#	Calculate the Amn coefficients in the Dirichlet-Neumann case.
		#
		#	input:	k_0	: wavenumber in medium 0
		#               a	: spheroid major semiaxis
		#               b	: spheroid minor semiaxis
		#               M	: Max size of 'm' index
		#               method	: method of calculation. '1' for xi >> 1 and '2' for xi ~ 1
		#               eta_inc	: incidence angle
		#		flag	: 'D' or 'N' for Dirichlet o neumann boundary conditions
		#
		#	output:		vector of coefficients
		#

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Type conversion
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		convert( AbstractFloat, k_0 ) ;
		convert( AbstractFloat, a ) ;
		convert( AbstractFloat, b ) ;
		convert( Int , M ) ;
		convert( Int , method ) ;
		convert( AbstractFloat, eta_inc ) ;

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Preliminaries
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		# Derived parameters
		const Size = round( Int64, (M+1)*(M+2)/2 );
		const d = 2*sqrt( a^2 - b^2 ) ;
		const xi_0 = 1/sqrt( 1 - (b/a)^2 ) ; # Prolate
		const c_0 = ( d/2 )*k_0 ;

		# Structure declaration
		const Amn = zeros( Complex{BigFloat}, Size, 1) ;
		const R1_0 = zeros( BigFloat, Size, 1) ;
		const R1p_0 = zeros( BigFloat, Size, 1) ;
		const R3_0 = zeros( Complex{BigFloat}, Size, 1) ;
		const R3p_0 = zeros( Complex{BigFloat}, Size, 1) ;

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Rmn calculation (Adelman-Duraiswami-Gumerov software)
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		for i = 1 : M + 1 # 'm' loop
			for j = i : M + 1 # 'n' loop
				indice = Index( M, i, j) ;
				pro_lambdamn_approx( c_0, i-1, j-1 ) ;
				run(`./pro_sphwv_S.R.sh $c_0 $(i-1) $(j-1) $eta_inc $xi_0`);
				R_0 = ReadFileToArrayBF( "Out_R.dat", ',', 0)  ;
				if (method == 2)# 'Method 2'
					R1_0[indice] = R_0[6] ;
					R1p_0[indice] = R_0[7] ;
					R3_0[indice] = complex( R_0[6], R_0[12] ) ;
					R3p_0[indice] = complex( R_0[7], R_0[13] ) ;
				else # 'Method 1'
					R1_0[indice] = R_0[4] ;
					R1p_0[indice] = R_0[5] ;
					R3_0[indice]= complex( R_0[4], R_0[10] ) ;
					R3p_0[indice]= complex( R_0[5], R_0[11]) ;
				end
			end
		end

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Amn calculation
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		if ( flag == "N" ) # Neumann
			for i = 1 : M + 1
				for j = i : M + 1
					indice = Index( M, i, j )  ;
					if ( R3p_0[indice] != 0 )
						Amn[indice] = - R1p_0[indice]/R3p_0[indice] ;
					end
				end
			end
		elseif ( flag == "D" )	# Dirichlet
			for i = 1 : M + 1
		 		for j = i : M + 1
		 			indice = Index( M, i, j ) ;
		 			if ( R3_0[indice] !=0 )
		 				Amn[indice] = - R1_0[indice]/R3_0[indice];
		 			end
		 		end
		 	end
		else
			println("ERROR: Boundary conditions not defined.")
		end

		return Amn ;
 	end

	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	Oblate Spheroid coefficientes (Dirichlet-Neumann)
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function Coeff_OblDN( k_0::Any, a::Any, b::Any, M::Any, method::Any, eta_inc::Any, flag::AbstractString  )
		#	Calculate the Amn coefficients in the Dirichlet-Neumann case.
		#
		#	input:	k_0	: wavenumber in medium 0
		#               a	: spheroid major semiaxis
		#               b	: spheroid minor semiaxis
		#               M	: Max size of 'm' index
		#               method	: method of calculation. '1' for xi >> 1 and '2' for xi ~ 1
		#               eta_inc	: incidence angle
		#		flag	: 'D' or 'N' for Dirichlet o neumann boundary conditions
		#
		#	output:	Amn	: Vector of coefficients
		#

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Type conversion
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		convert( AbstractFloat, k_0 ) ;
		convert( AbstractFloat, a ) ;
		convert( AbstractFloat, b ) ;
		convert( Int , M ) ;
		convert( Int , method ) ;
		convert( AbstractFloat, eta_inc ) ;

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Preliminaries
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		# Derived parameters
		const Size = round( Int64,(M+1)*(M+2)/2 );
		const d = 2*sqrt(a^2 - b^2) ;	# Oblate
		const xi_0 = 1/sqrt((a/b)^2-1) ;
		const c_0 = (d/2)*k_0 ;

		# Structure declaration
		const Amn = zeros(Complex{BigFloat}, Size, 1) ;
		const R1_0 = zeros(BigFloat, Size, 1) ;
		const R1p_0 = zeros(BigFloat, Size, 1) ;
		const R3_0 = zeros(Complex{BigFloat}, Size, 1) ;
		const R3p_0 = zeros(Complex{BigFloat}, Size, 1) ;

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Rmn calculation (Adelman-Duraiswami-Gumerov software)
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		for i = 1 : M + 1 # 'm' loop
			for j = i : M + 1 # 'n' loop
				indice = Index( M, i, j) ;
				obl_lambdamn_approx( c_0, i-1, j-1 ) ;
				run( `./obl_sphwv_S.R.sh $c_0 $(i-1) $(j-1) $eta_inc $xi_0` );
				R_0 = ReadFileToArrayBF( "Out_R.dat", ',', 0)  ;
				if (method == 2)# 'Method 2'
					R1_0[indice] = R_0[5] ;
					R1p_0[indice] = R_0[6] ;
					R3_0[indice] = complex( R_0[5], R_0[11] ) ;
					R3p_0[indice] = complex( R_0[6], R_0[12] ) ;
				else # 'Method 1'
					R1_0[indice] = R_0[3] ;
					R1p_0[indice] = R_0[4] ;
					R3_0[indice]= complex( R_0[3], R_0[9] ) ;
					R3p_0[indice]= complex( R_0[4], R_0[10]) ;
				end
			end
		end

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# 	Amn Calculation
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		if ( flag == "N" ) # Neumann
			for i = 1 : M + 1
				for j = i : M + 1
					indice = Index( M, i, j )  ;
					if ( R3p_0[indice] != 0 )
						Amn[indice] = - R1p_0[indice]/R3p_0[indice] ;
					end
				end
			end
		elseif ( flag == "D" ) # Dirichlet
			for i = 1 : M + 1
		 		for j = i : M + 1
		 			indice = Index( M, i, j ) ;
		 			if ( R3_0[indice] !=0 )
		 				Amn[indice] = - R1_0[indice]/R3_0[indice];
		 			end
		 		end
		 	end
		else
			println("ERROR: Boundary conditions not defined.")
		end

		return Amn ;
	end


text:
function sum_3_5(n::Int)
    s = 0
    for i = 3:3:n-1
        s += i
    end

    for i = 5:5:n-1
        if i% 3 != 0
            s += i
        end
    end
    return s
end
print(sum_3_5(1000))


text:
module VisJS

# package code goes here
export VisJSGraph, NodeType

type NodeType
  id::Int64
  label::ASCIIString
  color::ASCIIString
  group::ASCIIString
end

type VisJSGraph
    n::Array{NodeType,1}
    e
    divid
end

VisJSGraph(n,e) = VisJSGraph(n, e, randstring())

include("display.jl")
include("utils.jl")

if Pkg.installed("Escher") != nothing
       include("escher_integration.jl")
end

end # module


text:
function Ybus_current_kirchoff(inputs, V_r, V_i, I_injections_r, I_injections_i)
    I_bus = get_aux_arrays(inputs)[5]
    I_balance = get_aux_arrays(inputs)[6]
    Ybus = get_Ybus(inputs) + get_total_shunts(inputs)
    # Note: BLAS doesn't work because the the type of Vr and Vi is not Matrix of Complex
    LinearAlgebra.mul!(I_bus, Ybus, (V_r + V_i .* 1im))
    bus_count = length(I_bus)
    for n in eachindex(I_bus)
        I_balance[n] = I_injections_r[n] - real(I_bus[n])
        I_balance[n + bus_count] = I_injections_i[n] - imag(I_bus[n])
    end

    return I_balance
end


text:
module cindex

export parse_header, cu_type, ty_kind, name, spelling,
       value, children, cu_file, resolve_type, return_type,
       tokenize, pointee_type, typedef_type
export CLType, CLCursor, CXString, CXTypeKind, CursorList, TokenList

import Base.search
export search

export function_return_modifiers, function_arg_modifiers, cu_visitor_cb

import ..libclang

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

include("cindex/defs.jl")
include("cindex/types.jl")
include("cindex/base.jl")

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

# Main entry point for parsing
# Returns root CXCursor in the TranslationUnit for a given header
#
# Required argument:
#   "header.h"          header file to parse
#
# Optional (keyword) arguments:
#   index:              CXIndex pointer (pass to avoid re-allocation)
#   diagnostics:        Display Clang diagnostics
#   cplusplus:          Parse as C++
#   includes:           Vector of extra include directories to search
#   args:               Compiler switches as string array, eg: ["-x", "c++", "-fno-elide-type"]
#   flags:              Bitwise OR of TranslationUnitFlags
#
function parse_header(header::AbstractString;
                index                           = Union{},
                diagnostics::Bool               = true,
                cplusplus::Bool                 = false,
                args                            = [""],
                includes                        = String[],
                flags                           = TranslationUnit_Flags.None)
    if !isfile(header)
        error(header, " not found")

        # TODO: support parsing in-memory with CXUnsavedFile arg
        #       to _parseTranslationUnit.jj
        #unsaved_file = CXUnsavedFile()
    end
    if (index == Union{})
        index = idx_create(0, (diagnostics ? 1 : 0))
    end
    if (cplusplus)
        append!(args, ["-x", "c++"])
    end
    if (length(includes) > 0)
        args = vcat(args, [["-I",x] for x in includes]...)
    end

    tu = tu_parse(index, header, args, length(args),
                  C_NULL, 0, flags)
    if (tu == C_NULL)
        error("ParseTranslationUnit returned NULL; unable to create TranslationUnit for header $(header)")
    end

    return tu_cursor(tu)
end

###############################################################################
# Search function for CursorList
# Returns vector of CXCursors in CursorList matching predicate
#
# Required arguments:
#   CursorList      List to search
#   IsMatch(CXCursor)
#                   Predicate Function, accepting a CXCursor argument
#
function search(cl::Vector{CLCursor}, ismatch::Function)
    ret = CLCursor[]
    for cu in cl
        ismatch(cu) && push!(ret, cu)
    end
    ret
end
search(cu::CLCursor, ismatch::Function) = search(children(cu), ismatch)
search(cu::CLCursor, T::DataType) = search(cu, x->isa(x, T))
search(cu::CLCursor, name::String) = search(cu, x->(cindex.spelling(x) == name))

###############################################################################
# Extended search function
# Returns a Dict{ DataType => CLCursor

function matchchildren(cu::CLCursor, types::Array{DataType,1})
    ret = Dict{Any,Any}([(t, CLCursor[]) for t in types])
    for child in children(cu)
        for t in types
            isa(child, t) && push!(ret[t], child)
        end
    end
    return ret
end

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

# TODO: macro version should be more efficient.
anymatch(first, args...) = any([==(first, a) for a in args])

cu_type(c::CLCursor) = getCursorType(c)
ty_kind(c::CLType) = convert(Int, c.typ.kind)
name(c::CLCursor) = getCursorDisplayName(c)
spelling(c::CLType) = getTypeKindSpelling(convert(Int32, ty_kind(c)))
spelling(c::CLCursor) = getCursorSpelling(c)
is_null(c::CLCursor) = (Cursor_isNull(c) != 0)

function resolve_type(rt::CLType)
    # This helper attempts to work around some limitations of the
    # current libclang API.
    if ty_kind(rt) == cindex.TypeKind.Unexposed
        # try to resolve Unexposed type to cursor definition.
        rtdef_cu = cindex.getTypeDeclaration(rt)
        if (!is_null(rtdef_cu) && !isa(rtdef_cu, NoDeclFound))
            return cu_type(rtdef_cu)
        end
    end
    # otherwise, this will either be a builtin or unexposed
    # client needs to sort out.
    return rt
end

function return_type(c::Union{FunctionDecl, CXXMethod}, resolve::Bool)
    if (resolve)
        return resolve_type( getCursorResultType(c) )
    else
        return getCursorResultType(c)
    end
end
return_type(c::CLCursor) = return_type(c, true)

function pointee_type(t::Pointer)
    return cindex.getPointeeType(t)
end
pointee_type(cu::CLCursor) = error("pointee_type(CLCursor) is discontinued, please use pointee_type(cindex.cu_type(cu))")

function typedef_type(c::TypedefDecl)
    t = cindex.getTypedefDeclUnderlyingType(c)
    if isa(t, Unexposed)
        t = cu_type(children(c)[1])
    end
    return t
end

function value(c::EnumConstantDecl)
    t = cu_type(c)
    if isa(t, IntType) || isa(t, Long) || isa(t, LongLong)
            return getEnumConstantDeclValue(c)
    elseif isa(t, UInt) || isa(t, ULong) || isa(t, ULongLong)
            return getEnumConstantDeclUnsignedValue(c)
    end
    error("Unknown EnumConstantDecl type: ", t, " cursor: ", c)
end

function cu_children_visitor(cursor::CXCursor, parent::CXCursor, client_data::Ptr{Cvoid})
    list = unsafe_pointer_to_objref(client_data)
    push!(list, cursor)
    return Cuint(1) # CXChildVisit_Continue TODO: use enum
end

# TODO implement -- reduce allocations
#=
function cu_children_search(cursor::CXCursor, parent::CXCursor, client_data::Ptr{Cvoid})
    # expect ("target", store::Array{CXCursor,1}, findfirst::Bool)
    target,store,findfirst = unsafe_pointer_to_objref(client_data)
end
=#

function children(cu::CLCursor)
    # TODO: possible to use sizehint! here?
    list = CLCursor[]
    ccall( (:clang_visitChildren, libclang), Nothing,
           (CXCursor, Ptr{Cvoid}, Ptr{Cvoid}),
           cu, cu_visitor_cb, pointer_from_objref(list))
    list
end

# TODO add cu_location
function cu_file(cu::CLCursor)
    # TODO should return a struct or namedtuple
    # TODO turn this in to a normal wrapper
    loc = ccall( (:clang_getCursorLocation, libclang), CXSourceLocation,
                 (CXCursor,), cu)

    cxfile = Ref{CXFile}()
    line = Ref{Cuint}(); col = Ref{Cuint}(); offset = Ref{Cuint}()

    ccall( (:clang_getExpansionLocation, libclang), (Nothing),
            (CXSourceLocation, Ref{CXFile}, Ref{Cuint}, Ref{Cuint}, Ref{Cuint}),
            loc,               cxfile,       line,       col,        offset)

    getFileName(cxfile[])
end

################################################################################
# Tokenizer access
################################################################################

# Returns TokenList
function tokenize(cursor::CLCursor)
    tu = Cursor_getTranslationUnit(cursor)
    sourcerange = getCursorExtent(cursor)
    return cindex.tokenize(tu, sourcerange)
end

Base.firstindex(tl::TokenList) = 1
Base.lastindex(tl::TokenList) = length(tl)
Base.length(tl::TokenList) = tl.size
Base.iterate(tl::TokenList, state=1) = state > tl.size ? nothing : (tl[state], state+1)

function Base.getindex(tl::TokenList, i::Int)
    if (i < 1 || i > length(tl))
        # throw(BoundsError())
        @warn "TokenList miss at index ", i
        return Comment(string("TokenList miss at index ", i))
    end

    c = unsafe_load(tl.ptr, i)
    kind = c.int_data[1]
    spelling = cindex.getTokenSpelling(tl.tunit, c)

    if (kind == TokenKind.Punctuation)
        return Punctuation(spelling)
    elseif (kind == TokenKind.Keyword)
        return Keyword(spelling)
    elseif (kind == TokenKind.Identifier)
        return Identifier(spelling)
    elseif (kind == TokenKind.Literal)
        return Literal(spelling)
    elseif (kind == TokenKind.Comment)
        return Comment(spelling)
    end
    throw("Unknown TokenKind: $kind")
end

################################################################################
# Helper functions
################################################################################

# Retrieve function arguments for a given cursor
function_args(cursor::Union{FunctionDecl, CXXMethod}) = search(cursor, ParmDecl)

#returns a tuple with the default argument values for a C++ function
#only seems to work if cplusplus=true in parse_header
function function_arg_defaults(method::Union{cindex.CXXMethod, cindex.FunctionDecl, cindex.Constructor})
    defvals = Any[]
    for c in children(method)
        if isa(c,cindex.ParmDecl)

            #value in case default was not specified
            lit = nothing
            toks = tokenize(c)

            n = 1
            #get index of '='
            for tok in toks
                if isa(tok, cindex.Punctuation) && tok.text == "="
                    break
                end
                n += 1
            end

            ts = collect(toks)

            #default value is the one after '='
            if n<length(toks)
                lit = join(map(x->x.text, ts[n+1:end]), "")
            end

            #parse default value as Julia code
            pl = if (lit != nothing)
                    try
                        Meta.parse(lit)
                    catch err
                        @info "Error parsing function_default_arg value: $lit"
                        return tuple(Any[])
                    end
                 else
                     lit
                 end

            #for empty string, need to add quotes
            if pl == ""
                pl = "\"\""
            end
            push!(defvals, pl)
        end
    end
    return tuple(defvals...)
end

#returns an array of the function return type modifiers
function function_return_modifiers(f::Union{FunctionDecl, CXXMethod})
    modifs = Any[]
    for tok in tokenize(f)
        if isa(tok, cindex.Keyword)
            if tok.text in ["const"]
                push!(modifs, tok)
            end
        end
        isa(tok, cindex.Identifier) && tok.text==spelling(f) && break
    end
    return modifs
end

function function_arg_modifiers(p::ParmDecl)
    modifs = cindex.Keyword[]
    for tok in tokenize(p)

        #only read keywords
        isa(tok, cindex.Keyword) || continue

        #read up to variable identifier
        isa(tok, cindex.Identifier) && break

        if tok.text in ["const", "virtual"]
            push!(modifs, tok)
        end
    end
    return modifs
end

################################################################################
# Display overrides
################################################################################

Base.show(io::IO, tk::CLToken)  = print(io, typeof(tk), "(\"", tk.text, "\")")
Base.show(io::IO, ty::CLType)   = print(io, "CLType (", typeof(ty), ") ")
Base.show(io::IO, cu::CLCursor) = print(io, "CLCursor (", typeof(cu), ") ", name(cu))

###############################################################################
# Internal functions
###############################################################################

tu_init(hdrfile::Any) = tu_init(hdrfile, 0, false, 0)
function tu_init(hdrfile::Any, diagnostics, cpp::Bool, opts::Int)
    idx = idx_create(0,diagnostics)
    tu = tu_parse(idx, hdrfile, (cpp ? ["-x", "c++"] : [""]), opts)
    return tu
end

function tu_dispose(tu::CXTranslationUnit)
    ccall( (:clang_disposeTranslationUnit, libclang), Nothing, (Ptr{Cvoid},), tu)
end

function tu_cursor(tu::CXTranslationUnit)
    tu == C_NULL && error("Invalid TranslationUnit!")
    getTranslationUnitCursor(tu)
end

function tu_parse(CXIndex,
                  source_filename::AbstractString,
                  cl_args::Array{S,1},
                  num_clargs,
                  unsaved_files::CXUnsavedFile,
                  num_unsaved_files,
                  options) where S <: String

    ccall(  (:clang_parseTranslationUnit, libclang),
            CXTranslationUnit,
            (Ptr{Cvoid}, Ptr{UInt8}, Ptr{Ptr{UInt8}}, UInt32, Ptr{Cvoid}, UInt32, UInt32),
                CXIndex,
                source_filename,
                cl_args,
                num_clargs,
                unsaved_files,
                num_unsaved_files,
                options)
end

function idx_create(excludeDeclsFromPCH::Int, displayDiagnostics::Int)
    ccall(  (:clang_createIndex, libclang),
            CXTranslationUnit,
            (Int32, Int32),
                excludeDeclsFromPCH,
                displayDiagnostics)
end
idx_create() = idx_create(0,0)

function getFile(tu::CXTranslationUnit, file::String)
    ccall( (:clang_getFile, libclang),
            CXFile,
            (Ptr{Cvoid}, Ptr{UInt8}),
                tu,
                file)
end

function __init__()
    global cu_visitor_cb = @cfunction(cu_children_visitor, Cuint, (CXCursor, CXCursor, Ptr{Cvoid}))
end

end # module


text:
Base.show(io::IO, ::MIME"text/plain", ir::StaticIR) =
    print_ir(io, ir)

function print_ir(io::IO, ir::StaticIR)
    println(io, "== Static IR ==")
    args = join((string(arg.name) for arg in ir.arg_nodes), ", ")
    println(io, "Arguments: ($args)")
    for node in ir.nodes
        node in ir.arg_nodes && continue
        print(io, "  "); print_ir(io, node); println(io)
    end
    print(io, "  return $(ir.return_node.name)")
end

function print_ir(io::IO, node::TrainableParameterNode)
    print(io, "@param $(node.name)::$(node.typ)")
end

function print_ir(io::IO, node::JuliaNode)
    inputs = join((string(i.name) for i in node.inputs), ", ")
    print(io, "$(node.name) = $(node.fn)($inputs)")
end

function print_ir(io::IO, node::GenerativeFunctionCallNode)
    inputs = join((string(i.name) for i in node.inputs), ", ")
    gen_fn_name = ir_name(node.generative_function)
    print(io, "$(node.name) = @trace($(gen_fn_name)($inputs), :$(node.addr))")
end

function print_ir(io::IO, node::RandomChoiceNode)
    inputs = join((string(i.name) for i in node.inputs), ", ")
    print(io, "$(node.name) = @trace($(node.dist)($inputs), :$(node.addr))")
end

ir_name(fn::GenerativeFunction) = nameof(typeof(fn))
ir_name(fn::DynamicDSLFunction) = nameof(fn.julia_function)


text:
Random.seed!(17)
N = 100
net = CompleteBipartiteGraph(N, N, Network)
eprop!(net, "w", e -> rand())
vprop!(net, "b", v -> 1)
res_match = Matching.run_bp(net, maxiters=100, verbose=false)
res_bmatch = BMatching.run_bp(net, maxiters=100, verbose=false)
bmatchmap = [m[1] for m in res_bmatch.match]
@test all(bmatchmap .== res_match.match)

Random.seed!(17)
N = 100
net = CompleteGraph(N, Network)
eprop!(net, "w", e -> rand())
vprop!(net, "b", v -> rand(1:5))
res = BMatching.run_bp(net, maxiters=100, verbose=false)
@test res.num_violations <= 1 # = 0 on Julia 1.6
@test res.iters <= 100


text:
using DSP, Test
using DelimitedFiles: readdlm
using Polynomials: Polynomial

#######################################
#
#  http://www.mathworks.com.au/help/signal/ref/freqz.html
#  Example 1: Frequency response from transfer function
#
#  dlmwrite('freqz-eg1.txt',[w, abs(h)], 'delimiter', '\t', 'precision', '%.12f')
#
#######################################
@testset "freqresp from TF" begin
    # Matlab
    freqz_eg1_w_abs = readdlm(joinpath(dirname(@__FILE__), "data", "freqz-eg1.txt"),'\t')
    matlab_abs   = freqz_eg1_w_abs[:,2]

    # Julia
    b0 = 0.05634
    b1 = [1,  1]
    b2 = [1, -1.0166, 1]
    a1 = [1, -0.683]
    a2 = [1, -1.4461, 0.7957]
    b = b0*conv(b1, b2)
    a = conv(a1, a2)

    w     = range(0, stop=6.280045284537, length=2001)
    h     = freqresp(PolynomialRatio(b, a), w)
    h_abs = convert(Array{Float64}, abs.(h))
    h_ref = [b0 * Polynomial(b1)(z) * Polynomial(b2)(z) / (Polynomial(a1)(z) * Polynomial(a2)(z)) for z ∈ exp.(-im*w)]

    # Test
    @test h_abs ≈ matlab_abs
    @test h ≈ h_ref
end

#######################################
#
# freqresp with conversion to PolynomialRatio may lead to undesirable roundoff errors
# check by using the same poles and zeros such that they should cancel
#
#######################################
@testset "freqresp pole/zero cancellation" begin
    rs = [1-10.0^-n for n in 1:15]
    @test freqresp(ZeroPoleGain(rs, reverse(rs), 42.0), range(0, stop=2π, length=50)) ≈ fill(42., 50)

    biq = Biquad(42.0, 84.0*real(.999999*exp(im)), 42.0*.999999^2, 2.0*real(.999999*exp(im)), .999999^2)
    @test freqresp(biq, range(0, stop=2π, length=50)) ≈ fill(42.0, 50)

    pr = [1-10.0^-n for n in 1:10]
    zr = reverse(pr)
    bs = [Biquad(42.0, 84.0*real(r[1]*exp(im)), 42.0*r[1]^2, 2.0*real(r[2]*exp(im)), r[2]^2) for r in zip(zr, pr)]
    sos = SecondOrderSections(bs, 42.0^-9)
    @test freqresp(sos, range(0, stop=2π, length=50)) ≈ fill(42.0, 50)

    @test freqresp(ZeroPoleGain(ComplexF64[], ComplexF64[], 42.0), range(0, stop=2π, length=50)) == fill(42.0, 50)
    @test freqresp(SecondOrderSections(Biquad{:z,Float64}[], 42.0), range(0, stop=2π, length=50)) == fill(42.0, 50)
end

#######################################
#
#  Test frequency, phase, impulse and step response
#
#  Data from Matlab using b,a and a from above:
#  ir = impz(b,a, 512);
#  stepr = stepz(b,a, 512);
#  h = abs(freqz(b,a));
#  [phi,f] =  phasez(b, a);
#  all = [f, ir, stepr, h, phi];
#  dlmwrite('responses-eg1.txt',all, 'delimiter', '\t', 'precision', '%.12f')
#
#######################################
@testset "impresp, stepresp, freqresp, phaseresp" begin
    Hdelay = PolynomialRatio{:z}([0, 1], [1])
    W = range(0, stop=2π, length=100)
    for Tf in (PolynomialRatio, ZeroPoleGain, Biquad, SecondOrderSections)
        H = convert(Tf, Hdelay)
        @test freqresp(H, W) ≈ exp.(-W*im)
        @test phaseresp(H, W) ≈ -W
        @test grpdelay(H, W) ≈ fill(1.0, length(W))
        @test impresp(H, 10) == [0, 1, 0, 0, 0, 0, 0, 0, 0, 0]
        @test stepresp(H, 10) == [0, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    end

    matlab_resp = readdlm(joinpath(dirname(@__FILE__), "data", "responses-eg1.txt"),'\t')
    b0 = 0.05634
    b1 = [1,  1]
    b2 = [1, -1.0166, 1]
    a1 = [1, -0.683]
    a2 = [1, -1.4461, 0.7957]
    b = b0*conv(b1, b2)
    a = conv(a1, a2)
    df = PolynomialRatio(b, a)
    w = matlab_resp[:,1]

    #Impulse response
    impz_matlab = matlab_resp[:,2]
    @test impresp(df, 512) ≈ impz_matlab

    #Step response
    stepz_matlab = matlab_resp[:,3]
    @test stepresp(df, 512) ≈ stepz_matlab

    h_matlab = matlab_resp[:,4]
    @test abs.(freqresp(df, w)) ≈ h_matlab
    @test abs.(freqresp(SecondOrderSections(df), w)) ≈ h_matlab
    @test abs.(freqresp(ZeroPoleGain(df), w)) ≈ h_matlab

    phi_matlab = matlab_resp[:,5]
    @test phaseresp(df, w) ≈ phi_matlab

    # Test diffent versions of the functions
    H, w = freqresp(df)
    @test H == freqresp(df, w)
    @test minimum(w) ≤ 0                 # w should cover frequencies from 0 ...
    @test maximum(w) ≥ float(π)          # ... to π ...
    @test maximum(abs.(diff(w))) < 0.03  # and be sufficiently dense

    phi, w = phaseresp(df)
    @test phi == phaseresp(df, w)
    @test minimum(w) ≤ 0                 # w should cover frequencies from 0 ...
    @test maximum(w) ≥ float(π)          # ... to π ...
    @test maximum(abs.(diff(w))) < 0.03  # and be sufficiently dense

    @test cumsum(impresp(df)) ≈ stepresp(df)
end

#######################################
#
#  http://www.mathworks.com.au/help/signal/ref/freqs.html
#  Example 1: Frequency response from the transfer function
#
#  dlmwrite('freqs-eg1.txt',[w; mag; phasedeg]', 'delimiter', '\t', 'precision', '%.10f')
#
#######################################
@testset "freqresp(::FilterCoefficients{:s})" begin
    w = 10 .^ range(-1, stop=1, length=50)
    for (b, a, H) in (
            ([1, 0], [1], w*im), # s
            ([1], [1, 0], 1 ./ (w*im)), # 1/s
            ([0, 0, 1], [2, 10], 1 ./ (2*w*im .+ 10)), # 1/(2s + 10)
            ([1, 0, 1], [2, 10], (1 .- w.^2) ./ (2*w*im .+ 10)), # (s^2 + 1)/(2s + 10)
        )
        @test freqresp(PolynomialRatio{:s}(b, a), w) ≈ H
        @test freqresp(ZeroPoleGain(PolynomialRatio{:s}(b, a)), w) ≈ H
        if isone(a[1]) && lastindex(Polynomial(reverse(b))) <= lastindex(Polynomial(reverse(a)))
            @test freqresp(Biquad(PolynomialRatio{:s}(b, a)), w) ≈ H
        else
            # cannot represent this PolynomialRatio as a Biquad due to implied a0=1
            @test_broken freqresp(Biquad(PolynomialRatio{:s}(b, a)), w) ≈ H
        end
        if lastindex(Polynomial(reverse(b))) <= lastindex(Polynomial(reverse(a)))
            @test freqresp(SecondOrderSections(PolynomialRatio{:s}(b, a)), w) ≈ H
        else
            # cannot represent this PolynomialRatio as a Biquad due to implied a0=1
            @test_broken freqresp(SecondOrderSections(PolynomialRatio{:s}(b, a)), w) ≈ H
        end
    end

    # Julia
    a = [1.0, 0.4, 1.0]
    b = [0.2, 0.3, 1.0]
    w = 10 .^ range(-1, stop=1, length=50)

    h        = freqresp(PolynomialRatio{:s}(b, a), w)
    mag      = convert(Array{Float64}, abs.(h))
    phasedeg = (180/pi)*phaseresp(PolynomialRatio{:s}(b, a), w)

    # Matlab
    freqs_eg1_w_mag_phasedeg = readdlm(joinpath(dirname(@__FILE__), "data", "freqs-eg1.txt"),'\t')
    matlab_w        = freqs_eg1_w_mag_phasedeg[:,1]
    matlab_mag      = freqs_eg1_w_mag_phasedeg[:,2]
    matlab_phasedeg = freqs_eg1_w_mag_phasedeg[:,3]
    h_ref = [Polynomial(reverse(b))(s) / Polynomial(reverse(a))(s) for s ∈ im*w]

    # Test
    @test w ≈ matlab_w
    @test h ≈ h_ref
    @test mag ≈ matlab_mag
    @test phasedeg ≈ matlab_phasedeg

    @test h ≈ freqresp(ZeroPoleGain(PolynomialRatio{:s}(b, a)), w)
    @test h ≈ freqresp(SecondOrderSections(PolynomialRatio{:s}(b, a)), w)
end

# ######################################
#
#  Test grpdelay
#
#  Data from Matlab using b,a and a from above:
#  [gd, w] = grpdelay(b, a, 512)
#  all = [w gd]
#  dlmwrite('grpdelay_eg1.txt', all, 'delimiter', '\t', 'precision', '%.12f')
#
# ######################################
@testset "grpdelay" begin
    matlab_delay = readdlm(joinpath(dirname(@__FILE__), "data", "grpdelay_eg1.txt"),'\t')
    b0 = 0.05634
    b1 = [1,  1]
    b2 = [1, -1.0166, 1]
    a1 = [1, -0.683]
    a2 = [1, -1.4461, 0.7957]
    b = b0*conv(b1, b2)
    a = conv(a1, a2)
    df = PolynomialRatio(b, a)
    w = matlab_delay[:, 1]

    grpdelay_matlab = matlab_delay[:, 2]
    @test grpdelay(df, w) ≈ grpdelay_matlab

    # Test with IIR filters types I-IV
    tau, w = grpdelay(PolynomialRatio([1, 1, 1, 1, 1], [1]))
    @test tau ≈ fill(2.0, length(w))
    tau, w = grpdelay(PolynomialRatio([1, 1, 1, 1, 1, 1], [1]))
    @test tau ≈ fill(2.5, length(w))
    tau, w = grpdelay(PolynomialRatio([1, 0, -1], [1]))
    @test tau ≈ fill(1.0, length(w))
    tau, w = grpdelay(PolynomialRatio([1, -1], [1]))
    @test tau ≈ fill(0.5, length(w))
    @test minimum(w) ≤ 0                 # w should cover frequencies from 0 ...
    @test maximum(w) ≥ float(π)          # ... to π ...
    @test maximum(abs.(diff(w))) < 0.03  # and be sufficiently dense

    # group delay of analog filters
    w = matlab_delay[:, 1]
    @test grpdelay(PolynomialRatio{:s}([1, 2], [123]), w) ≈ 0.5 ./ (1 .+ w.^2/4)
    @test grpdelay(PolynomialRatio{:s}([123], [1, 2]), w) ≈ -0.5 ./ (1 .+ w.^2/4)
    @test grpdelay(PolynomialRatio{:s}([1, 2], [3, 9]), w) ≈ 0.5 ./ (1 .+ w.^2/4) - 1/3 ./ (1 .+ w.^2/9)
    tau, w = grpdelay(PolynomialRatio{:s}([1, 2], [123]))
    @test tau ≈ 0.5 ./ (1 .+ w.^2/4)
    @test minimum(w) < 2 < maximum(w) # position of zero covered
    tau, w = grpdelay(PolynomialRatio{:s}([123], [1, 2]))
    @test tau ≈ -0.5 ./ (1 .+ w.^2/4)
    @test minimum(w) < 2 < maximum(w) # position of pole covered
    tau, w = grpdelay(PolynomialRatio{:s}([1, 2], [3, 9]))
    @test tau ≈ 0.5 ./ (1 .+ w.^2/4) - 1/3 ./ (1 .+ w.^2/9)
    @test minimum(w) < 2 && 3 < maximum(w) # position of pole and zero covered
end


text:
# Auxiliary functions for natural cubic spline interpolation

# Stores parameters for fitted spline polynomials
mutable struct Spline{T}
    x::Vector{T} # (x, y) are the original points
    y::Vector{Float64}
    params::Vector{Float64} # stores parameters for polynomials. [ a1, b1, c1, d1, a2, b2, c2, d2 ...]

    function Spline(x::Vector{T}, y::Vector{Float64}, params::Vector{Float64}) where {T}
        polynms_count = length(x) - 1
        @assert length(params) == polynms_count * 4 "params length $(length(params)) does not conform to the expected number of polynomials ($(polynms_count))" # each polynomial has 4 parameters
        return new{T}(x, y, params)
    end
end

# Given a polinomial index, returns the indexer of the last parameter before the first parameter of the referenced polinomial
_base_param_index_(poly_index::Int) = (poly_index-1)*4

# Performs natural cubic spline interpolation
function splineint(s::Spline{T}, x_out::Number) where {T<:Number}
    local poly_index::Int = 1
    local base_idx::Int

    if x_out > s.x[end]
        # Extrapolation after last point
        poly_index = length(s.x) - 1
        base_idx = _base_param_index_(poly_index)
        return (x_out - s.x[end])*(s.params[base_idx + 2] + 2*s.params[base_idx + 3]*s.x[end] + 3*s.params[base_idx + 4]*(s.x[end]^2)) + s.y[end]
    elseif x_out < s.x[1]
        # Extrapolation before first point
        base_idx = 0
        return (x_out - s.x[1])*(s.params[base_idx + 2] + 2*s.params[base_idx + 3]*s.x[1] + 3*s.params[base_idx + 4]*(s.x[1]^2)) + s.y[1]
    else
        # Interplation
        while x_out > s.x[poly_index+1]
            poly_index += 1
        end

        base_idx = _base_param_index_(poly_index)

        #P1                       P2                    P3                     ...
        #1   2    3      4        5   6    7      8     9   10   11     12
        #a + bx + cx^2 + dx^3     a + bx + cx^2 + dx^3  a + bx + cx^2 + dx^3   ...
        return s.params[base_idx + 1] + s.params[base_idx + 2]*x_out +
            s.params[base_idx + 3]*(x_out^2) + s.params[base_idx + 4]*(x_out^3)
    end
end

# Performs natural cubic spline interpolation
function splineint(s::Spline{T1}, x_out::Vector{T2}) :: Vector{Float64} where {T1, T2}
    len = length(x_out)
    y_out = Vector{Float64}(undef, len)
    for i in 1:len
        y_out[i] = splineint(s, x_out[i])
    end
    return y_out
end

# Build a Spline object by fitting 3rd order polynomials around points (x_in, y_in)
function splinefit(x_in::Vector{T}, y_in::Vector{Float64}) :: Spline{T} where {T}
    #
    # TODO: optimize. See http://www.math.ntnu.no/emner/TMA4215/2008h/cubicsplines.pdf
    #
    points_count = length(x_in)
    @assert points_count == length(y_in) "x_in and y_in doesn't conform on sizes."
    matrix_n = 4*(points_count-1) # the main matrix is a square matrix matrix_n by matrix_n

    A = zeros(matrix_n, matrix_n)
    b = zeros(matrix_n)

    # Known values for polys

    # First Point
    A[1,1] = 1.0
    A[1,2] = x_in[1]
    A[1,3] = x_in[1]^2
    A[1,4] = x_in[1]^3
    b[1] = y_in[1]

    # Last Point
    base_idx = 4*(points_count - 2)
    A[2, base_idx + 1] = 1.0
    A[2, base_idx + 2] = x_in[points_count]
    A[2, base_idx + 3] = x_in[points_count]^2
    A[2, base_idx + 4] = x_in[points_count]^3
    b[2] = y_in[points_count]

    # Inner Points
    row = 3
    for i in 2:(points_count-1)
        # Connecting to left poly
        A[row, (i-2)*4 + 1] = 1.0
        A[row, (i-2)*4 + 2] = x_in[i]
        A[row, (i-2)*4 + 3] = x_in[i]^2
        A[row, (i-2)*4 + 4] = x_in[i]^3
        b[row] = y_in[i]

        row += 1

        # Connecting to right poly
        A[row, (i-1)*4 + 1] = 1.0
        A[row, (i-1)*4 + 2] = x_in[i]
        A[row, (i-1)*4 + 3] = x_in[i]^2
        A[row, (i-1)*4 + 4] = x_in[i]^3
        b[row] = y_in[i]

        row += 1

        # Conditions on first order derivatives
        A[row, (i-2)*4 + 2] = 1.0
        A[row, (i-2)*4 + 3] = 2.0*x_in[i]
        A[row, (i-2)*4 + 4] = 3.0*(x_in[i]^2)
        A[row, (i-1)*4 + 2] = -1.0
        A[row, (i-1)*4 + 3] = -2.0*x_in[i]
        A[row, (i-1)*4 + 4] = -3.0*(x_in[i]^2)

        row += 1

        # Conditions on second order derivatives
        A[row, (i-2)*4 + 3] = 2.0
        A[row, (i-2)*4 + 4] = 6.0 * x_in[i]
        A[row, (i-1)*4 + 3] = -2.0
        A[row, (i-1)*4 + 4] = -6.0 * x_in[i]

        row += 1
    end

    # Conditions for natural cubic spline
    A[row, 3] = 2.0
    A[row, 4] = 6.0 * x_in[1]

    row += 1

    A[row, (points_count-2)*4 + 3] = 2.0
    A[row, (points_count-2)*4 + 4] = 6.0 * x_in[points_count]

    return Spline(x_in, y_in, A \ b)
end


text:
function EG!(g::AbstractGraph, order=vertices(g))
    mg = construct_metagraph(g)
    origin_edges = edges(g)
    old_mv = 0
    for v in order
        mv = vertex_s2m(mg, v)
        sedges = getsaturate(mg, mv)
        add_edges!(mg, sedges)
        sedges::Vector{Vector{Int}} = [edge_m2s(mg, e[1], e[2]) for e in sedges]
        add_edges!(g, sedges)
        rem_vertex!(mg, mv)
        old_mv = mv
    end
    data = Dict("added_edges" => setdiff(edges(g), origin_edges))
    return data
end

function EG(g::AbstractGraph, order=vertices(g))
    cg = deepcopy(g)
    data = EG!(cg, order)
    return cg, data
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen