source: trunk/src/SDMath.cc@ 434

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

consolidate cursor selection so that when function
setCursorSlice is extended to handle multi-dimensional
cursors (but contiguous pixels, all of the computation
functions will operate correctly.

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