source: trunk/src/STFrequencies.cpp@ 909

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

fixed bug in print(): was using column instead of keyword to retrieve frame

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