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