#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 <iostream.h>

#include "vector.h"
#include "lava.h"
#include "grid.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.01;  // 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

/* How much should velocity be damped?  This constant determines the damping force. */

double damping = -0.03;   // damping constant (kilograms per second)

/* Thermal conductivity constant... not particularly well defined. */

double thermal_conductivity = 0.9; // thermal conductivity constant (dimensionaless)
double glass_conductivity = 0.99;

/* Specify the fraction of the volume that will be filled with lava.  The remainder will be
   filled with kerosene */

double lava_fraction = 0.50; 

/* There are two types of particles in the system: wax ("lava") and kerosene ("the other liquid").
   Sometimes it's instructive to draw the kerosene, although normally it's invisible/transparent.  */

bool draw_kerosene = false; 

/* We can specify how many samples are used when drawing a sphere.  */

int spheresteps = 4; 

/* Do we want OpenGL lighting, or just a flat display?  Lighting would be better, but it doesn't
   work very well right now.  Needs to be fixed. */

bool use_lighting = false;

/* The naive algorithm for particle interactions is just to compare every particle to every other
   particle.  Obviously this is O(n^2) with respect to the number of particles and rapidly becomes
   very slow.  Because only particles nearby in space interact, we don't need to do that if we have
   some idea which particles may be nearby.  If use_grid is true, then a 3D grid structure will be
   used to find the nearest neighbors, dramatically increasing the speed of the simulation.  There's
   really no reason to not use it. */

bool use_grid = true;

/* If thermal diffusion is on, then heat is transferred between nearby particles.  If thermal
   diffusion is off, then we just assume a static temperature gradient from the bottom to the
   top of the lamp.
*/

bool thermal_diffusion = true; 

/* If thermal diffusion is on, and a particle is within
   boundary_effect distance from a boundary, then its temperature will
   move towards the temperature of the boundary.
*/
double boundary_effect = radius * 2;

int verbosity = 0; // verbosity


double delta_t = 0.1;   // this is the time step used for Eulerian integration

bool color_code = true;
bool use_force_interaction = true;

bool volume_is_rectangular = true;

lava_t *lava = 0;        // this is the array of particles
grid_t *the_grid;        // this is a grid used to find nearest neighbors



// return a random number in the range [0,1]
// note: some operating systems (eg, Solaris) lie about RAND_MAX

double my_random(void) {
  return (double)rand()/RAND_MAX;
}

// for a cylindrical simulation volume, this function returns the radius as a function of height

double container_radius(double z) {
  return width/2.0;
}


// return the temperature at a given position
// not used if thermal diffusion is turned on

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

double boundary_temperature(vector position) {
  if (position.z < (height/10)) 
    return temperature_max;
  else
    return temperature_min;
}

// return the new temperature of a particle of lava after a given time interval
// not used if thermal diffusion is turned on

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


// This is the thermal diffusion interaction, used to calculate the change in
// temperature of particle "a" due to particle "b" in the time interval dt.
// distance is the separation of a and b, eg (*b - *a).length()

void thermal_interaction(lava_t *a, const lava_t *b, double dt, double distance) {
  
  double dT = b->temperature - a->temperature;
  double radius = (a->radius + b->radius)/2.0;

  // there's a problem, where multiplying by radius/distance will result in NaN. I think
  // this is because particles are getting too close together.  I'm not sure how to
  // prevent this.  I don't want the particles to be closer together than one radius.

  if (distance < 3*radius) 
    a->dT = a->dT + dt*dT*thermal_conductivity; // *radius/distance;
  
}

// calculate the buoyant force on a bit of lava.  this has not been adapted for use with
// thermal diffusion.

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;
  temperature_mean = temperature(lava.position);
  double delta_density = (lava.temperature - temperature_mean)/(temperature_max - temperature_min);
  double force = .5 * delta_density * lava.mass;
  temperature_mean = (temperature_max + temperature_min)/2.0;
  if (force < 0) {
      force *= 1.5;
  } else {
      force *= 1.1;
  }
  /*
  if (lava.temperature < temperature_mean) {
      force += .7 * delta_density * lava.mass;
  } else {
      force += .7 * delta_density * lava.mass;
  }
  */
  return vector(0.0,0.0,force);
}

int interactions = 0;
int successful_interactions;

// calculate the force on "a" due to "b", given their separation vector and the
// length of that vector.

vector interaction_force(lava_t a, lava_t b, vector separation, double distance) {
  // compute the force on a due to b

  interactions++;

  // 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 =  (a.mass + b.mass); 
  double abs_temp = (a.temperature < b.temperature) ?
      b.temperature - a.temperature :
      a.temperature - b.temperature;

  if (distance > 2*radius)
    return 0.0;

  successful_interactions++;

  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;
    if (abs_temp < (temperature_max - temperature_min) * 0.1)
	force_magnitude = -force_magnitude;

    return force_magnitude * separation;
    
  } else {
    // dissimilar materials

         force_magnitude = interaction_constant * distance;
    if (abs_temp < (temperature_max - temperature_min) * 0.1)
	force_magnitude = -force_magnitude;
    return force_magnitude * separation;
  }

  return vector();
}


// create all the wax and kerosene particles

lava_t *create_lava(int elements) {

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

  if (verbosity > 2) printf("creating lava elements (%d)\n", elements);

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

  for (int i=0; i<elements; i++) {
    if (i > lava_elements) {
      lava[i].position.z = (1.0 - lava_fraction) * height * my_random() + 
	lava_fraction*height;
      lava[i].type = 1; // other liquid
    } else {
      lava[i].position.z = lava_fraction * height * my_random();
      lava[i].type = 0; // lava
    }

    do {
      lava[i].position.x = width * my_random() - width/2.0;
      lava[i].position.y = width * my_random() - width/2.0;
    } while ((!volume_is_rectangular) && 
	     (sqrt(pow(lava[i].position.x,2) + pow(lava[i].position.y,2)) + radius >
	      container_radius(lava[i].position.z)));
    
    if (verbosity > 3) printf(" creating lava element %d type %d location (%f,%f,%f)\n",
			      i,
			      lava[i].type,
			      lava[i].position.x,
			      lava[i].position.y,
			      lava[i].position.z);

    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
    
  }

  // add the lava to the grid

  for (int i=0; i<elements; i++) {
//    printf("adding element %d to grid\n",i);
    the_grid->add(&lava[i]);
  }

  return lava;
}

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

// this is the main time step function

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

  // we have to create a new grid, into which we'll store the lava at its new positions
  // at the end, we'll replace the old grid with the new one, and delete the old one

  grid_t *new_grid = new grid_t(the_grid->grid_min_x, the_grid->grid_min_y, the_grid->grid_min_z,
			      the_grid->grid_max_x, the_grid->grid_max_y, the_grid->grid_max_z,
			      the_grid->grid_nx,    the_grid->grid_ny,    the_grid->grid_nz);

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

    // start off with no forces, and no temperature change

    lava[i].dT = 0.0;                  // temperature change
    lava[i].force = 0.0;               // net force

    
    lava[i].force = lava[i].force + buoyant_force(lava[i]);                // buoyant force
    lava[i].force = lava[i].force + damping * lava[i].velocity;            // drag

    if (use_grid) {
      int x,y,z;
      the_grid->whichcell(lava[i].position, &x, &y, &z);
     
      for (int ni=-1; ni<2; ni++) 
	for (int nj=-1; nj<2; nj++) 
	  for (int nk=-1; nk<2; nk++) {
	    
	    list<lava_t *> *cell_list = the_grid->cellptr(x+ni,y+nj,z+nk);
	    if (!cell_list) continue;
	  
	    list<lava_t *>::const_iterator pos;
	    
	    for (pos = cell_list->begin(); pos != cell_list->end(); ++pos) {

	      lava_t *a = &lava[i];
	      lava_t *b = *pos;

	      if (a==b) continue; // no self-interaction
	      if (a->type != b->type) continue;

	      // This loops through all neighborly interactions

	      // precompute the separation just to avoid duplication of effort

	      vector separation = b->position - a->position;
	      double distance = separation.length();

	      // force interaction
	      if (use_force_interaction)
		  lava[i].force = lava[i].force + interaction_force(lava[i],**pos,separation,distance);
	      
	      // thermal interaction
	      
	      if (thermal_diffusion) {
		thermal_interaction(a,b,delta_t,distance);
	      }
	    }
	  }
    } else {
      
      // brute force O(n^2) mechanism, if you don't want to use the grid method

      for (int j=0; j<elements; j++) 
	if (i!=j) {
	  vector separation = lava[j].position - lava[i].position;
	  double distance = separation.length();

	  lava[i].force = lava[i].force + interaction_force(lava[i],lava[j],separation,distance);

	  if (thermal_diffusion)
	    thermal_interaction(&lava[i],&lava[j],delta_t,distance);
	}
    }

    // calculate interaction with boundary.  if a particle is within a certain distance of
    // the boundary, we'll invent a virtual particle on the boundary and have it interact
    // with this one

    if (thermal_diffusion) {
	if (volume_is_rectangular) {
	    if (lava[i].position.x <= -width/2.0+boundary_effect ||
		lava[i].position.x >= width/2.0-boundary_effect)
	    {
		double dT = temperature(lava[i].position) -
		    lava[i].temperature;
		lava[i].dT -= dT * delta_t * glass_conductivity;
	    }
	    if (lava[i].position.y <= -width/2.0+boundary_effect ||
		lava[i].position.y >= width/2.0-boundary_effect)
	    {
		double dT = temperature(lava[i].position) -
		    lava[i].temperature;
		lava[i].dT -= dT * delta_t * glass_conductivity;
	    }
	} else {
	    double rho = sqrt(pow(lava[i].position.x,2) + pow(lava[i].position.y,2));
	    if (rho >= container_radius(lava[i].position.z) - boundary_effect) {
		double dT = temperature(lava[i].position) - lava[i].temperature;
		lava[i].dT -= dT * delta_t * glass_conductivity;
	    }
	    
	}
	if (lava[i].position.z <= 0.0+boundary_effect) {
	    double dT = temperature_max - lava[i].temperature;
	    lava[i].dT += dT * delta_t * thermal_conductivity;
	} else if (lava[i].position.z >= height-boundary_effect) {
	    double dT = temperature_min - lava[i].temperature;
	    lava[i].dT += dT * delta_t * thermal_conductivity;
	}
    }

    // for now we'll just do some simple thermal stuff

    /* THIS IS NOT HERE YET */

    // printf("total force on %d is %f Newtons\n",i,lava[i].force.length());
    
    // calculate the change in temperature    

    if (thermal_diffusion) {
	lava[i].temperature = lava[i].temperature + lava[i].dT;
	lava[i].temperature = clamp(lava[i].temperature, temperature_min, temperature_max);
	if (lava[i].type == 0) {
	    lava[i].radius = radius * (1 + 0.25 * (lava[i].temperature -
						   initial_temperature) /
				       (temperature_max -
					temperature_min));
	} else {
	    lava[i].radius = radius * (1 + 0.05 * (lava[i].temperature -
						   initial_temperature) /
				       (temperature_max -
					temperature_min));
	}
    } else {
	// assume a simple gradient
	lava[i].temperature = new_temperature(lava[i],delta_t);
    }

    // let delta T pass; assume constant acceleration for the duration

    // EULERIAN INTEGRATION

    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;  //we might want complete damping
    
    // we might want to add in some brownian motion (perturb the position)
    
    // enforce the limits of the container

    if (volume_is_rectangular) {
	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);
    } else {
	double theta = atan2(lava[i].position.y,lava[i].position.x);
	double rho = sqrt(pow(lava[i].position.x,2) + pow(lava[i].position.y,2));

	lava[i].position.z = clamp(lava[i].position.z,0.0+radius,height-radius);
	
	if (rho > container_radius(lava[i].position.z)-radius) { 
	    rho = container_radius(lava[i].position.z)-radius;
	    lava[i].position.x = rho * cos(theta);
	    lava[i].position.y = rho * sin(theta);
	}
    }

    // add this bit of lava to the new grid, at its new location

    if (use_grid) {
      new_grid->add(&lava[i]);     
    }

  }

  // replace the old grid with the new grid

  if (use_grid) {
    delete the_grid;
    the_grid = new_grid;
  }

}

void color_scale(float x, float *zr, float *zg, float *zb){
  float r,g,b;
  x *= 255.0;
  if (x > 128.0) {
    if (x > 192.0) {
      r = 255.0;
      g = (255.0 - x)*4.0;
      b = 0.0;
    } else  {
      r = (x - 128.0)*4.0;
      g = 255.0;
      b = 0.0;
    }
  } else {
    if (x > 64.0) {
      r = 0;
      g = 255.0;
      b = 255 - (x-64.0)*4;
    } else {
      r = 0;
      g = x*4;
      b = 255.0;
    }
  }
 
  *zr = r/255.0;
  *zg = g/255.0;
  *zb = b/255.0;
 
}
void display(void) {

  glPushMatrix();
  glClear(GL_COLOR_BUFFER_BIT);
  // glLoadIdentity();

  /*
  glScalef(view_sdepth,view_sdepth,view_sdepth); // spitzer
  glRotatef(-view_stheta,1.,0.,0.);              // spitzer
  glRotatef(view_sphi,0.,0.,1.);                 // spitzer

  */

  glRotatef(-90,1.0,0.0,0.0);
  glTranslatef(0.0,0.0,-10*height/2.0);
  glScalef(10.0,10.0,10.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: 
	if(color_code) {
	    float r,g,b;
	    float relative_temperature;
	    relative_temperature = (lava[i].temperature - temperature_min)/(temperature_max - temperature_min);
	    color_scale(relative_temperature,&r,&g,&b);
	    glColor3f(r,g,b);
	} else {
	    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();
  glPopMatrix();
  glutSwapBuffers();
}


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

void mouse_move(int x, int y) {       
  // TODO: Add crystal ball interface
}

void mouse(int button, int state, int x, int y) { 
  // TODO: Add crystal ball interface
}


void key(unsigned char k, int x, int y) {
  return;
}


void idle(void) {

  //printf("calculating...\n");
  interactions = 0;
  successful_interactions = 0;
  iterate(lava,elements,delta_t);
  //printf("%d interactions, %d successful (%d elements)\n",interactions,successful_interactions,elements);
  //printf("drawing...\n");
  glutPostRedisplay();
  return;
}

void rotate(float degree, float x, float y, float z) {
  cout << "Rotate " << degree << " " << x << " " << y << " " << z << endl;
}

void move(float x, float y, float z) {
  cout << "Translate " << x << " " << y << " " << z << endl;
}

/* for initial setting of world coordinates */
void setCamera(float x, float y, float z, float degree, float xx, float yy, float zz) {
  move(-x, -y, -z);
  rotate(degree, xx, yy, zz);
}

void light(char *type, float intensity) {
  cout << "LightSource \"" << type << "\" 3 \"intensity\" " << intensity << endl;
}

void color(float r, float g, float b) {
  cout << "Color [" << r << " " << g << " " << b << "]" << endl;
}

/* specify shader */
void surface(char *shader, char *paramters) {
  cout << "Surface \"" << shader << "\" " << paramters << endl;
}

void polygon(float x00, float x01, float x02, float x10, float x11,
	     float x12, float x20, float x21, float x22, float x30,
	     float x31, float x32)
{
  cout << "Polygon \"P\" [ "
    << x00 << " "
    << x01 << " "
    << x02 << " "
    << x10 << " "
    << x11 << " "
    << x12 << " "
    << x20 << " "
    << x21 << " "
    << x22 << " "
    << x30 << " "
    << x31 << " "
    << x32 << " "
    << "]" << endl;
}

void sphere(float radius) {
  cout << "Sphere " << radius << " -" << radius << " " << radius << " 360" << endl;
}

/* commons */
void  beginWorld()      { cout << "WorldBegin" << endl; }
void  endWorld()        { cout << "WorldEnd" << endl; }
void  beginTransform()  { cout << "TransformBegin" << endl; }
void  endTransform()    { cout << "TransformEnd" << endl; }
void  beginAttribute()  { cout << "AttributeBegin" << endl; }
void  endAttribute()    { cout << "AttributeEnd" << endl; }

double lampSize = 0.22 * 2.95;
double largeWidth = (0.1 + 7.0 / 80.0) * 2.95;
double smallWidth = 0.075 * 2.95;
double totalHeight = .725 * 2.75;
double topHeight = 0.4 * 2.7;
double upperHeight = 0.1 * 2.7;

void lamp_strip(float width, float length, char *shader) {
  surface(shader,"");
  polygon(0,0,0,
          width,0,0,
          width,length,0,
          0,length,0);
}

void lamp_cone(float degreeIncrement) {
  float degree;
  float radius, width, length;

  float l1,l2,angle;

  radius = largeWidth / 2.0;
  width = 2 * 3.141592654 * radius * degreeIncrement / 360.0;

  l1 = (largeWidth - smallWidth) / 2.0;
  l2 = (totalHeight - topHeight) / 2.0;
  length = sqrt(l1*l1 + l2*l2);

  angle = atan2(l1, l2) * 360.0 / 2.0 / 3.141592654;

  for (degree = 0; degree < 360; degree += degreeIncrement)
  {
    beginTransform();
      rotate(degree, 0,1,0);
      move(0,0,lampSize * radius);
      rotate(angle,-1,0,0);
      lamp_strip(lampSize * width, lampSize * length ,"plastic");
    endTransform();
  }  
}

void lamp_stand(float degreeIncrement) {
  lamp_cone(degreeIncrement);
}

void lamp_bottom(float degreeIncrement) {
  beginTransform();
    move(0,(totalHeight - topHeight) * lampSize,0);
    rotate(180,1,0,0);
    lamp_cone(degreeIncrement);
  endTransform();
}

void lamp_glass(float degreeIncrement) {
  float degree;
  float radius, width, length;

  float l1,l2,angle;

  radius = largeWidth / 2.0;
  width = 2 * 3.141592654 * radius * degreeIncrement / 360.0;

  l1 = ((topHeight - upperHeight) / topHeight) * (largeWidth - smallWidth) / 2.0;
  l2 = (topHeight - upperHeight);
  length = sqrt(l1*l1 + l2*l2);

  angle = atan2(l1, l2) * 360.0 / 2.0 / 3.141592654;

  beginTransform();
    move(0,(totalHeight - topHeight) * lampSize,0);

    for (degree = 0; degree < 360; degree += degreeIncrement)
    {
      beginTransform();
        rotate(degree, 0,1,0);
        move(0,0,lampSize * radius);
        rotate(angle,-1,0,0);
        lamp_strip(lampSize * width, lampSize * length, "myglass");
      endTransform();
    }  
  endTransform();
}

void lamp_top(float degreeIncrement) {
  float degree;
  float radius, width, length;

  float l1,l2,angle;

  radius = ((upperHeight / topHeight) * (largeWidth - smallWidth) + smallWidth) / 2.0;
  width = 2 * 3.141592654 * radius * degreeIncrement / 360.0;

  l1 = (upperHeight / topHeight) * (largeWidth - smallWidth) / 2.0;
  l2 = upperHeight;
  length = sqrt(l1*l1 + l2*l2);

  angle = atan2(l1, l2) * 360.0 / 2.0 / 3.141592654;

  beginTransform();
    move(0,(totalHeight - upperHeight) * lampSize,0);

    for (degree = 0; degree < 360; degree += degreeIncrement)
    {
      beginTransform();
        rotate(degree, 0,1,0);
        move(0,0,lampSize * radius);
        rotate(angle,-1,0,0);
        lamp_strip(lampSize * width, lampSize * length, "plastic");
      endTransform();
    }  
  endTransform();
}

void lamp() {
    float r = 0;
    float g = 0.5;
    float b = 0.75;

    float step = 1;
    float level = 0.85;
    
    /* lightsource within lamp */
    beginTransform();
    move(0,(totalHeight - topHeight) * lampSize / 2.0,0);
    light("pointlight", 140.0);
    endTransform();
    
    color(r,g,b);
    lamp_stand(step);
    lamp_bottom(step);
    lamp_top(step);
    
    color(level * r,level * g,level * b);
    lamp_glass(step);
}

void output_frame(lava_t *lava, int elements, int frameno) {
    // The Frame Stuff
    cout << "FrameBegin " << frameno << endl;
    cout << "Display \"frame" << frameno << ".tif\" \"tiff\" \"rgb\"" << endl;
    cout << "Format 200 400 1.0" << endl;
    cout << "Projection \"perspective\" \"fov\" 90" << endl;
    cout << "Translate 0 -0.3 0.5" << endl;

    // The Lights
    cout << "WorldBegin" << endl;
    cout << "LightSource \"ambientlight\" 1 \"intensity\" 0.5" << endl;
    cout << "TransformBegin" << endl;
    cout << "Translate 3 4 -8" << endl;
    cout << "LightSource \"pointlight\" 2 \"intensity\" 120" << endl;
    cout << "TransformEnd" << endl;

    // The Table
    cout << "AttributeBegin" << endl;
    cout << "Surface \"wood2\"" << endl;
    cout << "TransformBegin" << endl;
    cout << "Translate 0 -0.5 0" << endl;
    cout << "Polygon \"P\" [-2 0 -2 -2 0 2 2 0 2 2 0 -2]" << endl;
    cout << "TransformEnd" << endl;
    cout << "AttributeEnd" << endl;

    // The Lamp
    cout << "TransformBegin" << endl;
    cout << "Translate 0 -0.475 0" << endl;
    lamp();
    cout << "TransformEnd" << endl;

    // The Wax
    cout << "AttributeBegin" << endl;
    cout << "Color [1 1 0]" << endl;
    cout << "Surface \"myplastic\"" << endl;
    for(int i = 0; i < elements; i++) {
	cout << "TransformBegin" << endl;
	printf("Translate %f %f %f\n", lava[i].position.x, lava[i].position.z, lava[i].position.y);
	printf("Sphere %f %f %f 360\n", lava[i].radius, -lava[i].radius, lava[i].radius);
	cout << "TransformEnd" << endl;
    }

    // The End
    cout << "AttributeEnd" << endl;
    cout << "WorldEnd" << endl;
    cout << "FrameEnd" << endl;
}

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

  int seed = time(NULL);
  //printf("random seed is %d\n", seed);  
  srandom(seed);

  double total_volume, lava_volume, particle_volume;

  if (volume_is_rectangular) {
    total_volume = height*width*width;
  } else {
    // this is not right if container_radius is not constant
    total_volume = height*M_PI*container_radius(height/2.0)*container_radius(height/2.0);
  }

  lava_volume = lava_fraction*total_volume;
  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 the simulation  
  
  the_grid = new grid_t(-width/2.0,-width/2.0, 0.0,
			width/2.0, width/2.0, height,
			20, 20, 20);
  
  lava = create_lava(elements);
  
  if (argc == 2 && strcmp(argv[1], "-opengl") == 0) {
      
      // Initialize GLUT
      
      glutInit(&argc, argv);
      glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
      glutCreateWindow("Lava Lamp");
      glutDisplayFunc(display);
      glutReshapeFunc(resize);
      glutMouseFunc(mouse);
      glutMotionFunc(mouse_move);
      glutKeyboardFunc(key);
      glutIdleFunc(idle);
      
      // 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 light1_amb[] = { 0.1, 0.1, 0.1, 1.0 };
      GLfloat light1_dif[] = { 0.9, 0.9, 0.9, 1.0 };
      GLfloat light1_spec[] = { 1.0, 1.0, 1.0, 1.0 };
      GLfloat light1_pos[] = { 0.5, 0.5, 0.5, 0.0 };  
      
      if (use_lighting) {  
	  
	  // TODO: this is not right. 
	  
	  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_AMBIENT, light1_amb);
	  glLightfv(GL_LIGHT0, GL_DIFFUSE, light1_dif);
	  glLightfv(GL_LIGHT0, GL_SPECULAR, light1_spec);
	  glLightfv(GL_LIGHT0, GL_POSITION, light1_pos);
	  glEnable (GL_LIGHTING);
	  glEnable (GL_LIGHT0);
	  //glDepthFunc(GL_LESS);
	  //glEnable(GL_DEPTH_TEST);
      }
      
      // Switch to main loop
      glutMainLoop();
      free(lava);
      return 0;        
  }

  if (argc == 2) {
      int frameno = 1;
      cout << "Option \"searchpath\" \"shader\" \"/usr/local/bmrt/shaders:.:myshaders\"" << endl;
      while(frameno <= atoi(argv[1])) {
	  iterate(lava, elements, delta_t);
	  output_frame(lava, elements, frameno++);
      }
  }

  free(lava);
}