diff --git a/include/Utils/qem.h b/include/Utils/qem.h
index 823e8dcd0a42ca697879a495e0bdb491acbc7291..606fa249f8780f27dcd1c87061ea82c537c8104e 100644
--- a/include/Utils/qem.h
+++ b/include/Utils/qem.h
@@ -40,8 +40,6 @@
 // Eigen includes
 #include <Eigen/Dense>
 
-#define CONST_VAL -5212368.54127 // random value
-
 /*! \namespace CGoGN
  * \brief namespace for all elements composing the CGoGN library
  */
@@ -563,12 +561,22 @@ public:
 	static std::vector<VEC3> coefsFromTensors(Geom::Tensor3d* T) ;
 
 	/*!
-	 * \brief method to complete a symmetric tensor that was
-	 * only filled in its first half
+	 * \brief method to complete a symmetric matrix that was
+	 * only filled in its first half (line >= column)
 	 *
-	 * \param T the tensor to fill
+	 * \param M the matrix to fill
 	 */
-	static void fillTensor(Geom::Tensor3d& T) ;
+	static void fillSymmetricMatrix(Eigen::MatrixXd& M) ;
+
+	/*!
+	 * \brief method to rotate a tensor representing a polynomial light field
+	 *
+	 * \param T the tensor representing a polynomial light field
+	 * \param R the 3x3 matrix representing a rotation in (u,v,1)-space
+	 *
+	 * \return a new rotated tensor representing a polynomial light field.
+	 */
+	static Geom::Tensor3d rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R) ;
 
 private:
 	// Double computation is crucial for stability
@@ -595,7 +603,7 @@ private:
 	 */
 	bool optimize(std::vector<VEC3>& coefs) const ;
 
-	Geom::Tensor3d rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R) ;
+//	Geom::Tensor3d rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R) ;
 
 	/*!
 	 * \brief method to build a rotate matrix (rotation in tangent plane)
@@ -639,6 +647,9 @@ private:
 	 * \return the integral of product of monomes
 	 */
 	Eigen::MatrixXd buildIntegralMatrix_C(const REAL& alpha, unsigned int size) ;
+
+	Eigen::MatrixXd buildLowerLeftIntegralMatrix_A(const REAL& alpha, unsigned int size) ;
+	Eigen::MatrixXd buildLowerLeftIntegralMatrix_C(const REAL& alpha, unsigned int size) ;
 } ;
 
 } // Utils
diff --git a/include/Utils/qem.hpp b/include/Utils/qem.hpp
index e42b86735253527afc23239e79426091378cf1b3..536ec7fb5c86e55f960ea5767e2c75efe0264d93 100644
--- a/include/Utils/qem.hpp
+++ b/include/Utils/qem.hpp
@@ -344,18 +344,6 @@ QuadricHF<REAL>::QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REA
 	Eigen::MatrixXd integ_b = buildIntegralMatrix_B(alpha,nbcoefs) ; // Parameterized integral matrix b
 	Eigen::MatrixXd integ_c = buildIntegralMatrix_C(alpha,nbcoefs) ; // Parameterized integral matrix c
 
-	//		for (unsigned int i = 0 ; i < nbcoefs ; ++i)
-	//			for (unsigned int j = 0 ; j < nbcoefs ; ++j)
-	//				std::cout << "(" << i << "," << j << ")=" << m_A(i,j) << std::endl ;
-	//
-	//		for (unsigned int i = 0 ; i < nbcoefs ; ++i)
-	//			for (unsigned int j = 0 ; j < nbcoefs ; ++j)
-	//				std::cout << "(" << i << "," << j << ")=" << integ_b(i,j) << std::endl ;
-	//
-	//		for (unsigned int i = 0 ; i < nbcoefs ; ++i)
-	//			for (unsigned int j = 0 ; j < nbcoefs ; ++j)
-	//				std::cout << "(" << i << "," << j << ")=" << integ_c(i,j) << std::endl ;
-
 	for (unsigned int c = 0 ; c < 3 ; ++c)
 	{
 		Eigen::VectorXd v ;	v.resize(nbcoefs) ;
@@ -369,11 +357,7 @@ QuadricHF<REAL>::QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REA
 
 template <typename REAL>
 QuadricHF<REAL>::~QuadricHF()
-{
-	//delete m_A ;
-	//for (unsigned int c = 0 ; c < 3 ; ++c)
-	//	delete m_b[c] ;
-}
+{}
 
 template <typename REAL>
 void
@@ -462,17 +446,6 @@ QuadricHF<REAL>::operator /= (const REAL& v)
 	return *this ;
 }
 
-
-//	REAL
-// QuadricHF<REAL>::operator() (const VECNp& v) const
-//	{
-//		VECN hv ;
-//		for (unsigned int i = 0 ; i < 3 ; ++i)
-//			hv[i] = v[i] ;
-//
-//		return evaluate(v) ;
-//}
-
 template <typename REAL>
 REAL
 QuadricHF<REAL>::operator() (const std::vector<VEC3>& coefs) const
@@ -484,16 +457,6 @@ template <typename REAL>
 bool
 QuadricHF<REAL>::findOptimizedCoefs(std::vector<VEC3>& coefs)
 {
-	//		unsigned int nbcoefs = 6 ;
-	//		for (unsigned int i = 0 ; i < nbcoefs ; ++i)
-	//			for (unsigned int j = 0 ; j < nbcoefs ; ++j)
-	//				std::cout << "A(" << i << "," << j << ")=" << m_A(i,j) << std::endl ;
-	//
-	//		for (unsigned int i = 0 ; i < nbcoefs ; ++i)
-	//				std::cout << "b(" << i << ")=" << m_b[0][i] << std::endl ;
-	//
-	//		std::cout << "c=" << m_c[0] << std::endl ;
-
 	return optimize(coefs) ;
 }
 
@@ -523,7 +486,7 @@ QuadricHF<REAL>::optimize(std::vector<VEC3>& coefs) const
 {
 	coefs.resize(m_b[0].size()) ;
 
-	if (fabs(m_A.determinant()) < 1e-6 )
+	if (fabs(m_A.determinant()) < 1e-10 )
 		return false ;
 
 	Eigen::MatrixXd Ainv = m_A.inverse() ;
@@ -578,171 +541,183 @@ QuadricHF<REAL>::rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R)
 	return Tp ;
 }
 
+template <typename REAL>
+void
+QuadricHF<REAL>::fillSymmetricMatrix(Eigen::MatrixXd& M)
+{
+	assert(M.cols() == M.rows() || !"QuadricHF<REAL>::fillSymmetricMatrix: matrix to fill should be a square matrix") ;
+	for (unsigned int c = 0 ; c < M.cols() ; ++c)
+		for (unsigned int l = 0 ; l < c ; ++l)
+			M( l, c) = M (c, l) ;
+}
+
 template <typename REAL>
 Eigen::MatrixXd
 QuadricHF<REAL>::buildIntegralMatrix_A(const REAL& alpha, unsigned int size)
+{
+	Eigen::MatrixXd res = buildLowerLeftIntegralMatrix_A(alpha,size) ;
+	fillSymmetricMatrix(res) ;
+
+	return res ;
+}
+
+template <typename REAL>
+Eigen::MatrixXd
+QuadricHF<REAL>::buildIntegralMatrix_C(const REAL& alpha, unsigned int size)
+{
+	Eigen::MatrixXd res = buildLowerLeftIntegralMatrix_C(alpha,size) ;
+	fillSymmetricMatrix(res) ;
+
+	return res ;
+}
+
+
+template <typename REAL>
+Eigen::MatrixXd
+QuadricHF<REAL>::buildLowerLeftIntegralMatrix_A(const REAL& alpha, unsigned int size)
 {
 	Eigen::MatrixXd A(size,size) ;
 
-	A( 0 , 0 ) =  2*(M_PI-alpha) ;
-	A( 0 , 1 ) =  M_PI*sin(alpha)/2.0 ;
-	A( 0 , 2 ) =  0 ;
-	A( 0 , 3 ) =  0 ;
-	A( 0 , 4 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ;
-	A( 0 , 5 ) =  2.0*(M_PI-alpha)/3.0 ;
-	A( 1 , 0 ) =  M_PI*sin(alpha)/2.0 ;
-	A( 1 , 1 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ;
-	A( 1 , 2 ) =  0 ;
-	A( 1 , 3 ) =  0 ;
-	A( 1 , 4 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/8.0 ;
-	A( 1 , 5 ) =  M_PI*sin(alpha)/8.0 ;
-	A( 2 , 0 ) =  0 ;
-	A( 2 , 1 ) =  0 ;
-	A( 2 , 2 ) =  2.0*(M_PI-alpha)/3.0 ;
-	A( 2 , 3 ) =  M_PI*sin(alpha)/8.0 ;
-	A( 2 , 4 ) =  0 ;
-	A( 2 , 5 ) =  0 ;
-	A( 3 , 0 ) =  0 ;
-	A( 3 , 1 ) =  0 ;
-	A( 3 , 2 ) =  M_PI*sin(alpha)/8.0 ;
-	A( 3 , 3 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
-	A( 3 , 4 ) =  0 ;
-	A( 3 , 5 ) =  0 ;
-	A( 4 , 0 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ;
-	A( 4 , 1 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/8.0 ;
-	A( 4 , 2 ) =  0 ;
-	A( 4 , 3 ) =  0 ;
-	A( 4 , 4 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ;
-	A( 4 , 5 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
-	A( 5 , 0 ) =  2.0*(M_PI-alpha)/3.0 ;
-	A( 5 , 1 ) =  M_PI*sin(alpha)/8.0 ;
-	A( 5 , 2 ) =  0 ;
-	A( 5 , 3 ) =  0 ;
-	A( 5 , 4 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
-	A( 5 , 5 ) =  2.0*(M_PI-alpha)/5.0 ;
+	A( 0 , 0 ) = 2*(M_PI-alpha) ;
+
+	A( 1 , 0 ) = M_PI*sin(alpha)/2.0 ;
+	A( 1 , 1 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ;
+
+	A( 2 , 0 ) = 0 ;
+	A( 2 , 1 ) = 0 ;
+	A( 2 , 2 ) = 2.0*(M_PI-alpha)/3.0 ;
+
+	A( 3 , 0 ) = 0 ;
+	A( 3 , 1 ) = 0 ;
+	A( 3 , 2 ) = M_PI*sin(alpha)/8.0 ;
+	A( 3 , 3 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
+
+	A( 4 , 0 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ;
+	A( 4 , 1 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/8.0 ;
+	A( 4 , 2 ) = 0 ;
+	A( 4 , 3 ) = 0 ;
+	A( 4 , 4 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ;
+
+	A( 5 , 0 ) = 2.0*(M_PI-alpha)/3.0 ;
+	A( 5 , 1 ) = M_PI*sin(alpha)/8.0 ;
+	A( 5 , 2 ) = 0 ;
+	A( 5 , 3 ) = 0 ;
+	A( 5 , 4 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
+	A( 5 , 5 ) = 2.0*(M_PI-alpha)/5.0 ;
 
 	if (size < 7)
 		return A ;
 
-	A( 6 , 0 ) =  0 ;
-	A( 6 , 1 ) =  0 ;
-	A( 6 , 2 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
-	A( 6 , 3 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
-	A( 6 , 4 ) =  0 ;
-	A( 6 , 5 ) =  0 ;
-	A( 6 , 6 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ;
-	A( 6 , 7 ) =  0 ;
-	A( 6 , 8 ) =  0 ;
-	A( 6 , 9 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
-	A( 7 , 0 ) =  M_PI*sin(alpha)/8.0 ;
-	A( 7 , 1 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
-	A( 7 , 2 ) =  0 ;
-	A( 7 , 3 ) =  0 ;
-	A( 7 , 4 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
-	A( 7 , 5 ) =  M_PI*sin(alpha)/1.6E+1 ;
-	A( 7 , 6 ) =  0 ;
-	A( 7 , 7 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
-	A( 7 , 8 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ;
-	A( 7 , 9 ) =  0 ;
-	A( 8 , 0 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/8.0 ;
-	A( 8 , 1 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ;
-	A( 8 , 2 ) =  0 ;
-	A( 8 , 3 ) =  0 ;
-	A( 8 , 4 ) =  M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/4.8E+1 ;
-	A( 8 , 5 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
-	A( 8 , 6 ) =  0 ;
-	A( 8 , 7 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ;
-	A( 8 , 8 ) =  -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/2.1E+2 ;
-	A( 8 , 9 ) =  0 ;
-	A( 9 , 0 ) =  0 ;
-	A( 9 , 1 ) =  0 ;
-	A( 9 , 2 ) =  2.0*(M_PI-alpha)/5.0 ;
-	A( 9 , 3 ) =  M_PI*sin(alpha)/1.6E+1 ;
-	A( 9 , 4 ) =  0 ;
-	A( 9 , 5 ) =  0 ;
-	A( 9 , 6 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
-	A( 9 , 7 ) =  0 ;
-	A( 9 , 8 ) =  0 ;
-	A( 9 , 9 ) =  2.0*(M_PI-alpha)/7.0 ;
+	A( 6 , 0 ) = 0 ;
+	A( 6 , 1 ) = 0 ;
+	A( 6 , 2 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
+	A( 6 , 3 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
+	A( 6 , 4 ) = 0 ;
+	A( 6 , 5 ) = 0 ;
+	A( 6 , 6 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ;
+
+	A( 7 , 0 ) = M_PI*sin(alpha)/8.0 ;
+	A( 7 , 1 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
+	A( 7 , 2 ) = 0 ;
+	A( 7 , 3 ) = 0 ;
+	A( 7 , 4 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
+	A( 7 , 5 ) = M_PI*sin(alpha)/1.6E+1 ;
+	A( 7 , 6 ) = 0 ;
+	A( 7 , 7 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
+
+	A( 8 , 0 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/8.0 ;
+	A( 8 , 1 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ;
+	A( 8 , 2 ) = 0 ;
+	A( 8 , 3 ) = 0 ;
+	A( 8 , 4 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/4.8E+1 ;
+	A( 8 , 5 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
+	A( 8 , 6 ) = 0 ;
+	A( 8 , 7 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ;
+	A( 8 , 8 ) = -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/2.1E+2 ;
+
+	A( 9 , 0 ) = 0 ;
+	A( 9 , 1 ) = 0 ;
+	A( 9 , 2 ) = 2.0*(M_PI-alpha)/5.0 ;
+	A( 9 , 3 ) = M_PI*sin(alpha)/1.6E+1 ;
+	A( 9 , 4 ) = 0 ;
+	A( 9 , 5 ) = 0 ;
+	A( 9 , 6 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
+	A( 9 , 7 ) = 0 ;
+	A( 9 , 8 ) = 0 ;
+	A( 9 , 9 ) = 2.0*(M_PI-alpha)/7.0 ;
 
 	if (size < 11)
 		return A ;
 
-	A( 10 , 0 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
-	A( 10 , 1 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
-	A( 10 , 2 ) =  0 ;
-	A( 10 , 3 ) =  0 ;
-	A( 10 , 4 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ;
-	A( 10 , 5 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
-	A( 10 , 6 ) =  0 ;
-	A( 10 , 7 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ;
-	A( 10 , 8 ) =  M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/3.84E+2 ;
-	A( 10 , 9 ) =  0 ;
-	A( 10 , 10 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ;
-	A( 10 , 11 ) =  0 ;
-	A( 10 , 12 ) =  0 ;
-	A( 10 , 13 ) =  -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/1.89E+3 ;
-	A( 10 , 14 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/6.3E+1 ;
-	A( 11 , 0 ) =  0 ;
-	A( 11 , 1 ) =  0 ;
-	A( 11 , 2 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
-	A( 11 , 3 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ;
-	A( 11 , 4 ) =  0 ;
-	A( 11 , 5 ) =  0 ;
-	A( 11 , 6 ) =  M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/3.84E+2 ;
-	A( 11 , 7 ) =  0 ;
-	A( 11 , 8 ) =  0 ;
-	A( 11 , 9 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ;
-	A( 11 , 10 ) =  0 ;
-	A( 11 , 11 ) =  -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/1.89E+3 ;
-	A( 11 , 12 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ;
-	A( 11 , 13 ) =  0 ;
-	A( 11 , 14 ) =  0 ;
-	A( 12 , 0 ) =  0 ;
-	A( 12 , 1 ) =  0 ;
-	A( 12 , 2 ) =  M_PI*sin(alpha)/1.6E+1 ;
-	A( 12 , 3 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
-	A( 12 , 4 ) =  0 ;
-	A( 12 , 5 ) =  0 ;
-	A( 12 , 6 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ;
-	A( 12 , 7 ) =  0 ;
-	A( 12 , 8 ) =  0 ;
-	A( 12 , 9 ) =  5.0*M_PI*sin(alpha)/1.28E+2 ;
-	A( 12 , 10 ) =  0 ;
-	A( 12 , 11 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ;
-	A( 12 , 12 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/6.3E+1 ;
-	A( 12 , 13 ) =  0 ;
-	A( 12 , 14 ) =  0 ;
-	A( 13 , 0 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ;
-	A( 13 , 1 ) =  M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/4.8E+1 ;
-	A( 13 , 2 ) =  0 ;
-	A( 13 , 3 ) =  0 ;
-	A( 13 , 4 ) =  -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/2.1E+2 ;
-	A( 13 , 5 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ;
-	A( 13 , 6 ) =  0 ;
-	A( 13 , 7 ) =  M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/3.84E+2 ;
-	A( 13 , 8 ) =  -M_PI*(5*pow(sin(alpha),7)-21*pow(sin(alpha),5)+35*pow(sin(alpha),3)-35*sin(alpha))/1.28E+2 ;
-	A( 13 , 9 ) =  0 ;
-	A( 13 , 10 ) =  -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/1.89E+3 ;
-	A( 13 , 11 ) =  0 ;
-	A( 13 , 12 ) =  0 ;
-	A( 13 , 13 ) =  -(3*sin(8*alpha)+168*sin(4*alpha)-128*pow(sin(2*alpha),3)+768*sin(2*alpha)+840*alpha-840*M_PI)/3.78E+3 ;
-	A( 13 , 14 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ;
-	A( 14 , 0 ) =  2.0*(M_PI-alpha)/5.0 ;
-	A( 14 , 1 ) =  M_PI*sin(alpha)/1.6E+1 ;
-	A( 14 , 2 ) =  0 ;
-	A( 14 , 3 ) =  0 ;
-	A( 14 , 4 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
-	A( 14 , 5 ) =  2.0*(M_PI-alpha)/7.0 ;
-	A( 14 , 6 ) =  0 ;
-	A( 14 , 7 ) =  5.0*M_PI*sin(alpha)/1.28E+2 ;
-	A( 14 , 8 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ;
-	A( 14 , 9 ) =  0 ;
-	A( 14 , 10 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/6.3E+1 ;
-	A( 14 , 11 ) =  0 ;
-	A( 14 , 12 ) =  0 ;
-	A( 14 , 13 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ;
-	A( 14 , 14 ) =  2.0*(M_PI-alpha)/9.0 ;
+	A( 10 , 0 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
+	A( 10 , 1 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
+	A( 10 , 2 ) = 0 ;
+	A( 10 , 3 ) = 0 ;
+	A( 10 , 4 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ;
+	A( 10 , 5 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
+	A( 10 , 6 ) = 0 ;
+	A( 10 , 7 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ;
+	A( 10 , 8 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/3.84E+2 ;
+	A( 10 , 9 ) = 0 ;
+	A( 10 , 10 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ;
+
+	A( 11 , 0 ) = 0 ;
+	A( 11 , 1 ) = 0 ;
+	A( 11 , 2 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
+	A( 11 , 3 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ;
+	A( 11 , 4 ) = 0 ;
+	A( 11 , 5 ) = 0 ;
+	A( 11 , 6 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/3.84E+2 ;
+	A( 11 , 7 ) = 0 ;
+	A( 11 , 8 ) = 0 ;
+	A( 11 , 9 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ;
+	A( 11 , 10 ) = 0 ;
+	A( 11 , 11 ) = -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/1.89E+3 ;
+
+	A( 12 , 0 ) = 0 ;
+	A( 12 , 1 ) = 0 ;
+	A( 12 , 2 ) = M_PI*sin(alpha)/1.6E+1 ;
+	A( 12 , 3 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
+	A( 12 , 4 ) = 0 ;
+	A( 12 , 5 ) = 0 ;
+	A( 12 , 6 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ;
+	A( 12 , 7 ) = 0 ;
+	A( 12 , 8 ) = 0 ;
+	A( 12 , 9 ) = 5.0*M_PI*sin(alpha)/1.28E+2 ;
+	A( 12 , 10 ) = 0 ;
+	A( 12 , 11 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ;
+	A( 12 , 12 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/6.3E+1 ;
+
+	A( 13 , 0 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ;
+	A( 13 , 1 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/4.8E+1 ;
+	A( 13 , 2 ) = 0 ;
+	A( 13 , 3 ) = 0 ;
+	A( 13 , 4 ) = -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/2.1E+2 ;
+	A( 13 , 5 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ;
+	A( 13 , 6 ) = 0 ;
+	A( 13 , 7 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/3.84E+2 ;
+	A( 13 , 8 ) = -M_PI*(5*pow(sin(alpha),7)-21*pow(sin(alpha),5)+35*pow(sin(alpha),3)-35*sin(alpha))/1.28E+2 ;
+	A( 13 , 9 ) = 0 ;
+	A( 13 , 10 ) = -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/1.89E+3 ;
+	A( 13 , 11 ) = 0 ;
+	A( 13 , 12 ) = 0 ;
+	A( 13 , 13 ) = -(3*sin(8*alpha)+168*sin(4*alpha)-128*pow(sin(2*alpha),3)+768*sin(2*alpha)+840*alpha-840*M_PI)/3.78E+3 ;
+
+	A( 14 , 0 ) = 2.0*(M_PI-alpha)/5.0 ;
+	A( 14 , 1 ) = M_PI*sin(alpha)/1.6E+1 ;
+	A( 14 , 2 ) = 0 ;
+	A( 14 , 3 ) = 0 ;
+	A( 14 , 4 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
+	A( 14 , 5 ) = 2.0*(M_PI-alpha)/7.0 ;
+	A( 14 , 6 ) = 0 ;
+	A( 14 , 7 ) = 5.0*M_PI*sin(alpha)/1.28E+2 ;
+	A( 14 , 8 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ;
+	A( 14 , 9 ) = 0 ;
+	A( 14 , 10 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/6.3E+1 ;
+	A( 14 , 11 ) = 0 ;
+	A( 14 , 12 ) = 0 ;
+	A( 14 , 13 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ;
+	A( 14 , 14 ) = 2.0*(M_PI-alpha)/9.0 ;
 
 	return A ;
 }
@@ -753,362 +728,421 @@ QuadricHF<REAL>::buildIntegralMatrix_B(const REAL& alpha, unsigned int size)
 {
 	Eigen::MatrixXd B(size,size) ;
 
-	B( 0 , 0 ) =  2*(M_PI-alpha) ;
-	B( 0 , 1 ) =  M_PI*(sin(2*alpha)+sin(alpha))/2.0 ;
-	B( 0 , 2 ) =  0 ;
-	B( 0 , 3 ) =  0 ;
-	B( 0 , 4 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.0 ;
-	B( 0 , 5 ) =  2.0*(M_PI-alpha)/3.0 ;
-	B( 1 , 0 ) =  M_PI*sin(alpha)/2.0 ;
-	B( 1 , 1 ) =  4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/3.0 ;
-	B( 1 , 2 ) =  0 ;
-	B( 1 , 3 ) =  0 ;
-	B( 1 , 4 ) =  3.0*M_PI*((sin(5*alpha)+3*sin(3*alpha)+6*sin(alpha))/1.2E+1+sin( 2*alpha)/3.0)/8.0 ;
-	B( 1 , 5 ) =  M_PI*sin(alpha)/8.0 ;
-	B( 2 , 0 ) =  0 ;
-	B( 2 , 1 ) =  0 ;
-	B( 2 , 2 ) =  2.0*(M_PI-alpha)/3.0 ;
-	B( 2 , 3 ) =  M_PI*(sin(2*alpha)+sin(alpha))/8.0 ;
-	B( 2 , 4 ) =  0 ;
-	B( 2 , 5 ) =  0 ;
-	B( 3 , 0 ) =  0 ;
-	B( 3 , 1 ) =  0 ;
-	B( 3 , 2 ) =  M_PI*sin(alpha)/8.0 ;
-	B( 3 , 3 ) =  4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/1.5E+1 ;
-	B( 3 , 4 ) =  0 ;
-	B( 3 , 5 ) =  0 ;
-	B( 4 , 0 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ;
-	B( 4 , 1 ) =  3.0*M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0 )/8.0 ;
-	B( 4 , 2 ) =  0 ;
-	B( 4 , 3 ) =  0 ;
-	B( 4 , 4 ) =  1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/1.5E+1 ;
-	B( 4 , 5 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
-	B( 5 , 0 ) =  2.0*(M_PI-alpha)/3.0 ;
-	B( 5 , 1 ) =  M_PI*(sin(2*alpha)+sin(alpha))/8.0 ;
-	B( 5 , 2 ) =  0 ;
-	B( 5 , 3 ) =  0 ;
-	B( 5 , 4 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
-	B( 5 , 5 ) =  2.0*(M_PI-alpha)/5.0 ;
-
-	if (size < 7)
-		return B ;
-
-	B( 6 , 0 ) =  0 ;
-	B( 6 , 1 ) =  0 ;
-	B( 6 , 2 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
-	B( 6 , 3 ) =  M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0)/1.6E+1 ;
-	B( 6 , 4 ) =  0 ;
-	B( 6 , 5 ) =  0 ;
-	B( 6 , 6 ) =  1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/1.05E+2 ;
-	B( 6 , 7 ) =  0 ;
-	B( 6 , 8 ) =  0 ;
-	B( 6 , 9 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
-	B( 7 , 0 ) =  M_PI*sin(alpha)/8.0 ;
-	B( 7 , 1 ) =  4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/1.5E+1 ;
-	B( 7 , 2 ) =  0 ;
-	B( 7 , 3 ) =  0 ;
-	B( 7 , 4 ) =  M_PI*((sin(5*alpha)+3*sin(3*alpha)+6*sin(alpha))/1.2E+1+sin(2*alpha )/3.0)/1.6E+1 ;
-	B( 7 , 5 ) =  M_PI*sin(alpha)/1.6E+1 ;
-	B( 7 , 6 ) =  0 ;
-	B( 7 , 7 ) =  4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/3.5E+1 ;
-	B( 7 , 8 ) =  1.6E+1*((3*sin(3*alpha)+6*sin(alpha))/3.2E+1-(sin(7*alpha)+2*sin(5 *alpha)+6*sin(3*alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.05E+2 ;
-	B( 7 , 9 ) =  0 ;
-	B( 8 , 0 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/8.0 ;
-	B( 8 , 1 ) =  1.6E+1*(5.0*sin(alpha)/3.2E+1-(sin(5*alpha)+6*sin(3*alpha)+2*sin(alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.5E+1 ;
-	B( 8 , 2 ) =  0 ;
-	B( 8 , 3 ) =  0 ;
-	B( 8 , 4 ) =  5.0*M_PI*((3*sin(7*alpha)+15*sin(5*alpha)+55*sin(3*alpha)+75*sin (alpha))/2.4E+2+sin(2*alpha)/5.0)/1.6E+1 ;
-	B( 8 , 5 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
-	B( 8 , 6 ) =  0 ;
-	B( 8 , 7 ) =  1.6E+1*(5.0*sin(alpha)/3.2E+1-(sin(5*alpha)+6*sin(3*alpha)+2*sin(alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.05E+2 ;
-	B( 8 , 8 ) =  3.2E+1*((29*sin(3*alpha)+45*sin(alpha))/3.84E+2-(2*sin(9*alpha)+9*sin(7*alpha)+27*sin(5*alpha)+54*sin(3*alpha)+(12*alpha-12*M_PI)*cos(3*alpha)+18*sin(alpha)+(108*alpha-108*M_PI)*cos(alpha))/3.84E+2)/3.5E+1 ;
-	B( 8 , 9 ) =  0 ;
-	B( 9 , 0 ) =  0 ;
-	B( 9 , 1 ) =  0 ;
-	B( 9 , 2 ) =  2.0*(M_PI-alpha)/5.0 ;
-	B( 9 , 3 ) =  M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
-	B( 9 , 4 ) =  0 ;
-	B( 9 , 5 ) =  0 ;
-	B( 9 , 6 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
-	B( 9 , 7 ) =  0 ;
-	B( 9 , 8 ) =  0 ;
-	B( 9 , 9 ) =  2.0*(M_PI-alpha)/7.0 ;
-
-	if (size < 11)
-		return B ;
-
-	B( 10 , 0 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
-	B( 10 , 1 ) =  M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0)/1.6E+1 ;
-	B( 10 , 2 ) =  0 ;
-	B( 10 , 3 ) =  0 ;
-	B( 10 , 4 ) =  1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/1.05E+2 ;
-	B( 10 , 5 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
-	B( 10 , 6 ) =  0 ;
-	B( 10 , 7 ) =  3.0*M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0 )/1.28E+2 ;
-	B( 10 , 8 ) =  5.0*M_PI*((3*sin(8*alpha)+10*sin(6*alpha)+30*sin(4*alpha)+90*sin (2*alpha))/2.4E+2+(7*sin(3*alpha)+15*sin(alpha))/6.0E+1)/1.28E+2 ;
-	B( 10 , 9 ) =  0 ;
-	B( 10 , 10 ) =  1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/3.15E+2 ;
-	B( 10 , 11 ) =  0 ;
-	B( 10 , 12 ) =  0 ;
-	B( 10 , 13 ) =  3.2E+1*((7*sin(4*alpha)+30*sin(2*alpha))/1.92E+2-(sin(10*alpha)+3* sin(8*alpha)+9*sin(6*alpha)+24*sin(4*alpha)+18*sin(2*alpha)+(24 *alpha-24*M_PI)*cos(2*alpha)+36*alpha-36*M_PI)/1.92E+2)/3.15E+2 ;
-	B( 10 , 14 ) =  -(sin(2*alpha)+2*alpha-2*M_PI)/6.3E+1 ;
-	B( 11 , 0 ) =  0 ;
-	B( 11 , 1 ) =  0 ;
-	B( 11 , 2 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
-	B( 11 , 3 ) =  1.6E+1*(5.0*sin(alpha)/3.2E+1-(sin(5*alpha)+6*sin(3*alpha)+2*sin(alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.05E+2 ;
-	B( 11 , 4 ) =  0 ;
-	B( 11 , 5 ) =  0 ;
-	B( 11 , 6 ) =  5.0*M_PI*((3*sin(7*alpha)+15*sin(5*alpha)+55*sin(3*alpha)+75*sin (alpha))/2.4E+2+sin(2*alpha)/5.0)/1.28E+2 ;
-	B( 11 , 7 ) =  0 ;
-	B( 11 , 8 ) =  0 ;
-	B( 11 , 9 ) =  -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ;
-	B( 11 , 10 ) =  0 ;
-	B( 11 , 11 ) =  3.2E+1*((29*sin(3*alpha)+45*sin(alpha))/3.84E+2-(2*sin(9*alpha)+9*sin(7*alpha)+27*sin(5*alpha)+54*sin(3*alpha)+(12*alpha-12*M_PI)* cos(3*alpha)+18*sin(alpha)+(108*alpha-108*M_PI)*cos(alpha))/3.84E+2)/3.15E+2 ;
-	B( 11 , 12 ) =  1.6E+1*(5.0*sin(alpha)/3.2E+1-(sin(5*alpha)+6*sin(3*alpha)+2*sin(alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/3.15E+2 ;
-	B( 11 , 13 ) =  0 ;
-	B( 11 , 14 ) =  0 ;
-	B( 12 , 0 ) =  0 ;
-	B( 12 , 1 ) =  0 ;
-	B( 12 , 2 ) =  M_PI*sin(alpha)/1.6E+1 ;
-	B( 12 , 3 ) =  4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/3.5E+1 ;
-	B( 12 , 4 ) =  0 ;
-	B( 12 , 5 ) =  0 ;
-	B( 12 , 6 ) =  3.0*M_PI*((sin(5*alpha)+3*sin(3*alpha)+6*sin(alpha))/1.2E+1+sin( 2*alpha)/3.0)/1.28E+2 ;
-	B( 12 , 7 ) =  0 ;
-	B( 12 , 8 ) =  0 ;
-	B( 12 , 9 ) =  5.0*M_PI*sin(alpha)/1.28E+2 ;
-	B( 12 , 10 ) =  0 ;
-	B( 12 , 11 ) =  1.6E+1*((3*sin(3*alpha)+6*sin(alpha))/3.2E+1-(sin(7*alpha)+2*sin(5 *alpha)+6*sin(3*alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/3.15E+2 ;
-	B( 12 , 12 ) =  4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/6.3E+1 ;
-	B( 12 , 13 ) =  0 ;
-	B( 12 , 14 ) =  0 ;
-	B( 13 , 0 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ;
-	B( 13 , 1 ) =  5.0*M_PI*((3*sin(6*alpha)+20*sin(4*alpha)+95*sin(2*alpha))/2.4E+2 + sin(alpha)/5.0)/1.6E+1 ;
-	B( 13 , 2 ) =  0 ;
-	B( 13 , 3 ) =  0 ;
-	B( 13 , 4 ) =  3.2E+1*(1.1E+1*sin(2*alpha)/9.6E+1-(sin(8*alpha)+6*sin(6*alpha)+21 *sin(4*alpha)+24*sin(2*alpha)+(24*alpha-24*M_PI)*cos(2*alpha)+36 *alpha-36*M_PI)/1.92E+2)/3.5E+1 ;
-	B( 13 , 5 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ;
-	B( 13 , 6 ) =  0 ;
-	B( 13 , 7 ) =  5.0*M_PI*((3*sin(6*alpha)+20*sin(4*alpha)+95*sin(2*alpha))/2.4E+2+sin(alpha)/5.0)/1.28E+2 ;
-	B( 13 , 8 ) =  3.5E+1*M_PI*((5*sin(10*alpha)+28*sin(8*alpha)+91*sin(6*alpha)+280*sin(4*alpha)+630*sin(2*alpha))/2.24E+3+(13*sin(3*alpha)+21*sin(alpha))/1.4E+2)/1.28E+2 ;
-	B( 13 , 9 ) =  0 ;
-	B( 13 , 10 ) =  3.2E+1*(1.1E+1*sin(2*alpha)/9.6E+1-(sin(8*alpha)+6*sin(6*alpha)+21 *sin(4*alpha)+24*sin(2*alpha)+(24*alpha-24*M_PI)*cos(2*alpha)+36 *alpha-36*M_PI)/1.92E+2)/3.15E+2 ;
-	B( 13 , 11 ) =  0 ;
-	B( 13 , 12 ) =  0 ;
-	B( 13 , 13 ) =  2.56E+2*((103*sin(4*alpha)+352*sin(2*alpha))/3.072E+3-(3*sin(12*alpha)+16*sin(10*alpha)+52*sin(8*alpha)+144*sin(6*alpha)+324*sin(4*alpha)+(24*alpha-24*M_PI)*cos(4*alpha)+288*sin(2*alpha)+(384*alpha-384*M_PI)*cos(2*alpha)+432*alpha-432*M_PI)/3.072E+3)/3.15E+2 ;
-	B( 13 , 14 ) =  -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ;
-	B( 14 , 0 ) =  2.0*(M_PI-alpha)/5.0 ;
-	B( 14 , 1 ) =  M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
-	B( 14 , 2 ) =  0 ;
-	B( 14 , 3 ) =  0 ;
-	B( 14 , 4 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
-	B( 14 , 5 ) =  2.0*(M_PI-alpha)/7.0 ;
-	B( 14 , 6 ) =  0 ;
-	B( 14 , 7 ) =  5.0*M_PI*(sin(2*alpha)+sin(alpha))/1.28E+2 ;
-	B( 14 , 8 ) =  3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/1.28E+2 ;
-	B( 14 , 9 ) =  0 ;
-	B( 14 , 10 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/6.3E+1 ;
-	B( 14 , 11 ) =  0 ;
-	B( 14 , 12 ) =  0 ;
-	B( 14 , 13 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/3.15E+2 ;
-	B( 14 , 14 ) =  2.0*(M_PI-alpha)/9.0 ;
+	 B( 0 , 0 ) = 2*(M_PI-alpha) ;
+	 B( 0 , 1 ) = M_PI*(sin(2*alpha)+sin(alpha))/2.0 ;
+	 B( 0 , 2 ) = 0 ;
+	 B( 0 , 3 ) = 0 ;
+	 B( 0 , 4 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.0 ;
+	 B( 0 , 5 ) = 2.0*(M_PI-alpha)/3.0 ;
+
+	 B( 1 , 0 ) = M_PI*sin(alpha)/2.0 ;
+	 B( 1 , 1 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/3.0 ;
+	 B( 1 , 2 ) = 0 ;
+	 B( 1 , 3 ) = 0 ;
+	 B( 1 , 4 ) = 3.0*M_PI*((sin(5*alpha)+3*sin(3*alpha)+6*sin(alpha))/1.2E+1+sin( 2*alpha)/3.0)/8.0 ;
+	 B( 1 , 5 ) = M_PI*sin(alpha)/8.0 ;
+
+	 B( 2 , 0 ) = 0 ;
+	 B( 2 , 1 ) = 0 ;
+	 B( 2 , 2 ) = 2.0*(M_PI-alpha)/3.0 ;
+	 B( 2 , 3 ) = M_PI*(sin(2*alpha)+sin(alpha))/8.0 ;
+	 B( 2 , 4 ) = 0 ;
+	 B( 2 , 5 ) = 0 ;
+
+	 B( 3 , 0 ) = 0 ;
+	 B( 3 , 1 ) = 0 ;
+	 B( 3 , 2 ) = M_PI*sin(alpha)/8.0 ;
+	 B( 3 , 3 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/1.5E+1 ;
+	 B( 3 , 4 ) = 0 ;
+	 B( 3 , 5 ) = 0 ;
+
+	 B( 4 , 0 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ;
+	 B( 4 , 1 ) = 3.0*M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0 )/8.0 ;
+	 B( 4 , 2 ) = 0 ;
+	 B( 4 , 3 ) = 0 ;
+	 B( 4 , 4 ) = 1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/1.5E+1 ;
+	 B( 4 , 5 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
+
+	 B( 5 , 0 ) = 2.0*(M_PI-alpha)/3.0 ;
+	 B( 5 , 1 ) = M_PI*(sin(2*alpha)+sin(alpha))/8.0 ;
+	 B( 5 , 2 ) = 0 ;
+	 B( 5 , 3 ) = 0 ;
+	 B( 5 , 4 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
+	 B( 5 , 5 ) = 2.0*(M_PI-alpha)/5.0 ;
+
+	 if (size < 7)
+		 return B ;
+
+	 B( 0 , 6 ) = 0 ;
+	 B( 0 , 7 ) = M_PI*(sin(2*alpha)+sin(alpha))/8.0 ;
+	 B( 0 , 8 ) = 3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/8.0 ;
+	 B( 0 , 9 ) = 0 ;
+
+	 B( 1 , 6 ) = 0 ;
+	 B( 1 , 7 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/1.5E+1 ;
+	 B( 1 , 8 ) = 1.6E+1*((3*sin(3*alpha)+6*sin(alpha))/3.2E+1-(sin(7*alpha)+2*sin(5 *alpha)+6*sin(3*alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.5E+1 ;
+	 B( 1 , 9 ) = 0 ;
+
+	 B( 2 , 6 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
+	 B( 2 , 7 ) = 0 ;
+	 B( 2 , 8 ) = 0 ;
+	 B( 2 , 9 ) = 2.0*(M_PI-alpha)/5.0 ;
+
+	 B( 3 , 6 ) = M_PI*((sin(5*alpha)+3*sin(3*alpha)+6*sin(alpha))/1.2E+1+sin(2*alpha )/3.0)/1.6E+1 ;
+	 B( 3 , 7 ) = 0 ;
+	 B( 3 , 8 ) = 0 ;
+	 B( 3 , 9 ) = M_PI*sin(alpha)/1.6E+1 ;
+
+	 B( 4 , 6 ) = 0 ;
+	 B( 4 , 7 ) = M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0)/1.6E+1 ;
+	 B( 4 , 8 ) = 5.0*M_PI*((3*sin(8*alpha)+10*sin(6*alpha)+30*sin(4*alpha)+90*sin (2*alpha))/2.4E+2+(7*sin(3*alpha)+15*sin(alpha))/6.0E+1)/1.6E+1 ;
+	 B( 4 , 9 ) = 0 ;
+
+	 B( 5 , 6 ) = 0 ;
+	 B( 5 , 7 ) = M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
+	 B( 5 , 8 ) = M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin (alpha))/3.0)/1.6E+1 ;
+	 B( 5 , 9 ) = 0 ;
+
+	 B( 6 , 0 ) = 0 ;
+	 B( 6 , 1 ) = 0 ;
+	 B( 6 , 2 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
+	 B( 6 , 3 ) = M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0)/1.6E+1 ;
+	 B( 6 , 4 ) = 0 ;
+	 B( 6 , 5 ) = 0 ;
+	 B( 6 , 6 ) = 1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/1.05E+2 ;
+	 B( 6 , 7 ) = 0 ;
+	 B( 6 , 8 ) = 0 ;
+	 B( 6 , 9 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
+
+	 B( 7 , 0 ) = M_PI*sin(alpha)/8.0 ;
+	 B( 7 , 1 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/1.5E+1 ;
+	 B( 7 , 2 ) = 0 ;
+	 B( 7 , 3 ) = 0 ;
+	 B( 7 , 4 ) = M_PI*((sin(5*alpha)+3*sin(3*alpha)+6*sin(alpha))/1.2E+1+sin(2*alpha )/3.0)/1.6E+1 ;
+	 B( 7 , 5 ) = M_PI*sin(alpha)/1.6E+1 ;
+	 B( 7 , 6 ) = 0 ;
+	 B( 7 , 7 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/3.5E+1 ;
+	 B( 7 , 8 ) = 1.6E+1*((3*sin(3*alpha)+6*sin(alpha))/3.2E+1-(sin(7*alpha)+2*sin(5 *alpha)+6*sin(3*alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.05E+2 ;
+	 B( 7 , 9 ) = 0 ;
+
+	 B( 8 , 0 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/8.0 ;
+	 B( 8 , 1 ) = 1.6E+1*(5.0*sin(alpha)/3.2E+1-(sin(5*alpha)+6*sin(3*alpha)+2*sin(alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.5E+1 ;
+	 B( 8 , 2 ) = 0 ;
+	 B( 8 , 3 ) = 0 ;
+	 B( 8 , 4 ) = 5.0*M_PI*((3*sin(7*alpha)+15*sin(5*alpha)+55*sin(3*alpha)+75*sin (alpha))/2.4E+2+sin(2*alpha)/5.0)/1.6E+1 ;
+	 B( 8 , 5 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
+	 B( 8 , 6 ) = 0 ;
+	 B( 8 , 7 ) = 1.6E+1*(5.0*sin(alpha)/3.2E+1-(sin(5*alpha)+6*sin(3*alpha)+2*sin(alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.05E+2 ;
+	 B( 8 , 8 ) = 3.2E+1*((29*sin(3*alpha)+45*sin(alpha))/3.84E+2-(2*sin(9*alpha)+9*sin(7*alpha)+27*sin(5*alpha)+54*sin(3*alpha)+(12*alpha-12*M_PI)*cos(3*alpha)+18*sin(alpha)+(108*alpha-108*M_PI)*cos(alpha))/3.84E+2)/3.5E+1 ;
+	 B( 8 , 9 ) = 0 ;
+
+	 B( 9 , 0 ) = 0 ;
+	 B( 9 , 1 ) = 0 ;
+	 B( 9 , 2 ) = 2.0*(M_PI-alpha)/5.0 ;
+	 B( 9 , 3 ) = M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
+	 B( 9 , 4 ) = 0 ;
+	 B( 9 , 5 ) = 0 ;
+	 B( 9 , 6 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
+	 B( 9 , 7 ) = 0 ;
+	 B( 9 , 8 ) = 0 ;
+	 B( 9 , 9 ) = 2.0*(M_PI-alpha)/7.0 ;
+
+	 if (size < 10)
+		 return B ;
+
+	 B( 0 , 10 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
+	 B( 0 , 11 ) = 0 ;
+	 B( 0 , 12 ) = 0 ;
+	 B( 0 , 13 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.5E+1 ;
+	 B( 0 , 14 ) = 2.0*(M_PI-alpha)/5.0 ;
+
+	 B( 1 , 10 ) = M_PI*((sin(5*alpha)+3*sin(3*alpha)+6*sin(alpha))/1.2E+1+sin(2*alpha )/3.0)/1.6E+1 ;
+	 B( 1 , 11 ) = 0 ;
+	 B( 1 , 12 ) = 0 ;
+	 B( 1 , 13 ) = 5.0*M_PI*((3*sin(9*alpha)+5*sin(7*alpha)+20*sin(5*alpha)+60*sin(3*alpha)+90*sin(alpha))/2.4E+2+(sin(4*alpha)+10*sin(2*alpha))/3.0E+1)/1.6E+1 ;
+	 B( 1 , 14 ) = M_PI*sin(alpha)/1.6E+1 ;
+
+	 B( 2 , 10 ) = 0 ;
+	 B( 2 , 11 ) = M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin (alpha))/3.0)/1.6E+1 ;
+	 B( 2 , 12 ) = M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
+	 B( 2 , 13 ) = 0 ;
+	 B( 2 , 14 ) = 0 ;
+
+	 B( 3 , 10 ) = 0 ;
+	 B( 3 , 11 ) = 1.6E+1*((3*sin(3*alpha)+6*sin(alpha))/3.2E+1-(sin(7*alpha)+2*sin(5 *alpha)+6*sin(3*alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.05E+2 ;
+	 B( 3 , 12 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/3.5E+1 ;
+	 B( 3 , 13 ) = 0 ;
+	 B( 3 , 14 ) = 0 ;
+
+	 B( 4 , 10 ) = 1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/1.05E+2 ;
+	 B( 4 , 11 ) = 0 ;
+	 B( 4 , 12 ) = 0 ;
+	 B( 4 , 13 ) = 3.2E+1*((7*sin(4*alpha)+30*sin(2*alpha))/1.92E+2-(sin(10*alpha)+3* sin(8*alpha)+9*sin(6*alpha)+24*sin(4*alpha)+18*sin(2*alpha)+(24 *alpha-24*M_PI)*cos(2*alpha)+36*alpha-36*M_PI)/1.92E+2)/3.5E+1 ;
+	 B( 4 , 14 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
+
+	 B( 5 , 10 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
+	 B( 5 , 11 ) = 0 ;
+	 B( 5 , 12 ) = 0 ;
+	 B( 5 , 13 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.05E+2 ;
+	 B( 5 , 14 ) = 2.0*(M_PI-alpha)/7.0 ;
+
+	 B( 6 , 10 ) = 0 ;
+	 B( 6 , 11 ) = 5.0*M_PI*((3*sin(8*alpha)+10*sin(6*alpha)+30*sin(4*alpha)+90*sin (2*alpha))/2.4E+2+(7*sin(3*alpha)+15*sin(alpha))/6.0E+1)/1.28E+2 ;
+	 B( 6 , 12 ) = 3.0*M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0 )/1.28E+2 ;
+	 B( 6 , 13 ) = 0 ;
+	 B( 6 , 14 ) = 0 ;
+
+	 B( 7 , 10 ) = 3.0*M_PI*((sin(5*alpha)+3*sin(3*alpha)+6*sin(alpha))/1.2E+1+sin( 2*alpha)/3.0)/1.28E+2 ;
+	 B( 7 , 11 ) = 0 ;
+	 B( 7 , 12 ) = 0 ;
+	 B( 7 , 13 ) = 5.0*M_PI*((3*sin(9*alpha)+5*sin(7*alpha)+20*sin(5*alpha)+60*sin(3*alpha)+90*sin(alpha))/2.4E+2+(sin(4*alpha)+10*sin(2*alpha))/3.0E+1)/1.28E+2 ;
+	 B( 7 , 14 ) = 5.0*M_PI*sin(alpha)/1.28E+2 ;
+
+	 B( 8 , 10 ) = 5.0*M_PI*((3*sin(7*alpha)+15*sin(5*alpha)+55*sin(3*alpha)+75*sin (alpha))/2.4E+2+sin(2*alpha)/5.0)/1.28E+2 ;
+	 B( 8 , 11 ) = 0 ;
+	 B( 8 , 12 ) = 0 ;
+	 B( 8 , 13 ) = 3.5E+1*M_PI*((5*sin(11*alpha)+21*sin(9*alpha)+63*sin(7*alpha)+175*sin(5*alpha)+490*sin(3*alpha)+490*sin(alpha))/2.24E+3+(3*sin(4*alpha)+14*sin(2*alpha))/7.0E+1)/1.28E+2 ;
+	 B( 8 , 14 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ;
+
+	 B( 9 , 10 ) = 0 ;
+	 B( 9 , 11 ) = 3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/1.28E+2 ;
+	 B( 9 , 12 ) = 5.0*M_PI*(sin(2*alpha)+sin(alpha))/1.28E+2 ;
+	 B( 9 , 13 ) = 0 ;
+	 B( 9 , 14 ) = 0 ;
+
+	 B( 10 , 0 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;
+	 B( 10 , 1 ) = M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0)/1.6E+1 ;
+	 B( 10 , 2 ) = 0 ;
+	 B( 10 , 3 ) = 0 ;
+	 B( 10 , 4 ) = 1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/1.05E+2 ;
+	 B( 10 , 5 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ;
+	 B( 10 , 6 ) = 0 ;
+	 B( 10 , 7 ) = 3.0*M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0 )/1.28E+2 ;
+	 B( 10 , 8 ) = 5.0*M_PI*((3*sin(8*alpha)+10*sin(6*alpha)+30*sin(4*alpha)+90*sin (2*alpha))/2.4E+2+(7*sin(3*alpha)+15*sin(alpha))/6.0E+1)/1.28E+2 ;
+	 B( 10 , 9 ) = 0 ;
+	 B( 10 , 10 ) = 1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/3.15E+2 ;
+	 B( 10 , 11 ) = 0 ;
+	 B( 10 , 12 ) = 0 ;
+	 B( 10 , 13 ) = 3.2E+1*((7*sin(4*alpha)+30*sin(2*alpha))/1.92E+2-(sin(10*alpha)+3* sin(8*alpha)+9*sin(6*alpha)+24*sin(4*alpha)+18*sin(2*alpha)+(24 *alpha-24*M_PI)*cos(2*alpha)+36*alpha-36*M_PI)/1.92E+2)/3.15E+2 ;
+	 B( 10 , 14 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/6.3E+1 ;
+
+	 B( 11 , 0 ) = 0 ;
+	 B( 11 , 1 ) = 0 ;
+	 B( 11 , 2 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ;
+	 B( 11 , 3 ) = 1.6E+1*(5.0*sin(alpha)/3.2E+1-(sin(5*alpha)+6*sin(3*alpha)+2*sin(alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.05E+2 ;
+	 B( 11 , 4 ) = 0 ;
+	 B( 11 , 5 ) = 0 ;
+	 B( 11 , 6 ) = 5.0*M_PI*((3*sin(7*alpha)+15*sin(5*alpha)+55*sin(3*alpha)+75*sin (alpha))/2.4E+2+sin(2*alpha)/5.0)/1.28E+2 ;
+	 B( 11 , 7 ) = 0 ;
+	 B( 11 , 8 ) = 0 ;
+	 B( 11 , 9 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ;
+	 B( 11 , 10 ) = 0 ;
+	 B( 11 , 11 ) = 3.2E+1*((29*sin(3*alpha)+45*sin(alpha))/3.84E+2-(2*sin(9*alpha)+9*sin(7*alpha)+27*sin(5*alpha)+54*sin(3*alpha)+(12*alpha-12*M_PI)*cos(3*alpha)+18*sin(alpha)+(108*alpha-108*M_PI)*cos(alpha))/3.84E+2)/3.15E+2 ;
+	 B( 11 , 12 ) = 1.6E+1*(5.0*sin(alpha)/3.2E+1-(sin(5*alpha)+6*sin(3*alpha)+2*sin(alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/3.15E+2 ;
+	 B( 11 , 13 ) = 0 ;
+	 B( 11 , 14 ) = 0 ;
+
+	 B( 12 , 0 ) = 0 ;
+	 B( 12 , 1 ) = 0 ;
+	 B( 12 , 2 ) = M_PI*sin(alpha)/1.6E+1 ;
+	 B( 12 , 3 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/3.5E+1 ;
+	 B( 12 , 4 ) = 0 ;
+	 B( 12 , 5 ) = 0 ;
+	 B( 12 , 6 ) = 3.0*M_PI*((sin(5*alpha)+3*sin(3*alpha)+6*sin(alpha))/1.2E+1+sin( 2*alpha)/3.0)/1.28E+2 ;
+	 B( 12 , 7 ) = 0 ;
+	 B( 12 , 8 ) = 0 ;
+	 B( 12 , 9 ) = 5.0*M_PI*sin(alpha)/1.28E+2 ;
+	 B( 12 , 10 ) = 0 ;
+	 B( 12 , 11 ) = 1.6E+1*((3*sin(3*alpha)+6*sin(alpha))/3.2E+1-(sin(7*alpha)+2*sin(5 *alpha)+6*sin(3*alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/3.15E+2 ;
+	 B( 12 , 12 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/6.3E+1 ;
+	 B( 12 , 13 ) = 0 ;
+	 B( 12 , 14 ) = 0 ;
+
+	 B( 13 , 0 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ;
+	 B( 13 , 1 ) = 5.0*M_PI*((3*sin(6*alpha)+20*sin(4*alpha)+95*sin(2*alpha))/2.4E+2+sin(alpha)/5.0)/1.6E+1 ;
+	 B( 13 , 2 ) = 0 ;
+	 B( 13 , 3 ) = 0 ;
+	 B( 13 , 4 ) = 3.2E+1*(1.1E+1*sin(2*alpha)/9.6E+1-(sin(8*alpha)+6*sin(6*alpha)+21 *sin(4*alpha)+24*sin(2*alpha)+(24*alpha-24*M_PI)*cos(2*alpha)+36 *alpha-36*M_PI)/1.92E+2)/3.5E+1 ;
+	 B( 13 , 5 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ;
+	 B( 13 , 6 ) = 0 ;
+	 B( 13 , 7 ) = 5.0*M_PI*((3*sin(6*alpha)+20*sin(4*alpha)+95*sin(2*alpha))/2.4E+2+sin(alpha)/5.0)/1.28E+2 ;
+	 B( 13 , 8 ) = 3.5E+1*M_PI*((5*sin(10*alpha)+28*sin(8*alpha)+91*sin(6*alpha)+280*sin(4*alpha)+630*sin(2*alpha))/2.24E+3+(13*sin(3*alpha)+21*sin(alpha))/1.4E+2)/1.28E+2 ;
+	 B( 13 , 9 ) = 0 ;
+	 B( 13 , 10 ) = 3.2E+1*(1.1E+1*sin(2*alpha)/9.6E+1-(sin(8*alpha)+6*sin(6*alpha)+21 *sin(4*alpha)+24*sin(2*alpha)+(24*alpha-24*M_PI)*cos(2*alpha)+36 *alpha-36*M_PI)/1.92E+2)/3.15E+2 ;
+	 B( 13 , 11 ) = 0 ;
+	 B( 13 , 12 ) = 0 ;
+	 B( 13 , 13 ) = 2.56E+2*((103*sin(4*alpha)+352*sin(2*alpha))/3.072E+3-(3*sin(12*alpha)+16*sin(10*alpha)+52*sin(8*alpha)+144*sin(6*alpha)+324*sin(4*alpha)+(24*alpha-24*M_PI)*cos(4*alpha)+288*sin(2*alpha)+(384*alpha-384*M_PI)*cos(2*alpha)+432*alpha-432*M_PI)/3.072E+3)/3.15E+2 ;
+	 B( 13 , 14 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ;
+
+	 B( 14 , 0 ) = 2.0*(M_PI-alpha)/5.0 ;
+	 B( 14 , 1 ) = M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
+	 B( 14 , 2 ) = 0 ;
+	 B( 14 , 3 ) = 0 ;
+	 B( 14 , 4 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
+	 B( 14 , 5 ) = 2.0*(M_PI-alpha)/7.0 ;
+	 B( 14 , 6 ) = 0 ;
+	 B( 14 , 7 ) = 5.0*M_PI*(sin(2*alpha)+sin(alpha))/1.28E+2 ;
+	 B( 14 , 8 ) = 3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/1.28E+2 ;
+	 B( 14 , 9 ) = 0 ;
+	 B( 14 , 10 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/6.3E+1 ;
+	 B( 14 , 11 ) = 0 ;
+	 B( 14 , 12 ) = 0 ;
+	 B( 14 , 13 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/3.15E+2 ;
+	 B( 14 , 14 ) = 2.0*(M_PI-alpha)/9.0 ;
 
 	return B ;
 }
 
 template <typename REAL>
 Eigen::MatrixXd
-QuadricHF<REAL>::buildIntegralMatrix_C(const REAL& alpha, unsigned int size)
+QuadricHF<REAL>::buildLowerLeftIntegralMatrix_C(const REAL& alpha, unsigned int size)
 {
 	Eigen::MatrixXd C(size,size) ;
 
-	C( 0 , 0 ) =  2*(M_PI-alpha) ;
-	C( 0 , 1 ) =  M_PI*(sin(2*alpha)+sin(alpha))/2.0 ;
-	C( 0 , 2 ) =  0 ;
-	C( 0 , 3 ) =  0 ;
-	C( 0 , 4 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.0 ;
-	C( 0 , 5 ) =  2.0*(M_PI-alpha)/3.0 ;
-	C( 1 , 0 ) =  M_PI*(sin(2*alpha)+sin(alpha))/2.0 ;
-	C( 1 , 1 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.0 ;
-	C( 1 , 2 ) =  0 ;
-	C( 1 , 3 ) =  0 ;
-	C( 1 , 4 ) =  3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/8.0 ;
-	C( 1 , 5 ) =  M_PI*(sin(2*alpha)+sin(alpha))/8.0 ;
-	C( 2 , 0 ) =  0 ;
-	C( 2 , 1 ) =  0 ;
-	C( 2 , 2 ) =  2.0*(M_PI-alpha)/3.0 ;
-	C( 2 , 3 ) =  M_PI*(sin(2*alpha)+sin(alpha))/8.0 ;
-	C( 2 , 4 ) =  0 ;
-	C( 2 , 5 ) =  0 ;
-	C( 3 , 0 ) =  0 ;
-	C( 3 , 1 ) =  0 ;
-	C( 3 , 2 ) =  M_PI*(sin(2*alpha)+sin(alpha))/8.0 ;
-	C( 3 , 3 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
-	C( 3 , 4 ) =  0 ;
-	C( 3 , 5 ) =  0 ;
-	C( 4 , 0 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.0 ;
-	C( 4 , 1 ) =  3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/8.0 ;
-	C( 4 , 2 ) =  0 ;
-	C( 4 , 3 ) =  0 ;
-	C( 4 , 4 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.5E+1 ;
-	C( 4 , 5 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
-	C( 5 , 0 ) =  2.0*(M_PI-alpha)/3.0 ;
-	C( 5 , 1 ) =  M_PI*(sin(2*alpha)+sin(alpha))/8.0 ;
-	C( 5 , 2 ) =  0 ;
-	C( 5 , 3 ) =  0 ;
-	C( 5 , 4 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
-	C( 5 , 5 ) =  2.0*(M_PI-alpha)/5.0 ;
+	C( 0 , 0 ) = 2*(M_PI-alpha) ;
+
+	C( 1 , 0 ) = M_PI*(sin(2*alpha)+sin(alpha))/2.0 ;
+	C( 1 , 1 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.0 ;
+
+	C( 2 , 0 ) = 0 ;
+	C( 2 , 1 ) = 0 ;
+	C( 2 , 2 ) = 2.0*(M_PI-alpha)/3.0 ;
+
+	C( 3 , 0 ) = 0 ;
+	C( 3 , 1 ) = 0 ;
+	C( 3 , 2 ) = M_PI*(sin(2*alpha)+sin(alpha))/8.0 ;
+	C( 3 , 3 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
+
+	C( 4 , 0 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.0 ;
+	C( 4 , 1 ) = 3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/8.0 ;
+	C( 4 , 2 ) = 0 ;
+	C( 4 , 3 ) = 0 ;
+	C( 4 , 4 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.5E+1 ;
+
+	C( 5 , 0 ) = 2.0*(M_PI-alpha)/3.0 ;
+	C( 5 , 1 ) = M_PI*(sin(2*alpha)+sin(alpha))/8.0 ;
+	C( 5 , 2 ) = 0 ;
+	C( 5 , 3 ) = 0 ;
+	C( 5 , 4 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
+	C( 5 , 5 ) = 2.0*(M_PI-alpha)/5.0 ;
 
 	if (size < 7)
 		return C ;
 
-	C( 6 , 0 ) =  0 ;
-	C( 6 , 1 ) =  0 ;
-	C( 6 , 2 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
-	C( 6 , 3 ) =  M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin (alpha))/3.0)/1.6E+1 ;
-	C( 6 , 4 ) =  0 ;
-	C( 6 , 5 ) =  0 ;
-	C( 6 , 6 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.05E+2 ;
-	C( 6 , 7 ) =  0 ;
-	C( 6 , 8 ) =  0 ;
-	C( 6 , 9 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
-	C( 7 , 0 ) =  M_PI*(sin(2*alpha)+sin(alpha))/8.0 ;
-	C( 7 , 1 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
-	C( 7 , 2 ) =  0 ;
-	C( 7 , 3 ) =  0 ;
-	C( 7 , 4 ) =  M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin (alpha))/3.0)/1.6E+1 ;
-	C( 7 , 5 ) =  M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
-	C( 7 , 6 ) =  0 ;
-	C( 7 , 7 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
-	C( 7 , 8 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.05E+2 ;
-	C( 7 , 9 ) =  0 ;
-	C( 8 , 0 ) =  3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/8.0 ;
-	C( 8 , 1 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.5E+1 ;
-	C( 8 , 2 ) =  0 ;
-	C( 8 , 3 ) =  0 ;
-	C( 8 , 4 ) =  5.0*M_PI*((3*pow(sin(2*alpha),5)-10*pow(sin(2*alpha),3)+15*sin(2*alpha)) /1.5E+1+(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/1.5E+1 )/1.6E+1 ;
-	C( 8 , 5 ) =  M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin (alpha))/3.0)/1.6E+1 ;
-	C( 8 , 6 ) =  0 ;
-	C( 8 , 7 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.05E+2 ;
-	C( 8 , 8 ) =  3.2E+1*((9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha )/1.92E+2-(9*sin(8*alpha)-4*pow(sin(4*alpha),3)+48*sin(4*alpha)+120 *alpha-60*M_PI)/1.92E+2)/3.5E+1 ;
-	C( 8 , 9 ) =  0 ;
-	C( 9 , 0 ) =  0 ;
-	C( 9 , 1 ) =  0 ;
-	C( 9 , 2 ) =  2.0*(M_PI-alpha)/5.0 ;
-	C( 9 , 3 ) =  M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
-	C( 9 , 4 ) =  0 ;
-	C( 9 , 5 ) =  0 ;
-	C( 9 , 6 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
-	C( 9 , 7 ) =  0 ;
-	C( 9 , 8 ) =  0 ;
-	C( 9 , 9 ) =  2.0*(M_PI-alpha)/7.0 ;
+	C( 6 , 0 ) = 0 ;
+	C( 6 , 1 ) = 0 ;
+	C( 6 , 2 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
+	C( 6 , 3 ) = M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin (alpha))/3.0)/1.6E+1 ;
+	C( 6 , 4 ) = 0 ;
+	C( 6 , 5 ) = 0 ;
+	C( 6 , 6 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.05E+2 ;
+
+	C( 7 , 0 ) = M_PI*(sin(2*alpha)+sin(alpha))/8.0 ;
+	C( 7 , 1 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
+	C( 7 , 2 ) = 0 ;
+	C( 7 , 3 ) = 0 ;
+	C( 7 , 4 ) = M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin (alpha))/3.0)/1.6E+1 ;
+	C( 7 , 5 ) = M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
+	C( 7 , 6 ) = 0 ;
+	C( 7 , 7 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
+
+	C( 8 , 0 ) = 3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/8.0 ;
+	C( 8 , 1 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.5E+1 ;
+	C( 8 , 2 ) = 0 ;
+	C( 8 , 3 ) = 0 ;
+	C( 8 , 4 ) = 5.0*M_PI*((3*pow(sin(2*alpha),5)-10*pow(sin(2*alpha),3)+15*sin(2*alpha)) /1.5E+1+(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/1.5E+1 )/1.6E+1 ;
+	C( 8 , 5 ) = M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin (alpha))/3.0)/1.6E+1 ;
+	C( 8 , 6 ) = 0 ;
+	C( 8 , 7 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.05E+2 ;
+	C( 8 , 8 ) = 3.2E+1*((9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha )/1.92E+2-(9*sin(8*alpha)-4*pow(sin(4*alpha),3)+48*sin(4*alpha)+120 *alpha-60*M_PI)/1.92E+2)/3.5E+1 ;
+
+	C( 9 , 0 ) = 0 ;
+	C( 9 , 1 ) = 0 ;
+	C( 9 , 2 ) = 2.0*(M_PI-alpha)/5.0 ;
+	C( 9 , 3 ) = M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
+	C( 9 , 4 ) = 0 ;
+	C( 9 , 5 ) = 0 ;
+	C( 9 , 6 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
+	C( 9 , 7 ) = 0 ;
+	C( 9 , 8 ) = 0 ;
+	C( 9 , 9 ) = 2.0*(M_PI-alpha)/7.0 ;
 
 	if (size < 11)
 		return C ;
 
-	C( 10 , 0 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
-	C( 10 , 1 ) =  M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin (alpha))/3.0)/1.6E+1 ;
-	C( 10 , 2 ) =  0 ;
-	C( 10 , 3 ) =  0 ;
-	C( 10 , 4 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.05E+2 ;
-	C( 10 , 5 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
-	C( 10 , 6 ) =  0 ;
-	C( 10 , 7 ) =  3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/1.28E+2 ;
-	C( 10 , 8 ) =  5.0*M_PI*((3*pow(sin(2*alpha),5)-10*pow(sin(2*alpha),3)+15*sin(2*alpha)) /1.5E+1+(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/1.5E+1)/1.28E+2 ;
-	C( 10 , 9 ) =  0 ;
-	C( 10 , 10 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/3.15E+2 ;
-	C( 10 , 11 ) =  0 ;
-	C( 10 , 12 ) =  0 ;
-	C( 10 , 13 ) =  3.2E+1*((9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha )/1.92E+2-(9*sin(8*alpha)-4*pow(sin(4*alpha),3)+48*sin(4*alpha)+120 *alpha-60*M_PI)/1.92E+2)/3.15E+2 ;
-	C( 10 , 14 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/6.3E+1 ;
-	C( 11 , 0 ) =  0 ;
-	C( 11 , 1 ) =  0 ;
-	C( 11 , 2 ) =  M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin (alpha))/3.0)/1.6E+1 ;
-	C( 11 , 3 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.05E+2 ;
-	C( 11 , 4 ) =  0 ;
-	C( 11 , 5 ) =  0 ;
-	C( 11 , 6 ) =  5.0*M_PI*((3*pow(sin(2*alpha),5)-10*pow(sin(2*alpha),3)+15*sin(2*alpha)) /1.5E+1+(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/1.5E+1 )/1.28E+2 ;
-	C( 11 , 7 ) =  0 ;
-	C( 11 , 8 ) =  0 ;
-	C( 11 , 9 ) =  3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/1.28E+2 ;
-	C( 11 , 10 ) =  0 ;
-	C( 11 , 11 ) =  3.2E+1*((9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha )/1.92E+2-(9*sin(8*alpha)-4*pow(sin(4*alpha),3)+48*sin(4*alpha)+120 *alpha-60*M_PI)/1.92E+2)/3.15E+2 ;
-	C( 11 , 12 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/3.15E+2 ;
-	C( 11 , 13 ) =  0 ;
-	C( 11 , 14 ) =  0 ;
-	C( 12 , 0 ) =  0 ;
-	C( 12 , 1 ) =  0 ;
-	C( 12 , 2 ) =  M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
-	C( 12 , 3 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
-	C( 12 , 4 ) =  0 ;
-	C( 12 , 5 ) =  0 ;
-	C( 12 , 6 ) =  3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/1.28E+2 ;
-	C( 12 , 7 ) =  0 ;
-	C( 12 , 8 ) =  0 ;
-	C( 12 , 9 ) =  5.0*M_PI*(sin(2*alpha)+sin(alpha))/1.28E+2 ;
-	C( 12 , 10 ) =  0 ;
-	C( 12 , 11 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/3.15E+2 ;
-	C( 12 , 12 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/6.3E+1 ;
-	C( 12 , 13 ) =  0 ;
-	C( 12 , 14 ) =  0 ;
-	C( 13 , 0 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.5E+1 ;
-	C( 13 , 1 ) =  5.0*M_PI*((3*pow(sin(2*alpha),5)-10*pow(sin(2*alpha),3)+15*sin(2*alpha)) /1.5E+1+(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/1.5E+1 )/1.6E+1 ;
-	C( 13 , 2 ) =  0 ;
-	C( 13 , 3 ) =  0 ;
-	C( 13 , 4 ) =  3.2E+1*((9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha )/1.92E+2-(9*sin(8*alpha)-4*pow(sin(4*alpha),3)+48*sin(4*alpha)+120 *alpha-60*M_PI)/1.92E+2)/3.5E+1 ;
-	C( 13 , 5 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.05E+2 ;
-	C( 13 , 6 ) =  0 ;
-	C( 13 , 7 ) =  5.0*M_PI*((3*pow(sin(2*alpha),5)-10*pow(sin(2*alpha),3)+15*sin(2*alpha)) /1.5E+1+(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/1.5E+1 )/1.28E+2 ;
-	C( 13 , 8 ) =  3.5E+1*M_PI*(-(5*pow(sin(2*alpha),7)-21*pow(sin(2*alpha),5)+35*pow(sin(2*alpha),3)-35*sin(2*alpha))/3.5E+1-(5*pow(sin(alpha),7)-21*pow(sin(alpha),5)+35*pow(sin(alpha),3)-35*sin(alpha))/3.5E+1)/1.28E+2 ;
-	C( 13 , 9 ) =  0 ;
-	C( 13 , 10 ) =  3.2E+1*((9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha )/1.92E+2-(9*sin(8*alpha)-4*pow(sin(4*alpha),3)+48*sin(4*alpha)+120 *alpha-60*M_PI)/1.92E+2)/3.15E+2 ;
-	C( 13 , 11 ) =  0 ;
-	C( 13 , 12 ) =  0 ;
-	C( 13 , 13 ) =  2.56E+2*((3*sin(8*alpha)+168*sin(4*alpha)-128*pow(sin(2*alpha),3)+768* sin(2*alpha)+840*alpha)/3.072E+3-(3*sin(16*alpha)+168*sin(8*alpha)-128*pow(sin(4*alpha),3)+768*sin(4*alpha)+1680*alpha-840*M_PI)/3.072E+3)/3.15E+2 ;
-	C( 13 , 14 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/3.15E+2 ;
-	C( 14 , 0 ) =  2.0*(M_PI-alpha)/5.0 ;
-	C( 14 , 1 ) =  M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
-	C( 14 , 2 ) =  0 ;
-	C( 14 , 3 ) =  0 ;
-	C( 14 , 4 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
-	C( 14 , 5 ) =  2.0*(M_PI-alpha)/7.0 ;
-	C( 14 , 6 ) =  0 ;
-	C( 14 , 7 ) =  5.0*M_PI*(sin(2*alpha)+sin(alpha))/1.28E+2 ;
-	C( 14 , 8 ) =  3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/1.28E+2 ;
-	C( 14 , 9 ) =  0 ;
-	C( 14 , 10 ) =  4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/6.3E+1 ;
-	C( 14 , 11 ) =  0 ;
-	C( 14 , 12 ) =  0 ;
-	C( 14 , 13 ) =  1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/3.15E+2 ;
-	C( 14 , 14 ) =  2.0*(M_PI-alpha)/9.0 ;
+	C( 10 , 0 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ;
+	C( 10 , 1 ) = M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin (alpha))/3.0)/1.6E+1 ;
+	C( 10 , 2 ) = 0 ;
+	C( 10 , 3 ) = 0 ;
+	C( 10 , 4 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.05E+2 ;
+	C( 10 , 5 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
+	C( 10 , 6 ) = 0 ;
+	C( 10 , 7 ) = 3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/1.28E+2 ;
+	C( 10 , 8 ) = 5.0*M_PI*((3*pow(sin(2*alpha),5)-10*pow(sin(2*alpha),3)+15*sin(2*alpha)) /1.5E+1+(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/1.5E+1)/1.28E+2 ;
+	C( 10 , 9 ) = 0 ;
+	C( 10 , 10 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/3.15E+2 ;
+
+	C( 11 , 0 ) = 0 ;
+	C( 11 , 1 ) = 0 ;
+	C( 11 , 2 ) = M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin (alpha))/3.0)/1.6E+1 ;
+	C( 11 , 3 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.05E+2 ;
+	C( 11 , 4 ) = 0 ;
+	C( 11 , 5 ) = 0 ;
+	C( 11 , 6 ) = 5.0*M_PI*((3*pow(sin(2*alpha),5)-10*pow(sin(2*alpha),3)+15*sin(2*alpha)) /1.5E+1+(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/1.5E+1 )/1.28E+2 ;
+	C( 11 , 7 ) = 0 ;
+	C( 11 , 8 ) = 0 ;
+	C( 11 , 9 ) = 3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/1.28E+2 ;
+	C( 11 , 10 ) = 0 ;
+	C( 11 , 11 ) = 3.2E+1*((9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha )/1.92E+2-(9*sin(8*alpha)-4*pow(sin(4*alpha),3)+48*sin(4*alpha)+120 *alpha-60*M_PI)/1.92E+2)/3.15E+2 ;
+
+	C( 12 , 0 ) = 0 ;
+	C( 12 , 1 ) = 0 ;
+	C( 12 , 2 ) = M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
+	C( 12 , 3 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
+	C( 12 , 4 ) = 0 ;
+	C( 12 , 5 ) = 0 ;
+	C( 12 , 6 ) = 3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/1.28E+2 ;
+	C( 12 , 7 ) = 0 ;
+	C( 12 , 8 ) = 0 ;
+	C( 12 , 9 ) = 5.0*M_PI*(sin(2*alpha)+sin(alpha))/1.28E+2 ;
+	C( 12 , 10 ) = 0 ;
+	C( 12 , 11 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/3.15E+2 ;
+	C( 12 , 12 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/6.3E+1 ;
+
+	C( 13 , 0 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.5E+1 ;
+	C( 13 , 1 ) = 5.0*M_PI*((3*pow(sin(2*alpha),5)-10*pow(sin(2*alpha),3)+15*sin(2*alpha)) /1.5E+1+(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/1.5E+1 )/1.6E+1 ;
+	C( 13 , 2 ) = 0 ;
+	C( 13 , 3 ) = 0 ;
+	C( 13 , 4 ) = 3.2E+1*((9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha )/1.92E+2-(9*sin(8*alpha)-4*pow(sin(4*alpha),3)+48*sin(4*alpha)+120 *alpha-60*M_PI)/1.92E+2)/3.5E+1 ;
+	C( 13 , 5 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/1.05E+2 ;
+	C( 13 , 6 ) = 0 ;
+	C( 13 , 7 ) = 5.0*M_PI*((3*pow(sin(2*alpha),5)-10*pow(sin(2*alpha),3)+15*sin(2*alpha)) /1.5E+1+(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/1.5E+1 )/1.28E+2 ;
+	C( 13 , 8 ) = 3.5E+1*M_PI*(-(5*pow(sin(2*alpha),7)-21*pow(sin(2*alpha),5)+35*pow(sin(2*alpha),3)-35*sin(2*alpha))/3.5E+1-(5*pow(sin(alpha),7)-21*pow(sin(alpha),5)+35*pow(sin(alpha),3)-35*sin(alpha))/3.5E+1)/1.28E+2 ;
+	C( 13 , 9 ) = 0 ;
+	C( 13 , 10 ) = 3.2E+1*((9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha )/1.92E+2-(9*sin(8*alpha)-4*pow(sin(4*alpha),3)+48*sin(4*alpha)+120 *alpha-60*M_PI)/1.92E+2)/3.15E+2 ;
+	C( 13 , 11 ) = 0 ;
+	C( 13 , 12 ) = 0 ;
+	C( 13 , 13 ) = 2.56E+2*((3*sin(8*alpha)+168*sin(4*alpha)-128*pow(sin(2*alpha),3)+768* sin(2*alpha)+840*alpha)/3.072E+3-(3*sin(16*alpha)+168*sin(8*alpha)-128*pow(sin(4*alpha),3)+768*sin(4*alpha)+1680*alpha-840*M_PI)/3.072E+3)/3.15E+2 ;
+
+	C( 14 , 0 ) = 2.0*(M_PI-alpha)/5.0 ;
+	C( 14 , 1 ) = M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ;
+	C( 14 , 2 ) = 0 ;
+	C( 14 , 3 ) = 0 ;
+	C( 14 , 4 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ;
+	C( 14 , 5 ) = 2.0*(M_PI-alpha)/7.0 ;
+	C( 14 , 6 ) = 0 ;
+	C( 14 , 7 ) = 5.0*M_PI*(sin(2*alpha)+sin(alpha))/1.28E+2 ;
+	C( 14 , 8 ) = 3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/1.28E+2 ;
+	C( 14 , 9 ) = 0 ;
+	C( 14 , 10 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/6.3E+1 ;
+	C( 14 , 11 ) = 0 ;
+	C( 14 , 12 ) = 0 ;
+	C( 14 , 13 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/3.15E+2 ;
+	C( 14 , 14 ) = 2.0*(M_PI-alpha)/9.0 ;
 
 	return C ;
 }
 
-template <typename REAL>
-void
-QuadricHF<REAL>::fillTensor(Geom::Tensor3d& T)
-{
-	std::vector<unsigned int> p ;
-	p.resize(T.order(), 0) ;
-	do
-	{
-		//		for (unsigned int i = 0 ; i < p.size() ; ++i)
-			//			std::cout << p[i] << " " ;
-		//		std::cout << std::endl ;
-		if (T(p) == CONST_VAL)
-		{
-			std::vector<unsigned int> sorted_p = p ;
-			std::sort(sorted_p.begin(), sorted_p.begin() + T.order()) ;
-			//			std::cout << "sort: " ;
-			//			for (unsigned int i = 0 ; i < p.size() ; ++i)
-				//				std::cout << sorted_p[i] << " " ;
-			//			std::cout << std::endl ;
-			T(p) = T(sorted_p) ;
-		}
-	} while (Geom::Tensor3d::incremIndex(p)) ;
-}
-
 template <typename REAL>
 Geom::Tensor3d*
 QuadricHF<REAL>::tensorsFromCoefs(const std::vector<VEC3>& coefs)
@@ -1120,7 +1154,7 @@ QuadricHF<REAL>::tensorsFromCoefs(const std::vector<VEC3>& coefs)
 	for (unsigned int col = 0 ; col < 3 ; ++col)
 	{
 		A[col] = Geom::Tensor3d(degree) ;
-		A[col].setConst(CONST_VAL) ;
+		//A[col].setConst(CONST_VAL) ;
 	}
 
 	std::vector<unsigned int> index ;
@@ -1193,7 +1227,7 @@ QuadricHF<REAL>::tensorsFromCoefs(const std::vector<VEC3>& coefs)
 	}
 
 	for (unsigned int col = 0 ; col < 3 ; ++col)
-		fillTensor(A[col]) ;
+		A[col].completeSymmetricTensor() ;
 
 	return A ;
 }