c03b3e9d8ce0e63bd5ca4a38a8df49f15702d38a
[occt.git] / src / GccGeo / GccGeo_Circ2d2TanRad.gxx
1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15 #include <ElCLib.hxx>
16 #include <gp_Ax2d.hxx>
17 #include <gp_Circ2d.hxx>
18 #include <gp_Lin2d.hxx>
19 #include <Standard_NegativeValue.hxx>
20 #include <Standard_OutOfRange.hxx>
21 #include <StdFail_NotDone.hxx>
22 #include <TColStd_Array1OfReal.hxx>
23 #include <GccEnt_BadQualifier.hxx>
24 #include <IntRes2d_Domain.hxx>
25 #include <IntRes2d_IntersectionPoint.hxx>
26
27 // circulaire tant a une courbe et une droite ,de rayon donne
28 //==============================================================
29
30 //========================================================================
31 // On initialise WellDone a false.                                       +
32 // On recupere la courbe Cu2 et la droite L1.                            +
33 // On sort en erreur dans les cas ou la construction est impossible.     +
34 // On fait la parallele a Cu2 dans le bon sens.                          +
35 // On fait la parallele a L1 dans le bon sens.                           +
36 // On intersecte les paralleles ==> point de centre de la solution.      +
37 // On cree la solution qu on ajoute aux solutions deja trouvees.         +
38 // On remplit les champs.                                                +
39 //========================================================================
40
41 GccGeo_Circ2d2TanRad::
42 GccGeo_Circ2d2TanRad (const GccEnt_QualifiedLin&  Qualified1,
43                       const TheQCurve&            Qualified2,
44                       const Standard_Real         Radius    ,
45                       const Standard_Real         Tolerance ):
46
47 //========================================================================
48 // initialisation des champs.                                            +
49 //========================================================================
50
51 cirsol(1,16)   ,
52 qualifier1(1,16),
53 qualifier2(1,16),
54 TheSame1(1,16) ,
55 TheSame2(1,16) ,
56 pnttg1sol(1,16),
57 pnttg2sol(1,16),
58 par1sol(1,16)  ,
59 par2sol(1,16)  ,
60 pararg1(1,16)  ,
61 pararg2(1,16)  
62 {
63
64   //========================================================================
65   // Traitement.                                                           +
66   //========================================================================
67
68   Standard_Real Tol = Abs(Tolerance);
69   Standard_Real thefirst = -100000.;
70   Standard_Real thelast  =  100000.;
71   Standard_Real firstparam;
72   Standard_Real lastparam;
73   gp_Dir2d dirx(1.,0.);
74   TColStd_Array1OfReal cote1(1,2);
75   TColStd_Array1OfReal cote2(1,2);
76   Standard_Integer nbrcote1=0;
77   Standard_Integer nbrcote2=0;
78   WellDone = Standard_False;
79   NbrSol = 0;
80   if (!(Qualified1.IsEnclosed() ||
81     Qualified1.IsOutside() || Qualified1.IsUnqualified()) ||
82     !(Qualified2.IsEnclosed() || Qualified2.IsEnclosing() || 
83     Qualified2.IsOutside() || Qualified2.IsUnqualified())) {
84
85       GccEnt_BadQualifier::Raise();
86       return;
87   }
88   gp_Lin2d L1 = Qualified1.Qualified();
89   Standard_Real x1dir = (L1.Direction()).X();
90   Standard_Real y1dir = (L1.Direction()).Y();
91   Standard_Real lxloc = (L1.Location()).X();
92   Standard_Real lyloc = (L1.Location()).Y();
93   gp_Pnt2d origin1(lxloc,lyloc);
94   gp_Dir2d normL1(-y1dir,x1dir);
95   TheCurve Cu2= Qualified2.Qualified();
96   if (Radius < 0.0) { Standard_NegativeValue::Raise(); }
97   else {
98     if (Qualified1.IsEnclosed() && Qualified2.IsEnclosed()) {
99       //   =======================================================
100       nbrcote1 = 1;
101       nbrcote2 = 1;
102       cote1(1) = Radius;
103       cote2(1) = Radius;
104     }
105     else if(Qualified1.IsEnclosed() && Qualified2.IsOutside()) {
106       //   ==========================================================
107       nbrcote1 = 1;
108       nbrcote2 = 1;
109       cote1(1) = Radius;
110       cote2(1) = -Radius;
111     }
112     else if (Qualified1.IsOutside() && Qualified2.IsEnclosed()) {
113       //   ===========================================================
114       nbrcote1 = 1;
115       nbrcote2 = 1;
116       cote1(1) = -Radius;
117       cote2(1) = Radius;
118     }
119     else if(Qualified1.IsOutside() && Qualified2.IsOutside()) {
120       //   =========================================================
121       nbrcote1 = 1;
122       nbrcote2 = 1;
123       cote1(1) = -Radius;
124       cote2(1) = -Radius;
125     }
126     if(Qualified1.IsEnclosed() && Qualified2.IsUnqualified()) {
127       //   =========================================================
128       nbrcote1 = 1;
129       nbrcote2 = 2;
130       cote1(1) = Radius;
131       cote2(1) = Radius;
132       cote2(2) = -Radius;
133     }
134     if(Qualified1.IsUnqualified() && Qualified2.IsEnclosed()) {
135       //   =========================================================
136       nbrcote1 = 2;
137       nbrcote2 = 1;
138       cote1(1) = Radius;
139       cote1(2) = -Radius;
140       cote2(1) = Radius;
141     }
142     else if(Qualified1.IsOutside() && Qualified2.IsUnqualified()) {
143       //   =============================================================
144       nbrcote1 = 1;
145       nbrcote2 = 2;
146       cote1(1) = -Radius;
147       cote2(1) = Radius;
148       cote2(2) = -Radius;
149     }
150     if(Qualified1.IsUnqualified() && Qualified2.IsOutside()) {
151       //   ========================================================
152       nbrcote1 = 2;
153       nbrcote2 = 1;
154       cote1(1) = Radius;
155       cote1(2) = -Radius;
156       cote2(1) = -Radius;
157     }
158     else if(Qualified1.IsUnqualified() && Qualified2.IsUnqualified()) {
159       //   =================================================================
160       nbrcote1 = 2;
161       nbrcote2 = 2;
162       cote1(1) = Radius;
163       cote1(2) = -Radius;
164       cote2(1) = Radius;
165       cote2(2) = -Radius;
166     }
167     gp_Dir2d Dir(-y1dir,x1dir);
168     for (Standard_Integer jcote1 = 1 ; jcote1 <= nbrcote1 ; jcote1++) {
169       gp_Pnt2d Point(L1.Location().XY()+cote1(jcote1)*Dir.XY());
170       gp_Lin2d Line(Point,L1.Direction()); // ligne avec deport.
171       IntRes2d_Domain D1;
172       for (Standard_Integer jcote2 = 1 ; jcote2 <= nbrcote2 ; jcote2++) {
173         Handle(TheHParGenCurve) HCu2 = new TheHParGenCurve(Cu2);
174         TheParGenCurve C2(HCu2,cote2(jcote2));
175         firstparam = Max(TheCurvePGTool::FirstParameter(C2),thefirst);
176         lastparam  = Min(TheCurvePGTool::LastParameter(C2),thelast);
177         IntRes2d_Domain D2(TheCurvePGTool::Value(C2,firstparam),firstparam,Tol,
178           TheCurvePGTool::Value(C2,lastparam),lastparam,Tol);
179         TheIntConicCurve Intp(Line,D1,C2,D2,Tol,Tol);
180         if (Intp.IsDone()) {
181           if (!Intp.IsEmpty()) {
182             for (Standard_Integer i = 1 ; i <= Intp.NbPoints() ; i++) {
183               NbrSol++;
184               gp_Pnt2d Center(Intp.Point(i).Value());
185               cirsol(NbrSol) = gp_Circ2d(gp_Ax2d(Center,dirx),Radius);
186               //             =======================================================
187               gp_Dir2d dc1(origin1.XY()-Center.XY());
188               qualifier2(NbrSol) = Qualified2.Qualifier();
189               if (!Qualified1.IsUnqualified()) { 
190                 qualifier1(NbrSol) = Qualified1.Qualifier();
191               }
192               else if (dc1.Dot(normL1) > 0.0) {
193                 qualifier1(NbrSol) = GccEnt_outside;
194               }
195               else { qualifier1(NbrSol) = GccEnt_enclosed; }
196               TheSame1(NbrSol) = 0;
197               TheSame2(NbrSol) = 0;
198               pararg1(NbrSol) = Intp.Point(i).ParamOnFirst();
199               pararg2(NbrSol) = Intp.Point(i).ParamOnSecond();
200               pnttg1sol(NbrSol) = ElCLib::Value(pararg1(NbrSol),L1);
201               pnttg2sol(NbrSol) = TheTool::Value(Cu2,pararg2(NbrSol));
202               par1sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
203                 pnttg1sol(NbrSol));
204               par2sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
205                 pnttg2sol(NbrSol));
206             }
207           }
208           WellDone = Standard_True;
209         }
210       }
211     }
212   }
213 }
214
215 // circulaire tant a une courbe et un cercle ,de rayon donne
216 //=============================================================
217
218 //========================================================================
219 // On initialise WellDone a false.                                       +
220 // On recupere la courbe Cu2 et le cercle C1.                            +
221 // On sort en erreur dans les cas ou la construction est impossible.     +
222 // On fait la parallele a Cu2 dans le bon sens.                          +
223 // On fait la parallele a C1 dans le bon sens.                           +
224 // On intersecte les paralleles ==> point de centre de la solution.      +
225 // On cree la solution qu on ajoute aux solutions deja trouvees.         +
226 // On remplit les champs.                                                +
227 //========================================================================
228
229 GccGeo_Circ2d2TanRad::
230 GccGeo_Circ2d2TanRad (const GccEnt_QualifiedCirc& Qualified1,
231                       const TheQCurve&            Qualified2,
232                       const Standard_Real         Radius    ,
233                       const Standard_Real         Tolerance ):
234
235 //========================================================================
236 // initialisation des champs.                                            +
237 //========================================================================
238
239 cirsol(1,16)   ,
240 qualifier1(1,16),
241 qualifier2(1,16),
242 TheSame1(1,16) ,
243 TheSame2(1,16) ,
244 pnttg1sol(1,16),
245 pnttg2sol(1,16),
246 par1sol(1,16)  ,
247 par2sol(1,16)  ,
248 pararg1(1,16)  ,
249 pararg2(1,16)  
250 {
251
252   //========================================================================
253   // Traitement.                                                           +
254   //========================================================================
255
256   Standard_Real Tol = Abs(Tolerance);
257   Standard_Real thefirst = -100000.;
258   Standard_Real thelast  =  100000.;
259   Standard_Real firstparam;
260   Standard_Real lastparam;
261   gp_Dir2d dirx(1.,0.);
262   TColStd_Array1OfReal cote1(1,2);
263   TColStd_Array1OfReal cote2(1,2);
264   Standard_Integer nbrcote1=0;
265   Standard_Integer nbrcote2=0;
266   WellDone = Standard_False;
267   NbrSol = 0;
268   if (!(Qualified1.IsEnclosed() || Qualified1.IsEnclosing() || 
269     Qualified1.IsOutside() || Qualified1.IsUnqualified()) ||
270     !(Qualified2.IsEnclosed() || Qualified2.IsEnclosing() || 
271     Qualified2.IsOutside() || Qualified2.IsUnqualified())) {
272       GccEnt_BadQualifier::Raise();
273       return;
274   }
275   gp_Circ2d C1 = Qualified1.Qualified();
276   gp_Pnt2d center1(C1.Location());
277   TheCurve Cu2 = Qualified2.Qualified();
278   if (Radius < 0.0) { Standard_NegativeValue::Raise(); }
279   else {
280     if (Qualified1.IsEnclosed() && Qualified2.IsEnclosed()) {
281       //   =======================================================
282       nbrcote1 = 1;
283       nbrcote2 = 1;
284       cote1(1) = Radius;
285       cote2(1) = Radius;
286     }
287     else if(Qualified1.IsEnclosed() && Qualified2.IsOutside()) {
288       //   ==========================================================
289       nbrcote1 = 1;
290       nbrcote2 = 1;
291       cote1(1) = Radius;
292       cote2(1) = -Radius;
293     }
294     else if (Qualified1.IsOutside() && Qualified2.IsEnclosed()) {
295       //   ===========================================================
296       nbrcote1 = 1;
297       nbrcote2 = 1;
298       cote1(1) = -Radius;
299       cote2(1) = Radius;
300     }
301     else if(Qualified1.IsOutside() && Qualified2.IsOutside()) {
302       //   =========================================================
303       nbrcote1 = 1;
304       nbrcote2 = 1;
305       cote1(1) = -Radius;
306       cote2(1) = -Radius;
307     }
308     if(Qualified1.IsEnclosed() && Qualified2.IsUnqualified()) {
309       //   =========================================================
310       nbrcote1 = 1;
311       nbrcote2 = 2;
312       cote1(1) = Radius;
313       cote2(1) = Radius;
314       cote2(2) = -Radius;
315     }
316     if(Qualified1.IsUnqualified() && Qualified2.IsEnclosed()) {
317       //   =========================================================
318       nbrcote1 = 2;
319       nbrcote2 = 1;
320       cote1(1) = Radius;
321       cote1(2) = -Radius;
322       cote2(1) = Radius;
323     }
324     else if(Qualified1.IsOutside() && Qualified2.IsUnqualified()) {
325       //   =============================================================
326       nbrcote1 = 1;
327       nbrcote2 = 2;
328       cote1(1) = -Radius;
329       cote2(1) = Radius;
330       cote2(2) = -Radius;
331     }
332     if(Qualified1.IsUnqualified() && Qualified2.IsOutside()) {
333       //   ========================================================
334       nbrcote1 = 2;
335       nbrcote2 = 1;
336       cote1(1) = Radius;
337       cote1(2) = -Radius;
338       cote2(1) = -Radius;
339     }
340     else if(Qualified1.IsUnqualified() && Qualified2.IsUnqualified()) {
341       //   =================================================================
342       nbrcote1 = 2;
343       nbrcote2 = 2;
344       cote1(1) = Radius;
345       cote1(2) = -Radius;
346       cote2(1) = Radius;
347       cote2(2) = -Radius;
348     }
349     Standard_Real R1 = C1.Radius();
350     TheIntConicCurve Intp;
351     for (Standard_Integer jcote1 = 1 ; jcote1 <= nbrcote1 ; jcote1++) {
352       gp_Circ2d Circ(C1.XAxis(),R1+cote1(jcote1));
353       IntRes2d_Domain D1(ElCLib::Value(0.,Circ),   0.,Tol,
354         ElCLib::Value(2.*M_PI,Circ),2.*M_PI,Tol);
355       D1.SetEquivalentParameters(0.,2.*M_PI);
356       for (Standard_Integer jcote2 = 1 ; jcote2 <= nbrcote2 ; jcote2++) {
357         Handle(TheHParGenCurve) HCu2 = new TheHParGenCurve(Cu2);
358         TheParGenCurve C2(HCu2,cote2(jcote2));
359         firstparam = Max(TheCurvePGTool::FirstParameter(C2),thefirst);
360         lastparam  = Min(TheCurvePGTool::LastParameter(C2),thelast);
361         IntRes2d_Domain D2(TheCurvePGTool::Value(C2,firstparam),firstparam,Tol,
362           TheCurvePGTool::Value(C2,lastparam),lastparam,Tol);
363         Intp.Perform(Circ,D1,C2,D2,Tol,Tol);
364         if (Intp.IsDone()) {
365           if (!Intp.IsEmpty()) {
366             for (Standard_Integer i = 1 ; i <= Intp.NbPoints() ; i++) {
367               NbrSol++;
368               gp_Pnt2d Center(Intp.Point(i).Value());
369               cirsol(NbrSol) = gp_Circ2d(gp_Ax2d(Center,dirx),Radius);
370               //             =======================================================
371 #ifdef DEB
372               gp_Dir2d dir1(Center.XY()-center1.XY());
373 #else
374               Center.XY() ;
375               center1.XY() ;
376 #endif
377               Standard_Real distcc1 = Center.Distance(center1);
378               if (!Qualified1.IsUnqualified()) { 
379                 qualifier1(NbrSol) = Qualified1.Qualifier();
380               }
381               else if (Abs(distcc1+Radius-R1) < Tol) {
382                 qualifier1(NbrSol) = GccEnt_enclosed;
383               }
384               else if (Abs(distcc1-R1-Radius) < Tol) {
385                 qualifier1(NbrSol) = GccEnt_outside;
386               }
387               else { qualifier1(NbrSol) = GccEnt_enclosing; }
388               qualifier2(NbrSol) = Qualified2.Qualifier();
389               TheSame1(NbrSol) = 0;
390               TheSame2(NbrSol) = 0;
391               pararg1(NbrSol) = Intp.Point(i).ParamOnFirst();
392               pararg2(NbrSol) = Intp.Point(i).ParamOnSecond();
393               pnttg1sol(NbrSol) = ElCLib::Value(pararg1(NbrSol),C1);
394               pnttg2sol(NbrSol) = TheTool::Value(Cu2,pararg2(NbrSol));
395               par1sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
396                 pnttg1sol(NbrSol));
397               par2sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
398                 pnttg2sol(NbrSol));
399             }
400           }
401           WellDone = Standard_True;
402         }
403       }
404     }
405   }
406 }
407
408 // circulaire tant a une courbe et un point ,de rayon donne
409 //============================================================
410
411 //========================================================================
412 // On initialise WellDone a false.                                       +
413 // On recupere la courbe Cu1 et le point P2.                             +
414 // On sort en erreur dans les cas ou la construction est impossible.     +
415 // On fait la parallele a Cu1 dans le bon sens.                          +
416 // On fait la parallele a P2 dans le bon sens.                           +
417 // On intersecte les paralleles ==> point de centre de la solution.      +
418 // On cree la solution qu on ajoute aux solutions deja trouvees.         +
419 // On remplit les champs.                                                +
420 //========================================================================
421
422 GccGeo_Circ2d2TanRad::
423 GccGeo_Circ2d2TanRad (const TheQCurve&    Qualified1,
424                       const gp_Pnt2d&     Point2    ,
425                       const Standard_Real Radius    ,
426                       const Standard_Real Tolerance ):
427
428 //========================================================================
429 // initialisation des champs.                                            +
430 //========================================================================
431
432 cirsol(1,16)   ,
433 qualifier1(1,16),
434 qualifier2(1,16),
435 TheSame1(1,16) ,
436 TheSame2(1,16) ,
437 pnttg1sol(1,16),
438 pnttg2sol(1,16),
439 par1sol(1,16)  ,
440 par2sol(1,16)  ,
441 pararg1(1,16)  ,
442 pararg2(1,16)  
443 {
444
445   //========================================================================
446   // Traitement.                                                           +
447   //========================================================================
448
449   Standard_Real Tol = Abs(Tolerance);
450   Standard_Real thefirst = -100000.;
451   Standard_Real thelast  =  100000.;
452   Standard_Real firstparam;
453   Standard_Real lastparam;
454   gp_Dir2d dirx(1.,0.);
455   TColStd_Array1OfReal cote1(1,2);
456   Standard_Integer nbrcote1=0;
457   WellDone = Standard_False;
458   NbrSol = 0;
459   if (!(Qualified1.IsEnclosed() || Qualified1.IsEnclosing() || 
460     Qualified1.IsOutside() || Qualified1.IsUnqualified())) {
461       GccEnt_BadQualifier::Raise();
462       return;
463   }
464   TheCurve Cu1 = Qualified1.Qualified();
465   if (Radius < 0.0) { Standard_NegativeValue::Raise(); }
466   else {
467     if (Qualified1.IsEnclosed()) {
468       //    ===========================
469       nbrcote1 = 1;
470       cote1(1) = Radius;
471     }
472     else if(Qualified1.IsOutside()) {
473       //   ===============================
474       nbrcote1 = 1;
475       cote1(1) = -Radius;
476     }
477     else if(Qualified1.IsUnqualified()) {
478       //   ===================================
479       nbrcote1 = 2;
480       cote1(1) = Radius;
481       cote1(2) = -Radius;
482     }
483     gp_Circ2d Circ(gp_Ax2d(Point2,gp_Dir2d(1.,0.)),Radius);
484     IntRes2d_Domain D1(ElCLib::Value(0.,Circ),   0.,Tol,
485       ElCLib::Value(M_PI+M_PI,Circ),M_PI+M_PI,Tol);
486     D1.SetEquivalentParameters(0.,M_PI+M_PI);
487     TheIntConicCurve Intp;
488     for (Standard_Integer jcote1 = 1 ; jcote1 <= nbrcote1 ; jcote1++) {
489       Handle(TheHParGenCurve) HCu1 = new TheHParGenCurve(Cu1);
490       TheParGenCurve Cu2(HCu1,cote1(jcote1));
491       firstparam = Max(TheCurvePGTool::FirstParameter(Cu2),thefirst);
492       lastparam  = Min(TheCurvePGTool::LastParameter(Cu2),thelast);
493       IntRes2d_Domain D2(TheCurvePGTool::Value(Cu2,firstparam),firstparam,Tol,
494         TheCurvePGTool::Value(Cu2,lastparam),lastparam,Tol);
495       Intp.Perform(Circ,D1,Cu2,D2,Tol,Tol);
496       if (Intp.IsDone()) {
497         if (!Intp.IsEmpty()) {
498           for (Standard_Integer i = 1 ; i <= Intp.NbPoints() ; i++) {
499             NbrSol++;
500             gp_Pnt2d Center(Intp.Point(i).Value());
501             cirsol(NbrSol) = gp_Circ2d(gp_Ax2d(Center,dirx),Radius);
502             //           =======================================================
503             qualifier1(NbrSol) = Qualified1.Qualifier();
504             qualifier2(NbrSol) = GccEnt_noqualifier;
505             TheSame1(NbrSol) = 0;
506             TheSame2(NbrSol) = 0;
507             pararg1(NbrSol) = Intp.Point(i).ParamOnSecond();
508             pararg2(NbrSol) = 0.;
509             pnttg1sol(NbrSol) = TheTool::Value(Cu1,pararg1(NbrSol));
510             pnttg2sol(NbrSol) = Point2;
511             par1sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
512               pnttg1sol(NbrSol));
513             par2sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
514               pnttg2sol(NbrSol));
515           }
516         }
517         WellDone = Standard_True;
518       }
519     }
520   }
521 }
522
523 static void PrecRoot(const TheParGenCurve& theC1,
524                                  const TheParGenCurve& theC2,
525                                  const Standard_Real theU0,
526                                  const Standard_Real theV0,
527                                  const Standard_Real theUmin,
528                                  const Standard_Real theUmax,
529                                  const Standard_Real theVmin,
530                                  const Standard_Real theVmax,
531                                  Standard_Real& theUfinal,
532                                  Standard_Real& theVfinal)
533 {
534   const Standard_Real aInitStepU = (theUmax - theUmin)/2.0,
535                       aInitStepV = (theVmax - theVmin)/2.0;
536
537   Standard_Real aStepU = aInitStepU, aStepV = aInitStepV;
538
539   const Standard_Real aTol = Precision::PConfusion() * Precision::PConfusion();
540   const Standard_Integer aNbIterMax = 100;
541
542   gp_Pnt2d aP1, aP2;
543   gp_Vec2d aD1, aD2;
544
545   TheCurvePGTool::D1(theC1, theU0, aP1, aD1);
546   TheCurvePGTool::D1(theC2, theV0, aP2, aD2);
547
548   gp_Vec2d vP12(aP1.XY() - aP2.XY());
549
550   Standard_Real aU = theU0, aV = theV0;
551   theUfinal = theU0;
552   theVfinal = theV0;
553
554   Standard_Real aSQDistPrev = aP1.SquareDistance(aP2);
555   
556   Standard_Integer aNbIter = 1;
557
558   do
559   {
560     Standard_Real aDetH = aD1.Y()*aD2.X() - aD1.X()*aD2.Y();
561     if(aDetH == 0.0)
562       break;
563
564     aU += aStepU*(aD2.Y() * vP12.X() - aD2.X()*vP12.Y())/aDetH;
565     aV += aStepV*(aD1.Y() * vP12.X() - aD1.X()*vP12.Y())/aDetH;
566
567     if(Abs(aU - theUmin) > 1000.0)
568       //method diverges
569       return;
570
571     if(Abs(aU - theUmax) > 1000.0)
572       //method diverges
573       return;
574
575     if(Abs(aV - theVmin) > 1000.0)
576       //method diverges
577       return;
578
579     if(Abs(aV - theVmax) > 1000.0)
580       //method diverges
581       return;
582
583     TheCurvePGTool::D1(theC1, aU, aP1, aD1);
584     TheCurvePGTool::D1(theC2, aV, aP2, aD2);
585     const Standard_Real aSQDist = aP1.SquareDistance(aP2);
586
587     if(Precision::IsInfinite(aSQDist))
588       //method diverges
589       return;
590
591     vP12.SetXY(aP1.XY() - aP2.XY());
592
593     if(aSQDist < aSQDistPrev)
594     {
595       aSQDistPrev = aSQDist;
596       aStepU = aInitStepU;
597       aStepV = aInitStepV;
598       theUfinal = aU;
599       theVfinal = aV;
600     }
601     else
602     {
603       aStepU /= 2.0;
604       aStepV /= 2.0;
605     }
606   }
607   while((aNbIter++ < aNbIterMax) && ((aStepU > aTol) || (aStepV > aTol)));
608
609   Standard_Boolean isInBound = Standard_True;
610   if(theUfinal < theUmin)
611   {
612     aU = theUfinal;
613     aV = theVfinal;
614
615     theUfinal = theUmin;
616     isInBound = Standard_False;
617   }
618
619   if(theUfinal > theUmax)
620   {
621     aU = theUfinal;
622     aV = theVfinal;
623
624     theUfinal = theUmax;
625     isInBound = Standard_False;
626   }
627
628   if(!isInBound)
629   {
630     TheCurvePGTool::D1(theC1, aU, aP1, aD1);
631     TheCurvePGTool::D1(theC2, aV, aP2, aD2);
632     Standard_Real aV1 = (aD2.X() == 0.0) ? aV :((theUfinal - aU)*aD1.X() + aV*aD2.X() + (aP1.X() - aP2.X()))/aD2.X();
633     Standard_Real aV2 = (aD2.Y() == 0.0) ? aV :((theUfinal - aU)*aD1.Y() + aV*aD2.Y() + (aP1.Y() - aP2.Y()))/aD2.Y();
634
635     if(aV1 < theVmin)
636       aV1 = theVmin;
637
638     if(aV1 > theVmax)
639       aV1 = theVmax;
640
641     if(aV2 < theVmin)
642       aV2 = theVmin;
643
644     if(aV2 > theVmax)
645       aV2 = theVmax;
646
647     aP1 = TheCurvePGTool::Value(theC1,theUfinal);
648     aP2 = TheCurvePGTool::Value(theC2,aV1);
649
650     Standard_Real aSQ1 = aP1.SquareDistance(aP2);
651
652     aP2 = TheCurvePGTool::Value(theC2,aV2);
653     Standard_Real aSQ2 = aP1.SquareDistance(aP2);
654
655     if(aSQ1 < aSQ2)
656       theVfinal = aV1;
657     else
658       theVfinal = aV2;
659
660     return;
661   }
662
663   if(theVfinal < theVmin)
664   {
665     aU = theUfinal;
666     aV = theVfinal;
667
668     theVfinal = theVmin;
669     isInBound = Standard_False;
670   }
671
672   if(theVfinal > theVmax)
673   {
674     aU = theUfinal;
675     aV = theVfinal;
676
677     theVfinal = theVmax;
678     isInBound = Standard_False;
679   }
680
681   if(isInBound)
682     return;
683
684   TheCurvePGTool::D1(theC1, aU, aP1, aD1);
685   TheCurvePGTool::D1(theC2, aV, aP2, aD2);
686   Standard_Real aU1 = (aD1.X() == 0.0) ? aU :((theVfinal - aV)*aD2.X() + aU*aD1.X() + (aP2.X() - aP1.X()))/aD1.X();
687   Standard_Real aU2 = (aD1.Y() == 0.0) ? aU :((theVfinal - aV)*aD2.Y() + aU*aD1.Y() + (aP2.Y() - aP1.Y()))/aD1.Y();
688
689   if(aU1 < theUmin)
690     aU1 = theUmin;
691
692   if(aU1 > theUmax)
693     aU1 = theUmax;
694
695   if(aU2 < theUmin)
696     aU2 = theUmin;
697
698   if(aU2 > theUmax)
699     aU2 = theUmax;
700
701   aP2 = TheCurvePGTool::Value(theC2,theVfinal);
702   aP1 = TheCurvePGTool::Value(theC1,aU1);
703
704   Standard_Real aSQ1 = aP1.SquareDistance(aP2);
705
706   aP1 = TheCurvePGTool::Value(theC1,aU2);
707   Standard_Real aSQ2 = aP1.SquareDistance(aP2);
708
709   if(aSQ1 < aSQ2)
710     theUfinal = aU1;
711   else
712     theUfinal = aU2;  
713 }
714
715 // circulaire tant a deux courbes ,de rayon donne
716 //==================================================
717
718 //========================================================================
719 // On initialise WellDone a false.                                       +
720 // On recupere les courbes Cu1 et Cu2.                                   +
721 // On sort en erreur dans les cas ou la construction est impossible.     +
722 // On fait la parallele a Cu1 dans le bon sens.                          +
723 // On fait la parallele a Cu2 dans le bon sens.                          +
724 // On intersecte les paralleles ==> point de centre de la solution.      +
725 // On cree la solution qu on ajoute aux solutions deja trouvees.         +
726 // On remplit les champs.                                                +
727 //========================================================================
728
729 GccGeo_Circ2d2TanRad::
730 GccGeo_Circ2d2TanRad (const TheQCurve&    Qualified1,
731                       const TheQCurve&    Qualified2,
732                       const Standard_Real Radius    ,
733                       const Standard_Real Tolerance ):
734
735 //========================================================================
736 // initialisation des champs.                                            +
737 //========================================================================
738
739 cirsol(1,16)   ,
740 qualifier1(1,16),
741 qualifier2(1,16),
742 TheSame1(1,16) ,
743 TheSame2(1,16) ,
744 pnttg1sol(1,16),
745 pnttg2sol(1,16),
746 par1sol(1,16)  ,
747 par2sol(1,16)  ,
748 pararg1(1,16)  ,
749 pararg2(1,16)  
750 {
751
752   //========================================================================
753   // Traitement.                                                           +
754   //========================================================================
755
756   Standard_Real Tol = Abs(Tolerance);
757   Standard_Real thefirst = -100000.;
758   Standard_Real thelast  =  100000.;
759   Standard_Real firstparam;
760   Standard_Real lastparam;
761   gp_Dir2d dirx(1.,0.);
762   TColStd_Array1OfReal cote1(1,2);
763   TColStd_Array1OfReal cote2(1,2);
764   Standard_Integer nbrcote1=0;
765   Standard_Integer nbrcote2=0;
766   WellDone = Standard_False;
767   NbrSol = 0;
768   if (!(Qualified1.IsEnclosed() || Qualified1.IsEnclosing() || 
769     Qualified1.IsOutside() || Qualified1.IsUnqualified()) ||
770     !(Qualified2.IsEnclosed() || Qualified2.IsEnclosing() || 
771     Qualified2.IsOutside() || Qualified2.IsUnqualified())) {
772       GccEnt_BadQualifier::Raise();
773       return;
774   }
775   TheCurve Cu1 = Qualified1.Qualified();
776   TheCurve Cu2 = Qualified2.Qualified();
777   if (Radius < 0.0) { Standard_NegativeValue::Raise(); }
778   else {
779     if (Qualified1.IsEnclosed() && Qualified2.IsEnclosed()) {
780       //   =======================================================
781       nbrcote1 = 1;
782       nbrcote2 = 1;
783       cote1(1) = Radius;
784       cote2(1) = Radius;
785     }
786     else if(Qualified1.IsEnclosed() && Qualified2.IsOutside()) {
787       //   ==========================================================
788       nbrcote1 = 1;
789       nbrcote2 = 1;
790       cote1(1) = Radius;
791       cote2(1) = -Radius;
792     }
793     else if (Qualified1.IsOutside() && Qualified2.IsEnclosed()) {
794       //   ===========================================================
795       nbrcote1 = 1;
796       nbrcote2 = 1;
797       cote1(1) = -Radius;
798       cote2(1) = Radius;
799     }
800     else if(Qualified1.IsOutside() && Qualified2.IsOutside()) {
801       //   =========================================================
802       nbrcote1 = 1;
803       nbrcote2 = 1;
804       cote1(1) = -Radius;
805       cote2(1) = -Radius;
806     }
807     if(Qualified1.IsEnclosed() && Qualified2.IsUnqualified()) {
808       //   =========================================================
809       nbrcote1 = 1;
810       nbrcote2 = 2;
811       cote1(1) = Radius;
812       cote2(1) = Radius;
813       cote2(2) = -Radius;
814     }
815     if(Qualified1.IsUnqualified() && Qualified2.IsEnclosed()) {
816       //   =========================================================
817       nbrcote1 = 2;
818       nbrcote2 = 1;
819       cote1(1) = Radius;
820       cote1(2) = -Radius;
821       cote2(1) = Radius;
822     }
823     else if(Qualified1.IsOutside() && Qualified2.IsUnqualified()) {
824       //   =============================================================
825       nbrcote1 = 1;
826       nbrcote2 = 2;
827       cote1(1) = -Radius;
828       cote2(1) = Radius;
829       cote2(2) = -Radius;
830     }
831     if(Qualified1.IsUnqualified() && Qualified2.IsOutside()) {
832       //   ========================================================
833       nbrcote1 = 2;
834       nbrcote2 = 1;
835       cote1(1) = Radius;
836       cote1(2) = -Radius;
837       cote2(1) = -Radius;
838     }
839     else if(Qualified1.IsUnqualified() && Qualified2.IsUnqualified()) {
840       //   =================================================================
841       nbrcote1 = 2;
842       nbrcote2 = 2;
843       cote1(1) = Radius;
844       cote1(2) = -Radius;
845       cote2(1) = Radius;
846       cote2(2) = -Radius;
847     }
848     TheIntCurveCurve Intp;
849     for (Standard_Integer jcote1 = 1 ; jcote1 <= nbrcote1 ; jcote1++) {
850       Handle(TheHParGenCurve) HCu1 = new TheHParGenCurve(Cu1); 
851       TheParGenCurve C1(HCu1,cote1(jcote1));
852       firstparam = Max(TheCurvePGTool::FirstParameter(C1),thefirst);
853       lastparam  = Min(TheCurvePGTool::LastParameter(C1),thelast);
854 #ifdef DEB
855       IntRes2d_Domain D2(TheCurvePGTool::Value(C1,firstparam),firstparam,Tol,
856         TheCurvePGTool::Value(C1,lastparam),lastparam,Tol);
857 #else
858       TheCurvePGTool::Value(C1,firstparam);
859       TheCurvePGTool::Value(C1,lastparam);
860 #endif
861       for (Standard_Integer jcote2 = 1 ; jcote2 <= nbrcote2 ; jcote2++) {
862         Handle(TheHParGenCurve) HCu2 = new TheHParGenCurve(Cu2);
863         TheParGenCurve C2(HCu2,cote2(jcote2));
864         firstparam = Max(TheCurvePGTool::FirstParameter(C2),thefirst);
865         lastparam  = Min(TheCurvePGTool::LastParameter(C2),thelast);
866 #ifdef DEB
867         IntRes2d_Domain D2(TheCurvePGTool::Value(C2,firstparam),firstparam,Tol,
868           TheCurvePGTool::Value(C2,lastparam),lastparam,Tol);
869 #else
870         TheCurvePGTool::Value(C2,firstparam);
871         TheCurvePGTool::Value(C2,lastparam);
872 #endif
873         Intp.Perform(C1,C2,Tol,Tol);
874         if (Intp.IsDone()) {
875           if (!Intp.IsEmpty()) {
876             for (Standard_Integer i = 1 ; i <= Intp.NbPoints() ; i++)
877             {
878               Standard_Real aU0 = Intp.Point(i).ParamOnFirst();
879               Standard_Real aV0 = Intp.Point(i).ParamOnSecond();
880
881               Standard_Real aU1 = aU0-Precision::PApproximation();
882               Standard_Real aV1 = aV0-Precision::PApproximation();
883
884               Standard_Real aU2 = aU0+Precision::PApproximation();
885               Standard_Real aV2 = aV0+Precision::PApproximation();
886
887               gp_Pnt2d P11 = TheCurvePGTool::Value(C1,aU1);
888               gp_Pnt2d P12 = TheCurvePGTool::Value(C2,aV1);
889               gp_Pnt2d P21 = TheCurvePGTool::Value(C1,aU2);
890               gp_Pnt2d P22 = TheCurvePGTool::Value(C2,aV2);
891
892               Standard_Real aDist1112 = P11.Distance(P12);
893               Standard_Real aDist1122 = P11.Distance(P22);
894
895               Standard_Real aDist1221 = P12.Distance(P21);
896               Standard_Real aDist2122 = P21.Distance(P22);
897
898               if( Min(aDist1112, aDist1122) <= Precision::Approximation() &&
899                   Min(aDist1221, aDist2122) <= Precision::Approximation())
900               {
901                 PrecRoot(C1, C2, aU0, aV0,
902                          Max(TheCurvePGTool::FirstParameter(C1), aU0 - 10.0),
903                          Min(TheCurvePGTool::LastParameter(C1), aU0 + 10.0),
904                          Max(TheCurvePGTool::FirstParameter(C2), aV0 - 10.0),
905                          Min(TheCurvePGTool::LastParameter(C2), aV0 + 10.0),
906                          aU0, aV0);
907               }
908
909               NbrSol++;
910               gp_Pnt2d Center(TheCurvePGTool::Value(C1, aU0));
911               cirsol(NbrSol) = gp_Circ2d(gp_Ax2d(Center,dirx),Radius);
912               //             =======================================================
913               qualifier1(NbrSol) = Qualified1.Qualifier();
914               qualifier1(NbrSol) = Qualified1.Qualifier();
915               TheSame1(NbrSol) = 0;
916               TheSame2(NbrSol) = 0;
917               pararg1(NbrSol) = Intp.Point(i).ParamOnFirst();
918               pararg2(NbrSol) = Intp.Point(i).ParamOnSecond();
919               pnttg1sol(NbrSol) = TheTool::Value(Cu1,pararg1(NbrSol));
920               pnttg2sol(NbrSol) = TheTool::Value(Cu2,pararg2(NbrSol));
921               par1sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
922                 pnttg1sol(NbrSol));
923               par2sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
924                 pnttg2sol(NbrSol));
925             }
926           }
927
928           WellDone = Standard_True;
929         }
930       }
931     }
932   }
933 }
934
935 //=========================================================================
936
937 Standard_Boolean GccGeo_Circ2d2TanRad::
938 IsDone () const { return WellDone; }
939
940 Standard_Integer GccGeo_Circ2d2TanRad::
941 NbSolutions () const { return NbrSol; }
942
943 gp_Circ2d GccGeo_Circ2d2TanRad::
944 ThisSolution (const Standard_Integer Index) const 
945 {
946   if (!WellDone) { StdFail_NotDone::Raise(); }
947   if (Index <= 0 ||Index > NbrSol) { Standard_OutOfRange::Raise(); }
948   return cirsol(Index);
949 }
950
951 void GccGeo_Circ2d2TanRad::
952 WhichQualifier(const Standard_Integer Index   ,
953                GccEnt_Position& Qualif1 ,
954                GccEnt_Position& Qualif2 ) const
955 {
956   if (!WellDone) { StdFail_NotDone::Raise(); }
957   else if (Index <= 0 ||Index > NbrSol) { Standard_OutOfRange::Raise(); }
958   else {
959     Qualif1 = qualifier1(Index);
960     Qualif2 = qualifier2(Index);
961   }
962 }
963
964 void GccGeo_Circ2d2TanRad::
965 Tangency1 (const Standard_Integer Index,
966            Standard_Real&   ParSol,
967            Standard_Real&   ParArg,
968            gp_Pnt2d&        PntSol) const{
969              if (!WellDone) { StdFail_NotDone::Raise(); }
970              else if (Index <= 0 ||Index > NbrSol) { Standard_OutOfRange::Raise(); }
971              else {
972                if (TheSame1(Index) == 0) {
973                  ParSol = par1sol(Index);
974                  ParArg = pararg1(Index);
975                  PntSol = gp_Pnt2d(pnttg1sol(Index));
976                }
977                else { StdFail_NotDone::Raise(); }
978              }
979 }
980
981 void GccGeo_Circ2d2TanRad::
982 Tangency2 (const Standard_Integer Index,
983            Standard_Real&   ParSol,
984            Standard_Real&   ParArg,
985            gp_Pnt2d&        PntSol) const{
986              if (!WellDone) { StdFail_NotDone::Raise(); }
987              else if (Index <= 0 ||Index > NbrSol) { Standard_OutOfRange::Raise(); }
988              else {
989                if (TheSame2(Index) == 0) {
990                  ParSol = par2sol(Index);
991                  ParArg = pararg2(Index);
992                  PntSol = gp_Pnt2d(pnttg2sol(Index));
993                }
994                else { StdFail_NotDone::Raise(); }
995              }
996 }
997
998 Standard_Boolean GccGeo_Circ2d2TanRad::
999 IsTheSame1 (const Standard_Integer Index) const
1000 {
1001   if (!WellDone) { StdFail_NotDone::Raise(); }
1002   if (Index <= 0 ||Index > NbrSol) { Standard_OutOfRange::Raise(); }
1003
1004   if (TheSame1(Index) == 0) { return Standard_False; }
1005   return Standard_True;
1006 }
1007
1008 Standard_Boolean GccGeo_Circ2d2TanRad::
1009 IsTheSame2 (const Standard_Integer Index) const
1010 {
1011   if (!WellDone) { StdFail_NotDone::Raise(); }
1012   if (Index <= 0 ||Index > NbrSol) { Standard_OutOfRange::Raise(); }
1013
1014   if (TheSame2(Index) == 0) { return Standard_False; }
1015   return Standard_True;
1016 }
1017