#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <unistd.h>
#include <GL/glut.h>
#include <GL/gl.h>
#include <assert.h>
#include <time.h>
#include "vector.h"


double height = 0.30;  // height of the lava lamp [meters]
double width = 0.15;   // width and depth of the lamp [meters]

double radius = 0.02;  // radius of each differential fluid element [meters]
int elements = 0;      // number of fluid elements (will be chosen automatically)
double temperature_min = 10;  // temperature at the top of the lamp
double temperature_max = 30;  // temperature at the bottom of the lamp
double initial_temperature = 18; // initial temperature of fluid elements
double damping = 0;   // damping constant (kilograms per second)
double thermal_conductivity = 0.01; // thermal conductivity constant (dimensionaless)
double lava_fraction = 0.10;  // fraction of material that is "lava"
int draw_kerosene = 0; // draw in the kerosene material?
int spheresteps = 4; // how well will the spheres be drawn?

class lava_t {
public:
 double temperature;
 vector position;
 vector velocity;
 vector force;
 double radius;
 double mass;
 int type;
};

lava_t *lava = 0;


// return the temperature at a given position

double temperature(vector position) {
  // linear interpolation between temperature_min and temperature_max
  return temperature_min + ((height-position.z)/height)*(temperature_max - temperature_min);
}

double new_temperature(lava_t lava, double dt) {
 double dT = temperature(lava.position) - lava.temperature;
 return lava.temperature + dt*dT*thermal_conductivity;
}

vector buoyant_force(lava_t lava) {
 
  // assume the density of the material changes by +/- 10% at the temperature extremes
  // and is equal to unity at a mean temperature (middle of the gradient).  then th
  // force is equal to +/- 10% of the mass off the material times gravitational
  // acceleration
  if (lava.type == 1) return 0; // kerosene has no buoyant force in this model
  double temperature_mean = (temperature_max + temperature_min)/2.0;
  double delta_density = 0.10 * (lava.temperature - temperature_mean)/(temperature_max - temperature_min);
  double force = 10 * delta_density * lava.mass;
  return vector(0.0,0.0,force);
}

vector interaction_force(lava_t a, lava_t b) {
  // compute the force on a due to b
  
  vector separation = (b.position - a.position);
  double distance = separation.length();
  double radius = (a.radius + b.radius)/2.0;
  double force_magnitude = 0.0;
  double interaction_constant = 50 * (a.mass + b.mass); 

  if (distance > 2*radius)
    return 0.0;

  separation = separation.normalized();
  
  if (a.type == b.type) {
    // similar materials
    
    if (distance < 1.5*radius) {

      force_magnitude = interaction_constant * (radius - distance);
      
    } else {
      force_magnitude = interaction_constant * (distance - 2*radius);
    }
    if (force_magnitude < 0) force_magnitude *= 2;
    return force_magnitude * separation;
    
  } else {
    // dissimilar materials

    force_magnitude = interaction_constant * distance;
    return force_magnitude * separation;
  }

  return vector();
}

lava_t *create_lava(int elements) {

  lava_t *lava = (lava_t *)calloc(sizeof(lava_t), elements);

  int lava_elements = (int)(elements * lava_fraction);

  for (int i=0; i<elements; i++) {

    lava[i].position.x = width * ((double)random()/RAND_MAX) - width/2.0;
    lava[i].position.y = width * ((double)random()/RAND_MAX) - width/2.0;
    
    if (i > lava_elements) {
      lava[i].position.z = (1.0 - lava_fraction) * height * ((double)random()/RAND_MAX) + 
	lava_fraction*height;
      lava[i].type = 1; // other liquid
    } else {
      lava[i].position.z = lava_fraction * height * ((double)random()/RAND_MAX);
      lava[i].type = 0; // lava
    }

    lava[i].temperature = initial_temperature;             // celcius
    lava[i].force = 0.0;                                   // kilograms * meters / seconds
    lava[i].velocity = 0.0;                                // meters / seconds
    lava[i].radius = radius;                               // meters 
    lava[i].mass = 1000 * 1.50 * pow(lava[i].radius, 3.0); // kilograms
    
  }
  
  return lava;
}

double clamp(double x, double low, double high) {
  if (x<low) return low;
  if (x>high) return high;
  return x;
}

void iterate(lava_t *lava, int elements, double delta_t) {

  for (int i=0; i<elements; i++) {
    
    lava[i].force = buoyant_force(lava[i]);                                // buoyant force
    
    lava[i].force = lava[i].force + damping * lava[i].velocity;            // drag
    
    for (int j=0; j<elements; j++) 
      if (i!=j)
	lava[i].force = lava[i].force + interaction_force(lava[i],lava[j]);  // interaction
    
    //printf("total force on %d is %f Newtons\n",i,lava[i].force.length());
    
    // calculate the change in temperature
    
    lava[i].temperature = new_temperature(lava[i],delta_t);
    
    // let delta T pass; assume constant acceleration for the duration
    
    lava[i].position =            lava[i].position + 
      lava[i].velocity * delta_t + 
      (0.5 / lava[i].mass ) * lava[i].force * delta_t * delta_t;
    
    lava[i].velocity = lava[i].velocity + lava[i].force * (delta_t / lava[i].mass);
    lava[i].velocity = 0.0;
    
    // add in some brownian motion (perturb the position)
    
    // enforce the limits of the container
    lava[i].position.x = clamp(lava[i].position.x,-width/2.0+radius,width/2.0-radius);
    lava[i].position.y = clamp(lava[i].position.y,-width/2.0+radius,width/2.0-radius);      
    lava[i].position.z = clamp(lava[i].position.z,0.0+radius,height-radius);
    
    
  }
}



void display(void) {
 glClear(GL_COLOR_BUFFER_BIT);
 glLoadIdentity();

 glRotatef(-90,1.0,0.0,0.0);
 glTranslatef(0.0,0.0,-10*height/2.0);
 glScalef(10.0,10.0,10.0);
 //glRotatef(-20,1.0,1.0,1.0);
 // make the box

 double foo = width/2.0;

 glColor3f(40,40,90);

 for (int z=0; z<2; z++) {
   glBegin(GL_LINES);
   glVertex3f(-foo,-foo,height*z);
   glVertex3f(-foo, foo,height*z);
   glVertex3f(-foo, foo,height*z);
   glVertex3f( foo, foo,height*z);
   glVertex3f( foo, foo,height*z);
   glVertex3f( foo,-foo,height*z);
   glVertex3f( foo,-foo,height*z);
   glVertex3f(-foo,-foo,height*z);
   glEnd();
 }

 glBegin(GL_LINES);
 for (int a=0;a<2;a++)
   for (int b=0;b<2;b++) {
     glVertex3f((a?1:-1)*foo,(b?1:-1)*foo,0.0);
     glVertex3f((a?1:-1)*foo,(b?1:-1)*foo,height);
   }
 glEnd();


 // draw the lava bits
 for (int i=0; i<elements; i++) {
   switch (lava[i].type) {
   case 0: 
     glColor3f(1,.3,.3);
     break;
   case 1:
     glColor3f(0.2,0.5,0.5);     
     break;
   }

   if ((lava[i].type == 1) && (!draw_kerosene))
     continue;

   glPushMatrix();
   /* printf("drawing lava %d at (%f,%f,%f) at temperature %f (local temperature %f)\n",
	  i,
	  lava[i].position.x,
	  lava[i].position.y,
	  lava[i].position.z,
	  lava[i].temperature,
	  temperature(lava[i].position));
   */
   glTranslatef(lava[i].position.x,lava[i].position.y,lava[i].position.z);
   glutSolidSphere(lava[i].radius,spheresteps,spheresteps);   
   glPopMatrix();
 }
 glFinish();
 glutSwapBuffers();
}


void resize(int x,int y) {
 glViewport(0,0,x,y);
}

void mouse(int button,int state,int x,int y) {
  iterate(lava,elements,0.1);
  display();
}

void mouse_move(int x, int y) { 
}

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

  srandom(time(0));

  double total_volume = height*width*width;
  double lava_volume = lava_fraction*total_volume;
  double particle_volume = 1.5*M_PI*radius*radius*radius;

  printf(" total volume is %f m^3\n",total_volume);
  printf(" lava volume is %f m^3\n",lava_volume);
  printf(" particle volume is %f m^3\n",particle_volume);
  printf(" total number of elements is %d\n",(int)(total_volume/particle_volume));
  printf(" number of lava elements should be %d\n",(int)(lava_volume/particle_volume));
  
  elements = (int)(total_volume/particle_volume);

  // Initialize GLUT

  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
  glutCreateWindow("Lava Lamp");
  glutDisplayFunc(display);
  glutReshapeFunc(resize);
  glutMouseFunc(mouse);
  //  glutMotionFunc(mouse_move);

  // Initialize GL

  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glOrtho(-2.0,2.0,-2.0,2.0,-100,100);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();


  GLfloat specular[] = {1.0, 1.0, 1.0, 1.0 };
  GLfloat shininess[] = { 100.0 };
  GLfloat position[] = { 0.3, 0.3, 0.3, 0.0 };

  glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, specular);
  glMaterialfv(GL_FRONT_AND_BACK, GL_SHININESS, shininess);
  glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
  glEnable(GL_COLOR_MATERIAL);

  glLightfv(GL_LIGHT0, GL_POSITION, position);
  glEnable(GL_LIGHTING);
  glEnable(GL_LIGHT0);

  // initialize the simulatio  
  lava = create_lava(elements);

  // Switch to main loop
  glutMainLoop();
  
  free(lava);
  return 0;        
}