diff options
author | Franz Zotter <fzotter@users.sourceforge.net> | 2009-01-14 10:58:32 +0000 |
---|---|---|
committer | Franz Zotter <fzotter@users.sourceforge.net> | 2009-01-14 10:58:32 +0000 |
commit | 92b9deaf8d7c2a96c3977e1204341d6b6feb1fc1 (patch) | |
tree | 6e17422343f2a63bba88422d02e5d06fdf8d7baf /src/mtx_sh/sharmonics.c | |
parent | 96987170c95a18b779cd2bfc316d88d754db4b8e (diff) |
renamed [mtx_sh] to [mtx_spherical_harmonics].
svn path=/trunk/externals/iem/iemmatrix/; revision=10549
Diffstat (limited to 'src/mtx_sh/sharmonics.c')
-rw-r--r-- | src/mtx_sh/sharmonics.c | 157 |
1 files changed, 0 insertions, 157 deletions
diff --git a/src/mtx_sh/sharmonics.c b/src/mtx_sh/sharmonics.c deleted file mode 100644 index 8593182..0000000 --- a/src/mtx_sh/sharmonics.c +++ /dev/null @@ -1,157 +0,0 @@ -/* - * Recursive computation of (arbitrary degree) spherical harmonics, - * 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, 2008. - * - * This code is published under the Gnu General Public License, see - * "LICENSE.txt" - */ - -#include "mtx_sh/sharmonics.h" - -// HELPER ROUTINES - -// preliminarily writing normalized Legendre functions into the result -// Y_n^m(theta) = N_n^m * P_n^m(cos(theta)) -// ny0 and np0 denote where the position (n,m)=(n,0) is in the arrays -// ly0 and lp0 denote the starting position for one vertex in the arrays -// see below to find out how the data is arranged -static void sharmonics_initlegendrenormlzd(SHWorkSpace *ws) { - int n,m,ny0,np0; - int l,ly0,lp0; - const int pincr=(ws->nmax+1)*(ws->nmax+2)/2; - const int yincr=(ws->nmax+1)*(ws->nmax+1); - - for (n=0,ny0=0,np0=0; n<=ws->nmax; n++) { - for (m=0; m<=n; m++) { - ly0=0; - lp0=0; - for (l=0; l<ws->l; l++) { - ws->y[ly0+ny0+m] = ws->wn->n[np0+m] * ws->wl->p[lp0+np0+m]; - ws->y[ly0+ny0-m] = ws->y[ly0+ny0+m]; - ly0+=yincr; - lp0+=pincr; - } - } - ny0+=2*n+2; - np0+=n+1; - } -} - -// multiplying Chebyshev sin/cos to the preliminary result -// Y_n^m(phi,theta) = Y_n^m(theta) * T_m(phi) -// ny0 and nt0 denote where the position (n,m)=(n,0) or m=0 is in the arrays -// ly0 and lt0 denote the starting position for one vertex in the arrays -// see below to find out how the data is arranged -static void sharmonics_multcheby12(SHWorkSpace *ws) { - int n,m,ny0; - const int nt0=ws->nmax; - int l,ly0,lt0; - const int tincr=2*ws->nmax+1; - const int yincr=(ws->nmax+1)*(ws->nmax+1); - - for (n=0,ny0=0; n<=ws->nmax; n++) { - m=0; - ly0=0; - lt0=nt0; - for (l=0; l<ws->l; l++) { - ws->y[ly0+ny0+m]*= ws->wc->t[lt0+m]; - ly0+=yincr; - lt0+=tincr; - } - for (m=1; m<=n; m++) { - ly0=0; - lt0=nt0; - for (l=0; l<ws->l; l++) { - ws->y[ly0+ny0-m]*= -ws->wc->t[lt0-m]; - ws->y[ly0+ny0+m]*= ws->wc->t[lt0+m]; - ly0+=yincr; - lt0+=tincr; - } - } - ny0+=2*n+2; - } -} - - -/* MAIN PROGRAM. IMPORTANT EXPRESSIONS - - p... vector containing Legendre functions evaluated at the vector z=cos(theta) - structure [P_0^0(z1) P_1^0(z1) P_1^1(z1) P_2^0(z1) .... Pnmax^nmax(z1) - P_0^0(z2) P_1^0(z1) P_1^1(z2) P_2^0(z2) .... Pnmax^nmax(z2) - ... - P_0^0(zL) P_1^0(zL) P_1^1(zL) P_2^0(zL) .... Pnmax^nmax(zL)] - with length L X (nmax+1)*(nmax+2)/2 - - t... vector containing Chebyshev polynomials sin/cos evaluated at the vector phi - structure [T_-nmax(phi1) ... T_-1(phi1) T_0(phi1) T_1(phi1) ... T_nmax(phi1) - T_-nmax(phi2) ... T_-1(phi2) T_0(phi2) T_1(phi2) ... T_nmax(phi2) - ... - T_-nmax(phiL) ... T_-1(phiL) T_0(phiL) T_1(phiL) ... T_nmax(phiL)] - with length L X (2*nmax+1); negative indices are sine, positive ones - cosine terms - - norml ... vector containing normalization terms - structure [N_0^0 N_1^0 N_1^1 N_2^0 N_2^1 N_2^2 .... N_nmax^nmax] - with length (nmax+1)*(nmax+2)/2 - - - y ... THE RESULT: containing the spherical harmonics, with negative m for sine - positive m for cosine terms; p=(phi,theta) - structure [Y_0^0(p1) Y_1^-1(p1) Y_1^0(p1) Y_1^1 ... Y_nmax^nmax(p1) - Y_0^0(p2) Y_1^-1(p2) Y_1^0(p2) Y_1^1 ... Y_nmax^nmax(p2) - ... - Y_0^0(pL) Y_1^-1(pL) Y_1^0(pL) Y_1^1 ... Y_nmax^nmax(pL)] - with length L X (nmax+1)^2 - -*/ - -SHWorkSpace *sharmonics_alloc(size_t nmax, size_t l) { - SHWorkSpace *ws=0; - - if ((ws=(SHWorkSpace*)calloc(1,sizeof(SHWorkSpace)))!=0) { - ws->y=(double*)calloc(l*(nmax+1)*(nmax+1),sizeof(double)); - - ws->wl=(LegendreWorkSpace*)legendre_a_alloc(nmax,l); - ws->wc=(Cheby12WorkSpace*)chebyshev12_alloc(nmax,l); - ws->wn=(SHNorml*)sharmonics_normalization_new(nmax); - - if ((ws->y==0)||(ws->wl==0)||(ws->wc==0)||(ws->wn==0)) { - sharmonics_free(ws); - ws=0; - } - else { - ws->l=l; - ws->nmax=nmax; - } - - } - return ws; -} - -void sharmonics_free(SHWorkSpace *ws) { - if (ws!=0) { - legendre_a_free(ws->wl); - chebyshev12_free(ws->wc); - sharmonics_normalization_free(ws->wn); - free(ws); - } -} - -void sharmonics(double *phi, double *theta, SHWorkSpace *ws) { - - if ((ws!=0)&&(theta!=0)&&(phi!=0)) { - chebyshev12(phi,ws->wc); - legendre_a(theta,ws->wl); - - sharmonics_initlegendrenormlzd(ws); - sharmonics_multcheby12(ws); - } -} - - |