source: tags/asap-4.1.0/src/STFrequencies.cpp@ 2747

Last change on this file since 2747 was 2662, checked in by Malte Marquarding, 12 years ago

merge from trunk into release candidate

File size: 12.3 KB
RevLine 
[806]1//
2// C++ Implementation: STFrequencies
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#include <casa/iostream.h>
13#include <casa/iomanip.h>
14#include <casa/Exceptions/Error.h>
15#include <casa/Containers/RecordField.h>
16#include <casa/Arrays/IPosition.h>
[2662]17#include <casa/Logging/LogIO.h>
[806]18
19#include <tables/Tables/TableDesc.h>
20#include <tables/Tables/SetupNewTab.h>
21#include <tables/Tables/ScaColDesc.h>
22#include <tables/Tables/TableRecord.h>
23#include <tables/Tables/TableParse.h>
24#include <tables/Tables/TableRow.h>
25
26#include <coordinates/Coordinates/CoordinateSystem.h>
27#include <coordinates/Coordinates/CoordinateUtil.h>
28
29#include "STFrequencies.h"
30
31
32using namespace casa;
33
34namespace asap {
35
[847]36const String STFrequencies::name_ = "FREQUENCIES";
[806]37
[849]38STFrequencies::STFrequencies(const Scantable& parent) :
39 STSubTable(parent, name_)
[806]40{
41 setup();
42}
43
[1360]44STFrequencies::STFrequencies( casa::Table tab ) :
[856]45 STSubTable(tab, name_)
[849]46{
47 refpixCol_.attach(table_,"REFPIX");
48 refvalCol_.attach(table_,"REFVAL");
49 incrCol_.attach(table_,"INCREMENT");
[806]50
[849]51}
52
[806]53STFrequencies::~STFrequencies()
54{
55}
56
[1360]57STFrequencies & STFrequencies::operator=( const STFrequencies & other )
[849]58{
59 if ( this != &other ) {
60 static_cast<STSubTable&>(*this) = other;
61 refpixCol_.attach(table_,"REFPIX");
62 refvalCol_.attach(table_,"REFVAL");
63 incrCol_.attach(table_,"INCREMENT");
64 }
65 return *this;
66}
67
[806]68void STFrequencies::setup( )
69{
70 // add to base class table
71 table_.addColumn(ScalarColumnDesc<Double>("REFPIX"));
72 table_.addColumn(ScalarColumnDesc<Double>("REFVAL"));
73 table_.addColumn(ScalarColumnDesc<Double>("INCREMENT"));
74
[847]75 table_.rwKeywordSet().define("FRAME", String("TOPO"));
76 table_.rwKeywordSet().define("BASEFRAME", String("TOPO"));
[806]77 table_.rwKeywordSet().define("EQUINOX",String( "J2000"));
[941]78 table_.rwKeywordSet().define("UNIT", String(""));
[806]79 table_.rwKeywordSet().define("DOPPLER", String("RADIO"));
80
81 // new cached columns
82 refpixCol_.attach(table_,"REFPIX");
83 refvalCol_.attach(table_,"REFVAL");
84 incrCol_.attach(table_,"INCREMENT");
85}
86
87uInt STFrequencies::addEntry( Double refpix, Double refval, Double inc )
88{
89 // test if this already exists
90 Table result = table_( near(table_.col("REFVAL"), refval)
[2242]91 && near(table_.col("REFPIX"), refpix)
92 && near(table_.col("INCREMENT"), inc), 1 );
[806]93 uInt resultid = 0;
94 if ( result.nrow() > 0) {
95 ROScalarColumn<uInt> c(result, "ID");
96 c.get(0, resultid);
97
98 } else {
99 uInt rno = table_.nrow();
100 table_.addRow();
101 // get last assigned freq_id and increment
102 if ( rno > 0 ) {
103 idCol_.get(rno-1, resultid);
104 resultid++;
105 }
106 refpixCol_.put(rno, refpix);
107 refvalCol_.put(rno, refval);
108 incrCol_.put(rno, inc);
109 idCol_.put(rno, resultid);
110 }
111 return resultid;
112}
113
114
[830]115
116void STFrequencies::getEntry( Double& refpix, Double& refval, Double& inc,
117 uInt id )
118{
[2243]119 Table t = table_(table_.col("ID") == Int(id), 1 );
[830]120 if (t.nrow() == 0 ) {
121 throw(AipsError("STFrequencies::getEntry - freqID out of range"));
122 }
123 ROTableRow row(t);
124 // get first row - there should only be one matching id
125 const TableRecord& rec = row.get(0);
126 refpix = rec.asDouble("REFPIX");
127 refval = rec.asDouble("REFVAL");
128 inc = rec.asDouble("INCREMENT");
129}
130
[1819]131void STFrequencies::setEntry( Double refpix, Double refval, Double inc, uInt id )
132{
[2243]133 Table t = table_(table_.col("ID") == Int(id), 1 );
[1819]134 if (t.nrow() == 0 ) {
135 throw(AipsError("STFrequencies::getEntry - freqID out of range"));
136 }
137 for ( uInt i = 0 ; i < table_.nrow() ; i++ ) {
138 uInt fid ;
139 idCol_.get( i, fid ) ;
140 if ( fid == id ) {
141 refpixCol_.put( i, refpix ) ;
142 refvalCol_.put( i, refval ) ;
143 incrCol_.put( i, inc ) ;
144 }
145 }
146}
147
[847]148SpectralCoordinate STFrequencies::getSpectralCoordinate( uInt id ) const
[806]149{
[2243]150 Table t = table_(table_.col("ID") == Int(id), 1 );
[806]151
152 if (t.nrow() == 0 ) {
[847]153 throw(AipsError("STFrequencies::getSpectralCoordinate - ID out of range"));
[806]154 }
155
156 // get the data
157 ROTableRow row(t);
158 // get first row - there should only be one matching id
159 const TableRecord& rec = row.get(0);
[847]160 return SpectralCoordinate( getFrame(true), rec.asDouble("REFVAL"),
[806]161 rec.asDouble("INCREMENT"),
162 rec.asDouble("REFPIX"));
163}
164
[1819]165/**
[847]166SpectralCoordinate
[1360]167 STFrequencies::getSpectralCoordinate( const MDirection& md,
168 const MPosition& mp,
169 const MEpoch& me,
170 Double restfreq, uInt id ) const
[1819]171**/
172SpectralCoordinate
173 STFrequencies::getSpectralCoordinate( const MDirection& md,
174 const MPosition& mp,
175 const MEpoch& me,
176 Vector<Double> restfreq, uInt id ) const
[806]177{
[847]178 SpectralCoordinate spc = getSpectralCoordinate(id);
[1819]179 //spc.setRestFrequency(restfreq, True);
180 // for now just use the first rest frequency
181 if (restfreq.nelements()==0 ) {
182 restfreq.resize(1);
183 restfreq[0] = 0;
184 }
185 spc.setRestFrequency(restfreq[0], True);
[847]186 if ( !spc.setReferenceConversion(getFrame(), me, mp, md) ) {
187 throw(AipsError("Couldn't convert frequency frame."));
188 }
189 String unitstr = getUnitString();
190 if ( !unitstr.empty() ) {
191 Unit unitu(unitstr);
192 if ( unitu == Unit("Hz") ) {
[866]193 Vector<String> wau(1); wau = unitu.getName();
194 spc.setWorldAxisUnits(wau);
[847]195 } else {
196 spc.setVelocity(unitstr, getDoppler());
197 }
198 }
199 return spc;
200}
201
202
203void STFrequencies::rescale( Float factor, const std::string& mode )
204{
[806]205 TableRow row(table_);
206 TableRecord& outrec = row.record();
207 RecordFieldPtr<Double> rv(outrec, "REFVAL");
208 RecordFieldPtr<Double> rp(outrec, "REFPIX");
209 RecordFieldPtr<Double> inc(outrec, "INCREMENT");
210 for (uInt i=0; i<table_.nrow(); ++i) {
211
212 const TableRecord& rec = row.get(i);
213
[847]214 SpectralCoordinate sc ( getFrame(true), rec.asDouble("REFVAL"),
[806]215 rec.asDouble("INCREMENT"), rec.asDouble("REFPIX") );
216
217 SpectralCoordinate scout;
218 if (mode == "BIN") {
219 scout = binCsys(sc, Int(factor));
220 } else if (mode == "RESAMPLE") {
221 scout = resampleCsys(sc, factor);
222 }
223 *rv = scout.referenceValue()[0];
224 *rp = scout.referencePixel()[0];
225 *inc = scout.increment()[0];
226 row.put(i);
227 }
228}
229
230SpectralCoordinate STFrequencies::binCsys(const SpectralCoordinate& sc,
231 Int factor)
232{
233 CoordinateSystem csys;
234 csys.addCoordinate(sc);
235 IPosition factors(1, factor);
236 CoordinateSystem binnedcs =
237 CoordinateUtil::makeBinnedCoordinateSystem(factors, csys, False);
238 return binnedcs.spectralCoordinate(0);
239}
240
241SpectralCoordinate STFrequencies::resampleCsys(const SpectralCoordinate& sc,
242 Float width)
243{
244 Vector<Float> offset(1,0.0);
[1694]245 Vector<Float> factors(1,width);
[806]246 Vector<Int> newshape;
247 CoordinateSystem csys;
248 csys.addCoordinate(sc);
249 CoordinateSystem csys2 = csys.subImage(offset, factors, newshape);
250 return csys2.spectralCoordinate(0);
251}
252
253
[847]254MFrequency::Types STFrequencies::getFrame(bool base) const
[806]255{
256 // get the ref frame
[921]257 String rf;
258 if ( base )
259 rf = table_.keywordSet().asString("BASEFRAME");
260 else
261 rf = table_.keywordSet().asString("FRAME");
[806]262
263 // Create SpectralCoordinate (units Hz)
264 MFrequency::Types mft;
265 if (!MFrequency::getType(mft, rf)) {
266 ostringstream oss;
[2662]267 LogIO os( casa::LogOrigin( "STFrequencies", "getFrame") );
268 os << LogIO::WARN << "WARNING: Frequency type unknown assuming TOPO"
269 << LogIO::POST;
[806]270 mft = MFrequency::TOPO;
271 }
272
273 return mft;
274}
275
[1360]276std::string STFrequencies::getFrameString( bool base ) const
[847]277{
278 if ( base ) return table_.keywordSet().asString("BASEFRAME");
279 else return table_.keywordSet().asString("FRAME");
280}
281
[1360]282std::string STFrequencies::getUnitString( ) const
[847]283{
284 return table_.keywordSet().asString("UNIT");
285}
286
[1360]287Unit STFrequencies::getUnit( ) const
[847]288{
289 return Unit(table_.keywordSet().asString("UNIT"));
290}
291
[1360]292std::string STFrequencies::getDopplerString( ) const
[847]293{
294 return table_.keywordSet().asString("DOPPLER");
295}
296
[1360]297MDoppler::Types STFrequencies::getDoppler( ) const
[847]298{
299 String dpl = table_.keywordSet().asString("DOPPLER");
300
301 // Create SpectralCoordinate (units Hz)
302 MDoppler::Types mdt;
303 if (!MDoppler::getType(mdt, dpl)) {
304 throw(AipsError("Doppler type unknown"));
305 }
306 return mdt;
307}
308
[1375]309std::string STFrequencies::print( int id, Bool strip ) const
[806]310{
311 Table t;
312 ostringstream oss;
313 if ( id < 0 ) t = table_;
[2243]314 else t = table_(table_.col("ID") == Int(id), 1 );
[806]315 ROTableRow row(t);
316 for (uInt i=0; i<t.nrow(); ++i) {
317 const TableRecord& rec = row.get(i);
318 oss << setw(8)
[1375]319 << t.keywordSet().asString("FRAME") << setw(16) << setprecision(8)
320 << rec.asDouble("REFVAL") << setw(7)
321 << rec.asDouble("REFPIX")
[1694]322 << setw(15)
[1375]323 << rec.asDouble("INCREMENT");
[806]324 }
[1375]325 String outstr(oss);
326 if (strip) {
327 int f = outstr.find_first_not_of(' ');
328 int l = outstr.find_last_not_of(' ', outstr.size());
[1694]329 if (f < 0) {
330 f = 0;
[1375]331 }
[1694]332 if ( l < f || l < f ) {
[1375]333 l = outstr.size();
334 }
335 return outstr.substr(f,l);
336 }
337 return outstr;
[806]338}
339
340float STFrequencies::getRefFreq( uInt id, uInt channel )
341{
[2243]342 Table t = table_(table_.col("ID") == Int(id), 1 );
[806]343 if ( t.nrow() == 0 ) throw(AipsError("Selected Illegal frequency id"));
344 ROTableRow row(t);
345 const TableRecord& rec = row.get(0);
346 return (Double(channel/2) - rec.asDouble("REFPIX"))
347 * rec.asDouble("INCREMENT") + rec.asDouble("REFVAL");
348}
349
[1360]350bool STFrequencies::conformant( const STFrequencies& other ) const
[840]351{
352 const Record& r = table_.keywordSet();
353 const Record& ro = other.table_.keywordSet();
[847]354 return ( r.asString("FRAME") == ro.asString("FRAME") &&
[840]355 r.asString("EQUINOX") == ro.asString("EQUINOX") &&
356 r.asString("UNIT") == ro.asString("UNIT") &&
357 r.asString("DOPPLER") == ro.asString("DOPPLER")
358 );
359}
360
[1360]361std::vector< std::string > STFrequencies::getInfo( ) const
[847]362{
363 const Record& r = table_.keywordSet();
364 std::vector<std::string> out;
365 out.push_back(r.asString("UNIT"));
366 out.push_back(r.asString("FRAME"));
367 out.push_back(r.asString("DOPPLER"));
[866]368 return out;
[847]369}
370
[1360]371void STFrequencies::setInfo( const std::vector< std::string >& theinfo )
[847]372{
373 if ( theinfo.size() != 3 ) throw(AipsError("setInfo needs three parameters"));
[866]374 try {
375 setUnit(theinfo[0]);
376 setFrame(theinfo[1]);
377 setDoppler(theinfo[2]);
378 } catch (AipsError& e) {
379 throw(e);
380 }
381}
382
[1360]383void STFrequencies::setUnit( const std::string & unit )
[866]384{
385 if (unit == "" || unit == "pixel" || unit == "channel" ) {
386 table_.rwKeywordSet().define("UNIT", "");
[847]387 } else {
[866]388 Unit u(unit);
389 if ( u == Unit("km/s") || u == Unit("Hz") )
390 table_.rwKeywordSet().define("UNIT", unit);
391 else {
392 throw(AipsError("Illegal spectral unit."));
393 }
[847]394 }
395}
[866]396
[1360]397void STFrequencies::setFrame(MFrequency::Types frame, bool base )
[847]398{
[921]399 String f = MFrequency::showType(frame);
400 if (base)
401 table_.rwKeywordSet().define("BASEFRAME", f);
402 else
403 table_.rwKeywordSet().define("FRAME", f);
404
405}
406
[1360]407void STFrequencies::setFrame( const std::string & frame, bool base )
[921]408{
[847]409 MFrequency::Types mdr;
410 if (!MFrequency::getType(mdr, frame)) {
411 Int a,b;const uInt* c;
412 const String* valid = MFrequency::allMyTypes(a, b, c);
[866]413 Vector<String> ftypes(IPosition(1,a), valid);
414 ostringstream oss;
415 oss << String("Please specify a legal frequency type. Types are\n");
416 oss << ftypes;
417 String msg(oss);
418 throw(AipsError(msg));
[847]419 } else {
[921]420 if (base)
421 table_.rwKeywordSet().define("BASEFRAME", frame);
422 else
423 table_.rwKeywordSet().define("FRAME", frame);
[847]424 }
425}
426
[1360]427void STFrequencies::setDoppler( const std::string & doppler )
[866]428{
429 MDoppler::Types mdt;
430 if (!MDoppler::getType(mdt, doppler)) {
431 Int a,b;const uInt* c;
432 const String* valid = MDoppler::allMyTypes(a, b, c);
433 Vector<String> ftypes(IPosition(1,a), valid);
434 ostringstream oss;
435 oss << String("Please specify a legal doppler type. Types are\n");
436 oss << ftypes;
437 String msg(oss);
438 throw(AipsError(msg));
439 } else {
440 table_.rwKeywordSet().define("DOPPLER", doppler);
441 }
442}
[849]443
[1694]444void STFrequencies::shiftRefPix(int npix, uInt id)
[1360]445{
[2243]446 Table t = table_(table_.col("ID") == Int(id), 1 );
[1360]447 if ( t.nrow() == 0 ) throw(AipsError("Selected Illegal frequency id"));
448 ScalarColumn<Double> tcol(t, "REFPIX");
449 tcol.put(0, tcol(0)+Double(npix));
450}
[866]451
[806]452} // namespace
Note: See TracBrowser for help on using the repository browser.