aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--doc/mtx_spherical_radial-help.pd241
-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
6 files changed, 538 insertions, 0 deletions
diff --git a/doc/mtx_spherical_radial-help.pd b/doc/mtx_spherical_radial-help.pd
new file mode 100644
index 0000000..5b6afe6
--- /dev/null
+++ b/doc/mtx_spherical_radial-help.pd
@@ -0,0 +1,241 @@
+#N canvas 379 0 1021 852 10;
+#N canvas 0 0 450 300 (subpatch) 0;
+#X array j0 100 float 1;
+#A 0 0.998334 0.964227 0.888657 0.777588 0.639683 0.485483 0.32644
+0.173872 0.0379576 -0.0731249 -0.153847 -0.201696 -0.217232 -0.203847
+-0.167252 -0.114765 -0.054476 0.00562416 0.058443 0.0984482 0.122122
+0.128169 0.117474 0.0928211 0.0584223 0.0193205 -0.0192641 -0.0525636
+-0.0767977 -0.0895608 -0.0900205 -0.0789194 -0.0583864 -0.0315936 -0.00230552
+0.0256184 0.0487061 0.0642826 0.0707609 0.067784 0.0562077 0.0379347
+0.0156251 -0.00767335 -0.0289267 -0.045498 -0.0554624 -0.0578212 -0.0525892
+-0.0407518 -0.0240998 -0.00497058 0.0140721 0.0305726 0.0424939 0.0484624
+0.0479197 0.0411655 0.0292891 0.0140031 -0.002598 -0.018319 -0.03115
+-0.0395184 -0.0424783 -0.0398125 -0.0320383 -0.0203186 -0.0062929 0.00814866
+0.0211162 0.0309648 0.0365 0.03712 0.0328773 0.0244546 0.013059 0.000252624
+-0.0122595 -0.0228541 -0.0301981 -0.0334145 -0.0321848 -0.0267766 -0.0179956
+-0.00706707 0.00453248 0.0152732 0.0237724 0.0289707 0.0302626 0.0275654
+0.0213185 0.0124147 0.00207441 -0.00832043 -0.0174098 -0.0240321 -0.0273718
+-0.0270586;
+#X coords 0 1 99 -1 200 140 1;
+#X restore 69 219 graph;
+#N canvas 0 0 450 300 (subpatch) 0;
+#X array y0 100 float 1;
+#A 0 -9.95004 -1.91809 -0.810089 -0.304619 -0.00483966 0.181763 0.288657
+0.333212 0.328231 0.285574 0.21711 0.13461 0.0491928 -0.0293734 -0.0933163
+-0.137346 -0.158906 -0.15816 -0.137711 -0.102106 -0.0571689 -0.00924676
+0.0355409 0.0719789 0.0962818 0.106425 0.102249 0.0853441 0.0587414
+0.0264484 -0.00709537 -0.037585 -0.0613658 -0.0758406 -0.0797276 -0.0731431
+-0.0575093 -0.0353085 -0.00971747 0.0158235 0.0380541 0.0542939 0.0627506
+0.0627045 0.0545503 0.0396971 0.0203463 -0.000821204 -0.0210054 -0.037647
+-0.048744 -0.0530835 -0.0503669 -0.0412123 -0.027042 -0.00987069 0.00797327
+0.0241575 0.0366444 0.0439426 0.0452802 0.04068 0.0309308 0.0174615
+0.00213646 -0.0129974 -0.0259804 -0.0351907 -0.0395426 -0.0386117 -0.0326705
+-0.0226335 -0.00992048 0.00374029 0.016545 0.0268491 0.0333767 0.0353758
+0.0327023 0.0258234 0.0157413 0.00384979 -0.00825582 -0.0189914 -0.0269888
+-0.0312691 -0.0313604 -0.0273472 -0.0198471 -0.00991801 0.00108892
+0.0117122 0.0205722 0.026549 0.028923 0.0274601 0.0224315 0.0145684
+0.00495724 -0.00510983;
+#X coords 0 1 99 -1 200 140 1;
+#X restore 273 218 graph;
+#N canvas 0 0 450 300 (subpatch) 0;
+#X array j1 100 float 1;
+#A 0 0.0333 0.151926 0.258502 0.344765 0.404366 0.433434 0.430906 0.398561
+0.340773 0.264018 0.176172 0.0857012 0.000808178 -0.0713561 -0.125349
+-0.157887 -0.168057 -0.15727 -0.128968 -0.0881426 -0.0407019 0.00722327
+0.0499589 0.0828816 0.102861 0.108515 0.100244 0.0800755 0.051316 0.0180848
+-0.0152242 -0.0444836 -0.0663113 -0.0784363 -0.0799115 -0.0711576 -0.0538387
+-0.0305939 -0.00466337 0.0205417 0.0418694 0.0568064 0.063761 0.0622198
+0.0527642 0.0369499 0.0170698 -0.00416483 -0.0239835 -0.0399079 -0.0500544
+-0.0533485 -0.049631 -0.0396435 -0.0249017 -0.00747387 0.0103011 0.0261224
+0.0380184 0.0445884 0.0451623 0.0398627 0.0295634 0.0157541 0.000329767
+-0.0146647 -0.027302 -0.0360164 -0.0397945 -0.0382901 -0.0318491 -0.0214458
+-0.0085399 0.00512516 0.0177551 0.0277372 0.0338448 0.0353847 0.0322742
+0.0250353 0.0147129 0.00272587 -0.00932521 -0.0198704 -0.0275726 -0.0314957
+-0.0312168 -0.0268688 -0.0191109 -0.00903089 0.00200534 0.0125378 0.0212037
+0.0269128 0.0289832 0.0272214 0.0219371 0.0138931 0.00419583 -0.00585494
+;
+#X coords 0 1 99 -1 200 140 1;
+#X restore 69 379 graph;
+#N canvas 0 0 450 300 (subpatch) 0;
+#X array y1 100 float 1;
+#A 0 100.499 5.08199 1.86277 1.03198 0.642778 0.391258 0.200655 0.0486347
+-0.0704963 -0.157309 -0.211619 -0.234338 -0.228189 -0.197797 -0.149379
+-0.0901827 -0.0277821 0.0306545 0.0790444 0.112931 0.129831 0.129357
+0.113112 0.0843665 0.047579 0.00780913 -0.0299028 -0.0611178 -0.0824773
+-0.0920306 -0.0893798 -0.075634 -0.0531885 -0.0253627 0.00405364 0.031287
+0.0530402 0.0668722 0.071455 0.0666825 0.0536247 0.0343387 0.0115672
+-0.0116345 -0.032295 -0.047895 -0.0566644 -0.0577737 -0.0513997 -0.0386632
+-0.0214493 -0.00214039 0.016706 0.0326874 0.043856 0.0489506 0.0475324
+0.0400125 0.0275701 0.0119765 -0.00465167 -0.0201339 -0.0325078 -0.0402728
+-0.0425692 -0.0392682 -0.0309667 -0.0188886 -0.00470961 0.00967236
+0.0223871 0.0318329 0.0368752 0.0369804 0.0322683 0.0234795 0.0118627
+-0.000998855 -0.0134016 -0.0237446 -0.0307342 -0.033544 -0.0319104
+-0.0261532 -0.0171201 -0.00606465 0.00552617 0.0161298 0.0243871 0.0292745
+0.0302296 0.0272147 0.020709 0.0116366 0.00123572 -0.00910834 -0.0180468
+-0.0244415 -0.0275097 -0.0269179;
+#X coords 0 1 99 -1 200 140 1;
+#X restore 273 378 graph;
+#N canvas 0 0 450 300 (subpatch) 0;
+#X array j2 100 float 1;
+#A 0 0.000666191 0.0142423 0.0438713 0.0861799 0.136337 0.188585 0.236874
+0.275523 0.299836 0.306613 0.294482 0.264041 0.217772 0.159759 0.0952301
+0.0299879 -0.0302173 -0.0802925 -0.116323 -0.135954 -0.138587 -0.125385
+-0.0990795 -0.0636155 -0.0236694 0.0158917 0.0505546 0.0766421 0.0916826
+0.0946273 0.0858963 0.0672542 0.041536 0.0122611 -0.0168159 -0.0421625
+-0.0608785 -0.071014 -0.0717602 -0.0634945 -0.0476817 -0.0266473 -0.00325543
+0.019465 0.0387006 0.0521913 0.0584877 0.0570986 0.0485147 0.0341096
+0.0159345 -0.00356237 -0.0218585 -0.0366754 -0.0462567 -0.0495713 -0.0464185
+-0.037425 -0.0239387 -0.00783391 0.00874297 0.0236544 0.0350433 0.0415603
+0.0425204 0.0379701 0.0286599 0.015928 0.00151276 -0.0126817 -0.024833
+-0.0334324 -0.037469 -0.0365463 -0.0309168 -0.0214326 -0.00941995 0.00350276
+0.015641 0.0254441 0.0317013 0.0336896 0.0312552 0.0248197 0.0153124
+0.00403802 -0.00749991 -0.0177981 -0.0255479 -0.0298004 -0.0300804
+-0.0264389 -0.0194339 -0.0100484 0.000446893 0.0106636 0.0192785 0.0252034
+0.0277219 0.0265749;
+#X coords 0 1 99 -1 200 140 1;
+#X restore 69 549 graph;
+#N canvas 0 0 450 300 (subpatch) 0;
+#X array y2 100 float 1;
+#A 0 3024.91 34.6483 7.52992 2.89013 1.2384 0.426713 -0.0263457 -0.278374
+-0.398111 -0.424692 -0.386041 -0.305084 -0.201668 -0.0928367 0.00748753
+0.0889225 0.144905 0.172714 0.173186 0.150159 0.109688 0.0591149 0.0061069
+-0.0422499 -0.0802067 -0.103891 -0.111582 -0.103722 -0.0826649 -0.0522311
+-0.0171175 0.017751 0.0478499 0.0695893 0.0806976 0.0804173 0.0695012
+0.050022 0.0250285 -0.00189889 -0.0271343 -0.0474708 -0.0605065 -0.0649094
+-0.0605325 -0.0483731 -0.0303889 -0.00920143 0.0122733 0.031212 0.045245
+0.0527412 0.0529879 0.0462443 0.0336689 0.0171336 -0.00104611 -0.0184281
+-0.0327644 -0.0422855 -0.0459131 -0.0433748 -0.0352119 -0.0226813 -0.00756812
+0.00806366 0.0221484 0.0328881 0.0389768 0.0397568 0.0352832 0.0262963
+0.0141048 0.000398699 -0.012982 -0.024291 -0.0321012 -0.0354818 -0.0341065
+-0.0282798 -0.0188812 -0.00723461 0.00507498 0.0164157 0.0253228 0.0306859
+0.0318857 0.0288629 0.0221128 0.0126073 0.00165734 -0.00926688 -0.0187316
+-0.0255258 -0.0288155 -0.0282441 -0.0239688 -0.0166291 -0.00725295
+0.00288612;
+#X coords 0 1 99 -1 200 140 1;
+#X restore 273 548 graph;
+#N canvas 0 0 450 300 plot 0;
+#X obj 60 19 inlet;
+#X obj 44 198 mtx;
+#X obj 58 42 t a a;
+#X obj 39 70 mtx_size;
+#X obj 39 112 until;
+#X obj 39 92 t f b;
+#X obj 68 131 + 1;
+#X obj 39 131 i;
+#X msg 82 93 1;
+#X obj 43 152 t f f;
+#X obj 152 243 outlet;
+#X obj 44 222 list prepend 0;
+#X obj 45 243 outlet;
+#X obj 152 219 - 1;
+#X msg 44 175 col \$1;
+#X connect 0 0 2 0;
+#X connect 1 0 11 0;
+#X connect 2 0 3 0;
+#X connect 2 1 1 1;
+#X connect 3 1 5 0;
+#X connect 4 0 7 0;
+#X connect 5 0 4 0;
+#X connect 5 1 8 0;
+#X connect 6 0 7 1;
+#X connect 7 0 6 0;
+#X connect 7 0 9 0;
+#X connect 8 0 7 1;
+#X connect 9 0 14 0;
+#X connect 9 1 13 0;
+#X connect 11 0 12 0;
+#X connect 13 0 10 0;
+#X connect 14 0 1 0;
+#X restore 124 112 pd plot;
+#X obj 124 160 s;
+#X msg 162 135 symbol j\$1;
+#N canvas 0 0 450 300 plot 0;
+#X obj 60 19 inlet;
+#X obj 44 198 mtx;
+#X obj 58 42 t a a;
+#X obj 39 70 mtx_size;
+#X obj 39 112 until;
+#X obj 39 92 t f b;
+#X obj 68 131 + 1;
+#X obj 39 131 i;
+#X msg 82 93 1;
+#X obj 43 152 t f f;
+#X obj 152 243 outlet;
+#X obj 44 222 list prepend 0;
+#X obj 45 243 outlet;
+#X obj 152 219 - 1;
+#X msg 44 175 col \$1;
+#X connect 0 0 2 0;
+#X connect 1 0 11 0;
+#X connect 2 0 3 0;
+#X connect 2 1 1 1;
+#X connect 3 1 5 0;
+#X connect 4 0 7 0;
+#X connect 5 0 4 0;
+#X connect 5 1 8 0;
+#X connect 6 0 7 1;
+#X connect 7 0 6 0;
+#X connect 7 0 9 0;
+#X connect 8 0 7 1;
+#X connect 9 0 14 0;
+#X connect 9 1 13 0;
+#X connect 11 0 12 0;
+#X connect 13 0 10 0;
+#X connect 14 0 1 0;
+#X restore 245 114 pd plot;
+#X obj 245 162 s;
+#X msg 283 137 symbol y\$1;
+#N canvas 0 0 450 300 (subpatch) 0;
+#X array j3 100 float 1;
+#A 0 9.51852e-06 0.00095102 0.00527028 0.0150905 0.0317082 0.0553716
+0.085194 0.119212 0.154584 0.187911 0.215627 0.234435 0.241716 0.235869
+0.216543 0.184724 0.142677 0.0937343 0.0419586 -0.00827382 -0.0527338
+-0.0877843 -0.110761 -0.120243 -0.11619 -0.09992 -0.0739438 -0.0416654
+-0.00699329 0.0260988 0.0540063 0.0738778 0.0839026 0.083473 0.0732052
+0.0548197 0.0308987 0.00455261 -0.020964 -0.0426398 -0.0580521 -0.0656311
+-0.0648136 -0.0560715 -0.0408163 -0.0211929 0.000206524 0.020674 0.0377202
+0.0493699 0.0543867 0.0523989 0.0439155 0.0302337 0.0132523 -0.00478447
+-0.0215758 -0.0350539 -0.0436334 -0.0463949 -0.0431797 -0.0345861 -0.0218717
+-0.00677619 0.00871263 0.0226157 0.0332128 0.0392526 0.0400974 0.0357879
+0.027019 0.0150343 0.00145374 -0.0119425 -0.0234446 -0.031629 -0.0355328
+-0.0347652 -0.0295429 -0.0206482 -0.00931503 0.00293997 0.0145178 0.0239443
+0.0300561 0.0321429 0.0300285 0.0240813 0.015155 0.00446824 -0.00655985
+-0.016497 -0.0240825 -0.0283853 -0.0289184 -0.0256915 -0.0192001 -0.0103516
+-0.000340117 0.00951389;
+#X coords 0 1 99 -1 200 140 1;
+#X restore 69 709 graph;
+#N canvas 0 0 450 300 (subpatch) 0;
+#X array y3 100 float 1;
+#A 0 151145 366.834 43.4101 11.0361 3.31824 0.714766 -0.258057 -0.571765
+-0.587221 -0.468662 -0.301997 -0.135561 0.0035989 0.102198 0.15655
+0.169761 0.149492 0.106014 0.0504976 -0.00643979 -0.0558783 -0.0913754
+-0.109365 -0.10918 -0.0927436 -0.0639955 -0.0281465 0.0091362 0.042514
+0.0676428 0.0816513 0.0833923 0.0734539 0.0539494 0.028129 -0.000125397
+-0.0268511 -0.0485288 -0.0625167 -0.0673434 -0.0628338 -0.0500594 -0.031131
+-0.00886785 0.013607 0.0332908 0.047688 0.0551132 0.0548748 0.0473213
+0.0337506 0.0162001 -0.00285082 -0.0208224 -0.0353766 -0.0447137 -0.0477865
+-0.0444104 -0.0352552 -0.0217275 -0.00576019 0.0104582 0.0247791 0.0353731
+0.0409598 0.0409568 0.0355346 0.0255707 0.0125128 -0.00182793 -0.0155244
+-0.0267899 -0.0342077 -0.036906 -0.0346574 -0.0278903 -0.0176154 -0.00527729
+0.00744578 0.0188686 0.0275192 0.0323273 0.0327536 0.0288477 0.0212272
+0.0109833 -0.00047449 -0.0116095 -0.0209631 -0.0273442 -0.0299787 -0.0286024
+-0.0234837 -0.0153772 -0.00541359 0.0050563 0.0146438 0.0221049 0.0265009
+0.0273153;
+#X coords 0 1 99 -1 200 140 1;
+#X restore 273 708 graph;
+#X obj 365 54 hsl 128 15 1 300 0 0 empty empty empty -2 -8 0 10 -262144
+-1 -1 1500 1;
+#X msg 124 42 0.1 \$1 100;
+#X obj 125 87 mtx_spherical_radial h 3;
+#X obj 125 63 mtx_linspace;
+#X connect 6 0 7 0;
+#X connect 6 1 8 0;
+#X connect 8 0 7 1;
+#X connect 9 0 10 0;
+#X connect 9 1 11 0;
+#X connect 11 0 10 1;
+#X connect 14 0 15 0;
+#X connect 15 0 17 0;
+#X connect 16 0 6 0;
+#X connect 16 1 9 0;
+#X connect 17 0 16 0;
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-12-27 10:57:11 +0000