source: trunk/src/SDMath.cc@ 478

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

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