source: branches/alma/src/Scantable.cpp@ 1644

Last change on this file since 1644 was 1616, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-729, CAS-1147

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Use LogIO instead of std::cout and std::cerr.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 48.8 KB
Line 
1//
2// C++ Implementation: Scantable
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2005
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <map>
13#include <fstream>
14
15#include <casa/aips.h>
16#include <casa/iostream.h>
17#include <casa/iomanip.h>
18#include <casa/OS/Path.h>
19#include <casa/OS/File.h>
20#include <casa/Arrays/Array.h>
21#include <casa/Arrays/ArrayMath.h>
22#include <casa/Arrays/MaskArrMath.h>
23#include <casa/Arrays/ArrayLogical.h>
24#include <casa/Arrays/ArrayAccessor.h>
25#include <casa/Arrays/Vector.h>
26#include <casa/Arrays/VectorSTLIterator.h>
27#include <casa/Arrays/Slice.h>
28#include <casa/BasicMath/Math.h>
29#include <casa/BasicSL/Constants.h>
30#include <casa/Quanta/MVAngle.h>
31#include <casa/Containers/RecordField.h>
32#include <casa/Utilities/GenSort.h>
33#include <casa/Logging/LogIO.h>
34
35#include <tables/Tables/TableParse.h>
36#include <tables/Tables/TableDesc.h>
37#include <tables/Tables/TableCopy.h>
38#include <tables/Tables/SetupNewTab.h>
39#include <tables/Tables/ScaColDesc.h>
40#include <tables/Tables/ArrColDesc.h>
41#include <tables/Tables/TableRow.h>
42#include <tables/Tables/TableVector.h>
43#include <tables/Tables/TableIter.h>
44
45#include <tables/Tables/ExprNode.h>
46#include <tables/Tables/TableRecord.h>
47#include <casa/Quanta/MVTime.h>
48#include <casa/Quanta/MVAngle.h>
49#include <measures/Measures/MeasRef.h>
50#include <measures/Measures/MeasTable.h>
51// needed to avoid error in .tcc
52#include <measures/Measures/MCDirection.h>
53//
54#include <measures/Measures/MDirection.h>
55#include <measures/Measures/MFrequency.h>
56#include <measures/Measures/MEpoch.h>
57#include <measures/TableMeasures/TableMeasRefDesc.h>
58#include <measures/TableMeasures/TableMeasValueDesc.h>
59#include <measures/TableMeasures/TableMeasDesc.h>
60#include <measures/TableMeasures/ScalarMeasColumn.h>
61#include <coordinates/Coordinates/CoordinateUtil.h>
62
63#include "Scantable.h"
64#include "STPolLinear.h"
65#include "STPolCircular.h"
66#include "STPolStokes.h"
67#include "STAttr.h"
68#include "MathUtils.h"
69
70using namespace casa;
71
72namespace asap {
73
74std::map<std::string, STPol::STPolFactory *> Scantable::factories_;
75
76void Scantable::initFactories() {
77 if ( factories_.empty() ) {
78 Scantable::factories_["linear"] = &STPolLinear::myFactory;
79 Scantable::factories_["circular"] = &STPolCircular::myFactory;
80 Scantable::factories_["stokes"] = &STPolStokes::myFactory;
81 }
82}
83
84Scantable::Scantable(Table::TableType ttype) :
85 type_(ttype)
86{
87 initFactories();
88 setupMainTable();
89 freqTable_ = STFrequencies(*this);
90 table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_.table());
91 weatherTable_ = STWeather(*this);
92 table_.rwKeywordSet().defineTable("WEATHER", weatherTable_.table());
93 focusTable_ = STFocus(*this);
94 table_.rwKeywordSet().defineTable("FOCUS", focusTable_.table());
95 tcalTable_ = STTcal(*this);
96 table_.rwKeywordSet().defineTable("TCAL", tcalTable_.table());
97 moleculeTable_ = STMolecules(*this);
98 table_.rwKeywordSet().defineTable("MOLECULES", moleculeTable_.table());
99 historyTable_ = STHistory(*this);
100 table_.rwKeywordSet().defineTable("HISTORY", historyTable_.table());
101 fitTable_ = STFit(*this);
102 table_.rwKeywordSet().defineTable("FIT", fitTable_.table());
103 originalTable_ = table_;
104 attach();
105}
106
107Scantable::Scantable(const std::string& name, Table::TableType ttype) :
108 type_(ttype)
109{
110 initFactories();
111 Table tab(name, Table::Update);
112 uInt version = tab.keywordSet().asuInt("VERSION");
113 if (version != version_) {
114 throw(AipsError("Unsupported version of ASAP file."));
115 }
116 if ( type_ == Table::Memory ) {
117 table_ = tab.copyToMemoryTable(generateName());
118 } else {
119 table_ = tab;
120 }
121
122 attachSubtables();
123 originalTable_ = table_;
124 attach();
125}
126
127Scantable::Scantable( const Scantable& other, bool clear )
128{
129 // with or without data
130 String newname = String(generateName());
131 type_ = other.table_.tableType();
132 if ( other.table_.tableType() == Table::Memory ) {
133 if ( clear ) {
134 table_ = TableCopy::makeEmptyMemoryTable(newname,
135 other.table_, True);
136 } else
137 table_ = other.table_.copyToMemoryTable(newname);
138 } else {
139 other.table_.deepCopy(newname, Table::New, False,
140 other.table_.endianFormat(),
141 Bool(clear));
142 table_ = Table(newname, Table::Update);
143 table_.markForDelete();
144 }
145 /// @todo reindex SCANNO, recompute nbeam, nif, npol
146 if ( clear ) copySubtables(other);
147 attachSubtables();
148 originalTable_ = table_;
149 attach();
150}
151
152void Scantable::copySubtables(const Scantable& other) {
153 Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
154 TableCopy::copyRows(t, other.freqTable_.table());
155 t = table_.rwKeywordSet().asTable("FOCUS");
156 TableCopy::copyRows(t, other.focusTable_.table());
157 t = table_.rwKeywordSet().asTable("WEATHER");
158 TableCopy::copyRows(t, other.weatherTable_.table());
159 t = table_.rwKeywordSet().asTable("TCAL");
160 TableCopy::copyRows(t, other.tcalTable_.table());
161 t = table_.rwKeywordSet().asTable("MOLECULES");
162 TableCopy::copyRows(t, other.moleculeTable_.table());
163 t = table_.rwKeywordSet().asTable("HISTORY");
164 TableCopy::copyRows(t, other.historyTable_.table());
165 t = table_.rwKeywordSet().asTable("FIT");
166 TableCopy::copyRows(t, other.fitTable_.table());
167}
168
169void Scantable::attachSubtables()
170{
171 freqTable_ = STFrequencies(table_);
172 focusTable_ = STFocus(table_);
173 weatherTable_ = STWeather(table_);
174 tcalTable_ = STTcal(table_);
175 moleculeTable_ = STMolecules(table_);
176 historyTable_ = STHistory(table_);
177 fitTable_ = STFit(table_);
178}
179
180Scantable::~Scantable()
181{
182 //cout << "~Scantable() " << this << endl;
183}
184
185void Scantable::setupMainTable()
186{
187 TableDesc td("", "1", TableDesc::Scratch);
188 td.comment() = "An ASAP Scantable";
189 td.rwKeywordSet().define("VERSION", uInt(version_));
190
191 // n Cycles
192 td.addColumn(ScalarColumnDesc<uInt>("SCANNO"));
193 // new index every nBeam x nIF x nPol
194 td.addColumn(ScalarColumnDesc<uInt>("CYCLENO"));
195
196 td.addColumn(ScalarColumnDesc<uInt>("BEAMNO"));
197 td.addColumn(ScalarColumnDesc<uInt>("IFNO"));
198 // linear, circular, stokes
199 td.rwKeywordSet().define("POLTYPE", String("linear"));
200 td.addColumn(ScalarColumnDesc<uInt>("POLNO"));
201
202 td.addColumn(ScalarColumnDesc<uInt>("FREQ_ID"));
203 td.addColumn(ScalarColumnDesc<uInt>("MOLECULE_ID"));
204 td.addColumn(ScalarColumnDesc<Int>("REFBEAMNO"));
205
206 td.addColumn(ScalarColumnDesc<Double>("TIME"));
207 TableMeasRefDesc measRef(MEpoch::UTC); // UTC as default
208 TableMeasValueDesc measVal(td, "TIME");
209 TableMeasDesc<MEpoch> mepochCol(measVal, measRef);
210 mepochCol.write(td);
211
212 td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
213
214 td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
215 // Type of source (on=0, off=1, other=-1)
216 ScalarColumnDesc<Int> stypeColumn("SRCTYPE");
217 stypeColumn.setDefault(Int(-1));
218 td.addColumn(stypeColumn);
219 td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
220
221 //The actual Data Vectors
222 td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
223 td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
224 td.addColumn(ArrayColumnDesc<Float>("TSYS"));
225
226 td.addColumn(ArrayColumnDesc<Double>("DIRECTION",
227 IPosition(1,2),
228 ColumnDesc::Direct));
229 TableMeasRefDesc mdirRef(MDirection::J2000); // default
230 TableMeasValueDesc tmvdMDir(td, "DIRECTION");
231 // the TableMeasDesc gives the column a type
232 TableMeasDesc<MDirection> mdirCol(tmvdMDir, mdirRef);
233 // a uder set table type e.g. GALCTIC, B1950 ...
234 td.rwKeywordSet().define("DIRECTIONREF", String("J2000"));
235 // writing create the measure column
236 mdirCol.write(td);
237 td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
238 td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
239 td.addColumn(ScalarColumnDesc<Float>("PARANGLE"));
240 td.addColumn(ScalarColumnDesc<Float>("OPACITY"));
241
242 td.addColumn(ScalarColumnDesc<uInt>("TCAL_ID"));
243 ScalarColumnDesc<Int> fitColumn("FIT_ID");
244 fitColumn.setDefault(Int(-1));
245 td.addColumn(fitColumn);
246
247 td.addColumn(ScalarColumnDesc<uInt>("FOCUS_ID"));
248 td.addColumn(ScalarColumnDesc<uInt>("WEATHER_ID"));
249
250 // columns which just get dragged along, as they aren't used in asap
251 td.addColumn(ScalarColumnDesc<Double>("SRCVELOCITY"));
252 td.addColumn(ArrayColumnDesc<Double>("SRCPROPERMOTION"));
253 td.addColumn(ArrayColumnDesc<Double>("SRCDIRECTION"));
254 td.addColumn(ArrayColumnDesc<Double>("SCANRATE"));
255
256 td.rwKeywordSet().define("OBSMODE", String(""));
257
258 // Now create Table SetUp from the description.
259 SetupNewTable aNewTab(generateName(), td, Table::Scratch);
260 table_ = Table(aNewTab, type_, 0);
261 originalTable_ = table_;
262}
263
264
265void Scantable::attach()
266{
267 timeCol_.attach(table_, "TIME");
268 srcnCol_.attach(table_, "SRCNAME");
269 srctCol_.attach(table_, "SRCTYPE");
270 specCol_.attach(table_, "SPECTRA");
271 flagsCol_.attach(table_, "FLAGTRA");
272 tsysCol_.attach(table_, "TSYS");
273 cycleCol_.attach(table_,"CYCLENO");
274 scanCol_.attach(table_, "SCANNO");
275 beamCol_.attach(table_, "BEAMNO");
276 ifCol_.attach(table_, "IFNO");
277 polCol_.attach(table_, "POLNO");
278 integrCol_.attach(table_, "INTERVAL");
279 azCol_.attach(table_, "AZIMUTH");
280 elCol_.attach(table_, "ELEVATION");
281 dirCol_.attach(table_, "DIRECTION");
282 paraCol_.attach(table_, "PARANGLE");
283 fldnCol_.attach(table_, "FIELDNAME");
284 rbeamCol_.attach(table_, "REFBEAMNO");
285
286 mfitidCol_.attach(table_,"FIT_ID");
287 mfreqidCol_.attach(table_, "FREQ_ID");
288 mtcalidCol_.attach(table_, "TCAL_ID");
289 mfocusidCol_.attach(table_, "FOCUS_ID");
290 mmolidCol_.attach(table_, "MOLECULE_ID");
291}
292
293void Scantable::setHeader(const STHeader& sdh)
294{
295 table_.rwKeywordSet().define("nIF", sdh.nif);
296 table_.rwKeywordSet().define("nBeam", sdh.nbeam);
297 table_.rwKeywordSet().define("nPol", sdh.npol);
298 table_.rwKeywordSet().define("nChan", sdh.nchan);
299 table_.rwKeywordSet().define("Observer", sdh.observer);
300 table_.rwKeywordSet().define("Project", sdh.project);
301 table_.rwKeywordSet().define("Obstype", sdh.obstype);
302 table_.rwKeywordSet().define("AntennaName", sdh.antennaname);
303 table_.rwKeywordSet().define("AntennaPosition", sdh.antennaposition);
304 table_.rwKeywordSet().define("Equinox", sdh.equinox);
305 table_.rwKeywordSet().define("FreqRefFrame", sdh.freqref);
306 table_.rwKeywordSet().define("FreqRefVal", sdh.reffreq);
307 table_.rwKeywordSet().define("Bandwidth", sdh.bandwidth);
308 table_.rwKeywordSet().define("UTC", sdh.utc);
309 table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit);
310 table_.rwKeywordSet().define("Epoch", sdh.epoch);
311 table_.rwKeywordSet().define("POLTYPE", sdh.poltype);
312}
313
314STHeader Scantable::getHeader() const
315{
316 STHeader sdh;
317 table_.keywordSet().get("nBeam",sdh.nbeam);
318 table_.keywordSet().get("nIF",sdh.nif);
319 table_.keywordSet().get("nPol",sdh.npol);
320 table_.keywordSet().get("nChan",sdh.nchan);
321 table_.keywordSet().get("Observer", sdh.observer);
322 table_.keywordSet().get("Project", sdh.project);
323 table_.keywordSet().get("Obstype", sdh.obstype);
324 table_.keywordSet().get("AntennaName", sdh.antennaname);
325 table_.keywordSet().get("AntennaPosition", sdh.antennaposition);
326 table_.keywordSet().get("Equinox", sdh.equinox);
327 table_.keywordSet().get("FreqRefFrame", sdh.freqref);
328 table_.keywordSet().get("FreqRefVal", sdh.reffreq);
329 table_.keywordSet().get("Bandwidth", sdh.bandwidth);
330 table_.keywordSet().get("UTC", sdh.utc);
331 table_.keywordSet().get("FluxUnit", sdh.fluxunit);
332 table_.keywordSet().get("Epoch", sdh.epoch);
333 table_.keywordSet().get("POLTYPE", sdh.poltype);
334 return sdh;
335}
336
337void Scantable::setSourceType( int stype )
338{
339 if ( stype < 0 || stype > 1 )
340 throw(AipsError("Illegal sourcetype."));
341 TableVector<Int> tabvec(table_, "SRCTYPE");
342 tabvec = Int(stype);
343}
344
345bool Scantable::conformant( const Scantable& other )
346{
347 return this->getHeader().conformant(other.getHeader());
348}
349
350
351
352std::string Scantable::formatSec(Double x) const
353{
354 Double xcop = x;
355 MVTime mvt(xcop/24./3600.); // make days
356
357 if (x < 59.95)
358 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
359 else if (x < 3599.95)
360 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_H,7)+" ";
361 else {
362 ostringstream oss;
363 oss << setw(2) << std::right << setprecision(1) << mvt.hour();
364 oss << ":" << mvt.string(MVTime::TIME_CLEAN_NO_H,7) << " ";
365 return String(oss);
366 }
367};
368
369std::string Scantable::formatDirection(const MDirection& md) const
370{
371 Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue();
372 Int prec = 7;
373
374 MVAngle mvLon(t[0]);
375 String sLon = mvLon.string(MVAngle::TIME,prec);
376 uInt tp = md.getRef().getType();
377 if (tp == MDirection::GALACTIC ||
378 tp == MDirection::SUPERGAL ) {
379 sLon = mvLon(0.0).string(MVAngle::ANGLE_CLEAN,prec);
380 }
381 MVAngle mvLat(t[1]);
382 String sLat = mvLat.string(MVAngle::ANGLE+MVAngle::DIG2,prec);
383 return sLon + String(" ") + sLat;
384}
385
386
387std::string Scantable::getFluxUnit() const
388{
389 return table_.keywordSet().asString("FluxUnit");
390}
391
392void Scantable::setFluxUnit(const std::string& unit)
393{
394 String tmp(unit);
395 Unit tU(tmp);
396 if (tU==Unit("K") || tU==Unit("Jy")) {
397 table_.rwKeywordSet().define(String("FluxUnit"), tmp);
398 } else {
399 throw AipsError("Illegal unit - must be compatible with Jy or K");
400 }
401}
402
403void Scantable::setInstrument(const std::string& name)
404{
405 bool throwIt = true;
406 // create an Instrument to see if this is valid
407 STAttr::convertInstrument(name, throwIt);
408 String nameU(name);
409 nameU.upcase();
410 table_.rwKeywordSet().define(String("AntennaName"), nameU);
411}
412
413void Scantable::setFeedType(const std::string& feedtype)
414{
415 if ( Scantable::factories_.find(feedtype) == Scantable::factories_.end() ) {
416 std::string msg = "Illegal feed type "+ feedtype;
417 throw(casa::AipsError(msg));
418 }
419 table_.rwKeywordSet().define(String("POLTYPE"), feedtype);
420}
421
422MPosition Scantable::getAntennaPosition () const
423{
424 Vector<Double> antpos;
425 table_.keywordSet().get("AntennaPosition", antpos);
426 MVPosition mvpos(antpos(0),antpos(1),antpos(2));
427 return MPosition(mvpos);
428}
429
430void Scantable::makePersistent(const std::string& filename)
431{
432 String inname(filename);
433 Path path(inname);
434 /// @todo reindex SCANNO, recompute nbeam, nif, npol
435 inname = path.expandedName();
436 table_.deepCopy(inname, Table::New);
437}
438
439int Scantable::nbeam( int scanno ) const
440{
441 if ( scanno < 0 ) {
442 Int n;
443 table_.keywordSet().get("nBeam",n);
444 return int(n);
445 } else {
446 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
447 Table t = table_(table_.col("SCANNO") == scanno);
448 ROTableRow row(t);
449 const TableRecord& rec = row.get(0);
450 Table subt = t( t.col("IFNO") == Int(rec.asuInt("IFNO"))
451 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
452 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
453 ROTableVector<uInt> v(subt, "BEAMNO");
454 return int(v.nelements());
455 }
456 return 0;
457}
458
459int Scantable::nif( int scanno ) const
460{
461 if ( scanno < 0 ) {
462 Int n;
463 table_.keywordSet().get("nIF",n);
464 return int(n);
465 } else {
466 // take the first POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
467 Table t = table_(table_.col("SCANNO") == scanno);
468 ROTableRow row(t);
469 const TableRecord& rec = row.get(0);
470 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
471 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
472 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
473 if ( subt.nrow() == 0 ) return 0;
474 ROTableVector<uInt> v(subt, "IFNO");
475 return int(v.nelements());
476 }
477 return 0;
478}
479
480int Scantable::npol( int scanno ) const
481{
482 if ( scanno < 0 ) {
483 Int n;
484 table_.keywordSet().get("nPol",n);
485 return n;
486 } else {
487 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
488 Table t = table_(table_.col("SCANNO") == scanno);
489 ROTableRow row(t);
490 const TableRecord& rec = row.get(0);
491 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
492 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
493 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
494 if ( subt.nrow() == 0 ) return 0;
495 ROTableVector<uInt> v(subt, "POLNO");
496 return int(v.nelements());
497 }
498 return 0;
499}
500
501int Scantable::ncycle( int scanno ) const
502{
503 if ( scanno < 0 ) {
504 Block<String> cols(2);
505 cols[0] = "SCANNO";
506 cols[1] = "CYCLENO";
507 TableIterator it(table_, cols);
508 int n = 0;
509 while ( !it.pastEnd() ) {
510 ++n;
511 ++it;
512 }
513 return n;
514 } else {
515 Table t = table_(table_.col("SCANNO") == scanno);
516 ROTableRow row(t);
517 const TableRecord& rec = row.get(0);
518 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
519 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
520 && t.col("IFNO") == Int(rec.asuInt("IFNO")) );
521 if ( subt.nrow() == 0 ) return 0;
522 return int(subt.nrow());
523 }
524 return 0;
525}
526
527
528int Scantable::nrow( int scanno ) const
529{
530 return int(table_.nrow());
531}
532
533int Scantable::nchan( int ifno ) const
534{
535 if ( ifno < 0 ) {
536 Int n;
537 table_.keywordSet().get("nChan",n);
538 return int(n);
539 } else {
540 // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't
541 // vary with these
542 Table t = table_(table_.col("IFNO") == ifno);
543 if ( t.nrow() == 0 ) return 0;
544 ROArrayColumn<Float> v(t, "SPECTRA");
545 return v.shape(0)(0);
546 }
547 return 0;
548}
549
550int Scantable::nscan() const {
551 Vector<uInt> scannos(scanCol_.getColumn());
552 uInt nout = genSort( scannos, Sort::Ascending,
553 Sort::QuickSort|Sort::NoDuplicates );
554 return int(nout);
555}
556
557int Scantable::getChannels(int whichrow) const
558{
559 return specCol_.shape(whichrow)(0);
560}
561
562int Scantable::getBeam(int whichrow) const
563{
564 return beamCol_(whichrow);
565}
566
567std::vector<uint> Scantable::getNumbers(ScalarColumn<uInt>& col)
568{
569 Vector<uInt> nos(col.getColumn());
570 uInt n = genSort( nos, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates );
571 nos.resize(n, True);
572 std::vector<uint> stlout;
573 nos.tovector(stlout);
574 return stlout;
575}
576
577int Scantable::getIF(int whichrow) const
578{
579 return ifCol_(whichrow);
580}
581
582int Scantable::getPol(int whichrow) const
583{
584 return polCol_(whichrow);
585}
586
587std::string Scantable::formatTime(const MEpoch& me, bool showdate) const
588{
589 MVTime mvt(me.getValue());
590 if (showdate)
591 mvt.setFormat(MVTime::YMD);
592 else
593 mvt.setFormat(MVTime::TIME);
594 ostringstream oss;
595 oss << mvt;
596 return String(oss);
597}
598
599void Scantable::calculateAZEL()
600{
601 MPosition mp = getAntennaPosition();
602 MEpoch::ROScalarColumn timeCol(table_, "TIME");
603 ostringstream oss;
604 oss << "Computed azimuth/elevation using " << endl
605 << mp << endl;
606 for (Int i=0; i<nrow(); ++i) {
607 MEpoch me = timeCol(i);
608 MDirection md = getDirection(i);
609 oss << " Time: " << formatTime(me,False) << " Direction: " << formatDirection(md)
610 << endl << " => ";
611 MeasFrame frame(mp, me);
612 Vector<Double> azel =
613 MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
614 frame)
615 )().getAngle("rad").getValue();
616 azCol_.put(i,Float(azel[0]));
617 elCol_.put(i,Float(azel[1]));
618 oss << "azel: " << azel[0]/C::pi*180.0 << " "
619 << azel[1]/C::pi*180.0 << " (deg)" << endl;
620 }
621 pushLog(String(oss));
622}
623
624void Scantable::flag(const std::vector<bool>& msk, bool unflag)
625{
626 std::vector<bool>::const_iterator it;
627 uInt ntrue = 0;
628 for (it = msk.begin(); it != msk.end(); ++it) {
629 if ( *it ) {
630 ntrue++;
631 }
632 }
633 if ( selector_.empty() && (msk.size() == 0 || msk.size() == ntrue) )
634 throw(AipsError("Trying to flag whole scantable."));
635 if ( msk.size() == 0 ) {
636 uChar userflag = 1 << 7;
637 if ( unflag ) {
638 userflag = 0 << 7;
639 }
640 for ( uInt i=0; i<table_.nrow(); ++i) {
641 Vector<uChar> flgs = flagsCol_(i);
642 flgs = userflag;
643 flagsCol_.put(i, flgs);
644 }
645 return;
646 }
647 if ( int(msk.size()) != nchan() ) {
648 throw(AipsError("Mask has incorrect number of channels."));
649 }
650 for ( uInt i=0; i<table_.nrow(); ++i) {
651 Vector<uChar> flgs = flagsCol_(i);
652 if ( flgs.nelements() != msk.size() ) {
653 throw(AipsError("Mask has incorrect number of channels."
654 " Probably varying with IF. Please flag per IF"));
655 }
656 std::vector<bool>::const_iterator it;
657 uInt j = 0;
658 uChar userflag = 1 << 7;
659 if ( unflag ) {
660 userflag = 0 << 7;
661 }
662 for (it = msk.begin(); it != msk.end(); ++it) {
663 if ( *it ) {
664 flgs(j) = userflag;
665 }
666 ++j;
667 }
668 flagsCol_.put(i, flgs);
669 }
670}
671
672std::vector<bool> Scantable::getMask(int whichrow) const
673{
674 Vector<uChar> flags;
675 flagsCol_.get(uInt(whichrow), flags);
676 Vector<Bool> bflag(flags.shape());
677 convertArray(bflag, flags);
678 bflag = !bflag;
679 std::vector<bool> mask;
680 bflag.tovector(mask);
681 return mask;
682}
683
684std::vector<float> Scantable::getSpectrum( int whichrow,
685 const std::string& poltype ) const
686{
687 String ptype = poltype;
688 if (poltype == "" ) ptype = getPolType();
689 if ( whichrow < 0 || whichrow >= nrow() )
690 throw(AipsError("Illegal row number."));
691 std::vector<float> out;
692 Vector<Float> arr;
693 uInt requestedpol = polCol_(whichrow);
694 String basetype = getPolType();
695 if ( ptype == basetype ) {
696 specCol_.get(whichrow, arr);
697 } else {
698 CountedPtr<STPol> stpol(STPol::getPolClass(Scantable::factories_, basetype));
699 uInt row = uInt(whichrow);
700 stpol->setSpectra(getPolMatrix(row));
701 Float fang,fhand,parang;
702 fang = focusTable_.getTotalFeedAngle(mfocusidCol_(row));
703 fhand = focusTable_.getFeedHand(mfocusidCol_(row));
704 parang = paraCol_(row);
705 /// @todo re-enable this
706 // disable total feed angle to support paralactifying Caswell style
707 stpol->setPhaseCorrections(parang, -parang, fhand);
708 arr = stpol->getSpectrum(requestedpol, ptype);
709 }
710 if ( arr.nelements() == 0 )
711 pushLog("Not enough polarisations present to do the conversion.");
712 arr.tovector(out);
713 return out;
714}
715
716void Scantable::setSpectrum( const std::vector<float>& spec,
717 int whichrow )
718{
719 Vector<Float> spectrum(spec);
720 Vector<Float> arr;
721 specCol_.get(whichrow, arr);
722 if ( spectrum.nelements() != arr.nelements() )
723 throw AipsError("The spectrum has incorrect number of channels.");
724 specCol_.put(whichrow, spectrum);
725}
726
727
728String Scantable::generateName()
729{
730 return (File::newUniqueName("./","temp")).baseName();
731}
732
733const casa::Table& Scantable::table( ) const
734{
735 return table_;
736}
737
738casa::Table& Scantable::table( )
739{
740 return table_;
741}
742
743std::string Scantable::getPolType() const
744{
745 return table_.keywordSet().asString("POLTYPE");
746}
747
748void Scantable::unsetSelection()
749{
750 table_ = originalTable_;
751 attach();
752 selector_.reset();
753}
754
755void Scantable::setSelection( const STSelector& selection )
756{
757 Table tab = const_cast<STSelector&>(selection).apply(originalTable_);
758 if ( tab.nrow() == 0 ) {
759 throw(AipsError("Selection contains no data. Not applying it."));
760 }
761 table_ = tab;
762 attach();
763 selector_ = selection;
764}
765
766std::string Scantable::summary( bool verbose )
767{
768 // Format header info
769 ostringstream oss;
770 oss << endl;
771 oss << asap::SEPERATOR << endl;
772 oss << " Scan Table Summary" << endl;
773 oss << asap::SEPERATOR << endl;
774 oss.flags(std::ios_base::left);
775 oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
776 << setw(15) << "IFs:" << setw(4) << nif() << endl
777 << setw(15) << "Polarisations:" << setw(4) << npol()
778 << "(" << getPolType() << ")" << endl
779 << setw(15) << "Channels:" << setw(4) << nchan() << endl;
780 oss << endl;
781 String tmp;
782 oss << setw(15) << "Observer:"
783 << table_.keywordSet().asString("Observer") << endl;
784 oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
785 table_.keywordSet().get("Project", tmp);
786 oss << setw(15) << "Project:" << tmp << endl;
787 table_.keywordSet().get("Obstype", tmp);
788 oss << setw(15) << "Obs. Type:" << tmp << endl;
789 table_.keywordSet().get("AntennaName", tmp);
790 oss << setw(15) << "Antenna Name:" << tmp << endl;
791 table_.keywordSet().get("FluxUnit", tmp);
792 oss << setw(15) << "Flux Unit:" << tmp << endl;
793 //Vector<Double> vec(moleculeTable_.getRestFrequencies());
794 int nid = moleculeTable_.nrow();
795 Bool firstline = True;
796 oss << setw(15) << "Rest Freqs:";
797 for (int i=0; i<nid; i++) {
798 Table t = table_(table_.col("MOLECULE_ID") == i);
799 if (t.nrow() > 0) {
800 Vector<Double> vec(moleculeTable_.getRestFrequency(i));
801 if (vec.nelements() > 0) {
802 if (firstline) {
803 oss << setprecision(10) << vec << " [Hz]" << endl;
804 firstline=False;
805 }
806 else{
807 oss << setw(15)<<" " << setprecision(10) << vec << " [Hz]" << endl;
808 }
809 } else {
810 oss << "none" << endl;
811 }
812 }
813 }
814
815 oss << setw(15) << "Abcissa:" << getAbcissaLabel(0) << endl;
816 oss << selector_.print() << endl;
817 oss << endl;
818 // main table
819 String dirtype = "Position ("
820 + getDirectionRefString()
821 + ")";
822 oss << setw(5) << "Scan" << setw(15) << "Source"
823 << setw(10) << "Time" << setw(18) << "Integration" << endl;
824 oss << setw(5) << "" << setw(5) << "Beam" << setw(3) << "" << dirtype << endl;
825 oss << setw(10) << "" << setw(3) << "IF" << setw(6) << ""
826 << setw(8) << "Frame" << setw(16)
827 << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment" <<endl;
828 oss << asap::SEPERATOR << endl;
829 TableIterator iter(table_, "SCANNO");
830 while (!iter.pastEnd()) {
831 Table subt = iter.table();
832 ROTableRow row(subt);
833 MEpoch::ROScalarColumn timeCol(subt,"TIME");
834 const TableRecord& rec = row.get(0);
835 oss << setw(4) << std::right << rec.asuInt("SCANNO")
836 << std::left << setw(1) << ""
837 << setw(15) << rec.asString("SRCNAME")
838 << setw(10) << formatTime(timeCol(0), false);
839 // count the cycles in the scan
840 TableIterator cyciter(subt, "CYCLENO");
841 int nint = 0;
842 while (!cyciter.pastEnd()) {
843 ++nint;
844 ++cyciter;
845 }
846 oss << setw(3) << std::right << nint << setw(3) << " x " << std::left
847 << setw(6) << formatSec(rec.asFloat("INTERVAL")) << endl;
848
849 TableIterator biter(subt, "BEAMNO");
850 while (!biter.pastEnd()) {
851 Table bsubt = biter.table();
852 ROTableRow brow(bsubt);
853 const TableRecord& brec = brow.get(0);
854 uInt row0 = bsubt.rowNumbers(table_)[0];
855 oss << setw(5) << "" << setw(4) << std::right << brec.asuInt("BEAMNO")<< std::left;
856 oss << setw(4) << "" << formatDirection(getDirection(row0)) << endl;
857 TableIterator iiter(bsubt, "IFNO");
858 while (!iiter.pastEnd()) {
859 Table isubt = iiter.table();
860 ROTableRow irow(isubt);
861 const TableRecord& irec = irow.get(0);
862 oss << setw(10) << "";
863 oss << setw(3) << std::right << irec.asuInt("IFNO") << std::left
864 << setw(2) << "" << frequencies().print(irec.asuInt("FREQ_ID"))
865 << endl;
866
867 ++iiter;
868 }
869 ++biter;
870 }
871 ++iter;
872 }
873 /// @todo implement verbose mode
874 return String(oss);
875}
876
877std::string Scantable::getTime(int whichrow, bool showdate) const
878{
879 MEpoch::ROScalarColumn timeCol(table_, "TIME");
880 MEpoch me;
881 if (whichrow > -1) {
882 me = timeCol(uInt(whichrow));
883 } else {
884 Double tm;
885 table_.keywordSet().get("UTC",tm);
886 me = MEpoch(MVEpoch(tm));
887 }
888 return formatTime(me, showdate);
889}
890
891MEpoch Scantable::getEpoch(int whichrow) const
892{
893 if (whichrow > -1) {
894 return timeCol_(uInt(whichrow));
895 } else {
896 Double tm;
897 table_.keywordSet().get("UTC",tm);
898 return MEpoch(MVEpoch(tm));
899 }
900}
901
902std::string Scantable::getDirectionString(int whichrow) const
903{
904 return formatDirection(getDirection(uInt(whichrow)));
905}
906
907std::vector< double > Scantable::getAbcissa( int whichrow ) const
908{
909 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal row number"));
910 std::vector<double> stlout;
911 int nchan = specCol_(whichrow).nelements();
912 String us = freqTable_.getUnitString();
913 if ( us == "" || us == "pixel" || us == "channel" ) {
914 for (int i=0; i<nchan; ++i) {
915 stlout.push_back(double(i));
916 }
917 return stlout;
918 }
919
920 const MPosition& mp = getAntennaPosition();
921 const MDirection& md = getDirection(whichrow);
922 const MEpoch& me = timeCol_(whichrow);
923 //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
924 Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
925 SpectralCoordinate spc =
926 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
927 Vector<Double> pixel(nchan);
928 Vector<Double> world;
929 indgen(pixel);
930 if ( Unit(us) == Unit("Hz") ) {
931 for ( int i=0; i < nchan; ++i) {
932 Double world;
933 spc.toWorld(world, pixel[i]);
934 stlout.push_back(double(world));
935 }
936 } else if ( Unit(us) == Unit("km/s") ) {
937 Vector<Double> world;
938 spc.pixelToVelocity(world, pixel);
939 world.tovector(stlout);
940 }
941 return stlout;
942}
943void Scantable::setDirectionRefString( const std::string & refstr )
944{
945 MDirection::Types mdt;
946 if (refstr != "" && !MDirection::getType(mdt, refstr)) {
947 throw(AipsError("Illegal Direction frame."));
948 }
949 if ( refstr == "" ) {
950 String defaultstr = MDirection::showType(dirCol_.getMeasRef().getType());
951 table_.rwKeywordSet().define("DIRECTIONREF", defaultstr);
952 } else {
953 table_.rwKeywordSet().define("DIRECTIONREF", String(refstr));
954 }
955}
956
957std::string Scantable::getDirectionRefString( ) const
958{
959 return table_.keywordSet().asString("DIRECTIONREF");
960}
961
962MDirection Scantable::getDirection(int whichrow ) const
963{
964 String usertype = table_.keywordSet().asString("DIRECTIONREF");
965 String type = MDirection::showType(dirCol_.getMeasRef().getType());
966 if ( usertype != type ) {
967 MDirection::Types mdt;
968 if (!MDirection::getType(mdt, usertype)) {
969 throw(AipsError("Illegal Direction frame."));
970 }
971 return dirCol_.convert(uInt(whichrow), mdt);
972 } else {
973 return dirCol_(uInt(whichrow));
974 }
975}
976
977std::string Scantable::getAbcissaLabel( int whichrow ) const
978{
979 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal ro number"));
980 const MPosition& mp = getAntennaPosition();
981 const MDirection& md = getDirection(whichrow);
982 const MEpoch& me = timeCol_(whichrow);
983 //const Double& rf = mmolidCol_(whichrow);
984 const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
985 SpectralCoordinate spc =
986 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
987
988 String s = "Channel";
989 Unit u = Unit(freqTable_.getUnitString());
990 if (u == Unit("km/s")) {
991 s = CoordinateUtil::axisLabel(spc, 0, True,True, True);
992 } else if (u == Unit("Hz")) {
993 Vector<String> wau(1);wau = u.getName();
994 spc.setWorldAxisUnits(wau);
995 s = CoordinateUtil::axisLabel(spc, 0, True, True, False);
996 }
997 return s;
998
999}
1000
1001/**
1002void asap::Scantable::setRestFrequencies( double rf, const std::string& name,
1003 const std::string& unit )
1004**/
1005void Scantable::setRestFrequencies( vector<double> rf, const vector<std::string>& name,
1006 const std::string& unit )
1007
1008{
1009 ///@todo lookup in line table to fill in name and formattedname
1010 Unit u(unit);
1011 //Quantum<Double> urf(rf, u);
1012 Quantum<Vector<Double> >urf(rf, u);
1013 Vector<String> formattedname(0);
1014 //cerr<<"Scantable::setRestFrequnecies="<<urf<<endl;
1015
1016 //uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
1017 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), mathutil::toVectorString(name), formattedname);
1018 TableVector<uInt> tabvec(table_, "MOLECULE_ID");
1019 tabvec = id;
1020}
1021
1022/**
1023void asap::Scantable::setRestFrequencies( const std::string& name )
1024{
1025 throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
1026 ///@todo implement
1027}
1028**/
1029void Scantable::setRestFrequencies( const vector<std::string>& name )
1030{
1031 throw(AipsError("setRestFrequencies( const vector<std::string>& name ) NYI"));
1032 ///@todo implement
1033}
1034
1035std::vector< unsigned int > Scantable::rownumbers( ) const
1036{
1037 std::vector<unsigned int> stlout;
1038 Vector<uInt> vec = table_.rowNumbers();
1039 vec.tovector(stlout);
1040 return stlout;
1041}
1042
1043
1044Matrix<Float> Scantable::getPolMatrix( uInt whichrow ) const
1045{
1046 ROTableRow row(table_);
1047 const TableRecord& rec = row.get(whichrow);
1048 Table t =
1049 originalTable_( originalTable_.col("SCANNO") == Int(rec.asuInt("SCANNO"))
1050 && originalTable_.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
1051 && originalTable_.col("IFNO") == Int(rec.asuInt("IFNO"))
1052 && originalTable_.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
1053 ROArrayColumn<Float> speccol(t, "SPECTRA");
1054 return speccol.getColumn();
1055}
1056
1057std::vector< std::string > Scantable::columnNames( ) const
1058{
1059 Vector<String> vec = table_.tableDesc().columnNames();
1060 return mathutil::tovectorstring(vec);
1061}
1062
1063MEpoch::Types Scantable::getTimeReference( ) const
1064{
1065 return MEpoch::castType(timeCol_.getMeasRef().getType());
1066}
1067
1068void Scantable::addFit( const STFitEntry& fit, int row )
1069{
1070 //cout << mfitidCol_(uInt(row)) << endl;
1071 LogIO os( LogOrigin( "Scantable", "addFit()", WHERE ) ) ;
1072 os << mfitidCol_(uInt(row)) << LogIO::POST ;
1073 uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row)));
1074 mfitidCol_.put(uInt(row), id);
1075}
1076
1077void Scantable::shift(int npix)
1078{
1079 Vector<uInt> fids(mfreqidCol_.getColumn());
1080 genSort( fids, Sort::Ascending,
1081 Sort::QuickSort|Sort::NoDuplicates );
1082 for (uInt i=0; i<fids.nelements(); ++i) {
1083 frequencies().shiftRefPix(npix, i);
1084 }
1085}
1086
1087std::string asap::Scantable::getAntennaName() const
1088{
1089 String out;
1090 table_.keywordSet().get("AntennaName", out);
1091 return out;
1092}
1093
1094int asap::Scantable::checkScanInfo(const std::vector<int>& scanlist) const
1095{
1096 String tbpath;
1097 int ret = 0;
1098 if ( table_.keywordSet().isDefined("GBT_GO") ) {
1099 table_.keywordSet().get("GBT_GO", tbpath);
1100 Table t(tbpath,Table::Old);
1101 // check each scan if other scan of the pair exist
1102 int nscan = scanlist.size();
1103 for (int i = 0; i < nscan; i++) {
1104 Table subt = t( t.col("SCAN") == scanlist[i]+1 );
1105 if (subt.nrow()==0) {
1106 //cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
1107 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1108 os <<LogIO::WARN<<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<LogIO::POST;
1109 ret = 1;
1110 break;
1111 }
1112 ROTableRow row(subt);
1113 const TableRecord& rec = row.get(0);
1114 int scan1seqn = rec.asuInt("PROCSEQN");
1115 int laston1 = rec.asuInt("LASTON");
1116 if ( rec.asuInt("PROCSIZE")==2 ) {
1117 if ( i < nscan-1 ) {
1118 Table subt2 = t( t.col("SCAN") == scanlist[i+1]+1 );
1119 if ( subt2.nrow() == 0) {
1120 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1121
1122 //cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
1123 os<<LogIO::WARN<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<LogIO::POST;
1124 ret = 1;
1125 break;
1126 }
1127 ROTableRow row2(subt2);
1128 const TableRecord& rec2 = row2.get(0);
1129 int scan2seqn = rec2.asuInt("PROCSEQN");
1130 int laston2 = rec2.asuInt("LASTON");
1131 if (scan1seqn == 1 && scan2seqn == 2) {
1132 if (laston1 == laston2) {
1133 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1134 //cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1135 os<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
1136 i +=1;
1137 }
1138 else {
1139 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1140 //cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1141 os<<LogIO::WARN<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
1142 }
1143 }
1144 else if (scan1seqn==2 && scan2seqn == 1) {
1145 if (laston1 == laston2) {
1146 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1147 //cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
1148 os<<LogIO::WARN<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<LogIO::POST;
1149 ret = 1;
1150 break;
1151 }
1152 }
1153 else {
1154 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1155 //cerr<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
1156 os<<LogIO::WARN<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<LogIO::POST;
1157 ret = 1;
1158 break;
1159 }
1160 }
1161 }
1162 else {
1163 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1164 //cerr<<"The scan does not appear to be standard obsevation."<<endl;
1165 os<<LogIO::WARN<<"The scan does not appear to be standard obsevation."<<LogIO::POST;
1166 }
1167 //if ( i >= nscan ) break;
1168 }
1169 }
1170 else {
1171 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1172 //cerr<<"No reference to GBT_GO table."<<endl;
1173 os<<LogIO::WARN<<"No reference to GBT_GO table."<<LogIO::POST;
1174 ret = 1;
1175 }
1176 return ret;
1177}
1178
1179std::vector<double> asap::Scantable::getDirectionVector(int whichrow) const
1180{
1181 Vector<Double> Dir = dirCol_(whichrow).getAngle("rad").getValue();
1182 std::vector<double> dir;
1183 Dir.tovector(dir);
1184 return dir;
1185}
1186
1187void asap::Scantable::reshapeSpectrum( int nmin, int nmax )
1188 throw( casa::AipsError )
1189{
1190 // assumed that all rows have same nChan
1191 Vector<Float> arr = specCol_( 0 ) ;
1192 int nChan = arr.nelements() ;
1193
1194 // if nmin < 0 or nmax < 0, nothing to do
1195 if ( nmin < 0 ) {
1196 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1197 }
1198 if ( nmax < 0 ) {
1199 throw( casa::indexError<int>( nmax, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1200 }
1201
1202 // if nmin > nmax, exchange values
1203 if ( nmin > nmax ) {
1204 int tmp = nmax ;
1205 nmax = nmin ;
1206 nmin = tmp ;
1207 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1208 os << "Swap values. Applied range is ["
1209 << nmin << ", " << nmax << "]" << LogIO::POST ;
1210 }
1211
1212 // if nmin exceeds nChan, nothing to do
1213 if ( nmin >= nChan ) {
1214 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Specified minimum exceeds nChan." ) ) ;
1215 }
1216
1217 // if nmax exceeds nChan, reset nmax to nChan
1218 if ( nmax >= nChan ) {
1219 if ( nmin == 0 ) {
1220 // nothing to do
1221 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1222 os << "Whole range is selected. Nothing to do." << LogIO::POST ;
1223 return ;
1224 }
1225 else {
1226 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1227 os << "Specified maximum exceeds nChan. Applied range is ["
1228 << nmin << ", " << nChan-1 << "]." << LogIO::POST ;
1229 nmax = nChan - 1 ;
1230 }
1231 }
1232
1233 // reshape specCol_ and flagCol_
1234 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1235 reshapeSpectrum( nmin, nmax, irow ) ;
1236 }
1237
1238 // update FREQUENCIES subtable
1239 Double refpix ;
1240 Double refval ;
1241 Double increment ;
1242 int freqnrow = freqTable_.table().nrow() ;
1243 Vector<uInt> oldId( freqnrow ) ;
1244 Vector<uInt> newId( freqnrow ) ;
1245 for ( int irow = 0 ; irow < freqnrow ; irow++ ) {
1246 freqTable_.getEntry( refpix, refval, increment, irow ) ;
1247 /***
1248 * need to shift refpix to nmin
1249 * note that channel nmin in old index will be channel 0 in new one
1250 ***/
1251 refval = refval - ( refpix - nmin ) * increment ;
1252 refpix = 0 ;
1253 freqTable_.setEntry( refpix, refval, increment, irow ) ;
1254 }
1255
1256 // update nchan
1257 int newsize = nmax - nmin + 1 ;
1258 table_.rwKeywordSet().define( "nChan", newsize ) ;
1259
1260 // update bandwidth
1261 // assumed all spectra in the scantable have same bandwidth
1262 table_.rwKeywordSet().define( "Bandwidth", increment * newsize ) ;
1263
1264 return ;
1265}
1266
1267void asap::Scantable::reshapeSpectrum( int nmin, int nmax, int irow )
1268{
1269 // reshape specCol_ and flagCol_
1270 Vector<Float> oldspec = specCol_( irow ) ;
1271 Vector<uChar> oldflag = flagsCol_( irow ) ;
1272 uInt newsize = nmax - nmin + 1 ;
1273 specCol_.put( irow, oldspec( Slice( nmin, newsize, 1 ) ) ) ;
1274 flagsCol_.put( irow, oldflag( Slice( nmin, newsize, 1 ) ) ) ;
1275
1276 return ;
1277}
1278
1279void asap::Scantable::regridChannel( int nChan, double dnu )
1280{
1281 LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
1282 os << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << LogIO::POST ;
1283 // assumed that all rows have same nChan
1284 Vector<Float> arr = specCol_( 0 ) ;
1285 int oldsize = arr.nelements() ;
1286
1287 // if oldsize == nChan, nothing to do
1288 if ( oldsize == nChan ) {
1289 os << "Specified channel number is same as current one. Nothing to do." << LogIO::POST ;
1290 return ;
1291 }
1292
1293 // if oldChan < nChan, unphysical operation
1294 if ( oldsize < nChan ) {
1295 os << "Unphysical operation. Nothing to do." << LogIO::POST ;
1296 return ;
1297 }
1298
1299 // change channel number for specCol_ and flagCol_
1300 Vector<Float> newspec( nChan, 0 ) ;
1301 Vector<uChar> newflag( nChan, false ) ;
1302 vector<string> coordinfo = getCoordInfo() ;
1303 string oldinfo = coordinfo[0] ;
1304 coordinfo[0] = "Hz" ;
1305 setCoordInfo( coordinfo ) ;
1306 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1307 regridChannel( nChan, dnu, irow ) ;
1308 }
1309 coordinfo[0] = oldinfo ;
1310 setCoordInfo( coordinfo ) ;
1311
1312
1313 // NOTE: this method does not update metadata such as
1314 // FREQUENCIES subtable, nChan, Bandwidth, etc.
1315
1316 return ;
1317}
1318
1319void asap::Scantable::regridChannel( int nChan, double dnu, int irow )
1320{
1321 // logging
1322 //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ;
1323 //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ;
1324
1325 Vector<Float> oldspec = specCol_( irow ) ;
1326 Vector<uChar> oldflag = flagsCol_( irow ) ;
1327 Vector<Float> newspec( nChan, 0 ) ;
1328 Vector<uChar> newflag( nChan, false ) ;
1329
1330 // regrid
1331 vector<double> abcissa = getAbcissa( irow ) ;
1332 int oldsize = abcissa.size() ;
1333 double olddnu = abcissa[1] - abcissa[0] ;
1334 //int refChan = 0 ;
1335 //double frac = 0.0 ;
1336 //double wedge = 0.0 ;
1337 //double pile = 0.0 ;
1338 int ichan = 0 ;
1339 double wsum = 0.0 ;
1340 Vector<Float> z( nChan ) ;
1341 z[0] = abcissa[0] - 0.5 * olddnu + 0.5 * dnu ;
1342 for ( int ii = 1 ; ii < nChan ; ii++ )
1343 z[ii] = z[ii-1] + dnu ;
1344 Vector<Float> zi( nChan+1 ) ;
1345 Vector<Float> yi( oldsize + 1 ) ;
1346 zi[0] = z[0] - 0.5 * dnu ;
1347 zi[1] = z[0] + 0.5 * dnu ;
1348 for ( int ii = 2 ; ii < nChan ; ii++ )
1349 zi[ii] = zi[ii-1] + dnu ;
1350 zi[nChan] = z[nChan-1] + 0.5 * dnu ;
1351 yi[0] = abcissa[0] - 0.5 * olddnu ;
1352 yi[1] = abcissa[1] + 0.5 * olddnu ;
1353 for ( int ii = 2 ; ii < oldsize ; ii++ )
1354 yi[ii] = abcissa[ii-1] + olddnu ;
1355 yi[oldsize] = abcissa[oldsize-1] + 0.5 * olddnu ;
1356 if ( dnu > 0.0 ) {
1357 for ( int ii = 0 ; ii < nChan ; ii++ ) {
1358 double zl = zi[ii] ;
1359 double zr = zi[ii+1] ;
1360 for ( int j = ichan ; j < oldsize ; j++ ) {
1361 double yl = yi[j] ;
1362 double yr = yi[j+1] ;
1363 if ( yl <= zl ) {
1364 if ( yr <= zl ) {
1365 continue ;
1366 }
1367 else if ( yr <= zr ) {
1368 newspec[ii] += oldspec[j] * ( yr - zl ) ;
1369 newflag[ii] = newflag[ii] || oldflag[j] ;
1370 wsum += ( yr - zl ) ;
1371 }
1372 else {
1373 newspec[ii] += oldspec[j] * dnu ;
1374 newflag[ii] = newflag[ii] || oldflag[j] ;
1375 wsum += dnu ;
1376 ichan = j ;
1377 break ;
1378 }
1379 }
1380 else if ( yl < zr ) {
1381 if ( yr <= zr ) {
1382 newspec[ii] += oldspec[j] * ( yr - yl ) ;
1383 newflag[ii] = newflag[ii] || oldflag[j] ;
1384 wsum += ( yr - yl ) ;
1385 }
1386 else {
1387 newspec[ii] += oldspec[j] * ( zr - yl ) ;
1388 newflag[ii] = newflag[ii] || oldflag[j] ;
1389 wsum += ( zr - yl ) ;
1390 ichan = j ;
1391 break ;
1392 }
1393 }
1394 else {
1395 ichan = j - 1 ;
1396 break ;
1397 }
1398 }
1399 newspec[ii] /= wsum ;
1400 wsum = 0.0 ;
1401 }
1402 }
1403 else if ( dnu < 0.0 ) {
1404 for ( int ii = 0 ; ii < nChan ; ii++ ) {
1405 double zl = zi[ii] ;
1406 double zr = zi[ii+1] ;
1407 for ( int j = ichan ; j < oldsize ; j++ ) {
1408 double yl = yi[j] ;
1409 double yr = yi[j+1] ;
1410 if ( yl >= zl ) {
1411 if ( yr >= zl ) {
1412 continue ;
1413 }
1414 else if ( yr >= zr ) {
1415 newspec[ii] += oldspec[j] * abs( yr - zl ) ;
1416 newflag[ii] = newflag[ii] || oldflag[j] ;
1417 wsum += abs( yr - zl ) ;
1418 }
1419 else {
1420 newspec[ii] += oldspec[j] * abs( dnu ) ;
1421 newflag[ii] = newflag[ii] || oldflag[j] ;
1422 wsum += abs( dnu ) ;
1423 ichan = j ;
1424 break ;
1425 }
1426 }
1427 else if ( yl > zr ) {
1428 if ( yr >= zr ) {
1429 newspec[ii] += oldspec[j] * abs( yr - yl ) ;
1430 newflag[ii] = newflag[ii] || oldflag[j] ;
1431 wsum += abs( yr - yl ) ;
1432 }
1433 else {
1434 newspec[ii] += oldspec[j] * abs( zr - yl ) ;
1435 newflag[ii] = newflag[ii] || oldflag[j] ;
1436 wsum += abs( zr - yl ) ;
1437 ichan = j ;
1438 break ;
1439 }
1440 }
1441 else {
1442 ichan = j - 1 ;
1443 break ;
1444 }
1445 }
1446 newspec[ii] /= wsum ;
1447 wsum = 0.0 ;
1448 }
1449 }
1450// * ichan = 0
1451// ***/
1452// //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
1453// pile += dnu ;
1454// wedge = olddnu * ( refChan + 1 ) ;
1455// while ( wedge < pile ) {
1456// newspec[0] += olddnu * oldspec[refChan] ;
1457// newflag[0] = newflag[0] || oldflag[refChan] ;
1458// //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
1459// refChan++ ;
1460// wedge += olddnu ;
1461// wsum += olddnu ;
1462// //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
1463// }
1464// frac = ( wedge - pile ) / olddnu ;
1465// wsum += ( 1.0 - frac ) * olddnu ;
1466// newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
1467// newflag[0] = newflag[0] || oldflag[refChan] ;
1468// //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
1469// //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
1470// newspec[0] /= wsum ;
1471// //ofs << "newspec[0] = " << newspec[0] << endl ;
1472// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
1473
1474// /***
1475// * ichan = 1 - nChan-2
1476// ***/
1477// for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
1478// pile += dnu ;
1479// newspec[ichan] += frac * olddnu * oldspec[refChan] ;
1480// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
1481// //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
1482// refChan++ ;
1483// wedge += olddnu ;
1484// wsum = frac * olddnu ;
1485// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
1486// while ( wedge < pile ) {
1487// newspec[ichan] += olddnu * oldspec[refChan] ;
1488// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
1489// //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
1490// refChan++ ;
1491// wedge += olddnu ;
1492// wsum += olddnu ;
1493// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
1494// }
1495// frac = ( wedge - pile ) / olddnu ;
1496// wsum += ( 1.0 - frac ) * olddnu ;
1497// newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
1498// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
1499// //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
1500// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
1501// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
1502// newspec[ichan] /= wsum ;
1503// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
1504// }
1505
1506// /***
1507// * ichan = nChan-1
1508// ***/
1509// // NOTE: Assumed that all spectra have the same bandwidth
1510// pile += dnu ;
1511// newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
1512// newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
1513// //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
1514// refChan++ ;
1515// wedge += olddnu ;
1516// wsum = frac * olddnu ;
1517// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
1518// for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
1519// newspec[nChan-1] += olddnu * oldspec[jchan] ;
1520// newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
1521// wsum += olddnu ;
1522// //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
1523// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
1524// }
1525// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
1526// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
1527// newspec[nChan-1] /= wsum ;
1528// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
1529
1530// specCol_.put( irow, newspec ) ;
1531// flagsCol_.put( irow, newflag ) ;
1532
1533// // ofs.close() ;
1534
1535
1536 return ;
1537}
1538
1539}
1540//namespace asap
Note: See TracBrowser for help on using the repository browser.