/*
  IMPACT - Virtual Reality Workshop
    First implementation by John Henckel on Oct 7, 1993

  Acknowledgements:
  * This program was developed using Borland Turbo C++ 3.0.  The speed and
    versatility of the graphics of this program is a tribute to the talented
    programmers at Borland.  Also their online help is fabulous.
  * The equations for collision of irregularly shaped objects are taken from
    "Dynamics" by Pestel and Thomson, 1968, I haven't found such a detailed
    analysis of this problem in any other book.
*/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <graphics.h>
#include <conio.h>
#include "mouse.h"

// Macros and constants

#define MA 10       // array size
#define MB 50       // bodies
#define MX 4        // virtual x dimension size (origin lower-left)
#define MY 3        // virtual y size
#define EPS 0.00001 // small number
#define PI 3.1415926
#define frand(x) (x/1000.0*random(1000))
#define HOME "\x1B[;H"         // ansi escape code
#define rline(a,b,c,d) line((a)*xm/MX,ym-(b)*ym/MY,(c)*xm/MX,ym-(d)*ym/MY)

// Types

typedef float fa[MA];
typedef struct { int x,y; } pt;

typedef struct {     // Body:       (- means scratch field)
  int color;         //  color of the body
  float r;           //  radius (maximum)
  float m;           //  mass
  float i;           //  coefficient of moment of inertia ( = I/m )
  int n;             //  number of points in the shape (0=ball,1=??)
  fa sx,sy;          //  shape (center of mass is always origin)
  fa px,py;          //- corner points
  fa len;            //  length of each edge (from i to i+1)
  pt poly[2][MA+1];  //- polygons for draw and erase
  int pp;            //- current poly 0 or 1
  float vx,vy,w;     //  velocity
  float cx,cy,a;     //  center position
  float ocx,ocy,oa;  //- old center position
} body;

// Global Variables

int xm,ym;
body bod[MB];        // Array of bodies
int numb;

// Global options

float fric=0.0; // air friction
float grav=0.0;   // gravity
float scut=1.0;   // collision center adjustment
float rest=1.0;   // coefficient of restitution
int trace=0;      // erase image on redraw
int clear=0;      // clear screen
int mous=0;       // mouse present?
int delet=0;      // delete last body
int debug=0;
float covel=0;      // maximum contact velocity
float pene=0;
/*----------------------------------------------------------------------
   Redraw_all - erase all bodies and redraw them in the new position
*/
int redraw_all(int p)
{
  int b;

  if (clear) cleardevice();
  clear=0;
  for (b=0; b<numb; ++b) {
    setcolor(0);                         // black
    if (!trace) drawpoly(bod[b].n+1,(int*)bod[b].poly[1-p]);
    if (b<numb-delet) {
      setcolor(bod[b].color);
      drawpoly(bod[b].n+1,(int*)bod[b].poly[p]);
    }
  }
  numb -= delet;
  if (numb<0) numb=0;
  delet=0;

  setcolor(1);
  line(0,0,xm,0);
  line(xm,0,xm,ym);
  line(xm,ym,0,ym);
  line(0,ym,0,0);

  return 0;
}

/*----------------------------------------------------------------------
   computes new c,a  +=  new v,w
*/
int compute_new(void)
{
  int b;
  for (b=0; b<numb; ++b) {
    bod[b].cx += bod[b].vx;
    bod[b].cy += bod[b].vy;
    bod[b].a  += bod[b].w;
  }
  return 0;
}
/*----------------------------------------------------------------------
   sets px,py (corners) based on shape and cx,cy (center), a (angle)
*/
int compute_pxy(void)
{
  int i,b;
  float cx,cy,ca,sa;

  for (b=0; b<numb; ++b) {
    cx = bod[b].cx;
    cy = bod[b].cy;
    ca = cos(bod[b].a);
    sa = sin(bod[b].a);

    for (i=0; i<bod[b].n; ++i) {
      bod[b].px[i] = cx + bod[b].sx[i]*ca - bod[b].sy[i]*sa;
      bod[b].py[i] = cy + bod[b].sx[i]*sa + bod[b].sy[i]*ca;
    }
  }
  return 0;
}

/*----------------------------------------------------------------------
   set_old  - copy new c,a,v,w to old.
*/
int set_old(void)
{
  int b;
  for (b=0; b<numb; ++b) {
    bod[b].ocx = bod[b].cx; 
    bod[b].ocy = bod[b].cy; 
    bod[b].oa  = bod[b].a ; 
  }
  return 0;
}
/*----------------------------------------------------------------------
   copy px,py to poly[p],  adjust for scale.
   the range of pxy is [MX,MY], poly is [xm,ym] y inverted.
*/
int set_poly(int p)
{
  int i,b;

  for (b=0; b<numb; ++b) {
    for (i=0; i<bod[b].n; ++i) {
      bod[b].poly[p][i].x = bod[b].px[i]*xm/MX;
      bod[b].poly[p][i].y = (MY-bod[b].py[i])*ym/MY;
    }
    bod[b].poly[p][i].x = bod[b].poly[p][0].x;      // close the polygon
    bod[b].poly[p][i].y = bod[b].poly[p][0].y;
  }
  return 0;
}

/*----------------------------------------------------------------------
   draw mouse pointer
*/
int draw_mouse(int x, int y, int c)
{
  setcolor(c);
  line(x+5,y+5,x-5,y-5);
  line(x-5,y+5,x+5,y-5);
  return 0;
}


/*----------------------------------------------------------------------
   help - display help
*/
int help(void)
{
  setcolor(15);
  outtextxy(50,100,"This is a bouncing bodies simulation. ");
  outtextxy(50,110,"Version 1.1 (12/93) by John Henckel, "
                   "henck@rchvmv.vnet.ibm.com");
  outtextxy(50,120,"All comments and suggestions are welcome!");
  outtextxy(50,130,"Press a to add a new bodies");
  outtextxy(50,140,"Press d to delete the last body");
  outtextxy(50,150,"Press g to change gravitation constant [0.0]");
  outtextxy(50,160,"Press f to change air friction [0.0]");
  outtextxy(50,170,"Press r to change restitution coefficient [1.0]");
  outtextxy(50,180,"Press t to toggle trace mode [0]");
  outtextxy(50,190,"Press x to toggle debug mode [0]");
  outtextxy(50,200,"Press e for an energy report");
  outtextxy(50,210,"Press ? for help");
  outtextxy(50,220,"Press c to clear screen");
  outtextxy(50,230,"Press s to change body center adjustment factor [1]");
  outtextxy(50,240,"[] = default");
  if (mous) {
    outtextxy(50,250,"You can use the mouse to pull things.");
    outtextxy(50,260,"  button 1=pull, button 2=stop spin, both=pull ALL!");
  }
  outtextxy(50,280,"Press a key to continue");
  while (!kbhit());
  return clear=1;
}

/*----------------------------------------------------------------------
   msg - display a message on the message line
*/
int msg(int c,char *s)
{
  outtextxy(50,20+c*10,s);
  return clear=1;
}
int msgd(int c,char *s,int d)
{
  char x[80];
  sprintf(x,"%-30s %d",s,d);
  return msg(c,x);
}
int msgf(int c,char *s,float d)
{
  char x[80];
  sprintf(x,"%-30s %9.7f",s,d);
  return msg(c,x);
}


/*----------------------------------------------------------------------
   get_input  and redraw mouse pointer
*/
#define VV       0.1
int get_input(void)
{
  static int b=-1,x=0,y=0;
  static char s[320];
  int i,j,k,c,n;
  float d,t,x1,y1,a,m;

  if (!numb || b==-1) help();        // first time
  if (mous) {
    draw_mouse(x,y,0);
    b=get_mouse(&x,&y);
    draw_mouse(x,y,15-b);
  } else b=0;
  if (b) {
    x1 = 1.0*x/xm*MX;
    y1 = (1-1.0*y/ym)*MY;
    if (b<3) {
      d = 99999;                //  find closest body
      for (j=i=0; i<numb; ++i) {
        t = hypot(x1-bod[i].cx,y1-bod[i].cy);
        if (t < d) { d=t; j=i; }
      }
      if (d > EPS) {
        bod[j].vx = (x1-bod[j].cx)*VV;
        bod[j].vy = (y1-bod[j].cy)*VV;
      }
      if (b>1) bod[j].w = 0;           // stop spin
    }
    else {                                // attract everything
      for (j=0; j<numb; ++j) {
        bod[j].vx = (x1-bod[j].cx)*VV*0.1;
        bod[j].vy = (y1-bod[j].cy)*VV*0.1;
      }
    }
  }
  if (kbhit()) {
    while (kbhit()) b=getch();      // skip over typeahead and extended codes
    if (b==27 || b=='q') return 1;
    if (b=='g') {
      msg(0,"Enter gravitation constant*10000 (-10.0 to 10.0)");
      scanf("%f",&grav);
      grav /= 10000;
    }
    if (b=='f') {
      msg(0,"Enter air friction constant (0.0-100.0)");
      scanf("%f",&fric);
    }
    if (b=='r') {
      msg(0,"Enter coefficient of restitution (0.0 - 1.0)");
      scanf("%f",&rest);
      if (rest<-1) rest=-1;
      if (rest>2)  rest=2;
    }
    if (b=='v') {
      msg(0,"Enter minimum collision velocity (0.0 - 0.1)");
      scanf("%f",&covel);
    }
    if (b=='p') {
      msg(0,"Enter penetration acceleration*10000 (0.0 - 0.1)");
      scanf("%f",&pene);
      pene /= 10000;
    }
    if (b=='e') {
      t=a=d=m=0;
      for (i=0; i<numb; ++i) {
        t += 0.5*bod[i].m*(bod[i].vx*bod[i].vx + bod[i].vy*bod[i].vy);
        a += 0.5*bod[i].m*bod[i].i*bod[i].w*bod[i].w;
        d += bod[i].m*bod[i].cy*grav;
      }
      m=t+a+d;
      msgd(0,"Total number of bodies",numb);
      msgf(1,"Total translation energy",t);
      msgf(2,"Total rotation energy",a);
      msgf(3,"Total potential energy",d);
      msgf(4,"Total energy",m);
      msgf(6,"Gravitation",grav);
      msgf(7,"Friction",fric);
      msgf(8,"Restitution coefficient",rest);
      msgf(9,"Collision adjustment",scut);
      msgd(10,"Mouse detected",mous);
      msgd(11,"Debug mode",debug);
      while (!kbhit());
    }
    if (b=='s') {
      msg(0,"Enter collision center adjustment (0.0 - 1.0)");
      scanf("%f",&scut);
      if (scut<-1) scut=-1;
      if (scut>2)  scut=2;
    }
    if (b=='d') delet=1;
    if (b=='a' && numb<MB) {      // Add new body
      msg(0,"Enter data: Count[1] Sides[4] Aspect[1] Radius[.1] Color[0?]");
      msg(1,"            Noise[0] Mass[0] Moment[0] Pos[.5 .5 0] Vel[0 0 0]");
      msg(2,"[] = default.  Default mass and moment is based on radius, "
          "unit density.");
      msg(3,"For example, for five rectangles type: 5 4 2");
      n=1; j=4; d=.1; x1=1; c=0; y1=0.0001; m=0; a=0;   // defaults
      memset(&bod[numb],0,sizeof *bod);
      bod[numb].cx = bod[numb].cy = 0.5;
      do { gets(s); } while(*s < ' ');
      sscanf(s,"%d %d %f %f %d %f %f %f %f %f %f %f %f %f",
             &n,&j,&x1,&d,&c,&y1,&m,&a,
             &bod[numb].cx,&bod[numb].cy,&bod[numb].a,
             &bod[numb].vx,&bod[numb].vy,&bod[numb].w);
      if (d<EPS) d=0.1;
      if (d>2) d=2;
      if (m==0) m=PI*d*d*x1*100;
      if (a==0) a=d*d*x1;
      if (j<2) j=2;
      if (j>=MA) j=MA-1;
      if (n<0) n=0;
      if (n+numb > MB) n=MB-numb;
      if (y1>1) y1=1;

      // Set up the first new body
      bod[numb].color = c;
      bod[numb].r = d;
      bod[numb].m = m;
      bod[numb].i = a;

      for (i=numb; i<numb+n; ++i) {
        bod[i] = bod[numb];
        if (c==0) bod[i].color=9+random(6);
        if (i>numb)
          bod[i].cy = bod[i-1].cy + 2*d;
        if (bod[i].cy+d > MY) {
          bod[i].cy = d;
          bod[i].cx = bod[i-1].cx + 2*d*x1;
        }
        bod[i].n = j;
        for (k=0; k<j; ++k) {
          bod[i].sx[k] = d*(1-frand(y1))*cos(PI*2*(k+0.5+frand(y1))/j)*-x1;
          bod[i].sy[k] = d*(1-frand(y1))*sin(PI*2*(k+0.5+frand(y1))/j);
        }
        for (k=0; k<j; ++k) {        // wrong
          bod[i].len[k] = hypot(bod[i].sx[k],bod[i].sy[k]);
          if (bod[i].len[k] < EPS) bod[i].len[k] = EPS;
        }
        bod[i].w += frand(PI*y1);
      }
      numb += n;
    }
    if (b=='t') trace = 1-trace;
    if (b=='c') clear = 1;
    if (b=='x') debug = 1;
    if (b=='?') help();
  }
  return 0;
}
/*********      
      //-----------------------------------------------
      //  calculate, rx ry, velocity of i1 relative to b2
      x1 = px-bod[b1].ocx;
      y1 = py-bod[b1].ocy;
      x2 = px-bod[b2].ocx;
      y2 = py-bod[b2].ocy;
      rx = bod[b1].vx - bod[b2].vx - bod[b1].w*y1 + bod[b2].w*y2;
      ry = bod[b1].vy - bod[b2].vy + bod[b1].w*x1 - bod[b2].w*x2;
      //--------------------------------------------------
      //  calculate the previous location of i1, relative to b2
      qx = px - rx;
      qy = py - ry;
      //-----------------------------------------------------
      //  Find the edge of b2 farthest from q (using signed distance).
      //   (This is where we assume all bodies are convex.)
      i2=-1; d2=-99999;
      for (j=0; j<bod[b2].n; ++j) {
        d = dptl(qx,qy,b2,j);
        if (d>d2) { i2=j; d2=d; }
      }
      d1 = -dptl(px,py,b2,i2);
      if ((d2+d1) < EPS) d2 = d1+EPS;      // slow moving body!
      t = d1 / (d2+d1);    // t=how long ago did impact occur
      //---------------------------------------------------------
      //   Find the first impact that actually occurred
******/

/*----------------------------------------------------------------------
   distance from a point to a line (the edge of another body).
   result is negative inside (clockwise).
*/
float dptl(float x,float y,int b,int i)
{
  int t;
  float nx,ny;
  t = i+1;
  if (t==bod[b].n) t=0;
  nx = bod[b].py[i] - bod[b].py[t];
  ny = bod[b].px[t] - bod[b].px[i];    // n is normal vector
  x -= bod[b].px[i];
  y -= bod[b].py[i];
  nx = (x*nx + y*ny) / bod[b].len[i];
  return nx;
}

/*-----------------------------------------------------------------
   collision - test for a collision between two bodies, if so,
   return the point of impact, and unit vector normal to surface
   of impact (points into b2), and penetration distance.
*/
int collision(int b1,int b2,
              int *B1,int *B2,
              float *hx,float *hy,
              float *nx,float *ny,float *pen)
{
  float d1,d2,d,px,py;
  int i,j,i1,i2,t=0;

  d2 = -1;
  while (++t < 3) {                // repeat only two times
    for (i=0; i<bod[b1].n; ++i) {
      px = bod[b1].px[i];
      py = bod[b1].py[i];
      d1 = -99999;
      for (j=0; j<bod[b2].n; ++j) {
        d = dptl(px,py,b2,j);
        if (d > 0) break;             // not inside
        if (d > d1) { d1=d;  i1=j; }  // find nearest edge to px,py
      }
      if (j==bod[b2].n && d1<d2) {    // maximum penetration
        d2 = d1;
        i2 = i1+1;                // i2 = other end point of the edge
        if (i2==bod[b2].n) i2=0;
        *ny = (bod[b2].px[i1] - bod[b2].px[i2])/bod[b2].len[i1];
        *nx = (bod[b2].py[i2] - bod[b2].py[i1])/bod[b2].len[i1];
        *hx = px;        // + *nx*d1;
        *hy = py;        // + *ny*d1;
        *B1 = b1;
        *B2 = b2;
        *pen = -d1;
      }
    }
    i = b1; b1 = b2; b2 = i;            // swap b1 and b2
  }  
  if (d2 < 0) return 0;
  return 1;
}

/*----------------------------------------------------------------------
   get_accel  - get any changes to velocity from collisions, gravity, etc
*/
int get_accel(void)
{
  int t,t1,t2;
  int i1,b1,b2;
  static float x1,x2,y1,y2,v1,v2,c1,c2,w1,w2,m1,m2,o1,o2;
  static float nx,ny,z,a1,a2,a,b,c,hx,hy;

  // test if any body bumped into a wall

  for (b1=0; b1<numb; ++b1) {
    x1=MX; x2=0;
    for (i1=0; i1<bod[b1].n; ++i1) {
      if (x1 > bod[b1].px[i1]) {        // find leftmost
        x1=bod[b1].px[i1];
        y1=bod[b1].py[i1];
      }
      if (x2 < bod[b1].px[i1]) {        // find rightmost
        x2=bod[b1].px[i1];
        y2=bod[b1].py[i1];
      }
    }
    if (x2 > MX) { x1=x2-MX; y1=y2; }
    else if (x1 > 0) x1=0;
    if (x1) {
      y1 = bod[b1].ocy - y1;
      w1 = bod[b1].w;
      v1 = bod[b1].vx;
      a = y1*y1/bod[b1].i;
      c1 = ((a-rest)*v1 - y1*(1+rest)*w1)/(a+1);
      bod[b1].vx += c1-v1;
      bod[b1].w += y1*(c1 - v1)/bod[b1].i;
      bod[b1].cx -= scut*x1;
    }
    y1=MY; y2=0;
    for (i1=0; i1<bod[b1].n; ++i1) {
      if (y1 > bod[b1].py[i1]) {        // find bottommost
        x1=bod[b1].px[i1];
        y1=bod[b1].py[i1];
      }
      if (y2 < bod[b1].py[i1]) {        // find topmost
        x2=bod[b1].px[i1];
        y2=bod[b1].py[i1];
      }
    }
    if (y2 > MY) { y1=y2-MY; x1=x2; }
    else if (y1 > 0) y1=0;
    y2=y1;
    if (y2) {
      y1 = -bod[b1].ocx + x1;
      w1 = bod[b1].w;
      v1 = bod[b1].vy;
      a = y1*y1/bod[b1].i;
      c1 = ((a-rest)*v1 - y1*(1+rest)*w1)/(a+1);
      bod[b1].vy += c1-v1;
      bod[b1].w += y1*(c1 - v1)/bod[b1].i;
      bod[b1].cy -= scut*y2;
    }
  }
  //  For every pair of shapes, test if any of b1's corners are inside b2.

  for (t1=0; t1<numb; ++t1)
    for (t2=t1+1; t2<numb; ++t2)
      if (collision(t1,t2,&b1,&b2,&hx,&hy,&nx,&ny,&z)) {
        // Prepare to compute collision
        v1 = bod[b1].vx*nx + bod[b1].vy*ny;
        v2 = bod[b2].vx*nx + bod[b2].vy*ny;
        y1 = (hx-bod[b1].ocx)*ny - (hy-bod[b1].ocy)*nx;
        y2 = (hx-bod[b2].ocx)*ny - (hy-bod[b2].ocy)*nx;
        w1 = bod[b1].w;
        w2 = bod[b2].w;
        m1 = bod[b1].m;
        m2 = bod[b2].m;
        a1 = bod[b1].i;
        a2 = bod[b2].i;
        /*--------------------------------------------------------------------
           find change in velocity using "Dynamics" impact equations
        */
        a = (a1*m1*y2*y2 + a2*m2*y1*y1)/(a1*a2);
        b = m1 + m2;
        if (b < EPS) break;           // mass is too small
        c = y1*w1 - y2*w2;
        c1 = ((a + m1 - rest*m2)*v1 + (1+rest)*m2*(v2 - c)) / (a + b);
        c2 = ((a + m2 - rest*m1)*v2 + (1+rest)*m1*(v1 + c)) / (a + b);
        o1 = y1*(c1 - v1)/a1;
        o2 = y2*(c2 - v2)/a2;
        bod[b1].w += o1;
        bod[b2].w += o2;
        bod[b1].vx += nx*(c1 - v1);
        bod[b1].vy += ny*(c1 - v1);
        bod[b2].vx += nx*(c2 - v2);
        bod[b2].vy += ny*(c2 - v2);
        /*--------------------------------------------------------------------
           adjust the centers of the two bodies so they aren't overlapped,
           movement is inversely proportional to mass.
        */
        z *= scut/b;
        bod[b1].cx -= nx*z*m2;
        bod[b1].cy -= ny*z*m2;
        bod[b2].cx += nx*z*m1;
        bod[b2].cy += ny*z*m1;
      }

  //   Add air friction and gravity
  //   and check for centers out of range

  for (b1=0; b1<numb; ++b1) {
    z = hypot(bod[b1].vx,bod[b1].vy);
    z = 1/(1+z*fric);                // air friction
    bod[b1].vy *= z;
    bod[b1].vx *= z;
    bod[b1].vy -= grav;             // gravity

    t = 0;
    if (bod[b1].cx < 0) { t=1;  bod[b1].cx = 0; }
    if (bod[b1].cx >MX) { t=1;  bod[b1].cx = MX; }
    if (bod[b1].cy < 0) { t=1;  bod[b1].cy = 0; }
    if (bod[b1].cy >MY) { t=1;  bod[b1].cy = MY; }
    if (t) {
      bod[b1].vx *= 0.9;          // slow it down!
      bod[b1].vy *= 0.9;
      bod[b1].w *= 0.9;
    }
  }

  return 0;
}
/*----------------------------------------------------------------------
   Animate
   This is the main action loop of the program.
   On entry these must be set.
     bod, numb, xm, ym
   and graphics screen must be up and clear.
*/
int animate(void)
{
  int p=0;

  compute_pxy();           // compute corners based on new c,a
  set_poly(p);             // initialize old poly

  while (1) {
    p = 1-p;
    set_poly(p);           // set poly based on corners
    if (get_input()) break; // test mouse and keyboard
    redraw_all(p);         // erase and draw all bodies
    set_old();             // set old c,a,v,w to new ones.
    compute_new();         // new c,a = old c,a + old v,w
    compute_pxy();         // compute corners based on new c,a
    get_accel();           // test for any collisions, etc.
  }
  return 0;
}
/*----------------------------------------------------------------------
   M A I N
*/
int main(int argc, char **argv)
{
  int gm,gd=0;

  registerbgidriver(EGAVGA_driver);
  initgraph(&gd,&gm,"");
  srand(time(NULL));
  if (graphresult() != grOk) return
    printf("Unable to initialize VGA adapter");
  xm = getmaxx();
  ym = getmaxy();
  if (init_mouse()==-1) {
//    closegraph();
//    printf("Unable to initialize mouse.");
//    return 0;
//  }
//  else {
    mous = 1;
    hide_mouse();              // we make our own pointer.
  }
  numb = 0;
//  if (argc>1) sscanf(argv[1],"%d",&numb);
//  if (numb >= MB) numb=MB-1;

//  for (b=0; b<numb; ++b) random_quad(b,1.3);

/*
  numb=1;

  b=0;
  bod[b].n = 4;
  bod[b].sx[0]= t;   bod[b].sy[0]= t;
  bod[b].sx[1]= t;   bod[b].sy[1]=-t;
  bod[b].sx[2]=-t;   bod[b].sy[2]=-t;
  bod[b].sx[3]=-t;   bod[b].sy[3]= t;
  bod[b].color = random(5)+2;
  bod[b].m = 1;
  bod[b].r = 1;
  bod[b].i = 0.7;               //  I/m = r*r/2
  bod[b].cx = 1; bod[b].cy = 2;
  bod[b].a = 0;
  bod[b].vx = 0.1; bod[b].vy = 0;
  bod[b].w = 0;

  b=0;
  bod[b].n = 4;
  bod[b].sx[0]= t;   bod[b].sy[0]= t;
  bod[b].sx[1]= t;   bod[b].sy[1]=-t;
  bod[b].sx[2]=-t;   bod[b].sy[2]=-t;
  bod[b].sx[3]=-t;   bod[b].sy[3]= t;
  bod[b].color = random(5)+2;
  bod[b].m = 1;
  bod[b].r = 1;
  bod[b].i = 0.7;               //  I/m = r*r/2
  bod[b].cx = 3; bod[b].cy = 2;
  bod[b].a = 0.3;
  bod[b].vx = -0.1; bod[b].vy = 0;
  bod[b].w = 0;
*/
  animate();
  closegraph();
  return 0;
}