Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-8a6b96f
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALSimIMRSEOBNRROMUtilities.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2014, 2015 Michael Puerrer, John Veitch
3 * Reduced Order Model for SEOBNR
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/**
22 * \author Michael Puerrer
23 *
24 * \file
25 *
26 * \brief Auxiliary functions for SEOBNRv1/v2 reduced order modeling codes in
27 * LALSimIMRSEOBNRv2ROMDoubleSpin.c, LALSimIMRSEOBNRv1ROMDoubleSpin.c,
28 * LALSimIMRSEOBNRv2ROMEffectiveSpin.c, LALSimIMRSEOBNRv1ROMEffectiveSpin.c,
29 * LALSimIMRSEOBNRv2ChirpTime.c, LALSimIMRSEOBNRv2ROMDoubleSpinHI.c.
30 *
31 * Here I collect common auxiliary functions pertaining to reading
32 * data stored in gsl binary vectors and matrices, HDF5 files using the LAL interface,
33 * parameter space interpolation with B-splines, fitting to a cubic,
34 * a custom gsl error handler and adjustment of nearby parameter values.
35 */
36
37#include <stdio.h>
38#include <stdarg.h>
39#include <lal/XLALError.h>
40#include <stdbool.h>
41#include <gsl/gsl_math.h>
42#include <gsl/gsl_multifit.h>
43#include <gsl/gsl_interp.h>
44#include <gsl/gsl_fit.h>
46#include "gsl_bspline_basis.h"
47
48#ifdef LAL_HDF5_ENABLED
49#include <lal/H5FileIO.h>
50#endif
51
52UNUSED static int read_vector(const char dir[], const char fname[], gsl_vector *v);
53UNUSED static int read_matrix(const char dir[], const char fname[], gsl_matrix *m);
54/* SEOBNRv4HM_ROM functions */
55UNUSED static char* concatenate_strings(int count, ...);
56UNUSED static REAL8 sigmoid(REAL8 x);
57UNUSED static UINT4 blend(gsl_vector * freqs, gsl_vector * out_fun, REAL8 freq_1, REAL8 freq_2);
58UNUSED static UINT4 blend_functions(gsl_vector * freqs_out, gsl_vector * out_fun, gsl_vector * freq_in_1, gsl_vector * fun_in_1, gsl_vector * freq_in_2, gsl_vector * fun_in_2, REAL8 freq_1, REAL8 freq_2);
59UNUSED static UINT4 compute_i_max_LF_i_min_HF(INT8 * i_max_LF, INT8 * i_min_LF,gsl_vector * freqs_in_1,gsl_vector * freqs_in_2,REAL8 freq_1);
60UNUSED static REAL8 Get_omegaQNM_SEOBNRv4(REAL8 q, REAL8 chi1z, REAL8 chi2z, UINT4 l, UINT4 m);
61UNUSED static REAL8 Get_omegaQNM_SEOBNRv5(REAL8 q, REAL8 chi1z, REAL8 chi2z, UINT4 l, UINT4 m);
62UNUSED static UINT4 unwrap_phase(gsl_vector* phaseout,gsl_vector* phasein);
63UNUSED static UINT8 compute_i_at_f(gsl_vector * freq_array, REAL8 freq);
64UNUSED static UINT4 align_wfs_window(gsl_vector* f_array_1, gsl_vector* f_array_2, gsl_vector* phase_1, gsl_vector* phase_2, REAL8* Deltat, REAL8* Deltaphi, REAL8 f_align_start, REAL8 f_align_end);
65UNUSED static UINT4 align_wfs_window_from_22(gsl_vector* f_array_1, gsl_vector* f_array_2, gsl_vector* phase_1, gsl_vector* phase_2, REAL8 f_align_start, REAL8 f_align_end, REAL8 Deltat22, REAL8 Deltaphi22, INT4 modeM);
66
67#ifdef LAL_HDF5_ENABLED
68UNUSED static int CheckVectorFromHDF5(LALH5File *file, const char name[], const double *v, size_t n);
69UNUSED static int ReadHDF5RealVectorDataset(LALH5File *file, const char *name, gsl_vector **data);
70UNUSED static int ReadHDF5RealMatrixDataset(LALH5File *file, const char *name, gsl_matrix **data);
71UNUSED static int ReadHDF5LongVectorDataset(LALH5File *file, const char *name, gsl_vector_long **data);
72UNUSED static int ReadHDF5LongMatrixDataset(LALH5File *file, const char *name, gsl_matrix_long **data);
73UNUSED static void PrintInfoStringAttribute(LALH5File *file, const char attribute[]);
74#endif
75
77 gsl_vector *v,
78 REAL8 eta,
79 REAL8 chi1,
80 REAL8 chi2,
81 int ncy,
82 int ncz,
83 gsl_bspline_workspace *bwx,
84 gsl_bspline_workspace *bwy,
85 gsl_bspline_workspace *bwz
86);
87
89 gsl_vector *v,
90 REAL8 eta,
91 REAL8 chi,
92 int ncx,
93 int ncy,
94 gsl_bspline_workspace *bwx,
95 gsl_bspline_workspace *bwy
96);
97
98UNUSED static gsl_vector *Fit_cubic(const gsl_vector *xi, const gsl_vector *yi);
99
100UNUSED static bool approximately_equal(REAL8 x, REAL8 y, REAL8 epsilon);
101UNUSED static void nudge(REAL8 *x, REAL8 X, REAL8 epsilon);
102
103UNUSED static double SEOBNRROM_Ringdown_Mf_From_Mtot_q(
104 const double Mtot_sec,
105 const double q,
106 const double chi1,
107 const double chi2,
108 Approximant apx
109);
110UNUSED static double SEOBNRROM_Ringdown_Mf_From_Mtot_Eta(
111 const double Mtot_sec,
112 const double eta,
113 const double chi1,
114 const double chi2,
115 Approximant apx
116);
117
118// Helper functions to read gsl_vector and gsl_matrix data with error checking
119static int read_vector(const char dir[], const char fname[], gsl_vector *v) {
120 size_t size = strlen(dir) + strlen(fname) + 2;
121 char *path = XLALMalloc(size);
122 snprintf(path, size, "%s/%s", dir, fname);
123
124 FILE *f = fopen(path, "rb");
125 if (!f)
126 XLAL_ERROR(XLAL_EIO, "Could not find ROM data file at path `%s'", path);
127 int ret = gsl_vector_fread(f, v);
128 if (ret != 0)
129 XLAL_ERROR(XLAL_EIO, "Error reading data from `%s'", path);
130 fclose(f);
131
132 XLAL_PRINT_INFO("Sucessfully read data file `%s'", path);
133 XLALFree(path);
134 return(XLAL_SUCCESS);
135}
136
137static int read_matrix(const char dir[], const char fname[], gsl_matrix *m) {
138 size_t size = strlen(dir) + strlen(fname) + 2;
139 char *path = XLALMalloc(size);
140 snprintf(path, size, "%s/%s", dir, fname);
141
142 FILE *f = fopen(path, "rb");
143 if (!f)
144 XLAL_ERROR(XLAL_EIO, "Could not find ROM data file at path `%s'", path);
145 int ret = gsl_matrix_fread(f, m);
146 if (ret != 0)
147 XLAL_ERROR(XLAL_EIO, "Error reading data from `%s'", path);
148 fclose(f);
149
150 XLAL_PRINT_INFO("Sucessfully read data file `%s'", path);
151 XLALFree(path);
152 return(XLAL_SUCCESS);
153}
154
155#ifdef LAL_HDF5_ENABLED
156static int CheckVectorFromHDF5(LALH5File *file, const char name[], const double *v, size_t n) {
157 gsl_vector *temp = NULL;
158 ReadHDF5RealVectorDataset(file, name, &temp);
159
160 if (temp->size != n)
161 XLAL_ERROR(XLAL_EIO, "Length of data %s disagrees", name);
162
163 for (size_t i=0; i<temp->size; i++)
164 if (gsl_vector_get(temp, i) != v[i])
165 XLAL_ERROR(XLAL_EIO, "Data %s disagrees", name);
166 gsl_vector_free(temp);
167 return XLAL_SUCCESS;
168}
169
170static int ReadHDF5RealVectorDataset(LALH5File *file, const char *name, gsl_vector **data) {
171 LALH5Dataset *dset;
172 UINT4Vector *dimLength;
173 size_t n;
174
175 if (file == NULL || name == NULL || data == NULL)
177
178 dset = XLALH5DatasetRead(file, name);
179 if (dset == NULL)
181
183 XLALH5DatasetFree(dset);
184 XLAL_ERROR(XLAL_ETYPE, "Dataset `%s' is wrong type", name);
185 }
186
187 dimLength = XLALH5DatasetQueryDims(dset);
188 if (dimLength == NULL) {
189 XLALH5DatasetFree(dset);
191 }
192 if (dimLength->length != 1) {
193 XLALH5DatasetFree(dset);
194 XLAL_ERROR(XLAL_EDIMS, "Dataset `%s' must be 1-dimensional", name);
195 }
196
197 n = dimLength->data[0];
198 XLALDestroyUINT4Vector(dimLength);
199
200 if (*data == NULL) {
201 *data = gsl_vector_alloc(n);
202 if (*data == NULL) {
203 XLALH5DatasetFree(dset);
204 XLAL_ERROR(XLAL_ENOMEM, "gsl_vector_alloc(%zu) failed", n);
205 }
206 }
207 else if ((*data)->size != n) {
208 XLALH5DatasetFree(dset);
209 XLAL_ERROR(XLAL_EINVAL, "Expected gsl_vector `%s' of size %zu", name, n);
210 }
211
212 // Now read the data
213 if (XLALH5DatasetQueryData((*data)->data, dset) < 0) {
214 XLALH5DatasetFree(dset);
216 }
217
218 XLALH5DatasetFree(dset);
219 return 0;
220}
221
222static int ReadHDF5RealMatrixDataset(LALH5File *file, const char *name, gsl_matrix **data) {
223 LALH5Dataset *dset;
224 UINT4Vector *dimLength;
225 size_t n1, n2;
226
227 if (file == NULL || name == NULL || data == NULL)
229
230 dset = XLALH5DatasetRead(file, name);
231 if (dset == NULL)
233
235 XLALH5DatasetFree(dset);
236 XLAL_ERROR(XLAL_ETYPE, "Dataset `%s' is wrong type", name);
237 }
238
239 dimLength = XLALH5DatasetQueryDims(dset);
240 if (dimLength == NULL) {
241 XLALH5DatasetFree(dset);
243 }
244 if (dimLength->length != 2) {
245 XLALH5DatasetFree(dset);
246 XLAL_ERROR(XLAL_EDIMS, "Dataset `%s' must be 2-dimensional", name);
247 }
248
249 n1 = dimLength->data[0];
250 n2 = dimLength->data[1];
251 XLALDestroyUINT4Vector(dimLength);
252
253 if (*data == NULL) {
254 *data = gsl_matrix_alloc(n1, n2);
255 if (*data == NULL) {
256 XLALH5DatasetFree(dset);
257 XLAL_ERROR(XLAL_ENOMEM, "gsl_matrix_alloc(%zu, %zu) failed", n1, n2);
258 }
259 }
260 else if ((*data)->size1 != n1 || (*data)->size2 != n2) {
261 XLALH5DatasetFree(dset);
262 XLAL_ERROR(XLAL_EINVAL, "Expected gsl_matrix `%s' of size %zu x %zu", name, n1, n2);
263 }
264
265 // Now read the data
266 if (XLALH5DatasetQueryData((*data)->data, dset) < 0) {
267 XLALH5DatasetFree(dset);
269 }
270
271 XLALH5DatasetFree(dset);
272 return 0;
273}
274
275static int ReadHDF5LongVectorDataset(LALH5File *file, const char *name, gsl_vector_long **data) {
276 LALH5Dataset *dset;
277 UINT4Vector *dimLength;
278 size_t n;
279
280 if (file == NULL || name == NULL || data == NULL)
282
283 dset = XLALH5DatasetRead(file, name);
284 if (dset == NULL)
286
288 XLALH5DatasetFree(dset);
289 XLAL_ERROR(XLAL_ETYPE, "Dataset `%s' is wrong type", name);
290 }
291
292 dimLength = XLALH5DatasetQueryDims(dset);
293 if (dimLength == NULL) {
294 XLALH5DatasetFree(dset);
296 }
297 if (dimLength->length != 1) {
298 XLALH5DatasetFree(dset);
299 XLAL_ERROR(XLAL_EDIMS, "Dataset `%s' must be 1-dimensional", name);
300 }
301
302 n = dimLength->data[0];
303 XLALDestroyUINT4Vector(dimLength);
304
305 if (*data == NULL) {
306 *data = gsl_vector_long_alloc(n);
307 if (*data == NULL) {
308 XLALH5DatasetFree(dset);
309 XLAL_ERROR(XLAL_ENOMEM, "gsl_vector_long_alloc(%zu) failed", n);
310 }
311 }
312 else if ((*data)->size != n) {
313 XLALH5DatasetFree(dset);
314 XLAL_ERROR(XLAL_EINVAL, "Expected gsl_vector `%s' of size %zu", name, n);
315 }
316
317 // Now read the data
318 if (XLALH5DatasetQueryData((*data)->data, dset) < 0) {
319 XLALH5DatasetFree(dset);
321 }
322
323 XLALH5DatasetFree(dset);
324 return 0;
325}
326
327static int ReadHDF5LongMatrixDataset(LALH5File *file, const char *name, gsl_matrix_long **data) {
328 LALH5Dataset *dset;
329 UINT4Vector *dimLength;
330 size_t n1, n2;
331
332 if (file == NULL || name == NULL || data == NULL)
334
335 dset = XLALH5DatasetRead(file, name);
336 if (dset == NULL)
338
340 XLALH5DatasetFree(dset);
341 XLAL_ERROR(XLAL_ETYPE, "Dataset `%s' is wrong type", name);
342 }
343
344 dimLength = XLALH5DatasetQueryDims(dset);
345 if (dimLength == NULL) {
346 XLALH5DatasetFree(dset);
348 }
349 if (dimLength->length != 2) {
350 XLALH5DatasetFree(dset);
351 XLAL_ERROR(XLAL_EDIMS, "Dataset `%s' must be 2-dimensional", name);
352 }
353
354 n1 = dimLength->data[0];
355 n2 = dimLength->data[1];
356 XLALDestroyUINT4Vector(dimLength);
357
358 if (*data == NULL) {
359 *data = gsl_matrix_long_alloc(n1, n2);
360 if (*data == NULL) {
361 XLALH5DatasetFree(dset);
362 XLAL_ERROR(XLAL_ENOMEM, "gsl_matrix_long_alloc(%zu, %zu) failed", n1, n2);
363 }
364 }
365 else if ((*data)->size1 != n1 || (*data)->size2 != n2) {
366 XLALH5DatasetFree(dset);
367 XLAL_ERROR(XLAL_EINVAL, "Expected gsl_matrix_long `%s' of size %zu x %zu", name, n1, n2);
368 }
369
370 // Now read the data
371 if (XLALH5DatasetQueryData((*data)->data, dset) < 0) {
372 XLALH5DatasetFree(dset);
374 }
375
376 XLALH5DatasetFree(dset);
377 return 0;
378}
379
380static void PrintInfoStringAttribute(LALH5File *file, const char attribute[]) {
381 LALH5Generic gfile = {.file = file};
382 int len = XLALH5AttributeQueryStringValue(NULL, 0, gfile, attribute) + 1;
383 char *str = XLALMalloc(len);
384 XLALH5FileQueryStringAttributeValue(str, len, file, attribute);
385 XLALPrintInfo("%s:\n%s\n", attribute, str);
386 LALFree(str);
387}
388#endif
389
390// Helper function to perform tensor product spline interpolation with gsl
391// The gsl_vector v contains the ncx x ncy x ncz dimensional coefficient tensor in vector form
392// that should be interpolated and evaluated at position (eta,chi1,chi2).
394 gsl_vector *v,
395 REAL8 eta,
396 REAL8 chi1,
397 REAL8 chi2,
398 int ncy,
399 int ncz,
400 gsl_bspline_workspace *bwx,
401 gsl_bspline_workspace *bwy,
402 gsl_bspline_workspace *bwz
403) {
404 // Store nonzero cubic (order k=4) B-spline basis functions in the eta and chi directions.
405 gsl_vector *Bx4 = gsl_vector_alloc(4);
406 gsl_vector *By4 = gsl_vector_alloc(4);
407 gsl_vector *Bz4 = gsl_vector_alloc(4);
408
409 size_t isx, isy, isz; // first non-zero spline
410 // Evaluate all potentially nonzero cubic B-spline basis functions for
411 // positions (eta,chi) and store them in the vectors Bx4, By4, Bz4.
412 // Since the B-splines are of compact support we only need to store a small
413 // number of basis functions to avoid computing terms that would be zero anyway.
414 // https://www.gnu.org/software/gsl/manual/html_node/Overview-of-B_002dsplines.html#Overview-of-B_002dsplines
415 gsl_bspline_basis(eta, Bx4, &isx, bwx);
416 gsl_bspline_basis(chi1, By4, &isy, bwy);
417 gsl_bspline_basis(chi2, Bz4, &isz, bwz);
418
419 // Now compute coefficient at desired parameters (q,chi1,chi2)
420 // from C(eta,chi1,chi2) = c_ijk * Beta_i * Bchi1_j * Bchi2_k
421 // while summing over indices i,j,k where the B-splines are nonzero.
422 // Note: in the 2D case we were able to use gsl_matrix c = gsl_matrix_view_vector(&v, ncx, ncy).matrix
423 // to convert vector view of the coefficient matrix to a matrix view.
424 // However, since tensors are not supported in gsl, we have to do the indexing explicitly.
425 double sum = 0;
426 for (int i=0; i<4; i++)
427 for (int j=0; j<4; j++)
428 for (int k=0; k<4; k++) {
429 int ii = isx + i;
430 int jj = isy + j;
431 int kk = isz + k;
432 double cijk = gsl_vector_get(v, (ii*ncy + jj)*ncz + kk);
433 sum += cijk * gsl_vector_get(Bx4, i) * gsl_vector_get(By4, j) * gsl_vector_get(Bz4, k);
434 }
435
436 gsl_vector_free(Bx4);
437 gsl_vector_free(By4);
438 gsl_vector_free(Bz4);
439
440 return sum;
441}
442
443// Helper function to perform tensor product spline interpolation with gsl
444// The gsl_vector v contains the ncx x ncy dimensional coefficient matrix in vector form
445// that should be interpolated and evaluated at position (eta,chi).
447 gsl_vector *v,
448 REAL8 eta,
449 REAL8 chi,
450 int ncx,
451 int ncy,
452 gsl_bspline_workspace *bwx,
453 gsl_bspline_workspace *bwy
454) {
455 gsl_matrix c = gsl_matrix_view_vector(v, ncx, ncy).matrix; // Convert coefficient matrix from vector view to matrix view c_ij.
456
457 // Store nonzero cubic (order k=4) B-spline basis functions in the eta and chi directions.
458 gsl_vector *Bx4 = gsl_vector_alloc(4);
459 gsl_vector *By4 = gsl_vector_alloc(4);
460
461 REAL8 sum = 0;
462 size_t isx, isy; // first non-zero spline
463 // Evaluate all potentially nonzero cubic B-spline basis functions for positions (eta,chi) and stores them in the vectors Bx4, By4.
464 // Since the B-splines are of compact support we only need to store a small number of basis functions
465 // to avoid computing terms that would be zero anyway.
466 // https://www.gnu.org/software/gsl/manual/html_node/Overview-of-B_002dsplines.html#Overview-of-B_002dsplines
467 gsl_bspline_basis(eta, Bx4, &isx, bwx);
468 gsl_bspline_basis(chi, By4, &isy, bwy);
469
470 // Now compute coefficient at desired parameters (eta,chi) from C(eta,chi) = c_ij * Beta_i * Bchi_j
471 // summing over indices i,j where the B-splines are nonzero.
472 for (int i=0; i<4; i++)
473 for (int j=0; j<4; j++)
474 sum += gsl_matrix_get(&c, isx + i, isy + j) * gsl_vector_get(Bx4, i) * gsl_vector_get(By4, j);
475
476 gsl_vector_free(Bx4);
477 gsl_vector_free(By4);
478
479 return sum;
480}
481
482// Returns fitting coefficients for cubic y = c[0] + c[1]*x + c[2]*x**2 + c[3]*x**3
483static gsl_vector *Fit_cubic(const gsl_vector *xi, const gsl_vector *yi) {
484 const int n = xi->size; // how many data points are we fitting
485 const int p = 4; // order of fitting polynomial
486
487 gsl_matrix *X = gsl_matrix_alloc(n, p); // design matrix
488 gsl_vector *y = gsl_vector_alloc(n); // data vector
489 gsl_vector *c = gsl_vector_alloc(p); // coefficient vector
490 gsl_matrix *cov = gsl_matrix_alloc(p, p);
491 double chisq;
492
493 for (int i=0; i < n; i++) {
494 double xval = gsl_vector_get(xi, i);
495 double xval2 = xval*xval;
496 double xval3 = xval2*xval;
497 gsl_matrix_set(X, i, 0, 1.0);
498 gsl_matrix_set(X, i, 1, xval);
499 gsl_matrix_set(X, i, 2, xval2);
500 gsl_matrix_set(X, i, 3, xval3);
501
502 double yval = gsl_vector_get(yi, i);
503 gsl_vector_set (y, i, yval);
504 }
505
506 gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(n, p);
507 gsl_multifit_linear(X, y, c, cov, &chisq, work); // perform linear least-squares fit
508
509 // clean up
510 gsl_multifit_linear_free(work);
511 gsl_matrix_free(X);
512 gsl_matrix_free(cov);
513 gsl_vector_free(y);
514
515 return c;
516}
517
518// This function determines whether x and y are approximately equal to a relative accuracy epsilon.
519// Note that x and y are compared to relative accuracy, so this function is not suitable for testing whether a value is approximately zero.
520static bool approximately_equal(REAL8 x, REAL8 y, REAL8 epsilon) {
521 return !gsl_fcmp(x, y, epsilon);
522}
523
524// If x and X are approximately equal to relative accuracy epsilon then set x = X.
525// If X = 0 then use an absolute comparison.
526static void nudge(REAL8 *x, REAL8 X, REAL8 epsilon) {
527 if (X != 0.0) {
528 if (approximately_equal(*x, X, epsilon)) {
529 XLAL_PRINT_INFO("Nudging value %.15g to %.15g\n", *x, X);
530 *x = X;
531 }
532 }
533 else {
534 if (fabs(*x - X) < epsilon)
535 *x = X;
536 }
537}
538
539/* compute ringdown frequency for (2,2) mode from total mass and mass ratio */
541 const double Mtot_sec,
542 const double q,
543 const double chi1,
544 const double chi2,
545 Approximant apx
546)
547{
548 double Mtot_SI = Mtot_sec / LAL_MTSUN_SI * LAL_MSUN_SI;
549 double m1_SI = Mtot_SI * q / (1. + q);
550 double m2_SI = Mtot_SI * 1. / (1. + q);
551 return XLALSimInspiralGetFinalFreq(m1_SI, m2_SI, 0, 0, chi1, 0, 0, chi2, apx) * Mtot_sec;
552}
553
554/* compute ringdown frequency for (2,2) mode from total mass and eta */
556 const double Mtot_sec,
557 const double eta,
558 const double chi1,
559 const double chi2,
560 Approximant apx
561)
562{
563 double q = (1. + sqrt(1. - 4. * eta) - 2. * eta) / (2. * eta);
564 return SEOBNRROM_Ringdown_Mf_From_Mtot_q(Mtot_sec, q, chi1, chi2, apx);
565}
566
567/* SEOBNRv4HM_ROM functions */
568
569
570// Helper function to concatenate strings
571char* concatenate_strings(int count, ...)
572{
573 va_list ap;
574 int i;
575
576 // Find required length to store merged string
577 int len = 1; // room for NULL
578 va_start(ap, count);
579 for(i=0 ; i<count ; i++)
580 len += strlen(va_arg(ap, char*));
581 va_end(ap);
582
583 // Allocate memory to concat strings
584 char *merged = XLALCalloc(sizeof(char),len);
585 int null_pos = 0;
586
587 // Actually concatenate strings
588 va_start(ap, count);
589 for(i=0 ; i<count ; i++)
590 {
591 char *s = va_arg(ap, char*);
592 strcpy(merged+null_pos, s);
593 null_pos += strlen(s);
594 }
595 va_end(ap);
596
597 return merged;
598}
599
600/* Sigmoid function for hybridization */
602{
603 REAL8 exp_value;
604 REAL8 return_value;
605
606 /*** Exponential calculation ***/
607 exp_value = exp(-x);
608
609 /*** Final sigmoid value ***/
610 return_value = 1 / (1 + exp_value);
611
612 return return_value;
613}
614
615// Compute activation function in a window [freq_1, freq_2]
616UINT4 blend(gsl_vector * freqs, gsl_vector * out_fun, REAL8 freq_1, REAL8 freq_2){
617 for(unsigned int i = 0; i < freqs->size; i++){
618 if(freqs->data[i] <= freq_1){
619 out_fun->data[i] = 0.;
620 }
621 else if((freqs->data[i]>freq_1) && (freqs->data[i]<freq_2)){
622 out_fun->data[i] = sigmoid(- (freq_2 - freq_1)/(freqs->data[i]-freq_1) - (freq_2 - freq_1)/(freqs->data[i]-freq_2));
623 }
624 else{
625 out_fun->data[i] = 1.;
626 }
627 }
628 return XLAL_SUCCESS;
629}
630
631// Function to compute the indices associated to [freq_1] in two arrays
632UINT4 compute_i_max_LF_i_min_HF(INT8 * i_max_LF, INT8 * i_min_HF,gsl_vector * freqs_in_1,gsl_vector * freqs_in_2, REAL8 freq_1){
633
634 for(unsigned int i = 0; i < freqs_in_1->size; i++){
635 if(freqs_in_1->data[i]<freq_1){
636 *i_max_LF = i;
637 }
638 }
639 for(unsigned int i = freqs_in_2->size -1; i > 0; i--){
640 if(freqs_in_2->data[i]>=freq_1){
641 *i_min_HF = i;
642 }
643 }
644return XLAL_SUCCESS;
645}
646
647// Function to compute the index associated to [freq_1] in an array
648UINT8 compute_i_at_f(gsl_vector * freq_array, REAL8 freq){
649 UINT8 index = 0;
650 REAL8 diff = fabs(freq-freq_array->data[index]);
651 for(unsigned int i = 1; i < freq_array->size; i++){
652 if(fabs(freq-freq_array->data[i])< diff){
653 diff = fabs(freq-freq_array->data[i]);
654 index = i;
655 }
656 }
657 return index;
658}
659
660/* Function to blend LF and HF ROM functions in a window [freq_1, freq_2] */
661UINT4 blend_functions(gsl_vector * freqs_out, gsl_vector * out_fun, gsl_vector * freq_in_1, gsl_vector * fun_in_1, gsl_vector * freq_in_2, gsl_vector * fun_in_2, REAL8 freq_1, REAL8 freq_2){
662
663 gsl_vector *bld_2 = gsl_vector_calloc(out_fun->size);
664
665 /* Generate the window function */
666 UNUSED UINT4 ret = blend(freqs_out, bld_2,freq_1,freq_2);
667
668 gsl_interp_accel *acc = gsl_interp_accel_alloc ();
669 gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, freq_in_1->size);
670 gsl_spline_init (spline, freq_in_1->data, fun_in_1->data, freq_in_1->size);
671
672 unsigned int i = 0;
673 REAL8 inter_func;
674 /* Build the contribution to out_fun from fun_in_1 in the range [f_ini, freq_2] */
675 /* where f_ini is the starting frequency of the array freq_in_1 */
676 /* for freq > freq_2 the contribution of fun_in_1 to out_fun is null since 1 - bld_2 = 0 there*/
677 /* for freq < freq_1, bld_2 = 0, so I could have another easier while in the range [freq_1,freq_2]*/
678
679 while(freqs_out->data[i]<=freq_2){
680 /* Interpolate the original function on the new grid */
681 inter_func = gsl_spline_eval (spline, freqs_out->data[i], acc);
682 /* Multiply the original function by the window function */
683 out_fun->data[i] = inter_func*(1.-bld_2->data[i]);
684 i++;
685 }
686 gsl_spline_free (spline);
687 gsl_interp_accel_free(acc);
688
689 /* Build the contribution to out_fun from fun_in_2 in the range [freq_2, f_end] */
690 /* where f_end is the ending frequency of the array freq_in_2 */
691 /* for freq > freq_2, bld_2 = 1*/
692 spline = gsl_spline_alloc (gsl_interp_cspline, freq_in_2->size);
693 acc = gsl_interp_accel_alloc ();
694 gsl_spline_init (spline, freq_in_2->data, fun_in_2->data, freq_in_2->size);
695 i = freqs_out->size-1;
696 while(freqs_out->data[i]>freq_2){
697 inter_func = gsl_spline_eval (spline, freqs_out->data[i], acc);
698 out_fun->data[i] = inter_func;
699 i--;
700 }
701
702 /* Build the contribution to out_fun from fun_in_2 in the range [freq_1, freq_2] */
703 /* where f_end is the ending frequency of the array freq_in_2 */
704 while((freqs_out->data[i]>=freq_1)&&(freqs_out->data[i]<=freq_2)){
705 inter_func = gsl_spline_eval (spline, freqs_out->data[i], acc);
706 out_fun->data[i] += inter_func*bld_2->data[i];
707 /* Note that the += is important here */
708 i--;
709 }
710
711 /* Cleanup */
712 gsl_vector_free(bld_2);
713 gsl_spline_free (spline);
714 gsl_interp_accel_free(acc);
715
716 return XLAL_SUCCESS;
717}
718
719// Function to align two phases in a window [f_align_start,f_align_end]
720// Output phase shift, time shift from the linear fit of phi2 - phi1
721// Convention for alignment quantities:
722// phase_2 ~= phase_1 + Deltaphi + 2pi*f*Deltat
723UINT4 align_wfs_window(gsl_vector* f_array_1, gsl_vector* f_array_2, gsl_vector* phase_1, gsl_vector* phase_2, REAL8* Deltat, REAL8* Deltaphi, REAL8 f_align_start, REAL8 f_align_end){
724
725 /* Check frequency range */
726 if((f_align_start<f_array_1->data[0]) || (f_align_start<f_array_2->data[0]) || (f_align_end>f_array_1->data[f_array_1->size-1]) || (f_align_end>f_array_2->data[f_array_2->size-1])) {
727 XLAL_ERROR(XLAL_EINVAL, "Incompatible frequency ranges.");
728 }
729 /* Number of points to fit the phase difference */
730 /* arbitrary, suitable for linear fit */
731 UINT8 npt_fit = 10;
732 // Interpolate phase_1 phase on phase_2 frequency grid between [f_align_start,f_align_end] and compute phase difference
733 gsl_interp_accel *acc_1 = gsl_interp_accel_alloc();
734 gsl_spline *spline_1 = gsl_spline_alloc(gsl_interp_cspline, phase_1->size);
735 gsl_spline_init(spline_1, f_array_1->data, phase_1->data, phase_1->size);
736 gsl_interp_accel *acc_2 = gsl_interp_accel_alloc();
737 gsl_spline *spline_2 = gsl_spline_alloc(gsl_interp_cspline, phase_2->size);
738 gsl_spline_init(spline_2, f_array_2->data, phase_2->data, phase_2->size);
739 gsl_vector* delta_phi = gsl_vector_alloc(npt_fit);
740 gsl_vector* freq_delta_phi = gsl_vector_alloc(npt_fit);
741 // delta_phi is phi2-phi1
742 REAL8 f = 0.;
743 for(unsigned int i=0; i<freq_delta_phi->size; i++){
744 f = f_align_start + ((double) i)/(freq_delta_phi->size - 1) * (f_align_end-f_align_start);
745 freq_delta_phi->data[i] = f;
746 delta_phi->data[i] = gsl_spline_eval(spline_2, f, acc_2) - gsl_spline_eval(spline_1, f, acc_1);
747 }
748 REAL8 c0, c1, cov00, cov01, cov11, chisq;
749 gsl_fit_linear(freq_delta_phi->data, 1, delta_phi->data, 1, freq_delta_phi->size, &c0, &c1, &cov00, &cov01, &cov11, &chisq);
750 *Deltaphi = c0;
751 *Deltat = c1 / (2*LAL_PI);
752
753 // Remove the linear fit from phase_1
754 for(unsigned int i = 0; i < f_array_1->size; i++){
755 phase_1->data[i] = phase_1->data[i] + (c1*f_array_1->data[i] + c0);
756 }
757 /* Cleanup */
758 gsl_vector_free(delta_phi);
759 gsl_vector_free(freq_delta_phi);
760 gsl_spline_free(spline_1);
761 gsl_interp_accel_free(acc_1);
762 gsl_spline_free(spline_2);
763 gsl_interp_accel_free(acc_2);
764
765 return XLAL_SUCCESS;
766
767}
768
769// Useful debugging functions
770UNUSED static void test_save_gsl_vector(const char filename[], gsl_vector *v);
771UNUSED static void test_save_gsl_spline(const char filename[], gsl_spline *s);
772UNUSED static void test_save_cmplx_freq_series(COMPLEX16FrequencySeries *mode, const char filename[], bool save_real_part);
773
774/********************* Definitions begin here ********************/
775
776UNUSED static void test_save_gsl_vector(const char filename[], gsl_vector *v) {
777 FILE *fp = fopen(filename, "w");
778 gsl_vector_fprintf(fp, v, "%g");
779 fclose(fp);
780}
781
782UNUSED static void test_save_gsl_spline(const char filename[], gsl_spline *s) {
783 FILE *fp = fopen(filename, "w");
784 for (size_t i=0; i < s->size; i++) {
785 fprintf(fp, "%g\t%g\n", s->x[i], s->y[i]);
786 }
787 fclose(fp);
788}
789
792 const char filename[],
793 bool save_real_part)
794{
795 // Save either the real or imaginary part of a COMPLEX16FrequencySeries to a file.
796 gsl_vector *v = gsl_vector_alloc(mode->data->length);
797 double (*part_fun_ptr)(double complex);
798 if (save_real_part)
799 part_fun_ptr = creal;
800 else
801 part_fun_ptr = cimag;
802 for (unsigned int i=0; i<mode->data->length; i++)
803 gsl_vector_set(v, i, part_fun_ptr(mode->data->data[i]));
805 gsl_vector_free(v);
806}
807
808// Function to align two phases with a given 22-mode constant and slope
809// using the window [f_align_start,f_align_end] to determine an additional pi
810// (in general there is a pi ambiguity in m/2*Deltaphi_22)
811// Convention for alignment quantities:
812// phase2_22 ~= phase1_22 + Deltaphi22_align + 2pi*f*Deltat22_align
813UINT4 align_wfs_window_from_22(gsl_vector* f_array_1, gsl_vector* f_array_2, gsl_vector* phase_1, gsl_vector* phase_2, REAL8 f_align_start, REAL8 f_align_end, REAL8 Deltat22, REAL8 Deltaphi22, INT4 modeM){
814
815 /* Check frequency range */
816 if((f_align_start<f_array_1->data[0]) || (f_align_start<f_array_2->data[0]) || (f_align_end>f_array_1->data[f_array_1->size-1]) || (f_align_end>f_array_2->data[f_array_2->size-1])) {
817 XLAL_ERROR(XLAL_EINVAL, "Incompatible frequency ranges.");
818 }
819 /* Number of points to fit the phase difference */
820 /* arbitrary, suitable for linear fit */
821 UINT8 npt_fit = 10;
822 // Interpolate phase_1 phase on phase_2 frequency grid between [f_align_start,f_align_end] and compute phase difference
823 gsl_interp_accel *acc_1 = gsl_interp_accel_alloc();
824 gsl_spline *spline_1 = gsl_spline_alloc(gsl_interp_cspline, phase_1->size);
825 gsl_spline_init(spline_1, f_array_1->data, phase_1->data, phase_1->size);
826 gsl_interp_accel *acc_2 = gsl_interp_accel_alloc();
827 gsl_spline *spline_2 = gsl_spline_alloc(gsl_interp_cspline, phase_2->size);
828 gsl_spline_init(spline_2, f_array_2->data, phase_2->data, phase_2->size);
829 gsl_vector* delta_phi = gsl_vector_alloc(npt_fit);
830 gsl_vector* freq_delta_phi = gsl_vector_alloc(npt_fit);
831 // delta_phi is phi2-phi1, taking into account the constant and linear term from alignment to be added
832 REAL8 f = 0.;
833 for(unsigned int i=0; i<freq_delta_phi->size; i++){
834 f = f_align_start + ((double) i)/(freq_delta_phi->size - 1) * (f_align_end-f_align_start);
835 delta_phi->data[i] = gsl_spline_eval(spline_2, f, acc_2) - gsl_spline_eval(spline_1, f, acc_1) - (2*LAL_PI*Deltat22*f + modeM/2.*Deltaphi22);
836 }
837
838 // Determine a possible additional pi-shift (pi-ambiguity in m/2*Deltaphi22)
839 REAL8 av_phase_diff = 0.;
840 for(unsigned int i = 0; i < freq_delta_phi->size; i++){
841 av_phase_diff += delta_phi->data[i];
842 }
843 av_phase_diff = av_phase_diff / (delta_phi->size);
844
845 // Closest multiple of pi to eliminate the residual phase diff
846 REAL8 pishift = floor((av_phase_diff + LAL_PI/2) / LAL_PI) * LAL_PI;
847
848 // Remove the scaled linear fit and pi shift from phase_1
849 for(unsigned int i = 0; i < f_array_1->size; i++){
850 phase_1->data[i] = phase_1->data[i] + (2*LAL_PI*Deltat22*f_array_1->data[i] + modeM/2.*Deltaphi22 + pishift);
851 }
852
853 /* Cleanup */
854 gsl_vector_free(delta_phi);
855 gsl_vector_free(freq_delta_phi);
856 gsl_spline_free(spline_1);
857 gsl_interp_accel_free(acc_1);
858 gsl_spline_free(spline_2);
859 gsl_interp_accel_free(acc_2);
860
861 return XLAL_SUCCESS;
862
863}
864
865// Function to compute the QNM frequency
867 // Total mass M is not important here, we return the QNM frequencies in geometric units
868 REAL8 M = 100.;
869 REAL8 Ms = M * LAL_MTSUN_SI;
870 REAL8 m1 = M * q/(1+q);
871 REAL8 m2 = M * 1/(1+q);
872 Approximant SpinAlignedEOBapproximant = SEOBNRv4;
873 COMPLEX16Vector modefreqVec;
874 COMPLEX16 modeFreq;
875 modefreqVec.length = 1;
876 modefreqVec.data = &modeFreq;
877 REAL8 spin1[3] = {0., 0., chi1z};
878 REAL8 spin2[3] = {0., 0., chi2z};
879 // XLALSimIMREOBGenerateQNMFreqV2 is returning QNM frequencies in SI units, we multiply by Ms to convert it in geometric units
880 UNUSED UINT4 ret = XLALSimIMREOBGenerateQNMFreqV2(&modefreqVec, m1, m2, spin1, spin2, l, m, 1, SpinAlignedEOBapproximant);
881 return Ms * creal(modefreqVec.data[0]);
882}
883
884// Function to compute the SEOBNRv5 QNM frequency
886 // Total mass M is not important here, we return the QNM frequencies in geometric units
887 REAL8 M = 100.;
888 REAL8 Ms = M * LAL_MTSUN_SI;
889 REAL8 m1 = M * q/(1+q);
890 REAL8 m2 = M * 1/(1+q);
891 Approximant SpinAlignedEOBapproximant = SEOBNRv5_ROM; // SEOBNRv5 uses different final state fits compared to SEOBNRv4
892 COMPLEX16Vector modefreqVec;
893 COMPLEX16 modeFreq;
894 modefreqVec.length = 1;
895 modefreqVec.data = &modeFreq;
896 REAL8 spin1[3] = {0., 0., chi1z};
897 REAL8 spin2[3] = {0., 0., chi2z};
898 // XLALSimIMREOBGenerateQNMFreqV5 is returning QNM frequencies in SI units, we multiply by Ms to convert it in geometric units
899 UNUSED UINT4 ret = XLALSimIMREOBGenerateQNMFreqV5(&modefreqVec, m1, m2, spin1, spin2, l, m, 1, SpinAlignedEOBapproximant);
900 return Ms * creal(modefreqVec.data[0]);
901}
902
903/* Function to unwrap phases mod 2pi - acts on a REAL8Vector representing the phase */
904/* FIXME: for long vectors there are small differences with the numpy unwrap function - to be checked */
906 gsl_vector* phaseout, /* Output: unwrapped phase vector, already allocated - can be the same as phasein, in which case unwrapping is done in place */
907 gsl_vector* phasein) /* Input: phase vector */
908{
909 int N = phasein->size;
910 double* p = phasein->data;
911 double* pmod = (double*) malloc(sizeof(double) * N);
912 int* jumps = (int*) malloc(sizeof(int) * N);
913 int* cumul = (int*) malloc(sizeof(int) * N);
914
915 /* Compute phase mod 2pi (shifted to be between -pi and pi) */
916 for(int i=0; i<N; i++) {
917 pmod[i] = p[i] - floor((p[i] + LAL_PI) / (2*LAL_PI))*(2*LAL_PI);
918 }
919
920 /* Identify jumps */
921 jumps[0] = 0;
922 double d = 0.;
923 for(int i=1; i<N; i++) {
924 d = pmod[i] - pmod[i-1];
925 if(d<-LAL_PI) jumps[i] = 1;
926 else if(d>LAL_PI) jumps[i] = -1;
927 else jumps[i] = 0;
928 }
929
930 /* Cumulative of the jump sequence */
931 int c = 0;
932 cumul[0] = 0;
933 for(int i=1; i<N; i++) {
934 c += jumps[i];
935 cumul[i] = c;
936 }
937
938 /* Correct mod 2pi phase series by the number of 2pi factor given by the cumulative of the jumps */
939 double* pout = phaseout->data;
940 for(int i=0; i<N; i++) {
941 pout[i] = pmod[i] + 2*LAL_PI*cumul[i];
942 }
943
944 /* Cleanup */
945 free(pmod);
946 free(jumps);
947 free(cumul);
948
949 return XLAL_SUCCESS;
950}
struct tagLALH5File LALH5File
struct tagLALH5Dataset LALH5Dataset
#define LALFree(p)
INT4 XLALSimIMREOBGenerateQNMFreqV5(COMPLEX16Vector *modefreqs, const REAL8 mass1, const REAL8 mass2, const REAL8 spin1[3], const REAL8 spin2[3], UINT4 l, INT4 m, UINT4 nmodes, Approximant approximant)
These functions generate the quasinormal mode frequencies for a black hole ringdown,...
const double c1
const double c0
static UNUSED void test_save_gsl_vector(const char filename[], gsl_vector *v)
static UNUSED UINT8 compute_i_at_f(gsl_vector *freq_array, REAL8 freq)
static UNUSED gsl_vector * Fit_cubic(const gsl_vector *xi, const gsl_vector *yi)
static UNUSED void nudge(REAL8 *x, REAL8 X, REAL8 epsilon)
static UNUSED UINT4 align_wfs_window_from_22(gsl_vector *f_array_1, gsl_vector *f_array_2, gsl_vector *phase_1, gsl_vector *phase_2, REAL8 f_align_start, REAL8 f_align_end, REAL8 Deltat22, REAL8 Deltaphi22, INT4 modeM)
static UNUSED REAL8 Get_omegaQNM_SEOBNRv5(REAL8 q, REAL8 chi1z, REAL8 chi2z, UINT4 l, UINT4 m)
static UNUSED double SEOBNRROM_Ringdown_Mf_From_Mtot_q(const double Mtot_sec, const double q, const double chi1, const double chi2, Approximant apx)
static UNUSED void test_save_gsl_spline(const char filename[], gsl_spline *s)
static UNUSED UINT4 compute_i_max_LF_i_min_HF(INT8 *i_max_LF, INT8 *i_min_LF, gsl_vector *freqs_in_1, gsl_vector *freqs_in_2, REAL8 freq_1)
static UNUSED REAL8 Interpolate_Coefficent_Tensor(gsl_vector *v, REAL8 eta, REAL8 chi1, REAL8 chi2, int ncy, int ncz, gsl_bspline_workspace *bwx, gsl_bspline_workspace *bwy, gsl_bspline_workspace *bwz)
static UNUSED REAL8 Get_omegaQNM_SEOBNRv4(REAL8 q, REAL8 chi1z, REAL8 chi2z, UINT4 l, UINT4 m)
static UNUSED UINT4 unwrap_phase(gsl_vector *phaseout, gsl_vector *phasein)
static UNUSED int read_vector(const char dir[], const char fname[], gsl_vector *v)
static UNUSED void test_save_cmplx_freq_series(COMPLEX16FrequencySeries *mode, const char filename[], bool save_real_part)
static UNUSED UINT4 align_wfs_window(gsl_vector *f_array_1, gsl_vector *f_array_2, gsl_vector *phase_1, gsl_vector *phase_2, REAL8 *Deltat, REAL8 *Deltaphi, REAL8 f_align_start, REAL8 f_align_end)
static UNUSED UINT4 blend(gsl_vector *freqs, gsl_vector *out_fun, REAL8 freq_1, REAL8 freq_2)
static UNUSED double SEOBNRROM_Ringdown_Mf_From_Mtot_Eta(const double Mtot_sec, const double eta, const double chi1, const double chi2, Approximant apx)
static UNUSED char * concatenate_strings(int count,...)
static UNUSED REAL8 sigmoid(REAL8 x)
static UNUSED UINT4 blend_functions(gsl_vector *freqs_out, gsl_vector *out_fun, gsl_vector *freq_in_1, gsl_vector *fun_in_1, gsl_vector *freq_in_2, gsl_vector *fun_in_2, REAL8 freq_1, REAL8 freq_2)
static UNUSED int read_matrix(const char dir[], const char fname[], gsl_matrix *m)
static UNUSED REAL8 Interpolate_Coefficent_Matrix(gsl_vector *v, REAL8 eta, REAL8 chi, int ncx, int ncy, gsl_bspline_workspace *bwx, gsl_bspline_workspace *bwy)
static UNUSED bool approximately_equal(REAL8 x, REAL8 y, REAL8 epsilon)
double XLALSimInspiralGetFinalFreq(REAL8 m1, REAL8 m2, const REAL8 S1x, const REAL8 S1y, const REAL8 S1z, const REAL8 S2x, const REAL8 S2y, const REAL8 S2z, Approximant approximant)
Function that gives the default ending frequencies of the given approximant.
static vector d(const double L_norm, const double J_norm, const vector roots)
Internal function that returns the coefficients "d_0", "d_2" and "d_4" from 1703.03967 corresponding ...
#define c
const char * name
static REAL8 sum(REAL8 *arr1, REAL8 *arr2, int k)
#define fprintf
int s
Definition: bh_qnmode.c:137
int l
Definition: bh_qnmode.c:135
REAL8 M
Definition: bh_qnmode.c:133
double i
Definition: bh_ringdown.c:118
sigmaKerr data[0]
void XLALH5DatasetFree(LALH5Dataset UNUSED *dset)
int XLALH5AttributeQueryStringValue(char UNUSED *value, size_t UNUSED size, const LALH5Generic UNUSED object, const char UNUSED *key)
LALTYPECODE XLALH5DatasetQueryType(LALH5Dataset UNUSED *dset)
LALH5Dataset * XLALH5DatasetRead(LALH5File UNUSED *file, const char UNUSED *name)
int XLALH5DatasetQueryData(void UNUSED *data, LALH5Dataset UNUSED *dset)
int XLALH5FileQueryStringAttributeValue(char UNUSED *value, size_t UNUSED size, LALH5File UNUSED *file, const char UNUSED *key)
UINT4Vector * XLALH5DatasetQueryDims(LALH5Dataset UNUSED *dset)
#define LAL_MSUN_SI
#define LAL_PI
#define LAL_MTSUN_SI
uint64_t UINT8
double complex COMPLEX16
double REAL8
int64_t INT8
uint32_t UINT4
int32_t INT4
LAL_I8_TYPE_CODE
LAL_D_TYPE_CODE
void * XLALMalloc(size_t n)
void * XLALCalloc(size_t m, size_t n)
void XLALFree(void *p)
INT4 XLALSimIMREOBGenerateQNMFreqV2(COMPLEX16Vector *modefreqs, const REAL8 mass1, const REAL8 mass2, const REAL8 spin1[3], const REAL8 spin2[3], UINT4 l, INT4 m, UINT4 nmodes, Approximant approximant)
These functions generate the quasinormal mode frequencies for a black hole ringdown.
Approximant
Enum that specifies the PN approximant to be used in computing the waveform.
@ SEOBNRv4
Spin nonprecessing EOBNR model v4.
@ SEOBNRv5_ROM
Time domain, precessing phenomenological IMR waveform model with subdominant modes ([arXiv: 20XY....
static const INT4 m
static const INT4 q
void XLALDestroyUINT4Vector(UINT4Vector *vector)
#define XLAL_PRINT_INFO(...)
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_ERROR(...)
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFUNC
XLAL_ETYPE
XLAL_EIO
XLAL_EDIMS
XLAL_EINVAL
static int gsl_bspline_basis(const double x, gsl_vector *Bk, size_t *istart, gsl_bspline_workspace *w)
list y
path
filename
int N
p
x
COMPLEX16Sequence * data
COMPLEX16 * data
UINT4 * data
FILE * fp
LALH5File * file