[occt.git] / src / GeomFill / GeomFill_QuasiAngularConvertor.cxx

// Created on: 1997-08-06
// Created by: Philippe MANGIN
// Copyright (c) 1997-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.


#include <Convert_CompPolynomialToPoles.hxx>
#include <GeomFill_QuasiAngularConvertor.hxx>
#include <gp_Ax1.hxx>
#include <gp_Mat.hxx>
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>
#include <PLib.hxx>
#include <StdFail_NotDone.hxx>
#include <TColStd_Array1OfReal.hxx>
#include <TColStd_HArray1OfInteger.hxx>
#include <TColStd_HArray1OfReal.hxx>
#include <TColStd_HArray2OfReal.hxx>

#define NullAngle 1.e-6

// QuasiAngular is rational definition of Cos(theta(t) and sin(theta)
// on [-alpha, +alpha] with
//                     2      2
//                    U   -  V 
//  cos (theta(t)) = ----------
//                     2      2
//                    U   +  V 

//                      2 * U*V
//   sin (theta(t)) = ----------
//                      2      2
//                     U   +  V 
//                     3
//       V(t) = t + c t
//                     2
//       U(t) = 1 + b t
//            1
//       c = ---  + b   
//            3 
//            -1                     gamma
//       b =---------  +   -----------------------  
//                 2 
//            gamma         3*(tang gamma - gamma) 
//     with gamma = alpha / 2


GeomFill_QuasiAngularConvertor::GeomFill_QuasiAngularConvertor(): 
                    myinit(Standard_False),
		    B(1, 7, 1, 7), 
		    Px(1, 7), Py(1, 7), 
		    W(1,7), Vx(1, 7), Vy(1, 7), Vw(1,7)
{
}

Standard_Boolean GeomFill_QuasiAngularConvertor::Initialized() const
{
 return myinit;
}

void GeomFill_QuasiAngularConvertor::Init() 
{
  if (myinit) return; //On n'initialise qu'une fois
  Standard_Integer ii, jj, Ordre=7;
  Standard_Real terme;
  TColStd_Array1OfReal Coeffs(1, Ordre*Ordre), TrueInter(1,2), Inter(1,2);
  Handle(TColStd_HArray2OfReal) 
    Poles1d = new (TColStd_HArray2OfReal) (1, Ordre, 1, Ordre);

  //Calcul de B
  Inter.SetValue(1, -1);
  Inter.SetValue(2,  1);
  TrueInter.SetValue(1, -1);
  TrueInter.SetValue(2,  1);

  Coeffs.Init(0);
  for (ii=1; ii<=Ordre; ii++) {  Coeffs.SetValue(ii+(ii-1)*Ordre, 1); }

   //Convertion
  Convert_CompPolynomialToPoles 
	 AConverter(Ordre, Ordre-1,  Ordre-1,
		    Coeffs, 
		    Inter,
		    TrueInter);
  AConverter.Poles(Poles1d);

  for (jj=1; jj<=Ordre; jj++) {
      for (ii=1; ii<=Ordre; ii++) {
        terme = Poles1d->Value(ii,jj);
        if (Abs(terme-1) < 1.e-9) terme = 1 ; //petite retouche
        if (Abs(terme+1) < 1.e-9) terme = -1; 
	B(ii, jj) =  terme;
      }
  }
 
  // Init des polynomes
  Vx.Init(0);
  Vx(1) = 1;
  Vy.Init(0);
  Vy(2) = 2;
  Vw.Init(0);
  Vw(1) = 1;
  myinit = Standard_True;
}

void GeomFill_QuasiAngularConvertor::Section(const gp_Pnt& FirstPnt,
				       const gp_Pnt& Center,
				       const gp_Vec& Dir,
				       const Standard_Real Angle,
				       TColgp_Array1OfPnt& Poles,
				       TColStd_Array1OfReal& Weights) 
{ 
  Standard_Real b, b2, c, c2,tan_b;
  Standard_Integer ii;
  Standard_Real beta, beta2, beta3, beta4, beta5, beta6, wi;
  gp_XYZ aux;
  gp_Mat Rot;
  // Calcul de la transformation
  gp_Vec V1(Center, FirstPnt), V2;
  Rot.SetRotation(Dir.XYZ(), Angle/2);
  aux = V1.XYZ();
  aux *= Rot;
  V1.SetXYZ(aux);
  V2 =  Dir^V1;
 
  gp_Mat M(V1.X(), V2.X(), 0, 
	   V1.Y(), V2.Y(), 0,
	   V1.Z(), V2.Z(), 0);

  // Calcul des coeffs  -----------
  beta = Angle/4;
  beta2 = beta * beta;
  beta3 = beta * beta2;
  beta4 = beta2*beta2;
  beta5 = beta3*beta2;
  beta6 = beta3*beta3; 
 
  if  ((M_PI/2 - beta)> NullAngle) {
    if (Abs(beta) < NullAngle) {
     Standard_Real cf = 2.0/(3*5*7);
     b = - (0.2+cf*beta2) / (1+ 0.2*beta2);
//     b = beta5 / cf; 
    }
    else {
      tan_b = Tan(beta);
      b = - 1.0e0 / beta2;
      b += beta / (3*(tan_b - beta));
    }
  }
  else b = ((Standard_Real) -1)/beta2;
  c =  ((Standard_Real) 1)/ 3 + b;
  b2 = b*b;
  c2 = c*c;

  // X = U*U - V*V
  Vx(3) = beta2*(2*b - 1);
  Vx(5) = beta4*(b2 - 2*c); 
  Vx(7) = -beta6*c2;

  //Y = 2*U*V
  Vy(2) = 2*beta;
  Vy(4) = beta3*2*(c+b);
  Vy(6) = 2*beta5*b*c;

 // W = U*U + V*V
  Vw(3) = beta2*(1 + 2*b);
  Vw(5) = beta4*(2*c + b2); 
  Vw(7) = beta6*c2;   
    
// Calculs des poles
  Px.Multiply(B, Vx);
  Py.Multiply(B, Vy);
  W.Multiply(B, Vw);

// Transfo
  gp_XYZ pnt;
  for (ii=1; ii<=7; ii++) {
     wi = W(ii);
     pnt.SetCoord(Px(ii)/wi, Py(ii)/wi, 0);
     pnt *= M;
     pnt += Center.XYZ();
     Poles(ii).ChangeCoord() = pnt;
     Weights(ii) = wi;
  }
}

void GeomFill_QuasiAngularConvertor::Section(const gp_Pnt& FirstPnt,
				       const gp_Vec& DFirstPnt,
				       const gp_Pnt& Center,
				       const gp_Vec& DCenter,
				       const gp_Vec& Dir,
				       const gp_Vec& DDir,
				       const Standard_Real Angle,
				       const Standard_Real DAngle,
				       TColgp_Array1OfPnt& Poles,
				       TColgp_Array1OfVec& DPoles,
				       TColStd_Array1OfReal& Weights,
				       TColStd_Array1OfReal& DWeights) 
{
  Standard_Integer Ordre = 7;
  math_Vector DVx(1, Ordre), DVy(1, Ordre), DVw(1, Ordre),
              DPx(1, Ordre), DPy(1, Ordre), DW(1, Ordre);
  Standard_Real b, b2, c, c2,tan_b;
  Standard_Real bpr, dtan_b;
  Standard_Integer ii;
  Standard_Real beta, beta2, beta3, beta4, beta5, beta6, betaprim;
  gp_Vec V1(Center, FirstPnt), V1Prim, V2;

  // Calcul des  transformations
  gp_XYZ aux;
  Standard_Real Sina, Cosa;
  gp_Mat Rot, RotPrim, D, DPrim;
  // La rotation s'ecrit I +  sin(Ang) * D + (1. - cos(Ang)) * D*D
  // ou D est l'application x -> Dir ^ x
  Rot.SetRotation(Dir.XYZ(), Angle/2);
  // La derive s'ecrit donc :
  // AngPrim * (sin(Ang)*D*D + cos(Ang)*D) 
  // + sin(Ang)*DPrim  + (1. - cos(Ang)) *(DPrim*D + D*DPrim)
  Sina = Sin(Angle/2); 
  Cosa = Cos(Angle/2);
  D.SetCross(Dir.XYZ());
  DPrim.SetCross(DDir.XYZ());

  RotPrim =  (D.Powered(2)).Multiplied(Sina);
  RotPrim += D.Multiplied(Cosa);
  RotPrim *= DAngle/2;
  RotPrim += DPrim.Multiplied(Sina);
  RotPrim += ((DPrim.Multiplied(D)).Added(D.Multiplied(DPrim))).Multiplied(1-Cosa);
  
  aux = (DFirstPnt - DCenter).XYZ().Multiplied(Rot);
  aux += V1.XYZ().Multiplied(RotPrim);
  V1Prim.SetXYZ(aux);
  aux = V1.XYZ();
  aux.Multiply(Rot);
  V1.SetXYZ(aux);
  V2 =  Dir^V1;
  gp_Mat M (V1.X(), V2.X(), 0,
	    V1.Y(), V2.Y(), 0,
	    V1.Z(), V2.Z(), 0);
  V2 = (DDir.Crossed(V1)).Added(Dir.Crossed(V1Prim));
  gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
		V1Prim.Y(), V2.Y(), 0,
		V1Prim.Z(), V2.Z(), 0);

  // Calcul des constante  -----------
  beta = Angle/4;
  betaprim = DAngle/4;
  beta2 = beta * beta;
  beta3 = beta * beta2;
  beta4 = beta2*beta2;
  beta5 = beta3*beta2;
  beta6 = beta3*beta3;

  if (Abs(beta) < NullAngle) {
    // On calcul b par D.L
     Standard_Real cf = 2.0/(3*5*7);
     Standard_Real Num, Denom;
     Num = 0.2 + cf*beta2;
     Denom = 1+0.2*beta2;
     b = - Num/Denom;
     bpr = -2*beta*betaprim*(cf*Denom - 0.2*Num)/(Denom*Denom);
  }
  else {
    b = ((Standard_Real) -1)/beta2;
    bpr = (2*betaprim) / beta3;  
    if  ((M_PI/2 - beta)> NullAngle) {
      tan_b = Tan(beta);
      dtan_b = betaprim * (1 + tan_b*tan_b);
      b2 = tan_b - beta;
      b += beta / (3*b2);
      bpr += (betaprim*tan_b - beta*dtan_b) / (3*b2*b2);
    }
  }

  c =  ((Standard_Real) 1)/ 3 + b;
  b2 = b*b;
  c2 = c*c;


  // X = U*U - V*V
  Vx(3) = beta2*(2*b - 1);
  Vx(5) = beta4*(b2 - 2*c); 
  Vx(7) = -beta6*c2;
  DVx.Init(0);
  DVx(3) = 2*(beta*betaprim*(2*b - 1) + bpr*beta2);
  DVx(5) = 4*beta3*betaprim*(b2 - 2*c) + 2*beta4*bpr*(b-1); 
  DVx(7) = - 6*beta5*betaprim*c2 - 2*beta6*bpr*c;  

  //Y = 2*U*V
  Vy(2) = 2*beta;
  Vy(4) = beta3*2*(c+b);
  Vy(6) = 2*beta5*b*c;
  DVy.Init(0);
  DVy(2) = 2*betaprim;
  DVy(4) = 6*beta2*betaprim*(b+c) + 4*beta3*bpr;
  DVy(6) = 10*beta4*betaprim*b*c + 2*beta5*bpr*(b+c);

 // W = U*U + V*V
  Vw(3) = beta2*(1 + 2*b);
  Vw(5) = beta4*(2*c + b2); 
  Vw(7) = beta6*c2;   
  DVw.Init(0);
//  DVw(3) = 2*(beta*betaprim*(1 + 2*b) + beta2*bpr);
  DVw(3) = 2*beta*(betaprim*(1 + 2*b) + beta*bpr);
//  DVw(5) = 4*beta3*betaprim*(2*c + b2) + 2*beta4*bpr*(b+1);
  DVw(5) = 2*beta3*(2*betaprim*(2*c + b2) + beta*bpr*(b+1)); 
//  DVw(7) = 6*beta5*betaprim*c2 + 2*beta6*bpr*c;
  DVw(7) = 2*beta5*c*(3*betaprim*c + beta*bpr);

 
  // Calcul des poles
  Px.Multiply(B, Vx);
  Py.Multiply(B, Vy);
  W.Multiply(B, Vw);
  DPx.Multiply(B, DVx);
  DPy.Multiply(B, DVy);
  DW.Multiply(B, DVw);

  gp_XYZ P, DP;
  Standard_Real wi;

  for (ii=1; ii<=Ordre; ii++) {
    wi = W(ii);
    P.SetCoord(Px(ii)/wi, Py(ii)/wi, 0);
    DP.SetCoord(DPx(ii)/wi, DPy(ii)/wi, 0);
    DP -= (DW(ii)/wi)*P;

    Poles(ii).ChangeCoord() = M*P + Center.XYZ();
    P *= MPrim;
    DP *= M;
    aux.SetLinearForm(1, P, 1, DP, DCenter.XYZ());
    DPoles(ii).SetXYZ(aux);
    Weights(ii) = wi;
    DWeights(ii) = DW(ii);
  }
}

void GeomFill_QuasiAngularConvertor::Section(const gp_Pnt& FirstPnt,
				       const gp_Vec& DFirstPnt,
				       const gp_Vec& D2FirstPnt,
				       const gp_Pnt& Center,
				       const gp_Vec& DCenter,
				       const gp_Vec& D2Center,
				       const gp_Vec& Dir,
				       const gp_Vec& DDir,
				       const gp_Vec& D2Dir,
				       const Standard_Real Angle,
				       const Standard_Real DAngle,
				       const Standard_Real D2Angle,
				       TColgp_Array1OfPnt& Poles,
				       TColgp_Array1OfVec& DPoles,
				       TColgp_Array1OfVec& D2Poles,
				       TColStd_Array1OfReal& Weights,
				       TColStd_Array1OfReal& DWeights,
				       TColStd_Array1OfReal& D2Weights) 
{
  Standard_Integer Ordre = 7;
  math_Vector DVx(1, Ordre), DVy(1, Ordre),  DVw(1, Ordre),
              D2Vx(1, Ordre), D2Vy(1, Ordre),  D2Vw(1, Ordre);
  math_Vector DPx(1, Ordre), DPy(1, Ordre), DW(1, Ordre),
              D2Px(1, Ordre), D2Py(1, Ordre), D2W(1, Ordre);
 
  Standard_Integer ii;
  Standard_Real daux, b, b2, c, c2, bpr, bsc;
  gp_Vec V1(Center, FirstPnt), V1Prim, V1Secn, V2;

  // Calcul des  transformations
  gp_XYZ auxyz;
  Standard_Real Sina, Cosa;
  gp_Mat Rot, RotPrim, RotSecn, D, DPrim, DSecn, DDP, Maux;
  // La rotation s'ecrit I +  sin(Ang) * D + (1. - cos(Ang)) * D*D
  // ou D est l'application x -> Dir ^ x
  Rot.SetRotation(Dir.XYZ(), Angle/2);
  // La derive s'ecrit donc :
  // AngPrim * (sin(Ang)*D*D + cos(Ang)*D) 
  // + sin(Ang)*DPrim  + (1. - cos(Ang)) *(DPrim*D + D*DPrim)
  Sina = Sin(Angle/2); 
  Cosa = Cos(Angle/2);
  D.SetCross(Dir.XYZ());
  DPrim.SetCross(DDir.XYZ());
  DSecn.SetCross(D2Dir.XYZ());
  
  DDP = (DPrim.Multiplied(D)).Added(D.Multiplied(DPrim));
  RotPrim =  (D.Powered(2)).Multiplied(Sina);
  RotPrim += D.Multiplied(Cosa);
  RotPrim *= DAngle/2;
  RotPrim += DPrim.Multiplied(Sina);
  RotPrim += DDP.Multiplied(1-Cosa);

  RotSecn =  (D.Powered(2)).Multiplied(Sina);
  RotSecn += D.Multiplied(Cosa);
  RotSecn *= D2Angle/2;
  Maux = (D.Powered(2)).Multiplied(Cosa);
  Maux -= D.Multiplied(Sina);
  Maux *= DAngle/2;
  Maux += DDP.Multiplied(2*Sina);
  Maux += DPrim.Multiplied(2*Cosa);
  Maux *= DAngle/2;
  RotSecn += Maux;
  Maux = (DSecn.Multiplied(D)).Added(D.Multiplied(DSecn));
  Maux +=  (DPrim.Powered(2)).Multiplied(2);
  Maux *= 1 - Cosa;
  Maux += DSecn.Multiplied(Sina);
  RotSecn += Maux;
  
  V1Prim = DFirstPnt - DCenter;
  auxyz = (D2FirstPnt - D2Center).XYZ().Multiplied(Rot);
  auxyz += 2*(V1Prim.XYZ().Multiplied(RotPrim));
  auxyz += V1.XYZ().Multiplied(RotSecn);
  V1Secn.SetXYZ(auxyz);
  auxyz = V1Prim.XYZ().Multiplied(Rot);
  auxyz += V1.XYZ().Multiplied(RotPrim);
  V1Prim.SetXYZ(auxyz);
  auxyz = V1.XYZ();
  auxyz.Multiply(Rot);
  V1.SetXYZ(auxyz);
  V2 =  Dir^V1;

  gp_Mat M (V1.X(), V2.X(), 0,
	    V1.Y(), V2.Y(), 0,
	    V1.Z(), V2.Z(), 0);
  V2 = (DDir.Crossed(V1)).Added(Dir.Crossed(V1Prim));
  gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
		V1Prim.Y(), V2.Y(), 0,
		V1Prim.Z(), V2.Z(), 0);
 
  V2 = DDir^V1Prim;
  V2 *= 2;
  V2 += (D2Dir.Crossed(V1)).Added(Dir.Crossed(V1Secn));
  gp_Mat MSecn (V1Secn.X(), V2.X(), 0,
		V1Secn.Y(), V2.Y(), 0,
		V1Secn.Z(), V2.Z(), 0);

     
  // Calcul des coeff  -----------
  Standard_Real tan_b, dtan_b, d2tan_b;
  Standard_Real beta, beta2, beta3, beta4, beta5, beta6, betaprim, betasecn;
  Standard_Real betaprim2, bpr2;
  beta = Angle/4;
  betaprim = DAngle/4;
  betasecn = D2Angle/4;
  beta2 = beta * beta;
  beta3 = beta * beta2;
  beta4 = beta2*beta2;
  beta5 = beta3*beta2;
  beta6 = beta3*beta3;
  betaprim2 =  betaprim * betaprim;

 if (Abs(beta) < NullAngle) {
    // On calcul b par D.L
     Standard_Real cf =-2.0/21;
     Standard_Real Num, Denom, aux;
     Num = 0.2 + cf*beta2;
     Denom = 1+0.2*beta2;
     aux = (cf*Denom - 0.2*Num)/(Denom*Denom);
     b = - Num/Denom;
     bpr = -2*beta*betaprim*aux;
     bsc = 2*aux*(betaprim2 + beta*betasecn - 2*beta*betaprim2);
    
  }
  else { 
    b = ((Standard_Real) -1)/beta2;
    bpr = (2*betaprim) / beta3;
    bsc = (2*betasecn - 6*betaprim*(betaprim/beta)) / beta3; 
    if  ((M_PI/2 - beta)> NullAngle) {
      tan_b = Tan(beta);
      dtan_b = betaprim * (1 + tan_b*tan_b);
      d2tan_b = betasecn * (1 + tan_b*tan_b)
	+ 2*betaprim * tan_b * dtan_b;
      b2 = tan_b - beta;
      b += beta / (3*b2);
      Standard_Real aux =  betaprim*tan_b - beta*dtan_b;
      bpr += aux / (3*b2*b2);
      daux =  betasecn*tan_b - beta*d2tan_b;
      bsc += (daux - 2*aux*betaprim*tan_b*tan_b/b2)/(3*b2*b2);
    }
  }

  c =  ((Standard_Real) 1)/ 3 + b;
  b2 = b*b;
  c2 = c*c;
  bpr2 = bpr * bpr;


  // X = U*U - V*V
  Vx(3) = beta2*(2*b - 1);
  Vx(5) = beta4*(b2 - 2*c); 
  Vx(7) = -beta6*c2;
  DVx.Init(0);
  DVx(3) = 2*(beta*betaprim*(2*b - 1) + bpr*beta2);
  DVx(5) = 4*beta3*betaprim*(b2 - 2*c) + 2*beta4*bpr*(b-1); 
  DVx(7) = - 6*beta5*betaprim*c2 - 2*beta6*bpr*c;
  D2Vx.Init(0);  
  D2Vx(3) = 2*((betaprim2+beta*betasecn)*(2*b - 1) 
	       + 8*beta*betaprim*bpr 
	       + bsc*beta2);
  D2Vx(5) =  4*(b2 - 2*c)*(3*beta2*betaprim2 + beta3*betasecn)
             + 16*beta3*betaprim*bpr*(b-1)
             + 2*beta4*(bsc*(b-1)+bpr2);
  D2Vx(7) = - 6 * c2 * (5*beta4*betaprim2+beta5*betasecn)
            - 24*beta5*betaprim*bpr*c
	    - 2*beta6*(bsc*c + bpr2);

  //Y = 2*U*V
  Vy(2) = 2*beta;
  Vy(4) = beta3*2*(c+b);
  Vy(6) = 2*beta5*b*c;
  DVy.Init(0);
  DVy(2) = 2*betaprim;
  DVy(4) = 6*beta2*betaprim*(b+c) + 4*beta3*bpr;
  DVy(6) = 10*beta4*betaprim*b*c + 2*beta5*bpr*(b+c);
  D2Vy.Init(0);
  D2Vy(2) = 2*betasecn;
  D2Vy(4) = 6*(b+c)*(2*beta*betaprim2 + beta2*betasecn) 
          + 24*beta2*betaprim*bpr*(b+c)
	  + 4*beta3*bsc;
  D2Vy(6) = 10*b*c*(4*beta3*betaprim2 + beta4*betasecn)
          + 40 * beta4*betaprim*bpr*(b+c)
	  + 2*beta5*(bsc*(b+c)+ 2*bpr2);

 // W = U*U + V*V
  Vw(3) = beta2*(1 + 2*b);
  Vw(5) = beta4*(2*c + b2); 
  Vw(7) = beta6*c2;   
  DVw.Init(0);
  DVw(3) = 2*(beta*betaprim*(1 + 2*b) + beta2*bpr);
  DVw(5) = 4*beta3*betaprim*(2*c + b2) + 2*beta4*bpr*(b+1); 
  DVw(7) = 6*beta5*betaprim*c2 + 2*beta6*bpr*c;
  D2Vw.Init(0);
  D2Vw(3) = 2*((betaprim2+beta*betasecn)*(2*b + 1) 
	       + 8*beta*betaprim*bpr 
	       + bsc*beta2);
  D2Vw(5) = 4*(b2 + 2*c)*(3*beta2*betaprim2 + beta3*betasecn)
             + 16*beta3*betaprim*bpr*(b+11)
             + 2*beta4*(bsc*(b+1)+bpr2); 
  D2Vw(7) = 6 * c2 * (5*beta4*betaprim2+beta5*betasecn)
            + 24*beta5*betaprim*bpr*c
	    + 2*beta6*(bsc*c + bpr2);
 

  // Calcul des poles
  Px = B * Vx;
  Py = B * Vy;
  W.Multiply(B, Vw);
  DPx = B * DVx;
  DPy = B * DVy;
  DW.Multiply(B, DVw);
  D2Px = B * D2Vx;
  D2Py = B * D2Vy;
  D2W.Multiply(B, D2Vw);

  gp_XYZ P, DP, D2P;
  Standard_Real wi, dwi;
  for (ii=1; ii<=Ordre; ii++) {
     wi = W(ii);
     dwi = DW(ii);
     P.SetCoord(Px(ii)/wi, Py(ii)/wi, 0);
     DP.SetCoord(DPx(ii)/wi, DPy(ii)/wi, 0);
   
     D2P.SetCoord(D2Px(ii)/wi, D2Py(ii)/wi, 0);
     D2P -= 2*(dwi/wi)*DP;
     D2P += (2*Pow(dwi/wi, 2) - D2W(ii)/wi)*P;
     DP -= (DW(ii)/wi)*P;  

     Poles(ii).ChangeCoord() = M*P + Center.XYZ();
     auxyz.SetLinearForm(1, MPrim*P, 
			 1, M*DP, 
			 DCenter.XYZ());  
     DPoles(ii).SetXYZ(auxyz);
     P  *= MSecn;
     DP *= MPrim;
     D2P*= M;
     auxyz.SetLinearForm(1, P, 
			 2, DP, 
			 1, D2P,
			 D2Center.XYZ());
     D2Poles(ii).SetXYZ(auxyz);
     Weights(ii) = wi;
     DWeights(ii) = dwi; 
     D2Weights(ii) = D2W(ii);     
  }
}
Commit	Line	Data
b311480e	1	// Created on: 1997-08-06
	2	// Created by: Philippe MANGIN
	3	// Copyright (c) 1997-1999 Matra Datavision
973c2be1	4	// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e	5	//
973c2be1	6	// This file is part of Open CASCADE Technology software library.
b311480e	7	//
d5f74e42	8	// This library is free software; you can redistribute it and/or modify it under
d5f74e42	9	// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1	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.
b311480e	13	//
973c2be1	14	// Alternatively, this file may be used under the terms of Open CASCADE
973c2be1	15	// commercial license or contractual agreement.
7fd59977	16
7fd59977	17
7fd59977	18	#include <Convert_CompPolynomialToPoles.hxx>
42cf5bc1	19	#include <GeomFill_QuasiAngularConvertor.hxx>
	20	#include <gp_Ax1.hxx>
	21	#include <gp_Mat.hxx>
	22	#include <gp_Pnt.hxx>
	23	#include <gp_Vec.hxx>
	24	#include <PLib.hxx>
	25	#include <StdFail_NotDone.hxx>
7fd59977	26	#include <TColStd_Array1OfReal.hxx>
7fd59977	27	#include <TColStd_HArray1OfInteger.hxx>
7fd59977	28	#include <TColStd_HArray1OfReal.hxx>
42cf5bc1	29	#include <TColStd_HArray2OfReal.hxx>
7fd59977	30
	31	#define NullAngle 1.e-6
	32
	33	// QuasiAngular is rational definition of Cos(theta(t) and sin(theta)
	34	// on [-alpha, +alpha] with
	35	// 2 2
	36	// U - V
	37	// cos (theta(t)) = ----------
	38	// 2 2
	39	// U + V
7fd59977	40
	41	// 2 * U*V
	42	// sin (theta(t)) = ----------
	43	// 2 2
	44	// U + V
	45	// 3
	46	// V(t) = t + c t
	47	// 2
	48	// U(t) = 1 + b t
	49	// 1
	50	// c = --- + b
	51	// 3
7fd59977	52	// -1 gamma
	53	// b =--------- + -----------------------
	54	// 2
	55	// gamma 3*(tang gamma - gamma)
7fd59977	56	// with gamma = alpha / 2
	57
	58
	59	GeomFill_QuasiAngularConvertor::GeomFill_QuasiAngularConvertor():
	60	myinit(Standard_False),
	61	B(1, 7, 1, 7),
	62	Px(1, 7), Py(1, 7),
	63	W(1,7), Vx(1, 7), Vy(1, 7), Vw(1,7)
	64	{
	65	}
	66
	67	Standard_Boolean GeomFill_QuasiAngularConvertor::Initialized() const
	68	{
	69	return myinit;
	70	}
	71
	72	void GeomFill_QuasiAngularConvertor::Init()
	73	{
	74	if (myinit) return; //On n'initialise qu'une fois
	75	Standard_Integer ii, jj, Ordre=7;
	76	Standard_Real terme;
	77	TColStd_Array1OfReal Coeffs(1, Ordre*Ordre), TrueInter(1,2), Inter(1,2);
	78	Handle(TColStd_HArray2OfReal)
	79	Poles1d = new (TColStd_HArray2OfReal) (1, Ordre, 1, Ordre);
	80
	81	//Calcul de B
	82	Inter.SetValue(1, -1);
	83	Inter.SetValue(2, 1);
	84	TrueInter.SetValue(1, -1);
	85	TrueInter.SetValue(2, 1);
	86
	87	Coeffs.Init(0);
	88	for (ii=1; ii<=Ordre; ii++) { Coeffs.SetValue(ii+(ii-1)*Ordre, 1); }
	89
	90	//Convertion
	91	Convert_CompPolynomialToPoles
	92	AConverter(Ordre, Ordre-1, Ordre-1,
	93	Coeffs,
	94	Inter,
	95	TrueInter);
	96	AConverter.Poles(Poles1d);
	97
	98	for (jj=1; jj<=Ordre; jj++) {
	99	for (ii=1; ii<=Ordre; ii++) {
	100	terme = Poles1d->Value(ii,jj);
	101	if (Abs(terme-1) < 1.e-9) terme = 1 ; //petite retouche
	102	if (Abs(terme+1) < 1.e-9) terme = -1;
	103	B(ii, jj) = terme;
	104	}
	105	}
	106
	107	// Init des polynomes
	108	Vx.Init(0);
	109	Vx(1) = 1;
	110	Vy.Init(0);
	111	Vy(2) = 2;
	112	Vw.Init(0);
	113	Vw(1) = 1;
	114	myinit = Standard_True;
	115	}
	116
	117	void GeomFill_QuasiAngularConvertor::Section(const gp_Pnt& FirstPnt,
	118	const gp_Pnt& Center,
	119	const gp_Vec& Dir,
120	const Standard_Real Angle,
121	TColgp_Array1OfPnt& Poles,
122	TColStd_Array1OfReal& Weights)
123	{
124	Standard_Real b, b2, c, c2,tan_b;
125	Standard_Integer ii;
126	Standard_Real beta, beta2, beta3, beta4, beta5, beta6, wi;
127	gp_XYZ aux;
128	gp_Mat Rot;
129	// Calcul de la transformation
130	gp_Vec V1(Center, FirstPnt), V2;
131	Rot.SetRotation(Dir.XYZ(), Angle/2);
132	aux = V1.XYZ();
133	aux *= Rot;
134	V1.SetXYZ(aux);
135	V2 = Dir^V1;
136
137	gp_Mat M(V1.X(), V2.X(), 0,
138	V1.Y(), V2.Y(), 0,
139	V1.Z(), V2.Z(), 0);
7fd59977	140
	141	// Calcul des coeffs -----------
	142	beta = Angle/4;
	143	beta2 = beta * beta;
	144	beta3 = beta * beta2;
	145	beta4 = beta2*beta2;
	146	beta5 = beta3*beta2;
	147	beta6 = beta3*beta3;
	148
c6541a0c	149	if ((M_PI/2 - beta)> NullAngle) {
7fd59977	150	if (Abs(beta) < NullAngle) {
	151	Standard_Real cf = 2.0/(357);
	152	b = - (0.2+cfbeta2) / (1+ 0.2beta2);
	153	// b = beta5 / cf;
	154	}
	155	else {
	156	tan_b = Tan(beta);
	157	b = - 1.0e0 / beta2;
	158	b += beta / (3*(tan_b - beta));
	159	}
	160	}
	161	else b = ((Standard_Real) -1)/beta2;
	162	c = ((Standard_Real) 1)/ 3 + b;
	163	b2 = b*b;
	164	c2 = c*c;
	165
	166	// X = UU - VV
	167	Vx(3) = beta2(2b - 1);
	168	Vx(5) = beta4(b2 - 2c);
	169	Vx(7) = -beta6*c2;
	170
	171	//Y = 2UV
	172	Vy(2) = 2*beta;
	173	Vy(4) = beta32(c+b);
	174	Vy(6) = 2beta5b*c;
	175
	176	// W = UU + VV
	177	Vw(3) = beta2(1 + 2b);
	178	Vw(5) = beta4(2c + b2);
	179	Vw(7) = beta6*c2;
	180
	181	// Calculs des poles
	182	Px.Multiply(B, Vx);
	183	Py.Multiply(B, Vy);
	184	W.Multiply(B, Vw);
	185
	186	// Transfo
	187	gp_XYZ pnt;
	188	for (ii=1; ii<=7; ii++) {
	189	wi = W(ii);
	190	pnt.SetCoord(Px(ii)/wi, Py(ii)/wi, 0);
	191	pnt *= M;
	192	pnt += Center.XYZ();
	193	Poles(ii).ChangeCoord() = pnt;
	194	Weights(ii) = wi;
	195	}
	196	}
	197
	198	void GeomFill_QuasiAngularConvertor::Section(const gp_Pnt& FirstPnt,
	199	const gp_Vec& DFirstPnt,
	200	const gp_Pnt& Center,
	201	const gp_Vec& DCenter,
	202	const gp_Vec& Dir,
	203	const gp_Vec& DDir,
	204	const Standard_Real Angle,
	205	const Standard_Real DAngle,
	206	TColgp_Array1OfPnt& Poles,
	207	TColgp_Array1OfVec& DPoles,
	208	TColStd_Array1OfReal& Weights,
	209	TColStd_Array1OfReal& DWeights)
	210	{
	211	Standard_Integer Ordre = 7;
	212	math_Vector DVx(1, Ordre), DVy(1, Ordre), DVw(1, Ordre),
	213	DPx(1, Ordre), DPy(1, Ordre), DW(1, Ordre);
214	Standard_Real b, b2, c, c2,tan_b;
215	Standard_Real bpr, dtan_b;
216	Standard_Integer ii;
217	Standard_Real beta, beta2, beta3, beta4, beta5, beta6, betaprim;
218	gp_Vec V1(Center, FirstPnt), V1Prim, V2;
219
220	// Calcul des transformations
221	gp_XYZ aux;
222	Standard_Real Sina, Cosa;
223	gp_Mat Rot, RotPrim, D, DPrim;
224	// La rotation s'ecrit I + sin(Ang) * D + (1. - cos(Ang)) * D*D
225	// ou D est l'application x -> Dir ^ x
226	Rot.SetRotation(Dir.XYZ(), Angle/2);
227	// La derive s'ecrit donc :
228	// AngPrim * (sin(Ang)DD + cos(Ang)*D)
229	// + sin(Ang)DPrim + (1. - cos(Ang)) (DPrimD + DDPrim)
230	Sina = Sin(Angle/2);
231	Cosa = Cos(Angle/2);
232	D.SetCross(Dir.XYZ());
233	DPrim.SetCross(DDir.XYZ());
234
235	RotPrim = (D.Powered(2)).Multiplied(Sina);
236	RotPrim += D.Multiplied(Cosa);
237	RotPrim *= DAngle/2;
238	RotPrim += DPrim.Multiplied(Sina);
239	RotPrim += ((DPrim.Multiplied(D)).Added(D.Multiplied(DPrim))).Multiplied(1-Cosa);
240
241	aux = (DFirstPnt - DCenter).XYZ().Multiplied(Rot);
242	aux += V1.XYZ().Multiplied(RotPrim);
243	V1Prim.SetXYZ(aux);
244	aux = V1.XYZ();
245	aux.Multiply(Rot);
246	V1.SetXYZ(aux);
247	V2 = Dir^V1;
248	gp_Mat M (V1.X(), V2.X(), 0,
249	V1.Y(), V2.Y(), 0,
250	V1.Z(), V2.Z(), 0);
251	V2 = (DDir.Crossed(V1)).Added(Dir.Crossed(V1Prim));
252	gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
253	V1Prim.Y(), V2.Y(), 0,
254	V1Prim.Z(), V2.Z(), 0);
255
256	// Calcul des constante -----------
257	beta = Angle/4;
258	betaprim = DAngle/4;
259	beta2 = beta * beta;
260	beta3 = beta * beta2;
261	beta4 = beta2*beta2;
262	beta5 = beta3*beta2;
263	beta6 = beta3*beta3;
264
265	if (Abs(beta) < NullAngle) {
266	// On calcul b par D.L
267	Standard_Real cf = 2.0/(357);
268	Standard_Real Num, Denom;
269	Num = 0.2 + cf*beta2;
270	Denom = 1+0.2*beta2;
271	b = - Num/Denom;
272	bpr = -2betabetaprim(cfDenom - 0.2Num)/(DenomDenom);
273	}
274	else {
275	b = ((Standard_Real) -1)/beta2;
276	bpr = (2*betaprim) / beta3;
c6541a0c	277	if ((M_PI/2 - beta)> NullAngle) {
7fd59977	278	tan_b = Tan(beta);
	279	dtan_b = betaprim * (1 + tan_b*tan_b);
	280	b2 = tan_b - beta;
	281	b += beta / (3*b2);
	282	bpr += (betaprimtan_b - betadtan_b) / (3b2b2);
	283	}
	284	}
	285
	286	c = ((Standard_Real) 1)/ 3 + b;
	287	b2 = b*b;
	288	c2 = c*c;
	289
	290
	291	// X = UU - VV
	292	Vx(3) = beta2(2b - 1);
	293	Vx(5) = beta4(b2 - 2c);
	294	Vx(7) = -beta6*c2;
	295	DVx.Init(0);
	296	DVx(3) = 2(betabetaprim(2b - 1) + bpr*beta2);
	297	DVx(5) = 4beta3betaprim(b2 - 2c) + 2beta4bpr*(b-1);
	298	DVx(7) = - 6beta5betaprimc2 - 2beta6bprc;
	299
	300	//Y = 2UV
	301	Vy(2) = 2*beta;
	302	Vy(4) = beta32(c+b);
	303	Vy(6) = 2beta5b*c;
	304	DVy.Init(0);
	305	DVy(2) = 2*betaprim;
	306	DVy(4) = 6beta2betaprim(b+c) + 4beta3*bpr;
	307	DVy(6) = 10beta4betaprimbc + 2beta5bpr*(b+c);
	308
	309	// W = UU + VV
	310	Vw(3) = beta2(1 + 2b);
	311	Vw(5) = beta4(2c + b2);
	312	Vw(7) = beta6*c2;
	313	DVw.Init(0);
	314	// DVw(3) = 2(betabetaprim(1 + 2b) + beta2*bpr);
	315	DVw(3) = 2beta(betaprim(1 + 2b) + beta*bpr);
	316	// DVw(5) = 4beta3betaprim(2c + b2) + 2beta4bpr*(b+1);
	317	DVw(5) = 2beta3(2betaprim(2c + b2) + betabpr*(b+1));
	318	// DVw(7) = 6beta5betaprimc2 + 2beta6bprc;
	319	DVw(7) = 2beta5c(3betaprimc + betabpr);
	320
	321
	322
	323	// Calcul des poles
	324	Px.Multiply(B, Vx);
	325	Py.Multiply(B, Vy);
	326	W.Multiply(B, Vw);
	327	DPx.Multiply(B, DVx);
	328	DPy.Multiply(B, DVy);
	329	DW.Multiply(B, DVw);
	330
	331	gp_XYZ P, DP;
	332	Standard_Real wi;
	333
	334	for (ii=1; ii<=Ordre; ii++) {
	335	wi = W(ii);
	336	P.SetCoord(Px(ii)/wi, Py(ii)/wi, 0);
	337	DP.SetCoord(DPx(ii)/wi, DPy(ii)/wi, 0);
	338	DP -= (DW(ii)/wi)*P;
	339
	340	Poles(ii).ChangeCoord() = M*P + Center.XYZ();
	341	P *= MPrim;
342	DP *= M;
343	aux.SetLinearForm(1, P, 1, DP, DCenter.XYZ());
344	DPoles(ii).SetXYZ(aux);
345	Weights(ii) = wi;
346	DWeights(ii) = DW(ii);
347	}
348	}
349
350	void GeomFill_QuasiAngularConvertor::Section(const gp_Pnt& FirstPnt,
351	const gp_Vec& DFirstPnt,
352	const gp_Vec& D2FirstPnt,
353	const gp_Pnt& Center,
354	const gp_Vec& DCenter,
355	const gp_Vec& D2Center,
356	const gp_Vec& Dir,
357	const gp_Vec& DDir,
358	const gp_Vec& D2Dir,
359	const Standard_Real Angle,
360	const Standard_Real DAngle,
361	const Standard_Real D2Angle,
362	TColgp_Array1OfPnt& Poles,
363	TColgp_Array1OfVec& DPoles,
364	TColgp_Array1OfVec& D2Poles,
365	TColStd_Array1OfReal& Weights,
366	TColStd_Array1OfReal& DWeights,
367	TColStd_Array1OfReal& D2Weights)
368	{
369	Standard_Integer Ordre = 7;
370	math_Vector DVx(1, Ordre), DVy(1, Ordre), DVw(1, Ordre),
371	D2Vx(1, Ordre), D2Vy(1, Ordre), D2Vw(1, Ordre);
372	math_Vector DPx(1, Ordre), DPy(1, Ordre), DW(1, Ordre),
373	D2Px(1, Ordre), D2Py(1, Ordre), D2W(1, Ordre);
374
375	Standard_Integer ii;
51740958	376	Standard_Real daux, b, b2, c, c2, bpr, bsc;
7fd59977	377	gp_Vec V1(Center, FirstPnt), V1Prim, V1Secn, V2;
	378
	379	// Calcul des transformations
	380	gp_XYZ auxyz;
	381	Standard_Real Sina, Cosa;
	382	gp_Mat Rot, RotPrim, RotSecn, D, DPrim, DSecn, DDP, Maux;
	383	// La rotation s'ecrit I + sin(Ang) * D + (1. - cos(Ang)) * D*D
	384	// ou D est l'application x -> Dir ^ x
	385	Rot.SetRotation(Dir.XYZ(), Angle/2);
	386	// La derive s'ecrit donc :
	387	// AngPrim * (sin(Ang)DD + cos(Ang)*D)
	388	// + sin(Ang)DPrim + (1. - cos(Ang)) (DPrimD + DDPrim)
	389	Sina = Sin(Angle/2);
	390	Cosa = Cos(Angle/2);
	391	D.SetCross(Dir.XYZ());
	392	DPrim.SetCross(DDir.XYZ());
	393	DSecn.SetCross(D2Dir.XYZ());
	394
	395	DDP = (DPrim.Multiplied(D)).Added(D.Multiplied(DPrim));
	396	RotPrim = (D.Powered(2)).Multiplied(Sina);
	397	RotPrim += D.Multiplied(Cosa);
	398	RotPrim *= DAngle/2;
	399	RotPrim += DPrim.Multiplied(Sina);
	400	RotPrim += DDP.Multiplied(1-Cosa);
	401
	402	RotSecn = (D.Powered(2)).Multiplied(Sina);
	403	RotSecn += D.Multiplied(Cosa);
	404	RotSecn *= D2Angle/2;
	405	Maux = (D.Powered(2)).Multiplied(Cosa);
	406	Maux -= D.Multiplied(Sina);
	407	Maux *= DAngle/2;
	408	Maux += DDP.Multiplied(2*Sina);
	409	Maux += DPrim.Multiplied(2*Cosa);
	410	Maux *= DAngle/2;
	411	RotSecn += Maux;
	412	Maux = (DSecn.Multiplied(D)).Added(D.Multiplied(DSecn));
	413	Maux += (DPrim.Powered(2)).Multiplied(2);
	414	Maux *= 1 - Cosa;
	415	Maux += DSecn.Multiplied(Sina);
	416	RotSecn += Maux;
	417
	418	V1Prim = DFirstPnt - DCenter;
	419	auxyz = (D2FirstPnt - D2Center).XYZ().Multiplied(Rot);
	420	auxyz += 2*(V1Prim.XYZ().Multiplied(RotPrim));
	421	auxyz += V1.XYZ().Multiplied(RotSecn);
	422	V1Secn.SetXYZ(auxyz);
	423	auxyz = V1Prim.XYZ().Multiplied(Rot);
	424	auxyz += V1.XYZ().Multiplied(RotPrim);
	425	V1Prim.SetXYZ(auxyz);
	426	auxyz = V1.XYZ();
	427	auxyz.Multiply(Rot);
	428	V1.SetXYZ(auxyz);
	429	V2 = Dir^V1;
	430
	431	gp_Mat M (V1.X(), V2.X(), 0,
	432	V1.Y(), V2.Y(), 0,
	433	V1.Z(), V2.Z(), 0);
	434	V2 = (DDir.Crossed(V1)).Added(Dir.Crossed(V1Prim));
	435	gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
	436	V1Prim.Y(), V2.Y(), 0,
	437	V1Prim.Z(), V2.Z(), 0);
	438
	439	V2 = DDir^V1Prim;
	440	V2 *= 2;
441	V2 += (D2Dir.Crossed(V1)).Added(Dir.Crossed(V1Secn));
442	gp_Mat MSecn (V1Secn.X(), V2.X(), 0,
443	V1Secn.Y(), V2.Y(), 0,
444	V1Secn.Z(), V2.Z(), 0);
445
446
447	// Calcul des coeff -----------
448	Standard_Real tan_b, dtan_b, d2tan_b;
449	Standard_Real beta, beta2, beta3, beta4, beta5, beta6, betaprim, betasecn;
450	Standard_Real betaprim2, bpr2;
451	beta = Angle/4;
452	betaprim = DAngle/4;
453	betasecn = D2Angle/4;
454	beta2 = beta * beta;
455	beta3 = beta * beta2;
456	beta4 = beta2*beta2;
457	beta5 = beta3*beta2;
458	beta6 = beta3*beta3;
459	betaprim2 = betaprim * betaprim;
460
461	if (Abs(beta) < NullAngle) {
462	// On calcul b par D.L
463	Standard_Real cf =-2.0/21;
464	Standard_Real Num, Denom, aux;
465	Num = 0.2 + cf*beta2;
466	Denom = 1+0.2*beta2;
467	aux = (cfDenom - 0.2Num)/(Denom*Denom);
468	b = - Num/Denom;
469	bpr = -2betabetaprim*aux;
470	bsc = 2aux(betaprim2 + betabetasecn - 2beta*betaprim2);
471
472	}
473	else {
474	b = ((Standard_Real) -1)/beta2;
475	bpr = (2*betaprim) / beta3;
476	bsc = (2betasecn - 6betaprim*(betaprim/beta)) / beta3;
c6541a0c	477	if ((M_PI/2 - beta)> NullAngle) {
7fd59977	478	tan_b = Tan(beta);
	479	dtan_b = betaprim * (1 + tan_b*tan_b);
	480	d2tan_b = betasecn * (1 + tan_b*tan_b)
	481	+ 2betaprim tan_b * dtan_b;
	482	b2 = tan_b - beta;
	483	b += beta / (3*b2);
51740958	484	Standard_Real aux = betaprimtan_b - betadtan_b;
7fd59977	485	bpr += aux / (3b2b2);
	486	daux = betasecntan_b - betad2tan_b;
	487	bsc += (daux - 2auxbetaprimtan_btan_b/b2)/(3b2b2);
	488	}
	489	}
	490
	491	c = ((Standard_Real) 1)/ 3 + b;
	492	b2 = b*b;
	493	c2 = c*c;
	494	bpr2 = bpr * bpr;
	495
	496
	497	// X = UU - VV
	498	Vx(3) = beta2(2b - 1);
	499	Vx(5) = beta4(b2 - 2c);
	500	Vx(7) = -beta6*c2;
	501	DVx.Init(0);
	502	DVx(3) = 2(betabetaprim(2b - 1) + bpr*beta2);
	503	DVx(5) = 4beta3betaprim(b2 - 2c) + 2beta4bpr*(b-1);
	504	DVx(7) = - 6beta5betaprimc2 - 2beta6bprc;
	505	D2Vx.Init(0);
	506	D2Vx(3) = 2((betaprim2+betabetasecn)(2b - 1)
	507	+ 8betabetaprim*bpr
	508	+ bsc*beta2);
	509	D2Vx(5) = 4(b2 - 2c)(3beta2betaprim2 + beta3betasecn)
	510	+ 16beta3betaprimbpr(b-1)
	511	+ 2beta4(bsc*(b-1)+bpr2);
	512	D2Vx(7) = - 6 * c2 * (5beta4betaprim2+beta5*betasecn)
	513	- 24beta5betaprimbprc
	514	- 2beta6(bsc*c + bpr2);
	515
	516	//Y = 2UV
	517	Vy(2) = 2*beta;
	518	Vy(4) = beta32(c+b);
	519	Vy(6) = 2beta5b*c;
	520	DVy.Init(0);
	521	DVy(2) = 2*betaprim;
	522	DVy(4) = 6beta2betaprim(b+c) + 4beta3*bpr;
	523	DVy(6) = 10beta4betaprimbc + 2beta5bpr*(b+c);
	524	D2Vy.Init(0);
	525	D2Vy(2) = 2*betasecn;
	526	D2Vy(4) = 6(b+c)(2betabetaprim2 + beta2*betasecn)
	527	+ 24beta2betaprimbpr(b+c)
	528	+ 4beta3bsc;
	529	D2Vy(6) = 10bc(4beta3betaprim2 + beta4betasecn)
	530	+ 40 * beta4betaprimbpr*(b+c)
	531	+ 2beta5(bsc(b+c)+ 2bpr2);
	532
	533	// W = UU + VV
	534	Vw(3) = beta2(1 + 2b);
	535	Vw(5) = beta4(2c + b2);
	536	Vw(7) = beta6*c2;
	537	DVw.Init(0);
	538	DVw(3) = 2(betabetaprim(1 + 2b) + beta2*bpr);
	539	DVw(5) = 4beta3betaprim(2c + b2) + 2beta4bpr*(b+1);
	540	DVw(7) = 6beta5betaprimc2 + 2beta6bprc;
	541	D2Vw.Init(0);
	542	D2Vw(3) = 2((betaprim2+betabetasecn)(2b + 1)
	543	+ 8betabetaprim*bpr
	544	+ bsc*beta2);
	545	D2Vw(5) = 4(b2 + 2c)(3beta2betaprim2 + beta3betasecn)
	546	+ 16beta3betaprimbpr(b+11)
	547	+ 2beta4(bsc*(b+1)+bpr2);
	548	D2Vw(7) = 6 * c2 * (5beta4betaprim2+beta5*betasecn)
549	+ 24beta5betaprimbprc
550	+ 2beta6(bsc*c + bpr2);
551
552
553	// Calcul des poles
554	Px = B * Vx;
555	Py = B * Vy;
556	W.Multiply(B, Vw);
557	DPx = B * DVx;
558	DPy = B * DVy;
559	DW.Multiply(B, DVw);
560	D2Px = B * D2Vx;
561	D2Py = B * D2Vy;
562	D2W.Multiply(B, D2Vw);
563
564	gp_XYZ P, DP, D2P;
565	Standard_Real wi, dwi;
566	for (ii=1; ii<=Ordre; ii++) {
567	wi = W(ii);
568	dwi = DW(ii);
569	P.SetCoord(Px(ii)/wi, Py(ii)/wi, 0);
570	DP.SetCoord(DPx(ii)/wi, DPy(ii)/wi, 0);
571
572	D2P.SetCoord(D2Px(ii)/wi, D2Py(ii)/wi, 0);
573	D2P -= 2(dwi/wi)DP;
574	D2P += (2Pow(dwi/wi, 2) - D2W(ii)/wi)P;
575	DP -= (DW(ii)/wi)*P;
576
577	Poles(ii).ChangeCoord() = M*P + Center.XYZ();
578	auxyz.SetLinearForm(1, MPrim*P,
579	1, M*DP,
580	DCenter.XYZ());
581	DPoles(ii).SetXYZ(auxyz);
582	P *= MSecn;
583	DP *= MPrim;
584	D2P*= M;
585	auxyz.SetLinearForm(1, P,
586	2, DP,
587	1, D2P,
588	D2Center.XYZ());
589	D2Poles(ii).SetXYZ(auxyz);
590	Weights(ii) = wi;
591	DWeights(ii) = dwi;
592	D2Weights(ii) = D2W(ii);
593	}
594	}