#include <math.h>
#include "GR_BdryPatch2D.h"

// *******************************************************************
//
// BdryInterpolated2D
// -------
//
// Descendant of BdryPatch2D. This boundary patch is very different than
// other ones in the sense that it is not a patch itself, but rather a collection
// of other patches (cubic parametric ones in this case).
//
// *******************************************************************

BdryInterpolated2D::
BdryInterpolated2D(Bdry2D* pBdryObject, const int iBCL, const int iBCR,
		   const int iRL, const int iRR,
		   const List<int> LiThePoints) :
  BdryPatch2D(pBdryObject, iBCL, iBCR, iRL, iRR, LiThePoints[0],
	      LiThePoints[LiThePoints.iLength()-1]),
  LpBCP2D(), LiPoints(), LdOrderedRoots(),
  LdOrientationChangeAtRoots(), dOrientationChange(0),
  qUseLastCoeffs(false), iLastNewPoint(-1)
{
  // Constructor
  // Must have at least three points
  //	assert(LiThePoints.iLength() >= 3);
  // The verts HAVE to be valid..
  assert(LiThePoints[0] >= 0);
  assert(LiThePoints[LiThePoints.iLength()-1] >= 0);
  // Copy all the points
  LiPoints.vClear();
  LiPoints.vAppendList(LiThePoints);
  // Set boundary type
  eGeom2DType = eInterpolated;
  // We will have LiPoints.iLength()-1 splines...
  dMaxParam = LiPoints.iLength()-1;

  // Must now create all the cubic parametric boundary patches
  vComputeAllSplines();

  // Refresh...
  vRefresh();
}

void BdryInterpolated2D::
vComputeAllSplines()
{
  // We have at least three points in LiPoints
  // Compute the coefficients needed to initialize the splines.
  int iNP = LiPoints.iLength();
  // Set our arrays needed for this...
  double *adRX=NULL, *adRY=NULL, adPt0[2], adPt1[2], adPt2[2];

  if (iNP > 2) {
    adRX = new double[iNP-2];
    adRY = new double[iNP-2];

    // Assign the proper values to the rows...
    for (int i = 0; i < iNP-2; i++) {
      // Get the points...
      pB2D->vGetPoint(LiPoints[i],   adPt0);
      pB2D->vGetPoint(LiPoints[i+1], adPt1);
      pB2D->vGetPoint(LiPoints[i+2], adPt2);
      // Set the result column..
      adRX[i] = 3.0 * (adPt0[0] - 2.0*adPt1[0] + adPt2[0]);
      adRY[i] = 3.0 * (adPt0[1] - 2.0*adPt1[1] + adPt2[1]);
    }

    // Assume end B's are zero for now: have take them into account otherwise

    // Solve the tri-diagonal systems
    vFindBCoeffs(iNP-2, adRX);
    vFindBCoeffs(iNP-2, adRY);
  }

  // Now create the cubic parametric patches needed
  double dAx, dBx, dCx, dDx, dAy, dBy, dCy, dDy;
  double dBnx, dBnpx, dBny, dBnpy;
  double adPtn[2], adPtnp[2];

  // Calculate the coefficients and create the patches
  for (int i = 0; i < iNP-1; i++) {
    // Get the points
    pB2D->vGetPoint(LiPoints[i], adPtn);
    pB2D->vGetPoint(LiPoints[i+1], adPtnp);

    // Set the bn's and bn-1's
    // Must be careful with the extremities
    if (i == 0) {
      // First b's
      dBnx = 0.0;
      dBny = 0.0;
    }
    else {
      dBnx  = adRX[i-1];
      dBny  = adRY[i-1];
    }

    if (i == iNP-2) {
      // Last b's
      dBnpx = 0.0;
      dBnpy = 0.0;
    }
    else {
      dBnpx = adRX[i];
      dBnpy = adRY[i];
    }

    // Set the coefficients.

    dAx = (dBnpx - dBnx)/3.0;
    dBx = dBnx;
    dCx = (adPtnp[0]-adPtn[0]) - ((dBnpx + 2.0*dBnx)/3.0);
    dDx = adPtn[0];

    dAy = (dBnpy - dBny)/3.0;
    dBy = dBny;
    dCy = (adPtnp[1]-adPtn[1]) - ((dBnpy + 2.0*dBny)/3.0);
    dDy = adPtn[1];

    // Call the constructor...

    BdryCubicParam2D *pBCP =
      new BdryCubicParam2D(pB2D, iLeftBdryCond(), iRightBdryCond(),
			   iLeftRegion(), iRightRegion(),
			   LiPoints[i], LiPoints[i+1],
			   dAx, dBx, dCx, dDx, dAy, dBy, dCy, dDy);

    LpBCP2D.vAppendItem(pBCP);
  }

  if (iNP > 2) {
    delete [] adRX;
    delete [] adRY;
  }
}

void BdryInterpolated2D::
vFindBCoeffs(const int iNumRows, double *adR)
{
  double *adA, *adB, *adC;

  adA = new double[iNumRows];
  adB = new double[iNumRows];
  adC = new double[iNumRows];

  int i;
  // Fill the tridiagonal matrix
  for (i = 0; i < iNumRows; i++) {
    adA[i] = 1.0;
    adB[i] = 4.0;
    adC[i] = 1.0;
  }

  // Use Thomas algorithm to solve system
  for (i = 0; i < iNumRows-1; i++) {
    adC[i] /= adB[i];
    adR[i] /= adB[i];
    adB[i+1] -= adC[i]*adA[i+1];
    adR[i+1] -= adR[i]*adA[i+1];
  }
  adR[iNumRows-1] /= adB[iNumRows-1];

  for (i = iNumRows-2; i >= 0; i--) {
    adR[i] -= adR[i+1]*adC[i];
  }

  delete [] adA;
  delete [] adB;
  delete [] adC;
}

void BdryInterpolated2D::
vRefresh()
{
  dOrientationChange = 0.0;
  LdOrderedRoots.vClear();
  LdOrientationChangeAtRoots.vClear();
  // We have a list of splines -- many things need to be calculated.
  for (int i = 0; i < LpBCP2D.iLength(); i++) {
    BdryCubicParam2D *pBCP = LpBCP2D[i];
    // Make sure they are up to date..
    pBCP->vRefresh();

    double di = i;

    // Roots need to be updated too..
    int iNumRoots = pBCP->iNumOrderedRoots();

    // Insert the very first one...
    if (i == 0) {
      int iNB = LdOrderedRoots.iLength();
      LdOrderedRoots.vAddItem(0.0);
      if (LdOrderedRoots.iLength() > iNB)
	LdOrientationChangeAtRoots.vAppendItem(0.0);
    }

    for (int j = 1; j < iNumRoots+2; j++) {
      double dPatchRoot = pBCP->dGetOrderedRoot(j);
      double dPatchChange = pBCP->dGetOrientationChangeAtRoot(j);
      int iNB = LdOrderedRoots.iLength();
      LdOrderedRoots.vAddItem(di + dPatchRoot);
      if (LdOrderedRoots.iLength() > iNB)
	LdOrientationChangeAtRoots.vAppendItem(dOrientationChange + dPatchChange);
    }

    // Update orientation change
    dOrientationChange += pBCP->dGetOrientationChange();
  }
  /*
    vMessage(3, "Total orientation change for spline: %f\n", dOrientationChange);
    if (LdOrderedRoots.iLength() > 0) {
    vMessage(3, "Here is a list of the ordered roots\n");
    for (int i = 0; i < LdOrderedRoots.iLength(); i++ ) {
    vMessage(3, "  %d: %f -- %f\n", i, LdOrderedRoots[i], LdOrientationChangeAtRoots[i]);
    }
    }
  */
}

double BdryInterpolated2D::
dOrientationChangeAtParameter(const double dt) const
{
  int iSpline;
  double dRemaining;
  vGetSplineFromParameter(dt, &iSpline, &dRemaining);
  // Get the orientation at the beginning of the spline
  double dOrChangeAtStart;
  int iIndex = LdOrderedRoots.iSeek(iSpline);
  dOrChangeAtStart = LdOrientationChangeAtRoots[iIndex];
  return dOrChangeAtStart + LpBCP2D[iSpline]->dOrientationChangeAtParameter(dRemaining);
}


double BdryInterpolated2D::
dParameterAtOrientationChange(const double dChange, const double dStart) const
{
  // 1) Find between which roots to look at.
  // 2) Find the parametric curve to refer to...
  // 3) Ask the parametric curve to do the job for you!

  // Find the biggest root smaller (or equal) to dStart
  int iStart = 0;
  bool qFound = false;
  double dOrStart = 0.0;

  assert ( (dStart >= 0.0) && (dStart <= dMaxParam) );

  while (!qFound) {
    if (LdOrderedRoots[iStart] > dStart) {
      qFound = true;
      iStart--;
    }
    else {
      iStart++;
    }
  }

  dOrStart = dOrientationChangeAtParameter(dStart);

  //	vMessage(3, "%f is just less than %f\n", LdOrderedRoots[iStart], dStart);
  //	vMessage(3, "Orientation change at root is %f\n", dOrStart);

  // Now find the roots with an orientation change just smaller, and just bigger
  qFound = false;
  int iLower = iStart;
  int iUpper = iStart;

  while (!qFound) {
    iUpper++;
    if ( iFuzzyComp(LdOrientationChangeAtRoots[iUpper] - dOrStart,
		    dChange) >= 0) {
      qFound = true;
    }
    else {
      iLower++;
    }
  }

  //	vMessage(3, "Orientation change limits are %f and %f\n",
  //		LdOrientationChangeAtRoots[iLower], LdOrientationChangeAtRoots[iUpper]);
  //	vMessage(3, "For which the roots are %f and %f\n",
  //		LdOrderedRoots[iLower], LdOrderedRoots[iUpper]);

  int iSpline;
  double dRemaining;
  double dNewChange = (dOrStart + dChange) - LdOrientationChangeAtRoots[iLower];
  vGetSplineFromParameter(LdOrderedRoots[iLower], &iSpline, &dRemaining);

  return iSpline + LpBCP2D[iSpline]->dParameterAtOrientationChange(dNewChange, dRemaining);

}

double BdryInterpolated2D::
dCurvature(const double adLoc[2]) const
{
  // Get parameter for the point

  double dt = dParameterAtPoint(adLoc);

  // Find which spline we have to look at...
  double dParamLeft;
  int    iSpline;
  vGetSplineFromParameter(dt, &iSpline, &dParamLeft);

  // Call the proper function..
  return LpBCP2D[iSpline]->dCurvature(adLoc);
  // as simple as that..

}

double BdryInterpolated2D::
dParameterAtPoint(const double adLoc[2]) const
{
  // Let's try the following approach..
  // 1) Get the distance from all the cubic parametric patches
  // 2) Only look at the ones that report it to be very close...

  int iClosestPatch = -1;
  double dDist;
  double dSmallestDist = 1e20;

  //	vMessage(3, "Looking at (%f, %f)\n", adLoc[0], adLoc[1]);

  for (int iP = 0; iP < LpBCP2D.iLength(); iP++) {
    dDist = LpBCP2D[iP]->dDistFromBdry(adLoc);
    //		vMessage(3, "Patch: %d -- Distance: %3.12e\n", iP, dDist);

    if (iFuzzyComp(fabs(dDist), fabs(dSmallestDist)) < 0)	{
      //			vMessage(3, "That was the smallest yet\n");
      dSmallestDist = dDist;
      iClosestPatch = iP;
    }
  }

  //	vMessage(3, "Closest patch: %d -- Distance: %3.18e\n", iClosestPatch, dSmallestDist);

  assert(iClosestPatch >= 0);

  double dParaminPatch;

  LpBCP2D[iClosestPatch]->vParameterClosestToPoint(adLoc, &dParaminPatch);

  return iClosestPatch + dParaminPatch;
}

void   BdryInterpolated2D::
vGetSplineFromParameter(const double dParam, int *iSpline, double *dRemainingAL) const
{
  *iSpline = static_cast<int>(floor(dParam));
  *dRemainingAL = dParam - (*iSpline);

  if ( (iFuzzyComp(*dRemainingAL, 0.0) == 0) && (*iSpline == LpBCP2D.iLength()) ) {
    // We're at the very end of the interpolated spline...
    *iSpline = *iSpline - 1;
    *dRemainingAL = 1.0;
  }

}

double BdryInterpolated2D::dGeodesicDistance(const double dParamIn1,
					     const double dParamIn2) const
{
  double dParam1, dParam2;
  if (dParamIn1 < dParamIn2) {
    dParam1 = dParamIn1;
    dParam2 = dParamIn2;
  }
  else {
    dParam2 = dParamIn1;
    dParam1 = dParamIn2;
  }

  // Find out which splines these two fall into:
  double dParamExtra1, dParamExtra2;
  int    iSpline1,     iSpline2;

  assert(dParam1 <= dParam2);

  vGetSplineFromParameter(dParam1, &iSpline1, &dParamExtra1);
  vGetSplineFromParameter(dParam2, &iSpline2, &dParamExtra2);

  assert(iSpline1 <= iSpline2);
  assert(iFuzzyComp(dParamExtra1, 0.0) >= 0 &&
	 iFuzzyComp(dParamExtra1, 1.0) <= 0);
  assert(iFuzzyComp(dParamExtra2, 0.0) >= 0 &&
	 iFuzzyComp(dParamExtra2, 1.0) <= 0);

  if (iSpline1 == iSpline2) {
    return LpBCP2D[iSpline1]->dGeodesicDistance(dParamExtra1, dParamExtra2);
  }
  else {
    vMessage(4, "iSpline1 = %2d, iSpline2 = %2d\n", iSpline1, iSpline2);
    // Finish off this subspline
    double dRetVal = LpBCP2D[iSpline1]->dGeodesicDistance(dParamExtra1, 1.0);
    for (int ii = iSpline1 + 1; ii < iSpline2; ii++) {
      dRetVal += LpBCP2D[ii]->dGeodesicDistance(0.0, 1.0);
    }
    dRetVal += LpBCP2D[iSpline2]->dGeodesicDistance(0.0, dParamExtra2);
    return dRetVal;
  }
}

double BdryInterpolated2D::dSpeed(const double dParam) const
{
  // Find which spline we have to look at...
  double dParamLeft;
  int    iSpline;
  vGetSplineFromParameter(dParam, &iSpline, &dParamLeft);

  // Call the proper function..
  return LpBCP2D[iSpline]->dSpeed(dParamLeft);
  // as simple as that..
}

void   BdryInterpolated2D::
vPointAtParameter(const double dParam, double adLoc[2]) const
{
  // Find which spline we have to look at...
  double dParamLeft;
  int    iSpline;
  vGetSplineFromParameter(dParam, &iSpline, &dParamLeft);

  // Call the proper function..
  LpBCP2D[iSpline]->vPointAtParameterRatio(dParamLeft, adLoc);
  // as simple as that..
}


void   BdryInterpolated2D::
vPointAtParameterRatio(const double dRatio, double adLoc[2]) const
{
  // Find which spline we have to look at...
  double dParamLeft;
  int    iSpline;
  vGetSplineFromParameter(dRatio * dMaxParam, &iSpline, &dParamLeft);

  // Call the proper function..
  LpBCP2D[iSpline]->vPointAtParameterRatio(dParamLeft, adLoc);
  // as simple as that..
}


void   BdryInterpolated2D::
vUnitNormal(const double adLoc[], double adNorm[]) const
{
  double dt = dParameterAtPoint(adLoc);
  double dRatio = dt/dMaxParam;
  vNormalAtParameterRatio(dRatio, adNorm);
}

void   BdryInterpolated2D::
vNormalAtParameterRatio(const double dRatio, double adNorm[2]) const
{
  // Find which spline we have to look at...
  double dParamLeft;
  int    iSpline;
  vGetSplineFromParameter(dRatio * dMaxParam, &iSpline, &dParamLeft);

  // Call the proper function..
  LpBCP2D[iSpline]->vNormalAtParameterRatio(dParamLeft, adNorm);
  // as simple as that..

}

void   BdryInterpolated2D::
vTangentAtParameterRatio(const double dRatio, double a2dTangent[2]) {
  // Find which spline we have to look at...
  double dParamLeft;
  int    iSpline;
  vGetSplineFromParameter(dRatio * dMaxParam, &iSpline, &dParamLeft);

  // Call the proper function..
  LpBCP2D[iSpline]->vTangentAtParameterRatio(dParamLeft, a2dTangent);
  // as simple as that..

}


void   BdryInterpolated2D::
vAssignData(Bdry2D*, const int, const int, const int, const int,
	    const int, const int)
{
  assert(0);
}

void   BdryInterpolated2D::
vPointAtDistFromStart(const double dDist, double adLoc[2]) const
{
  // Find roots at a distance just smaller and just greater than
  // dDist...

  bool qFoundHigher = false;
  bool qFoundIt = false;
  int  iLower=1000;
  int  iHigher=-1000;
  int  i = -1;
  double dRoot;
  double dRootDist;
  int  iFC;

  double adStart[2];
  double adRootPt[2];

  vGetBeginPoint(adStart);

  do {
    i++;
    dRoot = LdOrderedRoots[i];
    vPointAtParameter(dRoot, adRootPt);
    dRootDist = dDIST2D(adStart, adRootPt);

    iFC = iFuzzyComp(dDist, dRootDist);

    if (iFC == 0) {
      // Right on it!
      qFoundIt = true;
    }

    if (iFC < 0) {
      iHigher = i;
      qFoundHigher = true;
    }
  } while (!qFoundHigher && !qFoundIt);

  if (qFoundIt) {
    adLoc[0] = adRootPt[0];
    adLoc[1] = adRootPt[1];
    return;
  }

  assert(iHigher > 0);
  iLower = iHigher - 1;

  double dLowDist, dHighDist, dNewDist;
  double dLow, dHigh, dNew, dRatio;

  dLow = LdOrderedRoots[iLower];
  dHigh = LdOrderedRoots[iHigher];

  vPointAtParameter(dLow, adRootPt);
  dLowDist = dDIST2D(adStart, adRootPt);
  vPointAtParameter(dHigh, adRootPt);
  dHighDist = dDIST2D(adStart, adRootPt);

  // Now have lower and higher limits... Zoom in on the point...
  while (!qFoundIt) {
    dRatio = fabs(dLowDist - dDist) / ( fabs(dLowDist - dDist) + fabs(dHighDist - dDist) );
    dNew = dLow + dRatio*(dHigh - dLow);

    vPointAtParameter(dNew, adRootPt);
    dNewDist = dDIST2D(adStart, adRootPt);

    iFC = iFuzzyComp(dDist, dNewDist);

    if (iFC < 0) {
      // New distance is greater than what it should be..
      dHigh = dNew;
      dHighDist = dNewDist;
    }

    if (iFC > 0) {
      // New distance is less than what it should be..
      dLow = dNew;
      dLowDist = dNewDist;
    }

    if (iFC == 0) {
      qFoundIt = true;
    }
  }

  adLoc[0] = adRootPt[0];
  adLoc[1] = adRootPt[1];

}

void   BdryInterpolated2D::
vPointAtDistFromEnd(const double dDist, double adLoc[2]) const
{
  // Find roots at a distance just smaller and just greater than
  // dDist...

  bool qFoundLower = false;
  bool qFoundIt = false;
  int  iLower=1000;
  int  iHigher=-1000;
  int  i = LdOrderedRoots.iLength();
  double dRoot;
  double dRootDist;
  int  iFC;

  double adEnd[2];
  double adRootPt[2];

  vGetEndPoint(adEnd);

  do {
    i--;
    dRoot = LdOrderedRoots[i];
    vPointAtParameter(dRoot, adRootPt);
    dRootDist = dDIST2D(adEnd, adRootPt);

    iFC = iFuzzyComp(dDist, dRootDist);

    if (iFC == 0) {
      // Right on it!
      qFoundIt = true;
    }

    if (iFC < 0) {
      iLower = i;
      qFoundLower = true;
    }
  } while (!qFoundLower && !qFoundIt);

  if (qFoundIt) {
    adLoc[0] = adRootPt[0];
    adLoc[1] = adRootPt[1];
    return;
  }

  assert(iLower < LdOrderedRoots.iLength()-1);
  iHigher = iLower + 1;

  double dLowDist, dHighDist, dNewDist;
  double dLow, dHigh, dNew, dRatio;

  dLow = LdOrderedRoots[iLower];
  dHigh = LdOrderedRoots[iHigher];

  vPointAtParameter(dLow, adRootPt);
  dLowDist = dDIST2D(adEnd, adRootPt);
  vPointAtParameter(dHigh, adRootPt);
  dHighDist = dDIST2D(adEnd, adRootPt);

  // Now have lower and higher limits... Zoom in on the point...
  while (!qFoundIt) {
    dRatio = fabs(dLowDist - dDist) / ( fabs(dLowDist - dDist) + fabs(dHighDist - dDist) );
    dNew = dLow + dRatio*(dHigh - dLow);

    vPointAtParameter(dNew, adRootPt);
    dNewDist = dDIST2D(adEnd, adRootPt);

    iFC = iFuzzyComp(dDist, dNewDist);

    if (iFC < 0) {
      // New distance is greater than what it should be..
      dLow = dNew;
      dLowDist = dNewDist;
    }

    if (iFC > 0) {
      // New distance is less than what it should be..
      dHigh = dNew;
      dHighDist = dNewDist;
    }

    if (iFC == 0) {
      qFoundIt = true;
    }
  }

  adLoc[0] = adRootPt[0];
  adLoc[1] = adRootPt[1];
}

void   BdryInterpolated2D::
vBdryPointClosestToPoint(const double /*adLoc*/[2],
			 double /*adPoint*/[2]) const
{
  vMessage(1, "CALLED EMPTY FUNCTION BdryInterpolated::vBdryPointClosestToPoint()!\n");
}

double  BdryInterpolated2D::
dDistFromBdry(const double adLoc[2], int * const piSide) const
{
  // Call all of the splines...
  // Return the smallest distance (and the side returned along with that
  // distance)

  int iP;
  double dNewSmallestDist = LARGE_DBL;
  double dD;
  int    iS, iR;
  int    iNewSmallestRoot = -1;
  int    iNewSmallestSide = SIDE_UNDETERMINED;
  BdryCubicParam2D *pBCP;

  int    *aiCross = new int   [LdOrderedRoots.iLength()];

  for (iR = 0; iR < LdOrderedRoots.iLength(); iR++) {
    // Get each root.
    // Calculate distance vector
    // Get normal too
    // Calculate dot product.
    double adPt[2];
    double adNorm[2];
    double adVec[2];

    vPointAtParameter(LdOrderedRoots[iR], adPt);
    vNormalAtParameterRatio(LdOrderedRoots[iR]/dMaxParam, adNorm);
    // Calculate vector
    adVec[0] = adPt[0] - adLoc[0];
    adVec[1] = adPt[1] - adLoc[1];
    // Cross product..
    double dCross = adVec[0]*adNorm[1] - adVec[1]*adNorm[0];
    double dDist  = dDIST2D(adPt, adLoc);

    if (iFuzzyComp(dDist, dNewSmallestDist) < 0) {
      dNewSmallestDist = dDist;
      iNewSmallestRoot = iR;
    }

    aiCross[iR] = iFuzzyComp(dCross, 0.0);
  }

  // Now find cubics that are candidates
  List<int> LiCubics;
  int iLast = aiCross[0];
  int iNow;

  for (iR = 1; iR < LdOrderedRoots.iLength(); iR++) {
    iNow = aiCross[iR];

    if (iNow != iLast) {
      // Intervals point to the first item
      double dDum;
      int iCubicLast, iCubicNow;
      vGetSplineFromParameter(LdOrderedRoots[iR-1], &iCubicLast, &dDum);
      vGetSplineFromParameter(LdOrderedRoots[iR], &iCubicNow, &dDum);
      LiCubics.vAddItem(iCubicLast);
      LiCubics.vAddItem(iCubicNow);
    }

    iLast = iNow;
  }

  if (LiCubics.iLength() > 0) {
    // Now only need to check the splines we isolated
    for (iP = 0; iP < LiCubics.iLength(); iP++) {
      pBCP = LpBCP2D[LiCubics[iP]];
      dD = pBCP->dDistFromBdry(adLoc, &iS);
      if (iFuzzyComp(fabs(dD), dNewSmallestDist) < 0) {
	dNewSmallestDist = dD;
	iNewSmallestSide = iS;
      }
    }
  }
  else {
    // Simply return the smallest distance then...
    int iSpline;
    double dDum;
    vGetSplineFromParameter(LdOrderedRoots[iNewSmallestRoot], &iSpline, &dDum);
    pBCP = LpBCP2D[iSpline];
    dNewSmallestDist = pBCP->dDistFromBdry(adLoc, &iNewSmallestSide);
  }
  if (piSide) *piSide = iNewSmallestSide;
  delete aiCross;
  return dNewSmallestDist;
}

void   BdryInterpolated2D::
vRatioBetweenTwoPoints(const double dRatio, const double adLoc1[2],
		       const double adLoc2[2], double adPoint[2]) const
{
  // - Find the parameter for the two points (dt1, dt2)
  // - Find the total orientation for each point (dOr1, dOr2)
  // - Find the point with orientation change dRatio*(dOr2-dOr1) starting from
  //   dOr1

  assert ( (dRatio >= 0.0) && (dRatio <= 1.0) );

  double dCorrectedRatio = dRatio;
  double dtdum; // dummy var
  // Get parameters
  double dt1 = dParameterAtPoint(adLoc1);
  //	vMessage(3, "dt1 is %f\n", dt1);
  double dt2 = dParameterAtPoint(adLoc2);
  //	vMessage(3, "dt2 is %f\n", dt2);
  // Also need to get the orientation changes there...
  double dOr1 = dOrientationChangeAtParameter(dt1);
  double dOr2 = dOrientationChangeAtParameter(dt2);

  // Order them properly
  if (dt1 > dt2) {
    dtdum = dt2;
    dt2 = dt1;
    dt1 = dtdum;
    dtdum = dOr2;
    dOr2 = dOr1;
    dOr1 = dtdum;
    dCorrectedRatio = 1-dCorrectedRatio;
  }

  // What if the spline is closed on itself?
  if (iBeginPointIndex() == iEndPointIndex() ) {
    if (iFuzzyComp(dt1, 0.0) == 0) {
      // Maybe dt1 should be dMaxParam instead..
      if (dOr2 > dOr1 + 1.0) {
	// It should have been the other way around..
	dt1 = dMaxParam;
	dOr1 = dOrientationChange;
      }
    }
    else {
      if (iFuzzyComp(dt2, dMaxParam) == 0) {
	// Maybe dt2 should be 0.0 instead...
	if (dOr1 < dOr2 - 1.0) {
	  dt2 = 0.0;
	  dOr2 = 0.0;
	}
      }
    }
  }

  // Re-Order them properly
  if (dt1 > dt2) {
    dtdum = dt2;
    dt2 = dt1;
    dt1 = dtdum;
    dtdum = dOr2;
    dOr2 = dOr1;
    dOr1 = dtdum;
    dCorrectedRatio = 1-dCorrectedRatio;
  }


  // Calculate new location
  if (fabs(dOr2-dOr1) < 0.01) {
    dtdum = dt1 + dCorrectedRatio*(dt2-dt1);
  }
  else {
    dtdum = dParameterAtOrientationChange(dCorrectedRatio*fabs(dOr2-dOr1), dt1);
  }
  // Get point
  assert( iFuzzyComp(dtdum, dt1) >= 0 &&
	  iFuzzyComp(dtdum, dt2) <= 0 );
  vPointAtParameter(dtdum, adPoint);
  // done
}

double BdryInterpolated2D::
dXMax()
{
  double dMax = -1E16;
  double dPatchMax;
  for (int i = 0; i < LpBCP2D.iLength(); i++) {
    dPatchMax = LpBCP2D[i]->dXMax();
    if (dPatchMax > dMax) dMax = dPatchMax;
  }

  return dMax;
}

double BdryInterpolated2D::
dYMax()
{
  double dMax = -1E16;
  double dPatchMax;
  for (int i = 0; i < LpBCP2D.iLength(); i++) {
    dPatchMax = LpBCP2D[i]->dYMax();
    if (dPatchMax > dMax) dMax = dPatchMax;
  }

  return dMax;
}

double BdryInterpolated2D::
dXMin()
{
  double dMin = 1E16;
  double dPatchMin;
  for (int i = 0; i < LpBCP2D.iLength(); i++) {
    dPatchMin = LpBCP2D[i]->dXMin();
    if (dPatchMin < dMin) dMin = dPatchMin;
  }

  return dMin;
}

double BdryInterpolated2D::
dYMin()
{
  double dMin = 1E16;
  double dPatchMin;
  for (int i = 0; i < LpBCP2D.iLength(); i++) {
    dPatchMin = LpBCP2D[i]->dYMin();
    if (dPatchMin < dMin) dMin = dPatchMin;
  }

  return dMin;
}

List<double>  BdryInterpolated2D::
LdOriginalDiscretization(
			 bool *qInsertVertPointers,
			 bool *qCloseLoop
			 )
{
  // Determine the original discretization of the patch.
  // Works like this:
  // - Determine the total (absolute) change in orientation of the patch
  // - Number of points to insert = (int)(total change in orientation)
  // - Find the proper locations, i.e. at every (change in orientation)/(numpoints+1)
  //   (that one will be interesting)

  int iNPts = 2*static_cast<int>(floor(dOrientationChange)) + 1;
  List<double> LdPointList;
  double adPt[2];
  double dPrevt = 0.0;
  double dNewt = 0.0;

  LdPointList.vClear();

  // In general, these are set this way...
  *qCloseLoop = false;
  *qInsertVertPointers = true;

  // Orientation change along each discretization segment
  double dChange = dOrientationChange/(iNPts + 1);

  // Find the other points in the middle
  for (int k = 1; k <= iNPts; k++) {
    // Get t
    dNewt = dParameterAtOrientationChange(dChange, dPrevt);
    // Find point
    vPointAtParameter(dNewt, adPt);
    // Insert in list
    //		vMessage(2, "t: %f, Or: %f (%f, %f)\n", dNewt, k*dChange, adPt[0], adPt[1]);
    //		vMessage(2, "But I find this when I try to find back:\n");
    LdPointList.vAppendItem(adPt[0]);
    LdPointList.vAppendItem(adPt[1]);
    // Make sure to update dPrevt
    dPrevt = dNewt;
  }
  /*
    if ( iBeginPointIndex() == iEndPointIndex() ) {
    // The curve loops on itself
    *qCloseLoop = true;
    }
  */
  return LdPointList;

}

void   BdryInterpolated2D::
vGetSplineCoefficients(const int iSpline,
		       double *dAx, double *dBx, double *dCx, double *dDx,
		       double *dAy, double *dBy, double *dCy, double *dDy)
{
  LpBCP2D[iSpline]->vGetAllCoefficients(dAx, dBx, dCx, dDx, dAy, dBy, dCy, dDy);
}

void   BdryInterpolated2D::
vSetSplineCoefficients(const int iSpline, const double dAx,
		       const double dBx, const double dCx, const double dDx,
		       const double dAy, const double dBy, const double dCy, const double dDy)
{
  LpBCP2D[iSpline]->vSetAllCoefficients(dAx, dBx, dCx, dDx, dAy, dBy, dCy, dDy);
  LpBCP2D[iSpline]->vRefresh();
}


void   BdryInterpolated2D::
vSplitSplineAtParameter(const int /*iSpline*/, const double dSplit,
			double *dAx1, double *dBx1, double *dCx1, double *dDx1,
			double *dAy1, double *dBy1, double *dCy1, double *dDy1,
			double *dAx2, double *dBx2, double *dCx2, double *dDx2,
			double *dAy2, double *dBy2, double *dCy2, double *dDy2 )
{
  // Splits the spline at parameter dSplit, and returns the new parameters
  // of the two resulting cubic parametric curves

  double dcAx, dcAy, dcBx, dcBy, dcCx, dcCy, dcDx, dcDy;

  // The coefficients to use are in the first set...
  dcAx = *dAx1; dcBx = *dBx1; dcCx = *dCx1; dcDx = *dDx1;
  dcAy = *dAy1; dcBy = *dBy1; dcCy = *dCy1; dcDy = *dDy1;

  // For the first part...

  double dtc, dtc2, dtc3;

  assert ( (dSplit >= 0.0) && (dSplit <= 1.0) );

  if (dSplit < 1.0) {

    dtc = dSplit;
    dtc2 = dtc*dtc;
    dtc3 = dtc2*dtc;

    *dAx1 = dcAx*dtc3; *dAy1 = dcAy*dtc3;
    *dBx1 = dcBx*dtc2; *dBy1 = dcBy*dtc2;
    *dCx1 = dcCx*dtc;  *dCy1 = dcCy*dtc;
    *dDx1 = dcDx;      *dDy1 = dcDy;

    // For the second part...

    *dAx2 = dcAx*(1.0 - 3.0*dtc + 3.0*dtc2 - dtc3);
    *dAy2 = dcAy*(1.0 - 3.0*dtc + 3.0*dtc2 - dtc3);
    *dBx2 = dcAx*(3.0*dtc - 6.0*dtc2 + 3.0*dtc3) + dcBx*(1.0 - 2.0*dtc + dtc2);
    *dBy2 = dcAy*(3.0*dtc - 6.0*dtc2 + 3.0*dtc3) + dcBy*(1.0 - 2.0*dtc + dtc2);
    *dCx2 = dcAx*(-3.0*dtc3 + 3.0*dtc2) + dcBx*(-2.0*dtc2 + 2.0*dtc) + dcCx*(1.0 - dtc);
    *dCy2 = dcAy*(-3.0*dtc3 + 3.0*dtc2) + dcBy*(-2.0*dtc2 + 2.0*dtc) + dcCy*(1.0 - dtc);
    *dDx2 = (dcAx*dtc3 + dcBx*dtc2 + dcCx*dtc + dcDx);
    *dDy2 = (dcAy*dtc3 + dcBy*dtc2 + dcCy*dtc + dcDy);
  }

}

double BdryInterpolated2D::
dCalculateNewShortSplinesChange() {

  // returns the number of splines and
  // adds the new control points to the
  // B2D object
  int iNumShortSplines = static_cast<int>(floor(dOrientationChange));

  // Orientation change along each discretization segment
  double dChange = dOrientationChange/(iNumShortSplines + 1);

  return dChange;
}

/*
  void   BdryInterpolated2D::
  vGetNewShortSpline(const double dStart, const double dChange, double *dEnd, List<int> *LiPoints, List<double> *LdCoeffs)
  {
  // Find in which spline dStart lies in
  int iStartSpline;
  double dStartRemaining;
  vGetSplineFromParameter(dStart, &iStartSpline, &dStartRemaining);

  if (iFuzzyComp(dStartRemaining, 0.0) == 0) {
  if (iStartSpline > 0) {
  iStartSpline--;
  dStartRemaining = 1.0;
  }
  }

  // Find the new location...
  double dNewt = dParameterAtOrientationChange(dChange, dStart);
  *dEnd = dNewt;

  int iEndSpline;
  double dEndRemaining;
  vGetSplineFromParameter(dNewt, &iEndSpline, &dEndRemaining);

  // Split!
  double dAx1, dBx1, dCx1, dDx1, dAy1, dBy1, dCy1, dDy1;
  double dAx2, dBx2, dCx2, dDx2, dAy2, dBy2, dCy2, dDy2;
  int    iNewPoint = 0;

  if ((iStartSpline == iEndSpline) && (qUseLastCoeffs)) {
  // This is a hint to replace the current params with
  // the given ones..
  iNewPoint = -1;
  dAx1 = dLastCoeffs[0][0]; dAy1 = dLastCoeffs[1][0];
  dBx1 = dLastCoeffs[0][1]; dBy1 = dLastCoeffs[1][1];
  dCx1 = dLastCoeffs[0][2]; dCy1 = dLastCoeffs[1][2];
  dDx1 = dLastCoeffs[0][3]; dDy1 = dLastCoeffs[1][3];
  }

  vSplitSplineAtParameter(dNewt, &iNewPoint,
  &dAx1, &dBx1, &dCx1, &dDx1, &dAy1, &dBy1, &dCy1, &dDy1,
  &dAx2, &dBx2, &dCx2, &dDx2, &dAy2, &dBy2, &dCy2, &dDy2);

  LiPoints->vClear();
  LdCoeffs->vClear();

  if (qUseLastCoeffs) {
  // Add point
  LiPoints->vAppendItem(iLastNewPoint);

  if (iStartSpline != iEndSpline) {
  // add coeffs too
  LdCoeffs->vAppendItem(dLastCoeffs[0][0]);
  LdCoeffs->vAppendItem(dLastCoeffs[0][1]);
  LdCoeffs->vAppendItem(dLastCoeffs[0][2]);
  LdCoeffs->vAppendItem(dLastCoeffs[0][3]);
  LdCoeffs->vAppendItem(dLastCoeffs[1][0]);
  LdCoeffs->vAppendItem(dLastCoeffs[1][1]);
  LdCoeffs->vAppendItem(dLastCoeffs[1][2]);
  LdCoeffs->vAppendItem(dLastCoeffs[1][3]);
  }

  // Add the begin points of remaining splines
  for (int i = iStartSpline+1; i <= iEndSpline; i++) {
  LiPoints->vAppendItem(LpBCP2D[i]->iBeginPointIndex());
  // Don't add the last coefficients
  if (i != iEndSpline) {
  double d1, d2, d3, d4, d5, d6, d7, d8;
  LpBCP2D[i]->vGetAllCoefficients(&d1, &d2, &d3, &d4, &d5, &d6, &d7, &d8);
  // Add the coefficients
  LdCoeffs->vAppendItem(d1);
  LdCoeffs->vAppendItem(d2);
  LdCoeffs->vAppendItem(d3);
  LdCoeffs->vAppendItem(d4);
  LdCoeffs->vAppendItem(d5);
  LdCoeffs->vAppendItem(d6);
  LdCoeffs->vAppendItem(d7);
  LdCoeffs->vAppendItem(d8);
  }
  }

  if (iFuzzyComp(dEndRemaining, 0.0) != 0) {
  // add the last coefficients.....
  if ( (iFuzzyComp(dNewt, dMaxParam) != 0) ||
  (iStartSpline == iEndSpline) ) {
  LdCoeffs->vAppendItem(dAx1);
  LdCoeffs->vAppendItem(dBx1);
  LdCoeffs->vAppendItem(dCx1);
  LdCoeffs->vAppendItem(dDx1);
  LdCoeffs->vAppendItem(dAy1);
  LdCoeffs->vAppendItem(dBy1);
  LdCoeffs->vAppendItem(dCy1);
  LdCoeffs->vAppendItem(dDy1);
  }

  LiPoints->vAppendItem(iNewPoint);
  }
  }
  else {
  // only when we start from 0.0 ...
  for (int i = iStartSpline; i <= iEndSpline; i++) {
  LiPoints->vAppendItem(LpBCP2D[i]->iBeginPointIndex());
  // Don't add the last ones...
  if (i != iEndSpline) {
  double d1, d2, d3, d4, d5, d6, d7, d8;
  LpBCP2D[i]->vGetAllCoefficients(&d1, &d2, &d3, &d4, &d5, &d6, &d7, &d8);
  // Add the coefficients
  LdCoeffs->vAppendItem(d1);
  LdCoeffs->vAppendItem(d2);
  LdCoeffs->vAppendItem(d3);
  LdCoeffs->vAppendItem(d4);
  LdCoeffs->vAppendItem(d5);
  LdCoeffs->vAppendItem(d6);
  LdCoeffs->vAppendItem(d7);
  LdCoeffs->vAppendItem(d8);
  }
  }

  if (iFuzzyComp(dEndRemaining, 0.0) != 0) {
  // add the last coefficients.....
  LdCoeffs->vAppendItem(dAx1);
  LdCoeffs->vAppendItem(dBx1);
  LdCoeffs->vAppendItem(dCx1);
  LdCoeffs->vAppendItem(dDx1);
  LdCoeffs->vAppendItem(dAy1);
  LdCoeffs->vAppendItem(dBy1);
  LdCoeffs->vAppendItem(dCy1);
  LdCoeffs->vAppendItem(dDy1);

  LiPoints->vAppendItem(iNewPoint);
  }
  }

  qUseLastCoeffs = true;
  iLastNewPoint = iNewPoint;
  dLastCoeffs[0][0] = dAx2; dLastCoeffs[1][0] = dAy2;
  dLastCoeffs[0][1] = dBx2; dLastCoeffs[1][1] = dBy2;
  dLastCoeffs[0][2] = dCx2; dLastCoeffs[1][2] = dCy2;
  dLastCoeffs[0][3] = dDx2; dLastCoeffs[1][3] = dDy2;
  }

*/

// Only true if the region is NOT to the right...
bool   BdryInterpolated2D::
qIsColinear(BdryPatch2D*, Vert*, int*)
{
  assert(0);
  return false;
}


BdryInterpolated2D::
~BdryInterpolated2D()
{
  // Destructor
  LiPoints.vClear();
  LpBCP2D.vClear();
}
