public class Bessel
{
  public static double modBessI0(double x)
  {
    double ax, ans;
    double y;
    if ((ax = Math.abs(x)) < 3.75)
    {
      y = x/3.75;
      y *= y;
      ans = 1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
       +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
    }
    else
    {
      y=3.75/ax;
      ans=(Math.exp(ax)/Math.sqrt(ax))*(0.39894228+y*(0.1328592e-1
      +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
      +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
      +y*0.392377e-2))))))));  
    }
    return ans;
  }

  public static double modBessI1(double x)
  {
    double ax,ans;
    double y; 
    if ((ax=Math.abs(x)) < 3.75) 
    { 
      y=x/3.75;
      y*=y;
      ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
      +y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
    } 
    else 
    {
      y=3.75/ax;
      ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1
      -y*0.420059e-2));
      ans=0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2
      +y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
      ans *= (Math.exp(ax)/Math.sqrt(ax));
    }
    return x < 0.0 ? -ans : ans;
  }

  public static double  modBessI(int n, double x)
  {
    final double ACC = 40.0;
    final double  BIGNO = 1.0e10;
    final double  BIGNI = 1.0e-10;

    if (n == 0)
      return modBessI0(x);
    if (n == 1)
      return modBessI1(x);

    int j;
    double bi,bim,bip,tox,ans;
    if (x == 0.0)
      return 0.0;
    else 
    {
      tox=2.0/Math.abs(x);
      bip=ans=0.0;
      bi=1.0;
      for (j=2*(n+(int) Math.sqrt(ACC*n));j>0;j--) 
      {
        bim=bip+j*tox*bi;
        bip=bi;
        bi=bim;
        if (Math.abs(bi) > BIGNO) 
        {
          ans *= BIGNI;
          bi *= BIGNI;
          bip *= BIGNI;
        }
        if (j == n) 
          ans=bip;
      }
      ans *= modBessI0(x)/bi; 
      return ( (x < 0.0) && ((n & 1) == 1) ) ? -ans : ans;
    }
  } 
}