projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0024511: Remove obsolete Image_PixelAddress.cxx
[occt.git] / src / IntAna / IntAna_QuadQuadGeo.cxx
CommitLineData
b311480e 1// Created on: 1992-08-06
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//
973c2be1 8// This library is free software; you can redistribute it and / or modify it
9// under the terms of the GNU Lesser General Public 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
b311480e 16
7fd59977 17//----------------------------------------------------------------------
18//-- Purposse: Geometric Intersection between two Natural Quadric
19//-- If the intersection is not a conic,
20//-- analytical methods must be called.
21//----------------------------------------------------------------------
22#ifndef DEB
23#define No_Standard_RangeError
24#define No_Standard_OutOfRange
25#endif
26
27#include <IntAna_QuadQuadGeo.ixx>
28
29#include <IntAna_IntConicQuad.hxx>
30#include <StdFail_NotDone.hxx>
31#include <Standard_DomainError.hxx>
32#include <Standard_OutOfRange.hxx>
33#include <math_DirectPolynomialRoots.hxx>
34
35#include <gp.hxx>
36#include <gp_Pln.hxx>
37#include <gp_Vec.hxx>
38#include <ElSLib.hxx>
39#include <ElCLib.hxx>
40
41#include <gp_Dir.hxx>
42#include <gp_XYZ.hxx>
43#include <gp_Pnt2d.hxx>
44#include <gp_Vec2d.hxx>
45#include <gp_Dir2d.hxx>
46
47
48static
49 gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D);
77088633 50static
51 void RefineDir(gp_Dir& aDir);
7fd59977 52
53//=======================================================================
54//class :
55//purpose : O p e r a t i o n s D i v e r s e s s u r d e s A x 1
56//=======================================================================
57class AxeOperator {
58 public:
59 AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2);
60
61 void Distance(Standard_Real& dist,
7eed5d29 62 Standard_Real& Param1,
63 Standard_Real& Param2);
7fd59977 64
65 gp_Pnt PtIntersect() {
66 return ptintersect;
67 }
68 Standard_Boolean Coplanar(void) {
69 return thecoplanar;
70 }
71 Standard_Boolean Same(void) {
72 return theparallel && (thedistance<myEPSILON_DISTANCE);
73 }
74 Standard_Real Distance(void) {
75 return thedistance ;
76 }
77 Standard_Boolean Intersect(void) {
78 return (thecoplanar && (!theparallel));
79 }
80 Standard_Boolean Parallel(void) {
81 return theparallel;
82 }
83 Standard_Boolean Normal(void) {
84 return thenormal;
85 }
86
87 protected:
88 Standard_Real Det33(const Standard_Real a11,
7eed5d29 89 const Standard_Real a12,
90 const Standard_Real a13,
91 const Standard_Real a21,
92 const Standard_Real a22,
93 const Standard_Real a23,
94 const Standard_Real a31,
95 const Standard_Real a32,
96 const Standard_Real a33) {
7fd59977 97 Standard_Real theReturn =
98 a11*(a22*a33-a32*a23) - a21*(a12*a33-a32*a13) + a31*(a12*a23-a22*a13) ;
99 return theReturn ;
100 }
101
102 private:
103 gp_Pnt ptintersect;
104 gp_Ax1 Axe1;
105 gp_Ax1 Axe2;
106 Standard_Real thedistance;
107 Standard_Boolean theparallel;
108 Standard_Boolean thecoplanar;
109 Standard_Boolean thenormal;
110 //
111 Standard_Real myEPSILON_DISTANCE;
112 Standard_Real myEPSILON_AXES_PARA;
113};
114
115//=======================================================================
116//function : AxeOperator::AxeOperator
117//purpose :
118//=======================================================================
119 AxeOperator::AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2)
120{
121 myEPSILON_DISTANCE=0.00000000000001;
122 myEPSILON_AXES_PARA=0.000000000001;
123 Axe1=A1;
124 Axe2=A2;
125 //---------------------------------------------------------------------
126 gp_Dir V1=Axe1.Direction();
127 gp_Dir V2=Axe2.Direction();
128 gp_Pnt P1=Axe1.Location();
129 gp_Pnt P2=Axe2.Location();
77088633 130 //
131 RefineDir(V1);
132 RefineDir(V2);
7fd59977 133 thecoplanar= Standard_False;
134 thenormal = Standard_False;
135
136 //--- check if the two axis are parallel
137 theparallel=V1.IsParallel(V2, myEPSILON_AXES_PARA);
138 //--- Distance between the two axis
139 gp_XYZ perp(A1.Direction().XYZ().Crossed(A2.Direction().XYZ()));
140 if (theparallel) {
141 gp_Lin L1(A1);
142 thedistance = L1.Distance(A2.Location());
143 }
144 else {
145 thedistance = Abs(gp_Vec(perp.Normalized()).Dot(gp_Vec(Axe1.Location(),
7eed5d29 146 Axe2.Location())));
7fd59977 147 }
148 //--- check if Axis are Coplanar
149 Standard_Real D33;
150 if(thedistance<myEPSILON_DISTANCE) {
151 D33=Det33(V1.X(),V1.Y(),V1.Z()
7eed5d29 152 ,V2.X(),V2.Y(),V2.Z()
153 ,P1.X()-P2.X(),P1.Y()-P2.Y(),P1.Z()-P2.Z());
7fd59977 154 if(Abs(D33)<=myEPSILON_DISTANCE) {
155 thecoplanar=Standard_True;
156 }
157 }
158 else {
159 thecoplanar=Standard_True;
160 thenormal=(V1.Dot(V2)==0.0)? Standard_True : Standard_False;
161 }
162 //--- check if the two axis are concurrent
163 if(thecoplanar && (!theparallel)) {
164 Standard_Real smx=P2.X() - P1.X();
165 Standard_Real smy=P2.Y() - P1.Y();
166 Standard_Real smz=P2.Z() - P1.Z();
167 Standard_Real Det1,Det2,Det3,A;
168 Det1=V1.Y() * V2.X() - V1.X() * V2.Y();
169 Det2=V1.Z() * V2.Y() - V1.Y() * V2.Z();
170 Det3=V1.Z() * V2.X() - V1.X() * V2.Z();
171
172 if((Det1!=0.0) && ((Abs(Det1) >= Abs(Det2))&&(Abs(Det1) >= Abs(Det3)))) {
173 A=(smy * V2.X() - smx * V2.Y())/Det1;
174 }
175 else if((Det2!=0.0)
7eed5d29 176 && ((Abs(Det2) >= Abs(Det1))
177 &&(Abs(Det2) >= Abs(Det3)))) {
7fd59977 178 A=(smz * V2.Y() - smy * V2.Z())/Det2;
179 }
180 else {
181 A=(smz * V2.X() - smx * V2.Z())/Det3;
182 }
183 ptintersect.SetCoord( P1.X() + A * V1.X()
7eed5d29 184 ,P1.Y() + A * V1.Y()
185 ,P1.Z() + A * V1.Z());
7fd59977 186 }
187 else {
188 ptintersect.SetCoord(0,0,0); //-- Pour eviter des FPE
189 }
190}
191//=======================================================================
192//function : Distance
193//purpose :
194//=======================================================================
195 void AxeOperator::Distance(Standard_Real& dist,Standard_Real& Param1,Standard_Real& Param2)
196 {
197 gp_Vec O1O2(Axe1.Location(),Axe2.Location()); //-----------------------------
198 gp_Dir U1 = Axe1.Direction(); //-- juste pour voir.
199 gp_Dir U2 = Axe2.Direction();
200
201 gp_Dir N = U1.Crossed(U2);
202 Standard_Real D = Det33(U1.X(),U2.X(),N.X(),
7eed5d29 203 U1.Y(),U2.Y(),N.Y(),
204 U1.Z(),U2.Z(),N.Z());
7fd59977 205 if(D) {
206 dist = Det33(U1.X(),U2.X(),O1O2.X(),
7eed5d29 207 U1.Y(),U2.Y(),O1O2.Y(),
208 U1.Z(),U2.Z(),O1O2.Z()) / D;
7fd59977 209 Param1 = Det33(O1O2.X(),U2.X(),N.X(),
7eed5d29 210 O1O2.Y(),U2.Y(),N.Y(),
211 O1O2.Z(),U2.Z(),N.Z()) / (-D);
7fd59977 212 //------------------------------------------------------------
213 //-- On resout P1 * Dir1 + P2 * Dir2 + d * N = O1O2
214 //-- soit : Segment perpendiculaire : O1+P1 D1
215 //-- O2-P2 D2
216 Param2 = Det33(U1.X(),O1O2.X(),N.X(),
7eed5d29 217 U1.Y(),O1O2.Y(),N.Y(),
218 U1.Z(),O1O2.Z(),N.Z()) / (D);
7fd59977 219 }
220}
221//=======================================================================
222//function : DirToAx2
223//purpose : returns a gp_Ax2 where D is the main direction
224//=======================================================================
225gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D)
226{
227 Standard_Real x=D.X(); Standard_Real ax=Abs(x);
228 Standard_Real y=D.Y(); Standard_Real ay=Abs(y);
229 Standard_Real z=D.Z(); Standard_Real az=Abs(z);
230 if( (ax==0.0) || ((ax<ay) && (ax<az)) ) {
231 return(gp_Ax2(P,D,gp_Dir(gp_Vec(0.0,-z,y))));
232 }
233 else if( (ay==0.0) || ((ay<ax) && (ay<az)) ) {
234 return(gp_Ax2(P,D,gp_Dir(gp_Vec(-z,0.0,x))));
235 }
236 else {
237 return(gp_Ax2(P,D,gp_Dir(gp_Vec(-y,x,0.0))));
238 }
239}
240//=======================================================================
241//function : IntAna_QuadQuadGeo
242//purpose : Empty constructor
243//=======================================================================
244 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(void)
245 : done(Standard_False),
246 nbint(0),
247 typeres(IntAna_Empty),
248 pt1(0,0,0),
249 pt2(0,0,0),
7eed5d29 250 pt3(0,0,0),
251 pt4(0,0,0),
7fd59977 252 param1(0),
253 param2(0),
7eed5d29 254 param3(0),
255 param4(0),
7fd59977 256 param1bis(0),
257 param2bis(0),
258 myCommonGen(Standard_False),
259 myPChar(0,0,0)
260{
261 InitTolerances();
262}
263//=======================================================================
264//function : InitTolerances
265//purpose :
266//=======================================================================
267 void IntAna_QuadQuadGeo::InitTolerances()
268{
269 myEPSILON_DISTANCE = 0.00000000000001;
270 myEPSILON_ANGLE_CONE = 0.000000000001;
271 myEPSILON_MINI_CIRCLE_RADIUS = 0.000000001;
272 myEPSILON_CYLINDER_DELTA_RADIUS = 0.0000000000001;
273 myEPSILON_CYLINDER_DELTA_DISTANCE= 0.0000001;
274 myEPSILON_AXES_PARA = 0.000000000001;
275}
276//=======================================================================
277//function : IntAna_QuadQuadGeo
278//purpose : Pln Pln
279//=======================================================================
280 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P1,
7eed5d29 281 const gp_Pln& P2,
282 const Standard_Real TolAng,
283 const Standard_Real Tol)
7fd59977 284: done(Standard_False),
285 nbint(0),
286 typeres(IntAna_Empty),
287 pt1(0,0,0),
288 pt2(0,0,0),
7eed5d29 289 pt3(0,0,0),
290 pt4(0,0,0),
7fd59977 291 param1(0),
292 param2(0),
7eed5d29 293 param3(0),
294 param4(0),
7fd59977 295 param1bis(0),
296 param2bis(0),
297 myCommonGen(Standard_False),
298 myPChar(0,0,0)
299{
300 InitTolerances();
301 Perform(P1,P2,TolAng,Tol);
302}
303//=======================================================================
304//function : Perform
305//purpose :
306//=======================================================================
307 void IntAna_QuadQuadGeo::Perform (const gp_Pln& P1,
7eed5d29 308 const gp_Pln& P2,
309 const Standard_Real TolAng,
310 const Standard_Real Tol)
7fd59977 311{
312 done=Standard_False;
313 //
314 param2bis=0.0;
315
316 Standard_Real A1 = 0., B1 = 0., C1 = 0., D1 = 0., A2 = 0., B2 = 0., C2 = 0., D2 = 0.;
317 P1.Coefficients(A1,B1,C1,D1);
318 P2.Coefficients(A2,B2,C2,D2);
319
320 gp_Vec vd(gp_Vec(A1,B1,C1).Crossed(gp_Vec(A2,B2,C2)));
321 Standard_Real dist1= A2*P1.Location().X() + B2*P1.Location().Y() + C2*P1.Location().Z() + D2;
322 Standard_Real dist2= A1*P2.Location().X() + B1*P2.Location().Y() + C1*P2.Location().Z() + D1;
323
324 if(vd.Magnitude() <=TolAng) {
325 // normalles are collinear - planes are same or parallel
326 typeres = (Abs(dist1) <= Tol && Abs(dist2) <= Tol) ? IntAna_Same : IntAna_Empty;
327 }
328 else {
329 Standard_Real denom=A1*A2 + B1*B2 + C1*C2;
330
331 Standard_Real denom2 = denom*denom;
332 Standard_Real ddenom = 1. - denom2;
4355f260 333 //denom = ( Abs(ddenom) <= 1.e-9 ) ? 1.e-9 : ddenom;
334 denom = ( Abs(ddenom) <= 1.e-16 ) ? 1.e-16 : ddenom;
7fd59977 335
336 Standard_Real par1 = dist1/denom;
337 Standard_Real par2 = -dist2/denom;
338
339 gp_Vec inter1(gp_Vec(A1,B1,C1).Crossed(vd));
340 gp_Vec inter2(gp_Vec(A2,B2,C2).Crossed(vd));
341
342 Standard_Real X1=P1.Location().X() + par1*inter1.X();
343 Standard_Real Y1=P1.Location().Y() + par1*inter1.Y();
344 Standard_Real Z1=P1.Location().Z() + par1*inter1.Z();
345 Standard_Real X2=P2.Location().X() + par2*inter2.X();
346 Standard_Real Y2=P2.Location().Y() + par2*inter2.Y();
347 Standard_Real Z2=P2.Location().Z() + par2*inter2.Z();
348
349 pt1=gp_Pnt((X1+X2)*0.5, (Y1+Y2)*0.5, (Z1+Z2)*0.5);
350 dir1 = gp_Dir(vd);
351 typeres = IntAna_Line;
352 nbint = 1;
353
354 }
355 done=Standard_True;
356}
357//=======================================================================
358//function : IntAna_QuadQuadGeo
359//purpose : Pln Cylinder
360//=======================================================================
361 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Pln& P
362 ,const gp_Cylinder& Cl
363 ,const Standard_Real Tolang
04cbc9d3 364 ,const Standard_Real Tol
365 ,const Standard_Real H)
7fd59977 366 : done(Standard_False),
367 nbint(0),
368 typeres(IntAna_Empty),
369 pt1(0,0,0),
370 pt2(0,0,0),
7eed5d29 371 pt3(0,0,0),
372 pt4(0,0,0),
7fd59977 373 param1(0),
374 param2(0),
7eed5d29 375 param3(0),
376 param4(0),
7fd59977 377 param1bis(0),
378 param2bis(0),
379 myCommonGen(Standard_False),
380 myPChar(0,0,0)
381{
382 InitTolerances();
04cbc9d3 383 Perform(P,Cl,Tolang,Tol,H);
7fd59977 384}
385//=======================================================================
386//function : Perform
387//purpose :
388//=======================================================================
389 void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
04cbc9d3 390 ,const gp_Cylinder& Cl
391 ,const Standard_Real Tolang
392 ,const Standard_Real Tol
393 ,const Standard_Real H)
7fd59977 394{
395 done = Standard_False;
396 Standard_Real dist,radius;
397 Standard_Real A,B,C,D;
398 Standard_Real X,Y,Z;
399 Standard_Real sint,cost,h;
400 gp_XYZ axex,axey,omega;
401
402
403 param2bis=0.0;
404 radius = Cl.Radius();
405
406 gp_Lin axec(Cl.Axis());
407 gp_XYZ normp(P.Axis().Direction().XYZ());
408
409 P.Coefficients(A,B,C,D);
410 axec.Location().Coord(X,Y,Z);
411 dist = A*X + B*Y + C*Z + D; // la distance axe/plan est evaluee a l origine de l axe.
412
413 Standard_Real tolang = Tolang;
414 Standard_Boolean newparams = Standard_False;
415
416 gp_Vec ldv( axec.Direction() );
417 gp_Vec npv( normp );
418 Standard_Real dA = Abs( ldv.Angle( npv ) );
c6541a0c 419 if( dA > (M_PI/4.) )
7fd59977 420 {
c6541a0c 421 Standard_Real dang = Abs( ldv.Angle( npv ) ) - M_PI/2.;
7fd59977 422 Standard_Real dangle = Abs( dang );
423 if( dangle > Tolang )
7eed5d29 424 {
425 Standard_Real sinda = Abs( Sin( dangle ) );
426 Standard_Real dif = Abs( sinda - Tol );
427 if( dif < Tol )
428 {
429 tolang = sinda * 2.;
430 newparams = Standard_True;
431 }
432 }
7fd59977 433 }
434
435 nbint = 0;
04cbc9d3 436 IntAna_IntConicQuad inter(axec,P,tolang,Tol,H);
7fd59977 437
438 if (inter.IsParallel()) {
439 // Le resultat de l intersection Plan-Cylindre est de type droite.
440 // il y a 1 ou 2 droites
441
442 typeres = IntAna_Line;
443 omega.SetCoord(X-dist*A,Y-dist*B,Z-dist*C);
444
445 if (Abs(Abs(dist)-radius) < Tol)
446 {
7eed5d29 447 nbint = 1;
448 pt1.SetXYZ(omega);
449
450 if( newparams )
451 {
452 gp_XYZ omegaXYZ(X,Y,Z);
453 gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
454 Standard_Real Xt,Yt,Zt,distt;
455 omegaXYZtrnsl.Coord(Xt,Yt,Zt);
456 distt = A*Xt + B*Yt + C*Zt + D;
457 gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A, omegaXYZtrnsl.Y()-distt*B, omegaXYZtrnsl.Z()-distt*C );
458 gp_Pnt ppt1;
459 ppt1.SetXYZ( omega1 );
460 gp_Vec vv1(pt1,ppt1);
461 gp_Dir dd1( vv1 );
462 dir1 = dd1;
463 }
464 else
465 dir1 = axec.Direction();
7fd59977 466 }
467 else if (Abs(dist) < radius)
468 {
7eed5d29 469 nbint = 2;
470 h = Sqrt(radius*radius - dist*dist);
471 axey = axec.Direction().XYZ().Crossed(normp); // axey est normalise
472
473 pt1.SetXYZ(omega - h*axey);
474 pt2.SetXYZ(omega + h*axey);
475
476 if( newparams )
477 {
478 gp_XYZ omegaXYZ(X,Y,Z);
479 gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
480 Standard_Real Xt,Yt,Zt,distt,ht;
481 omegaXYZtrnsl.Coord(Xt,Yt,Zt);
482 distt = A*Xt + B*Yt + C*Zt + D;
483 // ht = Sqrt(radius*radius - distt*distt);
484 Standard_Real anSqrtArg = radius*radius - distt*distt;
485 ht = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
486
487 gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A, omegaXYZtrnsl.Y()-distt*B, omegaXYZtrnsl.Z()-distt*C );
488 gp_Pnt ppt1,ppt2;
489 ppt1.SetXYZ( omega1 - ht*axey);
490 ppt2.SetXYZ( omega1 + ht*axey);
491 gp_Vec vv1(pt1,ppt1);
492 gp_Vec vv2(pt2,ppt2);
493 gp_Dir dd1( vv1 );
494 gp_Dir dd2( vv2 );
495 dir1 = dd1;
496 dir2 = dd2;
497 }
498 else
499 {
500 dir1 = axec.Direction();
501 dir2 = axec.Direction();
502 }
7fd59977 503 }
504 // else nbint = 0
505
506 // debug JAG : le nbint = 0 doit etre remplace par typeres = IntAna_Empty
507 // et ne pas etre seulement supprime...
508
509 else {
510 typeres = IntAna_Empty;
511 }
512 }
513 else { // Il y a un point d intersection. C est le centre du cercle
514 // ou de l ellipse solution.
515
516 nbint = 1;
517 axey = normp.Crossed(axec.Direction().XYZ());
518 sint = axey.Modulus();
519
520 pt1 = inter.Point(1);
521
522 if (sint < Tol/radius) {
523
524 // on construit un cercle avec comme axes X et Y ceux du cylindre
525 typeres = IntAna_Circle;
526
527 dir1 = axec.Direction(); // axe Z
528 dir2 = Cl.Position().XDirection();
529 param1 = radius;
530 }
531 else {
532
533 // on construit un ellipse
534 typeres = IntAna_Ellipse;
535 cost = Abs(axec.Direction().XYZ().Dot(normp));
536 axex = axey.Crossed(normp);
537
538 dir1.SetXYZ(normp); //Modif ds ce bloc
539 dir2.SetXYZ(axex);
540
541 param1 = radius/cost;
542 param1bis = radius;
543 }
544 }
788cbaf4 545
7fd59977 546 done = Standard_True;
547}
548//=======================================================================
549//function : IntAna_QuadQuadGeo
550//purpose : Pln Cone
551//=======================================================================
552 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
7eed5d29 553 const gp_Cone& Co,
554 const Standard_Real Tolang,
555 const Standard_Real Tol)
7fd59977 556: done(Standard_False),
557 nbint(0),
558 typeres(IntAna_Empty),
559 pt1(0,0,0),
560 pt2(0,0,0),
7eed5d29 561 pt3(0,0,0),
562 pt4(0,0,0),
7fd59977 563 param1(0),
564 param2(0),
7eed5d29 565 param3(0),
566 param4(0),
7fd59977 567 param1bis(0),
568 param2bis(0),
569 myCommonGen(Standard_False),
570 myPChar(0,0,0)
571{
572 InitTolerances();
573 Perform(P,Co,Tolang,Tol);
574}
575//=======================================================================
576//function : Perform
577//purpose :
578//=======================================================================
579 void IntAna_QuadQuadGeo::Perform(const gp_Pln& P,
7eed5d29 580 const gp_Cone& Co,
581 const Standard_Real Tolang,
582 const Standard_Real Tol)
7fd59977 583{
584
585 done = Standard_False;
586 nbint = 0;
587
588 Standard_Real A,B,C,D;
589 Standard_Real X,Y,Z;
590 Standard_Real dist,sint,cost,sina,cosa,angl,costa;
591 Standard_Real dh;
592 gp_XYZ axex,axey;
593
594 gp_Lin axec(Co.Axis());
595 P.Coefficients(A,B,C,D);
596 gp_Pnt apex(Co.Apex());
597
598 apex.Coord(X,Y,Z);
599 dist = A*X + B*Y + C*Z + D; // distance signee sommet du cone/ Plan
600
601 gp_XYZ normp = P.Axis().Direction().XYZ();
602 if(P.Direct()==Standard_False) { //-- lbr le 14 jan 97
603 normp.Reverse();
604 }
605
606 axey = normp.Crossed(Co.Axis().Direction().XYZ());
607 axex = axey.Crossed(normp);
608
609
610 angl = Co.SemiAngle();
611
612 cosa = Cos(angl);
613 sina = Abs(Sin(angl));
614
615
616 // Angle entre la normale au plan et l axe du cone, ramene entre 0. et PI/2.
617
618 sint = axey.Modulus();
619 cost = Abs(Co.Axis().Direction().XYZ().Dot(normp));
620
621 // Le calcul de costa permet de determiner si le plan contient
622 // un generatrice du cone : on calcul Sin((PI/2. - t) - angl)
623
624 costa = cost*cosa - sint*sina; // sin((PI/2 -t)-angl)=cos(t+angl)
625 // avec t ramene entre 0 et pi/2.
626
627 if (Abs(dist) < Tol) {
628 // on considere que le plan contient le sommet du cone.
629 // les solutions possibles sont donc : 1 point, 1 droite, 2 droites
630 // selon l inclinaison du plan.
631
632 if (Abs(costa) < Tolang) { // plan parallele a la generatrice
633 typeres = IntAna_Line;
634 nbint = 1;
635 gp_XYZ ptonaxe(apex.XYZ() + 10.*(Co.Axis().Direction().XYZ()));
636 // point sur l axe du cone cote z positif
637
638 dist = A*ptonaxe.X() + B*ptonaxe.Y() + C*ptonaxe.Z() + D;
639 ptonaxe = ptonaxe - dist*normp;
640 pt1 = apex;
641 dir1.SetXYZ(ptonaxe - pt1.XYZ());
642 }
643 else if (cost < sina) { // plan "interieur" au cone
644 typeres = IntAna_Line;
645 nbint = 2;
646 pt1 = apex;
647 pt2 = apex;
648 dh = Sqrt(sina*sina-cost*cost)/cosa;
649 dir1.SetXYZ(axex + dh*axey);
650 dir2.SetXYZ(axex - dh*axey);
651 }
652 else { // plan "exterieur" au cone
653 typeres = IntAna_Point;
654 nbint = 1;
655 pt1 = apex;
656 }
657 }
658 else {
659 // Solutions possibles : cercle, ellipse, parabole, hyperbole selon
660 // l inclinaison du plan.
661 Standard_Real deltacenter, distance;
662
663 if (cost < Tolang) {
664 // Le plan contient la direction de l axe du cone. La solution est
665 // l hyperbole
666 typeres = IntAna_Hyperbola;
667 nbint = 2;
668 pt1.SetXYZ(apex.XYZ()-dist*normp);
669 pt2 = pt1;
670 dir1=normp;
671 dir2.SetXYZ(axex);
672 param1 = param2 = Abs(dist/Tan(angl));
673 param1bis = param2bis = Abs(dist);
674 }
675 else {
676
677 IntAna_IntConicQuad inter(axec,P,Tolang); // on a necessairement 1 point.
678
679 gp_Pnt center(inter.Point(1));
680
681 // En fonction de la position de l intersection par rapport au sommet
682 // du cone, on change l axe x en -x et y en -y. Le parametre du sommet
683 // sur axec est negatif (voir definition du cone)
684
685 distance = apex.Distance(center);
686
687 if (inter.ParamOnConic(1) + Co.RefRadius()/Tan(angl) < 0.) {
7eed5d29 688 axex.Reverse();
689 axey.Reverse();
7fd59977 690 }
691
692 if (Abs(costa) < Tolang) { // plan parallele a une generatrice
7eed5d29 693 typeres = IntAna_Parabola;
694 nbint = 1;
695 deltacenter = distance/2./cosa;
696 axex.Normalize();
697 pt1.SetXYZ(center.XYZ()-deltacenter*axex);
698 dir1 = normp;
699 dir2.SetXYZ(axex);
700 param1 = deltacenter*sina*sina;
7fd59977 701 }
702 else if (sint < Tolang) { // plan perpendiculaire a l axe
7eed5d29 703 typeres = IntAna_Circle;
704 nbint = 1;
705 pt1 = center;
706 dir1 = Co.Position().Direction();
707 dir2 = Co.Position().XDirection();
708 param1 = apex.Distance(center)*Abs(Tan(angl));
7fd59977 709 }
710 else if (cost < sina ) {
7eed5d29 711 typeres = IntAna_Hyperbola;
712 nbint = 2;
713 axex.Normalize();
714
715 deltacenter = sint*sina*sina*distance/(sina*sina - cost*cost);
716 pt1.SetXYZ(center.XYZ() - deltacenter*axex);
717 pt2 = pt1;
718 dir1 = normp;
719 dir2.SetXYZ(axex);
720 param1 = param2 = cost*sina*cosa*distance /(sina*sina-cost*cost);
721 param1bis = param2bis = cost*sina*distance / Sqrt(sina*sina-cost*cost);
7fd59977 722
723 }
724 else { // on a alors cost > sina
7eed5d29 725 typeres = IntAna_Ellipse;
726 nbint = 1;
727 Standard_Real radius = cost*sina*cosa*distance/(cost*cost-sina*sina);
728 deltacenter = sint*sina*sina*distance/(cost*cost-sina*sina);
729 axex.Normalize();
730 pt1.SetXYZ(center.XYZ() + deltacenter*axex);
731 dir1 = normp;
732 dir2.SetXYZ(axex);
733 param1 = radius;
734 param1bis = cost*sina*distance/ Sqrt(cost*cost - sina*sina);
7fd59977 735 }
736 }
737 }
738
739 //-- On a du mal a gerer plus loin (Value ProjLib, Params ... )
740 //-- des hyperboles trop bizarres
741 //-- On retourne False -> Traitement par biparametree
742 static Standard_Real EllipseLimit = 1.0E+9; //OCC513(apo) 1000000
743 static Standard_Real HyperbolaLimit = 2.0E+6; //OCC537(apo) 50000
744 if(typeres==IntAna_Ellipse && nbint>=1) {
745 if(Abs(param1) > EllipseLimit || Abs(param1bis) > EllipseLimit) {
746 done=Standard_False;
747 return;
748 }
749 }
750 if(typeres==IntAna_Hyperbola && nbint>=2) {
751 if(Abs(param2) > HyperbolaLimit || Abs(param2bis) > HyperbolaLimit) {
752 done = Standard_False;
753 return;
754 }
755 }
756 if(typeres==IntAna_Hyperbola && nbint>=1) {
757 if(Abs(param1) > HyperbolaLimit || Abs(param1bis) > HyperbolaLimit) {
758 done=Standard_False;
759 return;
760 }
761 }
762
763 done = Standard_True;
764}
765
766//=======================================================================
767//function : IntAna_QuadQuadGeo
768//purpose : Pln Sphere
769//=======================================================================
770 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
7eed5d29 771 const gp_Sphere& S)
7fd59977 772: done(Standard_False),
773 nbint(0),
774 typeres(IntAna_Empty),
775 pt1(0,0,0),
776 pt2(0,0,0),
7eed5d29 777 pt3(0,0,0),
778 pt4(0,0,0),
7fd59977 779 param1(0),
780 param2(0),
7eed5d29 781 param3(0),
782 param4(0),
7fd59977 783 param1bis(0),
784 param2bis(0),
785 myCommonGen(Standard_False),
786 myPChar(0,0,0)
787{
788 InitTolerances();
789 Perform(P,S);
790}
791//=======================================================================
792//function : Perform
793//purpose :
794//=======================================================================
795 void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
7eed5d29 796 ,const gp_Sphere& S)
7fd59977 797{
798
799 done = Standard_False;
800 Standard_Real A,B,C,D,dist, radius;
801 Standard_Real X,Y,Z;
802
803 nbint = 0;
804// debug JAG : on met typeres = IntAna_Empty par defaut...
805 typeres = IntAna_Empty;
806
807 P.Coefficients(A,B,C,D);
808 S.Location().Coord(X,Y,Z);
809 radius = S.Radius();
810
811 dist = A * X + B * Y + C * Z + D;
812
813 if (Abs( Abs(dist) - radius) < Epsilon(radius)) {
814 // on a une seule solution : le point projection du centre de la sphere
815 // sur le plan
816 nbint = 1;
817 typeres = IntAna_Point;
818 pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
819 }
820 else if (Abs(dist) < radius) {
821 // on a un cercle solution
822 nbint = 1;
823 typeres = IntAna_Circle;
824 pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
825 dir1 = P.Axis().Direction();
826 if(P.Direct()==Standard_False) dir1.Reverse();
827 dir2 = P.Position().XDirection();
828 param1 = Sqrt(radius*radius - dist*dist);
829 }
830 param2bis=0.0; //-- pour eviter param2bis not used ....
831 done = Standard_True;
832}
833
834//=======================================================================
835//function : IntAna_QuadQuadGeo
836//purpose : Cylinder - Cylinder
837//=======================================================================
838 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl1,
7eed5d29 839 const gp_Cylinder& Cyl2,
840 const Standard_Real Tol)
7fd59977 841: done(Standard_False),
842 nbint(0),
843 typeres(IntAna_Empty),
844 pt1(0,0,0),
845 pt2(0,0,0),
7eed5d29 846 pt3(0,0,0),
847 pt4(0,0,0),
7fd59977 848 param1(0),
849 param2(0),
7eed5d29 850 param3(0),
851 param4(0),
7fd59977 852 param1bis(0),
853 param2bis(0),
854 myCommonGen(Standard_False),
855 myPChar(0,0,0)
856{
857 InitTolerances();
858 Perform(Cyl1,Cyl2,Tol);
859}
860//=======================================================================
861//function : Perform
862//purpose :
863//=======================================================================
864 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
7eed5d29 865 const gp_Cylinder& Cyl2,
866 const Standard_Real Tol)
7fd59977 867{
868 done=Standard_True;
869 //---------------------------- Parallel axes -------------------------
870 AxeOperator A1A2(Cyl1.Axis(),Cyl2.Axis());
871 Standard_Real R1=Cyl1.Radius();
872 Standard_Real R2=Cyl2.Radius();
873 Standard_Real RmR, RmR_Relative;
874 RmR=(R1>R2)? (R1-R2) : (R2-R1);
875 {
876 Standard_Real Rmax, Rmin;
877 Rmax=(R1>R2)? R1 : R2;
878 Rmin=(R1>R2)? R2 : R1;
879 RmR_Relative=RmR/Rmax;
880 }
881
882 Standard_Real DistA1A2=A1A2.Distance();
883
884 if(A1A2.Parallel()) {
885 if(DistA1A2<=Tol) {
886 if(RmR<=Tol) {
7eed5d29 887 typeres=IntAna_Same;
7fd59977 888 }
889 else {
7eed5d29 890 typeres=IntAna_Empty;
7fd59977 891 }
892 }
893 else { //-- DistA1A2 > Tol
894 gp_Pnt P1=Cyl1.Location();
895 gp_Pnt P2t=Cyl2.Location();
896 gp_Pnt P2;
897 //-- P2t is projected on the plane (P1,DirCylX,DirCylY)
898 gp_Dir DirCyl = Cyl1.Position().Direction();
899 Standard_Real ProjP2OnDirCyl1=gp_Vec(DirCyl).Dot(gp_Vec(P1,P2t));
900
901 P2.SetCoord( P2t.X() - ProjP2OnDirCyl1*DirCyl.X()
7eed5d29 902 ,P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y()
903 ,P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
7fd59977 904 //--
905 Standard_Real R1pR2=R1+R2;
906 if(DistA1A2>(R1pR2+Tol)) {
7eed5d29 907 typeres=IntAna_Empty;
908 nbint=0;
7fd59977 909 }
910 else if(DistA1A2>(R1pR2)) {
7eed5d29 911 //-- 1 Tangent line -------------------------------------OK
912 typeres=IntAna_Line;
913
914 nbint=1;
915 dir1=DirCyl;
916 Standard_Real R1_R1pR2=R1/R1pR2;
917 pt1.SetCoord( P1.X() + R1_R1pR2 * (P2.X()-P1.X())
918 ,P1.Y() + R1_R1pR2 * (P2.Y()-P1.Y())
919 ,P1.Z() + R1_R1pR2 * (P2.Z()-P1.Z()));
920
7fd59977 921 }
922 else if(DistA1A2>RmR) {
7eed5d29 923 //-- 2 lines ---------------------------------------------OK
924 typeres=IntAna_Line;
925 nbint=2;
926 dir1=DirCyl;
927 gp_Vec P1P2(P1,P2);
928 gp_Dir DirA1A2=gp_Dir(P1P2);
929 gp_Dir Ortho_dir1_P1P2 = dir1.Crossed(DirA1A2);
930 dir2=dir1;
931 Standard_Real Alpha=0.5*(R1*R1-R2*R2+DistA1A2*DistA1A2)/(DistA1A2);
932
933// Standard_Real Beta = Sqrt(R1*R1-Alpha*Alpha);
934 Standard_Real anSqrtArg = R1*R1-Alpha*Alpha;
935 Standard_Real Beta = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
936
937 if((Beta+Beta)<Tol) {
938 nbint=1;
939 pt1.SetCoord( P1.X() + Alpha*DirA1A2.X()
940 ,P1.Y() + Alpha*DirA1A2.Y()
941 ,P1.Z() + Alpha*DirA1A2.Z());
942 }
943 else {
944 pt1.SetCoord( P1.X() + Alpha*DirA1A2.X() + Beta*Ortho_dir1_P1P2.X()
945 ,P1.Y() + Alpha*DirA1A2.Y() + Beta*Ortho_dir1_P1P2.Y()
946 ,P1.Z() + Alpha*DirA1A2.Z() + Beta*Ortho_dir1_P1P2.Z() );
947
948 pt2.SetCoord( P1.X() + Alpha*DirA1A2.X() - Beta*Ortho_dir1_P1P2.X()
949 ,P1.Y() + Alpha*DirA1A2.Y() - Beta*Ortho_dir1_P1P2.Y()
950 ,P1.Z() + Alpha*DirA1A2.Z() - Beta*Ortho_dir1_P1P2.Z());
951 }
7fd59977 952 }
953 else if(DistA1A2>(RmR-Tol)) {
7eed5d29 954 //-- 1 Tangent ------------------------------------------OK
955 typeres=IntAna_Line;
956 nbint=1;
957 dir1=DirCyl;
958 Standard_Real R1_RmR=R1/RmR;
7fd59977 959
7eed5d29 960 if(R1 < R2) R1_RmR = -R1_RmR;
7fd59977 961
7eed5d29 962 pt1.SetCoord( P1.X() + R1_RmR * (P2.X()-P1.X())
963 ,P1.Y() + R1_RmR * (P2.Y()-P1.Y())
964 ,P1.Z() + R1_RmR * (P2.Z()-P1.Z()));
7fd59977 965 }
966 else {
7eed5d29 967 nbint=0;
968 typeres=IntAna_Empty;
7fd59977 969 }
970 }
971 }
972 else { //-- No Parallel Axis ---------------------------------OK
973 if((RmR_Relative<=myEPSILON_CYLINDER_DELTA_RADIUS)
974 && (DistA1A2 <= myEPSILON_CYLINDER_DELTA_DISTANCE)) {
975 //-- PI/2 between the two axis and Intersection
976 //-- and identical radius
977 typeres=IntAna_Ellipse;
978 nbint=2;
979 gp_Dir DirCyl1=Cyl1.Position().Direction();
980 gp_Dir DirCyl2=Cyl2.Position().Direction();
981 pt1=pt2=A1A2.PtIntersect();
982
983 Standard_Real A=DirCyl1.Angle(DirCyl2);
984 Standard_Real B;
c6541a0c 985 B=Abs(Sin(0.5*(M_PI-A)));
7fd59977 986 A=Abs(Sin(0.5*A));
987
988 if(A==0.0 || B==0.0) {
7eed5d29 989 typeres=IntAna_Same;
990 return;
7fd59977 991 }
992
993
994 gp_Vec dircyl1(DirCyl1);gp_Vec dircyl2(DirCyl2);
995 dir1 = gp_Dir(dircyl1.Added(dircyl2));
996 dir2 = gp_Dir(dircyl1.Subtracted(dircyl2));
7eed5d29 997
7fd59977 998 param2 = Cyl1.Radius() / A;
999 param1 = Cyl1.Radius() / B;
1000 param2bis= param1bis = Cyl1.Radius();
1001 if(param1 < param1bis) {
7eed5d29 1002 A=param1; param1=param1bis; param1bis=A;
7fd59977 1003 }
1004 if(param2 < param2bis) {
7eed5d29 1005 A=param2; param2=param2bis; param2bis=A;
7fd59977 1006 }
1007 }
1008 else {
1009 if(Abs(DistA1A2-Cyl1.Radius()-Cyl2.Radius())<Tol) {
7eed5d29 1010 typeres = IntAna_Point;
1011 Standard_Real d,p1,p2;
1012
1013 gp_Dir D1 = Cyl1.Axis().Direction();
1014 gp_Dir D2 = Cyl2.Axis().Direction();
1015 A1A2.Distance(d,p1,p2);
1016 gp_Pnt P = Cyl1.Axis().Location();
1017 gp_Pnt P1(P.X() - p1*D1.X(),
1018 P.Y() - p1*D1.Y(),
1019 P.Z() - p1*D1.Z());
1020 P = Cyl2.Axis().Location();
1021 gp_Pnt P2(P.X() - p2*D2.X(),
1022 P.Y() - p2*D2.Y(),
1023 P.Z() - p2*D2.Z());
1024 gp_Vec P1P2(P1,P2);
1025 D1=gp_Dir(P1P2);
1026 p1=Cyl1.Radius();
1027 pt1.SetCoord(P1.X() + p1*D1.X(),
1028 P1.Y() + p1*D1.Y(),
1029 P1.Z() + p1*D1.Z());
1030 nbint = 1;
7fd59977 1031 }
1032 else {
7eed5d29 1033 typeres=IntAna_NoGeometricSolution;
7fd59977 1034 }
1035 }
1036 }
1037}
1038//=======================================================================
1039//function : IntAna_QuadQuadGeo
1040//purpose : Cylinder - Cone
1041//=======================================================================
1042 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
7eed5d29 1043 const gp_Cone& Con,
1044 const Standard_Real Tol)
7fd59977 1045: done(Standard_False),
1046 nbint(0),
1047 typeres(IntAna_Empty),
1048 pt1(0,0,0),
1049 pt2(0,0,0),
7eed5d29 1050 pt3(0,0,0),
1051 pt4(0,0,0),
7fd59977 1052 param1(0),
1053 param2(0),
7eed5d29 1054 param3(0),
1055 param4(0),
7fd59977 1056 param1bis(0),
1057 param2bis(0),
1058 myCommonGen(Standard_False),
1059 myPChar(0,0,0)
1060{
1061 InitTolerances();
1062 Perform(Cyl,Con,Tol);
1063}
1064//=======================================================================
1065//function : Perform
1066//purpose :
1067//=======================================================================
1068 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
7eed5d29 1069 const gp_Cone& Con,
1070 const Standard_Real )
7fd59977 1071{
1072 done=Standard_True;
1073 AxeOperator A1A2(Cyl.Axis(),Con.Axis());
1074 if(A1A2.Same()) {
1075 gp_Pnt Pt=Con.Apex();
1076 Standard_Real dist=Cyl.Radius()/(Tan(Con.SemiAngle()));
1077 gp_Dir dir=Cyl.Position().Direction();
1078 pt1.SetCoord( Pt.X() + dist*dir.X()
7eed5d29 1079 ,Pt.Y() + dist*dir.Y()
1080 ,Pt.Z() + dist*dir.Z());
7fd59977 1081 pt2.SetCoord( Pt.X() - dist*dir.X()
7eed5d29 1082 ,Pt.Y() - dist*dir.Y()
1083 ,Pt.Z() - dist*dir.Z());
7fd59977 1084 dir1=dir2=dir;
1085 param1=param2=Cyl.Radius();
1086 nbint=2;
1087 typeres=IntAna_Circle;
1088
1089 }
1090 else {
1091 typeres=IntAna_NoGeometricSolution;
1092 }
1093}
1094//=======================================================================
1095//function :
1096//purpose : Cylinder - Sphere
1097//=======================================================================
1098 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
7eed5d29 1099 const gp_Sphere& Sph,
1100 const Standard_Real Tol)
7fd59977 1101: done(Standard_False),
1102 nbint(0),
1103 typeres(IntAna_Empty),
1104 pt1(0,0,0),
1105 pt2(0,0,0),
7eed5d29 1106 pt3(0,0,0),
1107 pt4(0,0,0),
7fd59977 1108 param1(0),
1109 param2(0),
7eed5d29 1110 param3(0),
1111 param4(0),
7fd59977 1112 param1bis(0),
1113 param2bis(0),
1114 myCommonGen(Standard_False),
1115 myPChar(0,0,0)
1116{
1117 InitTolerances();
1118 Perform(Cyl,Sph,Tol);
1119}
1120//=======================================================================
1121//function : Perform
1122//purpose :
1123//=======================================================================
1124 void IntAna_QuadQuadGeo::Perform( const gp_Cylinder& Cyl
7eed5d29 1125 ,const gp_Sphere& Sph
1126 ,const Standard_Real)
7fd59977 1127{
1128 done=Standard_True;
1129 gp_Pnt Pt=Sph.Location();
1130 AxeOperator A1A2(Cyl.Axis(),Sph.Position().Axis());
1131 if((A1A2.Intersect() && Pt.Distance(A1A2.PtIntersect())==0.0 )
1132 || (A1A2.Same())) {
1133 if(Sph.Radius() < Cyl.Radius()) {
1134 typeres = IntAna_Empty;
1135 }
1136 else {
1137 Standard_Real dist=Sqrt( Sph.Radius() * Sph.Radius() - Cyl.Radius() * Cyl.Radius() );
1138 gp_Dir dir=Cyl.Position().Direction();
1139 dir1 = dir2 = dir;
1140 typeres=IntAna_Circle;
1141 pt1.SetCoord( Pt.X() + dist*dir.X()
7eed5d29 1142 ,Pt.Y() + dist*dir.Y()
1143 ,Pt.Z() + dist*dir.Z());
7fd59977 1144 nbint=1;
1145 param1 = Cyl.Radius();
1146 if(dist>RealEpsilon()) {
7eed5d29 1147 pt2.SetCoord( Pt.X() - dist*dir.X()
1148 ,Pt.Y() - dist*dir.Y()
1149 ,Pt.Z() - dist*dir.Z());
1150 param2=Cyl.Radius();
1151 nbint=2;
7fd59977 1152 }
1153 }
1154 }
1155 else {
1156 typeres=IntAna_NoGeometricSolution;
1157 }
1158}
1159
1160//=======================================================================
1161//function : IntAna_QuadQuadGeo
1162//purpose : Cone - Cone
1163//=======================================================================
1164 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con1,
7eed5d29 1165 const gp_Cone& Con2,
1166 const Standard_Real Tol)
7fd59977 1167: done(Standard_False),
1168 nbint(0),
1169 typeres(IntAna_Empty),
1170 pt1(0,0,0),
1171 pt2(0,0,0),
7eed5d29 1172 pt3(0,0,0),
1173 pt4(0,0,0),
7fd59977 1174 param1(0),
1175 param2(0),
7eed5d29 1176 param3(0),
1177 param4(0),
7fd59977 1178 param1bis(0),
1179 param2bis(0),
1180 myCommonGen(Standard_False),
1181 myPChar(0,0,0)
1182{
1183 InitTolerances();
1184 Perform(Con1,Con2,Tol);
1185}
1186//
1187//=======================================================================
1188//function : Perform
1189//purpose :
1190//=======================================================================
1191 void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con1,
7eed5d29 1192 const gp_Cone& Con2,
1193 const Standard_Real Tol)
7fd59977 1194{
1195 done=Standard_True;
1196 //
1197 Standard_Real tg1, tg2, aDA1A2, aTol2;
1198 gp_Pnt aPApex1, aPApex2;
4bd102b8 1199
1200 Standard_Real TOL_APEX_CONF = 1.e-10;
1201
7fd59977 1202 //
1203 tg1=Tan(Con1.SemiAngle());
1204 tg2=Tan(Con2.SemiAngle());
1205
1206 if((tg1 * tg2) < 0.) {
1207 tg2 = -tg2;
1208 }
1209 //
7fd59977 1210 aTol2=Tol*Tol;
1211 aPApex1=Con1.Apex();
1212 aPApex2=Con2.Apex();
1213 aDA1A2=aPApex1.SquareDistance(aPApex2);
7fd59977 1214 //
1215 AxeOperator A1A2(Con1.Axis(),Con2.Axis());
1216 //
1217 // 1
1218 if(A1A2.Same()) {
1219 //-- two circles
1220 Standard_Real x;
1221 gp_Pnt P=Con1.Apex();
1222 gp_Dir D=Con1.Position().Direction();
1223 Standard_Real d=gp_Vec(D).Dot(gp_Vec(P,Con2.Apex()));
1224
1225 if(Abs(tg1-tg2)>myEPSILON_ANGLE_CONE) {
4bd102b8 1226 if (fabs(d) < TOL_APEX_CONF) {
7eed5d29 1227 typeres = IntAna_Point;
1228 nbint = 1;
1229 pt1 = P;
1230 return;
4bd102b8 1231 }
7fd59977 1232 x=(d*tg2)/(tg1+tg2);
1233 pt1.SetCoord( P.X() + x*D.X()
7eed5d29 1234 ,P.Y() + x*D.Y()
1235 ,P.Z() + x*D.Z());
7fd59977 1236 param1=Abs(x*tg1);
1237
1238 x=(d*tg2)/(tg2-tg1);
1239 pt2.SetCoord( P.X() + x*D.X()
7eed5d29 1240 ,P.Y() + x*D.Y()
1241 ,P.Z() + x*D.Z());
7fd59977 1242 param2=Abs(x*tg1);
1243 dir1 = dir2 = D;
1244 nbint=2;
1245 typeres=IntAna_Circle;
1246 }
1247 else {
4bd102b8 1248 if (fabs(d) < TOL_APEX_CONF) {
7eed5d29 1249 typeres=IntAna_Same;
7fd59977 1250 }
1251 else {
7eed5d29 1252 typeres=IntAna_Circle;
1253 nbint=1;
1254 x=d*0.5;
1255 pt1.SetCoord( P.X() + x*D.X()
1256 ,P.Y() + x*D.Y()
1257 ,P.Z() + x*D.Z());
1258 param1 = Abs(x * tg1);
1259 dir1 = D;
7fd59977 1260 }
1261 }
1262 } //-- fin A1A2.Same
1263 // 2
1264 else if((Abs(tg1-tg2)<myEPSILON_ANGLE_CONE) && (A1A2.Parallel())) {
1265 //-- voir AnVer12mai98
1266 Standard_Real DistA1A2=A1A2.Distance();
1267 gp_Dir DA1=Con1.Position().Direction();
1268 gp_Vec O1O2(Con1.Apex(),Con2.Apex());
b045e6a4 1269 gp_Dir O1O2n(O1O2); // normalization of the vector before projection
1270 Standard_Real O1O2_DA1=gp_Vec(DA1).Dot(gp_Vec(O1O2n));
1271
1272 gp_Vec O1_Proj_A2(O1O2n.X()-O1O2_DA1*DA1.X(),
7eed5d29 1273 O1O2n.Y()-O1O2_DA1*DA1.Y(),
1274 O1O2n.Z()-O1O2_DA1*DA1.Z());
7fd59977 1275 gp_Dir DB1=gp_Dir(O1_Proj_A2);
b045e6a4 1276
7fd59977 1277 Standard_Real yO1O2=O1O2.Dot(gp_Vec(DA1));
1278 Standard_Real ABSTG1 = Abs(tg1);
1279 Standard_Real X2 = (DistA1A2/ABSTG1 - yO1O2)*0.5;
1280 Standard_Real X1 = X2+yO1O2;
1281
1282 gp_Pnt P1(Con1.Apex().X() + X1*( DA1.X() + ABSTG1*DB1.X()),
7eed5d29 1283 Con1.Apex().Y() + X1*( DA1.Y() + ABSTG1*DB1.Y()),
1284 Con1.Apex().Z() + X1*( DA1.Z() + ABSTG1*DB1.Z()));
7fd59977 1285
1286 gp_Pnt MO1O2(0.5*(Con1.Apex().X()+Con2.Apex().X()),
7eed5d29 1287 0.5*(Con1.Apex().Y()+Con2.Apex().Y()),
1288 0.5*(Con1.Apex().Z()+Con2.Apex().Z()));
7fd59977 1289 gp_Vec P1MO1O2(P1,MO1O2);
1290
1291 gp_Dir DA1_X_DB1=DA1.Crossed(DB1);
1292 gp_Dir OrthoPln = DA1_X_DB1.Crossed(gp_Dir(P1MO1O2));
1293
1294 IntAna_QuadQuadGeo INTER_QUAD_PLN(gp_Pln(P1,OrthoPln),Con1,Tol,Tol);
1295 if(INTER_QUAD_PLN.IsDone()) {
1296 switch(INTER_QUAD_PLN.TypeInter()) {
7eed5d29 1297 case IntAna_Ellipse: {
1298 typeres=IntAna_Ellipse;
1299 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1300 pt1 = E.Location();
1301 dir1 = E.Position().Direction();
1302 dir2 = E.Position().XDirection();
1303 param1 = E.MajorRadius();
1304 param1bis = E.MinorRadius();
1305 nbint = 1;
1306 break;
7fd59977 1307 }
1308 case IntAna_Circle: {
7eed5d29 1309 typeres=IntAna_Circle;
1310 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1311 pt1 = C.Location();
1312 dir1 = C.Position().XDirection();
1313 dir2 = C.Position().YDirection();
1314 param1 = C.Radius();
1315 nbint = 1;
1316 break;
7fd59977 1317 }
1318 case IntAna_Hyperbola: {
7eed5d29 1319 typeres=IntAna_Hyperbola;
1320 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1321 pt1 = pt2 = H.Location();
1322 dir1 = H.Position().Direction();
1323 dir2 = H.Position().XDirection();
1324 param1 = param2 = H.MajorRadius();
1325 param1bis = param2bis = H.MinorRadius();
1326 nbint = 2;
1327 break;
7fd59977 1328 }
1329 case IntAna_Line: {
7eed5d29 1330 typeres=IntAna_Line;
1331 gp_Lin H=INTER_QUAD_PLN.Line(1);
1332 pt1 = pt2 = H.Location();
1333 dir1 = dir2 = H.Position().Direction();
1334 param1 = param2 = 0.0;
1335 param1bis = param2bis = 0.0;
1336 nbint = 2;
1337 break;
7fd59977 1338 }
1339 default:
7eed5d29 1340 typeres=IntAna_NoGeometricSolution;
7fd59977 1341 }
1342 }
1343 }// else if((Abs(tg1-tg2)<EPSILON_ANGLE_CONE) && (A1A2.Parallel()))
7fd59977 1344 // 3
1345 else if (aDA1A2<aTol2) {
7fd59977 1346 //
1347 // When apices are coinsided there can be 3 possible cases
1348 // 3.1 - empty solution (iRet=0)
1349 // 3.2 - one line when cone1 touches cone2 (iRet=1)
1350 // 3.3 - two lines when cone1 intersects cone2 (iRet=2)
1351 //
1352 Standard_Integer iRet;
1353 Standard_Real aGamma, aBeta1, aBeta2;
1354 Standard_Real aD1, aR1, aTgBeta1, aTgBeta2, aHalfPI;
1355 Standard_Real aCosGamma, aSinGamma, aDx, aR2, aRD2, aD2;
1356 gp_Pnt2d aP0, aPA1, aP1, aPA2;
1357 gp_Vec2d aVAx2;
1358 gp_Ax1 aAx1, aAx2;
1359 //
1360 // Preliminary analysis. Determination of iRet
1361 //
1362 iRet=0;
c6541a0c 1363 aHalfPI=0.5*M_PI;
7fd59977 1364 aD1=1.;
1365 aPA1.SetCoord(aD1, 0.);
1366 aP0.SetCoord(0., 0.);
1367 //
1368 aAx1=Con1.Axis();
1369 aAx2=Con2.Axis();
1370 aGamma=aAx1.Angle(aAx2);
1371 if (aGamma>aHalfPI){
c6541a0c 1372 aGamma=M_PI-aGamma;
7fd59977 1373 }
1374 aCosGamma=Cos(aGamma);
1375 aSinGamma=Sin(aGamma);
1376 //
1377 aBeta1=Con1.SemiAngle();
1378 aTgBeta1=Tan(aBeta1);
1379 aTgBeta1=Abs(aTgBeta1);
1380 //
1381 aBeta2=Con2.SemiAngle();
1382 aTgBeta2=Tan(aBeta2);
1383 aTgBeta2=Abs(aTgBeta2);
1384 //
1385 aR1=aD1*aTgBeta1;
1386 aP1.SetCoord(aD1, aR1);
1387 //
1388 // PA2
1389 aVAx2.SetCoord(aCosGamma, aSinGamma);
1390 gp_Dir2d aDAx2(aVAx2);
1391 gp_Lin2d aLAx2(aP0, aDAx2);
1392 //
1393 gp_Vec2d aV(aP0, aP1);
1394 aDx=aVAx2.Dot(aV);
1395 aPA2=aP0.Translated(aDx*aDAx2);
1396 //
1397 // aR2
1398 aDx=aPA2.Distance(aP0);
1399 aR2=aDx*aTgBeta2;
1400 //
1401 // aRD2
1402 aRD2=aPA2.Distance(aP1);
1403 //
1404 if (aRD2>(aR2+Tol)) {
1405 iRet=0;
7fd59977 1406 typeres=IntAna_Empty; //nothing
4101383e 1407 return;
7fd59977 1408 }
1409 //
1410 iRet=1; //touch case => 1 line
1411 if (aRD2<(aR2-Tol)) {
1412 iRet=2;//intersection => couple of lines
1413 }
1414 //
1415 // Finding the solution in 3D
1416 //
1417 Standard_Real aDa;
1418 gp_Pnt aQApex1, aQA1, aQA2, aQX, aQX1, aQX2;
1419 gp_Dir aD3Ax1, aD3Ax2;
1420 gp_Lin aLin;
1421 IntAna_QuadQuadGeo aIntr;
1422 //
1423 aQApex1=Con1.Apex();
1424 aD3Ax1=aAx1.Direction();
1425 aQA1.SetCoord(aQApex1.X()+aD1*aD3Ax1.X(),
7eed5d29 1426 aQApex1.Y()+aD1*aD3Ax1.Y(),
1427 aQApex1.Z()+aD1*aD3Ax1.Z());
7fd59977 1428 //
1429 aDx=aD3Ax1.Dot(aAx2.Direction());
1430 if (aDx<0.) {
1431 aAx2.Reverse();
1432 }
1433 aD3Ax2=aAx2.Direction();
1434 //
1435 aD2=aD1*sqrt((1.+aTgBeta1*aTgBeta1)/(1.+aTgBeta2*aTgBeta2));
1436 //
1437 aQA2.SetCoord(aQApex1.X()+aD2*aD3Ax2.X(),
7eed5d29 1438 aQApex1.Y()+aD2*aD3Ax2.Y(),
1439 aQApex1.Z()+aD2*aD3Ax2.Z());
7fd59977 1440 //
1441 gp_Pln aPln1(aQA1, aD3Ax1);
1442 gp_Pln aPln2(aQA2, aD3Ax2);
1443 //
1444 aIntr.Perform(aPln1, aPln2, Tol, Tol);
1445 if (!aIntr.IsDone()) {
1446 iRet=-1; // just in case. it must not be so
1447 typeres=IntAna_NoGeometricSolution;
1448 return;
1449 }
1450 //
1451 aLin=aIntr.Line(1);
1452 const gp_Dir& aDLin=aLin.Direction();
1453 gp_Vec aVLin(aDLin);
1454 gp_Pnt aOrig=aLin.Location();
1455 gp_Vec aVr(aQA1, aOrig);
1456 aDx=aVLin.Dot(aVr);
1457 aQX=aOrig.Translated(aDx*aVLin);
1458 //
1459 // Final part
1460 //
1461 typeres=IntAna_Line;
1462 //
1463 param1=0.;
1464 param2 =0.;
1465 param1bis=0.;
1466 param2bis=0.;
1467 //
1468 if (iRet==1) {
1469 // one line
1470 nbint=1;
1471 pt1=aQApex1;
1472 gp_Vec aVX(aQApex1, aQX);
1473 dir1=gp_Dir(aVX);
7fd59977 1474 }
1475
1476 else {//iRet=2
1477 // two lines
1478 nbint=2;
1479 aDa=aQA1.Distance(aQX);
1480 aDx=sqrt(aR1*aR1-aDa*aDa);
1481 aQX1=aQX.Translated(aDx*aVLin);
1482 aQX2=aQX.Translated(-aDx*aVLin);
1483 //
1484 pt1=aQApex1;
1485 pt2=aQApex1;
1486 gp_Vec aVX1(aQApex1, aQX1);
1487 dir1=gp_Dir(aVX1);
1488 gp_Vec aVX2(aQApex1, aQX2);
1489 dir2=gp_Dir(aVX2);
7fd59977 1490 }
1491 } //else if (aDA1A2<aTol2) {
7fd59977 1492 //Case when cones have common generatrix
1493 else if(A1A2.Intersect()) {
1494 //Check if apex of one cone belongs another one
1495 Standard_Real u, v, tol2 = Tol*Tol;
1496 ElSLib::Parameters(Con2, aPApex1, u, v);
1497 gp_Pnt p = ElSLib::Value(u, v, Con2);
1498 if(aPApex1.SquareDistance(p) > tol2) {
1499 typeres=IntAna_NoGeometricSolution;
1500 return;
1501 }
1502 //
1503 ElSLib::Parameters(Con1, aPApex2, u, v);
1504 p = ElSLib::Value(u, v, Con1);
1505 if(aPApex2.SquareDistance(p) > tol2) {
1506 typeres=IntAna_NoGeometricSolution;
1507 return;
1508 }
1509
1510 //Cones have a common generatrix passing through apexes
1511 myCommonGen = Standard_True;
1512
1513 //common generatrix of cones
1514 gp_Lin aGen(aPApex1, gp_Dir(gp_Vec(aPApex1, aPApex2)));
1515
1516 //Intersection point of axes
1517 gp_Pnt aPAxeInt = A1A2.PtIntersect();
1518
1519 //Characteristic point of intersection curve
1520 u = ElCLib::Parameter(aGen, aPAxeInt);
1521 myPChar = ElCLib::Value(u, aGen);
1522
1523
1524 //Other generatrixes of cones laying in maximal plane
c6541a0c
D		 1525 gp_Lin aGen1 = aGen.Rotated(Con1.Axis(), M_PI);
1526 gp_Lin aGen2 = aGen.Rotated(Con2.Axis(), M_PI);
7fd59977 1527 //
1528 //Intersection point of generatrixes
1529 gp_Dir aN; //solution plane normal
1530 gp_Dir aD1 = aGen1.Direction();
1531
1532 gp_Dir aD2(aD1.Crossed(aGen.Direction()));
1533
1534 if(aD1.IsParallel(aGen2.Direction(), Precision::Angular())) {
1535 aN = aD1.Crossed(aD2);
1536 }
1537 else if(aGen1.SquareDistance(aGen2) > tol2) {
1538 //Something wrong ???
1539 typeres=IntAna_NoGeometricSolution;
1540 return;
1541 }
1542 else {
1543 gp_Dir D1 = aGen1.Position().Direction();
1544 gp_Dir D2 = aGen2.Position().Direction();
1545 gp_Pnt O1 = aGen1.Location();
1546 gp_Pnt O2 = aGen2.Location();
1547 Standard_Real D1DotD2 = D1.Dot(D2);
1548 Standard_Real aSin = 1.-D1DotD2*D1DotD2;
1549 gp_Vec O1O2 (O1,O2);
1550 Standard_Real U2 = (D1.XYZ()*(O1O2.Dot(D1))-(O1O2.XYZ())).Dot(D2.XYZ());
1551 U2 /= aSin;
1552 gp_Pnt aPGint(ElCLib::Value(U2, aGen2));
1553
1554 aD1 = gp_Dir(gp_Vec(aPGint, myPChar));
1555 aN = aD1.Crossed(aD2);
1556 }
1557 //Plane that must contain intersection curves
1558 gp_Pln anIntPln(myPChar, aN);
1559
1560 IntAna_QuadQuadGeo INTER_QUAD_PLN(anIntPln,Con1,Tol,Tol);
1561
1562 if(INTER_QUAD_PLN.IsDone()) {
1563 switch(INTER_QUAD_PLN.TypeInter()) {
7eed5d29 1564 case IntAna_Ellipse: {
1565 typeres=IntAna_Ellipse;
1566 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1567 pt1 = E.Location();
1568 dir1 = E.Position().Direction();
1569 dir2 = E.Position().XDirection();
1570 param1 = E.MajorRadius();
1571 param1bis = E.MinorRadius();
1572 nbint = 1;
1573 break;
7fd59977 1574 }
1575 case IntAna_Circle: {
7eed5d29 1576 typeres=IntAna_Circle;
1577 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1578 pt1 = C.Location();
1579 dir1 = C.Position().XDirection();
1580 dir2 = C.Position().YDirection();
1581 param1 = C.Radius();
1582 nbint = 1;
1583 break;
7fd59977 1584 }
1585 case IntAna_Parabola: {
7eed5d29 1586 typeres=IntAna_Parabola;
1587 gp_Parab Prb=INTER_QUAD_PLN.Parabola(1);
1588 pt1 = Prb.Location();
1589 dir1 = Prb.Position().Direction();
1590 dir2 = Prb.Position().XDirection();
1591 param1 = Prb.Focal();
1592 nbint = 1;
1593 break;
7fd59977 1594 }
1595 case IntAna_Hyperbola: {
7eed5d29 1596 typeres=IntAna_Hyperbola;
1597 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1598 pt1 = pt2 = H.Location();
1599 dir1 = H.Position().Direction();
1600 dir2 = H.Position().XDirection();
1601 param1 = param2 = H.MajorRadius();
1602 param1bis = param2bis = H.MinorRadius();
1603 nbint = 2;
1604 break;
7fd59977 1605 }
1606 default:
7eed5d29 1607 typeres=IntAna_NoGeometricSolution;
7fd59977 1608 }
1609 }
1610 }
4101383e 1611
7fd59977 1612 else {
1613 typeres=IntAna_NoGeometricSolution;
1614 }
1615}
1616//=======================================================================
1617//function : IntAna_QuadQuadGeo
1618//purpose : Sphere - Cone
1619//=======================================================================
1620 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
7eed5d29 1621 const gp_Cone& Con,
1622 const Standard_Real Tol)
7fd59977 1623: done(Standard_False),
1624 nbint(0),
1625 typeres(IntAna_Empty),
1626 pt1(0,0,0),
1627 pt2(0,0,0),
7eed5d29 1628 pt3(0,0,0),
1629 pt4(0,0,0),
7fd59977 1630 param1(0),
1631 param2(0),
7eed5d29 1632 param3(0),
1633 param4(0),
7fd59977 1634 param1bis(0),
1635 param2bis(0),
1636 myCommonGen(Standard_False),
1637 myPChar(0,0,0)
1638{
1639 InitTolerances();
1640 Perform(Sph,Con,Tol);
1641}
1642//=======================================================================
1643//function : Perform
1644//purpose :
1645//=======================================================================
1646 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
7eed5d29 1647 const gp_Cone& Con,
1648 const Standard_Real)
7fd59977 1649{
77088633 1650
1651 //
7fd59977 1652 done=Standard_True;
77088633 1653 //
7fd59977 1654 AxeOperator A1A2(Con.Axis(),Sph.Position().Axis());
1655 gp_Pnt Pt=Sph.Location();
77088633 1656 //
7fd59977 1657 if((A1A2.Intersect() && (Pt.Distance(A1A2.PtIntersect())==0.0))
1658 || A1A2.Same()) {
1659 gp_Pnt ConApex= Con.Apex();
1660 Standard_Real dApexSphCenter=Pt.Distance(ConApex);
1661 gp_Dir ConDir;
1662 if(dApexSphCenter>RealEpsilon()) {
1663 ConDir = gp_Dir(gp_Vec(ConApex,Pt));
1664 }
1665 else {
1666 ConDir = Con.Position().Direction();
1667 }
1668
1669 Standard_Real Rad=Sph.Radius();
1670 Standard_Real tga=Tan(Con.SemiAngle());
1671
1672
1673 //-- 2 circles
1674 //-- x: Roots of (x**2 + y**2 = Rad**2)
1675 //-- tga = y / (x+dApexSphCenter)
1676 Standard_Real tgatga = tga * tga;
1677 math_DirectPolynomialRoots Eq( 1.0+tgatga
7eed5d29 1678 ,2.0*tgatga*dApexSphCenter
1679 ,-Rad*Rad + dApexSphCenter*dApexSphCenter*tgatga);
7fd59977 1680 if(Eq.IsDone()) {
1681 Standard_Integer nbsol=Eq.NbSolutions();
1682 if(nbsol==0) {
7eed5d29 1683 typeres=IntAna_Empty;
7fd59977 1684 }
1685 else {
7eed5d29 1686 typeres=IntAna_Circle;
1687 if(nbsol>=1) {
1688 Standard_Real x = Eq.Value(1);
1689 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1690 nbint=1;
1691 pt1.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1692 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1693 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1694 param1 = tga * dApexSphCenterpx;
1695 param1 = Abs(param1);
1696 dir1 = ConDir;
1697 if(param1<=myEPSILON_MINI_CIRCLE_RADIUS) {
1698 typeres=IntAna_PointAndCircle;
1699 param1=0.0;
1700 }
1701 }
1702 if(nbsol>=2) {
1703 Standard_Real x=Eq.Value(2);
1704 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1705 nbint=2;
1706 pt2.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1707 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1708 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1709 param2 = tga * dApexSphCenterpx;
1710 param2 = Abs(param2);
1711 dir2=ConDir;
1712 if(param2<=myEPSILON_MINI_CIRCLE_RADIUS) {
1713 typeres=IntAna_PointAndCircle;
1714 param2=0.0;
1715 }
1716 }
7fd59977 1717 }
1718 }
1719 else {
1720 done=Standard_False;
1721 }
1722 }
1723 else {
1724 typeres=IntAna_NoGeometricSolution;
1725 }
1726}
1727
1728//=======================================================================
1729//function : IntAna_QuadQuadGeo
1730//purpose : Sphere - Sphere
1731//=======================================================================
1732 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Sphere& Sph1
7eed5d29 1733 ,const gp_Sphere& Sph2
1734 ,const Standard_Real Tol)
7fd59977 1735: done(Standard_False),
1736 nbint(0),
1737 typeres(IntAna_Empty),
1738 pt1(0,0,0),
1739 pt2(0,0,0),
7eed5d29 1740 pt3(0,0,0),
1741 pt4(0,0,0),
7fd59977 1742 param1(0),
1743 param2(0),
7eed5d29 1744 param3(0),
1745 param4(0),
7fd59977 1746 param1bis(0),
1747 param2bis(0),
1748 myCommonGen(Standard_False),
1749 myPChar(0,0,0)
1750{
1751 InitTolerances();
1752 Perform(Sph1,Sph2,Tol);
1753}
1754//=======================================================================
1755//function : Perform
1756//purpose :
1757//=======================================================================
1758 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph1,
7eed5d29 1759 const gp_Sphere& Sph2,
1760 const Standard_Real Tol)
7fd59977 1761{
1762 done=Standard_True;
1763 gp_Pnt O1=Sph1.Location();
1764 gp_Pnt O2=Sph2.Location();
1765 Standard_Real dO1O2=O1.Distance(O2);
1766 Standard_Real R1=Sph1.Radius();
1767 Standard_Real R2=Sph2.Radius();
1768 Standard_Real Rmin,Rmax;
1769 typeres=IntAna_Empty;
1770 param2bis=0.0; //-- pour eviter param2bis not used ....
1771
1772 if(R1>R2) { Rmin=R2; Rmax=R1; } else { Rmin=R1; Rmax=R2; }
1773
1774 if(dO1O2<=Tol && (Abs(R1-R2) <= Tol)) {
1775 typeres = IntAna_Same;
1776 }
1777 else {
1778 if(dO1O2<=Tol) { return; }
1779 gp_Dir Dir=gp_Dir(gp_Vec(O1,O2));
1780 Standard_Real t = Rmax - dO1O2 - Rmin;
1781
1782 //----------------------------------------------------------------------
1783 //-- |----------------- R1 --------------------|
1784 //-- |----dO1O2-----|-----------R2----------|
1785 //-- --->--<-- t
1786 //--
1787 //-- |------ R1 ------|---------dO1O2----------|
1788 //-- |-------------------R2-----------------------|
1789 //-- --->--<-- t
1790 //----------------------------------------------------------------------
1791 if(t >= 0.0 && t <=Tol) {
1792 typeres = IntAna_Point;
1793 nbint = 1;
1794 Standard_Real t2;
1795 if(R1==Rmax) t2=(R1 + (R2 + dO1O2)) * 0.5;
1796 else t2=(-R1+(dO1O2-R2))*0.5;
7eed5d29 1797
7fd59977 1798 pt1.SetCoord( O1.X() + t2*Dir.X()
7eed5d29 1799 ,O1.Y() + t2*Dir.Y()
1800 ,O1.Z() + t2*Dir.Z());
7fd59977 1801 }
1802 else {
1803 //-----------------------------------------------------------------
1804 //-- |----------------- dO1O2 --------------------|
1805 //-- |----R1-----|-----------R2----------|-Tol-|
1806 //--
1807 //-- |----------------- Rmax --------------------|
1808 //-- |----Rmin----|-------dO1O2-------|-Tol-|
1809 //--
1810 //-----------------------------------------------------------------
1811 if((dO1O2 > (R1+R2+Tol)) || (Rmax > (dO1O2+Rmin+Tol))) {
7eed5d29 1812 typeres=IntAna_Empty;
7fd59977 1813 }
1814 else {
7eed5d29 1815 //---------------------------------------------------------------
1816 //--
1817 //--
1818 //---------------------------------------------------------------
1819 Standard_Real Alpha=0.5*(R1*R1-R2*R2+dO1O2*dO1O2)/(dO1O2);
1820 Standard_Real Beta = R1*R1-Alpha*Alpha;
1821 Beta = (Beta>0.0)? Sqrt(Beta) : 0.0;
1822
1823 if(Beta<= myEPSILON_MINI_CIRCLE_RADIUS) {
1824 typeres = IntAna_Point;
1825 Alpha = (R1 + (dO1O2 - R2)) * 0.5;
1826 }
1827 else {
1828 typeres = IntAna_Circle;
1829 dir1 = Dir;
1830 param1 = Beta;
1831 }
1832 pt1.SetCoord( O1.X() + Alpha*Dir.X()
1833 ,O1.Y() + Alpha*Dir.Y()
1834 ,O1.Z() + Alpha*Dir.Z());
1835
1836 nbint=1;
7fd59977 1837 }
1838 }
1839 }
1840}
7eed5d29 1841
1842//=======================================================================
1843//function : IntAna_QuadQuadGeo
1844//purpose : Plane - Torus
1845//=======================================================================
1846IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& Pln,
1847 const gp_Torus& Tor,
1848 const Standard_Real Tol)
1849: done(Standard_False),
1850 nbint(0),
1851 typeres(IntAna_Empty),
1852 pt1(0,0,0),
1853 pt2(0,0,0),
1854 pt3(0,0,0),
1855 pt4(0,0,0),
1856 param1(0),
1857 param2(0),
1858 param3(0),
1859 param4(0),
1860 param1bis(0),
1861 param2bis(0),
1862 myCommonGen(Standard_False),
1863 myPChar(0,0,0)
1864{
1865 InitTolerances();
1866 Perform(Pln,Tor,Tol);
1867}
1868//=======================================================================
1869//function : Perform
1870//purpose :
1871//=======================================================================
1872void IntAna_QuadQuadGeo::Perform(const gp_Pln& Pln,
1873 const gp_Torus& Tor,
1874 const Standard_Real Tol)
1875{
1876 done = Standard_True;
1877 //
1878 Standard_Real aRMin, aRMaj;
1879 //
1880 aRMin = Tor.MinorRadius();
1881 aRMaj = Tor.MajorRadius();
1882 if (aRMin >= aRMaj) {
1883 typeres = IntAna_NoGeometricSolution;
1884 return;
1885 }
1886 //
1887 const gp_Ax1 aPlnAx = Pln.Axis();
1888 const gp_Ax1 aTorAx = Tor.Axis();
1889 //
1890 Standard_Boolean bParallel, bNormal;
1891 //
1892 bParallel = aTorAx.IsParallel(aPlnAx, myEPSILON_AXES_PARA);
1893 bNormal = !bParallel ? aTorAx.IsNormal(aPlnAx, myEPSILON_AXES_PARA) : Standard_False;
1894 if (!bNormal && !bParallel) {
1895 typeres = IntAna_NoGeometricSolution;
1896 return;
1897 }
1898 //
1899 Standard_Real aDist;
1900 //
1901 gp_Pnt aTorLoc = aTorAx.Location();
1902 if (bParallel) {
1903 Standard_Real aDt, X, Y, Z, A, B, C, D;
1904 //
1905 Pln.Coefficients(A,B,C,D);
1906 aTorLoc.Coord(X,Y,Z);
1907 aDist = A*X + B*Y + C*Z + D;
1908 //
1909 if ((Abs(aDist) - aRMin) > Tol) {
1910 typeres=IntAna_Empty;
1911 return;
1912 }
1913 //
1914 typeres = IntAna_Circle;
1915 //
1916 pt1.SetCoord(X - aDist*A, Y - aDist*B, Z - aDist*C);
1917 aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
1918 param1 = aRMaj + aDt;
1919 dir1 = aTorAx.Direction();
1920 nbint = 1;
1921 if ((Abs(aDist) < aRMin) && (aDt > Tol)) {
1922 pt2 = pt1;
1923 param2 = aRMaj - aDt;
1924 dir2 = dir1;
1925 nbint = 2;
1926 }
1927 }
1928 //
1929 else {
1930 aDist = Pln.Distance(aTorLoc);
1931 if (aDist > myEPSILON_DISTANCE) {
1932 typeres = IntAna_NoGeometricSolution;
1933 return;
1934 }
1935 //
1936 typeres = IntAna_Circle;
1937 param2 = param1 = aRMin;
1938 dir2 = dir1 = aPlnAx.Direction();
1939 nbint = 2;
1940 //
1941 gp_Dir aDir = aTorAx.Direction()^dir1;
1942 pt1.SetXYZ(aTorLoc.XYZ() + aRMaj*aDir.XYZ());
1943 pt2.SetXYZ(aTorLoc.XYZ() - aRMaj*aDir.XYZ());
1944 }
1945}
1946
1947//=======================================================================
1948//function : IntAna_QuadQuadGeo
1949//purpose : Cylinder - Torus
1950//=======================================================================
1951IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1952 const gp_Torus& Tor,
1953 const Standard_Real Tol)
1954: done(Standard_False),
1955 nbint(0),
1956 typeres(IntAna_Empty),
1957 pt1(0,0,0),
1958 pt2(0,0,0),
1959 pt3(0,0,0),
1960 pt4(0,0,0),
1961 param1(0),
1962 param2(0),
1963 param3(0),
1964 param4(0),
1965 param1bis(0),
1966 param2bis(0),
1967 myCommonGen(Standard_False),
1968 myPChar(0,0,0)
1969{
1970 InitTolerances();
1971 Perform(Cyl,Tor,Tol);
1972}
1973//=======================================================================
1974//function : Perform
1975//purpose :
1976//=======================================================================
1977void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
1978 const gp_Torus& Tor,
1979 const Standard_Real Tol)
1980{
1981 done = Standard_True;
1982 //
1983 Standard_Real aRMin, aRMaj;
1984 //
1985 aRMin = Tor.MinorRadius();
1986 aRMaj = Tor.MajorRadius();
1987 if (aRMin >= aRMaj) {
1988 typeres = IntAna_NoGeometricSolution;
1989 return;
1990 }
1991 //
1992 const gp_Ax1 aCylAx = Cyl.Axis();
1993 const gp_Ax1 aTorAx = Tor.Axis();
1994 //
1995 const gp_Lin aLin(aTorAx);
1996 const gp_Pnt aLocCyl = Cyl.Location();
1997 //
1998 if (!aTorAx.IsParallel(aCylAx, myEPSILON_AXES_PARA) ||
1999 (aLin.Distance(aLocCyl) > myEPSILON_DISTANCE)) {
2000 typeres = IntAna_NoGeometricSolution;
2001 return;
2002 }
2003 //
2004 Standard_Real aRCyl;
2005 //
2006 aRCyl = Cyl.Radius();
2007 if (((aRCyl + Tol) < (aRMaj - aRMin)) || ((aRCyl - Tol) > (aRMaj + aRMin))) {
2008 typeres = IntAna_Empty;
2009 return;
2010 }
2011 //
2012 typeres = IntAna_Circle;
2013 //
2014 Standard_Real aDist = Sqrt(Abs(aRMin*aRMin - (aRCyl-aRMaj)*(aRCyl-aRMaj)));
2015 gp_XYZ aTorLoc = aTorAx.Location().XYZ();
2016 //
2017 dir1 = aTorAx.Direction();
2018 pt1.SetXYZ(aTorLoc + aDist*dir1.XYZ());
2019 param1 = aRCyl;
2020 nbint = 1;
2021 if ((aDist > Tol) && (aRCyl > (aRMaj - aRMin)) &&
2022 (aRCyl < (aRMaj + aRMin))) {
2023 dir2 = dir1;
2024 pt2.SetXYZ(aTorLoc - aDist*dir2.XYZ());
2025 param2 = param1;
2026 nbint = 2;
2027 }
2028}
2029
2030//=======================================================================
2031//function : IntAna_QuadQuadGeo
2032//purpose : Cone - Torus
2033//=======================================================================
2034IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con,
2035 const gp_Torus& Tor,
2036 const Standard_Real Tol)
2037: done(Standard_False),
2038 nbint(0),
2039 typeres(IntAna_Empty),
2040 pt1(0,0,0),
2041 pt2(0,0,0),
2042 pt3(0,0,0),
2043 pt4(0,0,0),
2044 param1(0),
2045 param2(0),
2046 param3(0),
2047 param4(0),
2048 param1bis(0),
2049 param2bis(0),
2050 myCommonGen(Standard_False),
2051 myPChar(0,0,0)
2052{
2053 InitTolerances();
2054 Perform(Con,Tor,Tol);
2055}
2056//=======================================================================
2057//function : Perform
2058//purpose :
2059//=======================================================================
2060void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con,
2061 const gp_Torus& Tor,
2062 const Standard_Real Tol)
2063{
2064 done = Standard_True;
2065 //
2066 Standard_Real aRMin, aRMaj;
2067 //
2068 aRMin = Tor.MinorRadius();
2069 aRMaj = Tor.MajorRadius();
2070 if (aRMin >= aRMaj) {
2071 typeres = IntAna_NoGeometricSolution;
2072 return;
2073 }
2074 //
2075 const gp_Ax1 aConAx = Con.Axis();
2076 const gp_Ax1 aTorAx = Tor.Axis();
2077 //
2078 const gp_Lin aLin(aTorAx);
2079 const gp_Pnt aConApex = Con.Apex();
2080 //
2081 if (!aTorAx.IsParallel(aConAx, myEPSILON_AXES_PARA) ||
2082 (aLin.Distance(aConApex) > myEPSILON_DISTANCE)) {
2083 typeres = IntAna_NoGeometricSolution;
2084 return;
2085 }
2086 //
2087 Standard_Real anAngle, aDist, aParam[4];
2088 Standard_Integer i;
2089 gp_Pnt aTorLoc, aPCT, aPN, aPt[4];
2090 gp_Dir aDir[4];
2091 //
2092 anAngle = Con.SemiAngle();
2093 aTorLoc = aTorAx.Location();
2094 //
2095 aPN.SetXYZ(aTorLoc.XYZ() + aRMaj*Tor.YAxis().Direction().XYZ());
2096 gp_Dir aDN (gp_Vec(aTorLoc, aPN));
2097 gp_Ax1 anAxCLRot(aConApex, aDN);
2098 gp_Lin aConL = aLin.Rotated(anAxCLRot, anAngle);
2099 gp_Dir aDL = aConL.Position().Direction();
2100 gp_Dir aXDir = Tor.XAxis().Direction();
2101 //
2102 typeres = IntAna_Empty;
2103 //
2104 for (i = 0; i < 2; ++i) {
2105 if (i) {
2106 aXDir.Reverse();
2107 }
2108 aPCT.SetXYZ(aTorLoc.XYZ() + aRMaj*aXDir.XYZ());
2109 //
2110 aDist = aConL.Distance(aPCT);
2111 if (aDist > aRMin+Tol) {
2112 continue;
2113 }
2114 //
2115 typeres = IntAna_Circle;
2116 //
2117 gp_XYZ aPh = aPCT.XYZ() - aDist*aConL.Normal(aPCT).Direction().XYZ();
2118 aDist = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
2119 //
2120 gp_Pnt aP;
2121 gp_XYZ aDVal = aDist*aDL.XYZ();
2122 aP.SetXYZ(aPh + aDVal);
2123 aParam[nbint] = aLin.Distance(aP);
2124 aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2125 aDir[nbint] = aTorAx.Direction();
2126 ++nbint;
2127 if ((aDist < aRMin) && (aDVal.Modulus() > Tol)) {
2128 aP.SetXYZ(aPh - aDVal);
2129 aParam[nbint] = aLin.Distance(aP);
2130 aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2131 aDir[nbint] = aDir[nbint-1];
2132 ++nbint;
2133 }
2134 }
2135 //
2136 for (i = 0; i < nbint; ++i) {
2137 switch (i) {
2138 case 0:{
2139 pt1 = aPt[i];
2140 param1 = aParam[i];
2141 dir1 = aDir[i];
2142 break;
2143 }
2144 case 1:{
2145 pt2 = aPt[i];
2146 param2 = aParam[i];
2147 dir2 = aDir[i];
2148 break;
2149 }
2150 case 2:{
2151 pt3 = aPt[i];
2152 param3 = aParam[i];
2153 dir3 = aDir[i];
2154 break;
2155 }
2156 case 3:{
2157 pt4 = aPt[i];
2158 param4 = aParam[i];
2159 dir4 = aDir[i];
2160 break;
2161 }
2162 default:
2163 break;
2164 }
2165 }
2166}
2167
2168//=======================================================================
2169//function : IntAna_QuadQuadGeo
2170//purpose : Sphere - Torus
2171//=======================================================================
2172IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
2173 const gp_Torus& Tor,
2174 const Standard_Real Tol)
2175: done(Standard_False),
2176 nbint(0),
2177 typeres(IntAna_Empty),
2178 pt1(0,0,0),
2179 pt2(0,0,0),
2180 pt3(0,0,0),
2181 pt4(0,0,0),
2182 param1(0),
2183 param2(0),
2184 param3(0),
2185 param4(0),
2186 param1bis(0),
2187 param2bis(0),
2188 myCommonGen(Standard_False),
2189 myPChar(0,0,0)
2190{
2191 InitTolerances();
2192 Perform(Sph,Tor,Tol);
2193}
2194//=======================================================================
2195//function : Perform
2196//purpose :
2197//=======================================================================
2198void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
2199 const gp_Torus& Tor,
2200 const Standard_Real Tol)
2201{
2202 done = Standard_True;
2203 //
2204 Standard_Real aRMin, aRMaj;
2205 //
2206 aRMin = Tor.MinorRadius();
2207 aRMaj = Tor.MajorRadius();
2208 if (aRMin >= aRMaj) {
2209 typeres = IntAna_NoGeometricSolution;
2210 return;
2211 }
2212 //
2213 const gp_Ax1 aTorAx = Tor.Axis();
2214 const gp_Lin aLin(aTorAx);
2215 const gp_Pnt aSphLoc = Sph.Location();
2216 //
2217 if (aLin.Distance(aSphLoc) > myEPSILON_DISTANCE) {
2218 typeres = IntAna_NoGeometricSolution;
2219 return;
2220 }
2221 //
2222 Standard_Real aRSph, aDist;
2223 gp_Pnt aTorLoc;
2224 //
2225 gp_Dir aXDir = Tor.XAxis().Direction();
2226 aTorLoc.SetXYZ(aTorAx.Location().XYZ() + aRMaj*aXDir.XYZ());
2227 aRSph = Sph.Radius();
2228 //
2229 gp_Vec aVec12(aTorLoc, aSphLoc);
2230 aDist = aVec12.Magnitude();
2231 if (((aDist - Tol) > (aRMin + aRSph)) ||
2232 ((aDist + Tol) < Abs(aRMin - aRSph))) {
2233 typeres = IntAna_Empty;
2234 return;
2235 }
2236 //
2237 typeres = IntAna_Circle;
2238 //
2239 Standard_Real anAlpha, aBeta;
2240 //
2241 anAlpha = 0.5*(aRMin*aRMin - aRSph*aRSph + aDist*aDist ) / aDist;
2242 aBeta = Sqrt(Abs(aRMin*aRMin - anAlpha*anAlpha));
2243 //
2244 gp_Dir aDir12(aVec12);
2245 gp_XYZ aPh = aTorLoc.XYZ() + anAlpha*aDir12.XYZ();
2246 gp_Dir aDC = Tor.YAxis().Direction()^aDir12;
2247 //
2248 gp_Pnt aP;
2249 gp_XYZ aDVal = aBeta*aDC.XYZ();
2250 aP.SetXYZ(aPh + aDVal);
2251 param1 = aLin.Distance(aP);
2252 pt1.SetXYZ(aP.XYZ() - param1*aXDir.XYZ());
2253 dir1 = aTorAx.Direction();
2254 nbint = 1;
2255 if ((aDist < (aRSph + aRMin)) && (aDist > Abs(aRSph - aRMin)) &&
2256 (aDVal.Modulus() > Tol)) {
2257 aP.SetXYZ(aPh - aDVal);
2258 param2 = aLin.Distance(aP);
2259 pt2.SetXYZ(aP.XYZ() - param2*aXDir.XYZ());
2260 dir2 = dir1;
2261 nbint = 2;
2262 }
2263}
2264
2265//=======================================================================
2266//function : IntAna_QuadQuadGeo
2267//purpose : Torus - Torus
2268//=======================================================================
2269IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Torus& Tor1,
2270 const gp_Torus& Tor2,
2271 const Standard_Real Tol)
2272: done(Standard_False),
2273 nbint(0),
2274 typeres(IntAna_Empty),
2275 pt1(0,0,0),
2276 pt2(0,0,0),
2277 pt3(0,0,0),
2278 pt4(0,0,0),
2279 param1(0),
2280 param2(0),
2281 param3(0),
2282 param4(0),
2283 param1bis(0),
2284 param2bis(0),
2285 myCommonGen(Standard_False),
2286 myPChar(0,0,0)
2287{
2288 InitTolerances();
2289 Perform(Tor1,Tor2,Tol);
2290}
2291//=======================================================================
2292//function : Perform
2293//purpose :
2294//=======================================================================
2295void IntAna_QuadQuadGeo::Perform(const gp_Torus& Tor1,
2296 const gp_Torus& Tor2,
2297 const Standard_Real Tol)
2298{
2299 done = Standard_True;
2300 //
2301 Standard_Real aRMin1, aRMin2, aRMaj1, aRMaj2;
2302 //
2303 aRMin1 = Tor1.MinorRadius();
2304 aRMaj1 = Tor1.MajorRadius();
2305 aRMin2 = Tor2.MinorRadius();
2306 aRMaj2 = Tor2.MajorRadius();
2307 if (aRMin1 >= aRMaj1 || aRMin2 >= aRMaj2) {
2308 typeres = IntAna_NoGeometricSolution;
2309 return;
2310 }
2311 //
2312 const gp_Ax1 anAx1 = Tor1.Axis();
2313 const gp_Ax1 anAx2 = Tor2.Axis();
2314 //
2315 gp_Lin aL1(anAx1);
2316 if (!anAx1.IsParallel(anAx2, myEPSILON_AXES_PARA) ||
2317 (aL1.Distance(anAx2.Location()) > myEPSILON_DISTANCE)) {
2318 typeres = IntAna_NoGeometricSolution;
2319 return;
2320 }
2321 //
2322 gp_Pnt aLoc1, aLoc2;
2323 //
2324 aLoc1 = anAx1.Location();
2325 aLoc2 = anAx2.Location();
2326 //
2327 if (aLoc1.IsEqual(aLoc2, Tol) &&
2328 (Abs(aRMin1 - aRMin2) <= Tol) &&
2329 (Abs(aRMaj1 - aRMaj2) <= Tol)) {
2330 typeres = IntAna_Same;
2331 return;
2332 }
2333 //
2334 Standard_Real aDist;
2335 gp_Pnt aP1, aP2;
2336 //
2337 gp_Dir aXDir1 = Tor1.XAxis().Direction();
2338 aP1.SetXYZ(aLoc1.XYZ() + aRMaj1*aXDir1.XYZ());
2339 aP2.SetXYZ(aLoc2.XYZ() + aRMaj2*aXDir1.XYZ());
2340 //
2341 gp_Vec aV12(aP1, aP2);
2342 aDist = aV12.Magnitude();
2343 if (((aDist - Tol) > (aRMin1 + aRMin2)) ||
2344 ((aDist + Tol) < Abs(aRMin1 - aRMin2))) {
2345 typeres = IntAna_Empty;
2346 return;
2347 }
2348 //
2349 typeres = IntAna_Circle;
2350 //
2351 Standard_Real anAlpha, aBeta;
2352 //
2353 anAlpha = 0.5*(aRMin1*aRMin1 - aRMin2*aRMin2 + aDist*aDist ) / aDist;
2354 aBeta = Sqrt(Abs(aRMin1*aRMin1 - anAlpha*anAlpha));
2355 //
2356 gp_Dir aDir12(aV12);
2357 gp_XYZ aPh = aP1.XYZ() + anAlpha*aDir12.XYZ();
2358 gp_Dir aDC = Tor1.YAxis().Direction()^aDir12;
2359 //
2360 gp_Pnt aP;
2361 gp_XYZ aDVal = aBeta*aDC.XYZ();
2362 aP.SetXYZ(aPh + aDVal);
2363 param1 = aL1.Distance(aP);
2364 pt1.SetXYZ(aP.XYZ() - param1*aXDir1.XYZ());
2365 dir1 = anAx1.Direction();
2366 nbint = 1;
2367 if ((aDist < (aRMin1 + aRMin2)) && (aDist > Abs(aRMin1 - aRMin2)) &&
2368 aDVal.Modulus() > Tol) {
2369 aP.SetXYZ(aPh - aDVal);
2370 param2 = aL1.Distance(aP);
2371 pt2.SetXYZ(aP.XYZ() - param2*aXDir1.XYZ());
2372 dir2 = dir1;
2373 nbint = 2;
2374 }
2375}
2376
7fd59977 2377//=======================================================================
2378//function : Point
2379//purpose : Returns a Point
2380//=======================================================================
2381 gp_Pnt IntAna_QuadQuadGeo::Point(const Standard_Integer n) const
2382{
2383 if(!done) { StdFail_NotDone::Raise(); }
2384 if(n>nbint || n<1) { Standard_DomainError::Raise(); }
2385 if(typeres==IntAna_PointAndCircle) {
2386 if(n!=1) { Standard_DomainError::Raise(); }
2387 if(param1==0.0) return(pt1);
2388 return(pt2);
2389 }
2390 else if(typeres==IntAna_Point) {
2391 if(n==1) return(pt1);
2392 return(pt2);
2393 }
2394
2395 // WNT (what can you expect from MicroSoft ?)
2396 return gp_Pnt(0,0,0);
2397}
2398//=======================================================================
2399//function : Line
2400//purpose : Returns a Line
2401//=======================================================================
2402 gp_Lin IntAna_QuadQuadGeo::Line(const Standard_Integer n) const
2403{
2404 if(!done) { StdFail_NotDone::Raise(); }
2405 if((n>nbint) || (n<1) || (typeres!=IntAna_Line)) {
2406 Standard_DomainError::Raise();
2407 }
2408 if(n==1) { return(gp_Lin(pt1,dir1)); }
2409 else { return(gp_Lin(pt2,dir2)); }
2410}
2411//=======================================================================
2412//function : Circle
2413//purpose : Returns a Circle
2414//=======================================================================
2415 gp_Circ IntAna_QuadQuadGeo::Circle(const Standard_Integer n) const
2416{
2417 if(!done) { StdFail_NotDone::Raise(); }
2418 if(typeres==IntAna_PointAndCircle) {
2419 if(n!=1) { Standard_DomainError::Raise(); }
2420 if(param2==0.0) return(gp_Circ(DirToAx2(pt1,dir1),param1));
2421 return(gp_Circ(DirToAx2(pt2,dir2),param2));
2422 }
2423 else if((n>nbint) || (n<1) || (typeres!=IntAna_Circle)) {
2424 Standard_DomainError::Raise();
2425 }
7eed5d29 2426 if (n==1) { return(gp_Circ(DirToAx2(pt1,dir1),param1));}
2427 else if (n==2) { return(gp_Circ(DirToAx2(pt2,dir2),param2));}
2428 else if (n==3) { return(gp_Circ(DirToAx2(pt3,dir3),param3));}
2429 else { return(gp_Circ(DirToAx2(pt4,dir4),param4));}
7fd59977 2430}
2431
2432//=======================================================================
2433//function : Ellipse
2434//purpose : Returns a Elips
2435//=======================================================================
2436 gp_Elips IntAna_QuadQuadGeo::Ellipse(const Standard_Integer n) const
2437{
2438 if(!done) { StdFail_NotDone::Raise(); }
2439 if((n>nbint) || (n<1) || (typeres!=IntAna_Ellipse)) {
2440 Standard_DomainError::Raise();
2441 }
2442
2443 if(n==1) {
2444 Standard_Real R1=param1, R2=param1bis, aTmp;
2445 if (R1<R2) {
2446 aTmp=R1; R1=R2; R2=aTmp;
2447 }
2448 gp_Ax2 anAx2(pt1, dir1 ,dir2);
2449 gp_Elips anElips (anAx2, R1, R2);
2450 return anElips;
2451 }
2452 else {
2453 Standard_Real R1=param2, R2=param2bis, aTmp;
2454 if (R1<R2) {
2455 aTmp=R1; R1=R2; R2=aTmp;
2456 }
2457 gp_Ax2 anAx2(pt2, dir2 ,dir1);
2458 gp_Elips anElips (anAx2, R1, R2);
2459 return anElips;
2460 }
2461}
2462//=======================================================================
2463//function : Parabola
2464//purpose : Returns a Parabola
2465//=======================================================================
2466 gp_Parab IntAna_QuadQuadGeo::Parabola(const Standard_Integer n) const
2467{
2468 if(!done) {
2469 StdFail_NotDone::Raise();
2470 }
2471 if (typeres!=IntAna_Parabola) {
2472 Standard_DomainError::Raise();
2473 }
2474 if((n>nbint) || (n!=1)) {
2475 Standard_OutOfRange::Raise();
2476 }
2477 return(gp_Parab(gp_Ax2( pt1
7eed5d29 2478 ,dir1
2479 ,dir2)
2480 ,param1));
7fd59977 2481}
2482//=======================================================================
2483//function : Hyperbola
2484//purpose : Returns a Hyperbola
2485//=======================================================================
2486 gp_Hypr IntAna_QuadQuadGeo::Hyperbola(const Standard_Integer n) const
2487{
2488 if(!done) {
2489 StdFail_NotDone::Raise();
2490 }
2491 if((n>nbint) || (n<1) || (typeres!=IntAna_Hyperbola)) {
2492 Standard_DomainError::Raise();
2493 }
2494 if(n==1) {
2495 return(gp_Hypr(gp_Ax2( pt1
7eed5d29 2496 ,dir1
2497 ,dir2)
2498 ,param1,param1bis));
7fd59977 2499 }
2500 else {
2501 return(gp_Hypr(gp_Ax2( pt2
7eed5d29 2502 ,dir1
2503 ,dir2.Reversed())
2504 ,param2,param2bis));
7fd59977 2505 }
2506}
7fd59977 2507//=======================================================================
2508//function : HasCommonGen
2509//purpose :
2510//=======================================================================
7fd59977 2511Standard_Boolean IntAna_QuadQuadGeo::HasCommonGen() const
2512{
2513 return myCommonGen;
2514}
7fd59977 2515//=======================================================================
2516//function : PChar
2517//purpose :
2518//=======================================================================
7fd59977 2519const gp_Pnt& IntAna_QuadQuadGeo::PChar() const
2520{
2521 return myPChar;
2522}
77088633 2523//=======================================================================
2524//function : RefineDir
2525//purpose :
2526//=======================================================================
2527void RefineDir(gp_Dir& aDir)
2528{
2529 Standard_Integer k, m, n;
2530 Standard_Real aC[3];
2531 //
2532 aDir.Coord(aC[0], aC[1], aC[2]);
2533 //
2534 m=0;
2535 n=0;
2536 for (k=0; k<3; ++k) {
2537 if (aC[k]==1. || aC[k]==-1.) {
2538 ++m;
2539 }
2540 else if (aC[k]!=0.) {
2541 ++n;
2542 }
2543 }
2544 //
2545 if (m && n) {
2546 Standard_Real aEps, aR1, aR2, aNum;
2547 //
2548 aEps=RealEpsilon();
2549 aR1=1.-aEps;
2550 aR2=1.+aEps;
2551 //
2552 for (k=0; k<3; ++k) {
2553 m=(aC[k]>0.);
2554 aNum=(m)? aC[k] : -aC[k];
2555 if (aNum>aR1 && aNum<aR2) {
7eed5d29 2556 if (m) {
2557 aC[k]=1.;
2558 }
2559 else {
2560 aC[k]=-1.;
2561 }
2562 //
2563 aC[(k+1)%3]=0.;
2564 aC[(k+2)%3]=0.;
2565 break;
77088633 2566 }
2567 }
2568 aDir.SetCoord(aC[0], aC[1], aC[2]);
2569 }
2570}
7fd59977 2571
2572
2573
Open CASCADE Technology repository