Heavy ion Analysis Libriares
Loading...
Searching...
No Matches
FemtoYlmMath.cxx
1/*
2 * FemtoYlmMathMath.cxx
3 *
4 * Created on: 25 sie 2022
5 * Author: Daniel Wielanek
6 * E-mail: daniel.wielanek@gmail.com
7 * Warsaw University of Technology, Faculty of Physics
8 */
9#include "FemtoYlmMath.h"
10
11
12#include <TMath.h>
13#include <cmath>
14#include <cstdlib>
15#include <iostream>
16
17namespace Hal {
18
19 FemtoYlmMath::FemtoYlmMath() { InitializeYlms(); }
20
21 FemtoYlmMath::~FemtoYlmMath() {}
22
23 FemtoYlmMath::FemtoYlmMath(const FemtoYlmMath& aYlm) : TObject(aYlm) { InitializeYlms(); }
24
25 FemtoYlmMath& FemtoYlmMath::operator=(const FemtoYlmMath& aYlm) {
26 if (this == &aYlm) return *this;
27
28 InitializeYlms();
29
30 return *this;
31 }
32
33 std::complex<double> FemtoYlmMath::Ceiphi(double phi) const { return std::complex<double>(cos(phi), sin(phi)); }
34
35 double FemtoYlmMath::Legendre(int ell, int em, double ctheta) const {
36 // Calculate a single Legendre value
37 // *** Warning - NOT optimal - calculated all Plms up to L !!!
38 double lbuf[36];
39 FemtoYlmMath::LegendreUpToYlm(ell, ctheta, lbuf);
40
41 return lbuf[fPlmshift[ell] - abs(em)];
42 }
43 double FemtoYlmMath::DeltaJ(double aJot1, double aJot2, double aJot) const {
44 if ((aJot1 + aJot2 - aJot) < 0) { return 0; }
45 if ((aJot1 - aJot2 + aJot) < 0) { return 0; }
46 if ((-aJot1 + aJot2 + aJot) < 0) { return 0; }
47 if ((aJot1 + aJot2 + aJot + 1) < 0) { return 0; }
48 double res = TMath::Sqrt(1.0 * fFactorials[lrint(aJot1 + aJot2 - aJot)] * fFactorials[lrint(aJot1 - aJot2 + aJot)]
49 * fFactorials[lrint(-aJot1 + aJot2 + aJot)] / fFactorials[lrint(aJot1 + aJot2 + aJot + 1)]);
50 return res;
51 }
52
53 double FemtoYlmMath::ClebschGordan(double aJot1, double aEm1, double aJot2, double aEm2, double aJot, double aEm) const {
54 int mint, maxt;
55 double cgc = 0.0;
56 int titer;
57 double coef;
58
59 maxt = lrint(aJot1 + aJot2 - aJot);
60 mint = 0;
61 if (lrint(aJot1 - aEm1) < maxt) maxt = lrint(aJot1 - aEm1);
62 if (lrint(aJot2 + aEm2) < maxt) maxt = lrint(aJot2 + aEm2);
63 if (lrint(-(aJot - aJot2 + aEm1)) > mint) mint = lrint(-(aJot - aJot2 + aEm1));
64 if (lrint(-(aJot - aJot1 - aEm2)) > mint) mint = lrint(-(aJot - aJot1 - aEm2));
65
66 for (titer = mint; titer <= maxt; titer++) {
67 coef = TMath::Power(-1, titer);
68 coef *= TMath::Sqrt((2 * aJot + 1) * fFactorials[lrint(aJot1 + aEm1)] * fFactorials[lrint(aJot1 - aEm1)]
69 * fFactorials[lrint(aJot2 + aEm2)] * fFactorials[lrint(aJot2 - aEm2)] * fFactorials[lrint(aJot + aEm)]
70 * fFactorials[lrint(aJot - aEm)]);
71 coef /= (fFactorials[titer] * fFactorials[lrint(aJot1 + aJot2 - aJot - titer)] * fFactorials[lrint(aJot1 - aEm1 - titer)]
72 * fFactorials[lrint(aJot2 + aEm2 - titer)] * fFactorials[lrint(aJot - aJot2 + aEm1 + titer)]
73 * fFactorials[lrint(aJot - aJot1 - aEm2 + titer)]);
74
75 cgc += coef;
76 }
77
78 cgc *= DeltaJ(aJot1, aJot2, aJot);
79
80 return cgc;
81 }
82
83 double FemtoYlmMath::WignerSymbol(double aJot1, double aEm1, double aJot2, double aEm2, double aJot, double aEm) const {
84 if (lrint(aEm1 + aEm2 + aEm) != 0.0) return 0.0;
85 double cge = ClebschGordan(aJot1, aEm1, aJot2, aEm2, aJot, -aEm);
86 if (lrint(abs(aJot1 - aJot2 - aEm)) % 2) cge *= -1.0;
87 cge /= sqrt(2 * aJot + 1);
88
89 if (cge == -0.0) cge = 0.0;
90
91 return cge;
92 }
93
94 std::complex<double> FemtoYlmMath::Ylm(int ell, int m, double theta, double phi) const {
95 // Calculate Ylm spherical input
96 double ctheta;
97 std::complex<double> answer;
98 std::complex<double> ci(0.0, 1.0);
99 ctheta = cos(theta);
100 answer = (fPrefactors[fPrefshift[ell] + m] * Legendre(ell, m, ctheta)) * Ceiphi(m * phi);
101
102 return answer;
103 }
104
105 std::complex<double> FemtoYlmMath::Ylm(int ell, int m, double x, double y, double z) const {
106 // Calculate Ylm cartesian input
107 std::complex<double> answer;
108 double ctheta, phi;
109 double r = sqrt(x * x + y * y + z * z);
110 if (r < 1e-10 || fabs(z) < 1e-10)
111 ctheta = 0.0;
112 else
113 ctheta = z / r;
114 phi = atan2(y, x);
115 answer = (fPrefactors[fPrefshift[ell] + m] * Legendre(ell, m, ctheta)) * Ceiphi(m * phi);
116
117 return answer;
118 }
119
120 double FemtoYlmMath::ReYlm(int ell, int m, double theta, double phi) const {
121 return real(FemtoYlmMath::Ylm(ell, m, theta, phi));
122 }
123
124 double FemtoYlmMath::ImYlm(int ell, int m, double theta, double phi) const {
125 return imag(FemtoYlmMath::Ylm(ell, m, theta, phi));
126 }
127
128 double FemtoYlmMath::ReYlm(int ell, int m, double x, double y, double z) const {
129 return real(FemtoYlmMath::Ylm(ell, m, x, y, z));
130 }
131
132 double FemtoYlmMath::ImYlm(int ell, int m, double x, double y, double z) const {
133 return imag(FemtoYlmMath::Ylm(ell, m, x, y, z));
134 }
135
136 void FemtoYlmMath::InitializeYlms() {
137 // Calculate prefactors for fast Ylm calculation
138
139 double oneoversqrtpi = 1.0 / TMath::Sqrt(TMath::Pi());
140
141
142 // l=0 prefactors
143 fPrefactors[0] = 0.5 * oneoversqrtpi;
144
145 // l=1 prefactors
146 fPrefactors[1] = 0.5 * sqrt(3.0 / 2.0) * oneoversqrtpi;
147 fPrefactors[2] = 0.5 * sqrt(3.0) * oneoversqrtpi;
148 fPrefactors[3] = -fPrefactors[1];
149
150 // l=2 prefactors
151 fPrefactors[4] = 0.25 * sqrt(15.0 / 2.0) * oneoversqrtpi;
152 fPrefactors[5] = 0.5 * sqrt(15.0 / 2.0) * oneoversqrtpi;
153 fPrefactors[6] = 0.25 * sqrt(5.0) * oneoversqrtpi;
154 fPrefactors[7] = -fPrefactors[5];
155 fPrefactors[8] = fPrefactors[4];
156
157 // l=3 prefactors
158 fPrefactors[9] = 0.125 * sqrt(35.0) * oneoversqrtpi;
159 fPrefactors[10] = 0.25 * sqrt(105.0 / 2.0) * oneoversqrtpi;
160 fPrefactors[11] = 0.125 * sqrt(21.0) * oneoversqrtpi;
161 fPrefactors[12] = 0.25 * sqrt(7.0) * oneoversqrtpi;
162 fPrefactors[13] = -fPrefactors[11];
163 fPrefactors[14] = fPrefactors[10];
164 fPrefactors[15] = -fPrefactors[9];
165
166 // l=4 prefactors
167 fPrefactors[16] = 3.0 / 16.0 * sqrt(35.0 / 2.0) * oneoversqrtpi;
168 fPrefactors[17] = 3.0 / 8.0 * sqrt(35.0) * oneoversqrtpi;
169 fPrefactors[18] = 3.0 / 8.0 * sqrt(5.0 / 2.0) * oneoversqrtpi;
170 fPrefactors[19] = 3.0 / 8.0 * sqrt(5.0) * oneoversqrtpi;
171 fPrefactors[20] = 3.0 / 16.0 * oneoversqrtpi;
172 fPrefactors[21] = -fPrefactors[19];
173 fPrefactors[22] = fPrefactors[18];
174 fPrefactors[23] = -fPrefactors[17];
175 fPrefactors[24] = fPrefactors[16];
176
177 // l=5 prefactors
178 fPrefactors[25] = 3.0 / 32.0 * sqrt(77.0) * oneoversqrtpi;
179 fPrefactors[26] = 3.0 / 16.0 * sqrt(385.0 / 2.0) * oneoversqrtpi;
180 fPrefactors[27] = 1.0 / 32.0 * sqrt(385.0) * oneoversqrtpi;
181 fPrefactors[28] = 1.0 / 8.0 * sqrt(1155.0 / 2.0) * oneoversqrtpi;
182 fPrefactors[29] = 1.0 / 16.0 * sqrt(165.0 / 2.0) * oneoversqrtpi;
183 fPrefactors[30] = 1.0 / 16.0 * sqrt(11.0) * oneoversqrtpi;
184 fPrefactors[31] = -fPrefactors[29];
185 fPrefactors[32] = fPrefactors[28];
186 fPrefactors[33] = -fPrefactors[27];
187 fPrefactors[34] = fPrefactors[26];
188 fPrefactors[35] = -fPrefactors[25];
189
190 fPrefshift[0] = 0;
191 fPrefshift[1] = 2;
192 fPrefshift[2] = 6;
193 fPrefshift[3] = 12;
194 fPrefshift[4] = 20;
195 fPrefshift[5] = 30;
196
197 fPlmshift[0] = 0;
198 fPlmshift[1] = 2;
199 fPlmshift[2] = 5;
200 fPlmshift[3] = 9;
201 fPlmshift[4] = 14;
202 fPlmshift[5] = 20;
203 fYlms = new std::complex<double>[36];
204
205 Double_t fac = 1;
206 fFactorials[0] = 1;
207 for (int iter = 1; iter < 24; iter++) {
208 fac *= iter;
209 fFactorials[iter] = fac;
210 }
211 }
212
213 void FemtoYlmMath::LegendreUpToYlm(int lmax, double ctheta, double* lbuf) const {
214 // Calculate a set of legendre polynomials up to a given l
215 // with spherical input
216 double sins[6];
217 double coss[6];
218 sins[0] = 0.0;
219 coss[0] = 1.0;
220 sins[1] = sqrt(1 - ctheta * ctheta);
221 coss[1] = ctheta;
222 for (int iter = 2; iter < 6; iter++) {
223 sins[iter] = sins[iter - 1] * sins[1];
224 coss[iter] = coss[iter - 1] * coss[1];
225 }
226
227 // Legendre polynomials l=0
228 lbuf[0] = 1.0;
229
230 // Legendre polynomials l=1
231 if (lmax > 0) {
232 lbuf[1] = sins[1];
233 lbuf[2] = coss[1];
234 }
235
236 // Legendre polynomials l=2
237 if (lmax > 1) {
238 lbuf[3] = sins[2];
239 lbuf[4] = sins[1] * coss[1];
240 lbuf[5] = 3 * coss[2] - 1;
241 }
242
243 // Legendre polynomials l=3
244 if (lmax > 2) {
245 lbuf[6] = sins[3];
246 lbuf[7] = sins[2] * coss[1];
247 lbuf[8] = (5 * coss[2] - 1) * sins[1];
248 lbuf[9] = 5 * coss[3] - 3 * coss[1];
249 }
250
251 // Legendre polynomials l=4
252 if (lmax > 3) {
253 lbuf[10] = sins[4];
254 lbuf[11] = sins[3] * coss[1];
255 lbuf[12] = (7 * coss[2] - 1) * sins[2];
256 lbuf[13] = (7 * coss[3] - 3 * coss[1]) * sins[1];
257 lbuf[14] = 35 * coss[4] - 30 * coss[2] + 3;
258 }
259
260 // Legendre polynomials l=5
261 if (lmax > 4) {
262 lbuf[15] = sins[5];
263 lbuf[16] = sins[4] * coss[1];
264 lbuf[17] = (9 * coss[2] - 1) * sins[3];
265 lbuf[18] = (3 * coss[3] - 1 * coss[1]) * sins[2];
266 lbuf[19] = (21 * coss[4] - 14 * coss[2] + 1) * sins[1];
267 lbuf[20] = 63 * coss[5] - 70 * coss[3] + 15 * coss[1];
268 }
269 }
270
271 std::complex<double>* FemtoYlmMath::YlmUpToL(int lmax, double x, double y, double z) {
272 // Calculate a set of Ylms up to a given l
273 // with cartesian input
274 double ctheta, phi;
275
276 double r = sqrt(x * x + y * y + z * z);
277 if (r < 1e-10 || fabs(z) < 1e-10)
278 ctheta = 0.0;
279 else
280 ctheta = z / r;
281 phi = atan2(y, x);
282
283 return YlmUpToL(lmax, ctheta, phi);
284 }
285
286 std::complex<double>* FemtoYlmMath::YlmUpToL(int lmax, double ctheta, double phi) {
287 // Calculate a set of Ylms up to a given l
288 // with spherical input
289 int lcur = 0;
290 double lpol;
291
292 double coss[6];
293 double sins[6];
294
295 double lbuf[21]; // was 36 but not used
296 LegendreUpToYlm(lmax, ctheta, lbuf);
297
298 for (int iter = 1; iter <= lmax; iter++) {
299 coss[iter - 1] = cos(iter * phi);
300 sins[iter - 1] = sin(iter * phi);
301 }
302 fYlms[lcur++] = fPrefactors[0] * lbuf[0] * std::complex<double>(1, 0);
303
304 for (int il = 1; il <= lmax; il++) {
305 // First im = 0
306 fYlms[lcur + il] = fPrefactors[fPrefshift[il]] * lbuf[fPlmshift[il]] * std::complex<double>(1.0, 0.0);
307 // Im != 0
308 for (int im = 1; im <= il; im++) {
309 lpol = lbuf[fPlmshift[il] - im];
310 fYlms[lcur + il - im] = fPrefactors[fPrefshift[il] - im] * lpol * std::complex<double>(coss[im - 1], -sins[im - 1]);
311 fYlms[lcur + il + im] = fPrefactors[fPrefshift[il] + im] * lpol * std::complex<double>(coss[im - 1], sins[im - 1]);
312 }
313 lcur += 2 * il + 1;
314 }
315 return fYlms;
316 }
317} // namespace Hal