text:
#=
# Format specific functions for writing data output files.
=#

function output_values_raw(vars, file)
    # do one variable/component at a time
    if typeof(vars) <: Array
        for vi=1:length(vars) # each variable
            output_values_raw(vars[vi], file);
        end
    else
        # Here vars is just one variable
        N = size(vars.values, 2); # number of points
        for ci=1:size(vars.values,1) # each component for that variable
            file.write(vars.values[ci,:]); # write the values
        end
    end

end

function output_values_csv(vars, file)
    if config.solver_type == FV
        grid = fv_grid;
    else
        grid = grid_data;
    end
    # columns will be like x, y, z, u1,...
    if typeof(vars) <: Array
        N = size(vars[1].values, 2); # number of points (must be same for all vars)
        if vars[1].location == CELL
            x = fv_info.cellCenters;
        else
            x = grid.allnodes;
        end

        # header
        line = "x";
        if dim > 1
            line *= ", y";
        end
        if dim > 2
            line *= ", z";
        end
        for vi=1:length(vars) # each variable
            if size(vars[vi].values,1) > 1
                for ci=1:size(vars[vi].values,1)
                    line *= ", "*string(vars[vi].symbol)*"_"*string(ci);
                end
            else
                line *= ", "*string(vars[vi].symbol);
            end
        end
        println(file, line);

        # data
        for i=1:N
            # x, y, z
            line = string(x[1,i]);
            if dim > 1
                line *= ", "*string(x[2,i]);
            end
            if dim > 2
                line *= ", "*string(x[3,i]);
            end

            # u1, u2, ...
            for vi=1:length(vars) # each variable
                for ci=1:size(vars[vi].values,1)
                    line *= ", "*string(vars[vi].values[ci,i]);
                end
            end
            println(file, line);
        end

    else
        N = size(vars.values, 2); # number of points
        if vars.location == CELL
            x = fv_info.cellCenters;
        else
            x = grid.allnodes;
        end
        dim = size(x,1);

        # header
        line = "x";
        if dim > 1
            line *= ", y";
        end
        if dim > 2
            line *= ", z";
        end
        if size(vars.values,1) > 1
            for ci=1:size(vars.values,1)
                line *= ", "*string(vars.symbol)*"_"*string(ci);
            end
        else
            line *= ", "*string(vars.symbol);
        end
        println(file, line);

        # data
        for i=1:N
            # x, y, z
            line = string(x[1,i]);
            if dim > 1
                line *= ", "*string(x[2,i]);
            end
            if dim > 2
                line *= ", "*string(x[3,i]);
            end

            # u1, u2, ...
            for ci=1:size(vars.values,1)
                line *= ", "*string(vars.values[ci,i]);
            end
            println(file, line);
        end
    end
end

function output_values_vtk(vars, file, ascii)
    if config.solver_type == FV
        grid = fv_grid;
        nel = grid.nel_owned; # Only write owned elements
    else
        grid = grid_data;
        nel = size(grid.loc2glb,2); # all elements
    end
    # Use the WriteVTK package.
    # First build MeshCell array
    if config.dimension == 1
        np_to_type = Dict([(2,3)]); # line
    elseif config.dimension == 2
        np_to_type = Dict([(3,5),(4,9)]); # triangle, quad
    else
        np_to_type = Dict([(4,10),(8,12)]); # tet, hex
    end

    cells = Array{WriteVTK.MeshCell,1}(undef, nel);
    # nodes = fill(-2e20, 3, size(grid.allnodes, 2));
    for ei=1:nel
        nvertex = 0;
        for vi=1:length(grid.glbvertex[:,ei])
            if grid.glbvertex[vi,ei] > 0
                nvertex = nvertex+1;
            else
                break;
            end
        end
        if nvertex < 2
            printerr("Found element with too few vertices when writing VTK output: el="*string(ei)*", vertexmap="*string(grid.glbvertex[:,ei]));
            return;
        end
        cell_type = np_to_type[nvertex];
        cells[ei] = WriteVTK.MeshCell(VTKCellType(cell_type), grid.glbvertex[1:nvertex,ei]);
        # if config.dimension == 1
        #     nodes[1, grid.glbvertex[1:nvertex,ei]] = grid.allnodes[:, grid.glbvertex[1:nvertex,ei]];
        #     nodes[2:3, grid.glbvertex[1:nvertex,ei]] = zeros(2, nvertex);
        # elseif config.dimension == 2
        #     nodes[1:2, grid.glbvertex[1:nvertex,ei]] = grid.allnodes[:, grid.glbvertex[1:nvertex,ei]];
        #     nodes[3, grid.glbvertex[1:nvertex,ei]] = zeros(1, nvertex);
        # else
        #     nodes[:, grid.glbvertex[1:nvertex,ei]] = grid.allnodes[:, grid.glbvertex[1:nvertex,ei]];
        # end
    end

    # Let's try all nodes. If it causes problems, uncomment above to remove ghosts
    if config.dimension == 1
        nodes = zeros(3, size(grid.allnodes, 2));
        nodes[1, :] = grid.allnodes[1, :];
    elseif config.dimension == 2
        nodes = zeros(3, size(grid.allnodes, 2));
        nodes[1:2, :] = grid.allnodes[1:2, :];
    else
        nodes = grid.allnodes;
    end

    # Initialize the file
    # For multiple partitions write a pvtk file
    if config.num_procs > 1
        # # remove ghost nodes (they should have coordinates -2e20)
        # nextind = 1;
        # for ni=1:size(nodes,2)
        #     if !(nodes[1,ni] == -2e20)
        #         nodes[:,nextind] = nodes[:,ni];
        #         nextind += 1;
        #     end
        # end
        # nodes = nodes[:,1:nextind-1];

        vtkfile = pvtk_grid(file, nodes, cells, ascii=ascii, part=config.proc_rank+1, nparts=config.num_procs);
        vtkfile["partition"] = fill(config.partition_index, nel);
    else
        vtkfile = vtk_grid(file, nodes, cells, ascii=ascii);
    end

    # Add values
    if typeof(vars) <: Array
        for vi=1:length(vars)
            if length(vars[vi].symvar) > 1
                for ci=1:length(vars[vi].symvar)
                    compname = string(vars[vi].symbol) * "_" * string(ci);
                    vtkfile[compname] = vars[vi].values[ci,:];
                end
            else
                vtkfile[string(vars[vi].symbol)] = vars[vi].values[:];
            end
        end
    else
        if length(vars.symvar) > 1
            for ci=1:length(vars.symvar)
                compname = string(vars.symbol) * "_" * string(ci);
                vtkfile[compname] = vars.values[ci,:];
            end
        else
            vtkfile[string(vars.symbol)] = vars.values[:];
        end
    end

    return vtk_save(vtkfile);
end

########################################################################################################
# My attempt at writing a vtu file.
function output_values_myvtu(var, file, ascii)
    if config.solver_type == FV
        grid = fv_grid;
    else
        grid = grid_data;
    end

    if !(typeof(var) <: Array)
        vars = [var];
    else
        vars = var;
    end

    cells_only = config.solver_type == FV; # If FV, only use cell data.
    num_points = size(grid.allnodes,2);
    num_cells = grid.nel_owned;
    if ascii
        format = "ascii";
    else
        format = "appended";
    end
    # For appended data
    appended = zeros(UInt64, 0);
    chunk_offset = 1;
    byte_offset = 0;

    # cell types: ??
    nodes_per_element = size(grid.glbvertex,1);
    cell_types = zeros(Int8, num_cells);
    cells_per_line = 16;
    if config.dimension == 1
        for ci=1:num_cells
            cell_types[ci] = 3;
        end
    elseif config.dimension == 2
        for ci=1:num_cells
            if nodes_per_element == 3
                cell_types[ci] = 5;
            elseif nodes_per_element == 4
                cell_types[ci] = 9;
            end
        end
    else
        for ci=1:num_cells
            if nodes_per_element == 4
                cell_types[ci] = 10;
            elseif nodes_per_element == 8
                cell_types[ci] = 12;
            end
        end
    end

    content = "";

    # header
    content *= "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n";
    content *= "<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" byte_order=\"LittleEndian\""*(ascii ? "" : " header_type=\"UInt64\"")*">\n";
    content *= "  <UnstructuredGrid>\n";
    content *= "    <Piece NumberOfPoints=\""*string(num_points)*"\" NumberOfCells=\""*string(num_cells)*"\">\n";

    # Points
    content *= "      <Points>\n";
    if ascii
        content *= "        <DataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"3\" format=\"ascii\">\n";
    else
        content *= "        <DataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"3\" format=\"appended\" offset=\"0\">\n";
        (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, num_points*3 * 8, UInt64);
    end
    for ni=1:num_points
        scoper = chunk_offset + byte_offset;
        if ascii
            content *= "          ";
            content *= string(grid.allnodes[1, ni]) * " ";
            content *= (config.dimension > 1) ? (string(grid.allnodes[2, ni]) * " ") : "0 ";
            content *= (config.dimension > 2) ? (string(grid.allnodes[3, ni]) * " ") : "0 ";
            content *= "\n";
        else
            tmp1 = grid.allnodes[1, ni];
            tmp2 = (config.dimension > 1) ? grid.allnodes[2, ni] : 0.0;
            tmp3 = (config.dimension > 2) ? grid.allnodes[3, ni] : 0.0;
            (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, tmp1, Float64);
            (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, tmp2, Float64);
            (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, tmp3, Float64);
        end
    end
    content *= "        </DataArray>\n";
    content *= "      </Points>\n";

    # Cells
    content *= "      <Cells>\n";
    if ascii
        content *= "        <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n";
    else
        content *= "        <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\""*string((chunk_offset-1)*8 + byte_offset)*"\">\n";
        (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, num_cells*nodes_per_element * 4, UInt64);
    end
    for ci=1:num_cells
        scoper = chunk_offset + byte_offset;
        if ascii
            content *= "          ";
            for ni=1:nodes_per_element
                content *= string(grid.glbvertex[ni, ci]-1) * " ";
            end
            content *= "\n";
        else
            for ni=1:nodes_per_element
                (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, grid.glbvertex[ni, ci]-1, Int32);
            end
            #println("cell "*string(ci)*" nn "*string(nodes_per_element)*" co "*string(chunk_offset)*" bo "*string(byte_offset))
        end

    end
    content *= "        </DataArray>\n";

    if ascii
        content *= "        <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n";
    else
        #println("doing offsets co "*string(chunk_offset)*" bo "*string(byte_offset)*" -> "*string((chunk_offset-1)*8 + byte_offset))
        content *= "        <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\""*string((chunk_offset-1)*8 + byte_offset)*"\">\n";
        (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, num_cells * 4, UInt64);
    end
    tmp = 0;
    offset = 0;
    for ci=1:num_cells
        scoper = chunk_offset + byte_offset;
        if ascii
            if (tmp == 0) content *= "          "; end
            offset += nodes_per_element;
            content *= string(offset) * " ";
            tmp += 1;
            if (tmp == cells_per_line)
                content *= "\n";
                tmp = 0;
            elseif ci == num_cells
                content *= "\n";
            end
        else
            offset += nodes_per_element;
            (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, offset, Int32);
            #println("cell "*string(ci)*" off "*string(offset)*" co "*string(chunk_offset)*" bo "*string(byte_offset))
        end
    end
    content *= "        </DataArray>\n";

    if ascii
        content *= "        <DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n";
    else
        #println("doing types co "*string(chunk_offset)*" bo "*string(byte_offset)*" -> "*string((chunk_offset-1)*8 + byte_offset))
        content *= "        <DataArray type=\"UInt8\" Name=\"types\" format=\"appended\" offset=\""*string((chunk_offset-1)*8 + byte_offset)*"\">\n";
        (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, num_cells, UInt64);
    end
    tmp = 0;
    for ci=1:num_cells
        scoper = chunk_offset + byte_offset;
        if ascii
            if (tmp == 0) content *= "          "; end
            content *= string(cell_types[ci]) * " ";
            tmp += 1;
            if (tmp == cells_per_line)
                content *= "\n";
                tmp = 0;
            elseif ci == num_cells
                content *= "\n";
            end
        else
            (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, cell_types[ci], Int8);
        end
    end
    content *= "        </DataArray>\n";
    content *= "      </Cells>\n";

    # Point data
    content *= "      <PointData>\n";
    if !cells_only
        for vi=1:length(vars)
            scoper = chunk_offset + byte_offset;
            comps = size(vars[vi].values,1);
            if ascii
                content *= "        <DataArray type=\"Float64\" Name=\""*string(vars[vi].symbol)*"\""*
                                " NumberOfComponents=\""*string(comps)*"\"" * " format=\"ascii\">\n";
            else
                content *= "        <DataArray type=\"Float64\" Name=\""*string(vars[vi].symbol)*"\""*
                                " NumberOfComponents=\""*string(comps)*"\"" * " format=\"appended\" offset=\""*string((chunk_offset-1)*8 + byte_offset)*"\">\n";
                (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, num_points*comps * 8, UInt64);
            end
            tmp = 0;
            for ni=1:num_points
                if ascii
                    if (tmp == 0) content *= "          "; end
                    if comps > 1
                        for d=1:comps
                            content *= string(vars[vi].values[d, ni]) * " ";
                        end
                        content *= "\n";
                    else
                        content *= string(vars[vi].values[1, ni]) * " ";
                        tmp += 1;
                        if (tmp >= cells_per_line/2)
                            content *= "\n";
                            tmp = 0;
                        elseif ni==num_points
                            content *= "\n";
                        end
                    end
                else
                    for d=1:comps
                        scoper = chunk_offset + byte_offset;
                        (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, vars[vi].values[d, ni], Float64);
                    end
                end
            end
            content *= "        </DataArray>\n";
        end
    end
    content *= "      </PointData>\n";

    # Cell data
    content *= "      <CellData>\n";
    if cells_only
        for vi=1:length(vars)
            scoper = chunk_offset + byte_offset;
            comps = size(vars[vi].values,1);
            if ascii
                content *= "        <DataArray type=\"Float64\" Name=\""*string(vars[vi].symbol)*
                                " NumberOfComponents=\""*string(comps)*"\"" * " format=\"ascii\">\n";
            else
                content *= "        <DataArray type=\"Float64\" Name=\""*string(vars[vi].symbol)*
                                " NumberOfComponents=\""*string(comps)*"\"" * " format=\"appended\" offset=\""*string((chunk_offset-1)*8 + byte_offset)*"\">\n";
                (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, num_cells*comps * 8, UInt64);
            end
            tmp = 0;
            for ni=1:num_cells
                if ascii
                    if (tmp == 0) content *= "          "; end
                    if comps > 1
                        for d=1:comps
                            content *= string(vars[vi].values[d, ni]) * " ";
                        end
                        content *= "\n";
                    else
                        content *= string(vars[vi].values[1, ni]) * " ";
                        tmp += 1;
                        if (tmp >= cells_per_line/2)
                            content *= "\n";
                            tmp = 0;
                        elseif ni==num_cells
                            content *= "\n";
                        end
                    end
                else
                    for d=1:comps
                        scoper = chunk_offset + byte_offset;
                        (chunk_offset, byte_offset) = add_to_appended_data!(appended, chunk_offset, byte_offset, vars[vi].values[d, ni], Float64);
                    end
                end
            end
            content *= "        </DataArray>\n";
        end
    end
    content *= "      </CellData>\n";

    content *= "    </Piece>\n";
    content *= "  </UnstructuredGrid>\n";

    # appended data if needed
    if !ascii
        content *= "  <AppendedData encoding=\"raw\">\n";
        content *= "_";
        print(file, content);
        for i=1:length(appended)
            write(file, htol(appended[i]));
        end

        content = "\n  </AppendedData>\n";
    end

    # ending
    content *= "</VTKFile>\n";

    println(file, content);
end

# adds one piece of data to an appended data set.
function add_to_appended_data!(data, chunk_offset, byte_offset, newdata, type)
    bytes = UInt64(0);
    tmp = UInt64(0);
    bc = sizeof(type);
    if type == Float64
        intdata = reinterpret(UInt64, newdata)[1];
    elseif type == Float32
        intdata = reinterpret(UInt32, Float32(newdata))[1];
    else
        intdata = newdata;
    end

    # insert the bytes of newdata in little endian
    for i=1:bc
        tmp = (intdata >>> (8*(i-1))) & 0xFF;
        tmp = tmp << (8*(i-1));
        bytes = bytes | tmp;
        tmp = UInt64(0);
    end

    if byte_offset > 0
        lowbc = min(8-byte_offset, bc);
        highbc = bc-lowbc;
        lowtmp = data[chunk_offset];
        hightmp = UInt64(0);

        lowtmp = lowtmp | (bytes << (8*(byte_offset)));

        hightmp = bytes >>> (8*(8-byte_offset));

        data[chunk_offset] = lowtmp;
        if highbc > 0
            push!(data, hightmp);
            return (chunk_offset + 1, highbc);
        else
            byte_offset = byte_offset + lowbc;
            if byte_offset == 8
                chunk_offset += 1;
                byte_offset = 0;
            end
            return (chunk_offset, byte_offset);
        end

    else
        push!(data, bytes);
        if bc == 8
            chunk_offset += 1;
            byte_offset = 0;
        else
            byte_offset = bc;
        end

        return (chunk_offset, byte_offset);
    end
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
### Utility Functions

function build_component_funcs(symex)
  funcs = Vector{Expr}(undef,0) # Get all of the functions for symbolic computation
  for (i,arg) in enumerate(symex.args)
    if i%2 == 0
      ex = arg.args[2]
      if (typeof(ex) <: Symbol) || (typeof(ex) <: Number)
        push!(funcs,:($ex*1))
      else # It's an expression, just push
        fix_ex = copy(ex)
        flip_mult!(fix_ex)
        push!(funcs,fix_ex)
      end
    end
  end
  funcs
end

function modelingtoolkitize_expr(ex::Expr,vars,curmod)
    names = [tosymbol(x) for x in vars]
    ex.head === :if && (ex = Expr(:call, ifelse, ex.args...))
    ex.head === :call || throw(ArgumentError("internal representation does not support non-call Expr"))
    op = ex.args[1] ∈ names ? vars[findfirst(x->ex.args[1] == tosymbol(x),vars)] : getproperty(curmod,ex.args[1]) # HACK
    return op((modelingtoolkitize_expr(x,vars,curmod) for x in ex.args[2:end])...)
end

function modelingtoolkitize_expr(ex::Sym,vars,curmod)
    ex
end

function modelingtoolkitize_expr(ex::Symbol,vars,curmod)
    names = tosymbol.(vars)
    idx = findfirst(x->ex == x, names)
    op = idx !== nothing ? vars[idx] : getproperty(curmod,ex)
    op
end

function modelingtoolkitize_expr(ex::Number,vars,curmod)
    ex
end

ext: jl


text:
module FCSFiles

using FileIO

include("type.jl")
include("utils.jl")
include("parse.jl")

export FlowSample

function load(f::File{format"FCS"})
    open(f) do io
        offsets = parse_header(io)

        text_mappings = parse_text(io, offsets[1], offsets[2])
        verify_text(text_mappings)

        return parse_data(io, offsets[3], offsets[4], text_mappings)
    end
end

end # module


text:
using OutlierDetection
using Documenter
using DocumenterMarkdown

DocMeta.setdocmeta!(OutlierDetection, :DocTestSetup, :(using OutlierDetection); recursive=true)
makedocs(;
    doctest = VERSION == v"1.6",
    sitename="OutlierDetection.jl",
    authors="David Muhr <muhrdavid@gmail.com> and contributors",
    repo="https://github.com/davnn/OutlierDetection.jl/blob/{commit}{path}#L{line}",
    format = Markdown(),
    clean=true
)


text:
for N in [Float64, Rational{Int}, Float32]
    # random empty set
    rand(EmptySet)

    E = EmptySet{N}(2)
    B = BallInf(ones(N, 2), N(1))

    # testing that the empty set is an absorbing element for the cartesian product
    @test B * E isa EmptySet && E * B isa EmptySet
    # testing the mathematical alias ×
    @test B × E isa EmptySet && E × B isa EmptySet

    cpa = CartesianProductArray([B, N(2) * B, N(3) * B])
    @test cpa * E isa EmptySet && E * cpa isa EmptySet
    @test cpa × E isa EmptySet && E × cpa isa EmptySet

    # testing cp of empty set with itself
    @test E * E == E

    # testing that the empty set is an absorbing element for the Minkowski sum
    @test B + E isa EmptySet && E + B isa EmptySet
    # testing the mathematical alias ⊕
    @test B ⊕ E isa EmptySet && E ⊕ B isa EmptySet

    msa = MinkowskiSumArray([B, N(2) * B, N(3) * B])
    @test msa + E isa EmptySet && E + msa isa EmptySet
    @test msa ⊕ E isa EmptySet && E ⊕ msa isa EmptySet

    # testing M-sum of empty set with itself
    @test E + E == E

    # testing that the emptyset is neutral for the convex hull
    @test CH(B, E) == B
    @test CH(E, B) == B

    # test convex hull of empty set with itself
    @test CH(E, E) == E

    # dim
    @test dim(E) == 2

    # support vector
    @test_throws ErrorException σ(N[0], E)

    # boundedness
    @test isbounded(E)

    # isuniversal
    res, w = isuniversal(E, true)
    @test !isuniversal(E) && !res && w ∉ E

    # membership
    @test N[0] ∉ E
    @test N[0, 0] ∉ E

    # subset
    @test E ⊆ B && ⊆(E, B, true)[1]
    subset, point = ⊆(B, E, true)
    @test B ⊈ E && !subset && point ∈ B
    @test E ⊆ E && ⊆(E, E, true)[1]

    # emptiness check
    @test isempty(E)

    # an_element/norm/radius/diameter functions
    @test_throws ErrorException an_element(E)
    @test_throws ErrorException norm(E)
    @test_throws ErrorException radius(E)
    @test_throws ErrorException diameter(E)

    # linear map of an empty set
    linear_map(ones(N, 2, 2), E) == E

    # translation
    @test translate(E, N[1, 2]) == E
end

# default Float64 constructor
@test EmptySet(2) == ∅(2) == EmptySet{Float64}(2)


text:
"""
    spar = splineparams(T, D1, Nseg, pcof)

Constructor for struct splineparams, which sets up the parameters for a regular B-spline function
(without carrier waves).

# Arguments
- `T:: Float64`: Duration of spline function
- `D1:: Int64`: Number of basis functions in each spline
- `Nseg:: Int64`:  Number of splines (real, imaginary, different ctrl func)
- `pcof:: Array{Float64, 1}`: Coefficient vector. Must have D1*Nseg elements

# External links
* [Spline Wavelet](https://en.wikipedia.org/wiki/Spline_wavelet#Quadratic_B-spline) on Wikipedia.
"""
struct splineparams
    T::Float64
    D1::Int64 # Number of coefficients per spline
    tcenter::Array{Float64,1}
    dtknot::Float64
    pcof::Array{Float64,1} # pcof should have D1*Nseg elements
    Nseg::Int64 # Number of segments (real, imaginary, different ctrl func)
    Ncoeff:: Int64 # Total number of coefficients

# new, simplified constructor
    function splineparams(T, D1, Nseg, pcof)
        dtknot = T/(D1 -2)
        tcenter = dtknot.*(collect(1:D1) .- 1.5)
        new(T, D1, tcenter, dtknot, pcof, Nseg, Nseg*D1)
    end

end

# bspline2: Evaluate quadratic bspline function
"""
    f = bspline2(t, splineparam, splinefunc)

Evaluate a B-spline function. See also the `splineparams` constructor.

# Arguments
- `t::Float64`: Evaluate spline at parameter t ∈ [0, param.T]
- `param::splineparams`: Parameters for the spline
- `splinefunc::Int64`: Spline function index ∈ [0, param.Nseg-1]
"""
@inline function bspline2(t::Float64, param::splineparams, splinefunc::Int64)
  f = 0.0

  dtknot = param.dtknot
  width = 3*dtknot

  offset = splinefunc*param.D1

  k = max.(3, ceil.(Int64,t./dtknot + 2)) # Unsure if this line does what it is supposed to
  k = min.(k, param.D1)

  # 1st segment of nurb k
  tc = param.tcenter[k]
  tau = (t .- tc)./width
  f = f + param.pcof[offset+k] * (9/8 .+ 4.5*tau + 4.5*tau^2) # test to remove square for extra speed

  # 2nd segment of nurb k-1
  tc = param.tcenter[k-1]
  tau = (t - tc)./width
  f = f .+ param.pcof[offset+k.-1] .* (0.75 - 9 *tau^2)

  # 3rd segment of nurb k-2
  tc = param.tcenter[k.-2]
  tau = (t .- tc)./width
  f = f + param.pcof[offset+k.-2] * (9/8 - 4.5*tau + 4.5*tau.^2)
end

"""
    g = gradbspline2(t, param, splinefunc)

Evaluate the gradient of a spline function with respect to all coefficient.
NOTE: the index of the spline functions are 0-based. For a set of
coupled controls, mod(`splinefunc`,2)=0 corresponds to ∇ p_j(t) and mod(`splinefunc`,2) = 1
corresponds to ∇ q_j(t), where j = div(splinefunc,2).

# Arguments
- `t::Float64`: Evaluate spline at parameter t ∈ [0, param.T]
- `param::splineparams`: Spline parameter object
- `splinefunc::Int64`: Spline function index ∈ [0, param.Nseg-1]
"""
@inline function gradbspline2(t::Float64,param::splineparams, splinefunc::Int64)

# NOTE: param.Nseg used to be '2'
  g = zeros(param.Nseg*param.D1) # real and imag parts for both f and g

  dtknot = param.dtknot
  width = 3*dtknot

  offset = splinefunc*param.D1

  k = max.(3, ceil.(Int64,t./dtknot .+ 2)) # t_knot(k-1) < t <= t_knot(k), but t=0 needs to give k=3
  k = min.(k, param.D1) # protect agains roundoff that sometimes makes t/dt > N_nurbs-2

  #1st segment of nurb k
  tc = param.tcenter[k]
  tau = (t .- tc)./width
  g[offset .+ k] = (9/8 .+ 4.5.*tau .+ 4.5.*tau.^2);

  #2nd segment of nurb k-1
  tc = param.tcenter[k.-1]
  tau = (t .- tc)./width
  g[offset .+ k.-1] = (0.75 .- 9 .*tau.^2)

  # 3rd segment og nurb k-2
  tc = param.tcenter[k.-2]
  tau = (t .- tc)./width
  g[offset .+ k.-2] = (9/8 .- 4.5.*tau .+ 4.5.*tau.^2);
  return g
end

"""
    bcpar = bcparams(T, D1, Ncoupled, Nunc, omega, pcof)

General constructor of `struct bcparams` for setting up B-splines with carrier waves.

    bcpar = bcparams(T, D1, omega, pcof)

Simplified constructor for the case when there are no uncoupled controls and `Ncoupled = size(omega,1)`.

# Arguments
- `T:: Float64`: Duration of spline function
- `D1:: Int64`: Number of basis functions in each segment
- `Ncoupled::Int64`: Number of coupled controls in the simulation
- `Nunc::Int64`: Number of uncoupled controls in the simulation
- `omega::Array{Float64,2}`: Carrier wave frequencies
- `pcof:: Array{Float64, 1}`: Coefficient vector. Must have D1*Nseg elements

# First dimensions of the `omega` array:
- Without uncoupled controls, `Nunc=0` and `size(omega,1) = Ncoupled`.
- With uncoupled controls, `Nunc > 0` and `size(omega,1) = Ncoupled + Nunc`.

# Second dimension of the `omega` array:
- `size(omega, 2) = Nfreq`

# Ordering of the `pcof` array:
First consider the case without uncoupled control functions, `Nunc = 0`:
Then the `pcof` array then has `2*Ncoupled*Nfreq*D1` elements.
Each `ctrl ∈ [1,Ncoupled]` and `freq ∈ [1,Nfreq]` corresponds to `D1` elements in
the `pcof` vector. For the case `Ncoupled = 2` and `Nfreq = 2`, the elements are ordered according to

| ctrl    | freq    | α_1  | α_2 |
| ------- | -------- | -------- | -------|
| 1      | 1        | 1:D1             | D1+1:2 D1 |
| 1      | 2        | 2 D1+1: 3 D1 | 3 D1+1:4 D1 |
| 2      | 1        | 4 D1+1: 5 D1 | 5 D1+1:6 D1 |
| 2      | 2        | 6 D1+1: 7 D1 | 7 D1+1: 8D1 |

If there are uncoupled controls, `Nunc > 0`, the `pcof` array should have `(2*Ncoupled + Nunc)*Nfreq*D1` elements.
The last `Nunc*Nfreq*D1` elements correspond to the uncoupled control functions and are ordered in a corresponding way.

# External links
* [Spline Wavelet](https://en.wikipedia.org/wiki/Spline_wavelet#Quadratic_B-spline) on Wikipedia.
"""
struct bcparams
    T ::Float64
    D1::Int64 # number of B-spline coefficients per control function
    om::Array{Float64,2} #Carrier wave frequencies [rad/s], size Nfreq
    tcenter::Array{Float64,1}
    dtknot::Float64
    pcof::Array{Float64,1} # coefficients for all 2*Ncoupled splines, size Ncoupled*D1*Nfreq*2 (*2 because of sin/cos)
    Nfreq::Int64 # Number of frequencies
    Ncoeff:: Int64 # Total number of coefficients
    Ncoupled::Int64 # Number of B-splines functions for the coupled ctrl Hamiltonians
    Nunc::Int64 # Number of B-spline functions  for the UNcoupled ctrl Hamiltonians

    # New constructor to allow defining number of symmetric Hamiltonian terms
    function bcparams(T::Float64, D1::Int64, Ncoupled::Int64, Nunc::Int64, omega::Array{Float64,2}, pcof::Array{Float64,1})
        dtknot = T/(D1 -2)
        tcenter = dtknot.*(collect(1:D1) .- 1.5)
        Nfreq = size(omega,2)
        nCoeff = Nfreq*D1*2*(Ncoupled + Nunc)
        if nCoeff != length(pcof)
            println("nCoeff = ", nCoeff, " Nfreq = ", Nfreq, " D1 = ", D1, " Ncoupled = ", Ncoupled, " Nunc = ", Nunc, " len(pcof) = ", length(pcof))
            throw(DimensionMismatch("Inconsistent number of coefficients and size of parameter vector (nCoeff ≠ length(pcof)."))
        end
        new(T, D1, omega, tcenter, dtknot, pcof, Nfreq, nCoeff, Ncoupled, Nunc)
    end

end

# simplified constructor (assumes no uncoupled terms)
function bcparams(T::Float64, D1::Int64, omega::Array{Float64,2}, pcof::Array{Float64,1})
  dtknot = T/(D1 -2)
  tcenter = dtknot.*(collect(1:D1) .- 1.5)
  Ncoupled = size(omega,1) # should check that Ncoupled >=1
  Nfreq = size(omega,2)
  Nunc = 0
  nCoeff = Nfreq*D1*2*Ncoupled
  if nCoeff != length(pcof)
    throw(DimensionMismatch("Inconsistent number of coefficients and size of parameter vector (nCoeff ≠ length(pcof)."))
  end
  bcparams(T, D1, Ncoupled, Nunc, omega, pcof)
end

"""
    f = bcarrier2(t, params, func)

Evaluate a B-spline function with carrier waves. See also the `bcparams` constructor.

# Arguments
- `t::Float64`: Evaluate spline at parameter t ∈ [0, param.T]
- `param::params`: Parameters for the spline
- `func::Int64`: Spline function index ∈ [0, param.Nseg-1]
"""
@inline function bcarrier2(t::Float64, params::bcparams, func::Int64)
    # for a single oscillator, func=0 corresponds to p(t) and func=1 to q(t)
    # in general, 0 <= func < 2*Ncoupled + Nunc

    # compute basic offset: func 0 and 1 use the same spline coefficients, but combined in a different way
    osc = div(func, 2) # osc is base 0; 0<= osc < Ncoupled
    q_func = func % 2 # q_func = 0 for p and q_func=1 for q

    f = 0.0 # initialize

    dtknot = params.dtknot
    width = 3*dtknot

    k = max.(3, ceil.(Int64,t./dtknot + 2)) # pick out the index of the last basis function corresponding to t
    k = min.(k, params.D1) #  Make sure we don't access outside the array

    if func < 2*(params.Ncoupled + params.Nunc)
        # Coupled and uncoupled controls
        @fastmath @inbounds @simd for freq in 1:params.Nfreq
            fbs1 = 0.0 # initialize
            fbs2 = 0.0 # initialize
            # offset in parameter array (osc = 0,1,2,...
            # Vary freq first, then osc
            offset1 = 2*osc*params.Nfreq*params.D1 + (freq-1)*2*params.D1
            offset2 = 2*osc*params.Nfreq*params.D1 + (freq-1)*2*params.D1 + params.D1

            # 1st segment of nurb k
            tc = params.tcenter[k]
            tau = (t .- tc)./width
            fbs1 += params.pcof[offset1+k] * (9/8 .+ 4.5*tau + 4.5*tau^2)
            fbs2 += params.pcof[offset2+k] * (9/8 .+ 4.5*tau + 4.5*tau^2)

            # 2nd segment of nurb k-1
            tc = params.tcenter[k-1]
            tau = (t - tc)./width
            fbs1 += params.pcof[offset1+k.-1] .* (0.75 - 9 *tau^2)
            fbs2 += params.pcof[offset2+k.-1] .* (0.75 - 9 *tau^2)

            # 3rd segment of nurb k-2
            tc = params.tcenter[k.-2]
            tau = (t .- tc)./width
            fbs1 += params.pcof[offset1+k-2] * (9/8 - 4.5*tau + 4.5*tau.^2)
            fbs2 += params.pcof[offset2+k-2] * (9/8 - 4.5*tau + 4.5*tau.^2)

            #    end # for carrier phase
            # p(t)
            if q_func==1
                f += fbs1 * sin(params.om[osc+1,freq]*t) + fbs2 * cos(params.om[osc+1,freq]*t) # q-func
            else
                f += fbs1 * cos(params.om[osc+1,freq]*t) - fbs2 * sin(params.om[osc+1,freq]*t) # p-func
            end
        end # for freq
    end # if
    # else

    #   # uncoupled control
    #   @fastmath @inbounds @simd for freq in 1:params.Nfreq
    #     fbs = 0.0 # initialize

    #     # offset to uncoupled parameters
    #     # offset = (2*params.Ncoupled)*params.D1*params.Nfreq + (fun - 2*params.Ncoupled)*params.D1*params.Nfreq  +  (freq-1)*params.D1
    #     offset = func*params.D1*params.Nfreq  +  (freq-1)*params.D1

    #     # 1st segment of nurb k
    #     if k <= params.D1  # could do the if statement outside the for-loop instead
    #       tc = params.tcenter[k]
    #       tau = (t .- tc)./width
    #       fbs += params.pcof[offset+k] * (9/8 .+ 4.5*tau + 4.5*tau^2)
    #     end

    #     # 2nd segment of nurb k-1
    #     if k >= 2 && k <= params.D1 + 1
    #       tc = params.tcenter[k-1]
    #       tau = (t - tc)./width
    #       fbs += params.pcof[offset+k.-1] .* (0.75 - 9 *tau^2)
    #     end

    #     # 3rd segment of nurb k-2
    #     if k >= 3 && k <= params.D1 + 2
    #       tc = params.tcenter[k.-2]
    #       tau = (t .- tc)./width
    #       fbs += params.pcof[offset+k-2] * (9/8 - 4.5*tau + 4.5*tau.^2)
    #     end

    #     # FMG: For now this would assume that every oscillator would have
    #     # an uncoupled term if uncoupled terms exist in the problem. Perhaps
    #     # update parameter struct to track terms if they are only present on some
    #     # oscillators?
    #     ind = func - 2*params.Ncoupled
    #     f += fbs * cos(params.om[ind+1,freq]*t) # spl is 0-based
    #   end # for

    return f
end

"""
    gradbcarrier2!(t, params, func, g) -> g
Evaluate the gradient of a control function with respect to all coefficient.

NOTE: the index of the control functions is 0-based. For a set of
coupled controls, mod(`func`,2)=0 corresponds to ∇ p_j(t) and mod(`func`,2) = 1
corresponds to ∇ q_j(t), where j = div(func,2).

# Arguments
- `t::Float64`: Evaluate spline at parameter t ∈ [0, param.T]
- `params::bcparams`: Parameters for the spline
- `func::Int64`: Control function index ∈ [0, param.Nseg-1]
- `g::Array{Float64,1}`: Preallocated array to store calculated gradient
"""
function gradbcarrier2!(t::Float64, params::bcparams, func::Int64, g::Array{Float64,1})

    # compute basic offset: func 0 and 1 use the same spline coefficients, but combined in a different way
    osc = div(func, 2)
    q_func = func % 2 # q_func = 0 for p and q_func=1 for q

    # allocate array for returning the results
    # g = zeros(length(params.pcof)) # cos and sin parts

    g .= 0.0

    dtknot = params.dtknot
    width = 3*dtknot

    k = max.(3, ceil.(Int64,t./dtknot .+ 2)) # t_knot(k-1) < t <= t_knot(k), but t=0 needs to give k=3
    k = min.(k, params.D1) # protect agains roundoff that sometimes makes t/dt > N_nurbs-2

    if func < 2*(params.Ncoupled + params.Nunc)
        @fastmath @inbounds @simd for freq in 1:params.Nfreq

            # offset in parameter array (osc = 0,1,2,...
            # Vary freq first, then osc
            offset1 = 2*osc*params.Nfreq*params.D1 + (freq-1)*2*params.D1
            offset2 = 2*osc*params.Nfreq*params.D1 + (freq-1)*2*params.D1 + params.D1

            #1st segment of nurb k
            tc = params.tcenter[k]
            tau = (t .- tc)./width
            bk = (9/8 .+ 4.5.*tau .+ 4.5.*tau.^2)
            if q_func==1
                g[offset1 .+ k] = bk * sin(params.om[osc+1,freq]*t)
                g[offset2 .+ k] = bk * cos(params.om[osc+1,freq]*t)
            else # p-func
                g[offset1 .+ k] = bk * cos(params.om[osc+1,freq]*t)
                g[offset2 .+ k] = -bk * sin(params.om[osc+1,freq]*t)
            end

            #2nd segment of nurb k-1
            tc = params.tcenter[k.-1]
            tau = (t .- tc)./width
            bk = (0.75 .- 9 .*tau.^2)
            if q_func==1
                g[offset1 .+ (k-1)] = bk * sin(params.om[osc+1,freq]*t)
                g[offset2 .+ (k-1)] = bk * cos(params.om[osc+1,freq]*t)
            else # p-func
                g[offset1 .+ (k-1)] = bk * cos(params.om[osc+1,freq]*t)
                g[offset2 .+ (k-1)] = -bk * sin(params.om[osc+1,freq]*t)
            end

            # 3rd segment og nurb k-2
            tc = params.tcenter[k.-2]
            tau = (t .- tc)./width
            bk = (9/8 .- 4.5.*tau .+ 4.5.*tau.^2)
            if q_func==1
                g[offset1 .+ (k-2)] = bk * sin(params.om[osc+1,freq]*t)
                g[offset2 .+ (k-2)] = bk * cos(params.om[osc+1,freq]*t)
            else # p-func
                g[offset1 .+ (k-2)] = bk * cos(params.om[osc+1,freq]*t)
                g[offset2 .+ (k-2)] = -bk * sin(params.om[osc+1,freq]*t)
            end

        end #for freq
    # else
    #     # uncoupled control case
    #   @fastmath @inbounds @simd for freq in 1:params.Nfreq

    #     # offset
    #     offset = func*params.D1*params.Nfreq  +  (freq-1)*params.D1
    #     ind = func - 2*params.Ncoupled

    #     #1st segment of nurb k
    #     if k <= params.D1  # could do the if statement outside the for-loop instead
    #       tc = params.tcenter[k]
    #       tau = (t .- tc)./width
    #       g[offset .+ k] = (9/8 .+ 4.5.*tau .+ 4.5.*tau.^2)*cos(params.om[ind+1,freq]*t)
    #     end

    #     #2nd segment of nurb k-1
    #     if k >= 2 && k <= params.D1 + 1
    #       tc = params.tcenter[k.-1]
    #       tau = (t .- tc)./width
    #       g[offset .+ k.-1] = (0.75 .- 9 .*tau.^2)*cos(params.om[ind+1,freq]*t)
    #     end

    #     # 3rd segment og nurb k-2
    #     if k >= 3 && k <= params.D1 + 2
    #       tc = params.tcenter[k.-2]
    #       tau = (t .- tc)./width
    #       g[offset .+ k.-2] = (9/8 .- 4.5.*tau .+ 4.5.*tau.^2)*cos(params.om[ind+1,freq]*t)
    #     end

    #   end #for freq

       end #if
 end


text:
"""
    fk425(r1950, d1950, dr1950, dd1950, p1950, v1950)

Convert B1950.0 FK4 star catalog data to J2000.0 FK5.

This function converts a star's catalog data from the old FK4
Bessel-Newcomb) system to the later IAU 1976 FK5 (Fricke) system.

### Given (all B1950.0, FK4) ###

- `r1950,d1950`: B1950.0 RA,Dec (rad)
- `dr1950,dd1950`: B1950.0 proper motions (rad/trop.yr)
- `p1950`: Parallax (arcsec)
- `v1950`: Radial velocity (km/s, +ve = moving away)

### Returned (all J2000.0, FK5) ###

- `r2000,d2000`: J2000.0 RA,Dec (rad)
- `dr2000,dd2000`: J2000.0 proper motions (rad/Jul.yr)
- `p2000`: Parallax (arcsec)
- `v2000`: Radial velocity (km/s, +ve = moving away)

### Notes ###

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

2. The conversion is somewhat complicated, for several reasons:

   - Change of standard epoch from B1950.0 to J2000.0.

   - An intermediate transition date of 1984 January 1.0 TT.

   - A change of precession model.

   - Change of time unit for proper motion (tropical to Julian).

   - FK4 positions include the E-terms of aberration, to simplify
     the hand computation of annual aberration.  FK5 positions
     assume a rigorous aberration computation based on the Earth's
     barycentric velocity.

   - The E-terms also affect proper motions, and in particular cause
     objects at large distances to exhibit fictitious proper
     motions.

   The algorithm is based on Smith et al. (1989) and Yallop et al.
   (1989), which presented a matrix method due to Standish (1982) as
   developed by Aoki et al. (1983), using Kinoshita's development of
   Andoyer's post-Newcomb precession.  The numerical constants from
   Seidelmann (1992) are used canonically.

3. Conversion from B1950.0 FK4 to J2000.0 FK5 only is provided for.
   Conversions for different epochs and equinoxes would require
   additional treatment for precession, proper motion and E-terms.

4. In the FK4 catalog the proper motions of stars within 10 degrees
   of the poles do not embody differential E-terms effects and
   should, strictly speaking, be handled in a different manner from
   stars outside these regions.  However, given the general lack of
   homogeneity of the star data available for routine astrometry,
   the difficulties of handling positions that may have been
   determined from astrometric fields spanning the polar and non-
   polar regions, the likelihood that the differential E-terms
   effect was not taken into account when allowing for proper motion
   in past astrometry, and the undesirability of a discontinuity in
   the algorithm, the decision has been made in this ERFA algorithm
   to include the effects of differential E-terms on the proper
   motions for all stars, whether polar or not.  At epoch J2000.0,
   and measuring "on the sky" rather than in terms of RA change, the
   errors resulting from this simplification are less than
   1 milliarcsecond in position and 1 milliarcsecond per century in
   proper motion.

### Called ###

- [`anp`](@ref): normalize angle into range 0 to 2pi
- [`pv2s`](@ref): pv-vector to spherical coordinates
- [`pdp`](@ref): scalar product of two p-vectors
- [`pvmpv`](@ref): pv-vector minus pv_vector
- [`pvppv`](@ref): pv-vector plus pv_vector
- [`s2pv`](@ref): spherical coordinates to pv-vector
- [`sxp`](@ref): multiply p-vector by scalar

### References ###

- Aoki, S. et al., 1983, "Conversion matrix of epoch B1950.0
  FK4-based positions of stars to epoch J2000.0 positions in
  accordance with the new IAU resolutions".  Astron.Astrophys.
  128, 263-267.

- Seidelmann, P.K. (ed), 1992, "Explanatory Supplement to the
  Astronomical Almanac", ISBN 0-935702-68-7.

- Smith, C.A. et al., 1989, "The transformation of astrometric
  catalog systems to the equinox J2000.0".  Astron.J. 97, 265.

- Standish, E.M., 1982, "Conversion of positions and proper motions
  from B1950.0 to the IAU system at J2000.0".  Astron.Astrophys.,
  115, 1, 20-22.

- Yallop, B.D. et al., 1989, "Transformation of mean star places
  from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space".
  Astron.J. 97, 274.
"""
function fk425(r1950, d1950, dr1950, dd1950, p1950, v1950)
    r2000 = Ref{Cdouble}()
    d2000 = Ref{Cdouble}()
    dr2000 = Ref{Cdouble}()
    dd2000 = Ref{Cdouble}()
    p2000 = Ref{Cdouble}()
    v2000 = Ref{Cdouble}()
    ccall((:eraFk425, liberfa), Cvoid,
          (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble,
           Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
           Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},),
          r1950, d1950, dr1950, dd1950, p1950, v1950,
          r2000, d2000, dr2000, dd2000, p2000, v2000)
    return r2000[], d2000[], dr2000[], dd2000[], p2000[], v2000[]
end

"""
    fk45z(r1950, d1950, bepoch)

Convert a B1950.0 FK4 star position to J2000.0 FK5, assuming zero
proper motion in the FK5 system.

This function converts a star's catalog data from the old FK4
(Bessel-Newcomb) system to the later IAU 1976 FK5 (Fricke) system,
in such a way that the FK5 proper motion is zero.  Because such a
star has, in general, a non-zero proper motion in the FK4 system,
the function requires the epoch at which the position in the FK4
system was determined.

### Given ###

- `r1950`, `d1950`: B1950.0 FK4 RA,Dec at epoch (rad)
- `bepoch`: Besselian epoch (e.g. 1979.3)

### Returned ###

- `r2000`, `d2000`: J2000.0 FK5 RA,Dec (rad)

### Notes ###

1. The epoch bepoch is strictly speaking Beselian, but if a
   Julian epoch is supplied the result will be affected only to a
   negligible extent.

2. The method is from Appendix 2 of Aoki et al. (1983), but using
   the constants of Seidelmann (1992).  See the function eraFk425
   for a general introduction to the FK4 to FK5 conversion.

3. Conversion from equinox B1950.0 FK4 to equinox J2000.0 FK5 only
   is provided for.  Conversions for different starting and/or
   ending epochs would require additional treatment for precession,
   proper motion and E-terms.

4. In the FK4 catalog the proper motions of stars within 10 degrees
   of the poles do not embody differential E-terms effects and
   should, strictly speaking, be handled in a different manner from
   stars outside these regions.  However, given the general lack of
   homogeneity of the star data available for routine astrometry,
   the difficulties of handling positions that may have been
   determined from astrometric fields spanning the polar and non-
   polar regions, the likelihood that the differential E-terms
   effect was not taken into account when allowing for proper motion
   in past astrometry, and the undesirability of a discontinuity in
   the algorithm, the decision has been made in this ERFA algorithm
   to include the effects of differential E-terms on the proper
   motions for all stars, whether polar or not.  At epoch 2000.0,
   and measuring "on the sky" rather than in terms of RA change, the
   errors resulting from this simplification are less than
   1 milliarcsecond in position and 1 milliarcsecond per century in
   proper motion.

### References ###

- Aoki, S. et al., 1983, "Conversion matrix of epoch B1950.0
  FK4-based positions of stars to epoch J2000.0 positions in
  accordance with the new IAU resolutions".  Astron.Astrophys.
  128, 263-267.

- Seidelmann, P.K. (ed), 1992, "Explanatory Supplement to the
  Astronomical Almanac", ISBN 0-935702-68-7.

### Called ###

- [`anp`](@ref): normalize angle into range 0 to 2pi
- [`c2s`](@ref): p-vector to spherical
- [`epb2jd`](@ref): Besselian epoch to Julian date
- [`epj`](@ref): Julian date to Julian epoch
- [`pdp`](@ref): scalar product of two p-vectors
- [`pmp`](@ref): p-vector minus p-vector
- [`ppsp`](@ref): p-vector plus scaled p-vector
- [`pvu`](@ref): update a pv-vector
- [`s2c`](@ref): spherical to p-vectors
"""
function fk45z(r1950, d1950, bepoch)
    r2000 = Ref{Cdouble}()
    d2000 = Ref{Cdouble}()
    ccall((:eraFk45z, liberfa), Cvoid,
          (Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
          r1950, d1950, bepoch, r2000, d2000)
    return r2000[], d2000[]
end

"""
    function fk524(r2000, d2000, dr2000, dd2000, p2000, v2000)

Convert J2000.0 FK5 star catalog data to B1950.0 FK4.

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

- `r2000`, `d2000`: J2000.0 RA,Dec (rad)
- `dr2000`, `dd2000`: J2000.0 proper motions (rad/Jul.yr)
- `p2000`: parallax (arcsec)
- `v2000`: radial velocity (km/s, +ve = moving away)

### Returned (all B1950.0, FK4) ###

- `r1950`, `d1950`: B1950.0 RA,Dec (rad)
- `dr1950`, `dd1950`: B1950.0 proper motions (rad/trop.yr)
- `p1950`: parallax (arcsec)
- `v1950`: radial velocity (km/s, +ve = moving away)

### Notes ###

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

2. The conversion is somewhat complicated, for several reasons:

   - Change of standard epoch from J2000.0 to B1950.0.

   - An intermediate transition date of 1984 January 1.0 TT.

   - A change of precession model.

   - Change of time unit for proper motion (Julian to tropical).

   - FK4 positions include the E-terms of aberration, to simplify
     the hand computation of annual aberration.  FK5 positions
     assume a rigorous aberration computation based on the Earth's
     barycentric velocity.

   - The E-terms also affect proper motions, and in particular cause
     objects at large distances to exhibit fictitious proper
     motions.

   The algorithm is based on Smith et al. (1989) and Yallop et al.
   (1989), which presented a matrix method due to Standish (1982) as
   developed by Aoki et al. (1983), using Kinoshita's development of
   Andoyer's post-Newcomb precession.  The numerical constants from
   Seidelmann (1992) are used canonically.

4. In the FK4 catalog the proper motions of stars within 10 degrees
   of the poles do not embody differential E-terms effects and
   should, strictly speaking, be handled in a different manner from
   stars outside these regions.  However, given the general lack of
   homogeneity of the star data available for routine astrometry,
   the difficulties of handling positions that may have been
   determined from astrometric fields spanning the polar and non-
   polar regions, the likelihood that the differential E-terms
   effect was not taken into account when allowing for proper motion
   in past astrometry, and the undesirability of a discontinuity in
   the algorithm, the decision has been made in this ERFA algorithm
   to include the effects of differential E-terms on the proper
   motions for all stars, whether polar or not.  At epoch J2000.0,
   and measuring "on the sky" rather than in terms of RA change, the
   errors resulting from this simplification are less than
   1 milliarcsecond in position and 1 milliarcsecond per century in
   proper motion.

### Called ###

- [`anp`](@ref): normalize angle into range 0 to 2pi
- [`pdp`](@ref): scalar product of two p-vectors
- [`pm`](@ref): modulus of p-vector
- [`pmp`](@ref): p-vector minus p-vector
- [`ppp`](@ref): p-vector plus p-vector
- [`pv2s`](@ref): pv-vector to spherical coordinates
- [`s2pv`](@ref): spherical coordinates to pv-vector
- [`sxp`](@ref): multiply p-vector by scalar

### References ###

- Aoki, S. et al., 1983, "Conversion matrix of epoch B1950.0
  FK4-based positions of stars to epoch J2000.0 positions in
  accordance with the new IAU resolutions".  Astron.Astrophys.
  128, 263-267.

- Seidelmann, P.K. (ed), 1992, "Explanatory Supplement to the
  Astronomical Almanac", ISBN 0-935702-68-7.

- Smith, C.A. et al., 1989, "The transformation of astrometric
  catalog systems to the equinox J2000.0".  Astron.J. 97, 265.

- Standish, E.M., 1982, "Conversion of positions and proper motions
  from B1950.0 to the IAU system at J2000.0".  Astron.Astrophys.,
  115, 1, 20-22.

- Yallop, B.D. et al., 1989, "Transformation of mean star places
  from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space".
  Astron.J. 97, 274.
"""
function fk524(r2000, d2000, dr2000, dd2000, p2000, v2000)
    r1950 = Ref{Cdouble}()
    d1950 = Ref{Cdouble}()
    dr1950 = Ref{Cdouble}()
    dd1950 = Ref{Cdouble}()
    p1950 = Ref{Cdouble}()
    v1950 = Ref{Cdouble}()
    ccall((:eraFk524, liberfa), Cvoid,
          (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble,
           Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
           Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
          r2000, d2000, dr2000, dd2000, p2000, v2000,
          r1950, d1950, dr1950, dd1950, p1950, v1950)
    return r1950[], d1950[], dr1950[], dd1950[], p1950[], v1950[]
end

"""
    fk54z(r2000, d2000, bepoch)

Convert a J2000.0 FK5 star position to B1950.0 FK4, assuming zero
proper motion in FK5 and parallax.

### Given ###

- `r2000`, `d2000`: J2000.0 FK5 RA,Dec (rad)
- `bepoch`: Besselian epoch (e.g. 1950.0)

### Returned ###

- `r1950`, `d1950`: B1950.0 FK4 RA,Dec (rad) at epoch BEPOCH
- `dr1950`, `dd1950`: B1950.0 FK4 proper motions (rad/trop.yr)

### Notes ###

1. In contrast to the eraFk524 function, here the FK5 proper
   motions, the parallax and the radial velocity are presumed zero.

2. This function converts a star position from the IAU 1976 FK5
   (Fricke) system to the former FK4 (Bessel-Newcomb) system, for
   cases such as distant radio sources where it is presumed there is
   zero parallax and no proper motion.  Because of the E-terms of
   aberration, such objects have (in general) non-zero proper motion
   in FK4, and the present function returns those fictitious proper
   motions.

3. Conversion from B1950.0 FK4 to J2000.0 FK5 only is provided for.
   Conversions involving other equinoxes would require additional
   treatment for precession.

4. The position returned by this function is in the B1950.0 FK4
   reference system but at Besselian epoch BEPOCH.  For comparison
   with catalogs the BEPOCH argument will frequently be 1950.0. (In
   this context the distinction between Besselian and Julian epoch
   is insignificant.)

5. The RA component of the returned (fictitious) proper motion is
   dRA/dt rather than cos(Dec)*dRA/dt.

### Called ###

- [`anp`](@ref): normalize angle into range 0 to 2pi
- [`c2s`](@ref): p-vector to spherical
- [`fk524`](@ref): FK4 to FK5
- [`s2c`](@ref): spherical to p-vector
"""
function fk54z(r2000, d2000, bepoch)
    r1950 = Ref{Cdouble}()
    d1950 = Ref{Cdouble}()
    dr1950 = Ref{Cdouble}()
    dd1950 = Ref{Cdouble}()
    ccall((:eraFk54z, liberfa), Cvoid,
          (Cdouble, Cdouble, Cdouble,
           Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
          r2000, d2000, bepoch, r1950, d1950, dr1950, dd1950)
    return r1950[], d1950[], dr1950[], dd1950[]
end

"""
    fk5hip()

FK5 to Hipparcos rotation and spin.

### Returned ###

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

### Notes ###

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

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

         P_Hipparcos = r5h x P_FK5

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

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

### Called ###

- [`rv2m`](@ref): r-vector to r-matrix

### Reference ###

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

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

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

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

### Given ###

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

### Returned ###

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

### Notes ###

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

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

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

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

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

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

5. See also [`fk52h`](@ref), [`h2fk5`](@ref), [`hfk5z`](@ref).

### Called ###

- [`s2c`](@ref): spherical coordinates to unit vector
- [`fk5hip`](@ref): FK5 to Hipparcos rotation and spin
- [`sxp`](@ref): multiply p-vector by scalar
- [`rv2m`](@ref): r-vector to r-matrix
- [`trxp`](@ref): product of transpose of r-matrix and p-vector
- [`pxp`](@ref): vector product of two p-vectors
- [`c2s`](@ref): p-vector to spherical
- [`anp`](@ref): normalize angle into range 0 to 2pi

### Reference ###

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

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

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

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

### Given ###

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

### Returned ###

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

### Notes ###

1. Naming the following points:

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

   the four Fukushima-Williams angles are as follows:

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

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

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

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

### Called ###

- [`fw2m`](@ref): F-W angles to r-matrix
- [`bpn2xy`](@ref): extract CIP X,Y coordinates from NPB matrix

### Reference ###

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

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

"""
    fad03(t)

Fundamental argument, IERS Conventions (2003):

mean elongation of the Moon from the Sun.

### Given ###

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

### Returned ###

- `D`, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

"""
fad03

"""
    fae03(t)

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

### Given ###

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

### Returned ###

- Mean longitude of Earth, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

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

"""
fae03

"""
    faf03(t)

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

### Given ###

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

### Returned ###

- `F`, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

"""
faf03

"""
    faju03(t)

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

### Given ###

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

### Returned ###

- Mean longitude of Jupiter, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

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

"""
faju03

"""
    fal03(t)

Fundamental argument, IERS Conventions (2003):

mean anomaly of the Moon.

### Given ###

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

### Returned ###

- `l`, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

"""
fal03

"""
    falp03(t)

Fundamental argument, IERS Conventions (2003):

mean anomaly of the Sun.

### Given ###

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

### Returned ###

- `l'`, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

"""
falp03

"""
    fama03(t)

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

### Given ###

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

### Returned ###

- Mean longitude of Mars, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

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

"""
fama03

"""
    fame03(t)

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

### Given ###

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

### Returned ###

- Mean longitude of Mercury, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

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

"""
fame03

"""
    fane03(t)

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

### Given ###

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

### Returned ###

- Mean longitude of Neptune, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

"""
fane03

"""
    faom03(t)

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

### Given ###

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

### Returned ###

- `Omega`, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

"""
faom03

"""
    fapa03(t)

Fundamental argument, IERS Conventions (2003):

general accumulated precession in longitude.

### Given ###

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

### Returned ###

- General precession in longitude, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

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

"""
fapa03

"""
    fasa03(t)

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

### Given ###

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

### Returned ###

- Mean longitude of Saturn, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

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

"""
fasa03

"""
    faur03(t)

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

### Given ###

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

### Returned  ###

- Mean longitude of Uranus, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

"""
faur03

"""
    fave03(t)

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

### Given ###

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

### Returned ###

- Mean longitude of Venus, radians (Note 2)

### Notes ###

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

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

### References ###

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

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

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

"""
fave03

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

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

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

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

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

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

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

### Notes ###

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

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

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

4. See also [`h2fk5`](@ref), [`fk5hz`](@ref), [`hfk5z`](@ref).

### Called ###

- [`starpv`](@ref): star catalog data to space motion pv-vector
- [`fk5hip`](@ref): FK5 to Hipparcos rotation and spin
- [`rxp`](@ref): product of r-matrix and p-vector
- [`pxp`](@ref): vector product of two p-vectors
- [`ppp`](@ref): p-vector plus p-vector
- [`pvstar`](@ref): space motion pv-vector to star catalog data

### Reference ###

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

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

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

Form rotation matrix given the Fukushima-Williams angles.

### Given ###

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

### Returned ###

- `r`: Rotation matrix

### Notes ###

1. Naming the following points:

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

   the four Fukushima-Williams angles are as follows:

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

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

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

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

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

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

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

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

### Called ###

- [`ir`](@ref): initialize r-matrix to identity
- [`rz`](@ref): rotate around Z-axis
- [`rx`](@ref): rotate around X-axis

### Reference ###

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

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


text:
@time begin
    @time using GaussianBeamLens
    @time using GR

    @time X = GaussianBeamLens.GaussianBeamLensPropagate(0, 25, 1000, 1000);
    # See https://gr-framework.org/workstations.html
    #ENV["GKSwstype"] = 41 # Use quick native Win32 GUI
    @time GR.heatmap(X)
    println("Total time:") #  7.397603 seconds
end

max_stars_repo_name: mkitti/GaussianBeamLens

max_issues_repo_name: mkitti/GaussianBeamLens

max_forks_repo_name: mkitti/GaussianBeamLens


text:
module MutableStringTestUTF8
using MutableStrings

const nstr = 10^5
const lstr = 6

function makestr()
    iob = IOBuffer()
    for idx in 1:nstr
        write(iob, randstring(lstr))
        write(iob, ' ')
    end
    takebuf_string(iob)
end

#mk_blank(x::MutableASCIIString, y::Range1{Int64}) = (x.data[y] = uint8(' '))
mk_blank(x) = ""
mk_space(x) = " "

function warmup()
    _mstr = MutableUTF8String(copy("abcdef ΑΒΓΔΕΖΗΘ αβγδεζηθ".data))
    l = length(_mstr)
    lowercase!(_mstr);
    @assert _mstr == "abcdef αβγδεζηθ αβγδεζηθ"
    #println(_mstr)
    uppercase!(_mstr);
    @assert _mstr == "ABCDEF ΑΒΓΔΕΖΗΘ ΑΒΓΔΕΖΗΘ"
    #println(_mstr)

    reverse!(_mstr);
    lcfirst!(_mstr);
    @assert _mstr[1] == 'θ'
    ucfirst!(_mstr);
    @assert _mstr[1] == 'Θ'

    @assert _mstr[end] == 'A'
    @assert length(_mstr) == l
    bl = length(_mstr.data)
    _mstr[end] = 'Δ'
    @assert _mstr[end] == 'Δ'
    @assert length(_mstr) == l
    @assert length(_mstr.data) > bl
    bl = length(_mstr.data)
    @assert _mstr[end-1] == 'B'
    _mstr[end-1] = 'β'
    @assert _mstr[end-2] == 'β'
    @assert _mstr == "ΘΗΖΕΔΓΒΑ ΘΗΖΕΔΓΒΑ FEDCβΔ"
end

function test_set_get()
    ms = MutableUTF8String("Lorem ipsum dolor sit amet, consectetur adipisicing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum.")
    @assert ms[1] == 'L'
    @assert ms[1:10] == "Lorem ipsu"
    @assert ms[2:10] == "orem ipsu"
    @assert typeof(ms[1:10]) == UTF8String
    ms[1] = 'X'
    @assert ms[1:10] == "Xorem ipsu"
    ms[1:10] = 'X'
    @assert ms[1:10] == "X"^10
    ms[1:10] = "Lorem ipsu"
    @assert ms[1:10] == "Lorem ipsu"

    replace!(ms, "Lorem", "O RUM")
    @assert ms[1:5] == "O RUM"
    replace!(ms, r"O RUM", "Lorem")
    @assert ms[1:10] == "Lorem ipsu"

    s = replace(ms, "Lorem", "O RUM")
    @assert ms[1:5] != "O RUM"
    @assert beginswith(s, "O RUM")

    replace!(ms, "adipisicing", "FEDCβΔ")
    @assert start(search(ms, "adipisicing")) == 0
    @assert start(search(ms, "FEDCβΔ")) > 0

    replace!(ms, r"FEDCβΔ", "adipisicing")
    @assert start(search(ms, "adipisicing")) > 0
    @assert start(search(ms, "FEDCβΔ")) == 0
end

function test_time()
    master_str = makestr();

    str = utf8(master_str);
    mstr = MutableUTF8String(str);
    @printf("%25s%30s%40s\n", "", "string", "mutable string")
    @printf("%25s%30s%40s\n", "length", "$(length(str))", "$(length(mstr))")
    @printf("%25s%20s%20s%20s%20s\n", "", "time", "bytes", "time", "bytes")

    gc(); _ret, t, b = @timed str = lowercase(str);
    gc(); _mret, mt, mb = @timed lowercase!(mstr);
    @printf("%25s%20s%20s%20s%20s\n", "lowercase", "$t", "$b", "$mt", "$mb")

    gc(); _ret, t, b = @timed str = reverse(str);
    gc(); _mret, mt, mb = @timed reverse!(mstr);
    @printf("%25s%20s%20s%20s%20s\n", "reverse", "$t", "$b", "$mt", "$mb")

end
end # module

MutableStringTestUTF8.warmup()
MutableStringTestUTF8.test_set_get()
MutableStringTestUTF8.test_time()


text:
# sec2delete = "# ╟─2b4a8e35-7cef-4eed-bad2-29b9290756b6
# # ╟─e488cf57-0948-48a8-9caf-c18405ecfdae
# # ╟─19e15dda-281e-4e3d-82e7-9b1ae3ae2fa1
# # ╠═665e56b1-d063-4d58-9ff4-61b9ef914641
# # ╠═3bcd446d-652c-4b4c-ac37-7ad70009a831
# # ╟─0febdbc8-3829-44e2-adfc-7b4656549a5b
# # ╠═623c3306-6fbb-4500-8675-a6c66dfac2aa
# # ╠═9fb0ddcd-ee66-42fa-ac14-5ea6db4c5b14
# # ╠═dad2eee8-166a-46f1-8f0b-1abf0a490eda
# # ╟─fdb9b72c-b107-42e1-844e-5ee33f626a13
# # ╠═716b3c84-d02b-4ea0-81d9-d679ceb53503
# # ╟─43e800b7-ef12-44e5-923c-e42323a92872
# # ╠═8d6d11f4-4329-4a7b-ae0e-5de7db67f88c
# # ╟─aee72572-eddb-4a74-810b-5cbdac7325dc
# # ╠═6cb8decc-d3ad-408f-96fa-2c667cf73c38
# # ╟─4e1538f1-a134-4cf0-913c-7f5e9f2864ab
# # ╠═71a4f0ca-f6db-486f-af45-7914ef71b90d
# # ╠═48c140c1-0e29-4dae-bcaf-d816bebcbdb7
# # ╟─e002ef9d-c03b-498c-b830-796f6bf4a5e3
# # ╟─dfd9a282-1b1d-428c-b306-8a44edbf0c3e
# # ╠═037f673c-5a05-4e90-84ee-f10d4769e5ca
# # ╠═378c2188-7cb3-44cd-87d6-7d35a2eb594c"
function sectiondeletePluto(fpath::String,newpath::String, sec2delete)
lines = open(fpath) do file
	readlines(file);
end

id2 = lines .== "# ╔═╡ Cell order:";
lenid2 = findall(id2);
lnco = length(lenid2)==1 ? lenid2[1] : error("Something wrong");

part1 = lines[1:lnco-1];
part2 = lines[lnco:end];
key1s = get_keys(part1, "# ╔═╡ ");
key2s = get_keys(part2, "(# ╠═|# ╟─)");

sort(key1s) == sort(key2s) ? nothing : error("Keys not matched.")

id0s = findall(occursin.(r"# ╔═╡ ",part1));
id1s = findall(occursin.(r"# ╔═╡ ",part1)) .- 1;
popfirst!(id1s);
push!(id1s,lnco-1);
ids1 = hcat(id0s, id1s);
ids2 = findall(occursin.(r"(# ╠═|# ╟─)",part2));

keys2delete = collect(eachmatch(Regex("\\w{8}-\\w{4}-\\w{4}-\\w{4}-\\w{6}"),sec2delete))
keys2delete = [mt.match for mt in keys2delete];

length(key1s) == length(id0s) ? nothing : error("Section keys and section line numbers mismatch.");

delete_lines!(key1s, part1, keys2delete,ids1);
delete_lines!(key2s, part2, keys2delete, ids2);

open(newpath,"w") do io
	for line in vcat(part1, part2)
		println(io, line);
	end
end

end

function delete_lines!(key1s, part1, keys2delete, ids)
	srng2delete_pt1 = [];
	for keyd in keys2delete
		push!(srng2delete_pt1, ids[occursin.(Regex(keyd), key1s), :]);
	end

	deleteat!(part1, vcat(collect.(map(x -> x[1]:x[end], srng2delete_pt1))...));
end

function get_keys(partn,prefix)
	match1 = match.(Regex("(?<=$prefix)\\w{8}-\\w{4}-\\w{4}-\\w{4}-\\w{6}"),partn);
	keys = [];
	for mt in match1
		!isnothing(mt) ? push!(keys,mt.match) : continue
	end
	return keys
end


text:
mutable struct Problem
  F::Firm

  N::Int
  S::Int
  p::Matrix{Float64}

  dtol::Float64
  maxiter::Int
end

function updatefirmτ!(Φ::Problem, τ::Vector{Float64})::Problem
  Φ.F.τ .= τ
  refreshR!(Φ.F, Φ.p)
  findpolicy!(Φ::Problem)
  return Φ
end

function findpolicy!(Φ::Problem)::Problem
  local rslices::Matrix{SubArray} #view of r in each state
  local vslices::Matrix{SubArray} #view of v in each state
  local πslices::Matrix{SubArray} #view of the policy in each state
  local mslices::Array{SubArray,3} #view of the pricing kernel in each state and each future state
  local iter::Int = 0
  local derr::Float64 = 1.0
  local vrow::Matrix{Float64} = Matrix{Float64}(undef,1,Φ.N) #captures updated values
  local πrow::Matrix{CartesianIndex} = Matrix{CartesianIndex}(undef,1,Φ.N) #captures updated policy
  local vold::Array{Float64,3} = deepcopy(Φ.F.V) #holds the previous value for comparison purposes

  #each view here maps all current states to a single future state
  rslices = Matrix{SubArray}(undef, Φ.S, Φ.S)
  mslices = Array{SubArray,3}(undef, Φ.S, Φ.S, Φ.S)
  vslices = Matrix{SubArray}(undef, Φ.S, Φ.S) #will a view for each lagged and current state
  πslices = Matrix{SubArray}(undef, Φ.S, Φ.S)

  #make our convenience views
  for sₗ ∈ 1:Φ.S, s ∈ 1:Φ.S
    rslices[sₗ, s] = view(Φ.F.R, :, :, sₗ, s)
    vslices[sₗ, s] = view(Φ.F.V, :, sₗ, s)
    πslices[sₗ, s] = view(Φ.F.Π, :, sₗ, s)
    mslices[sₗ, s,:] .= ((sₜ₁::Int)->view(Φ.F.M, :, :, sₗ, s, sₜ₁)).(1:Φ.S)
  end

  #these are useful for the BLAS libraries
  C::Matrix{Float64} = zeros(Φ.N, Φ.N)
  β::Float64 = 0.0
  N1::Vector{Float64} = ones(Φ.N)
  N²::Int = Φ.N*Φ.N

  while (derr > Φ.dtol) && (iter < Φ.maxiter)
    for sₗ ∈ 1:Φ.S, s ∈ 1:Φ.S #need to use vec since findmax returns column vectors

      vrow, πrow = findmax(rslices[sₗ, s] .+
        (Φ.p[1, s] .* mslices[sₗ, s,1] .* vslices[s, 1]) .+
        (Φ.p[2, s] .* mslices[sₗ, s,2] .* vslices[s, 2]), dims=1)

      #need to do this to fix the dimensionality of the returned values from findmax
      vslices[sₗ, s] .= vec(vrow)
      πslices[sₗ, s] .= vec(πrow)

    end

    derr = maximum((abs).(vold .- Φ.F.V))
    vold .= Φ.F.V

    iter += 1

    (iter % 500 == 0) && println("iter: $iter, derr: $derr")

  end

  (iter == Φ.maxiter) && @warn("Convergence failure iter: $iter, derr: $derr")

  return Φ
end

function Problem(F::Firm{NTuple{2,Float64}}, p::Matrix{Float64};
    dtol::Float64=DEF_D_TOL, maxiter::Int = DEF_MAX_ITER,
    verbose::Bool = length(F.ψ[1]) < 1000)::Problem

  N::Int = length(F.ψ[1])
  S::Int = size(p,1)

  Φ::Problem = Problem(F, N, S, p, dtol, maxiter)
  findpolicy!(Φ)

  #display(Φ.F.V)
  if verbose
    highpolicy::Matrix{Float64} = Matrix{Int}(undef, Φ.N, Φ.S)
    lowpolicy::Matrix{Float64} = Matrix{Int}(undef, Φ.N, Φ.S)

    highvalue::Matrix{Float64} = Matrix{Float64}(undef, Φ.N, Φ.S)
    lowvalue::Matrix{Float64} = Matrix{Float64}(undef, Φ.N, Φ.S)
    for s ∈ 1:Φ.S
      highpolicy[:,s] .= ((n)->(maximum(Φ.F.R[:,n,1,s]) > 0.0) ? (Φ.F.Π[n,1,s])[1] : -1.0).(1:Φ.N)
      lowpolicy[:,s] .= ((n)->(maximum(Φ.F.R[:,n,2,s]) > 0.0) ? (Φ.F.Π[n,2,s])[1] : -1.0).(1:Φ.N)

      highvalue[:,s] = ((n)->(maximum(Φ.F.R[:,n,1,s]) > 0.0) ? (Φ.F.V[n,1,s]) : -1.0).(1:Φ.N)
      lowvalue[:,s]  = ((n)->(maximum(Φ.F.R[:,n,2,s]) > 0.0) ? (Φ.F.V[n,2,s]) : -1.0).(1:Φ.N)

      highvalue[:,s] .= (v::Float64->round(v,digits=1)).(highvalue[:,s])
      lowvalue[:,s] .= (v::Float64->round(v,digits=1)).(lowvalue[:,s])
    end

    for n ∈ 1:Φ.N
      k::Float64 = round(Φ.F.ψ[1][n], digits=2)
      yh::Float64 = round(Φ.F.ψ[2][n], digits=2)

      ylgh::Float64 = getyl(Φ.F.ψ[2][n], Φ.F.ψ[1][n], Φ.p[1,1], Φ.F)
      ylgl::Float64 = getyl(Φ.F.ψ[2][n], Φ.F.ψ[1][n], Φ.p[1,2], Φ.F)

      ylgh = round(ylgh, digits=2)
      ylgl = round(ylgl, digits=2)

      #println("n: $n; k,yh: ($k, $yh); yl, v,p|sₗ=h: $ylgh, $(highvalue[n,:]), $(highpolicy[n,:])")
      println("n: $n; k,yh: ($k, $yh); yl, p|sₗ=h: $ylgh, $(highpolicy[n,:]); p|sₗ=l: $ylgl, $(lowpolicy[n,:] )")
    end
  end

  return Φ
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using TrajectoryOptimization
using BenchmarkTools
using Dates
using Logging
using LinearAlgebra
import TrajectoryOptimization: trim_entry
import BenchmarkTools: prettytime, prettymemory, prettydiff, prettypercent
using Formatting

const paramsfile = "benchmark/benchmark_params.json"
const histfile = "benchmark/benchmark_history.json"

include("benchmark_methods.jl")
# include("dubinscar_benchmarks.jl")
# include("pendulum_benchmarks.jl")
include("ilqr_benchmarks.jl")

loadparams!(SUITE, BenchmarkTools.load(paramsfile)[1])
# results = run(suite, verbose = true)


text:
export StrainReconParam, getB, getStrainReconParam

"""
type StrainReconParam

Data and temporary variables for strain reconstruction problem

Fields:

	1. d::Vector           - measurement vector (SNP data)
	2. m::Int              - number of SNP sites
	3. n::Int              - number of strains to be estimated
	4. p::Int              - number of categories (i.e., possible values for one SNP site)
	5. q::Int              - length of the measurement vector, q = m*(p-1)
	6. pwr                 - exponent in the likelihood function
	7. R                   - Weight matrix for the likelihood as the inverse square root of Gamma
	8. label               - label
	9. lambda              - regularization parameter (usually zero)
	10. B::Array{Bool}     - (p-1)*p^n x n binary matrix of all possible blocks of size (p-1) x n
	11.-14. Ci, ci, Ce, ce - (in)equalities to be used in optimization algorithms
"""
type StrainReconParam
	d::Vector{Float64}
	m::Int
	n::Int
	p::Int
	q::Int
	pwr::Float64
	R::SparseMatrixCSC{Float64}
	label::String
	lambda::Float64
	B::BitArray{2}
	Ci::AbstractArray{Float64}
	ci::Vector{Float64}
	Ce::AbstractArray{Float64}
	ce::Vector{Float64}
end

function Base.display(param::StrainReconParam)
	print("StrainReconParam(")
	(!isempty(param.label)) && print("label=$(param.label), ")
	print("length(d)=$(length(param.d)), ")
	print("nStrains=$(param.n))\n")
	print("nCategories=$(param.p))\n")
	print("nSites=$(param.m))\n")
end

"""
function getStrainReconParam(d,n;p,pwr,R,label,lambda,B,Ci,ci,Ce,ce)

Construct and return a StrainReconParam object.

Input: see StrainReconParam. The weight matrix R can also be given as a vector
(of length q) which is then converted to a diagonal matrix, or as a scalar if
the diagonal is constant.
"""
function getStrainReconParam(d::Vector{Float64},n::Int; p::Int=2, pwr::Float64=2.0,
	                     R=sdiag(ones(size(d,1))),label::String="",lambda::Float64=0.0,
	                     B::BitArray{2}=falses(0,0),
	                     Ci::AbstractArray{Float64}=zeros(0,n),ci::Vector{Float64}=zeros(0),
						 Ce::AbstractArray{Float64}=zeros(0,n),ce::Vector{Float64}=zeros(0))

	if length(d) % (p-1) != 0
		error("the length of data must be divisible by nCategories-1")
	else
		m = Int(length(d)/(p-1))
		q = m*(p-1)
	end

	if length(R) == q
		R = sdiag(R)
	elseif length(R) == 1
		R = R*sdiag(ones(q))
	end

	return StrainReconParam(d,m,n,p,q,pwr,R,label,lambda,B,Ci,ci,Ce,ce)
end

"""
function getConstraints(param::StrainReconParam)

constructs contraint matrices and stores them inside param. Recomputes them only if needed.
"""
function getConstraints(param::StrainReconParam)
	if isempty(param.Ci) || isempty(param.Ce) || isempty(param.ci) || isempty(param.ce) ||
		size(param.Ci,2)!=param.n || size(param.Ce,2)!=param.n

		# equality constraint
		param.Ce = ones(1,param.n)
		param.ce = ones(1)

		# lower bound
        Ci1 = eye(param.n)[end,:]
        ci1 = zeros(param.n)[end]

		# monotonicity
		Ci2 = - full(ddx(param.n-1))
		ci2 = zeros(param.n-1)
		param.Ci = [Ci1';Ci2]
		param.ci = [ci1;ci2]
	end
	return param.Ci, param.ci, param.Ce, param.ce
end

"""
function getB(param::StrainReconParam)

returns (p-1)*p^n x n block binary matrix (block size
is (p-1) x n), corresponding to all possible p-ary
barcodes of length n (all possible blocks of matrix M)

Matrix is only constructed when needed and then stored
in param.B

Input:

	param::StrainReconParam

Output:
	B::Array{Bool} - size(B)=((p-1)*p^param.n,param.n)
"""
function getB(param::StrainReconParam)
	if isempty(param.B) || size(param.B,1)!=(param.p-1)*param.p^param.n || size(param.B,2)!=param.n
		param.B = getB(param.n,param.p)
	end
	return param.B
end

function getB(n::Int,p::Int=2)
	# create all possible rows of n elements with integers 1...p
	intmat = zeros(Int,p^n,n)
	idx = 1
	for i = 0:n-1
		for j = 0:p^n-1
			intmat[idx] = mod(floor(Int,j/p^i),p)+1
			idx += 1
		end
	end

	# convert each row to a binary block with p-1 rows (the last row would be redundant)
	B = falses((p-1)*p^n,n)
	for k=1:size(intmat,1)
		fullblock = falses(p,n)
		fullblock[sub2ind(size(fullblock),intmat[k,:],1:n)] = true
		B[(k-1)*(p-1)+1 : k*(p-1),:] = fullblock[1:p-1,:]
	end

	# arrange so that falses go to the upper left corner
	B = flipdim(flipdim(B,1),2)

	return B
end


text:
module Diffusion

using QuantExQASM.NCU
using QuantExQASM.Oracle
using QuantExQASM.Circuit
using QuantExQASM.GateOps

function apply_h!(cct::Circuit.Circ, tgt, reg::Union{String, Nothing}=nothing)
    Circuit.add_gatecall!(cct, GateOps.hadamard(tgt, reg))
end

function apply_x!(c::Circuit.Circ, tgt, reg::Union{String, Nothing}=nothing)
    Circuit.add_gatecall!(c, GateOps.pauli_x(tgt, reg))
end

"""
    apply_diffusion(cct::Circuit.Circ, ctrl_indices::Vector, tgt_index, reg::Union{String, Nothing}=nothing)

Application of the Grover diffusion operator to marked register.
"""
function apply_diffusion(cct::Circuit.Circ, ctrl_indices::Vector, tgt_index, reg::Union{String, Nothing}=nothing)
    aux_idx = typeof(ctrl_indices)()

    for ctrl in vcat(ctrl_indices, tgt_index)
        apply_h!(cct, ctrl, reg)
        apply_x!(cct, ctrl, reg)
    end

    NCU.apply_ncu!(cct, ctrl_indices, aux_idx, tgt_index, GateOps.GateLabel(:z))

    for ctrl in vcat(ctrl_indices, tgt_index)
        apply_x!(cct, ctrl, reg)
        apply_h!(cct, ctrl, reg)
    end
end

"""
    apply_diffusion(cct::Circuit.Circ, ctrl_indices::Vector, aux_indices::Vector, tgt_index, reg::Union{String, Nothing}=nothing)

Application of the Grover diffusion operator to marked register. Uses additionally provided auxiliary qubits to reduce depth.
"""
function apply_diffusion(cct::Circuit.Circ, ctrl_indices::Vector, aux_indices::Vector, tgt_index, reg::Union{String, Nothing}=nothing)

    for ctrl in vcat(ctrl_indices, tgt_index)
        apply_h!(cct, ctrl, reg)
        apply_x!(cct, ctrl, reg)
    end

    NCU.apply_ncu!(cct, ctrl_indices, aux_indices, tgt_index, GateOps.GateLabel(:z))

    for ctrl in vcat(ctrl_indices, tgt_index)
        apply_x!(cct, ctrl, reg)
        apply_h!(cct, ctrl, reg)
    end
end

end