LALPulsar  6.1.0.1-89842e6
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 
38 static void LALComputeAM (LALStatus *, AMCoeffs *coe, LIGOTimeGPS *ts, AMCoeffsParams *params);
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  */
54 void 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,
132  HOUGHPeakGramVector *pgV,
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,
217  HOUGHPeakGramVector *pgV,
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  */
301 void 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 
474 void 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 
580 void 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 
633 void 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 
668 void 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 */
885  ASSERT (weightV->length == catalog->length, status, LALHOUGHH_ESZMM, LALHOUGHH_MSGESZMM);
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