should I define a differential equation’s callback in the main system such that it u[1] is bracketed

in [-pi, pi] range?

so yes I tried call back function defined GitHub - awage/Basins.jl on this example

this is code for replication

```
function rotor_1(du, u, params, t)
θ1,pθ1,θ2,pθ2 = u
k,l = params
r = ((2 + l^2 - 2*l*(cos(θ1)- cos(θ2)) - 2*cos(θ2-θ1))^1/2)
h = exp(-k*r) * (1+k*r) / r^3;
du[1] = pθ1;
du[2] = (h*(sin(θ1-θ2)+l*sin(θ1)))
du[3] = pθ2;
du[4] = (h*(sin(θ2-θ1)-l*sin(θ2)))
end
u0 = [-pi,0.,pi,0.]
params = [2.0,1.5]
# now use callback
# We have to define a callback to wrap the phase in [-π,π]
function affect!(integrator)
if integrator.u[1] < 0
integrator.u[1] += 2*π
else
integrator.u[1] -= 2*π
end
end
condition(u,t,integrator) = (integrator.u[1] <= -π || integrator.u[1] >= π)
cb = DiscreteCallback(condition,affect!)
diffeq = (alg = Vern7(), reltol = 1.0e-10, abstol = 1.0e-10, atol = 1.0e-10, rtol = 1.0e-10, callback=cb)
# coupled_ode
rotor_system = CoupledODEs(rotor_1, u0, params; diffeq,t0=0.0)
```

now let’s check the trajectory generated now

```
Points, t = trajectory(rotor_system, 50000.0; Ttr = 0.0, Δt =0.01)
maximum(Points[:,1])
# 3.4152245610206577
minimum(Points[:,1])
# -3.4250833268013063
maximum(Points[:,3])
# 3.1415926536003895
minimum(Points[:,3])
# -3.1270534882495484
```

so callback works but for θ1 I’m missing the bracket of limits (help needed) and now

energy checking

```
function hamil_cal(p)
θ1,pθ1,θ2,pθ2 = p
k, l = params
r = (sqrt(abs(2 + l^2 - 2*l*(cos(θ1)- cos(θ2)) - 2*cos(θ2-θ1))))
H = ((pθ1^2+pθ2^2)/2) + (exp(-k*r)/r)
return H
end
energy = Vector{Float64}(undef,0)
energy_abs = Vector{Float64}(undef,0)
function energy_check(Points)
global energy = Vector{Float64}(undef,0)
global energy_abs = Vector{Float64}(undef,0)
@showprogress for u in Points
push!(energy,hamil_cal(u))
push!(energy_abs,abs(hamil_cal(u)-hamil_cal(u0)))
end
@show maximum(energy_abs)
end
energy_check(Points)
# maximum(energy_abs) = 4.0548661650596785e-10
```

it’s seems to follow the system’s trend so reinitialize and find the poincare map

```
rotor_system = CoupledODEs(rotor_1, u0, params; diffeq,t0=0.0)
psos = poincaresos(rotor_system, (1,0.0), Tmax = 500000.0; u0=u0)
```

this gives me

```
julia> psos = poincaresos(rotor_system, (1,0.0), Tmax = 500000.0; u0=u0)
4-dimensional StateSpaceSet{Float64} with 1 points
NaN NaN NaN NaN
```

I have looked into poincaresos manpage Orbit Diagrams & PSOS · DynamicalSystems.jl

from what I see hyperplane is defined in this for plotting section of θ2,pθ2 whenever θ1 becomes 0.0

ofcourse If I look at the Points[:,1] it’s indeed crossing 0 so maybe poincaresos is not using callback defined, I don’t know how to proceed further. Is there any other way or package where I can do this?