diff --git a/src/Numerics.Tests/SpecialFunctionsTests/KelvinTests.cs b/src/Numerics.Tests/SpecialFunctionsTests/KelvinTests.cs
new file mode 100644
index 00000000..96857265
--- /dev/null
+++ b/src/Numerics.Tests/SpecialFunctionsTests/KelvinTests.cs
@@ -0,0 +1,135 @@
+using NUnit.Framework;
+using System;
+
+namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
+{
+ ///
+ /// Kelvin functions tests.
+ ///
+ [TestFixture, Category("Functions")]
+ public class KelvinTests
+ {
+ [Test]
+ public void KelvinBerApprox([Range(-8, 8, 0.25)] double x)
+ {
+ // Approx by Abramowitz/Stegun 9.11.1
+ Assert.AreEqual(Polynomial.Evaluate(x / 8.0,
+ 1.0,
+ 0.0, 0.0, 0.0, -64.0, 0.0,
+ 0.0, 0.0, 113.77777774, 0.0, 0.0,
+ 0.0, -32.36345652, 0.0, 0.0, 0.0,
+ 2.64191397, 0.0, 0.0, 0.0, -0.08349609,
+ 0.0, 0.0, 0.0, 0.00122552, 0.0,
+ 0.0, 0.0, -0.00000901), SpecialFunctions.KelvinBer(x), 1e-9);
+ }
+
+ [Test]
+ public void KelvinBeiApprox([Range(-8, 8, 0.25)] double x)
+ {
+ // Approx by Abramowitz/Stegun 9.11.2
+ Assert.AreEqual(Polynomial.Evaluate(x / 8.0,
+ 0.0,
+ 0.0, 16.0, 0.0, 0.0, 0.0,
+ -113.77777774, 0.0, 0.0, 0.0, 72.81777742,
+ 0.0, 0.0, 0.0, -10.56765779, 0.0,
+ 0.0, 0.0, 0.52185615, 0.0, 0.0,
+ 0.0, -0.01103667, 0.0, 0.0, 0.0,
+ 0.00011346),
+ SpecialFunctions.KelvinBei(x), 6e-9);
+ }
+
+ [Test]
+ public void KelvinKerApprox([Range(0.25, 8, 0.25)] double x)
+ {
+ // Approx by Abramowitz/Stegun 9.11.3
+ Assert.AreEqual(
+ Polynomial.Evaluate(x / 8.0,
+ -Math.Log(x / 2.0) * SpecialFunctions.KelvinBer(x) + SpecialFunctions.KelvinBei(x) * Constants.PiOver4 - 0.57721566,
+ 0.0, 0.0, 0.0, -59.05819744, 0.0,
+ 0.0, 0.0, 171.36272133, 0.0, 0.0,
+ 0.0, -60.60977451, 0.0, 0.0, 0.0,
+ 5.65539121, 0.0, 0.0, 0.0, -0.19636347,
+ 0.0, 0.0, 0.0, 0.00309699, 0.0,
+ 0.0, 0.0, -0.00002458),
+ SpecialFunctions.KelvinKer(x), 1e-8);
+ }
+
+ [Test]
+ public void KelvinKeiApprox([Range(0.25, 8, 0.25)] double x)
+ {
+ // Approx by Abramowitz/Stegun 9.11.4
+ Assert.AreEqual(
+ -Math.Log(x / 2.0) * SpecialFunctions.KelvinBei(x) - Constants.PiOver4 * SpecialFunctions.KelvinBer(x)
+ + Polynomial.Evaluate(x / 8.0,
+ 0.0,
+ 0.0, 6.76454936, 0.0, 0.0, 0.0,
+ -142.91827687, 0.0, 0.0, 0.0, 124.23569650,
+ 0.0, 0.0, 0.0, -21.30060904, 0.0,
+ 0.0, 0.0, 1.17509064, 0.0, 0.0,
+ 0.0, -0.02695875, 0.0, 0.0, 0.0,
+ 0.00029532),
+ SpecialFunctions.KelvinKei(x), 3e-9);
+ }
+
+ [Test]
+ public void KelvinBerPrimeApprox([Range(-8, 8, 0.25)] double x)
+ {
+ // Approx by Abramowitz/Stegun 9.11.5
+ Assert.AreEqual(x * Polynomial.Evaluate(x / 8.0,
+ 0.0,
+ 0.0, -4.0, 0.0, 0.0, 0.0,
+ 14.22222222, 0.0, 0.0, 0.0, -6.06814810,
+ 0.0, 0.0, 0.0, 0.66047849, 0.0,
+ 0.0, 0.0, -0.02609253, 0.0, 0.0,
+ 0.0, 0.00045957, 0.0, 0.0, 0.0,
+ -0.00000394), SpecialFunctions.KelvinBerPrime(x), 2.1e-8);
+ }
+
+ [Test]
+ public void KelvinBeiPrimeApprox([Range(-8, 8, 0.25)] double x)
+ {
+ // Approx by Abramowitz/Stegun 9.11.6
+ Assert.AreEqual(x * Polynomial.Evaluate(x / 8.0,
+ 0.5,
+ 0.0, 0.0, 0.0, -10.66666666, 0.0,
+ 0.0, 0.0, 11.37777772, 0.0, 0.0,
+ 0.0, -2.31167514, 0.0, 0.0, 0.0,
+ 0.14677204, 0.0, 0.0, 0.0, -0.00379386,
+ 0.0, 0.0, 0.0, 0.00004609),
+ SpecialFunctions.KelvinBeiPrime(x), 7e-8);
+ }
+
+ [Test]
+ public void KelvinKerPrimeApprox([Range(0.25, 8, 0.25)] double x)
+ {
+ // Approx by Abramowitz/Stegun 9.11.7
+ Assert.AreEqual(
+ -Math.Log(x / 2.0) * SpecialFunctions.KelvinBerPrime(x) - SpecialFunctions.KelvinBer(x) / x + Constants.PiOver4 * SpecialFunctions.KelvinBeiPrime(x)
+ + x * Polynomial.Evaluate(x / 8.0,
+ 0.0,
+ 0.0, -3.69113734, 0.0, 0.0, 0.0,
+ 21.42034017, 0.0, 0.0, 0.0, -11.36433272,
+ 0.0, 0.0, 0.0, 1.41384780, 0.0,
+ 0.0, 0.0, -0.06136358, 0.0, 0.0,
+ 0.0, 0.00116137, 0.0, 0.0, 0.0,
+ -0.00001075),
+ SpecialFunctions.KelvinKerPrime(x), 8e-8);
+ }
+
+ [Test]
+ public void KelvinKeiPrimeApprox([Range(0.25, 8, 0.25)] double x)
+ {
+ // Approx by Abramowitz/Stegun 9.11.8
+ Assert.AreEqual(
+ -Math.Log(x / 2.0) * SpecialFunctions.KelvinBeiPrime(x) - SpecialFunctions.KelvinBei(x) / x - Constants.PiOver4 * SpecialFunctions.KelvinBerPrime(x)
+ + x * Polynomial.Evaluate(x / 8.0,
+ 0.21139217,
+ 0.0, 0.0, 0.0, -13.39858846, 0.0,
+ 0.0, 0.0, 19.41182758, 0.0, 0.0,
+ 0.0, -4.65950823, 0.0, 0.0, 0.0,
+ 0.33049424, 0.0, 0.0, 0.0, -0.00926707,
+ 0.0, 0.0, 0.0, 0.00011997),
+ SpecialFunctions.KelvinKeiPrime(x), 7e-8);
+ }
+ }
+}
diff --git a/src/Numerics/SpecialFunctions/Kelvin.cs b/src/Numerics/SpecialFunctions/Kelvin.cs
new file mode 100644
index 00000000..9c6e467a
--- /dev/null
+++ b/src/Numerics/SpecialFunctions/Kelvin.cs
@@ -0,0 +1,264 @@
+using System;
+using System.Numerics;
+
+namespace MathNet.Numerics
+{
+ ///
+ /// This partial implementation of the SpecialFunctions class contains all methods related to the modified Bessel function.
+ ///
+ public static partial class SpecialFunctions
+ {
+ ///
+ /// Returns the Kelvin function of the first kind.
+ /// KelvinBe(nu, x) is given by BesselJ(0, j * sqrt(j) * x) where j = sqrt(-1).
+ /// KelvinBer(nu, x) and KelvinBei(nu, x) are the real and imaginary parts of the KelvinBe(nu, x)
+ ///
+ /// the order of the the Kelvin function.
+ /// The value to compute the Kelvin function of.
+ /// The Kelvin function of the first kind.
+ public static Complex KelvinBe(double nu, double x)
+ {
+ Complex ISqrtI = new Complex(-Constants.Sqrt1Over2, Constants.Sqrt1Over2); // j * sqrt(j) = (-1)^(3/4) = (-1 + j)/sqrt(2)
+ return BesselJ(nu, ISqrtI * x);
+ }
+
+ ///
+ /// Returns the Kelvin function ber.
+ /// KelvinBer(nu, x) is given by the real part of BesselJ(nu, j * sqrt(j) * x) where j = sqrt(-1).
+ ///
+ /// the order of the the Kelvin function.
+ /// The value to compute the Kelvin function of.
+ /// The Kelvin function ber.
+ public static double KelvinBer(double nu, double x)
+ {
+ return KelvinBe(nu, x).Real;
+ }
+
+ ///
+ /// Returns the Kelvin function ber.
+ /// KelvinBer(x) is given by the real part of BesselJ(0, j * sqrt(j) * x) where j = sqrt(-1).
+ /// KelvinBer(x) is equivalent to KelvinBer(0, x).
+ ///
+ /// The value to compute the Kelvin function of.
+ /// The Kelvin function ber.
+ public static double KelvinBer(double x)
+ {
+ return KelvinBe(0, x).Real;
+ }
+
+ ///
+ /// Returns the Kelvin function bei.
+ /// KelvinBei(nu, x) is given by the imaginary part of BesselJ(nu, j * sqrt(j) * x) where j = sqrt(-1).
+ ///
+ /// the order of the the Kelvin function.
+ /// The value to compute the Kelvin function of.
+ /// The Kelvin function bei.
+ public static double KelvinBei(double nu, double x)
+ {
+ return KelvinBe(nu, x).Imaginary;
+ }
+
+ ///
+ /// Returns the Kelvin function bei.
+ /// KelvinBei(x) is given by the imaginary part of BesselJ(0, j * sqrt(j) * x) where j = sqrt(-1).
+ /// KelvinBei(x) is equivalent to KelvinBei(0, x).
+ ///
+ /// The value to compute the Kelvin function of.
+ /// The Kelvin function bei.
+ public static double KelvinBei(double x)
+ {
+ return KelvinBe(0, x).Imaginary;
+ }
+
+ ///
+ /// Returns the derivative of the Kelvin function ber.
+ ///
+ /// The order of the Kelvin function.
+ /// The value to compute the derivative of the Kelvin function of.
+ /// the derivative of the Kelvin function ber
+ public static double KelvinBerPrime(double nu, double x)
+ {
+ const double inv2Sqrt2 = 0.35355339059327376220042218105242451964241796884424; // 1/(2 * sqrt(2))
+ return inv2Sqrt2 * (-KelvinBer(nu - 1, x) + KelvinBer(nu + 1, x) - KelvinBei(nu - 1, x) + KelvinBei(nu + 1, x));
+ }
+
+ ///
+ /// Returns the derivative of the Kelvin function ber.
+ ///
+ /// The value to compute the derivative of the Kelvin function of.
+ /// The derivative of the Kelvin function ber.
+ public static double KelvinBerPrime(double x)
+ {
+ return KelvinBerPrime(0, x);
+ }
+
+ ///
+ /// Returns the derivative of the Kelvin function bei.
+ ///
+ /// The order of the Kelvin function.
+ /// The value to compute the derivative of the Kelvin function of.
+ /// the derivative of the Kelvin function bei.
+ public static double KelvinBeiPrime(double nu, double x)
+ {
+ const double inv2Sqrt2 = 0.35355339059327376220042218105242451964241796884424; // 1/(2 * sqrt(2))
+ return inv2Sqrt2 * (KelvinBer(nu - 1, x) - KelvinBer(nu + 1, x) - KelvinBei(nu - 1, x) + KelvinBei(nu + 1, x));
+ }
+
+ ///
+ /// Returns the derivative of the Kelvin function bei.
+ ///
+ /// The value to compute the derivative of the Kelvin function of.
+ /// The derivative of the Kelvin function bei.
+ public static double KelvinBeiPrime(double x)
+ {
+ return KelvinBeiPrime(0, x);
+ }
+
+ ///
+ /// Returns the Kelvin function of the second kind
+ /// KelvinKe(nu, x) is given by Exp(-nu * pi * j / 2) * BesselK(nu, x * sqrt(j)) where j = sqrt(-1).
+ /// KelvinKer(nu, x) and KelvinKei(nu, x) are the real and imaginary parts of the KelvinBe(nu, x)
+ ///
+ /// The order of the Kelvin function.
+ /// The value to calculate the kelvin function of,
+ ///
+ public static Complex KelvinKe(double nu, double x)
+ {
+ Complex PiIOver2 = new Complex(0.0, Constants.PiOver2); // pi * I / 2
+ Complex SqrtI = new Complex(Constants.Sqrt1Over2, Constants.Sqrt1Over2); // sqrt(j) = (-1)^(1/4) = (1 + j)/sqrt(2)
+ return Complex.Exp(-nu * PiIOver2) * BesselK(nu, SqrtI * x);
+ }
+
+ ///
+ /// Returns the Kelvin function ker.
+ /// KelvinKer(nu, x) is given by the real part of Exp(-nu * pi * j / 2) * BesselK(nu, sqrt(j) * x) where j = sqrt(-1).
+ ///
+ /// the order of the the Kelvin function.
+ /// The non-negative real value to compute the Kelvin function of.
+ /// The Kelvin function ker.
+ public static double KelvinKer(double nu, double x)
+ {
+ if (x <= 0.0)
+ {
+ throw new ArithmeticException();
+ }
+
+ return KelvinKe(nu, x).Real;
+ }
+
+ ///
+ /// Returns the Kelvin function ker.
+ /// KelvinKer(x) is given by the real part of Exp(-nu * pi * j / 2) * BesselK(0, sqrt(j) * x) where j = sqrt(-1).
+ /// KelvinKer(x) is equivalent to KelvinKer(0, x).
+ ///
+ /// The non-negative real value to compute the Kelvin function of.
+ /// The Kelvin function ker.
+ public static double KelvinKer(double x)
+ {
+ if (x <= 0.0)
+ {
+ throw new ArithmeticException();
+ }
+
+ return KelvinKe(0, x).Real;
+ }
+
+ ///
+ /// Returns the Kelvin function kei.
+ /// KelvinKei(nu, x) is given by the imaginary part of Exp(-nu * pi * j / 2) * BesselK(nu, sqrt(j) * x) where j = sqrt(-1).
+ ///
+ /// the order of the the Kelvin function.
+ /// The non-negative real value to compute the Kelvin function of.
+ /// The Kelvin function kei.
+ public static double KelvinKei(double nu, double x)
+ {
+ if (x <= 0.0)
+ {
+ throw new ArithmeticException();
+ }
+
+ return KelvinKe(nu, x).Imaginary;
+ }
+
+ ///
+ /// Returns the Kelvin function kei.
+ /// KelvinKei(x) is given by the imaginary part of Exp(-nu * pi * j / 2) * BesselK(0, sqrt(j) * x) where j = sqrt(-1).
+ /// KelvinKei(x) is equivalent to KelvinKei(0, x).
+ ///
+ /// The non-negative real value to compute the Kelvin function of.
+ /// The Kelvin function kei.
+ public static double KelvinKei(double x)
+ {
+ if (x <= 0.0)
+ {
+ throw new ArithmeticException();
+ }
+
+ return KelvinKe(0, x).Imaginary;
+ }
+
+ ///
+ /// Returns the derivative of the Kelvin function ker.
+ ///
+ /// The order of the Kelvin function.
+ /// The non-negative real value to compute the derivative of the Kelvin function of.
+ /// The derivative of the Kelvin function ker.
+ public static double KelvinKerPrime(double nu, double x)
+ {
+ if (x <= 0.0)
+ {
+ throw new ArithmeticException();
+ }
+
+ const double inv2Sqrt2 = 0.35355339059327376220042218105242451964241796884424; // 1/(2 * sqrt(2))
+ return inv2Sqrt2 * (-KelvinKer(nu - 1, x) + KelvinKer(nu + 1, x) - KelvinKei(nu - 1, x) + KelvinKei(nu + 1, x));
+ }
+
+ ///
+ /// Returns the derivative of the Kelvin function ker.
+ ///
+ /// The value to compute the derivative of the Kelvin function of.
+ /// The derivative of the Kelvin function ker.
+ public static double KelvinKerPrime(double x)
+ {
+ if (x <= 0.0)
+ {
+ throw new ArithmeticException();
+ }
+
+ return KelvinKerPrime(0, x);
+ }
+
+ ///
+ /// Returns the derivative of the Kelvin function kei.
+ ///
+ /// The order of the Kelvin function.
+ /// The value to compute the derivative of the Kelvin function of.
+ /// The derivative of the Kelvin function kei.
+ public static double KelvinKeiPrime(double nu, double x)
+ {
+ if (x <= 0.0)
+ {
+ throw new ArithmeticException();
+ }
+
+ const double inv2Sqrt2 = 0.35355339059327376220042218105242451964241796884424; // 1/(2 * sqrt(2))
+ return inv2Sqrt2 * (KelvinKer(nu - 1, x) - KelvinKer(nu + 1, x) - KelvinKei(nu - 1, x) + KelvinKei(nu + 1, x));
+ }
+
+ ///
+ /// Returns the derivative of the Kelvin function kei.
+ ///
+ /// The value to compute the derivative of the Kelvin function of.
+ /// The derivative of the Kelvin function kei.
+ public static double KelvinKeiPrime(double x)
+ {
+ if (x <= 0.0)
+ {
+ throw new ArithmeticException();
+ }
+
+ return KelvinKeiPrime(0, x);
+ }
+ }
+}