// This includes a variable dt and seems to work.
// I changed the centering scheme to match vdt.html.
// This version is for the 64 sig bit Pentium.
// This has logic to compare two nearby paths thru phase space.
// The Runge Kutta scheme.
// Switching to D dimensions instead of 2.
// This code has the dubious ability to track masses in 4D thru
// an inverse square law field.
// Deal with 0 masses (test points)
// This hacked up version seems to show stability of L4 with mass ratio = 999.
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
typedef double ex; typedef double d;
// #define sqrt sqrtl
#define D 2
#define N 3
static ex s(ex x){return x*x;}
typedef struct{ex x, xd;} Dv;
typedef struct {Dv dv[D]; ex m;} pt;
typedef struct {pt pos[N]; ex t, dt;} ss;

static void spin(ss * gurg, ex tlim, ex del, char draw){
  char nm[] = {'t','r','a','i','l','A',0};
  FILE * out = draw?fopen (nm, "wb"):0;
  ex dt = .000001; int n=0; ex t = 0;
  pt z[] = {{{{.002,0},{0,.002}},7.992},
       {{{-1.998,0},{0,-1.998}},.008},
       {{{-.998,-sqrt(3)},{sqrt(3)+.001,-.998}},0}}; // Lagrange
  // {{{{del,0},{4,0}},3},{{{3,0},{0,0}},4},{{{0,0},{0,0}},5}}; //345
  // {{{{del,0},{4,0},{0,0.1}},3},{{{3,0},{0,0},{0,0}},4},{{{0,0},{0,0},{0,0}},5}}; // 345 3D
  // {{{{-1+del,0},{0,-.5}},1},{{{1,0},{0,.5}},1}}; // Circle
  // {{{{0,sqrt(2)},{del,0}},1},{{{0,-sqrt(2)},{0.5,0}},1}}; // parabolas
  void ri(){ex cg[D], M[D], L[D][D], tm=0;
     {int j=D; while(j--) {cg[j] = M[j] = 0; {int k=D; while(k--) L[j][k]=0;}}}
     {int j=N; while(j--){tm += z[j].m;
        {int k=D; while(k--) {cg[k] += z[j].dv[k].x*z[j].m;
           M[k] += z[j].dv[k].xd*z[j].m;}}
        {int k=D; while(k--){int l=D; while(l--) L[k][l] +=
             z[j].m*(z[j].dv[k].x*z[j].dv[l].xd -  z[j].dv[l].x*z[j].dv[k].xd);}}}}
         {void pDex(ex * tp, ex f, char* s){int j;
             for(j=0; j<D; ++j) printf("%20.14f", (d)(*(tp+j)*f));
             printf("<-%s\n", s);}
         pDex(cg, 1/tm, "cg"); pDex(M, 1, "M");
         {int j; for(j=0; j<D; ++j) pDex(&L[j][0], 1, "L");}}}
  ex e(pt* j, pt* k){ex S=0; int m=D; while(m--) S += s(j->dv[m].x - k->dv[m].x); 
           return - j->m*k->m/sqrt(S);}
  ex E(pt* j){ex S=0; int m=D; while(m--) S += s(j->dv[m].xd); return j->m*S/2;}
  void pr(char * m){int k;
    {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=%15.10f, dt=%6.1e, te = %14.11f\n", m, (d)t, (d)dt, (d)(ke+pe));
        printf("n=%d pe = %16.10f ke = %16.10f\n", n, (d)pe, (d)ke);}
    for(k=0; k<N; ++k) {pt p = z[k]; int n; for(n=0; n<D; ++n)
          printf("%13.10f %13.10f ", (d)p.dv[n].x, (d)p.dv[n].xd); printf("\n");}}
  void step(){typedef struct {Dv dv[N][D];} kt;
     char slow = 0, notfast=0;
     kt f(ex t, kt y, char dtc) {kt nv; {int n = N; while(n--) {int d=D; while(d--){
       nv.dv[n][d].x = y.dv[n][d].xd; nv.dv[n][d].xd = 0;}}}
       {int j=N; while(j--) {int k=j; while (k--) if((z[j].m != 0) || (z[k].m != 0))
       {ex dx[D]; ex s2=0;
         {int dd=D; while(dd--) s2 += s(dx[dd] = y.dv[j][dd].x - y.dv[k][dd].x);}
         {ex d = 1./(sqrt(s2)*s2);
         {int dd=D; while(dd--) {ex F = dx[dd]*d;
         nv.dv[j][dd].xd -= F*z[k].m; nv.dv[k][dd].xd += F*z[j].m;}}
         if(dtc){ex b = 0; {int dd=D; while(dd--) b +=s(z[j].dv[dd].xd - z[k].dv[dd].xd);}
           b*=s(dt); if(s2 < 6000000.*b) slow = 1; else if(s2 < 15000000.*b) notfast = 1;}}}}}
       return nv;}
  kt lc(ex sc, kt x, kt y){kt rt; {int j=N; while(j--) {int dd=D; while(dd--){
     rt.dv[j][dd] = (Dv){sc*x.dv[j][dd].x + y.dv[j][dd].x, sc*x.dv[j][dd].xd + y.dv[j][dd].xd};}}}
     return rt;}
  kt k0; {int j=N; while(j--) {int dd=D; while(dd--) k0.dv[j][dd] = z[j].dv[dd];}}
  {kt k1 = f(t, k0, 1),
      k2 = f(t + dt/2, lc(dt/2,k1,k0), 0),
      k3 = f(t + dt/2, lc(dt/2,k2,k0), 0),
      k4 = f(t + dt, lc(dt, k3, k0), 0),
      k5 = lc(1, lc(2,k2,k1), lc(2,k3,k4));
   {int j=N; while(j--) {int dd=D; while(dd--)
       {z[j].dv[dd].x += dt*k5.dv[j][dd].x/6; z[j].dv[dd].xd += dt*k5.dv[j][dd].xd/6;}}}
     t += dt; if(slow) dt *= .5; else if(!notfast) dt *= 2;}}
  void qt(int q, char r){
    if(draw) fclose(out); printf("qt %d\n", q); pr("fin"); ri(); if(r) exit(-3);}

  if(sizeof(z) != N*sizeof(z[0])) qt(4, 1); if(!out && draw) qt(6, 1);
  if(0) {float rnd(){return (float)rand()/RAND_MAX;}
     Dv rp(){float x, y, S; do {x = 2*rnd() - 1; y = 2*rnd() - 1;} while((S=x*x+y*y)>=1);
          {float L = sqrt(-2*log(S)/S); return (Dv){x*L, y*L};}}
     {int n=N; srand(45); while(n--) {int d = D; while(d--) z[n].dv[d] = rp();}}}
  pr("init"); ri();
  {ex cg[D], tm[D], tM=0; {int d = D; while(d--) cg[d] = tm[d] = 0;}
    {int n = N; while(n--) {int d = D; while(d--) {cg[d] += z[n].m*z[n].dv[d].x;
      tm[d] += z[n].m*z[n].dv[d].xd;} tM += z[n].m;}}
    {int n = N; while(n--) {int d = D; while(d--)
         {z[n].dv[d].x -= cg[d]/tM; z[n].dv[d].xd -= tm[d]/tM;}}}}
  pr("t0"); ri();
  {/*ML*/ typedef struct{short x, y;} px; int nct = -1; ex tretic = 0;
   px v[N]; int ol=200000; ex fst = 1000; ex cpnt = 1;
   float clip(float c){return c<-1?0:c<1?(c+1)*350:699;}
   float pcx(ex x) {return clip(30*x+30);}
   float pcy(ex x) {return clip(30*x-51.5);}
   {int n = N; while(n--) v[n] = (px){pcx(z[n].dv[0].x), pcy(z[n].dv[1].x)};}
   while(1) {/* Loop */ if(n>nct) {pr(">"); nct += 10000;}
      if(draw) {void ph(short s){putc(255&(s>>8), out); putc(255&s, out);}
         void shv(short * sc, float pq){float c = *sc;
            if (pq<c) --*sc; else if (pq>c+1.25) ++*sc;}
         int cx(float x, float y, int j){px a = v[j]; shv(&v[j].x, x); shv(&v[j].y, y);
            // if(j == 2) printf("xxs %e %e %d %d\n", (d)x, (d)y, v[j].x, v[j].y);
            return a.x != v[j].x || a.y != v[j].y;}
         if(t > fst)
           {fclose(out); ++nm[5]; out = fopen (nm, "wb"); if(!out) qt(5, 1); fst += 8;}
         {int k=N; float c = cos(t), s = sin(t);
         while(k--) {float x = z[k].dv[0].x, y = z[k].dv[1].x;
           if (cx(pcx(x*c + y*s), pcy(y*c - x*s), k))
           {putc(k, out); ph(v[k].x); ph(v[k].y); if(!--ol) qt(2, 1);}}}
         if(t>tretic) {int k=N; while(k--) {putc(3, out); ph(v[k].x); ph(v[k].y);}
             tretic +=2*3.1415926535;}}
      if(gurg && (t > cpnt)){ss grab(){ss s; {int j=N; while(j--) s.pos[j] = z[j];}
                s.t = t; s.dt = dt; return s;}
         ss tmp = grab(); dt = cpnt - t; step();
         *(gurg++) = grab(); cpnt += 1.;
         {int j=N; while(j--) z[j] = tmp.pos[j];} t = tmp.t; dt = tmp.dt;}
      if(t>tlim) {qt(7, 0); return;}
      step(); ++n;}}}
#define last 50
ss A[last], B[last];
int main(){
  if(0 /*MS*/){spin(A, last, 0, 1); spin(B, last, 1.e-11, 0);
  {int j; for(j=0; j<last; ++j) {ex mx = 0;
    {int k = N; while(k--){int dd=D; while(dd--) {
      mx += s(A[j].pos[k].dv[dd].x  - B[j].pos[k].dv[dd].x)
          + s(A[j].pos[k].dv[dd].xd - B[j].pos[k].dv[dd].xd);}}
    printf("%d %e\n", j, (d)sqrt(mx));}}}}
  else spin(0, 160, 0, 1);
  return 0;}