Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-5e288d3
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALSimNRSurrogateUtilities.c
Go to the documentation of this file.
1/* Copyright (C) 2018 Jonathan Blackman
2 * Utility functions that are useful for evaluating NR surrogates.
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
21#include <lal/LALConstants.h>
22
23
24/****************** Init and Destroy functions ****************/
25
26/**
27 * Initialize a MultiModalWaveforme. We also have the usual frame data, which is initially the
28m structure
29 */
31 MultiModalWaveform **wave, /**< Output; the MultiModalWaveform. */
32 int LMax, /**< (ell, m) modes with 2 <= ell <= LMax will be created. */
33 int n_times /**< Number of time samples per mode. */
34) {
35 if (!wave) exit(1);
36 if (LMax < 2) XLAL_ERROR_VOID(XLAL_FAILURE, "Got LMax=%d < 2!\n", LMax);
37 if (*wave) MultiModalWaveform_Destroy(*wave);
38 (*wave) = XLALCalloc(1, sizeof(MultiModalWaveform));
39
40 int n_modes = LMax*(LMax+2) - 3;
41
42 (*wave)->n_modes = n_modes;
43 (*wave)->lvals = XLALCalloc(n_modes, sizeof(int));
44 (*wave)->mvals = XLALCalloc(n_modes, sizeof(int));
45 (*wave)->n_times = n_times;
46 gsl_vector **modes_real_part = XLALMalloc(n_modes * sizeof(*modes_real_part));
47 gsl_vector **modes_imag_part = XLALMalloc(n_modes * sizeof(*modes_imag_part));
48 (*wave)->modes_real_part = modes_real_part;
49 (*wave)->modes_imag_part = modes_imag_part;
50
51 int ell=2;
52 int m=-2;
53 for (int i=0; i<n_modes; i++) {
54 (*wave)->lvals[i] = ell;
55 (*wave)->mvals[i] = m;
56 (*wave)->modes_real_part[i] = gsl_vector_calloc(n_times);
57 (*wave)->modes_imag_part[i] = gsl_vector_calloc(n_times);
58 m += 1;
59 if (m > ell) {
60 ell += 1;
61 m = -1 * ell;
62 }
63 }
64}
65
66/**
67 * Destroy a MultiModalWaveform structure; free all associated memory.
68 */
70 if (!wave) return;
71 for (int i=0; i<wave->n_modes; i++) {
72 if (wave->modes_real_part[i]) gsl_vector_free(wave->modes_real_part[i]);
73 if (wave->modes_imag_part[i]) gsl_vector_free(wave->modes_imag_part[i]);
74 }
75 free(wave->modes_real_part);
76 free(wave->modes_imag_part);
77 XLALFree(wave->lvals);
78 XLALFree(wave->mvals);
79 XLALFree(wave);
80}
81
82/**
83 * Initialize a ComplexPowers structure
84 */
86 ComplexPowers **cp, /**< The ComplexPowers structure to be initialized.*/
87 int LMax, /**< The maximum ell that will be needed by WignerDMatrices.*/
88 int n_times /**< The number of time samples.*/
89) {
90 // include powers from -2*LMax to 2*LMax inclusive
91 if (!cp) exit(1);
92 if (*cp) ComplexPowers_Destroy(*cp);
93 (*cp) = XLALCalloc(1, sizeof(ComplexPowers));
94
95 int n_entries = 4*LMax+1;
96 (*cp)->LMax = LMax;
97 (*cp)->n_entries = n_entries;
98
99 gsl_vector **real_part = XLALMalloc(n_entries * sizeof(*real_part));
100 gsl_vector **imag_part = XLALMalloc(n_entries * sizeof(*imag_part));
101 (*cp)->real_part = real_part;
102 (*cp)->imag_part = imag_part;
103
104 for (int i=0; i<n_entries; i++) {
105 (*cp)->real_part[i] = gsl_vector_calloc(n_times);
106 (*cp)->imag_part[i] = gsl_vector_calloc(n_times);
107 }
108}
109
110/**
111 * Destroy a ComplexPowers structure; free all associated memory.
112 */
114 if (!cp) return;
115 for (int i=0; i<cp->n_entries; i++) {
116 if (cp->real_part[i]) gsl_vector_free(cp->real_part[i]);
117 if (cp->imag_part[i]) gsl_vector_free(cp->imag_part[i]);
118 }
119 free(cp->real_part);
120 free(cp->imag_part);
121 XLALFree(cp);
122}
123
124/**
125 * Initialize a RealPowers structure
126 */
127static void RealPowers_Init(
128 RealPowers **rp, /**< The RealPowers structure to be initialized.*/
129 int LMax, /**< The maximum ell that will be needed by WignerDMatrices.*/
130 int n_times /**< The number of time samples.*/
131) {
132 // include powers from 0 to 2*LMax inclusive
133 if (!rp) exit(1);
134 if (*rp) RealPowers_Destroy(*rp);
135 (*rp) = XLALCalloc(1, sizeof(RealPowers));
136
137 int n_entries = 2*LMax+1;
138 (*rp)->LMax = LMax;
139 (*rp)->n_entries = n_entries;
140
141 gsl_vector **powers = XLALMalloc(n_entries * sizeof(*powers));
142 (*rp)->powers = powers;
143
144 for (int i=0; i<n_entries; i++) {
145 (*rp)->powers[i] = gsl_vector_calloc(n_times);
146 }
147}
148
149/**
150 * Destroy a RealPowers structure; free all associated memory.
151 */
153 if (!rp) return;
154 for (int i=0; i<rp->n_entries; i++) {
155 if (rp->powers[i]) gsl_vector_free(rp->powers[i]);
156 }
157 free(rp->powers);
158 XLALFree(rp);
159}
160
161/**
162 * Initialize a WignerDMatrices structure
163 */
165 WignerDMatrices **matrices, /**< Output; the WignerDMatrices.*/
166 int n_times, /**< The number of time samples for each (ell, m, m') component.*/
167 int LMax /**< The maximum ell to generate. */
168) {
169 if (!matrices) exit(1);
170 if (LMax < 2) XLAL_ERROR_VOID(XLAL_FAILURE, "Got LMax=%d < 2!\n", LMax);
171 if (*matrices) WignerDMatrices_Destroy(*matrices);
172 (*matrices) = XLALCalloc(1, sizeof(WignerDMatrices));
173
174 int n_entries = 0;
175 for (int ell=2; ell<=LMax; ell++) {
176 n_entries += (2*ell+1) * (2*ell+1);
177 }
178
179 (*matrices)->LMax = LMax;
180 (*matrices)->n_entries = n_entries;
181 (*matrices)->n_times = n_times;
182
183 gsl_vector **real_part = XLALMalloc(n_entries * sizeof(*real_part));
184 gsl_vector **imag_part = XLALMalloc(n_entries * sizeof(*imag_part));
185 (*matrices)->real_part = real_part;
186 (*matrices)->imag_part = imag_part;
187
188 for (int i=0; i<n_entries; i++) {
189 (*matrices)->real_part[i] = gsl_vector_calloc(n_times);
190 (*matrices)->imag_part[i] = gsl_vector_calloc(n_times);
191 }
192}
193
194/**
195 * Destroy a WignerDMatrices structure; free all associated memory.
196 */
198 if (!matrices) return;
199 for (int i=0; i<matrices->n_entries; i++) {
200 if (matrices->real_part[i]) gsl_vector_free(matrices->real_part[i]);
201 if (matrices->imag_part[i]) gsl_vector_free(matrices->imag_part[i]);
202 }
203 free(matrices->real_part);
204 free(matrices->imag_part);
205 XLALFree(matrices);
206}
207
208
209
210/**************** Utility functions **************************/
211
212/**
213 * Computes the index of the (ell, m, mp) component in a WignerDMatices structure
214 */
215static int WignerDMatrix_Index(int ell, int m, int mp) {
216 // Start of the (m, mp) matrix, which has size (2*ell + 1) X (2*ell + 1)
217 int i0 = (ell*(ell*ell*4 - 1))/3 - 10;
218 int res = i0 + (2*ell+1)*(ell+m) + (ell+mp);
219 return res;
220}
221
222static REAL8 factorial(int n) {
223 if (n <= 0) return 1.0;
224 return factorial(n-1) * n;
225}
226
227static REAL8 factorial_ratio(int n, int k) {
228 if (n <= k) return 1.0;
229 return factorial_ratio(n-1, k) * n;
230}
231
232static REAL8 binomial(int n, int k) {
233 return factorial_ratio(n, k) / factorial(n-k);
234}
235
236static REAL8 wigner_coef(int ell, int mp, int m) {
237 return sqrt(factorial(ell+m) * factorial(ell-m) / (factorial(ell+mp) * factorial(ell-mp)));
238}
239
240/**
241 * Multiplies z1(t) * z2(t) = (x1(t) + i*y1(t)) * (x2(t) + i*y2(t)),
242 * storing the result in x1 and y1.
243 */
245 gsl_vector *x1, /**< Real part of z1. */
246 gsl_vector *y1, /**< Imag part of z1. */
247 gsl_vector *x2, /**< Real part of z2. */
248 gsl_vector *y2, /**< Imag part of z2. */
249 gsl_vector *tmpx, /**< A vector for storing temporary results. */
250 gsl_vector *tmpy /**< A second vector for storing temporary results. */
251) {
252 gsl_vector_set_zero(tmpx);
253 gsl_vector_set_zero(tmpy);
254
255 // Store the temporary results x1*y2 and y1*y2
256 gsl_vector_add(tmpx, y1);
257 gsl_vector_mul(tmpx, y2);
258 gsl_vector_add(tmpy, x1);
259 gsl_vector_mul(tmpy, y2);
260
261 // Now we can safely transform x1->x1*x2 and y1->y1*x2
262 gsl_vector_mul(x1, x2);
263 gsl_vector_mul(y1, x2);
264
265 // Now we add in the temporary results
266 gsl_vector_sub(x1, tmpx);
267 gsl_vector_add(y1, tmpy);
268}
269
270
271/*************** Heavy lifting functions ********************/
272/*
273 * Given dy/dt evaluated as a function of y and t at 4 time nodes, compute
274 * dy = y_5 - y_4 using a 4th-order accurate Adams-Bashforth scheme.
275 * See https://dcc.ligo.org/T1800069 for details.
276 */
277static void NRSur_ab4_dy(
278 REAL8 *dy, /**< Result */
279 REAL8 *k1, /**< dy/dt evaluated at node 1 */
280 REAL8 *k2, /**< dy/dt evaluated at node 2 */
281 REAL8 *k3, /**< dy/dt evaluated at node 3 */
282 REAL8 *k4, /**< dy/dt evaluated at node 4 */
283 REAL8 t12, /**< t_2 - t_1 */
284 REAL8 t23, /**< t_3 - t_2 */
285 REAL8 t34, /**< t_4 - t_3 */
286 REAL8 t45, /**< t_5 - t_4 */
287 int dim /**< Vector dimension. Length of dy, k1, k2, k3, k4. */
288) {
289
290 // These are defined in https://dcc.ligo.org/T1800069, except we rename Dj to denomj.
291 REAL8 t13, t14, t24, denom1, denom2, denom3;
292 REAL8 A1, A2, A3, A4, B1, B2, B3, B4, C1, C2, C3, C4; // Matrix coefficients in the M matrix.
293 REAL8 dx_vec[4]; // Integral of the p vector from t4 to t5; [t45, t45^2/2, t45^3/3, t45^4/4]
294
295 // Temp variable for the coefficients multiplying 1, x, x^2, and x^3 (k * M in the DCC document)
296 REAL8 tmp_coef;
297
298 int i;
299
300 // Compute time intervals
301 t13 = t12 + t23;
302 t14 = t13 + t34;
303 t24 = t23 + t34;
304
305 // Compute denomenators
306 denom1 = t12 * t13 * t14;
307 denom2 = t12 * t23 * t24;
308 denom3 = t13 * t23 * t34;
309
310 // Compute dx_vec
311 dx_vec[0] = t45;
312 dx_vec[1] = t45 * t45 / 2.0;
313 dx_vec[2] = t45 * t45 * t45 / 3.0;
314 dx_vec[3] = dx_vec[1] * dx_vec[1];
315
316 // Column 1, which goes with dx_vec[1]:
317 A1 = -1.0 * t24 * t34 / denom1;
318 A2 = t14 * t34 / denom2;
319 A3 = -1.0 * t14 * t24 / denom3;
320 A4 = -1 * (A1 + A2 + A3);
321
322 // Column 2, which goes with dx_vec[2]:
323 B1 = -1.0 * (t24 + t34) / denom1;
324 B2 = (t14 + t34) / denom2;
325 B3 = -1.0 * (t14 + t24) / denom3;
326 B4 = -1 * (B1 + B2 + B3);
327
328 // Column 3, which goes with dx_vec[3]:
329 C1 = -1.0 / denom1;
330 C2 = 1.0 / denom2;
331 C3 = -1.0 / denom3;
332 C4 = -1 * (C1 + C2 + C3);
333
334 // Compute results
335 for (i=0; i<dim; i++) {
336 tmp_coef = k4[i] * dx_vec[0];
337 tmp_coef += (A1 * k1[i] + A2 * k2[i] + A3 * k3[i] + A4 * k4[i]) * dx_vec[1];
338 tmp_coef += (B1 * k1[i] + B2 * k2[i] + B3 * k3[i] + B4 * k4[i]) * dx_vec[2];
339 tmp_coef += (C1 * k1[i] + C2 * k2[i] + C3 * k3[i] + C4 * k4[i]) * dx_vec[3];
340 dy[i] = tmp_coef;
341 }
342}
343
344
345
346/**
347 * Computes all powers of z(t) = x(t) + i*y(t) needed for WignerDMatrices.
348 */
349static void ComplexPowers_Compute(ComplexPowers *cp, gsl_vector *x, gsl_vector *y) {
350 int i, j, power;
351
352 // z = x + iy
353 gsl_vector *tmp = gsl_vector_calloc(x->size);
354 gsl_vector *z_mag_sqr = gsl_vector_calloc(x->size);
355
356 // Set zero'th power
357 i = 2*cp->LMax;
358 gsl_vector_add_constant(cp->real_part[i], 1.0);
359
360 // Set first power
361 i += 1;
362 gsl_vector_add(cp->real_part[i], x);
363 gsl_vector_add(cp->imag_part[i], y);
364
365 // Compute positive powers
366 for (power=2; power <= 2*cp->LMax; power++) {
367 // Compute z^n = z^{n-1} * z. Currently, i indexes z^{n-1}.
368 // Re[z^{n-1]] * Re[z]
369 gsl_vector_add(tmp, cp->real_part[i]);
370 gsl_vector_mul(tmp, x);
371 gsl_vector_add(cp->real_part[i+1], tmp);
372 gsl_vector_set_zero(tmp);
373
374 // Re[z^{n-1]] * Im[z]
375 gsl_vector_add(tmp, cp->real_part[i]);
376 gsl_vector_mul(tmp, y);
377 gsl_vector_add(cp->imag_part[i+1], tmp);
378 gsl_vector_set_zero(tmp);
379
380 // Im[z^{n-1]] * Re[z]
381 gsl_vector_add(tmp, cp->imag_part[i]);
382 gsl_vector_mul(tmp, x);
383 gsl_vector_add(cp->imag_part[i+1], tmp);
384 gsl_vector_set_zero(tmp);
385
386 // Im[z^{n-1]] * Im[z]
387 gsl_vector_add(tmp, cp->imag_part[i]);
388 gsl_vector_mul(tmp, y);
389 gsl_vector_sub(cp->real_part[i+1], tmp);
390 gsl_vector_set_zero(tmp);
391
392 i += 1;
393 }
394
395 // Compute z^{-n} = (z^n)* / |(z^n)^2|
396 for (power=1; power <= 2*cp->LMax; power++) {
397 i = 2*cp->LMax + power;
398 j = 2*cp->LMax - power;
399
400 // Compute |(z^n)|^2
401 gsl_vector_add(tmp, cp->real_part[i]);
402 gsl_vector_mul(tmp, tmp);
403 gsl_vector_add(z_mag_sqr, tmp);
404 gsl_vector_set_zero(tmp);
405 gsl_vector_add(tmp, cp->imag_part[i]);
406 gsl_vector_mul(tmp, tmp);
407 gsl_vector_add(z_mag_sqr, tmp);
408 gsl_vector_set_zero(tmp);
409
410 // Set z^{-n}
411 gsl_vector_add(cp->real_part[j], cp->real_part[i]);
412 gsl_vector_div(cp->real_part[j], z_mag_sqr);
413 gsl_vector_sub(cp->imag_part[j], cp->imag_part[i]);
414 gsl_vector_div(cp->imag_part[j], z_mag_sqr);
415
416 gsl_vector_set_zero(z_mag_sqr);
417 }
418
419 gsl_vector_free(tmp);
420 gsl_vector_free(z_mag_sqr);
421
422}
423
424/**
425 * Computes all powers of x(t) needed for WignerDMatrices.
426 */
427static void RealPowers_Compute(RealPowers *rp, gsl_vector *x) {
428 int power;
429
430 gsl_vector *tmp = gsl_vector_calloc(x->size);
431
432 // Set zero'th power
433 gsl_vector_add_constant(rp->powers[0], 1.0);
434
435 // Set first power
436 gsl_vector_add(rp->powers[1], x);
437
438 // Compute positive powers
439 for (power=2; power <= 2*rp->LMax; power++) {
440 // Compute x^n = x^{n-1} * x.
441 gsl_vector_add(tmp, rp->powers[power-1]);
442 gsl_vector_mul(tmp, x);
443 gsl_vector_add(rp->powers[power], tmp);
444 gsl_vector_set_zero(tmp);
445 }
446
447 gsl_vector_free(tmp);
448}
449
450
451/**
452 * Transforms modes from the coorbital frame to the inertial frame.
453 * Note that this is a coorbital frame determined entirely from the waveform
454 * modes and *not* from the trajectories; it is a rotation of the coprecessing
455 * frame about the z-axis by an orbital phase computed from the waveform modes.
456 * See eqns 6-9 and the surrounding text of https://arxiv.org/abs/1705.07089
457 */
459 MultiModalWaveform *h, /**< Output. Dimensionless waveform modes.
460 Should be initialized already. */
461 MultiModalWaveform *h_coorb, /**< Coorbital frame waveform modes. */
462 gsl_vector **quat, /**< Coprecessing frame quaternions - 4 vectors. */
463 gsl_vector *orbphase /**< Orbital phase. */
464) {
465 int i, j, ell, m, mp;
466 int n_times = h_coorb->n_times;
467 int lmax = h_coorb->lvals[h_coorb->n_modes -1];
468 int lmin;
469 gsl_vector *cosmphi = gsl_vector_alloc(n_times);
470 gsl_vector *sinmphi = gsl_vector_alloc(n_times);
471 gsl_vector *tmp_vec = gsl_vector_alloc(n_times);
472
473 // First transform to the coprecessing frame:
474 // h^{\ell, m}_\mathrm{copr} = e^{-i m \varphi} h^{\ell, m}_\mathrm{coorb}.
475 // Loop over m first, since fixed m modes transform the same way.
476 MultiModalWaveform *h_copr = NULL;
477 MultiModalWaveform_Init(&h_copr, lmax, n_times);
478 for (m = -1*lmax; m <= lmax; m++) {
479 if (m==0) {
480 // No transformation
481 for (ell = 2; ell <= lmax; ell++) {
482 i = ell*(ell+1) - 4;
483 gsl_vector_add(h_copr->modes_real_part[i], h_coorb->modes_real_part[i]);
484 gsl_vector_add(h_copr->modes_imag_part[i], h_coorb->modes_imag_part[i]);
485 }
486 } else {
487 for (j=0; j<n_times; j++) {
488 gsl_vector_set(cosmphi, j, cos(m * gsl_vector_get(orbphase, j)));
489 gsl_vector_set(sinmphi, j, sin(m * gsl_vector_get(orbphase, j)));
490 }
491 lmin = abs(m);
492 if (lmin < 2) {
493 lmin = 2;
494 }
495 for (ell = lmin; ell<=lmax; ell++) {
496 i = ell*(ell+1) - 4 + m;
497 // compute and add combinations of {real, imag} and {cosmphi, sinmphi}
498 // real * cosmphi
499 gsl_vector_set_zero(tmp_vec);
500 gsl_vector_add(tmp_vec, h_coorb->modes_real_part[i]);
501 gsl_vector_mul(tmp_vec, cosmphi);
502 gsl_vector_add(h_copr->modes_real_part[i], tmp_vec);
503 // -1 * real * I * sinmphi
504 gsl_vector_set_zero(tmp_vec);
505 gsl_vector_add(tmp_vec, h_coorb->modes_real_part[i]);
506 gsl_vector_mul(tmp_vec, sinmphi);
507 gsl_vector_sub(h_copr->modes_imag_part[i], tmp_vec);
508 // I * imag * cosmphi
509 gsl_vector_set_zero(tmp_vec);
510 gsl_vector_add(tmp_vec, h_coorb->modes_imag_part[i]);
511 gsl_vector_mul(tmp_vec, cosmphi);
512 gsl_vector_add(h_copr->modes_imag_part[i], tmp_vec);
513 // imag * sinmphi
514 gsl_vector_set_zero(tmp_vec);
515 gsl_vector_add(tmp_vec, h_coorb->modes_imag_part[i]);
516 gsl_vector_mul(tmp_vec, sinmphi);
517 gsl_vector_add(h_copr->modes_real_part[i], tmp_vec);
518 }
519 }
520 }
521
522 // Now we transform the coprecessing modes to the intertial frame, using the quaternions.
523 WignerDMatrices *matrices = NULL;
524 WignerDMatrices_Init(&matrices, n_times, lmax);
525 WignerDMatrices_Compute(matrices, quat);
526 int matrix_index;
527 for (ell = 2; ell <= lmax; ell++) {
528 for (m= -1 * ell; m <= ell; m++) {
529 i = ell*(ell+1) - 4 + m;
530 for (mp = -1 * ell; mp <= ell; mp++) {
531 j = ell*(ell+1) - 4 + mp;
532 matrix_index = WignerDMatrix_Index(ell, m, mp);
533 // Re[h] * Re[D]
534 gsl_vector_set_zero(tmp_vec);
535 gsl_vector_add(tmp_vec, h_copr->modes_real_part[j]);
536 gsl_vector_mul(tmp_vec, matrices->real_part[matrix_index]);
537 gsl_vector_add(h->modes_real_part[i], tmp_vec);
538
539 // Re[h] * Im[D]
540 gsl_vector_set_zero(tmp_vec);
541 gsl_vector_add(tmp_vec, h_copr->modes_real_part[j]);
542 gsl_vector_mul(tmp_vec, matrices->imag_part[matrix_index]);
543 gsl_vector_add(h->modes_imag_part[i], tmp_vec);
544
545 // Im[h] * Re[D]
546 gsl_vector_set_zero(tmp_vec);
547 gsl_vector_add(tmp_vec, h_copr->modes_imag_part[j]);
548 gsl_vector_mul(tmp_vec, matrices->real_part[matrix_index]);
549 gsl_vector_add(h->modes_imag_part[i], tmp_vec);
550
551 // Im[h] * Im[D]
552 gsl_vector_set_zero(tmp_vec);
553 gsl_vector_add(tmp_vec, h_copr->modes_imag_part[j]);
554 gsl_vector_mul(tmp_vec, matrices->imag_part[matrix_index]);
555 gsl_vector_sub(h->modes_real_part[i], tmp_vec);
556 }
557 }
558 }
559 WignerDMatrices_Destroy(matrices);
561 gsl_vector_free(cosmphi);
562 gsl_vector_free(sinmphi);
563 gsl_vector_free(tmp_vec);
564}
565
566/**
567 * Computes all the Wigner D Matrices.
568 * We compute them all at once because a lot of the work is just processing quat.
569 * Parts of this function are adapted from GWFrames:
570 * https://github.com/moble/GWFrames
571 * written by Michael Boyle, based on his paper: http://arxiv.org/abs/1302.2919
572 * although we are working with the conjugate of quat.
573 */
575 WignerDMatrices *matrices, /**< The initialized WignerDMatrices which will contain the output.*/
576 gsl_vector **quat /**< The 4 quaternion components representing the coorbital frame.*/
577) {
578
579 int n_times = quat[0]->size;
580 int i, j, ell, m, mp, rho_min, rho_max, rho;
581 REAL8 tmp_re, tmp_im, tmp_abs_sqr, coef, c;
582 REAL8 eps_sqr = 1.0e-24;
583
584 // ra = q[0] - I*q[3], and rb = -q[2] - I*q[1]
585 // It's possible that |ra| = 0 or |rb| = 0, making it unsafe to divide by them.
586 // We keep track of which indices have this problem, use a safe-but-garbage value of quat=[0.5, 0.5, 0.5, 0.5] at those indices,
587 // and then fix them afterwards.
588 gsl_vector_int *index_types = gsl_vector_int_calloc(n_times);
589 gsl_vector *safe_ra_mag_sqr = gsl_vector_calloc(n_times);
590 gsl_vector *safe_rb_mag_sqr = gsl_vector_calloc(n_times);
591 gsl_vector *safe_ra_real = gsl_vector_calloc(n_times);
592 gsl_vector *safe_ra_imag = gsl_vector_calloc(n_times);
593 gsl_vector *safe_rb_real = gsl_vector_calloc(n_times);
594 gsl_vector *safe_rb_imag = gsl_vector_calloc(n_times);
595 gsl_vector_add(safe_ra_real, quat[0]);
596 gsl_vector_sub(safe_ra_imag, quat[3]);
597 gsl_vector_sub(safe_rb_real, quat[2]);
598 gsl_vector_sub(safe_rb_imag, quat[1]);
599 for (i=0; i<n_times; i++) {
600 tmp_re = gsl_vector_get(quat[2], i);
601 tmp_im = gsl_vector_get(quat[1], i);
602 tmp_abs_sqr = tmp_re*tmp_re + tmp_im*tmp_im;
603 if (tmp_abs_sqr < eps_sqr) {
604 gsl_vector_int_set(index_types, i, 2); // Use 0 for normal, 1 for ra small, 2 for ra ok but rb small
605 gsl_vector_set(safe_rb_mag_sqr, i, 0.5);
606 gsl_vector_set(safe_rb_real, i, 0.25);
607 gsl_vector_set(safe_rb_imag, i, 0.25);
608 } else {
609 gsl_vector_set(safe_rb_mag_sqr, i, tmp_abs_sqr);
610 }
611
612 tmp_re = gsl_vector_get(quat[0], i);
613 tmp_im = gsl_vector_get(quat[3], i);
614 tmp_abs_sqr = tmp_re*tmp_re + tmp_im*tmp_im;
615 if (tmp_abs_sqr < eps_sqr) {
616 gsl_vector_int_set(index_types, i, 1); // Use 0 for normal, 1 for ra small, 2 for ra ok but rb small
617 gsl_vector_set(safe_ra_mag_sqr, i, 0.5);
618 gsl_vector_set(safe_ra_real, i, 0.25);
619 gsl_vector_set(safe_ra_imag, i, 0.25);
620 } else {
621 gsl_vector_set(safe_ra_mag_sqr, i, tmp_abs_sqr);
622 }
623 }
624
625 // We will need various integer powers of ra, rb, |ra^2|, and |(rb/ra)^2| that will be used in many WignerD terms
626 gsl_vector *safe_abs_rb_over_ra_sqr = gsl_vector_calloc(n_times);
627 gsl_vector_add(safe_abs_rb_over_ra_sqr, safe_rb_mag_sqr);
628 gsl_vector_div(safe_abs_rb_over_ra_sqr, safe_ra_mag_sqr);
629
630 ComplexPowers *ra_powers = NULL;
631 ComplexPowers *rb_powers = NULL;
632 RealPowers *abs_ra_sqr_powers = NULL;
633 RealPowers *abs_rb_over_ra_sqr_powers = NULL;
634
635 ComplexPowers_Init(&ra_powers, matrices->LMax, n_times);
636 ComplexPowers_Init(&rb_powers, matrices->LMax, n_times);
637 RealPowers_Init(&abs_ra_sqr_powers, matrices->LMax, n_times);
638 RealPowers_Init(&abs_rb_over_ra_sqr_powers, matrices->LMax, n_times);
639
640
641 ComplexPowers_Compute(ra_powers, safe_ra_real, safe_ra_imag);
642 ComplexPowers_Compute(rb_powers, safe_rb_real, safe_rb_imag);
643 RealPowers_Compute(abs_ra_sqr_powers, safe_ra_mag_sqr);
644 RealPowers_Compute(abs_rb_over_ra_sqr_powers, safe_abs_rb_over_ra_sqr);
645
646
647 // Compute the matrices. We don't use safe_* anymore, so we can use them as temporary results
648 for (ell=2; ell <= matrices->LMax; ell++) {
649 for (m = -1*ell; m <= ell; m++) {
650 for (mp = -1*ell; mp <= ell; mp++) {
651 i = WignerDMatrix_Index(ell, m, mp);
652 coef = wigner_coef(ell, mp, m);
653 gsl_vector_set_zero(safe_ra_mag_sqr);
654
655 // factor = coef * (ra^(m+mp)) * (rb^(m-mp)) * (|ra^2|^(ell-m))
656 // The result is
657 // factor * \sum_{rho} c(rho) * (|(rb/ra)^2|^rho)
658 // where c(rho) = (-1)^rho * binom(ell+mp, rho) * binom(ell-mp, ell-rho-m)
659
660 // Store factor in matrices, using safe_rb_real and safe_rb_imag as temporary results
661 gsl_vector_add(matrices->real_part[i], ra_powers->real_part[2*matrices->LMax + m + mp]);
662 gsl_vector_add(matrices->imag_part[i], ra_powers->imag_part[2*matrices->LMax + m + mp]);
663 complex_vector_mult(matrices->real_part[i], matrices->imag_part[i],
664 rb_powers->real_part[2*matrices->LMax + m - mp], rb_powers->imag_part[2*matrices->LMax + m - mp],
665 safe_rb_real, safe_rb_imag);
666 gsl_vector_scale(matrices->real_part[i], coef);
667 gsl_vector_scale(matrices->imag_part[i], coef);
668 gsl_vector_mul(matrices->real_part[i], abs_ra_sqr_powers->powers[ell-m]);
669 gsl_vector_mul(matrices->imag_part[i], abs_ra_sqr_powers->powers[ell-m]);
670
671 // Compute the sum, storing it in safe_ra_mag_sqr
672 rho_min = (mp > m) ? (mp - m) : 0;
673 rho_max = (mp > -1*m) ? (ell-m) : (ell+mp);
674 for (rho=rho_min; rho <= rho_max; rho++) {
675 c = binomial(ell+mp, rho) * binomial(ell-mp, ell-rho-m);
676 if ((rho%2)==1) {
677 c *= -1;
678 }
679 // Store the temporary term in safe_rb_mag_sqr
680 gsl_vector_set_zero(safe_rb_mag_sqr);
681 gsl_vector_add(safe_rb_mag_sqr, abs_rb_over_ra_sqr_powers->powers[rho]);
682 gsl_vector_scale(safe_rb_mag_sqr, c);
683 gsl_vector_add(safe_ra_mag_sqr, safe_rb_mag_sqr);
684 }
685
686 // multiply by the (real) sum to get the result
687 gsl_vector_mul(matrices->real_part[i], safe_ra_mag_sqr);
688 gsl_vector_mul(matrices->imag_part[i], safe_ra_mag_sqr);
689
690 }
691 }
692 }
693
694 REAL8 zx, zy;
695 int k;
696 // Now fix the values at bad indices
697 for (j=0; j<n_times; j++) {
698 rho = gsl_vector_int_get(index_types, j);
699 if (rho != 0) {
700 for (ell=2; ell <= matrices->LMax; ell++) {
701 for (m=-1*ell; m <= ell; m++) {
702 for (mp = -1*ell; mp <= ell; mp++) {
703 i = WignerDMatrix_Index(ell, m, mp);
704 zx = 0.0;
705 zy = 0.0;
706 if (rho==1 && (mp == -m)) {
707 // z = (-1)^(ell+m+1) * rb^(2*m)
708 tmp_re = -1 * gsl_vector_get(quat[2], j);
709 tmp_im = -1 * gsl_vector_get(quat[1], j);
710 zx = 1.0;
711 zy = 0.0;
712 for (k=0; k < 2*m; k++) {
713 c = zx * tmp_im;
714 zx = zx*tmp_re - zy*tmp_im;
715 zy = zy*tmp_re + c;
716 }
717 if ((ell+m)%2 == 0) {
718 zx *= -1;
719 zy *= -1;
720 }
721 gsl_vector_set(matrices->real_part[i], j, zx);
722 gsl_vector_set(matrices->imag_part[i], j, zy);
723 } else if (rho==2 && (mp == m)) {
724 // z = ra^(2*m)
725 tmp_re = gsl_vector_get(quat[0], j);
726 tmp_im = -1 * gsl_vector_get(quat[3], j);
727 zx = 1.0;
728 zy = 0.0;
729 for (k=0; k < 2*m; k++) {
730 c = zx * tmp_im;
731 zx = zx*tmp_re - zy*tmp_im;
732 zy = zy*tmp_re + c;
733 }
734 gsl_vector_set(matrices->real_part[i], j, zx);
735 gsl_vector_set(matrices->imag_part[i], j, zy);
736 }
737 gsl_vector_set(matrices->real_part[i], j, zx);
738 gsl_vector_set(matrices->imag_part[i], j, zy);
739 }
740 }
741 }
742 }
743 }
744
745 // Cleanup
746 gsl_vector_int_free(index_types);
747 gsl_vector_free(safe_ra_mag_sqr);
748 gsl_vector_free(safe_rb_mag_sqr);
749 gsl_vector_free(safe_ra_real);
750 gsl_vector_free(safe_ra_imag);
751 gsl_vector_free(safe_rb_real);
752 gsl_vector_free(safe_rb_imag);
753 gsl_vector_free(safe_abs_rb_over_ra_sqr);
754 ComplexPowers_Destroy(ra_powers);
755 ComplexPowers_Destroy(rb_powers);
756 RealPowers_Destroy(abs_ra_sqr_powers);
757 RealPowers_Destroy(abs_rb_over_ra_sqr_powers);
758}
759
760
761/**
762 * Computes inverse, qInv, of a quaternion q such that q*qInv = 1.
763 */
764UNUSED static int quatInv(
765 REAL8 *qInv, /**< Output: Inverse of q. Should be initialized with size=4 */
766 REAL8 *q /**< Input quaternion. Should have size=4*/
767) {
768 const REAL8 normSqr = q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3];
769 qInv[0] = q[0]/normSqr;
770 for (int i=1; i<4; i++) {
771 qInv[i] = -q[i]/normSqr;
772 }
773
774 return XLAL_SUCCESS;
775}
776
777/**
778 * Multiplies two quaternions.
779 */
780UNUSED static int multiplyQuats(
781 REAL8 *prod, /**< Output: Product of q1 and q2. Should be initialized with size=4. */
782 REAL8 *q1, /**< First quaternion. Should have size=4. */
783 REAL8 *q2 /**< Second quaternion. Should have size=4. */
784) {
785 prod[0] = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3];
786 prod[1] = q1[2]*q2[3] - q2[2]*q1[3] + q1[0]*q2[1] + q2[0]*q1[1];
787 prod[2] = q1[3]*q2[1] - q2[3]*q1[1] + q1[0]*q2[2] + q2[0]*q1[2];
788 prod[3] = q1[1]*q2[2] - q2[1]*q1[2] + q1[0]*q2[3] + q2[0]*q1[3];
789 return XLAL_SUCCESS;
790}
791
792
793/**
794 * Given a coprecessing frame quaternion (quat), transforms a 3-vector (vec)
795 * from the coprecessing frame to the inertial frame.
796 * Reference: Eq.(A6) of https://arxiv.org/pdf/1302.2919.pdf
797 */
799 REAL8 *new_vec, /**< Output: Transformed vector in inertial frame. Should be initialized with size=3 */
800 REAL8 *quat, /**< Coprecessing frame quaternion. Should have size=4.*/
801 REAL8 *vec /**< Input vector in coprecessing frame. Should have size=3.*/
802) {
803
804 // quaternion representation of vec, such that the scalar component is zero
805 REAL8 vec_quat[4] = {0, vec[0], vec[1], vec[2]};
806
807 // q_prod = vec_quat * inv(quat)
808 REAL8 q_prod[4];
809 REAL8 qInv[4];
810 quatInv(qInv, quat);
811 multiplyQuats(q_prod, vec_quat, qInv);
812
813 // quaternion representation of output vec, such that the
814 // scalar component is zero
815 REAL8 new_vec_quat[4];
816
817 // new_vec_quat = quat * vec_quat * inv(quat)
818 multiplyQuats(new_vec_quat, quat, q_prod);
819 new_vec[0] = new_vec_quat[1];
820 new_vec[1] = new_vec_quat[2];
821 new_vec[2] = new_vec_quat[3];
822
823 return XLAL_SUCCESS;
824}
825
826/**
827 * Transforms a vector from the coprecessing frame to the inertial frame
828 * using the coprecessing frame quaternions.
829 */
831 gsl_vector **vec, /**< Input and Output: 3 time-dependent components of vec in the coprecessing frame */
832 gsl_vector **quat /**< 4 time-dependent components of coprecessing frame quaternions. */
833) {
834
835 REAL8 tmp_vec[3];
836 REAL8 tmp_new_vec[3];
837 REAL8 tmp_quat[4];
838
839 REAL8 *x = vec[0]->data;
840 REAL8 *y = vec[1]->data;
841 REAL8 *z = vec[2]->data;
842
843 REAL8 *q0 = quat[0]->data;
844 REAL8 *q1 = quat[1]->data;
845 REAL8 *q2 = quat[2]->data;
846 REAL8 *q3 = quat[3]->data;
847
848 int n = quat[0]->size;
849 for (int i=0; i<n; i++) {
850
851 tmp_vec[0] = x[i];
852 tmp_vec[1] = y[i];
853 tmp_vec[2] = z[i];
854
855 tmp_quat[0] = q0[i];
856 tmp_quat[1] = q1[i];
857 tmp_quat[2] = q2[i];
858 tmp_quat[3] = q3[i];
859
860 // update x,y,z, which also updates vec, with the transformed vector
861 quaternionTransformVector(tmp_new_vec, tmp_quat, tmp_vec);
862 x[i] = tmp_new_vec[0];
863 y[i] = tmp_new_vec[1];
864 z[i] = tmp_new_vec[2];
865 }
866
867 return XLAL_SUCCESS;
868}
869
870
871/**
872 * Wrapper to get dimensionless omegaMin/omegaRef from fMin/fRef.
873 */
874UNUSED static int get_dimless_omega(
875 REAL8 *omegaMin_dimless, /**< Output. omegaMin in units rad/M. */
876 REAL8 *omegaRef_dimless, /**< Output. omegaRef in units rad/M. */
877 const REAL8 fMin, /**< fMin in Hertz. */
878 const REAL8 fRef, /**< fRef in Hertz. */
879 const REAL8 Mtot_sec /**< Total mass in seconds. */
880) {
881
882 // Orbital angular frequency = 0.5*(wave angular frequency) = pi*(wave frequency)
883 *omegaMin_dimless = (Mtot_sec * fMin) * LAL_PI;
884
885 if (fRef == 0) {
886 // If fRef is 0, set it to fMin
887 *omegaRef_dimless = *omegaMin_dimless;
888 }
889 else {
890 *omegaRef_dimless = (Mtot_sec * fRef) * LAL_PI;
891 }
892
893 if (*omegaRef_dimless + 1e-13 < *omegaMin_dimless){
894 XLAL_ERROR(XLAL_EINVAL, "fRef cannot be lesser than fMin.");
895 }
896
897 return XLAL_SUCCESS;
898}
#define LMax
#define c
static REAL8 UNUSED q3(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
static REAL8 UNUSED C2(REAL8 e0, REAL8 f0)
static REAL8 UNUSED C1(REAL8 Mtotal)
static REAL8 UNUSED C4(REAL8 e0, REAL8 f0)
static REAL8 UNUSED q2(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
static REAL8 UNUSED q1(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
static REAL8 UNUSED C3(REAL8 e0, REAL8 f0)
static void WignerDMatrices_Compute(WignerDMatrices *matrices, gsl_vector **quat)
Computes all the Wigner D Matrices.
static int WignerDMatrix_Index(int ell, int m, int mp)
Computes the index of the (ell, m, mp) component in a WignerDMatices structure.
static void NRSur_ab4_dy(REAL8 *dy, REAL8 *k1, REAL8 *k2, REAL8 *k3, REAL8 *k4, REAL8 t12, REAL8 t23, REAL8 t34, REAL8 t45, int dim)
static UNUSED int quatInv(REAL8 *qInv, REAL8 *q)
Computes inverse, qInv, of a quaternion q such that q*qInv = 1.
static void ComplexPowers_Init(ComplexPowers **cp, int LMax, int n_times)
Initialize a ComplexPowers structure.
static void MultiModalWaveform_Destroy(MultiModalWaveform *wave)
Destroy a MultiModalWaveform structure; free all associated memory.
static void WignerDMatrices_Destroy(WignerDMatrices *matrices)
Destroy a WignerDMatrices structure; free all associated memory.
static void complex_vector_mult(gsl_vector *x1, gsl_vector *y1, gsl_vector *x2, gsl_vector *y2, gsl_vector *tmpx, gsl_vector *tmpy)
Multiplies z1(t) * z2(t) = (x1(t) + i*y1(t)) * (x2(t) + i*y2(t)), storing the result in x1 and y1.
static void RealPowers_Compute(RealPowers *rp, gsl_vector *x)
Computes all powers of x(t) needed for WignerDMatrices.
static void TransformModesCoorbitalToInertial(MultiModalWaveform *h, MultiModalWaveform *h_coorb, gsl_vector **quat, gsl_vector *orbphase)
Transforms modes from the coorbital frame to the inertial frame.
static void MultiModalWaveform_Init(MultiModalWaveform **wave, int LMax, int n_times)
Initialize a MultiModalWaveforme.
static void WignerDMatrices_Init(WignerDMatrices **matrices, int n_times, int LMax)
Initialize a WignerDMatrices structure.
static UNUSED int multiplyQuats(REAL8 *prod, REAL8 *q1, REAL8 *q2)
Multiplies two quaternions.
static void ComplexPowers_Compute(ComplexPowers *cp, gsl_vector *x, gsl_vector *y)
Computes all powers of z(t) = x(t) + i*y(t) needed for WignerDMatrices.
static UNUSED int get_dimless_omega(REAL8 *omegaMin_dimless, REAL8 *omegaRef_dimless, const REAL8 fMin, const REAL8 fRef, const REAL8 Mtot_sec)
Wrapper to get dimensionless omegaMin/omegaRef from fMin/fRef.
static void ComplexPowers_Destroy(ComplexPowers *cp)
Destroy a ComplexPowers structure; free all associated memory.
static UNUSED int transformTimeDependentVector(gsl_vector **vec, gsl_vector **quat)
Transforms a vector from the coprecessing frame to the inertial frame using the coprecessing frame qu...
static UNUSED int quaternionTransformVector(REAL8 *new_vec, REAL8 *quat, REAL8 *vec)
Given a coprecessing frame quaternion (quat), transforms a 3-vector (vec) from the coprecessing frame...
static void RealPowers_Destroy(RealPowers *rp)
Destroy a RealPowers structure; free all associated memory.
static REAL8 wigner_coef(int ell, int mp, int m)
static REAL8 factorial(int n)
static REAL8 factorial_ratio(int n, int k)
static void RealPowers_Init(RealPowers **rp, int LMax, int n_times)
Initialize a RealPowers structure.
static REAL8 binomial(int n, int k)
double i
Definition: bh_ringdown.c:118
double e
Definition: bh_ringdown.c:117
#define LAL_PI
double REAL8
void * XLALMalloc(size_t n)
void * XLALCalloc(size_t m, size_t n)
void XLALFree(void *p)
static const INT4 m
static const INT4 q
#define XLAL_ERROR_VOID(...)
#define XLAL_ERROR(...)
XLAL_SUCCESS
XLAL_EINVAL
XLAL_FAILURE
list y
x
Helper structure for computing WignerDMatrices, which require z(t)^n for many values of n.
int n_entries
Number of entries in **powers.
gsl_vector ** imag_part
The imag part of z^n.
int LMax
Maximum L; powers computed depend on LMax.
gsl_vector ** real_part
The real part of z^n.
Structure for a multi-modal waveform.
gsl_vector ** modes_real_part
The real part of each mode - n_modes vectors of length n_times.
int * mvals
The m values of each mode (n_modes entries)
int * lvals
The ell values of each mode (n_modes entries)
gsl_vector ** modes_imag_part
The imaginary part of each mode - n_modes vectors of length n_times.
int n_modes
Number of modes.
int n_times
Number of time samples in each mode.
Helper structure for computing WignerDMatrices, which require x(t)^n for many values of n.
gsl_vector ** powers
The data; x^n.
int LMax
Maximum L; powers computed depend on LMax.
int n_entries
Number of entries in **powers.
int LMax
Includes (ell, m) modes with 2 <= ell <= LMax.
gsl_vector ** imag_part
The imaginary part of each entry.
int n_entries
Number of (ell, m, m') entries.
gsl_vector ** real_part
The real part of each entry.