0031035: Coding - uninitialized class fields reported by Visual Studio Code Analysis
[occt.git] / src / BlendFunc / BlendFunc_ConstRad.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 09/09/1996 PMN Adde Nb(Intervalls), IsRationnal
18 //                         Optimisation, use of GetCircle
19 // Modified 20/02/1998 PMN Singular surfaces management
20
21 #include <Adaptor3d_HCurve.hxx>
22 #include <Adaptor3d_HSurface.hxx>
23 #include <Blend_Point.hxx>
24 #include <BlendFunc.hxx>
25 #include <BlendFunc_ConstRad.hxx>
26 #include <ElCLib.hxx>
27 #include <GeomFill.hxx>
28 #include <gp.hxx>
29 #include <gp_Ax1.hxx>
30 #include <gp_Circ.hxx>
31 #include <gp_Pnt.hxx>
32 #include <gp_Vec.hxx>
33 #include <gp_Vec2d.hxx>
34 #include <math_Gauss.hxx>
35 #include <math_Matrix.hxx>
36 #include <math_SVD.hxx>
37 #include <Precision.hxx>
38 #include <Standard_DomainError.hxx>
39 #include <Standard_NotImplemented.hxx>
40 #include <Standard_TypeDef.hxx>
41
42 #define Eps 1.e-15
43
44
45 //=======================================================================
46 //function : BlendFunc_ConstRad
47 //purpose  : 
48 //=======================================================================
49
50 BlendFunc_ConstRad::BlendFunc_ConstRad(const Handle(Adaptor3d_HSurface)& S1,
51                                        const Handle(Adaptor3d_HSurface)& S2,
52                                        const Handle(Adaptor3d_HCurve)& C)
53                                  :
54                   surf1(S1),surf2(S2),
55                   curv(C), tcurv(C),
56                   istangent(Standard_True), param(0.0),
57                   ray1(0.0), ray2(0.0),
58                   choix(0), xval(1, 4),
59                   E(1,4), DEDX(1,4,1,4),  DEDT(1,4), 
60                   D2EDX2(4,4,4),
61                   D2EDXDT(1,4,1,4), D2EDT2(1,4),
62                   maxang(RealFirst()), minang(RealLast()),
63                   distmin(RealLast()),
64                   mySShape(BlendFunc_Rational)
65
66 // Initialisaton of cash control variables.
67   tval = -9.876e100;
68   xval.Init(-9.876e100);
69   myXOrder = -1;
70   myTOrder = -1;
71 }
72
73 //=======================================================================
74 //function : NbEquations
75 //purpose  : 
76 //=======================================================================
77
78 Standard_Integer BlendFunc_ConstRad::NbEquations () const
79 {
80   return 4;
81 }
82
83
84 //=======================================================================
85 //function : Set
86 //purpose  : 
87 //=======================================================================
88
89 void BlendFunc_ConstRad::Set(const Standard_Real Radius, const Standard_Integer Choix)
90 {
91   choix = Choix;
92   switch (choix) {
93   case 1:
94   case 2:
95     {
96       ray1 = -Radius;
97       ray2 = -Radius;
98     }
99     break;
100   case 3:
101   case 4:
102     {
103       ray1 = Radius;
104       ray2 = -Radius;
105     }
106     break;
107   case 5:
108   case 6:
109     {
110       ray1 = Radius;
111       ray2 = Radius;
112     }
113     break;
114   case 7:
115   case 8:
116     {
117       ray1 = -Radius;
118       ray2 = Radius;
119     }
120     break;
121   default:
122     ray1 = ray2 = -Radius;
123   }
124 }
125
126 //=======================================================================
127 //function : Set
128 //purpose  : 
129 //=======================================================================
130
131 void BlendFunc_ConstRad::Set(const BlendFunc_SectionShape TypeSection)
132 {
133   mySShape = TypeSection;
134 }
135
136
137 //=======================================================================
138 //function : ComputeValues
139 //purpose  : OBLIGATORY passage for all calculations
140 //           This method manages positioning on Surfaces and Curves
141 //           Calculate the equations and their partial derivates
142 //           Stock certain intermediate results in fields to 
143 //           use in other methods.
144 //=======================================================================
145
146 Standard_Boolean BlendFunc_ConstRad::ComputeValues(const math_Vector& X,
147                                                    const Standard_Integer Order,
148                                                    const Standard_Boolean byParam,
149                                                    const Standard_Real Param)
150 {
151  // static declaration to avoid systematic reallocation
152  
153  static gp_Vec d3u1,d3v1,d3uuv1,d3uvv1,d3u2,d3v2,d3uuv2,d3uvv2; 
154  static gp_Vec d1gui, d2gui, d3gui;
155  static gp_Pnt ptgui;
156  static Standard_Real invnormtg, dinvnormtg;
157  Standard_Real T =  Param, aux;
158
159  // Case of implicite parameter
160  if ( !byParam) { T = param;}
161
162  // Is the work already done ?
163  Standard_Boolean myX_OK = (Order<=myXOrder) ;
164  for (Standard_Integer ii=1; ((ii<=X.Length()) && myX_OK); ii++) {
165    myX_OK = ( X(ii) == xval(ii) );
166  }
167
168  Standard_Boolean t_OK =( (T == tval) 
169                        && ((Order<=myTOrder)||(!byParam)) ); 
170
171  if (myX_OK && (t_OK) ) {
172    return Standard_True;
173  }
174
175  // Processing of t
176  if (!t_OK) {
177    tval = T;
178    if (byParam) { myTOrder = Order;}
179    else         { myTOrder = 0;}
180    //----- Positioning on the curve ----------------
181    switch (myTOrder) {
182    case 0 :
183      {
184        tcurv->D1(T, ptgui, d1gui);
185        nplan  = d1gui.Normalized();
186      }
187      break;
188
189    case 1 :
190      { 
191        tcurv->D2(T,ptgui,d1gui,d2gui);
192        nplan  = d1gui.Normalized();
193        invnormtg = ((Standard_Real) 1 ) / d1gui.Magnitude();
194        dnplan.SetLinearForm(invnormtg, d2gui,
195                             -invnormtg*(nplan.Dot(d2gui)), nplan);
196      break;
197      }
198    case 2 :
199      { 
200        tcurv->D3(T,ptgui,d1gui,d2gui,d3gui);
201        nplan  = d1gui.Normalized();
202        invnormtg = ((Standard_Real) 1 ) / d1gui.Magnitude();
203        dnplan.SetLinearForm(invnormtg, d2gui,
204                            -invnormtg*(nplan.Dot(d2gui)), nplan);
205        dinvnormtg = - nplan.Dot(d2gui)*invnormtg*invnormtg;
206        d2nplan.SetLinearForm(invnormtg,  d3gui,
207                              dinvnormtg, d2gui);
208        aux = dinvnormtg*(nplan.Dot(d2gui)) + invnormtg*( dnplan.Dot(d2gui) 
209                                                        + nplan.Dot(d3gui) ); 
210        d2nplan.SetLinearForm(-invnormtg*(nplan.Dot(d2gui)), dnplan,
211                              -aux, nplan, d2nplan );
212        break;
213    }
214    default:
215      return Standard_False;
216    }
217  }
218
219  // Processing of X
220  if (!myX_OK) {
221    xval = X;
222    myXOrder = Order;
223    //-------------- Positioning on surfaces -----------------
224    switch (myXOrder) {
225    case 0 :
226      {
227        surf1->D1(X(1),X(2),pts1,d1u1,d1v1);
228        nsurf1 = d1u1.Crossed(d1v1);
229        surf2->D1(X(3),X(4),pts2,d1u2,d1v2);
230        nsurf2 = d1u2.Crossed(d1v2);
231        break;
232      }
233    case 1 :
234      {
235        surf1->D2(X(1),X(2),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
236        nsurf1 = d1u1.Crossed(d1v1);
237        dns1u1 = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
238        dns1v1 = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
239        surf2->D2(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
240        nsurf2 = d1u2.Crossed(d1v2);
241        dns1u2 = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
242        dns1v2 = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
243        break;
244      }
245    case 2 :
246      {
247        surf1->D3(X(1),X(2),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1,d3u1,d3v1,d3uuv1,d3uvv1);
248        nsurf1 = d1u1.Crossed(d1v1);
249        dns1u1 = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
250        dns1v1 = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
251
252        surf2->D3(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2,d3u2,d3v2,d3uuv2,d3uvv2);
253        nsurf2 = d1u2.Crossed(d1v2);
254        dns1u2 = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
255        dns1v2 = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
256        break;
257      }
258    default:
259      return Standard_False;
260    }
261    // Case of degenerated surfaces
262    if (nsurf1.Magnitude() < Eps ) {
263      //gp_Vec normal;
264      gp_Pnt2d P(X(1), X(2)); 
265      if (Order == 0) BlendFunc::ComputeNormal(surf1, P, nsurf1);
266      else  BlendFunc::ComputeDNormal(surf1, P, nsurf1, dns1u1, dns1v1);
267    }
268   if (nsurf2.Magnitude() < Eps) {
269      //gp_Vec normal;
270      gp_Pnt2d P(X(3), X(4));
271      if (Order==0) BlendFunc::ComputeNormal(surf2, P, nsurf2);
272      else BlendFunc::ComputeDNormal(surf2, P, nsurf2, dns1u2, dns1v2);
273    }
274  }
275
276  // -------------------- Positioning of order 0 ---------------------
277  Standard_Real invnorm1, invnorm2, ndotns1, ndotns2, theD;
278  gp_Vec ncrossns1,ncrossns2,resul,temp;
279
280  theD =  - (nplan.XYZ().Dot(ptgui.XYZ()));
281
282  E(1) = (nplan.X() * (pts1.X() + pts2.X()) + 
283          nplan.Y() * (pts1.Y() + pts2.Y()) + 
284          nplan.Z() * (pts1.Z() + pts2.Z())) /2 + theD;
285
286  ncrossns1 = nplan.Crossed(nsurf1);
287  ncrossns2 = nplan.Crossed(nsurf2);
288
289  invnorm1 = ncrossns1.Magnitude();
290  invnorm2 = ncrossns2.Magnitude();
291
292  if (invnorm1 > Eps) invnorm1 = ((Standard_Real) 1) /invnorm1;
293  else {
294    invnorm1 = 1; // Unsatisfactory, but it is not necessary to crash
295 #ifdef OCCT_DEBUG
296    std::cout << " ConstRad : Surface singuliere " << std::endl;
297 #endif
298  }
299  if (invnorm2 > Eps) invnorm2 = ((Standard_Real) 1) /invnorm2;
300  else {
301    invnorm2 = 1; //  Unsatisfactory, but it is not necessary to crash
302 #ifdef OCCT_DEBUG
303    std::cout << " ConstRad : Surface singuliere " << std::endl;
304 #endif
305  }
306
307  ndotns1 = nplan.Dot(nsurf1);
308  ndotns2 = nplan.Dot(nsurf2);
309
310  temp.SetLinearForm(ndotns1,nplan,-1.,nsurf1);
311  temp.Multiply(invnorm1);
312
313  resul.SetLinearForm(ray1,temp,gp_Vec(pts2,pts1));
314  temp.SetLinearForm(ndotns2,nplan,-1.,nsurf2);
315  temp.Multiply(invnorm2);
316  resul.Subtract(ray2*temp);
317
318  E(2) = resul.X();
319  E(3) = resul.Y();
320  E(4) = resul.Z();
321
322  // -------------------- Positioning of order 1 ---------------------
323  if (Order >= 1) {
324    Standard_Real  grosterme, cube, carre;
325   
326
327    DEDX(1,1) = nplan.Dot(d1u1)/2;
328    DEDX(1,2) = nplan.Dot(d1v1)/2;
329    DEDX(1,3) = nplan.Dot(d1u2)/2;
330    DEDX(1,4) = nplan.Dot(d1v2)/2;
331
332    cube =invnorm1*invnorm1*invnorm1;
333    // Derived in relation to u1
334    grosterme = - ncrossns1.Dot(nplan.Crossed(dns1u1))*cube;
335    dndu1.SetLinearForm( grosterme*ndotns1
336                       + invnorm1*nplan.Dot(dns1u1), nplan,
337                       - grosterme, nsurf1,
338                       - invnorm1,  dns1u1 );
339
340    resul.SetLinearForm(ray1, dndu1, d1u1);
341    DEDX(2,1) = resul.X();
342    DEDX(3,1) = resul.Y();
343    DEDX(4,1) = resul.Z();
344
345    // Derived in relation to v1
346
347    grosterme = - ncrossns1.Dot(nplan.Crossed(dns1v1))*cube;
348    dndv1.SetLinearForm( grosterme*ndotns1
349                        +invnorm1*nplan.Dot(dns1v1), nplan,
350                        -grosterme, nsurf1,
351                        -invnorm1, dns1v1);
352    
353    resul.SetLinearForm(ray1, dndv1, d1v1);
354    DEDX(2,2) = resul.X();
355    DEDX(3,2) = resul.Y();
356    DEDX(4,2) = resul.Z();
357
358    cube = invnorm2*invnorm2*invnorm2;
359    // Derived in relation to u2
360    grosterme = - ncrossns2.Dot(nplan.Crossed(dns1u2))*cube;
361    dndu2.SetLinearForm( grosterme*ndotns2
362                        +invnorm2*nplan.Dot(dns1u2), nplan,
363                        -grosterme, nsurf2,
364                        -invnorm2, dns1u2);
365
366    resul.SetLinearForm(-ray2, dndu2, -1, d1u2);
367    DEDX(2,3) = resul.X();
368    DEDX(3,3) = resul.Y();
369    DEDX(4,3) = resul.Z();
370
371    // Derived in relation to v2
372    grosterme = -ncrossns2.Dot(nplan.Crossed(dns1v2))*cube;
373    dndv2.SetLinearForm( grosterme*ndotns2
374                        +invnorm2*nplan.Dot(dns1v2), nplan,
375                        -grosterme, nsurf2,
376                        -invnorm2 , dns1v2);
377
378    resul.SetLinearForm(-ray2,dndv2, -1, d1v2);   
379    DEDX(2,4) = resul.X();
380    DEDX(3,4) = resul.Y();
381    DEDX(4,4) = resul.Z();
382
383    if (byParam) {
384      temp.SetXYZ( (pts1.XYZ()+pts2.XYZ())/2 - ptgui.XYZ());
385      // Derived from n1 in relation to w     
386      grosterme = ncrossns1.Dot(dnplan.Crossed(nsurf1))*invnorm1*invnorm1;
387      dn1w.SetLinearForm((dnplan.Dot(nsurf1)-grosterme*ndotns1)*invnorm1, nplan,
388                          ndotns1*invnorm1,dnplan,
389                          grosterme*invnorm1,nsurf1);
390   
391      // Derivee from n2 in relation to w
392      grosterme = ncrossns2.Dot(dnplan.Crossed(nsurf2))*invnorm2*invnorm2;
393      dn2w.SetLinearForm((dnplan.Dot(nsurf2)-grosterme*ndotns2)*invnorm2,nplan,
394                         ndotns2*invnorm2,dnplan,
395                         grosterme*invnorm2,nsurf2);
396
397
398      DEDT(1) =  dnplan.Dot(temp) - 1./invnormtg ;
399      DEDT(2) =  ray1*dn1w.X() - ray2*dn2w.X();
400      DEDT(3) =  ray1*dn1w.Y() - ray2*dn2w.Y();
401      DEDT(4) =  ray1*dn1w.Z() - ray2*dn2w.Z();
402    }
403    // ------   Positioning of order 2  -----------------------------
404    if (Order == 2) {
405 //     gp_Vec d2ndu1,  d2ndu2, d2ndv1, d2ndv2, d2nduv1, d2nduv2;
406      gp_Vec d2ns1u1,  d2ns1u2, d2ns1v1, d2ns1v2, d2ns1uv1, d2ns1uv2;
407      Standard_Real uterm, vterm, smallterm, p1, p2, p12;
408      Standard_Real DPrim, DSecn;
409      D2EDX2.Init(0);
410
411      D2EDX2(1, 1, 1) = nplan.Dot(d2u1)/2;
412      D2EDX2(1, 2, 1) = D2EDX2(1, 1, 2) = nplan.Dot(d2uv1)/2;
413      D2EDX2(1, 2, 2) = nplan.Dot(d2v1)/2;
414
415      D2EDX2(1, 3, 3) = nplan.Dot(d2u2)/2;
416      D2EDX2(1, 4, 3) = D2EDX2(1, 3, 4) = nplan.Dot(d2uv2)/2;
417      D2EDX2(1, 4, 4) = nplan.Dot(d2v2)/2;
418                                // ================
419                                // ==  Surface 1 ==
420                                // ================
421      carre = invnorm1*invnorm1;
422      cube  = carre*invnorm1;
423      // Derived double compared to u1       
424        // Derived from the norm
425      d2ns1u1.SetLinearForm(1, d3u1.Crossed(d1v1),
426                            2, d2u1.Crossed(d2uv1),
427                            1, d1u1.Crossed(d3uuv1));
428      DPrim = ncrossns1.Dot(nplan.Crossed(dns1u1));
429      smallterm =  - 2*DPrim*cube;
430      DSecn = ncrossns1.Dot(nplan.Crossed(d2ns1u1))
431           + (nplan.Crossed(dns1u1)).SquareMagnitude();
432      grosterme  = (3*DPrim*DPrim*carre - DSecn) * cube;     
433
434      temp.SetLinearForm( grosterme*ndotns1, nplan,
435                         -grosterme        , nsurf1);     
436      p1 =  nplan.Dot(dns1u1);
437      p2 =  nplan.Dot(d2ns1u1);    
438      d2ndu1.SetLinearForm( invnorm1*p2 
439                           +smallterm*p1, nplan,
440                           -smallterm,    dns1u1,
441                           -invnorm1,     d2ns1u1);
442      d2ndu1 += temp;
443      resul.SetLinearForm(ray1, d2ndu1, d2u1);
444      D2EDX2(2,1,1) = resul.X();
445      D2EDX2(3,1,1) = resul.Y();
446      D2EDX2(4,1,1) = resul.Z();
447
448      // Derived double compared to u1, v1       
449        // Derived from the norm
450      d2ns1uv1 =  (d3uuv1.Crossed(d1v1))
451               +  (d2u1  .Crossed(d2v1))
452               +  (d1u1  .Crossed(d3uvv1));
453      uterm =  ncrossns1.Dot(nplan.Crossed(dns1u1));
454      vterm =  ncrossns1.Dot(nplan.Crossed(dns1v1));
455      DSecn = (nplan.Crossed(dns1v1)).Dot(nplan.Crossed(dns1u1))
456            +  ncrossns1.Dot(nplan.Crossed(d2ns1uv1));
457      grosterme  = (3*uterm*vterm*carre-DSecn)*cube;
458      uterm *= -cube; //and only now
459      vterm *= -cube;
460        
461      p1 = nplan.Dot(dns1u1);
462      p2 = nplan.Dot(dns1v1);
463      temp.SetLinearForm(  grosterme*ndotns1, nplan,
464                         - grosterme, nsurf1,
465                         - invnorm1, d2ns1uv1);
466      d2nduv1.SetLinearForm( invnorm1*nplan.Dot(d2ns1uv1)
467                           + uterm*p2 
468                           + vterm*p1, nplan,
469                           - uterm,    dns1v1,
470                           - vterm,    dns1u1);
471                         
472      d2nduv1 += temp;
473      resul.SetLinearForm(ray1, d2nduv1, d2uv1);
474
475      D2EDX2(2,2,1) = D2EDX2(2,1,2) = resul.X();
476      D2EDX2(3,2,1) = D2EDX2(3,1,2) = resul.Y();
477      D2EDX2(4,2,1) = D2EDX2(4,1,2) = resul.Z();    
478
479      // Derived double compared to v1       
480        // Derived from the norm
481      d2ns1v1.SetLinearForm(1, d1u1.Crossed(d3v1),
482                            2, d2uv1.Crossed(d2v1),
483                            1, d3uvv1.Crossed(d1v1));
484      DPrim = ncrossns1.Dot(nplan.Crossed(dns1v1));
485      smallterm =  - 2*DPrim*cube;
486      DSecn = ncrossns1.Dot(nplan.Crossed(d2ns1v1))
487           + (nplan.Crossed(dns1v1)).SquareMagnitude();
488      grosterme  = (3*DPrim*DPrim*carre - DSecn) * cube;
489        
490      p1 =  nplan.Dot(dns1v1);
491      p2 =  nplan.Dot(d2ns1v1);
492      temp.SetLinearForm( grosterme*ndotns1, nplan,
493                         -grosterme        , nsurf1);
494      d2ndv1.SetLinearForm( invnorm1*p2 
495                           +smallterm*p1, nplan,
496                           -smallterm,    dns1v1,
497                           -invnorm1,     d2ns1v1);
498      d2ndv1 += temp;
499      resul.SetLinearForm(ray1, d2ndv1, d2v1);
500
501      D2EDX2(2,2,2) = resul.X();
502      D2EDX2(3,2,2) = resul.Y();
503      D2EDX2(4,2,2) = resul.Z();
504                                // ================
505                                // ==  Surface 2 ==
506                                // ================
507      carre = invnorm2*invnorm2;
508      cube  = carre*invnorm2;
509      // Derived double compared to u2       
510        // Derived from the norm
511      d2ns1u2.SetLinearForm(1, d3u2.Crossed(d1v2),
512                            2, d2u2.Crossed(d2uv2),
513                            1, d1u2.Crossed(d3uuv2));
514      DPrim = ncrossns2.Dot(nplan.Crossed(dns1u2));
515      smallterm =  - 2*DPrim*cube;
516      DSecn = ncrossns2.Dot(nplan.Crossed(d2ns1u2))
517           + (nplan.Crossed(dns1u2)).SquareMagnitude();
518      grosterme  = (3*DPrim*DPrim*carre - DSecn) * cube;     
519
520      temp.SetLinearForm( grosterme*ndotns2, nplan,
521                         -grosterme        , nsurf2);     
522      p1 =  nplan.Dot(dns1u2);
523      p2 =  nplan.Dot(d2ns1u2);    
524      d2ndu2.SetLinearForm( invnorm2*p2 
525                           +smallterm*p1, nplan,
526                           -smallterm,    dns1u2,
527                           -invnorm2,     d2ns1u2);
528      d2ndu2 += temp;
529      resul.SetLinearForm(-ray2, d2ndu2, -1, d2u2);
530      D2EDX2(2,3,3) = resul.X();
531      D2EDX2(3,3,3) = resul.Y();
532      D2EDX2(4,3,3) = resul.Z();
533
534      // Derived double compared to u2, v2       
535        // Derived from the norm
536      d2ns1uv2 =  (d3uuv2.Crossed(d1v2))
537               +  (d2u2  .Crossed(d2v2))
538               +  (d1u2  .Crossed(d3uvv2));
539      uterm =  ncrossns2.Dot(nplan.Crossed(dns1u2));
540      vterm =  ncrossns2.Dot(nplan.Crossed(dns1v2));
541      DSecn = (nplan.Crossed(dns1v2)).Dot(nplan.Crossed(dns1u2))
542            +  ncrossns2.Dot(nplan.Crossed(d2ns1uv2));
543      grosterme  = (3*uterm*vterm*carre-DSecn)*cube;
544      uterm *= -cube; //and only now
545      vterm *= -cube;
546        
547      p1 = nplan.Dot(dns1u2);
548      p2 = nplan.Dot(dns1v2);
549      temp.SetLinearForm(  grosterme*ndotns2, nplan,
550                         - grosterme, nsurf2,
551                         - invnorm2, d2ns1uv2);
552      d2nduv2.SetLinearForm( invnorm2*nplan.Dot(d2ns1uv2)
553                           + uterm*p2 
554                           + vterm*p1, nplan,
555                           - uterm,    dns1v2,
556                           - vterm,    dns1u2);
557                         
558      d2nduv2 += temp;
559      resul.SetLinearForm(-ray2, d2nduv2, -1, d2uv2);
560
561      D2EDX2(2,4,3) = D2EDX2(2,3,4) = resul.X();
562      D2EDX2(3,4,3) = D2EDX2(3,3,4) = resul.Y();
563      D2EDX2(4,4,3) = D2EDX2(4,3,4) = resul.Z();    
564
565      // Derived double compared to v2       
566        // Derived from the norm
567      d2ns1v2.SetLinearForm(1, d1u2.Crossed(d3v2),
568                            2, d2uv2.Crossed(d2v2),
569                            1, d3uvv2.Crossed(d1v2));
570      DPrim = ncrossns2.Dot(nplan.Crossed(dns1v2));
571      smallterm =  - 2*DPrim*cube;
572      DSecn = ncrossns2.Dot(nplan.Crossed(d2ns1v2))
573           + (nplan.Crossed(dns1v2)).SquareMagnitude();
574      grosterme  = (3*DPrim*DPrim*carre - DSecn) * cube;
575        
576      p1 =  nplan.Dot(dns1v2);
577      p2 =  nplan.Dot(d2ns1v2);
578      temp.SetLinearForm( grosterme*ndotns2, nplan,
579                         -grosterme        , nsurf2);
580      d2ndv2.SetLinearForm( invnorm2*p2 
581                           +smallterm*p1, nplan,
582                           -smallterm,    dns1v2,
583                           -invnorm2,     d2ns1v2);
584      d2ndv2 += temp;
585      resul.SetLinearForm(-ray2, d2ndv2, -1, d2v2);
586
587      D2EDX2(2,4,4) = resul.X();
588      D2EDX2(3,4,4) = resul.Y();
589      D2EDX2(4,4,4) = resul.Z();
590
591      if (byParam) {
592        Standard_Real tterm;
593         //  ---------- Derivation double in t, X --------------------------
594        D2EDXDT(1,1) = dnplan.Dot(d1u1)/2;
595        D2EDXDT(1,2) = dnplan.Dot(d1v1)/2;
596        D2EDXDT(1,3) = dnplan.Dot(d1u2)/2;
597        D2EDXDT(1,4) = dnplan.Dot(d1v2)/2;
598
599        carre = invnorm1*invnorm1;
600        cube = carre*invnorm1;
601        //--> Derived compared to u1 and t
602        tterm =  ncrossns1.Dot(dnplan.Crossed(nsurf1));
603        smallterm  = - tterm*cube;
604        // Derived from the norm
605        uterm =  ncrossns1.Dot(nplan. Crossed(dns1u1));
606        DSecn = (nplan.Crossed(dns1u1)).Dot(dnplan.Crossed(nsurf1))
607               + ncrossns1.Dot(dnplan.Crossed(dns1u1));
608        grosterme = (3*uterm*tterm*carre - DSecn) * cube;
609        uterm *= -cube;
610
611        p1 = dnplan.Dot(nsurf1);
612        p2 = nplan. Dot(dns1u1);
613        p12 = dnplan.Dot(dns1u1);  
614      
615        d2ndtu1.SetLinearForm( invnorm1*p12 
616                             + smallterm*p2
617                             + uterm*p1
618                             + grosterme*ndotns1, nplan,
619                               invnorm1*p2 
620                             + uterm*ndotns1,     dnplan,
621                             - smallterm,         dns1u1);
622        d2ndtu1 -= grosterme*nsurf1;
623
624        resul = ray1*d2ndtu1;
625        D2EDXDT(2,1) = resul.X();
626        D2EDXDT(3,1) = resul.Y();
627        D2EDXDT(4,1) = resul.Z();
628
629        //--> Derived compared to v1 and t
630        // Derived from the norm
631        uterm =  ncrossns1.Dot(nplan. Crossed(dns1v1));
632        DSecn = (nplan. Crossed(dns1v1)).Dot(dnplan.Crossed(nsurf1))
633               + ncrossns1.Dot(dnplan.Crossed(dns1v1));
634        grosterme = (3*uterm*tterm*carre - DSecn) * cube;
635        uterm *= -cube;
636
637        p1 =  dnplan.Dot(nsurf1);
638        p2 =  nplan. Dot(dns1v1);
639        p12 = dnplan.Dot(dns1v1);       
640        d2ndtv1.SetLinearForm( invnorm1*p12 
641                            + uterm*p1 
642                            + smallterm*p2
643                            + grosterme*ndotns1, nplan,
644                              invnorm1*p2 
645                            + uterm*ndotns1,    dnplan,
646                            - smallterm     ,   dns1v1);
647        d2ndtv1 -= grosterme*nsurf1;
648
649        resul = ray1* d2ndtv1;
650        D2EDXDT(2,2) = resul.X();
651        D2EDXDT(3,2) = resul.Y();
652        D2EDXDT(4,2) = resul.Z();     
653
654        carre = invnorm2*invnorm2;
655        cube = carre*invnorm2;
656        //--> Derived compared to u2 and t
657        tterm =  ncrossns2.Dot(dnplan.Crossed(nsurf2));
658        smallterm = -tterm*cube;
659        // Derived from the norm
660        uterm =  ncrossns2.Dot(nplan. Crossed(dns1u2));
661        DSecn = (nplan. Crossed(dns1u2)).Dot(dnplan.Crossed(nsurf2))
662              + ncrossns2.Dot(dnplan.Crossed(dns1u2));
663        grosterme = (3*uterm*tterm*carre - DSecn) * cube;
664        uterm *= -cube;
665
666        p1 =  dnplan.Dot(nsurf2);
667        p2 =  nplan. Dot(dns1u2);
668        p12 = dnplan.Dot(dns1u2);
669        
670        d2ndtu2.SetLinearForm( invnorm2*p12
671                            + smallterm*p2
672                            + uterm*p1
673                            + grosterme*ndotns2, nplan,
674                              invnorm2*p2
675                            + uterm*ndotns2,    dnplan,
676                             -smallterm     ,   dns1u2);
677        d2ndtu2 -= grosterme*nsurf2;
678
679        resul = - ray2*d2ndtu2;
680        D2EDXDT(2,3) = resul.X();
681        D2EDXDT(3,3) = resul.Y();
682        D2EDXDT(4,3) = resul.Z();
683
684        //--> Derived compared to v2 and t
685        // Derived from the norm
686        uterm =  ncrossns2.Dot(nplan. Crossed(dns1v2));
687        DSecn = (nplan.Crossed(dns1v2)).Dot(dnplan.Crossed(nsurf2))
688              +  ncrossns2.Dot(dnplan.Crossed(dns1v2)); 
689        grosterme = (3*uterm*tterm*carre - DSecn) * cube;
690        uterm *= - cube;
691
692        p1 =  dnplan.Dot(nsurf2);
693        p2 =  nplan. Dot(dns1v2);
694        p12 = dnplan.Dot(dns1v2); 
695       
696        d2ndtv2.SetLinearForm( invnorm2*p12 
697                             + smallterm*p2
698                             + uterm*p1
699                             + grosterme*ndotns2, nplan,
700                               invnorm2*p2
701                             + uterm*ndotns2,     dnplan,
702                              -smallterm        , dns1v2);
703        d2ndtv2 -= grosterme*nsurf2;
704
705        resul = - ray2 * d2ndtv2;
706        D2EDXDT(2,4) = resul.X();
707        D2EDXDT(3,4) = resul.Y();
708        D2EDXDT(4,4) = resul.Z();     
709
710
711         //  ---------- Derivation double in t -----------------------------
712     // Derived from n1 compared to w
713        carre = invnorm1*invnorm1;
714        cube = carre*invnorm1;
715        // Derived from the norm
716        DPrim =  ncrossns1.Dot(dnplan.Crossed(nsurf1));
717        smallterm = - 2*DPrim*cube;
718        DSecn     = (dnplan.Crossed(nsurf1)).SquareMagnitude()
719                   + ncrossns1.Dot(d2nplan.Crossed(nsurf1));
720        grosterme = (3*DPrim*DPrim*carre - DSecn) * cube; 
721       
722        p1 =  dnplan. Dot(nsurf1);
723        p2 =  d2nplan.Dot(nsurf1);
724
725        temp.SetLinearForm( grosterme*ndotns1, nplan,
726                           -grosterme , nsurf1);
727        d2n1w.SetLinearForm(  smallterm*p1
728                            + invnorm1*p2,      nplan,
729                              smallterm*ndotns1
730                            + 2*invnorm1*p1,    dnplan,
731                              ndotns1*invnorm1,   d2nplan);
732        d2n1w += temp;
733   
734      // Derived from n2 compared to w
735        carre = invnorm2*invnorm2;
736        cube = carre*invnorm2;
737        // Derived from the norm
738        DPrim =  ncrossns2.Dot(dnplan.Crossed(nsurf2));
739        smallterm = - 2*DPrim*cube;
740        DSecn     = (dnplan.Crossed(nsurf2)).SquareMagnitude()
741                   + ncrossns2.Dot(d2nplan.Crossed(nsurf2));
742        grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
743    
744        p1 =  dnplan. Dot(nsurf2);
745        p2 =  d2nplan.Dot(nsurf2);
746
747        temp.SetLinearForm( grosterme*ndotns2, nplan,
748                           -grosterme ,        nsurf2);
749        d2n2w.SetLinearForm(  smallterm*p1
750                            + invnorm2*p2,        nplan,
751                              smallterm*ndotns2
752                            + 2*invnorm2*p1,     dnplan,
753                              ndotns2*invnorm2,  d2nplan);
754        d2n2w += temp;
755
756        temp.SetXYZ( (pts1.XYZ()+pts2.XYZ())/2 - ptgui.XYZ());
757        D2EDT2(1) = d2nplan.Dot(temp) - 2*dnplan.Dot(d1gui) - nplan.Dot(d2gui);
758        D2EDT2(2) = ray1*d2n1w.X() - ray2*d2n2w.X();
759        D2EDT2(3) = ray1*d2n1w.Y() - ray2*d2n2w.Y();
760        D2EDT2(4) = ray1*d2n1w.Z() - ray2*d2n2w.Z();
761      }
762    }
763  }
764  return Standard_True; 
765 }
766
767 //=======================================================================
768 //function : Set
769 //purpose  : 
770 //=======================================================================
771
772 void BlendFunc_ConstRad::Set(const Standard_Real Param)
773 {
774   param = Param;
775 }
776
777 //=======================================================================
778 //function : Set
779 //purpose  : Segmentation of the useful part of the curve
780 //           Precision is taken at random and small !?
781 //=======================================================================
782
783 void BlendFunc_ConstRad::Set(const Standard_Real First, const Standard_Real Last)
784 {
785   tcurv = curv->Trim(First, Last, 1.e-12);
786 }
787
788 //=======================================================================
789 //function : GetTolerance
790 //purpose  : 
791 //=======================================================================
792
793 void BlendFunc_ConstRad::GetTolerance(math_Vector& Tolerance, const Standard_Real Tol) const
794 {
795   Tolerance(1) = surf1->UResolution(Tol);
796   Tolerance(2) = surf1->VResolution(Tol);
797   Tolerance(3) = surf2->UResolution(Tol);
798   Tolerance(4) = surf2->VResolution(Tol);
799 }
800
801
802 //=======================================================================
803 //function : GetBounds
804 //purpose  : 
805 //=======================================================================
806
807 void BlendFunc_ConstRad::GetBounds(math_Vector& InfBound, math_Vector& SupBound) const
808 {
809   InfBound(1) = surf1->FirstUParameter();
810   InfBound(2) = surf1->FirstVParameter();
811   InfBound(3) = surf2->FirstUParameter();
812   InfBound(4) = surf2->FirstVParameter();
813   SupBound(1) = surf1->LastUParameter();
814   SupBound(2) = surf1->LastVParameter();
815   SupBound(3) = surf2->LastUParameter();
816   SupBound(4) = surf2->LastVParameter();
817
818   for(Standard_Integer i = 1; i <= 4; i++){
819     if(!Precision::IsInfinite(InfBound(i)) &&
820        !Precision::IsInfinite(SupBound(i))) {
821       Standard_Real range = (SupBound(i) - InfBound(i));
822       InfBound(i) -= range;
823       SupBound(i) += range;
824     }
825   }
826 }
827
828
829 //=======================================================================
830 //function : IsSolution
831 //purpose  : 
832 //=======================================================================
833
834 Standard_Boolean BlendFunc_ConstRad::IsSolution(const math_Vector& Sol, const Standard_Real Tol)
835 {
836   Standard_Real norm, Cosa, Sina, Angle;
837   Standard_Boolean Ok=Standard_True;
838
839   Ok = ComputeValues(Sol, 1, Standard_True, param);
840
841   if (Abs(E(1)) <= Tol &&
842       E(2)*E(2) + E(3)*E(3) +  E(4)*E(4) <= Tol*Tol) { 
843
844     // ns1, ns2 and  np are copied locally to avoid crushing the fields !
845     gp_Vec ns1,ns2,np;
846     ns1 = nsurf1;
847     ns2 = nsurf2;
848     np = nplan;
849   
850     norm = nplan.Crossed(ns1).Magnitude();
851     if (norm < Eps)  {
852       norm = 1; // Unsatisfactory, but it is not necessary to stop
853     }
854     ns1.SetLinearForm(nplan.Dot(ns1)/norm,nplan, -1./norm, ns1);
855
856     norm = nplan.Crossed(ns2).Magnitude();
857     if (norm < Eps)  {
858       norm = 1; // Unsatisfactory, but it is not necessary to stop
859     }
860     ns2.SetLinearForm(nplan.Dot(ns2)/norm,nplan, -1./norm, ns2);   
861
862     Standard_Real maxpiv = 1.e-9;
863     math_Vector controle(1,4),solution(1,4), tolerances(1,4);
864     GetTolerance(tolerances,Tol);
865
866     istangent = Standard_True;
867     math_Gauss Resol(DEDX,maxpiv);
868     if (Resol.IsDone()) {
869       Resol.Solve(-DEDT,solution);
870       istangent = Standard_False;
871       controle = DEDT.Added(DEDX.Multiplied(solution));
872       if(Abs(controle(1)) > tolerances(1) ||
873          Abs(controle(2)) > tolerances(2) ||
874          Abs(controle(3)) > tolerances(3) ||
875          Abs(controle(4)) > tolerances(4)){
876         istangent = Standard_True;
877       } 
878     }
879
880     if (istangent) {
881       math_SVD SingRS (DEDX);
882       if (SingRS.IsDone()) {
883         SingRS.Solve(-DEDT, solution, 1.e-6);
884         istangent = Standard_False;
885         controle = DEDT.Added(DEDX.Multiplied(solution));
886         if(Abs(controle(1)) > tolerances(1) ||
887            Abs(controle(2)) > tolerances(2) ||
888            Abs(controle(3)) > tolerances(3) ||
889            Abs(controle(4)) > tolerances(4)){
890 #ifdef OCCT_DEBUG
891         std::cout<<"Cheminement : echec calcul des derivees"<<std::endl;
892 #endif
893           istangent = Standard_True;
894         } 
895       }
896     }
897
898     if(!istangent){
899       tg1.SetLinearForm(solution(1),d1u1,solution(2),d1v1);
900       tg2.SetLinearForm(solution(3),d1u2,solution(4),d1v2);
901       tg12d.SetCoord(solution(1),solution(2));
902       tg22d.SetCoord(solution(3),solution(4));
903     }
904  
905   // update of maxang
906
907     if (ray1 > 0.) {
908       ns1.Reverse();
909     }
910     if (ray2 >0.) {
911       ns2.Reverse();
912     }
913     Cosa = ns1.Dot(ns2);
914     Sina = np.Dot(ns1.Crossed(ns2));
915     if (choix%2 != 0) {
916       Sina = -Sina;  //nplan is changed in -nplan
917     }
918
919     if(Cosa > 1.) {Cosa = 1.; Sina = 0.;}
920     Angle = ACos(Cosa);
921
922  // Reframing on  ]-pi/2, 3pi/2]
923     if (Sina <0.) {
924       if (Cosa > 0.) Angle = -Angle;
925       else           Angle =  2.*M_PI - Angle;
926     }
927
928 //    std::cout << "Angle : " <<Angle << std::endl;
929 //    if ((Angle < 0) || (Angle > M_PI)) {
930 //      std::cout << "t = " << param << std::endl;
931 //    }
932
933     if (Abs(Angle)>maxang) {maxang = Abs(Angle);}
934     if (Abs(Angle)<minang) {minang = Abs(Angle);}
935     distmin = Min( distmin, pts1.Distance(pts2));
936
937     return Ok;
938   }
939   istangent = Standard_True;
940   return Standard_False;
941 }
942
943 //=======================================================================
944 //function : GetMinimalDistance
945 //purpose  : 
946 //=======================================================================
947
948 Standard_Real BlendFunc_ConstRad::GetMinimalDistance() const
949 {
950   return distmin;
951 }
952
953 //=======================================================================
954 //function : Value
955 //purpose  : 
956 //=======================================================================
957
958 Standard_Boolean BlendFunc_ConstRad::Value(const math_Vector& X, math_Vector& F)
959 {
960   const Standard_Boolean Ok = ComputeValues(X, 0);
961   F = E;
962   return Ok;
963 }
964
965
966
967 //=======================================================================
968 //function : Derivatives
969 //purpose  : 
970 //=======================================================================
971
972 Standard_Boolean BlendFunc_ConstRad::Derivatives(const math_Vector& X, math_Matrix& D)
973 {
974   const Standard_Boolean Ok = ComputeValues(X, 1);
975   D = DEDX;
976   return Ok;
977 }
978
979 //=======================================================================
980 //function : Values
981 //purpose  : 
982 //=======================================================================
983
984 Standard_Boolean BlendFunc_ConstRad::Values(const math_Vector& X, math_Vector& F, math_Matrix& D)
985 {
986   const Standard_Boolean Ok = ComputeValues(X, 1);
987   F = E;
988   D = DEDX;
989   return Ok;
990 }  
991
992 //=======================================================================
993 //function : PointOnS1
994 //purpose  : 
995 //=======================================================================
996
997 const gp_Pnt& BlendFunc_ConstRad::PointOnS1 () const
998 {
999   return pts1;
1000 }
1001
1002
1003 //=======================================================================
1004 //function : PointOnS2
1005 //purpose  : 
1006 //=======================================================================
1007
1008 const gp_Pnt& BlendFunc_ConstRad::PointOnS2 () const
1009 {
1010   return pts2;
1011 }
1012
1013
1014 //=======================================================================
1015 //function : IsTangencyPoint
1016 //purpose  : 
1017 //=======================================================================
1018
1019 Standard_Boolean BlendFunc_ConstRad::IsTangencyPoint () const
1020 {
1021   return istangent;
1022 }
1023
1024
1025 //=======================================================================
1026 //function : TangentOnS1
1027 //purpose  : 
1028 //=======================================================================
1029
1030 const gp_Vec& BlendFunc_ConstRad::TangentOnS1 () const
1031 {
1032   if (istangent)
1033     throw Standard_DomainError("BlendFunc_ConstRad::TangentOnS1");
1034   return tg1;
1035 }
1036
1037
1038 //=======================================================================
1039 //function : TangentOnS2
1040 //purpose  : 
1041 //=======================================================================
1042
1043 const gp_Vec& BlendFunc_ConstRad::TangentOnS2 () const
1044 {
1045   if (istangent)
1046     throw Standard_DomainError("BlendFunc_ConstRad::TangentOnS2");
1047   return tg2;
1048 }
1049
1050
1051 //=======================================================================
1052 //function : Tangent2dOnS1
1053 //purpose  : 
1054 //=======================================================================
1055
1056 const gp_Vec2d& BlendFunc_ConstRad::Tangent2dOnS1 () const
1057 {
1058   if (istangent)
1059     throw Standard_DomainError("BlendFunc_ConstRad::Tangent2dOnS1");
1060   return tg12d;
1061 }
1062
1063
1064 //=======================================================================
1065 //function : Tangent2dOnS2
1066 //purpose  : 
1067 //=======================================================================
1068
1069 const gp_Vec2d& BlendFunc_ConstRad::Tangent2dOnS2 () const
1070 {
1071   if (istangent)
1072     throw Standard_DomainError("BlendFunc_ConstRad::Tangent2dOnS2");
1073   return tg22d;
1074 }
1075
1076
1077 //=======================================================================
1078 //function : Tangent
1079 //purpose  : 
1080 //=======================================================================
1081
1082 void BlendFunc_ConstRad::Tangent(const Standard_Real U1,
1083                                  const Standard_Real V1,
1084                                  const Standard_Real U2,
1085                                  const Standard_Real V2,
1086                                  gp_Vec& TgF,
1087                                  gp_Vec& TgL,
1088                                  gp_Vec& NmF,
1089                                  gp_Vec& NmL) const
1090 {
1091   gp_Pnt Center;
1092   gp_Vec ns1;
1093   Standard_Real  invnorm1;
1094
1095   if ((U1!=xval(1)) || (V1!=xval(2)) || 
1096       (U2!=xval(3)) || (V2!=xval(4))) {
1097     gp_Vec d1u,d1v;
1098     gp_Pnt bid;
1099     surf1->D1(U1,V1,bid,d1u,d1v);  
1100     NmF = ns1 = d1u.Crossed(d1v);
1101     surf2->D1(U2,V2,bid,d1u,d1v);
1102     NmL = d1u.Crossed(d1v);
1103   }
1104   else {
1105     NmF = ns1 = nsurf1;
1106     NmL = nsurf2;
1107   }
1108
1109   invnorm1 = nplan.Crossed(ns1).Magnitude();
1110   if (invnorm1 < Eps) invnorm1 = 1;
1111   else invnorm1 = 1. / invnorm1;
1112
1113   ns1.SetLinearForm(nplan.Dot(ns1)*invnorm1,nplan, -invnorm1,ns1);
1114   Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1115
1116   TgF = nplan.Crossed(gp_Vec(Center,pts1));
1117   TgL = nplan.Crossed(gp_Vec(Center,pts2));
1118   if (choix%2 == 1) {
1119     TgF.Reverse();
1120     TgL.Reverse();
1121   }
1122 }
1123
1124 //=======================================================================
1125 //function : TwistOnS1
1126 //purpose  : 
1127 //=======================================================================
1128
1129 Standard_Boolean BlendFunc_ConstRad::TwistOnS1() const
1130 {
1131   if (istangent)
1132     throw Standard_DomainError("BlendFunc_ConstRad::TwistOnS1");
1133   return tg1.Dot(nplan) < 0.;
1134 }
1135
1136 //=======================================================================
1137 //function : TwistOnS2
1138 //purpose  : 
1139 //=======================================================================
1140
1141 Standard_Boolean BlendFunc_ConstRad::TwistOnS2() const
1142 {
1143   if (istangent)
1144     throw Standard_DomainError("BlendFunc_ConstRad::TwistOnS2");
1145   return tg2.Dot(nplan) < 0.;
1146 }
1147
1148 //=======================================================================
1149 //function : Section
1150 //purpose  : 
1151 //=======================================================================
1152
1153 void BlendFunc_ConstRad::Section(const Standard_Real Param,
1154                                  const Standard_Real U1,
1155                                  const Standard_Real V1,
1156                                  const Standard_Real U2,
1157                                  const Standard_Real V2,
1158                                  Standard_Real& Pdeb,
1159                                  Standard_Real& Pfin,
1160                                  gp_Circ& C)
1161 {
1162   gp_Pnt Center;
1163   gp_Vec ns1,np;
1164
1165   math_Vector X(1,4);
1166   X(1) = U1; X(2) = V1; X(3) = U2; X(4) = V2; 
1167   Standard_Real prm = Param;
1168
1169   ComputeValues(X, 0, Standard_True, prm);
1170
1171   ns1 = nsurf1;
1172   np = nplan;
1173
1174   Standard_Real  norm1;
1175   norm1 = nplan.Crossed(ns1).Magnitude();
1176   if (norm1 < Eps)  {
1177     norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1178   }   
1179   ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1180   Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1181
1182 // ns1 is oriented from the center to pts1,
1183
1184   if (ray1 > 0.) {
1185     ns1.Reverse();
1186   }
1187   if (choix%2 != 0) {
1188     np.Reverse();
1189   }
1190   C.SetRadius(Abs(ray1));
1191   C.SetPosition(gp_Ax2(Center,np,ns1));
1192   Pdeb = 0.;
1193   Pfin = ElCLib::Parameter(C,pts2);
1194   // Test negative and almost null angles : Singular Case
1195   if (Pfin>1.5*M_PI) {
1196     np.Reverse();
1197     C.SetPosition(gp_Ax2(Center,np,ns1));
1198     Pfin = ElCLib::Parameter(C,pts2);
1199   }
1200   if (Pfin < Precision::PConfusion()) Pfin += Precision::PConfusion();
1201 }
1202
1203 //=======================================================================
1204 //function : IsRational
1205 //purpose  : 
1206 //=======================================================================
1207
1208 Standard_Boolean BlendFunc_ConstRad::IsRational () const
1209 {
1210  return  (mySShape==BlendFunc_Rational || mySShape==BlendFunc_QuasiAngular);
1211 }
1212
1213 //=======================================================================
1214 //function : GetSectionSize
1215 //purpose  : 
1216 //=======================================================================
1217 Standard_Real BlendFunc_ConstRad::GetSectionSize() const 
1218 {
1219   return maxang*Abs(ray1);
1220 }
1221
1222 //=======================================================================
1223 //function : GetMinimalWeight
1224 //purpose  : 
1225 //=======================================================================
1226 void BlendFunc_ConstRad::GetMinimalWeight(TColStd_Array1OfReal& Weigths) const 
1227 {
1228   BlendFunc::GetMinimalWeights(mySShape, myTConv, minang, maxang, Weigths );
1229   // It is supposed that it does not depend on the Radius!
1230 }
1231
1232 //=======================================================================
1233 //function : NbIntervals
1234 //purpose  : 
1235 //=======================================================================
1236
1237 Standard_Integer BlendFunc_ConstRad::NbIntervals (const GeomAbs_Shape S) const
1238 {
1239   return curv->NbIntervals(BlendFunc::NextShape(S));
1240 }
1241
1242
1243 //=======================================================================
1244 //function : Intervals
1245 //purpose  : 
1246 //=======================================================================
1247
1248 void BlendFunc_ConstRad::Intervals (TColStd_Array1OfReal& T,
1249                                     const GeomAbs_Shape S) const
1250 {
1251   curv->Intervals(T, BlendFunc::NextShape(S));
1252 }
1253
1254 //=======================================================================
1255 //function : GetShape
1256 //purpose  : 
1257 //=======================================================================
1258
1259 void BlendFunc_ConstRad::GetShape (Standard_Integer& NbPoles,
1260                                    Standard_Integer& NbKnots,
1261                                    Standard_Integer& Degree,
1262                                    Standard_Integer& NbPoles2d)
1263 {
1264   NbPoles2d = 2;
1265   BlendFunc::GetShape(mySShape,maxang,NbPoles,NbKnots,Degree,myTConv);
1266 }
1267
1268 //=======================================================================
1269 //function : GetTolerance
1270 //purpose  : Determine Tolerances used for approximations.
1271 //=======================================================================
1272 void BlendFunc_ConstRad::GetTolerance(const Standard_Real BoundTol, 
1273                                       const Standard_Real SurfTol, 
1274                                       const Standard_Real AngleTol, 
1275                                       math_Vector& Tol3d, 
1276                                       math_Vector& Tol1d) const
1277 {
1278  Standard_Integer low = Tol3d.Lower() , up=Tol3d.Upper();
1279  Standard_Real Tol;
1280  Tol= GeomFill::GetTolerance(myTConv, minang, Abs(ray1),
1281                              AngleTol, SurfTol);
1282  Tol1d.Init(SurfTol);
1283  Tol3d.Init(SurfTol);
1284  Tol3d(low+1) = Tol3d(up-1) = Min(Tol, SurfTol);
1285  Tol3d(low) = Tol3d(up) = Min(Tol, BoundTol);
1286     
1287 }
1288
1289 //=======================================================================
1290 //function : Knots
1291 //purpose  : 
1292 //=======================================================================
1293
1294 void BlendFunc_ConstRad::Knots(TColStd_Array1OfReal& TKnots)
1295 {
1296   GeomFill::Knots(myTConv,TKnots);
1297 }
1298
1299
1300 //=======================================================================
1301 //function : Mults
1302 //purpose  : 
1303 //=======================================================================
1304
1305 void BlendFunc_ConstRad::Mults(TColStd_Array1OfInteger& TMults)
1306 {
1307   GeomFill::Mults(myTConv,TMults);
1308 }
1309
1310
1311 //=======================================================================
1312 //function : Section
1313 //purpose  : 
1314 //=======================================================================
1315
1316 void BlendFunc_ConstRad::Section(const Blend_Point& P,
1317                                  TColgp_Array1OfPnt& Poles,
1318                                  TColgp_Array1OfPnt2d& Poles2d,
1319                                  TColStd_Array1OfReal& Weights)
1320 {
1321   gp_Pnt Center;
1322   gp_Vec ns1,ns2,np;
1323
1324   math_Vector X(1,4);
1325   Standard_Real prm = P.Parameter();
1326
1327   Standard_Integer low = Poles.Lower();
1328   Standard_Integer upp = Poles.Upper();
1329
1330   P.ParametersOnS1(X(1), X(2));
1331   P.ParametersOnS2(X(3), X(4));
1332
1333   ComputeValues(X, 0, Standard_True, prm);
1334   distmin = Min (distmin, pts1.Distance(pts2));
1335
1336   // ns1, ns2, np are copied locally to avoid crushing the fields !
1337   ns1 = nsurf1;
1338   ns2 = nsurf2;
1339   np = nplan;
1340
1341   
1342   Poles2d(Poles2d.Lower()).SetCoord(X(1), X(2));
1343   Poles2d(Poles2d.Upper()).SetCoord(X(3), X(4));
1344
1345   if (mySShape == BlendFunc_Linear) {
1346     Poles(low) = pts1;
1347     Poles(upp) = pts2;
1348     Weights(low) = 1.0;
1349     Weights(upp) = 1.0;
1350     return;
1351   }
1352
1353   Standard_Real  norm1, norm2;
1354   norm1 = nplan.Crossed(ns1).Magnitude();
1355   norm2 = nplan.Crossed(ns2).Magnitude();
1356   if (norm1 < Eps)  {
1357     norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1358   }
1359   if (norm2 < Eps) {
1360     norm2 = 1; // Unsatisfactory, but it is not necessary to stop
1361   }
1362
1363   ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1364   ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2,ns2);
1365
1366   Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1367
1368 // ns1 (resp. ns2) is oriented from center to pts1 (resp. pts2),
1369 // and the triedron ns1,ns2,nplan is made direct.
1370
1371   if (ray1 > 0.) {
1372     ns1.Reverse();
1373   }
1374   if (ray2 >0.) {
1375     ns2.Reverse();
1376   }
1377   if (choix%2 != 0) {
1378     np.Reverse();
1379   }
1380
1381   GeomFill::GetCircle(myTConv,
1382                       ns1, ns2, 
1383                       np, pts1, pts2,
1384                       Abs(ray1), Center, 
1385                       Poles, Weights);
1386 }
1387  
1388
1389 //=======================================================================
1390 //function : Section
1391 //purpose  : 
1392 //=======================================================================
1393
1394 Standard_Boolean BlendFunc_ConstRad::Section
1395   (const Blend_Point& P,
1396    TColgp_Array1OfPnt& Poles,
1397    TColgp_Array1OfVec& DPoles,
1398    TColgp_Array1OfPnt2d& Poles2d,
1399    TColgp_Array1OfVec2d& DPoles2d,
1400    TColStd_Array1OfReal& Weights,
1401    TColStd_Array1OfReal& DWeights)
1402 {
1403   gp_Vec ns1, ns2, np, dnp, dnorm1w, dnorm2w, tgc;
1404   Standard_Real norm1, norm2;
1405
1406   gp_Pnt Center;
1407   math_Vector sol(1,4), secmember(1,4);
1408
1409   Standard_Real prm = P.Parameter();
1410   Standard_Integer low = Poles.Lower();
1411   Standard_Integer upp = Poles.Upper();
1412   Standard_Boolean istgt = Standard_True;
1413
1414   P.ParametersOnS1(sol(1),sol(2));
1415   P.ParametersOnS2(sol(3),sol(4));
1416
1417   // Calculation of equations
1418   ComputeValues(sol, 1, Standard_True, prm);
1419   distmin = Min (distmin, pts1.Distance(pts2));
1420
1421   // ns1, ns2, np are copied locally to avoid crushing the fields !
1422   ns1 = nsurf1;
1423   ns2 = nsurf2;
1424   np = nplan;
1425   dnp = dnplan;
1426
1427   if ( ! pts1.IsEqual(pts2, 1.e-4)) {
1428
1429     // Calculation of derivates Processing Normal
1430     math_Gauss Resol(DEDX, 1.e-9);
1431
1432     if (Resol.IsDone()) {  
1433       Resol.Solve(-DEDT, secmember);
1434       istgt = Standard_False;
1435     }
1436   }
1437
1438  if (istgt) {
1439     math_SVD SingRS (DEDX);
1440     if (SingRS.IsDone()) {
1441       SingRS.Solve(-DEDT, secmember, 1.e-6);
1442       istgt = Standard_False;
1443       }
1444   }
1445
1446  if (!istgt) {
1447     tg1.SetLinearForm(secmember(1),d1u1, secmember(2),d1v1);
1448     tg2.SetLinearForm(secmember(3),d1u2, secmember(4),d1v2);
1449
1450     dnorm1w.SetLinearForm(secmember(1),dndu1, secmember(2),dndv1, dn1w);
1451     dnorm2w.SetLinearForm(secmember(3),dndu2, secmember(4),dndv2, dn2w);
1452   }
1453
1454
1455   // Tops 2d  
1456   Poles2d(Poles2d.Lower()).SetCoord(sol(1),sol(2));
1457   Poles2d(Poles2d.Upper()).SetCoord(sol(3),sol(4));
1458   if (!istgt) {
1459     DPoles2d(Poles2d.Lower()).SetCoord(secmember(1),secmember(2));
1460     DPoles2d(Poles2d.Upper()).SetCoord(secmember(3),secmember(4));
1461   }
1462
1463   // the linear case is processed...
1464   if (mySShape == BlendFunc_Linear) {
1465     Poles(low) = pts1;
1466     Poles(upp) = pts2;
1467     Weights(low) = 1.0;
1468     Weights(upp) = 1.0;
1469     if (!istgt) {
1470       DPoles(low) = tg1;
1471       DPoles(upp) = tg2;
1472       DWeights(low) = 0.0;
1473       DWeights(upp) = 0.0;
1474     }
1475     return (!istgt);
1476   }
1477
1478   // Case of the circle
1479   norm1 = nplan.Crossed(ns1).Magnitude();
1480   norm2 = nplan.Crossed(ns2).Magnitude();
1481   if (norm1 < Eps) {
1482     norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1483 #ifdef OCCT_DEBUG
1484     std::cout << " ConstRad : Surface singuliere " << std::endl;
1485 #endif
1486   }
1487   if (norm2 < Eps) {
1488    norm2 = 1; // Unsatisfactory, but it is not necessary to stop
1489 #ifdef OCCT_DEBUG
1490    std::cout << " ConstRad : Surface singuliere " << std::endl;
1491 #endif
1492  }
1493
1494   ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1495   ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2,ns2);
1496
1497   Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1498   if (!istgt) {
1499     tgc.SetLinearForm(ray1,dnorm1w,tg1); //  = tg1.Added(ray1*dn1w);
1500   }
1501
1502   // ns1 is oriented from the center to pts1, and ns2 from the center to pts2
1503   // and the trihedron ns1,ns2,nplan is made direct
1504
1505   if (ray1 > 0.) {
1506     ns1.Reverse();
1507     if (!istgt) {
1508       dnorm1w.Reverse();
1509     }
1510   }
1511   if (ray2 >0.) {
1512     ns2.Reverse();
1513     if (!istgt) {
1514       dnorm2w.Reverse();
1515     }
1516   }
1517   if (choix%2 != 0) {
1518     np.Reverse();
1519     dnp.Reverse();
1520   }
1521   
1522   if (!istgt) {
1523     return GeomFill::GetCircle(myTConv, 
1524                                ns1, ns2, 
1525                                dnorm1w, dnorm2w, 
1526                                np, dnp, 
1527                                pts1, pts2, 
1528                                tg1, tg2, 
1529                                Abs(ray1), 0, 
1530                                Center, tgc, 
1531                                Poles, 
1532                                DPoles,
1533                                Weights, 
1534                                DWeights); 
1535   }
1536   else {
1537     GeomFill::GetCircle(myTConv,
1538                         ns1, ns2, 
1539                         np, pts1, pts2,
1540                         Abs(ray1), Center, 
1541                         Poles, Weights);
1542     return Standard_False;
1543   }
1544 }
1545
1546 //=======================================================================
1547 //function : Section
1548 //purpose  : 
1549 //=======================================================================
1550
1551 Standard_Boolean BlendFunc_ConstRad::Section
1552   (const Blend_Point& P,
1553    TColgp_Array1OfPnt& Poles,
1554    TColgp_Array1OfVec& DPoles,
1555    TColgp_Array1OfVec& D2Poles,
1556    TColgp_Array1OfPnt2d& Poles2d,
1557    TColgp_Array1OfVec2d& DPoles2d,
1558    TColgp_Array1OfVec2d& D2Poles2d,
1559    TColStd_Array1OfReal& Weights,
1560    TColStd_Array1OfReal& DWeights,
1561    TColStd_Array1OfReal& D2Weights)
1562 {
1563   gp_Vec ns1,ns2, np, dnp, d2np, dnorm1w, dnorm2w, d2norm1w, d2norm2w;
1564   gp_Vec tgc, dtgc, dtg1, dtg2, temp, tempbis;
1565   Standard_Real norm1, norm2;
1566
1567   gp_Pnt Center;
1568   math_Vector X(1,4), sol(1,4), secmember(1,4);
1569   math_Matrix D2DXdSdt(1,4,1,4);
1570
1571   Standard_Real prm = P.Parameter();
1572   Standard_Integer low = Poles.Lower();
1573   Standard_Integer upp = Poles.Upper();
1574   Standard_Boolean istgt = Standard_True;
1575
1576   P.ParametersOnS1(X(1), X(2));
1577   P.ParametersOnS2(X(3), X(4));
1578
1579 /*  Pour debuger par des D.F
1580 #ifdef OCCT_DEBUG
1581   Standard_Real deltat = 1.e-7;
1582   if (prm==tcurv->LastParameter()){deltat *= -1;} //Pour les discont
1583   Standard_Real deltaX = 1.e-7;
1584   Standard_Real seuil = 1.e-3;
1585   Standard_Integer ii, jj;
1586   gp_Vec d_plan, d1, d2, pdiff;
1587   math_Matrix M(1,4,1,4), MDiff(1,4,1,4);
1588   math_Matrix Mu1(1,4,1,4), Mv1(1,4,1,4);
1589   math_Matrix Mu2(1,4,1,4), Mv2(1,4,1,4);  
1590   math_Vector V(1,4), VDiff(1,4),dx(1,4);
1591
1592   dx = X;
1593   dx(1)+=deltaX;
1594   ComputeValues(dx, 1, Standard_True, prm );
1595   Mu1 = DEDX;
1596
1597   dx = X;
1598   dx(2)+=deltaX;
1599   ComputeValues(dx, 1, Standard_True, prm);
1600   Mv1 = DEDX;
1601
1602   dx = X;
1603   dx(3)+=deltaX;
1604   ComputeValues(dx, 1, Standard_True, prm  );
1605   Mu2 = DEDX;
1606
1607   dx = X;
1608   dx(4)+=deltaX;
1609   ComputeValues(dx, 1,  Standard_True, prm );
1610   Mv2 = DEDX;
1611
1612   ComputeValues(X, 1, Standard_True, prm+deltat);
1613   M = DEDX;
1614   V = DEDT;
1615   d_plan = dnplan;
1616   d1 = dn1w;
1617   d2 = dn2w;
1618 # endif 
1619 */
1620
1621   // Calculation of equations
1622   ComputeValues(X, 2, Standard_True, prm);
1623   distmin = Min (distmin, pts1.Distance(pts2));
1624
1625 /*
1626 #ifdef OCCT_DEBUG
1627   MDiff = (M - DEDX)*(1/deltat);
1628   VDiff = (V - DEDT)*(1/deltat);
1629
1630   pdiff = (d_plan - dnplan)*(1/deltat);
1631   if ((pdiff-d2nplan).Magnitude() > seuil*(pdiff.Magnitude()+1))
1632     {
1633       std::cout << "d2nplan = (" << d2nplan.X() << ","<< d2nplan.Y() << ","<< d2nplan.Z() << ")"<<std::endl;
1634       std::cout << "Diff fi = (" << pdiff.X() << ","<<  pdiff.Y() << ","<<  pdiff.Z() << ")"<<std::endl;
1635     }
1636
1637   pdiff = (d1 - dn1w)*(1/deltat);
1638   if ((pdiff-d2n1w).Magnitude() > seuil*(pdiff.Magnitude()+1))
1639     {
1640       std::cout << "d2n1w   = (" << d2n1w.X() << ","<< d2n1w.Y() << ","<< d2n1w.Z() << ")"<<std::endl;
1641       std::cout << "Diff fi = (" << pdiff.X() << ","<<  pdiff.Y() << ","<<  pdiff.Z() << ")"<<std::endl;
1642     }
1643   pdiff = (d2 - dn2w)*(1/deltat);
1644   if ((pdiff-d2n2w).Magnitude() > seuil*(pdiff.Magnitude()+1))
1645     {
1646       std::cout << "d2n2w   = (" << d2n2w.X() << ","<< d2n2w.Y() << ","<< d2n2w.Z() << ")"<<std::endl;
1647       std::cout << "Diff fi = (" << pdiff.X() << ","<<  pdiff.Y() << ","<<  pdiff.Z() << ")"<<std::endl;
1648     }
1649  
1650
1651   for ( ii=1; ii<=4; ii++) {
1652     if (Abs(VDiff(ii)-D2EDT2(ii)) > seuil*(Abs(D2EDT2(ii))+1)) 
1653       {
1654           std::cout << "erreur sur D2EDT2 : "<< ii << std::endl;
1655           std::cout << D2EDT2(ii) << " D.F = " << VDiff(ii) << std::endl;
1656         } 
1657         for (jj=1; jj<=4; jj++) {
1658           if (Abs(MDiff(ii,jj)-D2EDXDT(ii, jj)) > 
1659               1.e-3*(Abs(D2EDXDT(ii, jj))+1.e-2)) 
1660               {
1661                 std::cout << "erreur sur D2EDXDT : "<< ii << " , " << jj << std::endl;
1662                 std::cout << D2EDXDT(ii,jj) << " D.F = " << MDiff(ii,jj) << std::endl;
1663               }
1664         }
1665   }
1666 // Test de D2EDX2 en u1
1667   MDiff = (Mu1 - DEDX)/deltaX;
1668   for (ii=1; ii<=4; ii++) {
1669     for (jj=1; jj<=4; jj++) {
1670       if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 1)) > 
1671           seuil*(Abs(D2EDX2(ii, jj, 1))+1)) 
1672         {
1673           std::cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 1 << std::endl;
1674           std::cout << D2EDX2(ii,jj, 1) << " D.F = " << MDiff(ii,jj) << std::endl;
1675         }
1676     }
1677   }
1678
1679 // Test de D2EDX2 en v1
1680   MDiff = (Mv1 - DEDX)/deltaX;
1681   for (ii=1; ii<=4; ii++) {
1682     for (jj=1; jj<=4; jj++) {
1683       if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 2)) > 
1684           seuil*(Abs(D2EDX2(ii, jj, 2))+1)) 
1685         {
1686           std::cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 2 << std::endl;
1687           std::cout << D2EDX2(ii,jj, 2) << " D.F = " << MDiff(ii,jj) << std::endl;
1688         }
1689     }
1690   }
1691 // Test de D2EDX2 en u2
1692   MDiff = (Mu2 - DEDX)/deltaX;
1693   for (ii=1; ii<=4; ii++) {
1694     for (jj=1; jj<=4; jj++) {
1695       if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 3)) > 
1696           seuil*(Abs(D2EDX2(ii, jj, 3))+1)) 
1697         {
1698           std::cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 3 << std::endl;
1699           std::cout << D2EDX2(ii,jj, 3) << " D.F = " << MDiff(ii,jj) << std::endl;
1700         }
1701     }
1702   }
1703
1704 // Test de D2EDX2 en v2
1705   MDiff = (Mv2 - DEDX)/deltaX;
1706   for (ii=1; ii<=4; ii++) {
1707     for (jj=1; jj<=4; jj++) {
1708       if (Abs(MDiff(ii,jj)-D2EDX2(ii, jj, 4)) > 
1709           seuil*(Abs(D2EDX2(ii, jj, 4))+1)) 
1710         {
1711           std::cout << "erreur sur D2EDX2 : "<< ii << " , " << jj << " , " << 4 << std::endl;
1712           std::cout << D2EDX2(ii,jj, 4) << " D.F = " << MDiff(ii,jj) << std::endl;
1713         }
1714     }
1715   }
1716
1717 #endif
1718 */
1719   // ns1, ns2, np are copied locally to avois crushing the fields !
1720   ns1 = nsurf1;
1721   ns2 = nsurf2;
1722   np  = nplan;
1723   dnp = dnplan;
1724   d2np = d2nplan;
1725
1726   // Calculation of derivatives
1727
1728
1729   if ( ! pts1.IsEqual(pts2, 1.e-4)) {
1730     math_Gauss Resol(DEDX, 1.e-9); // Precise tolerance !!!!! 
1731     // Calculation of derivatives Processing Normal
1732     if (Resol.IsDone()) {  
1733       Resol.Solve(-DEDT, sol);
1734       D2EDX2.Multiply(sol, D2DXdSdt);    
1735       secmember = - (D2EDT2 + (2*D2EDXDT+D2DXdSdt)*sol);
1736       Resol.Solve(secmember);
1737       istgt = Standard_False;
1738     }
1739   }
1740   
1741   if (istgt) {
1742     math_SVD SingRS (DEDX);
1743     math_Vector Vbis(1,4);
1744     if (SingRS.IsDone()) {
1745       SingRS.Solve(-DEDT, sol, 1.e-6);
1746       D2EDX2.Multiply(sol, D2DXdSdt);
1747       Vbis = - (D2EDT2 + (2*D2EDXDT+D2DXdSdt)*sol);      
1748       SingRS.Solve( Vbis, secmember, 1.e-6);
1749       istgt = Standard_False;
1750     }
1751   }
1752
1753   if (!istgt) { 
1754     tg1.SetLinearForm(sol(1),d1u1, sol(2),d1v1);
1755     tg2.SetLinearForm(sol(3),d1u2, sol(4),d1v2);
1756
1757     dnorm1w.SetLinearForm(sol(1),dndu1, sol(2),dndv1, dn1w);
1758     dnorm2w.SetLinearForm(sol(3),dndu2, sol(4),dndv2, dn2w);
1759     temp.SetLinearForm(sol(1)*sol(1), d2u1, 
1760                        2*sol(1)*sol(2), d2uv1,
1761                        sol(2)*sol(2), d2v1);
1762
1763     dtg1.SetLinearForm(secmember(1),d1u1, secmember(2),d1v1, temp);
1764
1765     temp.SetLinearForm(sol(3)*sol(3), d2u2, 
1766                        2*sol(3)*sol(4), d2uv2,
1767                        sol(4)*sol(4), d2v2);
1768     dtg2.SetLinearForm(secmember(3),d1u2, secmember(4),d1v2, temp);    
1769
1770     temp.SetLinearForm(sol(1)*sol(1), d2ndu1, 
1771                        2*sol(1)*sol(2), d2nduv1,
1772                        sol(2)*sol(2), d2ndv1);
1773
1774     tempbis.SetLinearForm(2*sol(1), d2ndtu1, 
1775                           2*sol(2), d2ndtv1,
1776                           d2n1w);
1777     temp+= tempbis;
1778     d2norm1w.SetLinearForm(secmember(1),dndu1, secmember(2),dndv1, temp);
1779
1780
1781     temp.SetLinearForm(sol(3)*sol(3),   d2ndu2, 
1782                        2*sol(3)*sol(4), d2nduv2,
1783                        sol(4)*sol(4),   d2ndv2);
1784     tempbis.SetLinearForm(2*sol(3), d2ndtu2, 
1785                           2*sol(4), d2ndtv2,
1786                           d2n2w);
1787     temp+= tempbis;
1788     d2norm2w.SetLinearForm(secmember(3),dndu2, secmember(4),dndv2, temp);
1789   }
1790
1791   // Tops 2d  
1792   Poles2d(Poles2d.Lower()).SetCoord(X(1),X(2));
1793   Poles2d(Poles2d.Upper()).SetCoord(X(3),X(4));
1794   if (!istgt) {
1795     DPoles2d(Poles2d.Lower()) .SetCoord(sol(1),sol(2));
1796     DPoles2d(Poles2d.Upper()) .SetCoord(sol(3),sol(4));
1797     D2Poles2d(Poles2d.Lower()).SetCoord(secmember(1), secmember(2));
1798     D2Poles2d(Poles2d.Upper()).SetCoord(secmember(3), secmember(4));    
1799   }
1800
1801   // linear case is processed...
1802   if (mySShape == BlendFunc_Linear) {
1803     Poles(low) = pts1;
1804     Poles(upp) = pts2;
1805     Weights(low) = 1.0;
1806     Weights(upp) = 1.0;
1807     if (!istgt) {
1808       DPoles(low) = tg1;
1809       DPoles(upp) = tg2;
1810       DPoles(low) = dtg1;
1811       DPoles(upp) = dtg2;
1812       DWeights(low) = 0.0;
1813       DWeights(upp) = 0.0;
1814       D2Weights(low) = 0.0;
1815       D2Weights(upp) = 0.0;
1816     }
1817     return (!istgt);
1818   }
1819
1820   // Case of circle
1821   norm1 = nplan.Crossed(ns1).Magnitude();
1822   norm2 = nplan.Crossed(ns2).Magnitude();
1823   if (norm1 < Eps) {
1824     norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1825 #ifdef OCCT_DEBUG
1826     std::cout << " ConstRad : Surface singuliere " << std::endl;
1827 #endif
1828   }
1829   if (norm2 < Eps) {
1830     norm2 = 1; // Unsatisfactory, but it is not necessary to stop
1831 #ifdef OCCT_DEBUG
1832     std::cout << " ConstRad : Surface singuliere " << std::endl;
1833 #endif
1834  }
1835
1836   ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1, ns1);
1837   ns2.SetLinearForm(nplan.Dot(ns2)/norm2,nplan, -1./norm2, ns2);
1838
1839   Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1840   if (!istgt) {
1841     tgc.SetLinearForm(ray1,dnorm1w,tg1); 
1842     dtgc.SetLinearForm(ray1, d2norm1w, dtg1);
1843   }
1844
1845   // ns1 is oriented from the center to pts1 and ns2 from the center to pts2
1846   // trihedron ns1,ns2,nplan is made direct
1847
1848   if (ray1 > 0.) {
1849     ns1.Reverse();
1850     if (!istgt) {
1851       dnorm1w.Reverse();
1852       d2norm1w.Reverse();
1853     }
1854   }
1855   if (ray2 >0.) {
1856     ns2.Reverse();
1857     if (!istgt) {
1858       dnorm2w.Reverse();
1859       d2norm2w.Reverse();
1860     }
1861   }
1862   if (choix%2 != 0) {
1863     np.Reverse();
1864     dnp.Reverse();
1865     d2np.Reverse();
1866   }
1867   
1868   if (!istgt) {
1869     return GeomFill::GetCircle(myTConv, 
1870                                ns1, ns2, 
1871                                dnorm1w, dnorm2w,
1872                                d2norm1w, d2norm2w, 
1873                                np, dnp, d2np,
1874                                pts1, pts2, 
1875                                tg1, tg2, 
1876                                dtg1, dtg2,
1877                                Abs(ray1), 0, 0, 
1878                                Center, tgc, dtgc, 
1879                                Poles, 
1880                                DPoles,
1881                                D2Poles,
1882                                Weights, 
1883                                DWeights,
1884                                D2Weights); 
1885   }
1886   else {
1887     GeomFill::GetCircle(myTConv,
1888                         ns1, ns2, 
1889                         nplan, pts1, pts2,
1890                         Abs(ray1), Center, 
1891                         Poles, Weights);
1892     return Standard_False;
1893   }
1894 }
1895
1896 //=======================================================================
1897 //function : AxeRot
1898 //purpose  : 
1899 //=======================================================================
1900
1901 gp_Ax1 BlendFunc_ConstRad::AxeRot (const Standard_Real Prm)
1902 {
1903   gp_Ax1 axrot;
1904   gp_Vec dirax, d1gui, d2gui, np, dnp;
1905   gp_Pnt oriax, ptgui;
1906
1907   curv->D2(Prm,ptgui,d1gui,d2gui);
1908
1909   Standard_Real normtg = d1gui.Magnitude();
1910   np  = d1gui.Normalized();
1911   dnp.SetLinearForm(1./normtg, d2gui,
1912                    -1./normtg*(np.Dot(d2gui)), np);
1913
1914   dirax = np.Crossed(dnp);
1915   if (dirax.Magnitude() >= gp::Resolution()) {
1916
1917     axrot.SetDirection(dirax);
1918   }
1919   else {
1920     axrot.SetDirection(np);  // To avoid stop
1921   }
1922   if (dnp.Magnitude() >= gp::Resolution()) {
1923     oriax.SetXYZ(ptgui.XYZ()+
1924                  (normtg/dnp.Magnitude())*dnp.Normalized().XYZ());
1925   }
1926   else {
1927     oriax.SetXYZ(ptgui.XYZ());
1928   }
1929   axrot.SetLocation(oriax);
1930   return axrot;
1931 }
1932
1933 void BlendFunc_ConstRad::Resolution(const Standard_Integer IC2d,
1934                                     const Standard_Real Tol,
1935                                     Standard_Real& TolU,
1936                                     Standard_Real& TolV) const
1937 {
1938   if(IC2d == 1){
1939     TolU = surf1->UResolution(Tol);
1940     TolV = surf1->VResolution(Tol);
1941   }
1942   else {
1943     TolU = surf2->UResolution(Tol);
1944     TolV = surf2->VResolution(Tol);
1945   }
1946 }