How to create an interactive Poincare section from the given code?

Question

Here we discussed Poincare section for perturbed string model taken from the paper Chaotic dynamics of a suspended string in a gravitational background with magnetic field. We try to reproduce some data from Chaos of Wilson Loop from String Motion near Black Hole Horizon.
Code for the Poincare section of a chaotic system is given by

    value = {22.50340313, 60.54615508, 15.70200865, 6.12469956, 12.24939912, -1.68535, 12.4194}; {K1, K2, K3, K4, K5, \[Omega]sq[0], \[Omega]sq[1]} = Rationalize[value, 10^(-30)];
(* The Lagrangian of the system *)
lagrangian = Sum[Derivative[1][c[n]][t]^2 - c[n][t]^2*\[Omega]sq[n], {n, {0, 1}}] + K1*c[0][t]^3 + K2*c[0][t]*c[1][t]^2 + K3*c[0][t]*Derivative[1][c[0]][t]^2 + 
    K4*c[0][t]*Derivative[1][c[1]][t]^2 + K5*Derivative[1][c[0]][t]*c[1][t]*Derivative[1][c[1]][t]; 
c[0][t_] := OverTilde[c][0][t] + \[Alpha]1*OverTilde[c][0][t]^2 + \[Alpha]2*OverTilde[c][1][t]^2; 
c[1][t_] := OverTilde[c][1][t] + \[Alpha]3*OverTilde[c][0][t]*OverTilde[c][1][t]; 
\[Alpha]1 = -2; \[Alpha]2 = -2^(-1); \[Alpha]3 = -1; 
n = Expand[lagrangian]; 
vars = {OverTilde[c][0][t], OverTilde[c][1][t], Derivative[1][OverTilde[c][0]][t], Derivative[1][OverTilde[c][1]][t]}; 
lagrangian = Normal[Series[n /. Thread[vars -> m*vars], {m, 0, 3}]] /. m -> 1; 
momentum[n_] := D[lagrangian, Derivative[1][OverTilde[c][n]][t]]


(* Getting Hamiltonian from the given Lagrangian )

 hamiltonian = Expand[Sum[momentum[n]Derivative[1][OverTilde[c][n]][t], {n, {0, 1}}] - lagrangian];
(* Finding the equations of motion *)
    eulerLagrange[lagrangian_, vars_, dvars_] := Thread[Table[D[D[lagrangian, dvar], t], {dvar, dvars}] - Table[D[lagrangian, var], {var, vars}] == ConstantArray[0, Length[vars]]]; 
    equationsOfMotion = eulerLagrange[lagrangian, {OverTilde[c][0][t], OverTilde[c][1][t]}, {Derivative[1][OverTilde[c][0]][t], Derivative[1][OverTilde[c][1]][t]}];
(* Finding the Poincare Section *)
    sol = Table[Block[{a, b, [Chi], d, habd}, habd = hamiltonian /. t -> 0 /. {OverTilde[c][0][0] -> a, Derivative[1][OverTilde[c][0]][0] -> b, OverTilde[c][1][0] -> [Chi], 
             Derivative[1][OverTilde[c][1]][0] -> d}; {a, b, [Chi]} = {x, 0.0001, 0.0001}; d = d /. FindRoot[habd == 10^(-5), {d, -x}]; 
          Reap[NDSolve[{equationsOfMotion, OverTilde[c][0][0] == a, Derivative[1][OverTilde[c][0]][0] == b, OverTilde[c][1][0] == [Chi], Derivative[1][OverTilde[c][1]][0] == d, 
             WhenEvent[OverTilde[c][1][t] == 0 && Derivative[1][OverTilde[c][1]][t] >= 0, Sow[{OverTilde[c][0][t], Derivative[1][OverTilde[c][0]][t]}]]}, 
            {OverTilde[c][0][t], OverTilde[c][1][t]}, {t, 0, 8000}, Method -> {"FixedStep", Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 10, 
                "ImplicitSolver" -> {"Newton", AccuracyGoal -> MachinePrecision, PrecisionGoal -> MachinePrecision, "IterationSafetyFactor" -> 1}}}, StartingStepSize -> 1/10]]], 
        {x, -0.002, -0.0001, 0.00005}]; 
    ListPlot[Table[Flatten[sol[[i]][[2]], 1], {i, Length[sol]}], PlotTheme -> "Detailed", PlotRange -> {{-0.002, 0}, {-0.002, 0.002}}, 
      FrameLabel -> {OverTilde[Subscript[c, 0]], OverDot[OverTilde[Subscript[c, 0]]]}, PlotStyle -> PointSize[0.005]]

How can I make it interactive using the ClickPoincarePlot2D resource function available in Wolfram function repository?

PS. The result of the execution of the code is as follows (Few warnings are omitted.). user64494

Answer 1

To display Poincare section using ResourceFunction["ClickPoincarePlot2D"] we need to define Hamiltonian and equations of motion in the Hamilton form as follows

h := Function[{S}, -((33707 S[[1]]^2)/20000) - (
   815338710539 S[[1]]^3)/51728115000 + S[[2]]^2 + (
   154040173 S[[1]] S[[2]]^2)/20000000 + (62097 S[[3]]^2)/5000 - (
   20799664898903 S[[1]] S[[3]]^2)/248503740000 + (
   103117489 S[[2]] S[[3]] S[[4]])/12500000 + S[[4]]^2 + (
   103117489 S[[1]] S[[4]]^2)/25000000]
eqs = {Derivative[1][x][t] == px[t], 
   Derivative[1][px][t] == -((154040173 px[t]^2)/20000000) - (
     103117489 py[t]^2)/25000000 + (33707 x[t])/10000 + (
     815338710539 x[t]^2)/17242705000 + (20799664898903 y[t]^2)/
     248503740000, Derivative[1][y][t] == py[t], 
   Derivative[1][py][t] == -((103117489 px[t] py[t])/12500000) - (
     62097 y[t])/2500 + (20799664898903 x[t] y[t])/124251870000};

Then evaluate

ResourceFunction["ClickPoincarePlot2D"][eqs, h, 10^-5, t, 6000, 
 y[t], {x[t], px[t]}, {PlotTheme -> "Detailed", 
  PlotRange -> {{-0.002, 0}, {-0.002, 0.002}}, Frame -> True, 
  FrameLabel -> {OverTilde[Subscript[c, 0]], 
    OverDot[OverTilde[Subscript[c, 0]]]}, 
  PlotStyle -> PointSize[0.005]}]

Appendix 1. To compute expression h from initial code we use

value = {22.50340313, 60.54615508, 15.70200865, 6.12469956, 
  12.24939912, -1.68535, 12.4194}; {K1, K2, K3, K4, 
  K5, \[Omega]sq[0], \[Omega]sq[1]} = Rationalize[value, 10^(-30)];
(The Lagrangian of the system)
lagrangian = 
  Sum[Derivative[1][c[n]][t]^2 - 
     c[n][t]^2[Omega]sq[n], {n, {0, 1}}] + K1c[0][t]^3 + 
   K2c[0][t]c[1][t]^2 + K3c[0][t]Derivative[1][c[0]][t]^2 + 
   K4c[0][t]Derivative[1][c[1]][t]^2 + 
   K5Derivative[1][c[0]][t]c[1][t]Derivative[1][c[1]][t];
c[0][t_] := 
  OverTilde[c][0][t] + [Alpha]1OverTilde[c][0][t]^2 + [Alpha]2*
    OverTilde[c][1][t]^2;
c[1][t_] := 
  OverTilde[c][1][t] + [Alpha]3OverTilde[c][0][t]OverTilde[c][1][t];
[Alpha]1 = -2; [Alpha]2 = -2^(-1); [Alpha]3 = -1;
n = Expand[lagrangian];
vars = {OverTilde[c][0][t], OverTilde[c][1][t], 
   Derivative[1][OverTilde[c][0]][t], 
   Derivative[1][OverTilde[c][1]][t]};
lagrangian = 
  Normal[Series[n /. Thread[vars -> m*vars], {m, 0, 3}]] /. m -> 1;
momentum[n_] := D[lagrangian, Derivative[1][OverTilde[c][n]][t]]
(Getting Hamiltonian from the given Lagrangian)
hamiltonian = 
  Expand[Sum[
     momentum[n]*Derivative[1][OverTilde[c][n]][t], {n, {0, 1}}] - 
    lagrangian];
rules = {OverTilde[c][0][t] -> S[[1]], Derivative[1][OverTilde[c][0]][t] -> S[[2]], OverTilde[c][1][t] -> S[[3]], 
    Derivative[1][OverTilde[c][1]][t] -> S[[4]]}; hamiltonian /. rules
(-((33707 S[[1]]^2)/20000) - (815338710539 S[[1]]^3)/51728115000 + S[[2]]^2 + (
 154040173 S[[1]] S[[2]]^2)/20000000 + (62097 S[[3]]^2)/5000 - (
 20799664898903 S[[1]] S[[3]]^2)/248503740000 + (
 103117489 S[[2]] S[[3]] S[[4]])/12500000 + S[[4]]^2 + (
 103117489 S[[1]] S[[4]]^2)/25000000 )

Last expression used in definition h:=Function[{S},...]. Equations of motion

rules1 = {OverTilde[c][0][t] -> x[t], Derivative[1][OverTilde[c][0]][t] -> px[t], OverTilde[c][1][t] -> y[t], 
    Derivative[1][OverTilde[c][1]][t] ->py[t]};H= hamiltonian /. rules1; 
eqs = {x'[t] == px[t], px'[t] == -D[H, x[t]], y'[t] == py[t],
   py'[t] == -D[H, y[t]]}
Out[]= {Derivative[1][x][t] == px[t], 
 Derivative[1][px][t] == -((154040173 px[t]^2)/20000000) - (
   103117489 py[t]^2)/25000000 + (33707 x[t])/10000 + (
   815338710539 x[t]^2)/17242705000 + (20799664898903 y[t]^2)/
   248503740000, Derivative[1][y][t] == py[t], 
 Derivative[1][py][t] == -((103117489 px[t] py[t])/12500000) - (
   62097 y[t])/2500 + (20799664898903 x[t] y[t])/124251870000}