1 ////////////////////////////////////////////////////////////////////////////////
3 // AliFemtoYlm - the class to calculate varous components of spherical //
6 // Authors: Adam Kisiel kisiel@mps.ohio-state.edu //
8 ////////////////////////////////////////////////////////////////////////////////
10 #include "AliFemtoYlm.h"
14 double *AliFemtoYlm::fgPrefactors = 0x0;
15 int *AliFemtoYlm::fgPrefshift = 0x0;
16 int *AliFemtoYlm::fgPlmshift = 0x0;
18 AliFemtoYlm::AliFemtoYlm() {
22 AliFemtoYlm::~AliFemtoYlm() {
29 AliFemtoYlm::AliFemtoYlm(const AliFemtoYlm& aYlm){
30 fgPrefshift = aYlm.fgPrefshift;
34 AliFemtoYlm& AliFemtoYlm::operator=(const AliFemtoYlm& aYlm){
43 std::complex<double> AliFemtoYlm::Ceiphi(double phi){
44 return std::complex<double>(cos(phi),sin(phi));
47 double AliFemtoYlm::Legendre(int ell, int em, double ctheta){
48 // Calculate a single Legendre value
49 // *** Warning - NOT optimal - calculated all Plms up to L !!!
51 AliFemtoYlm::LegendreUpToYlm(ell, ctheta, lbuf);
53 return lbuf[fgPlmshift[ell]-abs(em)];
56 std::complex<double> AliFemtoYlm::Ylm(int ell,int m,double theta,double phi){
57 // Calculate Ylm spherical input
59 std::complex<double> answer;
60 std::complex<double> ci(0.0,1.0);
62 answer = (fgPrefactors[fgPrefshift[ell]+m]*Legendre(ell,m,ctheta))*Ceiphi(m*phi);
67 std::complex<double> AliFemtoYlm::Ylm(int ell, int m, double x, double y, double z){
68 // Calculate Ylm cartesian input
69 std::complex<double> answer;
71 double r = sqrt(x*x+y*y+z*z);
72 if ( r < 1e-10 || fabs(z) < 1e-10 ) ctheta = 0.0;
75 answer = (fgPrefactors[fgPrefshift[ell]+m]*Legendre(ell,m,ctheta))*Ceiphi(m*phi);
80 double AliFemtoYlm::ReYlm(int ell, int m, double theta, double phi){
81 return real(AliFemtoYlm::Ylm(ell,m,theta,phi));
84 double AliFemtoYlm::ImYlm(int ell, int m, double theta, double phi){
85 return imag(AliFemtoYlm::Ylm(ell,m,theta,phi));
88 double AliFemtoYlm::ReYlm(int ell, int m, double x,double y,double z){
89 return real(AliFemtoYlm::Ylm(ell,m,x,y,z));
92 double AliFemtoYlm::ImYlm(int ell, int m, double x,double y,double z){
93 return imag(AliFemtoYlm::Ylm(ell,m,x,y,z));
96 void AliFemtoYlm::InitializeYlms()
98 // Calculate prefactors for fast Ylm calculation
100 double oneoversqrtpi = 1.0/TMath::Sqrt(TMath::Pi());
102 if(fgPrefactors!=NULL)
104 if(fgPrefshift!=NULL)
109 fgPrefactors = (double *) malloc(sizeof(double) * 36);
110 fgPrefshift = (int *) malloc(sizeof(int) * 6);
111 fgPlmshift = (int *) malloc(sizeof(int) * 6);
114 fgPrefactors[0] = 0.5*oneoversqrtpi;
117 fgPrefactors[1] = 0.5*sqrt(3.0/2.0)*oneoversqrtpi;
118 fgPrefactors[2] = 0.5*sqrt(3.0)*oneoversqrtpi;
119 fgPrefactors[3] = -fgPrefactors[1];
122 fgPrefactors[4] = 0.25*sqrt(15.0/2.0)*oneoversqrtpi;
123 fgPrefactors[5] = 0.5*sqrt(15.0/2.0)*oneoversqrtpi;
124 fgPrefactors[6] = 0.25*sqrt(5.0)*oneoversqrtpi;
125 fgPrefactors[7] = -fgPrefactors[5];
126 fgPrefactors[8] = fgPrefactors[4];
129 fgPrefactors[9] = 0.125*sqrt(35.0)*oneoversqrtpi;
130 fgPrefactors[10] = 0.25*sqrt(105.0/2.0)*oneoversqrtpi;
131 fgPrefactors[11] = 0.125*sqrt(21.0)*oneoversqrtpi;
132 fgPrefactors[12] = 0.25*sqrt(7.0)*oneoversqrtpi;
133 fgPrefactors[13] = -fgPrefactors[11];
134 fgPrefactors[14] = fgPrefactors[10];
135 fgPrefactors[15] = -fgPrefactors[9];
138 fgPrefactors[16] = 3.0/16.0*sqrt(35.0/2.0)*oneoversqrtpi;
139 fgPrefactors[17] = 3.0/8.0*sqrt(35.0)*oneoversqrtpi;
140 fgPrefactors[18] = 3.0/8.0*sqrt(5.0/2.0)*oneoversqrtpi;
141 fgPrefactors[19] = 3.0/8.0*sqrt(5.0)*oneoversqrtpi;
142 fgPrefactors[20] = 3.0/16.0*oneoversqrtpi;
143 fgPrefactors[21] = -fgPrefactors[19];
144 fgPrefactors[22] = fgPrefactors[18];
145 fgPrefactors[23] = -fgPrefactors[17];
146 fgPrefactors[24] = fgPrefactors[16];
149 fgPrefactors[25] = 3.0/32.0*sqrt(77.0)*oneoversqrtpi;
150 fgPrefactors[26] = 3.0/16.0*sqrt(385.0/2.0)*oneoversqrtpi;
151 fgPrefactors[27] = 1.0/32.0*sqrt(385.0)*oneoversqrtpi;
152 fgPrefactors[28] = 1.0/8.0*sqrt(1155.0/2.0)*oneoversqrtpi;
153 fgPrefactors[29] = 1.0/16.0*sqrt(165.0/2.0)*oneoversqrtpi;
154 fgPrefactors[30] = 1.0/16.0*sqrt(11.0)*oneoversqrtpi;
155 fgPrefactors[31] = -fgPrefactors[29];
156 fgPrefactors[32] = fgPrefactors[28];
157 fgPrefactors[33] = -fgPrefactors[27];
158 fgPrefactors[34] = fgPrefactors[26];
159 fgPrefactors[35] = -fgPrefactors[25];
176 void AliFemtoYlm::LegendreUpToYlm(int lmax, double ctheta, double *lbuf)
178 // Calculate a set of legendre polynomials up to a given l
179 // with spherical input
184 sins[1] = sqrt(1-ctheta*ctheta);
186 for (int iter=2; iter<6; iter++) {
187 sins[iter] = sins[iter-1]*sins[1];
188 coss[iter] = coss[iter-1]*coss[1];
191 // Legendre polynomials l=0
194 // Legendre polynomials l=1
200 // Legendre polynomials l=2
203 lbuf[4] = sins[1]*coss[1];
204 lbuf[5] = 3*coss[2]-1;
207 // Legendre polynomials l=3
210 lbuf[7] = sins[2]*coss[1];
211 lbuf[8] = (5*coss[2]-1)*sins[1];
212 lbuf[9] = 5*coss[3]-3*coss[1];
215 // Legendre polynomials l=4
218 lbuf[11] = sins[3]*coss[1];
219 lbuf[12] = (7*coss[2]-1)*sins[2];
220 lbuf[13] = (7*coss[3]-3*coss[1])*sins[1];
221 lbuf[14] = 35*coss[4]-30*coss[2]+3;
224 // Legendre polynomials l=5
227 lbuf[16] = sins[4]*coss[1];
228 lbuf[17] = (9*coss[2]-1)*sins[3];
229 lbuf[18] = (3*coss[3]-1*coss[1])*sins[2];
230 lbuf[19] = (21*coss[4]-14*coss[2]+1)*sins[1];
231 lbuf[20] = 63*coss[5]-70*coss[3]+15*coss[1];
235 void AliFemtoYlm::YlmUpToL(int lmax, double x, double y, double z, std::complex<double> *ylms)
237 // Calculate a set of Ylms up to a given l
238 // with cartesian input
241 double r = sqrt(x*x+y*y+z*z);
242 if ( r < 1e-10 || fabs(z) < 1e-10 ) ctheta = 0.0;
246 YlmUpToL(lmax, ctheta, phi, ylms);
250 void AliFemtoYlm::YlmUpToL(int lmax, double ctheta, double phi, std::complex<double> *ylms)
252 // Calculate a set of Ylms up to a given l
253 // with spherical input
261 LegendreUpToYlm(lmax, ctheta, lbuf);
263 for (int iter=1; iter<=lmax; iter++) {
264 coss[iter-1] = cos(iter*phi);
265 sins[iter-1] = sin(iter*phi);
267 ylms[lcur++] = fgPrefactors[0]*lbuf[0] * std::complex<double>(1,0);
269 for (int il = 1; il<=lmax; il++) {
271 ylms[lcur+il] = fgPrefactors[fgPrefshift[il]]*lbuf[fgPlmshift[il]]*std::complex<double>(1.0,0.0);
273 for (int im=1; im<=il; im++) {
274 lpol = lbuf[fgPlmshift[il]-im];
275 ylms[lcur+il-im] = fgPrefactors[fgPrefshift[il]-im]*lpol*std::complex<double>(coss[im-1],-sins[im-1]);
276 ylms[lcur+il+im] = fgPrefactors[fgPrefshift[il]+im]*lpol*std::complex<double>(coss[im-1],sins[im-1]);