[occt.git] / src / IntAna / IntAna_Curve.cxx

// Created on: 1992-06-30
// Created by: Laurent BUCHARD
// Copyright (c) 1992-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.

#ifndef OCCT_DEBUG
#define No_Standard_RangeError
#define No_Standard_OutOfRange
#endif

//----------------------------------------------------------------------
//-- Differents constructeurs sont proposes qui correspondent aux
//-- polynomes en Z :  
//--    A(Sin(Theta),Cos(Theta)) Z**2 
//--  + B(Sin(Theta),Cos(Theta)) Z
//--  + C(Sin(Theta),Cos(Theta))
//--
//-- Une Courbe est definie sur un domaine 
//--
//-- Value retourne le point de parametre U(Theta),V(Theta)
//--       ou V est la solution du polynome A V**2 + B V + C
//--       (Selon les cas, on prend V+ ou V-)
//--
//-- D1u   calcule le vecteur tangent a la courbe 
//--       et retourne le booleen Standard_False si ce calcul ne peut
//--       pas etre mene a bien.  
//----------------------------------------------------------------------

#include <algorithm>

#include <ElSLib.hxx>
#include <gp_Cone.hxx>
#include <gp_Cylinder.hxx>
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>
#include <gp_XYZ.hxx>
#include <IntAna_Curve.hxx>
#include <math_DirectPolynomialRoots.hxx>
#include <Precision.hxx>
#include <Standard_DomainError.hxx>

//=======================================================================
//function : IntAna_Curve
//purpose  : 
//=======================================================================
IntAna_Curve::IntAna_Curve()
{
  typequadric=GeomAbs_OtherSurface;
  firstbounded=Standard_False;
  lastbounded=Standard_False;
}
//=======================================================================
//function : SetConeQuadValues
//purpose  : Description de l intersection Cone Quadrique
//=======================================================================
void IntAna_Curve::SetConeQuadValues(const gp_Cone& Cone,
                                     const Standard_Real Qxx,
                                     const Standard_Real Qyy,
                                     const Standard_Real Qzz,
                                     const Standard_Real Qxy,
                                     const Standard_Real Qxz,
                                     const Standard_Real Qyz,
                                     const Standard_Real Qx,
                                     const Standard_Real Qy,
                                     const Standard_Real Qz,
                                     const Standard_Real Q1,
                                     const Standard_Real TOL,
                                     const Standard_Real DomInf,
                                     const Standard_Real DomSup,
                                     const Standard_Boolean twocurves,
                                     const Standard_Boolean takezpositive)
{
  
  Ax3        = Cone.Position();
  RCyl       = Cone.RefRadius();

  Angle      = Cone.SemiAngle();
  Standard_Real UnSurTgAngle = 1.0/(Tan(Cone.SemiAngle()));

  typequadric= GeomAbs_Cone;
  
  TwoCurves     = twocurves;         //-- deux  Z pour un meme parametre 
  TakeZPositive = takezpositive;     //-- Prendre sur la courbe le Z Positif
                                     //--   ( -B + Sqrt()) et non (-B - Sqrt())
  
  
  Z0Cte      = Q1;                   //-- Attention On a    Z?Cos Cos(t)
  Z0Sin      = 0.0;                  //-- et Non          2 Z?Cos Cos(t) !!!
  Z0Cos      = 0.0;                  //-- Ce pour tous les Parametres
  Z0CosCos   = 0.0;                  //--  ie pas de Coefficient 2  
  Z0SinSin   = 0.0;                  //--     devant les termes CS C S 
  Z0CosSin   = 0.0;
  
  Z1Cte      = 2.0*(UnSurTgAngle)*Qz;
  Z1Sin      = Qy+Qy;
  Z1Cos      = Qx+Qx;
  Z1CosCos   = 0.0;
  Z1SinSin   = 0.0;
  Z1CosSin   = 0.0;
  
  Z2Cte      = Qzz * UnSurTgAngle*UnSurTgAngle;
  Z2Sin      = (UnSurTgAngle+UnSurTgAngle)*Qyz;
  Z2Cos      = (UnSurTgAngle+UnSurTgAngle)*Qxz;
  Z2CosCos   = Qxx;
  Z2SinSin   = Qyy;
  Z2CosSin   = Qxy;

  Tolerance  = TOL;
  DomainInf = DomInf;
  DomainSup = DomSup;
  
  RestrictedInf = RestrictedSup = Standard_True;   //-- Le Domaine est Borne
  firstbounded = lastbounded = Standard_False;

  myFirstParameter = DomainInf;
  myLastParameter = (TwoCurves) ? DomainSup + DomainSup - DomainInf :
                                  DomainSup;
}

//=======================================================================
//function : SetCylinderQuadValues
//purpose  : Description de l intersection Cylindre Quadrique
//=======================================================================
void IntAna_Curve::SetCylinderQuadValues(const gp_Cylinder& Cyl,
                                         const Standard_Real Qxx,
                                         const Standard_Real Qyy,
                                         const Standard_Real Qzz,
                                         const Standard_Real Qxy,
                                         const Standard_Real Qxz,
                                         const Standard_Real Qyz,
                                         const Standard_Real Qx,
                                         const Standard_Real Qy,
                                         const Standard_Real Qz,
                                         const Standard_Real Q1,
                                         const Standard_Real TOL,
                                         const Standard_Real DomInf,
                                         const Standard_Real DomSup,
                                         const Standard_Boolean twocurves,
                                         const Standard_Boolean takezpositive)
{
  
  Ax3        = Cyl.Position();
  RCyl       = Cyl.Radius();
  typequadric= GeomAbs_Cylinder;
  
  TwoCurves     = twocurves;         //-- deux  Z pour un meme parametre 
  TakeZPositive = takezpositive;     //-- Prendre sur la courbe le Z Positif
  Standard_Real RCylmul2 = RCyl+RCyl;         //--   ( -B + Sqrt()) 

  Z0Cte      = Q1;
  Z0Sin      = RCylmul2*Qy;
  Z0Cos      = RCylmul2*Qx;
  Z0CosCos   = Qxx*RCyl*RCyl;
  Z0SinSin   = Qyy*RCyl*RCyl;
  Z0CosSin   = RCyl*RCyl*Qxy;

  Z1Cte      = Qz+Qz;
  Z1Sin      = RCylmul2*Qyz;
  Z1Cos      = RCylmul2*Qxz;
  Z1CosCos   = 0.0;
  Z1SinSin   = 0.0;
  Z1CosSin   = 0.0;

  Z2Cte      = Qzz;
  Z2Sin      = 0.0;
  Z2Cos      = 0.0;
  Z2CosCos   = 0.0;
  Z2SinSin   = 0.0;
  Z2CosSin   = 0.0;

  Tolerance  = TOL;
  DomainInf = DomInf;
  DomainSup = DomSup;
  
  RestrictedInf = RestrictedSup = Standard_True;
  firstbounded = lastbounded = Standard_False;

  myFirstParameter = DomainInf;
  myLastParameter  = (TwoCurves) ? DomainSup + DomainSup - DomainInf :
                                   DomainSup;
}

//=======================================================================
//function : IsOpen
//purpose  : 
//=======================================================================
Standard_Boolean IntAna_Curve::IsOpen() const
{
  return(RestrictedInf && RestrictedSup);
}

//=======================================================================
//function : Domain
//purpose  : 
//=======================================================================
void IntAna_Curve::Domain(Standard_Real& theFirst,
                          Standard_Real& theLast) const
{
  if (RestrictedInf && RestrictedSup)
  {
    theFirst = myFirstParameter;
    theLast = myLastParameter;
  }
  else
  {
    throw Standard_DomainError("IntAna_Curve::Domain");
  }
}
//=======================================================================
//function : IsConstant
//purpose  : 
//=======================================================================
Standard_Boolean IntAna_Curve::IsConstant() const
{
  //-- ???  Pas facile de decider a la seule vue des Param.
  return(Standard_False);
}

//=======================================================================
//function : IsFirstOpen
//purpose  : 
//=======================================================================
Standard_Boolean IntAna_Curve::IsFirstOpen() const
{
  return(firstbounded);
}

//=======================================================================
//function : IsLastOpen
//purpose  : 
//=======================================================================
Standard_Boolean IntAna_Curve::IsLastOpen() const
{
  return(lastbounded);
}
//=======================================================================
//function : SetIsFirstOpen
//purpose  : 
//=======================================================================
void IntAna_Curve::SetIsFirstOpen(const Standard_Boolean Flag)
{
  firstbounded = Flag;
}

//=======================================================================
//function : SetIsLastOpen
//purpose  : 
//=======================================================================
void IntAna_Curve::SetIsLastOpen(const Standard_Boolean Flag)
{
  lastbounded = Flag;
}

//=======================================================================
//function : InternalUVValue
//purpose  : 
//=======================================================================
void IntAna_Curve::InternalUVValue(const Standard_Real theta,
                                   Standard_Real& Param1,
                                   Standard_Real& Param2,
                                   Standard_Real& A,
                                   Standard_Real& B,
                                   Standard_Real& C,
                                   Standard_Real& cost,
                                   Standard_Real& sint,
                                   Standard_Real& SigneSqrtDis) const
{
  const Standard_Real aRelTolp = 1.0+Epsilon(1.0), aRelTolm = 1.0-Epsilon(1.0);
  
  // Infinitesimal step of increasing curve parameter. See comment below.
  const Standard_Real aDT = 100.0*Epsilon(DomainSup + DomainSup - DomainInf);

  Standard_Real Theta=theta;
  Standard_Boolean SecondSolution=Standard_False; 

  if ((Theta<DomainInf*aRelTolm) ||
      ((Theta>DomainSup*aRelTolp) && (!TwoCurves)) ||
      (Theta>(DomainSup + DomainSup - DomainInf)*aRelTolp))
  {
    SigneSqrtDis = 0.;
    throw Standard_DomainError("IntAna_Curve::Domain");
  }
  
  if (Abs(Theta - DomainSup) < aDT)
  {
    // Point of Null-discriminant.
    Theta = DomainSup;
  }
  else if (Theta>DomainSup)
  {
    Theta = DomainSup + DomainSup - Theta;
    SecondSolution=Standard_True; 
  }

  Param1=Theta;

  if(!TwoCurves) { 
    SecondSolution=TakeZPositive; 
  }
  //
  cost = Cos(Theta);
  sint = Sin(Theta);
  const Standard_Real aSin2t = Sin(Theta + Theta);
  const Standard_Real aCos2t = Cos(Theta + Theta);
  
  A=Z2Cte+sint*(Z2Sin+sint*Z2SinSin)+cost*(Z2Cos+cost*Z2CosCos)
    + Z2CosSin*aSin2t;
  
  const Standard_Real aDA = cost*Z2Sin - sint*Z2Cos + 
                            aSin2t*(Z2SinSin - Z2CosCos) + 
                            aCos2t*(Z2CosSin * Z2CosSin);

  B=Z1Cte+sint*(Z1Sin+sint*Z1SinSin)+cost*(Z1Cos+cost*Z1CosCos)
    + Z1CosSin*aSin2t;

  const Standard_Real aDB = Z1Sin*cost - Z1Cos*sint +
                            aSin2t*(Z1SinSin - Z1CosCos) +
                            aCos2t*(Z1CosSin + Z1CosSin);
  
  C=Z0Cte+sint*(Z0Sin+sint*Z0SinSin)+cost*(Z0Cos+cost*Z0CosCos)
    + Z0CosSin*aSin2t;
  
  const Standard_Real aDC = Z0Sin*cost - Z0Cos*sint +
                            aSin2t*(Z0SinSin - Z0CosCos) + 
                            aCos2t*(Z0CosSin + Z0CosSin);

  Standard_Real aDiscriminant = B*B-4.0*A*C;

  // We consider that infinitesimal dt = aDT.
  // Error of discriminant computation is equal to
  // (d(Disc)/dt)*dt, where 1st derivative d(Disc)/dt = 2*B*aDB - 4*(A*aDC + C*aDA).

  const Standard_Real aTolD = 2.0*aDT*Abs(B*aDB - 2.0*(A*aDC + C*aDA));
  
  if (aDiscriminant < aTolD)
    aDiscriminant = 0.0;

  if (Abs(A) <= Precision::PConfusion())
  {
    if (Abs(B) <= Precision::PConfusion())
    {
      Param2 = 0.0;
    }
    else
    {
      Param2 = -C / B;
    }
  }
  else
  {
    SigneSqrtDis = (SecondSolution) ? Sqrt(aDiscriminant) : -Sqrt(aDiscriminant);
    Param2 = (-B + SigneSqrtDis) / (A + A);
  }
}

//=======================================================================
//function : Value
//purpose  : 
//=======================================================================
gp_Pnt IntAna_Curve::Value(const Standard_Real theta)
{
  Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
  //
  A=0.0;  B=0.0;   C=0.0;
  U=0.0;  V=0.0;  
  sint=0.0;  cost=0.0;
  SigneSqrtDis=0.0;
  InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis); 
  //-- checked the parameter U and Raises Domain Error if Error
  return(InternalValue(U,V));
}
//=======================================================================
//function : D1u
//purpose  : 
//=======================================================================
Standard_Boolean IntAna_Curve::D1u(const Standard_Real theta,
                                   gp_Pnt& Pt,
                                   gp_Vec& Vec)
{
  //-- Pour detecter le cas ou le calcul est impossible 
  Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
  A=0.0;  B=0.0;   C=0.0;
  U=0.0;  V=0.0;  
  sint=0.0;  cost=0.0;
  //
  InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
  //
  Pt = Value(theta);
  if (Abs(A)<1.0e-7 || Abs(SigneSqrtDis)<1.0e-10) return(Standard_False);
  
  
  //-- Approximation de la derivee (mieux que le calcul mathematique!)
  Standard_Real dtheta = (DomainSup - DomainInf)*1.0e-6;
  Standard_Real theta2 = theta+dtheta;
  if ((theta2<DomainInf) || ((theta2>DomainSup) && (!TwoCurves))
      || (theta2>(DomainSup + DomainSup - DomainInf + 1.0e-14)))
  {
    dtheta = -dtheta;
    theta2 = theta+dtheta;
  }
  gp_Pnt P2 = Value(theta2);
  dtheta = 1.0/dtheta;
  Vec.SetCoord((P2.X()-Pt.X())*dtheta,
	       (P2.Y()-Pt.Y())*dtheta,
	       (P2.Z()-Pt.Z())*dtheta);
  
  return(Standard_True);
}  
//=======================================================================
//function : FindParameter
//purpose  : Projects P to the ALine. Returns the list of parameters as a results
//            of projection.
//           Sometimes aline can be self-intersected line (see bug #29807 where 
//            ALine goes through the cone apex).
//=======================================================================
void IntAna_Curve::FindParameter(const gp_Pnt& theP,
                                 TColStd_ListOfReal& theParams) const
{
  const Standard_Real aPIpPI = M_PI + M_PI,
                      anEpsAng = 1.e-8,
                      aSqTolPrecision=1.0e-8;
  Standard_Real aTheta = 0.0;
  //
  switch (typequadric)
  {
    case GeomAbs_Cylinder:
    {
      Standard_Real aZ;
      ElSLib::CylinderParameters(Ax3, RCyl, theP, aTheta, aZ);
    }
    break;

    case GeomAbs_Cone:
    {
      Standard_Real aZ;
      ElSLib::ConeParameters(Ax3, RCyl, Angle, theP, aTheta, aZ);
    }
    break;

    default:
      return;
  }
  //
  if (!firstbounded && (DomainInf > aTheta) && ((DomainInf - aTheta) <= anEpsAng))
  {
    aTheta = DomainInf;
  }
  else if (!lastbounded && (aTheta > DomainSup) && ((aTheta - DomainSup) <= anEpsAng))
  {
    aTheta = DomainSup;
  }
  //
  if (aTheta < DomainInf)
  {
    aTheta = aTheta + aPIpPI;
  }
  else if (aTheta > DomainSup)
  {
    aTheta = aTheta - aPIpPI;
  }

  const Standard_Integer aMaxPar = 5;
  Standard_Real aParams[aMaxPar] = {DomainInf, DomainSup, aTheta,
                                    (TwoCurves)? DomainSup + DomainSup - aTheta : RealLast(),
                                    (TwoCurves) ? DomainSup + DomainSup - DomainInf : RealLast()};

  std::sort(aParams, aParams + aMaxPar - 1);

  for (Standard_Integer i = 0; i < aMaxPar; i++)
  {
    if (aParams[i] > myLastParameter)
      break;
    
    if (aParams[i] < myFirstParameter)
      continue;

    if (i && (aParams[i] - aParams[i - 1]) < Precision::PConfusion())
      continue;

    Standard_Real U = 0.0, V= 0.0, 
                  A = 0.0, B = 0.0, C = 0.0,
                  sint = 0.0, cost = 0.0, SigneSqrtDis = 0.0;
    InternalUVValue(aParams[i], U, V, A, B, C, 
                    cost, sint, SigneSqrtDis);
    const gp_Pnt aP(InternalValue(U, V));
    if (aP.SquareDistance(theP) < aSqTolPrecision)
    {
      theParams.Append(aParams[i]);
    }
  }
}
//=======================================================================
//function : InternalValue
//purpose  : 
//=======================================================================
gp_Pnt IntAna_Curve::InternalValue(const Standard_Real U,
                                   const Standard_Real _V) const
{
  //-- cout<<" ["<<U<<","<<V<<"]";
  Standard_Real V = _V;
  if(V > 100000.0 )   {   V= 100000.0; }       
  if(V < -100000.0 )  {   V=-100000.0; }      

  switch(typequadric) {
  case GeomAbs_Cone:
    {
      //------------------------------------------------
      //-- Parametrage : X = V * Cos(U)              ---
      //--               Y = V * Sin(U)              ---
      //--               Z = (V-RCyl) / Tan(SemiAngle)--
      //------------------------------------------------ 
      //-- Angle Vaut Cone.SemiAngle()  
      return(ElSLib::ConeValue(U,(V-RCyl)/Sin(Angle),Ax3,RCyl,Angle));
    }
    break;
    
  case GeomAbs_Cylinder:
    return(ElSLib::CylinderValue(U,V,Ax3,RCyl)); 
  case GeomAbs_Sphere:
    return(ElSLib::SphereValue(U,V,Ax3,RCyl)); 
  default:
    return(gp_Pnt(0.0,0.0,0.0));
  }
}

//
//=======================================================================
//function : SetDomain
//purpose  : 
//=======================================================================
void IntAna_Curve::SetDomain(const Standard_Real theFirst,
                             const Standard_Real theLast)
{
  if (theLast <= theFirst)
  {
    throw Standard_DomainError("IntAna_Curve::Domain");
  }
  //
  myFirstParameter = theFirst;
  myLastParameter = theLast;
}
Commit	Line	Data
b311480e	1	// Created on: 1992-06-30
	2	// Created by: Laurent BUCHARD
	3	// Copyright (c) 1992-1999 Matra Datavision
973c2be1	4	// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e	5	//
973c2be1	6	// This file is part of Open CASCADE Technology software library.
b311480e	7	//
d5f74e42	8	// This library is free software; you can redistribute it and/or modify it under
d5f74e42	9	// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1	10	// by the Free Software Foundation, with special exception defined in the file
	11	// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
	12	// distribution for complete text of the license and disclaimer of any warranty.
b311480e	13	//
973c2be1	14	// Alternatively, this file may be used under the terms of Open CASCADE
973c2be1	15	// commercial license or contractual agreement.
7fd59977	16
0797d9d3	17	#ifndef OCCT_DEBUG
7fd59977	18	#define No_Standard_RangeError
	19	#define No_Standard_OutOfRange
	20	#endif
	21
	22	//----------------------------------------------------------------------
	23	//-- Differents constructeurs sont proposes qui correspondent aux
	24	//-- polynomes en Z :
	25	//-- A(Sin(Theta),Cos(Theta)) Z**2
	26	//-- + B(Sin(Theta),Cos(Theta)) Z
	27	//-- + C(Sin(Theta),Cos(Theta))
	28	//--
	29	//-- Une Courbe est definie sur un domaine
	30	//--
	31	//-- Value retourne le point de parametre U(Theta),V(Theta)
	32	//-- ou V est la solution du polynome A V**2 + B V + C
	33	//-- (Selon les cas, on prend V+ ou V-)
	34	//--
	35	//-- D1u calcule le vecteur tangent a la courbe
	36	//-- et retourne le booleen Standard_False si ce calcul ne peut
	37	//-- pas etre mene a bien.
	38	//----------------------------------------------------------------------
	39
3306fdd9	40	#include <algorithm>
3306fdd9	41
7fd59977	42	#include <ElSLib.hxx>
42cf5bc1	43	#include <gp_Cone.hxx>
	44	#include <gp_Cylinder.hxx>
	45	#include <gp_Pnt.hxx>
	46	#include <gp_Vec.hxx>
7fd59977	47	#include <gp_XYZ.hxx>
42cf5bc1	48	#include <IntAna_Curve.hxx>
	49	#include <math_DirectPolynomialRoots.hxx>
	50	#include <Precision.hxx>
	51	#include <Standard_DomainError.hxx>
7fd59977	52
	53	//=======================================================================
	54	//function : IntAna_Curve
	55	//purpose :
	56	//=======================================================================
3306fdd9	57	IntAna_Curve::IntAna_Curve()
7fd59977	58	{
	59	typequadric=GeomAbs_OtherSurface;
	60	firstbounded=Standard_False;
	61	lastbounded=Standard_False;
	62	}
	63	//=======================================================================
	64	//function : SetConeQuadValues
	65	//purpose : Description de l intersection Cone Quadrique
	66	//=======================================================================
3306fdd9	67	void IntAna_Curve::SetConeQuadValues(const gp_Cone& Cone,
	68	const Standard_Real Qxx,
	69	const Standard_Real Qyy,
	70	const Standard_Real Qzz,
	71	const Standard_Real Qxy,
	72	const Standard_Real Qxz,
	73	const Standard_Real Qyz,
	74	const Standard_Real Qx,
	75	const Standard_Real Qy,
	76	const Standard_Real Qz,
	77	const Standard_Real Q1,
	78	const Standard_Real TOL,
	79	const Standard_Real DomInf,
	80	const Standard_Real DomSup,
	81	const Standard_Boolean twocurves,
	82	const Standard_Boolean takezpositive)
7fd59977	83	{
	84
	85	Ax3 = Cone.Position();
	86	RCyl = Cone.RefRadius();
	87
	88	Angle = Cone.SemiAngle();
	89	Standard_Real UnSurTgAngle = 1.0/(Tan(Cone.SemiAngle()));
	90
	91	typequadric= GeomAbs_Cone;
	92
	93	TwoCurves = twocurves; //-- deux Z pour un meme parametre
	94	TakeZPositive = takezpositive; //-- Prendre sur la courbe le Z Positif
	95	//-- ( -B + Sqrt()) et non (-B - Sqrt())
	96
	97
	98	Z0Cte = Q1; //-- Attention On a Z?Cos Cos(t)
	99	Z0Sin = 0.0; //-- et Non 2 Z?Cos Cos(t) !!!
	100	Z0Cos = 0.0; //-- Ce pour tous les Parametres
	101	Z0CosCos = 0.0; //-- ie pas de Coefficient 2
	102	Z0SinSin = 0.0; //-- devant les termes CS C S
	103	Z0CosSin = 0.0;
	104
	105	Z1Cte = 2.0(UnSurTgAngle)Qz;
	106	Z1Sin = Qy+Qy;
	107	Z1Cos = Qx+Qx;
	108	Z1CosCos = 0.0;
	109	Z1SinSin = 0.0;
	110	Z1CosSin = 0.0;
	111
	112	Z2Cte = Qzz * UnSurTgAngle*UnSurTgAngle;
	113	Z2Sin = (UnSurTgAngle+UnSurTgAngle)*Qyz;
	114	Z2Cos = (UnSurTgAngle+UnSurTgAngle)*Qxz;
	115	Z2CosCos = Qxx;
	116	Z2SinSin = Qyy;
3306fdd9	117	Z2CosSin = Qxy;
7fd59977	118
7fd59977	119	Tolerance = TOL;
3306fdd9	120	DomainInf = DomInf;
3306fdd9	121	DomainSup = DomSup;
7fd59977	122
	123	RestrictedInf = RestrictedSup = Standard_True; //-- Le Domaine est Borne
	124	firstbounded = lastbounded = Standard_False;
3306fdd9	125
	126	myFirstParameter = DomainInf;
	127	myLastParameter = (TwoCurves) ? DomainSup + DomainSup - DomainInf :
	128	DomainSup;
7fd59977	129	}
	130
	131	//=======================================================================
	132	//function : SetCylinderQuadValues
	133	//purpose : Description de l intersection Cylindre Quadrique
	134	//=======================================================================
3306fdd9	135	void IntAna_Curve::SetCylinderQuadValues(const gp_Cylinder& Cyl,
	136	const Standard_Real Qxx,
	137	const Standard_Real Qyy,
	138	const Standard_Real Qzz,
	139	const Standard_Real Qxy,
	140	const Standard_Real Qxz,
	141	const Standard_Real Qyz,
	142	const Standard_Real Qx,
	143	const Standard_Real Qy,
	144	const Standard_Real Qz,
	145	const Standard_Real Q1,
	146	const Standard_Real TOL,
	147	const Standard_Real DomInf,
	148	const Standard_Real DomSup,
	149	const Standard_Boolean twocurves,
	150	const Standard_Boolean takezpositive)
7fd59977	151	{
	152
	153	Ax3 = Cyl.Position();
	154	RCyl = Cyl.Radius();
	155	typequadric= GeomAbs_Cylinder;
	156
	157	TwoCurves = twocurves; //-- deux Z pour un meme parametre
	158	TakeZPositive = takezpositive; //-- Prendre sur la courbe le Z Positif
	159	Standard_Real RCylmul2 = RCyl+RCyl; //-- ( -B + Sqrt())
	160
	161	Z0Cte = Q1;
	162	Z0Sin = RCylmul2*Qy;
	163	Z0Cos = RCylmul2*Qx;
	164	Z0CosCos = QxxRCylRCyl;
	165	Z0SinSin = QyyRCylRCyl;
3306fdd9	166	Z0CosSin = RCylRCylQxy;
7fd59977	167
	168	Z1Cte = Qz+Qz;
	169	Z1Sin = RCylmul2*Qyz;
	170	Z1Cos = RCylmul2*Qxz;
	171	Z1CosCos = 0.0;
	172	Z1SinSin = 0.0;
	173	Z1CosSin = 0.0;
	174
	175	Z2Cte = Qzz;
	176	Z2Sin = 0.0;
	177	Z2Cos = 0.0;
	178	Z2CosCos = 0.0;
	179	Z2SinSin = 0.0;
	180	Z2CosSin = 0.0;
	181
	182	Tolerance = TOL;
3306fdd9	183	DomainInf = DomInf;
3306fdd9	184	DomainSup = DomSup;
7fd59977	185
	186	RestrictedInf = RestrictedSup = Standard_True;
	187	firstbounded = lastbounded = Standard_False;
3306fdd9	188
	189	myFirstParameter = DomainInf;
	190	myLastParameter = (TwoCurves) ? DomainSup + DomainSup - DomainInf :
	191	DomainSup;
7fd59977	192	}
	193
	194	//=======================================================================
	195	//function : IsOpen
	196	//purpose :
	197	//=======================================================================
3306fdd9	198	Standard_Boolean IntAna_Curve::IsOpen() const
7fd59977	199	{
	200	return(RestrictedInf && RestrictedSup);
	201	}
	202
	203	//=======================================================================
	204	//function : Domain
	205	//purpose :
	206	//=======================================================================
3306fdd9	207	void IntAna_Curve::Domain(Standard_Real& theFirst,
3306fdd9	208	Standard_Real& theLast) const
7fd59977	209	{
3306fdd9	210	if (RestrictedInf && RestrictedSup)
	211	{
	212	theFirst = myFirstParameter;
	213	theLast = myLastParameter;
7fd59977	214	}
3306fdd9	215	else
3306fdd9	216	{
9775fa61	217	throw Standard_DomainError("IntAna_Curve::Domain");
7fd59977	218	}
	219	}
	220	//=======================================================================
	221	//function : IsConstant
	222	//purpose :
	223	//=======================================================================
3306fdd9	224	Standard_Boolean IntAna_Curve::IsConstant() const
7fd59977	225	{
	226	//-- ??? Pas facile de decider a la seule vue des Param.
	227	return(Standard_False);
	228	}
	229
	230	//=======================================================================
	231	//function : IsFirstOpen
	232	//purpose :
	233	//=======================================================================
3306fdd9	234	Standard_Boolean IntAna_Curve::IsFirstOpen() const
7fd59977	235	{
	236	return(firstbounded);
	237	}
	238
	239	//=======================================================================
	240	//function : IsLastOpen
	241	//purpose :
	242	//=======================================================================
3306fdd9	243	Standard_Boolean IntAna_Curve::IsLastOpen() const
7fd59977	244	{
	245	return(lastbounded);
	246	}
	247	//=======================================================================
	248	//function : SetIsFirstOpen
	249	//purpose :
	250	//=======================================================================
3306fdd9	251	void IntAna_Curve::SetIsFirstOpen(const Standard_Boolean Flag)
7fd59977	252	{
	253	firstbounded = Flag;
	254	}
	255
	256	//=======================================================================
	257	//function : SetIsLastOpen
	258	//purpose :
	259	//=======================================================================
3306fdd9	260	void IntAna_Curve::SetIsLastOpen(const Standard_Boolean Flag)
7fd59977	261	{
	262	lastbounded = Flag;
	263	}
	264
	265	//=======================================================================
	266	//function : InternalUVValue
	267	//purpose :
	268	//=======================================================================
3306fdd9	269	void IntAna_Curve::InternalUVValue(const Standard_Real theta,
	270	Standard_Real& Param1,
	271	Standard_Real& Param2,
	272	Standard_Real& A,
	273	Standard_Real& B,
	274	Standard_Real& C,
	275	Standard_Real& cost,
	276	Standard_Real& sint,
	277	Standard_Real& SigneSqrtDis) const
7fd59977	278	{
e2e0498b	279	const Standard_Real aRelTolp = 1.0+Epsilon(1.0), aRelTolm = 1.0-Epsilon(1.0);
3306fdd9	280
	281	// Infinitesimal step of increasing curve parameter. See comment below.
	282	const Standard_Real aDT = 100.0*Epsilon(DomainSup + DomainSup - DomainInf);
	283
7fd59977	284	Standard_Real Theta=theta;
	285	Standard_Boolean SecondSolution=Standard_False;
	286
3306fdd9	287	if ((Theta<DomainInf*aRelTolm) \|\|
	288	((Theta>DomainSup*aRelTolp) && (!TwoCurves)) \|\|
	289	(Theta>(DomainSup + DomainSup - DomainInf)*aRelTolp))
	290	{
7c65581d	291	SigneSqrtDis = 0.;
9775fa61	292	throw Standard_DomainError("IntAna_Curve::Domain");
7fd59977	293	}
7fd59977	294
3306fdd9	295	if (Abs(Theta - DomainSup) < aDT)
	296	{
	297	// Point of Null-discriminant.
	298	Theta = DomainSup;
	299	}
	300	else if (Theta>DomainSup)
	301	{
	302	Theta = DomainSup + DomainSup - Theta;
7fd59977	303	SecondSolution=Standard_True;
	304	}
	305
	306	Param1=Theta;
	307
	308	if(!TwoCurves) {
	309	SecondSolution=TakeZPositive;
	310	}
	311	//
	312	cost = Cos(Theta);
	313	sint = Sin(Theta);
3306fdd9	314	const Standard_Real aSin2t = Sin(Theta + Theta);
3306fdd9	315	const Standard_Real aCos2t = Cos(Theta + Theta);
7fd59977	316
7fd59977	317	A=Z2Cte+sint(Z2Sin+sintZ2SinSin)+cost(Z2Cos+costZ2CosCos)
3306fdd9	318	+ Z2CosSin*aSin2t;
7fd59977	319
3306fdd9	320	const Standard_Real aDA = costZ2Sin - sintZ2Cos +
	321	aSin2t*(Z2SinSin - Z2CosCos) +
	322	aCos2t(Z2CosSin Z2CosSin);
	323
7fd59977	324	B=Z1Cte+sint(Z1Sin+sintZ1SinSin)+cost(Z1Cos+costZ1CosCos)
3306fdd9	325	+ Z1CosSin*aSin2t;
	326
	327	const Standard_Real aDB = Z1Sincost - Z1Cossint +
	328	aSin2t*(Z1SinSin - Z1CosCos) +
	329	aCos2t*(Z1CosSin + Z1CosSin);
7fd59977	330
7fd59977	331	C=Z0Cte+sint(Z0Sin+sintZ0SinSin)+cost(Z0Cos+costZ0CosCos)
3306fdd9	332	+ Z0CosSin*aSin2t;
7fd59977	333
3306fdd9	334	const Standard_Real aDC = Z0Sincost - Z0Cossint +
	335	aSin2t*(Z0SinSin - Z0CosCos) +
	336	aCos2t*(Z0CosSin + Z0CosSin);
	337
	338	Standard_Real aDiscriminant = BB-4.0A*C;
7fd59977	339
3306fdd9	340	// We consider that infinitesimal dt = aDT.
	341	// Error of discriminant computation is equal to
	342	// (d(Disc)/dt)dt, where 1st derivative d(Disc)/dt = 2BaDB - 4(AaDC + CaDA).
	343
	344	const Standard_Real aTolD = 2.0aDTAbs(BaDB - 2.0(AaDC + CaDA));
7fd59977	345
3306fdd9	346	if (aDiscriminant < aTolD)
	347	aDiscriminant = 0.0;
	348
	349	if (Abs(A) <= Precision::PConfusion())
	350	{
	351	if (Abs(B) <= Precision::PConfusion())
	352	{
	353	Param2 = 0.0;
7fd59977	354	}
3306fdd9	355	else
	356	{
	357	Param2 = -C / B;
7fd59977	358	}
7fd59977	359	}
3306fdd9	360	else
	361	{
	362	SigneSqrtDis = (SecondSolution) ? Sqrt(aDiscriminant) : -Sqrt(aDiscriminant);
	363	Param2 = (-B + SigneSqrtDis) / (A + A);
7fd59977	364	}
	365	}
	366
	367	//=======================================================================
	368	//function : Value
	369	//purpose :
	370	//=======================================================================
3306fdd9	371	gp_Pnt IntAna_Curve::Value(const Standard_Real theta)
7fd59977	372	{
	373	Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
	374	//
	375	A=0.0; B=0.0; C=0.0;
	376	U=0.0; V=0.0;
	377	sint=0.0; cost=0.0;
	378	SigneSqrtDis=0.0;
	379	InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
	380	//-- checked the parameter U and Raises Domain Error if Error
	381	return(InternalValue(U,V));
	382	}
	383	//=======================================================================
	384	//function : D1u
	385	//purpose :
	386	//=======================================================================
3306fdd9	387	Standard_Boolean IntAna_Curve::D1u(const Standard_Real theta,
	388	gp_Pnt& Pt,
	389	gp_Vec& Vec)
7fd59977	390	{
	391	//-- Pour detecter le cas ou le calcul est impossible
	392	Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
	393	A=0.0; B=0.0; C=0.0;
	394	U=0.0; V=0.0;
	395	sint=0.0; cost=0.0;
	396	//
	397	InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
	398	//
	399	Pt = Value(theta);
3306fdd9	400	if (Abs(A)<1.0e-7 \|\| Abs(SigneSqrtDis)<1.0e-10) return(Standard_False);
7fd59977	401
	402
	403	//-- Approximation de la derivee (mieux que le calcul mathematique!)
3306fdd9	404	Standard_Real dtheta = (DomainSup - DomainInf)*1.0e-6;
7fd59977	405	Standard_Real theta2 = theta+dtheta;
3306fdd9	406	if ((theta2<DomainInf) \|\| ((theta2>DomainSup) && (!TwoCurves))
	407	\|\| (theta2>(DomainSup + DomainSup - DomainInf + 1.0e-14)))
	408	{
7fd59977	409	dtheta = -dtheta;
	410	theta2 = theta+dtheta;
	411	}
	412	gp_Pnt P2 = Value(theta2);
	413	dtheta = 1.0/dtheta;
	414	Vec.SetCoord((P2.X()-Pt.X())*dtheta,
	415	(P2.Y()-Pt.Y())*dtheta,
	416	(P2.Z()-Pt.Z())*dtheta);
	417
	418	return(Standard_True);
	419	}
	420	//=======================================================================
	421	//function : FindParameter
3306fdd9	422	//purpose : Projects P to the ALine. Returns the list of parameters as a results
	423	// of projection.
	424	// Sometimes aline can be self-intersected line (see bug #29807 where
	425	// ALine goes through the cone apex).
7fd59977	426	//=======================================================================
3306fdd9	427	void IntAna_Curve::FindParameter(const gp_Pnt& theP,
3306fdd9	428	TColStd_ListOfReal& theParams) const
7fd59977	429	{
3306fdd9	430	const Standard_Real aPIpPI = M_PI + M_PI,
	431	anEpsAng = 1.e-8,
	432	aSqTolPrecision=1.0e-8;
	433	Standard_Real aTheta = 0.0;
7fd59977	434	//
3306fdd9	435	switch (typequadric)
	436	{
	437	case GeomAbs_Cylinder:
7fd59977	438	{
3306fdd9	439	Standard_Real aZ;
3306fdd9	440	ElSLib::CylinderParameters(Ax3, RCyl, theP, aTheta, aZ);
7fd59977	441	}
7fd59977	442	break;
3306fdd9	443
3306fdd9	444	case GeomAbs_Cone:
7fd59977	445	{
3306fdd9	446	Standard_Real aZ;
3306fdd9	447	ElSLib::ConeParameters(Ax3, RCyl, Angle, theP, aTheta, aZ);
7fd59977	448	}
	449	break;
	450
3306fdd9	451	default:
3306fdd9	452	return;
7fd59977	453	}
7fd59977	454	//
3306fdd9	455	if (!firstbounded && (DomainInf > aTheta) && ((DomainInf - aTheta) <= anEpsAng))
	456	{
	457	aTheta = DomainInf;
7fd59977	458	}
3306fdd9	459	else if (!lastbounded && (aTheta > DomainSup) && ((aTheta - DomainSup) <= anEpsAng))
	460	{
	461	aTheta = DomainSup;
7fd59977	462	}
7fd59977	463	//
3306fdd9	464	if (aTheta < DomainInf)
	465	{
	466	aTheta = aTheta + aPIpPI;
7fd59977	467	}
3306fdd9	468	else if (aTheta > DomainSup)
	469	{
	470	aTheta = aTheta - aPIpPI;
7fd59977	471	}
7fd59977	472
3306fdd9	473	const Standard_Integer aMaxPar = 5;
	474	Standard_Real aParams[aMaxPar] = {DomainInf, DomainSup, aTheta,
	475	(TwoCurves)? DomainSup + DomainSup - aTheta : RealLast(),
	476	(TwoCurves) ? DomainSup + DomainSup - DomainInf : RealLast()};
7fd59977	477
3306fdd9	478	std::sort(aParams, aParams + aMaxPar - 1);
7fd59977	479
3306fdd9	480	for (Standard_Integer i = 0; i < aMaxPar; i++)
	481	{
	482	if (aParams[i] > myLastParameter)
	483	break;
	484
	485	if (aParams[i] < myFirstParameter)
	486	continue;
	487
	488	if (i && (aParams[i] - aParams[i - 1]) < Precision::PConfusion())
	489	continue;
	490
	491	Standard_Real U = 0.0, V= 0.0,
	492	A = 0.0, B = 0.0, C = 0.0,
	493	sint = 0.0, cost = 0.0, SigneSqrtDis = 0.0;
	494	InternalUVValue(aParams[i], U, V, A, B, C,
	495	cost, sint, SigneSqrtDis);
	496	const gp_Pnt aP(InternalValue(U, V));
	497	if (aP.SquareDistance(theP) < aSqTolPrecision)
	498	{
	499	theParams.Append(aParams[i]);
7fd59977	500	}
7fd59977	501	}
7fd59977	502	}
	503	//=======================================================================
	504	//function : InternalValue
	505	//purpose :
	506	//=======================================================================
3306fdd9	507	gp_Pnt IntAna_Curve::InternalValue(const Standard_Real U,
3306fdd9	508	const Standard_Real _V) const
7fd59977	509	{
	510	//-- cout<<" ["<<U<<","<<V<<"]";
	511	Standard_Real V = _V;
	512	if(V > 100000.0 ) { V= 100000.0; }
	513	if(V < -100000.0 ) { V=-100000.0; }
	514
	515	switch(typequadric) {
	516	case GeomAbs_Cone:
	517	{
	518	//------------------------------------------------
	519	//-- Parametrage : X = V * Cos(U) ---
	520	//-- Y = V * Sin(U) ---
	521	//-- Z = (V-RCyl) / Tan(SemiAngle)--
	522	//------------------------------------------------
	523	//-- Angle Vaut Cone.SemiAngle()
	524	return(ElSLib::ConeValue(U,(V-RCyl)/Sin(Angle),Ax3,RCyl,Angle));
	525	}
	526	break;
	527
	528	case GeomAbs_Cylinder:
	529	return(ElSLib::CylinderValue(U,V,Ax3,RCyl));
	530	case GeomAbs_Sphere:
	531	return(ElSLib::SphereValue(U,V,Ax3,RCyl));
	532	default:
	533	return(gp_Pnt(0.0,0.0,0.0));
	534	}
7fd59977	535	}
	536
	537	//
	538	//=======================================================================
	539	//function : SetDomain
	540	//purpose :
	541	//=======================================================================
3306fdd9	542	void IntAna_Curve::SetDomain(const Standard_Real theFirst,
3306fdd9	543	const Standard_Real theLast)
7fd59977	544	{
3306fdd9	545	if (theLast <= theFirst)
3306fdd9	546	{
9775fa61	547	throw Standard_DomainError("IntAna_Curve::Domain");
7fd59977	548	}
7fd59977	549	//
3306fdd9	550	myFirstParameter = theFirst;
3306fdd9	551	myLastParameter = theLast;
7fd59977	552	}