From 1a5c1fba740dd72835f41912c4bf33147376107e Mon Sep 17 00:00:00 2001
From: diluculo <diluculo@gmail.com>
Date: Tue, 6 Nov 2018 19:53:42 +0900
Subject: [PATCH] Special Functions: cleanup Bessel-type functions - Scale enum
 is added and replaced expScaled boolean flag. - Scaled Bessel functions are
 added. - HankelH1 and HankelH2    - allow only complex argument because they
 return usually complex values. - SphericalBesselJ and SphericalBesselY:    -
 scaling options are removed.

---
 .../SpecialFunctionsTests/BesselTests.cs      |   4 +-
 src/Numerics/Constants.cs                     |   3 +
 src/Numerics/SpecialFunctions/Airy.cs         | 214 ++++++++++------
 src/Numerics/SpecialFunctions/Bessel.cs       | 231 ++++++++++++------
 src/Numerics/SpecialFunctions/Hankel.cs       |  68 +++---
 src/Numerics/SpecialFunctions/Options.cs      |  23 ++
 .../SpecialFunctions/SphericalBessel.cs       | 126 +++++++---
 7 files changed, 454 insertions(+), 215 deletions(-)
 create mode 100644 src/Numerics/SpecialFunctions/Options.cs

diff --git a/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs b/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs
index 60ed1416..f8dd3451 100644
--- a/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs
+++ b/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs
@@ -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/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 1e69093c..c2f7936a 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 d9c15f7f..0bbb8945 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 3ad42386..6eaf4156 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/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 afa1a1b8..2858cc5a 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);
         }
     }
 }