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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
/*
Evaluates all associated legendre functions
at the angles theta up to the order nmax
using the three-term recurrence of the Legendre functions.
P has dimensions length(theta) x (nmax+1)(nmax+2)
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/legendre_a.h"
static void legendre_first_recurrence (double *sintheta, LegendreWorkSpace *wl) {
int n;
unsigned int l,l0;
int nmo0=0;
int n0=1;
const int incr=(wl->nmax+1)*(wl->nmax+2)/2;
// computes the legendre functions P_n^m(costheta) for m=n
// from P_0^0
for (n=1; n<=wl->nmax; n++) {
for (l=0,l0=0; l<wl->l; l++,l0+=incr) {
wl->p[l0+n0+n] = -(2*n-1) * wl->p[l0+nmo0+n-1] * sintheta[l];
}
nmo0=n0;
n0+=n+1;
}
}
static void legendre_second_recurrence (double *costheta, LegendreWorkSpace *wl) {
unsigned int m,n,l,l0;
int nmt0=-1;
int nmo0=0;
int n0=1;
const int incr=(wl->nmax+1)*(wl->nmax+2)/2;
// computes the Legendre functions P_n^m(costheta) from
// P_n^m with m=n
for (n=1; n<=wl->nmax; n++) {
for (m=0; m<n; m++) {
if (m<=n-2) {
for (l=0,l0=0; l<wl->l; l++,l0+=incr) {
wl->p[l0+n0+m] = (
(2*n-1) * costheta[l] * wl->p[l0+nmo0+m]
- (n+m-1) * wl->p[l0+nmt0+m]
) / (n-m);
}
}
else {
for (l=0,l0=0; l<wl->l; l++,l0+=incr) {
wl->p[l0+n0+m] = (
(2*n-1) * costheta[l] * wl->p[l0+nmo0+m]
) / (n-m);
}
}
}
nmt0=nmo0;
nmo0=n0;
n0+=n+1;
}
}
LegendreWorkSpace *legendre_a_alloc(const size_t nmax, const size_t l) {
LegendreWorkSpace *wl;
// memory allocation
if ((wl=(LegendreWorkSpace*)calloc(1,sizeof(LegendreWorkSpace)))!=0) {
wl->l=l;
wl->nmax=nmax;
if ((wl->p=(double*)calloc(l*(nmax+1)*(nmax+2)/2,sizeof(double)))==0) {
free(wl);
wl = 0;
}
}
return wl;
}
void legendre_a_free(LegendreWorkSpace *wl) {
if (wl!=0) {
free(wl->p);
free(wl);
}
}
void legendre_a(double *theta, LegendreWorkSpace *wl) {
unsigned int l,l0;
const int incr=(wl->nmax+1)*(wl->nmax+2)/2;
double *costheta;
double *sintheta;
// memory allocation
if ((wl!=0)&&(theta!=0)) {
if ((costheta=(double*)calloc(wl->l,sizeof(double)))==0) {
return;
}
if ((sintheta=(double*)calloc(wl->l,sizeof(double)))==0) {
free(costheta);
return;
}
// constants and initialization
for (l=0, l0=0; l<wl->l; l++, l0+=incr) {
costheta[l]=cos(theta[l]);
sintheta[l]=sin(theta[l]);
// initial value P_0^0=1
wl->p[l0]=1;
}
// recurrences evaluating all the Legendre functions
legendre_first_recurrence(sintheta,wl);
legendre_second_recurrence(costheta,wl);
free(sintheta);
free(costheta);
}
}
|
