src/GeomFill/GeomFill_QuasiAngularConvertor.cxx

   1 // File:        GeomFill_QuasiAngularConvertor.cxx
   2 // Created:     Wed Aug  6 09:31:38 1997
   3 // Author:      Philippe MANGIN
   4 //              <pmn@sgi29>
   5
   6
   7 #include <GeomFill_QuasiAngularConvertor.ixx>
   8
   9
  10 #include <PLib.hxx>
  11 #include <gp_Mat.hxx>
  12 #include <gp_Ax1.hxx>
  13
  14 #include <Convert_CompPolynomialToPoles.hxx>
  15
  16 #include <TColStd_Array1OfReal.hxx>
  17 #include <TColStd_HArray2OfReal.hxx>
  18 #include <TColStd_HArray1OfInteger.hxx>
  19 #include <TColStd_HArray1OfReal.hxx>
  20
  21 #define NullAngle 1.e-6
  22
  23 // QuasiAngular is rational definition of Cos(theta(t) and sin(theta)
  24 // on [-alpha, +alpha] with
  25 //                     2      2
  26 //                    U   -  V
  27 //  cos (theta(t)) = ----------
  28 //                     2      2
  29 //                    U   +  V
  30 //
  31
  32 //                      2 * U*V
  33 //   sin (theta(t)) = ----------
  34 //                      2      2
  35 //                     U   +  V
  36 //                     3
  37 //       V(t) = t + c t
  38 //                     2
  39 //       U(t) = 1 + b t
  40 //            1
  41 //       c = ---  + b
  42 //            3
  43 //
  44 //
  45 //            -1                     gamma
  46 //       b =---------  +   -----------------------
  47 //                 2
  48 //            gamma         3*(tang gamma - gamma)
  49 //
  50 //     with gamma = alpha / 2
  51
  52
  53 GeomFill_QuasiAngularConvertor::GeomFill_QuasiAngularConvertor():
  54                     myinit(Standard_False),
  55                     B(1, 7, 1, 7),
  56                     Px(1, 7), Py(1, 7),
  57                     W(1,7), Vx(1, 7), Vy(1, 7), Vw(1,7)
  58 {
  59 }
  60
  61 Standard_Boolean GeomFill_QuasiAngularConvertor::Initialized() const
  62 {
  63  return myinit;
  64 }
  65
  66 void GeomFill_QuasiAngularConvertor::Init()
  67 {
  68   if (myinit) return; //On n'initialise qu'une fois
  69   Standard_Integer ii, jj, Ordre=7;
  70   Standard_Real terme;
  71   TColStd_Array1OfReal Coeffs(1, Ordre*Ordre), TrueInter(1,2), Inter(1,2);
  72   Handle(TColStd_HArray2OfReal)
  73     Poles1d = new (TColStd_HArray2OfReal) (1, Ordre, 1, Ordre);
  74
  75   //Calcul de B
  76   Inter.SetValue(1, -1);
  77   Inter.SetValue(2,  1);
  78   TrueInter.SetValue(1, -1);
  79   TrueInter.SetValue(2,  1);
  80
  81   Coeffs.Init(0);
  82   for (ii=1; ii<=Ordre; ii++) {  Coeffs.SetValue(ii+(ii-1)*Ordre, 1); }
  83
  84    //Convertion
  85   Convert_CompPolynomialToPoles
  86          AConverter(Ordre, Ordre-1,  Ordre-1,
  87                     Coeffs,
  88                     Inter,
  89                     TrueInter);
  90   AConverter.Poles(Poles1d);
  91
  92   for (jj=1; jj<=Ordre; jj++) {
  93       for (ii=1; ii<=Ordre; ii++) {
  94         terme = Poles1d->Value(ii,jj);
  95         if (Abs(terme-1) < 1.e-9) terme = 1 ; //petite retouche
  96         if (Abs(terme+1) < 1.e-9) terme = -1;
  97         B(ii, jj) =  terme;
  98       }
  99   }
 100
 101   // Init des polynomes
 102   Vx.Init(0);
 103   Vx(1) = 1;
 104   Vy.Init(0);
 105   Vy(2) = 2;
 106   Vw.Init(0);
 107   Vw(1) = 1;
 108   myinit = Standard_True;
 109 }
 110
 111 void GeomFill_QuasiAngularConvertor::Section(const gp_Pnt& FirstPnt,
 112                                        const gp_Pnt& Center,
 113                                        const gp_Vec& Dir,
 114                                        const Standard_Real Angle,
 115                                        TColgp_Array1OfPnt& Poles,
 116                                        TColStd_Array1OfReal& Weights)
 117 {
 118   Standard_Real b, b2, c, c2,tan_b;
 119   Standard_Integer ii;
 120   Standard_Real beta, beta2, beta3, beta4, beta5, beta6, wi;
 121   gp_XYZ aux;
 122   gp_Mat Rot;
 123   // Calcul de la transformation
 124   gp_Vec V1(Center, FirstPnt), V2;
 125   Rot.SetRotation(Dir.XYZ(), Angle/2);
 126   aux = V1.XYZ();
 127   aux *= Rot;
 128   V1.SetXYZ(aux);
 129   V2 =  Dir^V1;
 130
 131   gp_Mat M(V1.X(), V2.X(), 0,
 132            V1.Y(), V2.Y(), 0,
 133            V1.Z(), V2.Z(), 0);
 134 #ifdef DEB
 135   Standard_Real r = FirstPnt.Distance(Center);
 136 #else
 137   FirstPnt.Distance(Center);
 138 #endif
 139
 140   // Calcul des coeffs  -----------
 141   beta = Angle/4;
 142   beta2 = beta * beta;
 143   beta3 = beta * beta2;
 144   beta4 = beta2*beta2;
 145   beta5 = beta3*beta2;
 146   beta6 = beta3*beta3;
 147
 148   if  ((PI/2 - beta)> NullAngle) {
 149     if (Abs(beta) < NullAngle) {
 150      Standard_Real cf = 2.0/(3*5*7);
 151      b = - (0.2+cf*beta2) / (1+ 0.2*beta2);
 152 //     b = beta5 / cf;
 153     }
 154     else {
 155       tan_b = Tan(beta);
 156       b = - 1.0e0 / beta2;
 157       b += beta / (3*(tan_b - beta));
 158     }
 159   }
 160   else b = ((Standard_Real) -1)/beta2;
 161   c =  ((Standard_Real) 1)/ 3 + b;
 162   b2 = b*b;
 163   c2 = c*c;
 164
 165   // X = U*U - V*V
 166   Vx(3) = beta2*(2*b - 1);
 167   Vx(5) = beta4*(b2 - 2*c);
 168   Vx(7) = -beta6*c2;
 169
 170   //Y = 2*U*V
 171   Vy(2) = 2*beta;
 172   Vy(4) = beta3*2*(c+b);
 173   Vy(6) = 2*beta5*b*c;
 174
 175  // W = U*U + V*V
 176   Vw(3) = beta2*(1 + 2*b);
 177   Vw(5) = beta4*(2*c + b2);
 178   Vw(7) = beta6*c2;
 179
 180 // Calculs des poles
 181   Px.Multiply(B, Vx);
 182   Py.Multiply(B, Vy);
 183   W.Multiply(B, Vw);
 184
 185 // Transfo
 186   gp_XYZ pnt;
 187   for (ii=1; ii<=7; ii++) {
 188      wi = W(ii);
 189      pnt.SetCoord(Px(ii)/wi, Py(ii)/wi, 0);
 190      pnt *= M;
 191      pnt += Center.XYZ();
 192      Poles(ii).ChangeCoord() = pnt;
 193      Weights(ii) = wi;
 194   }
 195 }
 196
 197 void GeomFill_QuasiAngularConvertor::Section(const gp_Pnt& FirstPnt,
 198                                        const gp_Vec& DFirstPnt,
 199                                        const gp_Pnt& Center,
 200                                        const gp_Vec& DCenter,
 201                                        const gp_Vec& Dir,
 202                                        const gp_Vec& DDir,
 203                                        const Standard_Real Angle,
 204                                        const Standard_Real DAngle,
 205                                        TColgp_Array1OfPnt& Poles,
 206                                        TColgp_Array1OfVec& DPoles,
 207                                        TColStd_Array1OfReal& Weights,
 208                                        TColStd_Array1OfReal& DWeights)
 209 {
 210   Standard_Integer Ordre = 7;
 211   math_Vector DVx(1, Ordre), DVy(1, Ordre), DVw(1, Ordre),
 212               DPx(1, Ordre), DPy(1, Ordre), DW(1, Ordre);
 213   Standard_Real b, b2, c, c2,tan_b;
 214   Standard_Real bpr, dtan_b;
 215   Standard_Integer ii;
 216   Standard_Real beta, beta2, beta3, beta4, beta5, beta6, betaprim;
 217   gp_Vec V1(Center, FirstPnt), V1Prim, V2;
 218
 219   // Calcul des  transformations
 220   gp_XYZ aux;
 221   Standard_Real Sina, Cosa;
 222   gp_Mat Rot, RotPrim, D, DPrim;
 223   // La rotation s'ecrit I +  sin(Ang) * D + (1. - cos(Ang)) * D*D
 224   // ou D est l'application x -> Dir ^ x
 225   Rot.SetRotation(Dir.XYZ(), Angle/2);
 226   // La derive s'ecrit donc :
 227   // AngPrim * (sin(Ang)*D*D + cos(Ang)*D)
 228   // + sin(Ang)*DPrim  + (1. - cos(Ang)) *(DPrim*D + D*DPrim)
 229   Sina = Sin(Angle/2);
 230   Cosa = Cos(Angle/2);
 231   D.SetCross(Dir.XYZ());
 232   DPrim.SetCross(DDir.XYZ());
 233
 234   RotPrim =  (D.Powered(2)).Multiplied(Sina);
 235   RotPrim += D.Multiplied(Cosa);
 236   RotPrim *= DAngle/2;
 237   RotPrim += DPrim.Multiplied(Sina);
 238   RotPrim += ((DPrim.Multiplied(D)).Added(D.Multiplied(DPrim))).Multiplied(1-Cosa);
 239
 240   aux = (DFirstPnt - DCenter).XYZ().Multiplied(Rot);
 241   aux += V1.XYZ().Multiplied(RotPrim);
 242   V1Prim.SetXYZ(aux);
 243   aux = V1.XYZ();
 244   aux.Multiply(Rot);
 245   V1.SetXYZ(aux);
 246   V2 =  Dir^V1;
 247   gp_Mat M (V1.X(), V2.X(), 0,
 248             V1.Y(), V2.Y(), 0,
 249             V1.Z(), V2.Z(), 0);
 250   V2 = (DDir.Crossed(V1)).Added(Dir.Crossed(V1Prim));
 251   gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
 252                 V1Prim.Y(), V2.Y(), 0,
 253                 V1Prim.Z(), V2.Z(), 0);
 254
 255   // Calcul des constante  -----------
 256   beta = Angle/4;
 257   betaprim = DAngle/4;
 258   beta2 = beta * beta;
 259   beta3 = beta * beta2;
 260   beta4 = beta2*beta2;
 261   beta5 = beta3*beta2;
 262   beta6 = beta3*beta3;
 263
 264   if (Abs(beta) < NullAngle) {
 265     // On calcul b par D.L
 266      Standard_Real cf = 2.0/(3*5*7);
 267      Standard_Real Num, Denom;
 268      Num = 0.2 + cf*beta2;
 269      Denom = 1+0.2*beta2;
 270      b = - Num/Denom;
 271      bpr = -2*beta*betaprim*(cf*Denom - 0.2*Num)/(Denom*Denom);
 272   }
 273   else {
 274     b = ((Standard_Real) -1)/beta2;
 275     bpr = (2*betaprim) / beta3;
 276     if  ((PI/2 - beta)> NullAngle) {
 277       tan_b = Tan(beta);
 278       dtan_b = betaprim * (1 + tan_b*tan_b);
 279       b2 = tan_b - beta;
 280       b += beta / (3*b2);
 281       bpr += (betaprim*tan_b - beta*dtan_b) / (3*b2*b2);
 282     }
 283   }
 284
 285   c =  ((Standard_Real) 1)/ 3 + b;
 286   b2 = b*b;
 287   c2 = c*c;
 288
 289
 290   // X = U*U - V*V
 291   Vx(3) = beta2*(2*b - 1);
 292   Vx(5) = beta4*(b2 - 2*c);
 293   Vx(7) = -beta6*c2;
 294   DVx.Init(0);
 295   DVx(3) = 2*(beta*betaprim*(2*b - 1) + bpr*beta2);
 296   DVx(5) = 4*beta3*betaprim*(b2 - 2*c) + 2*beta4*bpr*(b-1);
 297   DVx(7) = - 6*beta5*betaprim*c2 - 2*beta6*bpr*c;
 298
 299   //Y = 2*U*V
 300   Vy(2) = 2*beta;
 301   Vy(4) = beta3*2*(c+b);
 302   Vy(6) = 2*beta5*b*c;
 303   DVy.Init(0);
 304   DVy(2) = 2*betaprim;
 305   DVy(4) = 6*beta2*betaprim*(b+c) + 4*beta3*bpr;
 306   DVy(6) = 10*beta4*betaprim*b*c + 2*beta5*bpr*(b+c);
 307
 308  // W = U*U + V*V
 309   Vw(3) = beta2*(1 + 2*b);
 310   Vw(5) = beta4*(2*c + b2);
 311   Vw(7) = beta6*c2;
 312   DVw.Init(0);
 313 //  DVw(3) = 2*(beta*betaprim*(1 + 2*b) + beta2*bpr);
 314   DVw(3) = 2*beta*(betaprim*(1 + 2*b) + beta*bpr);
 315 //  DVw(5) = 4*beta3*betaprim*(2*c + b2) + 2*beta4*bpr*(b+1);
 316   DVw(5) = 2*beta3*(2*betaprim*(2*c + b2) + beta*bpr*(b+1));
 317 //  DVw(7) = 6*beta5*betaprim*c2 + 2*beta6*bpr*c;
 318   DVw(7) = 2*beta5*c*(3*betaprim*c + beta*bpr);
 319
 320
 321
 322   // Calcul des poles
 323   Px.Multiply(B, Vx);
 324   Py.Multiply(B, Vy);
 325   W.Multiply(B, Vw);
 326   DPx.Multiply(B, DVx);
 327   DPy.Multiply(B, DVy);
 328   DW.Multiply(B, DVw);
 329
 330   gp_XYZ P, DP;
 331   Standard_Real wi;
 332
 333   for (ii=1; ii<=Ordre; ii++) {
 334     wi = W(ii);
 335     P.SetCoord(Px(ii)/wi, Py(ii)/wi, 0);
 336     DP.SetCoord(DPx(ii)/wi, DPy(ii)/wi, 0);
 337     DP -= (DW(ii)/wi)*P;
 338
 339     Poles(ii).ChangeCoord() = M*P + Center.XYZ();
 340     P *= MPrim;
 341     DP *= M;
 342     aux.SetLinearForm(1, P, 1, DP, DCenter.XYZ());
 343     DPoles(ii).SetXYZ(aux);
 344     Weights(ii) = wi;
 345     DWeights(ii) = DW(ii);
 346   }
 347 }
 348
 349 void GeomFill_QuasiAngularConvertor::Section(const gp_Pnt& FirstPnt,
 350                                        const gp_Vec& DFirstPnt,
 351                                        const gp_Vec& D2FirstPnt,
 352                                        const gp_Pnt& Center,
 353                                        const gp_Vec& DCenter,
 354                                        const gp_Vec& D2Center,
 355                                        const gp_Vec& Dir,
 356                                        const gp_Vec& DDir,
 357                                        const gp_Vec& D2Dir,
 358                                        const Standard_Real Angle,
 359                                        const Standard_Real DAngle,
 360                                        const Standard_Real D2Angle,
 361                                        TColgp_Array1OfPnt& Poles,
 362                                        TColgp_Array1OfVec& DPoles,
 363                                        TColgp_Array1OfVec& D2Poles,
 364                                        TColStd_Array1OfReal& Weights,
 365                                        TColStd_Array1OfReal& DWeights,
 366                                        TColStd_Array1OfReal& D2Weights)
 367 {
 368   Standard_Integer Ordre = 7;
 369   math_Vector DVx(1, Ordre), DVy(1, Ordre),  DVw(1, Ordre),
 370               D2Vx(1, Ordre), D2Vy(1, Ordre),  D2Vw(1, Ordre);
 371   math_Vector DPx(1, Ordre), DPy(1, Ordre), DW(1, Ordre),
 372               D2Px(1, Ordre), D2Py(1, Ordre), D2W(1, Ordre);
 373
 374   Standard_Integer ii;
 375   Standard_Real aux, daux, b, b2, c, c2, bpr, bsc;
 376   gp_Vec V1(Center, FirstPnt), V1Prim, V1Secn, V2;
 377
 378   // Calcul des  transformations
 379   gp_XYZ auxyz;
 380   Standard_Real Sina, Cosa;
 381   gp_Mat Rot, RotPrim, RotSecn, D, DPrim, DSecn, DDP, Maux;
 382   // La rotation s'ecrit I +  sin(Ang) * D + (1. - cos(Ang)) * D*D
 383   // ou D est l'application x -> Dir ^ x
 384   Rot.SetRotation(Dir.XYZ(), Angle/2);
 385   // La derive s'ecrit donc :
 386   // AngPrim * (sin(Ang)*D*D + cos(Ang)*D)
 387   // + sin(Ang)*DPrim  + (1. - cos(Ang)) *(DPrim*D + D*DPrim)
 388   Sina = Sin(Angle/2);
 389   Cosa = Cos(Angle/2);
 390   D.SetCross(Dir.XYZ());
 391   DPrim.SetCross(DDir.XYZ());
 392   DSecn.SetCross(D2Dir.XYZ());
 393
 394   DDP = (DPrim.Multiplied(D)).Added(D.Multiplied(DPrim));
 395   RotPrim =  (D.Powered(2)).Multiplied(Sina);
 396   RotPrim += D.Multiplied(Cosa);
 397   RotPrim *= DAngle/2;
 398   RotPrim += DPrim.Multiplied(Sina);
 399   RotPrim += DDP.Multiplied(1-Cosa);
 400
 401   RotSecn =  (D.Powered(2)).Multiplied(Sina);
 402   RotSecn += D.Multiplied(Cosa);
 403   RotSecn *= D2Angle/2;
 404   Maux = (D.Powered(2)).Multiplied(Cosa);
 405   Maux -= D.Multiplied(Sina);
 406   Maux *= DAngle/2;
 407   Maux += DDP.Multiplied(2*Sina);
 408   Maux += DPrim.Multiplied(2*Cosa);
 409   Maux *= DAngle/2;
 410   RotSecn += Maux;
 411   Maux = (DSecn.Multiplied(D)).Added(D.Multiplied(DSecn));
 412   Maux +=  (DPrim.Powered(2)).Multiplied(2);
 413   Maux *= 1 - Cosa;
 414   Maux += DSecn.Multiplied(Sina);
 415   RotSecn += Maux;
 416
 417   V1Prim = DFirstPnt - DCenter;
 418   auxyz = (D2FirstPnt - D2Center).XYZ().Multiplied(Rot);
 419   auxyz += 2*(V1Prim.XYZ().Multiplied(RotPrim));
 420   auxyz += V1.XYZ().Multiplied(RotSecn);
 421   V1Secn.SetXYZ(auxyz);
 422   auxyz = V1Prim.XYZ().Multiplied(Rot);
 423   auxyz += V1.XYZ().Multiplied(RotPrim);
 424   V1Prim.SetXYZ(auxyz);
 425   auxyz = V1.XYZ();
 426   auxyz.Multiply(Rot);
 427   V1.SetXYZ(auxyz);
 428   V2 =  Dir^V1;
 429
 430   gp_Mat M (V1.X(), V2.X(), 0,
 431             V1.Y(), V2.Y(), 0,
 432             V1.Z(), V2.Z(), 0);
 433   V2 = (DDir.Crossed(V1)).Added(Dir.Crossed(V1Prim));
 434   gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
 435                 V1Prim.Y(), V2.Y(), 0,
 436                 V1Prim.Z(), V2.Z(), 0);
 437
 438   V2 = DDir^V1Prim;
 439   V2 *= 2;
 440   V2 += (D2Dir.Crossed(V1)).Added(Dir.Crossed(V1Secn));
 441   gp_Mat MSecn (V1Secn.X(), V2.X(), 0,
 442                 V1Secn.Y(), V2.Y(), 0,
 443                 V1Secn.Z(), V2.Z(), 0);
 444
 445
 446   // Calcul des coeff  -----------
 447   Standard_Real tan_b, dtan_b, d2tan_b;
 448   Standard_Real beta, beta2, beta3, beta4, beta5, beta6, betaprim, betasecn;
 449   Standard_Real betaprim2, bpr2;
 450   beta = Angle/4;
 451   betaprim = DAngle/4;
 452   betasecn = D2Angle/4;
 453   beta2 = beta * beta;
 454   beta3 = beta * beta2;
 455   beta4 = beta2*beta2;
 456   beta5 = beta3*beta2;
 457   beta6 = beta3*beta3;
 458   betaprim2 =  betaprim * betaprim;
 459
 460  if (Abs(beta) < NullAngle) {
 461     // On calcul b par D.L
 462      Standard_Real cf =-2.0/21;
 463      Standard_Real Num, Denom, aux;
 464      Num = 0.2 + cf*beta2;
 465      Denom = 1+0.2*beta2;
 466      aux = (cf*Denom - 0.2*Num)/(Denom*Denom);
 467      b = - Num/Denom;
 468      bpr = -2*beta*betaprim*aux;
 469      bsc = 2*aux*(betaprim2 + beta*betasecn - 2*beta*betaprim2);
 470
 471   }
 472   else {
 473     b = ((Standard_Real) -1)/beta2;
 474     bpr = (2*betaprim) / beta3;
 475     bsc = (2*betasecn - 6*betaprim*(betaprim/beta)) / beta3;
 476     if  ((PI/2 - beta)> NullAngle) {
 477       tan_b = Tan(beta);
 478       dtan_b = betaprim * (1 + tan_b*tan_b);
 479       d2tan_b = betasecn * (1 + tan_b*tan_b)
 480         + 2*betaprim * tan_b * dtan_b;
 481       b2 = tan_b - beta;
 482       b += beta / (3*b2);
 483       aux =  betaprim*tan_b - beta*dtan_b;
 484       bpr += aux / (3*b2*b2);
 485       daux =  betasecn*tan_b - beta*d2tan_b;
 486       bsc += (daux - 2*aux*betaprim*tan_b*tan_b/b2)/(3*b2*b2);
 487     }
 488   }
 489
 490   c =  ((Standard_Real) 1)/ 3 + b;
 491   b2 = b*b;
 492   c2 = c*c;
 493   bpr2 = bpr * bpr;
 494
 495
 496   // X = U*U - V*V
 497   Vx(3) = beta2*(2*b - 1);
 498   Vx(5) = beta4*(b2 - 2*c);
 499   Vx(7) = -beta6*c2;
 500   DVx.Init(0);
 501   DVx(3) = 2*(beta*betaprim*(2*b - 1) + bpr*beta2);
 502   DVx(5) = 4*beta3*betaprim*(b2 - 2*c) + 2*beta4*bpr*(b-1);
 503   DVx(7) = - 6*beta5*betaprim*c2 - 2*beta6*bpr*c;
 504   D2Vx.Init(0);
 505   D2Vx(3) = 2*((betaprim2+beta*betasecn)*(2*b - 1)
 506                + 8*beta*betaprim*bpr
 507                + bsc*beta2);
 508   D2Vx(5) =  4*(b2 - 2*c)*(3*beta2*betaprim2 + beta3*betasecn)
 509              + 16*beta3*betaprim*bpr*(b-1)
 510              + 2*beta4*(bsc*(b-1)+bpr2);
 511   D2Vx(7) = - 6 * c2 * (5*beta4*betaprim2+beta5*betasecn)
 512             - 24*beta5*betaprim*bpr*c
 513             - 2*beta6*(bsc*c + bpr2);
 514
 515   //Y = 2*U*V
 516   Vy(2) = 2*beta;
 517   Vy(4) = beta3*2*(c+b);
 518   Vy(6) = 2*beta5*b*c;
 519   DVy.Init(0);
 520   DVy(2) = 2*betaprim;
 521   DVy(4) = 6*beta2*betaprim*(b+c) + 4*beta3*bpr;
 522   DVy(6) = 10*beta4*betaprim*b*c + 2*beta5*bpr*(b+c);
 523   D2Vy.Init(0);
 524   D2Vy(2) = 2*betasecn;
 525   D2Vy(4) = 6*(b+c)*(2*beta*betaprim2 + beta2*betasecn)
 526           + 24*beta2*betaprim*bpr*(b+c)
 527           + 4*beta3*bsc;
 528   D2Vy(6) = 10*b*c*(4*beta3*betaprim2 + beta4*betasecn)
 529           + 40 * beta4*betaprim*bpr*(b+c)
 530           + 2*beta5*(bsc*(b+c)+ 2*bpr2);
 531
 532  // W = U*U + V*V
 533   Vw(3) = beta2*(1 + 2*b);
 534   Vw(5) = beta4*(2*c + b2);
 535   Vw(7) = beta6*c2;
 536   DVw.Init(0);
 537   DVw(3) = 2*(beta*betaprim*(1 + 2*b) + beta2*bpr);
 538   DVw(5) = 4*beta3*betaprim*(2*c + b2) + 2*beta4*bpr*(b+1);
 539   DVw(7) = 6*beta5*betaprim*c2 + 2*beta6*bpr*c;
 540   D2Vw.Init(0);
 541   D2Vw(3) = 2*((betaprim2+beta*betasecn)*(2*b + 1)
 542                + 8*beta*betaprim*bpr
 543                + bsc*beta2);
 544   D2Vw(5) = 4*(b2 + 2*c)*(3*beta2*betaprim2 + beta3*betasecn)
 545              + 16*beta3*betaprim*bpr*(b+11)
 546              + 2*beta4*(bsc*(b+1)+bpr2);
 547   D2Vw(7) = 6 * c2 * (5*beta4*betaprim2+beta5*betasecn)
 548             + 24*beta5*betaprim*bpr*c
 549             + 2*beta6*(bsc*c + bpr2);
 550
 551
 552   // Calcul des poles
 553   Px = B * Vx;
 554   Py = B * Vy;
 555   W.Multiply(B, Vw);
 556   DPx = B * DVx;
 557   DPy = B * DVy;
 558   DW.Multiply(B, DVw);
 559   D2Px = B * D2Vx;
 560   D2Py = B * D2Vy;
 561   D2W.Multiply(B, D2Vw);
 562
 563   gp_XYZ P, DP, D2P;
 564   Standard_Real wi, dwi;
 565   for (ii=1; ii<=Ordre; ii++) {
 566      wi = W(ii);
 567      dwi = DW(ii);
 568      P.SetCoord(Px(ii)/wi, Py(ii)/wi, 0);
 569      DP.SetCoord(DPx(ii)/wi, DPy(ii)/wi, 0);
 570
 571      D2P.SetCoord(D2Px(ii)/wi, D2Py(ii)/wi, 0);
 572      D2P -= 2*(dwi/wi)*DP;
 573      D2P += (2*Pow(dwi/wi, 2) - D2W(ii)/wi)*P;
 574      DP -= (DW(ii)/wi)*P;
 575
 576      Poles(ii).ChangeCoord() = M*P + Center.XYZ();
 577      auxyz.SetLinearForm(1, MPrim*P,
 578                          1, M*DP,
 579                          DCenter.XYZ());
 580      DPoles(ii).SetXYZ(auxyz);
 581      P  *= MSecn;
 582      DP *= MPrim;
 583      D2P*= M;
 584      auxyz.SetLinearForm(1, P,
 585                          2, DP,
 586                          1, D2P,
 587                          D2Center.XYZ());
 588      D2Poles(ii).SetXYZ(auxyz);
 589      Weights(ii) = wi;
 590      DWeights(ii) = dwi;
 591      D2Weights(ii) = D2W(ii);
 592   }
 593 }