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

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

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