]>
Commit | Line | Data |
---|---|---|
17b6dded | 1 | //////////////////////////////////////////////////////////////////////////////// |
2 | // // | |
3 | // AliFemtoYlm - the class to calculate varous components of spherical // | |
4 | // harmonics // | |
5 | // // | |
6 | // Authors: Adam Kisiel kisiel@mps.ohio-state.edu // | |
7 | // // | |
8 | //////////////////////////////////////////////////////////////////////////////// | |
9 | ||
10 | #include "AliFemtoYlm.h" | |
11 | #include <TMath.h> | |
12 | #include <iostream> | |
13 | ||
14 | double *AliFemtoYlm::fgPrefactors = 0x0; | |
15 | int *AliFemtoYlm::fgPrefshift = 0x0; | |
16 | int *AliFemtoYlm::fgPlmshift = 0x0; | |
17 | ||
18 | AliFemtoYlm::AliFemtoYlm() { | |
19 | InitializeYlms(); | |
20 | } | |
21 | ||
22 | AliFemtoYlm::~AliFemtoYlm() {} | |
23 | ||
24 | ||
25 | AliFemtoYlm::AliFemtoYlm(const AliFemtoYlm& aYlm){ | |
04328c66 | 26 | fgPrefshift = aYlm.fgPrefshift; |
17b6dded | 27 | InitializeYlms(); |
28 | } | |
29 | ||
30 | AliFemtoYlm& AliFemtoYlm::operator=(const AliFemtoYlm& aYlm){ | |
31 | if (this == &aYlm) | |
32 | return *this; | |
33 | ||
34 | InitializeYlms(); | |
35 | ||
36 | return *this; | |
37 | } | |
38 | ||
39 | std::complex<double> AliFemtoYlm::Ceiphi(double phi){ | |
40 | return std::complex<double>(cos(phi),sin(phi)); | |
41 | } | |
42 | ||
43 | double AliFemtoYlm::Legendre(int ell, int em, double ctheta){ | |
44 | // Calculate a single Legendre value | |
45 | // *** Warning - NOT optimal - calculated all Plms up to L !!! | |
46 | double lbuf[36]; | |
47 | AliFemtoYlm::LegendreUpToYlm(ell, ctheta, lbuf); | |
48 | ||
49 | return lbuf[fgPlmshift[ell]-abs(em)]; | |
50 | } | |
51 | ||
52 | std::complex<double> AliFemtoYlm::Ylm(int ell,int m,double theta,double phi){ | |
53 | // Calculate Ylm spherical input | |
54 | double ctheta; | |
55 | std::complex<double> answer; | |
56 | std::complex<double> ci(0.0,1.0); | |
57 | ctheta=cos(theta); | |
58 | answer = (fgPrefactors[fgPrefshift[ell]+m]*Legendre(ell,m,ctheta))*Ceiphi(m*phi); | |
59 | ||
60 | return answer; | |
61 | } | |
62 | ||
63 | std::complex<double> AliFemtoYlm::Ylm(int ell, int m, double x, double y, double z){ | |
64 | // Calculate Ylm cartesian input | |
65 | std::complex<double> answer; | |
66 | double ctheta,phi; | |
67 | double r = sqrt(x*x+y*y+z*z); | |
68 | if ( r < 1e-10 || fabs(z) < 1e-10 ) ctheta = 0.0; | |
69 | else ctheta=z/r; | |
70 | phi=atan2(y,x); | |
71 | answer = (fgPrefactors[fgPrefshift[ell]+m]*Legendre(ell,m,ctheta))*Ceiphi(m*phi); | |
72 | ||
73 | return answer; | |
74 | } | |
75 | ||
76 | double AliFemtoYlm::ReYlm(int ell, int m, double theta, double phi){ | |
77 | return real(AliFemtoYlm::Ylm(ell,m,theta,phi)); | |
78 | } | |
79 | ||
80 | double AliFemtoYlm::ImYlm(int ell, int m, double theta, double phi){ | |
81 | return imag(AliFemtoYlm::Ylm(ell,m,theta,phi)); | |
82 | } | |
83 | ||
84 | double AliFemtoYlm::ReYlm(int ell, int m, double x,double y,double z){ | |
85 | return real(AliFemtoYlm::Ylm(ell,m,x,y,z)); | |
86 | } | |
87 | ||
88 | double AliFemtoYlm::ImYlm(int ell, int m, double x,double y,double z){ | |
89 | return imag(AliFemtoYlm::Ylm(ell,m,x,y,z)); | |
90 | } | |
91 | ||
92 | void AliFemtoYlm::InitializeYlms() | |
93 | { | |
94 | // Calculate prefactors for fast Ylm calculation | |
95 | ||
96 | double oneoversqrtpi = 1.0/TMath::Sqrt(TMath::Pi()); | |
97 | ||
98 | fgPrefactors = (double *) malloc(sizeof(double) * 36); | |
99 | fgPrefshift = (int *) malloc(sizeof(int) * 6); | |
100 | fgPlmshift = (int *) malloc(sizeof(int) * 6); | |
101 | ||
102 | // l=0 prefactors | |
103 | fgPrefactors[0] = 0.5*oneoversqrtpi; | |
104 | ||
105 | // l=1 prefactors | |
106 | fgPrefactors[1] = 0.5*sqrt(3.0/2.0)*oneoversqrtpi; | |
107 | fgPrefactors[2] = 0.5*sqrt(3.0)*oneoversqrtpi; | |
108 | fgPrefactors[3] = -fgPrefactors[1]; | |
109 | ||
110 | // l=2 prefactors | |
111 | fgPrefactors[4] = 0.25*sqrt(15.0/2.0)*oneoversqrtpi; | |
112 | fgPrefactors[5] = 0.5*sqrt(15.0/2.0)*oneoversqrtpi; | |
113 | fgPrefactors[6] = 0.25*sqrt(5.0)*oneoversqrtpi; | |
114 | fgPrefactors[7] = -fgPrefactors[5]; | |
115 | fgPrefactors[8] = fgPrefactors[4]; | |
116 | ||
117 | // l=3 prefactors | |
118 | fgPrefactors[9] = 0.125*sqrt(35.0)*oneoversqrtpi; | |
119 | fgPrefactors[10] = 0.25*sqrt(105.0/2.0)*oneoversqrtpi; | |
120 | fgPrefactors[11] = 0.125*sqrt(21.0)*oneoversqrtpi; | |
121 | fgPrefactors[12] = 0.25*sqrt(7.0)*oneoversqrtpi; | |
122 | fgPrefactors[13] = -fgPrefactors[11]; | |
123 | fgPrefactors[14] = fgPrefactors[10]; | |
124 | fgPrefactors[15] = -fgPrefactors[9]; | |
125 | ||
126 | // l=4 prefactors | |
127 | fgPrefactors[16] = 3.0/16.0*sqrt(35.0/2.0)*oneoversqrtpi; | |
128 | fgPrefactors[17] = 3.0/8.0*sqrt(35.0)*oneoversqrtpi; | |
129 | fgPrefactors[18] = 3.0/8.0*sqrt(5.0/2.0)*oneoversqrtpi; | |
130 | fgPrefactors[19] = 3.0/8.0*sqrt(5.0)*oneoversqrtpi; | |
131 | fgPrefactors[20] = 3.0/16.0*oneoversqrtpi; | |
132 | fgPrefactors[21] = -fgPrefactors[19]; | |
133 | fgPrefactors[22] = fgPrefactors[18]; | |
134 | fgPrefactors[23] = -fgPrefactors[17]; | |
135 | fgPrefactors[24] = fgPrefactors[16]; | |
136 | ||
137 | // l=5 prefactors | |
138 | fgPrefactors[25] = 3.0/32.0*sqrt(77.0)*oneoversqrtpi; | |
139 | fgPrefactors[26] = 3.0/16.0*sqrt(385.0/2.0)*oneoversqrtpi; | |
140 | fgPrefactors[27] = 1.0/32.0*sqrt(385.0)*oneoversqrtpi; | |
141 | fgPrefactors[28] = 1.0/8.0*sqrt(1155.0/2.0)*oneoversqrtpi; | |
142 | fgPrefactors[29] = 1.0/16.0*sqrt(165.0/2.0)*oneoversqrtpi; | |
143 | fgPrefactors[30] = 1.0/16.0*sqrt(11.0)*oneoversqrtpi; | |
144 | fgPrefactors[31] = -fgPrefactors[29]; | |
145 | fgPrefactors[32] = fgPrefactors[28]; | |
146 | fgPrefactors[33] = -fgPrefactors[27]; | |
147 | fgPrefactors[34] = fgPrefactors[26]; | |
148 | fgPrefactors[35] = -fgPrefactors[25]; | |
149 | ||
150 | fgPrefshift[0] = 0; | |
151 | fgPrefshift[1] = 2; | |
152 | fgPrefshift[2] = 6; | |
153 | fgPrefshift[3] = 12; | |
154 | fgPrefshift[4] = 20; | |
155 | fgPrefshift[5] = 30; | |
156 | ||
157 | fgPlmshift[0] = 0; | |
158 | fgPlmshift[1] = 2; | |
159 | fgPlmshift[2] = 5; | |
160 | fgPlmshift[3] = 9; | |
161 | fgPlmshift[4] = 14; | |
162 | fgPlmshift[5] = 20; | |
163 | } | |
164 | ||
165 | void AliFemtoYlm::LegendreUpToYlm(int lmax, double ctheta, double *lbuf) | |
166 | { | |
167 | // Calculate a set of legendre polynomials up to a given l | |
168 | // with spherical input | |
169 | double sins[6]; | |
170 | double coss[6]; | |
171 | sins[0] = 0.0; | |
172 | coss[0] = 1.0; | |
173 | sins[1] = sqrt(1-ctheta*ctheta); | |
174 | coss[1] = ctheta; | |
175 | for (int iter=2; iter<6; iter++) { | |
176 | sins[iter] = sins[iter-1]*sins[1]; | |
177 | coss[iter] = coss[iter-1]*coss[1]; | |
178 | } | |
179 | ||
180 | // Legendre polynomials l=0 | |
181 | lbuf[0] = 1.0; | |
182 | ||
183 | // Legendre polynomials l=1 | |
184 | if (lmax>0) { | |
185 | lbuf[1] = sins[1]; | |
186 | lbuf[2] = coss[1]; | |
187 | } | |
188 | ||
189 | // Legendre polynomials l=2 | |
190 | if (lmax>1) { | |
191 | lbuf[3] = sins[2]; | |
192 | lbuf[4] = sins[1]*coss[1]; | |
193 | lbuf[5] = 3*coss[2]-1; | |
194 | } | |
195 | ||
196 | // Legendre polynomials l=3 | |
197 | if (lmax>2) { | |
198 | lbuf[6] = sins[3]; | |
199 | lbuf[7] = sins[2]*coss[1]; | |
200 | lbuf[8] = (5*coss[2]-1)*sins[1]; | |
201 | lbuf[9] = 5*coss[3]-3*coss[1]; | |
202 | } | |
203 | ||
204 | // Legendre polynomials l=4 | |
205 | if (lmax>3) { | |
206 | lbuf[10] = sins[4]; | |
207 | lbuf[11] = sins[3]*coss[1]; | |
208 | lbuf[12] = (7*coss[2]-1)*sins[2]; | |
209 | lbuf[13] = (7*coss[3]-3*coss[1])*sins[1]; | |
210 | lbuf[14] = 35*coss[4]-30*coss[2]+3; | |
211 | } | |
212 | ||
213 | // Legendre polynomials l=5 | |
214 | if (lmax>4) { | |
215 | lbuf[15] = sins[5]; | |
216 | lbuf[16] = sins[4]*coss[1]; | |
217 | lbuf[17] = (9*coss[2]-1)*sins[3]; | |
218 | lbuf[18] = (3*coss[3]-1*coss[1])*sins[2]; | |
219 | lbuf[19] = (21*coss[4]-14*coss[2]+1)*sins[1]; | |
220 | lbuf[20] = 63*coss[5]-70*coss[3]+15*coss[1]; | |
221 | } | |
222 | } | |
223 | ||
224 | void AliFemtoYlm::YlmUpToL(int lmax, double x, double y, double z, std::complex<double> *ylms) | |
225 | { | |
226 | // Calculate a set of Ylms up to a given l | |
227 | // with cartesian input | |
228 | double ctheta,phi; | |
229 | ||
230 | double r = sqrt(x*x+y*y+z*z); | |
231 | if ( r < 1e-10 || fabs(z) < 1e-10 ) ctheta = 0.0; | |
232 | else ctheta=z/r; | |
233 | phi=atan2(y,x); | |
234 | ||
235 | YlmUpToL(lmax, ctheta, phi, ylms); | |
236 | ||
237 | } | |
238 | ||
239 | void AliFemtoYlm::YlmUpToL(int lmax, double ctheta, double phi, std::complex<double> *ylms) | |
240 | { | |
241 | // Calculate a set of Ylms up to a given l | |
242 | // with spherical input | |
243 | int lcur = 0; | |
244 | double lpol; | |
245 | ||
246 | double coss[6]; | |
247 | double sins[6]; | |
248 | ||
249 | double lbuf[36]; | |
250 | LegendreUpToYlm(lmax, ctheta, lbuf); | |
251 | ||
252 | for (int iter=1; iter<=lmax; iter++) { | |
253 | coss[iter-1] = cos(iter*phi); | |
254 | sins[iter-1] = sin(iter*phi); | |
255 | } | |
256 | ylms[lcur++] = fgPrefactors[0]*lbuf[0] * std::complex<double>(1,0); | |
257 | ||
258 | for (int il = 1; il<=lmax; il++) { | |
259 | // First im = 0 | |
260 | ylms[lcur+il] = fgPrefactors[fgPrefshift[il]]*lbuf[fgPlmshift[il]]*std::complex<double>(1.0,0.0); | |
261 | // Im != 0 | |
262 | for (int im=1; im<=il; im++) { | |
263 | lpol = lbuf[fgPlmshift[il]-im]; | |
264 | ylms[lcur+il-im] = fgPrefactors[fgPrefshift[il]-im]*lpol*std::complex<double>(coss[im-1],-sins[im-1]); | |
265 | ylms[lcur+il+im] = fgPrefactors[fgPrefshift[il]+im]*lpol*std::complex<double>(coss[im-1],sins[im-1]); | |
266 | } | |
267 | lcur += 2*il + 1; | |
268 | } | |
269 | } |