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