I’m trying to plot an orbit as a surface in the extended phase space of the system

$ $ \dot{q}=p$ $ $ $ \dot{p}=-sin(q)(1+\epsilon sin(\theta))$ $ $ $ \dot{\theta}=\omega$ $

When $ \epsilon=0$ the dynamics can be portrayed on the $ p,q$ phase space (without the extension of the phase space)

` f[p_, q_, \[Theta]_] := p; g[p_, q_, \[Theta]_] := -Sin[ q] (1 + \[Epsilon] Sin[\[Theta]]) ; h[p_, q_, \[Theta]_] := \[Omega]; H[q_, p_] := 1/2 p^2 - Cos[q]; \[Delta] = 0; \[Epsilon] = 0.1; \[Omega] = 1; tf = 10; Show[ {ContourPlot[H[q, p] == 1, {q, -2 \[Pi], 2 \[Pi]}, {p, -6, 6}, Contours -> 100, ContourStyle -> {Red, Thick}], StreamPlot[{p, -Sin[q]}, {q, -2 \[Pi], 2 \[Pi]}, {p, -6, 6}, StreamPoints -> Fine, PlotRangePadding -> 0, ImageSize -> 500]}] `

where I’m interested in the dynamics on the red curve. d I wish the plot trajectories on the orbit (red curve) in the extended phase space such that $ \theta=mod(\theta,2\pi)$

When I attempt to plot a trajectory, I have a connecting lines between each section which was shifted with the $ mod$ function

` \[Delta] = 0; \[Epsilon] = 0.5; \[Omega] = 4; tf = 10; sol = NDSolve[{q'[t] == f[p[t], q[t], \[Theta][t]], p'[t] == g[p[t], q[t], \[Theta][t]], \[Theta]'[t] == h[p[t], q[t], \[Theta][t]], WhenEvent[Mod[\[Theta][t], 2 \[Pi]] == 0 , \[Theta][t] -> 0], p[0] == .1, q[0] == .1, \[Theta][0] == 0}, {p, q, \[Theta]}, {t, 0, tf}] Plot[Evaluate[\[Theta][t] /. sol], {t, 0, tf}] ParametricPlot3D[ Evaluate[{p[t], q[t], \[Theta][t]} /. sol], {t, 0, tf}, AspectRatio -> 1/2] `

at the end I would like to get a 3d plot of the surface with the trajectories that bounded to it, such as in the figure below