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_TrigonometricEquationFunction.hxx>
29 #include <math_DirectPolynomialRoots.hxx>
30 #include <Standard_OutOfRange.hxx>
31 #include <math_FunctionWithDerivative.hxx>
32 #include <math_NewtonFunctionRoot.hxx>
33 #include <Precision.hxx>
35 math_TrigonometricFunctionRoots::math_TrigonometricFunctionRoots
36 (const Standard_Real theD,
37 const Standard_Real theE,
38 const Standard_Real theInfBound,
39 const Standard_Real theSupBound)
42 InfiniteStatus(Standard_False),
45 const Standard_Real A(0.0), B(0.0), C(0.0);
46 Perform(A, B, C, theD, theE, theInfBound, theSupBound);
50 math_TrigonometricFunctionRoots::math_TrigonometricFunctionRoots
51 (const Standard_Real theC,
52 const Standard_Real theD,
53 const Standard_Real theE,
54 const Standard_Real theInfBound,
55 const Standard_Real theSupBound)
58 InfiniteStatus(Standard_False),
61 const Standard_Real A(0.0), B(0.0);
62 Perform(A, B, theC, theD, theE, theInfBound, theSupBound);
67 math_TrigonometricFunctionRoots::math_TrigonometricFunctionRoots
68 (const Standard_Real theA,
69 const Standard_Real theB,
70 const Standard_Real theC,
71 const Standard_Real theD,
72 const Standard_Real theE,
73 const Standard_Real theInfBound,
74 const Standard_Real theSupBound)
77 InfiniteStatus(Standard_False),
80 Perform(theA, theB, theC, theD, theE, theInfBound, theSupBound);
83 void math_TrigonometricFunctionRoots::Perform(const Standard_Real A,
84 const Standard_Real B,
85 const Standard_Real C,
86 const Standard_Real D,
87 const Standard_Real E,
88 const Standard_Real InfBound,
89 const Standard_Real SupBound) {
91 Standard_Integer i, j=0, k, l, NZer=0, Nit = 10;
92 Standard_Real Depi, Delta, Mod, AA, BB, CC, MyBorneInf;
93 Standard_Real Teta, X;
94 Standard_Real Eps, Tol1 = 1.e-15;
95 TColStd_Array1OfReal ko(1,5), Zer(1,4);
96 Standard_Boolean Flag4;
97 InfiniteStatus = Standard_False;
103 if (InfBound <= RealFirst() && SupBound >= RealLast()) {
108 else if (SupBound >= RealLast()) {
109 MyBorneInf = InfBound;
111 Mod = MyBorneInf/Depi;
113 else if (InfBound <= RealFirst()) {
114 MyBorneInf = SupBound - Depi;
116 Mod = MyBorneInf/Depi;
119 MyBorneInf = InfBound;
120 Delta = SupBound-InfBound;
122 if ((SupBound-InfBound) > Depi) { Delta = Depi;}
125 if ((Abs(A) <= Eps) && (Abs(B) <= Eps)) {
129 InfiniteStatus = Standard_True; // infinite de solutions.
138 // Equation du type d*sin(x) + e = 0
139 // =================================
147 Zer(2) = M_PI - Zer(1);
149 for (i = 1; i <= NZer; i++) {
150 if (Zer(i) <= -Eps) {
151 Zer(i) = Depi - Abs(Zer(i));
153 // On rend les solutions entre InfBound et SupBound:
154 // =================================================
155 Zer(i) += IntegerPart(Mod)*Depi;
156 X = Zer(i)-MyBorneInf;
157 if ((X > (-Epsilon(Delta))) && (X < Delta+ Epsilon(Delta))) {
165 else if (Abs(D) <= Eps) {
167 // Equation du premier degre de la forme c*cos(x) + e = 0
168 // ======================================================
178 for (i = 1; i <= NZer; i++) {
179 if (Zer(i) <= -Eps) {
180 Zer(i) = Depi - Abs(Zer(i));
182 // On rend les solutions entre InfBound et SupBound:
183 // =================================================
184 Zer(i) += IntegerPart(Mod)*2.*M_PI;
185 X = Zer(i)-MyBorneInf;
186 if ((X >= (-Epsilon(Delta))) && (X <= Delta+ Epsilon(Delta))) {
195 // Equation du second degre:
196 // =========================
201 math_DirectPolynomialRoots Resol(AA, BB, CC);
202 if (!Resol.IsDone()) {
203 Done = Standard_False;
206 else if(!Resol.InfiniteRoots()) {
207 NZer = Resol.NbSolutions();
208 for (i = 1; i <= NZer; i++) {
209 Zer(i) = Resol.Value(i);
212 else if (Resol.InfiniteRoots()) {
213 InfiniteStatus = Standard_True;
219 // Two additional analytical cases.
220 if ((Abs(A) <= Eps) &&
225 // 2 * B * sin * cos + D * sin = 0
231 if (Abs(AA) <= 1.0 + Precision::PConfusion())
247 Zer(4) = Depi - Zer(3);
252 for (i = 1; i <= NZer; i++)
254 if (Zer(i) <= MyBorneInf - Eps)
258 // On rend les solutions entre InfBound et SupBound:
259 // =================================================
260 Zer(i) += IntegerPart(Mod)*2.*M_PI;
261 X = Zer(i)-MyBorneInf;
262 if ((X >= (-Precision::PConfusion())) &&
263 (X <= Delta + Precision::PConfusion()))
265 if (Zer(i) < InfBound)
267 if (Zer(i) > SupBound)
277 // 2 * B * sin * cos + C * cos = 0
280 Zer(2) = M_PI * 3.0 / 2.0;
283 if (Abs(AA) <= 1.0 + Precision::PConfusion())
294 Zer(3) = M_PI * 3.0 / 2.0;
295 Zer(4) = M_PI * 3.0 / 2.0;
300 Zer(4) = M_PI - Zer(3);
305 for (i = 1; i <= NZer; i++)
307 if (Zer(i) <= MyBorneInf - Eps)
311 // On rend les solutions entre InfBound et SupBound:
312 // =================================================
313 Zer(i) += IntegerPart(Mod)*2.*M_PI;
314 X = Zer(i)-MyBorneInf;
315 if ((X >= (-Precision::PConfusion())) &&
316 (X <= Delta + Precision::PConfusion()))
318 if (Zer(i) < InfBound)
320 if (Zer(i) > SupBound)
330 // Equation du 4 ieme degre
331 // ========================
337 Standard_Boolean bko;
338 Standard_Integer nbko=0;
341 math_DirectPolynomialRoots Resol4(ko(1), ko(2), ko(3), ko(4), ko(5));
342 if (!Resol4.IsDone()) {
343 Done = Standard_False;
346 else if (!Resol4.InfiniteRoots()) {
347 NZer = Resol4.NbSolutions();
348 for (i = 1; i <= NZer; i++) {
349 Zer(i) = Resol4.Value(i);
352 else if (Resol4.InfiniteRoots()) {
353 InfiniteStatus = Standard_True;
356 Standard_Boolean triok;
359 for(i=1;i<NZer;i++) {
360 if(Zer(i)>Zer(i+1)) {
361 Standard_Real t=Zer(i);
364 triok=Standard_False;
368 while(triok==Standard_False);
370 for(i=1;i<NZer;i++) {
371 if(Abs(Zer(i+1)-Zer(i))<Eps) {
372 //-- est ce une racine double ou une erreur numerique ?
373 Standard_Real qw=Zer(i+1);
374 Standard_Real va=ko(4)+qw*(2.0*ko(3)+qw*(3.0*ko(2)+qw*(4.0*ko(1))));
375 //-- std::cout<<" Val Double ("<<qw<<")=("<<va<<")"<<std::endl;
381 // std::cout<<"Pb ds math_TrigonometricFunctionRoots CC="
382 // <<A<<" CS="<<B<<" C="<<C<<" S="<<D<<" Cte="<<E<<std::endl;
390 //-- Si il y a un coeff petit, on divise
404 // Verification des solutions par rapport aux bornes:
405 // ==================================================
406 Standard_Real SupmInfs100 = (SupBound-InfBound)*0.01;
408 for (i = 1; i <= NZer; i++) {
409 Teta = atan(Zer(i)); Teta+=Teta;
410 if (Zer(i) <= (-Eps)) {
411 Teta = Depi-Abs(Teta);
413 Teta += IntegerPart(Mod)*Depi;
414 if (Teta-MyBorneInf < 0) Teta += Depi;
416 X = Teta -MyBorneInf;
417 if ((X >= (-Epsilon(Delta))) && (X <= Delta+ Epsilon(Delta))) {
421 //OCC541(apo): Standard_Real TetaNewton=0;
422 Standard_Real TetaNewton = Teta;
423 math_TrigonometricEquationFunction MyF(A, B, C, D, E);
424 math_NewtonFunctionRoot Resol(MyF, X, Tol1, Eps, Nit);
425 if (Resol.IsDone()) {
426 TetaNewton = Resol.Root();
428 //-- lbr le 7 mars 97 (newton converge tres tres loin de la solution initilale)
429 Standard_Real DeltaNewton = TetaNewton-Teta;
430 if((DeltaNewton > SupmInfs100) || (DeltaNewton < -SupmInfs100)) {
431 //-- std::cout<<"\n Newton X0="<<Teta<<" -> "<<TetaNewton<<std::endl;
437 Flag4 = Standard_False;
439 for(k = 1; k <= NbSol; k++) {
440 //On met les valeurs par ordre croissant:
442 for (l = k; l <= NbSol; l++) {
448 Flag4 = Standard_True;
458 // Cas particulier de PI:
460 Standard_Integer startIndex = NbSol + 1;
461 for( Standard_Integer solIt = startIndex; solIt <= 4; solIt++) {
462 Teta = M_PI + IntegerPart(Mod)*2.0*M_PI;;
463 X = Teta - MyBorneInf;
464 if ((X >= (-Epsilon(Delta))) && (X <= Delta + Epsilon(Delta))) {
465 if (Abs(A-C+E) <= Eps) {
466 Flag4 = Standard_False;
467 for (k = 1; k <= NbSol; k++) {
470 Flag4 = Standard_True;
473 if ((solIt == startIndex) && (Abs(Teta-Sol(k)) <= Eps)) {
483 for (k = j; k <= NbSol; k++) {
497 void math_TrigonometricFunctionRoots::Dump(Standard_OStream& o) const
499 o << " math_TrigonometricFunctionRoots: \n";
503 else if (InfiniteStatus) {
504 o << " There is an infinity of roots\n";
506 else if (!InfiniteStatus) {
507 o << " Number of solutions = " << NbSol <<"\n";
508 for (Standard_Integer i = 1; i <= NbSol; i++) {
509 o << " Value number " << i << "= " << Sol(i) << "\n";