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
DriveHough.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2005 Badri Krishnan, Alicia Sintes
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 <lal/LALHough.h>
21
22/** \cond DONT_DOXYGEN */
23
24#ifdef __GNUC__
25#define UNUSED __attribute__ ((unused))
26#else
27#define UNUSED
28#endif
29
30/** \endcond */
31
32/*
33 * The functions that make up the guts of this module
34 */
35
36#define SQ(x) (x) * (x)
37
39
40/** \addtogroup LALHough_h */
41/** @{ */
42
43/**
44 * constructs the space of <tt>phmd</tt> <tt>PHMDVectorSequence *phmdVS</tt>, given a
45 * HOUGHPeakGramVector *pgV and HOUGHptfLUTVector *lutV.
46 * The minimum frequency bin present corresponds to <tt>phmdVS->fBinMin</tt> and the
47 * total number of different frequencies is
48 * <tt>phmdVS->nfSize</tt>. At this moment the \c fBinMin line corresponds to the
49 * first row of the cylinder and <tt>phmdVS->breakLine</tt> is set to zero.
50 * <tt>phmdVS->breakLine</tt>
51 * \f$ \in\, [0, \f$ \c nfSize) is <em>the pointer</em> which
52 * identifies the position of the \c fBinMin row in the circular-cylinder buffer.
53 */
54void LALHOUGHConstructSpacePHMD (LALStatus *status, /**< pointer to LALStatus structure */
55 PHMDVectorSequence *phmdVS, /**< Cylindrical buffer of PHMDs */
56 HOUGHPeakGramVector *pgV, /**< Vetor of peakgrams */
57 HOUGHptfLUTVector *lutV /**< vector of look up tables */)
58{
59
60 UINT4 k,j;
61 UINT4 nfSize; /* number of different frequencies */
62 UINT4 length; /* number of elements for each frequency */
63 UINT8 fBinMin; /* present minimum frequency bin */
64 UINT8 fBin; /* present frequency bin */
65
66
67 /* --------------------------------------------- */
70
71
72 /* Make sure the arguments are not NULL: */
76 /* ------------------------------------------- */
77
81 /* ------------------------------------------- */
82
83 /* Make sure there is no size mismatch */
84 ASSERT (pgV->length == lutV->length, status,
86 ASSERT (pgV->length == phmdVS->length, status,
88 /* ------------------------------------------- */
89
90 /* Make sure there are elements to be computed*/
93
94 /* at the beggining, the fBinMin line corresponds to the first row */
95 phmdVS->breakLine = 0; /* mark [0,nfSize) (of the circular buffer)
96 pointing to the starting of the fBinMin line */
97
98 length = phmdVS->length;
99 nfSize = phmdVS->nfSize;
100 fBinMin = phmdVS->fBinMin;
101 phmdVS->deltaF = lutV->lut[0].deltaF; /* frequency resolution */
102
103 for ( k=0; k<length; ++k ){
104
105 /* make sure all deltaF are consistent */
106 ASSERT (phmdVS->deltaF == lutV->lut[k].deltaF,
108
109 fBin = fBinMin;
110
111 for ( j=0; j< nfSize; ++j ){
112 phmdVS->phmd[ j*length+k ].fBin = fBin;
113
114 TRY( LALHOUGHPeak2PHMD(status->statusPtr,
115 &(phmdVS->phmd[ j*length+k ]),
116 &(lutV->lut[k]), &(pgV->pg[k]) ), status);
117 ++fBin;
118 }
119 }
120
121
123 /* normal exit */
124 RETURN (status);
125}
126
127/**
128 * This function updates the space of <tt>phmd</tt> increasing the frequency <tt>phmdVS->fBinMin</tt> by one.
129 */
131 PHMDVectorSequence *phmdVS,
133 HOUGHptfLUTVector *lutV)
134{
135 UINT4 k,breakLine;
136 UINT4 nfSize; /* number of different frequencies */
137 UINT4 length; /* number of elements for each frequency */
138 UINT8 fBinMin; /* minimum frequency bin */
139 UINT8 fBin; /* present frequency bin */
140
141
142 /* --------------------------------------------- */
145
146
147 /* Make sure the arguments are not NULL: */
151 /* ------------------------------------------- */
152
156 /* ------------------------------------------- */
157
158 /* Make sure there is no size mismatch */
159 ASSERT (pgV->length == lutV->length, status,
161 ASSERT (pgV->length == phmdVS->length, status,
163 /* ------------------------------------------- */
164
165 /* Make sure there are elements to be computed*/
168 /* ------------------------------------------- */
169
170
171 length = phmdVS->length;
172 nfSize = phmdVS->nfSize;
173
174 breakLine = phmdVS->breakLine; /* old Break Line */
175 fBinMin = phmdVS->fBinMin; /* initial frequency value */
176
177 /* Make sure initial breakLine is in [0,nfSize) */
178 ASSERT ( breakLine < nfSize, status, LALHOUGHH_EVAL, LALHOUGHH_MSGEVAL);
179
180 /* ------------------------------------------- */
181
182 /* Updating the space of PHMD increasing frequency */
183
184 fBin = fBinMin + nfSize;
185
186 for ( k=0; k<length; ++k ){
187 /* make sure all deltaF are consistent */
188 ASSERT (phmdVS->deltaF == lutV->lut[k].deltaF,
190
191 phmdVS->phmd[ breakLine*length+k ].fBin = fBin;
192 TRY( LALHOUGHPeak2PHMD(status->statusPtr,
193 &(phmdVS->phmd[ breakLine*length+k ]),
194 &(lutV->lut[k]), &(pgV->pg[k]) ), status);
195 }
196
197 /* Shift fBinMin and its mark */
198 ++phmdVS->fBinMin;
199
200 phmdVS->breakLine = (breakLine +1) % nfSize;
201 /* mark [0,nfSize) (of the circular buffer, modulus nfSize)
202 pointing to the starting of the new fBinMin line */
203
205 /* normal exit */
206 RETURN (status);
207}
208
209
210/**
211 * Function for shifting the cylindrical buffer of PHMDs down by one
212 * frequency bin -- the highest frequency bin is dropped and an
213 * extra frequency bin is added at the lowest frequency
214 */
216 PHMDVectorSequence *phmdVS,
218 HOUGHptfLUTVector *lutV)
219{
220 UINT4 k,breakLine;
221 UINT4 nfSize; /* number of different frequencies */
222 UINT4 length; /* number of elements for each frequency */
223
224 UINT8 fBin; /* present frequency bin */
225
226 /* --------------------------------------------- */
229
230
231 /* Make sure the arguments are not NULL: */
235 /* ------------------------------------------- */
236
240 /* ------------------------------------------- */
241
242 /* Make sure there is no size mismatch */
243 ASSERT (pgV->length == lutV->length, status,
245 ASSERT (pgV->length == phmdVS->length, status,
247 /* ------------------------------------------- */
248
249 /* Make sure there are elements to be computed*/
252 /* ------------------------------------------- */
253
254 length = phmdVS->length;
255 nfSize = phmdVS->nfSize;
256
257 breakLine = phmdVS->breakLine; /* old Break Line */
258
259 /* Make sure initial breakLine is in [0,nfSize) */
260 ASSERT ( breakLine < nfSize, status, LALHOUGHH_EVAL, LALHOUGHH_MSGEVAL);
261
262 /* ------------------------------------------- */
263
264 /* Updating the space of PHMD decreasing frequency */
265
266 /* Shift fBinMin and its mark */
267 fBin = --phmdVS->fBinMin; /* initial frequency value */
268
269 phmdVS->breakLine = (breakLine + nfSize- 1) % nfSize;
270 /* mark [0,nfSize) (of the circular buffer, modulus nfSize)
271 pointing to the starting of the new fBinMin line */
272
273 breakLine = phmdVS->breakLine; /* the new Break Line */
274
275 for ( k=0; k<length; ++k ){
276 /* make sure all deltaF are consistent */
277 ASSERT (phmdVS->deltaF == lutV->lut[k].deltaF,
279
280 phmdVS->phmd[ breakLine*length+k ].fBin = fBin;
281 TRY( LALHOUGHPeak2PHMD(status->statusPtr,
282 &(phmdVS->phmd[ breakLine*length+k ]),
283 &(lutV->lut[k]), &(pgV->pg[k]) ), status);
284 }
285
286
287
289 /* normal exit */
290 RETURN (status);
291}
292
293
294/**
295 * Given PHMDVectorSequence *phmdVS, the space of \c phmd, and
296 * UINT8FrequencyIndexVector *freqInd, a structure containing the frequency
297 * indices of the \c phmd at different time stamps that have to be combined
298 * to form a Hough map, the function LALHOUGHConstructHMT() produces the
299 * total Hough map.
300 */
301void LALHOUGHConstructHMT (LALStatus *status, /**< pointer to LALStatus structure */
302 HOUGHMapTotal *ht, /**< The output hough map */
303 UINT8FrequencyIndexVector *freqInd,/**< time-frequency trajectory */
304 PHMDVectorSequence *phmdVS /**< set of partial hough map derivatives */)
305{
306
307
308 UINT4 k,j;
309 UINT4 breakLine;
310 UINT4 nfSize; /* number of different frequencies */
311 UINT4 length; /* number of elements for each frequency */
312 UINT8 fBinMin; /* present minimum frequency bin */
313 INT8 fBin; /* present frequency bin */
314 UINT2 xSide,ySide;
315
316 HOUGHMapDeriv hd; /* the Hough map derivative */
317
318 /* --------------------------------------------- */
321
322 /* Make sure the arguments are not NULL: */
326 /* ------------------------------------------- */
327
330 /* ------------------------------------------- */
331
332 /* Make sure there is no size mismatch */
333 ASSERT (freqInd->length == phmdVS->length, status,
335 ASSERT (freqInd->deltaF == phmdVS->deltaF, status,
337 /* ------------------------------------------- */
338
339 /* Make sure there are elements */
342 /* ------------------------------------------- */
343
344 /* Make sure the ht map contains some pixels */
347
348 length = phmdVS->length;
349 nfSize = phmdVS->nfSize;
350
351 fBinMin = phmdVS->fBinMin; /* initial frequency value od the cilinder*/
352
353 breakLine = phmdVS->breakLine;
354
355 /* number of physical pixels */
356 xSide = ht->xSide;
357 ySide = ht->ySide;
358
359 /* Make sure initial breakLine is in [0,nfSize) */
360 ASSERT ( breakLine < nfSize, status, LALHOUGHH_EVAL, LALHOUGHH_MSGEVAL);
361
362 /* ------------------------------------------- */
363
364 /* Initializing hd map and memory allocation */
365 hd.xSide = xSide;
366 hd.ySide = ySide;
367 hd.map = (HoughDT *)LALMalloc(ySide*(xSide+1)*sizeof(HoughDT));
368 if (hd. map == NULL) {
370 }
371 /* ------------------------------------------- */
372
373
374 TRY( LALHOUGHInitializeHD(status->statusPtr, &hd), status);
375 for ( k=0; k<length; ++k ){
376 /* read the frequency index and make sure is in the proper interval*/
377 fBin =freqInd->data[k] -fBinMin;
378
379 ASSERT ( fBin < nfSize, status, LALHOUGHH_EVAL, LALHOUGHH_MSGEVAL);
381
382 /* find index */
383 j = (fBin + breakLine) % nfSize;
384
385 /* Add the corresponding PHMD to HD */
386 TRY( LALHOUGHAddPHMD2HD(status->statusPtr,
387 &hd, &(phmdVS->phmd[j*length+k]) ), status);
388 }
389
390 TRY( LALHOUGHIntegrHD2HT(status->statusPtr, ht, &hd), status);
391
392 /* Free memory and exit */
393 LALFree(hd.map);
394
396 /* normal exit */
397 RETURN (status);
398}
399
400
401
402/**
403 * This function computes the corresponding frequency bin of
404 * a \c phmd <tt>UINT8 *fBinMap</tt> for a given intrinsic
405 * frequency bin of a source <tt>UINT8 *f0Bin</tt>, and information regarding the
406 * time and the residual spin down parameters HOUGHResidualSpinPar *rs.
407 */
409 UINT8 *fBinMap,
410 UINT8 *f0Bin,
412{
413
414 UINT4 i;
415 INT4 shiftFBin;
416 REAL8 shiftF;
417
418 UINT4 spinOrder;
419 REAL8 *spinF;
420 REAL8 timeDiff; /* T(t)-T(t0) */
421 REAL8 timeDiffProd;
422 REAL8 deltaF; /* df=1/TCOH */
423 /* --------------------------------------------- */
426
427 /* Make sure the arguments are not NULL: */
431 /* ------------------------------------------- */
432
433 /* Make sure the Input/Output pointers are not the same */
434 ASSERT (fBinMap != f0Bin, status, LALHOUGHH_ESAME, LALHOUGHH_MSGESAME);
435
436 shiftFBin = 0;
437 shiftF = 0.0;
438
439 spinOrder = rs->spinRes.length;
440
441 if(spinOrder){
443 timeDiff = rs->timeDiff;
444 timeDiffProd = timeDiff;
445
446 deltaF = rs->deltaF;
447 spinF = rs->spinRes.data;
448
449 for (i=0; i<spinOrder; ++i ){
450 shiftF += spinF[i] * timeDiffProd;
451 timeDiffProd *= timeDiff;
452 }
453 shiftFBin = rint( shiftF/deltaF) ; /* positive or negative */
454 }
455
456 *fBinMap = *f0Bin + shiftFBin;
457
459 /* normal exit */
460 RETURN (status);
461}
462
463
464
465
466
467/**
468 * Calculates the total hough map for a given trajectory in the
469 * time-frequency plane and a set of partial hough map derivatives allowing
470 * each PHMD to have a different weight factor to account for varying
471 * sensitivity at different sky-locations.
472 */
473
474void LALHOUGHConstructHMT_W (LALStatus *status, /**< pointer to LALStatus structure */
475 HOUGHMapTotal *ht, /**< The output hough map */
476 UINT8FrequencyIndexVector *freqInd, /**< time-frequency trajectory */
477 PHMDVectorSequence *phmdVS /**< set of partial hough map derivatives */)
478{
479
480
481 UINT4 k,j;
482 UINT4 breakLine;
483 UINT4 nfSize; /* number of different frequencies */
484 UINT4 length; /* number of elements for each frequency */
485 UINT8 fBinMin; /* present minimum frequency bin */
486 INT8 fBin; /* present frequency bin */
487 UINT2 xSide,ySide;
488
489 HOUGHMapDeriv hd; /* the Hough map derivative */
490
491 /* --------------------------------------------- */
494
495 /* Make sure the arguments are not NULL: */
499 /* ------------------------------------------- */
500
503 /* ------------------------------------------- */
504
505 /* Make sure there is no size mismatch */
506 ASSERT (freqInd->length == phmdVS->length, status,
508 ASSERT (freqInd->deltaF == phmdVS->deltaF, status,
510 /* ------------------------------------------- */
511
512 /* Make sure there are elements */
515 /* ------------------------------------------- */
516
517 /* Make sure the ht map contains some pixels */
520
521 length = phmdVS->length;
522 nfSize = phmdVS->nfSize;
523
524 fBinMin = phmdVS->fBinMin; /* initial frequency value od the cilinder*/
525
526 breakLine = phmdVS->breakLine;
527
528 /* number of physical pixels */
529 xSide = ht->xSide;
530 ySide = ht->ySide;
531
532 /* Make sure initial breakLine is in [0,nfSize) */
533 ASSERT ( breakLine < nfSize, status, LALHOUGHH_EVAL, LALHOUGHH_MSGEVAL);
534
535 /* ------------------------------------------- */
536
537 /* Initializing hd map and memory allocation */
538 hd.xSide = xSide;
539 hd.ySide = ySide;
540 hd.map = (HoughDT *)LALMalloc(ySide*(xSide+1)*sizeof(HoughDT));
541 if (hd. map == NULL) {
543 }
544
545 /* ------------------------------------------- */
546
547 TRY( LALHOUGHInitializeHD(status->statusPtr, &hd), status);
548 for ( k=0; k<length; ++k ){
549 /* read the frequency index and make sure is in the proper interval*/
550 fBin =freqInd->data[k] -fBinMin;
551
552 ASSERT ( fBin < nfSize, status, LALHOUGHH_EVAL, LALHOUGHH_MSGEVAL);
554
555 /* find index */
556 j = (fBin + breakLine) % nfSize;
557
558 /* Add the corresponding PHMD to HD */
559 TRY( LALHOUGHAddPHMD2HD_W(status->statusPtr,
560 &hd, &(phmdVS->phmd[j*length+k]) ), status);
561 }
562
563 TRY( LALHOUGHIntegrHD2HT(status->statusPtr, ht, &hd), status);
564
565 /* Free memory and exit */
566 LALFree(hd.map);
567
569 /* normal exit */
570 RETURN (status);
571}
572
573
574
575/**
576 * Adds weight factors for set of partial hough map derivatives -- the
577 * weights must be calculated outside this function.
578 */
579
580void LALHOUGHWeighSpacePHMD (LALStatus *status, /**< pointer to LALStatus structure */
581 PHMDVectorSequence *phmdVS, /**< partial hough map derivatives */
582 REAL8Vector *weightV /**< vector of weights */)
583{
584
585 UINT4 k,j;
586 UINT4 nfSize; /* number of different frequencies */
587 UINT4 length; /* number of elements for each frequency */
588 /* --------------------------------------------- */
589
592
593
594 /* Make sure the arguments are not NULL: */
597 /* ------------------------------------------- */
598
601 /* ------------------------------------------- */
602
603 /* Make sure there is no size mismatch */
604 ASSERT (weightV->length == phmdVS->length, status,
606 /* ------------------------------------------- */
607
608 /* Make sure there are elements to be computed*/
611
612
613 length = phmdVS->length;
614 nfSize = phmdVS->nfSize;
615
616 /* weigh the phmds according to weightV */
617 for ( k=0; k<length; ++k ){
618 for ( j=0; j< nfSize; ++j ){
619 phmdVS->phmd[ j*length+k ].weight = (HoughDT)weightV->data[k];
620 }
621 }
622
623
625 /* normal exit */
626 RETURN (status);
627}
628
629
630
631/** Initializes weight factors to unity */
632
633void LALHOUGHInitializeWeights (LALStatus *status, /**< pointer to LALStatus structure */
634 REAL8Vector *weightV /**< vector of weights */)
635{
636
637 UINT4 j, length;
638
639 /* --------------------------------------------- */
642
643
644 /* Make sure the arguments are not NULL: */
647 /* ------------------------------------------- */
648
649 /* Make sure there are elements to be computed*/
651
652 length = weightV->length;
653
654 for (j=0; j<length; j++) {
655 weightV->data[j] = 1.0;
656 }
657
659 /* normal exit */
660 RETURN (status);
661}
662
663
664
665
666/** Normalizes weight factors so that their sum is N */
667
668void LALHOUGHNormalizeWeights (LALStatus *status, /**< pointer to LALStatus structure */
669 REAL8Vector *weightV /**< vector of weights */)
670{
671
672 UINT4 j, length;
673 REAL8 sum;
674
675 /* --------------------------------------------- */
678
679
680 /* Make sure the arguments are not NULL: */
683 /* ------------------------------------------- */
684
685 /* Make sure there are elements to be computed*/
687
688 length = weightV->length;
689
690 /* calculate sum of weights */
691 sum = 0.0;
692 for (j=0; j<length; j++) {
693 sum += weightV->data[j];
694 }
695
696 /* normalize weights */
697 for (j=0; j<length; j++) {
698 weightV->data[j] *= length/sum;
699 }
700
702 /* normal exit */
703 RETURN (status);
704}
705
706
707
708/**
709 * Computes weight factors arising from amplitude modulation -- it multiplies
710 * an existing weight vector
711 */
713 REAL8Vector *weightV,
714 LIGOTimeGPSVector *timeV,
717 REAL8 alpha,
718 REAL8 delta)
719{
720
721 UINT4 length, j;
722
723 /* amplitude modulation stuff */
724 REAL4Vector *aVec, *bVec;
725 REAL8 A, B, a, b;
727 AMCoeffsParams *amParams;
728 EarthState earth;
729 BarycenterInput baryinput;
730
731 /* --------------------------------------------- */
734
735
736 /* Make sure the arguments are not NULL: */
741
744 /* ------------------------------------------- */
745
746 /* Make sure there is no size mismatch */
748 /* ------------------------------------------- */
749
750 /* Make sure there are elements to be computed*/
752
753 length = timeV->length;
754
755
756 /* detector location */
757 baryinput.site.location[0] = detector->location[0]/LAL_C_SI;
758 baryinput.site.location[1] = detector->location[1]/LAL_C_SI;
759 baryinput.site.location[2] = detector->location[2]/LAL_C_SI;
760 baryinput.dInv = 0.e0;
761
762 /* alpha and delta must come from the skypatch */
763 /* for now set it to something arbitrary */
764 /* baryinput.alpha = 0.0; */
765 /* baryinput.delta = 0.0; */
766
767 baryinput.alpha = alpha;
768 baryinput.delta = delta;
769
770 /* Allocate space for amParams stucture */
771 /* Here, amParams->das is the Detector and Source info */
772 amParams = (AMCoeffsParams *)LALMalloc(sizeof(AMCoeffsParams));
773 if (amParams == NULL) {
775 }
776
777 amParams->das = (LALDetAndSource *)LALMalloc(sizeof(LALDetAndSource));
778 if (amParams->das == NULL) {
780 }
781
782 amParams->das->pSource = (LALSource *)LALMalloc(sizeof(LALSource));
783 if (amParams->das->pSource == NULL) {
785 }
786
787 /* Fill up AMCoeffsParams structure */
788 amParams->baryinput = &baryinput;
789 amParams->earth = &earth;
790 amParams->edat = edat;
791 amParams->das->pDetector = detector;
792 /* make sure alpha and delta are correct */
793 amParams->das->pSource->equatorialCoords.latitude = baryinput.delta;
794 amParams->das->pSource->equatorialCoords.longitude = baryinput.alpha;
795 amParams->das->pSource->orientation = 0.0;
797 amParams->polAngle = amParams->das->pSource->orientation ; /* These two have to be the same!!*/
798
799 /* allocate memory for a[i] and b[i] */
800 /* TRY( LALSCreateVector( status, &(amc.a), length), status); */
801 /* TRY( LALSCreateVector( status, &(amc.b), length), status); */
802
803 amc.a = NULL;
804 amc.a = (REAL4Vector *)LALMalloc(sizeof(REAL4Vector));
805 if (amc.a == NULL) {
807 }
808
809 amc.a->length = length;
810 amc.a->data = (REAL4 *)LALMalloc(length*sizeof(REAL4));
811 if (amc.a->data == NULL) {
813 }
814
815 amc.b = NULL;
816 amc.b = (REAL4Vector *)LALMalloc(sizeof(REAL4Vector));
817 if (amc.b == NULL) {
819 }
820
821 amc.b->length = length;
822 amc.b->data = (REAL4 *)LALMalloc(length*sizeof(REAL4));
823 if (amc.b->data == NULL) {
825 }
826
827 TRY (LALComputeAM( status->statusPtr, &amc, timeV->data, amParams), status);
828 aVec = amc.a; /* a and b are as defined in JKS */
829 bVec = amc.b;
830 A = amc.A; /* note A is twice average of a[i]^2 */
831 B = amc.B; /* note B is twice average of b[i]^2 */
832
833 for(j=0; j<length; j++){
834 a = aVec->data[j];
835 b = bVec->data[j];
836 weightV->data[j] *= 2.0*(a*a + b*b)/(A+B);
837 }
838
839 /* normalize weights */
840 TRY( LALHOUGHNormalizeWeights( status->statusPtr, weightV ), status);
841
842 /* TRY( LALSDestroyVector( status, &(amc.a)), status); */
843 /* TRY( LALSDestroyVector( status, &(amc.b)), status); */
844 LALFree(amc.a->data);
845 LALFree(amc.b->data);
846 LALFree(amc.a);
847 LALFree(amc.b);
848
849 LALFree(amParams->das->pSource);
850 LALFree(amParams->das);
851 LALFree(amParams);
852
854 /* normal exit */
855 RETURN (status);
856}
857
858
859
860
861/**
862 * Computes weight factors arising from amplitude modulation -- it multiplies
863 * an existing weight vector
864 */
866 REAL8Vector *weightV,
867 SFTCatalog *catalog,
869 REAL8 UNUSED alpha,
870 REAL8 UNUSED delta)
871{
872
873 /* --------------------------------------------- */
876
880
883
884 /* Make sure there is no size mismatch */
886
887 (void)weightV;
888 (void)catalog;
889 (void)edat;
890
892 /* normal exit */
893 RETURN (status);
894}
895
896
897/**
898 * Original antenna-pattern function by S Berukoff
899 */
901 AMCoeffs *coe,
902 LIGOTimeGPS *ts,
904{
905
906 REAL4 zeta; /* sine of angle between detector arms */
907 INT4 i; /* temporary loop index */
908 LALDetAMResponse response; /* output of LALComputeDetAMResponse */
909
910 REAL4 sumA2=0.0;
911 REAL4 sumB2=0.0;
912 REAL4 sumAB=0.0; /* variables to store scalar products */
913 INT4 length=coe->a->length; /* length of input time series */
914
915 REAL4 cos2psi;
916 REAL4 sin2psi; /* temp variables */
917
920
921 /* Must put an ASSERT checking that ts vec and coe vec are same length */
922
923 /* Compute the angle between detector arms, then the reciprocal */
924 {
925 LALFrDetector det = params->das->pDetector->frDetector;
926 zeta = 1.0/(sin(det.xArmAzimuthRadians - det.yArmAzimuthRadians));
927 if(params->das->pDetector->type == LALDETECTORTYPE_CYLBAR) zeta=1.0;
928 }
929
930 cos2psi = cos(2.0 * params->polAngle);
931 sin2psi = sin(2.0 * params->polAngle);
932
933 /* Note the length is the same for the b vector */
934 for(i=0; i<length; ++i)
935 {
936 REAL4 *a = coe->a->data;
937 REAL4 *b = coe->b->data;
938
939 /* Compute F_plus, F_cross */
940 LALComputeDetAMResponse(status->statusPtr, &response, params->das, &ts[i]);
941
942 /* Compute a, b from JKS eq 10,11
943 * a = zeta * (F_plus*cos(2\psi)-F_cross*sin(2\psi))
944 * b = zeta * (F_cross*cos(2\psi)+Fplus*sin(2\psi))
945 */
946 a[i] = zeta * (response.plus*cos2psi-response.cross*sin2psi);
947 b[i] = zeta * (response.cross*cos2psi+response.plus*sin2psi);
948
949 /* Compute scalar products */
950 sumA2 += SQ(a[i]); /* A */
951 sumB2 += SQ(b[i]); /* B */
952 sumAB += (a[i]) * (b[i]); /* C */
953 }
954
955 {
956 /* Normalization factor */
957 REAL8 L = 2.0/(REAL8)length;
958
959 /* Assign output values and normalise */
960 coe->A = L*sumA2;
961 coe->B = L*sumB2;
962 coe->C = L*sumAB;
963 coe->D = ( coe->A * coe->B - SQ(coe->C) );
964 /* protection against case when AB=C^2 */
965 if(coe->D == 0) coe->D=1.0e-9;
966 }
967
968 /* Normal exit */
969
971 RETURN(status);
972
973} /* LALComputeAM() */
974
975/** @} */
#define SQ(x)
Definition: DriveHough.c:36
REAL8 zeta
int j
int k
#define LALMalloc(n)
#define LALFree(p)
static double double delta
#define ABORT(statusptr, code, mesg)
#define TRY(func, statusptr)
#define ATTATCHSTATUSPTR(statusptr)
#define ASSERT(assertion, statusptr, code, mesg)
#define DETATCHSTATUSPTR(statusptr)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
#define L(i, j)
UINT2 A
Definition: SFTnaming.c:46
UINT2 B
Definition: SFTnaming.c:47
AMCoeffs amc
void LALComputeDetAMResponse(LALStatus *status, LALDetAMResponse *pResponse, const LALDetAndSource *pDetAndSrc, const LIGOTimeGPS *gps)
void LALHOUGHAddPHMD2HD_W(LALStatus *status, HOUGHMapDeriv *hd, HOUGHphmd *phmd)
Adds a hough map derivative into a total hough map derivative taking into account the weight of the p...
Definition: HoughMap.c:231
void LALHOUGHInitializeHD(LALStatus *status, HOUGHMapDeriv *hd)
This function initializes the Hough map derivative space HOUGHMapDeriv *hd to zero.
Definition: HoughMap.c:31
void LALHOUGHAddPHMD2HD(LALStatus *status, HOUGHMapDeriv *hd, HOUGHphmd *phmd)
Given an initial Hough map derivative HOUGHMapDeriv *hd and a representation of a phmd HOUGHphmd *phm...
Definition: HoughMap.c:123
void LALHOUGHIntegrHD2HT(LALStatus *status, HOUGHMapTotal *ht, HOUGHMapDeriv *hd)
This function constructs a total Hough map HOUGHMapTotal *ht from its derivative HOUGHMapDeriv *hd by...
Definition: HoughMap.c:353
uint64_t UINT8
double REAL8
#define XLAL_INIT_DECL(var,...)
int64_t INT8
uint16_t UINT2
uint32_t UINT4
int32_t INT4
float REAL4
LALDETECTORTYPE_CYLBAR
static void LALComputeAM(LALStatus *, AMCoeffs *coe, LIGOTimeGPS *ts, AMCoeffsParams *params)
Original antenna-pattern function by S Berukoff.
Definition: DriveHough.c:900
#define LALHOUGHH_EVAL
Definition: LALHough.h:116
void LALHOUGHConstructHMT_W(LALStatus *status, HOUGHMapTotal *ht, UINT8FrequencyIndexVector *freqInd, PHMDVectorSequence *phmdVS)
Calculates the total hough map for a given trajectory in the time-frequency plane and a set of partia...
Definition: DriveHough.c:474
#define LALHOUGHH_MSGESZMM
Definition: LALHough.h:121
#define LALHOUGHH_MSGESIZE
Definition: LALHough.h:120
#define LALHOUGHH_ESZMM
Definition: LALHough.h:112
#define LALHOUGHH_MSGEMEM
Definition: LALHough.h:126
#define LALHOUGHH_ESAME
Definition: LALHough.h:114
void LALHOUGHConstructSpacePHMD(LALStatus *status, PHMDVectorSequence *phmdVS, HOUGHPeakGramVector *pgV, HOUGHptfLUTVector *lutV)
constructs the space of phmd PHMDVectorSequence *phmdVS, given a HOUGHPeakGramVector *pgV and HOUGHpt...
Definition: DriveHough.c:54
void LALHOUGHupdateSpacePHMDdn(LALStatus *status, PHMDVectorSequence *phmdVS, HOUGHPeakGramVector *pgV, HOUGHptfLUTVector *lutV)
Function for shifting the cylindrical buffer of PHMDs down by one frequency bin – the highest frequen...
Definition: DriveHough.c:215
#define LALHOUGHH_MSGESAME
Definition: LALHough.h:123
void LALHOUGHNormalizeWeights(LALStatus *status, REAL8Vector *weightV)
Normalizes weight factors so that their sum is N.
Definition: DriveHough.c:668
#define LALHOUGHH_ENULL
Definition: LALHough.h:110
#define LALHOUGHH_MSGENULL
Definition: LALHough.h:119
void LALHOUGHComputeFBinMap(LALStatus *status, UINT8 *fBinMap, UINT8 *f0Bin, HOUGHResidualSpinPar *rs)
This function computes the corresponding frequency bin of a phmd UINT8 *fBinMap for a given intrinsic...
Definition: DriveHough.c:408
void LALHOUGHComputeAMWeights(LALStatus *status, REAL8Vector *weightV, LIGOTimeGPSVector *timeV, LALDetector *detector, EphemerisData *edat, REAL8 alpha, REAL8 delta)
Computes weight factors arising from amplitude modulation – it multiplies an existing weight vector.
Definition: DriveHough.c:712
void LALHOUGHWeighSpacePHMD(LALStatus *status, PHMDVectorSequence *phmdVS, REAL8Vector *weightV)
Adds weight factors for set of partial hough map derivatives – the weights must be calculated outside...
Definition: DriveHough.c:580
void LALHOUGHInitializeWeights(LALStatus *status, REAL8Vector *weightV)
Initializes weight factors to unity.
Definition: DriveHough.c:633
#define LALHOUGHH_MSGEVAL
Definition: LALHough.h:125
void LALHOUGHComputeMultiIFOAMWeights(LALStatus *status, REAL8Vector *weightV, SFTCatalog *catalog, EphemerisData *edat, REAL8 UNUSED alpha, REAL8 UNUSED delta)
Computes weight factors arising from amplitude modulation – it multiplies an existing weight vector.
Definition: DriveHough.c:865
void LALHOUGHupdateSpacePHMDup(LALStatus *status, PHMDVectorSequence *phmdVS, HOUGHPeakGramVector *pgV, HOUGHptfLUTVector *lutV)
This function updates the space of phmd increasing the frequency phmdVS->fBinMin by one.
Definition: DriveHough.c:130
#define LALHOUGHH_ESIZE
Definition: LALHough.h:111
#define LALHOUGHH_EMEM
Definition: LALHough.h:117
void LALHOUGHConstructHMT(LALStatus *status, HOUGHMapTotal *ht, UINT8FrequencyIndexVector *freqInd, PHMDVectorSequence *phmdVS)
Given PHMDVectorSequence *phmdVS, the space of phmd, and UINT8FrequencyIndexVector *freqInd,...
Definition: DriveHough.c:301
REAL8 HoughDT
Hough Map derivative pixel type.
Definition: PHMD.h:121
void LALHOUGHPeak2PHMD(LALStatus *status, HOUGHphmd *phmd, HOUGHptfLUT *lut, HOUGHPeakGram *pg)
Construction of Partial-Hough-Map-Derivatives (phmd) given a peak-gram and the look-up-table.
Definition: Peak2PHMD.c:61
static const INT4 a
COORDINATESYSTEM_EQUATORIAL
int deltaF
ts
double alpha
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
This structure contains the parameters for the routine.
Definition: LALComputeAM.h:75
LALDetAndSource * das
det and source information
Definition: LALComputeAM.h:79
REAL4 polAngle
polarization angle
Definition: LALComputeAM.h:81
EarthState * earth
from XLALBarycenter()
Definition: LALComputeAM.h:77
EphemerisData * edat
the ephemerides
Definition: LALComputeAM.h:78
BarycenterInput * baryinput
data from Barycentring routine
Definition: LALComputeAM.h:76
Basic input structure to LALBarycenter.c.
REAL8 alpha
Source right ascension in ICRS J2000 coords (radians).
REAL8 dInv
1/(distance to source), in 1/sec.
REAL8 delta
Source declination in ICRS J2000 coords (radians)
LALDetector site
detector site info.
Basic output structure of LALBarycenterEarth.c.
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
This structure stores the Hough map derivative.
Definition: HoughMap.h:121
UINT2 ySide
number of physical pixels in the y direction
Definition: HoughMap.h:123
UINT2 xSide
number of physical pixels in the x direction
Definition: HoughMap.h:122
HoughDT * map
the pixel count derivatives; the number of elements to allocate is ySide*(xSide+1)*
Definition: HoughMap.h:124
This structure stores the Hough map.
Definition: HoughMap.h:130
UINT2 ySide
number of physical pixels in the y direction
Definition: HoughMap.h:142
UINT2 xSide
number of physical pixels in the x direction
Definition: HoughMap.h:141
This structure contains a vector of peak-grams (for the different time stamps)
Definition: LALHough.h:167
UINT4 length
number of elements
Definition: LALHough.h:168
HOUGHPeakGram * pg
the Peakgrams
Definition: LALHough.h:169
This structure stores the residual spin-down parameters at a given time.
Definition: LALHough.h:199
REAL8 deltaF
Frequency resolution; df=1/TCOH.
Definition: LALHough.h:200
REAL8Vector spinRes
length: Maximum order of spdwn parameter *data: pointer to residual Spin parameter set fk
Definition: LALHough.h:202
REAL8 timeDiff
: time difference
Definition: LALHough.h:201
UINT8 fBin
Frequency bin of this partial map derivative.
Definition: PHMD.h:142
HoughDT weight
First column border, containing the edge effects when clipping on a finite patch.
Definition: PHMD.h:152
REAL8 deltaF
Frequency resolution df=1/TCOH, where 1/TCOH is the coherent integration time used in teh demodulatio...
Definition: LUT.h:249
This structure contains a vector of partial look up tables (for the different time stamps)
Definition: LALHough.h:173
HOUGHptfLUT * lut
the partial Look Up Tables
Definition: LALHough.h:175
UINT4 length
number of elements
Definition: LALHough.h:174
LALSource * pSource
const LALDetector * pDetector
REAL8 location[3]
REAL4 xArmAzimuthRadians
REAL4 yArmAzimuthRadians
SkyPosition equatorialCoords
REAL8 orientation
A vector of 'timestamps' of type LIGOTimeGPS.
Definition: SFTfileIO.h:188
LIGOTimeGPS * data
array of timestamps
Definition: SFTfileIO.h:193
UINT4 length
number of timestamps
Definition: SFTfileIO.h:192
This structure contains a vector sequence of partial-Hough maps derivatives (for different time stamp...
Definition: LALHough.h:189
UINT8 fBinMin
frequency index of smallest intrinsic frequency in circular buffer
Definition: LALHough.h:192
UINT4 breakLine
Mark [0, nfSize) (of the circular buffer) pointing to the starting of the fBinMin line.
Definition: LALHough.h:194
UINT4 nfSize
number of different frequencies
Definition: LALHough.h:190
REAL8 deltaF
frequency resolution
Definition: LALHough.h:193
HOUGHphmd * phmd
the partial Hough map derivatives
Definition: LALHough.h:195
UINT4 length
number of elements for each frequency
Definition: LALHough.h:191
REAL4 * data
REAL8 * data
An "SFT-catalogue": a vector of SFTdescriptors, as returned by XLALSFTdataFind()
Definition: SFTfileIO.h:238
SFTDescriptor * data
array of data-entries describing matched SFTs
Definition: SFTfileIO.h:243
UINT4 length
number of SFTs in catalog
Definition: SFTfileIO.h:242
REAL8 longitude
REAL8 latitude
CoordinateSystem system
This structure stores the frequency indexes of the partial-Hough map derivatives at different time st...
Definition: LALHough.h:149
UINT8 * data
the frequency indexes
Definition: LALHough.h:152
REAL8 deltaF
frequency resolution
Definition: LALHough.h:151
UINT4 length
number of elements
Definition: LALHough.h:150