Changeset 455 for trunk/src

Timestamp:

02/16/05 12:29:26 (20 years ago)

Author:

mar637

Message:

• added FITS subtable for fitting and relevant functions to set fits and put/get SDFitContainer
• cleanup

Location:

trunk/src

Files:

• SDMemTable.cc (modified) (70 diffs)
• SDMemTable.h (modified) (9 diffs)
• SDMemTableWrapper.h (modified) (2 diffs)

trunk/src/SDMemTable.cc

@@ -40 +40 @@
 #include <casa/Arrays/ArrayLogical.h>
 #include <casa/Arrays/ArrayAccessor.h>
+#include <casa/Arrays/VectorSTLIterator.h>
 #include <casa/Arrays/Vector.h>
 #include <casa/BasicMath/Math.h>
@@ -113 +114 @@
   polSel_(0)
 {
-  cout << exprs << endl;
   Table t = tableCommand(exprs,tab);
   if (t.nrow() == 0)
@@ -160 +160 @@
   TableDesc td("", "1", TableDesc::Scratch);
   td.comment() = "A SDMemTable";
-//
+  
   td.addColumn(ScalarColumnDesc<Double>("TIME"));
   td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
@@ -180 +180 @@
   td.addColumn(ScalarColumnDesc<Int>("REFBEAM"));
   td.addColumn(ArrayColumnDesc<String>("HISTORY"));
+  td.addColumn(ArrayColumnDesc<Int>("FITID"));
+
 
   // Now create Table SetUp from the description.
@@ -185 +187 @@
   SetupNewTable aNewTab("dummy", td, Table::New);
 
-// Bind the Stokes Virtual machine to the STOKES column
-// Because we don't know how many polarizations will
-// be in the data at this point, we must bind the
-// Virtual Engine regardless.  The STOKES column won't
-// be accessed if not appropriate (nPol=4)
+  // Bind the Stokes Virtual machine to the STOKES column Because we
+  // don't know how many polarizations will be in the data at this
+  // point, we must bind the Virtual Engine regardless.  The STOKES
+  // column won't be accessed if not appropriate (nPol=4)
 
    SDStokesEngine stokesEngine(String("STOKES"), String("SPECTRA"));
    aNewTab.bindColumn ("STOKES", stokesEngine);
 
-// Create Table
-
+   // Create Table
   table_ = Table(aNewTab, Table::Memory, 0);
-}
-
+  // add subtable 
+  TableDesc tdf("", "1", TableDesc::Scratch);
+  tdf.addColumn(ArrayColumnDesc<String>("FUNCTIONS"));
+  tdf.addColumn(ArrayColumnDesc<Int>("COMPONENTS"));
+  tdf.addColumn(ArrayColumnDesc<Double>("PARAMETERS"));
+  tdf.addColumn(ArrayColumnDesc<Bool>("PARMASK"));
+  SetupNewTable fittab("fits", tdf, Table::New);
+  Table fitTable(fittab, Table::Memory);
+  table_.rwKeywordSet().defineTable("FITS", fitTable);
+
+
+}
 
 void SDMemTable::attach()
@@ -221 +231 @@
   rbeamCol_.attach(table_, "REFBEAM");
   histCol_.attach(table_, "HISTORY");
+  fitCol_.attach(table_,"FITID");
 }
 
@@ -311 +322 @@
 
   std::vector<bool> mask;
-//
+
   Array<uChar> arr;
   flagsCol_.get(whichRow, arr);
-//
+
   ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
   aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
@@ -347 +358 @@
   Array<Float> arr;
   specCol_.get(whichRow, arr);
-  return getFloatSpectrum (arr);
-}
-
-std::vector<float> SDMemTable::getStokesSpectrum(Int whichRow, Bool doPol, Float paOffset) const 
-//
-// Gets
-//  doPol=False  : I,Q,U,V
-//  doPol=True   : I,P,PA,V   ; P = sqrt(Q**2+U**2), PA = 0.5*atan2(Q,U)
-//
+  return getFloatSpectrum(arr);
+}
+
+std::vector<float> SDMemTable::getStokesSpectrum(Int whichRow, Bool doPol, 
+                                                 Float paOffset) const 
+  // Gets
+  //  doPol=False  : I,Q,U,V
+  //  doPol=True   : I,P,PA,V   ; P = sqrt(Q**2+U**2), PA = 0.5*atan2(Q,U)
+  //
 {
   AlwaysAssert(asap::nAxes==4,AipsError);
@@ -363 +374 @@
   Array<Float> arr;
   stokesCol_.get(whichRow, arr);
-//
+
   if (doPol && (polSel_==1 || polSel_==2)) {       //   Q,U --> P, P.A.
 
-// Set current cursor location
-
+    // Set current cursor location
      const IPosition& shape = arr.shape();
      IPosition start, end;
      setCursorSlice (start, end, shape);
 
-// Get Q and U slices
-
-     Array<Float> Q = SDPolUtil::getStokesSlice (arr,start,end,"Q");
-     Array<Float> U = SDPolUtil::getStokesSlice (arr,start,end,"U");
-
-// Compute output 
+     // Get Q and U slices
+
+     Array<Float> Q = SDPolUtil::getStokesSlice(arr,start,end,"Q");
+     Array<Float> U = SDPolUtil::getStokesSlice(arr,start,end,"U");
+
+     // Compute output 
 
      Array<Float> out;
@@ -386 +396 @@
      }
 
-// Copy to output
+     // Copy to output
 
      IPosition vecShape(1,shape(asap::ChanAxis));
      Vector<Float> outV = out.reform(vecShape);
-     return convertVector(outV);
+     std::vector<float> stlout;
+     outV.tovector(stlout);
+     return stlout;
+
   } else {
 
-// Selects at the cursor location
-
-    return getFloatSpectrum (arr);
-  }
-}
-
-std::vector<float> SDMemTable::getCircularSpectrum(Int whichRow, Bool doRR) const 
-//
-// Gets
-//  RR = I + V
-//  LL = I - V
-//
+    // Selects at the cursor location
+    return getFloatSpectrum(arr);
+  }
+}
+
+std::vector<float> SDMemTable::getCircularSpectrum(Int whichRow, 
+                                                   Bool doRR) const 
+  // Gets
+  //  RR = I + V
+  //  LL = I - V
 {
   AlwaysAssert(asap::nAxes==4,AipsError);
@@ -413 +424 @@
   stokesCol_.get(whichRow, arr);
 
-// Set current cursor location
+  // Set current cursor location
 
   const IPosition& shape = arr.shape();
   IPosition start, end;
-  setCursorSlice (start, end, shape);
-
-// Get I and V slices
-
-  Array<Float> I = SDPolUtil::getStokesSlice (arr,start,end,"I");
-  Array<Float> V = SDPolUtil::getStokesSlice (arr,start,end,"V");
-
-// Compute output 
-
-  Array<Float> out = SDPolUtil::circularPolarizationFromStokes (I, V, doRR);
-
-// Copy to output
+  setCursorSlice(start, end, shape);
+
+  // Get I and V slices
+
+  Array<Float> I = SDPolUtil::getStokesSlice(arr,start,end,"I");
+  Array<Float> V = SDPolUtil::getStokesSlice(arr,start,end,"V");
+
+  // Compute output 
+
+  Array<Float> out = SDPolUtil::circularPolarizationFromStokes(I, V, doRR);
+
+  // Copy to output
 
   IPosition vecShape(1,shape(asap::ChanAxis));
   Vector<Float> outV = out.reform(vecShape);
-  return convertVector(outV);
+  std::vector<float> stlout;
+  outV.tovector(stlout);
+
+  return stlout;
 }
 
@@ -464 +478 @@
   dpl = theinfo[2];             // Doppler
   brfrm = theinfo[3];           // Base frame
-//
+
   Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
-//
+
   Vector<Double> rstf;
   t.keywordSet().get("RESTFREQS",rstf);
-//
+
   Bool canDo = True;
   Unit u1("km/s");Unit u2("Hz");
@@ -482 +496 @@
   }
   t.rwKeywordSet().define("UNIT", un);
-//
+
   MFrequency::Types mdr;
   if (!MFrequency::getType(mdr, rfrm)) {
@@ -493 +507 @@
     t.rwKeywordSet().define("REFFRAME",rfrm);
   }
-//
+
   MDoppler::Types dtype;
   dpl.upcase();
@@ -501 +515 @@
     t.rwKeywordSet().define("DOPPLER",dpl);
   }
-//
+
   if (!MFrequency::getType(mdr, brfrm)) {
      Int a,b;const uInt* c;
@@ -517 +531 @@
   std::vector<double> abc(nChan());
 
-// Get header units keyword
-
+  // Get header units keyword
   Table t = table_.keywordSet().asTable("FREQUENCIES");
   String sunit;
@@ -525 +538 @@
   Unit u(sunit);
 
-// Easy if just wanting pixels
-
+  // Easy if just wanting pixels
   if (sunit==String("pixel")) {
     // assume channels/pixels
@@ -534 +546 @@
       (*it) = Double(i++);
     }
-//
     return abc;
   }
 
-// Continue with km/s or Hz.  Get FreqIDs
-
+  // Continue with km/s or Hz.  Get FreqIDs
   Vector<uInt> freqIDs;
   freqidCol_.get(whichRow, freqIDs);
@@ -546 +556 @@
   uInt restFreqID = freqIDs(IFSel_);
 
-// Get SpectralCoordinate, set reference frame conversion,
-// velocity conversion, and rest freq state
+  // Get SpectralCoordinate, set reference frame conversion,
+  // velocity conversion, and rest freq state
 
   SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow);
   Vector<Double> pixel(nChan());
   indgen(pixel);
-//
+
   if (u == Unit("km/s")) {
      Vector<Double> world;
@@ -564 +574 @@
   } else if (u == Unit("Hz")) {
 
-// Set world axis units
-
+    // Set world axis units    
     Vector<String> wau(1); wau = u.getName();
     spc.setWorldAxisUnits(wau);
-//
+
     std::vector<double>::iterator it;
     Double tmp;
@@ -584 +593 @@
 {
   Table t = table_.keywordSet().asTable("FREQUENCIES");
-//
+
   String sunit;
   t.keywordSet().get("UNIT",sunit);
   if (sunit == "") sunit = "pixel";
   Unit u(sunit);
-//
+
   Vector<uInt> freqIDs;
   freqidCol_.get(whichRow, freqIDs);
@@ -596 +605 @@
   uInt restFreqID = freqIDs(IFSel_);
 
-// Get SpectralCoordinate, with frame, velocity, rest freq state set
-
+  // Get SpectralCoordinate, with frame, velocity, rest freq state set
   SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow);
-//
+
   String s = "Channel";
   if (u == Unit("km/s")) { 
@@ -606 +614 @@
     Vector<String> wau(1);wau = u.getName();
     spc.setWorldAxisUnits(wau);
-//
+
     s = CoordinateUtil::axisLabel(spc,0,True,True,False);
   }
@@ -616 +624 @@
   Array<Float> arr;
   specCol_.get(whichRow, arr);
-  if (spectrum.size() != arr.shape()(3)) {
+  if (spectrum.size() != arr.shape()(asap::ChanAxis)) {
     throw(AipsError("Attempting to set spectrum with incorrect length."));
   }
 
-// Setup accessors
-
+  // Setup accessors
   ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
   aa0.reset(aa0.begin(uInt(beamSel_)));                   // Beam selection
@@ -629 +636 @@
   aa2.reset(aa2.begin(uInt(polSel_)));                    // Pol selection
 
-// Iterate
-
+  // Iterate
   std::vector<float>::iterator it = spectrum.begin();
   for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
@@ -644 +650 @@
   specCol_.get(whichRow, arr);
 
-// Iterate and extract
-
+  // Iterate and extract
   spectrum.resize(arr.shape()(3));
   ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
@@ -653 +658 @@
   ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
   aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
-//
+
   ArrayAccessor<Float, Axis<asap::BeamAxis> > va(spectrum);
   for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
@@ -696 +701 @@
 */
 
-MaskedArray<Float> SDMemTable::rowAsMaskedArray(uInt whichRow, Bool toStokes) const 
-{
-
-// Get flags
-
+MaskedArray<Float> SDMemTable::rowAsMaskedArray(uInt whichRow, 
+                                                Bool toStokes) const 
+{
+  // Get flags
   Array<uChar> farr;
   flagsCol_.get(whichRow, farr);
 
-// Get data and convert mask to Bool
-
+  // Get data and convert mask to Bool
   Array<Float> arr;
   Array<Bool> mask;
   if (toStokes) {
      stokesCol_.get(whichRow, arr);
-//
+
      Array<Bool> tMask(farr.shape());
      convertArray(tMask, farr);
@@ -719 +722 @@
      convertArray(mask, farr);
   }
-//
+
   return MaskedArray<Float>(arr,!mask);
 }
@@ -728 +731 @@
   tsCol_.get(whichRow, arr);
   Float out;
-//
+
   IPosition ip(arr.shape());
   ip(asap::BeamAxis) = beamSel_;
@@ -734 +737 @@
   ip(asap::PolAxis) = polSel_;
   ip(asap::ChanAxis)=0;               // First channel only
-//
+
   out = arr(ip);
   return out;
@@ -762 +765 @@
 {
   MEpoch::Types met = getTimeReference();
-//
+
   Double obstime;
   timeCol_.get(whichRow,obstime);
@@ -793 +796 @@
   t.keywordSet().get("BASEREFFRAME",rf);
 
-// Create SpectralCoordinate (units Hz)
-
+  // Create SpectralCoordinate (units Hz)
   MFrequency::Types mft;
   if (!MFrequency::getType(mft, rf)) {
@@ -803 +805 @@
   rvc.get(freqID, rv);
   incc.get(freqID, inc);
-//
+
   SpectralCoordinate spec(mft,rv,inc,rp);
-//
   return spec;
 }
 
 
-SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID, uInt restFreqID,
+SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID, 
+                                                     uInt restFreqID,
                                                      uInt whichRow) const
 {
-
-// Create basic SC
-
+  
+  // Create basic SC
   SpectralCoordinate spec = getSpectralCoordinate (freqID);
-//
+
   Table t = table_.keywordSet().asTable("FREQUENCIES");
 
-// Set rest frequency
-
+  // Set rest frequency
   Vector<Double> restFreqIDs;
   t.keywordSet().get("RESTFREQS",restFreqIDs);
-  if (restFreqID < restFreqIDs.nelements()) {                  // SHould always be true
+  if (restFreqID < restFreqIDs.nelements()) {    // Should always be true
     spec.setRestFrequency(restFreqIDs(restFreqID),True);
   }
 
-// Set up frame conversion layer
-
+  // Set up frame conversion layer
   String frm;
   t.keywordSet().get("REFFRAME",frm);
@@ -835 +834 @@
   MFrequency::Types mtype;
   if (!MFrequency::getType(mtype, frm)) {
-    cerr << "Frequency type unknown assuming TOPO" << endl;       // SHould never happen
+    // Should never happen
+    cerr << "Frequency type unknown assuming TOPO" << endl;
     mtype = MFrequency::TOPO;
   }
 
-// Set reference frame conversion  (requires row)
-
+  // Set reference frame conversion  (requires row)
   MDirection dir = getDirection(whichRow);
   MEpoch epoch = getEpoch(whichRow);
   MPosition pos = getAntennaPosition();
-//
+
   if (!spec.setReferenceConversion(mtype,epoch,pos,dir)) {
     throw(AipsError("Couldn't convert frequency frame."));
   }
 
-// Now velocity conversion if appropriate
-
+  // Now velocity conversion if appropriate
   String unitStr;
   t.keywordSet().get("UNIT",unitStr);
-//
+
   String dpl;
   t.keywordSet().get("DOPPLER",dpl);
@@ -859 +857 @@
   MDoppler::getType(dtype, dpl);
 
-// Only set velocity unit if non-blank and non-Hz
-
+  // Only set velocity unit if non-blank and non-Hz
   if (!unitStr.empty()) {
      Unit unitU(unitStr);
@@ -868 +865 @@
      }
   }
-//
+
   return spec;
 }
@@ -937 +934 @@
   rfvec = q.getValue("Hz");
   aTable.rwKeywordSet().define("RESTFREQS", rfvec);
-  table_.rwKeywordSet().defineTable ("FREQUENCIES", aTable);
-  return True;
+  table_.rwKeywordSet().defineTable("FREQUENCIES", aTable);
+  return true;
+}
+
+bool SDMemTable::putSDFitTable(const SDFitTable& sdft)
+{
+  TableDesc td("", "1", TableDesc::Scratch);
+  td.addColumn(ArrayColumnDesc<String>("FUNCTIONS"));
+  td.addColumn(ArrayColumnDesc<Int>("COMPONENTS"));
+  td.addColumn(ArrayColumnDesc<Double>("PARAMETERS"));
+  td.addColumn(ArrayColumnDesc<Bool>("PARMASK"));
+  SetupNewTable aNewTab("fits", td, Table::New);
+  Table aTable(aNewTab, Table::Memory);
+  ArrayColumn<String> sc0(aTable, "FUNCTIONS");
+  ArrayColumn<Int> sc1(aTable, "COMPONENTS");
+  ArrayColumn<Double> sc2(aTable, "PARAMETERS");
+  ArrayColumn<Bool> sc3(aTable, "PARMASK");
+  for (uInt i; i<sdft.length(); ++i) {
+    const Vector<Double>& parms = sdft.getFitParameters(i);
+    const Vector<Bool>& parmask = sdft.getFitParameterMask(i);
+    const Vector<String>& funcs = sdft.getFitFunctions(i);
+    const Vector<Int>& comps = sdft.getFitComponents(i);
+    sc0.put(i,funcs);
+    sc1.put(i,comps);
+    sc3.put(i,parmask);
+    sc2.put(i,parms);
+  }
+  table_.rwKeywordSet().defineTable("FITS", aTable);
+  return true;
+}
+
+SDFitTable SDMemTable::getSDFitTable() const
+{
+  const Table& t = table_.keywordSet().asTable("FITS");
+  Vector<Double> parms;
+  Vector<Bool> parmask;
+  Vector<String> funcs;
+  Vector<Int> comps;  
+  ROArrayColumn<Double> parmsCol(t, "PARAMETERS");
+  ROArrayColumn<Bool> parmaskCol(t, "PARMASK");
+  ROArrayColumn<Int> compsCol(t, "COMPONENTS");
+  ROArrayColumn<String> funcsCol(t, "FUNCTIONS");
+  uInt n = t.nrow();
+  SDFitTable sdft(n);
+  for (uInt i=0; i<n; ++i) {
+    parmaskCol.get(i, parmask);
+    sdft.putFitParameterMask(i, parmask);
+    parmsCol.get(i, parms);
+    sdft.putFitParameters(i, parms);
+    funcsCol.get(i, funcs);
+    sdft.putFitFunctions(i, funcs);
+    compsCol.get(i, comps);
+    sdft.putFitComponents(i, comps);
+  }
+  return sdft;
+}
+
+
+void SDMemTable::addFit(uInt whichRow,
+                        const Vector<Double>& p, const Vector<Bool>& m,
+                        const Vector<String>& f, const Vector<Int>& c)
+{
+  if (whichRow >= nRow()) {
+    throw(AipsError("Specified row out of range"));
+  }
+  Table t = table_.keywordSet().asTable("FITS");
+  uInt nrow = t.nrow();  
+  t.addRow();
+  ArrayColumn<Double> parmsCol(t, "PARAMETERS");
+  ArrayColumn<Bool> parmaskCol(t, "PARMASK");
+  ArrayColumn<Int> compsCol(t, "COMPONENTS");
+  ArrayColumn<String> funcsCol(t, "FUNCTIONS");
+  parmsCol.put(nrow, p);
+  parmaskCol.put(nrow, m);
+  compsCol.put(nrow, c);
+  funcsCol.put(nrow, f);
+
+  Array<Int> fitarr;
+  fitCol_.get(whichRow, fitarr);
+
+  Array<Int> newarr;               // The new Array containing the fitid 
+  Int pos =-1;                     // The fitid position in the array 
+  if ( fitarr.nelements() == 0 ) { // no fits at all in this row
+    Array<Int> arr(IPosition(4,nBeam(),nIF(),nPol(),1));
+    arr = -1;
+    newarr.reference(arr);
+    pos = 0;      
+  } else {
+    IPosition shp = fitarr.shape();
+    IPosition start(4, beamSel_, IFSel_, polSel_,0);
+    IPosition end(4, beamSel_, IFSel_, polSel_, shp[3]-1);
+    // reform the output array slice to be of dim=1
+    Array<Int> tmp = (fitarr(start, end)).reform(IPosition(1,shp[3]));
+    const Vector<Int>& fits = tmp;
+    VectorSTLIterator<Int> it(fits);
+    Int i = 0;
+    while (it != fits.end()) {
+      if (*it == -1) {
+        pos = i;
+        break;
+      }
+      ++i;
+      ++it;
+    };
+  }
+  if (pos == -1) {
+      mathutil::extendLastArrayAxis(newarr,fitarr, -1);
+      pos = fitarr.shape()[3];     // the new element position
+  } else {
+    if (fitarr.nelements() > 0)
+      newarr = fitarr;
+  }
+  newarr(IPosition(4, beamSel_, IFSel_, polSel_, pos)) = Int(nrow); 
+  fitCol_.put(whichRow, newarr);
+
 }
 
@@ -946 +1056 @@
   SDFrequencyTable sdft;
 
-// Add refpix/refval/incr.  What are the units ? Hz I suppose
-// but it's nowhere described...
-
+  // Add refpix/refval/incr.  What are the units ? Hz I suppose
+  // but it's nowhere described...
   Vector<Double> refPix, refVal, incr;
   ScalarColumn<Double> refPixCol(t, "REFPIX");
@@ -956 +1065 @@
   refVal = refValCol.getColumn();
   incr = incrCol.getColumn();
-//
+
   uInt n = refPix.nelements();
   for (uInt i=0; i<n; i++) {
@@ -962 +1071 @@
   }
 
-// Frequency reference frame.  I don't know if this
-// is the correct frame.  It might be 'REFFRAME'
-// rather than 'BASEREFFRAME' ?
-
+  // Frequency reference frame.  I don't know if this
+  // is the correct frame.  It might be 'REFFRAME'
+  // rather than 'BASEREFFRAME' ?
   String baseFrame;
   t.keywordSet().get("BASEREFFRAME",baseFrame);
   sdft.setRefFrame(baseFrame);
 
-// Equinox
-
+  // Equinox  
   Float equinox;
   t.keywordSet().get("EQUINOX", equinox);
   sdft.setEquinox(equinox);
 
-// Rest Frequency
-
+  // Rest Frequency
   Vector<Double> restFreqs;
   t.keywordSet().get("RESTFREQS", restFreqs);
@@ -984 +1090 @@
   }
   sdft.setRestFrequencyUnit(String("Hz"));
-//
+
   return sdft;
 }
@@ -992 +1098 @@
   uInt rno = table_.nrow();
   table_.addRow();
-//
+
   timeCol_.put(rno, sdc.timestamp);
   srcnCol_.put(rno, sdc.sourcename);
@@ -1011 +1117 @@
   paraCol_.put(rno, sdc.parangle);
   histCol_.put(rno, sdc.getHistory());
-//
+  fitCol_.put(rno, sdc.getFitMap());
+
   return true;
 }
@@ -1031 +1138 @@
   sdc.tcal[0] = tc[0];sdc.tcal[1] = tc[1];
   tcaltCol_.get(whichRow, sdc.tcaltime);
-//
+
   Array<Float> spectrum;
   Array<Float> tsys;
@@ -1038 +1145 @@
   Array<Double> direction;
   Vector<String> histo;
-//
+  Array<Int> fits;
+  
   specCol_.get(whichRow, spectrum);
   sdc.putSpectrum(spectrum);
@@ -1053 +1161 @@
   histCol_.get(whichRow, histo);
   sdc.putHistory(histo);
+  fitCol_.get(whichRow, fits);
+  sdc.putFitMap(fits);
   return sdc;
 }
@@ -1168 +1278 @@
 {
   Bool throwIt = True;
-  Instrument ins = convertInstrument (name, throwIt);
+  Instrument ins = convertInstrument(name, throwIt);
   String nameU(name);
   nameU.upcase();
@@ -1176 +1286 @@
 std::string SDMemTable::summary(bool verbose) const  {
 
-// Format header info
-
+  // Format header info
   ostringstream oss;
   oss << endl;
@@ -1214 +1323 @@
       << "IF[" << getIF() << "] " << "Pol[" << getPol() << "]" << endl;
   oss << endl;
-//
+
   String dirtype ="Position ("+ MDirection::showType(getDirectionReference()) + ")";
   oss << setw(5) << "Scan"
@@ -1224 +1333 @@
   oss << "--------------------------------------------------------------------------------" << endl;
   
-// Generate list of scan start and end integrations
-
+  // Generate list of scan start and end integrations
   Vector<Int> scanIDs = scanCol_.getColumn();
   Vector<uInt> startInt, endInt;
   mathutil::scanBoundaries(startInt, endInt, scanIDs);
-//
+
   const uInt nScans = startInt.nelements();
   String name;
   Vector<uInt> freqIDs, listFQ;
   uInt nInt;
-//
+
   for (uInt i=0; i<nScans; i++) {
-
-// Get things from first integration of scan
-
-      String time = getTime(startInt(i),False);
-      String tInt = formatSec(Double(getInterval(startInt(i))));
-      String posit = formatDirection(getDirection(startInt(i),True));
-      srcnCol_.getScalar(startInt(i),name);
-
-// Find all the FreqIDs in this scan
-
-      listFQ.resize(0);      
-      for (uInt j=startInt(i); j<endInt(i)+1; j++) {
-         freqidCol_.get(j, freqIDs);
-         for (uInt k=0; k<freqIDs.nelements(); k++) {
-            mathutil::addEntry(listFQ, freqIDs(k));
-         }
+    // Get things from first integration of scan
+    String time = getTime(startInt(i),False);
+    String tInt = formatSec(Double(getInterval(startInt(i))));
+    String posit = formatDirection(getDirection(startInt(i),True));
+    srcnCol_.getScalar(startInt(i),name);
+
+    // Find all the FreqIDs in this scan
+    listFQ.resize(0);      
+    for (uInt j=startInt(i); j<endInt(i)+1; j++) {
+      freqidCol_.get(j, freqIDs);
+      for (uInt k=0; k<freqIDs.nelements(); k++) {
+        mathutil::addEntry(listFQ, freqIDs(k));
       }
-//
-      nInt = endInt(i) - startInt(i) + 1;
-      oss << setw(3) << std::right << i << std::left << setw(2) << "  "
-          << setw(15) << name
-          << setw(24) << posit
-          << setw(10) << time 
-          << setw(3) << std::right << nInt  << setw(3) << " x " << std::left
-          << setw(6) <<  tInt
-          << " " << listFQ << endl;
+    }
+
+    nInt = endInt(i) - startInt(i) + 1;
+    oss << setw(3) << std::right << i << std::left << setw(2) << "  "
+        << setw(15) << name
+        << setw(24) << posit
+        << setw(10) << time 
+        << setw(3) << std::right << nInt  << setw(3) << " x " << std::left
+        << setw(6) <<  tInt
+        << " " << listFQ << endl;
   }
   oss << endl;
-  oss << "Table contains " << table_.nrow() << " integration(s) in " << nScans << " scan(s)." << endl;
-
-// Frequency Table
-
+  oss << "Table contains " << table_.nrow() << " integration(s) in " 
+      << nScans << " scan(s)." << endl;
+
+  // Frequency Table
   if (verbose) {
     std::vector<string> info = getCoordInfo();
@@ -1292 +1397 @@
   return n;
 }
+
 Int SDMemTable::nIF() const {
   Int n;
@@ -1297 +1403 @@
   return n;
 }
+
 Int SDMemTable::nPol() const {
   Int n;
@@ -1302 +1409 @@
   return n;
 }
+
 Int SDMemTable::nChan() const {
   Int n;
@@ -1307 +1415 @@
   return n;
 }
+
 bool SDMemTable::appendHistory(const std::string& hist, int whichRow)
 {
@@ -1378 +1487 @@
     cerr  << "Unknown equinox using J2000" << endl;
   }
-//
+
   return mdr;
 }
@@ -1391 +1500 @@
     met = MEpoch::UTC;
   }
-//
+
   return met;
 }
@@ -1402 +1511 @@
    t.upcase();
 
-// The strings are what SDReader returns, after cunning interrogation
-// of station names... :-(
-
+   // The strings are what SDReader returns, after cunning interrogation
+   // of station names... :-(
    Instrument inst = asap::UNKNOWN;
    if (t==String("DSS-43")) {                
@@ -1426 +1534 @@
 }
 
-Bool SDMemTable::setRestFreqs (const Vector<Double>& restFreqsIn, 
-                                     const String& sUnit,
-                                     const vector<string>& lines,
-                                     const String& source,
-                                     Int whichIF)
+Bool SDMemTable::setRestFreqs(const Vector<Double>& restFreqsIn, 
+                              const String& sUnit,
+                              const vector<string>& lines,
+                              const String& source,
+                              Int whichIF)
 {
    const Int nIFs = nIF();
@@ -1437 +1545 @@
    }
 
-// FInd vector of restfrequencies
-// Double takes precedence over String
-
+   // Find vector of restfrequencies
+   // Double takes precedence over String
    Unit unit;
    Vector<Double> restFreqs;
@@ -1457 +1564 @@
             restFreqs[i] = lineFreq.getValue().getValue();          // Hz
          } else {
-            String s = String(lines[i]) + String(" is an unrecognized spectral line");
+            String s = String(lines[i]) + 
+              String(" is an unrecognized spectral line");
             throw(AipsError(s));
          }
@@ -1465 +1573 @@
    }
 
-// If multiple restfreqs, must be length nIF. In this
-// case we will just replace the rest frequencies
-// 
-
+   // If multiple restfreqs, must be length nIF. In this
+   // case we will just replace the rest frequencies
    const uInt nRestFreqs = restFreqs.nelements();
    Int idx = -1;
@@ -1474 +1580 @@
 
    if (nRestFreqs>1) {
-
-// Replace restFreqs, one per IF
-
+     // Replace restFreqs, one per IF
       if (nRestFreqs != nIFs) {
          throw (AipsError("Number of rest frequencies must be equal to the number of IFs"));
       }
       cout << "Replacing rest frequencies with given list, one per IF" << endl;
-//
       sdft.deleteRestFrequencies();
       for (uInt i=0; i<nRestFreqs; i++) {
@@ -1489 +1592 @@
    } else {
 
-// Add new rest freq
-
+     // Add new rest freq
       Quantum<Double> rf(restFreqs[0], unit);
       idx = sdft.addRestFrequency(rf.getValue("Hz"));
       cout << "Selecting given rest frequency" << endl;
    }
-
-// Replace
-
+   
+   // Replace
    Bool empty = source.empty();
    Bool ok = False;
@@ -1507 +1608 @@
          srcnCol_.get(i, srcName);
          restfreqidCol_.get(i,restFreqIDs);       
-// 
+
          if (idx==-1) {
-
-// Replace vector of restFreqs; one per IF. 
-// No selection possible
-
+           // Replace vector of restFreqs; one per IF.
+           // No selection possible
             for (uInt i=0; i<nIFs; i++) restFreqIDs[i] = i;
          } else {
-
-// Set RestFreqID for selected data
-
+           // Set RestFreqID for selected data
             if (empty || source==srcName) {
                if (whichIF<0) {
@@ -1526 +1623 @@
             }
          }
-//
          restfreqidCol_.put(i,restFreqIDs);       
       }
       ok = True;
    } else {
-      ok = False;
+     ok = False;
    }
-//
+
    return ok;
 }
@@ -1542 +1638 @@
    MFrequency lineFreq;
    Double freq;
-//
+
    cout.flags(std::ios_base::left);
    cout << "Line      Frequency (Hz)" << endl;
@@ -1549 +1645 @@
      MeasTable::Line(lineFreq, lines[i]);
      freq = lineFreq.getValue().getValue();          // Hz
-//
      cout << setw(11) << lines[i] << setprecision(10) << freq << endl;
    }
@@ -1577 +1672 @@
 
 void SDMemTable::rotateXYPhase (Float value, Bool doAll)
-//
-// phase in degrees
-// Applies to all Beams and IFs
-// Might want to optionally select on Beam/IF
-//
-{
-   if (nPol() != 4) {
-      throw(AipsError("You must have 4 polarizations to run this function"));
-   }
-//
-   IPosition start(asap::nAxes,0);
-   IPosition end(asap::nAxes);
-//
-   uInt nRow = specCol_.nrow();
-   Array<Float> data;
-   for (uInt i=0; i<nRow;++i) {
-      specCol_.get(i,data);
-      IPosition shape = data.shape();
-
-// Set slice
-
-      if (!doAll) {
-         setCursorSlice (start, end, shape);
-      } else {
-         end = shape-1;
-      }
-
-// Get polarization slice references
-
-      start(asap::PolAxis) = 2;
-      end(asap::PolAxis) = 2;
-      Array<Float> C3 = data(start,end);
-//
-      start(asap::PolAxis) = 3;
-      end(asap::PolAxis) = 3;
-      Array<Float> C4 = data(start,end);
-
-// Rotate
-
-      SDPolUtil::rotateXYPhase(C3, C4, value);
-
-// Put
-
-      specCol_.put(i,data);
-   }
-}
-
-
-void SDMemTable::setCursorSlice (IPosition& start, IPosition& end,
-                                 const IPosition& shape) const
-{
-   const uInt nDim = shape.nelements();
-   start.resize(nDim);
-   end.resize(nDim);
-//
-   start(asap::BeamAxis) = beamSel_;
-   end(asap::BeamAxis) = beamSel_;
-//
-   start(asap::IFAxis) = IFSel_;
-   end(asap::IFAxis) = IFSel_;
-//
-   start(asap::PolAxis) = polSel_;
-   end(asap::PolAxis) = polSel_;
-//
-   start(asap::ChanAxis) = 0;
-   end(asap::ChanAxis) = shape(asap::ChanAxis) - 1;
-}
-
-
-std::vector<float> SDMemTable::convertVector (const Vector<Float>& in) const
-{
-   std::vector<float> out(in.nelements());
-   for (uInt i=0; i<in.nelements(); i++) {
-      out[i] = in[i];
-   }
-   return out;
-}
-
-
-std::vector<float> SDMemTable::getFloatSpectrum (const Array<Float>& arr) const
-//
-// Get spectrum at cursor location
-//
-{
-
-// Setup accessors
-
+  // phase in degrees
+  // Applies to all Beams and IFs
+  // Might want to optionally select on Beam/IF
+{
+  if (nPol() != 4) {
+    throw(AipsError("You must have 4 polarizations to run this function"));
+  }
+  
+  IPosition start(asap::nAxes,0);
+  IPosition end(asap::nAxes);
+  
+  uInt nRow = specCol_.nrow();
+  Array<Float> data;
+  for (uInt i=0; i<nRow;++i) {
+    specCol_.get(i,data);
+    IPosition shape = data.shape();
+    
+    // Set slice
+    if (!doAll) {
+      setCursorSlice (start, end, shape);
+    } else {
+      end = shape-1;
+    }
+
+    // Get polarization slice references
+    start(asap::PolAxis) = 2;
+    end(asap::PolAxis) = 2;
+    Array<Float> C3 = data(start,end);
+    
+    start(asap::PolAxis) = 3;
+    end(asap::PolAxis) = 3;
+    Array<Float> C4 = data(start,end);
+    
+    // Rotate
+    SDPolUtil::rotateXYPhase(C3, C4, value);
+    
+    // Put
+    specCol_.put(i,data);
+  }
+}
+
+void SDMemTable::setCursorSlice(IPosition& start, IPosition& end,
+                                const IPosition& shape) const
+{
+  const uInt nDim = shape.nelements();
+  start.resize(nDim);
+  end.resize(nDim);
+  
+  start(asap::BeamAxis) = beamSel_;
+  end(asap::BeamAxis) = beamSel_;
+  start(asap::IFAxis) = IFSel_;
+  end(asap::IFAxis) = IFSel_;
+
+  start(asap::PolAxis) = polSel_;
+  end(asap::PolAxis) = polSel_;
+
+  start(asap::ChanAxis) = 0;
+  end(asap::ChanAxis) = shape(asap::ChanAxis) - 1;
+}
+
+
+std::vector<float> SDMemTable::getFloatSpectrum(const Array<Float>& arr) const
+  // Get spectrum at cursor location
+{
+
+  // Setup accessors
   ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
   aa0.reset(aa0.begin(uInt(beamSel_)));                    // Beam selection
@@ -1672 +1746 @@
   ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
   aa2.reset(aa2.begin(uInt(polSel_)));                     // Pol selection
-//
+  
   std::vector<float> spectrum;
   for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {

trunk/src/SDMemTable.h

@@ -36 +36 @@
 #include <vector>
 // AIPS++
- #include <casa/aips.h>
+#include <casa/aips.h>
 #include <casa/Arrays/MaskedArray.h>
 #include <casa/BasicSL/String.h>
@@ -49 +49 @@
 
 namespace asap {
-
+  
 class SDContainer;
 class SDHeader;
 class SDFrequencyTable;
-
+class SDFitTable;
 
 
@@ -95 +95 @@
   virtual std::vector<bool> getMask(casa::Int whichRow=0) const;
 
-
   // When we can handle input correlations being either circular or
-  // linear, we should probably have functions; getLinear, getCircular, getStokes
-  // since one can inter-convert between all three regardless of the input
-  // provided there are 4 polarizations. For now, we are assuming, if
-  // there are 4 polarizations, that we have linears.  getSpectrum 
-  // is then 'getLinear', getStokesSpectrum is 'getStokes' and getCircularSpectrum
-  // is 'getCircular'  
-
-
-  // Get Stokes at cursor location. One of either I,Q,U,V or I,P,PA,V (doPol=True)
-  // If the latter, you can add a PA offset (degrees)
+  // linear, we should probably have functions; getLinear,
+  // getCircular, getStokes since one can inter-convert between all
+  // three regardless of the input provided there are 4
+  // polarizations. For now, we are assuming, if there are 4
+  // polarizations, that we have linears.  getSpectrum is then
+  // 'getLinear', getStokesSpectrum is 'getStokes' and
+  // getCircularSpectrum is 'getCircular'
+
+
+  // Get Stokes at cursor location. One of either I,Q,U,V or I,P,PA,V
+  // (doPol=True) If the latter, you can add a PA offset (degrees)
   virtual std::vector<float> getStokesSpectrum(casa::Int whichRow=0, 
                                                casa::Bool doPol=casa::False,
@@ -130 +130 @@
   casa::MDirection getDirection(casa::Int whichRow=0,
                                 casa::Bool refBeam=casa::False) const; 
-//
+
   std::string getSourceName(casa::Int whichRow=0) const;
   double getInterval(casa::Int whichRow=0) const;
@@ -137 +137 @@
   virtual void setCoordInfo(std::vector<string> theinfo);
 
-// Set RestFreqID.  source="" and IF=-1 means select all
-  virtual casa::Bool setRestFreqs (const casa::Vector<casa::Double>& restFreqs, 
-                                   const casa::String& unit,
-                                   const std::vector<std::string>& lines,
-                                   const casa::String& source,
-                                   casa::Int whichIF=-1);
-
-// List lines
-   void spectralLines() const;
-
-// Get/Set flux unit
+  // Set RestFreqID.  source="" and IF=-1 means select all
+  virtual casa::Bool setRestFreqs(const casa::Vector<casa::Double>& restFreqs, 
+                                  const casa::String& unit,
+                                  const std::vector<std::string>& lines,
+                                  const casa::String& source,
+                                  casa::Int whichIF=-1);
+
+  // List lines
+  void spectralLines() const;
+
+  // Get/Set flux unit
   std::string getFluxUnit() const;
   void setFluxUnit (const std::string& unit);
-
-// Set Instrument
+  
+  // Set Instrument
   void setInstrument (const std::string& instrument);
 
@@ -243 +243 @@
   static Instrument convertInstrument(const casa::String& instrument,
                                       casa::Bool throwIt);
-
+  
+  bool putSDFitTable(const SDFitTable& sdft);
+  SDFitTable getSDFitTable() const;
+
+  void addFit(casa::uInt whichRow,
+              const casa::Vector<casa::Double>& p, 
+              const casa::Vector<casa::Bool>& m,
+              const casa::Vector<casa::String>& f, 
+              const casa::Vector<casa::Int>& c);
+  
 private:
   // utility func for nice printout
   casa::String formatSec(casa::Double x) const;
   casa::String formatDirection(const casa::MDirection& md) const;
-  std::vector<float> getFloatSpectrum (const casa::Array<casa::Float>& arr) const;
+  std::vector<float> getFloatSpectrum(const casa::Array<casa::Float>& arr) const;
   void setup();
   void attach();
@@ -254 +263 @@
 
   // Generate start and end for shape and current cursor selection
-  void setCursorSlice(casa::IPosition& start, casa::IPosition& end, const casa::IPosition& shape) const;
-
-  // Convert Vector to vector
-  std::vector<float> convertVector (const casa::Vector<casa::Float>& in) const;
-
+  void setCursorSlice(casa::IPosition& start, casa::IPosition& end, 
+                      const casa::IPosition& shape) const;
 
   // the current cursor into the array
@@ -266 +272 @@
   casa::Table table_;
 
-// Cached Columns to avoid reconstructing them for each row get/put
+  // Cached Columns to avoid reconstructing them for each row get/put
   casa::ScalarColumn<casa::Double> timeCol_, integrCol_;
   casa::ScalarColumn<casa::Float> azCol_, elCol_, paraCol_;
@@ -276 +282 @@
   casa::ArrayColumn<casa::uInt> freqidCol_, restfreqidCol_;
   casa::ArrayColumn<casa::String> histCol_;
+  casa::ArrayColumn<casa::Int> fitCol_;
   casa::ROArrayColumn<casa::Float> stokesCol_;
 };

trunk/src/SDMemTableWrapper.h

@@ -81 +81 @@
   }
 
-  std::vector<float> getStokesSpectrum(int whichRow=0, bool doPol=false, float paOff=0.0) const {
+  std::vector<float> getStokesSpectrum(int whichRow=0, bool doPol=false, 
+                                       float paOff=0.0) const {
     return table_->getStokesSpectrum(whichRow, doPol, paOff);
   }
 
-  std::vector<float> getCircularSpectrum(int whichRow=0, bool doRR=true) const {
+  std::vector<float> getCircularSpectrum(int whichRow=0, 
+                                         bool doRR=true) const {
     return table_->getCircularSpectrum(whichRow, doRR);
   }
@@ -187 +189 @@
   }
 
+  void addFit(int whichRow, const std::vector<double>& p,
+              const std::vector<bool>& m, const std::vector<string>& f,
+              const std::vector<int>& c) {
+    
+    casa::Vector<casa::Double> p2(p);
+    casa::Vector<casa::Bool> m2(m); 
+    casa::Vector<casa::String> f2(f.size());
+    casa::uInt i=0;
+    std::vector<std::string>::const_iterator it;
+    for (it=f.begin();it != f.end();++it) {
+      f2[i] = casa::String(*it);
+    }
+    casa::Vector<casa::Int> c2(c);
+    table_->addFit(casa::uInt(whichRow), p2,m2,f2,c2);
+  }
+
   void rotateXYPhase (float value, bool doAll) {
       table_->rotateXYPhase(value, doAll);