// This is intended to be the general N-body problem with prespecified time duration.
// It reports ke, pe & te.
// This includes a variable dt.
// This version is for the 64 sig bit Pentium (depending on the def of R).
// The Runge Kutta scheme.
// It is for D dimensions.
// This code has the dubious ability to track masses in 4D thru
// an inverse square law field.
// All of these tests report 0 as momentum, center of mass and angular momentum,
// bacause each point is either at the origin, or of negligible mass.
#include "h.h"
void prx(char * m, pt z[N]){
  R E(pt* j){R S=0; for(int m=0; m<D; ++m) S += s(j->dv[m].xd); return j->m*S/2;}
  R e(pt* j, pt* k){R S=0; for(int m=0; m<D; ++m) S += s(j->dv[m].x - k->dv[m].x); 
           return - j->m*k->m/sqrtl(S);}
 {R pe = 0, ke = 0;
    for(int j=0; j<N; ++j) {for(int k=0; k<j; ++k) pe += e(z+j, z+k); ke += E(z+j);}
    printf("%s pe=%16.12f ke=%16.12f te= %16.12f\n", m, (d)pe, (d)ke, (d)(ke+pe));}
 for(int k=0; k<N; ++k) {pt p = z[k]; for(int n=0; n<D; ++n)
          printf("%18.15Lf %18.15Lf ", p.dv[n].x, p.dv[n].xd); printf("\n");}}

void spin(pt z[N], R tlim, char draw, int phop){
  char nm[] = {'t','r','a','i','l','A',0};
  FILE * out = draw?fopen (nm, "wb"):0;
  R dt=0; int n=0; R t = 0;
  void ri(){R cg[D], M[D], L[D][D], tm=0;
     for(int j=0; j<D; ++j) {cg[j] = M[j] = 0; for(int k=0; k<D; ++k) L[j][k]=0;}
     for(int j=0; j<N; ++j) {tm += z[j].m;
        for(int k=0; k<D; ++k) {cg[k] += z[j].dv[k].x*z[j].m;
           M[k] += z[j].dv[k].xd*z[j].m;}
        for(int k=0; k<D; ++k) for(int l=0; l<D; ++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(R * tp, R f, char* s){
             for(int 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");
         for(int j=0; j<D; ++j) pDex(&L[j][0], 1, "L");}}

  void pr(char * m){if (phop < 0) return;
    printf("n=%d, t=%15.10f, dt=%7.2e\n", n, (d)t, (d)dt); prx(m, z);}

  int step(){typedef struct {Dv dv[N][D];} kt;
     {d dd = 1.e20; for(int j=0; j<N; ++j) for(int k=0; k<j; ++k)
       {d s2 = 0; for(int i=0; i<D; ++i) s2 += s(z[j].dv[i].x - z[k].dv[i].x);
          dd = dd<s2?dd:s2;}
       {d ds = sqrt(dd); dt = 0.0004*ds*sqrt(ds);}}
     char done = t+dt>tlim; if (done) dt = tlim-t;
     kt f(kt y) {kt nv; for(int n=0; n<N; ++n) for(int i=0; i<D; ++i){
       nv.dv[n][i].x = y.dv[n][i].xd; nv.dv[n][i].xd = 0;}
       for(int j=0; j<N; ++j) for(int k=0; k<j; ++k) {R dx[D]; R s2=0;
         for(int d=0; d<D; ++d) s2 += s(dx[d] = y.dv[j][d].x - y.dv[k][d].x);
         {R o = z[j].m*z[k].m/(sqrtl(s2)*s2);
         for(int d=0; d<D; ++d) {R F = dx[d]*o;
         nv.dv[j][d].xd -= F/z[j].m; nv.dv[k][d].xd += F/z[k].m;}}}
       return nv;}
     kt lc(R sc, kt x, kt y){kt rt; for(int j=0; j<N; ++j) for(int d=0; d<D; ++d) {
       rt.dv[j][d] = (Dv){sc*x.dv[j][d].x  + y.dv[j][d].x,
                          sc*x.dv[j][d].xd + y.dv[j][d].xd};}
       return rt;}
     kt k0; for(int j=0; j<N; ++j) for(int d=0; d<D; ++d) k0.dv[j][d] = z[j].dv[d];
     {kt k1 = f(k0), // for time t
        k2 = f(lc(dt/2,k1,k0)), // for time t+dt/2
        k3 = f(lc(dt/2,k2,k0)), // for time t+dt/2
        k4 = f(lc(dt, k3, k0)), // for time t+dt
        k5 = lc(1, lc(2,k2,k1), lc(2,k3,k4));
        for(int j=0; j<N; ++j) for(int e=0; e<D; ++e)
         {z[j].dv[e].x+=dt*k5.dv[j][e].x/6; z[j].dv[e].xd+=dt*k5.dv[j][e].xd/6;}
     t += dt;}
    return done;}

  void qt(int q){if(phop>0) printf("qt %d\n", q);
    if(draw) fclose(out);
    if(phop>0) {pr("fin"); ri();}}
  if(!out && draw) qt(6);
  pr("init"); if(phop>0) ri();
  {/*ML*/ typedef struct{short x, y;} px; int nct = -1; R tretic = 0;
   px v[N]; int ol=200000; R fst = 8;
   short pc(R x) {if(x<-2.5) return 0; if(x<17.5) return 100*x + 250; return 999;}
   int cx(px w, int j){px a = v[j]; v[j] = w; return a.x != w.x || a.y != w.y;}
   do {/* Loop */ if(n>nct) {pr(">"); nct += phop;}
      if(draw) {void ph(short s){putc(255&(s>>8), out); putc(255&s, out);}
         if(t > fst)
           {fclose(out); ++nm[5]; out = fopen (nm, "wb"); if(!out) qt(5); fst += 8;}
         for(int k=0; k<N; ++k) {if (cx((px){pc(z[k].dv[0].x), pc(z[k].dv[1].x)}, k))
           {putc(k, out); ph(v[k].x); ph(v[k].y); if(!--ol) qt(2);}}
         if(t>tretic) {for(int k=0; k<N; ++k) {putc(3, out); ph(v[k].x); ph(v[k].y);}
                       tretic +=1;}}
      ++n;} while(!step()); qt(7);}}