0

I need to use the bessel function of the first kind in unity3d (using c#). After trying to use the mathnet.numerics .dll in unity I found many errors apparently due to unity3d not supporting .NET 4. I decided to then use jetbrains dotPeek decompiler to find the bessel functions source.

However the bessel function gives the wrong output and to be honest the maths is above my paygrade so debugging it has been very difficult. Can anyone see any errors in what i have used.

for example, the output with x = 20 is 43558282.5571362 where is is supposed to be 0.167025

public static double BesselI0(double x)
{
    if (x < 0.0)
        x = -x;

    if (x <= 8.0)
    {
        double x1 = x / 2.0 - 2.0;
        return Math.Exp(x) * ChebyshevA(BesselI0A, x1);
    }
    else
    {
        double x1 = 32.0 / x - 2.0;
        return Math.Exp(x) * ChebyshevA(BesselI0B, x1) / Math.Sqrt(x);
    }   
}

internal static double ChebyshevA(double[] coefficients, double x)
{
    int num1 = 0;
    double[] numArray = coefficients;
    int index = num1;
    int num2 = 1;
    int num3 = index + num2;
    double num4 = numArray[index];
    double num5 = 0.0;
    int num6 = coefficients.Length - 1;
    double num7;
    do
    {
        num7 = num5;
        num5 = num4;
        num4 = x * num5 - num7 + coefficients[num3++];
    }
    while (--num6 > 0);
    return 0.5 * (num4 - num7);
}

private static readonly double[] BesselI0A = new double[30]
{
    -4.41534164647934E-18,
    3.33079451882224E-17,
    -2.43127984654795E-16,
    1.71539128555513E-15,
    -1.16853328779935E-14,
    7.67618549860494E-14,
    0.0 / 1.0,
    0.0 / 1.0,
    0.0 / 1.0,
    0.0 / 1.0,
    -5.18979560163526E-10,
    2.65982372468239E-09,
    -1.30002500998625E-08,
    6.04699502254192E-08,
    -2.67079385394061E-07,
    1.1173875391201E-06,
    -4.41673835845875E-06,
    1.64484480707289E-05,
    -5.7541950100821E-05,
    0.000188502885095842,
    -0.000576375574538582,
    0.00163947561694134,
    -0.00432430999505058,
    0.010546460394595,
    -0.0237374148058995,
    0.0493052842396707,
    -0.0949010970480476,
    0.171620901522209,
    -0.304682672343198,
    0.676795274409476
};

private static readonly double[] BesselI0B = new double[25]
{
    -7.23318048787475E-18,
    -4.83050448594418E-18,
    4.46562142029676E-17,
    3.46122286769746E-17,
    -2.82762398051658E-16,
    -3.42548561967722E-16,
    1.77256013305653E-15,
    3.81168066935262E-15,
    -9.55484669882831E-15,
    -4.15056934728722E-14,
    1.54008621752141E-14,
    0.0 / 1.0,
    0.0 / 1.0,
    0.0 / 1.0,
    0.0 / 1.0,
    0.0 / 1.0,
    0.0 / 1.0,
    4.94060238822497E-10,
    3.39623202570839E-09,
    2.26666899049818E-08,
    2.04891858946906E-07,
    2.89137052083476E-06,
    6.88975834691682E-05,
    0.00336911647825569,
    0.804490411014109
};
Mike Zboray
  • 39,828
  • 3
  • 90
  • 122
KMoore
  • 171
  • 1
  • 10
  • Without delving into specifics (admittedly my Calculus is beyond rusty), have you tried to verify that you have your order of operations correct? I've seen quite a few issues with mathematical functions returning unexpected results simply due to order of operations. – Leon Newswanger Jun 22 '14 at 19:04
  • I'll try verify that. thanks – KMoore Jun 22 '14 at 19:18
  • 1
    Are you sure that is the correct function? The Bessel functions of the first and second kinds are usually denote by J and Y, while the *modified* Bessel functions of the first and second kinds are usually denoted by I and K. – Mike Zboray Jun 22 '14 at 19:34
  • Also [Math.net](https://github.com/mathnet/mathnet-numerics) ([here too](https://www.nuget.org/packages/MathNet.Numerics/)) claims they support .NET 3.5. – Mike Zboray Jun 22 '14 at 19:40
  • Sorry all, I was looking through the documentation. The code was for the modified bessel function. Not the bessel function. Thanks – KMoore Jun 22 '14 at 20:02

1 Answers1

0

If you are trying to use Bessel function, not the modified Bessel function, you should look for BesselJ, not BesselI. However, I didn't find such in this lib.

alky
  • 1