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