/* navier_gl.c     navier_gl  geom.dat
 *            see file   navier_gl.txt   for more information
 *            simple numerical solution of boundary value Navier-Stokes
 *            water problems. Solved in three dimensions with time stepping.
 *            Tested using two viscosities, basically for air and water.
 *            Water is assumed incompressible and air must be open to
 *            the standard atmosphere. My problems of interest have water
 *            from 1 centimeter to 10 meters deep with surface open to the
 *            standard atmosphere.  Falling or partially submerged solid
 *            objects may be moving in air or water
 *            (controlled by users 'move_object()').  
 *            A vacuum may be created due to separation.
 *            The model holds the elementary volumes fixed in space.
 *            The fluid may move from one element to an adjacent element
 *            at any time step. A fraction "full" is maintained in order
 *            to maintain conservation of mass.
 *            In general the initial state will have potential energy
 *            and no kinetic energy. An object dropped into calm water
 *            or a dam bursting or a slouch gate opening converts potential
 *            energy into kinetic energy and after a long time results in
 *            a state of only potential energy. The final potential energy
 *            is expected to be less than the initial potential energy
 *            because friction will have caused loss of energy due to heating. 
 */

/* may be compiled and linked on windows using   cl /GX /ML /I. navier_gl.c
 * may be compiled and linked on Linux using   gcc -o navier_gl navier_gl.c -lm
 * sample 4 by 4 by 4 test case executed by  navier_gl
 * run your problem:   navier_gl  your_problem.dat
 */

/*  
 * numerical approximations may fail for density and/or viscosity more
 * than a factor of two from that of water.  This is not intended to be
 * a general fluid dynamics program.
 * The program should handle transient solutions of a collapsing tower
 * of water.
 */

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

#undef  abs
#define abs(x) ((x)>0?(x):(-(x)))
#undef  sign
#define sign(x) ((x)<0.0?(-1.0):((x)>0.0?(1.0):(0.0)))
#undef  min
#define min(x,y) ((x)<(y)?(x):(y))

/* only MKS, metric units are used: Meter,  Kilogram, Second
 *                             for: Length, Mass,     Time
 *
 *                                            dimensionality
 * distance in meter                          L
 * velocity in meter/second                   L/T
 * acceleration in meter/second^2             L/TT
 * area in meter^2                            LL
 * volume in meter^3                          LLL
 * force in kilogram*meter/second^2           ML/TT
 * pressure in kilogram/(meter*second^2)      M/LTT
 * density in kilogram/meter^3                M/LLL
 * viscosity in kilogram/(meter*second)       M/LT
 *
 */

/* change parameters to suit your needs.                                     */
/* these may be read from a file                                             */

/* debug controls */
static int    debug = 1;       /* debug level, 0 is no debug                 */
static FILE * debug_handle;    /* for file navier_gl.debug output            */
static int    box_debug = 0;   /* special debug of a specific boxes          */
static int    box_d[5][3];     /* up to 5 x,y,z box coordinates              */
static double max_force;       /* for checking                               */   
static double max_accel;
static double max_veloc;
static double max_press;

/* simulation controls */
static double t        = 0.0;  /* time now in seconds                        */
static double dt       = 0.1;  /* present value of delta time in seconds     */
static int reduce_dt   = 0;    /* automatic scaling of time step             */
static int increase_dt = 0;    /* automatic scaling of time step             */
static double dt_user  = 0.1;  /* user, maximum dt                           */
static double dt_min   = 0.001;/* minimum dt, min(dx,dy,dz)/water_sound      */ 
static double t_max    = 9.9;  /* prevent run away                           */
static int    run      =   0;  /* toggled by key 'r'                         */
static double max_err  = 0.0;  /* present maximum error, internal check      */
static double avg_err  = 0.0;  /* present iteration average error            */
static double max_pres = 0.0;  /* max pressure for scaling plot              */
static char dat_file_name[255]; /* save file name for 'clear'                */

/* physical constants */
static double g               = 9.80665;  /* meter/second^2  acceleration    */
static double air_density     = 1.223;    /* kilogram/meter^3  15.7 'C       */
static double air_pressure    = 101.32;   /* kilogram/(meter*second^2) sea   */
static double air_viscosity   = 182.7E-6; /* kilogram/(meter*second) 18 'C   */
static double water_density   = 998.2;    /* kilogram/meter^3  20 'C         */
static double water_viscosity = 0.001066; /* kilogram/(meter*second)  20 'C  */
static double water_pressure  = 9690.93;  /* killogram/(meter*second^2) at 1m*/
static double water_mass      = 998.2;    /* kilogram for 1 meter cube       */
static double water_force     = 9690.0;   /* kilogram*meter/second^2 at 1m   */
static double water_sound     = 1461.0;   /* meter/second at 19 'C           */
static double solid_density   = 100000.0; /* kilogram/meter^3 user set       */
static double wall_viscosity  = 0.002;    /* kilogram/(meter*second) user set*/
static double Pi = 3.141592653589793238462643383279502884197;


/* bounding box for cells or nodes or incremental volumes                    */
static int nx = 2; /* default, over written by file read                     */
static int ny = 2;
static int nz = 2;

/* dimension of every incremental volume=dx * dy * dz, z is vertical(gravity)*/
static double dx = 0.01;     /* width,  meter                                */
static double dy = 0.01;     /* side,   meter                                */
static double dz = 0.01;     /* height, meter                                */
static double dV = 0.000001; /* dx * dy * dz volume, meter**3                */
static double dM = 0.0009982;/* dV * water mass                              */
static double dAxy = 0.001;  /* dx * dy area of top and bottom, meter**2     */
static double dAxz = 0.001;  /* dx * dz area of front and back, meter**2     */
static double dAyz = 0.001;  /* dy * dz area of sides,  meter**2             */

/*  primary data structure, a cell or node or element or incremental volume  */

#define md 64                    /* max dimension, changes size of 3D array  */
struct cell { int    type;       /* type of cell 0=water,1=air,2=vac,3=s,4=f */
              int    init_id;    /* initial id xxxyyyzzz                     */
              double force[3];   /* force on centroid, F=ma                  */
              double accel[3];   /* old [bx] [by] [bz]  new [cx] [cy] [cz]   */
              double veloc[6];   /* velocity of centroid V=integral accel dt */
              double posit[6];   /* delta position within element            */
              double full;       /* fraction full = fraction of mass         */
              double press[6];   /* pressure on six sides, matches adjacent  */
              int    side[6];    /* type of side 0=water,1=air,2=vac,3=s,4=f */
              char   status[4];  /* status flags                             */
              double density;    /* kilograms per meter**3                   */
              double viscosity;  /* water or air or solid                    */
              double color;      /* ignored for air, used for water mixing   */
              double error;      /* local error check                        */
            } box[md][md][md],   /* problem geometry in x, y, z elements     */
              water_init,        /* default water element, twat=0            */
              air_init,          /* default air element, tair=1              */
              vac_init,          /* default vacuum element, tvac=2           */
              solid_init,        /* solid object element, tsol=3             */
              fixed_init;        /* fixed (boundary) element, tfix=4         */
              
const int bx = 0; /* subscript of force, accel, veloc, posit  */
const int by = 1;
const int bz = 2;
const int cx = 3; /* new computed veloc  */
const int cy = 4;
const int cz = 5;

const int xp = 0; /* subscript of pressure and side type */
const int xm = 1; /* x minus (smaller coordinate) side of box */
const int yp = 2; /* y positive (larger) side of box */
const int ym = 3;
const int zp = 4;
const int zm = 5; /* z minus side of box */

const int twat = 0; /* water element type and side */
const int tair = 1; /* air element type and side */
const int tvac = 2; /* vacumme element type and side */
const int tsol = 3; /* solid element type and side, movable */
const int tfix = 4; /* fixed (boundary) element type and side */
                    /* equivalent to x=-1 or x=nx etc. bounding box */

/* OpenGL controls */
static int surfacew;        /* surface display window handle */
static int velocityw;       /* velocity vector display window handle */
static int pressurew;       /* pressure display window handle */
static int forcew;          /* force vector display window handle */
static int controlw;        /* control window handle */
static int reportw;         /* report window handle */
static float theta = -16.0; /* object rotation about X axis */
static float chi   = 0.0;   /* object rotation about Y axis */
static float phi   = -22.0; /* object rotation about Z axis */
static int width   = 300;   /* window width  float scale -1.0 to 1.0 */
static int height  = 300;   /* window height float scale -1.0 to 1.0 */
static int cwidth  = 200;   /* control width float scale 0.0 to 200.0 */
static int grid_on  = 0;    /* display grid points */
static int cells_on = 0;    /* display cells */
static int vectors_on = 0;  /* display vectors */
static int xmouse;          /* mouse x */
static int ymouse;          /* mouse y corrected */
static float xstep;         /* graphic x grid size 1.0/nx */ 
static float ystep;         /* graphic y grid size 1.0/ny */
static float zstep;         /* graphic z grid size 1.0/ny */
static float offset=-0.5;   /* volume centered in unit sphere, center */
static int hide_empty = 1;  /* hide water elements full < hew */
static double hew = 0.02;   /* hide empty water amount */
static double max_pressure; /* for plotting, the average side pressure */
static double min_pressure; /* for plotting, 12 steps between  */
GLfloat colors[][3] = {{1.0,0.0,0.0},{1.0,0.5,0.0},{1.0,0.8,0.0},
                       {1.0,1.0,0.0},{0.8,0.9,0.0},{0.5,0.9,0.0},
                       {0.0,1.0,0.0},{0.0,1.0,0.6},
                       {0.0,1.0,1.0},{0.0,0.8,1.0},{0.0,0.5,1.0},
                       {0.0,0.0,1.0}};
static int ncolor = 12;


/* function prototypes of internal functions */

static void read_geom(char * file_name);
static void init_box_types();
static void DrawSurface();
static void DrawPressure();
static void DrawVelocity();
static void DrawForce();
static void set_bound();
static void move_object();
static void other_side(const int side, int* o_side, int* ix, int* iy, int* iz);
static void compute_pressure();
static void compute_flow();
static void update_flow();
static void update_new_old();
static void update_location();
static void check_error();
static void debug_box(); 
static void print_box(int x, int y, int z); 
static void compute();

static void draw_box(float x1, float y1, float z1,
                     float x2, float y2, float z2);
static void line_box(float x1, float y1, float z1,
                     float x2, float y2, float z2);
static void writeText(GLfloat x, GLfloat y, GLfloat z,
                      GLfloat size, char * msg);
static void draw_grid(void);
static void draw_cells(void);
static void draw_vectors(void);
static void display_surface(void);
static void display_velocity(void);
static void display_pressure(void);
static void display_force(void);
static void display_control(void);
static void char_val(char *s, float x);
static void display_report(void);
static void update_displays(void);   /* all */
static void myinit(void);
static void mouse(int btn, int state, int x, int y);
static void mouse_control(int btn, int state, int x, int y);
static int in_rect(int x, int y, int w, int h);
static void keys(unsigned char key, int x, int y);
static void reshape(int w, int h);
static void controlReshape(int w, int h);
    /* int  main(int argc, char * argv[]); */
static void finish(void);


/*************** read_geom *******************/
static void read_geom(char * file_name)
{
  FILE * inp;
  char line[80];
  int i, len, x, y, z;
  int first_box = 0;

  if(debug) fprintf(debug_handle, "read_geom running \n");  
  box_debug = 0;
  hide_empty = 0;
  strcpy(dat_file_name, file_name);
  inp = fopen(dat_file_name, "r");
  if(inp==NULL)
  {
    printf("can not open file %s \n", file_name);
    exit(1);
  }
  init_box_types();
  while(!feof(inp))
  {
    fgets(line, 79, inp);
    len = strlen(line);
    if(len<1) continue;
    line[len-1] = '\0';
    if(debug) fprintf(debug_handle, "line = %s \n", line);
    if(line[0]=='#') continue; /* comment */
    if(line[0]=='n'&&line[1]=='x')      sscanf(&line[2], "%d", &nx);
    else if(line[0]=='n'&&line[1]=='y') sscanf(&line[2], "%d", &ny);
    else if(line[0]=='n'&&line[1]=='z') sscanf(&line[2], "%d", &nz);
    else if(line[0]=='d'&&line[1]=='x') sscanf(&line[2], "%lf", &dx);
    else if(line[0]=='d'&&line[1]=='y') sscanf(&line[2], "%lf", &dy);
    else if(line[0]=='d'&&line[1]=='z') sscanf(&line[2], "%lf", &dz);
    else if(line[0]=='d'&&line[1]=='t') sscanf(&line[2], "%lf", &dt);
    else if(line[0]=='t'&&line[1]=='t') sscanf(&line[2], "%lf", &t_max);
    else if(line[0]=='d'&&line[1]=='b') sscanf(&line[2], "%d", &debug);
    else if(line[0]=='h'&&line[1]=='e') sscanf(&line[2], "%lf", &hew);
    else if(line[0]=='h'&&line[1]=='i') sscanf(&line[2], "%d", &hide_empty);
    else if(line[0]=='w'&&line[1]=='a')
    {
      sscanf(&line[2], "%d %d %d", &x, &y, &z);
      box[x][y][z] = water_init;
      box[x][y][z].init_id = 1000000*x + 1000*y +z;
    }
    else if(line[0]=='a'&&line[1]=='i')
    {
      sscanf(&line[2], "%d %d %d", &x, &y, &z);
      box[x][y][z] = air_init;
      box[x][y][z].init_id = 1000000*x + 1000*y +z;
    }
    else if(line[0]=='v'&&line[1]=='a')
    {
      sscanf(&line[2], "%d %d %d", &x, &y, &z);
      box[x][y][z] = vac_init;
      box[x][y][z].init_id = 1000000*x + 1000*y +z;
    }
    else if(line[0]=='s'&&line[1]=='o')
    {
      sscanf(&line[2], "%d %d %d", &x, &y, &z);
      box[x][y][z] = solid_init;
      box[x][y][z].init_id = 1000000*x + 1000*y +z;
    }
    else if(line[0]=='f'&&line[1]=='i')
    {
      sscanf(&line[2], "%d %d %d", &x, &y, &z);
      box[x][y][z] = fixed_init;
      box[x][y][z].init_id = 1000000*x + 1000*y +z;
    }
    else if(line[0]=='b'&&line[1]=='o')
    {
      sscanf(&line[2], "%d %d %d", &box_d[box_debug][0],
                                   &box_d[box_debug][1],
                                   &box_d[box_debug][2]);
      box_debug++;
    }
    line[0] = '\0'; /* avoid junk on last line */
  } /* end while, get next line */
  fclose(inp);
  if(debug) fprintf(debug_handle,
                    "Geometry nx=%d, ny=%d, nz=%d, dx=%g, dy=%g, dz=%g \n",
                    nx, ny, nz, dx, dy, dz);
         
  /* initialize dependent variables */
  t = 0.0;
  dV = dx * dy * dz;    /* volume of element in meter^3 */
  dM = dV * water_mass; /* mass of full water element in Kg */
  dAxy = dx * dy;       /* area of top and bottom side =4,5 */
  dAxz = dx * dz;       /* area of side = 2,3 */
  dAyz = dy * dz;       /* area of side = 0,1 */
  dt_min = min(dx, min(dy, dz))/water_sound; /* minimum dt */
  xstep = 1.0/(float)nx;
  ystep = 1.0/(float)ny;
  zstep = 1.0/(float)nz;
  offset = -0.5; /*  xx = offset + ix * xstep */
  if(debug) fprintf(debug_handle,
                    "Geometry dV=%g, dM=%g, dAxy=%g, dAxz=%g, dAyz=%g \n",
                    dV, dM, dAxy, dAxz, dAyz);
  if(debug) fprintf(debug_handle, "debug=%d, dt=%g, t_max=%g, dt_min=%g \n",
                    debug, dt, t_max, dt_min);
  if(debug) fprintf(debug_handle,
                    "g=%g, air_density=%g, air_pressure=%g, air_viscosity=%g \n",
                    g, air_density, air_pressure, air_viscosity);
  if(debug) fprintf(debug_handle,
                    "water_density=%g, water_viscosity=%g, water_pressure=%g \n",
                    water_density, water_viscosity, water_pressure);
  if(debug) fprintf(debug_handle,
                    "water_mass=%g, water_force=%g, water_sound=%g \n", 
                    water_mass, water_force, water_sound); 
  if(debug) fprintf(debug_handle,
                    "solid_density=%g, wall_viscosity=%g, Pi=%g \n", 
                    solid_density, wall_viscosity, Pi); 
  if(debug) fprintf(debug_handle, "xstep=%f, ystep=%f, zstep=%f, offset=%f \n",
                    xstep, ystep, zstep, offset);
  
  /* set up element bounds (sides types) */
  set_bound();
  compute_pressure();
  if(debug)
  {
    fprintf(debug_handle, "read_geom exiting \n");
    debug_box();
  }
} /* end read_geom */

/*************** init_box_types *******************/
static void init_box_types()
{
  int i, x, y, z;

  if(debug) fprintf(debug_handle, "init_box_types running \n");
  
  air_init.type = tair;     /* air cell */
  air_init.density = air_density;
  air_init.viscosity = 0.0; /* flag, not real air_viscosity */
  air_init.color = 0.0;
  air_init.error = 0.0;
  air_init.init_id = 0;
  for(i=0; i<3; i++){
    air_init.force[i] = 0.0;
    air_init.accel[i] = 0.0;
  }
  air_init.full = 0.0;
  for(i=0; i<6; i++){
    air_init.press[i] = air_pressure;
    air_init.veloc[i] = 0.0;
    air_init.posit[i] = 0.0;
    air_init.side[i]  = 0;
  }
  for(i=0; i<4; i++){
    air_init.status[i] = 0;
  }

  vac_init = air_init;        /* vacuum cell */
  vac_init.type = tvac;
  vac_init.density = 0.0;     /* no mass */
  vac_init.viscosity = 0.0;
  for(i=0; i<6; i++){
    vac_init.press[i] = 0.0;
  }

  water_init = air_init;      /* water cell */
  water_init.type = twat;
  water_init.full = 1.0;
  water_init.density = water_density;
  water_init.viscosity = water_viscosity;
  for(i=0; i<6; i++){
    water_init.press[i] = 0.0;
  }
  /* pressure, force, etc. computed later */
  
  solid_init = air_init;     /* solid cell */
  solid_init.type = tsol;
  solid_init.density = solid_density;
  solid_init.viscosity = wall_viscosity;
  for(i=0; i<6; i++){
    solid_init.press[i] = 0.0;
  }
  
  fixed_init = air_init;     /* fixed cell, boundary */
  fixed_init.type = tfix;
  fixed_init.density = 0.0;     /* no mass */
  fixed_init.viscosity = wall_viscosity;
  for(i=0; i<6; i++){
    fixed_init.press[i] = 0.0;
  }
   
  /* initialize total volume to air */
  for(x=0; x<md; x++){
  for(y=0; y<md; y++){
  for(z=0; z<md; z++){
    box[x][y][z] = air_init;
    box[x][y][z].init_id = 1000000*x + 1000*y +z;
  }
  }
  }  
  if(debug) fprintf(debug_handle, "init_box_types exiting \n");  
} /* end init_box_types */

static void DrawSurface()
{
  float xx, yy, yyt, zz;
  int x, y, z;

  if(debug>4) fprintf(debug_handle, "in DrawSurface \n");
  for(x=0; x<nx; x++){
  for(y=0; y<ny; y++){
  for(z=0; z<nz; z++){
    if(box[x][y][z].type==tair) continue; /* do not draw air */
    if(hide_empty==1 && box[x][y][z].type==twat && 
       box[x][y][z].full<hew) continue; /* do not draw almost empty */ 
    xx = offset+x*xstep; /* get drawing lower left corner */
    yy = offset+y*ystep;
    zz = offset+z*zstep;
    if(box[x][y][z].type==twat)
    {
      /* water boxes blue and outline */
      glColor3f(0.0, 0.0, 1.0);
      draw_box(xx, yy, zz, xx+xstep, yy+ystep, zz+box[x][y][z].full*zstep);
      glColor3f(0.0, 0.0, 0.0);
      line_box(xx, yy, zz, xx+xstep, yy+ystep, zz+box[x][y][z].full*zstep);
    }
    else if(box[x][y][z].type==tvac)
    {
      /* outline vacuum */
      glColor3f(0.0, 0.0, 0.0);
      line_box(xx, yy, zz, xx+xstep, yy+ystep, zz+zstep);
    }
    else if(box[x][y][z].type==tsol)
    {
      /* red box */
      glColor3f(1.0, 0.0, 0.0);
      draw_box(xx, yy, zz, xx+xstep, yy+ystep, zz+zstep);
    }
    else if(box[x][y][z].type==tfix)
    {
      /* fixed cell black */
      glColor3f(0.0, 0.0, 0.0);
      draw_box(xx, yy, zz, xx+xstep, yy+ystep, zz+zstep);
    }
  }
  }
  }
} /* end DrawSurface, one Surface draw of screen */


/*************** DrawPressure *********************/
static void DrawPressure()
{
  float xx, yy, zz;
  double inc_pressure, dif_pressure, avg_pressure;
  int x, y, z;
  int i;

  if(debug>4) fprintf(debug_handle, "in DrawPressure \n");
  dif_pressure = max_pressure-min_pressure;
  
  for(x=0; x<nx; x++){
  for(y=0; y<ny; y++){
  for(z=0; z<nz; z++){

    xx = offset+x*xstep;
    yy = offset+y*ystep;
    zz = offset+z*zstep;

    if(box[x][y][z].type != twat) continue; /* draw only water */
    /* draw pressure, colorcode */
    avg_pressure = (box[x][y][z].press[xp]+box[x][y][z].press[xm]+
                    box[x][y][z].press[yp]+box[x][y][z].press[ym]+
                    box[x][y][z].press[zp]+box[x][y][z].press[zm])/6.0;
    inc_pressure = (max_pressure-avg_pressure)/dif_pressure;
    i = (int)(inc_pressure*(double)ncolor); /* expect truncation towards zero */
    if(i>ncolor-1) i = ncolor-1;
    if(i<0) i = 0;
    glColor3fv(colors[i]);
    draw_box(xx, yy, zz, xx+xstep, yy+ystep, zz+zstep);
    line_box(xx, yy, zz, xx+xstep, yy+ystep, zz+zstep);
  }
  }
  }
} /* end DrawPressure, one Pressure draw of screen */


/*************** DrawVelocity *************************/
static void DrawVelocity()
{
  float xx, yy, zz;
  int x, y, z;
  double vx, vy, vz, vv;
  double angx, angy, angz;

  if(debug>4) fprintf(debug_handle, "in DrawVelocity \n");
  glColor3f(1.0, 0.0, 0.0);
  for(x=0; x<nx; x++){
  for(y=0; y<ny; y++){
  for(z=0; z<nz; z++){

    xx = offset+x*xstep;
    yy = offset+y*ystep;
    zz = offset+z*zstep;

    if(box[x][y][z].type != twat) continue; /* draw only water */

    vx = box[x][y][z].veloc[bx];
    vy = box[x][y][z].veloc[by];
    vz = box[x][y][z].veloc[bz];
    vv = sqrt(vx*vx + vy*vy + vz*vz);
    if(vv > 0.0 || debug>0)
    {
      angx = 360.0*atan2(vx, vv)/Pi;
      angy = 360.0*atan2(vy, vv)/Pi;
      angz = 360.0*atan2(vz, vv)/Pi;
      glPushMatrix();
      glTranslatef(xx, yy, zz);
      glRotatef(angx, 1.0, 0.0, 0.0);
      glRotatef(angy, 0.0, 1.0, 0.0);
      glRotatef(angz, 0.0, 0.0, 1.0);
      glutSolidCone(0.03, 0.1, 6, 6);
      glPopMatrix();
    }
  }
  }
  }
} /* end DrawVelocity, one Velocity draw of screen */


/*************** DrawForce *************************/
static void DrawForce()
{
  float xx, yy, zz;
  int x, y, z;
  double fx, fy, fz, ff;
  double angx, angy, angz;

  if(debug>4) fprintf(debug_handle, "in DrawForce \n");
  glColor3f(1.0, 0.0, 1.0);
  for(x=0; x<nx; x++){
  for(y=0; y<ny; y++){
  for(z=0; z<nz; z++){

    xx = offset+x*xstep;
    yy = offset+y*ystep;
    zz = offset+z*zstep;

    if(box[x][y][z].type != twat) continue; /* draw only water */

    fx = box[x][y][z].force[bx];
    fy = box[x][y][z].force[by];
    fz = box[x][y][z].force[bz];
    ff = sqrt(fx*fx + fy*fy + fz*fz);
    if(ff > 0.0 || debug>0)
    {
      angx = 360.0*atan2(fx, ff)/Pi;
      angy = 360.0*atan2(fy, ff)/Pi;
      angz = 360.0*atan2(fz, ff)/Pi;
      glPushMatrix();
      glTranslatef(xx, yy, zz);
      glRotatef(angx, 1.0, 0.0, 0.0);
      glRotatef(angy, 0.0, 1.0, 0.0);
      glRotatef(angz, 0.0, 0.0, 1.0);
      glutSolidCone(0.03, 0.1, 6, 6);
      glPopMatrix();
    }
  }
  }
  }
} /* end DrawForce, one Force draw of screen */

/******************* set_bound *******************/
static void set_bound()
{
  /* type twat=0 for water, tair=1 for air, tvac=2, tsol=3, tfix=4 */
  /* xp = 0  x plus side, subscript of pressure and side type
   * xm = 1  x minus (smaller) side
   * yp = 2
   * ym = 3
   * zp = 4
   * zm = 5  z minus side */
  int x, y, z, o_side, ix, iy, iz;

  if(debug) fprintf(debug_handle, "in set_bound \n");
  if(debug) fprintf(debug_handle, "xp=%d, xm=%d, yp=%d, ym=%d, zp=%d, zm=%d \n",
                    xp, xm, yp, ym, zp, zm);
  if(debug) fprintf(debug_handle, "twat=%d, tair=%d, tvac=%d, tsol=%d, tfix=%d\n",
                    twat, tair, tvac, tsol, tfix);
  for(x=0; x<nx; x++){
  for(y=0; y<ny; y++){
  for(z=0; z<nz; z++){
    if(box[x][y][z].type != twat) continue;
    /* set side type, sides and bottom are walls, top open to air */
    if(x==nx-1) box[x][y][z].side[xp] = tfix;
    else
    {
      other_side(xp, &o_side, &ix, &iy, &iz);
      box[x][y][z].side[xp] = box[x+ix][y+iy][z+iz].type;
    }
    if(x==0   ) box[x][y][z].side[xm] = tfix;
    else
    {
      other_side(xm, &o_side, &ix, &iy, &iz);
      box[x][y][z].side[xm] = box[x+ix][y+iy][z+iz].type;
    }
    if(y==ny-1) box[x][y][z].side[yp] = tfix;
    else
    {
      other_side(yp, &o_side, &ix, &iy, &iz);
      box[x][y][z].side[yp] = box[x+ix][y+iy][z+iz].type;
    }
    if(y==0   ) box[x][y][z].side[ym] = tfix;
    else
    {
      other_side(ym, &o_side, &ix, &iy, &iz);
      box[x][y][z].side[ym] = box[x+ix][y+iy][z+iz].type;
    }
    if(z==nz-1) box[x][y][z].side[zp] = tair; /* top open to air */
    else
    {
      other_side(zp, &o_side, &ix, &iy, &iz);
      box[x][y][z].side[zp] = box[x+ix][y+iy][z+iz].type;
    }
    if(z==0   ) box[x][y][z].side[zm] = tfix;
    else
    {
      other_side(zm, &o_side, &ix, &iy, &iz);
      box[x][y][z].side[zm] = box[x+ix][y+iy][z+iz].type;
    }        
  }
  }
  }
  if(debug)
  {
    fprintf(debug_handle, "exiting set_bound \n");
  }
} /* end set_bound */


/********************* move_object *******************/
static void move_object() /* for use to move an object through water */
{

  if(debug>1)printf("move_object, t=%g \n", t);
  if(debug>1)fprintf(debug_handle, "move_object, t=%g \n", t);

  if(debug>1)fprintf(debug_handle, "exiting move_object \n");
} /* end move_object */


/******************** other_side ******************************/
static void other_side(const int side, int* o_side, int* ix, int* iy, int* iz)
{
  /* input your side number */
  /*        o_side ix iy iz
   * xp = 0   1     1  0  0  x positive subscript of pressure and side
   * xm = 1   0    -1  0  0
   * yp = 2   3     0  1  0
   * ym = 3   2     0 -1  0
   * zp = 4   5     0  0  1
   * zm = 5   4     0  0 -1  z negative */
   int to_side[6]={ 1,  0,  3,  2,  5,  4};
   int tix[6]    ={ 1, -1,  0,  0,  0,  0};
   int tiy[6]    ={ 0,  0,  1, -1,  0,  0};
   int tiz[6]    ={ 0,  0,  0,  0,  1, -1};

   *o_side=to_side[side];      /* return neighbors side number */
   *ix    =tix[side];          /* return offsets from your x,y,z to neighbor */
   *iy    =tiy[side];
   *iz    =tiz[side];
} /* end other_side */

/************************* compute_pressure *************/
static void compute_pressure()
{
  /*        o_side ix iy iz
   * xp = 0   1     1  0  0  x positive subscript of pressure and side
   * xm = 1   0    -1  0  0
   * yp = 2   3     0  1  0
   * ym = 3   2     0 -1  0
   * zp = 4   5     0  0  1
   * zm = 5   4     0  0 -1  z negative */
  int x, y, z;
  int ix, iy, iz;
  double avg_press;
  int side, o_side;

  if(debug>1)fprintf(debug_handle, "compute_pressure, t=%g \n", t);
  max_pressure = 0.0;
  min_pressure = 1.0E6;
  
  max_pres=0.0; /* for scaling pressure graph */
  for(z=nz-1; z>=0; z--){ /* work from top down */
    for(y=0; y<ny; y++){
      for(x=0; x<nx; x++){
        /* set initial pressure on all six sides, assumes z=nz-1 open to air */
        if(box[x][y][z].type==twat) /* water element */
        {
          /* all pressure measured positive out of each side of element */
          /* walls and solid objects exert equal and opposite pressure  */
          if(z==nz-1) box[x][y][z].press[zp] = air_pressure;
          if(z<nz-1)
          {
            box[x][y][z].press[zp] = box[x][y][z+1].press[zm];
            if(box[x][y][z].full<1.0) box[x][y][z].press[zp] = air_pressure;
          }
          box[x][y][z].press[zm] = box[x][y][z].press[zp] +
                                   g * box[x][y][z].density * 
                                   box[x][y][z].full * dz;
          if(box[x][y][z].press[zm] > max_pres) max_pres=box[x][y][z].press[zm];
          avg_press = (box[x][y][z].press[zp] + box[x][y][z].press[zm]) / 2.0;
          if(avg_press>max_pressure) max_pressure = avg_press;
          if(avg_press<min_pressure) min_pressure = avg_press;
          box[x][y][z].press[xm] = avg_press;
          box[x][y][z].press[xp] = avg_press;
          box[x][y][z].press[ym] = avg_press;
          box[x][y][z].press[yp] = avg_press;
          /* dynamic pressure */
          box[x][y][z].press[xp] = box[x][y][z].press[xp] + 
                                   sign(box[x][y][z].veloc[bx])*0.5*box[x][y][z].density*
                                   box[x][y][z].veloc[bx]*box[x][y][z].veloc[bx];
          box[x][y][z].press[xm] = box[x][y][z].press[xm] - 
                                   sign(box[x][y][z].veloc[bx])*0.5*box[x][y][z].density*
                                   box[x][y][z].veloc[bx]*box[x][y][z].veloc[bx];
          box[x][y][z].press[yp] = box[x][y][z].press[yp] + 
                                   sign(box[x][y][z].veloc[by])*0.5*box[x][y][z].density*
                                   box[x][y][z].veloc[by]*box[x][y][z].veloc[by];
          box[x][y][z].press[ym] = box[x][y][z].press[ym] - 
                                   sign(box[x][y][z].veloc[by])*0.5*box[x][y][z].density*
                                   box[x][y][z].veloc[by]*box[x][y][z].veloc[by];
          box[x][y][z].press[zp] = box[x][y][z].press[zp] + 
                                   sign(box[x][y][z].veloc[bz])*0.5*box[x][y][z].density*
                                   box[x][y][z].veloc[bz]*box[x][y][z].veloc[bz];
          box[x][y][z].press[zm] = box[x][y][z].press[zm] - 
                                   sign(box[x][y][z].veloc[bz])*0.5*box[x][y][z].density*
                                   box[x][y][z].veloc[bz]*box[x][y][z].veloc[bz];
        }
      }
    }
  }
  if(debug>1)fprintf(debug_handle, "exiting compute_pressure \n");
  if(debug>2)debug_box();
} /* end compute_pressure */

/************************* compute_flow *************/
static void compute_flow()
{
  int x, y, z;
  double fx, fy, fz;        /* temporary forces kilogram meter/second^2 */
  double Vx, Vy, Vz;        /* temporary velocities meter/second        */
  double ax, ay, az;        /* temporary accelerations meter/second^2   */
  double dPx, dPy, dPz;     /* signed difference in pressure            */
  double dPxp, dPxm, dPyp, dPym, dPzp, dPzm; /* side to neighbor side pressure difference */
  double dVyxp, dVzxp, dVyxm, dVzxm; /* side to neighbor side velocity difference */
  double dVxyp, dVzyp, dVxym, dVzym; /* side to neighbor side velocity difference */
  double dVxzp, dVyzp, dVxzm, dVyzm; /* side to neighbor side velocity difference */
  double Pxp, Pxm, Pyp, Pym, Pzp, Pzm, TPA;  /* adjusted side pressures */
  double muxp, muxm, muyp, muym, muzp, muzm; /* viscosity at face */
  double Sxp, Sxm, Syp, Sym, Szp, Szm, ST; /* fraction of mass moving */
  
  if(debug>1)fprintf(debug_handle, "compute_flow, t=%g \n", t);
  for(x=0; x<nx; x++){
  for(y=0; y<ny; y++){
  for(z=0; z<nz; z++){
    if(box[x][y][z].type != twat) continue;
    /* get difference of pressure with neighbor */
    if(x==nx-1) dPxp = 0.0; /* wall */
    else        dPxp = box[x][y][z].press[xp]-box[x+1][y][z].press[xm];
    if(x==0   ) dPxm = 0.0; /* wall */
    else        dPxm = box[x][y][z].press[xm]-box[x-1][y][z].press[xp];
    if(y==ny-1) dPyp = 0.0; /* wall */
    else        dPyp = box[x][y][z].press[yp]-box[x][y+1][z].press[ym];
    if(y==0   ) dPym = 0.0; /* wall */
    else        dPym = box[x][y][z].press[ym]-box[x][y-1][z].press[yp];
    if(z==nz-1) dPzp = 0.0; /* air */
    else        dPzp = box[x][y][z].press[zp]-box[x][y][z+1].press[zm];
    if(z==0   ) dPzm = 0.0; /* wall */
    else        dPzm = box[x][y][z].press[zm]-box[x][y][z-1].press[zp];  
    /* sum and distribute pressure due to external pressure */
    TPA = (dPxp + dPxm + dPyp + dPym + dPzp + dPzm)/6.0;
    Pxp = box[x][y][z].press[xp] + TPA;
    Pxm = box[x][y][z].press[xm] + TPA;
    Pyp = box[x][y][z].press[yp] + TPA;
    Pym = box[x][y][z].press[ym] + TPA;
    Pzp = box[x][y][z].press[zp] + TPA;
    Pzm = box[x][y][z].press[zm] + TPA;
    /* compute force, pressure difference */
    fx = (Pxp-Pxm)*dAyz;
    fy = (Pyp-Pym)*dAxz;
    fz = (Pzp-(Pzm-g*box[x][y][z].density*box[x][y][z].full*dz))*dAxy;
    if(abs(fz)<1.0E-8) fz = 0.0;
    /* add force due to viscous friction */
    /* get velocity difference with neighbor */
    if(x==nx-1) dVyxp = box[x][y][z].veloc[by]; /* wall */
    else        dVyxp = box[x][y][z].veloc[by]-box[x+1][y][z].veloc[by];
    if(x==nx-1) dVzxp = box[x][y][z].veloc[bz]; /* wall */
    else        dVzxp = box[x][y][z].veloc[bz]-box[x+1][y][z].veloc[bz];
    if(x==0   ) dVyxm = box[x][y][z].veloc[by]; /* wall */
    else        dVyxm = box[x][y][z].veloc[by]-box[x-1][y][z].veloc[by];
    if(x==0   ) dVzxm = box[x][y][z].veloc[bz]; /* wall */
    else        dVzxm = box[x][y][z].veloc[bz]-box[x-1][y][z].veloc[bz];
    if(y==ny-1) dVxyp = box[x][y][z].veloc[bx]; /* wall */
    else        dVxyp = box[x][y][z].veloc[bx]-box[x][y+1][z].veloc[bx];
    if(y==ny-1) dVzyp = box[x][y][z].veloc[bz]; /* wall */
    else        dVzyp = box[x][y][z].veloc[bz]-box[x][y+1][z].veloc[bz];
    if(y==0   ) dVxym = box[x][y][z].veloc[bx]; /* wall */
    else        dVxym = box[x][y][z].veloc[bx]-box[x][y-1][z].veloc[bx];
    if(y==0   ) dVzym = box[x][y][z].veloc[bz]; /* wall */
    else        dVzym = box[x][y][z].veloc[bz]-box[x][y-1][z].veloc[bz];
    if(z==nz-1) dVxzp = box[x][y][z].veloc[bx]; /* air */
    else        dVxzp = box[x][y][z].veloc[bx]-box[x][y][z+1].veloc[bx];
    if(z==nz-1) dVyzp = box[x][y][z].veloc[by]; /* air */
    else        dVyzp = box[x][y][z].veloc[by]-box[x][y][z+1].veloc[by];
    if(z==0   ) dVxzm = box[x][y][z].veloc[bx]; /* wall */
    else        dVxzm = box[x][y][z].veloc[bx]-box[x][y][z-1].veloc[bx];  
    if(z==0   ) dVyzm = box[x][y][z].veloc[by]; /* wall */
    else        dVyzm = box[x][y][z].veloc[by]-box[x][y][z-1].veloc[by];  
    /* get viscosity of face with neighbor */
    if(x==nx-1) muxp = wall_viscosity; /* wall */
    else        muxp = box[x+1][y][z].viscosity;
    if(x==0   ) muxm = wall_viscosity; /* wall */
    else        muxm = box[x-1][y][z].viscosity;
    if(y==ny-1) muyp = wall_viscosity; /* wall */
    else        muyp = box[x][y+1][z].viscosity;
    if(y==0   ) muym = wall_viscosity; /* wall */
    else        muym = box[x][y-1][z].viscosity;
    if(z==nz-1) muzp = air_viscosity; /* air */
    else        muzp = box[x][y][z+1].viscosity;
    if(z==0   ) muzm = wall_viscosity; /* wall */
    else        muzm = box[x][y][z-1].viscosity; 
    Vx = box[x][y][z].veloc[bx];
    Vy = box[x][y][z].veloc[by];
    Vz = box[x][y][z].veloc[bz];
    fx = fx + muyp*(Vx-dVxyp)*dy + muym*(Vx-dVxym)*dy +
              muzp*(Vx-dVxzp)*dz + muzm*(Vx-dVxzm)*dz;
    fy = fy + muxp*(Vy-dVyxp)*dx + muxm*(Vy-dVyxm)*dx +
              muzp*(Vy-dVyzp)*dz + muzm*(Vy-dVyzm)*dz;
    fz = fz + muxp*(Vz-dVzxp)*dx + muxm*(Vz-dVzxm)*dx +
              muyp*(Vz-dVzyp)*dy + muym*(Vz-dVzym)*dy;
    box[x][y][z].force[bx] = fx;
    box[x][y][z].force[by] = fy;
    box[x][y][z].force[bz] = fz;
    ax = fx/dM; 
    ay = fy/dM;
    az = fz/dM; 
    box[x][y][z].accel[bx] = ax;
    box[x][y][z].accel[by] = ay;
    box[x][y][z].accel[bz] = az;
    box[x][y][z].veloc[cx] = box[x][y][z].veloc[bx] + ax * dt;
    box[x][y][z].veloc[cy] = box[x][y][z].veloc[by] + ay * dt;
    box[x][y][z].veloc[cz] = box[x][y][z].veloc[bz] + az * dt;
    Sxp = ( Vx*dt + 0.5*(dPxp*dAyz/dM)*dt*dt)/dx;
    Sxm = (-Vx*dt + 0.5*(dPxm*dAyz/dM)*dt*dt)/dx;
    Syp = ( Vy*dt + 0.5*(dPyp*dAxz/dM)*dt*dt)/dy;
    Sym = (-Vy*dt + 0.5*(dPym*dAxz/dM)*dt*dt)/dy;
    Szp = ( Vz*dt + 0.5*(dPzp*dAxy/dM)*dt*dt)/dz;
    Szm = (-Vz*dt + 0.5*(dPzm*dAxy/dM)*dt*dt)/dz;
    ST  = Sxp + Sxm + Syp + Sym + Szp + Szm;
    /* reduce dt if necessary */
    if(abs(ST)>0.25) printf("moving too much ST=%g \n", ST);
    if(abs(Sxp)>0.16) printf("moving too much Sxp=%g \n", Sxp);
    if(abs(Sxm)>0.16) printf("moving too much Sxm=%g \n", Sxm);
    if(abs(Syp)>0.16) printf("moving too much Syp=%g \n", Syp);
    if(abs(Sym)>0.16) printf("moving too much Sym=%g \n", Sym);
    if(abs(Szp)>0.16) printf("moving too much Szp=%g \n", Szp);
    if(abs(Szm)>0.16) printf("moving too much Szm=%g \n", Szm);
    
    box[x][y][z].posit[xp] = Sxp; 
    box[x][y][z].posit[xm] = Sxm; 
    box[x][y][z].posit[yp] = Syp;
    box[x][y][z].posit[ym] = Sym; 
    box[x][y][z].posit[zm] = Szm; 
    box[x][y][z].posit[zp] = Szp;
    if(debug>4)
    {
      fprintf(debug_handle, "flow for x=%d, y=%d, z=%d, TPA=%g \n", x, y, z, TPA);
      fprintf(debug_handle, "dPxp=%g, dPxm=%g, dPyp=%g \n", dPxp, dPxm, dPyp);
      fprintf(debug_handle, "dPym=%g, dPzp=%g, dPzm=%g \n", dPym, dPzp, dPzm);
      fprintf(debug_handle, "Pxp=%g, Pxm=%g, Pyp=%g \n", Pxp, Pxm, Pyp);
      fprintf(debug_handle, "Pym=%g, Pzp=%g, Pzm=%g \n", Pym, Pzp, Pzm);
      fprintf(debug_handle, "dVxyp=%g, dVxym=%g, dVzyp=%g, dVzym=%g \n",
                             dVxyp, dVxym, dVzyp, dVzym);
      fprintf(debug_handle, "dVyxp=%g, dVyxm=%g, dVzxp=%g, dVzxm=%g \n",
                             dVyxp, dVyxm, dVzxp, dVzxm);
      fprintf(debug_handle, "dVxzp=%g, dVxzm=%g, dVyzp=%g, dVyzm=%g \n",
                             dVxzp, dVxzm, dVyzp, dVyzm);
      fprintf(debug_handle, "fx=%g, fy=%g, fz=%g \n", fx, fy, fz);
      fprintf(debug_handle, "ax=%g, ay=%g, az=%g \n", ax, ay, az);
      fprintf(debug_handle, "Vx=%g, Vy=%g, Vz=%g \n", box[x][y][z].veloc[cx],
                            box[x][y][z].veloc[cy], box[x][y][z].veloc[cz]);
      fprintf(debug_handle, "Sxp=%g, Sxm=%g, Syp=%g \n", box[x][y][z].posit[xp],
                            box[x][y][z].posit[xm], box[x][y][z].posit[yp]);
      fprintf(debug_handle, "Sym=%g, Szp=%g, Szm=%g \n", box[x][y][z].posit[ym],
                            box[x][y][z].posit[zp], box[x][y][z].posit[zm]);
    }
  }
  }
  }
  if(debug>1)fprintf(debug_handle, "exiting compute_flow \n");
}


/****************** update_flow *********************/
static void update_flow()
{
  int x, y, z;
  double fx, fy, fz; /* forces   kilogram meter/second**2 */
  double dm;         /* mass of incremental volume        */
  double ax, ay, az; /* accelerations  meter/second**2    */
  double dVx, dVy, dVz;     /* special case difference in velocity */
  double dPx, dPy, dPz;     /* special case difference in pressure */
  double dVxm, dVym, dVzm;  /* special case difference in velocity */
  double dVxp, dVyp, dVzp;  /* special case difference in velocity */

  if(debug>1)printf("update_flow, t=%g \n", t);
  if(debug>1)fprintf(debug_handle, "update_flow, t=%g \n", t);

  for(z=0; z<nz; z++){
  for(y=0; y<ny; y++){
  for(x=0; x<nx; x++){
    if(box[x][y][z].type!=twat) continue;
    /* have a water cell */
    if(box[x][y][z].posit[xp]>0.0) /* push water out */
    {
      box[x][y][z].full = box[x][y][z].full - box[x][y][z].posit[xp];
      if(x<nx-1)
      {
        box[x+1][y][z].full = box[x+1][y][z].full + box[x][y][z].posit[xp];
        if(box[x+1][y][z].type==tair || box[x+1][y][z].type==tvac)
        {
          box[x+1][y][z].type=twat;
          box[x+1][y][z].viscosity = water_viscosity;
          box[x+1][y][z].density = water_density;
        }
      }
    }
    if(box[x][y][z].posit[xm]>0.0) /* push water out */
    {
      box[x][y][z].full = box[x][y][z].full - box[x][y][z].posit[xm];
      if(x>0)
      {
        box[x-1][y][z].full = box[x-1][y][z].full + box[x][y][z].posit[xm];
        if(box[x-1][y][z].type==tair || box[x-1][y][z].type==tvac)
        {
          box[x-1][y][z].type=twat;
          box[x-1][y][z].viscosity = water_viscosity;
          box[x-1][y][z].density = water_density;
        }
      }
    }
    if(box[x][y][z].posit[yp]>0.0) /* push water out */
    {
      box[x][y][z].full = box[x][y][z].full - box[x][y][z].posit[yp];
      if(y<ny-1)
      {
        box[x][y+1][z].full = box[x][y+1][z].full + box[x][y][z].posit[yp];
        if(box[x][y+1][z].type==tair || box[x][y+1][z].type==tvac)
        {
          box[x][y+1][z].type=twat;
          box[x][y+1][z].viscosity = water_viscosity;
          box[x][y+1][z].density = water_density;
        }
      }
    }
    if(box[x][y][z].posit[ym]>0.0) /* push water out */
    {
      box[x][y][z].full = box[x][y][z].full - box[x][y][z].posit[ym];
      if(y>0)
      {
        box[x][y-1][z].full = box[x][y-1][z].full + box[x][y][z].posit[ym];
        if(box[x][y-1][z].type==tair || box[x][y-1][z].type==tvac)
        {
          box[x][y-1][z].type=twat;
          box[x][y-1][z].viscosity = water_viscosity;
          box[x][y-1][z].density = water_density;
        }
      }
    }
    if(box[x][y][z].posit[zp]>0.0) /* push water out */
    {
      box[x][y][z].full = box[x][y][z].full - box[x][y][z].posit[zp];
      if(z<nz-1)
      {
        box[x][y][z+1].full = box[x][y][z+1].full + box[x][y][z].posit[zp];
        if(box[x][y][z+1].type==tair || box[x][y][z+1].type==tvac)
        {
          box[x][y][z+1].type=twat;
          box[x][y][z+1].viscosity = water_viscosity;
          box[x][y][z+1].density = water_density;
        }
      }
    }
    if(box[x][y][z].posit[zm]>0.0) /* push water out */
    {
      box[x][y][z].full = box[x][y][z].full - box[x][y][z].posit[zm];
      if(z>0)
      {
        box[x][y][z-1].full = box[x][y][z-1].full + box[x][y][z].posit[zm];
        if(box[x][y][z-1].type==tair || box[x][y][z-1].type==tvac)
        {
          box[x][y][z-1].type=twat;
          box[x][y][z-1].viscosity = water_viscosity;
          box[x][y][z-1].density = water_density;
        }
      }
    }
  }
  }
  }
  if(debug>1) fprintf(debug_handle, "exiting update_flow \n");
} /* end update_flow */


/*************************** update_new_old ******************/
static void update_new_old()
{
  /* old [bx] [by] [bz]  new [cx] [cy] [cz] */
  int x, y, z;

  if(debug>1)printf("update_new_old, t=%g \n", t);
  if(debug>1)fprintf(debug_handle, "update_new_old, t=%g \n", t);

  for(z=0; z<nz; z++)
  {
    for(y=0; y<ny; y++)
    {
      for(x=0; x<nx; x++)
      {
        if(box[x][y][z].type!=twat) continue;
        /* have a water cell */

        box[x][y][z].veloc[bx] = box[x][y][z].veloc[cx];
        box[x][y][z].veloc[by] = box[x][y][z].veloc[cy];
        box[x][y][z].veloc[bz] = box[x][y][z].veloc[cz];
      }
    }
  }
  if(debug>1)fprintf(debug_handle, "exiting update_old_new \n");
}


/*************************** update_location ***************/
static void update_location()
{
  if(debug>1)printf("update_location, t=%g \n", t);
  if(debug>1)fprintf(debug_handle, "update_location, t=%g \n", t);
  /* need to use 'full' where possible to move fluid */
  /* some may change air->water, water->air, etc */


  if(debug>1)fprintf(debug_handle, "exiting update_location \n");
} /* end update location */


/******************* check_error *********************/
static void check_error()
{
  int x, y, z;
  double err = 0.0;
  int count = 1;

  if(debug>1)printf("check_error, t=%g \n", t);
  if(debug>1)fprintf(debug_handle, "check_error, t=%g \n", t);
  for(z=0; z<nz; z++)
  {
    for(y=0; y<ny; y++)
    {
      for(x=0; x<nx; x++)
      {
      }
    }
  }
  avg_err=avg_err/(double)count;
  if(debug>1)fprintf(debug_handle, "exiting check_error max=%g, avg=%g\n",
                     max_err, avg_err);
} /* end check_error */


/************************ debug_box ********************/
static void debug_box()
{
  int x, y, z;

  fprintf(debug_handle, "t=%g, dt=%g, dx=%g, dy=%g, dz=%g \n",
          t, dt, dx, dy, dz);
  for(x=0; x<nx; x++){
  for(y=0; y<ny; y++){
  for(z=0; z<nz; z++){
    if(box[x][y][z].type==tair) continue; /* no debug for air */
    print_box(x, y, z);
  }
  }
  }  
} /* end debug_box */

/***************************** print box for debugging ******************/
static void print_box(int x, int y, int z)
{
    fprintf(debug_handle,
   "   [%d,%d,%d] type=%d, id=%d, density=%g, viscosity=%g, color=%g, err=%g \n",
   x, y, z, box[x][y][z].type, box[x][y][z].init_id, box[x][y][z].density,
   box[x][y][z].viscosity, box[x][y][z].color, box[x][y][z].error);
    fprintf(debug_handle,
            "      force x=%g, y=%g, z=%g \n",
            box[x][y][z].force[bx], box[x][y][z].force[by],
            box[x][y][z].force[bz]);
    if(abs(box[x][y][z].force[bx])>max_force) max_force=box[x][y][z].force[bx];
    if(abs(box[x][y][z].force[by])>max_force) max_force=box[x][y][z].force[by];
    if(abs(box[x][y][z].force[bz])>max_force) max_force=box[x][y][z].force[bz];
    fprintf(debug_handle,
            "      accel x=%g, y=%g, z=%g \n",
            box[x][y][z].accel[bx], box[x][y][z].accel[by],
            box[x][y][z].accel[bz]);
    if(abs(box[x][y][z].accel[bx])>max_accel) max_accel=box[x][y][z].accel[bx];
    if(abs(box[x][y][z].accel[by])>max_accel) max_accel=box[x][y][z].accel[by];
    if(abs(box[x][y][z].accel[bz])>max_accel) max_accel=box[x][y][z].accel[bz];
    fprintf(debug_handle,
            "      veloc old vx=%g, vy=%g, vz=%g \n",
            box[x][y][z].veloc[bx], box[x][y][z].veloc[by],
            box[x][y][z].veloc[bz]);
    if(abs(box[x][y][z].veloc[bx])>max_veloc) max_veloc=box[x][y][z].veloc[bx];
    if(abs(box[x][y][z].veloc[by])>max_veloc) max_veloc=box[x][y][z].veloc[by];
    if(abs(box[x][y][z].veloc[bz])>max_veloc) max_veloc=box[x][y][z].veloc[bz];
    fprintf(debug_handle,
            "      veloc new vx=%g, vy=%g, vz=%g \n",
            box[x][y][z].veloc[cx], box[x][y][z].veloc[cy],
            box[x][y][z].veloc[cz]);
    if(abs(box[x][y][z].veloc[cx])>max_veloc) max_veloc=box[x][y][z].veloc[cx];
    if(abs(box[x][y][z].veloc[cy])>max_veloc) max_veloc=box[x][y][z].veloc[cy];
    if(abs(box[x][y][z].veloc[cz])>max_veloc) max_veloc=box[x][y][z].veloc[cz];
    fprintf(debug_handle,
            "      posit Sxp=%g, Sxm=%g, Syp=%g, Sym=%g, Szp=%g, Szm=%g \n",
            box[x][y][z].posit[xp], box[x][y][z].posit[xm],
            box[x][y][z].posit[yp], box[x][y][z].posit[ym],
            box[x][y][z].posit[zp], box[x][y][z].posit[zm]);
    fprintf(debug_handle,
            "      full  =%g, status[0]=%d, [1]=%d, [2]=%d, [3]=%d \n",
            box[x][y][z].full, box[x][y][z].status[0], box[x][y][z].status[1],
            box[x][y][z].status[2],box[x][y][z].status[3]);
    fprintf(debug_handle,
            "      press xp=%g, xm=%g, yp=%g, ym=%g, zp=%g, zm=%g \n",
            box[x][y][z].press[xp], box[x][y][z].press[xm],
            box[x][y][z].press[yp], box[x][y][z].press[ym],
            box[x][y][z].press[zp], box[x][y][z].press[zm]);
    if(abs(box[x][y][z].press[xp])>max_press) max_press=box[x][y][z].press[xp];
    if(abs(box[x][y][z].press[xm])>max_press) max_press=box[x][y][z].press[xm];
    if(abs(box[x][y][z].press[yp])>max_press) max_press=box[x][y][z].press[yp];
    if(abs(box[x][y][z].press[ym])>max_press) max_press=box[x][y][z].press[ym];
    if(abs(box[x][y][z].press[zp])>max_press) max_press=box[x][y][z].press[zp];
    if(abs(box[x][y][z].press[zm])>max_press) max_press=box[x][y][z].press[zm];
    fprintf(debug_handle,
            "      side type xp=%d, xm=%d, yp=%d, ym=%d, zp=%d, zm=%d \n",
            box[x][y][z].side[xp], box[x][y][z].side[xm],
            box[x][y][z].side[yp], box[x][y][z].side[ym],
            box[x][y][z].side[zp], box[x][y][z].side[zm]);
}

/***************************** main compute step ******************/
static void compute()
{
  int i;
  
  printf("main compute step, t=%g \n", t);
  if(debug) fprintf(debug_handle, "main compute step, t=%g \n", t);

  do{
    compute_flow();
    if(reduce_dt)
    {
      dt = dt/2.0;
      if(debug>0)printf("reduced dt, now dt=%g \n", dt);
      reduce_dt = 0;
      continue;
    }
    if(increase_dt)
    {
      dt = dt*2.0;
      if(debug>0)printf("increased dt, now dt=%g \n", dt);
      increase_dt = 0;
      continue;
    }
    update_flow();
    update_location();
    move_object();
    compute_pressure();
    check_error();
    for(i=0; i<box_debug; i++) print_box(box_d[i][0], box_d[i][1], box_d[i][2]);
    update_new_old();
    if(debug) debug_box();
    update_displays();
    t = t + dt;
    if(t > t_max) {printf("t_max=%g reached\n", t_max); break;}
  }while(run);

  if(debug>5)printf("exiting main compute, now t=%g \n", t);
  if(debug>1)fprintf(debug_handle, "exiting main compute, now t=%g \n", t);
  return ;
} /* end compute */


static void draw_box(float x1, float y1, float z1, float x2, float y2, float z2)
{
  /* six faces: x1,y1,z1 x1,y1,z2 x2,y1,z2 x2,y1,z1  */
  /*            x1,y1,z1 x1,y2,z1 x2,y2,z1 x2,y1,z1  */
  /*            x1,y1,z1 x1,y2,z1 x1,y2,z2 x1,y1,z2  */
  /*            x2,y2,z2 x2,y2,z1 x1,y2,z1 x1,y2,z2  */
  /*            x2,y2,z2 x2,y1,z2 x1,y1,z2 x1,y2,z2  */
  /*            x2,y2,z2 x2,y1,z2 x2,y1,z1 x2,y2,z1  */
  glBegin(GL_POLYGON);    /* set color before call */
    glVertex3f(x1,y1,z1);
    glVertex3f(x1,y1,z2);
    glVertex3f(x2,y1,z2);
    glVertex3f(x2,y1,z1);
  glEnd();
  glBegin(GL_POLYGON);
    glVertex3f(x1,y1,z1);
    glVertex3f(x1,y2,z1);
    glVertex3f(x2,y2,z1);
    glVertex3f(x2,y1,z1);
  glEnd();
  glBegin(GL_POLYGON);
    glVertex3f(x1,y1,z1);
    glVertex3f(x1,y2,z1);
    glVertex3f(x1,y2,z2);
    glVertex3f(x1,y1,z2);
  glEnd();
  glBegin(GL_POLYGON);
    glVertex3f(x2,y2,z2);
    glVertex3f(x2,y2,z1);
    glVertex3f(x1,y2,z1);
    glVertex3f(x1,y2,z2);
  glEnd();
  glBegin(GL_POLYGON);
    glVertex3f(x2,y2,z2);
    glVertex3f(x2,y1,z2);
    glVertex3f(x1,y1,z2);
    glVertex3f(x1,y2,z2);
  glEnd();
  glBegin(GL_POLYGON);
    glVertex3f(x2,y2,z2);
    glVertex3f(x2,y1,z2);
    glVertex3f(x2,y1,z1);
    glVertex3f(x2,y2,z1);
  glEnd();
} /* end draw_box */

static void line_box(float x1, float y1, float z1, float x2, float y2, float z2)
{
  /* six faces: x1,y1,z1 x1,y1,z2 x2,y1,z2 x2,y1,z1  */
  /*            x1,y1,z1 x1,y2,z1 x2,y2,z1 x2,y1,z1  */
  /*            x1,y1,z1 x1,y2,z1 x1,y2,z2 x1,y1,z2  */
  /*            x2,y2,z2 x2,y2,z1 x1,y2,z1 x1,y2,z2  */
  /*            x2,y2,z2 x2,y1,z2 x1,y1,z2 x1,y2,z2  */
  /*            x2,y2,z2 x2,y1,z2 x2,y1,z1 x2,y2,z1  */
  glBegin(GL_LINE_LOOP);    /* set color before call */
    glVertex3f(x1,y1,z1);
    glVertex3f(x1,y1,z2);
    glVertex3f(x2,y1,z2);
    glVertex3f(x2,y1,z1);
  glEnd();
  glBegin(GL_LINE_LOOP);
    glVertex3f(x1,y1,z1);
    glVertex3f(x1,y2,z1);
    glVertex3f(x2,y2,z1);
    glVertex3f(x2,y1,z1);
  glEnd();
  glBegin(GL_LINE_LOOP);
    glVertex3f(x1,y1,z1);
    glVertex3f(x1,y2,z1);
    glVertex3f(x1,y2,z2);
    glVertex3f(x1,y1,z2);
  glEnd();
  glBegin(GL_LINE_LOOP);
    glVertex3f(x2,y2,z2);
    glVertex3f(x2,y2,z1);
    glVertex3f(x1,y2,z1);
    glVertex3f(x1,y2,z2);
  glEnd();
  glBegin(GL_LINE_LOOP);
    glVertex3f(x2,y2,z2);
    glVertex3f(x2,y1,z2);
    glVertex3f(x1,y1,z2);
    glVertex3f(x1,y2,z2);
  glEnd();
  glBegin(GL_LINE_LOOP);
    glVertex3f(x2,y2,z2);
    glVertex3f(x2,y1,z2);
    glVertex3f(x2,y1,z1);
    glVertex3f(x2,y2,z1);
  glEnd();
} /* end line_box */

static void writeText(GLfloat x, GLfloat y, GLfloat z, GLfloat size, char * msg)
{
  char * p;
  
  glPushMatrix();
    glLoadIdentity ();
    glEnable(GL_LINE_SMOOTH);
    glEnable(GL_BLEND);
    glTranslatef(x, y, z);
    glScalef(size/1000.0, size/1000.0, 0.0);
    for(p=msg; *p; p++)
      glutStrokeCharacter(GLUT_STROKE_ROMAN, *p);
    glDisable(GL_LINE_SMOOTH);
    glDisable(GL_BLEND);
  glPopMatrix();
} /* end writeText, may have to reset glEnable's after call */

static void draw_grid(void)
{
  int i, j, k;
  float x, y, z;
  
  glColor3f (0.0, 0.0, 1.0);
  glPointSize(2.0);
  for(i=1; i<nx; i++){
  for(j=1; j<ny; j++){
  for(k=1; k<nz; k++){
    x = (float)i * xstep -0.5;
    y = (float)j * ystep -0.5;
    z = (float)k * zstep -0.5;
    glBegin(GL_POINTS);
      glVertex3f(x, y, z);
    glEnd();
  }
  }
  }
} /* end draw_grid */

static void double_grid(void)
{
  int x, y, z, xx, yy, zz;
  
  if(nx*2>md || ny*2>md || nz*2>md) return;
  
  for(x=nx-1; x>=0; x--)
  {
    xx = 2*x;
    for(y=0; y<ny; y++)
    {
      for(z=0; z<nz; z++)
      {
        box[xx+1][y][z] = box[x][y][z];
        box[xx][y][z] = box[x][y][z];
      }
    }
  }
  nx = nx * 2;
  for(y=ny-1; y>=0; y--)
  {
    yy = 2*y;
    for(x=0; x<nx; x++)
    {
      for(z=0; z<nz; z++)
      {
        box[x][yy+1][z] = box[x][y][z];
        box[x][yy][z] = box[x][y][z];
      }
    }
  }
  ny = ny * 2;
  for(z=nz-1; z>=0; z--)
  {
    zz = 2*z;
    for(y=0; y<ny; y++)
    {
      for(x=0; x<nx; x++)
      {
        box[x][y][zz+1] = box[x][y][z];
        box[x][y][zz] = box[x][y][z];
      }
    }
  }
  nz = nz * 2;  
  /* re-initialize dependent variables */
  dx = dx/2.0;
  dy = dy/2.0;
  dz = dz/2.0;
  if(debug) fprintf(debug_handle,
                    "Geometry nx=%d, ny=%d, nz=%d, dx=%g, dy=%g, dz=%g \n",
                    nx, ny, nz, dx, dy, dz);
         
  dV = dx * dy * dz;    /* volume of element in meter^3 */
  dM = dV * water_mass; /* mass of full water element in Kg */
  dAxy = dx * dy;       /* area of top and bottom side =4,5 */
  dAxz = dx * dz;       /* area of side = 2,3 */
  dAyz = dy * dz;       /* area of side = 0,1 */
  dt_min = min(dx, min(dy, dz))/water_sound; /* minimum dt */
  if(debug) fprintf(debug_handle, "dV=%g, dM=%g, dt_min=%g \n",
                    dV, dM, dt_min);
  if(debug) fprintf(debug_handle, "dAxy=%g, dAxz=%g, dAyz=%g \n",
                    dAxy, dAxz, dAyz);
  xstep = 1.0/(float)nx;
  ystep = 1.0/(float)ny;
  zstep = 1.0/(float)nz;
  offset = -0.5; /*  xx = offset + ix * xstep */
  if(debug) fprintf(debug_handle, "xstep=%f, ystep=%f, zstep=%f, offset=%f \n",
                    xstep, ystep, zstep, offset);
} /* end double_grid */

static void draw_cells(void)
{
  int i, j, k;
  float x, y, z;
  
  glColor3f (1.0, 0.0, 0.0);
  glPointSize(2.0);
  glBegin(GL_LINES);
    for(i=0; i<nx+1; i++){
    for(j=0; j<ny+1; j++){
      x = (float)i * xstep -0.5;
      y = (float)j * ystep -0.5;
      z = (float)0 * zstep -0.5;
      glVertex3f(x, y, z);
      glVertex3f(x, y, z+1.0);
    }
    }
    for(i=0; i<nx+1; i++){
    for(k=0; k<nz+1; k++){
      x = (float)i * xstep -0.5;
      y = (float)0 * ystep -0.5;
      z = (float)k * zstep-0.5;
      glVertex3f(x, y, z);
      glVertex3f(x, y+1.0, z);
    }
    }
    for(j=0; j<ny+1; j++){
    for(k=0; k<nz+1; k++){
      x = (float)0 * xstep -0.5;
      y = (float)j * ystep -0.5;
      z = (float)k * zstep -0.5;
      glVertex3f(x, y, z);
      glVertex3f(x+1.0, y, z);
    }
    }
  glEnd();
} /* end draw_cells */

static void draw_vectors(void)
{
  int i, j, k;
  float x, y, z;
  float ang = 0.0;
  
  glColor3f (0.0, 1.0, .0);
    for(i=0; i<nx; i++){
    for(j=0; j<ny; j++){
    for(k=0; k<nz; k++){
      x = (float)i * xstep -0.5 + xstep/2.0;
      y = (float)j * ystep -0.5 + ystep/2.0;
      z = (float)k * zstep -0.5 + zstep/2.0;
      glPushMatrix();
      glTranslatef(x, y, z);
      glRotatef(ang, 1.0, 0.0, 0.0);
      ang = ang+3.0; /* arbitrary rotation */
      glutSolidCone(0.03, 0.1, 6, 6);
      glPopMatrix();
    }
    }
    }
} /* end draw_vectors */

static void display_surface(void)
{
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glLoadIdentity();
  glRotatef(theta, 1.0, 0.0, 0.0);
  glRotatef(chi,   0.0, 1.0, 0.0);
  glRotatef(phi,   0.0, 0.0, 1.0);
  if(grid_on) draw_grid();
  DrawSurface();
  glColor3f (0.0, 0.0, 0.0);
  line_box(-0.5, -0.5, -0.5, 0.5, 0.5, 0.5);
  glColor3f(1.0, 0.0, 0.0);
  glPointSize(3.0);
  glBegin(GL_POINTS);
    glVertex3f(-0.5, -0.5, 0.55);
  glEnd();
  glFlush();
} /* end display_surface */

static void display_velocity(void)
{
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glLoadIdentity();
  glRotatef(theta, 1.0, 0.0, 0.0);
  glRotatef(chi,   0.0, 1.0, 0.0);
  glRotatef(phi,   0.0, 0.0, 1.0);
  if(cells_on) draw_cells();
  DrawVelocity();
  glColor3f (0.0, 0.0, 0.0);
  line_box(-0.5, -0.5, -0.5, 0.5, 0.5, 0.5);
  glFlush();
} /* end display_velocity */

static void display_pressure(void)
{
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glLoadIdentity();
  glRotatef(theta, 1.0, 0.0, 0.0);
  glRotatef(chi,   0.0, 1.0, 0.0);
  glRotatef(phi,   0.0, 0.0, 1.0);
  if(cells_on) draw_cells();
  if(vectors_on) draw_vectors();
  if(grid_on) draw_grid();
  DrawPressure();
  glColor3f (0.0, 0.0, 0.0);
  line_box(-0.5, -0.5, -0.5, 0.5, 0.5, 0.5);
  glFlush();
} /* end display_pressure */

static void display_force(void)
{
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glLoadIdentity();
  glRotatef(theta, 1.0, 0.0, 0.0);
  glRotatef(chi,   0.0, 1.0, 0.0);
  glRotatef(phi,   0.0, 0.0, 1.0);
  if(vectors_on) draw_vectors();
  DrawForce();
  glColor3f (0.0, 0.0, 0.0);
  line_box(-0.5, -0.5, -0.5, 0.5, 0.5, 0.5);
  glFlush();
} /* end display_force */

static void display_control(void)
{
  float fh;

  fh = (float)height;
  glClear(GL_COLOR_BUFFER_BIT);
  glLoadIdentity();
  glColor3f(0.0, 0.0, 0.0);
  writeText(5.0, fh- 25.0, 0.0, 100.0, "key commands");
  writeText(5.0, fh- 45.0, 0.0, 100.0, "x,X,y,Y,z,Z rotate");
  writeText(5.0, fh- 65.0, 0.0, 100.0, "c  outlines cells");
  writeText(5.0, fh- 85.0, 0.0, 100.0, "v  vector direction");
  writeText(5.0, fh-105.0, 0.0, 100.0, "g  grid points");

  glColor3f(0.0, 1.0, 1.0);
  glBegin(GL_POLYGON);
    glVertex3f(0.0,0.0,0.0);
    glVertex3f((float)cwidth,0.0,0.0);
    glVertex3f((float)cwidth,20.0,0.0);
    glVertex3f(0.0,20.0,0.0);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  writeText(5.0, 3.0, 0.0, 150.0, "clear");

  glColor3f(1.0, 1.0, 0.0);
  glBegin(GL_POLYGON);
    glVertex3f(0.0,22.0,0.0);
    glVertex3f((float)cwidth,22.0,0.0);
    glVertex3f((float)cwidth,42.0,0.0);
    glVertex3f(0.0,42.0,0.0);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  writeText(5.0, 25.0, 0.0,150.0, "exit");

  glColor3f(0.0, 1.0, 1.0);
  glBegin(GL_POLYGON);
    glVertex3f(0.0,44.0,0.0);
    glVertex3f((float)cwidth,44.0,0.0);
    glVertex3f((float)cwidth,64.0,0.0);
    glVertex3f(0.0,64.0,0.0);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  writeText(5.0, 47.0, 0.0, 150.0, "run");

  glColor3f(1.0, 1.0, 0.0);
  glBegin(GL_POLYGON);
    glVertex3f(0.0,66.0,0.0);
    glVertex3f((float)cwidth,66.0,0.0);
    glVertex3f((float)cwidth,86.0,0.0);
    glVertex3f(0.0,86.0,0.0);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  writeText(5.0, 69.0, 0.0,150.0, "stop");

  glColor3f(0.0, 1.0, 1.0);
  glBegin(GL_POLYGON);
    glVertex3f(0.0,88.0,0.0);
    glVertex3f((float)cwidth,88.0,0.0);
    glVertex3f((float)cwidth,108.0,0.0);
    glVertex3f(0.0,108.0,0.0);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  writeText(5.0, 91.0, 0.0, 150.0, "step");

  glColor3f(1.0, 1.0, 0.0);
  glBegin(GL_POLYGON);
    glVertex3f(0.0,110.0,0.0);
    glVertex3f((float)cwidth,110.0,0.0);
    glVertex3f((float)cwidth,130.0,0.0);
    glVertex3f(0.0,130.0,0.0);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  writeText(5.0, 113.0, 0.0,150.0, "finer grid");

  glColor3f(0.0, 1.0, 1.0);
  glBegin(GL_POLYGON);
    glVertex3f(0.0,132.0,0.0);
    glVertex3f((float)cwidth,132.0,0.0);
    glVertex3f((float)cwidth,152.0,0.0);
    glVertex3f(0.0,152.0,0.0);
  glEnd();
  glColor3f(0.0, 0.0, 0.0);
  writeText(5.0, 135.0, 0.0,150.0, "smaller dt");

  glFlush();
} /* end display_control */

static void char_val(char *s, float x)
{
  sprintf(s, "%4.0f", x);
  *(s+4)=' ';
}

static void display_report(void)
{
  char t_val[]     = "T  =         seconds";
  char dt_val[]    = "dT =         seconds";
  char theta_val[] = "X rotation =      degrees";
  char chi_val[]   = "Y rotation =      degrees";
  char phi_val[]   = "Z rotation =      degrees";
  
  glClear(GL_COLOR_BUFFER_BIT);
  glLoadIdentity();
  glColor3f(0.0, 0.0, 0.0);
  sprintf(&t_val[5], "%7.3f", t);
  t_val[12]=' '; /* remove null from sprintf string */
  writeText(5.0, (float)(height-30), 0.0, 100.0, t_val);
  sprintf(&dt_val[5], "%7.3f", dt);
  dt_val[12]=' '; /* remove null from sprintf string */

  writeText(5.0, (float)(height-50), 0.0, 100.0, dt_val);
  char_val(&theta_val[13], theta);
  writeText(5.0, 45.0, 0.0, 100.0, theta_val);
  char_val(&chi_val[13], chi);
  writeText(5.0, 25.0, 0.0, 100.0, chi_val);
  char_val(&phi_val[13], phi);
  writeText(5.0,  5.0, 0.0, 100.0, phi_val);
  glFlush();
} /* end display_report */

static void update_displays(void)   /* all */
{
  glutSetWindow(surfacew);   /* draw surface */
  glutPostRedisplay(); 
  glutSetWindow(velocityw);  /* draw velocity */
  glutPostRedisplay();
  glutSetWindow(pressurew);  /* draw pressure */
  glutPostRedisplay(); 
  glutSetWindow(forcew);     /* draw force */
  glutPostRedisplay();
  glutSetWindow(controlw);   /* draw control window */
  glutPostRedisplay(); 
  glutSetWindow(reportw);    /* draw report window */
  glutPostRedisplay(); 
} /* end update_displays */

static void myinit(void) /* shared */
{
  glClearColor (1.0, 1.0, 1.0, 1.0);
  glColor3f (0.0, 0.0, 0.0);
  xstep = 1.0/(float)nx;
  ystep = 1.0/(float)ny;
  zstep = 1.0/(float)nz;
  offset = -0.5; /*  xx = offset + ix * xstep */
  if(debug>6)fprintf(debug_handle, "xstep=%f, ystep=%f, zstep=%f, offset=%f \n",
                     xstep, ystep, zstep, offset);
} /* end myinit */

static void mouse(int btn, int state, int x, int y) /* shared */
{
  if(btn==GLUT_LEFT_BUTTON && state==GLUT_DOWN)
  {
    theta = theta + 2.0;
    if(theta>360.0) theta = theta - 360.0;
    if(theta<-360.0) theta = theta + 360.0;
  }
  if(btn==GLUT_RIGHT_BUTTON && state==GLUT_DOWN)
  {
    phi = phi + 2.0;
    if(phi>360.0) phi = phi - 360.0;
    if(phi<-360.0) phi = phi + 360.0;
  }
  update_displays();
} /* end mouse */

static void mouse_control(int btn, int state, int x, int y)
{
  xmouse = x; /* global mouse coordinates */
  ymouse = height-y;
  if(btn==GLUT_LEFT_BUTTON && state==GLUT_DOWN)
  {
    if(in_rect(0, 0, cwidth, 20)) /* clear */
    {
      theta = -16.0;
      chi   = 0.0;
      phi   = -22.0;
      grid_on = 0;
      cells_on = 0;
      vectors_on = 0;
      run = 0;
      /* re-read input */
      read_geom(dat_file_name);
    }
    else if(in_rect(0, 22, cwidth, 20)) /* exit */
    {
      finish(); /* save debug info */
    }
    else if(in_rect(0, 44, cwidth, 20)) /* run */
    {
      run = 1;
      compute();
    }
    else if(in_rect(0, 66, cwidth, 20)) /* stop */
    {
      run = 0;
    }
    else if(in_rect(0, 88, cwidth, 20)) /* step */
    {
      run = 0;
      compute();
    }
    else if(in_rect(0, 110, cwidth, 20)) /* finer grid */
    {
      double_grid();
    }
    else if(in_rect(0, 132, cwidth, 20)) /* smaller dt */
    {
      dt = dt/2.0;
    }
    else if(in_rect(0, 195, cwidth, 20)) /* grid on */
    {
      grid_on = 1 - grid_on;
    }
    else if(in_rect(0, 215, cwidth, 20)) /* vectors on */
    {
      vectors_on = 1 - vectors_on;
    }
    else if(in_rect(0, 235, cwidth, 20)) /* cells on */
    {
      cells_on = 1 - cells_on;
    }
  }
  update_displays();
} /* end mouse_control */

static int in_rect(int x, int y, int w, int h)
{
  return xmouse>x && xmouse<x+w && ymouse>y && ymouse<y+h;
} /* end in_rect */

static void keys(unsigned char key, int x, int y)
{
  /* Use x, X, y, Y, z, Z to rotate about axis */
  if(key == 'x'){theta -= 2.0; if(theta<-360.0) theta = theta + 360.0;}                
  if(key == 'X'){theta += 2.0; if(theta>360.0) theta = theta - 360.0;}
  if(key == 'y'){chi   -= 2.0; if(chi<-360.0) chi = chi + 360.0;}
  if(key == 'Y'){chi   += 2.0; if(chi>360.0) chi = chi - 360.0;}
  if(key == 'z'){phi   -= 2.0; if(phi<-360.0) phi = phi + 360.0;}
  if(key == 'Z'){phi   += 2.0; if(phi>360.0) phi = phi - 360.0;}
  if(key == 'g') grid_on = 1 - grid_on;
  if(key == 'c') cells_on = 1 - cells_on;
  if(key == 'v') vectors_on = 1 - vectors_on;
  if(key == 'q') finish(); /* quit = exit */
  if(key == 'e') finish(); /* exit = quit */
  if(key == 'r') run = 1 - run; /* run */
  if(key == 's') { run = 0; compute();} ; /* step */
  update_displays();
} /* end keys */

static void reshape(int w, int h) /* shared */
{
  width = w;
  height = h;
  glViewport(0, 0, w, h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  if (w <= h) 
    glOrtho (-1.0, 1.0, -1.0*(GLfloat)h/(GLfloat)w, 
             1.0*(GLfloat)h/(GLfloat)w, -1.0, 1.0);
  else 
    glOrtho (-1.0*(GLfloat)w/(GLfloat)h, 
              1.0*(GLfloat)w/(GLfloat)h, -1.0, 1.0, -1.0, 1.0);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity ();
  glEnable(GL_DEPTH_TEST);
} /* end reshape */

static void controlReshape(int w, int h)
{
  glViewport(0, 0, w, h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluOrtho2D(0.0, (float)cwidth, 0.0, (float)height);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity ();
} /* end controlReshape */

/*********************************** main  navier_gl ********************/
int main(int argc, char* argv[])
{
  printf("navier_gl starting \n");
  debug_handle = fopen("navier_gl.debug", "w");
  fprintf(debug_handle, "navier_gl starting\n");
  if(argc>1) read_geom(argv[1]);
  else       read_geom("navier_gl.dat");
  
  glutInit(&argc,argv);
  glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB | GLUT_DEPTH);

  glutInitWindowPosition(25,50);
  glutInitWindowSize(cwidth,height);
  controlw=glutCreateWindow("control");
  myinit();
  glutDisplayFunc(display_control);
  glutReshapeFunc(controlReshape);
  glutMouseFunc(mouse_control);
  glutKeyboardFunc(keys);

  glutInitWindowPosition(25,50+height+34);
  glutInitWindowSize(cwidth,height);
  reportw=glutCreateWindow("report");
  myinit();
  glutDisplayFunc(display_report);
  glutReshapeFunc(controlReshape);
  glutMouseFunc(mouse);
  glutKeyboardFunc(keys);

  glutInitWindowPosition(25+cwidth+8,50);
  glutInitWindowSize(width,height);
  surfacew=glutCreateWindow("Surface");
  myinit();
  glutDisplayFunc(display_surface); 
  glutReshapeFunc(reshape);
  glutMouseFunc(mouse);
  glutKeyboardFunc(keys);
  /* glEnable(GL_DEPTH_TEST); */

  glutInitWindowPosition(25+cwidth+8+width+8,50);
  glutInitWindowSize(width,height);
  velocityw=glutCreateWindow("Velocity vectors");
  myinit();
  glutDisplayFunc(display_velocity);
  glutReshapeFunc(reshape);
  glutMouseFunc(mouse);
  glutKeyboardFunc(keys);

  glutInitWindowPosition(25+cwidth+8,50+height+34);
  glutInitWindowSize(width,height);
  pressurew=glutCreateWindow("Pressure");
  myinit();
  glutDisplayFunc(display_pressure); 
  glutReshapeFunc(reshape);
  glutMouseFunc(mouse);
  glutKeyboardFunc(keys);

  glutInitWindowPosition(25+cwidth+8+width+8,50+height+34);
  glutInitWindowSize(width,height);
  forcew=glutCreateWindow("Force vectors");
  myinit();
  glutDisplayFunc(display_force);
  glutReshapeFunc(reshape);
  glutMouseFunc(mouse);
  glutKeyboardFunc(keys);

  /* user must initiate "run" or "step" */
  
  printf("going into glutMainLoop \n"); 
  glutMainLoop();
  return 0; /* exit should be via "quit" to finish() */
} /* end main */

static void finish(void)
{
  if(debug)
  { 
    fprintf(debug_handle, "navier_gl finished\n");
    fflush(debug_handle);
    fclose(debug_handle);
  }
  printf("navier_gl finished\n");
  exit(0);
} /* end finish */