/*-------------------------------------------------------------------------
   Create a contour slice through a 3 vertex facet "p"
   Given the normal of the cutting plane "n" and a point on the plane "p0"
   Return
       0 if the contour plane doesn't cut the facet
       2 if it does cut the facet, the contour line segment is p1->p2
      -1 for an unexpected occurence
   If a vertex touches the contour plane nothing need to be drawn!?
   Note: the following has been written as a "stand along" piece of
   code that will work but is far from efficient....
*/
int ContourFacet(XYZ *p,XYZ n,XYZ p0,XYZ *p1,XYZ *p2)
{
   double A,B,C,D;
   double l;
   double side[3];

   /*
      Determine the equation of the plane as
      Ax + By + Cz + D = 0
   */
   l = sqrt(n.x*n.x + n.y*n.y + n.z*n.z);
   A = n.x / l;
   B = n.y / l;
   C = n.z / l;
   D = -(n.x*p0.x + n.y*p0.y + n.z*p0.z);

   /*
      Evaluate the equation of the plane for each vertex
      If side < 0 then it is on the side to be retained
      else it is to be clipped
   */
   side[0] = A*p[0].x + B*p[0].y + C*p[0].z + D;
   side[1] = A*p[1].x + B*p[1].y + C*p[1].z + D;
   side[2] = A*p[2].x + B*p[2].y + C*p[2].z + D;

   /* Are all the vertices on the same side */
   if (side[0] >= 0 && side[1] >= 0 && side[2] >= 0)
      return(0);
   if (side[0] <= 0 && side[1] <= 0 && side[2] <= 0)
      return(0);

   /* Is p0 the only point on a side by itself */
   if ((SIGN(side[0]) != SIGN(side[1])) && (SIGN(side[0]) != SIGN(side[2]))) {
      p1->x = p[0].x - side[0] * (p[2].x - p[0].x) / (side[2] - side[0]);
      p1->y = p[0].y - side[0] * (p[2].y - p[0].y) / (side[2] - side[0]);
      p1->z = p[0].z - side[0] * (p[2].z - p[0].z) / (side[2] - side[0]);
      p2->x = p[0].x - side[0] * (p[1].x - p[0].x) / (side[1] - side[0]);
      p2->y = p[0].y - side[0] * (p[1].y - p[0].y) / (side[1] - side[0]);
      p2->z = p[0].z - side[0] * (p[1].z - p[0].z) / (side[1] - side[0]);
      return(2);
   }

   /* Is p1 the only point on a side by itself */
   if ((SIGN(side[1]) != SIGN(side[0])) && (SIGN(side[1]) != SIGN(side[2]))) {
      p1->x = p[1].x - side[1] * (p[2].x - p[1].x) / (side[2] - side[1]);
      p1->y = p[1].y - side[1] * (p[2].y - p[1].y) / (side[2] - side[1]);
      p1->z = p[1].z - side[1] * (p[2].z - p[1].z) / (side[2] - side[1]);
      p2->x = p[1].x - side[1] * (p[0].x - p[1].x) / (side[0] - side[1]);
      p2->y = p[1].y - side[1] * (p[0].y - p[1].y) / (side[0] - side[1]);
      p2->z = p[1].z - side[1] * (p[0].z - p[1].z) / (side[0] - side[1]);
      return(2);
   }

   /* Is p2 the only point on a side by itself */
   if ((SIGN(side[2]) != SIGN(side[0])) && (SIGN(side[2]) != SIGN(side[1]))) {
      p1->x = p[2].x - side[2] * (p[0].x - p[2].x) / (side[0] - side[2]);
      p1->y = p[2].y - side[2] * (p[0].y - p[2].y) / (side[0] - side[2]);
      p1->z = p[2].z - side[2] * (p[0].z - p[2].z) / (side[0] - side[2]);
      p2->x = p[2].x - side[2] * (p[1].x - p[2].x) / (side[1] - side[2]);
      p2->y = p[2].y - side[2] * (p[1].y - p[2].y) / (side[1] - side[2]);
      p2->z = p[2].z - side[2] * (p[1].z - p[2].z) / (side[1] - side[2]);
      return(2);
   }

   /* Shouldn't get here */
   return(-1);
}