Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInspiral 5.0.3.1-6c6b863
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALInspiralBCVBank.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2007 Thomas Cokelaer
3*
4* This program is free software; you can redistribute it and/or modify
5* it under the terms of the GNU General Public License as published by
6* the Free Software Foundation; either version 2 of the License, or
7* (at your option) any later version.
8*
9* This program is distributed in the hope that it will be useful,
10* but WITHOUT ANY WARRANTY; without even the implied warranty of
11* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12* GNU General Public License for more details.
13*
14* You should have received a copy of the GNU General Public License
15* along with with program; see the file COPYING. If not, write to the
16* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17* MA 02110-1301 USA
18*/
19
20#include <stdio.h>
21#include <lal/LALInspiralBank.h>
22#include <lal/AVFactories.h>
23#include <lal/SeqFactories.h>
24#include <lal/LALStdio.h>
25#include <lal/FindRoot.h>
26
27#ifdef __GNUC__
28#define UNUSED __attribute__ ((unused))
29#else
30#define UNUSED
31#endif
32
33/**
34 * \defgroup LALInspiralBCVBank_c Module LALInspiralBCVBank.c
35 * \ingroup LALInspiralBank_h
36 * \brief BCV template bank
37 */
38/** @{ */
39
40/**
41 * \brief Lay a flat grid of BCV templates in the user specified range
42 * of the parameters \f$(\psi_0, \psi_3)\f$ in \c coarseIn structure
43 * \author Cokelaer, T
44 *
45 * ### Description ###
46 *
47 * Given the range of the parameters \f$(\psi_0, \psi_3),\f$
48 * the number of templates in the \c fCut direction,
49 * \e minimalMatch, noise spectral density, upper and
50 * lower frequency cutoffs (all in the input structure \c coarseIn)
51 * this routine outputs the list of templates in the BCV bank
52 * for the parameters \f$(\psi_0, \psi_3, f_{\mathrm{cut}}).\f$
53 *
54 * ### Algorithm ###
55 *
56 * A flat signal manifold is assumed and templates are laid
57 * uniform in the three dimensions. See below for an explanation
58 * of how templates are chosen in the \c fcut direction.
59 */
60void
62 LALStatus *status, /**< LAL status pointer */
63 InspiralTemplateList **list, /**< [out] an array containing the template bank parameters. */
64 INT4 *nlist, /**< [out] the number of templates in bank */
65 InspiralCoarseBankIn coarseIn /**< UNDOCUMENTED */
66 )
67
68{
69 INT4 j = 0;
70 INT4 nlistOld = 0;
71 /*do we really need static declaration here ? */
72 static InspiralBankParams bankParams;
73 static InspiralMetric metric;
75 static CreateVectorSequenceIn in;
76 static REAL4VectorSequence *tempList = NULL;
77
80
81 ASSERT( coarseIn.psi0Min > 0., status,
82 LALINSPIRALBANKH_ESIZE, LALINSPIRALBANKH_MSGESIZE );
83 ASSERT( coarseIn.psi0Max > coarseIn.psi0Min, status,
84 LALINSPIRALBANKH_ESIZE, LALINSPIRALBANKH_MSGESIZE );
85 ASSERT( coarseIn.psi3Min < 0., status,
86 LALINSPIRALBANKH_ESIZE, LALINSPIRALBANKH_MSGESIZE );
87 ASSERT( coarseIn.psi3Max > coarseIn.psi3Min, status,
88 LALINSPIRALBANKH_ESIZE, LALINSPIRALBANKH_MSGESIZE );
89 ASSERT( coarseIn.LowGM < coarseIn.HighGM, status,
90 LALINSPIRALBANKH_ESIZE, LALINSPIRALBANKH_MSGESIZE );
91
92 /* populate the param structure so as to call
93 * ComputeMetricBCV */
94 params.fLower = coarseIn.fLower;
95 params.fCutoff = coarseIn.fUpper;
96 params.alpha = coarseIn.alpha;
97 /* Get the BCV metric in psi0psi3 plane. */
99 &metric, &coarseIn.shf, &params );
101
102 /* print the metric if lalinfo is on*/
103 if ( lalDebugLevel & LALINFO )
104 {
105 REAL8 dx0 = sqrt( 2.L * (1.L-coarseIn.mmCoarse)/metric.g00 );
106 REAL8 dx1 = sqrt( 2.L * (1.L-coarseIn.mmCoarse)/metric.g11 );
107 LALPrintError( "G00=%e G01=%e G11=%e\n",
108 metric.G00, metric.G01, metric.G11 );
109 LALPrintError( "g00=%e g11=%e theta=%e\n",
110 metric.g00, metric.g11, metric.theta );
111 LALPrintError( "dp0=%e dp1=%e\n", dx0, dx1 );
112 }
113
114 /* We have the metric, which is constant. Now we need to place
115 * the templates in the parameter space which is define as follows by
116 * the psi0 and psi3 range:
117 * */
118 bankParams.metric = &metric;
119 bankParams.minimalMatch = coarseIn.mmCoarse;
120 bankParams.x0Min = coarseIn.psi0Min;
121 bankParams.x0Max = coarseIn.psi0Max;
122 bankParams.x1Min = coarseIn.psi3Min;
123 bankParams.x1Max = coarseIn.psi3Max;
124
125 /* Let us define a temporary list of templates. */
126 in.length = 1;
127 in.vectorLength = 2;
128 LALSCreateVectorSequence( status->statusPtr, &tempList, &in );
130
131 /* First we place templates in the psi0/psi3 plane.
132 *
133 * Historically we had two template banks. If gridSpacing is set to
134 * S2BCV, then, the code generates the bank used during S2. This bank
135 * uses a non-oriented square placement in psi0/psi3 plane and the
136 * fcut dimension placement is done BCVRegularFcutBank or BCVFCutBank
137 * function.
138 * If gridSpacing is not S3BCV, then we have the choice between a
139 * square or an hexagonal placement, oriented or not and fcut is placed
140 * using BankFcutS3S4.
141 * */
142 if (coarseIn.gridSpacing != S2BCV)
143 {
144 LALInspiralCreateFlatBankS3S4( status->statusPtr, tempList, &bankParams , coarseIn);
146 }
147 else
148 {
149 LALInspiralCreateFlatBank( status->statusPtr, tempList, &bankParams);
151 }
152
153 *nlist = tempList->length;
154 *list = (InspiralTemplateList *)
155 LALCalloc( *nlist, sizeof(InspiralTemplateList) );
156 if ( ! *list )
157 {
158 ABORT (status, LALINSPIRALBANKH_EMEM, LALINSPIRALBANKH_MSGEMEM);
159 }
160 /* We populate the output.
161 * */
162 for ( j = 0; j < *nlist; ++j )
163 {
164 (*list)[j].params.psi0 = (REAL8) tempList->data[2*j];
165 (*list)[j].params.psi3 = (REAL8) tempList->data[2*j+1];
166 (*list)[j].params.fLower = params.fLower;
167 (*list)[j].params.nStartPad = 0;
168 (*list)[j].params.nEndPad = 0;
169 (*list)[j].params.tSampling = coarseIn.tSampling;
170 (*list)[j].params.distance = 1.;
171 (*list)[j].params.signalAmplitude = 1.;
172 (*list)[j].params.approximant = BCV;
173 (*list)[j].params.massChoice = psi0Andpsi3;
174 (*list)[j].params.order = LAL_PNORDER_TWO;
175 (*list)[j].metric = metric;
176 (*list)[j].params.alpha = coarseIn.alpha;
177 }
178 nlistOld = *nlist;
179
180 /* Once the psi0/psi3 plane is populated, for each coordinate, we
181 * populate along the fcut dimension. Again, for historical reason,
182 * we call one of the following functions which slightly differs
183 * from each other (see documentation).
184 * */
185
186 /* If coarseIn.lowGM == - 1 then LowGM is unphysical. Hence, we use a
187 * Regular grid in cutoff frequency which is independant of LowGM or
188 * HighGM and which lays between Flower and Fsampling/2. If
189 * coarseIn.lowGM != -1 then, we populate between two frequencyies
190 * defines by low and high GM. (i.e lowGM = 6 means fLSO).
191 * */
192 if (coarseIn.gridSpacing != S2BCV)
193 {
195 list, nlist, coarseIn);
197 }
198 else if (coarseIn.LowGM == -1)
199 {
201 list, nlist, coarseIn);
203 }
204 else
205 {
206 LALInspiralBCVFcutBank( status->statusPtr,
207 list, nlist, coarseIn);
209 }
210
211 if ( lalDebugLevel & LALINFO )
212 {
214 "LALInspiralBCVBank: template numbers before %d and after %d calling LALInspiralBCVBank\n",
215 nlistOld, *nlist );
216 }
217
218 LALSDestroyVectorSequence( status->statusPtr, &tempList );
220
222 RETURN( status );
223}
224
225
226/**
227 * The code expects <tt>list->vectorLength=2</tt> and allocates just the
228 * requisite amount of memory to \c list and returns the number
229 * of grid points in <tt>list->length</tt>.
230 * The data points <tt>list->data[2j]</tt>, <tt>j=1,2, ... list->length</tt>,
231 * contain the \f$x_0\f$-coordinates of the grid and data points <tt>list->data[2j+1]</tt>,
232 * contain the \f$x_1\f$-coordinates of the grid.
233 *
234 * ### Description ###
235 *
236 * Given the \c metric and the \c minimalMatch this routine calls
237 * <tt>bank/LALInspiralUpdateParams</tt> to get the spacings in user coordinates (which are
238 * not necessarily the eigen-directions) and lays a uniform grid of templates in
239 * the range specified in (<tt>bankParams->x0Min</tt>, <tt>bankParams->x0Max</tt>) and
240 * (<tt>bankParams->x1Min</tt>, <tt>bankParams->x1Max</tt>).
241 *
242 * ### Algorithm ###
243 *
244 * The algorithm to lay templates is as follows: Given the increments \f$Dx_0\f$ and
245 * \f$Dx_1\f$ found from calling <tt>bank/LALInspiralUpdateParams</tt> lay a rectangular
246 * grid in the space of \f$(x_0, x_1).\f$
247 *
248 * <tt>
249 * <ul>
250 * <li> \f$x_1 = x_1^{\mathrm{min}}\f$
251 * <li> do while (\f$x_1 \le x_1^{\mathrm{max}}\f$)<br>
252 * <ul>
253 * <li> \f$x_0 = x_0^{\mathrm{min}}\f$
254 * <li> do while (\f$x_0 \le x_0^{\mathrm{max}}\f$)<br>
255 * <ul>
256 * <li> Add (\f$x_0, x_1\f$) to list
257 * <li> numTemplates++
258 * <li> Increment \f$x_0: \; x_0 = x_0 + Dx_0\f$
259 * </ul>
260 * <li> Increment \f$x_1: \; x_1 = x_1 + Dx_1\f$
261 * </ul>
262 * </ul>
263 * </tt>
264 */
265void
267 LALStatus *status, /**< LAL status pointer */
268 REAL4VectorSequence *list, /**< [out] an array containing the template bank parameters */
269 InspiralBankParams *bankParams /**< [in] It is necessary and sufficient to input
270 * the eigenvalues of the metric and the angle between the \f$x_0\f$ axis and the
271 * semi-major axis of the ambiguity ellipse, that is,
272 * <tt>bankParams.metric.g00, bankParams.metric.g11, bankParams.metric.theta</tt>,
273 * the minimal match, <tt>bankParams.minimalMatch</tt> and the range of the two
274 * coordinates over which templates must be chosen:
275 * (<tt>bankParams->x0Min</tt>, <tt>bankParams->x0Max</tt>) and
276 * (<tt>bankParams->x1Min</tt>, <tt>bankParams->x1Max</tt>)
277 */
278 )
279{
280 InspiralMetric *metric;
281 REAL8 minimalMatch;
282 REAL8 x0, x1;
283 UINT4 nlist = 0;
284
287
288 /* From the knowledge of the metric and the minimal match find the */
289 /* constant increments bankParams->dx0 and bankParmams->dx1 */
290 metric = bankParams->metric;
291 minimalMatch = bankParams->minimalMatch;
293 bankParams, *metric, minimalMatch );
295
296 /* Construct the template bank */
297 for (x1 = bankParams->x1Min; x1 <= bankParams->x1Max; x1 += bankParams->dx1)
298 {
299 for (x0 = bankParams->x0Min; x0 <= bankParams->x0Max; x0 += bankParams->dx0)
300 {
301 UINT4 ndx = 2 * nlist;
302 list->data = (REAL4 *) LALRealloc( list->data, (ndx+2) * sizeof(REAL4) );
303 if ( !list->data )
304 {
305 ABORT(status, LALINSPIRALBANKH_EMEM, LALINSPIRALBANKH_MSGEMEM);
306 }
307 list->data[ndx] = x0;
308 list->data[ndx + 1] = x1;
309 ++nlist;
310 }
311 }
312
313 list->length = nlist;
314
316 RETURN (status);
317}
318
319/**
320 * \brief Given a grid of templates with distinct values of \f$(\psi_0, \psi_3)\f$
321 * this routine returns a new grid in which every template has \c numFcutTemplates
322 * partners differing from one another in the ending frequency <tt>fendBCV</tt>.
323 *
324 * A call to this function should be preceeded by a call to LALInspiralCreateFlatBank()
325 * or a similar function, that gives a grid in \f$(\psi_0, \psi_3)\f$ space.
326 *
327 * ### Description ###
328 *
329 * A lattice of templates for BCV models should include,
330 * in addition to the values of \f$(\psi_0, \psi_3)\f$
331 * a range of \f$f_{\mathrm{cut}}\f$ -- the cutoff frequency.
332 * The right approach would be
333 * to compute the metric in the three-dimensional space of
334 * \f$(\psi_0, \psi_3, f_{\mathrm{cut}})\f$ and to choose templates as
335 * dictated by the metric. However, analytic computation of the
336 * metric has not been easy. Therefore, it has become necessary
337 * (at least for the time being) to make alternate choice of
338 * the cutoff frequencies.
339 *
340 * In this routine we implement a simple
341 * choice based on physical grounds: The post-Newtonian models
342 * predict an ending frequency that is larger than, but close to,
343 * the Schwarzschild last-stable orbit frequency
344 * \f$f_{\mathrm{lso}} = (6^{3/2} \pi M )^{-1}\f$ where \f$M\f$ is the total mass,
345 * while the effective one-body model has an ending frequency close
346 * to the light-ring, whose Schwarzschild value is
347 * \f$f_{\mathrm{lr}} = (3^{3/2} \pi M )^{-1}\f$. It is necessary to know
348 * the total mass of the system in both cases. However, not all
349 * pairs of \f$(\psi_0, \psi_3)\f$ can be inverted to get a positive
350 * \f$M\f$ but only when \f$\psi_0 > 0\f$ and \f$\psi_3 < 0\f$. Even then
351 * it is not guaranteed that the symmetric mass ratio will be
352 * less than \f$1/4,\f$ a necessary condition so that the component
353 * masses are found to be real. However, we do not demand that the
354 * symmetric mass ratio is less than a quarter. If the total mass
355 * is non-negative then we compute the \f$(f_{\mathrm{lso}}, f_{\mathrm{lr}})\f$
356 * and choose a user specified \c numFcutTemplates number of
357 * templates with their cutoff frequency <tt>list->fFinal</tt> defined
358 * uniformly spaced in the range \f$[f_{\mathrm{lso}},\ f_{\mathrm{lr}}]\f$.
359 *
360 * Furthermore, this routine discards all templates for which
361 * either the mass is not defined or, when defined, <tt>list->fFinal</tt> is
362 * smaller than the user defined lower frequency cutoff or larger
363 * than the Nyquist frequency of templates.
364 * Thus, the number of templates returned by this routine could
365 * be larger or fewer than the input number of templates.
366 *
367 * ### Algorithm ###
368 *
369 * Given \f$(\psi_0, \psi_3)\f$ one can solve for \f$(M, \eta)\f$ using:
370 * \f{equation}{
371 * M = \frac{-\psi_3}{16 \pi^2 \psi_0},\ \ \eta = \frac{3}{128 \psi_0 (\pi M)^{5/3}}.
372 * \f}
373 * Given the total mass compute the last stable orbit and light-ring frequencies using
374 * \f{equation}{
375 * f_{\mathrm{lso}} = (6^{3/2} \pi M)^{-1},\ \ f_{\mathrm{lr}} = (3^{3/2} \pi M)^{-1}.
376 * \f}
377 * Divide the range \f$(f_{\mathrm{lso}}, f_{\mathrm{lr}})\f$ so as to have \f$n_{\mathrm{cut}}= \mathtt{numFcutTemplates}\f$
378 * templates over this range:
379 * \f{equation}{
380 * df = f_{\mathrm{lr}} \frac {\left( 1 - 2^{-3/2} \right) }{ (n_{\mathrm{cut}} -1) }.
381 * \f}
382 * Next, choose templates at \f$f_k = f_\mathrm{lr} - k \times df,\f$ where \f$k=0, \ldots, n_\mathrm{cut}-1\f$.
383 * Note that by definition \f$f_0 = f_\mathrm{lr}\f$ and \f$f_{n_\mathrm{cut}-1} = f_\mathrm{lso}\f$;
384 * there are exatly \f$n_\mathrm{cut}\f$ templates in the range \f$(f_\mathrm{lso}, f_\mathrm{lr})\f$.
385 * We discard a template if either \f$M\f$ is not defined or if \f$f_\mathrm{cut}\f$ is smaller
386 * than the lower frequency cutoff specified in <tt>list[j]->fLower</tt>.
387 */
388void
390 LALStatus *status, /**< LAL status pointer */
391 InspiralTemplateList **list, /**< [out,in] an array initially containing the template
392 * bank with the values of <tt>list[j]->psi0, list[j]->psi3, list[j]->fLower,</tt> specified,
393 * is replaced on return with a re-sized array specifying also <tt>list->fFinal.</tt>
394 */
395 INT4 *NList, /**< [out,in] the number of templates in the Input bank is replaced by the number of templates in the output bank. */
396 InspiralCoarseBankIn coarseIn /**< UNDOCUMENTED */
397 )
398{
399 UINT4 nf; /* number of layers */
400 UINT4 nlist; /* number of final templates */
401 UINT4 j;
402 UINT4 ndx; /* temporary number of templates */
403 REAL8 frac; /* general variable*/
404 REAL8 fendBCV;
405 REAL4 LowGM;
406 REAL4 HighGM;
407
409
410 nf = coarseIn.numFcutTemplates;
411 ndx = nlist = *NList;
412
413 LowGM = coarseIn.LowGM;
414 HighGM = coarseIn.HighGM;
415
416 /* if we have only one layer, we don't need HighGM.
417 * And default value for LowGM is 3GM*/
418 if ( nf == 1 )
419 {
420 frac = 1;
421 }
422 else
423 {
424 frac = (1.L - 1.L/pow(HighGM/3., 1.5L)) / (nf-1.L);
425 }
426
427 /* for each psi0/psi3 pair, we generate the fcut layers */
428 for ( j = 0; j < nlist; ++j )
429 {
430 UINT4 valid = 0;
431 /* let us get the estimated params.fFinal*/
432 LALPSItoMasses(status, &((*list)[j].params), &valid , LowGM);
433
434 if ( valid )
435 {
436 UINT4 i;
437 REAL8 fMax;
438
439 fMax = (*list)[j].params.fFinal;
440 /* for each fcut layer */
441 for ( i = 0; i < nf; ++i )
442 {
443 fendBCV = fMax * (1.L - (REAL8) i * frac);
444
445 if ( (*list)[j].params.tSampling <= 0 )
446 {
447 ABORT( status, LALINSPIRALBANKH_ESIZE, LALINSPIRALBANKH_MSGESIZE );
448 }
449 /* the fFinal must be > flower and less than Nyquist.
450 * if not, no template is generated. */
451 if ( fendBCV > (*list)[j].params.fLower &&
452 fendBCV < (*list)[j].params.tSampling / 2.0 )
453 {
454 ++ndx;
455
456 *list = (InspiralTemplateList *)
457 LALRealloc( *list, ndx * sizeof(InspiralTemplateList) );
458 if ( ! *list )
459 {
460 ABORT( status, LALINSPIRALBANKH_EMEM, LALINSPIRALBANKH_MSGEMEM );
461 }
462 memset( *list + ndx - 1, 0, sizeof(InspiralTemplate) );
463 (*list)[ndx-1] = (*list)[j];
464 (*list)[ndx-1].params.fFinal = fendBCV;
465 (*list)[ndx-1].metric = (*list)[0].metric;
466 (*list)[ndx-1].nLayer = i;
467 }
468 }
469 }
470 }
471 /**/
472 for ( j = nlist; j < ndx; ++j )
473 {
474 (*list)[j-nlist] = (*list)[j];
475 }
476 /**/
477 *NList = ndx - nlist;
478 *list = LALRealloc( *list, *NList * sizeof(InspiralTemplateList) );
479 if ( ! *list )
480 {
481 ABORT( status, LALINSPIRALBANKH_EMEM, LALINSPIRALBANKH_MSGEMEM );
482 }
483
484 RETURN( status );
485}
486
487/** UNDOCUMENTED */
488void
492 UINT4 *valid,
493 REAL4 HighGM
494 )
495{
496
499
500 if ( params->psi0 <= 0.L || params->psi3 >= 0.L )
501 {
502 *valid = 0;
503 }
504 else
505 {
506 REAL8 totalMass;
507 REAL8 eta;
508 /*we estimate the total mass and then fFinal*/
509 params->totalMass = -params->psi3/(16.L*LAL_PI * LAL_PI * params->psi0);
510 eta = params->eta =
511 3.L/(128.L * params->psi0 * pow(LAL_PI*params->totalMass, (5./3.)));
512 totalMass = params->totalMass;
513 params->fFinal = 1.L/( LAL_PI * pow(HighGM, 1.5L) * params->totalMass );
514 params->totalMass /= LAL_MTSUN_SI;
515 *valid = 1;
516
517#if 0
518 if (params->eta > 0.25L)
519 {
520 *valid = 0;
521 }
522 else
523 {
526 *valid = 1;
527 }
528#endif
529
530 params->t0 = 5.0 / ( 256.0*eta*pow(totalMass, (5./3.)) *
531 pow(LAL_PI * params->fLower, (8./3.)));
532 params->t3 = LAL_PI/(8.0*eta*pow(totalMass, (2./3.)) *
533 pow(LAL_PI * params->fLower, (5./3.)));
534 }
536 RETURN (status);
537}
538
539/** UNDOCUMENTED */
540void
544 INT4 *NList,
545 InspiralCoarseBankIn coarseIn
546 )
547
548{
549 UINT4 nf; /* number of layers */
550 UINT4 nlist; /* number of final templates */
551 UINT4 j;
552 UINT4 ndx; /* temporary number of templates */
553 REAL8 frac; /* general variable*/
554 REAL8 fendBCV;
555 REAL4 LowGM;
556 REAL4 HighGM;
557
559
560 nf = coarseIn.numFcutTemplates;
561 ndx = nlist = *NList;
562 LowGM = 3.;
563 HighGM = coarseIn.HighGM;
564
565 /* for each template, we get the fcut layers*/
566 for ( j = 0; j < nlist; ++j )
567 {
568 UINT4 valid = 0;
570 &((*list)[j].params), &valid , LowGM);
571
572 if (valid)
573 {
574 UINT4 i;
575 REAL8 fMax;
576
577 fMax = (*list)[j].params.fFinal;
578 if ( (*list)[j].params.tSampling <= 0 )
579 {
580 ABORT( status, LALINSPIRALBANKH_ESIZE, LALINSPIRALBANKH_MSGESIZE );
581 }
582 /* the user might request only one layer */
583 if (nf == 1)
584 {
585 frac = 1;
586 }
587 else
588 {
589 frac = (1.L - 1.L/pow(HighGM/3., 1.5L)) / (nf-1.L);
590 }
591
592 /* sometimes since fMin is greater than the nyquist frequency, there
593 * is no template generated. This is not acceptable. We need at least
594 * one frequency at the nyquist frequency otherwise low masses
595 * systems are missed. */
596 if (((fMax * (1.L - (REAL4) (nf-1) * frac)) >= (*list)[j].params.tSampling/2.0))
597 {
598 fMax = (*list)[j].params.tSampling/2.0 - 1. ;
599 frac = -1;
600 }
601
602 /*Similarly, for high masses. */
603 /*if (((fMax * (1.L - (REAL4) (nf-1) * frac)) <= (*list)[j].params.fLower * 1.5))
604 {
605 fMax = (*list)[j].params.fLower * 1.5 ;
606 }
607 */
608 for (i=0; i<nf; i++)
609 {
610 fendBCV = fMax * (1.L - (REAL4) i * frac);
611 /* we search for valid expression of fendBCV and populate the bank */
612 if ( fendBCV >= (*list)[j].params.fLower * 1.5 &&
613 fendBCV < (*list)[j].params.tSampling / 2.0 )
614 {
615 ++ndx;
616 *list = (InspiralTemplateList *)
617 LALRealloc( *list, ndx * sizeof(InspiralTemplateList) );
618 if ( ! *list )
619 {
620 ABORT( status, LALINSPIRALBANKH_EMEM, LALINSPIRALBANKH_MSGEMEM );
621 }
622 memset( *list + ndx - 1, 0, sizeof(InspiralTemplate) );
623 (*list)[ndx-1] = (*list)[j];
624 (*list)[ndx-1].params.fFinal = fendBCV;
625 (*list)[ndx-1].metric = (*list)[0].metric;
626 (*list)[ndx-1].nLayer = i;
627 }
628 }
629 }
630 }
631 for ( j = nlist; j < ndx; ++j )
632 {
633 (*list)[j-nlist] = (*list)[j];
634 }
635 *NList = ndx - nlist;
636 *list = LALRealloc( *list, *NList * sizeof(InspiralTemplateList) );
637 if ( ! *list )
638 {
639 ABORT( status, LALINSPIRALBANKH_EMEM, LALINSPIRALBANKH_MSGEMEM );
640 }
641
642
643 RETURN( status );
644}
645
646/** UNDOCUMENTED */
647void
651 INT4 *NList,
652 InspiralCoarseBankIn coarseIn
653 )
654
655{
656 /* no restriction of physical masses.
657 * And regular layer of templates in the Frequency dimension */
658 UINT4 i,nf, nlist, j, ndx;
659 REAL8 fendBCV;
660
662
663 nf = coarseIn.numFcutTemplates;
664 ndx = nlist = *NList;
665
666 for ( j = 0; j < nlist; ++j )
667 {
668 for ( i = 1; i <=nf; ++i )
669 {
670 fendBCV = (*list)[j].params.fLower
671 + i * ((*list)[j].params.tSampling/2.0 - (*list)[j].params.fLower) / nf ;
672 ++ndx;
673
674 *list = (InspiralTemplateList *)
675 LALRealloc( *list, ndx * sizeof(InspiralTemplateList) );
676 if ( ! *list )
677 {
678 ABORT( status, LALINSPIRALBANKH_EMEM, LALINSPIRALBANKH_MSGEMEM );
679 }
680 memset( *list + ndx - 1, 0, sizeof(InspiralTemplate) );
681 (*list)[ndx-1] = (*list)[j];
682 (*list)[ndx-1].params.fFinal = fendBCV;
683 (*list)[ndx-1].metric = (*list)[0].metric;
684 (*list)[ndx-1].nLayer = i;
685 }
686 }
687
688
689 for ( j = nlist; j < ndx; ++j )
690 {
691 (*list)[j-nlist] = (*list)[j];
692 }
693
694 *NList = ndx - nlist;
695 *list = LALRealloc( *list, *NList * sizeof(InspiralTemplateList) );
696 if ( ! *list )
697 {
698 ABORT( status, LALINSPIRALBANKH_EMEM, LALINSPIRALBANKH_MSGEMEM );
699 }
700
701 RETURN( status );
702}
703
704/** UNDOCUMENTED */
705void
709 UINT4 *valid,
710 REAL4 lightring
711 )
712{
713
716
717 if ( params->psi0 <= 0.L || params->psi3 >= 0.L )
718 {
719 *valid = 0;
720 }
721 else
722 {
723 params->totalMass = -params->psi3/(16.L*LAL_PI * LAL_PI * params->psi0);
724 params->totalMass = params->totalMass * 2. ; /* The factor 2 is purely empiricail and
725 comes from simulaitons. ?It seems indeed
726 tyhat the relation between psi0 and psi3
727 which gives the total mass is not really
728 suitable. Ususally, the total mass is
729 twice as much as the estimated one.
730 */
731 params->fFinal = 1.L/( LAL_PI * pow(lightring, 1.5L) * params->totalMass );
732 params->totalMass /= LAL_MTSUN_SI; /* it it used later ? */
733
734 *valid = 1;
735 }
736
738 RETURN (status);
739}
740
741/** UNDOCUMENTED */
742void
746 InspiralBankParams *bankParams,
747 InspiralCoarseBankIn coarseIn
748 )
749
750{
751 InspiralMetric *metric;
752 REAL8 minimalMatch;
753 REAL8 x0, x1, dx1=0, dx0=0, x=0, y=0;
754 UINT4 nlist = 0;
755 INT4 layer = 1;
756 INT4 valid = -1;
757 REAL4 xp[8] = {3000, 40000, 100000, 300000, 550000, 550000, 250000,3000};
758 REAL4 yp[8] = {0, -3000, -3000, -2000, -1500, -300, -300, 0};
759 INT4 npol = 8;
760
761
764
765
766 /* From the knowledge of the metric and the minimal match
767 find the constant increments bankParams->dx0 and
768 bankParmams->dx1 */
769 metric = bankParams->metric;
770 minimalMatch = bankParams->minimalMatch;
771
772 switch (coarseIn.gridSpacing){
773 case HybridHexagonal:
774 case S2BCV:
775 case Hexagonal:
776 dx0 = sqrt(2.L * (1.L - minimalMatch)/metric->g00 );
777 dx1 = sqrt(2.L * (1.L - minimalMatch)/metric->g11 );
778 dx0 *=3./2./sqrt(2.);
779 dx1 *=sqrt(3./2.);
780 break;
781 case Square:
782 dx0 = sqrt(2.L * (1.L - minimalMatch)/metric->g00 );
783 dx1 = sqrt(2.L * (1.L - minimalMatch)/metric->g11 );
784 break;
787 bankParams, *metric, minimalMatch );
789 dx0 = bankParams->dx0 * 3./2./sqrt(2.);
790 dx1 = bankParams->dx1 * sqrt(3./2.);
791 break;
792
795 bankParams, *metric, minimalMatch );
797 dx0 = bankParams->dx0;
798 dx1 = bankParams->dx1;
799 break;
800 }
801
802
803 switch (coarseIn.gridSpacing)
804 {
805 case HybridHexagonal:
806 case S2BCV:
807 case Hexagonal:
809
810 /* x1==psi3 and x0==psi0 */
811 for (x1 = bankParams->x1Min -1e6; x1 <= bankParams->x1Max + 1e6; x1 += dx1)
812 {
813 layer++;
814 for (x0 = bankParams->x0Min - 1e6 +dx0/2.*(layer%2); x0 <= bankParams->x0Max+1e6; x0 += dx0 )
815 {
816 UINT4 ndx = 2 * nlist;
817 if ( coarseIn.gridSpacing == Hexagonal)
818 {
819 x = x0 *cos(metric->theta) + sin(metric->theta)* x1;
820 y = x0 *sin(metric->theta) - cos(metric->theta)* x1;
821 }
822 else
823 {
824 x = x0;
825 y = x1;
826 }
827
828 if ( (x > bankParams->x0Min -dx0/2.) && (y < bankParams->x1Max + dx1/2.) &&
829 (x < bankParams->x0Max +dx0/2.) && (y > bankParams->x1Min - dx1/2.))
830 {
831
832 if (coarseIn.insidePolygon == True)
833 {
834 LALInsidePolygon(status->statusPtr, xp, yp, npol, x, y, &valid);
835 }
836 else
837 {
838 LALExcludeTemplate(status->statusPtr, &valid, bankParams, x, y);
839 }
840
841
842 if (valid == 1)
843 {
844 list->data = (REAL4 *) LALRealloc( list->data, (ndx+2) * sizeof(REAL4) );
845 if ( !list->data )
846 {
847 ABORT(status, LALINSPIRALBANKH_EMEM, LALINSPIRALBANKH_MSGEMEM);
848 }
849 list->data[ndx] = x;
850 list->data[ndx + 1] = y;
851 ++nlist;
852 }
853 }
854 }
855 }
856 break;
857
858 case Square:
860
861 /* !! dx1 and dx0 are computed in a different way de[pending on the
862 value of BANKGRId */
863 for (x1 = bankParams->x1Min -1e6; x1 <= bankParams->x1Max + 1e6; x1 += dx1)
864 {
865 for (x0 = bankParams->x0Min - 1e6 ; x0 <= bankParams->x0Max+1e6; x0 += dx0 )
866 {
867 UINT4 ndx = 2 * nlist;
868
869 if (coarseIn.gridSpacing == Square)
870 {
871 x = x0 *cos(metric->theta) + sin(metric->theta)* x1 ;
872 y = x0 *sin(metric->theta) - cos(metric->theta)* x1;
873 }
874 else if (coarseIn.gridSpacing == SquareNotOriented)
875 {
876 x = x0;
877 y = x1;
878 }
879 if ( (x > bankParams->x0Min - dx0/2.) && (y < bankParams->x1Max + dx1/2.) &&
880 (x < bankParams->x0Max + dx0/2.) && (y > bankParams->x1Min - dx1/2.))
881
882 {
883
884 if (coarseIn.insidePolygon == True)
885 {
886 LALInsidePolygon(status->statusPtr, xp, yp, npol, x, y, &valid);
887 }
888 else
889 {
890 LALExcludeTemplate(status->statusPtr, &valid, bankParams, x, y);
891 }
892 if (valid)
893 {
894 list->data = (REAL4 *) LALRealloc( list->data, (ndx+2) * sizeof(REAL4) );
895 if ( !list->data )
896 {
897 ABORT(status, LALINSPIRALBANKH_EMEM, LALINSPIRALBANKH_MSGEMEM);
898 }
899 list->data[ndx] = x;
900 list->data[ndx + 1] = y;
901 ++nlist;
902 }
903 }
904 }
905 }
906 break;
907 }
908
909 list->length = nlist;
910
912 RETURN (status);
913}
914
915
916/**
917 * Thomas: 31 Aug 2006. This function is redundant with the polygon fit.
918 * It was design for BBH and therefore had tight boundary. For a more general
919 * purpose, I extended the range to generous values
920 */
921void
924 INT4 *valid,
925 InspiralBankParams UNUSED *bankParams,
926 REAL4 x,
927 REAL4 y)
928{
929 REAL4 psi0Int = 2500000.;
930 REAL4 psi3Int = -10000.;
931
934
935 if (x > psi0Int && y < psi3Int )
936 {
937 *valid = 0 ;
938 }
939 else
940 {
941 *valid = 1;
942 }
943
945 RETURN (status);
946}
947/** @} */
void LALInspiralParameterCalc(LALStatus *status, InspiralTemplate *params)
#define LALRealloc(p, n)
#define LALCalloc(m, n)
#define ABORT(statusptr, code, mesg)
#define CHECKSTATUSPTR(statusptr)
#define ATTATCHSTATUSPTR(statusptr)
#define ASSERT(assertion, statusptr, code, mesg)
#define DETATCHSTATUSPTR(statusptr)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
double i
void LALSDestroyVectorSequence(LALStatus *status, REAL4VectorSequence **vectorSequence)
void LALSCreateVectorSequence(LALStatus *status, REAL4VectorSequence **vectorSequence, CreateVectorSequenceIn *vSeqParams)
#define LAL_PI
#define LAL_MTSUN_SI
double REAL8
uint32_t UINT4
int32_t INT4
float REAL4
#define lalDebugLevel
LALINFO
int LALPrintError(const char *fmt,...)
@ psi0Andpsi3
BCV parameters and .
Definition: LALInspiral.h:187
void LALInspiralBCVRegularFcutBank(LALStatus *status, InspiralTemplateList **list, INT4 *NList, InspiralCoarseBankIn coarseIn)
UNDOCUMENTED.
void LALInspiralCreateBCVBank(LALStatus *status, InspiralTemplateList **list, INT4 *nlist, InspiralCoarseBankIn coarseIn)
Lay a flat grid of BCV templates in the user specified range of the parameters in coarseIn structure...
void LALPSItoMasses(LALStatus *status, InspiralTemplate *params, UINT4 *valid, REAL4 HighGM)
UNDOCUMENTED.
void LALInspiralBCVBankFcutS3S4(LALStatus *status, InspiralTemplateList **list, INT4 *NList, InspiralCoarseBankIn coarseIn)
UNDOCUMENTED.
void LALInspiralCreateFlatBank(LALStatus *status, REAL4VectorSequence *list, InspiralBankParams *bankParams)
The code expects list->vectorLength=2 and allocates just the requisite amount of memory to list and r...
void LALEmpiricalPSItoMassesConversion(LALStatus *status, InspiralTemplate *params, UINT4 *valid, REAL4 lightring)
UNDOCUMENTED.
void LALInspiralCreateFlatBankS3S4(LALStatus *status, REAL4VectorSequence *list, InspiralBankParams *bankParams, InspiralCoarseBankIn coarseIn)
UNDOCUMENTED.
void LALExcludeTemplate(LALStatus *status, INT4 *valid, InspiralBankParams UNUSED *bankParams, REAL4 x, REAL4 y)
Thomas: 31 Aug 2006.
void LALInspiralBCVFcutBank(LALStatus *status, InspiralTemplateList **list, INT4 *NList, InspiralCoarseBankIn coarseIn)
Given a grid of templates with distinct values of this routine returns a new grid in which every tem...
void LALInsidePolygon(LALStatus *status, REAL4 *inputx, REAL4 *inputy, INT4 n, REAL4 x, REAL4 y, INT4 *valid)
Module to check whether a point with coordinates (x,y) is inside a polygon defined by the vectors (vx...
void LALInspiralUpdateParams(LALStatus *status, InspiralBankParams *bankParams, InspiralMetric metric, REAL8 minimalMatch)
Function to update the parameters used in creating a coarse bank based on a square lattice.
void LALInspiralComputeMetricBCV(LALStatus *status, InspiralMetric *metric, REAL8FrequencySeries *psd, InspiralTemplate *params)
UNDOCUMENTED.
#define LALINSPIRALBANKH_ESIZE
Invalid input range.
#define LALINSPIRALBANKH_EMEM
Memory allocation failure.
@ SquareNotOriented
UNDOCUMENTED.
@ Square
UNDOCUMENTED.
@ S2BCV
UNDOCUMENTED.
@ HexagonalNotOriented
UNDOCUMENTED.
@ HybridHexagonal
UNDOCUMENTED.
@ Hexagonal
UNDOCUMENTED.
@ True
BCV
LAL_PNORDER_TWO
list y
x
This is a structure needed in the inner workings of the LALInspiralCreateCoarseBank code.
REAL8 x1Max
maximum value of the second coordinate as defined by the search region
REAL8 x1Min
minimum value of the second coordinate as defined by the search region
REAL8 x0Min
minimum value of the first coordinate as defined by the search region
REAL8 dx0
increment in the x0-direction
InspiralMetric * metric
pointer to the metric at the current location
REAL8 dx1
increment in the x1-direction
REAL8 x0Max
maximum value of the first coordinate as defined by the search region
REAL8 minimalMatch
UNDOCUMENTED.
Input for choosing a template bank.
REAL8 tSampling
Sampling rate.
GridSpacing gridSpacing
Type of gridspacing required.
REAL8FrequencySeries shf
Frequency series containing the PSD.
REAL8 psi3Min
minimum value of the parameter
InsidePolygonEnum insidePolygon
UNDOCUMENTED.
REAL8 fUpper
Upper frequency cutoff.
REAL8 fLower
Lower frequency cutoff.
REAL4 HighGM
UNDOCUMENTED.
UINT4 numFcutTemplates
Number of templates required in the fCut (upper cutoff) dimension and the value of upper and lower cu...
REAL8 psi3Max
maximum value of the parameter
REAL8 mmCoarse
Coarse grid minimal match.
REAL4 LowGM
UNDOCUMENTED.
REAL8 alpha
the BCV amplitude correction parameter
REAL8 psi0Min
minimum value of the parameter
REAL8 psi0Max
maximum value of the parameter
Structure to store metric at various points the signal manifold.
REAL8 g11
11-component of the diagonalised metric
REAL8 G01
01-component of the metric in coordinates
REAL8 g00
00-component of the diagonalised metric
REAL8 theta
Angle from tau0 to semi-major axis of the ellipse.
REAL8 G00
00-component of the metric in coordinates
REAL8 G11
11-component of the metric in coordinates
The inspiral waveform parameter structure containing information about the waveform to be generated.
Definition: LALInspiral.h:205
A grid of inspiral templates (ie a template list).
LALDict * params