#include <ctype.h>
#include <math.h>
#include <string.h>
#include "GR_events.h"
#include "GR_BdryPatch2D.h"
#include "GR_Bdry2D.h"
#include "GR_List.h"

// *******************************************************************
//
// BdryArc
// -------
//
// Another descendant of BdryPatch2D.
// It implements features for circular arcs (and circles).
//
// *******************************************************************

BdryArc::
BdryArc () : iCenterPt(-1), dRadius(0), dAngleStart(0), dAngleEnd(0)
{
  iBCLeft = iBCRight = iInvalidBC;
  iRegionLeft = iRegionRight = iInvalidRegion;
  iBeginPt = iEndPt = -1;
  eGeom2DType = eCircularArc;
}

// Constructor for circles...
BdryArc::
BdryArc (Bdry2D* pBdryObject, const int iBCL, const int iBCR,
	 const int iRL, const int iRR,
	 const double dRad, const int iC):
  BdryPatch2D (pBdryObject, iBCL, iBCR, iRL, iRR, -1, -1),
  iCenterPt(iC), dRadius(dRad), dAngleStart(0.), dAngleEnd(2*M_PI)

{
  // Set boundary type
  eGeom2DType = eCircle;

  // Set data

  assert2(dRad > 0, "BdryCircle: Radius has to be greater than zero");
  assert2(iC >= 0, "BdryCircle: Center point index must be greater than zero");
  pB2D->vGetPoint(iC, adCenter);
  
  // Calculate information about segment...
  vRefresh ();
}

 	// constructor for arcs
BdryArc::
BdryArc(Bdry2D* pBdryObject, const int iBCL, const int iBCR,
	const int iRL, const int iRR,
	const int iBPt, const int iEPt, const double dRad,
	const bool qLongArc):
  BdryPatch2D(pBdryObject, iBCL, iBCR, iRL, iRR, iBPt, iEPt),
  iCenterPt(-1), dRadius(dRad), dAngleStart(0), dAngleEnd(0)
{
		
  // Set boundary type

  if (qLongArc) {
    eGeom2DType = eLongArc;
  }
  else {
    eGeom2DType = eCircularArc;
  }
	
  // Set the radius
  assert2(dRad > 0, "BdryArc: Radius has to be greater than zero");	
	
  // verts are given.. use them
  assert (iBPt != -1);
  assert (iEPt != -1);		
  
  double d1[2], d2[2];
  double dDX, dDY;
  double dDist;
  double dAngle, dOrient;
  //	char strInfo[150];
  
  vGetBeginPoint(d1);
  vGetEndPoint(d2);
  
  dDX = d2[0] - d1[0];
  dDY = d2[1] - d1[1];
  
  dOrient = atan2(dDY, dDX);
  
  // Get distance
  dDist = dDIST2D(d1, d2);			
  
  // Build info phrase..
  //	sprintf(strInfo, "Distance between (%f, %f) and (%f, %f) -- %f -- is too large for radius %f",
  //					d1[0], d1[1], d2[0], d2[1], dDist, dRadius);
  
  // must have this, or else impossible to have a circle..		
  assert(iFuzzyComp(dDist, 2.0*dRadius) <= 0/*, strInfo*/);
  
  // The arc is always CCW from pVB to pVE
  // Find angle between line joining two points and center
  double dacos = dDist/(2.0*dRadius);
  if (dacos > 1.0) dacos = 1.0;
  
  dAngle = acos(dacos);
  
  // Long arcs have the center on the other side of the chord.
  if (qLongArc) dAngle *= -1;
  
  // Set center
  adCenter[0] = d1[0] + dRadius*cos(dOrient + dAngle);
  adCenter[1] = d1[1] + dRadius*sin(dOrient + dAngle);
  
  // Find the angles..
  dAngleStart = atan2(d1[1] - adCenter[1], d1[0] - adCenter[0]);
  dAngleEnd   = atan2(d2[1] - adCenter[1], d2[0] - adCenter[0]);
  
  
  
  //	if (dOrient < 0) {
  //		dAngleStart = -dAngleStart;
  //		dAngleEnd   = -dAngleEnd;
  //	}
  
  // Make them positive
  if (dAngleStart < 0) dAngleStart = dAngleStart + 2*M_PI;
  if (dAngleEnd < 0) dAngleEnd = dAngleEnd + 2*M_PI;	
  
  if (dAngleStart >= 2*M_PI) dAngleStart -= 2*M_PI;
  
  while (dAngleStart > dAngleEnd) {
    dAngleEnd += 2*M_PI;
  }
//   while (dAngleStart < dAngleEnd - 2*M_PI) {
//     dAngleEnd += 2*M_PI;
//   }
  
  //	vMessage(2, "BdryArc: dOrient: %f, dAngle: %f, dAngleStart: %f, dAngleEnd: %f, adCenter: (%f, %f)\n",
  //							dOrient, dAngle, dAngleStart, dAngleEnd, adCenter[0], adCenter[1]);
  //	vMessage(2, "BdryArc: d1 (%f, %f), d2 (%f, %f), dDist: %f, dRadius: %f\n",
  //							d1[0], d1[1], d2[0], d2[1], dDist, dRadius);							
  
  // calculate stuff
  vRefresh();
}

void BdryArc::
vAssignData (Bdry2D *pBdryObject, const int iBCL, const int iBCR,
	     const int iRL, const int iRR,
	     const int iBPt, const int iEPt)
{
  BdryPatch2D::vAssignData (pBdryObject, iBCL, iBCR, iRL, iRR, iBPt, iEPt);
  // Compute data..
  vRefresh ();
}


void BdryArc::
vComputeMaxParameter ()
{
	char strInfo[100];
  // Compute the length of the segment.
  // For a circular arc, it's pretty easy...
  dMaxParam = (dAngleEnd - dAngleStart)*dRadius;

  // build info string
  sprintf(strInfo, "Angles (Start: %f, End: %f) not measured counterclockwise\n", dAngleStart, dAngleEnd);
  assert2(dMaxParam > 0, strInfo);
}

void BdryArc::
vPointAtParameter(const double dParamIn, double adLoc[2]) const
{
  // Returns the point located at arclength dParamIn, starting
  // from dAngleStart...
  // dParamIn IS ALLOWED to be greater than the arc's real arclength.
  
  // Arc know their pVBegin/pVEnd, so this might be
  // worth a shot...
  
  if ( ((eGeom2DType == eCircularArc) || (eGeom2DType == eLongArc)) &&
       ( iFuzzyComp(dParamIn, 0.0) == 0) ) {
    vGetBeginPoint(adLoc);
    return;
  }
		
  if ( ((eGeom2DType == eCircularArc) || (eGeom2DType == eLongArc)) &&
       (iFuzzyComp(dParamIn, dMaxParam) == 0) ) {
    vGetEndPoint(adLoc);
    return;
  }
				
  // Get the angle spanned by the arclength specified...
  double dAngle = dParamIn/dRadius;
  
  // The location is pretty easy to find now...
  adLoc[0] = dRadius*cos(dAngleStart+dAngle) + adCenter[0];
  adLoc[1] = dRadius*sin(dAngleStart+dAngle) + adCenter[1];
  // done..	
}

void BdryArc::
vPointAtParameterRatio (const double dRatio, double adLoc[2]) const
{
  // Return the point located at the ratio between pVBegin and pVEnd..
  // Ratio must be between 0 and 1 ...
  assert ((iFuzzyComp(dRatio, 0.0) != -1) &&
	  (iFuzzyComp(dRatio, 1.0) != 1));
  double dRealRatio = dRatio;
  if (dRatio < 0.0) dRealRatio = 0;
  if (dRatio > 1.0) dRealRatio = 1;

  // Get the coords
  vPointAtParameter(dRealRatio*dMaxParam, adLoc);
}

void BdryArc::
vPointAtDistFromStart(const double dDist, double adLoc[2]) const
{
	// Returns the point that is located a LINEAR distance dDist from
	// the start of the patch.
	
	// The angle spanned by such length is the following..
	double dAngle = 2.0*asin( dDist/(2.0*dRadius) );
	// once the angle is know, the rest is easy...
	vPointAtParameter(dAngle*dRadius, adLoc);
}

void BdryArc::
vPointAtDistFromEnd(const double dDist, double adLoc[2]) const
{
	// The angle spanned by such length is the following..
	double dAngle = 2.0*asin( dDist/(2.0*dRadius) );
	
	// once the angle is know, the rest is easy...
	vPointAtParameter(dMaxParam - (dAngle*dRadius), adLoc);

}
double BdryArc::
dCurvature (const double[2]) const
{
  // Calculate the curvature of a circular arc
  // Actually, simple since the curvature is the inverse
  // of the radius, and we know the radius...

  return (1.0/dRadius);
}

double BdryArc::
dAngleClosestToPoint ( const double adLoc[2]) const {
	// Does the same as dParameterClosestToPoint, but returns the angle,
	// instead of the arclength...
	// The angle is relative to dAngleStart

  double adPtBegin[2], adPt90[2];
  double adNormBegin[2], adNorm90[2], adNorm[2];
  double dAngle, dAngle90;
  double dDot, dDot90;

	// Get the point at the beginning
	vPointAtParameter(0.0, adPtBegin);
	// Get the point at 90 degrees
	vPointAtParameter(dRadius*M_PI/2.0, adPt90);
		
	// Normal at beginning
	vUnitNormal(adPtBegin, adNormBegin);
	// Normal at 90 degrees
	vUnitNormal(adPt90, adNorm90);
	// Normal at point given by user
	vUnitNormal(adLoc, adNorm);
	
	// Calculate dot products
	dDot = dDOT2D(adNormBegin, adNorm);
	dDot90 = dDOT2D(adNorm90, adNorm);
	
	if (dDot > 1.0) dDot = 1.0;
	if (dDot < -1.0) dDot = -1.0;
	if (dDot90 > 1.0) dDot90 = 1.0;
	if (dDot90 < -1.0) dDot90 = -1.0;
	
	// Calculate the angles spanned by these dotproducts..
	dAngle = acos(dDot);
	dAngle90 = acos(dDot90);
	
	// Use the second dot product to determine real angle..!
	if (iFuzzyComp(dAngle90, M_PI/2.0) == 1) {
		// Angle was definitely greater than 90...
		dAngle = (2.0*M_PI) - dAngle;
	}
	else {
		// Angle was less than 90...
		// No change, I guess..
	}
		
	if (dAngle > 2.0*M_PI) {
		dAngle = dAngle - 2.0*M_PI;
	}
	
	// Angle should be positive, right..?
	assert(dAngle >= 0.0);

	return dAngle;		
	
}

double BdryArc::
dParameterClosestToPoint (const double adLoc[2]) const
{
  // Returns the arclength from pVBegin to the point on the arc
  // closest to the point given by adLoc[2]. The arclength might be negative
  // or greater than dParam, depending on the location of adLoc,
  // i.e. the location given by the arclength calculated in here is not
  // restricted to be between pVBegin and pVEnd.

  double dAngle = dAngleClosestToPoint(adLoc);
	
  // Arclength is easily calculated now..
  return dRadius*dAngle;
}

void BdryArc::
vBdryPointClosestToPoint (const double adLoc[2], double adPoint[2]) const
{
  // Returns the point ON THE CIRCULAR ARC closest to the point
  // given in adLoc[2].
  double dAngle, dAngleE, dAngleS;
  double dMagS, dMagE;

  dAngle = dAngleClosestToPoint (adLoc);
  // Must make this one absolute..
  dAngle += dAngleStart;

  // is it in between?

  if ( (dAngle > dAngleStart) && (dAngle < dAngleEnd) ) {
		// return the point right there...
		vPointAtParameter((dAngle-dAngleStart)*dRadius, adPoint);
  }
  else {

  	// Find out if it's closer to start or end..
  	dAngleE = dAngleEnd;
  	dAngleS = dAngleStart;

  	// Find differences...  	
  	dMagS = fabs(dAngle - dAngleS);	  	
  	dMagE = fabs(dAngle - dAngleE);
	  	
  	// Restrict results to shortest distance on circle..
	  if (dMagS > M_PI) dMagS = 2.0*M_PI - dMagS;
	  if (dMagE > M_PI) dMagE = 2.0*M_PI - dMagE;
	
	  if (dMagS < dMagE) {
			// it's closer to start...
			vPointAtParameter(0, adPoint);
	  }
	  else {
			// it's closer to end...
			vPointAtParameter(dMaxParam, adPoint);
	  }
	
	}

}

double BdryArc::
dDistFromBdry(const double adLoc[2], int * const piSide) const
{
  double adPoint[2], adNormal[2], dX, dY, dDot;

  // Find closest point location..
  vBdryPointClosestToPoint (adLoc, adPoint);
  // Get the normal there..
  // Points outward except if iRegionSide() returns REGION_SIDE_LEFT
  vUnitNormal(adPoint, adNormal);

  // Calculate distance..
  dX = adPoint[0] - adLoc[0];
  dY = adPoint[1] - adLoc[1];

  // Calculate if point is inside (where normal points) or outside
  // Use Dot product
  dDot = dX * adNormal[0] + dY * adNormal[1];

  double dDist = sqrt( dX*dX + dY*dY );

  // -ve dot means in direction of normal
  if (piSide) {
    if ((dDot <= 0.0)) {
      if (iRegionWhichSide() == REGION_SIDE_LEFT) {
	*piSide = SIDE_LEFT;
      } else {
	*piSide = SIDE_RIGHT;
      }
    }
    else {
      if (iRegionWhichSide() == REGION_SIDE_LEFT) {
	*piSide = SIDE_RIGHT;
      } else {
	*piSide = SIDE_LEFT;
      }
    }
  }
  return dDist;
}


double BdryArc::
dParameterAtPoint (const double adLoc[2]) const
{
  // Returns the arclength for the point...
  // This will actually return the arclength for the point
  // on the boundary closest to the point given.

  return dParameterClosestToPoint (adLoc);
}

void BdryArc::
vNormalAtParameterRatio(const double dRatio, double adNorm[2]) const
{

	// Find the point first
	double adPoint[2];
	
	vPointAtParameterRatio(dRatio, adPoint);
	vUnitNormal(adPoint, adNorm);

}

void BdryArc::
vUnitNormal (const double adLoc[], double adNorm[]) const
{
  // Returns the normal at location adLoc[].

  // adLoc is assumed to be on the arc/circle. If not, the projection
  // is used (i.e. no check is made to make sure it is on the arc)

  double dMag;

  // Calculate the "outward" normal

  adNorm[0] = adLoc[0] - adCenter[0];
  adNorm[1] = adLoc[1] - adCenter[1];
  dMag = sqrt(adNorm[0]*adNorm[0] + adNorm[1]*adNorm[1]);

  // Could be a problem if the CENTER was given as the location...
  // If that's the case, then point towards the beginning of the arc
  if (iFuzzyComp(dMag, 0.0) == 0) {
		// Magnitude is VERY small, it is very close, so point towards the
		// beginning
		adNorm[0] = cos(dAngleStart);
		adNorm[1] = sin(dAngleStart);
  }
  else {
	  // Make it unit normal

  	adNorm[0] /= dMag;
	  adNorm[1] /= dMag;	
	
	}

  // Now, make sure orientation is correct
	// Change is needed only if region is to the left, i.e. "inside"
	
  if (iRegionSide == REGION_SIDE_LEFT) {
		// Revert the normal...
		adNorm[0] = -adNorm[0];
		adNorm[1] = -adNorm[1];
  }

  // Done..
}

bool BdryArc::
qIsColinear (BdryPatch2D *, Vert *, int *)
{
  assert2(0, "qIsColinear not yet defined for arcs\n");
  return(false);
}

List<double> BdryArc::
LdOriginalDiscretization(bool *qInsertVertPointers, bool *qCloseLoop)
{
  // This returns the original discretization info
  // For a BdryArc, it is a little more complicated...
  // If it's an arc, tell it to insert vert pointers,
  // and to not close the loop
  // If it's a circle, there are no vert pointers,
  // and it must close the loop.
  
  double adPt[2];
  List<double> LdList;
	
  if (eGeom2DType == eCircularArc || eGeom2DType == eLongArc) {
    // it's an arc
    
    // Determine # of insertion points..
    int iNumVerts = static_cast<int>(floor(dMaxParam/dRadius));
    
    if (iNumVerts > 0) {
      int k;
      int iNumSegs = iNumVerts + 1;
      double dSteps = dMaxParam / iNumSegs;
      
      for (k = 1; k < iNumSegs; k++) {
	vPointAtParameter(dSteps*k, adPt);
	LdList.vAppendItem(adPt[0]);
	LdList.vAppendItem(adPt[1]);
      }
    } else {
      // no other things to insert..
    }
    
    *qInsertVertPointers = true;
    *qCloseLoop = false;
  }
  else {
    // it must be a circle
    assert(eGeom2DType == eCircle);
    
    // Determine the insertion points
    // For a circle, it's easy. It is split in six sections,
    // and we need six verts
    
    int kmax = 6;
    int k;
    
    for (k = 0; k < kmax; k++) {
      double dk = k;
      double dkmax = kmax;
      double dRatio = dk/dkmax;
      // Get point
      vPointAtParameterRatio(dRatio, adPt);
      // add to list
      LdList.vAppendItem(adPt[0]);
      LdList.vAppendItem(adPt[1]);
    }
    
    *qInsertVertPointers = false;
    *qCloseLoop = true;
  }
  
  return LdList;
}

void BdryArc::
vRatioBetweenTwoPoints(const double dRatio, const double adLoc1[2],
		       const double adLoc2[2], double adPoint[2])  const
{
	// Return the location of the point on the boundary that is located
	// at ratio dRatio between adLoc1[] and adLoc2[].
	
	// Both these points "should" be on the boundary...	
//	assert(iFuzzyComp(dDistFromBdry(adLoc1), 0.0) == 0);
//	assert(iFuzzyComp(dDistFromBdry(adLoc2), 0.0) == 0);
	
	double dAL1, dAL2, dALRatio;
	
  // Get first arclength
  dAL1 = dParameterAtPoint(adLoc1);
  // Get the second one
  dAL2 = dParameterAtPoint(adLoc2);
	
  if (eGeom2DType == BdryPatch2D::eCircle) {
  	// Special case for circle -- there are 2 midpoints!
  	bool qOneGreater = (dAL1 > dAL2);
		double dHalf = M_PI*dRadius;
		
		double dDiff = fabs(dAL1 - dAL2);
		
		if (dDiff > dHalf) {
			dDiff = 2.0*dHalf - dDiff;
			qOneGreater = !qOneGreater;
		}
		
		dDiff = dDiff*dRatio;
		
		// Choose the "shortest" midpoint		
		if (qOneGreater) {
			dALRatio = dAL1 - dDiff;
		} else {
			dALRatio = dAL1 + dDiff;
		}
		
		if (dALRatio < 0.0) dALRatio = 2*dHalf + dALRatio;
				
	}	else {
		// it's an arc..
		// Get the arclength we want
		
		// Check a few things
		if (iFuzzyComp(dAL1, dMaxParam) > 0) dAL1 -= 2.0*M_PI*dRadius;
		if (iFuzzyComp(dAL2, dMaxParam) > 0) dAL2 -= 2.0*M_PI*dRadius;
		
		dALRatio = dAL1 + dRatio*(dAL2 - dAL1);		
	}
	
	// Return the point
	vPointAtParameter(dALRatio, adPoint);	  										  										
}

void BdryArc::vGeodesicDistanceRatioBetweenTwoPoints(double dRatio, const double adLoc1[2], 
						     const double adLoc2[2], double adPoint[2]) const {
  vRatioBetweenTwoPoints(dRatio, adLoc1, adLoc2, adPoint);
}


double BdryArc::
dXMax() {
	// Returns the maximum X coordinate of the patch..
	
	// Maximum coordinate will be when angle is closer to 0 or 2Pi
	
	double adLoc[2];
	double adPt[2];
	
	adLoc[0] = adCenter[0] + dRadius;
	adLoc[1] = adCenter[1];
	
	vBdryPointClosestToPoint(adLoc, adPt);
	return adPt[0];
}

double BdryArc::
dYMax() {

	// Maximum coordinate will be when angle is closer to Pi/2
	
	double adLoc[2];
	double adPt[2];
	
	adLoc[0] = adCenter[0];
	adLoc[1] = adCenter[1] + dRadius;
	
	vBdryPointClosestToPoint(adLoc, adPt);
	return adPt[1];

}

double BdryArc::
dXMin() {
	// Minimum coordinate will be when angle is closer to Pi
	
	double adLoc[2];
	double adPt[2];
	
	adLoc[0] = adCenter[0] - dRadius;
	adLoc[1] = adCenter[1];
	
	vBdryPointClosestToPoint(adLoc, adPt);
	return adPt[0];

}

double BdryArc::
dYMin() {
	// Returns the minimum Y coordinate of the patch..
	// Minimum coordinate will be when angle is closer to 3Pi/2
	
	double adLoc[2];
	double adPt[2];
	
	adLoc[0] = adCenter[0];
	adLoc[1] = adCenter[1] - dRadius;
	
	vBdryPointClosestToPoint(adLoc, adPt);
	return adPt[1];

}

void BdryArc::
vRefresh ()
{
  // Make sure the pre-calculated items are up-to-date...

  // Make sure the length is computed properly...
  vComputeMaxParameter ();
}

BdryArc::~BdryArc ()
{
  // Destructor
}

