FormFactorPolygon now requests origin at center of mass.

AnisoPyramid broken.

Core/FormFactors/FormFactorAnisoPyramid.cpp

@@ -52,7 +52,9 @@ void FormFactorAnisoPyramid::onChange()
     double cot_alpha = MathFunctions::cot(m_alpha);
     if( !std::isfinite(cot_alpha) || cot_alpha<0 )
         throw Exceptions::OutOfBoundsException("pyramid angle alpha out of bounds");
-    if(cot_alpha*m_height > m_length || cot_alpha*m_height > m_width) {
+    double r = cot_alpha*2 * m_height / m_length;
+    double s = cot_alpha*2 * m_height / m_width;
+    if( r>1 || s>1 ) {
         std::ostringstream ostr;
         ostr << "FormFactorAnisoPyramid() -> Error in class initialization with parameters";
         ostr << " length:" << m_length;
@@ -64,21 +66,23 @@ void FormFactorAnisoPyramid::onChange()
     }
 
     double D = m_length/2;
-    double d = m_length/2 - m_height*cot_alpha;
+    double d = m_length/2 * (1-r);
     double W = m_width/2;
-    double w = m_width/2 - m_height*cot_alpha;
+    double w = m_width/2 * (1-s);
 
-    setPolyhedron( topology, 0, false, {
+    double zcom = m_height * ( .5 - 2*r/3 + r*r/4 ) / ( 1 - r + r*r/3 ); // center of mass
+
+    setPolyhedron( topology, -zcom, false, {
         // base:
-        { -D, -W, 0. },
-        {  D, -W, 0. },
-        {  D,  W, 0. },
-        { -D,  W, 0. },
+        { -D, -W, -zcom },
+        {  D, -W, -zcom },
+        {  D,  W, -zcom },
+        { -D,  W, -zcom },
         // top:
-        { -d, -w, m_height },
-        {  d, -w, m_height },
-        {  d,  w, m_height },
-        { -d,  w, m_height } } );
+        { -d, -w, m_height-zcom },
+        {  d, -w, m_height-zcom },
+        {  d,  w, m_height-zcom },
+        { -d,  w, m_height-zcom } } );
 }
 
 FormFactorAnisoPyramid* FormFactorAnisoPyramid::clone() const

Core/FormFactors/FormFactorCone6.cpp

@@ -47,38 +47,41 @@ void FormFactorCone6::onChange()
     double cot_alpha = MathFunctions::cot(m_alpha);
     if( !std::isfinite(cot_alpha) || cot_alpha<0 )
         throw Exceptions::OutOfBoundsException("pyramid angle alpha out of bounds");
-    if(cot_alpha*m_height > m_base_edge) {
+    double r = cot_alpha*2/sqrt(3) * m_height / m_base_edge; // L(top)/L(base)
+    if ( r > 1 ) {
         std::ostringstream ostr;
         ostr << "FormFactorCone6() -> Error in class initialization with parameters";
         ostr << " base_edge:" << m_base_edge;
         ostr << " height:" << m_height;
         ostr << " alpha[rad]:" << m_alpha << "\n\n";
-        ostr << "Check for 'height <= base_edge*tan(alpha)' failed.";
+        ostr << "Incompatible parameters\n";
         throw Exceptions::ClassInitializationException(ostr.str());
     }
 
     double a = m_base_edge;
     double as = a/2;
     double ac = a*sqrt(3)/2;
-    double b = m_base_edge - 2*m_height/sqrt(3)*cot_alpha;
+    double b = a * (1-r);
     double bs = b/2;
     double bc = b*sqrt(3)/2;
 
-    setPolyhedron( topology, 0, false, {
+    double zcom = m_height * ( .5 - 2*r/3 + r*r/4 ) / ( 1 - r + r*r/3 ); // center of mass
+
+    setPolyhedron( topology, -zcom, false, {
         // base:
-        {  a,   0., 0. },
-        {  as,  ac, 0. },
-        { -as,  ac, 0. },
-        { -a,   0., 0. },
-        { -as, -ac, 0. },
-        {  as, -ac, 0. },
+        {  a,   0., -zcom },
+        {  as,  ac, -zcom },
+        { -as,  ac, -zcom },
+        { -a,   0., -zcom },
+        { -as, -ac, -zcom },
+        {  as, -ac, -zcom },
         // top:
-        {  b,   0., m_height },
-        {  bs,  bc, m_height },
-        { -bs,  bc, m_height },
-        { -b,   0., m_height },
-        { -bs, -bc, m_height },
-        {  bs, -bc, m_height } } );
+        {  b,   0., m_height-zcom },
+        {  bs,  bc, m_height-zcom },
+        { -bs,  bc, m_height-zcom },
+        { -b,   0., m_height-zcom },
+        { -bs, -bc, m_height-zcom },
+        {  bs, -bc, m_height-zcom } } );
 }
 
 FormFactorCone6* FormFactorCone6::clone() const

Core/FormFactors/FormFactorPolyhedron.cpp

@@ -152,7 +152,7 @@ PolyhedralFace::PolyhedralFace( const std::vector<kvector_t>& V, bool _sym_S2 )
         }
         m_normal += ee.unit();
     }
-    m_normal /= NE;
+    m_normal = m_normal.unit();
     m_rperp = 0;
     for( size_t j=0; j<NV; ++j )
         m_rperp += V[j].dot(m_normal);
@@ -161,12 +161,16 @@ PolyhedralFace::PolyhedralFace( const std::vector<kvector_t>& V, bool _sym_S2 )
     for ( size_t j=1; j<NV; ++j )
         if( std::abs(V[j].dot(m_normal) - m_rperp) > 1e-14*m_radius_3d )
             throw std::runtime_error("Face is not planar");
-    // compute m_area
+    // compute area and center of gravity
     m_area = 0;
+    m_center = kvector_t();
     for ( size_t j=0; j<NV; ++j ) {
         size_t jj = (j+1)%NV;
-        m_area += m_normal.dot( V[j].cross( V[jj] ) ) / 2;
+        double wgt = m_normal.dot( V[j].cross( V[jj] ) ) / 2;
+        m_area += wgt;
+        m_center += wgt * ( m_rperp*m_normal+V[j]+V[jj] )/3;
     }
+    m_center /= 4*m_area;
     // only now deal with inversion symmetry
     if( sym_S2 ) {
         if( NE&1 )
@@ -385,9 +389,16 @@ void FormFactorPolyhedron::setPolyhedron(
 
     m_radius = 0;
     m_volume = 0;
+    kvector_t center;
     for( const PolyhedralFace& Gk: m_faces ) {
         m_radius = std::max( m_radius, Gk.radius3d() );
         m_volume += Gk.pyramidalVolume();
+        center += Gk.center() * Gk.pyramidalVolume();
+    }
+    center /= m_volume;
+    if( center.mag() > 1e-14*m_radius ) {
+        std::cerr << "center of mass: " << center << "\n";
+        throw std::runtime_error("Center of mass is not at origin");
     }
 
     if( m_sym_Ci ) {

Core/FormFactors/FormFactorPolyhedron.h

@@ -57,6 +57,7 @@ public:
     PolyhedralFace( const std::vector<kvector_t>& _V=std::vector<kvector_t>(), bool _sym_S2=false );
 
     double area() const { return m_area; }
+    kvector_t center() const { return m_center; }
     double pyramidalVolume() const { return m_rperp*m_area/3; }
     double radius3d() const { return m_radius_3d; }
     complex_t ff_n( int m, const cvector_t q ) const;
@@ -75,6 +76,7 @@ private:
     double m_rperp; //!< distance of this polygon's plane from the origin, along 'm_normal'
     double m_radius_2d; //!< radius of enclosing cylinder
     double m_radius_3d; //!< radius of enclosing sphere
+    kvector_t m_center; //!< center of mass
 
     void decompose_q( const cvector_t q, complex_t& qperp, cvector_t& qpa ) const;
     complex_t ff_n_core( int m, const cvector_t qpa ) const;

Core/FormFactors/FormFactorPyramid.cpp

@@ -49,7 +49,8 @@ void FormFactorPyramid::onChange()
     double cot_alpha = MathFunctions::cot(m_alpha);
     if( !std::isfinite(cot_alpha) || cot_alpha<0 )
         throw Exceptions::OutOfBoundsException("pyramid angle alpha out of bounds");
-    if(cot_alpha*m_height > m_base_edge) {
+    double r = cot_alpha * m_height / m_base_edge; // L(top)/L(base)
+    if ( r > 1 ) {
         std::ostringstream ostr;
         ostr << "FormFactorPyramid() -> Error in class initialization with parameters";
         ostr << " base_edge:" << m_base_edge;
@@ -60,19 +61,21 @@ void FormFactorPyramid::onChange()
     }
 
     double a = m_base_edge/2;
-    double b = m_base_edge/2 - m_height*cot_alpha;
+    double b = a * (1-r);
 
-    setPolyhedron( topology, 0, false, {
-        // base:
-        { -a, -a, 0. },
-        {  a, -a, 0. },
-        {  a,  a, 0. },
-        { -a,  a, 0. },
-        // top:
-        { -b, -b, m_height },
-        {  b, -b, m_height },
-        {  b,  b, m_height },
-        { -b,  b, m_height } } );
+    double zcom = m_height * ( .5 - 2*r/3 + r*r/4 ) / ( 1 - r + r*r/3 ); // center of mass
+
+    setPolyhedron( topology, -zcom, false, {
+            // base:
+            { -a, -a, -zcom },
+            {  a, -a, -zcom },
+            {  a,  a, -zcom },
+            { -a,  a, -zcom },
+            // top:
+            { -b, -b, m_height-zcom },
+            {  b, -b, m_height-zcom },
+            {  b,  b, m_height-zcom },
+            { -b,  b, m_height-zcom } } );
 }
 
 FormFactorPyramid* FormFactorPyramid::clone() const

Core/FormFactors/FormFactorTetrahedron.cpp

@@ -50,7 +50,8 @@ void FormFactorTetrahedron::onChange()
     double cot_alpha = MathFunctions::cot(m_alpha);
     if( !std::isfinite(cot_alpha) || cot_alpha<0 )
         throw Exceptions::OutOfBoundsException("pyramid angle alpha out of bounds");
-    if (cot_alpha * 2*std::sqrt(3.) * m_height > m_base_edge) {
+    double r = cot_alpha * 2*std::sqrt(3.) * m_height / m_base_edge; // L(top)/L(base)
+    if ( r > 1 ) {
         std::ostringstream ostr;
         ostr << "FormFactorTetrahedron() -> Error in class initialization with parameters ";
         ostr << " height:" << m_height;
@@ -64,20 +65,22 @@ void FormFactorTetrahedron::onChange()
     double as = a/2;
     double ac = a/sqrt(3)/2;
     double ah = a/sqrt(3);
-    double b = a - 2*sqrt(3)*m_height*cot_alpha;
+    double b = a * (1-r);
     double bs = b/2;
     double bc = b/sqrt(3)/2;
     double bh = b/sqrt(3);
 
-    setPolyhedron( topology, 0, false, {
+    double zcom = m_height * ( .5 - 2*r/3 + r*r/4 ) / ( 1 - r + r*r/3 ); // center of mass
+
+    setPolyhedron( topology, -zcom, false, {
             // base:
-            { -as, -ac, 0. },
-            {  as, -ac, 0. },
-            {  0.,  ah, 0. },
+            { -as, -ac, -zcom },
+            {  as, -ac, -zcom },
+            {  0.,  ah, -zcom },
             // top:
-            { -bs, -bc, m_height },
-            {  bs, -bc, m_height },
-            {  0.,  bh, m_height } } );
+            { -bs, -bc, m_height-zcom },
+            {  bs, -bc, m_height-zcom },
+            {  0.,  bh, m_height-zcom } } );
 }
 
 FormFactorTetrahedron* FormFactorTetrahedron::clone() const