projects / occt.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
0024142: Wrong section curve
[occt.git] / src / IntAna / IntAna_QuadQuadGeo.cxx
... / ...
CommitLineData
1// Created on: 1992-08-06
2// Created by: Laurent BUCHARD
3// Copyright (c) 1992-1999 Matra Datavision
4// Copyright (c) 1999-2012 OPEN CASCADE SAS
5//
6// The content of this file is subject to the Open CASCADE Technology Public
7// License Version 6.5 (the "License"). You may not use the content of this file
8// except in compliance with the License. Please obtain a copy of the License
9// at http://www.opencascade.org and read it completely before using this file.
10//
11// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13//
14// The Original Code and all software distributed under the License is
15// distributed on an "AS IS" basis, without warranty of any kind, and the
16// Initial Developer hereby disclaims all such warranties, including without
17// limitation, any warranties of merchantability, fitness for a particular
18// purpose or non-infringement. Please see the License for the specific terms
19// and conditions governing the rights and limitations under the License.
20
21//----------------------------------------------------------------------
22//-- Purposse: Geometric Intersection between two Natural Quadric
23//-- If the intersection is not a conic,
24//-- analytical methods must be called.
25//----------------------------------------------------------------------
26#ifndef DEB
27#define No_Standard_RangeError
28#define No_Standard_OutOfRange
29#endif
30
31#include <IntAna_QuadQuadGeo.ixx>
32
33#include <IntAna_IntConicQuad.hxx>
34#include <StdFail_NotDone.hxx>
35#include <Standard_DomainError.hxx>
36#include <Standard_OutOfRange.hxx>
37#include <math_DirectPolynomialRoots.hxx>
38
39#include <gp.hxx>
40#include <gp_Pln.hxx>
41#include <gp_Vec.hxx>
42#include <ElSLib.hxx>
43#include <ElCLib.hxx>
44
45#include <gp_Dir.hxx>
46#include <gp_XYZ.hxx>
47#include <gp_Pnt2d.hxx>
48#include <gp_Vec2d.hxx>
49#include <gp_Dir2d.hxx>
50
51
52static
53 gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D);
54static
55 void RefineDir(gp_Dir& aDir);
56
57//=======================================================================
58//class :
59//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
60//=======================================================================
61class AxeOperator {
62 public:
63 AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2);
64
65 void Distance(Standard_Real& dist,
66 Standard_Real& Param1,
67 Standard_Real& Param2);
68
69 gp_Pnt PtIntersect() {
70 return ptintersect;
71 }
72 Standard_Boolean Coplanar(void) {
73 return thecoplanar;
74 }
75 Standard_Boolean Same(void) {
76 return theparallel && (thedistance<myEPSILON_DISTANCE);
77 }
78 Standard_Real Distance(void) {
79 return thedistance ;
80 }
81 Standard_Boolean Intersect(void) {
82 return (thecoplanar && (!theparallel));
83 }
84 Standard_Boolean Parallel(void) {
85 return theparallel;
86 }
87 Standard_Boolean Normal(void) {
88 return thenormal;
89 }
90
91 protected:
92 Standard_Real Det33(const Standard_Real a11,
93 const Standard_Real a12,
94 const Standard_Real a13,
95 const Standard_Real a21,
96 const Standard_Real a22,
97 const Standard_Real a23,
98 const Standard_Real a31,
99 const Standard_Real a32,
100 const Standard_Real a33) {
101 Standard_Real theReturn =
102 a11*(a22*a33-a32*a23) - a21*(a12*a33-a32*a13) + a31*(a12*a23-a22*a13) ;
103 return theReturn ;
104 }
105
106 private:
107 gp_Pnt ptintersect;
108 gp_Ax1 Axe1;
109 gp_Ax1 Axe2;
110 Standard_Real thedistance;
111 Standard_Boolean theparallel;
112 Standard_Boolean thecoplanar;
113 Standard_Boolean thenormal;
114 //
115 Standard_Real myEPSILON_DISTANCE;
116 Standard_Real myEPSILON_AXES_PARA;
117};
118
119//=======================================================================
120//function : AxeOperator::AxeOperator
121//purpose :
122//=======================================================================
123 AxeOperator::AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2)
124{
125 myEPSILON_DISTANCE=0.00000000000001;
126 myEPSILON_AXES_PARA=0.000000000001;
127 Axe1=A1;
128 Axe2=A2;
129 //---------------------------------------------------------------------
130 gp_Dir V1=Axe1.Direction();
131 gp_Dir V2=Axe2.Direction();
132 gp_Pnt P1=Axe1.Location();
133 gp_Pnt P2=Axe2.Location();
134 //
135 RefineDir(V1);
136 RefineDir(V2);
137 thecoplanar= Standard_False;
138 thenormal = Standard_False;
139
140 //--- check if the two axis are parallel
141 theparallel=V1.IsParallel(V2, myEPSILON_AXES_PARA);
142 //--- Distance between the two axis
143 gp_XYZ perp(A1.Direction().XYZ().Crossed(A2.Direction().XYZ()));
144 if (theparallel) {
145 gp_Lin L1(A1);
146 thedistance = L1.Distance(A2.Location());
147 }
148 else {
149 thedistance = Abs(gp_Vec(perp.Normalized()).Dot(gp_Vec(Axe1.Location(),
150 Axe2.Location())));
151 }
152 //--- check if Axis are Coplanar
153 Standard_Real D33;
154 if(thedistance<myEPSILON_DISTANCE) {
155 D33=Det33(V1.X(),V1.Y(),V1.Z()
156 ,V2.X(),V2.Y(),V2.Z()
157 ,P1.X()-P2.X(),P1.Y()-P2.Y(),P1.Z()-P2.Z());
158 if(Abs(D33)<=myEPSILON_DISTANCE) {
159 thecoplanar=Standard_True;
160 }
161 }
162 else {
163 thecoplanar=Standard_True;
164 thenormal=(V1.Dot(V2)==0.0)? Standard_True : Standard_False;
165 }
166 //--- check if the two axis are concurrent
167 if(thecoplanar && (!theparallel)) {
168 Standard_Real smx=P2.X() - P1.X();
169 Standard_Real smy=P2.Y() - P1.Y();
170 Standard_Real smz=P2.Z() - P1.Z();
171 Standard_Real Det1,Det2,Det3,A;
172 Det1=V1.Y() * V2.X() - V1.X() * V2.Y();
173 Det2=V1.Z() * V2.Y() - V1.Y() * V2.Z();
174 Det3=V1.Z() * V2.X() - V1.X() * V2.Z();
175
176 if((Det1!=0.0) && ((Abs(Det1) >= Abs(Det2))&&(Abs(Det1) >= Abs(Det3)))) {
177 A=(smy * V2.X() - smx * V2.Y())/Det1;
178 }
179 else if((Det2!=0.0)
180 && ((Abs(Det2) >= Abs(Det1))
181 &&(Abs(Det2) >= Abs(Det3)))) {
182 A=(smz * V2.Y() - smy * V2.Z())/Det2;
183 }
184 else {
185 A=(smz * V2.X() - smx * V2.Z())/Det3;
186 }
187 ptintersect.SetCoord( P1.X() + A * V1.X()
188 ,P1.Y() + A * V1.Y()
189 ,P1.Z() + A * V1.Z());
190 }
191 else {
192 ptintersect.SetCoord(0,0,0); //-- Pour eviter des FPE
193 }
194}
195//=======================================================================
196//function : Distance
197//purpose :
198//=======================================================================
199 void AxeOperator::Distance(Standard_Real& dist,Standard_Real& Param1,Standard_Real& Param2)
200 {
201 gp_Vec O1O2(Axe1.Location(),Axe2.Location()); //-----------------------------
202 gp_Dir U1 = Axe1.Direction(); //-- juste pour voir.
203 gp_Dir U2 = Axe2.Direction();
204
205 gp_Dir N = U1.Crossed(U2);
206 Standard_Real D = Det33(U1.X(),U2.X(),N.X(),
207 U1.Y(),U2.Y(),N.Y(),
208 U1.Z(),U2.Z(),N.Z());
209 if(D) {
210 dist = Det33(U1.X(),U2.X(),O1O2.X(),
211 U1.Y(),U2.Y(),O1O2.Y(),
212 U1.Z(),U2.Z(),O1O2.Z()) / D;
213 Param1 = Det33(O1O2.X(),U2.X(),N.X(),
214 O1O2.Y(),U2.Y(),N.Y(),
215 O1O2.Z(),U2.Z(),N.Z()) / (-D);
216 //------------------------------------------------------------
217 //-- On resout P1 * Dir1 + P2 * Dir2 + d * N = O1O2
218 //-- soit : Segment perpendiculaire : O1+P1 D1
219 //-- O2-P2 D2
220 Param2 = Det33(U1.X(),O1O2.X(),N.X(),
221 U1.Y(),O1O2.Y(),N.Y(),
222 U1.Z(),O1O2.Z(),N.Z()) / (D);
223 }
224}
225//=======================================================================
226//function : DirToAx2
227//purpose : returns a gp_Ax2 where D is the main direction
228//=======================================================================
229gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D)
230{
231 Standard_Real x=D.X(); Standard_Real ax=Abs(x);
232 Standard_Real y=D.Y(); Standard_Real ay=Abs(y);
233 Standard_Real z=D.Z(); Standard_Real az=Abs(z);
234 if( (ax==0.0) || ((ax<ay) && (ax<az)) ) {
235 return(gp_Ax2(P,D,gp_Dir(gp_Vec(0.0,-z,y))));
236 }
237 else if( (ay==0.0) || ((ay<ax) && (ay<az)) ) {
238 return(gp_Ax2(P,D,gp_Dir(gp_Vec(-z,0.0,x))));
239 }
240 else {
241 return(gp_Ax2(P,D,gp_Dir(gp_Vec(-y,x,0.0))));
242 }
243}
244//=======================================================================
245//function : IntAna_QuadQuadGeo
246//purpose : Empty constructor
247//=======================================================================
248 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(void)
249 : done(Standard_False),
250 nbint(0),
251 typeres(IntAna_Empty),
252 pt1(0,0,0),
253 pt2(0,0,0),
254 param1(0),
255 param2(0),
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,
281 const gp_Pln& P2,
282 const Standard_Real TolAng,
283 const Standard_Real Tol)
284: done(Standard_False),
285 nbint(0),
286 typeres(IntAna_Empty),
287 pt1(0,0,0),
288 pt2(0,0,0),
289 param1(0),
290 param2(0),
291 param1bis(0),
292 param2bis(0),
293 myCommonGen(Standard_False),
294 myPChar(0,0,0)
295{
296 InitTolerances();
297 Perform(P1,P2,TolAng,Tol);
298}
299//=======================================================================
300//function : Perform
301//purpose :
302//=======================================================================
303 void IntAna_QuadQuadGeo::Perform (const gp_Pln& P1,
304 const gp_Pln& P2,
305 const Standard_Real TolAng,
306 const Standard_Real Tol)
307{
308 done=Standard_False;
309 //
310 param2bis=0.0;
311
312 Standard_Real A1 = 0., B1 = 0., C1 = 0., D1 = 0., A2 = 0., B2 = 0., C2 = 0., D2 = 0.;
313 P1.Coefficients(A1,B1,C1,D1);
314 P2.Coefficients(A2,B2,C2,D2);
315
316 gp_Vec vd(gp_Vec(A1,B1,C1).Crossed(gp_Vec(A2,B2,C2)));
317 Standard_Real dist1= A2*P1.Location().X() + B2*P1.Location().Y() + C2*P1.Location().Z() + D2;
318 Standard_Real dist2= A1*P2.Location().X() + B1*P2.Location().Y() + C1*P2.Location().Z() + D1;
319
320 if(vd.Magnitude() <=TolAng) {
321 // normalles are collinear - planes are same or parallel
322 typeres = (Abs(dist1) <= Tol && Abs(dist2) <= Tol) ? IntAna_Same : IntAna_Empty;
323 }
324 else {
325 Standard_Real denom=A1*A2 + B1*B2 + C1*C2;
326
327 Standard_Real denom2 = denom*denom;
328 Standard_Real ddenom = 1. - denom2;
329 //denom = ( Abs(ddenom) <= 1.e-9 ) ? 1.e-9 : ddenom;
330 denom = ( Abs(ddenom) <= 1.e-16 ) ? 1.e-16 : ddenom;
331
332 Standard_Real par1 = dist1/denom;
333 Standard_Real par2 = -dist2/denom;
334
335 gp_Vec inter1(gp_Vec(A1,B1,C1).Crossed(vd));
336 gp_Vec inter2(gp_Vec(A2,B2,C2).Crossed(vd));
337
338 Standard_Real X1=P1.Location().X() + par1*inter1.X();
339 Standard_Real Y1=P1.Location().Y() + par1*inter1.Y();
340 Standard_Real Z1=P1.Location().Z() + par1*inter1.Z();
341 Standard_Real X2=P2.Location().X() + par2*inter2.X();
342 Standard_Real Y2=P2.Location().Y() + par2*inter2.Y();
343 Standard_Real Z2=P2.Location().Z() + par2*inter2.Z();
344
345 pt1=gp_Pnt((X1+X2)*0.5, (Y1+Y2)*0.5, (Z1+Z2)*0.5);
346 dir1 = gp_Dir(vd);
347 typeres = IntAna_Line;
348 nbint = 1;
349
350 }
351 done=Standard_True;
352}
353//=======================================================================
354//function : IntAna_QuadQuadGeo
355//purpose : Pln Cylinder
356//=======================================================================
357 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Pln& P
358 ,const gp_Cylinder& Cl
359 ,const Standard_Real Tolang
360 ,const Standard_Real Tol
361 ,const Standard_Real H)
362 : done(Standard_False),
363 nbint(0),
364 typeres(IntAna_Empty),
365 pt1(0,0,0),
366 pt2(0,0,0),
367 param1(0),
368 param2(0),
369 param1bis(0),
370 param2bis(0),
371 myCommonGen(Standard_False),
372 myPChar(0,0,0)
373{
374 InitTolerances();
375 Perform(P,Cl,Tolang,Tol,H);
376}
377//=======================================================================
378//function : Perform
379//purpose :
380//=======================================================================
381 void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
382 ,const gp_Cylinder& Cl
383 ,const Standard_Real Tolang
384 ,const Standard_Real Tol
385 ,const Standard_Real H)
386{
387 done = Standard_False;
388 Standard_Real dist,radius;
389 Standard_Real A,B,C,D;
390 Standard_Real X,Y,Z;
391 Standard_Real sint,cost,h;
392 gp_XYZ axex,axey,omega;
393
394
395 param2bis=0.0;
396 radius = Cl.Radius();
397
398 gp_Lin axec(Cl.Axis());
399 gp_XYZ normp(P.Axis().Direction().XYZ());
400
401 P.Coefficients(A,B,C,D);
402 axec.Location().Coord(X,Y,Z);
403 dist = A*X + B*Y + C*Z + D; // la distance axe/plan est evaluee a l origine de l axe.
404
405 Standard_Real tolang = Tolang;
406 Standard_Boolean newparams = Standard_False;
407
408 gp_Vec ldv( axec.Direction() );
409 gp_Vec npv( normp );
410 Standard_Real dA = Abs( ldv.Angle( npv ) );
411 if( dA > (M_PI/4.) )
412 {
413 Standard_Real dang = Abs( ldv.Angle( npv ) ) - M_PI/2.;
414 Standard_Real dangle = Abs( dang );
415 if( dangle > Tolang )
416 {
417 Standard_Real sinda = Abs( Sin( dangle ) );
418 Standard_Real dif = Abs( sinda - Tol );
419 if( dif < Tol )
420 {
421 tolang = sinda * 2.;
422 newparams = Standard_True;
423 }
424 }
425 }
426
427 nbint = 0;
428 IntAna_IntConicQuad inter(axec,P,tolang,Tol,H);
429
430 if (inter.IsParallel()) {
431 // Le resultat de l intersection Plan-Cylindre est de type droite.
432 // il y a 1 ou 2 droites
433
434 typeres = IntAna_Line;
435 omega.SetCoord(X-dist*A,Y-dist*B,Z-dist*C);
436
437 if (Abs(Abs(dist)-radius) < Tol)
438 {
439 nbint = 1;
440 pt1.SetXYZ(omega);
441
442 if( newparams )
443 {
444 gp_XYZ omegaXYZ(X,Y,Z);
445 gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
446 Standard_Real Xt,Yt,Zt,distt;
447 omegaXYZtrnsl.Coord(Xt,Yt,Zt);
448 distt = A*Xt + B*Yt + C*Zt + D;
449 gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A, omegaXYZtrnsl.Y()-distt*B, omegaXYZtrnsl.Z()-distt*C );
450 gp_Pnt ppt1;
451 ppt1.SetXYZ( omega1 );
452 gp_Vec vv1(pt1,ppt1);
453 gp_Dir dd1( vv1 );
454 dir1 = dd1;
455 }
456 else
457 dir1 = axec.Direction();
458 }
459 else if (Abs(dist) < radius)
460 {
461 nbint = 2;
462 h = Sqrt(radius*radius - dist*dist);
463 axey = axec.Direction().XYZ().Crossed(normp); // axey est normalise
464
465 pt1.SetXYZ(omega - h*axey);
466 pt2.SetXYZ(omega + h*axey);
467
468 if( newparams )
469 {
470 gp_XYZ omegaXYZ(X,Y,Z);
471 gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
472 Standard_Real Xt,Yt,Zt,distt,ht;
473 omegaXYZtrnsl.Coord(Xt,Yt,Zt);
474 distt = A*Xt + B*Yt + C*Zt + D;
475 // ht = Sqrt(radius*radius - distt*distt);
476 Standard_Real anSqrtArg = radius*radius - distt*distt;
477 ht = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
478
479 gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A, omegaXYZtrnsl.Y()-distt*B, omegaXYZtrnsl.Z()-distt*C );
480 gp_Pnt ppt1,ppt2;
481 ppt1.SetXYZ( omega1 - ht*axey);
482 ppt2.SetXYZ( omega1 + ht*axey);
483 gp_Vec vv1(pt1,ppt1);
484 gp_Vec vv2(pt2,ppt2);
485 gp_Dir dd1( vv1 );
486 gp_Dir dd2( vv2 );
487 dir1 = dd1;
488 dir2 = dd2;
489 }
490 else
491 {
492 dir1 = axec.Direction();
493 dir2 = axec.Direction();
494 }
495 }
496 // else nbint = 0
497
498 // debug JAG : le nbint = 0 doit etre remplace par typeres = IntAna_Empty
499 // et ne pas etre seulement supprime...
500
501 else {
502 typeres = IntAna_Empty;
503 }
504 }
505 else { // Il y a un point d intersection. C est le centre du cercle
506 // ou de l ellipse solution.
507
508 nbint = 1;
509 axey = normp.Crossed(axec.Direction().XYZ());
510 sint = axey.Modulus();
511
512 pt1 = inter.Point(1);
513
514 if (sint < Tol/radius) {
515
516 // on construit un cercle avec comme axes X et Y ceux du cylindre
517 typeres = IntAna_Circle;
518
519 dir1 = axec.Direction(); // axe Z
520 dir2 = Cl.Position().XDirection();
521 param1 = radius;
522 }
523 else {
524
525 // on construit un ellipse
526 typeres = IntAna_Ellipse;
527 cost = Abs(axec.Direction().XYZ().Dot(normp));
528 axex = axey.Crossed(normp);
529
530 dir1.SetXYZ(normp); //Modif ds ce bloc
531 dir2.SetXYZ(axex);
532
533 param1 = radius/cost;
534 param1bis = radius;
535 }
536 }
537
538 done = Standard_True;
539}
540//=======================================================================
541//function : IntAna_QuadQuadGeo
542//purpose : Pln Cone
543//=======================================================================
544 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
545 const gp_Cone& Co,
546 const Standard_Real Tolang,
547 const Standard_Real Tol)
548: done(Standard_False),
549 nbint(0),
550 typeres(IntAna_Empty),
551 pt1(0,0,0),
552 pt2(0,0,0),
553 param1(0),
554 param2(0),
555 param1bis(0),
556 param2bis(0),
557 myCommonGen(Standard_False),
558 myPChar(0,0,0)
559{
560 InitTolerances();
561 Perform(P,Co,Tolang,Tol);
562}
563//=======================================================================
564//function : Perform
565//purpose :
566//=======================================================================
567 void IntAna_QuadQuadGeo::Perform(const gp_Pln& P,
568 const gp_Cone& Co,
569 const Standard_Real Tolang,
570 const Standard_Real Tol)
571{
572
573 done = Standard_False;
574 nbint = 0;
575
576 Standard_Real A,B,C,D;
577 Standard_Real X,Y,Z;
578 Standard_Real dist,sint,cost,sina,cosa,angl,costa;
579 Standard_Real dh;
580 gp_XYZ axex,axey;
581
582 gp_Lin axec(Co.Axis());
583 P.Coefficients(A,B,C,D);
584 gp_Pnt apex(Co.Apex());
585
586 apex.Coord(X,Y,Z);
587 dist = A*X + B*Y + C*Z + D; // distance signee sommet du cone/ Plan
588
589 gp_XYZ normp = P.Axis().Direction().XYZ();
590 if(P.Direct()==Standard_False) { //-- lbr le 14 jan 97
591 normp.Reverse();
592 }
593
594 axey = normp.Crossed(Co.Axis().Direction().XYZ());
595 axex = axey.Crossed(normp);
596
597
598 angl = Co.SemiAngle();
599
600 cosa = Cos(angl);
601 sina = Abs(Sin(angl));
602
603
604 // Angle entre la normale au plan et l axe du cone, ramene entre 0. et PI/2.
605
606 sint = axey.Modulus();
607 cost = Abs(Co.Axis().Direction().XYZ().Dot(normp));
608
609 // Le calcul de costa permet de determiner si le plan contient
610 // un generatrice du cone : on calcul Sin((PI/2. - t) - angl)
611
612 costa = cost*cosa - sint*sina; // sin((PI/2 -t)-angl)=cos(t+angl)
613 // avec t ramene entre 0 et pi/2.
614
615 if (Abs(dist) < Tol) {
616 // on considere que le plan contient le sommet du cone.
617 // les solutions possibles sont donc : 1 point, 1 droite, 2 droites
618 // selon l inclinaison du plan.
619
620 if (Abs(costa) < Tolang) { // plan parallele a la generatrice
621 typeres = IntAna_Line;
622 nbint = 1;
623 gp_XYZ ptonaxe(apex.XYZ() + 10.*(Co.Axis().Direction().XYZ()));
624 // point sur l axe du cone cote z positif
625
626 dist = A*ptonaxe.X() + B*ptonaxe.Y() + C*ptonaxe.Z() + D;
627 ptonaxe = ptonaxe - dist*normp;
628 pt1 = apex;
629 dir1.SetXYZ(ptonaxe - pt1.XYZ());
630 }
631 else if (cost < sina) { // plan "interieur" au cone
632 typeres = IntAna_Line;
633 nbint = 2;
634 pt1 = apex;
635 pt2 = apex;
636 dh = Sqrt(sina*sina-cost*cost)/cosa;
637 dir1.SetXYZ(axex + dh*axey);
638 dir2.SetXYZ(axex - dh*axey);
639 }
640 else { // plan "exterieur" au cone
641 typeres = IntAna_Point;
642 nbint = 1;
643 pt1 = apex;
644 }
645 }
646 else {
647 // Solutions possibles : cercle, ellipse, parabole, hyperbole selon
648 // l inclinaison du plan.
649 Standard_Real deltacenter, distance;
650
651 if (cost < Tolang) {
652 // Le plan contient la direction de l axe du cone. La solution est
653 // l hyperbole
654 typeres = IntAna_Hyperbola;
655 nbint = 2;
656 pt1.SetXYZ(apex.XYZ()-dist*normp);
657 pt2 = pt1;
658 dir1=normp;
659 dir2.SetXYZ(axex);
660 param1 = param2 = Abs(dist/Tan(angl));
661 param1bis = param2bis = Abs(dist);
662 }
663 else {
664
665 IntAna_IntConicQuad inter(axec,P,Tolang); // on a necessairement 1 point.
666
667 gp_Pnt center(inter.Point(1));
668
669 // En fonction de la position de l intersection par rapport au sommet
670 // du cone, on change l axe x en -x et y en -y. Le parametre du sommet
671 // sur axec est negatif (voir definition du cone)
672
673 distance = apex.Distance(center);
674
675 if (inter.ParamOnConic(1) + Co.RefRadius()/Tan(angl) < 0.) {
676 axex.Reverse();
677 axey.Reverse();
678 }
679
680 if (Abs(costa) < Tolang) { // plan parallele a une generatrice
681 typeres = IntAna_Parabola;
682 nbint = 1;
683 deltacenter = distance/2./cosa;
684 axex.Normalize();
685 pt1.SetXYZ(center.XYZ()-deltacenter*axex);
686 dir1 = normp;
687 dir2.SetXYZ(axex);
688 param1 = deltacenter*sina*sina;
689 }
690 else if (sint < Tolang) { // plan perpendiculaire a l axe
691 typeres = IntAna_Circle;
692 nbint = 1;
693 pt1 = center;
694 dir1 = Co.Position().Direction();
695 dir2 = Co.Position().XDirection();
696 param1 = apex.Distance(center)*Abs(Tan(angl));
697 }
698 else if (cost < sina ) {
699 typeres = IntAna_Hyperbola;
700 nbint = 2;
701 axex.Normalize();
702
703 deltacenter = sint*sina*sina*distance/(sina*sina - cost*cost);
704 pt1.SetXYZ(center.XYZ() - deltacenter*axex);
705 pt2 = pt1;
706 dir1 = normp;
707 dir2.SetXYZ(axex);
708 param1 = param2 = cost*sina*cosa*distance /(sina*sina-cost*cost);
709 param1bis = param2bis = cost*sina*distance / Sqrt(sina*sina-cost*cost);
710
711 }
712 else { // on a alors cost > sina
713 typeres = IntAna_Ellipse;
714 nbint = 1;
715 Standard_Real radius = cost*sina*cosa*distance/(cost*cost-sina*sina);
716 deltacenter = sint*sina*sina*distance/(cost*cost-sina*sina);
717 axex.Normalize();
718 pt1.SetXYZ(center.XYZ() + deltacenter*axex);
719 dir1 = normp;
720 dir2.SetXYZ(axex);
721 param1 = radius;
722 param1bis = cost*sina*distance/ Sqrt(cost*cost - sina*sina);
723 }
724 }
725 }
726
727 //-- On a du mal a gerer plus loin (Value ProjLib, Params ... )
728 //-- des hyperboles trop bizarres
729 //-- On retourne False -> Traitement par biparametree
730 static Standard_Real EllipseLimit = 1.0E+9; //OCC513(apo) 1000000
731 static Standard_Real HyperbolaLimit = 2.0E+6; //OCC537(apo) 50000
732 if(typeres==IntAna_Ellipse && nbint>=1) {
733 if(Abs(param1) > EllipseLimit || Abs(param1bis) > EllipseLimit) {
734 done=Standard_False;
735 return;
736 }
737 }
738 if(typeres==IntAna_Hyperbola && nbint>=2) {
739 if(Abs(param2) > HyperbolaLimit || Abs(param2bis) > HyperbolaLimit) {
740 done = Standard_False;
741 return;
742 }
743 }
744 if(typeres==IntAna_Hyperbola && nbint>=1) {
745 if(Abs(param1) > HyperbolaLimit || Abs(param1bis) > HyperbolaLimit) {
746 done=Standard_False;
747 return;
748 }
749 }
750
751 done = Standard_True;
752}
753
754//=======================================================================
755//function : IntAna_QuadQuadGeo
756//purpose : Pln Sphere
757//=======================================================================
758 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
759 const gp_Sphere& S)
760: done(Standard_False),
761 nbint(0),
762 typeres(IntAna_Empty),
763 pt1(0,0,0),
764 pt2(0,0,0),
765 param1(0),
766 param2(0),
767 param1bis(0),
768 param2bis(0),
769 myCommonGen(Standard_False),
770 myPChar(0,0,0)
771{
772 InitTolerances();
773 Perform(P,S);
774}
775//=======================================================================
776//function : Perform
777//purpose :
778//=======================================================================
779 void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
780 ,const gp_Sphere& S)
781{
782
783 done = Standard_False;
784 Standard_Real A,B,C,D,dist, radius;
785 Standard_Real X,Y,Z;
786
787 nbint = 0;
788// debug JAG : on met typeres = IntAna_Empty par defaut...
789 typeres = IntAna_Empty;
790
791 P.Coefficients(A,B,C,D);
792 S.Location().Coord(X,Y,Z);
793 radius = S.Radius();
794
795 dist = A * X + B * Y + C * Z + D;
796
797 if (Abs( Abs(dist) - radius) < Epsilon(radius)) {
798 // on a une seule solution : le point projection du centre de la sphere
799 // sur le plan
800 nbint = 1;
801 typeres = IntAna_Point;
802 pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
803 }
804 else if (Abs(dist) < radius) {
805 // on a un cercle solution
806 nbint = 1;
807 typeres = IntAna_Circle;
808 pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
809 dir1 = P.Axis().Direction();
810 if(P.Direct()==Standard_False) dir1.Reverse();
811 dir2 = P.Position().XDirection();
812 param1 = Sqrt(radius*radius - dist*dist);
813 }
814 param2bis=0.0; //-- pour eviter param2bis not used ....
815 done = Standard_True;
816}
817
818//=======================================================================
819//function : IntAna_QuadQuadGeo
820//purpose : Cylinder - Cylinder
821//=======================================================================
822 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl1,
823 const gp_Cylinder& Cyl2,
824 const Standard_Real Tol)
825: done(Standard_False),
826 nbint(0),
827 typeres(IntAna_Empty),
828 pt1(0,0,0),
829 pt2(0,0,0),
830 param1(0),
831 param2(0),
832 param1bis(0),
833 param2bis(0),
834 myCommonGen(Standard_False),
835 myPChar(0,0,0)
836{
837 InitTolerances();
838 Perform(Cyl1,Cyl2,Tol);
839}
840//=======================================================================
841//function : Perform
842//purpose :
843//=======================================================================
844 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
845 const gp_Cylinder& Cyl2,
846 const Standard_Real Tol)
847{
848 done=Standard_True;
849 //---------------------------- Parallel axes -------------------------
850 AxeOperator A1A2(Cyl1.Axis(),Cyl2.Axis());
851 Standard_Real R1=Cyl1.Radius();
852 Standard_Real R2=Cyl2.Radius();
853 Standard_Real RmR, RmR_Relative;
854 RmR=(R1>R2)? (R1-R2) : (R2-R1);
855 {
856 Standard_Real Rmax, Rmin;
857 Rmax=(R1>R2)? R1 : R2;
858 Rmin=(R1>R2)? R2 : R1;
859 RmR_Relative=RmR/Rmax;
860 }
861
862 Standard_Real DistA1A2=A1A2.Distance();
863
864 if(A1A2.Parallel()) {
865 if(DistA1A2<=Tol) {
866 if(RmR<=Tol) {
867 typeres=IntAna_Same;
868 }
869 else {
870 typeres=IntAna_Empty;
871 }
872 }
873 else { //-- DistA1A2 > Tol
874 gp_Pnt P1=Cyl1.Location();
875 gp_Pnt P2t=Cyl2.Location();
876 gp_Pnt P2;
877 //-- P2t is projected on the plane (P1,DirCylX,DirCylY)
878 gp_Dir DirCyl = Cyl1.Position().Direction();
879 Standard_Real ProjP2OnDirCyl1=gp_Vec(DirCyl).Dot(gp_Vec(P1,P2t));
880
881 P2.SetCoord( P2t.X() - ProjP2OnDirCyl1*DirCyl.X()
882 ,P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y()
883 ,P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
884 //--
885 Standard_Real R1pR2=R1+R2;
886 if(DistA1A2>(R1pR2+Tol)) {
887 typeres=IntAna_Empty;
888 nbint=0;
889 }
890 else if(DistA1A2>(R1pR2)) {
891 //-- 1 Tangent line -------------------------------------OK
892 typeres=IntAna_Line;
893
894 nbint=1;
895 dir1=DirCyl;
896 Standard_Real R1_R1pR2=R1/R1pR2;
897 pt1.SetCoord( P1.X() + R1_R1pR2 * (P2.X()-P1.X())
898 ,P1.Y() + R1_R1pR2 * (P2.Y()-P1.Y())
899 ,P1.Z() + R1_R1pR2 * (P2.Z()-P1.Z()));
900
901 }
902 else if(DistA1A2>RmR) {
903 //-- 2 lines ---------------------------------------------OK
904 typeres=IntAna_Line;
905 nbint=2;
906 dir1=DirCyl;
907 gp_Vec P1P2(P1,P2);
908 gp_Dir DirA1A2=gp_Dir(P1P2);
909 gp_Dir Ortho_dir1_P1P2 = dir1.Crossed(DirA1A2);
910 dir2=dir1;
911 Standard_Real Alpha=0.5*(R1*R1-R2*R2+DistA1A2*DistA1A2)/(DistA1A2);
912
913// Standard_Real Beta = Sqrt(R1*R1-Alpha*Alpha);
914 Standard_Real anSqrtArg = R1*R1-Alpha*Alpha;
915 Standard_Real Beta = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
916
917 if((Beta+Beta)<Tol) {
918 nbint=1;
919 pt1.SetCoord( P1.X() + Alpha*DirA1A2.X()
920 ,P1.Y() + Alpha*DirA1A2.Y()
921 ,P1.Z() + Alpha*DirA1A2.Z());
922 }
923 else {
924 pt1.SetCoord( P1.X() + Alpha*DirA1A2.X() + Beta*Ortho_dir1_P1P2.X()
925 ,P1.Y() + Alpha*DirA1A2.Y() + Beta*Ortho_dir1_P1P2.Y()
926 ,P1.Z() + Alpha*DirA1A2.Z() + Beta*Ortho_dir1_P1P2.Z() );
927
928 pt2.SetCoord( P1.X() + Alpha*DirA1A2.X() - Beta*Ortho_dir1_P1P2.X()
929 ,P1.Y() + Alpha*DirA1A2.Y() - Beta*Ortho_dir1_P1P2.Y()
930 ,P1.Z() + Alpha*DirA1A2.Z() - Beta*Ortho_dir1_P1P2.Z());
931 }
932 }
933 else if(DistA1A2>(RmR-Tol)) {
934 //-- 1 Tangent ------------------------------------------OK
935 typeres=IntAna_Line;
936 nbint=1;
937 dir1=DirCyl;
938 Standard_Real R1_RmR=R1/RmR;
939
940 if(R1 < R2) R1_RmR = -R1_RmR;
941
942 pt1.SetCoord( P1.X() + R1_RmR * (P2.X()-P1.X())
943 ,P1.Y() + R1_RmR * (P2.Y()-P1.Y())
944 ,P1.Z() + R1_RmR * (P2.Z()-P1.Z()));
945 }
946 else {
947 nbint=0;
948 typeres=IntAna_Empty;
949 }
950 }
951 }
952 else { //-- No Parallel Axis ---------------------------------OK
953 if((RmR_Relative<=myEPSILON_CYLINDER_DELTA_RADIUS)
954 && (DistA1A2 <= myEPSILON_CYLINDER_DELTA_DISTANCE)) {
955 //-- PI/2 between the two axis and Intersection
956 //-- and identical radius
957 typeres=IntAna_Ellipse;
958 nbint=2;
959 gp_Dir DirCyl1=Cyl1.Position().Direction();
960 gp_Dir DirCyl2=Cyl2.Position().Direction();
961 pt1=pt2=A1A2.PtIntersect();
962
963 Standard_Real A=DirCyl1.Angle(DirCyl2);
964 Standard_Real B;
965 B=Abs(Sin(0.5*(M_PI-A)));
966 A=Abs(Sin(0.5*A));
967
968 if(A==0.0 || B==0.0) {
969 typeres=IntAna_Same;
970 return;
971 }
972
973
974 gp_Vec dircyl1(DirCyl1);gp_Vec dircyl2(DirCyl2);
975 dir1 = gp_Dir(dircyl1.Added(dircyl2));
976 dir2 = gp_Dir(dircyl1.Subtracted(dircyl2));
977
978 param2 = Cyl1.Radius() / A;
979 param1 = Cyl1.Radius() / B;
980 param2bis= param1bis = Cyl1.Radius();
981 if(param1 < param1bis) {
982 A=param1; param1=param1bis; param1bis=A;
983 }
984 if(param2 < param2bis) {
985 A=param2; param2=param2bis; param2bis=A;
986 }
987 }
988 else {
989 if(Abs(DistA1A2-Cyl1.Radius()-Cyl2.Radius())<Tol) {
990 typeres = IntAna_Point;
991 Standard_Real d,p1,p2;
992
993 gp_Dir D1 = Cyl1.Axis().Direction();
994 gp_Dir D2 = Cyl2.Axis().Direction();
995 A1A2.Distance(d,p1,p2);
996 gp_Pnt P = Cyl1.Axis().Location();
997 gp_Pnt P1(P.X() - p1*D1.X(),
998 P.Y() - p1*D1.Y(),
999 P.Z() - p1*D1.Z());
1000 P = Cyl2.Axis().Location();
1001 gp_Pnt P2(P.X() - p2*D2.X(),
1002 P.Y() - p2*D2.Y(),
1003 P.Z() - p2*D2.Z());
1004 gp_Vec P1P2(P1,P2);
1005 D1=gp_Dir(P1P2);
1006 p1=Cyl1.Radius();
1007 pt1.SetCoord(P1.X() + p1*D1.X(),
1008 P1.Y() + p1*D1.Y(),
1009 P1.Z() + p1*D1.Z());
1010 nbint = 1;
1011 }
1012 else {
1013 typeres=IntAna_NoGeometricSolution;
1014 }
1015 }
1016 }
1017}
1018//=======================================================================
1019//function : IntAna_QuadQuadGeo
1020//purpose : Cylinder - Cone
1021//=======================================================================
1022 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1023 const gp_Cone& Con,
1024 const Standard_Real Tol)
1025: done(Standard_False),
1026 nbint(0),
1027 typeres(IntAna_Empty),
1028 pt1(0,0,0),
1029 pt2(0,0,0),
1030 param1(0),
1031 param2(0),
1032 param1bis(0),
1033 param2bis(0),
1034 myCommonGen(Standard_False),
1035 myPChar(0,0,0)
1036{
1037 InitTolerances();
1038 Perform(Cyl,Con,Tol);
1039}
1040//=======================================================================
1041//function : Perform
1042//purpose :
1043//=======================================================================
1044 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
1045 const gp_Cone& Con,
1046 const Standard_Real )
1047{
1048 done=Standard_True;
1049 AxeOperator A1A2(Cyl.Axis(),Con.Axis());
1050 if(A1A2.Same()) {
1051 gp_Pnt Pt=Con.Apex();
1052 Standard_Real dist=Cyl.Radius()/(Tan(Con.SemiAngle()));
1053 gp_Dir dir=Cyl.Position().Direction();
1054 pt1.SetCoord( Pt.X() + dist*dir.X()
1055 ,Pt.Y() + dist*dir.Y()
1056 ,Pt.Z() + dist*dir.Z());
1057 pt2.SetCoord( Pt.X() - dist*dir.X()
1058 ,Pt.Y() - dist*dir.Y()
1059 ,Pt.Z() - dist*dir.Z());
1060 dir1=dir2=dir;
1061 param1=param2=Cyl.Radius();
1062 nbint=2;
1063 typeres=IntAna_Circle;
1064
1065 }
1066 else {
1067 typeres=IntAna_NoGeometricSolution;
1068 }
1069}
1070//=======================================================================
1071//function :
1072//purpose : Cylinder - Sphere
1073//=======================================================================
1074 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1075 const gp_Sphere& Sph,
1076 const Standard_Real Tol)
1077: done(Standard_False),
1078 nbint(0),
1079 typeres(IntAna_Empty),
1080 pt1(0,0,0),
1081 pt2(0,0,0),
1082 param1(0),
1083 param2(0),
1084 param1bis(0),
1085 param2bis(0),
1086 myCommonGen(Standard_False),
1087 myPChar(0,0,0)
1088{
1089 InitTolerances();
1090 Perform(Cyl,Sph,Tol);
1091}
1092//=======================================================================
1093//function : Perform
1094//purpose :
1095//=======================================================================
1096 void IntAna_QuadQuadGeo::Perform( const gp_Cylinder& Cyl
1097 ,const gp_Sphere& Sph
1098 ,const Standard_Real)
1099{
1100 done=Standard_True;
1101 gp_Pnt Pt=Sph.Location();
1102 AxeOperator A1A2(Cyl.Axis(),Sph.Position().Axis());
1103 if((A1A2.Intersect() && Pt.Distance(A1A2.PtIntersect())==0.0 )
1104 || (A1A2.Same())) {
1105 if(Sph.Radius() < Cyl.Radius()) {
1106 typeres = IntAna_Empty;
1107 }
1108 else {
1109 Standard_Real dist=Sqrt( Sph.Radius() * Sph.Radius() - Cyl.Radius() * Cyl.Radius() );
1110 gp_Dir dir=Cyl.Position().Direction();
1111 dir1 = dir2 = dir;
1112 typeres=IntAna_Circle;
1113 pt1.SetCoord( Pt.X() + dist*dir.X()
1114 ,Pt.Y() + dist*dir.Y()
1115 ,Pt.Z() + dist*dir.Z());
1116 nbint=1;
1117 param1 = Cyl.Radius();
1118 if(dist>RealEpsilon()) {
1119 pt2.SetCoord( Pt.X() - dist*dir.X()
1120 ,Pt.Y() - dist*dir.Y()
1121 ,Pt.Z() - dist*dir.Z());
1122 param2=Cyl.Radius();
1123 nbint=2;
1124 }
1125 }
1126 }
1127 else {
1128 typeres=IntAna_NoGeometricSolution;
1129 }
1130}
1131
1132//=======================================================================
1133//function : IntAna_QuadQuadGeo
1134//purpose : Cone - Cone
1135//=======================================================================
1136 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con1,
1137 const gp_Cone& Con2,
1138 const Standard_Real Tol)
1139: done(Standard_False),
1140 nbint(0),
1141 typeres(IntAna_Empty),
1142 pt1(0,0,0),
1143 pt2(0,0,0),
1144 param1(0),
1145 param2(0),
1146 param1bis(0),
1147 param2bis(0),
1148 myCommonGen(Standard_False),
1149 myPChar(0,0,0)
1150{
1151 InitTolerances();
1152 Perform(Con1,Con2,Tol);
1153}
1154//
1155//=======================================================================
1156//function : Perform
1157//purpose :
1158//=======================================================================
1159 void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con1,
1160 const gp_Cone& Con2,
1161 const Standard_Real Tol)
1162{
1163 done=Standard_True;
1164 //
1165 Standard_Real tg1, tg2, aDA1A2, aTol2;
1166 gp_Pnt aPApex1, aPApex2;
1167
1168 Standard_Real TOL_APEX_CONF = 1.e-10;
1169
1170 //
1171 tg1=Tan(Con1.SemiAngle());
1172 tg2=Tan(Con2.SemiAngle());
1173
1174 if((tg1 * tg2) < 0.) {
1175 tg2 = -tg2;
1176 }
1177 //
1178 aTol2=Tol*Tol;
1179 aPApex1=Con1.Apex();
1180 aPApex2=Con2.Apex();
1181 aDA1A2=aPApex1.SquareDistance(aPApex2);
1182 //
1183 AxeOperator A1A2(Con1.Axis(),Con2.Axis());
1184 //
1185 // 1
1186 if(A1A2.Same()) {
1187 //-- two circles
1188 Standard_Real x;
1189 gp_Pnt P=Con1.Apex();
1190 gp_Dir D=Con1.Position().Direction();
1191 Standard_Real d=gp_Vec(D).Dot(gp_Vec(P,Con2.Apex()));
1192
1193 if(Abs(tg1-tg2)>myEPSILON_ANGLE_CONE) {
1194 if (fabs(d) < TOL_APEX_CONF) {
1195 typeres = IntAna_Point;
1196 nbint = 1;
1197 pt1 = P;
1198 return;
1199 }
1200 x=(d*tg2)/(tg1+tg2);
1201 pt1.SetCoord( P.X() + x*D.X()
1202 ,P.Y() + x*D.Y()
1203 ,P.Z() + x*D.Z());
1204 param1=Abs(x*tg1);
1205
1206 x=(d*tg2)/(tg2-tg1);
1207 pt2.SetCoord( P.X() + x*D.X()
1208 ,P.Y() + x*D.Y()
1209 ,P.Z() + x*D.Z());
1210 param2=Abs(x*tg1);
1211 dir1 = dir2 = D;
1212 nbint=2;
1213 typeres=IntAna_Circle;
1214 }
1215 else {
1216 if (fabs(d) < TOL_APEX_CONF) {
1217 typeres=IntAna_Same;
1218 }
1219 else {
1220 typeres=IntAna_Circle;
1221 nbint=1;
1222 x=d*0.5;
1223 pt1.SetCoord( P.X() + x*D.X()
1224 ,P.Y() + x*D.Y()
1225 ,P.Z() + x*D.Z());
1226 param1 = Abs(x * tg1);
1227 dir1 = D;
1228 }
1229 }
1230 } //-- fin A1A2.Same
1231 // 2
1232 else if((Abs(tg1-tg2)<myEPSILON_ANGLE_CONE) && (A1A2.Parallel())) {
1233 //-- voir AnVer12mai98
1234 Standard_Real DistA1A2=A1A2.Distance();
1235 gp_Dir DA1=Con1.Position().Direction();
1236 gp_Vec O1O2(Con1.Apex(),Con2.Apex());
1237 gp_Dir O1O2n(O1O2); // normalization of the vector before projection
1238 Standard_Real O1O2_DA1=gp_Vec(DA1).Dot(gp_Vec(O1O2n));
1239
1240 gp_Vec O1_Proj_A2(O1O2n.X()-O1O2_DA1*DA1.X(),
1241 O1O2n.Y()-O1O2_DA1*DA1.Y(),
1242 O1O2n.Z()-O1O2_DA1*DA1.Z());
1243 gp_Dir DB1=gp_Dir(O1_Proj_A2);
1244
1245 Standard_Real yO1O2=O1O2.Dot(gp_Vec(DA1));
1246 Standard_Real ABSTG1 = Abs(tg1);
1247 Standard_Real X2 = (DistA1A2/ABSTG1 - yO1O2)*0.5;
1248 Standard_Real X1 = X2+yO1O2;
1249
1250 gp_Pnt P1(Con1.Apex().X() + X1*( DA1.X() + ABSTG1*DB1.X()),
1251 Con1.Apex().Y() + X1*( DA1.Y() + ABSTG1*DB1.Y()),
1252 Con1.Apex().Z() + X1*( DA1.Z() + ABSTG1*DB1.Z()));
1253
1254 gp_Pnt MO1O2(0.5*(Con1.Apex().X()+Con2.Apex().X()),
1255 0.5*(Con1.Apex().Y()+Con2.Apex().Y()),
1256 0.5*(Con1.Apex().Z()+Con2.Apex().Z()));
1257 gp_Vec P1MO1O2(P1,MO1O2);
1258
1259 gp_Dir DA1_X_DB1=DA1.Crossed(DB1);
1260 gp_Dir OrthoPln = DA1_X_DB1.Crossed(gp_Dir(P1MO1O2));
1261
1262 IntAna_QuadQuadGeo INTER_QUAD_PLN(gp_Pln(P1,OrthoPln),Con1,Tol,Tol);
1263 if(INTER_QUAD_PLN.IsDone()) {
1264 switch(INTER_QUAD_PLN.TypeInter()) {
1265 case IntAna_Ellipse: {
1266 typeres=IntAna_Ellipse;
1267 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1268 pt1 = E.Location();
1269 dir1 = E.Position().Direction();
1270 dir2 = E.Position().XDirection();
1271 param1 = E.MajorRadius();
1272 param1bis = E.MinorRadius();
1273 nbint = 1;
1274 break;
1275 }
1276 case IntAna_Circle: {
1277 typeres=IntAna_Circle;
1278 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1279 pt1 = C.Location();
1280 dir1 = C.Position().XDirection();
1281 dir2 = C.Position().YDirection();
1282 param1 = C.Radius();
1283 nbint = 1;
1284 break;
1285 }
1286 case IntAna_Hyperbola: {
1287 typeres=IntAna_Hyperbola;
1288 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1289 pt1 = pt2 = H.Location();
1290 dir1 = H.Position().Direction();
1291 dir2 = H.Position().XDirection();
1292 param1 = param2 = H.MajorRadius();
1293 param1bis = param2bis = H.MinorRadius();
1294 nbint = 2;
1295 break;
1296 }
1297 case IntAna_Line: {
1298 typeres=IntAna_Line;
1299 gp_Lin H=INTER_QUAD_PLN.Line(1);
1300 pt1 = pt2 = H.Location();
1301 dir1 = dir2 = H.Position().Direction();
1302 param1 = param2 = 0.0;
1303 param1bis = param2bis = 0.0;
1304 nbint = 2;
1305 break;
1306 }
1307 default:
1308 typeres=IntAna_NoGeometricSolution;
1309 }
1310 }
1311 }// else if((Abs(tg1-tg2)<EPSILON_ANGLE_CONE) && (A1A2.Parallel()))
1312 // 3
1313 else if (aDA1A2<aTol2) {
1314 //
1315 // When apices are coinsided there can be 3 possible cases
1316 // 3.1 - empty solution (iRet=0)
1317 // 3.2 - one line when cone1 touches cone2 (iRet=1)
1318 // 3.3 - two lines when cone1 intersects cone2 (iRet=2)
1319 //
1320 Standard_Integer iRet;
1321 Standard_Real aGamma, aBeta1, aBeta2;
1322 Standard_Real aD1, aR1, aTgBeta1, aTgBeta2, aHalfPI;
1323 Standard_Real aCosGamma, aSinGamma, aDx, aR2, aRD2, aD2;
1324 gp_Pnt2d aP0, aPA1, aP1, aPA2;
1325 gp_Vec2d aVAx2;
1326 gp_Ax1 aAx1, aAx2;
1327 //
1328 // Preliminary analysis. Determination of iRet
1329 //
1330 iRet=0;
1331 aHalfPI=0.5*M_PI;
1332 aD1=1.;
1333 aPA1.SetCoord(aD1, 0.);
1334 aP0.SetCoord(0., 0.);
1335 //
1336 aAx1=Con1.Axis();
1337 aAx2=Con2.Axis();
1338 aGamma=aAx1.Angle(aAx2);
1339 if (aGamma>aHalfPI){
1340 aGamma=M_PI-aGamma;
1341 }
1342 aCosGamma=Cos(aGamma);
1343 aSinGamma=Sin(aGamma);
1344 //
1345 aBeta1=Con1.SemiAngle();
1346 aTgBeta1=Tan(aBeta1);
1347 aTgBeta1=Abs(aTgBeta1);
1348 //
1349 aBeta2=Con2.SemiAngle();
1350 aTgBeta2=Tan(aBeta2);
1351 aTgBeta2=Abs(aTgBeta2);
1352 //
1353 aR1=aD1*aTgBeta1;
1354 aP1.SetCoord(aD1, aR1);
1355 //
1356 // PA2
1357 aVAx2.SetCoord(aCosGamma, aSinGamma);
1358 gp_Dir2d aDAx2(aVAx2);
1359 gp_Lin2d aLAx2(aP0, aDAx2);
1360 //
1361 gp_Vec2d aV(aP0, aP1);
1362 aDx=aVAx2.Dot(aV);
1363 aPA2=aP0.Translated(aDx*aDAx2);
1364 //
1365 // aR2
1366 aDx=aPA2.Distance(aP0);
1367 aR2=aDx*aTgBeta2;
1368 //
1369 // aRD2
1370 aRD2=aPA2.Distance(aP1);
1371 //
1372 if (aRD2>(aR2+Tol)) {
1373 iRet=0;
1374 typeres=IntAna_Empty; //nothing
1375 return;
1376 }
1377 //
1378 iRet=1; //touch case => 1 line
1379 if (aRD2<(aR2-Tol)) {
1380 iRet=2;//intersection => couple of lines
1381 }
1382 //
1383 // Finding the solution in 3D
1384 //
1385 Standard_Real aDa;
1386 gp_Pnt aQApex1, aQA1, aQA2, aQX, aQX1, aQX2;
1387 gp_Dir aD3Ax1, aD3Ax2;
1388 gp_Lin aLin;
1389 IntAna_QuadQuadGeo aIntr;
1390 //
1391 aQApex1=Con1.Apex();
1392 aD3Ax1=aAx1.Direction();
1393 aQA1.SetCoord(aQApex1.X()+aD1*aD3Ax1.X(),
1394 aQApex1.Y()+aD1*aD3Ax1.Y(),
1395 aQApex1.Z()+aD1*aD3Ax1.Z());
1396 //
1397 aDx=aD3Ax1.Dot(aAx2.Direction());
1398 if (aDx<0.) {
1399 aAx2.Reverse();
1400 }
1401 aD3Ax2=aAx2.Direction();
1402 //
1403 aD2=aD1*sqrt((1.+aTgBeta1*aTgBeta1)/(1.+aTgBeta2*aTgBeta2));
1404 //
1405 aQA2.SetCoord(aQApex1.X()+aD2*aD3Ax2.X(),
1406 aQApex1.Y()+aD2*aD3Ax2.Y(),
1407 aQApex1.Z()+aD2*aD3Ax2.Z());
1408 //
1409 gp_Pln aPln1(aQA1, aD3Ax1);
1410 gp_Pln aPln2(aQA2, aD3Ax2);
1411 //
1412 aIntr.Perform(aPln1, aPln2, Tol, Tol);
1413 if (!aIntr.IsDone()) {
1414 iRet=-1; // just in case. it must not be so
1415 typeres=IntAna_NoGeometricSolution;
1416 return;
1417 }
1418 //
1419 aLin=aIntr.Line(1);
1420 const gp_Dir& aDLin=aLin.Direction();
1421 gp_Vec aVLin(aDLin);
1422 gp_Pnt aOrig=aLin.Location();
1423 gp_Vec aVr(aQA1, aOrig);
1424 aDx=aVLin.Dot(aVr);
1425 aQX=aOrig.Translated(aDx*aVLin);
1426 //
1427 // Final part
1428 //
1429 typeres=IntAna_Line;
1430 //
1431 param1=0.;
1432 param2 =0.;
1433 param1bis=0.;
1434 param2bis=0.;
1435 //
1436 if (iRet==1) {
1437 // one line
1438 nbint=1;
1439 pt1=aQApex1;
1440 gp_Vec aVX(aQApex1, aQX);
1441 dir1=gp_Dir(aVX);
1442 }
1443
1444 else {//iRet=2
1445 // two lines
1446 nbint=2;
1447 aDa=aQA1.Distance(aQX);
1448 aDx=sqrt(aR1*aR1-aDa*aDa);
1449 aQX1=aQX.Translated(aDx*aVLin);
1450 aQX2=aQX.Translated(-aDx*aVLin);
1451 //
1452 pt1=aQApex1;
1453 pt2=aQApex1;
1454 gp_Vec aVX1(aQApex1, aQX1);
1455 dir1=gp_Dir(aVX1);
1456 gp_Vec aVX2(aQApex1, aQX2);
1457 dir2=gp_Dir(aVX2);
1458 }
1459 } //else if (aDA1A2<aTol2) {
1460 //Case when cones have common generatrix
1461 else if(A1A2.Intersect()) {
1462 //Check if apex of one cone belongs another one
1463 Standard_Real u, v, tol2 = Tol*Tol;
1464 ElSLib::Parameters(Con2, aPApex1, u, v);
1465 gp_Pnt p = ElSLib::Value(u, v, Con2);
1466 if(aPApex1.SquareDistance(p) > tol2) {
1467 typeres=IntAna_NoGeometricSolution;
1468 return;
1469 }
1470 //
1471 ElSLib::Parameters(Con1, aPApex2, u, v);
1472 p = ElSLib::Value(u, v, Con1);
1473 if(aPApex2.SquareDistance(p) > tol2) {
1474 typeres=IntAna_NoGeometricSolution;
1475 return;
1476 }
1477
1478 //Cones have a common generatrix passing through apexes
1479 myCommonGen = Standard_True;
1480
1481 //common generatrix of cones
1482 gp_Lin aGen(aPApex1, gp_Dir(gp_Vec(aPApex1, aPApex2)));
1483
1484 //Intersection point of axes
1485 gp_Pnt aPAxeInt = A1A2.PtIntersect();
1486
1487 //Characteristic point of intersection curve
1488 u = ElCLib::Parameter(aGen, aPAxeInt);
1489 myPChar = ElCLib::Value(u, aGen);
1490
1491
1492 //Other generatrixes of cones laying in maximal plane
1493 gp_Lin aGen1 = aGen.Rotated(Con1.Axis(), M_PI);
1494 gp_Lin aGen2 = aGen.Rotated(Con2.Axis(), M_PI);
1495 //
1496 //Intersection point of generatrixes
1497 gp_Dir aN; //solution plane normal
1498 gp_Dir aD1 = aGen1.Direction();
1499
1500 gp_Dir aD2(aD1.Crossed(aGen.Direction()));
1501
1502 if(aD1.IsParallel(aGen2.Direction(), Precision::Angular())) {
1503 aN = aD1.Crossed(aD2);
1504 }
1505 else if(aGen1.SquareDistance(aGen2) > tol2) {
1506 //Something wrong ???
1507 typeres=IntAna_NoGeometricSolution;
1508 return;
1509 }
1510 else {
1511 gp_Dir D1 = aGen1.Position().Direction();
1512 gp_Dir D2 = aGen2.Position().Direction();
1513 gp_Pnt O1 = aGen1.Location();
1514 gp_Pnt O2 = aGen2.Location();
1515 Standard_Real D1DotD2 = D1.Dot(D2);
1516 Standard_Real aSin = 1.-D1DotD2*D1DotD2;
1517 gp_Vec O1O2 (O1,O2);
1518 Standard_Real U2 = (D1.XYZ()*(O1O2.Dot(D1))-(O1O2.XYZ())).Dot(D2.XYZ());
1519 U2 /= aSin;
1520 gp_Pnt aPGint(ElCLib::Value(U2, aGen2));
1521
1522 aD1 = gp_Dir(gp_Vec(aPGint, myPChar));
1523 aN = aD1.Crossed(aD2);
1524 }
1525 //Plane that must contain intersection curves
1526 gp_Pln anIntPln(myPChar, aN);
1527
1528 IntAna_QuadQuadGeo INTER_QUAD_PLN(anIntPln,Con1,Tol,Tol);
1529
1530 if(INTER_QUAD_PLN.IsDone()) {
1531 switch(INTER_QUAD_PLN.TypeInter()) {
1532 case IntAna_Ellipse: {
1533 typeres=IntAna_Ellipse;
1534 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1535 pt1 = E.Location();
1536 dir1 = E.Position().Direction();
1537 dir2 = E.Position().XDirection();
1538 param1 = E.MajorRadius();
1539 param1bis = E.MinorRadius();
1540 nbint = 1;
1541 break;
1542 }
1543 case IntAna_Circle: {
1544 typeres=IntAna_Circle;
1545 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1546 pt1 = C.Location();
1547 dir1 = C.Position().XDirection();
1548 dir2 = C.Position().YDirection();
1549 param1 = C.Radius();
1550 nbint = 1;
1551 break;
1552 }
1553 case IntAna_Parabola: {
1554 typeres=IntAna_Parabola;
1555 gp_Parab Prb=INTER_QUAD_PLN.Parabola(1);
1556 pt1 = Prb.Location();
1557 dir1 = Prb.Position().Direction();
1558 dir2 = Prb.Position().XDirection();
1559 param1 = Prb.Focal();
1560 nbint = 1;
1561 break;
1562 }
1563 case IntAna_Hyperbola: {
1564 typeres=IntAna_Hyperbola;
1565 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1566 pt1 = pt2 = H.Location();
1567 dir1 = H.Position().Direction();
1568 dir2 = H.Position().XDirection();
1569 param1 = param2 = H.MajorRadius();
1570 param1bis = param2bis = H.MinorRadius();
1571 nbint = 2;
1572 break;
1573 }
1574 default:
1575 typeres=IntAna_NoGeometricSolution;
1576 }
1577 }
1578 }
1579
1580 else {
1581 typeres=IntAna_NoGeometricSolution;
1582 }
1583}
1584//=======================================================================
1585//function : IntAna_QuadQuadGeo
1586//purpose : Sphere - Cone
1587//=======================================================================
1588 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
1589 const gp_Cone& Con,
1590 const Standard_Real Tol)
1591: done(Standard_False),
1592 nbint(0),
1593 typeres(IntAna_Empty),
1594 pt1(0,0,0),
1595 pt2(0,0,0),
1596 param1(0),
1597 param2(0),
1598 param1bis(0),
1599 param2bis(0),
1600 myCommonGen(Standard_False),
1601 myPChar(0,0,0)
1602{
1603 InitTolerances();
1604 Perform(Sph,Con,Tol);
1605}
1606//=======================================================================
1607//function : Perform
1608//purpose :
1609//=======================================================================
1610 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
1611 const gp_Cone& Con,
1612 const Standard_Real)
1613{
1614
1615 //
1616 done=Standard_True;
1617 //
1618 AxeOperator A1A2(Con.Axis(),Sph.Position().Axis());
1619 gp_Pnt Pt=Sph.Location();
1620 //
1621 if((A1A2.Intersect() && (Pt.Distance(A1A2.PtIntersect())==0.0))
1622 || A1A2.Same()) {
1623 gp_Pnt ConApex= Con.Apex();
1624 Standard_Real dApexSphCenter=Pt.Distance(ConApex);
1625 gp_Dir ConDir;
1626 if(dApexSphCenter>RealEpsilon()) {
1627 ConDir = gp_Dir(gp_Vec(ConApex,Pt));
1628 }
1629 else {
1630 ConDir = Con.Position().Direction();
1631 }
1632
1633 Standard_Real Rad=Sph.Radius();
1634 Standard_Real tga=Tan(Con.SemiAngle());
1635
1636
1637 //-- 2 circles
1638 //-- x: Roots of (x**2 + y**2 = Rad**2)
1639 //-- tga = y / (x+dApexSphCenter)
1640 Standard_Real tgatga = tga * tga;
1641 math_DirectPolynomialRoots Eq( 1.0+tgatga
1642 ,2.0*tgatga*dApexSphCenter
1643 ,-Rad*Rad + dApexSphCenter*dApexSphCenter*tgatga);
1644 if(Eq.IsDone()) {
1645 Standard_Integer nbsol=Eq.NbSolutions();
1646 if(nbsol==0) {
1647 typeres=IntAna_Empty;
1648 }
1649 else {
1650 typeres=IntAna_Circle;
1651 if(nbsol>=1) {
1652 Standard_Real x = Eq.Value(1);
1653 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1654 nbint=1;
1655 pt1.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1656 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1657 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1658 param1 = tga * dApexSphCenterpx;
1659 param1 = Abs(param1);
1660 dir1 = ConDir;
1661 if(param1<=myEPSILON_MINI_CIRCLE_RADIUS) {
1662 typeres=IntAna_PointAndCircle;
1663 param1=0.0;
1664 }
1665 }
1666 if(nbsol>=2) {
1667 Standard_Real x=Eq.Value(2);
1668 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1669 nbint=2;
1670 pt2.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1671 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1672 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1673 param2 = tga * dApexSphCenterpx;
1674 param2 = Abs(param2);
1675 dir2=ConDir;
1676 if(param2<=myEPSILON_MINI_CIRCLE_RADIUS) {
1677 typeres=IntAna_PointAndCircle;
1678 param2=0.0;
1679 }
1680 }
1681 }
1682 }
1683 else {
1684 done=Standard_False;
1685 }
1686 }
1687 else {
1688 typeres=IntAna_NoGeometricSolution;
1689 }
1690}
1691
1692//=======================================================================
1693//function : IntAna_QuadQuadGeo
1694//purpose : Sphere - Sphere
1695//=======================================================================
1696 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Sphere& Sph1
1697 ,const gp_Sphere& Sph2
1698 ,const Standard_Real Tol)
1699: done(Standard_False),
1700 nbint(0),
1701 typeres(IntAna_Empty),
1702 pt1(0,0,0),
1703 pt2(0,0,0),
1704 param1(0),
1705 param2(0),
1706 param1bis(0),
1707 param2bis(0),
1708 myCommonGen(Standard_False),
1709 myPChar(0,0,0)
1710{
1711 InitTolerances();
1712 Perform(Sph1,Sph2,Tol);
1713}
1714//=======================================================================
1715//function : Perform
1716//purpose :
1717//=======================================================================
1718 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph1,
1719 const gp_Sphere& Sph2,
1720 const Standard_Real Tol)
1721{
1722 done=Standard_True;
1723 gp_Pnt O1=Sph1.Location();
1724 gp_Pnt O2=Sph2.Location();
1725 Standard_Real dO1O2=O1.Distance(O2);
1726 Standard_Real R1=Sph1.Radius();
1727 Standard_Real R2=Sph2.Radius();
1728 Standard_Real Rmin,Rmax;
1729 typeres=IntAna_Empty;
1730 param2bis=0.0; //-- pour eviter param2bis not used ....
1731
1732 if(R1>R2) { Rmin=R2; Rmax=R1; } else { Rmin=R1; Rmax=R2; }
1733
1734 if(dO1O2<=Tol && (Abs(R1-R2) <= Tol)) {
1735 typeres = IntAna_Same;
1736 }
1737 else {
1738 if(dO1O2<=Tol) { return; }
1739 gp_Dir Dir=gp_Dir(gp_Vec(O1,O2));
1740 Standard_Real t = Rmax - dO1O2 - Rmin;
1741
1742 //----------------------------------------------------------------------
1743 //-- |----------------- R1 --------------------|
1744 //-- |----dO1O2-----|-----------R2----------|
1745 //-- --->--<-- t
1746 //--
1747 //-- |------ R1 ------|---------dO1O2----------|
1748 //-- |-------------------R2-----------------------|
1749 //-- --->--<-- t
1750 //----------------------------------------------------------------------
1751 if(t >= 0.0 && t <=Tol) {
1752 typeres = IntAna_Point;
1753 nbint = 1;
1754 Standard_Real t2;
1755 if(R1==Rmax) t2=(R1 + (R2 + dO1O2)) * 0.5;
1756 else t2=(-R1+(dO1O2-R2))*0.5;
1757
1758 pt1.SetCoord( O1.X() + t2*Dir.X()
1759 ,O1.Y() + t2*Dir.Y()
1760 ,O1.Z() + t2*Dir.Z());
1761 }
1762 else {
1763 //-----------------------------------------------------------------
1764 //-- |----------------- dO1O2 --------------------|
1765 //-- |----R1-----|-----------R2----------|-Tol-|
1766 //--
1767 //-- |----------------- Rmax --------------------|
1768 //-- |----Rmin----|-------dO1O2-------|-Tol-|
1769 //--
1770 //-----------------------------------------------------------------
1771 if((dO1O2 > (R1+R2+Tol)) || (Rmax > (dO1O2+Rmin+Tol))) {
1772 typeres=IntAna_Empty;
1773 }
1774 else {
1775 //---------------------------------------------------------------
1776 //--
1777 //--
1778 //---------------------------------------------------------------
1779 Standard_Real Alpha=0.5*(R1*R1-R2*R2+dO1O2*dO1O2)/(dO1O2);
1780 Standard_Real Beta = R1*R1-Alpha*Alpha;
1781 Beta = (Beta>0.0)? Sqrt(Beta) : 0.0;
1782
1783 if(Beta<= myEPSILON_MINI_CIRCLE_RADIUS) {
1784 typeres = IntAna_Point;
1785 Alpha = (R1 + (dO1O2 - R2)) * 0.5;
1786 }
1787 else {
1788 typeres = IntAna_Circle;
1789 dir1 = Dir;
1790 param1 = Beta;
1791 }
1792 pt1.SetCoord( O1.X() + Alpha*Dir.X()
1793 ,O1.Y() + Alpha*Dir.Y()
1794 ,O1.Z() + Alpha*Dir.Z());
1795
1796 nbint=1;
1797 }
1798 }
1799 }
1800}
1801//=======================================================================
1802//function : Point
1803//purpose : Returns a Point
1804//=======================================================================
1805 gp_Pnt IntAna_QuadQuadGeo::Point(const Standard_Integer n) const
1806{
1807 if(!done) { StdFail_NotDone::Raise(); }
1808 if(n>nbint || n<1) { Standard_DomainError::Raise(); }
1809 if(typeres==IntAna_PointAndCircle) {
1810 if(n!=1) { Standard_DomainError::Raise(); }
1811 if(param1==0.0) return(pt1);
1812 return(pt2);
1813 }
1814 else if(typeres==IntAna_Point) {
1815 if(n==1) return(pt1);
1816 return(pt2);
1817 }
1818
1819 // WNT (what can you expect from MicroSoft ?)
1820 return gp_Pnt(0,0,0);
1821}
1822//=======================================================================
1823//function : Line
1824//purpose : Returns a Line
1825//=======================================================================
1826 gp_Lin IntAna_QuadQuadGeo::Line(const Standard_Integer n) const
1827{
1828 if(!done) { StdFail_NotDone::Raise(); }
1829 if((n>nbint) || (n<1) || (typeres!=IntAna_Line)) {
1830 Standard_DomainError::Raise();
1831 }
1832 if(n==1) { return(gp_Lin(pt1,dir1)); }
1833 else { return(gp_Lin(pt2,dir2)); }
1834}
1835//=======================================================================
1836//function : Circle
1837//purpose : Returns a Circle
1838//=======================================================================
1839 gp_Circ IntAna_QuadQuadGeo::Circle(const Standard_Integer n) const
1840{
1841 if(!done) { StdFail_NotDone::Raise(); }
1842 if(typeres==IntAna_PointAndCircle) {
1843 if(n!=1) { Standard_DomainError::Raise(); }
1844 if(param2==0.0) return(gp_Circ(DirToAx2(pt1,dir1),param1));
1845 return(gp_Circ(DirToAx2(pt2,dir2),param2));
1846 }
1847 else if((n>nbint) || (n<1) || (typeres!=IntAna_Circle)) {
1848 Standard_DomainError::Raise();
1849 }
1850 if(n==1) { return(gp_Circ(DirToAx2(pt1,dir1),param1)); }
1851 else { return(gp_Circ(DirToAx2(pt2,dir2),param2)); }
1852}
1853
1854//=======================================================================
1855//function : Ellipse
1856//purpose : Returns a Elips
1857//=======================================================================
1858 gp_Elips IntAna_QuadQuadGeo::Ellipse(const Standard_Integer n) const
1859{
1860 if(!done) { StdFail_NotDone::Raise(); }
1861 if((n>nbint) || (n<1) || (typeres!=IntAna_Ellipse)) {
1862 Standard_DomainError::Raise();
1863 }
1864
1865 if(n==1) {
1866 Standard_Real R1=param1, R2=param1bis, aTmp;
1867 if (R1<R2) {
1868 aTmp=R1; R1=R2; R2=aTmp;
1869 }
1870 gp_Ax2 anAx2(pt1, dir1 ,dir2);
1871 gp_Elips anElips (anAx2, R1, R2);
1872 return anElips;
1873 }
1874 else {
1875 Standard_Real R1=param2, R2=param2bis, aTmp;
1876 if (R1<R2) {
1877 aTmp=R1; R1=R2; R2=aTmp;
1878 }
1879 gp_Ax2 anAx2(pt2, dir2 ,dir1);
1880 gp_Elips anElips (anAx2, R1, R2);
1881 return anElips;
1882 }
1883}
1884//=======================================================================
1885//function : Parabola
1886//purpose : Returns a Parabola
1887//=======================================================================
1888 gp_Parab IntAna_QuadQuadGeo::Parabola(const Standard_Integer n) const
1889{
1890 if(!done) {
1891 StdFail_NotDone::Raise();
1892 }
1893 if (typeres!=IntAna_Parabola) {
1894 Standard_DomainError::Raise();
1895 }
1896 if((n>nbint) || (n!=1)) {
1897 Standard_OutOfRange::Raise();
1898 }
1899 return(gp_Parab(gp_Ax2( pt1
1900 ,dir1
1901 ,dir2)
1902 ,param1));
1903}
1904//=======================================================================
1905//function : Hyperbola
1906//purpose : Returns a Hyperbola
1907//=======================================================================
1908 gp_Hypr IntAna_QuadQuadGeo::Hyperbola(const Standard_Integer n) const
1909{
1910 if(!done) {
1911 StdFail_NotDone::Raise();
1912 }
1913 if((n>nbint) || (n<1) || (typeres!=IntAna_Hyperbola)) {
1914 Standard_DomainError::Raise();
1915 }
1916 if(n==1) {
1917 return(gp_Hypr(gp_Ax2( pt1
1918 ,dir1
1919 ,dir2)
1920 ,param1,param1bis));
1921 }
1922 else {
1923 return(gp_Hypr(gp_Ax2( pt2
1924 ,dir1
1925 ,dir2.Reversed())
1926 ,param2,param2bis));
1927 }
1928}
1929//=======================================================================
1930//function : HasCommonGen
1931//purpose :
1932//=======================================================================
1933Standard_Boolean IntAna_QuadQuadGeo::HasCommonGen() const
1934{
1935 return myCommonGen;
1936}
1937//=======================================================================
1938//function : PChar
1939//purpose :
1940//=======================================================================
1941const gp_Pnt& IntAna_QuadQuadGeo::PChar() const
1942{
1943 return myPChar;
1944}
1945//=======================================================================
1946//function : RefineDir
1947//purpose :
1948//=======================================================================
1949void RefineDir(gp_Dir& aDir)
1950{
1951 Standard_Integer k, m, n;
1952 Standard_Real aC[3];
1953 //
1954 aDir.Coord(aC[0], aC[1], aC[2]);
1955 //
1956 m=0;
1957 n=0;
1958 for (k=0; k<3; ++k) {
1959 if (aC[k]==1. || aC[k]==-1.) {
1960 ++m;
1961 }
1962 else if (aC[k]!=0.) {
1963 ++n;
1964 }
1965 }
1966 //
1967 if (m && n) {
1968 Standard_Real aEps, aR1, aR2, aNum;
1969 //
1970 aEps=RealEpsilon();
1971 aR1=1.-aEps;
1972 aR2=1.+aEps;
1973 //
1974 for (k=0; k<3; ++k) {
1975 m=(aC[k]>0.);
1976 aNum=(m)? aC[k] : -aC[k];
1977 if (aNum>aR1 && aNum<aR2) {
1978 if (m) {
1979 aC[k]=1.;
1980 }
1981 else {
1982 aC[k]=-1.;
1983 }
1984 //
1985 aC[(k+1)%3]=0.;
1986 aC[(k+2)%3]=0.;
1987 break;
1988 }
1989 }
1990 aDir.SetCoord(aC[0], aC[1], aC[2]);
1991 }
1992}
1993
1994
1995
Open CASCADE Technology repository