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);
}
+
+ ///
+ /// Gets all the points (x values) where the derivative is 0
+ ///
+ /// An array of X values where the derivative of the spline is 0
+ public double[] GetHorizontalPoints()
+ {
+ List points = new List();
+ 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();
+ }
+
+ ///
+ /// Returns the minimum value for the spline within its domain.
+ ///
+ ///
+ 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;
+ }
}
}