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

Last change on this file since 2170 was 1757, checked in by Kana Sugimoto, 15 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

Note you also need to update casa/code/atnf.


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