diff --git a/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs b/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs
index 8304cc04..e919bf09 100644
--- a/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs
+++ b/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs
@@ -17,7 +17,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
         public void BesselJ0Approx([Range(-3, 3, 0.25)] double x)
         {
             // Approx by Abramowitz/Stegun 9.4.1
-            Assert.AreEqual(Evaluate.Polynomial(x / 3.0, 1.0, 0.0, -2.2499997, 0.0, 1.2656208, 0.0, -0.3163866, 0.0, 0.0444479, 0.0, -0.0039444, 0.0, 0.0002100), SpecialFunctions.BesselJ(0, x), 1e-7);
+            Assert.AreEqual(Polynomial.Evaluate(x / 3.0, 1.0, 0.0, -2.2499997, 0.0, 1.2656208, 0.0, -0.3163866, 0.0, 0.0444479, 0.0, -0.0039444, 0.0, 0.0002100), SpecialFunctions.BesselJ(0, x), 1e-7);
         }
 
         [TestCase(0, 0.0, 0.0, 1.0000000000000000000, 0.0000000000000000000, 14)]
@@ -57,7 +57,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
         {
             // Approx by Abramowitz/Stegun 9.4.2
             Assert.AreEqual(
-                Evaluate.Polynomial(x / 3.0, 2.0 / Math.PI * Math.Log(x / 2.0) * SpecialFunctions.BesselJ(0, x) + 0.36746691, 0.0, 0.60559366, 0.0, -0.74350384, 0.0, 0.25300117, 0.0, -0.04261214, 0.0, 0.00427916, 0.0, -0.00024846),
+                Polynomial.Evaluate(x / 3.0, 2.0 / Math.PI * Math.Log(x / 2.0) * SpecialFunctions.BesselJ(0, x) + 0.36746691, 0.0, 0.60559366, 0.0, -0.74350384, 0.0, 0.25300117, 0.0, -0.04261214, 0.0, 0.00427916, 0.0, -0.00024846),
                 SpecialFunctions.BesselY(0, x), 1e-7);
         }
 
@@ -91,7 +91,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
         public void BesselI0Approx([Range(-3.75, 3.75, 0.25)] double x)
         {
             // Approx by Abramowitz/Stegun 9.8.1
-            Assert.AreEqual(Evaluate.Polynomial(x / 3.75, 1.0, 0.0, 3.5156229, 0.0, 3.0899424, 0.0, 1.2067492, 0.0, 0.2659732, 0.0, 0.0360768, 0.0, 0.0045813), SpecialFunctions.BesselI(0, x), 1e-7);
+            Assert.AreEqual(Polynomial.Evaluate(x / 3.75, 1.0, 0.0, 3.5156229, 0.0, 3.0899424, 0.0, 1.2067492, 0.0, 0.2659732, 0.0, 0.0360768, 0.0, 0.0045813), SpecialFunctions.BesselI(0, x), 1e-7);
         }
 
         [TestCase(0.0, 1.0)]
@@ -111,7 +111,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
         public void BesselI1Approx([Range(-3.75, 3.75, 0.25)] double x)
         {
             // Approx by Abramowitz/Stegun 9.8.3
-            Assert.AreEqual(Evaluate.Polynomial(x / 3.75, 0.5, 0.0, 0.87890594, 0.0, 0.51498869, 0.0, 0.15084934, 0.0, 0.02658733, 0.0, 0.00301532, 0.0, 0.00032411) * x, SpecialFunctions.BesselI(1, x), 1e-8);
+            Assert.AreEqual(Polynomial.Evaluate(x / 3.75, 0.5, 0.0, 0.87890594, 0.0, 0.51498869, 0.0, 0.15084934, 0.0, 0.02658733, 0.0, 0.00301532, 0.0, 0.00032411) * x, SpecialFunctions.BesselI(1, x), 1e-8);
         }
 
         [TestCase(0.0, 0.0)]
@@ -177,7 +177,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
         public void BesselK0Approx([Range(0.20, 2.0, 0.20)] double x)
         {
             // Approx by Abramowitz/Stegun 9.8.5
-            Assert.AreEqual(Evaluate.Polynomial(x / 2.0, -Math.Log(x / 2.0) * SpecialFunctions.BesselI(0, x) - 0.57721566, 0.0, 0.42278420, 0.0, 0.23069756, 0.0, 0.03488590, 0.0, 0.00262698, 0.0, 0.00010750, 0.0, 0.00000740), SpecialFunctions.BesselK(0, x), 1e-8);
+            Assert.AreEqual(Polynomial.Evaluate(x / 2.0, -Math.Log(x / 2.0) * SpecialFunctions.BesselI(0, x) - 0.57721566, 0.0, 0.42278420, 0.0, 0.23069756, 0.0, 0.03488590, 0.0, 0.00262698, 0.0, 0.00010750, 0.0, 0.00000740), SpecialFunctions.BesselK(0, x), 1e-8);
         }
 
         [TestCase(1e-10, 23.14178244559887)]
@@ -196,7 +196,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
         public void BesselK1Approx([Range(0.20, 2.0, 0.20)] double x)
         {
             // Approx by Abramowitz/Stegun 9.8.7
-            Assert.AreEqual(Evaluate.Polynomial(x / 2.0, x * Math.Log(x / 2.0) * SpecialFunctions.BesselI(1, x) + 1.0, 0.0, 0.15443144, 0.0, -0.67278579, 0.0, -0.18156897, 0.0, -0.01919402, 0.0, -0.00110404, 0.0, -0.00004686), SpecialFunctions.BesselK(1, x) * x, 1e-8);
+            Assert.AreEqual(Polynomial.Evaluate(x / 2.0, x * Math.Log(x / 2.0) * SpecialFunctions.BesselI(1, x) + 1.0, 0.0, 0.15443144, 0.0, -0.67278579, 0.0, -0.18156897, 0.0, -0.01919402, 0.0, -0.00110404, 0.0, -0.00004686), SpecialFunctions.BesselK(1, x) * x, 1e-8);
         }
 
         [TestCase(1e-10, 1.0e+10)]
@@ -251,7 +251,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
         public void BesselIRatioExact(int n, double zr, double zi, double cyr, double cyi, int decimalPlaces)
         {
             var z = new Complex(zr, zi);
-            var actual = SpecialFunctions.BesselI(n + 1, z, true) / SpecialFunctions.BesselI(n, z, true);
+            var actual = SpecialFunctions.BesselI(n + 1, z, SpecialFunctions.Scale.Exponential) / SpecialFunctions.BesselI(n, z, SpecialFunctions.Scale.Exponential);
             AssertHelpers.AlmostEqualRelative(new Complex(cyr, cyi), actual, decimalPlaces);
         }
 
@@ -267,7 +267,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
         public void BesselKRatioExact(int n, double zr, double zi, double cyr, double cyi, int decimalPlaces)
         {
             var z = new Complex(zr, zi);
-            var actual = SpecialFunctions.BesselK(n + 1, z, true) / SpecialFunctions.BesselK(n, z, true);
+            var actual = SpecialFunctions.BesselK(n + 1, z, SpecialFunctions.Scale.Exponential) / SpecialFunctions.BesselK(n, z, SpecialFunctions.Scale.Exponential);
             AssertHelpers.AlmostEqualRelative(new Complex(cyr, cyi), actual, decimalPlaces);
         }
 
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
+{
+    /// <summary>
+    /// Kelvin functions tests.
+    /// </summary>
+    [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/Constants.cs b/src/Numerics/Constants.cs
index 5b6fb130..59931864 100644
--- a/src/Numerics/Constants.cs
+++ b/src/Numerics/Constants.cs
@@ -96,6 +96,9 @@ namespace MathNet.Numerics
         /// <summary>The number sqrt(2pi)</summary>
         public const double Sqrt2Pi = 2.5066282746310005024157652848110452530069867406099d;
 
+        /// <summary>The number sqrt(pi/2)</summary>
+        public const double SqrtPiOver2 = 1.2533141373155002512078826424055226265034933703050d;
+
         /// <summary>The number sqrt(2*pi*e)</summary>
         public const double Sqrt2PiE = 4.1327313541224929384693918842998526494455219169913d;
 
diff --git a/src/Numerics/SpecialFunctions/Airy.cs b/src/Numerics/SpecialFunctions/Airy.cs
index 5f60ba70..37cd7f46 100644
--- a/src/Numerics/SpecialFunctions/Airy.cs
+++ b/src/Numerics/SpecialFunctions/Airy.cs
@@ -8,121 +8,195 @@ namespace MathNet.Numerics
     public static partial class SpecialFunctions
     {
         /// <summary>
-        /// Airy function Ai(z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(zta) * Ai(z), where zta = (2/3) * z * Sqrt(z).
+        /// Returns the Airy function Ai.
+        /// <para>AiryAi(z) is a solution to the Airy equation, y'' - y * z = 0.</para>
+        /// <para>AiryAi(z, Scale.Exponential) returns Exp(zta) * AiryAi(z), where zta = (2/3) * z * Sqrt(z).</para>
         /// </summary>
         /// <param name="z">The value to compute the Airy function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
-        /// <returns></returns>
-        public static Complex AiryAi(Complex z, bool expScaled = false)
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The Airy function Ai.</returns>
+        public static Complex AiryAi(Complex z, Scale scale = Scale.Unity)
         {
-            return (expScaled) ? Amos.ScaledCairy(z) : Amos.Cairy(z);
+            return (scale == Scale.Exponential) ? Amos.ScaledCairy(z) : Amos.Cairy(z);
         }
 
         /// <summary>
-        /// Airy function Ai(z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(zta) * Ai(z), where zta = (2/3) * z * Sqrt(z).
+        /// Returns the exponentially scaled Airy function Ai.
+        /// <para>ScaledAiryAi(z) is given by Exp(zta) * AiryAi(z), where zta = (2/3) * z * Sqrt(z).</para>
         /// </summary>
         /// <param name="z">The value to compute the Airy function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
-        /// <returns></returns>
-        public static double AiryAi(double z, bool expScaled = false)
+        /// <returns>The exponentially scaled Airy function Ai.</returns>
+        public static Complex ScaledAiryAi(Complex z)
         {
-            if (expScaled)
-            {
-                return Amos.ScaledCairy(z);
-            }
-            else
-            {
-                return AiryAi(new Complex(z, 0), expScaled).Real;
-            }
+            return Amos.ScaledCairy(z);
         }
 
         /// <summary>
-        /// Derivative of the Airy function Ai.
-        /// <p/>
-        /// If expScaled is true, returns Exp(zta) * d/dz Ai(z), where zta = (2/3) * z * Sqrt(z).
+        /// Returns the Airy function Ai.
+        /// <para>AiryAi(z) is a solution to the Airy equation, y'' - y * z = 0.</para>
+        /// <para>AiryAi(z, Scale.Exponential) returns Exp(zta) * AiryAi(z), where zta = (2/3) * z * Sqrt(z).</para>
+        /// </summary>
+        /// <param name="z">The value to compute the Airy function of.</param>
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The Airy function Ai.</returns>
+        public static double AiryAi(double z, Scale scale = Scale.Unity)
+        {
+            return (scale == Scale.Exponential) ? Amos.ScaledCairy(z) : AiryAi(new Complex(z, 0), scale).Real;
+        }
+
+        /// <summary>
+        /// Returns the exponentially scaled Airy function Ai.
+        /// <para>ScaledAiryAi(z) is given by Exp(zta) * AiryAi(z), where zta = (2/3) * z * Sqrt(z).</para>
+        /// </summary>
+        /// <param name="z">The value to compute the Airy function of.</param>
+        /// <returns>The exponentially scaled Airy function Ai.</returns>
+        public static double ScaledAiryAi(double z)
+        {
+            return Amos.ScaledCairy(z);
+        }
+
+        /// <summary>
+        /// Returns the derivative of the Airy function Ai.
+        /// <para>AiryAiPrime(z) is defined as d/dz AiryAi(z).</para>
+        /// <para>AiryAiPrime(z, Scale.Exponential) returns Exp(zta) * AiryAiPrime(z), where zta = (2/3) * z * Sqrt(z).</para>
+        /// </summary>
+        /// <param name="z">The value to compute the derivative of the Airy function of.</param>
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The derivative of the Airy function Ai.</returns>
+        public static Complex AiryAiPrime(Complex z, Scale scale = Scale.Unity)
+        {
+            return (scale == Scale.Exponential) ? Amos.ScaledCairyPrime(z) : Amos.CairyPrime(z);
+        }
+
+        /// <summary>
+        /// Returns the exponentially scaled derivative of Airy function Ai
+        /// <para>ScaledAiryAiPrime(z) is given by Exp(zta) * AiryAiPrime(z), where zta = (2/3) * z * Sqrt(z).</para>
+        /// </summary>
+        /// <param name="z">The value to compute the derivative of the Airy function of.</param>
+        /// <returns>The exponentially scaled derivative of Airy function Ai.</returns>
+        public static Complex ScaledAiryAiPrime(Complex z)
+        {
+            return Amos.ScaledCairyPrime(z);
+        }
+
+        /// <summary>
+        /// Returns the derivative of the Airy function Ai.
+        /// <para>AiryAiPrime(z) is defined as d/dz AiryAi(z).</para>
+        /// <para>AiryAiPrime(z, Scale.Exponential) returns Exp(zta) * AiryAiPrime(z), where zta = (2/3) * z * Sqrt(z).</para>
         /// </summary>
         /// <param name="z">The value to compute the derivative of the Airy function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
-        /// <returns></returns>
-        public static Complex AiryAiPrime(Complex z, bool expScaled = false)
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The derivative of the Airy function Ai.</returns>
+        public static double AiryAiPrime(double z, Scale scale = Scale.Unity)
         {
-            return (expScaled) ? Amos.ScaledCairyPrime(z) : Amos.CairyPrime(z);
+            return (scale == Scale.Exponential) ? Amos.ScaledCairyPrime(z) : AiryAiPrime(new Complex(z, 0), scale).Real;
         }
 
         /// <summary>
-        /// Derivative of the Airy function Ai.
-        /// <p/>
-        /// If expScaled is true, returns Exp(zta) * d/dz Ai(z), where zta = (2/3) * z * Sqrt(z).
+        /// Returns the expoenntially scaled derivative of the Airy function Ai.
+        /// <para>ScaledAiryAiPrime(z) is given by Exp(zta) * AiryAiPrime(z), where zta = (2/3) * z * Sqrt(z).</para>
         /// </summary>
         /// <param name="z">The value to compute the derivative of the Airy function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
-        /// <returns></returns>
-        public static double AiryAiPrime(double z, bool expScaled = false)
+        /// <returns>The expoenntially scaled derivative of the Airy function Ai.</returns>
+        public static double ScaledAiryAiPrime(double z)
         {
-            if (expScaled)
-            {
-                return Amos.ScaledCairyPrime(z);
-            }
-            else
-            {
-                return AiryAiPrime(new Complex(z, 0), expScaled).Real;
-            }
+            return Amos.ScaledCairyPrime(z);
         }
 
         /// <summary>
-        /// Airy function Bi(z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-axzta) * Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real).
+        /// Returns the Airy function Bi.
+        /// <para>AiryBi(z) is a solution to the Airy equation, y'' - y * z = 0.</para>
+        /// <para>AiryBi(z, Scale.Exponential) returns Exp(-Abs(zta.Real)) * AiryBi(z) where zta = (2 / 3) * z * Sqrt(z).</para>
         /// </summary>
         /// <param name="z">The value to compute the Airy function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
-        /// <returns></returns>
-        public static Complex AiryBi(Complex z, bool expScaled = false)
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The Airy function Bi.</returns>
+        public static Complex AiryBi(Complex z, Scale scale = Scale.Unity)
         {
-            return (expScaled) ? Amos.ScaledCbiry(z) : Amos.Cbiry(z);
+            return (scale == Scale.Exponential) ? Amos.ScaledCbiry(z) : Amos.Cbiry(z);
         }
 
         /// <summary>
-        /// Airy function Bi(x).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-axzta) * Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real).
+        /// Returns the exponentially scaled Airy function Bi.
+        /// <para>ScaledAiryBi(z) is given by Exp(-Abs(zta.Real)) * AiryBi(z) where zta = (2 / 3) * z * Sqrt(z).</para>
         /// </summary>
         /// <param name="z">The value to compute the Airy function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
-        /// <returns></returns>
-        public static double AiryBi(double z, bool expScaled = false)
+        /// <returns>The exponentially scaled Airy function Bi(z).</returns>
+        public static Complex ScaledAiryBi(Complex z)
+        {
+            return Amos.ScaledCbiry(z);
+        }
+
+        /// <summary>
+        /// Returns the Airy function Bi.
+        /// <para>AiryBi(z) is a solution to the Airy equation, y'' - y * z = 0.</para>
+        /// <para>AiryBi(z, Scale.Exponential) returns Exp(-Abs(zta.Real)) * AiryBi(z) where zta = (2 / 3) * z * Sqrt(z).</para>
+        /// </summary>
+        /// <param name="z">The value to compute the Airy function of.</param>
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The Airy function Bi.</returns>
+        public static double AiryBi(double z, Scale scale = Scale.Unity)
+        {
+            return AiryBi(new Complex(z, 0), scale).Real;
+        }
+
+        /// <summary>
+        /// Returns the exponentially scaled Airy function Bi.
+        /// <para>ScaledAiryBi(z) is given by Exp(-Abs(zta.Real)) * AiryBi(z) where zta = (2 / 3) * z * Sqrt(z).</para>
+        /// </summary>
+        /// <param name="z">The value to compute the Airy function of.</param>
+        /// <returns>The exponentially scaled Airy function Bi.</returns>
+        public static double ScaledAiryBi(double z)
+        {
+            return AiryBi(new Complex(z, 0), Scale.Exponential).Real;
+        }
+
+        /// <summary>
+        /// Returns the derivative of the Airy function Bi.
+        /// <para>AiryBiPrime(z) is defined as d/dz AiryBi(z).</para>
+        /// <para>AiryBiPrime(z, Scale.Exponential) returns Exp(-Abs(zta.Real)) * AiryBiPrime(z) where zta = (2 / 3) * z * Sqrt(z).</para>
+        /// </summary>
+        /// <param name="z">The value to compute the derivative of the Airy function of.</param>
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The derivative of the Airy function Bi.</returns>
+        public static Complex AiryBiPrime(Complex z, Scale scale = Scale.Unity)
+        {
+            return (scale == Scale.Exponential) ? Amos.ScaledCbiryPrime(z) : Amos.CbiryPrime(z);
+        }
+
+        /// <summary>
+        /// Returns the exponentially scaled derivative of the Airy function Bi.
+        /// <para>ScaledAiryBiPrime(z) is given by Exp(-Abs(zta.Real)) * AiryBiPrime(z) where zta = (2 / 3) * z * Sqrt(z).</para>
+        /// </summary>
+        /// <param name="z">The value to compute the derivative of the Airy function of.</param>
+        /// <returns>The exponentially scaled derivative of the Airy function Bi.</returns>
+        public static Complex ScaledAiryBiPrime(Complex z)
         {
-            return AiryBi(new Complex(z, 0), expScaled).Real;
+            return Amos.ScaledCbiryPrime(z);
         }
 
         /// <summary>
-        /// Derivative of the Airy function Bi(z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-axzta) * d/dz Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real).
+        /// Returns the derivative of the Airy function Bi.
+        /// <para>AiryBiPrime(z) is defined as d/dz AiryBi(z).</para>
+        /// <para>AiryBiPrime(z, Scale.Exponential) returns Exp(-Abs(zta.Real)) * AiryBiPrime(z) where zta = (2 / 3) * z * Sqrt(z).</para>
         /// </summary>
         /// <param name="z">The value to compute the derivative of the Airy function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
-        /// <returns></returns>
-        public static Complex AiryBiPrime(Complex z, bool expScaled = false)
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The derivative of the Airy function Bi.</returns>
+        public static double AiryBiPrime(double z, Scale scale = Scale.Unity)
         {
-            return (expScaled) ? Amos.ScaledCbiryPrime(z) : Amos.CbiryPrime(z);
+            return AiryBiPrime(new Complex(z, 0), scale).Real;
         }
 
         /// <summary>
-        /// Derivative of the Airy function Bi(z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-axzta) * d/dz Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real).
+        /// Returns the exponentially scaled derivative of the Airy function Bi.
+        /// <para>ScaledAiryBiPrime(z) is given by Exp(-Abs(zta.Real)) * AiryBiPrime(z) where zta = (2 / 3) * z * Sqrt(z).</para>
         /// </summary>
         /// <param name="z">The value to compute the derivative of the Airy function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
-        /// <returns></returns>
-        public static double AiryBiPrime(double z, bool expScaled = false)
+        /// <returns>The exponentially scaled derivative of the Airy function Bi.</returns>
+        public static double ScaledAiryBiPrime(double z)
         {
-            return AiryBiPrime(new Complex(z, 0), expScaled).Real;
+            return AiryBiPrime(new Complex(z, 0), Scale.Exponential).Real;
         }
     }
 }
diff --git a/src/Numerics/SpecialFunctions/Bessel.cs b/src/Numerics/SpecialFunctions/Bessel.cs
index 83c3e85b..c80eaed5 100644
--- a/src/Numerics/SpecialFunctions/Bessel.cs
+++ b/src/Numerics/SpecialFunctions/Bessel.cs
@@ -8,122 +8,211 @@ namespace MathNet.Numerics
     public static partial class SpecialFunctions
     {
         /// <summary>
-        /// Bessel function of the first kind, J(v, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-Abs(y)) * J(v, z) where y = z.Imaginary.
+        /// Returns the Bessel function of the first kind.
+        /// <para>BesselJ(n, z) is a solution to the Bessel differential equation.</para>
+        /// <para>BesselJ(n, z, Scale.Exponential) returns Exp(-Abs(z.Imaginary)) * BesselJ(n, z).</para>
         /// </summary>
-        /// <param name="v">The order of the Bessel function</param>
+        /// <param name="n">The order of the Bessel function.</param>
         /// <param name="z">The value to compute the Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
-        /// <returns></returns>
-        public static Complex BesselJ(double v, Complex z, bool expScaled = false)
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The Bessel function of the first kind.</returns>
+        public static Complex BesselJ(double n, Complex z, Scale scale = Scale.Unity)
         {
-            return (expScaled) ? Amos.ScaledCbesj(v, z) : Amos.Cbesj(v, z);
+            return (scale == Scale.Exponential) ? Amos.ScaledCbesj(n, z) : Amos.Cbesj(n, z);
         }
 
         /// <summary>
-        /// Bessel function of the first kind, J(v, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-Abs(y)) * J(v, z) where y = z.Imaginary.
+        /// Returns the exponentially scaled Bessel function of the first kind.
+        /// <para>ScaledBesselJ(n, z) is given by Exp(-Abs(z.Imaginary)) * BesselJ(n, z).</para>
         /// </summary>
-        /// <param name="v">The order of the Bessel function</param>
+        /// <param name="n">The order of the Bessel function.</param>
         /// <param name="z">The value to compute the Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
-        /// <returns></returns>
-        public static double BesselJ(double v, double z, bool expScaled = false)
+        /// <returns>The exponentially scaled Bessel function of the first kind.</returns>
+        public static Complex ScaledBesselJ(double n, Complex z)
         {
-            return (expScaled) ? Amos.ScaledCbesj(v, z) : Amos.Cbesj(v, z);
+            return Amos.ScaledCbesj(n, z);
         }
 
         /// <summary>
-        /// Bessel function of the second kind, Y(v, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-Abs(y)) * Y(v, z) where y = z.Imaginary.
+        /// Returns the Bessel function of the first kind.
+        /// <para>BesselJ(n, z) is a solution to the Bessel differential equation.</para>
+        /// <para>BesselJ(n, z, Scale.Exponential) returns Exp(-Abs(z.Imaginary)) * J(n, z).</para>
         /// </summary>
-        /// <param name="v">The order of the Bessel function</param>
+        /// <param name="n">The order of the Bessel function.</param>
         /// <param name="z">The value to compute the Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
-        /// <returns></returns>
-        public static Complex BesselY(double v, Complex z, bool expScaled = false)
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The Bessel function of the first kind.</returns>
+        public static double BesselJ(double n, double z, Scale scale = Scale.Unity)
         {
-            return (expScaled) ? Amos.ScaledCbesy(v, z) : Amos.Cbesy(v, z);
+            return (scale == Scale.Exponential) ? Amos.ScaledCbesj(n, z) : Amos.Cbesj(n, z);
         }
 
         /// <summary>
-        /// Bessel function of the second kind, Y(v, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-Abs(y)) * Y(v, z) where y = z.Imaginary.
+        /// Returns the exponentially scaled Bessel function of the first kind.
+        /// <para>ScaledBesselJ(n, z) is given by Exp(-Abs(z.Imaginary)) * BesselJ(n, z).</para>
         /// </summary>
-        /// <param name="v">The order of the Bessel function</param>
+        /// <param name="n">The order of the Bessel function.</param>
         /// <param name="z">The value to compute the Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
-        /// <returns></returns>
-        public static double BesselY(double v, double z, bool expScaled = false)
+        /// <returns>The exponentially scaled Bessel function of the first kind.</returns>
+        public static double ScaledBesselJ(double n, double z)
         {
-            return (expScaled) ? Amos.ScaledCbesy(v, z) : Amos.Cbesy(v, z);
+            return Amos.ScaledCbesj(n, z);
         }
 
         /// <summary>
-        /// Modified Bessel function of the first kind, I(v, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-Abs(x)) * I(v, z) where x = z.Real.
+        /// Returns the Bessel function of the second kind.
+        /// <para>BesselY(n, z) is a solution to the Bessel differential equation.</para>
+        /// <para>BesselY(n, z, Scale.Exponential) returns Exp(-Abs(z.Imaginary)) * BesselY(n, z).</para>
         /// </summary>
-        /// <param name="v">The order of the Bessel function</param>
+        /// <param name="n">The order of the Bessel function.</param>
         /// <param name="z">The value to compute the Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
-        /// <returns></returns>
-        public static Complex BesselI(double v, Complex z, bool expScaled = false)
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The Bessel function of the second kind.</returns>
+        public static Complex BesselY(double n, Complex z, Scale scale = Scale.Unity)
         {
-            return (expScaled) ? Amos.ScaledCbesi(v, z) : Amos.Cbesi(v, z);
+            return (scale == Scale.Exponential) ? Amos.ScaledCbesy(n, z) : Amos.Cbesy(n, z);
         }
 
         /// <summary>
-        /// Modified Bessel function of the first kind, I(v, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-Abs(x)) * I(v, z) where x = z.Real.
+        /// Returns the exponentially scaled Bessel function of the second kind.
+        /// <para>ScaledBesselY(n, z) is given by Exp(-Abs(z.Imaginary)) * Y(n, z).</para>
         /// </summary>
-        /// <param name="v">The order of the Bessel function</param>
+        /// <param name="n">The order of the Bessel function.</param>
         /// <param name="z">The value to compute the Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
-        /// <returns></returns>
-        public static double BesselI(double v, double z, bool expScaled = false)
+        /// <returns>The exponentially scaled Bessel function of the second kind.</returns>
+        public static Complex ScaledBesselY(double n, Complex z)
         {
-            if (expScaled)
-            {
-                return Amos.ScaledCbesi(v, z);
-            }
-            else
-            {
-                return BesselI(v, new Complex(z, 0), expScaled).Real;
-            }
+            return Amos.ScaledCbesy(n, z);
         }
 
         /// <summary>
-        /// Modified Bessel function of the second kind, K(v, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(z) * K(v, z).
+        /// Returns the Bessel function of the second kind.
+        /// <para>BesselY(n, z) is a solution to the Bessel differential equation.</para>
+        /// <para>BesselY(n, z, Scale.Exponential) returns Exp(-Abs(z.Imaginary)) * BesselY(n, z).</para>
         /// </summary>
-        /// <param name="v">The order of the Bessel function</param>
+        /// <param name="n">The order of the Bessel function.</param>
         /// <param name="z">The value to compute the Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
-        /// <returns></returns>
-        public static Complex BesselK(double v, Complex z, bool expScaled = false)
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The Bessel function of the second kind.</returns>
+        public static double BesselY(double n, double z, Scale scale = Scale.Unity)
         {
-            return (expScaled) ? Amos.ScaledCbesk(v, z) : Amos.Cbesk(v, z);
+            return (scale == Scale.Exponential) ? Amos.ScaledCbesy(n, z) : Amos.Cbesy(n, z);
         }
 
         /// <summary>
-        /// Modified Bessel function of the second kind, K(v, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(z) * K(v, z).
+        /// Returns the exponentially scaled Bessel function of the second kind.
+        /// <para>ScaledBesselY(n, z) is given by Exp(-Abs(z.Imaginary)) * BesselY(n, z).</para>
         /// </summary>
-        /// <param name="v">The order of the Bessel function</param>
+        /// <param name="n">The order of the Bessel function.</param>
         /// <param name="z">The value to compute the Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
-        /// <returns></returns>
-        public static double BesselK(double v, double z, bool expScaled = false)
+        /// <returns>The exponentially scaled Bessel function of the second kind.</returns>
+        public static double ScaledBesselY(double n, double z)
         {
-            return (expScaled) ? Amos.ScaledCbesk(v, z) : Amos.Cbesk(v, z);
+            return Amos.ScaledCbesy(n, z);
+        }
+
+        /// <summary>
+        /// Returns the modified Bessel function of the first kind.
+        /// <para>BesselI(n, z) is a solution to the modified Bessel differential equation.</para>
+        /// <para>BesselI(n, z, Scale.Exponential) returns Exp(-Abs(z.Real)) * BesselI(n, z).</para>
+        /// </summary>
+        /// <param name="n">The order of the modified Bessel function.</param>
+        /// <param name="z">The value to compute the modified Bessel function of.</param>
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The modified Bessel function of the first kind.</returns>
+        public static Complex BesselI(double n, Complex z, Scale scale = Scale.Unity)
+        {
+            return (scale == Scale.Exponential) ? Amos.ScaledCbesi(n, z) : Amos.Cbesi(n, z);
+        }
+
+        /// <summary>
+        /// Returns the exponentially scaled modified Bessel function of the first kind.
+        /// <para>ScaledBesselI(n, z) is given by Exp(-Abs(z.Real)) * BesselI(n, z).</para>
+        /// </summary>
+        /// <param name="n">The order of the modified Bessel function.</param>
+        /// <param name="z">The value to compute the modified Bessel function of.</param>
+        /// <returns>The exponentially scaled modified Bessel function of the first kind.</returns>
+        public static Complex ScaledBesselI(double n, Complex z)
+        {
+            return Amos.ScaledCbesi(n, z);
+        }
+
+        /// <summary>
+        /// Returns the modified Bessel function of the first kind.
+        /// <para>BesselI(n, z) is a solution to the modified Bessel differential equation.</para>
+        /// <para>BesselI(n, z, Scale.Exponential) returns Exp(-Abs(z.Real)) * BesselI(n, z).</para>
+        /// </summary>
+        /// <param name="n">The order of the modified Bessel function.</param>
+        /// <param name="z">The value to compute the modified Bessel function of.</param>
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The modified Bessel function of the first kind.</returns>
+        public static double BesselI(double n, double z, Scale scale = Scale.Unity)
+        {
+            return (scale == Scale.Exponential) ? Amos.ScaledCbesi(n, z) : BesselI(n, new Complex(z, 0), scale).Real;
+        }
+
+        /// <summary>
+        /// Returns the exponentially scaled modified Bessel function of the first kind.
+        /// <para>ScaledBesselI(n, z) is given by Exp(-Abs(z.Real)) * BesselI(n, z).</para>
+        /// </summary>
+        /// <param name="n">The order of the modified Bessel function.</param>
+        /// <param name="z">The value to compute the modified Bessel function of.</param>
+        /// <returns>The exponentially scaled modified Bessel function of the first kind.</returns>
+        public static double ScaledBesselI(double n, double z)
+        {
+            return Amos.ScaledCbesi(n, z);
+        }
+
+        /// <summary>
+        /// Returns the modified Bessel function of the second kind.
+        /// <para>BesselK(n, z) is a solution to the modified Bessel differential equation.</para>
+        /// <para>BesselK(n, z, Scale.Exponential) returns Exp(z) * BesselK(n, z).</para>
+        /// </summary>
+        /// <param name="n">The order of the modified Bessel function.</param>
+        /// <param name="z">The value to compute the modified Bessel function of.</param>
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The modified Bessel function of the second kind.</returns>
+        public static Complex BesselK(double n, Complex z, Scale scale = Scale.Unity)
+        {
+            return (scale == Scale.Exponential) ? Amos.ScaledCbesk(n, z) : Amos.Cbesk(n, z);
+        }
+
+        /// <summary>
+        /// Returns the exponentially scaled modified Bessel function of the second kind.
+        /// <para>ScaledBesselK(n, z) is given by Exp(z) * BesselK(n, z).</para>
+        /// </summary>
+        /// <param name="n">The order of the modified Bessel function.</param>
+        /// <param name="z">The value to compute the modified Bessel function of.</param>
+        /// <returns>The exponentially scaled modified Bessel function of the second kind.</returns>
+        public static Complex ScaledBesselK(double n, Complex z)
+        {
+            return Amos.ScaledCbesk(n, z);
+        }
+
+        /// <summary>
+        /// Returns the modified Bessel function of the second kind.
+        /// <para>BesselK(n, z) is a solution to the modified Bessel differential equation.</para>
+        /// <para>BesselK(n, z, Scale.Exponential) returns Exp(z) * BesselK(n, z).</para>
+        /// </summary>
+        /// <param name="n">The order of the modified Bessel function.</param>
+        /// <param name="z">The value to compute the modified Bessel function of.</param>
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The modified Bessel function of the second kind.</returns>
+        public static double BesselK(double n, double z, Scale scale = Scale.Unity)
+        {
+            return (scale == Scale.Exponential) ? Amos.ScaledCbesk(n, z) : Amos.Cbesk(n, z);
+        }
+
+        /// <summary>
+        /// Returns the exponentially scaled modified Bessel function of the second kind.
+        /// <para>ScaledBesselK(n, z) is given by Exp(z) * BesselK(n, z).</para>
+        /// </summary>
+        /// <param name="n">The order of the modified Bessel function.</param>
+        /// <param name="z">The value to compute the modified Bessel function of.</param>
+        /// <returns>The exponentially scaled modified Bessel function of the second kind.</returns>
+        public static double ScaledBesselK(double n, double z)
+        {
+            return Amos.ScaledCbesk(n, z);
         }
     }
 }
diff --git a/src/Numerics/SpecialFunctions/Hankel.cs b/src/Numerics/SpecialFunctions/Hankel.cs
index 32943511..1f820c33 100644
--- a/src/Numerics/SpecialFunctions/Hankel.cs
+++ b/src/Numerics/SpecialFunctions/Hankel.cs
@@ -8,59 +8,55 @@ namespace MathNet.Numerics
     public static partial class SpecialFunctions
     {
         /// <summary>
-        /// Hankel function of the first kind, H1(n, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-z * j) * H1(n, z) where j = Sqrt(-1).
+        /// Returns the Hankel function of the first kind.
+        /// <para>HankelH1(n, z) is defined as BesselJ(n, z) + j * BesselY(n, z).</para>
+        /// <para>HankelH1(n, z, Scale.Exponential) returns Exp(-z * j) * HankelH1(n, z) where j = Sqrt(-1).</para>
         /// </summary>
-        /// <param name="n">The order of the Bessel function</param>
-        /// <param name="z">The value to compute the Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Hankel function</param>
-        /// <returns></returns>
-        public static Complex HankelH1(double n, Complex z, bool expScaled = false)
+        /// <param name="n">The order of the Hankel function.</param>
+        /// <param name="z">The value to compute the Hankel function of.</param>
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The Hankel function of the first kind.</returns>
+        public static Complex HankelH1(double n, Complex z, Scale scale = Scale.Unity)
         {
-            return (expScaled) ? Amos.ScaledCbesh1(n, z) : Amos.Cbesh1(n, z);
+            return (scale ==  Scale.Exponential) ? Amos.ScaledCbesh1(n, z) : Amos.Cbesh1(n, z);
         }
 
         /// <summary>
-        /// Hankel function of the first kind, H1(n, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-z * j) * H1(n, z) where j = Sqrt(-1).
+        /// Returns the exponentially scaled Hankel function of the first kind.
+        /// <para>ScaledHankelH1(n, z) is given by Exp(-z * j) * HankelH1(n, z) where j = Sqrt(-1).</para>
         /// </summary>
-        /// <param name="n">The order of the Bessel function</param>
-        /// <param name="z">The value to compute the Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Hankel function</param>
-        /// <returns></returns>
-        public static double HankelH1(double n, double z, bool expScaled = false)
+        /// <param name="n">The order of the Hankel function.</param>
+        /// <param name="z">The value to compute the Hankel function of.</param>
+        /// <returns>The exponentially scaled Hankel function of the first kind.</returns>
+        public static Complex ScaledHankelH1(double n, Complex z)
         {
-            return HankelH1(n, new Complex(z, 0), expScaled).Real;
+            return Amos.ScaledCbesh1(n, z);
         }
 
         /// <summary>
-        /// Hankel function of the second kind, H2(n, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(z * j) * H2(n, z) where j = Sqrt(-1).
+        /// Returns the Hankel function of the second kind.
+        /// <para>HankelH2(n, z) is defined as BesselJ(n, z) - j * BesselY(n, z).</para>
+        /// <para>HankelH2(n, z, Scale.Exponential) returns Exp(z * j) * HankelH2(n, z) where j = Sqrt(-1).</para>
         /// </summary>
-        /// <param name="n">The order of the Hankel function</param>
-        /// <param name="z">The value to compute the Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Hankel function</param>
-        /// <returns></returns>
-        public static Complex HankelH2(double n, Complex z, bool expScaled = false)
+        /// <param name="n">The order of the Hankel function.</param>
+        /// <param name="z">The value to compute the Hankel function of.</param>
+        /// <param name="scale">The option to set the scaling factor.</param>
+        /// <returns>The Hankel function of the second kind.</returns>
+        public static Complex HankelH2(double n, Complex z, Scale scale = Scale.Unity)
         {
-            return (expScaled) ? Amos.ScaledCbesh2(n, z) : Amos.Cbesh2(n, z);
+            return (scale == Scale.Exponential) ? Amos.ScaledCbesh2(n, z) : Amos.Cbesh2(n, z);
         }
 
         /// <summary>
-        /// Hankel function of the second kind, H2(n, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(z * j) * H2(n, z) where j = Sqrt(-1).
+        /// Returns the exponentially scaled Hankel function of the second kind.
+        /// <para>ScaledHankelH2(n, z) is given by Exp(z * j) * HankelH2(n, z) where j = Sqrt(-1).</para>
         /// </summary>
-        /// <param name="n">The order of the Bessel function</param>
-        /// <param name="z">The value to compute the Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled Hankel function</param>
-        /// <returns></returns>
-        public static double HankelH2(double n, double z, bool expScaled = false)
+        /// <param name="n">The order of the Hankel function.</param>
+        /// <param name="z">The value to compute the Hankel function of.</param>
+        /// <returns>The exponentially scaled Hankel function of the second kind.</returns>
+        public static Complex ScaledHankelH2(double n, Complex z)
         {
-            return HankelH2(n, new Complex(z, 0), expScaled).Real;
+            return Amos.ScaledCbesh2(n, z);
         }
     }
 }
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
+{
+    /// <summary>
+    /// This partial implementation of the SpecialFunctions class contains all methods related to the modified Bessel function.
+    /// </summary>
+    public static partial class SpecialFunctions
+    {
+        /// <summary>
+        /// Returns the Kelvin function of the first kind.
+        /// <para>KelvinBe(nu, x) is given by BesselJ(0, j * sqrt(j) * x) where j = sqrt(-1).</para>
+        /// <para>KelvinBer(nu, x) and KelvinBei(nu, x) are the real and imaginary parts of the KelvinBe(nu, x)</para>
+        /// </summary>
+        /// <param name="nu">the order of the the Kelvin function.</param>
+        /// <param name="x">The value to compute the Kelvin function of.</param>
+        /// <returns>The Kelvin function of the first kind.</returns>
+        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);
+        }
+
+        /// <summary>
+        /// Returns the Kelvin function ber.
+        /// <para>KelvinBer(nu, x) is given by the real part of BesselJ(nu, j * sqrt(j) * x) where j = sqrt(-1).</para>
+        /// </summary>
+        /// <param name="nu">the order of the the Kelvin function.</param>
+        /// <param name="x">The value to compute the Kelvin function of.</param>
+        /// <returns>The Kelvin function ber.</returns>
+        public static double KelvinBer(double nu, double x)
+        {
+            return KelvinBe(nu, x).Real;
+        }
+
+        /// <summary>
+        /// Returns the Kelvin function ber.
+        /// <para>KelvinBer(x) is given by the real part of BesselJ(0, j * sqrt(j) * x) where j = sqrt(-1).</para>
+        /// <para>KelvinBer(x) is equivalent to KelvinBer(0, x).</para>
+        /// </summary>
+        /// <param name="x">The value to compute the Kelvin function of.</param>
+        /// <returns>The Kelvin function ber.</returns>
+        public static double KelvinBer(double x)
+        {
+            return KelvinBe(0, x).Real;
+        }
+
+        /// <summary>
+        /// Returns the Kelvin function bei.
+        /// <para>KelvinBei(nu, x) is given by the imaginary part of BesselJ(nu, j * sqrt(j) * x) where j = sqrt(-1).</para>
+        /// </summary>
+        /// <param name="nu">the order of the the Kelvin function.</param>
+        /// <param name="x">The value to compute the Kelvin function of.</param>
+        /// <returns>The Kelvin function bei.</returns>
+        public static double KelvinBei(double nu, double x)
+        {
+            return KelvinBe(nu, x).Imaginary;
+        }
+
+        /// <summary>
+        /// Returns the Kelvin function bei.
+        /// <para>KelvinBei(x) is given by the imaginary part of BesselJ(0, j * sqrt(j) * x) where j = sqrt(-1).</para>
+        /// <para>KelvinBei(x) is equivalent to KelvinBei(0, x).</para>
+        /// </summary>
+        /// <param name="x">The value to compute the Kelvin function of.</param>
+        /// <returns>The Kelvin function bei.</returns>
+        public static double KelvinBei(double x)
+        {
+            return KelvinBe(0, x).Imaginary;
+        }
+
+        /// <summary>
+        /// Returns the derivative of the Kelvin function ber.
+        /// </summary>
+        /// <param name="nu">The order of the Kelvin function.</param>
+        /// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
+        /// <returns>the derivative of the Kelvin function ber</returns>
+        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));
+        }
+
+        /// <summary>
+        /// Returns the derivative of the Kelvin function ber.
+        /// </summary>
+        /// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
+        /// <returns>The derivative of the Kelvin function ber.</returns>
+        public static double KelvinBerPrime(double x)
+        {
+            return KelvinBerPrime(0, x);
+        }
+
+        /// <summary>
+        /// Returns the derivative of the Kelvin function bei.
+        /// </summary>
+        /// <param name="nu">The order of the Kelvin function.</param>
+        /// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
+        /// <returns>the derivative of the Kelvin function bei.</returns>
+        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));
+        }
+
+        /// <summary>
+        /// Returns the derivative of the Kelvin function bei.
+        /// </summary>
+        /// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
+        /// <returns>The derivative of the Kelvin function bei.</returns>
+        public static double KelvinBeiPrime(double x)
+        {
+            return KelvinBeiPrime(0, x);
+        }
+
+        /// <summary>
+        /// Returns the Kelvin function of the second kind
+        /// <para>KelvinKe(nu, x) is given by Exp(-nu * pi * j / 2) * BesselK(nu, x * sqrt(j)) where j = sqrt(-1).</para>
+        /// <para>KelvinKer(nu, x) and KelvinKei(nu, x) are the real and imaginary parts of the KelvinBe(nu, x)</para>
+        /// </summary>
+        /// <param name="nu">The order of the Kelvin function.</param>
+        /// <param name="x">The value to calculate the kelvin function of,</param>
+        /// <returns></returns>
+        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);
+        }
+
+        /// <summary>
+        /// Returns the Kelvin function ker.
+        /// <para>KelvinKer(nu, x) is given by the real part of Exp(-nu * pi * j / 2) * BesselK(nu, sqrt(j) * x) where j = sqrt(-1).</para>
+        /// </summary>
+        /// <param name="nu">the order of the the Kelvin function.</param>
+        /// <param name="x">The non-negative real value to compute the Kelvin function of.</param>
+        /// <returns>The Kelvin function ker.</returns>
+        public static double KelvinKer(double nu, double x)
+        {
+            if (x <= 0.0)
+            {
+                throw new ArithmeticException();
+            }
+
+            return KelvinKe(nu, x).Real;
+        }
+
+        /// <summary>
+        /// Returns the Kelvin function ker.
+        /// <para>KelvinKer(x) is given by the real part of Exp(-nu * pi * j / 2) * BesselK(0, sqrt(j) * x) where j = sqrt(-1).</para>
+        /// <para>KelvinKer(x) is equivalent to KelvinKer(0, x).</para>
+        /// </summary>
+        /// <param name="x">The non-negative real value to compute the Kelvin function of.</param>
+        /// <returns>The Kelvin function ker.</returns>
+        public static double KelvinKer(double x)
+        {
+            if (x <= 0.0)
+            {
+                throw new ArithmeticException();
+            }
+
+            return KelvinKe(0, x).Real;
+        }
+
+        /// <summary>
+        /// Returns the Kelvin function kei.
+        /// <para>KelvinKei(nu, x) is given by the imaginary part of Exp(-nu * pi * j / 2) * BesselK(nu, sqrt(j) * x) where j = sqrt(-1).</para>
+        /// </summary>
+        /// <param name="nu">the order of the the Kelvin function.</param>
+        /// <param name="x">The non-negative real value to compute the Kelvin function of.</param>
+        /// <returns>The Kelvin function kei.</returns>
+        public static double KelvinKei(double nu, double x)
+        {
+            if (x <= 0.0)
+            {
+                throw new ArithmeticException();
+            }
+
+            return KelvinKe(nu, x).Imaginary;
+        }
+
+        /// <summary>
+        /// Returns the Kelvin function kei.
+        /// <para>KelvinKei(x) is given by the imaginary part of Exp(-nu * pi * j / 2) * BesselK(0, sqrt(j) * x) where j = sqrt(-1).</para>
+        /// <para>KelvinKei(x) is equivalent to KelvinKei(0, x).</para>
+        /// </summary>
+        /// <param name="x">The non-negative real value to compute the Kelvin function of.</param>
+        /// <returns>The Kelvin function kei.</returns>
+        public static double KelvinKei(double x)
+        {
+            if (x <= 0.0)
+            {
+                throw new ArithmeticException();
+            }
+
+            return KelvinKe(0, x).Imaginary;
+        }
+
+        /// <summary>
+        /// Returns the derivative of the Kelvin function ker.
+        /// </summary>
+        /// <param name="nu">The order of the Kelvin function.</param>
+        /// <param name="x">The non-negative real value to compute the derivative of the Kelvin function of.</param>
+        /// <returns>The derivative of the Kelvin function ker.</returns>
+        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));
+        }
+
+        /// <summary>
+        /// Returns the derivative of the Kelvin function ker.
+        /// </summary>
+        /// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
+        /// <returns>The derivative of the Kelvin function ker.</returns>
+        public static double KelvinKerPrime(double x)
+        {
+            if (x <= 0.0)
+            {
+                throw new ArithmeticException();
+            }
+
+            return KelvinKerPrime(0, x);
+        }
+
+        /// <summary>
+        /// Returns the derivative of the Kelvin function kei.
+        /// </summary>
+        /// <param name="nu">The order of the Kelvin function.</param>
+        /// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
+        /// <returns>The derivative of the Kelvin function kei.</returns>
+        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));
+        }
+
+        /// <summary>
+        /// Returns the derivative of the Kelvin function kei.
+        /// </summary>
+        /// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
+        /// <returns>The derivative of the Kelvin function kei.</returns>
+        public static double KelvinKeiPrime(double x)
+        {
+            if (x <= 0.0)
+            {
+                throw new ArithmeticException();
+            }
+
+            return KelvinKeiPrime(0, x);
+        }
+    }
+}
diff --git a/src/Numerics/SpecialFunctions/Options.cs b/src/Numerics/SpecialFunctions/Options.cs
new file mode 100644
index 00000000..6342f0c7
--- /dev/null
+++ b/src/Numerics/SpecialFunctions/Options.cs
@@ -0,0 +1,23 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+namespace MathNet.Numerics
+{
+    public static partial class SpecialFunctions
+    {
+        public enum Scale
+        {
+            /// <summary>
+            /// For Bessel-related functions, no scaling factor is applied.
+            /// </summary>
+            Unity = 0,
+
+            /// <summary>
+            /// For Bessel-related functions, exponential scaling is applied.
+            /// </summary>
+            Exponential = 1
+        }
+    }
+}
diff --git a/src/Numerics/SpecialFunctions/SphericalBessel.cs b/src/Numerics/SpecialFunctions/SphericalBessel.cs
index 5bf43db1..160cfa5c 100644
--- a/src/Numerics/SpecialFunctions/SphericalBessel.cs
+++ b/src/Numerics/SpecialFunctions/SphericalBessel.cs
@@ -9,67 +9,121 @@ namespace MathNet.Numerics
     public static partial class SpecialFunctions
     {
         /// <summary>
-        /// Spherical Bessel function of the first kind, j(v, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-Abs(y)) * j(v, z) where y = z.Imaginary.
+        /// Returns the spherical Bessel function of the first kind.
+        /// <para>SphericalBesselJ(n, z) is given by Sqrt(pi/2) / Sqrt(z) * BesselJ(n + 1/2, z).</para>
         /// </summary>
-        /// <param name="v">The order of the spherical Bessel function</param>
+        /// <param name="n">The order of the spherical Bessel function.</param>
         /// <param name="z">The value to compute the spherical Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled spherical Bessel function</param>
-        /// <returns></returns>
-        public static Complex SphericalBesselJ(double v, Complex z, bool expScaled = false)
+        /// <returns>The spherical Bessel function of the first kind.</returns>
+        public static Complex SphericalBesselJ(double n, Complex z)
         {
-            const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
+            if (double.IsNaN(n) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
+            {
+                return new Complex(double.NaN, double.NaN);
+            }
+            
+            if (double.IsInfinity(z.Real))
+            {
+                return (z.Imaginary == 0) ? Complex.Zero : new Complex(double.PositiveInfinity, double.PositiveInfinity);
+            }
 
-            return rthpi * BesselJ(v + 0.5, z, expScaled) / Complex.Sqrt(z);
+            if (z.Real == 0 && z.Imaginary == 0)
+            {
+                return (n == 0) ? 1 : 0;
+            }
+
+            return Constants.SqrtPiOver2 * BesselJ(n + 0.5, z, Scale.Unity) / Complex.Sqrt(z);
         }
 
         /// <summary>
-        /// Spherical Bessel function of the first kind, j(v, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-Abs(y)) * j(v, z) where y = z.Imaginary.
+        /// Returns the spherical Bessel function of the first kind.
+        /// <para>SphericalBesselJ(n, z) is given by Sqrt(pi/2) / Sqrt(z) * BesselJ(n + 1/2, z).</para>
         /// </summary>
-        /// <param name="v">The order of the spherical Bessel function</param>
+        /// <param name="n">The order of the spherical Bessel function.</param>
         /// <param name="z">The value to compute the spherical Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled spherical Bessel function</param>
-        /// <returns></returns>
-        public static double SphericalBesselJ(double v, double z, bool expScaled = false)
+        /// <returns>The spherical Bessel function of the first kind.</returns>
+        public static double SphericalBesselJ(double n, double z)
         {
-            const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
+            if (double.IsNaN(n) || double.IsNaN(z))
+            {
+                return double.NaN;
+            }
+
+            if (n < 0)
+            {
+                return double.NaN;
+            }
+
+            if (double.IsInfinity(z))
+            {
+                return 0;
+            }
 
-            return rthpi * BesselJ(v + 0.5, z, expScaled) / Math.Sqrt(z);
+            if (z == 0)
+            {
+                return (n == 0) ? 1 : 0;
+            }
+
+            return Constants.SqrtPiOver2 * BesselJ(n + 0.5, z, Scale.Unity) / Math.Sqrt(z);
         }
 
         /// <summary>
-        /// Spherical Bessel function of the second kind, y(v, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-Abs(y)) * y(v, z) where y = z.Imaginary.
+        /// Returns the spherical Bessel function of the second kind.
+        /// <para>SphericalBesselY(n, z) is given by Sqrt(pi/2) / Sqrt(z) * BesselY(n + 1/2, z).</para>
         /// </summary>
-        /// <param name="v">The order of the spherical Bessel function</param>
+        /// <param name="n">The order of the spherical Bessel function.</param>
         /// <param name="z">The value to compute the spherical Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled spherical Bessel function</param>
-        /// <returns></returns>
-        public static Complex SphericalBesselY(double v, Complex z, bool expScaled = false)
+        /// <returns>The spherical Bessel function of the second kind.</returns>
+        public static Complex SphericalBesselY(double n, Complex z)
         {
-            const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
+            if (double.IsNaN(n) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
+            {
+                return new Complex(double.NaN, double.NaN);
+            }
+
+            if (double.IsInfinity(z.Real))
+            {
+                return (z.Imaginary == 0) ? Complex.Zero : new Complex(double.PositiveInfinity, double.PositiveInfinity);
+            }
 
-            return rthpi * BesselY(v + 0.5, z, expScaled) / Complex.Sqrt(z);
+            if (z.Real == 0 && z.Imaginary == 0)
+            {
+                return new Complex(double.NaN, double.NaN);
+            }
+
+            return Constants.SqrtPiOver2 * BesselY(n + 0.5, z, Scale.Unity) / Complex.Sqrt(z);
         }
 
         /// <summary>
-        /// Spherical Bessel function of the second kind, y(v, z).
-        /// <p/>
-        /// If expScaled is true, returns Exp(-Abs(y)) * y(v, z) where y = z.Imaginary.
+        /// Returns the spherical Bessel function of the second kind.
+        /// <para>SphericalBesselY(n, z) is given by Sqrt(pi/2) / Sqrt(z) * BesselY(n + 1/2, z).</para>
         /// </summary>
-        /// <param name="v">The order of the spherical Bessel function</param>
+        /// <param name="n">The order of the spherical Bessel function.</param>
         /// <param name="z">The value to compute the spherical Bessel function of.</param>
-        /// <param name="expScaled">If true, returns exponentially-scaled spherical Bessel function</param>
-        /// <returns></returns>
-        public static double SphericalBesselY(double v, double z, bool expScaled = false)
+        /// <returns>The spherical Bessel function of the second kind.</returns>
+        public static double SphericalBesselY(double n, double z)
         {
-            const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
+            if (double.IsNaN(n) || double.IsNaN(z))
+            {
+                return double.NaN;
+            }
+
+            if (n < 0)
+            {
+                return double.NaN;
+            }
+
+            if (double.IsInfinity(z))
+            {
+                return 0;
+            }
+
+            if (z == 0)
+            {
+                return double.NegativeInfinity;
+            }
 
-            return rthpi * BesselY(v + 0.5, z, expScaled) / Math.Sqrt(z);
+            return Constants.SqrtPiOver2 * BesselY(n + 0.5, z, Scale.Unity) / Math.Sqrt(z);
         }
     }
 }