source: trunk/src/SDMath.cc@ 307

Last change on this file since 307 was 304, checked in by kil064, 20 years ago

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