bessel.c Source File

00001 #include <sadie.h>
00002 
00003 /*----------------------------------------------------------------------------*/
00004 /*-General Information--------------------------------------------------------*/
00005 /*                                                                            */
00006 /*   This function computes the modified Bessel Function of the first         */
00007 /*   kind and zero order. Taken from the Numerical Recipes book.              */
00008 /*                                                                            */
00009 /*----------------------------------------------------------------------------*/
00010 /*-Interface Information------------------------------------------------------*/
00011 double mBesselF (
00012 double x        /*  I   Argument to compute Io(x).                            */
00013 /*----------------------------------------------------------------------------*/
00014 ) {
00015     double ax, bessel;
00016     double y;
00017 
00018     if ((ax=fabs(x)) < 3.75) {
00019         y = x/3.75;
00020         y *= y;
00021         bessel = 1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
00022               +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
00023     } else {
00024         y = 3.75/ax;
00025         bessel = (exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
00026               +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
00027               +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
00028               +y*0.392377e-2))))))));
00029     }
00030 
00031     return (bessel);
00032 }