Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-5e288d3
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALSimInspiralHGimri.c
Go to the documentation of this file.
1
2//====================================================================================================
3// Code to generate a quasi-circular intermediate mass-ratio inspiral trajectory and graviational
4// waveform following Huerta & Gair 2011 (1009.1985). The system's evolution is divided into four
5// stages. See Huerta & Gair 2011 for discussion of each.
6//
7// 1. Adiabatic inspiral
8// 2. Transition (see also Ori & Thorne 2000 -- gr-qc/0003032)
9// 3. Geodesic Plunge
10// 4. Quasi-Normal Ringdown
11//
12// System is evolved using 2PN angular momentum flux from Gair and Glampedakis 2006 (gr-qc/0510129)
13// with higher order fits to Teukolsky results. Azimuthal phase is evolved using the 2PN conservative
14// corrections discussed in Huerta & Gair 2008 (0812.4208), and the waveform is generated with the
15// flat space quadrupole formula (see Babak et al 2007 -- gr-qc/0607007).
16//
17// Tom Callister
18// thomas.a.callister@gmail.com
19// 2014
20//
21//====================================================================================================
22
23#include <math.h>
24#include <gsl/gsl_vector.h>
25#include <gsl/gsl_spline.h>
26#include <gsl/gsl_multiroots.h>
27#include <gsl/gsl_roots.h>
28
29#include <lal/LALConstants.h>
30#include <lal/LALDatatypes.h>
31#include <lal/TimeSeries.h>
32#include <lal/Units.h>
33#include <lal/Date.h>
34
35#define c ((REAL8)(LAL_C_SI))
36#define G ((REAL8)(LAL_G_SI))
37#define pi ((REAL8)(LAL_PI))
38#define pc ((REAL8)(LAL_PC_SI))
39#define Msun ((REAL8)(LAL_MSUN_SI))
40#define Msun_sec ((REAL8)(G*Msun/(c*c*c)))
41#define GPC_sec ((REAL8)(pc*pow(10.,9.)/c))
42
43//Function prototypes
44static REAL8 XLALHGimri_AngMomFlux(REAL8 q, REAL8 r, REAL8 nu);
45static INT4 XLALHGimri_PlusEquations(const gsl_vector * x, void *params, gsl_vector * f);
46static INT4 XLALHGimri_CrossEquations(const gsl_vector * z, void *params, gsl_vector * g);
47static REAL8 XLALHGimri_dFdr(REAL8 E, REAL8 Lz, REAL8 a, REAL8 r);
48static REAL8 XLALHGimri_initialP(REAL8 p0, void *params);
49static REAL8 XLALHGimri_dLzdr(REAL8 q, REAL8 r);
50static INT4 HGimri_start(REAL8 m, REAL8 M, REAL8 q, REAL8 D, REAL8 Sdotn, REAL8 phi0, REAL8 p0, REAL8Sequence *hplus, REAL8Sequence *hcross, REAL8 dt, UINT4 Npts);
51
52INT4 XLALHGimriGenerator(REAL8TimeSeries **hplus,REAL8TimeSeries **hcross,REAL8 phi0,REAL8 dt,REAL8 m1,REAL8 m2,REAL8 f_min,REAL8 r,REAL8 inc,REAL8 s1z);
53
54//Data type to hold current simulation regime
55enum stage {INSPIRAL, TRANSITION, PLUNGE, FINALPLUNGE};
56
57//Structure template used to match plunging waveform onto QNM ringdown at three points across the light-ring (LR)
58struct rparams {
59 REAL8 M; //Total binary mass (units seconds)
60 REAL8 Sdotn; //Cosine of inclination angle (relative to detector)
61 REAL8 final_mass; //Final BH mass (units seconds)
62 REAL8 final_q; //Final BH reduced spin
63 REAL8 w0, w1, w2; //Frequencies of n=0,1,2 overtones of the (l,m)=(2,+/-2) QNMs, real parts
64 REAL8 wi0, wi1, wi2; //Frequencies of n=0,1,2 overtones of the (l,m)=(2,+/-2) QNMs, imaginary parts
65 REAL8 ab_factor; //Amplitude factor
66 REAL8 hp_1, dhpdi_1; //hplus and its derivative at point 1
67 REAL8 hp_2, dhpdi_2; //hplus and its derivative at point 2
68 REAL8 hp_3, dhpdi_3; //hplus and its derivative at point 3
69 REAL8 hx_1, dhxdi_1; //hcross and its derivative at point 1
70 REAL8 hx_2, dhxdi_2; //hcross and its derivative at point 2
71 REAL8 hx_3, dhxdi_3; //hcross and its derivative at point 3
72 REAL8 dt; //Dimensionless time step
73 };
74
75//Structure template used in computing in initial radius
76struct pParams {
77 REAL8 m; //CO mass (seconds)
78 REAL8 M; //BH mass (seconds)
79 REAL8 q; //BH reduced spin
80 REAL8 f_min; //Desired initial GW frequency
81 };
82
83static REAL8 XLALHGimri_initialP(REAL8 p0, void *params) {
84
85 //====================================================
86 // Root-finding function used to find initial radius p
87 // corresponding to a GW frequency f_min.
88 //====================================================
89
90 //Read params
91 struct pParams *p = (struct pParams *) params;
92 REAL8 m = p->m;
93 REAL8 M = p->M;
94 REAL8 q = p->q;
95 REAL8 f_min = p->f_min;
96 REAL8 nu = (m*M)/pow(m+M,2.);
97
98 //Conservative corrections
99 REAL8 d0 = 1./8.;
100 REAL8 d1 = 1975./896.;
101 REAL8 d1_5 = -27.*pi/10.;
102 REAL8 l1_5 = -191./160.;
103 REAL8 d2 = 1152343./451584.;
104
105 //Match orbital frequency to twice the GW frequency f_min
106 return (1./(pow(p0,3/2.)+q))*(1. + nu*(d0+d1/p0+(d1_5+q*l1_5)*pow(p0,-1.5)+d2*pow(p0,-2.0))) - pi*(f_min)*(m+M)*Msun_sec;
107
108 }
109
110
111static REAL8 XLALHGimri_AngMomFlux(REAL8 q, REAL8 r, REAL8 nu) {
112
113 //===================================================================
114 // Angular momentum flux from Gair & Glampedakis 2006 (gr-qc/0510129).
115 // See Eq (45) for 2PN fluxes, and (57) for the higher order fits to
116 // circular Teukolsky data.
117 //===================================================================
118
119 REAL8 pref = -6.4/(pow(r,7./2.));
120 REAL8 rtr=sqrt(r);
121
122 //Cosine of orbital inclination. 1 if prograde, -1 if retrograde
123 REAL8 cosInc;
124 REAL8 sign;
125
126 if (q >= 0.0) {
127 cosInc = 1.;
128 sign = 1.;
129 }
130 else {
131 cosInc = -1.;
132 sign = -1.;
133 }
134
135 //Absolute value of spin
136 REAL8 q_abs = fabs(q);
137
138 //Higher order Teukolsky terms
139 REAL8 c1a,c1b,c1c;
140 REAL8 c2a,c2b,c2c;
141 REAL8 c3a,c3b,c3c;
142 REAL8 c4a,c4b,c4c;
143 REAL8 c5a,c5b,c5c;
144 REAL8 c6a,c6b,c6c;
145 REAL8 c7a,c7b,c7c;
146 REAL8 c8a,c8b,c8c;
147 REAL8 c9a,c9b,c9c;
148 REAL8 c10a,c10b,c10c;
149 REAL8 c11a,c11b,c11c;
150 REAL8 higherorderfit,c1,c2,correction,circbit;
151
152 c1a=-10.741956;
153 c1b=28.5942157;
154 c1c=-9.077378144;
155 c1=c1a+(c1b+c1c/rtr)/rtr;
156
157 c2a=-1.428362761;
158 c2b=10.70029768;
159 c2c=-33.70903016;
160 c2=(c2a+(c2b+c2c/rtr)/rtr);
161
162 c3a=-28.15174147;
163 c3b=60.9607071973;
164 c3c=40.99984205;
165
166 c4a=-0.348161211;
167 c4b=2.37258476;
168 c4c=-66.65840948;
169
170 c5a=-0.715392387;
171 c5b=3.21592568;
172 c5c=5.28887649;
173
174 c6a=-7.6103411;
175 c6b=128.87778309;
176 c6c=-475.4650442;
177
178 c7a=12.290783385;
179 c7b=-113.1250548;
180 c7c=306.11883292;
181
182 c8a=40.9258725;
183 c8b=-347.2713496;
184 c8c=886.50332051;
185
186 c9a=-25.48313727;
187 c9b=224.22721861;
188 c9c=-490.98212316;
189
190 c10a=-9.006337706;
191 c10b=91.17666278;
192 c10c=-297.001939215;
193
194 c11a=-0.64500047;
195 c11b=-5.13591989;
196 c11c=47.19818628;
197
198 correction = q_abs*(
199 736.2086781-283.9553066*rtr+q_abs*(-1325.1852209+483.266206498*rtr)+pow(q_abs,2.0)*(634.49936445-219.223848944*rtr)
200 +(82.07804475-25.82025864*rtr)+q_abs*(-904.16109275+301.477789146*rtr)+pow(q_abs,2.0)*(827.31891826-271.9659423*rtr)
201 )/(r*rtr);
202
203 higherorderfit = q_abs*c1 + pow(q_abs,3.)*c2
204 + cosInc*( c3a+(c3b+c3c/rtr)/rtr )
205 + pow(q_abs,2.)*cosInc*(c4a+(c4b+c4c/rtr)/rtr)
206 + pow(q_abs,4.)*cosInc*(c5a+(c5b+c5c/rtr)/rtr)
207 + q_abs*(c6a+(c6b+c6c/rtr)/rtr)
208 + pow(q_abs,3.)*(c7a+(c7b+c7c/rtr)/rtr)
209 + pow(q_abs,2.)*cosInc*(c8a+(c8b+c8c/rtr)/rtr)
210 + pow(q_abs,4.)*cosInc*(c9a+(c9b+c9c/rtr)/rtr)
211 + pow(q_abs,3.)*(c10a+(c10b+c10c/rtr)/rtr)
212 + pow(q_abs,4.)*cosInc*(c11a+(c11b+c11c/rtr)/rtr)
213 + correction;
214
215 //Bracketed flux terms from Eq (45) with higher order corrections
216 circbit = cosInc - (61./12.)*q_abs/(r*rtr) -1247.*cosInc/(336.*r) + 4.*pi*cosInc/(r*rtr) - 44711.*cosInc/(9072.*r*r)
217 + (33./16.)*q_abs*q_abs*cosInc/(r*r) + higherorderfit/(r*r*rtr);
218
219 //Return full flux. (sign) factor needed to guarantee that Ldot is negative at leading order,
220 //following our adopted spin convention.
221 return(sign*nu*pref*circbit);
222
223 }
224
225static INT4 XLALHGimri_PlusEquations(const gsl_vector * x, void *params, gsl_vector * f) {
226
227 //=============================================================================
228 // Used by gsl_multiroots to match plus-polarized plunge waveform onto ringdown.
229 // Ringdown waveform consists of (l,m)=(2,+/-2) QNMs, with n=0,1,2 overtones. The
230 // RD waveform and its time derivative is matched to plunge waveform at three
231 // points across the light ring, considering the n=0 fundamental tone at the first
232 // time step, the n=0,1 tones at the second, and the n=0,1,2 tones at the third
233 //=============================================================================
234
235 //See definition rparams above for structure information
236 REAL8 dt = ((struct rparams *) params)->dt;
237 REAL8 wi0 = ((struct rparams *) params)->wi0;
238 REAL8 w0 = ((struct rparams *) params)->w0;
239 REAL8 wi1 = ((struct rparams *) params)->wi1;
240 REAL8 w1 = ((struct rparams *) params)->w1;
241 REAL8 wi2 = ((struct rparams *) params)->wi2;
242 REAL8 w2 = ((struct rparams *) params)->w2;
243 REAL8 final_mass = ((struct rparams *) params)->final_mass;
244 REAL8 final_q = ((struct rparams *) params)->final_q;
245 REAL8 Sdotn = ((struct rparams *) params)->Sdotn;
246 REAL8 M = ((struct rparams *) params)->M;
247 REAL8 ab_factor = ((struct rparams *) params)->ab_factor;
248
249 REAL8 hp_1 = ((struct rparams *) params)->hp_1;
250 REAL8 dhpdi_1 = ((struct rparams *) params)->dhpdi_1;
251 REAL8 hp_2 = ((struct rparams *) params)->hp_2;
252 REAL8 dhpdi_2 = ((struct rparams *) params)->dhpdi_2;
253 REAL8 hp_3 = ((struct rparams *) params)->hp_3;
254 REAL8 dhpdi_3 = ((struct rparams *) params)->dhpdi_3;
255
256 //Leading amplitude factors
257 REAL8 aone = 1.+ 2.*Sdotn + Sdotn*Sdotn;
258 REAL8 atwo = 2.+ Sdotn - 4.*Sdotn*Sdotn - 3.*Sdotn*Sdotn*Sdotn;
259
260 //Unknown RD constants
261 const REAL8 a0n = gsl_vector_get (x, 0);
262 const REAL8 a0p = gsl_vector_get (x, 1);
263 const REAL8 a1n = gsl_vector_get (x, 2);
264 const REAL8 a1p = gsl_vector_get (x, 3);
265 const REAL8 a2n = gsl_vector_get (x, 4);
266 const REAL8 a2p = gsl_vector_get (x, 5);
267
268 //Functions of the form: (Ringdown Quantity) - (Plunge Quantity). We seek the roots of these functions
269
270 //h_plus at time step 1
271 const REAL8 yy0 = ab_factor*a0n*(aone-(2./9.)*atwo*w0*final_q) - hp_1;
272
273 //h_plus time derivative at time step 1
274 const REAL8 yy1 = ab_factor*(aone-(2./9.)*atwo*w0*final_q)*(-a0p*w0 - a0n*wi0) - dhpdi_1*final_mass/(M*dt);
275
276 //h_plus at time step 2
277 const REAL8 yy2 = (((aone - (2./9.)*atwo*w0*final_q)*(a0n*cos((M*dt*w0)/final_mass) - 1.*a0p*sin((M*dt*w0)/final_mass)))/exp((M*dt*wi0)/final_mass)
278 + ((aone - (2./9.)*atwo*w1*final_q)*(a1n*cos((M*dt*w1)/final_mass) - 1.*a1p*sin((M*dt*w1)/final_mass)))/exp((M*dt*wi1)/final_mass))
279 - hp_2/ab_factor;
280
281 //h_plus time derivative at time step 2
282 const REAL8 yy3 = (((-1.*(aone - (2./9.)*atwo*w0*final_q)*wi0*(a0n*cos((M*dt*w0)/final_mass) - 1.*a0p*sin((M*dt*w0)/final_mass)))/exp((M*dt*wi0)/final_mass)
283 + ((aone - (2./9.)*atwo*w0*final_q)*(-1.*a0p*w0*cos((M*dt*w0)/final_mass) - 1.*a0n*w0*sin((M*dt*w0)/final_mass)))/exp((M*dt*wi0)/final_mass)
284 - (1.*(aone - (2./9.)*atwo*w1*final_q)*wi1*(a1n*cos((M*dt*w1)/final_mass) - 1.*a1p*sin((M*dt*w1)/final_mass)))/exp((M*dt*wi1)/final_mass)
285 + ((aone - (2./9.)*atwo*w1*final_q)*(-1.*a1p*w1*cos((M*dt*w1)/final_mass) - 1.*a1n*w1*sin((M*dt*w1)/final_mass)))/exp((M*dt*wi1)/final_mass)))
286 - dhpdi_2*final_mass/(M*dt*ab_factor);
287
288 //h_plus at time step 3
289 const REAL8 yy4 = (((aone - (2./9.)*atwo*w0*final_q)*(a0n*cos((2.*M*dt*w0)/final_mass) - 1.*a0p*sin((2.*M*dt*w0)/final_mass)))/exp((2.*M*dt*wi0)/final_mass)
290 + ((aone - (2./9.)*atwo*w1*final_q)*(a1n*cos((2.*M*dt*w1)/final_mass) - 1.*a1p*sin((2.*M*dt*w1)/final_mass)))/exp((2.*M*dt*wi1)/final_mass)
291 + ((aone - (2./9.)*atwo*w2*final_q)*(a2n*cos((2.*M*dt*w2)/final_mass) - 1.*a2p*sin((2.*M*dt*w2)/final_mass)))/exp((2.*M*dt*wi2)/final_mass))
292 - hp_3/ab_factor;
293
294 //h_plus time derivative at time step 3
295 const REAL8 yy5 = (((-1.*(aone - (2./9.)*atwo*w0*final_q)*wi0*(a0n*cos((2.*M*dt*w0)/final_mass) - 1.*a0p*sin((2.*M*dt*w0)/final_mass)))/exp((2.*M*dt*wi0)/final_mass)
296 + ((aone - (2./9.)*atwo*w0*final_q)*(-1.*a0p*w0*cos((2.*M*dt*w0)/final_mass) - 1.*a0n*w0*sin((2.*M*dt*w0)/final_mass)))/exp((2.*M*dt*wi0)/final_mass)
297 - (1.*(aone - (2./9.)*atwo*w1*final_q)*wi1*(a1n*cos((2.*M*dt*w1)/final_mass) - 1.*a1p*sin((2.*M*dt*w1)/final_mass)))/exp((2.*M*dt*wi1)/final_mass)
298 + ((aone - (2./9.)*atwo*w1*final_q)*(-1.*a1p*w1*cos((2.*M*dt*w1)/final_mass) - 1.*a1n*w1*sin((2.*M*dt*w1)/final_mass)))/exp((2.*M*dt*wi1)/final_mass)
299 - (1.*(aone - (2./9.)*atwo*w2*final_q)*wi2*(a2n*cos((2.*M*dt*w2)/final_mass) - 1.*a2p*sin((2.*M*dt*w2)/final_mass)))/exp((2.*M*dt*wi2)/final_mass)
300 + ((aone - (2./9.)*atwo*w2*final_q)*(-1.*a2p*w2*cos((2.*M*dt*w2)/final_mass) - 1.*a2n*w2*sin((2.*M*dt*w2)/final_mass)))/exp((2.*M*dt*wi2)/final_mass)))
301 - dhpdi_3*final_mass/(M*dt*ab_factor);
302
303 gsl_vector_set (f, 0, yy0);
304 gsl_vector_set (f, 1, yy1);
305 gsl_vector_set (f, 2, yy2);
306 gsl_vector_set (f, 3, yy3);
307 gsl_vector_set (f, 4, yy4);
308 gsl_vector_set (f, 5, yy5);
309
310 return GSL_SUCCESS;
311 }
312
313
314static INT4 XLALHGimri_CrossEquations(const gsl_vector * z, void *params, gsl_vector * g) {
315
316 //=============================================================================
317 // Used by gsl_multiroots to match cross-polarized plunge waveform onto ringdown.
318 // Ringdown waveform consists of (l,m)=(2,+/-2) QNMs, with n=0,1,2 overtones. The
319 // RD waveform and its time derivative is matched to plunge waveform at three
320 // points across the light ring, considering the n=0 fundamental tone at the first
321 // time step, the n=0,1 tones at the second, and the n=0,1,2 tones at the third
322 //=============================================================================
323
324 //See definition rparams above for structure information
325 REAL8 dt = ((struct rparams *) params)->dt;
326 REAL8 wi0 = ((struct rparams *) params)->wi0;
327 REAL8 w0 = ((struct rparams *) params)->w0;
328 REAL8 wi1 = ((struct rparams *) params)->wi1;
329 REAL8 w1 = ((struct rparams *) params)->w1;
330 REAL8 wi2 = ((struct rparams *) params)->wi2;
331 REAL8 w2 = ((struct rparams *) params)->w2;
332 REAL8 final_mass = ((struct rparams *) params)->final_mass;
333 REAL8 final_q = ((struct rparams *) params)->final_q;
334 REAL8 Sdotn = ((struct rparams *) params)->Sdotn;
335 REAL8 M = ((struct rparams *) params)->M;
336 REAL8 ab_factor = ((struct rparams *) params)->ab_factor;
337 REAL8 hx_1 = ((struct rparams *) params)->hx_1;
338 REAL8 dhxdi_1 = ((struct rparams *) params)->dhxdi_1;
339 REAL8 hx_2 = ((struct rparams *) params)->hx_2;
340 REAL8 dhxdi_2 = ((struct rparams *) params)->dhxdi_2;
341 REAL8 hx_3 = ((struct rparams *) params)->hx_3;
342 REAL8 dhxdi_3 = ((struct rparams *) params)->dhxdi_3;
343
344 //Leading amplitude factors
345 REAL8 aone = 1.+ 2.*Sdotn + Sdotn*Sdotn;
346 REAL8 atwo = 2.+ Sdotn - 4.*Sdotn*Sdotn - 3.*Sdotn*Sdotn*Sdotn;
347
348 //Unknown RD constants
349 const REAL8 a0c = gsl_vector_get (z, 0);
350 const REAL8 a0cp = gsl_vector_get (z, 1);
351 const REAL8 a1c = gsl_vector_get (z, 2);
352 const REAL8 a1cp = gsl_vector_get (z, 3);
353 const REAL8 a2c = gsl_vector_get (z, 4);
354 const REAL8 a2cp = gsl_vector_get (z, 5);
355
356 //Functions of the form: (Ringdown Quantity) - (Plunge Quantity). We seek the roots of these functions
357
358 //h_cross at time step 1
359 const REAL8 yy0 = ab_factor*a0c*(aone-(2./9.)*atwo*w0*final_q) + hx_1;
360
361 //h_cross time derivative at time step 1
362 const REAL8 yy1 = ab_factor*(aone-(2./9.)*atwo*w0*final_q)*(a0cp*w0 - a0c*wi0) + dhxdi_1*final_mass/(M*dt);
363
364 //h_cross at time step 2
365 const REAL8 yy2 = (((aone - (2./9.)*atwo*w0*final_q)*(a0c*cos((M*dt*w0)/final_mass) + a0cp*sin((M*dt*w0)/final_mass)))/exp((M*dt*wi0)/final_mass)
366 + ((aone - (2./9.)*atwo*w1*final_q)*(a1c*cos((M*dt*w1)/final_mass) + a1cp*sin((M*dt*w1)/final_mass)))/exp((M*dt*wi1)/final_mass))
367 + hx_2/ab_factor;
368
369 //h_cross time derivative at time step 2
370 const REAL8 yy3 = (((-1.*(aone - (2./9.)*atwo*w0*final_q)*wi0*(a0c*cos((M*dt*w0)/final_mass) + a0cp*sin((M*dt*w0)/final_mass)))/exp((M*dt*wi0)/final_mass)
371 + ((aone - (2./9.)*atwo*w0*final_q)*(a0cp*w0*cos((M*dt*w0)/final_mass) - 1.*a0c*w0*sin((M*dt*w0)/final_mass)))/exp((M*dt*wi0)/final_mass)
372 - (1.*(aone - (2./9.)*atwo*w1*final_q)*wi1*(a1c*cos((M*dt*w1)/final_mass) + a1cp*sin((M*dt*w1)/final_mass)))/exp((M*dt*wi1)/final_mass)
373 + ((aone - (2./9.)*atwo*w1*final_q)*(a1cp*w1*cos((M*dt*w1)/final_mass) - 1.*a1c*w1*sin((M*dt*w1)/final_mass)))/exp((M*dt*wi1)/final_mass)))
374 + dhxdi_2*final_mass/(M*dt*ab_factor);
375
376 //h_cross at time step 3
377 const REAL8 yy4 = (((aone - (2./9.)*atwo*w0*final_q)*(a0c*cos((2.*M*dt*w0)/final_mass) + a0cp*sin((2.*M*dt*w0)/final_mass)))/exp((2.*M*dt*wi0)/final_mass)
378 + ((aone - (2./9.)*atwo*w1*final_q)*(a1c*cos((2.*M*dt*w1)/final_mass) + a1cp*sin((2.*M*dt*w1)/final_mass)))/exp((2.*M*dt*wi1)/final_mass)
379 + ((aone - (2./9.)*atwo*w2*final_q)*(a2c*cos((2.*M*dt*w2)/final_mass) + a2cp*sin((2.*M*dt*w2)/final_mass)))/exp((2.*M*dt*wi2)/final_mass))
380 + hx_3/ab_factor;
381
382 //h_cross time derivative at time step 3
383 const REAL8 yy5 = (((-1.*(aone - (2./9.)*atwo*w0*final_q)*wi0*(a0c*cos((2.*M*dt*w0)/final_mass) + a0cp*sin((2.*M*dt*w0)/final_mass)))/exp((2.*M*dt*wi0)/final_mass)
384 + ((aone - (2./9.)*atwo*w0*final_q)*(a0cp*w0*cos((2.*M*dt*w0)/final_mass) - 1.*a0c*w0*sin((2.*M*dt*w0)/final_mass)))/exp((2.*M*dt*wi0)/final_mass)
385 - (1.*(aone - (2./9.)*atwo*w1*final_q)*wi1*(a1c*cos((2.*M*dt*w1)/final_mass) + a1cp*sin((2.*M*dt*w1)/final_mass)))/exp((2.*M*dt*wi1)/final_mass)
386 + ((aone - (2./9.)*atwo*w1*final_q)*(a1cp*w1*cos((2.*M*dt*w1)/final_mass) - 1.*a1c*w1*sin((2.*M*dt*w1)/final_mass)))/exp((2.*M*dt*wi1)/final_mass)
387 - (1.*(aone - (2./9.)*atwo*w2*final_q)*wi2*(a2c*cos((2.*M*dt*w2)/final_mass) + a2cp*sin((2.*M*dt*w2)/final_mass)))/exp((2.*M*dt*wi2)/final_mass)
388 + ((aone - (2./9.)*atwo*w2*final_q)*(a2cp*w2*cos((2.*M*dt*w2)/final_mass) - 1.*a2c*w2*sin((2.*M*dt*w2)/final_mass)))/exp((2.*M*dt*wi2)/final_mass)))
389 + dhxdi_3*final_mass/(M*dt*ab_factor);
390
391 gsl_vector_set (g, 0, yy0);
392 gsl_vector_set (g, 1, yy1);
393 gsl_vector_set (g, 2, yy2);
394 gsl_vector_set (g, 3, yy3);
395 gsl_vector_set (g, 4, yy4);
396 gsl_vector_set (g, 5, yy5);
397
398 return GSL_SUCCESS;
399 }
400
401static REAL8 XLALHGimri_dLzdr(REAL8 q, REAL8 r) {
402
403 //==============================================
404 //Partial derivative (\partial L_z)/(\partial r)
405 //for circular geodesic orbits
406 //==============================================
407
408 //REAL8 sign;
409 //if (q >= 0.) sign = 1.;
410 //else sign = -1.;
411
412 REAL8 denom1 = r*r*r-3.*r*r+2.*q*pow(r,3./2.);
413 return ((2.*r-q/sqrt(r)-0.5*(r*r-2.*q*sqrt(r)+q*q)*(3.*r*r-6.*r+3.*q*sqrt(r))/denom1)/sqrt(denom1));
414
415 }
416
417static REAL8 XLALHGimri_dFdr(REAL8 E, REAL8 Lz, REAL8 a, REAL8 r) {
418
419 //==================================================
420 // Dimensionless (\partial F)/(\partial r), where F=R/(V_t)^2 and
421 // R and V_t are the Kerr radial and time potentials.
422 // Related to radial coordinate acceleration via
423 // (d^2 r)/(dt^2) = (1/2)dFdr
424 //==================================================
425
426 REAL8 R = pow(E*(a*a+r*r)-a*Lz,2.) - (r*r-2.*r+a*a)*((Lz-a*E)*(Lz-a*E) + r*r);
427 REAL8 Vt = a*(Lz-a*E) + ((r*r+a*a)/(r*r-2.*r+a*a))*(E*(r*r+a*a)-Lz*a);
428 REAL8 dRdr = 4.*E*r*(E*(a*a+r*r)-a*Lz) - 2.*(r-1.)*(pow(Lz-a*E,2.)+r*r) - 2.*r*(r*r-2.*r+a*a);
429 REAL8 dVtdr = 2.*r*(E*(r*r+a*a)-a*Lz)/(r*r-2.*r+a*a) + (r*r+a*a)*(2.*E*r)/(r*r-2.*r+a*a) - (r*r+a*a)*(E*(r*r+a*a)-a*Lz)*(2.*r-2.)/pow(r*r-2.*r+a*a,2.);
430
431 return ( dRdr/pow(Vt,2.) - 2.*R*dVtdr/pow(Vt,3.) );
432 }
433
434static INT4 HGimri_start(REAL8 m, REAL8 M, REAL8 q, REAL8 D, REAL8 Sdotn, REAL8 phi0, REAL8 p0,
435 REAL8Sequence *hplus, REAL8Sequence *hcross, REAL8 dt, UINT4 Npts) {
436
437 //====================================================================================================
438 // Function which evolves an inspiral trajectory and gravitational waveform for a circular intermediate
439 // mass-ratio inspiral, following Huerta & Gair 2011 (see file header above for more info).
440 //
441 // INPUTS:
442 //
443 // m: Compact object mass (in units of seconds)
444 // M: Central BH mass (seconds)
445 // q: Dimensionless reduced spin of the central BH (-1<q<1)
446 // D: Distance to source (in seconds)
447 // Sdotn: Cosine of the inclination angle between the detector line of sight and the orbital
448 // angular momentum vector of the IMRI
449 // phi0: Initial azimuthal angle
450 // p0: Initial orbital radius
451 // hplus: REAL8TimeSequence object to store plus-polarized GW
452 // hcross: REAL8TimeSequence object to store cross-polarized GW
453 // dt: Dimensionless integration time step
454 // Npts: The maximum number of steps to evolve
455 //
456 // NOTE:
457 //
458 // Within this function, masses and source distance are in units of seconds. Central black hole
459 // spin is expressed as the dimensionless reduced spin (-1<q<1). Orbital radius and time are dimensionless,
460 // with dimensions removed by multiplying by the appropriate factor of total mass (m+M).
461 //
462 //====================================================================================================
463
464 //Get intrinsic source parameters. Note: masses passed to HGimri_start() in units of seconds.
465
466 REAL8 mu; //Reduced mass in seconds
467 REAL8 nu; //Reduced mass ratio (mu/(m+M))
468 REAL8 q_abs; //Absolute value of spin q
469 REAL8 final_mass; //Final post-merger mass in units of seconds, from Huerta & Gair (1009.1985) Eq.40
470 REAL8 final_q; //Final post-merger spin, from Huerta & Gair (1009.1985) Eq. 41
471
472 mu = ((m*M)/(m+M));
473 nu = mu/(m+M);
474 q_abs = fabs(q);
475 final_mass = (M+m)*(1. + nu*(sqrt(8./9.)-1.) - 0.498*nu*nu);
476 final_q = q_abs - 0.129*q_abs*q_abs*nu - 0.384*q_abs*nu*nu - 2.686*q_abs*nu + 2.*sqrt(3.)*nu - 3.454*nu*nu + 2.353*nu*nu*nu;
477
478 //Prograde (sign = 1) or retrograde (sign = -1)
479
480 REAL8 sign;
481 if (q < 0.0)
482 sign = -1.0;
483 else
484 sign = 1.0;
485
486 //Conservative corrections as defined in Huerta & Gair (1009.1985v3) Eq.9
487
488 REAL8 d0 = 1./8.;
489 REAL8 d1 = 1975./896.;
490 REAL8 d1_5 = -27.*pi/10.;
491 REAL8 l1_5 = -191./160.;
492 REAL8 d2 = 1152343./451584.;
493
494 //Find ISCO radius and constants
495
496 REAL8 Z1,Z2;
497 REAL8 r_isco; //Radius at ISCO
498 REAL8 E_isco; //Energy at ISCO
499 REAL8 L_isco; //Angular momentum at ISCO
500 REAL8 omega_isco; //Orbital angular frequency at ISCO
501 REAL8 gamma_isco; //Lorentz factor d\tau/dt at ISCO
502
503 Z1 = 1. + pow((1.-q*q),1./3.)*(pow(1.+q, 1./3.)+pow(1.-q, 1./3.));
504 Z2 = sqrt(3.*q*q + Z1*Z1);
505 r_isco = 3. + Z2 - sign*sqrt((3. - Z1)*(3. + Z1 + 2.*Z2));
506 E_isco = (1.+q/pow(r_isco,1.5)-2./r_isco)/sqrt(1.+(2.*q)/pow(r_isco,1.5)-3./r_isco);
507 L_isco = sqrt(r_isco)*(1.+pow(q,2.)/pow(r_isco,2.)-(2.*q)/pow(r_isco,1.5))/sqrt(1.+(2.*q)/pow(r_isco,1.5)-3./r_isco);
508 omega_isco = (1./(pow(r_isco,3/2.)+q))*(1. + nu*(d0+d1/r_isco+(d1_5+q*l1_5)*pow(r_isco,-1.5)+d2*pow(r_isco,-2.0)));
509 gamma_isco = sqrt(1. - 3./r_isco + 2.*q/pow(r_isco,1.5))/(1. + q/pow(r_isco,1.5));
510
511 //Find dimensionless constants governing transition behavior (Ori & Thorne /gr-qc/0003032 Equations 3.9,3.18,3.19)
512
513 REAL8 a_isco;
514 REAL8 b_isco;
515 REAL8 k_isco;
516 REAL8 eps_dot;
517
518 eps_dot = (1.13197 + 0.0713235*q*q + (0.183783*q)/(-2.34484 + 2.15323*q));
519 a_isco = (3./pow(r_isco,6.)) * (r_isco*r_isco + 2.*( q*q*(E_isco*E_isco-1.) - L_isco*L_isco )*r_isco + 10.*pow(L_isco-q*E_isco,2.));
520 b_isco = (2./pow(r_isco,4.)) * ( (L_isco - q*q*E_isco*omega_isco)*r_isco - 3.*(L_isco - q*E_isco)*(1. - q*omega_isco) );
521 k_isco = (32./5.)*pow(omega_isco,7./3.)*eps_dot/gamma_isco;
522
523 //Calculate transition point, plunge point, and light ring radii
524
525 REAL8 r_match; //Radius at which inspiral is matched onto transition (Corresponds to T=-1 in Ori & Thorne Eq. 3.23)
526 REAL8 L_match; //Angular momentum at match
527 REAL8 r_plunge, E_plunge, L_plunge; //Conditions at the start of plunge. Corresponds to T=2 (using Ori & Thorne Eq. 3.25)
528 REAL8 r_lightring; //Radius of Light Ring
529
530 r_match = r_isco + sqrt(1.0*pow(nu*b_isco*k_isco,4./5.)/pow(a_isco,6./5.));
531 L_match = sqrt(r_match)*(1.+pow(q,2.)/pow(r_match,2.)-(2.*q)/pow(r_match,1.5))/sqrt(1.+(2.*q)/pow(r_match,1.5)-3./r_match);
532 r_plunge = r_isco + pow(nu*b_isco*k_isco,2./5.)*pow(a_isco,-3./5.)*(-6./pow(3.412-.75,2.));
533 E_plunge = E_isco - omega_isco*k_isco*pow(a_isco*b_isco*k_isco,-1./5.)*(3.412-0.75)*pow(nu,4./5.);
534 L_plunge = L_isco - k_isco*pow(a_isco*b_isco*k_isco,-1./5.)*(3.412-0.75)*pow(nu,4./5.);
535 r_lightring = 2.*(1.+cos((2./3.)*acos(-q)));
536
537 //Misc variable declaration
538
539 REAL8 dr; //Differential radial step corresponding to the time step dt
540 REAL8 d2phidt2; //Second coordinate time derivatives of phi
541 REAL8 d2rdt2; //Second coordinate time derivatives of radius
542 INT4 i_match=0.; //Index at start of transition.
543 INT4 i_lightring=0; //Index at which the CO reaches the light ring
544 INT4 i_finalplunge = 0; //Used to count the final 10 iterations inside the lightring. Plunge halted when i_finalplunge=10
545
546 REAL8 r; //Instantaneous radial coordinate
547 REAL8 r_old; //Radius at previous time step (used for sanity checks)
548 REAL8 drdt=0.; //Radial velocity
549 REAL8 phi; //Azimuthal position
550 REAL8 dphidt; //Orbital angular velocity
551 REAL8 dphidt_old; //Orbital angular velocity at previous time step (used to numerically compute d2phidt2)
552 REAL8 r_K2,drdt_K2, dphidt_K2,d2rdt2_K2=0.; //For Stage 2 of 4th-order Runge Kutta integration
553 REAL8 r_K3,drdt_K3,dphidt_K3,d2rdt2_K3; //For Stage 3 of RK4 integration
554 REAL8 r_K4,drdt_K4,dphidt_K4,d2rdt2_K4; //For stage 4 of RK4 integration
555
556 //Initialize coordinates
557 r = p0;
558 r_old = p0;
559 phi = phi0;
560 dphidt = (1./(pow(r,3/2.)+q))*(1. + nu*(d0+d1/r+(d1_5+q*l1_5)*pow(r,-1.5)+d2*pow(r,-2.0)));
561 enum stage currentStage = INSPIRAL;
562
563 //SANITY CHECK: Verify that we're starting outside transition. If not, abort.
564 if (r <= r_match) {
565 XLALPrintError("XLAL Error: Beginning inside Transition (Stage 2). Specify larger initial frequency f_min\n");
566 XLAL_ERROR(XLAL_EDOM);
567 }
568
569 //==========================================
570 // Iterate through trajectory
571 // Stage 1: Adiabatic Inspiral
572 // Stage 2: Transition Phase
573 // Stage 3: Plunge
574 //
575 // Evolve trajectory using an 4th order
576 // Runge-Kutta integrator
577 //==========================================
578
579 for (UINT4 i=0; i<Npts; i++) {
580
581 //Direct to the appropriate iterator
582 switch (currentStage) {
583
584 //=============================
585 // Stage 1: Adiabatic Inspiral
586 //=============================
587 case INSPIRAL:
588
589 //Dimensionless radial speed and angular velocity.
590 drdt = XLALHGimri_AngMomFlux(q,r,nu)/XLALHGimri_dLzdr(q,r);
591 dphidt = (1./(pow(r,3/2.)+q))*(1. + nu*(d0+d1/r+(d1_5+q*l1_5)*pow(r,-1.5)+d2*pow(r,-2.0)));
592
593 //RK4 Step 2
594 r_K2 = r + (dt/2.)*drdt;
595 drdt_K2 = XLALHGimri_AngMomFlux(q,r_K2,nu)/XLALHGimri_dLzdr(q,r_K2);
596 dphidt_K2 = (1./(pow(r_K2,3/2.)+q))*(1. + nu*(d0+d1/r_K2+(d1_5+q*l1_5)*pow(r_K2,-1.5)+d2*pow(r_K2,-2.0)));
597
598 //RK4 Step 3
599 r_K3 = r + (dt/2.)*drdt_K2;
600 drdt_K3 = XLALHGimri_AngMomFlux(q,r_K3,nu)/XLALHGimri_dLzdr(q,r_K3);
601 dphidt_K3 = (1./(pow(r_K3,3/2.)+q))*(1. + nu*(d0+d1/r_K3+(d1_5+q*l1_5)*pow(r_K3,-1.5)+d2*pow(r_K3,-2.0)));
602
603 //RK4 Step 4
604 r_K4 = r + (dt)*drdt_K3;
605 drdt_K4 = XLALHGimri_AngMomFlux(q,r_K4,nu)/XLALHGimri_dLzdr(q,r_K4);
606 dphidt_K4 = (1./(pow(r_K4,3/2.)+q))*(1. + nu*(d0+d1/r_K4+(d1_5+q*l1_5)*pow(r_K4,-1.5)+d2*pow(r_K4,-2.0)));
607
608 //Gravitational wave amplitudes, following Huerta & Gair (1009.1985) Eqs.14,15
609 hplus->data[i] = ((4.*mu*pow(dphidt*r,2.))/D) * ((1.+pow(Sdotn,2.))/2.) * cos(2.*phi);
610 hcross->data[i] = ((4.*mu*pow(dphidt*r,2.))/D) * Sdotn * sin(2.*phi);
611
612 //If we've hit the transition point, signal move into the transition regime.
613 dr = (1./6.)*dt*(drdt+2.*drdt_K2+2.*drdt_K3+drdt_K4);
614 if (r+dr < r_match) {
615
616 //Record transition start
617 i_match = i;
618 currentStage = TRANSITION;
619
620 }
621
622 //Update coordinates
623 r_old = r;
624 r += dr;
625 phi += (1./6.)*dt*(dphidt+2.*dphidt_K2+2.*dphidt_K3+dphidt_K4);
626
627 break;
628
629 //=============================
630 // Stage 2: Transition
631 //=============================
632 case TRANSITION:
633
634 dphidt_old = dphidt;
635
636 //Update acceleration using H&G (1009.1985) Eq.29. The gamma^2 factor makes this (to very
637 //good approximation) a coordinate t derivative. Directly compute d^2phi/dt^2.
638 d2rdt2 = gamma_isco*gamma_isco*(
639 - a_isco*(r-r_isco)*(r-r_isco)
640 + b_isco*(L_match-L_isco-nu*k_isco*gamma_isco*(i-i_match)*dt));
641 dphidt = (1./(pow(r,3/2.)+q))*(1. + nu*(d0+d1/r+(d1_5+q*l1_5)*pow(r,-1.5)+d2*pow(r,-2.0)));
642 d2phidt2 = (dphidt-dphidt_old)/(dt);
643
644 //Step 2
645 r_K2 = r+(dt/2.)*drdt;
646 drdt_K2 = drdt+(dt/2.)*d2rdt2;
647 dphidt_K2 = (1./(pow(r_K2,3/2.)+q))*(1. + nu*(d0+d1/r_K2+(d1_5+q*l1_5)*pow(r_K2,-1.5)+d2*pow(r_K2,-2.0)));
648 d2rdt2_K2 = gamma_isco*gamma_isco*(
649 - a_isco*(r_K2-r_isco)*(r_K2-r_isco)
650 + b_isco*(L_match-L_isco-nu*k_isco*gamma_isco*((i-i_match)*dt+(dt/2.))));
651
652 //Step 3
653 r_K3 = r+(dt/2.)*drdt_K2;
654 drdt_K3 = drdt+(dt/2.)*d2rdt2_K2;
655 dphidt_K3 = (1./(pow(r_K3,3/2.)+q))*(1. + nu*(d0+d1/r_K3+(d1_5+q*l1_5)*pow(r_K3,-1.5)+d2*pow(r_K3,-2.0)));
656 d2rdt2_K3 = gamma_isco*gamma_isco*(
657 - a_isco*(r_K3-r_isco)*(r_K3-r_isco)
658 + b_isco*(L_match-L_isco-nu*k_isco*gamma_isco*((i-i_match)*dt+(dt/2.))));
659
660 //Step 4
661 r_K4 = r+(dt)*drdt_K3;
662 drdt_K4 = drdt+(dt/2.)*d2rdt2_K3;
663 dphidt_K4 = (1./(pow(r_K4,3/2.)+q))*(1. + nu*(d0+d1/r_K4+(d1_5+q*l1_5)*pow(r_K4,-1.5)+d2*pow(r_K4,-2.0)));
664 d2rdt2_K4 = gamma_isco*gamma_isco*(
665 - a_isco*(r_K4-r_isco)*(r_K4-r_isco)
666 + b_isco*(L_match-L_isco-nu*k_isco*gamma_isco*((i-i_match)*dt+(dt))));
667
668 //GW amplitudes. Huerta & Gair (1009.1985) Eq. 37 and 38
669 hplus->data[i] = (mu/(2.*D))*(
670 (1. - 2.*(2.*Sdotn*Sdotn-1.)*pow(cos(phi),2.) - 3.*cos(2.*phi))*drdt*drdt
671 + (3. + (2.*Sdotn*Sdotn-1.))*(2.*cos(2.*phi)*dphidt*dphidt + sin(2.*phi)*d2phidt2)*r*r
672 + (4.*(3.+(2.*Sdotn*Sdotn-1.))*sin(2.*phi)*dphidt*drdt
673 + (1.-2.*(2.*Sdotn*Sdotn-1.)*pow(cos(phi),2.)-3.*cos(2.*phi))*d2rdt2)*r);
674
675 hcross->data[i] = (-2.*mu/D)*Sdotn*( sin(2.*(phi))*drdt*drdt
676 + r*r*(-2.*sin(2.*(phi))*dphidt*dphidt + cos(2.*( phi))*d2phidt2)
677 + r*(4.*cos(2.*(phi))*dphidt*drdt + sin(2.*(phi))*d2rdt2));
678
679 //Check if we've reached reached plunge
680 dr = (1./6.)*dt*(drdt+2.*drdt_K2+2.*drdt_K3+drdt_K4);
681 if (r+dr < r_plunge) {
682
683 //Recompute dphidt using geodesic equation to avoid discontinuous d^2phi/dt^2
684 dphidt = (L_plunge*(r-2.)+2.*E_plunge*q)/(E_plunge*(r*r*r+(2.+r)*q*q) - 2.*q*L_plunge);
685
686 //Record plunge start
687 currentStage = PLUNGE;
688
689 }
690
691 //Iterate
692 r_old = r;
693 r += dr;
694 phi += (1./6.)*dt*(dphidt+2.*dphidt_K2+2.*dphidt_K3+dphidt_K4);
695 drdt += (1./6.)*dt*(d2rdt2+2.*d2rdt2_K2+2.*d2rdt2_K3+d2rdt2_K4);
696
697 break;
698
699 //====================================
700 // Stage 3: Plunge
701 //====================================
702
703 //If inside LR, iterate i_finalplunge and move on to PLUNGE case below
704 case FINALPLUNGE:
705
706 //Update counter to track # of iterations inside LR
707 i_finalplunge ++;
708
709 //After 10 iterations inside light ring, stop integration.
710 if (i_finalplunge == 10) {
711
712 i = Npts;
713 break;
714
715 }
716#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
717 __attribute__ ((fallthrough));
718#endif
719
720 //In all cases...
721 case PLUNGE:
722
723 //RK1
724 dphidt_old = dphidt;
725 dphidt = (L_plunge*(r-2.)+2.*E_plunge*q)/(E_plunge*(r*r*r+(2.+r)*q*q) - 2.*q*L_plunge);
726 d2phidt2 = (dphidt-dphidt_old)/(dt);
727 d2rdt2 = (1./2.)*XLALHGimri_dFdr(E_plunge,L_plunge,q,r);
728
729 //RK2
730 r_K2 = r + (dt/2.)*drdt;
731 drdt_K2 = drdt + (dt/2.)*d2rdt2;
732 dphidt_K2 = (L_plunge*(r_K2-2.)+2.*E_plunge*q)/(E_plunge*(r_K2*r_K2*r_K2+(2.+r_K2)*q*q) - 2.*q*L_plunge);
733 d2rdt2 = (1./2.)*XLALHGimri_dFdr(E_plunge,L_plunge,q,r_K2);
734
735 //RK3
736 r_K3 = r + (dt/2.)*drdt_K2;
737 drdt_K3 = drdt + (dt/2.)*d2rdt2_K2;
738 dphidt_K3 = (L_plunge*(r_K3-2.)+2.*E_plunge*q)/(E_plunge*(r_K3*r_K3*r_K3+(2.+r_K3)*q*q) - 2.*q*L_plunge);
739 d2rdt2_K3 = (1./2.)*XLALHGimri_dFdr(E_plunge,L_plunge,q,r_K3);
740
741 //RK4
742 r_K4 = r + (dt)*drdt_K3;
743 drdt_K4 = drdt + (dt)*d2rdt2_K3;
744 dphidt_K4 = (L_plunge*(r_K4-2.)+2.*E_plunge*q)/(E_plunge*(r_K4*r_K4*r_K4+(2.+r_K4)*q*q) - 2.*q*L_plunge);
745 d2rdt2_K4 = (1./2.)*XLALHGimri_dFdr(E_plunge,L_plunge,q,r_K4);
746
747 //GW amplitudes
748 hplus->data[i] = (mu/(2.*D))*((1. - 2.*(2.*Sdotn*Sdotn-1.)*pow(cos(phi),2.) - 3.*cos(2.*phi))*drdt*drdt
749 + (3. + (2.*Sdotn*Sdotn-1.))*(2.*cos(2.*phi)*dphidt*dphidt + sin(2.*phi)*d2phidt2)*r*r
750 + (4.*(3.+(2.*Sdotn*Sdotn-1.))*sin(2.*phi)*dphidt*drdt
751 + (1.-2.*(2.*Sdotn*Sdotn-1.)*pow(cos(phi),2.)-3.*cos(2.*phi))*d2rdt2)*r);
752
753 hcross->data[i] = (-2.*mu/D)*Sdotn*( sin(2.*(phi))*drdt*drdt
754 + r*r*(-2.*sin(2.*(phi))*dphidt*dphidt + cos(2.*( phi))*d2phidt2)
755 + r*(4.*cos(2.*(phi))*dphidt*drdt + sin(2.*(phi))*d2rdt2));
756
757 //Update coordinates
758 r_old = r;
759 r += (1./6.)*dt*(drdt+2.*drdt_K2+2.*drdt_K3+drdt_K4);
760 phi += (1./6.)*dt*(dphidt+2.*dphidt_K2+2.*dphidt_K3+dphidt_K4);
761 drdt += (1./6.)*dt*(d2rdt2+2.*d2rdt2_K2+2.*d2rdt2_K3+d2rdt2_K4);
762
763 //check to see if we've passed light ring
764 if ((r < r_lightring) && (currentStage == PLUNGE)) {
765
766 //If so, initiate end of plunge
767 currentStage = FINALPLUNGE;
768 i_lightring = i;
769
770 }
771
772 break;
773
774 }
775
776 //===============
777 // SANITY CHECKS
778 //===============
779
780 //If r is negative or NaN, raise error.
781 if (r != r || r < 0.0) {
782 XLALPrintError("XLAL Error: Radius is negative or NaN. ");
783 if (currentStage == INSPIRAL) XLALPrintError("Halting INSPIRAL\n.");
784 if (currentStage == TRANSITION) XLALPrintError("Halting TRANSITION.\n");
785 if (currentStage == PLUNGE) XLALPrintError("Halting PLUNGE.\n");
786 if (currentStage == FINALPLUNGE) {
787 if (i_finalplunge <= 10)
788 XLALPrintError("Failed to complete 10 steps inside light ring. Use smaller time step dt. ");
789 XLALPrintError("Halting FINALPLUNGE.\n");
790 }
791 XLAL_ERROR(XLAL_EFAILED);
792 }
793
794 //If r has increased since the last time step, raise error.
795 if (r > r_old) {
796 XLALPrintError("XLAL Error: Increasing radius. ");
797 if (currentStage == INSPIRAL) XLALPrintError("Halting INSPIRAL\n.");
798 if (currentStage == TRANSITION) XLALPrintError("Halting TRANSITION.\n");
799 if (currentStage == PLUNGE) XLALPrintError("Halting PLUNGE.\n");
800 if (currentStage == FINALPLUNGE) {
801 if (i_finalplunge <= 10)
802 XLALPrintError("Failed to complete 10 steps inside light ring. Use smaller time step dt. ");
803 XLALPrintError("Halting FINALPLUNGE.\n");
804 }
805 XLAL_ERROR(XLAL_EFAILED);
806 }
807
808 //If we've reached the max array length and failed to plunge, either something has gone wrong numerically,
809 //or Npts is too small. In the latter case, choose smaller time steps or increase Npts in XLALHGimri_generator() below
810 if (i == Npts && currentStage != FINALPLUNGE) {
811 XLALPrintError("XLAL Error: Reached Npts before finishing plunge. Binary either failed to plunge, or Npts is too small.\n");
812 XLAL_ERROR(XLAL_EFAILED);
813 }
814
815 }
816
817
818 //=====================================================================
819 // Get ringdown frequencies by interpolating data presented in Table 2
820 // of Berti et al 2006 (gr-qc/0512160). See Huerta & Gair for details.
821 //=====================================================================
822
823 INT4 N = 11;
824
825 //Final q values
826 REAL8 fq[11] = {0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.98};
827
828 //Real and imaginary dimensionless frequencies for mode (l=2,m=2,n=0). From Berti et. al.
829 REAL8 w0_arr[11] = {0.3737, 0.3870, 0.4021, 0.4195, 0.4398, 0.4641, 0.4940, 0.5326, 0.5860, 0.6716, 0.8254};
830 REAL8 wi0_arr[11] = {0.0890, 0.0887, 0.0883, 0.0877, 0.0869, 0.0856, 0.0838, 0.0808, 0.0756, 0.0649, 0.0386};
831
832 //Real and imaginary dimensionless frequencies for mode (l=2,m=2,n=1). From Berti et. al.
833 REAL8 w1_arr[11] = {0.3467, 0.3619, 0.3790, 0.3984, 0.4208, 0.4474, 0.4798, 0.5212, 0.5779, 0.6677, 0.8249};
834 REAL8 wi1_arr[11] = {0.2739, 0.2725, 0.2705, 0.2680, 0.2647, 0.2602, 0.2538, 0.2442, 0.2281, 0.1953, 0.1159};
835
836 //Real and imaginary dimensionless frequencies for mode (l=2,m=2,n=2). From Berti et. al.
837 REAL8 w2_arr[11] = {0.3011, 0.3192, 0.3393, 0.3619, 0.3878, 0.4179, 0.4542, 0.4999, 0.5622, 0.6598, 0.8238};
838 REAL8 wi2_arr[11] = {0.4783, 0.4735, 0.4679, 0.4613, 0.4533, 0.4433, 0.4303, 0.4123, 0.3839, 0.3275, 0.1933};
839
840 //Set up interpolator
841 const gsl_interp_type *bn = gsl_interp_cspline;
842 gsl_interp_accel *acc = gsl_interp_accel_alloc();
843
844 //Define a bunch of interpolation splines
845 gsl_spline *spline0 = gsl_spline_alloc(bn,N);
846 gsl_spline *spline1 = gsl_spline_alloc(bn,N);
847 gsl_spline *spline2 = gsl_spline_alloc(bn,N);
848 gsl_spline *spline3 = gsl_spline_alloc(bn,N);
849 gsl_spline *spline4 = gsl_spline_alloc(bn,N);
850 gsl_spline *spline5 = gsl_spline_alloc(bn,N);
851
852 //Initialize splines
853 gsl_spline_init(spline0, fq, w0_arr, N);
854 gsl_spline_init(spline1, fq, wi0_arr, N);
855 gsl_spline_init(spline2, fq, w1_arr, N);
856 gsl_spline_init(spline3, fq, wi1_arr, N);
857 gsl_spline_init(spline4, fq, w2_arr, N);
858 gsl_spline_init(spline5, fq, wi2_arr, N);
859
860 //Get the real and imaginary frequencies for our final BH
861 REAL8 w0, wi0;
862 REAL8 w1, wi1;
863 REAL8 w2, wi2;
864
865 w0 = gsl_spline_eval(spline0, final_q, acc);
866 wi0 = gsl_spline_eval(spline1, final_q, acc);
867 w1 = gsl_spline_eval(spline2, final_q, acc);
868 wi1 = gsl_spline_eval(spline3, final_q, acc);
869 w2 = gsl_spline_eval(spline4, final_q, acc);
870 wi2 = gsl_spline_eval(spline5, final_q, acc);
871
872 //Free memory
873 gsl_spline_free(spline0);
874 gsl_spline_free(spline1);
875 gsl_spline_free(spline2);
876 gsl_spline_free(spline3);
877 gsl_spline_free(spline4);
878 gsl_spline_free(spline5);
879
880
881 //==================================================================
882 // Build interpolation functions for h_plus and h_cross at 20 time
883 // steps centered at r_lightring.
884 //==================================================================
885
886 //Define arrays to hold the final 20 values of hplus, hcross, and their derivatives during plunge.
887 //Note, derivatives are technically dh/di, not dh/dt. Get the latter using dh/dt = (dh/di)/dt
888 REAL8 i_interp[20];
889 REAL8 hp_interp[20];
890 REAL8 hx_interp[20];
891
892 //Unknown constant amplitudes in ringdown equation.
893 REAL8 a0n, a0p;
894 REAL8 a1n, a1p;
895 REAL8 a2n, a2p;
896 REAL8 a0c, a0cp;
897 REAL8 a1c, a1cp;
898 REAL8 a2c, a2cp;
899
900 //Read out the final twenty values of hplus[] and hcross[].
901 //The lightring bisects i_lightring (i_interp[10]) and i_lightring-1 (i_interp[9])
902 for (INT4 i=0; i<20; i++) {
903 i_interp[i] = i+(i_lightring-10);
904 hp_interp[i] = hplus->data[i+(i_lightring-10)]*D;
905 hx_interp[i] = hcross->data[i+(i_lightring-10)]*D;
906 }
907
908 //Make some splines!
909 gsl_spline *spline_hp = gsl_spline_alloc(gsl_interp_cspline, 20); //To interpolate across hplus
910 gsl_spline *spline_hx = gsl_spline_alloc(gsl_interp_cspline, 20); //To interpolate across hcross
911 gsl_spline_init(spline_hp, i_interp, hp_interp, 20);
912 gsl_spline_init(spline_hx, i_interp, hx_interp, 20);
913
914 //Define hp, hx, and their first, second, and third derivatives with respect to i at three steps across the light radius.
915 REAL8 hp_1, dhpdi_1;
916 REAL8 hx_1, dhxdi_1;
917 REAL8 hp_2, dhpdi_2;
918 REAL8 hx_2, dhxdi_2;
919 REAL8 hp_3, dhpdi_3;
920 REAL8 hx_3, dhxdi_3;
921
922 //Evaluate first point at i_interp[9] = i_lightring-1
923 hp_1 = gsl_spline_eval(spline_hp, i_interp[9], acc);
924 dhpdi_1 = gsl_spline_eval_deriv(spline_hp, i_interp[9], acc);
925 hx_1 = gsl_spline_eval(spline_hx, i_interp[9], acc);
926 dhxdi_1 = gsl_spline_eval_deriv(spline_hx, i_interp[9], acc);
927
928 //Evaluate second point at i_interp[10] = i_lightring
929 hp_2 = gsl_spline_eval(spline_hp, i_interp[10], acc);
930 dhpdi_2 = gsl_spline_eval_deriv(spline_hp, i_interp[10], acc);
931 hx_2 = gsl_spline_eval(spline_hx, i_interp[10], acc);
932 dhxdi_2 = gsl_spline_eval_deriv(spline_hx, i_interp[10], acc);
933
934 //Evaluate third point at i_interp[11] = i_lightring+1
935 hp_3 = gsl_spline_eval(spline_hp, i_interp[11], acc);
936 dhpdi_3 = gsl_spline_eval_deriv(spline_hp, i_interp[11], acc);
937 hx_3 = gsl_spline_eval(spline_hx, i_interp[11], acc);
938 dhxdi_3 = gsl_spline_eval_deriv(spline_hx, i_interp[11], acc);
939
940 //Free memory
941 gsl_spline_free(spline_hp);
942 gsl_spline_free(spline_hx);
943
944 //==============================================================
945 // Match plunging waveform onto QNM ringdown waveform by
946 // numerically finding the QNM amplitudes which match the
947 // waveform amplitudes and derivatives at the three points above.
948 //
949 // See comments in declarations of XLALHGimri_Plus/CrossEquations()
950 // for more detail.
951 // =============================================================
952
953 //Create an instance of rparams, populate
954 struct rparams p;
955 p.dt = dt;
956 p.M = (M+m);
957 p.Sdotn = Sdotn;
958 p.final_mass = final_mass;
959 p.final_q = final_q;
960 p.w0 = w0;
961 p.w1 = w1;
962 p.w2 = w2;
963 p.wi0 = wi0;
964 p.wi1 = wi1;
965 p.wi2 = wi2;
966 p.ab_factor = (1./8.)*sqrt(5./pi)*final_mass;
967 p.hp_1 = hp_1;
968 p.hp_2 = hp_2;
969 p.hp_3 = hp_3;
970 p.dhpdi_1 = dhpdi_1;
971 p.dhpdi_2 = dhpdi_2;
972 p.dhpdi_3 = dhpdi_3;
973 p.hx_1 = hx_1;
974 p.hx_2 = hx_2;
975 p.hx_3 = hx_3;
976 p.dhxdi_1 = dhxdi_1;
977 p.dhxdi_2 = dhxdi_2;
978 p.dhxdi_3 = dhxdi_3;
979
980 //Initialize root solvers
981 const gsl_multiroot_fsolver_type *T;
982 const gsl_multiroot_fsolver_type *Q;
983 gsl_multiroot_fsolver *plus_solver;
984 gsl_multiroot_fsolver *cross_solver;
985 INT4 statusp, statusc;
986 size_t iter = 0;
987 const size_t nmatch = 6;
988
989 //Define target functions to be XLALHGimri_Plus/CrossEquations
990 gsl_multiroot_function f_plus = {&XLALHGimri_PlusEquations, nmatch, &p};
991 gsl_multiroot_function f_cross = {&XLALHGimri_CrossEquations, nmatch, &p};
992
993 //Define two 6D vectors, and feed initial guesses
994 REAL8 x_init[6] = {-0.392254, 4.575194, -0.870431, 5.673678, 0.979356, -3.637467};
995 REAL8 z_init[6] = {-0.392254, 4.575194, -0.870431, 5.673678, 0.979356, -3.637467};
996 gsl_vector *xmr = gsl_vector_alloc(nmatch);
997 gsl_vector *z = gsl_vector_alloc(nmatch);
998
999 for (size_t i=0; i<nmatch; i++) {
1000 gsl_vector_set(xmr, i, x_init[i]);
1001 gsl_vector_set(z, i, z_init[i]);
1002 }
1003
1004 //Set up solvers
1005 T = gsl_multiroot_fsolver_hybrids;
1006 Q = gsl_multiroot_fsolver_hybrids;
1007 plus_solver = gsl_multiroot_fsolver_alloc(T, 6);
1008 cross_solver = gsl_multiroot_fsolver_alloc(Q, 6);
1009 gsl_multiroot_fsolver_set(plus_solver, &f_plus, xmr);
1010 gsl_multiroot_fsolver_set(cross_solver, &f_cross, z);
1011
1012 REAL8 fminval = 1.e-5;
1013 size_t MaxITS = 200000;
1014
1015 //Find root of XLALHGimri_PlusEquations
1016 do {
1017 iter++;
1018 statusp = gsl_multiroot_fsolver_iterate(plus_solver);
1019 if (statusp) break;
1020 statusp = gsl_multiroot_test_residual(plus_solver->f, fminval);
1021 }
1022 while (statusp == GSL_CONTINUE && iter < MaxITS);
1023
1024 //Find root of XLALHGimri_CrossEquations
1025 iter = 0;
1026 do {
1027 iter++;
1028 statusc = gsl_multiroot_fsolver_iterate(cross_solver);
1029 if (statusc) break;
1030 statusc = gsl_multiroot_test_residual(cross_solver->f, fminval);
1031 }
1032 while (statusc == GSL_CONTINUE && iter < MaxITS);
1033
1034 if (statusc != GSL_SUCCESS) {
1035 XLALPrintError("XLAL Error: Ringdown rootfinding failed.\n");
1036 XLAL_ERROR(XLAL_EFAILED);
1037 }
1038
1039 //Read out fit parameters
1040 a0n = gsl_vector_get(plus_solver->x,0);
1041 a0p = gsl_vector_get(plus_solver->x,1);
1042 a1n = gsl_vector_get(plus_solver->x,2);
1043 a1p = gsl_vector_get(plus_solver->x,3);
1044 a2n = gsl_vector_get(plus_solver->x,4);
1045 a2p = gsl_vector_get(plus_solver->x,5);
1046 a0c = gsl_vector_get(cross_solver->x,0);
1047 a0cp = gsl_vector_get(cross_solver->x,1);
1048 a1c = gsl_vector_get(cross_solver->x,2);
1049 a1cp = gsl_vector_get(cross_solver->x,3);
1050 a2c = gsl_vector_get(cross_solver->x,4);
1051 a2cp = gsl_vector_get(cross_solver->x,5);
1052
1053 //Free memory
1054 gsl_multiroot_fsolver_free(plus_solver);
1055 gsl_multiroot_fsolver_free(cross_solver);
1056 gsl_vector_free(xmr);
1057 gsl_vector_free(z);
1058
1059 //=====================================================================================
1060 // Recompute wave amplitudes from (i_lightring-1), the first matching point, onwards
1061 //=====================================================================================
1062
1063 REAL8 aone, atwo;
1064 REAL8 hp0, hp1, hp2;
1065 REAL8 hc0, hc1, hc2;
1066 REAL8 ringdown_amp;
1067 REAL8 hp_ringdown, hc_ringdown;
1068
1069 aone = 1.+ 2.*Sdotn + Sdotn*Sdotn;
1070 atwo = 2.+ Sdotn - 4.*Sdotn*Sdotn - 3.*Sdotn*Sdotn*Sdotn;
1071 ringdown_amp = (1./8.)*sqrt(5./pi)*(final_mass/D);
1072
1073 for (UINT4 i=i_lightring; i<Npts; i++) {
1074
1075 //NOTE: final_q and the ringdown frequencies wi are nondimensionalized by final_mass. Their direct product is therefore already
1076 //dimensionless. When multiplying the wi's by dt, however, factors of (M/final_mass), however, are needed to properly rescale the wi's.
1077 hp0 = exp(-(i-(i_lightring-1.))*dt*(wi0*(m+M)/final_mass))*( aone-(2./9.)*final_q*w0*atwo )*( a0n*cos(w0*(i-(i_lightring-1.))*(dt*(m+M)/final_mass)) - a0p*sin(w0*(i-(i_lightring-1.))*(dt*(m+M)/final_mass)));
1078 hp1 = exp(-(i-(i_lightring-1.))*dt*(wi1*(m+M)/final_mass))*( aone-(2./9.)*final_q*w1*atwo )*( a1n*cos(w1*(i-(i_lightring-1.))*(dt*(m+M)/final_mass)) - a1p*sin(w1*(i-(i_lightring-1.))*(dt*(m+M)/final_mass)));
1079 hp2 = exp(-(i-(i_lightring-1.))*dt*(wi2*(m+M)/final_mass))*( aone-(2./9.)*final_q*w2*atwo )*( a2n*cos(w2*(i-(i_lightring-1.))*(dt*(m+M)/final_mass)) - a2p*sin(w2*(i-(i_lightring-1.))*(dt*(m+M)/final_mass)));
1080 hc0 = exp(-(i-(i_lightring-1.))*dt*(wi0*(m+M)/final_mass))*( aone-(2./9.)*final_q*w0*atwo )*( a0c*cos(w0*(i-(i_lightring-1.))*(dt*(m+M)/final_mass)) + a0cp*sin(w0*(i-(i_lightring-1.))*(dt*(m+M)/final_mass)));
1081 hc1 = exp(-(i-(i_lightring-1.))*dt*(wi1*(m+M)/final_mass))*( aone-(2./9.)*final_q*w1*atwo )*( a1c*cos(w1*(i-(i_lightring-1.))*(dt*(m+M)/final_mass)) + a1cp*sin(w1*(i-(i_lightring-1.))*(dt*(m+M)/final_mass)));
1082 hc2 = exp(-(i-(i_lightring-1.))*dt*(wi2*(m+M)/final_mass))*( aone-(2./9.)*final_q*w2*atwo )*( a2c*cos(w2*(i-(i_lightring-1.))*(dt*(m+M)/final_mass)) + a2cp*sin(w2*(i-(i_lightring-1.))*(dt*(m+M)/final_mass)));
1083
1084 //Get hplus and hcross
1085 hp_ringdown = (hp0+hp1+hp2);
1086 hc_ringdown = -(hc0+hc1+hc2);
1087
1088 //Record wave amplitudes
1089 hplus->data[i] = ringdown_amp*hp_ringdown;
1090 hcross->data[i] = ringdown_amp*hc_ringdown;
1091
1092 //Halt when wave amplitudes become sufficiently low
1093 if (fabs(hplus->data[i]) < 1.e-30) {
1094
1095 hplus->length=i-1; //Redefine REAL8TimeSeries length to be the current index.
1096 hcross->length=i-1;
1097 i = Npts; //Halt.
1098
1099 }
1100
1101 }
1102
1103 return i_lightring;
1104
1105 }
1106
1107
1108/**
1109 * @addtogroup LALSimInspiralHGimri_c
1110 * @brief Routines for generating the Huerta-Gair Intermediate-Mass-Ratio Inspiral model.
1111 * @{
1112 */
1113
1114//===================================================
1115// Generator function for the Huerta-Gair IMRI model.
1116//===================================================
1117
1118INT4 XLALHGimriGenerator(
1119 REAL8TimeSeries **hplus,
1120 REAL8TimeSeries **hcross,
1121 REAL8 phi0, //Initial phi
1122 REAL8 dt, //Integration time step (in seconds)
1123 REAL8 m1, //BH mass (passed in kg)
1124 REAL8 m2, //CO mass (passed in kg)
1125 REAL8 f_min, //Initial frequency (passed in Hz)
1126 REAL8 r, //Distance to system (passed in meters)
1127 REAL8 inc, //Inclination angle between detector line of sight n and central BH spin axis
1128 REAL8 s1z //Central BH reduced spin
1129 ) {
1130
1131 //Parse input, adjusting dimensiones
1132 REAL8 q = s1z; //Rename spin
1133 REAL8 Sdotn = cos(inc); //Cosine(inclination)
1134 REAL8 M = m1/LAL_MSUN_SI; //Units of solar masses
1135 REAL8 m = m2/LAL_MSUN_SI; //Units of solar masses
1136 REAL8 dist = r/(LAL_PC_SI*1.e9); //Units of Gpc
1137
1138 //SANITY CHECKS
1139 if (f_min < 0.) {
1140 XLALPrintError("XLAL Error: Minimum frequency is %f. f_min must be greater than 0\n",f_min);
1141 XLAL_ERROR(XLAL_EDOM);
1142 }
1143
1144 if (q >= 1 || q<=-1) {
1145 XLALPrintError("XLAL Error: Magnitude of BH spin %f must be less than 1.\n",q);
1146 XLAL_ERROR(XLAL_EDOM);
1147 }
1148
1149 if (m/M > 1./10.) {
1150 XLALPrintWarning("XLAL Warning: Mass-ratio is %f. Code likely inaccurate for mass ratios m2/m1 > 1/10\n",m/M);
1151 }
1152
1153 if (dist <= 0) {
1154 XLALPrintError("XLAL Error: Distance %f Mpc to source must be greater than 0 Mpc.\n",dist);
1155 XLAL_ERROR(XLAL_EDOM);
1156 }
1157
1158 //==========================================================================
1159 //Numerically solve for initial radius corresponding to initial gravitational
1160 //wave frequency of f_min (or an initial orbital frequency of twice f_min).
1161 //==========================================================================
1162
1163 //Construct a gsl rootfinder.
1164 INT4 status;
1165 size_t iter=0, max_iter = 100;
1166 const gsl_root_fsolver_type *T;
1167 gsl_root_fsolver *s;
1168
1169 //Populate pParams object
1170 struct pParams p0Params;
1171 p0Params.m = m;
1172 p0Params.M = M;
1173 p0Params.q = q;
1174 p0Params.f_min = f_min;
1175
1176 //Define target function
1177 gsl_function F;
1178 F.function = &XLALHGimri_initialP;
1179 F.params = &p0Params;
1180
1181 //Initialize p0 and the search bounds.
1182 REAL8 p0 = 0.;
1183 REAL8 p_high = 500.;
1184 REAL8 p_low = 0.;
1185
1186 //Set a reasonable lower limit by equating p_low to the ISCO radius.
1187 if (q==0.) p_low = 6.;
1188 else {
1189 REAL8 Z1 = 1. + pow((1.-q*q),1./3.)*(pow(1.+q, 1./3.)+pow(1.-q, 1./3.));
1190 REAL8 Z2 = sqrt(3.*q*q + Z1*Z1);
1191 p_low = 3. + Z2 - (q/fabs(q))*sqrt((3. - Z1)*(3. + Z1 + 2.*Z2));
1192 }
1193
1194 //Find p0
1195 T = gsl_root_fsolver_brent;
1196 s = gsl_root_fsolver_alloc(T);
1197 gsl_root_fsolver_set(s, &F, p_low, p_high);
1198
1199 do {
1200 iter++;
1201 status = gsl_root_fsolver_iterate(s);
1202 p0 = gsl_root_fsolver_root(s);
1203 p_low = gsl_root_fsolver_x_lower(s);
1204 p_high = gsl_root_fsolver_x_upper(s);
1205 status = gsl_root_test_interval (p_low, p_high, 0, 0.001);
1206 }
1207 while (status == GSL_CONTINUE && iter < max_iter);
1208 gsl_root_fsolver_free(s);
1209
1210 if (status != GSL_SUCCESS) {
1211 XLALPrintError("XLAL Error: Rootfinding failed. Could not find inital radius ISCO < r < 500M satisfying specified initial frequency.\n");
1212 XLAL_ERROR(XLAL_EFAILED); //"generic failure"
1213 }
1214
1215 if (p0 >= 20.) {
1216 XLALPrintWarning("XLAL Warning: Beginning at large initial radius p0 = %fM. Code will require long runtime to compute waveform.\n",p0);
1217 }
1218
1219 //==============================================================
1220 //Create REAL8TimeSeries to hold waveform data, start simulation
1221 //==============================================================
1222
1223 REAL8 fDyne = 0.0;
1224 size_t length = 1000/dt;
1225 LIGOTimeGPS tc = LIGOTIMEGPSZERO;
1226 *hplus = XLALCreateREAL8TimeSeries("h_plus",&tc,fDyne,dt,&lalStrainUnit,length);
1227 *hcross = XLALCreateREAL8TimeSeries("h_cross",&tc,fDyne,dt,&lalStrainUnit,length);
1228 if (*hplus == NULL || *hcross == NULL)
1229 XLAL_ERROR(XLAL_EFUNC);
1230
1231 INT4 i_ref = 0;
1232 i_ref = HGimri_start(m*Msun_sec,M*Msun_sec,q,dist*GPC_sec,Sdotn,phi0,p0,(*hplus)->data,(*hcross)->data,dt/((m+M)*Msun_sec),length);
1233
1234 //Redefine reference epoch to beginning of ringdown
1235 XLALGPSAdd(&(*hplus)->epoch,-1.*i_ref*dt);
1236 XLALGPSAdd(&(*hcross)->epoch,-1.*i_ref*dt);
1237
1238 return 0;
1239
1240 }
1241
1242/** @} */
c1
const double c1
Definition: LALSimIMRLackeyTidal2013.c:51
c2
const double c2
Definition: LALSimIMRLackeyTidal2013.c:52
HGimri_start
static INT4 HGimri_start(REAL8 m, REAL8 M, REAL8 q, REAL8 D, REAL8 Sdotn, REAL8 phi0, REAL8 p0, REAL8Sequence *hplus, REAL8Sequence *hcross, REAL8 dt, UINT4 Npts)
Definition: LALSimInspiralHGimri.c:434
XLALHGimri_AngMomFlux
static REAL8 XLALHGimri_AngMomFlux(REAL8 q, REAL8 r, REAL8 nu)
Definition: LALSimInspiralHGimri.c:111
pi
#define pi
Definition: LALSimInspiralHGimri.c:37
XLALHGimri_CrossEquations
static INT4 XLALHGimri_CrossEquations(const gsl_vector *z, void *params, gsl_vector *g)
Definition: LALSimInspiralHGimri.c:314
XLALHGimri_dFdr
static REAL8 XLALHGimri_dFdr(REAL8 E, REAL8 Lz, REAL8 a, REAL8 r)
Definition: LALSimInspiralHGimri.c:417
XLALHGimri_dLzdr
static REAL8 XLALHGimri_dLzdr(REAL8 q, REAL8 r)
Definition: LALSimInspiralHGimri.c:401
XLALHGimri_PlusEquations
static INT4 XLALHGimri_PlusEquations(const gsl_vector *x, void *params, gsl_vector *f)
Definition: LALSimInspiralHGimri.c:225
GPC_sec
#define GPC_sec
Definition: LALSimInspiralHGimri.c:41
Msun_sec
#define Msun_sec
Definition: LALSimInspiralHGimri.c:40
INSPIRAL
@ INSPIRAL
Definition: LALSimInspiralHGimri.c:55
TRANSITION
@ TRANSITION
Definition: LALSimInspiralHGimri.c:55
FINALPLUNGE
@ FINALPLUNGE
Definition: LALSimInspiralHGimri.c:55
XLALHGimri_initialP
static REAL8 XLALHGimri_initialP(REAL8 p0, void *params)
Definition: LALSimInspiralHGimri.c:83
s
int s
Definition: bh_qnmode.c:137
M
REAL8 M
Definition: bh_qnmode.c:133
dt
double dt
Definition: bh_ringdown.c:113
i
double i
Definition: bh_ringdown.c:118
Q
const double Q
Definition: exact_derivatives-Hreal.c:33
w2
const double w2
Definition: exact_derivatives-HrealHeader.c:55
__attribute__
#define __attribute__(x)
LAL_MSUN_SI
#define LAL_MSUN_SI
LAL_PC_SI
#define LAL_PC_SI
LIGOTIMEGPSZERO
#define LIGOTIMEGPSZERO
REAL8
double REAL8
UINT4
uint32_t UINT4
INT4
int32_t INT4
XLALHGimriGenerator
INT4 XLALHGimriGenerator(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, REAL8 phi0, REAL8 dt, REAL8 m1, REAL8 m2, REAL8 f_min, REAL8 r, REAL8 inc, REAL8 s1z)
Definition: LALSimInspiralHGimri.c:1118
r
static const INT4 r
m
static const INT4 m
q
static const INT4 q
a
static const INT4 a
XLALCreateREAL8TimeSeries
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
lalStrainUnit
const LALUnit lalStrainUnit
XLAL_ERROR
#define XLAL_ERROR(...)
XLALPrintError
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLALPrintWarning
int int XLALPrintWarning(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_EFUNC
XLAL_EFUNC
XLAL_EDOM
XLAL_EDOM
XLAL_EFAILED
XLAL_EFAILED
XLALGPSAdd
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
mu
list mu
status
string status
N
int N
p
p
x
x
REAL8Vector::data
REAL8 * data
params
Definition: burst.c:245
f_min
double f_min
Definition: unicorn.c:22