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