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

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

Added swap_pol & invert_phase
  • 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 *= -1.0f;
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;
1486 tmp = C1;
1487 C1 = C2;
1488 C2 = C1;
1489
1490// Put
1491
1492 specCol.put(i,data);
1493 }
1494}
1495
1496// 'private' functions
1497
1498void SDMath::convertBrightnessUnits (SDMemTable* pTabOut, const SDMemTable& in,
1499 Bool toKelvin, Float cFac, Bool doAll)
1500{
1501
1502// Get header
1503
1504 SDHeader sh = in.getSDHeader();
1505 const uInt nChan = sh.nchan;
1506
1507// Get instrument
1508
1509 Bool throwIt = True;
1510 Instrument inst = SDAttr::convertInstrument(sh.antennaname, throwIt);
1511
1512// Get Diameter (m)
1513
1514 SDAttr sdAtt;
1515// Get epoch of first row
1516
1517 MEpoch dateObs = in.getEpoch(0);
1518
1519// Generate a Vector of correction factors. One per FreqID
1520
1521 SDFrequencyTable sdft = in.getSDFreqTable();
1522 Vector<uInt> freqIDs;
1523//
1524 Vector<Float> freqs(sdft.length());
1525 for (uInt i=0; i<sdft.length(); i++) {
1526 freqs(i) = (nChan/2 - sdft.referencePixel(i))*sdft.increment(i) + sdft.referenceValue(i);
1527 }
1528//
1529 Vector<Float> JyPerK = sdAtt.JyPerK(inst, dateObs, freqs);
1530 ostringstream oss;
1531 oss << "Jy/K = " << JyPerK;
1532 pushLog(String(oss));
1533 Vector<Float> factors = cFac * JyPerK;
1534 if (toKelvin) factors = Float(1.0) / factors;
1535
1536// Get data slice bounds
1537
1538 IPosition start, end;
1539 setCursorSlice (start, end, doAll, in);
1540 const uInt ifAxis = in.getIF();
1541
1542// Iteration axes
1543
1544 IPosition axes(asap::nAxes-1,0);
1545 for (uInt i=0,j=0; i<asap::nAxes; i++) {
1546 if (i!=asap::IFAxis) {
1547 axes(j++) = i;
1548 }
1549 }
1550
1551// Loop over rows and apply correction factor
1552
1553 Float factor = 1.0;
1554 const uInt axis = asap::ChanAxis;
1555 for (uInt i=0; i < in.nRow(); ++i) {
1556
1557// Get data
1558
1559 MaskedArray<Float> dataIn = in.rowAsMaskedArray(i);
1560 Array<Float>& values = dataIn.getRWArray(); // Ref to dataIn
1561 Array<Float> values2 = values(start,end); // Ref to values to dataIn
1562
1563// Get SDCOntainer
1564
1565 SDContainer sc = in.getSDContainer(i);
1566
1567// Get FreqIDs
1568
1569 freqIDs = sc.getFreqMap();
1570
1571// Now the conversion factor depends only upon frequency
1572// So we need to iterate through by IF only giving
1573// us BEAM/POL/CHAN cubes
1574
1575 ArrayIterator<Float> itIn(values2, axes);
1576 uInt ax = 0;
1577 while (!itIn.pastEnd()) {
1578 itIn.array() *= factors(freqIDs(ax)); // Writes back to dataIn
1579 itIn.next();
1580 }
1581
1582// Write out
1583
1584 putDataInSDC(sc, dataIn.getArray(), dataIn.getMask());
1585//
1586 pTabOut->putSDContainer(sc);
1587 }
1588}
1589
1590
1591
1592SDMemTable* SDMath::frequencyAlign(const SDMemTable& in,
1593 MFrequency::Types freqSystem,
1594 const String& refTime,
1595 const String& methodStr,
1596 Bool perFreqID)
1597{
1598// Get Header
1599
1600 SDHeader sh = in.getSDHeader();
1601 const uInt nChan = sh.nchan;
1602 const uInt nRows = in.nRow();
1603 const uInt nIF = sh.nif;
1604
1605// Get Table reference
1606
1607 const Table& tabIn = in.table();
1608
1609// Get Columns from Table
1610
1611 ROScalarColumn<Double> mjdCol(tabIn, "TIME");
1612 ROScalarColumn<String> srcCol(tabIn, "SRCNAME");
1613 ROArrayColumn<uInt> fqIDCol(tabIn, "FREQID");
1614 Vector<Double> times = mjdCol.getColumn();
1615
1616// Generate DataDesc table
1617
1618 Matrix<uInt> ddIdx;
1619 SDDataDesc dDesc;
1620 generateDataDescTable(ddIdx, dDesc, nIF, in, tabIn, srcCol,
1621 fqIDCol, perFreqID);
1622
1623// Get reference Epoch to time of first row or given String
1624
1625 Unit DAY(String("d"));
1626 MEpoch::Ref epochRef(in.getTimeReference());
1627 MEpoch refEpoch;
1628 if (refTime.length()>0) {
1629 refEpoch = epochFromString(refTime, in.getTimeReference());
1630 } else {
1631 refEpoch = in.getEpoch(0);
1632 }
1633 ostringstream oss;
1634 oss << "Aligned at reference Epoch " << formatEpoch(refEpoch) << " in frame " << MFrequency::showType(freqSystem);
1635 pushLog(String(oss));
1636
1637 // Get Reference Position
1638
1639 MPosition refPos = in.getAntennaPosition();
1640
1641 // Create FrequencyAligner Block. One FA for each possible
1642 // source/freqID (perFreqID=True) or source/IF (perFreqID=False)
1643 // combination
1644
1645 PtrBlock<FrequencyAligner<Float>* > a(dDesc.length());
1646 generateFrequencyAligners(a, dDesc, in, nChan, freqSystem, refPos,
1647 refEpoch, perFreqID);
1648
1649 // Generate and fill output Frequency Table. WHen perFreqID=True,
1650 // there is one output FreqID for each entry in the SDDataDesc
1651 // table. However, in perFreqID=False mode, there may be some
1652 // degeneracy, so we need a little translation map
1653
1654 SDFrequencyTable freqTabOut = in.getSDFreqTable();
1655 freqTabOut.setLength(0);
1656 Vector<String> units(1);
1657 units = String("Hz");
1658 Bool linear=True;
1659 //
1660 Vector<uInt> ddFQTrans(dDesc.length(),0);
1661 for (uInt i=0; i<dDesc.length(); i++) {
1662
1663 // Get Aligned SC in Hz
1664
1665 SpectralCoordinate sC = a[i]->alignedSpectralCoordinate(linear);
1666 sC.setWorldAxisUnits(units);
1667
1668 // Add FreqID
1669
1670 uInt idx = freqTabOut.addFrequency(sC.referencePixel()[0],
1671 sC.referenceValue()[0],
1672 sC.increment()[0]);
1673 // output FreqID = ddFQTrans(ddIdx)
1674 ddFQTrans(i) = idx;
1675 }
1676
1677 // Interpolation method
1678
1679 InterpolateArray1D<Double,Float>::InterpolationMethod interp;
1680 convertInterpString(interp, methodStr);
1681
1682 // New output Table
1683
1684 //pushLog("Create output table");
1685 SDMemTable* pTabOut = new SDMemTable(in,True);
1686 pTabOut->putSDFreqTable(freqTabOut);
1687
1688 // Loop over rows in Table
1689
1690 Bool extrapolate=False;
1691 const IPosition polChanAxes(2, asap::PolAxis, asap::ChanAxis);
1692 Bool useCachedAbcissa = False;
1693 Bool first = True;
1694 Bool ok;
1695 Vector<Float> yOut;
1696 Vector<Bool> maskOut;
1697 Vector<uInt> freqID(nIF);
1698 uInt ifIdx, faIdx;
1699 Vector<Double> xIn;
1700 //
1701 for (uInt iRow=0; iRow<nRows; ++iRow) {
1702// if (iRow%10==0) {
1703// ostringstream oss;
1704// oss << "Processing row " << iRow;
1705// pushLog(String(oss));
1706// }
1707
1708 // Get EPoch
1709
1710 Quantum<Double> tQ2(times[iRow],DAY);
1711 MVEpoch mv2(tQ2);
1712 MEpoch epoch(mv2, epochRef);
1713
1714 // Get copy of data
1715
1716 const MaskedArray<Float>& mArrIn(in.rowAsMaskedArray(iRow));
1717 Array<Float> values = mArrIn.getArray();
1718 Array<Bool> mask = mArrIn.getMask();
1719
1720 // For each row, the Frequency abcissa will be the same
1721 // regardless of polarization. For all other axes (IF and BEAM)
1722 // the abcissa will change. So we iterate through the data by
1723 // pol-chan planes to mimimize the work. Probably won't work for
1724 // multiple beams at this point.
1725
1726 ArrayIterator<Float> itValuesPlane(values, polChanAxes);
1727 ArrayIterator<Bool> itMaskPlane(mask, polChanAxes);
1728 while (!itValuesPlane.pastEnd()) {
1729
1730 // Find the IF index and then the FA PtrBlock index
1731
1732 const IPosition& pos = itValuesPlane.pos();
1733 ifIdx = pos(asap::IFAxis);
1734 faIdx = ddIdx(iRow,ifIdx);
1735
1736 // Generate abcissa for perIF. Could cache this in a Matrix on
1737 // a per scan basis. Pretty expensive doing it for every row.
1738
1739 if (!perFreqID) {
1740 xIn.resize(nChan);
1741 uInt fqID = dDesc.secID(ddIdx(iRow,ifIdx));
1742 SpectralCoordinate sC = in.getSpectralCoordinate(fqID);
1743 Double w;
1744 for (uInt i=0; i<nChan; i++) {
1745 sC.toWorld(w,Double(i));
1746 xIn[i] = w;
1747 }
1748 }
1749
1750 VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1);
1751 VectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
1752
1753 // Iterate through the plane by vector and align
1754
1755 first = True;
1756 useCachedAbcissa=False;
1757 while (!itValuesVec.pastEnd()) {
1758 if (perFreqID) {
1759 ok = a[faIdx]->align (yOut, maskOut, itValuesVec.vector(),
1760 itMaskVec.vector(), epoch, useCachedAbcissa,
1761 interp, extrapolate);
1762 } else {
1763 ok = a[faIdx]->align (yOut, maskOut, xIn, itValuesVec.vector(),
1764 itMaskVec.vector(), epoch, useCachedAbcissa,
1765 interp, extrapolate);
1766 }
1767 //
1768 itValuesVec.vector() = yOut;
1769 itMaskVec.vector() = maskOut;
1770 //
1771 itValuesVec.next();
1772 itMaskVec.next();
1773 //
1774 if (first) {
1775 useCachedAbcissa = True;
1776 first = False;
1777 }
1778 }
1779 //
1780 itValuesPlane.next();
1781 itMaskPlane.next();
1782 }
1783
1784 // Create SDContainer and put back
1785
1786 SDContainer sc = in.getSDContainer(iRow);
1787 putDataInSDC(sc, values, mask);
1788
1789 // Set output FreqIDs
1790
1791 for (uInt i=0; i<nIF; i++) {
1792 uInt idx = ddIdx(iRow,i); // Index into SDDataDesc table
1793 freqID(i) = ddFQTrans(idx); // FreqID in output FQ table
1794 }
1795 sc.putFreqMap(freqID);
1796 //
1797 pTabOut->putSDContainer(sc);
1798 }
1799
1800 // Now we must set the base and extra frames to the input frame
1801 std::vector<string> info = pTabOut->getCoordInfo();
1802 info[1] = MFrequency::showType(freqSystem); // Conversion frame
1803 info[3] = info[1]; // Base frame
1804 pTabOut->setCoordInfo(info);
1805
1806 // Clean up PointerBlock
1807 for (uInt i=0; i<a.nelements(); i++) delete a[i];
1808
1809 return pTabOut;
1810}
1811
1812
1813SDMemTable* SDMath::frequencySwitch(const SDMemTable& in)
1814{
1815 if (in.nIF() != 2) {
1816 throw(AipsError("nIF != 2 "));
1817 }
1818 Bool clear = True;
1819 SDMemTable* pTabOut = new SDMemTable(in, clear);
1820 const Table& tabIn = in.table();
1821
1822 // Shape of input and output data
1823 const IPosition& shapeIn = in.rowAsMaskedArray(0).shape();
1824
1825 const uInt nRows = in.nRow();
1826 AlwaysAssert(asap::nAxes==4,AipsError);
1827
1828 ROArrayColumn<Float> tSysCol;
1829 Array<Float> tsys;
1830 tSysCol.attach(tabIn,"TSYS");
1831
1832 for (uInt iRow=0; iRow<nRows; iRow++) {
1833 // Get data for this row
1834 MaskedArray<Float> marr(in.rowAsMaskedArray(iRow));
1835 tSysCol.get(iRow, tsys);
1836
1837 // whole Array for IF 0
1838 IPosition start(asap::nAxes,0);
1839 IPosition end = shapeIn-1;
1840 end(asap::IFAxis) = 0;
1841
1842 MaskedArray<Float> on = marr(start,end);
1843 Array<Float> ton = tsys(start,end);
1844 // Make a copy as "src" is a refrence which is manipulated.
1845 // oncopy is needed for the inverse quotient
1846 MaskedArray<Float> oncopy = on.copy();
1847
1848 // whole Array for IF 1
1849 start(asap::IFAxis) = 1;
1850 end(asap::IFAxis) = 1;
1851
1852 MaskedArray<Float> off = marr(start,end);
1853 Array<Float> toff = tsys(start,end);
1854
1855 on /= off; on -= 1.0f;
1856 on *= ton;
1857 off /= oncopy; off -= 1.0f;
1858 off *= toff;
1859
1860 SDContainer sc = in.getSDContainer(iRow);
1861 putDataInSDC(sc, marr.getArray(), marr.getMask());
1862 pTabOut->putSDContainer(sc);
1863 }
1864 return pTabOut;
1865}
1866
1867void SDMath::fillSDC(SDContainer& sc,
1868 const Array<Bool>& mask,
1869 const Array<Float>& data,
1870 const Array<Float>& tSys,
1871 Int scanID, Double timeStamp,
1872 Double interval, const String& sourceName,
1873 const Vector<uInt>& freqID)
1874{
1875// Data and mask
1876
1877 putDataInSDC(sc, data, mask);
1878
1879// TSys
1880
1881 sc.putTsys(tSys);
1882
1883// Time things
1884
1885 sc.timestamp = timeStamp;
1886 sc.interval = interval;
1887 sc.scanid = scanID;
1888//
1889 sc.sourcename = sourceName;
1890 sc.putFreqMap(freqID);
1891}
1892
1893void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum,
1894 MaskedArray<Float>& sum, Array<Float>& sumSq,
1895 Array<Float>& nPts, Array<Float>& tSysSum,
1896 Array<Float>& tSysSqSum,
1897 const Array<Float>& tSys, const Array<Float>& nInc,
1898 const Vector<Bool>& mask, Double time, Double interval,
1899 const std::vector<CountedPtr<SDMemTable> >& in,
1900 uInt iTab, uInt iRow, uInt axis,
1901 uInt nAxesSub, Bool useMask,
1902 WeightType wtType)
1903{
1904
1905// Get data
1906
1907 MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
1908 Array<Float>& valuesIn = dataIn.getRWArray(); // writable reference
1909 const Array<Bool>& maskIn = dataIn.getMask(); // RO reference
1910//
1911 if (wtType==NONE) {
1912 const MaskedArray<Float> n(nInc,dataIn.getMask());
1913 nPts += n; // Only accumulates where mask==T
1914 } else if (wtType==TINT) {
1915
1916// We are weighting the data by integration time.
1917
1918 valuesIn *= Float(interval);
1919
1920 } else if (wtType==VAR) {
1921
1922// We are going to average the data, weighted by the noise for each pol, beam and IF.
1923// So therefore we need to iterate through by spectrum (axis 3)
1924
1925 VectorIterator<Float> itData(valuesIn, axis);
1926 ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
1927 Float fac = 1.0;
1928 IPosition pos(nAxesSub,0);
1929//
1930 while (!itData.pastEnd()) {
1931
1932// Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
1933
1934 if (useMask) {
1935 MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
1936 fac = 1.0/variance(tmp);
1937 } else {
1938 MaskedArray<Float> tmp(itData.vector(),itMask.vector());
1939 fac = 1.0/variance(tmp);
1940 }
1941
1942// Scale data
1943
1944 itData.vector() *= fac; // Writes back into 'dataIn'
1945//
1946// Accumulate variance per if/pol/beam averaged over spectrum
1947// This method to get pos2 from itData.pos() is only valid
1948// because the spectral axis is the last one (so we can just
1949// copy the first nAXesSub positions out)
1950
1951 pos = itData.pos().getFirst(nAxesSub);
1952 sumSq(pos) += fac;
1953//
1954 itData.next();
1955 itMask.next();
1956 }
1957 } else if (wtType==TSYS || wtType==TINTSYS) {
1958
1959// We are going to average the data, weighted by 1/Tsys**2 for each pol, beam and IF.
1960// So therefore we need to iterate through by spectrum (axis 3). Although
1961// Tsys is stored as a vector of length nChan, the values are replicated.
1962// We will take a short cut and just use the value from the first channel
1963// for now.
1964//
1965 VectorIterator<Float> itData(valuesIn, axis);
1966 ReadOnlyVectorIterator<Float> itTSys(tSys, axis);
1967 IPosition pos(nAxesSub,0);
1968//
1969 Float fac = 1.0;
1970 if (wtType==TINTSYS) fac *= interval;
1971 while (!itData.pastEnd()) {
1972 Float t = itTSys.vector()[0];
1973 fac *= 1.0/t/t;
1974
1975// Scale data
1976
1977 itData.vector() *= fac; // Writes back into 'dataIn'
1978//
1979// Accumulate Tsys per if/pol/beam averaged over spectrum
1980// This method to get pos2 from itData.pos() is only valid
1981// because the spectral axis is the last one (so we can just
1982// copy the first nAXesSub positions out)
1983
1984 pos = itData.pos().getFirst(nAxesSub);
1985 tSysSqSum(pos) += fac;
1986//
1987 itData.next();
1988 itTSys.next();
1989 }
1990 }
1991
1992// Accumulate sum of (possibly scaled) data
1993
1994 sum += dataIn;
1995
1996// Accumulate Tsys, time, and interval
1997
1998 tSysSum += tSys;
1999 timeSum += time;
2000 intSum += interval;
2001 nAccum += 1;
2002}
2003
2004
2005void SDMath::normalize(MaskedArray<Float>& sum,
2006 const Array<Float>& sumSq,
2007 const Array<Float>& tSysSqSum,
2008 const Array<Float>& nPts,
2009 Double intSum,
2010 WeightType wtType, Int axis,
2011 Int nAxesSub)
2012{
2013 IPosition pos2(nAxesSub,0);
2014//
2015 if (wtType==NONE) {
2016
2017// We just average by the number of points accumulated.
2018// We need to make a MA out of nPts so that no divide by
2019// zeros occur
2020
2021 MaskedArray<Float> t(nPts, (nPts>Float(0.0)));
2022 sum /= t;
2023 } else if (wtType==TINT) {
2024
2025// Average by sum of Tint
2026
2027 sum /= Float(intSum);
2028 } else if (wtType==VAR) {
2029
2030// Normalize each spectrum by sum(1/var) where the variance
2031// is worked out for each spectrum
2032
2033 Array<Float>& data = sum.getRWArray();
2034 VectorIterator<Float> itData(data, axis);
2035 while (!itData.pastEnd()) {
2036 pos2 = itData.pos().getFirst(nAxesSub);
2037 itData.vector() /= sumSq(pos2);
2038 itData.next();
2039 }
2040 } else if (wtType==TSYS || wtType==TINTSYS) {
2041
2042// Normalize each spectrum by sum(1/Tsys**2) (TSYS) or
2043// sum(Tint/Tsys**2) (TINTSYS) where the pseudo
2044// replication over channel for Tsys has been dropped.
2045
2046 Array<Float>& data = sum.getRWArray();
2047 VectorIterator<Float> itData(data, axis);
2048 while (!itData.pastEnd()) {
2049 pos2 = itData.pos().getFirst(nAxesSub);
2050 itData.vector() /= tSysSqSum(pos2);
2051 itData.next();
2052 }
2053 }
2054}
2055
2056
2057
2058
2059void SDMath::setCursorSlice (IPosition& start, IPosition& end, Bool doAll, const SDMemTable& in) const
2060{
2061 const uInt nDim = asap::nAxes;
2062 DebugAssert(nDim==4,AipsError);
2063//
2064 start.resize(nDim);
2065 end.resize(nDim);
2066 if (doAll) {
2067 start = 0;
2068 end(0) = in.nBeam()-1;
2069 end(1) = in.nIF()-1;
2070 end(2) = in.nPol()-1;
2071 end(3) = in.nChan()-1;
2072 } else {
2073 start(0) = in.getBeam();
2074 end(0) = start(0);
2075//
2076 start(1) = in.getIF();
2077 end(1) = start(1);
2078//
2079 start(2) = in.getPol();
2080 end(2) = start(2);
2081//
2082 start(3) = 0;
2083 end(3) = in.nChan()-1;
2084 }
2085}
2086
2087
2088void SDMath::convertWeightString(WeightType& wtType, const String& weightStr,
2089 Bool listType)
2090{
2091 String tStr(weightStr);
2092 tStr.upcase();
2093 String msg;
2094 if (tStr.contains(String("NONE"))) {
2095 wtType = NONE;
2096 msg = String("Weighting type selected : None");
2097 } else if (tStr.contains(String("VAR"))) {
2098 wtType = VAR;
2099 msg = String("Weighting type selected : Variance");
2100 } else if (tStr.contains(String("TINTSYS"))) {
2101 wtType = TINTSYS;
2102 msg = String("Weighting type selected : Tint&Tsys");
2103 } else if (tStr.contains(String("TINT"))) {
2104 wtType = TINT;
2105 msg = String("Weighting type selected : Tint");
2106 } else if (tStr.contains(String("TSYS"))) {
2107 wtType = TSYS;
2108 msg = String("Weighting type selected : Tsys");
2109 } else {
2110 msg = String("Weighting type selected : None");
2111 throw(AipsError("Unrecognized weighting type"));
2112 }
2113//
2114 if (listType) pushLog(msg);
2115}
2116
2117
2118void SDMath::convertInterpString(casa::InterpolateArray1D<Double,Float>::InterpolationMethod& type,
2119 const casa::String& interp)
2120{
2121 String tStr(interp);
2122 tStr.upcase();
2123 if (tStr.contains(String("NEAR"))) {
2124 type = InterpolateArray1D<Double,Float>::nearestNeighbour;
2125 } else if (tStr.contains(String("LIN"))) {
2126 type = InterpolateArray1D<Double,Float>::linear;
2127 } else if (tStr.contains(String("CUB"))) {
2128 type = InterpolateArray1D<Double,Float>::cubic;
2129 } else if (tStr.contains(String("SPL"))) {
2130 type = InterpolateArray1D<Double,Float>::spline;
2131 } else {
2132 throw(AipsError("Unrecognized interpolation type"));
2133 }
2134}
2135
2136void SDMath::putDataInSDC(SDContainer& sc, const Array<Float>& data,
2137 const Array<Bool>& mask)
2138{
2139 sc.putSpectrum(data);
2140//
2141 Array<uChar> outflags(data.shape());
2142 convertArray(outflags,!mask);
2143 sc.putFlags(outflags);
2144}
2145
2146Table SDMath::readAsciiFile(const String& fileName) const
2147{
2148 String formatString;
2149 Table tbl = readAsciiTable (formatString, Table::Memory, fileName, "", "", False);
2150 return tbl;
2151}
2152
2153
2154
2155void SDMath::scaleFromAsciiTable(SDMemTable* pTabOut,
2156 const SDMemTable& in, const String& fileName,
2157 const String& col0, const String& col1,
2158 const String& methodStr, Bool doAll,
2159 const Vector<Float>& xOut, Bool doTSys)
2160{
2161
2162// Read gain-elevation ascii file data into a Table.
2163
2164 Table geTable = readAsciiFile (fileName);
2165//
2166 scaleFromTable (pTabOut, in, geTable, col0, col1, methodStr, doAll, xOut, doTSys);
2167}
2168
2169void SDMath::scaleFromTable(SDMemTable* pTabOut, const SDMemTable& in,
2170 const Table& tTable, const String& col0,
2171 const String& col1,
2172 const String& methodStr, Bool doAll,
2173 const Vector<Float>& xOut, Bool doTsys)
2174{
2175
2176// Get data from Table
2177
2178 ROScalarColumn<Float> geElCol(tTable, col0);
2179 ROScalarColumn<Float> geFacCol(tTable, col1);
2180 Vector<Float> xIn = geElCol.getColumn();
2181 Vector<Float> yIn = geFacCol.getColumn();
2182 Vector<Bool> maskIn(xIn.nelements(),True);
2183
2184// Interpolate (and extrapolate) with desired method
2185
2186 InterpolateArray1D<Double,Float>::InterpolationMethod method;
2187 convertInterpString(method, methodStr);
2188 Int intMethod(method);
2189//
2190 Vector<Float> yOut;
2191 Vector<Bool> maskOut;
2192 InterpolateArray1D<Float,Float>::interpolate(yOut, maskOut, xOut,
2193 xIn, yIn, maskIn, intMethod,
2194 True, True);
2195// Apply
2196
2197 scaleByVector(pTabOut, in, doAll, Float(1.0)/yOut, doTsys);
2198}
2199
2200
2201void SDMath::scaleByVector(SDMemTable* pTabOut, const SDMemTable& in,
2202 Bool doAll, const Vector<Float>& factor,
2203 Bool doTSys)
2204{
2205
2206// Set up data slice
2207
2208 IPosition start, end;
2209 setCursorSlice (start, end, doAll, in);
2210
2211// Get Tsys column
2212
2213 const Table& tIn = in.table();
2214 ArrayColumn<Float> tSysCol(tIn, "TSYS");
2215 Array<Float> tSys;
2216
2217// Loop over rows and apply correction factor
2218
2219 const uInt axis = asap::ChanAxis;
2220 for (uInt i=0; i < in.nRow(); ++i) {
2221
2222// Get data
2223
2224 MaskedArray<Float> dataIn(in.rowAsMaskedArray(i));
2225 MaskedArray<Float> dataIn2 = dataIn(start,end); // reference to dataIn
2226//
2227 if (doTSys) {
2228 tSysCol.get(i, tSys);
2229 Array<Float> tSys2 = tSys(start,end) * factor[i];
2230 tSysCol.put(i, tSys);
2231 }
2232
2233// Apply factor
2234
2235 dataIn2 *= factor[i];
2236
2237// Write out
2238
2239 SDContainer sc = in.getSDContainer(i);
2240 putDataInSDC(sc, dataIn.getArray(), dataIn.getMask());
2241//
2242 pTabOut->putSDContainer(sc);
2243 }
2244}
2245
2246
2247
2248
2249void SDMath::generateDataDescTable (Matrix<uInt>& ddIdx,
2250 SDDataDesc& dDesc,
2251 uInt nIF,
2252 const SDMemTable& in,
2253 const Table& tabIn,
2254 const ROScalarColumn<String>& srcCol,
2255 const ROArrayColumn<uInt>& fqIDCol,
2256 Bool perFreqID)
2257{
2258 const uInt nRows = tabIn.nrow();
2259 ddIdx.resize(nRows,nIF);
2260//
2261 String srcName;
2262 Vector<uInt> freqIDs;
2263 for (uInt iRow=0; iRow<nRows; iRow++) {
2264 srcCol.get(iRow, srcName);
2265 fqIDCol.get(iRow, freqIDs);
2266 const MDirection& dir = in.getDirection(iRow);
2267//
2268 if (perFreqID) {
2269
2270// One entry per source/freqID pair
2271
2272 for (uInt iIF=0; iIF<nIF; iIF++) {
2273 ddIdx(iRow,iIF) = dDesc.addEntry(srcName, freqIDs[iIF], dir, 0);
2274 }
2275 } else {
2276
2277// One entry per source/IF pair. Hang onto the FreqID as well
2278
2279 for (uInt iIF=0; iIF<nIF; iIF++) {
2280 ddIdx(iRow,iIF) = dDesc.addEntry(srcName, iIF, dir, freqIDs[iIF]);
2281 }
2282 }
2283 }
2284}
2285
2286
2287
2288
2289
2290MEpoch SDMath::epochFromString(const String& str, MEpoch::Types timeRef)
2291{
2292 Quantum<Double> qt;
2293 if (MVTime::read(qt,str)) {
2294 MVEpoch mv(qt);
2295 MEpoch me(mv, timeRef);
2296 return me;
2297 } else {
2298 throw(AipsError("Invalid format for Epoch string"));
2299 }
2300}
2301
2302
2303String SDMath::formatEpoch(const MEpoch& epoch) const
2304{
2305 MVTime mvt(epoch.getValue());
2306 return mvt.string(MVTime::YMD) + String(" (") + epoch.getRefString() + String(")");
2307}
2308
2309
2310
2311void SDMath::generateFrequencyAligners(PtrBlock<FrequencyAligner<Float>* >& a,
2312 const SDDataDesc& dDesc,
2313 const SDMemTable& in, uInt nChan,
2314 MFrequency::Types system,
2315 const MPosition& refPos,
2316 const MEpoch& refEpoch,
2317 Bool perFreqID)
2318{
2319 for (uInt i=0; i<dDesc.length(); i++) {
2320 uInt ID = dDesc.ID(i);
2321 uInt secID = dDesc.secID(i);
2322 const MDirection& refDir = dDesc.secDir(i);
2323//
2324 if (perFreqID) {
2325
2326// One aligner per source/FreqID pair.
2327
2328 SpectralCoordinate sC = in.getSpectralCoordinate(ID);
2329 a[i] = new FrequencyAligner<Float>(sC, nChan, refEpoch, refDir, refPos, system);
2330 } else {
2331
2332// One aligner per source/IF pair. But we still need the FreqID to
2333// get the right SC. Hence the messing about with the secondary ID
2334
2335 SpectralCoordinate sC = in.getSpectralCoordinate(secID);
2336 a[i] = new FrequencyAligner<Float>(sC, nChan, refEpoch, refDir, refPos, system);
2337 }
2338 }
2339}
2340
2341Vector<uInt> SDMath::getRowRange(const SDMemTable& in) const
2342{
2343 Vector<uInt> range(2);
2344 range[0] = 0;
2345 range[1] = in.nRow()-1;
2346 return range;
2347}
2348
2349
2350Bool SDMath::rowInRange(uInt i, const Vector<uInt>& range) const
2351{
2352 return (i>=range[0] && i<=range[1]);
2353}
2354
Note: See TracBrowser for help on using the repository browser.