src/BlendFunc/BlendFunc_ConstRad.cxx

   1 // Created on: 1993-12-02
   2 // Created by: Jacques GOUSSARD
   3 // Copyright (c) 1993-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 // Modified 09/09/1996 PMN Adde Nb(Intervalls), IsRationnal
  18 //                         Optimisation, use of GetCircle
  19 // Modified 20/02/1998 PMN Singular surfaces management
  20
  21 #include <BlendFunc_ConstRad.ixx>
  22
  23 #include <math_Gauss.hxx>
  24 #include <math_SVD.hxx>
  25 #include <gp.hxx>
  26 #include <BlendFunc.hxx>
  27 #include <GeomFill.hxx>
  28 #include <Standard_TypeDef.hxx>
  29 #include <Standard_DomainError.hxx>
  30 #include <Standard_NotImplemented.hxx>
  31 #include <ElCLib.hxx>
  32 #include <Precision.hxx>
  33
  34 #define Eps 1.e-15
  35
  36
  37 //=======================================================================
  38 //function : BlendFunc_ConstRad
  39 //purpose  :
  40 //=======================================================================
  41
  42 BlendFunc_ConstRad::BlendFunc_ConstRad(const Handle(Adaptor3d_HSurface)& S1,
  43                                        const Handle(Adaptor3d_HSurface)& S2,
  44                                        const Handle(Adaptor3d_HCurve)& C)
  45                                  :
  46                                   surf1(S1),surf2(S2),
  47                                   curv(C), tcurv(C),
  48                                   istangent(Standard_True),
  49                                   xval(1,4),
  50                                   E(1,4), DEDX(1,4,1,4),  DEDT(1,4),
  51                                   D2EDX2(4,4,4),
  52                                   D2EDXDT(1,4,1,4), D2EDT2(1,4),
  53                                   maxang(RealFirst()), minang(RealLast()),
  54                                   distmin(RealLast()),
  55                                   mySShape(BlendFunc_Rational)
  56 {
  57 // Initialisaton of cash control variables.
  58   tval = -9.876e100;
  59   xval.Init(-9.876e100);
  60   myXOrder = -1;
  61   myTOrder = -1;
  62 }
  63
  64 //=======================================================================
  65 //function : NbEquations
  66 //purpose  :
  67 //=======================================================================
  68
  69 Standard_Integer BlendFunc_ConstRad::NbEquations () const
  70 {
  71   return 4;
  72 }
  73
  74
  75 //=======================================================================
  76 //function : Set
  77 //purpose  :
  78 //=======================================================================
  79
  80 void BlendFunc_ConstRad::Set(const Standard_Real Radius, const Standard_Integer Choix)
  81 {
  82   choix = Choix;
  83   switch (choix) {
  84   case 1:
  85   case 2:
  86     {
  87       ray1 = -Radius;
  88       ray2 = -Radius;
  89     }
  90     break;
  91   case 3:
  92   case 4:
  93     {
  94       ray1 = Radius;
  95       ray2 = -Radius;
  96     }
  97     break;
  98   case 5:
  99   case 6:
 100     {
 101       ray1 = Radius;
 102       ray2 = Radius;
 103     }
 104     break;
 105   case 7:
 106   case 8:
 107     {
 108       ray1 = -Radius;
 109       ray2 = Radius;
 110     }
 111     break;
 112   default:
 113     ray1 = ray2 = -Radius;
 114   }
 115 }
 116
 117 //=======================================================================
 118 //function : Set
 119 //purpose  :
 120 //=======================================================================
 121
 122 void BlendFunc_ConstRad::Set(const BlendFunc_SectionShape TypeSection)
 123 {
 124   mySShape = TypeSection;
 125 }
 126
 127
 128 //=======================================================================
 129 //function : ComputeValues
 130 //purpose  : OBLIGATORY passage for all calculations
 131 //           This method manages positioning on Surfaces and Curves
 132 //           Calculate the equations and their partial derivates
 133 //           Stock certain intermediate results in fields to
 134 //           use in other methods.
 135 //=======================================================================
 136
 137 Standard_Boolean BlendFunc_ConstRad::ComputeValues(const math_Vector& X,
 138                                                    const Standard_Integer Order,
 139                                                    const Standard_Boolean byParam,
 140                                                    const Standard_Real Param)
 141 {
 142  // static declaration to avoid systematic reallocation
 143
 144  static gp_Vec d3u1,d3v1,d3uuv1,d3uvv1,d3u2,d3v2,d3uuv2,d3uvv2;
 145  static gp_Vec d1gui, d2gui, d3gui;
 146  static gp_Pnt ptgui;
 147  static Standard_Real invnormtg, dinvnormtg;
 148  Standard_Real T =  Param, aux;
 149
 150  // Case of implicite parameter
 151  if ( !byParam) { T = param;}
 152
 153  // Is the work already done ?
 154  Standard_Boolean myX_OK = (Order<=myXOrder) ;
 155  for (Standard_Integer ii=1; ((ii<=X.Length()) && myX_OK); ii++) {
 156    myX_OK = ( X(ii) == xval(ii) );
 157  }
 158
 159  Standard_Boolean t_OK =( (T == tval)
 160                        && ((Order<=myTOrder)||(!byParam)) );
 161
 162  if (myX_OK && (t_OK) ) {
 163    return Standard_True;
 164  }
 165
 166  // Processing of t
 167  if (!t_OK) {
 168    tval = T;
 169    if (byParam) { myTOrder = Order;}
 170    else         { myTOrder = 0;}
 171    //----- Positioning on the curve ----------------
 172    switch (myTOrder) {
 173    case 0 :
 174      {
 175        tcurv->D1(T, ptgui, d1gui);
 176        nplan  = d1gui.Normalized();
 177      }
 178      break;
 179
 180    case 1 :
 181      {
 182        tcurv->D2(T,ptgui,d1gui,d2gui);
 183        nplan  = d1gui.Normalized();
 184        invnormtg = ((Standard_Real) 1 ) / d1gui.Magnitude();
 185        dnplan.SetLinearForm(invnormtg, d2gui,
 186                             -invnormtg*(nplan.Dot(d2gui)), nplan);
 187      break;
 188      }
 189    case 2 :
 190      {
 191        tcurv->D3(T,ptgui,d1gui,d2gui,d3gui);
 192        nplan  = d1gui.Normalized();
 193        invnormtg = ((Standard_Real) 1 ) / d1gui.Magnitude();
 194        dnplan.SetLinearForm(invnormtg, d2gui,
 195                            -invnormtg*(nplan.Dot(d2gui)), nplan);
 196        dinvnormtg = - nplan.Dot(d2gui)*invnormtg*invnormtg;
 197        d2nplan.SetLinearForm(invnormtg,  d3gui,
 198                              dinvnormtg, d2gui);
 199        aux = dinvnormtg*(nplan.Dot(d2gui)) + invnormtg*( dnplan.Dot(d2gui)
 200                                                        + nplan.Dot(d3gui) );
 201        d2nplan.SetLinearForm(-invnormtg*(nplan.Dot(d2gui)), dnplan,
 202                              -aux, nplan, d2nplan );
 203        break;
 204    }
 205    default:
 206      return Standard_False;
 207    }
 208  }
 209
 210  // Processing of X
 211  if (!myX_OK) {
 212    xval = X;
 213    myXOrder = Order;
 214    //-------------- Positioning on surfaces -----------------
 215    switch (myXOrder) {
 216    case 0 :
 217      {
 218        surf1->D1(X(1),X(2),pts1,d1u1,d1v1);
 219        nsurf1 = d1u1.Crossed(d1v1);
 220        surf2->D1(X(3),X(4),pts2,d1u2,d1v2);
 221        nsurf2 = d1u2.Crossed(d1v2);
 222        break;
 223      }
 224    case 1 :
 225      {
 226        surf1->D2(X(1),X(2),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1);
 227        nsurf1 = d1u1.Crossed(d1v1);
 228        dns1u1 = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
 229        dns1v1 = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
 230        surf2->D2(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2);
 231        nsurf2 = d1u2.Crossed(d1v2);
 232        dns1u2 = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
 233        dns1v2 = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
 234        break;
 235      }
 236    case 2 :
 237      {
 238        surf1->D3(X(1),X(2),pts1,d1u1,d1v1,d2u1,d2v1,d2uv1,d3u1,d3v1,d3uuv1,d3uvv1);
 239        nsurf1 = d1u1.Crossed(d1v1);
 240        dns1u1 = d2u1.Crossed(d1v1).Added(d1u1.Crossed(d2uv1));
 241        dns1v1 = d2uv1.Crossed(d1v1).Added(d1u1.Crossed(d2v1));
 242
 243        surf2->D3(X(3),X(4),pts2,d1u2,d1v2,d2u2,d2v2,d2uv2,d3u2,d3v2,d3uuv2,d3uvv2);
 244        nsurf2 = d1u2.Crossed(d1v2);
 245        dns1u2 = d2u2.Crossed(d1v2).Added(d1u2.Crossed(d2uv2));
 246        dns1v2 = d2uv2.Crossed(d1v2).Added(d1u2.Crossed(d2v2));
 247        break;
 248      }
 249    default:
 250      return Standard_False;
 251    }
 252    // Case of degenerated surfaces
 253    if (nsurf1.Magnitude() < Eps ) {
 254      //gp_Vec normal;
 255      gp_Pnt2d P(X(1), X(2));
 256      if (Order == 0) BlendFunc::ComputeNormal(surf1, P, nsurf1);
 257      else  BlendFunc::ComputeDNormal(surf1, P, nsurf1, dns1u1, dns1v1);
 258    }
 259   if (nsurf2.Magnitude() < Eps) {
 260      //gp_Vec normal;
 261      gp_Pnt2d P(X(3), X(4));
 262      if (Order==0) BlendFunc::ComputeNormal(surf2, P, nsurf2);
 263      else BlendFunc::ComputeDNormal(surf2, P, nsurf2, dns1u2, dns1v2);
 264    }
 265  }
 266
 267  // -------------------- Positioning of order 0 ---------------------
 268  Standard_Real invnorm1, invnorm2, ndotns1, ndotns2, theD;
 269  gp_Vec ncrossns1,ncrossns2,resul,temp;
 270
 271  theD =  - (nplan.XYZ().Dot(ptgui.XYZ()));
 272
 273  E(1) = (nplan.X() * (pts1.X() + pts2.X()) +
 274          nplan.Y() * (pts1.Y() + pts2.Y()) +
 275          nplan.Z() * (pts1.Z() + pts2.Z())) /2 + theD;
 276
 277  ncrossns1 = nplan.Crossed(nsurf1);
 278  ncrossns2 = nplan.Crossed(nsurf2);
 279
 280  invnorm1 = ncrossns1.Magnitude();
 281  invnorm2 = ncrossns2.Magnitude();
 282
 283  if (invnorm1 > Eps) invnorm1 = ((Standard_Real) 1) /invnorm1;
 284  else {
 285    invnorm1 = 1; // Unsatisfactory, but it is not necessary to crash
 286 #ifdef OCCT_DEBUG
 287    cout << " ConstRad : Surface singuliere " << endl;
 288 #endif
 289  }
 290  if (invnorm2 > Eps) invnorm2 = ((Standard_Real) 1) /invnorm2;
 291  else {
 292    invnorm2 = 1; //  Unsatisfactory, but it is not necessary to crash
 293 #ifdef OCCT_DEBUG
 294    cout << " ConstRad : Surface singuliere " << endl;
 295 #endif
 296  }
 297
 298  ndotns1 = nplan.Dot(nsurf1);
 299  ndotns2 = nplan.Dot(nsurf2);
 300
 301  temp.SetLinearForm(ndotns1,nplan,-1.,nsurf1);
 302  temp.Multiply(invnorm1);
 303
 304  resul.SetLinearForm(ray1,temp,gp_Vec(pts2,pts1));
 305  temp.SetLinearForm(ndotns2,nplan,-1.,nsurf2);
 306  temp.Multiply(invnorm2);
 307  resul.Subtract(ray2*temp);
 308
 309  E(2) = resul.X();
 310  E(3) = resul.Y();
 311  E(4) = resul.Z();
 312
 313  // -------------------- Positioning of order 1 ---------------------
 314  if (Order >= 1) {
 315    Standard_Real  grosterme, cube, carre;
 316
 317
 318    DEDX(1,1) = nplan.Dot(d1u1)/2;
 319    DEDX(1,2) = nplan.Dot(d1v1)/2;
 320    DEDX(1,3) = nplan.Dot(d1u2)/2;
 321    DEDX(1,4) = nplan.Dot(d1v2)/2;
 322
 323    cube =invnorm1*invnorm1*invnorm1;
 324    // Derived in relation to u1
 325    grosterme = - ncrossns1.Dot(nplan.Crossed(dns1u1))*cube;
 326    dndu1.SetLinearForm( grosterme*ndotns1
 327                       + invnorm1*nplan.Dot(dns1u1), nplan,
 328                       - grosterme, nsurf1,
 329                       - invnorm1,  dns1u1 );
 330
 331    resul.SetLinearForm(ray1, dndu1, d1u1);
 332    DEDX(2,1) = resul.X();
 333    DEDX(3,1) = resul.Y();
 334    DEDX(4,1) = resul.Z();
 335
 336    // Derived in relation to v1
 337
 338    grosterme = - ncrossns1.Dot(nplan.Crossed(dns1v1))*cube;
 339    dndv1.SetLinearForm( grosterme*ndotns1
 340                        +invnorm1*nplan.Dot(dns1v1), nplan,
 341                        -grosterme, nsurf1,
 342                        -invnorm1, dns1v1);
 343
 344    resul.SetLinearForm(ray1, dndv1, d1v1);
 345    DEDX(2,2) = resul.X();
 346    DEDX(3,2) = resul.Y();
 347    DEDX(4,2) = resul.Z();
 348
 349    cube = invnorm2*invnorm2*invnorm2;
 350    // Derived in relation to u2
 351    grosterme = - ncrossns2.Dot(nplan.Crossed(dns1u2))*cube;
 352    dndu2.SetLinearForm( grosterme*ndotns2
 353                        +invnorm2*nplan.Dot(dns1u2), nplan,
 354                        -grosterme, nsurf2,
 355                        -invnorm2, dns1u2);
 356
 357    resul.SetLinearForm(-ray2, dndu2, -1, d1u2);
 358    DEDX(2,3) = resul.X();
 359    DEDX(3,3) = resul.Y();
 360    DEDX(4,3) = resul.Z();
 361
 362    // Derived in relation to v2
 363    grosterme = -ncrossns2.Dot(nplan.Crossed(dns1v2))*cube;
 364    dndv2.SetLinearForm( grosterme*ndotns2
 365                        +invnorm2*nplan.Dot(dns1v2), nplan,
 366                        -grosterme, nsurf2,
 367                        -invnorm2 , dns1v2);
 368
 369    resul.SetLinearForm(-ray2,dndv2, -1, d1v2);
 370    DEDX(2,4) = resul.X();
 371    DEDX(3,4) = resul.Y();
 372    DEDX(4,4) = resul.Z();
 373
 374    if (byParam) {
 375      temp.SetXYZ( (pts1.XYZ()+pts2.XYZ())/2 - ptgui.XYZ());
 376      // Derived from n1 in relation to w
 377      grosterme = ncrossns1.Dot(dnplan.Crossed(nsurf1))*invnorm1*invnorm1;
 378      dn1w.SetLinearForm((dnplan.Dot(nsurf1)-grosterme*ndotns1)*invnorm1, nplan,
 379                          ndotns1*invnorm1,dnplan,
 380                          grosterme*invnorm1,nsurf1);
 381
 382      // Derivee from n2 in relation to w
 383      grosterme = ncrossns2.Dot(dnplan.Crossed(nsurf2))*invnorm2*invnorm2;
 384      dn2w.SetLinearForm((dnplan.Dot(nsurf2)-grosterme*ndotns2)*invnorm2,nplan,
 385                         ndotns2*invnorm2,dnplan,
 386                         grosterme*invnorm2,nsurf2);
 387
 388
 389      DEDT(1) =  dnplan.Dot(temp) - 1./invnormtg ;
 390      DEDT(2) =  ray1*dn1w.X() - ray2*dn2w.X();
 391      DEDT(3) =  ray1*dn1w.Y() - ray2*dn2w.Y();
 392      DEDT(4) =  ray1*dn1w.Z() - ray2*dn2w.Z();
 393    }
 394    // ------   Positioning of order 2  -----------------------------
 395    if (Order == 2) {
 396 //     gp_Vec d2ndu1,  d2ndu2, d2ndv1, d2ndv2, d2nduv1, d2nduv2;
 397      gp_Vec d2ns1u1,  d2ns1u2, d2ns1v1, d2ns1v2, d2ns1uv1, d2ns1uv2;
 398      Standard_Real uterm, vterm, smallterm, p1, p2, p12;
 399      Standard_Real DPrim, DSecn;
 400      D2EDX2.Init(0);
 401
 402      D2EDX2(1, 1, 1) = nplan.Dot(d2u1)/2;
 403      D2EDX2(1, 2, 1) = D2EDX2(1, 1, 2) = nplan.Dot(d2uv1)/2;
 404      D2EDX2(1, 2, 2) = nplan.Dot(d2v1)/2;
 405
 406      D2EDX2(1, 3, 3) = nplan.Dot(d2u2)/2;
 407      D2EDX2(1, 4, 3) = D2EDX2(1, 3, 4) = nplan.Dot(d2uv2)/2;
 408      D2EDX2(1, 4, 4) = nplan.Dot(d2v2)/2;
 409                                // ================
 410                                // ==  Surface 1 ==
 411                                // ================
 412      carre = invnorm1*invnorm1;
 413      cube  = carre*invnorm1;
 414      // Derived double compared to u1
 415        // Derived from the norm
 416      d2ns1u1.SetLinearForm(1, d3u1.Crossed(d1v1),
 417                            2, d2u1.Crossed(d2uv1),
 418                            1, d1u1.Crossed(d3uuv1));
 419      DPrim = ncrossns1.Dot(nplan.Crossed(dns1u1));
 420      smallterm =  - 2*DPrim*cube;
 421      DSecn = ncrossns1.Dot(nplan.Crossed(d2ns1u1))
 422           + (nplan.Crossed(dns1u1)).SquareMagnitude();
 423      grosterme  = (3*DPrim*DPrim*carre - DSecn) * cube;
 424
 425      temp.SetLinearForm( grosterme*ndotns1, nplan,
 426                         -grosterme        , nsurf1);
 427      p1 =  nplan.Dot(dns1u1);
 428      p2 =  nplan.Dot(d2ns1u1);
 429      d2ndu1.SetLinearForm( invnorm1*p2
 430                           +smallterm*p1, nplan,
 431                           -smallterm,    dns1u1,
 432                           -invnorm1,     d2ns1u1);
 433      d2ndu1 += temp;
 434      resul.SetLinearForm(ray1, d2ndu1, d2u1);
 435      D2EDX2(2,1,1) = resul.X();
 436      D2EDX2(3,1,1) = resul.Y();
 437      D2EDX2(4,1,1) = resul.Z();
 438
 439      // Derived double compared to u1, v1
 440        // Derived from the norm
 441      d2ns1uv1 =  (d3uuv1.Crossed(d1v1))
 442               +  (d2u1  .Crossed(d2v1))
 443               +  (d1u1  .Crossed(d3uvv1));
 444      uterm =  ncrossns1.Dot(nplan.Crossed(dns1u1));
 445      vterm =  ncrossns1.Dot(nplan.Crossed(dns1v1));
 446      DSecn = (nplan.Crossed(dns1v1)).Dot(nplan.Crossed(dns1u1))
 447            +  ncrossns1.Dot(nplan.Crossed(d2ns1uv1));
 448      grosterme  = (3*uterm*vterm*carre-DSecn)*cube;
 449      uterm *= -cube; //and only now
 450      vterm *= -cube;
 451
 452      p1 = nplan.Dot(dns1u1);
 453      p2 = nplan.Dot(dns1v1);
 454      temp.SetLinearForm(  grosterme*ndotns1, nplan,
 455                         - grosterme, nsurf1,
 456                         - invnorm1, d2ns1uv1);
 457      d2nduv1.SetLinearForm( invnorm1*nplan.Dot(d2ns1uv1)
 458                           + uterm*p2
 459                           + vterm*p1, nplan,
 460                           - uterm,    dns1v1,
 461                           - vterm,    dns1u1);
 462
 463      d2nduv1 += temp;
 464      resul.SetLinearForm(ray1, d2nduv1, d2uv1);
 465
 466      D2EDX2(2,2,1) = D2EDX2(2,1,2) = resul.X();
 467      D2EDX2(3,2,1) = D2EDX2(3,1,2) = resul.Y();
 468      D2EDX2(4,2,1) = D2EDX2(4,1,2) = resul.Z();
 469
 470      // Derived double compared to v1
 471        // Derived from the norm
 472      d2ns1v1.SetLinearForm(1, d1u1.Crossed(d3v1),
 473                            2, d2uv1.Crossed(d2v1),
 474                            1, d3uvv1.Crossed(d1v1));
 475      DPrim = ncrossns1.Dot(nplan.Crossed(dns1v1));
 476      smallterm =  - 2*DPrim*cube;
 477      DSecn = ncrossns1.Dot(nplan.Crossed(d2ns1v1))
 478           + (nplan.Crossed(dns1v1)).SquareMagnitude();
 479      grosterme  = (3*DPrim*DPrim*carre - DSecn) * cube;
 480
 481      p1 =  nplan.Dot(dns1v1);
 482      p2 =  nplan.Dot(d2ns1v1);
 483      temp.SetLinearForm( grosterme*ndotns1, nplan,
 484                         -grosterme        , nsurf1);
 485      d2ndv1.SetLinearForm( invnorm1*p2
 486                           +smallterm*p1, nplan,
 487                           -smallterm,    dns1v1,
 488                           -invnorm1,     d2ns1v1);
 489      d2ndv1 += temp;
 490      resul.SetLinearForm(ray1, d2ndv1, d2v1);
 491
 492      D2EDX2(2,2,2) = resul.X();
 493      D2EDX2(3,2,2) = resul.Y();
 494      D2EDX2(4,2,2) = resul.Z();
 495                                // ================
 496                                // ==  Surface 2 ==
 497                                // ================
 498      carre = invnorm2*invnorm2;
 499      cube  = carre*invnorm2;
 500      // Derived double compared to u2
 501        // Derived from the norm
 502      d2ns1u2.SetLinearForm(1, d3u2.Crossed(d1v2),
 503                            2, d2u2.Crossed(d2uv2),
 504                            1, d1u2.Crossed(d3uuv2));
 505      DPrim = ncrossns2.Dot(nplan.Crossed(dns1u2));
 506      smallterm =  - 2*DPrim*cube;
 507      DSecn = ncrossns2.Dot(nplan.Crossed(d2ns1u2))
 508           + (nplan.Crossed(dns1u2)).SquareMagnitude();
 509      grosterme  = (3*DPrim*DPrim*carre - DSecn) * cube;
 510
 511      temp.SetLinearForm( grosterme*ndotns2, nplan,
 512                         -grosterme        , nsurf2);
 513      p1 =  nplan.Dot(dns1u2);
 514      p2 =  nplan.Dot(d2ns1u2);
 515      d2ndu2.SetLinearForm( invnorm2*p2
 516                           +smallterm*p1, nplan,
 517                           -smallterm,    dns1u2,
 518                           -invnorm2,     d2ns1u2);
 519      d2ndu2 += temp;
 520      resul.SetLinearForm(-ray2, d2ndu2, -1, d2u2);
 521      D2EDX2(2,3,3) = resul.X();
 522      D2EDX2(3,3,3) = resul.Y();
 523      D2EDX2(4,3,3) = resul.Z();
 524
 525      // Derived double compared to u2, v2
 526        // Derived from the norm
 527      d2ns1uv2 =  (d3uuv2.Crossed(d1v2))
 528               +  (d2u2  .Crossed(d2v2))
 529               +  (d1u2  .Crossed(d3uvv2));
 530      uterm =  ncrossns2.Dot(nplan.Crossed(dns1u2));
 531      vterm =  ncrossns2.Dot(nplan.Crossed(dns1v2));
 532      DSecn = (nplan.Crossed(dns1v2)).Dot(nplan.Crossed(dns1u2))
 533            +  ncrossns2.Dot(nplan.Crossed(d2ns1uv2));
 534      grosterme  = (3*uterm*vterm*carre-DSecn)*cube;
 535      uterm *= -cube; //and only now
 536      vterm *= -cube;
 537
 538      p1 = nplan.Dot(dns1u2);
 539      p2 = nplan.Dot(dns1v2);
 540      temp.SetLinearForm(  grosterme*ndotns2, nplan,
 541                         - grosterme, nsurf2,
 542                         - invnorm2, d2ns1uv2);
 543      d2nduv2.SetLinearForm( invnorm2*nplan.Dot(d2ns1uv2)
 544                           + uterm*p2
 545                           + vterm*p1, nplan,
 546                           - uterm,    dns1v2,
 547                           - vterm,    dns1u2);
 548
 549      d2nduv2 += temp;
 550      resul.SetLinearForm(-ray2, d2nduv2, -1, d2uv2);
 551
 552      D2EDX2(2,4,3) = D2EDX2(2,3,4) = resul.X();
 553      D2EDX2(3,4,3) = D2EDX2(3,3,4) = resul.Y();
 554      D2EDX2(4,4,3) = D2EDX2(4,3,4) = resul.Z();
 555
 556      // Derived double compared to v2
 557        // Derived from the norm
 558      d2ns1v2.SetLinearForm(1, d1u2.Crossed(d3v2),
 559                            2, d2uv2.Crossed(d2v2),
 560                            1, d3uvv2.Crossed(d1v2));
 561      DPrim = ncrossns2.Dot(nplan.Crossed(dns1v2));
 562      smallterm =  - 2*DPrim*cube;
 563      DSecn = ncrossns2.Dot(nplan.Crossed(d2ns1v2))
 564           + (nplan.Crossed(dns1v2)).SquareMagnitude();
 565      grosterme  = (3*DPrim*DPrim*carre - DSecn) * cube;
 566
 567      p1 =  nplan.Dot(dns1v2);
 568      p2 =  nplan.Dot(d2ns1v2);
 569      temp.SetLinearForm( grosterme*ndotns2, nplan,
 570                         -grosterme        , nsurf2);
 571      d2ndv2.SetLinearForm( invnorm2*p2
 572                           +smallterm*p1, nplan,
 573                           -smallterm,    dns1v2,
 574                           -invnorm2,     d2ns1v2);
 575      d2ndv2 += temp;
 576      resul.SetLinearForm(-ray2, d2ndv2, -1, d2v2);
 577
 578      D2EDX2(2,4,4) = resul.X();
 579      D2EDX2(3,4,4) = resul.Y();
 580      D2EDX2(4,4,4) = resul.Z();
 581
 582      if (byParam) {
 583        Standard_Real tterm;
 584         //  ---------- Derivation double in t, X --------------------------
 585        D2EDXDT(1,1) = dnplan.Dot(d1u1)/2;
 586        D2EDXDT(1,2) = dnplan.Dot(d1v1)/2;
 587        D2EDXDT(1,3) = dnplan.Dot(d1u2)/2;
 588        D2EDXDT(1,4) = dnplan.Dot(d1v2)/2;
 589
 590        carre = invnorm1*invnorm1;
 591        cube = carre*invnorm1;
 592        //--> Derived compared to u1 and t
 593        tterm =  ncrossns1.Dot(dnplan.Crossed(nsurf1));
 594        smallterm  = - tterm*cube;
 595        // Derived from the norm
 596        uterm =  ncrossns1.Dot(nplan. Crossed(dns1u1));
 597        DSecn = (nplan.Crossed(dns1u1)).Dot(dnplan.Crossed(nsurf1))
 598               + ncrossns1.Dot(dnplan.Crossed(dns1u1));
 599        grosterme = (3*uterm*tterm*carre - DSecn) * cube;
 600        uterm *= -cube;
 601
 602        p1 = dnplan.Dot(nsurf1);
 603        p2 = nplan. Dot(dns1u1);
 604        p12 = dnplan.Dot(dns1u1);
 605
 606        d2ndtu1.SetLinearForm( invnorm1*p12
 607                             + smallterm*p2
 608                             + uterm*p1
 609                             + grosterme*ndotns1, nplan,
 610                               invnorm1*p2
 611                             + uterm*ndotns1,     dnplan,
 612                             - smallterm,         dns1u1);
 613        d2ndtu1 -= grosterme*nsurf1;
 614
 615        resul = ray1*d2ndtu1;
 616        D2EDXDT(2,1) = resul.X();
 617        D2EDXDT(3,1) = resul.Y();
 618        D2EDXDT(4,1) = resul.Z();
 619
 620        //--> Derived compared to v1 and t
 621        // Derived from the norm
 622        uterm =  ncrossns1.Dot(nplan. Crossed(dns1v1));
 623        DSecn = (nplan. Crossed(dns1v1)).Dot(dnplan.Crossed(nsurf1))
 624               + ncrossns1.Dot(dnplan.Crossed(dns1v1));
 625        grosterme = (3*uterm*tterm*carre - DSecn) * cube;
 626        uterm *= -cube;
 627
 628        p1 =  dnplan.Dot(nsurf1);
 629        p2 =  nplan. Dot(dns1v1);
 630        p12 = dnplan.Dot(dns1v1);
 631        d2ndtv1.SetLinearForm( invnorm1*p12
 632                            + uterm*p1
 633                            + smallterm*p2
 634                            + grosterme*ndotns1, nplan,
 635                              invnorm1*p2
 636                            + uterm*ndotns1,    dnplan,
 637                            - smallterm     ,   dns1v1);
 638        d2ndtv1 -= grosterme*nsurf1;
 639
 640        resul = ray1* d2ndtv1;
 641        D2EDXDT(2,2) = resul.X();
 642        D2EDXDT(3,2) = resul.Y();
 643        D2EDXDT(4,2) = resul.Z();
 644
 645        carre = invnorm2*invnorm2;
 646        cube = carre*invnorm2;
 647        //--> Derived compared to u2 and t
 648        tterm =  ncrossns2.Dot(dnplan.Crossed(nsurf2));
 649        smallterm = -tterm*cube;
 650        // Derived from the norm
 651        uterm =  ncrossns2.Dot(nplan. Crossed(dns1u2));
 652        DSecn = (nplan. Crossed(dns1u2)).Dot(dnplan.Crossed(nsurf2))
 653              + ncrossns2.Dot(dnplan.Crossed(dns1u2));
 654        grosterme = (3*uterm*tterm*carre - DSecn) * cube;
 655        uterm *= -cube;
 656
 657        p1 =  dnplan.Dot(nsurf2);
 658        p2 =  nplan. Dot(dns1u2);
 659        p12 = dnplan.Dot(dns1u2);
 660
 661        d2ndtu2.SetLinearForm( invnorm2*p12
 662                            + smallterm*p2
 663                            + uterm*p1
 664                            + grosterme*ndotns2, nplan,
 665                              invnorm2*p2
 666                            + uterm*ndotns2,    dnplan,
 667                             -smallterm     ,   dns1u2);
 668        d2ndtu2 -= grosterme*nsurf2;
 669
 670        resul = - ray2*d2ndtu2;
 671        D2EDXDT(2,3) = resul.X();
 672        D2EDXDT(3,3) = resul.Y();
 673        D2EDXDT(4,3) = resul.Z();
 674
 675        //--> Derived compared to v2 and t
 676        // Derived from the norm
 677        uterm =  ncrossns2.Dot(nplan. Crossed(dns1v2));
 678        DSecn = (nplan.Crossed(dns1v2)).Dot(dnplan.Crossed(nsurf2))
 679              +  ncrossns2.Dot(dnplan.Crossed(dns1v2));
 680        grosterme = (3*uterm*tterm*carre - DSecn) * cube;
 681        uterm *= - cube;
 682
 683        p1 =  dnplan.Dot(nsurf2);
 684        p2 =  nplan. Dot(dns1v2);
 685        p12 = dnplan.Dot(dns1v2);
 686
 687        d2ndtv2.SetLinearForm( invnorm2*p12
 688                             + smallterm*p2
 689                             + uterm*p1
 690                             + grosterme*ndotns2, nplan,
 691                               invnorm2*p2
 692                             + uterm*ndotns2,     dnplan,
 693                              -smallterm        , dns1v2);
 694        d2ndtv2 -= grosterme*nsurf2;
 695
 696        resul = - ray2 * d2ndtv2;
 697        D2EDXDT(2,4) = resul.X();
 698        D2EDXDT(3,4) = resul.Y();
 699        D2EDXDT(4,4) = resul.Z();
 700
 701
 702         //  ---------- Derivation double in t -----------------------------
 703     // Derived from n1 compared to w
 704        carre = invnorm1*invnorm1;
 705        cube = carre*invnorm1;
 706        // Derived from the norm
 707        DPrim =  ncrossns1.Dot(dnplan.Crossed(nsurf1));
 708        smallterm = - 2*DPrim*cube;
 709        DSecn     = (dnplan.Crossed(nsurf1)).SquareMagnitude()
 710                   + ncrossns1.Dot(d2nplan.Crossed(nsurf1));
 711        grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
 712
 713        p1 =  dnplan. Dot(nsurf1);
 714        p2 =  d2nplan.Dot(nsurf1);
 715
 716        temp.SetLinearForm( grosterme*ndotns1, nplan,
 717                           -grosterme , nsurf1);
 718        d2n1w.SetLinearForm(  smallterm*p1
 719                            + invnorm1*p2,      nplan,
 720                              smallterm*ndotns1
 721                            + 2*invnorm1*p1,    dnplan,
 722                              ndotns1*invnorm1,   d2nplan);
 723        d2n1w += temp;
 724
 725      // Derived from n2 compared to w
 726        carre = invnorm2*invnorm2;
 727        cube = carre*invnorm2;
 728        // Derived from the norm
 729        DPrim =  ncrossns2.Dot(dnplan.Crossed(nsurf2));
 730        smallterm = - 2*DPrim*cube;
 731        DSecn     = (dnplan.Crossed(nsurf2)).SquareMagnitude()
 732                   + ncrossns2.Dot(d2nplan.Crossed(nsurf2));
 733        grosterme = (3*DPrim*DPrim*carre - DSecn) * cube;
 734
 735        p1 =  dnplan. Dot(nsurf2);
 736        p2 =  d2nplan.Dot(nsurf2);
 737
 738        temp.SetLinearForm( grosterme*ndotns2, nplan,
 739                           -grosterme ,        nsurf2);
 740        d2n2w.SetLinearForm(  smallterm*p1
 741                            + invnorm2*p2,        nplan,
 742                              smallterm*ndotns2
 743                            + 2*invnorm2*p1,     dnplan,
 744                              ndotns2*invnorm2,  d2nplan);
 745        d2n2w += temp;
 746
 747        temp.SetXYZ( (pts1.XYZ()+pts2.XYZ())/2 - ptgui.XYZ());
 748        D2EDT2(1) = d2nplan.Dot(temp) - 2*dnplan.Dot(d1gui) - nplan.Dot(d2gui);
 749        D2EDT2(2) = ray1*d2n1w.X() - ray2*d2n2w.X();
 750        D2EDT2(3) = ray1*d2n1w.Y() - ray2*d2n2w.Y();
 751        D2EDT2(4) = ray1*d2n1w.Z() - ray2*d2n2w.Z();
 752      }
 753    }
 754  }
 755  return Standard_True;
 756 }
 757
 758 //=======================================================================
 759 //function : Set
 760 //purpose  :
 761 //=======================================================================
 762
 763 void BlendFunc_ConstRad::Set(const Standard_Real Param)
 764 {
 765   param = Param;
 766 }
 767
 768 //=======================================================================
 769 //function : Set
 770 //purpose  : Segmentation of the useful part of the curve
 771 //           Precision is taken at random and small !?
 772 //=======================================================================
 773
 774 void BlendFunc_ConstRad::Set(const Standard_Real First, const Standard_Real Last)
 775 {
 776   tcurv = curv->Trim(First, Last, 1.e-12);
 777 }
 778
 779 //=======================================================================
 780 //function : GetTolerance
 781 //purpose  :
 782 //=======================================================================
 783
 784 void BlendFunc_ConstRad::GetTolerance(math_Vector& Tolerance, const Standard_Real Tol) const
 785 {
 786   Tolerance(1) = surf1->UResolution(Tol);
 787   Tolerance(2) = surf1->VResolution(Tol);
 788   Tolerance(3) = surf2->UResolution(Tol);
 789   Tolerance(4) = surf2->VResolution(Tol);
 790 }
 791
 792
 793 //=======================================================================
 794 //function : GetBounds
 795 //purpose  :
 796 //=======================================================================
 797
 798 void BlendFunc_ConstRad::GetBounds(math_Vector& InfBound, math_Vector& SupBound) const
 799 {
 800   InfBound(1) = surf1->FirstUParameter();
 801   InfBound(2) = surf1->FirstVParameter();
 802   InfBound(3) = surf2->FirstUParameter();
 803   InfBound(4) = surf2->FirstVParameter();
 804   SupBound(1) = surf1->LastUParameter();
 805   SupBound(2) = surf1->LastVParameter();
 806   SupBound(3) = surf2->LastUParameter();
 807   SupBound(4) = surf2->LastVParameter();
 808
 809   for(Standard_Integer i = 1; i <= 4; i++){
 810     if(!Precision::IsInfinite(InfBound(i)) &&
 811        !Precision::IsInfinite(SupBound(i))) {
 812       Standard_Real range = (SupBound(i) - InfBound(i));
 813       InfBound(i) -= range;
 814       SupBound(i) += range;
 815     }
 816   }
 817 }
 818
 819
 820 //=======================================================================
 821 //function : IsSolution
 822 //purpose  :
 823 //=======================================================================
 824
 825 Standard_Boolean BlendFunc_ConstRad::IsSolution(const math_Vector& Sol, const Standard_Real Tol)
 826 {
 827   Standard_Real norm, Cosa, Sina, Angle;
 828   Standard_Boolean Ok=Standard_True;
 829
 830   Ok = ComputeValues(Sol, 1, Standard_True, param);
 831
 832   if (Abs(E(1)) <= Tol &&
 833       E(2)*E(2) + E(3)*E(3) +  E(4)*E(4) <= Tol*Tol) {
 834
 835     // ns1, ns2 and  np are copied locally to avoid crushing the fields !
 836     gp_Vec ns1,ns2,np;
 837     ns1 = nsurf1;
 838     ns2 = nsurf2;
 839     np = nplan;
 840
 841     norm = nplan.Crossed(ns1).Magnitude();
 842     if (norm < Eps)  {
 843       norm = 1; // Unsatisfactory, but it is not necessary to stop
 844     }
 845     ns1.SetLinearForm(nplan.Dot(ns1)/norm,nplan, -1./norm, ns1);
 846
 847     norm = nplan.Crossed(ns2).Magnitude();
 848     if (norm < Eps)  {
 849       norm = 1; // Unsatisfactory, but it is not necessary to stop
 850     }
 851     ns2.SetLinearForm(nplan.Dot(ns2)/norm,nplan, -1./norm, ns2);
 852
 853     Standard_Real maxpiv = 1.e-9;
 854     math_Vector controle(1,4),solution(1,4), tolerances(1,4);
 855     GetTolerance(tolerances,Tol);
 856
 857     istangent = Standard_True;
 858     math_Gauss Resol(DEDX,maxpiv);
 859     if (Resol.IsDone()) {
 860       Resol.Solve(-DEDT,solution);
 861       istangent = Standard_False;
 862       controle = DEDT.Added(DEDX.Multiplied(solution));
 863       if(Abs(controle(1)) > tolerances(1) ||
 864          Abs(controle(2)) > tolerances(2) ||
 865          Abs(controle(3)) > tolerances(3) ||
 866          Abs(controle(4)) > tolerances(4)){
 867         istangent = Standard_True;
 868       }
 869     }
 870
 871     if (istangent) {
 872       math_SVD SingRS (DEDX);
 873       if (SingRS.IsDone()) {
 874         SingRS.Solve(-DEDT, solution, 1.e-6);
 875         istangent = Standard_False;
 876         controle = DEDT.Added(DEDX.Multiplied(solution));
 877         if(Abs(controle(1)) > tolerances(1) ||
 878            Abs(controle(2)) > tolerances(2) ||
 879            Abs(controle(3)) > tolerances(3) ||
 880            Abs(controle(4)) > tolerances(4)){
 881 #ifdef OCCT_DEBUG
 882         cout<<"Cheminement : echec calcul des derivees"<<endl;
 883 #endif
 884           istangent = Standard_True;
 885         }
 886       }
 887     }
 888
 889     if(!istangent){
 890       tg1.SetLinearForm(solution(1),d1u1,solution(2),d1v1);
 891       tg2.SetLinearForm(solution(3),d1u2,solution(4),d1v2);
 892       tg12d.SetCoord(solution(1),solution(2));
 893       tg22d.SetCoord(solution(3),solution(4));
 894     }
 895
 896   // update of maxang
 897
 898     if (ray1 > 0.) {
 899       ns1.Reverse();
 900     }
 901     if (ray2 >0.) {
 902       ns2.Reverse();
 903     }
 904     Cosa = ns1.Dot(ns2);
 905     Sina = np.Dot(ns1.Crossed(ns2));
 906     if (choix%2 != 0) {
 907       Sina = -Sina;  //nplan is changed in -nplan
 908     }
 909
 910     if(Cosa > 1.) {Cosa = 1.; Sina = 0.;}
 911     Angle = ACos(Cosa);
 912
 913  // Reframing on  ]-pi/2, 3pi/2]
 914     if (Sina <0.) {
 915       if (Cosa > 0.) Angle = -Angle;
 916       else           Angle =  2.*M_PI - Angle;
 917     }
 918
 919 //    cout << "Angle : " <<Angle << endl;
 920 //    if ((Angle < 0) || (Angle > M_PI)) {
 921 //      cout << "t = " << param << endl;
 922 //    }
 923
 924     if (Abs(Angle)>maxang) {maxang = Abs(Angle);}
 925     if (Abs(Angle)<minang) {minang = Abs(Angle);}
 926     distmin = Min( distmin, pts1.Distance(pts2));
 927
 928     return Ok;
 929   }
 930   istangent = Standard_True;
 931   return Standard_False;
 932 }
 933
 934 //=======================================================================
 935 //function : GetMinimalDistance
 936 //purpose  :
 937 //=======================================================================
 938
 939 Standard_Real BlendFunc_ConstRad::GetMinimalDistance() const
 940 {
 941   return distmin;
 942 }
 943
 944 //=======================================================================
 945 //function : Value
 946 //purpose  :
 947 //=======================================================================
 948
 949 Standard_Boolean BlendFunc_ConstRad::Value(const math_Vector& X, math_Vector& F)
 950 {
 951   const Standard_Boolean Ok = ComputeValues(X, 0);
 952   F = E;
 953   return Ok;
 954 }
 955
 956
 957
 958 //=======================================================================
 959 //function : Derivatives
 960 //purpose  :
 961 //=======================================================================
 962
 963 Standard_Boolean BlendFunc_ConstRad::Derivatives(const math_Vector& X, math_Matrix& D)
 964 {
 965   const Standard_Boolean Ok = ComputeValues(X, 1);
 966   D = DEDX;
 967   return Ok;
 968 }
 969
 970 //=======================================================================
 971 //function : Values
 972 //purpose  :
 973 //=======================================================================
 974
 975 Standard_Boolean BlendFunc_ConstRad::Values(const math_Vector& X, math_Vector& F, math_Matrix& D)
 976 {
 977   const Standard_Boolean Ok = ComputeValues(X, 1);
 978   F = E;
 979   D = DEDX;
 980   return Ok;
 981 }
 982
 983 //=======================================================================
 984 //function : PointOnS1
 985 //purpose  :
 986 //=======================================================================
 987
 988 const gp_Pnt& BlendFunc_ConstRad::PointOnS1 () const
 989 {
 990   return pts1;
 991 }
 992
 993
 994 //=======================================================================
 995 //function : PointOnS2
 996 //purpose  :
 997 //=======================================================================
 998
 999 const gp_Pnt& BlendFunc_ConstRad::PointOnS2 () const
1000 {
1001   return pts2;
1002 }
1003
1004
1005 //=======================================================================
1006 //function : IsTangencyPoint
1007 //purpose  :
1008 //=======================================================================
1009
1010 Standard_Boolean BlendFunc_ConstRad::IsTangencyPoint () const
1011 {
1012   return istangent;
1013 }
1014
1015
1016 //=======================================================================
1017 //function : TangentOnS1
1018 //purpose  :
1019 //=======================================================================
1020
1021 const gp_Vec& BlendFunc_ConstRad::TangentOnS1 () const
1022 {
1023   if (istangent)
1024     Standard_DomainError::Raise("BlendFunc_ConstRad::TangentOnS1");
1025   return tg1;
1026 }
1027
1028
1029 //=======================================================================
1030 //function : TangentOnS2
1031 //purpose  :
1032 //=======================================================================
1033
1034 const gp_Vec& BlendFunc_ConstRad::TangentOnS2 () const
1035 {
1036   if (istangent)
1037     Standard_DomainError::Raise("BlendFunc_ConstRad::TangentOnS2");
1038   return tg2;
1039 }
1040
1041
1042 //=======================================================================
1043 //function : Tangent2dOnS1
1044 //purpose  :
1045 //=======================================================================
1046
1047 const gp_Vec2d& BlendFunc_ConstRad::Tangent2dOnS1 () const
1048 {
1049   if (istangent)
1050     Standard_DomainError::Raise("BlendFunc_ConstRad::Tangent2dOnS1");
1051   return tg12d;
1052 }
1053
1054
1055 //=======================================================================
1056 //function : Tangent2dOnS2
1057 //purpose  :
1058 //=======================================================================
1059
1060 const gp_Vec2d& BlendFunc_ConstRad::Tangent2dOnS2 () const
1061 {
1062   if (istangent)
1063     Standard_DomainError::Raise("BlendFunc_ConstRad::Tangent2dOnS2");
1064   return tg22d;
1065 }
1066
1067
1068 //=======================================================================
1069 //function : Tangent
1070 //purpose  :
1071 //=======================================================================
1072
1073 void BlendFunc_ConstRad::Tangent(const Standard_Real U1,
1074                                  const Standard_Real V1,
1075                                  const Standard_Real U2,
1076                                  const Standard_Real V2,
1077                                  gp_Vec& TgF,
1078                                  gp_Vec& TgL,
1079                                  gp_Vec& NmF,
1080                                  gp_Vec& NmL) const
1081 {
1082   gp_Pnt Center;
1083   gp_Vec ns1;
1084   Standard_Real  invnorm1;
1085
1086   if ((U1!=xval(1)) || (V1!=xval(2)) ||
1087       (U2!=xval(3)) || (V2!=xval(4))) {
1088     gp_Vec d1u,d1v;
1089     gp_Pnt bid;
1090     surf1->D1(U1,V1,bid,d1u,d1v);
1091     NmF = ns1 = d1u.Crossed(d1v);
1092     surf2->D1(U2,V2,bid,d1u,d1v);
1093     NmL = d1u.Crossed(d1v);
1094   }
1095   else {
1096     NmF = ns1 = nsurf1;
1097     NmL = nsurf2;
1098   }
1099
1100   invnorm1 = nplan.Crossed(ns1).Magnitude();
1101   if (invnorm1 < Eps) invnorm1 = 1;
1102   else invnorm1 = 1. / invnorm1;
1103
1104   ns1.SetLinearForm(nplan.Dot(ns1)*invnorm1,nplan, -invnorm1,ns1);
1105   Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1106
1107   TgF = nplan.Crossed(gp_Vec(Center,pts1));
1108   TgL = nplan.Crossed(gp_Vec(Center,pts2));
1109   if (choix%2 == 1) {
1110     TgF.Reverse();
1111     TgL.Reverse();
1112   }
1113 }
1114
1115 //=======================================================================
1116 //function : TwistOnS1
1117 //purpose  :
1118 //=======================================================================
1119
1120 Standard_Boolean BlendFunc_ConstRad::TwistOnS1() const
1121 {
1122   if (istangent)
1123     Standard_DomainError::Raise("BlendFunc_ConstRad::TwistOnS1");
1124   return tg1.Dot(nplan) < 0.;
1125 }
1126
1127 //=======================================================================
1128 //function : TwistOnS2
1129 //purpose  :
1130 //=======================================================================
1131
1132 Standard_Boolean BlendFunc_ConstRad::TwistOnS2() const
1133 {
1134   if (istangent)
1135     Standard_DomainError::Raise("BlendFunc_ConstRad::TwistOnS2");
1136   return tg2.Dot(nplan) < 0.;
1137 }
1138
1139 //=======================================================================
1140 //function : Section
1141 //purpose  :
1142 //=======================================================================
1143
1144 void BlendFunc_ConstRad::Section(const Standard_Real Param,
1145                                  const Standard_Real U1,
1146                                  const Standard_Real V1,
1147                                  const Standard_Real U2,
1148                                  const Standard_Real V2,
1149                                  Standard_Real& Pdeb,
1150                                  Standard_Real& Pfin,
1151                                  gp_Circ& C)
1152 {
1153   gp_Pnt Center;
1154   gp_Vec ns1,np;
1155
1156   math_Vector X(1,4);
1157   X(1) = U1; X(2) = V1; X(3) = U2; X(4) = V2;
1158   Standard_Real prm = Param;
1159
1160   ComputeValues(X, 0, Standard_True, prm);
1161
1162   ns1 = nsurf1;
1163   np = nplan;
1164
1165   Standard_Real  norm1;
1166   norm1 = nplan.Crossed(ns1).Magnitude();
1167   if (norm1 < Eps)  {
1168     norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1169   }
1170   ns1.SetLinearForm(nplan.Dot(ns1)/norm1,nplan, -1./norm1,ns1);
1171   Center.SetXYZ(pts1.XYZ()+ray1*ns1.XYZ());
1172
1173 // ns1 is oriented from the center to pts1,
1174
1175   if (ray1 > 0.) {
1176     ns1.Reverse();
1177   }
1178   if (choix%2 != 0) {
1179     np.Reverse();
1180   }
1181   C.SetRadius(Abs(ray1));
1182   C.SetPosition(gp_Ax2(Center,np,ns1));
1183   Pdeb = 0.;
1184   Pfin = ElCLib::Parameter(C,pts2);
1185   // Test negative and almost null angles : Singular Case
1186   if (Pfin>1.5*M_PI) {
1187     np.Reverse();
1188     C.SetPosition(gp_Ax2(Center,np,ns1));
1189     Pfin = ElCLib::Parameter(C,pts2);
1190   }
1191   if (Pfin < Precision::PConfusion()) Pfin += Precision::PConfusion();
1192 }
1193
1194 //=======================================================================
1195 //function : IsRational
1196 //purpose  :
1197 //=======================================================================
1198
1199 Standard_Boolean BlendFunc_ConstRad::IsRational () const
1200 {
1201  return  (mySShape==BlendFunc_Rational || mySShape==BlendFunc_QuasiAngular);
1202 }
1203
1204 //=======================================================================
1205 //function : GetSectionSize
1206 //purpose  :
1207 //=======================================================================
1208 Standard_Real BlendFunc_ConstRad::GetSectionSize() const
1209 {
1210   return maxang*Abs(ray1);
1211 }
1212
1213 //=======================================================================
1214 //function : GetMinimalWeight
1215 //purpose  :
1216 //=======================================================================
1217 void BlendFunc_ConstRad::GetMinimalWeight(TColStd_Array1OfReal& Weigths) const
1218 {
1219   BlendFunc::GetMinimalWeights(mySShape, myTConv, minang, maxang, Weigths );
1220   // It is supposed that it does not depend on the Radius!
1221 }
1222
1223 //=======================================================================
1224 //function : NbIntervals
1225 //purpose  :
1226 //=======================================================================
1227
1228 Standard_Integer BlendFunc_ConstRad::NbIntervals (const GeomAbs_Shape S) const
1229 {
1230   return curv->NbIntervals(BlendFunc::NextShape(S));
1231 }
1232
1233
1234 //=======================================================================
1235 //function : Intervals
1236 //purpose  :
1237 //=======================================================================
1238
1239 void BlendFunc_ConstRad::Intervals (TColStd_Array1OfReal& T,
1240                                     const GeomAbs_Shape S) const
1241 {
1242   curv->Intervals(T, BlendFunc::NextShape(S));
1243 }
1244
1245 //=======================================================================
1246 //function : GetShape
1247 //purpose  :
1248 //=======================================================================
1249
1250 void BlendFunc_ConstRad::GetShape (Standard_Integer& NbPoles,
1251                                    Standard_Integer& NbKnots,
1252                                    Standard_Integer& Degree,
1253                                    Standard_Integer& NbPoles2d)
1254 {
1255   NbPoles2d = 2;
1256   BlendFunc::GetShape(mySShape,maxang,NbPoles,NbKnots,Degree,myTConv);
1257 }
1258
1259 //=======================================================================
1260 //function : GetTolerance
1261 //purpose  : Determine Tolerances used for approximations.
1262 //=======================================================================
1263 void BlendFunc_ConstRad::GetTolerance(const Standard_Real BoundTol,
1264                                       const Standard_Real SurfTol,
1265                                       const Standard_Real AngleTol,
1266                                       math_Vector& Tol3d,
1267                                       math_Vector& Tol1d) const
1268 {
1269  Standard_Integer low = Tol3d.Lower() , up=Tol3d.Upper();
1270  Standard_Real Tol;
1271  Tol= GeomFill::GetTolerance(myTConv, minang, Abs(ray1),
1272                              AngleTol, SurfTol);
1273  Tol1d.Init(SurfTol);
1274  Tol3d.Init(SurfTol);
1275  Tol3d(low+1) = Tol3d(up-1) = Min(Tol, SurfTol);
1276  Tol3d(low) = Tol3d(up) = Min(Tol, BoundTol);
1277
1278 }
1279
1280 //=======================================================================
1281 //function : Knots
1282 //purpose  :
1283 //=======================================================================
1284
1285 void BlendFunc_ConstRad::Knots(TColStd_Array1OfReal& TKnots)
1286 {
1287   GeomFill::Knots(myTConv,TKnots);
1288 }
1289
1290
1291 //=======================================================================
1292 //function : Mults
1293 //purpose  :
1294 //=======================================================================
1295
1296 void BlendFunc_ConstRad::Mults(TColStd_Array1OfInteger& TMults)
1297 {
1298   GeomFill::Mults(myTConv,TMults);
1299 }
1300
1301
1302 //=======================================================================
1303 //function : Section
1304 //purpose  :
1305 //=======================================================================
1306
1307 void BlendFunc_ConstRad::Section(const Blend_Point& P,
1308                                  TColgp_Array1OfPnt& Poles,
1309                                  TColgp_Array1OfPnt2d& Poles2d,
1310                                  TColStd_Array1OfReal& Weights)
1311 {
1312   gp_Pnt Center;
1313   gp_Vec ns1,ns2,np;
1314
1315   math_Vector X(1,4);
1316   Standard_Real prm = P.Parameter();
1317
1318   Standard_Integer low = Poles.Lower();
1319   Standard_Integer upp = Poles.Upper();
1320
1321   P.ParametersOnS1(X(1), X(2));
1322   P.ParametersOnS2(X(3), X(4));
1323
1324   ComputeValues(X, 0, Standard_True, prm);
1325   distmin = Min (distmin, pts1.Distance(pts2));
1326
1327   // ns1, ns2, np are copied locally to avoid crushing the fields !
1328   ns1 = nsurf1;
1329   ns2 = nsurf2;
1330   np = nplan;
1331
1332
1333   Poles2d(Poles2d.Lower()).SetCoord(X(1), X(2));
1334   Poles2d(Poles2d.Upper()).SetCoord(X(3), X(4));
1335
1336   if (mySShape == BlendFunc_Linear) {
1337     Poles(low) = pts1;
1338     Poles(upp) = pts2;
1339     Weights(low) = 1.0;
1340     Weights(upp) = 1.0;
1341     return;
1342   }
1343
1344   Standard_Real  norm1, norm2;
1345   norm1 = nplan.Crossed(ns1).Magnitude();
1346   norm2 = nplan.Crossed(ns2).Magnitude();
1347   if (norm1 < Eps)  {
1348     norm1 = 1; // Unsatisfactory, but it is not necessary to stop
1349   }
1350   if (norm2 < Eps) {
1351     norm2 = 1; // Unsatisfactory, but it is not necessary to stop
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 #ifdef OCCT_DEBUG
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 #ifdef OCCT_DEBUG
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 #ifdef OCCT_DEBUG
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 #ifdef OCCT_DEBUG
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 #ifdef OCCT_DEBUG
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 #ifdef OCCT_DEBUG
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 }