Thanks for sharing the information, it's just amazing (-.-)

These pictures are really great, but I think Arnol'd ( the A of KAM ) would shake his head, reading this thread. Especially because these pictures aren't tori, but maybe they are Poincare sections. I don't know.

A KAM tori is really easy to make in any programming language. Start with a unit area that maps to the surface of a torus, and identify each dimension with a phase space angle of uncoupled harmonic oscillators. Third and fourth phase space dimensions can be constant conserved energies or one conserved energy and one conserved angular momentum. Time evolution has constant angular velocity in both phase spaces, so the graph is just a line in the unit area. Adding in the topology you have a line that (probably) goes dense on the surface of a Torus (or in the unit area). With angular velocity ratio $\sqrt{2}$:

Torus = {r Cos[\[Phi]], r Sin[\[Phi]], .5 Sin[\[Theta]]} /. {r -> (1 - 0.5 Cos[\[Theta]])};
Trajectory = {r Cos[\[Phi]], r Sin[\[Phi]], .6 Sin[\[Theta]]
} /. {r -> (1 - 0.6 Cos[\[Theta]])} /. {\[Theta] -> Sqrt[2] \[Phi]};
Show[
 ParametricPlot3D[Trajectory, {\[Phi], 0, 100 Pi}],
 ParametricPlot3D[Torus, {\[Phi], 0, 2 Pi}, {\[Theta], 0, 2 Pi}],
 Boxed -> False,
 Axes -> False
 ]
Graphics[Point@Trajectory[[1 ;; 2]] /. \[Theta] -> # 2 Pi & /@  Range[100]]

1. You can add a small r^4 perturbation, but then you need to change parameterization slightly. The advantage of r^4 is that angular momentum and energy are still constants. Using r^4 perturbation, it's possible to see precession, which is the one remarkable difference from the multi-periodic, harmonic example.

2. If you add a perturbation with D4 symmetry, you lose conservation of momentum, but still have a three dimensional conserved energy surface, where I think it's possible for chaos to happen. Then you will get some pictures like the above by taking Poincare sections, and it will be more intelligible and less plagiarism. Note: I do not recommend triangular Henon Heiles, as odd powers potentials are more difficult to work with than strictly even.

So my suggestion is basically don't jump to stage 3 or 4 before you go through the basics. Ignore this advice at the risk of delusion, confusion, etc.

If you have never done phase space geometry before, my article Plane Pendulum and Beyond is a beginner introduction, where you can pick up some of the basic ideas of symplectic geometry without needing to worry about what Darboux's theorem looks like on a banquet table, and does it apply to dessert?

Best Regards.

I found another issue with the compiled function implementation: i had already assigned a numerical value to "angle" in my notebook and was not aware of it, but when you do a clean restart and this is not the case then the compiled function will not build correctly. In consequence the angle variable has also to be put in the interface to the compiled function. So below is a corrected version and I checked it with a clean restart.

f1 = Compile[{{xy, _Real, 1}}, 
   Module[{cs = Cos[xy[[3]]], sn = Sin[xy[[3]]], x = xy[[1]], 
     y = xy[[2]]}, {N[x cs + (x^2 - y) sn], 
     N[x sn - (x^2 - y) cs], xy[[3]]}]];

kamtorus2[ angle_, step_, ende_, points_] := 
 Module[{ x, xold, y, orbit},
  For[orbit = 0, orbit < ende, orbit += step,
   pts = AppendTo[pts, 
      NestWhileList[f1, {orbit/3., orbit/3., angle}, 
       Total[Abs[#]] < 2. 10^38 &, 1, points ]];
   ]]

pts = {};
Timing[kamtorus2[1.3, 0.05, 1.5, 200]]
ListPlot[Flatten[pts, 1][[All, 1 ;; 2]], AspectRatio -> 1]

Now it should compile and run smoothly and the output reacts to the setting of angle

Here is a more Mathematica style implementation (which runs faster):

kamtorus2[ angle_, step_, ende_, points_] := 
 Module[{ x, xold, y, orbit},
  For[orbit = 0, orbit < ende, orbit += step,
   pts = AppendTo[pts, 
      NestWhileList[f1, orbit/3. {1, 1}, Total[Abs[#]] < 2. 10^38 &, 
       1, points ]];
   ]]

pts = {};
Timing[kamtorus2[1.3, 0.05, 1.5, 200]]
ListPlot[Flatten[pts, 1], AspectRatio -> 1]

Below is a somewhat naive implementation.

kamtorus[ angle_, step_, ende_, points_] := 
 Module[{ x, xold, y, orbit},
  For[orbit = 0, orbit < ende, orbit += step,
   x = orbit/3.;
   y = x;
   ptno = 0;
   While[Abs[x] < 1. 10^38 && Abs[y] < 1. 10^38 &&  ptno < points,
    xold = x;
    x = N[ x        Cos[angle] +     (x^2 - y) Sin[angle]];
    y =  N[ xold Sin[angle] - (xold^2 - y) Cos[angle]];
    pts = AppendTo[pts, {x, y}];
    ptno++
    ];
   ]
  ]

pts = {};
Timing[kamtorus[1.3, 0.05, 1.5, 200]]
ListPlot[pts, AspectRatio -> 1]