source: trunk/src/SDMemTable.cc @ 465

Last change on this file since 465 was 465, checked in by mar637, 19 years ago

Added SDFitTable to handle fits and expose them to python vi the sdfit class.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 49.9 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDMemTable.cc: A MemoryTable container for single dish integrations
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
32#include <map>
33
34#include <casa/aips.h>
35#include <casa/iostream.h>
36#include <casa/iomanip.h>
37#include <casa/Arrays/Array.h>
38#include <casa/Arrays/ArrayMath.h>
39#include <casa/Arrays/MaskArrMath.h>
40#include <casa/Arrays/ArrayLogical.h>
41#include <casa/Arrays/ArrayAccessor.h>
42#include <casa/Arrays/VectorSTLIterator.h>
43#include <casa/Arrays/Vector.h>
44#include <casa/BasicMath/Math.h>
45#include <casa/Quanta/MVAngle.h>
46
47#include <tables/Tables/TableParse.h>
48#include <tables/Tables/TableDesc.h>
49#include <tables/Tables/SetupNewTab.h>
50#include <tables/Tables/ScaColDesc.h>
51#include <tables/Tables/ArrColDesc.h>
52
53#include <tables/Tables/ExprNode.h>
54#include <tables/Tables/TableRecord.h>
55#include <measures/Measures/MFrequency.h>
56#include <measures/Measures/MeasTable.h>
57#include <coordinates/Coordinates/CoordinateUtil.h>
58#include <casa/Quanta/MVTime.h>
59#include <casa/Quanta/MVAngle.h>
60
61#include "SDDefs.h"
62#include "SDContainer.h"
63#include "MathUtils.h"
64#include "SDPol.h"
65
66#include "SDMemTable.h"
67
68using namespace casa;
69using namespace asap;
70
71SDMemTable::SDMemTable() :
72  IFSel_(0),
73  beamSel_(0),
74  polSel_(0)
75{
76  setup();
77  attach();
78}
79
80SDMemTable::SDMemTable(const std::string& name) :
81  IFSel_(0),
82  beamSel_(0),
83  polSel_(0)
84{
85  Table tab(name);
86  table_ = tab.copyToMemoryTable("dummy");
87  //cerr << "hello from C SDMemTable @ " << this << endl;
88  attach();
89}
90
91SDMemTable::SDMemTable(const SDMemTable& other, Bool clear)
92{
93  IFSel_= other.IFSel_;
94  beamSel_= other.beamSel_;
95  polSel_= other.polSel_;
96  chanMask_ = other.chanMask_;
97  table_ = other.table_.copyToMemoryTable(String("dummy"));
98  // clear all rows()
99  if (clear) {
100    table_.removeRow(this->table_.rowNumbers());
101  } else {
102    IFSel_ = other.IFSel_;
103    beamSel_ = other.beamSel_;
104    polSel_ = other.polSel_;
105  }
106//
107  attach();
108  //cerr << "hello from CC SDMemTable @ " << this << endl;
109}
110
111SDMemTable::SDMemTable(const Table& tab, const std::string& exprs) :
112  IFSel_(0),
113  beamSel_(0),
114  polSel_(0)
115{
116  Table t = tableCommand(exprs,tab);
117  if (t.nrow() == 0)
118      throw(AipsError("Query unsuccessful."));
119  table_ = t.copyToMemoryTable("dummy");
120  attach();
121  renumber();
122}
123
124SDMemTable::~SDMemTable()
125{
126  //cerr << "goodbye from SDMemTable @ " << this << endl;
127}
128
129SDMemTable SDMemTable::getScan(Int scanID) const
130{
131  String cond("SELECT * from $1 WHERE SCANID == ");
132  cond += String::toString(scanID);
133  return SDMemTable(table_, cond);
134}
135
136SDMemTable &SDMemTable::operator=(const SDMemTable& other)
137{
138  if (this != &other) {
139     IFSel_= other.IFSel_;
140     beamSel_= other.beamSel_;
141     polSel_= other.polSel_;
142     chanMask_.resize(0);
143     chanMask_ = other.chanMask_;
144     table_ = other.table_.copyToMemoryTable(String("dummy"));
145     attach();
146  }
147  //cerr << "hello from ASS SDMemTable @ " << this << endl;
148  return *this;
149}
150
151SDMemTable SDMemTable::getSource(const std::string& source) const
152{
153  String cond("SELECT * from $1 WHERE SRCNAME == ");
154  cond += source;
155  return SDMemTable(table_, cond);
156}
157
158void SDMemTable::setup()
159{
160  TableDesc td("", "1", TableDesc::Scratch);
161  td.comment() = "A SDMemTable";
162 
163  td.addColumn(ScalarColumnDesc<Double>("TIME"));
164  td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
165  td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
166  td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
167  td.addColumn(ArrayColumnDesc<Float>("TSYS"));
168  td.addColumn(ArrayColumnDesc<Float>("STOKES"));
169  td.addColumn(ScalarColumnDesc<Int>("SCANID"));
170  td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
171  td.addColumn(ArrayColumnDesc<uInt>("FREQID"));
172  td.addColumn(ArrayColumnDesc<uInt>("RESTFREQID"));
173  td.addColumn(ArrayColumnDesc<Double>("DIRECTION"));
174  td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
175  td.addColumn(ScalarColumnDesc<String>("TCALTIME"));
176  td.addColumn(ArrayColumnDesc<Float>("TCAL"));
177  td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
178  td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
179  td.addColumn(ScalarColumnDesc<Float>("PARANGLE"));
180  td.addColumn(ScalarColumnDesc<Int>("REFBEAM"));
181  td.addColumn(ArrayColumnDesc<String>("HISTORY"));
182  td.addColumn(ArrayColumnDesc<Int>("FITID"));
183
184
185  // Now create Table SetUp from the description.
186
187  SetupNewTable aNewTab("dummy", td, Table::New);
188
189  // Bind the Stokes Virtual machine to the STOKES column Because we
190  // don't know how many polarizations will be in the data at this
191  // point, we must bind the Virtual Engine regardless.  The STOKES
192  // column won't be accessed if not appropriate (nPol=4)
193
194
195   SDStokesEngine::registerClass();
196   SDStokesEngine stokesEngine(String("STOKES"), String("SPECTRA"));
197   aNewTab.bindColumn ("STOKES", stokesEngine);
198
199   // Create Table
200  table_ = Table(aNewTab, Table::Memory, 0);
201  // add subtable
202  TableDesc tdf("", "1", TableDesc::Scratch);
203  tdf.addColumn(ArrayColumnDesc<String>("FUNCTIONS"));
204  tdf.addColumn(ArrayColumnDesc<Int>("COMPONENTS"));
205  tdf.addColumn(ArrayColumnDesc<Double>("PARAMETERS"));
206  tdf.addColumn(ArrayColumnDesc<Bool>("PARMASK"));
207  tdf.addColumn(ArrayColumnDesc<String>("FRAMEINFO"));
208  SetupNewTable fittab("fits", tdf, Table::New);
209  Table fitTable(fittab, Table::Memory);
210  table_.rwKeywordSet().defineTable("FITS", fitTable);
211
212
213}
214
215void SDMemTable::attach()
216{
217  timeCol_.attach(table_, "TIME");
218  srcnCol_.attach(table_, "SRCNAME");
219  specCol_.attach(table_, "SPECTRA");
220  flagsCol_.attach(table_, "FLAGTRA");
221  tsCol_.attach(table_, "TSYS");
222  stokesCol_.attach(table_, "STOKES");
223  scanCol_.attach(table_, "SCANID");
224  integrCol_.attach(table_, "INTERVAL");
225  freqidCol_.attach(table_, "FREQID");
226  restfreqidCol_.attach(table_, "RESTFREQID");
227  dirCol_.attach(table_, "DIRECTION");
228  fldnCol_.attach(table_, "FIELDNAME");
229  tcaltCol_.attach(table_, "TCALTIME");
230  tcalCol_.attach(table_, "TCAL");
231  azCol_.attach(table_, "AZIMUTH");
232  elCol_.attach(table_, "ELEVATION");
233  paraCol_.attach(table_, "PARANGLE");
234  rbeamCol_.attach(table_, "REFBEAM");
235  histCol_.attach(table_, "HISTORY");
236  fitCol_.attach(table_,"FITID");
237}
238
239
240std::string SDMemTable::getSourceName(Int whichRow) const
241{
242  String name;
243  srcnCol_.get(whichRow, name);
244  return name;
245}
246
247std::string SDMemTable::getTime(Int whichRow, Bool showDate) const
248{
249  Double tm;
250  if (whichRow > -1) {
251    timeCol_.get(whichRow, tm);
252  } else {
253    table_.keywordSet().get("UTC",tm);
254  }
255  MVTime mvt(tm);
256  if (showDate)
257    mvt.setFormat(MVTime::YMD);
258  else
259    mvt.setFormat(MVTime::TIME);
260  ostringstream oss;
261  oss << mvt;
262  return String(oss);
263}
264
265double SDMemTable::getInterval(Int whichRow) const
266{
267  Double intval;
268  integrCol_.get(whichRow, intval);
269  return intval;
270}
271
272bool SDMemTable::setIF(Int whichIF)
273{
274  if ( whichIF >= 0 && whichIF < nIF()) {
275    IFSel_ = whichIF;
276    return true;
277  }
278  return false;
279}
280
281bool SDMemTable::setBeam(Int whichBeam)
282{
283  if ( whichBeam >= 0 && whichBeam < nBeam()) {
284    beamSel_ = whichBeam;
285    return true;
286  }
287  return false;
288}
289
290bool SDMemTable::setPol(Int whichPol)
291{
292  if ( whichPol >= 0 && whichPol < nPol()) {
293    polSel_ = whichPol;
294    return true;
295  }
296  return false;
297}
298
299void SDMemTable::resetCursor()
300{
301   polSel_ = 0;
302   IFSel_ = 0;
303   beamSel_ = 0;
304}
305
306bool SDMemTable::setMask(std::vector<int> whichChans)
307{
308  std::vector<int>::iterator it;
309  uInt n = flagsCol_.shape(0)(3);
310  if (whichChans.empty()) {
311    chanMask_ = std::vector<bool>(n,true);
312    return true;     
313  }
314  chanMask_.resize(n,true);
315  for (it = whichChans.begin(); it != whichChans.end(); ++it) {
316    if (*it < n) {
317      chanMask_[*it] = false;
318    }
319  }
320  return true;
321}
322
323std::vector<bool> SDMemTable::getMask(Int whichRow) const
324{
325
326  std::vector<bool> mask;
327
328  Array<uChar> arr;
329  flagsCol_.get(whichRow, arr);
330
331  ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
332  aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
333  ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
334  aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
335  ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
336  aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
337
338  Bool useUserMask = ( chanMask_.size() == arr.shape()(3) );
339
340  std::vector<bool> tmp;
341  tmp = chanMask_; // WHY the fxxx do I have to make a copy here
342  std::vector<bool>::iterator miter;
343  miter = tmp.begin();
344
345  for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
346    bool out =!static_cast<bool>(*i);
347    if (useUserMask) {
348      out = out && (*miter);
349      miter++;
350    }
351    mask.push_back(out);
352  }
353  return mask;
354}
355
356
357
358std::vector<float> SDMemTable::getSpectrum(Int whichRow) const
359{
360
361  Array<Float> arr;
362  specCol_.get(whichRow, arr);
363  return getFloatSpectrum(arr);
364}
365
366std::vector<float> SDMemTable::getStokesSpectrum(Int whichRow, Bool doPol,
367                                                 Float paOffset) const
368  // Gets
369  //  doPol=False  : I,Q,U,V
370  //  doPol=True   : I,P,PA,V   ; P = sqrt(Q**2+U**2), PA = 0.5*atan2(Q,U)
371  //
372{
373  AlwaysAssert(asap::nAxes==4,AipsError);
374  if (nPol()!=1 && nPol()!=2 && nPol()!=4) {
375     throw (AipsError("You must have 1,2 or 4 polarizations to get the Stokes parameters"));
376  }
377  Array<Float> arr;
378  stokesCol_.get(whichRow, arr);
379
380  if (doPol && (polSel_==1 || polSel_==2)) {       //   Q,U --> P, P.A.
381
382    // Set current cursor location
383     const IPosition& shape = arr.shape();
384     IPosition start, end;
385     setCursorSlice (start, end, shape);
386
387     // Get Q and U slices
388
389     Array<Float> Q = SDPolUtil::getStokesSlice(arr,start,end,"Q");
390     Array<Float> U = SDPolUtil::getStokesSlice(arr,start,end,"U");
391
392     // Compute output
393
394     Array<Float> out;
395     if (polSel_==1) {                                        // P
396        out = SDPolUtil::polarizedIntensity(Q,U);
397     } else if (polSel_==2) {                                 // P.A.
398        out = SDPolUtil::positionAngle(Q,U) + paOffset;
399     }
400
401     // Copy to output
402
403     IPosition vecShape(1,shape(asap::ChanAxis));
404     Vector<Float> outV = out.reform(vecShape);
405     std::vector<float> stlout;
406     outV.tovector(stlout);
407     return stlout;
408
409  } else {
410
411    // Selects at the cursor location
412    return getFloatSpectrum(arr);
413  }
414}
415
416std::vector<float> SDMemTable::getCircularSpectrum(Int whichRow,
417                                                   Bool doRR) const
418  // Gets
419  //  RR = I + V
420  //  LL = I - V
421{
422  AlwaysAssert(asap::nAxes==4,AipsError);
423  if (nPol()!=4) {
424     throw (AipsError("You must have 4 polarizations to get RR or LL"));
425  }
426  Array<Float> arr;
427  stokesCol_.get(whichRow, arr);
428
429  // Set current cursor location
430
431  const IPosition& shape = arr.shape();
432  IPosition start, end;
433  setCursorSlice(start, end, shape);
434
435  // Get I and V slices
436
437  Array<Float> I = SDPolUtil::getStokesSlice(arr,start,end,"I");
438  Array<Float> V = SDPolUtil::getStokesSlice(arr,start,end,"V");
439
440  // Compute output
441
442  Array<Float> out = SDPolUtil::circularPolarizationFromStokes(I, V, doRR);
443
444  // Copy to output
445
446  IPosition vecShape(1,shape(asap::ChanAxis));
447  Vector<Float> outV = out.reform(vecShape);
448  std::vector<float> stlout;
449  outV.tovector(stlout);
450
451  return stlout;
452}
453
454
455std::vector<string> SDMemTable::getCoordInfo() const
456{
457  String un;
458  Table t = table_.keywordSet().asTable("FREQUENCIES");
459  String sunit;
460  t.keywordSet().get("UNIT",sunit);
461  String dpl;
462  t.keywordSet().get("DOPPLER",dpl);
463  if (dpl == "") dpl = "RADIO";
464  String rfrm;
465  t.keywordSet().get("REFFRAME",rfrm);
466  std::vector<string> inf;
467  inf.push_back(sunit);
468  inf.push_back(rfrm);
469  inf.push_back(dpl);
470  t.keywordSet().get("BASEREFFRAME",rfrm);
471  inf.push_back(rfrm);
472  return inf;
473}
474
475void SDMemTable::setCoordInfo(std::vector<string> theinfo)
476{
477  std::vector<string>::iterator it;
478  String un,rfrm, brfrm,dpl;
479  un = theinfo[0];              // Abcissa unit
480  rfrm = theinfo[1];            // Active (or conversion) frame
481  dpl = theinfo[2];             // Doppler
482  brfrm = theinfo[3];           // Base frame
483
484  Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
485
486  Vector<Double> rstf;
487  t.keywordSet().get("RESTFREQS",rstf);
488
489  Bool canDo = True;
490  Unit u1("km/s");Unit u2("Hz");
491  if (Unit(un) == u1) {
492    Vector<Double> rstf;
493    t.keywordSet().get("RESTFREQS",rstf);
494    if (rstf.nelements() == 0) {
495        throw(AipsError("Can't set unit to km/s if no restfrequencies are specified"));
496    }
497  } else if (Unit(un) != u2 && un != "") {
498        throw(AipsError("Unit not conformant with Spectral Coordinates"));
499  }
500  t.rwKeywordSet().define("UNIT", un);
501
502  MFrequency::Types mdr;
503  if (!MFrequency::getType(mdr, rfrm)) {
504   
505    Int a,b;const uInt* c;
506    const String* valid = MFrequency::allMyTypes(a, b, c);
507    String pfix = "Please specify a legal frame type. Types are\n";
508    throw(AipsError(pfix+(*valid)));
509  } else {
510    t.rwKeywordSet().define("REFFRAME",rfrm);
511  }
512
513  MDoppler::Types dtype;
514  dpl.upcase();
515  if (!MDoppler::getType(dtype, dpl)) {
516    throw(AipsError("Doppler type unknown"));
517  } else {
518    t.rwKeywordSet().define("DOPPLER",dpl);
519  }
520
521  if (!MFrequency::getType(mdr, brfrm)) {
522     Int a,b;const uInt* c;
523     const String* valid = MFrequency::allMyTypes(a, b, c);
524     String pfix = "Please specify a legal frame type. Types are\n";
525     throw(AipsError(pfix+(*valid)));
526   } else {
527    t.rwKeywordSet().define("BASEREFFRAME",brfrm);
528   }
529}
530
531
532std::vector<double> SDMemTable::getAbcissa(Int whichRow) const
533{
534  std::vector<double> abc(nChan());
535
536  // Get header units keyword
537  Table t = table_.keywordSet().asTable("FREQUENCIES");
538  String sunit;
539  t.keywordSet().get("UNIT",sunit);
540  if (sunit == "") sunit = "pixel";
541  Unit u(sunit);
542
543  // Easy if just wanting pixels
544  if (sunit==String("pixel")) {
545    // assume channels/pixels
546    std::vector<double>::iterator it;
547    uInt i=0;
548    for (it = abc.begin(); it != abc.end(); ++it) {
549      (*it) = Double(i++);
550    }
551    return abc;
552  }
553
554  // Continue with km/s or Hz.  Get FreqIDs
555  Vector<uInt> freqIDs;
556  freqidCol_.get(whichRow, freqIDs);
557  uInt freqID = freqIDs(IFSel_);
558  restfreqidCol_.get(whichRow, freqIDs);
559  uInt restFreqID = freqIDs(IFSel_);
560
561  // Get SpectralCoordinate, set reference frame conversion,
562  // velocity conversion, and rest freq state
563
564  SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow);
565  Vector<Double> pixel(nChan());
566  indgen(pixel);
567
568  if (u == Unit("km/s")) {
569     Vector<Double> world;
570     spc.pixelToVelocity(world,pixel);
571     std::vector<double>::iterator it;
572     uInt i = 0;
573     for (it = abc.begin(); it != abc.end(); ++it) {
574       (*it) = world[i];
575       i++;
576     }
577  } else if (u == Unit("Hz")) {
578
579    // Set world axis units   
580    Vector<String> wau(1); wau = u.getName();
581    spc.setWorldAxisUnits(wau);
582
583    std::vector<double>::iterator it;
584    Double tmp;
585    uInt i = 0;
586    for (it = abc.begin(); it != abc.end(); ++it) {
587      spc.toWorld(tmp,pixel[i]);
588      (*it) = tmp;
589      i++;
590    }
591  }
592  return abc;
593}
594
595std::string SDMemTable::getAbcissaString(Int whichRow) const
596{
597  Table t = table_.keywordSet().asTable("FREQUENCIES");
598
599  String sunit;
600  t.keywordSet().get("UNIT",sunit);
601  if (sunit == "") sunit = "pixel";
602  Unit u(sunit);
603
604  Vector<uInt> freqIDs;
605  freqidCol_.get(whichRow, freqIDs);
606  uInt freqID = freqIDs(IFSel_); 
607  restfreqidCol_.get(whichRow, freqIDs);
608  uInt restFreqID = freqIDs(IFSel_);
609
610  // Get SpectralCoordinate, with frame, velocity, rest freq state set
611  SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow);
612
613  String s = "Channel";
614  if (u == Unit("km/s")) {
615    s = CoordinateUtil::axisLabel(spc,0,True,True,True);
616  } else if (u == Unit("Hz")) {
617    Vector<String> wau(1);wau = u.getName();
618    spc.setWorldAxisUnits(wau);
619
620    s = CoordinateUtil::axisLabel(spc,0,True,True,False);
621  }
622  return s;
623}
624
625void SDMemTable::setSpectrum(std::vector<float> spectrum, int whichRow)
626{
627  Array<Float> arr;
628  specCol_.get(whichRow, arr);
629  if (spectrum.size() != arr.shape()(asap::ChanAxis)) {
630    throw(AipsError("Attempting to set spectrum with incorrect length."));
631  }
632
633  // Setup accessors
634  ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
635  aa0.reset(aa0.begin(uInt(beamSel_)));                   // Beam selection
636  ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
637  aa1.reset(aa1.begin(uInt(IFSel_)));                     // IF selection
638  ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
639  aa2.reset(aa2.begin(uInt(polSel_)));                    // Pol selection
640
641  // Iterate
642  std::vector<float>::iterator it = spectrum.begin();
643  for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
644    (*i) = Float(*it);
645    it++;
646  }
647  specCol_.put(whichRow, arr);
648}
649
650void SDMemTable::getSpectrum(Vector<Float>& spectrum, Int whichRow) const
651{
652  Array<Float> arr;
653  specCol_.get(whichRow, arr);
654
655  // Iterate and extract
656  spectrum.resize(arr.shape()(3));
657  ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
658  aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
659  ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
660  aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
661  ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
662  aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
663
664  ArrayAccessor<Float, Axis<asap::BeamAxis> > va(spectrum);
665  for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
666    (*va) = (*i);
667    va++;
668  }
669}
670
671
672/*
673void SDMemTable::getMask(Vector<Bool>& mask, Int whichRow) const {
674  Array<uChar> arr;
675  flagsCol_.get(whichRow, arr);
676  mask.resize(arr.shape()(3));
677
678  ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
679  aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
680  ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
681  aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
682  ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
683  aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
684
685  Bool useUserMask = ( chanMask_.size() == arr.shape()(3) );
686
687  ArrayAccessor<Bool, Axis<asap::BeamAxis> > va(mask);
688  std::vector<bool> tmp;
689  tmp = chanMask_; // WHY the fxxx do I have to make a copy here. The
690                   // iterator should work on chanMask_??
691  std::vector<bool>::iterator miter;
692  miter = tmp.begin();
693
694  for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
695    bool out =!static_cast<bool>(*i);
696    if (useUserMask) {
697      out = out && (*miter);
698      miter++;
699    }
700    (*va) = out;
701    va++;
702  }
703}
704*/
705
706MaskedArray<Float> SDMemTable::rowAsMaskedArray(uInt whichRow,
707                                                Bool toStokes) const
708{
709  // Get flags
710  Array<uChar> farr;
711  flagsCol_.get(whichRow, farr);
712
713  // Get data and convert mask to Bool
714  Array<Float> arr;
715  Array<Bool> mask;
716  if (toStokes) {
717     stokesCol_.get(whichRow, arr);
718
719     Array<Bool> tMask(farr.shape());
720     convertArray(tMask, farr);
721     mask = SDPolUtil::stokesData (tMask, True);
722  } else {
723     specCol_.get(whichRow, arr);
724     mask.resize(farr.shape());
725     convertArray(mask, farr);
726  }
727
728  return MaskedArray<Float>(arr,!mask);
729}
730
731Float SDMemTable::getTsys(Int whichRow) const
732{
733  Array<Float> arr;
734  tsCol_.get(whichRow, arr);
735  Float out;
736
737  IPosition ip(arr.shape());
738  ip(asap::BeamAxis) = beamSel_;
739  ip(asap::IFAxis) = IFSel_;
740  ip(asap::PolAxis) = polSel_;
741  ip(asap::ChanAxis)=0;               // First channel only
742
743  out = arr(ip);
744  return out;
745}
746
747MDirection SDMemTable::getDirection(Int whichRow, Bool refBeam) const
748{
749  MDirection::Types mdr = getDirectionReference();
750  Array<Double> posit;
751  dirCol_.get(whichRow,posit);
752  Vector<Double> wpos(2);
753  Int rb;
754  rbeamCol_.get(whichRow,rb);
755  wpos[0] = posit(IPosition(2,beamSel_,0));
756  wpos[1] = posit(IPosition(2,beamSel_,1));
757  if (refBeam && rb > -1) {  // use refBeam instead if it exists
758    wpos[0] = posit(IPosition(2,rb,0));
759    wpos[1] = posit(IPosition(2,rb,1));
760  }
761
762  Quantum<Double> lon(wpos[0],Unit(String("rad")));
763  Quantum<Double> lat(wpos[1],Unit(String("rad")));
764  return MDirection(lon, lat, mdr);
765}
766
767MEpoch SDMemTable::getEpoch(Int whichRow) const
768{
769  MEpoch::Types met = getTimeReference();
770
771  Double obstime;
772  timeCol_.get(whichRow,obstime);
773  MVEpoch tm2(Quantum<Double>(obstime, Unit(String("d"))));
774  return MEpoch(tm2, met);
775}
776
777MPosition SDMemTable::getAntennaPosition () const
778{
779  Vector<Double> antpos;
780  table_.keywordSet().get("AntennaPosition", antpos);
781  MVPosition mvpos(antpos(0),antpos(1),antpos(2));
782  return MPosition(mvpos);
783}
784
785
786SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID) const
787{
788 
789  Table t = table_.keywordSet().asTable("FREQUENCIES");
790  if (freqID> t.nrow() ) {
791    throw(AipsError("SDMemTable::getSpectralCoordinate - freqID out of range"));
792  }
793
794  Double rp,rv,inc;
795  String rf;
796  ROScalarColumn<Double> rpc(t, "REFPIX");
797  ROScalarColumn<Double> rvc(t, "REFVAL");
798  ROScalarColumn<Double> incc(t, "INCREMENT");
799  t.keywordSet().get("BASEREFFRAME",rf);
800
801  // Create SpectralCoordinate (units Hz)
802  MFrequency::Types mft;
803  if (!MFrequency::getType(mft, rf)) {
804    cerr << "Frequency type unknown assuming TOPO" << endl;
805    mft = MFrequency::TOPO;
806  }
807  rpc.get(freqID, rp);
808  rvc.get(freqID, rv);
809  incc.get(freqID, inc);
810
811  SpectralCoordinate spec(mft,rv,inc,rp);
812  return spec;
813}
814
815
816SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID,
817                                                     uInt restFreqID,
818                                                     uInt whichRow) const
819{
820 
821  // Create basic SC
822  SpectralCoordinate spec = getSpectralCoordinate (freqID);
823
824  Table t = table_.keywordSet().asTable("FREQUENCIES");
825
826  // Set rest frequency
827  Vector<Double> restFreqIDs;
828  t.keywordSet().get("RESTFREQS",restFreqIDs);
829  if (restFreqID < restFreqIDs.nelements()) {    // Should always be true
830    spec.setRestFrequency(restFreqIDs(restFreqID),True);
831  }
832
833  // Set up frame conversion layer
834  String frm;
835  t.keywordSet().get("REFFRAME",frm);
836  if (frm == "") frm = "TOPO";
837  MFrequency::Types mtype;
838  if (!MFrequency::getType(mtype, frm)) {
839    // Should never happen
840    cerr << "Frequency type unknown assuming TOPO" << endl;
841    mtype = MFrequency::TOPO;
842  }
843
844  // Set reference frame conversion  (requires row)
845  MDirection dir = getDirection(whichRow);
846  MEpoch epoch = getEpoch(whichRow);
847  MPosition pos = getAntennaPosition();
848
849  if (!spec.setReferenceConversion(mtype,epoch,pos,dir)) {
850    throw(AipsError("Couldn't convert frequency frame."));
851  }
852
853  // Now velocity conversion if appropriate
854  String unitStr;
855  t.keywordSet().get("UNIT",unitStr);
856
857  String dpl;
858  t.keywordSet().get("DOPPLER",dpl);
859  MDoppler::Types dtype;
860  MDoppler::getType(dtype, dpl);
861
862  // Only set velocity unit if non-blank and non-Hz
863  if (!unitStr.empty()) {
864     Unit unitU(unitStr);
865     if (unitU==Unit("Hz")) {
866     } else {
867        spec.setVelocity(unitStr, dtype);
868     }
869  }
870
871  return spec;
872}
873
874
875Bool SDMemTable::setCoordinate(const SpectralCoordinate& speccord,
876                               uInt freqID) {
877  Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
878  if (freqID > t.nrow() ) {
879    throw(AipsError("SDMemTable::setCoordinate - coord no out of range"));
880  }
881  ScalarColumn<Double> rpc(t, "REFPIX");
882  ScalarColumn<Double> rvc(t, "REFVAL");
883  ScalarColumn<Double> incc(t, "INCREMENT");
884
885  rpc.put(freqID, speccord.referencePixel()[0]);
886  rvc.put(freqID, speccord.referenceValue()[0]);
887  incc.put(freqID, speccord.increment()[0]);
888
889  return True;
890}
891
892Int SDMemTable::nCoordinates() const
893{
894  return table_.keywordSet().asTable("FREQUENCIES").nrow();
895}
896
897
898std::vector<double> SDMemTable::getRestFreqs() const
899{
900  Table t = table_.keywordSet().asTable("FREQUENCIES");
901  Vector<Double> tvec;
902  t.keywordSet().get("RESTFREQS",tvec);
903  std::vector<double> stlout;
904  tvec.tovector(stlout);
905  return stlout; 
906}
907
908bool SDMemTable::putSDFreqTable(const SDFrequencyTable& sdft)
909{
910  TableDesc td("", "1", TableDesc::Scratch);
911  td.addColumn(ScalarColumnDesc<Double>("REFPIX"));
912  td.addColumn(ScalarColumnDesc<Double>("REFVAL"));
913  td.addColumn(ScalarColumnDesc<Double>("INCREMENT"));
914  SetupNewTable aNewTab("freqs", td, Table::New);
915  Table aTable (aNewTab, Table::Memory, sdft.length());
916  ScalarColumn<Double> sc0(aTable, "REFPIX");
917  ScalarColumn<Double> sc1(aTable, "REFVAL");
918  ScalarColumn<Double> sc2(aTable, "INCREMENT");
919  for (uInt i=0; i < sdft.length(); ++i) {
920    sc0.put(i,sdft.referencePixel(i));
921    sc1.put(i,sdft.referenceValue(i));
922    sc2.put(i,sdft.increment(i));
923  }
924  String rf = sdft.refFrame();
925  if (rf.contains("TOPO")) rf = "TOPO";
926
927  aTable.rwKeywordSet().define("BASEREFFRAME", rf);
928  aTable.rwKeywordSet().define("REFFRAME", rf);
929  aTable.rwKeywordSet().define("EQUINOX", sdft.equinox());
930  aTable.rwKeywordSet().define("UNIT", String(""));
931  aTable.rwKeywordSet().define("DOPPLER", String("RADIO"));
932  Vector<Double> rfvec;
933  String rfunit;
934  sdft.restFrequencies(rfvec,rfunit);
935  Quantum<Vector<Double> > q(rfvec, rfunit);
936  rfvec.resize();
937  rfvec = q.getValue("Hz");
938  aTable.rwKeywordSet().define("RESTFREQS", rfvec);
939  table_.rwKeywordSet().defineTable("FREQUENCIES", aTable);
940  return true;
941}
942
943bool SDMemTable::putSDFitTable(const SDFitTable& sdft)
944{
945  TableDesc td("", "1", TableDesc::Scratch);
946  td.addColumn(ArrayColumnDesc<String>("FUNCTIONS"));
947  td.addColumn(ArrayColumnDesc<Int>("COMPONENTS"));
948  td.addColumn(ArrayColumnDesc<Double>("PARAMETERS"));
949  td.addColumn(ArrayColumnDesc<Bool>("PARMASK"));
950  td.addColumn(ArrayColumnDesc<String>("FRAMEINFO"));
951  SetupNewTable aNewTab("fits", td, Table::New);
952  Table aTable(aNewTab, Table::Memory);
953  ArrayColumn<String> sc0(aTable, "FUNCTIONS");
954  ArrayColumn<Int> sc1(aTable, "COMPONENTS");
955  ArrayColumn<Double> sc2(aTable, "PARAMETERS");
956  ArrayColumn<Bool> sc3(aTable, "PARMASK");
957  ArrayColumn<String> sc4(aTable, "FRAMEINFO");
958  for (uInt i; i<sdft.length(); ++i) {
959    const Vector<Double>& parms = sdft.getParameters(i);
960    const Vector<Bool>& parmask = sdft.getParameterMask(i);
961    const Vector<String>& funcs = sdft.getFunctions(i);
962    const Vector<Int>& comps = sdft.getComponents(i);
963    const Vector<String>& finfo = sdft.getFrameInfo(i);
964    sc0.put(i,funcs);
965    sc1.put(i,comps);
966    sc3.put(i,parmask);
967    sc2.put(i,parms);
968    sc4.put(i,finfo);
969  }
970  table_.rwKeywordSet().defineTable("FITS", aTable);
971  return true;
972}
973
974SDFitTable SDMemTable::getSDFitTable() const
975{
976  const Table& t = table_.keywordSet().asTable("FITS");
977  Vector<Double> parms;
978  Vector<Bool> parmask;
979  Vector<String> funcs;
980  Vector<String> finfo;
981  Vector<Int> comps; 
982  ROArrayColumn<Double> parmsCol(t, "PARAMETERS");
983  ROArrayColumn<Bool> parmaskCol(t, "PARMASK");
984  ROArrayColumn<Int> compsCol(t, "COMPONENTS");
985  ROArrayColumn<String> funcsCol(t, "FUNCTIONS");
986  ROArrayColumn<String> finfoCol(t, "FRAMEINFO");
987  uInt n = t.nrow();
988  SDFitTable sdft;
989  for (uInt i=0; i<n; ++i) {
990    parmaskCol.get(i, parmask);
991    parmsCol.get(i, parms);
992    funcsCol.get(i, funcs);
993    compsCol.get(i, comps);
994    finfoCol.get(i, finfo);
995    sdft.addFit(parms, parmask, funcs, comps, finfo);
996  }
997  return sdft;
998}
999
1000SDFitTable SDMemTable::getSDFitTable(uInt whichRow) const {
1001  Array<Int> fitid;
1002  fitCol_.get(whichRow, fitid);
1003  if (fitid.nelements() == 0) return SDFitTable();
1004
1005  IPosition shp = fitid.shape();
1006  IPosition start(4, beamSel_, IFSel_, polSel_,0);
1007  IPosition end(4, beamSel_, IFSel_, polSel_, shp[3]-1);
1008
1009  // reform the output array slice to be of dim=1
1010  Vector<Int> tmp = (fitid(start, end)).reform(IPosition(1,shp[3]));
1011
1012  const Table& t = table_.keywordSet().asTable("FITS");
1013  Vector<Double> parms;
1014  Vector<Bool> parmask;
1015  Vector<String> funcs;
1016  Vector<String> finfo;
1017  Vector<Int> comps;   
1018  ROArrayColumn<Double> parmsCol(t, "PARAMETERS");
1019  ROArrayColumn<Bool> parmaskCol(t, "PARMASK");
1020  ROArrayColumn<Int> compsCol(t, "COMPONENTS");
1021  ROArrayColumn<String> funcsCol(t, "FUNCTIONS");
1022  ROArrayColumn<String> finfoCol(t, "FRAMEINFO");
1023  if (t.nrow() == 0) return SDFitTable();
1024  SDFitTable sdft;
1025  Int k=-1;
1026  for (uInt i=0; i< tmp.nelements(); ++i) {
1027    k = tmp[i];
1028    if ( k > -1 && k < t.nrow() ) {
1029      parms.resize();
1030      parmsCol.get(k, parms);
1031      parmask.resize();
1032      parmaskCol.get(k, parmask);
1033      funcs.resize();
1034      funcsCol.get(k, funcs);
1035      comps.resize();
1036      compsCol.get(k, comps);
1037      finfo.resize();
1038      finfoCol.get(k, finfo);
1039      sdft.addFit(parms, parmask, funcs, comps, finfo);
1040    }
1041  }
1042  return sdft;
1043}
1044
1045void SDMemTable::addFit(uInt whichRow,
1046                        const Vector<Double>& p, const Vector<Bool>& m,
1047                        const Vector<String>& f, const Vector<Int>& c)
1048{
1049  if (whichRow >= nRow()) {
1050    throw(AipsError("Specified row out of range"));
1051  }
1052  Table t = table_.keywordSet().asTable("FITS");
1053  uInt nrow = t.nrow(); 
1054  t.addRow();
1055  ArrayColumn<Double> parmsCol(t, "PARAMETERS");
1056  ArrayColumn<Bool> parmaskCol(t, "PARMASK");
1057  ArrayColumn<Int> compsCol(t, "COMPONENTS");
1058  ArrayColumn<String> funcsCol(t, "FUNCTIONS");
1059  ArrayColumn<String> finfoCol(t, "FRAMEINFO");
1060  parmsCol.put(nrow, p);
1061  parmaskCol.put(nrow, m);
1062  compsCol.put(nrow, c);
1063  funcsCol.put(nrow, f);
1064  Vector<String> fi = mathutil::toVectorString(getCoordInfo());
1065  finfoCol.put(nrow, fi);
1066
1067  Array<Int> fitarr;
1068  fitCol_.get(whichRow, fitarr);
1069
1070  Array<Int> newarr;               // The new Array containing the fitid
1071  Int pos =-1;                     // The fitid position in the array
1072  if ( fitarr.nelements() == 0 ) { // no fits at all in this row
1073    Array<Int> arr(IPosition(4,nBeam(),nIF(),nPol(),1));
1074    arr = -1;
1075    newarr.reference(arr);
1076    pos = 0;     
1077  } else {
1078    IPosition shp = fitarr.shape();
1079    IPosition start(4, beamSel_, IFSel_, polSel_,0);
1080    IPosition end(4, beamSel_, IFSel_, polSel_, shp[3]-1);
1081    // reform the output array slice to be of dim=1
1082    Array<Int> tmp = (fitarr(start, end)).reform(IPosition(1,shp[3]));
1083    const Vector<Int>& fits = tmp;
1084    VectorSTLIterator<Int> it(fits);
1085    Int i = 0;
1086    while (it != fits.end()) {
1087      if (*it == -1) {
1088        pos = i;
1089        break;
1090      }
1091      ++i;
1092      ++it;
1093    };
1094  }
1095  if (pos == -1) {
1096      mathutil::extendLastArrayAxis(newarr,fitarr, -1);
1097      pos = fitarr.shape()[3];     // the new element position
1098  } else {
1099    if (fitarr.nelements() > 0)
1100      newarr = fitarr;
1101  }
1102  newarr(IPosition(4, beamSel_, IFSel_, polSel_, pos)) = Int(nrow);
1103  fitCol_.put(whichRow, newarr);
1104
1105}
1106
1107SDFrequencyTable SDMemTable::getSDFreqTable() const
1108{
1109  const Table& t = table_.keywordSet().asTable("FREQUENCIES");
1110  SDFrequencyTable sdft;
1111
1112  // Add refpix/refval/incr.  What are the units ? Hz I suppose
1113  // but it's nowhere described...
1114  Vector<Double> refPix, refVal, incr;
1115  ScalarColumn<Double> refPixCol(t, "REFPIX");
1116  ScalarColumn<Double> refValCol(t, "REFVAL");
1117  ScalarColumn<Double> incrCol(t, "INCREMENT");
1118  refPix = refPixCol.getColumn();
1119  refVal = refValCol.getColumn();
1120  incr = incrCol.getColumn();
1121
1122  uInt n = refPix.nelements();
1123  for (uInt i=0; i<n; i++) {
1124     sdft.addFrequency(refPix[i], refVal[i], incr[i]);
1125  }
1126
1127  // Frequency reference frame.  I don't know if this
1128  // is the correct frame.  It might be 'REFFRAME'
1129  // rather than 'BASEREFFRAME' ?
1130  String baseFrame;
1131  t.keywordSet().get("BASEREFFRAME",baseFrame);
1132  sdft.setRefFrame(baseFrame);
1133
1134  // Equinox 
1135  Float equinox;
1136  t.keywordSet().get("EQUINOX", equinox);
1137  sdft.setEquinox(equinox);
1138
1139  // Rest Frequency
1140  Vector<Double> restFreqs;
1141  t.keywordSet().get("RESTFREQS", restFreqs);
1142  for (uInt i=0; i<restFreqs.nelements(); i++) {
1143     sdft.addRestFrequency(restFreqs[i]);
1144  }
1145  sdft.setRestFrequencyUnit(String("Hz"));
1146
1147  return sdft;
1148}
1149
1150bool SDMemTable::putSDContainer(const SDContainer& sdc)
1151{
1152  uInt rno = table_.nrow();
1153  table_.addRow();
1154
1155  timeCol_.put(rno, sdc.timestamp);
1156  srcnCol_.put(rno, sdc.sourcename);
1157  fldnCol_.put(rno, sdc.fieldname);
1158  specCol_.put(rno, sdc.getSpectrum());
1159  flagsCol_.put(rno, sdc.getFlags());
1160  tsCol_.put(rno, sdc.getTsys());
1161  scanCol_.put(rno, sdc.scanid);
1162  integrCol_.put(rno, sdc.interval);
1163  freqidCol_.put(rno, sdc.getFreqMap());
1164  restfreqidCol_.put(rno, sdc.getRestFreqMap());
1165  dirCol_.put(rno, sdc.getDirection());
1166  rbeamCol_.put(rno, sdc.refbeam);
1167  tcalCol_.put(rno, sdc.tcal);
1168  tcaltCol_.put(rno, sdc.tcaltime);
1169  azCol_.put(rno, sdc.azimuth);
1170  elCol_.put(rno, sdc.elevation);
1171  paraCol_.put(rno, sdc.parangle);
1172  histCol_.put(rno, sdc.getHistory());
1173  fitCol_.put(rno, sdc.getFitMap());
1174
1175  return true;
1176}
1177
1178SDContainer SDMemTable::getSDContainer(uInt whichRow) const
1179{
1180  SDContainer sdc(nBeam(),nIF(),nPol(),nChan());
1181  timeCol_.get(whichRow, sdc.timestamp);
1182  srcnCol_.get(whichRow, sdc.sourcename);
1183  integrCol_.get(whichRow, sdc.interval);
1184  scanCol_.get(whichRow, sdc.scanid);
1185  fldnCol_.get(whichRow, sdc.fieldname);
1186  rbeamCol_.get(whichRow, sdc.refbeam);
1187  azCol_.get(whichRow, sdc.azimuth);
1188  elCol_.get(whichRow, sdc.elevation);
1189  paraCol_.get(whichRow, sdc.parangle);
1190  Vector<Float> tc;
1191  tcalCol_.get(whichRow, tc);
1192  sdc.tcal[0] = tc[0];sdc.tcal[1] = tc[1];
1193  tcaltCol_.get(whichRow, sdc.tcaltime);
1194
1195  Array<Float> spectrum;
1196  Array<Float> tsys;
1197  Array<uChar> flagtrum;
1198  Vector<uInt> fmap;
1199  Array<Double> direction;
1200  Vector<String> histo;
1201  Array<Int> fits;
1202 
1203  specCol_.get(whichRow, spectrum);
1204  sdc.putSpectrum(spectrum);
1205  flagsCol_.get(whichRow, flagtrum);
1206  sdc.putFlags(flagtrum);
1207  tsCol_.get(whichRow, tsys);
1208  sdc.putTsys(tsys);
1209  freqidCol_.get(whichRow, fmap);
1210  sdc.putFreqMap(fmap);
1211  restfreqidCol_.get(whichRow, fmap);
1212  sdc.putRestFreqMap(fmap);
1213  dirCol_.get(whichRow, direction);
1214  sdc.putDirection(direction);
1215  histCol_.get(whichRow, histo);
1216  sdc.putHistory(histo);
1217  fitCol_.get(whichRow, fits);
1218  sdc.putFitMap(fits);
1219  return sdc;
1220}
1221
1222bool SDMemTable::putSDHeader(const SDHeader& sdh)
1223{
1224  table_.rwKeywordSet().define("nIF", sdh.nif);
1225  table_.rwKeywordSet().define("nBeam", sdh.nbeam);
1226  table_.rwKeywordSet().define("nPol", sdh.npol);
1227  table_.rwKeywordSet().define("nChan", sdh.nchan);
1228  table_.rwKeywordSet().define("Observer", sdh.observer);
1229  table_.rwKeywordSet().define("Project", sdh.project);
1230  table_.rwKeywordSet().define("Obstype", sdh.obstype);
1231  table_.rwKeywordSet().define("AntennaName", sdh.antennaname);
1232  table_.rwKeywordSet().define("AntennaPosition", sdh.antennaposition);
1233  table_.rwKeywordSet().define("Equinox", sdh.equinox);
1234  table_.rwKeywordSet().define("FreqRefFrame", sdh.freqref);
1235  table_.rwKeywordSet().define("FreqRefVal", sdh.reffreq);
1236  table_.rwKeywordSet().define("Bandwidth", sdh.bandwidth);
1237  table_.rwKeywordSet().define("UTC", sdh.utc);
1238  table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit);
1239  table_.rwKeywordSet().define("Epoch", sdh.epoch);
1240  return true;
1241}
1242
1243SDHeader SDMemTable::getSDHeader() const
1244{
1245  SDHeader sdh;
1246  table_.keywordSet().get("nBeam",sdh.nbeam);
1247  table_.keywordSet().get("nIF",sdh.nif);
1248  table_.keywordSet().get("nPol",sdh.npol);
1249  table_.keywordSet().get("nChan",sdh.nchan);
1250  table_.keywordSet().get("Observer", sdh.observer);
1251  table_.keywordSet().get("Project", sdh.project);
1252  table_.keywordSet().get("Obstype", sdh.obstype);
1253  table_.keywordSet().get("AntennaName", sdh.antennaname);
1254  table_.keywordSet().get("AntennaPosition", sdh.antennaposition);
1255  table_.keywordSet().get("Equinox", sdh.equinox);
1256  table_.keywordSet().get("FreqRefFrame", sdh.freqref);
1257  table_.keywordSet().get("FreqRefVal", sdh.reffreq);
1258  table_.keywordSet().get("Bandwidth", sdh.bandwidth);
1259  table_.keywordSet().get("UTC", sdh.utc);
1260  table_.keywordSet().get("FluxUnit", sdh.fluxunit);
1261  table_.keywordSet().get("Epoch", sdh.epoch);
1262  return sdh;
1263}
1264void SDMemTable::makePersistent(const std::string& filename)
1265{
1266  table_.deepCopy(filename,Table::New);
1267}
1268
1269Int SDMemTable::nScan() const {
1270  Int n = 0;
1271  Int previous = -1;Int current=0;
1272  for (uInt i=0; i< scanCol_.nrow();i++) {
1273    scanCol_.getScalar(i,current);
1274    if (previous != current) {
1275      previous = current;
1276      n++;
1277    }
1278  }
1279  return n;
1280}
1281
1282String SDMemTable::formatSec(Double x) const
1283{
1284  Double xcop = x;
1285  MVTime mvt(xcop/24./3600.);  // make days
1286
1287  if (x < 59.95)
1288    return  String("      ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
1289  else if (x < 3599.95)
1290    return String("   ") + mvt.string(MVTime::TIME_CLEAN_NO_H,7)+" ";
1291  else {
1292    ostringstream oss;
1293    oss << setw(2) << std::right << setprecision(1) << mvt.hour();
1294    oss << ":" << mvt.string(MVTime::TIME_CLEAN_NO_H,7) << " ";
1295    return String(oss);
1296  }   
1297};
1298
1299String SDMemTable::formatDirection(const MDirection& md) const
1300{
1301  Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue();
1302  Int prec = 7;
1303
1304  MVAngle mvLon(t[0]);
1305  String sLon = mvLon.string(MVAngle::TIME,prec);
1306  MVAngle mvLat(t[1]);
1307  String sLat = mvLat.string(MVAngle::ANGLE+MVAngle::DIG2,prec);
1308  return sLon + String(" ") + sLat;
1309}
1310
1311
1312std::string SDMemTable::getFluxUnit() const
1313{
1314  String tmp;
1315  table_.keywordSet().get("FluxUnit", tmp);
1316  return tmp;
1317}
1318
1319void SDMemTable::setFluxUnit(const std::string& unit)
1320{
1321  String tmp(unit);
1322  Unit tU(tmp);
1323  if (tU==Unit("K") || tU==Unit("Jy")) {
1324     table_.rwKeywordSet().define(String("FluxUnit"), tmp);
1325  } else {
1326     throw AipsError("Illegal unit - must be compatible with Jy or K");
1327  }
1328}
1329
1330
1331void SDMemTable::setInstrument(const std::string& name)
1332{
1333  Bool throwIt = True;
1334  Instrument ins = convertInstrument(name, throwIt);
1335  String nameU(name);
1336  nameU.upcase();
1337  table_.rwKeywordSet().define(String("AntennaName"), nameU);
1338}
1339
1340std::string SDMemTable::summary(bool verbose) const  {
1341
1342  // Format header info
1343  ostringstream oss;
1344  oss << endl;
1345  oss << "--------------------------------------------------------------------------------" << endl;
1346  oss << " Scan Table Summary" << endl;
1347  oss << "--------------------------------------------------------------------------------" << endl;
1348  oss.flags(std::ios_base::left);
1349  oss << setw(15) << "Beams:" << setw(4) << nBeam() << endl
1350      << setw(15) << "IFs:" << setw(4) << nIF() << endl
1351      << setw(15) << "Polarisations:" << setw(4) << nPol() << endl
1352      << setw(15) << "Channels:"  << setw(4) << nChan() << endl;
1353  oss << endl;
1354  String tmp;
1355  table_.keywordSet().get("Observer", tmp);
1356  oss << setw(15) << "Observer:" << tmp << endl;
1357  oss << setw(15) << "Obs Date:" << getTime(-1,True) << endl;
1358  table_.keywordSet().get("Project", tmp);
1359  oss << setw(15) << "Project:" << tmp << endl;
1360  table_.keywordSet().get("Obstype", tmp);
1361  oss << setw(15) << "Obs. Type:" << tmp << endl;
1362  table_.keywordSet().get("AntennaName", tmp);
1363  oss << setw(15) << "Antenna Name:" << tmp << endl;
1364  table_.keywordSet().get("FluxUnit", tmp);
1365  oss << setw(15) << "Flux Unit:" << tmp << endl;
1366  Table t = table_.keywordSet().asTable("FREQUENCIES");
1367  Vector<Double> vec;
1368  t.keywordSet().get("RESTFREQS",vec);
1369  oss << setw(15) << "Rest Freqs:";
1370  if (vec.nelements() > 0) {
1371      oss << setprecision(0) << vec << " [Hz]" << endl;
1372  } else {
1373      oss << "None set" << endl;
1374  }
1375  oss << setw(15) << "Abcissa:" << getAbcissaString() << endl;
1376  oss << setw(15) << "Cursor:" << "Beam[" << getBeam() << "] "
1377      << "IF[" << getIF() << "] " << "Pol[" << getPol() << "]" << endl;
1378  oss << endl;
1379
1380  String dirtype ="Position ("+ MDirection::showType(getDirectionReference()) + ")";
1381  oss << setw(5) << "Scan"
1382      << setw(15) << "Source"
1383      << setw(24) << dirtype
1384      << setw(10) << "Time"
1385      << setw(18) << "Integration"
1386      << setw(7) << "FreqIDs" << endl;
1387  oss << "--------------------------------------------------------------------------------" << endl;
1388 
1389  // Generate list of scan start and end integrations
1390  Vector<Int> scanIDs = scanCol_.getColumn();
1391  Vector<uInt> startInt, endInt;
1392  mathutil::scanBoundaries(startInt, endInt, scanIDs);
1393
1394  const uInt nScans = startInt.nelements();
1395  String name;
1396  Vector<uInt> freqIDs, listFQ;
1397  uInt nInt;
1398
1399  for (uInt i=0; i<nScans; i++) {
1400    // Get things from first integration of scan
1401    String time = getTime(startInt(i),False);
1402    String tInt = formatSec(Double(getInterval(startInt(i))));
1403    String posit = formatDirection(getDirection(startInt(i),True));
1404    srcnCol_.getScalar(startInt(i),name);
1405
1406    // Find all the FreqIDs in this scan
1407    listFQ.resize(0);     
1408    for (uInt j=startInt(i); j<endInt(i)+1; j++) {
1409      freqidCol_.get(j, freqIDs);
1410      for (uInt k=0; k<freqIDs.nelements(); k++) {
1411        mathutil::addEntry(listFQ, freqIDs(k));
1412      }
1413    }
1414
1415    nInt = endInt(i) - startInt(i) + 1;
1416    oss << setw(3) << std::right << i << std::left << setw(2) << "  "
1417        << setw(15) << name
1418        << setw(24) << posit
1419        << setw(10) << time
1420        << setw(3) << std::right << nInt  << setw(3) << " x " << std::left
1421        << setw(6) <<  tInt
1422        << " " << listFQ << endl;
1423  }
1424  oss << endl;
1425  oss << "Table contains " << table_.nrow() << " integration(s) in "
1426      << nScans << " scan(s)." << endl;
1427
1428  // Frequency Table
1429  if (verbose) {
1430    std::vector<string> info = getCoordInfo();
1431    SDFrequencyTable sdft = getSDFreqTable();
1432    oss << endl << endl;
1433    oss << "FreqID  Frame   RefFreq(Hz)     RefPix   Increment(Hz)" << endl;
1434    oss << "--------------------------------------------------------------------------------" << endl;
1435    for (uInt i=0; i<sdft.length(); i++) {
1436      oss << setw(8) << i << setw(8)
1437          << info[3] << setw(16) << setprecision(8)
1438          << sdft.referenceValue(i) << setw(10)
1439          << sdft.referencePixel(i) << setw(12)
1440          << sdft.increment(i) << endl;
1441    }
1442    oss << "--------------------------------------------------------------------------------" << endl;
1443  }
1444  return String(oss);
1445}
1446
1447Int SDMemTable::nBeam() const
1448{
1449  Int n;
1450  table_.keywordSet().get("nBeam",n);
1451  return n;
1452}
1453
1454Int SDMemTable::nIF() const {
1455  Int n;
1456  table_.keywordSet().get("nIF",n);
1457  return n;
1458}
1459
1460Int SDMemTable::nPol() const {
1461  Int n;
1462  table_.keywordSet().get("nPol",n);
1463  return n;
1464}
1465
1466Int SDMemTable::nChan() const {
1467  Int n;
1468  table_.keywordSet().get("nChan",n);
1469  return n;
1470}
1471
1472bool SDMemTable::appendHistory(const std::string& hist, int whichRow)
1473{
1474  Vector<String> history;
1475  histCol_.get(whichRow, history);
1476  history.resize(history.nelements()+1,True);
1477  history[history.nelements()-1] = hist;
1478  histCol_.put(whichRow, history);
1479}
1480
1481std::vector<std::string> SDMemTable::history(int whichRow) const
1482{
1483  Vector<String> history;
1484  histCol_.get(whichRow, history);
1485  std::vector<std::string> stlout = mathutil::tovectorstring(history);
1486  return stlout;
1487}
1488/*
1489void SDMemTable::maskChannels(const std::vector<Int>& whichChans ) {
1490
1491  std::vector<int>::iterator it;
1492  ArrayAccessor<uChar, Axis<asap::PolAxis> > j(flags_);
1493  for (it = whichChans.begin(); it != whichChans.end(); it++) {
1494    j.reset(j.begin(uInt(*it)));
1495    for (ArrayAccessor<uChar, Axis<asap::BeamAxis> > i(j); i != i.end(); ++i) {
1496      for (ArrayAccessor<uChar, Axis<asap::IFAxis> > ii(i); ii != ii.end(); ++ii) {
1497        for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > iii(ii);
1498             iii != iii.end(); ++iii) {
1499          (*iii) =
1500        }
1501      }
1502    }
1503  }
1504
1505}
1506*/
1507void SDMemTable::flag(int whichRow)
1508{
1509  Array<uChar> arr;
1510  flagsCol_.get(whichRow, arr);
1511
1512  ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
1513  aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
1514  ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
1515  aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
1516  ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
1517  aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
1518
1519  for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
1520    (*i) = uChar(True);
1521  }
1522
1523  flagsCol_.put(whichRow, arr);
1524}
1525
1526MDirection::Types SDMemTable::getDirectionReference() const
1527{ 
1528  Float eq;
1529  table_.keywordSet().get("Equinox",eq);
1530  std::map<float,string> mp;
1531  mp[2000.0] = "J2000";
1532  mp[1950.0] = "B1950";
1533  MDirection::Types mdr;
1534  if (!MDirection::getType(mdr, mp[eq])) {   
1535    mdr = MDirection::J2000;
1536    cerr  << "Unknown equinox using J2000" << endl;
1537  }
1538
1539  return mdr;
1540}
1541
1542MEpoch::Types SDMemTable::getTimeReference() const
1543{
1544  MEpoch::Types met;
1545  String ep;
1546  table_.keywordSet().get("Epoch",ep);
1547  if (!MEpoch::getType(met, ep)) {
1548    cerr << "Epoch type unknown - using UTC" << endl;
1549    met = MEpoch::UTC;
1550  }
1551
1552  return met;
1553}
1554
1555
1556Instrument SDMemTable::convertInstrument(const String& instrument,
1557                                         Bool throwIt)
1558{
1559   String t(instrument);
1560   t.upcase();
1561
1562   // The strings are what SDReader returns, after cunning interrogation
1563   // of station names... :-(
1564   Instrument inst = asap::UNKNOWN;
1565   if (t==String("DSS-43")) {               
1566      inst = TIDBINBILLA;
1567   } else if (t==String("ATPKSMB")) {
1568      inst = ATPKSMB;
1569   } else if (t==String("ATPKSHOH")) {
1570      inst = ATPKSHOH;
1571   } else if (t==String("ATMOPRA")) {
1572      inst = ATMOPRA;
1573   } else if (t==String("CEDUNA")) {
1574      inst = CEDUNA;
1575   } else if (t==String("HOBART")) {
1576      inst = HOBART;
1577   } else {
1578     if (throwIt) {
1579       throw AipsError("Unrecognized instrument - use function scan.set_instrument to set");
1580     }
1581   }
1582   return inst;
1583}
1584
1585Bool SDMemTable::setRestFreqs(const Vector<Double>& restFreqsIn,
1586                              const String& sUnit,
1587                              const vector<string>& lines,
1588                              const String& source,
1589                              Int whichIF)
1590{
1591   const Int nIFs = nIF();
1592   if (whichIF>=nIFs) {
1593      throw(AipsError("Illegal IF index"));
1594   }
1595
1596   // Find vector of restfrequencies
1597   // Double takes precedence over String
1598   Unit unit;
1599   Vector<Double> restFreqs;
1600   if (restFreqsIn.nelements()>0) {
1601      restFreqs.resize(restFreqsIn.nelements());
1602      restFreqs = restFreqsIn;
1603      unit = Unit(sUnit);
1604   } else if (lines.size()>0) {
1605      const uInt nLines = lines.size();
1606      unit = Unit("Hz");
1607      restFreqs.resize(nLines);
1608      MFrequency lineFreq;
1609      for (uInt i=0; i<nLines; i++) {
1610         String tS(lines[i]);
1611         tS.upcase();
1612         if (MeasTable::Line(lineFreq, tS)) {
1613            restFreqs[i] = lineFreq.getValue().getValue();          // Hz
1614         } else {
1615            String s = String(lines[i]) +
1616              String(" is an unrecognized spectral line");
1617            throw(AipsError(s));
1618         }
1619      }
1620   } else {
1621      throw(AipsError("You have not specified any rest frequencies or lines"));
1622   }
1623
1624   // If multiple restfreqs, must be length nIF. In this
1625   // case we will just replace the rest frequencies
1626   const uInt nRestFreqs = restFreqs.nelements();
1627   Int idx = -1;
1628   SDFrequencyTable sdft = getSDFreqTable();
1629
1630   if (nRestFreqs>1) {
1631     // Replace restFreqs, one per IF
1632      if (nRestFreqs != nIFs) {
1633         throw (AipsError("Number of rest frequencies must be equal to the number of IFs"));
1634      }
1635      cout << "Replacing rest frequencies with given list, one per IF" << endl;
1636      sdft.deleteRestFrequencies();
1637      for (uInt i=0; i<nRestFreqs; i++) {
1638         Quantum<Double> rf(restFreqs[i], unit);
1639         sdft.addRestFrequency(rf.getValue("Hz"));
1640      }
1641   } else {
1642
1643     // Add new rest freq
1644      Quantum<Double> rf(restFreqs[0], unit);
1645      idx = sdft.addRestFrequency(rf.getValue("Hz"));
1646      cout << "Selecting given rest frequency" << endl;
1647   }
1648   
1649   // Replace
1650   Bool empty = source.empty();
1651   Bool ok = False;
1652   if (putSDFreqTable(sdft)) {
1653      const uInt nRow = table_.nrow();
1654      String srcName;
1655      Vector<uInt> restFreqIDs;
1656      for (uInt i=0; i<nRow; i++) {
1657         srcnCol_.get(i, srcName);
1658         restfreqidCol_.get(i,restFreqIDs);       
1659
1660         if (idx==-1) {
1661           // Replace vector of restFreqs; one per IF.
1662           // No selection possible
1663            for (uInt i=0; i<nIFs; i++) restFreqIDs[i] = i;
1664         } else {
1665           // Set RestFreqID for selected data
1666            if (empty || source==srcName) {
1667               if (whichIF<0) {
1668                  restFreqIDs = idx;
1669               } else {             
1670                  restFreqIDs[whichIF] = idx;
1671               }
1672            }
1673         }
1674         restfreqidCol_.put(i,restFreqIDs);       
1675      }
1676      ok = True;
1677   } else {
1678     ok = False;
1679   }
1680
1681   return ok;
1682}
1683
1684void SDMemTable::spectralLines() const
1685{
1686   Vector<String> lines = MeasTable::Lines();
1687   MFrequency lineFreq;
1688   Double freq;
1689
1690   cout.flags(std::ios_base::left);
1691   cout << "Line      Frequency (Hz)" << endl;
1692   cout << "-----------------------" << endl;
1693   for (uInt i=0; i<lines.nelements(); i++) {
1694     MeasTable::Line(lineFreq, lines[i]);
1695     freq = lineFreq.getValue().getValue();          // Hz
1696     cout << setw(11) << lines[i] << setprecision(10) << freq << endl;
1697   }
1698}
1699
1700void SDMemTable::renumber()
1701{
1702  uInt nRow = scanCol_.nrow();
1703  Int newscanid = 0;
1704  Int cIdx;// the current scanid
1705  // get the first scanid
1706  scanCol_.getScalar(0,cIdx);
1707  Int pIdx = cIdx;// the scanid of the previous row
1708  for (uInt i=0; i<nRow;++i) {
1709    scanCol_.getScalar(i,cIdx);
1710    if (pIdx == cIdx) {
1711      // renumber
1712      scanCol_.put(i,newscanid);
1713    } else {
1714      ++newscanid;
1715      pIdx = cIdx;   // store scanid
1716      --i;           // don't increment next loop
1717    }
1718  }
1719}
1720
1721
1722void SDMemTable::setCursorSlice(IPosition& start, IPosition& end,
1723                                const IPosition& shape) const
1724{
1725  const uInt nDim = shape.nelements();
1726  start.resize(nDim);
1727  end.resize(nDim);
1728 
1729  start(asap::BeamAxis) = beamSel_;
1730  end(asap::BeamAxis) = beamSel_;
1731  start(asap::IFAxis) = IFSel_;
1732  end(asap::IFAxis) = IFSel_;
1733
1734  start(asap::PolAxis) = polSel_;
1735  end(asap::PolAxis) = polSel_;
1736
1737  start(asap::ChanAxis) = 0;
1738  end(asap::ChanAxis) = shape(asap::ChanAxis) - 1;
1739}
1740
1741
1742std::vector<float> SDMemTable::getFloatSpectrum(const Array<Float>& arr) const
1743  // Get spectrum at cursor location
1744{
1745
1746  // Setup accessors
1747  ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
1748  aa0.reset(aa0.begin(uInt(beamSel_)));                    // Beam selection
1749
1750  ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
1751  aa1.reset(aa1.begin(uInt(IFSel_)));                      // IF selection
1752
1753  ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
1754  aa2.reset(aa2.begin(uInt(polSel_)));                     // Pol selection
1755 
1756  std::vector<float> spectrum;
1757  for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
1758    spectrum.push_back(*i);
1759  }
1760  return spectrum;
1761}
1762
Note: See TracBrowser for help on using the repository browser.