source: trunk/src/SDMath.cc@ 277

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

add reference time arg. to function 'VelocityAlignment'

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