Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInspiral 5.0.3.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALInspiralParameterCalc.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2007 David Churches, Duncan Brown, Jolien Creighton, B.S. Sathyaprakash, Anand Sengupta, Thomas Cokelaer
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 Sathyaprakash, B. S.
22 * \file
23 * \ingroup LALInspiral_h
24 *
25 * \brief Given a pair of masses (or other equivalent parameters) compute
26 * related chirp parameters.
27 *
28 * ### Prototypes ###
29 *
30 * <tt>XLALInspiralParameterCalc()</tt>
31 * <ul>
32 * <li>\c params: Input/Output, given a pair of binary parameters and a lower
33 * frequency cutoff, other equivalent parameters are computed by this function.</li>
34 * </ul>
35 *
36 * ### Description ###
37 *
38 * The code takes as its input <tt>params->fLower</tt> in Hz and
39 * a pair of masses (in units of \f$M_\odot\f$) or chirptimes (in seconds measured from <tt>params->fLower</tt>)
40 * and computes all the other {\em mass} parameters in the \c params structure.
41 * Users choice of input pair of {\em masses} should be specified by appropriately setting
42 * the variable <tt>params->massChoice</tt> as described in the Table below:
43 *
44 * <table align="center" class="doxtable">
45 * <caption align="top" style="text-align: left; font-weight: normal;">
46 * Table I. For a given <tt>params->massChoice</tt> in column 1 the user should specify the
47 * parameters as in column 2, in units as in column 3. Column 4 gives the conventional meaning
48 * of the parameters. Chirp times are measured from a lower frequency cutoff given
49 * in <tt>params->fLower.</tt></caption>
50 * <tr><th><tt>params->massChoice</tt></th><th>User should set</th><th>in units</th><th>which means</th></tr>
51 * <tr><td>m1Andm2</td><td>(<tt>mass1, mass2</tt>)</td><td>\f$(M_\odot, M_\odot)\f$</td><td>\f$(m_1,m_2)\f$</td></tr>
52 * <tr><td>totalMassAndEta</td><td>(<tt>totalmass, eta</tt>)</td><td>\f$(M_\odot, 0 < \eta \le 1/4)\f$</td><td>\f$(m, \eta)\f$</td></tr>
53 * <tr><td>totalMassAndMu</td><td>(<tt>totalmass, mu</tt>)</td><td>\f$(M_\odot, M_\odot)\f$</td><td>\f$(m, \mu)\f$</td></tr>
54 * <tr><td>t02</td><td>(<tt>t0, t2</tt>)</td><td>(sec, sec)</td><td>\f$(\tau_0, \tau_2)\f$</td></tr>
55 * <tr><td>t03</td><td>(<tt>t0, t3</tt>)</td><td>(sec, sec)</td><td>\f$(\tau_0, \tau_3)\f$</td></tr>
56 * <tr><td>t04</td><td>(<tt>t0, t4</tt>)</td><td>(sec, sec)</td><td>\f$(\tau_0, \tau_4)\f$</td></tr>
57 * </table>
58 *
59 * If \c massChoice is not set properly an error condition will occur and
60 * the function is aborted with status code defined by
61 * #LALINSPIRALH_EMASSCHOICE in \ref LALInspiral.h.
62 * In the above list \f$m_{1}\f$ and \f$m_{2}\f$ are the masses of
63 * the two compact objects, \f$m=m_{1}+m_{2}\f$ is the total
64 * mass, \f$\eta=m_{1}m_{2}/(m_{1}+m_{2})^{2}\f$ is the
65 * symmetric mass ratio, \f$\mu=m_{1}m_{2}/(m_{1}+m_{2})\f$ is
66 * the reduced mass and \f$\tau\f$'s are the chirptimes
67 * defined in terms of \f$f_{a}\f$=\c fLower by:
68 * \f{eqnarray}{
69 * \tau_{0} = \frac{5}{256 \eta m^{5/3} (\pi f_{a})^{8/3}}, \ \ \
70 * \tau_{2} = \frac{(3715 + 4620 \eta)}{64512 \eta m (\pi f_{a})^{2}}, \ \ \
71 * \tau_{3} = \frac{\pi}{8 \eta m^{2/3} (\pi f_{a})^{5/3}} \\
72 * \tau_{4} = \frac{5}{128 \eta m^{1/3} (\pi f_{a})^{4/3}} \left[ \frac{3058673}{1016064} +
73 * \frac{5429}{1008} \eta + \frac{617}{144} \eta^{2} \right],\ \ \
74 * \tau_5 = \frac {5}{256\eta f_a} \left (\frac {7729}{252} + \eta \right ).
75 * \f}
76 * %% Beyond 2.5 PN order, chirp times do not have an
77 * %% explicit expression in terms of the masses and \f$f_a.\f$
78 * Whichever pair of parameters is given to the function as an input, the function
79 * calculates the rest. Apart from the various masses and chirptimes the function
80 * also calculates the chirp mass \f$\mathcal{M}=(\mu^{3} m^{2})^{1/5}\f$ and
81 * the total chirp time \f$\tau_C\f$ consistent with the approximation chosen:
82 *
83 * <table align="center" class="doxtable">
84 * <caption align="top" style="text-align: left; font-weight: normal;">
85 * Table II: \f$t_C\f$ will be set according to the PN order chosen in <tt>params->approximant.</tt>
86 * </caption>
87 * <tr><th></th><th>Newtonian</th><th>onePN</th><th>onePointFivePN</th><th>twoPN</th><th>twoPointFivePN</th></tr>
88 * <tr><td>\f$\tau_C\f$</td><td>\f$\tau_0\f$</td><td>\f$\tau_0 + \tau_2\f$</td><td>\f$\tau_0 + \tau_2-\tau_3\f$
89 * </td><td>\f$\tau_0 + \tau_2-\tau_3 + \tau_4\f$</td><td>\f$\tau_0 + \tau_2-\tau_3 + \tau_4 - \tau_5\f$</td></tr>
90 * </table>
91 *
92 * ### Algorithm ###
93 *
94 * Root finding by bisection method is used to solve for mass ratio \f$\eta\f$ when
95 * chirptimes \f$(\tau_0,\, \tau_2)\f$ or \f$(\tau_0,\, \tau_4)\f$ is input.
96 *
97 * ### Uses ###
98 *
99 * When appropriate this function calls:
100 * <code>
101 * XLALDBisectionFindRoot()
102 * XLALEtaTau02()
103 * XLALEtaTau04()
104 * </code>
105 *
106 * ### Notes ###
107 *
108 */
109
110
111
112#include <math.h>
113#include <lal/LALInspiral.h>
114#include <lal/FindRoot.h>
115
116void
120 )
121{
122 XLAL_PRINT_DEPRECATION_WARNING("XLALInspiralParameterCalc");
123
126
128 if (xlalErrno)
130
132 RETURN(status);
133}
134
135int
138 )
139{
140 REAL8 m1, m2, totalMass, eta, mu, piFl, etamin, tiny, ieta;
141 REAL8 x1, x2, A0, A2, A3, A4, B2, B4, C4,v,tN;
142 REAL8 theta = -11831.L/9240.L;
143 REAL8 lambda = -1987.L/3080.L;
144 static REAL8 oneby4;
145 void *pars;
146 REAL8 (*rootfunction)(REAL8, void *);
147 REAL8 xmin, xmax, xacc;
148 EtaTau02In Tau2In;
149 EtaTau04In Tau4In;
150
151 if (params == NULL)
153 if ((INT4)params->massChoice < 0)
155 if ((INT4)params->massChoice > 15)
157
158 totalMass = 0.0;
159 ieta = params->ieta;
160 ieta = 1.;
161 oneby4 = 1./4.;
162 etamin = 1.e-10;
163 tiny = 1.e-10;
164 piFl = LAL_PI * params->fLower;
165
166 switch(params->massChoice)
167 {
168 case massesAndSpin:
169 /*case spinOnly:*/
170 case minmaxTotalMass:
171 case m1Andm2:
172 case fixedMasses:
173
174 if (params->mass1 <= 0)
176 if (params->mass2 <= 0)
178
179 m1 = params->mass1;
180 m2 = params->mass2;
181 params->totalMass = totalMass = m1+m2;
182 params->eta = eta = m1*m2/(totalMass*totalMass);
183 if (params->eta > oneby4) {
184 params->eta -= tiny;
185 }
186 params->mu = mu = m1*m2/totalMass;
187 params->chirpMass = pow(mu,0.6)*pow(totalMass,0.4);
188 params->psi0 = 3./128./params->eta
189 * 1. * pow((LAL_PI * params->totalMass * LAL_MTSUN_SI),-5./3.) ;
190 params->psi3 = -3./128./params->eta
191 * (16 * LAL_PI) * pow((LAL_PI * params->totalMass * LAL_MTSUN_SI),-2./3.);
192
193 break;
194
195 case totalMassAndEta:
196 case totalMassUAndEta:
197
198 if (params->totalMass <= 0)
200 if (params->eta <= 0)
202
203 if (params->eta > oneby4) {
204 params->eta -= tiny;
205 }
206 if (params->eta > oneby4)
208
209 totalMass = params->totalMass;
210 eta = params->eta;
211 if (eta <= oneby4) {
212 params->mass1 = 0.5*totalMass * ( 1.L + sqrt(1.L - 4.L*eta));
213 params->mass2 = 0.5*totalMass * ( 1.L - sqrt(1.L - 4.L*eta));
214 }
215 params->mu = eta*totalMass;
216 params->chirpMass = pow(eta,0.6)*totalMass;
217
218 break;
219
220 case totalMassAndMu:
221
222 if (params->totalMass <= 0)
224 if (params->mu <= 0)
226 if (params->mu >= params->totalMass)
228
229 totalMass = params->totalMass;
230 mu = params->mu;
231 eta = (params->mu)/totalMass;
232 if (eta > oneby4) {
233 eta -= tiny;
234 }
235 params->eta = eta;
236 if (eta <= oneby4) {
237 params->mass1 = 0.5*totalMass * ( 1.L + sqrt(1.L - 4.L*eta));
238 params->mass2 = 0.5*totalMass * ( 1.L - sqrt(1.L - 4.L*eta));
239 }
240 params->chirpMass = pow(eta,0.6)*totalMass;
241 params->mu = eta*totalMass;
242
243 break;
244
245 case t02:
246
247 if (params->t0 <= 0)
249 if (params->t2 <= 0)
251
252 A0 = 5./ pow(piFl, (8./3.))/256.;
253 A2 = 3715.0/(64512.0*pow(piFl,2.0));
254 B2 = 4620.0/3715 * ieta;
255 Tau2In.t2 = params->t2;
256 Tau2In.A2 = A2 * pow(params->t0/A0, 0.6);
257 Tau2In.B2 = B2;
258
259 pars = (void *) &Tau2In;
260 rootfunction = &XLALEtaTau02;
261 xmax = oneby4+tiny;
262 xmin = etamin;
263 xacc = 1.e-8;
264 x1 = XLALEtaTau02(xmax, pars);
267 x2 = XLALEtaTau02(xmin, pars);
270
271 if (x1*x2 > 0) {
272 params->eta = 0.;
273 return XLAL_SUCCESS;
274 } else {
275 eta = XLALDBisectionFindRoot(rootfunction, xmin, xmax, xacc, pars);
276 if (XLAL_IS_REAL8_FAIL_NAN(eta))
278 }
279 if (eta > oneby4) {
280 eta-=tiny;
281 }
282 params->eta = eta;
283 totalMass = pow(A0/(eta*params->t0), 0.6);
284 totalMass = params->totalMass = totalMass/LAL_MTSUN_SI;
285 if (eta <= oneby4) {
286 params->mass1 = 0.5*totalMass * ( 1.L + sqrt(1.L - 4.L*eta));
287 params->mass2 = 0.5*totalMass * ( 1.L - sqrt(1.L - 4.L*eta));
288 }
289 params->chirpMass = pow(eta,0.6)*totalMass;
290 params->mu = eta*totalMass;
291
292 break;
293
294 case t03:
295
296 if (params->t0 <= 0)
298 if (params->t3 <= 0)
300
301 A0 = 5./ pow(piFl, (8./3.))/256.;
302 A3 = LAL_PI / pow(piFl, (5./3.))/8.;
303 totalMass = A0 * params->t3/(A3 * params->t0);
304 eta = A0/(params->t0 * pow(totalMass, (5./3.)));
305
306 if (eta > oneby4) {
307 eta-=tiny;
308 }
309 params->eta = eta;
310 totalMass = params->totalMass = totalMass/LAL_MTSUN_SI;
311 if (eta <= oneby4) {
312 params->mass1 = 0.5*totalMass * ( 1.L + sqrt(1.L - 4.L*eta));
313 params->mass2 = 0.5*totalMass * ( 1.L - sqrt(1.L - 4.L*eta));
314 }
315 params->chirpMass = pow(eta,0.6)*totalMass;
316 params->mu = eta*totalMass;
317
318 break;
319
320 case t04:
321
322 if (params->t0 <= 0)
324 if (params->t4 <= 0)
326
327 A0 = 5./(256. * pow(piFl, (8./3.)));
328 A4 = 5./(128.0 * pow(piFl,(4./3.))) * 3058673./1016064.;
329 B4 = 5429./1008 * 1016064./3058673. * ieta;
330 C4 = 617./144. * 1016064./3058673. * ieta;
331 Tau4In.t4 = params->t4;
332 Tau4In.A4 = A4 * pow(params->t0/A0, 0.2);
333 Tau4In.B4 = B4;
334 Tau4In.C4 = C4;
335
336 pars = (void *) &Tau4In;
337 rootfunction = &XLALEtaTau04;
338 xmax = oneby4+tiny;
339 xmin = etamin;
340 xacc = 1.e-8;
341 x1 = XLALEtaTau04(xmax, pars);
344 x2 = XLALEtaTau04(xmin, pars);
347
348 if (x1*x2 > 0) {
349 params->eta = 0.;
350 return XLAL_SUCCESS;
351 } else {
352 eta = XLALDBisectionFindRoot(rootfunction, xmin, xmax, xacc, pars);
353 if (XLAL_IS_REAL8_FAIL_NAN(eta))
355 }
356 if (eta > oneby4) {
357 eta-=tiny;
358 }
359 params->eta = eta;
360 totalMass = pow(A0/(eta*params->t0), 0.6);
361 totalMass = params->totalMass = totalMass/LAL_MTSUN_SI;
362 if (eta <= oneby4) {
363 params->mass1 = 0.5*totalMass * ( 1.L + sqrt(1.L - 4.L*eta));
364 params->mass2 = 0.5*totalMass * ( 1.L - sqrt(1.L - 4.L*eta));
365 }
366 params->chirpMass = pow(eta,0.6)*totalMass;
367 params->mu = eta*totalMass;
368
369 break;
370
371
372 case psi0Andpsi3:
373 if (params->psi0 > 0 && params->psi3 < 0)
374 {
375 params->totalMass = totalMass = -params->psi3/(16.L * LAL_PI * LAL_PI * params->psi0)/LAL_MTSUN_SI;
376 params->eta = eta = 3.L/(128.L * params->psi0 * pow (LAL_PI * totalMass*LAL_MTSUN_SI, (5./3.)));
377
378 /* if eta < 1/4 amd M > 0 then physical values*/
379 if (eta <= oneby4)
380 {
381 params->mass1 = 0.5*totalMass * ( 1.L + sqrt(1.L - 4.L*eta));
382 params->mass2 = 0.5*totalMass * ( 1.L - sqrt(1.L - 4.L*eta));
383 params->mu = eta*totalMass;
384 params->chirpMass = pow(eta,0.6)*totalMass;
385 }
386 }
387 else
388 {
389 params->eta = 0.;
390 return XLAL_SUCCESS;
391 }
392 break;
393
394 default:
395 XLALPrintError("XLAL Error - %s: Improper choice for massChoice\n", __func__);
397 break;
398 }
399
400 if (params->eta > oneby4) {
401 params->eta-=tiny;
402 }
403 totalMass = totalMass*LAL_MTSUN_SI;
404
405 /* Should use the coefficients from LALInspiraSetup.c to avoid errors.
406 * */
407 v = pow(piFl * totalMass, 1.L/3.L);
408 tN = 5.L/256.L / eta * totalMass / pow(v,8.L);
409
410 params->t0 = 5.0L/(256.0L*eta*pow(totalMass,(5./3.))*pow(piFl,(8./3.)));
411 params->t2 = (3715.0L + (4620.0L*ieta*eta))/(64512.0*eta*totalMass*pow(piFl,2.0));
412 params->t3 = LAL_PI/(8.0*eta*pow(totalMass,(2./3.))*pow(piFl,(5./3.)));
413 params->t4 = (5.0/(128.0*eta*pow(totalMass,(1./3.))*pow(piFl,(4./3.))))
414 * (3058673./1016064. + 5429.*ieta*eta/1008. +617.*ieta*eta*eta/144.);
415 params->t5 = -5.*(7729./252. - 13./3.*ieta*eta)/(256.*eta*params->fLower);
416 /* This is a ddraft. t6 and t7 need to be checked propely*/
417 params->t6 = -10052469856691./23471078400. + 128./3.*LAL_PI*LAL_PI
418 +(15335597827.L/15240960.L-451.L/12.L*LAL_PI*LAL_PI+352./3.*theta-2464.L/9.L*lambda)*ieta*eta
419 +6848.L/105.L* LAL_GAMMA
420 -15211.L/1728.L*ieta*eta*eta+25565.L/1296.L*eta*eta*eta*ieta;
421 params->t6 = tN * (params->t6 + 6848.L/105.L*log(4.*v)) * pow(v,6);
422 params->t7 = (-15419335.L/127008.L-75703.L/756.L*ieta*eta+14809.L/378.L*ieta*eta*eta) * LAL_PI * tN * pow(v,7);
423
424 params->psi0 = 3.L/(128.L * eta * pow(LAL_PI * totalMass, (5./3.)));
425 params->psi3 = -3.L * LAL_PI/(8.L * eta * pow(LAL_PI * totalMass, (2./3.)));
426
427 switch (params->order) {
428
430 case LAL_PNORDER_HALF:
431 params->t2=0.0;
432/* params->t3=0.0;*/
433 params->t4=0.0;
434 params->t5=0.0;
435 params->t6=0.0;
436 params->t7=0.0;
437 params->tC = params->t0;
438 break;
439
440 case LAL_PNORDER_ONE:
441 params->t3=0.0;
442 params->t4=0.0;
443 params->t5=0.0;
444 params->t6=0.0;
445 params->t7=0.0;
446 params->tC = params->t0 + params->t2;
447 break;
448
450 params->t4=0.0;
451 params->t5=0.0;
452 params->t6=0.0;
453 params->t7=0.0;
454 params->tC = params->t0 + params->t2 - params->t3;
455 break;
456
457 case LAL_PNORDER_TWO:
458 params->t5=0.0;
459 params->t6=0.0;
460 params->t7=0.0;
461 params->tC = params->t0 + params->t2 - params->t3 + params->t4;
462 break;
463
465 params->t6 = 0.0;
466 params->t7 = 0.0;
467 params->tC = params->t0 + params->t2 - params->t3 + params->t4 - params->t5;
468 break;
469
471 /*check the initialisation and then comment the next line. For now we
472 * set t6=0*/
473 params->t6 = 0;
474 params->t7 = 0.0;
475 params->tC = params->t0 + params->t2 - params->t3 + params->t4 - params->t5 + params->t6;
476 break;
477
479 default:
480 /*check the initialisation and then comment the next line. For now we
481 * set t6=0 and t7=0*/
482 params->t6 = 0;
483 params->t7 = 0.0;
484 params->tC = params->t0 + params->t2 - params->t3 + params->t4 - params->t5 + params->t6 - params->t7;
485 break;
486 }
487
488 return XLAL_SUCCESS;
489
490}
REAL8 XLALEtaTau02(REAL8 eta, void *p)
Definition: LALEtaTau02.c:88
REAL8 XLALEtaTau04(REAL8 eta, void *p)
Definition: LALEtaTau04.c:89
int XLALInspiralParameterCalc(InspiralTemplate *params)
void LALInspiralParameterCalc(LALStatus *status, InspiralTemplate *params)
#define tiny
#define ATTATCHSTATUSPTR(statusptr)
#define DETATCHSTATUSPTR(statusptr)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
#define ABORTXLAL(sp)
double theta
REAL8 XLALDBisectionFindRoot(REAL8(*y)(REAL8, void *), REAL8 xmin, REAL8 xmax, REAL8 xacc, void *params)
#define LAL_PI
#define LAL_MTSUN_SI
#define LAL_GAMMA
double REAL8
int32_t INT4
@ massesAndSpin
UNDOCUMENTED.
Definition: LALInspiral.h:193
@ totalMassAndMu
total mass and reduced mass
Definition: LALInspiral.h:182
@ fixedMasses
The two masses are given by the input parameter structure.
Definition: LALInspiral.h:190
@ totalMassAndEta
total mass and symmetric mass ratio
Definition: LALInspiral.h:180
@ psi0Andpsi3
BCV parameters and .
Definition: LALInspiral.h:187
@ t03
chirptimes and , and
Definition: LALInspiral.h:185
@ t02
chirptimes and
Definition: LALInspiral.h:184
@ totalMassUAndEta
total mass and eta but uniform distribution in totalMass
Definition: LALInspiral.h:181
@ m1Andm2
component masses
Definition: LALInspiral.h:179
@ minmaxTotalMass
UNDOCUMENTED.
Definition: LALInspiral.h:194
@ t04
chirptimes and
Definition: LALInspiral.h:186
LAL_PNORDER_TWO_POINT_FIVE
LAL_PNORDER_THREE
LAL_PNORDER_TWO
LAL_PNORDER_ONE
LAL_PNORDER_THREE_POINT_FIVE
LAL_PNORDER_HALF
LAL_PNORDER_ONE_POINT_FIVE
LAL_PNORDER_NEWTONIAN
#define xlalErrno
#define XLAL_ERROR(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_PRINT_DEPRECATION_WARNING(replacement)
#define XLAL_IS_REAL8_FAIL_NAN(val)
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFUNC
XLAL_EDOM
XLAL_EINVAL
list mu
These are the input structures needed to solve for the mass ratio given the chirptimes or .
Definition: LALInspiral.h:150
The inspiral waveform parameter structure containing information about the waveform to be generated.
Definition: LALInspiral.h:205