source: trunk/src/STFrequencies.cpp@ 1014

Last change on this file since 1014 was 941, checked in by mar637, 19 years ago

fixed the output alignment of summary
File size: 10.9 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(""));
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;
228 if ( base )
229 rf = table_.keywordSet().asString("BASEFRAME");
230 else
231 rf = table_.keywordSet().asString("FRAME");
232
233 // Create SpectralCoordinate (units Hz)
234 MFrequency::Types mft;
235 if (!MFrequency::getType(mft, rf)) {
236 ostringstream oss;
237 pushLog("WARNING: Frequency type unknown assuming TOPO");
238 mft = MFrequency::TOPO;
239 }
240
241 return mft;
242}
243
244std::string asap::STFrequencies::getFrameString( bool base ) const
245{
246 if ( base ) return table_.keywordSet().asString("BASEFRAME");
247 else return table_.keywordSet().asString("FRAME");
248}
249
250std::string asap::STFrequencies::getUnitString( ) const
251{
252 return table_.keywordSet().asString("UNIT");
253}
254
255Unit asap::STFrequencies::getUnit( ) const
256{
257 return Unit(table_.keywordSet().asString("UNIT"));
258}
259
260std::string asap::STFrequencies::getDopplerString( ) const
261{
262 return table_.keywordSet().asString("DOPPLER");
263}
264
265MDoppler::Types asap::STFrequencies::getDoppler( ) const
266{
267 String dpl = table_.keywordSet().asString("DOPPLER");
268
269 // Create SpectralCoordinate (units Hz)
270 MDoppler::Types mdt;
271 if (!MDoppler::getType(mdt, dpl)) {
272 throw(AipsError("Doppler type unknown"));
273 }
274 return mdt;
275}
276
277std::string asap::STFrequencies::print( int id )
278{
279 Table t;
280 ostringstream oss;
281 if ( id < 0 ) t = table_;
282 else t = table_(table_.col("ID") == Int(id) );
283 ROTableRow row(t);
284 for (uInt i=0; i<t.nrow(); ++i) {
285 const TableRecord& rec = row.get(i);
286 oss << setw(8)
287 << t.keywordSet().asString("FRAME") << setw(16) << setprecision(8)
288 << rec.asDouble("REFVAL") << setw(7)
289 << rec.asDouble("REFPIX") << setw(12)
290 << rec.asDouble("INCREMENT") << endl;
291 }
292 return String(oss);
293}
294
295float STFrequencies::getRefFreq( uInt id, uInt channel )
296{
297 Table t = table_(table_.col("ID") == Int(id) );
298 if ( t.nrow() == 0 ) throw(AipsError("Selected Illegal frequency id"));
299 ROTableRow row(t);
300 const TableRecord& rec = row.get(0);
301 return (Double(channel/2) - rec.asDouble("REFPIX"))
302 * rec.asDouble("INCREMENT") + rec.asDouble("REFVAL");
303}
304
305bool asap::STFrequencies::conformant( const STFrequencies& other ) const
306{
307 const Record& r = table_.keywordSet();
308 const Record& ro = other.table_.keywordSet();
309 return ( r.asString("FRAME") == ro.asString("FRAME") &&
310 r.asString("EQUINOX") == ro.asString("EQUINOX") &&
311 r.asString("UNIT") == ro.asString("UNIT") &&
312 r.asString("DOPPLER") == ro.asString("DOPPLER")
313 );
314}
315
316std::vector< std::string > asap::STFrequencies::getInfo( ) const
317{
318 const Record& r = table_.keywordSet();
319 std::vector<std::string> out;
320 out.push_back(r.asString("UNIT"));
321 out.push_back(r.asString("FRAME"));
322 out.push_back(r.asString("DOPPLER"));
323 return out;
324}
325
326void asap::STFrequencies::setInfo( const std::vector< std::string >& theinfo )
327{
328 if ( theinfo.size() != 3 ) throw(AipsError("setInfo needs three parameters"));
329 try {
330 setUnit(theinfo[0]);
331 setFrame(theinfo[1]);
332 setDoppler(theinfo[2]);
333 } catch (AipsError& e) {
334 throw(e);
335 }
336}
337
338void asap::STFrequencies::setUnit( const std::string & unit )
339{
340 if (unit == "" || unit == "pixel" || unit == "channel" ) {
341 table_.rwKeywordSet().define("UNIT", "");
342 } else {
343 Unit u(unit);
344 if ( u == Unit("km/s") || u == Unit("Hz") )
345 table_.rwKeywordSet().define("UNIT", unit);
346 else {
347 throw(AipsError("Illegal spectral unit."));
348 }
349 }
350}
351
352void asap::STFrequencies::setFrame(MFrequency::Types frame, bool base )
353{
354 String f = MFrequency::showType(frame);
355 if (base)
356 table_.rwKeywordSet().define("BASEFRAME", f);
357 else
358 table_.rwKeywordSet().define("FRAME", f);
359
360}
361
362void asap::STFrequencies::setFrame( const std::string & frame, bool base )
363{
364 MFrequency::Types mdr;
365 if (!MFrequency::getType(mdr, frame)) {
366 Int a,b;const uInt* c;
367 const String* valid = MFrequency::allMyTypes(a, b, c);
368 Vector<String> ftypes(IPosition(1,a), valid);
369 ostringstream oss;
370 oss << String("Please specify a legal frequency type. Types are\n");
371 oss << ftypes;
372 String msg(oss);
373 throw(AipsError(msg));
374 } else {
375 if (base)
376 table_.rwKeywordSet().define("BASEFRAME", frame);
377 else
378 table_.rwKeywordSet().define("FRAME", frame);
379 }
380}
381
382void asap::STFrequencies::setDoppler( const std::string & doppler )
383{
384 MDoppler::Types mdt;
385 if (!MDoppler::getType(mdt, doppler)) {
386 Int a,b;const uInt* c;
387 const String* valid = MDoppler::allMyTypes(a, b, c);
388 Vector<String> ftypes(IPosition(1,a), valid);
389 ostringstream oss;
390 oss << String("Please specify a legal doppler type. Types are\n");
391 oss << ftypes;
392 String msg(oss);
393 throw(AipsError(msg));
394 } else {
395 table_.rwKeywordSet().define("DOPPLER", doppler);
396 }
397}
398
399
400} // namespace
Note: See TracBrowser for help on using the repository browser.