# This does the 3 body 3-4-5 configuration.
  This version is the simplest.
  This version updates x too often. #

MODE FL = LONG LONG REAL;
OP S = (FL x)FL: x UP 2;
OP SR = (REAL x)REAL: sqrt(x);
OP SR = (LONG REAL x)LONG REAL: long sqrt(x);
OP SR = (LONG LONG REAL x)LONG LONG REAL: long long sqrt(x);

MODE PT = STRUCT(FL x, y, xd, yd, m, mid);
[3]PT z := ((0,4,0,0,3,0),(3,0,0,0,4,0),(0,0,0,0,5,0));
FL dt = .000001; FL t := 0;
PROC pr = (STRING m)VOID: (print((t, " ", m));
   FOR k TO 3 DO PT p = z[k]; print(newline);
      print((x OF p, y OF p, xd OF p/dt, yd OF p/dt)) OD;
      print(newline));
PROC conf = (REF PT j, k)VOID: (FL dx = x OF j - x OF k, dy = y OF j - y OF k,
    s2 = S dx + S dy, d = m OF j*m OF k/((SR s2)*s2), fx = dx*d, fy = dy*d;
    xd OF j -:= fx* mid OF j; yd OF j -:= fy* mid OF j;
    xd OF k -:= fx* mid OF k; yd OF k -:= fy* mid OF k;
    x OF j +:= xd OF j; y OF j +:= yd OF j;
    x OF k +:= xd OF k; y OF k +:= yd OF k);
PROC step = VOID:(conf(z[1], z[2]); conf(z[1], z[3]); conf(z[2], z[3]); t +:= dt);
FOR j TO 3 DO mid OF z[j] := S dt / m OF z[j] OD;
  pr("init"); step; pr("once");
  TO 10 DO TO 1000 DO step OD; pr("Chug") OD