#include "main.h"

/*
**	Global Variables
*/

double eps_dist = 0.00001;
double dot_tol = 0.9;

Vector ll = {-1.5, -1.5, -1.5};
Vector ur = {1.5, 1.5, 1.5};
double nx = 2, ny = 2, nz = 2;


main(argc, argv)        
     char **argv;
{
  vol_scan(ll.x, ll.y, ll.z, ur.x, ur.y, ur.z
	   ,(ur.x-ll.x)/nx, (ur.y-ll.y)/ny, (ur.z-ll.z)/nz);
}


vol_scan(x0, y0, z0, x1, y1, z1, xinc, yinc, zinc)
	double x0, y0, z0, x1, y1, z1, xinc, yinc, zinc;
{
  Vertex v[8];
  int k, side, hit;
  double x, y, z;

  for (z = z0; z <= z1; z += zinc) {
    for (y = y0; y <= y1; y += yinc) {
      for (x = x0; x <= x1; x += xinc) {
	for (k = 0, hit = FALSE; k < 8; k++) {
	  v[k].p.x = (k & 01) ? x : x - xinc;
	  v[k].p.y = (k & 02) ? y : y - yinc;
	  v[k].p.z = (k & 04) ? z : z - zinc;

	  v[k].fval = f_value(v[k].p);
	  if (k == 0) 
	    side = sign(v[0].fval);
	  else if (side != sign(v[k].fval))
	    hit = TRUE;
	}
	if (hit) { 
	  simplex(&(v[0]), &(v[1]) ,&(v[3]) , &(v[7]));
	  simplex(&(v[0]), &(v[5]) ,&(v[1]) , &(v[7]));
	  simplex(&(v[0]), &(v[3]) ,&(v[2]) , &(v[7]));
	  simplex(&(v[0]), &(v[2]) ,&(v[6]) , &(v[7]));
	  simplex(&(v[0]), &(v[4]) ,&(v[5]) , &(v[7]));
	  simplex(&(v[0]), &(v[6]) ,&(v[4]) , &(v[7]));
	}
      }
    }
  }
}


double f_value(p)
     Vector p;
{
  double r2 = 1.0;

  return (sqrt(p.x*p.x + p.y*p.y + p.z*p.z) - r2);
}


Vector f_ngrad(p)
     Vector p;
{
  return g3_scalev(1.0 / g3_norm(p), p);
}


int sign(v)
     double v;
{
  return (v < 0.0)? -1.0 : 1.0;
}