aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/iemmatrix_sources.c1
-rw-r--r--src/iemmatrix_sources.h1
-rw-r--r--src/mtx_spherical_harmonics/sph_radial.c82
-rw-r--r--src/mtx_spherical_harmonics/sph_radial.h26
-rw-r--r--src/mtx_spherical_radial.c187
5 files changed, 297 insertions, 0 deletions
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 <math.h>
+
+#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<n;k++) {
+ y[k+1] = -y[k-1] + y[k]/x * (2*k+1);
+ }
+}
+
+static void radialDiffRecurrence (double x, double *y1, double *yd, int n) {
+ int k;
+ for (k=0;k<n;k++) {
+ yd[k] = y1[k]/x * n - y1[k+1];
+ }
+}
+
+void sphBessel (double x, double *y, int n) { //TODO: small values!
+ if (y==0)
+ return;
+ if (n>=0)
+ y[0] = (x<EPS)?1.0:sin(x)/x;
+ if (n>=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<EPS)?1.0:sin(x)/x) - y[0]/x;
+ radialRecurrence (x,y,n);
+}
+
+void sphBesselDiff (double x, double *y, int n) {
+ double *y1;
+ if (n<0)
+ return;
+ if ((y1 = (double*)calloc(n+2,sizeof(double)))==0)
+ return;
+ sphBessel (x,y1,n+1);
+ radialDiffRecurrence (x,y1,y,n);
+ free(y1);
+}
+
+void sphNeumannDiff (double x, double *y, int n) {
+ double *y1;
+ if (n<0)
+ return;
+ if ((y1 = (double*)calloc(n+2,sizeof(double)))==0)
+ return;
+ sphNeumann (x,y,n+1);
+ radialDiffRecurrence (x,y1,y,n);
+ free(y1);
+}
diff --git a/src/mtx_spherical_harmonics/sph_radial.h b/src/mtx_spherical_harmonics/sph_radial.h
new file mode 100644
index 0000000..43661ce
--- /dev/null
+++ b/src/mtx_spherical_harmonics/sph_radial.h
@@ -0,0 +1,26 @@
+/*
+ * 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"
+ */
+
+#ifndef __SPH_RADIAL_H__
+#define __SPH_RADIAL_H__
+
+void sphBessel (double x, double *y, int n);
+
+void sphNeumann (double x, double *y, int n);
+
+void sphBesselDiff (double x, double *y, int n);
+
+void sphNeumannDiff (double x, double *y, int n);
+
+#endif // __SPH_RADIAL_H__
diff --git a/src/mtx_spherical_radial.c b/src/mtx_spherical_radial.c
new file mode 100644
index 0000000..5554a6f
--- /dev/null
+++ b/src/mtx_spherical_radial.c
@@ -0,0 +1,187 @@
+/*
+ * iemmatrix
+ *
+ * objects for manipulating simple matrices
+ * mostly refering to matlab/octave matrix functions
+ * this functions depends on the GNU scientific library
+ *
+ * Copyright (c) 2009, Franz Zotter
+ * IEM, Graz, Austria
+ *
+ * For information on usage and redistribution, and for a DISCLAIMER OF ALL
+ * WARRANTIES, see the file, "LICENSE.txt," in this distribution.
+ *
+ */
+
+#include "iemmatrix.h"
+#include <stdlib.h>
+#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_size<size)
+ post("mtx_spherical_radial: sparse matrix not yet supported: use \"mtx_check\"");
+ else if ((rows!=1)||(columns<1))
+ post("mtx_spherical_radial: 1 X L matrix expected with kr and h vector, but got more rows/no entries");
+ else {
+ if (x->l!=columns) {
+ deleteMTXSphdata(x);
+ x->l=columns;
+ allocMTXSphdata(x);
+ }
+ for (n=0;n<x->l;n++) {
+ x->kr[n]=(double) atom_getfloat(argv+n);
+ }
+
+ if (x->h_re!=0)
+ for (n=0,ofs=0;n<x->l;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;n<x->l;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;n<x->l*(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;n<x->l*(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();
+}
+
+
generated by cgit v1.2.1 (git 2.18.0) at 2024-04-28 04:32:23 +0000