source: trunk/src/STSideBandSep.cpp @ 3106

Last change on this file since 3106 was 3106, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/No?

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...


Check-in asap modifications from Jim regarding casacore namespace conversion.
File size: 51.0 KB
Line 
1// C++ Interface: STSideBandSep
2//
3// Description:
4//    A class to invoke sideband separation of Scantable
5//
6// Author: Kana Sugimoto <kana.sugi@nao.ac.jp>, (C) 2012
7//
8// Copyright: See COPYING file that comes with this distribution
9//
10//
11
12// asap
13#include "STGrid.h"
14#include "STMath.h"
15#include "MathUtils.h"
16#include "STSideBandSep.h"
17
18// STL
19#include <ctype.h>
20
21// cascore
22#include <casa/OS/File.h>
23#include <casa/Logging/LogIO.h>
24#include <casa/Quanta/QuantumHolder.h>
25
26#include <measures/Measures/MFrequency.h>
27#include <measures/Measures/MCFrequency.h>
28
29#include <tables/Tables/TableRow.h>
30#include <tables/Tables/TableRecord.h>
31#include <tables/Tables/TableVector.h>
32
33using namespace std ;
34using namespace casacore ;
35using namespace asap ;
36
37// #ifndef KS_DEBUG
38// #define KS_DEBUG
39// #endif
40
41namespace asap {
42
43// constructors
44STSideBandSep::STSideBandSep(const vector<string> &names)
45{
46  LogIO os(LogOrigin("STSideBandSep","STSideBandSep()", WHERE));
47  os << "Setting scantable names to process." << LogIO::POST ;
48  // Set file names
49  ntable_ = names.size();
50  infileList_.resize(ntable_);
51  for (unsigned int i = 0; i < ntable_; i++){
52    if (!checkFile(names[i], "d"))
53      throw( AipsError("File does not exist") );
54    infileList_[i] = names[i];
55  }
56  intabList_.resize(0);
57
58  init();
59
60  {// Summary
61    os << ntable_ << " files are set: [";
62    for (unsigned int i = 0; i < ntable_; i++) {
63      os << " '" << infileList_[i] << "' ";
64      if (i != ntable_-1) os << ",";
65    }
66    os << "] " << LogIO::POST;
67  }
68};
69
70STSideBandSep::STSideBandSep(const vector<ScantableWrapper> &tables)
71{
72  LogIO os(LogOrigin("STSideBandSep","STSideBandSep()", WHERE));
73  os << "Setting list of scantables to process." << LogIO::POST ;
74  // Set file names
75  ntable_ = tables.size();
76  intabList_.resize(ntable_);
77  for (unsigned int i = 0; i < ntable_; i++){
78    intabList_[i] = tables[i].getCP();
79  }
80  infileList_.resize(0);
81
82  init();
83  tp_ = intabList_[0]->table().tableType();
84
85  os << ntable_ << " tables are set." << LogIO::POST;
86};
87
88STSideBandSep::~STSideBandSep()
89{
90};
91
92void STSideBandSep::init()
93{
94  // frequency setup
95  sigIfno_= -1;
96  ftol_ = -1;
97  solFrame_ = MFrequency::N_Types;
98  // shifts
99  initshift();
100  // direction tolerance
101  xtol_ = ytol_ = 9.69627e-6; // 2arcsec
102  // solution parameters
103  otherside_ = false;
104  doboth_ = false;
105  rejlimit_ = 0.2;
106  // LO1 values
107  lo1Freq_ = -1;
108  loTime_ = -1;
109  loDir_ = "";
110  // Default LO frame is TOPO
111  loFrame_ = MFrequency::TOPO;
112  // scantable storage
113  tp_ = Table::Memory;
114};
115
116void STSideBandSep::initshift()
117{
118  // shifts
119  nshift_ = 0;
120  nchan_ = 0;
121  sigShift_.resize(0);
122  imgShift_.resize(0);
123  tableList_.resize(0);
124};
125
126void STSideBandSep::setFrequency(const int ifno,
127                                 const string freqtol,
128                                 const string frame)
129{
130  LogIO os(LogOrigin("STSideBandSep","setFrequency()", WHERE));
131
132  initshift();
133
134  // IFNO
135  sigIfno_ = ifno;
136
137  // Frequency tolerance
138  Quantum<Double> qftol;
139  readQuantity(qftol, String(freqtol));
140  if (!qftol.getUnit().empty()){
141    // make sure the quantity is frequency
142    if (qftol.getFullUnit().getValue() != Unit("Hz").getValue())
143      throw( AipsError("Invalid quantity for frequency tolerance.") );
144    qftol.convert("Hz");
145  }
146  ftol_ = qftol;
147
148  // Frequency Frame
149  if (!frame.empty()){
150    MFrequency::Types mft;
151    if (!MFrequency::getType(mft, frame))
152      throw( AipsError("Invalid frame type.") );
153    solFrame_ = mft;
154  } else {
155    solFrame_ = MFrequency::N_Types;
156  }
157
158  {// Summary
159    const String sframe = ( (solFrame_ == MFrequency::N_Types) ?
160                            "table frame" :
161                            MFrequency::showType(solFrame_) );
162    os << "Frequency setup to search IF group: "
163       << "IFNO of table[0] = " << sigIfno_
164       << " , Freq tolerance = " << ftol_.getValue() << " [ "
165       << (ftol_.getUnit().empty() ? "channel" : ftol_.getUnit() )
166       << " ] (in " << sframe <<")" << LogIO::POST;
167  }
168};
169
170
171void STSideBandSep::setDirTolerance(const vector<string> dirtol)
172{
173  LogIO os(LogOrigin("STSideBandSep","setDirTolerance()", WHERE));
174  Quantum<Double> qcell;
175  if ( (dirtol.size() == 1) && !dirtol[0].empty() ) {
176    readQuantity(qcell, String(dirtol[0]));
177    if (qcell.getFullUnit().getValue() == Unit("rad").getValue())
178      xtol_ = ytol_ = qcell.getValue("rad");
179    else
180      throw( AipsError("Invalid unit for direction tolerance.") );
181  }
182  else if (dirtol.size() > 1) {
183    if ( dirtol[0].empty() && dirtol[1].empty() )
184      throw( AipsError("Direction tolerance is empty.") );
185    if ( !dirtol[0].empty() ) {
186      readQuantity(qcell, String(dirtol[0]));
187      if (qcell.getFullUnit().getValue() == Unit("rad").getValue())
188        xtol_ = qcell.getValue("rad");
189      else
190        throw( AipsError("Invalid unit for direction tolerance.") );
191    }
192    if ( !dirtol[1].empty() ) {
193      readQuantity(qcell, String(dirtol[1]));
194      if (qcell.getFullUnit().getValue() == Unit("rad").getValue())
195        ytol_ = qcell.getValue("rad");
196      else
197        throw( AipsError("Invalid unit for direction tolerance.") );
198    }
199    else {
200      ytol_ = xtol_;
201    }
202  }
203  else throw( AipsError("Invalid direction tolerance.") );
204
205  os << "Direction tolerance: ( "
206     << xtol_ << " , " << ytol_ << " ) [rad]" << LogIO::POST;
207};
208
209void STSideBandSep::setShift(const vector<double> &shift)
210{
211  LogIO os(LogOrigin("STSideBandSep","setShift()", WHERE));
212  imgShift_.resize(shift.size());
213  for (unsigned int i = 0; i < shift.size(); i++)
214    imgShift_[i] = shift[i];
215
216  if (imgShift_.size() == 0) {
217    os << "Channel shifts are cleared." << LogIO::POST;
218  } else {
219    os << "Channel shifts of image sideband are set: ( ";
220    for (unsigned int i = 0; i < imgShift_.size(); i++) {
221      os << imgShift_[i];
222      if (i != imgShift_.size()-1) os << " , ";
223    }
224    os << " ) [channels]" << LogIO::POST;
225  }
226};
227
228void STSideBandSep::setThreshold(const double limit)
229{
230  LogIO os(LogOrigin("STSideBandSep","setThreshold()", WHERE));
231  if (limit < 0)
232    throw( AipsError("Rejection limit should be a positive number.") );
233
234  rejlimit_ = limit;
235  os << "Rejection limit is set to " << rejlimit_ << LogIO::POST;
236};
237
238void STSideBandSep::separate(string outname)
239{
240  LogIO os(LogOrigin("STSideBandSep","separate()", WHERE));
241  if (outname.empty())
242    outname = "sbseparated.asap";
243
244  // Set up a goup of IFNOs in the list of scantables within
245  // the frequency tolerance and make them a list.
246  nshift_ = setupShift();
247  if (nshift_ < 2)
248    throw( AipsError("At least 2 IFs are necessary for convolution.") );
249  // Grid scantable and generate output tables
250  ScantableWrapper gridst = gridTable();
251  sigTab_p = gridst.getCP();
252  if (doboth_)
253    imgTab_p = gridst.copy().getCP();
254  vector<unsigned int> remRowIds;
255  remRowIds.resize(0);
256  Matrix<float> specMat(nchan_, nshift_);
257  Matrix<bool> flagMat(nchan_, nshift_);
258  vector<float> sigSpec(nchan_), imgSpec(nchan_);
259  Vector<bool> flagVec(nchan_);
260  vector<uInt> tabIdvec;
261
262  //Generate FFTServer
263  fftsf.resize(IPosition(1, nchan_), FFTEnums::REALTOCOMPLEX);
264  fftsi.resize(IPosition(1, nchan_), FFTEnums::COMPLEXTOREAL);
265
266  /// Loop over sigTab_p and separate sideband
267  for (int irow = 0; irow < sigTab_p->nrow(); irow++){
268    tabIdvec.resize(0);
269    const int polId = sigTab_p->getPol(irow);
270    const int beamId = sigTab_p->getBeam(irow);
271    const vector<double> dir = sigTab_p->getDirectionVector(irow);
272    // Get a set of spectra to solve
273    if (!getSpectraToSolve(polId, beamId, dir[0], dir[1],
274                           specMat, flagMat, tabIdvec)){
275      remRowIds.push_back(irow);
276#ifdef KS_DEBUG
277      cout << "no matching row found. skipping row = " << irow << endl;
278#endif
279      continue;
280    }
281    // Solve signal sideband
282    sigSpec = solve(specMat, tabIdvec, true);
283    sigTab_p->setSpectrum(sigSpec, irow);
284    if (sigTab_p->isAllChannelsFlagged(irow)){
285      // unflag the spectrum since there should be some valid data
286      sigTab_p->flagRow(vector<uInt>(irow), true);
287      // need to unflag whole channels anyway
288      sigTab_p->flag(irow, vector<bool>(), true);
289    }
290    // apply channel flag
291    flagVec = collapseFlag(flagMat, tabIdvec, true);
292    //boolVec = !boolVec; // flag
293    vector<bool> tmpflag;
294    flagVec.tovector(tmpflag);
295    sigTab_p->flag(irow, tmpflag, false);
296
297    // Solve image sideband
298    if (doboth_) {
299      imgSpec = solve(specMat, tabIdvec, false);
300      imgTab_p->setSpectrum(imgSpec, irow);
301      if (imgTab_p->isAllChannelsFlagged(irow)){
302        // unflag the spectrum since there should be some valid data
303        imgTab_p->flagRow(vector<uInt>(irow), true);
304        // need to unflag whole channels anyway
305        imgTab_p->flag(irow, vector<bool>(), true);
306      }
307      // apply channel flag
308      flagVec = collapseFlag(flagMat, tabIdvec, false);
309      //boolVec = !boolVec; // flag
310      flagVec.tovector(tmpflag);
311      imgTab_p->flag(irow, tmpflag, false);
312    }
313  } // end of row loop
314
315  // Remove or flag rows without relevant data from gridded tables
316  if (remRowIds.size() > 0) {
317    const size_t nrem = remRowIds.size();
318    if ( sigTab_p->table().canRemoveRow() ) {
319      sigTab_p->table().removeRow(remRowIds);
320      os << "Removing " << nrem << " rows from the signal band table"
321         << LogIO::POST;
322    } else {
323      sigTab_p->flagRow(remRowIds, false);
324      os << "Cannot remove rows from the signal band table. Flagging "
325         << nrem << " rows" << LogIO::POST;
326    }
327
328    if (doboth_) {
329      if ( imgTab_p->table().canRemoveRow() ) {
330        imgTab_p->table().removeRow(remRowIds);
331        os << "Removing " << nrem << " rows from the image band table"
332           << LogIO::POST;
333      } else {
334        imgTab_p->flagRow(remRowIds, false);
335        os << "Cannot remove rows from the image band table. Flagging "
336           << nrem << " rows" << LogIO::POST;
337      }
338    }
339  }
340
341  // Finally, save tables on disk
342  if (outname.size() ==0)
343    outname = "sbseparated.asap";
344  const string sigName = outname + ".signalband";
345  os << "Saving SIGNAL sideband table: " << sigName << LogIO::POST;
346  sigTab_p->makePersistent(sigName);
347  if (doboth_) {
348    solveImageFrequency();
349    const string imgName = outname + ".imageband";
350    os << "Saving IMAGE sideband table: " << sigName << LogIO::POST;
351    imgTab_p->makePersistent(imgName);
352  }
353
354};
355
356unsigned int STSideBandSep::setupShift()
357{
358  LogIO os(LogOrigin("STSideBandSep","setupShift()", WHERE));
359  if (infileList_.size() == 0 && intabList_.size() == 0)
360    throw( AipsError("No scantable has been set. Set a list of scantables first.") );
361
362  const bool byname = (intabList_.size() == 0);
363  // Make sure sigIfno_ exists in the first table.
364  CountedPtr<Scantable> stab;
365  vector<string> coordsav;
366  vector<string> coordtmp(3);
367  os << "Checking IFNO in the first table." << LogIO::POST;
368  if (byname) {
369    if (!checkFile(infileList_[0], "d"))
370      os << LogIO::SEVERE << "Could not find scantable '" << infileList_[0]
371         << "'" << LogIO::EXCEPTION;
372    stab = CountedPtr<Scantable>(new Scantable(infileList_[0]));
373  } else {
374    stab = intabList_[0];
375  }
376  if (sigIfno_ < 0) {
377    sigIfno_ = (int) stab->getIF(0);
378    os << "IFNO to process has not been set. Using the first IF = "
379       << sigIfno_ << LogIO::POST;
380  }
381
382  unsigned int basench;
383  double basech0, baseinc, ftolval, inctolval;
384  coordsav = stab->getCoordInfo();
385  const string stfframe = coordsav[1];
386  coordtmp[0] = "Hz";
387  coordtmp[1] = ( (solFrame_ == MFrequency::N_Types) ?
388                  stfframe :
389                  MFrequency::showType(solFrame_) );
390  coordtmp[2] = coordsav[2];
391  stab->setCoordInfo(coordtmp);
392  if (!getFreqInfo(stab, (unsigned int) sigIfno_, basech0, baseinc, basench)) {
393    os << LogIO::SEVERE << "No data with IFNO=" << sigIfno_
394       << " found in the first table." << LogIO::EXCEPTION;
395  }
396  else {
397    os << "Found IFNO = " << sigIfno_
398       << " in the first table." << LogIO::POST;
399#ifdef KS_DEBUG
400    cout << "Frequencies in the first table: freq0 = " << basech0 << ", incr = " << baseinc << ", nch = " << basench << endl;
401#endif
402  }
403  if (ftol_.getUnit().empty()) {
404    // tolerance in unit of channels
405    ftolval = abs(ftol_.getValue() * baseinc);
406  }
407  else {
408    ftolval = abs(ftol_.getValue("Hz"));
409  }
410  inctolval = abs(baseinc/(double) basench);
411  const string poltype0 = stab->getPolType();
412#ifdef KS_DEBUG
413  cout << "Looking for pairs that satisfies freqtol = " << ftolval << ", inctol = " << inctolval << endl;
414#endif
415
416  // Initialize shift values
417  initshift();
418
419//   const bool setImg = ( doboth_ && (imgShift_.size() == 0) );
420  const bool setImg = (imgShift_.size() == 0);
421  // Select IFs
422  for (unsigned int itab = 0; itab < ntable_; itab++ ){
423    os << "Table " << itab << LogIO::POST;
424    if (itab > 0) {
425      if (byname) {
426        if (!checkFile(infileList_[itab], "d"))
427          os << LogIO::SEVERE << "Could not find scantable '"
428             << infileList_[itab] << "'" << LogIO::EXCEPTION;
429        stab = CountedPtr<Scantable>(new Scantable(infileList_[itab]));
430      } else {
431        stab = intabList_[itab];
432      }
433      //POLTYPE should be the same.
434      if (stab->getPolType() != poltype0 ) {
435        os << LogIO::WARN << "POLTYPE differs from the first table."
436           << " Skipping the table" << LogIO::POST;
437        continue;
438      }
439      // Multi beam data may not handled properly
440      if (stab->nbeam() > 1)
441        os <<  LogIO::WARN << "Table contains multiple beams. "
442           << "It may not be handled properly."  << LogIO::POST;
443
444      coordsav = stab->getCoordInfo();
445      coordtmp[2] = coordsav[2];
446      stab->setCoordInfo(coordtmp);
447    }
448    bool selected = false;
449    vector<uint> ifnos = stab->getIFNos();
450    vector<uint>::iterator iter;
451    const STSelector& basesel = stab->getSelection();
452    for (iter = ifnos.begin(); iter != ifnos.end(); iter++){
453      unsigned int nch;
454      double freq0, incr;
455      if ( getFreqInfo(stab, *iter, freq0, incr, nch) ){
456#ifdef KS_DEBUG
457        cout << "Checking IF=" << *iter << " int Table " << itab <<": freq0 = " << freq0 << ", incr = " << incr << ", nch = " << nch << endl;
458        cout << "- number of channels: " << nch << endl;
459        cout << "- chan0 difference: " << abs(freq0-basech0) << " (threshold: " << ftolval << ")" << endl;
460        cout << "- inc difference: " << abs(incr-baseinc) << " (threshold: " << inctolval << ")" << endl;
461#endif
462        if ( (nch == basench) && (abs(freq0-basech0) < ftolval)
463             && (abs(incr-baseinc) < inctolval) ){
464          //Found
465          STSelector sel(basesel);
466          sel.setIFs(vector<int>(1,(int) *iter));
467          stab->setSelection(sel);
468          CountedPtr<Scantable> seltab = ( new Scantable(*stab, false) );
469          stab->setSelection(basesel);
470          seltab->setCoordInfo(coordsav);
471          const double chShift = (freq0 - basech0) / baseinc;
472          tableList_.push_back(seltab);
473          sigShift_.push_back(-chShift);
474          if (setImg)
475            imgShift_.push_back(chShift);
476
477          selected = true;
478          os << "- IF" << *iter << " selected: sideband shift = "
479             << chShift << " channels" << LogIO::POST;
480        }
481      }
482    } // ifno loop
483    stab->setCoordInfo(coordsav);
484    if (!selected)
485      os << LogIO::WARN << "No data selected in Table "
486         << itab << LogIO::POST;
487  } // table loop
488  nchan_ = basench;
489
490  os << "Total number of IFs selected = " << tableList_.size()
491     << LogIO::POST;
492  if ( setImg && (sigShift_.size() != imgShift_.size()) ){
493      os << LogIO::SEVERE
494         << "User defined channel shift of image sideband has "
495         << imgShift_.size()
496         << " elements, while selected IFNOs are " << sigShift_.size()
497         << "\nThe frequency tolerance (freqtol) may be too small."
498         << LogIO::EXCEPTION;
499  }
500
501  return tableList_.size();
502};
503
504bool STSideBandSep::getFreqInfo(const CountedPtr<Scantable> &stab,
505                                const unsigned int &ifno,
506                                double &freq0, double &incr,
507                                unsigned int &nchan)
508{
509    vector<uint> ifnos = stab->getIFNos();
510    bool found = false;
511    vector<uint>::iterator iter;
512    for (iter = ifnos.begin(); iter != ifnos.end(); iter++){
513      if (*iter == ifno) {
514        found = true;
515        break;
516      }
517    }
518    if (!found)
519      return false;
520
521    const STSelector& basesel = stab->getSelection();
522    STSelector sel(basesel);
523    sel.setIFs(vector<int>(1,(int) ifno));
524    stab->setSelection(sel);
525    vector<double> freqs;
526    freqs = stab->getAbcissa(0);
527    freq0 = freqs[0];
528    incr = freqs[1] - freqs[0];
529    nchan = freqs.size();
530    stab->setSelection(basesel);
531    return true;
532};
533
534ScantableWrapper STSideBandSep::gridTable()
535{
536  LogIO os(LogOrigin("STSideBandSep","gridTable()", WHERE));
537  if (tableList_.size() == 0)
538    throw( AipsError("Internal error. No scantable has been set to grid.") );
539  Double xmin, xmax, ymin, ymax;
540  mapExtent(tableList_, xmin, xmax, ymin, ymax);
541  const Double centx = 0.5 * (xmin + xmax);
542  const Double centy = 0.5 * (ymin + ymax);
543  const int nx = max(1, (int) ceil( (xmax - xmin) * cos(centy) /xtol_ ) );
544  const int ny = max(1, (int) ceil( (ymax - ymin) / ytol_ ) );
545
546  string scellx, scelly;
547  {
548    ostringstream oss;
549    oss << xtol_ << "rad" ;
550    scellx = oss.str();
551  }
552  {
553    ostringstream oss;
554    oss << ytol_ << "rad" ;
555    scelly = oss.str();
556  }
557
558  ScantableWrapper stab0;
559  if (intabList_.size() > 0)
560    stab0 = ScantableWrapper(intabList_[0]);
561  else
562    stab0 = ScantableWrapper(infileList_[0]);
563
564  string scenter;
565  {
566    ostringstream oss;
567    oss << stab0.getCP()->getDirectionRefString() << " "
568        << centx << "rad" << " " << centy << "rad";
569    scenter = oss.str();
570  }
571 
572  STGrid2 gridder = STGrid2(stab0);
573  gridder.setIF(sigIfno_);
574  gridder.defineImage(nx, ny, scellx, scelly, scenter);
575  //  gridder.setFunc("box", 1); // convsupport=1 fails
576  gridder.setFunc("box");
577  gridder.setWeight("uniform");
578#ifdef KS_DEBUG
579  cout << "Grid parameter summary: " << endl;
580  cout << "- IF = " << sigIfno_ << endl;
581  cout << "- center = " << scenter << "\n"
582       << "- npix = (" << nx << ", " << ny << ")\n"
583       << "- cell = (" << scellx << ", " << scelly << ")" << endl;
584#endif
585  gridder.grid();
586  const int itp = (tp_ == Table::Memory ? 0 : 1);
587  ScantableWrapper gtab = gridder.getResultAsScantable(itp);
588  // WORKAROUND : Shift TIME for proper pointing resolution in future imaging.
589  shiftTimeInGriddedST(gtab.getCP());
590  return gtab;
591};
592
593void STSideBandSep::mapExtent(vector< CountedPtr<Scantable> > &tablist,
594                              Double &xmin, Double &xmax,
595                              Double &ymin, Double &ymax)
596{
597  ROArrayColumn<Double> dirCol_;
598  dirCol_.attach( tablist[0]->table(), "DIRECTION" );
599  Matrix<Double> direction = dirCol_.getColumn();
600  Vector<Double> ra( direction.row(0) );
601  mathutil::rotateRA(ra);
602  minMax( xmin, xmax, ra );
603  minMax( ymin, ymax, direction.row(1) );
604  Double amin, amax, bmin, bmax;
605  const uInt ntab = tablist.size();
606  for ( uInt i = 1 ; i < ntab ; i++ ) {
607    dirCol_.attach( tablist[i]->table(), "DIRECTION" );
608    direction.assign( dirCol_.getColumn() );
609    ra.assign( direction.row(0) );
610    mathutil::rotateRA(ra);
611    minMax( amin, amax, ra );
612    minMax( bmin, bmax, direction.row(1) );
613    xmin = min(xmin, amin);
614    xmax = max(xmax, amax);
615    ymin = min(ymin, bmin);
616    ymax = max(ymax, bmax);
617  }
618};
619
620// STGrid sets the identical time for all rows in scantable
621// which is reasonable thing to do in position based averaging.
622// However, this prevents CASA from finding proper pointing
623// per spectra once the gridded scantable is converted to
624// measurement set (MS). It is because MS does not
625// have ability to store per spectra pointing information.
626// MS stores pointing information in a subtable, POINTING,
627// with corresponding TIME when an antenna pointed the direction.
628// The pointing direction corresponding to a spectra is resolved
629// in MS by interpolating DIRECTION in POINTING subtable in TIME
630// the spectra is observed. If there are multiple match,
631// the first match is adopted. Therefore, gridded table (whose TIME
632// is set to a single value) is misunderstood in MS that all data
633// come from a single pointing.
634// The function workarounds this defect by artificially shifting
635// TIME by INTERVAL in each row.
636void STSideBandSep::shiftTimeInGriddedST(const CountedPtr<Scantable> &stab)
637{
638  LogIO os(LogOrigin("STSideBandSep", "shiftTimeInGriddedST()", WHERE));
639  // Gridded table usually has an IF and a BEAM.
640  {
641    std::vector<uint> bmnos = stab->getBeamNos();
642    if (bmnos.size() > 1)
643      throw( AipsError("Multiple BEAMNOs found in the scantable. This may not a gridded table") );
644    std::vector<uint> ifnos = stab->getIFNos();
645    if (ifnos.size() > 1)
646      throw( AipsError("Multiple IFNOs found in the scantable. This may not a gridded table") );
647  }
648  // Rows in gridded table usually sorted by DIRECTION
649  const Table& tab = stab->table();
650  ScalarColumn<Double> mjdCol( tab, "TIME");
651  ROScalarColumn<Double> intCol( tab, "INTERVAL");
652  ROArrayColumn<Double> dirCol( tab, "DIRECTION");
653  Matrix<Double> direction = dirCol.getColumn();
654  Vector<Double> ra( direction.row(0) );
655  Vector<Double> dec( direction.row(1) );
656  Double prevTime, prevInt, prevRA(ra[0]), prevDec(dec[0]);
657  mjdCol.get(0, prevTime);
658  intCol.get(0, prevInt);
659  Double currInt, currRA, currDec;
660  Double dx(xtol_*0.95), dy(ytol_*0.95);
661  Double secToDay(1./24./3600.);
662  for (int irow = 0; irow < stab->nrow(); ++irow){
663    currRA = ra[irow];
664    currDec = dec[irow];
665    if ((prevRA+dx-currRA)*(currRA-prevRA+dx)>=0 &&
666        (prevDec+dy-currDec)*(currDec-prevDec+dy)>=0) {
667      // the same time stamp as the previous row
668      mjdCol.put(irow, prevTime);
669      // remember the longest interval
670      intCol.get(irow, currInt);
671      if (currInt > prevInt) prevInt = currInt;
672    } else {
673      // a new direction. need to set new time stamp.
674      prevTime += prevInt*secToDay;
675      mjdCol.put(irow, prevTime);
676      // new interval and direction
677      intCol.get(irow, prevInt);
678      prevRA = currRA;
679      prevDec = currDec;
680    }
681  }
682};
683
684// bool STSideBandSep::getSpectraToSolve(const int polId, const int beamId,
685//                                    const double dirX, const double dirY,
686//                                    Matrix<float> &specMat, vector<uInt> &tabIdvec)
687bool STSideBandSep::getSpectraToSolve(const int polId, const int beamId,
688                                      const double dirX, const double dirY,
689                                      Matrix<float> &specMat,
690                                      Matrix<bool> &flagMat,
691                                      vector<uInt> &tabIdvec)
692{
693  LogIO os(LogOrigin("STSideBandSep","getSpectraToSolve()", WHERE));
694
695  tabIdvec.resize(0);
696  specMat.resize(nchan_, nshift_);
697  Vector<float> spec;
698  Vector<bool> boolVec;
699  uInt nspec = 0;
700  STMath stm(false); // insitu has no effect for average.
701  for (uInt itab = 0 ; itab < nshift_ ; itab++) {
702    CountedPtr<Scantable> currtab_p = tableList_[itab];
703    // Selection by POLNO and BEAMNO
704    const STSelector& basesel = currtab_p->getSelection();
705    STSelector sel(basesel);
706    sel.setPolarizations(vector<int>(1, polId));
707    sel.setBeams(vector<int>(1, beamId));
708    try {
709      currtab_p->setSelection(sel);
710    } catch (...) {
711#ifdef KS_DEBUG
712      cout << "Table " << itab << " - No spectrum found. skipping the table."
713           << endl;
714#endif
715      continue;
716    }
717    // Selection by direction;
718    vector<int> selrow(0);
719    vector<double> currDir(2, 0.);
720    const int nselrow = currtab_p->nrow();
721    for (int irow = 0 ; irow < nselrow ; irow++) {
722      currDir = currtab_p->getDirectionVector(irow);
723      if ( (abs(currDir[0]-dirX) > xtol_) ||
724           (abs(currDir[1]-dirY) > ytol_) )
725        continue;
726      // within direction tolerance
727      selrow.push_back(irow);
728    } // end of row loop
729
730    if (selrow.size() < 1) {
731      currtab_p->setSelection(basesel);
732
733#ifdef KS_DEBUG
734      cout << "Table " << itab << " - No spectrum found. skipping the table."
735           << endl;
736#endif
737
738      continue;
739    }
740
741    // At least a spectrum is selected in this table
742    CountedPtr<Scantable> seltab_p = ( new Scantable(*currtab_p, false) );
743    currtab_p->setSelection(basesel);
744    STSelector sel2(seltab_p->getSelection());
745    sel2.setRows(selrow);
746    seltab_p->setSelection(sel2);
747    CountedPtr<Scantable> avetab_p;
748    if (seltab_p->nrow() > 1) {
749      // STMath::average also merges FLAGTRA and FLAGROW
750      avetab_p = stm.average(vector< CountedPtr<Scantable> >(1, seltab_p),
751                             vector<bool>(), "TINTSYS", "NONE");
752#ifdef KS_DEBUG
753      cout << "Table " << itab << " - more than a spectrum is selected. averaging rows..."
754           << endl;
755#endif
756      if (avetab_p->nrow() > 1)
757        throw( AipsError("Averaged table has more than a row. Somethigs went wrong.") );
758    } else {
759      avetab_p = seltab_p;
760    }
761    // Check FLAGTRA and FLAGROW if there's any valid channel in the spectrum
762    if (avetab_p->getFlagRow(0) || avetab_p->isAllChannelsFlagged(0)) {
763#ifdef KS_DEBUG
764      cout << "Table " << itab << " - All data are flagged. skipping the table."
765           << endl;
766#endif
767      continue;
768    }
769    // Interpolate flagged channels of the spectrum.
770    Vector<Float> tmpSpec = avetab_p->getSpectrum(0);
771    // Mask is true if the data is valid (!flag)
772    vector<bool> mask = avetab_p->getMask(0);
773    mathutil::doZeroOrderInterpolation(tmpSpec, mask);
774    spec.reference(specMat.column(nspec));
775    spec = tmpSpec;
776    boolVec.reference(flagMat.column(nspec));
777    boolVec = mask; // cast std::vector to casa::Vector
778    boolVec = !boolVec;
779    tabIdvec.push_back((uInt) itab);
780    nspec++;
781    //Liberate from reference
782    spec.unique();
783    boolVec.unique();
784  } // end of table loop
785  // Check the number of selected spectra and resize matrix.
786  if (nspec != nshift_){
787    //shiftSpecmat.resize(nchan_, nspec, true);
788    specMat.resize(nchan_, nspec, true);
789    flagMat.resize(nchan_, nspec, true);
790#ifdef KS_DEBUG
791      cout << "Could not find corresponding rows in some tables."
792           << endl;
793      cout << "Number of spectra selected = " << nspec << endl;
794#endif
795  }
796  if (nspec < 2) {
797#ifdef KS_DEBUG
798      cout << "At least 2 spectra are necessary for convolution"
799           << endl;
800#endif
801      return false;
802  }
803  return true;
804};
805
806
807Vector<bool> STSideBandSep::collapseFlag(const Matrix<bool> &flagMat,
808                                         const vector<uInt> &tabIdvec,
809                                         const bool signal)
810{
811  LogIO os(LogOrigin("STSideBandSep","collapseFlag()", WHERE));
812  if (tabIdvec.size() == 0)
813    throw(AipsError("Internal error. Table index is not defined."));
814  if (flagMat.ncolumn() != tabIdvec.size())
815    throw(AipsError("Internal error. The row number of input matrix is not conformant."));
816  if (flagMat.nrow() != nchan_)
817    throw(AipsError("Internal error. The channel size of input matrix is not conformant."));
818 
819  const size_t nspec = tabIdvec.size();
820  vector<double> *thisShift;
821  if (signal == otherside_) {
822    // (solve signal && solveother = T) OR (solve image && solveother = F)
823    thisShift = &imgShift_;
824  } else {
825    // (solve signal && solveother = F) OR (solve image && solveother = T)
826    thisShift =  &sigShift_;
827 }
828
829  Vector<bool> outflag(nchan_, false);
830  double tempshift;
831  Vector<bool> shiftvec(nchan_, false);
832  Vector<bool> accflag(nchan_, false);
833  uInt shiftId;
834  for (uInt i = 0 ; i < nspec; ++i) {
835    shiftId = tabIdvec[i];
836    tempshift = - thisShift->at(shiftId);
837    shiftFlag(flagMat.column(i), tempshift, shiftvec);
838    // Now accumulate Flag
839    for (uInt j = 0 ; j < nchan_ ; ++j)
840      accflag[j] |= shiftvec[j];
841  }
842  outflag = accflag;
843  // Shift back Flag
844  //cout << "Shifting FLAG back to " << thisShift->at(0) << " channels" << endl;
845  //shiftFlag(accflag, thisShift->at(0), outflag);
846
847  return outflag;
848}
849
850
851vector<float> STSideBandSep::solve(const Matrix<float> &specmat,
852                                   const vector<uInt> &tabIdvec,
853                                   const bool signal)
854{
855  LogIO os(LogOrigin("STSideBandSep","solve()", WHERE));
856  if (tabIdvec.size() == 0)
857    throw(AipsError("Internal error. Table index is not defined."));
858  if (specmat.ncolumn() != tabIdvec.size())
859    throw(AipsError("Internal error. The row number of input matrix is not conformant."));
860  if (specmat.nrow() != nchan_)
861    throw(AipsError("Internal error. The channel size of input matrix is not conformant."));
862 
863
864#ifdef KS_DEBUG
865  cout << "Solving " << (signal ? "SIGNAL" : "IMAGE") << " sideband."
866     << endl;
867#endif
868
869  const size_t nspec = tabIdvec.size();
870  vector<double> *thisShift, *otherShift;
871  if (signal == otherside_) {
872    // (solve signal && solveother = T) OR (solve image && solveother = F)
873    thisShift = &imgShift_;
874    otherShift = &sigShift_;
875#ifdef KS_DEBUG
876    cout << "Image sideband will be deconvolved." << endl;
877#endif
878  } else {
879    // (solve signal && solveother = F) OR (solve image && solveother = T)
880    thisShift =  &sigShift_;
881    otherShift = &imgShift_;
882#ifdef KS_DEBUG
883    cout << "Signal sideband will be deconvolved." << endl;
884#endif
885 }
886
887  vector<double> spshift(nspec);
888  Matrix<float> shiftSpecmat(nchan_, nspec, 0.);
889  double tempshift;
890  Vector<float> shiftspvec;
891  uInt shiftId;
892  for (uInt i = 0 ; i < nspec; i++) {
893    shiftId = tabIdvec[i];
894    spshift[i] = otherShift->at(shiftId) - thisShift->at(shiftId);
895    tempshift = - thisShift->at(shiftId);
896    shiftspvec.reference(shiftSpecmat.column(i));
897    shiftSpectrum(specmat.column(i), tempshift, shiftspvec);
898  }
899
900  Matrix<float> convmat(nchan_, nspec*(nspec-1)/2, 0.);
901  vector<float> thisvec(nchan_, 0.);
902
903  float minval, maxval;
904  minMax(minval, maxval, shiftSpecmat);
905#ifdef KS_DEBUG
906  cout << "Max/Min of input Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
907#endif
908
909#ifdef KS_DEBUG
910  cout << "starting deconvolution" << endl;
911#endif
912  deconvolve(shiftSpecmat, spshift, rejlimit_, convmat);
913#ifdef KS_DEBUG
914  cout << "finished deconvolution" << endl;
915#endif
916
917  minMax(minval, maxval, convmat);
918#ifdef KS_DEBUG
919  cout << "Max/Min of output Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
920#endif
921
922  aggregateMat(convmat, thisvec);
923
924  if (!otherside_) return thisvec;
925
926  // subtract from the other side band.
927  vector<float> othervec(nchan_, 0.);
928  subtractFromOther(shiftSpecmat, thisvec, spshift, othervec);
929  return othervec;
930};
931
932
933void STSideBandSep::shiftSpectrum(const Vector<float> &invec,
934                                  double shift,
935                                  Vector<float> &outvec)
936{
937  LogIO os(LogOrigin("STSideBandSep","shiftSpectrum()", WHERE));
938  if (invec.size() != nchan_)
939    throw(AipsError("Internal error. The length of input vector differs from nchan_"));
940  if (outvec.size() != nchan_)
941    throw(AipsError("Internal error. The length of output vector differs from nchan_"));
942
943#ifdef KS_DEBUG
944  cout << "Start shifting spectrum for " << shift << " channels" << endl;
945#endif
946
947  // tweak shift to be in 0 ~ nchan_-1
948  if ( fabs(shift) > nchan_ ) shift = fmod(shift, nchan_);
949  if (shift < 0.) shift += nchan_;
950  double rweight = fmod(shift, 1.);
951  if (rweight < 0.) rweight += 1.;
952  double lweight = 1. - rweight;
953  uInt lchan, rchan;
954
955  outvec = 0;
956  for (uInt i = 0 ; i < nchan_ ; i++) {
957    lchan = uInt( floor( fmod( (i + shift), nchan_ ) ) );
958    if (lchan < 0.) lchan += nchan_;
959    rchan = ( (lchan + 1) % nchan_ );
960    outvec(lchan) += invec(i) * lweight;
961    outvec(rchan) += invec(i) * rweight;
962#ifdef KS_DEBUG
963    if (i == 2350 || i== 2930) {
964      cout << "Channel= " << i << " of input vector: " << endl;
965      cout << "L channel = " << lchan << endl;
966      cout << "R channel = " << rchan << endl;
967      cout << "L weight = " << lweight << endl;
968      cout << "R weight = " << rweight << endl;
969    }
970#endif
971  }
972};
973
974
975void STSideBandSep::shiftFlag(const Vector<bool> &invec,
976                                  double shift,
977                                  Vector<bool> &outvec)
978{
979  LogIO os(LogOrigin("STSideBandSep","shiftFlag()", WHERE));
980  if (invec.size() != nchan_)
981    throw(AipsError("Internal error. The length of input vector differs from nchan_"));
982  if (outvec.size() != nchan_)
983    throw(AipsError("Internal error. The length of output vector differs from nchan_"));
984
985#ifdef KS_DEBUG
986  cout << "Start shifting flag for " << shift << "channels" << endl;
987#endif
988
989  // shift is almost integer think it as int.
990  // tolerance should be in 0 - 1
991  double tolerance = 0.01;
992  // tweak shift to be in 0 ~ nchan_-1
993  if ( fabs(shift) > nchan_ ) shift = fmod(shift, nchan_);
994  if (shift < 0.) shift += nchan_;
995  double rweight = fmod(shift, 1.);
996  bool ruse(true), luse(true);
997  if (rweight < 0.) rweight += 1.;
998  if (rweight < tolerance){
999    // the shift is almost lchan
1000    ruse = false;
1001    luse = true;
1002  }
1003  if (rweight > 1-tolerance){
1004    // the shift is almost rchan
1005    ruse = true;
1006    luse = false;
1007  }
1008  uInt lchan, rchan;
1009
1010  outvec = false;
1011  for (uInt i = 0 ; i < nchan_ ; i++) {
1012    lchan = uInt( floor( fmod( (i + shift), nchan_ ) ) );
1013    if (lchan < 0.) lchan += nchan_;
1014    rchan = ( (lchan + 1) % nchan_ );
1015    outvec(lchan) |= (invec(i) && luse);
1016    outvec(rchan) |= (invec(i) && ruse);
1017#ifdef KS_DEBUG
1018    if (i == 2350 || i == 2930) {
1019      cout << "Channel= " << i << " of input vector: " << endl;
1020      cout << "L channel = " << lchan << endl;
1021      cout << "R channel = " << rchan << endl;
1022      cout << "L channel will be " << (luse ? "used" : "ignored") << endl;
1023      cout << "R channel will be " << (ruse ? "used" : "ignored") << endl;
1024    }
1025#endif
1026  }
1027};
1028
1029
1030void STSideBandSep::deconvolve(Matrix<float> &specmat,
1031                               const vector<double> shiftvec,
1032                               const double threshold,
1033                               Matrix<float> &outmat)
1034{
1035  LogIO os(LogOrigin("STSideBandSep","deconvolve()", WHERE));
1036  if (specmat.nrow() != nchan_)
1037    throw(AipsError("Internal error. The length of input matrix differs from nchan_"));
1038  if (specmat.ncolumn() != shiftvec.size())
1039    throw(AipsError("Internal error. The number of input shifts and spectrum  differs."));
1040
1041#ifdef KS_DEBUG
1042  float minval, maxval;
1043#endif
1044#ifdef KS_DEBUG
1045  minMax(minval, maxval, specmat);
1046  cout << "Max/Min of input Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
1047#endif
1048
1049  uInt ninsp = shiftvec.size();
1050  outmat.resize(nchan_, ninsp*(ninsp-1)/2, 0.);
1051  Matrix<Complex> fftspmat(nchan_/2+1, ninsp, 0.);
1052  Vector<float> rvecref(nchan_, 0.);
1053  Vector<Complex> cvecref(nchan_/2+1, 0.);
1054  uInt icol = 0;
1055  unsigned int nreject = 0;
1056
1057#ifdef KS_DEBUG
1058  cout << "Starting initial FFT. The number of input spectra = " << ninsp << endl;
1059  cout << "out matrix has ncolumn = " << outmat.ncolumn() << endl;
1060#endif
1061
1062  for (uInt isp = 0 ; isp < ninsp ; isp++) {
1063    rvecref.reference( specmat.column(isp) );
1064    cvecref.reference( fftspmat.column(isp) );
1065
1066#ifdef KS_DEBUG
1067    minMax(minval, maxval, rvecref);
1068    cout << "Max/Min of inv FFTed Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
1069#endif
1070
1071    fftsf.fft0(cvecref, rvecref, true);
1072
1073#ifdef KS_DEBUG
1074    double maxr=cvecref[0].real(), minr=cvecref[0].real(),
1075      maxi=cvecref[0].imag(), mini=cvecref[0].imag();
1076    for (uInt i = 1; i<cvecref.size();i++){
1077      maxr = max(maxr, cvecref[i].real());
1078      maxi = max(maxi, cvecref[i].imag());
1079      minr = min(minr, cvecref[i].real());
1080      mini = min(mini, cvecref[i].imag());
1081    }
1082    cout << "Max/Min of inv FFTed Matrix (size=" << cvecref.size() << ") = (max: " << maxr << " + " << maxi << "j , min: " << minr << " + " << mini << "j)" << endl;
1083#endif
1084  }
1085
1086  //Liberate from reference
1087  rvecref.unique();
1088
1089  Vector<Complex> cspec(nchan_/2+1, 0.);
1090  const double PI = 6.0 * asin(0.5);
1091  const double nchani = 1. / (float) nchan_;
1092  const Complex trans(0., 1.);
1093#ifdef KS_DEBUG
1094  cout << "starting actual deconvolution" << endl;
1095#endif
1096  for (uInt j = 0 ; j < ninsp ; j++) {
1097    for (uInt k = j+1 ; k < ninsp ; k++) {
1098      const double dx = (shiftvec[k] - shiftvec[j]) * 2. * PI * nchani;
1099
1100#ifdef KS_DEBUG
1101      cout << "icol = " << icol << endl;
1102#endif
1103
1104      for (uInt ichan = 0 ; ichan < cspec.size() ; ichan++){
1105        cspec[ichan] = ( fftspmat(ichan, j) + fftspmat(ichan, k) )*0.5;
1106        double phase = dx*ichan;
1107        if ( fabs( sin(phase) ) > threshold){
1108          cspec[ichan] += ( fftspmat(ichan, j) - fftspmat(ichan, k) ) * 0.5
1109            * trans * sin(phase) / ( 1. - cos(phase) );
1110        } else {
1111          nreject++;
1112        }
1113      } // end of channel loop
1114
1115#ifdef KS_DEBUG
1116      cout << "done calculation of cspec" << endl;
1117#endif
1118
1119      Vector<Float> rspec;
1120      rspec.reference( outmat.column(icol) );
1121
1122#ifdef KS_DEBUG
1123      cout << "Starting inverse FFT. icol = " << icol << endl;
1124      //cout << "- size of complex vector = " << cspec.size() << endl;
1125      double maxr=cspec[0].real(), minr=cspec[0].real(),
1126        maxi=cspec[0].imag(), mini=cspec[0].imag();
1127      for (uInt i = 1; i<cspec.size();i++){
1128        maxr = max(maxr, cspec[i].real());
1129        maxi = max(maxi, cspec[i].imag());
1130        minr = min(minr, cspec[i].real());
1131        mini = min(mini, cspec[i].imag());
1132      }
1133      cout << "Max/Min of conv vector (size=" << cspec.size() << ") = (max: " << maxr << " + " << maxi << "j , min: " << minr << " + " << mini << "j)" << endl;
1134#endif
1135
1136      fftsi.fft0(rspec, cspec, false);
1137
1138#ifdef KS_DEBUG
1139      //cout << "- size of inversed real vector = " << rspec.size() << endl;
1140      minMax(minval, maxval, rspec);
1141      cout << "Max/Min of inv FFTed Vector (size=" << rspec.size() << ") = (max: " << maxval << ", min: " << minval << ")" << endl;
1142      //cout << "Done inverse FFT. icol = " << icol << endl;
1143#endif
1144
1145      icol++;
1146    }
1147  }
1148
1149#ifdef KS_DEBUG
1150  minMax(minval, maxval, outmat);
1151  cout << "Max/Min of inv FFTed Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
1152#endif
1153
1154  os << "Threshold = " << threshold << ", Rejected channels = " << nreject << endl;
1155};
1156
1157
1158void STSideBandSep::aggregateMat(Matrix<float> &inmat,
1159                                 vector<float> &outvec)
1160{
1161  LogIO os(LogOrigin("STSideBandSep","aggregateMat()", WHERE));
1162  if (inmat.nrow() != nchan_)
1163    throw(AipsError("Internal error. The row numbers of input matrix differs from nchan_"));
1164//   if (outvec.size() != nchan_)
1165//     throw(AipsError("Internal error. The size of output vector should be equal to nchan_"));
1166
1167  os << "Averaging " << inmat.ncolumn() << " spectra in the input matrix."
1168     << LogIO::POST;
1169
1170  const uInt nspec = inmat.ncolumn();
1171  const double scale = 1./( (double) nspec );
1172  // initialize values with 0
1173  outvec.assign(nchan_, 0);
1174  for (uInt isp = 0 ; isp < nspec ; isp++) {
1175    for (uInt ich = 0 ; ich < nchan_ ; ich++) {
1176      outvec[ich] += inmat(ich, isp);
1177    }
1178  }
1179
1180  vector<float>::iterator iter;
1181  for (iter = outvec.begin(); iter != outvec.end(); iter++){
1182    *iter *= scale;
1183  }
1184};
1185
1186void STSideBandSep::subtractFromOther(const Matrix<float> &shiftmat,
1187                                      const vector<float> &invec,
1188                                      const vector<double> &shift,
1189                                      vector<float> &outvec)
1190{
1191  LogIO os(LogOrigin("STSideBandSep","subtractFromOther()", WHERE));
1192  if (shiftmat.nrow() != nchan_)
1193    throw(AipsError("Internal error. The row numbers of input matrix differs from nchan_"));
1194  if (invec.size() != nchan_)
1195    throw(AipsError("Internal error. The length of input vector should be nchan_"));
1196  if (shift.size() != shiftmat.ncolumn())
1197    throw(AipsError("Internal error. The column numbers of input matrix != the number of elements in shift"));
1198
1199  const uInt nspec = shiftmat.ncolumn();
1200  Vector<float> subsp(nchan_, 0.), shiftsub;
1201  Matrix<float> submat(nchan_, nspec, 0.);
1202  vector<float>::iterator iter;
1203  for (uInt isp = 0 ; isp < nspec ; isp++) {
1204    for (uInt ich = 0; ich < nchan_ ; ich++) {
1205      subsp(ich) = shiftmat(ich, isp) - invec[ich];
1206    }
1207    shiftsub.reference(submat.column(isp));
1208    shiftSpectrum(subsp, shift[isp], shiftsub);
1209  }
1210
1211  aggregateMat(submat, outvec);
1212};
1213
1214
1215void STSideBandSep::setLO1(const string lo1, const string frame,
1216                           const double reftime, const string refdir)
1217{
1218  Quantum<Double> qfreq;
1219  readQuantity(qfreq, String(lo1));
1220  lo1Freq_ = qfreq.getValue("Hz");
1221  MFrequency::getType(loFrame_, frame);
1222  loTime_ = reftime;
1223  loDir_ = refdir;
1224
1225#ifdef KS_DEBUG
1226  cout << "STSideBandSep::setLO1" << endl;
1227  if (lo1Freq_ > 0.)
1228    cout << "lo1 = " << lo1Freq_ << " [Hz] (" << frame << ")" << endl;
1229  if (loTime_ > 0.)
1230    cout << "ref time = " << loTime_ << " [day]" << endl;
1231  if (!loDir_.empty())
1232    cout << "ref direction = " << loDir_ << " [day]" << endl;
1233#endif
1234};
1235
1236void STSideBandSep::setLO1Root(string name)
1237{
1238   LogIO os(LogOrigin("STSideBandSep","setLO1Root()", WHERE));
1239   os << "Searching for '" << name << "'..." << LogIO::POST;
1240  // Check for existance of the file
1241  if (!checkFile(name)) {
1242     throw(AipsError("File does not exist"));
1243  }
1244  if (name[(name.size()-1)] == '/')
1245    name = name.substr(0,(name.size()-2));
1246
1247  if (checkFile(name+"/Receiver.xml", "file") &&
1248      checkFile(name+"/SpectralWindow.xml", "file")){
1249    os << "Found '" << name << "/Receiver.xml' ... got an ASDM name." << LogIO::POST;
1250    asdmName_ = name;
1251  } else if (checkFile(name+"/ASDM_RECEIVER") &&
1252             checkFile(name+"/ASDM_SPECTRALWINDOW")){
1253    os << "Found '" << name << "/ASDM_RECEIVER' ... got a Table name." << LogIO::POST;
1254    asisName_ = name;
1255  } else {
1256    throw(AipsError("Invalid file name. Set an MS or ASDM name."));
1257  }
1258
1259#ifdef KS_DEBUG
1260  cout << "STSideBandSep::setLO1Root" << endl;
1261  if (!asdmName_.empty())
1262    cout << "asdm name = " << asdmName_ << endl;
1263  if (!asisName_.empty())
1264    cout << "MS name = " << asisName_ << endl;
1265#endif
1266};
1267
1268
1269void STSideBandSep::solveImageFrequency()
1270{
1271  LogIO os(LogOrigin("STSideBandSep","solveImageFreqency()", WHERE));
1272  os << "Start calculating frequencies of image side band" << LogIO::POST;
1273
1274  if (imgTab_p.null())
1275    throw AipsError("STSideBandSep::solveImageFreqency - an image side band scantable should be set first");
1276
1277  // Convert frequency REFVAL to the value in frame of LO.
1278  // The code assumes that imgTab_p has only an IF and only a FREQ_ID
1279  // is associated to an IFNO
1280  // TODO: More complete Procedure would be
1281  // 1. Get freq IDs associated to sigIfno_
1282  // 2. Get freq information of the freq IDs
1283  // 3. For each freqIDs, get freq infromation in TOPO and an LO1
1284  //    frequency and calculate image band frequencies.
1285  STFrequencies freqTab_ = imgTab_p->frequencies();
1286  // get the base frame of table
1287  const MFrequency::Types tabframe = freqTab_.getFrame(true);
1288  TableVector<uInt> freqIdVec( imgTab_p->table(), "FREQ_ID" );
1289  // assuming single freqID per IFNO
1290  uInt freqid = freqIdVec(0);
1291  int nChan = imgTab_p->nchan(imgTab_p->getIF(0));
1292  double refpix, refval, increment ;
1293  freqTab_.getEntry(refpix, refval, increment, freqid);
1294  //MFrequency sigrefval = MFrequency(MVFrequency(refval),tabframe);
1295  // get freq infromation of sigIfno_ in loFrame_
1296  const MPosition mp = imgTab_p->getAntennaPosition();
1297  MEpoch me;
1298  MDirection md;
1299  if (loTime_ < 0.)
1300    me = imgTab_p->getEpoch(-1);
1301  else
1302    me = MEpoch(MVEpoch(loTime_));
1303  if (loDir_.empty()) {
1304    ArrayColumn<Double> srcdirCol_;
1305    srcdirCol_.attach(imgTab_p->table(), "SRCDIRECTION");
1306    // Assuming J2000 and SRCDIRECTION in unit of rad
1307    Quantum<Double> srcra = Quantum<Double>(srcdirCol_(0)(IPosition(1,0)), "rad");
1308    Quantum<Double> srcdec = Quantum<Double>(srcdirCol_(0)(IPosition(1,1)), "rad");
1309    md = MDirection(srcra, srcdec, MDirection::J2000);
1310    //imgTab_p->getDirection(0);
1311  } else {
1312    // parse direction string
1313    string::size_type pos0 = loDir_.find(" ");
1314   
1315    if (pos0 == string::npos) {
1316      throw AipsError("bad string format in LO1 direction");
1317    }
1318    string::size_type pos1 = loDir_.find(" ", pos0+1);
1319    String sepoch, sra, sdec;
1320    if (pos1 != string::npos) {
1321      sepoch = loDir_.substr(0, pos0);
1322      sra = loDir_.substr(pos0+1, pos1-pos0);
1323      sdec = loDir_.substr(pos1+1);
1324    }
1325    MDirection::Types epoch;
1326    MDirection::getType(epoch, sepoch);
1327    QuantumHolder qh ;
1328    String err ;
1329    qh.fromString( err, sra);
1330    Quantum<Double> ra = qh.asQuantumDouble() ;
1331    qh.fromString( err, sdec ) ;
1332    Quantum<Double> dec = qh.asQuantumDouble() ;
1333    //md = MDirection(ra.getValue("rad"), dec.getValue("rad"),epoch);
1334    md = MDirection(ra, dec, epoch);
1335  }
1336  MeasFrame mframe( me, mp, md );
1337  MFrequency::Convert tobframe(loFrame_, MFrequency::Ref(tabframe, mframe));
1338  MFrequency::Convert toloframe(tabframe, MFrequency::Ref(loFrame_, mframe));
1339  // Convert refval to loFrame_
1340  double sigrefval;
1341  if (tabframe == loFrame_)
1342    sigrefval = refval;
1343  else
1344    sigrefval = toloframe(Quantum<Double>(refval, "Hz")).get("Hz").getValue();
1345
1346  // Check for the availability of LO1
1347  if (lo1Freq_ > 0.) {
1348    os << "Using user defined LO1 frequency." << LogIO::POST;
1349  } else if (!asisName_.empty()) {
1350      // MS name is set.
1351    os << "Using user defined MS (asis): " << asisName_ << LogIO::POST;
1352    if (!getLo1FromAsisTab(asisName_, sigrefval, refpix, increment, nChan)) {
1353      throw AipsError("Failed to get LO1 frequency from MS");
1354    }
1355  } else if (!asdmName_.empty()) {
1356      // ASDM name is set.
1357    os << "Using user defined ASDM: " << asdmName_ << LogIO::POST;
1358    if (!getLo1FromAsdm(asdmName_, sigrefval, refpix, increment, nChan)) {
1359      throw AipsError("Failed to get LO1 frequency from ASDM");
1360    }
1361  } else {
1362    // Try getting ASDM name from scantable header
1363    os << "Try getting information from scantable header" << LogIO::POST;
1364    if (!getLo1FromScanTab(tableList_[0], sigrefval, refpix, increment, nChan)) {
1365      //throw AipsError("Failed to get LO1 frequency from asis table");
1366      os << LogIO::WARN << "Failed to get LO1 frequency using information in scantable." << LogIO::POST;
1367      os << LogIO::WARN << "Could not fill frequency information of IMAGE sideband properly." << LogIO::POST;
1368      os << LogIO::WARN << "Storing values of SIGNAL sideband in FREQUENCIES table" << LogIO::POST;
1369      return;
1370    }
1371  }
1372
1373  // LO1 should now be ready.
1374  if (lo1Freq_ < 0.)
1375    throw(AipsError("Got negative LO1 Frequency"));
1376
1377  // Print summary
1378  {
1379    // LO1
1380    Vector<Double> dirvec = md.getAngle(Unit(String("rad"))).getValue();
1381    os << "[LO1 settings]" << LogIO::POST;
1382    os << "- Frequency: " << lo1Freq_ << " [Hz] ("
1383       << MFrequency::showType(loFrame_) << ")" << LogIO::POST;
1384    os << "- Reference time: " << me.get(Unit(String("d"))).getValue()
1385       << " [day]" << LogIO::POST;
1386    os << "- Reference direction: [" << dirvec[0] << ", " << dirvec[1]
1387       << "] (" << md.getRefString() << ") " << LogIO::POST;
1388
1389    // signal sideband
1390    os << "[Signal side band]" << LogIO::POST;
1391    os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqid << ")"
1392       << LogIO::POST;
1393    os << "- Reference value: " << refval << " [Hz] ("
1394       << MFrequency::showType(tabframe) << ") = "
1395       << sigrefval << " [Hz] (" <<  MFrequency::showType(loFrame_)
1396       << ")" << LogIO::POST;
1397    os << "- Reference pixel: " << refpix  << LogIO::POST;
1398    os << "- Increment: " << increment << " [Hz]" << LogIO::POST;
1399  }
1400
1401  // Calculate image band incr and refval in loFrame_
1402  Double imgincr = -increment;
1403  Double imgrefval = 2 * lo1Freq_ - sigrefval;
1404  Double imgrefval_tab = imgrefval;
1405  // Convert imgrefval back to table base frame
1406  if (tabframe != loFrame_)
1407    imgrefval = tobframe(Quantum<Double>(imgrefval, "Hz")).get("Hz").getValue();
1408  // Set new frequencies table
1409  uInt fIDnew = freqTab_.addEntry(refpix, imgrefval, imgincr);
1410  // Update FREQ_ID in table.
1411  freqIdVec = fIDnew;
1412
1413  // Print summary (Image sideband)
1414  {
1415    os << "[Image side band]" << LogIO::POST;
1416    os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqIdVec(0)
1417       << ")" << LogIO::POST;
1418    os << "- Reference value: " << imgrefval << " [Hz] ("
1419       << MFrequency::showType(tabframe) << ") = "
1420       << imgrefval_tab << " [Hz] (" <<  MFrequency::showType(loFrame_)
1421       << ")" << LogIO::POST;
1422    os << "- Reference pixel: " << refpix  << LogIO::POST;
1423    os << "- Increment: " << imgincr << " [Hz]" << LogIO::POST;
1424  }
1425};
1426
1427Bool STSideBandSep::checkFile(const string name, string type)
1428{
1429  File file(name);
1430  if (!file.exists()){
1431    return false;
1432  } else if (type.empty()) {
1433    return true;
1434  } else {
1435    // Check for file type
1436    switch (std::tolower(type[0])) {
1437    case 'f':
1438      return file.isRegular(True);
1439    case 'd':
1440      return file.isDirectory(True);
1441    case 's':
1442      return file.isSymLink();
1443    default:
1444      throw AipsError("Invalid file type. Available types are 'file', 'directory', and 'symlink'.");
1445    }
1446  }
1447};
1448
1449bool STSideBandSep::getLo1FromAsdm(const string asdmname,
1450                                   const double /*refval*/,
1451                                   const double /*refpix*/,
1452                                   const double /*increment*/,
1453                                   const int /*nChan*/)
1454{
1455  // Check for relevant tables.
1456  string spwname = asdmname + "/SpectralWindow.xml";
1457  string recname = asdmname + "/Receiver.xml";
1458  if (!checkFile(spwname) || !checkFile(recname)) {
1459    throw(AipsError("Could not find subtables in ASDM"));
1460  }
1461
1462  return false;
1463
1464};
1465
1466
1467bool STSideBandSep::getLo1FromScanTab(CountedPtr< Scantable > &scantab,
1468                                      const double refval,
1469                                      const double refpix,
1470                                      const double increment,
1471                                      const int nChan)
1472{
1473  LogIO os(LogOrigin("STSideBandSep","getLo1FromScanTab()", WHERE));
1474  // Check for relevant tables.
1475  const TableRecord &rec = scantab->table().keywordSet() ;
1476  String spwname, recname;
1477  if (rec.isDefined("ASDM_SPECTRALWINDOW") && rec.isDefined("ASDM_RECEIVER")){
1478    spwname = rec.asString("ASDM_SPECTRALWINDOW");
1479    recname = rec.asString("ASDM_RECEIVER");
1480  }
1481  else {
1482    // keywords are not there
1483    os << LogIO::WARN
1484       << "Could not find necessary table names in scantable header."
1485       << LogIO::POST;
1486    return false;
1487  }
1488  if (!checkFile(spwname,"directory") || !checkFile(recname,"directory")) {
1489    throw(AipsError("Could not find relevant subtables in MS"));
1490  }
1491  // Get root MS name
1492  string msname;
1493  const String recsuff = "/ASDM_RECEIVER";
1494  String::size_type pos;
1495  pos = recname.size()-recsuff.size();
1496  if (recname.substr(pos) == recsuff)
1497    msname = recname.substr(0, pos);
1498  else
1499    throw(AipsError("Internal error in parsing table name from a scantable keyword."));
1500
1501  if (!checkFile(msname))
1502    throw(AipsError("Internal error in parsing MS name from a scantable keyword."));
1503
1504  return getLo1FromAsisTab(msname, refval, refpix, increment, nChan);
1505
1506};
1507
1508bool STSideBandSep::getLo1FromAsisTab(const string msname,
1509                                      const double refval,
1510                                      const double refpix,
1511                                      const double increment,
1512                                      const int nChan)
1513{
1514  LogIO os(LogOrigin("STSideBandSep","getLo1FromAsisTab()", WHERE));
1515  os << "Searching an LO1 frequency in '" << msname << "'" << LogIO::POST;
1516  // Check for relevant tables.
1517  const string spwname = msname + "/ASDM_SPECTRALWINDOW";
1518  const string recname = msname + "/ASDM_RECEIVER";
1519  if (!checkFile(spwname,"directory") || !checkFile(recname,"directory")) {
1520    throw(AipsError("Could not find relevant tables in MS"));
1521  }
1522
1523  Table spwtab_ = Table(spwname);
1524  String asdmSpw;
1525  ROTableRow spwrow(spwtab_);
1526  const Double rtol = 0.01;
1527  for (uInt idx = 0; idx < spwtab_.nrow(); idx++){
1528    const TableRecord& rec = spwrow.get(idx);
1529    // Compare nchan
1530    if (rec.asInt("numChan") != (Int) nChan)
1531      continue;
1532    // Compare increment
1533    Double asdminc;
1534    Array<Double> incarr = rec.asArrayDouble("chanWidthArray");
1535    if (incarr.size() > 0)
1536      asdminc = incarr(IPosition(1, (uInt) refpix));
1537    else
1538      asdminc = rec.asDouble("chanWidth");
1539    if (abs(asdminc - abs(increment)) > rtol * abs(increment))
1540      continue;
1541    // Compare refval
1542    Double asdmrefv;
1543    Array<Double> refvarr = rec.asArrayDouble("chanFreqArray");
1544    if (refvarr.size() > 0){
1545      const uInt iref = (uInt) refpix;
1546      const Double ratio = refpix - (Double) iref;
1547      asdmrefv = refvarr(IPosition(1, iref))*(1.-ratio)
1548        + refvarr(IPosition(1,iref+1))*ratio;
1549    }
1550    else {
1551      const Double ch0 = rec.asDouble("chanFreqStart");
1552      const Double chstep = rec.asDouble("chanFreqStep");
1553      asdmrefv = ch0 + chstep * refpix;
1554    }
1555    if (abs(asdmrefv - refval) < 0.5*abs(asdminc)){
1556      asdmSpw = rec.asString("spectralWindowId");
1557      break;
1558    }
1559  }
1560
1561  if (asdmSpw.empty()){
1562    os << LogIO::WARN << "Could not find relevant SPW ID in " << spwname << LogIO::POST;
1563    return false;
1564  }
1565  else {
1566    os << asdmSpw << " in " << spwname
1567       << " matches the freqeuncies of signal side band." << LogIO::POST;
1568  }
1569
1570  Table rectab_ = Table(recname);
1571  ROTableRow recrow(rectab_);
1572  for (uInt idx = 0; idx < rectab_.nrow(); idx++){
1573    const TableRecord& rec = recrow.get(idx);
1574    if (rec.asString("spectralWindowId") == asdmSpw){
1575      const Array<Double> loarr = rec.asArrayDouble("freqLO");
1576      lo1Freq_ = loarr(IPosition(1,0));
1577      os << "Found LO1 Frequency in " << recname << ": "
1578         << lo1Freq_ << " [Hz]" << LogIO::POST;
1579      return true;
1580    }
1581  }
1582  os << LogIO::WARN << "Could not find " << asdmSpw << " in " << recname
1583     << LogIO::POST;
1584  return false;
1585};
1586
1587} //namespace asap
Note: See TracBrowser for help on using the repository browser.