From bb1e00e474a38dc609abb7ff4e4f637cbd9b3754 Mon Sep 17 00:00:00 2001
From: Richard Brice <37087370+RickBrice@users.noreply.github.com>
Date: Tue, 18 Feb 2025 12:43:11 -0800
Subject: [PATCH] Fixes problems mapping vertical alignment segments business
 logic to geometry

---
 src/ifcparse/IfcAlignmentHelper.cpp | 39 ++++++++++++++++++++++-------
 1 file changed, 30 insertions(+), 9 deletions(-)

diff --git a/src/ifcparse/IfcAlignmentHelper.cpp b/src/ifcparse/IfcAlignmentHelper.cpp
index 8e26375ab5..f60feacc3e 100644
--- a/src/ifcparse/IfcAlignmentHelper.cpp
+++ b/src/ifcparse/IfcAlignmentHelper.cpp
@@ -30,6 +30,7 @@
 // @todo use std::numbers::pi when upgrading to C++ 20
 static const double PI = boost::math::constants::pi<double>();
 
+#include <boost/math/quadrature/trapezoidal.hpp>
 
 #ifdef HAS_SCHEMA_4x3_add2
 
@@ -243,14 +244,14 @@ std::tuple<typename aggregate_of<Ifc4x3_add2::IfcObjectDefinition>::ptr, typenam
         // back gradient
         auto dxBG = xPVI - xPBG;
         auto dyBG = yPVI - yPBG;
-        auto start_slope = atan2(dyBG, dxBG);
+        auto start_slope = tan(atan2(dyBG,dxBG));
 
         // forward gradient
         point_iter++;
         std::tie(xPFG, yPFG) = *point_iter;
         auto dxFG = xPFG - xPVI;
         auto dyFG = yPFG - yPVI;
-        auto end_slope = atan2(dyFG, dxFG);
+        auto end_slope = tan(atan2(dyFG,dxFG));
 
         double xEVC = xPVI + length / 2;
         double yEVC = yPVI + end_slope * length / 2;
@@ -292,8 +293,9 @@ std::tuple<typename aggregate_of<Ifc4x3_add2::IfcObjectDefinition>::ptr, typenam
     // create last tangent run
     auto dx = xPVI - xPBG;
     auto dy = yPVI - yPBG;
-    auto slope = atan2(dy, dx);
-    auto gradient_length = sqrt(dx * dx + dy * dy);
+    auto slope = tan(atan2(dy,dx));
+    auto gradient_length = dx;
+
     file.addDoublet<Ifc4x3_add2::IfcCartesianPoint>(xPBG, yPBG);
     auto design_parameters = new Ifc4x3_add2::IfcAlignmentVerticalSegment(boost::none, boost::none, xPBG, gradient_length, yPBG, slope, slope, boost::none, Ifc4x3_add2::IfcAlignmentVerticalSegmentTypeEnum::IfcAlignmentVerticalSegmentType_CONSTANTGRADIENT);
     auto alignment_segment = new Ifc4x3_add2::IfcAlignmentSegment(IfcParse::IfcGlobalId(), nullptr, boost::none, boost::none, boost::none, nullptr, nullptr, design_parameters);
@@ -688,13 +690,18 @@ std::pair<Ifc4x3_add2::IfcCurveSegment*, Ifc4x3_add2::IfcCurveSegment*> mapAlign
             new Ifc4x3_add2::IfcCartesianPoint(std::vector<double>({0, 0})),
             new Ifc4x3_add2::IfcVector(new Ifc4x3_add2::IfcDirection(std::vector<double>{1, 0}), 1.0));
 
+        // IfcCurveSegment.SegmentLength is the length of the curve segment, not the horizontal length.
+        auto dx = cos(atan(start_gradient));
+        auto dy = sin(atan(start_gradient));
+        auto segment_curve_length = horizontal_length / dx;
+
         auto curve_segment = new Ifc4x3_add2::IfcCurveSegment(
             Ifc4x3_add2::IfcTransitionCode::IfcTransitionCode_CONTSAMEGRADIENT,
             new Ifc4x3_add2::IfcAxis2Placement2D(
                 new Ifc4x3_add2::IfcCartesianPoint({start_distance_along, start_height}),
-                new Ifc4x3_add2::IfcDirection({sqrt(1.0 - start_gradient * start_gradient), start_gradient})),
+                new Ifc4x3_add2::IfcDirection({dx,dy})),
             new Ifc4x3_add2::IfcLengthMeasure(0.0), // start
-            new Ifc4x3_add2::IfcLengthMeasure(horizontal_length),
+            new Ifc4x3_add2::IfcLengthMeasure(segment_curve_length),
             parent_curve);
 
         result.first = curve_segment;
@@ -710,11 +717,22 @@ std::pair<Ifc4x3_add2::IfcCurveSegment*, Ifc4x3_add2::IfcCurveSegment*> mapAlign
             std::vector<double>{A, B, C},
             boost::none);
 
+        // IfcCurveSegment.SegmentLength is the length of the curve segment, not the horizontal length.
+        // The curve length is calculated by integrating the differential curve length equation sqrt(1 + (dy/dx)^2) from 0 to horizontal_length.
+        // y = A + Bx + Cx^2
+        // dy/dx = B + 2Cx
+        auto dx = cos(atan(start_gradient));
+        auto dy = sin(atan(start_gradient));
+        auto curve_length_fn = [B, C](double x) { return sqrt(1 + pow(B + C * x, 2)); };
+        auto segment_curve_length = boost::math::quadrature::trapezoidal(curve_length_fn, 0.0, horizontal_length);
+
         auto curve_segment = new Ifc4x3_add2::IfcCurveSegment(
             Ifc4x3_add2::IfcTransitionCode::IfcTransitionCode_CONTSAMEGRADIENT,
-            new Ifc4x3_add2::IfcAxis2Placement2D(new Ifc4x3_add2::IfcCartesianPoint({start_distance_along, start_height}), new Ifc4x3_add2::IfcDirection({sqrt(1.0 - start_gradient * start_gradient), start_gradient})),
+            new Ifc4x3_add2::IfcAxis2Placement2D(
+               new Ifc4x3_add2::IfcCartesianPoint({start_distance_along, start_height}), 
+               new Ifc4x3_add2::IfcDirection({dx,dy})),
             new Ifc4x3_add2::IfcLengthMeasure(0.0),
-            new Ifc4x3_add2::IfcLengthMeasure(horizontal_length),
+            new Ifc4x3_add2::IfcLengthMeasure(segment_curve_length),
             parent_curve);
 
         result.first = curve_segment;
@@ -735,11 +753,14 @@ std::pair<Ifc4x3_add2::IfcCurveSegment*, Ifc4x3_add2::IfcCurveSegment*> mapAlign
                                                  new Ifc4x3_add2::IfcDirection(std::vector<double>{1, 0})),
             radius);
 
+
+        auto segment_curve_length = radius * fabs(end_angle - start_angle);
+
         Ifc4x3_add2::IfcCurveSegment* curve_segment = new Ifc4x3_add2::IfcCurveSegment(
             Ifc4x3_add2::IfcTransitionCode::IfcTransitionCode_CONTSAMEGRADIENT,
             new Ifc4x3_add2::IfcAxis2Placement2D(new Ifc4x3_add2::IfcCartesianPoint({start_distance_along, start_height}), new Ifc4x3_add2::IfcDirection({1.0, 0.0})),
             new Ifc4x3_add2::IfcLengthMeasure(0.0),
-            new Ifc4x3_add2::IfcLengthMeasure(horizontal_length),
+            new Ifc4x3_add2::IfcLengthMeasure(segment_curve_length),
             parent_curve);
 
         result.first = curve_segment;