source: trunk/src/SDMath.cc@ 428

Last change on this file since 428 was 414, checked in by mar637, 20 years ago

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