"""
    SimpleTimeSteppers

Implements a few simple time integration schemes useful for testing the accuracy of spatial discretisations.
"""
module SimpleTimeSteppers

export rk4
export rkn4
export centered_2nd_order_fd

"""
    rk4(F, vₙ, tₙ, Δtₙ)

Solves ``vₜ = F(v,t)`` using the standard 4th order RK-method.
"""
function rk4(F, vₙ, tₙ, Δtₙ)
    k₁ = F(vₙ,tₙ)
    k₂ = F(vₙ + Δtₙ/2*k₁, tₙ + Δtₙ/2)
    k₃ = F(vₙ + Δtₙ/2*k₂, tₙ + Δtₙ/2)
    k₄ = F(vₙ + Δtₙ  *k₃, tₙ + Δtₙ  )

    vₙ₊₁ = vₙ + (1/6)*(k₁+2*(k₂+k₃)+k₄)*Δtₙ

    return vₙ₊₁
end

"""
    rkn4(F, vₙ, v̇ₙ, tₙ, Δtₙ)

Solves ``vₜₜ = F(v,v̇,t)`` using the Runge-Kutta-Nyström method based on the
standard 4th order RK-method.
"""
function rkn4(F, vₙ, v̇ₙ, tₙ, Δtₙ)
    k₁ = F(vₙ,                          v̇ₙ,            tₙ        )
    k₂ = F(vₙ + Δtₙ/2*v̇ₙ,               v̇ₙ + Δtₙ/2*k₁, tₙ + Δtₙ/2)
    k₃ = F(vₙ + Δtₙ/2*v̇ₙ + Δtₙ^2/4*k₁,  v̇ₙ + Δtₙ/2*k₂, tₙ + Δtₙ/2)
    k₄ = F(vₙ + Δtₙ  *v̇ₙ + Δtₙ^2/2*k₂,  v̇ₙ + Δtₙ  *k₃, tₙ + Δtₙ  )

    vₙ₊₁ = vₙ + Δtₙ*v̇ₙ + Δtₙ^2*(1/6)*(k₁ + k₂ + k₃);
    v̇ₙ₊₁ = v̇ₙ + Δtₙ*(k₁ + 2*k₂ + 2*k₃ + k₄)/6;

    return vₙ₊₁, v̇ₙ₊₁
end

"""
    centered_2nd_order_fd(F, vₙ, vₙ₋₁, Δtₙ)

Solves ``vₜₜ = F(v,t)`` using the 2nd order discretization
``(vₙ₊₁ - 2vₙ + vₙ₋₁)/Δt² = F(vₙ,tₙ)``.
"""
function centered_2nd_order_fd(F, vₙ, vₙ₋₁, tₙ, Δt)
    vₙ₊₁  = Δt^2*F(vₙ,tₙ) + 2vₙ - vₙ₋₁

    return vₙ₊₁, vₙ
end

# TODO: Rewrite without allocations


# TODO: Add euler forward and try a local error estimate for testing
end # module