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