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