source: trunk/src/SDMath.cc@ 404

Last change on this file since 404 was 397, checked in by kil064, 20 years ago

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