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
LALSimIMRPhenomC_internals.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2012 Prayush Kumar, Frank Ohme
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/* The paper refered to here as the Main paper, is Phys. Rev. D 82, 064016 (2010)
22 * */
23
24
26
27/**
28 *
29 * Internal functions
30 *
31 * */
32
33/*********************************************************************/
34/* Compute PN coefficients for non-precessing binaries */
35/* Ref. Eq.(A3)-(A5) of http://arxiv.org/pdf/1005.3306v3.pdf */
36/* */
37/*********************************************************************/
39 const REAL8 m1, /**< mass of companion 1 (solar masses) */
40 const REAL8 m2, /**< mass of companion 2 (solar masses) */
41 const REAL8 chi, /**< Reduced spin of the binary, defined in the main paper */
42 LALDict *extraParams) /**< linked list containing the extra testing GR parameters */
43{
44
47 memset(p, 0, sizeof(BBHPhenomCParams));
48
49 /* calculate the total mass and symmetric mass ratio */
50 const REAL8 M = m1 + m2;
51 const REAL8 eta = m1 * m2 / (M * M);
52 const REAL8 piM = M * LAL_PI * LAL_MTSUN_SI;
53 const REAL8 eta2 = eta*eta;
54 const REAL8 chisum = chi + chi;
55 const REAL8 chiprod = chi * chi;
56 const REAL8 chi2 = chi * chi;
57
58 // Store the total Mass of the system in seconds
59 p->m_sec = M * LAL_MTSUN_SI;
60 p->piM = piM;
61
62 /* Calculate the PN phasing terms */
63 p->pfaN = 3.0/(128.0 * eta);
64 p->pfa1 = 0.0;
65 p->pfa2 = (3715./756.) + (55.*eta/9.0);
66 p->pfa3 = -16.0*LAL_PI + (113./3.)*chi - 38.*eta*chisum/3.;
67 p->pfa4 = (152.93365/5.08032) - 50.*chi2 + eta*(271.45/5.04 + 1.25*chiprod) +
68 3085.*eta2/72.;
69 p->pfa5 = LAL_PI*(386.45/7.56 - 65.*eta/9.) -
70 chi*(735.505/2.268 + 130.*eta/9.) + chisum*(1285.0*eta/8.1 + 170.*eta2/9.) -
71 10.*chi2*chi/3. + 10.*eta*chi*chiprod;
72 p->pfa6 = 11583.231236531/4.694215680 - 640.0*LAL_PI*LAL_PI/3. -
73 6848.0*LAL_GAMMA/21. - 684.8*log(64.)/6.3 +
74 eta*(2255.*LAL_PI*LAL_PI/12. - 15737.765635/3.048192) +
75 76.055*eta2/1.728 - (127.825*eta2*eta/1.296) +
76 2920.*LAL_PI*chi/3. - (175. - 1490.*eta)*chi2/3. -
77 (1120.*LAL_PI/3. - 1085.*chi/3.)*eta*chisum +
78 (269.45*eta/3.36 - 2365.*eta2/6.)*chiprod;
79
80 p->pfa6log = -6848./63.;
81
82 p->pfa7 = LAL_PI*(770.96675/2.54016 + 378.515*eta/1.512 - 740.45*eta2/7.56) -
83 chi*(20373.952415/3.048192 + 1509.35*eta/2.24 - 5786.95*eta2/4.32) +
84 chisum*(4862.041225*eta/1.524096 + 1189.775*eta2/1.008 - 717.05*eta2*eta/2.16 - 830.*eta*chi2/3. + 35.*eta2*chiprod/3.) -
85 560.*LAL_PI*chi2 + 20.*LAL_PI*eta*chiprod +
86 chi2*chi*(945.55/1.68 - 85.*eta) + chi*chiprod*(396.65*eta/1.68 + 255.*eta2);
87
88 p->pfaN*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRPhi1(extraParams));
90 p->pfa2*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRDChi2(extraParams));
91 p->pfa3*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRDChi3(extraParams));
92 p->pfa4*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRDChi4(extraParams));
93 p->pfa5*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRDChi5(extraParams));
94 p->pfa6*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRDChi6(extraParams));
95 p->pfa6log*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRDChi6L(extraParams));
96 p->pfa7*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRDChi7(extraParams));
97
98 /* Coefficients to calculate xdot, that comes in the fourier amplitude */
99 p->xdotaN = 64.*eta/5.;
100 p->xdota2 = -7.43/3.36 - 11.*eta/4.;
101 p->xdota3 = 4.*LAL_PI - 11.3*chi/1.2 + 19.*eta*chisum/6.;
102 p->xdota4 = 3.4103/1.8144 + 5*chi2 + eta*(13.661/2.016 - chiprod/8.) + 5.9*eta2/1.8;
103 p->xdota5 = -LAL_PI*(41.59/6.72 + 189.*eta/8.) - chi*(31.571/1.008 - 116.5*eta/2.4) +
104 chisum*(21.863*eta/1.008 - 79.*eta2/6.) - 3*chi*chi2/4. +
105 9.*eta*chi*chiprod/4.;
106 p->xdota6 = 164.47322263/1.39708800 - 17.12*LAL_GAMMA/1.05 +
107 16.*LAL_PI*LAL_PI/3 - 8.56*log(16.)/1.05 +
108 eta*(45.1*LAL_PI*LAL_PI/4.8 - 561.98689/2.17728) +
109 5.41*eta2/8.96 - 5.605*eta*eta2/2.592 - 80.*LAL_PI*chi/3. +
110 eta*chisum*(20.*LAL_PI/3. - 113.5*chi/3.6) +
111 chi2*(64.153/1.008 - 45.7*eta/3.6) -
112 chiprod*(7.87*eta/1.44 - 30.37*eta2/1.44);
113
114 p->xdota6log = -856./105.;
115
116 p->xdota7 = -LAL_PI*(4.415/4.032 - 358.675*eta/6.048 - 91.495*eta2/1.512) -
117 chi*(252.9407/2.7216 - 845.827*eta/6.048 + 415.51*eta2/8.64) +
118 chisum*(158.0239*eta/5.4432 - 451.597*eta2/6.048 + 20.45*eta2*eta/4.32 + 107.*eta*chi2/6. - 5.*eta2*chiprod/24.) +
119 12.*LAL_PI*chi2 - chi2*chi*(150.5/2.4 + eta/8.) +
120 chi*chiprod*(10.1*eta/2.4 + 3.*eta2/8.);
121
122 /* Coefficients to compute the time-domain amplitude, which also enters the
123 * fourier amplitude. */
124 p->AN = 8.*eta*sqrt(LAL_PI/5.);
125 p->A2 = (-107. + 55.*eta)/42.;
126 p->A3 = 2.*LAL_PI - 4.*chi/3. + 2.*eta*chisum/3.;
127 p->A4 = -2.173/1.512 - eta*(10.69/2.16 - 2.*chiprod) + 2.047*eta2/1.512;
128 p->A5 = -10.7*LAL_PI/2.1 + eta*(3.4*LAL_PI/2.1);
129
130 p->A5imag = -24.*eta;
131
132 p->A6 = 270.27409/6.46800 - 8.56*LAL_GAMMA/1.05 +
133 2.*LAL_PI*LAL_PI/3. +
134 eta*(4.1*LAL_PI*LAL_PI/9.6 - 27.8185/3.3264) -
135 20.261*eta2/2.772 + 11.4635*eta*eta2/9.9792 -
136 4.28*log(16.)/1.05;
137
138 p->A6log = -428./105.;
139
140 p->A6imag = 4.28*LAL_PI/1.05;
141
142 return p;
143}
144
145/*********************************************************************/
146/* Compute phenomenological parameters for non-precessing binaries */
147/* Ref. Eq.(5.14) of http://arxiv.org/pdf/1005.3306v3.pdf */
148/* and Table II of the same paper. */
149/* */
150/*********************************************************************/
151
153 const REAL8 m1, /**< mass of companion 1 (solar masses) */
154 const REAL8 m2, /**< mass of companion 2 (solar masses) */
155 const REAL8 chi, /**< Reduced spin of the binary, defined in the main paper */
156 LALDict *extraParams) /**< linked list containing the extra testing GR parameters */
157{
158 BBHPhenomCParams *p = NULL;
159 p = ComputeIMRPhenomCParamsSPA( m1, m2, chi, extraParams );
160 if( !p )
162
163 const REAL8 M = m1 + m2;
164 const REAL8 eta = m1 * m2 / (M * M);
165
166 const REAL8 z101 = -2.417e-03;
167 const REAL8 z102 = -1.093e-03;
168 const REAL8 z111 = -1.917e-02;
169 const REAL8 z110 = 7.267e-02;
170 const REAL8 z120 = -2.504e-01;
171
172 const REAL8 z201 = 5.962e-01;
173 const REAL8 z202 = -5.600e-02;
174 const REAL8 z211 = 1.520e-01;
175 const REAL8 z210 = -2.970e+00;
176 const REAL8 z220 = 1.312e+01;
177
178 const REAL8 z301 = -3.283e+01;
179 const REAL8 z302 = 8.859e+00;
180 const REAL8 z311 = 2.931e+01;
181 const REAL8 z310 = 7.954e+01;
182 const REAL8 z320 = -4.349e+02;
183
184 const REAL8 z401 = 1.619e+02;
185 const REAL8 z402 = -4.702e+01;
186 const REAL8 z411 = -1.751e+02;
187 const REAL8 z410 = -3.225e+02;
188 const REAL8 z420 = 1.587e+03;
189
190 const REAL8 z501 = -6.320e+02;
191 const REAL8 z502 = 2.463e+02;
192 const REAL8 z511 = 1.048e+03;
193 const REAL8 z510 = 3.355e+02;
194 const REAL8 z520 = -5.115e+03;
195
196 const REAL8 z601 = -4.809e+01;
197 const REAL8 z602 = -3.643e+02;
198 const REAL8 z611 = -5.215e+02;
199 const REAL8 z610 = 1.870e+03;
200 const REAL8 z620 = 7.354e+02;
201
202 const REAL8 z701 = 4.149e+00;
203 const REAL8 z702 = -4.070e+00;
204 const REAL8 z711 = -8.752e+01;
205 const REAL8 z710 = -4.897e+01;
206 const REAL8 z720 = 6.665e+02;
207
208 const REAL8 z801 = -5.472e-02;
209 const REAL8 z802 = 2.094e-02;
210 const REAL8 z811 = 3.554e-01;
211 const REAL8 z810 = 1.151e-01;
212 const REAL8 z820 = 9.640e-01;
213
214 const REAL8 z901 = -1.235e+00;
215 const REAL8 z902 = 3.423e-01;
216 const REAL8 z911 = 6.062e+00;
217 const REAL8 z910 = 5.949e+00;
218 const REAL8 z920 = -1.069e+01;
219
220 /* Calculate alphas, gamma, deltas from Table II and Eq 5.14 of Main paper */
221 REAL8 eta2 = eta*eta;
222 REAL8 chi2 = chi*chi;
223 REAL8 etachi = eta * chi;
224
225 p->a1 = z101 * chi + z102 * chi2 + z111 * eta * chi + z110 * eta + z120 * eta2;
226 p->a2 = z201 * chi + z202 * chi2 + z211 * eta * chi + z210 * eta + z220 * eta2;
227 p->a3 = z301 * chi + z302 * chi2 + z311 * eta * chi + z310 * eta + z320 * eta2;
228 p->a4 = z401 * chi + z402 * chi2 + z411 * eta * chi + z410 * eta + z420 * eta2;
229 p->a5 = z501 * chi + z502 * chi2 + z511 * eta * chi + z510 * eta + z520 * eta2;
230 p->a6 = z601 * chi + z602 * chi2 + z611 * eta * chi + z610 * eta + z620 * eta2;
231
232 p->g1 = z701 * chi + z702 * chi2 + z711 * eta * chi + z710 * eta + z720 * eta2;
233
234 p->del1 = z801 * chi + z802 * chi2 + z811 * eta * chi + z810 * eta + z820 * eta2;
235 p->del2 = z901 * chi + z902 * chi2 + z911 * eta * chi + z910 * eta + z920 * eta2;
236
237 if (extraParams!=NULL)
238 {
239 p->a1*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRDXi1(extraParams));
240 p->a2*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRDXi2(extraParams));
241 p->a3*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRDXi3(extraParams));
242 p->a4*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRDXi4(extraParams));
243 p->a5*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRDXi5(extraParams));
244 p->a6*=(1.0+XLALSimInspiralWaveformParamsLookupNonGRDXi6(extraParams));
245 }
246
247 /* Get the Spin of the final BH */
248 REAL8 s4 = -0.129;
249 REAL8 s5 = -0.384;
250 REAL8 t0 = -2.686;
251 REAL8 t2 = -3.454;
252 REAL8 t3 = 2.353;
253 REAL8 finspin = chi + s4*chi*etachi + s5*etachi*eta + t0*etachi + 2.*sqrt(3.)*eta +
254 t2*eta2 + t3*eta2*eta;
255
256 if( fabs(finspin) > 1.0 )
258
259 p->afin = finspin;
260
261 /* Get the Ringdown frequency */
262 REAL8 prefac = (1./(2.*LAL_PI)) * LAL_C_SI * LAL_C_SI * LAL_C_SI / (LAL_G_SI * M * LAL_MSUN_SI);
263 REAL8 k1 = 1.5251;
264 REAL8 k2 = -1.1568;
265 REAL8 k3 = 0.1292;
266 //printf("fRD prefac = %12.18f\n", prefac);
267
268 p->fRingDown = (prefac * (k1 + k2 * pow(1. - fabs(finspin), k3)));
269 p->MfRingDown = p->m_sec * p->fRingDown;
270
271 /* Get the quality factor of ring-fown, using Eq (5.6) of Main paper */
272 p->Qual = (0.7000 + (1.4187 * pow(1.0 - fabs(finspin), -0.4990)) );;
273
274 /* Get the transition frequencies, at which the model switches phase and
275 * amplitude prescriptions, as used in Eq.(5.9), (5.13) of the Main paper */
276 p->f1 = 0.1 * p->fRingDown;
277 p->f2 = p->fRingDown;
278 p->f0 = 0.98 * p->fRingDown;
279 p->d1 = 0.005;
280 p->d2 = 0.005;
281 p->d0 = 0.015;
282
283 /* Get the coefficients beta1, beta2, defined in Eq 5.7 of the main paper */
284 REAL8 Mfrd = p->MfRingDown;
285
286 p->b2 = ((-5./3.)* p->a1 * pow(Mfrd,(-8./3.)) - p->a2/(Mfrd*Mfrd) -
287 (p->a3/3.)*pow(Mfrd,(-4./3.)) + (2./3.)* p->a5 * pow(Mfrd,(-1./3.)) + p->a6)/eta;
288
289 REAL8 psiPMrd = IMRPhenomCGeneratePhasePM( p->fRingDown, eta, p );
290
291 p->b1 = psiPMrd - (p->b2 * Mfrd);
292
293 /* Taking the upper cut-off frequency as 0.15M */
294 //p->fCut = (1.7086 * eta * eta - 0.26592 * eta + 0.28236) / p->piM;
295 p->fCut = 0.15 / p->m_sec;
296
297 return p;
298
299}
300
301/*********************************************************************/
302/* The following function return the hyperbolic-Tan+ windows used */
303/* in Eq.(5.9), (5.13) of the Main paper */
304/*********************************************************************/
305static REAL8 wPlus( const REAL8 f, const REAL8 f0, const REAL8 d, const BBHPhenomCParams *params )
306{
307
308 REAL8 Mf = params->m_sec * f;
309 REAL8 Mf0 = params->m_sec * f0;
310
311 return ( 0.5 * (1. + tanh(4.*(Mf - Mf0)/d) ) );
312
313}
314
315/*********************************************************************/
316/* The following function return the hyperbolic-Tan- windows used */
317/* in Eq.(5.9), (5.13) of the Main paper */
318/*********************************************************************/
319static REAL8 wMinus( const REAL8 f, const REAL8 f0, const REAL8 d, const BBHPhenomCParams *params )
320{
321
322 REAL8 Mf = params->m_sec * f;
323 REAL8 Mf0 = params->m_sec * f0;
324
325 return ( 0.5 * (1. - tanh(4.*(Mf - Mf0)/d) ) );
326
327}
328
329/*********************************************************************/
330/* The following function return the closest higher power of 2 */
331/*********************************************************************/
332static size_t NextPow2_PC(const size_t n) {
333 return 1 << (size_t) ceil(log2(n));
334}
335
336/**
337 *
338 * main functions
339 *
340 */
341
342/***********************************************************************************/
343/* The following private function generates the Pre-Merger phase, as defined in */
344/* Eq. (5.3), (5.9) of the Main paper. */
345/***********************************************************************************/
346
348 REAL8 f, /**< frequency (Hz) */
349 REAL8 eta, /**< dimensionless mass-ratio */
350 const BBHPhenomCParams *params /**< pointer to Object storing coefficients and constants */
351 )
352{
353 REAL8 Mf = params->m_sec * f;
354 REAL8 w = pow( Mf, (1./3.) );
355 REAL8 w2 = w*w;
356 REAL8 w3 = w2*w;
357 REAL8 w5 = w3*w2;
358
359 REAL8 phasing = (params->a1/w5) + (params->a2/w3) + (params->a3/w) + params->a4 +
360 (params->a5*w2) +(params->a6*w3);
361 phasing /= eta;
362
363 return phasing;
364}
365
366
367/***********************************************************************************/
368/* The following private function generates the complete amplitude and phase of */
369/* PhenomC waveform, at a given frequency. */
370/* Eq. (5.3), (5.9) of the Main paper. */
371/***********************************************************************************/
372
374 REAL8 *amplitude, /**< pointer to memory for phenomC amplitude */
375 REAL8 *phasing, /**< pointer to memory for phenomC phase */
376 REAL8 f, /**< frequency (Hz) */
377 REAL8 eta, /**< dimensionless mass-ratio */
378 const BBHPhenomCParams *params /**< pointer to Object storing coefficients and constants */
379 )
380{
381 *amplitude = 0.0;
382 *phasing = 0.0;
383
384 /* Get the phase */
385 REAL8 v = cbrt(params->piM * f);
386 REAL8 Mf = params->m_sec * f;
387
388 if( v >= 1.0 )
390
391 REAL8 v2 = v*v;
392 REAL8 v3 = v*v2;
393 REAL8 v4 = v2*v2;
394 REAL8 v5 = v3*v2;
395 REAL8 v6 = v3*v3;
396 REAL8 v7 = v4*v3;
397 REAL8 v10 = v5*v5;
398
399 /* SPA part of the phase */
400 REAL8 phSPA = 1. + params->pfa1 * v + params->pfa2 * v2 + params->pfa3 * v3 + params->pfa4 * v4 +
401 (1. + log(v3)) * params->pfa5 * v5 + (params->pfa6 + params->pfa6log * log(v3))*v6 +
402 params->pfa7 * v7;
403 phSPA *= (params->pfaN / v5);
404
405 // Taking t0 = phi0 = 0
406 phSPA -= (LAL_PI / 4.);
407
408 REAL8 w = cbrt( Mf );
409 REAL8 w2 = w*w;
410 REAL8 w3 = w2*w;
411 REAL8 w5 = w3*w2;
412
413 /* The Pre-Merger (PM) phase */
414 REAL8 phPM = (params->a1/w5) + (params->a2/w3) + (params->a3/w) + params->a4 +
415 (params->a5*w2) +(params->a6*w3);
416 phPM /= eta;
417
418 /* Ring-down phase */
419 REAL8 phRD = params->b1 + params->b2 * params->m_sec * f;
420
421 REAL8 wPlusf1 = wPlus( f, params->f1, params->d1, params );
422 REAL8 wPlusf2 = wPlus( f, params->f2, params->d2, params );
423 REAL8 wMinusf1 = wMinus( f, params->f1, params->d1, params );
424 REAL8 wMinusf2 = wMinus( f, params->f2, params->d2, params );
425
426 *phasing = phSPA * wMinusf1 + phPM * wPlusf1 * wMinusf2 + phRD * wPlusf2;
427
428 /* Get the amplitude */
429 REAL8 xdot = 1. + params->xdota2*v2 + params->xdota3*v3 + params->xdota4*v4 +
430 params->xdota5*v5 + (params->xdota6 + params->xdota6log*log(v2))*v6 +
431 params->xdota7 * v7;
432 xdot *= (params->xdotaN * v10);
433
434 if( xdot < 0.0 && f < params->f1 )
435 {
436 XLALPrintError("omegaDot < 0, while frequency is below SPA-PM matching freq.");
438 }
439
440 REAL8 aPM = 0.0;
441
442 /* Following Emma's code, take only the absolute value of omegaDot, when
443 * computing the amplitude */
444 REAL8 omgdot = 1.5*v*xdot;
445 REAL8 ampfac = sqrt(fabs(LAL_PI/omgdot));
446
447 /* Get the real and imaginary part of the PM amplitude */
448 REAL8 AmpPMre = ampfac * params->AN * v2 * (1. + params->A2*v2 + params->A3*v3 +
449 params->A4*v4 + params->A5*v5 + (params->A6 + params->A6log*log(v2))*v6);
450 REAL8 AmpPMim = ampfac * params->AN * v2 * (params->A5imag * v5 + params->A6imag * v6);
451
452 /* Following Emma's code, we take the absolute part of the complex SPA
453 * amplitude, and keep that as the amplitude */
454 aPM = sqrt( AmpPMre * AmpPMre + AmpPMim * AmpPMim );
455
456 aPM += (params->g1 * pow(Mf,(5./6.)));
457
458 /* The Ring-down aamplitude */
459 REAL8 Mfrd = params->MfRingDown;
460
461 /* From Emma's code, use sigma = fRD * del2 / Qual */
462 REAL8 sig = Mfrd * params->del2 / params->Qual;
463 REAL8 sig2 = sig*sig;
464 REAL8 L = sig2 / ((Mf - Mfrd)*(Mf - Mfrd) + sig2/4.);
465
466 REAL8 aRD = params->del1 * L * pow(Mf, (-7./6.));
467
468 REAL8 wPlusf0 = wPlus( f, params->f0, params->d0, params );
469 REAL8 wMinusf0 = wMinus( f, params->f0, params->d0, params );
470
471 *amplitude = - (aPM * wMinusf0 + aRD * wPlusf0);
472
473 return XLAL_SUCCESS;
474}
475
476
477
478
479/*********************************************************************************/
480/* The following functions can be used to generate the individual components the */
481/* PhenomC waveform, both amplitude and phase. */
482/*********************************************************************************/
483
484//static size_t find_instant_freq(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 target, const size_t start);
485//static size_t find_peak_amp(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc);
486//static int apply_phase_shift(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 shift);
487//static int apply_inclination(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 inclination);
488
489/*
490static REAL8 IMRPhenomCGeneratePhaseSPA( REAL8 f, const BBHPhenomCParams *params );
491static REAL8 IMRPhenomCGeneratePhaseRD( REAL8 f, const BBHPhenomCParams *params );
492static REAL8 IMRPhenomCGenerateAmplitudePM( REAL8 f, const BBHPhenomCParams *params );
493static REAL8 IMRPhenomCGenerateAmplitudeRD( REAL8 f, const BBHPhenomCParams *params );
494*/
495
496/*
497static REAL8 IMRPhenomCGeneratePhaseSPA( REAL8 f, const BBHPhenomCParams *params )
498{
499 REAL8 v = cbrt(params->piM * f);
500
501 if( v >= 1.0 )
502 XLAL_ERROR(XLAL_EDOM);
503
504 REAL8 v2 = v*v;
505 REAL8 v3 = v*v2;
506 REAL8 v4 = v2*v2;
507 REAL8 v5 = v3*v2;
508 REAL8 v6 = v3*v3;
509 REAL8 v7 = v4*v3;
510
511 REAL8 phasing = 1. + params->pfa1 * v + params->pfa2 * v2 + params->pfa3 * v3 + params->pfa4 * v4 +
512 (1. + log(v3)) * params->pfa5 * v5 + (params->pfa6 + params->pfa6log * log(v3))*v6 +
513 params->pfa7 * v7;
514 phasing *= (params->pfaN / v5);
515
516 // Taking t0 = phi0 = 0
517 phasing -= (LAL_PI / 4.);
518
519 return phasing;
520}
521*/
522
523/*
524static REAL8 IMRPhenomCGeneratePhaseRD( REAL8 f, const BBHPhenomCParams *params )
525{
526 return (params->b1 + params->b2 * params->m_sec * f);
527}
528*/
529
530#if 0
531static REAL8 IMRPhenomCGenerateAmplitudePM( REAL8 f, const BBHPhenomCParams *params )
532{
533 REAL8 v = cbrt(params->piM * f);
534 REAL8 Mf = params->m_sec * f;
535
536 if( v >= 1.0 )
538
539 REAL8 v2 = v*v;
540 REAL8 v3 = v*v2;
541 REAL8 v4 = v2*v2;
542 REAL8 v5 = v3*v2;
543 REAL8 v6 = v3*v3;
544 REAL8 v7 = v4*v3;
545 REAL8 v10 = v5*v5;
546
547 REAL8 xdot = 1. + params->xdota2*v2 + params->xdota3*v3 + params->xdota4*v4 +
548 params->xdota5*v5 + (params->xdota6 + params->xdota6log*log(v2))*v6 +
549 params->xdota7 * v7;
550 xdot *= (params->xdotaN * v10);
551
552 if( xdot < 0.0 && f < params->fRingDown )
553 {
554 XLALPrintError("omegaDot < 0, while frequency is below RingDown freq.");
556 }
557
558 REAL8 amplitude = 0.0;
559 if( xdot > 0.0 )
560 {
561 REAL8 omgdot = 1.5*v*xdot;
562 REAL8 ampfac = sqrt(LAL_PI/omgdot);
563
564 /* Get the real and imaginary part of the PM amplitude */
565 REAL8 AmpPMre = ampfac * params->AN * v2 * (1. + params->A2*v2 + params->A3*v3 +
566 params->A4*v4 + params->A5*v5 + (params->A6 + params->A6log*log(v2))*v6);
567 REAL8 AmpPMim = ampfac * params->AN * v2 * (params->A5imag * v5 + params->A6imag * v6);
568
569 /* Following Emma's code, we take the absolute part of the complex SPA
570 * amplitude, and keep that as the amplitude */
571 amplitude = sqrt( AmpPMre * AmpPMre + AmpPMim * AmpPMim );
572 }
573
574 amplitude += (params->g1 * pow(Mf,(5./6.)));
575
576 return amplitude;
577}
578
579static REAL8 IMRPhenomCGenerateAmplitudeRD( REAL8 f, const BBHPhenomCParams *params )
580{
581 REAL8 Mf = params->m_sec * f;
582 REAL8 Mfrd = params->MfRingDown;
583
584 /* From Emma's code, use sigma = fRD * del2 / Qual */
585 REAL8 sig = Mfrd * params->del2 / params->Qual;
586 REAL8 sig2 = sig*sig;
587 REAL8 L = sig2 / ((Mf - Mfrd)*(Mf - Mfrd) + sig2/4.);
588
589 REAL8 amplitude = params->del1 * L * pow(Mf, (-7./6.));
590
591 return amplitude;
592}
593#endif
594
595#if 0
596/* To be put in the function IMRPhenomCGenerateFD
597 *
598 * */
599/* Get the phase */
600 /*
601 phSPA = IMRPhenomCGeneratePhaseSPA( f, params );
602 phPM = IMRPhenomCGeneratePhasePM( f, eta, params );
603 phRD = IMRPhenomCGeneratePhaseRD( f, params );
604
605 wPlusf1 = wPlus( f, params->f1, params->d1, params );
606 wPlusf2 = wPlus( f, params->f2, params->d2, params );
607 wMinusf1 = wMinus( f, params->f1, params->d1, params );
608 wMinusf2 = wMinus( f, params->f2, params->d2, params );
609
610 phPhenomC = phi0 + phSPA * wMinusf1 + phPM * wPlusf1 * wMinusf2 + phRD * wPlusf2;
611 */
612 /* Get the amplitude */
613 /*
614 aPM = IMRPhenomCGenerateAmplitudePM( f, params );
615 aRD = IMRPhenomCGenerateAmplitudeRD( f, params );
616
617 wPlusf0 = wPlus( f, params->f0, params->d0, params );
618 wMinusf0 = wMinus( f, params->f0, params->d0, params );
619
620 aPhenomC = aPM * wMinusf0 + aRD * wPlusf0;
621 */
622#endif
static size_t NextPow2_PC(const size_t n)
static int IMRPhenomCGenerateAmpPhase(REAL8 *amplitude, REAL8 *phasing, REAL8 f, REAL8 eta, const BBHPhenomCParams *params)
static REAL8 IMRPhenomCGeneratePhasePM(REAL8 f, REAL8 eta, const BBHPhenomCParams *params)
main functions
static BBHPhenomCParams * ComputeIMRPhenomCParamsSPA(const REAL8 m1, const REAL8 m2, const REAL8 chi, LALDict *extraParams)
Internal functions.
static REAL8 wPlus(const REAL8 f, const REAL8 f0, const REAL8 d, const BBHPhenomCParams *params)
static REAL8 wMinus(const REAL8 f, const REAL8 f0, const REAL8 d, const BBHPhenomCParams *params)
static BBHPhenomCParams * ComputeIMRPhenomCParams(const REAL8 m1, const REAL8 m2, const REAL8 chi, LALDict *extraParams)
static vector d(const double L_norm, const double J_norm, const vector roots)
Internal function that returns the coefficients "d_0", "d_2" and "d_4" from 1703.03967 corresponding ...
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi5(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi6(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDXi3(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi6L(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDXi5(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi1(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi4(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi3(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi7(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDXi1(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDXi2(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDXi6(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRPhi1(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDXi4(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi2(LALDict *params)
REAL8 M
Definition: bh_qnmode.c:133
const double w2
const double w
#define LAL_C_SI
#define LAL_MSUN_SI
#define LAL_PI
#define LAL_MTSUN_SI
#define LAL_GAMMA
#define LAL_G_SI
double REAL8
void * XLALMalloc(size_t n)
#define XLAL_ERROR_NULL(...)
#define XLAL_ERROR(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EDOM
p
Definition: burst.c:245