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

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

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