// system includes
#include <math.h>
#include <iostream.h>
#include <stdlib.h>

// openGL includes
#include <GL/gl.h>
#include <GL/glu.h>

// gprims utility library: need definition of max_flt, etc.
#include "gprims/vector.H"

// local includes
#include "drawutils.H"
#include "object.H"


//////////////////////////////////////////////////////////////////////////
//
// GeomObject -- base virtual class
//
//////////////////////////////////////////////////////////////////////////
GeomObject::GeomObject( void ) 
{
  _pos[0] = 0; _pos[1] = 0; _pos[2] = 0;
  randomColor ( &_col[0] );  
  _constraint = NULL;
}

void       
GeomObject::setPosition( const Point3& newPos ) 
{  
  _pos = newPos;
}

 
//////////////////////////////////////////////////////////////////////////
//
// RandomMeshscape  (Jon Heiner & Marc Lebovitz)
//
// A Mesh "landscape" of width, height & depth of size on plane of constant z.
// It has nribs rows & columns of randomly (??) generated height.
//
//////////////////////////////////////////////////////////////////////////

// all functions with jjj-debug have not been changed from the XYRect functions.
// thus they will ... not work!

float
RandomMeshscape::hypotenuse(float x, float y) {
  float mag = powf(x,2) + powf(y,2);
  return sqrt(mag);
}


RandomMeshscape::RandomMeshscape( float s, int r ) :
  _size( s ),
  _nribs( r )
{
  debug = 1;
  debug_checkheight = 0;
  debug_checkbox = 1;
  debug_intersect_triangle = 0;
  debug_intersectWithRay = 1;
  debug_setPixel = 0;

  boxes = new int **[_nribs];
  // set boxes to 0's
  for (int q=0; q<_nribs; q++) {
    boxes[q] = new int *[_nribs];
    for (int r=0; r<_nribs; r++) {
      boxes[q][r] = new int[5]; // last 2 are for tri1/tri2 intersection
      boxes[q][r][0] = 0;
      boxes[q][r][1] = 0;
      boxes[q][r][2] = 0;
    }
  }

  float rnd, dist, val;
  const float rndness = .75;
  const float mtnness = .25;
  // create the depth matrix
  _depth = new float *[ r+1 ];           // depth is a 2d array of heights
  _inc = (float) s / (float) r;
  srand(13);                            // random number seed

  for( int i=0; i<r+1; i++) {
    _depth[i] = new float [ r+1 ];

    for( int j=0; j<r+1; j++) {          // bias depth to objects closer to the center
      rnd = (float) 1.0/(log10f(rand()));
      dist = 1.0f / hypotenuse((float) i-(float) r/2.0f + 0.5f, (float) j-(float) r/2.0f + 0.5f);

      val = (float) (0.5 * s * (mtnness*dist + rndness*rnd) / (mtnness + rndness));
      if (val > s) val = s;            // handle cases of infinite depth
      _depth[i][j] = val;
    }
  }
  _squaremins = new float *[ r ];
  _squaremaxs = new float *[ r ];
  float *a = new float;
  float *b = new float;
  float *c = new float;
  float *d = new float;
  float max,min;
  
  for( i=0; i<r; i++) {
    _squaremaxs[i] = new float [r];
    _squaremins[i] = new float [r];
    
    for( int j=0; j<r; j++) {          // bias depth to objects closer to the center
      getSquareHeights(i,j, a, b, c, d);
      if(*a>*b) { 
	max = *a; 
	min = *b;
	  }
      else {
	min = *a; 
	max= *b;
      }
      if (*c > max) max = *c;
      if (*c < min) min = *c;
      if (*d > max) max = *d;
      if (*d < min) min = *d;
      
      _squaremins[i][j] = min;
      _squaremaxs[i][j] = max;
    }
  }
}

/* assigns vals of 1st float[3] to be those of 2nd */
void
RandomMeshscape::assign_3vec(float* dest, float* src){
  dest[0] = src[0];
  dest[1] = src[1];
  dest[2] = src[2];
}

/* returns 1 if pt1 is before or at pt2 along ray from eye->far */
int
RandomMeshscape::first(float* pt1, float* pt2, const Point3& eye, const Point3& far) {
  // we know pt1,pt2 are on the ray, so just check with x values
  float t1,t2;
  float a = far[0] - eye[0];
  if (a!=0) { // use x's
    float b = eye[0];
    t1 = (pt1[0] - b) / a;
    t2 = (pt2[0] - b) / a;
  }
  else {
    float c = far[1] - eye[1];
    if (c!=0) { // use y's
      float d = eye[1];
      t1 = (pt1[1] - d) / c;
      t2 = (pt2[1] - d) / c;
    }
    else {
      float e = far[2] - eye[2];
      if (e!=0) { // use z's
	float f = eye[2];
	t1 = (pt1[2] - f) / e;
	t2 = (pt2[2] - f) / e;
      }
      else { // sanity check
	cout << "error in RandomMeshscape::first: ";
	cout << "eye, far are the same point" << endl;
      }
    }
  }
  if (t1<=t2) return 1;
  else return 0;
}

/* returns the magnitude of vector vec */
float
RandomMeshscape::magnitude(float* vec) {
  float mag = powf(vec[0],2) + powf(vec[1],2) + powf(vec[2],2);
  mag = sqrtf(mag);
  return mag;
}

/* normalizes the vector vec
   -> modifies argument */
void
RandomMeshscape::normalize(float* vec) {
  float mag = magnitude(vec);
  vec[0] = vec[0] / mag;
  vec[1] = vec[1] / mag;
  vec[2] = vec[2] / mag;

  // sanity check
  mag = magnitude(vec) - 1.0F;
  if (abs(mag) > 0.01) cout << "error in normalize" << endl;
}

/* cross product of 3-vecs */
void
RandomMeshscape::cross_product(float* v1, float* v2, float* cross) {
  cross[0] = v1[1]*v2[2] - v2[1]*v1[2];
  cross[1] = v1[2]*v2[0] - v2[2]*v1[0];
  cross[2] = v1[0]*v2[1] - v2[0]*v1[1];
}

/* returns the area of the triangle defined by
   clockwise points: pt1,pt2,pt3 */
// not being used!
float
RandomMeshscape::area(float* pt1, float* pt2, float* pt3) {
  float vecA[3];
  vecA[0] = pt2[0] - pt1[0];
  vecA[1] = pt2[1] - pt1[1];
  vecA[2] = pt2[2] - pt1[2];
  float vecB[3];
  vecB[0] = pt3[0] - pt1[0];
  vecB[1] = pt3[1] - pt1[1];
  vecB[2] = pt3[2] - pt1[2];
  float cross[3];
  cross_product(vecA, vecB, cross);
  float mag = magnitude(cross);
  mag = mag/2;
  return mag;
}

/* returns area of 2D triangle a,b,c
   where points are listed counterclockwise */
float
RandomMeshscape::area2(float* a, float* b, float* c) {
  float area;
  area = a[0] * b[1] - a[1] * b[0] +
    a[1] * c[0] - a[0] * c[1] +
    b[0] * c[1] - c[0] * b[1];
  return area;
}

/* returns 1 if pt3 is left of line pt1->pt2
   returns 1 if pt3 is on line pt1->pt2
   else returns 0
   note: is solving in 2D (x,y) because of nature of height field */
int
RandomMeshscape::left_on(float* pt1, float* pt2, float* pt3, int dir) {
  float a[2];
  float b[2];
  float c[2];
  if (dir == 0) {
    a[0] = pt1[1];
    a[1] = pt1[2];
    b[0] = pt2[1];
    b[1] = pt2[2];
    c[0] = pt3[1];
    c[1] = pt3[2];
  }
  if (dir == 1) {
    a[0] = pt1[0];
    a[1] = pt1[2];
    b[0] = pt2[0];
    b[1] = pt2[2];
    c[0] = pt3[0];
    c[1] = pt3[2];
  }
  if (dir == 2) {
    a[0] = pt1[0];
    a[1] = pt1[1];
    b[0] = pt2[0];
    b[1] = pt2[1];
    c[0] = pt3[0];
    c[1] = pt3[1];
  }

  return (area2(a,b,c) >= 0);
}

/* if intersectPt is inside triangle pt1,pt2,pt3, returns 1
   else returns 0
   dir gives the general direction of eye/far ray...
   0->i, 1->j, 2->k
   project triangle problem onto plane of other 2 directions
*/
int
RandomMeshscape::inside_triangle(float* pt1,float* pt2,float* pt3,float* intersectPt, int dir) {

   if ( left_on(pt1,pt2,intersectPt,dir) &&
	left_on(pt2,pt3,intersectPt,dir) &&
	left_on(pt3,pt1,intersectPt,dir) )
    return 1;
  else return 0;
}

/* if eye->far ray intersects plane defined by pt1,pt2,pt3, sets
   intersectPt to the point of intersection and returns 1
   else returns 0 */
int
RandomMeshscape::intersect_plane(float* pt1, float* pt2, float* pt3, float* intersectPt,
		const Point3& e, const Point3& f) {
  // find plane A,B,C,D through pt1,pt2,pt3
  // line1 = pt2 - pt1
  float line1[3];
  line1[0] = pt2[0] - pt1[0];
  line1[1] = pt2[1] - pt1[1];
  line1[2] = pt2[2] - pt1[2];
  // line2 = pt3 - pt1
  float line2[3];
  line2[0] = pt3[0] - pt1[0];
  line2[1] = pt3[1] - pt1[1];
  line2[2] = pt3[2] - pt1[2];
  // normal A,B,C to plane = line1xline2
  float normal[3];
  cross_product(line1,line2,normal);
  normalize(normal);
  // to find D, plug pt1 into equation Ax+By+Cz+D=0 and solve for D
  float A,B,C,D;
  A=normal[0];
  B=normal[1];
  C=normal[2];
  D = -1 * ( A*pt1[0] + B*pt1[1] + C*pt1[2] );
  // sanity checks on plane equation
  float check = A*pt2[0] + B*pt2[1] + C*pt2[2] + D;
  if (abs(check) > 0.01) cout << "error in plane equation...pt2" << endl;
  check = A*pt3[0] + B*pt3[1] + C*pt3[2] + D;
  if (abs(check) > 0.01) cout << "error in plane equation...pt3" << endl;

  // intersect plane A,B,C,D with line through e,f
  float num = A*e[0] + B*e[1] + C*e[2] + D;
  float denom = A*(e[0] - f[0]) + B*(e[1] - f[1]) + C*(e[2] - f[2]);
  if (denom==0) return 0; // line || to plane or inside it -> 0 or inf solns
  float t = num/denom;
  // sanity check
  if (t<0) { // pt of intersect before beginning of ray...should never happen!
    cout << "you should never read this! ...see intersect_plane" << endl;
    return 0;
  }

  // line: p = e + t*(f-e)
  intersectPt[0] = e[0] + t*(f[0] - e[0]);
  intersectPt[1] = e[1] + t*(f[1] - e[1]);
  intersectPt[2] = e[2] + t*(f[2] - e[2]);

  // sanity check
  check = A*intersectPt[0] + B*intersectPt[1] + C*intersectPt[2] + D;
  if (abs(check) > 0.01) cout << "error in plane equation...intersectPt not on plane" << endl;


  return 1;
}

/* find major direction (i,j,or k) of eye/far ray
   returns: 0 for i, 1 for j, 2 for k */
int
RandomMeshscape::find_major_dir(const Point3& eye, const Point3& far) {
  float a,c,e;
  a=far[0]-eye[0]; c=far[1]-eye[1]; e=far[2]-eye[2];
  int max = 0;
  float maxval = a;
  if (c>maxval) {
    max = 1;
    maxval = c;
  }
  if (e>maxval) {
    max = 2;
    maxval = e;
  }
  return max;
}

/* if eye->far ray intersects the triangle defined by pt1,pt2,pt3,
   returns 1 and gives the point of intersection in intersectPt
   else returns 0 */
int
RandomMeshscape::intersect_triangle(float* pt1, float* pt2, float* pt3, float* intersectPt,
				    const Point3& eye, const Point3& far) {
  // remember to convert e/f to object space by subtracting _pos!
  if (intersect_plane(pt1,pt2,pt3,intersectPt,eye,far)) {
    if (debug_intersect_triangle) cout << "inside plane" << endl;
    int dir = find_major_dir(eye,far);
    // make copy of intersectPt st inside_triangle doesn't alter it
    //  is this necessary?
    float ipt[3];
    assign_3vec(ipt, intersectPt);
    if (inside_triangle(pt1,pt2,pt3,ipt, dir)) {
      if (debug_intersect_triangle) cout << "inside triangle" << endl;
      return 1;
    }
    else return 0;
  }
  else return 0;
}

/* if eye/far ray intersects box, returns 1 and intersectPt is set to 1st point of intersection along ray
 else returns 0 */
int
RandomMeshscape::checkbox(int x, int y, float* intersectPt, const Point3& eye, const Point3& far) {
  int tri1 = 0;
  int tri2 = 0;
  float pt1[3];
  float pt2[3];
  float bleft[3]; 
  float bright[3];
  float tleft[3];
  float tright[3];

  if (checkheight(x,y,eye,far)) { // if the ray isn't way out of z-range
    boxes[x][y][1] = 1;
    if (debug_checkbox) cout << "within height bounds" << endl;
    getSquareVertices(x,y,bleft,bright,tleft,tright); // ?? does jon's code work? ??
    if (intersect_triangle(bleft, bright, tleft, pt1, eye, far)) { // must order triangle counterclockwise
      if (debug_checkbox) cout << "intersects triangle 1" << endl;
      tri1 = 1;
      };
    if (intersect_triangle(bright, tright, tleft, pt2, eye, far)) { // must order triangle counterclockwise
      if (debug_checkbox) cout << "intersects triangle 2" << endl;
      tri2 = 1;
      };
  }
  else return 0; // the ray was way out of z-range
  if (!tri1 && !tri2) return 0; // intersects neither
  if (tri1 && !tri2) assign_3vec(intersectPt,pt1);
  if (tri2 && !tri1) assign_3vec(intersectPt,pt2);
  if (tri1 && tri2) { // find the first pt along the ray and return it
    if (first(pt1,pt2,eye,far)) assign_3vec(intersectPt,pt1);
    else assign_3vec(intersectPt,pt2);
  }
  boxes[x][y][2] = 1;
  return 1;
}

/* returns 1 if eye->far ray is within bounding box of square x,y
   else returns 0 */
int
RandomMeshscape::checkheight(int i, int j, const Point3& eye, const Point3& far) {
  // first, let's parameterize eye/far ray
  // L=eye+(far-eye)t
  //   x=at+b, y=ct+d, z=et+f
  float a,b,c,d,e,f;
  b=eye[0]; d=eye[1]; f=eye[2];
  a=far[0]-eye[0]; c=far[1]-eye[1]; e=far[2]-eye[2];

  // second, let's get the bounds for the box
  float xl,xr,yb,yt;
  xl = getSquareXLeft(i);
  xr = getSquareXRight(i);
  yb = getSquareYBottom(j);
  yt = getSquareYTop(j);

  if (debug_checkheight) {
    cout << "xl: " << xl << endl;
    cout << "xr: " << xr << endl;
    cout << "yb: " << yb << endl;
    cout << "yt: " << yt << endl;
  }

  // sanity checks
  if (yb >= yt) cout << "error: RandomMeshscape::checkheight: yb>=yt" << endl;
  if (xl >= xr) cout << "error: RandomMeshscape::checkheight: xl>=xr" << endl;

  /*            third, let's check where ray intersects box edges                 */
  float t,x,y,z1,z2;
  int found = 0;

  if (a!=0) { // otherwise is vertical line
    // check if e/f intersects on left edge of box
    t = (xl - b) / a;
    y = c*t + d;
    if (debug_checkheight) {
    cout << "x: " << xl << endl;
    cout << "y: " << y << endl;
    cout << "t: " << t << endl;
    }

    // sanity checks
    if ((y>yt) && (y<yb)) cout << "error in checkheight: ((y>yt) && (y<yb))" << endl;
    // debug checks
    if (debug_checkheight) {
      if (y>yt) cout << "left edge intersect: y too high" << endl;
      if (y<yb) cout << "left edge intersect: y too low" << endl;
    }

    if ( (y<=yt) && (y>=yb) ) {
      // x,y within box at this t
      if (debug_checkheight) cout << "intersected left edge" << endl;
      if (t<0) if (debug_checkheight) cout << "error: intersection of ray with box before eye!" << endl;
      z1 = e*t + f;
      if (debug_checkheight) cout << "z1: " << z1 << endl;
      found++;
    }
    
    // check if e/f intersects on right edge of box
    t = (xr - b) / a;
    y = c*t + d;
    if (debug_checkheight) {
      cout << "x: " << xr << endl;
      cout << "y: " << y << endl;
      cout << "t: " << t << endl;
    }

    // sanity checks
    if ((y>yt) && (y<yb)) cout << "error in checkheight: ((y>yt) && (y<yb))" << endl;
    // debug checks
    if (debug_checkheight) {
      if (y>yt) cout << "left edge intersect: y too high" << endl;
      if (y<yb) cout << "left edge intersect: y too low" << endl;
    }

    if ( (y<=yt) && (y>=yb) ) {
      // x,y within box at this t
      if (debug_checkheight) cout << "intersected right edge" << endl;
      if (t<0) if (debug_checkheight) cout << "error: intersection of ray with box before eye!" << endl;
      if (found == 0) z1 = e*t + f;
      if (found == 1) z2 = e*t + f;
      if (debug_checkheight) {
	cout << "z1: " << z1 << endl;
	cout << "z2: " << z2 << endl;
      }
      found++;
    }
  }

  if (c!=0) { // otherwise is horizontal line {
    if (found < 2) {
      // check if e/f intersects on top edge of box
      t = (yt - d) / c;
      x = a*t + b;
      if (debug_checkheight) {
	cout << "x: " << x << endl;
	cout << "y: " << yt << endl;
	cout << "t: " << t << endl;
      }

    // sanity checks
    if ((x>xr) && (x<xl)) cout << "error in checkheight: ((x>xr) && (x<xl))" << endl;
    // debug checks
    if (debug_checkheight) {
      if (x>xr) cout << "left edge intersect: x too high" << endl;
      if (x<xl) cout << "left edge intersect: x too low" << endl;
    }

      if ( (x<=xr) && (x>=xl) ) {
	// x,y within box at this t
	if (debug_checkheight) cout << "intersected top edge" << endl;
	if (t<0) if (debug_checkheight) cout << "error: intersection of ray with box before eye!" << endl;
	if (found == 0) z1 = e*t + f;
	if (found == 1) z2 = e*t + f;
	if (debug_checkheight) {
	cout << "z1: " << z1 << endl;
	cout << "z2: " << z2 << endl;
	}
	found++;
      }
    }
    
    if (found < 2) {
      // check if e/f intersects on bottom edge of box
      t = (yb - d) / c;
      x = a*t + b;
      if (debug_checkheight) {
      cout << "x: " << x << endl;
      cout << "y: " << yb << endl;
      cout << "t: " << t << endl;
      }

    // sanity checks
    if ((x>xr) && (x<xl)) cout << "error in checkheight: ((x>xr) && (x<xl))" << endl;
    // debug checks
    if (debug_checkheight) {
    if (x>xr) cout << "left edge intersect: x too high" << endl;
    if (x<xl) cout << "left edge intersect: x too low" << endl;
    }

      if ( (x<=xr) && (x>=xl) ) {
	// x,y within box at this t
	if (debug_checkheight) cout << "intersected bottom edge" << endl;
	if (t<0) if (debug_checkheight) cout << "error: intersection of ray with box before eye!" << endl;
	z2 = e*t + f;
	if (debug_checkheight) cout << "z2: " << z2 << endl;
	found++;
      }
    }
  }

  // sanity check
  if ((found != 2) && (found != 0)) cout << "error: must intersect twice or not at all in checkHeight" << endl;


  // fourth, let's find max/min z's of the 2 ray intersections with the box edges
  float maxz,minz;
  if (z1>z2) {
    maxz = z1;
    minz = z2;
  }
  else {
    maxz = z2;
    minz = z1;
  }

  if (debug_checkheight) {
  cout << "maxz: " << maxz << endl;
  cout << "minz: " << minz << endl;
  cout << "square min: " << _squaremins[i][j] << endl;
  cout << "square max: " << _squaremaxs[i][j] << endl;
  }

  // fifth (finally), let's compare ray intersection heights with box heights
  // if box heights totally above ray, no intersection
  if (_squaremins[i][j] > maxz) return 0;
  // if box heights totally below ray, no intersection
  if (_squaremaxs[i][j] < minz) return 0;
  // intersection may exist
  return 1;
}

/* returns square where line y=mx+b is at value x=0 and x=width of mesh
   whether square is validSquare or not */
void
RandomMeshscape::getBoundingSquares(int* square1, int* square2, float m, float b)
{
  float y;
  y = b;
  square1[0] = 0;
  square1[1] = floor((y/_size)* (float)_nribs);
  y = m*_size + b;
  square2[0] = _nribs - 1;
  square2[1] = floor((y/_size)* (float)_nribs);
}

/* returns true if square x,y is a valid part of the mesh */
int
RandomMeshscape::validSquare(int x,int y) {
  // x,y can have values 0 --> (_nribs-1)
  if ( (x < _nribs) &&
       (x >= 0) &&
       (y < _nribs) &&
       (y >= 0) ) return 1;
  else return 0;
}



// jjj-debug
int        
RandomMeshscape::intersectWithRay ( const Point3& eye, const Point3& far, 
			   Point3& ipt ) {
  _intersect = 0;

  if (debug) {
    // set boxes to 0's
    for (int q=0; q<_nribs; q++) {
      for (int r=0; r<_nribs; r++) {
	boxes[q][r][0] = 0;
	boxes[q][r][1] = 0;
	boxes[q][r][2] = 0;
      }
    }
  }

  float intersectPt[3];
  float newPt[3]; // new intersection point
  int do_bres = 1;

  // translate eye,far to object space
  float eyer[3];
  float farr[3];
  eyer[0] = eye[0] - _pos[0];
  eyer[1] = eye[1] - _pos[1];
  eyer[2] = eye[2] - _pos[2];
  farr[0] = far[0] - _pos[0];
  farr[1] = far[1] - _pos[1];
  farr[2] = far[2] - _pos[2]; 

  if (debug_intersectWithRay) {
    cout << "in world coordinates: " << endl;
    cout << "  eye: " << eye[0] << " " << eye[1] << " " << eye[2] << endl;
    cout << "  far: " << far[0] << " " << far[1] << " " << far[2] << endl;
    cout << "  mesh pos: " << _pos[0] << " " << _pos[1] << " " << _pos[2] << endl;

    float v0,v1,v2;
    v0 = far[0] - eye[0];
    v1 = far[1] - eye[1];
    v2 = far[2] - eye[2];

    eyex = eye[0];
    eyey = eye[1];
    eyez = eye[2];
    farx = far[0] + 40*v0;
    fary = far[1] + 40*v1;
    farz = far[2] + 40*v2;
  }

  if (debug_intersectWithRay) cout << "eye/far line info (projected into z=0 plane): " << endl;
  float A,B,C,m,b;
  m=0.0;b=0.0; // default vals
  A = eyer[1] - farr[1];
  B = farr[0] - eyer[0];
  C = eyer[0]*farr[1] - eyer[1]*farr[0];
  if (B==0) { // vertical line
    if (debug_intersectWithRay) cout << "  vertical line" << endl;
    // sanity check
    if (A==0) cout << "vert line error in intersectWithRay" << endl;
    int xv = floor(( (-C/A) /_size) * _nribs);
    for (int yv = 0; yv < _nribs; yv++) {
      setPixel(xv, yv, ipt, newPt, eyer, farr);
      if (debug && validSquare(xv,yv)) boxes[xv][yv][0] = 1;
    }
    do_bres = 0;
  }
  if ((A!=0) && (B!=0)) { // do bresenham
    if (debug_intersectWithRay) cout << "  candidate for bresenham" << endl;
    m = -A/B;
    b = -C/B;
    if (debug_intersectWithRay) cout << "  2D slope,intercept: " << m << " " << b << endl;
  }
  
  int square1[2];
  int square2[2];
  getBoundingSquares(square1,square2,m,b);  // can i display to debug?
  int x1 = square1[0];
  int y1 = square1[1];
  int x2 = square2[0];
  int y2 = square2[1];
  bx1 = x1;
  by1 = y1;
  bx2 = x2;
  by2 = y2;

  if (y1 == y2) { // horizontal enough to avoid bresenham
    if (debug_intersectWithRay) cout << "  horizontal enough to avoid bresenham" << endl;
    // sanity check
    if (B==0) cout << "horiz line error in intersectWithRay" << endl;
    int yh = floor(( (-C/B) /_size) * _nribs);
    for (int xh = 0; xh < _nribs; xh++) {
      setPixel(xh, yh, ipt, newPt, eyer, farr);
      if (debug && validSquare(xh,yh)) boxes[xh][yh][0] = yh;
    }
    do_bres = 0;
  }

  // sanity check
  if ((x1 == x2) && (y1 == y2)) cout << "error in bounding squares ...see intersectWithRay" << endl;

  if (debug_intersectWithRay) {
    cout << "line through boxes info in 2D:" << endl;
    cout << "  endpt1: " << x1 << " " << y1 << endl;
    cout << "  endpt2: " << x2 << " " << y2 << endl;
    cout << "  m: " << m << endl;
    cout << "  b: " << b << endl;
    segx0=0;
    segy0=m*(segx0)+b;
    segx1=_size;
    segy1=m*(segx1)+b;

    // sanity checks on eye/far 2D projection
    float test;
    test = A*eye[0] + B*eye[1] +C;
    if (abs(test) > 0.01) cout << "  eye/far 2D projection wrong...A,B,C don't agree with eye" << endl;
    test = A*far[0] + B*far[1] +C;
    if (abs(test) > 0.01) cout << "  eye/far 2D projection wrong...A,B,C don't agree with far" << endl;
    test = m*eye[0] + b - eye[1];
    if (abs(test) > 0.01) cout << "  eye/far 2D projection wrong...m,b don't agree with eye" << endl;
    test = m*far[0] + b - far[1];
    if (abs(test) > 0.01) cout << "  eye/far 2D projection wrong...m,b don't agree with far" << endl;
    if ((m==0.0)&&(b==0.0)) cout << "m and b were never set" << endl;
  }

  if (do_bres) {
    do_bresenham(m, b, x1, y1, x2, y2, ipt, newPt, eyer, farr);
  }

  // translate intersection point out of object space
  ipt[0] = ipt[0] + _pos[0];
  ipt[1] = ipt[1] + _pos[1];
  ipt[2] = ipt[2] + _pos[2];

  if (_intersect) {
    if (debug) cout << "ray intersects" << endl;
    return 1;
  }
  if (debug) cout << "ray does not intersect" << endl;
  return 0;
}

/* runs augmented bresenham and computes intersection on chosen boxes */
void
RandomMeshscape::do_bresenham(float m, float b, int x1, int y1, int x2, int y2, float* ipt, float* newPt, float* eyer, float* farr) {
  int x, y, dx, dy, npix;
  int xinc, yinc, xup, yup; // advance step; error step
  int e, einc, emax;        // error; err inc; err limit
  
  dx = x2-x1;  dy = y2-y1;  // x,y deltas and sign bits
  xinc = (dx >= 0) ? 1 : -1;
  yinc = (dy >= 0) ? 1 : -1;
  
  if ( (xinc * dx) > (yinc * dy) ) { // loop across x
    npix = (dx * xinc) + 1;          // note:  non-neg
    emax = dx * xinc;
    einc = (dy * yinc) << 1;
    xup = 0;
    yup = yinc;
    yinc = 0; // adv, err steps
  }
  else {                             // loop across y
    npix = (dy * yinc) + 1;          // note:  non-neg
    emax = dy * yinc;
    einc = (dx * xinc) << 1;
    yup = 0;
    xup = xinc;
    xinc = 0; // adv, err steps
  }
  
  x = x1;
  y = y1;
  e = 0;    // error initially zero
  while ( npix-- > 0 ) {
    setPixel(x, y, ipt, newPt, eyer, farr);
    if (debug && validSquare(x,y)) boxes[x][y][0] = 1;
    x += xinc;
    y += yinc; 
    e += einc; // advance step
    if ( e > emax ) {
      // check alternate possibility
      setPixel(x, y, ipt, newPt, eyer, farr);
      if (debug && validSquare(x,y)) boxes[x][y][0] = 1;

      x += xup;
      y += yup;
      e -= emax << 1; // error step
    }
    else { // check alternate possibility
      x+=xup;
      y+=yup;
      e+=einc;

      setPixel(x, y, ipt, newPt, eyer, farr);
      if (debug && validSquare(x,y)) boxes[x][y][0] = 1;

      x-=xup;
      y-=yup;
      e-=einc;
    }
  }
}


/* determines and computes intersection of ray eyer/farr with box x,y of heightfield */
void
RandomMeshscape::setPixel(int x, int y, float* ipt, float* newPt, float* eyer, float* farr) {
      if (debug_setPixel) cout << "found a square: " << x << " " << y << endl;
      if (validSquare(x,y)) {
	if (debug_setPixel) cout << "the square is valid" << endl;
	
	if (checkbox(x,y,newPt,eyer,farr )) {
	  if (debug_setPixel) cout << "checkbox confirms intersection" << endl;
	  _intersect = 1;
	  if (!_intersect) { // first intersection
	    if (debug_intersectWithRay) cout << "found first intersection" << endl;
	    assign_3vec(ipt, newPt);
	    }
	  }
	  else { // multiple intersection
	    if (debug_setPixel) cout << "found a multiple intersection" << endl;
	    if (first (newPt, ipt, eyer, farr)) {
	      assign_3vec(ipt, newPt);
	  }
	}
      }
}

// jjj-debug
// int nearestToRay ( Point3D eye, Point3D far, Point3D cpt )
//
//   returns 1 iff nearest point could be determined, 0 otherwise
//
int
RandomMeshscape::nearestToRay ( const Point3& eye, 
		       const Point3& far, 
		       Point3& npt ) {
  return 0;
}

// jjj-debug
//
// int nearestToPoint( Point3D p, Point3D cpt )
//
//   returns 1 iff nearest point could be determined, 0 otherwise
int
RandomMeshscape::nearestToPoint ( const Point3& p, Point3& cpt ) {
  return 1;
}

// jjj-debug
//
// int normalAtPoint( Point3D p, Point3D n )
//
//   returns 1 iff normal could be determined, 0 otherwise
int
RandomMeshscape::normalAtPoint ( const Point3&, Vector3& n ) {
  return 1;
}

//
// return axial bounding box of object
//
void
RandomMeshscape::bounds( Point3& min, Point3& max ) {
  
  Vector3 r ( 0.5 * _size,  0.5 * _size, 0.5 * _size );
  min  =  _pos - r;
  max  =  _pos + r;
  
}

/* draws triangle in red
   of box x,y
   tri1: bl,br,tl
   tri2: br,tl,tr */
void
RandomMeshscape::drawTriangle(int x, int y, int tri) {
    float lx,by,rx,ty;
    lx = getSquareXLeft(x);
    rx = getSquareXRight(x);
    by = getSquareYBottom(y);
    ty = getSquareYTop(y);

    if (tri==1) {
      glBegin( GL_TRIANGLES);
      glColor3f( 1.0F, 0.0F, 0.0F );
      glVertex3f(lx, by, _depth[x][y]);
      glVertex3f(rx, by, _depth[x+1][y]);
      glVertex3f(lx, ty, _depth[x][y+1]);
      glEnd();
    }
    if (tri==2) {
      glBegin( GL_TRIANGLES);
      glColor3f( 1.0F, 0.0F, 0.0F );
      glVertex3f(rx, by, _depth[x+1][y]);
      glVertex3f(lx, ty, _depth[x][y+1]);
      glVertex3f(rx, ty, _depth[x+1][y+1]);
      glEnd();
    }
}

/* draws outline of box x,y at z=0 */
void
RandomMeshscape::drawbox(int x, int y, int col, float z1, float z2){
    float lx,by,rx,ty;
    lx = getSquareXLeft(x);
    rx = getSquareXRight(x);
    by = getSquareYBottom(y);
    ty = getSquareYTop(y);
    glBegin(GL_LINES);
    glColor3fv( &_col[col] );                         // set generic color information
    glVertex3f(lx,by,z1);
    glVertex3f(lx,ty,z1);
    glEnd();
    glBegin(GL_LINES);
    glColor3fv( &_col[col] );                         // set generic color information
    glVertex3f(lx,ty,z1);
    glVertex3f(rx,ty,z1);
    glEnd();
    glBegin(GL_LINES);
    glColor3fv( &_col[col] );                         // set generic color information
    glVertex3f(rx,ty,z1);
    glVertex3f(rx,by,z1);
    glEnd();
    glBegin(GL_LINES);
    glColor3fv( &_col[col] );                         // set generic color information
    glVertex3f(rx,by,z1);
    glVertex3f(lx,by,z1);
    glEnd();

    glBegin(GL_LINES);
    glColor3fv( &_col[col] );                         // set generic color information
    glVertex3f(lx,by,z2);
    glVertex3f(lx,ty,z2);
    glEnd();
    glBegin(GL_LINES);
    glColor3fv( &_col[col] );                         // set generic color information
    glVertex3f(lx,ty,z2);
    glVertex3f(rx,ty,z2);
    glEnd();
    glBegin(GL_LINES);
    glColor3fv( &_col[col] );                         // set generic color information
    glVertex3f(rx,ty,z2);
    glVertex3f(rx,by,z2);
    glEnd();
    glBegin(GL_LINES);
    glColor3fv( &_col[col] );                         // set generic color information
    glVertex3f(rx,by,z2);
    glVertex3f(lx,by,z2);
    glEnd();


    glBegin(GL_LINES);
    glColor3fv( &_col[col] );                         // set generic color information
    glVertex3f(lx,by,z1);
    glVertex3f(lx,by,z2);
    glEnd();
    glBegin(GL_LINES);
    glColor3fv( &_col[col] );                         // set generic color information
    glVertex3f(lx,ty,z1);
    glVertex3f(lx,ty,z2);
    glEnd();
    glBegin(GL_LINES);
    glColor3fv( &_col[col] );                         // set generic color information
    glVertex3f(rx,ty,z1);
    glVertex3f(rx,ty,z2);
    glEnd();
    glBegin(GL_LINES);
    glColor3fv( &_col[col] );                         // set generic color information
    glVertex3f(rx,by,z1);
    glVertex3f(rx,by,z2);
    glEnd();

}

void       
RandomMeshscape::draw( void ) 
{
  if (debug) {
    // eye/far proj onto z=0
    glBegin(GL_LINES);
    glColor3fv( &_col[5] );                         // set generic color information
    glVertex3f(segx0,segy0,0.0);
    glVertex3f(segx1,segy1,0.0);
    glEnd();

    // eye/far in world coords
    glBegin(GL_LINES);
    glColor3fv( &_col[4] );                         // set generic color information
    glVertex3f(eyex,eyey,eyez);
    glVertex3f(farx,fary,farz);
    glEnd();

    // draw markers for projected eye/far boxes
    for (int xpos = 0; xpos < _nribs; xpos++) {
      for (int ypos = 0; ypos < _nribs; ypos++) {
	if (boxes[xpos][ypos][0]) drawbox(xpos,ypos,5, -0.05,.4);
	if (boxes[xpos][ypos][1]) drawbox(xpos,ypos,4, -0.1,.45);
	if (boxes[xpos][ypos][2]) drawbox(xpos,ypos,1, -0.15,.5);
      }
    }

    // draw intersected triangles

  } // if debug

  glPushMatrix();                                 // save
  glTranslatef( _pos[0], _pos[1], _pos[2] );      // move
  
  float xpos , ypos;
  xpos = 0.0F;

  for( int i=0; i<(_nribs-1) ; i++ ) {
    ypos = 0.0F;
    glBegin( GL_TRIANGLE_STRIP);                    // create triangle strips from depth matrix
    glColor3fv( &_col[0] );                         // set generic color information

    for( int j=0; j<_nribs ; j++ ) {
      glVertex3f( xpos,     ypos, _depth[i][j] );
      glVertex3f( xpos+_inc, ypos, _depth[i+1][j] );
      ypos += _inc;
    }
    glEnd();
    xpos += _inc;
  }
  glPopMatrix();                                 // restore
}

// draw bounding box around meshscape
void 
RandomMeshscape::drawBounds( void ){
  glPushMatrix();                         // save
  glTranslatef( _pos[0], _pos[1], _pos[2] );
  
  glLineWidth( 2 ) ;
  glColor3f( 1, 0, 0 );                  // red
  drawBox( -0.5f*_size , 1.05f*_size,          // x coords 
	   -0.5f*_size , 1.05f*_size,          // y coords 
	   -_size*0.5F, _size*1.05F,          // z coords
	   GL_LINE_LOOP);
  glLineWidth( 1 ) ;

  glPopMatrix();                         // restore
}

//
// set position to point implicated by ray ^ constraint
//
int 
RandomMeshscape::update( const Point3& eye, const Point3& far ){  
  
  if ( _constraint ) {
    Point3 intersectionPt;
    _constraint->nearestToRay( eye, far, intersectionPt );
    
    _pos = intersectionPt;    // set position of meshscape
    
    return 1;                 // object was modified; force redraw
  }
  else {
    cout << "RandomMeshscape::update(), nothing to constrain with" << endl;
    return 0;
  }
}

// given a pointer to four floats corresponding to each of the
// top left/right & bottom left/right corners of square (indexing starts
// at 0 and goes to nribs - 1) returns the heights of the square's corners
// at those points.
// BE CAREFUL--for speed, there is no bounds check on the xpos/ypos index's
void       
RandomMeshscape::getSquareHeights ( int xpos, int ypos, 
				    float * bleft , float * bright , 
				    float * tleft , float * tright ) {
  *bleft = _depth[xpos][ypos];
  *bright = _depth[xpos][ypos+1];
  *tleft = _depth[xpos+1][ypos];
  *tright = _depth[xpos+1][ypos+1];
}

// given a pointer to four floats corresponding to each of the
// top left/right & bottom left/right corners of square (indexing starts
// at 0 and goes to nribs - 1) returns the vertices of the square's corners
// at those points.
// BE CAREFUL--for speed, there is no bounds check on the xpos/ypos index's
void       
RandomMeshscape::getSquareVertices ( int xpos, int ypos, 
				    float * bleft , float * bright , 
				    float * tleft , float * tright ) {
  bleft[0] = getSquareXLeft(xpos);
  bleft[1] = getSquareYBottom(ypos);
  bleft[2] = _depth[xpos][ypos];

  bright[0] = getSquareXRight(xpos);
  bright[1] = getSquareYBottom(ypos);
  bright[2] = _depth[xpos][ypos+1];

  tleft[0] = getSquareXLeft(xpos);
  tleft[1] = getSquareYTop(ypos);
  tleft[2] = _depth[xpos+1][ypos];

  tright[0] = getSquareXRight(xpos);
  tright[1] = getSquareYTop(ypos);
  tright[2] = _depth[xpos+1][ypos+1];
}

float      
RandomMeshscape::getSquareYTop(int ypos) {
  float width = _size / _nribs;
  float yt = ((float)ypos + 1.0F) * width;
  return yt;
  //  return ((-_size/2.0F)+((ypos+1.0F)*_inc));
}
float      
RandomMeshscape::getSquareYBottom(int ypos) {
  float width = _size / _nribs;
  float yb = (float)ypos * width;
  return yb;
  //  return ((-_size/2.0F)+(ypos*_inc));
}
float     
RandomMeshscape::getSquareXLeft(int xpos) {
  float width = _size / _nribs;
  float xl = (float)xpos * width;
  return xl;
  //  return ((-_size/2.0F)+(xpos*_inc));
}
float    
RandomMeshscape::getSquareXRight(int xpos) {
  float width = _size / _nribs;
  float xr = ((float)xpos + 1.0F) * width;
  return xr;
  //  return ((-_size/2.0F)+((xpos+1.0F)*_inc));
}

//////////////////////////////////////////////////////////////////////////
//
// Sphere
//
//////////////////////////////////////////////////////////////////////////

Sphere::Sphere( float radius, int nlunes, int nbands ):
  _radius( radius ),
  _nlunes( nlunes ),
  _nbands( nbands ),
  _sphereQuadric ( NULL )
{
}

// int nearestToRay( Point3D eye, Point3D far, Point3D cpt )
//
//   returns 1 iff nearest point could be determined, 0 otherwise
int
Sphere::nearestToRay ( const Point3& eye, const Point3& far, Point3& cpt ) {
  
  // if ray intersects sphere, resulting point is nearest
  if ( intersectWithRay ( eye, far, cpt ) )
    return 1;

  // cout << "No direct intersection...checking nearest" << endl;
  /* find parametric eqn for eye-far ray */
  Point3 p1,p2;

  p2 = eye; 
  p1 = far - p2;

  /* minimize distance of point on eye-far ray to center of sphere
     minimize: D = (at + g - x)2 + (bt + h - y)2 + (ct + i - z)2
     so solve: 0 = 2a(at+g-x) + 2b(bt+h-y) + 2c(ct+i-z)
     assumes eye at infinity!!!! */
  float num,denom;
  
  num =  (p1 * (p2 - _pos)).sum();
  denom = powf( p1[0], 2 ) + powf( p1[1], 2 ) + powf( p1[2], 2 );
  float t = ( -1 * num ) / denom;
  
  Point3 startPt = _pos;
  Point3 endPt = ( p1 * t ) + p2;
  
  /* find where ray from center of sphere to computed nearest
     point on input ray intersects surface of sphere */
  if ( intersectWithRay( startPt, endPt, cpt ) )
    return 1;
  
  // otherwise, error
  cout << "error in Sphere::nearestToRay()" << endl;
  return 0;
}

// int nearestToPoint( Point3 p, Point3 cpt )
//
//   returns 1 iff nearest point could be determined, 0 otherwise
int
Sphere::nearestToPoint ( const Point3 &p, Point3& cpt ) {

  // if query point is sphere center, query is ill-defined
  
  
  if ( p == _pos ) 
    return 0;
  
  // otherwise, nearest point is along ray from center to query point
  Point3 ep;
  ep = _pos;

  if ( intersectWithRay( ep, p, cpt ) )
    return 1;

  // otherwise, error
  cout << "error in Sphere::nearestToPoint()" << endl;
  return 0;
}

//
// int normalAtPoint( Point3D p, Point3D n )
//
//   returns 1 iff normal could be determined, 0 otherwise
int
Sphere::normalAtPoint ( const Point3& p, Vector3& n ) {
  
  // if query point is sphere center, query is ill-defined
  if ( p == _pos )
    return 0;

  // otherwise, normal (gradient) is directed away from center
  n = p - _pos;

  // unitize it
  n.normalize();
  
  return 1;
}

int        
Sphere::intersectWithRay ( const Point3& eye, 
			   const Point3& far, 
			   Point3& ipt ) { 
  
  /* constrain addition to sphere
     line through e,f:
     x=at+g, y=bt+h, z=ct+i
     e will be at t=0, f at t=1 */
  /* also, translate e,f into sphere's local coordinate system */
  
  Point3 p1,p2;

  // translate eye, far to the object's position. 
  p2 = eye - _pos; 
  p1 = (far - _pos) - p2;

  //for coeffs of quadratic when plug parametric eqns into sphere constraint 
  float A,B,C;

  // Point3's pow method raises each element of the vector to the
  // scalar specifed.
  A = (p1*p1).sum();
  B = ( 2* (p1*p2) ).sum();
  C =  ( p2*p2 ).sum()  -  powf( _radius, 2 );
  
  
  // if t exists,
  // this will give 2 answers for t, so pick closest one to e
  //  and find corresponding x,y,z
  if (powf(B,2) >= (4*A*C)) {
    // chose smallest positive t
    float t = ( -B - sqrtf( powf( B, 2 ) - ( 4 * A * C ) ) ) / ( 2 * A );
    if ( t < 0 )
      t = ( -B + sqrtf( powf( B, 2 ) - ( 4 * A * C ) ) ) / ( 2* A );
    
    Point3 pt3 = ( p1 * t ) + p2;
    
    // return computed intersection point
    ipt = pt3;
    return 1;
  } else {
    // no real root, so no intersection
    return 0;
  }  
}

void
Sphere::bounds( Point3& min, Point3& max ) {
  
  Vector3 r ( _radius, _radius, _radius );
  min  =  _pos - r;
  max  =  _pos + r;
  
}

// draw sphere
void       
Sphere::draw( void ) {
  
  // save GL matrices
  glPushMatrix();
  
  // position sphere
  glTranslatef( _pos[0], _pos[1], _pos[2] );
  if ( !_sphereQuadric ) // lazy alloc
	_sphereQuadric =  gluNewQuadric();
  glColor3f( 1.0F,0.0F,0.0F );
  gluSphere( _sphereQuadric, _radius, _nlunes, _nbands );

  // restore
  glPopMatrix();

}

// draw bounding box around sphere
void 
Sphere::drawBounds( void ){
  float halfSide = ( _radius * 1.25F );

  // save
  glPushMatrix();
  glTranslatef( _pos[0], _pos[1], _pos[2] );
  
  glLineWidth( 2 ) ;

  glColor3f( 1, 0, 0 ); // red
  drawBox( -halfSide, halfSide,  // x coords 
		   -halfSide, halfSide,  // y coords
		   -halfSide, halfSide,  // z coords
		   GL_LINE_LOOP);

  glLineWidth( 1 ) ;

  // restore
  glPopMatrix();
}

int
Sphere::update( const Point3& eye, const Point3& far ){
  
  // if object is constrained, move sphere to closest point on base object
  if ( _constraint ) {
    // compute intersection point
    Point3 intersectionPt;
    _constraint->nearestToRay( eye, far, intersectionPt );
    
    // set position of sphere to computed intersection point
    _pos = intersectionPt;
    
    // modified, request redraw
    return 1;
  }
  else {
    cout << "Sphere::update(), nothing to constrain with" << endl;
    return 0;
  }
}

 
//////////////////////////////////////////////////////////////////////////
//
// XYRect
//
// A rectangle of dx == width, dy == height on plane of constant z.
//
//////////////////////////////////////////////////////////////////////////

XYRect::XYRect( float w, float h ) :
  _width( w ),
  _height( h )
{
}

int        
XYRect::intersectWithRay ( const Point3& eye, const Point3& far, 
			   Point3& ipt ) {
  Vector3 v;
  
  if( eye[2] != far[2] ) { 
    
    // find intersection point of ray with plane of constant Z
    float ratio = (eye[2] - _pos[2]) / (eye[2] - far[2]);
    
    v = eye + ( ratio * ( far - eye ) ) ;
    
    // case in which eye/far intersects XYRect plane, 
    // but is not within the bounds of the XYRect
    float halfWidth = _width / 2.0F;    
    float halfHeight = _height / 2.0F;  
    if(    ( v[0] < _pos[0] - halfWidth )
	   || ( v[0] > _pos[0] + halfWidth )
	   || ( v[1] < _pos[1] - halfHeight )
	   || ( v[1] > _pos[1] + halfHeight ) )
      return 0;

    // return computed intersection
    ipt = v;
    
    return 1;
  } else {
    cout << "Eye/Far parallel to z==0 plane!\n" << endl;
    return 0;
  }
}

// int nearestToRay ( Point3D eye, Point3D far, Point3D cpt )
//
//   returns 1 iff nearest point could be determined, 0 otherwise
//
int
XYRect::nearestToRay ( const Point3& eye, 
		       const Point3& far, 
		       Point3& npt ) {
  
  if( eye[2] != far[2] ) { 
    
    // find intersection point of ray with plane of constant Z
    float ratio = (eye[2] - _pos[2]) / (eye[2] - far[2]);
    npt = eye + (ratio * ( far - eye )) ;
    
    // case in which eye/far intersects XYRect plane, 
    // but is not within the bounds of the XYRect
    float halfWidth = _width / 2.0F;    
    float halfHeight = _height / 2.0F;  
    if ( npt[0] < _pos[0] - halfWidth ) npt[0] = _pos[0] - halfWidth;
    if ( npt[0] > _pos[0] + halfWidth ) npt[0] = _pos[0] + halfWidth;
    if ( npt[1] < _pos[1] - halfHeight ) npt[1] = _pos[1] - halfHeight;
    if ( npt[1] > _pos[1] + halfHeight ) npt[1] = _pos[1] + halfHeight;

    return 1;
  } else {
    cout << "Eye/Far parallel to z==0 plane!\n" << endl;
    return 0;
  }
}

//
// int nearestToPoint( Point3D p, Point3D cpt )
//
//   returns 1 iff nearest point could be determined, 0 otherwise
int
XYRect::nearestToPoint ( const Point3& p, Point3& cpt ) {
  
  // simply project perpendicular to plane
  cpt[0] = p[0];
  cpt[1] = p[1];
  cpt[2] = _pos[2];
  
  // clamp to dimensions of object
  float halfWidth = _width / 2.0F;    
  float halfHeight = _height / 2.0F;    
  if ( cpt[0] < -halfWidth )
    cpt[0] = -halfWidth;
  if ( cpt[1] < -halfHeight )
    cpt[1] = -halfHeight;
  if ( cpt[0] > halfWidth )
    cpt[0] = halfWidth;
  if ( cpt[1] > halfHeight )
    cpt[1] = halfHeight;

  return 1;
}

//
// int normalAtPoint( Point3D p, Point3D n )
//
//   returns 1 iff normal could be determined, 0 otherwise
int
XYRect::normalAtPoint ( const Point3&, Vector3& n ) {

  // normal is just zhat
  n = Point3( 0.F, 0.F, 1.F );
  return 1;
}

//
// return axial bounding box of object
//

void
XYRect::bounds( Point3& min, Point3& max ) {
  
  Vector3 r ( 0.5 * _width,  0.5 *_height, 0 );
  min  =  _pos - r;
  max  =  _pos + r;
  
}

void       
XYRect::draw( void ) 
{  
  // save
  glPushMatrix();  
  
  glTranslatef( _pos[0], _pos[1], _pos[2] );
  
  glBegin( GL_POLYGON );
  glColor3fv( &_col[0] );
  float halfWidth  = _width / 2.0F;
  float halfHeight = _height / 2.0F;
  glVertex2f( -halfWidth, halfHeight );
  glVertex2f( -halfWidth, -halfHeight );
  glVertex2f( halfWidth, -halfHeight ); 
  glVertex2f( halfWidth, halfHeight );
  glEnd();

  // restore
  glPopMatrix();  
}

void       
XYRect::drawBounds( void ) 
{  
  // save
  glPushMatrix();  

  glTranslatef( _pos[0], _pos[1], _pos[2] );
  
  glBegin( GL_LINE_LOOP );
  glColor3fv( &_col[0] );
  float halfWidth  = _width / 2.0F;
  float halfHeight = _height / 2.0F;
  glVertex2f( -halfWidth, halfHeight );
  glVertex2f( -halfWidth, -halfHeight );
  glVertex2f( halfWidth, -halfHeight ); 
  glVertex2f( halfWidth, halfHeight );
  glEnd();

  // restore
  glPopMatrix();  
}

//
// set position to point implicated by ray ^ constraint
//
int 
XYRect::update( const Point3& eye, const Point3& far ){  
  
  if ( _constraint ) {
    Point3 intersectionPt;
    _constraint->nearestToRay( eye, far, intersectionPt );
    
    // set position of rect
    _pos = intersectionPt;
    
    // object was modified; force redraw
    return 1;
  }
  else {
    cout << "XYRect::update(), nothing to constrain with" << endl;
    return 0;
  }
}