Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-8a6b96f
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
XLALComputeAMTest.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010 Reinhard Prix
3 * Copyright (C) 2006 John Whelan, Reinhard Prix
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with with program; see the file COPYING. If not, write to the
17 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
18 * MA 02110-1301 USA
19 */
20
21#include <config.h>
22
23#include <math.h>
24#include <sys/times.h>
25
26#include <gsl/gsl_math.h>
27
28#include <lal/ComputeFstat.h>
29#include <lal/LALgetopt.h>
30#include <lal/LALBarycenter.h>
31#include <lal/LALInitBarycenter.h>
32#include <lal/AVFactories.h>
33#include <lal/SinCosLUT.h>
34
35/**
36 * \author Reinhard Prix, John Whelan
37 * \file
38 * \ingroup LALComputeAM_h
39 *
40 * \brief Test for XLALComputeAMCoeffs() and XLALComputeMultiAMCoeffs() by
41 * comparison with the old LAL functions old_LALGetAMCoeffs() and old_LALGetMultiAMCoeffs().
42 *
43 * Note, we run a comparison only for the 2-IFO multiAM functions XLALComputeMultiAMCoeffs()
44 * comparing it to old_LALGetMultiAMCoeffs() [combined with XLALWeightMultiAMCoeffs()],
45 * as this excercises the 1-IFO functions as well.
46 *
47 * Sky-location is picked at random each time, which allows a minimal
48 * Monte-Carlo validation by simply running this script repeatedly.
49 *
50 */
51
52/* ----- internal prototypes ---------- */
53static int XLALCompareMultiAMCoeffs( MultiAMCoeffs *multiAM1, MultiAMCoeffs *multiAM2, REAL8 tolerance );
54static void old_LALGetMultiAMCoeffs( LALStatus *, MultiAMCoeffs **multiAMcoef, const MultiDetectorStateSeries *multiDetStates, SkyPosition pos );
55static void old_LALNewGetAMCoeffs( LALStatus *, AMCoeffs *coeffs, const DetectorStateSeries *DetectorStates, SkyPosition skypos );
56
57#define LALCOMPUTEAMH_ENULL 7
58#define LALCOMPUTEAMH_EINPUT 8
59#define LALCOMPUTEAMH_ENONULL 9
60#define LALCOMPUTEAMH_EMEM 10
61
62#define LALCOMPUTEAMH_MSGENULL "Arguments contained an unexpected null pointer"
63#define LALCOMPUTEAMH_MSGEINPUT "Invalid input"
64#define LALCOMPUTEAMH_MSGENONULL "Output pointer is non-NULL"
65#define LALCOMPUTEAMH_MSGEMEM "Out of memory. Bad."
66
67/* ----- function definitions ---------- */
68int main( int argc, char *argv[] )
69{
70 int opt; /* Command-line option. */
71
72 UINT4 numIFOs = 2;
73 const CHAR *sites[2] = {"H1", "V1"};
74
75 UINT4 startTime = 714180733;
76 UINT4 duration = 180000; /* 50 hours */
77 UINT4 Tsft = 1800; /* assume 30min SFTs */
78
79 REAL8 tolerance = 2e-6; /* same algorithm, should be basically identical results */
80
81 char earthEphem[] = TEST_PKG_DATA_DIR "earth00-40-DE405.dat.gz";
82 char sunEphem[] = TEST_PKG_DATA_DIR "sun00-40-DE405.dat.gz";
83
84 UINT4 numChecks = 1; /* Number of times to check */
85
86 /* read user input */
87
88 while ( ( opt = LALgetopt( argc, argv, "n:qv:" ) ) != -1 ) {
89 switch ( opt ) {
90 case 'v': /* set lalDebugLevel */
91 break;
92 case 'n': /* number of times to check */
93 numChecks = atoi( LALoptarg );
94 break;
95 }
96 }
97
98 /* init random-generator */
99 struct tms buf;
100 UINT4 seed = times( &buf );
101 srand( seed );
102 XLALPrintInfo( "%s: seed used = %d\n", __func__, seed );
103
104 /* ----- init ephemeris ----- */
106 if ( ( edat = XLALInitBarycenter( earthEphem, sunEphem ) ) == NULL ) {
107 XLALPrintError( "%s: XLALInitBarycenter() failed with xlalErrno = %d\n", __func__, xlalErrno );
108 return XLAL_EFAILED;
109 }
110
111 /* ----- init detector info ---------- */
112 UINT4 X;
113 MultiLALDetector multiDet;
114 multiDet.length = numIFOs;
115 for ( X = 0; X < numIFOs; X ++ ) {
116 const LALDetector *site;
117 if ( ( site = XLALGetSiteInfo( sites[X] ) ) == NULL ) {
118 XLALPrintError( "%s: Failed to get site-info for detector '%s'\n", __func__, sites[X] );
119 return XLAL_EFAILED;
120 }
121 multiDet.sites[X] = ( *site ); /* copy! */
122 }
123
124 /* ----- init multi-IFO timestamps vector ---------- */
125 UINT4 numSteps = ( UINT4 ) ceil( duration / Tsft );
127 if ( ( multiTS = XLALCalloc( 1, sizeof( *multiTS ) ) ) == NULL ) {
129 }
130 multiTS->length = numIFOs;
131 if ( ( multiTS->data = XLALCalloc( numIFOs, sizeof( *multiTS->data ) ) ) == NULL ) {
133 }
134 for ( X = 0; X < numIFOs; X ++ ) {
135 if ( ( multiTS->data[X] = XLALCreateTimestampVector( numSteps ) ) == NULL ) {
136 XLALPrintError( "%s: XLALCreateTimestampVector(%d) failed.\n", __func__, numSteps );
137 return XLAL_EFAILED;
138 }
139 multiTS->data[X]->deltaT = Tsft;
140
141 UINT4 i;
142 for ( i = 0; i < numSteps; i ++ ) {
143 UINT4 ti = startTime + i * Tsft;
144 multiTS->data[X]->data[i].gpsSeconds = ti;
145 multiTS->data[X]->data[i].gpsNanoSeconds = 0;
146 } /* for i < numSteps */
147
148 } /* for X < numIFOs */
149
150 /* ---------- compute multi-detector states -------------------- */
151 MultiDetectorStateSeries *multiDetStates;
152 if ( ( multiDetStates = XLALGetMultiDetectorStates( multiTS, &multiDet, edat, 0.5 * Tsft ) ) == NULL ) {
153 XLALPrintError( "%s: XLALGetMultiDetectorStates() failed.\n", __func__ );
154 return XLAL_EFAILED;
155 }
158
159 /* ========== MAIN LOOP: N-trials of comparisons XLAL <--> LAL multiAM functions ========== */
160 while ( numChecks-- ) {
162
163 /* ----- pick skyposition at random ----- */
164 SkyPosition skypos;
165 skypos.longitude = LAL_TWOPI * ( 1.0 * rand() / ( RAND_MAX + 1.0 ) ); /* uniform in [0, 2pi) */
166 skypos.latitude = LAL_PI_2 - acos( 1 - 2.0 * rand() / RAND_MAX ); /* sin(delta) uniform in [-1,1] */
168
169 MultiNoiseWeights *weights = NULL; /* for now we only deal with unit-weights case */
170
171 /* ----- compute multiAM using old LAL function ----- */
172 MultiAMCoeffs *multiAM_LAL = NULL;
173 old_LALGetMultiAMCoeffs( &status, &multiAM_LAL, multiDetStates, skypos );
174 if ( status.statusCode ) {
175 XLALPrintError( "%s: old_LALGetMultiAMCoeffs() failed with statusCode = %d : %s\n", __func__, status.statusCode, status.statusDescription );
176 return XLAL_EFAILED;
177 }
178 if ( XLALWeightMultiAMCoeffs( multiAM_LAL, weights ) != XLAL_SUCCESS ) {
179 XLALPrintError( "%s: XLALWeightMultiAMCoeffs() failed with xlalErrno = %d\n", __func__, xlalErrno );
180 return XLAL_EFAILED;
181 }
182
183
184 /* ----- compute multiAM using XLAL function ----- */
185 MultiAMCoeffs *multiAM_XLAL;
186 if ( ( multiAM_XLAL = XLALComputeMultiAMCoeffs( multiDetStates, weights, skypos ) ) == NULL ) {
187 XLALPrintError( "%s: XLALComputeMultiAMCoeffs() failed with xlalErrno = %d\n", __func__, xlalErrno );
188 return XLAL_EFAILED;
189 }
190
191 /* now run comparison */
192 if ( XLALCompareMultiAMCoeffs( multiAM_XLAL, multiAM_LAL, tolerance ) != XLAL_SUCCESS ) {
193 XLALPrintError( "%s: comparison between multiAM_XLAL and multiAM_LAL failed.\n", __func__ );
194 return XLAL_EFAILED;
195 }
196
197 /* free memory created inside this loop */
198 XLALDestroyMultiAMCoeffs( multiAM_LAL );
199 XLALDestroyMultiAMCoeffs( multiAM_XLAL );
200
201 } /* for numChecks */
202
203 /* we're done: free memory */
204 XLALDestroyMultiDetectorStateSeries( multiDetStates );
205
207
208 printf( "OK. All tests passed successfully\n\n" );
209
210 return 0; /* OK */
211
212} /* main() */
213
214/*
215 * Comparison function for two multiAM vectors, return success or failure for given tolerance.
216 * we compare avg() and max of |a1_i - a2_i|^2 and |b1_i - b2_i|^2 respectively,
217 * and error in |A1 - A2|, |B1 - B2|, |C1 - C2|.
218 * These numbers are typically ~ O(1), so we simply compare these absolute errors to the tolerance.
219 *
220 */
221int
223{
224 /* check input */
225 if ( !multiAM1 || !multiAM2 || tolerance <= 0 ) {
226 XLALPrintError( "%s: invalid NULL input or non-positive tolerance\n", __func__ );
228 }
229
230 UINT4 numDet = multiAM1->length;
231 if ( numDet != multiAM2->length ) {
232 XLALPrintError( "%s: number of detectors differ multiAM1 = %d, multiAM2 = %d\n", __func__, multiAM1->length, multiAM2->length );
234 }
235
236 UINT4 X;
237 REAL8 maxerr_ab = 0, avgerr_ab = 0;
238 UINT4 numTerms = 0;
239 for ( X = 0; X < numDet; X ++ ) {
240 UINT4 numSteps = multiAM1->data[X]->a->length;
241 if ( numSteps != multiAM2->data[X]->a->length ) {
242 XLALPrintError( "%s: number of timesteps differ multiAM1[%d]->a = %d, multiAM2[%d]->a = %d\n", __func__, X, multiAM1->data[X]->a->length, X, multiAM2->data[X]->a->length );
244 }
245 if ( numSteps != multiAM1->data[X]->b->length || numSteps != multiAM2->data[X]->b->length ) {
246 XLALPrintError( "%s: number of timesteps differ multiAM1[%d]->b = %d, multiAM2[%d]->b = %d\n", __func__, X, multiAM1->data[X]->b->length, X, multiAM2->data[X]->b->length );
248 }
249
250 UINT4 i;
251 for ( i = 0; i < numSteps; i ++ ) {
252 REAL8 err_a = fabs( multiAM1->data[X]->a->data[i] - multiAM2->data[X]->a->data[i] );
253 REAL8 err_b = fabs( multiAM1->data[X]->b->data[i] - multiAM2->data[X]->b->data[i] );
254
255 if ( err_a > maxerr_ab ) {
256 maxerr_ab = err_a;
257 }
258 if ( err_b > maxerr_ab ) {
259 maxerr_ab = err_b;
260 }
261
262 avgerr_ab += err_a + err_b;
263 numTerms += 1;
264
265 } /* for i < numSteps */
266
267 } /* for X < numDet */
268
269 avgerr_ab /= ( 2.0 * numTerms );
270
271 /* now compute absolute maximal error in AntennaPattern matrix terms */
272 AntennaPatternMatrix *Mmunu1 = &multiAM1->Mmunu;
273 AntennaPatternMatrix *Mmunu2 = &multiAM2->Mmunu;
274
275 REAL8 maxxerr_Ad = 0, maxxerr_Bd = 0, maxxerr_Cd = 0, maxxerr_Dd = 0;
276 REAL8 err;
277
278 err = fabs( Mmunu1->Ad - Mmunu2->Ad );
279 if ( err > maxxerr_Ad ) {
280 maxxerr_Ad = err;
281 }
282 err = fabs( Mmunu1->Bd - Mmunu2->Bd );
283 if ( err > maxxerr_Bd ) {
284 maxxerr_Bd = err;
285 }
286 err = fabs( Mmunu1->Cd - Mmunu2->Cd );
287 if ( err > maxxerr_Cd ) {
288 maxxerr_Cd = err;
289 }
290 err = fabs( Mmunu1->Dd - Mmunu2->Dd );
291 if ( err > maxxerr_Dd ) {
292 maxxerr_Dd = err;
293 }
294
295 UINT4 failed = 0;
296 /* special treatment for Sinv_Tsft: independent of AM-functions, should agree to within numerics */
297 double eps = 1e-15;
298 if ( gsl_fcmp( Mmunu1->Sinv_Tsft, Mmunu2->Sinv_Tsft, eps ) ) {
299 XLALPrintError( "%s: Sinv_Tsft differs by more than %g relative error\n", __func__, eps );
300 failed ++;
301 } else {
302 XLALPrintInfo( "%s: Sinv_Tsft 1 = %g, 2 = %g, 1-2 = %g\n", __func__, Mmunu1->Sinv_Tsft, Mmunu2->Sinv_Tsft, Mmunu1->Sinv_Tsft - Mmunu2->Sinv_Tsft );
303 }
304
305 /* ----- compare matrix elements A,B,C -------------------- */
306 /* Mmunu = {A,B,C} are sums of N terms a^2, b^2, a*b, each with small numerical errors
307 * we therefore need to relax the tolerance from a,b,c by a factor of sqrt(N):
308 */
309 REAL8 tolerance_Mmunu = tolerance * sqrt( 1.0 * numTerms );
310 XLALPrintError( "%s: tolerances tol(a,b,c)=%g, tol(Mmunu) = %g\n", __func__, tolerance, tolerance_Mmunu );
311
312 if ( maxxerr_Ad > tolerance_Mmunu ) {
313 XLALPrintError( "%s: maximal difference in Ad is %g, which exceeds the tolerance %g\n", __func__, maxxerr_Ad, tolerance_Mmunu );
314 failed ++;
315 } else {
316 XLALPrintInfo( "%s: maxxerr_Ad = %g (< %g)\n", __func__, maxxerr_Ad, tolerance_Mmunu );
317 }
318
319 if ( maxxerr_Bd > tolerance_Mmunu ) {
320 XLALPrintError( "%s: maximal difference in Bd is %g, which exceeds the tolerance %g\n", __func__, maxxerr_Bd, tolerance_Mmunu );
321 failed ++;
322 } else {
323 XLALPrintInfo( "%s: maxxerr_Bd = %g (< %g)\n", __func__, maxxerr_Bd, tolerance_Mmunu );
324 }
325
326 if ( maxxerr_Cd > tolerance_Mmunu ) {
327 XLALPrintError( "%s: maximal difference in Cd is %g, which exceeds the tolerance %g\n", __func__, maxxerr_Cd, tolerance_Mmunu );
328 failed ++;
329 } else {
330 XLALPrintInfo( "%s: maxxerr_Cd = %g (< %g)\n", __func__, maxxerr_Cd, tolerance_Mmunu );
331 }
332
333 /* matrix can be quite ill-conditioned, so the error on Dd is very hard to constrain.
334 * in principle for random-parameters this can go arbitrarly badly, so we don't use
335 * any constraints in this test to avoid spurious test-failures
336 */
337 XLALPrintError( "%s: maxxerr_Dd = %g (%g)\n", __func__, maxxerr_Dd, tolerance_Mmunu );
338
339 /* ----- compare individual a,b,c errors -------------------- */
340 if ( maxerr_ab > tolerance ) {
341 XLALPrintError( "%s: maximal difference in {a, b} coefficients is %g, which exceeds the tolerance %g\n", __func__, maxerr_ab, tolerance );
342 failed ++;
343 } else {
344 XLALPrintInfo( "%s: maxerr_ab = %g (< %g)\n", __func__, maxerr_ab, tolerance );
345 }
346
347 if ( avgerr_ab > tolerance ) {
348 XLALPrintError( "%s: average difference in {a, b} coefficients is %g, which exceeds the tolerance %g\n", __func__, avgerr_ab, tolerance );
349 failed ++;
350 } else {
351 XLALPrintInfo( "%s: avgerr_ab = %g (< %g)\n", __func__, avgerr_ab, tolerance );
352 }
353
354 return failed;
355
356} /* XLALCompareMultiAMCoeffs() */
357
358/*
359 * Multi-IFO version of old_LALGetAMCoeffs().
360 * Get all antenna-pattern coefficients for all input detector-series.
361 *
362 * NOTE: contrary to old_LALGetAMCoeffs(), this functions *allocates* the output-vector,
363 * use XLALDestroyMultiAMCoeffs() to free this.
364 */
365void
366old_LALGetMultiAMCoeffs( LALStatus *status, /*< [in/out] LAL status structure pointer */
367 MultiAMCoeffs **multiAMcoef, /*< [out] AM-coefficients for all input detector-state series */
368 const MultiDetectorStateSeries *multiDetStates, /*< [in] detector-states at timestamps t_i */
369 SkyPosition skypos /*< source sky-position [in equatorial coords!] */
370 )
371{
373 MultiAMCoeffs *ret = NULL;
374
377
378 /* check input */
384
385 numDetectors = multiDetStates->length;
386
387 if ( ( ret = LALCalloc( 1, sizeof( *ret ) ) ) == NULL ) {
389 }
390 ret->length = numDetectors;
391 if ( ( ret->data = LALCalloc( numDetectors, sizeof( *ret->data ) ) ) == NULL ) {
392 LALFree( ret );
394 }
395
396 for ( X = 0; X < numDetectors; X ++ ) {
397 AMCoeffs *amcoeX = NULL;
398 UINT4 numStepsX = multiDetStates->data[X]->length;
399
400 ret->data[X] = LALCalloc( 1, sizeof( *( ret->data[X] ) ) );
401 amcoeX = ret->data[X];
402 amcoeX->a = XLALCreateREAL4Vector( numStepsX );
403 if ( ( amcoeX->b = XLALCreateREAL4Vector( numStepsX ) ) == NULL ) {
404 LALPrintError( "\nOut of memory!\n\n" );
405 goto failed;
406 }
407
408 old_LALNewGetAMCoeffs( status->statusPtr, amcoeX, multiDetStates->data[X], skypos );
409 if ( status->statusPtr->statusCode ) {
410 LALPrintError( "\nCall to old_LALNewGetAMCoeffs() has failed ... \n\n" );
411 goto failed;
412 }
413
414 } /* for X < numDetectors */
415
416 goto success;
417
418failed:
419 /* free all memory allocated so far */
421 ABORT( status, -1, "old_LALGetMultiAMCoeffs() failed" );
422
423success:
424 ( *multiAMcoef ) = ret;
425
427 RETURN( status );
428
429} /* old_LALGetMultiAMCoeffs() */
430
431/*
432 * Compute the 'amplitude coefficients' \f$ a(t)\sin\zeta \f$ ,
433 * \f$ b(t)\sin\zeta \f$ as defined in \cite JKS98 for a series of
434 * timestamps.
435 *
436 * The input consists of the DetectorState-timeseries, which contains
437 * the detector-info and the LMST's corresponding to the different times.
438 *
439 * In order to allow re-using the output-structure AMCoeffs for subsequent
440 * calls, we require the REAL4Vectors a and b to be allocated already and
441 * to have the same length as the DetectoStates-timeseries.
442 *
443 * \note This is an alternative implementation to both LALComputeAM()
444 * and old_LALGetAMCoeffs(), which uses the geometrical definition of
445 * \f$ a\sin\zeta \f$ and \f$ b\sin\zeta \f$ as detector response
446 * coefficients in a preferred polarization basis. (It is thereby
447 * more general than the JKS expressions and could be used e.g., with
448 * the response tensor of a bar detector with no further modification
449 * needed.)
450 */
451void
452old_LALNewGetAMCoeffs( LALStatus *status, /*< [in/out] LAL status structure pointer */
453 AMCoeffs *coeffs, /*< [out] amplitude-coeffs {a(t_i), b(t_i)} */
454 const DetectorStateSeries *DetectorStates,/*< timeseries of detector states */
455 SkyPosition skypos /*< {alpha,delta} of the source */
456 )
457{
459 REAL4 sin1delta, cos1delta;
460 REAL4 sin1alpha, cos1alpha;
461
462 REAL4 xi1, xi2;
463 REAL4 eta1, eta2, eta3;
464 REAL4 norm;
465 UINT4 i, numSteps;
466
468
469 /*---------- check input ---------- */
471
472 numSteps = DetectorStates->length;
473
474 /* require the coeffients-vectors to be allocated and consistent with timestamps */
477 ASSERT( ( coeffs->a->length == numSteps ) && ( coeffs->b->length == numSteps ), status,
479
480 /* require sky-pos to be in equatorial coordinates */
482 SKYCOORDINATESH_ESYS, SKYCOORDINATESH_MSGESYS );
483
484 /*---------- We write components of xi and eta vectors in SSB-fixed coords */
485 alpha = skypos.longitude;
486 delta = skypos.latitude;
487
488 if ( XLALSinCosLUT( &sin1delta, &cos1delta, delta ) != XLAL_SUCCESS ) {
489 ABORT( status->statusPtr, LAL_EXLAL, "XLALSinCosLUT (&sin1delta, &cos1delta, delta ) failed" );
490 }
491
492 if ( XLALSinCosLUT( &sin1alpha, &cos1alpha, alpha ) != XLAL_SUCCESS ) {
493 ABORT( status->statusPtr, LAL_EXLAL, "XLALSinCosLUT (&sin1alpha, &cos1alpha, alpha ) failed" );
494 }
495
496 // see Eq.(17) in CFSv2 notes (version v3):
497 // https://dcc.ligo.org/cgi-bin/private/DocDB/ShowDocument?docid=1665&version=3
498 xi1 = sin1alpha;
499 xi2 = -cos1alpha;
500 eta1 = -sin1delta * cos1alpha;
501 eta2 = -sin1delta * sin1alpha;
502 eta3 = cos1delta;
503
504 /*---------- Compute the a(t_i) and b(t_i) ---------- */
505 coeffs->A = 0;
506 coeffs->B = 0;
507 coeffs->C = 0;
508 coeffs->D = 0;
509 for ( i = 0; i < numSteps; i++ ) {
510 REAL4 ai, bi;
511
512 SymmTensor3 *d = &( DetectorStates->data[i].detT );
513
514 ai = d->d11 * ( xi1 * xi1 - eta1 * eta1 )
515 + 2 * d->d12 * ( xi1 * xi2 - eta1 * eta2 )
516 - 2 * d->d13 * eta1 * eta3
517 + d->d22 * ( xi2 * xi2 - eta2 * eta2 )
518 - 2 * d->d23 * eta2 * eta3
519 - d->d33 * eta3 * eta3;
520
521 bi = d->d11 * 2 * xi1 * eta1
522 + 2 * d->d12 * ( xi1 * eta2 + xi2 * eta1 )
523 + 2 * d->d13 * xi1 * eta3
524 + d->d22 * 2 * xi2 * eta2
525 + 2 * d->d23 * xi2 * eta3;
526
527 /*
528 printf("xi = (%f,%f)\n",xi1,xi2);
529 printf("eta = (%f,%f,%f)\n",eta1,eta2,eta3);
530 printf("d = (%f %f %f\n",d->d11,d->d12,d->d13);
531 printf(" %f %f %f\n",d->d12,d->d22,d->d23);
532 printf(" %f %f %f)\n",d->d13,d->d23,d->d33);
533 */
534
535 coeffs->a->data[i] = ai;
536 coeffs->b->data[i] = bi;
537
538 /* sum A, B, C on the fly */
539 coeffs->A += ai * ai;
540 coeffs->B += bi * bi;
541 coeffs->C += ai * bi;
542
543 } /* for i < numSteps */
544
545 /* finish calculation of A,B,C, D */
546 norm = 2.0f / numSteps;
547 coeffs->A *= norm;
548 coeffs->B *= norm;
549 coeffs->C *= norm;
550
551 coeffs->D = coeffs->A * coeffs->B - coeffs->C * coeffs->C;
552
553 RETURN( status );
554
555} /* old_LALNewGetAMCoeffs() */
const REAL8 eps
#define __func__
log an I/O error, i.e.
void LALCheckMemoryLeaks(void)
#define LALCalloc(m, n)
#define LALFree(p)
static double double delta
#define ABORT(statusptr, code, mesg)
#define LAL_EXLAL
#define ATTATCHSTATUSPTR(statusptr)
#define ASSERT(assertion, statusptr, code, mesg)
#define DETATCHSTATUSPTR(statusptr)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
int LALgetopt(int argc, char *const *argv, const char *optstring)
char * LALoptarg
int main(int argc, char *argv[])
#define LALCOMPUTEAMH_EINPUT
#define LALCOMPUTEAMH_MSGENONULL
#define LALCOMPUTEAMH_EMEM
#define LALCOMPUTEAMH_MSGEMEM
#define LALCOMPUTEAMH_ENULL
#define LALCOMPUTEAMH_ENONULL
#define LALCOMPUTEAMH_MSGENULL
static void old_LALNewGetAMCoeffs(LALStatus *, AMCoeffs *coeffs, const DetectorStateSeries *DetectorStates, SkyPosition skypos)
#define LALCOMPUTEAMH_MSGEINPUT
static void old_LALGetMultiAMCoeffs(LALStatus *, MultiAMCoeffs **multiAMcoef, const MultiDetectorStateSeries *multiDetStates, SkyPosition pos)
static int XLALCompareMultiAMCoeffs(MultiAMCoeffs *multiAM1, MultiAMCoeffs *multiAM2, REAL8 tolerance)
double e
const double xi2
void XLALDestroyMultiDetectorStateSeries(MultiDetectorStateSeries *mdetStates)
Helper function to get rid of a multi-IFO DetectorStateSeries Note, this is "NULL-robust" in the sens...
MultiDetectorStateSeries * XLALGetMultiDetectorStates(const MultiLIGOTimeGPSVector *multiTS, const MultiLALDetector *multiIFO, const EphemerisData *edat, REAL8 tOffset)
Get the detector-time series for the given MultiLIGOTimeGPSVector.
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
void XLALDestroyEphemerisData(EphemerisData *edat)
Destructor for EphemerisData struct, NULL robust.
void XLALDestroyMultiAMCoeffs(MultiAMCoeffs *multiAMcoef)
Destroy a MultiAMCoeffs structure.
Definition: LALComputeAM.c:469
int XLALWeightMultiAMCoeffs(MultiAMCoeffs *multiAMcoef, const MultiNoiseWeights *multiWeights)
Replace AM-coeffs by weighted AM-coeffs, i.e.
Definition: LALComputeAM.c:188
MultiAMCoeffs * XLALComputeMultiAMCoeffs(const MultiDetectorStateSeries *multiDetStates, const MultiNoiseWeights *multiWeights, SkyPosition skypos)
Multi-IFO version of XLALComputeAMCoeffs().
Definition: LALComputeAM.c:379
#define LAL_PI_2
#define LAL_TWOPI
double REAL8
#define XLAL_INIT_DECL(var,...)
char CHAR
uint32_t UINT4
float REAL4
int LALPrintError(const char *fmt,...)
void * XLALCalloc(size_t m, size_t n)
const LALDetector * XLALGetSiteInfo(const CHAR *name)
Find the site geometry-information 'LALDetector' for given a detector name (or prefix).
Definition: SFTnaming.c:218
LIGOTimeGPSVector * XLALCreateTimestampVector(UINT4 len)
Allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:47
void XLALDestroyMultiTimestamps(MultiLIGOTimeGPSVector *multiTS)
Destroy a MultiLIGOTimeGPSVector timestamps vector.
int XLALSinCosLUT(REAL4 *sinx, REAL4 *cosx, REAL8 x)
Calculate sin(x) and cos(x) to roughly 1e-7 precision using a lookup-table and Taylor-expansion.
Definition: SinCosLUT.c:83
#define SKYCOORDINATESH_ESYS
COORDINATESYSTEM_EQUATORIAL
REAL4Vector * XLALCreateREAL4Vector(UINT4 length)
#define xlalErrno
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_ERROR(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EINVAL
XLAL_EFAILED
list weights
double alpha
size_t numDetectors
This structure contains the per-SFT (weighted) antenna-pattern functions , with the SFT-index,...
Definition: LALComputeAM.h:63
REAL4 B
summed antenna-pattern matrix coefficient:
Definition: LALComputeAM.h:67
REAL4Vector * b
(weighted) per-SFT antenna-pattern function
Definition: LALComputeAM.h:65
REAL4 A
summed antenna-pattern matrix coefficient:
Definition: LALComputeAM.h:66
REAL4 C
summed antenna-pattern matrix coefficient:
Definition: LALComputeAM.h:68
REAL4Vector * a
(weighted) per-SFT antenna-pattern function
Definition: LALComputeAM.h:64
REAL4 D
determinant
Definition: LALComputeAM.h:69
Struct holding the "antenna-pattern" matrix , in terms of the multi-detector scalar product.
Definition: LALComputeAM.h:127
REAL8 Sinv_Tsft
normalization-factor (using single-sided PSD!)
Definition: LALComputeAM.h:133
REAL4 Dd
determinant factor , such that
Definition: LALComputeAM.h:132
SymmTensor3 detT
Detector-tensor components in SSB-fixed, Cartesian coordinates.
Timeseries of DetectorState's, representing the detector-info at different timestamps.
DetectorState * data
array of DetectorState entries
UINT4 length
total number of entries
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
Multi-IFO container for antenna-pattern coefficients and atenna-pattern matrix .
Definition: LALComputeAM.h:137
UINT4 length
number of IFOs
Definition: LALComputeAM.h:141
AMCoeffs ** data
noise-weighted AM-coeffs , and
Definition: LALComputeAM.h:142
AntennaPatternMatrix Mmunu
antenna-pattern matrix
Definition: LALComputeAM.h:143
Multi-IFO time-series of DetectorStates.
UINT4 length
number of detectors
DetectorStateSeries ** data
vector of pointers to DetectorStateSeries
array of detectors definitions 'LALDetector'
UINT4 length
number of detectors
LALDetector sites[PULSAR_MAX_DETECTORS]
array of site information
A collection of (multi-IFO) LIGOTimeGPSVector time-stamps vectors.
Definition: SFTfileIO.h:198
One noise-weight (number) per SFT (therefore indexed over IFOs and SFTs.
Definition: PSDutils.h:71
REAL4 * data
REAL8 longitude
REAL8 latitude
CoordinateSystem system
A symmetric 3x3 tensor (such as detector-tensors), storing only the upper triangle.
enum @4 site