LALSimulation  5.4.0.1-fe68b98
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
44 static REAL8 XLALHGimri_AngMomFlux(REAL8 q, REAL8 r, REAL8 nu);
45 static INT4 XLALHGimri_PlusEquations(const gsl_vector * x, void *params, gsl_vector * f);
46 static INT4 XLALHGimri_CrossEquations(const gsl_vector * z, void *params, gsl_vector * g);
47 static REAL8 XLALHGimri_dFdr(REAL8 E, REAL8 Lz, REAL8 a, REAL8 r);
48 static REAL8 XLALHGimri_initialP(REAL8 p0, void *params);
49 static REAL8 XLALHGimri_dLzdr(REAL8 q, REAL8 r);
50 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);
51 
52 INT4 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
55 enum 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)
58 struct 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
76 struct 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 
83 static 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 
111 static 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 
225 static 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 
314 static 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 
401 static 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 
417 static 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 
434 static 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 
1118 INT4 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)
x
list x
p
list p
mu
list mu
status
string status
N
int N
REAL8Vector::data
REAL8 * data
params
Definition: burst.c:245
f_min
double f_min
Definition: unicorn.c:22