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

Last change on this file since 1538 was 1451, checked in by TakTsutsumi, 16 years ago

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: No

What Interface Changed:

Test Programs: scantable.summary()

Put in Release Notes: No

Description: Added missing endl in summary


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