casa6

#include <cassert>

#include <memory>

#include <iostream>

​

#include <casacore/casa/Arrays/ArrayMath.h>

#include <casacore/casa/Arrays/Vector.h>

#include <casacore/casa/Utilities/Assert.h>

#include <casacore/casa/Logging.h>

#include <casacore/measures/Measures/MDirection.h>

#include <casacore/measures/TableMeasures/ScalarMeasColumn.h>

#include <casacore/ms/MSSel/MSSelection.h>

#include <casacore/tables/TaQL/ExprNode.h>

​

#include <synthesis/Utilities/PointingDirectionCalculator.h>

#include <synthesis/Utilities/SingleDishBeamUtil.h>

​

using namespace std;

using namespace casacore;

​

#define _ORIGIN LogOrigin("SingleDishBeamUtil", __func__, WHERE)

​

namespace casa { //# NAMESPACE CASA - BEGIN

​

SingleDishBeamUtil::SingleDishBeamUtil(const MeasurementSet &ms,

                       const String &referenceFrame,

                       const String &movingSource,

                       const String &pointingColumn,

                       const String &antenna)

       : referenceFrame_(referenceFrame), movingSource_(movingSource),

     pointingColumn_(pointingColumn), antSel_(antenna)

{

  ms_ = new MeasurementSet(ms);

  directionUnit_ = Unit("rad");

}

​

Bool SingleDishBeamUtil::getMapPointings(Matrix<Double> &pointingList) {

    try {

        PointingDirectionCalculator calc(*ms_);

​

        calc.setDirectionColumn(pointingColumn_);

    calc.selectData(antSel_);

        calc.setFrame(referenceFrame_);

        MDirection::Types refType = MDirection::J2000; // any non planet value

        Bool status = False;

        status = MDirection::getType(refType, movingSource_);

        Bool doMovingSourceCorrection = (status == True &&

                MDirection::N_Types < refType &&

                refType < MDirection::N_Planets);

        if (doMovingSourceCorrection) {

            calc.setMovingSource(movingSource_);

        }

        calc.setDirectionListMatrixShape(PointingDirectionCalculator::COLUMN_MAJOR);

​

        pointingList = calc.getDirection();

    directionUnit_ = Unit("rad");

        Vector<Double> longitude = pointingList.column(0);

        Vector<Double> latitude = pointingList.column(1);

​

    if (longitude.size() < 2) return True; // no need for boundary check.

​

        // Diagnose if longitude values are divided by periodic boundary surface

        // (+-pi or 0, 2pi)

        // In this case, mean of longitude should be around 0 (+-pi) or pi (0,2pi)