final

.gitignore

@@ -1,4 +1,5 @@
 *Manifest.toml
+!examples/Manifest.toml
 
 docs/build/
 docs/site/

Project.toml

@@ -4,8 +4,11 @@ authors = ["Michael Kraus <[email protected]>"]
 version = "0.4.0"
 
 [deps]
+ChebyshevTransforms = "2752694e-a3a0-42bd-9f17-5f3242da1ce6"
+CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]

docs/make.jl

@@ -8,14 +8,6 @@ DocMeta.setdocmeta!(
     end
 )
 
-DocMeta.setdocmeta!(
-    PoincareInvariants.SecondChebyshevPlans.PaduaTransforms,
-    :DocTestSetup,
-    quote
-        using PoincareInvariants.SecondChebyshevPlans.PaduaTransforms
-    end
-)
-
 DocMeta.setdocmeta!(
     PoincareInvariants.CanonicalSymplecticForms,
     :DocTestSetup,
@@ -29,18 +21,18 @@ makedocs(
     modules=[PoincareInvariants],
     pages = [
         "Home" => "index.md",
-        "First Poincaré Invariants" => "first_poincare_invariants.md",
-        "Second Poincaré Invariants" => [
-            "Guide" => "second_poincare_invariants.md",
-            "Chebyshev Implementation" => "chebyshev_implementation.md",
-            "Padua Transforms" => "padua_transforms.md"
+        "Theory" => "theory.md",
+        "Guides" => [
+            "Pendulum" => "guides/pendulum.md",
+            "Charged Particle" => "guides/charged_particle.md",
+            "DifferentialEquations.jl" => "guides/diffeq.md",
+            "Plans" => "guides/plans.md",
         ],
-        "Canonical Symplectic Structures" => "canonical_symplectic_structures.md",
         "Reference" => "reference.md"
     ],
     format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true")
 )
 
-deploydocs(
-    repo   = "github.com/JuliaGNI/PoincareInvariants.jl"
-)
+# deploydocs(
+#     repo   = "github.com/JuliaGNI/PoincareInvariants.jl"
+# )

docs/src/guides/charged_particle.md

@@ -0,0 +1,185 @@
+# Charged Particle
+
+In the second part of the tutorial we will treat a non-canonical system and show how to [`compute!`](@ref) Poincaré integral invariants for it.
+The system consists of a 2D charged particle subject to electromagnetic fields.
+The equations of motion are
+
+```math
+\dot{x} = v,\; \dot{v} = E(x) + v × B(x)
+```
+
+where $x$ is the position vector, $v$ is the velocity, $E$ and $B$ are the electric and magnetic field and the mass and charge of the particle is set to one. For our example we will use a constant magnetic field of strength $10$ pointing in the z-direction and an electric field given by $E(x, y) = (-x, -y^3)$.
+
+## Integration
+
+Again, we'll quickly implement our own integrators here, namely the [Runge-Kutta method (RK4)](https://en.wikipedia.org/wiki/Runge–Kutta_methods#The_Runge–Kutta_method) and a [second order symplectic splitting method for charged particles in electromagnetic fields](https://doi.org/10.1016/j.physleta.2016.12.031).
+
+```julia
+using PoincareInvariants
+using StaticArrays
+
+E(x, y) = (-x, -y^3)
+
+A(x, y) = 5 .* (-y, x)
+B(x, y) = 10.0
+Bx(x, y, Δx) = 10 * Δx  # integral in x direction from x to x + Δx
+By(x, y, Δy) = 10 * Δy  # integral in y direction from y to y + Δy
+
+struct Split2 end
+
+function ϕx((x, y, vx, vy), dt)
+    Δx = vx * dt
+    return (x + Δx, y, vx, vy - Bx(x, y, Δx))
+end
+
+function ϕy((x, y, vx, vy), dt)
+    Δy = vy * dt
+    return (x, y + Δy, vx + By(x, y, Δy), vy)
+end
+
+function ϕE((x, y, vx, vy), dt)
+    (Δvx, Δvy) = dt .* E(x, y)
+    return (x, y, vx + Δvx, vy + Δvy)
+end
+
+function timestep(z, dt, ::Split2)
+    hdt = 0.5 * dt
+    z = ϕx(z, hdt)
+    z = ϕy(z, hdt)
+    z = ϕE(z, dt)
+    z = ϕy(z, hdt)
+    z = ϕx(z, hdt)
+    return z
+end
+
+struct RK4 end
+
+function zdot((x, y, vx, vy))
+    ex, ey = E(x, y); b = B(x, y)
+    (vx, vy, ex + b * vy, ey - b * vx)
+end
+
+function timestep(z, dt, ::RK4)
+    hdt = 0.5 .* dt
+    k1 = zdot(z)
+    k2 = zdot(z .+ hdt .* k1)
+    k3 = zdot(z .+ hdt .* k2)
+    k4 = zdot(z .+  dt .* k3)
+    return z .+ dt .* (k1 .+ 2 .* k2 .+ 2 .* k3 .+ k4) .* (1/6)
+end
+
+"""
+    integrate(z0, dt, nsteps, nt, method)
+
+start at `z0` and integrate the equations of motion using `method`.
+Returns the timeseries as a vector of tuples. `nt` points are saved,
+`nsteps` steps are taken from saved point to saved point and `dt` is the
+size of each time step.
+"""
+function integrate(z0, dt, nsteps, nt, method)
+    out = Vector{NTuple{4, Float64}}(undef, nt)
+    z = out[1] = z0
+    for i in 2:nt
+        for _ in 1:nsteps
+            z = timestep(z, dt, method)
+        end
+        out[i] = z
+    end
+    return out
+end
+
+integrate(mat::AbstractMatrix, dt, nsteps, nt, method) = map(eachrow(mat)) do r
+    integrate((r[1], r[2], r[3], r[4]), dt, nsteps, nt, method)
+end
+```
+
+## Computing the Invariants
+
+Having gotten integration out of the way, we can move onto the invariants. For this system, the first invariant is given by
+
+```math
+I_{1} = \int_{\gamma} v(x) + A(x) \, dx
+```
+
+where $A(x)$ is the magnetic vector potential as a function of the position vector $x$. The second invariant is given by
+
+```math
+I_{2} = \int_{S} \omega_{ij} (q) \, dz^{i} \, dz^{j}
+```
+
+with the two-form $\omega$ given by
+
+```math
+\begin{pmatrix}
+ 0 & B & -1 &  0 \\
+-B & 0 &  0 & -1 \\
+ 1 & 0 &  0 &  0 \\
+ 0 & 1 &  0 &  0
+\end{pmatrix}
+```
+
+In code, we can define these forms as follows.
+
+```julia
+function oneform((x, y, vx, vy), ::Real, ::Any)
+    p = (vx, vy) .+ A(x, y)
+    @SVector [p[1], p[2], 0, 0]
+end
+
+function twoform(z, ::Real, ::Any)
+    b = B(z[1], z[2])
+    @SMatrix [ 0  b -1  0;
+              -b  0  0 -1;
+               1  0  0  0;
+               0  1  0  0]
+end
+```
+
+Forms always have the signature `form(z, t, p)`, where `z` is the phasespace position, `t` is a time and `p` is an arbitrary parameter.
+To use these forms to compute our invariants, we must create the setup objects `FirstPI` and `SecondPI` and initialise the curve and surface we want in phasespace.
+
+```julia
+pi1 = FirstPI{Float64, 4}(oneform, 1_000)
+pi2 = SecondPI{Float64, 4}(twoform, 10_000)
+
+I1 = 0
+pnts1 = getpoints(pi1) do θ
+    10 .* (0, 0, sinpi(2θ), cospi(2θ))
+end
+
+I2 = 1_000
+pnts2 = getpoints(pi2) do x, y
+    10 .* (y - 0.5, x - 0.5, 0, 0)
+end
+```
+
+The curve is just a point in space, so the first integral invariant evaluates to zero.
+The second invariant is equal to $1000$ since, it is equal to the magnetic field of strength $10$ integrated over a $10$ by $10$ square.
+We can confirm this by computing the inital invariants.
+
+```julia
+@assert isapprox(compute!(pi1, pnts1, 52.3, "optional parameter"), I1; atol=10^(-15))
+@assert isapprox(compute!(pi2, pnts2), I2; atol=10^(-11))
+```
+
+By default `p` is `nothing` and `t` is zero, but as shown in the first line we can also explicitly supply a time and parameter if the forms we defined earlier had required them. [`compute!`](@ref) used on a time series accepts an iterable of times and an arbitrary parameter.
+
+```julia
+times = range(0.0; step=0.05 * 50, length=5)
+
+series1 = integrate(pnts1, 0.05, 50, 5, Split2())
+@assert all(compute!(pi1, series1, times, ("optional parameters", 3.7, 42))) do I
+    abs(I - I1) < 10^(-13)
+end
+
+series2 = integrate(pnts2, 0.05, 50, 5, Split2())
+@assert all(compute!(pi2, series2)) do I
+    abs(I - I2) < 10^(-11)
+end
+```
+
+Again, we can plot the results. The top two plots show the increasingly distorted curve and surface over time (projected onto x and y position components), while the bottom two plots show the error in the invariant over time for the two integration algorithms.
+
+![](particle.png)
+
+We see that RK4 does not preserve the non-canonical invariant while the splitting method does.

docs/src/guides/diffeq.md

@@ -0,0 +1,56 @@
+# Compute Invariants With DifferentialEquations.jl
+
+Before reading this, it's worth familiarising yourself with the [Ensemble Simulations](https://diffeq.sciml.ai/stable/features/ensemble/#ensemble) section of the DifferentialEquations.jl docs.
+
+To compute invariants with DifferentialEquations.jl, we'll need to create a `PIEnsembleProblem`. This works as follows
+
+```Julia
+PIEnsembleProblem(init, prob, pinv;
+    prob_func = (prob,i,repeat)->prob,
+    output_func = (sol,i) -> (sol,false),
+    reduction = (u,data,I)->(append!(u,data),false),
+    u_init = [], safetycopy = false
+)
+```
+
+`init` is the phasespace curve or surface parameterisation, which determines the initial conditions;
+`prob` is a problem from DifferentialEquations.jl, which specifies the time span and differential equation;
+`pinv` is an AbstractPoincareInvariant;
+`prob_func` allows the user to `remake` the problem for each trajectory in the ensemble;
+all further arguments are exactly like the `EnsembleProblem` type from DifferentialEquations.jl. You need not add a `prob_func` of your own to have the correct initial conditions. These will be entered automaticaly according ot the given setup objet and parameterisation. However, by default, the type of the initial condition is `Vector{T}`. If you want to use an `ArrayPartion` or a `StaticArray` or other type to represent a point in phase space, you should specify a `prob_func` to do the conversion.
+
+All in all, this can look like so:
+
+```Julia
+using OrdinaryDiffEq
+using RecursiveArrayTools: ArrayPartition
+using PoincareInvariants
+
+dt = 0.1
+prob = SecondOrderODEProblem((p, θ, params, t) -> [-sin(θ[1])], 0.0, 0.0, (0.0, 2.0))
+pf(prob, i, repeat) = remake(prob; u0 = ArrayPartition((prob.u0[1:1], prob.u0[2:2])))
+
+pi1 = CanonicalFirstPI{Float64, 2}(1_000)
+ens_prob = PIEnsembleProblem(ϕ -> (sinpi(2ϕ), 3 * cospi(2ϕ)), prob, pi1; prob_func=pf)
+```
+
+Once, we have our `PIEnsembleProblem`, we can call solve on it, just like the `EnsembleProblem` it wraps. The general pattern is
+
+```Julia
+solve(prob::PIEnsembleProblem, alg, ensemblealg; kwargs...)
+```
+
+Most arguments are simply passed to the solve function called on the wrapped `EnsembleProblem`.
+The `trajectories` keyword argument is already set as `getpointnum(pinv)` and the `adaptive` keyword is set to false by default. Otherwise the behaviour is identical. For our example, we would have
+
+```Julia
+sol = solve(ens_prob, SymplecticEuler(), EnsembleSerial(); dt=dt)
+```
+
+Finally, we can call `compute!` on the result, which returns a `Vector{T}` of invariant values for every saved timestep. You must make sure all trajectories save at the same time steps.
+
+```Julia
+compute!(pi1, sol[, p])
+```
+
+The times and points are taken from the solution. An optional parameter `p` may be given, which is passed to the differential form.