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