using LinearAlgebra # for norm, dot, and cross
using PyPlot        # for plotting
using SymPy         # for symbolic math
using Roots         # to find zeros numerically
using QuadGK        # to integerate numerically
using ForwardDiff   # for derivatives

uvec(r) = r / norm(r) # unit vector

D(r::Function, n=1) = n > 1 ? D(D(r),n-1) : t -> ForwardDiff.derivative(r, float(t))
Base.adjoint(r::Function) = D(r)

xs_ys(vs) = (A=hcat(vs...); Tuple([A[i,:] for i in eachindex(vs[1])]))
xs_ys(v,vs...) = xs_ys([v, vs...])
xs_ys(r::Function, a, b, n=100) = xs_ys(r.(range(a, stop=b, length=n)))

function arrow!(plt::Plots.Plot, p, v; kwargs...)
  if length(p) == 2
     quiver!(plt, xs_ys([p])..., quiver=Tuple(xs_ys([v])); kwargs...)
  elseif length(p) == 3
    # 3d quiver needs support
    # https://github.com/JuliaPlots/Plots.jl/issues/319#issue-159652535
    # headless arrow instead
    plot!(plt, xs_ys(p, p+v)...; kwargs...)
	end
end
arrow!(p,v;kwargs...) = arrow!(Plots.current(), p, v; kwargs...)



























using SymPy
@vars t real=true