source: trunk/src/SDMath.cc@ 791

Last change on this file since 791 was 780, checked in by mar637, 19 years ago

merge from Release12

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