[mtx_bessel] extended to gsl/gsl_sf_bessel.h support

svn path=/trunk/externals/iem/iemmatrix/; revision=14939

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
 }