[occt.git] / src / Extrema / Extrema_ExtPElC.cxx

// Copyright (c) 1995-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 <ElCLib.hxx>
#include <Extrema_ExtPElC.hxx>
#include <Extrema_POnCurv.hxx>
#include <gp_Circ.hxx>
#include <gp_Elips.hxx>
#include <gp_Hypr.hxx>
#include <gp_Lin.hxx>
#include <gp_Parab.hxx>
#include <gp_Pnt.hxx>
#include <math_DirectPolynomialRoots.hxx>
#include <math_TrigonometricFunctionRoots.hxx>
#include <Precision.hxx>
#include <Standard_NotImplemented.hxx>
#include <Standard_OutOfRange.hxx>
#include <StdFail_NotDone.hxx>

//=======================================================================
//function : Extrema_ExtPElC
//purpose  : 
//=======================================================================
Extrema_ExtPElC::Extrema_ExtPElC()
{
  myDone = Standard_False;
  myNbExt = 0;

  for (Standard_Integer i = 0; i < 4; i++)
  {
    mySqDist[i] = RealLast();
    myIsMin[i] = Standard_False;
  }
}

//=======================================================================
//function : Extrema_ExtPElC
//purpose  : 
//=======================================================================
Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P, 
				  const gp_Lin&       L,
				  const Standard_Real Tol,
				  const Standard_Real Uinf, 
				  const Standard_Real Usup)
{
  Perform(P, L, Tol, Uinf, Usup);
}
//=======================================================================
//function : Perform
//purpose  : 
//=======================================================================
void Extrema_ExtPElC::Perform(const gp_Pnt&       P, 
			      const gp_Lin&       L,
			      const Standard_Real Tol,
			      const Standard_Real Uinf, 
			      const Standard_Real Usup)
{
  myDone = Standard_False;
  myNbExt = 0;
  gp_Vec V1 (L.Direction());
  gp_Pnt OR = L.Location();
  gp_Vec V(OR, P);
  Standard_Real Mydist = V1.Dot(V);
  if ((Mydist >= Uinf-Tol) && 
      (Mydist <= Usup+Tol)){ 
    
    gp_Pnt MyP = OR.Translated(Mydist*V1);
    Extrema_POnCurv MyPOnCurve(Mydist, MyP);
    mySqDist[0] = P.SquareDistance(MyP);
    myPoint[0] = MyPOnCurve;
    myIsMin[0] = Standard_True;
    myNbExt = 1;
    myDone = Standard_True;
  }

}


Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P,
				  const gp_Circ&      C,
				  const Standard_Real Tol,
				  const Standard_Real Uinf,
				  const Standard_Real Usup)
{
  Perform(P, C, Tol, Uinf, Usup);
}

void Extrema_ExtPElC::Perform(const gp_Pnt&       P, 
			      const gp_Circ&      C,
			      const Standard_Real Tol,
			      const Standard_Real Uinf,
			      const Standard_Real Usup)
/*-----------------------------------------------------------------------------
Function:
  Find values of parameter u such as:
   - dist(P,C(u)) pass by an extrema,
   - Uinf <= u <= Usup.

Method:
  Pass 3 stages:
  1- Projection of point P in the plane of the circle,
  2- Calculation of u solutions in [0.,2.*M_PI]:
      Let Pp, the projected point and 
           O, the center of the circle;
     2 cases:
     - if Pp is mixed with 0, there is an infinite number of solutions;
       IsDone() renvoie Standard_False.
     - otherwise, 2 points are solutions for the complete circle:
       . Us1 = angle(OPp,OX) corresponds to the minimum,
       . let Us2 = ( Us1 + M_PI if Us1 < M_PI,
                    ( Us1 - M_PI otherwise;
         Us2 corresponds to the maximum.
  3- Calculate the extrema in [Uinf,Usup].
-----------------------------------------------------------------------------*/
{
  myDone = Standard_False;
  myNbExt = 0;

// 1- Projection of the point P in the plane of circle -> Pp ...

  gp_Pnt O = C.Location();
  gp_Vec Axe (C.Axis().Direction());
  gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
  gp_Pnt Pp = P.Translated(Trsl);

// 2- Calculate u solutions in [0.,2.*PI] ...

  gp_Vec OPp (O,Pp);
  if (OPp.Magnitude() < Tol) { return; }
  Standard_Real Usol[2];
  Usol[0] = C.XAxis().Direction().AngleWithRef(OPp,Axe); // -M_PI<U1<M_PI

  const Standard_Real aAngTol = Precision::Angular();
  if ( Usol[0] + M_PI < aAngTol )
    Usol[0] = -M_PI;
  else if ( Usol[0] - M_PI > -aAngTol )
    Usol[0] = M_PI; 

  Usol[1] = Usol[0] + M_PI;

  Standard_Real myuinf = Uinf;
  //Standard_Real TolU = Tol*C.Radius();
  Standard_Real TolU, aR;
  aR=C.Radius();
  TolU=Precision::Infinite();
  if (aR > gp::Resolution()) {
    TolU= Tol/aR;
  }
  //
  ElCLib::AdjustPeriodic(Uinf, Uinf+2*M_PI, TolU, myuinf, Usol[0]);
  ElCLib::AdjustPeriodic(Uinf, Uinf+2*M_PI, TolU, myuinf, Usol[1]);
  if (((Usol[0]-2*M_PI-Uinf) < TolU) && ((Usol[0]-2*M_PI-Uinf) > -TolU)) Usol[0] = Uinf;
  if (((Usol[1]-2*M_PI-Uinf) < TolU) && ((Usol[1]-2*M_PI-Uinf) > -TolU)) Usol[1] = Uinf;


// 3- Calculate extrema in [Umin,Umax] ...

  gp_Pnt Cu;
  Standard_Real Us;
  for (Standard_Integer NoSol = 0; NoSol <= 1; NoSol++) {
    Us = Usol[NoSol];
    if (((Uinf-Us) < TolU) && ((Us-Usup) < TolU)) {
      Cu = ElCLib::Value(Us,C);
      mySqDist[myNbExt] = Cu.SquareDistance(P);
      myIsMin[myNbExt] = (NoSol == 0);
      myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
      myNbExt++;
    }
  }
  myDone = Standard_True;
}
//=============================================================================

Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P,
				  const gp_Elips&     C,
				  const Standard_Real Tol,
				  const Standard_Real Uinf, 
				  const Standard_Real Usup)
{
  Perform(P, C, Tol, Uinf, Usup);
}


void Extrema_ExtPElC::Perform (const gp_Pnt&       P, 
			       const gp_Elips&     C,
			       const Standard_Real Tol,
			       const Standard_Real Uinf, 
			       const Standard_Real Usup)
/*-----------------------------------------------------------------------------
Function:
  Find values of parameter u so that:
   - dist(P,C(u)) passes by an extremum,
   - Uinf <= u <= Usup.

Method:
  Takes 2 steps:
  1- Projection of point P in the plane of the ellipse,
  2- Calculation of the solutions:
     Let Pp, the projected point; find values u so that:
       (C(u)-Pp).C'(u) = 0. (1)
     Let Cos = cos(u) and Sin = sin(u),
          C(u) = (A*Cos,B*Sin) and Pp = (X,Y);
     Then, (1) <=> (A*Cos-X,B*Sin-Y).(-A*Sin,B*Cos) = 0.
                    (B**2-A**2)*Cos*Sin - B*Y*Cos + A*X*Sin = 0.
     Use algorithm math_TrigonometricFunctionRoots to solve this equation.
-----------------------------------------------------------------------------*/
{
  myDone = Standard_False;
  myNbExt = 0;

// 1- Projection of point P in the plane of the ellipse -> Pp ...

  gp_Pnt O = C.Location();
  gp_Vec Axe (C.Axis().Direction());
  gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
  gp_Pnt Pp = P.Translated(Trsl);

// 2- Calculation of solutions ...

  Standard_Integer NoSol, NbSol;
  Standard_Real A = C.MajorRadius();
  Standard_Real B = C.MinorRadius();
  gp_Vec OPp (O,Pp);
  Standard_Real OPpMagn = OPp.Magnitude();
  if (OPpMagn < Tol) { if (Abs(A-B) < Tol) { return; } }
  Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
  Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
  //  Standard_Real Y = Sqrt(OPpMagn*OPpMagn-X*X);

  Standard_Real ko2 = (B*B-A*A)/2., ko3 = -B*Y, ko4 = A*X;
  if(Abs(ko3) < 1.e-16*Max(Abs(ko2), Abs(ko3))) ko3 = 0.0;

//  math_TrigonometricFunctionRoots Sol(0.,(B*B-A*A)/2.,-B*Y,A*X,0.,Uinf,Usup);
  math_TrigonometricFunctionRoots Sol(0.,ko2, ko3, ko4, 0.,Uinf,Usup);

  if (!Sol.IsDone()) { return; }
  gp_Pnt Cu;
  Standard_Real Us;
  NbSol = Sol.NbSolutions();
  for (NoSol = 1; NoSol <= NbSol; NoSol++) {
    Us = Sol.Value(NoSol);
    Cu = ElCLib::Value(Us,C);
    mySqDist[myNbExt] = Cu.SquareDistance(P);
    myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
    Cu = ElCLib::Value(Us + 0.1, C);
    myIsMin[myNbExt] = mySqDist[myNbExt] < Cu.SquareDistance(P);
    myNbExt++;
  }
  myDone = Standard_True;
}
//=============================================================================

Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P, 
				  const gp_Hypr&      C,
				  const Standard_Real Tol,
				  const Standard_Real Uinf,
				  const Standard_Real Usup)
{
  Perform(P, C, Tol, Uinf, Usup);
}


void Extrema_ExtPElC::Perform(const gp_Pnt&       P, 
			      const gp_Hypr&      C,
			      const Standard_Real Tol,
			      const Standard_Real Uinf,
			      const Standard_Real Usup)
/*-----------------------------------------------------------------------------
Function:
  Find values of parameter u so that:
   - dist(P,C(u)) passes by an extremum,
   - Uinf <= u <= Usup.

Method:
  Takes 2 steps:
  1- Projection of point P in the plane of the hyperbola,
  2- Calculation of solutions:
      Let Pp, le point projete; on recherche les valeurs u telles que:
       (C(u)-Pp).C'(u) = 0. (1)
      Let R and r be the radiuses of the hyperbola,
          Chu = Cosh(u) and Shu = Sinh(u),
          C(u) = (R*Chu,r*Shu) and Pp = (X,Y);
     Then, (1) <=> (R*Chu-X,r*Shu-Y).(R*Shu,r*Chu) = 0.
                    (R**2+r**2)*Chu*Shu - X*R*Shu - Y*r*Chu = 0. (2)
     Let v = e**u;
     Then, by using Chu = (e**u+e**(-u))/2. and Sh = (e**u-e**(-u)))/2.
            (2) <=> ((R**2+r**2)/4.) * (v**2-v**(-2)) -
	            ((X*R+Y*r)/2.) * v +
	            ((X*R-Y*r)/2.) * v**(-1) = 0.
      (2)* v**2	<=> ((R**2+r**2)/4.) * v**4 -
                     ((X*R+Y*r)/2.)  * v**3 +
		     ((X*R-Y*r)/2.)  * v    -
		    ((R**2+r**2)/4.)        = 0.
  Use algorithm math_DirectPolynomialRoots to solve this equation by v.
-----------------------------------------------------------------------------*/
{
  myDone = Standard_False;
  myNbExt = 0;

// 1- Projection of point P in the plane of hyperbola -> Pp ...

  gp_Pnt O = C.Location();
  gp_Vec Axe (C.Axis().Direction());
  gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
  gp_Pnt Pp = P.Translated(Trsl);

// 2- Calculation of solutions ...

  Standard_Real Tol2 = Tol * Tol;
  Standard_Real R = C.MajorRadius();
  Standard_Real r = C.MinorRadius();
  gp_Vec OPp (O,Pp);
  Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
  Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));

  Standard_Real C1 = (R*R+r*r)/4.;
  math_DirectPolynomialRoots Sol(C1,-(X*R+Y*r)/2.,0.,(X*R-Y*r)/2.,-C1);
  if (!Sol.IsDone()) { return; }
  gp_Pnt Cu;
  Standard_Real Us, Vs;
  Standard_Integer NbSol = Sol.NbSolutions();
  Standard_Boolean DejaEnr;
  Standard_Integer NoExt;
  gp_Pnt TbExt[4];
  for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
    Vs = Sol.Value(NoSol);
    if (Vs > 0.) {
      Us = Log(Vs);
      if ((Us >= Uinf) && (Us <= Usup)) {
	Cu = ElCLib::Value(Us,C);
	DejaEnr = Standard_False;
	for (NoExt = 0; NoExt < myNbExt; NoExt++) {
	  if (TbExt[NoExt].SquareDistance(Cu) < Tol2) {
	    DejaEnr = Standard_True;
	    break;
	  }
	}
	if (!DejaEnr) {
	  TbExt[myNbExt] = Cu;
	  mySqDist[myNbExt] = Cu.SquareDistance(P);
	  myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
	  myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
	  myNbExt++;
	}
      } // if ((Us >= Uinf) && (Us <= Usup))
    } // if (Vs > 0.)
  } // for (Standard_Integer NoSol = 1; ...
  myDone = Standard_True;
}
//=============================================================================

Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P,
				  const gp_Parab&     C,
				  const Standard_Real Tol,
				  const Standard_Real Uinf,
				  const Standard_Real Usup)
{
  Perform(P, C, Tol, Uinf, Usup);
}


void Extrema_ExtPElC::Perform(const gp_Pnt&       P, 
			      const gp_Parab&     C,
//			      const Standard_Real Tol,
			      const Standard_Real ,
			      const Standard_Real Uinf,
			      const Standard_Real Usup)
/*-----------------------------------------------------------------------------
Function:
  Find values of parameter u so that:
   - dist(P,C(u)) pass by an extremum,
   - Uinf <= u <= Usup.

Method:
  Takes 2 steps:
  1- Projection of point P in the plane of the parabola,
  2- Calculation of solutions:
     Let Pp, the projected point; find values u so that:
       (C(u)-Pp).C'(u) = 0. (1)
     Let F the focus of the parabola,
          C(u) = ((u*u)/(4.*F),u) and Pp = (X,Y);
     Alors, (1) <=> ((u*u)/(4.*F)-X,u-Y).(u/(2.*F),1) = 0.
                    (1./(4.*F)) * U**3 + (2.*F-X) * U - 2*F*Y = 0.
    Use algorithm math_DirectPolynomialRoots to solve this equation by U.
-----------------------------------------------------------------------------*/
{
  myDone = Standard_False;
  myNbExt = 0;

// 1- Projection of point P in the plane of the parabola -> Pp ...

  gp_Pnt O = C.Location();
  gp_Vec Axe (C.Axis().Direction());
  gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
  gp_Pnt Pp = P.Translated(Trsl);

// 2- Calculation of solutions ...

  Standard_Real F = C.Focal();
  gp_Vec OPp (O,Pp);
  Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
//  Standard_Real Y = Sqrt(OPpMagn*OPpMagn-X*X);
  Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
  math_DirectPolynomialRoots Sol(1./(4.*F),0.,2.*F-X,-2.*F*Y);
  if (!Sol.IsDone()) { return; }
  gp_Pnt Cu;
  Standard_Real Us;
  Standard_Integer NbSol = Sol.NbSolutions();
  Standard_Boolean DejaEnr;
  Standard_Integer NoExt;
  gp_Pnt TbExt[3];
  for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
    Us = Sol.Value(NoSol);
    if ((Us >= Uinf) && (Us <= Usup)) {
      Cu = ElCLib::Value(Us,C);
      DejaEnr = Standard_False;
      for (NoExt = 0; NoExt < myNbExt; NoExt++) {
        if (TbExt[NoExt].SquareDistance(Cu) < Precision::SquareConfusion()) {
	  DejaEnr = Standard_True;
	  break;
	}
      }
      if (!DejaEnr) {
	TbExt[myNbExt] = Cu;
	mySqDist[myNbExt] = Cu.SquareDistance(P);
	myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
	myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
	myNbExt++;
      }
    } // if ((Us >= Uinf) && (Us <= Usup))
  } // for (Standard_Integer NoSol = 1; ...
  myDone = Standard_True;
}
//=============================================================================

Standard_Boolean Extrema_ExtPElC::IsDone () const { return myDone; }
//=============================================================================

Standard_Integer Extrema_ExtPElC::NbExt () const
{
  if (!IsDone()) { throw StdFail_NotDone(); }
  return myNbExt;
}
//=============================================================================

Standard_Real Extrema_ExtPElC::SquareDistance (const Standard_Integer N) const
{
  if ((N < 1) || (N > NbExt())) { throw Standard_OutOfRange(); }
  return mySqDist[N-1];
}
//=============================================================================

Standard_Boolean Extrema_ExtPElC::IsMin (const Standard_Integer N) const
{
  if ((N < 1) || (N > NbExt())) { throw Standard_OutOfRange(); }
  return myIsMin[N-1];
}
//=============================================================================

const Extrema_POnCurv& Extrema_ExtPElC::Point (const Standard_Integer N) const
{
  if ((N < 1) || (N > NbExt())) { throw Standard_OutOfRange(); }
  return myPoint[N-1];
}
//=============================================================================
Commit	Line	Data
b311480e	1	// Copyright (c) 1995-1999 Matra Datavision
973c2be1	2	// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e	3	//
973c2be1	4	// This file is part of Open CASCADE Technology software library.
b311480e	5	//
d5f74e42	6	// This library is free software; you can redistribute it and/or modify it under
d5f74e42	7	// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1	8	// by the Free Software Foundation, with special exception defined in the file
	9	// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
	10	// distribution for complete text of the license and disclaimer of any warranty.
b311480e	11	//
973c2be1	12	// Alternatively, this file may be used under the terms of Open CASCADE
973c2be1	13	// commercial license or contractual agreement.
b311480e	14
42cf5bc1	15
	16	#include <ElCLib.hxx>
	17	#include <Extrema_ExtPElC.hxx>
	18	#include <Extrema_POnCurv.hxx>
	19	#include <gp_Circ.hxx>
	20	#include <gp_Elips.hxx>
	21	#include <gp_Hypr.hxx>
	22	#include <gp_Lin.hxx>
	23	#include <gp_Parab.hxx>
	24	#include <gp_Pnt.hxx>
7fd59977	25	#include <math_DirectPolynomialRoots.hxx>
7fd59977	26	#include <math_TrigonometricFunctionRoots.hxx>
7fd59977	27	#include <Precision.hxx>
42cf5bc1	28	#include <Standard_NotImplemented.hxx>
	29	#include <Standard_OutOfRange.hxx>
	30	#include <StdFail_NotDone.hxx>
7fd59977	31
42ff8f5b	32	//=======================================================================
	33	//function : Extrema_ExtPElC
	34	//purpose :
	35	//=======================================================================
638ad7f3	36	Extrema_ExtPElC::Extrema_ExtPElC()
	37	{
	38	myDone = Standard_False;
	39	myNbExt = 0;
	40
	41	for (Standard_Integer i = 0; i < 4; i++)
	42	{
	43	mySqDist[i] = RealLast();
	44	myIsMin[i] = Standard_False;
	45	}
	46	}
7fd59977	47
42ff8f5b	48	//=======================================================================
	49	//function : Extrema_ExtPElC
	50	//purpose :
	51	//=======================================================================
7fd59977	52	Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
	53	const gp_Lin& L,
	54	const Standard_Real Tol,
	55	const Standard_Real Uinf,
	56	const Standard_Real Usup)
	57	{
	58	Perform(P, L, Tol, Uinf, Usup);
	59	}
42ff8f5b	60	//=======================================================================
	61	//function : Perform
	62	//purpose :
	63	//=======================================================================
7fd59977	64	void Extrema_ExtPElC::Perform(const gp_Pnt& P,
	65	const gp_Lin& L,
	66	const Standard_Real Tol,
	67	const Standard_Real Uinf,
	68	const Standard_Real Usup)
	69	{
	70	myDone = Standard_False;
	71	myNbExt = 0;
42ff8f5b	72	gp_Vec V1 (L.Direction());
7fd59977	73	gp_Pnt OR = L.Location();
	74	gp_Vec V(OR, P);
	75	Standard_Real Mydist = V1.Dot(V);
	76	if ((Mydist >= Uinf-Tol) &&
	77	(Mydist <= Usup+Tol)){
	78
	79	gp_Pnt MyP = OR.Translated(Mydist*V1);
	80	Extrema_POnCurv MyPOnCurve(Mydist, MyP);
	81	mySqDist[0] = P.SquareDistance(MyP);
	82	myPoint[0] = MyPOnCurve;
	83	myIsMin[0] = Standard_True;
	84	myNbExt = 1;
	85	myDone = Standard_True;
	86	}
	87
	88	}
	89
	90
	91
	92	Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
	93	const gp_Circ& C,
	94	const Standard_Real Tol,
	95	const Standard_Real Uinf,
	96	const Standard_Real Usup)
	97	{
	98	Perform(P, C, Tol, Uinf, Usup);
	99	}
	100
	101	void Extrema_ExtPElC::Perform(const gp_Pnt& P,
	102	const gp_Circ& C,
	103	const Standard_Real Tol,
	104	const Standard_Real Uinf,
	105	const Standard_Real Usup)
	106	/*-----------------------------------------------------------------------------
0d969553 Y	107	Function:
	108	Find values of parameter u such as:
	109	- dist(P,C(u)) pass by an extrema,
7fd59977	110	- Uinf <= u <= Usup.
7fd59977	111
0d969553 Y	112	Method:
	113	Pass 3 stages:
	114	1- Projection of point P in the plane of the circle,
c6541a0c	115	2- Calculation of u solutions in [0.,2.*M_PI]:
0d969553 Y	116	Let Pp, the projected point and
	117	O, the center of the circle;
	118	2 cases:
	119	- if Pp is mixed with 0, there is an infinite number of solutions;
7fd59977	120	IsDone() renvoie Standard_False.
0d969553 Y	121	- otherwise, 2 points are solutions for the complete circle:
0d969553 Y	122	. Us1 = angle(OPp,OX) corresponds to the minimum,
c6541a0c D	123	. let Us2 = ( Us1 + M_PI if Us1 < M_PI,
c6541a0c D	124	( Us1 - M_PI otherwise;
0d969553 Y	125	Us2 corresponds to the maximum.
0d969553 Y	126	3- Calculate the extrema in [Uinf,Usup].
7fd59977	127	-----------------------------------------------------------------------------*/
	128	{
	129	myDone = Standard_False;
	130	myNbExt = 0;
	131
0d969553	132	// 1- Projection of the point P in the plane of circle -> Pp ...
7fd59977	133
	134	gp_Pnt O = C.Location();
	135	gp_Vec Axe (C.Axis().Direction());
	136	gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
	137	gp_Pnt Pp = P.Translated(Trsl);
	138
0d969553	139	// 2- Calculate u solutions in [0.,2.*PI] ...
7fd59977	140
	141	gp_Vec OPp (O,Pp);
	142	if (OPp.Magnitude() < Tol) { return; }
	143	Standard_Real Usol[2];
c6541a0c	144	Usol[0] = C.XAxis().Direction().AngleWithRef(OPp,Axe); // -M_PI<U1<M_PI
c764e804	145
	146	const Standard_Real aAngTol = Precision::Angular();
	147	if ( Usol[0] + M_PI < aAngTol )
	148	Usol[0] = -M_PI;
	149	else if ( Usol[0] - M_PI > -aAngTol )
	150	Usol[0] = M_PI;
	151
c6541a0c	152	Usol[1] = Usol[0] + M_PI;
7fd59977	153
7fd59977	154	Standard_Real myuinf = Uinf;
7fd59977	155	//Standard_Real TolU = Tol*C.Radius();
	156	Standard_Real TolU, aR;
	157	aR=C.Radius();
	158	TolU=Precision::Infinite();
	159	if (aR > gp::Resolution()) {
	160	TolU= Tol/aR;
	161	}
42ff8f5b	162	//
c6541a0c D	163	ElCLib::AdjustPeriodic(Uinf, Uinf+2*M_PI, TolU, myuinf, Usol[0]);
	164	ElCLib::AdjustPeriodic(Uinf, Uinf+2*M_PI, TolU, myuinf, Usol[1]);
	165	if (((Usol[0]-2M_PI-Uinf) < TolU) && ((Usol[0]-2M_PI-Uinf) > -TolU)) Usol[0] = Uinf;
	166	if (((Usol[1]-2M_PI-Uinf) < TolU) && ((Usol[1]-2M_PI-Uinf) > -TolU)) Usol[1] = Uinf;
7fd59977	167
7fd59977	168
0d969553	169	// 3- Calculate extrema in [Umin,Umax] ...
7fd59977	170
	171	gp_Pnt Cu;
	172	Standard_Real Us;
	173	for (Standard_Integer NoSol = 0; NoSol <= 1; NoSol++) {
	174	Us = Usol[NoSol];
	175	if (((Uinf-Us) < TolU) && ((Us-Usup) < TolU)) {
	176	Cu = ElCLib::Value(Us,C);
	177	mySqDist[myNbExt] = Cu.SquareDistance(P);
	178	myIsMin[myNbExt] = (NoSol == 0);
	179	myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
	180	myNbExt++;
	181	}
	182	}
	183	myDone = Standard_True;
	184	}
	185	//=============================================================================
	186
	187	Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
	188	const gp_Elips& C,
	189	const Standard_Real Tol,
	190	const Standard_Real Uinf,
	191	const Standard_Real Usup)
	192	{
	193	Perform(P, C, Tol, Uinf, Usup);
	194	}
	195
	196
	197
	198	void Extrema_ExtPElC::Perform (const gp_Pnt& P,
	199	const gp_Elips& C,
	200	const Standard_Real Tol,
	201	const Standard_Real Uinf,
	202	const Standard_Real Usup)
	203	/*-----------------------------------------------------------------------------
0d969553 Y	204	Function:
	205	Find values of parameter u so that:
	206	- dist(P,C(u)) passes by an extremum,
7fd59977	207	- Uinf <= u <= Usup.
7fd59977	208
0d969553 Y	209	Method:
	210	Takes 2 steps:
	211	1- Projection of point P in the plane of the ellipse,
	212	2- Calculation of the solutions:
	213	Let Pp, the projected point; find values u so that:
7fd59977	214	(C(u)-Pp).C'(u) = 0. (1)
0d969553 Y	215	Let Cos = cos(u) and Sin = sin(u),
	216	C(u) = (ACos,BSin) and Pp = (X,Y);
	217	Then, (1) <=> (ACos-X,BSin-Y).(-ASin,BCos) = 0.
7fd59977	218	(B2-A2)CosSin - BYCos + AXSin = 0.
0d969553	219	Use algorithm math_TrigonometricFunctionRoots to solve this equation.
7fd59977	220	-----------------------------------------------------------------------------*/
	221	{
	222	myDone = Standard_False;
	223	myNbExt = 0;
	224
0d969553	225	// 1- Projection of point P in the plane of the ellipse -> Pp ...
7fd59977	226
	227	gp_Pnt O = C.Location();
	228	gp_Vec Axe (C.Axis().Direction());
	229	gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
	230	gp_Pnt Pp = P.Translated(Trsl);
	231
0d969553	232	// 2- Calculation of solutions ...
7fd59977	233
	234	Standard_Integer NoSol, NbSol;
	235	Standard_Real A = C.MajorRadius();
	236	Standard_Real B = C.MinorRadius();
	237	gp_Vec OPp (O,Pp);
	238	Standard_Real OPpMagn = OPp.Magnitude();
	239	if (OPpMagn < Tol) { if (Abs(A-B) < Tol) { return; } }
	240	Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
	241	Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
	242	// Standard_Real Y = Sqrt(OPpMagnOPpMagn-XX);
	243
	244	Standard_Real ko2 = (BB-AA)/2., ko3 = -BY, ko4 = AX;
	245	if(Abs(ko3) < 1.e-16*Max(Abs(ko2), Abs(ko3))) ko3 = 0.0;
	246
	247	// math_TrigonometricFunctionRoots Sol(0.,(BB-AA)/2.,-BY,AX,0.,Uinf,Usup);
	248	math_TrigonometricFunctionRoots Sol(0.,ko2, ko3, ko4, 0.,Uinf,Usup);
	249
	250	if (!Sol.IsDone()) { return; }
	251	gp_Pnt Cu;
	252	Standard_Real Us;
	253	NbSol = Sol.NbSolutions();
	254	for (NoSol = 1; NoSol <= NbSol; NoSol++) {
	255	Us = Sol.Value(NoSol);
	256	Cu = ElCLib::Value(Us,C);
	257	mySqDist[myNbExt] = Cu.SquareDistance(P);
7fd59977	258	myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
a20400b3	259	Cu = ElCLib::Value(Us + 0.1, C);
a20400b3	260	myIsMin[myNbExt] = mySqDist[myNbExt] < Cu.SquareDistance(P);
7fd59977	261	myNbExt++;
	262	}
	263	myDone = Standard_True;
	264	}
	265	//=============================================================================
	266
	267	Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
	268	const gp_Hypr& C,
	269	const Standard_Real Tol,
	270	const Standard_Real Uinf,
	271	const Standard_Real Usup)
	272	{
	273	Perform(P, C, Tol, Uinf, Usup);
	274	}
	275
	276
	277	void Extrema_ExtPElC::Perform(const gp_Pnt& P,
	278	const gp_Hypr& C,
	279	const Standard_Real Tol,
	280	const Standard_Real Uinf,
	281	const Standard_Real Usup)
	282	/*-----------------------------------------------------------------------------
0d969553 Y	283	Function:
	284	Find values of parameter u so that:
	285	- dist(P,C(u)) passes by an extremum,
7fd59977	286	- Uinf <= u <= Usup.
7fd59977	287
0d969553 Y	288	Method:
	289	Takes 2 steps:
	290	1- Projection of point P in the plane of the hyperbola,
	291	2- Calculation of solutions:
	292	Let Pp, le point projete; on recherche les valeurs u telles que:
7fd59977	293	(C(u)-Pp).C'(u) = 0. (1)
0d969553 Y	294	Let R and r be the radiuses of the hyperbola,
	295	Chu = Cosh(u) and Shu = Sinh(u),
	296	C(u) = (RChu,rShu) and Pp = (X,Y);
	297	Then, (1) <=> (RChu-X,rShu-Y).(RShu,rChu) = 0.
7fd59977	298	(R2+r2)ChuShu - XRShu - YrChu = 0. (2)
0d969553 Y	299	Let v = e**u;
0d969553 Y	300	Then, by using Chu = (eu+e(-u))/2. and Sh = (eu-e(-u)))/2.
7fd59977	301	(2) <=> ((R2+r2)/4.) * (v2-v(-2)) -
	302	((XR+Yr)/2.) * v +
	303	((XR-Yr)/2.) * v**(-1) = 0.
	304	(2)* v2 <=> ((R2+r*2)/4.) v**4 -
	305	((XR+Yr)/2.) * v**3 +
	306	((XR-Yr)/2.) * v -
	307	((R2+r2)/4.) = 0.
0d969553	308	Use algorithm math_DirectPolynomialRoots to solve this equation by v.
7fd59977	309	-----------------------------------------------------------------------------*/
	310	{
	311	myDone = Standard_False;
	312	myNbExt = 0;
	313
0d969553	314	// 1- Projection of point P in the plane of hyperbola -> Pp ...
7fd59977	315
	316	gp_Pnt O = C.Location();
	317	gp_Vec Axe (C.Axis().Direction());
	318	gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
	319	gp_Pnt Pp = P.Translated(Trsl);
	320
0d969553	321	// 2- Calculation of solutions ...
7fd59977	322
	323	Standard_Real Tol2 = Tol * Tol;
	324	Standard_Real R = C.MajorRadius();
	325	Standard_Real r = C.MinorRadius();
	326	gp_Vec OPp (O,Pp);
7fd59977	327	Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
	328	Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
	329
	330	Standard_Real C1 = (RR+rr)/4.;
	331	math_DirectPolynomialRoots Sol(C1,-(XR+Yr)/2.,0.,(XR-Yr)/2.,-C1);
	332	if (!Sol.IsDone()) { return; }
	333	gp_Pnt Cu;
	334	Standard_Real Us, Vs;
	335	Standard_Integer NbSol = Sol.NbSolutions();
	336	Standard_Boolean DejaEnr;
	337	Standard_Integer NoExt;
	338	gp_Pnt TbExt[4];
	339	for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
	340	Vs = Sol.Value(NoSol);
	341	if (Vs > 0.) {
	342	Us = Log(Vs);
	343	if ((Us >= Uinf) && (Us <= Usup)) {
	344	Cu = ElCLib::Value(Us,C);
	345	DejaEnr = Standard_False;
	346	for (NoExt = 0; NoExt < myNbExt; NoExt++) {
	347	if (TbExt[NoExt].SquareDistance(Cu) < Tol2) {
	348	DejaEnr = Standard_True;
	349	break;
	350	}
	351	}
	352	if (!DejaEnr) {
	353	TbExt[myNbExt] = Cu;
	354	mySqDist[myNbExt] = Cu.SquareDistance(P);
	355	myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
	356	myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
	357	myNbExt++;
	358	}
	359	} // if ((Us >= Uinf) && (Us <= Usup))
	360	} // if (Vs > 0.)
	361	} // for (Standard_Integer NoSol = 1; ...
	362	myDone = Standard_True;
	363	}
	364	//=============================================================================
	365
	366	Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
	367	const gp_Parab& C,
	368	const Standard_Real Tol,
	369	const Standard_Real Uinf,
	370	const Standard_Real Usup)
	371	{
	372	Perform(P, C, Tol, Uinf, Usup);
	373	}
	374
	375
	376	void Extrema_ExtPElC::Perform(const gp_Pnt& P,
	377	const gp_Parab& C,
	378	// const Standard_Real Tol,
	379	const Standard_Real ,
	380	const Standard_Real Uinf,
	381	const Standard_Real Usup)
	382	/*-----------------------------------------------------------------------------
0d969553 Y	383	Function:
	384	Find values of parameter u so that:
	385	- dist(P,C(u)) pass by an extremum,
7fd59977	386	- Uinf <= u <= Usup.
7fd59977	387
0d969553 Y	388	Method:
	389	Takes 2 steps:
	390	1- Projection of point P in the plane of the parabola,
	391	2- Calculation of solutions:
	392	Let Pp, the projected point; find values u so that:
7fd59977	393	(C(u)-Pp).C'(u) = 0. (1)
0d969553 Y	394	Let F the focus of the parabola,
0d969553 Y	395	C(u) = ((uu)/(4.F),u) and Pp = (X,Y);
7fd59977	396	Alors, (1) <=> ((uu)/(4.F)-X,u-Y).(u/(2.*F),1) = 0.
7fd59977	397	(1./(4.F)) U*3 + (2.F-X) * U - 2FY = 0.
0d969553	398	Use algorithm math_DirectPolynomialRoots to solve this equation by U.
7fd59977	399	-----------------------------------------------------------------------------*/
	400	{
	401	myDone = Standard_False;
	402	myNbExt = 0;
	403
0d969553	404	// 1- Projection of point P in the plane of the parabola -> Pp ...
7fd59977	405
	406	gp_Pnt O = C.Location();
	407	gp_Vec Axe (C.Axis().Direction());
	408	gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
	409	gp_Pnt Pp = P.Translated(Trsl);
	410
0d969553	411	// 2- Calculation of solutions ...
7fd59977	412
	413	Standard_Real F = C.Focal();
	414	gp_Vec OPp (O,Pp);
7fd59977	415	Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
	416	// Standard_Real Y = Sqrt(OPpMagnOPpMagn-XX);
	417	Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
	418	math_DirectPolynomialRoots Sol(1./(4.F),0.,2.F-X,-2.FY);
	419	if (!Sol.IsDone()) { return; }
	420	gp_Pnt Cu;
	421	Standard_Real Us;
	422	Standard_Integer NbSol = Sol.NbSolutions();
	423	Standard_Boolean DejaEnr;
	424	Standard_Integer NoExt;
	425	gp_Pnt TbExt[3];
	426	for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
	427	Us = Sol.Value(NoSol);
	428	if ((Us >= Uinf) && (Us <= Usup)) {
	429	Cu = ElCLib::Value(Us,C);
	430	DejaEnr = Standard_False;
	431	for (NoExt = 0; NoExt < myNbExt; NoExt++) {
638ad7f3	432	if (TbExt[NoExt].SquareDistance(Cu) < Precision::SquareConfusion()) {
7fd59977	433	DejaEnr = Standard_True;
	434	break;
	435	}
	436	}
	437	if (!DejaEnr) {
	438	TbExt[myNbExt] = Cu;
	439	mySqDist[myNbExt] = Cu.SquareDistance(P);
	440	myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
	441	myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
	442	myNbExt++;
	443	}
	444	} // if ((Us >= Uinf) && (Us <= Usup))
	445	} // for (Standard_Integer NoSol = 1; ...
	446	myDone = Standard_True;
	447	}
	448	//=============================================================================
	449
	450	Standard_Boolean Extrema_ExtPElC::IsDone () const { return myDone; }
	451	//=============================================================================
	452
	453	Standard_Integer Extrema_ExtPElC::NbExt () const
	454	{
9775fa61	455	if (!IsDone()) { throw StdFail_NotDone(); }
7fd59977	456	return myNbExt;
	457	}
	458	//=============================================================================
	459
	460	Standard_Real Extrema_ExtPElC::SquareDistance (const Standard_Integer N) const
	461	{
9775fa61	462	if ((N < 1) \|\| (N > NbExt())) { throw Standard_OutOfRange(); }
7fd59977	463	return mySqDist[N-1];
	464	}
	465	//=============================================================================
	466
	467	Standard_Boolean Extrema_ExtPElC::IsMin (const Standard_Integer N) const
	468	{
9775fa61	469	if ((N < 1) \|\| (N > NbExt())) { throw Standard_OutOfRange(); }
7fd59977	470	return myIsMin[N-1];
	471	}
	472	//=============================================================================
	473
5d99f2c8	474	const Extrema_POnCurv& Extrema_ExtPElC::Point (const Standard_Integer N) const
7fd59977	475	{
9775fa61	476	if ((N < 1) \|\| (N > NbExt())) { throw Standard_OutOfRange(); }
7fd59977	477	return myPoint[N-1];
	478	}
	479	//=============================================================================