/* This file contains the functions to calculate */
/* spherical Bessel functions. */

#include <bessel.h>

double bessel_spher(int l, double x)
{
     double s = 0.0, s0, s1, s_limit, l_d, term1, term2, double_fact;
     double *bes;
     int i, l_it;
     
     s0 = sin(x) / x;
     s1 = ( sin(x) / pow(x,2.) - (cos(x) / x) );
     
     if (l == 0)
        s = s0;
     if (l == 1)
        s = s1;

     if (l > 1)
     {
     bes = (double *) malloc( (size_t) ((l + 1) * sizeof(double)) );
       
     bes[0] = s0;
     bes[1] = s1;
     
     for(l_it = 2; l_it <= l; l_it++)
     {
     double_fact = 1.0;
     l_d = (double) l_it; 
     for (i = 1; i <= (2 * l_it + 1); i += 2)
        double_fact *= i;  
     s_limit = pow(x, l_d) / double_fact; 

     if ( s_limit < pow(10., - 8.) )
        bes[l_it] = s_limit;
     
     else
     { 
       term1 = (2.0 * l_d - 1.0) / x * bes[l_it - 1];
       term2 = bes[l_it - 2];
       bes[l_it] = term1 - term2; 
     }
     }
     s = bes[l];
     free( (FREE_ARG) (bes));
     
     } 
     
     return s;
}