2 // Created: Mon Aug 28 16:32:43 1995
3 // Author: Laurent BOURESCHE
5 // Modified: 28/02/1996 by PMN : HermiteCoefficients added
6 // Modified: 18/06/1996 by PMN : NULL reference.
7 // Modified: 19/02/1997 by JCT : EvalPoly2Var added
9 #define No_Standard_RangeError
10 #define No_Standard_OutOfRange
13 #include <math_Matrix.hxx>
14 #include <math_Gauss.hxx>
15 #include <Standard_ConstructionError.hxx>
16 #include <GeomAbs_Shape.hxx>
18 // To convert points array into Real ..
19 // *********************************
21 #define Dimension_gen 2
22 #define Array1OfPoints TColgp_Array1OfPnt2d
23 #define Point gp_Pnt2d
25 #include <PLib_ChangeDim.gxx>
31 #define Dimension_gen 3
32 #define Array1OfPoints TColgp_Array1OfPnt
35 #include <PLib_ChangeDim.gxx>
41 #include <math_Gauss.hxx>
44 void PLib::InternalBinomial(const Standard_Integer N,
45 Standard_Integer& maxbinom,
46 Standard_Address& binom)
49 Standard_Integer i,im1,ip1,id2,md2,md3,j,k;
50 Standard_Integer np1 = N + 1;
51 Standard_Integer** nwbin = new Standard_Integer* [np1];
54 for (i = 0; i <= maxbinom; i++) nwbin[i] =
55 ((Standard_Integer**)binom)[i];
56 delete [] ((Standard_Integer**)binom);
60 nwbin[0] = new Standard_Integer [1];
65 for (i = maxbinom + 1; i < np1; i++) {
72 ((Standard_Integer**)binom)[i] = new Standard_Integer [ip1];
74 for (j = 0; j < id2; j++) {
75 ((Standard_Integer**)binom)[i][j] =
76 k + ((Standard_Integer**)binom)[im1][j];
77 k = ((Standard_Integer**)binom)[im1][j];
80 if (j > md2) j = im1 - j;
81 ((Standard_Integer**)binom)[i][id2] =
82 k + ((Standard_Integer**)binom)[im1][j];
84 for (j = ip1 - md3; j < ip1; j++) {
85 ((Standard_Integer**)binom)[i][j] =
86 ((Standard_Integer**)binom)[i][i - j];
93 static Standard_Integer storage_size = 0 ;
94 static Standard_Real *derivative_storage= NULL;
95 static Standard_Integer binomial_size = 0;
96 static Standard_Real *binomial_array = NULL;
98 static void LocalArray(const Standard_Integer newsize,
99 Standard_Integer& size,
102 // update a local array
103 if (newsize > size) {
104 if (*arr) delete [] *arr;
106 *arr = new Standard_Real [size];
110 //=======================================================================
111 //function : RationalDerivative
113 //=======================================================================
115 void PLib::RationalDerivative(const Standard_Integer Degree,
116 const Standard_Integer DerivativeRequest,
117 const Standard_Integer Dimension,
119 Standard_Real& RDers,
120 const Standard_Boolean All)
123 // Our purpose is to compute f = (u/v) derivated N = DerivativeRequest times
126 // Let C(N,P) be the binomial
131 // u = SUM C (q,p) f v
138 // (q) ( (q) (p) (q-p) )
139 // f = (1/v) ( u - SUM C (q,p) f v )
143 // make arrays for the binomial since computing it each time could raise a performance
145 // As oppose to the method below the <Der> array is organized in the following
148 // u (1) u (2) .... u (Dimension) v (1)
151 // u (1) u (2) .... u (Dimension) v (1)
153 // ............................................
155 // (Degree) (Degree) (Degree) (Degree)
156 // u (1) u (2) .... u (Dimension) v (1)
159 Standard_Real Inverse;
160 Standard_Real *PolesArray = &Ders;
161 Standard_Real *RationalArray = &RDers;
162 Standard_Real Factor ;
163 Standard_Integer ii, Index, OtherIndex, Index1, Index2, jj;
164 if (Dimension == 3) {
165 Standard_Integer DeRequest1 = DerivativeRequest + 1;
166 Standard_Integer MinDegRequ = DerivativeRequest;
167 if (MinDegRequ > Degree) MinDegRequ = Degree;
168 if (DeRequest1 > binomial_size) {
169 if (binomial_size > 0) {
170 delete [] binomial_array;
172 binomial_array = new Standard_Real [DeRequest1];
173 binomial_size = DeRequest1;
176 for (ii = 0 ; ii < DeRequest1 ; ii++) {
177 binomial_array[ii] = 1.0e0 ;
180 Standard_Integer DimDeRequ1 = (DeRequest1 << 1) + DeRequest1;
181 if (storage_size < DimDeRequ1) {
182 if (storage_size > 0)
183 delete [] derivative_storage ;
184 derivative_storage = new Standard_Real [DimDeRequ1];
185 storage_size = DimDeRequ1;
187 RationalArray = derivative_storage ;
190 Inverse = 1.0e0 / PolesArray[3] ;
195 for (ii = 0 ; ii <= MinDegRequ ; ii++) {
198 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
199 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
200 RationalArray[Index] = PolesArray[OtherIndex];
204 for (jj = ii - 1 ; jj >= 0 ; jj--) {
205 Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
206 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
207 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
208 RationalArray[Index] -= Factor * RationalArray[Index1];
213 for (jj = ii ; jj >= 1 ; jj--) {
214 binomial_array[jj] += binomial_array[jj - 1] ;
216 RationalArray[Index] *= Inverse ; Index++;
217 RationalArray[Index] *= Inverse ; Index++;
218 RationalArray[Index] *= Inverse ; Index++;
221 for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
224 RationalArray[Index] = 0.0e0; Index++;
225 RationalArray[Index] = 0.0e0; Index++;
226 RationalArray[Index] = 0.0e0;
229 for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
230 Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
231 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
232 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
233 RationalArray[Index] -= Factor * RationalArray[Index1];
238 for (jj = ii ; jj >= 1 ; jj--) {
239 binomial_array[jj] += binomial_array[jj - 1] ;
241 RationalArray[Index] *= Inverse; Index++;
242 RationalArray[Index] *= Inverse; Index++;
243 RationalArray[Index] *= Inverse; Index++;
247 RationalArray = &RDers ;
248 Standard_Integer DimDeRequ = (DerivativeRequest << 1) + DerivativeRequest;
249 RationalArray[0] = derivative_storage[DimDeRequ]; DimDeRequ++;
250 RationalArray[1] = derivative_storage[DimDeRequ]; DimDeRequ++;
251 RationalArray[2] = derivative_storage[DimDeRequ];
256 Standard_Integer Dimension1 = Dimension + 1;
257 Standard_Integer Dimension2 = Dimension << 1;
258 Standard_Integer DeRequest1 = DerivativeRequest + 1;
259 Standard_Integer MinDegRequ = DerivativeRequest;
260 if (MinDegRequ > Degree) MinDegRequ = Degree;
261 if (DeRequest1 > binomial_size) {
262 if (binomial_size > 0) {
263 delete [] binomial_array;
265 binomial_array = new Standard_Real [DeRequest1];
266 binomial_size = DeRequest1;
269 for (ii = 0 ; ii < DeRequest1 ; ii++) {
270 binomial_array[ii] = 1.0e0 ;
273 Standard_Integer DimDeRequ1 = Dimension * DeRequest1;
274 if (storage_size < DimDeRequ1) {
275 if (storage_size > 0)
276 delete [] derivative_storage ;
277 derivative_storage = new Standard_Real [DimDeRequ1];
278 storage_size = DimDeRequ1;
280 RationalArray = derivative_storage ;
283 Inverse = 1.0e0 / PolesArray[Dimension] ;
285 Index2 = - Dimension2;
288 for (ii = 0 ; ii <= MinDegRequ ; ii++) {
292 for (kk = 0 ; kk < Dimension ; kk++) {
293 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
298 for (jj = ii - 1 ; jj >= 0 ; jj--) {
299 Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
301 for (kk = 0 ; kk < Dimension ; kk++) {
302 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
305 Index1 -= Dimension2 ;
308 for (jj = ii ; jj >= 1 ; jj--) {
309 binomial_array[jj] += binomial_array[jj - 1] ;
312 for (kk = 0 ; kk < Dimension ; kk++) {
313 RationalArray[Index] *= Inverse ; Index++;
317 for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
321 for (kk = 0 ; kk < Dimension ; kk++) {
322 RationalArray[Index] = 0.0e0 ; Index++;
326 for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
327 Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
329 for (kk = 0 ; kk < Dimension ; kk++) {
330 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
333 Index1 -= Dimension2 ;
336 for (jj = ii ; jj >= 1 ; jj--) {
337 binomial_array[jj] += binomial_array[jj - 1] ;
340 for (kk = 0 ; kk < Dimension ; kk++) {
341 RationalArray[Index] *= Inverse; Index++;
346 RationalArray = &RDers ;
347 Standard_Integer DimDeRequ = Dimension * DerivativeRequest;
349 for (kk = 0 ; kk < Dimension ; kk++) {
350 RationalArray[kk] = derivative_storage[DimDeRequ]; DimDeRequ++;
356 //=======================================================================
357 //function : RationalDerivatives
358 //purpose : Uses Homogeneous Poles Derivatives and Deivatives of Weights
359 //=======================================================================
361 void PLib::RationalDerivatives(const Standard_Integer DerivativeRequest,
362 const Standard_Integer Dimension,
363 Standard_Real& PolesDerivates,
364 // must be an array with
365 // (DerivativeRequest + 1) * Dimension slots
366 Standard_Real& WeightsDerivates,
367 // must be an array with
368 // (DerivativeRequest + 1) slots
369 Standard_Real& RationalDerivates)
372 // Our purpose is to compute f = (u/v) derivated N times
375 // Let C(N,P) be the binomial
380 // u = SUM C (q,p) f v
387 // (q) ( (q) (p) (q-p) )
388 // f = (1/v) ( u - SUM C (q,p) f v )
392 // make arrays for the binomial since computing it each time could
393 // raize a performance issue
395 Standard_Real Inverse;
396 Standard_Real *PolesArray = &PolesDerivates;
397 Standard_Real *WeightsArray = &WeightsDerivates;
398 Standard_Real *RationalArray = &RationalDerivates;
399 Standard_Real Factor ;
401 Standard_Integer ii, Index, Index1, Index2, jj;
402 Standard_Integer DeRequest1 = DerivativeRequest + 1;
404 if (DeRequest1 > binomial_size) {
405 if (binomial_size > 0) {
406 delete [] binomial_array;
408 binomial_array = new Standard_Real [DeRequest1];
409 binomial_size = DeRequest1;
412 for (ii = 0 ; ii < DeRequest1 ; ii++) {
413 binomial_array[ii] = 1.0e0 ;
415 Inverse = 1.0e0 / WeightsArray[0] ;
416 if (Dimension == 3) {
420 for (ii = 0 ; ii < DeRequest1 ; ii++) {
423 RationalArray[Index] = PolesArray[Index] ; Index++;
424 RationalArray[Index] = PolesArray[Index] ; Index++;
425 RationalArray[Index] = PolesArray[Index] ;
428 for (jj = ii - 1 ; jj >= 0 ; jj--) {
429 Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
430 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
431 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
432 RationalArray[Index] -= Factor * RationalArray[Index1];
437 for (jj = ii ; jj >= 1 ; jj--) {
438 binomial_array[jj] += binomial_array[jj - 1] ;
440 RationalArray[Index] *= Inverse ; Index++;
441 RationalArray[Index] *= Inverse ; Index++;
442 RationalArray[Index] *= Inverse ; Index++;
447 Standard_Integer Dimension2 = Dimension << 1;
449 Index2 = - Dimension2;
451 for (ii = 0 ; ii < DeRequest1 ; ii++) {
455 for (kk = 0 ; kk < Dimension ; kk++) {
456 RationalArray[Index] = PolesArray[Index]; Index++;
460 for (jj = ii - 1 ; jj >= 0 ; jj--) {
461 Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
463 for (kk = 0 ; kk < Dimension ; kk++) {
464 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
467 Index1 -= Dimension2;
470 for (jj = ii ; jj >= 1 ; jj--) {
471 binomial_array[jj] += binomial_array[jj - 1] ;
474 for (kk = 0 ; kk < Dimension ; kk++) {
475 RationalArray[Index] *= Inverse ; Index++;
481 //=======================================================================
482 //function : This evaluates a polynomial and its derivatives
483 //purpose : up to the requested order
484 //=======================================================================
486 void PLib::EvalPolynomial(const Standard_Real Par,
487 const Standard_Integer DerivativeRequest,
488 const Standard_Integer Degree,
489 const Standard_Integer Dimension,
490 Standard_Real& PolynomialCoeff,
491 Standard_Real& Results)
493 // the polynomial coefficients are assumed to be stored as follows :
495 // [0] [Dimension -1] X coefficient
497 // [Dimension] [Dimension + Dimension -1] X coefficient
499 // [2 * Dimension] [2 * Dimension + Dimension-1] X coefficient
501 // ...................................................
505 // [d * Dimension] [d * Dimension + Dimension-1] X coefficient
507 // where d is the Degree
510 Standard_Integer DegreeDimension = Degree * Dimension;
513 Standard_Real *RA = &Results ;
514 Standard_Real *PA = &PolynomialCoeff ;
515 Standard_Real *tmpRA = RA;
516 Standard_Real *tmpPA = PA + DegreeDimension;
522 if (DerivativeRequest > 0 ) {
524 Standard_Real *valRA;
525 Standard_Integer ii, LocalRequest;
526 Standard_Integer Index1, Index2;
527 Standard_Integer MaxIndex1, MaxIndex2;
528 if (DerivativeRequest < Degree) {
529 LocalRequest = DerivativeRequest;
530 MaxIndex2 = MaxIndex1 = LocalRequest;
533 LocalRequest = Degree;
534 MaxIndex2 = MaxIndex1 = Degree;
538 for (ii = 1; ii <= LocalRequest; ii++) {
539 *tmpRA = 0.0e0; tmpRA ++ ;
542 for (jj = Degree ; jj > 0 ; jj--) {
547 for (ii = LocalRequest ; ii > 0 ; ii--) {
549 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
554 *valRA = Par * (*valRA) + (*tmpPA);
559 for (jj = Degree ; jj > 0 ; jj--) {
561 *tmpRA = Par * (*tmpRA) + (*tmpPA);
567 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
568 *tmpRA = *tmpPA; tmpRA++;
570 if (DerivativeRequest > 0 ) {
571 Standard_Real *valRA;
572 Standard_Integer ii, LocalRequest;
573 Standard_Integer Index1, Index2;
574 Standard_Integer MaxIndex1, MaxIndex2;
575 if (DerivativeRequest < Degree) {
576 LocalRequest = DerivativeRequest;
577 MaxIndex2 = MaxIndex1 = LocalRequest << 1;
580 LocalRequest = Degree;
581 MaxIndex2 = MaxIndex1 = DegreeDimension;
585 for (ii = 1; ii <= LocalRequest; ii++) {
586 *tmpRA = 0.0e0; tmpRA++;
587 *tmpRA = 0.0e0; tmpRA++;
590 for (jj = Degree ; jj > 0 ; jj--) {
596 for (ii = LocalRequest ; ii > 0 ; ii--) {
598 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
603 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
605 Index1 = MaxIndex1 + 1;
606 Index2 = MaxIndex2 + 1;
608 for (ii = LocalRequest ; ii > 0 ; ii--) {
610 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
615 *valRA = Par * (*valRA) + (*tmpPA);
622 for (jj = Degree ; jj > 0 ; jj--) {
625 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
626 *tmpRA = Par * (*tmpRA) + (*tmpPA);
633 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
634 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
635 *tmpRA = *tmpPA; tmpRA++;
637 if (DerivativeRequest > 0 ) {
638 Standard_Real *valRA;
639 Standard_Integer ii, LocalRequest;
640 Standard_Integer Index1, Index2;
641 Standard_Integer MaxIndex1, MaxIndex2;
642 if (DerivativeRequest < Degree) {
643 LocalRequest = DerivativeRequest;
644 MaxIndex2 = MaxIndex1 = (LocalRequest << 1) + LocalRequest;
647 LocalRequest = Degree;
648 MaxIndex2 = MaxIndex1 = DegreeDimension;
652 for (ii = 1; ii <= LocalRequest; ii++) {
653 *tmpRA = 0.0e0; tmpRA++;
654 *tmpRA = 0.0e0; tmpRA++;
655 *tmpRA = 0.0e0; tmpRA++;
658 for (jj = Degree ; jj > 0 ; jj--) {
664 for (ii = LocalRequest ; ii > 0 ; ii--) {
666 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
671 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
673 Index1 = MaxIndex1 + 1;
674 Index2 = MaxIndex2 + 1;
676 for (ii = LocalRequest ; ii > 0 ; ii--) {
678 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
683 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
685 Index1 = MaxIndex1 + 2;
686 Index2 = MaxIndex2 + 2;
688 for (ii = LocalRequest ; ii > 0 ; ii--) {
690 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
695 *valRA = Par * (*valRA) + (*tmpPA);
702 for (jj = Degree ; jj > 0 ; jj--) {
705 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
706 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
707 *tmpRA = Par * (*tmpRA) + (*tmpPA);
714 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
715 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
716 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
718 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
719 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
720 *tmpRA = *tmpPA; tmpRA++;
722 if (DerivativeRequest > 0 ) {
723 Standard_Real *valRA;
724 Standard_Integer ii, LocalRequest;
725 Standard_Integer Index1, Index2;
726 Standard_Integer MaxIndex1, MaxIndex2;
727 if (DerivativeRequest < Degree) {
728 LocalRequest = DerivativeRequest;
729 MaxIndex2 = MaxIndex1 = (LocalRequest << 2) + (LocalRequest << 1);
732 LocalRequest = Degree;
733 MaxIndex2 = MaxIndex1 = DegreeDimension;
737 for (ii = 1; ii <= LocalRequest; ii++) {
738 *tmpRA = 0.0e0; tmpRA++;
739 *tmpRA = 0.0e0; tmpRA++;
740 *tmpRA = 0.0e0; tmpRA++;
742 *tmpRA = 0.0e0; tmpRA++;
743 *tmpRA = 0.0e0; tmpRA++;
744 *tmpRA = 0.0e0; tmpRA++;
747 for (jj = Degree ; jj > 0 ; jj--) {
753 for (ii = LocalRequest ; ii > 0 ; ii--) {
755 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
760 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
762 Index1 = MaxIndex1 + 1;
763 Index2 = MaxIndex2 + 1;
765 for (ii = LocalRequest ; ii > 0 ; ii--) {
767 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
772 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
774 Index1 = MaxIndex1 + 2;
775 Index2 = MaxIndex2 + 2;
777 for (ii = LocalRequest ; ii > 0 ; ii--) {
779 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
784 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
786 Index1 = MaxIndex1 + 3;
787 Index2 = MaxIndex2 + 3;
789 for (ii = LocalRequest ; ii > 0 ; ii--) {
791 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
796 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
798 Index1 = MaxIndex1 + 4;
799 Index2 = MaxIndex2 + 4;
801 for (ii = LocalRequest ; ii > 0 ; ii--) {
803 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
808 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
810 Index1 = MaxIndex1 + 5;
811 Index2 = MaxIndex2 + 5;
813 for (ii = LocalRequest ; ii > 0 ; ii--) {
815 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
820 *valRA = Par * (*valRA) + (*tmpPA);
827 for (jj = Degree ; jj > 0 ; jj--) {
831 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
832 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
833 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
835 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
836 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
837 *tmpRA = Par * (*tmpRA) + (*tmpPA);
844 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
845 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
846 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
848 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
849 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
850 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
852 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
853 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
854 *tmpRA = *tmpPA; tmpRA++;
856 if (DerivativeRequest > 0 ) {
857 Standard_Real *valRA;
858 Standard_Integer ii, LocalRequest;
859 Standard_Integer Index1, Index2;
860 Standard_Integer MaxIndex1, MaxIndex2;
861 if (DerivativeRequest < Degree) {
862 LocalRequest = DerivativeRequest;
863 MaxIndex2 = MaxIndex1 = (LocalRequest << 3) + LocalRequest;
866 LocalRequest = Degree;
867 MaxIndex2 = MaxIndex1 = DegreeDimension;
871 for (ii = 1; ii <= LocalRequest; ii++) {
872 *tmpRA = 0.0e0; tmpRA++;
873 *tmpRA = 0.0e0; tmpRA++;
874 *tmpRA = 0.0e0; tmpRA++;
876 *tmpRA = 0.0e0; tmpRA++;
877 *tmpRA = 0.0e0; tmpRA++;
878 *tmpRA = 0.0e0; tmpRA++;
880 *tmpRA = 0.0e0; tmpRA++;
881 *tmpRA = 0.0e0; tmpRA++;
882 *tmpRA = 0.0e0; tmpRA++;
885 for (jj = Degree ; jj > 0 ; jj--) {
891 for (ii = LocalRequest ; ii > 0 ; ii--) {
893 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
898 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
900 Index1 = MaxIndex1 + 1;
901 Index2 = MaxIndex2 + 1;
903 for (ii = LocalRequest ; ii > 0 ; ii--) {
905 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
910 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
912 Index1 = MaxIndex1 + 2;
913 Index2 = MaxIndex2 + 2;
915 for (ii = LocalRequest ; ii > 0 ; ii--) {
917 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
922 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
924 Index1 = MaxIndex1 + 3;
925 Index2 = MaxIndex2 + 3;
927 for (ii = LocalRequest ; ii > 0 ; ii--) {
929 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
934 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
936 Index1 = MaxIndex1 + 4;
937 Index2 = MaxIndex2 + 4;
939 for (ii = LocalRequest ; ii > 0 ; ii--) {
941 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
946 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
948 Index1 = MaxIndex1 + 5;
949 Index2 = MaxIndex2 + 5;
951 for (ii = LocalRequest ; ii > 0 ; ii--) {
953 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
958 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
960 Index1 = MaxIndex1 + 6;
961 Index2 = MaxIndex2 + 6;
963 for (ii = LocalRequest ; ii > 0 ; ii--) {
965 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
970 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
972 Index1 = MaxIndex1 + 7;
973 Index2 = MaxIndex2 + 7;
975 for (ii = LocalRequest ; ii > 0 ; ii--) {
977 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
982 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
984 Index1 = MaxIndex1 + 8;
985 Index2 = MaxIndex2 + 8;
987 for (ii = LocalRequest ; ii > 0 ; ii--) {
989 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
994 *valRA = Par * (*valRA) + (*tmpPA);
1001 for (jj = Degree ; jj > 0 ; jj--) {
1005 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1006 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1007 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1009 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1010 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1011 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1013 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1014 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1015 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1022 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1023 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1024 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1026 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1027 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1028 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1030 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1031 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1032 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1034 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1035 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1036 *tmpRA = *tmpPA; tmpRA++;
1038 if (DerivativeRequest > 0 ) {
1039 Standard_Real *valRA;
1040 Standard_Integer ii, LocalRequest;
1041 Standard_Integer Index1, Index2;
1042 Standard_Integer MaxIndex1, MaxIndex2;
1043 if (DerivativeRequest < Degree) {
1044 LocalRequest = DerivativeRequest;
1045 MaxIndex2 = MaxIndex1 = (LocalRequest << 3) + (LocalRequest << 2);
1048 LocalRequest = Degree;
1049 MaxIndex2 = MaxIndex1 = DegreeDimension;
1053 for (ii = 1; ii <= LocalRequest; ii++) {
1054 *tmpRA = 0.0e0; tmpRA++;
1055 *tmpRA = 0.0e0; tmpRA++;
1056 *tmpRA = 0.0e0; tmpRA++;
1058 *tmpRA = 0.0e0; tmpRA++;
1059 *tmpRA = 0.0e0; tmpRA++;
1060 *tmpRA = 0.0e0; tmpRA++;
1062 *tmpRA = 0.0e0; tmpRA++;
1063 *tmpRA = 0.0e0; tmpRA++;
1064 *tmpRA = 0.0e0; tmpRA++;
1066 *tmpRA = 0.0e0; tmpRA++;
1067 *tmpRA = 0.0e0; tmpRA++;
1068 *tmpRA = 0.0e0; tmpRA++;
1071 for (jj = Degree ; jj > 0 ; jj--) {
1077 for (ii = LocalRequest ; ii > 0 ; ii--) {
1078 valRA = &RA[Index1];
1079 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1083 valRA = &RA[Index1];
1084 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1086 Index1 = MaxIndex1 + 1;
1087 Index2 = MaxIndex2 + 1;
1089 for (ii = LocalRequest ; ii > 0 ; ii--) {
1090 valRA = &RA[Index1];
1091 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1095 valRA = &RA[Index1];
1096 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1098 Index1 = MaxIndex1 + 2;
1099 Index2 = MaxIndex2 + 2;
1101 for (ii = LocalRequest ; ii > 0 ; ii--) {
1102 valRA = &RA[Index1];
1103 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1107 valRA = &RA[Index1];
1108 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1110 Index1 = MaxIndex1 + 3;
1111 Index2 = MaxIndex2 + 3;
1113 for (ii = LocalRequest ; ii > 0 ; ii--) {
1114 valRA = &RA[Index1];
1115 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1119 valRA = &RA[Index1];
1120 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1122 Index1 = MaxIndex1 + 4;
1123 Index2 = MaxIndex2 + 4;
1125 for (ii = LocalRequest ; ii > 0 ; ii--) {
1126 valRA = &RA[Index1];
1127 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1131 valRA = &RA[Index1];
1132 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1134 Index1 = MaxIndex1 + 5;
1135 Index2 = MaxIndex2 + 5;
1137 for (ii = LocalRequest ; ii > 0 ; ii--) {
1138 valRA = &RA[Index1];
1139 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1143 valRA = &RA[Index1];
1144 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1146 Index1 = MaxIndex1 + 6;
1147 Index2 = MaxIndex2 + 6;
1149 for (ii = LocalRequest ; ii > 0 ; ii--) {
1150 valRA = &RA[Index1];
1151 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1155 valRA = &RA[Index1];
1156 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1158 Index1 = MaxIndex1 + 7;
1159 Index2 = MaxIndex2 + 7;
1161 for (ii = LocalRequest ; ii > 0 ; ii--) {
1162 valRA = &RA[Index1];
1163 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1167 valRA = &RA[Index1];
1168 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1170 Index1 = MaxIndex1 + 8;
1171 Index2 = MaxIndex2 + 8;
1173 for (ii = LocalRequest ; ii > 0 ; ii--) {
1174 valRA = &RA[Index1];
1175 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1179 valRA = &RA[Index1];
1180 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1182 Index1 = MaxIndex1 + 9;
1183 Index2 = MaxIndex2 + 9;
1185 for (ii = LocalRequest ; ii > 0 ; ii--) {
1186 valRA = &RA[Index1];
1187 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1191 valRA = &RA[Index1];
1192 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1194 Index1 = MaxIndex1 + 10;
1195 Index2 = MaxIndex2 + 10;
1197 for (ii = LocalRequest ; ii > 0 ; ii--) {
1198 valRA = &RA[Index1];
1199 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1203 valRA = &RA[Index1];
1204 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1206 Index1 = MaxIndex1 + 11;
1207 Index2 = MaxIndex2 + 11;
1209 for (ii = LocalRequest ; ii > 0 ; ii--) {
1210 valRA = &RA[Index1];
1211 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1215 valRA = &RA[Index1];
1216 *valRA = Par * (*valRA) + (*tmpPA);
1223 for (jj = Degree ; jj > 0 ; jj--) {
1227 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1228 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1229 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1231 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1232 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1233 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1235 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1236 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1237 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1239 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1240 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1241 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1249 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1250 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1251 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1253 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1254 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1255 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1257 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1258 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1259 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1261 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1262 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1263 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1265 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1266 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1267 *tmpRA = *tmpPA; tmpRA++;
1269 if (DerivativeRequest > 0 ) {
1270 Standard_Real *valRA;
1271 Standard_Integer ii, LocalRequest;
1272 Standard_Integer Index1, Index2;
1273 Standard_Integer MaxIndex1, MaxIndex2;
1274 if (DerivativeRequest < Degree) {
1275 LocalRequest = DerivativeRequest;
1276 MaxIndex2 = MaxIndex1 = (LocalRequest << 4) - LocalRequest;
1279 LocalRequest = Degree;
1280 MaxIndex2 = MaxIndex1 = DegreeDimension;
1284 for (ii = 1; ii <= LocalRequest; ii++) {
1285 *tmpRA = 0.0e0; tmpRA++;
1286 *tmpRA = 0.0e0; tmpRA++;
1287 *tmpRA = 0.0e0; tmpRA++;
1289 *tmpRA = 0.0e0; tmpRA++;
1290 *tmpRA = 0.0e0; tmpRA++;
1291 *tmpRA = 0.0e0; tmpRA++;
1293 *tmpRA = 0.0e0; tmpRA++;
1294 *tmpRA = 0.0e0; tmpRA++;
1295 *tmpRA = 0.0e0; tmpRA++;
1297 *tmpRA = 0.0e0; tmpRA++;
1298 *tmpRA = 0.0e0; tmpRA++;
1299 *tmpRA = 0.0e0; tmpRA++;
1301 *tmpRA = 0.0e0; tmpRA++;
1302 *tmpRA = 0.0e0; tmpRA++;
1303 *tmpRA = 0.0e0; tmpRA++;
1306 for (jj = Degree ; jj > 0 ; jj--) {
1312 for (ii = LocalRequest ; ii > 0 ; ii--) {
1313 valRA = &RA[Index1];
1314 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1318 valRA = &RA[Index1];
1319 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1321 Index1 = MaxIndex1 + 1;
1322 Index2 = MaxIndex2 + 1;
1324 for (ii = LocalRequest ; ii > 0 ; ii--) {
1325 valRA = &RA[Index1];
1326 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1330 valRA = &RA[Index1];
1331 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1333 Index1 = MaxIndex1 + 2;
1334 Index2 = MaxIndex2 + 2;
1336 for (ii = LocalRequest ; ii > 0 ; ii--) {
1337 valRA = &RA[Index1];
1338 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1342 valRA = &RA[Index1];
1343 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1345 Index1 = MaxIndex1 + 3;
1346 Index2 = MaxIndex2 + 3;
1348 for (ii = LocalRequest ; ii > 0 ; ii--) {
1349 valRA = &RA[Index1];
1350 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1354 valRA = &RA[Index1];
1355 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1357 Index1 = MaxIndex1 + 4;
1358 Index2 = MaxIndex2 + 4;
1360 for (ii = LocalRequest ; ii > 0 ; ii--) {
1361 valRA = &RA[Index1];
1362 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1366 valRA = &RA[Index1];
1367 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1369 Index1 = MaxIndex1 + 5;
1370 Index2 = MaxIndex2 + 5;
1372 for (ii = LocalRequest ; ii > 0 ; ii--) {
1373 valRA = &RA[Index1];
1374 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1378 valRA = &RA[Index1];
1379 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1381 Index1 = MaxIndex1 + 6;
1382 Index2 = MaxIndex2 + 6;
1384 for (ii = LocalRequest ; ii > 0 ; ii--) {
1385 valRA = &RA[Index1];
1386 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1390 valRA = &RA[Index1];
1391 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1393 Index1 = MaxIndex1 + 7;
1394 Index2 = MaxIndex2 + 7;
1396 for (ii = LocalRequest ; ii > 0 ; ii--) {
1397 valRA = &RA[Index1];
1398 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1402 valRA = &RA[Index1];
1403 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1405 Index1 = MaxIndex1 + 8;
1406 Index2 = MaxIndex2 + 8;
1408 for (ii = LocalRequest ; ii > 0 ; ii--) {
1409 valRA = &RA[Index1];
1410 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1414 valRA = &RA[Index1];
1415 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1417 Index1 = MaxIndex1 + 9;
1418 Index2 = MaxIndex2 + 9;
1420 for (ii = LocalRequest ; ii > 0 ; ii--) {
1421 valRA = &RA[Index1];
1422 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1426 valRA = &RA[Index1];
1427 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1429 Index1 = MaxIndex1 + 10;
1430 Index2 = MaxIndex2 + 10;
1432 for (ii = LocalRequest ; ii > 0 ; ii--) {
1433 valRA = &RA[Index1];
1434 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1438 valRA = &RA[Index1];
1439 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1441 Index1 = MaxIndex1 + 11;
1442 Index2 = MaxIndex2 + 11;
1444 for (ii = LocalRequest ; ii > 0 ; ii--) {
1445 valRA = &RA[Index1];
1446 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1450 valRA = &RA[Index1];
1451 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1453 Index1 = MaxIndex1 + 12;
1454 Index2 = MaxIndex2 + 12;
1456 for (ii = LocalRequest ; ii > 0 ; ii--) {
1457 valRA = &RA[Index1];
1458 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1462 valRA = &RA[Index1];
1463 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1465 Index1 = MaxIndex1 + 13;
1466 Index2 = MaxIndex2 + 13;
1468 for (ii = LocalRequest ; ii > 0 ; ii--) {
1469 valRA = &RA[Index1];
1470 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1474 valRA = &RA[Index1];
1475 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1477 Index1 = MaxIndex1 + 14;
1478 Index2 = MaxIndex2 + 14;
1480 for (ii = LocalRequest ; ii > 0 ; ii--) {
1481 valRA = &RA[Index1];
1482 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1486 valRA = &RA[Index1];
1487 *valRA = Par * (*valRA) + (*tmpPA);
1494 for (jj = Degree ; jj > 0 ; jj--) {
1498 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1499 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1500 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1502 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1503 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1504 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1506 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1507 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1508 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1510 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1511 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1512 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1514 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1515 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1516 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1523 Standard_Integer kk ;
1524 for ( kk = 0 ; kk < Dimension ; kk++) {
1525 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1528 if (DerivativeRequest > 0 ) {
1529 Standard_Real *valRA;
1530 Standard_Integer ii, LocalRequest;
1531 Standard_Integer Index1, Index2;
1532 Standard_Integer MaxIndex1, MaxIndex2;
1533 if (DerivativeRequest < Degree) {
1534 LocalRequest = DerivativeRequest;
1535 MaxIndex2 = MaxIndex1 = LocalRequest * Dimension;
1538 LocalRequest = Degree;
1539 MaxIndex2 = MaxIndex1 = DegreeDimension;
1541 MaxIndex2 -= Dimension;
1543 for (ii = 1; ii <= MaxIndex1; ii++) {
1544 *tmpRA = 0.0e0; tmpRA++;
1547 for (jj = Degree ; jj > 0 ; jj--) {
1550 for (kk = 0 ; kk < Dimension ; kk++) {
1551 Index1 = MaxIndex1 + kk;
1552 Index2 = MaxIndex2 + kk;
1554 for (ii = LocalRequest ; ii > 0 ; ii--) {
1555 valRA = &RA[Index1];
1556 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1557 Index1 -= Dimension;
1558 Index2 -= Dimension;
1560 valRA = &RA[Index1];
1561 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1568 for (jj = Degree ; jj > 0 ; jj--) {
1572 for (kk = 0 ; kk < Dimension ; kk++) {
1573 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1582 //=======================================================================
1583 //function : This evaluates a polynomial without derivative
1585 //=======================================================================
1587 void PLib::NoDerivativeEvalPolynomial(const Standard_Real Par,
1588 const Standard_Integer Degree,
1589 const Standard_Integer Dimension,
1590 const Standard_Integer DegreeDimension,
1591 Standard_Real& PolynomialCoeff,
1592 Standard_Real& Results)
1594 Standard_Integer jj;
1595 Standard_Real *RA = &Results ;
1596 Standard_Real *PA = &PolynomialCoeff ;
1597 Standard_Real *tmpRA = RA;
1598 Standard_Real *tmpPA = PA + DegreeDimension;
1600 switch (Dimension) {
1605 for (jj = Degree ; jj > 0 ; jj--) {
1608 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1613 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1617 for (jj = Degree ; jj > 0 ; jj--) {
1621 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1622 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1628 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1629 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1633 for (jj = Degree ; jj > 0 ; jj--) {
1637 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1638 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1639 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1645 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1646 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1647 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1649 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1650 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1654 for (jj = Degree ; jj > 0 ; jj--) {
1658 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1659 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1660 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1662 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1663 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1664 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1670 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1671 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1672 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1674 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1675 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1676 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1678 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1679 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1683 for (jj = Degree ; jj > 0 ; jj--) {
1687 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1688 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1689 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1691 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1692 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1693 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1695 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1696 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1697 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1703 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1704 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1705 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1707 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1708 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1709 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1711 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1712 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1713 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1715 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1716 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1720 for (jj = Degree ; jj > 0 ; jj--) {
1724 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1725 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1726 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1728 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1729 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1730 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1732 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1733 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1734 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1736 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1737 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1738 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1744 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1745 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1746 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1748 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1749 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1750 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1752 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1753 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1754 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1756 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1757 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1758 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1760 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1761 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1765 for (jj = Degree ; jj > 0 ; jj--) {
1769 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1770 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1771 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1773 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1774 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1775 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1777 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1778 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1779 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1781 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1782 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1783 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1785 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1786 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1787 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1793 Standard_Integer kk ;
1794 for ( kk = 0 ; kk < Dimension ; kk++) {
1795 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1799 for (jj = Degree ; jj > 0 ; jj--) {
1803 for (kk = 0 ; kk < Dimension ; kk++) {
1804 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1812 //=======================================================================
1813 //function : This evaluates a polynomial of 2 variables
1814 //purpose : or its derivative at the requested orders
1815 //=======================================================================
1817 void PLib::EvalPoly2Var(const Standard_Real UParameter,
1818 const Standard_Real VParameter,
1819 const Standard_Integer UDerivativeRequest,
1820 const Standard_Integer VDerivativeRequest,
1821 const Standard_Integer UDegree,
1822 const Standard_Integer VDegree,
1823 const Standard_Integer Dimension,
1824 Standard_Real& PolynomialCoeff,
1825 Standard_Real& Results)
1827 // the polynomial coefficients are assumed to be stored as follows :
1829 // [0] [Dimension -1] U V coefficient
1831 // [Dimension] [Dimension + Dimension -1] U V coefficient
1833 // [2 * Dimension] [2 * Dimension + Dimension-1] U V coefficient
1835 // ...................................................
1839 // [m * Dimension] [m * Dimension + Dimension-1] U V coefficient
1841 // where m = UDegree
1844 // [(m+1) * Dimension] [(m+1) * Dimension + Dimension-1] U V coefficient
1846 // ...................................................
1849 // [2*m * Dimension] [2*m * Dimension + Dimension-1] U V coefficient
1851 // ...................................................
1854 // [m*n * Dimension] [m*n * Dimension + Dimension-1] U V coefficient
1856 // where n = VDegree
1858 Standard_Integer Udim = (VDegree+1)*Dimension,
1859 index = Udim*UDerivativeRequest;
1860 TColStd_Array1OfReal Curve(1, Udim*(UDerivativeRequest+1));
1861 TColStd_Array1OfReal Point(1, Dimension*(VDerivativeRequest+1));
1862 Standard_Real * Result = (Standard_Real *) &Curve.ChangeValue(1);
1863 Standard_Real * Digit = (Standard_Real *) &Point.ChangeValue(1);
1864 Standard_Real * ResultArray ;
1865 ResultArray = &Results ;
1867 PLib::EvalPolynomial(UParameter,
1874 PLib::EvalPolynomial(VParameter,
1881 index = Dimension*VDerivativeRequest;
1883 for (Standard_Integer i=0;i<Dimension;i++) {
1884 ResultArray[i] = Digit[index+i];
1889 static Standard_Integer storage_divided = 0 ;
1890 static Standard_Real *divided_differences_array = NULL;
1892 //=======================================================================
1893 //function : This evaluates the lagrange polynomial and its derivatives
1894 //purpose : up to the requested order that interpolates a series of
1895 //points of dimension <Dimension> with given assigned parameters
1896 //=======================================================================
1899 PLib::EvalLagrange(const Standard_Real Parameter,
1900 const Standard_Integer DerivativeRequest,
1901 const Standard_Integer Degree,
1902 const Standard_Integer Dimension,
1903 Standard_Real& Values,
1904 Standard_Real& Parameters,
1905 Standard_Real& Results)
1908 // the points are assumed to be stored as follows in the Values array :
1910 // [0] [Dimension -1] first point coefficients
1912 // [Dimension] [Dimension + Dimension -1] second point coefficients
1914 // [2 * Dimension] [2 * Dimension + Dimension-1] third point coefficients
1916 // ...................................................
1920 // [d * Dimension] [d * Dimension + Dimension-1] d + 1 point coefficients
1922 // where d is the Degree
1924 // The ParameterArray stores the parameter value assign to each point in
1925 // order described above, that is
1926 // [0] is assign to first point
1927 // [1] is assign to second point
1929 Standard_Integer ii, jj, kk, Index, Index1, ReturnCode=0;
1930 Standard_Integer local_request = DerivativeRequest;
1931 Standard_Real *ParameterArray;
1932 Standard_Real difference;
1933 Standard_Real *PointsArray;
1934 Standard_Real *ResultArray ;
1936 PointsArray = &Values ;
1937 ParameterArray = &Parameters ;
1938 ResultArray = &Results ;
1939 if (local_request >= Degree) {
1940 local_request = Degree ;
1942 LocalArray((Degree + 1) * Dimension,
1943 storage_divided, ÷d_differences_array) ;
1945 // Build the divided differences array
1948 for (ii = 0 ; ii < (Degree + 1) * Dimension ; ii++) {
1949 divided_differences_array[ii] = PointsArray[ii] ;
1952 for (ii = Degree ; ii >= 0 ; ii--) {
1954 for (jj = Degree ; jj > Degree - ii ; jj--) {
1955 Index = jj * Dimension ;
1956 Index1 = Index - Dimension ;
1958 for (kk = 0 ; kk < Dimension ; kk++) {
1959 divided_differences_array[Index + kk] -=
1960 divided_differences_array[Index1 + kk] ;
1963 ParameterArray[jj] - ParameterArray[jj - Degree -1 + ii] ;
1964 if (Abs(difference) < RealSmall()) {
1968 difference = 1.0e0 / difference ;
1970 for (kk = 0 ; kk < Dimension ; kk++) {
1971 divided_differences_array[Index + kk] *= difference ;
1977 // Evaluate the divided difference array polynomial which expresses as
1979 // P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
1980 // + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
1982 // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
1985 Index = Degree * Dimension ;
1987 for (kk = 0 ; kk < Dimension ; kk++) {
1988 ResultArray[kk] = divided_differences_array[Index + kk] ;
1991 for (ii = Dimension ; ii < (local_request + 1) * Dimension ; ii++) {
1992 ResultArray[ii] = 0.0e0 ;
1995 for (ii = Degree ; ii >= 1 ; ii--) {
1996 difference = Parameter - ParameterArray[ii - 1] ;
1998 for (jj = local_request ; jj > 0 ; jj--) {
1999 Index = jj * Dimension ;
2000 Index1 = Index - Dimension ;
2002 for (kk = 0 ; kk < Dimension ; kk++) {
2003 ResultArray[Index + kk] *= difference ;
2004 ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj ;
2007 Index = (ii -1) * Dimension ;
2009 for (kk = 0 ; kk < Dimension ; kk++) {
2010 ResultArray[kk] *= difference ;
2011 ResultArray[kk] += divided_differences_array[Index+kk] ;
2015 return (ReturnCode) ;
2018 //=======================================================================
2019 //function : This evaluates the hermite polynomial and its derivatives
2020 //purpose : up to the requested order that interpolates a series of
2021 //points of dimension <Dimension> with given assigned parameters
2022 //=======================================================================
2024 Standard_Integer PLib::EvalCubicHermite
2025 (const Standard_Real Parameter,
2026 const Standard_Integer DerivativeRequest,
2027 const Standard_Integer Dimension,
2028 Standard_Real& Values,
2029 Standard_Real& Derivatives,
2030 Standard_Real& theParameters,
2031 Standard_Real& Results)
2034 // the points are assumed to be stored as follows in the Values array :
2036 // [0] [Dimension -1] first point coefficients
2038 // [Dimension] [Dimension + Dimension -1] last point coefficients
2041 // the derivatives are assumed to be stored as follows in
2042 // the Derivatives array :
2044 // [0] [Dimension -1] first point coefficients
2046 // [Dimension] [Dimension + Dimension -1] last point coefficients
2048 // The ParameterArray stores the parameter value assign to each point in
2049 // order described above, that is
2050 // [0] is assign to first point
2051 // [1] is assign to last point
2053 Standard_Integer ii, jj, kk, pp, Index, Index1, Degree, ReturnCode;
2054 Standard_Integer local_request = DerivativeRequest ;
2058 Standard_Real ParametersArray[4];
2059 Standard_Real difference;
2060 Standard_Real inverse;
2061 Standard_Real *FirstLast;
2062 Standard_Real *PointsArray;
2063 Standard_Real *DerivativesArray;
2064 Standard_Real *ResultArray ;
2066 DerivativesArray = &Derivatives ;
2067 PointsArray = &Values ;
2068 FirstLast = &theParameters ;
2069 ResultArray = &Results ;
2070 if (local_request >= Degree) {
2071 local_request = Degree ;
2073 LocalArray((Degree + 1) * Dimension,
2074 storage_size, ÷d_differences_array) ;
2076 for (ii = 0, jj = 0 ; ii < 2 ; ii++, jj+= 2) {
2077 ParametersArray[jj] =
2078 ParametersArray[jj+1] = FirstLast[ii] ;
2081 // Build the divided differences array
2084 // initialise it at the stage 2 of the building algorithm
2085 // for devided differences
2087 inverse = FirstLast[1] - FirstLast[0] ;
2088 inverse = 1.0e0 / inverse ;
2090 for (ii = 0, jj = Dimension, kk = 2 * Dimension, pp = 3 * Dimension ;
2092 ii++, jj++, kk++, pp++) {
2093 divided_differences_array[ii] = PointsArray[ii] ;
2094 divided_differences_array[kk] = inverse *
2095 (PointsArray[jj] - PointsArray[ii]) ;
2096 divided_differences_array[jj] = DerivativesArray[ii] ;
2097 divided_differences_array[pp] = DerivativesArray[jj] ;
2100 for (ii = 1 ; ii <= Degree ; ii++) {
2102 for (jj = Degree ; jj >= ii+1 ; jj--) {
2103 Index = jj * Dimension ;
2104 Index1 = Index - Dimension ;
2106 for (kk = 0 ; kk < Dimension ; kk++) {
2107 divided_differences_array[Index + kk] -=
2108 divided_differences_array[Index1 + kk] ;
2111 for (kk = 0 ; kk < Dimension ; kk++) {
2112 divided_differences_array[Index + kk] *= inverse ;
2118 // Evaluate the divided difference array polynomial which expresses as
2120 // P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
2121 // + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
2123 // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
2126 Index = Degree * Dimension ;
2128 for (kk = 0 ; kk < Dimension ; kk++) {
2129 ResultArray[kk] = divided_differences_array[Index + kk] ;
2132 for (ii = Dimension ; ii < (local_request + 1) * Dimension ; ii++) {
2133 ResultArray[ii] = 0.0e0 ;
2136 for (ii = Degree ; ii >= 1 ; ii--) {
2137 difference = Parameter - ParametersArray[ii - 1] ;
2139 for (jj = local_request ; jj > 0 ; jj--) {
2140 Index = jj * Dimension ;
2141 Index1 = Index - Dimension ;
2143 for (kk = 0 ; kk < Dimension ; kk++) {
2144 ResultArray[Index + kk] *= difference ;
2145 ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj;
2148 Index = (ii -1) * Dimension ;
2150 for (kk = 0 ; kk < Dimension ; kk++) {
2151 ResultArray[kk] *= difference ;
2152 ResultArray[kk] += divided_differences_array[Index+kk] ;
2156 return (ReturnCode) ;
2159 //=======================================================================
2160 //function : HermiteCoefficients
2161 //purpose : calcul des polynomes d'Hermite
2162 //=======================================================================
2164 Standard_Boolean PLib::HermiteCoefficients(const Standard_Real FirstParameter,
2165 const Standard_Real LastParameter,
2166 const Standard_Integer FirstOrder,
2167 const Standard_Integer LastOrder,
2168 math_Matrix& MatrixCoefs)
2170 Standard_Integer NbCoeff = FirstOrder + LastOrder + 2, Ordre[2];
2171 Standard_Integer ii, jj, pp, cote, iof=0;
2172 Standard_Real Prod, TBorne = FirstParameter;
2173 math_Vector Coeff(1,NbCoeff), B(1, NbCoeff, 0.0);
2174 math_Matrix MAT(1,NbCoeff, 1,NbCoeff, 0.0);
2176 // Test de validites
2178 if ((FirstOrder < 0) || (LastOrder < 0)) return Standard_False;
2179 Standard_Real D1 = fabs(FirstParameter), D2 = fabs(LastParameter);
2180 if (D1 > 100 || D2 > 100) return Standard_False;
2182 if (D2 < 0.01) return Standard_False;
2183 if (fabs(LastParameter - FirstParameter) / D2 < 0.01) return Standard_False;
2185 // Calcul de la matrice a inverser (MAT)
2187 Ordre[0] = FirstOrder+1;
2188 Ordre[1] = LastOrder+1;
2190 for (cote=0; cote<=1; cote++) {
2193 for (pp=1; pp<=Ordre[cote]; pp++) {
2197 for (jj=pp; jj<=NbCoeff; jj++) {
2198 // tout se passe dans les 3 lignes suivantes
2199 MAT(ii, jj) = Coeff(jj) * Prod;
2200 Coeff(jj) *= jj - pp;
2204 TBorne = LastParameter;
2208 // resolution du systemes
2209 math_Gauss ResolCoeff(MAT, 1.0e-10);
2210 if (!ResolCoeff.IsDone()) return Standard_False;
2212 for (ii=1; ii<=NbCoeff; ii++) {
2214 ResolCoeff.Solve(B, Coeff);
2215 MatrixCoefs.SetRow( ii, Coeff);
2218 return Standard_True;
2221 //=======================================================================
2222 //function : CoefficientsPoles
2224 //=======================================================================
2226 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt& Coefs,
2227 const TColStd_Array1OfReal& WCoefs,
2228 TColgp_Array1OfPnt& Poles,
2229 TColStd_Array1OfReal& Weights)
2231 TColStd_Array1OfReal tempC(1,3*Coefs.Length());
2232 PLib::SetPoles(Coefs,tempC);
2233 TColStd_Array1OfReal tempP(1,3*Poles.Length());
2234 PLib::SetPoles(Coefs,tempP);
2235 PLib::CoefficientsPoles(3,tempC,WCoefs,tempP,Weights);
2236 PLib::GetPoles(tempP,Poles);
2239 //=======================================================================
2240 //function : CoefficientsPoles
2242 //=======================================================================
2244 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt2d& Coefs,
2245 const TColStd_Array1OfReal& WCoefs,
2246 TColgp_Array1OfPnt2d& Poles,
2247 TColStd_Array1OfReal& Weights)
2249 TColStd_Array1OfReal tempC(1,2*Coefs.Length());
2250 PLib::SetPoles(Coefs,tempC);
2251 TColStd_Array1OfReal tempP(1,2*Poles.Length());
2252 PLib::SetPoles(Coefs,tempP);
2253 PLib::CoefficientsPoles(2,tempC,WCoefs,tempP,Weights);
2254 PLib::GetPoles(tempP,Poles);
2257 //=======================================================================
2258 //function : CoefficientsPoles
2260 //=======================================================================
2262 void PLib::CoefficientsPoles (const TColStd_Array1OfReal& Coefs,
2263 const TColStd_Array1OfReal& WCoefs,
2264 TColStd_Array1OfReal& Poles,
2265 TColStd_Array1OfReal& Weights)
2267 PLib::CoefficientsPoles(1,Coefs,WCoefs,Poles,Weights);
2270 //=======================================================================
2271 //function : CoefficientsPoles
2273 //=======================================================================
2275 void PLib::CoefficientsPoles (const Standard_Integer dim,
2276 const TColStd_Array1OfReal& Coefs,
2277 const TColStd_Array1OfReal& WCoefs,
2278 TColStd_Array1OfReal& Poles,
2279 TColStd_Array1OfReal& Weights)
2281 Standard_Boolean rat = &WCoefs != NULL;
2282 Standard_Integer loc = Coefs.Lower();
2283 Standard_Integer lop = Poles.Lower();
2284 Standard_Integer lowc=0;
2285 Standard_Integer lowp=0;
2286 Standard_Integer upc = Coefs.Upper();
2287 Standard_Integer upp = Poles.Upper();
2288 Standard_Integer upwc=0;
2289 Standard_Integer upwp=0;
2290 Standard_Integer reflen = Coefs.Length()/dim;
2291 Standard_Integer i,j,k;
2294 lowc = WCoefs.Lower(); lowp = Weights.Lower();
2295 upwc = WCoefs.Upper(); upwp = Weights.Upper();
2298 for (i = 0; i < dim; i++){
2299 Poles (lop + i) = Coefs (loc + i);
2300 Poles (upp - i) = Coefs (upc - i);
2303 Weights (lowp) = WCoefs (lowc);
2304 Weights (upwp) = WCoefs (upwc);
2308 PLib::Binomial(reflen - 1);
2310 for (i = 2; i < reflen; i++ ) {
2311 Cnp = PLib::Bin(reflen - 1, i - 1);
2312 if (rat) Weights (lowp + i - 1) = WCoefs (lowc + i - 1) / Cnp;
2314 for(j = 0; j < dim; j++){
2315 Poles(lop + dim * (i-1) + j)= Coefs(loc + dim * (i-1) + j) / Cnp;
2319 for (i = 1; i <= reflen - 1; i++) {
2321 for (j = reflen - 1; j >= i; j--) {
2322 if (rat) Weights (lowp + j) += Weights (lowp + j -1);
2324 for(k = 0; k < dim; k++){
2325 Poles(lop + dim * j + k) += Poles(lop + dim * (j - 1) + k);
2331 for (i = 1; i <= reflen; i++) {
2333 for(j = 0; j < dim; j++){
2334 Poles(lop + dim * (i-1) + j) /= Weights(lowp + i -1);
2340 //=======================================================================
2341 //function : Trimming
2343 //=======================================================================
2345 void PLib::Trimming(const Standard_Real U1,
2346 const Standard_Real U2,
2347 TColgp_Array1OfPnt& Coefs,
2348 TColStd_Array1OfReal& WCoefs)
2350 TColStd_Array1OfReal temp(1,3*Coefs.Length());
2351 PLib::SetPoles(Coefs,temp);
2352 PLib::Trimming(U1,U2,3,temp,WCoefs);
2353 PLib::GetPoles(temp,Coefs);
2356 //=======================================================================
2357 //function : Trimming
2359 //=======================================================================
2361 void PLib::Trimming(const Standard_Real U1,
2362 const Standard_Real U2,
2363 TColgp_Array1OfPnt2d& Coefs,
2364 TColStd_Array1OfReal& WCoefs)
2366 TColStd_Array1OfReal temp(1,2*Coefs.Length());
2367 PLib::SetPoles(Coefs,temp);
2368 PLib::Trimming(U1,U2,2,temp,WCoefs);
2369 PLib::GetPoles(temp,Coefs);
2372 //=======================================================================
2373 //function : Trimming
2375 //=======================================================================
2377 void PLib::Trimming(const Standard_Real U1,
2378 const Standard_Real U2,
2379 TColStd_Array1OfReal& Coefs,
2380 TColStd_Array1OfReal& WCoefs)
2382 PLib::Trimming(U1,U2,1,Coefs,WCoefs);
2385 //=======================================================================
2386 //function : Trimming
2388 //=======================================================================
2390 void PLib::Trimming(const Standard_Real U1,
2391 const Standard_Real U2,
2392 const Standard_Integer dim,
2393 TColStd_Array1OfReal& Coefs,
2394 TColStd_Array1OfReal& WCoefs)
2398 // on fait le changement de variable v = (u-U1) / (U2-U1)
2399 // on exprime u = f(v) que l'on remplace dans l'expression polynomiale
2400 // decomposee sous la forme du schema iteratif de horner.
2402 Standard_Real lsp = U2 - U1;
2403 Standard_Integer indc, indw=0;
2404 Standard_Integer upc = Coefs.Upper() - dim + 1, upw=0;
2405 Standard_Integer len = Coefs.Length()/dim;
2406 Standard_Boolean rat = &WCoefs != NULL;
2409 if(len != WCoefs.Length())
2410 Standard_Failure::Raise("PLib::Trimming : nbcoefs/dim != nbweights !!!");
2411 upw = WCoefs.Upper();
2415 for (Standard_Integer i = 1; i <= len; i++) {
2416 Standard_Integer j ;
2417 indc = upc - dim*(i-1);
2418 if (rat) indw = upw - i + 1;
2419 //calcul du coefficient de degre le plus faible a l'iteration i
2421 for( j = 0; j < dim; j++){
2422 Coefs(indc - dim + j) += U1 * Coefs(indc + j);
2424 if (rat) WCoefs(indw - 1) += U1 * WCoefs(indw);
2426 //calcul des coefficients intermediaires :
2431 for(Standard_Integer k = 0; k < dim; k++){
2432 Coefs(indc - dim + k) =
2433 U1 * Coefs(indc + k) + lsp * Coefs(indc - dim + k);
2437 WCoefs(indw - 1) = U1 * WCoefs(indw) + lsp * WCoefs(indw - 1);
2441 //calcul du coefficient de degre le plus eleve :
2443 for(j = 0; j < dim; j++){
2444 Coefs(upc + j) *= lsp;
2446 if (rat) WCoefs(upw) *= lsp;
2450 //=======================================================================
2451 //function : CoefficientsPoles
2453 // Modified: 21/10/1996 by PMN : PolesCoefficient (PRO5852).
2454 // on ne bidouille plus les u et v c'est a l'appelant de savoir ce qu'il
2455 // fait avec BuildCache ou plus simplement d'utiliser PolesCoefficients
2456 //=======================================================================
2458 void PLib::CoefficientsPoles (const TColgp_Array2OfPnt& Coefs,
2459 const TColStd_Array2OfReal& WCoefs,
2460 TColgp_Array2OfPnt& Poles,
2461 TColStd_Array2OfReal& Weights)
2463 Standard_Boolean rat = (&WCoefs != NULL);
2464 Standard_Integer LowerRow = Poles.LowerRow();
2465 Standard_Integer UpperRow = Poles.UpperRow();
2466 Standard_Integer LowerCol = Poles.LowerCol();
2467 Standard_Integer UpperCol = Poles.UpperCol();
2468 Standard_Integer ColLength = Poles.ColLength();
2469 Standard_Integer RowLength = Poles.RowLength();
2471 // Bidouille pour retablir u et v pour les coefs calcules
2473 // Standard_Boolean inv = Standard_False; //ColLength != Coefs.ColLength();
2475 Standard_Integer Row, Col;
2476 Standard_Real W, Cnp;
2478 Standard_Integer I1, I2;
2479 Standard_Integer NPoleu , NPolev;
2481 PLib::Binomial(RowLength - 1);
2483 for (NPoleu = LowerRow; NPoleu <= UpperRow; NPoleu++){
2484 Poles (NPoleu, LowerCol) = Coefs (NPoleu, LowerCol);
2486 Weights (NPoleu, LowerCol) = WCoefs (NPoleu, LowerCol);
2489 for (Col = LowerCol + 1; Col <= UpperCol - 1; Col++) {
2490 Cnp = PLib::Bin(RowLength - 1,Col - LowerCol);
2491 Temp = Coefs (NPoleu, Col).XYZ();
2493 Poles (NPoleu, Col).SetXYZ (Temp);
2495 Weights (NPoleu, Col) = WCoefs (NPoleu, Col) / Cnp;
2498 Poles (NPoleu, UpperCol) = Coefs (NPoleu, UpperCol);
2500 Weights (NPoleu, UpperCol) = WCoefs (NPoleu, UpperCol);
2503 for (I1 = 1; I1 <= RowLength - 1; I1++) {
2505 for (I2 = UpperCol; I2 >= LowerCol + I1; I2--) {
2507 (Poles (NPoleu, I2).XYZ(), Poles (NPoleu, I2-1).XYZ());
2508 Poles (NPoleu, I2).SetXYZ (Temp);
2509 if (rat) Weights(NPoleu, I2) += Weights(NPoleu, I2-1);
2513 PLib::Binomial(ColLength - 1);
2515 for (NPolev = LowerCol; NPolev <= UpperCol; NPolev++){
2517 for (Row = LowerRow + 1; Row <= UpperRow - 1; Row++) {
2518 Cnp = PLib::Bin(ColLength - 1,Row - LowerRow);
2519 Temp = Poles (Row, NPolev).XYZ();
2521 Poles (Row, NPolev).SetXYZ (Temp);
2522 if (rat) Weights(Row, NPolev) /= Cnp;
2525 for (I1 = 1; I1 <= ColLength - 1; I1++) {
2527 for (I2 = UpperRow; I2 >= LowerRow + I1; I2--) {
2529 (Poles (I2, NPolev).XYZ(), Poles (I2-1, NPolev).XYZ());
2530 Poles (I2, NPolev).SetXYZ (Temp);
2531 if (rat) Weights(I2, NPolev) += Weights(I2-1, NPolev);
2537 for (Row = LowerRow; Row <= UpperRow; Row++) {
2539 for (Col = LowerCol; Col <= UpperCol; Col++) {
2540 W = Weights (Row, Col);
2541 Temp = Poles(Row, Col).XYZ();
2543 Poles(Row, Col).SetXYZ (Temp);
2549 //=======================================================================
2550 //function : UTrimming
2552 //=======================================================================
2554 void PLib::UTrimming(const Standard_Real U1,
2555 const Standard_Real U2,
2556 TColgp_Array2OfPnt& Coeffs,
2557 TColStd_Array2OfReal& WCoeffs)
2559 Standard_Boolean rat = &WCoeffs != NULL;
2560 Standard_Integer lr = Coeffs.LowerRow();
2561 Standard_Integer ur = Coeffs.UpperRow();
2562 Standard_Integer lc = Coeffs.LowerCol();
2563 Standard_Integer uc = Coeffs.UpperCol();
2564 TColgp_Array1OfPnt Temp (lr,ur);
2565 TColStd_Array1OfReal Temw (lr,ur);
2567 for (Standard_Integer icol = lc; icol <= uc; icol++) {
2568 Standard_Integer irow ;
2569 for ( irow = lr; irow <= ur; irow++) {
2570 Temp (irow) = Coeffs (irow, icol);
2571 if (rat) Temw (irow) = WCoeffs (irow, icol);
2573 if (rat) PLib::Trimming (U1, U2, Temp, Temw);
2574 else PLib::Trimming (U1, U2, Temp, PLib::NoWeights());
2576 for (irow = lr; irow <= ur; irow++) {
2577 Coeffs (irow, icol) = Temp (irow);
2578 if (rat) WCoeffs (irow, icol) = Temw (irow);
2583 //=======================================================================
2584 //function : VTrimming
2586 //=======================================================================
2588 void PLib::VTrimming(const Standard_Real V1,
2589 const Standard_Real V2,
2590 TColgp_Array2OfPnt& Coeffs,
2591 TColStd_Array2OfReal& WCoeffs)
2593 Standard_Boolean rat = &WCoeffs != NULL;
2594 Standard_Integer lr = Coeffs.LowerRow();
2595 Standard_Integer ur = Coeffs.UpperRow();
2596 Standard_Integer lc = Coeffs.LowerCol();
2597 Standard_Integer uc = Coeffs.UpperCol();
2598 TColgp_Array1OfPnt Temp (lc,uc);
2599 TColStd_Array1OfReal Temw (lc,uc);
2601 for (Standard_Integer irow = lr; irow <= ur; irow++) {
2602 Standard_Integer icol ;
2603 for ( icol = lc; icol <= uc; icol++) {
2604 Temp (icol) = Coeffs (irow, icol);
2605 if (rat) Temw (icol) = WCoeffs (irow, icol);
2607 if (rat) PLib::Trimming (V1, V2, Temp, Temw);
2608 else PLib::Trimming (V1, V2, Temp, PLib::NoWeights());
2610 for (icol = lc; icol <= uc; icol++) {
2611 Coeffs (irow, icol) = Temp (icol);
2612 if (rat) WCoeffs (irow, icol) = Temw (icol);
2617 //=======================================================================
2618 //function : HermiteInterpolate
2620 //=======================================================================
2622 Standard_Boolean PLib::HermiteInterpolate
2623 (const Standard_Integer Dimension,
2624 const Standard_Real FirstParameter,
2625 const Standard_Real LastParameter,
2626 const Standard_Integer FirstOrder,
2627 const Standard_Integer LastOrder,
2628 const TColStd_Array2OfReal& FirstConstr,
2629 const TColStd_Array2OfReal& LastConstr,
2630 TColStd_Array1OfReal& Coefficients)
2632 Standard_Real Pattern[3][6];
2634 // portage HP : il faut les initialiser 1 par 1
2655 math_Matrix A(0,FirstOrder+LastOrder+1, 0,FirstOrder+LastOrder+1);
2656 // The initialisation of the matrix A
2657 Standard_Integer irow ;
2658 for ( irow=0; irow<=FirstOrder; irow++) {
2659 Standard_Real FirstVal = 1.;
2661 for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
2662 A(irow,icol) = Pattern[irow][icol]*FirstVal;
2663 if (irow <= icol) FirstVal *= FirstParameter;
2667 for (irow=0; irow<=LastOrder; irow++) {
2668 Standard_Real LastVal = 1.;
2670 for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
2671 A(irow+FirstOrder+1,icol) = Pattern[irow][icol]*LastVal;
2672 if (irow <= icol) LastVal *= LastParameter;
2676 // The filled matrix A for FirstOrder=LastOrder=2 is:
2678 // 1 FP FP**2 FP**3 FP**4 FP**5
2679 // 0 1 2*FP 3*FP**2 4*FP**3 5*FP**4 FP - FirstParameter
2680 // 0 0 2 6*FP 12*FP**2 20*FP**3
2681 // 1 LP LP**2 LP**3 LP**4 LP**5
2682 // 0 1 2*LP 3*LP**2 4*LP**3 5*LP**4 LP - LastParameter
2683 // 0 0 2 6*LP 12*LP**2 20*LP**3
2685 // If FirstOrder or LastOrder <=2 then some rows and columns are missing.
2687 // If FirstOrder=1 then 3th row and 6th column are missing
2688 // If FirstOrder=LastOrder=0 then 2,3,5,6th rows and 3,4,5,6th columns are missing
2690 math_Gauss Equations(A);
2691 // cout << "A=" << A << endl;
2693 for (Standard_Integer idim=1; idim<=Dimension; idim++) {
2694 // cout << "idim=" << idim << endl;
2696 math_Vector B(0,FirstOrder+LastOrder+1);
2697 Standard_Integer icol ;
2698 for ( icol=0; icol<=FirstOrder; icol++)
2699 B(icol) = FirstConstr(idim,icol);
2701 for (icol=0; icol<=LastOrder; icol++)
2702 B(FirstOrder+1+icol) = LastConstr(idim,icol);
2703 // cout << "B=" << B << endl;
2705 // The solving of equations system A * X = B. Then B = X
2707 // cout << "After Solving" << endl << "B=" << B << endl;
2709 if (Equations.IsDone()==Standard_False) return Standard_False;
2711 // the filling of the Coefficients
2713 for (icol=0; icol<=FirstOrder+LastOrder+1; icol++)
2714 Coefficients(Dimension*icol+idim-1) = B(icol);
2716 return Standard_True;
2719 //=======================================================================
2720 //function : JacobiParameters
2722 //=======================================================================
2724 void PLib::JacobiParameters(const GeomAbs_Shape ConstraintOrder,
2725 const Standard_Integer MaxDegree,
2726 const Standard_Integer Code,
2727 Standard_Integer& NbGaussPoints,
2728 Standard_Integer& WorkDegree)
2730 // ConstraintOrder: Ordre de contrainte aux extremites :
2731 // C0 = contraintes de passage aux bornes;
2732 // C1 = C0 + contraintes de derivees 1eres;
2733 // C2 = C1 + contraintes de derivees 2ndes.
2734 // MaxDegree: Nombre maxi de coeff de la "courbe" polynomiale
2735 // d' approximation (doit etre superieur ou egal a
2736 // 2*NivConstr+2 et inferieur ou egal a 50).
2737 // Code: Code d' init. des parametres de discretisation.
2738 // (choix de NBPNTS et de NDGJAC de MAPF1C,MAPFXC).
2739 // = -5 Calcul tres rapide mais peu precis (8pts)
2740 // = -4 ' ' ' ' ' ' (10pts)
2741 // = -3 ' ' ' ' ' ' (15pts)
2742 // = -2 ' ' ' ' ' ' (20pts)
2743 // = -1 ' ' ' ' ' ' (25pts)
2744 // = 1 calcul rapide avec precision moyenne (30pts).
2745 // = 2 calcul rapide avec meilleure precision (40pts).
2746 // = 3 calcul un peu plus lent avec bonne precision (50 pts).
2747 // = 4 calcul lent avec la meilleure precision possible
2750 // The possible values of NbGaussPoints
2752 const Standard_Integer NDEG8=8, NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
2753 NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
2755 Standard_Integer NivConstr=0;
2756 switch (ConstraintOrder) {
2757 case GeomAbs_C0: NivConstr = 0; break;
2758 case GeomAbs_C1: NivConstr = 1; break;
2759 case GeomAbs_C2: NivConstr = 2; break;
2761 Standard_ConstructionError::Raise("Invalid ConstraintOrder");
2763 if (MaxDegree < 2*NivConstr+1)
2764 Standard_ConstructionError::Raise("Invalid MaxDegree");
2767 WorkDegree = MaxDegree + 9;
2769 WorkDegree = MaxDegree + 6;
2771 //---> Nbre mini de points necessaires.
2772 Standard_Integer IPMIN=0;
2773 if (WorkDegree < NDEG8)
2775 else if (WorkDegree < NDEG10)
2777 else if (WorkDegree < NDEG15)
2779 else if (WorkDegree < NDEG20)
2781 else if (WorkDegree < NDEG25)
2783 else if (WorkDegree < NDEG30)
2785 else if (WorkDegree < NDEG40)
2787 else if (WorkDegree < NDEG50)
2789 else if (WorkDegree < NDEG61)
2792 Standard_ConstructionError::Raise("Invalid MaxDegree");
2793 // ---> Nbre de points voulus.
2794 Standard_Integer IWANT=0;
2796 case -5: IWANT=NDEG8; break;
2797 case -4: IWANT=NDEG10; break;
2798 case -3: IWANT=NDEG15; break;
2799 case -2: IWANT=NDEG20; break;
2800 case -1: IWANT=NDEG25; break;
2801 case 1: IWANT=NDEG30; break;
2802 case 2: IWANT=NDEG40; break;
2803 case 3: IWANT=NDEG50; break;
2804 case 4: IWANT=NDEG61; break;
2806 Standard_ConstructionError::Raise("Invalid Code");
2808 //--> NbGaussPoints est le nombre de points de discretisation de la fonction,
2809 // il ne peut prendre que les valeurs 8,10,15,20,25,30,40,50 ou 61.
2810 // NbGaussPoints doit etre superieur strictement a WorkDegree.
2811 NbGaussPoints = Max(IPMIN,IWANT);
2812 // NbGaussPoints +=2;
2815 //=======================================================================
2816 //function : NivConstr
2817 //purpose : translates from GeomAbs_Shape to Integer
2818 //=======================================================================
2820 Standard_Integer PLib::NivConstr(const GeomAbs_Shape ConstraintOrder)
2822 Standard_Integer NivConstr=0;
2823 switch (ConstraintOrder) {
2824 case GeomAbs_C0: NivConstr = 0; break;
2825 case GeomAbs_C1: NivConstr = 1; break;
2826 case GeomAbs_C2: NivConstr = 2; break;
2828 Standard_ConstructionError::Raise("Invalid ConstraintOrder");
2833 //=======================================================================
2834 //function : ConstraintOrder
2835 //purpose : translates from Integer to GeomAbs_Shape
2836 //=======================================================================
2838 GeomAbs_Shape PLib::ConstraintOrder(const Standard_Integer NivConstr)
2840 GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
2841 switch (NivConstr) {
2842 case 0: ConstraintOrder = GeomAbs_C0; break;
2843 case 1: ConstraintOrder = GeomAbs_C1; break;
2844 case 2: ConstraintOrder = GeomAbs_C2; break;
2846 Standard_ConstructionError::Raise("Invalid NivConstr");
2848 return ConstraintOrder;
2851 //=======================================================================
2852 //function : EvalLength
2854 //=======================================================================
2856 void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
2857 Standard_Real& PolynomialCoeff,
2858 const Standard_Real U1, const Standard_Real U2,
2859 Standard_Real& Length)
2861 Standard_Integer i,j,idim, degdim;
2862 Standard_Real C1,C2,Sum,Tran,X1,X2,Der1,Der2,D1,D2,DD;
2864 Standard_Real *PolynomialArray = &PolynomialCoeff ;
2866 Standard_Integer NbGaussPoints = 4 * Min((Degree/4)+1,10);
2868 math_Vector GaussPoints(1,NbGaussPoints);
2869 math::GaussPoints(NbGaussPoints,GaussPoints);
2871 math_Vector GaussWeights(1,NbGaussPoints);
2872 math::GaussWeights(NbGaussPoints,GaussWeights);
2874 C1 = (U2 + U1) / 2.;
2875 C2 = (U2 - U1) / 2.;
2877 //-----------------------------------------------------------
2878 //****** Integration - Boucle sur les intervalles de GAUSS **
2879 //-----------------------------------------------------------
2883 for (j=1; j<=NbGaussPoints/2; j++) {
2884 // Integration en tenant compte de la symetrie
2885 Tran = C2 * GaussPoints(j);
2889 //****** Derivation sur la dimension de l'espace **
2891 degdim = Degree*Dimension;
2893 for (idim=0; idim<Dimension; idim++) {
2894 D1 = D2 = Degree * PolynomialArray [idim + degdim];
2895 for (i=Degree-1; i>=1; i--) {
2896 DD = i * PolynomialArray [idim + i*Dimension];
2904 //****** Integration **
2906 Sum += GaussWeights(j) * C2 * (Sqrt(Der1) + Sqrt(Der2));
2908 //****** Fin de boucle dur les intervalles de GAUSS **
2914 //=======================================================================
2915 //function : EvalLength
2917 //=======================================================================
2919 void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
2920 Standard_Real& PolynomialCoeff,
2921 const Standard_Real U1, const Standard_Real U2,
2922 const Standard_Real Tol,
2923 Standard_Real& Length, Standard_Real& Error)
2926 Standard_Integer NbSubInt = 1, // Current number of subintervals
2927 MaxNbIter = 13, // Max number of iterations
2928 NbIter = 1; // Current number of iterations
2929 Standard_Real dU,OldLen,LenI;
2931 PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1,U2, Length);
2937 dU = (U2-U1)/NbSubInt;
2938 for (i=1; i<=NbSubInt; i++) {
2939 PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1+(i-1)*dU,U1+i*dU, LenI);
2943 Error = Abs(OldLen-Length);
2945 while (Error > Tol && NbIter <= MaxNbIter);