added scaled option to the Bessel functions in order to calculate I(n+1, z) / I(n, z) for large z.

src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs

@@ -241,5 +241,36 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
         }
 
         #endregion
+
+        #region Ratio of Bessel I : I(n + 1, z) / I(n, z)
+
+        [TestCase(0, 1, 1, 0.57495795977223997079, 0.35054769385125934266, 14)]
+        [TestCase(0, 10, 10, 0.97503660846868673171, 0.025655591609138262390, 14)]
+        [TestCase(0, 100, 100, 0.99750003174335729017, 0.0025062812447888833561, 12)]
+        [TestCase(0, 1E6, 1E6, 0.99999975000000000003, 2.5000006250003125000E-7, 8)] // cf. BesselI(0, 1e6 + 1e6 j) = 7.4E434290 - 6.9E434290 j and BesselI(1, 1e6 + 1e6 j) = 7.4E434290 - 6.9E434290 j
+        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);
+            AssertHelpers.AlmostEqualRelative(new Complex(cyr, cyi), actual, decimalPlaces);
+        }
+
+        #endregion
+
+        #region Ratio of Bessel K : K(n + 1, z) / K(n, z)
+
+        [TestCase(0, 1E-6, 1E-6, 38803.844378426554737, -34562.243460439510287, 14)]
+        [TestCase(0, 1, 1, 1.2397012468882428039, -0.20950920530836317445, 14)]
+        [TestCase(0, 10, 10, 1.0249731393042359657, -0.024405853995209668702, 14)]
+        [TestCase(0, 100, 100, 1.0024999692332050665, -0.0024937812450508461147, 12)]
+        [TestCase(0, 1E6, 1E6, 1.0000002500000000000, -2.4999993750003125000E-7, 9)]
+        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);
+            AssertHelpers.AlmostEqualRelative(new Complex(cyr, cyi), actual, decimalPlaces);
+        }
+
+        #endregion
     }
 }

src/Numerics/SpecialFunctions/Airy.cs

@@ -12,87 +12,127 @@ namespace MathNet.Numerics
     public static partial class SpecialFunctions
     {
         /// <summary>
-        /// Airy function Ai
+        /// Airy function Ai(z).
+        /// <p/>
+        /// If expScaled is true, returns Exp(zta) * Ai(z), where zta = (2/3) * z * Sqrt(z).
         /// </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)
+        public static Complex AiryAi(Complex z, bool expScaled = false)
         {
             var amos = new AmosWrapper();
-            return amos.Cairy(z);
+            return (expScaled) ? amos.ScaledCairy(z) : amos.Cairy(z);
         }
 
         /// <summary>
-        /// Airy function Ai
+        /// Airy function Ai(z).
+        /// <p/>
+        /// If expScaled is true, returns Exp(zta) * Ai(z), where zta = (2/3) * z * Sqrt(z).
         /// </summary>
-        /// <param name="x">The value to compute the Airy function of.</param>
+        /// <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 x)
+        public static double AiryAi(double z, bool expScaled = false)
         {
-            return AiryAi(new Complex(x, 0)).Real;
+            if (expScaled)
+            {
+                var amos = new AmosWrapper();
+                return amos.ScaledCairy(z);
+            }
+            else
+            {
+                return AiryAi(new Complex(z, 0), expScaled).Real;
+            }
         }
 
         /// <summary>
-        /// Derivative of the Airy function Ai
+        /// 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).
         /// </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)
+        public static Complex AiryAiPrime(Complex z, bool expScaled = false)
         {
             var amos = new AmosWrapper();
-            return amos.CairyPrime(z);
+            return (expScaled) ? amos.ScaledCairyPrime(z) : amos.CairyPrime(z);
         }
 
         /// <summary>
-        /// Derivative of the Airy function Ai
+        /// 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).
         /// </summary>
-        /// <param name="x">The value to compute the derivative of the Airy function of.</param>
+        /// <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 x)
+        public static double AiryAiPrime(double z, bool expScaled = false)
         {
-            return AiryAiPrime(new Complex(x, 0)).Real;
+            if (expScaled)
+            {
+                var amos = new AmosWrapper();
+                return amos.ScaledCairyPrime(z);
+            }
+            else
+            {
+                return AiryAiPrime(new Complex(z, 0), expScaled).Real;
+            }
         }
 
         /// <summary>
-        /// Airy function Bi(z)
+        /// 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).
         /// </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)
+        public static Complex AiryBi(Complex z, bool expScaled = false)
         {
             var amos = new AmosWrapper();
-            return amos.Cbiry(z);
+            return (expScaled) ? amos.ScaledCbiry(z) : amos.Cbiry(z);
         }
 
         /// <summary>
-        /// Airy function Bi(x)
+        /// 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).
         /// </summary>
-        /// <param name="x">The value to compute the Airy function of.</param>
+        /// <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 x)
+        public static double AiryBi(double z, bool expScaled = false)
         {
-            return AiryBi(new Complex(x, 0)).Real;
+            return AiryBi(new Complex(z, 0), expScaled).Real;
         }
 
         /// <summary>
-        /// Derivative of the Airy function Bi(z)
+        /// 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).
         /// </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)
+        public static Complex AiryBiPrime(Complex z, bool expScaled = false)
         {
             var amos = new AmosWrapper();
-            return amos.CbiryPrime(z);
+            return (expScaled) ? amos.ScaledCbiryPrime(z) : amos.CbiryPrime(z);
         }
 
         /// <summary>
-        /// Derivative of the Airy function Bi(z)
+        /// 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).
         /// </summary>
-        /// <param name="x">The value to compute the derivative of the Airy function of.</param>
+        /// <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 x)
+        public static double AiryBiPrime(double z, bool expScaled = false)
         {
-            return AiryBiPrime(new Complex(x, 0)).Real;
+            return AiryBiPrime(new Complex(z, 0), expScaled).Real;
         }
     }
 }

src/Numerics/SpecialFunctions/Amos/AmosWrapper.cs

@@ -14,6 +14,7 @@ namespace MathNet.Numerics
         {
             #region AiryAi
 
+            // Returns Ai(z)
             public Complex Cairy(Complex z)
             {
                 int id = 0;
@@ -29,6 +30,44 @@ namespace MathNet.Numerics
                 return new Complex(air, aii);
             }
 
+            // Returns Exp(zta) * Ai(z) where zta = (2/3) * z * Sqrt(z)
+            public Complex ScaledCairy(Complex z)
+            {
+                int id = 0;
+                int kode = 2;
+                int nz = 0;
+                int ierr = 0;
+
+                double air = double.NaN;
+                double aii = double.NaN;
+
+                var amos = new AmosHelper();
+                amos.zairy(z.Real, z.Imaginary, id, kode, ref air, ref aii, ref nz, ref ierr);
+                return new Complex(air, aii);
+            }
+
+            // Returns Exp(zta) * Ai(z) where zta = (2/3) * z * Sqrt(z)
+            public double ScaledCairy(double z)
+            {
+                if (z < 0)
+                {
+                    return double.NaN;
+                }
+
+                int id = 0;
+                int kode = 2;
+                int nz = 0;
+                int ierr = 0;
+
+                double air = double.NaN;
+                double aii = double.NaN;
+
+                var amos = new AmosHelper();
+                amos.zairy(z, 0.0, id, kode, ref air, ref aii, ref nz, ref ierr);
+                return air;
+            }
+
+            // Returns d/dz Ai(z)
             public Complex CairyPrime(Complex z)
             {
                 int id = 1;
@@ -44,10 +83,48 @@ namespace MathNet.Numerics
                 return new Complex(air, aii);
             }
 
+            // Returns Exp(zta) * d/dz Ai(z) where zta = (2/3) * z * Sqrt(z)
+            public Complex ScaledCairyPrime(Complex z)
+            {
+                int id = 1;
+                int kode = 2;
+                int nz = 0;
+                int ierr = 0;
+
+                double air = double.NaN;
+                double aii = double.NaN;
+
+                var amos = new AmosHelper();
+                amos.zairy(z.Real, z.Imaginary, id, kode, ref air, ref aii, ref nz, ref ierr);
+                return new Complex(air, aii);
+            }
+
+            // Returns Exp(zta) * d/dz Ai(z) where zta = (2/3) * z * Sqrt(z)
+            public double ScaledCairyPrime(double z)
+            {
+                if (z < 0)
+                {
+                    return double.NaN;
+                }
+
+                int id = 1;
+                int kode = 2;
+                int nz = 0;
+                int ierr = 0;
+
+                double air = double.NaN;
+                double aii = double.NaN;
+
+                var amos = new AmosHelper();
+                amos.zairy(z, 0.0, id, kode, ref air, ref aii, ref nz, ref ierr);
+                return air;
+            }
+
             #endregion
 
             #region AiryBi
 
+            // Returns Bi(z)
             public Complex Cbiry(Complex z)
             {
                 int id = 0;
@@ -63,6 +140,23 @@ namespace MathNet.Numerics
                 return new Complex(bir, bii);
             }
 
+            // Returns Exp(-axzta) * Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real)
+            public Complex ScaledCbiry(Complex z)
+            {
+                int id = 0;
+                int kode = 2;
+                int nz = 0;
+                int ierr = 0;
+
+                double bir = double.NaN;
+                double bii = double.NaN;
+
+                var amos = new AmosHelper();
+                amos.zbiry(z.Real, z.Imaginary, id, kode, ref bir, ref bii, ref nz, ref ierr);
+                return new Complex(bir, bii);
+            }
+
+            // Returns d/dz Bi(z)
             public Complex CbiryPrime(Complex z)
             {
                 int id = 1;
@@ -78,10 +172,27 @@ namespace MathNet.Numerics
                 return new Complex(bipr, bipi);
             }
 
+            // Returns Exp(-axzta) * d/dz Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real)
+            public Complex ScaledCbiryPrime(Complex z)
+            {
+                int id = 1;
+                int kode = 2;
+                int nz = 0;
+                int ierr = 0;
+
+                double bipr = double.NaN;
+                double bipi = double.NaN;
+
+                var amos = new AmosHelper();
+                amos.zbiry(z.Real, z.Imaginary, id, kode, ref bipr, ref bipi, ref nz, ref ierr);
+                return new Complex(bipr, bipi);
+            }
+
             #endregion
 
             #region BesselJ
 
+            // Return J(v, z)
             public Complex Cbesj(double v, Complex z)
             {
                 if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
@@ -115,7 +226,7 @@ namespace MathNet.Numerics
                 if (ierr == 2)
                 {
                     //overflow
-                    cyj = CbesjScaled(v, z);
+                    cyj = ScaledCbesj(v, z);
                     cyj = new Complex(cyj.Real * double.PositiveInfinity, cyj.Imaginary * double.PositiveInfinity);
                 }
 
@@ -145,7 +256,19 @@ namespace MathNet.Numerics
                 return cyj;
             }
 
-            private Complex CbesjScaled(double v, Complex z)
+            // Return J(v, z)
+            public double Cbesj(double v, double z)
+            {
+                if (z < 0 && v != (int)v)
+                {
+                    return double.NaN;
+                }
+
+                return Cbesj(v, new Complex(z, 0)).Real;                 
+            }
+
+            // Return Exp(-Abs(y)) * J(v, z) where y = z.Imaginary
+            public Complex ScaledCbesj(double v, Complex z)
             {
                 if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
                 {
@@ -201,10 +324,22 @@ namespace MathNet.Numerics
                 return cyj;
             }
 
+            // Return Exp(-Abs(y)) * J(v, z) where y = z.Imaginary
+            public double ScaledCbesj(double v, double z)
+            {
+                if (z < 0 && v != (int)v)
+                {
+                    return double.NaN;
+                }
+
+                return ScaledCbesj(v, new Complex(z, 0)).Real;
+            }
+
             #endregion
 
             #region BesselY
 
+            // Return Y(v, z)
             public Complex Cbesy(double v, Complex z)
             {
                 if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
@@ -279,7 +414,8 @@ namespace MathNet.Numerics
                 return cyy;
             }
 
-            public double CbesyReal(double v, double x)
+            // Return Y(v, z)
+            public double Cbesy(double v, double x)
             {
                 if (x < 0.0)
                 {
@@ -290,10 +426,87 @@ namespace MathNet.Numerics
                 return Cbesy(v, z).Real;
             }
 
+            // Return Exp(-Abs(y)) * Y(v, z) where y = z.Imaginary
+            public Complex ScaledCbesy(double v, Complex z)
+            {
+                if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
+                {
+                    return new Complex(double.NaN, double.NaN);
+                }
+
+                int sign = 1;
+                if (v < 0)
+                {
+                    v = -v;
+                    sign = -1;
+                }
+
+                int n = 1;
+                int kode = 2;
+                int nz = 0;
+                int ierr = 0;
+
+                double[] cyyr = new double[n];
+                double[] cyyi = new double[n];
+                double[] cwrkr = new double[n];
+                double[] cwrki = new double[n];
+                for (int i = 0; i < n; i++)
+                {
+                    cyyr[i] = double.NaN;
+                    cyyi[i] = double.NaN;
+                    cwrkr[i] = double.NaN;
+                    cwrki[i] = double.NaN;
+                }
+
+                var amos = new AmosHelper();
+                amos.zbesy(z.Real, z.Imaginary, v, kode, n, cyyr, cyyi, ref nz, cwrkr, cwrki, ref ierr);
+                Complex cyy = new Complex(cyyr[0], cyyi[0]);
+
+                if (ierr == 2)
+                {
+                    if (z.Real >= 0 && z.Imaginary == 0)
+                    {
+                        //overflow
+                        cyy = new Complex(double.PositiveInfinity, 0);
+                    }
+                }
+
+                if (sign == -1)
+                {
+                    if (!ReflectJY(ref cyy, v))
+                    {
+                        double[] cyjr = new double[n];
+                        double[] cyji = new double[n];
+                        for (int i = 0; i < n; i++)
+                        {
+                            cyjr[i] = double.NaN;
+                            cyji[i] = double.NaN;
+                        }
+
+                        amos.zbesj(z.Real, z.Imaginary, v, kode, n, cyjr, cyji, ref nz, ref ierr);
+                        Complex cyj = new Complex(cyjr[0], cyji[0]);
+                        cyy = RotateJY(cyy, cyj, -v);
+                    }
+                }
+                return cyy;
+            }
+
+            // Return Exp(-Abs(y)) * Y(v, z) where y = z.Imaginary
+            public double ScaledCbesy(double v, double x)
+            {
+                if (x < 0)
+                {
+                    return double.NaN;
+                }
+
+                return ScaledCbesy(v, new Complex(x, 0)).Real;
+            }
+
             #endregion
 
             #region BesselI
 
+            // Returns I(v, z)
             public Complex Cbesi(double v, Complex z)
             {
                 if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
@@ -337,7 +550,7 @@ namespace MathNet.Numerics
                     }
                     else
                     {
-                        cyi = CbesiScaled(v * sign, z);
+                        cyi = ScaledCbesi(v * sign, z);
                         cyi = new Complex(cyi.Real * double.PositiveInfinity, cyi.Imaginary * double.PositiveInfinity);
                     }
                 }
@@ -358,7 +571,8 @@ namespace MathNet.Numerics
                 return cyi;
             }
 
-            private Complex CbesiScaled(double v, Complex z)
+            // Return Exp(-Abs(x)) * I(v, z) where x = z.Real
+            public Complex ScaledCbesi(double v, Complex z)
             {
                 if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
                 {
@@ -412,10 +626,22 @@ namespace MathNet.Numerics
                 return cyi;
             }
 
+            // Return Exp(-Abs(x)) * I(v, z) where x = z.Real
+            public double ScaledCbesi(double v, double x)
+            {
+                if (v != Math.Floor(v) && x < 0)
+                {
+                    return double.NaN;
+                }
+
+                return ScaledCbesi(v, new Complex(x, 0)).Real;
+            }
+
             #endregion
 
             #region BesselK
 
+            // Returns K(v, z)
             public Complex Cbesk(double v, Complex z)
             {
                 if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
@@ -463,8 +689,9 @@ namespace MathNet.Numerics
 
                 return cyk;
             }
-            
-            public double CbeskReal(double v, double z)
+
+            // Returns K(v, z)
+            public double Cbesk(double v, double z)
             {
                 if (z < 0)
                 {
@@ -489,10 +716,73 @@ namespace MathNet.Numerics
                 }
             }
 
+            // returns Exp(z) * K(v, z)
+            public Complex ScaledCbesk(double v, Complex z)
+            {
+                if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
+                {
+                    return new Complex(double.NaN, double.NaN);
+                }
+
+                if (v < 0)
+                {
+                    //K_v == K_{-v} even for non-integer v
+                    v = -v;
+                }
+
+                int n = 1;
+                int kode = 2;
+                int nz = 0;
+                int ierr = 0;
+
+                double[] cykr = new double[n];
+                double[] cyki = new double[n];
+                for (int i = 0; i < n; i++)
+                {
+                    cykr[i] = double.NaN;
+                    cyki[i] = double.NaN;
+                }
+
+                var amos = new AmosHelper();
+
+                amos.zbesk(z.Real, z.Imaginary, v, kode, n, cykr, cyki, ref nz, ref ierr);
+                Complex cyk = new Complex(cykr[0], cyki[0]);
+
+                if (ierr == 2)
+                {
+                    if (z.Real >= 0 && z.Imaginary == 0)
+                    {
+                        //overflow
+                        cyk = new Complex(double.PositiveInfinity, 0);
+                    }
+                }
+
+                return cyk;
+            }
+
+            // returns Exp(z) * K(v, z)
+            public double ScaledCbesk(double v, double z)
+            {
+                if (z < 0)
+                {
+                    return double.NaN;
+                }
+                else if (z == 0)
+                {
+                    return double.PositiveInfinity;
+                }
+                else
+                {
+                    Complex w = new Complex(z, 0);
+                    return ScaledCbesk(v, w).Real;
+                }
+            }
+
             #endregion
 
             #region HankelH1
 
+            // Returns H1(v, z)
             public Complex Cbesh1(double v, Complex z)
             {
                 if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
@@ -533,10 +823,52 @@ namespace MathNet.Numerics
                 return cyh;
             }
 
+            // Returns Exp(-z * j) * H1(n, z) where j = Sqrt(-1)
+            public Complex ScaledCbesh1(double v, Complex z)
+            {
+                if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
+                {
+                    return new Complex(double.NaN, double.NaN); ;
+                }
+
+                int n = 1;
+                int kode = 2;
+                int m = 1;
+                int nz = 0;
+                int ierr = 0;
+
+                double[] cyhr = new double[n];
+                double[] cyhi = new double[n];
+                for (int i = 0; i < n; i++)
+                {
+                    cyhr[i] = double.NaN;
+                    cyhi[i] = double.NaN;
+                }
+
+                int sign = 1;
+                if (v < 0)
+                {
+                    v = -v;
+                    sign = -1;
+                }
+
+                var amos = new AmosHelper();
+                amos.zbesh(z.Real, z.Imaginary, v, kode, m, n, cyhr, cyhi, ref nz, ref ierr);
+                Complex cyh = new Complex(cyhr[0], cyhi[0]);
+
+                if (sign == -1)
+                {
+                    cyh = Rotate(cyh, v);
+                }
+
+                return cyh;
+            }
+
             #endregion
 
             #region HankelH2
 
+            // Returns H2(v, z)
             public Complex Cbesh2(double v, Complex z)
             {
                 if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
@@ -581,6 +913,51 @@ namespace MathNet.Numerics
                 return cyh;
             }
 
+            // Returns Exp(z * j) * H2(n, z) where j = Sqrt(-1)
+            public Complex ScaledCbesh2(double v, Complex z)
+            {
+                if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
+                {
+                    return new Complex(double.NaN, double.NaN);
+                }
+
+                if (v == 0 && z.Real == 0 && z.Imaginary == 0)
+                {
+                    return new Complex(double.NaN, double.NaN); // ComplexInfinity
+                }
+
+                int n = 1;
+                int kode = 2;
+                int m = 2;
+                int nz = 0;
+                int ierr = 0;
+
+                double[] cyhr = new double[n];
+                double[] cyhi = new double[n];
+                for (int i = 0; i < n; i++)
+                {
+                    cyhr[i] = double.NaN;
+                    cyhi[i] = double.NaN;
+                }
+
+                int sign = 1;
+                if (v < 0)
+                {
+                    v = -v;
+                    sign = -1;
+                }
+
+                var amos = new AmosHelper();
+                amos.zbesh(z.Real, z.Imaginary, v, kode, m, n, cyhr, cyhi, ref nz, ref ierr);
+                Complex cyh = new Complex(cyhr[0], cyhi[0]);
+
+                if (sign == -1)
+                {
+                    cyh = Rotate(cyh, -v);
+                }
+                return cyh;
+            }
+
             #endregion
 
             #region utilities

src/Numerics/SpecialFunctions/Bessel.cs

@@ -12,97 +12,130 @@ namespace MathNet.Numerics
     public static partial class SpecialFunctions
     {
         /// <summary>
-        /// Bessel function of the first kind
+        /// 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.
         /// </summary>
-        /// <param name="n">The order of the Bessel function</param>
+        /// <param name="v">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 n, Complex z)
+        public static Complex BesselJ(double v, Complex z, bool expScaled = false)
         {
             var amos = new AmosWrapper();
-            return amos.Cbesj(n, z);
+            return (expScaled) ? amos.ScaledCbesj(v, z) : amos.Cbesj(v, z);
         }
 
         /// <summary>
-        /// Bessel function of the first kind
+        /// 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.
         /// </summary>
-        /// <param name="n">The order of the Bessel function</param>
-        /// <param name="x">The value to compute the Bessel function of.</param>
+        /// <param name="v">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 n, double x)
+        public static double BesselJ(double v, double z, bool expScaled = false)
         {
-            return BesselJ(n, new Complex(x, 0)).Real;
+            var amos = new AmosWrapper();
+            return (expScaled) ? amos.ScaledCbesj(v, z) : amos.Cbesj(v, z);
         }
 
         /// <summary>
-        /// Bessel function of the second kind
+        /// 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.
         /// </summary>
-        /// <param name="n">The order of the Bessel function</param>
+        /// <param name="v">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 n, Complex z)
+        public static Complex BesselY(double v, Complex z, bool expScaled = false)
         {
             var amos = new AmosWrapper();
-            return amos.Cbesy(n, z);
+            return (expScaled) ? amos.ScaledCbesy(v, z) : amos.Cbesy(v, z);
         }
 
         /// <summary>
-        /// Bessel function of the second kind
+        /// 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.
         /// </summary>
-        /// <param name="n">The order of the Bessel function</param>
-        /// <param name="x">The value to compute the Bessel function of.</param>
+        /// <param name="v">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 n, double x)
+        public static double BesselY(double v, double z, bool expScaled = false)
         {
             var amos = new AmosWrapper();
-            return amos.CbesyReal(n, x);
+            return (expScaled) ? amos.ScaledCbesy(v, z) : amos.Cbesy(v, z);
         }
 
         /// <summary>
-        /// Modified Bessel function of the first kind, of order n 
+        /// 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.
         /// </summary>
-        /// <param name="n">The order of the Bessel function</param>
+        /// <param name="v">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 n, Complex z)
+        public static Complex BesselI(double v, Complex z, bool expScaled = false)
         {
             var amos = new AmosWrapper();
-            return amos.Cbesi(n, z);
+            return (expScaled) ? amos.ScaledCbesi(v, z) : amos.Cbesi(v, z);
         }
 
         /// <summary>
-        /// Modified Bessel function of the first kind, of order n  
+        /// 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.
         /// </summary>
-        /// <param name="n">The order of the Bessel function</param>
-        /// <param name="x">The value to compute the Bessel function of.</param>
+        /// <param name="v">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 n, double x)
+        public static double BesselI(double v, double z, bool expScaled = false)
         {
-            return BesselI(n, new Complex(x, 0)).Real;
+            if (expScaled)
+            {
+                var amos = new AmosWrapper();
+                return amos.ScaledCbesi(v, z);
+            }
+            else
+            {
+                return BesselI(v, new Complex(z, 0), expScaled).Real;
+            }
         }
 
         /// <summary>
-        /// Modified Bessel function of the second kind
+        /// Modified Bessel function of the second kind, K(v, z).
+        /// <p/>
+        /// If expScaled is true, returns Exp(z) * K(v, z).
         /// </summary>
-        /// <param name="n">The order of the Bessel function</param>
+        /// <param name="v">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 n, Complex z)
+        public static Complex BesselK(double v, Complex z, bool expScaled = false)
         {
             var amos = new AmosWrapper();
-            return amos.Cbesk(n, z);
+            return (expScaled) ? amos.ScaledCbesk(v, z) : amos.Cbesk(v, z);
         }
 
         /// <summary>
-        /// Modified Bessel function of the second kind
+        /// Modified Bessel function of the second kind, K(v, z).
+        /// <p/>
+        /// If expScaled is true, returns Exp(z) * K(v, z).
         /// </summary>
-        /// <param name="n">The order of the Bessel function</param>
-        /// <param name="x">The value to compute the Bessel function of.</param>
+        /// <param name="v">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 n, double x)
+        public static double BesselK(double v, double z, bool expScaled = false)
         {
             var amos = new AmosWrapper();
-            return amos.CbeskReal(n, x);
+            return (expScaled) ? amos.ScaledCbesk(v, z) : amos.Cbesk(v, z);
         }
     }
 }

src/Numerics/SpecialFunctions/Hankel.cs

@@ -12,49 +12,61 @@ namespace MathNet.Numerics
     public static partial class SpecialFunctions
     {
         /// <summary>
-        /// Hankel function of the first kind
+        /// 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).
         /// </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)
+        public static Complex HankelH1(double n, Complex z, bool expScaled = false)
         {
             var amos = new AmosWrapper();
-            return amos.Cbesh1(n, z);
+            return (expScaled) ? amos.ScaledCbesh1(n, z) : amos.Cbesh1(n, z);
         }
 
         /// <summary>
-        /// Hankel function of the first kind
+        /// 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).
         /// </summary>
         /// <param name="n">The order of the Bessel function</param>
-        /// <param name="x">The value to compute the Bessel function of.</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 x)
+        public static double HankelH1(double n, double z, bool expScaled = false)
         {
-            return HankelH1(n, new Complex(x, 0)).Real;
+            return HankelH1(n, new Complex(z, 0), expScaled).Real;
         }
 
         /// <summary>
-        /// Hankel function of the second kind
+        /// 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).
         /// </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)
+        public static Complex HankelH2(double n, Complex z, bool expScaled = false)
         {
             var amos = new AmosWrapper();
-            return amos.Cbesh2(n, z);
+            return (expScaled) ? amos.ScaledCbesh2(n, z) :amos.Cbesh2(n, z);
         }
 
         /// <summary>
-        /// Hankel function of the second kind
+        /// 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).
         /// </summary>
         /// <param name="n">The order of the Bessel function</param>
-        /// <param name="x">The value to compute the Bessel function of.</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 x)
+        public static double HankelH2(double n, double z, bool expScaled = false)
         {
-            return HankelH2(n, new Complex(x, 0)).Real;
+            return HankelH2(n, new Complex(z, 0), expScaled).Real;
         }
     }
 }

src/Numerics/SpecialFunctions/SphericalBessel.cs

@@ -12,55 +12,67 @@ namespace MathNet.Numerics
     public static partial class SpecialFunctions
     {
         /// <summary>
-        /// Spherical Bessel function of the first kind
+        /// 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.
         /// </summary>
-        /// <param name="n">The order of the spherical Bessel function</param>
+        /// <param name="v">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 n, Complex z)
+        public static Complex SphericalBesselJ(double v, Complex z, bool expScaled = false)
         {
             const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
 
-            return rthpi * BesselJ(n + 0.5, z) / Complex.Sqrt(z);
+            return rthpi * BesselJ(v + 0.5, z, expScaled) / Complex.Sqrt(z);
         }
 
         /// <summary>
-        /// Spherical Bessel function of the first kind
+        /// 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.
         /// </summary>
-        /// <param name="n">The order of the spherical Bessel function</param>
-        /// <param name="x">The value to compute the spherical Bessel function of.</param>
+        /// <param name="v">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 n, double x)
+        public static double SphericalBesselJ(double v, double z, bool expScaled = false)
         {
             const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
 
-            return rthpi * BesselJ(n + 0.5, x) / Math.Sqrt(x);
+            return rthpi * BesselJ(v + 0.5, z, expScaled) / Math.Sqrt(z);
         }
 
         /// <summary>
-        /// Spherical Bessel function of the second kind
+        /// 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.
         /// </summary>
-        /// <param name="n">The order of the spherical Bessel function</param>
+        /// <param name="v">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 n, Complex z)
+        public static Complex SphericalBesselY(double v, Complex z, bool expScaled = false)
         {
             const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
 
-            return rthpi * BesselY(n + 0.5, z) / Complex.Sqrt(z);
+            return rthpi * BesselY(v + 0.5, z, expScaled) / Complex.Sqrt(z);
         }
 
         /// <summary>
-        /// Spherical Bessel function of the second kind
+        /// 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.
         /// </summary>
-        /// <param name="n">The order of the spherical Bessel function</param>
-        /// <param name="x">The value to compute the spherical Bessel function of.</param>
+        /// <param name="v">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 n, double x)
+        public static double SphericalBesselY(double v, double z, bool expScaled = false)
         {
             const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
 
-            return rthpi * BesselY(n + 0.5, x) / Math.Sqrt(x);
+            return rthpi * BesselY(v + 0.5, z, expScaled) / Math.Sqrt(z);
         }
     }
 }