LAL  7.5.0.1-ec27e42
SphericalHarmonics.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2007 S.Fairhurst, B. Krishnan, L.Santamaria, C. Robinson,
3  * C. Pankow
4  *
5  * This program is free software; you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation; either version 2 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with with program; see the file COPYING. If not, write to the
17  * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
18  * MA 02110-1301 USA
19  */
20 
21 #include <lal/SphericalHarmonics.h>
22 #include <lal/LALError.h>
23 #include <lal/XLALGSL.h>
24 
25 #include <gsl/gsl_sf_legendre.h>
26 #include <gsl/gsl_sf_gamma.h>
27 
28 /**
29  * Computes the (s)Y(l,m) spin-weighted spherical harmonic.
30  *
31  * From somewhere ....
32  *
33  * See also:
34  * Implements Equations (II.9)-(II.13) of
35  * D. A. Brown, S. Fairhurst, B. Krishnan, R. A. Mercer, R. K. Kopparapu,
36  * L. Santamaria, and J. T. Whelan,
37  * "Data formats for numerical relativity waves",
38  * arXiv:0709.0093v1 (2007).
39  *
40  * Currently only supports s=-2, l=2,3,4,5,6,7,8 modes.
41  */
42 COMPLEX16 XLALSpinWeightedSphericalHarmonic(
43  REAL8 theta, /**< polar angle (rad) */
44  REAL8 phi, /**< azimuthal angle (rad) */
45  int s, /**< spin weight */
46  int l, /**< mode number l */
47  int m /**< mode number m */
48  )
49 {
50  REAL8 fac;
51  COMPLEX16 ans;
52 
53  /* sanity checks ... */
54  if ( l < abs(s) )
55  {
56  XLALPrintError("XLAL Error - %s: Invalid mode s=%d, l=%d, m=%d - require |s| <= l\n", __func__, s, l, m );
57  XLAL_ERROR_VAL(0, XLAL_EINVAL);
58  }
59  if ( l < abs(m) )
60  {
61  XLALPrintError("XLAL Error - %s: Invalid mode s=%d, l=%d, m=%d - require |m| <= l\n", __func__, s, l, m );
62  XLAL_ERROR_VAL(0, XLAL_EINVAL);
63  }
64 
65  if ( s == -2 )
66  {
67  if ( l == 2 )
68  {
69  switch ( m )
70  {
71  case -2:
72  fac = sqrt( 5.0 / ( 64.0 * LAL_PI ) ) * ( 1.0 - cos( theta ))*( 1.0 - cos( theta ));
73  break;
74  case -1:
75  fac = sqrt( 5.0 / ( 16.0 * LAL_PI ) ) * sin( theta )*( 1.0 - cos( theta ));
76  break;
77 
78  case 0:
79  fac = sqrt( 15.0 / ( 32.0 * LAL_PI ) ) * sin( theta )*sin( theta );
80  break;
81 
82  case 1:
83  fac = sqrt( 5.0 / ( 16.0 * LAL_PI ) ) * sin( theta )*( 1.0 + cos( theta ));
84  break;
85 
86  case 2:
87  fac = sqrt( 5.0 / ( 64.0 * LAL_PI ) ) * ( 1.0 + cos( theta ))*( 1.0 + cos( theta ));
88  break;
89  default:
90  XLALPrintError("XLAL Error - %s: Invalid mode s=%d, l=%d, m=%d - require |m| <= l\n", __func__, s, l, m );
91  XLAL_ERROR_VAL(0, XLAL_EINVAL);
92  break;
93  } /* switch (m) */
94  } /* l==2*/
95  else if ( l == 3 )
96  {
97  switch ( m )
98  {
99  case -3:
100  fac = sqrt(21.0/(2.0*LAL_PI))*cos(theta/2.0)*pow(sin(theta/2.0),5.0);
101  break;
102  case -2:
103  fac = sqrt(7.0/(4.0*LAL_PI))*(2.0 + 3.0*cos(theta))*pow(sin(theta/2.0),4.0);
104  break;
105  case -1:
106  fac = sqrt(35.0/(2.0*LAL_PI))*(sin(theta) + 4.0*sin(2.0*theta) - 3.0*sin(3.0*theta))/32.0;
107  break;
108  case 0:
109  fac = (sqrt(105.0/(2.0*LAL_PI))*cos(theta)*pow(sin(theta),2.0))/4.0;
110  break;
111  case 1:
112  fac = -sqrt(35.0/(2.0*LAL_PI))*(sin(theta) - 4.0*sin(2.0*theta) - 3.0*sin(3.0*theta))/32.0;
113  break;
114 
115  case 2:
116  fac = sqrt(7.0/LAL_PI)*pow(cos(theta/2.0),4.0)*(-2.0 + 3.0*cos(theta))/2.0;
117  break;
118 
119  case 3:
120  fac = -sqrt(21.0/(2.0*LAL_PI))*pow(cos(theta/2.0),5.0)*sin(theta/2.0);
121  break;
122 
123  default:
124  XLALPrintError("XLAL Error - %s: Invalid mode s=%d, l=%d, m=%d - require |m| <= l\n", __func__, s, l, m );
125  XLAL_ERROR_VAL(0, XLAL_EINVAL);
126  break;
127  }
128  } /* l==3 */
129  else if ( l == 4 )
130  {
131  switch ( m )
132  {
133  case -4:
134  fac = 3.0*sqrt(7.0/LAL_PI)*pow(cos(theta/2.0),2.0)*pow(sin(theta/2.0),6.0);
135  break;
136  case -3:
137  fac = 3.0*sqrt(7.0/(2.0*LAL_PI))*cos(theta/2.0)*(1.0 + 2.0*cos(theta))*pow(sin(theta/2.0),5.0);
138  break;
139 
140  case -2:
141  fac = (3.0*(9.0 + 14.0*cos(theta) + 7.0*cos(2.0*theta))*pow(sin(theta/2.0),4.0))/(4.0*sqrt(LAL_PI));
142  break;
143  case -1:
144  fac = (3.0*(3.0*sin(theta) + 2.0*sin(2.0*theta) + 7.0*sin(3.0*theta) - 7.0*sin(4.0*theta)))/(32.0*sqrt(2.0*LAL_PI));
145  break;
146  case 0:
147  fac = (3.0*sqrt(5.0/(2.0*LAL_PI))*(5.0 + 7.0*cos(2.0*theta))*pow(sin(theta),2.0))/16.0;
148  break;
149  case 1:
150  fac = (3.0*(3.0*sin(theta) - 2.0*sin(2.0*theta) + 7.0*sin(3.0*theta) + 7.0*sin(4.0*theta)))/(32.0*sqrt(2.0*LAL_PI));
151  break;
152  case 2:
153  fac = (3.0*pow(cos(theta/2.0),4.0)*(9.0 - 14.0*cos(theta) + 7.0*cos(2.0*theta)))/(4.0*sqrt(LAL_PI));
154  break;
155  case 3:
156  fac = -3.0*sqrt(7.0/(2.0*LAL_PI))*pow(cos(theta/2.0),5.0)*(-1.0 + 2.0*cos(theta))*sin(theta/2.0);
157  break;
158  case 4:
159  fac = 3.0*sqrt(7.0/LAL_PI)*pow(cos(theta/2.0),6.0)*pow(sin(theta/2.0),2.0);
160  break;
161  default:
162  XLALPrintError("XLAL Error - %s: Invalid mode s=%d, l=%d, m=%d - require |m| <= l\n", __func__, s, l, m );
163  XLAL_ERROR_VAL(0, XLAL_EINVAL);
164  break;
165  }
166  } /* l==4 */
167  else if ( l == 5 )
168  {
169  switch ( m )
170  {
171  case -5:
172  fac = sqrt(330.0/LAL_PI)*pow(cos(theta/2.0),3.0)*pow(sin(theta/2.0),7.0);
173  break;
174  case -4:
175  fac = sqrt(33.0/LAL_PI)*pow(cos(theta/2.0),2.0)*(2.0 + 5.0*cos(theta))*pow(sin(theta/2.0),6.0);
176  break;
177  case -3:
178  fac = (sqrt(33.0/(2.0*LAL_PI))*cos(theta/2.0)*(17.0 + 24.0*cos(theta) + 15.0*cos(2.0*theta))*pow(sin(theta/2.0),5.0))/4.0;
179  break;
180  case -2:
181  fac = (sqrt(11.0/LAL_PI)*(32.0 + 57.0*cos(theta) + 36.0*cos(2.0*theta) + 15.0*cos(3.0*theta))*pow(sin(theta/2.0),4.0))/8.0;
182  break;
183  case -1:
184  fac = (sqrt(77.0/LAL_PI)*(2.0*sin(theta) + 8.0*sin(2.0*theta) + 3.0*sin(3.0*theta) + 12.0*sin(4.0*theta) - 15.0*sin(5.0*theta)))/256.0;
185  break;
186  case 0:
187  fac = (sqrt(1155.0/(2.0*LAL_PI))*(5.0*cos(theta) + 3.0*cos(3.0*theta))*pow(sin(theta),2.0))/32.0;
188  break;
189  case 1:
190  fac = sqrt(77.0/LAL_PI)*(-2.0*sin(theta) + 8.0*sin(2.0*theta) - 3.0*sin(3.0*theta) + 12.0*sin(4.0*theta) + 15.0*sin(5.0*theta))/256.0;
191  break;
192  case 2:
193  fac = sqrt(11.0/LAL_PI)*pow(cos(theta/2.0),4.0)*(-32.0 + 57.0*cos(theta) - 36.0*cos(2.0*theta) + 15.0*cos(3.0*theta))/8.0;
194  break;
195  case 3:
196  fac = -sqrt(33.0/(2.0*LAL_PI))*pow(cos(theta/2.0),5.0)*(17.0 - 24.0*cos(theta) + 15.0*cos(2.0*theta))*sin(theta/2.0)/4.0;
197  break;
198  case 4:
199  fac = sqrt(33.0/LAL_PI)*pow(cos(theta/2.0),6.0)*(-2.0 + 5.0*cos(theta))*pow(sin(theta/2.0),2.0);
200  break;
201  case 5:
202  fac = -sqrt(330.0/LAL_PI)*pow(cos(theta/2.0),7.0)*pow(sin(theta/2.0),3.0);
203  break;
204  default:
205  XLALPrintError("XLAL Error - %s: Invalid mode s=%d, l=%d, m=%d - require |m| <= l\n", __func__, s, l, m );
206  XLAL_ERROR_VAL(0, XLAL_EINVAL);
207  break;
208  }
209  } /* l==5 */
210  else if ( l == 6 )
211  {
212  switch ( m )
213  {
214  case -6:
215  fac = (3.*sqrt(715./LAL_PI)*pow(cos(theta/2.0),4)*pow(sin(theta/2.0),8))/2.0;
216  break;
217  case -5:
218  fac = (sqrt(2145./LAL_PI)*pow(cos(theta/2.0),3)*(1. + 3.*cos(theta))*pow(sin(theta/2.0),7))/2.0;
219  break;
220  case -4:
221  fac = (sqrt(195./(2.0*LAL_PI))*pow(cos(theta/2.0),2)*(35. + 44.*cos(theta)
222  + 33.*cos(2.*theta))*pow(sin(theta/2.0),6))/8.0;
223  break;
224  case -3:
225  fac = (3.*sqrt(13./LAL_PI)*cos(theta/2.0)*(98. + 185.*cos(theta) + 110.*cos(2*theta)
226  + 55.*cos(3.*theta))*pow(sin(theta/2.0),5))/32.0;
227  break;
228  case -2:
229  fac = (sqrt(13./LAL_PI)*(1709. + 3096.*cos(theta) + 2340.*cos(2.*theta) + 1320.*cos(3.*theta)
230  + 495.*cos(4.*theta))*pow(sin(theta/2.0),4))/256.0;
231  break;
232  case -1:
233  fac = (sqrt(65./(2.0*LAL_PI))*cos(theta/2.0)*(161. + 252.*cos(theta) + 252.*cos(2.*theta)
234  + 132.*cos(3.*theta) + 99.*cos(4.*theta))*pow(sin(theta/2.0),3))/64.0;
235  break;
236  case 0:
237  fac = (sqrt(1365./LAL_PI)*(35. + 60.*cos(2.*theta) + 33.*cos(4.*theta))*pow(sin(theta),2))/512.0;
238  break;
239  case 1:
240  fac = (sqrt(65./(2.0*LAL_PI))*pow(cos(theta/2.0),3)*(161. - 252.*cos(theta) + 252.*cos(2.*theta)
241  - 132.*cos(3.*theta) + 99.*cos(4.*theta))*sin(theta/2.0))/64.0;
242  break;
243  case 2:
244  fac = (sqrt(13./LAL_PI)*pow(cos(theta/2.0),4)*(1709. - 3096.*cos(theta) + 2340.*cos(2.*theta)
245  - 1320*cos(3*theta) + 495*cos(4*theta)))/256.0;
246  break;
247  case 3:
248  fac = (-3.*sqrt(13./LAL_PI)*pow(cos(theta/2.0),5)*(-98. + 185.*cos(theta) - 110.*cos(2*theta)
249  + 55.*cos(3.*theta))*sin(theta/2.0))/32.0;
250  break;
251  case 4:
252  fac = (sqrt(195./(2.0*LAL_PI))*pow(cos(theta/2.0),6)*(35. - 44.*cos(theta)
253  + 33.*cos(2*theta))*pow(sin(theta/2.0),2))/8.0;
254  break;
255  case 5:
256  fac = -(sqrt(2145./LAL_PI)*pow(cos(theta/2.0),7)*(-1. + 3.*cos(theta))*pow(sin(theta/2.0),3))/2.0;
257  break;
258  case 6:
259  fac = (3.*sqrt(715./LAL_PI)*pow(cos(theta/2.0),8)*pow(sin(theta/2.0),4))/2.0;
260  break;
261  default:
262  XLALPrintError("XLAL Error - %s: Invalid mode s=%d, l=%d, m=%d - require |m| <= l\n", __func__, s, l, m );
263  XLAL_ERROR_VAL(0, XLAL_EINVAL);
264  break;
265  }
266  } /* l==6 */
267  else if ( l == 7 )
268  {
269  switch ( m )
270  {
271  case -7:
272  fac = sqrt(15015./(2.0*LAL_PI))*pow(cos(theta/2.0),5)*pow(sin(theta/2.0),9);
273  break;
274  case -6:
275  fac = (sqrt(2145./LAL_PI)*pow(cos(theta/2.0),4)*(2. + 7.*cos(theta))*pow(sin(theta/2.0),8))/2.0;
276  break;
277  case -5:
278  fac = (sqrt(165./(2.0*LAL_PI))*pow(cos(theta/2.0),3)*(93. + 104.*cos(theta)
279  + 91.*cos(2.*theta))*pow(sin(theta/2.0),7))/8.0;
280  break;
281  case -4:
282  fac = (sqrt(165./(2.0*LAL_PI))*pow(cos(theta/2.0),2)*(140. + 285.*cos(theta)
283  + 156.*cos(2.*theta) + 91.*cos(3.*theta))*pow(sin(theta/2.0),6))/16.0;
284  break;
285  case -3:
286  fac = (sqrt(15./(2.0*LAL_PI))*cos(theta/2.0)*(3115. + 5456.*cos(theta) + 4268.*cos(2.*theta)
287  + 2288.*cos(3.*theta) + 1001.*cos(4.*theta))*pow(sin(theta/2.0),5))/128.0;
288  break;
289  case -2:
290  fac = (sqrt(15./LAL_PI)*(5220. + 9810.*cos(theta) + 7920.*cos(2.*theta) + 5445.*cos(3.*theta)
291  + 2860.*cos(4.*theta) + 1001.*cos(5.*theta))*pow(sin(theta/2.0),4))/512.0;
292  break;
293  case -1:
294  fac = (3.*sqrt(5./(2.0*LAL_PI))*cos(theta/2.0)*(1890. + 4130.*cos(theta) + 3080.*cos(2.*theta)
295  + 2805.*cos(3.*theta) + 1430.*cos(4.*theta) + 1001.*cos(5*theta))*pow(sin(theta/2.0),3))/512.0;
296  break;
297  case 0:
298  fac = (3.*sqrt(35./LAL_PI)*cos(theta)*(109. + 132.*cos(2.*theta)
299  + 143.*cos(4.*theta))*pow(sin(theta),2))/512.0;
300  break;
301  case 1:
302  fac = (3.*sqrt(5./(2.0*LAL_PI))*pow(cos(theta/2.0),3)*(-1890. + 4130.*cos(theta) - 3080.*cos(2.*theta)
303  + 2805.*cos(3.*theta) - 1430.*cos(4.*theta) + 1001.*cos(5.*theta))*sin(theta/2.0))/512.0;
304  break;
305  case 2:
306  fac = (sqrt(15./LAL_PI)*pow(cos(theta/2.0),4)*(-5220. + 9810.*cos(theta) - 7920.*cos(2.*theta)
307  + 5445.*cos(3.*theta) - 2860.*cos(4.*theta) + 1001.*cos(5.*theta)))/512.0;
308  break;
309  case 3:
310  fac = -(sqrt(15./(2.0*LAL_PI))*pow(cos(theta/2.0),5)*(3115. - 5456.*cos(theta) + 4268.*cos(2.*theta)
311  - 2288.*cos(3.*theta) + 1001.*cos(4.*theta))*sin(theta/2.0))/128.0;
312  break;
313  case 4:
314  fac = (sqrt(165./(2.0*LAL_PI))*pow(cos(theta/2.0),6)*(-140. + 285.*cos(theta) - 156.*cos(2*theta)
315  + 91.*cos(3.*theta))*pow(sin(theta/2.0),2))/16.0;
316  break;
317  case 5:
318  fac = -(sqrt(165./(2.0*LAL_PI))*pow(cos(theta/2.0),7)*(93. - 104.*cos(theta)
319  + 91.*cos(2.*theta))*pow(sin(theta/2.0),3))/8.0;
320  break;
321  case 6:
322  fac = (sqrt(2145./LAL_PI)*pow(cos(theta/2.0),8)*(-2. + 7.*cos(theta))*pow(sin(theta/2.0),4))/2.0;
323  break;
324  case 7:
325  fac = -(sqrt(15015./(2.0*LAL_PI))*pow(cos(theta/2.0),9)*pow(sin(theta/2.0),5));
326  break;
327  default:
328  XLALPrintError("XLAL Error - %s: Invalid mode s=%d, l=%d, m=%d - require |m| <= l\n", __func__, s, l, m );
329  XLAL_ERROR_VAL(0, XLAL_EINVAL);
330  break;
331  }
332  } /* l==7 */
333  else if ( l == 8 )
334  {
335  switch ( m )
336  {
337  case -8:
338  fac = sqrt(34034./LAL_PI)*pow(cos(theta/2.0),6)*pow(sin(theta/2.0),10);
339  break;
340  case -7:
341  fac = sqrt(17017./(2.0*LAL_PI))*pow(cos(theta/2.0),5)*(1. + 4.*cos(theta))*pow(sin(theta/2.0),9);
342  break;
343  case -6:
344  fac = sqrt(255255./LAL_PI)*pow(cos(theta/2.0),4)*(1. + 2.*cos(theta))
345  *sin(LAL_PI/4.0 - theta/2.0)*sin(LAL_PI/4.0 + theta/2.0)*pow(sin(theta/2.0),8);
346  break;
347  case -5:
348  fac = (sqrt(12155./(2.0*LAL_PI))*pow(cos(theta/2.0),3)*(19. + 42.*cos(theta)
349  + 21.*cos(2.*theta) + 14.*cos(3.*theta))*pow(sin(theta/2.0),7))/8.0;
350  break;
351  case -4:
352  fac = (sqrt(935./(2.0*LAL_PI))*pow(cos(theta/2.0),2)*(265. + 442.*cos(theta) + 364.*cos(2.*theta)
353  + 182.*cos(3.*theta) + 91.*cos(4.*theta))*pow(sin(theta/2.0),6))/32.0;
354  break;
355  case -3:
356  fac = (sqrt(561./(2.0*LAL_PI))*cos(theta/2.0)*(869. + 1660.*cos(theta) + 1300.*cos(2.*theta)
357  + 910.*cos(3.*theta) + 455.*cos(4.*theta) + 182.*cos(5.*theta))*pow(sin(theta/2.0),5))/128.0;
358  break;
359  case -2:
360  fac = (sqrt(17./LAL_PI)*(7626. + 14454.*cos(theta) + 12375.*cos(2.*theta) + 9295.*cos(3.*theta)
361  + 6006.*cos(4.*theta) + 3003.*cos(5.*theta) + 1001.*cos(6.*theta))*pow(sin(theta/2.0),4))/512.0;
362  break;
363  case -1:
364  fac = (sqrt(595./(2.0*LAL_PI))*cos(theta/2.0)*(798. + 1386.*cos(theta) + 1386.*cos(2.*theta)
365  + 1001.*cos(3.*theta) + 858.*cos(4.*theta) + 429.*cos(5.*theta) + 286.*cos(6.*theta))*pow(sin(theta/2.0),3))/512.0;
366  break;
367  case 0:
368  fac = (3.*sqrt(595./LAL_PI)*(210. + 385.*cos(2.*theta) + 286.*cos(4.*theta)
369  + 143.*cos(6.*theta))*pow(sin(theta),2))/4096.0;
370  break;
371  case 1:
372  fac = (sqrt(595./(2.0*LAL_PI))*pow(cos(theta/2.0),3)*(798. - 1386.*cos(theta) + 1386.*cos(2.*theta)
373  - 1001.*cos(3.*theta) + 858.*cos(4.*theta) - 429.*cos(5.*theta) + 286.*cos(6.*theta))*sin(theta/2.0))/512.0;
374  break;
375  case 2:
376  fac = (sqrt(17./LAL_PI)*pow(cos(theta/2.0),4)*(7626. - 14454.*cos(theta) + 12375.*cos(2.*theta)
377  - 9295.*cos(3.*theta) + 6006.*cos(4.*theta) - 3003.*cos(5.*theta) + 1001.*cos(6.*theta)))/512.0;
378  break;
379  case 3:
380  fac = -(sqrt(561./(2.0*LAL_PI))*pow(cos(theta/2.0),5)*(-869. + 1660.*cos(theta) - 1300.*cos(2.*theta)
381  + 910.*cos(3.*theta) - 455.*cos(4.*theta) + 182.*cos(5.*theta))*sin(theta/2.0))/128.0;
382  break;
383  case 4:
384  fac = (sqrt(935./(2.0*LAL_PI))*pow(cos(theta/2.0),6)*(265. - 442.*cos(theta) + 364.*cos(2.*theta)
385  - 182.*cos(3.*theta) + 91.*cos(4.*theta))*pow(sin(theta/2.0),2))/32.0;
386  break;
387  case 5:
388  fac = -(sqrt(12155./(2.0*LAL_PI))*pow(cos(theta/2.0),7)*(-19. + 42.*cos(theta) - 21.*cos(2.*theta)
389  + 14.*cos(3.*theta))*pow(sin(theta/2.0),3))/8.0;
390  break;
391  case 6:
392  fac = sqrt(255255./LAL_PI)*pow(cos(theta/2.0),8)*(-1. + 2.*cos(theta))*sin(LAL_PI/4.0 - theta/2.0)
393  *sin(LAL_PI/4.0 + theta/2.0)*pow(sin(theta/2.0),4);
394  break;
395  case 7:
396  fac = -(sqrt(17017./(2.0*LAL_PI))*pow(cos(theta/2.0),9)*(-1. + 4.*cos(theta))*pow(sin(theta/2.0),5));
397  break;
398  case 8:
399  fac = sqrt(34034./LAL_PI)*pow(cos(theta/2.0),10)*pow(sin(theta/2.0),6);
400  break;
401  default:
402  XLALPrintError("XLAL Error - %s: Invalid mode s=%d, l=%d, m=%d - require |m| <= l\n", __func__, s, l, m );
403  XLAL_ERROR_VAL(0, XLAL_EINVAL);
404  break;
405  }
406  } /* l==8 */
407  else
408  {
409  XLALPrintError("XLAL Error - %s: Unsupported mode l=%d (only l in [2,8] implemented)\n", __func__, l);
410  XLAL_ERROR_VAL(0, XLAL_EINVAL);
411  }
412  }
413  else
414  {
415  XLALPrintError("XLAL Error - %s: Unsupported mode s=%d (only s=-2 implemented)\n", __func__, s);
416  XLAL_ERROR_VAL(0, XLAL_EINVAL);
417  }
418  if (m)
419  ans = cpolar(1.0, m*phi) * fac;
420  else
421  ans = fac;
422  return ans;
423 }
424 
425 
426 /**
427  * Computes the scalar spherical harmonic \f$ Y_{lm}(\theta, \phi) \f$.
428  */
429 int
430 XLALScalarSphericalHarmonic(
431  COMPLEX16 *y, /**< output */
432  UINT4 l, /**< value of l */
433  INT4 m, /**< value of m */
434  REAL8 theta, /**< angle theta */
435  REAL8 phi /**< angle phi */
436  )
437 {
438 
439  int gslStatus;
440  gsl_sf_result pLm;
441 
442  INT4 absM = abs( m );
443 
444  if ( absM > (INT4) l )
445  {
446  XLAL_ERROR( XLAL_EINVAL );
447  }
448 
449  /* For some reason GSL will not take negative m */
450  /* We will have to use the relation between sph harmonics of +ve and -ve m */
451  XLAL_CALLGSL( gslStatus = gsl_sf_legendre_sphPlm_e((INT4)l, absM, cos(theta), &pLm ) );
452  if (gslStatus != GSL_SUCCESS)
453  {
454  XLALPrintError("Error in GSL function\n" );
455  XLAL_ERROR( XLAL_EFUNC );
456  }
457 
458  /* Compute the values for the spherical harmonic */
459  *y = cpolar(pLm.val, m * phi);
460 
461  /* If m is negative, perform some jiggery-pokery */
462  if ( m < 0 && absM % 2 == 1 )
463  {
464  *y = - *y;
465  }
466 
467  return XLAL_SUCCESS;
468 }
469 
470 /**
471  * Computes the spin 2 weighted spherical harmonic. This function is now
472  * deprecated and will be removed soon. All calls should be replaced with
473  * calls to XLALSpinWeightedSphericalHarmonic().
474  */
475 INT4 XLALSphHarm ( COMPLEX16 *out, /**< output */
476  UINT4 L, /**< value of L */
477  INT4 M, /**< value of M */
478  REAL4 theta, /**< angle with respect to the z axis */
479  REAL4 phi /**< angle with respect to the x axis */
480  )
481 {
482 
483  XLAL_PRINT_DEPRECATION_WARNING("XLALSpinWeightedSphericalHarmonic");
484 
485  *out = XLALSpinWeightedSphericalHarmonic( theta, phi, -2, L, M );
486  if ( xlalErrno )
487  {
488  XLAL_ERROR( XLAL_EFUNC );
489  }
490 
491  return XLAL_SUCCESS;
492 }
493 
494 /**
495  * Computes the n-th Jacobi polynomial for polynomial weights alpha and beta.
496  * The implementation here is only valid for real x -- enforced by the argument
497  * type. An extension to complex values would require evaluation of several
498  * gamma functions.
499  *
500  * See http://en.wikipedia.org/wiki/Jacobi_polynomials
501  */
502 double XLALJacobiPolynomial(int n, int alpha, int beta, double x){
503  double f1 = (x-1)/2.0, f2 = (x+1)/2.0;
504  int s=0;
505  double sum=0, val=0;
506  if( n == 0 ) return 1.0;
507  for( s=0; n-s >= 0; s++ ){
508  val=1.0;
509  val *= gsl_sf_choose( n+alpha, s );
510  val *= gsl_sf_choose( n+beta, n-s );
511  if( n-s != 0 ) val *= pow( f1, n-s );
512  if( s != 0 ) val*= pow( f2, s );
513 
514  sum += val;
515  }
516  return sum;
517 }
518 
519 /**
520  * Computes the 'little' d Wigner matrix for the Euler angle beta. Single angle
521  * small d transform with major index 'l' and minor index transition from m to
522  * mp.
523  *
524  * Uses a slightly unconventional method since the intuitive version by Wigner
525  * is less suitable to algorthmic development.
526  *
527  * See http://en.wikipedia.org/wiki/Wigner_D-matrix#Wigner_.28small.29_d-matrix
528  */
529 #define MIN(a,b) ((a) < (b) ? (a) : (b))
530 double XLALWignerdMatrix(
531  int l, /**< mode number l */
532  int mp, /**< mode number m' */
533  int m, /**< mode number m */
534  double beta /**< euler angle (rad) */
535  )
536 {
537 
538  int k = MIN( l+m, MIN( l-m, MIN( l+mp, l-mp )));
539  double a=0, lam=0;
540  if(k == l+m){
541  a = mp-m;
542  lam = mp-m;
543  } else if(k == l-m) {
544  a = m-mp;
545  lam = 0;
546  } else if(k == l+mp) {
547  a = m-mp;
548  lam = 0;
549  } else if(k == l-mp) {
550  a = mp-m;
551  lam = mp-m;
552  }
553 
554  int b = 2*l-2*k-a;
555  double pref = pow(-1, lam) * sqrt(gsl_sf_choose( 2*l-k, k+a )) / sqrt(gsl_sf_choose( k+b, b ));
556 
557  return pref * pow(sin(beta/2.0), a) * pow( cos(beta/2.0), b) * XLALJacobiPolynomial(k, a, b, cos(beta));
558 
559 }
560 
561 /**
562  * Computes the full Wigner D matrix for the Euler angle alpha, beta, and gamma
563  * with major index 'l' and minor index transition from m to mp.
564  *
565  * Uses a slightly unconventional method since the intuitive version by Wigner
566  * is less suitable to algorthmic development.
567  *
568  * See http://en.wikipedia.org/wiki/Wigner_D-matrix
569  *
570  * Currently only supports the modes which are implemented for the spin
571  * weighted spherical harmonics.
572  */
573 COMPLEX16 XLALWignerDMatrix(
574  int l, /**< mode number l */
575  int mp, /**< mode number m' */
576  int m, /**< mode number m */
577  double alpha, /**< euler angle (rad) */
578  double beta, /**< euler angle (rad) */
579  double gam /**< euler angle (rad) */
580  )
581 {
582  return cexp( -(1.0I)*mp*alpha ) *
583  XLALWignerdMatrix( l, mp, m, beta ) *
584  cexp( -(1.0I)*m*gam );
585 }
MIN
#define MIN(a, b)
Computes the 'little' d Wigner matrix for the Euler angle beta.
Definition: SphericalHarmonics.c:529
XLALPrintError
int XLALPrintError(const char *fmt,...)
Definition: XLALError.c:68
XLAL_CALLGSL
#define XLAL_CALLGSL(statement)
Definition: XLALGSL.h:33
LAL_PI
#define LAL_PI
Archimedes's constant, pi.
Definition: LALConstants.h:179
cpolar
#define cpolar(r, th)
Construct a COMPLEX16 from polar modulus and argument.
Definition: LALAtomicDatatypes.h:125
COMPLEX16
double complex COMPLEX16
Double-precision floating-point complex number (16 bytes total)
Definition: LALAtomicDatatypes.h:120
REAL8
double REAL8
Double precision real floating-point number (8 bytes).
Definition: LALAtomicDatatypes.h:103
UINT4
uint32_t UINT4
Four-byte unsigned integer.
Definition: LALAtomicDatatypes.h:59
INT4
int32_t INT4
Four-byte signed integer.
Definition: LALAtomicDatatypes.h:57
REAL4
float REAL4
Single precision real floating-point number (4 bytes).
Definition: LALAtomicDatatypes.h:102
m
static const INT4 m
Definition: Random.c:80
a
static const INT4 a
Definition: Random.c:79
XLALWignerdMatrix
double XLALWignerdMatrix(int l, int mp, int m, double beta)
Definition: SphericalHarmonics.c:530
XLALWignerDMatrix
COMPLEX16 XLALWignerDMatrix(int l, int mp, int m, double alpha, double beta, double gam)
Computes the full Wigner D matrix for the Euler angle alpha, beta, and gamma with major index 'l' and...
Definition: SphericalHarmonics.c:573
XLALJacobiPolynomial
double XLALJacobiPolynomial(int n, int alpha, int beta, double x)
Computes the n-th Jacobi polynomial for polynomial weights alpha and beta.
Definition: SphericalHarmonics.c:502
XLALSpinWeightedSphericalHarmonic
COMPLEX16 XLALSpinWeightedSphericalHarmonic(REAL8 theta, REAL8 phi, int s, int l, int m)
Computes the (s)Y(l,m) spin-weighted spherical harmonic.
Definition: SphericalHarmonics.c:42
XLALScalarSphericalHarmonic
int XLALScalarSphericalHarmonic(COMPLEX16 *y, UINT4 l, INT4 m, REAL8 theta, REAL8 phi)
Computes the scalar spherical harmonic .
Definition: SphericalHarmonics.c:430
XLALSphHarm
INT4 XLALSphHarm(COMPLEX16 *out, UINT4 L, INT4 M, REAL4 theta, REAL4 phi)
Computes the spin 2 weighted spherical harmonic.
Definition: SphericalHarmonics.c:475
XLAL_ERROR_VAL
#define XLAL_ERROR_VAL(val,...)
Macro to invoke the XLALError() function and return with code val (it should not really be used itsel...
Definition: XLALError.h:687
xlalErrno
#define xlalErrno
Modifiable lvalue containing the XLAL error number.
Definition: XLALError.h:571
XLAL_ERROR
#define XLAL_ERROR(...)
Macro to invoke a failure from a XLAL routine returning an integer.
Definition: XLALError.h:700
XLAL_PRINT_DEPRECATION_WARNING
#define XLAL_PRINT_DEPRECATION_WARNING(replacement)
Prints a deprecation warning at the "warning" verbosity level.
Definition: XLALError.h:228
XLAL_SUCCESS
@ XLAL_SUCCESS
Success return value (not an error number)
Definition: XLALError.h:401
XLAL_EFUNC
@ XLAL_EFUNC
Internal function call failed bit: "or" this with existing error number.
Definition: XLALError.h:462
XLAL_EINVAL
@ XLAL_EINVAL
Invalid argument.
Definition: XLALError.h:409
test_iterutils.x
x
Definition: test_iterutils.py:51