src/BlendFunc/BlendFunc_ConstRad.cxx

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