source: trunk/src/STFrequencies.cpp@ 2748

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

Ticket #199: Excised Logger::pushLog; now everything is using LogIO

File size: 12.3 KB
Line 
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>
17#include <casa/Logging/LogIO.h>
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
36const String STFrequencies::name_ = "FREQUENCIES";
37
38STFrequencies::STFrequencies(const Scantable& parent) :
39 STSubTable(parent, name_)
40{
41 setup();
42}
43
44STFrequencies::STFrequencies( casa::Table tab ) :
45 STSubTable(tab, name_)
46{
47 refpixCol_.attach(table_,"REFPIX");
48 refvalCol_.attach(table_,"REFVAL");
49 incrCol_.attach(table_,"INCREMENT");
50
51}
52
53STFrequencies::~STFrequencies()
54{
55}
56
57STFrequencies & STFrequencies::operator=( const STFrequencies & other )
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
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
75 table_.rwKeywordSet().define("FRAME", String("TOPO"));
76 table_.rwKeywordSet().define("BASEFRAME", String("TOPO"));
77 table_.rwKeywordSet().define("EQUINOX",String( "J2000"));
78 table_.rwKeywordSet().define("UNIT", String(""));
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)
91 && near(table_.col("REFPIX"), refpix)
92 && near(table_.col("INCREMENT"), inc), 1 );
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
115
116void STFrequencies::getEntry( Double& refpix, Double& refval, Double& inc,
117 uInt id )
118{
119 Table t = table_(table_.col("ID") == Int(id), 1 );
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
131void STFrequencies::setEntry( Double refpix, Double refval, Double inc, uInt id )
132{
133 Table t = table_(table_.col("ID") == Int(id), 1 );
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
148SpectralCoordinate STFrequencies::getSpectralCoordinate( uInt id ) const
149{
150 Table t = table_(table_.col("ID") == Int(id), 1 );
151
152 if (t.nrow() == 0 ) {
153 throw(AipsError("STFrequencies::getSpectralCoordinate - ID out of range"));
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);
160 return SpectralCoordinate( getFrame(true), rec.asDouble("REFVAL"),
161 rec.asDouble("INCREMENT"),
162 rec.asDouble("REFPIX"));
163}
164
165/**
166SpectralCoordinate
167 STFrequencies::getSpectralCoordinate( const MDirection& md,
168 const MPosition& mp,
169 const MEpoch& me,
170 Double restfreq, uInt id ) const
171**/
172SpectralCoordinate
173 STFrequencies::getSpectralCoordinate( const MDirection& md,
174 const MPosition& mp,
175 const MEpoch& me,
176 Vector<Double> restfreq, uInt id ) const
177{
178 SpectralCoordinate spc = getSpectralCoordinate(id);
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);
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") ) {
193 Vector<String> wau(1); wau = unitu.getName();
194 spc.setWorldAxisUnits(wau);
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{
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
214 SpectralCoordinate sc ( getFrame(true), rec.asDouble("REFVAL"),
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);
245 Vector<Float> factors(1,width);
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
254MFrequency::Types STFrequencies::getFrame(bool base) const
255{
256 // get the ref frame
257 String rf;
258 if ( base )
259 rf = table_.keywordSet().asString("BASEFRAME");
260 else
261 rf = table_.keywordSet().asString("FRAME");
262
263 // Create SpectralCoordinate (units Hz)
264 MFrequency::Types mft;
265 if (!MFrequency::getType(mft, rf)) {
266 ostringstream oss;
267 LogIO os( casa::LogOrigin( "STFrequencies", "getFrame") );
268 os << LogIO::WARN << "WARNING: Frequency type unknown assuming TOPO"
269 << LogIO::POST;
270 mft = MFrequency::TOPO;
271 }
272
273 return mft;
274}
275
276std::string STFrequencies::getFrameString( bool base ) const
277{
278 if ( base ) return table_.keywordSet().asString("BASEFRAME");
279 else return table_.keywordSet().asString("FRAME");
280}
281
282std::string STFrequencies::getUnitString( ) const
283{
284 return table_.keywordSet().asString("UNIT");
285}
286
287Unit STFrequencies::getUnit( ) const
288{
289 return Unit(table_.keywordSet().asString("UNIT"));
290}
291
292std::string STFrequencies::getDopplerString( ) const
293{
294 return table_.keywordSet().asString("DOPPLER");
295}
296
297MDoppler::Types STFrequencies::getDoppler( ) const
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
309std::string STFrequencies::print( int id, Bool strip ) const
310{
311 Table t;
312 ostringstream oss;
313 if ( id < 0 ) t = table_;
314 else t = table_(table_.col("ID") == Int(id), 1 );
315 ROTableRow row(t);
316 for (uInt i=0; i<t.nrow(); ++i) {
317 const TableRecord& rec = row.get(i);
318 oss << setw(8)
319 << t.keywordSet().asString("FRAME") << setw(16) << setprecision(8)
320 << rec.asDouble("REFVAL") << setw(7)
321 << rec.asDouble("REFPIX")
322 << setw(15)
323 << rec.asDouble("INCREMENT");
324 }
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());
329 if (f < 0) {
330 f = 0;
331 }
332 if ( l < f || l < f ) {
333 l = outstr.size();
334 }
335 return outstr.substr(f,l);
336 }
337 return outstr;
338}
339
340float STFrequencies::getRefFreq( uInt id, uInt channel )
341{
342 Table t = table_(table_.col("ID") == Int(id), 1 );
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
350bool STFrequencies::conformant( const STFrequencies& other ) const
351{
352 const Record& r = table_.keywordSet();
353 const Record& ro = other.table_.keywordSet();
354 return ( r.asString("FRAME") == ro.asString("FRAME") &&
355 r.asString("EQUINOX") == ro.asString("EQUINOX") &&
356 r.asString("UNIT") == ro.asString("UNIT") &&
357 r.asString("DOPPLER") == ro.asString("DOPPLER")
358 );
359}
360
361std::vector< std::string > STFrequencies::getInfo( ) const
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"));
368 return out;
369}
370
371void STFrequencies::setInfo( const std::vector< std::string >& theinfo )
372{
373 if ( theinfo.size() != 3 ) throw(AipsError("setInfo needs three parameters"));
374 try {
375 setUnit(theinfo[0]);
376 setFrame(theinfo[1]);
377 setDoppler(theinfo[2]);
378 } catch (AipsError& e) {
379 throw(e);
380 }
381}
382
383void STFrequencies::setUnit( const std::string & unit )
384{
385 if (unit == "" || unit == "pixel" || unit == "channel" ) {
386 table_.rwKeywordSet().define("UNIT", "");
387 } else {
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 }
394 }
395}
396
397void STFrequencies::setFrame(MFrequency::Types frame, bool base )
398{
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
407void STFrequencies::setFrame( const std::string & frame, bool base )
408{
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);
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));
419 } else {
420 if (base)
421 table_.rwKeywordSet().define("BASEFRAME", frame);
422 else
423 table_.rwKeywordSet().define("FRAME", frame);
424 }
425}
426
427void STFrequencies::setDoppler( const std::string & doppler )
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}
443
444void STFrequencies::shiftRefPix(int npix, uInt id)
445{
446 Table t = table_(table_.col("ID") == Int(id), 1 );
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}
451
452} // namespace
Note: See TracBrowser for help on using the repository browser.