src/GeomFill/GeomFill_QuasiAngularConvertor.cxx

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