diff options
| -rw-r--r-- | src/config.h.in | 6 | ||||
| -rw-r--r-- | src/configure.ac | 5 | ||||
| -rw-r--r-- | src/mtx_bessel.c | 25 |
3 files changed, 29 insertions, 7 deletions
diff --git a/src/config.h.in b/src/config.h.in index ea56ead..c86dafa 100644 --- a/src/config.h.in +++ b/src/config.h.in @@ -9,9 +9,15 @@ /* can gsl compute complex eigenvalues? */ #undef HAVE_GSL_EIGEN_NONSYMM +/* can gsl compute Bessel functions? */ +#undef HAVE_GSL_BESSEL + /* Define to 1 if you have the <gsl/gsl_linalg.h> header file. */ #undef HAVE_GSL_GSL_LINALG_H +/* Define to 1 if you have the <gsl/gsl_sf_bessel.h> header file. */ +#undef HAVE_GSL_GSL_LINALG_H + /* Define to 1 if you have the <inttypes.h> header file. */ #undef HAVE_INTTYPES_H diff --git a/src/configure.ac b/src/configure.ac index b6285f9..2a0a7be 100644 --- a/src/configure.ac +++ b/src/configure.ac @@ -69,9 +69,10 @@ dnl AC_CHECK_LIB([m], [cos]) AC_CHECK_LIB([gslcblas], [cblas_dgemm]) AC_CHECK_LIB([gsl], [gsl_blas_dgemm]) AC_CHECK_LIB([gsl], [gsl_eigen_nonsymm],AC_DEFINE(HAVE_GSL_EIGEN_NONSYMM,1,[can gsl compute complex eigenvalues?])) +AC_CHECK_LIB([gsl], [gsl_sf_bessel_Jn],AC_DEFINE(HAVE_GSL_BESSEL,1,[can gsl compute Bessel functions?])) dnl for math.h Bessel/Neumann functions -AC_CHECK_LIB([m], [jnf],AC_DEFINE(HAVE_MATH_BESSEL,1,[can math compute Bessel functions?])) +AC_CHECK_LIB([m], [jn],AC_DEFINE(HAVE_MATH_BESSEL,1,[can math compute Bessel functions?])) if test "x$with_pd" != "x"; then @@ -100,7 +101,7 @@ AC_CHECK_LIB(pd, nullfn) dnl Checks for header files. AC_HEADER_STDC -AC_CHECK_HEADERS(stdlib.h stdio.h string.h math.h time.h sys/time.h fftw3.h sndfile.h gsl/gsl_linalg.h) +AC_CHECK_HEADERS(stdlib.h stdio.h string.h math.h time.h sys/time.h fftw3.h sndfile.h gsl/gsl_linalg.h gsl/gsl_sf_bessel.h) dnl Checks for typedefs, structures, and compiler characteristics. AC_HEADER_TIME diff --git a/src/mtx_bessel.c b/src/mtx_bessel.c index f17342f..458100c 100644 --- a/src/mtx_bessel.c +++ b/src/mtx_bessel.c @@ -16,6 +16,9 @@ #include "iemmatrix.h" #include "math.h" #include <stdlib.h> +#ifdef HAVE_GSL_BESSEL +#include <gsl/gsl_sf_bessel.h> +#endif static t_class *mtx_bessel_class; @@ -121,7 +124,7 @@ static void mTXBesselMatrix (MTXBessel *x, t_symbol *s, int n,m,ofs; -#ifdef HAVE_MATH_BESSEL +#if defined HAVE_MATH_BESSEL || defined HAVE_GSL_BESSEL /* size check */ if (!size) @@ -139,16 +142,28 @@ static void mTXBesselMatrix (MTXBessel *x, t_symbol *s, for (n=0;n<x->l;n++) { x->kr[n]=(double) atom_getfloat(argv+n); } - + +#ifdef HAVE_GSL_BESSEL if (x->h_re!=0) for (m=0;m<x->l;m++) for (n=0;n<x->nmax+1;n++) - x->h_re[n+m*(x->nmax+1)]=jnf(n,x->kr[m]); + x->h_re[n+m*(x->nmax+1)]=gsl_sf_bessel_Jn(n,x->kr[m]); if (x->h_im!=0) for (m=0;m<x->l;m++) for (n=0;n<x->nmax+1;n++) - x->h_im[n+m*(x->nmax+1)]=ynf(n,x->kr[m]); + x->h_im[n+m*(x->nmax+1)]=gsl_sf_bessel_Yn(n,x->kr[m]); +#else + if (x->h_re!=0) + for (m=0;m<x->l;m++) + for (n=0;n<x->nmax+1;n++) + x->h_re[n+m*(x->nmax+1)]=jn(n,x->kr[m]); + + if (x->h_im!=0) + for (m=0;m<x->l;m++) + for (n=0;n<x->nmax+1;n++) + x->h_im[n+m*(x->nmax+1)]=yn(n,x->kr[m]); +#endif if (x->h_re!=0) { @@ -169,7 +184,7 @@ static void mTXBesselMatrix (MTXBessel *x, t_symbol *s, } #else - post("mtx_bessel: implementation requires math.h that implements jnf and ynf Bessel functions"); + post("mtx_bessel: implementation requires math.h that implements jn and yn Bessel functions or gsl_sf_bessel.h"); #endif } |