1 // Created on: 1992-06-29
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.
19 #include <Standard_Stream.hxx>
22 #define No_Standard_RangeError
23 #define No_Standard_OutOfRange
26 //======================================================================
27 //== I n t e r s e c t i o n C O N E Q U A D R I Q U E
28 //== C Y L I N D R E Q U A D R I Q U E
29 //======================================================================
33 #include <gp_Cone.hxx>
34 #include <gp_Cylinder.hxx>
36 #include <IntAna_Curve.hxx>
37 #include <IntAna_IntQuadQuad.hxx>
38 #include <IntAna_Quadric.hxx>
39 #include <math_TrigonometricFunctionRoots.hxx>
40 #include <Standard_DomainError.hxx>
41 #include <Standard_OutOfRange.hxx>
42 #include <StdFail_NotDone.hxx>
44 //=======================================================================
45 //class : TrigonometricRoots
46 //purpose: Classe Interne (Donne des racines classees d un polynome trigo)
47 //======================================================================
48 class TrigonometricRoots {
51 Standard_Real Roots[4];
52 Standard_Boolean done;
53 Standard_Integer NbRoots;
54 Standard_Boolean infinite_roots;
57 TrigonometricRoots(const Standard_Real CC,
58 const Standard_Real SC,
59 const Standard_Real C,
60 const Standard_Real S,
61 const Standard_Real Cte,
62 const Standard_Real Binf,
63 const Standard_Real Bsup);
65 Standard_Boolean IsDone() {
69 Standard_Boolean IsARoot(Standard_Real u) {
71 Standard_Real aEps=RealEpsilon();
72 Standard_Real PIpPI = M_PI + M_PI;
74 for(i=0 ; i<NbRoots; ++i) {
75 if(Abs(u - Roots[i])<=aEps) {
78 if(Abs(u - Roots[i]-PIpPI)<=aEps) {
82 return Standard_False;
85 Standard_Integer NbSolutions() {
87 throw StdFail_NotDone();
92 Standard_Boolean InfiniteRoots() {
94 throw StdFail_NotDone();
96 return infinite_roots;
99 Standard_Real Value(const Standard_Integer n) {
100 if((!done)||(n>NbRoots)) {
101 throw StdFail_NotDone();
106 //=======================================================================
107 //function : TrigonometricRoots::TrigonometricRoots
109 //=======================================================================
110 TrigonometricRoots::TrigonometricRoots(const Standard_Real CC,
111 const Standard_Real SC,
112 const Standard_Real C,
113 const Standard_Real S,
114 const Standard_Real Cte,
115 const Standard_Real Binf,
116 const Standard_Real Bsup)
118 Standard_Integer i, j, SvNbRoots;
119 Standard_Boolean Triee;
120 Standard_Real PIpPI = M_PI + M_PI;
124 //-- F= AA*CN*CN+2*BB*CN*SN+CC*CN+DD*SN+EE;
125 math_TrigonometricFunctionRoots MTFR(CC, SC, C, S, Cte, Binf, Bsup);
131 if(MTFR.InfiniteRoots()) {
132 infinite_roots=Standard_True;
136 NbRoots=MTFR.NbSolutions();
138 for(i=0; i<NbRoots; ++i) {
139 Roots[i]=MTFR.Value(i+1);
148 //-- La recherche directe donne n importe quoi.
150 for(i=0; i<SvNbRoots; ++i) {
152 Standard_Real co=cos(Roots[i]);
153 Standard_Real si=sin(Roots[i]);
154 y=co*(CC*co + (SC+SC)*si + C) + S*si + Cte;
157 return; //-- le 1er avril 98
163 for(i=1,j=0; i<SvNbRoots; ++i,++j) {
164 if(Roots[i]<Roots[j]) {
165 Triee=Standard_False;
166 Standard_Real t=Roots[i];
174 infinite_roots=Standard_False;
176 if(!NbRoots) {//--!!!!! Detection du cas Pol = Cte ( 1e-50 ) !!!!
177 if((Abs(CC) + Abs(SC) + Abs(C) + Abs(S)) < 1e-10) {
178 if(Abs(Cte) < 1e-10) {
179 infinite_roots=Standard_True;
184 //=======================================================================
185 //class : MyTrigonometricFunction
187 // Classe interne : Donne Value et Derivative sur un polynome
189 //======================================================================
190 class MyTrigonometricFunction {
193 Standard_Real CC,SS,SC,S,C,Cte;
197 MyTrigonometricFunction(const Standard_Real xCC,
198 const Standard_Real xSS,
199 const Standard_Real xSC,
200 const Standard_Real xC,
201 const Standard_Real xS,
202 const Standard_Real xCte) {
211 Standard_Real Value(const Standard_Real& U) {
212 Standard_Real sinus, cosinus, aRet;
216 aRet= CC*cosinus*cosinus +
218 2.0*(sinus*(SC*cosinus+S)+cosinus*C)+
224 Standard_Real Derivative(const Standard_Real& U) {
225 Standard_Real sinus, cosinus;
230 return(2.0*((sinus*cosinus)*(SS-CC)
233 +SC*(cosinus*cosinus-sinus*sinus)));
238 //=======================================================================
239 //function : IntAna_IntQuadQuad::IntAna_IntQuadQuad
240 //purpose : C o n s t r u c t e u r v i d e
241 //=======================================================================
242 IntAna_IntQuadQuad::IntAna_IntQuadQuad(void) {
244 identical = Standard_False;
248 myEpsilon=0.00000001;
249 myEpsilonCoeffPolyNull=0.00000001;
251 //=======================================================================
252 //function : IntAna_IntQuadQuad::IntAna_IntQuadQuad
253 //purpose : I n t e r s e c t i o n C y l i n d r e Q u a d r i q u e
254 //=======================================================================
255 IntAna_IntQuadQuad::IntAna_IntQuadQuad(const gp_Cylinder& Cyl,
256 const IntAna_Quadric& Quad,
257 const Standard_Real Tol) {
259 myEpsilon=0.00000001;
260 myEpsilonCoeffPolyNull=0.00000001;
261 Perform(Cyl,Quad,Tol);
263 //=======================================================================
265 //purpose : I n t e r s e c t i o n C y l i n d r e Q u a d r i q u e
266 //=======================================================================
267 void IntAna_IntQuadQuad::Perform(const gp_Cylinder& Cyl,
268 const IntAna_Quadric& Quad,
271 done = Standard_True;
272 identical= Standard_False;
276 Standard_Boolean UN_SEUL_Z_PAR_THETA, DEUX_Z_PAR_THETA,
277 Z_POSITIF, Z_INDIFFERENT, Z_NEGATIF;
279 UN_SEUL_Z_PAR_THETA=Standard_False;
280 DEUX_Z_PAR_THETA=Standard_True;
281 Z_POSITIF=Standard_True;
282 Z_INDIFFERENT=Standard_True;
283 Z_NEGATIF=Standard_False;
285 Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, aRealEpsilon, RCyl, R2;
286 Standard_Real PIpPI = M_PI + M_PI;
288 for(Standard_Integer raz = 0 ; raz < myNbMaxCurves ; raz++) {
289 previouscurve[raz] = nextcurve[raz] = 0;
293 aRealEpsilon=RealEpsilon();
294 //----------------------------------------------------------------------
295 //-- Coefficients de la quadrique :
297 //-- Qxx x + Qyy y + Qzz z + 2 ( Qxy x y + Qxz x z + Qyz y z )
298 //-- + 2 ( Qx x + Qy y + Qz z ) + Q1
299 //----------------------------------------------------------------------
300 Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1);
302 //----------------------------------------------------------------------
303 //-- Ecriture des Coefficients de la Quadrique dans le repere liee
305 //-- Cyl.Position() -> Ax2
306 //----------------------------------------------------------------------
307 Quad.NewCoefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,Cyl.Position());
309 //----------------------------------------------------------------------
310 //-- Parametrage du Cylindre Cyl :
311 //-- X = Rcyl * Cos(Theta)
312 //-- Y = Rcyl * Sin(Theta)
314 //----------------------------------------------------------------------
315 //-- Donne une Equation de la forme :
316 //-- F(Sin(Theta),Cos(Theta),ZCyl) = 0
317 //-- (Equation du second degre en ZCyl)
318 //-- ZCyl**2 CoeffZ2(Theta) + ZCyl CoeffZ1(Theta) + CoeffZ0(Theta)
319 //----------------------------------------------------------------------
320 //-- CoeffZ0 = Q1 + 2*Qx*RCyl*Cos[Theta] + Qxx*RCyl^2*Cos[Theta]^2
321 //-- 2*RCyl*Sin[Theta]* ( Qy + Qxy*RCyl*Cos[Theta]);
322 //-- Qyy*RCyl^2*Sin[Theta]^2;
323 //-- CoeffZ1 =2.0 * ( Qz + RCyl*(Qxz*Cos[Theta] + Sin[Theta]*Qyz)) ;
325 //-- !!!! Attention , si on norme sur Qzz pour detecter le cas 1er degre
326 //----------------------------------------------------------------------
327 //-- On Cherche Les Racines en Theta du discriminant de cette equation :
328 //-- DIS(Theta) = C_1 + C_SS Sin()**2 + C_CC Cos()**2 + 2 C_SC Sin() Cos()
329 //-- + 2 C_S Sin() + 2 C_C Cos()
331 //-- Si Qzz = 0 Alors On Resout Z=Fct(Theta) sur le domaine de Theta
332 //----------------------------------------------------------------------
334 if(Abs(Qzz)<myEpsilonCoeffPolyNull) {
335 done=Standard_False; //-- le 12 mars 98
339 //----------------------------------------------------------------------
340 //--- Cas ou F(Z,Theta) est du second degre en Z ----
341 //----------------------------------------------------------------------
344 Standard_Real C_1 = Qz*Qz - Qzz*Q1;
345 Standard_Real C_SS = R2 * ( Qyz*Qyz - Qyy*Qzz );
346 Standard_Real C_CC = R2 * ( Qxz*Qxz - Qxx*Qzz );
347 Standard_Real C_S = RCyl * ( Qyz*Qz - Qy*Qzz );
348 Standard_Real C_C = RCyl * ( Qxz*Qz - Qx*Qzz );
349 Standard_Real C_SC = R2 * ( Qxz*Qyz - Qxy*Qzz );
351 MyTrigonometricFunction MTF(C_CC,C_SS,C_SC,C_C,C_S,C_1);
352 TrigonometricRoots PolDIS(C_CC-C_SS,
356 C_1+C_SS, 0., PIpPI);
357 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
358 if(PolDIS.IsDone()==Standard_False) {
362 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
363 if(PolDIS.InfiniteRoots()) {
364 TheCurve[0].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
368 TheCurve[1].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
376 //---------------------------------------------------------------
377 //-- La Recherche des Zero de DIS a reussi
378 //---------------------------------------------------------------
379 Standard_Integer nbsolDIS=PolDIS.NbSolutions();
381 //--------------------------------------------------------------
382 //-- Le Discriminant a un signe constant :
384 //-- Si Positif ---> 2 Courbes
385 //-- Sinon ---> Pas de solution
386 //--------------------------------------------------------------
387 if(MTF.Value(M_PI) >= -aRealEpsilon) {
389 TheCurve[0].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
393 TheCurve[1].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
401 //------------------------------------------------------------
402 //-- Le Discriminant est toujours Negatif
403 //------------------------------------------------------------
408 //------------------------------------------------------------
409 //-- Le Discriminant a des racines
410 //-- (Le Discriminant est une fonction periodique donc ... )
411 //------------------------------------------------------------
413 //------------------------------------------------------------
414 //-- Point de Tangence pour ce Theta Solution
415 //-- Si Signe de Discriminant >=0 pour tout Theta
417 //-- Courbe Solution pour la partie Positive
418 //-- entre les 2 racines ( Ici Tout le domaine )
419 //-- Sinon Seulement un point Tangent
420 //------------------------------------------------------------
421 if(MTF.Value(PolDIS.Value(1)+M_PI) >= -aRealEpsilon ) {
422 //------------------------------------------------------------
423 //-- On a Un Point de Tangence + Une Courbe Solution
424 //------------------------------------------------------------
425 TheCurve[0].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
429 TheCurve[1].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
437 //------------------------------------------------------------
438 //-- On a simplement un Point de tangence
439 //------------------------------------------------------------
440 //--Standard_Real Theta = PolDIS(1);
441 //--Standard_Real SecPar= -0.5 * MTFZ1.Value(Theta) / MTFZ2.Value(Theta);
442 //--Thepoints[Nbpoints] = ElSLib::CylinderValue(Theta,SecPar,Cyl.Position(),Cyl.Radius());
448 //------------------------------------------------------------
449 //-- On detecte : Les racines double
450 //-- : Les Zones Positives [Modulo 2 PI]
451 //-- Les Courbes solutions seront obtenues pour les valeurs
452 //-- de Theta ou Discriminant(Theta) > 0 (>=0? en limite)
453 //-- Par la resolution de l equation implicite avec Theta fixe
454 //------------------------------------------------------------
455 //-- Si tout est solution, Alors on est sur une iso
456 //-- Ce cas devrait etre traite en amont
457 //------------------------------------------------------------
458 //-- On construit la fonction permettant connaissant un Theta
459 //-- de calculer les Z+ et Z- Correspondants.
460 //-------------------------------------------------------------
462 //-------------------------------------------------------------
463 //-- Calcul des Intervalles ou Discriminant >=0
464 //-- On a au plus nbsol intervalles ( en fait 2 )
465 //-- (1,2) (2,3) .. (nbsol,1+2PI)
466 //-------------------------------------------------------------
468 Standard_Real Theta1, Theta2, Theta3, autrepar, qwet;
469 Standard_Boolean UnPtTg = Standard_False;
473 for(i=1; i<=nbsolDIS ; i++) {
474 Theta1=PolDIS.Value(i);
475 Theta2=(i<nbsolDIS)? PolDIS.Value(i+1) : (PolDIS.Value(1)+PIpPI);
476 //----------------------------------------------------------------
477 //-- On detecte les racines doubles
478 //-- Si il n y a que 2 racines alors on prend tout l interval
479 //----------------------------------------------------------------
480 if(Abs(Theta2-Theta1)<=aRealEpsilon) {
481 UnPtTg = Standard_True;
487 qwet=MTF.Value(autrepar);
489 TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
490 myEpsilon,Theta1,Theta1+PIpPI,
494 TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
495 myEpsilon,Theta1,Theta1+PIpPI,
504 for(i=1; UnPtTg == (Standard_False) && (i<=nbsolDIS) ; i++) {
505 Theta1=PolDIS.Value(i);
506 Theta2=(i<nbsolDIS)? PolDIS.Value(i+1) : (PolDIS.Value(1)+PIpPI);
507 //----------------------------------------------------------------
508 //-- On detecte les racines doubles
509 //-- Si il n y a que 2 racines alors on prend tout l interval
510 //----------------------------------------------------------------
511 if(Abs(Theta2-Theta1)<=1e-12) {
512 //-- cout<<"\n####### Un Point de Tangence en Theta = "<<Theta1<<endl;
515 else { //-- On evite les pbs numeriques (Tout est du meme signe entre les racines)
516 qwet=MTF.Value(0.5*(Theta1+Theta2))
517 +MTF.Value(0.4*Theta1+0.6*Theta2)
518 +MTF.Value(0.6*Theta1+0.4*Theta2);
520 //------------------------------------------------------------
521 //-- On est positif a partir de Theta1
522 //-- L intervalle Theta1,Theta2 est retenu
523 //------------------------------------------------------------
525 //-- le 8 octobre 1997 :
526 //-- PB: Racine en 0 pi/2 pi/2 et 2pi
527 //-- On cree 2 courbes : 0 -> pi/2 2zpartheta
528 //-- pi/2 -> 2pi 2zpartheta
530 //-- la courbe 0 pi/2 passe par le double pt de tangence pi/2
531 //-- il faut donc couper cette courbe en 2
533 Theta3 = ((i+1)<nbsolDIS)? PolDIS.Value(i+2) : (PolDIS.Value(1)+PIpPI);
535 if((Theta3-Theta2)<5.e-8) {
537 TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
538 myEpsilon,Theta1,Theta2,
542 TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
543 myEpsilon,Theta1,Theta2,
549 TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
550 myEpsilon,Theta1,Theta2,
556 }//else { //-- On evite les pbs numerique ...
557 } //for(i=1; UnPtTg == (Standard_False) && (i<=nbsolDIS) ; i++) {
565 //=======================================================================
566 //function :IntAna_IntQuadQuad::IntAna_IntQuadQuad
568 //=======================================================================
569 IntAna_IntQuadQuad::IntAna_IntQuadQuad(const gp_Cone& Cone,
570 const IntAna_Quadric& Quad,
571 const Standard_Real Tol)
574 myEpsilon=0.00000001;
575 myEpsilonCoeffPolyNull=0.00000001;
576 Perform(Cone,Quad,Tol);
578 //=======================================================================
581 //=======================================================================
582 void IntAna_IntQuadQuad::Perform(const gp_Cone& Cone,
583 const IntAna_Quadric& Quad,
587 Standard_Boolean UN_SEUL_Z_PAR_THETA,
588 Z_POSITIF, Z_NEGATIF;
590 UN_SEUL_Z_PAR_THETA=Standard_False;
591 Z_POSITIF=Standard_True;
592 Z_NEGATIF=Standard_False;
595 Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1;
596 Standard_Real Theta1, Theta2, TgAngle;
597 Standard_Real PIpPI = M_PI + M_PI;
600 identical = Standard_False;
604 for(i=0 ; i<myNbMaxCurves ; ++i) {
609 Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1);
611 gp_Ax3 tAx3(Cone.Position());
612 tAx3.SetLocation(Cone.Apex());
613 Quad.NewCoefficients(Qxx,Qyy,Qzz,
618 TgAngle=1./(Tan(Cone.SemiAngle()));
620 // The parametrization of the Cone
622 // x= z*tan(beta)*cos(t)
623 // y= z*tan(beta)*sin(t)
626 // The intersection curves are defined by the equation
629 // f(z,t)= A(t)*z + B(t)*z + C(t)=0
632 // 1. Try to solve A(t)=0 -> PolZ2
634 Standard_Integer nbsol, nbsol1, nbsolZ2;
635 Standard_Real Z2CC, Z2SS, Z2Cte, Z2SC, Z2C, Z2S;
636 Standard_Real Z1CC, Z1SS, Z1Cte, Z1SC, Z1C, Z1S;
637 Standard_Real C_1,C_SS, C_CC, C_S, C_C, C_SC;
641 Z2Cte= Qzz * TgAngle*TgAngle;
646 TrigonometricRoots PolZ2(Z2CC-Z2SS,Z2SC,Z2C+Z2C,Z2S+Z2S,Z2Cte+Z2SS,0.,PIpPI);
647 if(!PolZ2.IsDone()) {
652 //MyTrigonometricFunction MTF2(Z2CC,Z2SS,Z2SC,Z2C,Z2S,Z2Cte);
653 nbsolZ2 = PolZ2.NbSolutions();
655 // 2. Try to solve B(t)=0 -> PolZ1
657 Z1Cte = 2.*(TgAngle)*Qz;
664 TrigonometricRoots PolZ1(Z1CC-Z1SS,Z1SC, Z1C+Z1C,Z1S+Z1S, Z1Cte+Z1SS,0.,PIpPI);
665 if(!PolZ1.IsDone()) {
669 MyTrigonometricFunction MTFZ1(Z1CC,Z1SS,Z1SC,Z1C,Z1S,Z1Cte);
671 nbsol1=PolZ1.NbSolutions();
672 if(PolZ2.InfiniteRoots()) { // i.e A(t)=0 for any t
673 if(!PolZ1.InfiniteRoots()) {
675 //-- B(t)*z + C(t) = 0 avec C(t) != 0
676 TheCurve[0].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
685 for(ii=1; ii<= nbsol1 ; ++ii) {
686 Standard_Real Theta=PolZ1.Value(ii);
687 if(Abs(MTFZ1.Value(Theta))<=myEpsilon) {
688 //-- Une droite Solution Z= -INF -> +INF pour Theta
689 //-- cout<<"######## Droite Solution Pour Theta = "<<Theta<<endl;
692 //-- cout<<"\n#### _+_+_+_+_+ CAS A(t) Z + B = 0 avec A et B ==0 "<<endl;
699 if(Abs(Q1)<=myEpsilon) {
704 //-- Pas de Solutions
713 //-- f(z,t)=A(t)*z + B(t)*z + C(t)=0 avec A n est pas toujours nul
716 // 3. Try to explore s.c. Discriminant: D(t)=B(t)-4*A(t)*C(t) => Pol
718 // The function D(t) is just a formula that has sense for quadratic
720 // For cases when A(t)=0 (say at t=ti, t=tj. etc) the equation
723 // f(z,t)=B(t)*z + C(t)=0, where B(t)!=0,
725 // and the D(t) is nonsense for it.
727 C_1 = TgAngle*TgAngle * ( Qz * Qz - Qzz * Q1);
728 C_SS = Qy * Qy - Qyy * Q1;
729 C_CC = Qx * Qx - Qxx * Q1;
730 C_S = TgAngle*( Qy * Qz - Qyz * Q1);
731 C_C = TgAngle*( Qx * Qz - Qxz * Q1);
732 C_SC = Qx * Qy - Qxy * Q1;
734 TrigonometricRoots Pol(C_CC-C_SS, C_SC, C_C+C_C,C_S+C_S, C_1+C_SS,0.,PIpPI);
740 nbsol=Pol.NbSolutions();
741 MyTrigonometricFunction MTF(C_CC,C_SS,C_SC,C_C,C_S,C_1);
743 if(Pol.InfiniteRoots()) {
745 // f(z,t)= (z(t)-zo(t)) = 0 Discriminant(t)=0 pour tout t
747 TheCurve[0].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
749 UN_SEUL_Z_PAR_THETA, Z_POSITIF);
750 TheCurve[1].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
752 UN_SEUL_Z_PAR_THETA, Z_NEGATIF);
758 // f(z,t)=A(t)*z + B(t)*z + C(t) Discriminant(t) != 0
760 if(!nbsol && (MTF.Value(M_PI)<0.) ) {
761 //-- Discriminant signe constant negatif
766 //-- On Traite le cas : Discriminant(t) > 0 pour tout t en simulant 1
767 // racine double en 0
768 Standard_Boolean DiscriminantConstantPositif = Standard_False;
771 DiscriminantConstantPositif = Standard_True;
773 //----------------------------------------------------------------------
774 //-- Le discriminant admet au moins une racine ( -> Point de Tangence )
775 //-- ou est constant positif.
776 //----------------------------------------------------------------------
777 for(i=1; i<=nbsol; ++i) {
778 if(DiscriminantConstantPositif) {
780 Theta2 = PIpPI-myEpsilon;
784 Theta2=(i<nbsol)? Pol.Value(i+1) : (Pol.Value(1)+PIpPI);
787 if(Abs(Theta2-Theta1)<=myEpsilon){
792 Standard_Real qwet =MTF.Value(0.5*(Theta1+Theta2))+
793 MTF.Value(0.4*Theta1+0.6*Theta2)+
794 MTF.Value(0.6*Theta1+0.4*Theta2);
798 //---------------------------------------------------------------------
799 //-- On a des Solutions entre Theta1 et Theta 2
800 //---------------------------------------------------------------------
802 //---------------------------------------------------------------------
803 //-- On Subdivise si necessaire Theta1-->Theta2
804 //-- en Theta1-->t1 t1--->t2 ... tn--->Theta2
805 //-- ou t1,t2,...tn sont des racines de A(t)
807 //-- Seule la courbe Z_NEGATIF est affectee
808 //----------------------------------------------------------------------
809 Standard_Boolean RacinesdePolZ2DansTheta1Theta2;
813 //nbsolZ2=PolZ2.NbSolutions();
814 RacinesdePolZ2DansTheta1Theta2=Standard_False;
815 for(i2=1; i2<=nbsolZ2 && !RacinesdePolZ2DansTheta1Theta2; ++i2) {
817 if(r>Theta1 && r<Theta2) {
818 RacinesdePolZ2DansTheta1Theta2=Standard_True;
822 if(r>Theta1 && r<Theta2){
823 RacinesdePolZ2DansTheta1Theta2=Standard_True;
828 if(!RacinesdePolZ2DansTheta1Theta2) {
829 //--------------------------------------------------------------------
830 //-- Pas de Branches Infinies
831 TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,myEpsilon,
833 UN_SEUL_Z_PAR_THETA,Z_POSITIF);
835 TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,myEpsilon,
843 Standard_Boolean NoChanges;
844 Standard_Real NewMin, NewMax, to;
848 NoChanges=Standard_True;
850 for(i2=1; i2<= (nbsolZ2+nbsolZ2) ; ++i2) {
852 to=PolZ2.Value(i2-nbsolZ2) + PIpPI;
858 // to is value of the parameter t where A(t)=0, i.e.
859 // f(z,to)=B(to)*z + C(to)=0, B(to)!=0.
861 // z=-C(to)/B(to) is one and only
862 // solution inspite of the fact that D(t)>0 for any value of t.
864 if(to<NewMax && to>NewMin) {
865 //-----------------------------------------------------------------
866 //-- On coupe au moins une fois le domaine Theta1 Theta2
867 //-----------------------------------------------------------------
868 NoChanges=Standard_False;
869 TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
871 UN_SEUL_Z_PAR_THETA, Z_NEGATIF);
874 TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
876 UN_SEUL_Z_PAR_THETA, Z_POSITIF);
877 //------------------------------------------------------------
879 //-- Si B < 0 Alors Infini sur Z+
880 //-- Sinon Infini sur Z-
881 //------------------------------------------------------------
882 if(PolZ2.IsARoot(NewMin)) {
883 if(MTFZ1.Value(NewMin) < 0.) {
884 TheCurve[NbCurves].SetIsFirstOpen(Standard_True);
887 TheCurve[NbCurves-1].SetIsFirstOpen(Standard_True);
890 if(MTFZ1.Value(to) < 0.) {
891 TheCurve[NbCurves].SetIsLastOpen(Standard_True);
894 TheCurve[NbCurves-1].SetIsLastOpen(Standard_True);
896 //------------------------------------------------------------
899 }//if(to<NewMax && to>NewMin)
900 }// for(i2=1; i2<= (nbsolZ2+nbsolZ2) ; ++i2)
903 TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
905 UN_SEUL_Z_PAR_THETA, Z_NEGATIF);
907 TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
909 UN_SEUL_Z_PAR_THETA, Z_POSITIF);
910 //------------------------------------------------------------
912 //-- Si B < 0 Alors Infini sur Z+
913 //-- Sinon Infini sur Z-
914 //------------------------------------------------------------
915 if(PolZ2.IsARoot(Theta1)) {
916 if(MTFZ1.Value(Theta1) < 0.) {
917 TheCurve[NbCurves].SetIsFirstOpen(Standard_True);
920 TheCurve[NbCurves-1].SetIsFirstOpen(Standard_True);
923 if(PolZ2.IsARoot(Theta2)) {
924 if(MTFZ1.Value(Theta2) < 0.) {
925 TheCurve[NbCurves].SetIsLastOpen(Standard_True);
928 TheCurve[NbCurves-1].SetIsLastOpen(Standard_True);
931 //------------------------------------------------------------
935 TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
937 UN_SEUL_Z_PAR_THETA, Z_NEGATIF);
939 TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
941 UN_SEUL_Z_PAR_THETA, Z_POSITIF);
942 //------------------------------------------------------------
944 //-- Si B < 0 Alors Infini sur Z+
945 //-- Sinon Infini sur Z-
946 //------------------------------------------------------------
947 if(PolZ2.IsARoot(NewMin)) {
948 if(MTFZ1.Value(NewMin) < 0.) {
949 TheCurve[NbCurves].SetIsFirstOpen(Standard_True);
952 TheCurve[NbCurves-1].SetIsFirstOpen(Standard_True);
955 if(PolZ2.IsARoot(Theta2)) {
956 if(MTFZ1.Value(Theta2) < 0.) {
957 TheCurve[NbCurves].SetIsLastOpen(Standard_True);
960 TheCurve[NbCurves-1].SetIsLastOpen(Standard_True);
963 //------------------------------------------------------------
968 }//for(i=1; i<=nbsol ; i++) {
970 InternalSetNextAndPrevious();
972 //=======================================================================
973 //function :InternalSetNextAndPrevious
975 //-- Raccordement des courbes bout a bout
976 //-- - Utilisation du champ : IsFirstOpen
978 //-- Pas de verification si ces champs retournent Vrai.
981 //-- 1: Test de Fin(C1) = Debut(C2) ->N(C1) et P(C2)
982 //-- 2: Debut(C1) = Fin(C2) ->P(C1) et N(C2)
983 //=======================================================================
984 void IntAna_IntQuadQuad::InternalSetNextAndPrevious()
986 Standard_Boolean NotLastOpenC2, NotFirstOpenC2;
987 Standard_Integer c1,c2;
988 Standard_Real aEps, aEPSILON_DISTANCE;
989 Standard_Real DInfC1,DSupC1, DInfC2,DSupC2;
992 aEPSILON_DISTANCE=0.0000000001;
994 for(c1=0; c1<NbCurves; c1++) {
996 previouscurve[c1] = 0;
999 for(c1=0; c1 < NbCurves; c1++) {
1000 TheCurve[c1].Domain(DInfC1,DSupC1);
1002 for(c2 = 0; (c2 < NbCurves) && (c2!=c1) ; c2++) {
1004 NotLastOpenC2 = ! TheCurve[c2].IsLastOpen();
1005 NotFirstOpenC2 = ! TheCurve[c2].IsFirstOpen();
1006 TheCurve[c2].Domain(DInfC2,DSupC2);
1007 if(TheCurve[c1].IsFirstOpen() == Standard_False) {
1009 if(Abs(DInfC1-DSupC2)<=aEps &&
1010 (TheCurve[c1].Value(DInfC1)
1011 .Distance(TheCurve[c2].Value(DSupC2))<aEPSILON_DISTANCE)) {
1012 previouscurve[c1] = c2+1;
1013 nextcurve[c2] = c1+1;
1016 if(NotFirstOpenC2) {
1017 if(Abs(DInfC1-DInfC2)<=aEps
1018 && (TheCurve[c1].Value(DInfC1)
1019 .Distance(TheCurve[c2].Value(DInfC2))<aEPSILON_DISTANCE)) {
1020 previouscurve[c1] = -(c2+1);
1021 previouscurve[c2] = -(c1+1);
1025 if(TheCurve[c1].IsLastOpen() == Standard_False) {
1027 if(Abs(DSupC1-DSupC2)<=aEps
1028 && (TheCurve[c1].Value(DSupC1)
1029 .Distance(TheCurve[c2].Value(DSupC2))<aEPSILON_DISTANCE)) {
1031 nextcurve[c1] = -(c2+1);
1032 nextcurve[c2] = -(c1+1);
1035 if(NotFirstOpenC2) {
1036 if(Abs(DSupC1-DInfC2)<=aEps
1037 && (TheCurve[c1].Value(DSupC1)
1038 .Distance(TheCurve[c2].Value(DInfC2))<aEPSILON_DISTANCE)) {
1039 nextcurve[c1] = c2+1;
1040 previouscurve[c2] = c1+1;
1047 //=======================================================================
1048 //function :HasPreviousCurve
1050 //=======================================================================
1051 Standard_Boolean IntAna_IntQuadQuad::HasPreviousCurve(const Standard_Integer I) const
1054 throw StdFail_NotDone("IntQuadQuad Not done");
1057 throw Standard_DomainError("IntQuadQuad identical");
1059 if((I>NbCurves)||(I<=0)) {
1060 throw Standard_OutOfRange("Incorrect Curve Number 'HasPrevious Curve'");
1062 if(previouscurve[I-1]) {
1063 return Standard_True;
1065 return Standard_False;
1067 //=======================================================================
1068 //function :HasNextCurve
1070 //=======================================================================
1071 Standard_Boolean IntAna_IntQuadQuad::HasNextCurve(const Standard_Integer I) const
1074 throw StdFail_NotDone("IntQuadQuad Not done");
1077 throw Standard_DomainError("IntQuadQuad identical");
1079 if((I>NbCurves)||(I<=0)) {
1080 throw Standard_OutOfRange("Incorrect Curve Number 'HasNextCurve'");
1082 if(nextcurve[I-1]) {
1083 return Standard_True;
1085 return(Standard_False);
1087 //=======================================================================
1088 //function :PreviousCurve
1090 //=======================================================================
1091 Standard_Integer IntAna_IntQuadQuad::PreviousCurve (const Standard_Integer I,
1092 Standard_Boolean& Opposite) const
1094 if(HasPreviousCurve(I)) {
1095 if(previouscurve[I-1]>0) {
1096 Opposite = Standard_False;
1097 return(previouscurve[I-1]);
1100 Opposite = Standard_True;
1101 return( - previouscurve[I-1]);
1105 throw Standard_DomainError("Incorrect Curve Number 'PreviousCurve'");
1108 //=======================================================================
1109 //function :NextCurve
1111 //=======================================================================
1112 Standard_Integer IntAna_IntQuadQuad::NextCurve (const Standard_Integer I,
1113 Standard_Boolean& Opposite) const
1115 if(HasNextCurve(I)) {
1116 if(nextcurve[I]>0) {
1117 Opposite = Standard_False;
1118 return(nextcurve[I-1]);
1121 Opposite = Standard_True;
1122 return( - nextcurve[I-1]);
1126 throw Standard_DomainError("Incorrect Curve Number 'NextCurve'");
1129 //=======================================================================
1132 //=======================================================================
1133 const IntAna_Curve& IntAna_IntQuadQuad::Curve(const Standard_Integer i) const
1136 throw StdFail_NotDone("IntQuadQuad Not done");
1139 throw Standard_DomainError("IntQuadQuad identical");
1141 if((i <= 0) || (i>NbCurves)) {
1142 throw Standard_OutOfRange("Incorrect Curve Number");
1144 return TheCurve[i-1];
1146 //=======================================================================
1149 //=======================================================================
1150 const gp_Pnt& IntAna_IntQuadQuad::Point (const Standard_Integer i) const
1153 throw StdFail_NotDone("IntQuadQuad Not done");
1156 throw Standard_DomainError("IntQuadQuad identical");
1158 if((i <= 0) || (i>Nbpoints)) {
1159 throw Standard_OutOfRange("Incorrect Point Number");
1161 return Thepoints[i-1];
1163 //=======================================================================
1164 //function :Parameters
1166 //=======================================================================
1167 void IntAna_IntQuadQuad::Parameters (const Standard_Integer, //i,
1169 Standard_Real& ) const
1171 cout << "IntAna_IntQuadQuad::Parameters(...) is not yet implemented" << endl;
1174 /*********************************************************************************
1176 Mathematica Pour Cone Quadrique
1177 In[6]:= y[r_,t_]:=r*Sin[t]
1179 In[7]:= x[r_,t_]:=r*Cos[t]
1181 In[8]:= z[r_,t_]:=r*d
1183 In[14]:= Quad[x_,y_,z_]:=Qxx x x + Qyy y y + Qzz z z + 2 (Qxy x y + Qxz x z + Qyz y z + Qx x + Qy y + Qz z ) + Q1
1185 In[15]:= Quad[x[r,t],y[r,t],z[r,t]]
1188 Out[15]= Q1 + d Qzz r + Qxx r Cos[t] + Qyy r Sin[t] +
1191 > 2 (d Qz r + Qx r Cos[t] + d Qxz r Cos[t] + Qy r Sin[t] +
1194 > d Qyz r Sin[t] + Qxy r Cos[t] Sin[t])
1200 In[17]:= Collect[QQ,r]
1203 Out[17]= Q1 + r (2 d Qz + 2 Qx Cos[t] + 2 Qy Sin[t]) +
1206 > r (d Qzz + 2 d Qxz Cos[t] + Qxx Cos[t] + 2 d Qyz Sin[t] +
1209 > 2 Qxy Cos[t] Sin[t] + Qyy Sin[t] )
1210 ********************************************************************************/
1213 //**********************************************************************
1214 //*** C O N E - Q U A D R I Q U E ***
1216 //*** R ( d Qzz + 2 d Qxz Cos[t] + Qxx Cos[t] + 2 d Qyz Sin[t] + ***
1218 //*** 2 Qxy Cos[t] Sin[t] + Qyy Sin[t] ) ***
1220 //*** + R ( 2 d Qz + 2 Qx Cos[t] + 2 Qy Sin[t] ) ***
1223 //**********************************************************************
1224 //FortranForm= ( DIS = QQ1 QQ1 - 4 QQ0 QQ2 ) / 4
1225 // - d**2*Qz**2 - d**2*Qzz*Q1 + (Qx**2 - Qxx*Q1)*Cos(t)**2 +
1226 // - (2*d*Qy*Qz - 2*d*Qyz*Q1)*Sin(t) + (Qy**2 - Qyy*Q1)*Sin(t)**2 +
1227 // - Cos(t)*(2*d*Qx*Qz - 2*d*Qxz*Q1 + (2*Qx*Qy - 2*Qxy*Q1)*Sin(t))
1228 //**********************************************************************
1229 //modified by NIZNHY-PKV Fri Dec 2 10:56:03 2005f
1232 void DumpCurve(const Standard_Integer aIndex,
1234 //=======================================================================
1235 //function : DumpCurve
1237 //=======================================================================
1238 void DumpCurve(const Standard_Integer aIndex,
1241 Standard_Boolean bIsOpen, bIsConstant, bIsFirstOpen, bIsLastOpen;
1242 Standard_Integer i, aNb;
1243 Standard_Real aT1, aT2, aT, dT;
1246 aC.Domain(aT1, aT2);
1247 bIsOpen=aC.IsOpen();
1248 bIsConstant=aC.IsConstant();
1249 bIsFirstOpen=aC.IsFirstOpen();
1250 bIsLastOpen=aC.IsLastOpen();
1253 printf(" * IntAna_Curve #%d*\n", aIndex);
1254 printf(" Domain: [ %lf, %lf ]\n", aT1, aT2);
1255 printf(" IsOpen=%d\n", bIsOpen);
1256 printf(" IsConstant=%d\n", bIsConstant);
1257 printf(" IsFirstOpen =%d\n", bIsFirstOpen);
1258 printf(" IsLastOpen =%d\n", bIsLastOpen);
1261 dT=(aT2-aT1)/(aNb-1);
1262 for (i=0; i<aNb; ++i) {
1269 printf("point p%d_%d %lf %lf %lf\n",
1270 aIndex, i, aP.X(), aP.Y(), aP.Z());
1272 printf(" ---- end of curve ----\n");
1275 //modified by NIZNHY-PKV Fri Dec 2 10:42:16 2005t