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