/* * 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_spherical_harmonics/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) { unsigned int n,m,ny0,np0; unsigned 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; ll; 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 normalized 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) { unsigned int n,m,ny0; const int nt0=ws->nmax; unsigned 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; ll; 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; ll; 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); } }