0024057: Eliminate compiler warning C4100 in MSVC++ with warning level 4
[occt.git] / src / BlendFunc / BlendFunc_CSCircular.cxx
1 // Created on: 1995-01-04
2 // Created by: Jacques GOUSSARD
3 // Copyright (c) 1995-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
5 //
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
10 //
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13 //
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
20
21
22 // Modified 10/09/1996 PMN Ajout de (Nb)Intervalles, IsRationnal 
23 //                        + Utilisation de GeomFill::GetCircle dans Section.
24 // Modified 23/06/1997 PMN : Pb de division par 0.
25
26 #include <BlendFunc_CSCircular.ixx>
27
28 #include <math_Gauss.hxx>
29
30 #include <ElCLib.hxx>
31 #include <gp.hxx>
32 #include <BlendFunc.hxx>
33 #include <GeomFill.hxx>
34 #include <Standard_NotImplemented.hxx>
35 #include <Standard_DomainError.hxx>
36 #include <Precision.hxx>
37
38 #define Eps 1.e-15
39
40 //=======================================================================
41 //function : BlendFunc_CSCircular
42 //purpose  : 
43 //=======================================================================
44
45 BlendFunc_CSCircular::BlendFunc_CSCircular(const Handle(Adaptor3d_HSurface)& S,
46                                            const Handle(Adaptor3d_HCurve)& C,
47                                            const Handle(Adaptor3d_HCurve)& CGuide,
48                                            const Handle(Law_Function)& L) :
49        surf(S),curv(C),guide(CGuide),istangent(Standard_True),law(L),
50        //prmc, dprmc, istangent, ray, choix, normtg,
51        maxang(RealFirst()),minang(RealLast()),mySShape(BlendFunc_Rational)
52        //myTConv
53 {
54   law = L;
55 }
56
57
58 //=======================================================================
59 //function : NbVariables
60 //purpose  : 
61 //=======================================================================
62
63 Standard_Integer BlendFunc_CSCircular::NbVariables () const
64 {
65   return 2; 
66 }
67
68
69 //=======================================================================
70 //function : NbEquations
71 //purpose  : 
72 //=======================================================================
73
74 Standard_Integer BlendFunc_CSCircular::NbEquations () const
75 {
76   return 2; 
77 }
78
79
80 //=======================================================================
81 //function : Set
82 //purpose  : 
83 //=======================================================================
84
85 void BlendFunc_CSCircular::Set(const Standard_Real Radius, const Standard_Integer Choix)
86 {
87   choix = Choix;
88   switch (Choix)
89   {
90     case 3 :
91     case 4 :
92       ray = Abs(Radius);
93       break;
94     default :
95       ray = -Abs(Radius);
96       break;
97   }
98 }
99
100
101 //=======================================================================
102 //function : Set
103 //purpose  : 
104 //=======================================================================
105
106 void BlendFunc_CSCircular::Set(const BlendFunc_SectionShape TypeSection)
107 {
108   mySShape = TypeSection;
109 }
110
111
112 //=======================================================================
113 //function : Set
114 //purpose  : 
115 //=======================================================================
116
117 void BlendFunc_CSCircular::Set(const Standard_Real Param)
118 {
119   gp_Pnt ptgui;
120   guide->D2(Param,ptgui,d1gui,d2gui);
121   law->D1(Param,prmc,dprmc);
122
123   normtg = d1gui.Magnitude();
124   nplan  = d1gui.Normalized();
125 }
126
127 //=======================================================================
128 //function : Set
129 //purpose  : 
130 //=======================================================================
131
132 void BlendFunc_CSCircular::Set(const Standard_Real, const Standard_Real)
133 {
134   Standard_NotImplemented::Raise("BlendFunc_CSCircular::Set");
135 }
136
137 //=======================================================================
138 //function : GetTolerance
139 //purpose  : 
140 //=======================================================================
141
142 void BlendFunc_CSCircular::GetTolerance(math_Vector& Tolerance, const Standard_Real Tol) const
143 {
144   Tolerance(1) = surf->UResolution(Tol);
145   Tolerance(2) = surf->VResolution(Tol);
146 }
147
148
149 //=======================================================================
150 //function : GetBounds
151 //purpose  : 
152 //=======================================================================
153
154 void BlendFunc_CSCircular::GetBounds(math_Vector& InfBound, math_Vector& SupBound) const
155 {
156   InfBound(1) = surf->FirstUParameter();
157   InfBound(2) = surf->FirstVParameter();
158   SupBound(1) = surf->LastUParameter();
159   SupBound(2) = surf->LastVParameter();
160
161   if(!Precision::IsInfinite(InfBound(1)) &&
162      !Precision::IsInfinite(SupBound(1))) {
163     const Standard_Real range = (SupBound(1) - InfBound(1));
164     InfBound(1) -= range;
165     SupBound(1) += range;
166   }
167   if(!Precision::IsInfinite(InfBound(2)) &&
168      !Precision::IsInfinite(SupBound(2))) {
169     const Standard_Real range = (SupBound(2) - InfBound(2));
170     InfBound(2) -= range;
171     SupBound(2) += range;
172   }
173 }
174
175
176 //=======================================================================
177 //function : IsSolution
178 //purpose  : 
179 //=======================================================================
180
181 Standard_Boolean BlendFunc_CSCircular::IsSolution(const math_Vector& Sol, const Standard_Real Tol)
182 {
183   math_Vector valsol(1,2),secmember(1,2);
184   math_Matrix gradsol(1,2,1,2);
185
186   gp_Vec dnplan,d1u1,d1v1,d1c,d2c,temp,ns,ncrossns,resul,nc;
187   Standard_Real norm,ndotns,grosterme;
188   Standard_Real Cosa,Sina,Angle;
189
190   Values(Sol,valsol,gradsol);
191
192   if (Abs(valsol(1)) <= Tol &&
193       Abs(valsol(2)) <= Tol*Tol) {
194
195     // Calcul des tangentes
196
197     pt2d = gp_Pnt2d(Sol(1),Sol(2));
198
199     surf->D1(Sol(1),Sol(2),pts,d1u1,d1v1);
200     curv->D2(prmc,ptc,d1c,d2c);
201
202     dnplan.SetLinearForm(1./normtg,d2gui,
203                          -1./normtg*(nplan.Dot(d2gui)),nplan);
204
205     ns = d1u1.Crossed(d1v1);
206     ncrossns = nplan.Crossed(ns);
207     ndotns = nplan.Dot(ns);
208     norm = ncrossns.Magnitude();
209     if (norm < Eps) {
210       norm = 1.;
211 //#if DEB
212 //      cout << "CSCircular : Surface singuliere !" << endl;
213 //#endif
214     }
215
216     temp.SetXYZ(pts.XYZ() - ptc.XYZ());
217     secmember(1) = dprmc*(nplan.Dot(d1c)) - dnplan.Dot(temp);    
218
219     grosterme = ncrossns.Dot(dnplan.Crossed(ns))/norm/norm;
220     temp.SetLinearForm(ray/norm*(dnplan.Dot(ns)-grosterme*ndotns),nplan,
221                        ray*ndotns/norm,dnplan,
222                        ray*grosterme/norm,ns);
223     temp -= dprmc*d1c;
224
225     ns.SetLinearForm(ndotns/norm,nplan, -1./norm,ns);
226     resul.SetLinearForm(ray,ns,gp_Vec(ptc,pts));
227
228     secmember(2) = -2.*(resul.Dot(temp));
229
230     math_Gauss Resol(gradsol);
231     if (Resol.IsDone()) {
232
233       Resol.Solve(secmember);
234
235       tgs.SetLinearForm(secmember(1),d1u1,secmember(2),d1v1);
236       tgc = dprmc*d1c;
237       tg2d.SetCoord(secmember(1),secmember(2));
238       istangent = Standard_False;
239     }
240     else {
241       istangent = Standard_True;
242     }
243     // mise a jour de maxang
244
245     if(ray > 0.) ns.Reverse();
246     nc = -resul.Normalized();
247
248     Cosa = ns.Dot(nc);
249     Sina = nplan.Dot(ns.Crossed(nc));
250     if (choix%2 != 0) {
251       Sina = -Sina;  //nplan est change en -nplan
252     }
253
254     Angle = ACos(Cosa);
255     if (Sina <0.) {
256       Angle = 2.*M_PI - Angle;
257     }
258
259     if (Angle>maxang) {maxang = Angle;}
260     if (Angle<minang) {minang = Angle;}
261
262     return Standard_True;
263   }
264   istangent = Standard_True;
265   return Standard_False;
266 }
267
268
269 //=======================================================================
270 //function : Value
271 //purpose  : 
272 //=======================================================================
273
274 Standard_Boolean BlendFunc_CSCircular::Value(const math_Vector& X, math_Vector& F)
275 {
276   gp_Vec d1u1,d1v1,d1c;
277
278   surf->D1(X(1),X(2),pts,d1u1,d1v1);
279   curv->D1(prmc,ptc,d1c);
280
281   F(1) = nplan.XYZ().Dot(pts.XYZ()-ptc.XYZ());
282
283   gp_Vec ns = d1u1.Crossed(d1v1);
284   Standard_Real norm = nplan.Crossed(ns).Magnitude();
285   if (norm < Eps) {
286     norm = 1.;
287 //#if DEB
288 //    cout << "CSCircular : Surface singuliere !" << endl;
289 //#endif
290   }
291
292   ns.SetLinearForm(nplan.Dot(ns)/norm,nplan, -1./norm,ns);
293   gp_Vec vref;
294   vref.SetLinearForm(ray,ns,gp_Vec(ptc,pts));
295
296   F(2) = vref.SquareMagnitude() - ray*ray;
297
298   pt2d = gp_Pnt2d(X(1),X(2));
299   return Standard_True;
300 }
301
302
303 //=======================================================================
304 //function : Derivatives
305 //purpose  : 
306 //=======================================================================
307
308 Standard_Boolean BlendFunc_CSCircular::Derivatives(const math_Vector& X, math_Matrix& D)
309 {
310   gp_Vec d1u1,d1v1,d2u1,d2v1,d2uv1,d1c;
311   gp_Vec ns,ncrossns,resul,temp,nsov,vref;
312
313   Standard_Real norm,ndotns,grosterme;
314
315   surf->D2(X(1),X(2),pts,d1u1,d1v1,d2u1,d2v1,d2uv1);
316   curv->D1(prmc,ptc,d1c);
317
318   D(1,1) = nplan.Dot(d1u1);
319   D(1,2) = nplan.Dot(d1v1);
320
321   ns = d1u1.Crossed(d1v1);
322   ncrossns = nplan.Crossed(ns);
323   norm = ncrossns.Magnitude();
324   if (norm < Eps) {
325     norm = 1.;
326 //#if DEB
327 //    cout << "CSCircular : Surface singuliere !" << endl;
328 //#endif
329   }
330
331   ndotns = nplan.Dot(ns);
332
333   nsov.SetLinearForm(nplan.Dot(ns)/norm,nplan, -1./norm,ns);
334   vref.SetLinearForm(ray,nsov,gp_Vec(ptc,pts));
335
336   // Derivee par rapport a u de Ps + ray*ns
337   temp = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
338   grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
339   resul.SetLinearForm(-ray/norm*(grosterme*ndotns-nplan.Dot(temp)),nplan,
340                       ray*grosterme/norm,ns,
341                       -ray/norm,temp,
342                       d1u1);
343
344   D(2,1) = 2.*(resul.Dot(vref));
345
346   // Derivee par rapport a v
347   temp = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
348   grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
349   resul.SetLinearForm(-ray/norm*(grosterme*ndotns-nplan.Dot(temp)),nplan,
350                       ray*grosterme/norm,ns,
351                       -ray/norm,temp,
352                       d1v1);
353
354   D(2,2) = 2.*(resul.Dot(vref));
355
356   pt2d = gp_Pnt2d(X(1),X(2));
357   return Standard_True;
358 }
359
360
361 //=======================================================================
362 //function : Values
363 //purpose  : 
364 //=======================================================================
365
366 Standard_Boolean BlendFunc_CSCircular::Values(const math_Vector& X, math_Vector& F, math_Matrix& D)
367 {
368   gp_Vec d1u1,d1v1,d1c;
369   gp_Vec d2u1,d2v1,d2uv1;
370   gp_Vec ns,ncrossns,resul,temp,vref,nsov;
371
372   Standard_Real norm,ndotns,grosterme;
373
374   surf->D2(X(1),X(2),pts,d1u1,d1v1,d2u1,d2v1,d2uv1);
375   curv->D1(prmc,ptc,d1c);
376
377   ns = d1u1.Crossed(d1v1);
378   ncrossns = nplan.Crossed(ns);
379   norm = ncrossns.Magnitude();
380   if (norm < Eps) {
381     norm = 1.;
382 //#if DEB
383 //    cout << "CSCircular : Surface singuliere !" << endl;
384 //#endif
385   }
386
387   ndotns = nplan.Dot(ns);
388   nsov.SetLinearForm(ndotns/norm,nplan,-1./norm,ns);
389   vref.SetLinearForm(ray,nsov,gp_Vec(ptc,pts));
390
391   F(1) = nplan.XYZ().Dot(pts.XYZ()-ptc.XYZ());
392   F(2) = vref.SquareMagnitude() - ray*ray;
393
394   D(1,1) = nplan.Dot(d1u1);
395   D(1,2) = nplan.Dot(d1v1);
396
397    // Derivee par rapport a u
398   temp = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
399   grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
400   resul.SetLinearForm(-ray/norm*(grosterme*ndotns-nplan.Dot(temp)),nplan,
401                       ray*grosterme/norm,ns,
402                       -ray/norm,temp,
403                       d1u1);
404
405   D(2,1) = 2.*(resul.Dot(vref));
406
407   // Derivee par rapport a v
408   temp = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
409   grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
410   resul.SetLinearForm(-ray/norm*(grosterme*ndotns-nplan.Dot(temp)),nplan,
411                       ray*grosterme/norm,ns,
412                       -ray/norm,temp,
413                       d1v1);
414
415   D(2,2) = 2.*(resul.Dot(vref));
416
417   pt2d = gp_Pnt2d(X(1),X(2));
418   return Standard_True;
419 }
420
421
422 //=======================================================================
423 //function : PointOnS
424 //purpose  : 
425 //=======================================================================
426
427 const gp_Pnt& BlendFunc_CSCircular::PointOnS () const
428 {
429   return pts;
430 }
431
432 //=======================================================================
433 //function : PointOnC
434 //purpose  : 
435 //=======================================================================
436
437 const gp_Pnt& BlendFunc_CSCircular::PointOnC () const
438 {
439   return ptc;
440 }
441
442 //=======================================================================
443 //function : Pnt2d
444 //purpose  : 
445 //=======================================================================
446
447 const gp_Pnt2d& BlendFunc_CSCircular::Pnt2d () const
448 {
449   return pt2d;
450 }
451
452 //=======================================================================
453 //function : ParameterOnC
454 //purpose  : 
455 //=======================================================================
456
457 Standard_Real BlendFunc_CSCircular::ParameterOnC () const
458 {
459   return prmc;
460 }
461
462
463 //=======================================================================
464 //function : IsTangencyPoint
465 //purpose  : 
466 //=======================================================================
467
468 Standard_Boolean BlendFunc_CSCircular::IsTangencyPoint () const
469 {
470   return istangent;
471 }
472
473 //=======================================================================
474 //function : TangentOnS
475 //purpose  : 
476 //=======================================================================
477
478 const gp_Vec& BlendFunc_CSCircular::TangentOnS () const
479 {
480   if (istangent)
481     Standard_DomainError::Raise("BlendFunc_CSCircular::TangentOnS");
482   return tgs;
483 }
484
485 //=======================================================================
486 //function : TangentOnC
487 //purpose  : 
488 //=======================================================================
489
490 const gp_Vec& BlendFunc_CSCircular::TangentOnC () const
491 {
492   if (istangent)
493     Standard_DomainError::Raise("BlendFunc_CSCircular::TangentOnC");
494   return tgc;
495 }
496
497 //=======================================================================
498 //function : Tangent2d
499 //purpose  : 
500 //=======================================================================
501
502 const gp_Vec2d& BlendFunc_CSCircular::Tangent2d () const
503 {
504   if (istangent)
505     Standard_DomainError::Raise("BlendFunc_CSCircular::Tangent2d");
506   return tg2d;
507 }
508
509
510 //=======================================================================
511 //function : Tangent
512 //purpose  : 
513 //=======================================================================
514
515 void BlendFunc_CSCircular::Tangent(const Standard_Real U, const Standard_Real V,
516                                    gp_Vec& TgS, gp_Vec& NmS) const
517 {
518   gp_Pnt bid;
519   gp_Vec d1u,d1v,ns;
520   surf->D1(U,V,bid,d1u,d1v);
521   NmS = ns = d1u.Crossed(d1v);
522   
523   const Standard_Real norm = nplan.Crossed(ns).Magnitude();
524   ns.SetLinearForm(nplan.Dot(ns)/norm,nplan, -1./norm,ns);
525   if(ray > 0.) ns.Reverse();
526   TgS = nplan.Crossed(ns);
527   if (choix%2 == 1)
528     TgS.Reverse();
529 }
530
531
532 //=======================================================================
533 //function : Section
534 //purpose  : 
535 //=======================================================================
536
537 void BlendFunc_CSCircular::Section(const Standard_Real Param,
538                                    const Standard_Real U,
539                                    const Standard_Real V,
540                                    const Standard_Real W,
541                                    Standard_Real& Pdeb,
542                                    Standard_Real& Pfin,
543                                    gp_Circ& C)
544 {
545   gp_Vec d1u1,d1v1;
546   gp_Vec ns;//,temp;
547   Standard_Real norm;
548   gp_Pnt Center;
549   gp_Pnt ptgui;
550
551   guide->D1(Param,ptgui,d1gui);
552   nplan  = d1gui.Normalized();
553
554   surf->D1(U,V,pts,d1u1,d1v1);
555   ptc = curv->Value(W);
556
557   ns = d1u1.Crossed(d1v1);
558   norm = nplan.Crossed(ns).Magnitude();
559   ns.SetLinearForm(nplan.Dot(ns)/norm,nplan, -1./norm,ns);
560
561   Center.SetXYZ(pts.XYZ()+ray*ns.XYZ());
562   C.SetRadius(Abs(ray));
563
564   if(ray > 0.) ns.Reverse();
565
566   if (choix%2 == 0) {
567     C.SetPosition(gp_Ax2(Center,nplan,ns));
568   }
569   else {
570     C.SetPosition(gp_Ax2(Center,nplan.Reversed(),ns));
571   }
572   Pdeb = 0.;
573   Pfin = ElCLib::Parameter(C,ptc);
574 }
575
576
577 Standard_Boolean BlendFunc_CSCircular::Section(const Blend_Point& P, 
578                                            TColgp_Array1OfPnt& Poles, 
579                                            TColgp_Array1OfVec& DPoles, 
580                                            TColgp_Array1OfVec& D2Poles, 
581                                            TColgp_Array1OfPnt2d& Poles2d, 
582                                            TColgp_Array1OfVec2d& DPoles2d, 
583                                            TColgp_Array1OfVec2d& D2Poles2d, 
584                                            TColStd_Array1OfReal& Weigths, 
585                                            TColStd_Array1OfReal& DWeigths, 
586                                            TColStd_Array1OfReal& D2Weigths)
587 {
588   return Blend_CSFunction::Section(P,Poles,DPoles,D2Poles,Poles2d,DPoles2d,D2Poles2d,Weigths,DWeigths,D2Weigths);
589 }
590
591 //=======================================================================
592 //function : GetSection
593 //purpose  : 
594 //=======================================================================
595
596 Standard_Boolean BlendFunc_CSCircular::GetSection(const Standard_Real Param,
597                                                   const Standard_Real U,
598                                                   const Standard_Real V,
599                                                   const Standard_Real /*W*/,
600                                                   TColgp_Array1OfPnt& tabP,
601                                                   TColgp_Array1OfVec& tabV)
602 {
603   Standard_Integer NbPoint=tabP.Length();
604   if (NbPoint != tabV.Length() || NbPoint < 2) {Standard_RangeError::Raise();}
605
606   Standard_Integer i, lowp = tabP.Lower(), lowv = tabV.Lower();
607
608   gp_Vec d1u1,d1v1,d2u1,d2v1,d2uv1,d1c,d2c; //,d1u2,d1v2;
609   gp_Vec ns,dnplan,dnw,dn2w,ncrn,dncrn,ns2;
610   gp_Vec ncrossns,resul;
611   gp_Vec resulu,resulv,temp;
612
613   Standard_Real norm,ndotns,grosterme;
614   Standard_Real lambda,Cosa,Sina;
615   Standard_Real Angle = 0.,Dangle = 0.;
616   math_Vector sol(1,2),valsol(1,2),secmember(1,2);
617   math_Matrix gradsol(1,2,1,2);
618
619   Set(Param);
620   dnplan.SetLinearForm(1./normtg,d2gui,
621                        -1./normtg*(nplan.Dot(d2gui)),nplan);
622
623   curv->D2(prmc,ptc,d1c,d2c);
624   surf->D2(U,V,pts,d1u1,d1v1,d2u1,d2v1,d2uv1);
625
626   ns = d1u1.Crossed(d1v1);
627   ncrossns = nplan.Crossed(ns);
628   ndotns = nplan.Dot(ns);
629   norm = ncrossns.Magnitude();
630
631   temp.SetXYZ(pts.XYZ() - ptc.XYZ());
632   secmember(1) = dprmc*(nplan.Dot(d1c)) - dnplan.Dot(temp);    
633
634   ns2.SetLinearForm(ndotns/norm,nplan, -1./norm,ns);
635
636   // Derivee de n1 par rapport a w (param sur ligne guide)
637
638   grosterme = ncrossns.Dot(dnplan.Crossed(ns))/norm/norm;
639   dnw.SetLinearForm((dnplan.Dot(ns)-grosterme*ndotns)/norm,nplan,
640                      ndotns/norm,dnplan,
641                      grosterme/norm,ns);
642
643   temp.SetLinearForm(ray,dnw,-dprmc,d1c);
644   resul.SetLinearForm(ray,ns2,gp_Vec(ptc,pts));
645
646   secmember(2) = -2.*(resul.Dot(temp));
647
648   sol(1) = U; sol(2) = V;
649   Values(sol,valsol,gradsol);
650   math_Gauss Resol(gradsol);
651
652   if (Resol.IsDone()) {
653     
654     Resol.Solve(secmember);
655     
656     tgs.SetLinearForm(secmember(1),d1u1,secmember(2),d1v1);
657     tgc = dprmc*d1c;
658
659     // Derivee de n1 par rapport a u1
660     temp = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
661     grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
662     resulu.SetLinearForm(-(grosterme*ndotns-nplan.Dot(temp))/norm,nplan,
663                          grosterme/norm,ns,
664                          -1./norm,temp);
665
666     // Derivee de n1 par rapport a v1
667     temp = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
668     grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
669     resulv.SetLinearForm(-(grosterme*ndotns-nplan.Dot(temp))/norm,nplan,
670                          grosterme/norm,ns,
671                          -1./norm,temp);
672
673
674     dnw.SetLinearForm(secmember(1),resulu,secmember(2),resulv,dnw);
675     ns = ns2;
676     dn2w.SetLinearForm(ray, dnw, -1., tgc, tgs);
677     norm = resul.Magnitude();
678     dn2w.Divide(norm);
679     ns2 = -resul.Normalized();
680     dn2w.SetLinearForm(ns2.Dot(dn2w),ns2,-1.,dn2w);
681
682     if(ray > 0.) {
683       ns.Reverse();
684       dnw.Reverse();
685     }
686
687     if (choix%2 != 0) {
688       nplan.Reverse();
689       dnplan.Reverse();
690     }
691
692
693     tabP(lowp) = pts;
694     tabP(lowp+NbPoint-1)  = ptc;
695
696     tabV(lowv) = tgs;
697     tabV(lowv+NbPoint-1)  = tgc;
698
699     if (NbPoint >2) {
700
701       Cosa = ns.Dot(ns2);
702       Sina = nplan.Dot(ns.Crossed(ns2));
703       Angle = ACos(Cosa);
704       if (Sina <0.) {
705         Angle = 2.*M_PI - Angle;
706       }
707       Dangle = -(dnw.Dot(ns2) + ns.Dot(dn2w))/Sina;
708       ncrn = nplan.Crossed(ns);
709       dncrn = dnplan.Crossed(ns).Added(nplan.Crossed(dnw));
710     }
711
712     for (i=2; i <= NbPoint-1; i++) {
713       lambda = (Standard_Real)(i-1)/(Standard_Real)(NbPoint-1);
714       Cosa = Cos(lambda*Angle);
715       Sina = Sin(lambda*Angle);
716       tabP(lowp+i-1).SetXYZ(pts.XYZ() 
717                      +Abs(ray)*((Cosa-1)*ns.XYZ() + Sina*ncrn.XYZ()));
718
719       temp.SetLinearForm(-Sina,ns,Cosa,ncrn);
720       temp.Multiply(lambda*Dangle);
721       temp.Add(((Cosa-1)*dnw).Added(Sina*dncrn));
722       temp.Multiply(Abs(ray));
723       temp.Add(tgs);
724       tabV(lowv+i-1)= temp;
725     } 
726     return Standard_True;
727   }
728   return Standard_False;
729 }
730
731 //=======================================================================
732 //function : IsRational
733 //purpose  : 
734 //=======================================================================
735
736 Standard_Boolean BlendFunc_CSCircular::IsRational() const
737 {
738   return (mySShape==BlendFunc_Rational || mySShape==BlendFunc_QuasiAngular);
739 }
740
741 //=======================================================================
742 //function : GetSectionSize
743 //purpose  : 
744 //=======================================================================
745
746 Standard_Real BlendFunc_CSCircular::GetSectionSize() const 
747 {
748   return maxang*Abs(ray);
749 }
750
751 //=======================================================================
752 //function : GetMinimalWeight
753 //purpose  : 
754 //=======================================================================
755 void BlendFunc_CSCircular::GetMinimalWeight(TColStd_Array1OfReal& Weigths) const 
756 {
757   BlendFunc::GetMinimalWeights(mySShape, myTConv,minang,maxang,Weigths);
758   // On suppose que cela ne depend pas du Rayon! 
759 }
760
761 //=======================================================================
762 //function : NbIntervals
763 //purpose  : 
764 //=======================================================================
765
766 Standard_Integer BlendFunc_CSCircular::NbIntervals (const GeomAbs_Shape S) const
767 {
768   return curv->NbIntervals(BlendFunc::NextShape(S));
769 }
770
771
772 //=======================================================================
773 //function : Intervals
774 //purpose  : 
775 //=======================================================================
776
777 void BlendFunc_CSCircular::Intervals (TColStd_Array1OfReal& T, const GeomAbs_Shape S) const
778 {
779   curv->Intervals(T, BlendFunc::NextShape(S));
780 }
781
782 //=======================================================================
783 //function : GetShape
784 //purpose  : 
785 //=======================================================================
786
787 void BlendFunc_CSCircular::GetShape (Standard_Integer& NbPoles,
788                                      Standard_Integer& NbKnots,
789                                      Standard_Integer& Degree,
790                                      Standard_Integer& NbPoles2d)
791 {
792   NbPoles2d = 1;
793   BlendFunc::GetShape(mySShape,maxang,NbPoles,NbKnots,Degree,myTConv);
794 }
795
796
797 //=======================================================================
798 //function : GetTolerance
799 //purpose  : Determine les Tolerances a utiliser dans les approximations.
800 //=======================================================================
801 void BlendFunc_CSCircular::GetTolerance(const Standard_Real BoundTol, 
802                                       const Standard_Real SurfTol, 
803                                       const Standard_Real AngleTol, 
804                                       math_Vector& Tol3d, 
805                                       math_Vector& Tol1d) const
806 {
807   const Standard_Integer low = Tol3d.Lower();
808   const Standard_Integer up = Tol3d.Upper();
809   const Standard_Real Tol = GeomFill::GetTolerance(myTConv, minang, ray, AngleTol, SurfTol);
810   Tol1d.Init(SurfTol);
811   Tol3d.Init(SurfTol);
812   Tol3d(low+1) = Tol3d(up-1) = Min( Tol, SurfTol);
813   Tol3d(low) = Tol3d(up) = Min( Tol, BoundTol);
814 }
815
816
817 //=======================================================================
818 //function : Knots
819 //purpose  : 
820 //=======================================================================
821
822 void BlendFunc_CSCircular::Knots(TColStd_Array1OfReal& TKnots)
823 {
824   GeomFill::Knots(myTConv,TKnots);
825 }
826
827 //=======================================================================
828 //function : Mults
829 //purpose  : 
830 //=======================================================================
831
832 void BlendFunc_CSCircular::Mults(TColStd_Array1OfInteger& TMults)
833 {
834   GeomFill::Mults(myTConv,TMults);
835 }
836
837
838 //=======================================================================
839 //function : Section
840 //purpose  : 
841 //=======================================================================
842
843 void BlendFunc_CSCircular::Section(const Blend_Point& P,
844                                    TColgp_Array1OfPnt& Poles,
845                                    TColgp_Array1OfPnt2d& Poles2d,
846                                    TColStd_Array1OfReal& Weights)
847 {
848   gp_Vec d1u1,d1v1;//,d1;
849   gp_Vec ns,ns2;//,temp,np2;
850   gp_Pnt Center;
851
852   Standard_Real norm,u1,v1;
853
854   Standard_Real prm = P.Parameter();
855   Standard_Integer low = Poles.Lower();
856   Standard_Integer upp = Poles.Upper();
857
858   Set(prm);
859
860   P.ParametersOnS(u1,v1);
861   surf->D1(u1,v1,pts,d1u1,d1v1);
862   ptc = curv->Value(prmc);
863
864   Poles2d(Poles2d.Lower()).SetCoord(u1,v1);
865
866   // Cas Linear
867   if (mySShape == BlendFunc_Linear) {
868     Poles(low) = pts;
869     Poles(upp) = ptc;
870     Weights(low) = 1.0;
871     Weights(upp) = 1.0;
872     return;
873   }
874
875   ns = d1u1.Crossed(d1v1);
876   norm = nplan.Crossed(ns).Magnitude();
877
878   ns.SetLinearForm(nplan.Dot(ns)/norm,nplan, -1./norm,ns);
879
880   Center.SetXYZ(pts.XYZ()+ray*ns.XYZ());
881
882   ns2 = gp_Vec(Center,ptc).Normalized();
883   if(ray > 0.) ns.Reverse();
884
885   if (choix%2 != 0) {
886     nplan.Reverse();
887   }
888
889   GeomFill::GetCircle(myTConv,
890                       ns, ns2, 
891                       nplan, pts, ptc,
892                       Abs(ray), Center, 
893                       Poles, Weights);
894 }
895  
896
897 //=======================================================================
898 //function : Section
899 //purpose  : 
900 //=======================================================================
901
902 Standard_Boolean BlendFunc_CSCircular::Section
903   (const Blend_Point& P,
904    TColgp_Array1OfPnt& Poles,
905    TColgp_Array1OfVec& DPoles,
906    TColgp_Array1OfPnt2d& Poles2d,
907    TColgp_Array1OfVec2d& DPoles2d,
908    TColStd_Array1OfReal& Weights,
909    TColStd_Array1OfReal& DWeights)
910 {
911   gp_Vec d1u1,d1v1,d2u1,d2v1,d2uv1,d1,d2;
912   gp_Vec ns,ns2,dnplan,dnw,dn2w; //,np2,dnp2;
913   gp_Vec ncrossns;;
914   gp_Vec resulu,resulv,temp,tgct,resul;
915
916   gp_Pnt Center;
917
918   Standard_Real norm,ndotns,grosterme;;
919
920   math_Vector sol(1,2),valsol(1,2),secmember(1,2);
921   math_Matrix gradsol(1,2,1,2);
922
923   Standard_Real prm = P.Parameter();
924   Standard_Integer low = Poles.Lower();
925   Standard_Integer upp = Poles.Upper();
926   Standard_Boolean istgt;
927
928   Set(prm);
929   dnplan.SetLinearForm(1./normtg,d2gui,
930                        -1./normtg*(nplan.Dot(d2gui)),nplan);
931
932   curv->D2(prmc,ptc,d1,d2);
933   P.ParametersOnS(sol(1),sol(2));
934   surf->D2(sol(1),sol(2),pts,d1u1,d1v1,d2u1,d2v1,d2uv1);
935
936   ns = d1u1.Crossed(d1v1);
937   ncrossns = nplan.Crossed(ns);
938   ndotns = nplan.Dot(ns);
939   norm = ncrossns.Magnitude();
940
941   ns2.SetLinearForm(ndotns/norm,nplan, -1./norm,ns);
942   temp.SetXYZ(pts.XYZ() - ptc.XYZ());
943
944   secmember(1) = dprmc*(nplan.Dot(d1)) - dnplan.Dot(temp);    
945   
946   // Derivee de n1 par rapport a w
947
948   grosterme = ncrossns.Dot(dnplan.Crossed(ns))/norm/norm;
949   dnw.SetLinearForm((dnplan.Dot(ns)-grosterme*ndotns)/norm,nplan,
950                      ndotns/norm,dnplan,
951                      grosterme/norm,ns);
952
953   temp.SetLinearForm(ray,dnw,-dprmc,d1);
954   resul.SetLinearForm(ray,ns2,gp_Vec(ptc,pts));
955
956   secmember(2) = -2.*(resul.Dot(temp));
957
958   Values(sol,valsol,gradsol);
959   math_Gauss Resol(gradsol);
960
961   if (Resol.IsDone()) {
962     
963     Resol.Solve(secmember);
964     
965     tgs.SetLinearForm(secmember(1),d1u1,secmember(2),d1v1);
966     tgc = dprmc*d1;
967
968     // Derivee de n1 par rapport a u1
969     temp = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
970     grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
971     resulu.SetLinearForm(-(grosterme*ndotns-nplan.Dot(temp))/norm,nplan,
972                          grosterme/norm,ns,
973                          -1./norm,temp);
974
975     // Derivee de n1 par rapport a v1
976     temp = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
977     grosterme = ncrossns.Dot(nplan.Crossed(temp))/norm/norm;
978     resulv.SetLinearForm(-(grosterme*ndotns-nplan.Dot(temp))/norm,nplan,
979                          grosterme/norm,ns,
980                          -1./norm,temp);
981
982
983     dnw.SetLinearForm(secmember(1),resulu,secmember(2),resulv,dnw);
984     ns = ns2;
985     dn2w.SetLinearForm(ray, dnw,-1., tgc, tgs);
986     norm = resul.Magnitude();
987     dn2w.Divide(norm);
988     ns2 = -resul.Normalized();
989     dn2w.SetLinearForm(ns2.Dot(dn2w),ns2,-1.,dn2w);
990
991     istgt = Standard_False;
992   }
993   else {
994     ns.SetLinearForm(ndotns/norm,nplan, -1./norm,ns);
995     ns2 = -resul.Normalized();
996     istgt = Standard_True;
997   }
998
999   // Les poles 2d
1000   
1001   Poles2d(Poles2d.Lower()).SetCoord(sol(1),sol(2));
1002   if (!istgt) {
1003     DPoles2d(Poles2d.Lower()).SetCoord(secmember(1),secmember(2));
1004   }
1005   
1006   if (mySShape == BlendFunc_Linear) {
1007     Poles(low) = pts;
1008     Poles(upp) = ptc;
1009     Weights(low) = 1.0;
1010     Weights(upp) = 1.0;
1011     if (!istgt) {
1012       DPoles(low) = tgs;
1013       DPoles(upp) = tgc;
1014       DWeights(low) = 0.0;
1015       DWeights(upp) = 0.0;
1016     }
1017     return (!istgt);
1018   }
1019
1020
1021   // Cas du cercle
1022   Center.SetXYZ(pts.XYZ()+ray*ns.XYZ());
1023   if (!istgt) {
1024     tgct = tgs.Added(ray*dnw);
1025   }
1026
1027   if(ray > 0.) {
1028     ns.Reverse();
1029     if(!istgt) { dnw.Reverse(); }
1030   }
1031
1032   if (choix%2 != 0) {
1033     nplan.Reverse();
1034     dnplan.Reverse();
1035   }
1036
1037   if (!istgt) {
1038     return GeomFill::GetCircle(myTConv, 
1039                                ns, ns2, 
1040                                dnw, dn2w, 
1041                                nplan, dnplan, 
1042                                pts, ptc, 
1043                                tgs, tgc, 
1044                                Abs(ray), 0, 
1045                                Center, tgct, 
1046                                Poles, 
1047                                DPoles,
1048                                Weights, 
1049                                DWeights); 
1050   }
1051   else {
1052     GeomFill::GetCircle(myTConv,
1053                         ns, ns2, 
1054                         nplan, pts, ptc,
1055                         Abs(ray), Center, 
1056                         Poles, Weights);
1057     return Standard_False;
1058   }
1059 }
1060
1061 void BlendFunc_CSCircular::Resolution(const Standard_Integer , const Standard_Real Tol,
1062                                       Standard_Real& TolU, Standard_Real& TolV) const
1063 {
1064   TolU = surf->UResolution(Tol);
1065   TolV = surf->VResolution(Tol);
1066 }