source: trunk/src/SDMath.cc@ 511

Last change on this file since 511 was 502, checked in by kil064, 20 years ago

add function rotateLinPolPhase
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 56.1 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDMath.cc: A collection of single dish mathematical operations
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: Malte.Marquarding@csiro.au
23//# Postal address: Malte Marquarding,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
29//# $Id:
30//#---------------------------------------------------------------------------
31#include <vector>
32
33#include <casa/aips.h>
34#include <casa/iostream.h>
35#include <casa/iomanip.h>
36#include <casa/BasicSL/String.h>
37#include <casa/Arrays/IPosition.h>
38#include <casa/Arrays/Array.h>
39#include <casa/Arrays/ArrayIter.h>
40#include <casa/Arrays/VectorIter.h>
41#include <casa/Arrays/ArrayMath.h>
42#include <casa/Arrays/ArrayLogical.h>
43#include <casa/Arrays/MaskedArray.h>
44#include <casa/Arrays/MaskArrMath.h>
45#include <casa/Arrays/MaskArrLogi.h>
46#include <casa/Arrays/Matrix.h>
47#include <casa/BasicMath/Math.h>
48#include <casa/Containers/Block.h>
49#include <casa/Exceptions.h>
50#include <casa/Quanta/Quantum.h>
51#include <casa/Quanta/Unit.h>
52#include <casa/Quanta/MVEpoch.h>
53#include <casa/Quanta/MVTime.h>
54#include <casa/Utilities/Assert.h>
55
56#include <coordinates/Coordinates/SpectralCoordinate.h>
57#include <coordinates/Coordinates/CoordinateSystem.h>
58#include <coordinates/Coordinates/CoordinateUtil.h>
59#include <coordinates/Coordinates/FrequencyAligner.h>
60
61#include <lattices/Lattices/LatticeUtilities.h>
62#include <lattices/Lattices/RebinLattice.h>
63
64#include <measures/Measures/MEpoch.h>
65#include <measures/Measures/MDirection.h>
66#include <measures/Measures/MPosition.h>
67
68#include <scimath/Mathematics/VectorKernel.h>
69#include <scimath/Mathematics/Convolver.h>
70#include <scimath/Mathematics/InterpolateArray1D.h>
71#include <scimath/Functionals/Polynomial.h>
72
73#include <tables/Tables/Table.h>
74#include <tables/Tables/ScalarColumn.h>
75#include <tables/Tables/ArrayColumn.h>
76#include <tables/Tables/ReadAsciiTable.h>
77
78#include "MathUtils.h"
79#include "SDDefs.h"
80#include "SDAttr.h"
81#include "SDContainer.h"
82#include "SDMemTable.h"
83
84#include "SDMath.h"
85#include "SDPol.h"
86
87using namespace casa;
88using namespace asap;
89
90
91SDMath::SDMath()
92{;}
93
94SDMath::SDMath(const SDMath& other)
95{
96
97// No state
98
99}
100
101SDMath& SDMath::operator=(const SDMath& other)
102{
103 if (this != &other) {
104// No state
105 }
106 return *this;
107}
108
109SDMath::~SDMath()
110{;}
111
112
113
114SDMemTable* SDMath::frequencyAlignment(const SDMemTable& in,
115 const String& refTime,
116 const String& method,
117 Bool perFreqID) const
118{
119// Get frame info from Table
120
121 std::vector<std::string> info = in.getCoordInfo();
122
123// Parse frequency system
124
125 String systemStr(info[1]);
126 String baseSystemStr(info[3]);
127 if (baseSystemStr==systemStr) {
128 throw(AipsError("You have not set a frequency frame different from the initial - use function set_freqframe"));
129 }
130//
131 MFrequency::Types freqSystem;
132 MFrequency::getType(freqSystem, systemStr);
133
134// Do it
135
136 return frequencyAlign(in, freqSystem, refTime, method, perFreqID);
137}
138
139
140
141CountedPtr<SDMemTable> SDMath::average(const Block<CountedPtr<SDMemTable> >& in,
142 const Vector<Bool>& mask, Bool scanAv,
143 const String& weightStr, Bool alignFreq) const
144//
145// Weighted averaging of spectra from one or more Tables.
146//
147{
148
149// Convert weight type
150
151 WeightType wtType = NONE;
152 convertWeightString(wtType, weightStr);
153
154// Create output Table by cloning from the first table
155
156 SDMemTable* pTabOut = new SDMemTable(*in[0],True);
157 if (in.nelements() > 1) {
158 for (uInt i=1; i < in.nelements(); ++i) {
159 pTabOut->appendToHistoryTable(in[i]->getHistoryTable());
160 }
161 }
162// Setup
163
164 IPosition shp = in[0]->rowAsMaskedArray(0).shape(); // Must not change
165 Array<Float> arr(shp);
166 Array<Bool> barr(shp);
167 const Bool useMask = (mask.nelements() == shp(asap::ChanAxis));
168
169// Columns from Tables
170
171 ROArrayColumn<Float> tSysCol;
172 ROScalarColumn<Double> mjdCol;
173 ROScalarColumn<String> srcNameCol;
174 ROScalarColumn<Double> intCol;
175 ROArrayColumn<uInt> fqIDCol;
176 ROScalarColumn<Int> scanIDCol;
177
178// Create accumulation MaskedArray. We accumulate for each channel,if,pol,beam
179// Note that the mask of the accumulation array will ALWAYS remain ALL True.
180// The MA is only used so that when data which is masked Bad is added to it,
181// that data does not contribute.
182
183 Array<Float> zero(shp);
184 zero=0.0;
185 Array<Bool> good(shp);
186 good = True;
187 MaskedArray<Float> sum(zero,good);
188
189// Counter arrays
190
191 Array<Float> nPts(shp); // Number of points
192 nPts = 0.0;
193 Array<Float> nInc(shp); // Increment
194 nInc = 1.0;
195
196// Create accumulation Array for variance. We accumulate for
197// each if,pol,beam, but average over channel. So we need
198// a shape with one less axis dropping channels.
199
200 const uInt nAxesSub = shp.nelements() - 1;
201 IPosition shp2(nAxesSub);
202 for (uInt i=0,j=0; i<(nAxesSub+1); i++) {
203 if (i!=asap::ChanAxis) {
204 shp2(j) = shp(i);
205 j++;
206 }
207 }
208 Array<Float> sumSq(shp2);
209 sumSq = 0.0;
210 IPosition pos2(nAxesSub,0); // For indexing
211
212// Time-related accumulators
213
214 Double time;
215 Double timeSum = 0.0;
216 Double intSum = 0.0;
217 Double interval = 0.0;
218
219// To get the right shape for the Tsys accumulator we need to
220// access a column from the first table. The shape of this
221// array must not change
222
223 Array<Float> tSysSum;
224 {
225 const Table& tabIn = in[0]->table();
226 tSysCol.attach(tabIn,"TSYS");
227 tSysSum.resize(tSysCol.shape(0));
228 }
229 tSysSum =0.0;
230 Array<Float> tSys;
231
232// Scan and row tracking
233
234 Int oldScanID = 0;
235 Int outScanID = 0;
236 Int scanID = 0;
237 Int rowStart = 0;
238 Int nAccum = 0;
239 Int tableStart = 0;
240
241// Source and FreqID
242
243 String sourceName, oldSourceName, sourceNameStart;
244 Vector<uInt> freqID, freqIDStart, oldFreqID;
245
246// Loop over tables
247
248 Float fac = 1.0;
249 const uInt nTables = in.nelements();
250 for (uInt iTab=0; iTab<nTables; iTab++) {
251
252// Should check that the frequency tables don't change if doing FreqAlignment
253
254// Attach columns to Table
255
256 const Table& tabIn = in[iTab]->table();
257 tSysCol.attach(tabIn, "TSYS");
258 mjdCol.attach(tabIn, "TIME");
259 srcNameCol.attach(tabIn, "SRCNAME");
260 intCol.attach(tabIn, "INTERVAL");
261 fqIDCol.attach(tabIn, "FREQID");
262 scanIDCol.attach(tabIn, "SCANID");
263
264// Find list of start/end rows for each scan
265
266// Loop over rows in Table
267
268 const uInt nRows = in[iTab]->nRow();
269 for (uInt iRow=0; iRow<nRows; iRow++) {
270
271// Check conformance
272
273 IPosition shp2 = in[iTab]->rowAsMaskedArray(iRow).shape();
274 if (!shp.isEqual(shp2)) {
275 throw (AipsError("Shapes for all rows must be the same"));
276 }
277
278// If we are not doing scan averages, make checks for source and
279// frequency setup and warn if averaging across them
280
281 scanIDCol.getScalar(iRow, scanID);
282
283// Get quantities from columns
284
285 srcNameCol.getScalar(iRow, sourceName);
286 mjdCol.get(iRow, time);
287 tSysCol.get(iRow, tSys);
288 intCol.get(iRow, interval);
289 fqIDCol.get(iRow, freqID);
290
291// Initialize first source and freqID
292
293 if (iRow==0 && iTab==0) {
294 sourceNameStart = sourceName;
295 freqIDStart = freqID;
296 }
297
298// If we are doing scan averages, see if we are at the end of an
299// accumulation period (scan). We must check soutce names too,
300// since we might have two tables with one scan each but different
301// source names; we shouldn't average different sources together
302
303 if (scanAv && ( (scanID != oldScanID) ||
304 (iRow==0 && iTab>0 && sourceName!=oldSourceName))) {
305
306// Normalize data in 'sum' accumulation array according to weighting scheme
307
308 normalize(sum, sumSq, nPts, wtType, asap::ChanAxis, nAxesSub);
309
310// Get ScanContainer for the first row of this averaged Scan
311
312 SDContainer scOut = in[iTab]->getSDContainer(rowStart);
313
314// Fill scan container. The source and freqID come from the
315// first row of the first table that went into this average (
316// should be the same for all rows in the scan average)
317
318 Float nR(nAccum);
319 fillSDC(scOut, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
320 timeSum/nR, intSum, sourceNameStart, freqIDStart);
321
322// Write container out to Table
323
324 pTabOut->putSDContainer(scOut);
325
326// Reset accumulators
327
328 sum = 0.0;
329 sumSq = 0.0;
330 nAccum = 0;
331//
332 tSysSum =0.0;
333 timeSum = 0.0;
334 intSum = 0.0;
335 nPts = 0.0;
336
337// Increment
338
339 rowStart = iRow; // First row for next accumulation
340 tableStart = iTab; // First table for next accumulation
341 sourceNameStart = sourceName; // First source name for next accumulation
342 freqIDStart = freqID; // First FreqID for next accumulation
343//
344 oldScanID = scanID;
345 outScanID += 1; // Scan ID for next accumulation period
346 }
347
348// Accumulate
349
350 accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts, tSysSum,
351 tSys, nInc, mask, time, interval, in, iTab, iRow, asap::ChanAxis,
352 nAxesSub, useMask, wtType);
353//
354 oldSourceName = sourceName;
355 oldFreqID = freqID;
356 }
357 }
358
359// OK at this point we have accumulation data which is either
360// - accumulated from all tables into one row
361// or
362// - accumulated from the last scan average
363//
364// Normalize data in 'sum' accumulation array according to weighting scheme
365
366 normalize(sum, sumSq, nPts, wtType, asap::ChanAxis, nAxesSub);
367
368// Create and fill container. The container we clone will be from
369// the last Table and the first row that went into the current
370// accumulation. It probably doesn't matter that much really...
371
372 Float nR(nAccum);
373 SDContainer scOut = in[tableStart]->getSDContainer(rowStart);
374 fillSDC(scOut, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
375 timeSum/nR, intSum, sourceNameStart, freqIDStart);
376 pTabOut->putSDContainer(scOut);
377 pTabOut->resetCursor();
378//
379 return CountedPtr<SDMemTable>(pTabOut);
380}
381
382
383
384CountedPtr<SDMemTable> SDMath::binaryOperate(const CountedPtr<SDMemTable>&
385 left,
386 const CountedPtr<SDMemTable>&
387 right,
388 const String& op, Bool preserve,
389 Bool doTSys) const
390{
391
392// Check operator
393
394 String op2(op);
395 op2.upcase();
396 uInt what = 0;
397 if (op2=="ADD") {
398 what = 0;
399 } else if (op2=="SUB") {
400 what = 1;
401 } else if (op2=="MUL") {
402 what = 2;
403 } else if (op2=="DIV") {
404 what = 3;
405 } else if (op2=="QUOTIENT") {
406 what = 4;
407 doTSys = True;
408 } else {
409 throw( AipsError("Unrecognized operation"));
410 }
411
412// Check rows
413
414 const uInt nRowLeft = left->nRow();
415 const uInt nRowRight = right->nRow();
416 Bool ok = (nRowRight==1&&nRowLeft>0) ||
417 (nRowRight>=1&&nRowLeft==nRowRight);
418 if (!ok) {
419 throw (AipsError("The right Scan Table can have one row or the same number of rows as the left Scan Table"));
420 }
421
422// Input Tables
423
424 const Table& tLeft = left->table();
425 const Table& tRight = right->table();
426
427// TSys columns
428
429 ROArrayColumn<Float> tSysLeftCol, tSysRightCol;
430 if (doTSys) {
431 tSysLeftCol.attach(tLeft, "TSYS");
432 tSysRightCol.attach(tRight, "TSYS");
433 }
434
435// First row for right
436
437 Array<Float> tSysLeftArr, tSysRightArr;
438 if (doTSys) tSysRightCol.get(0, tSysRightArr);
439 MaskedArray<Float>* pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(0));
440 IPosition shpRight = pMRight->shape();
441
442// Output Table cloned from left
443
444 SDMemTable* pTabOut = new SDMemTable(*left, True);
445 pTabOut->appendToHistoryTable(right->getHistoryTable());
446// Loop over rows
447
448 for (uInt i=0; i<nRowLeft; i++) {
449
450// Get data
451
452 MaskedArray<Float> mLeft(left->rowAsMaskedArray(i));
453 IPosition shpLeft = mLeft.shape();
454 if (doTSys) tSysLeftCol.get(i, tSysLeftArr);
455//
456 if (nRowRight>1) {
457 delete pMRight;
458 pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(i));
459 shpRight = pMRight->shape();
460 if (doTSys) tSysRightCol.get(i, tSysRightArr);
461 }
462//
463 if (!shpRight.isEqual(shpLeft)) {
464 throw(AipsError("left and right scan tables are not conformant"));
465 }
466 if (doTSys) {
467 if (!tSysRightArr.shape().isEqual(tSysRightArr.shape())) {
468 throw(AipsError("left and right Tsys data are not conformant"));
469 }
470 if (!shpRight.isEqual(tSysRightArr.shape())) {
471 throw(AipsError("left and right scan tables are not conformant"));
472 }
473 }
474
475// Make container
476
477 SDContainer sc = left->getSDContainer(i);
478
479// Operate on data and TSys
480
481 if (what==0) {
482 MaskedArray<Float> tmp = mLeft + *pMRight;
483 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
484 if (doTSys) sc.putTsys(tSysLeftArr+tSysRightArr);
485 } else if (what==1) {
486 MaskedArray<Float> tmp = mLeft - *pMRight;
487 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
488 if (doTSys) sc.putTsys(tSysLeftArr-tSysRightArr);
489 } else if (what==2) {
490 MaskedArray<Float> tmp = mLeft * *pMRight;
491 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
492 if (doTSys) sc.putTsys(tSysLeftArr*tSysRightArr);
493 } else if (what==3) {
494 MaskedArray<Float> tmp = mLeft / *pMRight;
495 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
496 if (doTSys) sc.putTsys(tSysLeftArr/tSysRightArr);
497 } else if (what==4) {
498 if (preserve) {
499 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) -
500 tSysRightArr;
501 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
502 } else {
503 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) -
504 tSysLeftArr;
505 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
506 }
507 sc.putTsys(tSysRightArr);
508 }
509
510// Put new row in output Table
511
512 pTabOut->putSDContainer(sc);
513 }
514 if (pMRight) delete pMRight;
515 pTabOut->resetCursor();
516
517 return CountedPtr<SDMemTable>(pTabOut);
518}
519
520
521
522std::vector<float> SDMath::statistic(const CountedPtr<SDMemTable>& in,
523 const Vector<Bool>& mask,
524 const String& which, Int row) const
525//
526// Perhaps iteration over pol/beam/if should be in here
527// and inside the nrow iteration ?
528//
529{
530 const uInt nRow = in->nRow();
531
532// Specify cursor location
533
534 IPosition start, end;
535 Bool doAll = False;
536 setCursorSlice (start, end, doAll, *in);
537
538// Loop over rows
539
540 const uInt nEl = mask.nelements();
541 uInt iStart = 0;
542 uInt iEnd = in->nRow()-1;
543//
544 if (row>=0) {
545 iStart = row;
546 iEnd = row;
547 }
548//
549 std::vector<float> result(iEnd-iStart+1);
550 for (uInt ii=iStart; ii <= iEnd; ++ii) {
551
552// Get row and deconstruct
553
554 MaskedArray<Float> dataIn = (in->rowAsMaskedArray(ii))(start,end);
555 Array<Float> v = dataIn.getArray().nonDegenerate();
556 Array<Bool> m = dataIn.getMask().nonDegenerate();
557
558// Access desired piece of data
559
560// Array<Float> v((arr(start,end)).nonDegenerate());
561// Array<Bool> m((barr(start,end)).nonDegenerate());
562
563// Apply OTF mask
564
565 MaskedArray<Float> tmp;
566 if (m.nelements()==nEl) {
567 tmp.setData(v,m&&mask);
568 } else {
569 tmp.setData(v,m);
570 }
571
572// Get statistic
573
574 result[ii-iStart] = mathutil::statistics(which, tmp);
575 }
576//
577 return result;
578}
579
580
581SDMemTable* SDMath::bin(const SDMemTable& in, Int width) const
582{
583 SDHeader sh = in.getSDHeader();
584 SDMemTable* pTabOut = new SDMemTable(in, True);
585
586// Bin up SpectralCoordinates
587
588 IPosition factors(1);
589 factors(0) = width;
590 for (uInt j=0; j<in.nCoordinates(); ++j) {
591 CoordinateSystem cSys;
592 cSys.addCoordinate(in.getSpectralCoordinate(j));
593 CoordinateSystem cSysBin =
594 CoordinateUtil::makeBinnedCoordinateSystem(factors, cSys, False);
595//
596 SpectralCoordinate sCBin = cSysBin.spectralCoordinate(0);
597 pTabOut->setCoordinate(sCBin, j);
598 }
599
600// Use RebinLattice to find shape
601
602 IPosition shapeIn(1,sh.nchan);
603 IPosition shapeOut = RebinLattice<Float>::rebinShape(shapeIn, factors);
604 sh.nchan = shapeOut(0);
605 pTabOut->putSDHeader(sh);
606
607// Loop over rows and bin along channel axis
608
609 for (uInt i=0; i < in.nRow(); ++i) {
610 SDContainer sc = in.getSDContainer(i);
611//
612 Array<Float> tSys(sc.getTsys()); // Get it out before sc changes shape
613
614// Bin up spectrum
615
616 MaskedArray<Float> marr(in.rowAsMaskedArray(i));
617 MaskedArray<Float> marrout;
618 LatticeUtilities::bin(marrout, marr, asap::ChanAxis, width);
619
620// Put back the binned data and flags
621
622 IPosition ip2 = marrout.shape();
623 sc.resize(ip2);
624//
625 putDataInSDC(sc, marrout.getArray(), marrout.getMask());
626
627// Bin up Tsys.
628
629 Array<Bool> allGood(tSys.shape(),True);
630 MaskedArray<Float> tSysIn(tSys, allGood, True);
631//
632 MaskedArray<Float> tSysOut;
633 LatticeUtilities::bin(tSysOut, tSysIn, asap::ChanAxis, width);
634 sc.putTsys(tSysOut.getArray());
635//
636 pTabOut->putSDContainer(sc);
637 }
638 return pTabOut;
639}
640
641SDMemTable* SDMath::resample(const SDMemTable& in, const String& methodStr,
642 Float width) const
643//
644// Should add the possibility of width being specified in km/s. This means
645// that for each freqID (SpectralCoordinate) we will need to convert to an
646// average channel width (say at the reference pixel). Then we would need
647// to be careful to make sure each spectrum (of different freqID)
648// is the same length.
649//
650{
651 Bool doVel = False;
652 if (doVel) {
653 for (uInt j=0; j<in.nCoordinates(); ++j) {
654 SpectralCoordinate sC = in.getSpectralCoordinate(j);
655 }
656 }
657
658// Interpolation method
659
660 InterpolateArray1D<Double,Float>::InterpolationMethod interp;
661 convertInterpString(interp, methodStr);
662 Int interpMethod(interp);
663
664// Make output table
665
666 SDMemTable* pTabOut = new SDMemTable(in, True);
667
668// Resample SpectralCoordinates (one per freqID)
669
670 const uInt nCoord = in.nCoordinates();
671 Vector<Float> offset(1,0.0);
672 Vector<Float> factors(1,1.0/width);
673 Vector<Int> newShape;
674 for (uInt j=0; j<in.nCoordinates(); ++j) {
675 CoordinateSystem cSys;
676 cSys.addCoordinate(in.getSpectralCoordinate(j));
677 CoordinateSystem cSys2 = cSys.subImage(offset, factors, newShape);
678 SpectralCoordinate sC = cSys2.spectralCoordinate(0);
679//
680 pTabOut->setCoordinate(sC, j);
681 }
682
683// Get header
684
685 SDHeader sh = in.getSDHeader();
686
687// Generate resampling vectors
688
689 const uInt nChanIn = sh.nchan;
690 Vector<Float> xIn(nChanIn);
691 indgen(xIn);
692//
693 Int fac = Int(nChanIn/width);
694 Vector<Float> xOut(fac+10); // 10 to be safe - resize later
695 uInt i = 0;
696 Float x = 0.0;
697 Bool more = True;
698 while (more) {
699 xOut(i) = x;
700//
701 i++;
702 x += width;
703 if (x>nChanIn-1) more = False;
704 }
705 const uInt nChanOut = i;
706 xOut.resize(nChanOut,True);
707//
708 IPosition shapeIn(in.rowAsMaskedArray(0).shape());
709 sh.nchan = nChanOut;
710 pTabOut->putSDHeader(sh);
711
712// Loop over rows and resample along channel axis
713
714 Array<Float> valuesOut;
715 Array<Bool> maskOut;
716 Array<Float> tSysOut;
717 Array<Bool> tSysMaskIn(shapeIn,True);
718 Array<Bool> tSysMaskOut;
719 for (uInt i=0; i < in.nRow(); ++i) {
720
721// Get container
722
723 SDContainer sc = in.getSDContainer(i);
724
725// Get data and Tsys
726
727 const Array<Float>& tSysIn = sc.getTsys();
728 const MaskedArray<Float>& dataIn(in.rowAsMaskedArray(i));
729 Array<Float> valuesIn = dataIn.getArray();
730 Array<Bool> maskIn = dataIn.getMask();
731
732// Interpolate data
733
734 InterpolateArray1D<Float,Float>::interpolate(valuesOut, maskOut, xOut,
735 xIn, valuesIn, maskIn,
736 interpMethod, True, True);
737 sc.resize(valuesOut.shape());
738 putDataInSDC(sc, valuesOut, maskOut);
739
740// Interpolate TSys
741
742 InterpolateArray1D<Float,Float>::interpolate(tSysOut, tSysMaskOut, xOut,
743 xIn, tSysIn, tSysMaskIn,
744 interpMethod, True, True);
745 sc.putTsys(tSysOut);
746
747// Put container in output
748
749 pTabOut->putSDContainer(sc);
750 }
751//
752 return pTabOut;
753}
754
755SDMemTable* SDMath::unaryOperate(const SDMemTable& in, Float val, Bool doAll,
756 uInt what, Bool doTSys) const
757//
758// what = 0 Multiply
759// 1 Add
760{
761 SDMemTable* pOut = new SDMemTable(in,False);
762 const Table& tOut = pOut->table();
763 ArrayColumn<Float> specCol(tOut,"SPECTRA");
764 ArrayColumn<Float> tSysCol(tOut,"TSYS");
765 Array<Float> tSysArr;
766
767// Get data slice bounds
768
769 IPosition start, end;
770 setCursorSlice (start, end, doAll, in);
771//
772 for (uInt i=0; i<tOut.nrow(); i++) {
773
774// Modify data
775
776 MaskedArray<Float> dataIn(pOut->rowAsMaskedArray(i));
777 MaskedArray<Float> dataIn2 = dataIn(start,end); // Reference
778 if (what==0) {
779 dataIn2 *= val;
780 } else if (what==1) {
781 dataIn2 += val;
782 }
783 specCol.put(i, dataIn.getArray());
784
785// Modify Tsys
786
787 if (doTSys) {
788 tSysCol.get(i, tSysArr);
789 Array<Float> tSysArr2 = tSysArr(start,end); // Reference
790 if (what==0) {
791 tSysArr2 *= val;
792 } else if (what==1) {
793 tSysArr2 += val;
794 }
795 tSysCol.put(i, tSysArr);
796 }
797 }
798//
799 return pOut;
800}
801
802SDMemTable* SDMath::averagePol(const SDMemTable& in, const Vector<Bool>& mask,
803 const String& weightStr) const
804//
805// Average all polarizations together, weighted by variance
806//
807{
808 WeightType wtType = NONE;
809 convertWeightString(wtType, weightStr);
810
811 const uInt nRows = in.nRow();
812
813// Create output Table and reshape number of polarizations
814
815 Bool clear=True;
816 SDMemTable* pTabOut = new SDMemTable(in, clear);
817 SDHeader header = pTabOut->getSDHeader();
818 header.npol = 1;
819 pTabOut->putSDHeader(header);
820
821// Shape of input and output data
822
823 const IPosition& shapeIn = in.rowAsMaskedArray(0).shape();
824 IPosition shapeOut(shapeIn);
825 shapeOut(asap::PolAxis) = 1; // Average all polarizations
826 if (shapeIn(asap::PolAxis)==1) {
827 throw(AipsError("The input has only one polarisation"));
828 }
829//
830 const uInt nChan = shapeIn(asap::ChanAxis);
831 const IPosition vecShapeOut(4,1,1,1,nChan); // A multi-dim form of a Vector shape
832 IPosition start(4), end(4);
833
834// Output arrays
835
836 Array<Float> outData(shapeOut, 0.0);
837 Array<Bool> outMask(shapeOut, True);
838 const IPosition axes(2, asap::PolAxis, asap::ChanAxis); // pol-channel plane
839//
840 const Bool useMask = (mask.nelements() == shapeIn(asap::ChanAxis));
841
842// Loop over rows
843
844 for (uInt iRow=0; iRow<nRows; iRow++) {
845
846// Get data for this row
847
848 MaskedArray<Float> marr(in.rowAsMaskedArray(iRow));
849 Array<Float>& arr = marr.getRWArray();
850 const Array<Bool>& barr = marr.getMask();
851
852// Make iterators to iterate by pol-channel planes
853
854 ReadOnlyArrayIterator<Float> itDataPlane(arr, axes);
855 ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes);
856
857// Accumulations
858
859 Float fac = 1.0;
860 Vector<Float> vecSum(nChan,0.0);
861
862// Iterate through data by pol-channel planes
863
864 while (!itDataPlane.pastEnd()) {
865
866// Iterate through plane by polarization and accumulate Vectors
867
868 Vector<Float> t1(nChan); t1 = 0.0;
869 Vector<Bool> t2(nChan); t2 = True;
870 MaskedArray<Float> vecSum(t1,t2);
871 Float norm = 0.0;
872 {
873 ReadOnlyVectorIterator<Float> itDataVec(itDataPlane.array(), 1);
874 ReadOnlyVectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
875 while (!itDataVec.pastEnd()) {
876
877// Create MA of data & mask (optionally including OTF mask) and get variance for this spectrum
878
879 if (useMask) {
880 const MaskedArray<Float> spec(itDataVec.vector(),mask&&itMaskVec.vector());
881 if (wtType==VAR) fac = 1.0 / variance(spec);
882 } else {
883 const MaskedArray<Float> spec(itDataVec.vector(),itMaskVec.vector());
884 if (wtType==VAR) fac = 1.0 / variance(spec);
885 }
886
887// Normalize spectrum (without OTF mask) and accumulate
888
889 const MaskedArray<Float> spec(fac*itDataVec.vector(), itMaskVec.vector());
890 vecSum += spec;
891 norm += fac;
892
893// Next
894
895 itDataVec.next();
896 itMaskVec.next();
897 }
898 }
899
900// Normalize summed spectrum
901
902 vecSum /= norm;
903
904// FInd position in input data array. We are iterating by pol-channel
905// plane so all that will change is beam and IF and that's what we want.
906
907 IPosition pos = itDataPlane.pos();
908
909// Write out data. This is a bit messy. We have to reform the Vector
910// accumulator into an Array of shape (1,1,1,nChan)
911
912 start = pos;
913 end = pos;
914 end(asap::ChanAxis) = nChan-1;
915 outData(start,end) = vecSum.getArray().reform(vecShapeOut);
916 outMask(start,end) = vecSum.getMask().reform(vecShapeOut);
917
918// Step to next beam/IF combination
919
920 itDataPlane.next();
921 itMaskPlane.next();
922 }
923
924// Generate output container and write it to output table
925
926 SDContainer sc = in.getSDContainer();
927 sc.resize(shapeOut);
928//
929 putDataInSDC(sc, outData, outMask);
930 pTabOut->putSDContainer(sc);
931 }
932
933// Set polarization cursor to 0
934
935 pTabOut->setPol(0);
936//
937 return pTabOut;
938}
939
940
941SDMemTable* SDMath::smooth(const SDMemTable& in,
942 const casa::String& kernelType,
943 casa::Float width, Bool doAll) const
944//
945// Should smooth TSys as well
946//
947{
948
949// Number of channels
950
951 const uInt nChan = in.nChan();
952
953// Generate Kernel
954
955 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernelType);
956 Vector<Float> kernel = VectorKernel::make(type, width, nChan, True, False);
957
958// Generate Convolver
959
960 IPosition shape(1,nChan);
961 Convolver<Float> conv(kernel, shape);
962
963// New Table
964
965 SDMemTable* pTabOut = new SDMemTable(in,True);
966
967// Output Vectors
968
969 Vector<Float> valuesOut(nChan);
970 Vector<Bool> maskOut(nChan);
971
972// Get data slice bounds
973
974 IPosition start, end;
975 setCursorSlice (start, end, doAll, in);
976
977// Loop over rows in Table
978
979 for (uInt ri=0; ri < in.nRow(); ++ri) {
980
981// Get slice of data
982
983 MaskedArray<Float> dataIn = in.rowAsMaskedArray(ri);
984
985// Deconstruct and get slices which reference these arrays
986
987 Array<Float> valuesIn = dataIn.getArray();
988 Array<Bool> maskIn = dataIn.getMask();
989//
990 Array<Float> valuesIn2 = valuesIn(start,end); // ref to valuesIn
991 Array<Bool> maskIn2 = maskIn(start,end);
992
993// Iterate through by spectra
994
995 VectorIterator<Float> itValues(valuesIn2, asap::ChanAxis);
996 VectorIterator<Bool> itMask(maskIn2, asap::ChanAxis);
997 while (!itValues.pastEnd()) {
998
999// Smooth
1000
1001 if (kernelType==VectorKernel::HANNING) {
1002 mathutil::hanning(valuesOut, maskOut, itValues.vector(), itMask.vector());
1003 itMask.vector() = maskOut;
1004 } else {
1005 mathutil::replaceMaskByZero(itValues.vector(), itMask.vector());
1006 conv.linearConv(valuesOut, itValues.vector());
1007 }
1008//
1009 itValues.vector() = valuesOut;
1010//
1011 itValues.next();
1012 itMask.next();
1013 }
1014
1015// Create and put back
1016
1017 SDContainer sc = in.getSDContainer(ri);
1018 putDataInSDC(sc, valuesIn, maskIn);
1019//
1020 pTabOut->putSDContainer(sc);
1021 }
1022//
1023 return pTabOut;
1024}
1025
1026
1027
1028SDMemTable* SDMath::convertFlux(const SDMemTable& in, Float D, Float etaAp,
1029 Float JyPerK, Bool doAll) const
1030//
1031// etaAp = aperture efficiency (-1 means find)
1032// D = geometric diameter (m) (-1 means find)
1033// JyPerK
1034//
1035{
1036 SDHeader sh = in.getSDHeader();
1037 SDMemTable* pTabOut = new SDMemTable(in, True);
1038
1039// Find out how to convert values into Jy and K (e.g. units might be mJy or mK)
1040// Also automatically find out what we are converting to according to the
1041// flux unit
1042
1043 Unit fluxUnit(sh.fluxunit);
1044 Unit K(String("K"));
1045 Unit JY(String("Jy"));
1046//
1047 Bool toKelvin = True;
1048 Double cFac = 1.0;
1049 if (fluxUnit==JY) {
1050 cout << "Converting to K" << endl;
1051//
1052 Quantum<Double> t(1.0,fluxUnit);
1053 Quantum<Double> t2 = t.get(JY);
1054 cFac = (t2 / t).getValue(); // value to Jy
1055//
1056 toKelvin = True;
1057 sh.fluxunit = "K";
1058 } else if (fluxUnit==K) {
1059 cout << "Converting to Jy" << endl;
1060//
1061 Quantum<Double> t(1.0,fluxUnit);
1062 Quantum<Double> t2 = t.get(K);
1063 cFac = (t2 / t).getValue(); // value to K
1064//
1065 toKelvin = False;
1066 sh.fluxunit = "Jy";
1067 } else {
1068 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
1069 }
1070 pTabOut->putSDHeader(sh);
1071
1072// Make sure input values are converted to either Jy or K first...
1073
1074 Float factor = cFac;
1075
1076// Select method
1077
1078 if (JyPerK>0.0) {
1079 factor *= JyPerK;
1080 if (toKelvin) factor = 1.0 / JyPerK;
1081//
1082 cout << "Jy/K = " << JyPerK << endl;
1083 Vector<Float> factors(in.nRow(), factor);
1084 scaleByVector(pTabOut, in, doAll, factors, False);
1085 } else if (etaAp>0.0) {
1086 Bool throwIt = True;
1087 Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt);
1088 SDAttr sda;
1089 if (D < 0) D = sda.diameter(inst);
1090 Float JyPerK = SDAttr::findJyPerK (etaAp,D);
1091 cout << "Jy/K = " << JyPerK << endl;
1092 factor *= JyPerK;
1093 if (toKelvin) {
1094 factor = 1.0 / factor;
1095 }
1096//
1097 Vector<Float> factors(in.nRow(), factor);
1098 scaleByVector(pTabOut, in, doAll, factors, False);
1099 } else {
1100
1101// OK now we must deal with automatic look up of values.
1102// We must also deal with the fact that the factors need
1103// to be computed per IF and may be different and may
1104// change per integration.
1105
1106 cout << "Looking up conversion factors" << endl;
1107 convertBrightnessUnits (pTabOut, in, toKelvin, cFac, doAll);
1108 }
1109//
1110 return pTabOut;
1111}
1112
1113
1114
1115
1116
1117SDMemTable* SDMath::gainElevation(const SDMemTable& in,
1118 const Vector<Float>& coeffs,
1119 const String& fileName,
1120 const String& methodStr, Bool doAll) const
1121{
1122
1123// Get header and clone output table
1124
1125 SDHeader sh = in.getSDHeader();
1126 SDMemTable* pTabOut = new SDMemTable(in, True);
1127
1128// Get elevation data from SDMemTable and convert to degrees
1129
1130 const Table& tab = in.table();
1131 ROScalarColumn<Float> elev(tab, "ELEVATION");
1132 Vector<Float> x = elev.getColumn();
1133 x *= Float(180 / C::pi); // Degrees
1134//
1135 const uInt nC = coeffs.nelements();
1136 if (fileName.length()>0 && nC>0) {
1137 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
1138 }
1139
1140// Correct
1141
1142 if (nC>0 || fileName.length()==0) {
1143
1144// Find instrument
1145
1146 Bool throwIt = True;
1147 Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt);
1148
1149// Set polynomial
1150
1151 Polynomial<Float>* pPoly = 0;
1152 Vector<Float> coeff;
1153 String msg;
1154 if (nC>0) {
1155 pPoly = new Polynomial<Float>(nC);
1156 coeff = coeffs;
1157 msg = String("user");
1158 } else {
1159 SDAttr sdAttr;
1160 coeff = sdAttr.gainElevationPoly(inst);
1161 pPoly = new Polynomial<Float>(3);
1162 msg = String("built in");
1163 }
1164//
1165 if (coeff.nelements()>0) {
1166 pPoly->setCoefficients(coeff);
1167 } else {
1168 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
1169 }
1170//
1171 cout << "Making polynomial correction with " << msg << " coefficients" << endl;
1172 const uInt nRow = in.nRow();
1173 Vector<Float> factor(nRow);
1174 for (uInt i=0; i<nRow; i++) {
1175 factor[i] = 1.0 / (*pPoly)(x[i]);
1176 }
1177 delete pPoly;
1178//
1179 scaleByVector (pTabOut, in, doAll, factor, True);
1180 } else {
1181
1182// Indicate which columns to read from ascii file
1183
1184 String col0("ELEVATION");
1185 String col1("FACTOR");
1186
1187// Read and correct
1188
1189 cout << "Making correction from ascii Table" << endl;
1190 scaleFromAsciiTable (pTabOut, in, fileName, col0, col1,
1191 methodStr, doAll, x, True);
1192 }
1193//
1194 return pTabOut;
1195}
1196
1197
1198
1199SDMemTable* SDMath::opacity(const SDMemTable& in, Float tau, Bool doAll) const
1200{
1201
1202// Get header and clone output table
1203
1204 SDHeader sh = in.getSDHeader();
1205 SDMemTable* pTabOut = new SDMemTable(in, True);
1206
1207// Get elevation data from SDMemTable and convert to degrees
1208
1209 const Table& tab = in.table();
1210 ROScalarColumn<Float> elev(tab, "ELEVATION");
1211 Vector<Float> zDist = elev.getColumn();
1212 zDist = Float(C::pi_2) - zDist;
1213
1214// Generate correction factor
1215
1216 const uInt nRow = in.nRow();
1217 Vector<Float> factor(nRow);
1218 Vector<Float> factor2(nRow);
1219 for (uInt i=0; i<nRow; i++) {
1220 factor[i] = exp(tau)/cos(zDist[i]);
1221 }
1222
1223// Correct
1224
1225 scaleByVector (pTabOut, in, doAll, factor, True);
1226//
1227 return pTabOut;
1228}
1229
1230
1231void SDMath::rotateXYPhase(SDMemTable& in, Float value, Bool doAll)
1232//
1233// phase in degrees
1234// Applies to all Beams and IFs
1235// Might want to optionally select on Beam/IF
1236//
1237{
1238 if (in.nPol() != 4) {
1239 throw(AipsError("You must have 4 polarizations to run this function"));
1240 }
1241//
1242 const Table& tabIn = in.table();
1243 ArrayColumn<Float> specCol(tabIn,"SPECTRA");
1244 IPosition start(asap::nAxes,0);
1245 IPosition end(asap::nAxes);
1246
1247// Set cursor slice. Assumes shape the same for all rows
1248
1249 setCursorSlice (start, end, doAll, in);
1250 IPosition start3(start);
1251 start3(asap::PolAxis) = 2; // Real(XY)
1252 IPosition end3(end);
1253 end3(asap::PolAxis) = 2;
1254//
1255 IPosition start4(start);
1256 start4(asap::PolAxis) = 3; // Imag (XY)
1257 IPosition end4(end);
1258 end4(asap::PolAxis) = 3;
1259//
1260 uInt nRow = in.nRow();
1261 Array<Float> data;
1262 for (uInt i=0; i<nRow;++i) {
1263 specCol.get(i,data);
1264 IPosition shape = data.shape();
1265
1266// Get polarization slice references
1267
1268 Array<Float> C3 = data(start3,end3);
1269 Array<Float> C4 = data(start4,end4);
1270
1271// Rotate
1272
1273 SDPolUtil::rotatePhase(C3, C4, value);
1274
1275// Put
1276
1277 specCol.put(i,data);
1278 }
1279}
1280
1281
1282
1283void SDMath::rotateLinPolPhase(SDMemTable& in, Float value, Bool doAll)
1284//
1285// phase in degrees
1286// Applies to all Beams and IFs
1287// Might want to optionally select on Beam/IF
1288//
1289{
1290 if (in.nPol() != 4) {
1291 throw(AipsError("You must have 4 polarizations to run this function"));
1292 }
1293//
1294 const Table& tabIn = in.table();
1295 ArrayColumn<Float> specCol(tabIn,"SPECTRA");
1296 ROArrayColumn<Float> stokesCol(tabIn,"STOKES");
1297 IPosition start(asap::nAxes,0);
1298 IPosition end(asap::nAxes);
1299
1300// Set cursor slice. Assumes shape the same for all rows
1301
1302 setCursorSlice (start, end, doAll, in);
1303//
1304 IPosition start1(start);
1305 start1(asap::PolAxis) = 0; // C1 (XX)
1306 IPosition end1(end);
1307 end1(asap::PolAxis) = 0;
1308//
1309 IPosition start2(start);
1310 start2(asap::PolAxis) = 1; // C2 (YY)
1311 IPosition end2(end);
1312 end2(asap::PolAxis) = 1;
1313//
1314 IPosition start3(start);
1315 start3(asap::PolAxis) = 2; // C3 ( Real(XY) )
1316 IPosition end3(end);
1317 end3(asap::PolAxis) = 2;
1318//
1319 IPosition startI(start);
1320 startI(asap::PolAxis) = 0; // I
1321 IPosition endI(end);
1322 endI(asap::PolAxis) = 0;
1323//
1324 IPosition startQ(start);
1325 startQ(asap::PolAxis) = 1; // Q
1326 IPosition endQ(end);
1327 endQ(asap::PolAxis) = 1;
1328//
1329 IPosition startU(start);
1330 startU(asap::PolAxis) = 2; // U
1331 IPosition endU(end);
1332 endU(asap::PolAxis) = 2;
1333
1334//
1335 uInt nRow = in.nRow();
1336 Array<Float> data, stokes;
1337 for (uInt i=0; i<nRow;++i) {
1338 specCol.get(i,data);
1339 stokesCol.get(i,stokes);
1340 IPosition shape = data.shape();
1341
1342// Get linear polarization slice references
1343
1344 Array<Float> C1 = data(start1,end1);
1345 Array<Float> C2 = data(start2,end2);
1346 Array<Float> C3 = data(start3,end3);
1347
1348// Get STokes slice references
1349
1350 Array<Float> I = stokes(startI,endI);
1351 Array<Float> Q = stokes(startQ,endQ);
1352 Array<Float> U = stokes(startU,endU);
1353
1354// Rotate
1355
1356 SDPolUtil::rotateLinPolPhase(C1, C2, C3, I, Q, U, value);
1357
1358// Put
1359
1360 specCol.put(i,data);
1361 }
1362}
1363
1364// 'private' functions
1365
1366void SDMath::convertBrightnessUnits (SDMemTable* pTabOut, const SDMemTable& in,
1367 Bool toKelvin, Float cFac, Bool doAll) const
1368{
1369
1370// Get header
1371
1372 SDHeader sh = in.getSDHeader();
1373 const uInt nChan = sh.nchan;
1374
1375// Get instrument
1376
1377 Bool throwIt = True;
1378 Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt);
1379
1380// Get Diameter (m)
1381
1382 SDAttr sdAtt;
1383
1384// Get epoch of first row
1385
1386 MEpoch dateObs = in.getEpoch(0);
1387
1388// Generate a Vector of correction factors. One per FreqID
1389
1390 SDFrequencyTable sdft = in.getSDFreqTable();
1391 Vector<uInt> freqIDs;
1392//
1393 Vector<Float> freqs(sdft.length());
1394 for (uInt i=0; i<sdft.length(); i++) {
1395 freqs(i) = (nChan/2 - sdft.referencePixel(i))*sdft.increment(i) + sdft.referenceValue(i);
1396 }
1397//
1398 Vector<Float> JyPerK = sdAtt.JyPerK(inst, dateObs, freqs);
1399 cout << "Jy/K = " << JyPerK << endl;
1400 Vector<Float> factors = cFac * JyPerK;
1401 if (toKelvin) factors = Float(1.0) / factors;
1402
1403// Get data slice bounds
1404
1405 IPosition start, end;
1406 setCursorSlice (start, end, doAll, in);
1407 const uInt ifAxis = in.getIF();
1408
1409// Iteration axes
1410
1411 IPosition axes(asap::nAxes-1,0);
1412 for (uInt i=0,j=0; i<asap::nAxes; i++) {
1413 if (i!=asap::IFAxis) {
1414 axes(j++) = i;
1415 }
1416 }
1417
1418// Loop over rows and apply correction factor
1419
1420 Float factor = 1.0;
1421 const uInt axis = asap::ChanAxis;
1422 for (uInt i=0; i < in.nRow(); ++i) {
1423
1424// Get data
1425
1426 MaskedArray<Float> dataIn = in.rowAsMaskedArray(i);
1427 Array<Float>& values = dataIn.getRWArray(); // Ref to dataIn
1428 Array<Float> values2 = values(start,end); // Ref to values to dataIn
1429
1430// Get SDCOntainer
1431
1432 SDContainer sc = in.getSDContainer(i);
1433
1434// Get FreqIDs
1435
1436 freqIDs = sc.getFreqMap();
1437
1438// Now the conversion factor depends only upon frequency
1439// So we need to iterate through by IF only giving
1440// us BEAM/POL/CHAN cubes
1441
1442 ArrayIterator<Float> itIn(values2, axes);
1443 uInt ax = 0;
1444 while (!itIn.pastEnd()) {
1445 itIn.array() *= factors(freqIDs(ax)); // Writes back to dataIn
1446 itIn.next();
1447 }
1448
1449// Write out
1450
1451 putDataInSDC(sc, dataIn.getArray(), dataIn.getMask());
1452//
1453 pTabOut->putSDContainer(sc);
1454 }
1455}
1456
1457
1458
1459SDMemTable* SDMath::frequencyAlign (const SDMemTable& in,
1460 MFrequency::Types freqSystem,
1461 const String& refTime,
1462 const String& methodStr,
1463 Bool perFreqID) const
1464{
1465// Get Header
1466
1467 SDHeader sh = in.getSDHeader();
1468 const uInt nChan = sh.nchan;
1469 const uInt nRows = in.nRow();
1470 const uInt nIF = sh.nif;
1471
1472// Get Table reference
1473
1474 const Table& tabIn = in.table();
1475
1476// Get Columns from Table
1477
1478 ROScalarColumn<Double> mjdCol(tabIn, "TIME");
1479 ROScalarColumn<String> srcCol(tabIn, "SRCNAME");
1480 ROArrayColumn<uInt> fqIDCol(tabIn, "FREQID");
1481 Vector<Double> times = mjdCol.getColumn();
1482
1483// Generate DataDesc table
1484
1485 Matrix<uInt> ddIdx;
1486 SDDataDesc dDesc;
1487 generateDataDescTable (ddIdx, dDesc, nIF, in, tabIn, srcCol, fqIDCol, perFreqID);
1488
1489// Get reference Epoch to time of first row or given String
1490
1491 Unit DAY(String("d"));
1492 MEpoch::Ref epochRef(in.getTimeReference());
1493 MEpoch refEpoch;
1494 if (refTime.length()>0) {
1495 refEpoch = epochFromString(refTime, in.getTimeReference());
1496 } else {
1497 refEpoch = in.getEpoch(0);
1498 }
1499 cout << "Aligning at reference Epoch " << formatEpoch(refEpoch)
1500 << " in frame " << MFrequency::showType(freqSystem) << endl;
1501
1502// Get Reference Position
1503
1504 MPosition refPos = in.getAntennaPosition();
1505
1506// Create FrequencyAligner Block. One FA for each possible
1507// source/freqID (perFreqID=True) or source/IF (perFreqID=False) combination
1508
1509 PtrBlock<FrequencyAligner<Float>* > a(dDesc.length());
1510 generateFrequencyAligners (a, dDesc, in, nChan, freqSystem, refPos,
1511 refEpoch, perFreqID);
1512
1513// Generate and fill output Frequency Table. WHen perFreqID=True, there is one output FreqID
1514// for each entry in the SDDataDesc table. However, in perFreqID=False mode, there may be
1515// some degeneracy, so we need a little translation map
1516
1517 SDFrequencyTable freqTabOut = in.getSDFreqTable();
1518 freqTabOut.setLength(0);
1519 Vector<String> units(1);
1520 units = String("Hz");
1521 Bool linear=True;
1522//
1523 Vector<uInt> ddFQTrans(dDesc.length(),0);
1524 for (uInt i=0; i<dDesc.length(); i++) {
1525
1526// Get Aligned SC in Hz
1527
1528 SpectralCoordinate sC = a[i]->alignedSpectralCoordinate(linear);
1529 sC.setWorldAxisUnits(units);
1530
1531// Add FreqID
1532
1533 uInt idx = freqTabOut.addFrequency(sC.referencePixel()[0],
1534 sC.referenceValue()[0],
1535 sC.increment()[0]);
1536 ddFQTrans(i) = idx; // output FreqID = ddFQTrans(ddIdx)
1537 }
1538
1539// Interpolation method
1540
1541 InterpolateArray1D<Double,Float>::InterpolationMethod interp;
1542 convertInterpString(interp, methodStr);
1543
1544// New output Table
1545
1546 cout << "Create output table" << endl;
1547 SDMemTable* pTabOut = new SDMemTable(in,True);
1548 pTabOut->putSDFreqTable(freqTabOut);
1549
1550// Loop over rows in Table
1551
1552 Bool extrapolate=False;
1553 const IPosition polChanAxes(2, asap::PolAxis, asap::ChanAxis);
1554 Bool useCachedAbcissa = False;
1555 Bool first = True;
1556 Bool ok;
1557 Vector<Float> yOut;
1558 Vector<Bool> maskOut;
1559 Vector<uInt> freqID(nIF);
1560 uInt ifIdx, faIdx;
1561 Vector<Double> xIn;
1562//
1563 for (uInt iRow=0; iRow<nRows; ++iRow) {
1564 if (iRow%10==0) {
1565 cout << "Processing row " << iRow << endl;
1566 }
1567
1568// Get EPoch
1569
1570 Quantum<Double> tQ2(times[iRow],DAY);
1571 MVEpoch mv2(tQ2);
1572 MEpoch epoch(mv2, epochRef);
1573
1574// Get copy of data
1575
1576 const MaskedArray<Float>& mArrIn(in.rowAsMaskedArray(iRow));
1577 Array<Float> values = mArrIn.getArray();
1578 Array<Bool> mask = mArrIn.getMask();
1579
1580// For each row, the Frequency abcissa will be the same regardless
1581// of polarization. For all other axes (IF and BEAM) the abcissa
1582// will change. So we iterate through the data by pol-chan planes
1583// to mimimize the work. Probably won't work for multiple beams
1584// at this point.
1585
1586 ArrayIterator<Float> itValuesPlane(values, polChanAxes);
1587 ArrayIterator<Bool> itMaskPlane(mask, polChanAxes);
1588 while (!itValuesPlane.pastEnd()) {
1589
1590// Find the IF index and then the FA PtrBlock index
1591
1592 const IPosition& pos = itValuesPlane.pos();
1593 ifIdx = pos(asap::IFAxis);
1594 faIdx = ddIdx(iRow,ifIdx);
1595
1596// Generate abcissa for perIF. Could cache this in a Matrix
1597// on a per scan basis. Pretty expensive doing it for every row.
1598
1599 if (!perFreqID) {
1600 xIn.resize(nChan);
1601 uInt fqID = dDesc.secID(ddIdx(iRow,ifIdx));
1602 SpectralCoordinate sC = in.getSpectralCoordinate(fqID);
1603 Double w;
1604 for (uInt i=0; i<nChan; i++) {
1605 sC.toWorld(w,Double(i));
1606 xIn[i] = w;
1607 }
1608 }
1609//
1610 VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1);
1611 VectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
1612
1613// Iterate through the plane by vector and align
1614
1615 first = True;
1616 useCachedAbcissa=False;
1617 while (!itValuesVec.pastEnd()) {
1618 if (perFreqID) {
1619 ok = a[faIdx]->align (yOut, maskOut, itValuesVec.vector(),
1620 itMaskVec.vector(), epoch, useCachedAbcissa,
1621 interp, extrapolate);
1622 } else {
1623 ok = a[faIdx]->align (yOut, maskOut, xIn, itValuesVec.vector(),
1624 itMaskVec.vector(), epoch, useCachedAbcissa,
1625 interp, extrapolate);
1626 }
1627//
1628 itValuesVec.vector() = yOut;
1629 itMaskVec.vector() = maskOut;
1630//
1631 itValuesVec.next();
1632 itMaskVec.next();
1633//
1634 if (first) {
1635 useCachedAbcissa = True;
1636 first = False;
1637 }
1638 }
1639//
1640 itValuesPlane.next();
1641 itMaskPlane.next();
1642 }
1643
1644// Create SDContainer and put back
1645
1646 SDContainer sc = in.getSDContainer(iRow);
1647 putDataInSDC(sc, values, mask);
1648
1649// Set output FreqIDs
1650
1651 for (uInt i=0; i<nIF; i++) {
1652 uInt idx = ddIdx(iRow,i); // Index into SDDataDesc table
1653 freqID(i) = ddFQTrans(idx); // FreqID in output FQ table
1654 }
1655 sc.putFreqMap(freqID);
1656//
1657 pTabOut->putSDContainer(sc);
1658 }
1659
1660// Now we must set the base and extra frames to the
1661// input frame
1662
1663 std::vector<string> info = pTabOut->getCoordInfo();
1664 info[1] = MFrequency::showType(freqSystem); // Conversion frame
1665 info[3] = info[1]; // Base frame
1666 pTabOut->setCoordInfo(info);
1667
1668// Clean up PointerBlock
1669
1670 for (uInt i=0; i<a.nelements(); i++) delete a[i];
1671//
1672 return pTabOut;
1673}
1674
1675
1676void SDMath::fillSDC(SDContainer& sc,
1677 const Array<Bool>& mask,
1678 const Array<Float>& data,
1679 const Array<Float>& tSys,
1680 Int scanID, Double timeStamp,
1681 Double interval, const String& sourceName,
1682 const Vector<uInt>& freqID) const
1683{
1684// Data and mask
1685
1686 putDataInSDC(sc, data, mask);
1687
1688// TSys
1689
1690 sc.putTsys(tSys);
1691
1692// Time things
1693
1694 sc.timestamp = timeStamp;
1695 sc.interval = interval;
1696 sc.scanid = scanID;
1697//
1698 sc.sourcename = sourceName;
1699 sc.putFreqMap(freqID);
1700}
1701
1702void SDMath::normalize(MaskedArray<Float>& sum,
1703 const Array<Float>& sumSq,
1704 const Array<Float>& nPts,
1705 WeightType wtType, Int axis,
1706 Int nAxesSub) const
1707{
1708 IPosition pos2(nAxesSub,0);
1709//
1710 if (wtType==NONE) {
1711
1712// We just average by the number of points accumulated.
1713// We need to make a MA out of nPts so that no divide by
1714// zeros occur
1715
1716 MaskedArray<Float> t(nPts, (nPts>Float(0.0)));
1717 sum /= t;
1718 } else if (wtType==VAR) {
1719
1720// Normalize each spectrum by sum(1/var) where the variance
1721// is worked out for each spectrum
1722
1723 Array<Float>& data = sum.getRWArray();
1724 VectorIterator<Float> itData(data, axis);
1725 while (!itData.pastEnd()) {
1726 pos2 = itData.pos().getFirst(nAxesSub);
1727 itData.vector() /= sumSq(pos2);
1728 itData.next();
1729 }
1730 } else if (wtType==TSYS) {
1731 }
1732}
1733
1734
1735void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum,
1736 MaskedArray<Float>& sum, Array<Float>& sumSq,
1737 Array<Float>& nPts, Array<Float>& tSysSum,
1738 const Array<Float>& tSys, const Array<Float>& nInc,
1739 const Vector<Bool>& mask, Double time, Double interval,
1740 const Block<CountedPtr<SDMemTable> >& in,
1741 uInt iTab, uInt iRow, uInt axis,
1742 uInt nAxesSub, Bool useMask,
1743 WeightType wtType) const
1744{
1745
1746// Get data
1747
1748 MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
1749 Array<Float>& valuesIn = dataIn.getRWArray(); // writable reference
1750 const Array<Bool>& maskIn = dataIn.getMask(); // RO reference
1751//
1752 if (wtType==NONE) {
1753 const MaskedArray<Float> n(nInc,dataIn.getMask());
1754 nPts += n; // Only accumulates where mask==T
1755 } else if (wtType==VAR) {
1756
1757// We are going to average the data, weighted by the noise for each pol, beam and IF.
1758// So therefore we need to iterate through by spectrum (axis 3)
1759
1760 VectorIterator<Float> itData(valuesIn, axis);
1761 ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
1762 Float fac = 1.0;
1763 IPosition pos(nAxesSub,0);
1764//
1765 while (!itData.pastEnd()) {
1766
1767// Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
1768
1769 if (useMask) {
1770 MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
1771 fac = 1.0/variance(tmp);
1772 } else {
1773 MaskedArray<Float> tmp(itData.vector(),itMask.vector());
1774 fac = 1.0/variance(tmp);
1775 }
1776
1777// Scale data
1778
1779 itData.vector() *= fac; // Writes back into 'dataIn'
1780//
1781// Accumulate variance per if/pol/beam averaged over spectrum
1782// This method to get pos2 from itData.pos() is only valid
1783// because the spectral axis is the last one (so we can just
1784// copy the first nAXesSub positions out)
1785
1786 pos = itData.pos().getFirst(nAxesSub);
1787 sumSq(pos) += fac;
1788//
1789 itData.next();
1790 itMask.next();
1791 }
1792 } else if (wtType==TSYS) {
1793 }
1794
1795// Accumulate sum of (possibly scaled) data
1796
1797 sum += dataIn;
1798
1799// Accumulate Tsys, time, and interval
1800
1801 tSysSum += tSys;
1802 timeSum += time;
1803 intSum += interval;
1804 nAccum += 1;
1805}
1806
1807
1808void SDMath::setCursorSlice (IPosition& start, IPosition& end, Bool doAll, const SDMemTable& in) const
1809{
1810 const uInt nDim = asap::nAxes;
1811 DebugAssert(nDim==4,AipsError);
1812//
1813 start.resize(nDim);
1814 end.resize(nDim);
1815 if (doAll) {
1816 start = 0;
1817 end(0) = in.nBeam()-1;
1818 end(1) = in.nIF()-1;
1819 end(2) = in.nPol()-1;
1820 end(3) = in.nChan()-1;
1821 } else {
1822 start(0) = in.getBeam();
1823 end(0) = start(0);
1824//
1825 start(1) = in.getIF();
1826 end(1) = start(1);
1827//
1828 start(2) = in.getPol();
1829 end(2) = start(2);
1830//
1831 start(3) = 0;
1832 end(3) = in.nChan()-1;
1833 }
1834}
1835
1836
1837void SDMath::convertWeightString(WeightType& wtType, const String& weightStr) const
1838{
1839 String tStr(weightStr);
1840 tStr.upcase();
1841 if (tStr.contains(String("NONE"))) {
1842 wtType = NONE;
1843 } else if (tStr.contains(String("VAR"))) {
1844 wtType = VAR;
1845 } else if (tStr.contains(String("TSYS"))) {
1846 wtType = TSYS;
1847 throw(AipsError("T_sys weighting not yet implemented"));
1848 } else {
1849 throw(AipsError("Unrecognized weighting type"));
1850 }
1851}
1852
1853
1854void SDMath::convertInterpString(casa::InterpolateArray1D<Double,Float>::InterpolationMethod& type,
1855 const casa::String& interp) const
1856{
1857 String tStr(interp);
1858 tStr.upcase();
1859 if (tStr.contains(String("NEAR"))) {
1860 type = InterpolateArray1D<Double,Float>::nearestNeighbour;
1861 } else if (tStr.contains(String("LIN"))) {
1862 type = InterpolateArray1D<Double,Float>::linear;
1863 } else if (tStr.contains(String("CUB"))) {
1864 type = InterpolateArray1D<Double,Float>::cubic;
1865 } else if (tStr.contains(String("SPL"))) {
1866 type = InterpolateArray1D<Double,Float>::spline;
1867 } else {
1868 throw(AipsError("Unrecognized interpolation type"));
1869 }
1870}
1871
1872void SDMath::putDataInSDC(SDContainer& sc, const Array<Float>& data,
1873 const Array<Bool>& mask) const
1874{
1875 sc.putSpectrum(data);
1876//
1877 Array<uChar> outflags(data.shape());
1878 convertArray(outflags,!mask);
1879 sc.putFlags(outflags);
1880}
1881
1882Table SDMath::readAsciiFile (const String& fileName) const
1883{
1884 String formatString;
1885 Table tbl = readAsciiTable (formatString, Table::Memory, fileName, "", "", False);
1886 return tbl;
1887}
1888
1889
1890
1891void SDMath::scaleFromAsciiTable(SDMemTable* pTabOut,
1892 const SDMemTable& in, const String& fileName,
1893 const String& col0, const String& col1,
1894 const String& methodStr, Bool doAll,
1895 const Vector<Float>& xOut, Bool doTSys) const
1896{
1897
1898// Read gain-elevation ascii file data into a Table.
1899
1900 Table geTable = readAsciiFile (fileName);
1901//
1902 scaleFromTable (pTabOut, in, geTable, col0, col1, methodStr, doAll, xOut, doTSys);
1903}
1904
1905void SDMath::scaleFromTable(SDMemTable* pTabOut, const SDMemTable& in,
1906 const Table& tTable, const String& col0,
1907 const String& col1,
1908 const String& methodStr, Bool doAll,
1909 const Vector<Float>& xOut, Bool doTsys) const
1910{
1911
1912// Get data from Table
1913
1914 ROScalarColumn<Float> geElCol(tTable, col0);
1915 ROScalarColumn<Float> geFacCol(tTable, col1);
1916 Vector<Float> xIn = geElCol.getColumn();
1917 Vector<Float> yIn = geFacCol.getColumn();
1918 Vector<Bool> maskIn(xIn.nelements(),True);
1919
1920// Interpolate (and extrapolate) with desired method
1921
1922 InterpolateArray1D<Double,Float>::InterpolationMethod method;
1923 convertInterpString(method, methodStr);
1924 Int intMethod(method);
1925//
1926 Vector<Float> yOut;
1927 Vector<Bool> maskOut;
1928 InterpolateArray1D<Float,Float>::interpolate(yOut, maskOut, xOut,
1929 xIn, yIn, maskIn, intMethod,
1930 True, True);
1931// Apply
1932
1933 scaleByVector(pTabOut, in, doAll, Float(1.0)/yOut, doTsys);
1934}
1935
1936
1937void SDMath::scaleByVector(SDMemTable* pTabOut, const SDMemTable& in,
1938 Bool doAll, const Vector<Float>& factor,
1939 Bool doTSys) const
1940{
1941
1942// Set up data slice
1943
1944 IPosition start, end;
1945 setCursorSlice (start, end, doAll, in);
1946
1947// Get Tsys column
1948
1949 const Table& tIn = in.table();
1950 ArrayColumn<Float> tSysCol(tIn, "TSYS");
1951 Array<Float> tSys;
1952
1953// Loop over rows and apply correction factor
1954
1955 const uInt axis = asap::ChanAxis;
1956 for (uInt i=0; i < in.nRow(); ++i) {
1957
1958// Get data
1959
1960 MaskedArray<Float> dataIn(in.rowAsMaskedArray(i));
1961 MaskedArray<Float> dataIn2 = dataIn(start,end); // reference to dataIn
1962//
1963 if (doTSys) {
1964 tSysCol.get(i, tSys);
1965 Array<Float> tSys2 = tSys(start,end) * factor[i];
1966 tSysCol.put(i, tSys);
1967 }
1968
1969// Apply factor
1970
1971 dataIn2 *= factor[i];
1972
1973// Write out
1974
1975 SDContainer sc = in.getSDContainer(i);
1976 putDataInSDC(sc, dataIn.getArray(), dataIn.getMask());
1977//
1978 pTabOut->putSDContainer(sc);
1979 }
1980}
1981
1982
1983
1984
1985void SDMath::generateDataDescTable (Matrix<uInt>& ddIdx,
1986 SDDataDesc& dDesc,
1987 uInt nIF,
1988 const SDMemTable& in,
1989 const Table& tabIn,
1990 const ROScalarColumn<String>& srcCol,
1991 const ROArrayColumn<uInt>& fqIDCol,
1992 Bool perFreqID) const
1993{
1994 const uInt nRows = tabIn.nrow();
1995 ddIdx.resize(nRows,nIF);
1996//
1997 String srcName;
1998 Vector<uInt> freqIDs;
1999 for (uInt iRow=0; iRow<nRows; iRow++) {
2000 srcCol.get(iRow, srcName);
2001 fqIDCol.get(iRow, freqIDs);
2002 const MDirection& dir = in.getDirection(iRow);
2003//
2004 if (perFreqID) {
2005
2006// One entry per source/freqID pair
2007
2008 for (uInt iIF=0; iIF<nIF; iIF++) {
2009 ddIdx(iRow,iIF) = dDesc.addEntry(srcName, freqIDs[iIF], dir, 0);
2010 }
2011 } else {
2012
2013// One entry per source/IF pair. Hang onto the FreqID as well
2014
2015 for (uInt iIF=0; iIF<nIF; iIF++) {
2016 ddIdx(iRow,iIF) = dDesc.addEntry(srcName, iIF, dir, freqIDs[iIF]);
2017 }
2018 }
2019 }
2020}
2021
2022
2023
2024
2025
2026MEpoch SDMath::epochFromString (const String& str, MEpoch::Types timeRef) const
2027{
2028 Quantum<Double> qt;
2029 if (MVTime::read(qt,str)) {
2030 MVEpoch mv(qt);
2031 MEpoch me(mv, timeRef);
2032 return me;
2033 } else {
2034 throw(AipsError("Invalid format for Epoch string"));
2035 }
2036}
2037
2038
2039String SDMath::formatEpoch(const MEpoch& epoch) const
2040{
2041 MVTime mvt(epoch.getValue());
2042 return mvt.string(MVTime::YMD) + String(" (") + epoch.getRefString() + String(")");
2043}
2044
2045
2046
2047void SDMath::generateFrequencyAligners (PtrBlock<FrequencyAligner<Float>* >& a,
2048 const SDDataDesc& dDesc,
2049 const SDMemTable& in, uInt nChan,
2050 MFrequency::Types system,
2051 const MPosition& refPos,
2052 const MEpoch& refEpoch,
2053 Bool perFreqID) const
2054{
2055 for (uInt i=0; i<dDesc.length(); i++) {
2056 uInt ID = dDesc.ID(i);
2057 uInt secID = dDesc.secID(i);
2058 const MDirection& refDir = dDesc.secDir(i);
2059//
2060 if (perFreqID) {
2061
2062// One aligner per source/FreqID pair.
2063
2064 SpectralCoordinate sC = in.getSpectralCoordinate(ID);
2065 a[i] = new FrequencyAligner<Float>(sC, nChan, refEpoch, refDir, refPos, system);
2066 } else {
2067
2068// One aligner per source/IF pair. But we still need the FreqID to
2069// get the right SC. Hence the messing about with the secondary ID
2070
2071 SpectralCoordinate sC = in.getSpectralCoordinate(secID);
2072 a[i] = new FrequencyAligner<Float>(sC, nChan, refEpoch, refDir, refPos, system);
2073 }
2074 }
2075}
2076
2077Vector<uInt> SDMath::getRowRange (const SDMemTable& in) const
2078{
2079 Vector<uInt> range(2);
2080 range[0] = 0;
2081 range[1] = in.nRow()-1;
2082 return range;
2083}
2084
2085
2086Bool SDMath::rowInRange (uInt i, const Vector<uInt>& range) const
2087{
2088 return (i>=range[0] && i<=range[1]);
2089}
Note: See TracBrowser for help on using the repository browser.