1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
18 // Implementation de la classe resolvant les equations en cosinus-sinus.
19 // Equation de la forme a*cos(x)*cos(x)+2*b*cos(x)*sin(x)+c*cos(x)+d*sin(x)+e
22 #define No_Standard_RangeError
23 #define No_Standard_OutOfRange
24 #define No_Standard_DimensionError
27 #include <math_TrigonometricFunctionRoots.hxx>
28 #include <math_DirectPolynomialRoots.hxx>
29 #include <Standard_OutOfRange.hxx>
30 #include <math_FunctionWithDerivative.hxx>
31 #include <math_NewtonFunctionRoot.hxx>
32 #include <Precision.hxx>
35 class MyTrigoFunction: public math_FunctionWithDerivative
44 MyTrigoFunction(const Standard_Real A,
45 const Standard_Real B,
46 const Standard_Real C,
47 const Standard_Real D,
48 const Standard_Real E)
49 : AA(A), BB(B), CC(C), DD(D), EE(E)
53 Standard_Boolean Value(const Standard_Real X, Standard_Real& F);
54 Standard_Boolean Derivative(const Standard_Real X, Standard_Real& D);
55 Standard_Boolean Values(const Standard_Real X, Standard_Real& F, Standard_Real& D);
58 Standard_Boolean MyTrigoFunction::Value(const Standard_Real X, Standard_Real& F) {
59 Standard_Real CN= cos(X), SN = sin(X);
60 //-- F= AA*CN*CN+2*BB*CN*SN+CC*CN+DD*SN+EE;
61 F=CN*(AA*CN + (BB+BB)*SN + CC) + DD*SN + EE;
65 Standard_Boolean MyTrigoFunction::Derivative(const Standard_Real X, Standard_Real& D) {
66 Standard_Real CN= Cos(X), SN = Sin(X);
67 //-- D = -2*AA*CN*SN+2*BB*(CN*CN-SN*SN)-CC*SN+DD*CN;
68 D = -AA*CN*SN + BB*(CN*CN-SN*SN);
74 Standard_Boolean MyTrigoFunction::Values(const Standard_Real X, Standard_Real& F, Standard_Real& D) {
75 Standard_Real CN= Cos(X), SN = Sin(X);
76 //-- F= AA*CN*CN+2*BB*CN*SN+CC*CN+DD*SN+EE;
77 //-- D = -2*AA*CN*SN+2*BB*(CN*CN-SN*SN)-CC*SN+DD*CN;
78 Standard_Real AACN = AA*CN;
79 Standard_Real BBSN = BB*SN;
81 F = AACN*CN + BBSN*(CN+CN) + CC*CN + DD*SN + EE;
82 D = -AACN*SN + BB*(CN*CN+SN*SN);
89 math_TrigonometricFunctionRoots::math_TrigonometricFunctionRoots
90 (const Standard_Real theD,
91 const Standard_Real theE,
92 const Standard_Real theInfBound,
93 const Standard_Real theSupBound)
96 InfiniteStatus(Standard_False),
99 const Standard_Real A(0.0), B(0.0), C(0.0);
100 Perform(A, B, C, theD, theE, theInfBound, theSupBound);
104 math_TrigonometricFunctionRoots::math_TrigonometricFunctionRoots
105 (const Standard_Real theC,
106 const Standard_Real theD,
107 const Standard_Real theE,
108 const Standard_Real theInfBound,
109 const Standard_Real theSupBound)
112 InfiniteStatus(Standard_False),
113 Done (Standard_False)
115 const Standard_Real A(0.0), B(0.0);
116 Perform(A, B, theC, theD, theE, theInfBound, theSupBound);
121 math_TrigonometricFunctionRoots::math_TrigonometricFunctionRoots
122 (const Standard_Real theA,
123 const Standard_Real theB,
124 const Standard_Real theC,
125 const Standard_Real theD,
126 const Standard_Real theE,
127 const Standard_Real theInfBound,
128 const Standard_Real theSupBound)
131 InfiniteStatus(Standard_False),
132 Done (Standard_False)
134 Perform(theA, theB, theC, theD, theE, theInfBound, theSupBound);
137 void math_TrigonometricFunctionRoots::Perform(const Standard_Real A,
138 const Standard_Real B,
139 const Standard_Real C,
140 const Standard_Real D,
141 const Standard_Real E,
142 const Standard_Real InfBound,
143 const Standard_Real SupBound) {
145 Standard_Integer i, j=0, k, l, NZer=0, Nit = 10;
146 Standard_Real Depi, Delta, Mod, AA, BB, CC, MyBorneInf;
147 Standard_Real Teta, X;
148 Standard_Real Eps, Tol1 = 1.e-15;
149 TColStd_Array1OfReal ko(1,5), Zer(1,4);
150 Standard_Boolean Flag4;
151 InfiniteStatus = Standard_False;
152 Done = Standard_True;
157 if (InfBound <= RealFirst() && SupBound >= RealLast()) {
162 else if (SupBound >= RealLast()) {
163 MyBorneInf = InfBound;
165 Mod = MyBorneInf/Depi;
167 else if (InfBound <= RealFirst()) {
168 MyBorneInf = SupBound - Depi;
170 Mod = MyBorneInf/Depi;
173 MyBorneInf = InfBound;
174 Delta = SupBound-InfBound;
176 if ((SupBound-InfBound) > Depi) { Delta = Depi;}
179 if ((Abs(A) <= Eps) && (Abs(B) <= Eps)) {
183 InfiniteStatus = Standard_True; // infinite de solutions.
192 // Equation du type d*sin(x) + e = 0
193 // =================================
201 Zer(2) = M_PI - Zer(1);
203 for (i = 1; i <= NZer; i++) {
204 if (Zer(i) <= -Eps) {
205 Zer(i) = Depi - Abs(Zer(i));
207 // On rend les solutions entre InfBound et SupBound:
208 // =================================================
209 Zer(i) += IntegerPart(Mod)*Depi;
210 X = Zer(i)-MyBorneInf;
211 if ((X > (-Epsilon(Delta))) && (X < Delta+ Epsilon(Delta))) {
219 else if (Abs(D) <= Eps) {
221 // Equation du premier degre de la forme c*cos(x) + e = 0
222 // ======================================================
232 for (i = 1; i <= NZer; i++) {
233 if (Zer(i) <= -Eps) {
234 Zer(i) = Depi - Abs(Zer(i));
236 // On rend les solutions entre InfBound et SupBound:
237 // =================================================
238 Zer(i) += IntegerPart(Mod)*2.*M_PI;
239 X = Zer(i)-MyBorneInf;
240 if ((X >= (-Epsilon(Delta))) && (X <= Delta+ Epsilon(Delta))) {
249 // Equation du second degre:
250 // =========================
255 math_DirectPolynomialRoots Resol(AA, BB, CC);
256 if (!Resol.IsDone()) {
257 Done = Standard_False;
260 else if(!Resol.InfiniteRoots()) {
261 NZer = Resol.NbSolutions();
262 for (i = 1; i <= NZer; i++) {
263 Zer(i) = Resol.Value(i);
266 else if (Resol.InfiniteRoots()) {
267 InfiniteStatus = Standard_True;
273 // Two additional analytical cases.
274 if ((Abs(A) <= Eps) &&
279 // 2 * B * sin * cos + D * sin = 0
285 if (Abs(AA) <= 1.0 + Precision::PConfusion())
301 Zer(4) = Depi - Zer(3);
306 for (i = 1; i <= NZer; i++)
308 if (Zer(i) <= MyBorneInf - Eps)
312 // On rend les solutions entre InfBound et SupBound:
313 // =================================================
314 Zer(i) += IntegerPart(Mod)*2.*M_PI;
315 X = Zer(i)-MyBorneInf;
316 if ((X >= (-Precision::PConfusion())) &&
317 (X <= Delta + Precision::PConfusion()))
319 if (Zer(i) < InfBound)
321 if (Zer(i) > SupBound)
331 // 2 * B * sin * cos + C * cos = 0
334 Zer(2) = M_PI * 3.0 / 2.0;
337 if (Abs(AA) <= 1.0 + Precision::PConfusion())
348 Zer(3) = M_PI * 3.0 / 2.0;
349 Zer(4) = M_PI * 3.0 / 2.0;
354 Zer(4) = M_PI - Zer(3);
359 for (i = 1; i <= NZer; i++)
361 if (Zer(i) <= MyBorneInf - Eps)
365 // On rend les solutions entre InfBound et SupBound:
366 // =================================================
367 Zer(i) += IntegerPart(Mod)*2.*M_PI;
368 X = Zer(i)-MyBorneInf;
369 if ((X >= (-Precision::PConfusion())) &&
370 (X <= Delta + Precision::PConfusion()))
372 if (Zer(i) < InfBound)
374 if (Zer(i) > SupBound)
384 // Equation du 4 ieme degre
385 // ========================
391 Standard_Boolean bko;
392 Standard_Integer nbko=0;
395 math_DirectPolynomialRoots Resol4(ko(1), ko(2), ko(3), ko(4), ko(5));
396 if (!Resol4.IsDone()) {
397 Done = Standard_False;
400 else if (!Resol4.InfiniteRoots()) {
401 NZer = Resol4.NbSolutions();
402 for (i = 1; i <= NZer; i++) {
403 Zer(i) = Resol4.Value(i);
406 else if (Resol4.InfiniteRoots()) {
407 InfiniteStatus = Standard_True;
410 Standard_Boolean triok;
413 for(i=1;i<NZer;i++) {
414 if(Zer(i)>Zer(i+1)) {
415 Standard_Real t=Zer(i);
418 triok=Standard_False;
422 while(triok==Standard_False);
424 for(i=1;i<NZer;i++) {
425 if(Abs(Zer(i+1)-Zer(i))<Eps) {
426 //-- est ce une racine double ou une erreur numerique ?
427 Standard_Real qw=Zer(i+1);
428 Standard_Real va=ko(4)+qw*(2.0*ko(3)+qw*(3.0*ko(2)+qw*(4.0*ko(1))));
429 //-- cout<<" Val Double ("<<qw<<")=("<<va<<")"<<endl;
435 // cout<<"Pb ds math_TrigonometricFunctionRoots CC="
436 // <<A<<" CS="<<B<<" C="<<C<<" S="<<D<<" Cte="<<E<<endl;
444 //-- Si il y a un coeff petit, on divise
458 // Verification des solutions par rapport aux bornes:
459 // ==================================================
460 Standard_Real SupmInfs100 = (SupBound-InfBound)*0.01;
462 for (i = 1; i <= NZer; i++) {
463 Teta = atan(Zer(i)); Teta+=Teta;
464 if (Zer(i) <= (-Eps)) {
465 Teta = Depi-Abs(Teta);
467 Teta += IntegerPart(Mod)*Depi;
468 if (Teta-MyBorneInf < 0) Teta += Depi;
470 X = Teta -MyBorneInf;
471 if ((X >= (-Epsilon(Delta))) && (X <= Delta+ Epsilon(Delta))) {
475 //OCC541(apo): Standard_Real TetaNewton=0;
476 Standard_Real TetaNewton = Teta;
477 MyTrigoFunction MyF(A, B, C, D, E);
478 math_NewtonFunctionRoot Resol(MyF, X, Tol1, Eps, Nit);
479 if (Resol.IsDone()) {
480 TetaNewton = Resol.Root();
482 //-- lbr le 7 mars 97 (newton converge tres tres loin de la solution initilale)
483 Standard_Real DeltaNewton = TetaNewton-Teta;
484 if((DeltaNewton > SupmInfs100) || (DeltaNewton < -SupmInfs100)) {
485 //-- cout<<"\n Newton X0="<<Teta<<" -> "<<TetaNewton<<endl;
491 Flag4 = Standard_False;
493 for(k = 1; k <= NbSol; k++) {
494 //On met les valeurs par ordre croissant:
496 for (l = k; l <= NbSol; l++) {
502 Flag4 = Standard_True;
512 // Cas particulier de PI:
514 Standard_Integer startIndex = NbSol + 1;
515 for( Standard_Integer solIt = startIndex; solIt <= 4; solIt++) {
516 Teta = M_PI + IntegerPart(Mod)*2.0*M_PI;;
517 X = Teta - MyBorneInf;
518 if ((X >= (-Epsilon(Delta))) && (X <= Delta + Epsilon(Delta))) {
519 if (Abs(A-C+E) <= Eps) {
520 Flag4 = Standard_False;
521 for (k = 1; k <= NbSol; k++) {
524 Flag4 = Standard_True;
527 if ((solIt == startIndex) && (Abs(Teta-Sol(k)) <= Eps)) {
537 for (k = j; k <= NbSol; k++) {
551 void math_TrigonometricFunctionRoots::Dump(Standard_OStream& o) const
553 o << " math_TrigonometricFunctionRoots: \n";
557 else if (InfiniteStatus) {
558 o << " There is an infinity of roots\n";
560 else if (!InfiniteStatus) {
561 o << " Number of solutions = " << NbSol <<"\n";
562 for (Standard_Integer i = 1; i <= NbSol; i++) {
563 o << " Value number " << i << "= " << Sol(i) << "\n";