/* body3.c  three body gravitational attraction                        */
/*          F = G m1 m2 / d^2  F, A, V, S have x,y,z components        */
/*          F = m2 v^2 / d     d^2 = (x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2 */
/*          F = m A            acceleration along force components     */
/*          V = V + A dt       velocity update at delta time  dt       */
/*          S = S + V dt       position update at delta time  dt       */

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <GL/glut.h>

/* define number of masses NM, add initial values */
#define NM 3
static GLfloat theta[] = {0.0,0.0,0.0}; /* display angle */
GLUquadricObj *obj;
typedef struct { double x; double y; double z; } xyz;
      /* subscript is mass index */
static xyz S[NM] = {{0.0,0.0,0.0}, {1.0,0.0,0.0}, {3.0,0.0,0.0}};
static xyz V[NM]; /* initialized to orbit */
static xyz A[NM]; /* computed */
static xyz F[NM]; /* computed */ 
static GLfloat R[NM] = {0.125, 0.0625, 0.03125}; /* radius of mass */
static double M[NM]  = {1.0, 0.125, 0.00625}; /* normalized mass */
static GLfloat color[NM][3] = {{1.0,0.0,0.0}, {0.0,1.0,0.0}, {0.0,0.0,1.0}};
static double G = 2.0;  /* normalized gravitational constant */
static double dt = 0.001;
static debug = 0;


static void printstring(void *font, char *string)
{
   int len, i;

   len = (int) strlen(string);
   for (i = 0; i < len; i++)
      glutBitmapCharacter(font, string[i]);
} /* end printstring */

static void display(void)
{
  int i;
  
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glLoadIdentity();
  glRotatef(theta[0], 1.0, 0.0, 0.0);
  glRotatef(theta[1], 0.0, 1.0, 0.0);
  glRotatef(theta[2], 0.0, 0.0, 1.0);

  for(i=0; i<NM; i++)
  {
    glColor3fv(color[i]);
    glTranslatef(S[i].x, S[i].y, S[i].z);
    glutSolidSphere(R[i], 20, 20);
    glTranslatef(-S[i].x, -S[i].y, -S[i].z);
  }
  glColor3f(0.0, 0.0, 0.0);
  glRasterPos2i(100, 550);
  printstring(GLUT_BITMAP_HELVETICA_18,
  "slower: s or down, faster: f or up, reset: r");

  glutSwapBuffers(); 
} /* end display */

static void physics(void)
{
  int i, j;
  double fg, fd;
  double Dijx, Dijy, Dijz;
  double d2ij;
  
  /* Idle callback, compute force, acceleration, update velocity and position */
  for(i=0; i<NM; i++)
  {
    F[i].x = 0.0;
    F[i].y = 0.0;
    F[i].z = 0.0;
    if(debug) printf("old S%d=%f,%f,%f, V%d=%f,%f,%f\n",
                     i, S[i].x, S[i].y, S[i].z, i, V[i].x, V[i].y, V[i].z);
  }
  for(i=0; i<NM-1; i++)
  {
    for(j=i+1; j<NM; j++)
    {
      Dijx = S[i].x - S[j].x; 
      Dijy = S[i].y - S[j].y; 
      Dijz = S[i].z - S[j].z;
      d2ij = Dijx*Dijx + Dijy*Dijy + Dijz*Dijz;
      fg = G * M[i] * M[j] / d2ij;
      fd = sqrt(d2ij);
      if(debug) printf("i=%d, j=%d, fg=%f, fd=%f, d2ij=%f\n",
                       i, j, fg, fd, d2ij);
      F[i].x = F[i].x - fg*Dijx/fd;
      F[i].y = F[i].y - fg*Dijy/fd;
      F[i].z = F[i].z - fg*Dijz/fd;
      F[j].x = F[j].x + fg*Dijx/fd;
      F[j].y = F[j].y + fg*Dijy/fd;
      F[j].z = F[j].z + fg*Dijz/fd;
    }
  }
  for(i=0; i<NM; i++)
  {
    if(debug) printf("F%d = %f, %f, %f\n", i, F[i].x, F[i].y, F[i].z);
    A[i].x = F[i].x/M[i];
    A[i].y = F[i].y/M[i];
    A[i].z = F[i].z/M[i];
    if(debug)printf("A%d = %f, %f, %f\n", i, A[i].x, A[i].y, A[i].z);
    V[i].x = V[i].x + A[i].x*dt;
    V[i].y = V[i].y + A[i].y*dt;
    V[i].z = V[i].z + A[i].z*dt;
    S[i].x = S[i].x + V[i].x*dt;
    S[i].y = S[i].y + V[i].y*dt;
    S[i].z = S[i].z + V[i].z
    *dt;
  }
  glutPostRedisplay();
} /* end physics */

static void init(void)
{
  S[0].x = 0.0;
  S[0].y = 0.0;
  S[0].z = 0.0;
  S[1].x = 1.0;
  S[1].y = 0.0;
  S[1].z = 0.0;
  S[2].x = 3.0;
  S[2].y = 0.0;
  S[2].z = 0.0;
  V[0].x = 0.0;
  V[0].y = 0.0;
  V[0].z = 0.0;
  V[1].x = 0.0;
  V[1].y = 1.0;
  V[1].z = 0.0;
  V[2].x = 0.0;
  V[2].y = 0.0;
  V[2].z = 0.0;
  physics();
} /* end init */

static void mouse(int btn, int state, int x, int y)
{

} /* end mouse */

static void keyboard(unsigned char key, int x, int y)
{
  switch (key)
  {
    case 's':
      dt = dt*0.7071;
      break;
    case 'f':
      dt = dt*1.414;
      break;
    case 'r':
      init();
      break;
   }
} /* end keyboard */

static void special(int k, int x, int y)
{
  switch(k)
  {
    case GLUT_KEY_LEFT:
      break;
    case GLUT_KEY_RIGHT:
      break;
    case GLUT_KEY_DOWN:
      dt = dt*0.7071;
      break;
    case GLUT_KEY_UP:
      dt = dt*1.414;
      break;
  }
} /* end special */

static void myReshape(int w, int h)
{
  glViewport(0, 0, w, h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glOrtho(-5.0, 5.0, -5.0, 5.0, -10.0, 10.0);
  glMatrixMode(GL_MODELVIEW);
} /* end myReshape */

int main(int argc, char *argv[])
{
  /* need both double buffering and z buffer */
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
  glutInitWindowSize(600, 600);
  glutCreateWindow(argv[0]);
  glutReshapeFunc(myReshape);
  glutDisplayFunc(display);
  glutIdleFunc(physics);
  glutMouseFunc(mouse);
  glutKeyboardFunc(keyboard);
  glutSpecialFunc(special);
  glEnable(GL_DEPTH_TEST);
  obj = gluNewQuadric();
  /* gluQuadricDrawStyle(obj, GLU_LINE); */
  init();
  
  glClearColor(1.0,1.0,1.0,1.0);
  glColor3f(1.0,1.0,1.0);
  glutMainLoop();
  return 0;
} /* end main of body3.c */