| 1 | // Created on: 1995-08-28 |
| 2 | // Created by: Laurent BOURESCHE |
| 3 | // Copyright (c) 1995-1999 Matra Datavision |
| 4 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
| 5 | // |
| 6 | // This file is part of Open CASCADE Technology software library. |
| 7 | // |
| 8 | // This library is free software; you can redistribute it and / or modify it |
| 9 | // under the terms of the GNU Lesser General Public version 2.1 as published |
| 10 | // by the Free Software Foundation, with special exception defined in the file |
| 11 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
| 12 | // distribution for complete text of the license and disclaimer of any warranty. |
| 13 | // |
| 14 | // Alternatively, this file may be used under the terms of Open CASCADE |
| 15 | // commercial license or contractual agreement. |
| 16 | |
| 17 | // Modified: 28/02/1996 by PMN : HermiteCoefficients added |
| 18 | // Modified: 18/06/1996 by PMN : NULL reference. |
| 19 | // Modified: 19/02/1997 by JCT : EvalPoly2Var added |
| 20 | |
| 21 | #include <PLib.ixx> |
| 22 | #include <NCollection_LocalArray.hxx> |
| 23 | #include <math_Matrix.hxx> |
| 24 | #include <math_Gauss.hxx> |
| 25 | #include <Standard_ConstructionError.hxx> |
| 26 | #include <GeomAbs_Shape.hxx> |
| 27 | |
| 28 | #include <math_Gauss.hxx> |
| 29 | #include <math.hxx> |
| 30 | |
| 31 | // To convert points array into Real .. |
| 32 | // ********************************* |
| 33 | |
| 34 | //======================================================================= |
| 35 | //function : SetPoles |
| 36 | //purpose : |
| 37 | //======================================================================= |
| 38 | |
| 39 | void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles, |
| 40 | TColStd_Array1OfReal& FP) |
| 41 | { |
| 42 | Standard_Integer j = FP .Lower(); |
| 43 | Standard_Integer PLower = Poles.Lower(); |
| 44 | Standard_Integer PUpper = Poles.Upper(); |
| 45 | |
| 46 | for (Standard_Integer i = PLower; i <= PUpper; i++) { |
| 47 | const gp_Pnt2d& P = Poles(i); |
| 48 | FP(j) = P.Coord(1); j++; |
| 49 | FP(j) = P.Coord(2); j++; |
| 50 | } |
| 51 | } |
| 52 | |
| 53 | //======================================================================= |
| 54 | //function : SetPoles |
| 55 | //purpose : |
| 56 | //======================================================================= |
| 57 | |
| 58 | void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles, |
| 59 | const TColStd_Array1OfReal& Weights, |
| 60 | TColStd_Array1OfReal& FP) |
| 61 | { |
| 62 | Standard_Integer j = FP .Lower(); |
| 63 | Standard_Integer PLower = Poles.Lower(); |
| 64 | Standard_Integer PUpper = Poles.Upper(); |
| 65 | |
| 66 | for (Standard_Integer i = PLower; i <= PUpper; i++) { |
| 67 | Standard_Real w = Weights(i); |
| 68 | const gp_Pnt2d& P = Poles(i); |
| 69 | FP(j) = P.Coord(1) * w; j++; |
| 70 | FP(j) = P.Coord(2) * w; j++; |
| 71 | FP(j) = w; j++; |
| 72 | } |
| 73 | } |
| 74 | |
| 75 | //======================================================================= |
| 76 | //function : GetPoles |
| 77 | //purpose : |
| 78 | //======================================================================= |
| 79 | |
| 80 | void PLib::GetPoles(const TColStd_Array1OfReal& FP, |
| 81 | TColgp_Array1OfPnt2d& Poles) |
| 82 | { |
| 83 | Standard_Integer j = FP .Lower(); |
| 84 | Standard_Integer PLower = Poles.Lower(); |
| 85 | Standard_Integer PUpper = Poles.Upper(); |
| 86 | |
| 87 | for (Standard_Integer i = PLower; i <= PUpper; i++) { |
| 88 | gp_Pnt2d& P = Poles(i); |
| 89 | P.SetCoord(1,FP(j)); j++; |
| 90 | P.SetCoord(2,FP(j)); j++; |
| 91 | } |
| 92 | } |
| 93 | |
| 94 | //======================================================================= |
| 95 | //function : GetPoles |
| 96 | //purpose : |
| 97 | //======================================================================= |
| 98 | |
| 99 | void PLib::GetPoles(const TColStd_Array1OfReal& FP, |
| 100 | TColgp_Array1OfPnt2d& Poles, |
| 101 | TColStd_Array1OfReal& Weights) |
| 102 | { |
| 103 | Standard_Integer j = FP .Lower(); |
| 104 | Standard_Integer PLower = Poles.Lower(); |
| 105 | Standard_Integer PUpper = Poles.Upper(); |
| 106 | |
| 107 | for (Standard_Integer i = PLower; i <= PUpper; i++) { |
| 108 | Standard_Real w = FP(j + 2); |
| 109 | Weights(i) = w; |
| 110 | gp_Pnt2d& P = Poles(i); |
| 111 | P.SetCoord(1,FP(j) / w); j++; |
| 112 | P.SetCoord(2,FP(j) / w); j++; |
| 113 | j++; |
| 114 | } |
| 115 | } |
| 116 | |
| 117 | //======================================================================= |
| 118 | //function : SetPoles |
| 119 | //purpose : |
| 120 | //======================================================================= |
| 121 | |
| 122 | void PLib::SetPoles(const TColgp_Array1OfPnt& Poles, |
| 123 | TColStd_Array1OfReal& FP) |
| 124 | { |
| 125 | Standard_Integer j = FP .Lower(); |
| 126 | Standard_Integer PLower = Poles.Lower(); |
| 127 | Standard_Integer PUpper = Poles.Upper(); |
| 128 | |
| 129 | for (Standard_Integer i = PLower; i <= PUpper; i++) { |
| 130 | const gp_Pnt& P = Poles(i); |
| 131 | FP(j) = P.Coord(1); j++; |
| 132 | FP(j) = P.Coord(2); j++; |
| 133 | FP(j) = P.Coord(3); j++; |
| 134 | } |
| 135 | } |
| 136 | |
| 137 | //======================================================================= |
| 138 | //function : SetPoles |
| 139 | //purpose : |
| 140 | //======================================================================= |
| 141 | |
| 142 | void PLib::SetPoles(const TColgp_Array1OfPnt& Poles, |
| 143 | const TColStd_Array1OfReal& Weights, |
| 144 | TColStd_Array1OfReal& FP) |
| 145 | { |
| 146 | Standard_Integer j = FP .Lower(); |
| 147 | Standard_Integer PLower = Poles.Lower(); |
| 148 | Standard_Integer PUpper = Poles.Upper(); |
| 149 | |
| 150 | for (Standard_Integer i = PLower; i <= PUpper; i++) { |
| 151 | Standard_Real w = Weights(i); |
| 152 | const gp_Pnt& P = Poles(i); |
| 153 | FP(j) = P.Coord(1) * w; j++; |
| 154 | FP(j) = P.Coord(2) * w; j++; |
| 155 | FP(j) = P.Coord(3) * w; j++; |
| 156 | FP(j) = w; j++; |
| 157 | } |
| 158 | } |
| 159 | |
| 160 | //======================================================================= |
| 161 | //function : GetPoles |
| 162 | //purpose : |
| 163 | //======================================================================= |
| 164 | |
| 165 | void PLib::GetPoles(const TColStd_Array1OfReal& FP, |
| 166 | TColgp_Array1OfPnt& Poles) |
| 167 | { |
| 168 | Standard_Integer j = FP .Lower(); |
| 169 | Standard_Integer PLower = Poles.Lower(); |
| 170 | Standard_Integer PUpper = Poles.Upper(); |
| 171 | |
| 172 | for (Standard_Integer i = PLower; i <= PUpper; i++) { |
| 173 | gp_Pnt& P = Poles(i); |
| 174 | P.SetCoord(1,FP(j)); j++; |
| 175 | P.SetCoord(2,FP(j)); j++; |
| 176 | P.SetCoord(3,FP(j)); j++; |
| 177 | } |
| 178 | } |
| 179 | |
| 180 | //======================================================================= |
| 181 | //function : GetPoles |
| 182 | //purpose : |
| 183 | //======================================================================= |
| 184 | |
| 185 | void PLib::GetPoles(const TColStd_Array1OfReal& FP, |
| 186 | TColgp_Array1OfPnt& Poles, |
| 187 | TColStd_Array1OfReal& Weights) |
| 188 | { |
| 189 | Standard_Integer j = FP .Lower(); |
| 190 | Standard_Integer PLower = Poles.Lower(); |
| 191 | Standard_Integer PUpper = Poles.Upper(); |
| 192 | |
| 193 | for (Standard_Integer i = PLower; i <= PUpper; i++) { |
| 194 | Standard_Real w = FP(j + 3); |
| 195 | Weights(i) = w; |
| 196 | gp_Pnt& P = Poles(i); |
| 197 | P.SetCoord(1,FP(j) / w); j++; |
| 198 | P.SetCoord(2,FP(j) / w); j++; |
| 199 | P.SetCoord(3,FP(j) / w); j++; |
| 200 | j++; |
| 201 | } |
| 202 | } |
| 203 | |
| 204 | // specialized allocator |
| 205 | namespace |
| 206 | { |
| 207 | |
| 208 | class BinomAllocator |
| 209 | { |
| 210 | public: |
| 211 | |
| 212 | //! Main constructor |
| 213 | BinomAllocator (const Standard_Integer theMaxBinom) |
| 214 | : myBinom (NULL), |
| 215 | myMaxBinom (theMaxBinom) |
| 216 | { |
| 217 | Standard_Integer i, im1, ip1, id2, md2, md3, j, k; |
| 218 | Standard_Integer np1 = myMaxBinom + 1; |
| 219 | myBinom = new Standard_Integer*[np1]; |
| 220 | myBinom[0] = new Standard_Integer[1]; |
| 221 | myBinom[0][0] = 1; |
| 222 | for (i = 1; i < np1; ++i) |
| 223 | { |
| 224 | im1 = i - 1; |
| 225 | ip1 = i + 1; |
| 226 | id2 = i >> 1; |
| 227 | md2 = im1 >> 1; |
| 228 | md3 = ip1 >> 1; |
| 229 | k = 0; |
| 230 | myBinom[i] = new Standard_Integer[ip1]; |
| 231 | |
| 232 | for (j = 0; j < id2; ++j) |
| 233 | { |
| 234 | myBinom[i][j] = k + myBinom[im1][j]; |
| 235 | k = myBinom[im1][j]; |
| 236 | } |
| 237 | j = id2; |
| 238 | if (j > md2) j = im1 - j; |
| 239 | myBinom[i][id2] = k + myBinom[im1][j]; |
| 240 | |
| 241 | for (j = ip1 - md3; j < ip1; j++) |
| 242 | { |
| 243 | myBinom[i][j] = myBinom[i][i - j]; |
| 244 | } |
| 245 | } |
| 246 | } |
| 247 | |
| 248 | //! Destructor |
| 249 | ~BinomAllocator() |
| 250 | { |
| 251 | // free memory |
| 252 | for (Standard_Integer i = 0; i <= myMaxBinom; ++i) |
| 253 | { |
| 254 | delete[] myBinom[i]; |
| 255 | } |
| 256 | delete[] myBinom; |
| 257 | } |
| 258 | |
| 259 | Standard_Real Value (const Standard_Integer N, |
| 260 | const Standard_Integer P) const |
| 261 | { |
| 262 | Standard_OutOfRange_Raise_if (N > myMaxBinom, |
| 263 | "PLib, BinomAllocator: requested degree is greater than maximum supported"); |
| 264 | return Standard_Real (myBinom[N][P]); |
| 265 | } |
| 266 | |
| 267 | private: |
| 268 | Standard_Integer** myBinom; |
| 269 | Standard_Integer myMaxBinom; |
| 270 | |
| 271 | }; |
| 272 | |
| 273 | // we do not call BSplCLib here to avoid Cyclic dependency detection by WOK |
| 274 | //static BinomAllocator THE_BINOM (BSplCLib::MaxDegree() + 1); |
| 275 | static BinomAllocator THE_BINOM (25 + 1); |
| 276 | } |
| 277 | |
| 278 | //======================================================================= |
| 279 | //function : Bin |
| 280 | //purpose : |
| 281 | //======================================================================= |
| 282 | |
| 283 | Standard_Real PLib::Bin(const Standard_Integer N, |
| 284 | const Standard_Integer P) |
| 285 | { |
| 286 | return THE_BINOM.Value (N, P); |
| 287 | } |
| 288 | |
| 289 | //======================================================================= |
| 290 | //function : RationalDerivative |
| 291 | //purpose : |
| 292 | //======================================================================= |
| 293 | |
| 294 | void PLib::RationalDerivative(const Standard_Integer Degree, |
| 295 | const Standard_Integer DerivativeRequest, |
| 296 | const Standard_Integer Dimension, |
| 297 | Standard_Real& Ders, |
| 298 | Standard_Real& RDers, |
| 299 | const Standard_Boolean All) |
| 300 | { |
| 301 | // |
| 302 | // Our purpose is to compute f = (u/v) derivated N = DerivativeRequest times |
| 303 | // |
| 304 | // We Write u = fv |
| 305 | // Let C(N,P) be the binomial |
| 306 | // |
| 307 | // then we have |
| 308 | // |
| 309 | // (q) (p) (q-p) |
| 310 | // u = SUM C (q,p) f v |
| 311 | // p = 0 to q |
| 312 | // |
| 313 | // |
| 314 | // Therefore |
| 315 | // |
| 316 | // |
| 317 | // (q) ( (q) (p) (q-p) ) |
| 318 | // f = (1/v) ( u - SUM C (q,p) f v ) |
| 319 | // ( p = 0 to q-1 ) |
| 320 | // |
| 321 | // |
| 322 | // make arrays for the binomial since computing it each time could raise a performance |
| 323 | // issue |
| 324 | // As oppose to the method below the <Der> array is organized in the following |
| 325 | // fashion : |
| 326 | // |
| 327 | // u (1) u (2) .... u (Dimension) v (1) |
| 328 | // |
| 329 | // (1) (1) (1) (1) |
| 330 | // u (1) u (2) .... u (Dimension) v (1) |
| 331 | // |
| 332 | // ............................................ |
| 333 | // |
| 334 | // (Degree) (Degree) (Degree) (Degree) |
| 335 | // u (1) u (2) .... u (Dimension) v (1) |
| 336 | // |
| 337 | // |
| 338 | Standard_Real Inverse; |
| 339 | Standard_Real *PolesArray = &Ders; |
| 340 | Standard_Real *RationalArray = &RDers; |
| 341 | Standard_Real Factor ; |
| 342 | Standard_Integer ii, Index, OtherIndex, Index1, Index2, jj; |
| 343 | NCollection_LocalArray<Standard_Real> binomial_array; |
| 344 | NCollection_LocalArray<Standard_Real> derivative_storage; |
| 345 | if (Dimension == 3) { |
| 346 | Standard_Integer DeRequest1 = DerivativeRequest + 1; |
| 347 | Standard_Integer MinDegRequ = DerivativeRequest; |
| 348 | if (MinDegRequ > Degree) MinDegRequ = Degree; |
| 349 | binomial_array.Allocate (DeRequest1); |
| 350 | |
| 351 | for (ii = 0 ; ii < DeRequest1 ; ii++) { |
| 352 | binomial_array[ii] = 1.0e0 ; |
| 353 | } |
| 354 | if (!All) { |
| 355 | Standard_Integer DimDeRequ1 = (DeRequest1 << 1) + DeRequest1; |
| 356 | derivative_storage.Allocate (DimDeRequ1); |
| 357 | RationalArray = derivative_storage ; |
| 358 | } |
| 359 | |
| 360 | Inverse = 1.0e0 / PolesArray[3] ; |
| 361 | Index = 0 ; |
| 362 | Index2 = - 6; |
| 363 | OtherIndex = 0 ; |
| 364 | |
| 365 | for (ii = 0 ; ii <= MinDegRequ ; ii++) { |
| 366 | Index2 += 3; |
| 367 | Index1 = Index2; |
| 368 | RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++; |
| 369 | RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++; |
| 370 | RationalArray[Index] = PolesArray[OtherIndex]; |
| 371 | Index -= 2; |
| 372 | OtherIndex += 2; |
| 373 | |
| 374 | for (jj = ii - 1 ; jj >= 0 ; jj--) { |
| 375 | Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3]; |
| 376 | RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++; |
| 377 | RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++; |
| 378 | RationalArray[Index] -= Factor * RationalArray[Index1]; |
| 379 | Index -= 2; |
| 380 | Index1 -= 5; |
| 381 | } |
| 382 | |
| 383 | for (jj = ii ; jj >= 1 ; jj--) { |
| 384 | binomial_array[jj] += binomial_array[jj - 1] ; |
| 385 | } |
| 386 | RationalArray[Index] *= Inverse ; Index++; |
| 387 | RationalArray[Index] *= Inverse ; Index++; |
| 388 | RationalArray[Index] *= Inverse ; Index++; |
| 389 | } |
| 390 | |
| 391 | for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){ |
| 392 | Index2 += 3; |
| 393 | Index1 = Index2; |
| 394 | RationalArray[Index] = 0.0e0; Index++; |
| 395 | RationalArray[Index] = 0.0e0; Index++; |
| 396 | RationalArray[Index] = 0.0e0; |
| 397 | Index -= 2; |
| 398 | |
| 399 | for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) { |
| 400 | Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3]; |
| 401 | RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++; |
| 402 | RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++; |
| 403 | RationalArray[Index] -= Factor * RationalArray[Index1]; |
| 404 | Index -= 2; |
| 405 | Index1 -= 5; |
| 406 | } |
| 407 | |
| 408 | for (jj = ii ; jj >= 1 ; jj--) { |
| 409 | binomial_array[jj] += binomial_array[jj - 1] ; |
| 410 | } |
| 411 | RationalArray[Index] *= Inverse; Index++; |
| 412 | RationalArray[Index] *= Inverse; Index++; |
| 413 | RationalArray[Index] *= Inverse; Index++; |
| 414 | } |
| 415 | |
| 416 | if (!All) { |
| 417 | RationalArray = &RDers ; |
| 418 | Standard_Integer DimDeRequ = (DerivativeRequest << 1) + DerivativeRequest; |
| 419 | RationalArray[0] = derivative_storage[DimDeRequ]; DimDeRequ++; |
| 420 | RationalArray[1] = derivative_storage[DimDeRequ]; DimDeRequ++; |
| 421 | RationalArray[2] = derivative_storage[DimDeRequ]; |
| 422 | } |
| 423 | } |
| 424 | else { |
| 425 | Standard_Integer kk; |
| 426 | Standard_Integer Dimension1 = Dimension + 1; |
| 427 | Standard_Integer Dimension2 = Dimension << 1; |
| 428 | Standard_Integer DeRequest1 = DerivativeRequest + 1; |
| 429 | Standard_Integer MinDegRequ = DerivativeRequest; |
| 430 | if (MinDegRequ > Degree) MinDegRequ = Degree; |
| 431 | binomial_array.Allocate (DeRequest1); |
| 432 | |
| 433 | for (ii = 0 ; ii < DeRequest1 ; ii++) { |
| 434 | binomial_array[ii] = 1.0e0 ; |
| 435 | } |
| 436 | if (!All) { |
| 437 | Standard_Integer DimDeRequ1 = Dimension * DeRequest1; |
| 438 | derivative_storage.Allocate (DimDeRequ1); |
| 439 | RationalArray = derivative_storage ; |
| 440 | } |
| 441 | |
| 442 | Inverse = 1.0e0 / PolesArray[Dimension] ; |
| 443 | Index = 0 ; |
| 444 | Index2 = - Dimension2; |
| 445 | OtherIndex = 0 ; |
| 446 | |
| 447 | for (ii = 0 ; ii <= MinDegRequ ; ii++) { |
| 448 | Index2 += Dimension; |
| 449 | Index1 = Index2; |
| 450 | |
| 451 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 452 | RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++; |
| 453 | } |
| 454 | Index -= Dimension; |
| 455 | OtherIndex ++;; |
| 456 | |
| 457 | for (jj = ii - 1 ; jj >= 0 ; jj--) { |
| 458 | Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension]; |
| 459 | |
| 460 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 461 | RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++; |
| 462 | } |
| 463 | Index -= Dimension ; |
| 464 | Index1 -= Dimension2 ; |
| 465 | } |
| 466 | |
| 467 | for (jj = ii ; jj >= 1 ; jj--) { |
| 468 | binomial_array[jj] += binomial_array[jj - 1] ; |
| 469 | } |
| 470 | |
| 471 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 472 | RationalArray[Index] *= Inverse ; Index++; |
| 473 | } |
| 474 | } |
| 475 | |
| 476 | for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){ |
| 477 | Index2 += Dimension; |
| 478 | Index1 = Index2; |
| 479 | |
| 480 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 481 | RationalArray[Index] = 0.0e0 ; Index++; |
| 482 | } |
| 483 | Index -= Dimension; |
| 484 | |
| 485 | for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) { |
| 486 | Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension]; |
| 487 | |
| 488 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 489 | RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++; |
| 490 | } |
| 491 | Index -= Dimension ; |
| 492 | Index1 -= Dimension2 ; |
| 493 | } |
| 494 | |
| 495 | for (jj = ii ; jj >= 1 ; jj--) { |
| 496 | binomial_array[jj] += binomial_array[jj - 1] ; |
| 497 | } |
| 498 | |
| 499 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 500 | RationalArray[Index] *= Inverse; Index++; |
| 501 | } |
| 502 | } |
| 503 | |
| 504 | if (!All) { |
| 505 | RationalArray = &RDers ; |
| 506 | Standard_Integer DimDeRequ = Dimension * DerivativeRequest; |
| 507 | |
| 508 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 509 | RationalArray[kk] = derivative_storage[DimDeRequ]; DimDeRequ++; |
| 510 | } |
| 511 | } |
| 512 | } |
| 513 | } |
| 514 | |
| 515 | //======================================================================= |
| 516 | //function : RationalDerivatives |
| 517 | //purpose : Uses Homogeneous Poles Derivatives and Deivatives of Weights |
| 518 | //======================================================================= |
| 519 | |
| 520 | void PLib::RationalDerivatives(const Standard_Integer DerivativeRequest, |
| 521 | const Standard_Integer Dimension, |
| 522 | Standard_Real& PolesDerivates, |
| 523 | // must be an array with |
| 524 | // (DerivativeRequest + 1) * Dimension slots |
| 525 | Standard_Real& WeightsDerivates, |
| 526 | // must be an array with |
| 527 | // (DerivativeRequest + 1) slots |
| 528 | Standard_Real& RationalDerivates) |
| 529 | { |
| 530 | // |
| 531 | // Our purpose is to compute f = (u/v) derivated N times |
| 532 | // |
| 533 | // We Write u = fv |
| 534 | // Let C(N,P) be the binomial |
| 535 | // |
| 536 | // then we have |
| 537 | // |
| 538 | // (q) (p) (q-p) |
| 539 | // u = SUM C (q,p) f v |
| 540 | // p = 0 to q |
| 541 | // |
| 542 | // |
| 543 | // Therefore |
| 544 | // |
| 545 | // |
| 546 | // (q) ( (q) (p) (q-p) ) |
| 547 | // f = (1/v) ( u - SUM C (q,p) f v ) |
| 548 | // ( p = 0 to q-1 ) |
| 549 | // |
| 550 | // |
| 551 | // make arrays for the binomial since computing it each time could |
| 552 | // raize a performance issue |
| 553 | // |
| 554 | Standard_Real Inverse; |
| 555 | Standard_Real *PolesArray = &PolesDerivates; |
| 556 | Standard_Real *WeightsArray = &WeightsDerivates; |
| 557 | Standard_Real *RationalArray = &RationalDerivates; |
| 558 | Standard_Real Factor ; |
| 559 | |
| 560 | Standard_Integer ii, Index, Index1, Index2, jj; |
| 561 | Standard_Integer DeRequest1 = DerivativeRequest + 1; |
| 562 | |
| 563 | NCollection_LocalArray<Standard_Real> binomial_array (DeRequest1); |
| 564 | NCollection_LocalArray<Standard_Real> derivative_storage; |
| 565 | |
| 566 | for (ii = 0 ; ii < DeRequest1 ; ii++) { |
| 567 | binomial_array[ii] = 1.0e0 ; |
| 568 | } |
| 569 | Inverse = 1.0e0 / WeightsArray[0] ; |
| 570 | if (Dimension == 3) { |
| 571 | Index = 0 ; |
| 572 | Index2 = - 6 ; |
| 573 | |
| 574 | for (ii = 0 ; ii < DeRequest1 ; ii++) { |
| 575 | Index2 += 3; |
| 576 | Index1 = Index2; |
| 577 | RationalArray[Index] = PolesArray[Index] ; Index++; |
| 578 | RationalArray[Index] = PolesArray[Index] ; Index++; |
| 579 | RationalArray[Index] = PolesArray[Index] ; |
| 580 | Index -= 2; |
| 581 | |
| 582 | for (jj = ii - 1 ; jj >= 0 ; jj--) { |
| 583 | Factor = binomial_array[jj] * WeightsArray[ii - jj] ; |
| 584 | RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++; |
| 585 | RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++; |
| 586 | RationalArray[Index] -= Factor * RationalArray[Index1]; |
| 587 | Index -= 2; |
| 588 | Index1 -= 5; |
| 589 | } |
| 590 | |
| 591 | for (jj = ii ; jj >= 1 ; jj--) { |
| 592 | binomial_array[jj] += binomial_array[jj - 1] ; |
| 593 | } |
| 594 | RationalArray[Index] *= Inverse ; Index++; |
| 595 | RationalArray[Index] *= Inverse ; Index++; |
| 596 | RationalArray[Index] *= Inverse ; Index++; |
| 597 | } |
| 598 | } |
| 599 | else { |
| 600 | Standard_Integer kk; |
| 601 | Standard_Integer Dimension2 = Dimension << 1; |
| 602 | Index = 0 ; |
| 603 | Index2 = - Dimension2; |
| 604 | |
| 605 | for (ii = 0 ; ii < DeRequest1 ; ii++) { |
| 606 | Index2 += Dimension; |
| 607 | Index1 = Index2; |
| 608 | |
| 609 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 610 | RationalArray[Index] = PolesArray[Index]; Index++; |
| 611 | } |
| 612 | Index -= Dimension; |
| 613 | |
| 614 | for (jj = ii - 1 ; jj >= 0 ; jj--) { |
| 615 | Factor = binomial_array[jj] * WeightsArray[ii - jj] ; |
| 616 | |
| 617 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 618 | RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++; |
| 619 | } |
| 620 | Index -= Dimension; |
| 621 | Index1 -= Dimension2; |
| 622 | } |
| 623 | |
| 624 | for (jj = ii ; jj >= 1 ; jj--) { |
| 625 | binomial_array[jj] += binomial_array[jj - 1] ; |
| 626 | } |
| 627 | |
| 628 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 629 | RationalArray[Index] *= Inverse ; Index++; |
| 630 | } |
| 631 | } |
| 632 | } |
| 633 | } |
| 634 | |
| 635 | //======================================================================= |
| 636 | //function : This evaluates a polynomial and its derivatives |
| 637 | //purpose : up to the requested order |
| 638 | //======================================================================= |
| 639 | |
| 640 | void PLib::EvalPolynomial(const Standard_Real Par, |
| 641 | const Standard_Integer DerivativeRequest, |
| 642 | const Standard_Integer Degree, |
| 643 | const Standard_Integer Dimension, |
| 644 | Standard_Real& PolynomialCoeff, |
| 645 | Standard_Real& Results) |
| 646 | // |
| 647 | // the polynomial coefficients are assumed to be stored as follows : |
| 648 | // 0 |
| 649 | // [0] [Dimension -1] X coefficient |
| 650 | // 1 |
| 651 | // [Dimension] [Dimension + Dimension -1] X coefficient |
| 652 | // 2 |
| 653 | // [2 * Dimension] [2 * Dimension + Dimension-1] X coefficient |
| 654 | // |
| 655 | // ................................................... |
| 656 | // |
| 657 | // |
| 658 | // d |
| 659 | // [d * Dimension] [d * Dimension + Dimension-1] X coefficient |
| 660 | // |
| 661 | // where d is the Degree |
| 662 | // |
| 663 | { |
| 664 | Standard_Integer DegreeDimension = Degree * Dimension; |
| 665 | |
| 666 | Standard_Integer jj; |
| 667 | Standard_Real *RA = &Results ; |
| 668 | Standard_Real *PA = &PolynomialCoeff ; |
| 669 | Standard_Real *tmpRA = RA; |
| 670 | Standard_Real *tmpPA = PA + DegreeDimension; |
| 671 | |
| 672 | switch (Dimension) { |
| 673 | |
| 674 | case 1 : { |
| 675 | *tmpRA = *tmpPA; |
| 676 | if (DerivativeRequest > 0 ) { |
| 677 | tmpRA++ ; |
| 678 | Standard_Real *valRA; |
| 679 | Standard_Integer ii, LocalRequest; |
| 680 | Standard_Integer Index1, Index2; |
| 681 | Standard_Integer MaxIndex1, MaxIndex2; |
| 682 | if (DerivativeRequest < Degree) { |
| 683 | LocalRequest = DerivativeRequest; |
| 684 | MaxIndex2 = MaxIndex1 = LocalRequest; |
| 685 | } |
| 686 | else { |
| 687 | LocalRequest = Degree; |
| 688 | MaxIndex2 = MaxIndex1 = Degree; |
| 689 | } |
| 690 | MaxIndex2 --; |
| 691 | |
| 692 | for (ii = 1; ii <= LocalRequest; ii++) { |
| 693 | *tmpRA = 0.0e0; tmpRA ++ ; |
| 694 | } |
| 695 | |
| 696 | for (jj = Degree ; jj > 0 ; jj--) { |
| 697 | tmpPA --; |
| 698 | Index1 = MaxIndex1; |
| 699 | Index2 = MaxIndex2; |
| 700 | |
| 701 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 702 | valRA = &RA[Index1]; |
| 703 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 704 | Index1 --; |
| 705 | Index2 --; |
| 706 | } |
| 707 | valRA = &RA[Index1]; |
| 708 | *valRA = Par * (*valRA) + (*tmpPA); |
| 709 | } |
| 710 | } |
| 711 | else { |
| 712 | |
| 713 | for (jj = Degree ; jj > 0 ; jj--) { |
| 714 | tmpPA --; |
| 715 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 716 | } |
| 717 | } |
| 718 | break; |
| 719 | } |
| 720 | case 2 : { |
| 721 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 722 | *tmpRA = *tmpPA; tmpRA++; |
| 723 | tmpPA --; |
| 724 | if (DerivativeRequest > 0 ) { |
| 725 | Standard_Real *valRA; |
| 726 | Standard_Integer ii, LocalRequest; |
| 727 | Standard_Integer Index1, Index2; |
| 728 | Standard_Integer MaxIndex1, MaxIndex2; |
| 729 | if (DerivativeRequest < Degree) { |
| 730 | LocalRequest = DerivativeRequest; |
| 731 | MaxIndex2 = MaxIndex1 = LocalRequest << 1; |
| 732 | } |
| 733 | else { |
| 734 | LocalRequest = Degree; |
| 735 | MaxIndex2 = MaxIndex1 = DegreeDimension; |
| 736 | } |
| 737 | MaxIndex2 -= 2; |
| 738 | |
| 739 | for (ii = 1; ii <= LocalRequest; ii++) { |
| 740 | *tmpRA = 0.0e0; tmpRA++; |
| 741 | *tmpRA = 0.0e0; tmpRA++; |
| 742 | } |
| 743 | |
| 744 | for (jj = Degree ; jj > 0 ; jj--) { |
| 745 | tmpPA -= 2; |
| 746 | |
| 747 | Index1 = MaxIndex1; |
| 748 | Index2 = MaxIndex2; |
| 749 | |
| 750 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 751 | valRA = &RA[Index1]; |
| 752 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 753 | Index1 -= 2; |
| 754 | Index2 -= 2; |
| 755 | } |
| 756 | valRA = &RA[Index1]; |
| 757 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 758 | |
| 759 | Index1 = MaxIndex1 + 1; |
| 760 | Index2 = MaxIndex2 + 1; |
| 761 | |
| 762 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 763 | valRA = &RA[Index1]; |
| 764 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 765 | Index1 -= 2; |
| 766 | Index2 -= 2; |
| 767 | } |
| 768 | valRA = &RA[Index1]; |
| 769 | *valRA = Par * (*valRA) + (*tmpPA); |
| 770 | |
| 771 | tmpPA --; |
| 772 | } |
| 773 | } |
| 774 | else { |
| 775 | |
| 776 | for (jj = Degree ; jj > 0 ; jj--) { |
| 777 | tmpPA -= 2; |
| 778 | tmpRA = RA; |
| 779 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 780 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 781 | tmpPA --; |
| 782 | } |
| 783 | } |
| 784 | break; |
| 785 | } |
| 786 | case 3 : { |
| 787 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 788 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 789 | *tmpRA = *tmpPA; tmpRA++; |
| 790 | tmpPA -= 2; |
| 791 | if (DerivativeRequest > 0 ) { |
| 792 | Standard_Real *valRA; |
| 793 | Standard_Integer ii, LocalRequest; |
| 794 | Standard_Integer Index1, Index2; |
| 795 | Standard_Integer MaxIndex1, MaxIndex2; |
| 796 | if (DerivativeRequest < Degree) { |
| 797 | LocalRequest = DerivativeRequest; |
| 798 | MaxIndex2 = MaxIndex1 = (LocalRequest << 1) + LocalRequest; |
| 799 | } |
| 800 | else { |
| 801 | LocalRequest = Degree; |
| 802 | MaxIndex2 = MaxIndex1 = DegreeDimension; |
| 803 | } |
| 804 | MaxIndex2 -= 3; |
| 805 | |
| 806 | for (ii = 1; ii <= LocalRequest; ii++) { |
| 807 | *tmpRA = 0.0e0; tmpRA++; |
| 808 | *tmpRA = 0.0e0; tmpRA++; |
| 809 | *tmpRA = 0.0e0; tmpRA++; |
| 810 | } |
| 811 | |
| 812 | for (jj = Degree ; jj > 0 ; jj--) { |
| 813 | tmpPA -= 3; |
| 814 | |
| 815 | Index1 = MaxIndex1; |
| 816 | Index2 = MaxIndex2; |
| 817 | |
| 818 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 819 | valRA = &RA[Index1]; |
| 820 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 821 | Index1 -= 3; |
| 822 | Index2 -= 3; |
| 823 | } |
| 824 | valRA = &RA[Index1]; |
| 825 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 826 | |
| 827 | Index1 = MaxIndex1 + 1; |
| 828 | Index2 = MaxIndex2 + 1; |
| 829 | |
| 830 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 831 | valRA = &RA[Index1]; |
| 832 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 833 | Index1 -= 3; |
| 834 | Index2 -= 3; |
| 835 | } |
| 836 | valRA = &RA[Index1]; |
| 837 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 838 | |
| 839 | Index1 = MaxIndex1 + 2; |
| 840 | Index2 = MaxIndex2 + 2; |
| 841 | |
| 842 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 843 | valRA = &RA[Index1]; |
| 844 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 845 | Index1 -= 3; |
| 846 | Index2 -= 3; |
| 847 | } |
| 848 | valRA = &RA[Index1]; |
| 849 | *valRA = Par * (*valRA) + (*tmpPA); |
| 850 | |
| 851 | tmpPA -= 2; |
| 852 | } |
| 853 | } |
| 854 | else { |
| 855 | |
| 856 | for (jj = Degree ; jj > 0 ; jj--) { |
| 857 | tmpPA -= 3; |
| 858 | tmpRA = RA; |
| 859 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 860 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 861 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 862 | tmpPA -= 2; |
| 863 | } |
| 864 | } |
| 865 | break; |
| 866 | } |
| 867 | case 6 : { |
| 868 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 869 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 870 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 871 | |
| 872 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 873 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 874 | *tmpRA = *tmpPA; tmpRA++; |
| 875 | tmpPA -= 5; |
| 876 | if (DerivativeRequest > 0 ) { |
| 877 | Standard_Real *valRA; |
| 878 | Standard_Integer ii, LocalRequest; |
| 879 | Standard_Integer Index1, Index2; |
| 880 | Standard_Integer MaxIndex1, MaxIndex2; |
| 881 | if (DerivativeRequest < Degree) { |
| 882 | LocalRequest = DerivativeRequest; |
| 883 | MaxIndex2 = MaxIndex1 = (LocalRequest << 2) + (LocalRequest << 1); |
| 884 | } |
| 885 | else { |
| 886 | LocalRequest = Degree; |
| 887 | MaxIndex2 = MaxIndex1 = DegreeDimension; |
| 888 | } |
| 889 | MaxIndex2 -= 6; |
| 890 | |
| 891 | for (ii = 1; ii <= LocalRequest; ii++) { |
| 892 | *tmpRA = 0.0e0; tmpRA++; |
| 893 | *tmpRA = 0.0e0; tmpRA++; |
| 894 | *tmpRA = 0.0e0; tmpRA++; |
| 895 | |
| 896 | *tmpRA = 0.0e0; tmpRA++; |
| 897 | *tmpRA = 0.0e0; tmpRA++; |
| 898 | *tmpRA = 0.0e0; tmpRA++; |
| 899 | } |
| 900 | |
| 901 | for (jj = Degree ; jj > 0 ; jj--) { |
| 902 | tmpPA -= 6; |
| 903 | |
| 904 | Index1 = MaxIndex1; |
| 905 | Index2 = MaxIndex2; |
| 906 | |
| 907 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 908 | valRA = &RA[Index1]; |
| 909 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 910 | Index1 -= 6; |
| 911 | Index2 -= 6; |
| 912 | } |
| 913 | valRA = &RA[Index1]; |
| 914 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 915 | |
| 916 | Index1 = MaxIndex1 + 1; |
| 917 | Index2 = MaxIndex2 + 1; |
| 918 | |
| 919 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 920 | valRA = &RA[Index1]; |
| 921 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 922 | Index1 -= 6; |
| 923 | Index2 -= 6; |
| 924 | } |
| 925 | valRA = &RA[Index1]; |
| 926 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 927 | |
| 928 | Index1 = MaxIndex1 + 2; |
| 929 | Index2 = MaxIndex2 + 2; |
| 930 | |
| 931 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 932 | valRA = &RA[Index1]; |
| 933 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 934 | Index1 -= 6; |
| 935 | Index2 -= 6; |
| 936 | } |
| 937 | valRA = &RA[Index1]; |
| 938 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 939 | |
| 940 | Index1 = MaxIndex1 + 3; |
| 941 | Index2 = MaxIndex2 + 3; |
| 942 | |
| 943 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 944 | valRA = &RA[Index1]; |
| 945 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 946 | Index1 -= 6; |
| 947 | Index2 -= 6; |
| 948 | } |
| 949 | valRA = &RA[Index1]; |
| 950 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 951 | |
| 952 | Index1 = MaxIndex1 + 4; |
| 953 | Index2 = MaxIndex2 + 4; |
| 954 | |
| 955 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 956 | valRA = &RA[Index1]; |
| 957 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 958 | Index1 -= 6; |
| 959 | Index2 -= 6; |
| 960 | } |
| 961 | valRA = &RA[Index1]; |
| 962 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 963 | |
| 964 | Index1 = MaxIndex1 + 5; |
| 965 | Index2 = MaxIndex2 + 5; |
| 966 | |
| 967 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 968 | valRA = &RA[Index1]; |
| 969 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 970 | Index1 -= 6; |
| 971 | Index2 -= 6; |
| 972 | } |
| 973 | valRA = &RA[Index1]; |
| 974 | *valRA = Par * (*valRA) + (*tmpPA); |
| 975 | |
| 976 | tmpPA -= 5; |
| 977 | } |
| 978 | } |
| 979 | else { |
| 980 | |
| 981 | for (jj = Degree ; jj > 0 ; jj--) { |
| 982 | tmpPA -= 6; |
| 983 | tmpRA = RA; |
| 984 | |
| 985 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 986 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 987 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 988 | |
| 989 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 990 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 991 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 992 | tmpPA -= 5; |
| 993 | } |
| 994 | } |
| 995 | break; |
| 996 | } |
| 997 | case 9 : { |
| 998 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 999 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1000 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1001 | |
| 1002 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1003 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1004 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1005 | |
| 1006 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1007 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1008 | *tmpRA = *tmpPA; tmpRA++; |
| 1009 | tmpPA -= 8; |
| 1010 | if (DerivativeRequest > 0 ) { |
| 1011 | Standard_Real *valRA; |
| 1012 | Standard_Integer ii, LocalRequest; |
| 1013 | Standard_Integer Index1, Index2; |
| 1014 | Standard_Integer MaxIndex1, MaxIndex2; |
| 1015 | if (DerivativeRequest < Degree) { |
| 1016 | LocalRequest = DerivativeRequest; |
| 1017 | MaxIndex2 = MaxIndex1 = (LocalRequest << 3) + LocalRequest; |
| 1018 | } |
| 1019 | else { |
| 1020 | LocalRequest = Degree; |
| 1021 | MaxIndex2 = MaxIndex1 = DegreeDimension; |
| 1022 | } |
| 1023 | MaxIndex2 -= 9; |
| 1024 | |
| 1025 | for (ii = 1; ii <= LocalRequest; ii++) { |
| 1026 | *tmpRA = 0.0e0; tmpRA++; |
| 1027 | *tmpRA = 0.0e0; tmpRA++; |
| 1028 | *tmpRA = 0.0e0; tmpRA++; |
| 1029 | |
| 1030 | *tmpRA = 0.0e0; tmpRA++; |
| 1031 | *tmpRA = 0.0e0; tmpRA++; |
| 1032 | *tmpRA = 0.0e0; tmpRA++; |
| 1033 | |
| 1034 | *tmpRA = 0.0e0; tmpRA++; |
| 1035 | *tmpRA = 0.0e0; tmpRA++; |
| 1036 | *tmpRA = 0.0e0; tmpRA++; |
| 1037 | } |
| 1038 | |
| 1039 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1040 | tmpPA -= 9; |
| 1041 | |
| 1042 | Index1 = MaxIndex1; |
| 1043 | Index2 = MaxIndex2; |
| 1044 | |
| 1045 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1046 | valRA = &RA[Index1]; |
| 1047 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1048 | Index1 -= 9; |
| 1049 | Index2 -= 9; |
| 1050 | } |
| 1051 | valRA = &RA[Index1]; |
| 1052 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1053 | |
| 1054 | Index1 = MaxIndex1 + 1; |
| 1055 | Index2 = MaxIndex2 + 1; |
| 1056 | |
| 1057 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1058 | valRA = &RA[Index1]; |
| 1059 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1060 | Index1 -= 9; |
| 1061 | Index2 -= 9; |
| 1062 | } |
| 1063 | valRA = &RA[Index1]; |
| 1064 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1065 | |
| 1066 | Index1 = MaxIndex1 + 2; |
| 1067 | Index2 = MaxIndex2 + 2; |
| 1068 | |
| 1069 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1070 | valRA = &RA[Index1]; |
| 1071 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1072 | Index1 -= 9; |
| 1073 | Index2 -= 9; |
| 1074 | } |
| 1075 | valRA = &RA[Index1]; |
| 1076 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1077 | |
| 1078 | Index1 = MaxIndex1 + 3; |
| 1079 | Index2 = MaxIndex2 + 3; |
| 1080 | |
| 1081 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1082 | valRA = &RA[Index1]; |
| 1083 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1084 | Index1 -= 9; |
| 1085 | Index2 -= 9; |
| 1086 | } |
| 1087 | valRA = &RA[Index1]; |
| 1088 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1089 | |
| 1090 | Index1 = MaxIndex1 + 4; |
| 1091 | Index2 = MaxIndex2 + 4; |
| 1092 | |
| 1093 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1094 | valRA = &RA[Index1]; |
| 1095 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1096 | Index1 -= 9; |
| 1097 | Index2 -= 9; |
| 1098 | } |
| 1099 | valRA = &RA[Index1]; |
| 1100 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1101 | |
| 1102 | Index1 = MaxIndex1 + 5; |
| 1103 | Index2 = MaxIndex2 + 5; |
| 1104 | |
| 1105 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1106 | valRA = &RA[Index1]; |
| 1107 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1108 | Index1 -= 9; |
| 1109 | Index2 -= 9; |
| 1110 | } |
| 1111 | valRA = &RA[Index1]; |
| 1112 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1113 | |
| 1114 | Index1 = MaxIndex1 + 6; |
| 1115 | Index2 = MaxIndex2 + 6; |
| 1116 | |
| 1117 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1118 | valRA = &RA[Index1]; |
| 1119 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1120 | Index1 -= 9; |
| 1121 | Index2 -= 9; |
| 1122 | } |
| 1123 | valRA = &RA[Index1]; |
| 1124 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1125 | |
| 1126 | Index1 = MaxIndex1 + 7; |
| 1127 | Index2 = MaxIndex2 + 7; |
| 1128 | |
| 1129 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1130 | valRA = &RA[Index1]; |
| 1131 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1132 | Index1 -= 9; |
| 1133 | Index2 -= 9; |
| 1134 | } |
| 1135 | valRA = &RA[Index1]; |
| 1136 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1137 | |
| 1138 | Index1 = MaxIndex1 + 8; |
| 1139 | Index2 = MaxIndex2 + 8; |
| 1140 | |
| 1141 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1142 | valRA = &RA[Index1]; |
| 1143 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1144 | Index1 -= 9; |
| 1145 | Index2 -= 9; |
| 1146 | } |
| 1147 | valRA = &RA[Index1]; |
| 1148 | *valRA = Par * (*valRA) + (*tmpPA); |
| 1149 | |
| 1150 | tmpPA -= 8; |
| 1151 | } |
| 1152 | } |
| 1153 | else { |
| 1154 | |
| 1155 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1156 | tmpPA -= 9; |
| 1157 | tmpRA = RA; |
| 1158 | |
| 1159 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1160 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1161 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1162 | |
| 1163 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1164 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1165 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1166 | |
| 1167 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1168 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1169 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 1170 | tmpPA -= 8; |
| 1171 | } |
| 1172 | } |
| 1173 | break; |
| 1174 | } |
| 1175 | case 12 : { |
| 1176 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1177 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1178 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1179 | |
| 1180 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1181 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1182 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1183 | |
| 1184 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1185 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1186 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1187 | |
| 1188 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1189 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1190 | *tmpRA = *tmpPA; tmpRA++; |
| 1191 | tmpPA -= 11; |
| 1192 | if (DerivativeRequest > 0 ) { |
| 1193 | Standard_Real *valRA; |
| 1194 | Standard_Integer ii, LocalRequest; |
| 1195 | Standard_Integer Index1, Index2; |
| 1196 | Standard_Integer MaxIndex1, MaxIndex2; |
| 1197 | if (DerivativeRequest < Degree) { |
| 1198 | LocalRequest = DerivativeRequest; |
| 1199 | MaxIndex2 = MaxIndex1 = (LocalRequest << 3) + (LocalRequest << 2); |
| 1200 | } |
| 1201 | else { |
| 1202 | LocalRequest = Degree; |
| 1203 | MaxIndex2 = MaxIndex1 = DegreeDimension; |
| 1204 | } |
| 1205 | MaxIndex2 -= 12; |
| 1206 | |
| 1207 | for (ii = 1; ii <= LocalRequest; ii++) { |
| 1208 | *tmpRA = 0.0e0; tmpRA++; |
| 1209 | *tmpRA = 0.0e0; tmpRA++; |
| 1210 | *tmpRA = 0.0e0; tmpRA++; |
| 1211 | |
| 1212 | *tmpRA = 0.0e0; tmpRA++; |
| 1213 | *tmpRA = 0.0e0; tmpRA++; |
| 1214 | *tmpRA = 0.0e0; tmpRA++; |
| 1215 | |
| 1216 | *tmpRA = 0.0e0; tmpRA++; |
| 1217 | *tmpRA = 0.0e0; tmpRA++; |
| 1218 | *tmpRA = 0.0e0; tmpRA++; |
| 1219 | |
| 1220 | *tmpRA = 0.0e0; tmpRA++; |
| 1221 | *tmpRA = 0.0e0; tmpRA++; |
| 1222 | *tmpRA = 0.0e0; tmpRA++; |
| 1223 | } |
| 1224 | |
| 1225 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1226 | tmpPA -= 12; |
| 1227 | |
| 1228 | Index1 = MaxIndex1; |
| 1229 | Index2 = MaxIndex2; |
| 1230 | |
| 1231 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1232 | valRA = &RA[Index1]; |
| 1233 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1234 | Index1 -= 12; |
| 1235 | Index2 -= 12; |
| 1236 | } |
| 1237 | valRA = &RA[Index1]; |
| 1238 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1239 | |
| 1240 | Index1 = MaxIndex1 + 1; |
| 1241 | Index2 = MaxIndex2 + 1; |
| 1242 | |
| 1243 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1244 | valRA = &RA[Index1]; |
| 1245 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1246 | Index1 -= 12; |
| 1247 | Index2 -= 12; |
| 1248 | } |
| 1249 | valRA = &RA[Index1]; |
| 1250 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1251 | |
| 1252 | Index1 = MaxIndex1 + 2; |
| 1253 | Index2 = MaxIndex2 + 2; |
| 1254 | |
| 1255 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1256 | valRA = &RA[Index1]; |
| 1257 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1258 | Index1 -= 12; |
| 1259 | Index2 -= 12; |
| 1260 | } |
| 1261 | valRA = &RA[Index1]; |
| 1262 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1263 | |
| 1264 | Index1 = MaxIndex1 + 3; |
| 1265 | Index2 = MaxIndex2 + 3; |
| 1266 | |
| 1267 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1268 | valRA = &RA[Index1]; |
| 1269 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1270 | Index1 -= 12; |
| 1271 | Index2 -= 12; |
| 1272 | } |
| 1273 | valRA = &RA[Index1]; |
| 1274 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1275 | |
| 1276 | Index1 = MaxIndex1 + 4; |
| 1277 | Index2 = MaxIndex2 + 4; |
| 1278 | |
| 1279 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1280 | valRA = &RA[Index1]; |
| 1281 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1282 | Index1 -= 12; |
| 1283 | Index2 -= 12; |
| 1284 | } |
| 1285 | valRA = &RA[Index1]; |
| 1286 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1287 | |
| 1288 | Index1 = MaxIndex1 + 5; |
| 1289 | Index2 = MaxIndex2 + 5; |
| 1290 | |
| 1291 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1292 | valRA = &RA[Index1]; |
| 1293 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1294 | Index1 -= 12; |
| 1295 | Index2 -= 12; |
| 1296 | } |
| 1297 | valRA = &RA[Index1]; |
| 1298 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1299 | |
| 1300 | Index1 = MaxIndex1 + 6; |
| 1301 | Index2 = MaxIndex2 + 6; |
| 1302 | |
| 1303 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1304 | valRA = &RA[Index1]; |
| 1305 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1306 | Index1 -= 12; |
| 1307 | Index2 -= 12; |
| 1308 | } |
| 1309 | valRA = &RA[Index1]; |
| 1310 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1311 | |
| 1312 | Index1 = MaxIndex1 + 7; |
| 1313 | Index2 = MaxIndex2 + 7; |
| 1314 | |
| 1315 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1316 | valRA = &RA[Index1]; |
| 1317 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1318 | Index1 -= 12; |
| 1319 | Index2 -= 12; |
| 1320 | } |
| 1321 | valRA = &RA[Index1]; |
| 1322 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1323 | |
| 1324 | Index1 = MaxIndex1 + 8; |
| 1325 | Index2 = MaxIndex2 + 8; |
| 1326 | |
| 1327 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1328 | valRA = &RA[Index1]; |
| 1329 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1330 | Index1 -= 12; |
| 1331 | Index2 -= 12; |
| 1332 | } |
| 1333 | valRA = &RA[Index1]; |
| 1334 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1335 | |
| 1336 | Index1 = MaxIndex1 + 9; |
| 1337 | Index2 = MaxIndex2 + 9; |
| 1338 | |
| 1339 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1340 | valRA = &RA[Index1]; |
| 1341 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1342 | Index1 -= 12; |
| 1343 | Index2 -= 12; |
| 1344 | } |
| 1345 | valRA = &RA[Index1]; |
| 1346 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1347 | |
| 1348 | Index1 = MaxIndex1 + 10; |
| 1349 | Index2 = MaxIndex2 + 10; |
| 1350 | |
| 1351 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1352 | valRA = &RA[Index1]; |
| 1353 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1354 | Index1 -= 12; |
| 1355 | Index2 -= 12; |
| 1356 | } |
| 1357 | valRA = &RA[Index1]; |
| 1358 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1359 | |
| 1360 | Index1 = MaxIndex1 + 11; |
| 1361 | Index2 = MaxIndex2 + 11; |
| 1362 | |
| 1363 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1364 | valRA = &RA[Index1]; |
| 1365 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1366 | Index1 -= 12; |
| 1367 | Index2 -= 12; |
| 1368 | } |
| 1369 | valRA = &RA[Index1]; |
| 1370 | *valRA = Par * (*valRA) + (*tmpPA); |
| 1371 | |
| 1372 | tmpPA -= 11; |
| 1373 | } |
| 1374 | } |
| 1375 | else { |
| 1376 | |
| 1377 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1378 | tmpPA -= 12; |
| 1379 | tmpRA = RA; |
| 1380 | |
| 1381 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1382 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1383 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1384 | |
| 1385 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1386 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1387 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1388 | |
| 1389 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1390 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1391 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1392 | |
| 1393 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1394 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1395 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 1396 | tmpPA -= 11; |
| 1397 | } |
| 1398 | } |
| 1399 | break; |
| 1400 | } |
| 1401 | case 15 : { |
| 1402 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1403 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1404 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1405 | |
| 1406 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1407 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1408 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1409 | |
| 1410 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1411 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1412 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1413 | |
| 1414 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1415 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1416 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1417 | |
| 1418 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1419 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1420 | *tmpRA = *tmpPA; tmpRA++; |
| 1421 | tmpPA -= 14; |
| 1422 | if (DerivativeRequest > 0 ) { |
| 1423 | Standard_Real *valRA; |
| 1424 | Standard_Integer ii, LocalRequest; |
| 1425 | Standard_Integer Index1, Index2; |
| 1426 | Standard_Integer MaxIndex1, MaxIndex2; |
| 1427 | if (DerivativeRequest < Degree) { |
| 1428 | LocalRequest = DerivativeRequest; |
| 1429 | MaxIndex2 = MaxIndex1 = (LocalRequest << 4) - LocalRequest; |
| 1430 | } |
| 1431 | else { |
| 1432 | LocalRequest = Degree; |
| 1433 | MaxIndex2 = MaxIndex1 = DegreeDimension; |
| 1434 | } |
| 1435 | MaxIndex2 -= 15; |
| 1436 | |
| 1437 | for (ii = 1; ii <= LocalRequest; ii++) { |
| 1438 | *tmpRA = 0.0e0; tmpRA++; |
| 1439 | *tmpRA = 0.0e0; tmpRA++; |
| 1440 | *tmpRA = 0.0e0; tmpRA++; |
| 1441 | |
| 1442 | *tmpRA = 0.0e0; tmpRA++; |
| 1443 | *tmpRA = 0.0e0; tmpRA++; |
| 1444 | *tmpRA = 0.0e0; tmpRA++; |
| 1445 | |
| 1446 | *tmpRA = 0.0e0; tmpRA++; |
| 1447 | *tmpRA = 0.0e0; tmpRA++; |
| 1448 | *tmpRA = 0.0e0; tmpRA++; |
| 1449 | |
| 1450 | *tmpRA = 0.0e0; tmpRA++; |
| 1451 | *tmpRA = 0.0e0; tmpRA++; |
| 1452 | *tmpRA = 0.0e0; tmpRA++; |
| 1453 | |
| 1454 | *tmpRA = 0.0e0; tmpRA++; |
| 1455 | *tmpRA = 0.0e0; tmpRA++; |
| 1456 | *tmpRA = 0.0e0; tmpRA++; |
| 1457 | } |
| 1458 | |
| 1459 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1460 | tmpPA -= 15; |
| 1461 | |
| 1462 | Index1 = MaxIndex1; |
| 1463 | Index2 = MaxIndex2; |
| 1464 | |
| 1465 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1466 | valRA = &RA[Index1]; |
| 1467 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1468 | Index1 -= 15; |
| 1469 | Index2 -= 15; |
| 1470 | } |
| 1471 | valRA = &RA[Index1]; |
| 1472 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1473 | |
| 1474 | Index1 = MaxIndex1 + 1; |
| 1475 | Index2 = MaxIndex2 + 1; |
| 1476 | |
| 1477 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1478 | valRA = &RA[Index1]; |
| 1479 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1480 | Index1 -= 15; |
| 1481 | Index2 -= 15; |
| 1482 | } |
| 1483 | valRA = &RA[Index1]; |
| 1484 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1485 | |
| 1486 | Index1 = MaxIndex1 + 2; |
| 1487 | Index2 = MaxIndex2 + 2; |
| 1488 | |
| 1489 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1490 | valRA = &RA[Index1]; |
| 1491 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1492 | Index1 -= 15; |
| 1493 | Index2 -= 15; |
| 1494 | } |
| 1495 | valRA = &RA[Index1]; |
| 1496 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1497 | |
| 1498 | Index1 = MaxIndex1 + 3; |
| 1499 | Index2 = MaxIndex2 + 3; |
| 1500 | |
| 1501 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1502 | valRA = &RA[Index1]; |
| 1503 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1504 | Index1 -= 15; |
| 1505 | Index2 -= 15; |
| 1506 | } |
| 1507 | valRA = &RA[Index1]; |
| 1508 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1509 | |
| 1510 | Index1 = MaxIndex1 + 4; |
| 1511 | Index2 = MaxIndex2 + 4; |
| 1512 | |
| 1513 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1514 | valRA = &RA[Index1]; |
| 1515 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1516 | Index1 -= 15; |
| 1517 | Index2 -= 15; |
| 1518 | } |
| 1519 | valRA = &RA[Index1]; |
| 1520 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1521 | |
| 1522 | Index1 = MaxIndex1 + 5; |
| 1523 | Index2 = MaxIndex2 + 5; |
| 1524 | |
| 1525 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1526 | valRA = &RA[Index1]; |
| 1527 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1528 | Index1 -= 15; |
| 1529 | Index2 -= 15; |
| 1530 | } |
| 1531 | valRA = &RA[Index1]; |
| 1532 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1533 | |
| 1534 | Index1 = MaxIndex1 + 6; |
| 1535 | Index2 = MaxIndex2 + 6; |
| 1536 | |
| 1537 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1538 | valRA = &RA[Index1]; |
| 1539 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1540 | Index1 -= 15; |
| 1541 | Index2 -= 15; |
| 1542 | } |
| 1543 | valRA = &RA[Index1]; |
| 1544 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1545 | |
| 1546 | Index1 = MaxIndex1 + 7; |
| 1547 | Index2 = MaxIndex2 + 7; |
| 1548 | |
| 1549 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1550 | valRA = &RA[Index1]; |
| 1551 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1552 | Index1 -= 15; |
| 1553 | Index2 -= 15; |
| 1554 | } |
| 1555 | valRA = &RA[Index1]; |
| 1556 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1557 | |
| 1558 | Index1 = MaxIndex1 + 8; |
| 1559 | Index2 = MaxIndex2 + 8; |
| 1560 | |
| 1561 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1562 | valRA = &RA[Index1]; |
| 1563 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1564 | Index1 -= 15; |
| 1565 | Index2 -= 15; |
| 1566 | } |
| 1567 | valRA = &RA[Index1]; |
| 1568 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1569 | |
| 1570 | Index1 = MaxIndex1 + 9; |
| 1571 | Index2 = MaxIndex2 + 9; |
| 1572 | |
| 1573 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1574 | valRA = &RA[Index1]; |
| 1575 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1576 | Index1 -= 15; |
| 1577 | Index2 -= 15; |
| 1578 | } |
| 1579 | valRA = &RA[Index1]; |
| 1580 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1581 | |
| 1582 | Index1 = MaxIndex1 + 10; |
| 1583 | Index2 = MaxIndex2 + 10; |
| 1584 | |
| 1585 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1586 | valRA = &RA[Index1]; |
| 1587 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1588 | Index1 -= 15; |
| 1589 | Index2 -= 15; |
| 1590 | } |
| 1591 | valRA = &RA[Index1]; |
| 1592 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1593 | |
| 1594 | Index1 = MaxIndex1 + 11; |
| 1595 | Index2 = MaxIndex2 + 11; |
| 1596 | |
| 1597 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1598 | valRA = &RA[Index1]; |
| 1599 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1600 | Index1 -= 15; |
| 1601 | Index2 -= 15; |
| 1602 | } |
| 1603 | valRA = &RA[Index1]; |
| 1604 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1605 | |
| 1606 | Index1 = MaxIndex1 + 12; |
| 1607 | Index2 = MaxIndex2 + 12; |
| 1608 | |
| 1609 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1610 | valRA = &RA[Index1]; |
| 1611 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1612 | Index1 -= 15; |
| 1613 | Index2 -= 15; |
| 1614 | } |
| 1615 | valRA = &RA[Index1]; |
| 1616 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1617 | |
| 1618 | Index1 = MaxIndex1 + 13; |
| 1619 | Index2 = MaxIndex2 + 13; |
| 1620 | |
| 1621 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1622 | valRA = &RA[Index1]; |
| 1623 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1624 | Index1 -= 15; |
| 1625 | Index2 -= 15; |
| 1626 | } |
| 1627 | valRA = &RA[Index1]; |
| 1628 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1629 | |
| 1630 | Index1 = MaxIndex1 + 14; |
| 1631 | Index2 = MaxIndex2 + 14; |
| 1632 | |
| 1633 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1634 | valRA = &RA[Index1]; |
| 1635 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1636 | Index1 -= 15; |
| 1637 | Index2 -= 15; |
| 1638 | } |
| 1639 | valRA = &RA[Index1]; |
| 1640 | *valRA = Par * (*valRA) + (*tmpPA); |
| 1641 | |
| 1642 | tmpPA -= 14; |
| 1643 | } |
| 1644 | } |
| 1645 | else { |
| 1646 | |
| 1647 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1648 | tmpPA -= 15; |
| 1649 | tmpRA = RA; |
| 1650 | |
| 1651 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1652 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1653 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1654 | |
| 1655 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1656 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1657 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1658 | |
| 1659 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1660 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1661 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1662 | |
| 1663 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1664 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1665 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1666 | |
| 1667 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1668 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1669 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 1670 | tmpPA -= 14; |
| 1671 | } |
| 1672 | } |
| 1673 | break; |
| 1674 | } |
| 1675 | default : { |
| 1676 | Standard_Integer kk ; |
| 1677 | for ( kk = 0 ; kk < Dimension ; kk++) { |
| 1678 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1679 | } |
| 1680 | tmpPA -= Dimension; |
| 1681 | if (DerivativeRequest > 0 ) { |
| 1682 | Standard_Real *valRA; |
| 1683 | Standard_Integer ii, LocalRequest; |
| 1684 | Standard_Integer Index1, Index2; |
| 1685 | Standard_Integer MaxIndex1, MaxIndex2; |
| 1686 | if (DerivativeRequest < Degree) { |
| 1687 | LocalRequest = DerivativeRequest; |
| 1688 | MaxIndex2 = MaxIndex1 = LocalRequest * Dimension; |
| 1689 | } |
| 1690 | else { |
| 1691 | LocalRequest = Degree; |
| 1692 | MaxIndex2 = MaxIndex1 = DegreeDimension; |
| 1693 | } |
| 1694 | MaxIndex2 -= Dimension; |
| 1695 | |
| 1696 | for (ii = 1; ii <= MaxIndex1; ii++) { |
| 1697 | *tmpRA = 0.0e0; tmpRA++; |
| 1698 | } |
| 1699 | |
| 1700 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1701 | tmpPA -= Dimension; |
| 1702 | |
| 1703 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 1704 | Index1 = MaxIndex1 + kk; |
| 1705 | Index2 = MaxIndex2 + kk; |
| 1706 | |
| 1707 | for (ii = LocalRequest ; ii > 0 ; ii--) { |
| 1708 | valRA = &RA[Index1]; |
| 1709 | *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; |
| 1710 | Index1 -= Dimension; |
| 1711 | Index2 -= Dimension; |
| 1712 | } |
| 1713 | valRA = &RA[Index1]; |
| 1714 | *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; |
| 1715 | } |
| 1716 | tmpPA -= Dimension; |
| 1717 | } |
| 1718 | } |
| 1719 | else { |
| 1720 | |
| 1721 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1722 | tmpPA -= Dimension; |
| 1723 | tmpRA = RA; |
| 1724 | |
| 1725 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 1726 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1727 | } |
| 1728 | tmpPA -= Dimension; |
| 1729 | } |
| 1730 | } |
| 1731 | } |
| 1732 | } |
| 1733 | } |
| 1734 | |
| 1735 | //======================================================================= |
| 1736 | //function : This evaluates a polynomial without derivative |
| 1737 | //purpose : |
| 1738 | //======================================================================= |
| 1739 | |
| 1740 | void PLib::NoDerivativeEvalPolynomial(const Standard_Real Par, |
| 1741 | const Standard_Integer Degree, |
| 1742 | const Standard_Integer Dimension, |
| 1743 | const Standard_Integer DegreeDimension, |
| 1744 | Standard_Real& PolynomialCoeff, |
| 1745 | Standard_Real& Results) |
| 1746 | { |
| 1747 | Standard_Integer jj; |
| 1748 | Standard_Real *RA = &Results ; |
| 1749 | Standard_Real *PA = &PolynomialCoeff ; |
| 1750 | Standard_Real *tmpRA = RA; |
| 1751 | Standard_Real *tmpPA = PA + DegreeDimension; |
| 1752 | |
| 1753 | switch (Dimension) { |
| 1754 | |
| 1755 | case 1 : { |
| 1756 | *tmpRA = *tmpPA; |
| 1757 | |
| 1758 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1759 | tmpPA--; |
| 1760 | |
| 1761 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 1762 | } |
| 1763 | break; |
| 1764 | } |
| 1765 | case 2 : { |
| 1766 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1767 | *tmpRA = *tmpPA; |
| 1768 | tmpPA--; |
| 1769 | |
| 1770 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1771 | tmpPA -= 2; |
| 1772 | tmpRA = RA; |
| 1773 | |
| 1774 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1775 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 1776 | tmpPA--; |
| 1777 | } |
| 1778 | break; |
| 1779 | } |
| 1780 | case 3 : { |
| 1781 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1782 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1783 | *tmpRA = *tmpPA; |
| 1784 | tmpPA -= 2; |
| 1785 | |
| 1786 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1787 | tmpPA -= 3; |
| 1788 | tmpRA = RA; |
| 1789 | |
| 1790 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1791 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1792 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 1793 | tmpPA -= 2; |
| 1794 | } |
| 1795 | break; |
| 1796 | } |
| 1797 | case 6 : { |
| 1798 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1799 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1800 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1801 | |
| 1802 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1803 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1804 | *tmpRA = *tmpPA; |
| 1805 | tmpPA -= 5; |
| 1806 | |
| 1807 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1808 | tmpPA -= 6; |
| 1809 | tmpRA = RA; |
| 1810 | |
| 1811 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1812 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1813 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1814 | |
| 1815 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1816 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1817 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 1818 | tmpPA -= 5; |
| 1819 | } |
| 1820 | break; |
| 1821 | } |
| 1822 | case 9 : { |
| 1823 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1824 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1825 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1826 | |
| 1827 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1828 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1829 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1830 | |
| 1831 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1832 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1833 | *tmpRA = *tmpPA; |
| 1834 | tmpPA -= 8; |
| 1835 | |
| 1836 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1837 | tmpPA -= 9; |
| 1838 | tmpRA = RA; |
| 1839 | |
| 1840 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1841 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1842 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1843 | |
| 1844 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1845 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1846 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1847 | |
| 1848 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1849 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1850 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 1851 | tmpPA -= 8; |
| 1852 | } |
| 1853 | break; |
| 1854 | } |
| 1855 | case 12 : { |
| 1856 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1857 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1858 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1859 | |
| 1860 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1861 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1862 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1863 | |
| 1864 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1865 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1866 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1867 | |
| 1868 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1869 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1870 | *tmpRA = *tmpPA; |
| 1871 | tmpPA -= 11; |
| 1872 | |
| 1873 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1874 | tmpPA -= 12; |
| 1875 | tmpRA = RA; |
| 1876 | |
| 1877 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1878 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1879 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1880 | |
| 1881 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1882 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1883 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1884 | |
| 1885 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1886 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1887 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1888 | |
| 1889 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1890 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1891 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 1892 | tmpPA -= 11; |
| 1893 | } |
| 1894 | break; |
| 1895 | } |
| 1896 | case 15 : { |
| 1897 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1898 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1899 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1900 | |
| 1901 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1902 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1903 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1904 | |
| 1905 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1906 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1907 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1908 | |
| 1909 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1910 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1911 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1912 | |
| 1913 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1914 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1915 | *tmpRA = *tmpPA; |
| 1916 | tmpPA -= 14; |
| 1917 | |
| 1918 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1919 | tmpPA -= 15; |
| 1920 | tmpRA = RA; |
| 1921 | |
| 1922 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1923 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1924 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1925 | |
| 1926 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1927 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1928 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1929 | |
| 1930 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1931 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1932 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1933 | |
| 1934 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1935 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1936 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1937 | |
| 1938 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1939 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1940 | *tmpRA = Par * (*tmpRA) + (*tmpPA); |
| 1941 | tmpPA -= 14; |
| 1942 | } |
| 1943 | break; |
| 1944 | } |
| 1945 | default : { |
| 1946 | Standard_Integer kk ; |
| 1947 | for ( kk = 0 ; kk < Dimension ; kk++) { |
| 1948 | *tmpRA = *tmpPA; tmpRA++; tmpPA++; |
| 1949 | } |
| 1950 | tmpPA -= Dimension; |
| 1951 | |
| 1952 | for (jj = Degree ; jj > 0 ; jj--) { |
| 1953 | tmpPA -= Dimension; |
| 1954 | tmpRA = RA; |
| 1955 | |
| 1956 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 1957 | *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; |
| 1958 | } |
| 1959 | tmpPA -= Dimension; |
| 1960 | } |
| 1961 | } |
| 1962 | } |
| 1963 | } |
| 1964 | |
| 1965 | //======================================================================= |
| 1966 | //function : This evaluates a polynomial of 2 variables |
| 1967 | //purpose : or its derivative at the requested orders |
| 1968 | //======================================================================= |
| 1969 | |
| 1970 | void PLib::EvalPoly2Var(const Standard_Real UParameter, |
| 1971 | const Standard_Real VParameter, |
| 1972 | const Standard_Integer UDerivativeRequest, |
| 1973 | const Standard_Integer VDerivativeRequest, |
| 1974 | const Standard_Integer UDegree, |
| 1975 | const Standard_Integer VDegree, |
| 1976 | const Standard_Integer Dimension, |
| 1977 | Standard_Real& PolynomialCoeff, |
| 1978 | Standard_Real& Results) |
| 1979 | // |
| 1980 | // the polynomial coefficients are assumed to be stored as follows : |
| 1981 | // 0 0 |
| 1982 | // [0] [Dimension -1] U V coefficient |
| 1983 | // 1 0 |
| 1984 | // [Dimension] [Dimension + Dimension -1] U V coefficient |
| 1985 | // 2 0 |
| 1986 | // [2 * Dimension] [2 * Dimension + Dimension-1] U V coefficient |
| 1987 | // |
| 1988 | // ................................................... |
| 1989 | // |
| 1990 | // |
| 1991 | // m 0 |
| 1992 | // [m * Dimension] [m * Dimension + Dimension-1] U V coefficient |
| 1993 | // |
| 1994 | // where m = UDegree |
| 1995 | // |
| 1996 | // 0 1 |
| 1997 | // [(m+1) * Dimension] [(m+1) * Dimension + Dimension-1] U V coefficient |
| 1998 | // |
| 1999 | // ................................................... |
| 2000 | // |
| 2001 | // m 1 |
| 2002 | // [2*m * Dimension] [2*m * Dimension + Dimension-1] U V coefficient |
| 2003 | // |
| 2004 | // ................................................... |
| 2005 | // |
| 2006 | // m n |
| 2007 | // [m*n * Dimension] [m*n * Dimension + Dimension-1] U V coefficient |
| 2008 | // |
| 2009 | // where n = VDegree |
| 2010 | { |
| 2011 | Standard_Integer Udim = (VDegree+1)*Dimension, |
| 2012 | index = Udim*UDerivativeRequest; |
| 2013 | TColStd_Array1OfReal Curve(1, Udim*(UDerivativeRequest+1)); |
| 2014 | TColStd_Array1OfReal Point(1, Dimension*(VDerivativeRequest+1)); |
| 2015 | Standard_Real * Result = (Standard_Real *) &Curve.ChangeValue(1); |
| 2016 | Standard_Real * Digit = (Standard_Real *) &Point.ChangeValue(1); |
| 2017 | Standard_Real * ResultArray ; |
| 2018 | ResultArray = &Results ; |
| 2019 | |
| 2020 | PLib::EvalPolynomial(UParameter, |
| 2021 | UDerivativeRequest, |
| 2022 | UDegree, |
| 2023 | Udim, |
| 2024 | PolynomialCoeff, |
| 2025 | Result[0]); |
| 2026 | |
| 2027 | PLib::EvalPolynomial(VParameter, |
| 2028 | VDerivativeRequest, |
| 2029 | VDegree, |
| 2030 | Dimension, |
| 2031 | Result[index], |
| 2032 | Digit[0]); |
| 2033 | |
| 2034 | index = Dimension*VDerivativeRequest; |
| 2035 | |
| 2036 | for (Standard_Integer i=0;i<Dimension;i++) { |
| 2037 | ResultArray[i] = Digit[index+i]; |
| 2038 | } |
| 2039 | } |
| 2040 | |
| 2041 | |
| 2042 | |
| 2043 | //======================================================================= |
| 2044 | //function : This evaluates the lagrange polynomial and its derivatives |
| 2045 | //purpose : up to the requested order that interpolates a series of |
| 2046 | //points of dimension <Dimension> with given assigned parameters |
| 2047 | //======================================================================= |
| 2048 | |
| 2049 | Standard_Integer |
| 2050 | PLib::EvalLagrange(const Standard_Real Parameter, |
| 2051 | const Standard_Integer DerivativeRequest, |
| 2052 | const Standard_Integer Degree, |
| 2053 | const Standard_Integer Dimension, |
| 2054 | Standard_Real& Values, |
| 2055 | Standard_Real& Parameters, |
| 2056 | Standard_Real& Results) |
| 2057 | { |
| 2058 | // |
| 2059 | // the points are assumed to be stored as follows in the Values array : |
| 2060 | // |
| 2061 | // [0] [Dimension -1] first point coefficients |
| 2062 | // |
| 2063 | // [Dimension] [Dimension + Dimension -1] second point coefficients |
| 2064 | // |
| 2065 | // [2 * Dimension] [2 * Dimension + Dimension-1] third point coefficients |
| 2066 | // |
| 2067 | // ................................................... |
| 2068 | // |
| 2069 | // |
| 2070 | // |
| 2071 | // [d * Dimension] [d * Dimension + Dimension-1] d + 1 point coefficients |
| 2072 | // |
| 2073 | // where d is the Degree |
| 2074 | // |
| 2075 | // The ParameterArray stores the parameter value assign to each point in |
| 2076 | // order described above, that is |
| 2077 | // [0] is assign to first point |
| 2078 | // [1] is assign to second point |
| 2079 | // |
| 2080 | Standard_Integer ii, jj, kk, Index, Index1, ReturnCode=0; |
| 2081 | Standard_Integer local_request = DerivativeRequest; |
| 2082 | Standard_Real *ParameterArray; |
| 2083 | Standard_Real difference; |
| 2084 | Standard_Real *PointsArray; |
| 2085 | Standard_Real *ResultArray ; |
| 2086 | |
| 2087 | PointsArray = &Values ; |
| 2088 | ParameterArray = &Parameters ; |
| 2089 | ResultArray = &Results ; |
| 2090 | if (local_request >= Degree) { |
| 2091 | local_request = Degree ; |
| 2092 | } |
| 2093 | NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension); |
| 2094 | // |
| 2095 | // Build the divided differences array |
| 2096 | // |
| 2097 | |
| 2098 | for (ii = 0 ; ii < (Degree + 1) * Dimension ; ii++) { |
| 2099 | divided_differences_array[ii] = PointsArray[ii] ; |
| 2100 | } |
| 2101 | |
| 2102 | for (ii = Degree ; ii >= 0 ; ii--) { |
| 2103 | |
| 2104 | for (jj = Degree ; jj > Degree - ii ; jj--) { |
| 2105 | Index = jj * Dimension ; |
| 2106 | Index1 = Index - Dimension ; |
| 2107 | |
| 2108 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 2109 | divided_differences_array[Index + kk] -= |
| 2110 | divided_differences_array[Index1 + kk] ; |
| 2111 | } |
| 2112 | difference = |
| 2113 | ParameterArray[jj] - ParameterArray[jj - Degree -1 + ii] ; |
| 2114 | if (Abs(difference) < RealSmall()) { |
| 2115 | ReturnCode = 1 ; |
| 2116 | goto FINISH ; |
| 2117 | } |
| 2118 | difference = 1.0e0 / difference ; |
| 2119 | |
| 2120 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 2121 | divided_differences_array[Index + kk] *= difference ; |
| 2122 | } |
| 2123 | } |
| 2124 | } |
| 2125 | // |
| 2126 | // |
| 2127 | // Evaluate the divided difference array polynomial which expresses as |
| 2128 | // |
| 2129 | // P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ... |
| 2130 | // + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P |
| 2131 | // |
| 2132 | // The ith slot in the divided_differences_array is [t1,t2,...,ti+1] |
| 2133 | // |
| 2134 | // |
| 2135 | Index = Degree * Dimension ; |
| 2136 | |
| 2137 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 2138 | ResultArray[kk] = divided_differences_array[Index + kk] ; |
| 2139 | } |
| 2140 | |
| 2141 | for (ii = Dimension ; ii < (local_request + 1) * Dimension ; ii++) { |
| 2142 | ResultArray[ii] = 0.0e0 ; |
| 2143 | } |
| 2144 | |
| 2145 | for (ii = Degree ; ii >= 1 ; ii--) { |
| 2146 | difference = Parameter - ParameterArray[ii - 1] ; |
| 2147 | |
| 2148 | for (jj = local_request ; jj > 0 ; jj--) { |
| 2149 | Index = jj * Dimension ; |
| 2150 | Index1 = Index - Dimension ; |
| 2151 | |
| 2152 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 2153 | ResultArray[Index + kk] *= difference ; |
| 2154 | ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj ; |
| 2155 | } |
| 2156 | } |
| 2157 | Index = (ii -1) * Dimension ; |
| 2158 | |
| 2159 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 2160 | ResultArray[kk] *= difference ; |
| 2161 | ResultArray[kk] += divided_differences_array[Index+kk] ; |
| 2162 | } |
| 2163 | } |
| 2164 | FINISH : |
| 2165 | return (ReturnCode) ; |
| 2166 | } |
| 2167 | |
| 2168 | //======================================================================= |
| 2169 | //function : This evaluates the hermite polynomial and its derivatives |
| 2170 | //purpose : up to the requested order that interpolates a series of |
| 2171 | //points of dimension <Dimension> with given assigned parameters |
| 2172 | //======================================================================= |
| 2173 | |
| 2174 | Standard_Integer PLib::EvalCubicHermite |
| 2175 | (const Standard_Real Parameter, |
| 2176 | const Standard_Integer DerivativeRequest, |
| 2177 | const Standard_Integer Dimension, |
| 2178 | Standard_Real& Values, |
| 2179 | Standard_Real& Derivatives, |
| 2180 | Standard_Real& theParameters, |
| 2181 | Standard_Real& Results) |
| 2182 | { |
| 2183 | // |
| 2184 | // the points are assumed to be stored as follows in the Values array : |
| 2185 | // |
| 2186 | // [0] [Dimension -1] first point coefficients |
| 2187 | // |
| 2188 | // [Dimension] [Dimension + Dimension -1] last point coefficients |
| 2189 | // |
| 2190 | // |
| 2191 | // the derivatives are assumed to be stored as follows in |
| 2192 | // the Derivatives array : |
| 2193 | // |
| 2194 | // [0] [Dimension -1] first point coefficients |
| 2195 | // |
| 2196 | // [Dimension] [Dimension + Dimension -1] last point coefficients |
| 2197 | // |
| 2198 | // The ParameterArray stores the parameter value assign to each point in |
| 2199 | // order described above, that is |
| 2200 | // [0] is assign to first point |
| 2201 | // [1] is assign to last point |
| 2202 | // |
| 2203 | Standard_Integer ii, jj, kk, pp, Index, Index1, Degree, ReturnCode; |
| 2204 | Standard_Integer local_request = DerivativeRequest ; |
| 2205 | |
| 2206 | ReturnCode = 0 ; |
| 2207 | Degree = 3 ; |
| 2208 | Standard_Real ParametersArray[4]; |
| 2209 | Standard_Real difference; |
| 2210 | Standard_Real inverse; |
| 2211 | Standard_Real *FirstLast; |
| 2212 | Standard_Real *PointsArray; |
| 2213 | Standard_Real *DerivativesArray; |
| 2214 | Standard_Real *ResultArray ; |
| 2215 | |
| 2216 | DerivativesArray = &Derivatives ; |
| 2217 | PointsArray = &Values ; |
| 2218 | FirstLast = &theParameters ; |
| 2219 | ResultArray = &Results ; |
| 2220 | if (local_request >= Degree) { |
| 2221 | local_request = Degree ; |
| 2222 | } |
| 2223 | NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension); |
| 2224 | |
| 2225 | for (ii = 0, jj = 0 ; ii < 2 ; ii++, jj+= 2) { |
| 2226 | ParametersArray[jj] = |
| 2227 | ParametersArray[jj+1] = FirstLast[ii] ; |
| 2228 | } |
| 2229 | // |
| 2230 | // Build the divided differences array |
| 2231 | // |
| 2232 | // |
| 2233 | // initialise it at the stage 2 of the building algorithm |
| 2234 | // for devided differences |
| 2235 | // |
| 2236 | inverse = FirstLast[1] - FirstLast[0] ; |
| 2237 | inverse = 1.0e0 / inverse ; |
| 2238 | |
| 2239 | for (ii = 0, jj = Dimension, kk = 2 * Dimension, pp = 3 * Dimension ; |
| 2240 | ii < Dimension ; |
| 2241 | ii++, jj++, kk++, pp++) { |
| 2242 | divided_differences_array[ii] = PointsArray[ii] ; |
| 2243 | divided_differences_array[kk] = inverse * |
| 2244 | (PointsArray[jj] - PointsArray[ii]) ; |
| 2245 | divided_differences_array[jj] = DerivativesArray[ii] ; |
| 2246 | divided_differences_array[pp] = DerivativesArray[jj] ; |
| 2247 | } |
| 2248 | |
| 2249 | for (ii = 1 ; ii <= Degree ; ii++) { |
| 2250 | |
| 2251 | for (jj = Degree ; jj >= ii+1 ; jj--) { |
| 2252 | Index = jj * Dimension ; |
| 2253 | Index1 = Index - Dimension ; |
| 2254 | |
| 2255 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 2256 | divided_differences_array[Index + kk] -= |
| 2257 | divided_differences_array[Index1 + kk] ; |
| 2258 | } |
| 2259 | |
| 2260 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 2261 | divided_differences_array[Index + kk] *= inverse ; |
| 2262 | } |
| 2263 | } |
| 2264 | } |
| 2265 | // |
| 2266 | // |
| 2267 | // Evaluate the divided difference array polynomial which expresses as |
| 2268 | // |
| 2269 | // P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ... |
| 2270 | // + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P |
| 2271 | // |
| 2272 | // The ith slot in the divided_differences_array is [t1,t2,...,ti+1] |
| 2273 | // |
| 2274 | // |
| 2275 | Index = Degree * Dimension ; |
| 2276 | |
| 2277 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 2278 | ResultArray[kk] = divided_differences_array[Index + kk] ; |
| 2279 | } |
| 2280 | |
| 2281 | for (ii = Dimension ; ii < (local_request + 1) * Dimension ; ii++) { |
| 2282 | ResultArray[ii] = 0.0e0 ; |
| 2283 | } |
| 2284 | |
| 2285 | for (ii = Degree ; ii >= 1 ; ii--) { |
| 2286 | difference = Parameter - ParametersArray[ii - 1] ; |
| 2287 | |
| 2288 | for (jj = local_request ; jj > 0 ; jj--) { |
| 2289 | Index = jj * Dimension ; |
| 2290 | Index1 = Index - Dimension ; |
| 2291 | |
| 2292 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 2293 | ResultArray[Index + kk] *= difference ; |
| 2294 | ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj; |
| 2295 | } |
| 2296 | } |
| 2297 | Index = (ii -1) * Dimension ; |
| 2298 | |
| 2299 | for (kk = 0 ; kk < Dimension ; kk++) { |
| 2300 | ResultArray[kk] *= difference ; |
| 2301 | ResultArray[kk] += divided_differences_array[Index+kk] ; |
| 2302 | } |
| 2303 | } |
| 2304 | // FINISH : |
| 2305 | return (ReturnCode) ; |
| 2306 | } |
| 2307 | |
| 2308 | //======================================================================= |
| 2309 | //function : HermiteCoefficients |
| 2310 | //purpose : calcul des polynomes d'Hermite |
| 2311 | //======================================================================= |
| 2312 | |
| 2313 | Standard_Boolean PLib::HermiteCoefficients(const Standard_Real FirstParameter, |
| 2314 | const Standard_Real LastParameter, |
| 2315 | const Standard_Integer FirstOrder, |
| 2316 | const Standard_Integer LastOrder, |
| 2317 | math_Matrix& MatrixCoefs) |
| 2318 | { |
| 2319 | Standard_Integer NbCoeff = FirstOrder + LastOrder + 2, Ordre[2]; |
| 2320 | Standard_Integer ii, jj, pp, cote, iof=0; |
| 2321 | Standard_Real Prod, TBorne = FirstParameter; |
| 2322 | math_Vector Coeff(1,NbCoeff), B(1, NbCoeff, 0.0); |
| 2323 | math_Matrix MAT(1,NbCoeff, 1,NbCoeff, 0.0); |
| 2324 | |
| 2325 | // Test de validites |
| 2326 | |
| 2327 | if ((FirstOrder < 0) || (LastOrder < 0)) return Standard_False; |
| 2328 | Standard_Real D1 = fabs(FirstParameter), D2 = fabs(LastParameter); |
| 2329 | if (D1 > 100 || D2 > 100) return Standard_False; |
| 2330 | D2 += D1; |
| 2331 | if (D2 < 0.01) return Standard_False; |
| 2332 | if (fabs(LastParameter - FirstParameter) / D2 < 0.01) return Standard_False; |
| 2333 | |
| 2334 | // Calcul de la matrice a inverser (MAT) |
| 2335 | |
| 2336 | Ordre[0] = FirstOrder+1; |
| 2337 | Ordre[1] = LastOrder+1; |
| 2338 | |
| 2339 | for (cote=0; cote<=1; cote++) { |
| 2340 | Coeff.Init(1); |
| 2341 | |
| 2342 | for (pp=1; pp<=Ordre[cote]; pp++) { |
| 2343 | ii = pp + iof; |
| 2344 | Prod = 1; |
| 2345 | |
| 2346 | for (jj=pp; jj<=NbCoeff; jj++) { |
| 2347 | // tout se passe dans les 3 lignes suivantes |
| 2348 | MAT(ii, jj) = Coeff(jj) * Prod; |
| 2349 | Coeff(jj) *= jj - pp; |
| 2350 | Prod *= TBorne; |
| 2351 | } |
| 2352 | } |
| 2353 | TBorne = LastParameter; |
| 2354 | iof = Ordre[0]; |
| 2355 | } |
| 2356 | |
| 2357 | // resolution du systemes |
| 2358 | math_Gauss ResolCoeff(MAT, 1.0e-10); |
| 2359 | if (!ResolCoeff.IsDone()) return Standard_False; |
| 2360 | |
| 2361 | for (ii=1; ii<=NbCoeff; ii++) { |
| 2362 | B(ii) = 1; |
| 2363 | ResolCoeff.Solve(B, Coeff); |
| 2364 | MatrixCoefs.SetRow( ii, Coeff); |
| 2365 | B(ii) = 0; |
| 2366 | } |
| 2367 | return Standard_True; |
| 2368 | } |
| 2369 | |
| 2370 | //======================================================================= |
| 2371 | //function : CoefficientsPoles |
| 2372 | //purpose : |
| 2373 | //======================================================================= |
| 2374 | |
| 2375 | void PLib::CoefficientsPoles (const TColgp_Array1OfPnt& Coefs, |
| 2376 | const TColStd_Array1OfReal& WCoefs, |
| 2377 | TColgp_Array1OfPnt& Poles, |
| 2378 | TColStd_Array1OfReal& Weights) |
| 2379 | { |
| 2380 | TColStd_Array1OfReal tempC(1,3*Coefs.Length()); |
| 2381 | PLib::SetPoles(Coefs,tempC); |
| 2382 | TColStd_Array1OfReal tempP(1,3*Poles.Length()); |
| 2383 | PLib::SetPoles(Coefs,tempP); |
| 2384 | PLib::CoefficientsPoles(3,tempC,WCoefs,tempP,Weights); |
| 2385 | PLib::GetPoles(tempP,Poles); |
| 2386 | } |
| 2387 | |
| 2388 | //======================================================================= |
| 2389 | //function : CoefficientsPoles |
| 2390 | //purpose : |
| 2391 | //======================================================================= |
| 2392 | |
| 2393 | void PLib::CoefficientsPoles (const TColgp_Array1OfPnt2d& Coefs, |
| 2394 | const TColStd_Array1OfReal& WCoefs, |
| 2395 | TColgp_Array1OfPnt2d& Poles, |
| 2396 | TColStd_Array1OfReal& Weights) |
| 2397 | { |
| 2398 | TColStd_Array1OfReal tempC(1,2*Coefs.Length()); |
| 2399 | PLib::SetPoles(Coefs,tempC); |
| 2400 | TColStd_Array1OfReal tempP(1,2*Poles.Length()); |
| 2401 | PLib::SetPoles(Coefs,tempP); |
| 2402 | PLib::CoefficientsPoles(2,tempC,WCoefs,tempP,Weights); |
| 2403 | PLib::GetPoles(tempP,Poles); |
| 2404 | } |
| 2405 | |
| 2406 | //======================================================================= |
| 2407 | //function : CoefficientsPoles |
| 2408 | //purpose : |
| 2409 | //======================================================================= |
| 2410 | |
| 2411 | void PLib::CoefficientsPoles (const TColStd_Array1OfReal& Coefs, |
| 2412 | const TColStd_Array1OfReal& WCoefs, |
| 2413 | TColStd_Array1OfReal& Poles, |
| 2414 | TColStd_Array1OfReal& Weights) |
| 2415 | { |
| 2416 | PLib::CoefficientsPoles(1,Coefs,WCoefs,Poles,Weights); |
| 2417 | } |
| 2418 | |
| 2419 | //======================================================================= |
| 2420 | //function : CoefficientsPoles |
| 2421 | //purpose : |
| 2422 | //======================================================================= |
| 2423 | |
| 2424 | void PLib::CoefficientsPoles (const Standard_Integer dim, |
| 2425 | const TColStd_Array1OfReal& Coefs, |
| 2426 | const TColStd_Array1OfReal& WCoefs, |
| 2427 | TColStd_Array1OfReal& Poles, |
| 2428 | TColStd_Array1OfReal& Weights) |
| 2429 | { |
| 2430 | Standard_Boolean rat = &WCoefs != NULL; |
| 2431 | Standard_Integer loc = Coefs.Lower(); |
| 2432 | Standard_Integer lop = Poles.Lower(); |
| 2433 | Standard_Integer lowc=0; |
| 2434 | Standard_Integer lowp=0; |
| 2435 | Standard_Integer upc = Coefs.Upper(); |
| 2436 | Standard_Integer upp = Poles.Upper(); |
| 2437 | Standard_Integer upwc=0; |
| 2438 | Standard_Integer upwp=0; |
| 2439 | Standard_Integer reflen = Coefs.Length()/dim; |
| 2440 | Standard_Integer i,j,k; |
| 2441 | //Les Extremites. |
| 2442 | if (rat) { |
| 2443 | lowc = WCoefs.Lower(); lowp = Weights.Lower(); |
| 2444 | upwc = WCoefs.Upper(); upwp = Weights.Upper(); |
| 2445 | } |
| 2446 | |
| 2447 | for (i = 0; i < dim; i++){ |
| 2448 | Poles (lop + i) = Coefs (loc + i); |
| 2449 | Poles (upp - i) = Coefs (upc - i); |
| 2450 | } |
| 2451 | if (rat) { |
| 2452 | Weights (lowp) = WCoefs (lowc); |
| 2453 | Weights (upwp) = WCoefs (upwc); |
| 2454 | } |
| 2455 | |
| 2456 | Standard_Real Cnp; |
| 2457 | for (i = 2; i < reflen; i++ ) { |
| 2458 | Cnp = PLib::Bin(reflen - 1, i - 1); |
| 2459 | if (rat) Weights (lowp + i - 1) = WCoefs (lowc + i - 1) / Cnp; |
| 2460 | |
| 2461 | for(j = 0; j < dim; j++){ |
| 2462 | Poles(lop + dim * (i-1) + j)= Coefs(loc + dim * (i-1) + j) / Cnp; |
| 2463 | } |
| 2464 | } |
| 2465 | |
| 2466 | for (i = 1; i <= reflen - 1; i++) { |
| 2467 | |
| 2468 | for (j = reflen - 1; j >= i; j--) { |
| 2469 | if (rat) Weights (lowp + j) += Weights (lowp + j -1); |
| 2470 | |
| 2471 | for(k = 0; k < dim; k++){ |
| 2472 | Poles(lop + dim * j + k) += Poles(lop + dim * (j - 1) + k); |
| 2473 | } |
| 2474 | } |
| 2475 | } |
| 2476 | if (rat) { |
| 2477 | |
| 2478 | for (i = 1; i <= reflen; i++) { |
| 2479 | |
| 2480 | for(j = 0; j < dim; j++){ |
| 2481 | Poles(lop + dim * (i-1) + j) /= Weights(lowp + i -1); |
| 2482 | } |
| 2483 | } |
| 2484 | } |
| 2485 | } |
| 2486 | |
| 2487 | //======================================================================= |
| 2488 | //function : Trimming |
| 2489 | //purpose : |
| 2490 | //======================================================================= |
| 2491 | |
| 2492 | void PLib::Trimming(const Standard_Real U1, |
| 2493 | const Standard_Real U2, |
| 2494 | TColgp_Array1OfPnt& Coefs, |
| 2495 | TColStd_Array1OfReal& WCoefs) |
| 2496 | { |
| 2497 | TColStd_Array1OfReal temp(1,3*Coefs.Length()); |
| 2498 | PLib::SetPoles(Coefs,temp); |
| 2499 | PLib::Trimming(U1,U2,3,temp,WCoefs); |
| 2500 | PLib::GetPoles(temp,Coefs); |
| 2501 | } |
| 2502 | |
| 2503 | //======================================================================= |
| 2504 | //function : Trimming |
| 2505 | //purpose : |
| 2506 | //======================================================================= |
| 2507 | |
| 2508 | void PLib::Trimming(const Standard_Real U1, |
| 2509 | const Standard_Real U2, |
| 2510 | TColgp_Array1OfPnt2d& Coefs, |
| 2511 | TColStd_Array1OfReal& WCoefs) |
| 2512 | { |
| 2513 | TColStd_Array1OfReal temp(1,2*Coefs.Length()); |
| 2514 | PLib::SetPoles(Coefs,temp); |
| 2515 | PLib::Trimming(U1,U2,2,temp,WCoefs); |
| 2516 | PLib::GetPoles(temp,Coefs); |
| 2517 | } |
| 2518 | |
| 2519 | //======================================================================= |
| 2520 | //function : Trimming |
| 2521 | //purpose : |
| 2522 | //======================================================================= |
| 2523 | |
| 2524 | void PLib::Trimming(const Standard_Real U1, |
| 2525 | const Standard_Real U2, |
| 2526 | TColStd_Array1OfReal& Coefs, |
| 2527 | TColStd_Array1OfReal& WCoefs) |
| 2528 | { |
| 2529 | PLib::Trimming(U1,U2,1,Coefs,WCoefs); |
| 2530 | } |
| 2531 | |
| 2532 | //======================================================================= |
| 2533 | //function : Trimming |
| 2534 | //purpose : |
| 2535 | //======================================================================= |
| 2536 | |
| 2537 | void PLib::Trimming(const Standard_Real U1, |
| 2538 | const Standard_Real U2, |
| 2539 | const Standard_Integer dim, |
| 2540 | TColStd_Array1OfReal& Coefs, |
| 2541 | TColStd_Array1OfReal& WCoefs) |
| 2542 | { |
| 2543 | |
| 2544 | // principe : |
| 2545 | // on fait le changement de variable v = (u-U1) / (U2-U1) |
| 2546 | // on exprime u = f(v) que l'on remplace dans l'expression polynomiale |
| 2547 | // decomposee sous la forme du schema iteratif de horner. |
| 2548 | |
| 2549 | Standard_Real lsp = U2 - U1; |
| 2550 | Standard_Integer indc, indw=0; |
| 2551 | Standard_Integer upc = Coefs.Upper() - dim + 1, upw=0; |
| 2552 | Standard_Integer len = Coefs.Length()/dim; |
| 2553 | Standard_Boolean rat = &WCoefs != NULL; |
| 2554 | |
| 2555 | if (rat) { |
| 2556 | if(len != WCoefs.Length()) |
| 2557 | Standard_Failure::Raise("PLib::Trimming : nbcoefs/dim != nbweights !!!"); |
| 2558 | upw = WCoefs.Upper(); |
| 2559 | } |
| 2560 | len --; |
| 2561 | |
| 2562 | for (Standard_Integer i = 1; i <= len; i++) { |
| 2563 | Standard_Integer j ; |
| 2564 | indc = upc - dim*(i-1); |
| 2565 | if (rat) indw = upw - i + 1; |
| 2566 | //calcul du coefficient de degre le plus faible a l'iteration i |
| 2567 | |
| 2568 | for( j = 0; j < dim; j++){ |
| 2569 | Coefs(indc - dim + j) += U1 * Coefs(indc + j); |
| 2570 | } |
| 2571 | if (rat) WCoefs(indw - 1) += U1 * WCoefs(indw); |
| 2572 | |
| 2573 | //calcul des coefficients intermediaires : |
| 2574 | |
| 2575 | while (indc < upc){ |
| 2576 | indc += dim; |
| 2577 | |
| 2578 | for(Standard_Integer k = 0; k < dim; k++){ |
| 2579 | Coefs(indc - dim + k) = |
| 2580 | U1 * Coefs(indc + k) + lsp * Coefs(indc - dim + k); |
| 2581 | } |
| 2582 | if (rat) { |
| 2583 | indw ++; |
| 2584 | WCoefs(indw - 1) = U1 * WCoefs(indw) + lsp * WCoefs(indw - 1); |
| 2585 | } |
| 2586 | } |
| 2587 | |
| 2588 | //calcul du coefficient de degre le plus eleve : |
| 2589 | |
| 2590 | for(j = 0; j < dim; j++){ |
| 2591 | Coefs(upc + j) *= lsp; |
| 2592 | } |
| 2593 | if (rat) WCoefs(upw) *= lsp; |
| 2594 | } |
| 2595 | } |
| 2596 | |
| 2597 | //======================================================================= |
| 2598 | //function : CoefficientsPoles |
| 2599 | //purpose : |
| 2600 | // Modified: 21/10/1996 by PMN : PolesCoefficient (PRO5852). |
| 2601 | // on ne bidouille plus les u et v c'est a l'appelant de savoir ce qu'il |
| 2602 | // fait avec BuildCache ou plus simplement d'utiliser PolesCoefficients |
| 2603 | //======================================================================= |
| 2604 | |
| 2605 | void PLib::CoefficientsPoles (const TColgp_Array2OfPnt& Coefs, |
| 2606 | const TColStd_Array2OfReal& WCoefs, |
| 2607 | TColgp_Array2OfPnt& Poles, |
| 2608 | TColStd_Array2OfReal& Weights) |
| 2609 | { |
| 2610 | Standard_Boolean rat = (&WCoefs != NULL); |
| 2611 | Standard_Integer LowerRow = Poles.LowerRow(); |
| 2612 | Standard_Integer UpperRow = Poles.UpperRow(); |
| 2613 | Standard_Integer LowerCol = Poles.LowerCol(); |
| 2614 | Standard_Integer UpperCol = Poles.UpperCol(); |
| 2615 | Standard_Integer ColLength = Poles.ColLength(); |
| 2616 | Standard_Integer RowLength = Poles.RowLength(); |
| 2617 | |
| 2618 | // Bidouille pour retablir u et v pour les coefs calcules |
| 2619 | // par buildcache |
| 2620 | // Standard_Boolean inv = Standard_False; //ColLength != Coefs.ColLength(); |
| 2621 | |
| 2622 | Standard_Integer Row, Col; |
| 2623 | Standard_Real W, Cnp; |
| 2624 | |
| 2625 | Standard_Integer I1, I2; |
| 2626 | Standard_Integer NPoleu , NPolev; |
| 2627 | gp_XYZ Temp; |
| 2628 | |
| 2629 | for (NPoleu = LowerRow; NPoleu <= UpperRow; NPoleu++){ |
| 2630 | Poles (NPoleu, LowerCol) = Coefs (NPoleu, LowerCol); |
| 2631 | if (rat) { |
| 2632 | Weights (NPoleu, LowerCol) = WCoefs (NPoleu, LowerCol); |
| 2633 | } |
| 2634 | |
| 2635 | for (Col = LowerCol + 1; Col <= UpperCol - 1; Col++) { |
| 2636 | Cnp = PLib::Bin(RowLength - 1,Col - LowerCol); |
| 2637 | Temp = Coefs (NPoleu, Col).XYZ(); |
| 2638 | Temp.Divide (Cnp); |
| 2639 | Poles (NPoleu, Col).SetXYZ (Temp); |
| 2640 | if (rat) { |
| 2641 | Weights (NPoleu, Col) = WCoefs (NPoleu, Col) / Cnp; |
| 2642 | } |
| 2643 | } |
| 2644 | Poles (NPoleu, UpperCol) = Coefs (NPoleu, UpperCol); |
| 2645 | if (rat) { |
| 2646 | Weights (NPoleu, UpperCol) = WCoefs (NPoleu, UpperCol); |
| 2647 | } |
| 2648 | |
| 2649 | for (I1 = 1; I1 <= RowLength - 1; I1++) { |
| 2650 | |
| 2651 | for (I2 = UpperCol; I2 >= LowerCol + I1; I2--) { |
| 2652 | Temp.SetLinearForm |
| 2653 | (Poles (NPoleu, I2).XYZ(), Poles (NPoleu, I2-1).XYZ()); |
| 2654 | Poles (NPoleu, I2).SetXYZ (Temp); |
| 2655 | if (rat) Weights(NPoleu, I2) += Weights(NPoleu, I2-1); |
| 2656 | } |
| 2657 | } |
| 2658 | } |
| 2659 | |
| 2660 | for (NPolev = LowerCol; NPolev <= UpperCol; NPolev++){ |
| 2661 | |
| 2662 | for (Row = LowerRow + 1; Row <= UpperRow - 1; Row++) { |
| 2663 | Cnp = PLib::Bin(ColLength - 1,Row - LowerRow); |
| 2664 | Temp = Poles (Row, NPolev).XYZ(); |
| 2665 | Temp.Divide (Cnp); |
| 2666 | Poles (Row, NPolev).SetXYZ (Temp); |
| 2667 | if (rat) Weights(Row, NPolev) /= Cnp; |
| 2668 | } |
| 2669 | |
| 2670 | for (I1 = 1; I1 <= ColLength - 1; I1++) { |
| 2671 | |
| 2672 | for (I2 = UpperRow; I2 >= LowerRow + I1; I2--) { |
| 2673 | Temp.SetLinearForm |
| 2674 | (Poles (I2, NPolev).XYZ(), Poles (I2-1, NPolev).XYZ()); |
| 2675 | Poles (I2, NPolev).SetXYZ (Temp); |
| 2676 | if (rat) Weights(I2, NPolev) += Weights(I2-1, NPolev); |
| 2677 | } |
| 2678 | } |
| 2679 | } |
| 2680 | if (rat) { |
| 2681 | |
| 2682 | for (Row = LowerRow; Row <= UpperRow; Row++) { |
| 2683 | |
| 2684 | for (Col = LowerCol; Col <= UpperCol; Col++) { |
| 2685 | W = Weights (Row, Col); |
| 2686 | Temp = Poles(Row, Col).XYZ(); |
| 2687 | Temp.Divide (W); |
| 2688 | Poles(Row, Col).SetXYZ (Temp); |
| 2689 | } |
| 2690 | } |
| 2691 | } |
| 2692 | } |
| 2693 | |
| 2694 | //======================================================================= |
| 2695 | //function : UTrimming |
| 2696 | //purpose : |
| 2697 | //======================================================================= |
| 2698 | |
| 2699 | void PLib::UTrimming(const Standard_Real U1, |
| 2700 | const Standard_Real U2, |
| 2701 | TColgp_Array2OfPnt& Coeffs, |
| 2702 | TColStd_Array2OfReal& WCoeffs) |
| 2703 | { |
| 2704 | Standard_Boolean rat = &WCoeffs != NULL; |
| 2705 | Standard_Integer lr = Coeffs.LowerRow(); |
| 2706 | Standard_Integer ur = Coeffs.UpperRow(); |
| 2707 | Standard_Integer lc = Coeffs.LowerCol(); |
| 2708 | Standard_Integer uc = Coeffs.UpperCol(); |
| 2709 | TColgp_Array1OfPnt Temp (lr,ur); |
| 2710 | TColStd_Array1OfReal Temw (lr,ur); |
| 2711 | |
| 2712 | for (Standard_Integer icol = lc; icol <= uc; icol++) { |
| 2713 | Standard_Integer irow ; |
| 2714 | for ( irow = lr; irow <= ur; irow++) { |
| 2715 | Temp (irow) = Coeffs (irow, icol); |
| 2716 | if (rat) Temw (irow) = WCoeffs (irow, icol); |
| 2717 | } |
| 2718 | if (rat) PLib::Trimming (U1, U2, Temp, Temw); |
| 2719 | else PLib::Trimming (U1, U2, Temp, PLib::NoWeights()); |
| 2720 | |
| 2721 | for (irow = lr; irow <= ur; irow++) { |
| 2722 | Coeffs (irow, icol) = Temp (irow); |
| 2723 | if (rat) WCoeffs (irow, icol) = Temw (irow); |
| 2724 | } |
| 2725 | } |
| 2726 | } |
| 2727 | |
| 2728 | //======================================================================= |
| 2729 | //function : VTrimming |
| 2730 | //purpose : |
| 2731 | //======================================================================= |
| 2732 | |
| 2733 | void PLib::VTrimming(const Standard_Real V1, |
| 2734 | const Standard_Real V2, |
| 2735 | TColgp_Array2OfPnt& Coeffs, |
| 2736 | TColStd_Array2OfReal& WCoeffs) |
| 2737 | { |
| 2738 | Standard_Boolean rat = &WCoeffs != NULL; |
| 2739 | Standard_Integer lr = Coeffs.LowerRow(); |
| 2740 | Standard_Integer ur = Coeffs.UpperRow(); |
| 2741 | Standard_Integer lc = Coeffs.LowerCol(); |
| 2742 | Standard_Integer uc = Coeffs.UpperCol(); |
| 2743 | TColgp_Array1OfPnt Temp (lc,uc); |
| 2744 | TColStd_Array1OfReal Temw (lc,uc); |
| 2745 | |
| 2746 | for (Standard_Integer irow = lr; irow <= ur; irow++) { |
| 2747 | Standard_Integer icol ; |
| 2748 | for ( icol = lc; icol <= uc; icol++) { |
| 2749 | Temp (icol) = Coeffs (irow, icol); |
| 2750 | if (rat) Temw (icol) = WCoeffs (irow, icol); |
| 2751 | } |
| 2752 | if (rat) PLib::Trimming (V1, V2, Temp, Temw); |
| 2753 | else PLib::Trimming (V1, V2, Temp, PLib::NoWeights()); |
| 2754 | |
| 2755 | for (icol = lc; icol <= uc; icol++) { |
| 2756 | Coeffs (irow, icol) = Temp (icol); |
| 2757 | if (rat) WCoeffs (irow, icol) = Temw (icol); |
| 2758 | } |
| 2759 | } |
| 2760 | } |
| 2761 | |
| 2762 | //======================================================================= |
| 2763 | //function : HermiteInterpolate |
| 2764 | //purpose : |
| 2765 | //======================================================================= |
| 2766 | |
| 2767 | Standard_Boolean PLib::HermiteInterpolate |
| 2768 | (const Standard_Integer Dimension, |
| 2769 | const Standard_Real FirstParameter, |
| 2770 | const Standard_Real LastParameter, |
| 2771 | const Standard_Integer FirstOrder, |
| 2772 | const Standard_Integer LastOrder, |
| 2773 | const TColStd_Array2OfReal& FirstConstr, |
| 2774 | const TColStd_Array2OfReal& LastConstr, |
| 2775 | TColStd_Array1OfReal& Coefficients) |
| 2776 | { |
| 2777 | Standard_Real Pattern[3][6]; |
| 2778 | |
| 2779 | // portage HP : il faut les initialiser 1 par 1 |
| 2780 | |
| 2781 | Pattern[0][0] = 1; |
| 2782 | Pattern[0][1] = 1; |
| 2783 | Pattern[0][2] = 1; |
| 2784 | Pattern[0][3] = 1; |
| 2785 | Pattern[0][4] = 1; |
| 2786 | Pattern[0][5] = 1; |
| 2787 | Pattern[1][0] = 0; |
| 2788 | Pattern[1][1] = 1; |
| 2789 | Pattern[1][2] = 2; |
| 2790 | Pattern[1][3] = 3; |
| 2791 | Pattern[1][4] = 4; |
| 2792 | Pattern[1][5] = 5; |
| 2793 | Pattern[2][0] = 0; |
| 2794 | Pattern[2][1] = 0; |
| 2795 | Pattern[2][2] = 2; |
| 2796 | Pattern[2][3] = 6; |
| 2797 | Pattern[2][4] = 12; |
| 2798 | Pattern[2][5] = 20; |
| 2799 | |
| 2800 | math_Matrix A(0,FirstOrder+LastOrder+1, 0,FirstOrder+LastOrder+1); |
| 2801 | // The initialisation of the matrix A |
| 2802 | Standard_Integer irow ; |
| 2803 | for ( irow=0; irow<=FirstOrder; irow++) { |
| 2804 | Standard_Real FirstVal = 1.; |
| 2805 | |
| 2806 | for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) { |
| 2807 | A(irow,icol) = Pattern[irow][icol]*FirstVal; |
| 2808 | if (irow <= icol) FirstVal *= FirstParameter; |
| 2809 | } |
| 2810 | } |
| 2811 | |
| 2812 | for (irow=0; irow<=LastOrder; irow++) { |
| 2813 | Standard_Real LastVal = 1.; |
| 2814 | |
| 2815 | for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) { |
| 2816 | A(irow+FirstOrder+1,icol) = Pattern[irow][icol]*LastVal; |
| 2817 | if (irow <= icol) LastVal *= LastParameter; |
| 2818 | } |
| 2819 | } |
| 2820 | // |
| 2821 | // The filled matrix A for FirstOrder=LastOrder=2 is: |
| 2822 | // |
| 2823 | // 1 FP FP**2 FP**3 FP**4 FP**5 |
| 2824 | // 0 1 2*FP 3*FP**2 4*FP**3 5*FP**4 FP - FirstParameter |
| 2825 | // 0 0 2 6*FP 12*FP**2 20*FP**3 |
| 2826 | // 1 LP LP**2 LP**3 LP**4 LP**5 |
| 2827 | // 0 1 2*LP 3*LP**2 4*LP**3 5*LP**4 LP - LastParameter |
| 2828 | // 0 0 2 6*LP 12*LP**2 20*LP**3 |
| 2829 | // |
| 2830 | // If FirstOrder or LastOrder <=2 then some rows and columns are missing. |
| 2831 | // For example: |
| 2832 | // If FirstOrder=1 then 3th row and 6th column are missing |
| 2833 | // If FirstOrder=LastOrder=0 then 2,3,5,6th rows and 3,4,5,6th columns are missing |
| 2834 | |
| 2835 | math_Gauss Equations(A); |
| 2836 | // cout << "A=" << A << endl; |
| 2837 | |
| 2838 | for (Standard_Integer idim=1; idim<=Dimension; idim++) { |
| 2839 | // cout << "idim=" << idim << endl; |
| 2840 | |
| 2841 | math_Vector B(0,FirstOrder+LastOrder+1); |
| 2842 | Standard_Integer icol ; |
| 2843 | for ( icol=0; icol<=FirstOrder; icol++) |
| 2844 | B(icol) = FirstConstr(idim,icol); |
| 2845 | |
| 2846 | for (icol=0; icol<=LastOrder; icol++) |
| 2847 | B(FirstOrder+1+icol) = LastConstr(idim,icol); |
| 2848 | // cout << "B=" << B << endl; |
| 2849 | |
| 2850 | // The solving of equations system A * X = B. Then B = X |
| 2851 | Equations.Solve(B); |
| 2852 | // cout << "After Solving" << endl << "B=" << B << endl; |
| 2853 | |
| 2854 | if (Equations.IsDone()==Standard_False) return Standard_False; |
| 2855 | |
| 2856 | // the filling of the Coefficients |
| 2857 | |
| 2858 | for (icol=0; icol<=FirstOrder+LastOrder+1; icol++) |
| 2859 | Coefficients(Dimension*icol+idim-1) = B(icol); |
| 2860 | } |
| 2861 | return Standard_True; |
| 2862 | } |
| 2863 | |
| 2864 | //======================================================================= |
| 2865 | //function : JacobiParameters |
| 2866 | //purpose : |
| 2867 | //======================================================================= |
| 2868 | |
| 2869 | void PLib::JacobiParameters(const GeomAbs_Shape ConstraintOrder, |
| 2870 | const Standard_Integer MaxDegree, |
| 2871 | const Standard_Integer Code, |
| 2872 | Standard_Integer& NbGaussPoints, |
| 2873 | Standard_Integer& WorkDegree) |
| 2874 | { |
| 2875 | // ConstraintOrder: Ordre de contrainte aux extremites : |
| 2876 | // C0 = contraintes de passage aux bornes; |
| 2877 | // C1 = C0 + contraintes de derivees 1eres; |
| 2878 | // C2 = C1 + contraintes de derivees 2ndes. |
| 2879 | // MaxDegree: Nombre maxi de coeff de la "courbe" polynomiale |
| 2880 | // d' approximation (doit etre superieur ou egal a |
| 2881 | // 2*NivConstr+2 et inferieur ou egal a 50). |
| 2882 | // Code: Code d' init. des parametres de discretisation. |
| 2883 | // (choix de NBPNTS et de NDGJAC de MAPF1C,MAPFXC). |
| 2884 | // = -5 Calcul tres rapide mais peu precis (8pts) |
| 2885 | // = -4 ' ' ' ' ' ' (10pts) |
| 2886 | // = -3 ' ' ' ' ' ' (15pts) |
| 2887 | // = -2 ' ' ' ' ' ' (20pts) |
| 2888 | // = -1 ' ' ' ' ' ' (25pts) |
| 2889 | // = 1 calcul rapide avec precision moyenne (30pts). |
| 2890 | // = 2 calcul rapide avec meilleure precision (40pts). |
| 2891 | // = 3 calcul un peu plus lent avec bonne precision (50 pts). |
| 2892 | // = 4 calcul lent avec la meilleure precision possible |
| 2893 | // (61pts). |
| 2894 | |
| 2895 | // The possible values of NbGaussPoints |
| 2896 | |
| 2897 | const Standard_Integer NDEG8=8, NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25, |
| 2898 | NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61; |
| 2899 | |
| 2900 | Standard_Integer NivConstr=0; |
| 2901 | switch (ConstraintOrder) { |
| 2902 | case GeomAbs_C0: NivConstr = 0; break; |
| 2903 | case GeomAbs_C1: NivConstr = 1; break; |
| 2904 | case GeomAbs_C2: NivConstr = 2; break; |
| 2905 | default: |
| 2906 | Standard_ConstructionError::Raise("Invalid ConstraintOrder"); |
| 2907 | } |
| 2908 | if (MaxDegree < 2*NivConstr+1) |
| 2909 | Standard_ConstructionError::Raise("Invalid MaxDegree"); |
| 2910 | |
| 2911 | if (Code >= 1) |
| 2912 | WorkDegree = MaxDegree + 9; |
| 2913 | else |
| 2914 | WorkDegree = MaxDegree + 6; |
| 2915 | |
| 2916 | //---> Nbre mini de points necessaires. |
| 2917 | Standard_Integer IPMIN=0; |
| 2918 | if (WorkDegree < NDEG8) |
| 2919 | IPMIN=NDEG8; |
| 2920 | else if (WorkDegree < NDEG10) |
| 2921 | IPMIN=NDEG10; |
| 2922 | else if (WorkDegree < NDEG15) |
| 2923 | IPMIN=NDEG15; |
| 2924 | else if (WorkDegree < NDEG20) |
| 2925 | IPMIN=NDEG20; |
| 2926 | else if (WorkDegree < NDEG25) |
| 2927 | IPMIN=NDEG25; |
| 2928 | else if (WorkDegree < NDEG30) |
| 2929 | IPMIN=NDEG30; |
| 2930 | else if (WorkDegree < NDEG40) |
| 2931 | IPMIN=NDEG40; |
| 2932 | else if (WorkDegree < NDEG50) |
| 2933 | IPMIN=NDEG50; |
| 2934 | else if (WorkDegree < NDEG61) |
| 2935 | IPMIN=NDEG61; |
| 2936 | else |
| 2937 | Standard_ConstructionError::Raise("Invalid MaxDegree"); |
| 2938 | // ---> Nbre de points voulus. |
| 2939 | Standard_Integer IWANT=0; |
| 2940 | switch (Code) { |
| 2941 | case -5: IWANT=NDEG8; break; |
| 2942 | case -4: IWANT=NDEG10; break; |
| 2943 | case -3: IWANT=NDEG15; break; |
| 2944 | case -2: IWANT=NDEG20; break; |
| 2945 | case -1: IWANT=NDEG25; break; |
| 2946 | case 1: IWANT=NDEG30; break; |
| 2947 | case 2: IWANT=NDEG40; break; |
| 2948 | case 3: IWANT=NDEG50; break; |
| 2949 | case 4: IWANT=NDEG61; break; |
| 2950 | default: |
| 2951 | Standard_ConstructionError::Raise("Invalid Code"); |
| 2952 | } |
| 2953 | //--> NbGaussPoints est le nombre de points de discretisation de la fonction, |
| 2954 | // il ne peut prendre que les valeurs 8,10,15,20,25,30,40,50 ou 61. |
| 2955 | // NbGaussPoints doit etre superieur strictement a WorkDegree. |
| 2956 | NbGaussPoints = Max(IPMIN,IWANT); |
| 2957 | // NbGaussPoints +=2; |
| 2958 | } |
| 2959 | |
| 2960 | //======================================================================= |
| 2961 | //function : NivConstr |
| 2962 | //purpose : translates from GeomAbs_Shape to Integer |
| 2963 | //======================================================================= |
| 2964 | |
| 2965 | Standard_Integer PLib::NivConstr(const GeomAbs_Shape ConstraintOrder) |
| 2966 | { |
| 2967 | Standard_Integer NivConstr=0; |
| 2968 | switch (ConstraintOrder) { |
| 2969 | case GeomAbs_C0: NivConstr = 0; break; |
| 2970 | case GeomAbs_C1: NivConstr = 1; break; |
| 2971 | case GeomAbs_C2: NivConstr = 2; break; |
| 2972 | default: |
| 2973 | Standard_ConstructionError::Raise("Invalid ConstraintOrder"); |
| 2974 | } |
| 2975 | return NivConstr; |
| 2976 | } |
| 2977 | |
| 2978 | //======================================================================= |
| 2979 | //function : ConstraintOrder |
| 2980 | //purpose : translates from Integer to GeomAbs_Shape |
| 2981 | //======================================================================= |
| 2982 | |
| 2983 | GeomAbs_Shape PLib::ConstraintOrder(const Standard_Integer NivConstr) |
| 2984 | { |
| 2985 | GeomAbs_Shape ConstraintOrder=GeomAbs_C0; |
| 2986 | switch (NivConstr) { |
| 2987 | case 0: ConstraintOrder = GeomAbs_C0; break; |
| 2988 | case 1: ConstraintOrder = GeomAbs_C1; break; |
| 2989 | case 2: ConstraintOrder = GeomAbs_C2; break; |
| 2990 | default: |
| 2991 | Standard_ConstructionError::Raise("Invalid NivConstr"); |
| 2992 | } |
| 2993 | return ConstraintOrder; |
| 2994 | } |
| 2995 | |
| 2996 | //======================================================================= |
| 2997 | //function : EvalLength |
| 2998 | //purpose : |
| 2999 | //======================================================================= |
| 3000 | |
| 3001 | void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension, |
| 3002 | Standard_Real& PolynomialCoeff, |
| 3003 | const Standard_Real U1, const Standard_Real U2, |
| 3004 | Standard_Real& Length) |
| 3005 | { |
| 3006 | Standard_Integer i,j,idim, degdim; |
| 3007 | Standard_Real C1,C2,Sum,Tran,X1,X2,Der1,Der2,D1,D2,DD; |
| 3008 | |
| 3009 | Standard_Real *PolynomialArray = &PolynomialCoeff ; |
| 3010 | |
| 3011 | Standard_Integer NbGaussPoints = 4 * Min((Degree/4)+1,10); |
| 3012 | |
| 3013 | math_Vector GaussPoints(1,NbGaussPoints); |
| 3014 | math::GaussPoints(NbGaussPoints,GaussPoints); |
| 3015 | |
| 3016 | math_Vector GaussWeights(1,NbGaussPoints); |
| 3017 | math::GaussWeights(NbGaussPoints,GaussWeights); |
| 3018 | |
| 3019 | C1 = (U2 + U1) / 2.; |
| 3020 | C2 = (U2 - U1) / 2.; |
| 3021 | |
| 3022 | //----------------------------------------------------------- |
| 3023 | //****** Integration - Boucle sur les intervalles de GAUSS ** |
| 3024 | //----------------------------------------------------------- |
| 3025 | |
| 3026 | Sum = 0; |
| 3027 | |
| 3028 | for (j=1; j<=NbGaussPoints/2; j++) { |
| 3029 | // Integration en tenant compte de la symetrie |
| 3030 | Tran = C2 * GaussPoints(j); |
| 3031 | X1 = C1 + Tran; |
| 3032 | X2 = C1 - Tran; |
| 3033 | |
| 3034 | //****** Derivation sur la dimension de l'espace ** |
| 3035 | |
| 3036 | degdim = Degree*Dimension; |
| 3037 | Der1 = Der2 = 0.; |
| 3038 | for (idim=0; idim<Dimension; idim++) { |
| 3039 | D1 = D2 = Degree * PolynomialArray [idim + degdim]; |
| 3040 | for (i=Degree-1; i>=1; i--) { |
| 3041 | DD = i * PolynomialArray [idim + i*Dimension]; |
| 3042 | D1 = D1 * X1 + DD; |
| 3043 | D2 = D2 * X2 + DD; |
| 3044 | } |
| 3045 | Der1 += D1 * D1; |
| 3046 | Der2 += D2 * D2; |
| 3047 | } |
| 3048 | |
| 3049 | //****** Integration ** |
| 3050 | |
| 3051 | Sum += GaussWeights(j) * C2 * (Sqrt(Der1) + Sqrt(Der2)); |
| 3052 | |
| 3053 | //****** Fin de boucle dur les intervalles de GAUSS ** |
| 3054 | } |
| 3055 | Length = Sum; |
| 3056 | } |
| 3057 | |
| 3058 | |
| 3059 | //======================================================================= |
| 3060 | //function : EvalLength |
| 3061 | //purpose : |
| 3062 | //======================================================================= |
| 3063 | |
| 3064 | void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension, |
| 3065 | Standard_Real& PolynomialCoeff, |
| 3066 | const Standard_Real U1, const Standard_Real U2, |
| 3067 | const Standard_Real Tol, |
| 3068 | Standard_Real& Length, Standard_Real& Error) |
| 3069 | { |
| 3070 | Standard_Integer i; |
| 3071 | Standard_Integer NbSubInt = 1, // Current number of subintervals |
| 3072 | MaxNbIter = 13, // Max number of iterations |
| 3073 | NbIter = 1; // Current number of iterations |
| 3074 | Standard_Real dU,OldLen,LenI; |
| 3075 | |
| 3076 | PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1,U2, Length); |
| 3077 | |
| 3078 | do { |
| 3079 | OldLen = Length; |
| 3080 | Length = 0.; |
| 3081 | NbSubInt *= 2; |
| 3082 | dU = (U2-U1)/NbSubInt; |
| 3083 | for (i=1; i<=NbSubInt; i++) { |
| 3084 | PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1+(i-1)*dU,U1+i*dU, LenI); |
| 3085 | Length += LenI; |
| 3086 | } |
| 3087 | NbIter++; |
| 3088 | Error = Abs(OldLen-Length); |
| 3089 | } |
| 3090 | while (Error > Tol && NbIter <= MaxNbIter); |
| 3091 | } |