使用Mathematica绘制双摆系统的Poincaré section

双摆模型的方程为（系数做过一定的伸展之后）：

目标：绘制当alpha=0和alpha'>0时的Poincare截面（由点(beta,l_beta)组成）。

abc = {a'[t] == 2*(la[t] - (1 + Cos[b[t]])*lb[t])/(3 - Cos[2*b[t]]), 
   la'[t] == -2*Sin[a[t]] - Sin[a[t] + b[t]], 
   b'[t] == 
    2*(-(1 + Cos[b[t]])*la[t] + (3 + 2*Cos[b[t]])*lb[t])/(3 - 
        Cos[2*b[t]]), 
   lb'[t] == -Sin[a[t] + b[t]] - 
     2*Sin[b[t]]*((la[t] - lb[t])*lb[t])/(3 - Cos[2*b[t]]) + 
     2*Sin[2*b[t]]*(la[t]^2 - 
         2*(1 + Cos[b[t]])*la[t]*lb[t] + (3 + 2*Cos[b[t]])*
          lb[t]^2)/(3 - Cos[2*b[t]])^2};

psect[{a0_, la0_, b0_, lb0_}] := 
 Reap[NDSolve[{abc, a[0] == a0, la[0] == la0, b[0] == b0, 
     lb[0] == lb0, 
     WhenEvent[a[t] == 0 && -2*Sin[a[t]] - Sin[a[t] + b[t]] < 0, 
      Sow[{b[t], lb[t]}]], 
     WhenEvent[a[t] == 0 && -2*Sin[a[t]] - Sin[a[t] + b[t]] > 0, 
      Sow[{b[t], lb[t]}]]}, {a[t], b[t], la[t], lb[t]}, {t, 0, 1000}, 
    MaxSteps -> Infinity]][[-1, 1]]

abcdata = 
  Map[psect, {{0.01, 0.01, 0.01, 0}, {0.01, 0.01, 0.01, 0.01}, {0.01, 
     0.01, 0.01, -0.01}}];

ListPlot[abcdata, ImageSize -> Medium, AspectRatio -> 1.15]

ps = psect /@ Table[{0.01, 0.005, 0.01, i}, {i, -.01, .01, .00125}];
ListPlot[ps, ImageSize -> Medium,AspectRatio->1.1]