LALSimulation  5.4.0.1-fe68b98
LALSimIMRPhenomTPHM_FrameRotations.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2020 Hector Estelles
3  *
4  * This program is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation; either version 2 of the License, or
7  * (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with with program; see the file COPYING. If not, write to the
16  * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17  * MA 02110-1301 USA
18  */
19 
20 
21 /**
22  * \author Hector Estelles
23  */
24 
25 #include <math.h>
26 
27 /* LAL Header Files */
28 #include <lal/LALSimIMR.h>
29 #include <lal/LALConstants.h>
30 #include <lal/Date.h>
31 #include <lal/Units.h>
32 
33 /* GSL Header Files */
34 #include <gsl/gsl_linalg.h>
35 #include <gsl/gsl_complex.h>
36 #include <gsl/gsl_complex_math.h>
37 #include <gsl/gsl_roots.h>
38 
39 #include <lal/XLALGSL.h>
40 #include <gsl/gsl_errno.h>
41 #include <gsl/gsl_spline.h>
42 
43 
44 #include <gsl/gsl_sf_legendre.h>
45 #include <gsl/gsl_sf_gamma.h>
46 
47 #include <lal/LALSimInspiralPrecess.h>
48 
49 
50 /* Routines to perform frame rotations for a set of Spin-Weighted Spherical Harmonic modes. */
51 
52 /****************************************/
53 /********* FUNCTION DEFINITION **********/
54 /****************************************/
55 
56 /* Struct for storing recurrent squared roots in the Wigner coefficients */
57 typedef struct tagPhenomT_precomputed_sqrt
58 {
59  REAL8 sqrt2, sqrt2half, sqrt3, sqrt5, sqrt6, sqrt7, sqrt10, sqrt14, sqrt15, sqrt21, sqrt30, sqrt35, sqrt70, sqrt210;
60 } PhenomT_precomputed_sqrt;
61 
62 /* Struct to store Wignerd coefficients and powers of the exponentials of the precessing angles */
63 typedef struct tagPhenomTPWignerStruct
64 {
65  REAL8 wignerdL2[5][5];
66  REAL8 wignerdL3[7][7];
67  REAL8 wignerdL4[9][9];
68  REAL8 wignerdL5[11][11];
69 
70  COMPLEX16 expAlphaPowers[11];
71  COMPLEX16 expGammaPowers[11];
72 } PhenomTPWignerStruct;
73 
74 void IMRPhenomTPHM_SetPrecomputedSqrt(PhenomT_precomputed_sqrt *SQRT);
75 
76 /* Function to set the PhenomTPWignerStruct given some Euler angles.
77 Sign of beta has to be passed independently since this information is lost in cosBeta.
78 globalRot flag specifies how many coefficients have to be computed, since for the time dependent
79 rotation between the co-precessing and the J-frame not all are needed. */
80 void IMRPhenomTPHM_SetWignerDStruct( PhenomTPWignerStruct *wS, /**< Wigner struct */
81  PhenomT_precomputed_sqrt *SQRT, /**< Precomputed squared roots */
82  REAL8 cosBeta, /**< cosinus of opening angle beta */
83  REAL8 alpha, /**< precessing angle alpha */
84  REAL8 gamma, /**< precessing angle gamma */
85  INT4 LMAX, /**< Maximum l mode number */
86  INT4 sign, /**< Sign of beta */
87  UINT4 globalRot /**< Compute coefficients for global or local rotation */
88  );
89 
90 /* Function to compute the WignerDMatrix. Analogous to XLALWignerDMatrix (which employs XLALWignerdMatrix), but Wignerd coefficients
91 here are read from the precomputed struct, in order to improve computational speed. */
92 COMPLEX16 PhenomTWignerDMatrix(
93  int l, /**< mode number l */
94  int mp, /**< mode number m' */
95  int m, /**< mode number m */
96  PhenomTPWignerStruct *wS /**< Wigner Struct, which contains Wignerd elements and powers of alpha and gamma */
97  );
98 
99 /* Function to rotate a given set of Spin-Weighted Spherical Harmonic mode time series with an time-dependent Euler rotation specified by three Euler angles.*/
100 int PhenomTPHM_RotateModes(
101  SphHarmTimeSeries* h_lm, /**< spherical harmonic decomposed modes, modified in place */
102  REAL8TimeSeries* alpha, /**< alpha Euler angle time series */
103  REAL8TimeSeries* cosbeta, /**< beta Euler angle time series */
104  REAL8TimeSeries* gam, /**< gamma Euler angle time series */
105  PhenomT_precomputed_sqrt *SQRT /**< precomputed squared root factors */
106 );
107 
108 /* Function to rotate a given set of Spin-Weighted Spherical Harmonic mode time series with an global Euler rotation specified by three Euler angles.*/
109 int PhenomTPHM_RotateModes_Global(
110  SphHarmTimeSeries* h_lm, /**< spherical harmonic decomposed modes, modified in place */
111  REAL8 alpha, /**< alpha Euler angle time series */
112  REAL8 cosbeta, /**< beta Euler angle time series */
113  REAL8 gam, /**< gamma Euler angle time series */
114  size_t length, /**< Length of waveform modes */
115  PhenomT_precomputed_sqrt *SQRT /**< precomputed squared root factors */
116 );
117 
118 /*********************************************/
119 /********* WIGNER-D MATRIX ROUTINES **********/
120 /*********************************************/
121 
122 /* Function to precompute and store squared root factors for the Wigner coefficients */
123 void IMRPhenomTPHM_SetPrecomputedSqrt(PhenomT_precomputed_sqrt *SQRT)
124 {
125  SQRT->sqrt2 = sqrt(2.0);
126  SQRT->sqrt2half = sqrt(2.5);
127  SQRT->sqrt3 = sqrt(3.0);
128  SQRT->sqrt5 = sqrt(5.0);
129  SQRT->sqrt6 = sqrt(6.0);
130  SQRT->sqrt7 = sqrt(7.0);
131  SQRT->sqrt10 = sqrt(10.0);
132  SQRT->sqrt14 = sqrt(14.0);
133  SQRT->sqrt15 = sqrt(15.0);
134  SQRT->sqrt21 = sqrt(21.0);
135  SQRT->sqrt30 = sqrt(30.0);
136  SQRT->sqrt35 = sqrt(35.0);
137  SQRT->sqrt70 = sqrt(70.0);
138  SQRT->sqrt210 = sqrt(210.0);
139 }
140 
141 /* Function to compute the WignerDMatrix. Analogous to XLALWignerDMatrix (which employs XLALWignerdMatrix), but Wignerd coefficients
142 here are read from the precomputed struct, in order to improve computational speed.
143 
144 Format of the Winerd elements corresponds to dlm[idx] with idx from l-m' to l+m' and the corresponding element is given
145 the function WignerD[{l,m',m},0,-beta,0] in Mathematica 12. See eq. 5 of https://dcc.ligo.org/DocDB/0159/T1900080/008/STmodes.pdf for
146 the relation between Mathematica's WignerD and XLALWignerDMatrix. */
147 void IMRPhenomTPHM_SetWignerDStruct(PhenomTPWignerStruct *wS, PhenomT_precomputed_sqrt *SQRT, REAL8 cosBeta, REAL8 alpha, REAL8 gamma, INT4 LMAX, INT4 sign, UINT4 globalRot)
148 {
149  /*cos(beta/2) and sin(beta/2) powers*/
150 
151  REAL8 cBetah = sqrt(0.5*fabs(1 + cosBeta));
152  REAL8 sBetah = sign*sqrt(0.5*fabs(1 - cosBeta));
153 
154  REAL8 cBetah2 = cBetah * cBetah;
155  REAL8 cBetah3 = cBetah * cBetah2;
156  REAL8 cBetah4 = cBetah * cBetah3;
157 
158  REAL8 sBetah2 = sBetah * sBetah;
159  REAL8 sBetah3 = sBetah * sBetah2;
160  REAL8 sBetah4 = sBetah * sBetah3;
161 
162  REAL8 C2mS2 = cBetah2 - sBetah2;
163 
164  REAL8 cBetah5, cBetah6, cBetah7, cBetah8, cBetah9, cBetah10;
165  REAL8 sBetah5, sBetah6, sBetah7, sBetah8, sBetah9, sBetah10;
166 
167  cBetah5 = 0.0; cBetah6 = 0.0; cBetah7 = 0.0; cBetah8 = 0.0; cBetah9 = 0.0; cBetah10 = 0.0;
168  sBetah5 = 0.0; sBetah6 = 0.0; sBetah7 = 0.0; sBetah8 = 0.0; sBetah9 = 0.0; sBetah10 = 0.0;
169 
170 
171  /* L=2 */
172 
173  REAL8 d22[5] = {sBetah4, 2.0*cBetah*sBetah3, SQRT->sqrt6*sBetah2*cBetah2, 2.0*cBetah3*sBetah, cBetah4};
174  REAL8 d2m2[5] = {d22[4], -d22[3], d22[2], -d22[1], d22[0]};
175 
176  REAL8 d21[5] = {2.0*cBetah*sBetah3, 3.0*cBetah2*sBetah2 - sBetah4, SQRT->sqrt6*(cBetah3*sBetah - cBetah*sBetah3), cBetah2*(cBetah2 - 3.0*sBetah2), -2.0*cBetah3*sBetah};
177  REAL8 d2m1[5] = {-d21[4], d21[3], -d21[2], d21[1], -d21[0]};
178 
179  for(UINT4 i=0; i<5; i++){
180  wS->wignerdL2[0][i] = d2m2[i];
181  wS->wignerdL2[1][i] = d2m1[i];
182  wS->wignerdL2[2][i] = 0.0;
183  wS->wignerdL2[3][i] = d21[i];
184  wS->wignerdL2[4][i] = d22[i];
185  }
186 
187  switch(LMAX)
188  {
189 
190  case 3:
191  {
192  cBetah5 = cBetah * cBetah4;
193  cBetah6 = cBetah * cBetah5;
194 
195  sBetah5 = sBetah * sBetah4;
196  sBetah6 = sBetah * sBetah5;
197 
198  REAL8 d33[7] = {sBetah6, SQRT->sqrt6*cBetah*sBetah5, SQRT->sqrt15*cBetah2*sBetah4, 2.0*SQRT->sqrt5*cBetah3*sBetah3, SQRT->sqrt15*cBetah4*sBetah2, SQRT->sqrt6*cBetah5*sBetah, cBetah6};
199  REAL8 d3m3[7] = {d33[6], -d33[5], d33[4], -d33[3], d33[2], -d33[1], d33[0]};
200 
201  for(UINT4 i=0; i<7; i++){
202  wS->wignerdL3[0][i] = d3m3[i];
203  wS->wignerdL3[1][i] = 0.0;
204  wS->wignerdL3[2][i] = 0.0;
205  wS->wignerdL3[3][i] = 0.0;
206  wS->wignerdL3[4][i] = 0.0;
207  wS->wignerdL3[5][i] = 0.0;
208  wS->wignerdL3[6][i] = d33[i];
209  }
210 
211  break;
212  }
213 
214 
215  case 4:
216  {
217  cBetah5 = cBetah * cBetah4;
218  cBetah6 = cBetah * cBetah5;
219  cBetah7 = cBetah * cBetah6;
220  cBetah8 = cBetah * cBetah7;
221 
222  sBetah5 = sBetah * sBetah4;
223  sBetah6 = sBetah * sBetah5;
224  sBetah7 = sBetah * sBetah6;
225  sBetah8 = sBetah * sBetah7;
226 
227  REAL8 d33[7] = {sBetah6, SQRT->sqrt6*cBetah*sBetah5, SQRT->sqrt15*cBetah2*sBetah4, 2.0*SQRT->sqrt5*cBetah3*sBetah3, SQRT->sqrt15*cBetah4*sBetah2, SQRT->sqrt6*cBetah5*sBetah, cBetah6};
228  REAL8 d3m3[7] = {d33[6], -d33[5], d33[4], -d33[3], d33[2], -d33[1], d33[0]};
229 
230  REAL8 d44[9] = {sBetah8, 2.0*SQRT->sqrt2*cBetah*sBetah7, 2.0*SQRT->sqrt7*cBetah2*sBetah6, 2.0*SQRT->sqrt14*cBetah3*sBetah5, SQRT->sqrt70*cBetah4*sBetah4, 2.0*SQRT->sqrt14*cBetah5*sBetah3, 2.0*SQRT->sqrt7*cBetah6*sBetah2, 2.0*sqrt(2)*cBetah7*sBetah, cBetah8};
231  REAL8 d4m4[9] = {d44[8], -d44[7], d44[6], -d44[5], d44[4], -d44[3], d44[2], -d44[1], d44[0]};
232 
233  for(UINT4 i=0; i<7; i++){
234  wS->wignerdL3[0][i] = d3m3[i];
235  wS->wignerdL3[1][i] = 0.0;
236  wS->wignerdL3[2][i] = 0.0;
237  wS->wignerdL3[3][i] = 0.0;
238  wS->wignerdL3[4][i] = 0.0;
239  wS->wignerdL3[5][i] = 0.0;
240  wS->wignerdL3[6][i] = d33[i];
241  }
242 
243  for(UINT4 i=0; i<9; i++){
244  wS->wignerdL4[0][i] = d4m4[i];
245  wS->wignerdL4[1][i] = 0.0;
246  wS->wignerdL4[2][i] = 0.0;
247  wS->wignerdL4[3][i] = 0.0;
248  wS->wignerdL4[4][i] = 0.0;
249  wS->wignerdL4[5][i] = 0.0;
250  wS->wignerdL4[6][i] = 0.0;
251  wS->wignerdL4[7][i] = 0.0;
252  wS->wignerdL4[8][i] = d44[i];
253  }
254 
255  break;
256  }
257 
258  case 5:
259  {
260  cBetah5 = cBetah * cBetah4;
261  cBetah6 = cBetah * cBetah5;
262  cBetah7 = cBetah * cBetah6;
263  cBetah8 = cBetah * cBetah7;
264  cBetah9 = cBetah * cBetah8;
265  cBetah10 = cBetah * cBetah9;
266 
267  sBetah5 = sBetah * sBetah4;
268  sBetah6 = sBetah * sBetah5;
269  sBetah7 = sBetah * sBetah6;
270  sBetah8 = sBetah * sBetah7;
271  sBetah9 = sBetah * sBetah8;
272  sBetah10 = sBetah * sBetah9;
273 
274  REAL8 d33[7] = {sBetah6, SQRT->sqrt6*cBetah*sBetah5, SQRT->sqrt15*cBetah2*sBetah4, 2.0*SQRT->sqrt5*cBetah3*sBetah3, SQRT->sqrt15*cBetah4*sBetah2, SQRT->sqrt6*cBetah5*sBetah, cBetah6};
275  REAL8 d3m3[7] = {d33[6], -d33[5], d33[4], -d33[3], d33[2], -d33[1], d33[0]};
276 
277  REAL8 d44[9] = {sBetah8, 2.0*SQRT->sqrt2*cBetah*sBetah7, 2.0*SQRT->sqrt7*cBetah2*sBetah6, 2.0*SQRT->sqrt14*cBetah3*sBetah5, SQRT->sqrt70*cBetah4*sBetah4, 2.0*SQRT->sqrt14*cBetah5*sBetah3, 2.0*SQRT->sqrt7*cBetah6*sBetah2, 2.0*sqrt(2)*cBetah7*sBetah, cBetah8};
278  REAL8 d4m4[9] = {d44[8], -d44[7], d44[6], -d44[5], d44[4], -d44[3], d44[2], -d44[1], d44[0]};
279 
280  REAL8 d55[11] = {sBetah10, SQRT->sqrt10*cBetah*sBetah9, 3*SQRT->sqrt5*cBetah2*sBetah8, 2*SQRT->sqrt30*cBetah3*sBetah7, SQRT->sqrt210*cBetah4*sBetah6, 6.0*SQRT->sqrt7*cBetah5*sBetah5,\
281  SQRT->sqrt210*cBetah6*sBetah4, 2*SQRT->sqrt30*cBetah7*sBetah3, 3*SQRT->sqrt5*cBetah8*sBetah2, SQRT->sqrt10*cBetah9*sBetah, cBetah10};
282  REAL8 d5m5[11] = {d55[10], -d55[9], d55[8], -d55[7], d55[6], -d55[5], d55[4], -d55[3], d55[2], -d55[1], d55[0]};
283 
284  for(UINT4 i=0; i<7; i++){
285  wS->wignerdL3[0][i] = d3m3[i];
286  wS->wignerdL3[1][i] = 0.0;
287  wS->wignerdL3[2][i] = 0.0;
288  wS->wignerdL3[3][i] = 0.0;
289  wS->wignerdL3[4][i] = 0.0;
290  wS->wignerdL3[5][i] = 0.0;
291  wS->wignerdL3[6][i] = d33[i];
292  }
293 
294  for(UINT4 i=0; i<9; i++){
295  wS->wignerdL4[0][i] = d4m4[i];
296  wS->wignerdL4[1][i] = 0.0;
297  wS->wignerdL4[2][i] = 0.0;
298  wS->wignerdL4[3][i] = 0.0;
299  wS->wignerdL4[4][i] = 0.0;
300  wS->wignerdL4[5][i] = 0.0;
301  wS->wignerdL4[6][i] = 0.0;
302  wS->wignerdL4[7][i] = 0.0;
303  wS->wignerdL4[8][i] = d44[i];
304  }
305 
306  for(UINT4 i=0; i<11; i++){
307  wS->wignerdL5[0][i] = d5m5[i];
308  wS->wignerdL5[1][i] = 0.0;
309  wS->wignerdL5[2][i] = 0.0;
310  wS->wignerdL5[3][i] = 0.0;
311  wS->wignerdL5[4][i] = 0.0;
312  wS->wignerdL5[5][i] = 0.0;
313  wS->wignerdL5[6][i] = 0.0;
314  wS->wignerdL5[7][i] = 0.0;
315  wS->wignerdL5[8][i] = 0.0;
316  wS->wignerdL5[9][i] = 0.0;
317  wS->wignerdL5[10][i] = d55[i];
318  }
319 
320  break;
321  }
322  }
323 
324 
325  /* For performing the global rotation between the J-frame and L0-frame, all coefficients for a given L have to be computed,
326  since the list of modes in the J-frame contains all the modes for a given L */
327  if(globalRot==1)
328  {
329  UNUSED REAL8 sinBeta = sign*sqrt(fabs(1.0 - cosBeta*cosBeta));
330 
331  UNUSED REAL8 cos2Beta = cosBeta*cosBeta - sinBeta*sinBeta;
332  UNUSED REAL8 cos3Beta = cosBeta*(2.0*cos2Beta - 1.0);
333  UNUSED REAL8 cos4Beta = pow(sinBeta,4) + pow(cosBeta,4) - 6.0*sinBeta*sinBeta*cosBeta*cosBeta;
334 
335  switch(LMAX)
336  {
337  case 2:
338  {
339  REAL8 d20[5] = {SQRT->sqrt6*cBetah2*sBetah2 , SQRT->sqrt6*cBetah*sBetah*C2mS2 , 0.25*(1 + 3*(-4*cBetah2*sBetah2 + pow(C2mS2,2))) , -SQRT->sqrt6*cBetah*sBetah*C2mS2 , SQRT->sqrt6*cBetah2*sBetah2};
340 
341  for(UINT4 i=0; i<5; i++){
342  wS->wignerdL2[2][i] = d20[i];
343  }
344  break;
345  }
346 
347  case 3:
348  {
349  REAL8 d20[5] = {SQRT->sqrt6*cBetah2*sBetah2 , SQRT->sqrt6*cBetah*sBetah*C2mS2 , 0.25*(1 + 3*(-4*cBetah2*sBetah2 + pow(C2mS2,2))) , -SQRT->sqrt6*cBetah*sBetah*C2mS2 , SQRT->sqrt6*cBetah2*sBetah2};
350 
351  REAL8 d32[7] = {SQRT->sqrt6*cBetah*sBetah5, sBetah4*(5.0*cBetah2 - sBetah2), SQRT->sqrt10*sBetah3*(2.0*cBetah3 - cBetah*sBetah2), SQRT->sqrt30*cBetah2*(cBetah2 - sBetah2)*sBetah2, SQRT->sqrt10*cBetah3*(cBetah2*sBetah - 2.0*sBetah3), cBetah4*(cBetah2 - 5.0*sBetah2), -SQRT->sqrt6*cBetah5*sBetah};
352  REAL8 d3m2[7] = {-d32[6],d32[5],-d32[4],d32[3],-d32[2],d32[1],-d32[0]};
353 
354  REAL8 d31[7] = {SQRT->sqrt15*cBetah2*sBetah4, SQRT->sqrt2half*cBetah*sBetah3*(1 + 3.0*C2mS2), 0.125*sBetah2*(13.0 + 20.0*C2mS2 + 15.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)), 0.25*SQRT->sqrt3*cBetah*sBetah*(3.0 + 5.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)),\
355  0.125*cBetah2*(13.0 - 20.0*C2mS2 + 15.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)), -SQRT->sqrt2half*cBetah3*sBetah*(-1.0 + 3.0*C2mS2), SQRT->sqrt15*cBetah4*sBetah2};
356  REAL8 d3m1[7] = {d31[6], -d31[5], d31[4], -d31[3], d31[2], -d31[1], d31[0]};
357 
358  REAL8 d30[7] = {2.0*SQRT->sqrt5*cBetah3*sBetah3, SQRT->sqrt30*cBetah2*sBetah2*C2mS2, 0.25*SQRT->sqrt3*cBetah*sBetah*(3.0 + 5.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)), 0.125*(5.0*cos3Beta + 3*C2mS2),\
359  -0.25*SQRT->sqrt3*cBetah*sBetah*(3.0 + 5.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)), SQRT->sqrt30*cBetah2*sBetah2*C2mS2, -2.0*SQRT->sqrt5*cBetah3*sBetah3};
360 
361  for(UINT4 i=0; i<5; i++){
362  wS->wignerdL2[2][i] = d20[i];
363  }
364  for(UINT4 i=0; i<7; i++){
365  wS->wignerdL3[1][i] = d3m2[i];
366  wS->wignerdL3[2][i] = d3m1[i];
367  wS->wignerdL3[3][i] = d30[i];
368  wS->wignerdL3[4][i] = d31[i];
369  wS->wignerdL3[5][i] = d32[i];
370  }
371  break;
372  }
373 
374  case 4:
375  {
376  REAL8 d20[5] = {SQRT->sqrt6*cBetah2*sBetah2 , SQRT->sqrt6*cBetah*sBetah*C2mS2 , 0.25*(1 + 3*(-4*cBetah2*sBetah2 + pow(C2mS2,2))) , -SQRT->sqrt6*cBetah*sBetah*C2mS2 , SQRT->sqrt6*cBetah2*sBetah2};
377 
378  REAL8 d32[7] = {SQRT->sqrt6*cBetah*sBetah5, sBetah4*(5.0*cBetah2 - sBetah2), SQRT->sqrt10*sBetah3*(2.0*cBetah3 - cBetah*sBetah2), SQRT->sqrt30*cBetah2*(cBetah2 - sBetah2)*sBetah2, SQRT->sqrt10*cBetah3*(cBetah2*sBetah - 2.0*sBetah3), cBetah4*(cBetah2 - 5.0*sBetah2), -SQRT->sqrt6*cBetah5*sBetah};
379  REAL8 d3m2[7] = {-d32[6],d32[5],-d32[4],d32[3],-d32[2],d32[1],-d32[0]};
380 
381  REAL8 d31[7] = {SQRT->sqrt15*cBetah2*sBetah4, SQRT->sqrt2half*cBetah*sBetah3*(1 + 3.0*C2mS2), 0.125*sBetah2*(13.0 + 20.0*C2mS2 + 15.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)), 0.25*SQRT->sqrt3*cBetah*sBetah*(3.0 + 5.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)),\
382  0.125*cBetah2*(13.0 - 20.0*C2mS2 + 15.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)), -SQRT->sqrt2half*cBetah3*sBetah*(-1.0 + 3.0*C2mS2), SQRT->sqrt15*cBetah4*sBetah2};
383  REAL8 d3m1[7] = {d31[6], -d31[5], d31[4], -d31[3], d31[2], -d31[1], d31[0]};
384 
385  REAL8 d30[7] = {2.0*SQRT->sqrt5*cBetah3*sBetah3, SQRT->sqrt30*cBetah2*sBetah2*C2mS2, 0.25*SQRT->sqrt3*cBetah*sBetah*(3.0 + 5.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)), 0.125*(5.0*cos3Beta + 3*C2mS2),\
386  -0.25*SQRT->sqrt3*cBetah*sBetah*(3.0 + 5.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)), SQRT->sqrt30*cBetah2*sBetah2*C2mS2, -2.0*SQRT->sqrt5*cBetah3*sBetah3};
387 
388  REAL8 d43[9] = {2*SQRT->sqrt2*cBetah*sBetah7, 7*cBetah2*sBetah6-sBetah8, SQRT->sqrt14*(3*cBetah3*sBetah5-cBetah*sBetah7), SQRT->sqrt7*(5*cBetah4*sBetah4-3*cBetah2*sBetah6),\
389  2*5.916079783099616*(cBetah5*sBetah3-cBetah3*sBetah5), SQRT->sqrt7*(3*cBetah6*sBetah2-5*cBetah4*sBetah4), SQRT->sqrt14*(cBetah7*sBetah-3*cBetah5*sBetah3), cBetah8-7*cBetah6*sBetah2, -2.*SQRT->sqrt2*cBetah7*sBetah};
390  REAL8 d4m3[9] = {-d43[8], d43[7], -d43[6], d43[5],-d43[4],d43[3], -d43[2], d43[1], -d43[0]};
391 
392  REAL8 d42[9] = {2*SQRT->sqrt7*cBetah2*sBetah6, SQRT->sqrt14*cBetah*sBetah5*(1.0 +2.0*C2mS2), sBetah4*(1.0 + 7.0*C2mS2 + 7.0*C2mS2*C2mS2), 0.5*SQRT->sqrt2*cBetah*sBetah3*(6.0 + 7.0*cos2Beta + 7.0*C2mS2),\
393  0.5*SQRT->sqrt2half*cBetah2*(5.0 + 7.0*cos2Beta)*sBetah2, 0.5*SQRT->sqrt2*cBetah3*sBetah*(6.0 + 7.0*cos2Beta - 7.0*C2mS2), cBetah4*(1.0 - 7.0*C2mS2 + 7.0*C2mS2*C2mS2), -SQRT->sqrt14*cBetah5*sBetah*(-1.0 +2.0*C2mS2), 2*SQRT->sqrt7*cBetah6*sBetah2};
394  REAL8 d4m2[9] = {d42[8], -d42[7], d42[6], -d42[5], d42[4], -d42[3], d42[2], -d42[1], d42[0]};
395 
396  REAL8 d41[9] = {2*SQRT->sqrt14*cBetah3*sBetah5, SQRT->sqrt7*cBetah2*sBetah4*(1.0 + 4.0*C2mS2), 0.5*SQRT->sqrt2*cBetah*sBetah3*(6.0 + 7.0*cos2Beta + 7.0*C2mS2), 0.125*sBetah2*(15.0 +21.0*cos2Beta + 14.0*cos3Beta + 30.0*C2mS2),\
397  0.125*SQRT->sqrt5*cBetah*sBetah*(7.0*cos3Beta + 9.0*C2mS2), 0.125*cBetah2*(-15.0 + 30*cosBeta - 21.0*cos2Beta + 14.0*cos3Beta), 0.5*SQRT->sqrt2*cBetah3*sBetah*(-6.0 - 7.0*cos2Beta + 7.0*C2mS2),\
398  SQRT->sqrt7*cBetah4*sBetah2*(-1.0 + 4.0*C2mS2), -2*SQRT->sqrt14*cBetah5*sBetah3};
399  REAL8 d4m1[9] = {-d41[8], d41[7], -d41[6], d41[5], -d41[4], d41[3], -d41[2], d41[1], -d41[0]};
400 
401  REAL8 d40[9] = {SQRT->sqrt70*cBetah4*sBetah4, 2*SQRT->sqrt35*cBetah3*sBetah3*C2mS2, 0.5*SQRT->sqrt2half*cBetah2*(5. + 7.*cos2Beta)*sBetah2, 0.125*SQRT->sqrt5*cBetah*sBetah*(7.*cos3Beta + 9.*C2mS2),\
402  0.015625*(9 + 20.*cos2Beta + 35.*cos4Beta), -0.125*SQRT->sqrt5*cBetah*sBetah*(7.*cos3Beta + 9.*C2mS2), 0.5*SQRT->sqrt2half*cBetah2*(5. + 7.*cos2Beta)*sBetah2, -2.*SQRT->sqrt35*cBetah3*sBetah3*C2mS2, SQRT->sqrt70*cBetah4*sBetah4};
403 
404  for(UINT4 i=0; i<5; i++){
405  wS->wignerdL2[2][i] = d20[i];
406  }
407  for(UINT4 i=0; i<7; i++){
408  wS->wignerdL3[1][i] = d3m2[i];
409  wS->wignerdL3[2][i] = d3m1[i];
410  wS->wignerdL3[3][i] = d30[i];
411  wS->wignerdL3[4][i] = d31[i];
412  wS->wignerdL3[5][i] = d32[i];
413  }
414  for(UINT4 i=0; i<9; i++){
415  wS->wignerdL4[1][i] = d4m3[i];
416  wS->wignerdL4[2][i] = d4m2[i];
417  wS->wignerdL4[3][i] = d4m1[i];
418  wS->wignerdL4[4][i] = d40[i];
419  wS->wignerdL4[5][i] = d41[i];
420  wS->wignerdL4[6][i] = d42[i];
421  wS->wignerdL4[7][i] = d43[i];
422  }
423  break;
424  }
425 
426  case 5:
427  {
428  REAL8 d20[5] = {SQRT->sqrt6*cBetah2*sBetah2 , SQRT->sqrt6*cBetah*sBetah*C2mS2 , 0.25*(1 + 3*(-4*cBetah2*sBetah2 + pow(C2mS2,2))) , -SQRT->sqrt6*cBetah*sBetah*C2mS2 , SQRT->sqrt6*cBetah2*sBetah2};
429 
430  REAL8 d32[7] = {SQRT->sqrt6*cBetah*sBetah5, sBetah4*(5.0*cBetah2 - sBetah2), SQRT->sqrt10*sBetah3*(2.0*cBetah3 - cBetah*sBetah2), SQRT->sqrt30*cBetah2*(cBetah2 - sBetah2)*sBetah2, SQRT->sqrt10*cBetah3*(cBetah2*sBetah - 2.0*sBetah3), cBetah4*(cBetah2 - 5.0*sBetah2), -SQRT->sqrt6*cBetah5*sBetah};
431  REAL8 d3m2[7] = {-d32[6],d32[5],-d32[4],d32[3],-d32[2],d32[1],-d32[0]};
432 
433  REAL8 d31[7] = {SQRT->sqrt15*cBetah2*sBetah4, SQRT->sqrt2half*cBetah*sBetah3*(1 + 3.0*C2mS2), 0.125*sBetah2*(13.0 + 20.0*C2mS2 + 15.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)), 0.25*SQRT->sqrt3*cBetah*sBetah*(3.0 + 5.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)),\
434  0.125*cBetah2*(13.0 - 20.0*C2mS2 + 15.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)), -SQRT->sqrt2half*cBetah3*sBetah*(-1.0 + 3.0*C2mS2), SQRT->sqrt15*cBetah4*sBetah2};
435  REAL8 d3m1[7] = {d31[6], -d31[5], d31[4], -d31[3], d31[2], -d31[1], d31[0]};
436 
437  REAL8 d30[7] = {2.0*SQRT->sqrt5*cBetah3*sBetah3, SQRT->sqrt30*cBetah2*sBetah2*C2mS2, 0.25*SQRT->sqrt3*cBetah*sBetah*(3.0 + 5.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)), 0.125*(5.0*cos3Beta + 3*C2mS2),\
438  -0.25*SQRT->sqrt3*cBetah*sBetah*(3.0 + 5.0*(C2mS2*C2mS2 - 4.0*cBetah2*sBetah2)), SQRT->sqrt30*cBetah2*sBetah2*C2mS2, -2.0*SQRT->sqrt5*cBetah3*sBetah3};
439 
440  REAL8 d43[9] = {2*SQRT->sqrt2*cBetah*sBetah7, 7*cBetah2*sBetah6-sBetah8, SQRT->sqrt14*(3*cBetah3*sBetah5-cBetah*sBetah7), SQRT->sqrt7*(5*cBetah4*sBetah4-3*cBetah2*sBetah6),\
441  2*5.916079783099616*(cBetah5*sBetah3-cBetah3*sBetah5), SQRT->sqrt7*(3*cBetah6*sBetah2-5*cBetah4*sBetah4), SQRT->sqrt14*(cBetah7*sBetah-3*cBetah5*sBetah3), cBetah8-7*cBetah6*sBetah2, -2.*SQRT->sqrt2*cBetah7*sBetah};
442  REAL8 d4m3[9] = {-d43[8], d43[7], -d43[6], d43[5],-d43[4],d43[3], -d43[2], d43[1], -d43[0]};
443 
444  REAL8 d42[9] = {2*SQRT->sqrt7*cBetah2*sBetah6, SQRT->sqrt14*cBetah*sBetah5*(1.0 +2.0*C2mS2), sBetah4*(1.0 + 7.0*C2mS2 + 7.0*C2mS2*C2mS2), 0.5*SQRT->sqrt2*cBetah*sBetah3*(6.0 + 7.0*cos2Beta + 7.0*C2mS2),\
445  0.5*SQRT->sqrt2half*cBetah2*(5.0 + 7.0*cos2Beta)*sBetah2, 0.5*SQRT->sqrt2*cBetah3*sBetah*(6.0 + 7.0*cos2Beta - 7.0*C2mS2), cBetah4*(1.0 - 7.0*C2mS2 + 7.0*C2mS2*C2mS2), -SQRT->sqrt14*cBetah5*sBetah*(-1.0 +2.0*C2mS2), 2*SQRT->sqrt7*cBetah6*sBetah2};
446  REAL8 d4m2[9] = {d42[8], -d42[7], d42[6], -d42[5], d42[4], -d42[3], d42[2], -d42[1], d42[0]};
447 
448  REAL8 d41[9] = {2*SQRT->sqrt14*cBetah3*sBetah5, SQRT->sqrt7*cBetah2*sBetah4*(1.0 + 4.0*C2mS2), 0.5*SQRT->sqrt2*cBetah*sBetah3*(6.0 + 7.0*cos2Beta + 7.0*C2mS2), 0.125*sBetah2*(15.0 +21.0*cos2Beta + 14.0*cos3Beta + 30.0*C2mS2),\
449  0.125*SQRT->sqrt5*cBetah*sBetah*(7.0*cos3Beta + 9.0*C2mS2), 0.125*cBetah2*(-15.0 + 30*cosBeta - 21.0*cos2Beta + 14.0*cos3Beta), 0.5*SQRT->sqrt2*cBetah3*sBetah*(-6.0 - 7.0*cos2Beta + 7.0*C2mS2),\
450  SQRT->sqrt7*cBetah4*sBetah2*(-1.0 + 4.0*C2mS2), -2*SQRT->sqrt14*cBetah5*sBetah3};
451  REAL8 d4m1[9] = {-d41[8], d41[7], -d41[6], d41[5], -d41[4], d41[3], -d41[2], d41[1], -d41[0]};
452 
453  REAL8 d40[9] = {SQRT->sqrt70*cBetah4*sBetah4, 2*SQRT->sqrt35*cBetah3*sBetah3*C2mS2, 0.5*SQRT->sqrt2half*cBetah2*(5. + 7.*cos2Beta)*sBetah2, 0.125*SQRT->sqrt5*cBetah*sBetah*(7.*cos3Beta + 9.*C2mS2),\
454  0.015625*(9 + 20.*cos2Beta + 35.*cos4Beta), -0.125*SQRT->sqrt5*cBetah*sBetah*(7.*cos3Beta + 9.*C2mS2), 0.5*SQRT->sqrt2half*cBetah2*(5. + 7.*cos2Beta)*sBetah2, -2.*SQRT->sqrt35*cBetah3*sBetah3*C2mS2, SQRT->sqrt70*cBetah4*sBetah4};
455 
456  REAL8 d54[11] = {SQRT->sqrt10*cBetah*sBetah9, sBetah8*(4.0 + 5.0*C2mS2), (3./sqrt(2))*cBetah*sBetah7*(3.0 +5.0*C2mS2), 2.0*SQRT->sqrt3*cBetah2*sBetah6*(2.0 + 5.0*C2mS2), SQRT->sqrt21*cBetah3*sBetah5*(1.0 + 5.0*C2mS2),\
457  3.0*SQRT->sqrt70*cBetah4*sBetah4*C2mS2, SQRT->sqrt21*cBetah5*sBetah3*(-1.0 + 5.0*C2mS2), 2.0*SQRT->sqrt3*cBetah6*sBetah2*(-2.0 + 5.0*C2mS2), (3./sqrt(2))*cBetah7*sBetah*(-3.0 +5.0*C2mS2), cBetah8*(-4.0 + 5.0*C2mS2), -SQRT->sqrt10*cBetah9*sBetah};
458  REAL8 d5m4[11] = {-d54[10], d54[9], -d54[8], d54[7], -d54[6], d54[5], -d54[4], d54[3], -d54[2], d54[1], -d54[0]};
459 
460  REAL8 d53[11] = {3.0*SQRT->sqrt5*cBetah2*sBetah8, (3.0/SQRT->sqrt2)*cBetah*sBetah7*(3.0+5.0*C2mS2), 0.25*(13.0 + 54.0*C2mS2 + 45.0*C2mS2*C2mS2)*sBetah6, sqrt(1.5)*(1.0+12.0*C2mS2+15.0*C2mS2*C2mS2)*cBetah*sBetah5, 0.5*sqrt(10.5)*(-1.0+6.0*C2mS2+15.0*C2mS2*C2mS2)*cBetah2*sBetah4,\
461  0.25*SQRT->sqrt35*(7.0 + 9.0*cos2Beta)*cBetah3*sBetah3, 0.5*sqrt(10.5)*(-1.0-6.0*C2mS2+15.0*C2mS2*C2mS2)*cBetah4*sBetah2, sqrt(1.5)*(1.0-12.0*C2mS2+15.0*C2mS2*C2mS2)*cBetah5*sBetah, 0.25*(13.0 - 54.0*C2mS2 + 45.0*C2mS2*C2mS2)*cBetah6, (3.0/SQRT->sqrt2)*cBetah7*sBetah*(3.0-5.0*C2mS2), 3.0*SQRT->sqrt5*cBetah8*sBetah2};
462  REAL8 d5m3[11] = {d53[10], -d53[9], d53[8], -d53[7], d53[6], -d53[5], d53[4], -d53[3], d53[2], -d53[1], d53[0]};
463 
464  REAL8 d52[11] = {2*SQRT->sqrt30*cBetah3*sBetah7, 2.0*SQRT->sqrt3*(2.0+5.0*C2mS2)*cBetah2*sBetah6, sqrt(1.5)*(1.0+12.0*C2mS2+15.0*C2mS2*C2mS2)*cBetah*sBetah5, (-1.0+3.0*C2mS2+18.0*C2mS2*C2mS2+15.0*C2mS2*C2mS2*C2mS2)*sBetah4,\
465  0.5*SQRT->sqrt7*(-1.0-3.0*C2mS2+9.0*C2mS2*C2mS2+15.0*C2mS2*C2mS2*C2mS2)*cBetah*sBetah3, 0.5*sqrt(52.5)*C2mS2*cBetah2*sBetah2*(1.0+3.0*cos2Beta), 0.5*SQRT->sqrt7*(1.0-3.0*C2mS2-9.0*C2mS2*C2mS2+15.0*C2mS2*C2mS2*C2mS2)*cBetah3*sBetah,\
466  (1.0+3.0*C2mS2-18.0*C2mS2*C2mS2+15.0*C2mS2*C2mS2*C2mS2)*cBetah4, -sqrt(1.5)*(1.0-12.0*C2mS2+15.0*C2mS2*C2mS2)*cBetah5*sBetah, 2.0*SQRT->sqrt3*(-2.0+5.0*C2mS2)*cBetah6*sBetah2, -2*SQRT->sqrt30*cBetah7*sBetah3};
467  REAL8 d5m2[11] = {-d52[10], d52[9], -d52[8], d52[7], -d52[6], d52[5], -d52[4], d52[3], -d52[2], d52[1], -d52[0]};
468 
469  REAL8 d51[11] = {SQRT->sqrt210*cBetah4*sBetah6, SQRT->sqrt21*(1.0+5.0*C2mS2)*cBetah3*sBetah5, 0.5*sqrt(10.5)*(-1.0+6.0*C2mS2+15.0*C2mS2*C2mS2)*cBetah2*sBetah4,\
470  0.5*SQRT->sqrt7*(-1.0-3.0*C2mS2+9.0*C2mS2*C2mS2+15.0*C2mS2*C2mS2*C2mS2)*cBetah*sBetah3, 0.125*(1.0-28.0*C2mS2-42.0*C2mS2*C2mS2+84.0*C2mS2*C2mS2*C2mS2+105.0*C2mS2*C2mS2*C2mS2*C2mS2)*sBetah2,\
471  sqrt(7.5)/32.0*cBetah*sBetah*(15.0 + 28.0*cos2Beta + 21.0*cos4Beta), 0.125*(1.0+28.0*C2mS2-42.0*C2mS2*C2mS2-84.0*C2mS2*C2mS2*C2mS2+105.0*C2mS2*C2mS2*C2mS2*C2mS2)*cBetah2,\
472  -0.5*SQRT->sqrt7*(1.0-3.0*C2mS2-9.0*C2mS2*C2mS2+15.0*C2mS2*C2mS2*C2mS2)*cBetah3*sBetah, 0.5*sqrt(10.5)*(-1.0-6.0*C2mS2+15.0*C2mS2*C2mS2)*cBetah4*sBetah2,\
473  -SQRT->sqrt21*(-1.0+5.0*C2mS2)*cBetah5*sBetah3, SQRT->sqrt210*cBetah6*sBetah4};
474  REAL8 d5m1[11] = {d51[10], -d51[9], d51[8], -d51[7], d51[6], -d51[5], d51[4], -d51[3], d51[2], -d51[1], d51[0]};
475 
476  REAL8 d50[11] = {6.0*SQRT->sqrt7*cBetah5*sBetah5, 3.0*SQRT->sqrt70*C2mS2*cBetah4*sBetah4, 0.25*SQRT->sqrt35*cBetah3*sBetah3*(7.0+9.0*cos2Beta), 0.5*sqrt(52.2)*C2mS2*cBetah2*sBetah2*(1.0+3.0*cos2Beta),\
477  sqrt(7.5)/32.0*cBetah*sBetah*(15.0+28.0*cos2Beta+21.0*cos4Beta), 0.125*C2mS2*(15.0-70.0*C2mS2*C2mS2+63.0*C2mS2*C2mS2*C2mS2*C2mS2), -sqrt(7.5)/32.0*cBetah*sBetah*(15.0+28.0*cos2Beta+21.0*cos4Beta),\
478  0.5*sqrt(52.2)*C2mS2*cBetah2*sBetah2*(1.0+3.0*cos2Beta), -0.25*SQRT->sqrt35*cBetah3*sBetah3*(7.0+9.0*cos2Beta), 3.0*SQRT->sqrt70*C2mS2*cBetah4*sBetah4, -6.0*SQRT->sqrt7*cBetah5*sBetah5};
479 
480  for(UINT4 i=0; i<5; i++){
481  wS->wignerdL2[2][i] = d20[i];
482  }
483  for(UINT4 i=0; i<7; i++){
484  wS->wignerdL3[1][i] = d3m2[i];
485  wS->wignerdL3[2][i] = d3m1[i];
486  wS->wignerdL3[3][i] = d30[i];
487  wS->wignerdL3[4][i] = d31[i];
488  wS->wignerdL3[5][i] = d32[i];
489  }
490  for(UINT4 i=0; i<9; i++){
491  wS->wignerdL4[1][i] = d4m3[i];
492  wS->wignerdL4[2][i] = d4m2[i];
493  wS->wignerdL4[3][i] = d4m1[i];
494  wS->wignerdL4[4][i] = d40[i];
495  wS->wignerdL4[5][i] = d41[i];
496  wS->wignerdL4[6][i] = d42[i];
497  wS->wignerdL4[7][i] = d43[i];
498  }
499  for(UINT4 i=0; i<11; i++){
500  wS->wignerdL5[1][i] = d5m4[i];
501  wS->wignerdL5[2][i] = d5m3[i];
502  wS->wignerdL5[3][i] = d5m2[i];
503  wS->wignerdL5[4][i] = d5m1[i];
504  wS->wignerdL5[5][i] = d50[i];
505  wS->wignerdL5[6][i] = d51[i];
506  wS->wignerdL5[7][i] = d52[i];
507  wS->wignerdL5[8][i] = d53[i];
508  wS->wignerdL5[9][i] = d54[i];
509  }
510  break;
511  }
512  }
513  }
514 
515  /* expAlpha powers */
516 
517  COMPLEX16 expAlpha = cexp(I*alpha);
518  COMPLEX16 expAlpham = 1./expAlpha;
519 
520  wS->expAlphaPowers[4] = expAlpha;
521  wS->expAlphaPowers[3] = expAlpha*expAlpha;
522  wS->expAlphaPowers[2] = expAlpha*wS->expAlphaPowers[3];
523  wS->expAlphaPowers[1] = expAlpha*wS->expAlphaPowers[2];
524  wS->expAlphaPowers[0] = expAlpha*wS->expAlphaPowers[1];
525  wS->expAlphaPowers[5] = 1.0;
526  wS->expAlphaPowers[6] = expAlpham;
527  wS->expAlphaPowers[7] = expAlpham*expAlpham;
528  wS->expAlphaPowers[8] = expAlpham*wS->expAlphaPowers[7];
529  wS->expAlphaPowers[9] = expAlpham*wS->expAlphaPowers[8];
530  wS->expAlphaPowers[10] = expAlpham*wS->expAlphaPowers[9];
531 
532  /* expGamma powers */
533 
534  COMPLEX16 expGamma = cexp(I*gamma);
535  COMPLEX16 expGammam = 1./expGamma;
536 
537  wS->expGammaPowers[4] = expGamma;
538  wS->expGammaPowers[3] = expGamma*expGamma;
539  wS->expGammaPowers[2] = expGamma*wS->expGammaPowers[3];
540  wS->expGammaPowers[1] = expGamma*wS->expGammaPowers[2];
541  wS->expGammaPowers[0] = expGamma*wS->expGammaPowers[1];
542  wS->expGammaPowers[5] = 1.0;
543  wS->expGammaPowers[6] = expGammam;
544  wS->expGammaPowers[7] = expGammam*expGammam;
545  wS->expGammaPowers[8] = expGammam*wS->expGammaPowers[7];
546  wS->expGammaPowers[9] = expGammam*wS->expGammaPowers[8];
547  wS->expGammaPowers[10] = expGammam*wS->expGammaPowers[9];
548 }
549 
550 /* Function to compute the WignerDMatrix. Analogous to XLALWignerDMatrix (which employs XLALWignerdMatrix), but Wignerd coefficients
551 here are read from the precomputed struct, in order to improve computational speed. */
552 COMPLEX16 PhenomTWignerDMatrix(
553  int l, /**< mode number l */
554  int mp, /**< mode number m' */
555  int m, /**< mode number m */
556  PhenomTPWignerStruct *wS /**< Wigner Struct, which contains Wignerd elements and powers of alpha and gamma */
557  )
558 {
559  COMPLEX16 wignerd = 0.0;
560 
561  switch(l) // For a given l, m, m' it returns the corresponding Wigner-d matrix coefficient, precomputed in the Wigner Struct
562  {
563  case 2:
564  wignerd = (COMPLEX16)wS->wignerdL2[2+mp][2+m];
565  break;
566  case 3:
567  wignerd = (COMPLEX16)wS->wignerdL3[3+mp][3+m];
568  break;
569  case 4:
570  wignerd = (COMPLEX16)wS->wignerdL4[4+mp][4+m];
571  break;
572  case 5:
573  wignerd = (COMPLEX16)wS->wignerdL5[5+mp][5+m];
574  break;
575 
576  }
577 
578  // Construct the Wigner-D matrix coefficient adding the exponentials of alpha and gamma for the corresponding l, m, m'
579  COMPLEX16 WignerD = wS->expAlphaPowers[mp+5]*wignerd*wS->expGammaPowers[m+5];
580 
581  return WignerD;
582 }
583 
584 /* Function to rotate a given set of Spin-Weighted Spherical Harmonic mode time series with an time-dependent Euler rotation specified by three Euler angles.
585 Adaptation of XLALSimInspiralPrecessionRotateModes in SimInspiralPrecess.c , only difference is that at each time step the WignerDStruct is set and PhenomTWignerDMatrix
586 is then employed. */
587 int PhenomTPHM_RotateModes(
588  SphHarmTimeSeries* h_lm, /**< spherical harmonic decomposed modes, modified in place */
589  REAL8TimeSeries* alpha, /**< alpha Euler angle time series */
590  REAL8TimeSeries* cosbeta, /**< beta Euler angle time series */
591  REAL8TimeSeries* gam, /**< gamma Euler angle time series */
592  PhenomT_precomputed_sqrt *SQRT /**< precomputed squared root factors */
593 ){
594 
595  unsigned int i;
596  int l, lmax, m, mp;
597  lmax = XLALSphHarmTimeSeriesGetMaxL( h_lm );
598 
599  // Temporary holding variables
600  complex double *x_lm = XLALCalloc( 2*lmax+1, sizeof(complex double) );
601  COMPLEX16TimeSeries **h_xx = XLALCalloc( 2*lmax+1, sizeof(COMPLEX16TimeSeries) );
602 
603  for(i=0; i<alpha->data->length; i++){
604  PhenomTPWignerStruct *wStruct;
605  wStruct = XLALMalloc(sizeof(PhenomTPWignerStruct));
606  IMRPhenomTPHM_SetWignerDStruct(wStruct, SQRT, cosbeta->data->data[i], alpha->data->data[i], gam->data->data[i], lmax, 1, 0);
607  for(l=2; l<=lmax; l++){
608  for(m=0; m<2*l+1; m++){
609  h_xx[m] = XLALSphHarmTimeSeriesGetMode(h_lm, l, m-l);
610  if( !h_xx[m] ){
611  x_lm[m] = 0;
612  } else {
613  x_lm[m] = h_xx[m]->data->data[i];
614  h_xx[m]->data->data[i] = 0;
615  }
616  }
617 
618  for(m=0; m<2*l+1; m++){
619  for(mp=0; mp<2*l+1; mp++){
620  if( !h_xx[m] ) continue;
621  if(!(creal(h_xx[m]->data->data[i])==0 && creal(x_lm[mp])==0)) {
622  h_xx[m]->data->data[i] +=
623  x_lm[mp] * PhenomTWignerDMatrix( l, mp-l, m-l, wStruct );
624  }
625  }
626  }
627  }
628  LALFree(wStruct);
629  }
630 
631 
632  XLALFree( x_lm );
633  XLALFree( h_xx );
634  return XLAL_SUCCESS;
635 }
636 
637 /* Function to rotate a given set of Spin-Weighted Spherical Harmonic mode time series with an time-dependent Euler rotation specified by three Euler angles.
638 Main difference with previous function is that for a fixed rotation, Wigner-D matrix only have to be computed once. */
639 int PhenomTPHM_RotateModes_Global(
640  SphHarmTimeSeries* h_lm, /**< spherical harmonic decomposed modes, modified in place */
641  REAL8 alpha, /**< alpha Euler angle time series */
642  REAL8 cosbeta, /**< beta Euler angle time series */
643  REAL8 gam, /**< gamma Euler angle time series */
644  size_t length, /**< Length of waveform modes */
645  PhenomT_precomputed_sqrt *SQRT /**< precomputed squared root factors */
646 ){
647 
648  unsigned int i;
649  int l, lmax, m, mp;
650  lmax = XLALSphHarmTimeSeriesGetMaxL( h_lm );
651  // Temporary holding variables
652  complex double *x_lm = XLALCalloc( 2*lmax+1, sizeof(complex double) );
653  COMPLEX16TimeSeries **h_xx = XLALCalloc( 2*lmax+1, sizeof(COMPLEX16TimeSeries) );
654 
655  PhenomTPWignerStruct *wStruct;
656  wStruct = XLALMalloc(sizeof(PhenomTPWignerStruct));
657  IMRPhenomTPHM_SetWignerDStruct(wStruct, SQRT, cosbeta, alpha, gam, lmax, -1, 1);
658 
659  COMPLEX16 wigner[lmax-1][2*lmax+1][2*lmax+1];
660 
661  for(l=2; l<=lmax; l++){
662  for(m=0; m<2*l+1; m++){
663  for(mp=0; mp<2*l+1; mp++){
664  wigner[l-2][mp][m] = PhenomTWignerDMatrix( l, mp-l, m-l, wStruct );
665  }
666  }
667  }
668 
669  for(i=0; i<length; i++){
670  for(l=2; l<=lmax; l++){
671  for(m=0; m<2*l+1; m++){
672  h_xx[m] = XLALSphHarmTimeSeriesGetMode(h_lm, l, m-l);
673  if( !h_xx[m] ){
674  x_lm[m] = 0;
675  } else {
676  x_lm[m] = h_xx[m]->data->data[i];
677  h_xx[m]->data->data[i] = 0;
678  }
679  }
680 
681  for(m=0; m<2*l+1; m++){
682  for(mp=0; mp<2*l+1; mp++){
683  if( !h_xx[m] ) continue;
684  if(!(creal(h_xx[m]->data->data[i])==0 && creal(x_lm[mp])==0)) {
685  h_xx[m]->data->data[i] += x_lm[mp] * wigner[l-2][mp][m];
686  }
687  }
688  }
689  }
690  }
691  XLALFree(wStruct);
692  XLALFree( x_lm );
693  XLALFree( h_xx );
694  return XLAL_SUCCESS;
695 }
LALFree
#define LALFree(p)
PhenomTPHM_RotateModes_Global
int PhenomTPHM_RotateModes_Global(SphHarmTimeSeries *h_lm, REAL8 alpha, REAL8 cosbeta, REAL8 gam, size_t length, PhenomT_precomputed_sqrt *SQRT)
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:639
IMRPhenomTPHM_SetPrecomputedSqrt
void IMRPhenomTPHM_SetPrecomputedSqrt(PhenomT_precomputed_sqrt *SQRT)
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:123
PhenomTPHM_RotateModes
int PhenomTPHM_RotateModes(SphHarmTimeSeries *h_lm, REAL8TimeSeries *alpha, REAL8TimeSeries *cosbeta, REAL8TimeSeries *gam, PhenomT_precomputed_sqrt *SQRT)
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:587
PhenomTWignerDMatrix
COMPLEX16 PhenomTWignerDMatrix(int l, int mp, int m, PhenomTPWignerStruct *wS)
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:552
IMRPhenomTPHM_SetWignerDStruct
void IMRPhenomTPHM_SetWignerDStruct(PhenomTPWignerStruct *wS, PhenomT_precomputed_sqrt *SQRT, REAL8 cosBeta, REAL8 alpha, REAL8 gamma, INT4 LMAX, INT4 sign, UINT4 globalRot)
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:147
l
int l
Definition: bh_qnmode.c:135
i
double i
Definition: bh_ringdown.c:118
data
sigmaKerr data[0]
Definition: exact_derivatives-HrealHeader.c:4
COMPLEX16
double complex COMPLEX16
REAL8
double REAL8
UINT4
uint32_t UINT4
INT4
int32_t INT4
XLALCalloc
void * XLALCalloc(size_t m, size_t n)
XLALMalloc
void * XLALMalloc(size_t n)
XLALFree
void XLALFree(void *p)
XLALSphHarmTimeSeriesGetMode
COMPLEX16TimeSeries * XLALSphHarmTimeSeriesGetMode(SphHarmTimeSeries *ts, UINT4 l, INT4 m)
Get the time series of a waveform's (l,m) spherical harmonic mode from a SphHarmTimeSeries linked lis...
Definition: LALSimSphHarmSeries.c:145
XLALSphHarmTimeSeriesGetMaxL
UINT4 XLALSphHarmTimeSeriesGetMaxL(SphHarmTimeSeries *ts)
Get the largest l index of any mode in the SphHarmTimeSeries linked list.
Definition: LALSimSphHarmSeries.c:164
m
static const INT4 m
XLAL_SUCCESS
XLAL_SUCCESS
alpha
double alpha
Definition: sgwb.c:106
COMPLEX16TimeSeries::data
COMPLEX16Sequence * data
COMPLEX16Vector::data
COMPLEX16 * data
PhenomT_precomputed_sqrt
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:58
PhenomT_precomputed_sqrt::sqrt5
REAL8 sqrt5
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
PhenomT_precomputed_sqrt::sqrt7
REAL8 sqrt7
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
PhenomT_precomputed_sqrt::sqrt6
REAL8 sqrt6
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
PhenomT_precomputed_sqrt::sqrt210
REAL8 sqrt210
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
PhenomT_precomputed_sqrt::sqrt2half
REAL8 sqrt2half
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
PhenomT_precomputed_sqrt::sqrt2
REAL8 sqrt2
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
PhenomT_precomputed_sqrt::sqrt14
REAL8 sqrt14
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
PhenomT_precomputed_sqrt::sqrt15
REAL8 sqrt15
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
PhenomT_precomputed_sqrt::sqrt35
REAL8 sqrt35
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
PhenomT_precomputed_sqrt::sqrt10
REAL8 sqrt10
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
PhenomT_precomputed_sqrt::sqrt70
REAL8 sqrt70
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
PhenomT_precomputed_sqrt::sqrt30
REAL8 sqrt30
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
PhenomT_precomputed_sqrt::sqrt3
REAL8 sqrt3
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
PhenomT_precomputed_sqrt::sqrt21
REAL8 sqrt21
Definition: LALSimIMRPhenomTPHM_FrameRotations.c:59
REAL8TimeSeries::data
REAL8Sequence * data
REAL8Vector::data
REAL8 * data
SphHarmTimeSeries
Structure to carry a collection of spherical harmonic modes in COMPLEX16 time series.
Definition: LALSimSphHarmSeries.h:50