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