/* This file contains the routines to calculate 
   modified Bessel functions. */

#include <bessel_mod.h>

/* The function bessI_0(x) returns the modified bessel function I_0(x) 
   for positive real arguments x, 
   using the polynomial approximations of Abramowitz and Stegun, p.378. */

double bessI_0(double x)
{
        double ax, ans;
        double y;

        if ((ax=fabs(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=(exp(ax)/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;
}


/* The function bessK_0(x) returns the modified bessel function K_0(x) 
   for positive real arguments x, 
   using the polynomial approximations of Abramowitz and Stegun, p.378. */

double bessK_0(double x)
{
        double y,ans;

        if (x <= 2.0) {
                y=x*x/4.0;
                ans=(-log(x/2.0)*bessI_0(x))+(-0.57721566+y*(0.42278420
                        +y*(0.23069756+y*(0.3488590e-1+y*(0.262698e-2
                        +y*(0.10750e-3+y*0.74e-5))))));
        } else {
                y=2.0/x;
                ans=(exp(-x)/sqrt(x))*(1.25331414+y*(-0.7832358e-1
                        +y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2
                        +y*(-0.251540e-2+y*0.53208e-3))))));
        }
        return ans;
}