source: branches/Release12/src/SDMath.cc@ 1437

Last change on this file since 1437 was 799, checked in by phi196, 19 years ago

Bug in swappol, minor indentation reformat

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