#include <math.h>
#include "GR_misc.h"
#include "GR_Geometry.h"
#include "GR_Vertex.h"
#include "GR_Face.h"
#include "GR_Cell.h"
#include "GR_AdaptPred.h"
#include "GR_Vec.h"

//@ Compute the vector normal for a set of three verts
void
vNormal2D (const double adA[2], const double adB[2],
	   const double adC[2], double *const pdRes)
{
  double adTempB[] =
  {adB[0] - adA[0],
   adB[1] - adA[1]};
  double adTempC[] =
  {adC[0] - adA[0],
   adC[1] - adA[1]};
  *pdRes = adTempB[0] * adTempC[1] - adTempB[1] * adTempC[0];
}

//@ Compute the vector normal for a set of two verts
// This is always rotated 90 degrees counterclockwise from the vector
// running from vertex A to vertex B.
void
vNormal (const double adA[2], const double adB[2], double adRes[2])
{
  adRes[0] = adA[1] - adB[1];
  adRes[1] = adB[0] - adA[0];
}

//@ Determine whether a vertex lies inside, on or outside a triangle
int
iIsInTri (const Vert * const pV, const Vert * const pVertA,
	  const Vert * const pVertB, const Vert * const pVertC)
  // For pV to be inside the tri formed by pVertA-C, the orientation
  // has to be right wrt every face.  The return value is -1 if the
  // vert lies outside the tri, 0 if it's colinear with all three
  // faces (impossible for valid tris), 1 if it coincides with a vert
  // in the tri, 2 if it lies on an edge, and 3 if it lies strictly
  // inside. 
{
  assert (pV->qValid () && pV->iSpaceDimen () == 2);
  assert (pVertA->qValid () && pVertA->iSpaceDimen () == 2);
  assert (pVertB->qValid () && pVertB->iSpaceDimen () == 2);
  assert (pVertC->qValid () && pVertC->iSpaceDimen () == 2);
  int iRetVal = 0;
  int iCheck = iOrient2D (pVertA, pVertB, pV);
  if (iCheck == -1)
    return -1;
  iRetVal += iCheck;

  iCheck = iOrient2D (pVertB, pVertC, pV);
  if (iCheck == -1)
    return -1;
  iRetVal += iCheck;

  iCheck = iOrient2D (pVertC, pVertA, pV);
  if (iCheck == -1)
    return -1;
  iRetVal += iCheck;

  return (iRetVal);
}


//@ 2D orientation test
int
iOrient2D (const Vert * const pVertA, const Vert * const pVertB,
	   const Vert * const pVertC)
  // Computes the orientation of three verts in 2D
{
  assert (pVertA->qValid ());
  assert (pVertB->qValid ());
  assert (pVertC->qValid ());
  assert (pVertA->iSpaceDimen () == 2);
  assert (pVertB->iSpaceDimen () == 2);
  assert (pVertC->iSpaceDimen () == 2);
  return iOrient2D (pVertA->adCoords (), pVertB->adCoords (),
		    pVertC->adCoords ());
}

int
iOrient2D (const double adA[2], const double adB[2],
	   const double adC[2])
{
  if (qAdaptPred) {
    double dResult = orient2d_shew (adA, adB, adC);
    if (dResult > 0)
      return 1;
    else if (dResult < 0)
      return -1;
    else
      return 0;
  }
  else {
    double xa = adA[0];
    double ya = adA[1];

    double xb = adB[0];
    double yb = adB[1];

    double xc = adC[0];
    double yc = adC[1];

    // Currently no exact arithmetic
    double dDet = ((xb - xa) * (yc - ya) - (xc - xa) * (yb - ya));
    // Scale the determinant by the mean of the magnitudes of vectors:
    double dScale = (dDIST2D (adA, adB) + dDIST2D (adB, adC)
		     + dDIST2D (adC, adA)) / 3;
    dDet /= (dScale * dScale);
    if (dDet > 1.e-10)
      return 1;
    else if (dDet < -1.e-10)
      return -1;
    else
      return 0;
  }
}

//@ Incircle test
int
iIncircle (const Vert * const pVertA, const Vert * const pVertB,
	   const Vert * const pVertC, const Vert * const pVertD)
  // Determines whether pVertD is inside or outside the circumcircle
  // of the tri formed by pVertA-C, using as nearly exact arithmetic
  // as required.  Returns 1 if pVertD is inside the
  // circumcircle of the tri, -1 if outside, 0 if on. 
{
  assert (pVertA->qValid ());
  assert (pVertB->qValid ());
  assert (pVertC->qValid ());
  assert (pVertD->qValid ());
  return iIncircle(pVertA->adCoords(), pVertB->adCoords(),
		   pVertC->adCoords(), pVertD->adCoords());
}

int
iIncircle (const double adA[2], const double adB[2],
	   const double adC[2], const double adD[2])
  // Determines whether pVertD is inside or outside the circumcircle
  // of the tri formed by pVertA-C, using as nearly exact arithmetic
  // as required.  Returns 1 if pVertD is inside the
  // circumcircle of the tri, -1 if outside, 0 if on. 
{
  if (qAdaptPred) {
    double dResult = incircle_shew (adA, adB, adC, adD);
    if (dResult > 0)
      return 1;
    else if (dResult < 0)
      return -1;
    else
      return 0;
  }
  else {
    double xa = adA[0];
    double ya = adA[1];
    double wa2 = xa * xa + ya * ya;

    double xb = adB[0];
    double yb = adB[1];
    double wb2 = xb * xb + yb * yb;

    double xc = adC[0];
    double yc = adC[1];
    double wc2 = xc * xc + yc * yc;

    double xd = adD[0];
    double yd = adD[1];
    double wd2 = xd * xd + yd * yd;

    // Currently no exact arithmetic
    double dDet1 = ((xb - xa) * ((yc - ya) * (wd2 - wa2) - (yd - ya) * (wc2 - wa2))
		    + (xc - xa) * ((yd - ya) * (wb2 - wa2) - (yb - ya) * (wd2 - wa2))
		    + (xd - xa) * ((yb - ya) * (wc2 - wa2) - (yc - ya) * (wb2 - wa2))
		    );
    int iOrient = iOrient2D(adA, adB, adC);
    double dProd = dDet1 * iOrient;
    // Scale by the mean of the magnitudes of vectors:
    double dScale = (dDIST2D (adB, adA) +
		     dDIST2D (adC, adA) +
		     dDIST2D (adB, adC)) / 3;
    // Unclear exactly how to scale this...
    dProd /= (dScale * dScale * dScale * dScale);
    if (dProd > 1.e-10)
      return -1;
    else if (dProd < -1.e-10)
      return 1;
    else
      return 0;
  }
}

//@ End of file
//@ End of file
