- Timestamp:
- 02/16/06 12:02:18 (19 years ago)
- Location:
- trunk/src
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/STFiller.cpp
r804 r805 1 //#--------------------------------------------------------------------------- 2 //# SDReader.cc: A class to read single dish spectra from SDFITS, RPFITS 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 //#--------------------------------------------------------------------------- 1 // 2 // C++ Implementation: STFiller 3 // 4 // Description: 5 // 6 // 7 // Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006 8 // 9 // Copyright: See COPYING file that comes with this distribution 10 // 11 // 31 12 #include <casa/iostream.h> 32 13 #include <casa/iomanip.h> … … 36 17 #include <casa/OS/File.h> 37 18 #include <casa/Quanta/Unit.h> 19 #include <casa/Arrays/ArrayMath.h> 20 #include <casa/Utilities/Regex.h> 21 22 #include <casa/Containers/RecordField.h> 23 24 #include <tables/Tables/TableRow.h> 25 38 26 #include <atnf/PKSIO/PKSreader.h> 39 27 40 #include "SDReader.h"41 28 #include "SDDefs.h" 42 29 #include "SDAttr.h" 43 30 31 #include "STFiller.h" 32 44 33 using namespace casa; 45 using namespace asap; 46 47 SDReader::SDReader() : 34 35 namespace asap { 36 37 STFiller::STFiller() : 48 38 reader_(0), 49 39 header_(0), 50 frequencies_(0), 51 table_(new SDMemTable()), 52 haveXPol_(False) 53 { 54 cursor_ = 0; 55 } 56 SDReader::SDReader(const std::string& filename, 57 int whichIF, int whichBeam) : 40 table_(0) 41 { 42 } 43 44 STFiller::STFiller( CountedPtr< Scantable > stbl ) : 58 45 reader_(0), 59 46 header_(0), 60 frequencies_(0), 61 table_(new SDMemTable()), 62 haveXPol_(False) 63 { 64 cursor_ = 0; 65 open(filename, whichIF, whichBeam); 66 } 67 68 SDReader::SDReader(CountedPtr<SDMemTable> tbl) : 47 table_(stbl) 48 { 49 } 50 51 STFiller::STFiller(const std::string & filename, 52 int whichIF, int whichBeam ) : 69 53 reader_(0), 70 54 header_(0), 71 table_(tbl), 72 haveXPol_(False) 73 { 74 cursor_ = 0; 75 } 76 77 SDReader::~SDReader() { 78 if (reader_) delete reader_; 79 if (header_) delete header_; 80 if (frequencies_) delete frequencies_; 81 } 82 83 void SDReader::reset() { 84 cursor_ = 0; 85 table_ = new SDMemTable(); 86 open(filename_,ifOffset_, beamOffset_); 87 } 88 89 90 void SDReader::close() { 91 //cerr << "disabled" << endl; 92 } 93 94 void SDReader::open(const std::string& filename, 95 int whichIF, int whichBeam) { 96 if (reader_) delete reader_; reader_ = 0; 97 Bool haveBase, haveSpectra; 55 table_(0) 56 { 57 open(filename, whichIF, whichBeam); 58 } 59 60 STFiller::~STFiller() 61 { 62 close(); 63 } 64 65 void STFiller::open( const std::string & filename, 66 int whichIF, int whichBeam ) 67 { 68 if (table_.null()) table_ = new Scantable(); 69 if (reader_) { delete reader_; reader_ = 0; } 70 Bool haveBase, haveSpectra; 98 71 99 72 String inName(filename); … … 102 75 103 76 File file(inName); 104 if ( !file.exists()) {77 if ( !file.exists() ) { 105 78 throw(AipsError("File does not exist")); 106 79 } … … 110 83 String format; 111 84 Vector<Bool> beams; 112 if ( (reader_ = getPKSreader(inName, 0, False, format, beams, nIF_,85 if ( (reader_ = getPKSreader(inName, 0, 0, format, beams, nIF_, 113 86 nChan_, nPol_, haveBase, haveSpectra, 114 haveXPol_)) == 0 ) {87 haveXPol_)) == 0 ) { 115 88 throw(AipsError("Creation of PKSreader failed")); 116 89 } … … 124 97 nBeam_ = beams.nelements(); 125 98 // Get basic parameters. 126 if ( haveXPol_) {99 if ( haveXPol_ ) { 127 100 pushLog("Cross polarization present"); 128 101 nPol_ += 2; // Convert Complex -> 2 Floats … … 145 118 delete reader_; 146 119 reader_ = 0; 120 delete header_; 121 header_ = 0; 147 122 throw(AipsError("Failed to get header.")); 148 return;149 123 } 150 124 if ((header_->obstype).matches("*SW*")) { … … 153 127 "setting # of IFs = 1 "); 154 128 nIF_ = 1; 129 header_->obstype = String("fswitch"); 155 130 } 156 131 … … 166 141 header_->nif = nIF_; 167 142 header_->epoch = "UTC"; 168 143 // *** header_->frequnit = "Hz" 169 144 // Apply selection criteria. 170 145 … … 176 151 if (whichIF>=0) { 177 152 if (whichIF>=0 && whichIF<nIF_) { 178 179 180 181 182 153 IFsel = False; 154 IFsel(whichIF) = True; 155 header_->nif = 1; 156 nIF_ = 1; 157 ifOffset_ = whichIF; 183 158 } else { 184 throw(AipsError("Illegal IF selection")); 159 delete reader_; 160 reader_ = 0; 161 delete header_; 162 header_ = 0; 163 throw(AipsError("Illegal IF selection")); 185 164 } 186 165 } … … 188 167 beamOffset_ = 0; 189 168 if (whichBeam>=0) { 190 if (whichBeam>=0 && whichBeam<nBeam_) { 191 beamSel = False; 192 beamSel(whichBeam) = True; 193 header_->nbeam = 1; 194 nBeam_ = 1; 195 beamOffset_ = whichBeam; 196 } else { 197 throw(AipsError("Illegal Beam selection")); 198 } 169 if (whichBeam>=0 && whichBeam<nBeam_) { 170 beamSel = False; 171 beamSel(whichBeam) = True; 172 header_->nbeam = 1; 173 nBeam_ = 1; 174 beamOffset_ = whichBeam; 175 } else { 176 delete reader_; 177 reader_ = 0; 178 delete header_; 179 header_ = 0; 180 throw(AipsError("Illegal Beam selection")); 181 } 199 182 } 200 183 Vector<Int> start(nIF_, 1); … … 202 185 reader_->select(beamSel, IFsel, start, end, ref, True, haveXPol_); 203 186 table_->putSDHeader(*header_); 204 205 if (frequencies_) delete frequencies_; 206 frequencies_ = new SDFrequencyTable(); 207 frequencies_->setRefFrame(header_->freqref); 208 frequencies_->setBaseRefFrame(header_->freqref); 209 frequencies_->setRestFrequencyUnit("Hz"); 210 frequencies_->setEquinox(header_->equinox); 211 } 212 213 int SDReader::read(const std::vector<int>& seq) { 187 } 188 189 void STFiller::close( ) 190 { 191 delete reader_;reader_=0; 192 delete header_;header_=0; 193 table_ = 0; 194 } 195 196 int asap::STFiller::read( ) 197 { 214 198 int status = 0; 215 199 … … 226 210 Complex xCalFctr; 227 211 Vector<Complex> xPol; 228 uInt n = seq.size(); 229 230 231 uInt stepsize = header_->nif*header_->nbeam; 232 uInt seqi = 0; 233 Bool getAll = False; 234 235 if (seq[n-1] == -1) getAll = True; 236 while ( (cursor_ <= seq[n-1]) || getAll) { 237 // only needs to be create if in seq 238 SDContainer sc(header_->nbeam,header_->nif,header_->npol,header_->nchan); 239 // iterate over one correlator integration cycle = nBeam*nIF 240 for (uInt row=0; row < stepsize; row++) { 241 // stepsize as well 242 // spectra(nChan,nPol)!!! 243 status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName, 244 srcName, srcDir, srcPM, srcVel, IFno, refFreq, 245 bandwidth, freqInc, restFreq, tcal, tcalTime, 246 azimuth, elevation, parAngle, focusAxi, 247 focusTan, focusRot, temperature, pressure, 248 humidity, windSpeed, windAz, refBeam, 249 beamNo, direction, scanRate, 250 tsys, sigma, calFctr, baseLin, baseSub, 251 spectra, flagtra, xCalFctr, xPol); 252 253 // Make sure beam/IF numbers are 0-relative - dealing with 254 // possible IF or Beam selection 255 beamNo = beamNo - beamOffset_ - 1; 256 IFno = IFno - ifOffset_ - 1; 257 258 if (status) { 259 if (status == -1) { 260 // EOF. 261 if (row > 0 && row < stepsize-1) 262 pushLog("incomplete scan data.\n Probably means not all Beams/IFs/Pols within a scan are present."); 263 264 // flush frequency table 265 table_->putSDFreqTable(*frequencies_); 266 return status; 267 } 268 } 269 // if in the given list 270 if (cursor_ == seq[seqi] || getAll) { 271 // add integration cycle 272 if (row==0) { 273 //add invariant info: scanNo, mjd, interval, fieldName, 274 //srcName, azimuth, elevation, parAngle, focusAxi, focusTan, 275 //focusRot, temperature, pressure, humidity, windSpeed, 276 //windAz srcDir, srcPM, srcVel 277 sc.timestamp = mjd; 278 sc.interval = interval; 279 sc.sourcename = srcName; 280 sc.fieldname = fieldName; 281 sc.azimuth = azimuth; 282 sc.elevation = elevation; 283 } 284 // add specific info 285 // refPix = nChan/2+1 in 1-rel Integer arith.! 286 Int refPix = header_->nchan/2; // 0-rel 287 uInt freqID = frequencies_->addFrequency(refPix, refFreq, freqInc); 288 uInt restFreqID = frequencies_->addRestFrequency(restFreq); 289 290 sc.setFrequencyMap(freqID, IFno); 291 sc.setRestFrequencyMap(restFreqID, IFno); 292 293 sc.tcal[0] = tcal[0];sc.tcal[1] = tcal[1]; 294 sc.tcaltime = tcalTime; 295 sc.parangle = parAngle; 296 sc.refbeam = -1; //nbeams == 1 297 if (nBeam_ > 1) // circumvent a bug "asap0000" in read which 298 // returns a random refbema number on multiple 299 // reads 300 sc.refbeam = refBeam-1;//make it 0-based; 301 sc.scanid = scanNo-1;//make it 0-based 302 if (haveXPol_) { 303 sc.setSpectrum(spectra, xPol, beamNo, IFno); 304 sc.setFlags(flagtra, beamNo, IFno, True); 305 } else { 306 sc.setSpectrum(spectra, beamNo, IFno); 307 sc.setFlags(flagtra, beamNo, IFno, False); 308 } 309 sc.setTsys(tsys, beamNo, IFno, haveXPol_); 310 sc.setDirection(direction, beamNo); 311 } 312 } 313 314 if (cursor_ == seq[seqi] || getAll) { 315 // insert container into table/list 316 table_->putSDContainer(sc); 317 seqi++;// next in list 318 } 319 cursor_++;// increment position in file 320 321 } 322 table_->putSDFreqTable(*frequencies_); 212 while ( status == 0 ) { 213 status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName, 214 srcName, srcDir, srcPM, srcVel, IFno, refFreq, 215 bandwidth, freqInc, restFreq, tcal, tcalTime, 216 azimuth, elevation, parAngle, focusAxi, 217 focusTan, focusRot, temperature, pressure, 218 humidity, windSpeed, windAz, refBeam, 219 beamNo, direction, scanRate, 220 tsys, sigma, calFctr, baseLin, baseSub, 221 spectra, flagtra, xCalFctr, xPol); 222 if ( status != 0 ) break; 223 TableRow row(table_->table()); 224 TableRecord& rec = row.record(); 225 RecordFieldPtr<uInt> scanoCol(rec, "SCANNO"); 226 *scanoCol = scanNo-1; 227 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO"); 228 *cyclenoCol = cycleNo-1; 229 RecordFieldPtr<Double> mjdCol(rec, "TIME"); 230 *mjdCol = mjd; 231 RecordFieldPtr<Double> intCol(rec, "INTERVAL"); 232 *intCol = interval; 233 RecordFieldPtr<String> srcnCol(rec, "SRCNAME"); 234 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE"); 235 // try to auto-identify if it is on or off. 236 Regex rx(".*_[e|w|R]$"); 237 Regex rx2("_[e|w|R|S]$"); 238 Int match = srcName.matches(rx); 239 if (match) { 240 *srcnCol = srcName; 241 } else { 242 *srcnCol = srcName.before(rx2); 243 } 244 //*srcnCol = srcName;//.before(rx2); 245 *srctCol = match; 246 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO"); 247 *beamCol = beamNo-beamOffset_-1; 248 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO"); 249 Int rb = -1; 250 if (nBeam_ > 1 ) rb = refBeam-1; 251 *rbCol = rb; 252 RecordFieldPtr<uInt> ifCol(rec, "IFNO"); 253 *ifCol = IFno-ifOffset_- 1; 254 uInt id; 255 /// @todo this has to change when nchan isn't global anymore 256 id = table_->frequencies().addEntry(Double(header_->nchan/2), 257 refFreq, freqInc); 258 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID"); 259 *mfreqidCol = id; 260 261 id = table_->molecules().addEntry(restFreq); 262 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID"); 263 *molidCol = id; 264 265 id = table_->tcal().addEntry(tcalTime, tcal); 266 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID"); 267 *mcalidCol = id; 268 id = table_->weather().addEntry(temperature, pressure, humidity, 269 windSpeed, windAz); 270 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID"); 271 *mweatheridCol = id; 272 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID"); 273 id = table_->focus().addEntry(focusAxi, focusTan, focusRot); 274 *mfocusidCol = id; 275 RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION"); 276 *dirCol = direction; 277 RecordFieldPtr<Double> azCol(rec, "AZIMUTH"); 278 *azCol = azimuth; 279 RecordFieldPtr<Double> elCol(rec, "ELEVATION"); 280 *elCol = elevation; 281 282 RecordFieldPtr<Float> parCol(rec, "PARANGLE"); 283 *parCol = parAngle; 284 285 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA"); 286 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA"); 287 RecordFieldPtr< uInt > polnoCol(rec, "POLNO"); 288 289 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS"); 290 // Turn the (nchan,npol) matrix and possible complex xPol vector 291 // into 2-4 rows in the scantable 292 Vector<Float> tsysvec(1); 293 // Why is spectra.ncolumn() == 3 for haveXPol_ == True 294 uInt npol = (spectra.ncolumn()==1 ? 1: 2); 295 for ( uInt i=0; i< npol; ++i ) { 296 tsysvec = tsys(i); 297 *tsysCol = tsysvec; 298 *polnoCol = i; 299 300 *specCol = spectra.column(i); 301 *flagCol = flagtra.column(i); 302 table_->table().addRow(); 303 row.put(table_->table().nrow()-1, rec); 304 } 305 if ( haveXPol_ ) { 306 // no tsys given for xpol, so emulate it 307 tsysvec = sqrt(tsys[0]*tsys[1]); 308 *tsysCol = tsysvec; 309 // add real part of cross pol 310 *polnoCol = 2; 311 Vector<Float> r(real(xPol)); 312 *specCol = r; 313 // make up flags from linears 314 /// @fixme this has to be a bitwise or of both pols 315 *flagCol = flagtra.column(0);// | flagtra.column(1); 316 table_->table().addRow(); 317 row.put(table_->table().nrow()-1, rec); 318 // ad imaginary part of cross pol 319 *polnoCol = 3; 320 Vector<Float> im(imag(xPol)); 321 *specCol = im; 322 table_->table().addRow(); 323 row.put(table_->table().nrow()-1, rec); 324 } 325 } 326 if (status > 0) { 327 close(); 328 throw(AipsError("Reading error occured, data possibly corrupted.")); 329 } 330 323 331 return status; 324 332 } 333 334 }//namespace asap -
trunk/src/STFiller.h
r804 r805 1 //#--------------------------------------------------------------------------- 2 //# SDReader.h: A class to read single dish spectra from SDFITS, RPFITS 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 #ifndef SDREADER_H 32 #define SDREADER_H 1 // 2 // C++ Interface: STFiller 3 // 4 // Description: 5 // 6 // 7 // Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006 8 // 9 // Copyright: See COPYING file that comes with this distribution 10 // 11 // 12 #ifndef STFILLER_H 13 #define STFILLER_H 33 14 34 15 #include <vector> … … 41 22 #include <casa/Arrays/Vector.h> 42 23 43 #include "S DMemTable.h"24 #include "Scantable.h" 44 25 #include "SDContainer.h" 45 26 #include "SDLog.h" … … 49 30 namespace asap { 50 31 51 class SDReader : public SDLog { 32 /** 33 This class fills a Scantable from external data formats using the PKSReader class. 34 35 @author Malte Marquarding 36 @date 2006/01/16 37 @version 2.0a 38 */ 39 class STFiller : public SDLog { 52 40 public: 53 SDReader();54 SDReader(const std::string& filename,55 int whichIF=-1, int whichBeam=-1);56 SDReader(casa::CountedPtr<SDMemTable> tbl);57 virtual ~SDReader();58 41 59 void open(const std::string& filename, 60 int whichIF=-1, 61 int whichBeam=-1); 62 void close(); 63 int read(const std::vector<int>& seq); 42 STFiller(); 64 43 65 casa::CountedPtr<SDMemTable> getTable() const { return table_;} 66 67 void reset(); 68 69 std::vector<int> pseudoHeader() const { 70 std::vector<int> v; 71 v.push_back(nBeam_);v.push_back(nIF_); 72 v.push_back(nPol_);v.push_back(nChan_); 73 return v; 74 } 44 STFiller(casa::CountedPtr< Scantable > stbl); 75 45 76 46 77 protected: 78 47 /** 48 * A constructor for a filler with associated input file 49 * @param filename the input file (rpf,sdfite or ms) 50 * @param whichIF read a specific IF only (default -1 means all IFs) 51 * @param whichBeam read a specific beam only (default -1 means all beams) 52 */ 53 STFiller( const std::string& filename, int whichIF=-1, 54 int whichBeam=-1 ); 55 56 ~STFiller(); 57 58 /** 59 * associate the Filler with a file on disk 60 * @param filename the input file (rpf,sdfite or ms) 61 * @param whichIF read a specific IF only (default -1 means all IFs) 62 * @param whichBeam read a specific beam only (default -1 means all beams) 63 * @exception AipsError Creation of PKSreader failed 64 */ 65 void open( const std::string& filename, int whichIF=-1, int whichBeam=-1 ); 66 67 /** 68 * detatch from file 69 */ 70 void close( ); 71 72 /** 73 * Read in "rows" from the source file attached with open() 74 * @return a status flag passed on by PKSreader 75 */ 76 int read( ); 77 78 casa::CountedPtr<Scantable> getTable() const { return table_;} 79 79 80 private: 80 casa::Int nBeam_,nIF_,nPol_,nChan_; 81 PKSreader* reader_; 81 82 PKSreader* reader_; 82 83 SDHeader* header_; 83 SDFrequencyTable* frequencies_;84 casa::CountedPtr<SDMemTable> table_;85 84 casa::String filename_; 86 casa:: uInt cursor_;87 casa:: Double timestamp_;88 casa::uInt beamOffset_, ifOffset_;85 casa::CountedPtr< Scantable > table_; 86 casa::Int nIF_, nBeam_, nPol_, nChan_; 87 casa::uInt ifOffset_, beamOffset_; 89 88 casa::Bool haveXPol_; 90 89 }; 91 90 92 }// namespace 91 } // namespace 92 93 93 #endif -
trunk/src/STMath.cpp
r804 r805 1 //#--------------------------------------------------------------------------- 2 //# SDMath.cc: A collection of single dish mathematical operations 3 //#--------------------------------------------------------------------------- 4 //# Copyright (C) 2004 5 //# ATNF 6 //# 7 //# This program is free software; you can redistribute it and/or modify it 8 //# under the terms of the GNU General Public License as published by the Free 9 //# Software Foundation; either version 2 of the License, or (at your option) 10 //# any later version. 11 //# 12 //# This program is distributed in the hope that it will be useful, but 13 //# WITHOUT ANY WARRANTY; without even the implied warranty of 14 //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 15 //# Public License for more details. 16 //# 17 //# You should have received a copy of the GNU General Public License along 18 //# with this program; if not, write to the Free Software Foundation, Inc., 19 //# 675 Massachusetts Ave, Cambridge, MA 02139, USA. 20 //# 21 //# Correspondence concerning this software should be addressed as follows: 22 //# Internet email: Malte.Marquarding@csiro.au 23 //# Postal address: Malte Marquarding, 24 //# Australia Telescope National Facility, 25 //# P.O. Box 76, 26 //# Epping, NSW, 2121, 27 //# AUSTRALIA 28 //# 29 //# $Id: 30 //#--------------------------------------------------------------------------- 31 #include <vector> 32 33 34 #include <casa/aips.h> 35 #include <casa/iostream.h> 36 #include <casa/sstream.h> 1 // 2 // C++ Implementation: STMath 3 // 4 // Description: 5 // 6 // 7 // Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006 8 // 9 // Copyright: See COPYING file that comes with this distribution 10 // 11 // 12 37 13 #include <casa/iomanip.h> 14 #include <casa/Exceptions/Error.h> 15 #include <casa/Containers/Block.h> 38 16 #include <casa/BasicSL/String.h> 39 #include <casa/Arrays/IPosition.h> 40 #include <casa/Arrays/Array.h> 41 #include <casa/Arrays/ArrayIter.h> 42 #include <casa/Arrays/VectorIter.h> 17 #include <tables/Tables/TableIter.h> 18 #include <tables/Tables/TableCopy.h> 19 #include <casa/Arrays/MaskArrLogi.h> 20 #include <casa/Arrays/MaskArrMath.h> 21 #include <casa/Arrays/ArrayLogical.h> 43 22 #include <casa/Arrays/ArrayMath.h> 44 #include <casa/Arrays/ArrayLogical.h> 45 #include <casa/Arrays/MaskedArray.h> 46 #include <casa/Arrays/MaskArrMath.h> 47 #include <casa/Arrays/MaskArrLogi.h> 48 #include <casa/Arrays/Matrix.h> 49 #include <casa/BasicMath/Math.h> 50 #include <casa/Exceptions.h> 51 #include <casa/Quanta/Quantum.h> 52 #include <casa/Quanta/Unit.h> 53 #include <casa/Quanta/MVEpoch.h> 54 #include <casa/Quanta/MVTime.h> 55 #include <casa/Utilities/Assert.h> 56 57 #include <coordinates/Coordinates/SpectralCoordinate.h> 58 #include <coordinates/Coordinates/CoordinateSystem.h> 59 #include <coordinates/Coordinates/CoordinateUtil.h> 60 #include <coordinates/Coordinates/FrequencyAligner.h> 23 #include <casa/Containers/RecordField.h> 24 #include <tables/Tables/TableRow.h> 25 #include <tables/Tables/TableVector.h> 26 #include <tables/Tables/ExprNode.h> 27 #include <tables/Tables/TableRecord.h> 28 #include <tables/Tables/ReadAsciiTable.h> 61 29 62 30 #include <lattices/Lattices/LatticeUtilities.h> 63 #include <lattices/Lattices/RebinLattice.h>64 65 #include <measures/Measures/MEpoch.h>66 #include <measures/Measures/MDirection.h>67 #include <measures/Measures/MPosition.h>68 31 69 32 #include <scimath/Mathematics/VectorKernel.h> 70 33 #include <scimath/Mathematics/Convolver.h> 71 #include <scimath/Mathematics/InterpolateArray1D.h>72 34 #include <scimath/Functionals/Polynomial.h> 73 35 74 #include <tables/Tables/Table.h>75 #include <tables/Tables/ScalarColumn.h>76 #include <tables/Tables/ArrayColumn.h>77 #include <tables/Tables/ReadAsciiTable.h>78 79 36 #include "MathUtils.h" 80 #include " SDDefs.h"37 #include "RowAccumulator.h" 81 38 #include "SDAttr.h" 82 #include "SDContainer.h" 83 #include "SDMemTable.h" 84 85 #include "SDMath.h" 86 #include "SDPol.h" 39 #include "STMath.h" 87 40 88 41 using namespace casa; 42 89 43 using namespace asap; 90 44 91 92 SDMath::SDMath() 93 { 94 } 95 96 SDMath::SDMath(const SDMath& other) 97 { 98 99 // No state 100 101 } 102 103 SDMath& SDMath::operator=(const SDMath& other) 104 { 105 if (this != &other) { 106 // No state 107 } 108 return *this; 109 } 110 111 SDMath::~SDMath() 112 {;} 113 114 115 SDMemTable* SDMath::frequencyAlignment(const SDMemTable& in, 116 const String& refTime, 117 const String& method, 118 Bool perFreqID) 119 { 120 // Get frame info from Table 121 std::vector<std::string> info = in.getCoordInfo(); 122 123 // Parse frequency system 124 String systemStr(info[1]); 125 String baseSystemStr(info[3]); 126 if (baseSystemStr==systemStr) { 127 throw(AipsError("You have not set a frequency frame different from the initial - use function set_freqframe")); 128 } 129 130 MFrequency::Types freqSystem; 131 MFrequency::getType(freqSystem, systemStr); 132 133 return frequencyAlign(in, freqSystem, refTime, method, perFreqID); 134 } 135 136 137 138 CountedPtr<SDMemTable> 139 SDMath::average(const std::vector<CountedPtr<SDMemTable> >& in, 140 const Vector<Bool>& mask, Bool scanAv, 141 const String& weightStr, Bool alignFreq) 142 // Weighted averaging of spectra from one or more Tables. 143 { 144 // Convert weight type 145 WeightType wtType = NONE; 146 convertWeightString(wtType, weightStr, True); 147 148 // Create output Table by cloning from the first table 149 SDMemTable* pTabOut = new SDMemTable(*in[0],True); 150 if (in.size() > 1) { 151 for (uInt i=1; i < in.size(); ++i) { 152 pTabOut->appendToHistoryTable(in[i]->getHistoryTable()); 153 } 154 } 155 // Setup 156 IPosition shp = in[0]->rowAsMaskedArray(0).shape(); // Must not change 157 Array<Float> arr(shp); 158 Array<Bool> barr(shp); 159 const Bool useMask = (mask.nelements() == shp(asap::ChanAxis)); 160 161 // Columns from Tables 162 ROArrayColumn<Float> tSysCol; 163 ROScalarColumn<Double> mjdCol; 164 ROScalarColumn<String> srcNameCol; 165 ROScalarColumn<Double> intCol; 166 ROArrayColumn<uInt> fqIDCol; 45 STMath::STMath(bool insitu) : 46 insitu_(insitu) 47 { 48 } 49 50 51 STMath::~STMath() 52 { 53 } 54 55 CountedPtr<Scantable> 56 STMath::average( const std::vector<CountedPtr<Scantable> >& in, 57 const Vector<Bool>& mask, 58 const std::string& weight, 59 const std::string& avmode, 60 bool alignfreq) 61 { 62 if ( avmode == "SCAN" && in.size() != 1 ) 63 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables")); 64 WeightType wtype = stringToWeight(weight); 65 // output 66 // clone as this is non insitu 67 bool insitu = insitu_; 68 setInsitu(false); 69 CountedPtr< Scantable > out = getScantable(in[0], true); 70 setInsitu(insitu); 71 72 Table& tout = out->table(); 73 74 /// @todo check if all scantables are conformant 75 76 ArrayColumn<Float> specColOut(tout,"SPECTRA"); 77 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA"); 78 ArrayColumn<Float> tsysColOut(tout,"TSYS"); 79 ScalarColumn<Double> mjdColOut(tout,"TIME"); 80 ScalarColumn<Double> intColOut(tout,"INTERVAL"); 81 82 // set up the output table rows. These are based on the structure of the 83 // FIRST scantabel in the vector 84 const Table& baset = in[0]->table(); 85 86 Block<String> cols(3); 87 cols[0] = String("BEAMNO"); 88 cols[1] = String("IFNO"); 89 cols[2] = String("POLNO"); 90 if ( avmode == "SOURCE" ) { 91 cols.resize(4); 92 cols[3] = String("SRCNAME"); 93 } 94 if ( avmode == "SCAN" && in.size() == 1) { 95 cols.resize(4); 96 cols[3] = String("SCANNO"); 97 } 98 uInt outrowCount = 0; 99 TableIterator iter(baset, cols); 100 while (!iter.pastEnd()) { 101 Table subt = iter.table(); 102 // copy the first row of this selection into the new table 103 tout.addRow(); 104 TableCopy::copyRows(tout, subt, outrowCount, 0, 1); 105 ++outrowCount; 106 ++iter; 107 } 108 RowAccumulator acc(wtype); 109 acc.setUserMask(mask); 110 ROTableRow row(tout); 111 ROArrayColumn<Float> specCol, tsysCol; 112 ROArrayColumn<uChar> flagCol; 113 ROScalarColumn<Double> mjdCol, intCol; 167 114 ROScalarColumn<Int> scanIDCol; 168 115 169 // Create accumulation MaskedArray. We accumulate for each 170 // channel,if,pol,beam Note that the mask of the accumulation array 171 // will ALWAYS remain ALL True. The MA is only used so that when 172 // data which is masked Bad is added to it, that data does not 173 // contribute. 174 175 Array<Float> zero(shp); 176 zero=0.0; 177 Array<Bool> good(shp); 178 good = True; 179 MaskedArray<Float> sum(zero,good); 180 181 // Counter arrays 182 Array<Float> nPts(shp); // Number of points 183 nPts = 0.0; 184 Array<Float> nInc(shp); // Increment 185 nInc = 1.0; 186 187 // Create accumulation Array for variance. We accumulate for each 188 // if,pol,beam, but average over channel. So we need a shape with 189 // one less axis dropping channels. 190 const uInt nAxesSub = shp.nelements() - 1; 191 IPosition shp2(nAxesSub); 192 for (uInt i=0,j=0; i<(nAxesSub+1); i++) { 193 if (i!=asap::ChanAxis) { 194 shp2(j) = shp(i); 195 j++; 196 } 197 } 198 Array<Float> sumSq(shp2); 199 sumSq = 0.0; 200 IPosition pos2(nAxesSub,0); // For indexing 201 202 // Time-related accumulators 203 Double time; 204 Double timeSum = 0.0; 205 Double intSum = 0.0; 206 Double interval = 0.0; 207 208 // To get the right shape for the Tsys accumulator we need to access 209 // a column from the first table. The shape of this array must not 210 // change. Note however that since the TSysSqSum array is used in a 211 // normalization process, and that I ignore the channel axis 212 // replication of values for now, it loses a dimension 213 214 Array<Float> tSysSum, tSysSqSum; 215 { 216 const Table& tabIn = in[0]->table(); 217 tSysCol.attach(tabIn,"TSYS"); 218 tSysSum.resize(tSysCol.shape(0)); 219 tSysSqSum.resize(shp2); 220 } 221 tSysSum = 0.0; 222 tSysSqSum = 0.0; 223 Array<Float> tSys; 224 225 // Scan and row tracking 226 Int oldScanID = 0; 227 Int outScanID = 0; 228 Int scanID = 0; 229 Int rowStart = 0; 230 Int nAccum = 0; 231 Int tableStart = 0; 232 233 // Source and FreqID 234 String sourceName, oldSourceName, sourceNameStart; 235 Vector<uInt> freqID, freqIDStart, oldFreqID; 236 237 // Loop over tables 238 Float fac = 1.0; 239 const uInt nTables = in.size(); 240 for (uInt iTab=0; iTab<nTables; iTab++) { 241 242 // Should check that the frequency tables don't change if doing 243 // FreqAlignment 244 245 // Attach columns to Table 246 const Table& tabIn = in[iTab]->table(); 247 tSysCol.attach(tabIn, "TSYS"); 248 mjdCol.attach(tabIn, "TIME"); 249 srcNameCol.attach(tabIn, "SRCNAME"); 250 intCol.attach(tabIn, "INTERVAL"); 251 fqIDCol.attach(tabIn, "FREQID"); 252 scanIDCol.attach(tabIn, "SCANID"); 253 254 // Loop over rows in Table 255 const uInt nRows = in[iTab]->nRow(); 256 for (uInt iRow=0; iRow<nRows; iRow++) { 257 // Check conformance 258 IPosition shp2 = in[iTab]->rowAsMaskedArray(iRow).shape(); 259 if (!shp.isEqual(shp2)) { 260 delete pTabOut; 261 throw (AipsError("Shapes for all rows must be the same")); 116 for (uInt i=0; i < tout.nrow(); ++i) { 117 for ( int j=0; j < in.size(); ++j ) { 118 const Table& tin = in[j]->table(); 119 const TableRecord& rec = row.get(i); 120 ROScalarColumn<Double> tmp(tin, "TIME"); 121 Double td;tmp.get(0,td); 122 Table basesubt = tin(tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO")) 123 && tin.col("IFNO") == Int(rec.asuInt("IFNO")) 124 && tin.col("POLNO") == Int(rec.asuInt("POLNO")) ); 125 Table subt; 126 if ( avmode == "SOURCE") { 127 subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") ); 128 } else if (avmode == "SCAN") { 129 subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) ); 130 } else { 131 subt = basesubt; 132 } 133 specCol.attach(subt,"SPECTRA"); 134 flagCol.attach(subt,"FLAGTRA"); 135 tsysCol.attach(subt,"TSYS"); 136 intCol.attach(subt,"INTERVAL"); 137 mjdCol.attach(subt,"TIME"); 138 Vector<Float> spec,tsys; 139 Vector<uChar> flag; 140 Double inter,time; 141 for (uInt k = 0; k < subt.nrow(); ++k ) { 142 flagCol.get(k, flag); 143 Vector<Bool> bflag(flag.shape()); 144 convertArray(bflag, flag); 145 if ( allEQ(bflag, True) ) { 146 continue;//don't accumulate 262 147 } 263 264 // If we are not doing scan averages, make checks for source 265 // and frequency setup and warn if averaging across them 266 scanIDCol.getScalar(iRow, scanID); 267 268 // Get quantities from columns 269 srcNameCol.getScalar(iRow, sourceName); 270 mjdCol.get(iRow, time); 271 tSysCol.get(iRow, tSys); 272 intCol.get(iRow, interval); 273 fqIDCol.get(iRow, freqID); 274 275 // Initialize first source and freqID 276 if (iRow==0 && iTab==0) { 277 sourceNameStart = sourceName; 278 freqIDStart = freqID; 279 } 280 281 // If we are doing scan averages, see if we are at the end of 282 // an accumulation period (scan). We must check soutce names 283 // too, since we might have two tables with one scan each but 284 // different source names; we shouldn't average different 285 // sources together 286 if (scanAv && ( (scanID != oldScanID) || 287 (iRow==0 && iTab>0 && sourceName!=oldSourceName))) { 288 289 // Normalize data in 'sum' accumulation array according to 290 // weighting scheme 291 normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, 292 asap::ChanAxis, nAxesSub); 293 294 // Get ScanContainer for the first row of this averaged Scan 295 SDContainer scOut = in[iTab]->getSDContainer(rowStart); 296 297 // Fill scan container. The source and freqID come from the 298 // first row of the first table that went into this average 299 // ( should be the same for all rows in the scan average) 300 301 Float nR(nAccum); 302 fillSDC(scOut, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID, 303 timeSum/nR, intSum, sourceNameStart, freqIDStart); 304 305 // Write container out to Table 306 pTabOut->putSDContainer(scOut); 307 308 // Reset accumulators 309 sum = 0.0; 310 sumSq = 0.0; 311 nAccum = 0; 312 313 tSysSum =0.0; 314 tSysSqSum =0.0; 315 timeSum = 0.0; 316 intSum = 0.0; 317 nPts = 0.0; 318 319 // Increment 320 rowStart = iRow; // First row for next accumulation 321 tableStart = iTab; // First table for next accumulation 322 sourceNameStart = sourceName; // First source name for next 323 // accumulation 324 freqIDStart = freqID; // First FreqID for next accumulation 325 326 oldScanID = scanID; 327 outScanID += 1; // Scan ID for next 328 // accumulation period 329 } 330 331 // Accumulate 332 accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts, 333 tSysSum, tSysSqSum, tSys, 334 nInc, mask, time, interval, in, iTab, iRow, asap::ChanAxis, 335 nAxesSub, useMask, wtType); 336 oldSourceName = sourceName; 337 oldFreqID = freqID; 338 } 339 } 340 341 // OK at this point we have accumulation data which is either 342 // - accumulated from all tables into one row 343 // or 344 // - accumulated from the last scan average 345 // 346 // Normalize data in 'sum' accumulation array according to weighting 347 // scheme 348 349 normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, 350 asap::ChanAxis, nAxesSub); 351 352 // Create and fill container. The container we clone will be from 353 // the last Table and the first row that went into the current 354 // accumulation. It probably doesn't matter that much really... 355 Float nR(nAccum); 356 SDContainer scOut = in[tableStart]->getSDContainer(rowStart); 357 fillSDC(scOut, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID, 358 timeSum/nR, intSum, sourceNameStart, freqIDStart); 359 pTabOut->putSDContainer(scOut); 360 pTabOut->resetCursor(); 361 362 return CountedPtr<SDMemTable>(pTabOut); 363 } 364 365 366 367 CountedPtr<SDMemTable> 368 SDMath::binaryOperate(const CountedPtr<SDMemTable>& left, 369 const CountedPtr<SDMemTable>& right, 370 const String& op, Bool preserve, Bool doTSys) 371 { 372 373 // Check operator 374 String op2(op); 375 op2.upcase(); 376 uInt what = 0; 377 if (op2=="ADD") { 378 what = 0; 379 } else if (op2=="SUB") { 380 what = 1; 381 } else if (op2=="MUL") { 382 what = 2; 383 } else if (op2=="DIV") { 384 what = 3; 385 } else if (op2=="QUOTIENT") { 386 what = 4; 387 doTSys = True; 388 } else { 389 throw( AipsError("Unrecognized operation")); 390 } 391 392 // Check rows 393 const uInt nRowLeft = left->nRow(); 394 const uInt nRowRight = right->nRow(); 395 Bool ok = (nRowRight==1 && nRowLeft>0) || 396 (nRowRight>=1 && nRowLeft==nRowRight); 397 if (!ok) { 398 throw (AipsError("The right Scan Table can have one row or the same number of rows as the left Scan Table")); 399 } 400 401 // Input Tables 402 const Table& tLeft = left->table(); 403 const Table& tRight = right->table(); 404 405 // TSys columns 406 ROArrayColumn<Float> tSysLeftCol, tSysRightCol; 407 if (doTSys) { 408 tSysLeftCol.attach(tLeft, "TSYS"); 409 tSysRightCol.attach(tRight, "TSYS"); 410 } 411 412 // First row for right 413 Array<Float> tSysLeftArr, tSysRightArr; 414 if (doTSys) tSysRightCol.get(0, tSysRightArr); 415 MaskedArray<Float>* pMRight = 416 new MaskedArray<Float>(right->rowAsMaskedArray(0)); 417 418 IPosition shpRight = pMRight->shape(); 419 420 // Output Table cloned from left 421 SDMemTable* pTabOut = new SDMemTable(*left, True); 422 pTabOut->appendToHistoryTable(right->getHistoryTable()); 423 424 // Loop over rows 425 for (uInt i=0; i<nRowLeft; i++) { 426 427 // Get data 428 MaskedArray<Float> mLeft(left->rowAsMaskedArray(i)); 429 IPosition shpLeft = mLeft.shape(); 430 if (doTSys) tSysLeftCol.get(i, tSysLeftArr); 431 432 if (nRowRight>1) { 433 delete pMRight; 434 pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(i)); 435 shpRight = pMRight->shape(); 436 if (doTSys) tSysRightCol.get(i, tSysRightArr); 437 } 438 439 if (!shpRight.isEqual(shpLeft)) { 440 delete pTabOut; 441 delete pMRight; 442 throw(AipsError("left and right scan tables are not conformant")); 443 } 444 if (doTSys) { 445 if (!tSysRightArr.shape().isEqual(tSysRightArr.shape())) { 446 delete pTabOut; 447 delete pMRight; 448 throw(AipsError("left and right Tsys data are not conformant")); 148 specCol.get(k, spec); 149 tsysCol.get(k, tsys); 150 intCol.get(k, inter); 151 mjdCol.get(k, time); 152 // spectrum has to be added last to enable weighting by the other values 153 acc.add(spec, !bflag, tsys, inter, time); 449 154 } 450 if (!shpRight.isEqual(tSysRightArr.shape())) { 451 delete pTabOut; 452 delete pMRight; 453 throw(AipsError("left and right scan tables are not conformant")); 155 } 156 //write out 157 specColOut.put(i, acc.getSpectrum()); 158 const Vector<Bool>& msk = acc.getMask(); 159 Vector<uChar> flg(msk.shape()); 160 convertArray(flg, !msk); 161 flagColOut.put(i, flg); 162 tsysColOut.put(i, acc.getTsys()); 163 intColOut.put(i, acc.getInterval()); 164 mjdColOut.put(i, acc.getTime()); 165 acc.reset(); 166 } 167 return out; 168 } 169 170 CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in, 171 bool droprows) 172 { 173 if (insitu_) return in; 174 else { 175 // clone 176 Scantable* tabp = new Scantable(*in, Bool(droprows)); 177 return CountedPtr<Scantable>(tabp); 178 } 179 } 180 181 CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in, 182 float val, 183 const std::string& mode, 184 bool tsys ) 185 { 186 // modes are "ADD" and "MUL" 187 CountedPtr< Scantable > out = getScantable(in, false); 188 Table& tab = out->table(); 189 ArrayColumn<Float> specCol(tab,"SPECTRA"); 190 ArrayColumn<Float> tsysCol(tab,"TSYS"); 191 for (uInt i=0; i<tab.nrow(); ++i) { 192 Vector<Float> spec; 193 Vector<Float> ts; 194 specCol.get(i, spec); 195 tsysCol.get(i, ts); 196 if (mode == "MUL") { 197 spec *= val; 198 specCol.put(i, spec); 199 if ( tsys ) { 200 ts *= val; 201 tsysCol.put(i, ts); 454 202 } 455 } 456 457 // Make container 458 SDContainer sc = left->getSDContainer(i); 459 460 // Operate on data and TSys 461 if (what==0) { 462 MaskedArray<Float> tmp = mLeft + *pMRight; 463 putDataInSDC(sc, tmp.getArray(), tmp.getMask()); 464 if (doTSys) sc.putTsys(tSysLeftArr+tSysRightArr); 465 } else if (what==1) { 466 MaskedArray<Float> tmp = mLeft - *pMRight; 467 putDataInSDC(sc, tmp.getArray(), tmp.getMask()); 468 if (doTSys) sc.putTsys(tSysLeftArr-tSysRightArr); 469 } else if (what==2) { 470 MaskedArray<Float> tmp = mLeft * *pMRight; 471 putDataInSDC(sc, tmp.getArray(), tmp.getMask()); 472 if (doTSys) sc.putTsys(tSysLeftArr*tSysRightArr); 473 } else if (what==3) { 474 MaskedArray<Float> tmp = mLeft / *pMRight; 475 putDataInSDC(sc, tmp.getArray(), tmp.getMask()); 476 if (doTSys) sc.putTsys(tSysLeftArr/tSysRightArr); 477 } else if (what==4) { 478 if (preserve) { 479 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) - 480 tSysRightArr; 481 putDataInSDC(sc, tmp.getArray(), tmp.getMask()); 482 } else { 483 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) - 484 tSysLeftArr; 485 putDataInSDC(sc, tmp.getArray(), tmp.getMask()); 486 } 487 sc.putTsys(tSysRightArr); 488 } 489 490 // Put new row in output Table 491 pTabOut->putSDContainer(sc); 492 } 493 if (pMRight) delete pMRight; 494 pTabOut->resetCursor(); 495 496 return CountedPtr<SDMemTable>(pTabOut); 497 } 498 499 500 std::vector<float> SDMath::statistic(const CountedPtr<SDMemTable>& in, 501 const Vector<Bool>& mask, 502 const String& which, Int row) const 503 // 504 // Perhaps iteration over pol/beam/if should be in here 505 // and inside the nrow iteration ? 506 // 507 { 508 const uInt nRow = in->nRow(); 509 510 // Specify cursor location 511 512 IPosition start, end; 513 Bool doAll = False; 514 setCursorSlice(start, end, doAll, *in); 515 516 // Loop over rows 517 518 const uInt nEl = mask.nelements(); 519 uInt iStart = 0; 520 uInt iEnd = in->nRow()-1; 521 // 522 if (row>=0) { 523 iStart = row; 524 iEnd = row; 525 } 526 // 527 std::vector<float> result(iEnd-iStart+1); 528 for (uInt ii=iStart; ii <= iEnd; ++ii) { 529 530 // Get row and deconstruct 531 532 MaskedArray<Float> dataIn = (in->rowAsMaskedArray(ii))(start,end); 533 Array<Float> v = dataIn.getArray().nonDegenerate(); 534 Array<Bool> m = dataIn.getMask().nonDegenerate(); 535 536 // Access desired piece of data 537 538 // Array<Float> v((arr(start,end)).nonDegenerate()); 539 // Array<Bool> m((barr(start,end)).nonDegenerate()); 540 541 // Apply OTF mask 542 543 MaskedArray<Float> tmp; 544 if (m.nelements()==nEl) { 545 tmp.setData(v,m&&mask); 546 } else { 547 tmp.setData(v,m); 548 } 549 550 // Get statistic 551 552 result[ii-iStart] = mathutil::statistics(which, tmp); 553 } 554 // 555 return result; 556 } 557 558 559 SDMemTable* SDMath::bin(const SDMemTable& in, Int width) 560 { 561 SDHeader sh = in.getSDHeader(); 562 SDMemTable* pTabOut = new SDMemTable(in, True); 563 564 // Bin up SpectralCoordinates 565 566 IPosition factors(1); 567 factors(0) = width; 568 for (uInt j=0; j<in.nCoordinates(); ++j) { 569 CoordinateSystem cSys; 570 cSys.addCoordinate(in.getSpectralCoordinate(j)); 571 CoordinateSystem cSysBin = 572 CoordinateUtil::makeBinnedCoordinateSystem(factors, cSys, False); 573 // 574 SpectralCoordinate sCBin = cSysBin.spectralCoordinate(0); 575 pTabOut->setCoordinate(sCBin, j); 576 } 577 578 // Use RebinLattice to find shape 579 580 IPosition shapeIn(1,sh.nchan); 581 IPosition shapeOut = RebinLattice<Float>::rebinShape(shapeIn, factors); 582 sh.nchan = shapeOut(0); 583 pTabOut->putSDHeader(sh); 584 585 // Loop over rows and bin along channel axis 586 587 for (uInt i=0; i < in.nRow(); ++i) { 588 SDContainer sc = in.getSDContainer(i); 589 // 590 Array<Float> tSys(sc.getTsys()); // Get it out before sc changes shape 591 592 // Bin up spectrum 593 594 MaskedArray<Float> marr(in.rowAsMaskedArray(i)); 595 MaskedArray<Float> marrout; 596 LatticeUtilities::bin(marrout, marr, asap::ChanAxis, width); 597 598 // Put back the binned data and flags 599 600 IPosition ip2 = marrout.shape(); 601 sc.resize(ip2); 602 // 603 putDataInSDC(sc, marrout.getArray(), marrout.getMask()); 604 605 // Bin up Tsys. 606 607 Array<Bool> allGood(tSys.shape(),True); 608 MaskedArray<Float> tSysIn(tSys, allGood, True); 609 // 610 MaskedArray<Float> tSysOut; 611 LatticeUtilities::bin(tSysOut, tSysIn, asap::ChanAxis, width); 612 sc.putTsys(tSysOut.getArray()); 613 // 614 pTabOut->putSDContainer(sc); 615 } 616 return pTabOut; 617 } 618 619 SDMemTable* SDMath::resample(const SDMemTable& in, const String& methodStr, 620 Float width) 203 } else if ( mode == "ADD" ) { 204 spec += val; 205 specCol.put(i, spec); 206 if ( tsys ) { 207 ts += val; 208 tsysCol.put(i, ts); 209 } 210 } 211 } 212 return out; 213 } 214 215 MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s, 216 const Vector<uChar>& f) 217 { 218 Vector<Bool> mask; 219 mask.resize(f.shape()); 220 convertArray(mask, f); 221 return MaskedArray<Float>(s,!mask); 222 } 223 224 Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma) 225 { 226 const Vector<Bool>& m = ma.getMask(); 227 Vector<uChar> flags(m.shape()); 228 convertArray(flags, !m); 229 return flags; 230 } 231 232 CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable >& in, 233 const std::string & mode, 234 bool preserve ) 235 { 236 /// @todo make other modes available 237 /// modes should be "nearest", "pair" 238 // make this operation non insitu 239 const Table& tin = in->table(); 240 Table ons = tin(tin.col("SRCTYPE") == Int(0)); 241 Table offs = tin(tin.col("SRCTYPE") == Int(1)); 242 if ( offs.nrow() == 0 ) 243 throw(AipsError("No 'off' scans present.")); 244 // put all "on" scans into output table 245 246 bool insitu = insitu_; 247 setInsitu(false); 248 CountedPtr< Scantable > out = getScantable(in, true); 249 setInsitu(insitu); 250 Table& tout = out->table(); 251 252 TableCopy::copyRows(tout, ons); 253 TableRow row(tout); 254 ROScalarColumn<Double> offtimeCol(offs, "TIME"); 255 256 ArrayColumn<Float> outspecCol(tout, "SPECTRA"); 257 ROArrayColumn<Float> outtsysCol(tout, "TSYS"); 258 ArrayColumn<uChar> outflagCol(tout, "FLAGTRA"); 259 260 for (uInt i=0; i < tout.nrow(); ++i) { 261 const TableRecord& rec = row.get(i); 262 Double ontime = rec.asDouble("TIME"); 263 ROScalarColumn<Double> offtimeCol(offs, "TIME"); 264 Double mindeltat = min(abs(offtimeCol.getColumn() - ontime)); 265 Table sel = offs( abs(offs.col("TIME")-ontime) <= mindeltat 266 && offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO")) 267 && offs.col("IFNO") == Int(rec.asuInt("IFNO")) 268 && offs.col("POLNO") == Int(rec.asuInt("POLNO")) ); 269 270 TableRow offrow(sel); 271 const TableRecord& offrec = offrow.get(0);//should only be one row 272 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA"); 273 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS"); 274 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA"); 275 /// @fixme this assumes tsys is a scalar not vector 276 Float tsysoffscalar = (*tsysoff)(IPosition(1,0)); 277 Vector<Float> specon, tsyson; 278 outtsysCol.get(i, tsyson); 279 outspecCol.get(i, specon); 280 Vector<uChar> flagon; 281 outflagCol.get(i, flagon); 282 MaskedArray<Float> mon = maskedArray(specon, flagon); 283 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff); 284 MaskedArray<Float> quot = (tsysoffscalar * mon / moff); 285 if (preserve) { 286 quot -= tsysoffscalar; 287 } else { 288 quot -= tsyson[0]; 289 } 290 outspecCol.put(i, quot.getArray()); 291 outflagCol.put(i, flagsFromMA(quot)); 292 } 293 return out; 294 } 295 296 CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in ) 297 { 298 // make copy or reference 299 CountedPtr< Scantable > out = getScantable(in, false); 300 Table& tout = out->table(); 301 Block<String> cols(3); 302 cols[0] = String("SCANNO"); 303 cols[1] = String("BEAMNO"); 304 cols[2] = String("POLNO"); 305 TableIterator iter(tout, cols); 306 while (!iter.pastEnd()) { 307 Table subt = iter.table(); 308 // this should leave us with two rows for the two IFs....if not ignore 309 if (subt.nrow() != 2 ) { 310 continue; 311 } 312 ArrayColumn<Float> specCol(tout, "SPECTRA"); 313 ArrayColumn<Float> tsysCol(tout, "TSYS"); 314 ArrayColumn<uChar> flagCol(tout, "FLAGTRA"); 315 Vector<Float> onspec,offspec, ontsys, offtsys; 316 Vector<uChar> onflag, offflag; 317 tsysCol.get(0, ontsys); tsysCol.get(1, offtsys); 318 specCol.get(0, onspec); specCol.get(1, offspec); 319 flagCol.get(0, onflag); flagCol.get(1, offflag); 320 MaskedArray<Float> on = maskedArray(onspec, onflag); 321 MaskedArray<Float> off = maskedArray(offspec, offflag); 322 MaskedArray<Float> oncopy = on.copy(); 323 324 on /= off; on -= 1.0f; 325 on *= ontsys[0]; 326 off /= oncopy; off -= 1.0f; 327 off *= offtsys[0]; 328 specCol.put(0, on.getArray()); 329 const Vector<Bool>& m0 = on.getMask(); 330 Vector<uChar> flags0(m0.shape()); 331 convertArray(flags0, !m0); 332 flagCol.put(0, flags0); 333 334 specCol.put(1, off.getArray()); 335 const Vector<Bool>& m1 = off.getMask(); 336 Vector<uChar> flags1(m1.shape()); 337 convertArray(flags1, !m1); 338 flagCol.put(1, flags1); 339 340 } 341 342 return out; 343 } 344 345 std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in, 346 const std::vector< bool > & mask, 347 const std::string& which ) 348 { 349 350 Vector<Bool> m(mask); 351 const Table& tab = in->table(); 352 ROArrayColumn<Float> specCol(tab, "SPECTRA"); 353 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA"); 354 std::vector<float> out; 355 for (uInt i=0; i < tab.nrow(); ++i ) { 356 Vector<Float> spec; specCol.get(i, spec); 357 MaskedArray<Float> ma = maskedArray(spec, flagCol(i)); 358 float outstat; 359 if ( spec.nelements() == m.nelements() ) { 360 outstat = mathutil::statistics(which, ma(m)); 361 } else { 362 outstat = mathutil::statistics(which, ma); 363 } 364 out.push_back(outstat); 365 } 366 return out; 367 } 368 369 CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in, 370 int width ) 371 { 372 if ( !in->selection().empty() ) throw(AipsError("Can't bin subset of the data.")); 373 CountedPtr< Scantable > out = getScantable(in, false); 374 Table& tout = out->table(); 375 out->frequencies().rescale(width, "BIN"); 376 ArrayColumn<Float> specCol(tout, "SPECTRA"); 377 ArrayColumn<uChar> flagCol(tout, "FLAGTRA"); 378 for (uInt i=0; i < tout.nrow(); ++i ) { 379 MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i)); 380 MaskedArray<Float> maout; 381 LatticeUtilities::bin(maout, main, 0, Int(width)); 382 /// @todo implement channel based tsys binning 383 specCol.put(i, maout.getArray()); 384 flagCol.put(i, flagsFromMA(maout)); 385 // take only the first binned spectrum's length for the deprecated 386 // global header item nChan 387 if (i==0) tout.rwKeywordSet().define(String("nChan"), 388 Int(maout.getArray().nelements())); 389 } 390 return out; 391 } 392 393 CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in, 394 const std::string& method, 395 float width ) 621 396 // 622 397 // Should add the possibility of width being specified in km/s. This means … … 627 402 // 628 403 { 629 Bool doVel = False;630 if (doVel) {631 for (uInt j=0; j<in.nCoordinates(); ++j) {632 SpectralCoordinate sC = in.getSpectralCoordinate(j);633 }634 }635 636 // Interpolation method637 638 404 InterpolateArray1D<Double,Float>::InterpolationMethod interp; 639 convertInterpString(interp, methodStr); 640 Int interpMethod(interp); 641 642 // Make output table 643 644 SDMemTable* pTabOut = new SDMemTable(in, True); 405 Int interpMethod(stringToIMethod(method)); 406 407 CountedPtr< Scantable > out = getScantable(in, false); 408 Table& tout = out->table(); 645 409 646 410 // Resample SpectralCoordinates (one per freqID) 647 648 const uInt nCoord = in.nCoordinates(); 649 Vector<Float> offset(1,0.0); 650 Vector<Float> factors(1,1.0/width); 651 Vector<Int> newShape; 652 for (uInt j=0; j<in.nCoordinates(); ++j) { 653 CoordinateSystem cSys; 654 cSys.addCoordinate(in.getSpectralCoordinate(j)); 655 CoordinateSystem cSys2 = cSys.subImage(offset, factors, newShape); 656 SpectralCoordinate sC = cSys2.spectralCoordinate(0); 657 // 658 pTabOut->setCoordinate(sC, j); 659 } 660 661 // Get header 662 663 SDHeader sh = in.getSDHeader(); 664 665 // Generate resampling vectors 666 667 const uInt nChanIn = sh.nchan; 668 Vector<Float> xIn(nChanIn); 669 indgen(xIn); 670 // 671 Int fac = Int(nChanIn/width); 672 Vector<Float> xOut(fac+10); // 10 to be safe - resize later 673 uInt i = 0; 674 Float x = 0.0; 675 Bool more = True; 676 while (more) { 677 xOut(i) = x; 678 // 679 i++; 680 x += width; 681 if (x>nChanIn-1) more = False; 682 } 683 const uInt nChanOut = i; 684 xOut.resize(nChanOut,True); 685 // 686 IPosition shapeIn(in.rowAsMaskedArray(0).shape()); 687 sh.nchan = nChanOut; 688 pTabOut->putSDHeader(sh); 689 690 // Loop over rows and resample along channel axis 691 692 Array<Float> valuesOut; 693 Array<Bool> maskOut; 694 Array<Float> tSysOut; 695 Array<Bool> tSysMaskIn(shapeIn,True); 696 Array<Bool> tSysMaskOut; 697 for (uInt i=0; i < in.nRow(); ++i) { 698 699 // Get container 700 701 SDContainer sc = in.getSDContainer(i); 702 703 // Get data and Tsys 704 705 const Array<Float>& tSysIn = sc.getTsys(); 706 const MaskedArray<Float>& dataIn(in.rowAsMaskedArray(i)); 707 Array<Float> valuesIn = dataIn.getArray(); 708 Array<Bool> maskIn = dataIn.getMask(); 709 710 // Interpolate data 711 712 InterpolateArray1D<Float,Float>::interpolate(valuesOut, maskOut, xOut, 713 xIn, valuesIn, maskIn, 714 interpMethod, True, True); 715 sc.resize(valuesOut.shape()); 716 putDataInSDC(sc, valuesOut, maskOut); 717 718 // Interpolate TSys 719 720 InterpolateArray1D<Float,Float>::interpolate(tSysOut, tSysMaskOut, xOut, 721 xIn, tSysIn, tSysMaskIn, 722 interpMethod, True, True); 723 sc.putTsys(tSysOut); 724 725 // Put container in output 726 727 pTabOut->putSDContainer(sc); 728 } 729 // 730 return pTabOut; 731 } 732 733 SDMemTable* SDMath::unaryOperate(const SDMemTable& in, Float val, Bool doAll, 734 uInt what, Bool doTSys) 735 // 736 // what = 0 Multiply 737 // 1 Add 738 { 739 SDMemTable* pOut = new SDMemTable(in,False); 740 const Table& tOut = pOut->table(); 741 ArrayColumn<Float> specCol(tOut,"SPECTRA"); 742 ArrayColumn<Float> tSysCol(tOut,"TSYS"); 743 Array<Float> tSysArr; 744 745 // Get data slice bounds 746 747 IPosition start, end; 748 setCursorSlice (start, end, doAll, in); 749 // 750 for (uInt i=0; i<tOut.nrow(); i++) { 751 752 // Modify data 753 754 MaskedArray<Float> dataIn(pOut->rowAsMaskedArray(i)); 755 MaskedArray<Float> dataIn2 = dataIn(start,end); // Reference 756 if (what==0) { 757 dataIn2 *= val; 758 } else if (what==1) { 759 dataIn2 += val; 760 } 761 specCol.put(i, dataIn.getArray()); 762 763 // Modify Tsys 764 765 if (doTSys) { 766 tSysCol.get(i, tSysArr); 767 Array<Float> tSysArr2 = tSysArr(start,end); // Reference 768 if (what==0) { 769 tSysArr2 *= val; 770 } else if (what==1) { 771 tSysArr2 += val; 772 } 773 tSysCol.put(i, tSysArr); 774 } 775 } 776 // 777 return pOut; 778 } 779 780 SDMemTable* SDMath::averagePol(const SDMemTable& in, const Vector<Bool>& mask, 781 const String& weightStr) 782 // 783 // Average all polarizations together, weighted by variance 784 // 785 { 786 WeightType wtType = NONE; 787 convertWeightString(wtType, weightStr, True); 788 789 // Create output Table and reshape number of polarizations 790 791 Bool clear=True; 792 SDMemTable* pTabOut = new SDMemTable(in, clear); 793 SDHeader header = pTabOut->getSDHeader(); 794 header.npol = 1; 795 pTabOut->putSDHeader(header); 796 // 797 const Table& tabIn = in.table(); 798 799 // Shape of input and output data 800 801 const IPosition& shapeIn = in.rowAsMaskedArray(0).shape(); 802 IPosition shapeOut(shapeIn); 803 shapeOut(asap::PolAxis) = 1; // Average all polarizations 804 if (shapeIn(asap::PolAxis)==1) { 805 delete pTabOut; 806 throw(AipsError("The input has only one polarisation")); 807 } 808 // 809 const uInt nRows = in.nRow(); 810 const uInt nChan = shapeIn(asap::ChanAxis); 811 AlwaysAssert(asap::nAxes==4,AipsError); 812 const IPosition vecShapeOut(4,1,1,1,nChan); // A multi-dim form of a Vector shape 813 IPosition start(4), end(4); 814 815 // Output arrays 816 817 Array<Float> outData(shapeOut, 0.0); 818 Array<Bool> outMask(shapeOut, True); 819 const IPosition axes(2, asap::PolAxis, asap::ChanAxis); // pol-channel plane 820 821 // Attach Tsys column if needed 822 823 ROArrayColumn<Float> tSysCol; 824 Array<Float> tSys; 825 if (wtType==TSYS) { 826 tSysCol.attach(tabIn,"TSYS"); 827 } 828 // 829 const Bool useMask = (mask.nelements() == shapeIn(asap::ChanAxis)); 830 831 // Loop over rows 832 833 for (uInt iRow=0; iRow<nRows; iRow++) { 834 835 // Get data for this row 836 837 MaskedArray<Float> marr(in.rowAsMaskedArray(iRow)); 838 Array<Float>& arr = marr.getRWArray(); 839 const Array<Bool>& barr = marr.getMask(); 840 841 // Get Tsys 842 843 if (wtType==TSYS) { 844 tSysCol.get(iRow,tSys); 845 } 846 847 // Make iterators to iterate by pol-channel planes 848 // The tSys array is empty unless wtType=TSYS so only 849 // access the iterator is that is the case 850 851 ReadOnlyArrayIterator<Float> itDataPlane(arr, axes); 852 ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes); 853 ReadOnlyArrayIterator<Float>* pItTsysPlane = 0; 854 if (wtType==TSYS) 855 pItTsysPlane = new ReadOnlyArrayIterator<Float>(tSys, axes); 856 857 // Accumulations 858 859 Float fac = 1.0; 860 Vector<Float> vecSum(nChan,0.0); 861 862 // Iterate through data by pol-channel planes 863 864 while (!itDataPlane.pastEnd()) { 865 866 // Iterate through plane by polarization and accumulate Vectors 867 868 Vector<Float> t1(nChan); t1 = 0.0; 869 Vector<Bool> t2(nChan); t2 = True; 870 Float tSys = 0.0; 871 MaskedArray<Float> vecSum(t1,t2); 872 Float norm = 0.0; 873 { 874 ReadOnlyVectorIterator<Float> itDataVec(itDataPlane.array(), 1); 875 ReadOnlyVectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1); 876 // 877 ReadOnlyVectorIterator<Float>* pItTsysVec = 0; 878 if (wtType==TSYS) { 879 pItTsysVec = 880 new ReadOnlyVectorIterator<Float>(pItTsysPlane->array(), 1); 881 } 882 // 883 while (!itDataVec.pastEnd()) { 884 885 // Create MA of data & mask (optionally including OTF mask) and get variance for this spectrum 886 887 if (useMask) { 888 const MaskedArray<Float> spec(itDataVec.vector(), 889 mask&&itMaskVec.vector()); 890 if (wtType==VAR) { 891 fac = 1.0 / variance(spec); 892 } else if (wtType==TSYS) { 893 tSys = pItTsysVec->vector()[0]; // Drop pseudo channel dependency 894 fac = 1.0 / tSys / tSys; 895 } 896 } else { 897 const MaskedArray<Float> spec(itDataVec.vector(), 898 itMaskVec.vector()); 899 if (wtType==VAR) { 900 fac = 1.0 / variance(spec); 901 } else if (wtType==TSYS) { 902 tSys = pItTsysVec->vector()[0]; // Drop pseudo channel dependency 903 fac = 1.0 / tSys / tSys; 904 } 905 } 906 907 // Normalize spectrum (without OTF mask) and accumulate 908 909 const MaskedArray<Float> spec(fac*itDataVec.vector(), 910 itMaskVec.vector()); 911 vecSum += spec; 912 norm += fac; 913 914 // Next 915 916 itDataVec.next(); 917 itMaskVec.next(); 918 if (wtType==TSYS) pItTsysVec->next(); 919 } 920 921 // Clean up 922 923 if (pItTsysVec) { 924 delete pItTsysVec; 925 pItTsysVec = 0; 926 } 927 } 928 929 // Normalize summed spectrum 930 931 vecSum /= norm; 932 933 // FInd position in input data array. We are iterating by pol-channel 934 // plane so all that will change is beam and IF and that's what we want. 935 936 IPosition pos = itDataPlane.pos(); 937 938 // Write out data. This is a bit messy. We have to reform the Vector 939 // accumulator into an Array of shape (1,1,1,nChan) 940 941 start = pos; 942 end = pos; 943 end(asap::ChanAxis) = nChan-1; 944 outData(start,end) = vecSum.getArray().reform(vecShapeOut); 945 outMask(start,end) = vecSum.getMask().reform(vecShapeOut); 946 947 // Step to next beam/IF combination 948 949 itDataPlane.next(); 950 itMaskPlane.next(); 951 if (wtType==TSYS) pItTsysPlane->next(); 952 } 953 954 // Generate output container and write it to output table 955 956 SDContainer sc = in.getSDContainer(); 957 sc.resize(shapeOut); 958 // 959 putDataInSDC(sc, outData, outMask); 960 pTabOut->putSDContainer(sc); 961 // 962 if (wtType==TSYS) { 963 delete pItTsysPlane; 964 pItTsysPlane = 0; 965 } 966 } 967 968 // Set polarization cursor to 0 969 970 pTabOut->setPol(0); 971 // 972 return pTabOut; 973 } 974 975 976 SDMemTable* SDMath::smooth(const SDMemTable& in, 977 const casa::String& kernelType, 978 casa::Float width, Bool doAll) 979 // 980 // Should smooth TSys as well 981 // 982 { 983 984 // Number of channels 985 const uInt nChan = in.nChan(); 986 987 // Generate Kernel 988 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernelType); 989 Vector<Float> kernel = VectorKernel::make(type, width, nChan, True, False); 990 991 // Generate Convolver 992 IPosition shape(1,nChan); 993 Convolver<Float> conv(kernel, shape); 994 995 // New Table 996 SDMemTable* pTabOut = new SDMemTable(in,True); 997 998 // Output Vectors 999 Vector<Float> valuesOut(nChan); 1000 Vector<Bool> maskOut(nChan); 1001 1002 // Get data slice bounds 1003 IPosition start, end; 1004 setCursorSlice (start, end, doAll, in); 1005 1006 // Loop over rows in Table 1007 for (uInt ri=0; ri < in.nRow(); ++ri) { 1008 1009 // Get slice of data 1010 MaskedArray<Float> dataIn = in.rowAsMaskedArray(ri); 1011 1012 // Deconstruct and get slices which reference these arrays 1013 Array<Float> valuesIn = dataIn.getArray(); 1014 Array<Bool> maskIn = dataIn.getMask(); 1015 1016 Array<Float> valuesIn2 = valuesIn(start,end); // ref to valuesIn 1017 Array<Bool> maskIn2 = maskIn(start,end); 1018 1019 // Iterate through by spectra 1020 VectorIterator<Float> itValues(valuesIn2, asap::ChanAxis); 1021 VectorIterator<Bool> itMask(maskIn2, asap::ChanAxis); 1022 while (!itValues.pastEnd()) { 1023 1024 // Smooth 1025 if (kernelType==VectorKernel::HANNING) { 1026 mathutil::hanning(valuesOut, maskOut, itValues.vector(), 1027 itMask.vector()); 1028 itMask.vector() = maskOut; 1029 } else { 1030 mathutil::replaceMaskByZero(itValues.vector(), itMask.vector()); 1031 conv.linearConv(valuesOut, itValues.vector()); 1032 } 1033 1034 itValues.vector() = valuesOut; 1035 itValues.next(); 1036 itMask.next(); 1037 } 1038 1039 // Create and put back 1040 SDContainer sc = in.getSDContainer(ri); 1041 putDataInSDC(sc, valuesIn, maskIn); 1042 1043 pTabOut->putSDContainer(sc); 1044 } 1045 1046 return pTabOut; 1047 } 1048 1049 1050 1051 SDMemTable* SDMath::convertFlux(const SDMemTable& in, Float D, Float etaAp, 1052 Float JyPerK, Bool doAll) 1053 // 1054 // etaAp = aperture efficiency (-1 means find) 1055 // D = geometric diameter (m) (-1 means find) 1056 // JyPerK 1057 // 1058 { 1059 SDHeader sh = in.getSDHeader(); 1060 SDMemTable* pTabOut = new SDMemTable(in, True); 1061 1062 // Find out how to convert values into Jy and K (e.g. units might be 1063 // mJy or mK) Also automatically find out what we are converting to 1064 // according to the flux unit 1065 Unit fluxUnit(sh.fluxunit); 411 out->frequencies().rescale(width, "RESAMPLE"); 412 TableIterator iter(tout, "IFNO"); 413 TableRow row(tout); 414 while ( !iter.pastEnd() ) { 415 Table tab = iter.table(); 416 ArrayColumn<Float> specCol(tab, "SPECTRA"); 417 //ArrayColumn<Float> tsysCol(tout, "TSYS"); 418 ArrayColumn<uChar> flagCol(tab, "FLAGTRA"); 419 Vector<Float> spec; 420 Vector<uChar> flag; 421 specCol.get(0,spec); // the number of channels should be constant per IF 422 uInt nChanIn = spec.nelements(); 423 Vector<Float> xIn(nChanIn); indgen(xIn); 424 Int fac = Int(nChanIn/width); 425 Vector<Float> xOut(fac+10); // 10 to be safe - resize later 426 uInt k = 0; 427 Float x = 0.0; 428 while (x < Float(nChanIn) ) { 429 xOut(k) = x; 430 k++; 431 x += width; 432 } 433 uInt nChanOut = k; 434 xOut.resize(nChanOut, True); 435 // process all rows for this IFNO 436 Vector<Float> specOut; 437 Vector<Bool> maskOut; 438 Vector<uChar> flagOut; 439 for (uInt i=0; i < tab.nrow(); ++i) { 440 specCol.get(i, spec); 441 flagCol.get(i, flag); 442 Vector<Bool> mask(flag.nelements()); 443 convertArray(mask, flag); 444 445 IPosition shapeIn(spec.shape()); 446 //sh.nchan = nChanOut; 447 InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut, 448 xIn, spec, mask, 449 interpMethod, True, True); 450 /// @todo do the same for channel based Tsys 451 flagOut.resize(maskOut.nelements()); 452 convertArray(flagOut, maskOut); 453 specCol.put(i, specOut); 454 flagCol.put(i, flagOut); 455 } 456 ++iter; 457 } 458 459 return out; 460 } 461 462 STMath::imethod STMath::stringToIMethod(const std::string& in) 463 { 464 static STMath::imap lookup; 465 466 // initialize the lookup table if necessary 467 if ( lookup.empty() ) { 468 lookup["NEAR"] = InterpolateArray1D<Double,Float>::nearestNeighbour; 469 lookup["LIN"] = InterpolateArray1D<Double,Float>::linear; 470 lookup["CUB"] = InterpolateArray1D<Double,Float>::cubic; 471 lookup["SPL"] = InterpolateArray1D<Double,Float>::spline; 472 } 473 474 STMath::imap::const_iterator iter = lookup.find(in); 475 476 if ( lookup.end() == iter ) { 477 std::string message = in; 478 message += " is not a valid interpolation mode"; 479 throw(AipsError(message)); 480 } 481 return iter->second; 482 } 483 484 WeightType STMath::stringToWeight(const std::string& in) 485 { 486 static std::map<std::string, WeightType> lookup; 487 488 // initialize the lookup table if necessary 489 if ( lookup.empty() ) { 490 lookup["NONE"] = asap::NONE; 491 lookup["TINT"] = asap::TINT; 492 lookup["TINTSYS"] = asap::TINTSYS; 493 lookup["TSYS"] = asap::TSYS; 494 lookup["VAR"] = asap::VAR; 495 } 496 497 std::map<std::string, WeightType>::const_iterator iter = lookup.find(in); 498 499 if ( lookup.end() == iter ) { 500 std::string message = in; 501 message += " is not a valid weighting mode"; 502 throw(AipsError(message)); 503 } 504 return iter->second; 505 } 506 507 CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in, 508 const Vector< Float > & coeffs, 509 const std::string & filename, 510 const std::string& method) 511 { 512 // Get elevation data from Scantable and convert to degrees 513 CountedPtr< Scantable > out = getScantable(in, false); 514 Table& tab = in->table(); 515 ROScalarColumn<Float> elev(tab, "ELEVATION"); 516 Vector<Float> x = elev.getColumn(); 517 x *= Float(180 / C::pi); // Degrees 518 519 const uInt nc = coeffs.nelements(); 520 if ( filename.length() > 0 && nc > 0 ) { 521 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both")); 522 } 523 524 // Correct 525 if ( nc > 0 || filename.length() == 0 ) { 526 // Find instrument 527 Bool throwit = True; 528 Instrument inst = 529 SDAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), 530 throwit); 531 532 // Set polynomial 533 Polynomial<Float>* ppoly = 0; 534 Vector<Float> coeff; 535 String msg; 536 if ( nc > 0 ) { 537 ppoly = new Polynomial<Float>(nc); 538 coeff = coeffs; 539 msg = String("user"); 540 } else { 541 SDAttr sdAttr; 542 coeff = sdAttr.gainElevationPoly(inst); 543 ppoly = new Polynomial<Float>(3); 544 msg = String("built in"); 545 } 546 547 if ( coeff.nelements() > 0 ) { 548 ppoly->setCoefficients(coeff); 549 } else { 550 delete ppoly; 551 throw(AipsError("There is no known gain-elevation polynomial known for this instrument")); 552 } 553 ostringstream oss; 554 oss << "Making polynomial correction with " << msg << " coefficients:" << endl; 555 oss << " " << coeff; 556 pushLog(String(oss)); 557 const uInt nrow = tab.nrow(); 558 Vector<Float> factor(nrow); 559 for ( uInt i=0; i < nrow; ++i ) { 560 factor[i] = 1.0 / (*ppoly)(x[i]); 561 } 562 delete ppoly; 563 scaleByVector(tab, factor, true); 564 565 } else { 566 // Read and correct 567 pushLog("Making correction from ascii Table"); 568 scaleFromAsciiTable(tab, filename, method, x, true); 569 } 570 return out; 571 } 572 573 void STMath::scaleFromAsciiTable(Table& in, const std::string& filename, 574 const std::string& method, 575 const Vector<Float>& xout, bool dotsys) 576 { 577 578 // Read gain-elevation ascii file data into a Table. 579 580 String formatString; 581 Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False); 582 scaleFromTable(in, tbl, method, xout, dotsys); 583 } 584 585 void STMath::scaleFromTable(Table& in, 586 const Table& table, 587 const std::string& method, 588 const Vector<Float>& xout, bool dotsys) 589 { 590 591 ROScalarColumn<Float> geElCol(table, "ELEVATION"); 592 ROScalarColumn<Float> geFacCol(table, "FACTOR"); 593 Vector<Float> xin = geElCol.getColumn(); 594 Vector<Float> yin = geFacCol.getColumn(); 595 Vector<Bool> maskin(xin.nelements(),True); 596 597 // Interpolate (and extrapolate) with desired method 598 599 //InterpolateArray1D<Double,Float>::InterpolationMethod method; 600 Int intmethod(stringToIMethod(method)); 601 602 Vector<Float> yout; 603 Vector<Bool> maskout; 604 InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout, 605 xin, yin, maskin, intmethod, 606 True, True); 607 608 scaleByVector(in, Float(1.0)/yout, dotsys); 609 } 610 611 void STMath::scaleByVector( Table& in, 612 const Vector< Float >& factor, 613 bool dotsys ) 614 { 615 uInt nrow = in.nrow(); 616 if ( factor.nelements() != nrow ) { 617 throw(AipsError("factors.nelements() != table.nelements()")); 618 } 619 ArrayColumn<Float> specCol(in, "SPECTRA"); 620 ArrayColumn<uChar> flagCol(in, "FLAGTRA"); 621 ArrayColumn<Float> tsysCol(in, "TSYS"); 622 for (uInt i=0; i < nrow; ++i) { 623 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i)); 624 ma *= factor[i]; 625 specCol.put(i, ma.getArray()); 626 flagCol.put(i, flagsFromMA(ma)); 627 if ( dotsys ) { 628 Vector<Float> tsys; 629 tsysCol.get(i, tsys); 630 tsys *= factor[i]; 631 specCol.put(i,tsys); 632 } 633 } 634 } 635 636 CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in, 637 float d, float etaap, 638 float jyperk ) 639 { 640 CountedPtr< Scantable > out = getScantable(in, false); 641 Table& tab = in->table(); 642 Unit fluxUnit(tab.keywordSet().asString("FluxUnit")); 1066 643 Unit K(String("K")); 1067 644 Unit JY(String("Jy")); 1068 645 1069 Bool toKelvin = True;1070 Double c Fac = 1.0;1071 1072 if ( fluxUnit==JY) {646 bool tokelvin = true; 647 Double cfac = 1.0; 648 649 if ( fluxUnit == JY ) { 1073 650 pushLog("Converting to K"); 1074 651 Quantum<Double> t(1.0,fluxUnit); 1075 652 Quantum<Double> t2 = t.get(JY); 1076 c Fac = (t2 / t).getValue(); // value to Jy1077 1078 to Kelvin = True;1079 sh.fluxunit = "K";1080 } else if ( fluxUnit==K) {653 cfac = (t2 / t).getValue(); // value to Jy 654 655 tokelvin = true; 656 out->setFluxUnit("K"); 657 } else if ( fluxUnit == K ) { 1081 658 pushLog("Converting to Jy"); 1082 659 Quantum<Double> t(1.0,fluxUnit); 1083 660 Quantum<Double> t2 = t.get(K); 1084 c Fac = (t2 / t).getValue(); // value to K1085 1086 to Kelvin = False;1087 sh.fluxunit = "Jy";661 cfac = (t2 / t).getValue(); // value to K 662 663 tokelvin = false; 664 out->setFluxUnit("Jy"); 1088 665 } else { 1089 666 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K")); 1090 667 } 1091 1092 pTabOut->putSDHeader(sh);1093 1094 668 // Make sure input values are converted to either Jy or K first... 1095 Float factor = c Fac;669 Float factor = cfac; 1096 670 1097 671 // Select method 1098 if ( JyPerK>0.0) {1099 factor *= JyPerK;1100 if ( toKelvin) factor = 1.0 / JyPerK;672 if (jyperk > 0.0) { 673 factor *= jyperk; 674 if ( tokelvin ) factor = 1.0 / jyperk; 1101 675 ostringstream oss; 1102 oss << "Jy/K = " << JyPerK;676 oss << "Jy/K = " << jyperk; 1103 677 pushLog(String(oss)); 1104 Vector<Float> factors( in.nRow(), factor);1105 scaleByVector( pTabOut, in, doAll, factors, False);1106 } else if ( etaAp>0.0) {1107 Bool throwIt = True;1108 Instrument inst = SDAttr::convertInstrument(sh.antennaname, throwIt);678 Vector<Float> factors(tab.nrow(), factor); 679 scaleByVector(tab,factors, false); 680 } else if ( etaap > 0.0) { 681 Instrument inst = 682 SDAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), True); 1109 683 SDAttr sda; 1110 if ( D < 0) D= sda.diameter(inst);1111 Float JyPerK = SDAttr::findJyPerK(etaAp,D);684 if (d < 0) d = sda.diameter(inst); 685 Float jyPerk = SDAttr::findJyPerK(etaap, d); 1112 686 ostringstream oss; 1113 oss << "Jy/K = " << JyPerK;687 oss << "Jy/K = " << jyperk; 1114 688 pushLog(String(oss)); 1115 factor *= JyPerK;1116 if ( toKelvin) {689 factor *= jyperk; 690 if ( tokelvin ) { 1117 691 factor = 1.0 / factor; 1118 692 } 1119 1120 Vector<Float> factors(in.nRow(), factor); 1121 scaleByVector(pTabOut, in, doAll, factors, False); 693 Vector<Float> factors(tab.nrow(), factor); 694 scaleByVector(tab, factors, False); 1122 695 } else { 1123 696 … … 1128 701 1129 702 pushLog("Looking up conversion factors"); 1130 convertBrightnessUnits (pTabOut, in, toKelvin, cFac, doAll); 1131 } 1132 return pTabOut; 1133 } 1134 1135 1136 SDMemTable* SDMath::gainElevation(const SDMemTable& in, 1137 const Vector<Float>& coeffs, 1138 const String& fileName, 1139 const String& methodStr, Bool doAll) 1140 { 1141 1142 // Get header and clone output table 1143 SDHeader sh = in.getSDHeader(); 1144 SDMemTable* pTabOut = new SDMemTable(in, True); 1145 1146 // Get elevation data from SDMemTable and convert to degrees 1147 const Table& tab = in.table(); 703 convertBrightnessUnits(out, tokelvin, cfac); 704 } 705 706 return out; 707 } 708 709 void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in, 710 bool tokelvin, float cfac ) 711 { 712 Table& table = in->table(); 713 Instrument inst = 714 SDAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True); 715 TableIterator iter(table, "FREQ_ID"); 716 STFrequencies stfreqs = in->frequencies(); 717 SDAttr sdAtt; 718 while (!iter.pastEnd()) { 719 Table tab = iter.table(); 720 ArrayColumn<Float> specCol(tab, "SPECTRA"); 721 ArrayColumn<uChar> flagCol(tab, "FLAGTRA"); 722 ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID"); 723 MEpoch::ROScalarColumn timeCol(tab, "TIME"); 724 725 uInt freqid; freqidCol.get(0, freqid); 726 Vector<Float> tmpspec; specCol.get(0, tmpspec); 727 // SDAttr.JyPerK has a Vector interface... change sometime. 728 Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements())); 729 for ( uInt i=0; i<tab.nrow(); ++i) { 730 Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0]; 731 Float factor = cfac * jyperk; 732 if ( tokelvin ) factor = Float(1.0) / factor; 733 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i)); 734 ma *= factor; 735 specCol.put(i, ma.getArray()); 736 flagCol.put(i, flagsFromMA(ma)); 737 } 738 } 739 } 740 741 CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in, 742 float tau ) 743 { 744 CountedPtr< Scantable > out = getScantable(in, false); 745 Table& tab = in->table(); 1148 746 ROScalarColumn<Float> elev(tab, "ELEVATION"); 1149 Vector<Float> x = elev.getColumn(); 1150 x *= Float(180 / C::pi); // Degrees 1151 1152 const uInt nC = coeffs.nelements(); 1153 if (fileName.length()>0 && nC>0) { 1154 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both")); 1155 } 1156 1157 // Correct 1158 if (nC>0 || fileName.length()==0) { 1159 // Find instrument 1160 Bool throwIt = True; 1161 Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt); 1162 1163 // Set polynomial 1164 Polynomial<Float>* pPoly = 0; 1165 Vector<Float> coeff; 1166 String msg; 1167 if (nC>0) { 1168 pPoly = new Polynomial<Float>(nC); 1169 coeff = coeffs; 1170 msg = String("user"); 1171 } else { 1172 SDAttr sdAttr; 1173 coeff = sdAttr.gainElevationPoly(inst); 1174 pPoly = new Polynomial<Float>(3); 1175 msg = String("built in"); 1176 } 1177 1178 if (coeff.nelements()>0) { 1179 pPoly->setCoefficients(coeff); 1180 } else { 1181 delete pPoly; 1182 throw(AipsError("There is no known gain-elevation polynomial known for this instrument")); 1183 } 1184 ostringstream oss; 1185 oss << "Making polynomial correction with " << msg << " coefficients:" << endl; 1186 oss << " " << coeff; 1187 pushLog(String(oss)); 1188 const uInt nRow = in.nRow(); 1189 Vector<Float> factor(nRow); 1190 for (uInt i=0; i<nRow; i++) { 1191 factor[i] = 1.0 / (*pPoly)(x[i]); 1192 } 1193 delete pPoly; 1194 scaleByVector (pTabOut, in, doAll, factor, True); 1195 1196 } else { 1197 1198 // Indicate which columns to read from ascii file 1199 String col0("ELEVATION"); 1200 String col1("FACTOR"); 1201 1202 // Read and correct 1203 1204 pushLog("Making correction from ascii Table"); 1205 scaleFromAsciiTable (pTabOut, in, fileName, col0, col1, 1206 methodStr, doAll, x, True); 1207 } 1208 1209 return pTabOut; 1210 } 1211 1212 1213 SDMemTable* SDMath::opacity(const SDMemTable& in, Float tau, Bool doAll) 1214 { 1215 1216 // Get header and clone output table 1217 1218 SDHeader sh = in.getSDHeader(); 1219 SDMemTable* pTabOut = new SDMemTable(in, True); 1220 1221 // Get elevation data from SDMemTable and convert to degrees 1222 1223 const Table& tab = in.table(); 1224 ROScalarColumn<Float> elev(tab, "ELEVATION"); 1225 Vector<Float> zDist = elev.getColumn(); 1226 zDist = Float(C::pi_2) - zDist; 1227 1228 // Generate correction factor 1229 1230 const uInt nRow = in.nRow(); 1231 Vector<Float> factor(nRow); 1232 Vector<Float> factor2(nRow); 1233 for (uInt i=0; i<nRow; i++) { 1234 factor[i] = exp(tau)/cos(zDist[i]); 1235 } 1236 1237 // Correct 1238 1239 scaleByVector (pTabOut, in, doAll, factor, True); 1240 1241 return pTabOut; 1242 } 1243 1244 1245 void SDMath::rotateXYPhase(SDMemTable& in, Float value, Bool doAll) 1246 // 1247 // phase in degrees 1248 // assumes linear correlations 1249 // 1250 { 1251 if (in.nPol() != 4) { 1252 throw(AipsError("You must have 4 polarizations to run this function")); 1253 } 1254 1255 SDHeader sh = in.getSDHeader(); 1256 Instrument inst = SDAttr::convertInstrument(sh.antennaname, False); 1257 SDAttr sdAtt; 1258 if (sdAtt.feedPolType(inst) != LINEAR) { 1259 throw(AipsError("Only linear polarizations are supported")); 1260 } 1261 // 1262 const Table& tabIn = in.table(); 1263 ArrayColumn<Float> specCol(tabIn,"SPECTRA"); 1264 IPosition start(asap::nAxes,0); 1265 IPosition end(asap::nAxes); 1266 1267 // Set cursor slice. Assumes shape the same for all rows 1268 1269 setCursorSlice (start, end, doAll, in); 1270 IPosition start3(start); 1271 start3(asap::PolAxis) = 2; // Real(XY) 1272 IPosition end3(end); 1273 end3(asap::PolAxis) = 2; 1274 // 1275 IPosition start4(start); 1276 start4(asap::PolAxis) = 3; // Imag (XY) 1277 IPosition end4(end); 1278 end4(asap::PolAxis) = 3; 1279 // 1280 uInt nRow = in.nRow(); 1281 Array<Float> data; 1282 for (uInt i=0; i<nRow;++i) { 1283 specCol.get(i,data); 1284 IPosition shape = data.shape(); 1285 1286 // Get polarization slice references 1287 Array<Float> C3 = data(start3,end3); 1288 Array<Float> C4 = data(start4,end4); 1289 1290 // Rotate 1291 SDPolUtil::rotatePhase(C3, C4, value); 1292 1293 // Put 1294 specCol.put(i,data); 1295 } 1296 } 1297 1298 1299 void SDMath::rotateLinPolPhase(SDMemTable& in, Float value, Bool doAll) 1300 // 1301 // phase in degrees 1302 // assumes linear correlations 1303 // 1304 { 1305 if (in.nPol() != 4) { 1306 throw(AipsError("You must have 4 polarizations to run this function")); 1307 } 1308 // 1309 SDHeader sh = in.getSDHeader(); 1310 Instrument inst = SDAttr::convertInstrument(sh.antennaname, False); 1311 SDAttr sdAtt; 1312 if (sdAtt.feedPolType(inst) != LINEAR) { 1313 throw(AipsError("Only linear polarizations are supported")); 1314 } 1315 // 1316 const Table& tabIn = in.table(); 1317 ArrayColumn<Float> specCol(tabIn,"SPECTRA"); 1318 ROArrayColumn<Float> stokesCol(tabIn,"STOKES"); 1319 IPosition start(asap::nAxes,0); 1320 IPosition end(asap::nAxes); 1321 1322 // Set cursor slice. Assumes shape the same for all rows 1323 1324 setCursorSlice (start, end, doAll, in); 1325 // 1326 IPosition start1(start); 1327 start1(asap::PolAxis) = 0; // C1 (XX) 1328 IPosition end1(end); 1329 end1(asap::PolAxis) = 0; 1330 // 1331 IPosition start2(start); 1332 start2(asap::PolAxis) = 1; // C2 (YY) 1333 IPosition end2(end); 1334 end2(asap::PolAxis) = 1; 1335 // 1336 IPosition start3(start); 1337 start3(asap::PolAxis) = 2; // C3 ( Real(XY) ) 1338 IPosition end3(end); 1339 end3(asap::PolAxis) = 2; 1340 // 1341 IPosition startI(start); 1342 startI(asap::PolAxis) = 0; // I 1343 IPosition endI(end); 1344 endI(asap::PolAxis) = 0; 1345 // 1346 IPosition startQ(start); 1347 startQ(asap::PolAxis) = 1; // Q 1348 IPosition endQ(end); 1349 endQ(asap::PolAxis) = 1; 1350 // 1351 IPosition startU(start); 1352 startU(asap::PolAxis) = 2; // U 1353 IPosition endU(end); 1354 endU(asap::PolAxis) = 2; 1355 1356 // 1357 uInt nRow = in.nRow(); 1358 Array<Float> data, stokes; 1359 for (uInt i=0; i<nRow;++i) { 1360 specCol.get(i,data); 1361 stokesCol.get(i,stokes); 1362 IPosition shape = data.shape(); 1363 1364 // Get linear polarization slice references 1365 1366 Array<Float> C1 = data(start1,end1); 1367 Array<Float> C2 = data(start2,end2); 1368 Array<Float> C3 = data(start3,end3); 1369 1370 // Get STokes slice references 1371 1372 Array<Float> I = stokes(startI,endI); 1373 Array<Float> Q = stokes(startQ,endQ); 1374 Array<Float> U = stokes(startU,endU); 1375 1376 // Rotate 1377 1378 SDPolUtil::rotateLinPolPhase(C1, C2, C3, I, Q, U, value); 1379 1380 // Put 1381 1382 specCol.put(i,data); 1383 } 1384 } 1385 1386 // 'private' functions 1387 1388 void SDMath::convertBrightnessUnits (SDMemTable* pTabOut, const SDMemTable& in, 1389 Bool toKelvin, Float cFac, Bool doAll) 1390 { 1391 1392 // Get header 1393 1394 SDHeader sh = in.getSDHeader(); 1395 const uInt nChan = sh.nchan; 1396 1397 // Get instrument 1398 1399 Bool throwIt = True; 1400 Instrument inst = SDAttr::convertInstrument(sh.antennaname, throwIt); 1401 1402 // Get Diameter (m) 1403 1404 SDAttr sdAtt; 1405 // Get epoch of first row 1406 1407 MEpoch dateObs = in.getEpoch(0); 1408 1409 // Generate a Vector of correction factors. One per FreqID 1410 1411 SDFrequencyTable sdft = in.getSDFreqTable(); 1412 Vector<uInt> freqIDs; 1413 // 1414 Vector<Float> freqs(sdft.length()); 1415 for (uInt i=0; i<sdft.length(); i++) { 1416 freqs(i) = (nChan/2 - sdft.referencePixel(i))*sdft.increment(i) + sdft.referenceValue(i); 1417 } 1418 // 1419 Vector<Float> JyPerK = sdAtt.JyPerK(inst, dateObs, freqs); 1420 ostringstream oss; 1421 oss << "Jy/K = " << JyPerK; 1422 pushLog(String(oss)); 1423 Vector<Float> factors = cFac * JyPerK; 1424 if (toKelvin) factors = Float(1.0) / factors; 1425 1426 // Get data slice bounds 1427 1428 IPosition start, end; 1429 setCursorSlice (start, end, doAll, in); 1430 const uInt ifAxis = in.getIF(); 1431 1432 // Iteration axes 1433 1434 IPosition axes(asap::nAxes-1,0); 1435 for (uInt i=0,j=0; i<asap::nAxes; i++) { 1436 if (i!=asap::IFAxis) { 1437 axes(j++) = i; 747 ArrayColumn<Float> specCol(tab, "SPECTRA"); 748 ArrayColumn<uChar> flagCol(tab, "FLAGTRA"); 749 for ( uInt i=0; i<tab.nrow(); ++i) { 750 Float zdist = Float(C::pi_2) - elev(i); 751 Float factor = exp(tau)/cos(zdist); 752 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i)); 753 ma *= factor; 754 specCol.put(i, ma.getArray()); 755 flagCol.put(i, flagsFromMA(ma)); 756 } 757 return out; 758 } 759 760 CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in, 761 const std::string& kernel, float width ) 762 { 763 CountedPtr< Scantable > out = getScantable(in, false); 764 Table& table = in->table(); 765 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel); 766 // same IFNO should have same no of channels 767 // this saves overhead 768 TableIterator iter(table, "IFNO"); 769 while (!iter.pastEnd()) { 770 Table tab = iter.table(); 771 ArrayColumn<Float> specCol(tab, "SPECTRA"); 772 ArrayColumn<uChar> flagCol(tab, "FLAGTRA"); 773 Vector<Float> tmpspec; specCol.get(0, tmpspec); 774 uInt nchan = tmpspec.nelements(); 775 Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False); 776 Convolver<Float> conv(kvec, IPosition(1,nchan)); 777 Vector<Float> spec; 778 Vector<uChar> flag; 779 for ( uInt i=0; i<tab.nrow(); ++i) { 780 specCol.get(i, spec); 781 flagCol.get(i, flag); 782 Vector<Bool> mask(flag.nelements()); 783 convertArray(mask, flag); 784 Vector<Float> specout; 785 if ( type == VectorKernel::HANNING ) { 786 Vector<Bool> maskout; 787 mathutil::hanning(specout, maskout, spec , mask); 788 convertArray(flag, maskout); 789 flagCol.put(i, flag); 790 } else { 791 mathutil::replaceMaskByZero(specout, mask); 792 conv.linearConv(specout, spec); 1438 793 } 1439 } 1440 1441 // Loop over rows and apply correction factor 1442 1443 Float factor = 1.0; 1444 const uInt axis = asap::ChanAxis; 1445 for (uInt i=0; i < in.nRow(); ++i) { 1446 1447 // Get data 1448 1449 MaskedArray<Float> dataIn = in.rowAsMaskedArray(i); 1450 Array<Float>& values = dataIn.getRWArray(); // Ref to dataIn 1451 Array<Float> values2 = values(start,end); // Ref to values to dataIn 1452 1453 // Get SDCOntainer 1454 1455 SDContainer sc = in.getSDContainer(i); 1456 1457 // Get FreqIDs 1458 1459 freqIDs = sc.getFreqMap(); 1460 1461 // Now the conversion factor depends only upon frequency 1462 // So we need to iterate through by IF only giving 1463 // us BEAM/POL/CHAN cubes 1464 1465 ArrayIterator<Float> itIn(values2, axes); 1466 uInt ax = 0; 1467 while (!itIn.pastEnd()) { 1468 itIn.array() *= factors(freqIDs(ax)); // Writes back to dataIn 1469 itIn.next(); 1470 } 1471 1472 // Write out 1473 1474 putDataInSDC(sc, dataIn.getArray(), dataIn.getMask()); 1475 // 1476 pTabOut->putSDContainer(sc); 1477 } 1478 } 1479 1480 1481 1482 SDMemTable* SDMath::frequencyAlign(const SDMemTable& in, 1483 MFrequency::Types freqSystem, 1484 const String& refTime, 1485 const String& methodStr, 1486 Bool perFreqID) 1487 { 1488 // Get Header 1489 1490 SDHeader sh = in.getSDHeader(); 1491 const uInt nChan = sh.nchan; 1492 const uInt nRows = in.nRow(); 1493 const uInt nIF = sh.nif; 1494 1495 // Get Table reference 1496 1497 const Table& tabIn = in.table(); 1498 1499 // Get Columns from Table 1500 1501 ROScalarColumn<Double> mjdCol(tabIn, "TIME"); 1502 ROScalarColumn<String> srcCol(tabIn, "SRCNAME"); 1503 ROArrayColumn<uInt> fqIDCol(tabIn, "FREQID"); 1504 Vector<Double> times = mjdCol.getColumn(); 1505 1506 // Generate DataDesc table 1507 1508 Matrix<uInt> ddIdx; 1509 SDDataDesc dDesc; 1510 generateDataDescTable(ddIdx, dDesc, nIF, in, tabIn, srcCol, 1511 fqIDCol, perFreqID); 1512 1513 // Get reference Epoch to time of first row or given String 1514 1515 Unit DAY(String("d")); 1516 MEpoch::Ref epochRef(in.getTimeReference()); 1517 MEpoch refEpoch; 1518 if (refTime.length()>0) { 1519 refEpoch = epochFromString(refTime, in.getTimeReference()); 1520 } else { 1521 refEpoch = in.getEpoch(0); 1522 } 1523 ostringstream oss; 1524 oss << "Aligned at reference Epoch " << formatEpoch(refEpoch) << " in frame " << MFrequency::showType(freqSystem); 1525 pushLog(String(oss)); 1526 1527 // Get Reference Position 1528 1529 MPosition refPos = in.getAntennaPosition(); 1530 1531 // Create FrequencyAligner Block. One FA for each possible 1532 // source/freqID (perFreqID=True) or source/IF (perFreqID=False) 1533 // combination 1534 1535 PtrBlock<FrequencyAligner<Float>* > a(dDesc.length()); 1536 generateFrequencyAligners(a, dDesc, in, nChan, freqSystem, refPos, 1537 refEpoch, perFreqID); 1538 1539 // Generate and fill output Frequency Table. WHen perFreqID=True, 1540 // there is one output FreqID for each entry in the SDDataDesc 1541 // table. However, in perFreqID=False mode, there may be some 1542 // degeneracy, so we need a little translation map 1543 1544 SDFrequencyTable freqTabOut = in.getSDFreqTable(); 1545 freqTabOut.setLength(0); 1546 Vector<String> units(1); 1547 units = String("Hz"); 1548 Bool linear=True; 1549 // 1550 Vector<uInt> ddFQTrans(dDesc.length(),0); 1551 for (uInt i=0; i<dDesc.length(); i++) { 1552 1553 // Get Aligned SC in Hz 1554 1555 SpectralCoordinate sC = a[i]->alignedSpectralCoordinate(linear); 1556 sC.setWorldAxisUnits(units); 1557 1558 // Add FreqID 1559 1560 uInt idx = freqTabOut.addFrequency(sC.referencePixel()[0], 1561 sC.referenceValue()[0], 1562 sC.increment()[0]); 1563 // output FreqID = ddFQTrans(ddIdx) 1564 ddFQTrans(i) = idx; 1565 } 1566 1567 // Interpolation method 1568 1569 InterpolateArray1D<Double,Float>::InterpolationMethod interp; 1570 convertInterpString(interp, methodStr); 1571 1572 // New output Table 1573 1574 //pushLog("Create output table"); 1575 SDMemTable* pTabOut = new SDMemTable(in,True); 1576 pTabOut->putSDFreqTable(freqTabOut); 1577 1578 // Loop over rows in Table 1579 1580 Bool extrapolate=False; 1581 const IPosition polChanAxes(2, asap::PolAxis, asap::ChanAxis); 1582 Bool useCachedAbcissa = False; 1583 Bool first = True; 1584 Bool ok; 1585 Vector<Float> yOut; 1586 Vector<Bool> maskOut; 1587 Vector<uInt> freqID(nIF); 1588 uInt ifIdx, faIdx; 1589 Vector<Double> xIn; 1590 // 1591 for (uInt iRow=0; iRow<nRows; ++iRow) { 1592 // if (iRow%10==0) { 1593 // ostringstream oss; 1594 // oss << "Processing row " << iRow; 1595 // pushLog(String(oss)); 1596 // } 1597 1598 // Get EPoch 1599 1600 Quantum<Double> tQ2(times[iRow],DAY); 1601 MVEpoch mv2(tQ2); 1602 MEpoch epoch(mv2, epochRef); 1603 1604 // Get copy of data 1605 1606 const MaskedArray<Float>& mArrIn(in.rowAsMaskedArray(iRow)); 1607 Array<Float> values = mArrIn.getArray(); 1608 Array<Bool> mask = mArrIn.getMask(); 1609 1610 // For each row, the Frequency abcissa will be the same 1611 // regardless of polarization. For all other axes (IF and BEAM) 1612 // the abcissa will change. So we iterate through the data by 1613 // pol-chan planes to mimimize the work. Probably won't work for 1614 // multiple beams at this point. 1615 1616 ArrayIterator<Float> itValuesPlane(values, polChanAxes); 1617 ArrayIterator<Bool> itMaskPlane(mask, polChanAxes); 1618 while (!itValuesPlane.pastEnd()) { 1619 1620 // Find the IF index and then the FA PtrBlock index 1621 1622 const IPosition& pos = itValuesPlane.pos(); 1623 ifIdx = pos(asap::IFAxis); 1624 faIdx = ddIdx(iRow,ifIdx); 1625 1626 // Generate abcissa for perIF. Could cache this in a Matrix on 1627 // a per scan basis. Pretty expensive doing it for every row. 1628 1629 if (!perFreqID) { 1630 xIn.resize(nChan); 1631 uInt fqID = dDesc.secID(ddIdx(iRow,ifIdx)); 1632 SpectralCoordinate sC = in.getSpectralCoordinate(fqID); 1633 Double w; 1634 for (uInt i=0; i<nChan; i++) { 1635 sC.toWorld(w,Double(i)); 1636 xIn[i] = w; 1637 } 1638 } 1639 1640 VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1); 1641 VectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1); 1642 1643 // Iterate through the plane by vector and align 1644 1645 first = True; 1646 useCachedAbcissa=False; 1647 while (!itValuesVec.pastEnd()) { 1648 if (perFreqID) { 1649 ok = a[faIdx]->align (yOut, maskOut, itValuesVec.vector(), 1650 itMaskVec.vector(), epoch, useCachedAbcissa, 1651 interp, extrapolate); 1652 } else { 1653 ok = a[faIdx]->align (yOut, maskOut, xIn, itValuesVec.vector(), 1654 itMaskVec.vector(), epoch, useCachedAbcissa, 1655 interp, extrapolate); 1656 } 1657 // 1658 itValuesVec.vector() = yOut; 1659 itMaskVec.vector() = maskOut; 1660 // 1661 itValuesVec.next(); 1662 itMaskVec.next(); 1663 // 1664 if (first) { 1665 useCachedAbcissa = True; 1666 first = False; 1667 } 1668 } 1669 // 1670 itValuesPlane.next(); 1671 itMaskPlane.next(); 1672 } 1673 1674 // Create SDContainer and put back 1675 1676 SDContainer sc = in.getSDContainer(iRow); 1677 putDataInSDC(sc, values, mask); 1678 1679 // Set output FreqIDs 1680 1681 for (uInt i=0; i<nIF; i++) { 1682 uInt idx = ddIdx(iRow,i); // Index into SDDataDesc table 1683 freqID(i) = ddFQTrans(idx); // FreqID in output FQ table 1684 } 1685 sc.putFreqMap(freqID); 1686 // 1687 pTabOut->putSDContainer(sc); 1688 } 1689 1690 // Now we must set the base and extra frames to the input frame 1691 std::vector<string> info = pTabOut->getCoordInfo(); 1692 info[1] = MFrequency::showType(freqSystem); // Conversion frame 1693 info[3] = info[1]; // Base frame 1694 pTabOut->setCoordInfo(info); 1695 1696 // Clean up PointerBlock 1697 for (uInt i=0; i<a.nelements(); i++) delete a[i]; 1698 1699 return pTabOut; 1700 } 1701 1702 1703 SDMemTable* SDMath::frequencySwitch(const SDMemTable& in) 1704 { 1705 if (in.nIF() != 2) { 1706 throw(AipsError("nIF != 2 ")); 1707 } 1708 Bool clear = True; 1709 SDMemTable* pTabOut = new SDMemTable(in, clear); 1710 const Table& tabIn = in.table(); 1711 1712 // Shape of input and output data 1713 const IPosition& shapeIn = in.rowAsMaskedArray(0).shape(); 1714 1715 const uInt nRows = in.nRow(); 1716 AlwaysAssert(asap::nAxes==4,AipsError); 1717 1718 ROArrayColumn<Float> tSysCol; 1719 Array<Float> tsys; 1720 tSysCol.attach(tabIn,"TSYS"); 1721 1722 for (uInt iRow=0; iRow<nRows; iRow++) { 1723 // Get data for this row 1724 MaskedArray<Float> marr(in.rowAsMaskedArray(iRow)); 1725 tSysCol.get(iRow, tsys); 1726 1727 // whole Array for IF 0 1728 IPosition start(asap::nAxes,0); 1729 IPosition end = shapeIn-1; 1730 end(asap::IFAxis) = 0; 1731 1732 MaskedArray<Float> on = marr(start,end); 1733 Array<Float> ton = tsys(start,end); 1734 // Make a copy as "src" is a refrence which is manipulated. 1735 // oncopy is needed for the inverse quotient 1736 MaskedArray<Float> oncopy = on.copy(); 1737 1738 // whole Array for IF 1 1739 start(asap::IFAxis) = 1; 1740 end(asap::IFAxis) = 1; 1741 1742 MaskedArray<Float> off = marr(start,end); 1743 Array<Float> toff = tsys(start,end); 1744 1745 on /= off; on -= 1.0f; 1746 on *= ton; 1747 off /= oncopy; off -= 1.0f; 1748 off *= toff; 1749 1750 SDContainer sc = in.getSDContainer(iRow); 1751 putDataInSDC(sc, marr.getArray(), marr.getMask()); 1752 pTabOut->putSDContainer(sc); 1753 } 1754 return pTabOut; 1755 } 1756 1757 void SDMath::fillSDC(SDContainer& sc, 1758 const Array<Bool>& mask, 1759 const Array<Float>& data, 1760 const Array<Float>& tSys, 1761 Int scanID, Double timeStamp, 1762 Double interval, const String& sourceName, 1763 const Vector<uInt>& freqID) 1764 { 1765 // Data and mask 1766 1767 putDataInSDC(sc, data, mask); 1768 1769 // TSys 1770 1771 sc.putTsys(tSys); 1772 1773 // Time things 1774 1775 sc.timestamp = timeStamp; 1776 sc.interval = interval; 1777 sc.scanid = scanID; 1778 // 1779 sc.sourcename = sourceName; 1780 sc.putFreqMap(freqID); 1781 } 1782 1783 void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum, 1784 MaskedArray<Float>& sum, Array<Float>& sumSq, 1785 Array<Float>& nPts, Array<Float>& tSysSum, 1786 Array<Float>& tSysSqSum, 1787 const Array<Float>& tSys, const Array<Float>& nInc, 1788 const Vector<Bool>& mask, Double time, Double interval, 1789 const std::vector<CountedPtr<SDMemTable> >& in, 1790 uInt iTab, uInt iRow, uInt axis, 1791 uInt nAxesSub, Bool useMask, 1792 WeightType wtType) 1793 { 1794 1795 // Get data 1796 1797 MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow)); 1798 Array<Float>& valuesIn = dataIn.getRWArray(); // writable reference 1799 const Array<Bool>& maskIn = dataIn.getMask(); // RO reference 1800 // 1801 if (wtType==NONE) { 1802 const MaskedArray<Float> n(nInc,dataIn.getMask()); 1803 nPts += n; // Only accumulates where mask==T 1804 } else if (wtType==TINT) { 1805 1806 // We are weighting the data by integration time. 1807 1808 valuesIn *= Float(interval); 1809 1810 } else if (wtType==VAR) { 1811 1812 // We are going to average the data, weighted by the noise for each pol, beam and IF. 1813 // So therefore we need to iterate through by spectrum (axis 3) 1814 1815 VectorIterator<Float> itData(valuesIn, axis); 1816 ReadOnlyVectorIterator<Bool> itMask(maskIn, axis); 1817 Float fac = 1.0; 1818 IPosition pos(nAxesSub,0); 1819 // 1820 while (!itData.pastEnd()) { 1821 1822 // Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor 1823 1824 if (useMask) { 1825 MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector()); 1826 fac = 1.0/variance(tmp); 1827 } else { 1828 MaskedArray<Float> tmp(itData.vector(),itMask.vector()); 1829 fac = 1.0/variance(tmp); 1830 } 1831 1832 // Scale data 1833 1834 itData.vector() *= fac; // Writes back into 'dataIn' 1835 // 1836 // Accumulate variance per if/pol/beam averaged over spectrum 1837 // This method to get pos2 from itData.pos() is only valid 1838 // because the spectral axis is the last one (so we can just 1839 // copy the first nAXesSub positions out) 1840 1841 pos = itData.pos().getFirst(nAxesSub); 1842 sumSq(pos) += fac; 1843 // 1844 itData.next(); 1845 itMask.next(); 1846 } 1847 } else if (wtType==TSYS || wtType==TINTSYS) { 1848 1849 // We are going to average the data, weighted by 1/Tsys**2 for each pol, beam and IF. 1850 // So therefore we need to iterate through by spectrum (axis 3). Although 1851 // Tsys is stored as a vector of length nChan, the values are replicated. 1852 // We will take a short cut and just use the value from the first channel 1853 // for now. 1854 // 1855 VectorIterator<Float> itData(valuesIn, axis); 1856 ReadOnlyVectorIterator<Float> itTSys(tSys, axis); 1857 IPosition pos(nAxesSub,0); 1858 // 1859 Float fac = 1.0; 1860 if (wtType==TINTSYS) fac *= interval; 1861 while (!itData.pastEnd()) { 1862 Float t = itTSys.vector()[0]; 1863 fac *= 1.0/t/t; 1864 1865 // Scale data 1866 1867 itData.vector() *= fac; // Writes back into 'dataIn' 1868 // 1869 // Accumulate Tsys per if/pol/beam averaged over spectrum 1870 // This method to get pos2 from itData.pos() is only valid 1871 // because the spectral axis is the last one (so we can just 1872 // copy the first nAXesSub positions out) 1873 1874 pos = itData.pos().getFirst(nAxesSub); 1875 tSysSqSum(pos) += fac; 1876 // 1877 itData.next(); 1878 itTSys.next(); 1879 } 1880 } 1881 1882 // Accumulate sum of (possibly scaled) data 1883 1884 sum += dataIn; 1885 1886 // Accumulate Tsys, time, and interval 1887 1888 tSysSum += tSys; 1889 timeSum += time; 1890 intSum += interval; 1891 nAccum += 1; 1892 } 1893 1894 1895 void SDMath::normalize(MaskedArray<Float>& sum, 1896 const Array<Float>& sumSq, 1897 const Array<Float>& tSysSqSum, 1898 const Array<Float>& nPts, 1899 Double intSum, 1900 WeightType wtType, Int axis, 1901 Int nAxesSub) 1902 { 1903 IPosition pos2(nAxesSub,0); 1904 // 1905 if (wtType==NONE) { 1906 1907 // We just average by the number of points accumulated. 1908 // We need to make a MA out of nPts so that no divide by 1909 // zeros occur 1910 1911 MaskedArray<Float> t(nPts, (nPts>Float(0.0))); 1912 sum /= t; 1913 } else if (wtType==TINT) { 1914 1915 // Average by sum of Tint 1916 1917 sum /= Float(intSum); 1918 } else if (wtType==VAR) { 1919 1920 // Normalize each spectrum by sum(1/var) where the variance 1921 // is worked out for each spectrum 1922 1923 Array<Float>& data = sum.getRWArray(); 1924 VectorIterator<Float> itData(data, axis); 1925 while (!itData.pastEnd()) { 1926 pos2 = itData.pos().getFirst(nAxesSub); 1927 itData.vector() /= sumSq(pos2); 1928 itData.next(); 1929 } 1930 } else if (wtType==TSYS || wtType==TINTSYS) { 1931 1932 // Normalize each spectrum by sum(1/Tsys**2) (TSYS) or 1933 // sum(Tint/Tsys**2) (TINTSYS) where the pseudo 1934 // replication over channel for Tsys has been dropped. 1935 1936 Array<Float>& data = sum.getRWArray(); 1937 VectorIterator<Float> itData(data, axis); 1938 while (!itData.pastEnd()) { 1939 pos2 = itData.pos().getFirst(nAxesSub); 1940 itData.vector() /= tSysSqSum(pos2); 1941 itData.next(); 1942 } 1943 } 1944 } 1945 1946 1947 1948 1949 void SDMath::setCursorSlice (IPosition& start, IPosition& end, Bool doAll, const SDMemTable& in) const 1950 { 1951 const uInt nDim = asap::nAxes; 1952 DebugAssert(nDim==4,AipsError); 1953 // 1954 start.resize(nDim); 1955 end.resize(nDim); 1956 if (doAll) { 1957 start = 0; 1958 end(0) = in.nBeam()-1; 1959 end(1) = in.nIF()-1; 1960 end(2) = in.nPol()-1; 1961 end(3) = in.nChan()-1; 1962 } else { 1963 start(0) = in.getBeam(); 1964 end(0) = start(0); 1965 // 1966 start(1) = in.getIF(); 1967 end(1) = start(1); 1968 // 1969 start(2) = in.getPol(); 1970 end(2) = start(2); 1971 // 1972 start(3) = 0; 1973 end(3) = in.nChan()-1; 1974 } 1975 } 1976 1977 1978 void SDMath::convertWeightString(WeightType& wtType, const String& weightStr, 1979 Bool listType) 1980 { 1981 String tStr(weightStr); 1982 tStr.upcase(); 1983 String msg; 1984 if (tStr.contains(String("NONE"))) { 1985 wtType = NONE; 1986 msg = String("Weighting type selected : None"); 1987 } else if (tStr.contains(String("VAR"))) { 1988 wtType = VAR; 1989 msg = String("Weighting type selected : Variance"); 1990 } else if (tStr.contains(String("TINTSYS"))) { 1991 wtType = TINTSYS; 1992 msg = String("Weighting type selected : Tint&Tsys"); 1993 } else if (tStr.contains(String("TINT"))) { 1994 wtType = TINT; 1995 msg = String("Weighting type selected : Tint"); 1996 } else if (tStr.contains(String("TSYS"))) { 1997 wtType = TSYS; 1998 msg = String("Weighting type selected : Tsys"); 1999 } else { 2000 msg = String("Weighting type selected : None"); 2001 throw(AipsError("Unrecognized weighting type")); 2002 } 2003 // 2004 if (listType) pushLog(msg); 2005 } 2006 2007 2008 void SDMath::convertInterpString(casa::InterpolateArray1D<Double,Float>::InterpolationMethod& type, 2009 const casa::String& interp) 2010 { 2011 String tStr(interp); 2012 tStr.upcase(); 2013 if (tStr.contains(String("NEAR"))) { 2014 type = InterpolateArray1D<Double,Float>::nearestNeighbour; 2015 } else if (tStr.contains(String("LIN"))) { 2016 type = InterpolateArray1D<Double,Float>::linear; 2017 } else if (tStr.contains(String("CUB"))) { 2018 type = InterpolateArray1D<Double,Float>::cubic; 2019 } else if (tStr.contains(String("SPL"))) { 2020 type = InterpolateArray1D<Double,Float>::spline; 2021 } else { 2022 throw(AipsError("Unrecognized interpolation type")); 2023 } 2024 } 2025 2026 void SDMath::putDataInSDC(SDContainer& sc, const Array<Float>& data, 2027 const Array<Bool>& mask) 2028 { 2029 sc.putSpectrum(data); 2030 // 2031 Array<uChar> outflags(data.shape()); 2032 convertArray(outflags,!mask); 2033 sc.putFlags(outflags); 2034 } 2035 2036 Table SDMath::readAsciiFile(const String& fileName) const 2037 { 2038 String formatString; 2039 Table tbl = readAsciiTable (formatString, Table::Memory, fileName, "", "", False); 2040 return tbl; 2041 } 2042 2043 2044 2045 void SDMath::scaleFromAsciiTable(SDMemTable* pTabOut, 2046 const SDMemTable& in, const String& fileName, 2047 const String& col0, const String& col1, 2048 const String& methodStr, Bool doAll, 2049 const Vector<Float>& xOut, Bool doTSys) 2050 { 2051 2052 // Read gain-elevation ascii file data into a Table. 2053 2054 Table geTable = readAsciiFile (fileName); 2055 // 2056 scaleFromTable (pTabOut, in, geTable, col0, col1, methodStr, doAll, xOut, doTSys); 2057 } 2058 2059 void SDMath::scaleFromTable(SDMemTable* pTabOut, const SDMemTable& in, 2060 const Table& tTable, const String& col0, 2061 const String& col1, 2062 const String& methodStr, Bool doAll, 2063 const Vector<Float>& xOut, Bool doTsys) 2064 { 2065 2066 // Get data from Table 2067 2068 ROScalarColumn<Float> geElCol(tTable, col0); 2069 ROScalarColumn<Float> geFacCol(tTable, col1); 2070 Vector<Float> xIn = geElCol.getColumn(); 2071 Vector<Float> yIn = geFacCol.getColumn(); 2072 Vector<Bool> maskIn(xIn.nelements(),True); 2073 2074 // Interpolate (and extrapolate) with desired method 2075 2076 InterpolateArray1D<Double,Float>::InterpolationMethod method; 2077 convertInterpString(method, methodStr); 2078 Int intMethod(method); 2079 // 2080 Vector<Float> yOut; 2081 Vector<Bool> maskOut; 2082 InterpolateArray1D<Float,Float>::interpolate(yOut, maskOut, xOut, 2083 xIn, yIn, maskIn, intMethod, 2084 True, True); 2085 // Apply 2086 2087 scaleByVector(pTabOut, in, doAll, Float(1.0)/yOut, doTsys); 2088 } 2089 2090 2091 void SDMath::scaleByVector(SDMemTable* pTabOut, const SDMemTable& in, 2092 Bool doAll, const Vector<Float>& factor, 2093 Bool doTSys) 2094 { 2095 2096 // Set up data slice 2097 2098 IPosition start, end; 2099 setCursorSlice (start, end, doAll, in); 2100 2101 // Get Tsys column 2102 2103 const Table& tIn = in.table(); 2104 ArrayColumn<Float> tSysCol(tIn, "TSYS"); 2105 Array<Float> tSys; 2106 2107 // Loop over rows and apply correction factor 2108 2109 const uInt axis = asap::ChanAxis; 2110 for (uInt i=0; i < in.nRow(); ++i) { 2111 2112 // Get data 2113 2114 MaskedArray<Float> dataIn(in.rowAsMaskedArray(i)); 2115 MaskedArray<Float> dataIn2 = dataIn(start,end); // reference to dataIn 2116 // 2117 if (doTSys) { 2118 tSysCol.get(i, tSys); 2119 Array<Float> tSys2 = tSys(start,end) * factor[i]; 2120 tSysCol.put(i, tSys); 2121 } 2122 2123 // Apply factor 2124 2125 dataIn2 *= factor[i]; 2126 2127 // Write out 2128 2129 SDContainer sc = in.getSDContainer(i); 2130 putDataInSDC(sc, dataIn.getArray(), dataIn.getMask()); 2131 // 2132 pTabOut->putSDContainer(sc); 2133 } 2134 } 2135 2136 2137 2138 2139 void SDMath::generateDataDescTable (Matrix<uInt>& ddIdx, 2140 SDDataDesc& dDesc, 2141 uInt nIF, 2142 const SDMemTable& in, 2143 const Table& tabIn, 2144 const ROScalarColumn<String>& srcCol, 2145 const ROArrayColumn<uInt>& fqIDCol, 2146 Bool perFreqID) 2147 { 2148 const uInt nRows = tabIn.nrow(); 2149 ddIdx.resize(nRows,nIF); 2150 // 2151 String srcName; 2152 Vector<uInt> freqIDs; 2153 for (uInt iRow=0; iRow<nRows; iRow++) { 2154 srcCol.get(iRow, srcName); 2155 fqIDCol.get(iRow, freqIDs); 2156 const MDirection& dir = in.getDirection(iRow); 2157 // 2158 if (perFreqID) { 2159 2160 // One entry per source/freqID pair 2161 2162 for (uInt iIF=0; iIF<nIF; iIF++) { 2163 ddIdx(iRow,iIF) = dDesc.addEntry(srcName, freqIDs[iIF], dir, 0); 2164 } 2165 } else { 2166 2167 // One entry per source/IF pair. Hang onto the FreqID as well 2168 2169 for (uInt iIF=0; iIF<nIF; iIF++) { 2170 ddIdx(iRow,iIF) = dDesc.addEntry(srcName, iIF, dir, freqIDs[iIF]); 2171 } 2172 } 2173 } 2174 } 2175 2176 2177 2178 2179 2180 MEpoch SDMath::epochFromString(const String& str, MEpoch::Types timeRef) 2181 { 2182 Quantum<Double> qt; 2183 if (MVTime::read(qt,str)) { 2184 MVEpoch mv(qt); 2185 MEpoch me(mv, timeRef); 2186 return me; 2187 } else { 2188 throw(AipsError("Invalid format for Epoch string")); 2189 } 2190 } 2191 2192 2193 String SDMath::formatEpoch(const MEpoch& epoch) const 2194 { 2195 MVTime mvt(epoch.getValue()); 2196 return mvt.string(MVTime::YMD) + String(" (") + epoch.getRefString() + String(")"); 2197 } 2198 2199 2200 2201 void SDMath::generateFrequencyAligners(PtrBlock<FrequencyAligner<Float>* >& a, 2202 const SDDataDesc& dDesc, 2203 const SDMemTable& in, uInt nChan, 2204 MFrequency::Types system, 2205 const MPosition& refPos, 2206 const MEpoch& refEpoch, 2207 Bool perFreqID) 2208 { 2209 for (uInt i=0; i<dDesc.length(); i++) { 2210 uInt ID = dDesc.ID(i); 2211 uInt secID = dDesc.secID(i); 2212 const MDirection& refDir = dDesc.secDir(i); 2213 // 2214 if (perFreqID) { 2215 2216 // One aligner per source/FreqID pair. 2217 2218 SpectralCoordinate sC = in.getSpectralCoordinate(ID); 2219 a[i] = new FrequencyAligner<Float>(sC, nChan, refEpoch, refDir, refPos, system); 2220 } else { 2221 2222 // One aligner per source/IF pair. But we still need the FreqID to 2223 // get the right SC. Hence the messing about with the secondary ID 2224 2225 SpectralCoordinate sC = in.getSpectralCoordinate(secID); 2226 a[i] = new FrequencyAligner<Float>(sC, nChan, refEpoch, refDir, refPos, system); 2227 } 2228 } 2229 } 2230 2231 Vector<uInt> SDMath::getRowRange(const SDMemTable& in) const 2232 { 2233 Vector<uInt> range(2); 2234 range[0] = 0; 2235 range[1] = in.nRow()-1; 2236 return range; 2237 } 2238 2239 2240 Bool SDMath::rowInRange(uInt i, const Vector<uInt>& range) const 2241 { 2242 return (i>=range[0] && i<=range[1]); 2243 } 2244 794 specCol.put(i, specout); 795 } 796 } 797 return out; 798 } -
trunk/src/STMath.h
r804 r805 1 //#--------------------------------------------------------------------------- 2 //# SDMath.h: A collection of single dish mathematical operations 3 //#--------------------------------------------------------------------------- 4 //# Copyright (C) 2004 5 //# ATNF 6 //# 7 //# This program is free software; you can redistribute it and/or modify it 8 //# under the terms of the GNU General Public License as published by the Free 9 //# Software Foundation; either version 2 of the License, or (at your option) 10 //# any later version. 11 //# 12 //# This program is distributed in the hope that it will be useful, but 13 //# WITHOUT ANY WARRANTY; without even the implied warranty of 14 //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 15 //# Public License for more details. 16 //# 17 //# You should have received a copy of the GNU General Public License along 18 //# with this program; if not, write to the Free Software Foundation, Inc., 19 //# 675 Massachusetts Ave, Cambridge, MA 02139, USA. 20 //# 21 //# Correspondence concerning this software should be addressed as follows: 22 //# Internet email: Malte.Marquarding@csiro.au 23 //# Postal address: Malte Marquarding, 24 //# Australia Telescope National Facility, 25 //# P.O. Box 76, 26 //# Epping, NSW, 2121, 27 //# AUSTRALIA 28 //# 29 //# $Id: 30 //#--------------------------------------------------------------------------- 31 #ifndef SDMATH_H 32 #define SDMATH_H 1 // 2 // C++ Interface: STMath 3 // 4 // Description: 5 // 6 // 7 // Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006 8 // 9 // Copyright: See COPYING file that comes with this distribution 10 // 11 // 12 #ifndef ASAPSTMATH_H 13 #define ASAPSTMATH_H 33 14 34 15 #include <string> 35 #include <vector> 16 #include <map> 17 36 18 #include <casa/aips.h> 37 19 #include <casa/Utilities/CountedPtr.h> 38 #include <coordinates/Coordinates/FrequencyAligner.h> 20 #include <casa/BasicSL/String.h> 21 #include <casa/Arrays/Vector.h> 22 #include <scimath/Mathematics/InterpolateArray1D.h> 39 23 24 #include "Scantable.h" 40 25 #include "SDDefs.h" 41 26 #include "SDLog.h" 42 27 43 class casa::Table;44 class casa::MEpoch;45 class casa::MPosition;46 template<class T> class casa::PtrBlock;47 template<class T> class casa::Matrix;48 template<class T> class casa::ROScalarColumn;49 template<class T> class casa::ROArrayColumn;50 51 28 namespace asap { 52 29 53 class SDMemTable; 54 class SDDataDesc; 55 56 class SDMath : private SDLog { 30 /** 31 Mathmatical operations on Scantable objects 57 32 33 @author Malte Marquarding 34 */ 35 class STMath : private SDLog { 58 36 public: 59 60 // Default constructor61 SDMath();62 63 // Copy Constructor (copy semantics)64 SDMath(const SDMath& other);65 66 // Assignment (copy semantics)67 SDMath &operator=(const SDMath& other);68 69 // Destructor70 ~SDMath();71 72 // Binary Table operators. op=ADD, SUB, MUL, DIV, QUOTIENT73 casa::CountedPtr<SDMemTable> binaryOperate(const casa::CountedPtr<SDMemTable>& left,74 const casa::CountedPtr<SDMemTable>& right,75 const casa::String& op, casa::Bool preserve,76 casa::Bool tSys);77 78 // Average in time79 casa::CountedPtr<SDMemTable> average(const std::vector<casa::CountedPtr<SDMemTable> >& in,80 const casa::Vector<casa::Bool>& mask,81 casa::Bool scanAverage,82 const casa::String& weightStr,83 casa::Bool align=casa::False);84 85 // Statistics. If row<0, all rows are done otherwise, just the86 // specified row.87 std::vector<float> statistic(const casa::CountedPtr<SDMemTable>& in,88 const casa::Vector<casa::Bool>& mask,89 const casa::String& which, casa::Int row) const;90 91 // Bin up spectra92 SDMemTable* bin(const SDMemTable& in, casa::Int width);93 37 94 // Resample spectra 95 SDMemTable* resample(const SDMemTable& in, const casa::String& method, 96 casa::Float factor); 38 typedef casa::InterpolateArray1D<casa::Double, 39 casa::Float>::InterpolationMethod imethod; 97 40 98 // Smooth 99 SDMemTable* smooth(const SDMemTable& in, const casa::String& kernel, 100 casa::Float width, casa::Bool doAll); 41 typedef std::map<std::string, imethod> imap; 101 42 102 // Flux conversion between Jansky and Kelvin 103 SDMemTable* convertFlux(const SDMemTable& in, casa::Float D, 104 casa::Float etaAp, 105 casa::Float JyPerK, casa::Bool doAll); 43 STMath(bool insitu=true); 106 44 107 // Gain-elevation correction 108 SDMemTable* gainElevation(const SDMemTable& in, 109 const casa::Vector<casa::Float>& coeffs, 110 const casa::String& fileName, 111 const casa::String& method, 112 casa::Bool doAll); 45 ~STMath(); 113 46 114 // Frequency Alignment 115 SDMemTable* frequencyAlignment(const SDMemTable& in, const casa::String& refTime, 116 const casa::String& method, casa::Bool perFreqID); 47 /** 48 * set the @attr insitu attribute 49 * @param b 50 */ 51 bool insitu() const { return insitu_;}; 52 void setInsitu(bool b) { insitu_ = b; }; 117 53 118 // Opacity correction 119 SDMemTable* opacity(const SDMemTable& in, casa::Float tau, casa::Bool doAll); 54 casa::CountedPtr<Scantable> 55 average( const std::vector<casa::CountedPtr<Scantable> >& in, 56 const casa::Vector<casa::Bool>& mask, 57 const std::string& weight, 58 const std::string& avmode = "", 59 bool alignfreq = false ); 120 60 121 // Simple unary mathematical operations. what=0 (mul) or 1 (add) 122 SDMemTable* unaryOperate(const SDMemTable& in, casa::Float offset,123 casa::Bool doAll, casa::uInt what, casa::Bool tSys);61 casa::CountedPtr<Scantable> 62 unaryOperate( const casa::CountedPtr<Scantable>& in, float val, 63 const std::string& mode, bool tsys=false ); 124 64 125 // Average polarizations. 126 SDMemTable* averagePol(const SDMemTable& in, 127 const casa::Vector<casa::Bool>& mask, 128 const casa::String& wtStr); 65 casa::CountedPtr<Scantable> quotient( const casa::CountedPtr<Scantable>& in, 66 const std::string& mode = "NEAREST", 67 bool preserve = true ); 129 68 130 SDMemTable* frequencySwitch(const SDMemTable& in);131 69 casa::CountedPtr<Scantable> 70 freqSwitch( const casa::CountedPtr<Scantable>& in ); 132 71 133 // Rotate XY phase. Value in degrees. 134 void rotateXYPhase(SDMemTable& in, casa::Float value, casa::Bool doAll); 72 std::vector<float> statistic(const casa::CountedPtr<Scantable>& in, 73 const std::vector<bool>& mask, 74 const std::string& which); 135 75 136 // Rotate Q & U by operating on the raw correlations.Value in degrees. 137 void rotateLinPolPhase(SDMemTable& in, casa::Float value, casa::Bool doAll); 76 casa::CountedPtr<Scantable> bin( const casa::CountedPtr<Scantable>& in, 77 int width=5); 78 casa::CountedPtr<Scantable> 79 resample(const casa::CountedPtr<Scantable>& in, 80 const std::string& method, float width); 138 81 139 140 private: 82 casa::CountedPtr<Scantable> 83 smooth(const casa::CountedPtr<Scantable>& in, const std::string& kernel, 84 float width); 141 85 142 // Weighting type for time averaging 86 casa::CountedPtr<Scantable> 87 gainElevation(const casa::CountedPtr<Scantable>& in, 88 const casa::Vector<casa::Float>& coeffs, 89 const std::string& fileName, 90 const std::string& method); 91 casa::CountedPtr<Scantable> 92 convertFlux(const casa::CountedPtr<Scantable>& in, float d, 93 float etaap, float jyperk); 143 94 144 enum WeightType {NONE=0,VAR,TSYS,TINT,TINTSYS,nWeightTypes}; 95 casa::CountedPtr<Scantable> opacity(const casa::CountedPtr<Scantable>& in, 96 float tau); 145 97 146 // Function to use accumulate data during time averaging 98 /// @todo frequency alignment 147 99 148 void accumulate(casa::Double& timeSum, casa::Double& intSum, 149 casa::Int& nAccum, 150 casa::MaskedArray<casa::Float>& sum, 151 casa::Array<casa::Float>& sumSq, 152 casa::Array<casa::Float>& nPts, 153 casa::Array<casa::Float>& tSysSum, 154 casa::Array<casa::Float>& tSysSqSum, 155 const casa::Array<casa::Float>& tSys, 156 const casa::Array<casa::Float>& nInc, 157 const casa::Vector<casa::Bool>& mask, 158 casa::Double time, casa::Double interval, 159 const std::vector<casa::CountedPtr<SDMemTable> >& in, 160 casa::uInt iTab, casa::uInt iRow, 161 casa::uInt axis, casa::uInt nAxesSub, 162 casa::Bool useMask, WeightType wtType); 100 private: 101 static imethod stringToIMethod(const std::string& in); 102 static WeightType stringToWeight(const std::string& in); 163 103 164 // Work out conversion factor for converting Jy<->K per IF per row and apply 165 void convertBrightnessUnits(SDMemTable* pTabOut, const SDMemTable& in,166 casa::Bool toKelvin, casa::Float sFac, casa::Bool doAll);104 void scaleByVector(casa::Table& in, 105 const casa::Vector<casa::Float>& factor, 106 bool dotsys); 167 107 168 // Convert weight string to enum value 108 void scaleFromAsciiTable(casa::Table& in, const std::string& filename, 109 const std::string& method, 110 const casa::Vector<casa::Float>& xout, 111 bool dotsys); 169 112 170 void convertWeightString(WeightType& wt, const casa::String& weightStr, 171 casa::Bool listType);; 113 void scaleFromTable(casa::Table& in, const casa::Table& table, 114 const std::string& method, 115 const casa::Vector<casa::Float>& xout, bool dotsys); 172 116 173 // Convert interpolation type string 174 // void convertInterpString(casa::Int& type, const casa::String& interp); 175 void convertInterpString(casa::InterpolateArray1D<casa::Double,casa::Float>::InterpolationMethod& method, 176 const casa::String& interp); 117 void convertBrightnessUnits(casa::CountedPtr<Scantable>& in, 118 bool tokelvin, float cfac); 177 119 178 // Scale data with values from an ascii Table 179 void scaleFromAsciiTable(SDMemTable* pTabOut, const SDMemTable& in, 180 const casa::String& fileName, 181 const casa::String& col0, const casa::String& col1, 182 const casa::String& methodStr, casa::Bool doAll, 183 const casa::Vector<casa::Float>& xOut, casa::Bool doTSys); 120 casa::CountedPtr< Scantable > 121 getScantable(const casa::CountedPtr< Scantable >& in, bool droprows); 184 122 185 // Scale data with values from a Table 186 void scaleFromTable(SDMemTable* pTabOut, const SDMemTable& in, const casa::Table& tTable,187 const casa::String& col0, const casa::String& col1,188 const casa::String& methodStr, casa::Bool doAll,189 const casa::Vector<casa::Float>& xOut, casa::Bool doTSys);123 casa::MaskedArray<casa::Float> 124 maskedArray( const casa::Vector<casa::Float>& s, 125 const casa::Vector<casa::uChar>& f ); 126 casa::Vector<casa::uChar> 127 flagsFromMA(const casa::MaskedArray<casa::Float>& ma); 190 128 191 // Scale data and optionally TSys by values in a Vector 192 void scaleByVector(SDMemTable* pTabOut, const SDMemTable& in, 193 casa::Bool doAll, const casa::Vector<casa::Float>& factor, 194 casa::Bool doTSys); 129 bool insitu_; 130 }; 195 131 196 // Convert time String to Epoch 197 casa::MEpoch epochFromString(const casa::String& str, casa::MEpoch::Types timeRef); 198 199 // Function to fill Scan Container when averaging in time 200 201 void fillSDC(SDContainer& sc, const casa::Array<casa::Bool>& mask, 202 const casa::Array<casa::Float>& data, 203 const casa::Array<casa::Float>& tSys, 204 casa::Int scanID, casa::Double timeStamp, 205 casa::Double interval, const casa::String& sourceName, 206 const casa::Vector<casa::uInt>& freqID); 207 208 // Format EPoch 209 casa::String formatEpoch(const casa::MEpoch& epoch) const; 210 211 // Align in Frequency 212 SDMemTable* frequencyAlign(const SDMemTable& in, 213 casa::MFrequency::Types system, 214 const casa::String& timeRef, 215 const casa::String& method, 216 casa::Bool perIF); 217 218 // Generate frequency aligners 219 void generateFrequencyAligners(casa::PtrBlock<casa::FrequencyAligner<casa::Float>* >& a, 220 const SDDataDesc& dDesc, 221 const SDMemTable& in, casa::uInt nChan, 222 casa::MFrequency::Types system, 223 const casa::MPosition& refPos, 224 const casa::MEpoch& refEpoch, 225 casa::Bool perFreqID); 226 227 // Generate data description table(combines source and freqID); 228 void generateDataDescTable(casa::Matrix<casa::uInt>& ddIdx, 229 SDDataDesc& dDesc, 230 casa::uInt nIF, 231 const SDMemTable& in, 232 const casa::Table& tabIn, 233 const casa::ROScalarColumn<casa::String>& srcCol, 234 const casa::ROArrayColumn<casa::uInt>& fqIDCol, 235 casa::Bool perFreqID); 236 237 // Get row range from SDMemTable state 238 casa::Vector<casa::uInt> getRowRange(const SDMemTable& in) const; 239 240 // Is row in the row range ? 241 casa::Bool rowInRange(casa::uInt i, const casa::Vector<casa::uInt>& range) const; 242 243 // Set slice to cursor or all axes 244 void setCursorSlice(casa::IPosition& start, casa::IPosition& end, 245 casa::Bool doAll, const SDMemTable& in) const; 246 247 // Function to normalize data when averaging in time 248 249 void normalize(casa::MaskedArray<casa::Float>& data, 250 const casa::Array<casa::Float>& sumSq, 251 const casa::Array<casa::Float>& tSysSumSq, 252 const casa::Array<casa::Float>& nPts, 253 casa::Double intSum, 254 WeightType wtType, casa::Int axis, casa::Int nAxes); 255 256 // Put the data and mask into the SDContainer 257 void putDataInSDC(SDContainer& sc, const casa::Array<casa::Float>& data, 258 const casa::Array<casa::Bool>& mask); 259 260 // Read ascii file into a Table 261 262 casa::Table readAsciiFile(const casa::String& fileName) const; 263 264 }; // class 265 266 } // namespace 267 132 } 268 133 #endif -
trunk/src/Scantable.cpp
r804 r805 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 1 // 2 // C++ Implementation: Scantable 3 // 4 // Description: 5 // 6 // 7 // Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2005 8 // 9 // Copyright: See COPYING file that comes with this distribution 10 // 11 // 32 12 #include <map> 33 13 … … 35 15 #include <casa/iostream.h> 36 16 #include <casa/iomanip.h> 17 #include <casa/OS/Path.h> 18 #include <casa/OS/File.h> 37 19 #include <casa/Arrays/Array.h> 38 20 #include <casa/Arrays/ArrayMath.h> … … 45 27 #include <casa/BasicSL/Constants.h> 46 28 #include <casa/Quanta/MVAngle.h> 29 #include <casa/Containers/RecordField.h> 47 30 48 31 #include <tables/Tables/TableParse.h> … … 52 35 #include <tables/Tables/ScaColDesc.h> 53 36 #include <tables/Tables/ArrColDesc.h> 37 #include <tables/Tables/TableRow.h> 38 #include <tables/Tables/TableVector.h> 39 #include <tables/Tables/TableIter.h> 54 40 55 41 #include <tables/Tables/ExprNode.h> 56 42 #include <tables/Tables/TableRecord.h> 57 43 #include <measures/Measures/MFrequency.h> 44 #include <measures/Measures/MEpoch.h> 58 45 #include <measures/Measures/MeasTable.h> 46 #include <measures/Measures/MeasRef.h> 47 #include <measures/TableMeasures/TableMeasRefDesc.h> 48 #include <measures/TableMeasures/TableMeasValueDesc.h> 49 #include <measures/TableMeasures/TableMeasDesc.h> 50 #include <measures/TableMeasures/ScalarMeasColumn.h> 59 51 #include <coordinates/Coordinates/CoordinateUtil.h> 60 52 #include <casa/Quanta/MVTime.h> 61 53 #include <casa/Quanta/MVAngle.h> 62 54 63 #include "SDDefs.h" 64 #include "SDContainer.h" 65 #include "MathUtils.h" 66 #include "SDPol.h" 55 #include "Scantable.h" 67 56 #include "SDAttr.h" 68 57 69 #include "SDMemTable.h"70 71 58 using namespace casa; 72 using namespace asap; 73 74 SDMemTable::SDMemTable() : 75 IFSel_(0), 76 beamSel_(0), 77 polSel_(0) 78 { 79 setup(); 59 60 namespace asap { 61 62 Scantable::Scantable(Table::TableType ttype) : 63 type_(ttype), 64 freqTable_(ttype), 65 focusTable_(ttype), 66 weatherTable_(ttype), 67 tcalTable_(ttype), 68 moleculeTable_(ttype) 69 { 70 setupMainTable(); 71 table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_.table()); 72 table_.rwKeywordSet().defineTable("WEATHER", weatherTable_.table()); 73 table_.rwKeywordSet().defineTable("FOCUS", focusTable_.table()); 74 table_.rwKeywordSet().defineTable("TCAL", tcalTable_.table()); 75 table_.rwKeywordSet().defineTable("MOLECULES", moleculeTable_.table()); 76 setupHistoryTable(); 77 setupFitTable(); 78 79 originalTable_ = table_; 80 80 attach(); 81 81 } 82 82 83 SDMemTable::SDMemTable(const std::string& name) : 84 IFSel_(0), 85 beamSel_(0), 86 polSel_(0) 83 Scantable::Scantable(const std::string& name, Table::TableType ttype) : 84 type_(ttype), 85 freqTable_(ttype) 87 86 { 88 87 Table tab(name); … … 92 91 throw(AipsError("Unsupported version of ASAP file.")); 93 92 } 94 table_ = tab.copyToMemoryTable("dummy"); 93 if ( type_ == Table::Memory ) 94 table_ = tab.copyToMemoryTable("dummy"); 95 else 96 table_ = tab; 97 98 originalTable_ = table_; 95 99 attach(); 96 100 } 97 101 98 SDMemTable::SDMemTable(const SDMemTable& other, Bool clear) 99 { 100 table_ = other.table_.copyToMemoryTable(String("dummy")); 101 // clear all rows()102 103 Scantable::Scantable( const Scantable& other, bool clear ) 104 { 105 // with or without data 102 106 if (clear) { 103 table_.removeRow(this->table_.rowNumbers()); 104 IFSel_= 0; 105 beamSel_= 0; 106 polSel_= 0; 107 } else { 108 IFSel_= other.IFSel_; 109 beamSel_= other.beamSel_; 110 polSel_= other.polSel_; 111 } 112 107 table_ = TableCopy::makeEmptyMemoryTable(String("dummy"), other.table_, 108 True); 109 } else { 110 table_ = other.table_.copyToMemoryTable(String("dummy")); 111 } 112 113 originalTable_ = table_; 113 114 attach(); 114 115 } 115 116 116 SDMemTable::SDMemTable(const Table& tab, const std::string& exprs) : 117 IFSel_(0), 118 beamSel_(0), 119 polSel_(0) 120 { 121 Table t = tableCommand(exprs,tab); 122 if (t.nrow() == 0) 123 throw(AipsError("Query unsuccessful.")); 124 table_ = t.copyToMemoryTable("dummy"); 125 attach(); 126 renumber(); 127 } 128 129 SDMemTable::~SDMemTable() 130 { 131 //cerr << "goodbye from SDMemTable @ " << this << endl; 132 } 133 134 SDMemTable SDMemTable::getScan(Int scanID) const 135 { 136 String cond("SELECT * from $1 WHERE SCANID == "); 137 cond += String::toString(scanID); 138 return SDMemTable(table_, cond); 139 } 140 141 SDMemTable &SDMemTable::operator=(const SDMemTable& other) 142 { 143 if (this != &other) { 144 IFSel_= other.IFSel_; 145 beamSel_= other.beamSel_; 146 polSel_= other.polSel_; 147 table_ = other.table_.copyToMemoryTable(String("dummy")); 148 attach(); 149 } 150 return *this; 151 } 152 153 SDMemTable SDMemTable::getSource(const std::string& source) const 154 { 155 String cond("SELECT * from $1 WHERE SRCNAME == "); 156 cond += source; 157 return SDMemTable(table_, cond); 158 } 159 160 void SDMemTable::setup() 117 Scantable::~Scantable() 118 { 119 cout << "~Scantable() " << this << endl; 120 } 121 122 void Scantable::setupMainTable() 161 123 { 162 124 TableDesc td("", "1", TableDesc::Scratch); 163 td.comment() = "A SDMemTable";125 td.comment() = "An ASAP Scantable"; 164 126 td.rwKeywordSet().define("VERSION", Int(version_)); 165 127 128 // n Cycles 129 td.addColumn(ScalarColumnDesc<uInt>("SCANNO")); 130 // new index every nBeam x nIF x nPol 131 td.addColumn(ScalarColumnDesc<uInt>("CYCLENO")); 132 133 td.addColumn(ScalarColumnDesc<uInt>("BEAMNO")); 134 td.addColumn(ScalarColumnDesc<uInt>("IFNO")); 135 td.rwKeywordSet().define("POLTYPE", String("linear")); 136 td.addColumn(ScalarColumnDesc<uInt>("POLNO")); 137 138 td.addColumn(ScalarColumnDesc<uInt>("FREQ_ID")); 139 td.addColumn(ScalarColumnDesc<uInt>("MOLECULE_ID")); 140 // linear, circular, stokes [I Q U V], stokes1 [I Plinear Pangle V] 141 td.addColumn(ScalarColumnDesc<Int>("REFBEAMNO")); 142 166 143 td.addColumn(ScalarColumnDesc<Double>("TIME")); 144 TableMeasRefDesc measRef(MEpoch::UTC); // UTC as default 145 TableMeasValueDesc measVal(td, "TIME"); 146 TableMeasDesc<MEpoch> mepochCol(measVal, measRef); 147 mepochCol.write(td); 148 149 td.addColumn(ScalarColumnDesc<Double>("INTERVAL")); 150 167 151 td.addColumn(ScalarColumnDesc<String>("SRCNAME")); 152 // Type of source (on=0, off=1, other=-1) 153 td.addColumn(ScalarColumnDesc<Int>("SRCTYPE", Int(-1))); 154 td.addColumn(ScalarColumnDesc<String>("FIELDNAME")); 155 156 //The actual Data Vectors 168 157 td.addColumn(ArrayColumnDesc<Float>("SPECTRA")); 169 158 td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA")); 170 159 td.addColumn(ArrayColumnDesc<Float>("TSYS")); 171 td.addColumn(ArrayColumnDesc<Float>("STOKES")); 172 td.addColumn(ScalarColumnDesc<Int>("SCANID")); 173 td.addColumn(ScalarColumnDesc<Double>("INTERVAL")); 174 td.addColumn(ArrayColumnDesc<uInt>("FREQID")); 175 td.addColumn(ArrayColumnDesc<uInt>("RESTFREQID")); 176 td.addColumn(ArrayColumnDesc<Double>("DIRECTION")); 177 td.addColumn(ScalarColumnDesc<String>("FIELDNAME")); 178 td.addColumn(ScalarColumnDesc<String>("TCALTIME")); 179 td.addColumn(ArrayColumnDesc<Float>("TCAL")); 180 td.addColumn(ScalarColumnDesc<Float>("AZIMUTH")); 181 td.addColumn(ScalarColumnDesc<Float>("ELEVATION")); 160 161 td.addColumn(ArrayColumnDesc<Double>("DIRECTION", 162 IPosition(1,2), 163 ColumnDesc::Direct)); 164 TableMeasRefDesc mdirRef(MDirection::J2000); // default 165 TableMeasValueDesc tmvdMDir(td, "DIRECTION"); 166 // the TableMeasDesc gives the column a type 167 TableMeasDesc<MDirection> mdirCol(tmvdMDir, mdirRef); 168 // writing create the measure column 169 mdirCol.write(td); 170 td.addColumn(ScalarColumnDesc<Double>("AZIMUTH")); 171 td.addColumn(ScalarColumnDesc<Double>("ELEVATION")); 182 172 td.addColumn(ScalarColumnDesc<Float>("PARANGLE")); 183 td.addColumn(ScalarColumnDesc<Int>("REFBEAM")); 184 td.addColumn(ArrayColumnDesc<Int>("FITID")); 173 174 td.addColumn(ScalarColumnDesc<uInt>("TCAL_ID")); 175 td.addColumn(ScalarColumnDesc<uInt>("FIT_ID")); 176 177 td.addColumn(ScalarColumnDesc<uInt>("FOCUS_ID")); 178 td.addColumn(ScalarColumnDesc<uInt>("WEATHER_ID")); 179 180 td.rwKeywordSet().define("OBSMODE", String("")); 185 181 186 182 // Now create Table SetUp from the description. 187 183 SetupNewTable aNewTab("dummy", td, Table::New); 188 189 // Bind the Stokes Virtual machine to the STOKES column Because we190 // don't know how many polarizations will be in the data at this191 // point, we must bind the Virtual Engine regardless. The STOKES192 // column won't be accessed if not appropriate (nPol=4)193 SDStokesEngine::registerClass();194 SDStokesEngine stokesEngine(String("STOKES"), String("SPECTRA"));195 aNewTab.bindColumn("STOKES", stokesEngine);196 197 // Create Table198 184 table_ = Table(aNewTab, Table::Memory, 0); 199 // add subtable 200 TableDesc tdf("", "1", TableDesc::Scratch); 201 tdf.addColumn(ArrayColumnDesc<String>("FUNCTIONS")); 202 tdf.addColumn(ArrayColumnDesc<Int>("COMPONENTS")); 203 tdf.addColumn(ArrayColumnDesc<Double>("PARAMETERS")); 204 tdf.addColumn(ArrayColumnDesc<Bool>("PARMASK")); 205 tdf.addColumn(ArrayColumnDesc<String>("FRAMEINFO")); 206 SetupNewTable fittab("fits", tdf, Table::New); 207 Table fitTable(fittab, Table::Memory); 208 table_.rwKeywordSet().defineTable("FITS", fitTable); 209 185 186 originalTable_ = table_; 187 188 } 189 190 void asap::Scantable::setupHistoryTable( ) 191 { 210 192 TableDesc tdh("", "1", TableDesc::Scratch); 211 193 tdh.addColumn(ScalarColumnDesc<String>("ITEM")); 212 SetupNewTable histtab("hist ", tdh, Table::New);194 SetupNewTable histtab("history", tdh, Table::New); 213 195 Table histTable(histtab, Table::Memory); 214 196 table_.rwKeywordSet().defineTable("HISTORY", histTable); 215 197 } 216 198 217 void SDMemTable::attach() 218 { 219 timeCol_.attach(table_, "TIME"); 220 srcnCol_.attach(table_, "SRCNAME"); 221 specCol_.attach(table_, "SPECTRA"); 222 flagsCol_.attach(table_, "FLAGTRA"); 223 tsCol_.attach(table_, "TSYS"); 224 stokesCol_.attach(table_, "STOKES"); 225 scanCol_.attach(table_, "SCANID"); 226 integrCol_.attach(table_, "INTERVAL"); 227 freqidCol_.attach(table_, "FREQID"); 228 restfreqidCol_.attach(table_, "RESTFREQID"); 229 dirCol_.attach(table_, "DIRECTION"); 230 fldnCol_.attach(table_, "FIELDNAME"); 231 tcaltCol_.attach(table_, "TCALTIME"); 232 tcalCol_.attach(table_, "TCAL"); 233 azCol_.attach(table_, "AZIMUTH"); 234 elCol_.attach(table_, "ELEVATION"); 235 paraCol_.attach(table_, "PARANGLE"); 236 rbeamCol_.attach(table_, "REFBEAM"); 237 fitCol_.attach(table_,"FITID"); 238 } 239 240 241 std::string SDMemTable::getSourceName(Int whichRow) const 242 { 243 String name; 244 srcnCol_.get(whichRow, name); 245 return name; 246 } 247 248 float SDMemTable::getElevation(Int whichRow) const 249 { 250 float elevation; 251 elCol_.get(whichRow, elevation); 252 return elevation; 253 } 254 255 float SDMemTable::getAzimuth(Int whichRow) const 256 { 257 float azimuth; 258 azCol_.get(whichRow, azimuth); 259 return azimuth; 260 } 261 262 float SDMemTable::getParAngle(Int whichRow) const 263 { 264 float parangle; 265 paraCol_.get(whichRow, parangle); 266 return parangle; 267 } 268 269 std::string SDMemTable::getTime(Int whichRow, Bool showDate) const 270 { 271 Double tm; 272 if (whichRow > -1) { 273 timeCol_.get(whichRow, tm); 274 } else { 275 table_.keywordSet().get("UTC",tm); 276 } 277 MVTime mvt(tm); 278 if (showDate) 279 mvt.setFormat(MVTime::YMD); 280 else 281 mvt.setFormat(MVTime::TIME); 282 ostringstream oss; 283 oss << mvt; 284 return String(oss); 285 } 286 287 double SDMemTable::getInterval(Int whichRow) const 288 { 289 Double intval; 290 integrCol_.get(whichRow, intval); 291 return intval; 292 } 293 294 bool SDMemTable::setIF(Int whichIF) 295 { 296 if ( whichIF >= 0 && whichIF < nIF()) { 297 IFSel_ = whichIF; 298 return true; 299 } 300 return false; 301 } 302 303 bool SDMemTable::setBeam(Int whichBeam) 304 { 305 if ( whichBeam >= 0 && whichBeam < nBeam()) { 306 beamSel_ = whichBeam; 307 return true; 308 } 309 return false; 310 } 311 312 bool SDMemTable::setPol(Int whichPol) 313 { 314 if ( whichPol >= 0 && whichPol < nPol()) { 315 polSel_ = whichPol; 316 return true; 317 } 318 return false; 319 } 320 321 void SDMemTable::resetCursor() 322 { 323 polSel_ = 0; 324 IFSel_ = 0; 325 beamSel_ = 0; 326 } 327 328 std::vector<bool> SDMemTable::getMask(Int whichRow) const 329 { 330 331 std::vector<bool> mask; 332 333 Array<uChar> arr; 334 flagsCol_.get(whichRow, arr); 335 336 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr); 337 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam 338 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0); 339 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF 340 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1); 341 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol 342 343 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) { 344 bool out =!static_cast<bool>(*i); 345 mask.push_back(out); 346 } 347 return mask; 348 } 349 350 351 352 std::vector<float> SDMemTable::getSpectrum(Int whichRow) const 353 { 354 Array<Float> arr; 355 specCol_.get(whichRow, arr); 356 return getFloatSpectrum(arr); 357 } 358 359 360 int SDMemTable::nStokes() const 361 { 362 return stokesCol_.shape(0).nelements(); // All rows same shape 363 } 364 365 366 std::vector<float> SDMemTable::getStokesSpectrum(Int whichRow, 367 Bool doPol) const 368 // 369 // Gets one STokes parameter depending on cursor polSel location 370 // doPol=False : I,Q,U,V 371 // doPol=True : I,P,PA,V ; P = sqrt(Q**2+U**2), PA = 0.5*atan2(Q,U) 372 // 373 { 374 AlwaysAssert(asap::nAxes==4,AipsError); 375 if (nPol()!=1 && nPol()!=2 && nPol()!=4) { 376 throw (AipsError("You must have 1,2 or 4 polarizations to get the Stokes parameters")); 377 } 378 379 // For full conversion we are only supporting linears at the moment 380 381 if (nPol() > 2) { 382 String antName; 383 table_.keywordSet().get("AntennaName", antName); 384 Instrument inst = SDAttr::convertInstrument (antName, True); 385 SDAttr sdAtt; 386 if (sdAtt.feedPolType(inst) != LINEAR) { 387 throw(AipsError("Only linear polarizations are supported")); 388 } 389 } 390 391 Array<Float> arr; 392 stokesCol_.get(whichRow, arr); 393 394 if (doPol && (polSel_==1 || polSel_==2)) { // Q,U --> P, P.A. 395 396 // Set current cursor location 397 398 const IPosition& shape = arr.shape(); 399 IPosition start, end; 400 getCursorSlice(start, end, shape); 401 402 // Get Q and U slices 403 404 Array<Float> Q = SDPolUtil::getStokesSlice(arr,start,end,"Q"); 405 Array<Float> U = SDPolUtil::getStokesSlice(arr,start,end,"U"); 406 407 // Compute output 408 409 Array<Float> out; 410 if (polSel_==1) { // P 411 out = SDPolUtil::polarizedIntensity(Q,U); 412 } else if (polSel_==2) { // P.A. 413 out = SDPolUtil::positionAngle(Q,U); 414 } 415 416 // Copy to output 417 418 IPosition vecShape(1,shape(asap::ChanAxis)); 419 Vector<Float> outV = out.reform(vecShape); 420 std::vector<float> stlout; 421 outV.tovector(stlout); 422 return stlout; 423 424 } else { 425 // Selects at the cursor location 426 return getFloatSpectrum(arr); 427 } 428 } 429 430 std::string SDMemTable::getPolarizationLabel(Bool linear, Bool stokes, 431 Bool linPol, Int polIdx) const 432 { 433 uInt idx = polSel_; 434 if (polIdx >=0) idx = polIdx; 435 return SDPolUtil::polarizationLabel(idx, linear, stokes, linPol); 436 } 437 438 439 440 std::vector<float> SDMemTable::stokesToPolSpectrum(Int whichRow, 441 Bool toLinear, 442 Int polIdx) const 443 // 444 // polIdx 445 // 0:3 -> RR,LL,Real(RL),Imag(RL) 446 // XX,YY,Real(XY),Image(XY) 447 // 448 // Gets only 449 // RR = I + V 450 // LL = I - V 451 // at the moment 452 // 453 { 454 AlwaysAssert(asap::nAxes==4,AipsError); 455 if (nStokes()!=4) { 456 throw (AipsError("You must have 4 Stokes to convert to linear or circular")); 457 } 458 // 459 Array<Float> arr, out; 460 stokesCol_.get(whichRow, arr); 461 462 // Set current cursor location 463 464 const IPosition& shape = arr.shape(); 465 IPosition start, end; 466 getCursorSlice(start, end, shape); 467 468 // Get the slice 469 470 if (toLinear) { 471 throw(AipsError("Conversion to linears not yet supported")); 472 } else { 473 uInt selection = polSel_; 474 if (polIdx > -1) selection = polIdx; 475 Bool doRR = (selection==0); 476 if (selection>1) { 477 throw(AipsError("Only conversion to RR & LL is currently supported")); 478 } 479 480 // Get I and V slices 481 Array<Float> I = SDPolUtil::getStokesSlice(arr,start,end,"I"); 482 Array<Float> V = SDPolUtil::getStokesSlice(arr,start,end,"V"); 483 484 // Compute output 485 out = SDPolUtil::circularPolarizationFromStokes(I, V, doRR); 486 } 487 488 // Copy to output 489 IPosition vecShape(1,shape(asap::ChanAxis)); 490 Vector<Float> outV = out.reform(vecShape); 491 std::vector<float> stlout; 492 outV.tovector(stlout); 493 // 494 return stlout; 495 } 496 497 498 499 500 Array<Float> SDMemTable::getStokesSpectrum(Int whichRow, Int iBeam, Int iIF) const 501 { 502 503 // Get data 504 505 Array<Float> arr; 506 stokesCol_.get(whichRow, arr); 507 508 // Set current cursor location and overwrite polarization axis 509 510 const IPosition& shape = arr.shape(); 511 IPosition start(shape.nelements(),0); 512 IPosition end(shape-1); 513 if (iBeam!=-1) { 514 start(asap::BeamAxis) = iBeam; 515 end(asap::BeamAxis) = iBeam; 516 } 517 if (iIF!=-1) { 518 start(asap::IFAxis) = iIF; 519 end(asap::IFAxis) = iIF; 520 } 521 522 // Get slice 523 524 return arr(start,end); 525 } 526 527 528 std::vector<string> SDMemTable::getCoordInfo() const 529 { 530 String un; 531 Table t = table_.keywordSet().asTable("FREQUENCIES"); 532 String sunit; 533 t.keywordSet().get("UNIT",sunit); 534 String dpl; 535 t.keywordSet().get("DOPPLER",dpl); 536 if (dpl == "") dpl = "RADIO"; 537 String rfrm; 538 t.keywordSet().get("REFFRAME",rfrm); 539 std::vector<string> inf; 540 inf.push_back(sunit); 541 inf.push_back(rfrm); 542 inf.push_back(dpl); 543 t.keywordSet().get("BASEREFFRAME",rfrm); 544 inf.push_back(rfrm); 545 return inf; 546 } 547 548 void SDMemTable::setCoordInfo(std::vector<string> theinfo) 549 { 550 std::vector<string>::iterator it; 551 String un,rfrm, brfrm,dpl; 552 un = theinfo[0]; // Abcissa unit 553 rfrm = theinfo[1]; // Active (or conversion) frame 554 dpl = theinfo[2]; // Doppler 555 brfrm = theinfo[3]; // Base frame 556 Table t = table_.rwKeywordSet().asTable("FREQUENCIES"); 557 558 Vector<Double> rstf; 559 t.keywordSet().get("RESTFREQS",rstf); 560 561 Bool canDo = True; 562 Unit u1("km/s");Unit u2("Hz"); 563 if (Unit(un) == u1) { 564 Vector<Double> rstf; 565 t.keywordSet().get("RESTFREQS",rstf); 566 if (rstf.nelements() == 0) { 567 throw(AipsError("Can't set unit to km/s if no restfrequencies are specified")); 568 } 569 } else if (Unit(un) != u2 && un != "") { 570 throw(AipsError("Unit not conformant with Spectral Coordinates")); 571 } 572 t.rwKeywordSet().define("UNIT", un); 573 574 MFrequency::Types mdr; 575 if (!MFrequency::getType(mdr, rfrm)) { 576 577 Int a,b;const uInt* c; 578 const String* valid = MFrequency::allMyTypes(a, b, c); 579 String pfix = "Please specify a legal frame type. Types are\n"; 580 throw(AipsError(pfix+(*valid))); 581 } else { 582 t.rwKeywordSet().define("REFFRAME",rfrm); 583 } 584 585 MDoppler::Types dtype; 586 dpl.upcase(); 587 if (!MDoppler::getType(dtype, dpl)) { 588 throw(AipsError("Doppler type unknown")); 589 } else { 590 t.rwKeywordSet().define("DOPPLER",dpl); 591 } 592 593 if (!MFrequency::getType(mdr, brfrm)) { 594 Int a,b;const uInt* c; 595 const String* valid = MFrequency::allMyTypes(a, b, c); 596 String pfix = "Please specify a legal frame type. Types are\n"; 597 throw(AipsError(pfix+(*valid))); 598 } else { 599 t.rwKeywordSet().define("BASEREFFRAME",brfrm); 600 } 601 } 602 603 604 std::vector<double> SDMemTable::getAbcissa(Int whichRow) const 605 { 606 std::vector<double> abc(nChan()); 607 608 // Get header units keyword 609 Table t = table_.keywordSet().asTable("FREQUENCIES"); 610 String sunit; 611 t.keywordSet().get("UNIT",sunit); 612 if (sunit == "") sunit = "pixel"; 613 Unit u(sunit); 614 615 // Easy if just wanting pixels 616 if (sunit==String("pixel")) { 617 // assume channels/pixels 618 std::vector<double>::iterator it; 619 uInt i=0; 620 for (it = abc.begin(); it != abc.end(); ++it) { 621 (*it) = Double(i++); 622 } 623 return abc; 624 } 625 626 // Continue with km/s or Hz. Get FreqIDs 627 Vector<uInt> freqIDs; 628 freqidCol_.get(whichRow, freqIDs); 629 uInt freqID = freqIDs(IFSel_); 630 restfreqidCol_.get(whichRow, freqIDs); 631 uInt restFreqID = freqIDs(IFSel_); 632 633 // Get SpectralCoordinate, set reference frame conversion, 634 // velocity conversion, and rest freq state 635 636 SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow); 637 Vector<Double> pixel(nChan()); 638 indgen(pixel); 639 640 if (u == Unit("km/s")) { 641 Vector<Double> world; 642 spc.pixelToVelocity(world,pixel); 643 std::vector<double>::iterator it; 644 uInt i = 0; 645 for (it = abc.begin(); it != abc.end(); ++it) { 646 (*it) = world[i]; 647 i++; 648 } 649 } else if (u == Unit("Hz")) { 650 651 // Set world axis units 652 Vector<String> wau(1); wau = u.getName(); 653 spc.setWorldAxisUnits(wau); 654 655 std::vector<double>::iterator it; 656 Double tmp; 657 uInt i = 0; 658 for (it = abc.begin(); it != abc.end(); ++it) { 659 spc.toWorld(tmp,pixel[i]); 660 (*it) = tmp; 661 i++; 662 } 663 } 664 return abc; 665 } 666 667 std::string SDMemTable::getAbcissaString(Int whichRow) const 668 { 669 Table t = table_.keywordSet().asTable("FREQUENCIES"); 670 671 String sunit; 672 t.keywordSet().get("UNIT",sunit); 673 if (sunit == "") sunit = "pixel"; 674 Unit u(sunit); 675 676 Vector<uInt> freqIDs; 677 freqidCol_.get(whichRow, freqIDs); 678 uInt freqID = freqIDs(IFSel_); 679 restfreqidCol_.get(whichRow, freqIDs); 680 uInt restFreqID = freqIDs(IFSel_); 681 682 // Get SpectralCoordinate, with frame, velocity, rest freq state set 683 SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow); 684 685 String s = "Channel"; 686 if (u == Unit("km/s")) { 687 s = CoordinateUtil::axisLabel(spc,0,True,True,True); 688 } else if (u == Unit("Hz")) { 689 Vector<String> wau(1);wau = u.getName(); 690 spc.setWorldAxisUnits(wau); 691 692 s = CoordinateUtil::axisLabel(spc,0,True,True,False); 693 } 694 return s; 695 } 696 697 void SDMemTable::setSpectrum(std::vector<float> spectrum, int whichRow) 698 { 699 Array<Float> arr; 700 specCol_.get(whichRow, arr); 701 if (spectrum.size() != arr.shape()(asap::ChanAxis)) { 702 throw(AipsError("Attempting to set spectrum with incorrect length.")); 703 } 704 705 // Setup accessors 706 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr); 707 aa0.reset(aa0.begin(uInt(beamSel_))); // Beam selection 708 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0); 709 aa1.reset(aa1.begin(uInt(IFSel_))); // IF selection 710 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1); 711 aa2.reset(aa2.begin(uInt(polSel_))); // Pol selection 712 713 // Iterate 714 std::vector<float>::iterator it = spectrum.begin(); 715 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) { 716 (*i) = Float(*it); 717 it++; 718 } 719 specCol_.put(whichRow, arr); 720 } 721 722 void SDMemTable::getSpectrum(Vector<Float>& spectrum, Int whichRow) const 723 { 724 Array<Float> arr; 725 specCol_.get(whichRow, arr); 726 727 // Iterate and extract 728 spectrum.resize(arr.shape()(3)); 729 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr); 730 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam 731 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0); 732 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF 733 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1); 734 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol 735 736 ArrayAccessor<Float, Axis<asap::BeamAxis> > va(spectrum); 737 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) { 738 (*va) = (*i); 739 va++; 740 } 741 } 742 743 744 /* 745 void SDMemTable::getMask(Vector<Bool>& mask, Int whichRow) const { 746 Array<uChar> arr; 747 flagsCol_.get(whichRow, arr); 748 mask.resize(arr.shape()(3)); 749 750 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr); 751 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam 752 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0); 753 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF 754 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1); 755 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol 756 757 Bool useUserMask = ( chanMask_.size() == arr.shape()(3) ); 758 759 ArrayAccessor<Bool, Axis<asap::BeamAxis> > va(mask); 760 std::vector<bool> tmp; 761 tmp = chanMask_; // WHY the fxxx do I have to make a copy here. The 762 // iterator should work on chanMask_?? 763 std::vector<bool>::iterator miter; 764 miter = tmp.begin(); 765 766 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) { 767 bool out =!static_cast<bool>(*i); 768 if (useUserMask) { 769 out = out && (*miter); 770 miter++; 771 } 772 (*va) = out; 773 va++; 774 } 775 } 776 */ 777 778 MaskedArray<Float> SDMemTable::rowAsMaskedArray(uInt whichRow, 779 Bool toStokes) const 780 { 781 // Get flags 782 Array<uChar> farr; 783 flagsCol_.get(whichRow, farr); 784 785 // Get data and convert mask to Bool 786 Array<Float> arr; 787 Array<Bool> mask; 788 if (toStokes) { 789 stokesCol_.get(whichRow, arr); 790 791 Array<Bool> tMask(farr.shape()); 792 convertArray(tMask, farr); 793 mask = SDPolUtil::stokesData (tMask, True); 794 } else { 795 specCol_.get(whichRow, arr); 796 mask.resize(farr.shape()); 797 convertArray(mask, farr); 798 } 799 800 return MaskedArray<Float>(arr,!mask); 801 } 802 803 Float SDMemTable::getTsys(Int whichRow) const 804 { 805 Array<Float> arr; 806 tsCol_.get(whichRow, arr); 807 Float out; 808 809 IPosition ip(arr.shape()); 810 ip(asap::BeamAxis) = beamSel_; 811 ip(asap::IFAxis) = IFSel_; 812 ip(asap::PolAxis) = polSel_; 813 ip(asap::ChanAxis)=0; // First channel only 814 815 out = arr(ip); 816 return out; 817 } 818 819 MDirection SDMemTable::getDirection(Int whichRow, Bool refBeam) const 820 { 821 MDirection::Types mdr = getDirectionReference(); 822 Array<Double> posit; 823 dirCol_.get(whichRow,posit); 824 Vector<Double> wpos(2); 825 Int rb; 826 rbeamCol_.get(whichRow,rb); 827 wpos[0] = posit(IPosition(2,beamSel_,0)); 828 wpos[1] = posit(IPosition(2,beamSel_,1)); 829 if (refBeam && rb > -1) { // use refBeam instead if it exists 830 wpos[0] = posit(IPosition(2,rb,0)); 831 wpos[1] = posit(IPosition(2,rb,1)); 832 } 833 834 Quantum<Double> lon(wpos[0],Unit(String("rad"))); 835 Quantum<Double> lat(wpos[1],Unit(String("rad"))); 836 return MDirection(lon, lat, mdr); 837 } 838 839 MEpoch SDMemTable::getEpoch(Int whichRow) const 840 { 841 MEpoch::Types met = getTimeReference(); 842 843 Double obstime; 844 timeCol_.get(whichRow,obstime); 845 MVEpoch tm2(Quantum<Double>(obstime, Unit(String("d")))); 846 return MEpoch(tm2, met); 847 } 848 849 MPosition SDMemTable::getAntennaPosition () const 850 { 851 Vector<Double> antpos; 852 table_.keywordSet().get("AntennaPosition", antpos); 853 MVPosition mvpos(antpos(0),antpos(1),antpos(2)); 854 return MPosition(mvpos); 855 } 856 857 858 SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID) const 859 { 860 861 Table t = table_.keywordSet().asTable("FREQUENCIES"); 862 if (freqID> t.nrow() ) { 863 throw(AipsError("SDMemTable::getSpectralCoordinate - freqID out of range")); 864 } 865 866 Double rp,rv,inc; 867 String rf; 868 ROScalarColumn<Double> rpc(t, "REFPIX"); 869 ROScalarColumn<Double> rvc(t, "REFVAL"); 870 ROScalarColumn<Double> incc(t, "INCREMENT"); 871 t.keywordSet().get("BASEREFFRAME",rf); 872 873 // Create SpectralCoordinate (units Hz) 874 MFrequency::Types mft; 875 if (!MFrequency::getType(mft, rf)) { 876 ostringstream oss; 877 pushLog("WARNING: Frequency type unknown assuming TOPO"); 878 mft = MFrequency::TOPO; 879 } 880 rpc.get(freqID, rp); 881 rvc.get(freqID, rv); 882 incc.get(freqID, inc); 883 884 SpectralCoordinate spec(mft,rv,inc,rp); 885 return spec; 886 } 887 888 889 SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID, 890 uInt restFreqID, 891 uInt whichRow) const 892 { 893 894 // Create basic SC 895 SpectralCoordinate spec = getSpectralCoordinate (freqID); 896 897 Table t = table_.keywordSet().asTable("FREQUENCIES"); 898 899 // Set rest frequency 900 Vector<Double> restFreqIDs; 901 t.keywordSet().get("RESTFREQS",restFreqIDs); 902 if (restFreqID < restFreqIDs.nelements()) { // Should always be true 903 spec.setRestFrequency(restFreqIDs(restFreqID),True); 904 } 905 906 // Set up frame conversion layer 907 String frm; 908 t.keywordSet().get("REFFRAME",frm); 909 if (frm == "") frm = "TOPO"; 910 MFrequency::Types mtype; 911 if (!MFrequency::getType(mtype, frm)) { 912 // Should never happen 913 pushLog("WARNING: Frequency type unknown assuming TOPO"); 914 mtype = MFrequency::TOPO; 915 } 916 917 // Set reference frame conversion (requires row) 918 MDirection dir = getDirection(whichRow); 919 MEpoch epoch = getEpoch(whichRow); 920 MPosition pos = getAntennaPosition(); 921 922 if (!spec.setReferenceConversion(mtype,epoch,pos,dir)) { 923 throw(AipsError("Couldn't convert frequency frame.")); 924 } 925 926 // Now velocity conversion if appropriate 927 String unitStr; 928 t.keywordSet().get("UNIT",unitStr); 929 930 String dpl; 931 t.keywordSet().get("DOPPLER",dpl); 932 MDoppler::Types dtype; 933 MDoppler::getType(dtype, dpl); 934 935 // Only set velocity unit if non-blank and non-Hz 936 if (!unitStr.empty()) { 937 Unit unitU(unitStr); 938 if (unitU==Unit("Hz")) { 939 } else { 940 spec.setVelocity(unitStr, dtype); 941 } 942 } 943 944 return spec; 945 } 946 947 948 Bool SDMemTable::setCoordinate(const SpectralCoordinate& speccord, 949 uInt freqID) { 950 Table t = table_.rwKeywordSet().asTable("FREQUENCIES"); 951 if (freqID > t.nrow() ) { 952 throw(AipsError("SDMemTable::setCoordinate - coord no out of range")); 953 } 954 ScalarColumn<Double> rpc(t, "REFPIX"); 955 ScalarColumn<Double> rvc(t, "REFVAL"); 956 ScalarColumn<Double> incc(t, "INCREMENT"); 957 958 rpc.put(freqID, speccord.referencePixel()[0]); 959 rvc.put(freqID, speccord.referenceValue()[0]); 960 incc.put(freqID, speccord.increment()[0]); 961 962 return True; 963 } 964 965 Int SDMemTable::nCoordinates() const 966 { 967 return table_.keywordSet().asTable("FREQUENCIES").nrow(); 968 } 969 970 971 std::vector<double> SDMemTable::getRestFreqs() const 972 { 973 Table t = table_.keywordSet().asTable("FREQUENCIES"); 974 Vector<Double> tvec; 975 t.keywordSet().get("RESTFREQS",tvec); 976 std::vector<double> stlout; 977 tvec.tovector(stlout); 978 return stlout; 979 } 980 981 bool SDMemTable::putSDFreqTable(const SDFrequencyTable& sdft) 199 void Scantable::setupFitTable() 982 200 { 983 201 TableDesc td("", "1", TableDesc::Scratch); 984 td.addColumn(ScalarColumnDesc<Double>("REFPIX")); 985 td.addColumn(ScalarColumnDesc<Double>("REFVAL")); 986 td.addColumn(ScalarColumnDesc<Double>("INCREMENT")); 987 SetupNewTable aNewTab("freqs", td, Table::New); 988 Table aTable (aNewTab, Table::Memory, sdft.length()); 989 ScalarColumn<Double> sc0(aTable, "REFPIX"); 990 ScalarColumn<Double> sc1(aTable, "REFVAL"); 991 ScalarColumn<Double> sc2(aTable, "INCREMENT"); 992 for (uInt i=0; i < sdft.length(); ++i) { 993 sc0.put(i,sdft.referencePixel(i)); 994 sc1.put(i,sdft.referenceValue(i)); 995 sc2.put(i,sdft.increment(i)); 996 } 997 String rf = sdft.refFrame(); 998 if (rf.contains("TOPO")) rf = "TOPO"; 999 String brf = sdft.baseRefFrame(); 1000 if (brf.contains("TOPO")) brf = "TOPO"; 1001 1002 aTable.rwKeywordSet().define("BASEREFFRAME", brf); 1003 aTable.rwKeywordSet().define("REFFRAME", rf); 1004 aTable.rwKeywordSet().define("EQUINOX", sdft.equinox()); 1005 aTable.rwKeywordSet().define("UNIT", sdft.unit()); 1006 aTable.rwKeywordSet().define("DOPPLER", String("RADIO")); 1007 Vector<Double> rfvec; 1008 String rfunit; 1009 sdft.restFrequencies(rfvec,rfunit); 1010 Quantum<Vector<Double> > q(rfvec, rfunit); 1011 rfvec.resize(); 1012 rfvec = q.getValue("Hz"); 1013 aTable.rwKeywordSet().define("RESTFREQS", rfvec); 1014 table_.rwKeywordSet().defineTable("FREQUENCIES", aTable); 1015 return true; 1016 } 1017 1018 bool SDMemTable::putSDFitTable(const SDFitTable& sdft) 1019 { 1020 TableDesc td("", "1", TableDesc::Scratch); 202 td.addColumn(ScalarColumnDesc<uInt>("FIT_ID")); 1021 203 td.addColumn(ArrayColumnDesc<String>("FUNCTIONS")); 1022 204 td.addColumn(ArrayColumnDesc<Int>("COMPONENTS")); … … 1026 208 SetupNewTable aNewTab("fits", td, Table::New); 1027 209 Table aTable(aNewTab, Table::Memory); 1028 ArrayColumn<String> sc0(aTable, "FUNCTIONS");1029 ArrayColumn<Int> sc1(aTable, "COMPONENTS");1030 ArrayColumn<Double> sc2(aTable, "PARAMETERS");1031 ArrayColumn<Bool> sc3(aTable, "PARMASK");1032 ArrayColumn<String> sc4(aTable, "FRAMEINFO");1033 for (uInt i; i<sdft.length(); ++i) {1034 const Vector<Double>& parms = sdft.getParameters(i);1035 const Vector<Bool>& parmask = sdft.getParameterMask(i);1036 const Vector<String>& funcs = sdft.getFunctions(i);1037 const Vector<Int>& comps = sdft.getComponents(i);1038 const Vector<String>& finfo = sdft.getFrameInfo(i);1039 sc0.put(i,funcs);1040 sc1.put(i,comps);1041 sc3.put(i,parmask);1042 sc2.put(i,parms);1043 sc4.put(i,finfo);1044 }1045 210 table_.rwKeywordSet().defineTable("FITS", aTable); 1046 return true; 1047 } 1048 1049 SDFitTable SDMemTable::getSDFitTable() const 1050 { 1051 const Table& t = table_.keywordSet().asTable("FITS"); 1052 Vector<Double> parms; 1053 Vector<Bool> parmask; 1054 Vector<String> funcs; 1055 Vector<String> finfo; 1056 Vector<Int> comps; 1057 ROArrayColumn<Double> parmsCol(t, "PARAMETERS"); 1058 ROArrayColumn<Bool> parmaskCol(t, "PARMASK"); 1059 ROArrayColumn<Int> compsCol(t, "COMPONENTS"); 1060 ROArrayColumn<String> funcsCol(t, "FUNCTIONS"); 1061 ROArrayColumn<String> finfoCol(t, "FRAMEINFO"); 1062 uInt n = t.nrow(); 1063 SDFitTable sdft; 1064 for (uInt i=0; i<n; ++i) { 1065 parmaskCol.get(i, parmask); 1066 parmsCol.get(i, parms); 1067 funcsCol.get(i, funcs); 1068 compsCol.get(i, comps); 1069 finfoCol.get(i, finfo); 1070 sdft.addFit(parms, parmask, funcs, comps, finfo); 1071 } 1072 return sdft; 1073 } 1074 1075 SDFitTable SDMemTable::getSDFitTable(uInt whichRow) const { 1076 const Table& t = table_.keywordSet().asTable("FITS"); 1077 if (t.nrow() == 0 || whichRow >= t.nrow()) return SDFitTable(); 1078 Array<Int> fitid; 1079 fitCol_.get(whichRow, fitid); 1080 if (fitid.nelements() == 0) return SDFitTable(); 1081 1082 IPosition shp = fitid.shape(); 1083 IPosition start(4, beamSel_, IFSel_, polSel_,0); 1084 IPosition end(4, beamSel_, IFSel_, polSel_, shp[3]-1); 1085 1086 // reform the output array slice to be of dim=1 1087 Vector<Int> tmp = (fitid(start, end)).reform(IPosition(1,shp[3])); 1088 1089 Vector<Double> parms; 1090 Vector<Bool> parmask; 1091 Vector<String> funcs; 1092 Vector<String> finfo; 1093 Vector<Int> comps; 1094 ROArrayColumn<Double> parmsCol(t, "PARAMETERS"); 1095 ROArrayColumn<Bool> parmaskCol(t, "PARMASK"); 1096 ROArrayColumn<Int> compsCol(t, "COMPONENTS"); 1097 ROArrayColumn<String> funcsCol(t, "FUNCTIONS"); 1098 ROArrayColumn<String> finfoCol(t, "FRAMEINFO"); 1099 SDFitTable sdft; 1100 Int k=-1; 1101 for (uInt i=0; i< tmp.nelements(); ++i) { 1102 k = tmp[i]; 1103 if ( k > -1 && k < t.nrow() ) { 1104 parms.resize(); 1105 parmsCol.get(k, parms); 1106 parmask.resize(); 1107 parmaskCol.get(k, parmask); 1108 funcs.resize(); 1109 funcsCol.get(k, funcs); 1110 comps.resize(); 1111 compsCol.get(k, comps); 1112 finfo.resize(); 1113 finfoCol.get(k, finfo); 1114 sdft.addFit(parms, parmask, funcs, comps, finfo); 1115 } 1116 } 1117 return sdft; 1118 } 1119 1120 void SDMemTable::addFit(uInt whichRow, 1121 const Vector<Double>& p, const Vector<Bool>& m, 1122 const Vector<String>& f, const Vector<Int>& c) 1123 { 1124 if (whichRow >= nRow()) { 1125 throw(AipsError("Specified row out of range")); 1126 } 1127 Table t = table_.keywordSet().asTable("FITS"); 1128 uInt nrow = t.nrow(); 1129 t.addRow(); 1130 ArrayColumn<Double> parmsCol(t, "PARAMETERS"); 1131 ArrayColumn<Bool> parmaskCol(t, "PARMASK"); 1132 ArrayColumn<Int> compsCol(t, "COMPONENTS"); 1133 ArrayColumn<String> funcsCol(t, "FUNCTIONS"); 1134 ArrayColumn<String> finfoCol(t, "FRAMEINFO"); 1135 parmsCol.put(nrow, p); 1136 parmaskCol.put(nrow, m); 1137 compsCol.put(nrow, c); 1138 funcsCol.put(nrow, f); 1139 Vector<String> fi = mathutil::toVectorString(getCoordInfo()); 1140 finfoCol.put(nrow, fi); 1141 1142 Array<Int> fitarr; 1143 fitCol_.get(whichRow, fitarr); 1144 1145 Array<Int> newarr; // The new Array containing the fitid 1146 Int pos =-1; // The fitid position in the array 1147 if ( fitarr.nelements() == 0 ) { // no fits at all in this row 1148 Array<Int> arr(IPosition(4,nBeam(),nIF(),nPol(),1)); 1149 arr = -1; 1150 newarr.reference(arr); 1151 pos = 0; 1152 } else { 1153 IPosition shp = fitarr.shape(); 1154 IPosition start(4, beamSel_, IFSel_, polSel_,0); 1155 IPosition end(4, beamSel_, IFSel_, polSel_, shp[3]-1); 1156 // reform the output array slice to be of dim=1 1157 Array<Int> tmp = (fitarr(start, end)).reform(IPosition(1,shp[3])); 1158 const Vector<Int>& fits = tmp; 1159 VectorSTLIterator<Int> it(fits); 1160 Int i = 0; 1161 while (it != fits.end()) { 1162 if (*it == -1) { 1163 pos = i; 1164 break; 1165 } 1166 ++i; 1167 ++it; 1168 }; 1169 } 1170 if (pos == -1) { 1171 mathutil::extendLastArrayAxis(newarr,fitarr, -1); 1172 pos = fitarr.shape()[3]; // the new element position 1173 } else { 1174 if (fitarr.nelements() > 0) 1175 newarr = fitarr; 1176 } 1177 newarr(IPosition(4, beamSel_, IFSel_, polSel_, pos)) = Int(nrow); 1178 fitCol_.put(whichRow, newarr); 1179 1180 } 1181 1182 SDFrequencyTable SDMemTable::getSDFreqTable() const 1183 { 1184 const Table& t = table_.keywordSet().asTable("FREQUENCIES"); 1185 SDFrequencyTable sdft; 1186 1187 // Add refpix/refval/incr. What are the units ? Hz I suppose 1188 // but it's nowhere described... 1189 Vector<Double> refPix, refVal, incr; 1190 ScalarColumn<Double> refPixCol(t, "REFPIX"); 1191 ScalarColumn<Double> refValCol(t, "REFVAL"); 1192 ScalarColumn<Double> incrCol(t, "INCREMENT"); 1193 refPix = refPixCol.getColumn(); 1194 refVal = refValCol.getColumn(); 1195 incr = incrCol.getColumn(); 1196 1197 uInt n = refPix.nelements(); 1198 for (uInt i=0; i<n; i++) { 1199 sdft.addFrequency(refPix[i], refVal[i], incr[i]); 1200 } 1201 1202 // Frequency reference frame. I don't know if this 1203 // is the correct frame. It might be 'REFFRAME' 1204 // rather than 'BASEREFFRAME' ? 1205 String frame; 1206 t.keywordSet().get("REFFRAME",frame); 1207 sdft.setRefFrame(frame); 1208 t.keywordSet().get("BASEREFFRAME",frame); 1209 sdft.setBaseRefFrame(frame); 1210 1211 // Equinox 1212 Float equinox; 1213 t.keywordSet().get("EQUINOX", equinox); 1214 sdft.setEquinox(equinox); 1215 1216 String unit; 1217 t.keywordSet().get("UNIT", unit); 1218 sdft.setUnit(unit); 1219 1220 // Rest Frequency 1221 Vector<Double> restFreqs; 1222 t.keywordSet().get("RESTFREQS", restFreqs); 1223 for (uInt i=0; i<restFreqs.nelements(); i++) { 1224 sdft.addRestFrequency(restFreqs[i]); 1225 } 1226 sdft.setRestFrequencyUnit(String("Hz")); 1227 1228 return sdft; 1229 } 1230 1231 bool SDMemTable::putSDContainer(const SDContainer& sdc) 1232 { 1233 uInt rno = table_.nrow(); 1234 table_.addRow(); 1235 1236 timeCol_.put(rno, sdc.timestamp); 1237 srcnCol_.put(rno, sdc.sourcename); 1238 fldnCol_.put(rno, sdc.fieldname); 1239 specCol_.put(rno, sdc.getSpectrum()); 1240 flagsCol_.put(rno, sdc.getFlags()); 1241 tsCol_.put(rno, sdc.getTsys()); 1242 scanCol_.put(rno, sdc.scanid); 1243 integrCol_.put(rno, sdc.interval); 1244 freqidCol_.put(rno, sdc.getFreqMap()); 1245 restfreqidCol_.put(rno, sdc.getRestFreqMap()); 1246 dirCol_.put(rno, sdc.getDirection()); 1247 rbeamCol_.put(rno, sdc.refbeam); 1248 tcalCol_.put(rno, sdc.tcal); 1249 tcaltCol_.put(rno, sdc.tcaltime); 1250 azCol_.put(rno, sdc.azimuth); 1251 elCol_.put(rno, sdc.elevation); 1252 paraCol_.put(rno, sdc.parangle); 1253 fitCol_.put(rno, sdc.getFitMap()); 1254 return true; 1255 } 1256 1257 SDContainer SDMemTable::getSDContainer(uInt whichRow) const 1258 { 1259 SDContainer sdc(nBeam(),nIF(),nPol(),nChan()); 1260 timeCol_.get(whichRow, sdc.timestamp); 1261 srcnCol_.get(whichRow, sdc.sourcename); 1262 integrCol_.get(whichRow, sdc.interval); 1263 scanCol_.get(whichRow, sdc.scanid); 1264 fldnCol_.get(whichRow, sdc.fieldname); 1265 rbeamCol_.get(whichRow, sdc.refbeam); 1266 azCol_.get(whichRow, sdc.azimuth); 1267 elCol_.get(whichRow, sdc.elevation); 1268 paraCol_.get(whichRow, sdc.parangle); 1269 Vector<Float> tc; 1270 tcalCol_.get(whichRow, tc); 1271 sdc.tcal[0] = tc[0];sdc.tcal[1] = tc[1]; 1272 tcaltCol_.get(whichRow, sdc.tcaltime); 1273 1274 Array<Float> spectrum; 1275 Array<Float> tsys; 1276 Array<uChar> flagtrum; 1277 Vector<uInt> fmap; 1278 Array<Double> direction; 1279 Array<Int> fits; 1280 1281 specCol_.get(whichRow, spectrum); 1282 sdc.putSpectrum(spectrum); 1283 flagsCol_.get(whichRow, flagtrum); 1284 sdc.putFlags(flagtrum); 1285 tsCol_.get(whichRow, tsys); 1286 sdc.putTsys(tsys); 1287 freqidCol_.get(whichRow, fmap); 1288 sdc.putFreqMap(fmap); 1289 restfreqidCol_.get(whichRow, fmap); 1290 sdc.putRestFreqMap(fmap); 1291 dirCol_.get(whichRow, direction); 1292 sdc.putDirection(direction); 1293 fitCol_.get(whichRow, fits); 1294 sdc.putFitMap(fits); 1295 return sdc; 1296 } 1297 1298 bool SDMemTable::putSDHeader(const SDHeader& sdh) 211 } 212 213 void Scantable::attach() 214 { 215 historyTable_ = table_.keywordSet().asTable("HISTORY"); 216 fitTable_ = table_.keywordSet().asTable("FITS"); 217 218 timeCol_.attach(table_, "TIME"); 219 srcnCol_.attach(table_, "SRCNAME"); 220 specCol_.attach(table_, "SPECTRA"); 221 flagsCol_.attach(table_, "FLAGTRA"); 222 tsCol_.attach(table_, "TSYS"); 223 cycleCol_.attach(table_,"CYCLENO"); 224 scanCol_.attach(table_, "SCANNO"); 225 beamCol_.attach(table_, "BEAMNO"); 226 integrCol_.attach(table_, "INTERVAL"); 227 azCol_.attach(table_, "AZIMUTH"); 228 elCol_.attach(table_, "ELEVATION"); 229 dirCol_.attach(table_, "DIRECTION"); 230 paraCol_.attach(table_, "PARANGLE"); 231 fldnCol_.attach(table_, "FIELDNAME"); 232 rbeamCol_.attach(table_, "REFBEAMNO"); 233 234 mfitidCol_.attach(table_,"FIT_ID"); 235 fitidCol_.attach(fitTable_,"FIT_ID"); 236 237 mfreqidCol_.attach(table_, "FREQ_ID"); 238 239 mtcalidCol_.attach(table_, "TCAL_ID"); 240 241 mfocusidCol_.attach(table_, "FOCUS_ID"); 242 243 mmolidCol_.attach(table_, "MOLECULE_ID"); 244 } 245 246 void Scantable::putSDHeader(const SDHeader& sdh) 1299 247 { 1300 248 table_.rwKeywordSet().define("nIF", sdh.nif); … … 1314 262 table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit); 1315 263 table_.rwKeywordSet().define("Epoch", sdh.epoch); 1316 return true; 1317 } 1318 1319 SDHeader SDMemTable::getSDHeader() const 264 } 265 266 SDHeader Scantable::getSDHeader() const 1320 267 { 1321 268 SDHeader sdh; … … 1338 285 return sdh; 1339 286 } 1340 void SDMemTable::makePersistent(const std::string& filename) 1341 { 1342 table_.deepCopy(filename,Table::New); 1343 1344 } 1345 1346 Int SDMemTable::nScan() const { 1347 Int n = 0; 1348 Int previous = -1;Int current=0; 287 288 int Scantable::rowToScanIndex( int therow ) 289 { 290 int therealrow = -1; 291 292 return therealrow; 293 } 294 295 int Scantable::nScan() const { 296 int n = 0; 297 int previous = -1; int current = 0; 1349 298 for (uInt i=0; i< scanCol_.nrow();i++) { 1350 299 scanCol_.getScalar(i,current); … … 1357 306 } 1358 307 1359 String SDMemTable::formatSec(Double x) const308 std::string Scantable::formatSec(Double x) const 1360 309 { 1361 310 Double xcop = x; … … 1374 323 }; 1375 324 1376 String SDMemTable::formatDirection(const MDirection& md) const325 std::string Scantable::formatDirection(const MDirection& md) const 1377 326 { 1378 327 Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue(); … … 1387 336 1388 337 1389 std::string S DMemTable::getFluxUnit() const338 std::string Scantable::getFluxUnit() const 1390 339 { 1391 340 String tmp; … … 1394 343 } 1395 344 1396 void S DMemTable::setFluxUnit(const std::string& unit)345 void Scantable::setFluxUnit(const std::string& unit) 1397 346 { 1398 347 String tmp(unit); … … 1405 354 } 1406 355 1407 1408 void SDMemTable::setInstrument(const std::string& name) 1409 { 1410 Bool throwIt = True; 356 void Scantable::setInstrument(const std::string& name) 357 { 358 bool throwIt = true; 1411 359 Instrument ins = SDAttr::convertInstrument(name, throwIt); 1412 360 String nameU(name); … … 1415 363 } 1416 364 1417 std::string SDMemTable::summary(bool verbose) const { 1418 1419 // Format header info 1420 ostringstream oss; 1421 oss << endl; 1422 oss << "--------------------------------------------------------------------------------" << endl; 1423 oss << " Scan Table Summary" << endl; 1424 oss << "--------------------------------------------------------------------------------" << endl; 1425 oss.flags(std::ios_base::left); 1426 oss << setw(15) << "Beams:" << setw(4) << nBeam() << endl 1427 << setw(15) << "IFs:" << setw(4) << nIF() << endl 1428 << setw(15) << "Polarisations:" << setw(4) << nPol() << endl 1429 << setw(15) << "Channels:" << setw(4) << nChan() << endl; 1430 oss << endl; 1431 String tmp; 1432 table_.keywordSet().get("Observer", tmp); 1433 oss << setw(15) << "Observer:" << tmp << endl; 1434 oss << setw(15) << "Obs Date:" << getTime(-1,True) << endl; 1435 table_.keywordSet().get("Project", tmp); 1436 oss << setw(15) << "Project:" << tmp << endl; 1437 table_.keywordSet().get("Obstype", tmp); 1438 oss << setw(15) << "Obs. Type:" << tmp << endl; 1439 table_.keywordSet().get("AntennaName", tmp); 1440 oss << setw(15) << "Antenna Name:" << tmp << endl; 1441 table_.keywordSet().get("FluxUnit", tmp); 1442 oss << setw(15) << "Flux Unit:" << tmp << endl; 1443 Table t = table_.keywordSet().asTable("FREQUENCIES"); 1444 Vector<Double> vec; 1445 t.keywordSet().get("RESTFREQS",vec); 1446 oss << setw(15) << "Rest Freqs:"; 1447 if (vec.nelements() > 0) { 1448 oss << setprecision(10) << vec << " [Hz]" << endl; 1449 } else { 1450 oss << "None set" << endl; 1451 } 1452 oss << setw(15) << "Abcissa:" << getAbcissaString() << endl; 1453 oss << setw(15) << "Cursor:" << "Beam[" << getBeam() << "] " 1454 << "IF[" << getIF() << "] " << "Pol[" << getPol() << "]" << endl; 1455 oss << endl; 1456 1457 String dirtype ="Position ("+ MDirection::showType(getDirectionReference()) + ")"; 1458 oss << setw(5) << "Scan" 1459 << setw(15) << "Source" 1460 << setw(24) << dirtype 1461 << setw(10) << "Time" 1462 << setw(18) << "Integration" 1463 << setw(7) << "FreqIDs" << endl; 1464 oss << "--------------------------------------------------------------------------------" << endl; 1465 1466 // Generate list of scan start and end integrations 1467 Vector<Int> scanIDs = scanCol_.getColumn(); 1468 Vector<uInt> startInt, endInt; 1469 mathutil::scanBoundaries(startInt, endInt, scanIDs); 1470 1471 const uInt nScans = startInt.nelements(); 1472 String name; 1473 Vector<uInt> freqIDs, listFQ; 1474 Vector<uInt> restFreqIDs, listRestFQ; 1475 uInt nInt; 1476 1477 for (uInt i=0; i<nScans; i++) { 1478 // Get things from first integration of scan 1479 String time = getTime(startInt(i),False); 1480 String tInt = formatSec(Double(getInterval(startInt(i)))); 1481 String posit = formatDirection(getDirection(startInt(i),True)); 1482 srcnCol_.getScalar(startInt(i),name); 1483 1484 // Find all the FreqIDs in this scan 1485 listFQ.resize(0); 1486 listRestFQ.resize(0); 1487 for (uInt j=startInt(i); j<endInt(i)+1; j++) { 1488 freqidCol_.get(j, freqIDs); 1489 for (uInt k=0; k<freqIDs.nelements(); k++) { 1490 mathutil::addEntry(listFQ, freqIDs(k)); 1491 } 1492 // 1493 restfreqidCol_.get(j, restFreqIDs); 1494 for (uInt k=0; k<restFreqIDs.nelements(); k++) { 1495 mathutil::addEntry(listRestFQ, restFreqIDs(k)); 1496 } 1497 } 1498 1499 nInt = endInt(i) - startInt(i) + 1; 1500 oss << setw(3) << std::right << i << std::left << setw(2) << " " 1501 << setw(15) << name 1502 << setw(24) << posit 1503 << setw(10) << time 1504 << setw(3) << std::right << nInt << setw(3) << " x " << std::left 1505 << setw(6) << tInt 1506 << " " << listFQ << " " << listRestFQ << endl; 1507 } 1508 oss << endl; 1509 oss << "Table contains " << table_.nrow() << " integration(s) in " 1510 << nScans << " scan(s)." << endl; 1511 1512 // Frequency Table 1513 if (verbose) { 1514 std::vector<string> info = getCoordInfo(); 1515 SDFrequencyTable sdft = getSDFreqTable(); 1516 oss << endl << endl; 1517 oss << "FreqID Frame RefFreq(Hz) RefPix Increment(Hz)" << endl; 1518 oss << "--------------------------------------------------------------------------------" << endl; 1519 for (uInt i=0; i<sdft.length(); i++) { 1520 oss << setw(8) << i << setw(8) 1521 << info[3] << setw(16) << setprecision(8) 1522 << sdft.referenceValue(i) << setw(10) 1523 << sdft.referencePixel(i) << setw(12) 1524 << sdft.increment(i) << endl; 1525 } 1526 oss << "--------------------------------------------------------------------------------" << endl; 1527 } 1528 return String(oss); 1529 } 1530 /* 1531 std::string SDMemTable::scanSummary(const std::vector<int>& whichScans) { 1532 ostringstream oss; 1533 Vector<Int> scanIDs = scanCol_.getColumn(); 1534 Vector<uInt> startInt, endInt; 1535 mathutil::scanBoundaries(startInt, endInt, scanIDs); 1536 const uInt nScans = startInt.nelements(); 1537 std::vector<int>::const_iterator it(whichScans); 1538 return String(oss); 1539 } 1540 */ 1541 Int SDMemTable::nBeam() const 1542 { 1543 Int n; 1544 table_.keywordSet().get("nBeam",n); 1545 return n; 1546 } 1547 1548 Int SDMemTable::nIF() const { 1549 Int n; 1550 table_.keywordSet().get("nIF",n); 1551 return n; 1552 } 1553 1554 Int SDMemTable::nPol() const { 1555 Int n; 1556 table_.keywordSet().get("nPol",n); 1557 return n; 1558 } 1559 1560 Int SDMemTable::nChan() const { 1561 Int n; 1562 table_.keywordSet().get("nChan",n); 1563 return n; 1564 } 1565 1566 Table SDMemTable::getHistoryTable() const 365 MPosition Scantable::getAntennaPosition () const 366 { 367 Vector<Double> antpos; 368 table_.keywordSet().get("AntennaPosition", antpos); 369 MVPosition mvpos(antpos(0),antpos(1),antpos(2)); 370 return MPosition(mvpos); 371 } 372 373 void Scantable::makePersistent(const std::string& filename) 374 { 375 String inname(filename); 376 Path path(inname); 377 inname = path.expandedName(); 378 table_.deepCopy(inname, Table::New); 379 } 380 381 int asap::Scantable::nbeam( int scanno ) const 382 { 383 if ( scanno < 0 ) { 384 Int n; 385 table_.keywordSet().get("nBeam",n); 386 return int(n); 387 } else { 388 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these 389 Table tab = table_(table_.col("SCANNO") == scanno 390 && table_.col("POLNO") == 0 391 && table_.col("IFNO") == 0 392 && table_.col("CYCLENO") == 0 ); 393 ROTableVector<uInt> v(tab, "BEAMNO"); 394 return int(v.nelements()); 395 } 396 return 0; 397 } 398 399 int asap::Scantable::nif( int scanno ) const 400 { 401 if ( scanno < 0 ) { 402 Int n; 403 table_.keywordSet().get("nIF",n); 404 return int(n); 405 } else { 406 // take the first POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these 407 Table tab = table_(table_.col("SCANNO") == scanno 408 && table_.col("POLNO") == 0 409 && table_.col("BEAMNO") == 0 410 && table_.col("CYCLENO") == 0 ); 411 ROTableVector<uInt> v(tab, "IFNO"); 412 return int(v.nelements()); 413 } 414 return 0; 415 } 416 417 int asap::Scantable::npol( int scanno ) const 418 { 419 if ( scanno < 0 ) { 420 Int n; 421 table_.keywordSet().get("nPol",n); 422 return n; 423 } else { 424 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these 425 Table tab = table_(table_.col("SCANNO") == scanno 426 && table_.col("IFNO") == 0 427 && table_.col("BEAMNO") == 0 428 && table_.col("CYCLENO") == 0 ); 429 ROTableVector<uInt> v(tab, "POLNO"); 430 return int(v.nelements()); 431 } 432 return 0; 433 } 434 435 int asap::Scantable::nrow( int scanno ) const 436 { 437 if ( scanno < 0 ) { 438 return int(table_.nrow()); 439 } else { 440 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these 441 Table tab = table_(table_.col("SCANNO") == scanno 442 && table_.col("BEAMNO") == 0 443 && table_.col("IFNO") == 0 444 && table_.col("POLNO") == 0 ); 445 return int(tab.nrow()); 446 } 447 return 0; 448 } 449 450 451 int asap::Scantable::nchan( int scanno, int ifno ) const 452 { 453 if ( scanno < 0 && ifno < 0 ) { 454 Int n; 455 table_.keywordSet().get("nChan",n); 456 return int(n); 457 } else { 458 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these 459 Table tab = table_(table_.col("SCANNO") == scanno 460 && table_.col("IFNO") == ifno 461 && table_.col("BEAMNO") == 0 462 && table_.col("POLNO") == 0 463 && table_.col("CYCLENO") == 0 ); 464 ROArrayColumn<Float> v(tab, "SPECTRA"); 465 cout << v.shape(0) << endl; 466 return 0; 467 } 468 return 0; 469 } 470 471 Table Scantable::getHistoryTable() const 1567 472 { 1568 473 return table_.keywordSet().asTable("HISTORY"); 1569 474 } 1570 475 1571 void S DMemTable::appendToHistoryTable(const Table& otherHist)476 void Scantable::appendToHistoryTable(const Table& otherHist) 1572 477 { 1573 478 Table t = table_.rwKeywordSet().asTable("HISTORY"); 1574 const String sep = "--------------------------------------------------------------------------------"; 1575 addHistory( sep);479 480 addHistory(asap::SEPERATOR); 1576 481 TableCopy::copyRows(t, otherHist, t.nrow(), 0, otherHist.nrow()); 1577 addHistory( sep);1578 } 1579 1580 void S DMemTable::addHistory(const std::string& hist)482 addHistory(asap::SEPERATOR); 483 } 484 485 void Scantable::addHistory(const std::string& hist) 1581 486 { 1582 487 Table t = table_.rwKeywordSet().asTable("HISTORY"); … … 1587 492 } 1588 493 1589 std::vector<std::string> S DMemTable::getHistory() const494 std::vector<std::string> Scantable::getHistory() const 1590 495 { 1591 496 Vector<String> history; … … 1602 507 } 1603 508 1604 1605 /* 1606 void SDMemTable::maskChannels(const std::vector<Int>& whichChans ) { 1607 1608 std::vector<int>::iterator it; 1609 ArrayAccessor<uChar, Axis<asap::PolAxis> > j(flags_); 1610 for (it = whichChans.begin(); it != whichChans.end(); it++) { 1611 j.reset(j.begin(uInt(*it))); 1612 for (ArrayAccessor<uChar, Axis<asap::BeamAxis> > i(j); i != i.end(); ++i) { 1613 for (ArrayAccessor<uChar, Axis<asap::IFAxis> > ii(i); ii != ii.end(); ++ii) { 1614 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > iii(ii); 1615 iii != iii.end(); ++iii) { 1616 (*iii) = 1617 } 1618 } 1619 } 1620 } 1621 1622 } 1623 */ 1624 void SDMemTable::flag(int whichRow) 1625 { 1626 Array<uChar> arr; 1627 flagsCol_.get(whichRow, arr); 1628 1629 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr); 1630 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam 1631 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0); 1632 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF 1633 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1); 1634 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol 1635 1636 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) { 1637 (*i) = uChar(True); 1638 } 1639 1640 flagsCol_.put(whichRow, arr); 1641 } 1642 1643 MDirection::Types SDMemTable::getDirectionReference() const 1644 { 1645 Float eq; 1646 table_.keywordSet().get("Equinox",eq); 1647 std::map<float,string> mp; 1648 mp[2000.0] = "J2000"; 1649 mp[1950.0] = "B1950"; 1650 MDirection::Types mdr; 1651 if (!MDirection::getType(mdr, mp[eq])) { 1652 mdr = MDirection::J2000; 1653 pushLog("WARNING: Unknown equinox using J2000"); 1654 } 1655 1656 return mdr; 1657 } 1658 1659 MEpoch::Types SDMemTable::getTimeReference() const 1660 { 1661 MEpoch::Types met; 1662 String ep; 1663 table_.keywordSet().get("Epoch",ep); 1664 if (!MEpoch::getType(met, ep)) { 1665 pushLog("WARNING: Epoch type unknown - using UTC"); 1666 met = MEpoch::UTC; 1667 } 1668 1669 return met; 1670 } 1671 1672 1673 Bool SDMemTable::setRestFreqs(const Vector<Double>& restFreqsIn, 1674 const String& sUnit, 1675 const vector<string>& lines, 1676 const String& source, 1677 Int whichIF) 1678 { 1679 const Int nIFs = nIF(); 1680 if (whichIF>=nIFs) { 1681 throw(AipsError("Illegal IF index")); 1682 } 1683 1684 // Find vector of restfrequencies 1685 // Double takes precedence over String 1686 Unit unit; 1687 Vector<Double> restFreqs; 1688 if (restFreqsIn.nelements()>0) { 1689 restFreqs.resize(restFreqsIn.nelements()); 1690 restFreqs = restFreqsIn; 1691 unit = Unit(sUnit); 1692 } else if (lines.size()>0) { 1693 const uInt nLines = lines.size(); 1694 unit = Unit("Hz"); 1695 restFreqs.resize(nLines); 1696 MFrequency lineFreq; 1697 for (uInt i=0; i<nLines; i++) { 1698 String tS(lines[i]); 1699 tS.upcase(); 1700 if (MeasTable::Line(lineFreq, tS)) { 1701 restFreqs[i] = lineFreq.getValue().getValue(); // Hz 1702 } else { 1703 String s = String(lines[i]) + 1704 String(" is an unrecognized spectral line"); 1705 throw(AipsError(s)); 1706 } 1707 } 1708 } else { 1709 throw(AipsError("You have not specified any rest frequencies or lines")); 1710 } 1711 1712 // If multiple restfreqs, must be length nIF. In this 1713 // case we will just replace the rest frequencies 1714 // We can't disinguish scalar and vector of length 1 1715 const uInt nRestFreqs = restFreqs.nelements(); 1716 Int idx = -1; 1717 SDFrequencyTable sdft = getSDFreqTable(); 1718 1719 ostringstream oss; 1720 if (nRestFreqs>1) { 1721 // Replace restFreqs, one per IF 1722 if (nRestFreqs != nIFs) { 1723 throw (AipsError("Number of rest frequencies must be equal to the number of IFs")); 1724 } 1725 ostringstream oss; 1726 oss << "Replaced rest frequencies, one per IF, with given list : " << restFreqs; 1727 sdft.deleteRestFrequencies(); 1728 for (uInt i=0; i<nRestFreqs; i++) { 1729 Quantum<Double> rf(restFreqs[i], unit); 1730 sdft.addRestFrequency(rf.getValue("Hz")); 1731 } 1732 } else { 1733 1734 // Add new rest freq 1735 Quantum<Double> rf(restFreqs[0], unit); 1736 idx = sdft.addRestFrequency(rf.getValue("Hz")); 1737 if (whichIF>=0) { 1738 oss << "Selected given rest frequency (" << restFreqs[0] << ") for IF " << whichIF << endl; 1739 } else { 1740 oss << "Selected given rest frequency (" << restFreqs[0] << ") for all IFs" << endl; 1741 } 1742 } 1743 pushLog(String(oss)); 1744 // Replace 1745 Bool empty = source.empty(); 1746 Bool ok = False; 1747 if (putSDFreqTable(sdft)) { 1748 const uInt nRow = table_.nrow(); 1749 String srcName; 1750 Vector<uInt> restFreqIDs; 1751 for (uInt i=0; i<nRow; i++) { 1752 srcnCol_.get(i, srcName); 1753 restfreqidCol_.get(i,restFreqIDs); 1754 if (idx==-1) { 1755 // Replace vector of restFreqs; one per IF. 1756 // No selection possible 1757 for (uInt i=0; i<nIFs; i++) restFreqIDs[i] = i; 1758 } else { 1759 // Set RestFreqID for selected data 1760 if (empty || source==srcName) { 1761 if (whichIF<0) { 1762 restFreqIDs = idx; 1763 } else { 1764 restFreqIDs[whichIF] = idx; 1765 } 1766 } 1767 } 1768 restfreqidCol_.put(i,restFreqIDs); 1769 } 1770 ok = True; 1771 } else { 1772 ok = False; 1773 } 1774 1775 return ok; 1776 } 1777 1778 std::string SDMemTable::spectralLines() const 1779 { 1780 Vector<String> lines = MeasTable::Lines(); 1781 MFrequency lineFreq; 1782 Double freq; 1783 ostringstream oss; 1784 1785 oss.flags(std::ios_base::left); 1786 oss << "Line Frequency (Hz)" << endl; 1787 oss << "-----------------------" << endl; 1788 for (uInt i=0; i<lines.nelements(); i++) { 1789 MeasTable::Line(lineFreq, lines[i]); 1790 freq = lineFreq.getValue().getValue(); // Hz 1791 oss << setw(11) << lines[i] << setprecision(10) << freq << endl; 1792 } 1793 return String(oss); 1794 } 1795 1796 void SDMemTable::renumber() 1797 { 1798 uInt nRow = scanCol_.nrow(); 1799 Int newscanid = 0; 1800 Int cIdx;// the current scanid 1801 // get the first scanid 1802 scanCol_.getScalar(0,cIdx); 1803 Int pIdx = cIdx;// the scanid of the previous row 1804 for (uInt i=0; i<nRow;++i) { 1805 scanCol_.getScalar(i,cIdx); 1806 if (pIdx == cIdx) { 1807 // renumber 1808 scanCol_.put(i,newscanid); 1809 } else { 1810 ++newscanid; 1811 pIdx = cIdx; // store scanid 1812 --i; // don't increment next loop 1813 } 1814 } 1815 } 1816 1817 1818 void SDMemTable::getCursorSlice(IPosition& start, IPosition& end, 1819 const IPosition& shape) const 1820 { 1821 const uInt nDim = shape.nelements(); 1822 start.resize(nDim); 1823 end.resize(nDim); 1824 1825 start(asap::BeamAxis) = beamSel_; 1826 end(asap::BeamAxis) = beamSel_; 1827 start(asap::IFAxis) = IFSel_; 1828 end(asap::IFAxis) = IFSel_; 1829 1830 start(asap::PolAxis) = polSel_; 1831 end(asap::PolAxis) = polSel_; 1832 1833 start(asap::ChanAxis) = 0; 1834 end(asap::ChanAxis) = shape(asap::ChanAxis) - 1; 1835 } 1836 1837 1838 std::vector<float> SDMemTable::getFloatSpectrum(const Array<Float>& arr) const 1839 // Get spectrum at cursor location 1840 { 1841 1842 // Setup accessors 1843 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr); 1844 aa0.reset(aa0.begin(uInt(beamSel_))); // Beam selection 1845 1846 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0); 1847 aa1.reset(aa1.begin(uInt(IFSel_))); // IF selection 1848 1849 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1); 1850 aa2.reset(aa2.begin(uInt(polSel_))); // Pol selection 1851 1852 std::vector<float> spectrum; 1853 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) { 1854 spectrum.push_back(*i); 1855 } 1856 return spectrum; 1857 } 1858 1859 void SDMemTable::calculateAZEL() 509 std::string Scantable::formatTime(const MEpoch& me, bool showdate) const 510 { 511 MVTime mvt(me.getValue()); 512 if (showdate) 513 mvt.setFormat(MVTime::YMD); 514 else 515 mvt.setFormat(MVTime::TIME); 516 ostringstream oss; 517 oss << mvt; 518 return String(oss); 519 } 520 521 void Scantable::calculateAZEL() 1860 522 { 1861 523 MPosition mp = getAntennaPosition(); 524 MEpoch::ROScalarColumn timeCol(table_, "TIME"); 1862 525 ostringstream oss; 1863 526 oss << "Computed azimuth/elevation using " << endl 1864 527 << mp << endl; 1865 for (uInt i=0; i<nRow();++i) { 1866 MEpoch me = getEpoch(i); 1867 MDirection md = getDirection(i,False); 1868 oss << " Time: " << getTime(i,False) << " Direction: " << formatDirection(md) 528 for (uInt i=0; i<nrow(); ++i) { 529 MEpoch me = timeCol(i); 530 MDirection md = dirCol_(i); 531 dirCol_.get(i,md); 532 oss << " Time: " << formatTime(me,False) << " Direction: " << formatDirection(md) 1869 533 << endl << " => "; 1870 534 MeasFrame frame(mp, me); … … 1880 544 pushLog(String(oss)); 1881 545 } 546 547 double Scantable::getInterval(int whichrow) const 548 { 549 if (whichrow < 0) return 0.0; 550 Double intval; 551 integrCol_.get(Int(whichrow), intval); 552 return intval; 553 } 554 555 std::vector<bool> Scantable::getMask(int whichrow) const 556 { 557 Vector<uChar> flags; 558 flagsCol_.get(uInt(whichrow), flags); 559 Vector<Bool> bflag(flags.shape()); 560 convertArray(bflag, flags); 561 bflag = !bflag; 562 std::vector<bool> mask; 563 bflag.tovector(mask); 564 return mask; 565 } 566 567 std::vector<float> Scantable::getSpectrum(int whichrow) const 568 { 569 Vector<Float> arr; 570 specCol_.get(whichrow, arr); 571 std::vector<float> out; 572 arr.tovector(out); 573 return out; 574 } 575 576 String Scantable::generateName() 577 { 578 return (File::newUniqueName("./","temp")).baseName(); 579 } 580 581 const casa::Table& Scantable::table( ) const 582 { 583 return table_; 584 } 585 586 casa::Table& Scantable::table( ) 587 { 588 return table_; 589 } 590 591 std::string Scantable::getPolarizationLabel(bool linear, bool stokes, 592 bool linpol, int polidx) const 593 { 594 uInt idx = 0; 595 if (polidx >=0) idx = polidx; 596 return ""; 597 //return SDPolUtil::polarizationLabel(idx, linear, stokes, linpol); 598 } 599 600 void Scantable::unsetSelection() 601 { 602 table_ = originalTable_; 603 selector_.reset(); 604 } 605 606 void Scantable::setSelection( const STSelector& selection ) 607 { 608 Table tab = const_cast<STSelector&>(selection).apply(originalTable_); 609 if ( tab.nrow() == 0 ) { 610 throw(AipsError("Selection contains no data. Not applying it.")); 611 } 612 table_ = tab; 613 selector_ = selection; 614 } 615 616 std::string asap::Scantable::summary( bool verbose ) 617 { 618 // Format header info 619 ostringstream oss; 620 oss << endl; 621 oss << asap::SEPERATOR << endl; 622 oss << " Scan Table Summary" << endl; 623 oss << asap::SEPERATOR << endl; 624 oss.flags(std::ios_base::left); 625 oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl 626 << setw(15) << "IFs:" << setw(4) << nif() << endl 627 << setw(15) << "Polarisations:" << setw(4) << npol() << endl 628 << setw(15) << "Channels:" << setw(4) << nchan() << endl; 629 oss << endl; 630 String tmp; 631 //table_.keywordSet().get("Observer", tmp); 632 oss << setw(15) << "Observer:" << table_.keywordSet().asString("Observer") << endl; 633 oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl; 634 table_.keywordSet().get("Project", tmp); 635 oss << setw(15) << "Project:" << tmp << endl; 636 table_.keywordSet().get("Obstype", tmp); 637 oss << setw(15) << "Obs. Type:" << tmp << endl; 638 table_.keywordSet().get("AntennaName", tmp); 639 oss << setw(15) << "Antenna Name:" << tmp << endl; 640 table_.keywordSet().get("FluxUnit", tmp); 641 oss << setw(15) << "Flux Unit:" << tmp << endl; 642 Vector<Float> vec; 643 oss << setw(15) << "Rest Freqs:"; 644 if (vec.nelements() > 0) { 645 oss << setprecision(10) << vec << " [Hz]" << endl; 646 } else { 647 oss << "none" << endl; 648 } 649 oss << setw(15) << "Abcissa:" << "channel" << endl; 650 oss << selector_.print() << endl; 651 oss << endl; 652 // main table 653 String dirtype = "Position (" 654 + MDirection::showType(dirCol_.getMeasRef().getType()) 655 + ")"; 656 oss << setw(5) << "Scan" 657 << setw(15) << "Source" 658 // << setw(24) << dirtype 659 << setw(10) << "Time" 660 << setw(18) << "Integration" << endl 661 << setw(5) << "" << setw(10) << "Beam" << dirtype << endl 662 << setw(15) << "" << setw(5) << "IF" 663 << setw(8) << "Frame" << setw(16) 664 << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment" <<endl; 665 oss << asap::SEPERATOR << endl; 666 TableIterator iter(table_, "SCANNO"); 667 while (!iter.pastEnd()) { 668 Table subt = iter.table(); 669 ROTableRow row(subt); 670 MEpoch::ROScalarColumn timeCol(subt,"TIME"); 671 const TableRecord& rec = row.get(0); 672 oss << setw(4) << std::right << rec.asuInt("SCANNO") 673 << std::left << setw(1) << "" 674 << setw(15) << rec.asString("SRCNAME") 675 << setw(10) << formatTime(timeCol(0), false); 676 // count the cycles in the scan 677 TableIterator cyciter(subt, "CYCLENO"); 678 int nint = 0; 679 while (!cyciter.pastEnd()) { 680 ++nint; 681 ++cyciter; 682 } 683 oss << setw(3) << std::right << nint << setw(3) << " x " << std::left 684 << setw(6) << formatSec(rec.asFloat("INTERVAL")) << endl; 685 686 TableIterator biter(subt, "BEAMNO"); 687 while (!biter.pastEnd()) { 688 Table bsubt = biter.table(); 689 ROTableRow brow(bsubt); 690 MDirection::ROScalarColumn bdirCol(bsubt,"DIRECTION"); 691 const TableRecord& brec = brow.get(0); 692 oss << setw(6) << "" << setw(10) << brec.asuInt("BEAMNO"); 693 oss << setw(24) << formatDirection(bdirCol(0)) << endl; 694 TableIterator iiter(bsubt, "IFNO"); 695 while (!iiter.pastEnd()) { 696 Table isubt = iiter.table(); 697 ROTableRow irow(isubt); 698 const TableRecord& irec = irow.get(0); 699 oss << std::right <<setw(8) << "" << std::left << irec.asuInt("IFNO"); 700 oss << frequencies().print(irec.asuInt("FREQ_ID")); 701 702 ++iiter; 703 } 704 ++biter; 705 } 706 ++iter; 707 } 708 /// @todo implement verbose mode 709 return String(oss); 710 } 711 712 std::string Scantable::getTime(int whichrow, bool showdate) const 713 { 714 MEpoch::ROScalarColumn timeCol(table_, "TIME"); 715 MEpoch me; 716 if (whichrow > -1) { 717 me = timeCol(uInt(whichrow)); 718 } else { 719 Double tm; 720 table_.keywordSet().get("UTC",tm); 721 me = MEpoch(MVEpoch(tm)); 722 } 723 return formatTime(me, showdate); 724 } 725 726 }//namespace asap
Note:
See TracChangeset
for help on using the changeset viewer.