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
LALSimIMRPhenomHM.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2017 Sebastian Khan, Francesco Pannarale, Lionel London
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#include <lal/LALSimIMR.h>
21#include <lal/SphericalHarmonics.h>
22#include <lal/Date.h>
23#include <lal/Units.h>
24
25#include "LALSimIMRPhenomHM.h"
28#include "LALSimRingdownCW.h"
30
31/*
32 * Phase shift due to leading order complex amplitude
33 * [L.Blancet, arXiv:1310.1528 (Sec. 9.5)]
34 * "Spherical hrmonic modes for numerical relativity"
35 */
36/* List of phase shifts: the index is the azimuthal number m */
37static const double cShift[7] = {0.0,
38 LAL_PI_2 /* i shift */,
39 0.0,
40 -LAL_PI_2 /* -i shift */,
41 LAL_PI /* 1 shift */,
42 LAL_PI_2 /* -1 shift */,
43 0.0};
44
45/**
46 * read in a LALDict.
47 * If ModeArray in LALDict is NULL then create a ModrArray
48 * with the default modes in PhenomHM.
49 * If ModeArray is not NULL then use the modes supplied by user.
50 */
52 LALDict *extraParams)
53{
54
55 /* setup ModeArray */
56 if (extraParams == NULL)
57 extraParams = XLALCreateDict();
58 LALValue *ModeArray = XLALSimInspiralWaveformParamsLookupModeArray(extraParams);
59 if (ModeArray == NULL)
60 { /* Default behaviour */
61 XLAL_PRINT_INFO("Using default modes for PhenomHM.\n");
63 /* Only need to define the positive m modes/
64 * The negative m modes are automatically added.
65 */
72 // XLALSimInspiralModeArrayPrintModes(ModeArray);
73 /* Don't forget to insert ModeArray back into extraParams. */
75 }
76 else
77 {
78 XLAL_PRINT_INFO("Using custom modes for PhenomHM.\n");
79 }
80
81 XLALDestroyValue(ModeArray);
82
83 return extraParams;
84}
85
86/**
87 * Reads in a ModeArray and checks that it is valid.
88 * may only contain the modes in the model
89 * i.e., 22, 21, 33, 32, 44, 43
90 * Only checks upto ell=8 though.
91 */
92static int IMRPhenomHM_check_mode_array(LALValue *ModeArray)
93{
94 // if no 22,21,33,32,44,43 mode and active -> error
95 // these modes are not in the model
96 for (INT4 ell = 2; ell <= 8; ell++)
97 {
98 for (INT4 mm = -ell; mm < ell + 1; mm++)
99 {
100 if (ell == 2 && mm == 2)
101 {
102 continue;
103 }
104 else if (ell == 2 && mm == 1)
105 {
106 continue;
107 }
108 else if (ell == 3 && mm == 3)
109 {
110 continue;
111 }
112 else if (ell == 3 && mm == 2)
113 {
114 continue;
115 }
116 else if (ell == 4 && mm == 4)
117 {
118 continue;
119 }
120 else if (ell == 4 && mm == 3)
121 {
122 continue;
123 }
124
125 if (XLALSimInspiralModeArrayIsModeActive(ModeArray, ell, mm) == 1)
126 {
127 XLAL_ERROR(XLAL_EFUNC, "(%i,%i) mode in ModeArray but model does not include this!\n", ell, mm);
128 }
129 }
130 }
131 return XLAL_SUCCESS;
132}
133
134/**
135 *
136 */
138{
139 XLAL_CHECK(0 != p, XLAL_EFAULT, "p is NULL");
140 XLAL_CHECK(!(number < 0), XLAL_EDOM, "number must be non-negative");
141
142 // consider changing pow(x,1/6.0) to cbrt(x) and sqrt(x) - might be faster
143 p->itself = number;
144 p->sqrt = sqrt(number);
145 p->sixth = cbrt(p->sqrt);
146 p->m_sixth = 1.0 / p->sixth;
147 p->third = p->sixth * p->sixth;
148 p->two_thirds = p->third * p->third;
149 p->four_thirds = number * p->third;
150 p->five_thirds = number * p->two_thirds;
151 p->two = number * number;
152 p->seven_thirds = p->third * p->two;
153 p->eight_thirds = p->two_thirds * p->two;
154 p->m_seven_sixths = p->m_sixth / number;
155 p->m_five_sixths = p->m_seven_sixths * p->third;
156 p->m_sqrt = 1. / p->sqrt;
157
158 return XLAL_SUCCESS;
159}
160
162{
163 XLAL_CHECK(0 != p, XLAL_EFAULT, "p is NULL");
164 XLAL_CHECK(number >= 0, XLAL_EDOM, "number must be non-negative");
165
166 // consider changing pow(x,1/6.0) to cbrt(x) and sqrt(x) - might be faster
167 double sixth = pow(number, 1.0 / 6.0);
168 p->third = sixth * sixth;
169 p->two_thirds = p->third * p->third;
170 p->four_thirds = number * p->third;
171 p->five_thirds = p->four_thirds * p->third;
172 p->two = number * number;
173 p->seven_thirds = p->third * p->two;
174 p->eight_thirds = p->two_thirds * p->two;
175 p->inv = 1. / number;
176 double m_sixth = 1.0 / sixth;
177 p->m_seven_sixths = p->inv * m_sixth;
178 p->m_third = m_sixth * m_sixth;
179 p->m_two_thirds = p->m_third * p->m_third;
180 p->m_five_thirds = p->inv * p->m_two_thirds;
181
182 return XLAL_SUCCESS;
183}
184
185/**
186 * returns the real and imag parts of the complex ringdown frequency
187 * for the (l,m) mode.
188 */
190 REAL8 *fringdown,
191 REAL8 *fdamp,
192 UINT4 ell,
193 INT4 mm,
194 REAL8 finalmass,
195 REAL8 finalspin)
196{
197 const REAL8 inv2Pi = 0.5 / LAL_PI;
198 complex double ZZ;
199 ZZ = SimRingdownCW_CW07102016(SimRingdownCW_KAPPA(finalspin, ell, mm), ell, mm, 0);
200 const REAL8 Mf_RD_tmp = inv2Pi * creal(ZZ); /* GW ringdown frequency, converted from angular frequency */
201 *fringdown = Mf_RD_tmp / finalmass; /* scale by predicted final mass */
202 /* lm mode ringdown damping time (imaginary part of ringdown), geometric units */
203 const REAL8 f_DAMP_tmp = inv2Pi * cimag(ZZ); /* this is the 1./tau in the complex QNM */
204 *fdamp = f_DAMP_tmp / finalmass; /* scale by predicted final mass */
205
206 return XLAL_SUCCESS;
207}
208
209/**
210 * helper function to easily check if the
211 * input frequency sequence is uniformly space
212 * or a user defined set of discrete frequencies.
213 */
215 REAL8Sequence *freqs,
216 REAL8 deltaF)
217{
218 UINT4 freq_is_uniform = 0;
219 if ((freqs->length == 2) && (deltaF > 0.))
220 {
221 freq_is_uniform = 1;
222 }
223 else if ((freqs->length != 2) && (deltaF <= 0.))
224 {
225 freq_is_uniform = 0;
226 }
227
228 return freq_is_uniform;
229}
230
231/**
232 * derive frequency variables for PhenomHM based on input.
233 * used to set the index on arrays where we have non-zero values.
234 */
236 PhenomHMFrequencyBoundsStorage *p, /**< [out] PhenomHMFrequencyBoundsStorage struct */
237 REAL8Sequence *freqs, /**< Input list of GW frequencies [Hz] */
238 REAL8 Mtot, /**< total mass in solar masses */
239 REAL8 deltaF, /**< frequency spacing */
240 REAL8 f_ref_in /**< reference GW frequency */
241)
242{
243 p->deltaF = deltaF;
244 /* determine how to populate frequency sequence */
245 /* if len(freqs_in) == 2 and deltaF > 0. then
246 * f_min = freqs_in[0]
247 * f_max = freqs_in[1]
248 * else if len(freqs_in) != 2 and deltaF <= 0. then
249 * user has given an arbitrary set of frequencies to evaluate the model at.
250 */
251
252 p->freq_is_uniform = IMRPhenomHM_is_freq_uniform(freqs, p->deltaF);
253
254 if (p->freq_is_uniform == 1)
255 { /* This case we use regularly spaced frequencies */
256 p->f_min = freqs->data[0];
257 p->f_max = freqs->data[1];
258
259 /* If p->f_max == 0. Then we default to the ending frequency
260 * for PhenomHM
261 */
262 if (p->f_max == 0.)
263 {
266 }
267 /* we only need to evaluate the phase from
268 * f_min to f_max with a spacing of deltaF
269 */
270 p->npts = PhenomInternal_NextPow2(p->f_max / p->deltaF) + 1;
271 p->ind_min = (size_t)ceil(p->f_min / p->deltaF);
272 p->ind_max = (size_t)ceil(p->f_max / p->deltaF);
273 XLAL_CHECK((p->ind_max <= p->npts) && (p->ind_min <= p->ind_max), XLAL_EDOM, "minimum freq index %zu and maximum freq index %zu do not fulfill 0<=ind_min<=ind_max<=npts=%zu.", p->ind_min, p->ind_max, p->npts);
274 }
275 else if (p->freq_is_uniform == 0)
276 { /* This case we possibly use irregularly spaced frequencies */
277 /* Check that the frequencies are always increasing */
278 for (UINT4 i = 0; i < freqs->length - 1; i++)
279 {
281 freqs->data[i] - freqs->data[i + 1] < 0.,
283 "custom frequencies must be increasing.");
284 }
285
286 XLAL_PRINT_INFO("Using custom frequency input.\n");
287 p->f_min = freqs->data[0];
288 p->f_max = freqs->data[freqs->length - 1]; /* Last element */
289
290 p->npts = freqs->length;
291 p->ind_min = 0;
292 p->ind_max = p->npts;
293 }
294 else
295 { /* Throw an informative error. */
296 XLAL_PRINT_ERROR("Input sequence of frequencies and deltaF is not \
297 compatible.\nSpecify a f_min and f_max by using a REAL8Sequence of length = 2 \
298 along with a deltaF > 0.\
299 \nIf you want to supply an arbitrary list of frequencies to evaluate the with \
300 then supply those frequencies using a REAL8Sequence and also set deltaF <= 0.");
301 }
302
303 /* Fix default behaviour for f_ref */
304 /* If f_ref = 0. then set f_ref = f_min */
305 p->f_ref = f_ref_in;
306 if (p->f_ref == 0.)
307 {
308 p->f_ref = p->f_min;
309 }
310
311 return XLAL_SUCCESS;
312}
313
314/**
315 * Precompute a bunch of PhenomHM related quantities and store them filling in a
316 * PhenomHMStorage variable
317 */
320 const REAL8 m1_SI,
321 const REAL8 m2_SI,
322 const REAL8 chi1x,
323 const REAL8 chi1y,
324 const REAL8 chi1z,
325 const REAL8 chi2x,
326 const REAL8 chi2y,
327 const REAL8 chi2z,
328 REAL8Sequence *freqs,
329 const REAL8 deltaF,
330 const REAL8 f_ref,
331 const REAL8 phiRef)
332{
333 int retcode;
334 XLAL_CHECK(0 != p, XLAL_EFAULT, "p is NULL");
335
336 p->m1 = m1_SI / LAL_MSUN_SI;
337 p->m2 = m2_SI / LAL_MSUN_SI;
338 p->m1_SI = m1_SI;
339 p->m2_SI = m2_SI;
340 p->Mtot = p->m1 + p->m2;
341 p->eta = p->m1 * p->m2 / (p->Mtot * p->Mtot);
342 p->chi1x = chi1x;
343 p->chi1y = chi1y;
344 p->chi1z = chi1z;
345 p->chi2x = chi2x;
346 p->chi2y = chi2y;
347 p->chi2z = chi2z;
348 p->phiRef = phiRef;
349 p->deltaF = deltaF;
350 p->freqs = freqs;
351
352 if (p->eta > 0.25)
353 PhenomInternal_nudge(&(p->eta), 0.25, 1e-6);
354 if (p->eta > 0.25 || p->eta < 0.0)
355 XLAL_ERROR(XLAL_EDOM, "Unphysical eta. Must be between 0. and 0.25\n");
356 if (p->eta < MAX_ALLOWED_ETA)
357 XLAL_PRINT_WARNING("Warning: The model is not calibrated for mass-ratios above 20\n");
358
359 retcode = 0;
361 &(p->m1),
362 &(p->m2),
363 &(p->chi1x),
364 &(p->chi1y),
365 &(p->chi1z),
366 &(p->chi2x),
367 &(p->chi2y),
368 &(p->chi2z));
370 XLAL_SUCCESS == retcode,
372 "PhenomInternal_AlignedSpinEnforcePrimaryIsm1 failed");
373
374
375 p->chip = XLALSimPhenomUtilsChiP(p->m1, p->m2, p->chi1x, p->chi1y, p->chi2x, p->chi2y);
376
377 /* sanity checks on frequencies */
379 retcode = 0;
381 &pHMFS,
382 p->freqs,
383 p->Mtot,
384 p->deltaF,
385 f_ref);
387 XLAL_SUCCESS == retcode,
389 "init_IMRPhenomHMGet_FrequencyBounds_storage failed");
390
391 /* redundent storage */
392 p->f_min = pHMFS.f_min;
393 p->f_max = pHMFS.f_max;
394 p->f_ref = pHMFS.f_ref;
395 p->freq_is_uniform = pHMFS.freq_is_uniform;
396 p->npts = pHMFS.npts;
397 p->ind_min = pHMFS.ind_min;
398 p->ind_max = pHMFS.ind_max;
399
400 p->Mf_ref = XLALSimPhenomUtilsHztoMf(p->f_ref, p->Mtot);
401
402 p->finmass = XLALSimPhenomUtilsIMRPhenomDFinalMass(p->m1, p->m2, p->chi1z, p->chi2z);
403 p->finspin = XLALSimPhenomUtilsPhenomPv2FinalSpin(p->m1, p->m2, p->chi1z, p->chi2z, p->chip);
404
405 if (p->finspin > 1.0)
406 XLAL_ERROR(XLAL_EDOM, "PhenomD fring function: final spin > 1.0 not supported\n");
407
408 /* populate the ringdown frequency array */
409 /* If you want to model a new mode then you have to add it here. */
410 /* (l,m) = (2,2) */
412 &p->PhenomHMfring[2][2],
413 &p->PhenomHMfdamp[2][2],
414 2, 2,
415 p->finmass, p->finspin);
416
417 /* (l,m) = (2,1) */
419 &p->PhenomHMfring[2][1],
420 &p->PhenomHMfdamp[2][1],
421 2, 1,
422 p->finmass, p->finspin);
423
424 /* (l,m) = (3,3) */
426 &p->PhenomHMfring[3][3],
427 &p->PhenomHMfdamp[3][3],
428 3, 3,
429 p->finmass, p->finspin);
430
431 /* (l,m) = (3,2) */
433 &p->PhenomHMfring[3][2],
434 &p->PhenomHMfdamp[3][2],
435 3, 2,
436 p->finmass, p->finspin);
437
438 /* (l,m) = (4,4) */
440 &p->PhenomHMfring[4][4],
441 &p->PhenomHMfdamp[4][4],
442 4, 4,
443 p->finmass, p->finspin);
444
445 /* (l,m) = (4,3) */
447 &p->PhenomHMfring[4][3],
448 &p->PhenomHMfdamp[4][3],
449 4, 3,
450 p->finmass, p->finspin);
451
452 p->Mf_RD_22 = p->PhenomHMfring[2][2];
453 p->Mf_DM_22 = p->PhenomHMfdamp[2][2];
454
455 /* (l,m) = (2,2) */
456 int ell, mm;
457 ell = 2;
458 mm = 2;
459 p->Rholm[ell][mm] = p->Mf_RD_22 / p->PhenomHMfring[ell][mm];
460 p->Taulm[ell][mm] = p->PhenomHMfdamp[ell][mm] / p->Mf_DM_22;
461 /* (l,m) = (2,1) */
462 ell = 2;
463 mm = 1;
464 p->Rholm[ell][mm] = p->Mf_RD_22 / p->PhenomHMfring[ell][mm];
465 p->Taulm[ell][mm] = p->PhenomHMfdamp[ell][mm] / p->Mf_DM_22;
466 /* (l,m) = (3,3) */
467 ell = 3;
468 mm = 3;
469 p->Rholm[ell][mm] = p->Mf_RD_22 / p->PhenomHMfring[ell][mm];
470 p->Taulm[ell][mm] = p->PhenomHMfdamp[ell][mm] / p->Mf_DM_22;
471 /* (l,m) = (3,2) */
472 ell = 3;
473 mm = 2;
474 p->Rholm[ell][mm] = p->Mf_RD_22 / p->PhenomHMfring[ell][mm];
475 p->Taulm[ell][mm] = p->PhenomHMfdamp[ell][mm] / p->Mf_DM_22;
476 /* (l,m) = (4,4) */
477 ell = 4;
478 mm = 4;
479 p->Rholm[ell][mm] = p->Mf_RD_22 / p->PhenomHMfring[ell][mm];
480 p->Taulm[ell][mm] = p->PhenomHMfdamp[ell][mm] / p->Mf_DM_22;
481 /* (l,m) = (4,3) */
482 ell = 4;
483 mm = 3;
484 p->Rholm[ell][mm] = p->Mf_RD_22 / p->PhenomHMfring[ell][mm];
485 p->Taulm[ell][mm] = p->PhenomHMfdamp[ell][mm] / p->Mf_DM_22;
486
487 return XLAL_SUCCESS;
488};
489
490/**
491 * domain mapping function - ringdown
492 */
494 REAL8 Mf,
495 REAL8 Mf_RD_22,
496 REAL8 Mf_RD_lm,
497 const INT4 AmpFlag,
498 const INT4 ell,
499 const INT4 mm,
500 PhenomHMStorage *pHM)
501{
502 double ans = 0.0;
503 if (AmpFlag == AmpFlagTrue)
504 {
505 /* For amplitude */
506 ans = Mf - Mf_RD_lm + Mf_RD_22; /*Used for the Amplitude as an approx fix for post merger powerlaw slope */
507 }
508 else
509 {
510 /* For phase */
511 REAL8 Rholm = pHM->Rholm[ell][mm];
512 ans = Rholm * Mf; /* Used for the Phase */
513 }
514
515 return ans;
516}
517
518/**
519 * mathematica function Ti
520 * domain mapping function - inspiral
521 */
522double IMRPhenomHMTi(REAL8 Mf, const INT4 mm)
523{
524 return 2.0 * Mf / mm;
525}
526
527/**
528 * helper function for IMRPhenomHMFreqDomainMap
529 */
531 double *Am,
532 double *Bm,
533 const INT4 mm,
534 REAL8 fi,
535 REAL8 fr,
536 REAL8 Mf_RD_22,
537 REAL8 Mf_RD_lm,
538 const INT4 AmpFlag,
539 const INT4 ell,
540 PhenomHMStorage *pHM)
541{
542 REAL8 Trd = IMRPhenomHMTrd(fr, Mf_RD_22, Mf_RD_lm, AmpFlag, ell, mm, pHM);
543 REAL8 Ti = IMRPhenomHMTi(fi, mm);
544
545 //Am = ( Trd[fr]-Ti[fi] )/( fr - fi );
546 *Am = (Trd - Ti) / (fr - fi);
547
548 //Bm = Ti[fi] - fi*Am;
549 *Bm = Ti - fi * (*Am);
550
551 return XLAL_SUCCESS;
552}
553
554/**
555 * helper function for IMRPhenomHMFreqDomainMap
556 */
558 REAL8 *a,
559 REAL8 *b,
560 REAL8 flm,
561 REAL8 fi,
562 REAL8 fr,
563 REAL8 Ai,
564 REAL8 Bi,
565 REAL8 Am,
566 REAL8 Bm,
567 REAL8 Ar,
568 REAL8 Br)
569{
570 // Define function to output map params used depending on
571 if (flm > fi)
572 {
573 if (flm > fr)
574 {
575 *a = Ar;
576 *b = Br;
577 }
578 else
579 {
580 *a = Am;
581 *b = Bm;
582 }
583 }
584 else
585 {
586 *a = Ai;
587 *b = Bi;
588 };
589 return XLAL_SUCCESS;
590}
591
592/**
593 * helper function for IMRPhenomHMFreqDomainMap
594 */
596 REAL8 *a, /**< [Out] */
597 REAL8 *b, /**< [Out] */
598 REAL8 *fi, /**< [Out] */
599 REAL8 *fr, /**< [Out] */
600 REAL8 *f1, /**< [Out] */
601 const REAL8 flm, /**< input waveform frequency */
602 const INT4 ell, /**< spherical harmonics ell mode */
603 const INT4 mm, /**< spherical harmonics m mode */
604 PhenomHMStorage *pHM, /**< Stores quantities in order to calculate them only once */
605 const int AmpFlag /**< is ==1 then computes for amplitude, if ==0 then computes for phase */
606)
607{
608
609 /*check output points are NULL*/
610 XLAL_CHECK(a != NULL, XLAL_EFAULT);
611 XLAL_CHECK(b != NULL, XLAL_EFAULT);
612 XLAL_CHECK(fi != NULL, XLAL_EFAULT);
613 XLAL_CHECK(fr != NULL, XLAL_EFAULT);
614 XLAL_CHECK(f1 != NULL, XLAL_EFAULT);
615
616 /* Account for different f1 definition between PhenomD Amplitude and Phase derivative models */
617 REAL8 Mf_1_22 = 0.; /* initalise variable */
618 if (AmpFlag == AmpFlagTrue)
619 {
620 /* For amplitude */
621 Mf_1_22 = AMP_fJoin_INS; /* inspiral joining frequency from PhenomD [amplitude model], for the 22 mode */
622 }
623 else
624 {
625 /* For phase */
626 Mf_1_22 = PHI_fJoin_INS; /* inspiral joining frequency from PhenomD [phase model], for the 22 mode */
627 }
628
629 *f1 = Mf_1_22;
630
631 REAL8 Mf_RD_22 = pHM->Mf_RD_22;
632 REAL8 Mf_RD_lm = pHM->PhenomHMfring[ell][mm];
633
634 // Define a ratio of QNM frequencies to be used for scaling various quantities
635 REAL8 Rholm = pHM->Rholm[ell][mm];
636
637 // Given experiments with the l!=m modes, it appears that the QNM scaling rather than the PN scaling may be optimal for mapping f1
638 REAL8 Mf_1_lm = Mf_1_22 / Rholm;
639
640 /* Define transition frequencies */
641 *fi = Mf_1_lm;
642 *fr = Mf_RD_lm;
643
644 /*Define the slope and intercepts of the linear transformation used*/
645 REAL8 Ai = 2.0 / mm;
646 REAL8 Bi = 0.0;
647 REAL8 Am;
648 REAL8 Bm;
649 IMRPhenomHMSlopeAmAndBm(&Am, &Bm, mm, *fi, *fr, Mf_RD_22, Mf_RD_lm, AmpFlag, ell, pHM);
650
651 REAL8 Ar = 1.0;
652 REAL8 Br = 0.0;
653 if (AmpFlag == AmpFlagTrue)
654 {
655 /* For amplitude */
656 Br = -Mf_RD_lm + Mf_RD_22;
657 }
658 else
659 {
660 /* For phase */
661 Ar = Rholm;
662 }
663
664 /* Define function to output map params used depending on */
665 int ret = IMRPhenomHMMapParams(a, b, flm, *fi, *fr, Ai, Bi, Am, Bm, Ar, Br);
666 if (ret != XLAL_SUCCESS)
667 {
668 XLALPrintError("XLAL Error - IMRPhenomHMMapParams failed in IMRPhenomHMFreqDomainMapParams (1)\n");
670 }
671
672 return XLAL_SUCCESS;
673}
674
675/**
676 * IMRPhenomHMFreqDomainMap
677 * Input waveform frequency in Geometric units (Mflm)
678 * and computes what frequency this corresponds
679 * to scaled to the 22 mode.
680 */
682 REAL8 Mflm,
683 const INT4 ell,
684 const INT4 mm,
685 PhenomHMStorage *pHM,
686 const int AmpFlag)
687{
688
689 /* Mflm here has the same meaning as Mf_wf in XLALSimIMRPhenomHMFreqDomainMapHM (old deleted function). */
690 REAL8 a = 0.;
691 REAL8 b = 0.;
692 /* Following variables not used in this funciton but are returned in IMRPhenomHMFreqDomainMapParams */
693 REAL8 fi = 0.;
694 REAL8 fr = 0.;
695 REAL8 f1 = 0.;
696 int ret = IMRPhenomHMFreqDomainMapParams(&a, &b, &fi, &fr, &f1, Mflm, ell, mm, pHM, AmpFlag);
697 if (ret != XLAL_SUCCESS)
698 {
699 XLALPrintError("XLAL Error - IMRPhenomHMFreqDomainMapParams failed in IMRPhenomHMFreqDomainMap\n");
701 }
702 REAL8 Mf22 = a * Mflm + b;
703 return Mf22;
704}
705
707 HMPhasePreComp *q, /**< [out] HMPhasePreComp struct */
708 const INT4 ell, /**< ell spherical harmonic number */
709 const INT4 mm, /**< m spherical harmonic number */
710 PhenomHMStorage *pHM, /**< PhenomHMStorage struct */
711 UNUSED LALDict *extraParams /**< LALDict strcut */
712)
713{
714 REAL8 ai = 0.0;
715 REAL8 bi = 0.0;
716 REAL8 am = 0.0;
717 REAL8 bm = 0.0;
718 REAL8 ar = 0.0;
719 REAL8 br = 0.0;
720 REAL8 fi = 0.0;
721 REAL8 f1 = 0.0;
722 REAL8 fr = 0.0;
723
724 const INT4 AmpFlag = 0;
725
726 /* NOTE: As long as Mfshit isn't >= fr then the value of the shift is arbitrary. */
727 const REAL8 Mfshift = 0.0001;
728
729 int ret = IMRPhenomHMFreqDomainMapParams(&ai, &bi, &fi, &fr, &f1, Mfshift, ell, mm, pHM, AmpFlag);
730 if (ret != XLAL_SUCCESS)
731 {
732 XLALPrintError("XLAL Error - IMRPhenomHMFreqDomainMapParams failed in IMRPhenomHMFreqDomainMapParams - inspiral\n");
734 }
735 q->ai = ai;
736 q->bi = bi;
737
738 ret = IMRPhenomHMFreqDomainMapParams(&am, &bm, &fi, &fr, &f1, fi + Mfshift, ell, mm, pHM, AmpFlag);
739 if (ret != XLAL_SUCCESS)
740 {
741 XLALPrintError("XLAL Error - IMRPhenomHMFreqDomainMapParams failed in IMRPhenomHMFreqDomainMapParams - intermediate\n");
743 }
744 q->am = am;
745 q->bm = bm;
746
747 ret = IMRPhenomHMFreqDomainMapParams(&ar, &br, &fi, &fr, &f1, fr + Mfshift, ell, mm, pHM, AmpFlag);
748 if (ret != XLAL_SUCCESS)
749 {
750 XLALPrintError("XLAL Error - IMRPhenomHMFreqDomainMapParams failed in IMRPhenomHMFreqDomainMapParams - merger-ringdown\n");
752 }
753
754 q->ar = ar;
755 q->br = br;
756
757 q->fi = fi;
758 q->fr = fr;
759
760 double Rholm = pHM->Rholm[ell][mm];
761 double Taulm = pHM->Taulm[ell][mm];
762
763 PhenDAmpAndPhasePreComp pDPreComp;
765 &pDPreComp,
766 pHM->m1,
767 pHM->m2,
768 pHM->chi1x,
769 pHM->chi1y,
770 pHM->chi1z,
771 pHM->chi2x,
772 pHM->chi2y,
773 pHM->chi2z,
774 Rholm,
775 Taulm,
776 extraParams);
777 if (ret != XLAL_SUCCESS)
778 {
779 XLALPrintError("XLAL Error - IMRPhenomDSetupAmpAndPhaseCoefficients failed\n");
781 }
782
783 REAL8 PhDBMf = am * fi + bm;
784 q->PhDBconst = IMRPhenomDPhase_OneFrequency(PhDBMf, pDPreComp, Rholm, Taulm) / am;
785
786 REAL8 PhDCMf = ar * fr + br;
787 q->PhDCconst = IMRPhenomDPhase_OneFrequency(PhDCMf, pDPreComp, Rholm, Taulm) / ar;
788
789 REAL8 PhDBAMf = ai * fi + bi;
790 q->PhDBAterm = IMRPhenomDPhase_OneFrequency(PhDBAMf, pDPreComp, Rholm, Taulm) / ai;
791 return XLAL_SUCCESS;
792}
793
794/**
795 * Define function for FD PN amplitudes
796 */
798 REAL8 fM,
799 INT4 l,
800 INT4 m,
801 REAL8 M1,
802 REAL8 M2,
803 REAL8 X1z,
804 REAL8 X2z)
805{
806
807 // LLondon 2017
808
809 // Define effective intinsic parameters
810 COMPLEX16 Hlm = 0;
811 REAL8 M_INPUT = M1 + M2;
812 M1 = M1 / (M_INPUT);
813 M2 = M2 / (M_INPUT);
814 REAL8 M = M1 + M2;
815 REAL8 eta = M1 * M2 / (M * M);
816 REAL8 delta = sqrt(1.0 - 4 * eta);
817 REAL8 Xs = 0.5 * (X1z + X2z);
818 REAL8 Xa = 0.5 * (X1z - X2z);
819 COMPLEX16 ans = 0;
820
821 // Define PN parameter and realed powers
822 REAL8 v = pow(M * 2.0 * LAL_PI * fM / m, 1.0 / 3.0);
823 REAL8 v2 = v * v;
824 REAL8 v3 = v * v2;
825
826 // Define Leading Order Ampitude for each supported multipole
827 if (l == 2 && m == 2)
828 {
829 // (l,m) = (2,2)
830 // THIS IS LEADING ORDER
831 Hlm = 1.0;
832 }
833 else if (l == 2 && m == 1)
834 {
835 // (l,m) = (2,1)
836 // SPIN TERMS ADDED
837
838 // UP TO 4PN
839 REAL8 v4 = v * v3;
840 Hlm = (sqrt(2.0) / 3.0) * \
841 ( \
842 v * delta - v2 * 1.5 * (Xa + delta * Xs) + \
843 v3 * delta * ((335.0 / 672.0) + (eta * 117.0 / 56.0)
844 ) \
845 + \
846 v4 * \
847 ( \
848 Xa * (3427.0 / 1344 - eta * 2101.0 / 336) + \
849 delta * Xs * (3427.0 / 1344 - eta * 965 / 336) + \
850 delta * (-I * 0.5 - LAL_PI - 2 * I * 0.69314718056) \
851 )
852 );
853 }
854 else if (l == 3 && m == 3)
855 {
856 // (l,m) = (3,3)
857 // THIS IS LEADING ORDER
858 Hlm = 0.75 * sqrt(5.0 / 7.0) * (v * delta);
859 }
860 else if (l == 3 && m == 2)
861 {
862 // (l,m) = (3,2)
863 // NO SPIN TERMS to avoid roots
864 Hlm = (1.0 / 3.0) * sqrt(5.0 / 7.0) * (v2 * (1.0 - 3.0 * eta));
865 }
866 else if (l == 4 && m == 4)
867 {
868 // (l,m) = (4,4)
869 // THIS IS LEADING ORDER
870 Hlm = (4.0 / 9.0) * sqrt(10.0 / 7.0) * v2 * (1.0 - 3.0 * eta);
871 }
872 else if (l == 4 && m == 3)
873 {
874 // (l,m) = (4,3)
875 // NO SPIN TERMS TO ADD AT DESIRED ORDER
876 Hlm = 0.75 * sqrt(3.0 / 35.0) * v3 * delta * (1.0 - 2.0 * eta);
877 }
878 else
879 {
880 XLALPrintError("XLAL Error - requested ell = %i and m = %i mode not available, check documentation for available modes\n", l, m);
882 }
883 // Compute the final PN Amplitude at Leading Order in fM
884 ans = M * M * LAL_PI * sqrt(eta * 2.0 / 3) * pow(v, -3.5) * cabs(Hlm);
885
886 return ans;
887}
888
889/**
890 * @addtogroup LALSimIMRPhenom_c
891 * @{
892 *
893 * @name Routines for IMR Phenomenological Model "HM"
894 * @{
895 *
896 * @author Sebastian Khan, Francesco Pannarale, Lionel London
897 *
898 * @brief C code for IMRPhenomHM phenomenological waveform model.
899 *
900 * Inspiral-merger and ringdown phenomenological, frequecny domain
901 * waveform model for binary black holes systems.
902 * Models not only the dominant (l,|m|) = (2,2) modes
903 * but also some of the sub-domant modes too.
904 * Model described in PhysRevLett.120.161102/1708.00404.
905 * The model is based on IMRPhenomD (\cite Husa:2015iqa, \cite Khan:2015jqa)
906 *
907 * @note The higher mode information was not calibrated to Numerical Relativity
908 * simulation therefore the calibration range is inherited from PhenomD.
909 *
910 * @attention The model is usable outside this parameter range,
911 * and in tests to date gives sensible physical results,
912 * but conclusive statements on the physical fidelity of
913 * the model for these parameters await comparisons against further
914 * numerical-relativity simulations. For more information, see the review wiki
915 * under https://git.ligo.org/waveforms/reviews/phenomhm/wikis/home
916 * Also a technical document in the DCC https://dcc.ligo.org/LIGO-T1800295
917 */
918
919/**
920 * Returns h+ and hx in the frequency domain.
921 *
922 * This function can be called in the usual sense
923 * where you supply a f_min, f_max and deltaF.
924 * This is the case when deltaF > 0.
925 * If f_max = 0. then the default ending frequnecy is used.
926 * or you can also supply a custom set of discrete
927 * frequency points with which to evaluate the waveform.
928 * To do this you must call this function with
929 * deltaF <= 0.
930 *
931 */
933 UNUSED COMPLEX16FrequencySeries **hptilde, /**< [out] Frequency-domain waveform h+ */
934 UNUSED COMPLEX16FrequencySeries **hctilde, /**< [out] Frequency-domain waveform hx */
935 UNUSED REAL8Sequence *freqs, /**< Frequency points at which to evaluate the waveform (Hz) */
936 UNUSED REAL8 m1_SI, /**< mass of companion 1 (kg) */
937 UNUSED REAL8 m2_SI, /**< mass of companion 2 (kg) */
938 UNUSED REAL8 chi1z, /**< z-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
939 UNUSED REAL8 chi2z, /**< z-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
940 UNUSED const REAL8 distance, /**< distance of source (m) */
941 UNUSED const REAL8 inclination, /**< inclination of source (rad) */
942 UNUSED const REAL8 phiRef, /**< reference orbital phase (rad) */
943 UNUSED const REAL8 deltaF, /**< Sampling frequency (Hz). To use arbitrary frequency points set deltaF <= 0. */
944 UNUSED REAL8 f_ref, /**< Reference frequency */
945 UNUSED LALDict *extraParams /**<linked list containing the extra testing GR parameters */
946)
947{
948 /* define and init return code for this function */
949 int retcode;
950
951 /* sanity checks on input parameters: check pointers, etc. */
952 /* NOTE: a lot of checks are done in the function
953 * XLALSimIMRPhenomHMGethlmModes because that can also be used
954 * as a standalone function. It gets called through IMRPhenomHMCore
955 * so to avoid doubling up on checks alot of the checks are done in
956 * XLALSimIMRPhenomHMGethlmModes.
957 */
958 XLAL_CHECK(NULL != hptilde, XLAL_EFAULT);
959 XLAL_CHECK(NULL != hctilde, XLAL_EFAULT);
960 XLAL_CHECK(*hptilde == NULL, XLAL_EFAULT);
961 XLAL_CHECK(*hctilde == NULL, XLAL_EFAULT);
962 XLAL_CHECK(distance > 0, XLAL_EDOM, "distance must be positive.\n");
963
964 /* main: evaluate model at given frequencies */
965 retcode = 0;
966 retcode = IMRPhenomHMCore(
967 hptilde,
968 hctilde,
969 freqs,
970 m1_SI,
971 m2_SI,
972 chi1z,
973 chi2z,
974 distance,
975 inclination,
976 phiRef,
977 deltaF,
978 f_ref,
979 extraParams);
980 XLAL_CHECK(retcode == XLAL_SUCCESS,
981 XLAL_EFUNC, "IMRPhenomHMCore failed in XLALSimIMRPhenomHM.");
982
983 /* cleanup */
984 /* XLALDestroy and XLALFree any pointers. */
985
986 return XLAL_SUCCESS;
987}
988
989/** @} */
990/** @} */
991
992/**
993 * internal function that returns h+ and hx.
994 * Inside this function the my bulk of the work is done
995 * like the loop over frequencies.
996 */
998 UNUSED COMPLEX16FrequencySeries **hptilde, /**< [out] Frequency domain h+ GW strain */
999 UNUSED COMPLEX16FrequencySeries **hctilde, /**< [out] Frequency domain hx GW strain */
1000 REAL8Sequence *freqs, /**< GW frequecny list [Hz] */
1001 REAL8 m1_SI, /**< primary mass [kg] */
1002 REAL8 m2_SI, /**< secondary mass [kg] */
1003 REAL8 chi1z, /**< aligned spin of primary */
1004 REAL8 chi2z, /**< aligned spin of secondary */
1005 const REAL8 distance, /**< distance [m] */
1006 const REAL8 inclination, /**< inclination angle */
1007 const REAL8 phiRef, /**< orbital phase at f_ref */
1008 const REAL8 deltaF, /**< frequency spacing */
1009 REAL8 f_ref, /**< reference GW frequency */
1010 LALDict *extraParams /**< LALDict struct */
1011)
1012{
1013 int retcode;
1014
1015 /* Use an auxiliar laldict to not overwrite the input argument */
1016 LALDict *extraParams_aux;
1017
1018 /* setup mode array */
1019 if (extraParams == NULL)
1020 {
1021 extraParams_aux = XLALCreateDict();
1022 }
1023 else{
1024 extraParams_aux = XLALDictDuplicate(extraParams);
1025 }
1026 extraParams_aux = IMRPhenomHM_setup_mode_array(extraParams_aux);
1027
1028 /* evaluate all hlm modes */
1030 *hlms = NULL;
1031 retcode = 0;
1033 hlms,
1034 freqs,
1035 m1_SI,
1036 m2_SI,
1037 0.,
1038 0.,
1039 chi1z,
1040 0.,
1041 0.,
1042 chi2z,
1043 phiRef,
1044 deltaF,
1045 f_ref,
1046 extraParams_aux);
1047 XLAL_CHECK(XLAL_SUCCESS == retcode,
1048 XLAL_EFUNC, "XLALSimIMRPhenomHMGethlmModes failed");
1049
1050 /* compute the frequency bounds */
1051 const REAL8 Mtot = (m1_SI + m2_SI) / LAL_MSUN_SI;
1054 retcode = 0;
1056 pHMFS,
1057 freqs,
1058 Mtot,
1059 deltaF,
1060 f_ref);
1061 XLAL_CHECK(XLAL_SUCCESS == retcode,
1062 XLAL_EFUNC, "init_IMRPhenomHMGet_FrequencyBounds_storage failed");
1063
1064 /* now we have generated all hlm modes we need to
1065 * multiply them with the Ylm's and sum them.
1066 */
1067
1068 LIGOTimeGPS tC = LIGOTIMEGPSZERO; // = {0, 0}
1069 if (pHMFS->freq_is_uniform == 1)
1070 { /* 1. uniformly spaced */
1071 XLAL_PRINT_INFO("freq_is_uniform = True\n");
1072 /* coalesce at t=0 */
1073 /* Shift by overall length in time */
1074 XLAL_CHECK(
1075 XLALGPSAdd(&tC, -1. / deltaF),
1076 XLAL_EFUNC,
1077 "Failed to shift coalescence time to t=0,\
1078tried to apply shift of -1.0/deltaF with deltaF=%g.",
1079 deltaF);
1080 } /* else if 2. i.e. not uniformly spaced then we don't shift. */
1081
1082 /* Allocate hptilde and hctilde */
1083 *hptilde = XLALCreateCOMPLEX16FrequencySeries("hptilde: FD waveform", &tC, 0.0, deltaF, &lalStrainUnit, pHMFS->npts);
1084 if (!(hptilde))
1086 memset((*hptilde)->data->data, 0, pHMFS->npts * sizeof(COMPLEX16));
1087 XLALUnitMultiply(&(*hptilde)->sampleUnits, &(*hptilde)->sampleUnits, &lalSecondUnit);
1088
1089 *hctilde = XLALCreateCOMPLEX16FrequencySeries("hctilde: FD waveform", &tC, 0.0, deltaF, &lalStrainUnit, pHMFS->npts);
1090 if (!(hctilde))
1092 memset((*hctilde)->data->data, 0, pHMFS->npts * sizeof(COMPLEX16));
1093 XLALUnitMultiply(&(*hctilde)->sampleUnits, &(*hctilde)->sampleUnits, &lalSecondUnit);
1094
1095 /* Adding the modes to form hplus, hcross
1096 * - use of a function that copies XLALSimAddMode but for Fourier domain structures */
1097 INT4 sym; /* sym will decide whether to add the -m mode (when equatorial symmetry is present) */
1098
1099 /* setup ModeArray */
1100 LALValue *ModeArray = XLALSimInspiralWaveformParamsLookupModeArray(extraParams_aux);
1101
1102 /* loop over modes */
1103 /* at this point ModeArray should contain the list of modes
1104 * and therefore if NULL then something is wrong and abort.
1105 */
1106 if (ModeArray == NULL)
1107 {
1108 XLAL_ERROR(XLAL_EDOM, "ModeArray is NULL when it shouldn't be. Aborting.\n");
1109 }
1110 for (UINT4 ell = 2; ell < L_MAX_PLUS_1; ell++)
1111 {
1112 for (INT4 mm = 1; mm < (INT4)ell + 1; mm++)
1113 { /* loop over only positive m is intentional. negative m added automatically */
1114 /* first check if (l,m) mode is 'activated' in the ModeArray */
1115 /* if activated then generate the mode, else skip this mode. */
1116 if (XLALSimInspiralModeArrayIsModeActive(ModeArray, ell, mm) != 1)
1117 { /* skip mode */
1118 continue;
1119 } /* else: generate mode */
1120
1122 if (!(hlm))
1124
1125 /* We test for hypothetical m=0 modes */
1126 if (mm == 0)
1127 {
1128 sym = 0;
1129 }
1130 else
1131 {
1132 sym = 1;
1133 }
1134 PhenomInternal_IMRPhenomHMFDAddMode(*hptilde, *hctilde, hlm, inclination, 0., ell, mm, sym); /* The phase \Phi is set to 0 - assumes phiRef is defined as half the phase of the 22 mode h22 */
1135 }
1136 }
1137
1139 XLALFree(hlms);
1140
1141 /* Compute the amplitude pre-factor */
1142 const REAL8 amp0 = XLALSimPhenomUtilsFDamp0(Mtot, distance);
1143 #pragma omp parallel for
1144 for (size_t i = pHMFS->ind_min; i < pHMFS->ind_max; i++)
1145 {
1146 ((*hptilde)->data->data)[i] = ((*hptilde)->data->data)[i] * amp0;
1147 ((*hctilde)->data->data)[i] = ((*hctilde)->data->data)[i] * amp0;
1148 }
1149
1150 /* cleanup */
1151 XLALDestroyValue(ModeArray);
1152 LALFree(pHMFS);
1153
1154 XLALDestroyDict(extraParams_aux);
1155
1156 return XLAL_SUCCESS;
1157}
1158
1159/**
1160 * XLAL function that returns
1161 * a SphHarmFrequencySeries object
1162 * containing all the hlm modes
1163 * requested.
1164 * These have the correct relative phases between modes.
1165 * Note this has a similar interface to XLALSimIMRPhenomHM
1166 * because it is designed so it can be used independently.
1167 */
1169 UNUSED SphHarmFrequencySeries **hlms, /**< [out] SphHarmFrequencySeries struct containing hlm modes */
1170 UNUSED REAL8Sequence *freqs, /**< frequency sequency in Hz */
1171 UNUSED REAL8 m1_SI, /**< primary mass [kg] */
1172 UNUSED REAL8 m2_SI, /**< secondary mass [kg] */
1173 UNUSED REAL8 chi1x, /**< x-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
1174 UNUSED REAL8 chi1y, /**< y-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
1175 UNUSED REAL8 chi1z, /**< z-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
1176 UNUSED REAL8 chi2x, /**< x-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
1177 UNUSED REAL8 chi2y, /**< y-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
1178 UNUSED REAL8 chi2z, /**< z-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
1179 UNUSED const REAL8 phiRef, /**< orbital phase at f_ref */
1180 UNUSED const REAL8 deltaF, /**< frequency spacing */
1181 UNUSED REAL8 f_ref, /**< reference GW frequency */
1182 UNUSED LALDict *extraParams /**< LALDict struct */
1183)
1184{
1185 UNUSED int retcode;
1186
1187 /* Use an auxiliar laldict to not overwrite the input argument */
1188 LALDict *extraParams_aux;
1189
1190 /* sanity checks on input parameters: check pointers, etc. */
1191
1192 /* Check inputs for sanity */
1193 XLAL_CHECK(m1_SI > 0, XLAL_EDOM, "m1 must be positive.\n");
1194 XLAL_CHECK(m2_SI > 0, XLAL_EDOM, "m2 must be positive.\n");
1195 XLAL_CHECK(fabs(chi1z) <= 1.0, XLAL_EDOM, "Aligned spin chi1z=%g \
1196must be <= 1 in magnitude!\n", chi1z);
1197 XLAL_CHECK(fabs(chi2z) <= 1.0, XLAL_EDOM, "Aligned spin chi2z=%g \
1198must be <= 1 in magnitude!\n", chi2z);
1199 XLAL_CHECK(f_ref >= 0, XLAL_EDOM, "Reference frequency must be \
1200positive.\n");
1201
1202 /* setup ModeArray */
1203 if (extraParams == NULL){
1204 extraParams_aux = XLALCreateDict();
1205 }
1206 else{
1207 extraParams_aux = XLALDictDuplicate(extraParams);
1208 }
1209 extraParams_aux = IMRPhenomHM_setup_mode_array(extraParams_aux);
1210 LALValue *ModeArray = XLALSimInspiralWaveformParamsLookupModeArray(extraParams_aux);
1211 int rcode = IMRPhenomHM_check_mode_array(ModeArray);
1212 XLAL_CHECK(XLAL_SUCCESS == rcode, rcode, "IMRPhenomHM_check_mode_array failed");
1213
1214 /* setup frequency sequency */
1215 REAL8Sequence *amps = NULL;
1216 REAL8Sequence *phases = NULL;
1217 REAL8Sequence *freqs_geom = NULL; /* freqs is in geometric units */
1218
1219 LIGOTimeGPS tC = LIGOTIMEGPSZERO; // = {0, 0}
1220
1221 /* setup PhenomHM model storage struct / structs */
1222 /* Compute quantities/parameters related to PhenomD only once and store them */
1223 PhenomHMStorage *pHM;
1224 pHM = XLALMalloc(sizeof(PhenomHMStorage));
1225 retcode = 0;
1226 retcode = init_PhenomHM_Storage(
1227 pHM,
1228 m1_SI,
1229 m2_SI,
1230 chi1x,
1231 chi1y,
1232 chi1z,
1233 chi2x,
1234 chi2y,
1235 chi2z,
1236 freqs,
1237 deltaF,
1238 f_ref,
1239 phiRef);
1240 XLAL_CHECK(XLAL_SUCCESS == retcode, XLAL_EFUNC, "init_PhenomHM_Storage \
1241failed");
1242
1243 /* Two possibilities */
1244 if (pHM->freq_is_uniform == 1)
1245 { /* 1. uniformly spaced */
1246 XLAL_PRINT_INFO("freq_is_uniform = True\n");
1247
1248 freqs = XLALCreateREAL8Sequence(pHM->npts);
1249 phases = XLALCreateREAL8Sequence(pHM->npts);
1250 amps = XLALCreateREAL8Sequence(pHM->npts);
1251
1252 for (size_t i = 0; i < pHM->npts; i++)
1253 { /* populate the frequency unitformly from zero - this is the standard
1254 convention we use when generating waveforms in LAL. */
1255 freqs->data[i] = i * pHM->deltaF; /* This is in Hz */
1256 phases->data[i] = 0; /* initalise all phases to zero. */
1257 amps->data[i] = 0; /* initalise all amps to zero. */
1258 }
1259 /* coalesce at t=0 */
1260 XLAL_CHECK(
1261 XLALGPSAdd(&tC, -1. / pHM->deltaF),
1262 XLAL_EFUNC,
1263 "Failed to shift coalescence time to t=0,\
1264tried to apply shift of -1.0/deltaF with deltaF=%g.",
1265 pHM->deltaF);
1266 }
1267 else if (pHM->freq_is_uniform == 0)
1268 { /* 2. arbitrarily space */
1269 XLAL_PRINT_INFO("freq_is_uniform = False\n");
1270 freqs = pHM->freqs; /* This is in Hz */
1271 phases = XLALCreateREAL8Sequence(freqs->length);
1272 amps = XLALCreateREAL8Sequence(freqs->length);
1273 for (size_t i = 0; i < pHM->npts; i++)
1274 {
1275 phases->data[i] = 0; /* initalise all phases to zero. */
1276 amps->data[i] = 0; /* initalise all phases to zero. */
1277 }
1278 }
1279 else
1280 {
1281 XLAL_ERROR(XLAL_EDOM, "freq_is_uniform = %i and should be either 0 or 1.", pHM->freq_is_uniform);
1282 }
1283
1284 /* PhenomD functions take geometric frequencies */
1285 freqs_geom = XLALCreateREAL8Sequence(pHM->npts);
1286 for (size_t i = 0; i < pHM->npts; i++)
1287 {
1288 freqs_geom->data[i] = XLALSimPhenomUtilsHztoMf(freqs->data[i], pHM->Mtot); /* initalise all phases to zero. */
1289 }
1290
1291 PhenDAmpAndPhasePreComp pDPreComp22;
1293 &pDPreComp22,
1294 pHM->m1,
1295 pHM->m2,
1296 pHM->chi1x,
1297 pHM->chi1y,
1298 pHM->chi1z,
1299 pHM->chi2x,
1300 pHM->chi2y,
1301 pHM->chi2z,
1302 pHM->Rholm[2][2],
1303 pHM->Taulm[2][2],
1304 extraParams_aux);
1305 if (retcode != XLAL_SUCCESS)
1306 {
1307 XLALPrintError("XLAL Error - IMRPhenomDSetupAmpAndPhaseCoefficients failed\n");
1309 }
1310
1311 /* compute the reference phase shift need to align the waveform so that
1312 the phase is equal to phiRef at the reference frequency f_ref. */
1313 /* the phase shift is computed by evaluating the phase of the
1314 (l,m)=(2,2) mode.
1315 phi0 is the correction we need to add to each mode. */
1316 REAL8 phi_22_at_f_ref = IMRPhenomDPhase_OneFrequency(pHM->Mf_ref, pDPreComp22,
1317 1.0, 1.0);
1318 REAL8 phi0 = 0.5 * phi_22_at_f_ref + phiRef;
1319
1320 /* loop over modes */
1321 /* at this point ModeArray should contain the list of modes
1322 * and therefore if NULL then something is wrong and abort.
1323 */
1324 if (ModeArray == NULL)
1325 {
1326 XLAL_ERROR(XLAL_EDOM, "ModeArray is NULL when it shouldn't be. Aborting.\n");
1327 }
1328 for (UINT4 ell = 2; ell < L_MAX_PLUS_1; ell++)
1329 {
1330 for (INT4 mm = 1; mm < (INT4)ell + 1; mm++)
1331 { /* loop over only positive m is intentional. negative m added automatically */
1332 /* first check if (l,m) mode is 'activated' in the ModeArray */
1333 /* if activated then generate the mode, else skip this mode. */
1334 if (XLALSimInspiralModeArrayIsModeActive(ModeArray, ell, mm) != 1)
1335 { /* skip mode */
1336 XLAL_PRINT_INFO("SKIPPING ell = %i mm = %i\n", ell, mm);
1337 continue;
1338 } /* else: generate mode */
1339 XLAL_PRINT_INFO("generateing ell = %i mm = %i\n", ell, mm);
1340
1341 COMPLEX16FrequencySeries *hlm = XLALCreateCOMPLEX16FrequencySeries("hlm: FD waveform", &tC, 0.0, pHM->deltaF, &lalStrainUnit, pHM->npts);
1342 memset(hlm->data->data, 0, pHM->npts * sizeof(COMPLEX16));
1344 retcode = 0;
1345 retcode = IMRPhenomHMEvaluateOnehlmMode(&hlm,
1346 amps, phases,
1347 freqs_geom,
1348 pHM,
1349 ell, mm,
1350 phi0,
1351 extraParams_aux);
1352 XLAL_CHECK(XLAL_SUCCESS == retcode,
1353 XLAL_EFUNC, "IMRPhenomHMEvaluateOnehlmMode failed");
1354
1355 *hlms = XLALSphHarmFrequencySeriesAddMode(*hlms, hlm, ell, mm);
1356
1358 }
1359 }
1360
1361 /* cleanup */
1362 XLALDestroyREAL8Sequence(freqs_geom);
1363 XLALDestroyValue(ModeArray);
1364
1365 if (pHM->freq_is_uniform == 1)
1366 { /* 1. uniformly spaced */
1370 }
1371 else if (pHM->freq_is_uniform == 0)
1372 { /* 2. arbitrarily space */
1375 }
1376 else
1377 {
1378 XLAL_ERROR(XLAL_EDOM, "freq_is_uniform = %i and should be either 0 or 1.", pHM->freq_is_uniform);
1379 }
1380
1381 LALFree(pHM);
1382
1383 XLALDestroyDict(extraParams_aux);
1384
1385 return XLAL_SUCCESS;
1386}
1387
1388/**
1389 * Function to compute the one hlm mode.
1390 * Note this is not static so that IMRPhenomPv3HM
1391 * can also use this function
1392 */
1394 UNUSED COMPLEX16FrequencySeries **hlm, /**< [out] One hlm mode */
1395 UNUSED REAL8Sequence *amps, /**< amplitude frequency sequence */
1396 UNUSED REAL8Sequence *phases, /**< phase frequency sequence */
1397 UNUSED REAL8Sequence *freqs_geom, /**< dimensionless frequency sequence */
1398 UNUSED PhenomHMStorage *pHM, /**< PhenomHMStorage struct */
1399 UNUSED UINT4 ell, /**< ell spherical harmonic number */
1400 UNUSED INT4 mm, /**< m spherical harmonic number */
1401 UNUSED REAL8 phi0, /**< phase shift needed to align waveform to phiRef at f_ref. */
1402 UNUSED LALDict *extraParams /**< LALDict struct */
1403)
1404{
1405 int retcode;
1406
1407 /* generate phase */
1408 retcode = 0;
1409 retcode = IMRPhenomHMPhase(
1410 phases,
1411 freqs_geom,
1412 pHM,
1413 ell, mm,
1414 extraParams);
1415 XLAL_CHECK(XLAL_SUCCESS == retcode,
1416 XLAL_EFUNC, "IMRPhenomHMPhase failed");
1417
1418 /* generate amplitude */
1419 retcode = 0;
1420 retcode = IMRPhenomHMAmplitude(
1421 amps,
1422 freqs_geom,
1423 pHM,
1424 ell, mm,
1425 extraParams);
1426 XLAL_CHECK(XLAL_SUCCESS == retcode,
1427 XLAL_EFUNC, "IMRPhenomHMAmplitude failed");
1428
1429 /* compute time shift */
1431 pHM->eta, pHM->chi1z, pHM->chi2z,
1432 pHM->finspin, extraParams);
1433
1434 REAL8 phase_term1 = 0.0;
1435 REAL8 phase_term2 = 0.0;
1436 REAL8 Mf = 0.0; /* geometric frequency */
1437 /* combine together to make hlm */
1438 //loop over hlm COMPLEX16FrequencySeries
1439 for (size_t i = pHM->ind_min; i < pHM->ind_max; i++)
1440 {
1441 Mf = freqs_geom->data[i];
1442 phase_term1 = - t0 * (Mf - pHM->Mf_ref);
1443 phase_term2 = phases->data[i] - (mm * phi0);
1444 ((*hlm)->data->data)[i] = amps->data[i] * cexp(-I * (phase_term1 + phase_term2));
1445 }
1446
1447 /* cleanup */
1448
1449 return XLAL_SUCCESS;
1450}
1451
1452/**
1453 * returns IMRPhenomHM amplitude evaluated at a set of input frequencies
1454 * for the l,m mode
1455 */
1457 UNUSED REAL8Sequence *amps, /**< [out] amplitude frequency sequence */
1458 UNUSED REAL8Sequence *freqs_geom, /**< dimensionless frequency sequence */
1459 UNUSED PhenomHMStorage *pHM, /**< PhenomHMStorage struct */
1460 UNUSED UINT4 ell, /**< ell spherical harmonic number */
1461 UNUSED INT4 mm, /**< m spherical harmonic number */
1462 UNUSED LALDict *extraParams /**< LALDict struct */
1463)
1464{
1465 int retcode;
1466
1467 /* scale input frequencies according to PhenomHM model */
1468 /* LL: Map the input domain (frequencies) for this ell mm multipole
1469 to those appropirate for the ell=|mm| multipole */
1470 REAL8Sequence *freqs_amp = XLALCreateREAL8Sequence(freqs_geom->length);
1471 for (UINT4 i = 0; i < freqs_amp->length; i++)
1472 {
1473 freqs_amp->data[i] = IMRPhenomHMFreqDomainMap(
1474 freqs_geom->data[i], ell, mm, pHM, AmpFlagTrue);
1475 }
1476
1477 /* LL: Compute the PhenomD Amplitude at the mapped l=m=2 fequencies */
1478 retcode = 0;
1480 amps,
1481 freqs_amp,
1482 pHM->ind_min, pHM->ind_max,
1483 pHM->m1, pHM->m2,
1484 pHM->chi1x,
1485 pHM->chi1y,
1486 pHM->chi1z,
1487 pHM->chi2x,
1488 pHM->chi2y,
1489 pHM->chi2z);
1490 XLAL_CHECK(XLAL_SUCCESS == retcode,
1491 XLAL_EFUNC, "IMRPhenomDAmpFrequencySequence failed");
1492
1493 /*
1494 LL: Here we map the ampliude's range using two steps:
1495 (1) We divide by the leading order l=m=2 behavior, and then
1496 scale in the expected PN behavior for the multipole of interest.
1497 NOTE that this step is done at the mapped frequencies,
1498 which results in smooth behavior despite the sharp featured of the domain map.
1499 There are other (perhaps more intuitive) options for mapping the amplitudes,
1500 but these do not have the desired smooth features.
1501 (2) An additional scaling is needed to recover the desired PN ampitude.
1502 This is needed becuase only frequencies appropriate for the dominant
1503 quadrupole have been used thusly, so the current answer does not
1504 conform to PN expectations for inspiral.
1505 This is trikier than described here, so please give it a deeper think.
1506 */
1507
1508 int status_in_for = XLAL_SUCCESS;
1509 for (UINT4 i = 0; i < freqs_amp->length; i++)
1510 {
1511 PhenomHMUsefulPowers powers_of_freq_amp;
1512 status_in_for = PhenomHM_init_useful_powers(
1513 &powers_of_freq_amp, freqs_amp->data[i]);
1514 if (XLAL_SUCCESS != status_in_for)
1515 {
1516 XLALPrintError("PhenomHM_init_useful_powers failed for Mf, status_in_for=%d", status_in_for);
1517 retcode = status_in_for;
1518 }
1519 //new
1520
1521 /* LL: Calculate the corrective factor for step #2 */
1522
1523 double beta_term1 = IMRPhenomHMOnePointFiveSpinPN(
1524 freqs_geom->data[i],
1525 ell,
1526 mm,
1527 pHM->m1,
1528 pHM->m2,
1529 pHM->chi1z,
1530 pHM->chi2z);
1531
1532 double beta=0.;
1533 double beta_term2=0.;
1534 double HMamp_term1=1.;
1535 double HMamp_term2=1.;
1536 //HACK to fix equal black hole case producing NaNs.
1537 //More elegant solution needed.
1538 if (beta_term1 == 0.){
1539 beta = 0.;
1540 } else {
1541
1542 beta_term2 = IMRPhenomHMOnePointFiveSpinPN(2.0 * freqs_geom->data[i] / mm, ell, mm, pHM->m1, pHM->m2, pHM->chi1z, pHM->chi2z);
1543 beta = beta_term1 / beta_term2;
1544
1545 /* LL: Apply steps #1 and #2 */
1546 HMamp_term1 = IMRPhenomHMOnePointFiveSpinPN(
1547 freqs_amp->data[i],
1548 ell,
1549 mm,
1550 pHM->m1,
1551 pHM->m2,
1552 pHM->chi1z,
1553 pHM->chi2z);
1554 HMamp_term2 = IMRPhenomHMOnePointFiveSpinPN(freqs_amp->data[i], 2, 2, pHM->m1, pHM->m2, 0.0, 0.0);
1555
1556 }
1557
1558 //HMamp is computed here
1559 amps->data[i] *= beta * HMamp_term1 / HMamp_term2;
1560 }
1561
1562 /* cleanup */
1563 XLALDestroyREAL8Sequence(freqs_amp);
1564
1565 return XLAL_SUCCESS;
1566}
1567/**
1568 * returns IMRPhenomHM phase evaluated at a set of input frequencies
1569 * for the l,m mode
1570 */
1572 UNUSED REAL8Sequence *phases, /**< [out] phase frequency sequence */
1573 UNUSED REAL8Sequence *freqs_geom, /**< dimensionless frequency sequence */
1574 UNUSED PhenomHMStorage *pHM, /**< PhenomHMStorage struct */
1575 UNUSED UINT4 ell, /**< ell spherical harmonic number */
1576 UNUSED INT4 mm, /**< m spherical harmonic number */
1577 UNUSED LALDict *extraParams /**< LALDict struct */
1578)
1579{
1580 int retcode;
1581
1583 retcode = 0;
1584 retcode = IMRPhenomHMPhasePreComp(&q, ell, mm, pHM, extraParams);
1585 if (retcode != XLAL_SUCCESS)
1586 {
1587 XLALPrintError("XLAL Error - IMRPhenomHMPhasePreComp failed\n");
1589 }
1590 REAL8 Rholm = pHM->Rholm[ell][mm];
1591 REAL8 Taulm = pHM->Taulm[ell][mm];
1592
1593 PhenDAmpAndPhasePreComp pDPreComp;
1595 &pDPreComp,
1596 pHM->m1,
1597 pHM->m2,
1598 pHM->chi1x,
1599 pHM->chi1y,
1600 pHM->chi1z,
1601 pHM->chi2x,
1602 pHM->chi2y,
1603 pHM->chi2z,
1604 Rholm,
1605 Taulm,
1606 extraParams);
1607 if (retcode != XLAL_SUCCESS)
1608 {
1609 XLALPrintError("XLAL Error - IMRPhenomDSetupAmpAndPhaseCoefficients failed\n");
1611 }
1612
1613 REAL8 Mf_wf = 0.0;
1614 REAL8 Mf = 0.0;
1615 REAL8 Mfr = 0.0;
1616 REAL8 tmpphaseC = 0.0;
1617 for (UINT4 i = pHM->ind_min; i < pHM->ind_max; i++)
1618 {
1619 /* Add complex phase shift depending on 'm' mode */
1620 phases->data[i] = cShift[mm];
1621 Mf_wf = freqs_geom->data[i];
1622 // This if ladder is in the mathematica function HMPhase. PhenomHMDev.nb
1623 if (!(Mf_wf > q.fi))
1624 { /* in mathematica -> IMRPhenDPhaseA */
1625 Mf = q.ai * Mf_wf + q.bi;
1626 phases->data[i] += IMRPhenomDPhase_OneFrequency(Mf, pDPreComp, Rholm, Taulm) / q.ai;
1627 }
1628 else if (!(Mf_wf > q.fr))
1629 { /* in mathematica -> IMRPhenDPhaseB */
1630 Mf = q.am * Mf_wf + q.bm;
1631 phases->data[i] += IMRPhenomDPhase_OneFrequency(Mf, pDPreComp, Rholm, Taulm) / q.am - q.PhDBconst + q.PhDBAterm;
1632 }
1633 else if ((Mf_wf > q.fr))
1634 { /* in mathematica -> IMRPhenDPhaseC */
1635 Mfr = q.am * q.fr + q.bm;
1636 tmpphaseC = IMRPhenomDPhase_OneFrequency(Mfr, pDPreComp, Rholm, Taulm) / q.am - q.PhDBconst + q.PhDBAterm;
1637 Mf = q.ar * Mf_wf + q.br;
1638 phases->data[i] += IMRPhenomDPhase_OneFrequency(Mf, pDPreComp, Rholm, Taulm) / q.ar - q.PhDCconst + tmpphaseC;
1639 }
1640 else
1641 {
1642 XLALPrintError("XLAL_ERROR - should not get here - in function IMRPhenomHMPhase");
1644 }
1645 }
1646
1647 return XLAL_SUCCESS;
1648
1649
1650}
void XLALDestroyDict(LALDict *dict)
LALDict * XLALDictDuplicate(const LALDict *orig)
LALDict * XLALCreateDict(void)
#define LALFree(p)
UNUSED REAL8 IMRPhenomDPhase_OneFrequency(REAL8 Mf, PhenDAmpAndPhasePreComp pD, REAL8 Rholm, REAL8 Taulm)
Function to return the phenomD phase using the IMRPhenomDSetupAmpAndPhaseCoefficients struct.
int IMRPhenomDAmpFrequencySequence(REAL8Sequence *amps, REAL8Sequence *freqs, size_t ind_min, size_t ind_max, REAL8 m1, REAL8 m2, REAL8 chi1x, REAL8 chi1y, REAL8 chi1z, REAL8 chi2x, REAL8 chi2y, REAL8 chi2z)
Helper function used in PhenomHM and PhenomPv3HM Returns the phenomD amplitude.
REAL8 IMRPhenomDComputet0(REAL8 eta, REAL8 chi1z, REAL8 chi2z, REAL8 finspin, LALDict *extraParams)
computes the time shift as the approximate time of the peak of the 22 mode.
int IMRPhenomDSetupAmpAndPhaseCoefficients(PhenDAmpAndPhasePreComp *pDPreComp, REAL8 m1, REAL8 m2, REAL8 chi1x, REAL8 chi1y, REAL8 chi1z, REAL8 chi2x, REAL8 chi2y, REAL8 chi2z, const REAL8 Rholm, const REAL8 Taulm, LALDict *extraParams)
Function to compute the amplitude and phase coefficients for PhenomD Used to optimise the calls to IM...
#define PHI_fJoin_INS
Dimensionless frequency (Mf) at which the inspiral phase switches to the intermediate phase.
#define AMP_fJoin_INS
Dimensionless frequency (Mf) at which the inspiral amplitude switches to the intermediate amplitude.
Internal function for IMRPhenomD phenomenological waveform model. See LALSimIMRPhenom....
static double fdamp(double eta, double chi1, double chi2, double finspin)
fdamp is the complex part of the ringdown frequency 1508.07250 figure 9
int PhenomHM_init_useful_powers(PhenomHMUsefulPowers *p, REAL8 number)
must be called before the first usage of *p
double IMRPhenomHMTi(REAL8 Mf, const INT4 mm)
mathematica function Ti domain mapping function - inspiral
COMPLEX16 IMRPhenomHMOnePointFiveSpinPN(REAL8 fM, INT4 l, INT4 m, REAL8 M1, REAL8 M2, REAL8 X1z, REAL8 X2z)
Define function for FD PN amplitudes.
static int IMRPhenomHM_check_mode_array(LALValue *ModeArray)
Reads in a ModeArray and checks that it is valid.
static int init_PhenomHM_Storage(PhenomHMStorage *p, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi1x, const REAL8 chi1y, const REAL8 chi1z, const REAL8 chi2x, const REAL8 chi2y, const REAL8 chi2z, REAL8Sequence *freqs, const REAL8 deltaF, const REAL8 f_ref, const REAL8 phiRef)
Precompute a bunch of PhenomHM related quantities and store them filling in a PhenomHMStorage variabl...
int IMRPhenomHMEvaluateOnehlmMode(UNUSED COMPLEX16FrequencySeries **hlm, UNUSED REAL8Sequence *amps, UNUSED REAL8Sequence *phases, UNUSED REAL8Sequence *freqs_geom, UNUSED PhenomHMStorage *pHM, UNUSED UINT4 ell, UNUSED INT4 mm, UNUSED REAL8 phi0, UNUSED LALDict *extraParams)
Function to compute the one hlm mode.
int IMRPhenomHMFreqDomainMapParams(REAL8 *a, REAL8 *b, REAL8 *fi, REAL8 *fr, REAL8 *f1, const REAL8 flm, const INT4 ell, const INT4 mm, PhenomHMStorage *pHM, const int AmpFlag)
helper function for IMRPhenomHMFreqDomainMap
int init_IMRPhenomHMGet_FrequencyBounds_storage(PhenomHMFrequencyBoundsStorage *p, REAL8Sequence *freqs, REAL8 Mtot, REAL8 deltaF, REAL8 f_ref_in)
derive frequency variables for PhenomHM based on input.
static const double cShift[7]
double IMRPhenomHMFreqDomainMap(REAL8 Mflm, const INT4 ell, const INT4 mm, PhenomHMStorage *pHM, const int AmpFlag)
IMRPhenomHMFreqDomainMap Input waveform frequency in Geometric units (Mflm) and computes what frequen...
int IMRPhenomHMSlopeAmAndBm(double *Am, double *Bm, const INT4 mm, REAL8 fi, REAL8 fr, REAL8 Mf_RD_22, REAL8 Mf_RD_lm, const INT4 AmpFlag, const INT4 ell, PhenomHMStorage *pHM)
helper function for IMRPhenomHMFreqDomainMap
int IMRPhenomHMMapParams(REAL8 *a, REAL8 *b, REAL8 flm, REAL8 fi, REAL8 fr, REAL8 Ai, REAL8 Bi, REAL8 Am, REAL8 Bm, REAL8 Ar, REAL8 Br)
helper function for IMRPhenomHMFreqDomainMap
int IMRPhenomHMPhasePreComp(HMPhasePreComp *q, const INT4 ell, const INT4 mm, PhenomHMStorage *pHM, UNUSED LALDict *extraParams)
int IMRPhenomHMGetRingdownFrequency(REAL8 *fringdown, REAL8 *fdamp, UINT4 ell, INT4 mm, REAL8 finalmass, REAL8 finalspin)
returns the real and imag parts of the complex ringdown frequency for the (l,m) mode.
int IMRPhenomHMPhase(UNUSED REAL8Sequence *phases, UNUSED REAL8Sequence *freqs_geom, UNUSED PhenomHMStorage *pHM, UNUSED UINT4 ell, UNUSED INT4 mm, UNUSED LALDict *extraParams)
returns IMRPhenomHM phase evaluated at a set of input frequencies for the l,m mode
LALDict * IMRPhenomHM_setup_mode_array(LALDict *extraParams)
read in a LALDict.
int PhenomHM_init_useful_mf_powers(PhenomHMUsefulMfPowers *p, REAL8 number)
must be called before the first usage of *p
int IMRPhenomHMAmplitude(UNUSED REAL8Sequence *amps, UNUSED REAL8Sequence *freqs_geom, UNUSED PhenomHMStorage *pHM, UNUSED UINT4 ell, UNUSED INT4 mm, UNUSED LALDict *extraParams)
returns IMRPhenomHM amplitude evaluated at a set of input frequencies for the l,m mode
int XLALSimIMRPhenomHMGethlmModes(UNUSED SphHarmFrequencySeries **hlms, UNUSED REAL8Sequence *freqs, UNUSED REAL8 m1_SI, UNUSED REAL8 m2_SI, UNUSED REAL8 chi1x, UNUSED REAL8 chi1y, UNUSED REAL8 chi1z, UNUSED REAL8 chi2x, UNUSED REAL8 chi2y, UNUSED REAL8 chi2z, UNUSED const REAL8 phiRef, UNUSED const REAL8 deltaF, UNUSED REAL8 f_ref, UNUSED LALDict *extraParams)
XLAL function that returns a SphHarmFrequencySeries object containing all the hlm modes requested.
UINT4 IMRPhenomHM_is_freq_uniform(REAL8Sequence *freqs, REAL8 deltaF)
helper function to easily check if the input frequency sequence is uniformly space or a user defined ...
int IMRPhenomHMCore(UNUSED COMPLEX16FrequencySeries **hptilde, UNUSED COMPLEX16FrequencySeries **hctilde, REAL8Sequence *freqs, REAL8 m1_SI, REAL8 m2_SI, REAL8 chi1z, REAL8 chi2z, const REAL8 distance, const REAL8 inclination, const REAL8 phiRef, const REAL8 deltaF, REAL8 f_ref, LALDict *extraParams)
internal function that returns h+ and hx.
double IMRPhenomHMTrd(REAL8 Mf, REAL8 Mf_RD_22, REAL8 Mf_RD_lm, const INT4 AmpFlag, const INT4 ell, const INT4 mm, PhenomHMStorage *pHM)
domain mapping function - ringdown
#define L_MAX_PLUS_1
Highest ell multipole PhenomHM models + 1.
#define PHENOMHM_DEFAULT_MF_MAX
Dimensionless frequency (Mf) at which define the end of the waveform.
#define MAX_ALLOWED_ETA
eta is the symmetric mass-ratio.
#define AmpFlagTrue
int PhenomInternal_IMRPhenomHMFDAddMode(COMPLEX16FrequencySeries *hptilde, COMPLEX16FrequencySeries *hctilde, COMPLEX16FrequencySeries *hlmtilde, REAL8 theta, REAL8 phi, INT4 l, INT4 m, INT4 sym)
helper function to multiple hlm with Ylm.
size_t PhenomInternal_NextPow2(const size_t n)
Return the closest higher power of 2.
int PhenomInternal_PrecessingSpinEnforcePrimaryIsm1(REAL8 *m1, REAL8 *m2, REAL8 *chi1x, REAL8 *chi1y, REAL8 *chi1z, REAL8 *chi2x, REAL8 *chi2y, REAL8 *chi2z)
Given m1 with spins (chi1x, chi1y, chi1z) and m2 with spins (chi2x,chi2y,chi2z).
void PhenomInternal_nudge(REAL8 *x, REAL8 X, REAL8 epsilon)
If x and X are approximately equal to relative accuracy epsilon then set x = X.
static double double delta
double XLALSimPhenomUtilsFDamp0(REAL8 Mtot_Msun, REAL8 distance)
compute the frequency domain amplitude pre-factor
double XLALSimPhenomUtilsIMRPhenomDFinalMass(REAL8 m1, REAL8 m2, REAL8 chi1z, REAL8 chi2z)
Helper function used in PhenomHM and PhenomPv3HM Returns the final mass from the fit used in PhenomD.
REAL8 XLALSimPhenomUtilsPhenomPv2FinalSpin(const REAL8 m1, const REAL8 m2, const REAL8 chi1_l, const REAL8 chi2_l, const REAL8 chip)
Wrapper for final-spin formula based on:
REAL8 XLALSimPhenomUtilsChiP(const REAL8 m1, const REAL8 m2, const REAL8 s1x, const REAL8 s1y, const REAL8 s2x, const REAL8 s2y)
Function to compute the effective precession parameter chi_p (1408.1810)
double XLALSimPhenomUtilsMftoHz(REAL8 Mf, REAL8 Mtot_Msun)
Convert from geometric frequency to frequency in Hz.
double XLALSimPhenomUtilsHztoMf(REAL8 fHz, REAL8 Mtot_Msun)
Convert from frequency in Hz to geometric frequency.
static double beta(const double a, const double b, const sysq *system)
Internal function that computes the spin-orbit couplings.
LALValue * XLALSimInspiralWaveformParamsLookupModeArray(LALDict *params)
int XLALSimInspiralWaveformParamsInsertModeArray(LALDict *params, LALValue *value)
double SimRingdownCW_KAPPA(double jf, int l, int m)
complex double SimRingdownCW_CW07102016(double kappa, int l, int input_m, int n)
void XLALDestroyValue(LALValue *value)
int l
Definition: bh_qnmode.c:135
REAL8 M
Definition: bh_qnmode.c:133
double i
Definition: bh_ringdown.c:118
double e
Definition: bh_ringdown.c:117
COMPLEX16FrequencySeries * XLALCreateCOMPLEX16FrequencySeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaF, const LALUnit *sampleUnits, size_t length)
void XLALDestroyCOMPLEX16FrequencySeries(COMPLEX16FrequencySeries *series)
#define LAL_PI_2
#define LAL_MSUN_SI
#define LAL_PI
#define LIGOTIMEGPSZERO
double complex COMPLEX16
double REAL8
uint32_t UINT4
int32_t INT4
void * XLALMalloc(size_t n)
void XLALFree(void *p)
UNUSED int XLALSimIMRPhenomHM(UNUSED COMPLEX16FrequencySeries **hptilde, UNUSED COMPLEX16FrequencySeries **hctilde, UNUSED REAL8Sequence *freqs, UNUSED REAL8 m1_SI, UNUSED REAL8 m2_SI, UNUSED REAL8 chi1z, UNUSED REAL8 chi2z, UNUSED const REAL8 distance, UNUSED const REAL8 inclination, UNUSED const REAL8 phiRef, UNUSED const REAL8 deltaF, UNUSED REAL8 f_ref, UNUSED LALDict *extraParams)
Returns h+ and hx in the frequency domain.
int XLALSimInspiralModeArrayIsModeActive(LALValue *modes, unsigned l, int m)
LALValue * XLALSimInspiralModeArrayActivateMode(LALValue *modes, unsigned l, int m)
LALValue * XLALSimInspiralCreateModeArray(void)
Create a LALValue pointer to store the mode array.
COMPLEX16FrequencySeries * XLALSphHarmFrequencySeriesGetMode(SphHarmFrequencySeries *ts, UINT4 l, INT4 m)
Get the time series of a waveform's (l,m) spherical harmonic mode from a SphHarmFrequencySeries linke...
SphHarmFrequencySeries * XLALSphHarmFrequencySeriesAddMode(SphHarmFrequencySeries *appended, const COMPLEX16FrequencySeries *inmode, UINT4 l, INT4 m)
Prepend a node to a linked list of SphHarmFrequencySeries, or create a new head.
void XLALDestroySphHarmFrequencySeries(SphHarmFrequencySeries *ts)
Delete list from current pointer to the end of the list.
static const INT4 m
static const INT4 q
static const INT4 a
void XLALDestroyREAL8Sequence(REAL8Sequence *sequence)
REAL8Sequence * XLALCreateREAL8Sequence(size_t length)
const LALUnit lalStrainUnit
const LALUnit lalSecondUnit
LALUnit * XLALUnitMultiply(LALUnit *output, const LALUnit *unit1, const LALUnit *unit2)
#define XLAL_PRINT_INFO(...)
#define XLAL_ERROR(...)
#define XLAL_CHECK(assertion,...)
#define XLAL_PRINT_WARNING(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_PRINT_ERROR(...)
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFUNC
XLAL_EDOM
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
p
COMPLEX16Sequence * data
COMPLEX16 * data
Structure holding Higher Mode Phase pre-computations.
A struct the store the amplitude and phase structs for phenomD.
Structure storing pre-determined quantities that describe the frequency array and tells you over whic...
size_t npts
number of points in waveform array
size_t ind_min
start index containing non-zero values
UINT4 freq_is_uniform
If = 1 then assume uniform spaced, If = 0 then assume arbitrarily spaced.
size_t ind_max
end index containing non-zero values
Structure storing pre-determined quantities complying to the conventions of the PhenomHM model.
REAL8 Rholm[L_MAX_PLUS_1][L_MAX_PLUS_1]
ratio of (2,2) mode to (l,m) mode ringdown frequency
REAL8 Taulm[L_MAX_PLUS_1][L_MAX_PLUS_1]
ratio of (l,m) mode to (2,2) mode damping time
REAL8 chi1x
x-component of the dimensionless spin of object 1 w.r.t.
REAL8 chi2y
y-component of the dimensionless spin of object 2 w.r.t.
size_t npts
number of points in waveform array
REAL8Sequence * freqs
REAL8 chi2z
z-component of the dimensionless spin of object 2 w.r.t.
REAL8 m2
mass of lighter body in solar masses
REAL8 chi2x
x-component of the dimensionless spin of object 2 w.r.t.
REAL8 chi1z
z-component of the dimensionless spin of object 1 w.r.t.
REAL8 Mtot
total mass in solar masses
REAL8 m1
mass of larger body in solar masses
REAL8 PhenomHMfring[L_MAX_PLUS_1][L_MAX_PLUS_1]
REAL8 Mf_ref
reference frequnecy in geometric units
REAL8 chi1y
y-component of the dimensionless spin of object 1 w.r.t.
UINT4 freq_is_uniform
If = 1 then assume uniform spaced, If = 0 then assume arbitrarily spaced.
Useful powers of Mf: 1/6, 1/3, 2/3, 4/3, 5/3, 2, 7/3, 8/3, -7/6, -5/6, -1/2, -1/6,...
useful powers in GW waveforms: 1/6, 1/3, 2/3, 4/3, 5/3, 2, 7/3, 8/3, -1, -1/6, -7/6,...
REAL8 * data