0031035: Coding - uninitialized class fields reported by Visual Studio Code Analysis
[occt.git] / src / IntAna / IntAna_IntQuadQuad.cxx
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
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
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.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 #include <stdio.h>
18
19 #include <Standard_Stream.hxx>
20
21 #ifndef OCCT_DEBUG
22 #define No_Standard_RangeError
23 #define No_Standard_OutOfRange
24 #endif
25
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 //======================================================================
30
31 #include <ElSLib.hxx>
32 #include <gp_Ax2.hxx>
33 #include <gp_Ax3.hxx>
34 #include <gp_Cone.hxx>
35 #include <gp_Cylinder.hxx>
36 #include <gp_Pnt.hxx>
37 #include <IntAna_Curve.hxx>
38 #include <IntAna_IntQuadQuad.hxx>
39 #include <IntAna_Quadric.hxx>
40 #include <math_TrigonometricFunctionRoots.hxx>
41 #include <Standard_DomainError.hxx>
42 #include <Standard_OutOfRange.hxx>
43 #include <StdFail_NotDone.hxx>
44
45 //=======================================================================
46 //function : AddSpecialPoints
47 //purpose  : Sometimes the boundaries theTheta1 and theTheta2 are
48 //            computed with some inaccuracy. At that, some special points
49 //            (cone apex or sphere pole(s)), which are true intersection
50 //            points lie out of the domain [theTheta1, theTheta2] of the ALine.
51 //           This function corrects these boundaries to make them be included 
52 //            in the domain of the ALine.
53 //           Parameters Theta1 and Theta2 must be initialized
54 //            before calling this function.
55 //=======================================================================
56 template <class gpSmth>
57 static void AddSpecialPoints(const IntAna_Quadric& theQuad,
58                              const gpSmth& theGpObj,
59                              Standard_Real& theTheta1,
60                              Standard_Real& theTheta2)
61 {
62   const Standard_Real aPeriod = M_PI + M_PI;
63   const NCollection_List<gp_Pnt> &aLSP = theQuad.SpecialPoints();
64
65   if (aLSP.IsEmpty())
66     return;
67
68   Standard_Real aU = 0.0, aV = 0.0;
69   Standard_Real aMaxDelta = 0.0;
70   for (NCollection_List<gp_Pnt>::Iterator anItr(aLSP); anItr.More(); anItr.Next())
71   {
72     const gp_Pnt &aPt = anItr.Value();
73     ElSLib::Parameters(theGpObj, aPt, aU, aV);
74     const gp_Pnt aPProj(ElSLib::Value(aU, aV, theGpObj));
75
76     if (aPt.SquareDistance(aPProj) > Precision::SquareConfusion())
77     {
78       // aPt is not an intersection point
79       continue;
80     }
81
82     Standard_Real aDelta1 = Min(aU - theTheta1, 0.0),
83                   aDelta2 = Max(aU - theTheta2, 0.0);
84
85     if (aDelta1 < -M_PI)
86     {
87       // Must be aDelta1 = Min(aU - theTheta1 + aPeriod, 0.0).
88       // But aU - theTheta1 + aPeriod >= 0 always.
89       aDelta1 = 0.0;
90     }
91
92     if (aDelta2 > M_PI)
93     {
94       // Must be aDelta2 = Max(aU - theTheta2 - aPeriod, 0.0).
95       // But aU - theTheta2 - aPeriod <= 0 always.
96       aDelta2 = 0.0;
97     }
98
99     const Standard_Real aDelta = Max(-aDelta1, aDelta2);
100     aMaxDelta = Max(aMaxDelta, aDelta);
101   }
102
103   if(aMaxDelta != 0.0)
104   {
105     theTheta1 -= aMaxDelta;
106     theTheta2 += aMaxDelta;
107     if ((theTheta2 - theTheta1) > aPeriod)
108     {
109       theTheta2 = theTheta1 + aPeriod;
110     }
111   }
112 }
113
114 //=======================================================================
115 //class : TrigonometricRoots
116 //purpose: Classe Interne (Donne des racines classees d un polynome trigo)
117 //======================================================================
118 class TrigonometricRoots {
119
120  private:
121   Standard_Real Roots[4];
122   Standard_Boolean done;
123   Standard_Integer NbRoots;
124   Standard_Boolean infinite_roots;
125
126  public:
127   TrigonometricRoots(const Standard_Real CC,
128                      const Standard_Real SC,
129                      const Standard_Real C,
130                      const Standard_Real S,
131                      const Standard_Real Cte,
132                      const Standard_Real Binf,
133                      const Standard_Real Bsup);
134   //IsDone
135   Standard_Boolean IsDone() { 
136     return done; 
137   }
138   //IsARoot
139   Standard_Boolean IsARoot(Standard_Real u) {
140     Standard_Integer i;
141     Standard_Real aEps=RealEpsilon();
142     Standard_Real PIpPI = M_PI + M_PI;
143     //
144     for(i=0 ; i<NbRoots; ++i) {
145       if(Abs(u - Roots[i])<=aEps) {
146         return Standard_True;
147       }
148       if(Abs(u - Roots[i]-PIpPI)<=aEps) {
149         return Standard_True;
150       }
151     }
152     return Standard_False;
153   }
154   //NbSolutions
155   Standard_Integer NbSolutions() { 
156     if(!done) {
157       throw StdFail_NotDone();
158     }
159     return NbRoots; 
160   }
161   //InfiniteRoots
162   Standard_Boolean InfiniteRoots() { 
163     if(!done) {
164       throw StdFail_NotDone();
165     }    
166     return infinite_roots; 
167   }
168   //Value
169   Standard_Real Value(const Standard_Integer n) {
170     if((!done)||(n>NbRoots)) {
171       throw StdFail_NotDone();
172     }    
173     return Roots[n-1];
174   }
175 }; 
176 //=======================================================================
177 //function : TrigonometricRoots::TrigonometricRoots
178 //purpose  : 
179 //=======================================================================
180 TrigonometricRoots::TrigonometricRoots(const Standard_Real CC,
181                                        const Standard_Real SC,
182                                        const Standard_Real C,
183                                        const Standard_Real S,
184                                        const Standard_Real Cte,
185                                        const Standard_Real Binf,
186                                        const Standard_Real Bsup)
187 : infinite_roots(Standard_False)
188 {
189   Standard_Integer i, j, SvNbRoots;
190   Standard_Boolean Triee;
191   Standard_Real PIpPI = M_PI + M_PI;
192   //
193   done=Standard_False;
194   //
195   //-- F= AA*CN*CN+2*BB*CN*SN+CC*CN+DD*SN+EE;
196   math_TrigonometricFunctionRoots MTFR(CC, SC, C, S, Cte, Binf, Bsup); 
197   if(!MTFR.IsDone()) {
198     return;
199   }
200   //
201   done=Standard_True;
202   if(MTFR.InfiniteRoots()) {
203     infinite_roots=Standard_True;
204     return;
205   }
206   //
207   NbRoots=MTFR.NbSolutions();
208   //
209   for(i=0; i<NbRoots; ++i) {
210     Roots[i]=MTFR.Value(i+1);
211     if(Roots[i]<0.){
212       Roots[i]+=PIpPI;
213     }
214     if(Roots[i]>PIpPI) {
215       Roots[i]-=PIpPI;
216     }
217   }
218   //
219   //-- La recherche directe donne n importe quoi. 
220   SvNbRoots=NbRoots;
221   for(i=0; i<SvNbRoots; ++i) {
222     Standard_Real y;
223     Standard_Real co=cos(Roots[i]);
224     Standard_Real si=sin(Roots[i]);
225     y=co*(CC*co + (SC+SC)*si + C) + S*si + Cte;
226     if(Abs(y)>1e-8) {
227       done=Standard_False; 
228       return; //-- le 1er avril 98 
229     }
230   }
231   //
232   do {
233     Triee=Standard_True;
234     for(i=1,j=0; i<SvNbRoots; ++i,++j) {
235       if(Roots[i]<Roots[j]) {
236         Triee=Standard_False;
237         Standard_Real t=Roots[i]; 
238         Roots[i]=Roots[j]; 
239         Roots[j]=t;
240       }
241     }
242   }
243   while(!Triee);
244   //
245   infinite_roots=Standard_False;
246   //
247   if(!NbRoots) {//--!!!!! Detection du cas Pol = Cte ( 1e-50 ) !!!!
248     if((Abs(CC) + Abs(SC) + Abs(C) + Abs(S)) < 1e-10) {
249       if(Abs(Cte) < 1e-10)  {
250         infinite_roots=Standard_True;
251       }
252     }
253   }
254 }
255 //=======================================================================
256 //class    : MyTrigonometricFunction
257 //purpose  : 
258 //  Classe interne : Donne Value et Derivative sur un polynome 
259 //                   trigonometrique
260 //======================================================================
261 class MyTrigonometricFunction {
262  
263  private:
264   Standard_Real CC,SS,SC,S,C,Cte;
265  
266  public:
267   //
268   MyTrigonometricFunction(const Standard_Real xCC,
269                           const Standard_Real xSS,
270                           const Standard_Real xSC,
271                           const Standard_Real xC,
272                           const Standard_Real xS,
273                           const Standard_Real xCte) {
274     CC=xCC; 
275     SS=xSS; 
276     SC=xSC; 
277     S=xS; 
278     C=xC; 
279     Cte=xCte;
280   }
281  
282   Standard_Real Value(const Standard_Real& U) {
283     Standard_Real sinus, cosinus, aRet;
284     //
285     sinus=sin(U);
286     cosinus=cos(U);
287     aRet= CC*cosinus*cosinus + 
288           SS*sinus*sinus +
289           2.0*(sinus*(SC*cosinus+S)+cosinus*C)+
290           Cte;
291     //
292     return aRet;
293   }  
294   //
295   Standard_Real Derivative(const Standard_Real& U) {
296     Standard_Real sinus, cosinus;
297     //
298     sinus=sin(U);
299     cosinus=cos(U);
300     //
301     return(2.0*((sinus*cosinus)*(SS-CC)
302                 +S*cosinus
303                 -C*sinus
304                 +SC*(cosinus*cosinus-sinus*sinus)));
305   }
306 };
307
308 //////////
309 //=======================================================================
310 //function : IntAna_IntQuadQuad::IntAna_IntQuadQuad
311 //purpose  : C o n s t r u c t e u r    v i d e   
312 //=======================================================================
313 IntAna_IntQuadQuad::IntAna_IntQuadQuad(void) {
314   done=Standard_False;
315   identical = Standard_False;
316   NbCurves=0;
317   Nbpoints=0;
318   myNbMaxCurves=12;
319   myEpsilon=0.00000001;
320   myEpsilonCoeffPolyNull=0.00000001;
321   memset (nextcurve, 0, sizeof (nextcurve));
322   memset (previouscurve, 0, sizeof (previouscurve));
323 }
324 //=======================================================================
325 //function : IntAna_IntQuadQuad::IntAna_IntQuadQuad
326 //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 
327 //=======================================================================
328 IntAna_IntQuadQuad::IntAna_IntQuadQuad(const gp_Cylinder& Cyl,
329                                        const IntAna_Quadric& Quad,
330                                        const Standard_Real Tol) {
331   myNbMaxCurves=12;
332   myEpsilon=0.00000001;
333   myEpsilonCoeffPolyNull=0.00000001;
334   Perform(Cyl,Quad,Tol);
335 }
336 //=======================================================================
337 //function : Perform
338 //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 
339 //=======================================================================
340 void IntAna_IntQuadQuad::Perform(const gp_Cylinder& Cyl,
341                                  const IntAna_Quadric& Quad,
342                                  const Standard_Real) 
343 {
344   done = Standard_True;
345   identical= Standard_False;
346   NbCurves=0;
347   Nbpoints=0;
348   //
349   Standard_Boolean UN_SEUL_Z_PAR_THETA, DEUX_Z_PAR_THETA, 
350                    Z_POSITIF, Z_INDIFFERENT, Z_NEGATIF;
351   //
352   UN_SEUL_Z_PAR_THETA=Standard_False;
353   DEUX_Z_PAR_THETA=Standard_True;
354   Z_POSITIF=Standard_True;
355   Z_INDIFFERENT=Standard_True;
356   Z_NEGATIF=Standard_False;
357   //
358   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, aRealEpsilon, RCyl, R2;
359   Standard_Real PIpPI = M_PI + M_PI;
360   //
361   for(Standard_Integer raz = 0 ; raz < myNbMaxCurves ; raz++) {
362     previouscurve[raz] = nextcurve[raz] = 0;
363   }
364   //
365   RCyl=Cyl.Radius();
366   aRealEpsilon=RealEpsilon();
367   //----------------------------------------------------------------------
368   //-- Coefficients de la quadrique : 
369   //--      2       2       2
370   //-- Qxx x + Qyy y + Qzz z + 2 ( Qxy x y + Qxz x z + Qyz y z )
371   //-- + 2 ( Qx x + Qy y + Qz z ) + Q1
372   //----------------------------------------------------------------------
373   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1);
374   
375   //----------------------------------------------------------------------
376   //-- Ecriture des Coefficients de la Quadrique dans le repere liee 
377   //-- au Cylindre 
378   //--                 Cyl.Position() -> Ax2
379   //----------------------------------------------------------------------
380   Quad.NewCoefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,Cyl.Position());
381   
382   //----------------------------------------------------------------------
383   //-- Parametrage du Cylindre Cyl : 
384   //--     X = Rcyl * Cos(Theta)
385   //--     Y = Rcyl * Sin(Theta) 
386   //--     Z = Z
387   //----------------------------------------------------------------------
388   //-- Donne une Equation de la forme :
389   //--   F(Sin(Theta),Cos(Theta),ZCyl) = 0
390   //--   (Equation du second degre en ZCyl) 
391   //--    ZCyl**2  CoeffZ2(Theta) + ZCyl CoeffZ1(Theta) + CoeffZ0(Theta)
392   //----------------------------------------------------------------------
393   //-- CoeffZ0 = Q1 + 2*Qx*RCyl*Cos[Theta] + Qxx*RCyl^2*Cos[Theta]^2
394   //--           2*RCyl*Sin[Theta]* ( Qy + Qxy*RCyl*Cos[Theta]);
395   //--           Qyy*RCyl^2*Sin[Theta]^2;
396   //-- CoeffZ1 =2.0 * ( Qz + RCyl*(Qxz*Cos[Theta] + Sin[Theta]*Qyz)) ;
397   //-- CoeffZ2 =  Qzz;
398   //-- !!!! Attention , si on norme sur Qzz pour detecter le cas 1er degre 
399   //----------------------------------------------------------------------
400   //-- On Cherche Les Racines en Theta du discriminant de cette equation :
401   //-- DIS(Theta) = C_1 + C_SS Sin()**2 + C_CC Cos()**2 + 2 C_SC Sin() Cos()
402   //--                  + 2 C_S  Sin()    + 2 C_C Cos()
403   //--
404   //-- Si Qzz = 0   Alors  On Resout Z=Fct(Theta)  sur le domaine de Theta
405   //----------------------------------------------------------------------
406
407   if(Abs(Qzz)<myEpsilonCoeffPolyNull) {
408     done=Standard_False; //-- le 12 mars 98    
409     return;
410   }  
411   else { //#0
412     //----------------------------------------------------------------------
413     //--- Cas  ou  F(Z,Theta) est du second degre en Z                  ----
414     //----------------------------------------------------------------------
415     R2   = RCyl*RCyl;
416     
417     Standard_Real C_1  = Qz*Qz - Qzz*Q1;
418     Standard_Real C_SS = R2    * ( Qyz*Qyz - Qyy*Qzz );
419     Standard_Real C_CC = R2    * ( Qxz*Qxz - Qxx*Qzz );
420     Standard_Real C_S  = RCyl  * ( Qyz*Qz  - Qy*Qzz  );
421     Standard_Real C_C  = RCyl  * ( Qxz*Qz  - Qx*Qzz  );
422     Standard_Real C_SC = R2    * ( Qxz*Qyz - Qxy*Qzz );
423     //
424     MyTrigonometricFunction MTF(C_CC,C_SS,C_SC,C_C,C_S,C_1);
425     TrigonometricRoots PolDIS(C_CC-C_SS,
426                               C_SC,
427                               C_C+C_C,
428                               C_S+C_S,
429                               C_1+C_SS, 0., PIpPI);
430     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
431     if(PolDIS.IsDone()==Standard_False) {
432       done=Standard_False;
433       return;
434     }
435     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
436     if(PolDIS.InfiniteRoots()) {   
437       TheCurve[0].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
438                                         myEpsilon,0.,PIpPI,
439                                         UN_SEUL_Z_PAR_THETA,
440                                         Z_POSITIF);
441       TheCurve[1].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
442                                         myEpsilon,0.,PIpPI,
443                                         UN_SEUL_Z_PAR_THETA,
444                                         Z_NEGATIF);
445       NbCurves=2;
446     }
447
448     else { //#1
449       //---------------------------------------------------------------
450       //-- La Recherche des Zero de DIS a reussi
451       //---------------------------------------------------------------
452       Standard_Integer nbsolDIS=PolDIS.NbSolutions();
453       if(nbsolDIS==0) {
454         //--------------------------------------------------------------
455         //-- Le Discriminant a un signe constant : 
456         //-- 
457         //-- Si Positif  ---> 2 Courbes
458         //-- Sinon       ---> Pas de solution
459         //--------------------------------------------------------------
460         if(MTF.Value(M_PI) >= -aRealEpsilon) {
461
462           TheCurve[0].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
463                                             myEpsilon,0.0,PIpPI,
464                                             UN_SEUL_Z_PAR_THETA,
465                                             Z_POSITIF);
466           TheCurve[1].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
467                                             myEpsilon,0.0,PIpPI,
468                                             UN_SEUL_Z_PAR_THETA,
469                                             Z_NEGATIF);
470           
471           NbCurves=2;
472         }
473         else {
474           //------------------------------------------------------------
475           //-- Le Discriminant est toujours Negatif
476           //------------------------------------------------------------
477           NbCurves=0;
478         }
479       }
480       else { //#2
481         //------------------------------------------------------------
482         //-- Le Discriminant a des racines 
483         //-- (Le Discriminant est une fonction periodique donc ... )
484         //------------------------------------------------------------
485         if( nbsolDIS==1 ) {
486           //------------------------------------------------------------
487           //-- Point de Tangence pour ce Theta Solution
488           //-- Si Signe de Discriminant >=0 pour tout Theta
489           //--     Alors  
490           //--       Courbe Solution pour la partie Positive
491           //--       entre les 2 racines ( Ici Tout le domaine )
492           //-- Sinon Seulement un point Tangent
493           //------------------------------------------------------------
494           if(MTF.Value(PolDIS.Value(1)+M_PI) >= -aRealEpsilon ) {
495             //------------------------------------------------------------
496             //-- On a Un Point de Tangence + Une Courbe Solution
497             //------------------------------------------------------------
498             TheCurve[0].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
499                                               myEpsilon,0.0,PIpPI,
500                                               UN_SEUL_Z_PAR_THETA,
501                                               Z_POSITIF);
502             TheCurve[1].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
503                                               myEpsilon,0.0,PIpPI,
504                                               UN_SEUL_Z_PAR_THETA,
505                                               Z_NEGATIF);
506             
507             NbCurves=2;
508           }
509           else {
510             //------------------------------------------------------------
511             //-- On a simplement un Point de tangence
512             //------------------------------------------------------------
513             //--Standard_Real Theta = PolDIS(1);
514             //--Standard_Real SecPar= -0.5 * MTFZ1.Value(Theta) / MTFZ2.Value(Theta);
515             //--Thepoints[Nbpoints] = ElSLib::CylinderValue(Theta,SecPar,Cyl.Position(),Cyl.Radius());
516             //--Nbpoints++;
517             NbCurves=0;
518           }
519         }
520         else {  // #3
521           //------------------------------------------------------------
522           //-- On detecte   :  Les racines double         
523           //--              :  Les Zones Positives [Modulo 2 PI]
524           //-- Les Courbes solutions seront obtenues pour les valeurs
525           //-- de Theta ou Discriminant(Theta) > 0  (>=0? en limite)
526           //-- Par la resolution de l equation implicite avec Theta fixe
527           //------------------------------------------------------------
528           //-- Si tout est solution, Alors on est sur une iso 
529           //-- Ce cas devrait etre traite en amont
530           //------------------------------------------------------------
531           //-- On construit la fonction permettant connaissant un Theta
532           //-- de calculer les Z+ et Z- Correspondants.
533           //-------------------------------------------------------------
534
535           //-------------------------------------------------------------
536           //-- Calcul des Intervalles ou Discriminant >=0
537           //-- On a au plus nbsol intervalles ( en fait 2 )
538           //--  (1,2) (2,3) .. (nbsol,1+2PI)
539           //-------------------------------------------------------------
540           Standard_Integer i;
541           Standard_Real Theta1, Theta2, Theta3, autrepar, qwet;
542           Standard_Boolean UnPtTg = Standard_False;
543           //
544           NbCurves=0;
545           if(nbsolDIS == 2) { 
546             for(i=1; i<=nbsolDIS ; i++) {
547               Theta1=PolDIS.Value(i);
548               Theta2=(i<nbsolDIS)? PolDIS.Value(i+1)  :  (PolDIS.Value(1)+PIpPI);
549               //----------------------------------------------------------------
550               //-- On detecte les racines doubles 
551               //-- Si il n y a que 2 racines alors on prend tout l interval
552               //----------------------------------------------------------------
553               if(Abs(Theta2-Theta1)<=aRealEpsilon) {
554                 UnPtTg = Standard_True;
555                 autrepar=Theta1-0.1; 
556                 if(autrepar<0.) {
557                   autrepar=Theta1+0.1;
558                 }
559                 //
560                 qwet=MTF.Value(autrepar);
561                 if(qwet>=0.) { 
562                   Standard_Real aParam = Theta1 + PIpPI;
563                   AddSpecialPoints(Quad, Cyl, Theta1, aParam);
564                   TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
565                                                            myEpsilon,Theta1,aParam,
566                                                            UN_SEUL_Z_PAR_THETA,
567                                                            Z_POSITIF);
568                   NbCurves++;
569                   TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
570                                                            myEpsilon,Theta1,aParam,
571                                                            UN_SEUL_Z_PAR_THETA,
572                                                            Z_NEGATIF);
573                   NbCurves++;
574                 }
575               }
576             }
577           }
578           
579           for(i=1; UnPtTg == (Standard_False) && (i<=nbsolDIS) ; i++) {
580             Theta1=PolDIS.Value(i);
581             Theta2=(i<nbsolDIS)? PolDIS.Value(i+1) : (PolDIS.Value(1)+PIpPI);
582             //----------------------------------------------------------------
583             //-- On detecte les racines doubles 
584             //-- Si il n y a que 2 racines alors on prend tout l interval
585             //----------------------------------------------------------------
586             if(Abs(Theta2-Theta1)<=1e-12) {
587               //-- std::cout<<"\n####### Un Point de Tangence en Theta = "<<Theta1<<std::endl;
588               //-- i++;
589             }
590             else {  //-- On evite les pbs numeriques (Tout est du meme signe entre les racines) 
591               qwet=MTF.Value(0.5*(Theta1+Theta2))
592                     +MTF.Value(0.4*Theta1+0.6*Theta2)
593                       +MTF.Value(0.6*Theta1+0.4*Theta2);
594               if(qwet >= 0.) {
595                 //------------------------------------------------------------
596                 //-- On est positif a partir de Theta1
597                 //-- L intervalle Theta1,Theta2 est retenu
598                 //------------------------------------------------------------
599                 
600                 //-- le 8 octobre 1997 : 
601                 //-- PB: Racine en 0    pi/2 pi/2  et 2pi
602                 //-- On cree 2 courbes : 0    -> pi/2    2zpartheta
603                 //--                     pi/2 -> 2pi     2zpartheta
604                 //-- 
605                 //-- la courbe 0 pi/2   passe par le double pt de tangence pi/2
606                 //-- il faut donc couper cette courbe en 2
607                 //--
608                 Theta3 = ((i+1)<nbsolDIS)? PolDIS.Value(i+2)  :  (PolDIS.Value(1)+PIpPI);
609                 //ft
610                 if((Theta3-Theta2)<5.e-8) { 
611                 //
612                   AddSpecialPoints(Quad, Cyl, Theta1, Theta2);
613                   TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
614                                                            myEpsilon,Theta1,Theta2,
615                                                            UN_SEUL_Z_PAR_THETA,
616                                                            Z_POSITIF);
617                   NbCurves++;
618                   TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
619                                                            myEpsilon,Theta1,Theta2,
620                                                            UN_SEUL_Z_PAR_THETA,
621                                                            Z_NEGATIF);
622                   NbCurves++;
623                 }
624                 else { 
625                   AddSpecialPoints(Quad, Cyl, Theta1, Theta2);
626                   TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
627                                                            myEpsilon,Theta1,Theta2,
628                                                            DEUX_Z_PAR_THETA,
629                                                            Z_INDIFFERENT);
630                   NbCurves++;
631                 }
632               }//if(qwet >= 0.)
633             }//else {  //-- On evite les pbs numerique ...      
634           } //for(i=1; UnPtTg == (Standard_False) && (i<=nbsolDIS) ; i++) {    
635         }//else { // #3
636       }//else { //#2
637     }//else { //#1
638     
639   }//else { //#0
640 }
641
642 //=======================================================================
643 //function :IntAna_IntQuadQuad::IntAna_IntQuadQuad
644 //purpose  : 
645 //=======================================================================
646 IntAna_IntQuadQuad::IntAna_IntQuadQuad(const gp_Cone& Cone,
647                                        const IntAna_Quadric& Quad,
648                                        const Standard_Real Tol) 
649
650   myNbMaxCurves=12;
651   myEpsilon=0.00000001;
652   myEpsilonCoeffPolyNull=0.00000001;
653   Perform(Cone,Quad,Tol);
654 }
655 //=======================================================================
656 //function :Perform
657 //purpose  : 
658 //=======================================================================
659 void IntAna_IntQuadQuad::Perform(const gp_Cone& Cone,
660                                  const IntAna_Quadric& Quad,
661                                  const Standard_Real)  
662
663   //
664   Standard_Boolean UN_SEUL_Z_PAR_THETA,
665                    Z_POSITIF, Z_NEGATIF;
666   //
667   UN_SEUL_Z_PAR_THETA=Standard_False;
668   Z_POSITIF=Standard_True;
669   Z_NEGATIF=Standard_False;
670   //
671   Standard_Integer i;
672   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1;
673   Standard_Real Theta1, Theta2, TgAngle;
674   Standard_Real PIpPI = M_PI + M_PI;
675   //
676   done=Standard_True;
677   identical = Standard_False;
678   NbCurves=0;
679   Nbpoints=0;
680   //
681   for(i=0 ; i<myNbMaxCurves ; ++i) {
682     previouscurve[i]=0;
683     nextcurve[i]=0;
684   }
685   //
686   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1);
687   //
688   gp_Ax3 tAx3(Cone.Position());
689   tAx3.SetLocation(Cone.Apex());
690   Quad.NewCoefficients(Qxx,Qyy,Qzz,
691                        Qxy,Qxz,Qyz,
692                        Qx,Qy,Qz,Q1,
693                        tAx3);
694   //
695   TgAngle=1./(Tan(Cone.SemiAngle()));
696   //
697   // The parametrization of the Cone
698   // 
699   // x= z*tan(beta)*cos(t)
700   // y= z*tan(beta)*sin(t)
701   // z=z
702   // 
703   // The intersection curves are defined by the equation
704   //
705   //                        2
706   //          f(z,t)= A(t)*z + B(t)*z + C(t)=0
707   //
708   //   
709   // 1. Try to solve A(t)=0 -> PolZ2
710   //
711   Standard_Integer nbsol, nbsol1, nbsolZ2;
712   Standard_Real Z2CC, Z2SS, Z2Cte, Z2SC, Z2C, Z2S;
713   Standard_Real Z1CC, Z1SS, Z1Cte, Z1SC, Z1C, Z1S; 
714   Standard_Real C_1,C_SS, C_CC, C_S, C_C, C_SC;
715   //
716   Z2CC = Qxx;
717   Z2SS = Qyy;
718   Z2Cte= Qzz * TgAngle*TgAngle;
719   Z2SC = Qxy;
720   Z2C  = (TgAngle)*Qxz;
721   Z2S  = (TgAngle)*Qyz;
722   //
723   TrigonometricRoots PolZ2(Z2CC-Z2SS,Z2SC,Z2C+Z2C,Z2S+Z2S,Z2Cte+Z2SS,0.,PIpPI);
724   if(!PolZ2.IsDone()) {
725     done=!done;
726     return;
727   } 
728   //
729   //MyTrigonometricFunction MTF2(Z2CC,Z2SS,Z2SC,Z2C,Z2S,Z2Cte);
730   nbsolZ2 = PolZ2.NbSolutions();
731   //
732   // 2. Try to solve B(t)=0  -> PolZ1
733   //
734   Z1Cte = 2.*(TgAngle)*Qz;
735   Z1SS  = 0.;
736   Z1CC  = 0.;
737   Z1S   = Qy;
738   Z1C   = Qx;
739   Z1SC  = 0.;
740   //
741   TrigonometricRoots PolZ1(Z1CC-Z1SS,Z1SC, Z1C+Z1C,Z1S+Z1S, Z1Cte+Z1SS,0.,PIpPI);
742   if(!PolZ1.IsDone()) {
743     done=!done;
744     return;
745   }
746   MyTrigonometricFunction MTFZ1(Z1CC,Z1SS,Z1SC,Z1C,Z1S,Z1Cte);
747   //
748   nbsol1=PolZ1.NbSolutions();
749   if(PolZ2.InfiniteRoots()) { // i.e A(t)=0 for any t
750     if(!PolZ1.InfiniteRoots()) {    
751       if(!nbsol1) {
752         //-- B(t)*z + C(t) = 0    avec C(t) != 0
753         TheCurve[0].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
754                                       myEpsilon,0.,PIpPI,
755                                       UN_SEUL_Z_PAR_THETA,
756                                       Z_POSITIF);
757         NbCurves=1;
758       }
759       else {
760         /*
761         Standard_Integer ii;
762         for(ii=1; ii<= nbsol1 ; ++ii) {
763           Standard_Real Theta=PolZ1.Value(ii);
764           if(Abs(MTFZ1.Value(Theta))<=myEpsilon) {
765             //-- Une droite Solution  Z=  -INF -> +INF  pour Theta
766             //-- std::cout<<"######## Droite Solution Pour Theta = "<<Theta<<std::endl;
767           }
768           else {
769             //-- std::cout<<"\n#### _+_+_+_+_+ CAS A(t) Z + B = 0   avec A et B ==0 "<<std::endl;
770           }
771         }
772         */
773       }
774     }
775     else {
776       if(Abs(Q1)<=myEpsilon) {
777         done=!done;
778         return;
779       }
780       else {
781         //-- Pas de Solutions 
782       }
783     }
784     return;     
785   }
786   //
787   //else { //#5 
788   //
789   //                2
790   //-- f(z,t)=A(t)*z + B(t)*z + C(t)=0   avec A n est pas toujours nul
791   //
792   //                                              2
793   // 3. Try to explore s.c. Discriminant:  D(t)=B(t)-4*A(t)*C(t) => Pol
794   //
795   // The function D(t) is just a formula that has sense for quadratic
796   // equation above.
797   // For cases when A(t)=0 (say at t=ti, t=tj. etc) the equation
798   // will be
799   //
800   //     f(z,t)=B(t)*z + C(t)=0, where B(t)!=0, 
801   //
802   // and the D(t) is nonsense for it.
803   //
804   C_1  = TgAngle*TgAngle * ( Qz * Qz - Qzz * Q1);
805   C_SS = Qy * Qy - Qyy * Q1;
806   C_CC = Qx * Qx - Qxx * Q1; 
807   C_S  = TgAngle*( Qy * Qz - Qyz * Q1);
808   C_C  = TgAngle*( Qx * Qz - Qxz * Q1);
809   C_SC = Qx * Qy - Qxy * Q1;
810   //
811   TrigonometricRoots Pol(C_CC-C_SS, C_SC, C_C+C_C,C_S+C_S, C_1+C_SS,0.,PIpPI);
812   if(!Pol.IsDone()) {
813     done=!done;
814     return;
815   }
816   //
817   nbsol=Pol.NbSolutions();
818   MyTrigonometricFunction MTF(C_CC,C_SS,C_SC,C_C,C_S,C_1);
819   //
820   if(Pol.InfiniteRoots()) { 
821     //                             2
822     //         f(z,t)= (z(t)-zo(t)) = 0       Discriminant(t)=0 pour tout t 
823     //
824     TheCurve[0].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
825                                   0.,PIpPI,
826                                   UN_SEUL_Z_PAR_THETA, Z_POSITIF);
827     TheCurve[1].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
828                                   0.,PIpPI,
829                                   UN_SEUL_Z_PAR_THETA, Z_NEGATIF);
830     NbCurves=2;
831     return;
832   }
833   //else {//#4
834   //                     2
835   //        f(z,t)=A(t)*z + B(t)*z + C(t)      Discriminant(t) != 0 
836   //
837   if(!nbsol && (MTF.Value(M_PI)<0.) ) {
838     //-- Discriminant signe constant negatif
839     return;
840   }
841   //else {//# 3
842   //
843   //-- On Traite le cas : Discriminant(t) > 0 pour tout t en simulant 1 
844   //   racine double en 0
845   Standard_Boolean DiscriminantConstantPositif = Standard_False;
846   if(!nbsol) { 
847     nbsol=1;
848     DiscriminantConstantPositif = Standard_True;
849   }
850   //----------------------------------------------------------------------
851   //-- Le discriminant admet au moins une racine ( -> Point de Tangence )
852   //-- ou est constant positif.
853   //----------------------------------------------------------------------
854   for(i=1; i<=nbsol; ++i) {
855     if(DiscriminantConstantPositif) {
856       Theta1 = 0.;
857       Theta2 = PIpPI-myEpsilon; 
858     }
859     else {
860       Theta1=Pol.Value(i);
861       Theta2=(i<nbsol)? Pol.Value(i+1) : (Pol.Value(1)+PIpPI);
862     }
863     //
864     if(Abs(Theta2-Theta1)<=myEpsilon){
865       done=Standard_False; 
866       return;// !!! pkv
867     }
868     //else {// #2
869     Standard_Real qwet =MTF.Value(0.5*(Theta1+Theta2))+
870                         MTF.Value(0.4*Theta1+0.6*Theta2)+
871                         MTF.Value(0.6*Theta1+0.4*Theta2);
872     if(qwet < 0.) {
873       continue;
874     }
875     //---------------------------------------------------------------------
876     //-- On a des Solutions entre Theta1 et Theta 2 
877     //---------------------------------------------------------------------
878     
879     //---------------------------------------------------------------------
880     //-- On Subdivise si necessaire Theta1-->Theta2
881     //-- en Theta1-->t1   t1--->t2   ...  tn--->Theta2
882     //--  ou t1,t2,...tn sont des racines de A(t)
883     //--
884     //-- Seule la courbe Z_NEGATIF est affectee
885     //----------------------------------------------------------------------
886     Standard_Boolean RacinesdePolZ2DansTheta1Theta2;
887     Standard_Integer i2;
888     Standard_Real r;
889     //
890     //nbsolZ2=PolZ2.NbSolutions();
891     RacinesdePolZ2DansTheta1Theta2=Standard_False;
892     for(i2=1; i2<=nbsolZ2 && !RacinesdePolZ2DansTheta1Theta2; ++i2) {
893       r=PolZ2.Value(i2);
894       if(r>Theta1 && r<Theta2) {
895         RacinesdePolZ2DansTheta1Theta2=Standard_True;
896       }
897       else { 
898         r+=PIpPI;
899         if(r>Theta1 && r<Theta2){ 
900           RacinesdePolZ2DansTheta1Theta2=Standard_True;
901         }
902       }
903     }
904     //
905     if(!RacinesdePolZ2DansTheta1Theta2) {
906       //--------------------------------------------------------------------
907       //-- Pas de Branches Infinies 
908       TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,myEpsilon,
909                                            Theta1,Theta2,
910                                            UN_SEUL_Z_PAR_THETA,Z_POSITIF);
911       NbCurves++;
912       TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,myEpsilon,
913                                            Theta1,Theta2,
914                                            UN_SEUL_Z_PAR_THETA,
915                                            Z_NEGATIF);
916       NbCurves++;
917     }
918     
919     else { //#1
920       Standard_Boolean NoChanges;
921       Standard_Real NewMin, NewMax, to;
922       //
923       NewMin=Theta1;
924       NewMax=Theta2;
925       NoChanges=Standard_True;
926       //
927       for(i2=1; i2<= (nbsolZ2+nbsolZ2) ; ++i2) {
928         if(i2>nbsolZ2) {
929           to=PolZ2.Value(i2-nbsolZ2) + PIpPI; 
930         }
931         else {
932           to=PolZ2.Value(i2);
933         }
934         //
935         // to is value of the parameter t where A(t)=0, i.e. 
936         // f(z,to)=B(to)*z + C(to)=0, B(to)!=0. 
937         //                   
938         // z=-C(to)/B(to) is one and only
939         // solution inspite of the fact that D(t)>0 for any value of t.
940         //
941         if(to<NewMax && to>NewMin) {
942           //-----------------------------------------------------------------
943           //-- On coupe au moins une fois le domaine Theta1 Theta2
944           //-----------------------------------------------------------------
945           NoChanges=Standard_False;
946           TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
947                                                NewMin,to,
948                                                UN_SEUL_Z_PAR_THETA, Z_NEGATIF);
949           //
950           NbCurves++;
951           TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
952                                                NewMin,to,
953                                                UN_SEUL_Z_PAR_THETA, Z_POSITIF);
954           //------------------------------------------------------------
955           //-- A == 0    
956           //-- Si B < 0    Alors Infini sur   Z+  
957           //-- Sinon             Infini sur   Z-
958           //------------------------------------------------------------
959           if(PolZ2.IsARoot(NewMin)) {
960             if(MTFZ1.Value(NewMin) < 0.) {
961               TheCurve[NbCurves].SetIsFirstOpen(Standard_True);
962             }
963             else {
964               TheCurve[NbCurves-1].SetIsFirstOpen(Standard_True);
965             }
966           }
967           if(MTFZ1.Value(to) < 0.) {
968             TheCurve[NbCurves].SetIsLastOpen(Standard_True);
969           }
970           else {
971             TheCurve[NbCurves-1].SetIsLastOpen(Standard_True);
972           }
973           //------------------------------------------------------------
974           NbCurves++;
975           NewMin=to;
976         }//if(to<NewMax && to>NewMin)
977       }// for(i2=1; i2<= (nbsolZ2+nbsolZ2) ; ++i2)
978       //
979       if(NoChanges) {
980         TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,  myEpsilon,
981                                              Theta1,Theta2,
982                                              UN_SEUL_Z_PAR_THETA, Z_NEGATIF);
983         NbCurves++;
984         TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
985                                              Theta1,Theta2,
986                                              UN_SEUL_Z_PAR_THETA, Z_POSITIF);
987         //------------------------------------------------------------
988         //-- A == 0    
989         //-- Si B < 0    Alors Infini sur   Z+  
990         //-- Sinon             Infini sur   Z-
991         //------------------------------------------------------------
992         if(PolZ2.IsARoot(Theta1)) {
993           if(MTFZ1.Value(Theta1) < 0.) {
994             TheCurve[NbCurves].SetIsFirstOpen(Standard_True);
995           }
996           else {
997             TheCurve[NbCurves-1].SetIsFirstOpen(Standard_True);
998           }
999         }
1000         if(PolZ2.IsARoot(Theta2)) {
1001           if(MTFZ1.Value(Theta2) < 0.) {
1002             TheCurve[NbCurves].SetIsLastOpen(Standard_True);
1003           }
1004           else {
1005             TheCurve[NbCurves-1].SetIsLastOpen(Standard_True);
1006           }
1007         }
1008         //------------------------------------------------------------
1009         NbCurves++;
1010       }//if(NoChanges) {
1011       else {// #0
1012         TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
1013                                              NewMin,Theta2,
1014                                              UN_SEUL_Z_PAR_THETA, Z_NEGATIF);
1015         NbCurves++;
1016         TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
1017                                              NewMin,Theta2,
1018                                              UN_SEUL_Z_PAR_THETA, Z_POSITIF);
1019         //------------------------------------------------------------
1020         //-- A == 0    
1021         //-- Si B < 0    Alors Infini sur   Z+  
1022         //-- Sinon             Infini sur   Z-
1023         //------------------------------------------------------------
1024         if(PolZ2.IsARoot(NewMin)) {
1025           if(MTFZ1.Value(NewMin) < 0.) {
1026             TheCurve[NbCurves].SetIsFirstOpen(Standard_True);
1027           }
1028           else {
1029             TheCurve[NbCurves-1].SetIsFirstOpen(Standard_True);
1030           }
1031         }
1032         if(PolZ2.IsARoot(Theta2)) {
1033           if(MTFZ1.Value(Theta2) < 0.) {
1034             TheCurve[NbCurves].SetIsLastOpen(Standard_True);
1035           }
1036           else {
1037             TheCurve[NbCurves-1].SetIsLastOpen(Standard_True);
1038           }
1039         }
1040         //------------------------------------------------------------
1041         
1042         NbCurves++;
1043       }//else {// #0
1044     }//else { //#1
1045   }//for(i=1; i<=nbsol ; i++) {
1046   //}//else { //#5
1047   InternalSetNextAndPrevious();
1048 }
1049 //=======================================================================
1050 //function :InternalSetNextAndPrevious
1051 //purpose  : 
1052 //-- Raccordement des courbes bout a bout 
1053 //--    - Utilisation du champ : IsFirstOpen 
1054 //--    -                        IsLastOpen
1055 //-- Pas de verification si ces champs retournent Vrai.
1056 //--
1057 //--
1058 //-- 1: Test de      Fin(C1)    = Debut(C2)   ->N(C1) et P(C2)
1059 //-- 2:              Debut(C1)  = Fin(C2)     ->P(C1) et N(C2)
1060 //=======================================================================
1061 void IntAna_IntQuadQuad::InternalSetNextAndPrevious() 
1062 {
1063   Standard_Boolean NotLastOpenC2, NotFirstOpenC2;
1064   Standard_Integer c1,c2;
1065   Standard_Real aEps, aEPSILON_DISTANCE;
1066   Standard_Real DInfC1,DSupC1, DInfC2,DSupC2;
1067   //
1068   aEps=0.0000001;
1069   aEPSILON_DISTANCE=0.0000000001;
1070   //
1071   for(c1=0; c1<NbCurves; c1++) {
1072     nextcurve[c1] =0;
1073     previouscurve[c1] = 0;
1074   }
1075   //
1076   for(c1=0; c1 < NbCurves; c1++) {
1077     TheCurve[c1].Domain(DInfC1,DSupC1);
1078     
1079     for(c2 = 0; (c2 < NbCurves) && (c2!=c1) ; c2++) {
1080       
1081       NotLastOpenC2  = ! TheCurve[c2].IsLastOpen();
1082       NotFirstOpenC2 = ! TheCurve[c2].IsFirstOpen();
1083       TheCurve[c2].Domain(DInfC2,DSupC2);
1084       if(TheCurve[c1].IsFirstOpen() == Standard_False) {
1085         if(NotLastOpenC2) {
1086           if(Abs(DInfC1-DSupC2)<=aEps && 
1087              (TheCurve[c1].Value(DInfC1)
1088               .Distance(TheCurve[c2].Value(DSupC2))<aEPSILON_DISTANCE)) {
1089             previouscurve[c1] = c2+1;
1090             nextcurve[c2]     = c1+1;
1091           }
1092         }
1093         if(NotFirstOpenC2) {
1094           if(Abs(DInfC1-DInfC2)<=aEps
1095              && (TheCurve[c1].Value(DInfC1)
1096                  .Distance(TheCurve[c2].Value(DInfC2))<aEPSILON_DISTANCE)) {
1097             previouscurve[c1] = -(c2+1);
1098             previouscurve[c2] = -(c1+1);
1099           }
1100         }
1101       }
1102       if(TheCurve[c1].IsLastOpen() == Standard_False) {
1103         if(NotLastOpenC2) {
1104           if(Abs(DSupC1-DSupC2)<=aEps
1105              && (TheCurve[c1].Value(DSupC1)
1106                  .Distance(TheCurve[c2].Value(DSupC2))<aEPSILON_DISTANCE)) {
1107             
1108             nextcurve[c1] = -(c2+1);
1109             nextcurve[c2] = -(c1+1);
1110           }
1111         }
1112         if(NotFirstOpenC2) {
1113           if(Abs(DSupC1-DInfC2)<=aEps 
1114              && (TheCurve[c1].Value(DSupC1) 
1115                  .Distance(TheCurve[c2].Value(DInfC2))<aEPSILON_DISTANCE)) { 
1116             nextcurve[c1]     = c2+1;
1117             previouscurve[c2] = c1+1;
1118           }
1119         }
1120       }
1121     }
1122   }
1123 }
1124 //=======================================================================
1125 //function :HasPreviousCurve
1126 //purpose  : 
1127 //=======================================================================    
1128 Standard_Boolean IntAna_IntQuadQuad::HasPreviousCurve(const Standard_Integer I) const
1129 {
1130   if(!done) {
1131     throw StdFail_NotDone("IntQuadQuad Not done");
1132   }  
1133   if (identical) {
1134     throw Standard_DomainError("IntQuadQuad identical");
1135   }
1136   if((I>NbCurves)||(I<=0)) {
1137     throw Standard_OutOfRange("Incorrect Curve Number 'HasPrevious Curve'");
1138   }
1139   if(previouscurve[I-1]) {
1140     return Standard_True;
1141   }
1142   return Standard_False;
1143 }
1144 //=======================================================================
1145 //function :HasNextCurve
1146 //purpose  : 
1147 //=======================================================================    
1148 Standard_Boolean IntAna_IntQuadQuad::HasNextCurve(const Standard_Integer I) const
1149 {
1150   if(!done) {
1151     throw StdFail_NotDone("IntQuadQuad Not done");
1152   }  
1153   if (identical) {
1154     throw Standard_DomainError("IntQuadQuad identical");
1155   }
1156   if((I>NbCurves)||(I<=0)) {
1157     throw Standard_OutOfRange("Incorrect Curve Number 'HasNextCurve'");
1158   }
1159   if(nextcurve[I-1]) {
1160     return Standard_True;
1161   }
1162   return(Standard_False);
1163 }
1164 //=======================================================================
1165 //function :PreviousCurve
1166 //purpose  : 
1167 //=======================================================================     
1168 Standard_Integer IntAna_IntQuadQuad::PreviousCurve  (const Standard_Integer I,
1169                                                      Standard_Boolean& theOpposite) const
1170 {
1171   if(HasPreviousCurve(I)) {
1172     if(previouscurve[I-1]>0) {
1173       theOpposite = Standard_False;
1174       return(previouscurve[I-1]);
1175     }
1176     else {
1177       theOpposite = Standard_True;
1178       return( - previouscurve[I-1]);
1179     }
1180   }
1181   else {
1182     throw Standard_DomainError("Incorrect Curve Number 'PreviousCurve'");
1183   }
1184 }
1185 //=======================================================================
1186 //function :NextCurve
1187 //purpose  : 
1188 //=======================================================================    
1189 Standard_Integer IntAna_IntQuadQuad::NextCurve (const Standard_Integer I,
1190                                                 Standard_Boolean& theOpposite) const
1191 {
1192   if(HasNextCurve(I)) {
1193     if(nextcurve[I]>0) {
1194       theOpposite = Standard_False;
1195       return(nextcurve[I-1]);
1196     }
1197     else {
1198       theOpposite = Standard_True;
1199       return( - nextcurve[I-1]);
1200     }
1201   }
1202   else {
1203     throw Standard_DomainError("Incorrect Curve Number 'NextCurve'");
1204   }
1205 }
1206 //=======================================================================
1207 //function :Curve
1208 //purpose  : 
1209 //=======================================================================       
1210 const IntAna_Curve& IntAna_IntQuadQuad::Curve(const Standard_Integer i) const
1211 {
1212   if(!done) {
1213     throw StdFail_NotDone("IntQuadQuad Not done");
1214   }
1215   if (identical) {
1216     throw Standard_DomainError("IntQuadQuad identical");
1217   }
1218   if((i <= 0) || (i>NbCurves)) {
1219     throw Standard_OutOfRange("Incorrect Curve Number");
1220   }
1221   return TheCurve[i-1];
1222 }
1223 //=======================================================================
1224 //function :Point
1225 //purpose  : 
1226 //=======================================================================    
1227 const gp_Pnt& IntAna_IntQuadQuad::Point (const Standard_Integer i) const
1228 {
1229   if(!done) {
1230     throw StdFail_NotDone("IntQuadQuad Not done");
1231   }  
1232   if (identical) {
1233     throw Standard_DomainError("IntQuadQuad identical");
1234   }
1235   if((i <= 0) || (i>Nbpoints)) {
1236     throw Standard_OutOfRange("Incorrect Point Number");
1237   }
1238   return Thepoints[i-1];
1239 }
1240 //=======================================================================
1241 //function :Parameters
1242 //purpose  : 
1243 //=======================================================================    
1244 void IntAna_IntQuadQuad::Parameters (const Standard_Integer, //i, 
1245                                      Standard_Real& , 
1246                                      Standard_Real& ) const
1247 {
1248   std::cout << "IntAna_IntQuadQuad::Parameters(...) is not yet implemented" << std::endl;
1249 }
1250
1251 /*********************************************************************************
1252
1253 Mathematica Pour Cone Quadrique 
1254 In[6]:= y[r_,t_]:=r*Sin[t]
1255
1256 In[7]:= x[r_,t_]:=r*Cos[t]
1257
1258 In[8]:= z[r_,t_]:=r*d
1259
1260 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
1261
1262 In[15]:= Quad[x[r,t],y[r,t],z[r,t]]
1263
1264                2      2        2       2        2       2
1265 Out[15]= Q1 + d  Qzz r  + Qxx r  Cos[t]  + Qyy r  Sin[t]  + 
1266  
1267                                       2
1268 >    2 (d Qz r + Qx r Cos[t] + d Qxz r  Cos[t] + Qy r Sin[t] + 
1269  
1270                2               2
1271 >       d Qyz r  Sin[t] + Qxy r  Cos[t] Sin[t])
1272
1273 In[16]:= QQ=%
1274
1275
1276
1277 In[17]:= Collect[QQ,r]
1278 Collect[QQ,r]
1279
1280 Out[17]= Q1 + r (2 d Qz + 2 Qx Cos[t] + 2 Qy Sin[t]) + 
1281  
1282       2   2                                  2
1283 >    r  (d  Qzz + 2 d Qxz Cos[t] + Qxx Cos[t]  + 2 d Qyz Sin[t] + 
1284  
1285                                         2
1286 >       2 Qxy Cos[t] Sin[t] + Qyy Sin[t] )
1287 ********************************************************************************/
1288
1289
1290   //**********************************************************************
1291   //***                   C O N E   - Q U A D R I Q U E                ***
1292   //***   2    2                                  2                    ***
1293   //***  R  ( d  Qzz + 2 d Qxz Cos[t] + Qxx Cos[t]  + 2 d Qyz Sin[t] + ***
1294   //***                                                                ***
1295   //***       2 Qxy Cos[t] Sin[t] + Qyy Sin[t] )                       ***
1296   //***                                                                ***
1297   //*** + R  ( 2 d Qz + 2 Qx Cos[t] + 2 Qy Sin[t] )                    ***
1298   //***                                                                ***
1299   //*** + Q1                                                           ***
1300   //**********************************************************************
1301   //FortranForm= ( DIS = QQ1 QQ1 - 4 QQ0 QQ2  ) / 4 
1302   //   -  d**2*Qz**2 - d**2*Qzz*Q1 + (Qx**2 - Qxx*Q1)*Cos(t)**2 + 
1303   //   -   (2*d*Qy*Qz - 2*d*Qyz*Q1)*Sin(t) + (Qy**2 - Qyy*Q1)*Sin(t)**2 + 
1304   //   -   Cos(t)*(2*d*Qx*Qz - 2*d*Qxz*Q1 + (2*Qx*Qy - 2*Qxy*Q1)*Sin(t))
1305   //**********************************************************************
1306 //modified by NIZNHY-PKV Fri Dec  2 10:56:03 2005f
1307 /*
1308 static
1309   void DumpCurve(const Standard_Integer aIndex,
1310                  IntAna_Curve& aC);
1311 //=======================================================================
1312 //function : DumpCurve
1313 //purpose  : 
1314 //=======================================================================
1315 void DumpCurve(const Standard_Integer aIndex,
1316                IntAna_Curve& aC)
1317 {
1318   Standard_Boolean bIsOpen, bIsConstant, bIsFirstOpen, bIsLastOpen;
1319   Standard_Integer i, aNb;
1320   Standard_Real aT1, aT2, aT, dT;
1321   gp_Pnt aP;
1322   //
1323   aC.Domain(aT1, aT2);
1324   bIsOpen=aC.IsOpen();
1325   bIsConstant=aC.IsConstant();
1326   bIsFirstOpen=aC.IsFirstOpen();
1327   bIsLastOpen=aC.IsLastOpen();
1328   //
1329   printf("\n");
1330   printf(" * IntAna_Curve #%d*\n", aIndex);
1331   printf(" Domain: [ %lf, %lf ]\n", aT1, aT2);
1332   printf(" IsOpen=%d\n", bIsOpen);
1333   printf(" IsConstant=%d\n", bIsConstant);
1334   printf(" IsFirstOpen =%d\n", bIsFirstOpen);
1335   printf(" IsLastOpen =%d\n", bIsLastOpen);
1336   //
1337   aNb=11;
1338   dT=(aT2-aT1)/(aNb-1);
1339   for (i=0; i<aNb; ++i) {
1340     aT=aT1+i*dT;
1341     if (i==(aNb-1)){
1342       aT=aT2;
1343     }
1344     //
1345     aP=aC.Value(aT);
1346     printf("point p%d_%d %lf %lf %lf\n", 
1347            aIndex, i, aP.X(), aP.Y(), aP.Z());
1348   }
1349   printf(" ---- end of curve ----\n");
1350 }
1351 */
1352 //modified by NIZNHY-PKV Fri Dec  2 10:42:16 2005t