source: trunk/src/SDMath.cc@ 268

Last change on this file since 268 was 267, checked in by kil064, 20 years ago

split velocitryAlign into towo functins

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