1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
|
/*
Evaluates all circular harmonics
at the angles phi up to the order nmax.
using the recurrence for the Chebyshev
polynomials of the first and second kind
T has the dimensions length(phi) x 2nmax+1
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/chebyshev12.h"
Cheby12WorkSpace *chebyshev12_alloc(const size_t nmax, const size_t l) {
Cheby12WorkSpace *wc;
// memory allocation
if ((wc=(Cheby12WorkSpace*)calloc(1,sizeof(Cheby12WorkSpace)))!=0) {
wc->l=l;
wc->nmax=nmax;
if ((wc->t=(double*)calloc(l*(2*nmax+1),sizeof(double)))==0) {
free(wc);
return 0;
}
return wc;
}
return 0;
}
void chebyshev12_free(Cheby12WorkSpace *wc) {
if (wc!=0) {
free(wc->t);
free(wc);
}
}
void chebyshev12(double *phi, Cheby12WorkSpace *wc) {
int l,l0,n;
const int incr=2*wc->nmax+1;
double *cosphi;
double *sinphi;
// memory allocation
if ((wc!=0)&&(phi!=0)) {
if ((cosphi=(double*)calloc(wc->l,sizeof(double)))==0) {
return;
}
if ((sinphi=(double*)calloc(wc->l,sizeof(double)))==0) {
free(cosphi);
return;
}
// constants and initialization
for (l=0, l0=wc->nmax; l<wc->l; l++, l0+=incr) {
cosphi[l]=cos(phi[l]);
sinphi[l]=sin(phi[l]);
// initial value T_0=1
wc->t[l0]=1;
wc->t[l0+1]=cosphi[l];
wc->t[l0-1]=sinphi[l];
}
// recurrence for n>1
for (n=2; n<=wc->nmax; n++) {
for (l=0, l0=wc->nmax; l<wc->l; l++, l0+=incr) {
wc->t[l0+n]=cosphi[l]* wc->t[l0+n-1] - sinphi[l]* wc->t[l0-n+1];
wc->t[l0-n]=sinphi[l]* wc->t[l0+n-1] + cosphi[l]* wc->t[l0-n+1];
}
}
free(cosphi);
free(sinphi);
}
}
|