source: trunk/src/SDMath.cc@ 298

Last change on this file since 298 was 294, checked in by kil064, 20 years ago

move some code into function generateVelocityAligners
add arg. doTSys to binaryOperate and unaryOperate
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 45.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//
379 return CountedPtr<SDMemTable>(pTabOut);
380}
381
382
383
384CountedPtr<SDMemTable> SDMath::binaryOperate (const CountedPtr<SDMemTable>& left,
385 const CountedPtr<SDMemTable>& right,
386 const String& op, Bool preserve,
387 Bool doTSys) const
388{
389
390// Check operator
391
392 String op2(op);
393 op2.upcase();
394 uInt what = 0;
395 if (op2=="ADD") {
396 what = 0;
397 } else if (op2=="SUB") {
398 what = 1;
399 } else if (op2=="MUL") {
400 what = 2;
401 } else if (op2=="DIV") {
402 what = 3;
403 } else if (op2=="QUOTIENT") {
404 what = 4;
405 doTSys = True;
406 } else {
407 throw( AipsError("Unrecognized operation"));
408 }
409
410// Check rows
411
412 const uInt nRowLeft = left->nRow();
413 const uInt nRowRight = right->nRow();
414 Bool ok = (nRowRight==1&&nRowLeft>0) ||
415 (nRowRight>=1&&nRowLeft==nRowRight);
416 if (!ok) {
417 throw (AipsError("The right Scan Table can have one row or the same number of rows as the left Scan Table"));
418 }
419
420// Input Tables
421
422 const Table& tLeft = left->table();
423 const Table& tRight = right->table();
424
425// TSys columns
426
427 ROArrayColumn<Float> tSysLeftCol, tSysRightCol;
428 if (doTSys) {
429 tSysLeftCol.attach(tLeft, "TSYS");
430 tSysRightCol.attach(tRight, "TSYS");
431 }
432
433// First row for right
434
435 Array<Float> tSysLeftArr, tSysRightArr;
436 if (doTSys) tSysRightCol.get(0, tSysRightArr);
437 MaskedArray<Float>* pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(0));
438 IPosition shpRight = pMRight->shape();
439
440// Output Table cloned from left
441
442 SDMemTable* pTabOut = new SDMemTable(*left, True);
443
444// Loop over rows
445
446 for (uInt i=0; i<nRowLeft; i++) {
447
448// Get data
449
450 MaskedArray<Float> mLeft(left->rowAsMaskedArray(i));
451 IPosition shpLeft = mLeft.shape();
452 if (doTSys) tSysLeftCol.get(i, tSysLeftArr);
453//
454 if (nRowRight>1) {
455 delete pMRight;
456 pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(i));
457 shpRight = pMRight->shape();
458 if (doTSys) tSysRightCol.get(i, tSysRightArr);
459 }
460//
461 if (!shpRight.isEqual(shpLeft)) {
462 throw(AipsError("left and right scan tables are not conformant"));
463 }
464 if (doTSys) {
465 if (!tSysRightArr.shape().isEqual(tSysRightArr.shape())) {
466 throw(AipsError("left and right Tsys data are not conformant"));
467 }
468 if (!shpRight.isEqual(tSysRightArr.shape())) {
469 throw(AipsError("left and right scan tables are not conformant"));
470 }
471 }
472
473// Make container
474
475 SDContainer sc = left->getSDContainer(i);
476
477// Operate on data and TSys
478
479 if (what==0) {
480 MaskedArray<Float> tmp = mLeft + *pMRight;
481 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
482 if (doTSys) sc.putTsys(tSysLeftArr+tSysRightArr);
483 } else if (what==1) {
484 MaskedArray<Float> tmp = mLeft - *pMRight;
485 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
486 if (doTSys) sc.putTsys(tSysLeftArr-tSysRightArr);
487 } else if (what==2) {
488 MaskedArray<Float> tmp = mLeft * *pMRight;
489 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
490 if (doTSys) sc.putTsys(tSysLeftArr*tSysRightArr);
491 } else if (what==3) {
492 MaskedArray<Float> tmp = mLeft / *pMRight;
493 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
494 if (doTSys) sc.putTsys(tSysLeftArr/tSysRightArr);
495 } else if (what==4) {
496 if (preserve) {
497 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) - tSysRightArr;
498 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
499 } else {
500 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) - tSysLeftArr;
501 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
502 }
503 sc.putTsys(tSysRightArr);
504 }
505
506// Put new row in output Table
507
508 pTabOut->putSDContainer(sc);
509 }
510 if (pMRight) delete pMRight;
511//
512 return CountedPtr<SDMemTable>(pTabOut);
513}
514
515
516
517std::vector<float> SDMath::statistic(const CountedPtr<SDMemTable>& in,
518 const Vector<Bool>& mask,
519 const String& which, Int row) const
520//
521// Perhaps iteration over pol/beam/if should be in here
522// and inside the nrow iteration ?
523//
524{
525 const uInt nRow = in->nRow();
526
527// Specify cursor location
528
529 IPosition start, end;
530 getCursorLocation(start, end, *in);
531
532// Loop over rows
533
534 const uInt nEl = mask.nelements();
535 uInt iStart = 0;
536 uInt iEnd = in->nRow()-1;
537//
538 if (row>=0) {
539 iStart = row;
540 iEnd = row;
541 }
542//
543 std::vector<float> result(iEnd-iStart+1);
544 for (uInt ii=iStart; ii <= iEnd; ++ii) {
545
546// Get row and deconstruct
547
548 MaskedArray<Float> marr(in->rowAsMaskedArray(ii));
549 Array<Float> arr = marr.getArray();
550 Array<Bool> barr = marr.getMask();
551
552// Access desired piece of data
553
554 Array<Float> v((arr(start,end)).nonDegenerate());
555 Array<Bool> m((barr(start,end)).nonDegenerate());
556
557// Apply OTF mask
558
559 MaskedArray<Float> tmp;
560 if (m.nelements()==nEl) {
561 tmp.setData(v,m&&mask);
562 } else {
563 tmp.setData(v,m);
564 }
565
566// Get statistic
567
568 result[ii-iStart] = mathutil::statistics(which, tmp);
569 }
570//
571 return result;
572}
573
574
575SDMemTable* SDMath::bin(const SDMemTable& in, Int width) const
576{
577 SDHeader sh = in.getSDHeader();
578 SDMemTable* pTabOut = new SDMemTable(in, True);
579
580// Bin up SpectralCoordinates
581
582 IPosition factors(1);
583 factors(0) = width;
584 for (uInt j=0; j<in.nCoordinates(); ++j) {
585 CoordinateSystem cSys;
586 cSys.addCoordinate(in.getSpectralCoordinate(j));
587 CoordinateSystem cSysBin =
588 CoordinateUtil::makeBinnedCoordinateSystem(factors, cSys, False);
589//
590 SpectralCoordinate sCBin = cSysBin.spectralCoordinate(0);
591 pTabOut->setCoordinate(sCBin, j);
592 }
593
594// Use RebinLattice to find shape
595
596 IPosition shapeIn(1,sh.nchan);
597 IPosition shapeOut = RebinLattice<Float>::rebinShape(shapeIn, factors);
598 sh.nchan = shapeOut(0);
599 pTabOut->putSDHeader(sh);
600
601
602// Loop over rows and bin along channel axis
603
604 for (uInt i=0; i < in.nRow(); ++i) {
605 SDContainer sc = in.getSDContainer(i);
606//
607 Array<Float> tSys(sc.getTsys()); // Get it out before sc changes shape
608
609// Bin up spectrum
610
611 MaskedArray<Float> marr(in.rowAsMaskedArray(i));
612 MaskedArray<Float> marrout;
613 LatticeUtilities::bin(marrout, marr, asap::ChanAxis, width);
614
615// Put back the binned data and flags
616
617 IPosition ip2 = marrout.shape();
618 sc.resize(ip2);
619//
620 putDataInSDC(sc, marrout.getArray(), marrout.getMask());
621
622// Bin up Tsys.
623
624 Array<Bool> allGood(tSys.shape(),True);
625 MaskedArray<Float> tSysIn(tSys, allGood, True);
626//
627 MaskedArray<Float> tSysOut;
628 LatticeUtilities::bin(tSysOut, tSysIn, asap::ChanAxis, width);
629 sc.putTsys(tSysOut.getArray());
630//
631 pTabOut->putSDContainer(sc);
632 }
633 return pTabOut;
634}
635
636SDMemTable* SDMath::unaryOperate(const SDMemTable& in, Float val, Bool doAll,
637 uInt what, Bool doTSys) const
638//
639// what = 0 Multiply
640// 1 Add
641{
642 SDMemTable* pOut = new SDMemTable(in,False);
643 const Table& tOut = pOut->table();
644 ArrayColumn<Float> specCol(tOut,"SPECTRA");
645 ArrayColumn<Float> tSysCol(tOut,"TSYS");
646 Array<Float> tSysArr;
647//
648 if (doAll) {
649 for (uInt i=0; i < tOut.nrow(); i++) {
650
651// Modify data
652
653 MaskedArray<Float> dataIn(pOut->rowAsMaskedArray(i));
654 if (what==0) {
655 dataIn *= val;
656 } else if (what==1) {
657 dataIn += val;
658 }
659 specCol.put(i, dataIn.getArray());
660
661// Modify Tsys
662
663 if (doTSys) {
664 tSysCol.get(i, tSysArr);
665 if (what==0) {
666 tSysArr *= val;
667 } else if (what==1) {
668 tSysArr += val;
669 }
670 tSysCol.put(i, tSysArr);
671 }
672 }
673 } else {
674
675// Get cursor location
676
677 IPosition start, end;
678 getCursorLocation(start, end, in);
679//
680 for (uInt i=0; i < tOut.nrow(); i++) {
681
682// Modify data
683
684 MaskedArray<Float> dataIn(pOut->rowAsMaskedArray(i));
685 MaskedArray<Float> dataIn2 = dataIn(start,end); // Reference
686 if (what==0) {
687 dataIn2 *= val;
688 } else if (what==1) {
689 dataIn2 += val;
690 }
691 specCol.put(i, dataIn.getArray());
692
693// Modify Tsys
694
695 if (doTSys) {
696 tSysCol.get(i, tSysArr);
697 Array<Float> tSysArr2 = tSysArr(start,end); // Reference
698 if (what==0) {
699 tSysArr2 *= val;
700 } else if (what==1) {
701 tSysArr2 += val;
702 }
703 tSysCol.put(i, tSysArr);
704 }
705 }
706 }
707//
708 return pOut;
709}
710
711
712
713SDMemTable* SDMath::averagePol(const SDMemTable& in, const Vector<Bool>& mask) const
714//
715// Average all polarizations together, weighted by variance
716//
717{
718// WeightType wtType = NONE;
719// convertWeightString(wtType, weight);
720
721 const uInt nRows = in.nRow();
722
723// Create output Table and reshape number of polarizations
724
725 Bool clear=True;
726 SDMemTable* pTabOut = new SDMemTable(in, clear);
727 SDHeader header = pTabOut->getSDHeader();
728 header.npol = 1;
729 pTabOut->putSDHeader(header);
730
731// Shape of input and output data
732
733 const IPosition& shapeIn = in.rowAsMaskedArray(0u, False).shape();
734 IPosition shapeOut(shapeIn);
735 shapeOut(asap::PolAxis) = 1; // Average all polarizations
736//
737 const uInt nChan = shapeIn(asap::ChanAxis);
738 const IPosition vecShapeOut(4,1,1,1,nChan); // A multi-dim form of a Vector shape
739 IPosition start(4), end(4);
740
741// Output arrays
742
743 Array<Float> outData(shapeOut, 0.0);
744 Array<Bool> outMask(shapeOut, True);
745 const IPosition axes(2, asap::PolAxis, asap::ChanAxis); // pol-channel plane
746//
747 const Bool useMask = (mask.nelements() == shapeIn(asap::ChanAxis));
748
749// Loop over rows
750
751 for (uInt iRow=0; iRow<nRows; iRow++) {
752
753// Get data for this row
754
755 MaskedArray<Float> marr(in.rowAsMaskedArray(iRow));
756 Array<Float>& arr = marr.getRWArray();
757 const Array<Bool>& barr = marr.getMask();
758
759// Make iterators to iterate by pol-channel planes
760
761 ReadOnlyArrayIterator<Float> itDataPlane(arr, axes);
762 ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes);
763
764// Accumulations
765
766 Float fac = 1.0;
767 Vector<Float> vecSum(nChan,0.0);
768
769// Iterate through data by pol-channel planes
770
771 while (!itDataPlane.pastEnd()) {
772
773// Iterate through plane by polarization and accumulate Vectors
774
775 Vector<Float> t1(nChan); t1 = 0.0;
776 Vector<Bool> t2(nChan); t2 = True;
777 MaskedArray<Float> vecSum(t1,t2);
778 Float varSum = 0.0;
779 {
780 ReadOnlyVectorIterator<Float> itDataVec(itDataPlane.array(), 1);
781 ReadOnlyVectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
782 while (!itDataVec.pastEnd()) {
783
784// Create MA of data & mask (optionally including OTF mask) and get variance
785
786 if (useMask) {
787 const MaskedArray<Float> spec(itDataVec.vector(),mask&&itMaskVec.vector());
788 fac = 1.0 / variance(spec);
789 } else {
790 const MaskedArray<Float> spec(itDataVec.vector(),itMaskVec.vector());
791 fac = 1.0 / variance(spec);
792 }
793
794// Normalize spectrum (without OTF mask) and accumulate
795
796 const MaskedArray<Float> spec(fac*itDataVec.vector(), itMaskVec.vector());
797 vecSum += spec;
798 varSum += fac;
799
800// Next
801
802 itDataVec.next();
803 itMaskVec.next();
804 }
805 }
806
807// Normalize summed spectrum
808
809 vecSum /= varSum;
810
811// FInd position in input data array. We are iterating by pol-channel
812// plane so all that will change is beam and IF and that's what we want.
813
814 IPosition pos = itDataPlane.pos();
815
816// Write out data. This is a bit messy. We have to reform the Vector
817// accumulator into an Array of shape (1,1,1,nChan)
818
819 start = pos;
820 end = pos;
821 end(asap::ChanAxis) = nChan-1;
822 outData(start,end) = vecSum.getArray().reform(vecShapeOut);
823 outMask(start,end) = vecSum.getMask().reform(vecShapeOut);
824
825// Step to next beam/IF combination
826
827 itDataPlane.next();
828 itMaskPlane.next();
829 }
830
831// Generate output container and write it to output table
832
833 SDContainer sc = in.getSDContainer();
834 sc.resize(shapeOut);
835//
836 putDataInSDC(sc, outData, outMask);
837 pTabOut->putSDContainer(sc);
838 }
839//
840 return pTabOut;
841}
842
843
844SDMemTable* SDMath::smooth(const SDMemTable& in,
845 const casa::String& kernelType,
846 casa::Float width, Bool doAll) const
847{
848
849// Number of channels
850
851 const uInt chanAxis = asap::ChanAxis; // Spectral axis
852 SDHeader sh = in.getSDHeader();
853 const uInt nChan = sh.nchan;
854
855// Generate Kernel
856
857 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernelType);
858 Vector<Float> kernel = VectorKernel::make(type, width, nChan, True, False);
859
860// Generate Convolver
861
862 IPosition shape(1,nChan);
863 Convolver<Float> conv(kernel, shape);
864
865// New Table
866
867 SDMemTable* pTabOut = new SDMemTable(in,True);
868
869// Get cursor location
870
871 IPosition start, end;
872 getCursorLocation(start, end, in);
873//
874 IPosition shapeOut(4,1);
875
876// Output Vectors
877
878 Vector<Float> valuesOut(nChan);
879 Vector<Bool> maskOut(nChan);
880
881// Loop over rows in Table
882
883 for (uInt ri=0; ri < in.nRow(); ++ri) {
884
885// Get copy of data
886
887 const MaskedArray<Float>& dataIn(in.rowAsMaskedArray(ri));
888 AlwaysAssert(dataIn.shape()(asap::ChanAxis)==nChan, AipsError);
889//
890 Array<Float> valuesIn = dataIn.getArray();
891 Array<Bool> maskIn = dataIn.getMask();
892
893// Branch depending on whether we smooth all locations or just
894// those pointed at by the current selection cursor
895
896 if (doAll) {
897 uInt axis = asap::ChanAxis;
898 VectorIterator<Float> itValues(valuesIn, axis);
899 VectorIterator<Bool> itMask(maskIn, axis);
900 while (!itValues.pastEnd()) {
901
902// Smooth
903 if (kernelType==VectorKernel::HANNING) {
904 mathutil::hanning(valuesOut, maskOut, itValues.vector(), itMask.vector());
905 itMask.vector() = maskOut;
906 } else {
907 mathutil::replaceMaskByZero(itValues.vector(), itMask.vector());
908 conv.linearConv(valuesOut, itValues.vector());
909 }
910//
911 itValues.vector() = valuesOut;
912//
913 itValues.next();
914 itMask.next();
915 }
916 } else {
917
918// Set multi-dim Vector shape
919
920 shapeOut(asap::ChanAxis) = valuesIn.shape()(chanAxis);
921
922// Stuff about with shapes so that we don't have conformance run-time errors
923
924 Vector<Float> valuesIn2 = valuesIn(start,end).nonDegenerate();
925 Vector<Bool> maskIn2 = maskIn(start,end).nonDegenerate();
926
927// Smooth
928
929 if (kernelType==VectorKernel::HANNING) {
930 mathutil::hanning(valuesOut, maskOut, valuesIn2, maskIn2);
931 maskIn(start,end) = maskOut.reform(shapeOut);
932 } else {
933 mathutil::replaceMaskByZero(valuesIn2, maskIn2);
934 conv.linearConv(valuesOut, valuesIn2);
935 }
936//
937 valuesIn(start,end) = valuesOut.reform(shapeOut);
938 }
939
940// Create and put back
941
942 SDContainer sc = in.getSDContainer(ri);
943 putDataInSDC(sc, valuesIn, maskIn);
944//
945 pTabOut->putSDContainer(sc);
946 }
947//
948 return pTabOut;
949}
950
951
952
953SDMemTable* SDMath::convertFlux (const SDMemTable& in, Float a, Float eta, Bool doAll) const
954//
955// As it is, this function could be implemented with 'simpleOperate'
956// However, I anticipate that eventually we will look the conversion
957// values up in a Table and apply them in a frequency dependent way,
958// so I have implemented it fully here
959//
960{
961 SDHeader sh = in.getSDHeader();
962 SDMemTable* pTabOut = new SDMemTable(in, True);
963
964// FInd out how to convert values into Jy and K (e.g. units might be mJy or mK)
965// Also automatically find out what we are converting to according to the
966// flux unit
967
968 Unit fluxUnit(sh.fluxunit);
969 Unit K(String("K"));
970 Unit JY(String("Jy"));
971//
972 Bool toKelvin = True;
973 Double inFac = 1.0;
974 if (fluxUnit==JY) {
975 cerr << "Converting to K" << endl;
976//
977 Quantum<Double> t(1.0,fluxUnit);
978 Quantum<Double> t2 = t.get(JY);
979 inFac = (t2 / t).getValue();
980//
981 toKelvin = True;
982 sh.fluxunit = "K";
983 } else if (fluxUnit==K) {
984 cerr << "Converting to Jy" << endl;
985//
986 Quantum<Double> t(1.0,fluxUnit);
987 Quantum<Double> t2 = t.get(K);
988 inFac = (t2 / t).getValue();
989//
990 toKelvin = False;
991 sh.fluxunit = "Jy";
992 } else {
993 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
994 }
995 pTabOut->putSDHeader(sh);
996
997// Compute conversion factor. 'a' and 'eta' are really frequency, time and
998// telescope dependent and should be looked// up in a table
999
1000 Float factor = 2.0 * inFac * 1.0e-7 * 1.0e26 *
1001 QC::k.getValue(Unit(String("erg/K"))) / a / eta;
1002 if (toKelvin) {
1003 factor = 1.0 / factor;
1004 }
1005 cerr << "Applying conversion factor = " << factor << endl;
1006
1007// Generate correction vector. Apply same factor regardless
1008// of beam/pol/IF. This will need to change somewhen.
1009
1010 Vector<Float> factors(in.nRow(), factor);
1011
1012// Correct
1013
1014 correctFromVector (pTabOut, in, doAll, factors);
1015//
1016 return pTabOut;
1017}
1018
1019
1020SDMemTable* SDMath::gainElevation (const SDMemTable& in, const Vector<Float>& coeffs,
1021 const String& fileName,
1022 const String& methodStr, Bool doAll) const
1023{
1024
1025// Get header and clone output table
1026
1027 SDHeader sh = in.getSDHeader();
1028 SDMemTable* pTabOut = new SDMemTable(in, True);
1029
1030// Get elevation data from SDMemTable and convert to degrees
1031
1032 const Table& tab = in.table();
1033 ROScalarColumn<Float> elev(tab, "ELEVATION");
1034 Vector<Float> x = elev.getColumn();
1035 x *= Float(180 / C::pi);
1036//
1037 const uInt nC = coeffs.nelements();
1038 if (fileName.length()>0 && nC>0) {
1039 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
1040 }
1041
1042// Correct
1043
1044 if (nC>0 || fileName.length()==0) {
1045
1046// Find instrument
1047
1048 Bool throwIt = True;
1049 Instrument inst = SDMemTable::convertInstrument (sh.antennaname, throwIt);
1050
1051// Set polynomial
1052
1053 Polynomial<Float>* pPoly = 0;
1054 Vector<Float> coeff;
1055 String msg;
1056 if (nC>0) {
1057 pPoly = new Polynomial<Float>(nC);
1058 coeff = coeffs;
1059 msg = String("user");
1060 } else {
1061 if (inst==ATPKSMB) {
1062 } else if (inst==ATPKSHOH) {
1063 } else if (inst==TIDBINBILLA) {
1064 pPoly = new Polynomial<Float>(3);
1065 coeff.resize(3);
1066 coeff(0) = 3.58788e-1;
1067 coeff(1) = 2.87243e-2;
1068 coeff(2) = -3.219093e-4;
1069 } else if (inst==ATMOPRA) {
1070 } else {
1071 }
1072 msg = String("built in");
1073 }
1074//
1075 if (coeff.nelements()>0) {
1076 pPoly->setCoefficients(coeff);
1077 } else {
1078 throw(AipsError("There is no known gain-el polynomial known for this instrument"));
1079 }
1080//
1081 cerr << "Making polynomial correction with " << msg << " coefficients" << endl;
1082 const uInt nRow = in.nRow();
1083 Vector<Float> factor(nRow);
1084 for (uInt i=0; i<nRow; i++) {
1085 factor[i] = (*pPoly)(x[i]);
1086 }
1087 delete pPoly;
1088//
1089 correctFromVector (pTabOut, in, doAll, factor);
1090 } else {
1091
1092// Indicate which columns to read from ascii file
1093
1094 String col0("ELEVATION");
1095 String col1("FACTOR");
1096
1097// Read and correct
1098
1099 cerr << "Making correction from ascii Table" << endl;
1100 correctFromAsciiTable (pTabOut, in, fileName, col0, col1,
1101 methodStr, doAll, x);
1102 }
1103//
1104 return pTabOut;
1105}
1106
1107
1108
1109SDMemTable* SDMath::opacity (const SDMemTable& in, Float tau, Bool doAll) const
1110{
1111
1112// Get header and clone output table
1113
1114 SDHeader sh = in.getSDHeader();
1115 SDMemTable* pTabOut = new SDMemTable(in, True);
1116
1117// Get elevation data from SDMemTable and convert to degrees
1118
1119 const Table& tab = in.table();
1120 ROScalarColumn<Float> elev(tab, "ELEVATION");
1121 Vector<Float> zDist = elev.getColumn();
1122 zDist = Float(C::pi_2) - zDist;
1123
1124// Generate correction factor
1125
1126 const uInt nRow = in.nRow();
1127 Vector<Float> factor(nRow);
1128 Vector<Float> factor2(nRow);
1129 for (uInt i=0; i<nRow; i++) {
1130 factor[i] = exp(tau)/cos(zDist[i]);
1131 }
1132
1133// Correct
1134
1135 correctFromVector (pTabOut, in, doAll, factor);
1136//
1137 return pTabOut;
1138}
1139
1140
1141
1142
1143// 'private' functions
1144
1145SDMemTable* SDMath::velocityAlign (const SDMemTable& in,
1146 MFrequency::Types velSystem,
1147 const String& velUnit,
1148 MDoppler::Types doppler,
1149 const String& refTime) const
1150{
1151// Get Header
1152
1153 SDHeader sh = in.getSDHeader();
1154 const uInt nChan = sh.nchan;
1155 const uInt nRows = in.nRow();
1156
1157// Get Table reference
1158
1159 const Table& tabIn = in.table();
1160
1161// Get Columns from Table
1162
1163 ROScalarColumn<Double> mjdCol(tabIn, "TIME");
1164 ROScalarColumn<String> srcCol(tabIn, "SRCNAME");
1165 ROArrayColumn<uInt> fqIDCol(tabIn, "FREQID");
1166//
1167 Vector<Double> times = mjdCol.getColumn();
1168 Vector<String> srcNames = srcCol.getColumn();
1169 Vector<uInt> freqID;
1170
1171// Generate Source table
1172
1173 Vector<String> srcTab;
1174 Vector<uInt> srcIdx, firstRow;
1175 generateSourceTable (srcTab, srcIdx, firstRow, srcNames);
1176 const uInt nSrcTab = srcTab.nelements();
1177 cerr << "Found " << srcTab.nelements() << " sources to align " << endl;
1178
1179// Get reference Epoch to time of first row or given String
1180
1181 Unit DAY(String("d"));
1182 MEpoch::Ref epochRef(in.getTimeReference());
1183 MEpoch refEpoch;
1184 if (refTime.length()>0) {
1185 refEpoch = epochFromString(refTime, in.getTimeReference());
1186 } else {
1187 refEpoch = in.getEpoch(0);
1188 }
1189 cerr << "Aligning at reference Epoch " << formatEpoch(refEpoch) << endl;
1190
1191// Get Reference Position
1192
1193 MPosition refPos = in.getAntennaPosition();
1194
1195// Get Frequency Table
1196
1197 SDFrequencyTable fTab = in.getSDFreqTable();
1198 const uInt nFreqIDs = fTab.length();
1199
1200// Create VelocityAligner Block. One VA for each possible
1201// source/freqID combination
1202
1203 PtrBlock<VelocityAligner<Float>* > vA(nFreqIDs*nSrcTab);
1204 generateVelocityAligners (vA, in, nChan, nFreqIDs, nSrcTab, firstRow,
1205 velSystem, velUnit, doppler, refPos, refEpoch);
1206
1207// New output Table
1208
1209 SDMemTable* pTabOut = new SDMemTable(in,True);
1210
1211// Loop over rows in Table
1212
1213 const IPosition polChanAxes(2, asap::PolAxis, asap::ChanAxis);
1214 VelocityAligner<Float>::Method method = VelocityAligner<Float>::LINEAR;
1215 Bool extrapolate=False;
1216 Bool useCachedAbcissa = False;
1217 Bool first = True;
1218 Bool ok;
1219 Vector<Float> yOut;
1220 Vector<Bool> maskOut;
1221 uInt ifIdx, vaIdx;
1222//
1223 for (uInt iRow=0; iRow<nRows; ++iRow) {
1224 if (iRow%10==0) {
1225 cerr << "Processing row " << iRow << endl;
1226 }
1227
1228// Get EPoch
1229
1230 Quantum<Double> tQ2(times[iRow],DAY);
1231 MVEpoch mv2(tQ2);
1232 MEpoch epoch(mv2, epochRef);
1233
1234// Get FreqID vector. One freqID per IF
1235
1236 fqIDCol.get(iRow, freqID);
1237
1238// Get copy of data
1239
1240 const MaskedArray<Float>& mArrIn(in.rowAsMaskedArray(iRow));
1241 Array<Float> values = mArrIn.getArray();
1242 Array<Bool> mask = mArrIn.getMask();
1243
1244// cerr << "values in = " << values(IPosition(4,0,0,0,0),IPosition(4,0,0,0,9)) << endl;
1245
1246// For each row, the Velocity abcissa will be the same regardless
1247// of polarization. For all other axes (IF and BEAM) the abcissa
1248// will change. So we iterate through the data by pol-chan planes
1249// to mimimize the work. At this point, I think the Direction
1250// is stored as the same for each beam. DOn't know where the
1251// offsets are or what to do about them right now. For now
1252// all beams get same position and velocoity abcissa.
1253
1254 ArrayIterator<Float> itValuesPlane(values, polChanAxes);
1255 ArrayIterator<Bool> itMaskPlane(mask, polChanAxes);
1256 while (!itValuesPlane.pastEnd()) {
1257
1258// Find the IF index and then the VA PtrBlock index
1259
1260 const IPosition& pos = itValuesPlane.pos();
1261 ifIdx = pos(asap::IFAxis);
1262 vaIdx = (srcIdx[iRow]*nFreqIDs) + freqID[ifIdx];
1263//
1264 VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1);
1265 VectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
1266//
1267 first = True;
1268 useCachedAbcissa=False;
1269 while (!itValuesVec.pastEnd()) {
1270 ok = vA[vaIdx]->align (yOut, maskOut, itValuesVec.vector(),
1271 itMaskVec.vector(), epoch, useCachedAbcissa,
1272 method, extrapolate);
1273 itValuesVec.vector() = yOut;
1274 itMaskVec.vector() = maskOut;
1275//
1276 itValuesVec.next();
1277 itMaskVec.next();
1278//
1279 if (first) {
1280 useCachedAbcissa = True;
1281 first = False;
1282 }
1283 }
1284//
1285 itValuesPlane.next();
1286 itMaskPlane.next();
1287 }
1288
1289// cerr << "values out = " << values(IPosition(4,0,0,0,0),IPosition(4,0,0,0,9)) << endl;
1290
1291// Create and put back
1292
1293 SDContainer sc = in.getSDContainer(iRow);
1294 putDataInSDC(sc, values, mask);
1295//
1296 pTabOut->putSDContainer(sc);
1297 }
1298
1299// Clean up PointerBlock
1300
1301 for (uInt i=0; i<vA.nelements(); i++) delete vA[i];
1302//
1303 return pTabOut;
1304}
1305
1306
1307void SDMath::fillSDC(SDContainer& sc,
1308 const Array<Bool>& mask,
1309 const Array<Float>& data,
1310 const Array<Float>& tSys,
1311 Int scanID, Double timeStamp,
1312 Double interval, const String& sourceName,
1313 const Vector<uInt>& freqID) const
1314{
1315// Data and mask
1316
1317 putDataInSDC(sc, data, mask);
1318
1319// TSys
1320
1321 sc.putTsys(tSys);
1322
1323// Time things
1324
1325 sc.timestamp = timeStamp;
1326 sc.interval = interval;
1327 sc.scanid = scanID;
1328//
1329 sc.sourcename = sourceName;
1330 sc.putFreqMap(freqID);
1331}
1332
1333void SDMath::normalize(MaskedArray<Float>& sum,
1334 const Array<Float>& sumSq,
1335 const Array<Float>& nPts,
1336 WeightType wtType, Int axis,
1337 Int nAxesSub) const
1338{
1339 IPosition pos2(nAxesSub,0);
1340//
1341 if (wtType==NONE) {
1342
1343// We just average by the number of points accumulated.
1344// We need to make a MA out of nPts so that no divide by
1345// zeros occur
1346
1347 MaskedArray<Float> t(nPts, (nPts>Float(0.0)));
1348 sum /= t;
1349 } else if (wtType==VAR) {
1350
1351// Normalize each spectrum by sum(1/var) where the variance
1352// is worked out for each spectrum
1353
1354 Array<Float>& data = sum.getRWArray();
1355 VectorIterator<Float> itData(data, axis);
1356 while (!itData.pastEnd()) {
1357 pos2 = itData.pos().getFirst(nAxesSub);
1358 itData.vector() /= sumSq(pos2);
1359 itData.next();
1360 }
1361 } else if (wtType==TSYS) {
1362 }
1363}
1364
1365
1366void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum,
1367 MaskedArray<Float>& sum, Array<Float>& sumSq,
1368 Array<Float>& nPts, Array<Float>& tSysSum,
1369 const Array<Float>& tSys, const Array<Float>& nInc,
1370 const Vector<Bool>& mask, Double time, Double interval,
1371 const Block<CountedPtr<SDMemTable> >& in,
1372 uInt iTab, uInt iRow, uInt axis,
1373 uInt nAxesSub, Bool useMask,
1374 WeightType wtType) const
1375{
1376
1377// Get data
1378
1379 MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
1380 Array<Float>& valuesIn = dataIn.getRWArray(); // writable reference
1381 const Array<Bool>& maskIn = dataIn.getMask(); // RO reference
1382//
1383 if (wtType==NONE) {
1384 const MaskedArray<Float> n(nInc,dataIn.getMask());
1385 nPts += n; // Only accumulates where mask==T
1386 } else if (wtType==VAR) {
1387
1388// We are going to average the data, weighted by the noise for each pol, beam and IF.
1389// So therefore we need to iterate through by spectrum (axis 3)
1390
1391 VectorIterator<Float> itData(valuesIn, axis);
1392 ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
1393 Float fac = 1.0;
1394 IPosition pos(nAxesSub,0);
1395//
1396 while (!itData.pastEnd()) {
1397
1398// Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
1399
1400 if (useMask) {
1401 MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
1402 fac = 1.0/variance(tmp);
1403 } else {
1404 MaskedArray<Float> tmp(itData.vector(),itMask.vector());
1405 fac = 1.0/variance(tmp);
1406 }
1407
1408// Scale data
1409
1410 itData.vector() *= fac; // Writes back into 'dataIn'
1411//
1412// Accumulate variance per if/pol/beam averaged over spectrum
1413// This method to get pos2 from itData.pos() is only valid
1414// because the spectral axis is the last one (so we can just
1415// copy the first nAXesSub positions out)
1416
1417 pos = itData.pos().getFirst(nAxesSub);
1418 sumSq(pos) += fac;
1419//
1420 itData.next();
1421 itMask.next();
1422 }
1423 } else if (wtType==TSYS) {
1424 }
1425
1426// Accumulate sum of (possibly scaled) data
1427
1428 sum += dataIn;
1429
1430// Accumulate Tsys, time, and interval
1431
1432 tSysSum += tSys;
1433 timeSum += time;
1434 intSum += interval;
1435 nAccum += 1;
1436}
1437
1438
1439
1440
1441void SDMath::getCursorLocation(IPosition& start, IPosition& end,
1442 const SDMemTable& in) const
1443{
1444 const uInt nDim = 4;
1445 const uInt i = in.getBeam();
1446 const uInt j = in.getIF();
1447 const uInt k = in.getPol();
1448 const uInt n = in.nChan();
1449//
1450 start.resize(nDim);
1451 start(0) = i;
1452 start(1) = j;
1453 start(2) = k;
1454 start(3) = 0;
1455//
1456 end.resize(nDim);
1457 end(0) = i;
1458 end(1) = j;
1459 end(2) = k;
1460 end(3) = n-1;
1461}
1462
1463
1464void SDMath::convertWeightString(WeightType& wtType, const String& weightStr) const
1465{
1466 String tStr(weightStr);
1467 tStr.upcase();
1468 if (tStr.contains(String("NONE"))) {
1469 wtType = NONE;
1470 } else if (tStr.contains(String("VAR"))) {
1471 wtType = VAR;
1472 } else if (tStr.contains(String("TSYS"))) {
1473 wtType = TSYS;
1474 throw(AipsError("T_sys weighting not yet implemented"));
1475 } else {
1476 throw(AipsError("Unrecognized weighting type"));
1477 }
1478}
1479
1480void SDMath::convertInterpString(Int& type, const String& interp) const
1481{
1482 String tStr(interp);
1483 tStr.upcase();
1484 if (tStr.contains(String("NEAR"))) {
1485 type = InterpolateArray1D<Float,Float>::nearestNeighbour;
1486 } else if (tStr.contains(String("LIN"))) {
1487 type = InterpolateArray1D<Float,Float>::linear;
1488 } else if (tStr.contains(String("CUB"))) {
1489 type = InterpolateArray1D<Float,Float>::cubic;
1490 } else if (tStr.contains(String("SPL"))) {
1491 type = InterpolateArray1D<Float,Float>::spline;
1492 } else {
1493 throw(AipsError("Unrecognized interpolation type"));
1494 }
1495}
1496
1497void SDMath::putDataInSDC(SDContainer& sc, const Array<Float>& data,
1498 const Array<Bool>& mask) const
1499{
1500 sc.putSpectrum(data);
1501//
1502 Array<uChar> outflags(data.shape());
1503 convertArray(outflags,!mask);
1504 sc.putFlags(outflags);
1505}
1506
1507Table SDMath::readAsciiFile (const String& fileName) const
1508{
1509 String formatString;
1510 Table tbl = readAsciiTable (formatString, Table::Memory, fileName, "", "", False);
1511 return tbl;
1512}
1513
1514
1515
1516void SDMath::correctFromAsciiTable(SDMemTable* pTabOut,
1517 const SDMemTable& in, const String& fileName,
1518 const String& col0, const String& col1,
1519 const String& methodStr, Bool doAll,
1520 const Vector<Float>& xOut) const
1521{
1522
1523// Read gain-elevation ascii file data into a Table.
1524
1525 Table geTable = readAsciiFile (fileName);
1526//
1527 correctFromTable (pTabOut, in, geTable, col0, col1, methodStr, doAll, xOut);
1528}
1529
1530void SDMath::correctFromTable(SDMemTable* pTabOut, const SDMemTable& in,
1531 const Table& tTable, const String& col0,
1532 const String& col1,
1533 const String& methodStr, Bool doAll,
1534 const Vector<Float>& xOut) const
1535{
1536
1537// Get data from Table
1538
1539 ROScalarColumn<Float> geElCol(tTable, col0);
1540 ROScalarColumn<Float> geFacCol(tTable, col1);
1541 Vector<Float> xIn = geElCol.getColumn();
1542 Vector<Float> yIn = geFacCol.getColumn();
1543 Vector<Bool> maskIn(xIn.nelements(),True);
1544
1545// Interpolate (and extrapolate) with desired method
1546
1547 Int method = 0;
1548 convertInterpString(method, methodStr);
1549//
1550 Vector<Float> yOut;
1551 Vector<Bool> maskOut;
1552 InterpolateArray1D<Float,Float>::interpolate(yOut, maskOut, xOut,
1553 xIn, yIn, maskIn, method,
1554 True, True);
1555// Apply
1556
1557 correctFromVector (pTabOut, in, doAll, yOut);
1558}
1559
1560
1561void SDMath::correctFromVector (SDMemTable* pTabOut, const SDMemTable& in,
1562 Bool doAll, const Vector<Float>& factor) const
1563{
1564
1565// For operations only on specified cursor location
1566
1567 IPosition start, end;
1568 getCursorLocation(start, end, in);
1569
1570// Loop over rows and apply correction factor
1571
1572 const uInt axis = asap::ChanAxis;
1573 for (uInt i=0; i < in.nRow(); ++i) {
1574
1575// Get data
1576
1577 MaskedArray<Float> dataIn(in.rowAsMaskedArray(i));
1578
1579// Apply factor
1580
1581 if (doAll) {
1582 dataIn *= factor[i];
1583 } else {
1584 MaskedArray<Float> dataIn2 = dataIn(start,end); // reference
1585 dataIn2 *= factor[i];
1586 }
1587
1588// Write out
1589
1590 SDContainer sc = in.getSDContainer(i);
1591 putDataInSDC(sc, dataIn.getArray(), dataIn.getMask());
1592//
1593 pTabOut->putSDContainer(sc);
1594 }
1595}
1596
1597
1598void SDMath::generateSourceTable (Vector<String>& srcTab,
1599 Vector<uInt>& srcIdx,
1600 Vector<uInt>& firstRow,
1601 const Vector<String>& srcNames) const
1602//
1603// This algorithm assumes that if there are multiple beams
1604// that the source names are diffent. Oterwise we would need
1605// to look atthe direction for each beam...
1606//
1607{
1608 const uInt nRow = srcNames.nelements();
1609 srcTab.resize(0);
1610 srcIdx.resize(nRow);
1611 firstRow.resize(0);
1612//
1613 uInt nSrc = 0;
1614 for (uInt i=0; i<nRow; i++) {
1615 String srcName = srcNames[i];
1616
1617// Do we have this source already ?
1618
1619 Int idx = -1;
1620 if (nSrc>0) {
1621 for (uInt j=0; j<nSrc; j++) {
1622 if (srcName==srcTab[j]) {
1623 idx = j;
1624 break;
1625 }
1626 }
1627 }
1628
1629// Add new entry if not found
1630
1631 if (idx==-1) {
1632 nSrc++;
1633 srcTab.resize(nSrc,True);
1634 srcTab(nSrc-1) = srcName;
1635 idx = nSrc-1;
1636//
1637 firstRow.resize(nSrc,True);
1638 firstRow(nSrc-1) = i; // First row for which this source occurs
1639 }
1640
1641// Set index for this row
1642
1643 srcIdx[i] = idx;
1644 }
1645}
1646
1647MEpoch SDMath::epochFromString (const String& str, MEpoch::Types timeRef) const
1648{
1649 Quantum<Double> qt;
1650 if (MVTime::read(qt,str)) {
1651 MVEpoch mv(qt);
1652 MEpoch me(mv, timeRef);
1653 return me;
1654 } else {
1655 throw(AipsError("Invalid format for Epoch string"));
1656 }
1657}
1658
1659
1660String SDMath::formatEpoch(const MEpoch& epoch) const
1661{
1662 MVTime mvt(epoch.getValue());
1663 return mvt.string(MVTime::YMD) + String(" (") + epoch.getRefString() + String(")");
1664}
1665
1666
1667void SDMath::generateVelocityAligners (PtrBlock<VelocityAligner<Float>* >& vA,
1668 const SDMemTable& in, uInt nChan,
1669 uInt nFreqIDs, uInt nSrcTab,
1670 const Vector<uInt>& firstRow,
1671 MFrequency::Types velSystem,
1672 const String& velUnit,
1673 MDoppler::Types doppler,
1674 const MPosition& refPos,
1675 const MEpoch& refEpoch) const
1676{
1677 for (uInt fqID=0; fqID<nFreqIDs; fqID++) {
1678 SpectralCoordinate sC = in.getSpectralCoordinate(fqID);
1679 for (uInt iSrc=0; iSrc<nSrcTab; iSrc++) {
1680 MDirection refDir = in.getDirection(firstRow[iSrc]);
1681 uInt idx = (iSrc*nFreqIDs) + fqID;
1682 vA[idx] = new VelocityAligner<Float>(sC, nChan, refEpoch, refDir, refPos,
1683 velUnit, doppler, velSystem);
1684 }
1685 }
1686}
1687
Note: See TracBrowser for help on using the repository browser.