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