//# Correspondence concerning AIPS++ should be addressed as follows:
//# Internet email: aips2-request@nrao.edu.
//# Postal address: AIPS++ Project Office
//# National Radio Astronomy Observatory
//# Charlottesville, VA 22903-2475 USA
//----------------------------------------------------------------------------
#include <measures/Measures/Muvw.h>
#include <measures/Measures/MCBaseline.h>
#include <msvis/MSVis/ViImplementation2.h>
#include <synthesis/CalTables/CTIter.h>
#include <tables/Tables/ScalarColumn.h>
#include <casa/OS/Timer.h>
using namespace casacore;
namespace casa { //# NAMESPACE CASA - BEGIN
//----------------------------------------------------------------------------
ROCTIter::ROCTIter(NewCalTable tab, const Block<String>& sortcol) :
sortCols_(sortcol.begin( ),sortcol.end( )),
ti_=new TableIterator(tab,sortcol);
// Attach initial accessors:
// cout << "CTIter sort columns (" << sortCols_.nelements() << "): " << sortCols_ << endl;
// If SPW a sort column, then
singleSpw_=anyEQ(sortCols_,String("SPECTRAL_WINDOW_ID"));
// Initialize MSDerivedValues and MEpoch
msd_ = new MSDerivedValues();
msd_->setAntennas(calCol_.antenna());
epoch_ = calCol_.timeMeas()(0);
// Initialize antenna positions for uvw
cout << "singleSpw_ = " << boolalpha << singleSpw_ << endl;
cout << "inct_->nrow() = " << inct_->nrow() << endl;
cout << "this->nrow() = " << this->nrow() << endl;
cout << "iROCTMainCols_->spwId() = " << iROCTMainCols_->spwId().getColumn() << endl;
cout << "iROCTMainCols_->spwId()(0) = " << iROCTMainCols_->spwId()(0) << endl;
cout << "thisSpw() = " << this->thisSpw() << endl;
//----------------------------------------------------------------------------
if (ti_!=NULL) delete ti_;
if (iROCTMainCols_!=NULL) delete iROCTMainCols_;
if (inct_!=NULL) delete inct_;
if (msd_!=NULL) delete msd_;
// Advance the TableIterator
// attach accessors to new iteration
// Update phase center and uvw if needed
Double time = thisTime();
void ROCTIter::updatePhaseCenter() {
// Set MSDerivedValues phase center when changed
if (field != lastfield_) {
if ((time != lasttime_) || (field != lastfield_)) {
if (field != lastfield_) {
Double time = thisTime();
MDirection phaseCenter = calCol_.field().phaseDirMeas(field, time);
msd_->setFieldCenter(phaseCenter);
void ROCTIter::updatePhaseCenter() {
// Set MSDerivedValues phase center for field
phaseCenter_ = calCol_.field().phaseDirMeas(thisField(), thisTime());
msd_->setFieldCenter(phaseCenter_);
// Advance the TableIterator
Double ROCTIter::thisTime() const { return iROCTMainCols_->time()(0); };
Vector<Double> ROCTIter::time() const { return iROCTMainCols_->time().getColumn(); };
void ROCTIter::time(Vector<Double>& v) const { iROCTMainCols_->time().getColumn(v); };
Int ROCTIter::thisField() const { return iROCTMainCols_->fieldId()(0); };
casacore::Double ROCTIter::hourang(casacore::Double time) const {
return vi::ViImplementation2::hourangCalculate(time, *msd_, epoch_);
casacore::Float ROCTIter::parang0(casacore::Double time) const {
return vi::ViImplementation2::parang0Calculate(time, *msd_, epoch_);
casacore::Matrix<casacore::Double> ROCTIter::uvw() const {
casacore::Vector<casacore::Int> ant1 = antenna1();
casacore::Vector<casacore::Int> ant2 = antenna2();
auto nbaseline = ant1.size();
casacore::Matrix<casacore::Double> uvw(3, nbaseline);
for (uInt i = 0; i < ant1.size(); ++i) {
uvw.column(i) = antennaUVW_[ant2(i)] - antennaUVW_[ant1(i)];
void ROCTIter::initUVW() {
// Calculate relative positions of antennas
nAnt_ = calCol_.antenna().nrow();
auto antPosMeas = calCol_.antenna().positionMeas();
refAntPos_ = antPosMeas(0); // use first antenna for reference
auto refAntPosValue = refAntPos_.getValue();
// Set up baseline and uvw types
MBaseline::getType(baseline_type_, MPosition::showType(refAntPos_.getRef().getType()));
MBaseline::getType(phasedir_type_, MDirection::showType(phaseCenter_.getRef().getType()));
mvbaselines_.resize(nAnt_);
for (Int ant = 0; ant < nAnt_; ++ant) {
// MVBaselines are basically xyz Vectors, not Measures
mvbaselines_[ant] = MVBaseline(refAntPosValue, antPosMeas(ant).getValue());
void ROCTIter::calculateAntennaUVW() {
// Set antennaUVW_ for current iteration when field or time changes
MEpoch epoch = iROCTMainCols_->timeMeas()(0);
MeasFrame measFrame(refAntPos_, epoch, phaseCenter_);
MBaseline::Ref baselineRef(baseline_type_);
baselineMeas.set(mvbaseline, baselineRef);
baselineMeas.getRefPtr()->set(measFrame);
// Conversion engine to phasedir type
MBaseline::Ref uvwRef(phasedir_type_);
MBaseline::Convert baselineConv(baselineMeas, uvwRef);
// WSRT convention: phase opposite to VLA (l increases toward increasing RA)
bool wsrtConvention = calCol_.observation().telescopeName()(thisObs()) == "WSRT";
antennaUVW_.resize(nAnt_);
for (int i = 0; i < nAnt_; ++i) {
baselineMeas.set(mvbaselines_[i], baselineRef);
MBaseline baselineOutFrame = baselineConv(baselineMeas);
MVuvw uvwOutFrame(baselineOutFrame.getValue(), phaseCenter_.getValue());
antennaUVW_[i] = -uvwOutFrame.getValue();
antennaUVW_[i] = uvwOutFrame.getValue();
void ROCTIter::attach() {
if (iROCTMainCols_!=NULL) delete iROCTMainCols_;
if (inct_!=NULL) delete inct_;
inct_= new NewCalTable(ti_->table());
iROCTMainCols_ = new ROCTMainColumns(*inct_);
CTIter::CTIter(NewCalTable tab, const Block<String>& sortcol) :