Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALSimIMRNRHybSur3dq8.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2018 Vijay Varma
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with with program; see the file COPYING. If not, write to the
16 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17 * MA 02110-1301 USA
18 */
19
20/**
21 * \author Vijay Varma
22 *
23 * \file
24 *
25 * \brief C code for NRHybSur3dq8 waveform model, an NR-hybrid surrogate model
26 * for aligned-spin BBH.
27 *
28 * The binary data file (NRHybSur3dq8_lal_v1.0.h5) is available at:
29 * https://git.ligo.org/waveforms/software/lalsuite-waveform-data
30 * and on Zenodo at: https://zenodo.org/records/14999310.
31 * Get the lalsuite-waveform-data repo or put the data into a location in your
32 * LAL_DATA_PATH.
33 * The data is also available on CIT at /home/lalsimulation_data and via CVMFS
34 * at /cvmfs/shared.storage.igwn.org/igwn/shared/auxiliary/obs_sci/cbc/waveform/lalsimulation_data
35 * Make sure the files are in your LAL_DATA_PATH.
36 *
37 * **Paper**: https://arxiv.org/abs/1812.07865
38 *
39 * **Parameter ranges**:
40 *
41 * q = [1, 10.1] and \f$\chi_{1z}, \chi_{2z}\f$ = [-0.81, 0.81]
42 * or
43 * q = [1, 9.1] and \f$\chi_{1z}, \chi_{2z}\f$ = [-0.91, 0.91]
44 *
45 * modes: \f$ \ell \leq 4, m \geq 0 \f$, and (5,5), but not (4,1) or (4,0).
46 * m<0 modes are determined from the m \f$\geq0\f$ modes.
47 *
48 * \f$M \geq 2.25 M_{\odot} \f$, for fstart=20Hz, for all modes.
49 *
50 * **Training parameter ranges**:
51 *
52 * q = [1, 8]
53 *
54 * \f$\chi_{1z}, \chi_{2z}\f$ = [-0.8, 0.8]
55 *
56 * But extrapolates reasonably to the above mass ratios and spins.
57 */
58
59
60#ifdef __GNUC__
61#define UNUSED __attribute__ ((unused))
62#else
63#define UNUSED
64#endif
65
68
69
70#include <libgen.h>
71
72#include <lal/FileIO.h>
73#include <lal/SphericalHarmonics.h>
74#include <lal/H5FileIO.h>
75
76
77#ifdef LAL_PTHREAD_LOCK
78#include <pthread.h>
79#endif
80
81#ifdef LAL_PTHREAD_LOCK
82static pthread_once_t NRHybSur3dq8_is_initialized = PTHREAD_ONCE_INIT;
83#endif
84
85/**
86 * Global surrogate data.
87 * This data will be loaded at most once. Any executable which calls
88 * NRHybSur3dq8_Init_LALDATA directly or by calling any XLAL NRHybSur3dq8
89 * function will have a memory leak according to valgrind, because we never
90 * free this memory.
91 */
93
94//*************************************************************************/
95//************************* function definitions **************************/
96//*************************************************************************/
97
98/**
99 * Helper function to check if the NRHybSur3dq8 model has been initialized.
100 */
101static bool NRHybSur3dq8_IsSetup(void) {
103 return true;
104 else
105 return false;
106}
107
108
109/**
110 * Surrogate initialization.
111 *
112 * This needs to be called once, before __lalsim_NRHybSur3dq8_data is used. It
113 * finds the H5 data file with the NRHybSur3dq8 data and loads the surrogate.
114 * Can be called multiple times, will just return true if surrogate is already
115 * loaded.
116 */
117static void NRHybSur3dq8_Init_LALDATA(void) {
118
119 if (NRHybSur3dq8_IsSetup()) return;
120
122 if (path==NULL){
124 "Unable to resolve data file '%s' in $LAL_DATA_PATH.\n"
125 "Note: LALSuite versions >= 7.25 require data files that are publicly available at:\n"
126 "https://git.ligo.org/waveforms/software/lalsuite-waveform-data\n"
127 "and on Zenodo at: https://zenodo.org/records/14999310.\n"
128 "For earlier LALSuite versions, use the files in lalsuite-extra, available at:\n"
129 "https://git.ligo.org/lscsoft/lalsuite-extra\n",
131 }
132
133
134 char *dir = dirname(path);
135 const UINT4 size = strlen(dir) + strlen(NRHybSur3dq8_DATAFILE) + 2;
136 char *file_path = XLALMalloc(size);
137 snprintf(file_path, size, "%s/%s", dir, NRHybSur3dq8_DATAFILE);
138
139 #if LAL_HDF5_ENABLED
140 LALH5File *file = XLALH5FileOpen(file_path, "r");
141 if (file==NULL) {
143 "Unable to load data file %s in $LAL_DATA_PATH."
144 " File may be corrupted.\n", NRHybSur3dq8_DATAFILE);
145 }
146
147 int ret = NRHybSur_Init(&__lalsim_NRHybSur3dq8_data, file);
148
149 ret |= ROM_check_canonical_file_basename(file,NRHybSur3dq8_DATAFILE,
150 "CANONICAL_FILE_BASENAME");
151
153 if (ret != XLAL_SUCCESS) {
154 XLAL_ERROR_VOID(XLAL_FAILURE, "Failure loading data from %s\n",
155 file_path);
156 }
157 #else
158 XLAL_ERROR_VOID(XLAL_FAILURE, "Failure loading data from %s\n. HDF5 support is not enabled.\n", file_path);
159 #endif
160
161 XLALFree(path);
162 XLALFree(file_path);
163}
164
165
166/**
167 * Map from mass ratio and spins to surrogate fit parameters.
168 *
169 * The fit parameters are \f$[log_e(q), \hat{\chi}, \chi_a]\f$.
170 * \f$\hat{\chi}\f$ is defined in Eq.(46) of arxiv:1812.07865.
171 * \f$\chi_a = (\chi_{1z} - \chi_{2z})/2 \f$.
172 * These are described in Sec.VI.C.3 of arxiv:1812.07865.
173 */
175 gsl_vector* fit_params, /**< Output: mapped fit parameters. */
176 const REAL8 q, /**< Mass ratio m1 / m2 >= 1. */
177 const REAL8 chi1z, /**< Dimless spin of heavier BH. */
178 const REAL8 chi2z /**< Dimless spin of lighter BH. */
179) {
180
181 const REAL8 eta = q/(1.+q)/(1.+q);
182 const REAL8 chi_wtAvg = (q*chi1z+chi2z)/(1.+q);
183 const REAL8 chiHat
184 = (chi_wtAvg - 38.*eta/113.*(chi1z + chi2z))/(1. - 76.*eta/113.);
185 const REAL8 chi_a = (chi1z - chi2z)/2.;
186
187 XLAL_CHECK((fit_params != NULL) && (fit_params->size == 3), XLAL_EDIMS,
188 "NRHybSur3dq8_fitParams(): size of fit_params should be 3, not %zu.\n",
189 fit_params->size);
190
191 gsl_vector_set(fit_params, 0, log(q));
192 gsl_vector_set(fit_params, 1, chiHat);
193 gsl_vector_set(fit_params, 2, chi_a);
194
195 return XLAL_SUCCESS;
196}
197
198
199/**
200 * This is the core function of the NRHybSur3dq8 model.
201 *
202 * It evaluates all required waveform modes. For each mode, it evaluates each
203 * waveform data piece. The different data pieces are described in Sec.VI.B of
204 * arxiv:1812.07865.
205 * For the (2,2) mode the data pieces are the amplitude and phase. Note that we
206 * model the phase residual but add back the 0PN term at evaluation time.
207 * For all other modes the data pieces are the real and imaginary parts of the
208 * strain in the coorbital frame.
209 *
210 * The reference point is the time at which the (2,2) mode frequency equals
211 * fRef. If fRef=0, sets fRef=fMin. We set the orbital phase to 0 at the
212 * reference point. The orbital phase is obtained as \f$\phi_{22}/2\f$, so this
213 * leaves a pi ambiguity. But the surrogate data is already aligned such that
214 * the heavier BH is on the +ve x-axis at t=-1000M. See Sec.VI.A.4 of
215 * arxiv:1812.07865. This resolves the pi ambiguity. This means that after the
216 * realignment, the orbital phase at reference frequency fRef is 0, or Bh1 is
217 * on the +ve x-axis. Note that this is alignment is done using only waveform
218 * quantities, so this doesn't necessarily correspond to the (gauge dependent)
219 * NR Bh trajectories. The modes are returned in this frame, which agrees with
220 * the LAL convention. When combining the modes to get the polarizations, the
221 * Ylms should be evaluated at (inclination, pi/2 - phiRef), following the LAL
222 * convention.
223 *
224 * Only uses data at (2,2) mode frequencies >= fMin. This determines the start
225 * time. The start time, along with the step size deltaT, is used to determine
226 * the output_times. Uses cubic spline interpolation to interpolate from the
227 * surrogate's time array to output_times.
228 *
229 * NOTE: If mass ratio q<1, the labels of the two BHs are swapped internally.
230 *
231 * **Returned values**:
232 *
233 * phi_22: contains the phase of the (2,2) mode. This is always evaluated as
234 * this is required for other modes as well to transform from coorbital
235 * frame to inertial frame.
236 *
237 * evaluated_mode_dps: Contains all other data pieces. This is a list of size
238 * num_modes_incl, the number of modes to include. For each mode this
239 * contains the amplitude, and real and imaginary parts of the coorbital
240 * frame strain. For the (2,2) mode only the amplitude is required. For
241 * all other modes only the coorbital frame strain is required. So, we
242 * evaluate only the required data pieces of each mode. The amplitude and
243 * coorbital frame strain is in units of rh/M and needs to be rescaled to
244 * physical units.
245 *
246 * epoch: Initial time value, w.r.t. the peak (t=0 at the peak) of the total
247 * waveform amplitude, as defined in Eq.38 of arxiv:1812.07865.
248 */
250 gsl_vector **phi_22, /**< Output: phase of (2,2) mode. */
251 EvaluatedDataPieces **evaluated_mode_dps, /**< Output: All other data
252 pieces. */
253 LIGOTimeGPS *epoch, /**< Output: Initial time value, where t=0 is at
254 the peak of the total waveform amplitude. */
255 const REAL8 deltaT, /**< Sampling interval (s). */
256 const REAL8 fMin, /**< Start GW frequency (Hz). */
257 const REAL8 fRef, /**< Reference GW frequency (Hz). */
258 REAL8 q, /**< Mass ratio m1/m2. */
259 const REAL8 Mtot_sec, /**< Total mass in geometric units (s). */
260 REAL8 chi1z, /**< Dimensionless spin of Bh1. */
261 REAL8 chi2z, /**< Dimensionless spin of Bh2. */
262 LALValue* ModeArray, /**< Container for (ell, m) modes to generate. */
263 LALDict* LALparams /**< Dict with extra parameters */
264) {
265
266 const NRHybSurData *NR_hybsur_data = &__lalsim_NRHybSur3dq8_data;
267
268 REAL8 init_orbphase = 0; // In the LAL convention the larger BH should
269 // be on the +ve x-axis at fRef, this is done
270 // by setting orbphase = 0 at fRef.
271 // For the surrogate Bh1 is defined to be the one with the larger mass,
272 // Bh2 with the smaller mass. Therefore, if q<1, we swap the labels of the
273 // two Bhs
274 if (q < 1) {
275 q = 1./q;
276 REAL8 tmp = chi1z;
277 chi1z = chi2z;
278 chi2z = tmp;
279
280 // The above swap generates the requested waveform but rotated by
281 // pi about the z-axis, we undo this by adding pi to init_orbphase.
282 init_orbphase += LAL_PI;
283 }
284
285 const char *param_validity = "This model is valid for q <= 9.1 & "
286 "|chi1z,chi2z| <= 0.91, or q <= 10.1 & |chi1z,chi2z| <= 0.81";
287
288 // By default we do not allow unlimited_extrapolation
289 UINT4 unlim_extrap = 0;
290 if (LALparams != NULL &&
291 XLALDictContains(LALparams, "unlimited_extrapolation")) {
292 // Unless the user asks for it
293 unlim_extrap
294 = XLALDictLookupUINT4Value(LALparams, "unlimited_extrapolation");
295 }
296
297 // Sanity checks and warnings
298 if ((q > 10.1) && (unlim_extrap == 0)) {
300 "Too much extrapolation in mass ratio; q=%0.4f > 10.1\n%s\n", q,
301 param_validity);
302 }
303 if (q > 8.1) {
305 "Extrapolating outside training range q=%0.4f > 8.1\n", q);
306 }
307 if ((fabs(chi1z) > 0.91) && (unlim_extrap == 0)) {
309 "Too much extrapolation; |chi1z|=%0.4f > 0.91\n%s\n", fabs(chi1z),
310 param_validity);
311 }
312 if ((fabs(chi2z) > 0.91) && (unlim_extrap == 0)) {
314 "Too much extrapolation; |chi2z|=%0.4f > 0.91\n%s\n", fabs(chi2z),
315 param_validity);
316 }
317 if (fabs(chi1z) > 0.81) {
318 if ((q > 9.1) && (unlim_extrap == 0)) {
320 "Too much extrapolation; q=%0.4f > 9.1 & |chi1z|=%.04f"
321 " >0.81\n%s\n", q, fabs(chi1z), param_validity);
322 }
324 "Extrapolating outside training range |chi1z|=%0.4f > 0.81\n",
325 fabs(chi1z));
326 }
327 if (fabs(chi2z) > 0.81) {
328 if ((q > 9.1) && (unlim_extrap == 0)) {
330 "Too much extrapolation; q=%0.4f > 9.1 & |chi2z|=%.04f"
331 " >0.81\n%s\n", q, fabs(chi2z), param_validity);
332 }
334 "Extrapolating outside training range |chi2z|=%0.4f > 0.81\n",
335 fabs(chi2z));
336 }
337
338 // Get dimensionless start and reference frequencies. Note that cases where
339 // fRef<fMin, deltaT is too small, and fRef or fMin are too
340 // small/large (i.e. outside of the surrogate range) are handled in
341 // NRHybSur_eval_phase_22().
342
343 // dimensionless start angular frequency of (2,2) mode in rad/M
344 const REAL8 omegaM_22_min = 2*LAL_PI*fMin*Mtot_sec;
345
346 // dimensionless reference angular frequency of (2,2) mode in rad/M
347 REAL8 omegaM_22_ref;
348 if (fRef == 0) {
349 // If fRef is 0, set it to fMin
350 omegaM_22_ref = omegaM_22_min;
351 }
352 else {
353 omegaM_22_ref = 2*LAL_PI*fRef*Mtot_sec;
354 }
355
356 // dimensionless time step size
357 const REAL8 deltaTOverM = deltaT/Mtot_sec;
358
359 // Compute fit_params (initialize with dimension of the surrogate)
360 gsl_vector *fit_params = gsl_vector_alloc(NR_hybsur_data->params_dim);
361 int ret = NRHybSur3dq8_fitParams(fit_params, q, chi1z, chi2z);
362 if(ret != XLAL_SUCCESS) {
363 XLAL_ERROR(XLAL_EFUNC, "Failed to evaluate fit_params.");
364 }
365
366 // Training data parameters, required for GPR fits
367 const gsl_matrix *x_train = NR_hybsur_data->x_train;
368
369 // sparse time array of the surrogate
370 const gsl_vector *domain = NR_hybsur_data->domain;
371
372 // assign size to dummy_worker and dummy_dp
373 gsl_vector *dummy_worker = gsl_vector_alloc(NR_hybsur_data->params_dim);
374 gsl_vector *dummy_dp = gsl_vector_alloc(domain->size);
375
376 gsl_vector *output_times = NULL;
377
378 // symmetric mass ratio
379 const REAL8 eta = q/(1.+q)/(1.+q);
380
381 // Evaluate phase of (2,2) mode
382 // The phase of the (2,2,) mode is always evaluated as this is
383 // needed for the transformation from coorbital to inertial frame for all
384 // modes
385 ret = NRHybSur_eval_phase_22(phi_22, &output_times, eta, fit_params,
386 omegaM_22_min, deltaTOverM, init_orbphase, omegaM_22_ref, dummy_dp,
387 x_train, dummy_worker, NR_hybsur_data);
388 if(ret != XLAL_SUCCESS) {
390 "Failed to evaluate phi_22 data piece");
391 }
392
393 // set epoch to initial time. Note that t=0 is at the peak of the total
394 // waveform amplitude, as defined in Eq.38 of arxiv:1812.07865.
395 const REAL8 t0 = gsl_vector_get(output_times, 0);
396 XLALGPSAdd(epoch, Mtot_sec * t0);
397
398 // Evaluate other data pieces for required modes
399 const gsl_matrix_long *mode_list = NR_hybsur_data->mode_list;
400 const UINT4 num_modes_modeled = NR_hybsur_data->num_modes_modeled;
401 UINT4 incl_mode_idx = 0; // This tracks the output modes
402 for (UINT4 mode_idx = 0; mode_idx < num_modes_modeled; mode_idx++){
403
404 const UINT4 ell = gsl_matrix_long_get(mode_list, mode_idx, 0);
405 const UINT4 m = gsl_matrix_long_get(mode_list, mode_idx, 1);
406
407 // Evaluate a mode only if it is required
408 if (XLALSimInspiralModeArrayIsModeActive(ModeArray, ell, m) == 1) {
409
410 const ModeDataPieces *data_pieces
411 = NR_hybsur_data->mode_data_pieces[mode_idx];
412
413 if((ell != data_pieces->ell) || (m != data_pieces->m)){
414 XLAL_ERROR(XLAL_EDATA, "Modes do not agree");
415 }
416
417 evaluated_mode_dps[incl_mode_idx]
420
422 &evaluated_mode_dps[incl_mode_idx], ell, m,
423 data_pieces, output_times, fit_params, dummy_dp,
424 x_train, dummy_worker, NR_hybsur_data);
425 if(ret != XLAL_SUCCESS) {
427 "Failed to evaluate (%u, %u) mode", ell, m);
428 }
429
430 incl_mode_idx += 1;
431 }
432 }
433
434 gsl_vector_free(fit_params);
435 gsl_vector_free(dummy_dp);
436 gsl_vector_free(dummy_worker);
437 gsl_vector_free(output_times);
438
439 return XLAL_SUCCESS;
440}
441
442
443
444
445
446
447//************************ main functions ***************************/
448
449/**
450 * Reference: arxiv:1812.07865.
451 *
452 * Returns the plus and cross polarizations of NRHybSur3dq8 waveform model.
453 *
454 * The reference point is the time at which the (2,2) mode frequency equals
455 * fRef. If fRef=0, sets fRef=fMin. We set the orbital phase to 0 at the
456 * reference point. The orbital phase is obtained as \f$\phi_{22}/2\f$, so this
457 * leaves a pi ambiguity. But the surrogate data is already aligned such that
458 * the heavier BH is on the +ve x-axis at t=-1000M. See Sec.VI.A.4 of
459 * arxiv:1812.07865. This resolves the pi ambiguity. This means that after the
460 * realignment, the orbital phase at reference frequency fRef is 0, or Bh1 is
461 * on the +ve x-axis. Note that this is alignment is done using only waveform
462 * quantities, so this doesn't necessarily correspond to the (gauge dependent)
463 * NR Bh trajectories. The polarizations of the waveform are returned in the sky
464 * of this frame at (inclination, pi/2 - phiRef). This agrees with the LAL
465 * convention.
466 *
467 * Only uses data at (2,2) mode frequencies >= fMin. This determines the start
468 * time. The start time, along with the step size deltaT, is used to determine
469 * the output_times. Uses cubic spline interpolation to interpolate from the
470 * surrogate's time array to output_times.
471 *
472 * By default, uses all available modes of the surrogate, that is \f$ \ell \leq
473 * 4, m \geq 0 \f$, and (5,5), but not (4,1) or (4,0). For m>0 modes, the
474 * contribution from the m<0 counterparts is automatically added.
475 *
476 * If a custom ModeArray is given, only those modes are used. Note that it only
477 * looks for m>=0 modes in ModeArray, and will ignore m<0 modes even if present.
478 * The m<0 modes automatically get added.
479 *
480 * This surrogate model is trained on the following range.
481 * q <= 8, |chi_1|, |chi_2| <= 0.8
482 * If you want a guarantee of accuracy you should stick to the above ranges.
483 *
484 * We allow extrapolation to the following ranges, but with a warning:
485 * q = [1, 10.1] and \f$\chi_{1z}, \chi_{2z}\f$ = [-0.81, 0.81]
486 * or
487 * q = [1, 9.1] and \f$\chi_{1z}, \chi_{2z}\f$ = [-0.91, 0.91]
488 * We expect the model to be reasonable when extrapolated to these ranges.
489 *
490 * Beyond the above ranges, we raise an error. If you want to extrapolate
491 * beyond these limits you can specify unlimited_extrapolation = 1 in your
492 * dictParams as follows:
493 * # USE AT YOUR OWN RISK!!
494 * lal.DictInsertUINT4Value(dictParams, "unlimited_extrapolation", 1)
495 */
497 REAL8TimeSeries **hplus, /**<Output: \f$h_+\f$ polarization. */
498 REAL8TimeSeries **hcross, /**<Output: \f$h_{\times}\f$ polarization.*/
499 REAL8 phiRef, /**< azimuthal angle for Ylms */
500 REAL8 inclination, /**< Inclination angle. */
501 REAL8 deltaT, /**< Sampling interval (s). */
502 REAL8 m1, /**< Mass of Bh1 (kg). */
503 REAL8 m2, /**< Mass of Bh2 (kg). */
504 REAL8 distance, /**< Distance of source (m). */
505 REAL8 fMin, /**< Start GW frequency (Hz). */
506 REAL8 fRef, /**< Reference GW frequency (Hz). */
507 REAL8 chi1z, /**< Dimensionless spin of Bh1. */
508 REAL8 chi2z, /**< Dimensionless spin of Bh2. */
509 LALDict* LALparams /**< Dict with extra parameters */
510)
511{
512#ifdef LAL_PTHREAD_LOCK
513 (void) pthread_once(&NRHybSur3dq8_is_initialized, NRHybSur3dq8_Init_LALDATA);
514#else
516#endif
517
518 // Loaded surrogate data
519 const NRHybSurData *NR_hybsur_data = &__lalsim_NRHybSur3dq8_data;
520
521 if (NR_hybsur_data->setup != 1){
522 XLAL_ERROR(XLAL_FAILURE, "Surrogate data is not loaded.");
523 }
524
525 // If ModeArray is not specified, use all available modes
526 LALValue* ModeArray
528 if (ModeArray == NULL) {
529 ModeArray = XLALSimInspiralCreateModeArray();
530 NRHybSur_set_default_modes(ModeArray, NR_hybsur_data);
531 }
532
533 // Make sure we didn't request any unavailable modes, and get number of
534 // modes to include
535 UINT4 num_modes_incl, max_ell;
536 int ret = NRHybSur_check_mode_array(&num_modes_incl, &max_ell, ModeArray,
537 NR_hybsur_data);
538 if (ret != XLAL_SUCCESS) {
539 XLAL_ERROR(XLAL_EFUNC, "Inappropriate ModeArray.");
540 }
541
542 // Check Nyquist frequency is larger than the ringdow frequency
543 // of the (max_ell,max_ell,0) QNM overtone, where max_ell is the
544 // highest ell index included. Raises warning only, not error.
545 ret = NRHybSur_sanity_check_sample_rate(deltaT, m1, m2, chi1z, chi2z,
546 max_ell);
547 if (ret != XLAL_SUCCESS) {
548 XLAL_ERROR(XLAL_EFUNC, "check_sample_rate failed.");
549 }
550
551 // Total mass in seconds
552 const REAL8 Mtot_sec = (m1 + m2)/LAL_MSUN_SI * LAL_MTSUN_SI;
553
554 // mass ratio
555 const REAL8 q = m1/m2;
556
557 // rescale factor for the amplitude based on the distance and total mass
558 const REAL8 a0 = Mtot_sec * LAL_C_SI / distance;
559
560 // Initialize vector to hold phase of the (2,2,) mode
561 gsl_vector *phi_22 = NULL;
562
563 // Initialize data pieces for amplitude of (2,2) mode, and coorbital
564 // frame strain for other modes.
565 EvaluatedDataPieces **evaluated_mode_dps
566 = XLALMalloc(num_modes_incl * sizeof(*evaluated_mode_dps));
567
569
570 // Evaluate the surrogate
571 ret = NRHybSur3dq8_core(&phi_22, evaluated_mode_dps,
572 &epoch, deltaT, fMin, fRef, q, Mtot_sec, chi1z, chi2z,
573 ModeArray, LALparams);
574 if(ret != XLAL_SUCCESS) {
575 XLAL_ERROR(XLAL_EFUNC, "Surrogate evaluation failed.");
576 }
577
578 const UINT4 num_times = phi_22->size;
579
580 // initialize the hp and hc vectors to 0
581 REAL8TimeSeries *hPlusTS =
583 num_times);
584 REAL8TimeSeries *hCrossTS =
586 num_times);
587 memset(hPlusTS->data->data, 0, hPlusTS->data->length*sizeof(REAL8));
588 memset(hCrossTS->data->data, 0, hCrossTS->data->length*sizeof(REAL8));
589
590 // Sum up all required modes
591 COMPLEX16 swsh_coef;
592 COMPLEX16 swsh_coef_negM = 0; // some compilers complain if this is
593 // not initialized here
594 COMPLEX16 hlm;
595 COMPLEX16 hcomplex;
596 int pre_factor;
597 UINT4 ell, m;
598 for (UINT4 mode_idx = 0; mode_idx < num_modes_incl; mode_idx++){
599
600 EvaluatedDataPieces *this_mode_eval_dp = evaluated_mode_dps[mode_idx];
601 ell = this_mode_eval_dp->ell;
602 m = this_mode_eval_dp->m;
603
604 /* See Harald Pfeiffer, T18002260-v1 for frame diagram. The surrogate
605 * frame has its z direction aligned with the orbital angular momentum,
606 * which agrees with the lowercase source frame z in the diagram. The
607 * surrogate x direction, however, points along the line of ascending
608 * nodes (the funny omega with circles on the ends). The uppercase Z,
609 * which is the direction along which we want to evaluate the waveform,
610 * is always in the surrogate frame (y, z) plane. Z is rotated from z
611 * towards the +ive y surrogate axis, so we should always evaluate at
612 * (inclination, pi/2-phiRef). */
613 swsh_coef = XLALSpinWeightedSphericalHarmonic(inclination,
614 0.5 * LAL_PI - phiRef, -2, ell, m);
615
616 if (m > 0) {
617 // Note: Each time a new mode is added to the strain, the value of
618 // swsh_coef_negM persists. When m>0 the old value is overwritten,
619 // while if m=0 it has a meaningless value from the previous
620 // iteration. This is not a bug, but can be confusing, so just
621 // adding this note for anyone that is debugging this code.
622 swsh_coef_negM = XLALSpinWeightedSphericalHarmonic(inclination,
623 0.5 * LAL_PI - phiRef, -2, ell, -m);
624 }
625
626 for (UINT4 j=0; j < num_times; j++) {
627
628 // The (2,2) phase is required for all modes; for the other modes,
629 // we use it to transform from the coorbital frame to the inertial
630 // frame
631 const REAL8 tmp_phi_22 = gsl_vector_get(phi_22, j);
632
633 if (ell == 2 && m ==2){
634 const REAL8 tmp_Amp_22 = gsl_vector_get(
635 this_mode_eval_dp->ampl_eval, j);
636 // (2,2) mode strain in the inertial frame
637 hlm = tmp_Amp_22 * (cos(tmp_phi_22) - I*sin(tmp_phi_22));
638
639 } else {
640
641 // set the real and imaginary parts (in coorbital frame) to
642 // zero by default
643 REAL8 coorb_hlm_re = 0.0;
644 REAL8 coorb_hlm_im = 0.0;
645
646 // for m=0, l=odd, the real part is zero. For all other
647 // cases get the real part
648 if (m != 0 || ell % 2 == 0) {
649 coorb_hlm_re = gsl_vector_get(
650 this_mode_eval_dp->coorb_re_eval, j);
651 }
652
653 // for m=0, l=even, the imaginary part is zero. For all other
654 // cases get the imaginary part
655 if (m != 0 || ell % 2 == 1) {
656 coorb_hlm_im = gsl_vector_get(
657 this_mode_eval_dp->coorb_im_eval, j);
658 }
659
660 // strain in the inertial frame for all non (2,2) modes.
661 // See Eq.(39) of arxiv:1812.07865.for definition of coorbital
662 // frame strain
663 hlm = (coorb_hlm_re + I * coorb_hlm_im)
664 * (cos(m*tmp_phi_22/2.) - I*sin(m*tmp_phi_22/2.));
665 }
666
667 hcomplex = hlm * swsh_coef;
668 if (m > 0) {
669 // For m>0 modes, also add contribution from m<0 modes as
670 // (-1)**l * conj(hlm) * Y_{ell, -m}
671 if (ell % 2 == 0){
672 pre_factor = 1;
673 }
674 else {
675 pre_factor = -1;
676 }
677 hcomplex += pre_factor * conj(hlm) * swsh_coef_negM;
678 }
679
680 // hcomplex = h_+ - i*h_x
681 hPlusTS->data->data[j] += creal(hcomplex) * a0;
682 hCrossTS->data->data[j] -= cimag(hcomplex) * a0;
683 }
684 }
685
686 // Point the output pointers to the relevant time series and return
687 (*hplus) = hPlusTS;
688 (*hcross) = hCrossTS;
689
690 // Cleanup
691 NRHybSur_DestroyEvaluatedDataPieces(phi_22, evaluated_mode_dps,
692 num_modes_incl);
693 if(ModeArray) {
694 XLALDestroyValue(ModeArray);
695 }
696
697
698 return XLAL_SUCCESS;
699}
700
701
702
703/**
704 * Reference: arxiv:1812.07865.
705 *
706 * Returns the spin-weighted spherical harmonic modes of NRHybSur3dq8 waveform
707 * model.
708 *
709 * The reference point is the time at which the (2,2) mode frequency equals
710 * fRef. If fRef=0, sets fRef=fMin. We set the orbital phase to 0 at the
711 * reference point. The orbital phase is obtained as \f$\phi_{22}/2\f$, so this
712 * leaves a pi ambiguity. But the surrogate data is already aligned such that
713 * the heavier BH is on the +ve x-axis at t=-1000M. See Sec.VI.A.4 of
714 * arxiv:1812.07865. This resolves the pi ambiguity. This means that after the
715 * realignment, the orbital phase at reference frequency fRef is 0, or Bh1 is
716 * on the +ve x-axis. Note that this is alignment is done using only waveform
717 * quantities, so this doesn't necessarily correspond to the (gauge dependent)
718 * NR Bh trajectories. The modes are returned in this frame, which agrees with
719 * the LAL convention. When combining the modes to get the polarizations, the
720 * Ylms should be evaluated at (inclination, pi/2 - phiRef), following the LAL
721 * convention.
722 *
723 * Only uses data at (2,2) mode frequencies >= fMin. This determines the start
724 * time. The start time, along with the step size deltaT, is used to determine
725 * the output_times. Uses cubic spline interpolation to interpolate from the
726 * surrogate's time array to output_times.
727 *
728 * By default, returns all available modes of the surrogate, that is \f$ \ell
729 * \leq 4, m \geq 0 \f$, and (5,5), but not (4,1) or (4,0).
730 *
731 * If a custom ModeArray is given, only those modes are used. Note that it only
732 * looks for m>=0 modes in ModeArray, and will ignore m<0 modes even if present.
733 *
734 * This surrogate model is trained on the following range.
735 * q <= 8, |chi_1|, |chi_2| <= 0.8
736 * If you want a guarantee of accuracy you should stick to the above ranges.
737 *
738 * We allow extrapolation to the following ranges, but with a warning:
739 * q = [1, 10.1] and \f$\chi_{1z}, \chi_{2z}\f$ = [-0.81, 0.81]
740 * or
741 * q = [1, 9.1] and \f$\chi_{1z}, \chi_{2z}\f$ = [-0.91, 0.91]
742 * We expect the model to be reasonable when extrapolated to these ranges.
743 *
744 * Beyond the above ranges, we raise an error. If you want to extrapolate
745 * beyond these limits you can specify unlimited_extrapolation = 1 in your
746 * dictParams as follows:
747 * # USE AT YOUR OWN RISK!!
748 * lal.DictInsertUINT4Value(dictParams, "unlimited_extrapolation", 1)
749 */
751 REAL8 deltaT, /**< Sampling interval (s). */
752 REAL8 m1, /**< Mass of Bh1 (kg). */
753 REAL8 m2, /**< Mass of Bh2 (kg). */
754 REAL8 chi1z, /**< Dimensionless spin of Bh1. */
755 REAL8 chi2z, /**< Dimensionless spin of Bh2. */
756 REAL8 fMin, /**< Start GW frequency (Hz). */
757 REAL8 fRef, /**< Reference GW frequency (Hz). */
758 REAL8 distance, /**< Distance of source (m). */
759 LALDict* LALparams /**< Dict with extra parameters */
760) {
761#ifdef LAL_PTHREAD_LOCK
762 (void) pthread_once(&NRHybSur3dq8_is_initialized, NRHybSur3dq8_Init_LALDATA);
763#else
765#endif
766
767 // Loaded surrogate data
768 const NRHybSurData *NR_hybsur_data = &__lalsim_NRHybSur3dq8_data;
769
770 if (NR_hybsur_data->setup != 1){
771 XLAL_ERROR_NULL(XLAL_FAILURE, "Surrogate data is not loaded.");
772 }
773
774 // If ModeArray is not specified, use all available modes
775 LALValue* ModeArray
777 if (ModeArray == NULL) {
778 ModeArray = XLALSimInspiralCreateModeArray();
779 NRHybSur_set_default_modes(ModeArray, NR_hybsur_data);
780 }
781
782 // Make sure we didn't request any unavailable modes, and get number of
783 // modes to include
784 UINT4 num_modes_incl, max_ell;
785 int ret = NRHybSur_check_mode_array(&num_modes_incl, &max_ell, ModeArray,
786 NR_hybsur_data);
787 if (ret != XLAL_SUCCESS) {
788 XLAL_ERROR_NULL(XLAL_EFUNC, "Inappropriate ModeArray.");
789 }
790
791 // Check Nyquist frequency is larger than the ringdow frequency
792 // of the (max_ell,max_ell,0) QNM overtone, where max_ell is the
793 // highest ell index included. Raises warning only, not error.
794 ret = NRHybSur_sanity_check_sample_rate(deltaT, m1, m2, chi1z, chi2z,
795 max_ell);
796 if (ret != XLAL_SUCCESS) {
797 XLAL_ERROR_NULL(XLAL_EFUNC, "check_sample_rate failed.");
798 }
799
800 // Total mass in seconds
801 const REAL8 Mtot_sec = (m1 + m2)/LAL_MSUN_SI * LAL_MTSUN_SI;
802
803 // mass ratio
804 const REAL8 q = m1/m2;
805
806 // rescale factor for the amplitude based on the distance and total mass
807 const REAL8 a0 = Mtot_sec * LAL_C_SI / distance;
808
809 // Initialize vector to hold phase of the (2,2,) mode
810 gsl_vector *phi_22 = NULL;
811
812 // Initialize data pieces for amplitude of (2,2) mode, and coorbital
813 // frame strain for other modes.
814 EvaluatedDataPieces **evaluated_mode_dps
815 = XLALMalloc(num_modes_incl * sizeof(*evaluated_mode_dps));
816
818
819 // Evaluate the surrogate
820 ret = NRHybSur3dq8_core(&phi_22, evaluated_mode_dps,
821 &epoch, deltaT, fMin, fRef, q, Mtot_sec, chi1z, chi2z,
822 ModeArray, LALparams);
823 if(ret != XLAL_SUCCESS) {
824 XLAL_ERROR_NULL(XLAL_EFUNC, "Surrogate evaluation failed.");
825 }
826
827 const UINT4 num_times = phi_22->size;
828
829 // Evaluate all required modes
830 SphHarmTimeSeries *hlms = NULL;
831 COMPLEX16TimeSeries *tmp_mode;
832 UINT4 ell, m;
833 char mode_name[32];
834 for (UINT4 mode_idx = 0; mode_idx < num_modes_incl; mode_idx++){
835
836 EvaluatedDataPieces *this_mode_eval_dp = evaluated_mode_dps[mode_idx];
837 ell = this_mode_eval_dp->ell;
838 m = this_mode_eval_dp->m;
839
840 snprintf(mode_name, sizeof(mode_name), "(%u, %u) mode", ell, m);
841 tmp_mode = XLALCreateCOMPLEX16TimeSeries(mode_name, &epoch, 0.0,
842 deltaT, &lalStrainUnit, num_times);
843
844 for (UINT4 j=0; j < num_times; j++) {
845
846 const REAL8 tmp_phi_22 = gsl_vector_get(phi_22, j);
847
848 if (ell == 2 && m ==2){
849 const REAL8 tmp_Amp_22 = gsl_vector_get(
850 this_mode_eval_dp->ampl_eval, j);
851
852 // (2,2) mode strain in the inertial frame
853 tmp_mode->data->data[j]
854 = a0* tmp_Amp_22 * (cos(tmp_phi_22) - I*sin(tmp_phi_22));
855
856 } else {
857
858 // set the real and imaginary parts (in coorbital frame) to
859 // zero by default
860 REAL8 coorb_hlm_re = 0.0;
861 REAL8 coorb_hlm_im = 0.0;
862
863 // for m=0, l=odd, the real part is zero. For all other
864 // cases get the real part
865 if (m != 0 || ell % 2 == 0) {
866 coorb_hlm_re = gsl_vector_get(
867 this_mode_eval_dp->coorb_re_eval, j);
868 }
869
870 // for m=0, l=even, the imaginary part is zero. For all other
871 // cases get the imaginary part
872 if (m != 0 || ell % 2 == 1) {
873 coorb_hlm_im = gsl_vector_get(
874 this_mode_eval_dp->coorb_im_eval, j);
875 }
876
877 // strain in the inertial frame for all non (2,2) modes.
878 // See Eq.(39) of arxiv:1812.07865.for definition of coorbital
879 // frame strain
880 tmp_mode->data->data[j]
881 = a0 * (coorb_hlm_re + I * coorb_hlm_im)
882 * (cos(m*tmp_phi_22/2.) - I*sin(m*tmp_phi_22/2.));
883 }
884 }
885
886 hlms = XLALSphHarmTimeSeriesAddMode(hlms, tmp_mode, ell, m);
888 }
889
890 // Cleanup
891 NRHybSur_DestroyEvaluatedDataPieces(phi_22, evaluated_mode_dps,
892 num_modes_incl);
893 if(ModeArray) {
894 XLALDestroyValue(ModeArray);
895 }
896
897 return hlms;
898}
struct tagLALH5File LALH5File
int XLALDictContains(const LALDict *dict, const char *key)
UINT4 XLALDictLookupUINT4Value(const LALDict *dict, const char *key)
Auxiliary functions related to HDF5 waveform data files.
SphHarmTimeSeries * XLALSimIMRNRHybSur3dq8Modes(REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 chi1z, REAL8 chi2z, REAL8 fMin, REAL8 fRef, REAL8 distance, LALDict *LALparams)
Reference: arxiv:1812.07865.
int NRHybSur3dq8_fitParams(gsl_vector *fit_params, const REAL8 q, const REAL8 chi1z, const REAL8 chi2z)
Map from mass ratio and spins to surrogate fit parameters.
static bool NRHybSur3dq8_IsSetup(void)
Helper function to check if the NRHybSur3dq8 model has been initialized.
int NRHybSur3dq8_core(gsl_vector **phi_22, EvaluatedDataPieces **evaluated_mode_dps, LIGOTimeGPS *epoch, const REAL8 deltaT, const REAL8 fMin, const REAL8 fRef, REAL8 q, const REAL8 Mtot_sec, REAL8 chi1z, REAL8 chi2z, LALValue *ModeArray, LALDict *LALparams)
This is the core function of the NRHybSur3dq8 model.
static void NRHybSur3dq8_Init_LALDATA(void)
Surrogate initialization.
INT4 XLALSimIMRNRHybSur3dq8Polarizations(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, REAL8 phiRef, REAL8 inclination, REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 distance, REAL8 fMin, REAL8 fRef, REAL8 chi1z, REAL8 chi2z, LALDict *LALparams)
Reference: arxiv:1812.07865.
static NRHybSurData __lalsim_NRHybSur3dq8_data
Global surrogate data.
C code for NRHybSur3dq8 waveform model, an NR-hybrid surrogate model for aligned-spin BBH.
static const char NRHybSur3dq8_DATAFILE[]
LALValue * XLALSimInspiralWaveformParamsLookupModeArray(LALDict *params)
int NRHybSur_check_mode_array(UINT4 *num_modes_incl, UINT4 *max_ell, LALValue *ModeArray, const NRHybSurData *NR_hybsur_data)
Sanity checks on ModeArray.
void NRHybSur_set_default_modes(LALValue *ModeArray, const NRHybSurData *NR_hybsur_data)
Activates all modes of an NRHybSur model.
int NRHybSur_eval_phase_22(gsl_vector **phi_22, gsl_vector **output_times, const REAL8 eta, const gsl_vector *fit_params, const REAL8 omegaM_22_min, const REAL8 deltaTOverM, const REAL8 phiRef, const REAL8 omegaM_22_ref, gsl_vector *dummy_dp, const gsl_matrix *x_train, gsl_vector *dummy_worker, const NRHybSurData *NR_hybsur_data)
Evaluate the phase of the (2,2) mode.
int NRHybSur_eval_mode_data_pieces(EvaluatedDataPieces **this_mode_eval_dp, const UINT4 ell, const UINT4 m, const ModeDataPieces *data_pieces, const gsl_vector *output_times, const gsl_vector *fit_params, gsl_vector *dummy_dp, const gsl_matrix *x_train, gsl_vector *dummy_worker, const NRHybSurData *NR_hybsur_data)
Evaluate waveform data pieces of a single mode.
int NRHybSur_sanity_check_sample_rate(REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 chi1z, REAL8 chi2z, UINT4 max_ell)
Sanity check (warning only, not error) that the sample rate is high enough to capture the ringdown fr...
void NRHybSur_DestroyEvaluatedDataPieces(gsl_vector *phi_22, EvaluatedDataPieces **evaluated_mode_dps, const UINT4 num_modes_incl)
Destroy phi_22 and an EvaluatedDataPieces structure.
void XLALDestroyValue(LALValue *value)
#define XLAL_FILE_RESOLVE_PATH(fname)
void XLALH5FileClose(LALH5File UNUSED *file)
LALH5File * XLALH5FileOpen(const char UNUSED *path, const char UNUSED *mode)
#define LAL_C_SI
#define LAL_MSUN_SI
#define LAL_PI
#define LAL_MTSUN_SI
#define LIGOTIMEGPSZERO
double complex COMPLEX16
double REAL8
uint32_t UINT4
int32_t INT4
void * XLALMalloc(size_t n)
void XLALFree(void *p)
int XLALSimInspiralModeArrayIsModeActive(LALValue *modes, unsigned l, int m)
LALValue * XLALSimInspiralCreateModeArray(void)
Create a LALValue pointer to store the mode array.
SphHarmTimeSeries * XLALSphHarmTimeSeriesAddMode(SphHarmTimeSeries *appended, const COMPLEX16TimeSeries *inmode, UINT4 l, INT4 m)
Prepend a node to a linked list of SphHarmTimeSeries, or create a new head.
static const INT4 m
static const INT4 q
COMPLEX16 XLALSpinWeightedSphericalHarmonic(REAL8 theta, REAL8 phi, int s, int l, int m)
void XLALDestroyCOMPLEX16TimeSeries(COMPLEX16TimeSeries *series)
COMPLEX16TimeSeries * XLALCreateCOMPLEX16TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
const LALUnit lalStrainUnit
#define XLAL_ERROR_VOID(...)
#define XLAL_ERROR_NULL(...)
#define XLAL_ERROR(...)
#define XLAL_CHECK(assertion,...)
#define XLAL_PRINT_WARNING(...)
XLAL_EDATA
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EDOM
XLAL_EIO
XLAL_EDIMS
XLAL_FAILURE
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
path
COMPLEX16Sequence * data
COMPLEX16 * data
NRHybSur evaluated data for a single mode.
gsl_vector * coorb_re_eval
Coorbital frame real part evaluation.
gsl_vector * ampl_eval
Amplitude data piece evaluation.
gsl_vector * coorb_im_eval
Coorbital frame imag part evaluation.
NRHybSur data pieces of a single mode.
UINT4 m
Mode m value.
UINT4 ell
Mode value.
NRHybSur surrogate data for all modes, to be loaded from a h5 file.
ModeDataPieces ** mode_data_pieces
Data pieces of all modes, same order as mode_list.
gsl_vector * domain
Common time samples for the surrogate.
REAL8 params_dim
Dimensions of the model.
UINT4 setup
Indicates if NRHybSur has been initialized.
gsl_matrix * x_train
Training set parameters, needed for GPR fits.
UINT4 num_modes_modeled
Number of modeled modes.
gsl_matrix_long * mode_list
List of modeled modes, first element is always (2,2).
REAL8Sequence * data
REAL8 * data
Structure to carry a collection of spherical harmonic modes in COMPLEX16 time series.
LIGOTimeGPS epoch
Definition: unicorn.c:20
double deltaT
Definition: unicorn.c:24