Line data Source code
1 :
2 : // -----------------------------------------------------------------------------
3 :
4 : /*
5 :
6 : CalStatsDerived.cc
7 :
8 : Description:
9 : ------------
10 : This file contains the member functions of the classes derived from CalStats.
11 :
12 : Classes:
13 : --------
14 : CalStatsReal - This class feeds real data to the CalStats base class.
15 : CalStatsAmp - This class converts complex data to amplitudes and initializes
16 : the CalStats base class.
17 : CalStatsPhase - This class converts complex data to phases and initializes the
18 : CalStats base class.
19 :
20 : Inhertited classes:
21 : -------------------
22 : CalStats - This class calculates statistics on CASA caltables.
23 :
24 : Modification history:
25 : ---------------------
26 : 2011 Nov 15 - Nick Elias, NRAO
27 : Initial version created with classes CalStatsAmp and
28 : CalStatsPhase.
29 : 2011 Dec 11 - Nick Elias, NRAO
30 : Class CalStatsReal added.
31 : 2012 Jan 25 - Nick Elias, NRAO
32 : Logging capability added. Error checking added.
33 :
34 : */
35 :
36 : // -----------------------------------------------------------------------------
37 : // Includes
38 : // -----------------------------------------------------------------------------
39 :
40 : #include <calanalysis/CalAnalysis/CalStatsDerived.h>
41 : #include <casacore/casa/Utilities/GenSort.h>
42 :
43 : // -----------------------------------------------------------------------------
44 : // Start of casa namespace
45 : // -----------------------------------------------------------------------------
46 :
47 : using namespace casacore;
48 : namespace casa {
49 :
50 : const uInt CalStatsPhase::NUM_ITER_UNWRAP = 50;
51 : const Double CalStatsPhase::NEW_RANGE_FACTOR = 5.0;
52 :
53 : // -----------------------------------------------------------------------------
54 : // Start of CalStatsReal class
55 : // -----------------------------------------------------------------------------
56 :
57 : /*
58 :
59 : CalStatsReal
60 :
61 : Description:
62 : ------------
63 : This class feeds real data to the CalStats base class.
64 :
65 : Inhertited classes:
66 : -------------------
67 : CalStats - This class calculates statistics of new CASA caltables.
68 :
69 : Class public member functions:
70 : ------------------------------
71 : CalStatsReal - This constructor feeds real data to the CalStats base class.
72 : ~CalStatsReal - This destructor deallocates the internal memory of an instance.
73 :
74 : 2011 Dec 11 - Nick Elias, NRAO
75 : Initial version created with public member functions
76 : CalStatsReal() and ~CalStatsReal().
77 :
78 : */
79 :
80 : // -----------------------------------------------------------------------------
81 : // Start of CalStatsReal public member functions
82 : // -----------------------------------------------------------------------------
83 :
84 : /*
85 :
86 : CalStatsReal::CalStatsReal
87 :
88 : Description:
89 : ------------
90 : This class feeds real data to the CalStats base class..
91 :
92 : NB: The FEED axis is always included as an iteration axis by default because one
93 : cannot perform a fit along it. The other iteration axis is defined by the user.
94 :
95 : NB: The default CalStats constructor is called first as the default, then the
96 : standard one is called at the end.
97 :
98 : Inputs:
99 : -------
100 : oValue - This reference to a Cube<Double> instance contains the values.
101 : oValueErr - This reference to a Cube<Double> instance contains the value
102 : errors.
103 : oFlag - This reference to a Cube<Bool> instance contains the flags.
104 : oFeed - This reference to a Vector<String> instance is the feed
105 : abscissae.
106 : oFrequency - This reference to a Vector<Double> instance is the frequency
107 : abscissae.
108 : oTime - This reference to a Vector<Double> instance is the time
109 : abscissae.
110 : eAxisIterUserID - This reference to a CalStats::AXIS enum contains either the
111 : FREQUENCY or TIME iteration axes (user defined).
112 :
113 : Outputs:
114 : --------
115 : None.
116 :
117 : Modification history:
118 : ---------------------
119 : 2011 Dec 11 - Nick Elias, NRAO
120 : Initial version.
121 : 2012 Jan 25 - Nick Elias, NRAO
122 : Error checking added.
123 :
124 : */
125 :
126 : // -----------------------------------------------------------------------------
127 :
128 0 : CalStatsReal::CalStatsReal( const Cube<Double>& oValue,
129 : const Cube<Double>& oValueErr, const Cube<Bool>& oFlag,
130 : const Vector<String>& oFeed, const Vector<Double>& oFrequency,
131 0 : const Vector<Double>& oTime, const CalStats::AXIS& eAxisIterUserID )
132 0 : : CalStats() {
133 :
134 : // Create an instance of the CalStats base class constructor and copy its
135 : // state to this CalStatsReal instance
136 :
137 : CalStats* poCS;
138 :
139 : try {
140 0 : poCS = new CalStats( oValue, oValueErr, oFlag, oFeed, oFrequency, oTime,
141 0 : eAxisIterUserID );
142 : }
143 0 : catch ( AipsError oAE ) {
144 0 : throw( oAE );
145 0 : }
146 :
147 0 : oAxisIterID = IPosition( poCS->axisIterID() );
148 0 : eAxisNonIterID = poCS->axisNonIterID();
149 :
150 0 : oAxisIterFeed = Vector<String>( poCS->axisIterFeed() );
151 0 : oAxisIterUser = Vector<Double>( poCS->axisIterUser() );
152 0 : oAxisNonIter = Vector<Double>( poCS->axisNonIter() );
153 :
154 0 : oStatsShape = IPosition( poCS->statsShape() );
155 :
156 0 : poValue = new Cube<Double>( poCS->value() );
157 0 : poValueErr = new Cube<Double>( poCS->valueErr() );
158 0 : poFlag = new Cube<Bool>( poCS->flag() );
159 :
160 0 : poValueIter = new ArrayIterator<Double>( *poValue, oAxisIterID, false );
161 0 : poValueIter->reset();
162 :
163 0 : poValueErrIter = new ArrayIterator<Double>( *poValueErr, oAxisIterID, false );
164 0 : poValueErrIter->reset();
165 :
166 0 : poFlagIter = new ArrayIterator<Bool>( *poFlag, oAxisIterID, false );
167 0 : poFlagIter->reset();
168 :
169 0 : delete poCS;
170 :
171 :
172 : // Return
173 :
174 0 : return;
175 :
176 0 : }
177 :
178 : // -----------------------------------------------------------------------------
179 :
180 : /*
181 :
182 : CalStatsReal::~CalStatsReal
183 :
184 : Description:
185 : ------------
186 : This destructor deallocates the internal memory of an instance.
187 :
188 : Inputs:
189 : -------
190 : None.
191 :
192 : Outputs:
193 : --------
194 : None.
195 :
196 : Modification history:
197 : ---------------------
198 : 2011 Dec 11 - Nick Elias, NRAO
199 : Initial version.
200 :
201 : */
202 :
203 : // -----------------------------------------------------------------------------
204 :
205 0 : CalStatsReal::~CalStatsReal( void ) {}
206 :
207 : // -----------------------------------------------------------------------------
208 : // End of CalStatsReal public member functions
209 : // -----------------------------------------------------------------------------
210 :
211 : // -----------------------------------------------------------------------------
212 : // End of CalStatsReal class
213 : // -----------------------------------------------------------------------------
214 :
215 : // -----------------------------------------------------------------------------
216 : // Start of CalStatsAmp class
217 : // -----------------------------------------------------------------------------
218 :
219 : /*
220 :
221 : CalStatsAmp
222 :
223 : Description:
224 : ------------
225 : This class converts complex data to amplitudes and initializes the CalStats base
226 : class.
227 :
228 : Inhertited classes:
229 : -------------------
230 : CalStats - This class calculates statistics of new CASA caltables.
231 :
232 : CalStatsAmp public member functions:
233 : ------------------------------------
234 : CalStatsAmp - This generic constructor converts complex data to amplitudes and
235 : initializes the CalStats base class.
236 : ~CalStatsAmp - This destructor deallocates the internal memory of an instance.
237 :
238 : CalStatsAmp public static member functions:
239 : -------------------------------------------
240 : norm - This member function normalizes the amplitudes their errors.
241 :
242 : Modification history:
243 : ---------------------
244 : 2011 Nov 15 - Nick Elias, NRAO
245 : Initial version created with public member functions are
246 : CalStatsAmp() and ~CalStatsAmp(); and public static member
247 : function is norm().
248 :
249 : */
250 :
251 : // -----------------------------------------------------------------------------
252 : // Start of CalStatsAmp public member functions
253 : // -----------------------------------------------------------------------------
254 :
255 : /*
256 :
257 : CalStatsAmp::CalStatsAmp
258 :
259 : Description:
260 : ------------
261 : This class converts complex data and their errors to amplitudes and their errors
262 : and initializes the CalStats base class.
263 :
264 : NB: The FEED axis is always included as an iteration axis by default because one
265 : cannot perform a fit along it. The other iteration axis is defined by the user.
266 :
267 : NB: The default CalStats constructor is called first as the default, then the
268 : standard one is called at the end.
269 :
270 : NB: The are no input data real-imaginary cross correlations available, so the
271 : amplitude errors do not include them.
272 :
273 : Inputs:
274 : -------
275 : oValue - This reference to a Cube<DComplex> instance contains the
276 : values.
277 : oValueErr - This reference to a Cube<Double> instance contains the value
278 : errors (real and imaginary parts).
279 : oFlag - This reference to a Cube<Bool> instance contains the flags.
280 : oFeed - This reference to a Vector<String> instance is the feed
281 : abscissae.
282 : oFrequency - This reference to a Vector<Double> instance is the frequency
283 : abscissae.
284 : oTime - This reference to a Vector<Double> instance is the time
285 : abscissae.
286 : eAxisIterUserID - This reference to a CalStats::AXIS enum contains either the
287 : FREQUENCY or TIME iteration axes (user defined).
288 : bNorm - This reference to a Bool variable contains the normalization
289 : flag (true = normalize, false = don't normalize).
290 :
291 : Outputs:
292 : --------
293 : None.
294 :
295 : Modification history:
296 : ---------------------
297 : 2011 Nov 15 - Nick Elias, NRAO
298 : Initial version.
299 : 2012 Jan 25 - Nick Elias, NRAO
300 : Error checking added.
301 : 2012 Feb 15 - Nick Elias, NRAO
302 : Value error input parameter changed from DComplex to Double.
303 : 2012 Mar 28 - Nick Elias, NRAO
304 : Changed the normalization code so that vectors are iteratively fed
305 : to CalStatsAmp::norm().
306 : 2012 Mar 29 - Nick Elias, NRAO
307 : Member function can now normalize along the frequency and time
308 : axes.
309 :
310 : */
311 :
312 : // -----------------------------------------------------------------------------
313 :
314 0 : CalStatsAmp::CalStatsAmp( const Cube<DComplex>& oValue,
315 : const Cube<Double>& oValueErr, const Cube<Bool>& oFlag,
316 : const Vector<String>& oFeed, const Vector<Double>& oFrequency,
317 : const Vector<Double>& oTime, const CalStats::AXIS& eAxisIterUserID,
318 0 : const Bool& bNorm ) : CalStats() {
319 :
320 : // Calculate the amplitudes and their errors
321 :
322 0 : Cube<Double> oAmp( amplitude( oValue ) );
323 :
324 0 : Cube<Double> oAmpErr = oValueErr.copy();
325 :
326 :
327 : // Create the iterators. The input flag cube is copied since it cannot be
328 : // modified.
329 :
330 0 : IPosition oIterShape( 2, (ssize_t) CalStats::FEED, eAxisIterUserID );
331 :
332 0 : Cube<Bool> oFlagCopy( oFlag.copy() );
333 0 : ArrayIterator<Bool> oFlagIter( oFlagCopy, oIterShape, false );
334 :
335 0 : ArrayIterator<Double> oAmpIter( oAmp, oIterShape, false );
336 0 : ArrayIterator<Double> oAmpErrIter( oAmpErr, oIterShape, false );
337 :
338 :
339 : // If selected, normalize the amplitudes and their errors
340 :
341 0 : while ( bNorm && !oAmpIter.pastEnd() ) {
342 :
343 0 : uInt uiNumAbs = 0;
344 0 : if ( eAxisIterUserID == CalStats::TIME ) {
345 0 : uiNumAbs = oFrequency.nelements();
346 : } else {
347 0 : uiNumAbs = oTime.nelements();
348 : }
349 :
350 0 : if ( uiNumAbs <= 1 ) {
351 0 : oFlagIter.next(); oAmpIter.next(); oAmpErrIter.next();
352 0 : continue;
353 : }
354 :
355 0 : IPosition oShape( 1, uiNumAbs );
356 :
357 0 : Vector<Bool> oFlagV( oFlagIter.array().copy().reform(oShape) );
358 0 : Vector<Double> oAmpV( oAmpIter.array().copy().reform(oShape) );
359 0 : Vector<Double> oAmpErrV( oAmpErrIter.array().copy().reform(oShape) );
360 :
361 0 : norm( oAmpV, oAmpErrV, oFlagV );
362 :
363 0 : oFlagIter.array() = oFlagV;
364 0 : oAmpIter.array() = oAmpV;
365 0 : oAmpErrIter.array() = oAmpErrV;
366 :
367 0 : oFlagIter.next(); oAmpIter.next(); oAmpErrIter.next();
368 :
369 0 : }
370 :
371 :
372 : // Create an instance of the CalStats base class constructor and copy its
373 : // state to this CalStatsAmp instance
374 :
375 : CalStats* poCS;
376 :
377 : try {
378 0 : poCS = new CalStats( oAmp, oAmpErr, oFlagCopy, oFeed, oFrequency, oTime,
379 0 : eAxisIterUserID );
380 : }
381 0 : catch ( AipsError oAE ) {
382 0 : throw( oAE );
383 0 : }
384 :
385 0 : oAxisIterID = IPosition( poCS->axisIterID() );
386 0 : eAxisNonIterID = poCS->axisNonIterID();
387 :
388 0 : oAxisIterFeed = Vector<String>( poCS->axisIterFeed() );
389 0 : oAxisIterUser = Vector<Double>( poCS->axisIterUser() );
390 0 : oAxisNonIter = Vector<Double>( poCS->axisNonIter() );
391 :
392 0 : oStatsShape = IPosition( poCS->statsShape() );
393 :
394 0 : poValue = new Cube<Double>( poCS->value() );
395 0 : poValueErr = new Cube<Double>( poCS->valueErr() );
396 0 : poFlag = new Cube<Bool>( poCS->flag() );
397 :
398 0 : poValueIter = new ArrayIterator<Double>( *poValue, oAxisIterID, false );
399 0 : poValueIter->reset();
400 :
401 0 : poValueErrIter = new ArrayIterator<Double>( *poValueErr, oAxisIterID, false );
402 0 : poValueErrIter->reset();
403 :
404 0 : poFlagIter = new ArrayIterator<Bool>( *poFlag, oAxisIterID, false );
405 0 : poFlagIter->reset();
406 :
407 0 : delete poCS;
408 :
409 :
410 : // Return
411 :
412 0 : return;
413 :
414 0 : }
415 :
416 : // -----------------------------------------------------------------------------
417 :
418 : /*
419 :
420 : CalStatsAmp::~CalStatsAmp
421 :
422 : Description:
423 : ------------
424 : This destructor deallocates the internal memory of an instance.
425 :
426 : Inputs:
427 : -------
428 : None.
429 :
430 : Outputs:
431 : --------
432 : None.
433 :
434 : Modification history:
435 : ---------------------
436 : 2011 Nov 15 - Nick Elias, NRAO
437 : Initial version.
438 :
439 : */
440 :
441 : // -----------------------------------------------------------------------------
442 :
443 0 : CalStatsAmp::~CalStatsAmp( void ) {}
444 :
445 : // -----------------------------------------------------------------------------
446 : // End of CalStatsAmp public member functions
447 : // -----------------------------------------------------------------------------
448 :
449 : // -----------------------------------------------------------------------------
450 : // Start of CalStatsAmp public static member functions
451 : // -----------------------------------------------------------------------------
452 :
453 : /*
454 :
455 : CalStatsAmp::norm
456 :
457 : Description:
458 : ------------
459 : This member function normalizes the amplitudes and their errors.
460 :
461 : NB: The normalization is applied only along the frequency axis.
462 :
463 : NB: The normalization is applied only when the number of unflagged frequencies
464 : is greater than 1.
465 :
466 : NB: All flags corresponding to amplitudes less then 1.0E-08 times the peak
467 : amplitude (along the FREQUENCY axis) are updated to true.
468 :
469 : Inputs:
470 : -------
471 : oAmp - This reference to a Vector<Double> instance contains the unnormalized
472 : amplitudes.
473 : oAmpErr - This reference to a Vector<Double> instance contains the unnormalized
474 : amplitude errors.
475 : oFlag - This reference to a Vector<Bool> instance contains the flags.
476 :
477 : Outputs:
478 : --------
479 : oAmp - This reference to a Vector<Double> instance contains the normalized
480 : amplitudes.
481 : oAmpErr - This reference to a Vector<Double> instance contains the normalized
482 : amplitude errors.
483 : oFlag - This reference to a Vector<Bool> instance contains the flags.
484 :
485 : Modification history:
486 : ---------------------
487 : 2011 Dec 11 - Nick Elias, NRAO
488 : Initial version.
489 : 2012 Jan 25 - Nick Elias, NRAO
490 : Logging capability added.
491 : 2012 Mar 28 - Nick Elias, NRAO
492 : Eliminated iteration and made input amplitude, amplitude error,
493 : and flag cubes into vectors. Flags are now actually used.
494 :
495 : */
496 :
497 : // -----------------------------------------------------------------------------
498 :
499 0 : void CalStatsAmp::norm( Vector<Double>& oAmp, Vector<Double>& oAmpErr,
500 : Vector<Bool>& oFlag ) {
501 :
502 : // Eliminate the flagged amplitudes and their errors
503 :
504 0 : MaskedArray<Double> oAmpM( oAmp, !oFlag );
505 0 : Vector<Double> oAmpC( oAmpM.getCompressedArray() );
506 :
507 0 : MaskedArray<Double> oAmpErrM( oAmpErr, !oFlag );
508 0 : Vector<Double> oAmpErrC( oAmpErrM.getCompressedArray() );
509 :
510 :
511 : // Return if the length of the abscissa is unity
512 :
513 0 : uInt uiNumAbsC = oAmpC.nelements();
514 :
515 0 : if ( uiNumAbsC <= 1 ) {
516 0 : LogIO log( LogOrigin( "CalStatsAmp", "norm", WHERE ) );
517 : // An entire scan for a basline being flagged can cause this
518 : // This is not entirly uncommon so this was switched from WARN -> NORMAL
519 : log << LogIO::NORMAL
520 : << "Abscissa has a dimension <= 1, no normalization"
521 0 : << LogIO::POST;
522 0 : return;
523 0 : }
524 :
525 :
526 : // Normalize the amplitudes and their errors along the abscissa. Flag all low
527 : // amplitudes.
528 :
529 0 : Double dAmpMax = max( oAmpC );
530 :
531 0 : Vector<Bool> oFlagLow( oAmp <= (Double) 1.0E-08*dAmpMax );
532 0 : oFlag = oFlag || oFlagLow;
533 :
534 0 : uInt uiNumAbs = oAmp.nelements();
535 :
536 0 : for ( uInt a=0; a<uiNumAbs; a++ ) {
537 0 : if ( !oFlag[a] ) {
538 0 : oAmp[a] /= dAmpMax;
539 0 : oAmpErr[a] /= dAmpMax;
540 : }
541 : }
542 :
543 :
544 : // Return
545 :
546 0 : return;
547 :
548 0 : }
549 :
550 : // -----------------------------------------------------------------------------
551 : // End of CalStatsAmp public static member functions
552 : // -----------------------------------------------------------------------------
553 :
554 : // -----------------------------------------------------------------------------
555 : // End of CalStatsAmp class
556 : // -----------------------------------------------------------------------------
557 :
558 : // -----------------------------------------------------------------------------
559 : // Start of CalStatsPhase class
560 : // -----------------------------------------------------------------------------
561 :
562 : /*
563 :
564 : CalStatsPhase
565 :
566 : Description:
567 : ------------
568 : This class converts complex data to phases and initializes the CalStats base
569 : class.
570 :
571 : Inhertited classes:
572 : -------------------
573 : CalStats - This class calculates statistics of new CASA caltables.
574 :
575 : Class public member functions:
576 : ------------------------------
577 : CalStatsPhase - This generic constructor converts complex data to amplitudes
578 : and initializes the CalStats base class.
579 : ~CalStatsPhase - This destructor deallocates the internal memory of an instance.
580 :
581 : CalStatsPhase public static member functions:
582 : ---------------------------------------------
583 : unwrapGD - This member function unwraps the phases along the frequency axis
584 : with respect to the group delay.
585 : unwrapSimple - This member function performs a simple unwrapping procedure for
586 : both frequency and temporal abscissae.
587 :
588 : CalStatsPhase private static member functions:
589 : ----------------------------------------------
590 : fringePacket2 - This member function forms the squared-amplitude fringe packet.
591 :
592 : CalStatsPhase templated private static member functions:
593 : --------------------------------------------------------
594 : maxLocation - This member function finds the abscissa corresponding to the peak
595 : value of an ordinate vector.
596 :
597 : Modification history:
598 : ---------------------
599 : 2011 Nov 15 - Nick Elias, NRAO
600 : Initial version created with public member functions are
601 : CalStatsPhase() and ~CalStatsPhase(); and public static member
602 : function is unwrap().
603 : 2012 Mar 27 - Nick Elias, NRAO
604 : Private static member functions fringePacket2() and maxLocation()
605 : added. Private static member variables NUM_ITER_UNWRAP and
606 : NEW_RANGE_FACTOR added.
607 : 2012 Mar 30 - Nick Elias, NRAO
608 : Public static member function unwrap() renamed to unwrapGD().
609 : Public static member function unwrapSimple() added.
610 :
611 : */
612 :
613 : // -----------------------------------------------------------------------------
614 : // Start of CalStatsPhase public member functions
615 : // -----------------------------------------------------------------------------
616 :
617 : /*
618 :
619 : CalStatsPhase::CalStatsPhase
620 :
621 : Description:
622 : ------------
623 : This class converts complex data and their errors to phases and their errors and
624 : initializes the CalStats base class.
625 :
626 : NB: The FEED axis is always included as an iteration axis by default because one
627 : cannot perform a fit along it. The other iteration axis is defined by the user.
628 :
629 : NB: All flags corresponding to amplitudes less then 1.0E-08 times the peak
630 : amplitude are updated to true.
631 :
632 : NB: The default CalStats constructor is called first as the default, then the
633 : standard one is called at the end.
634 :
635 : NB: The are no input data real-imaginary cross correlations available, so the
636 : phase errors do not include them.
637 :
638 : NB: For unwrapping along the time axis (iteration axis CalStats::FREQUENCY), the
639 : unwrapSimple() function is used. For unwrapping along the frequency axis
640 : (iteration axis CalStats::TIME), dJumpMax==0.0 leads to unwrapGD() while
641 : dJumpMax!=0.0 leads to unwrapSimple().
642 :
643 : Inputs:
644 : -------
645 : oValue - This reference to a Cube<DComplex> instance contains the
646 : values.
647 : oValueErr - This reference to a Cube<Double> instance contains the value
648 : errors (real and imaginary parts).
649 : oFlag - This reference to a Cube<Bool> instance contains the flags.
650 : oFeed - This reference to a Vector<String> instance is the feed
651 : abscissae.
652 : oFrequency - This reference to a Vector<Double> instance is the frequency
653 : abscissae.
654 : oTime - This reference to a Vector<Double> instance is the time
655 : abscissae.
656 : eAxisIterUserID - This reference to a CalStats::AXIS enum contains either the
657 : CalStats::FREQUENCY or CalStats::TIME iteration axes (user
658 : defined).
659 : bUnwrap - This reference to a Bool variable contains the unwrapping flag
660 : (true = unwrap, false = don't unwrap).
661 : dJumpMax - This reference to a Double variable contains the maximum
662 : deviation from +/- M_PI for adjacent points to be unwrapped
663 : by +/- 2.0*M_PI (in radians). This parameter is always used
664 : when the specified iteration axis is CalStats::FREQUENCY
665 : (unwrapping along the CalStats::TIME axis). If the specified
666 : iteration axis is CalStats::TIME (unwrapping along the
667 : CalStats::FREQUENCY axis), this parameter selects the type of
668 : unwrapping (dJumpMax==0.0 --> group-delay unwrapping, dJumpMax
669 : != 0.0 --> simple unwrapping).
670 :
671 : Outputs:
672 : --------
673 : None.
674 :
675 : Modification history:
676 : ---------------------
677 : 2011 Nov 15 - Nick Elias, NRAO
678 : Initial version.
679 : 2011 Dec 11 - Nick Elias, NRAO
680 : Flag updates according to low amplitude was added.
681 : 2012 Jan 25 - Nick Elias, NRAO
682 : Error checking added.
683 : 2012 Feb 15 - Nick Elias, NRAO
684 : Value error input parameter changed from DComplex to Double.
685 : 2012 Mar 28 - Nick Elias, NRAO
686 : Changed the unwrapping code so that vectors are iteratively fed to
687 : CalStatsPhase::unwrapSimple() or CalStatsPhase::unwrapGD().
688 : 2012 Apr 01 - Nick Elias, NRAO
689 : Input parameter dJumpMax added.
690 :
691 : */
692 :
693 : // -----------------------------------------------------------------------------
694 :
695 0 : CalStatsPhase::CalStatsPhase( const Cube<DComplex>& oValue,
696 : const Cube<Double>& oValueErr, const Cube<Bool>& oFlag,
697 : const Vector<String>& oFeed, const Vector<Double>& oFrequency,
698 : const Vector<Double>& oTime, const CalStats::AXIS& eAxisIterUserID,
699 0 : const Bool& bUnwrap, const Double& dJumpMax ) : CalStats() {
700 :
701 : // Calculate the phases and the initial phase error cube (set to 0.0)
702 :
703 0 : Cube<Double> oPhase( phase(oValue) );
704 :
705 0 : Cube<Double> oPhaseErr( oPhase.shape() );
706 :
707 :
708 : // Create the iterators. The input flag cube is copied since it cannot be
709 : // modified.
710 :
711 0 : IPosition oIterShape( 2, (ssize_t) CalStats::FEED, eAxisIterUserID );
712 :
713 0 : ReadOnlyArrayIterator<DComplex> oValueIter( oValue, oIterShape, false );
714 0 : ReadOnlyArrayIterator<Double> oValueErrIter( oValueErr, oIterShape, false );
715 :
716 0 : Cube<Bool> oFlagCopy( oFlag.copy() );
717 0 : ArrayIterator<Bool> oFlagIter( oFlagCopy, oIterShape, false );
718 :
719 0 : ArrayIterator<Double> oPhaseIter( oPhase, oIterShape, false );
720 0 : ArrayIterator<Double> oPhaseErrIter( oPhaseErr, oIterShape, false );
721 :
722 :
723 : // If selected, unwrap the phases
724 :
725 0 : while ( bUnwrap && !oPhaseIter.pastEnd() ) {
726 :
727 0 : uInt uiNumAbs = 0;
728 0 : if ( eAxisIterUserID == CalStats::TIME ) {
729 0 : uiNumAbs = oFrequency.nelements();
730 : } else {
731 0 : uiNumAbs = oTime.nelements();
732 : }
733 :
734 0 : if ( uiNumAbs <= 1 ) {
735 0 : oFlagIter.next(); oPhaseIter.next();
736 0 : continue;
737 : }
738 :
739 0 : IPosition oShape( 1, uiNumAbs );
740 :
741 0 : Vector<Bool> oFlagV( oFlagIter.array().copy().reform(oShape) );
742 0 : Vector<Double> oPhaseV( oPhaseIter.array().copy().reform(oShape) );
743 :
744 0 : if ( eAxisIterUserID == CalStats::TIME ) {
745 0 : if ( dJumpMax == 0.0 ) {
746 0 : unwrapGD( oPhaseV, oFrequency, oFlagV );
747 : } else {
748 0 : unwrapSimple( oPhaseV, dJumpMax, oFlagV );
749 : }
750 : } else {
751 0 : unwrapSimple( oPhaseV, dJumpMax, oFlagV );
752 : }
753 :
754 0 : oPhaseIter.array() = oPhaseV;
755 :
756 0 : oFlagIter.next(); oPhaseIter.next();
757 :
758 0 : }
759 :
760 :
761 : // Reset the iterators
762 :
763 0 : oFlagIter.reset(); oPhaseIter.next();
764 :
765 :
766 : // Calculate the phase errors. They require dividing by amplitude. Set the
767 : // flags according to low amplitudes (corresponding to phase errors greater
768 : // than M_PI). All flagged phase errors are set to M_PI. The iteration axes
769 : // are FEED and TIME. A FREQUENCY for loop is located inside the FEED/TIME
770 : // iteration loop to set the errors.
771 :
772 0 : while ( !oValueIter.pastEnd() ) {
773 :
774 0 : uInt uiNumAbs = oValueIter.array().nelements();
775 0 : IPosition oShape( 1, uiNumAbs );
776 :
777 0 : Vector<Double> oPhaseErrV( oShape );
778 :
779 0 : Vector<DComplex> oValueV( oValueIter.array().copy().reform(oShape) );
780 0 : Vector<Double> oValueErrV( oValueErrIter.array().copy().reform(oShape) );
781 :
782 0 : Vector<Double> oAmpV( amplitude(oValueV) );
783 :
784 0 : Vector<Bool> oFlagV( oFlagIter.array().copy().reform(oShape) );
785 0 : Vector<Bool> oFlagLow( oAmpV <= oValueErrV/((Double) M_PI) );
786 :
787 0 : oFlagV = oFlagV || oFlagLow;
788 0 : oFlagIter.array() = oFlagV;
789 :
790 0 : for ( uInt a=0; a<uiNumAbs; a++ ) {
791 0 : if ( !oFlagV[a] ) {
792 0 : oPhaseErrV[a] = oValueErrV[a] / oAmpV[a];
793 : } else {
794 0 : oPhaseErrV[a] = (Double) M_PI;
795 : }
796 : }
797 :
798 0 : oPhaseErrIter.array() = oPhaseErrV;
799 :
800 0 : oValueIter.next(); oValueErrIter.next(); oFlagIter.next();
801 0 : oPhaseErrIter.next();
802 :
803 0 : }
804 :
805 :
806 : // Create an instance of the CalStats base class constructor and copy its
807 : // state to this CalStatsPhase instance
808 :
809 : CalStats* poCS;
810 :
811 : try {
812 0 : poCS = new CalStats( oPhase, oPhaseErr, oFlagCopy, oFeed, oFrequency, oTime,
813 0 : eAxisIterUserID );
814 : }
815 0 : catch ( AipsError oAE ) {
816 0 : throw( oAE );
817 0 : }
818 :
819 0 : oAxisIterID = IPosition( poCS->axisIterID() );
820 0 : eAxisNonIterID = poCS->axisNonIterID();
821 :
822 0 : oAxisIterFeed = Vector<String>( poCS->axisIterFeed() );
823 0 : oAxisIterUser = Vector<Double>( poCS->axisIterUser() );
824 0 : oAxisNonIter = Vector<Double>( poCS->axisNonIter() );
825 :
826 0 : oStatsShape = IPosition( poCS->statsShape() );
827 :
828 0 : poValue = new Cube<Double>( poCS->value() );
829 0 : poValueErr = new Cube<Double>( poCS->valueErr() );
830 0 : poFlag = new Cube<Bool>( poCS->flag() );
831 :
832 0 : poValueIter = new ArrayIterator<Double>( *poValue, oAxisIterID, false );
833 0 : poValueIter->reset();
834 :
835 0 : poValueErrIter = new ArrayIterator<Double>( *poValueErr, oAxisIterID, false );
836 0 : poValueErrIter->reset();
837 :
838 0 : poFlagIter = new ArrayIterator<Bool>( *poFlag, oAxisIterID, false );
839 0 : poFlagIter->reset();
840 :
841 0 : delete poCS;
842 :
843 :
844 : // Return
845 :
846 0 : return;
847 :
848 0 : }
849 :
850 : // -----------------------------------------------------------------------------
851 :
852 : /*
853 :
854 : CalStatsPhase::~CalStatsPhase
855 :
856 : Description:
857 : ------------
858 : This destructor deallocates the internal memory of an instance.
859 :
860 : Inputs:
861 : -------
862 : None.
863 :
864 : Outputs:
865 : --------
866 : None.
867 :
868 : Modification history:
869 : ---------------------
870 : 2011 Nov 15 - Nick Elias, NRAO
871 : Initial version.
872 :
873 : */
874 :
875 : // -----------------------------------------------------------------------------
876 :
877 0 : CalStatsPhase::~CalStatsPhase( void ) {}
878 :
879 : // -----------------------------------------------------------------------------
880 : // End of CalStatsPhase public member functions
881 : // -----------------------------------------------------------------------------
882 :
883 : // -----------------------------------------------------------------------------
884 : // Start of CalStatsPhase public static member functions
885 : // -----------------------------------------------------------------------------
886 :
887 : /*
888 :
889 : CalStatsPhase::unwrapSimple
890 :
891 : Description:
892 : ------------
893 : This member function performs a simple unwrapping procedure for both frequency
894 : and temporal abscissae.
895 :
896 : Algorithm:
897 : ----------
898 : * If the first point is a little less than +M_PI and the second point is a
899 : little more than -M_PI, add 2.0*M_PI to the second point.
900 : * If the first point is a little more than -M_PI and the second point is a
901 : ittle less than +M_PI, subtract 2.0*M_PI from the second point.
902 : * "A little more" means within dJumpMax of +/- M_PI.
903 : * Flagged data are ignored, which could be a problem if a lot of them occur
904 : sequentially.
905 :
906 : Inputs:
907 : -------
908 : oPhase - This reference to a Vector<Double> instance contains the wrapped
909 : phases.
910 : dJumpMax - This reference to a Double variable contains the maximum deviation
911 : from +/- M_PI for adjacent points to be unwrapped by +/- 2.0*M_PI (in
912 : radians).
913 : oFlag - This reference to a Vector<Bool> instance contains the flags.
914 :
915 : Outputs:
916 : --------
917 : oPhase - This reference to a Vector<Double> instance contains the unwrapped
918 : phases.
919 :
920 : Modification history:
921 : ---------------------
922 : 2012 Mar 30 - Nick Elias, NRAO
923 : Initial version.
924 :
925 : */
926 :
927 : // -----------------------------------------------------------------------------
928 :
929 0 : void CalStatsPhase::unwrapSimple( Vector<Double>& oPhase,
930 : const Double& dJumpMax, const Vector<Bool>& oFlag ) {
931 :
932 : // Initialize the number of elements and make an immutable copy of the input
933 : // phase vector
934 :
935 0 : uInt uiNumAbs = oPhase.nelements();
936 :
937 0 : Vector<Double> oPhaseIn( oPhase.copy() );
938 :
939 :
940 : // Perform the simple unwrapping. The unwrapping occurs if subsequent points
941 : // are both within dJumpMax of M_PI or -M_PI. Flagged data are ignored.
942 :
943 0 : for ( uInt a=1; a<uiNumAbs; a++ ) {
944 :
945 0 : if ( oFlag[a-1] ) continue;
946 :
947 0 : uInt a2 = 0;
948 0 : for ( a2=a; a2<uiNumAbs; a2++ ) {
949 0 : if ( !oFlag[a2] ) break;
950 : }
951 0 : if ( a2 >= uiNumAbs ) return;
952 :
953 : // The first point is a little less than +M_PI and the second point is a
954 : // little more than -M_PI, add 2.0*M_PI to the second point
955 0 : if ( M_PI-oPhaseIn[a-1] <= dJumpMax && M_PI+oPhaseIn[a2] <= dJumpMax ) {
956 0 : for ( uInt a3=a2; a3<uiNumAbs; a3++ ) oPhase[a3] += 2.0 * M_PI;
957 : }
958 :
959 : // The first point is a little more than -M_PI and the second point is a
960 : // little less than +M_PI, subtract 2.0*M_PI from the second point
961 0 : if ( M_PI+oPhaseIn[a-1] <= dJumpMax && M_PI-oPhaseIn[a2] <= dJumpMax ) {
962 0 : for ( uInt a3=a2; a3<uiNumAbs; a3++ ) oPhase[a3] -= 2.0 * M_PI;
963 : }
964 :
965 : }
966 :
967 :
968 : // Return
969 :
970 0 : return;
971 :
972 0 : }
973 :
974 : // -----------------------------------------------------------------------------
975 :
976 : /*
977 :
978 : CalStatsPhase::unwrapGD
979 :
980 : Description:
981 : ------------
982 : This member function unwraps the phases along the frequency axis with respect to
983 : the group delay.
984 :
985 : NB: The unwrapping is applied only along the frequency axis. The unwrapping is
986 : applied only when the number of unflagged frequencies is greater than 1 (setting
987 : the bar very low). If the number of unflagged frequencies is small, you may
988 : get crappy results.
989 :
990 : NB: This algorithm should work better with low-S/N data compared to simple
991 : unwrapping, but it will not work if the S/N is abysmally low. To understand
992 : why this statement is true, consider an analogy with BPOLY. BPOLY is used when
993 : the S/N per channel is low. The fit uses all bandpass data at once to maximize
994 : the available S/N. This phase unwrapping algorithm creates a fringe packet,
995 : using all of the bandpass data at once, to estimate the delay before unwrapping.
996 :
997 : NB: If the true delay is outside of the unaliased search window, the unwrapping
998 : aliases into the search window. This scenario is OK if only continuous phase
999 : plotting is required.
1000 :
1001 : NB: If the dispersion or other wiggles are large, this technique may not work
1002 : well.
1003 :
1004 : Algorithm:
1005 : ----------
1006 : * Sort the input frequencies and form the differences. Sort the frequency
1007 : differences. The desired frequency increment is the minimum value (except for
1008 : zero, which arises datasets from spectral windows with the same frequencies;
1009 : choose the next largest value).
1010 : * Determine the maximum frequency.
1011 : * For the first iteration:
1012 : - Calculate the initial time increment using the maximum frequency (Nyquist
1013 : criterion).
1014 : - Calculate the initial time range using the frequency increment.
1015 : - Form the initial time vector. The edges of the unaliased search window are
1016 : +/- 0.5 times the initial time range (centered around zero group delay).
1017 : - Form the initial squared-amplitude fringe packet.
1018 : - Find the peak of the initial squared-amplitude fringe packet. The time
1019 : corresponding to this peak is the initial group delay estimate.
1020 : * For subsequent iterations:
1021 : - Set the new time range to NEW_RANGE_FACTOR times the previous time
1022 : increment.
1023 : - Set the new time increment to the new time range divided by the number of
1024 : times (this number is the same as the initial iteration and doesn't change).
1025 : - Form the new time vector. The edges of the search window are +/- times the
1026 : new time range (centered around the previous group delay estimate).
1027 : - Form the new squared-amplitude fringe packet.
1028 : - Find the peak of the new squared-amplitude fringe packet. The time
1029 : corresponding to this peak is the new group delay estimate.
1030 : - Has the group delay converged? If not repeat for another iteration,
1031 : otherwise stop iterating. The maximum number of iterations is
1032 : NUM_ITER_UNWRAP.
1033 : * Form the phase versus frequency using the group delay.
1034 : * Subtract the group delay phases from the input phases to form the differential
1035 : delay.
1036 : * Subtract the first differential delay from all of the differential delays.
1037 : * Starting at the first frequency, search for phase jumps. For each one found,
1038 : add/subtract 2.0*M_PI*N (N=the number of phase jumps; it can be positive or
1039 : negative) for all subsequenct frequencies. Flagged data are not unwrapped.
1040 :
1041 : Inputs:
1042 : -------
1043 : oPhase - This reference to a Vector<Double> instance contains the wrapped
1044 : phases.
1045 : oFrequency - This reference to a Vector<Double> instance is the frequency
1046 : abscissa.
1047 : oFlag - This reference to a Vector<Bool> instance contains the flags.
1048 :
1049 : Outputs:
1050 : --------
1051 : oPhase - This reference to a Vector<Double> instance contains the unwrapped
1052 : phases.
1053 :
1054 : Modification history:
1055 : ---------------------
1056 : 2011 Nov 15 - Nick Elias, NRAO
1057 : Initial stub version.
1058 : 2012 Jan 25 - Nick Elias, NRAO
1059 : Logging capability added.
1060 : 2012 Mar 27 - Nick Elias, NRAO
1061 : Initial working version.
1062 : 2012 Mar 28 - Nick Elias, NRAO
1063 : Eliminated iteration and made input phase and flag cubes into
1064 : vectors. Flags are now actually used.
1065 :
1066 : */
1067 :
1068 : // -----------------------------------------------------------------------------
1069 :
1070 0 : void CalStatsPhase::unwrapGD( Vector<Double>& oPhase,
1071 : const Vector<Double>& oFrequency, const Vector<Bool>& oFlag ) {
1072 :
1073 : // Eliminate the flagged phases and frequencies
1074 :
1075 0 : MaskedArray<Double> oPhaseM( oPhase, !oFlag );
1076 0 : Vector<Double> oPhaseC( oPhaseM.getCompressedArray() );
1077 :
1078 0 : MaskedArray<Double> oFrequencyM( oFrequency, !oFlag );
1079 0 : Vector<Double> oFrequencyC( oFrequencyM.getCompressedArray() );
1080 :
1081 :
1082 : // Return if the length of the frequency axis is unity
1083 :
1084 0 : uInt uiNumFrequencyC = oFrequencyC.nelements();
1085 :
1086 0 : if ( uiNumFrequencyC <= 1 ) {
1087 0 : LogIO log( LogOrigin( "CalStatsPhase", "unwrap", WHERE ) );
1088 : log << LogIO::WARN
1089 : << "Frequency axis has a dimension <= 1, no unwrapping performed"
1090 0 : << LogIO::POST;
1091 0 : return;
1092 0 : }
1093 :
1094 :
1095 : // Calculate the minimum frequency increment. If there are overlapping
1096 : // spectral windows, some of the frequency increments may be zero if the
1097 : // duplicate channels are not flagged. Make sure to ignore all zeros.
1098 :
1099 0 : Vector<Double> oFreqDeltaC( uiNumFrequencyC-1 );
1100 :
1101 0 : for ( uInt f=1; f<uiNumFrequencyC; f++ ) {
1102 0 : oFreqDeltaC[f-1] = oFrequencyC[f] - oFrequencyC[f-1];
1103 : }
1104 :
1105 0 : Sort::Order eOrder = Sort::Ascending;
1106 0 : Int iOptions = Sort::QuickSort;
1107 0 : GenSort<Double>::sort( oFreqDeltaC, eOrder, (int) iOptions );
1108 :
1109 0 : uInt fD = 0;
1110 0 : for ( fD=0; fD<uiNumFrequencyC-1; fD++ ) {
1111 0 : if ( oFreqDeltaC[fD] != 0.0 ) break;
1112 : }
1113 :
1114 0 : if ( fD >= uiNumFrequencyC-1 ) {
1115 0 : LogIO log( LogOrigin( "CalStatsPhase", "unwrap", WHERE ) );
1116 0 : log << "Something is very wrong with the frequencies" << LogIO::EXCEPTION;
1117 0 : }
1118 :
1119 0 : Double dFreqDeltaC = oFreqDeltaC[fD];
1120 :
1121 :
1122 : // Calculate the maximum frequency
1123 :
1124 0 : Double dFreqMaxC = max( oFrequencyC );
1125 :
1126 :
1127 : // Calculate the initial time increment and range
1128 :
1129 0 : Double dTimeDelta = 0.5 / dFreqMaxC;
1130 :
1131 0 : Double dTimeRange = 1.0 / dFreqDeltaC;
1132 :
1133 :
1134 : // Initialize the time vector. The center of the time range is zero group
1135 : // delay.
1136 :
1137 0 : uInt uiNumTime = (uInt) fabs( dTimeRange/dTimeDelta + 1.0 );
1138 :
1139 0 : Vector<Double> oTime( uiNumTime );
1140 0 : indgen<Double>( oTime, 0.0-0.5*dTimeRange, dTimeDelta );
1141 :
1142 :
1143 : // Form the squared magnitude of the initial fringe packet and find the group
1144 : // delay
1145 :
1146 0 : Vector<Double> oFringePacket( fringePacket2( oPhaseC, oFrequencyC, oTime ) );
1147 :
1148 0 : Double dGroupDelay = maxLocation<Double>( oTime, oFringePacket );
1149 :
1150 :
1151 : // Initialize the number of iterations
1152 :
1153 0 : uInt uiNumIter = 0;
1154 :
1155 :
1156 : // Iteratively improve the group delay estimate
1157 :
1158 0 : while ( uiNumIter < NUM_ITER_UNWRAP ) {
1159 :
1160 0 : dTimeRange = NEW_RANGE_FACTOR * dTimeDelta;
1161 0 : Double dTimeDeltaNew = dTimeRange / (uiNumTime+1);
1162 :
1163 0 : indgen<Double>( oTime, dGroupDelay-0.5*dTimeDelta, dTimeDeltaNew );
1164 :
1165 0 : oFringePacket = fringePacket2( oPhaseC, oFrequencyC, oTime );
1166 0 : Double dGroupDelayNew = maxLocation<Double>( oTime, oFringePacket );
1167 :
1168 0 : if ( fabs(dGroupDelayNew-dGroupDelay)/dTimeDelta < 1.0E-08 ) break;
1169 :
1170 0 : dTimeDelta = dTimeDeltaNew;
1171 0 : dGroupDelay = dGroupDelayNew;
1172 :
1173 0 : uiNumIter++;
1174 :
1175 : }
1176 :
1177 :
1178 : // If the number of iterations has been exceeded, print a warning and set
1179 : // the group delay to the middle of the time range
1180 :
1181 0 : if ( uiNumIter > NUM_ITER_UNWRAP ) {
1182 :
1183 0 : dGroupDelay = mean( oTime );
1184 :
1185 0 : LogIO log( LogOrigin( "CalStatsPhase", "unwrap", WHERE ) );
1186 : log << LogIO::WARN
1187 : << "Number of iterations exceeded for group delay calculation\n"
1188 : << "Using the mean time from the last iteration"
1189 0 : << LogIO::POST;
1190 :
1191 0 : }
1192 :
1193 :
1194 : // Calculate the differential phases, which are the input phases minus the
1195 : // group delay. Make sure to subtract the initial value from all values.
1196 :
1197 0 : uInt uiNumFrequency = oFrequency.nelements();
1198 :
1199 0 : Vector<Double> oDiffPhase( oPhase - 2.0*M_PI*oFrequency*dGroupDelay );
1200 0 : oDiffPhase -= oDiffPhase[0];
1201 :
1202 :
1203 : // If phase wraps are found, unwrap them. Flagged data are not unwrapped.
1204 :
1205 0 : for ( uInt f=1; f<uiNumFrequency; f++ ) {
1206 :
1207 0 : if ( oFlag[f-1] ) continue;
1208 :
1209 0 : uInt f2 = 0;
1210 0 : for ( f2=f; f2<uiNumFrequency; f2++ ) if ( !oFlag[f2] ) break;
1211 :
1212 0 : if ( f2 >= uiNumFrequency ) break;
1213 :
1214 0 : Double dDiffPhase = oDiffPhase[f2] - oDiffPhase[f-1];
1215 :
1216 0 : Int N = 0;
1217 0 : while ( dDiffPhase + 2.0*M_PI*N <= -M_PI ) N++;
1218 0 : while ( dDiffPhase + 2.0*M_PI*N > M_PI ) N--;
1219 :
1220 0 : if ( N == 0 ) continue;
1221 :
1222 0 : for ( uInt f3=f2; f3<uiNumFrequency; f3++ ) oPhase[f3] += 2.0 * M_PI * N;
1223 :
1224 : }
1225 :
1226 :
1227 : // Return
1228 :
1229 0 : return;
1230 :
1231 0 : }
1232 :
1233 : // -----------------------------------------------------------------------------
1234 : // End of CalStatsPhase public static member functions
1235 : // -----------------------------------------------------------------------------
1236 :
1237 : // -----------------------------------------------------------------------------
1238 : // Start of CalStatsPhase private static member functions
1239 : // -----------------------------------------------------------------------------
1240 :
1241 : /*
1242 :
1243 : Description:
1244 : ------------
1245 : This member function forms the squared-amplitude fringe packet.
1246 :
1247 : Inputs:
1248 : -------
1249 : oPhase - This Vector<Double>() instance contains the wrapped phases.
1250 : oFrequency - This Vector<Double>() instance contains the frequencies.
1251 : oTime - This Vector<Double>() instance contains the times.
1252 :
1253 : Outputs:
1254 : --------
1255 : The reference to the Vector<Double>() instance containing the squared-amplitude
1256 : fringe packet, returned via the function value.
1257 :
1258 : Modification history:
1259 : ---------------------
1260 : 2012 Mar 27 - Nick Elias, NRAO
1261 : Initial version.
1262 :
1263 : */
1264 :
1265 : // -----------------------------------------------------------------------------
1266 :
1267 0 : Vector<Double>& CalStatsPhase::fringePacket2( const Vector<Double>& oPhase,
1268 : const Vector<Double>& oFrequency, const Vector<Double>& oTime ) {
1269 :
1270 : // Initialize the pointer to the squared-amplitude fringe packet
1271 :
1272 0 : uInt uiNumTime = oTime.nelements();
1273 :
1274 : Vector<Double>* poFringePacket;
1275 0 : poFringePacket = new Vector<Double>( uiNumTime, 0.0 );
1276 :
1277 :
1278 : // Calculate the complex fringe packet for all times using all frequencies and
1279 : // phases
1280 :
1281 0 : uInt uiNumFrequency = oFrequency.nelements();
1282 :
1283 0 : Vector<DComplex> oFringePacketC( uiNumTime, DComplex( 0.0, 0.0 ) );
1284 :
1285 0 : for ( uInt f=0; f<uiNumFrequency; f++ ) {
1286 0 : Vector<DComplex> oFringePacketCF( uiNumTime, DComplex( 0.0, 0.0 ) );
1287 0 : setReal( oFringePacketCF, cos( 2.0*M_PI*oFrequency[f]*oTime - oPhase[f] ) );
1288 0 : setImag( oFringePacketCF, sin( 2.0*M_PI*oFrequency[f]*oTime - oPhase[f] ) );
1289 0 : operator+=( oFringePacketC, oFringePacketCF );
1290 0 : }
1291 :
1292 :
1293 : // Calculate and return the reference to the squared-amplitude fringe packet
1294 :
1295 0 : *poFringePacket = square( amplitude( oFringePacketC ) );
1296 0 : *poFringePacket /= (Double) uiNumFrequency * uiNumFrequency;
1297 :
1298 0 : return( *poFringePacket );
1299 :
1300 0 : }
1301 :
1302 : // -----------------------------------------------------------------------------
1303 : // End of CalStatsPhase private static member functions
1304 : // -----------------------------------------------------------------------------
1305 :
1306 : // -----------------------------------------------------------------------------
1307 : // End of CalStatsPhase class
1308 : // -----------------------------------------------------------------------------
1309 :
1310 : };
1311 :
1312 : // -----------------------------------------------------------------------------
1313 : // End of casa namespace
1314 : // -----------------------------------------------------------------------------
|
