1 // Created on: 1992-06-30
2 // Created by: Laurent BUCHARD
3 // Copyright (c) 1992-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
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.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
18 #define No_Standard_RangeError
19 #define No_Standard_OutOfRange
22 //----------------------------------------------------------------------
23 //-- Differents constructeurs sont proposes qui correspondent aux
25 //-- A(Sin(Theta),Cos(Theta)) Z**2
26 //-- + B(Sin(Theta),Cos(Theta)) Z
27 //-- + C(Sin(Theta),Cos(Theta))
29 //-- Une Courbe est definie sur un domaine
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-)
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 //----------------------------------------------------------------------
43 #include <gp_Cone.hxx>
44 #include <gp_Cylinder.hxx>
48 #include <IntAna_Curve.hxx>
49 #include <math_DirectPolynomialRoots.hxx>
50 #include <Precision.hxx>
51 #include <Standard_DomainError.hxx>
53 //=======================================================================
54 //function : IntAna_Curve
56 //=======================================================================
57 IntAna_Curve::IntAna_Curve()
59 typequadric=GeomAbs_OtherSurface;
60 firstbounded=Standard_False;
61 lastbounded=Standard_False;
63 //=======================================================================
64 //function : SetConeQuadValues
65 //purpose : Description de l intersection Cone Quadrique
66 //=======================================================================
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)
85 Ax3 = Cone.Position();
86 RCyl = Cone.RefRadius();
88 Angle = Cone.SemiAngle();
89 Standard_Real UnSurTgAngle = 1.0/(Tan(Cone.SemiAngle()));
91 typequadric= GeomAbs_Cone;
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())
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
105 Z1Cte = 2.0*(UnSurTgAngle)*Qz;
112 Z2Cte = Qzz * UnSurTgAngle*UnSurTgAngle;
113 Z2Sin = (UnSurTgAngle+UnSurTgAngle)*Qyz;
114 Z2Cos = (UnSurTgAngle+UnSurTgAngle)*Qxz;
123 RestrictedInf = RestrictedSup = Standard_True; //-- Le Domaine est Borne
124 firstbounded = lastbounded = Standard_False;
126 myFirstParameter = DomainInf;
127 myLastParameter = (TwoCurves) ? DomainSup + DomainSup - DomainInf :
131 //=======================================================================
132 //function : SetCylinderQuadValues
133 //purpose : Description de l intersection Cylindre Quadrique
134 //=======================================================================
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)
153 Ax3 = Cyl.Position();
155 typequadric= GeomAbs_Cylinder;
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())
164 Z0CosCos = Qxx*RCyl*RCyl;
165 Z0SinSin = Qyy*RCyl*RCyl;
166 Z0CosSin = RCyl*RCyl*Qxy;
169 Z1Sin = RCylmul2*Qyz;
170 Z1Cos = RCylmul2*Qxz;
186 RestrictedInf = RestrictedSup = Standard_True;
187 firstbounded = lastbounded = Standard_False;
189 myFirstParameter = DomainInf;
190 myLastParameter = (TwoCurves) ? DomainSup + DomainSup - DomainInf :
194 //=======================================================================
197 //=======================================================================
198 Standard_Boolean IntAna_Curve::IsOpen() const
200 return(RestrictedInf && RestrictedSup);
203 //=======================================================================
206 //=======================================================================
207 void IntAna_Curve::Domain(Standard_Real& theFirst,
208 Standard_Real& theLast) const
210 if (RestrictedInf && RestrictedSup)
212 theFirst = myFirstParameter;
213 theLast = myLastParameter;
217 throw Standard_DomainError("IntAna_Curve::Domain");
220 //=======================================================================
221 //function : IsConstant
223 //=======================================================================
224 Standard_Boolean IntAna_Curve::IsConstant() const
226 //-- ??? Pas facile de decider a la seule vue des Param.
227 return(Standard_False);
230 //=======================================================================
231 //function : IsFirstOpen
233 //=======================================================================
234 Standard_Boolean IntAna_Curve::IsFirstOpen() const
236 return(firstbounded);
239 //=======================================================================
240 //function : IsLastOpen
242 //=======================================================================
243 Standard_Boolean IntAna_Curve::IsLastOpen() const
247 //=======================================================================
248 //function : SetIsFirstOpen
250 //=======================================================================
251 void IntAna_Curve::SetIsFirstOpen(const Standard_Boolean Flag)
256 //=======================================================================
257 //function : SetIsLastOpen
259 //=======================================================================
260 void IntAna_Curve::SetIsLastOpen(const Standard_Boolean Flag)
265 //=======================================================================
266 //function : InternalUVValue
268 //=======================================================================
269 void IntAna_Curve::InternalUVValue(const Standard_Real theta,
270 Standard_Real& Param1,
271 Standard_Real& Param2,
277 Standard_Real& SigneSqrtDis) const
279 const Standard_Real aRelTolp = 1.0+Epsilon(1.0), aRelTolm = 1.0-Epsilon(1.0);
281 // Infinitesimal step of increasing curve parameter. See comment below.
282 const Standard_Real aDT = 100.0*Epsilon(DomainSup + DomainSup - DomainInf);
284 Standard_Real Theta=theta;
285 Standard_Boolean SecondSolution=Standard_False;
287 if ((Theta<DomainInf*aRelTolm) ||
288 ((Theta>DomainSup*aRelTolp) && (!TwoCurves)) ||
289 (Theta>(DomainSup + DomainSup - DomainInf)*aRelTolp))
292 throw Standard_DomainError("IntAna_Curve::Domain");
295 if (Abs(Theta - DomainSup) < aDT)
297 // Point of Null-discriminant.
300 else if (Theta>DomainSup)
302 Theta = DomainSup + DomainSup - Theta;
303 SecondSolution=Standard_True;
309 SecondSolution=TakeZPositive;
314 const Standard_Real aSin2t = Sin(Theta + Theta);
315 const Standard_Real aCos2t = Cos(Theta + Theta);
317 A=Z2Cte+sint*(Z2Sin+sint*Z2SinSin)+cost*(Z2Cos+cost*Z2CosCos)
320 const Standard_Real aDA = cost*Z2Sin - sint*Z2Cos +
321 aSin2t*(Z2SinSin - Z2CosCos) +
322 aCos2t*(Z2CosSin * Z2CosSin);
324 B=Z1Cte+sint*(Z1Sin+sint*Z1SinSin)+cost*(Z1Cos+cost*Z1CosCos)
327 const Standard_Real aDB = Z1Sin*cost - Z1Cos*sint +
328 aSin2t*(Z1SinSin - Z1CosCos) +
329 aCos2t*(Z1CosSin + Z1CosSin);
331 C=Z0Cte+sint*(Z0Sin+sint*Z0SinSin)+cost*(Z0Cos+cost*Z0CosCos)
334 const Standard_Real aDC = Z0Sin*cost - Z0Cos*sint +
335 aSin2t*(Z0SinSin - Z0CosCos) +
336 aCos2t*(Z0CosSin + Z0CosSin);
338 Standard_Real aDiscriminant = B*B-4.0*A*C;
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 = 2*B*aDB - 4*(A*aDC + C*aDA).
344 const Standard_Real aTolD = 2.0*aDT*Abs(B*aDB - 2.0*(A*aDC + C*aDA));
346 if (aDiscriminant < aTolD)
349 if (Abs(A) <= Precision::PConfusion())
351 if (Abs(B) <= Precision::PConfusion())
362 SigneSqrtDis = (SecondSolution) ? Sqrt(aDiscriminant) : -Sqrt(aDiscriminant);
363 Param2 = (-B + SigneSqrtDis) / (A + A);
367 //=======================================================================
370 //=======================================================================
371 gp_Pnt IntAna_Curve::Value(const Standard_Real theta)
373 Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
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));
383 //=======================================================================
386 //=======================================================================
387 Standard_Boolean IntAna_Curve::D1u(const Standard_Real theta,
391 //-- Pour detecter le cas ou le calcul est impossible
392 Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
397 InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
400 if (Abs(A)<1.0e-7 || Abs(SigneSqrtDis)<1.0e-10) return(Standard_False);
403 //-- Approximation de la derivee (mieux que le calcul mathematique!)
404 Standard_Real dtheta = (DomainSup - DomainInf)*1.0e-6;
405 Standard_Real theta2 = theta+dtheta;
406 if ((theta2<DomainInf) || ((theta2>DomainSup) && (!TwoCurves))
407 || (theta2>(DomainSup + DomainSup - DomainInf + 1.0e-14)))
410 theta2 = theta+dtheta;
412 gp_Pnt P2 = Value(theta2);
414 Vec.SetCoord((P2.X()-Pt.X())*dtheta,
415 (P2.Y()-Pt.Y())*dtheta,
416 (P2.Z()-Pt.Z())*dtheta);
418 return(Standard_True);
420 //=======================================================================
421 //function : FindParameter
422 //purpose : Projects P to the ALine. Returns the list of parameters as a results
424 // Sometimes aline can be self-intersected line (see bug #29807 where
425 // ALine goes through the cone apex).
426 //=======================================================================
427 void IntAna_Curve::FindParameter(const gp_Pnt& theP,
428 TColStd_ListOfReal& theParams) const
430 const Standard_Real aPIpPI = M_PI + M_PI,
432 aSqTolPrecision=1.0e-8;
433 Standard_Real aTheta = 0.0;
437 case GeomAbs_Cylinder:
440 ElSLib::CylinderParameters(Ax3, RCyl, theP, aTheta, aZ);
447 ElSLib::ConeParameters(Ax3, RCyl, Angle, theP, aTheta, aZ);
455 if (!firstbounded && (DomainInf > aTheta) && ((DomainInf - aTheta) <= anEpsAng))
459 else if (!lastbounded && (aTheta > DomainSup) && ((aTheta - DomainSup) <= anEpsAng))
464 if (aTheta < DomainInf)
466 aTheta = aTheta + aPIpPI;
468 else if (aTheta > DomainSup)
470 aTheta = aTheta - aPIpPI;
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()};
478 std::sort(aParams, aParams + aMaxPar - 1);
480 for (Standard_Integer i = 0; i < aMaxPar; i++)
482 if (aParams[i] > myLastParameter)
485 if (aParams[i] < myFirstParameter)
488 if (i && (aParams[i] - aParams[i - 1]) < Precision::PConfusion())
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)
499 theParams.Append(aParams[i]);
503 //=======================================================================
504 //function : InternalValue
506 //=======================================================================
507 gp_Pnt IntAna_Curve::InternalValue(const Standard_Real U,
508 const Standard_Real _V) const
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; }
515 switch(typequadric) {
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));
528 case GeomAbs_Cylinder:
529 return(ElSLib::CylinderValue(U,V,Ax3,RCyl));
531 return(ElSLib::SphereValue(U,V,Ax3,RCyl));
533 return(gp_Pnt(0.0,0.0,0.0));
538 //=======================================================================
539 //function : SetDomain
541 //=======================================================================
542 void IntAna_Curve::SetDomain(const Standard_Real theFirst,
543 const Standard_Real theLast)
545 if (theLast <= theFirst)
547 throw Standard_DomainError("IntAna_Curve::Domain");
550 myFirstParameter = theFirst;
551 myLastParameter = theLast;