0022922: Clean up warnings on uninitialized / unused variables
[occt.git] / src / GeomFill / GeomFill_QuasiAngularConvertor.cxx
CommitLineData
7fd59977 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
53GeomFill_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
61Standard_Boolean GeomFill_QuasiAngularConvertor::Initialized() const
62{
63 return myinit;
64}
65
66void 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
111void 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);
7fd59977 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
c6541a0c 143 if ((M_PI/2 - beta)> NullAngle) {
7fd59977 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
192void 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;
c6541a0c 271 if ((M_PI/2 - beta)> NullAngle) {
7fd59977 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
344void 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;
c6541a0c 471 if ((M_PI/2 - beta)> NullAngle) {
7fd59977 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}