OCCT

// Created on: 1993-03-17
// Created by: Laurent BUCHARD
// Copyright (c) 1993-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#include <GeomAbs_SurfaceType.hxx>
#include <IntSurf_PntOn2S.hxx>
#include <Precision.hxx>
#include <StdFail_NotDone.hxx>
#include <NCollection_Array1.hxx>
#include <math_FunctionSetRoot.hxx>

//=================================================================================================
// function : IsSingular
// purpose  :    Returns TRUE if vectors theDU || theDV or if at least one
//            of them has null-magnitude.
//              theSqLinTol is square of linear tolerance.
//              theAngTol is angular tolerance.
//=================================================================================================

static bool IsSingular(const gp_Vec& theDU,
                       const gp_Vec& theDV,
                       const double  theSqLinTol,
                       const double  theAngTol)
{
  gp_Vec aDU(theDU), aDV(theDV);

  const double aSqMagnDU = aDU.SquareMagnitude(), aSqMagnDV = aDV.SquareMagnitude();

  if (aSqMagnDU < theSqLinTol)
    return true;

  aDU.Divide(sqrt(aSqMagnDU));

  if (aSqMagnDV < theSqLinTol)
    return true;

  aDV.Divide(sqrt(aSqMagnDV));

  // Here aDU and aDV vectors have magnitude 1.0.

  if (aDU.Crossed(aDV).SquareMagnitude() < theAngTol * theAngTol)
    return true;

  return false;
}

//=================================================================================================
// function : SingularProcessing
// purpose  :    Computes 2D-representation (in UV-coordinates) of
//            theTg3D vector on the surface in case when
//            theDU.Crossed(theDV).Magnitude() == 0.0. Stores result in
//            theTg2D variable.
//              theDU and theDV are vectors of 1st derivative
//            (with respect to U and V variables correspondingly).
//              If theIsTo3DTgCompute == TRUE then theTg3D has not been
//            defined yet (it should be computed).
//              theLinTol is SQUARE of the tolerance.
//
// Algorithm:
//              Condition
//                  Tg=theDU*theTg2D.X()+theDV*theTg2D.Y()
//            has to be satisfied strictly.
//              More over, vector Tg has to be NORMALIZED
//            (if theIsTo3DTgCompute == TRUE then new computed vector will
//            always have magnitude 1.0).
//=================================================================================================

static bool SingularProcessing(const gp_Vec& theDU,
                               const gp_Vec& theDV,
                               const bool    theIsTo3DTgCompute,
                               const double  theLinTol,
                               const double  theAngTol,
                               gp_Vec&       theTg3D,
                               gp_Vec2d&     theTg2D)
{
  // Attention: @ \sin theAngTol \approx theAngTol @ (for cross-product)

  // Really, vector theTg3D has to be normalized (if theIsTo3DTgCompute == FALSE).
  const double aSQTan = theTg3D.SquareMagnitude();

  const double aSqMagnDU = theDU.SquareMagnitude(), aSqMagnDV = theDV.SquareMagnitude();

  // There are some reasons of singularity

  // 1.
  if ((aSqMagnDU < theLinTol) && (aSqMagnDV < theLinTol))
  {
    // For future, this case can be processed as same as in case of
    // osculating surfaces (expanding in Taylor series). Here,
    // we return only.

    return false;
  }

  // 2.
  if (aSqMagnDU < theLinTol)
  {
    // In this case, theTg3D vector will be parallel with theDV.
    // Its true direction shall be precised later (the algorithm is
    // based on array of Walking-points).

    if (theIsTo3DTgCompute)
    {
      // theTg3D will be normalized. Its magnitude is
      const double aTgMagn = 1.0;

      const double aNorm = sqrt(aSqMagnDV);
      theTg3D            = theDV.Divided(aNorm);
      theTg2D.SetCoord(0.0, aTgMagn / aNorm);
    }
    else
    {
      // theTg3D is already defined.
      // Here we check only, if this tangent is parallel to theDV.

      if (theDV.Crossed(theTg3D).SquareMagnitude() < theAngTol * theAngTol * aSqMagnDV * aSQTan)
      {
        // theTg3D is parallel to theDV

        // Use sign "+" if theTg3D and theDV are codirectional
        // and sign "-" if opposite
        const double aDP = theTg3D.Dot(theDV);
        theTg2D.SetCoord(0.0, std::copysign(sqrt(aSQTan / aSqMagnDV), aDP));
      }
      else
      {
        // theTg3D is not parallel to theDV
        // It is abnormal

        return false;
      }
    }

    return true;
  }

  // 3.
  if (aSqMagnDV < theLinTol)
  {
    // In this case, theTg3D vector will be parallel with theDU.
    // Its true direction shall be precised later (the algorithm is
    // based on array of Walking-points).

    if (theIsTo3DTgCompute)
    {
      // theTg3D will be normalized. Its magnitude is
      const double aTgMagn = 1.0;

      const double aNorm = sqrt(aSqMagnDU);
      theTg3D            = theDU.Divided(aNorm);
      theTg2D.SetCoord(aTgMagn / aNorm, 0.0);
    }
    else
    {
      // theTg3D is already defined.
      // Here we check only, if this tangent is parallel to theDU.

      if (theDU.Crossed(theTg3D).SquareMagnitude() < theAngTol * theAngTol * aSqMagnDU * aSQTan)
      {
        // theTg3D is parallel to theDU

        // Use sign "+" if theTg3D and theDU are codirectional
        // and sign "-" if opposite
        const double aDP = theTg3D.Dot(theDU);
        theTg2D.SetCoord(std::copysign(sqrt(aSQTan / aSqMagnDU), aDP), 0.0);
      }
      else
      {
        // theTg3D is not parallel to theDU
        // It is abnormal

        return false;
      }
    }

    return true;
  }

  // 4. If aSqMagnDU > 0.0 && aSqMagnDV > 0.0 but theDV || theDU.

  const double aLenU = sqrt(aSqMagnDU), aLenV = sqrt(aSqMagnDV);

  // aLenSum > 0.0 definitely
  const double aLenSum = aLenU + aLenV;

  if (theDV.Dot(theDU) > 0.0)
  {
    // Vectors theDV and theDU are codirectional.

    if (theIsTo3DTgCompute)
    {
      theTg2D.SetCoord(1.0 / aLenSum, 1.0 / aLenSum);
      theTg3D = theDU * theTg2D.X() + theDV * theTg2D.Y();
    }
    else
    {
      // theTg3D is already defined.
      // Here we check only, if this tangent is parallel to theDU
      //(and theDV together).

      if (theDU.Crossed(theTg3D).SquareMagnitude() < theAngTol * theAngTol * aSqMagnDU * aSQTan)
      {
        // theTg3D is parallel to theDU

        const double aDP    = theTg3D.Dot(theDU);
        const double aLenTg = std::copysign(sqrt(aSQTan), aDP);
        theTg2D.SetCoord(aLenTg / aLenSum, aLenTg / aLenSum);
      }
      else
      {
        // theTg3D is not parallel to theDU
        // It is abnormal

        return false;
      }
    }
  }
  else
  {
    // Vectors theDV and theDU are opposite.

    if (theIsTo3DTgCompute)
    {
      // Here we chose theDU as direction of theTg3D.
      // True direction shall be precised later (the algorithm is
      // based on array of Walking-points).

      theTg2D.SetCoord(1.0 / aLenSum, -1.0 / aLenSum);
      theTg3D = theDU * theTg2D.X() + theDV * theTg2D.Y();
    }
    else
    {
      // theTg3D is already defined.
      // Here we check only, if this tangent is parallel to theDU
      //(and theDV together).

      if (theDU.Crossed(theTg3D).SquareMagnitude() < theAngTol * theAngTol * aSqMagnDU * aSQTan)
      {
        // theTg3D is parallel to theDU

        const double aDP    = theTg3D.Dot(theDU);
        const double aLenTg = std::copysign(sqrt(aSQTan), aDP);
        theTg2D.SetCoord(aLenTg / aLenSum, -aLenTg / aLenSum);
      }
      else
      {
        // theTg3D is not parallel to theDU
        // It is abnormal

        return false;
      }
    }
  }

  return true;
}

//=================================================================================================
// function : NonSingularProcessing
// purpose  :    Computes 2D-representation (in UV-coordinates) of
//            theTg3D vector on the surface in case when
//            theDU.Crossed(theDV).Magnitude() > 0.0. Stores result in
//            theTg2D variable.
//              theDU and theDV are vectors of 1st derivative
//            (with respect to U and V variables correspondingly).
//              theLinTol is SQUARE of the tolerance.
//
// Algorithm:
//              Condition
//                  Tg=theDU*theTg2D.X()+theDV*theTg2D.Y()
//            has to be satisfied strictly.
//              More over, vector Tg has always to be NORMALIZED.
//=================================================================================================

static bool NonSingularProcessing(const gp_Vec& theDU,
                                  const gp_Vec& theDV,
                                  const gp_Vec& theTg3D,
                                  const double  theLinTol,
                                  const double  theAngTol,
                                  gp_Vec2d&     theTg2D)
{
  const gp_Vec aNormal = theDU.Crossed(theDV);
  const double aSQMagn = aNormal.SquareMagnitude();

  if (IsSingular(theDU, theDV, theLinTol, theAngTol))
  {
    gp_Vec aTg(theTg3D);
    return SingularProcessing(theDU, theDV, false, theLinTol, theAngTol, aTg, theTg2D);
  }

  // If @\vec{T}=\vec{A}*U+\vec{B}*V@ then

  //  \left\{\begin{matrix}
  //  \vec{A} \times \vec{T} = (\vec{A} \times \vec{B})*V
  //  \vec{B} \times \vec{T} = (\vec{B} \times \vec{A})*U
  //  \end{matrix}\right.

  // From here, values of U and V can be found very easily
  //(if @\left \| \vec{A} \times \vec{B} \right \| > 0.0 @,
  // else it is singular case).

  const gp_Vec aTgU(theTg3D.Crossed(theDU)), aTgV(theTg3D.Crossed(theDV));
  const double aDeltaU = aTgV.SquareMagnitude() / aSQMagn;
  const double aDeltaV = aTgU.SquareMagnitude() / aSQMagn;

  theTg2D.SetCoord(std::copysign(sqrt(aDeltaU), aTgV.Dot(aNormal)),
                   -std::copysign(sqrt(aDeltaV), aTgU.Dot(aNormal)));

  return true;
}

//=================================================================================================

ApproxInt_ImpPrmSvSurfaces::ApproxInt_ImpPrmSvSurfaces(const TheISurface& ISurf,
                                                       const ThePSurface& PSurf)
    : MyIsTangent(false),
      MyHasBeenComputed(false),
      MyIsTangentbis(false),
      MyHasBeenComputedbis(false),
      MyImplicitFirst(true),
      MyZerImpFunc(PSurf, ISurf)
{
  SetUseSolver(true);
}

//=================================================================================================

ApproxInt_ImpPrmSvSurfaces::ApproxInt_ImpPrmSvSurfaces(const ThePSurface& PSurf,
                                                       const TheISurface& ISurf)
    : MyIsTangent(false),
      MyHasBeenComputed(false),
      MyIsTangentbis(false),
      MyHasBeenComputedbis(false),
      MyImplicitFirst(false),
      MyZerImpFunc(PSurf, ISurf)
{
  SetUseSolver(true);
}

//=================================================================================================

void ApproxInt_ImpPrmSvSurfaces::Pnt(const double u1,
                                     const double v1,
                                     const double u2,
                                     const double v2,
                                     gp_Pnt&      P)
{
  gp_Pnt   aP;
  gp_Vec   aT;
  gp_Vec2d aTS1, aTS2;
  double   tu1 = u1;
  double   tu2 = u2;
  double   tv1 = v1;
  double   tv2 = v2;
  this->Compute(tu1, tv1, tu2, tv2, aP, aT, aTS1, aTS2);
  P = MyPnt;
}

//=================================================================================================

bool ApproxInt_ImpPrmSvSurfaces::Tangency(const double u1,
                                          const double v1,
                                          const double u2,
                                          const double v2,
                                          gp_Vec&      T)
{
  gp_Pnt   aP;
  gp_Vec   aT;
  gp_Vec2d aTS1, aTS2;
  double   tu1 = u1;
  double   tu2 = u2;
  double   tv1 = v1;
  double   tv2 = v2;
  bool     t   = this->Compute(tu1, tv1, tu2, tv2, aP, aT, aTS1, aTS2);
  T            = MyTg;
  return (t);
}

//=================================================================================================

bool ApproxInt_ImpPrmSvSurfaces::TangencyOnSurf1(const double u1,
                                                 const double v1,
                                                 const double u2,
                                                 const double v2,
                                                 gp_Vec2d&    T)
{
  gp_Pnt   aP;
  gp_Vec   aT;
  gp_Vec2d aTS1, aTS2;
  double   tu1 = u1;
  double   tu2 = u2;
  double   tv1 = v1;
  double   tv2 = v2;
  bool     t   = this->Compute(tu1, tv1, tu2, tv2, aP, aT, aTS1, aTS2);
  T            = MyTguv1;
  return (t);
}

//=================================================================================================

bool ApproxInt_ImpPrmSvSurfaces::TangencyOnSurf2(const double u1,
                                                 const double v1,
                                                 const double u2,
                                                 const double v2,
                                                 gp_Vec2d&    T)
{
  gp_Pnt   aP;
  gp_Vec   aT;
  gp_Vec2d aTS1, aTS2;
  double   tu1 = u1;
  double   tu2 = u2;
  double   tv1 = v1;
  double   tv2 = v2;
  bool     t   = this->Compute(tu1, tv1, tu2, tv2, aP, aT, aTS1, aTS2);
  T            = MyTguv2;
  return (t);
}

//=================================================================================================
// function : Compute
// purpose  :    Computes point on curve, 3D and 2D-tangents of a curve and
//            parameters on the surfaces.
//=================================================================================================

bool ApproxInt_ImpPrmSvSurfaces::Compute(double&   u1,
                                         double&   v1,
                                         double&   u2,
                                         double&   v2,
                                         gp_Pnt&   P,
                                         gp_Vec&   Tg,
                                         gp_Vec2d& Tguv1,
                                         gp_Vec2d& Tguv2)
{
  const IntSurf_Quadric& aQSurf  = MyZerImpFunc.ISurface();
  const ThePSurface&     aPSurf  = MyZerImpFunc.PSurface();
  gp_Vec2d&              aQuadTg = MyImplicitFirst ? Tguv1 : Tguv2;
  gp_Vec2d&              aPrmTg  = MyImplicitFirst ? Tguv2 : Tguv1;

  // for square
  constexpr double aNullValue = Precision::Approximation() * Precision::Approximation(),
                   anAngTol   = Precision::Angular();

  double tu1 = u1;
  double tu2 = u2;
  double tv1 = v1;
  double tv2 = v2;

  if (MyHasBeenComputed)
  {
    if ((MyParOnS1.X() == u1) && (MyParOnS1.Y() == v1) && (MyParOnS2.X() == u2)
        && (MyParOnS2.Y() == v2))
    {
      return (MyIsTangent);
    }
    else if (MyHasBeenComputedbis == false)
    {
      MyTgbis              = MyTg;
      MyTguv1bis           = MyTguv1;
      MyTguv2bis           = MyTguv2;
      MyPntbis             = MyPnt;
      MyParOnS1bis         = MyParOnS1;
      MyParOnS2bis         = MyParOnS2;
      MyIsTangentbis       = MyIsTangent;
      MyHasBeenComputedbis = MyHasBeenComputed;
    }
  }

  if (MyHasBeenComputedbis)
  {
    if ((MyParOnS1bis.X() == u1) && (MyParOnS1bis.Y() == v1) && (MyParOnS2bis.X() == u2)
        && (MyParOnS2bis.Y() == v2))
    {

      gp_Vec   TV(MyTg);
      gp_Vec2d TV1(MyTguv1);
      gp_Vec2d TV2(MyTguv2);
      gp_Pnt   TP(MyPnt);
      gp_Pnt2d TP1(MyParOnS1);
      gp_Pnt2d TP2(MyParOnS2);
      bool     TB = MyIsTangent;

      MyTg        = MyTgbis;
      MyTguv1     = MyTguv1bis;
      MyTguv2     = MyTguv2bis;
      MyPnt       = MyPntbis;
      MyParOnS1   = MyParOnS1bis;
      MyParOnS2   = MyParOnS2bis;
      MyIsTangent = MyIsTangentbis;

      MyTgbis        = TV;
      MyTguv1bis     = TV1;
      MyTguv2bis     = TV2;
      MyPntbis       = TP;
      MyParOnS1bis   = TP1;
      MyParOnS2bis   = TP2;
      MyIsTangentbis = TB;

      return (MyIsTangent);
    }
  }

  math_Vector X(1, 2);
  math_Vector BornInf(1, 2), BornSup(1, 2), Tolerance(1, 2);
  //--- ThePSurfaceTool::GetResolution(aPSurf,Tolerance(1),Tolerance(2));
  Tolerance(1) = 1.0e-8;
  Tolerance(2) = 1.0e-8;
  double binfu, bsupu, binfv, bsupv;
  binfu               = ThePSurfaceTool::FirstUParameter(aPSurf);
  binfv               = ThePSurfaceTool::FirstVParameter(aPSurf);
  bsupu               = ThePSurfaceTool::LastUParameter(aPSurf);
  bsupv               = ThePSurfaceTool::LastVParameter(aPSurf);
  BornInf(1)          = binfu;
  BornSup(1)          = bsupu;
  BornInf(2)          = binfv;
  BornSup(2)          = bsupv;
  double TranslationU = 0., TranslationV = 0.;

  if (!FillInitialVectorOfSolution(u1,
                                   v1,
                                   u2,
                                   v2,
                                   binfu,
                                   bsupu,
                                   binfv,
                                   bsupv,
                                   X,
                                   TranslationU,
                                   TranslationV))
  {
    MyIsTangent = MyIsTangentbis = false;
    MyHasBeenComputed = MyHasBeenComputedbis = false;
    return (false);
  }

  bool   aRsnldIsDone = false;
  double PourTesterU  = X(1);
  double PourTesterV  = X(2);
  if (GetUseSolver())
  {
    math_FunctionSetRoot Rsnld(MyZerImpFunc);
    Rsnld.SetTolerance(Tolerance);
    Rsnld.Perform(MyZerImpFunc, X, BornInf, BornSup);
    aRsnldIsDone = Rsnld.IsDone();
    if (aRsnldIsDone)
      Rsnld.Root(X);
  }
  if (aRsnldIsDone || !GetUseSolver())
  {
    MyHasBeenComputed = true;

    double DistAvantApresU = std::abs(PourTesterU - X(1));
    double DistAvantApresV = std::abs(PourTesterV - X(2));

    MyPnt = P = ThePSurfaceTool::Value(aPSurf, X(1), X(2));

    if ((DistAvantApresV <= 0.001) && (DistAvantApresU <= 0.001))
    {
      gp_Vec aD1uPrm, aD1vPrm;
      gp_Vec aD1uQuad, aD1vQuad;

      if (MyImplicitFirst)
      {
        u2 = X(1) - TranslationU;
        v2 = X(2) - TranslationV;

        if (aQSurf.TypeQuadric() != GeomAbs_Plane)
        {
          while (u1 - tu1 > M_PI)
            u1 -= M_PI + M_PI;
          while (tu1 - u1 > M_PI)
            u1 += M_PI + M_PI;
        }

        MyParOnS1.SetCoord(tu1, tv1);
        MyParOnS2.SetCoord(tu2, tv2);

        gp_Pnt aP2;

        ThePSurfaceTool::D1(aPSurf, X(1), X(2), P, aD1uPrm, aD1vPrm);
        aQSurf.D1(u1, v1, aP2, aD1uQuad, aD1vQuad);

        // Middle-point of P-P2 segment
        P.BaryCenter(1.0, aP2, 1.0);
      }
      else
      {
        u1 = X(1) - TranslationU;
        v1 = X(2) - TranslationV;
        // aQSurf.Parameters(P, u2, v2);
        if (aQSurf.TypeQuadric() != GeomAbs_Plane)
        {
          while (u2 - tu2 > M_PI)
            u2 -= M_PI + M_PI;
          while (tu2 - u2 > M_PI)
            u2 += M_PI + M_PI;
        }

        MyParOnS1.SetCoord(tu1, tv1);
        MyParOnS2.SetCoord(tu2, tu2);

        gp_Pnt aP2;
        ThePSurfaceTool::D1(aPSurf, X(1), X(2), P, aD1uPrm, aD1vPrm);

        aQSurf.D1(u2, v2, aP2, aD1uQuad, aD1vQuad);

        // Middle-point of P-P2 segment
        P.BaryCenter(1.0, aP2, 1.0);
      }

      MyPnt = P;

      // Normals to the surfaces
      gp_Vec aNormalPrm(aD1uPrm.Crossed(aD1vPrm)), aNormalImp(aQSurf.Normale(MyPnt));

      const double aSQMagnPrm = aNormalPrm.SquareMagnitude(),
                   aSQMagnImp = aNormalImp.SquareMagnitude();

      bool isPrmSingular = false, isImpSingular = false;

      if (IsSingular(aD1uPrm, aD1vPrm, aNullValue, anAngTol))
      {
        isPrmSingular = true;
        if (!SingularProcessing(aD1uPrm, aD1vPrm, true, aNullValue, anAngTol, Tg, aPrmTg))
        {
          MyIsTangent       = false;
          MyHasBeenComputed = MyHasBeenComputedbis = false;
          return false;
        }

        MyTg = Tg;
      }
      else
      {
        aNormalPrm.Divide(sqrt(aSQMagnPrm));
      }

      // Analogically for implicit surface
      if (aSQMagnImp < aNullValue)
      {
        isImpSingular = true;

        if (!SingularProcessing(aD1uQuad,
                                aD1vQuad,
                                !isPrmSingular,
                                aNullValue,
                                anAngTol,
                                Tg,
                                aQuadTg))
        {
          MyIsTangent       = false;
          MyHasBeenComputed = MyHasBeenComputedbis = false;
          return false;
        }

        MyTg = Tg;
      }
      else
      {
        aNormalImp.Divide(sqrt(aSQMagnImp));
      }

      if (isImpSingular && isPrmSingular)
      {
        // All is OK. All abnormal cases were processed above.

        MyTguv1 = Tguv1;
        MyTguv2 = Tguv2;

        MyIsTangent = true;
        return MyIsTangent;
      }
      else if (!(isImpSingular || isPrmSingular))
      {
        // Processing pure non-singular case
        //(3D- and 2D-tangents are still not defined)

        // Ask to pay attention to the fact that here
        // aNormalImp and aNormalPrm are normalized.
        // Therefore, @ \left \| \vec{Tg} \right \| = 0.0 @
        // if and only if (aNormalImp || aNormalPrm).
        Tg = aNormalImp.Crossed(aNormalPrm);
      }

      const double aSQMagnTg = Tg.SquareMagnitude();

      if (aSQMagnTg < aNullValue)
      {
        MyIsTangent       = false;
        MyHasBeenComputed = MyHasBeenComputedbis = false;
        return false;
      }

      // Normalize Tg vector
      Tg.Divide(sqrt(aSQMagnTg));
      MyTg = Tg;

      if (!isPrmSingular)
      {
        // If isPrmSingular==TRUE then aPrmTg has already been computed.

        if (!NonSingularProcessing(aD1uPrm, aD1vPrm, Tg, aNullValue, anAngTol, aPrmTg))
        {
          MyIsTangent       = false;
          MyHasBeenComputed = MyHasBeenComputedbis = false;
          return false;
        }
      }

      if (!isImpSingular)
      {
        // If isImpSingular==TRUE then aQuadTg has already been computed.

        if (!NonSingularProcessing(aD1uQuad, aD1vQuad, Tg, aNullValue, anAngTol, aQuadTg))
        {
          MyIsTangent       = false;
          MyHasBeenComputed = MyHasBeenComputedbis = false;
          return false;
        }
      }

      MyTguv1 = Tguv1;
      MyTguv2 = Tguv2;

      MyIsTangent = true;

#ifdef OCCT_DEBUG
      // cout << "+++++++++++++++++  ApproxInt_ImpPrmSvSurfaces::Compute(...)  ++++++++++" << endl;
      // printf( "P2d_1(%+10.20f, %+10.20f); P2d_2(%+10.20f, %+10.20f)\n"
      //         "P(%+10.20f, %+10.20f, %+10.20f);\n"
      //         "Tg = {%+10.20f, %+10.20f, %+10.20f};\n"
      //         "Tguv1 = {%+10.20f, %+10.20f};\n"
      //         "Tguv2 = {%+10.20f, %+10.20f}\n",
      //         u1, v1, u2, v2,
      //         P.X(), P.Y(), P.Z(),
      //         Tg.X(), Tg.Y(), Tg.Z(),
      //         Tguv1.X(), Tguv1.Y(), Tguv2.X(), Tguv2.Y());
      // cout << "-----------------------------------------------------------------------" << endl;
#endif

      return true;
    }
    else
    {
      //-- cout<<" ApproxInt_ImpImpSvSurfaces.gxx : Distance apres recadrage Trop Grande "<<endl;

      MyIsTangent       = false;
      MyHasBeenComputed = MyHasBeenComputedbis = false;
      return false;
    }
  }
  else
  {
    MyIsTangent       = false;
    MyHasBeenComputed = MyHasBeenComputedbis = false;
    return false;
  }
}

//=================================================================================================
// function : SeekPoint
// purpose  :    Computes point on curve and
//            parameters on the surfaces.
//=================================================================================================

bool ApproxInt_ImpPrmSvSurfaces::SeekPoint(const double     u1,
                                           const double     v1,
                                           const double     u2,
                                           const double     v2,
                                           IntSurf_PntOn2S& Point)
{
  const IntSurf_Quadric& aQSurf = MyZerImpFunc.ISurface();
  const ThePSurface&     aPSurf = MyZerImpFunc.PSurface();

  math_Vector X(1, 2);
  math_Vector BornInf(1, 2), BornSup(1, 2), Tolerance(1, 2);
  //--- ThePSurfaceTool::GetResolution(aPSurf,Tolerance(1),Tolerance(2));
  Tolerance(1) = 1.0e-8;
  Tolerance(2) = 1.0e-8;
  double binfu, bsupu, binfv, bsupv;
  binfu               = ThePSurfaceTool::FirstUParameter(aPSurf);
  binfv               = ThePSurfaceTool::FirstVParameter(aPSurf);
  bsupu               = ThePSurfaceTool::LastUParameter(aPSurf);
  bsupv               = ThePSurfaceTool::LastVParameter(aPSurf);
  BornInf(1)          = binfu;
  BornSup(1)          = bsupu;
  BornInf(2)          = binfv;
  BornSup(2)          = bsupv;
  double TranslationU = 0., TranslationV = 0.;

  if (!FillInitialVectorOfSolution(u1,
                                   v1,
                                   u2,
                                   v2,
                                   binfu,
                                   bsupu,
                                   binfv,
                                   bsupv,
                                   X,
                                   TranslationU,
                                   TranslationV))
    return false;

  double NewU1, NewV1, NewU2, NewV2;

  math_FunctionSetRoot Rsnld(MyZerImpFunc);
  Rsnld.SetTolerance(Tolerance);
  Rsnld.Perform(MyZerImpFunc, X, BornInf, BornSup);
  if (Rsnld.IsDone())
  {
    MyHasBeenComputed = true;
    Rsnld.Root(X);

    MyPnt = ThePSurfaceTool::Value(aPSurf, X(1), X(2));

    if (MyImplicitFirst)
    {
      NewU2 = X(1) - TranslationU;
      NewV2 = X(2) - TranslationV;

      aQSurf.Parameters(MyPnt, NewU1, NewV1);
      // adjust U
      if (aQSurf.TypeQuadric() != GeomAbs_Plane)
      {
        double sign = (NewU1 > u1) ? -1 : 1;
        while (std::abs(u1 - NewU1) > M_PI)
          NewU1 += sign * (M_PI + M_PI);
      }
    }
    else
    {
      NewU1 = X(1) - TranslationU;
      NewV1 = X(2) - TranslationV;

      aQSurf.Parameters(MyPnt, NewU2, NewV2);
      // adjust U
      if (aQSurf.TypeQuadric() != GeomAbs_Plane)
      {
        double sign = (NewU2 > u2) ? -1 : 1;
        while (std::abs(u2 - NewU2) > M_PI)
          NewU2 += sign * (M_PI + M_PI);
      }
    }
  }
  else
    return false;

  Point.SetValue(MyPnt, NewU1, NewV1, NewU2, NewV2);
  return true;
}

//=================================================================================================

bool ApproxInt_ImpPrmSvSurfaces::FillInitialVectorOfSolution(const double u1,
                                                             const double v1,
                                                             const double u2,
                                                             const double v2,
                                                             const double binfu,
                                                             const double bsupu,
                                                             const double binfv,
                                                             const double bsupv,
                                                             math_Vector& X,
                                                             double&      TranslationU,
                                                             double&      TranslationV)
{
  const ThePSurface& aPSurf = MyZerImpFunc.PSurface();

  math_Vector F(1, 1);

  TranslationU = 0.0;
  TranslationV = 0.0;

  if (MyImplicitFirst)
  {
    if (u2 < binfu - 0.0000000001)
    {
      if (ThePSurfaceTool::IsUPeriodic(aPSurf))
      {
        double d = ThePSurfaceTool::UPeriod(aPSurf);
        do
        {
          TranslationU += d;
        } while (u2 + TranslationU < binfu);
      }
      else
        return (false);
    }
    else if (u2 > bsupu + 0.0000000001)
    {
      if (ThePSurfaceTool::IsUPeriodic(aPSurf))
      {
        double d = ThePSurfaceTool::UPeriod(aPSurf);
        do
        {
          TranslationU -= d;
        } while (u2 + TranslationU > bsupu);
      }
      else
        return (false);
    }
    if (v2 < binfv - 0.0000000001)
    {
      if (ThePSurfaceTool::IsVPeriodic(aPSurf))
      {
        double d = ThePSurfaceTool::VPeriod(aPSurf);
        do
        {
          TranslationV += d;
        } while (v2 + TranslationV < binfv);
      }
      else
        return (false);
    }
    else if (v2 > bsupv + 0.0000000001)
    {
      if (ThePSurfaceTool::IsVPeriodic(aPSurf))
      {
        double d = ThePSurfaceTool::VPeriod(aPSurf);
        do
        {
          TranslationV -= d;
        } while (v2 + TranslationV > bsupv);
      }
      else
        return (false);
    }
    X(1) = u2 + TranslationU;
    X(2) = v2 + TranslationV;
  }
  else
  {
    if (u1 < binfu - 0.0000000001)
    {
      if (ThePSurfaceTool::IsUPeriodic(aPSurf))
      {
        double d = ThePSurfaceTool::UPeriod(aPSurf);
        do
        {
          TranslationU += d;
        } while (u1 + TranslationU < binfu);
      }
      else
        return (false);
    }
    else if (u1 > bsupu + 0.0000000001)
    {
      if (ThePSurfaceTool::IsUPeriodic(aPSurf))
      {
        double d = ThePSurfaceTool::UPeriod(aPSurf);
        do
        {
          TranslationU -= d;
        } while (u1 + TranslationU > bsupu);
      }
      else
        return (false);
    }
    if (v1 < binfv - 0.0000000001)
    {
      if (ThePSurfaceTool::IsVPeriodic(aPSurf))
      {
        double d = ThePSurfaceTool::VPeriod(aPSurf);
        do
        {
          TranslationV += d;
        } while (v1 + TranslationV < binfv);
      }
      else
        return (false);
    }
    else if (v1 > bsupv + 0.0000000001)
    {
      if (ThePSurfaceTool::IsVPeriodic(aPSurf))
      {
        double d = ThePSurfaceTool::VPeriod(aPSurf);
        do
        {
          TranslationV -= d;
        } while (v1 + TranslationV > bsupv);
      }
      else
        return (false);
    }
    X(1) = u1 + TranslationU;
    X(2) = v1 + TranslationV;
  }

  //=================================================================================================

  // Make a small step from boundaries in order to avoid
  // finding "outboundaried" solution (Rsnld -> NotDone).
  if (GetUseSolver())
  {
    double du = std::max(Precision::Confusion(),
                         ThePSurfaceTool::UResolution(aPSurf, Precision::Confusion()));
    double dv = std::max(Precision::Confusion(),
                         ThePSurfaceTool::VResolution(aPSurf, Precision::Confusion()));
    if (X(1) - 0.0000000001 <= binfu)
      X(1) = X(1) + du;
    if (X(1) + 0.0000000001 >= bsupu)
      X(1) = X(1) - du;
    if (X(2) - 0.0000000001 <= binfv)
      X(2) = X(2) + dv;
    if (X(2) + 0.0000000001 >= bsupv)
      X(2) = X(2) - dv;
  }

  return true;
}