using System; using Complex = System.Numerics.Complex; namespace IStation.Numerics { // References: // [1] https://github.com/scipy/scipy/blob/master/scipy/special/amos_wrappers.c public static partial class SpecialFunctions { private static class Amos { #region AiryAi // Returns Ai(z) public static Complex Cairy(Complex z) { int id = 0; int kode = 1; int nz = 0; int ierr = 0; double air = double.NaN; double aii = double.NaN; AmosHelper.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 static Complex ScaledCairy(Complex z) { int id = 0; int kode = 2; int nz = 0; int ierr = 0; double air = double.NaN; double aii = double.NaN; AmosHelper.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 static 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; AmosHelper.zairy(z, 0.0, id, kode, ref air, ref aii, ref nz, ref ierr); return air; } // Returns d/dz Ai(z) public static Complex CairyPrime(Complex z) { int id = 1; int kode = 1; int nz = 0; int ierr = 0; double air = double.NaN; double aii = double.NaN; AmosHelper.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 static Complex ScaledCairyPrime(Complex z) { int id = 1; int kode = 2; int nz = 0; int ierr = 0; double air = double.NaN; double aii = double.NaN; AmosHelper.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 static 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; AmosHelper.zairy(z, 0.0, id, kode, ref air, ref aii, ref nz, ref ierr); return air; } #endregion #region AiryBi // Returns Bi(z) public static Complex Cbiry(Complex z) { int id = 0; int kode = 1; int nz = 0; int ierr = 0; double bir = double.NaN; double bii = double.NaN; AmosHelper.zbiry(z.Real, z.Imaginary, id, kode, ref bir, ref bii, ref nz, ref ierr); return new Complex(bir, bii); } // Returns Exp(-axzta) * Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real) public static Complex ScaledCbiry(Complex z) { int id = 0; int kode = 2; int nz = 0; int ierr = 0; double bir = double.NaN; double bii = double.NaN; AmosHelper.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 static Complex CbiryPrime(Complex z) { int id = 1; int kode = 1; int nz = 0; int ierr = 0; double bipr = double.NaN; double bipi = double.NaN; AmosHelper.zbiry(z.Real, z.Imaginary, id, kode, ref bipr, ref bipi, ref nz, ref ierr); 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 static Complex ScaledCbiryPrime(Complex z) { int id = 1; int kode = 2; int nz = 0; int ierr = 0; double bipr = double.NaN; double bipi = double.NaN; AmosHelper.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 static Complex Cbesj(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 = 1; int nz = 0; int ierr = 0; 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; } AmosHelper.zbesj(z.Real, z.Imaginary, v, kode, n, cyjr, cyji, ref nz, ref ierr); Complex cyj = new Complex(cyjr[0], cyji[0]); if (ierr == 2) { //overflow cyj = ScaledCbesj(v, z); cyj = new Complex(cyj.Real * double.PositiveInfinity, cyj.Imaginary * double.PositiveInfinity); } if (sign == -1) { if (!ReflectJY(ref cyj, v)) { 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; } AmosHelper.zbesy(z.Real, z.Imaginary, v, kode, n, cyyr, cyyi, ref nz, cwrkr, cwrki, ref ierr); Complex cyy = new Complex(cyyr[0], cyyi[0]); cyj = RotateJY(cyj, cyy, v); } } return cyj; } // Return J(v, z) public static 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 static Complex ScaledCbesj(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[] cyjr = new double[n]; double[] cyji = new double[n]; for (int i = 0; i < n; i++) { cyjr[i] = double.NaN; cyji[i] = double.NaN; } AmosHelper.zbesj(z.Real, z.Imaginary, v, kode, n, cyjr, cyji, ref nz, ref ierr); Complex cyj = new Complex(cyjr[0], cyji[0]); if (sign == -1) { if (!ReflectJY(ref cyj, v)) { double[] cyyr = new double[n]; double[] cyyi = new double[n]; double[] cworkr = new double[n]; double[] cworki = new double[n]; for (int i = 0; i < n; i++) { cyyr[i] = double.NaN; cyyi[i] = double.NaN; cworkr[i] = double.NaN; cworki[i] = double.NaN; } AmosHelper.zbesy(z.Real, z.Imaginary, v, kode, n, cyyr, cyyi, ref nz, cworkr, cworki, ref ierr); Complex cyy = new Complex(cyyr[0], cyyi[0]); cyj = RotateJY(cyj, cyy, v); } } return cyj; } // Return Exp(-Abs(y)) * J(v, z) where y = z.Imaginary public static 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 static Complex Cbesy(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 = 1; int nz = 0; int ierr = 0; Complex cyy; if (z.Real == 0 && z.Imaginary == 0) { //overflow cyy = new Complex(double.NegativeInfinity, 0); } else { double[] cyyr = new double[n]; double[] cyyi = new double[n]; double[] cworkr = new double[n]; double[] cworki = new double[n]; for (int i = 0; i < n; i++) { cyyr[i] = double.NaN; cyyi[i] = double.NaN; cworkr[i] = double.NaN; cworki[i] = double.NaN; } AmosHelper.zbesy(z.Real, z.Imaginary, v, kode, n, cyyr, cyyi, ref nz, cworkr, cworki, ref ierr); cyy = new Complex(cyyr[0], cyyi[0]); if (ierr == 2) { if (z.Real >= 0 && z.Imaginary == 0) { //overflow cyy = new Complex(double.NegativeInfinity, 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; } AmosHelper.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 Y(v, z) public static double Cbesy(double v, double x) { if (x < 0.0) { return double.NaN; } Complex z = new Complex(x, 0.0); return Cbesy(v, z).Real; } // Return Exp(-Abs(y)) * Y(v, z) where y = z.Imaginary public static 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; } AmosHelper.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; } AmosHelper.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 static 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 static Complex Cbesi(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 = 1; int nz = 0; int ierr = 0; double[] cyir = new double[n]; double[] cyii = new double[n]; for (int i = 0; i < n; i++) { cyir[i] = double.NaN; cyii[i] = double.NaN; } AmosHelper.zbesi(z.Real, z.Imaginary, v, kode, n, cyir, cyii, ref nz, ref ierr); Complex cyi = new Complex(cyir[0], cyii[0]); if (ierr == 2) { //overflow if (z.Imaginary == 0 && (z.Real >= 0 || v == Math.Floor(v))) { if (z.Real < 0 && v / 2 != Math.Floor(v / 2)) cyi = new Complex(double.NegativeInfinity, 0); else cyi = new Complex(double.PositiveInfinity, 0); } else { cyi = ScaledCbesi(v * sign, z); cyi = new Complex(cyi.Real * double.PositiveInfinity, cyi.Imaginary * double.PositiveInfinity); } } if (sign == -1) { if (!ReflectI(cyi, v)) { double[] cykr = new double[n]; double[] cyki = new double[n]; AmosHelper.zbesk(z.Real, z.Imaginary, v, kode, n, cykr, cyki, ref nz, ref ierr); Complex cyk = new Complex(cykr[0], cyki[0]); cyi = RotateI(cyi, cyk, v); } } return cyi; } // Return Exp(-Abs(x)) * I(v, z) where x = z.Real public static Complex ScaledCbesi(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[] cyir = new double[n]; double[] cyii = new double[n]; for (int i = 0; i < n; i++) { cyir[i] = double.NaN; cyii[i] = double.NaN; } AmosHelper.zbesi(z.Real, z.Imaginary, v, kode, n, cyir, cyii, ref nz, ref ierr); Complex cyi = new Complex(cyir[0], cyii[0]); if (sign == -1) { if (!ReflectI(cyi, v)) { double[] cykr = new double[n]; double[] cyki = new double[n]; AmosHelper.zbesk(z.Real, z.Imaginary, v, kode, n, cykr, cyki, ref nz, ref ierr); Complex cyk = new Complex(cykr[0], cyki[0]); //adjust scaling to match zbesi cyk = Rotate(cyk, -z.Imaginary / Math.PI); if (z.Real > 0) { cyk = new Complex(cyk.Real * Math.Exp(-2 * z.Real), cyk.Imaginary * Math.Exp(-2 * z.Real)); } //v -> -v cyi = RotateI(cyi, cyk, v); } } return cyi; } // Return Exp(-Abs(x)) * I(v, z) where x = z.Real public static 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 static Complex Cbesk(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 = 1; 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; } AmosHelper.zbesk(z.Real, z.Imaginary, v, kode, n, cykr, cyki, ref nz, ref ierr); Complex cyk = new Complex(cykr[0], cyki[0]); if (ierr == 1) { if (z.Real == 0.0 && z.Imaginary == 0.0) { cyk = new Complex(double.PositiveInfinity, 0); } } else if (ierr == 2) { if (z.Real >= 0 && z.Imaginary == 0) { //overflow cyk = new Complex(double.PositiveInfinity, 0); } } return cyk; } // Returns K(v, z) public static double Cbesk(double v, double z) { if (z < 0) { return double.NaN; } else if (z == 0) { return double.PositiveInfinity; } else if (z > 710 * (1 + Math.Abs(v))) { // Underflow. See uniform expansion https://dlmf.nist.gov/10.41 // This condition is not a strict bound (it can underflow earlier), // rather, we are here working around a restriction in AMOS. return 0; } else { Complex w = new Complex(z, 0); return Cbesk(v, w).Real; } } // returns Exp(z) * K(v, z) public static 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; } AmosHelper.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 static 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 static Complex Cbesh1(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 = 1; 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; } AmosHelper.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; } // Returns Exp(-z * j) * H1(n, z) where j = Sqrt(-1) public static 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; } AmosHelper.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 static Complex Cbesh2(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 = 1; 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; } AmosHelper.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; } // Returns Exp(z * j) * H2(n, z) where j = Sqrt(-1) public static 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; } AmosHelper.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 private static double SinPi(double x) { if (Math.Floor(x) == x && Math.Abs(x) < 1.0e14) { //Return 0 when at exact zero, as long as the floating point number is //small enough to distinguish integer points from other points. return 0; } return Math.Sin(Math.PI * x); } private static double CosPi(double x) { if (Math.Floor(x + 0.5) == x + 0.5 && Math.Abs(x) < 1.0E14) { //Return 0 when at exact zero, as long as the floating point number is //small enough to distinguish integer points from other points. return 0; } return Math.Cos(Math.PI * x); } private static Complex Rotate(Complex z, double v) { double c = CosPi(v); double s = SinPi(v); return new Complex(z.Real * c - z.Imaginary * s, z.Real * s + z.Imaginary * c); } private static Complex RotateJY(Complex j, Complex y, double v) { double c = CosPi(v); double s = SinPi(v); return new Complex(j.Real * c - y.Real * s, j.Imaginary * c - y.Imaginary * s); } private static bool ReflectJY(ref Complex jy, double v) { //NB: Y_v may be huge near negative integers -- so handle exact // integers carefully if (v != Math.Floor(v)) { return false; } int i = (int)(v - 16384.0 * Math.Floor(v / 16384.0)); if (i % 2 == 1) { jy = new Complex(-jy.Real, -jy.Imaginary); } return true; } private static bool ReflectI(Complex ik, double v) { if (v != Math.Floor(v)) { return false; } return true; //I is symmetric for integer v } private static Complex RotateI(Complex i, Complex k, double v) { double s = Math.Sin(v * Math.PI) * (2.0 / Math.PI); return new Complex(i.Real + s * k.Real, i.Imaginary + s * k.Imaginary); } #endregion } } }