#include <GL/glut.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <sys/types.h>
#include <stdio.h>
#include <assert.h>
#include <string.h> /* memset */
#include <unistd.h> /* close */
#define NUM_PARTICLES    100000
#define PAUSE              0
#define NORMALIZE_SPEED    1
#define QUIT               2
GLfloat  light0Amb[4] =  { 1.0, 0.6, 0.2, 1.0 };
GLfloat  light0Dif[4] =  { 1.0, 0.6, 0.2, 1.0 };
GLfloat  light0Spec[4] = { 0.0, 0.0, 0.0, 1.0 };
GLfloat  light0Pos[4] =  { 0.0, 0.0, 0.0, 1.0 };
GLfloat  light1Amb[4] =  { 0.0, 0.0, 0.0, 1.0 };
GLfloat  light1Dif[4] =  { 1.0, 1.0, 1.0, 1.0 };
GLfloat  light1Spec[4] = { 1.0, 1.0, 1.0, 1.0 };
GLfloat  light1Pos[4] =  { 0.0, 5.0, 5.0, 0.0 };
GLfloat  materialAmb[4] = { 0.25, 0.22, 0.26, 1.0 };
GLfloat  materialDif[4] = { 0.63, 0.57, 0.60, 1.0 };
GLfloat  materialSpec[4] = { 0.99, 0.91, 0.81, 1.0 };
GLfloat  materialShininess = 27.8;

struct particleData {
  float   position[3];
  float   speed[3];
  float   color[3];
};
typedef struct particleData    particleData;

typedef struct constants {
  int nframes;
  int npframe;
  float h;
  float dt;
  float rho0;
  float k;
  float mu;
  float g;
} constants;
void default_args(constants* args) {
  args->nframes = 400;
  args->npframe = 100;
  args->dt = 1e-4;
  args->h = 5e-2;
  args->rho0 = 1000;
  args->k = 1e3;
  args->mu = 0.1;
  args->g = 9.8;
}

typedef struct flow {
  int n;
  float mass;
  float* __restrict__ rho;
  float* __restrict__ x;
  float* __restrict__ vh;
  float* __restrict__ v;
  float* __restrict__ a;
} flow;
flow* alloc_currFlow(int n) {
    flow* s = (flow*) calloc(1, sizeof(flow));
    s->n   =  n;
    s->rho =  (float*) calloc(  n, sizeof(float));
    s->x =    (float*) calloc(3*n, sizeof(float));
    s->vh =   (float*) calloc(3*n, sizeof(float));
    s->v =    (float*) calloc(3*n, sizeof(float));
    s->a =    (float*) calloc(3*n, sizeof(float));
    return s;
}
void resetFlow(flow* s) {
    free(s->a);
    free(s->v);
    free(s->vh);
    free(s->x);
    free(s->rho);
    free(s);
}

typedef int (*fillRegion)(float, float, float);

particleData particles[NUM_PARTICLES];
static constants args;
static flow* currFlow;
static int nframes;
static int npframe;
static float dt;
static int n;

int fuel = 0;
float angle = 0.1;
int  wantNormalize = 1;
int  wantPause = 0;
int box_indicator(float x, float y, float z) {
  return (x < 0.5) && (y < 0.5) && (z < 0.5) ;
}
int sph_indicator(float x, float y, float z) {
  float dx = (x-0.5);
  float dy = (y-0.5);
  float dz = (y-0.3);
  float r3 = dx*dx + dy*dy + dz*dz;
  return (r3 < 0.25*0.25*0.25);
}

flow* place_particles(constants* param, fillRegion fillNum) {
  float h = param->h;
  float cellSize = h/1.3;

  int count = 0;
  float x;
  float y;
  float z;
  for (x = 0; x < 1; x += cellSize)
    for (y = 0; y < 1; y += cellSize)
	  for (z = 0; z < 1; z += cellSize)
       count += fillNum(x,y,z);

  flow* s = alloc_currFlow(count);
  int p = 0;
    for (x = 0; x < 1; x += cellSize) {
      for (y = 0; y < 1; y += cellSize) {
	    for (z = 0; z < 1; z += cellSize) {
          if (fillNum(x,y,z)) {
            s->x[2*p+0] = x;
            s->x[2*p+1] = y;
            s->v[2*p+0] = 0;
            s->v[2*p+1] = 0;
            ++p;
	      }
        }
      }
    }
  return s;
}

void compute_density(flow* s, constants* args) {
  int i = 0;
  int j = 0;
  int n = s->n;
  float* __restrict__ rho = s->rho;
  const float* __restrict__ x = s->x;

  float h = args->h;
  float h3 = h*h*h;
  float h12 = ( h3*h3 )*( h3*h3 );
  float C = 4 * s->mass / M_PI / h12;

  memset(rho, 0, n*sizeof(float));
  for (i = 0; i < n; ++i) {
    rho[i] += 4 * s->mass / M_PI / h3;
    for (j = i+1; j < n; ++j) {
      float dx = x[2*i+0]-x[2*j+0];
      float dy = x[2*i+1]-x[2*j+1];
      float dz = x[2*i+2]-x[2*j+2];
      float r3 = dx*dx + dy*dy + dz*dz;
      float sub = h3-r3;
      if (sub > 0) {
        float rho_ij = C*sub*sub*sub;
        rho[i] += rho_ij;
        rho[j] += rho_ij;
      }
    }
  }
}

void compute_accel(flow* currFlow, constants* args) {
  int i;
  int j;

  const float h = args->h;
  const float rho0 = args->rho0;
  const float k = args->k;
  const float mu = args->mu;
  const float g = args->g;
  const float mass = currFlow->mass;
  const float h3 = h*h*h;
  
  const float* __restrict__ rho = currFlow->rho;
  const float* __restrict__ x = currFlow->x;
  const float* __restrict__ v = currFlow->v;
  float* __restrict__ a = currFlow->a;
  int n = currFlow->n;
  
  compute_density(currFlow, args);
  
  for (i = 0; i < n; ++i) {
    a[2*i+0] = 0;
    a[2*i+1] = 0;
    a[2*i+2] = -g; //gravity
  }
  
  float C0 = mass / M_PI / ( (h3)*(h3) );
  float Cp = 15*k;
  float Cv = -40*mu;
  
  for (i = 0; i < n; ++i) {
    const float rhoi = rho[i];
     for (j = i+1; j < n; ++j) {
      float dx = x[2*i+0]-x[2*j+0];
      float dy = x[2*i+1]-x[2*j+1];
      float dz = x[2*i+2]-x[2*j+2];
      float r3 = dx*dx + dy*dy + dz*dz;
      if (r3 < h3) {
        const float rhoj = rho[j];
        float q = cbrt(r3)/h;
        float u = 1-q;
        float w0 = C0 * u/rhoi/rhoj;
        float wp = w0 * Cp * (rhoi+rhoj-2*rho0) * u/q;
        float wv = w0 * Cv;
        float dvx = v[2*i+0]-v[2*j+0];
        float dvy = v[2*i+1]-v[2*j+1];
        float dvz = v[2*i+2]-v[2*j+2];
        a[2*i+0] += (wp*dx + wv*dvx);
        a[2*i+1] += (wp*dy + wv*dvy);
        a[2*i+2] += (wp*dz + wv*dvz);
        a[2*j+0] -= (wp*dx + wv*dvx);
        a[2*j+1] -= (wp*dy + wv*dvy);
        a[2*i+2] -= (wp*dz + wv*dvz);
      }
    }
  }
}

void damp_reflect(int which, float barrier, float* x, float* v, float* vh) {
  const float DAMP = 0.75;
  if (v[which] == 0)
    return;
  float tbounce = (x[which]-barrier)/v[which];
  x[0] -= v[0]*(1-DAMP)*tbounce;
  x[1] -= v[1]*(1-DAMP)*tbounce;
  x[2] -= v[2]*(1-DAMP)*tbounce;
  x[which] = 2*barrier-x[which];
  v[which] = -v[which];
  vh[which] = -vh[which];
  v[0] *= DAMP; vh[0] *= DAMP;
  v[1] *= DAMP; vh[1] *= DAMP;
  v[2] *= DAMP; vh[2] *= DAMP;
}

void reflect_bc(flow* s) {
  int i;
  const float XMIN = 0.0;
  const float XMAX = 1.0;
  const float YMIN = 0.0;
  const float YMAX = 1.0;
  const float ZMIN = 0.0;
  const float ZMAX = 1.0;
  float* __restrict__ vh = s->vh;
  float* __restrict__ v = s->v;
  float* __restrict__ x = s->x;
  int n = s->n;
  for (i = 0; i < n; ++i, x += 2, v += 2, vh += 2) {
    if (x[0] < XMIN) damp_reflect(0, XMIN, x, v, vh);
    if (x[0] > XMAX) damp_reflect(0, XMAX, x, v, vh);
    if (x[1] < YMIN) damp_reflect(1, YMIN, x, v, vh);
    if (x[1] > YMAX) damp_reflect(1, YMAX, x, v, vh);
    if (x[2] < ZMIN) damp_reflect(2, ZMIN, x, v, vh);
    if (x[2] > ZMAX) damp_reflect(2, ZMAX, x, v, vh);
  }
}

void leapfrog_step(flow* s, double dt) {
  int i;
  const float* __restrict__ a = s->a;
  float* __restrict__ vh = s->vh;
  float* __restrict__ v = s->v;
  float* __restrict__ x = s->x;
  int n = s->n;
  for (i = 0; i < 2*n; ++i) vh[i] += a[i] * dt;
  for (i = 0; i < 2*n; ++i) v[i] = vh[i] + a[i] * dt / 2;
  for (i = 0; i < 2*n; ++i) x[i] += vh[i] * dt;
  reflect_bc(s);
}

void leapfrog_start(flow* s, double dt) {
  int i;
  const float* __restrict__ a = s->a;
  float* __restrict__ vh = s->vh;
  float* __restrict__ v = s->v;
  float* __restrict__ x = s->x;
  int n = s->n;
  for (i = 0; i < 2*n; ++i) {
    vh[i] = v[i] + a[i] * dt / 2;}
  for (i = 0; i < 2*n; ++i){
    v[i] += a[i] * dt;}
  for (i = 0; i < 2*n; ++i){
    x[i] += vh[i] * dt;}
  reflect_bc(s);
}

void normalize_mass(flow* s, constants* param) {
  int i;
  s->mass = 1;
  compute_density(s, param);
  float rho0 = param->rho0;
  float rho2s = 0;
  float rhos = 0;
  for (i = 0; i < s->n; ++i) {
    rho2s += (s->rho[i])*(s->rho[i]);
    rhos += s->rho[i];
  }
  s->mass *= ( rho0*rhos / rho2s );
}

flow* init_particles(constants* param) {
  flow* s = place_particles(param, box_indicator);
  normalize_mass(s, param);
  return s;
}

void newSpeed (float dest[3]) {
  float    x;
  float    y;
  float    z;
  float    len;

  x = (2.0 * ((GLfloat) rand ()) / ((GLfloat) RAND_MAX)) - 1.0;
  y = (2.0 * ((GLfloat) rand ()) / ((GLfloat) RAND_MAX)) - 1.0;
  z = (2.0 * ((GLfloat) rand ()) / ((GLfloat) RAND_MAX)) - 1.0;

  if (wantNormalize)
    {
      len = sqrt (x * x + y * y + z * z);

      if (len)
	{
	  x = x / len;
	  y = y / len;
	  z = z / len;
	}
    }

  x = x * 0.2;
  y = y * 0.2;
  z = z * 0.2;

  dest[0] = x;
  dest[1] = y;
  dest[2] = z;
}

void newExplosion (void) {
  int    i;

  for (i = 0; i < NUM_PARTICLES; i++)
    {
      particles[i].position[0] = 0.0;
      particles[i].position[1] = 0.0;
      particles[i].position[2] = 0.0;

      particles[i].color[0] = 1.0;
      particles[i].color[1] = 1.0;
      particles[i].color[2] = 0.5;

      newSpeed (particles[i].speed);
    }
  fuel = 100;
}

void display (void) {
  int    i;
  glClear (GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glLoadIdentity ();
  glTranslatef (0.0, 0.0, -10.0);
  glRotatef (angle, 0.0, 1.0, 0.0);

  if (fuel > 0)
    {
      glPushMatrix ();

      glDisable (GL_LIGHTING);
      glDisable (GL_DEPTH_TEST);

      glBegin (GL_POINTS);

      for (i = 0; i < NUM_PARTICLES; i++)
	{
	  glColor3fv (particles[i].color);
	  glVertex3fv (particles[i].position);
	}

      glEnd ();

      glPopMatrix ();

      glEnable (GL_LIGHTING);
      glEnable (GL_LIGHT0);
      glEnable (GL_DEPTH_TEST);

      glNormal3f (0.0, 0.0, 1.0);
    }

  glutSwapBuffers ();
}

void keyboard (unsigned char key, int x, int y) {
  switch (key)
    {
    case ' ':
    newExplosion ();
    printf("Space Pressed!\n");
      break;

    case 27:
      exit (0);
      break;

    case 'p':
      wantPause = 1 - wantPause;
      break;
    }
}

void check_currFlow(flow* s) {
  int i;
  for (i = 0; i < s->n; ++i) {
    float xi = s->x[2*i+0];
    float yi = s->x[2*i+1];
    float zi = s->x[2*i+2];
  }
}

void prepareSim (void) {
  default_args(&args);
  currFlow = init_particles(&args);
  nframes = args.nframes;
  npframe = args.npframe;
  dt = args.dt;
  n = currFlow->n;
  compute_accel(currFlow, &args);
  leapfrog_start(currFlow, dt);
  check_currFlow(currFlow);
}

void newPosition (void) {
  int i;
  if (!wantPause) {
    if (fuel > 0) {
  int i;
  int frame;
  printf("newPosition Running!!\n");
  /*
  for (frame = 1; frame < nframes; ++frame) {
    for (i = 0; i < npframe; ++i) {
      compute_accel(currFlow, &args);
      leapfrog_step(currFlow, dt);
      check_currFlow(currFlow);
    }
  }
*/
    
		
		
      if (fuel > 0) {
	  for (i = 0; i < NUM_PARTICLES; i++)
	    {
	      particles[i].position[0] += particles[i].speed[0];
	      particles[i].position[1] += particles[i].speed[1];
	      particles[i].position[2] += particles[i].speed[2];

	      particles[i].color[0] -= 1.0 / 500.0;
	      if (particles[i].color[0] < 0.0)
		{
		  particles[i].color[0] = 0.0;
		}

	      particles[i].color[1] -= 1.0 / 100.0;
	      if (particles[i].color[1] < 0.0)
		{
		  particles[i].color[1] = 0.0;
		}

	      particles[i].color[2] -= 1.0 / 50.0;
	      if (particles[i].color[2] < 0.0)
		{
		  particles[i].color[2] = 0.0;
		}
	    }
}
	  --fuel;
	}
 angle += 0.4;
    }

  glutPostRedisplay ();
}

void reshape (int w, int h) {
  glViewport (0.0, 0.0, (GLfloat) w, (GLfloat) h);

  glMatrixMode (GL_PROJECTION);
  glLoadIdentity ();
  gluPerspective (45.0, (GLfloat) w / (GLfloat) h, 0.1, 100.0);

  glMatrixMode (GL_MODELVIEW);
}

void menuSelect (int value) {
  switch (value)
    {
    case PAUSE:
      wantPause = 1 - wantPause;
      break;

    case NORMALIZE_SPEED:
      wantNormalize = 1 - wantNormalize;
      break;

    case QUIT:
      exit (0);
      break;
    }
}

int main (int argc, char   *argv[]) {
  prepareSim();

  glutInit (&argc, argv);
  glutInitDisplayMode (GLUT_DOUBLE | GLUT_DEPTH | GLUT_RGB);
  glutCreateWindow ("Explosion");
  glutKeyboardFunc (keyboard);

  glutIdleFunc (newPosition); 
  glutDisplayFunc (display);
  glutReshapeFunc (reshape);
  srand (time (NULL));
  glEnable (GL_LIGHT0);
  glEnable (GL_LIGHT1);
  glLightfv (GL_LIGHT0, GL_AMBIENT, light0Amb);
  glLightfv (GL_LIGHT0, GL_DIFFUSE, light0Dif);
  glLightfv (GL_LIGHT0, GL_SPECULAR, light0Spec);
  glLightfv (GL_LIGHT0, GL_POSITION, light0Pos);
  glLightfv (GL_LIGHT1, GL_AMBIENT, light1Amb);
  glLightfv (GL_LIGHT1, GL_DIFFUSE, light1Dif);
  glLightfv (GL_LIGHT1, GL_SPECULAR, light1Spec);
  glLightfv (GL_LIGHT1, GL_POSITION, light1Pos);
  glLightModelf (GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
  glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT, materialAmb);
  glMaterialfv (GL_FRONT_AND_BACK, GL_DIFFUSE, materialDif);
  glMaterialfv (GL_FRONT_AND_BACK, GL_SPECULAR, materialSpec);
  glMaterialf (GL_FRONT_AND_BACK, GL_SHININESS, materialShininess);
  glEnable (GL_NORMALIZE);
  glutCreateMenu (menuSelect);
  glutAddMenuEntry ("Pause", PAUSE);
  glutAddMenuEntry ("Toggle normalized speed vectors", NORMALIZE_SPEED);
  glutAddMenuEntry ("Quit", QUIT);
  glutAttachMenu (GLUT_RIGHT_BUTTON);

  glutMainLoop ();
  resetFlow(currFlow);
  return 0;
}