// system includes
#include <math.h>
#include <iostream.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!

RandomMeshscape::RandomMeshscape( float s, int r ) :
  _size( s ),
  _nribs( r )
{
  cout << "Heiner thinks he's the shit" << endl;
  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/(flog10(rand()));
      dist = 1.0f / hypotf((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;
    }
  }
  cout << "Marc owes Heiner $$$" << endl;
  _squaremins = new float *[ r ];
  _squaremaxs = new float *[ r ];
  float *a,*b,*c,*d;
  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;
    }
  }*/
  cout << "NOT!!" << endl;
}


/* returns 1 if pt1 is before or at pt2 along ray from eye->far */
int
RandomMeshscape::first(float* pt1, float* pt2, const Point3& e, const Point3& f) {
  // we know pt1,pt2 are on the ray, so just check with x values
  float a,b,t1,t2;
  b = e[0];
  a = f[0] - e[0];
  t1 = (pt1[0] - b) / a;
  t2 = (pt2[0] - b) / a;
  if (t1<=t2) return 1;
  else return 0;
}

/* 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 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_of(float* pt1, float* pt2, float* pt3) {
  // line = pt1 x pt2
  float cross[3];
  cross_product(pt1,pt2,cross);
  // if p.l>=0, is leftof or on
  float pdotl;
  pdotl = cross[0]*pt3[0] + cross[1]*pt3[1] + cross[2]*pt3[2];
  if (pdotl>=0) {
    cout << "left_of found" << endl;
    return 1;
  }
  else return 0;
}

/* if intersectPt is inside triangle pt1,pt2,pt3, returns 1
   else returns 0 */
int
RandomMeshscape::inside_triangle(float* pt1,float* pt2,float* pt3,float* intersectPt) {
  if ( left_of(pt1,pt2,intersectPt) &&
       left_of(pt2,pt3,intersectPt) &&
       left_of(pt3,pt1,intersectPt) )
    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 = pt1xpt2
  float line1[3];
  cross_product(pt1,pt2,line1);
  // line2 = pt1xpt3
  float line2[3];
  cross_product(pt1,pt3,line2);
  // normal A,B,C to plane = line1xline2
  float normal[3];
  cross_product(line1,line2,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] );
  // 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;
  if (t<0) { // pt of intersect before beginning of ray...should never happen!
    cout << "you should never read this!  i'm inside 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]);
  return 1;
}

/* 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)) {
    cout << "inside plane" << endl;
    if (inside_triangle(pt1,pt2,pt3,intersectPt)) {
      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
    cout << "within height bounds" << endl;
    getSquareVertices(x,y,bleft,bright,tleft,tright);
    if (intersect_triangle(bleft, bright, tleft, pt1, eye, far)) {
      cout << "intersects triangle 1" << endl;
      tri1 = 1;
      };
    if (intersect_triangle(bright, tright, tleft, pt2, eye, far)) {
      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) intersectPt = pt1;
  if (tri2 && !tri1) intersectPt = pt2;
  if (tri1 && tri2) { // find the first pt along the ray and return it
    if (first(pt1,pt2,eye,far)) intersectPt = pt1;
    else intersectPt = pt2;
  }
  return 1;
}

/* returns 1 if eye->far ray is within bounding box of square x,y
   else returns 0 */
int
RandomMeshscape::checkheight(int x, int y, const Point3& eye, const Point3& far) {
  //  _squaremins[i][j];
  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] = ceil(y/_nribs);
  y = m*_size + b;
  square2[0] = ceil(_size/_nribs);
  square2[1] = ceil(y/_nribs);
}

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



// jjj-debug
int        
RandomMeshscape::intersectWithRay ( const Point3& eye, const Point3& far, 
			   Point3& ipt ) {
  int intersect = 0;
  float intersectPt[3];

  float A,B,C,m,b;
  A = eye[1] - far[1];
  B = far[0] - eye[0];
  C = eye[0]*far[1] - eye[1]*far[0];
  if (B==0) { // vertical line
    cout << "fucking vertical line, mother fucker" << endl;
  }
  else { // do bresenham
    cout << "happy joy, non-vertical fuck" << endl;
    float m,b;
    m = -A/B;
    b = -C/B;
  }
  
  int square1[2];
  int square2[2];
  getBoundingSquares(square1,square2,m,b);
  int x1 = square1[0];
  int y1 = square1[1];
  int x2 = square2[0];
  int y2 = square2[1];
  float newPt[3]; // new intersection point
  
  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 ) {
    cout << "found a square" << endl;
    if (validSquare(x,y)) {
      cout << "the square is valid" << endl;
      if (checkbox(x,y,newPt,eye,far )) {
	cout << "checkbox confirms intersection" << endl;
	if (intersect)
	  cout << "found first intersection" << endl;
	  if (first (newPt, intersectPt, eye, far)) {
	    ipt[0] = newPt[0];
	    ipt[1] = newPt[1];
	    ipt[2] = newPt[2];
	  }
	else {
	  cout << "found a multiple intersection" << endl;
	  intersectPt[0] = newPt[0];
	  intersectPt[1] = newPt[1];
	  intersectPt[2] = newPt[2];
	  ipt[0] = newPt[0];
	  ipt[1] = newPt[1];
	  ipt[2] = newPt[2];
	}
	intersect = 1;
      }
    }
    x += xinc;
    y += yinc; 
    e += einc; // advance step
    if ( e > emax ) {
      x += xup;
      y += yup;
      e -= emax << 1; // error step
    }
  }
  
  if (intersect) {
    cout << "ray intersects" << endl;
    return 1;
  }
  cout << "ray does not intersect" << endl;
  return 0;
}

// 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;
  
}
/*  The old way....changed to make our relationship to the globabl coordinate system easier to deal with...
void       
RandomMeshscape::draw( void ) 
{  
  glPushMatrix();                                 // save
  glTranslatef( _pos[0], _pos[1], _pos[2] );      // move
  
  float xpos , ypos;
  xpos = - _size / 2.0F;

  for( int i=0; i<(_nribs-1) ; i++ ) {
    ypos = - _size / 2.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 ){
  float halfSide = ( _size * 0.55F );

  glPushMatrix();                         // save
  glTranslatef( _pos[0], _pos[1], _pos[2] );
  
  glLineWidth( 2 ) ;
  glColor3f( 1, 0, 0 );                  // red
  drawBox( -halfSide, halfSide,          // x coords 
		   -halfSide, halfSide,  // y coords
		   -halfSide+(_size*0.5F), halfSide+(_size*0.5F),  // z coords
		   GL_LINE_LOOP);
  glLineWidth( 1 ) ;

  glPopMatrix();                         // restore
}

*/
void       
RandomMeshscape::draw( void ) 
{  
  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) {
  return ((-_size/2.0F)+((ypos+1.0F)*_inc));
}
float      
RandomMeshscape::getSquareYBottom(int ypos) {
  return ((-_size/2.0F)+(ypos*_inc));
}
float     
RandomMeshscape::getSquareXLeft(int xpos) {
  return ((-_size/2.0F)+(xpos*_inc));
}
float    
RandomMeshscape::getSquareXRight(int xpos) {
  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();
  glColor3fv( &_col[0] );
  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;
  }
}