Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-5e288d3
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LISAspecifics.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2006, 2007 Reinhard Prix
3 * Copyright (C) 2006, 2007, 2008 John T. Whelan
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/**
22 * \author R. Prix, J. T. Whelan
23 * \file
24 * \brief
25 * LISA-specific implementations for Fstat/continuous-wave searches on LISA TDI observables.
26 *
27 */
28
29/*---------- INCLUDES ----------*/
30#include <math.h>
31#include <string.h>
32
33#include <lal/LISAspecifics.h>
34#include <lal/LALError.h>
35#include <lal/SinCosLUT.h>
36
37/*---------- local DEFINES ----------*/
38#define TRUE (1==1)
39#define FALSE (1==0)
40
41/*----- Macros ----- */
42/** Simple Euklidean scalar product for two 3-dim vectors in cartesian coords */
43#define SCALAR(u,v) ((u)[0]*(v)[0] + (u)[1]*(v)[1] + (u)[2]*(v)[2])
44#define SQ(x) ( (x) * (x) )
45#define GPS2REAL8(gps) (1.0 * (gps).gpsSeconds + 1.e-9 * (gps).gpsNanoSeconds )
46#define LAL_SQRT1_3 0.5773502691896257645091487805019575 /**< 1/sqrt(3) */
47/*----- SWITCHES -----*/
48
49/*---------- internal types ----------*/
50
51/*---------- Global variables ----------*/
52
53/*---------- internal prototypes ----------*/
54int XLALgetLISAtwoArmRAAIFO( CmplxDetectorTensor *detT, const DetectorArm *armA, const DetectorArm *armB, const FreqSkypos_t *freq_skypos );
55
56static REAL4 safe_sinc( REAL4 x );
57
58/*==================== FUNCTION DEFINITIONS ====================*/
59
60
61/**
62 * Return true if 'det' is a LISA \em LALDetector struct
63 */
64BOOLEAN
65XLALisLISAdetector( const LALDetector *det )
66{
67 return ( det != NULL ) && ( strncmp( det->frDetector.name, "LISA TDI ", 9 ) == 0 );
68}
69
70/**
71 * Set up the \em LALDetector structs representing LISA X, Y, Z TDI observables.
72 * Detectors will be registered with prefixes beginning with 'prefixLetter'.
73 */
74int
75XLALregisterLISAdetectors( const CHAR prefixLetter )
76{
77
78 for ( CHAR channelNum = '1'; channelNum <= '9'; ++channelNum ) {
79 LALDetector XLAL_INIT_DECL( Detector1 );
80
81 switch ( channelNum ) {
82 case '1':
83 strcpy( Detector1.frDetector.name, "LISA TDI X" );
84 break;
85 case '2':
86 strcpy( Detector1.frDetector.name, "LISA TDI Y" );
87 break;
88 case '3':
89 strcpy( Detector1.frDetector.name, "LISA TDI Z" );
90 break;
91 case '4':
92 strcpy( Detector1.frDetector.name, "LISA TDI Y-Z" );
93 break;
94 case '5':
95 strcpy( Detector1.frDetector.name, "LISA TDI Z-X" );
96 break;
97 case '6':
98 strcpy( Detector1.frDetector.name, "LISA TDI X-Y" );
99 break;
100 case '7':
101 strcpy( Detector1.frDetector.name, "LISA TDI A" );
102 break;
103 case '8':
104 strcpy( Detector1.frDetector.name, "LISA TDI E" );
105 break;
106 case '9':
107 strcpy( Detector1.frDetector.name, "LISA TDI T" );
108 break;
109 default:
110 XLAL_ERROR( XLAL_EINVAL, "Illegal LISA TDI index '%c': must be one of {'1', '2', '3', '4', '5', '6', '7', '8', '9'}.\n\n", channelNum );
111 break;
112 } /* switch (channelNum) */
113
114 Detector1.frDetector.prefix[0] = prefixLetter;
115 Detector1.frDetector.prefix[1] = channelNum;
116
117 /* fill frDetector with dummy numbers: meaningless for LISA */
118 Detector1.frDetector.vertexLongitudeRadians = 0;
119 Detector1.frDetector.vertexLatitudeRadians = 0;
120 Detector1.frDetector.vertexElevation = 0;
121 Detector1.frDetector.xArmAltitudeRadians = 0;
122 Detector1.frDetector.xArmAzimuthRadians = 0;
123 Detector1.frDetector.yArmAltitudeRadians = 0;
124 Detector1.frDetector.yArmAzimuthRadians = 0;
125
126 Detector1.type = LALDETECTORTYPE_ABSENT;
127
128 /* however: need to be careful to put some non-zero numbers for location
129 * otherwise XLALBarycenter() will spit out NaNs ...
130 */
131 Detector1.location[0] = 1;
132 Detector1.location[1] = 1;
133 Detector1.location[2] = 1;
134
135 XLAL_CHECK( XLALRegisterSpecialCWDetector( &Detector1 ) == XLAL_SUCCESS, XLAL_EFUNC );
136
137 }
138
139 return XLAL_SUCCESS;
140
141} /* XLALregisterLISAdetectors() */
142
143
144/**
145 * Precompute the arm-geometry for LISA, which is later used
146 * for assembling the RAA detector-tensor (which depends on
147 * frequency and skyposition
148 */
149int
150XLALprecomputeLISAarms( DetectorState *detState )
151{
152 REAL4 x[3], y[3], z[3]; /* position x,y,z of spacecraft n */
153
154 enum { SC1 = 0, SC2, SC3 }; /* translate TDI-conventions to C-indexing */
155 enum { I1 = 0, I2, I3 }; /* convenience indices */
156
157 UINT4 send_l[3] = { SC3, SC1, SC2 }; /* canonical sender-spacecraft 's' for arm 'l' */
158 UINT4 rec_l [3] = { SC2, SC3, SC1 }; /* canonical receiver-spacecraft 'r' for arm 'l' */
159
160 UINT4 arm;
161
162 if ( !detState ) {
163 return -1;
164 }
165
166 { /* determine space-craft positions [MLDC Challenge 1] */
167 REAL4 e = 0.00965; /* LISA eccentricity */
168 REAL4 ae = LAL_AU_SI * e / LAL_C_SI; /* in units of time */
169 REAL4 kappa = 0, lambda = 0; /* MLDC default */
170 REAL4 sin_alpha, cos_alpha;
171
172 REAL8 Om = LAL_TWOPI / LAL_YRSID_SI;
173 REAL8 alpha_t;
174 UINT4 i;
175 REAL8 tGPS = GPS2REAL8( detState->tGPS ) - LISA_TIME_ORIGIN;
176
177 alpha_t = Om * tGPS + kappa;
178 XLAL_CHECK( XLALSinCosLUT( &sin_alpha, &cos_alpha, alpha_t ) == XLAL_SUCCESS, XLAL_EFUNC );
179
180 for ( i = 0; i < 3; i ++ ) {
181 REAL4 beta = 2.0 * i * LAL_PI / 3.0 + lambda; /* relative orbital phase in constellation */
182 REAL4 sin_beta, cos_beta;
183 XLAL_CHECK( XLALSinCosLUT( &sin_beta, &cos_beta, beta ) == XLAL_SUCCESS, XLAL_EFUNC );
184
185 x[i] = detState->earthState.posNow[0] + ae * ( sin_alpha * cos_alpha * sin_beta - ( 1.0f + SQ( sin_alpha ) ) * cos_beta );
186 y[i] = detState->earthState.posNow[1] + ae * ( sin_alpha * cos_alpha * cos_beta - ( 1.0f + SQ( cos_alpha ) ) * sin_beta );
187 z[i] = detState->earthState.posNow[2] - sqrt( 3.0f ) * ae * ( cos_alpha * cos_beta + sin_alpha * sin_beta );
188 } /* for i=[0:2] */
189
190 } /* determine spacecraft positions */
191
192 for ( arm = 0; arm < 3; arm ++ ) {
193 UINT4 send, rec;
194 REAL4 n[3], L_c, invL;
195
196 /* get canonical (ie following TDI conventions) senders and receivers this arm */
197 send = send_l[arm];
198 rec = rec_l [arm];
199
200 /* get un-normalized arm-vectors first */
201 n[0] = x[rec] - x[send];
202 n[1] = y[rec] - y[send];
203 n[2] = z[rec] - z[send];
204
205 /* get armlength in seconds, and normalize */
206 L_c = sqrt( SCALAR( n, n ) );
207 detState->detArms[arm].armlength_c = L_c;
208 invL = 1.0f / L_c;
209
210 n[0] *= invL;
211 n[1] *= invL;
212 n[2] *= invL;
213
214 /* pre-compute the "basis tensor" n x n for this arm */
215 XLALTensorSquareVector3( &( detState->detArms[arm].basisT ), n );
216
217 /* store arm unit-vector */
218 detState->detArms[arm].n[0] = n[0];
219 detState->detArms[arm].n[1] = n[1];
220 detState->detArms[arm].n[2] = n[2];
221
222 } /* for arm = 0, 1, 2 */
223
224 return 0;
225
226} /* XLALprecomputeLISAarms() */
227
228
229/* Construct the long-wavelength-limit (LWL) detector tensor for LISA, given the prefix
230 * and the precomputed detector arms.
231 *
232 * RETURN 0 = OK, -1 = ERROR
233 */
234int
235XLALgetLISADetectorTensorLWL( SymmTensor3 *detT, /**< [out]: LISA LWL detector-tensor */
236 const Detector3Arms detArms, /**< [in] precomputed detector-arms */
237 CHAR channelNum ) /**< channel-number (as a char) '1', '2', '3' .. */
238{
239 LISAarmT armA = 0, armB = 0;
240 CHAR chan1 = 0;
241 CHAR chan2 = 0;
242 SymmTensor3 detT1, detT2;
243
244 if ( !detT ) {
245 return -1;
246 }
247
248 /* we distinuish (currently) 3 different TDI observables: X (channel=1), Y (channel=2), Z (channel=3) */
249 switch ( channelNum ) {
250 case '1': /* TDI observable 'X' */
251 armA = LISA_ARM2;
252 armB = LISA_ARM3;
253 break;
254 case '2': /* TDI observable 'Y' */
255 armA = LISA_ARM3;
256 armB = LISA_ARM1;
257 break;
258 case '3': /* TDI observable 'Z' */
259 armA = LISA_ARM1;
260 armB = LISA_ARM2;
261 break;
262 case '4':
263 chan1 = '2';
264 chan2 = '3';
265 break;
266 case '5':
267 chan1 = '3';
268 chan2 = '1';
269 break;
270 case '6':
271 chan1 = '1';
272 chan2 = '2';
273 break;
274 case '7':
275 case '8':
276 case '9':
277 chan2 = 'X';
278 break;
279 default: /* unknown */
280 XLALPrintError( "\nInvalid channel-number '%c' for LISA \n\n", channelNum );
281 xlalErrno = XLAL_EINVAL;
282 return -1;
283 break;
284 } /* switch channel[1] */
285
286 if ( chan1 == 0 ) {
287 if ( chan2 == 0 ) {
288 XLALSubtractSymmTensor3s( detT, &( detArms[armA].basisT ), &( detArms[armB].basisT ) );
289 XLALScaleSymmTensor3( detT, detT, 0.5f ); /* (1/2)*(nA x nA - nB x nB ) */
290 } else {
291 REAL8 weight1, weight2, weight3;
292 SymmTensor3 detT3;
293
294 if ( XLALgetLISADetectorTensorLWL( &detT1, detArms, '1' ) != 0 ) {
295 XLALPrintError( "\nXLALgetLISADetectorTensorLWL() failed !\n\n" );
296 xlalErrno = XLAL_EINVAL;
297 return -1;
298 }
299 if ( XLALgetLISADetectorTensorLWL( &detT2, detArms, '2' ) != 0 ) {
300 XLALPrintError( "\nXLALgetLISADetectorTensorLWL() failed !\n\n" );
301 xlalErrno = XLAL_EINVAL;
302 return -1;
303 }
304 if ( XLALgetLISADetectorTensorLWL( &detT3, detArms, '3' ) != 0 ) {
305 XLALPrintError( "\nXLALgetLISADetectorTensorLWL() failed !\n\n" );
306 xlalErrno = XLAL_EINVAL;
307 return -1;
308 }
309 switch ( channelNum ) {
310 case '7':
311 weight1 = ( 2.0 / 3.0 );
312 weight2 = ( -1.0 / 3.0 );
313 weight3 = ( -1.0 / 3.0 );
314 break;
315 case '8':
316 weight1 = 0.0;
317 weight2 = ( -1.0 * LAL_SQRT1_3 );
318 weight3 = LAL_SQRT1_3;
319 break;
320 case '9':
321 /* Note that in the LWL this will give detT = 0 */
322 weight1 = ( LAL_SQRT2 / 3.0 );
323 weight2 = ( LAL_SQRT2 / 3.0 );
324 weight3 = ( LAL_SQRT2 / 3.0 );
325 break;
326 default:
327 /* should never get to default */
328 XLALPrintError( "\nInvalid channel-number '%c' for LISA \n\n", channelNum );
329 xlalErrno = XLAL_EINVAL;
330 return -1;
331 break;
332 }
333
334 /* Note that the following in-place operations *are* safe */
335
336 if ( XLALScaleSymmTensor3( &detT1, &detT1, weight1 ) != 0 ) {
337 XLALPrintError( "\nXLALScaleSymmTensor3() failed !\n\n" );
338 xlalErrno = XLAL_EINVAL;
339 return -1;
340 }
341 if ( XLALScaleSymmTensor3( &detT2, &detT2, weight2 ) != 0 ) {
342 XLALPrintError( "\nXLALScaleSymmTensor3() failed !\n\n" );
343 xlalErrno = XLAL_EINVAL;
344 return -1;
345 }
346 if ( XLALScaleSymmTensor3( &detT3, &detT3, weight3 ) != 0 ) {
347 XLALPrintError( "\nXLALScaleSymmTensor3() failed !\n\n" );
348 xlalErrno = XLAL_EINVAL;
349 return -1;
350 }
351 if ( XLALAddSymmTensor3s( &detT2, &detT2, &detT3 ) != 0 ) {
352 XLALPrintError( "\nXLALSumSymmTensor3s() failed !\n\n" );
353 xlalErrno = XLAL_EINVAL;
354 return -1;
355 }
356 if ( XLALAddSymmTensor3s( detT, &detT1, &detT2 ) != 0 ) {
357 XLALPrintError( "\nXLALSumSymmTensor3s() failed !\n\n" );
358 xlalErrno = XLAL_EINVAL;
359 return -1;
360 }
361 } /* if (chan2 != 0) */
362 } else {
363 if ( XLALgetLISADetectorTensorLWL( &detT1, detArms, chan1 ) != 0 ) {
364 XLALPrintError( "\nXLALgetLISADetectorTensorLWL() failed !\n\n" );
365 xlalErrno = XLAL_EINVAL;
366 return -1;
367 }
368 if ( XLALgetLISADetectorTensorLWL( &detT2, detArms, chan2 ) != 0 ) {
369 XLALPrintError( "\nXLALgetLISADetectorTensorLWL() failed !\n\n" );
370 xlalErrno = XLAL_EINVAL;
371 return -1;
372 }
373
374 if ( XLALSubtractSymmTensor3s( detT, &detT1, &detT2 ) != 0 ) {
375 XLALPrintError( "\nXLALSubtractSymmTensor3s() failed !\n\n" );
376 xlalErrno = XLAL_EINVAL;
377 return -1;
378 } /* d_X - d_Y etc */
379 } /* multi-channel "detector" such as 'X-Y' etc */
380
381 return 0;
382
383} /* XLALgetLISADetectorTensorLWL() */
384
385
386/* Construct the rigid-adiabatic-approximation (RAA) detector tensor for LISA, given the prefix,
387 * the precomputed arms, and the Frequency and skyposition.
388 *
389 * RETURN 0 = OK, -1 = ERROR
390 */
391int
392XLALgetLISADetectorTensorRAA( CmplxDetectorTensor *detT, /**< [out]: LISA LWL detector-tensor */
393 const Detector3Arms detArms, /**< [in] precomputed detector-arms */
394 CHAR channelNum, /**< [in] channel-number (as a char) '1', '2', '3' .. */
395 const FreqSkypos_t *freq_skypos /**< [in] precompute frequency and skypos info */
396 )
397{
398 LISAarmT armA = 0, armB = 0;
399 CHAR chan1 = 0;
400 CHAR chan2 = 0;
401 CmplxDetectorTensor detT1, detT2;
402
403 if ( !detT ) {
404 return -1;
405 }
406
407 /* we distinuish (currently) 3 different TDI observables: X (channel=1), Y (channel=2), Z (channel=3) */
408 switch ( channelNum ) {
409 case '1': /* TDI observable 'X' */
410 armA = LISA_ARM2;
411 armB = LISA_ARM3;
412 break;
413 case '2': /* TDI observable 'Y' */
414 armA = LISA_ARM3;
415 armB = LISA_ARM1;
416 break;
417 case '3': /* TDI observable 'Z' */
418 armA = LISA_ARM1;
419 armB = LISA_ARM2;
420 break;
421 case '4':
422 chan1 = '2';
423 chan2 = '3';
424 break;
425 case '5':
426 chan1 = '3';
427 chan2 = '1';
428 break;
429 case '6':
430 chan1 = '1';
431 chan2 = '2';
432 break;
433 case '7':
434 case '8':
435 case '9':
436 chan2 = 'X';
437 break;
438 default: /* unknown */
439 XLALPrintError( "\nInvalid channel-number '%c' for LISA \n\n", channelNum );
440 xlalErrno = XLAL_EINVAL;
441 return -1;
442 break;
443 } /* switch channel[1] */
444
445 if ( chan1 == 0 ) {
446 if ( chan2 == 0 ) {
447 if ( XLALgetLISAtwoArmRAAIFO( detT, &( detArms[armA] ),
448 &( detArms[armB] ), freq_skypos )
449 != 0 ) {
450 XLALPrintError( "\nXLALgetLISAtwoArmRAAIFO() failed !\n\n" );
451 xlalErrno = XLAL_EINVAL;
452 return -1;
453 }
454 } else {
455 REAL8 weight1, weight2, weight3;
456 CmplxDetectorTensor detT3;
457
458 if ( XLALgetLISADetectorTensorRAA( &detT1, detArms, '1',
459 freq_skypos ) != 0 ) {
460 XLALPrintError( "\nXLALgetLISADetectorTensorRAA() failed !\n\n" );
461 xlalErrno = XLAL_EINVAL;
462 return -1;
463 }
464 if ( XLALgetLISADetectorTensorRAA( &detT2, detArms, '2',
465 freq_skypos ) != 0 ) {
466 XLALPrintError( "\nXLALgetLISADetectorTensorRAA() failed !\n\n" );
467 xlalErrno = XLAL_EINVAL;
468 return -1;
469 }
470 if ( XLALgetLISADetectorTensorRAA( &detT3, detArms, '3',
471 freq_skypos ) != 0 ) {
472 XLALPrintError( "\nXLALgetLISADetectorTensorRAA() failed !\n\n" );
473 xlalErrno = XLAL_EINVAL;
474 return -1;
475 }
476 switch ( channelNum ) {
477 case '7':
478 weight1 = ( 2.0 / 3.0 );
479 weight2 = ( -1.0 / 3.0 );
480 weight3 = ( -1.0 / 3.0 );
481 break;
482 case '8':
483 weight1 = 0.0;
484 weight2 = ( -1.0 * LAL_SQRT1_3 );
485 weight3 = LAL_SQRT1_3;
486 break;
487 case '9':
488 /* Note that in the RAA this will give detT != 0 */
489 weight1 = ( LAL_SQRT2 / 3.0 );
490 weight2 = ( LAL_SQRT2 / 3.0 );
491 weight3 = ( LAL_SQRT2 / 3.0 );
492 break;
493 default:
494 /* should never get to default */
495 XLALPrintError( "\nInvalid channel-number '%c' for LISA \n\n", channelNum );
496 xlalErrno = XLAL_EINVAL;
497 return -1;
498 break;
499 }
500
501 /* Note that the following in-place operations *are* safe */
502
503 if ( XLALScaleSymmTensor3( &detT1.re, &detT1.re, weight1 ) != 0 ) {
504 XLALPrintError( "\nXLALScaleSymmTensor3() failed !\n\n" );
505 xlalErrno = XLAL_EINVAL;
506 return -1;
507 }
508 if ( XLALScaleSymmTensor3( &detT1.im, &detT1.im, weight1 ) != 0 ) {
509 XLALPrintError( "\nXLALScaleSymmTensor3() failed !\n\n" );
510 xlalErrno = XLAL_EINVAL;
511 return -1;
512 }
513 if ( XLALScaleSymmTensor3( &detT2.re, &detT2.re, weight2 ) != 0 ) {
514 XLALPrintError( "\nXLALScaleSymmTensor3() failed !\n\n" );
515 xlalErrno = XLAL_EINVAL;
516 return -1;
517 }
518 if ( XLALScaleSymmTensor3( &detT2.im, &detT2.im, weight2 ) != 0 ) {
519 XLALPrintError( "\nXLALScaleSymmTensor3() failed !\n\n" );
520 xlalErrno = XLAL_EINVAL;
521 return -1;
522 }
523 if ( XLALScaleSymmTensor3( &detT3.re, &detT3.re, weight3 ) != 0 ) {
524 XLALPrintError( "\nXLALScaleSymmTensor3() failed !\n\n" );
525 xlalErrno = XLAL_EINVAL;
526 return -1;
527 }
528 if ( XLALScaleSymmTensor3( &detT3.im, &detT3.im, weight3 ) != 0 ) {
529 XLALPrintError( "\nXLALScaleSymmTensor3() failed !\n\n" );
530 xlalErrno = XLAL_EINVAL;
531 return -1;
532 }
533 if ( XLALAddSymmTensor3s( &detT2.re, &detT2.re, &detT3.re ) != 0 ) {
534 XLALPrintError( "\nXLALSumSymmTensor3s() failed !\n\n" );
535 xlalErrno = XLAL_EINVAL;
536 return -1;
537 }
538 if ( XLALAddSymmTensor3s( &detT2.im, &detT2.im, &detT3.im ) != 0 ) {
539 XLALPrintError( "\nXLALSumSymmTensor3s() failed !\n\n" );
540 xlalErrno = XLAL_EINVAL;
541 return -1;
542 }
543 if ( XLALAddSymmTensor3s( &detT->re, &detT1.re, &detT2.re ) != 0 ) {
544 XLALPrintError( "\nXLALSumSymmTensor3s() failed !\n\n" );
545 xlalErrno = XLAL_EINVAL;
546 return -1;
547 }
548 if ( XLALAddSymmTensor3s( &detT->im, &detT1.im, &detT2.im ) != 0 ) {
549 XLALPrintError( "\nXLALSumSymmTensor3s() failed !\n\n" );
550 xlalErrno = XLAL_EINVAL;
551 return -1;
552 }
553 } /* if (chan2 != 0) */
554 } else {
555 if ( XLALgetLISADetectorTensorRAA( &detT1, detArms, chan1, freq_skypos ) != 0 ) {
556 XLALPrintError( "\nXLALgetLISADetectorTensorRAA() failed !\n\n" );
557 xlalErrno = XLAL_EINVAL;
558 return -1;
559 }
560 if ( XLALgetLISADetectorTensorRAA( &detT2, detArms, chan2, freq_skypos ) != 0 ) {
561 XLALPrintError( "\nXLALgetLISADetectorTensorRAA() failed !\n\n" );
562 xlalErrno = XLAL_EINVAL;
563 return -1;
564 }
565 XLALSubtractSymmTensor3s( &detT->re, &detT1.re, &detT2.re );
566 XLALSubtractSymmTensor3s( &detT->im, &detT1.im, &detT2.im );
567 }
568
569 return 0;
570
571} /* XLALgetCmplxLISADetectorTensor() */
572
573
574/**
575 * return a rigid-adiabatic-approximation (RAA) two-arm IFO detector tensor for LISA, given 'armA' and 'armB'
576 * This implements LISA RAA-tensor using spacecraft-orbits described by Eq.(2.1) in LISA-MLCD
577 * 'challenge1.pdf' document, see http://astrogravs.nasa.gov/docs/mldc/round1.html
578 *
579 * RETURN 0 = OK, -1 = ERROR
580 *
581 * \note: armA=0 means "arm 1" in TDI notation, therefore the enums LISA_ARM1, LISA_ARM2, LISA_ARM3
582 * should be used to avoid confusion.
583 */
584int
585XLALgetLISAtwoArmRAAIFO( CmplxDetectorTensor *detT, /**< [out]: two-arm IFO detector-tensor */
586 const DetectorArm *detArmA, /**< [in] precomputed detector-arm 'A' */
587 const DetectorArm *detArmB, /**< [in] precomputed detector-arm 'B' */
588 const FreqSkypos_t *freq_skypos /**< [in] precomputed frequency and skypos info */
589 )
590{
591 REAL4 L_c = 0.5f * ( detArmA->armlength_c + detArmB->armlength_c );
592 REAL4 pifL_c = LAL_PI * freq_skypos->Freq * L_c;
593 REAL4 pifL_3c = pifL_c / 3.0f;
594 REAL4 kdotnA, kdotnB;
595 REAL4 eta, sinpha, cospha;
596 REAL4 sinc_eta;
597 COMPLEX8 coeffAA, coeffBB;
598
599 if ( !detT || !detArmA || !detArmB || !freq_skypos ) {
600 return -1;
601 }
602
603 kdotnA = - SCALAR( freq_skypos->skyposV, ( detArmA->n ) );
604 kdotnB = - SCALAR( freq_skypos->skyposV, ( detArmB->n ) );
605
606 /* ----- calculate complex coefficient in front of basis-tensor (1/2)*(nA x nA) ----- */
607
608 /* first term */
609 XLAL_CHECK( XLALSinCosLUT( &sinpha, &cospha, pifL_3c * ( 3.0f - ( kdotnA + 2.0f * kdotnB ) ) ) == XLAL_SUCCESS, XLAL_EFUNC );
610
611 eta = pifL_c * ( 1.0f + kdotnA );
612 sinc_eta = safe_sinc( eta );
613
614 coeffAA = crectf( 0.25f * cospha * sinc_eta, 0.25f * sinpha * sinc_eta );
615
616 /* second term */
617 XLAL_CHECK( XLALSinCosLUT( &sinpha, &cospha, pifL_3c * ( - 3.0f - ( kdotnA + 2.0f * kdotnB ) ) ) == XLAL_SUCCESS, XLAL_EFUNC );
618
619 eta = pifL_c * ( 1.0f - kdotnA );
620 sinc_eta = safe_sinc( eta );
621
622 coeffAA += crectf( 0.25f * cospha * sinc_eta, 0.25f * sinpha * sinc_eta );
623
624 /* ----- calculate coefficient in front of (1/2)*(nB x nB) */
625
626 /* first term */
627 XLAL_CHECK( XLALSinCosLUT( &sinpha, &cospha, pifL_3c * ( 3.0f + ( 2.0f * kdotnA + kdotnB ) ) ) == XLAL_SUCCESS, XLAL_EFUNC );
628
629 eta = pifL_c * ( 1.0f - kdotnB );
630 sinc_eta = safe_sinc( eta );
631
632 coeffBB = crectf( 0.25f * cospha * sinc_eta, 0.25f * sinpha * sinc_eta );
633
634 /* second term */
635 XLAL_CHECK( XLALSinCosLUT( &sinpha, &cospha, pifL_3c * ( - 3.0f + ( 2.0f * kdotnA + kdotnB ) ) ) == XLAL_SUCCESS, XLAL_EFUNC );
636
637 eta = pifL_c * ( 1.0f + kdotnB );
638 sinc_eta = safe_sinc( eta );
639
640 coeffBB += crectf( 0.25f * cospha * sinc_eta, 0.25f * sinpha * sinc_eta );
641
642 /* now we can express the detector tensor in the rigid adiabatic approximation (RAA):
643 * detT = coeffAA * basisA - coeffBB * basisB
644 */
645 {
646 SymmTensor3 tmpA, tmpB;
647 XLALScaleSymmTensor3( &tmpA, &detArmA->basisT, crealf( coeffAA ) );
648 XLALScaleSymmTensor3( &tmpB, &detArmB->basisT, crealf( coeffBB ) );
649 XLALSubtractSymmTensor3s( &detT->re, &tmpA, &tmpB );
650
651 XLALScaleSymmTensor3( &tmpA, &detArmA->basisT, cimagf( coeffAA ) );
652 XLALScaleSymmTensor3( &tmpB, &detArmB->basisT, cimagf( coeffBB ) );
653 XLALSubtractSymmTensor3s( &detT->im, &tmpA, &tmpB );
654 }
655
656 return 0;
657
658} /* XLALgetLISAtwoArmRAAIFO() */
659
660#define SINC_SAFETY 1e-5
661/**
662 * Unnormalized sinc(x) = sin(x) / x. Correctly handle the limit x->0
663 * where sinc(x) = 1
664 */
665static REAL4 safe_sinc( REAL4 x )
666{
667 REAL4 sinx, cosx;
668 if ( ( x > SINC_SAFETY ) || ( x < -SINC_SAFETY ) ) {
669 XLAL_CHECK( XLALSinCosLUT( &sinx, &cosx, x ) == XLAL_SUCCESS, XLAL_EFUNC );
670 return ( sinx / x );
671 } else {
672 return 1.0f;
673 }
674
675} /* sinc () */
XLALregisterLISAdetectors
int XLALregisterLISAdetectors(const CHAR prefixLetter)
Set up the LALDetector structs representing LISA X, Y, Z TDI observables.
Definition: LISAspecifics.c:75
XLALgetLISAtwoArmRAAIFO
int XLALgetLISAtwoArmRAAIFO(CmplxDetectorTensor *detT, const DetectorArm *armA, const DetectorArm *armB, const FreqSkypos_t *freq_skypos)
return a rigid-adiabatic-approximation (RAA) two-arm IFO detector tensor for LISA,...
Definition: LISAspecifics.c:585
safe_sinc
static REAL4 safe_sinc(REAL4 x)
Unnormalized sinc(x) = sin(x) / x.
Definition: LISAspecifics.c:665
XLALprecomputeLISAarms
int XLALprecomputeLISAarms(DetectorState *detState)
Precompute the arm-geometry for LISA, which is later used for assembling the RAA detector-tensor (whi...
Definition: LISAspecifics.c:150
SINC_SAFETY
#define SINC_SAFETY
Definition: LISAspecifics.c:660
XLALgetLISADetectorTensorRAA
int XLALgetLISADetectorTensorRAA(CmplxDetectorTensor *detT, const Detector3Arms detArms, CHAR channelNum, const FreqSkypos_t *freq_skypos)
Definition: LISAspecifics.c:392
LAL_SQRT1_3
#define LAL_SQRT1_3
1/sqrt(3)
Definition: LISAspecifics.c:46
XLALisLISAdetector
BOOLEAN XLALisLISAdetector(const LALDetector *det)
Return true if 'det' is a LISA LALDetector struct.
Definition: LISAspecifics.c:65
SCALAR
#define SCALAR(u, v)
Simple Euklidean scalar product for two 3-dim vectors in cartesian coords.
Definition: LISAspecifics.c:43
SQ
#define SQ(x)
Definition: LISAspecifics.c:44
XLALgetLISADetectorTensorLWL
int XLALgetLISADetectorTensorLWL(SymmTensor3 *detT, const Detector3Arms detArms, CHAR channelNum)
Definition: LISAspecifics.c:235
GPS2REAL8
#define GPS2REAL8(gps)
Definition: LISAspecifics.c:45
LISAarmT
LISAarmT
Translate TDI arm indices to C-indexing.
Definition: LISAspecifics.h:49
LISA_ARM2
@ LISA_ARM2
Definition: LISAspecifics.h:51
LISA_ARM1
@ LISA_ARM1
Definition: LISAspecifics.h:50
LISA_ARM3
@ LISA_ARM3
Definition: LISAspecifics.h:52
LISA_TIME_ORIGIN
#define LISA_TIME_ORIGIN
Definition: LISAspecifics.h:45
e
double e
XLALTensorSquareVector3
int XLALTensorSquareVector3(SymmTensor3 *vxv, REAL4 v[3])
Compute the "squared-tensor" v x v for given vector v, the result is returned in a "detectorTensor" s...
Definition: DetectorStates.c:140
XLALAddSymmTensor3s
int XLALAddSymmTensor3s(SymmTensor3 *sum, const SymmTensor3 *aT, const SymmTensor3 *bT)
Convenience function for adding two SymmTensor3s: aT + bT NOTE: it is safe to have sum point to the s...
Definition: DetectorStates.c:187
XLALSubtractSymmTensor3s
int XLALSubtractSymmTensor3s(SymmTensor3 *diff, const SymmTensor3 *aT, const SymmTensor3 *bT)
Convenience function for subtracting two SymmTensor3s: aT - bT NOTE: it is safe to have diff point to...
Definition: DetectorStates.c:211
XLALScaleSymmTensor3
int XLALScaleSymmTensor3(SymmTensor3 *mult, const SymmTensor3 *aT, REAL4 factor)
Convenience function for multiplying a SymmTensor3 by a scalar factor.
Definition: DetectorStates.c:235
Detector3Arms
DetectorArm Detector3Arms[3]
used to allow functions some type/size checking
Definition: DetectorStates.h:102
LAL_PI
#define LAL_PI
LAL_YRSID_SI
#define LAL_YRSID_SI
LAL_TWOPI
#define LAL_TWOPI
LAL_SQRT2
#define LAL_SQRT2
BOOLEAN
unsigned char BOOLEAN
REAL8
double REAL8
XLAL_INIT_DECL
#define XLAL_INIT_DECL(var,...)
crectf
#define crectf(re, im)
CHAR
char CHAR
UINT4
uint32_t UINT4
COMPLEX8
float complex COMPLEX8
REAL4
float REAL4
LALDETECTORTYPE_ABSENT
LALDETECTORTYPE_ABSENT
XLALRegisterSpecialCWDetector
int XLALRegisterSpecialCWDetector(const LALDetector *specialDetector)
Register a special detector for use with CW codes.
Definition: SFTnaming.c:73
XLALSinCosLUT
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
xlalErrno
#define xlalErrno
XLAL_ERROR
#define XLAL_ERROR(...)
XLAL_CHECK
#define XLAL_CHECK(assertion,...)
XLALPrintError
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EFUNC
XLAL_EINVAL
XLAL_EINVAL
y
list y
eta
eta
n
n
lalpulsar_knope.x
x
Definition: lalpulsar_knope.py:124
make_frame_cache.i
int i
Definition: make_frame_cache.py:77
test_ssbtodetector.LAL_AU_SI
int LAL_AU_SI
Definition: test_ssbtodetector.py:11
test_ssbtodetector.LAL_C_SI
int LAL_C_SI
Definition: test_ssbtodetector.py:9
CmplxDetectorTensor
The 'detector tensor' for a GW-detector: symmetric 3x3 matrix, storing only the upper triangle.
Definition: LISAspecifics.h:61
CmplxDetectorTensor::re
SymmTensor3 re
tensor holding real-parts of all components
Definition: LISAspecifics.h:62
CmplxDetectorTensor::im
SymmTensor3 im
tensor holding imaginary-parts of all components
Definition: LISAspecifics.h:63
DetectorArm
Struct containing pre-computed quantites describing a single detector arm: unit-vector along detector...
Definition: DetectorStates.h:96
DetectorArm::basisT
SymmTensor3 basisT
arm "basis-tensor" (n x n)
Definition: DetectorStates.h:98
DetectorArm::n
REAL4 n[3]
unit vector pointing along this arm
Definition: DetectorStates.h:97
DetectorArm::armlength_c
REAL4 armlength_c
armlengths in seconds L / c
Definition: DetectorStates.h:99
DetectorState
State-info about position, velocity and LMST of a detector together with corresponding EarthState.
Definition: DetectorStates.h:127
DetectorState::earthState
EarthState earthState
EarthState information.
Definition: DetectorStates.h:132
DetectorState::detArms
Detector3Arms detArms
include up to three arms to allow describing LISA
Definition: DetectorStates.h:133
DetectorState::tGPS
LIGOTimeGPS tGPS
GPS timestamps corresponding to this entry.
Definition: DetectorStates.h:128
EarthState::posNow
REAL8 posNow[3]
Cartesian coords of Earth's center at tgps, extrapolated from JPL DE405 ephemeris; units= sec.
Definition: LALBarycenter.h:171
FreqSkypos_t
Convenience container for precomputed pi f L/c and skyposition vector.
Definition: LISAspecifics.h:69
FreqSkypos_t::skyposV
REAL8 skyposV[3]
unit vector pointing to skyposition of source
Definition: LISAspecifics.h:71
FreqSkypos_t::Freq
REAL4 Freq
signal frequency
Definition: LISAspecifics.h:70
LALDetector::frDetector
LALFrDetector frDetector
LALFrDetector::name
CHAR name[LALNameLength]
SymmTensor3
A symmetric 3x3 tensor (such as detector-tensors), storing only the upper triangle.
Definition: DetectorStates.h:67