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

// Copyright (c) 1995-1999 Matra Datavision
// Copyright (c) 1999-2012 OPEN CASCADE SAS
//
// The content of this file is subject to the Open CASCADE Technology Public
// License Version 6.5 (the "License"). You may not use the content of this file
// except in compliance with the License. Please obtain a copy of the License
// at http://www.opencascade.org and read it completely before using this file.
//
// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
//
// The Original Code and all software distributed under the License is
// distributed on an "AS IS" basis, without warranty of any kind, and the
// Initial Developer hereby disclaims all such warranties, including without
// limitation, any warranties of merchantability, fitness for a particular
// purpose or non-infringement. Please see the License for the specific terms
// and conditions governing the rights and limitations under the License.

#include <Extrema_ExtPElC.ixx>
#include <StdFail_NotDone.hxx>
#include <math_DirectPolynomialRoots.hxx>
#include <math_TrigonometricFunctionRoots.hxx>
#include <ElCLib.hxx>
#include <Standard_OutOfRange.hxx>
#include <Standard_NotImplemented.hxx>
#include <Precision.hxx>

//=======================================================================
//function : Extrema_ExtPElC
//purpose  : 
//=======================================================================
Extrema_ExtPElC::Extrema_ExtPElC () { myDone = 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
  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()) { StdFail_NotDone::Raise(); }
  return myNbExt;
}
//=============================================================================

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

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

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