Loading [MathJax]/extensions/TeX/AMSsymbols.js
LAL 7.7.0.1-3a66518
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
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 */
42COMPLEX16 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 */
429int
430XLALScalarSphericalHarmonic(
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 */
475INT4 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 */
502double 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))
530double 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 */
573COMPLEX16 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:50