// This does the 3 body 3-4-5 configuration.
// It reports ke, pe & te.
// This includes a variable dt and seems to work.
// I changed the centering scheme to match vdt.html.
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
typedef double ex;
#define N 3
static ex s(ex x){return x*x;}
typedef struct {ex x, y, xd, yd, m, mid;} pt;
//  xd & yd include factor of dt. m is mass. mid is dt^2/mass.

int main(){
  FILE * out = fopen ("trails", "wb");
  ex dt = .000001;
  pt z[] = {{0,4,0,0,3,0},{3,0,0,0,4,0},{0,0,0,0,5,0}}; // 345
  // {{-1,0,0,-.5,1,0},{1,0,0,.5,1,0}}; // Circle
  // {{0,0,sqrt(2),0,1,0},{0,0.5,-sqrt(2),0,1,0}}; // parabolas
  void bh(ex a){int j=N; 
    while(j--) {z[j].x += a*z[j].xd; z[j].y += a*z[j].yd;}}
  int n=0; ex t = 0; char slow, notfast;
  ex e(pt* j, pt* k){return - j->m*k->m/sqrt(s(j->x - k->x)+ s(j->y - k->y));}
  ex E(pt* j){return j->m*(s(j->xd/dt) + s(j->yd/dt))/2;}
  void pr(char * m){int k; bh(0.5);
    {ex pe = 0, ke = 0;
      {int j=N; while(j--) {int k=j; while (k--) pe += e(z+j, z+k); ke += E(z+j);}}
        printf("%s t=%12.7f, dt=%6.1e, te = %16.11f\n", m, t, dt, ke+pe);
        printf("n=%d pe = %13.7f ke = %13.7f\n", n, pe, ke);}
    for(k=0; k<N; ++k) {pt p = z[k];
        printf("%10.7f %10.7f %10.7f %10.7f\n", p.x, p.y, p.xd/dt, p.yd/dt);}
    bh(-0.5);}
  void conf(pt* j, pt* k){
    ex dx = j->x - k->x, dy = j->y - k->y, s2 = s(dx) + s(dy);
    ex d = j->m*k->m/(sqrt(s2)*s2), fx = dx*d, fy = dy*d;
    j->xd -= fx*j->mid; j->yd -= fy*j->mid; k->xd += fx*k->mid; k->yd += fy*k->mid;
    {ex b = (s(j->xd - k->xd) + s(j->yd - k->yd));
      if(s2 < 600000000*b) slow = 1; else if(s2 < 1500000000*b) notfast = 1;}}
  void step(){t += dt; slow=notfast=0;
     {int j=N; while(j--) {z[j].x += z[j].xd; z[j].y += z[j].yd;}}
     {int j=N; while(j--) {int k=j; while (k--) conf(z+j, z+k);}}
     {void adt(ex dtf){ex x = (1-dtf)/2, d2=s(dtf); dt *= dtf;
       {int j=N; while(j--){pt*q = z+j;
       q->x += x*q->xd; q->y += x*q->yd; q->xd *= dtf; q->yd *= dtf; q->mid *= d2;}}}
     if(slow) adt(.5); else if(!notfast) adt(2);}}
  void qt(int q){fclose(out); printf("qt %d\n", q); pr("fin"); exit(-3);}

  if(sizeof(z) != N*sizeof(z[0])) qt(4);
  {int j=N; while(j--) {z[j].xd *= dt; z[j].yd *= dt; z[j].mid = s(dt)/z[j].m;}}
  bh(-0.5); pr("init");
  {short v[2*N]; int ol=200000;
   short f(ex x) {if(x<-2.5) return 0; if(x<17.5) return 100*x + 250; return 999;}
   int cx(pt w, int j){short a = v[j], b = v[j+1], A = f(w.x), B = f(w.y);
             v[j] = A; v[j+1] = B; return a != A || b != B;}
   while(1) {step(); ++n;
      if(!(n%1000000)) pr(">");
      if(0) {int k=N; while(k--) if (cx(z[k], 2*k))
         {putc(k, out); fwrite(v+2*k, 4, 1, out); if(!--ol) qt(2);}}
          if(t>40) qt(1);}}}