src/Convert/Convert_PolynomialCosAndSin.cxx

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