From a89fcd09ac3ca39d72d7b4e3e66cb0890672107d Mon Sep 17 00:00:00 2001
From: Eduardo Rojas <eduardo.rojas@earthsoft.com>
Date: Wed, 15 Dec 2021 11:42:19 -0300
Subject: [PATCH] Added function go get points where cubic spline derivative is
 0 (GetHorizontalPoints) and also to get the maximum value of the cubic spline
 function within its domain (GetMaxPoints)

---
 src/Numerics/Interpolation/CubicSpline.cs | 70 +++++++++++++++++++++++
 1 file changed, 70 insertions(+)

diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs
index 1830da86..f3446667 100644
--- a/src/Numerics/Interpolation/CubicSpline.cs
+++ b/src/Numerics/Interpolation/CubicSpline.cs
@@ -623,5 +623,75 @@ namespace MathNet.Numerics.Interpolation
 
             return Math.Min(Math.Max(index, 0), _x.Length - 2);
         }
+
+        /// <summary>
+        /// Gets all the points (x values) where the derivative is 0
+        /// </summary>
+        /// <returns>An array of X values where the derivative of the spline is 0</returns>
+        public double[] GetHorizontalPoints()
+        {
+            List<double> points = new List<double>();
+            for (int index = 0; index < _x.Length-1; index++)
+            {
+                double a = 6 * _c3[index]; //derive ax^3 and multiply by 2
+                double b = 2 * _c2[index]; //derive bx^2
+                double c = _c1[index];//derive cx
+                double d = b * b - 2 * a * c;
+                //first check if a is 0, if so its a linear function, this happens with quadratic condition
+                if (a.AlmostEqual(0))
+                {
+                double x = _x[index] - c / b;
+                //check if the result is in the domain
+                if (_x[index] <= x && x <= _x[index + 1]) points.Add(x);
+                }
+                else if (d.AlmostEqual(0))//its a quadratic with a single solution
+                {
+                double x = _x[index] - b / a;
+                if (_x[index] <= x && x <= _x[index + 1]) points.Add(x);
+                }
+                else if (d > 0)//only has a solution if d is greater than 0
+                {
+                d = (double)System.Math.Sqrt(d);
+                //apply quadratic equation
+                double x1 = _x[index]+ (-b + d) / a;
+                double x2 = _x[index]+ (-b - d) / a;
+                //Add any solution points that fall within the domain to the list
+                if ((_x[index] <= x1) && (x1 <= _x[index + 1])) points.Add(x1);
+                if ((_x[index] <= x2) && (x2 <= _x[index + 1])) points.Add(x2);
+                }
+            }
+            return points.ToArray();
+        }
+
+        /// <summary>
+        /// Returns the minimum value for the spline within its domain.
+        /// </summary>
+        /// <returns></returns>
+        public double GetMaxPoint()
+        {
+            double max = double.MinValue;
+            double x=0;
+            //go through the functions points to check if one of them has a higher value
+            foreach (double p in _x)
+            {
+                double y = Interpolate(p);
+                if (y > max)
+                {
+                max = y;
+                x = p;
+                }
+            }
+            //go through the inflexion, local minimums and local maximums
+            foreach (double p in GetHorizontalPoints())
+            {
+                double y = Interpolate(p);
+                if (y > max)
+                {
+                max = y;
+                x = p;
+                }
+            }
+            return x;
+        }
     }
 }