From 7b2be41d85f28bde6dca44faf60e6bf64dabe27f Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (o)" <j.wuttke@fz-juelich.de>
Date: Mon, 12 Dec 2022 18:09:10 +0100
Subject: [PATCH] Pyramid2 ditto

---
 Sample/HardParticle/Pyramid2.cpp | 24 +++++++++++-------------
 1 file changed, 11 insertions(+), 13 deletions(-)

diff --git a/Sample/HardParticle/Pyramid2.cpp b/Sample/HardParticle/Pyramid2.cpp
index 3b42525201a..3ad3ad53e24 100644
--- a/Sample/HardParticle/Pyramid2.cpp
+++ b/Sample/HardParticle/Pyramid2.cpp
@@ -49,23 +49,21 @@ std::string Pyramid2::validate() const
     if (!errs.empty())
         return jointError(errs);
 
-    const double cot_alpha = Math::cot(m_alpha);
-    double r = cot_alpha * 2 * m_height / m_length;
-    double s = cot_alpha * 2 * m_height / m_width;
-    if (r > 1)
-        errs.push_back("parameters violate condition 2 * ctg(alpha) <= m_height / m_length");
-    if (s > 1)
-        errs.push_back("parameters violate condition 2 * ctg(alpha) <= m_height / m_width");
+    const double D = m_length / 2;
+    const double W = m_width / 2;
+    const double R = m_height * Math::cot(m_alpha);
+    if (R > D)
+        errs.push_back("parameters violate condition H * ctg(alpha) <= m_length / 2");
+    if (R > W)
+        errs.push_back("parameters violate condition H * ctg(alpha) <= m_width / 2");
     if (!errs.empty())
         return jointError(errs);
 
-    double D = m_length / 2;
-    double d = m_length / 2 * (1 - r);
-    double W = m_width / 2;
-    double w = m_width / 2 * (1 - s);
-
     // center of mass
-    double zcom = m_height * 0.5 * (6 - 4 * (r + s) + 3 * r * s) / (6 - 3 * (r + s) + 2 * r * s);
+    double zcom = m_height / 2 * (6 * W * D - 4 * (W + D) * R + 3 * R * R)
+                  / (6 * W * D - 3 * (W + D) * R + 2 * R * R);
+    const double d = D - R;
+    const double w = W - R;
 
     setPolyhedron(topology, -zcom,
                   {// base:
-- 
GitLab