source: trunk/src/SDMath.cc @ 502

Last change on this file since 502 was 502, checked in by kil064, 19 years ago

add function rotateLinPolPhase

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