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
 135   // Calcul des coeffs  -----------
 136   beta = Angle/4;
 137   beta2 = beta * beta;
 138   beta3 = beta * beta2;
 139   beta4 = beta2*beta2;
 140   beta5 = beta3*beta2;
 141   beta6 = beta3*beta3;
 142
 143   if  ((M_PI/2 - beta)> NullAngle) {
 144     if (Abs(beta) < NullAngle) {
 145      Standard_Real cf = 2.0/(3*5*7);
 146      b = - (0.2+cf*beta2) / (1+ 0.2*beta2);
 147 //     b = beta5 / cf;
 148     }
 149     else {
 150       tan_b = Tan(beta);
 151       b = - 1.0e0 / beta2;
 152       b += beta / (3*(tan_b - beta));
 153     }
 154   }
 155   else b = ((Standard_Real) -1)/beta2;
 156   c =  ((Standard_Real) 1)/ 3 + b;
 157   b2 = b*b;
 158   c2 = c*c;
 159
 160   // X = U*U - V*V
 161   Vx(3) = beta2*(2*b - 1);
 162   Vx(5) = beta4*(b2 - 2*c);
 163   Vx(7) = -beta6*c2;
 164
 165   //Y = 2*U*V
 166   Vy(2) = 2*beta;
 167   Vy(4) = beta3*2*(c+b);
 168   Vy(6) = 2*beta5*b*c;
 169
 170  // W = U*U + V*V
 171   Vw(3) = beta2*(1 + 2*b);
 172   Vw(5) = beta4*(2*c + b2);
 173   Vw(7) = beta6*c2;
 174
 175 // Calculs des poles
 176   Px.Multiply(B, Vx);
 177   Py.Multiply(B, Vy);
 178   W.Multiply(B, Vw);
 179
 180 // Transfo
 181   gp_XYZ pnt;
 182   for (ii=1; ii<=7; ii++) {
 183      wi = W(ii);
 184      pnt.SetCoord(Px(ii)/wi, Py(ii)/wi, 0);
 185      pnt *= M;
 186      pnt += Center.XYZ();
 187      Poles(ii).ChangeCoord() = pnt;
 188      Weights(ii) = wi;
 189   }
 190 }
 191
 192 void GeomFill_QuasiAngularConvertor::Section(const gp_Pnt& FirstPnt,
 193                                        const gp_Vec& DFirstPnt,
 194                                        const gp_Pnt& Center,
 195                                        const gp_Vec& DCenter,
 196                                        const gp_Vec& Dir,
 197                                        const gp_Vec& DDir,
 198                                        const Standard_Real Angle,
 199                                        const Standard_Real DAngle,
 200                                        TColgp_Array1OfPnt& Poles,
 201                                        TColgp_Array1OfVec& DPoles,
 202                                        TColStd_Array1OfReal& Weights,
 203                                        TColStd_Array1OfReal& DWeights)
 204 {
 205   Standard_Integer Ordre = 7;
 206   math_Vector DVx(1, Ordre), DVy(1, Ordre), DVw(1, Ordre),
 207               DPx(1, Ordre), DPy(1, Ordre), DW(1, Ordre);
 208   Standard_Real b, b2, c, c2,tan_b;
 209   Standard_Real bpr, dtan_b;
 210   Standard_Integer ii;
 211   Standard_Real beta, beta2, beta3, beta4, beta5, beta6, betaprim;
 212   gp_Vec V1(Center, FirstPnt), V1Prim, V2;
 213
 214   // Calcul des  transformations
 215   gp_XYZ aux;
 216   Standard_Real Sina, Cosa;
 217   gp_Mat Rot, RotPrim, D, DPrim;
 218   // La rotation s'ecrit I +  sin(Ang) * D + (1. - cos(Ang)) * D*D
 219   // ou D est l'application x -> Dir ^ x
 220   Rot.SetRotation(Dir.XYZ(), Angle/2);
 221   // La derive s'ecrit donc :
 222   // AngPrim * (sin(Ang)*D*D + cos(Ang)*D)
 223   // + sin(Ang)*DPrim  + (1. - cos(Ang)) *(DPrim*D + D*DPrim)
 224   Sina = Sin(Angle/2);
 225   Cosa = Cos(Angle/2);
 226   D.SetCross(Dir.XYZ());
 227   DPrim.SetCross(DDir.XYZ());
 228
 229   RotPrim =  (D.Powered(2)).Multiplied(Sina);
 230   RotPrim += D.Multiplied(Cosa);
 231   RotPrim *= DAngle/2;
 232   RotPrim += DPrim.Multiplied(Sina);
 233   RotPrim += ((DPrim.Multiplied(D)).Added(D.Multiplied(DPrim))).Multiplied(1-Cosa);
 234
 235   aux = (DFirstPnt - DCenter).XYZ().Multiplied(Rot);
 236   aux += V1.XYZ().Multiplied(RotPrim);
 237   V1Prim.SetXYZ(aux);
 238   aux = V1.XYZ();
 239   aux.Multiply(Rot);
 240   V1.SetXYZ(aux);
 241   V2 =  Dir^V1;
 242   gp_Mat M (V1.X(), V2.X(), 0,
 243             V1.Y(), V2.Y(), 0,
 244             V1.Z(), V2.Z(), 0);
 245   V2 = (DDir.Crossed(V1)).Added(Dir.Crossed(V1Prim));
 246   gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
 247                 V1Prim.Y(), V2.Y(), 0,
 248                 V1Prim.Z(), V2.Z(), 0);
 249
 250   // Calcul des constante  -----------
 251   beta = Angle/4;
 252   betaprim = DAngle/4;
 253   beta2 = beta * beta;
 254   beta3 = beta * beta2;
 255   beta4 = beta2*beta2;
 256   beta5 = beta3*beta2;
 257   beta6 = beta3*beta3;
 258
 259   if (Abs(beta) < NullAngle) {
 260     // On calcul b par D.L
 261      Standard_Real cf = 2.0/(3*5*7);
 262      Standard_Real Num, Denom;
 263      Num = 0.2 + cf*beta2;
 264      Denom = 1+0.2*beta2;
 265      b = - Num/Denom;
 266      bpr = -2*beta*betaprim*(cf*Denom - 0.2*Num)/(Denom*Denom);
 267   }
 268   else {
 269     b = ((Standard_Real) -1)/beta2;
 270     bpr = (2*betaprim) / beta3;
 271     if  ((M_PI/2 - beta)> NullAngle) {
 272       tan_b = Tan(beta);
 273       dtan_b = betaprim * (1 + tan_b*tan_b);
 274       b2 = tan_b - beta;
 275       b += beta / (3*b2);
 276       bpr += (betaprim*tan_b - beta*dtan_b) / (3*b2*b2);
 277     }
 278   }
 279
 280   c =  ((Standard_Real) 1)/ 3 + b;
 281   b2 = b*b;
 282   c2 = c*c;
 283
 284
 285   // X = U*U - V*V
 286   Vx(3) = beta2*(2*b - 1);
 287   Vx(5) = beta4*(b2 - 2*c);
 288   Vx(7) = -beta6*c2;
 289   DVx.Init(0);
 290   DVx(3) = 2*(beta*betaprim*(2*b - 1) + bpr*beta2);
 291   DVx(5) = 4*beta3*betaprim*(b2 - 2*c) + 2*beta4*bpr*(b-1);
 292   DVx(7) = - 6*beta5*betaprim*c2 - 2*beta6*bpr*c;
 293
 294   //Y = 2*U*V
 295   Vy(2) = 2*beta;
 296   Vy(4) = beta3*2*(c+b);
 297   Vy(6) = 2*beta5*b*c;
 298   DVy.Init(0);
 299   DVy(2) = 2*betaprim;
 300   DVy(4) = 6*beta2*betaprim*(b+c) + 4*beta3*bpr;
 301   DVy(6) = 10*beta4*betaprim*b*c + 2*beta5*bpr*(b+c);
 302
 303  // W = U*U + V*V
 304   Vw(3) = beta2*(1 + 2*b);
 305   Vw(5) = beta4*(2*c + b2);
 306   Vw(7) = beta6*c2;
 307   DVw.Init(0);
 308 //  DVw(3) = 2*(beta*betaprim*(1 + 2*b) + beta2*bpr);
 309   DVw(3) = 2*beta*(betaprim*(1 + 2*b) + beta*bpr);
 310 //  DVw(5) = 4*beta3*betaprim*(2*c + b2) + 2*beta4*bpr*(b+1);
 311   DVw(5) = 2*beta3*(2*betaprim*(2*c + b2) + beta*bpr*(b+1));
 312 //  DVw(7) = 6*beta5*betaprim*c2 + 2*beta6*bpr*c;
 313   DVw(7) = 2*beta5*c*(3*betaprim*c + beta*bpr);
 314
 315
 316
 317   // Calcul des poles
 318   Px.Multiply(B, Vx);
 319   Py.Multiply(B, Vy);
 320   W.Multiply(B, Vw);
 321   DPx.Multiply(B, DVx);
 322   DPy.Multiply(B, DVy);
 323   DW.Multiply(B, DVw);
 324
 325   gp_XYZ P, DP;
 326   Standard_Real wi;
 327
 328   for (ii=1; ii<=Ordre; ii++) {
 329     wi = W(ii);
 330     P.SetCoord(Px(ii)/wi, Py(ii)/wi, 0);
 331     DP.SetCoord(DPx(ii)/wi, DPy(ii)/wi, 0);
 332     DP -= (DW(ii)/wi)*P;
 333
 334     Poles(ii).ChangeCoord() = M*P + Center.XYZ();
 335     P *= MPrim;
 336     DP *= M;
 337     aux.SetLinearForm(1, P, 1, DP, DCenter.XYZ());
 338     DPoles(ii).SetXYZ(aux);
 339     Weights(ii) = wi;
 340     DWeights(ii) = DW(ii);
 341   }
 342 }
 343
 344 void GeomFill_QuasiAngularConvertor::Section(const gp_Pnt& FirstPnt,
 345                                        const gp_Vec& DFirstPnt,
 346                                        const gp_Vec& D2FirstPnt,
 347                                        const gp_Pnt& Center,
 348                                        const gp_Vec& DCenter,
 349                                        const gp_Vec& D2Center,
 350                                        const gp_Vec& Dir,
 351                                        const gp_Vec& DDir,
 352                                        const gp_Vec& D2Dir,
 353                                        const Standard_Real Angle,
 354                                        const Standard_Real DAngle,
 355                                        const Standard_Real D2Angle,
 356                                        TColgp_Array1OfPnt& Poles,
 357                                        TColgp_Array1OfVec& DPoles,
 358                                        TColgp_Array1OfVec& D2Poles,
 359                                        TColStd_Array1OfReal& Weights,
 360                                        TColStd_Array1OfReal& DWeights,
 361                                        TColStd_Array1OfReal& D2Weights)
 362 {
 363   Standard_Integer Ordre = 7;
 364   math_Vector DVx(1, Ordre), DVy(1, Ordre),  DVw(1, Ordre),
 365               D2Vx(1, Ordre), D2Vy(1, Ordre),  D2Vw(1, Ordre);
 366   math_Vector DPx(1, Ordre), DPy(1, Ordre), DW(1, Ordre),
 367               D2Px(1, Ordre), D2Py(1, Ordre), D2W(1, Ordre);
 368
 369   Standard_Integer ii;
 370   Standard_Real aux, daux, b, b2, c, c2, bpr, bsc;
 371   gp_Vec V1(Center, FirstPnt), V1Prim, V1Secn, V2;
 372
 373   // Calcul des  transformations
 374   gp_XYZ auxyz;
 375   Standard_Real Sina, Cosa;
 376   gp_Mat Rot, RotPrim, RotSecn, D, DPrim, DSecn, DDP, Maux;
 377   // La rotation s'ecrit I +  sin(Ang) * D + (1. - cos(Ang)) * D*D
 378   // ou D est l'application x -> Dir ^ x
 379   Rot.SetRotation(Dir.XYZ(), Angle/2);
 380   // La derive s'ecrit donc :
 381   // AngPrim * (sin(Ang)*D*D + cos(Ang)*D)
 382   // + sin(Ang)*DPrim  + (1. - cos(Ang)) *(DPrim*D + D*DPrim)
 383   Sina = Sin(Angle/2);
 384   Cosa = Cos(Angle/2);
 385   D.SetCross(Dir.XYZ());
 386   DPrim.SetCross(DDir.XYZ());
 387   DSecn.SetCross(D2Dir.XYZ());
 388
 389   DDP = (DPrim.Multiplied(D)).Added(D.Multiplied(DPrim));
 390   RotPrim =  (D.Powered(2)).Multiplied(Sina);
 391   RotPrim += D.Multiplied(Cosa);
 392   RotPrim *= DAngle/2;
 393   RotPrim += DPrim.Multiplied(Sina);
 394   RotPrim += DDP.Multiplied(1-Cosa);
 395
 396   RotSecn =  (D.Powered(2)).Multiplied(Sina);
 397   RotSecn += D.Multiplied(Cosa);
 398   RotSecn *= D2Angle/2;
 399   Maux = (D.Powered(2)).Multiplied(Cosa);
 400   Maux -= D.Multiplied(Sina);
 401   Maux *= DAngle/2;
 402   Maux += DDP.Multiplied(2*Sina);
 403   Maux += DPrim.Multiplied(2*Cosa);
 404   Maux *= DAngle/2;
 405   RotSecn += Maux;
 406   Maux = (DSecn.Multiplied(D)).Added(D.Multiplied(DSecn));
 407   Maux +=  (DPrim.Powered(2)).Multiplied(2);
 408   Maux *= 1 - Cosa;
 409   Maux += DSecn.Multiplied(Sina);
 410   RotSecn += Maux;
 411
 412   V1Prim = DFirstPnt - DCenter;
 413   auxyz = (D2FirstPnt - D2Center).XYZ().Multiplied(Rot);
 414   auxyz += 2*(V1Prim.XYZ().Multiplied(RotPrim));
 415   auxyz += V1.XYZ().Multiplied(RotSecn);
 416   V1Secn.SetXYZ(auxyz);
 417   auxyz = V1Prim.XYZ().Multiplied(Rot);
 418   auxyz += V1.XYZ().Multiplied(RotPrim);
 419   V1Prim.SetXYZ(auxyz);
 420   auxyz = V1.XYZ();
 421   auxyz.Multiply(Rot);
 422   V1.SetXYZ(auxyz);
 423   V2 =  Dir^V1;
 424
 425   gp_Mat M (V1.X(), V2.X(), 0,
 426             V1.Y(), V2.Y(), 0,
 427             V1.Z(), V2.Z(), 0);
 428   V2 = (DDir.Crossed(V1)).Added(Dir.Crossed(V1Prim));
 429   gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
 430                 V1Prim.Y(), V2.Y(), 0,
 431                 V1Prim.Z(), V2.Z(), 0);
 432
 433   V2 = DDir^V1Prim;
 434   V2 *= 2;
 435   V2 += (D2Dir.Crossed(V1)).Added(Dir.Crossed(V1Secn));
 436   gp_Mat MSecn (V1Secn.X(), V2.X(), 0,
 437                 V1Secn.Y(), V2.Y(), 0,
 438                 V1Secn.Z(), V2.Z(), 0);
 439
 440
 441   // Calcul des coeff  -----------
 442   Standard_Real tan_b, dtan_b, d2tan_b;
 443   Standard_Real beta, beta2, beta3, beta4, beta5, beta6, betaprim, betasecn;
 444   Standard_Real betaprim2, bpr2;
 445   beta = Angle/4;
 446   betaprim = DAngle/4;
 447   betasecn = D2Angle/4;
 448   beta2 = beta * beta;
 449   beta3 = beta * beta2;
 450   beta4 = beta2*beta2;
 451   beta5 = beta3*beta2;
 452   beta6 = beta3*beta3;
 453   betaprim2 =  betaprim * betaprim;
 454
 455  if (Abs(beta) < NullAngle) {
 456     // On calcul b par D.L
 457      Standard_Real cf =-2.0/21;
 458      Standard_Real Num, Denom, aux;
 459      Num = 0.2 + cf*beta2;
 460      Denom = 1+0.2*beta2;
 461      aux = (cf*Denom - 0.2*Num)/(Denom*Denom);
 462      b = - Num/Denom;
 463      bpr = -2*beta*betaprim*aux;
 464      bsc = 2*aux*(betaprim2 + beta*betasecn - 2*beta*betaprim2);
 465
 466   }
 467   else {
 468     b = ((Standard_Real) -1)/beta2;
 469     bpr = (2*betaprim) / beta3;
 470     bsc = (2*betasecn - 6*betaprim*(betaprim/beta)) / beta3;
 471     if  ((M_PI/2 - beta)> NullAngle) {
 472       tan_b = Tan(beta);
 473       dtan_b = betaprim * (1 + tan_b*tan_b);
 474       d2tan_b = betasecn * (1 + tan_b*tan_b)
 475         + 2*betaprim * tan_b * dtan_b;
 476       b2 = tan_b - beta;
 477       b += beta / (3*b2);
 478       aux =  betaprim*tan_b - beta*dtan_b;
 479       bpr += aux / (3*b2*b2);
 480       daux =  betasecn*tan_b - beta*d2tan_b;
 481       bsc += (daux - 2*aux*betaprim*tan_b*tan_b/b2)/(3*b2*b2);
 482     }
 483   }
 484
 485   c =  ((Standard_Real) 1)/ 3 + b;
 486   b2 = b*b;
 487   c2 = c*c;
 488   bpr2 = bpr * bpr;
 489
 490
 491   // X = U*U - V*V
 492   Vx(3) = beta2*(2*b - 1);
 493   Vx(5) = beta4*(b2 - 2*c);
 494   Vx(7) = -beta6*c2;
 495   DVx.Init(0);
 496   DVx(3) = 2*(beta*betaprim*(2*b - 1) + bpr*beta2);
 497   DVx(5) = 4*beta3*betaprim*(b2 - 2*c) + 2*beta4*bpr*(b-1);
 498   DVx(7) = - 6*beta5*betaprim*c2 - 2*beta6*bpr*c;
 499   D2Vx.Init(0);
 500   D2Vx(3) = 2*((betaprim2+beta*betasecn)*(2*b - 1)
 501                + 8*beta*betaprim*bpr
 502                + bsc*beta2);
 503   D2Vx(5) =  4*(b2 - 2*c)*(3*beta2*betaprim2 + beta3*betasecn)
 504              + 16*beta3*betaprim*bpr*(b-1)
 505              + 2*beta4*(bsc*(b-1)+bpr2);
 506   D2Vx(7) = - 6 * c2 * (5*beta4*betaprim2+beta5*betasecn)
 507             - 24*beta5*betaprim*bpr*c
 508             - 2*beta6*(bsc*c + bpr2);
 509
 510   //Y = 2*U*V
 511   Vy(2) = 2*beta;
 512   Vy(4) = beta3*2*(c+b);
 513   Vy(6) = 2*beta5*b*c;
 514   DVy.Init(0);
 515   DVy(2) = 2*betaprim;
 516   DVy(4) = 6*beta2*betaprim*(b+c) + 4*beta3*bpr;
 517   DVy(6) = 10*beta4*betaprim*b*c + 2*beta5*bpr*(b+c);
 518   D2Vy.Init(0);
 519   D2Vy(2) = 2*betasecn;
 520   D2Vy(4) = 6*(b+c)*(2*beta*betaprim2 + beta2*betasecn)
 521           + 24*beta2*betaprim*bpr*(b+c)
 522           + 4*beta3*bsc;
 523   D2Vy(6) = 10*b*c*(4*beta3*betaprim2 + beta4*betasecn)
 524           + 40 * beta4*betaprim*bpr*(b+c)
 525           + 2*beta5*(bsc*(b+c)+ 2*bpr2);
 526
 527  // W = U*U + V*V
 528   Vw(3) = beta2*(1 + 2*b);
 529   Vw(5) = beta4*(2*c + b2);
 530   Vw(7) = beta6*c2;
 531   DVw.Init(0);
 532   DVw(3) = 2*(beta*betaprim*(1 + 2*b) + beta2*bpr);
 533   DVw(5) = 4*beta3*betaprim*(2*c + b2) + 2*beta4*bpr*(b+1);
 534   DVw(7) = 6*beta5*betaprim*c2 + 2*beta6*bpr*c;
 535   D2Vw.Init(0);
 536   D2Vw(3) = 2*((betaprim2+beta*betasecn)*(2*b + 1)
 537                + 8*beta*betaprim*bpr
 538                + bsc*beta2);
 539   D2Vw(5) = 4*(b2 + 2*c)*(3*beta2*betaprim2 + beta3*betasecn)
 540              + 16*beta3*betaprim*bpr*(b+11)
 541              + 2*beta4*(bsc*(b+1)+bpr2);
 542   D2Vw(7) = 6 * c2 * (5*beta4*betaprim2+beta5*betasecn)
 543             + 24*beta5*betaprim*bpr*c
 544             + 2*beta6*(bsc*c + bpr2);
 545
 546
 547   // Calcul des poles
 548   Px = B * Vx;
 549   Py = B * Vy;
 550   W.Multiply(B, Vw);
 551   DPx = B * DVx;
 552   DPy = B * DVy;
 553   DW.Multiply(B, DVw);
 554   D2Px = B * D2Vx;
 555   D2Py = B * D2Vy;
 556   D2W.Multiply(B, D2Vw);
 557
 558   gp_XYZ P, DP, D2P;
 559   Standard_Real wi, dwi;
 560   for (ii=1; ii<=Ordre; ii++) {
 561      wi = W(ii);
 562      dwi = DW(ii);
 563      P.SetCoord(Px(ii)/wi, Py(ii)/wi, 0);
 564      DP.SetCoord(DPx(ii)/wi, DPy(ii)/wi, 0);
 565
 566      D2P.SetCoord(D2Px(ii)/wi, D2Py(ii)/wi, 0);
 567      D2P -= 2*(dwi/wi)*DP;
 568      D2P += (2*Pow(dwi/wi, 2) - D2W(ii)/wi)*P;
 569      DP -= (DW(ii)/wi)*P;
 570
 571      Poles(ii).ChangeCoord() = M*P + Center.XYZ();
 572      auxyz.SetLinearForm(1, MPrim*P,
 573                          1, M*DP,
 574                          DCenter.XYZ());
 575      DPoles(ii).SetXYZ(auxyz);
 576      P  *= MSecn;
 577      DP *= MPrim;
 578      D2P*= M;
 579      auxyz.SetLinearForm(1, P,
 580                          2, DP,
 581                          1, D2P,
 582                          D2Center.XYZ());
 583      D2Poles(ii).SetXYZ(auxyz);
 584      Weights(ii) = wi;
 585      DWeights(ii) = dwi;
 586      D2Weights(ii) = D2W(ii);
 587   }
 588 }