From 8b5c07f8055d6888459fd662c85f706ef84e0471 Mon Sep 17 00:00:00 2001 From: Franz Zotter Date: Sun, 11 Oct 2009 21:22:05 +0000 Subject: added [mtx_spherical_radial], the radial functions of the spherical base-solutions svn path=/trunk/externals/iem/iemmatrix/; revision=12577 --- src/iemmatrix_sources.c | 1 + src/iemmatrix_sources.h | 1 + src/mtx_spherical_harmonics/sph_radial.c | 82 ++++++++++++++ src/mtx_spherical_harmonics/sph_radial.h | 26 +++++ src/mtx_spherical_radial.c | 187 +++++++++++++++++++++++++++++++ 5 files changed, 297 insertions(+) create mode 100644 src/mtx_spherical_harmonics/sph_radial.c create mode 100644 src/mtx_spherical_harmonics/sph_radial.h create mode 100644 src/mtx_spherical_radial.c (limited to 'src') diff --git a/src/iemmatrix_sources.c b/src/iemmatrix_sources.c index a246ec1..5e05345 100644 --- a/src/iemmatrix_sources.c +++ b/src/iemmatrix_sources.c @@ -86,6 +86,7 @@ void iemmatrix_sources_setup(void) iemtx_sndfileread_setup(); /* mtx_sndfileread.c */ iemtx_sort_setup(); /* mtx_sort.c */ iemtx_spherical_harmonics_setup(); /* mtx_spherical_harmonics.c */ + iemtx_spherical_radial_setup(); /* mtx_spherical_radial.c */ iemtx_sub_setup(); /* mtx_sub.c */ iemtx_sum_setup(); /* mtx_sum.c */ iemtx_svd_setup(); /* mtx_svd.c */ diff --git a/src/iemmatrix_sources.h b/src/iemmatrix_sources.h index cda8fc5..3e56eba 100644 --- a/src/iemmatrix_sources.h +++ b/src/iemmatrix_sources.h @@ -84,6 +84,7 @@ void iemtx_slice_setup(void); /* mtx_slice.c */ void iemtx_sndfileread_setup(void); /* mtx_sndfileread.c */ void iemtx_sort_setup(void); /* mtx_sort.c */ void iemtx_spherical_harmonics_setup(void); /* mtx_spherical_harmonics.c */ +void iemtx_spherical_radial_setup(void); /* mtx_spherical_radial.c */ void iemtx_sub_setup(void); /* mtx_sub.c */ void iemtx_sum_setup(void); /* mtx_sum.c */ void iemtx_svd_setup(void); /* mtx_svd.c */ diff --git a/src/mtx_spherical_harmonics/sph_radial.c b/src/mtx_spherical_harmonics/sph_radial.c new file mode 100644 index 0000000..ec3536e --- /dev/null +++ b/src/mtx_spherical_harmonics/sph_radial.c @@ -0,0 +1,82 @@ +/* + * Recursive computation of (arbitrary degree) spherical Bessel/Neumann/Hankel functions, + * according to Gumerov and Duraiswami, + * "The Fast Multipole Methods for the Helmholtz Equation in Three Dimensions", + * Elsevier, 2005. + * + * Implementation by Franz Zotter, Institute of Electronic Music and Acoustics + * (IEM), University of Music and Dramatic Arts (KUG), Graz, Austria + * http://iem.at/Members/zotter, 2007. + * + * This code is published under the Gnu General Public License, see + * "LICENSE.txt" + */ + +#include + +#include "sph_radial.h" + +#define EPS 1e-10 + +static void radialRecurrence (double x, double *y, int n); + +// the two recurrences for higher n: +// by now no numeric stabilization for the bessel function is performed + +static void radialRecurrence (double x, double *y, int n) { + int k; + for (k=1;k=0) + y[0] = (x=1) + y[1] = -cos(x)/x + y[0]/x; + post("before %f %f %f",y[0],y[1],y[2]); + radialRecurrence (x,y,n); + post("after %f %f %f",y[0],y[1],y[2]); +} + +void sphNeumann (double x, double *y, int n) { + if (y==0) + return; + if (n>=0) + y[0] = -cos(x)/x; + if (n>=1) + y[1] = ((x +#include "mtx_spherical_harmonics/sph_radial.c" + +static t_class *mtx_spherical_radial_class; + +typedef struct _MTXSph_ MTXSph; +struct _MTXSph_ +{ + t_object x_obj; + t_outlet *list_h_re_out; + t_outlet *list_h_im_out; + t_atom *list_h_re; + t_atom *list_h_im; + double *kr; + double *h_re; + double *h_im; + size_t nmax; + size_t l; +}; + +static void allocMTXSphdata (MTXSph *x) +{ + x->kr=(double*)calloc(x->l,sizeof(double)); + if (x->list_h_re_out!=0) { + x->list_h_re=(t_atom*)calloc(x->l*(x->nmax+1)+2,sizeof(t_atom)); + x->h_re=(double*)calloc(x->l*(x->nmax+1),sizeof(double)); + } + if (x->list_h_im_out!=0) { + x->list_h_im=(t_atom*)calloc(x->l*(x->nmax+1)+2,sizeof(t_atom)); + x->h_im=(double*)calloc(x->l*(x->nmax+1),sizeof(double)); + } +} + +static void deleteMTXSphdata (MTXSph *x) +{ + if (x->kr!=0) + free(x->kr); + if (x->h_re!=0) + free(x->h_re); + if (x->h_im!=0) + free(x->h_im); + if (x->list_h_re!=0) + free(x->list_h_re); + if (x->list_h_im!=0) + free(x->list_h_im); + x->list_h_re=0; + x->list_h_im=0; + x->h_re=0; + x->h_im=0; + x->kr=0; +} + +static void *newMTXSph (t_symbol *s, int argc, t_atom *argv) +{ + int nmax; + char whichfunction = 'j'; + t_symbol *fsym; + MTXSph *x = (MTXSph *) pd_new (mtx_spherical_radial_class); + x->list_h_re = 0; + x->list_h_im = 0; + x->list_h_im_out = 0; + x->list_h_re_out = 0; + x->kr = 0; + x->h_re = 0; + x->h_im = 0; + x->l=0; + fsym=atom_getsymbol(argv); + if (fsym->s_name!=0) + whichfunction=fsym->s_name[0]; + post("%s",fsym->s_name); + switch (whichfunction) { + default: + case 'j': + post("j"); + x->list_h_re_out = outlet_new (&x->x_obj, gensym("matrix")); + break; + case 'h': + post("h"); + x->list_h_re_out = outlet_new (&x->x_obj, gensym("matrix")); + case 'y': + post("h or y"); + x->list_h_im_out = outlet_new (&x->x_obj, gensym("matrix")); + } + nmax=(int) atom_getfloat(argv+1); + if (nmax<0) + nmax=0; + x->nmax=nmax; + + return ((void *) x); +} + +static void mTXSphBang (MTXSph *x) +{ + if (x->list_h_im!=0) { + outlet_anything(x->list_h_im_out, gensym("matrix"), x->l*(x->nmax+1)+2, x->list_h_im); + } + if (x->list_h_re!=0) { + outlet_anything(x->list_h_re_out, gensym("matrix"), x->l*(x->nmax+1)+2, x->list_h_re); + } +} + +static void mTXSphMatrix (MTXSph *x, t_symbol *s, + int argc, t_atom *argv) +{ + int rows = atom_getint (argv++); + int columns = atom_getint (argv++); + int size = rows * columns; + int in_size = argc-2; + int n,ofs; + + + /* size check */ + if (!size) + post("mtx_spherical_radial: invalid dimensions"); + else if (in_sizel!=columns) { + deleteMTXSphdata(x); + x->l=columns; + allocMTXSphdata(x); + } + for (n=0;nl;n++) { + x->kr[n]=(double) atom_getfloat(argv+n); + } + + if (x->h_re!=0) + for (n=0,ofs=0;nl;n++,ofs+=x->nmax+1) + sphBessel(x->kr[n], x->h_re+ofs, x->nmax); + + if (x->h_im!=0) + for (n=0,ofs=0;nl;n++,ofs+=x->nmax+1) + sphNeumann(x->kr[n], x->h_im+ofs, x->nmax); + + if (x->h_re!=0) { + SETFLOAT(x->list_h_re+1,(float)(x->nmax+1)); + SETFLOAT(x->list_h_re,(float)x->l); + for (n=0;nl*(x->nmax+1);n++) + SETFLOAT(x->list_h_re+n+2,(float)x->h_re[n]); + } + + if (x->h_im!=0) { + SETFLOAT(x->list_h_im+1,(float)(x->nmax+1)); + SETFLOAT(x->list_h_im,(float)x->l); + for (n=0;nl*(x->nmax+1);n++) + SETFLOAT(x->list_h_im+n+2,(float)x->h_im[n]); + } + + mTXSphBang(x); + } +} + +void mtx_spherical_radial_setup (void) +{ + mtx_spherical_radial_class = class_new + (gensym("mtx_spherical_radial"), + (t_newmethod) newMTXSph, + (t_method) deleteMTXSphdata, + sizeof (MTXSph), + CLASS_DEFAULT, A_GIMME, 0); + class_addbang (mtx_spherical_radial_class, (t_method) mTXSphBang); + class_addmethod (mtx_spherical_radial_class, (t_method) mTXSphMatrix, gensym("matrix"), A_GIMME,0); +} + +void iemtx_spherical_radial_setup(void){ + mtx_spherical_radial_setup(); +} + + -- cgit v1.2.1