// galaga.c (starble galaxy)
// simulation of n-body galaxy
// with initial random velocity distribution
// and collisions like marbles with u=.1
// will add other models later
#include <stdio.h>
#include <conio.h>
#include <dos.h>
#include <mem.h>
#include <stdlib.h>
#include <math.h>
#include <values.h>
#include <time.h>
#include <process.h>

#define X_OFFS  160.0
#define Y_OFFS  100.0
#define X_RES   320
#define Y_RES   200
#define X_RESM  319
#define Y_RESM  199
#define VIDEO 0xA000
#define PUTPIXS(x,y,color) canvas[x+((y)<<8)+((y)<<6)] = color
#define GETPIXS(x,y) canvas[x+((y)<<8)+((y)<<6)]
#define M_PIO2 1.57079632679489661923
#define M_PIT2 6.28318530717958647692
#define LOOPSPTICK 18.0 // central pick of loops per clock tick for option to scale moves speeds by LOOPSPTICK/prevloopcts
#define GRAV   3.50e-6 //some gravitational constant
#define MAXNSTARS 100

unsigned char far *SCREEN=MK_FP(VIDEO,0);
unsigned char far *canvas;

typedef struct
{
  float x;
  float y;
  float z;
} vector ;

typedef struct
{
  vector r; //position
  vector v; //velocity
  float m; //mass
  float rad; //radius
  vector w; //angular velocity
  int c; //color
} star ;

star stars[MAXNSTARS], oldstars[MAXNSTARS], saved[MAXNSTARS];
// oldstars is for erase, with some wasted fields

int nstars=20;
int savedstars;
float zthr=50.0;

typedef struct
{
  int x;
  int y;
  int c;
} pixel ;

void vga_mode(int mode)
{
  union REGS regs;

  regs.x.ax = mode;
  int86(0x10, &regs, &regs);
}

void clear_screen(void)
{
  _fmemset(canvas, 0, (unsigned)64000);
}

void vgapalset(void)
{
  int i,j,k;
  int vgapal[256][3]; //custom vga palette

  {
  k=0;
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[i+16*j][j]= (4*i<64) ? (4*i) : 63;
    vgapal[i+16*j][(j+1)%3]=0;
    vgapal[i+16*j][(j+2)%3]=0;}
   k+=16*3; //48
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[k+i+16*j][j]= (4*i<64) ? (4*i) : 63;
    vgapal[k+i+16*j][(j+1)%3]= (4*i<64) ? (4*i) : 63;
    vgapal[k+i+16*j][(j+2)%3]=0;}
   k+=16*3; //96
  for (i=0;i<16;i++) for (j=0;j<3;j++)
    vgapal[k+i][j]= (4*i<64) ? (4*i) : 63;
   k+=16;   //112
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[k+i+16*j][j]= (4*i<64) ? (4*i) : 63;
    vgapal[k+i+16*j][(j+1)%3]= 2*i;
    vgapal[k+i+16*j][(j+2)%3]= 2*i;}
   k+=16*3; //164
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[k+i+16*j][j]= (4*i<64) ? (4*i) : 63;
    vgapal[k+i+16*j][(j+1)%3]= (4*i<64) ? (4*i) : 63;
    vgapal[k+i+16*j][(j+2)%3]= 2*i;}
   k+=16*3; //208
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[k+i+16*j][j]= (4*i<64) ? (4*i) : 63;
    vgapal[k+i+16*j][(j+1)%3]= 1.3*i;
    vgapal[k+i+16*j][(j+2)%3]= 1.3*i;}
   k+=16*3; //256
  for(j=0;j<3;j++) {vgapal[0][j]=0;}
  }
  //outp(0x03C6,0xff);
  outp(0x03c8, 0);
  for(i=0;i<256;i++)
  {
    for(j=0;j<3;j++)
    {
    outp(0x03c9, vgapal[i][j]);
    }
  }
}

void pauser(int *pause)
{
  while (!kbhit()) ;
  getch();
  *pause=0;
}

void steppauser(int *steppause)
{
  int ch;
  while (!kbhit());
  ch=getch();
  if(ch!=' ') *steppause=0;
}

int yres(int y)
{
  if (y<0) return 0;
  if (y>=200) return 199;
  return y;
}

int imod( int a, int b )
//integer mod, no negative returns
{
  a = a % b ; 
  if (a<0) a += b;
  return a;
}

float twopimod( float angle )
{
  float answer;
  answer = fmod( angle , M_PIT2 );
  if ( answer<0.0 ) answer += M_PIT2;
  return answer;
}

float pimod( float angle )
{
  float answer;
  answer = fmod( angle , M_PI );
  if ( answer<0.0 ) answer += M_PI;
  return answer;
}

float nearfmodmult( float quant, float grain )
{ //quantizes quantity quant to nearest multiple of grain 
  //this should work, but test this function
  float flt;
  flt = fmod(quant,grain) ;
  if(flt==0.0) return quant;
  else
  {
    if(quant>0.0)
    {
      if ( 2.0*flt>=grain )           // fabs(grain) ) ; assume positive grain
        return quant+grain-flt;
      else
        return quant-flt;
    }
    if(quant<0.0)
    {
      if ( -2.0*flt>=grain )           // fabs(grain) ) ; assume positive grain
        return quant-grain-flt;
      else
        return quant-flt;
    }
  }
}

pixel screenproj(vector r)
{
// projection to screen
//#define ZTHR 50
#define XF 0.833333333333
  pixel p;
  if (r.z<zthr)
  {
    p.y = Y_OFFS - r.y/(zthr-r.z) * Y_OFFS;
    p.x = r.x / (zthr-r.z) * X_OFFS*XF + X_OFFS;
    //p.c = 8*zthr / (zthr-r.z)  - 15
    p.c = 8.0/zthr * (r.z-zthr);
    if(p.x<0)p.x=0; else if(p.x>X_RESM)p.x=X_RESM;
    if(p.y<0)p.y=0; else if(p.y>Y_RESM)p.y=Y_RESM;
    if(p.c<-15)p.c=-15; 
  }
  else
  { p.x=-1; p.y=-1; p.c=0; }
  return p;
}

vector vecttimes(float a, vector b)
{
  b.x *= a ;
  b.y *= a ;
  b.z *= a ;
  return b ;
}

vector vectsum(vector a, vector b)
{
  a.x += b.x ;
  a.y += b.y ;
  a.z += b.z ;
  return a ;
}

vector vectdiff(vector a, vector b)
{
  a.x -= b.x ;
  a.y -= b.y ;
  a.z -= b.z ;
  return a ;
}

float vectdot(vector a, vector b)
{
  a.x  = a.x*b.x;
  a.x += a.y*b.y;
  a.x += a.z*b.z;
  return a.x ;
}

vector vectcross(vector a, vector b)
{
  vector c ;
  c.x = a.y * b.z - b.y * a.z ;
  c.y = a.z * b.x - b.z * a.x ;
  c.z = a.x * b.y - b.x * a.y ;
  return c ;
}

float vabs2(vector a)
{
  return vectdot(a,a);
}

float vabs(vector a)
{
  return sqrt(vectdot(a,a));
}

vector unitvect(vector b)
{
  float flt;
  flt = vectdot(b,b);
  if (flt!=0.0)
    return vecttimes(1.0/sqrt(flt),b);
  else
  {
    b.x=b.y=b.z=0.0;
    return b;
  }
}

vector vectproj(vector a, vector b)
{ //returns vector of projection of a onto b
  b = unitvect(b);
  return vecttimes( vectdot(a,b), b) ;
}

float metric(vector a, vector b)
// distance between two points
{
  vector r;
  r = vectdiff(a,b);
  return sqrt(vectdot(r,r));
}

void centerdists(void)
{
  int i;
  vector ar={0.0,0.0,0.0},av={0.0,0.0,0.0};
  for (i=0;i<nstars;i++)
  {
    ar.x+=stars[i].r.x;
    ar.y+=stars[i].r.y;
    ar.z+=stars[i].r.z;
    av.x+=stars[i].v.x;
    av.y+=stars[i].v.y;
    av.z+=stars[i].v.z;
  }  
  ar.x/=nstars;
  ar.y/=nstars;
  ar.z/=nstars;
  av.x/=nstars;
  av.y/=nstars;
  av.z/=nstars;
  for (i=0;i<nstars;i++)
  {
    stars[i].r.x -= ar.x;
    stars[i].r.y -= ar.y;
    stars[i].r.z -= ar.z;
    stars[i].v.x -= av.x;
    stars[i].v.y -= av.y;
    stars[i].v.z -= av.z;
  }  
}

void init_stars()
{
//  vector r; //position
//  vector v; //velocity
//  float m; //mass
//  float rad; //radius
//  vector w; //angular velocity
//  int c; //color
// stars[nstars];
#define RBOX 20.0
#define VBOX 0.01
#define WBOX 0.001
  int i;
  for (i=0;i<nstars;i++)
  {
    stars[i].r.x = rand()*2.0*RBOX/RAND_MAX-RBOX;
    stars[i].r.y = rand()*2.0*RBOX/RAND_MAX-RBOX;
    stars[i].r.z = rand()*2.0*RBOX/RAND_MAX-RBOX;
    stars[i].v.x = rand()*2.0*VBOX/RAND_MAX-VBOX;
    stars[i].v.y = rand()*2.0*VBOX/RAND_MAX-VBOX;
    stars[i].v.z = rand()*2.0*VBOX/RAND_MAX-VBOX;
    stars[i].m   = 0.08;
    stars[i].rad = 1.0;
    stars[i].w.x = rand()*2.0*WBOX/RAND_MAX-WBOX;
    stars[i].w.y = rand()*2.0*WBOX/RAND_MAX-WBOX;
    stars[i].w.z = rand()*2.0*WBOX/RAND_MAX-WBOX;
    stars[i].c   = (rand()%6)*16+15;
  }
  centerdists();
}

void starserasecopy(void)
{
  int i;
  for (i=0;i<nstars;i++)
    oldstars[i]=stars[i];
}

void savecopy(void)
{
  int i;
  for (i=0;i<nstars;i++)
    saved[i]=stars[i];
  savedstars=nstars;
}

void restoresaved(void)
{
  int i;
  nstars=savedstars;
  for (i=0;i<nstars;i++)
    stars[i]=saved[i];
}

void bresenline(int x0, int y0, int x1, int y1, int c)
{
  int steep;
  int deltax;
  int deltay;
  int error;
  int ystep;
  int x,y,t;
  steep = ( abs(y1 - y0) > abs(x1 - x0) ) ;
  if (steep) {
    t=x0; x0=y0; y0=t;
    t=x1; x1=y1; y1=t;  }
  if (x0 > x1) {
    t=x0; x0=x1; x1=t;
    t=y0; y0=y1; y1=t;  }
  deltax = x1 - x0;
  deltay = abs(y1 - y0);
  error = 0;
  y = y0;
  if (y0 < y1) ystep = 1; else ystep = -1;
  for(x=x0;x<=x1;x++)
  { 
    if (steep) { PUTPIXS(y,x,c); }
    else { PUTPIXS(x,y,c); }
    error = error + deltay;
    if ( error<<1 >= deltax )  {
      y = y + ystep;
      error = error - deltax; }
  }
}

void erasestar(int i)
// 
{
//  vector r; //position
//  vector rad; //radius
//  oldstars[nstars];
  pixel p;
  vector sp,spp;
  float ra=oldstars[i].rad;
  sp = oldstars[i].r;
  spp = sp; 
  spp.z += ra;
  p = screenproj( spp );
  if(p.x==-1)return;
  PUTPIXS(p.x,p.y,0);
  spp = sp;
  spp.x += ra;
  p = screenproj( spp );
  PUTPIXS(p.x,p.y,0);
  spp = sp;
  spp.x -= ra;
  p = screenproj( spp );
  PUTPIXS(p.x,p.y,0);
  spp = sp; 
  spp.y += ra;
  p = screenproj( spp );
  PUTPIXS(p.x,p.y,0);
  spp = sp; 
  spp.y -= ra;
  p = screenproj( spp );
  PUTPIXS(p.x,p.y,0);
}

void drawstar(int i, int erase)
// 
{
//  vector r; //position
//  vector rad; //radius
//  int c; //color
//  stars[nstars];
  pixel p;
  vector sp,spp;
  float ra=stars[i].rad;
  sp = stars[i].r;
  spp = sp; 
  spp.z += ra;
  p = screenproj( spp );
  if(p.x==-1)return;
  p.c += stars[i].c;
  if(erase) p.c=0;
  PUTPIXS(p.x,p.y,p.c);
  spp = sp;
  spp.x += ra;
  p = screenproj( spp );
  p.c += stars[i].c;
  if(erase) p.c=0;
  PUTPIXS(p.x,p.y,p.c);
  spp = sp;
  spp.x -= ra;
  p = screenproj( spp );
  p.c += stars[i].c;
  if(erase) p.c=0;
  PUTPIXS(p.x,p.y,p.c);
  spp = sp; 
  spp.y += ra;
  p = screenproj( spp );
  p.c += stars[i].c;
  if(erase) p.c=0;
  PUTPIXS(p.x,p.y,p.c);
  spp = sp; 
  spp.y -= ra;
  p = screenproj( spp );
  p.c += stars[i].c;
  if(erase) p.c=0;
  PUTPIXS(p.x,p.y,p.c);
}

void erasestars()
{
  int i;
  for(i=0;i<nstars;i++)
    erasestar(i);
}

void drawstars(int erase)
{
  int i;
  for(i=0;i<nstars;i++)
    drawstar(i,erase);
}

void movestar(int s, float dt)
{
  stars[s].r = vectsum( stars[s].r, vecttimes( dt, stars[s].v ) ) ;
  // what if they also had some reference central angle / surface position to move by stars[s].w?
  // I could take the surface position and move it by surface velocity w X (surface point) with centripetal accelaration;
  // is there a better way to do that?
  // I could rotate the surface position by w dt rotation matrices.  What about non-commutativity?
  // Could have w dt shift euler angles, and that would be more of a specification of orientation than one surface point.
  // Quaternions?  Yeah, quaternions.
}

void movestars(float dt)
{
  int i;
  for(i=0;i<nstars;i++)
    movestar(i,dt);
}

float ballmoment(float m, float r)
{
  return 0.4*m*r*r;
}

int autocls( clock_t tock, int autoclear )
{
  static int lastclear=0;
  if (lastclear==0) lastclear = tock;
  if (tock<lastclear) lastclear = tock;
  if ( tock > lastclear + (64>>(autoclear-1)) ) //~4 seconds /2^(autoclear-1)
    { clear_screen(); lastclear = tock; return 1;}
  else return 0;
}

/*void p1p2pushout4(float dist)
{  // from sphrspwa
  float th1,th2; //collision line bearings, theta
  float halfovrlp,midovrlp1,midovrlp2;
  float p12radii;
  float loni,lati;
  float fraction;

  p12radii = (p1rad+p2rad)*M_PI/Y_RES;
  halfovrlp = ( p12radii - dist ) *0.5 ;
  midovrlp1 = ( p12radii - dist )*p2rad/(p1rad+p2rad) ;
  midovrlp2 = p12radii - dist - midovrlp1 ;

  fraction = ( p1rad*M_PI/Y_RES - halfovrlp ) / dist ;
  intermediatepoint(p1long, p1lat, p2long, p2lat, fraction, &loni, &lati);
  th1 = M_PI + angtofrom(loni, lati, p1long, p1lat) ;
  th2 = M_PI + angtofrom(loni, lati, p2long, p2lat) ;

  surfacerun(p1long,p1lat,th1,midovrlp1,&p1long,&p1lat);
  surfacerun(p2long,p2lat,th2,midovrlp2,&p2long,&p2lat);
}*/

void starblebounce(int j, int k, float uf)
//collision response treating stars as marbles
{
  vector disp;
  float dist, flt;
  vector cl; //collision line
  float cldot1, cldot2;
  star s1=stars[j],s2=stars[k];
  vector v1=s1.v, v2=s2.v;
  vector v1icl,v2icl; // collision line initial velocities
  float  v1iclm,v2iclm; // collision line initial velocities magnitudes
  vector v1fcl,v2fcl; // collision line final velocities
  float  v1fclm,v2fclm; // collision line final velocities magnitudes
  vector v1ici,v2ici; // collision transverse initial velocities
  vector v1fci,v2fci; // collision transverse final velocities
  float m1=s1.m, m2=s2.m;
  disp = vectdiff(s2.r,s1.r);
  dist = vabs(disp);  if(dist==0.0) return;
  cl = vecttimes( 1.0/dist , disp ); // collision line unit vector from 1 to 2
  v1iclm = vectdot( v1 , cl );
  v2iclm = vectdot( v2 , cl );
  if (v1iclm-v2iclm<0.0) return; //not colliding
  v1icl = vecttimes( v1iclm , cl );
  v2icl = vecttimes( v2iclm , cl );
  v1ici = vectdiff( v1 , v1icl );
  v2ici = vectdiff( v2 , v2icl );
  v1fclm = ( (m1-m2) * v1iclm  +  2.0 * m2 * v2iclm ) / ( m1 + m2 ) ;
  v2fclm = ( (m2-m1) * v2iclm  +  2.0 * m1 * v1iclm ) / ( m1 + m2 ) ;
  v1fcl = vecttimes( v1fclm , cl );
  v2fcl = vecttimes( v2fclm , cl );
  if ( uf!=0.0 )
  {
    vector cp; //collision point position
    float rad1=s1.rad, rad2=s2.rad;
    vector v1cp,v2cp; // collision-transverse velocities at the contact point
    vector rv1cp; // contact point relative collision-transverse velocity
    vector rv1cpuv,rv2cpuv; // '' unit vectors
    vector w1i=s1.w, w2i=s2.w, w1f, w2f;
    float cldv1,cldv2; //collision velocity delta for friction normal
    float i1,i2; // moments of inertia
    int overfric1,overfric2;
    float td1=1.0, td2=1.0;
    vector cpr1,cpr2; //contact point from ball centers
    cp = vectsum( vecttimes( 0.5*rad1 - 0.5*rad2 + 0.5*dist , cl ), s1.r ); //
    i1 = ballmoment(m1,rad1);
    i2 = ballmoment(m2,rad2);
     cpr1 = vectdiff(cp,s1.r);  //
     cpr2 = vectdiff(cp,s2.r);
     //cpr1 = vecttimes( rad1, unitvect( vectdiff(cp,s1.r) ) );  //
     //cpr2 = vecttimes( rad2, unitvect( vectdiff(cp,s2.r) ) );
    cldv1 = vabs(vectdiff(v1fcl,v1icl));
    cldv2 = vabs(vectdiff(v2fcl,v2icl));
    v1cp = vectsum( v1ici, vectcross( w1i, cpr1 ) );
    v2cp = vectsum( v2ici, vectcross( w2i, cpr2 ) );
    rv1cp = vectdiff( v1cp, v2cp );
    //rv2cp = vectdiff( v2cp, v1cp );
    rv1cpuv = unitvect(rv1cp);
    rv2cpuv = vecttimes(-1.0,rv1cpuv);
    v1fci = vectsum( v1ici, vecttimes( -uf*cldv1, rv1cpuv ) );
    v2fci = vectsum( v2ici, vecttimes( -uf*cldv2, rv2cpuv ) );
    w1f = vectsum( w1i, vectcross( cpr1, vecttimes( -uf*cldv1*m1/i1, rv1cpuv ) ) );
    w2f = vectsum( w2i, vectcross( cpr2, vecttimes( -uf*cldv2*m2/i2, rv2cpuv ) ) );
    if( vectdot( vectsum( v1fci, vectcross( w1f, cpr1) ), rv1cpuv)
        - vectdot( v2cp, rv1cpuv) < 0.0 ) overfric1=1; else overfric1=0;
    if( vectdot( vectsum( v2fci, vectcross( w2f, cpr2) ), rv2cpuv)
        - vectdot( v1cp, rv2cpuv) < 0.0 ) overfric2=1; else overfric2=0;
    if( overfric1 || overfric2 )
    {
      if ( overfric1 )
        td1 = vectdot(rv1cp,rv1cpuv) / ( uf * cldv1 * (1.0 + m1*vabs2(cpr1)/i1) );
      if ( overfric2 )
        td2 = vectdot(rv1cp,rv1cpuv) / ( uf * cldv2 * (1.0 + m2*vabs2(cpr2)/i2) );
      if(td1>1.0 || td2>1.0){ gotoxy(20,1); printf("%4f %4f",td1,td2); }//I hope this never happens
      if ( td2<td1 ) td1=td2 ;
      v1fci = vectsum( v1ici, vecttimes( -uf*cldv1*td1, rv1cpuv ) );
      v2fci = vectsum( v2ici, vecttimes( -uf*cldv2*td1, rv2cpuv ) );
      w1f = vectsum( w1i, vectcross( cpr1, vecttimes( -uf*cldv1*m1/i1*td1, rv1cpuv ) ) );
      w2f = vectsum( w2i, vectcross( cpr2, vecttimes( -uf*cldv2*m2/i2*td1, rv2cpuv ) ) );
    }
    if ( m1*vabs2(v1fci)+i1*vabs2(w1f)+m2*vabs2(v2fci)+i2*vabs2(w2f)>m1*vabs2(v1ici)+i1*vabs2(w1i)+m2*vabs2(v2ici)+i2*vabs2(w2i) )
    {
      float tlke,a,b;
      b = -2.0*m1*uf*cldv1*td1*vectdot(v1ici,rv1cpuv) -2.0*vabs(cpr1)*uf*cldv1*m1*td1*vectdot(w1i,rv1cpuv)
          -2.0*m2*uf*cldv2*td1*vectdot(v2ici,rv2cpuv) -2.0*vabs(cpr2)*uf*cldv2*m2*td1*vectdot(w2i,rv2cpuv);
      a =   m1*pow(uf*cldv1*td1,2) + i1*pow(vabs(cpr1)*uf*cldv1*m1/i1*td1,2)
          + m2*pow(uf*cldv2*td2,2) + i2*pow(vabs(cpr2)*uf*cldv2*m2/i2*td2,2);
      tlke = -b/a;
      if(tlke<0.0)tlke=0.0;
      v1fci = vectsum( v1ici, vecttimes( -uf*cldv1*td1*tlke, rv1cpuv ) );
      v2fci = vectsum( v2ici, vecttimes( -uf*cldv2*td1*tlke, rv2cpuv ) );
      w1f = vectsum( w1i, vectcross( cpr1, vecttimes( -uf*cldv1*m1/i1*td1*tlke, rv1cpuv ) ) );
      w2f = vectsum( w2i, vectcross( cpr2, vecttimes( -uf*cldv2*m2/i2*td1*tlke, rv2cpuv ) ) );
    }
    stars[j].w=w1f;
    stars[k].w=w2f;
    stars[j].v=vectsum(v1fcl,v1fci);
    stars[k].v=vectsum(v2fcl,v2fci);
  }
  else
  {
    stars[j].v=vectsum(v1fcl,v1ici);
    stars[k].v=vectsum(v2fcl,v2ici);
  }
}

/*
void forwardafterbackout(,,float a, float dt)
{ //forward moves after collision bounce after back-out of pre-collision penetration
}

float backoutprecollpenet(,,, float dt)
{ //backs out the pre-collision penetration 

}
*/

void starblecollisions(float dt)
{
  //static int colliding[((nstars-1)*nstars)/2];
  float radii;
  int j,k;
  for (j=0;j<nstars-1;j++)
    for (k=j+1;k<nstars;k++)
    {
      radii = stars[j].rad + stars[k].rad ;
      if ( vabs2(vectdiff(stars[k].r,stars[j].r)) < radii*radii )
      {
        //float backout; 
        //backout=backoutprecollpenet(,,,dist,dt)
        starblebounce(j,k,0.2);
        //forwardafterbackout(,,float a, float dt)
      }
    }
//  if (!colliding[(j+1)*(nstars-1 + nstars-1-j)/2 - nstars+1+j  +   k-j-1])
//if (!colliding[(j+1)*(nstars+nstars-2-j)/2 - nstars + k])
//nstars=6         1                      5      6      1..5    0..4
//                 2                      9      6      2..5    5..8
//                 3                      12     6      3..5    9..11
//  if (!colliding[(j)*(nstars+nstars-1-j)/2  +         k-j-1])
//nstars=6         0                      0             0..4    0..4
//                 1                      5             0..3    5..8
//                 2                      9             0..2    9..11
//
//what about sectorizing the space somehow and checking for collisions only for starbles in the same sector?
//
//what about multiple collisions and
// backing up time to the before the penetrating motion
// or only doing the collision bounce for the objects that have passed into each other the most for objects with more than one collision
//? need to think about this more
}

void starbreak(void)
{
// take a star and shoot it at 0,0,0
  static int corn=1;
  stars[0].r.x = -2.0*corn*RBOX;
  stars[0].r.y = -2.0*corn*RBOX;
  stars[0].r.z = 0;
  stars[0].v.x = 0.5*corn*VBOX;
  stars[0].v.y = 0.5*corn*VBOX;
  stars[0].v.z = 0;
  corn*=-1;
}

void starbreak2(void)
{
// make a star and shoot it at 0,0,0
  static int corn=1;
  if (nstars==MAXNSTARS) return;
  stars[nstars].r.x = -2.0*corn*RBOX;
  stars[nstars].r.y = -2.0*corn*RBOX;
  stars[nstars].r.z = 0;
  stars[nstars].v.x = 0.5*corn*VBOX;
  stars[nstars].v.y = 0.5*corn*VBOX;
  stars[nstars].v.z = 0;
  corn*=-1;
  stars[nstars].m   = 0.08;
  stars[nstars].rad = 1.0;
  stars[nstars].w.x = rand()*2.0*WBOX/RAND_MAX-WBOX;
  stars[nstars].w.y = rand()*2.0*WBOX/RAND_MAX-WBOX;
  stars[nstars].w.z = rand()*2.0*WBOX/RAND_MAX-WBOX;
  stars[nstars].c   = (rand()%6)*16+15;
  nstars+=1;
}

void orbitshot(void)
{
// make a star and shoot it at 0,0,something, possibly into an orbit
  static int corn=1;
  if (nstars==MAXNSTARS) return;
  stars[nstars].r.x = -2.0*corn*RBOX;
  stars[nstars].r.y = -2.0*corn*RBOX;
  stars[nstars].r.z = corn*RBOX;
  stars[nstars].v.x = nstars/10.0*corn*VBOX;
  stars[nstars].v.y = nstars/10.0*corn*VBOX;
  stars[nstars].v.z = 0;
  corn*=-1;
  stars[nstars].m   = 0.08;
  stars[nstars].rad = 1.0;
  stars[nstars].w.x = rand()*2.0*WBOX/RAND_MAX-WBOX;
  stars[nstars].w.y = rand()*2.0*WBOX/RAND_MAX-WBOX;
  stars[nstars].w.z = rand()*2.0*WBOX/RAND_MAX-WBOX;
  stars[nstars].c   = (rand()%6)*16+15;
  nstars+=1;
}

void gravity(int j, int k, float dt)
{
  vector fg;
  vector disp;
  float dist2;
  vector v1=stars[j].v, v2=stars[k].v;
  float m1=stars[j].m, m2=stars[k].m;
  disp = vectdiff(stars[k].r,stars[j].r);
  dist2 = vectdot(disp,disp);
  fg = vecttimes( GRAV*m1*m2/dist2, vecttimes( sqrt(dist2), disp ) ) ;
  v1 = vectsum( v1, vecttimes( dt/m1, fg ) );
  v2 = vectdiff(v2, vecttimes( dt/m2, fg ) );
  stars[j].v=v1;
  stars[k].v=v2;
}

void gravitate(float dt)
{
  int j,k;
  for (j=0;j<nstars-1;j++)
    for (k=j+1;k<nstars;k++)
      gravity(j,k,dt);
}

float gravitypot(int j, int k)
{  //measure potential energy for physdisplay
  float dist;
  float m1=stars[j].m, m2=stars[k].m;
  dist = vabs( vectdiff(stars[k].r,stars[j].r) );
  return -GRAV*m1*m2/dist ;
}

float gravitatepot(void)
{  //measure potential energy for physdisplay
  int j,k;
  float gpe=0.0;
  for (j=0;j<nstars-1;j++)
    for (k=j+1;k<nstars;k++)
      gpe += gravitypot(j,k);
  return gpe;
}

void physdisplay(int d)
{
  int i;
  vector tp={0.0,0.0,0.0},tl={0.0,0.0,0.0};
  float te=0.0;
  for (i=0;i<nstars;i++)
  {
    vector p;
    float m;
    m = stars[i].m;
    p = vecttimes( m, stars[i].v);
    tp = vectsum( tp , p );
    te += 0.5*vectdot(p,p)/m;
    tl = vectsum( tl , vectsum( vecttimes( ballmoment(m,stars[i].rad),stars[i].w ), vecttimes(m,vectcross(stars[i].r,p) ) ) );
  }
  if(d==2)
  {  
    te += gravitatepot();
    gotoxy(1,23);
    printf("te: %+9.5f",te);
  }
  gotoxy(1,24);
  printf("tp: %+9.5f %+9.5f %+9.5f",tp.x,tp.y,tp.z);
  gotoxy(1,25);
  printf("tl: %+9.5f %+9.5f %+9.5f",tl.x,tl.y,tl.z);
}

int main(void)
{
  int quit=0;
  int ch, i;
  int delnum=0;
  int delnummem=1;
  int intextscr=0;
  int pause=0;
  int autoclear=0;
  int erasetrails=1;
  int kollisions=1;
  unsigned long int loopctr=0;// loop counts so far // at 300 loops per second, would take almost half a year to flip
  unsigned long int prevtickloopct=0;
  int prevloopcts=16; // how many loops in the last clock tick
  int smoothedloopcts[16]={16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16}; // smoothed previous loop counts; earlier method (15* + new)/16 seemed odd
  int smoothedlci=1; //index into smoothedloopcts[]
  clock_t start, tick, tock, gtick, keytick, dptick;
  int physdisp=0;
  int blankphysdisp=0;
  int fpsreadout=0;
  int steppause=0;
  int movestiming=0;  // option to scale moves speeds by LOOPSPTICK/prevloopcts
  unsigned long int physdisploopct=0;

  start = clock();
  dptick = keytick = gtick = clock();

  canvas = SCREEN;
  vga_mode(0x13);
  vgapalset();
  clear_screen();

  init_stars();
  starserasecopy();
  savecopy();

  tick = clock();
  while (!quit)
  {
    delay(delnum);

    loopctr++;
    if ( ( tock=clock() ) > tick ) // how long does it take for clock to flip a u long integer? 7 years.
    {
      prevloopcts = loopctr - prevtickloopct;  // (int)(loopctr - prevtickloopct);
      ch=0; for(i=1;i<16;i++) ch+=smoothedloopcts[i];
      smoothedloopcts[0] = (ch+prevloopcts)>>4; smoothedloopcts[smoothedlci]=prevloopcts;
      smoothedlci = ( smoothedlci + 1 )%15 + 1;
      prevtickloopct = loopctr;
      tick = tock;
    } //system clock() ticks ~18.2 times a second

    if (pause) { pauser(&pause); }

    if(kbhit())
    {
      ch=getch();
      switch( ch )
      {
        case 0x1b:
        { 
          if(!intextscr) quit=1; 
          else
          {
            intextscr = 0;
            vga_mode(0x13);
            vgapalset();
          }
        }
        break;
        case '-':
        case '`': clear_screen();
        break;
        case '~': vga_mode(0x13);  vgapalset();
        break;
        case ' ': pause=1;
        break;
        case '|': steppause=1; //hit ' ' to step in steppause
        break;
        case '%': delnum++;
        break;
        case '5': {if (delnum!=0) {delnummem=delnum; delnum=0;}
                      else { delnum=delnummem; }  }
        break;
        case 'T': 
        case 't': erasetrails=!erasetrails;
        break;
        case 'A': zthr*=0.99; if(zthr<1.0)zthr=1.0;
        break;
        case 'a': zthr+=3.0;
        break;
        case '[': physdisp=(physdisp+1)%3;
        break;
        case 'k': kollisions=!kollisions;
        break;
        case 'b': starbreak(); //remove a starble to a distance
        break;
        case 'B': starbreak2(); //drop in another starble
        break;
        case 'o': orbitshot(); //throw a starble
        break;
        case 'r': nstars--; if(nstars==0)nstars=1;
        break;
        case 'c': centerdists(); //recenter the galaga
        break;
        case 'H': savecopy();
        break;
        case 'h': restoresaved();
        break;
        case '\\': fpsreadout=!fpsreadout;
        break;
        case '=':
        case 'I': nstars=20;
        case 'i': init_stars(); //random distribution of velocities and positions, centered
        break;
      } // end switch( ch )
    }

    if( autoclear && autocls(tock,autoclear) && !intextscr ) gtick=0;

    if ( !pause )
    {
/*      float dt=1.0, flt;  int dts=1, mtdts=1, h;
      flt=;
      for(i=0;i<MAX_SHOTS;i++)
      {
        if(shots[i].owner!=-1)
        {
          flt=max(flt,shots[i].v);
        }
      }
      if( flt>2.0*vinc ) { dts=ceil(flt/vinc); dt=1.0/dts; }
      if( movestiming )
      {
        flt = (float)LOOPSPTICK/prevloopcts;
        mtdts=ceil(flt); dt*=flt/mtdts;
      }
      //{ gotoxy(30,1); printf("%3d %3d",mtdts,dts); }
      for (h=0;h<mtdts;h++)
      for (i=0;i<dts;i++)  */
      {
        float dt=5.0;
        gravitate(dt);
        movestars(dt);
        if(kollisions)starblecollisions(dt);
      }
    }

    if(!intextscr && tock>gtick)
    { 
      gtick=tock;
      if(erasetrails) 
      {
        erasestars();
        starserasecopy();
      }
      drawstars(0);
    }
    if( fpsreadout ) { gotoxy(2,1); printf("%4d",smoothedloopcts[0]*18); } //*18.2); }
    if( physdisp && (loopctr-physdisploopct>prevloopcts || pause) )
    {
      physdisploopct = loopctr;
      physdisplay(physdisp);
    }
    if(steppause) {steppauser(&steppause);}    //
  }

  tock = clock();
  if(!intextscr) vga_mode(0x03);
  printf("Time elapsed: %.1f seconds.\n", (tock - start) / CLK_TCK);
  printf("In 1/18.2 second clock ticks: %d.\n", (tock - start) );
  pauser(&pause);
  return 0;
}