0022922: Clean up warnings on uninitialized / unused variables
[occt.git] / src / BRepBlend / BRepBlend_SurfRstEvolRad.cxx
1 // File:        BRepBlend_SurfRstEvolRad.cxx
2 // Created:     Mon Jul 28 17:50:33 1997
3 // Author:      Jerome LEMONIER
4 //              <jlr@sgi64>
5
6
7 #include <BRepBlend_SurfRstEvolRad.ixx>
8 #include <math_Gauss.hxx>
9 #include <math_SVD.hxx>
10
11 #include <ElCLib.hxx>
12 #include <gp.hxx>
13 #include <BlendFunc.hxx>
14 #include <GeomFill.hxx>
15 #include <TColStd_SequenceOfReal.hxx>
16 #include <Standard_DomainError.hxx>
17 #include <Standard_NotImplemented.hxx>
18 #include <Precision.hxx>
19
20 #define Eps 1.e-15
21
22 static void t3dto2d(Standard_Real& a,
23                     Standard_Real& b,
24                     const gp_Vec& A,
25                     const gp_Vec& B,
26                     const gp_Vec& C)
27 {
28   Standard_Real AB = A.Dot(B);
29   Standard_Real AC = A.Dot(C);
30   Standard_Real BC = B.Dot(C);
31   Standard_Real BB = B.Dot(B);
32   Standard_Real CC = C.Dot(C);
33   Standard_Real deno = (BB*CC-BC*BC);
34   a = (AB*CC-AC*BC)/deno;
35   b = (AC*BB-AB*BC)/deno;
36 }
37
38 static void FusionneIntervalles(const TColStd_Array1OfReal& I1,
39                          const TColStd_Array1OfReal& I2,
40                          TColStd_SequenceOfReal& Seq)
41 {
42   Standard_Integer ind1=1, ind2=1;
43   Standard_Real    Epspar = Precision::PConfusion()*0.99;
44   // it is supposed that positioning works with PConfusion()/2
45   Standard_Real    v1, v2;
46 // Initialisation : IND1 and IND2 point at the first element
47 // of each of 2 tables to be processed. INDS points at the last
48 // element created by TABSOR
49
50
51 //--- TABSOR is filled by parsing TABLE1 and TABLE2 simultaneously ---
52 //------------------ and eliminating multiple occurrencies ------------
53
54  while ((ind1<=I1.Upper()) && (ind2<=I2.Upper())) {
55       v1 = I1(ind1);
56       v2 = I2(ind2);
57       if (Abs(v1-v2)<= Epspar) {
58 // Here the elements of I1 and I2 fit.
59          Seq.Append((v1+v2)/2);
60          ind1++;
61          ind2++;
62        }
63       else if (v1 < v2) {
64         // Here the element of I1 fits.
65          Seq.Append(v1);
66          ind1++;
67        }
68       else {
69 // Here the element of TABLE2 fits.
70          Seq.Append(v2);
71          ind2++;
72        }
73     }
74
75   if (ind1>I1.Upper()) { 
76 //----- Here I1 is exhausted, completed with the end of TABLE2 -------
77
78     for (; ind2<=I2.Upper(); ind2++) {
79       Seq.Append(I2(ind2));
80     }
81   }
82
83   if (ind2>I2.Upper()) { 
84 //----- Here I2 is exhausted, completed with the end of I1 -------
85
86     for (; ind1<=I1.Upper(); ind1++) {
87       Seq.Append(I1(ind1));
88     }
89   }
90 }
91
92 //=======================================================================
93 //function : BRepBlend_SurfRstEvolRad
94 //purpose  : Contructor
95 //=======================================================================
96 BRepBlend_SurfRstEvolRad::BRepBlend_SurfRstEvolRad
97 (const Handle(Adaptor3d_HSurface)& Surf,
98  const Handle(Adaptor3d_HSurface)& SurfRst,
99  const Handle(Adaptor2d_HCurve2d)& Rst,
100  const Handle(Adaptor3d_HCurve)& CGuide,
101  const Handle(Law_Function)& Evol):
102  surf(Surf), surfrst(SurfRst), 
103  rst(Rst), cons(Rst,SurfRst), 
104  guide(CGuide), tguide(CGuide),
105  istangent(Standard_True), 
106  maxang(RealFirst()), minang(RealLast()),
107  distmin(RealLast()),
108  mySShape(BlendFunc_Rational)
109
110 { tevol=Evol;
111   fevol=Evol;
112 }
113
114 //=======================================================================
115 //function : 
116 //purpose  : 
117 //=======================================================================
118 Standard_Integer BRepBlend_SurfRstEvolRad::NbVariables() const
119 {  
120   return 3;
121 }
122
123 //=======================================================================
124 //function : 
125 //purpose  : 
126 //=======================================================================
127 Standard_Integer BRepBlend_SurfRstEvolRad::NbEquations() const
128 {
129   return 3;
130 }
131
132 //=======================================================================
133 //function : 
134 //purpose  : 
135 //=======================================================================
136 Standard_Boolean BRepBlend_SurfRstEvolRad::Value
137 (const math_Vector& X,
138  math_Vector& F) 
139 {
140   gp_Vec d1u1,d1v1,ns,vref;
141   Standard_Real norm;
142
143   surf->D1(X(1),X(2),pts,d1u1,d1v1);
144   ptrst = cons.Value(X(3));
145   
146   F(1) = nplan.XYZ().Dot(pts.XYZ()) + theD;
147   
148   F(2) = nplan.XYZ().Dot(ptrst.XYZ()) + theD;
149   
150   ns = d1u1.Crossed(d1v1);
151   norm = nplan.Crossed(ns).Magnitude();
152   ns.SetLinearForm(nplan.Dot(ns)/norm,nplan, -1./norm,ns);
153   vref.SetLinearForm(ray,ns,gp_Vec(ptrst,pts));
154   F(3) = vref.SquareMagnitude() - ray*ray;
155   return Standard_True;
156 }
157
158 //=======================================================================
159 //function : Derivatives
160 //purpose  : 
161 //=======================================================================
162 Standard_Boolean BRepBlend_SurfRstEvolRad::Derivatives
163 (const math_Vector& X,
164  math_Matrix& D) 
165 {  gp_Vec d1u1,d1v1,d2u1,d2v1,d2uv1,d1;
166   gp_Vec ns,ncrossns,resul,temp, vref;
167   
168   Standard_Real norm,ndotns,grosterme;
169   
170   surf->D2(X(1),X(2),pts,d1u1,d1v1,d2u1,d2v1,d2uv1);
171   cons.D1(X(3),ptrst,d1);
172   
173   D(1,1) = nplan.Dot(d1u1);
174   D(1,2) = nplan.Dot(d1v1);
175   D(1,3) = 0.;
176   
177   D(2,1) = 0.;
178   D(2,2) = 0.;
179   D(2,3) = nplan.Dot(d1);
180   
181   
182   ns = d1u1.Crossed(d1v1);
183   ncrossns = nplan.Crossed(ns);
184   norm = ncrossns.Magnitude();
185   ndotns = nplan.Dot(ns);
186   
187   vref.SetLinearForm(ndotns,nplan,-1.,ns);
188   vref.Divide(norm);
189   vref.SetLinearForm(ray,vref,gp_Vec(ptrst,pts));
190   
191   // Derivative corresponding to u1
192   temp = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
193   grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
194   resul.SetLinearForm(-ray/norm*(grosterme*ndotns-nplan.Dot(temp)),nplan,
195                       ray*grosterme/norm,ns,
196                       -ray/norm,temp,
197                       d1u1);
198   
199   D(3,1) = 2.*(resul.Dot(vref));
200   
201   
202   // Derivative corresponding to v1
203   temp = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
204   grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
205   resul.SetLinearForm(-ray/norm*(grosterme*ndotns-nplan.Dot(temp)),nplan,
206                       ray*grosterme/norm,ns,
207                       -ray/norm,temp,
208                       d1v1);
209   
210   D(3,2) = 2.*(resul.Dot(vref));
211   
212   D(3,3) = -2.*(d1.Dot(vref));
213   
214   return Standard_True;
215
216 }
217
218 //=======================================================================
219 //function : 
220 //purpose  : 
221 //=======================================================================
222 Standard_Boolean BRepBlend_SurfRstEvolRad::Values
223 (const math_Vector& X,
224  math_Vector& F,
225  math_Matrix& D) 
226 {
227   gp_Vec d1u1,d1v1,d1;
228   gp_Vec d2u1,d2v1,d2uv1;
229   gp_Vec ns,ncrossns,resul,temp,vref;
230   
231   Standard_Real norm,ndotns,grosterme;
232
233   surf->D2(X(1),X(2),pts,d1u1,d1v1,d2u1,d2v1,d2uv1);
234   cons.D1(X(3),ptrst,d1);
235   
236   F(1) = nplan.XYZ().Dot(pts.XYZ()) + theD;
237   F(2) = nplan.XYZ().Dot(ptrst.XYZ()) + theD;
238   
239   D(1,1) = nplan.Dot(d1u1);
240   D(1,2) = nplan.Dot(d1v1);
241   D(1,3) = 0.;
242   
243   D(2,1) = 0.;
244   D(2,2) = 0.;
245   D(2,3) = nplan.Dot(d1);
246   
247   
248   ns = d1u1.Crossed(d1v1);
249   ncrossns = nplan.Crossed(ns);
250   norm = ncrossns.Magnitude();
251   ndotns = nplan.Dot(ns);
252   
253   vref.SetLinearForm(ndotns,nplan,-1.,ns);
254   vref.Divide(norm);
255   vref.SetLinearForm(ray,vref,gp_Vec(ptrst,pts));
256   
257   F(3) = vref.SquareMagnitude() - ray*ray;
258   
259   
260   // Derivative corresponding to u1
261   temp = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
262   grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
263   resul.SetLinearForm(-ray/norm*(grosterme*ndotns-nplan.Dot(temp)),nplan,
264                       ray*grosterme/norm,ns,
265                       -ray/norm,temp,
266                       d1u1);
267   
268   D(3,1) = 2.*(resul.Dot(vref));
269   
270   
271   // Derivative corresponding to v1
272   temp = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
273   grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
274   resul.SetLinearForm(-ray/norm*(grosterme*ndotns-nplan.Dot(temp)),nplan,
275                       ray*grosterme/norm,ns,
276                       -ray/norm,temp,
277                       d1v1);
278   
279   D(3,2) = 2.*(resul.Dot(vref));
280   
281   D(3,3) = -2.*(d1.Dot(vref));
282   
283   return Standard_True;
284 }
285
286 //=======================================================================
287 //function : 
288 //purpose  : 
289 //=======================================================================
290 void BRepBlend_SurfRstEvolRad::Set
291 (const Handle(Adaptor3d_HSurface)& SurfRef,
292 const Handle(Adaptor2d_HCurve2d)& RstRef) 
293 {
294   surfref = SurfRef;
295   rstref = RstRef;
296 }
297
298 //=======================================================================
299 //function : 
300 //purpose  : 
301 //=======================================================================
302 void BRepBlend_SurfRstEvolRad::Set(const Standard_Real Param) 
303 {
304   d1gui = gp_Vec(0.,0.,0.);
305   nplan = gp_Vec(0.,0.,0.);
306   tguide->D2(Param,ptgui,d1gui,d2gui);
307   normtg = d1gui.Magnitude();
308   nplan.SetXYZ(d1gui.Normalized().XYZ());
309   gp_XYZ nplanXYZ(nplan.XYZ());
310   gp_XYZ ptguiXYZ(ptgui.XYZ());
311   theD =  nplanXYZ.Dot(ptguiXYZ)  ;
312   theD = theD  * (-1.) ;
313   tevol->D1(Param,ray,dray);
314   ray=sg1*ray;
315   dray=sg1*dray;
316 }
317
318 //=======================================================================
319 //function : 
320 //purpose  : Segments the curve in its useful part.
321 //           Precision is taken arbitrary small !?
322 //=======================================================================
323  void BRepBlend_SurfRstEvolRad::Set
324 (const Standard_Real First,
325 const Standard_Real Last) 
326
327   tguide = guide->Trim(First,Last,1.e-12);
328   tevol  = fevol->Trim(First,Last,1.e-12);
329 }
330
331 //=======================================================================
332 //function : 
333 //purpose  : 
334 //=======================================================================
335  void BRepBlend_SurfRstEvolRad::GetTolerance
336 (math_Vector& Tolerance,
337 const Standard_Real Tol) const
338 {
339   Tolerance(1) = surf->UResolution(Tol);
340   Tolerance(2) = surf->VResolution(Tol);
341   Tolerance(3) = cons.Resolution(Tol);
342 }
343
344 //=======================================================================
345 //function : 
346 //purpose  : 
347 //=======================================================================
348  void BRepBlend_SurfRstEvolRad::GetBounds
349 (math_Vector& InfBound,
350 math_Vector& SupBound) const
351 {
352   InfBound(1) = surf->FirstUParameter();
353   InfBound(2) = surf->FirstVParameter();
354   InfBound(3) = cons.FirstParameter();
355   SupBound(1) = surf->LastUParameter();
356   SupBound(2) = surf->LastVParameter();
357   SupBound(3) = cons.LastParameter();
358   
359   if(!Precision::IsInfinite(InfBound(1)) &&
360      !Precision::IsInfinite(SupBound(1))) {
361     Standard_Real range = (SupBound(1) - InfBound(1));
362     InfBound(1) -= range;
363     SupBound(1) += range;
364   }
365   if(!Precision::IsInfinite(InfBound(2)) &&
366      !Precision::IsInfinite(SupBound(2))) {
367     Standard_Real range = (SupBound(2) - InfBound(2));
368     InfBound(2) -= range;
369     SupBound(2) += range;
370   }
371 }
372
373 //=======================================================================
374 //function : 
375 //purpose  : 
376 //=======================================================================
377  Standard_Boolean BRepBlend_SurfRstEvolRad::IsSolution
378 (const math_Vector& Sol,
379 const Standard_Real Tol) 
380 {
381   math_Vector valsol(1,3),secmember(1,3);
382   math_Matrix gradsol(1,3,1,3);
383   
384   gp_Vec dnplan,d1u1,d1v1,d1urst,d1vrst,d1,temp,ns,ns2,ncrossns,resul;
385   gp_Pnt bid;
386   Standard_Real norm,ndotns,grosterme;
387   Standard_Real Cosa,Sina,Angle;
388   
389   Values(Sol,valsol,gradsol);
390   if (Abs(valsol(1)) <= Tol &&
391       Abs(valsol(2)) <= Tol &&
392       Abs(valsol(3)) <= 2*Tol*Abs(ray) ) {
393     
394     // Calculation of tangents
395     
396     pt2ds  = gp_Pnt2d(Sol(1),Sol(2));
397     prmrst = Sol(3);
398     pt2drst = rst->Value(prmrst);
399     surf->D1(Sol(1),Sol(2),pts,d1u1,d1v1);
400     cons.D1(Sol(3),ptrst,d1);
401     dnplan.SetLinearForm(1./normtg,d2gui,
402                          -1./normtg*(nplan.Dot(d2gui)),nplan);
403     
404     temp.SetXYZ(pts.XYZ() - ptgui.XYZ());
405     secmember(1) = normtg - dnplan.Dot(temp);
406     
407     temp.SetXYZ(ptrst.XYZ() - ptgui.XYZ());
408     secmember(2) = normtg - dnplan.Dot(temp);
409     
410     ns = d1u1.Crossed(d1v1);
411     ncrossns = nplan.Crossed(ns);
412     ndotns = nplan.Dot(ns);
413     norm = ncrossns.Magnitude();
414     
415     grosterme = ncrossns.Dot(dnplan.Crossed(ns))/norm/norm;
416     gp_Vec dnw;
417     dnw.SetLinearForm((dnplan.Dot(ns)-grosterme*ndotns)/norm,nplan,
418                       ndotns/norm,dnplan,
419                       grosterme/norm,ns);
420     
421     ns.SetLinearForm(ndotns/norm,nplan, -1./norm,ns);
422     resul.SetLinearForm(ray, ns, gp_Vec(ptrst,pts));
423     
424     secmember(3) = -2.*ray*(dnw.Dot(resul)) - 2.*dray*(ns.Dot(resul)) + 2.*ray*dray;
425     math_Gauss Resol(gradsol);
426     if (Resol.IsDone()) {
427       Resol.Solve(secmember);
428       istangent = Standard_False;
429     }
430     else {
431       math_SVD SingRS (gradsol);
432       if (SingRS.IsDone()) {
433         math_Vector DEDT(1,3);
434         DEDT = secmember;
435         SingRS.Solve(DEDT, secmember, 1.e-6);
436         istangent = Standard_False;
437       }
438       else istangent = Standard_True;
439     }      
440     
441     if (!istangent) {  
442       tgs.SetLinearForm(secmember(1),d1u1,secmember(2),d1v1);
443       tgrst = secmember(3)*d1;
444       tg2ds.SetCoord(secmember(1),secmember(2));
445       surfrst->D1(pt2drst.X(),pt2drst.Y(),bid,d1urst,d1vrst);
446       Standard_Real a, b;
447       t3dto2d(a,b,tgrst,d1urst,d1vrst);
448       tg2drst.SetCoord(a,b);
449       istangent = Standard_False;
450     }
451     else {
452       istangent = Standard_True;
453     }
454     // update of maxang
455     if(ray>0.) ns.Reverse();
456     ns2 = -resul.Normalized();
457     
458     Cosa = ns.Dot(ns2);
459     Sina = nplan.Dot(ns.Crossed(ns2));
460     if (choix%2 != 0) {
461       Sina = -Sina;  //nplan is changed into -nplan
462     }
463     
464     Angle = ACos(Cosa);
465     if (Sina <0.) {
466       Angle = 2.*M_PI - Angle;
467     }
468     
469     if (Angle>maxang) {maxang = Angle;}
470     if (Angle<minang) {minang = Angle;}
471     distmin = Min( distmin, pts.Distance(ptrst));
472     
473     return Standard_True;
474   }
475   istangent = Standard_True;
476   return Standard_False;
477 }
478
479
480 //=======================================================================
481 //function : GetMinimalDistance
482 //purpose  : 
483 //=======================================================================
484
485 Standard_Real BRepBlend_SurfRstEvolRad::GetMinimalDistance() const
486 {
487   return distmin;
488 }
489
490 //=======================================================================
491 //function : 
492 //purpose  : 
493 //=======================================================================
494 const gp_Pnt& BRepBlend_SurfRstEvolRad::PointOnS() const
495 {
496   return pts;
497 }
498
499 //=======================================================================
500 //function : 
501 //purpose  : 
502 //=======================================================================
503 const gp_Pnt& BRepBlend_SurfRstEvolRad::PointOnRst() const
504 {
505   return ptrst ;
506 }
507
508 //=======================================================================
509 //function : 
510 //purpose  : 
511 //=======================================================================
512 const gp_Pnt2d& BRepBlend_SurfRstEvolRad::Pnt2dOnS() const
513 {
514   return pt2ds;
515 }
516
517 //=======================================================================
518 //function : 
519 //purpose  : 
520 //=======================================================================
521 const gp_Pnt2d& BRepBlend_SurfRstEvolRad::Pnt2dOnRst() const
522 {
523   return pt2drst;
524 }
525
526 //=======================================================================
527 //function : 
528 //purpose  : 
529 //=======================================================================
530  Standard_Real BRepBlend_SurfRstEvolRad::ParameterOnRst() const
531 {
532   return prmrst;
533 }
534
535 //=======================================================================
536 //function : 
537 //purpose  : 
538 //=======================================================================
539  Standard_Boolean BRepBlend_SurfRstEvolRad::IsTangencyPoint() const
540 {
541   return istangent;
542 }
543
544 //=======================================================================
545 //function : 
546 //purpose  : 
547 //=======================================================================
548 const gp_Vec& BRepBlend_SurfRstEvolRad::TangentOnS() const
549 {
550   if (istangent) {Standard_DomainError::Raise();}
551   return tgs;
552 }
553
554 //=======================================================================
555 //function : 
556 //purpose  : 
557 //=======================================================================
558 const gp_Vec2d& BRepBlend_SurfRstEvolRad::Tangent2dOnS() const
559 {
560   if (istangent) {Standard_DomainError::Raise();}
561   return tg2ds;
562 }
563
564 //=======================================================================
565 //function : 
566 //purpose  : 
567 //=======================================================================
568 const gp_Vec& BRepBlend_SurfRstEvolRad::TangentOnRst() const
569 {
570   if (istangent) {Standard_DomainError::Raise();}
571   return tgrst;
572 }
573
574 //=======================================================================
575 //function : 
576 //purpose  : 
577 //=======================================================================
578 const gp_Vec2d& BRepBlend_SurfRstEvolRad::Tangent2dOnRst() const
579 {
580   if (istangent) {Standard_DomainError::Raise();}
581   return tg2drst;
582 }
583
584 //=======================================================================
585 //function : 
586 //purpose  : 
587 //=======================================================================
588  Standard_Boolean BRepBlend_SurfRstEvolRad::Decroch
589 (const math_Vector& Sol,
590 gp_Vec& NS,
591 gp_Vec& TgS) const
592 {
593   gp_Vec TgRst, NRst, NRstInPlane, NSInPlane;
594   gp_Pnt bid,Center;
595   gp_Vec d1u,d1v;
596   Standard_Real norm,unsurnorm;
597   
598   surf->D1(Sol(1),Sol(2),bid,d1u,d1v);
599   NS = NSInPlane = d1u.Crossed(d1v);
600   
601   norm = nplan.Crossed(NS).Magnitude();
602   unsurnorm = 1./norm;
603   NSInPlane.SetLinearForm(nplan.Dot(NS)*unsurnorm,nplan,-unsurnorm,NS);
604   
605   Center.SetXYZ(bid.XYZ()+ray*NSInPlane.XYZ());
606   if(choix>2) NSInPlane.Reverse();
607   TgS = nplan.Crossed(gp_Vec(Center,bid));
608   if (choix%2 == 1) {
609     TgS.Reverse();
610   }
611   Standard_Real u,v;
612   rstref->Value(Sol(3)).Coord(u,v);
613   surfref->D1(u,v,bid,d1u,d1v);
614   NRst = d1u.Crossed(d1v);
615   norm = nplan.Crossed(NRst).Magnitude();
616   unsurnorm = 1./norm;
617   NRstInPlane.SetLinearForm(nplan.Dot(NRst)*unsurnorm,nplan,-unsurnorm,NRst);
618   gp_Vec centptrst(Center,bid);
619   if(centptrst.Dot(NRstInPlane) < 0.) NRstInPlane.Reverse();
620   TgRst = nplan.Crossed(centptrst);
621   if (choix%2 == 1) {
622     TgRst.Reverse();
623   }
624
625   Standard_Real dot, NT = NRstInPlane.Magnitude();
626   NT *= TgRst.Magnitude();
627   if (Abs(NT) < 1.e-7) {
628     return Standard_False; // Singularity or Incoherence.
629   }
630   dot = NRstInPlane.Dot(TgRst);
631   dot /= NT;
632          
633   return (dot < 1.e-10);
634 }
635
636 //=======================================================================
637 //function : 
638 //purpose  : 
639 //=======================================================================
640  void BRepBlend_SurfRstEvolRad::Set (const Standard_Integer Choix) 
641 {
642   choix = Choix;
643   switch (choix) {
644   case 1 :
645   case 2 :
646     sg1 = -1;
647     break;
648   case 3 :
649   case 4 :
650     sg1 = 1;
651     break;
652   default :
653     sg1 = -1;
654     break;
655   }
656 }
657
658 //=======================================================================
659 //function : 
660 //purpose  : 
661 //=======================================================================
662  void BRepBlend_SurfRstEvolRad::Set(const BlendFunc_SectionShape TypeSection) 
663 {
664   mySShape = TypeSection;
665 }
666
667 //=======================================================================
668 //function : 
669 //purpose  : 
670 //=======================================================================
671  void BRepBlend_SurfRstEvolRad::Section
672 (const Standard_Real Param,
673 const Standard_Real U,
674 const Standard_Real V,
675 const Standard_Real W,
676 Standard_Real& Pdeb,
677 Standard_Real& Pfin,
678 gp_Circ& C) 
679 {
680   gp_Vec d1u1,d1v1;
681   gp_Vec ns, np;
682   Standard_Real norm;
683   gp_Pnt Center;
684   
685   tguide->D1(Param,ptgui,d1gui);
686   np  = d1gui.Normalized();
687   ray = sg1*tevol->Value(Param);
688   
689   surf->D1(U,V,pts,d1u1,d1v1);
690   ptrst = cons.Value(W);
691   
692   ns = d1u1.Crossed(d1v1);
693   
694   norm = nplan.Crossed(ns).Magnitude();
695   ns.SetLinearForm(nplan.Dot(ns)/norm,nplan, -1./norm,ns);
696   Center.SetXYZ(pts.XYZ()+ray*ns.XYZ());
697   C.SetRadius(Abs(ray));
698
699   if (ray > 0) {
700     ns.Reverse();
701   }    
702   if (choix%2 != 0) {
703     np.Reverse();
704   }
705   C.SetPosition(gp_Ax2(Center,np,ns));
706
707   Pdeb = 0.; //ElCLib::Parameter(C,pts);
708   Pfin = ElCLib::Parameter(C,ptrst);
709
710   // Test negative and almost null angles : Single Case
711   if (Pfin>1.5*M_PI) {
712     np.Reverse();
713     C.SetPosition(gp_Ax2(Center,np,ns));
714     Pfin = ElCLib::Parameter(C,ptrst);
715   }
716   if (Pfin < Precision::PConfusion()) Pfin += Precision::PConfusion();  
717 }
718
719 //=======================================================================
720 //function : 
721 //purpose  : 
722 //=======================================================================
723  Standard_Boolean BRepBlend_SurfRstEvolRad::IsRational() const
724 {
725   return  (mySShape==BlendFunc_Rational || mySShape==BlendFunc_QuasiAngular);
726 }
727
728 //=======================================================================
729 //function : 
730 //purpose  : 
731 //=======================================================================
732  Standard_Real BRepBlend_SurfRstEvolRad::GetSectionSize() const
733 {
734   return maxang*Abs(ray);
735 }
736
737 //=======================================================================
738 //function : 
739 //purpose  : 
740 //=======================================================================
741  void BRepBlend_SurfRstEvolRad::GetMinimalWeight(TColStd_Array1OfReal& Weigths) const
742 {
743   BlendFunc::GetMinimalWeights(mySShape, myTConv, minang, maxang, Weigths );
744   // It is supposed that it does not depend on the Radius! 
745 }
746
747 //=======================================================================
748 //function : 
749 //purpose  : 
750 //=======================================================================
751  Standard_Integer BRepBlend_SurfRstEvolRad::NbIntervals(const GeomAbs_Shape S) const
752 {
753   Standard_Integer Nb_Int_Courbe, Nb_Int_Loi;
754   Nb_Int_Courbe =  guide->NbIntervals(BlendFunc::NextShape(S));
755   Nb_Int_Loi    =  fevol->NbIntervals(S);
756
757   if  (Nb_Int_Loi==1) {
758     return Nb_Int_Courbe;
759   }
760
761   TColStd_Array1OfReal IntC(1, Nb_Int_Courbe+1);
762   TColStd_Array1OfReal IntL(1, Nb_Int_Loi+1);
763   TColStd_SequenceOfReal    Inter;
764   guide->Intervals(IntC, BlendFunc::NextShape(S));
765   fevol->Intervals(IntL, S);
766
767   FusionneIntervalles( IntC, IntL, Inter);
768   return Inter.Length()-1;
769 }
770
771 //=======================================================================
772 //function : 
773 //purpose  : 
774 //=======================================================================
775  void BRepBlend_SurfRstEvolRad::Intervals(TColStd_Array1OfReal& T,
776                                           const GeomAbs_Shape S) const
777 {
778   Standard_Integer Nb_Int_Courbe, Nb_Int_Loi;  
779   Nb_Int_Courbe =  guide->NbIntervals(BlendFunc::NextShape(S));
780   Nb_Int_Loi    =  fevol->NbIntervals(S);
781
782   if  (Nb_Int_Loi==1) {
783     guide->Intervals(T, BlendFunc::NextShape(S));
784   }
785   else {
786     TColStd_Array1OfReal IntC(1, Nb_Int_Courbe+1);
787     TColStd_Array1OfReal IntL(1, Nb_Int_Loi+1);
788     TColStd_SequenceOfReal    Inter;
789     guide->Intervals(IntC, BlendFunc::NextShape(S));
790     fevol->Intervals(IntL, S);
791
792     FusionneIntervalles( IntC, IntL, Inter);
793     for (Standard_Integer ii=1; ii<=Inter.Length(); ii++) {
794       T(ii) = Inter(ii);
795     }
796   } 
797 }
798
799
800 //=======================================================================
801 //function : 
802 //purpose  : 
803 //=======================================================================
804  void BRepBlend_SurfRstEvolRad::GetShape
805 (Standard_Integer& NbPoles,
806 Standard_Integer& NbKnots,
807 Standard_Integer& Degree,
808 Standard_Integer& NbPoles2d) 
809 {
810   NbPoles2d = 2;
811   BlendFunc::GetShape(mySShape,maxang,NbPoles,NbKnots,Degree,myTConv);
812 }
813
814 //=======================================================================
815 //function : 
816 //purpose  : 
817 //=======================================================================
818  void BRepBlend_SurfRstEvolRad::GetTolerance
819 (const Standard_Real BoundTol,
820 const Standard_Real SurfTol,
821 const Standard_Real AngleTol,
822 math_Vector& Tol3d,
823 math_Vector& Tol1d) const
824 {
825   Standard_Integer low = Tol3d.Lower() , up=Tol3d.Upper();
826   Standard_Real Tol;
827   Tol= GeomFill::GetTolerance(myTConv, minang, Abs(ray), 
828                                AngleTol, SurfTol);
829   Tol1d.Init(SurfTol);
830   Tol3d.Init(SurfTol);
831   Tol3d(low+1) = Tol3d(up-1) = Min( Tol, SurfTol);
832   Tol3d(low) = Tol3d(up) = Min( Tol, BoundTol);
833 }
834
835 //=======================================================================
836 //function : 
837 //purpose  : 
838 //=======================================================================
839  void BRepBlend_SurfRstEvolRad::Knots(TColStd_Array1OfReal& TKnots) 
840 {
841   GeomFill::Knots(myTConv,TKnots);
842 }
843
844 //=======================================================================
845 //function : 
846 //purpose  : 
847 //=======================================================================
848  void BRepBlend_SurfRstEvolRad::Mults(TColStd_Array1OfInteger& TMults) 
849 {
850   GeomFill::Mults(myTConv,TMults);
851 }
852
853 //=======================================================================
854 //function : 
855 //purpose  : 
856 //=======================================================================
857  Standard_Boolean BRepBlend_SurfRstEvolRad::Section
858 (const Blend_Point& P,
859 TColgp_Array1OfPnt& Poles,
860 TColgp_Array1OfVec& DPoles,
861 TColgp_Array1OfPnt2d& Poles2d,
862 TColgp_Array1OfVec2d& DPoles2d,
863 TColStd_Array1OfReal& Weigths,
864 TColStd_Array1OfReal& DWeigths) 
865 {
866   
867   gp_Vec d1u1,d1v1,d2u1,d2v1,d2uv1,d1;
868   gp_Vec ns,ns2,dnplan,dnw,dn2w;//,np2,dnp2;
869   gp_Vec ncrossns;;
870   gp_Vec resulu,resulv,temp,tgct,resul;
871   gp_Vec d1urst,d1vrst;
872   gp_Pnt Center,bid;
873   
874   Standard_Real norm,ndotns,grosterme,dray;
875   
876   math_Vector sol(1,3),valsol(1,3),secmember(1,3);
877   math_Matrix gradsol(1,3,1,3);
878   
879   Standard_Real prm = P.Parameter(),rayprim;
880   Standard_Integer low = Poles.Lower();
881   Standard_Integer upp = Poles.Upper();
882   Standard_Boolean istgt;
883   
884   tguide->D2(prm,ptgui,d1gui,d2gui);
885
886   tevol->D1(prm,ray,dray);
887   ray=sg1*ray;
888   dray=sg1*dray;
889   normtg = d1gui.Magnitude();
890   nplan  = d1gui.Normalized();
891   dnplan.SetLinearForm(1./normtg,d2gui,
892                        -1./normtg*(nplan.Dot(d2gui)),nplan);
893   
894   P.ParametersOnS(sol(1),sol(2));
895   sol(3) = prmrst = P.ParameterOnC();
896   pt2drst = rst->Value(prmrst);
897   
898   Values(sol,valsol,gradsol);
899   
900   surf->D2(sol(1),sol(2),pts,d1u1,d1v1,d2u1,d2v1,d2uv1);
901   cons.D1(sol(3),ptrst,d1);
902   
903   temp.SetXYZ(pts.XYZ()- ptgui.XYZ());
904   secmember(1) = normtg - dnplan.Dot(temp);
905   
906   temp.SetXYZ(ptrst.XYZ()- ptgui.XYZ());
907   secmember(2) = normtg - dnplan.Dot(temp);
908   
909   ns = d1u1.Crossed(d1v1);
910   ncrossns = nplan.Crossed(ns);
911   ndotns = nplan.Dot(ns);
912   norm = ncrossns.Magnitude();
913   if (norm < Eps)  {
914     norm = 1; // Not enough, but it is not necessary to stop
915 #if DEB
916     cout << " SurfRstEvolRad : Surface single " << endl;
917 #endif
918   }
919   
920   // Derivative of n1 corresponding to w
921   
922   grosterme = ncrossns.Dot(dnplan.Crossed(ns))/norm/norm;
923   dnw.SetLinearForm((dnplan.Dot(ns)-grosterme*ndotns)/norm,nplan,
924                     ndotns/norm,dnplan,
925                     grosterme/norm,ns);
926   
927   temp.SetLinearForm(ndotns/norm,nplan, -1./norm,ns);
928   resul.SetLinearForm(ray,temp,gp_Vec(ptrst,pts));
929
930   //secmember(3) = -2.*ray*(dnw.Dot(resul)); // jag 950105 il manquait ray
931   secmember(3) = -2.*ray*(dnw.Dot(resul)) - 2.*dray*(temp.Dot(resul)) + 2.*ray*dray;
932   math_Gauss Resol(gradsol);
933   
934   if (Resol.IsDone()) {
935     Resol.Solve(secmember);
936     istgt = Standard_False;
937   }
938   else {    
939     math_SVD SingRS (gradsol);
940     if (SingRS.IsDone()) {
941       math_Vector DEDT(1,3);
942       DEDT = secmember;
943       SingRS.Solve(DEDT, secmember, 1.e-6);
944       istgt = Standard_False;
945     }
946     else istgt = Standard_True;
947   }
948
949   if (!istgt) {
950
951     tgs.SetLinearForm(secmember(1),d1u1,secmember(2),d1v1);
952     tgrst = secmember(3)*d1;
953     // Derivative of n1 corresponding to u1
954     temp = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
955     grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
956     resulu.SetLinearForm(-(grosterme*ndotns-nplan.Dot(temp))/norm,nplan,
957                          grosterme/norm,ns,
958                          -1./norm,temp);
959     
960     // Derivative of n1 corresponding to v1
961     temp = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
962     grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
963     resulv.SetLinearForm(-(grosterme*ndotns-nplan.Dot(temp))/norm,nplan,
964                          grosterme/norm,ns,
965                          -1./norm,temp);
966     
967     
968     dnw.SetLinearForm(secmember(1),resulu,secmember(2),resulv,dnw);
969     ns.SetLinearForm(ndotns/norm,nplan, -1./norm,ns);
970     
971     dn2w.SetLinearForm(ray, dnw, -1., tgrst, tgs);
972     dn2w.SetLinearForm(dray,ns,dn2w);
973     norm = resul.Magnitude();
974     dn2w.Divide(norm);
975     ns2 = -resul.Normalized();
976     dn2w.SetLinearForm(ns2.Dot(dn2w),ns2,-1.,dn2w);
977     
978     istgt = Standard_False;
979   }
980   else {
981     ns.SetLinearForm(ndotns/norm,nplan, -1./norm,ns);
982     ns2 = -resul.Normalized();
983     istgt = Standard_True;
984   }
985   
986   // Tops 2D
987   
988   Poles2d(Poles2d.Lower()).SetCoord(sol(1),sol(2));
989   Poles2d(Poles2d.Upper()).SetCoord(pt2drst.X(),pt2drst.Y());
990   if (!istgt) {
991     DPoles2d(Poles2d.Lower()).SetCoord(secmember(1),secmember(2));
992     surfrst->D1(pt2drst.X(),pt2drst.Y(),bid,d1urst,d1vrst);
993     Standard_Real a, b;
994     t3dto2d(a,b,tgrst,d1urst,d1vrst);
995     DPoles2d(Poles2d.Upper()).SetCoord(a,b);
996   }
997   
998   // Linear Case
999   if (mySShape == BlendFunc_Linear) {
1000     Poles(low) = pts;
1001     Poles(upp) = ptrst;
1002     Weigths(low) = 1.0;
1003     Weigths(upp) = 1.0;
1004     if (!istgt) {
1005       DPoles(low) = tgs;
1006       DPoles(upp) = tgrst;
1007       DWeigths(low) = 0.0;
1008       DWeigths(upp) = 0.0;
1009     }
1010     return (!istgt);
1011   }
1012   
1013   // Case of the circle
1014   Center.SetXYZ(pts.XYZ()+ray*ns.XYZ());
1015   if (!istgt) {
1016     tgct.SetLinearForm(ray,dnw,dray,ns,tgs);
1017   }
1018   
1019   if (ray > 0.) {
1020     ns.Reverse();
1021     if (!istgt) {
1022       dnw.Reverse();
1023     }
1024   }
1025   if (choix%2 != 0) {
1026     nplan.Reverse();
1027     dnplan.Reverse();
1028   }
1029   if (!istgt) {
1030     if (ray < 0.) { // to avoid Abs(dray) some lines below
1031       rayprim = -dray;
1032     }
1033     else rayprim = dray;
1034
1035     return GeomFill::GetCircle(myTConv, 
1036                                ns, ns2, 
1037                                dnw, dn2w, 
1038                                nplan, dnplan, 
1039                                pts, ptrst, 
1040                                tgs, tgrst, 
1041                                Abs(ray),rayprim , 
1042                                Center, tgct, 
1043                                Poles, 
1044                                DPoles,
1045                                Weigths, 
1046                                DWeigths); 
1047   }
1048   else {
1049     GeomFill::GetCircle(myTConv,
1050                         ns, ns2, 
1051                         nplan, pts, ptrst,
1052                         Abs(ray), Center, 
1053                         Poles, Weigths);
1054     return Standard_False;
1055   }
1056 }
1057
1058 //=======================================================================
1059 //function : 
1060 //purpose  : 
1061 //=======================================================================
1062  Standard_Boolean BRepBlend_SurfRstEvolRad::Section
1063 (const Blend_Point& P,
1064 TColgp_Array1OfPnt& Poles,
1065 TColgp_Array1OfVec& DPoles,
1066 TColgp_Array1OfVec& D2Poles,
1067 TColgp_Array1OfPnt2d& Poles2d,
1068 TColgp_Array1OfVec2d& DPoles2d,
1069 TColgp_Array1OfVec2d& D2Poles2d,
1070 TColStd_Array1OfReal& Weigths,
1071 TColStd_Array1OfReal& DWeigths,
1072 TColStd_Array1OfReal& D2Weigths) 
1073 {
1074     return Standard_False;
1075 }
1076
1077 //=======================================================================
1078 //function : 
1079 //purpose  : 
1080 //=======================================================================
1081  void BRepBlend_SurfRstEvolRad::Section
1082 (const Blend_Point& P,
1083 TColgp_Array1OfPnt& Poles,
1084 TColgp_Array1OfPnt2d& Poles2d,
1085 TColStd_Array1OfReal& Weigths) 
1086 {
1087   gp_Vec d1u1,d1v1;//,d1;
1088   gp_Vec ns,ns2;//,temp,np2;
1089   gp_Pnt Center;
1090   
1091   Standard_Real norm,u1,v1,w;
1092   
1093   Standard_Real prm = P.Parameter();
1094   Standard_Integer low = Poles.Lower();
1095   Standard_Integer upp = Poles.Upper();
1096   
1097   tguide->D1(prm,ptgui,d1gui);
1098   ray = tevol->Value(prm);
1099   ray=sg1*ray;
1100   nplan  = d1gui.Normalized();
1101   
1102   P.ParametersOnS(u1,v1);
1103   w = P.ParameterOnC(); //jlr : point on curve not on surface
1104   gp_Pnt2d  pt2d = rst->Value(w);
1105
1106   surf->D1(u1,v1,pts,d1u1,d1v1);
1107   ptrst = cons.Value(w);
1108
1109   distmin = Min (distmin, pts.Distance(ptrst));
1110   
1111   Poles2d(Poles2d.Lower()).SetCoord(u1,v1);
1112   Poles2d(Poles2d.Upper()).SetCoord(pt2d.X(),pt2d.Y());
1113   
1114   // Linear case
1115   if (mySShape == BlendFunc_Linear) {
1116     Poles(low) = pts;
1117     Poles(upp) = ptrst;
1118     Weigths(low) = 1.0;
1119     Weigths(upp) = 1.0;
1120     return;
1121   }
1122   
1123   ns = d1u1.Crossed(d1v1);
1124   norm = nplan.Crossed(ns).Magnitude();
1125   
1126   ns.SetLinearForm(nplan.Dot(ns)/norm,nplan, -1./norm,ns);
1127   
1128   Center.SetXYZ(pts.XYZ()+ray*ns.XYZ());
1129   
1130   ns2 = gp_Vec(Center,ptrst).Normalized();
1131   if(ray>0) ns.Reverse();
1132   if (choix%2 != 0) {
1133     nplan.Reverse();
1134   }
1135   
1136   GeomFill::GetCircle(myTConv,
1137                       ns, ns2, 
1138                       nplan, pts, ptrst,
1139                       Abs(ray), Center, 
1140                       Poles, Weigths);
1141 }
1142
1143 void BRepBlend_SurfRstEvolRad::Resolution(const Standard_Integer IC2d,
1144                                           const Standard_Real Tol,
1145                                           Standard_Real& TolU,
1146                                           Standard_Real& TolV) const
1147 {
1148   if(IC2d == 1){
1149     TolU = surf->UResolution(Tol);
1150     TolV = surf->VResolution(Tol);
1151   }
1152   else {
1153     TolU = surfrst->UResolution(Tol);
1154     TolV = surfrst->VResolution(Tol);
1155   }
1156 }
1157
1158