0026937: Eliminate NO_CXX_EXCEPTION macro support
[occt.git] / src / Convert / Convert_PolynomialCosAndSin.cxx
1 // Created on: 1995-10-10
2 // Created by: Jacques GOUSSARD
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 under
9 // the terms of the GNU Lesser General Public License 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 #include <Convert_PolynomialCosAndSin.hxx>
18
19 #include <TColgp_Array1OfPnt2d.hxx>
20 #include <gp_Trsf2d.hxx>
21 #include <gp_Pnt2d.hxx>
22 #include <gp_Vec2d.hxx>
23 #include <gp_XY.hxx>
24
25 #include <gp.hxx>
26 #include <Precision.hxx>
27 #include <PLib.hxx>
28 #include <BSplCLib.hxx>
29
30 #include <Standard_ConstructionError.hxx>
31
32 static Standard_Real Locate(const Standard_Real Angfin,
33                             const TColgp_Array1OfPnt2d& TPoles,
34                             const Standard_Real Umin,
35                             const Standard_Real Umax)
36 {
37   Standard_Real umin = Umin;
38   Standard_Real umax = Umax;
39   Standard_Real Ptol = Precision::Angular();
40   Standard_Real Utol = Precision::PConfusion();
41   while (Abs(umax-umin)>= Utol) {
42     Standard_Real ptest = (umax+umin)/2.;
43     gp_Pnt2d valP;
44     BSplCLib::D0(ptest,TPoles,BSplCLib::NoWeights(),valP);
45     Standard_Real theta = ATan2(valP.Y(),valP.X());
46     if (theta < 0.) {
47       theta +=2.*M_PI;
48     }
49     if (Abs(theta - Angfin) < Ptol) {
50       return ptest;
51     }
52     if (theta < Angfin) {
53       umin = ptest;
54     }
55     else if (theta > Angfin) {
56       umax = ptest;
57     }
58   }
59   return (umin+umax)/2.;
60 }
61
62
63 void BuildPolynomialCosAndSin
64   (const Standard_Real UFirst,
65    const Standard_Real ULast,
66    const Standard_Integer num_poles,
67    Handle(TColStd_HArray1OfReal)& CosNumeratorPtr,
68    Handle(TColStd_HArray1OfReal)& SinNumeratorPtr,
69    Handle(TColStd_HArray1OfReal)& DenominatorPtr)
70 {
71
72   Standard_Real  Delta,
73   locUFirst,
74 //  locULast,
75 //  temp_value,
76   t_min,
77   t_max,
78   trim_min,
79   trim_max,
80   middle,
81   Angle,
82   PI2 = 2*M_PI ;
83   Standard_Integer ii, degree = num_poles -1 ;
84   locUFirst = UFirst ;
85
86   // Return UFirst in [-2PI; 2PI]
87   // to make rotations without risk
88   while (locUFirst > PI2) {
89     locUFirst -= PI2;
90   }
91   while (locUFirst < - PI2) {
92     locUFirst += PI2;
93   }
94
95 // Return to the arc [0, Delta]
96   Delta = ULast - UFirst;
97   middle =  0.5e0 * Delta ; 
98   
99   // coincide the required bisector of the angular sector with 
100   // axis -Ox definition of the circle in Bezier of degree 7 so that 
101   // parametre 1/2 of Bezier was exactly a point of the bissectrice 
102   // of the required angular sector.
103   //
104   Angle = middle - M_PI ;
105   //
106   // Circle of radius 1. See Euclid
107   //
108
109   TColgp_Array1OfPnt2d TPoles(1,8),
110   NewTPoles(1,8) ;
111   TPoles(1).SetCoord(1.,0.);
112   TPoles(2).SetCoord(1.,1.013854);
113   TPoles(3).SetCoord(-0.199043,1.871905);
114   TPoles(4).SetCoord(-1.937729,1.057323);
115   TPoles(5).SetCoord(-1.937729,-1.057323);
116   TPoles(6).SetCoord(-0.199043,-1.871905);
117   TPoles(7).SetCoord(1.,-1.013854);
118   TPoles(8).SetCoord(1.,0.);
119   gp_Trsf2d T;
120   T.SetRotation(gp::Origin2d(),Angle);
121   for (ii=1; ii<=num_poles; ii++) {
122     TPoles(ii).Transform(T);
123   }
124
125
126   t_min = 1.0e0 - (Delta * 1.3e0 / M_PI) ;
127   t_min *= 0.5e0 ;
128   t_min = Max(t_min,0.0e0) ;
129   t_max = 1.0e0 + (Delta * 1.3e0 / M_PI) ;
130   t_max *= 0.5e0 ;
131   t_max = Min(t_max,1.0e0) ;
132   trim_max = Locate(Delta,
133                     TPoles,
134                     t_min,
135                     t_max);
136   //
137   // as Bezier is symmetric correspondingly to the bissector 
138   // of the angular sector ...
139   
140   trim_min = 1.0e0 - trim_max ;
141   //
142   Standard_Real knot_array[2] ;
143   Standard_Integer  mults_array[2] ; 
144   knot_array[0] = 0.0e0 ;
145   knot_array[1] = 1.0e0 ;
146   mults_array[0] = degree + 1 ;
147   mults_array[1] = degree + 1 ;
148   
149   TColStd_Array1OfReal  the_knots(knot_array[0],1,2),
150   the_new_knots(knot_array[0],1,2);
151   TColStd_Array1OfInteger the_mults(mults_array[0],1,2),
152   the_new_mults(mults_array[0],1,2) ;
153   
154   BSplCLib::Trimming(degree,
155                      Standard_False,
156                      the_knots,
157                      the_mults,
158                      TPoles,
159                      BSplCLib::NoWeights(),
160                      trim_min,
161                      trim_max,
162                      the_new_knots,
163                      the_new_mults,
164                      NewTPoles,
165                      BSplCLib::NoWeights());
166
167   // readjustment is obviously redundant
168   Standard_Real SinD = Sin(Delta), CosD = Cos(Delta);
169   gp_Pnt2d Pdeb(1., 0.);
170   gp_Pnt2d Pfin(CosD, SinD);
171
172   Standard_Real dtg = NewTPoles(1).Distance(NewTPoles(2));
173   NewTPoles(1) = Pdeb;
174   gp_XY theXY(0.,dtg);
175   Pdeb.ChangeCoord() += theXY;
176   NewTPoles(2) = Pdeb;
177   
178   // readjustment to Euclid
179   dtg = NewTPoles(num_poles).Distance(NewTPoles(num_poles-1));
180   NewTPoles(num_poles) = Pfin;
181   theXY.SetCoord(dtg*SinD,-dtg*CosD);
182   Pfin.ChangeCoord() += theXY;
183   NewTPoles(num_poles-1) = Pfin;
184
185   // Rotation to return to the arc [LocUFirst, LocUFirst+Delta]
186   T.SetRotation(gp::Origin2d(), locUFirst);
187   for (ii=1; ii<=num_poles; ii++) {
188     NewTPoles(ii).Transform(T);
189   }  
190   
191   for (ii=1; ii<=num_poles; ii++) {
192     CosNumeratorPtr->SetValue(ii,NewTPoles(ii).X());
193     SinNumeratorPtr->SetValue(ii,NewTPoles(ii).Y());
194     DenominatorPtr->SetValue(ii,1.);
195   }
196 }
197
198 /*
199 void BuildHermitePolynomialCosAndSin
200   (const Standard_Real UFirst,
201    const Standard_Real ULast,
202    const Standard_Integer num_poles,
203    Handle(TColStd_HArray1OfReal)& CosNumeratorPtr,
204    Handle(TColStd_HArray1OfReal)& SinNumeratorPtr,
205    Handle(TColStd_HArray1OfReal)& DenominatorPtr)
206 {
207
208   if (num_poles%2 != 0) {
209     throw Standard_ConstructionError();
210   }
211   Standard_Integer ii;
212   Standard_Integer ordre_deriv = num_poles/2;
213   Standard_Real ang = ULast - UFirst;
214   Standard_Real Cd = Cos(UFirst);
215   Standard_Real Sd = Sin(UFirst);
216   Standard_Real Cf = Cos(ULast);
217   Standard_Real Sf = Sin(ULast);
218   
219   Standard_Integer Degree = num_poles-1;
220   TColStd_Array1OfReal FlatKnots(1,2*num_poles);
221   TColStd_Array1OfReal Parameters(1,num_poles);
222   TColStd_Array1OfInteger ContactOrderArray(1,num_poles);
223   TColgp_Array1OfPnt2d Poles(1,num_poles);
224   TColgp_Array1OfPnt2d TPoles(1,num_poles);
225  
226
227   for (ii=1; ii<=num_poles; ii++) {
228     FlatKnots(ii) = 0.;
229     FlatKnots(ii+num_poles) = 1.;
230   }
231
232   Standard_Real coef = 1.;
233   Standard_Real xd,yd,xf,yf; 
234
235   for (ii=1; ii<=ordre_deriv; ii++) {
236     Parameters(ii) = 0.;
237     Parameters(ii+ordre_deriv) = 1.;
238
239     ContactOrderArray(ii) = ContactOrderArray(num_poles-ii+1) = ii-1;
240
241     switch ((ii-1)%4) {
242     case 0:
243       {
244         xd = Cd*coef;
245         yd = Sd*coef;
246         xf = Cf*coef;
247         yf = Sf*coef;
248       }
249       break;
250     case 1:
251       {
252         xd = -Sd*coef;
253         yd =  Cd*coef;
254         xf = -Sf*coef;
255         yf =  Cf*coef;
256       }
257       break;
258     case 2:
259       {
260         xd = -Cd*coef;
261         yd = -Sd*coef;
262         xf = -Cf*coef;
263         yf = -Sf*coef;
264       }
265       break;
266     case 3:
267       {
268         xd =  Sd*coef;
269         yd = -Cd*coef;
270         xf =  Sf*coef;
271         yf = -Cf*coef;
272       }
273       break;
274     }
275
276     Poles(ii).SetX(xd);
277     Poles(ii).SetY(yd);
278     Poles(num_poles-ii+1).SetX(xf);
279     Poles(num_poles-ii+1).SetY(yf);
280
281     coef *= ang;
282   }
283
284   Standard_Integer InversionPb;
285   BSplCLib::Interpolate(Degree,FlatKnots,Parameters,
286                         ContactOrderArray,Poles,InversionPb);
287
288   if (InversionPb !=0) {
289     throw Standard_ConstructionError();
290   }
291   for (ii=1; ii<=num_poles; ii++) {
292     CosNumeratorPtr->SetValue(ii,Poles(ii).X());
293     SinNumeratorPtr->SetValue(ii,Poles(ii).Y());
294     DenominatorPtr->SetValue(ii,1.);
295   }
296
297 }
298 */