source: branches/alma/src/STFrequencies.cpp@ 1553

Last change on this file since 1553 was 1446, checked in by TakTsutsumi, 16 years ago

Merged recent updates (since 2007) from nrao-asap

File size: 12.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
130void STFrequencies::setEntry( Double refpix, Double refval, Double inc, uInt id )
131{
132 Table t = table_(table_.col("ID") == Int(id) );
133 if (t.nrow() == 0 ) {
134 throw(AipsError("STFrequencies::getEntry - freqID out of range"));
135 }
136 for ( uInt i = 0 ; i < table_.nrow() ; i++ ) {
137 uInt fid ;
138 idCol_.get( i, fid ) ;
139 if ( fid == id ) {
140 refpixCol_.put( i, refpix ) ;
141 refvalCol_.put( i, refval ) ;
142 incrCol_.put( i, inc ) ;
143 }
144 }
145}
146
147SpectralCoordinate STFrequencies::getSpectralCoordinate( uInt id ) const
148{
149 Table t = table_(table_.col("ID") == Int(id) );
150
151 if (t.nrow() == 0 ) {
152 throw(AipsError("STFrequencies::getSpectralCoordinate - ID out of range"));
153 }
154
155 // get the data
156 ROTableRow row(t);
157 // get first row - there should only be one matching id
158 const TableRecord& rec = row.get(0);
159
160 return SpectralCoordinate( getFrame(true), rec.asDouble("REFVAL"),
161 rec.asDouble("INCREMENT"),
162 rec.asDouble("REFPIX"));
163}
164
165/**
166SpectralCoordinate
167 asap::STFrequencies::getSpectralCoordinate( const MDirection& md,
168 const MPosition& mp,
169 const MEpoch& me,
170 Double restfreq, uInt id ) const
171**/
172SpectralCoordinate
173 asap::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,1.0/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 pushLog("WARNING: Frequency type unknown assuming TOPO");
268 mft = MFrequency::TOPO;
269 }
270
271 return mft;
272}
273
274std::string STFrequencies::getFrameString( bool base ) const
275{
276 if ( base ) return table_.keywordSet().asString("BASEFRAME");
277 else return table_.keywordSet().asString("FRAME");
278}
279
280std::string STFrequencies::getUnitString( ) const
281{
282 return table_.keywordSet().asString("UNIT");
283}
284
285Unit STFrequencies::getUnit( ) const
286{
287 return Unit(table_.keywordSet().asString("UNIT"));
288}
289
290std::string STFrequencies::getDopplerString( ) const
291{
292 return table_.keywordSet().asString("DOPPLER");
293}
294
295MDoppler::Types STFrequencies::getDoppler( ) const
296{
297 String dpl = table_.keywordSet().asString("DOPPLER");
298
299 // Create SpectralCoordinate (units Hz)
300 MDoppler::Types mdt;
301 if (!MDoppler::getType(mdt, dpl)) {
302 throw(AipsError("Doppler type unknown"));
303 }
304 return mdt;
305}
306
307std::string STFrequencies::print( int id, Bool strip ) const
308{
309 Table t;
310 ostringstream oss;
311 if ( id < 0 ) t = table_;
312 else t = table_(table_.col("ID") == Int(id) );
313 ROTableRow row(t);
314 for (uInt i=0; i<t.nrow(); ++i) {
315 const TableRecord& rec = row.get(i);
316 oss << setw(8)
317 << t.keywordSet().asString("FRAME") << setw(16) << setprecision(8)
318 << rec.asDouble("REFVAL") << setw(7)
319 << rec.asDouble("REFPIX")
320 << setw(12)
321 << rec.asDouble("INCREMENT");
322 }
323 String outstr(oss);
324 if (strip) {
325 int f = outstr.find_first_not_of(' ');
326 int l = outstr.find_last_not_of(' ', outstr.size());
327 if (f < 0) {
328 f = 0;
329 }
330 if ( l < f || l < f ) {
331 l = outstr.size();
332 }
333 return outstr.substr(f,l);
334 }
335 return outstr;
336}
337
338float STFrequencies::getRefFreq( uInt id, uInt channel )
339{
340 Table t = table_(table_.col("ID") == Int(id) );
341 if ( t.nrow() == 0 ) throw(AipsError("Selected Illegal frequency id"));
342 ROTableRow row(t);
343 const TableRecord& rec = row.get(0);
344 return (Double(channel/2) - rec.asDouble("REFPIX"))
345 * rec.asDouble("INCREMENT") + rec.asDouble("REFVAL");
346}
347
348bool STFrequencies::conformant( const STFrequencies& other ) const
349{
350 const Record& r = table_.keywordSet();
351 const Record& ro = other.table_.keywordSet();
352 return ( r.asString("FRAME") == ro.asString("FRAME") &&
353 r.asString("EQUINOX") == ro.asString("EQUINOX") &&
354 r.asString("UNIT") == ro.asString("UNIT") &&
355 r.asString("DOPPLER") == ro.asString("DOPPLER")
356 );
357}
358
359std::vector< std::string > STFrequencies::getInfo( ) const
360{
361 const Record& r = table_.keywordSet();
362 std::vector<std::string> out;
363 out.push_back(r.asString("UNIT"));
364 out.push_back(r.asString("FRAME"));
365 out.push_back(r.asString("DOPPLER"));
366 return out;
367}
368
369void STFrequencies::setInfo( const std::vector< std::string >& theinfo )
370{
371 if ( theinfo.size() != 3 ) throw(AipsError("setInfo needs three parameters"));
372 try {
373 setUnit(theinfo[0]);
374 setFrame(theinfo[1]);
375 setDoppler(theinfo[2]);
376 } catch (AipsError& e) {
377 throw(e);
378 }
379}
380
381void STFrequencies::setUnit( const std::string & unit )
382{
383 if (unit == "" || unit == "pixel" || unit == "channel" ) {
384 table_.rwKeywordSet().define("UNIT", "");
385 } else {
386 Unit u(unit);
387 if ( u == Unit("km/s") || u == Unit("Hz") )
388 table_.rwKeywordSet().define("UNIT", unit);
389 else {
390 throw(AipsError("Illegal spectral unit."));
391 }
392 }
393}
394
395void STFrequencies::setFrame(MFrequency::Types frame, bool base )
396{
397 String f = MFrequency::showType(frame);
398 if (base)
399 table_.rwKeywordSet().define("BASEFRAME", f);
400 else
401 table_.rwKeywordSet().define("FRAME", f);
402
403}
404
405void STFrequencies::setFrame( const std::string & frame, bool base )
406{
407 MFrequency::Types mdr;
408 if (!MFrequency::getType(mdr, frame)) {
409 Int a,b;const uInt* c;
410 const String* valid = MFrequency::allMyTypes(a, b, c);
411 Vector<String> ftypes(IPosition(1,a), valid);
412 ostringstream oss;
413 oss << String("Please specify a legal frequency type. Types are\n");
414 oss << ftypes;
415 String msg(oss);
416 throw(AipsError(msg));
417 } else {
418 if (base)
419 table_.rwKeywordSet().define("BASEFRAME", frame);
420 else
421 table_.rwKeywordSet().define("FRAME", frame);
422 }
423}
424
425void STFrequencies::setDoppler( const std::string & doppler )
426{
427 MDoppler::Types mdt;
428 if (!MDoppler::getType(mdt, doppler)) {
429 Int a,b;const uInt* c;
430 const String* valid = MDoppler::allMyTypes(a, b, c);
431 Vector<String> ftypes(IPosition(1,a), valid);
432 ostringstream oss;
433 oss << String("Please specify a legal doppler type. Types are\n");
434 oss << ftypes;
435 String msg(oss);
436 throw(AipsError(msg));
437 } else {
438 table_.rwKeywordSet().define("DOPPLER", doppler);
439 }
440}
441
442void STFrequencies::shiftRefPix(int npix, uInt id)
443{
444 Table t = table_(table_.col("ID") == Int(id) );
445 if ( t.nrow() == 0 ) throw(AipsError("Selected Illegal frequency id"));
446 ScalarColumn<Double> tcol(t, "REFPIX");
447 tcol.put(0, tcol(0)+Double(npix));
448}
449
450} // namespace
Note: See TracBrowser for help on using the repository browser.