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>
AddSpecialPoints(const IntAna_Quadric & theQuad,const gpSmth & theGpObj,Standard_Real & theTheta1,Standard_Real & theTheta2)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
IsDone()135 Standard_Boolean IsDone() {
136 return done;
137 }
138 //IsARoot
IsARoot(Standard_Real u)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
NbSolutions()155 Standard_Integer NbSolutions() {
156 if(!done) {
157 throw StdFail_NotDone();
158 }
159 return NbRoots;
160 }
161 //InfiniteRoots
InfiniteRoots()162 Standard_Boolean InfiniteRoots() {
163 if(!done) {
164 throw StdFail_NotDone();
165 }
166 return infinite_roots;
167 }
168 //Value
Value(const Standard_Integer n)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 //=======================================================================
TrigonometricRoots(const Standard_Real CC,const Standard_Real SC,const Standard_Real C,const Standard_Real S,const Standard_Real Cte,const Standard_Real Binf,const Standard_Real Bsup)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 //
MyTrigonometricFunction(const Standard_Real xCC,const Standard_Real xSS,const Standard_Real xSC,const Standard_Real xC,const Standard_Real xS,const Standard_Real xCte)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
Value(const Standard_Real & U)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 //
Derivative(const Standard_Real & U)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 //=======================================================================
IntAna_IntQuadQuad(void)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 //=======================================================================
IntAna_IntQuadQuad(const gp_Cylinder & Cyl,const IntAna_Quadric & Quad,const Standard_Real Tol)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 //=======================================================================
Perform(const gp_Cylinder & Cyl,const IntAna_Quadric & Quad,const Standard_Real)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 //=======================================================================
IntAna_IntQuadQuad(const gp_Cone & Cone,const IntAna_Quadric & Quad,const Standard_Real Tol)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 //=======================================================================
Perform(const gp_Cone & Cone,const IntAna_Quadric & Quad,const Standard_Real)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 in spite 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 //=======================================================================
InternalSetNextAndPrevious()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 //=======================================================================
HasPreviousCurve(const Standard_Integer I) const1128 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 //=======================================================================
HasNextCurve(const Standard_Integer I) const1148 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 //=======================================================================
PreviousCurve(const Standard_Integer I,Standard_Boolean & theOpposite) const1168 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 //=======================================================================
NextCurve(const Standard_Integer I,Standard_Boolean & theOpposite) const1189 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 //=======================================================================
Curve(const Standard_Integer i) const1210 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 //=======================================================================
Point(const Standard_Integer i) const1227 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 //=======================================================================
Parameters(const Standard_Integer,Standard_Real &,Standard_Real &) const1244 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
1353