// This does the 3 body 3-4-5 configuration.
// It reports ke, pe & te.
// This includes a variable dt and seems to work roughly.
// I changed the centering scheme to match vdt.html and now its broke!
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 3
static double s(double x){return x*x;}
typedef struct {double 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");
  double 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(double a){int j=N; 
    while(j--) {z[j].x += a*z[j].xd; z[j].y += a*z[j].yd;}}
  double t = 0;
  char slow, notfast;
  double e(pt* j, pt* k){return - j->m*k->m/
              sqrt(s((j->x) - (k->x))+ s((j->y) - (k->y)));}
  double E(pt* j){return j->m*(s(j->xd/dt) + s(j->yd/dt))/2;}
  int n=0;
  void pr(char * m){int k; bh(0.5);
    {double 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){
   double dx = j->x - k->x, dy = j->y - k->y, s2 = s(dx) + s(dy);
   double 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;
   {double b = (s(j->xd - k->xd) + s(j->yd - k->yd));
     if(s2 < 6000000*b) slow = 1;
     else if(s2 < 15000000*b) notfast = 1;}}
   void ap(pt * j){j->x += j->xd; j->y += j->yd;}
   void step(){t += dt; slow=notfast=0;
      {int j=N; while(j--) ap(z+j);}
      {int j=N; while(j--) {int k=j; while (k--) conf(z+j, z+k);}}
      {void adt(double dtf){double 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 t=%12.8f\n", q, t); 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=75000;
   short f(double 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%10000)) pr(">");
      {int k, n = 0; for(k=0; k<N; ++k) if (cx(z[k], 2*k)) n=1;
         if(n) {fwrite(v, 4, N, out); if(!--ol) qt(2); if(t>93) qt(1);}}}}}
  // Events at (-12.8166666) 1.8794, (-12.8166666) 15.829, (-12.775013) 29.88, (-12.77494017)