text:
using LinearAlgebra
using Printf

abstract type TimePropagationMethod end

COUNT_B = 0

mutable struct EquidistantTimeStepper
    method::TimePropagationMethod
    psi::WaveFunction
    t0::Real
    dt::Real
    steps::Int
    function EquidistantTimeStepper(method::TimePropagationMethod, psi::WaveFunction,
        t0::Real, dt::Real, steps::Int)
        global COUNT_B = 0
        new(method, psi, t0, dt, steps)
    end
end

#default initializer
function initialize!(method::TimePropagationMethod, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int)
    0 # iteration  starts with step 0
end

#default finalizer
function finalize!(method::TimePropagationMethod, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int, step::Int)
end

function Base.iterate(tsi::EquidistantTimeStepper, step=0)
    if step == 0
        set_propagate_time_together_with_A!(tsi.psi.m, true)
        set_time!(tsi.psi, tsi.t0)
        step = initialize!(tsi.method, tsi.psi, tsi.t0, tsi.dt, tsi.steps)
    end

    if step >= tsi.steps
        finalize!(tsi.method, tsi.psi, tsi.t0, tsi.dt, tsi.steps, step)
        return nothing
    end

    step!(tsi.method, tsi.psi, tsi.t0, tsi.dt, tsi.steps, step)
    return (step+1, tsi), step+1
end

function global_orders(method::TimePropagationMethod,
                       psi::WaveFunction, reference_solution::WaveFunction,
                       t0::Real, tend::Real, dt::Real; rows=8)
    @assert psi.m==reference_solution.m
    tab = zeros(Float64, rows, 3)

    wf_save_initial_value = clone(psi)
    copy!(wf_save_initial_value, psi)

    steps = Int(floor((tend-t0)/dt))
    dt1 = dt
    err_old = 0.0
    println("             dt         err           C      p      B calls")
    println("-----------------------------------------------------------")
    for row=1:rows
        for t in EquidistantTimeStepper(method, psi, t0, dt1, steps) end
        err = distance(psi, reference_solution)
        if (row==1)
            @Printf.printf("%3i%12.3e%12.3e                   %13i\n", row, Float64(dt1), Float64(err), Int64(COUNT_B))
            tab[row,1] = dt1
            tab[row,2] = err
            tab[row,3] = 0
        else
            p = log(err_old/err)/log(2.0)
            C = err/dt1^p
            @Printf.printf("%3i%12.3e%12.3e%12.3e%7.2f%13i\n", row, Float64(dt1), Float64(err),
                                                       Float64(C), Float64(p), Int64(COUNT_B))
            tab[row,1] = dt1
            tab[row,2] = err
            tab[row,3] = p
        end
        err_old = err
        dt1 = 0.5*dt1
        steps = 2*steps
        copy!(psi,wf_save_initial_value)
    end
end

#######################################
#Splitting methods
#######################################

mutable struct SplittingMethod <: TimePropagationMethod
    s::Int
    a::Vector{Float64}
    b::Vector{Float64}
    function SplittingMethod(a::Vector{Float64}, b::Vector{Float64})
         s = length(a)
         @assert s==length(b)
         new(s, a, b)
    end
end

function step!(m::SplittingMethod, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int, step::Int)
    for j = 1:m.s
        if m.a[j]!=0.0
            propagate_A!(psi, m.a[j]*dt)
        end
        if m.b[j]!=0.0
            propagate_B!(psi, m.b[j]*dt)
            global COUNT_B += 1
        end
    end
end

# Splitting with RK4 for B step (alternatively lower order RK2 for cheaper lower order methods (strang))
mutable struct SplittingRK4BMethod <: TimePropagationMethod
    s::Int
    a::Vector{Float64}
    b::Vector{Float64}
    secondorder::Bool
    G #::WaveFunction
    U #::Vector{WaveFunction}
    function SplittingRK4BMethod(a::Vector{Float64}, b::Vector{Float64}; secondorder::Bool=false)
         s = length(a)
         @assert s==length(b)
         new(s, a, b, secondorder)
    end
end

function initialize!(method::SplittingRK4BMethod, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int)
    method.G = wave_function(psi.m)
    method.U = WaveFunction[wave_function(psi.m) for j=1:4]
    0 # iteration  starts with step 0
end

function finalize!(method::SplittingRK4BMethod, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int, step::Int)
    method.G = nothing
    method.U = nothing
end

function step!(m::SplittingRK4BMethod, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int, step::Int)
    for j = 1:m.s
        if m.a[j]!=0.0
            propagate_A!(psi, m.a[j]*dt)
        end
        if m.b[j]!=0.0
            if !m.secondorder
                RK4_stepB!(psi, m.b[j]*dt, G=m.G, U=m.U)
            else
                RK2_stepB!(psi, m.b[j]*dt, G=m.G, U=m.U)
            end
        end
    end
end

function RK4_stepB!(psi::WaveFunction, dt::Number;
    # classical explicit Runge-Kutta method for propagating B by dt, using time at t+dt as implied
    # by autonomization for splitting methods.
    G::WaveFunction=wave_function(psi.m),
    U::Vector{WaveFunction}=[wave_function(psi.m) for j=1:4])

    t = get_time(psi)
    set_time!(psi, t+dt)
    for j=1:4
        copy!(U[j], psi)
        set_time!(U[j], t)
    end
    set_time!(psi, t)

    gen_rhs!(G, U[1])
    axpy!(U[2], G, 0.5*dt)
    axpy!(psi, G, 1/6*dt)

    gen_rhs!(G, U[2])
    axpy!(U[3], G, 0.5*dt)
    axpy!(psi, G, 1/3*dt)

    gen_rhs!(G, U[3])
    axpy!(U[4], G, dt)
    axpy!(psi, G, 1/3*dt)

    gen_rhs!(G, U[4])
    axpy!(psi, G, 1/6*dt)
    psi
end

function RK2_stepB!(psi::WaveFunction, dt::Number;
    # Second order Runge-Kutta method for propagating B by dt, using time at t+dt as implied
    # by autonomization for splitting methods.
    G::WaveFunction=wave_function(psi.m),
    U::Vector{WaveFunction}=[wave_function(psi.m) for j=1:2])

    t = get_time(psi)
    set_time!(psi, t+dt)
    for j=1:2
        copy!(U[j], psi)
        set_time!(U[j], t)
    end
    set_time!(psi, t)

    gen_rhs!(G, U[1])
    axpy!(U[2], G, dt)
    axpy!(psi, G, 1/2*dt)

    gen_rhs!(G, U[2])
    axpy!(psi, G, 1/2*dt)

    psi
end

#######################################
#Composition methods
#######################################
function gen_rhs!(rhs::WaveFunction, psi::WaveFunction)
    set!(rhs, 0)
    add_apply_B!(psi,rhs)
    global COUNT_B += 1
end

function get_coeffs_composition(g::Vector{Float64})
    n = length(g)
    a = zeros(n+1)
    b = vcat(g, 0)
    a[1] = g[1]/2
    a[n+1] = a[1]
    for j=2:div(n,2)+1
        a[j] = (g[j-1]+g[j])/2
        a[n-j+2] = a[j]
    end
    a,b
end

mutable struct CompositionMethod <: TimePropagationMethod
    s::Int
    c::Vector{Float64}
    a::Vector{Float64}
    b::Vector{Float64}
    iters::Int
    psi0 #::WaveFunction
    k1 #::WaveFunction
    k2 #::WaveFunction
    function CompositionMethod(c::Vector{Float64}, iters::Int)
         s = length(c)+1
         a,b = get_coeffs_composition(c)
         new(s, c, a, b, iters, nothing, nothing, nothing)
    end
end

function initialize!(method::CompositionMethod, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int)
    method.k1 = wave_function(psi.m)
    method.k2 = wave_function(psi.m)
    method.psi0 = wave_function(psi.m)
    0 # iteration  starts with step 0
end

function finalize!(method::CompositionMethod, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int, step::Int)
    method.k1 = nothing
    method.k2 = nothing
    method.psi0 = nothing
end

function step!(m::CompositionMethod, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int, step::Int)
    for j = 1:m.s
        if m.a[j]!=0.0
            propagate_A!(psi, m.a[j]*dt)
        end
        if m.b[j]!=0.0
            if m.iters>0
                copy!(m.psi0, psi)
            end

            #Euler step:
            gen_rhs!(m.k1, psi)
            axpy!(psi, m.k1, m.b[j]*dt)

            #fixed point iteration for implicit midpoint rule
            for it=1:m.iters
                copy!(m.k1, m.psi0)
                scale!(m.k1, 0.5)
                axpy!(m.k1, psi, 0.5)
                copy!(psi, m.psi0)
                set_time!(m.k1, get_time(psi))
                gen_rhs!(m.k2, m.k1)
                axpy!(psi, m.k2, m.b[j]*dt)
            end
        end
    end
end

#######################################
#(Exponential) Runge-Kutta
#######################################

mutable struct ExponentialRungeKutta <: TimePropagationMethod
    scheme::Int
    s::Int
    iters::Int
    G #::WaveFunction
    U #::Vector{WaveFunction}
    function ExponentialRungeKutta(scheme::Symbol=:krogstad; iters::Int=(scheme==:gauss_lawson ? 4 : 0))
        scheme1 = 42
        s = 4
        if scheme==:rk4
            scheme1=4
            s = 4
        elseif scheme==:erk2
            scheme1=39
            s = 2
        elseif scheme==:erk2a
            scheme1=40
            s = 2
        elseif scheme==:lawson
            scheme1=14
            s = 4
        elseif scheme==:krogstad
            scheme1=42
            s = 4
    #    elseif scheme==:cox_mathews
    #        scheme1=41
    #        s = 4
        elseif scheme==:strehmel_weiner
            scheme1=43
            s = 4
        elseif scheme==:gauss_lawson
            scheme1=24
            s = 3
        else
            warn("scheme not known, I use :krogstad")
        end
        new(scheme1, s, iters, nothing, nothing)
    end
end

function initialize!(method::ExponentialRungeKutta, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int)
    method.G = wave_function(psi.m)
    method.U = WaveFunction[wave_function(psi.m) for j=1:method.s]
    0 # iteration  starts with step 0
end

function finalize!(method::ExponentialRungeKutta, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int, step::Int)
    method.G = nothing
    method.U = nothing
end

function ERK2_step!(psi::WaveFunction, dt::Number;
    #eq. (2.39) in Hochbruck/Ostermann, c2=1
    G::WaveFunction=wave_function(psi.m),
    U::Vector{WaveFunction}=[wave_function(psi.m) for j=1:2])

    s = 2
    c = [0.0, 1.0]*dt
    t = get_time(psi)
    for j=1:s
        copy!(U[j], psi)
        set_time!(U[j], t+c[j])
    end
    set_time!(psi, t+dt)

    gen_rhs!(G, U[1])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[2], c[2], 1, dt)
    add_phi_A!(G, psi, dt, 1, dt)
    add_phi_A!(G, psi, dt, 2, -dt)

    gen_rhs!(G, U[2])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, psi, dt, 2, dt)
end

function ERK2A_step!(psi::WaveFunction, dt::Number;
    #eq. (2.40) in Hochbruck/Ostermann, c2=0.5
    G::WaveFunction=wave_function(psi.m),
    U::Vector{WaveFunction}=[wave_function(psi.m) for j=1:2])

    s = 2
    c = [0.0, 0.5]*dt
    t = get_time(psi)
    for j=1:s
        copy!(U[j], psi)
        set_time!(U[j], t+c[j])
    end
    set_time!(psi, t+dt)

    gen_rhs!(G, U[1])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[2], c[2], 1, 0.5*dt)

    gen_rhs!(G, U[2])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, psi, dt, 1, dt)
end

function ERK4_step_krogstad!(psi::WaveFunction, dt::Number;
    #Krogstad, eq. (2.42) in Hochbruck/Ostermann
    G::WaveFunction=wave_function(psi.m),
    U::Vector{WaveFunction}=[wave_function(psi.m) for j=1:4])
    s = 4
    c = [0.0, 0.5, 0.5, 1.0]*dt
    t = get_time(psi)
    for j=1:s
        copy!(U[j], psi)
        set_time!(U[j], t+c[j])
    end
    set_time!(psi, t+dt)

    gen_rhs!(G, U[1])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[2], c[2], 1, 0.5*dt)
    add_phi_A!(G, U[3], c[3], 1, 0.5*dt)
    add_phi_A!(G, U[3], c[3], 2, -1.0*dt)
    add_phi_A!(G, U[4], c[4], 1, 1.0*dt)
    add_phi_A!(G, U[4], c[4], 2, -2.0*dt)
    add_phi_A!(G, psi, dt, 1, 1.0*dt)
    add_phi_A!(G, psi, dt, 2, -3.0*dt)
    add_phi_A!(G, psi, dt, 3, 4.0*dt)

    gen_rhs!(G, U[2])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[3], c[3], 2, 1.0*dt)
    add_phi_A!(G, psi, dt, 2, 2.0*dt)
    add_phi_A!(G, psi, dt, 3, -4.0*dt)

    gen_rhs!(G, U[3])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[4], c[4], 2, 2.0*dt)
    add_phi_A!(G, psi, dt, 2, 2.0*dt)
    add_phi_A!(G, psi, dt, 3, -4.0*dt)

    gen_rhs!(G, U[4])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, psi, dt, 2, -1.0*dt)
    add_phi_A!(G, psi, dt, 3, 4.0*dt)
    psi
end

function RK4_step!(psi::WaveFunction, dt::Number;
    #classical explicit Runge-Kutta method (NOT exponential!)
    G::WaveFunction=wave_function(psi.m),
    U::Vector{WaveFunction}=[wave_function(psi.m) for j=1:4])
    s = 4
    c = [0.0, 0.5, 0.5, 1.0]*dt
    t = get_time(psi)
    for j=1:s
        copy!(U[j], psi)
        set_time!(U[j], t+c[j])
    end
    set_time!(psi, t+dt)

    gen_rhs!(G, U[1])
    add_apply_A!(U[1], G, 1.0)
    axpy!(U[2], G, 0.5*dt)
    axpy!(psi, G, 1/6*dt)

    gen_rhs!(G, U[2])
    add_apply_A!(U[2], G, 1.0)
    axpy!(U[3], G, 0.5*dt)
    axpy!(psi, G, 1/3*dt)

    gen_rhs!(G, U[3])
    add_apply_A!(U[3], G, 1.0)
    axpy!(U[4], G, dt)
    axpy!(psi, G, 1/3*dt)

    gen_rhs!(G, U[4])
    add_apply_A!(U[4], G, 1.0)
    axpy!(psi, G, 1/6*dt)
    psi
end

function ERK4_step_lawson!(psi::WaveFunction, dt::Number;
    #Lawson Runge-Kutta method
    G::WaveFunction=wave_function(psi.m),
    U::Vector{WaveFunction}=[wave_function(psi.m) for j=1:4])

    s = 4
    c = [0.0, 0.5, 0.5, 1.0]*dt
    t = get_time(psi)
    copy!(U[1], psi)
    copy!(U[2], U[1])
    propagate_A!(U[2], 0.5*dt)
    copy!(U[3], U[2])
    copy!(U[4], U[3])
    propagate_A!(U[4], 0.5*dt)
    copy!(psi, U[4])
    for j=1:s #seems to be not necessary: propagate_A! already yields right time...
        set_time!(U[j], t+c[j])
    end
    set_time!(psi, t+dt)

    gen_rhs!(G, U[1])
    propagate_A!(G, 0.5*dt)
    axpy!(U[2], G, 0.5*dt)
    propagate_A!(G, 0.5*dt)
    axpy!(psi, G, 1/6*dt)

    gen_rhs!(G, U[2])
    axpy!(U[3], G, 0.5*dt)
    propagate_A!(G, 0.5*dt)
    axpy!(psi, G, 1/3*dt)

    gen_rhs!(G, U[3])
    propagate_A!(G, 0.5*dt)
    axpy!(U[4], G, dt)
    axpy!(psi, G, 1/3*dt)

    gen_rhs!(G, U[4])
    axpy!(psi, G, 1/6*dt)
    psi
end

function ERK4_step_gauss_lawson!(psi::WaveFunction, dt::Number;
    iters::Int=4,
    #Lawson Runge-Kutta method
    G::WaveFunction=wave_function(psi.m),
    U::Vector{WaveFunction}=[wave_function(psi.m) for j=1:3])

    c1 = 1/2-sqrt(3)/6
    c2 = 1/2+sqrt(3)/6
    a11 = 1/4
    a12 = 1/4-sqrt(3)/6
    a21 = 1/4+sqrt(3)/6
    a22 = 1/4
    G1 = G
    G2 = U[3]
    copy!(U[1], psi)
    copy!(U[2], psi)
    propagate_A!(U[1], c1*dt)
    propagate_A!(U[2], c2*dt)

    gen_rhs!(G1, U[1])
    gen_rhs!(G2, U[2])

    for iter = 1:iters #fixed point iteration
        if iter>1
            copy!(U[1], psi)
            copy!(U[2], psi)
            propagate_A!(U[1], c1*dt)
            propagate_A!(U[2], c2*dt)
        end

        axpy!(U[1], G1, a11*dt)
        axpy!(U[2], G2, a22*dt)

        propagate_A!(G1, (c2-c1)*dt)
        propagate_A!(G2, (c1-c2)*dt)

        axpy!(U[1], G2, a12*dt)
        axpy!(U[2], G1, a21*dt)

        gen_rhs!(G1, U[1])
        gen_rhs!(G2, U[2])
    end

    propagate_A!(psi, dt)

    propagate_A!(G1, (1.0-c1)*dt)
    propagate_A!(G2, (1.0-c2)*dt)

    axpy!(psi, G1, 0.5*dt)
    axpy!(psi, G2, 0.5*dt)
end

function ERK4_step_strehmel_weiner!(psi::WaveFunction, dt::Number;
    #Strehmel/Weiner, eq. (2.43) in Hochbruck/Ostermann
    G::WaveFunction=wave_function(psi.m),
    U::Vector{WaveFunction}=[wave_function(psi.m) for j=1:4])

    s = 4
    c = [0.0, 0.5, 0.5, 1.0]*dt
    t = get_time(psi)
    for j=1:s
        copy!(U[j], psi)
        set_time!(U[j], t+c[j])
    end
    set_time!(psi, t+dt)

    gen_rhs!(G, U[1])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[2], c[2], 1, 0.5*dt)
    add_phi_A!(G, U[3], c[3], 1, 0.5*dt)
    add_phi_A!(G, U[3], c[3], 2, -0.5*dt)
    add_phi_A!(G, U[4], c[4], 1, 1.0*dt)
    add_phi_A!(G, U[4], c[4], 2, -2.0*dt)
    add_phi_A!(G, psi, dt, 1, 1.0*dt)
    add_phi_A!(G, psi, dt, 2, -3.0*dt)
    add_phi_A!(G, psi, dt, 3, 4.0*dt)

    gen_rhs!(G, U[2])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[3], c[3], 2, 0.5*dt)
    add_phi_A!(G, U[4], c[4], 2, -2.0*dt)

    gen_rhs!(G, U[3])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, U[4], c[4], 2, 4.0*dt)
    add_phi_A!(G, psi, dt, 2, 4.0*dt)
    add_phi_A!(G, psi, dt, 3, -8.0*dt)

    gen_rhs!(G, U[4])
    add_apply_A!(U[1], G, 1.0)
    add_phi_A!(G, psi, dt, 2, -1.0*dt)
    add_phi_A!(G, psi, dt, 3, 4.0*dt)
    psi
end

function step!(m::ExponentialRungeKutta, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int, step::Int)
    if m.scheme==4
        RK4_step!(psi, dt, G=m.G, U=m.U)
    elseif m.scheme==14
        ERK4_step_lawson!(psi, dt, G=m.G, U=m.U)
    elseif m.scheme==24
        ERK4_step_gauss_lawson!(psi, dt, G=m.G, U=m.U, iters=m.iters)
    elseif m.scheme==39
        ERK2_step!(psi, dt, G=m.G, U=m.U)
    elseif m.scheme==40
        ERK2A_step!(psi, dt, G=m.G, U=m.U)
    elseif m.scheme==42
        ERK4_step_krogstad!(psi, dt, G=m.G, U=m.U)
    elseif m.scheme==43
        ERK4_step_strehmel_weiner!(psi, dt, G=m.G, U=m.U)
    elseif m.scheme==41
        ERK4_step_cox_mathews!(psi, dt, G=m.G, U=m.U)
    end
 end

#######################################
# Exponential multistep
#######################################

struct QuadratureRule
    c::Vector{Float64} # nodes normed to interval [0,1]
    b::Vector{Float64} # weights
end

function gen_interpolation_matrix(x::Vector{Float64}, z::Vector{Float64})
    n = length(x)
    m = length(z)
    L = zeros(n, m)
    for j=1:m
        for i=1:n
            L[i,j] = prod(Float64[(z[j]-x[k])/(x[i]-x[k]) for k=1:i-1])*prod(Float64[(z[j]-x[k])/(x[i]-x[k]) for k=i+1:n])
        end
    end
    L
end

mutable struct ExponentialMultistep <: TimePropagationMethod
    N::Int
    N1::Int
    iters::Int
    ptr::Int
    version::Int
    C1::Union{Matrix{Float64}, Vector{Float64}}
    C2::Union{Matrix{Float64}, Vector{Float64}}
    a::Vector{Float64}
    b::Vector{Float64}
    final_iteration::Bool
    combine_first::Bool
    rhs_back #::Vector{WaveFunction} # storage for previous solution values
    psi0 #::WaveFunction
    acc #::WaveFunction
    starting_method #::TimePropagationMethod
    starting_subdivision::Int # subdivision of dt for calculation of starting values

    function ExponentialMultistep(N1::Int; version::Int=1, iters::Int=0,
        final_iteration::Bool=false, combine_first::Bool=true, starting_method::Union{Nothing,TimePropagationMethod}=nothing,
        quadrature::Union{Nothing,QuadratureRule}=nothing, starting_subdivision::Int=1)
        N = iters==0 ? N1 : N1+1
        @assert version>=0 && version<=2
        a = Float64[]
        b = Float64[]
        if quadrature!=nothing
            version=3
            a = vcat(quadrature.c,1)
            a = vcat(a[1], a[2:end]-a[1:end-1])
            b = vcat(quadrature.b,0)
            C1 = gen_interpolation_matrix(collect((-N1+1.0):0.0), quadrature.c)
            C2 = gen_interpolation_matrix(collect((-N1+1.0):1.0), quadrature.c)
        elseif version==1
            C1 = Matrix{Float64}(inv(Rational{Int}[n^m//factorial(m) for n=-N1+1:0, m=0:N1-1]))
            C2 = Matrix{Float64}(inv(Rational{Int}[n^m//factorial(m) for n=-N1+1:1, m=0:N1]))
        else
            C1 = Vector{Float64}(Rational{Int}[n^m for m=0:N1-1, n=-N1+1:0]\Rational{Int}[1//(m+1) for m=0:N1-1])
            C2 = Vector{Float64}(Rational{Int}[n^m for m=0:N1,   n=-N1+1:1]\Rational{Int}[1//(m+1) for m=0:N1])
        end
        ptr = N
        new(N, N1, iters, ptr, version, C1, C2, a, b, final_iteration, combine_first, nothing, nothing, nothing, starting_method, starting_subdivision)
    end
end

function gen_exponential_multistep_starting_values(psi::WaveFunction, dt::Number, N::Int;
          final_iteration=false, combine_first::Bool=false)
    acc = wave_function(psi.m)
    psi0 = wave_function(psi.m)
    rhs_back = WaveFunction[wave_function(psi.m) for j=1:N]
    gen_rhs!(rhs_back[1], psi)
    t0 = get_time(psi)
    copy!(psi0, psi)
    for K=2:N
        for R=1:(final_iteration&&K==N ? 2 : 1)
            copy!(psi, psi0)
            for J=2:K
                C  = Matrix{Float64}(inv(Rational{Int}[k^m//factorial(m) for k=-J+2:K-J, m=0:K-2]))
                if combine_first
                    copy!(acc, rhs_back[J-1])
                    add_apply_A!(psi, acc, 1.0)
                    add_phi_A!(acc, psi, dt, 1, dt)
                    set_time!(psi, get_time(psi)+dt)
                else
                    propagate_A!(psi, dt)
                end
                for m=(combine_first ? 2 : 1):K-1
                    set!(acc, 0.0)
                    for k=1:K-1
                        axpy!(acc, rhs_back[k], C[m,k])
                    end
                    add_phi_A!(acc, psi, dt, m, dt)
                end
                gen_rhs!(rhs_back[J], psi)
            end
        end
    end
    rhs_back
end

function gen_exponential_multistep2_starting_values(psi::WaveFunction, dt::Number, N::Int;
          final_iteration=false)
    psi0 = wave_function(psi.m)
    rhs_back = WaveFunction[wave_function(psi.m) for j=1:N]
    gen_rhs!(rhs_back[1], psi)
    t0 = get_time(psi)
    copy!(psi0, psi)
    for K=2:N
        for R=1:(final_iteration && K==N ? 2 : 1)
            copy!(psi, psi0)
            for k=1:K-1
                propagate_A!(rhs_back[k], -(K-2)*dt)
            end
            for J=2:K
                C = Vector{Float64}(Rational{Int}[n^m for m=0:K-2, n=-J+2:K-J]\Rational{Int}[1//(m+1) for m=0:K-2])
                for k=1:K-1
                   axpy!(psi, rhs_back[k], C[k]*dt)
                end
                propagate_A!(psi, dt)
                gen_rhs!(rhs_back[J], psi)
                propagate_A!(rhs_back[J], -dt)
                for k=1:K
                    propagate_A!(rhs_back[k], dt)
                end
            end
        end
    end
    rhs_back
end

function gen_extrapolated_wf!(psi::WaveFunction, L::Vector{Float64}, psi_back::Vector{WaveFunction}, ptr::Int;
                              order::Int=length(psi_back))
    n = length(psi_back)
    set!(psi, 0)
    for j=1:order
        k = mod(j-order+ptr-1, n)+1
        axpy!(psi, psi_back[k], L[j])
    end
end

function gen_exponential_multistep3_starting_values(psi::WaveFunction, dt::Number, N::Int,
    a::Vector{Float64}, b::Vector{Float64};final_iteration=false)
    rhs = wave_function(psi.m)
    s = length(a)
    aa = cumsum(a)
    psi0 = wave_function(psi.m)
    rhs_back = WaveFunction[wave_function(psi.m) for j=1:N]
    copy!(psi0, psi)
    gen_rhs!(rhs_back[1], psi)
    t0 = get_time(psi)
    for K=2:N
        for R=1:(final_iteration && K==N ? 2 : 1)
            copy!(psi, psi0)
            for J=2:K
                L = gen_interpolation_matrix(collect(0.0:(K-2.0)), J-2.0+aa) # TODO: check order->order+/-1
                for j = 1:s
                    if a[j]!=0.0
                        propagate_A!(psi, a[j]*dt)
                    end
                    if b[j]!=0.0
                        gen_extrapolated_wf!(rhs, L[:,j], rhs_back[1:K-1], 2*K-2) # TODO: check order->order+/-1
                        axpy!(psi, rhs, b[j]*dt)
                    end
                end
                gen_rhs!(rhs_back[J], psi)
            end
        end
    end
    rhs_back
end

function initialize!(method::ExponentialMultistep, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int)
    if method.starting_method!=nothing
        method.rhs_back = WaveFunction[wave_function(psi.m) for j=1:method.N]
        gen_rhs!(method.rhs_back[1], psi)
        k=1
        for I in EquidistantTimeStepper(method.starting_method, psi, t0, dt/method.starting_subdivision, (method.N1-1)*method.starting_subdivision)
            if I[1] % method.starting_subdivision == 0
                k += 1
                gen_rhs!(method.rhs_back[k], psi)
            end
        end
        if method.version==2
            for k=0:method.N1-1
                propagate_A!(method.rhs_back[method.N1-k], dt*k)
            end
        end
    else
        if method.version==1
            method.rhs_back = gen_exponential_multistep_starting_values(psi, dt, method.N1,
                           final_iteration=method.final_iteration, combine_first=method.combine_first)
        elseif method.version==2
            method.rhs_back = gen_exponential_multistep2_starting_values(psi, dt, method.N1,
                           final_iteration=method.final_iteration)
        elseif method.version==3
            method.rhs_back = gen_exponential_multistep3_starting_values(psi, dt, method.N1, method.a, method.b,
                           final_iteration=method.final_iteration)
        end
        if method.iters>=1
            push!(method.rhs_back, wave_function(psi.m))
        end
    end
    if method.iters>=1
        method.psi0 = wave_function(psi.m)
    end
    if method.version==1 || method.version==3
        method.acc = wave_function(psi.m)
    end
    method.ptr = method.N1
    method.N1-1
end

function finalize!(method::ExponentialMultistep, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int, step::Int)
    method.psi0 = nothing
    method.acc = nothing
    method.rhs_back = nothing
end

function step!(m::ExponentialMultistep, psi::WaveFunction,
         t0::Real, dt::Real, steps::Int, step::Int)

    if m.iters>=1
        copy!(m.psi0, psi)
    end

    #predictor
    if m.version==1
        if m.combine_first
            k1 = mod(m.ptr-1, m.N)+1
            copy!(m.acc, m.rhs_back[k1])
            add_apply_A!(psi, m.acc, 1.0)
            add_phi_A!(m.acc, psi, dt, 1, dt)
            set_time!(psi, get_time(psi)+dt)
        else
            propagate_A!(psi, dt)
        end
        for n=(m.combine_first ? 2 : 1):m.N1
            set!(m.acc, 0.0)
            for k=1:m.N1
                k1 = mod(k-m.N1+m.ptr-1, m.N)+1
                axpy!(m.acc, m.rhs_back[k1], m.C1[n,k])
            end
            add_phi_A!(m.acc, psi, dt, n, dt)
        end
    elseif m.version==3
        for j = 1:length(m.a)
            if m.a[j]!=0.0
                propagate_A!(psi, m.a[j]*dt)
            end
            if m.b[j]!=0.0
                gen_extrapolated_wf!(m.acc, m.C1[:,j], m.rhs_back, m.ptr, order=m.N1)
                axpy!(psi, m.acc, m.b[j]*dt)
            end
        end
    elseif m.version==2
        for k=1:m.N1
            k1 = mod(k-m.N1+m.ptr-1, m.N)+1
            axpy!(psi, m.rhs_back[k1], m.C1[k]*dt)
        end
        propagate_A!(psi, dt)
    end

    m.ptr = mod(m.ptr, m.N) + 1
    gen_rhs!(m.rhs_back[m.ptr], psi)
    if m.version==2 && m.iters>=1
        propagate_A!(m.rhs_back[m.ptr], -dt)
    end

    for iter=1:m.iters
        copy!(psi, m.psi0)

        #corrector
        if m.version==1
            if m.combine_first
                k1 = mod(m.ptr-2, m.N)+1
                copy!(m.acc, m.rhs_back[k1])
                add_apply_A!(psi, m.acc, 1.0)
                add_phi_A!(m.acc, psi, dt, 1, dt)
                set_time!(psi, get_time(psi)+dt)
            else
                propagate_A!(psi, dt)
            end
            for n=(m.combine_first ? 2 : 1):m.N1+1
                set!(m.acc, 0.0)
                for k=1:m.N1+1
                    k1 = mod(k-m.N1+m.ptr-2, m.N)+1
                    axpy!(m.acc, m.rhs_back[k1], m.C2[n,k])
                end
                add_phi_A!(m.acc, psi, dt, n, dt)
            end
        elseif m.version==3
            for j = 1:length(m.a)
                if m.a[j]!=0.0
                    propagate_A!(psi, m.a[j]*dt)
                end
                if m.b[j]!=0.0
                    gen_extrapolated_wf!(m.acc, m.C2[:,j], m.rhs_back, m.ptr, order=m.N1+1)
                    axpy!(psi, m.acc, m.b[j]*dt)
                end
            end
        elseif m.version==2
            for k=1:m.N1+1
                k1 = mod(k-m.N1+m.ptr-2, m.N)+1
                axpy!(psi, m.rhs_back[k1], m.C2[k]*dt)
            end
            propagate_A!(psi, dt)
        end
        gen_rhs!(m.rhs_back[m.ptr], psi)
    end

    if m.version==2
        for k=1:m.N
            if k==m.ptr continue end
            propagate_A!(m.rhs_back[k], dt)
        end
    end
end

###############################################################################
# Adaptive Propagators
###############################################################################

abstract type AdaptiveTimePropagationMethod end

mutable struct AdaptiveTimeStepper
    method::AdaptiveTimePropagationMethod

    psi::WaveFunction
    t0::Real
    tend::Real
    tol::Real
    dt::Real
end

#default initializer
function initialize!(method::AdaptiveTimePropagationMethod, psi::WaveFunction,
         t0::Real, tend::Real, tol::Real, dt::Real)
    t0 # iteration  starts with step 0
end

#default finalizer
function finalize!(method::AdaptiveTimePropagationMethod, psi::WaveFunction,
         t0::Real, tend::Real, tol::Real, dt::Real, t::Real)
end

function Base.iterate(tsi::AdaptiveTimeStepper, t::Real=0)
    # TODO hack, tend <= 0 cannot be used
    if t == 0
        set_propagate_time_together_with_A!(tsi.psi.m, true)
        set_time!(tsi.psi, tsi.t0)
        step = initialize!(tsi.method, tsi.psi, tsi.t0, tsi.tend, tsi.tol, tsi.dt)
    end

    if t >= tsi.tend
        finalize!(tsi.method, tsi.psi, tsi.t0, tsi.tend, tsi.tol, tsi.dt, t)
        return nothing
    end

    (tnew, dtnew) = step!(tsi.method, tsi.psi, tsi.t0, tsi.tend, tsi.tol, tsi.dt, t)
    tsi.dt = dtnew
    return (tnew, tsi), tnew
end

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

function solve_vander_trans(x::Vector, b::Vector)
    # Algorithm 4.6.2 from Golub/van Loan
    n = length(x)
    for k=1:n-1
        for i=n:-1:k+1
            b[i] -= x[k]*b[i-1]
        end
    end
    for k=n-1:-1:1
        for i=k+1:n
            b[i] /= (x[i]-x[i-k])
        end
        for i=k:n-1
            b[i] -= b[i+1]
        end
    end
end

#c = Float64[dt^m/(m+1) for m=0:N1-1]
#solve_vander_trans(t_back, c)
#c

mutable struct AdaptiveAdamsLawson <: AdaptiveTimePropagationMethod
    N::Int
    N1::Int
    ptr::Int
    t_back::Vector{Float64}
    rhs_back #::Vector{WaveFunction} # storage for previous solution values
    psi0 #::WaveFunction
    psi1 #::WaveFunction
    acc  #::WaveFunction
    starting_method #::TimePropagationMethod
    bootstrap_mode::Bool
    N1_final::Int

    combine_first::Bool
    version::Int

    function AdaptiveAdamsLawson(N1::Int; starting_method::Union{Nothing,TimePropagationMethod}=nothing,
                                combine_first::Bool=true, version::Int=2)
        new(N1+1, N1, N1, Float64[], nothing, nothing, nothing, nothing, starting_method, starting_method==nothing, N1, combine_first, version)
    end
end

function initialize!(m::AdaptiveAdamsLawson, psi::WaveFunction,
         t0::Real, tend::Real, tol::Real, dt::Real)
    m.bootstrap_mode = m.starting_method==nothing
    N_final = m.N1_final+1
    m.t_back = zeros(Float64, N_final)
    m.rhs_back = WaveFunction[wave_function(psi.m) for j=1:N_final]
    m.psi0 = wave_function(psi.m)
    m.psi1 = wave_function(psi.m)
    if m.version==1
        m.acc = wave_function(psi.m)
    end
    gen_rhs!(m.rhs_back[1], psi)
    if m.bootstrap_mode
        m.N = 2
        m.N1 = 1
    else
        k=1
        for I in EquidistantTimeStepper(m.starting_method, psi, t0, dt, m.N1-1)
            k += 1
            gen_rhs!(m.rhs_back[k], psi)
        end
        m.t_back[:] = Float64[dt*k for k=0:m.N1]
        if m.version==2 #Lawson
            for k=0:m.N1-1
                propagate_A!(m.rhs_back[m.N1-k], dt*k)
            end
        end
    end
    m.ptr = m.N1
    (m.N1-1)*dt
end

function finalize!(method::AdaptiveAdamsLawson, psi::WaveFunction,
         t0::Real, tend::Real, tol::Real, dt::Real, t::Real)
    method.psi0 = nothing
    method.psi1 = nothing
    method.acc = nothing
    method.rhs_back = nothing
end

function step!(m::AdaptiveAdamsLawson, psi::WaveFunction,
         t0::Real, tend::Real, tol::Real, dt::Real, t::Real)
    facmin = 0.25
    facmax = 4.0
    fac = 0.9

    copy!(m.psi0, psi)
    copy!(m.psi1, psi)

    ptr0 = m.ptr
    dt0 = dt

    err = 2.0 #error/t
    while err>=1.0
        dt = min(dt, tend-t)
        dt0 = dt

        #predictor
        tt = Float64[m.t_back[mod(k-m.N1+m.ptr-1, m.N)+1] for k=1:m.N1]
        tt = tt .- tt[end]
        if m.version==2 # Lawson
            C = Float64[dt^m/(m+1) for m=0:m.N1-1]
            solve_vander_trans(tt, C)
            for k=1:m.N1
                k1 = mod(k-m.N1+m.ptr-1, m.N)+1
                axpy!(m.psi1, m.rhs_back[k1], C[k]*dt)
            end
            propagate_A!(m.psi1, dt)
        elseif m.version==1
            C = inv([(t/dt)^m/factorial(m) for t in tt, m=0:m.N1-1])
            if m.combine_first
                k1 = mod(m.ptr-1, m.N)+1
                copy!(m.acc, m.rhs_back[k1])
                add_apply_A!(m.psi1, m.acc, 1.0)
                add_phi_A!(m.acc, m.psi1, dt, 1, dt)
                set_time!(m.psi1, get_time(m.psi1)+dt)
            else
                propagate_A!(m.psi1, dt)
            end
            for n=(m.combine_first ? 2 : 1):m.N1
                set!(m.acc, 0.0)
                for k=1:m.N1
                    k1 = mod(k-m.N1+m.ptr-1, m.N)+1
                    axpy!(m.acc, m.rhs_back[k1], C[n,k])
                end
                add_phi_A!(m.acc, m.psi1, dt, n, dt)
            end
        end

        if m.bootstrap_mode
            m.ptr += 1
        else
            m.ptr = mod(m.ptr, m.N) + 1
        end
        gen_rhs!(m.rhs_back[m.ptr], m.psi1)
        m.t_back[m.ptr] = t+dt
        if m.version==2
            propagate_A!(m.rhs_back[m.ptr], -dt)
        end

        #corrector
        tt = Float64[m.t_back[mod(k-m.N1+m.ptr-2, m.N)+1] for k=1:m.N1+1]
        tt = tt .- tt[end-1]
        if m.version==2 # Lawson
            C = Float64[dt^m/(m+1) for m=0:m.N1]
            solve_vander_trans(tt, C)
            for k=1:m.N1+1
                k1 = mod(k-m.N1+m.ptr-2, m.N)+1
                axpy!(psi, m.rhs_back[k1], C[k]*dt)
            end
            propagate_A!(psi, dt)
        elseif m.version==1
            C = inv([(t/dt)^m/factorial(m) for t in tt, m=0:m.N1])
            if m.combine_first
                k1 = mod(m.ptr-2, m.N)+1
                copy!(m.acc, m.rhs_back[k1])
                add_apply_A!(psi, m.acc, 1.0)
                add_phi_A!(m.acc, psi, dt, 1, dt)
                set_time!(psi, get_time(psi)+dt)
            else
                propagate_A!(psi, dt)
            end
            for n=(m.combine_first ? 2 : 1):m.N1+1
                set!(m.acc, 0.0)
                for k=1:m.N1+1
                    k1 = mod(k-m.N1+m.ptr-2, m.N)+1
                    axpy!(m.acc, m.rhs_back[k1], C[n,k])
                end
                add_phi_A!(m.acc, psi, dt, n, dt)
            end

        end

        err = distance(psi, m.psi1)/tol
        dt = dt*min(facmax, max(facmin, fac*(1.0/err)^(1.0/(m.N1+1)))) # CHECK m.N+1=p+1 !!!

        if err>=1.0
            copy!(m.psi1, m.psi0)
            copy!(psi, m.psi0)
            m.ptr = ptr0
            @Printf.printf("t=%17.9e  err=%17.8e  dt=%17.8e  rejected...\n",
                    Float64(t), Float64(err), Float64(dt))
        end
    end

    gen_rhs!(m.rhs_back[m.ptr], psi)
    t += dt0
    m.t_back[m.ptr] = t

    if m.version==2
        for k=1:m.N
            if k==m.ptr continue end
            propagate_A!(m.rhs_back[k], dt0)
        end
    end
    if m.bootstrap_mode
        if m.N1_final>m.N1
            m.N1 += 1
            m.N += 1
        else
            m.bootstrap_mode = false
        end
    end
    (t, dt)
end


text:
using DataFrames
using LinearAlgebra
using MixedModels
using MixedModelsSim
using PooledArrays
using Tables
using Test

@testset "factorproduct" begin
	items = (item = nlevels(50, 'I'), category = repeat(["ingroup", "outgroup"], inner=25))
	df = columntable(factorproduct(items, (subj = nlevels(100),)))
	@test length(df) == 3
	@test length(first(df)) == 5000
end

@testset "levels" begin
    levs = nlevels(100, 'k')
    @test length(levs) == 100
    @test length(unique(levs)) == 100
    @test first(levs) == "k001"
end

@testset "nlevstbl" begin
    subj = DataFrame(nlevstbl(:subj, 20))
    @test nrow(subj) == 20
    @test isone(ncol(subj))
    @test propertynames(subj) == [:subj]
    item = DataFrame(nlevstbl(:item, 9, :lev => ["low", "medium", "high"]))
    @test nrow(item) == 9
    @test ncol(item) == 2
    @test propertynames(item) == [:item, :lev]
    @test first(first(subj)) == "S01"
    df = crossjoin(subj, item)
    @test nrow(df) == 180
    @test ncol(df) == 3
    @test_throws ArgumentError nlevstbl(:baditem, 9, :nlev => ["low", "high"])
    twofac = nlevstbl(:item, 18, :level => ["low", "medium", "high"], :blur => ["N", "Y"])
    @test length(twofac) == 3
    @test length(first(twofac)) == 18
    @test_throws ArgumentError nlevstbl(:item, 18, :level => ["low", "high"], :blur => ["N", "Y"])

    item = nlevstbl(:subj, 128, :trt => ["C", "T"])
    @test item.subj isa PooledVector{String, Int16, Vector{Int16}}
    @test item.trt isa PooledVector{String, Int8, Vector{Int8}}
    @test length(item.subj) == length(item.trt)
end

@testset "pooled!" begin
    subject = (subj = ["S1","S2","S3","S4","S5"], age=["O","O","O","Y","Y"]);

    @test_throws ArgumentError pooled!(subject)

    df = pooled!(DataFrame(subject))

    @test df.subj isa PooledArray
    @test df.age isa PooledArray
end

@testset "flatlowertri" begin
    l = LowerTriangular([1 0 0; 2 3 0; 4 5 6])
    @test all(flatlowertri(l) .== [1, 2, 4, 3, 5, 6])
end

@testset "update!" begin
    @testset "LMM" begin
        fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)),
                  MixedModels.dataset(:sleepstudy); progress=false)
        update!(fm1; θ=[1.0, 0.0, 1.0])

        @test all(values(first(fm1.σs)) .== fm1.σ)
        @test only(VarCorr(fm1).σρ.subj.ρ) == 0.0
    end

    @testset "GLMM" begin
        cbpp = MixedModels.dataset(:cbpp)
        gm1 = fit(MixedModel, @formula((incid/hsz) ~ 1 + period + (1|herd)), cbpp, Binomial();
                wts=cbpp.hsz, fast=true, progress=false)

        β = repeat([-1.], 4)
        θ = [0.5]

        β₀ = copy(fixef(gm1))
        update!(gm1; θ=θ)
        @test all(values(first(gm1.σs)) .== θ)
        @test all(gm1.β .== β₀)

        refit!(gm1, fast=false, progress=false)
        β₀ = copy(fixef(gm1))
        update!(gm1; θ=θ)
        @test all(values(first(gm1.σs)) .== θ)
        @test all(gm1.β .== β₀)
    end

end

@testset "create_re|θ" begin
    fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)),
              MixedModels.dataset(:sleepstudy); progress=false)

    σs = values(first(fm1.σs)) ./ fm1.σ
    ρ = only(first(VarCorr(fm1).σρ).ρ)
    corr = [1 ρ; ρ 1]

    subjre = create_re(σs...; corrmat=corr)

    @test all(only(fm1.λ) .≈ subjre)
    @test create_re(σs...) == diagm([σs...])

    fm1unfitted = LinearMixedModel(@formula(reaction ~ 1 + days + (1 + days|subj)),
                                   MixedModels.dataset(:sleepstudy))

    @test_logs((:warn,
                "Specifying the random effects by position instead of name is deprecated"),
                update!(fm1unfitted, subjre))

    update!(fm1unfitted, subjre)

    vcu, vc = VarCorr(fm1unfitted), VarCorr(fm1)

    @test all(values(vcu.σρ.subj.σ) .≈ values(vc.σρ.subj.σ))
    @test all(values(vcu.σρ.subj.ρ) .≈ values(vc.σρ.subj.ρ))

    @test_throws ArgumentError update!(fm1unfitted; subjitem=subjre, θ=fm1.θ)
    @test_throws ArgumentError update!(fm1unfitted)
    @test_throws ArgumentError update!(fm1unfitted; item=subjre)

    fmkb = fit(MixedModel,
               @formula(rt_trunc ~ 1 + spkr + prec + load + (1|subj) + (1 + prec|item)),
               MixedModels.dataset(:kb07); progress=false)

    fmkb2 = LinearMixedModel(@formula(rt_trunc ~ 1 + spkr + prec + load + (1|subj) + (1 + prec|item)),
                            MixedModels.dataset(:kb07))
    update!(fmkb2; subj=fmkb.reterms[2].λ, item=fmkb.reterms[1].λ)

    vcu, vc = VarCorr(fmkb2), VarCorr(fmkb2)

    @test all(values(vcu.σρ.subj.σ) .≈ values(vc.σρ.subj.σ))
    @test all(values(vcu.σρ.subj.ρ) .≈ values(vc.σρ.subj.ρ))
    @test all(values(vcu.σρ.item.σ) .≈ values(vc.σρ.item.σ))
    @test all(values(vcu.σρ.item.ρ) .≈ values(vc.σρ.item.ρ))

end

ext: jl


text:
struct Nuclide
  charge::Float64
  mass::Float64
  charge_over_mass::Float64
end
Nuclide(q, m) = Nuclide(q, m, q / m)

charge(n::Nuclide) = n.charge
mass(n::Nuclide) = n.mass
charge_mass_ratio(n::Nuclide) = n.charge_over_mass

abstract type AbstractParticle end

mutable struct Particle{S<:AbstractShape} <: AbstractParticle
  nuclide::Nuclide
  basis::BasisFunction{S}
  velocity::Float64
end
function Particle(n::Nuclide, shape::AbstractShape=DeltaFunctionShape();
    kwargstuple...) where {S<:AbstractShape}
  kwargs = Dict(kwargstuple)
  haskey(kwargs, :x) && (kwargs[:position] = kwargs[:x])
  haskey(kwargs, :v) && (kwargs[:velocity] = kwargs[:v])
  haskey(kwargs, :w) && (kwargs[:weight] = kwargs[:w])
  position = get(kwargs, :position, 0.0)
  velocity = get(kwargs, :velocity, 0.0)
  weight = get(kwargs, :weight, 0.0)
  basis = BasisFunction(shape, position, weight)
  return Particle(n, basis, velocity)
end

velocity(p::Particle) = p.velocity
basis(p::Particle) = p.basis
Base.position(p::Particle) = centre(p)

for op ∈ (:weight, :lower, :upper, :width, :centre)
  @eval $op(p::Particle) = $op(basis(p))
end
for op ∈ (:mass, :charge, :charge_mass_ratio)
  @eval $op(p::Particle) = $op(p.nuclide)
end
for op ∈ (:-, :+)
  @eval function Base.$op(a::Particle{S}, b::Particle{S}) where {S<:AbstractShape}
    @assert a.nuclude == b.nuclude
    return Particle(a.nuclude, $op(a.basis, b.basis), $op(a.velocity, b.velocity))
  end
end

for op ∈ (:*, :/)
  @eval function Base.$op(a::Particle, x::Number)
    return Particle(a.nuclude, $op(a.basis, x), $op(a.velocity, x))
  end
end

function pushposition!(p::Particle, dt, bc::AbstractBC)
  translate!(p.basis, p.velocity * dt, bc)
  return p
end
function pushvelocity!(p::Particle, E::Number, dt)
  p.velocity += charge_mass_ratio(p) * E * dt
  return p
end
import Base.push!
Base.push!(p::Particle, E, dt, bc) = (pushposition!(p, dt, bc); pushvelocity!(p, E, dt); p)

function pushvelocity!(p::Particle, f::AbstractField, dt)
  E = electricfield(f, p)
  @assert isfinite(E)
  return pushvelocity!(p, E, dt)
end


text:
function  ueberpruefeGeradeoderPara(PX1,PX2,theta,Geradengleichung,a,b,ParameterdarstellungX1,ParameterdarstellungX2)

t = 0.6221;
Z = 3;
S = 1;
x = 0.7291;
Wahl = [];

if isnan(ParameterdarstellungX1(t,Z,S)) && isnan(ParameterdarstellungX1(t,Z,S)) && !isnan(Geradengleichung(Z,S,x))
    if Geradengleichung(Z,S,x) == tand(theta[S])*x+PX2[Z,S]-tand(theta[S])*PX1[Z,S]
        println("Ihr habt die Geradengleichung korrekt aufgestellt!")
        Wahl = 5.4; # 5.4 steht für Gerade
    elseif Geradengleichung(Z,S,x) == cosd(theta[S])*x+PX2[Z,S]-cosd(theta[S])*PX1[Z,S]
        println("Schaut euch noch einmal die Definitionen von Kosinus, Sinus und Tangens an und überprüft eure Gleichung.")
        Wahl = [];
    elseif Geradengleichung(Z,S,x) == sind(theta(S))*x+PX2(Z,S)-sind(theta(S))*PX1(Z,S)
        println("Schaut euch noch einmal die Definitionen von Kosinus, Sinus und Tangens an und überprüft eure Gleichung.")
        Wahl = [];
    else
        println("Eure Gleichung ist noch nicht korrekt.")
        Wahl = [];
    end
elseif !isnan(ParameterdarstellungX1(t,Z,S)) && !isnan(ParameterdarstellungX1(t,Z,S)) && isnan(Geradengleichung(Z,S,x))
   if ParameterdarstellungX1(t,Z,S) == PX1[Z,S]+t*cosd(theta[S])
        if ParameterdarstellungX2(t,Z,S) == PX2[Z,S]+t*sind(theta[S])
            println("Eure Parameterdarstellung für die Strahlen ist korrekt!");
            Wahl = 3.7; # 1 steht für Parameterdarstellung
        elseif ParameterdarstellungX2(t,Z,S) == PX2[Z,S]+t*cosd(theta[S])
            println("Der Teil für die X1 Koordinate ist korrekt. Schaut euch für die X2 Koordinaten noch einmal an, wie der Kosinus und Sinus definiert sind.")
            Wahl = [];
        else
            println("Der Teil für die X1 Koordinate ist korrekt. Überprüft erneut eure Formel für die X2 Koordinate.")
            Wahl = [];
        end
    elseif ParameterdarstellungX1(t,Z,S) == PX1[Z,S]+t*sind(theta[S])
        if ParameterdarstellungX2(t,Z,S) == PX2[Z,S]+t*sind(theta[S])
            println("Der Teil für die X2 Koordinate ist korrekt. Schaut euch für die X1 Koordinaten noch einmal an, wie der Kosinus und Sinus definiert sind.")
            Wahl = [];
        else
            println("Schaut euch erneut an, wie der Sinus und Kosinus definiert sind und überprüft erneut eure Formeln.")
            Wahl = [];
        end
    else
        println("Schaut euch erneut an, wie der Sinus und Kosinus definiert sind und überprüft erneut eure Formeln.")
        Wahl = [];
   end
else
    println("Ihr müsst euch für eine der beiden Darstellungsarten entscheiden.")
    Wahl = [];
end
return Wahl;
end


text:
#The LLRBTree data structure, that allows us to acces keys in log2 time.
#Translated from code by Robert Sedgewick
module LLRBTrees

import Base: haskey, getindex, setindex!, delete!, minimum, maximum, push!, show, replace

export TreeNode, LLRBTree, getnodefrom_key, getmaxdepth, orderedpairs, levelorder,
         populate!, isleftleaf, isrightleaf, ischildleaf, isleaf, isred

type TreeNode{K,V}
    key:: Nullable{K}
    value:: Nullable{V}
    isRed:: Bool
    left:: TreeNode{K,V}
    right:: TreeNode{K,V}

    TreeNode() = method_exists(isless, (K,K) ) ? new(Nullable{K}(), Nullable{V}(), false) : error("keys must be comparable")
    TreeNode(key::K, value::V, isRed::Bool, left::TreeNode{K,V}, right::TreeNode{K,V}) =
    method_exists(isless, (K,K) ) ? new(key, value, isRed, left, right): error("keys must be comparable")
    TreeNode(key::K, value::V, isRed::Bool) = method_exists(isless, (K,K) ) ? new(key, value, isRed) : error("keys must be comparable")

end
# New nodes are always red so no method that explicitly assigns the color is neccessary
TreeNode{K,V}(key::K,value::V,left::TreeNode{K,V},right::TreeNode{K,V}) = TreeNode(key,value,true,left,right)
TreeNode{K,V}(key::K,value::V) = TreeNode{K,V}(key, value, true)
TreeNode{K,V}(key::K,value::V, isRed::Bool) = TreeNode{K,V}(key, value, isRed)

# To be able to have different instances of LLRB trees, we create a tree type
# that characterizes a tree by its root
type LLRBTree{K,V}
    root:: TreeNode{K,V}

    LLRBTree()=new()
    LLRBTree(node:: TreeNode{K,V}) = (tree = new(); tree.root=node; tree.root.isRed=false; tree)
    LLRBTree(key::K, value::V) = (tree=new(TreeNode(key,value)); tree.root.isRed=false; tree)
end

LLRBTree{K,V}(key::K,value::V) = LLRBTree{K,V}(key, value)
LLRBTree{K,V}(node::TreeNode{K,V}) = LLRBTree{K,V}(node)

function isleaf{K,V}(node::TreeNode{K,V})
    isnull( node.key )
end

function ischildleaf{K,V}(node::TreeNode{K,V}, left::Bool)
    if(left)
        if(!isdefined(node, :left))
            return true
        else
            return isnull(node.left.key)
        end
    else
        if(!isdefined(node, :right) )
            return true
        else
            return isnull(node.right.key)
        end
    end
end

isleftleaf{K,V}(node::TreeNode{K,V})= ischildleaf(node, true)
isrightleaf{K,V}(node::TreeNode{K,V}) = ischildleaf(node, false)

function isred{K,V}(node::TreeNode{K,V}, dirs::Symbol...)

    steps = length(dirs)
    isred = node.isRed
    foundleaf = 0
    foundundef = 0

    for i in 1:steps
        if(dirs[i]!= :left && dirs[i]!= :right)
            Base.error("Directions must be :left or :right")
        end

    end

    for i in 1:steps

        dirbool = dirs[i]== :left ? true : false

        if(isdefined(node, dirs[i]))
            node = dirbool ? node.left : node.right
            isred = node.isRed

            foundleaf = isnull(node.key) ? i : foundleaf
        else
            isred = false
            foundundef = i
            break
        end
    end

    overflow = (foundundef < steps && foundundef!=0)
    overflow |= (foundleaf < steps && foundleaf!=0)

    return  overflow ? Base.error("Directions move beyond leaves") : isred

end

function getmaxdepth{K,V}(node::TreeNode{K,V}, depth=1)
      #We do a caching up of the number of recursions with the use of parameters and returns
      depthleft=depthright=depth
    if !isleftleaf(node)
          #("left in at $depth")
          depthleft = getmaxdepth(node.left, depth+1)
      end
    if !isrightleaf(node)
          #("right in at $depth")
          depthright = getmaxdepth(node.right, depth+1)
      end
      #("out at $depth")
      #([depthleft,depthright])
      return maximum([depthleft,depthright])
  end
  getmaxdepth(tree::LLRBTree) = getmaxdepth(tree.root)

function replace(original::ASCIIString, start::Int64, fin::Int64, string::ASCIIString)
    left = original[1:start]
    right = original[fin:end]
    return left*string*right
end

function show{K,V}(io::IO, tree::LLRBTree{K,V})

    if(isdefined(tree, :root))
        totalheight = getmaxdepth(tree)
        height = getmaxdepth(tree) < 5 ? totalheight : 5
        width = 2^(height +2)

        len = width * height * 2 + 2 |> Int
        stringtree = ""
        for i = 1:len
            stringtree*= i < len - 2 && i % width == 0 ? "\n" : " "
        end

        if(!isnull(tree.root.key))
            stringtree = showR(stringtree, Int(width / 2), 1, Int(width / 4), width, tree.root, " ")
        else
            stringtree = typeof(tree)
        end
    else
            stringtree = typeof(tree)
    end

    println(io, stringtree)
end

function showR(stringtree::ASCIIString, center::Int, row::Int, indent::Int, width::Int, node::TreeNode,
            edge::ASCIIString)

    if(indent%2 == 0)

        if (!isleftleaf(node))
            stringtree = showR(stringtree, center - indent, row + 2, Int(indent/2), width, node.left, " /");
        end

        valuestring = string(node.value.value)
        valuestring = node.isRed ? valuestring*"_r" : valuestring

        idx1 = Int( row * width + center - div( (length(valuestring) + 1) , 2 ) )
        idx2 = idx1 + length(valuestring) |> Int
        idx3 = idx1 - width |> Int
        if (idx2 < length(stringtree))
            stringtree = replace(stringtree, idx1, idx2 +1, valuestring)
            stringtree = replace(stringtree, idx3, idx3 + 3, edge)
        end

        if (!isrightleaf(node))
            stringtree = showR(stringtree, center + indent, row + 2, Int(indent/2), width, node.right, "\\ ")
        end
    end

    return stringtree
end

# The function that returns the node that corresponds to a certain key, below a root-node
# it returns a TreeLeaf when the root node is empty
function getnodefrom_key{K,V}(part::TreeNode{K,V}, key::K)
  # The recursive search for the key
    if get(part.key) == key
        return part
    elseif key < get(part.key)
        if( !isleftleaf(part) )
            getnodefrom_key(part.left, key)
        else
            return TreeNode{K,V}()
        end
    else
        if(!isrightleaf(part))
            getnodefrom_key(part.right, key)
        else
          return TreeNode{K,V}()
        end
    end
end
getnodefrom_key{K,V}(tree::LLRBTree{K,V}, key::K) = getnodefrom_key(tree.root, key)

haskey{K,V}(tree::LLRBTree{K,V}, key::K) =  !isleaf( getnodefrom_key(tree.root, key) )
haskey{K,V}(node::TreeNode{K,V}, key::K) = !isleaf( getnodefrom_key(node, key) )

# A getindex only that the indexes are keys and we get the value of the first node with the given key
function getindex{K,V}(node::TreeNode{K,V}, key::K)
    part::TreeNode{K,V} = getnodefrom_key(node,key)
    # Throw a KeyException if the key is node.left present

    if isleaf( part )
        throw(KeyError(key))
    else
        return part.value
    end
end

getindex{K,V}(tree::LLRBTree{K,V}, key::K) = getindex(tree.root, key)

function setindex!{K,V}(tree::LLRBTree{K,V}, value::V, key::K)
    if(isdefined(tree, :root))
        part::TreeNode{K,V} = getnodefrom_key(tree.root,key)
    else
        part = TreeNode{K,V}()
    end

    if isleaf( part )
        push!(tree, key, value)
    else
        part.value= value
    end

    tree
end

function rotateleft{K,V}(node::TreeNode{K,V})

    if isrightleaf(node)
        error("The right part of the given node must be a node")
    end
    son=node.right

    if !isleftleaf(son)
        node.right=son.left
    else
        node.right=TreeNode{K,V}()
    end
    son.left=node

    son.isRed=node.isRed
    node.isRed=true
    return son
end

function rotateright{K,V}(node::TreeNode{K,V})

    if isleftleaf(node)
        error("The left part of the given node must be a node")
    end
    son=node.left

    if !isrightleaf(son)
        node.left=son.right
    else
        node.left=TreeNode{K,V}()
    end
    son.right=node

    son.isRed=node.isRed
    node.isRed=true
    return son
end

function flipcolor!{K,V}(node::TreeNode{K,V})

    node.isRed = !node.isRed
    if !isleftleaf(node)
        node.left.isRed = !node.left.isRed
    end
    if !isrightleaf(node)
        node.right.isRed = !node.right.isRed
    end

end

function add_node{K,V}(node::TreeNode{K,V}, key::K, value::V)

    if isleaf(node)
        #("Add the node only when a leaf is found  ",tree)
        return TreeNode(key, value)
    end

    #put flipcolor here to get 2-3-4 trees

    if get(node.key) < key
        #("went right")
        node.right = !isrightleaf(node) ? node.right: TreeNode{K,V}()
        node.right = add_node(node.right, key, value)
    elseif get(node.key) > key
        #("went left")
        node.left = !isleftleaf(node) ? node.left: TreeNode{K,V}()
        node.left = add_node(node.left, key, value)
    else
        #("node was same")
        node.value = value
    end

    if ( isred(node, :right) && !isred(node, :left))
        #("Rotate left to prevent red node on right ", node)
        node=rotateleft(node)
    end

    if ( isred(node, :left) && isred(node, :left, :left) )
        #("Rotate right to prevent consecutive red nodes ", node)
        node=rotateright(node)
    end
    if isred(node, :left) && isred(node, :right)
        #("Split node with two red children")
        flipcolor!(node)
    end

    return node
end
# For trees and nodes as parameters; different name so the other one doesn't get used
push!{K,V}(tree::LLRBTree{K,V}, key::K, value::V) = (tree.root = isdefined(tree, :root)
? add_node(tree.root, key, value): TreeNode(key, value); tree.root.isRed=false; tree)

push!{K,V}(tree::LLRBTree{K,V}, node::TreeNode{K,V}) = push!(tree, node.key.value, node.value.value)

function populate!{K,V}(tree::LLRBTree{K,V}, N::Int64)
    keys = []

    for i in 1:N
        r::Int64 = rand(1:N*2+100)
        push!(tree, r, r )

        #rep_index = indexin([r], keys)
        #while rep_index != [0]
        #    deleteat!(keys, rep_index)
        #    rep_index = indexin([r], keys)
        #end

        push!(keys, r)
    end
    return keys
end

function moveredleft{K,V}(node::TreeNode{K,V})
    flipcolor!(node)

    if !isrightleaf(node)
        if ( isred(node, :right, :left) )
            node.right = rotateright(node.right)
            node = rotateleft(node)
            flipcolor!(node)
        end
    end
    return node
end

function moveredright{K,V}(node::TreeNode{K,V})
    flipcolor!(node)

    if ( isred(node, :left, :left) )
        node= rotateright(node)
        flipcolor!(node)
    end
    return node
end

function minimum{K,V}(node::TreeNode{K,V})
    while !isleftleaf(node)
        node=node.left
    end
    return node
end
minimum(tree::LLRBTree) = minimum(tree.root)

function maximum{K,V}(node::TreeNode{K,V})
    while !isrightleaf(node)
        node=node.right
    end
    return node
end
maximum(tree::LLRBTree) = maximum(tree.root)

function fixup{K,V}(node::TreeNode{K,V})

    if isred(node, :right)
        node=rotateleft(node)
    end

    if(!isred(node, :left) && isred(node, :right) )
        node=rotateleft(node)
    end

    if(isred(node, :left) && isred(node, :left, :left) )
        node=rotateright(node)
    end

    if (isred(node, :left) && isred(node, :right) )
        flipcolor!(node)
    end

    return node
end

function deletemin{K,V}(node::TreeNode{K,V})

    if(!isleftleaf(node))

        if(!isred(node, :left) && !isred(node, :left, :left) )
            moveredleft(node)
        end

        node.left = deletemin(node.left)

        return fixup(node)
    else

        return TreeNode{K,V}()
    end
end

function delete_node{K,V}(node::TreeNode{K,V}, key::K)
    if key<node.key.value

        if !isleftleaf(node)

            if(!isred(node, :left) && !isred(node, :left, :left) )
                node=moveredleft(node)
            end

            node.left=delete_node(node.left, key)

        else
            Base.warn_once("Key not found:", key)
            return fixup(node)
        end

    else
        if(isred(node, :left))
            node=rotateright(node)
        end

        if (key==node.key.value && isrightleaf(node))
            return TreeNode{K,V}()
        end

        if !isrightleaf(node)

            if(!isred(node, :right) && !isred(node, :right, :left) )
                node=moveredright(node)
            end

            if key==node.key.value
                min = minimum(node.right)
                node.key=min.key
                node.value=min.value
                node.right=deletemin(node.right)
            else
                node.right=delete_node(node.right,key)
            end

        else
            Base.warn_once("Key not found:", key)
            return fixup(node)
        end
    end
    return fixup(node)
end

function delete!{K}(tree::LLRBTree, key::K)

    if( isdefined(tree, :root))
        if(! isleaf(tree.root) )
            tree.root = delete_node(tree.root, key)

            tree.root.isRed = false
        else
            Base.warn_once("Root cannot be leaf")
        end
    else
        Base.warn_once("Root cannot be leaf")
    end

    tree
end
delete!{K,V}(tree::LLRBTree, node::TreeNode{K,V}) = delete!(tree, node.key.value)

function levelorder{K,V}(tree::LLRBTree{K,V})
    queue = []
    push!(queue, tree.root)
    println("levelorder:  ")
    count = 20
    last20 = ones(1:20)

    while (!isempty(queue))
        node = shift!(queue)

        #Only print the first 20 values
        if(count>0)
            print(node.value.value,"  ")
        elseif count == 0
            print("...  ")
        end

        #Keep the last 20 entries to print later
        shift!(last20)
        #pop!(last20)
        push!(last20, node.value.value)
        #println(size(last20,1))

        if( !isleftleaf(node) )
            push!(queue, node.left)
        end
        if( !isrightleaf(node) )
            push!(queue, node.right)
        end

        count -= 1
    end

    if count < -19
        for i in 0:19
            print(last20[end-i]," ")
        end
    end
end

#Return an ordered list with recursion
function inorder{K,V}(node::TreeNode{K,V}, list=Any[])
  if !isleaf(node)

      if !isleftleaf(node)
          #("Go down to the left first")
          list = inorder(node.left, list)[2]
      end

      #("Push to the list when leftmost")
      push!(list, [node.key.value, node.value.value])

      if !isrightleaf(node)
          #("Go right if no left")
          list = inorder(node.right, list)[2]
      end
  end
  #("Go up;  ",node,";  ",list)
  return node, list
end
orderedpairs{K,V}(tree::LLRBTree{K,V}) = inorder(tree.root)[2]

end


text:
# Copyright (c) 2021 California Institute of Technology (“Caltech”).
# U.S. Government sponsorship acknowledged.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# - Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# - Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# - Neither the name of Caltech nor its operating division, the Jet Propulsion
#   Laboratory, nor the names of its contributors may be used to endorse or
#   promote products derived from this software without specific prior written
#   permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

######################### Radial Functions #########################3

function radial_greens_function(L, r1, r2)
    if r1 <= r2
        return r2 / (2*L+1) * (r1/r2)^L
    else
        return r2 / (2*L+1) * (r2/r1)^(L+1)
    end
end

function get_radial_greens_matrix(L, ri, Δr)
    G = fill(NaN, length(ri), length(ri))
    for j=1:length(ri), i=j:length(ri)
        G[i,j] = G[j,i] = (ri[i] / ri[j]) * radial_greens_function(L, ri[i], ri[j])
    end
    G .*= Δr
    return Hermitian(G)
end

function get_radial_cryofunks(L, ri, Δr)
    B = diagm(@. √Δr * ri)
    G = get_radial_greens_matrix(L, ri, Δr)
    e = eigen(G)
    k = @. 1 / √e.values
    p = sortperm(k)
    k = k[p]
    Z = e.vectors[:,p]
    g = inv(B) * Z
    ginv = Z' * B
    return k, g, ginv
end

# vim: set sw=4 et sts=4 :


text:
## log
# 2d
Base.log(R::Angle2d) = Angle2dGenerator(R.theta)
Base.log(R::RotMatrix{2}) = RotMatrixGenerator(log(Angle2d(R)))
#= We can define log for Rotation{2} like this,
but the subtypes of Rotation{2} are only Angle2d and RotMatrix{2},
so we don't need this defnition. =#
# Base.log(R::Rotation{2}) = log(Angle2d(R))

# 3d
Base.log(R::RotationVec) = RotationVecGenerator(R.sx,R.sy,R.sz)
Base.log(R::RotMatrix{3}) = RotMatrixGenerator(log(RotationVec(R)))
Base.log(R::Rotation{3}) = log(RotationVec(R))

# General dimensions
function Base.log(R::RotMatrix{N}) where N
    # This will be faster when log(::SMatrix) is implemented in StaticArrays.jl
    @static if VERSION < v"1.7"
        # This if block is related to this PR.
        # https://github.com/JuliaLang/julia/pull/40573
        S = SMatrix{N,N}(real(log(Matrix(R))))
    else
        S = SMatrix{N,N}(log(Matrix(R)))
    end
    RotMatrixGenerator((S-S')/2)
end

## exp
# 2d
Base.exp(R::Angle2dGenerator) = Angle2d(R.v)
Base.exp(R::RotMatrixGenerator{2}) = RotMatrix(exp(Angle2dGenerator(R)))
# Same as log(R::Rotation{2}), this definition is not necessary until someone add a new subtype of RotationGenerator{2}.
# Base.exp(R::RotationGenerator{2}) = exp(Angle2dGenerator(R))

# 3d
Base.exp(R::RotationVecGenerator) = RotationVec(R.x,R.y,R.z)
Base.exp(R::RotMatrixGenerator{3}) = RotMatrix(exp(RotationVecGenerator(R)))
# Same as log(R::Rotation{2}), this definition is not necessary until someone add a new subtype of RotationGenerator{3}.
# Base.exp(R::RotationGenerator{3}) = exp(RotationVecGenerator(R))

# General dimensions
Base.exp(R::RotMatrixGenerator{N}) where N = RotMatrix(exp(SMatrix(R)))


text:
# An array type for colorized overlays of grayscale images
"""
```
A = Overlay(channels, colors)
A = Overlay(channels, colors, clim)
A = Overlay(channels, colors, mapi)
```

Create an `Overlay` array from grayscale channels.  `channels = (channel1,
channel2, ...)`, `colors` is a vector or tuple of `Color`s, and `clim` is a
vector or tuple of min/max values, e.g., `clim = ((min1,max1),(min2,max2),...)`.
Alternatively, you can supply a list of `MapInfo` objects.

See also: `OverlayImage`.
"""
immutable Overlay{T,N,NC,AT<:Tuple{Vararg{AbstractArray}},MITypes<:Tuple{Vararg{MapInfo}}} <: AbstractArray{RGB{T},N}
    channels::AT   # this holds the grayscale arrays
    colors::NTuple{NC,RGB{T}}
    mapi::MITypes

    function Overlay(channels::Tuple{Vararg{AbstractArray}}, colors, mapi::Tuple{Vararg{MapInfo}})
        length(channels) == NC || error("Number of channels must match number of colors")
        length(mapi) == NC   || error("Number of ConvertInfo objects must match number of colors")
        for i = 2:NC
            size(channels[i]) == size(channels[1]) || error("All arrays must have the same size")
        end
        new(channels, colors, mapi)
    end
end
Overlay(channels::Tuple{Vararg{AbstractArray}}, colors::AbstractVector, mapi::Tuple{Vararg{MapInfo}}) =
    Overlay(channels,tuple(colors...),mapi)
Overlay{NC,T}(channels::Tuple{Vararg{AbstractArray}}, colors::NTuple{NC,RGB{T}}, mapi::Tuple{Vararg{MapInfo}}) =
    Overlay{T,ndims(channels[1]),NC,typeof(channels),typeof(mapi)}(channels,colors,mapi)

function Overlay(channels::Tuple{Vararg{AbstractArray}}, colors,
                 clim = ntuple(i->(zero(eltype(channels[i])), one(eltype(channels[i]))), length(channels)))
    n = length(channels)
    for i = 1:n
        if length(clim[i]) != 2
            error("clim must be a 2-vector")
        end
    end
    mapi = ntuple(i->ScaleMinMax(Float32, channels[i], clim[i][1], clim[i][2]), n)
    Overlay(channels, colors, mapi)
end

# Returns the overlay as an image, if possible
"`OverlayImage` is identical to `Overlay`, except that it returns an Image."
function OverlayImage(channels::Tuple{Vararg{AbstractArray}}, colors::Tuple{Vararg{Colorant}},
            arg = ntuple(i->(zero(eltype(channels[i])), one(eltype(channels[i]))), length(channels)))
    ovr = Overlay(channels, colors, arg)
    local prop
    haveprop = false
    for i = 1:length(channels)
        if isa(channels[i], AbstractImage)
            prop = copy(properties(channels[i]))
            haveprop = true
            break
        end
    end
    if !haveprop
        prop = properties(channels[1])
    end
    haskey(prop, "colorspace") && delete!(prop, "colorspace")
    haskey(prop, "colordim")   && delete!(prop, "colordim")
    haskey(prop, "limits")     && delete!(prop, "limits")
    Image(ovr, prop)
end

for NC = 1:3
    NCm = NC-1
    for K = 1:4
        indexargs = Symbol[symbol(string("i_",d)) for d = 1:K]
        sigargs = Expr[:($a::Integer) for a in indexargs]
        @eval begin
            function getindex{T,N,AT,MITypes}(O::Overlay{T,N,$NC,AT,MITypes}, $(sigargs...))
                # one of them needs a bounds-check
                out = map(O.mapi[$NC], getindex(O.channels[$NC], $(indexargs...))) * O.colors[$NC]
                @inbounds @nexprs $NCm c->begin
                    out += map(O.mapi[c], getindex(O.channels[c], $(indexargs...))) * O.colors[c]
                end
                clamp01(eltype(O), out)
            end
        end
    end
end

setindex!(O::Overlay, val, I::Real...) = error("Overlays are read-only. Convert to Image{RGB} to adjust values.")

#### Other Overlay support functions ####
eltype{T}(o::Overlay{T}) = RGB{T}
length(o::Overlay) = isempty(o.channels) ? 0 : length(o.channels[1])
size(o::Overlay) = isempty(o.channels) ? (0,) : size(o.channels[1])
size(o::Overlay, d::Integer) = isempty(o.channels) ? 0 : size(o.channels[1],d)
nchannels(o::Overlay) = length(o.channels)

similar(o::Overlay) = Array(eltype(o), size(o))
similar(o::Overlay, ::NTuple{0}) = Array(eltype(o), size(o))
similar{T}(o::Overlay, ::Type{T}) = Array(T, size(o))
similar{T}(o::Overlay, ::Type{T}, sz::Int64) = Array(T, sz)
similar{T}(o::Overlay, ::Type{T}, sz::Int64...) = Array(T, sz)
similar{T}(o::Overlay, ::Type{T}, sz) = Array(T, sz)

showcompact(io::IO, o::Overlay) = print(io, summary(o), " with colors ", o.colors)


text:
## A diagonal block is used to identify the shifts and to
## find the eigenvalues for the deflated matrix

## For the QFactorization (Descending) case provide a 2 x 2 view of the full matrix, A[j:j+1,j:j+1]
## For the Twisted case, provide a view of A[j:J+1, j:J+1] that is
## approximate unless the rotators at j-1 and J+1 are both identity rotators (or implicitly so)

##
## For QRFactorization, we exploit fact that QF is Hessenberg and
## R  is upper  triangular to simplify the matrix multiplication:
function diagonal_block!(A, QF::QFactorization, RF, j, J)

    i, k = j-1, j+1
    qji, qjj, qjk, qkj, qkk = QF[j,i], QF[j,j], QF[j,k], QF[k,j], QF[k,k]

    rjj, rjk = RF[j,j], RF[j,k]
    rkk = RF[k,k]
    if i >= 1
        rij, rik = RF[i,j], RF[i,k]
    else
        rij,  rik = 0*rkk, 0*rkk
    end

    ## This is matrix multiplication of a Hessenberg matrix times a upper triangular
    ## (triu(Q,-1) * triu(R))[k-1:k, k-1:k]
    ## could be just sum(QF[i,l] * RF[l,j] for l in i-1:j)

    A[1,1] = qji * rij + qjj * rjj
    A[1,2] = qji * rik + qjj * rjk + qjk * rkk
    A[2,1] = qkj * rjj
    A[2,2] = qkj * rjk + qkk * rkk

    return nothing
end

## Twisted case
function diagonal_block!(A, QF::QFactorizationTwisted, RF, j, k)
    approx_diagonal_block!(A, QF, RF, j, k)
end

## For twisted we use thee approximate diagonal block
## this is exact when the the j-1 and  k+1  rotators are deflated or j=1 or k=end
function approx_diagonal_block!(A, QF, RF, j, k) where {T, S}
    R = RF
    D = QF.D
    Q = QF.Q
    k = min(k, length(Q))
    # make 1:(k-j+1) identity
    A .= zero(eltype(A))
    if j == 1
        _approx_diagonal_block!(A, Q, D, R, j, k)
    else
        _approx_diagonal_block!(A, Q, D, R, j-1, k)
        # shift up and left by 1
        m,n = size(A)
        for j in 1:n-1  # across columns, then down rows
            for i in 1:m-1
                A[i,j] = A[i+1, j+1]
            end
        end
    end

    nothing
end

function _approx_diagonal_block!(A, Q, D, R, j, k)

    S = eltype(A)
    for i in j:k+1
        for ii in j:(i-1)
            A[i-j+1,ii-j+1] = zero(S)
        end
        for ii in i:k+1
            A[i-j+1,ii-j+1] = D[i] *  R[i,ii]
        end
    end

    # we get first one
    Tw = view(Q, j:k)
    for i in reverse(position_vector_indices(Tw.pv))  # allocates here
        U = Tw.x[i]
        V = Rotator(vals(U)..., idx(U) - j + 1)
        lmul!(V, A)
    end

    return  nothing
end

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

# [a11 - l a12; a21 a22] -> l^2 -2 * (tr(A)/2) l + det(A)
# so we use b = tr(A)/2 for qdrtc routing
## function eigen_values(state::AbstractQRFactorizationState)

##     # this allocates:
##     # e1, e2 = eigvals(state.A)
##     # return real(e1), imag(e1), real(e2), imag(e2)

##     a11, a12 = state.A[1,1], state.A[1,2]
##     a21, a22 = state.A[2,1], state.A[2,2]

##     eigen_values(a11, a12, a21, a22)
## end

# return tuple
function eigen_values(A::AbstractArray{T, N}) where {T, N}
    eigen_values(A[1,1], A[1,2], A[2,1], A[2,2])
end

## function eigen_values(A::AbstractArray{S, N}) where {S <: Complex, N}
##     eigen_values(A[1,1], A[1,2], A[2,1], A[2,2])
## end

function  eigen_values(a11::T, a12::T, a21::T, a22::T) where {T <: Real}
    b = (a11 + a22) / 2
    c = a11 * a22 - a12 * a21

    e1r, e1i, e2r, e2i = qdrtc(b,c) #qdrtc(one(b), b, c)
    return  complex(e1r, e1i), complex(e2r, e2i)

end

# from `modified_quadratic.f90`
function  eigen_values(a11::S, a12::S, a21::S, a22::S) where {S <: Complex}

    tr = a11 + a22
    detm = a11 * a22 - a21 * a12
    disc::S = sqrt(tr * tr - 4.0 * detm)

    u::S = abs(tr + disc) > abs(tr - disc) ? tr + disc : tr - disc
    if iszero(u)
        return zero(S), zero(S)
    else
        e1::S = u / 2.0
        e2::S = detm / e1
        return e1, e2
    end

end

#  this  makes some  algorithm type stable.
_eigvals(A::AbstractMatrix{T}) where {T <: Real} = eigvals(A)::Vector{Complex{T}}
_eigvals(A::AbstractMatrix{S}) where {S <: Complex} =  eigvals(A)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test

n= 10 #Size of matrix to test
srand(1)

@testset "Interconversion between special matrix types" begin
    a = [1.0:n;]
    A = Diagonal(a)
    @testset for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Matrix]
       @test Matrix(convert(newtype, A)) == Matrix(A)
       @test Matrix(convert(newtype, Diagonal(GenericArray(a)))) == Matrix(A)
    end

    @testset for isupper in (true, false)
        A = Bidiagonal(a, [1.0:n-1;], ifelse(isupper, :U, :L))
        for newtype in [Bidiagonal, Tridiagonal, Matrix]
           @test Matrix(convert(newtype, A)) == Matrix(A)
           @test Matrix(newtype(A)) == Matrix(A)
        end
        @test_throws ArgumentError convert(SymTridiagonal, A)
        tritype = isupper ? UpperTriangular : LowerTriangular
        @test Matrix(tritype(A)) == Matrix(A)

        A = Bidiagonal(a, zeros(n-1), ifelse(isupper, :U, :L)) #morally Diagonal
        for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Matrix]
           @test Matrix(convert(newtype, A)) == Matrix(A)
           @test Matrix(newtype(A)) == Matrix(A)
        end
        @test Matrix(tritype(A)) == Matrix(A)
    end

    A = SymTridiagonal(a, [1.0:n-1;])
    for newtype in [Tridiagonal, Matrix]
       @test Matrix(convert(newtype, A)) == Matrix(A)
    end
    for newtype in [Diagonal, Bidiagonal]
       @test_throws ArgumentError convert(newtype,A)
    end
    A = SymTridiagonal(a, zeros(n-1))
    @test Matrix(convert(Bidiagonal,A)) == Matrix(A)

    A = Tridiagonal(zeros(n-1), [1.0:n;], zeros(n-1)) #morally Diagonal
    for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Matrix]
       @test Matrix(convert(newtype, A)) == Matrix(A)
    end
    A = Tridiagonal(ones(n-1), [1.0:n;], ones(n-1)) #not morally Diagonal
    for newtype in [SymTridiagonal, Matrix]
       @test Matrix(convert(newtype, A)) == Matrix(A)
    end
    for newtype in [Diagonal, Bidiagonal]
        @test_throws ArgumentError convert(newtype,A)
    end
    A = Tridiagonal(zeros(n-1), [1.0:n;], ones(n-1)) #not morally Diagonal
    @test Matrix(convert(Bidiagonal, A)) == Matrix(A)
    A = UpperTriangular(Tridiagonal(zeros(n-1), [1.0:n;], ones(n-1)))
    @test Matrix(convert(Bidiagonal, A)) == Matrix(A)
    A = Tridiagonal(ones(n-1), [1.0:n;], zeros(n-1)) #not morally Diagonal
    @test Matrix(convert(Bidiagonal, A)) == Matrix(A)
    A = LowerTriangular(Tridiagonal(ones(n-1), [1.0:n;], zeros(n-1)))
    @test Matrix(convert(Bidiagonal, A)) == Matrix(A)
    @test_throws ArgumentError convert(SymTridiagonal,A)

    A = LowerTriangular(Matrix(Diagonal(a))) #morally Diagonal
    for newtype in [Diagonal, Bidiagonal, SymTridiagonal, LowerTriangular, Matrix]
        @test Matrix(convert(newtype, A)) == Matrix(A)
    end
    A = UpperTriangular(Matrix(Diagonal(a))) #morally Diagonal
    for newtype in [Diagonal, Bidiagonal, SymTridiagonal, UpperTriangular, Matrix]
        @test Matrix(convert(newtype, A)) == Matrix(A)
    end
    A = UpperTriangular(triu(rand(n,n)))
    for newtype in [Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal]
        @test_throws ArgumentError convert(newtype,A)
    end
end

@testset "Binary ops among special types" begin
    a=[1.0:n;]
    A=Diagonal(a)
    Spectypes = [Diagonal, Bidiagonal, Tridiagonal, Matrix]
    for (idx, type1) in enumerate(Spectypes)
        for type2 in Spectypes
           B = convert(type1,A)
           C = convert(type2,A)
           @test Matrix(B + C) ≈ Matrix(A + A)
           @test Matrix(B - C) ≈ Matrix(A - A)
       end
    end
    B = SymTridiagonal(a, ones(n-1))
    for Spectype in [Diagonal, Bidiagonal, Tridiagonal, Matrix]
        @test Matrix(B + convert(Spectype,A)) ≈ Matrix(B + A)
        @test Matrix(convert(Spectype,A) + B) ≈ Matrix(B + A)
        @test Matrix(B - convert(Spectype,A)) ≈ Matrix(B - A)
        @test Matrix(convert(Spectype,A) - B) ≈ Matrix(A - B)
    end

    C = rand(n,n)
    for TriType in [Base.LinAlg.UnitLowerTriangular, Base.LinAlg.UnitUpperTriangular, UpperTriangular, LowerTriangular]
        D = TriType(C)
        for Spectype in [Diagonal, Bidiagonal, Tridiagonal, Matrix]
            @test Matrix(D + convert(Spectype,A)) ≈ Matrix(D + A)
            @test Matrix(convert(Spectype,A) + D) ≈ Matrix(A + D)
            @test Matrix(D - convert(Spectype,A)) ≈ Matrix(D - A)
            @test Matrix(convert(Spectype,A) - D) ≈ Matrix(A - D)
        end
    end
end

@testset "Triangular Types and QR" begin
    for typ in [UpperTriangular,LowerTriangular,Base.LinAlg.UnitUpperTriangular,Base.LinAlg.UnitLowerTriangular]
        a = rand(n,n)
        atri = typ(a)
        b = rand(n,n)
        qrb = qrfact(b,Val(true))
        @test Base.LinAlg.A_mul_Bc(atri,qrb[:Q]) ≈ Matrix(atri) * qrb[:Q]'
        @test Base.LinAlg.A_mul_Bc!(copy(atri),qrb[:Q]) ≈ Matrix(atri) * qrb[:Q]'
        qrb = qrfact(b,Val(false))
        @test Base.LinAlg.A_mul_Bc(atri,qrb[:Q]) ≈ Matrix(atri) * qrb[:Q]'
        @test Base.LinAlg.A_mul_Bc!(copy(atri),qrb[:Q]) ≈ Matrix(atri) * qrb[:Q]'
    end
end

# should all yield sparse arrays
@testset "concatenations of combinations of special and other matrix types" begin
    N = 4
    # Test concatenating pairwise combinations of special matrices
    diagmat = Diagonal(ones(N))
    bidiagmat = Bidiagonal(ones(N), ones(N-1), :U)
    tridiagmat = Tridiagonal(ones(N-1), ones(N), ones(N-1))
    symtridiagmat = SymTridiagonal(ones(N), ones(N-1))
    specialmats = (diagmat, bidiagmat, tridiagmat, symtridiagmat)
    for specialmata in specialmats, specialmatb in specialmats
        @test issparse(hcat(specialmata, specialmatb))
        @test issparse(vcat(specialmata, specialmatb))
        @test issparse(hvcat((1,1), specialmata, specialmatb))
        @test issparse(cat((1,2), specialmata, specialmatb))
    end
    # Test concatenating pairwise combinations of special matrices with sparse matrices,
    # dense matrices, or dense vectors
    densevec = ones(N)
    densemat = diagm(ones(N))
    spmat = spdiagm(ones(N))
    for specialmat in specialmats
        # --> Tests applicable only to pairs of matrices
        for othermat in (spmat, densemat)
            @test issparse(vcat(specialmat, othermat))
            @test issparse(vcat(othermat, specialmat))
        end
        # --> Tests applicable also to pairs including vectors
        for specialmat in specialmats, othermatorvec in (spmat, densemat, densevec)
            @test issparse(hcat(specialmat, othermatorvec))
            @test issparse(hcat(othermatorvec, specialmat))
            @test issparse(hvcat((2,), specialmat, othermatorvec))
            @test issparse(hvcat((2,), othermatorvec, specialmat))
            @test issparse(cat((1,2), specialmat, othermatorvec))
            @test issparse(cat((1,2), othermatorvec, specialmat))
        end
    end
end

# Test that concatenations of annotated sparse/special matrix types with other matrix
# types yield sparse arrays, and that the code which effects that does not make concatenations
# strictly involving un/annotated dense matrices yield sparse arrays
#
# TODO: As with the associated code, these tests should be moved to a more appropriate
# location, particularly some future equivalent of base/linalg/special.jl dedicated to
# intereactions between a broader set of matrix types
@testset "concatenations of annotated types" begin
    N = 4
    # The tested annotation types
    testfull = Bool(parse(Int,(get(ENV, "JULIA_TESTFULL", "0"))))
    utriannotations = (UpperTriangular, Base.LinAlg.UnitUpperTriangular)
    ltriannotations = (LowerTriangular, Base.LinAlg.UnitLowerTriangular)
    triannotations = (utriannotations..., ltriannotations...)
    symannotations = (Symmetric, Hermitian)
    annotations = testfull ? (triannotations..., symannotations...) : (LowerTriangular, Symmetric)
    # Concatenations involving these types, un/annotated, should yield sparse arrays
    spvec = spzeros(N)
    spmat = speye(N)
    diagmat = Diagonal(ones(N))
    bidiagmat = Bidiagonal(ones(N), ones(N-1), :U)
    tridiagmat = Tridiagonal(ones(N-1), ones(N), ones(N-1))
    symtridiagmat = SymTridiagonal(ones(N), ones(N-1))
    sparseconcatmats = testfull ? (spmat, diagmat, bidiagmat, tridiagmat, symtridiagmat) : (spmat, diagmat)
    # Concatenations involving strictly these types, un/annotated, should yield dense arrays
    densevec = ones(N)
    densemat = ones(N, N)
    # Annotated collections
    annodmats = [annot(densemat) for annot in annotations]
    annospcmats = [annot(spcmat) for annot in annotations, spcmat in sparseconcatmats]
    # Test that concatenations of pairwise combinations of annotated sparse/special
    # yield sparse matrices
    for annospcmata in annospcmats, annospcmatb in annospcmats
        @test issparse(vcat(annospcmata, annospcmatb))
        @test issparse(hcat(annospcmata, annospcmatb))
        @test issparse(hvcat((2,), annospcmata, annospcmatb))
        @test issparse(cat((1,2), annospcmata, annospcmatb))
    end
    # Test that concatenations of pairwise combinations of annotated sparse/special
    # matrices and other matrix/vector types yield sparse matrices
    for annospcmat in annospcmats
        # --> Tests applicable to pairs including only matrices
        for othermat in (densemat, annodmats..., sparseconcatmats...)
            @test issparse(vcat(annospcmat, othermat))
            @test issparse(vcat(othermat, annospcmat))
        end
        # --> Tests applicable to pairs including other vectors or matrices
        for other in (spvec, densevec, densemat, annodmats..., sparseconcatmats...)
            @test issparse(hcat(annospcmat, other))
            @test issparse(hcat(other, annospcmat))
            @test issparse(hvcat((2,), annospcmat, other))
            @test issparse(hvcat((2,), other, annospcmat))
            @test issparse(cat((1,2), annospcmat, other))
            @test issparse(cat((1,2), other, annospcmat))
        end
    end
    # The preceding tests should cover multi-way combinations of those types, but for good
    # measure test a few multi-way combinations involving those types
    @test issparse(vcat(spmat, densemat, annospcmats[1], annodmats[2]))
    @test issparse(vcat(densemat, spmat, annodmats[1], annospcmats[2]))
    @test issparse(hcat(spvec, annodmats[1], annospcmats[3], densevec, diagmat))
    @test issparse(hcat(annodmats[2], annospcmats[4], spvec, densevec, diagmat))
    @test issparse(hvcat((5,), diagmat, densevec, spvec, annodmats[1], annospcmats[1]))
    @test issparse(hvcat((5,), spvec, annodmats[2], diagmat, densevec, annospcmats[2]))
    @test issparse(cat((1,2), annodmats[1], diagmat, annospcmats[3], densevec, spvec))
    @test issparse(cat((1,2), spvec, diagmat, densevec, annospcmats[4], annodmats[2]))
    # Test that concatenations strictly involving un/annotated dense matrices/vectors
    # yield dense arrays
    for densemata in (densemat, annodmats...)
        # --> Tests applicable to pairs including only matrices
        for densematb in (densemat, annodmats...)
            @test !issparse(vcat(densemata, densematb))
            @test !issparse(vcat(densematb, densemata))
        end
        # --> Tests applicable to pairs including vectors or matrices
        for otherdense in (densevec, densemat, annodmats...)
            @test !issparse(hcat(densemata, otherdense))
            @test !issparse(hcat(otherdense, densemata))
            @test !issparse(hvcat((2,), densemata, otherdense))
            @test !issparse(hvcat((2,), otherdense, densemata))
            @test !issparse(cat((1,2), densemata, otherdense))
            @test !issparse(cat((1,2), otherdense, densemata))
        end
    end
end
@testset "vcat of Vectors with SparseVectors should yield SparseVector (#22225)" begin
    @test isa((@inferred vcat(Float64[], spzeros(1))), SparseVector)
end

max_stars_repo_name: anowacki/julia

max_issues_repo_name: anowacki/julia

max_forks_repo_name: anowacki/julia


text:
"""
A shape where anything inside could cause a disturbance at the listener position
from a point source wavefront starting at the listener. Also everything inside
has a positive `x` coordinate. It is equivalent to a segment of a circle.

More precisely, if the listener is at (l_x,l_y) then the interior of the shape
is defined as
sqrt((x-l_x)^2+(y-l_y)^2)<time and x>0
"""
struct TimeOfFlightFromPoint{T <: AbstractFloat} <: Shape{T,2}
    listener_position::Vector{T}
    time::T
end

name(shape::TimeOfFlightFromPoint) = "Time of flight from point source"

function volume(shape::TimeOfFlightFromPoint)
    θ = 2acos(-shape.listener_position[1] / shape.time)
    return shape.time^2 * (θ - sin(θ))/2
end

import Base.issubset
function issubset(circle::Sphere{T,2}, shape::TimeOfFlightFromPoint) where T
    (origin(circle)[1] - circle.radius) > 0 &&
    norm(origin(circle) - shape.listener_position) < (shape.time - 2circle.radius)
end

function bounding_box(shape::TimeOfFlightFromPoint{T}) where T <: AbstractFloat
    box_height = 2sqrt(shape.time^2 - shape.listener_position[1]^2)
    box_width = max(shape.time + shape.listener_position[1], zero(T))
    return Box([SVector(zero(T), -box_height/2), SVector(box_width, box_height/2)])
end

function boundary_functions(shape::TimeOfFlightFromPoint)

    function x(t)
        check_boundary_coord_range(t)
        if t <= 1//2
            θ = acos(-shape.listener_position[1] / shape.time)
            return shape.time*cos(θ*(4t-1)) + shape.listener_position[1]
        else
            return zero(t)
        end
    end

    function y(t)
        check_boundary_coord_range(t)
        if t <= 1//2
            θ = acos(-shape.listener_position[1] / shape.time)
            return shape.time*sin(θ*(4t-1)) + shape.listener_position[2]
        else
            return 2*(3//4-t)*2*sqrt(shape.time^2 - shape.listener_position[1]^2)
        end
    end

    return x, y
end