source: branches/asap-3.x/src/STFrequencies.cpp@ 2676

Last change on this file since 2676 was 1694, checked in by Malte Marquarding, 15 years ago

Ticket #172: inverse scaling factor was used in coordinate resampling. Also added no of channels to each IF in scantable.summary and made all get{IF/Beam/Pol}No functions const

File size: 11.2 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
43STFrequencies::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 & 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 return SpectralCoordinate( getFrame(true), rec.asDouble("REFVAL"),
143 rec.asDouble("INCREMENT"),
144 rec.asDouble("REFPIX"));
145}
146
147SpectralCoordinate
148 STFrequencies::getSpectralCoordinate( const MDirection& md,
149 const MPosition& mp,
150 const MEpoch& me,
151 Double restfreq, uInt id ) const
152{
153 SpectralCoordinate spc = getSpectralCoordinate(id);
154 spc.setRestFrequency(restfreq, True);
155 if ( !spc.setReferenceConversion(getFrame(), me, mp, md) ) {
156 throw(AipsError("Couldn't convert frequency frame."));
157 }
158 String unitstr = getUnitString();
159 if ( !unitstr.empty() ) {
160 Unit unitu(unitstr);
161 if ( unitu == Unit("Hz") ) {
162 Vector<String> wau(1); wau = unitu.getName();
163 spc.setWorldAxisUnits(wau);
164 } else {
165 spc.setVelocity(unitstr, getDoppler());
166 }
167 }
168 return spc;
169}
170
171
172void STFrequencies::rescale( Float factor, const std::string& mode )
173{
174 TableRow row(table_);
175 TableRecord& outrec = row.record();
176 RecordFieldPtr<Double> rv(outrec, "REFVAL");
177 RecordFieldPtr<Double> rp(outrec, "REFPIX");
178 RecordFieldPtr<Double> inc(outrec, "INCREMENT");
179 for (uInt i=0; i<table_.nrow(); ++i) {
180
181 const TableRecord& rec = row.get(i);
182
183 SpectralCoordinate sc ( getFrame(true), rec.asDouble("REFVAL"),
184 rec.asDouble("INCREMENT"), rec.asDouble("REFPIX") );
185
186 SpectralCoordinate scout;
187 if (mode == "BIN") {
188 scout = binCsys(sc, Int(factor));
189 } else if (mode == "RESAMPLE") {
190 scout = resampleCsys(sc, factor);
191 }
192 *rv = scout.referenceValue()[0];
193 *rp = scout.referencePixel()[0];
194 *inc = scout.increment()[0];
195 row.put(i);
196 }
197}
198
199SpectralCoordinate STFrequencies::binCsys(const SpectralCoordinate& sc,
200 Int factor)
201{
202 CoordinateSystem csys;
203 csys.addCoordinate(sc);
204 IPosition factors(1, factor);
205 CoordinateSystem binnedcs =
206 CoordinateUtil::makeBinnedCoordinateSystem(factors, csys, False);
207 return binnedcs.spectralCoordinate(0);
208}
209
210SpectralCoordinate STFrequencies::resampleCsys(const SpectralCoordinate& sc,
211 Float width)
212{
213 Vector<Float> offset(1,0.0);
214 Vector<Float> factors(1,width);
215 Vector<Int> newshape;
216 CoordinateSystem csys;
217 csys.addCoordinate(sc);
218 CoordinateSystem csys2 = csys.subImage(offset, factors, newshape);
219 return csys2.spectralCoordinate(0);
220}
221
222
223MFrequency::Types STFrequencies::getFrame(bool base) const
224{
225 // get the ref frame
226 String rf;
227 if ( base )
228 rf = table_.keywordSet().asString("BASEFRAME");
229 else
230 rf = table_.keywordSet().asString("FRAME");
231
232 // Create SpectralCoordinate (units Hz)
233 MFrequency::Types mft;
234 if (!MFrequency::getType(mft, rf)) {
235 ostringstream oss;
236 pushLog("WARNING: Frequency type unknown assuming TOPO");
237 mft = MFrequency::TOPO;
238 }
239
240 return mft;
241}
242
243std::string STFrequencies::getFrameString( bool base ) const
244{
245 if ( base ) return table_.keywordSet().asString("BASEFRAME");
246 else return table_.keywordSet().asString("FRAME");
247}
248
249std::string STFrequencies::getUnitString( ) const
250{
251 return table_.keywordSet().asString("UNIT");
252}
253
254Unit STFrequencies::getUnit( ) const
255{
256 return Unit(table_.keywordSet().asString("UNIT"));
257}
258
259std::string STFrequencies::getDopplerString( ) const
260{
261 return table_.keywordSet().asString("DOPPLER");
262}
263
264MDoppler::Types STFrequencies::getDoppler( ) const
265{
266 String dpl = table_.keywordSet().asString("DOPPLER");
267
268 // Create SpectralCoordinate (units Hz)
269 MDoppler::Types mdt;
270 if (!MDoppler::getType(mdt, dpl)) {
271 throw(AipsError("Doppler type unknown"));
272 }
273 return mdt;
274}
275
276std::string STFrequencies::print( int id, Bool strip ) const
277{
278 Table t;
279 ostringstream oss;
280 if ( id < 0 ) t = table_;
281 else t = table_(table_.col("ID") == Int(id) );
282 ROTableRow row(t);
283 for (uInt i=0; i<t.nrow(); ++i) {
284 const TableRecord& rec = row.get(i);
285 oss << setw(8)
286 << t.keywordSet().asString("FRAME") << setw(16) << setprecision(8)
287 << rec.asDouble("REFVAL") << setw(7)
288 << rec.asDouble("REFPIX")
289 << setw(15)
290 << rec.asDouble("INCREMENT");
291 }
292 String outstr(oss);
293 if (strip) {
294 int f = outstr.find_first_not_of(' ');
295 int l = outstr.find_last_not_of(' ', outstr.size());
296 if (f < 0) {
297 f = 0;
298 }
299 if ( l < f || l < f ) {
300 l = outstr.size();
301 }
302 return outstr.substr(f,l);
303 }
304 return outstr;
305}
306
307float STFrequencies::getRefFreq( uInt id, uInt channel )
308{
309 Table t = table_(table_.col("ID") == Int(id) );
310 if ( t.nrow() == 0 ) throw(AipsError("Selected Illegal frequency id"));
311 ROTableRow row(t);
312 const TableRecord& rec = row.get(0);
313 return (Double(channel/2) - rec.asDouble("REFPIX"))
314 * rec.asDouble("INCREMENT") + rec.asDouble("REFVAL");
315}
316
317bool STFrequencies::conformant( const STFrequencies& other ) const
318{
319 const Record& r = table_.keywordSet();
320 const Record& ro = other.table_.keywordSet();
321 return ( r.asString("FRAME") == ro.asString("FRAME") &&
322 r.asString("EQUINOX") == ro.asString("EQUINOX") &&
323 r.asString("UNIT") == ro.asString("UNIT") &&
324 r.asString("DOPPLER") == ro.asString("DOPPLER")
325 );
326}
327
328std::vector< std::string > STFrequencies::getInfo( ) const
329{
330 const Record& r = table_.keywordSet();
331 std::vector<std::string> out;
332 out.push_back(r.asString("UNIT"));
333 out.push_back(r.asString("FRAME"));
334 out.push_back(r.asString("DOPPLER"));
335 return out;
336}
337
338void STFrequencies::setInfo( const std::vector< std::string >& theinfo )
339{
340 if ( theinfo.size() != 3 ) throw(AipsError("setInfo needs three parameters"));
341 try {
342 setUnit(theinfo[0]);
343 setFrame(theinfo[1]);
344 setDoppler(theinfo[2]);
345 } catch (AipsError& e) {
346 throw(e);
347 }
348}
349
350void STFrequencies::setUnit( const std::string & unit )
351{
352 if (unit == "" || unit == "pixel" || unit == "channel" ) {
353 table_.rwKeywordSet().define("UNIT", "");
354 } else {
355 Unit u(unit);
356 if ( u == Unit("km/s") || u == Unit("Hz") )
357 table_.rwKeywordSet().define("UNIT", unit);
358 else {
359 throw(AipsError("Illegal spectral unit."));
360 }
361 }
362}
363
364void STFrequencies::setFrame(MFrequency::Types frame, bool base )
365{
366 String f = MFrequency::showType(frame);
367 if (base)
368 table_.rwKeywordSet().define("BASEFRAME", f);
369 else
370 table_.rwKeywordSet().define("FRAME", f);
371
372}
373
374void STFrequencies::setFrame( const std::string & frame, bool base )
375{
376 MFrequency::Types mdr;
377 if (!MFrequency::getType(mdr, frame)) {
378 Int a,b;const uInt* c;
379 const String* valid = MFrequency::allMyTypes(a, b, c);
380 Vector<String> ftypes(IPosition(1,a), valid);
381 ostringstream oss;
382 oss << String("Please specify a legal frequency type. Types are\n");
383 oss << ftypes;
384 String msg(oss);
385 throw(AipsError(msg));
386 } else {
387 if (base)
388 table_.rwKeywordSet().define("BASEFRAME", frame);
389 else
390 table_.rwKeywordSet().define("FRAME", frame);
391 }
392}
393
394void STFrequencies::setDoppler( const std::string & doppler )
395{
396 MDoppler::Types mdt;
397 if (!MDoppler::getType(mdt, doppler)) {
398 Int a,b;const uInt* c;
399 const String* valid = MDoppler::allMyTypes(a, b, c);
400 Vector<String> ftypes(IPosition(1,a), valid);
401 ostringstream oss;
402 oss << String("Please specify a legal doppler type. Types are\n");
403 oss << ftypes;
404 String msg(oss);
405 throw(AipsError(msg));
406 } else {
407 table_.rwKeywordSet().define("DOPPLER", doppler);
408 }
409}
410
411void STFrequencies::shiftRefPix(int npix, uInt id)
412{
413 Table t = table_(table_.col("ID") == Int(id) );
414 if ( t.nrow() == 0 ) throw(AipsError("Selected Illegal frequency id"));
415 ScalarColumn<Double> tcol(t, "REFPIX");
416 tcol.put(0, tcol(0)+Double(npix));
417}
418
419} // namespace
Note: See TracBrowser for help on using the repository browser.