I am having problems getting event detection to work properly in DifferentialEquations.jl. I have reduced the code as far as possible:
using DifferentialEquations using PyPlot function evaluate!(t, y, δy) rm = norm(y[1:3]) δy[1:3] = y[4:6] δy[4:6] = -mu .* y[1:3] ./ rm^3 end function condition(t, y, integrator) r = y[1:3] v = y[4:6] rm = norm(r) vm = norm(v) ecc = norm(((vm^2 - mu/rm)*r - v * (r ⋅ v))/mu) sma = -mu/(2*(vm^2/2 - mu/rm)) ese = (r ⋅ v) / sqrt(mu * sma) ece = rm * vm^2 / mu - 1.0 E = atan2(ese, ece) ano = 2*atan2(sqrt(1 + ecc)*sin(E/2), sqrt(1 - ecc)*cos(E/2)) push!(ts, t) push!(as, ano) ano end function affect!(integrator) push!(te, integrator.t) end const mu = 398600.4418 ts = Vector() as = Vector() te = Vector() t1 = 5562.095583019713 r = [ 6068279.27, -1692843.94, -2516619.18, ]/1000 v = [ -660.415582, 5495.938726, -5303.093233, ]/1000 y0 = [r; v] prob = ODEProblem( evaluate!, y0, (0.0, t1), callback=ContinuousCallback(condition, affect!), ) res = solve(prob, Vern9()) idx = sortperm(ts) t = ts[idx] ano = as[idx] plot(t, ano) for t in te plot([te, te], [-pi, pi], "r") end grid()
This is the output:
The blue line is the output of the
condition function. I would expect the callback to only be triggered at the upcrossing between 4000 and 5000 seconds. But somehow it also gets triggered at the downcrossing. Am I doing something wrong?