1 // File:        GeomFill_QuasiAngularConvertor.cxx
2 // Created:     Wed Aug  6 09:31:38 1997
3 // Author:      Philippe MANGIN
4 //              <pmn@sgi29>
7 #include <GeomFill_QuasiAngularConvertor.ixx>
10 #include <PLib.hxx>
11 #include <gp_Mat.hxx>
12 #include <gp_Ax1.hxx>
14 #include <Convert_CompPolynomialToPoles.hxx>
16 #include <TColStd_Array1OfReal.hxx>
17 #include <TColStd_HArray2OfReal.hxx>
18 #include <TColStd_HArray1OfInteger.hxx>
19 #include <TColStd_HArray1OfReal.hxx>
21 #define NullAngle 1.e-6
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 //
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
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 }
61 Standard_Boolean GeomFill_QuasiAngularConvertor::Initialized() const
62 {
63  return myinit;
64 }
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);
75   //Calcul de B
76   Inter.SetValue(1, -1);
77   Inter.SetValue(2,  1);
78   TrueInter.SetValue(1, -1);
79   TrueInter.SetValue(2,  1);
81   Coeffs.Init(0);
82   for (ii=1; ii<=Ordre; ii++) {  Coeffs.SetValue(ii+(ii-1)*Ordre, 1); }
84    //Convertion
85   Convert_CompPolynomialToPoles
86          AConverter(Ordre, Ordre-1,  Ordre-1,
87                     Coeffs,
88                     Inter,
89                     TrueInter);
90   AConverter.Poles(Poles1d);
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   }
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 }
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;
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
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;
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;
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;
170   //Y = 2*U*V
171   Vy(2) = 2*beta;
172   Vy(4) = beta3*2*(c+b);
173   Vy(6) = 2*beta5*b*c;
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;
180 // Calculs des poles
181   Px.Multiply(B, Vx);
182   Py.Multiply(B, Vy);
183   W.Multiply(B, Vw);
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 }
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;
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());
234   RotPrim =  (D.Powered(2)).Multiplied(Sina);
235   RotPrim += D.Multiplied(Cosa);
236   RotPrim *= DAngle/2;
237   RotPrim += DPrim.Multiplied(Sina);
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);
251   gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
252                 V1Prim.Y(), V2.Y(), 0,
253                 V1Prim.Z(), V2.Z(), 0);
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;
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   }
285   c =  ((Standard_Real) 1)/ 3 + b;
286   b2 = b*b;
287   c2 = c*c;
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;
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);
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);
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);
330   gp_XYZ P, DP;
331   Standard_Real wi;
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;
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 }
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);
374   Standard_Integer ii;
375   Standard_Real aux, daux, b, b2, c, c2, bpr, bsc;
376   gp_Vec V1(Center, FirstPnt), V1Prim, V1Secn, V2;
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());
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);
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;
412   Maux +=  (DPrim.Powered(2)).Multiplied(2);
413   Maux *= 1 - Cosa;
414   Maux += DSecn.Multiplied(Sina);
415   RotSecn += Maux;
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;
430   gp_Mat M (V1.X(), V2.X(), 0,
431             V1.Y(), V2.Y(), 0,
432             V1.Z(), V2.Z(), 0);
434   gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
435                 V1Prim.Y(), V2.Y(), 0,
436                 V1Prim.Z(), V2.Z(), 0);
438   V2 = DDir^V1Prim;
439   V2 *= 2;
441   gp_Mat MSecn (V1Secn.X(), V2.X(), 0,
442                 V1Secn.Y(), V2.Y(), 0,
443                 V1Secn.Z(), V2.Z(), 0);
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;
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);
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   }
490   c =  ((Standard_Real) 1)/ 3 + b;
491   b2 = b*b;
492   c2 = c*c;
493   bpr2 = bpr * bpr;
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);
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);
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);
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);
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);
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;
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 }