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 |
|
---|
31 | using namespace casa;
|
---|
32 |
|
---|
33 | namespace asap {
|
---|
34 |
|
---|
35 | const String STFrequencies::name_ = "FREQUENCIES";
|
---|
36 |
|
---|
37 | STFrequencies::STFrequencies(const Scantable& parent) :
|
---|
38 | STSubTable(parent, name_)
|
---|
39 | {
|
---|
40 | setup();
|
---|
41 | }
|
---|
42 |
|
---|
43 | 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 |
|
---|
52 | STFrequencies::~STFrequencies()
|
---|
53 | {
|
---|
54 | }
|
---|
55 |
|
---|
56 | STFrequencies & 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 |
|
---|
67 | void 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 |
|
---|
86 | uInt 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 |
|
---|
115 | void 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 |
|
---|
130 | void 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 |
|
---|
147 | SpectralCoordinate 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 | /**
|
---|
165 | SpectralCoordinate
|
---|
166 | STFrequencies::getSpectralCoordinate( const MDirection& md,
|
---|
167 | const MPosition& mp,
|
---|
168 | const MEpoch& me,
|
---|
169 | Double restfreq, uInt id ) const
|
---|
170 | **/
|
---|
171 | SpectralCoordinate
|
---|
172 | asap::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 |
|
---|
202 | void 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 |
|
---|
229 | SpectralCoordinate 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 |
|
---|
240 | SpectralCoordinate STFrequencies::resampleCsys(const SpectralCoordinate& sc,
|
---|
241 | Float width)
|
---|
242 | {
|
---|
243 | Vector<Float> offset(1,0.0);
|
---|
244 | Vector<Float> factors(1,1.0/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 |
|
---|
253 | MFrequency::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 |
|
---|
273 | std::string STFrequencies::getFrameString( bool base ) const
|
---|
274 | {
|
---|
275 | if ( base ) return table_.keywordSet().asString("BASEFRAME");
|
---|
276 | else return table_.keywordSet().asString("FRAME");
|
---|
277 | }
|
---|
278 |
|
---|
279 | std::string STFrequencies::getUnitString( ) const
|
---|
280 | {
|
---|
281 | return table_.keywordSet().asString("UNIT");
|
---|
282 | }
|
---|
283 |
|
---|
284 | Unit STFrequencies::getUnit( ) const
|
---|
285 | {
|
---|
286 | return Unit(table_.keywordSet().asString("UNIT"));
|
---|
287 | }
|
---|
288 |
|
---|
289 | std::string STFrequencies::getDopplerString( ) const
|
---|
290 | {
|
---|
291 | return table_.keywordSet().asString("DOPPLER");
|
---|
292 | }
|
---|
293 |
|
---|
294 | MDoppler::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 |
|
---|
306 | std::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(12)
|
---|
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 |
|
---|
337 | float 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 |
|
---|
347 | bool 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 |
|
---|
358 | std::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 |
|
---|
368 | void 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 |
|
---|
380 | void 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 |
|
---|
394 | void 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 |
|
---|
404 | void 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 |
|
---|
424 | void 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 |
|
---|
441 | void 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
|
---|