casa6

// Based on code/alma/apps/UVWCoords

​

#include <casacore/casa/Logging/LogIO.h>

#include <msvis/MSVis/MSUVWGenerator.h>

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

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

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

#include <casacore/ms/MeasurementSets/MSColumns.h>

#include <casacore/ms/MeasurementSets/MSAntennaColumns.h>

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

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

#include <map>

​

using namespace casacore;

namespace casa {

​

// The UvwCoords ctor has lines for the antennas, antenna offsets, and station

// positions.  This ctor assumes they're present in msc_p if present at all.

  MSUVWGenerator::MSUVWGenerator(MSColumns &msc_ref, const MBaseline::Types bltype,

                 const Muvw::Types) :

  msc_p(msc_ref),                       

  bl_csys_p(MBaseline::Ref(bltype)), // MBaseline::J2000, ITRF, etc.

  antColumns_p(msc_p.antenna()),

  antPositions_p(antColumns_p.positionMeas()),

  antOffset_p(antColumns_p.offsetMeas()),

  refpos_p(antPositions_p(0)),  // We use the first antenna for the reference

  feedOffset_p(msc_p.feed().positionMeas())

{

  // It seems that using a String is the only safe way to go from say,

  // MPosition::ITRF to MBaseline::ITRF.

  MBaseline::getType(refposref_p,

                     MPosition::showType(refpos_p.getRef().getType()));

​

  fill_bl_an(bl_an_p);      

}

​

MSUVWGenerator::~MSUVWGenerator(){

}

​

void MSUVWGenerator::fill_bl_an(Vector<MVBaseline>& bl_an_p)

{

  nant_p = antPositions_p.table().nrow();

  logSink() << LogIO::DEBUG1 << "nant_p: " << nant_p << LogIO::POST;

​

  Double max_baseline = -1.0;

  Double bl_len;

  const ScalarColumn<Double>& antDiams = antColumns_p.dishDiameter();

  Double smallestDiam = antDiams(0);

  Double secondSmallestDiam = antDiams(0);

  bl_an_p.resize(nant_p);

  for(uInt an = 0; an < nant_p; ++an){

    // MVBaselines are basically xyz Vectors, not Measures.

    bl_an_p[an] = MVBaseline(refpos_p.getValue(), antPositions_p(an).getValue());

​

    // MVBaseline has functions to return the length, but Manhattan distances

    // are good enough for this, and faster than a sqrt.

    Vector<Double> bluvw(bl_an_p[an].getValue());

    bl_len = fabs(bluvw[0]) + fabs(bluvw[1]) + fabs(bluvw[2]);

​

    if(bl_len > max_baseline)

      max_baseline = bl_len;

​

    if(antDiams(an) < secondSmallestDiam){

      if(antDiams(an) < smallestDiam){

    secondSmallestDiam = smallestDiam;

    smallestDiam = antDiams(an);

      }

      else

    secondSmallestDiam = antDiams(an);

    }

  }

​

  // Setting timeRes_p to 0.0025 * the time for a 1 radian phase change on the

  // longest baseline at 2x the primary beamwidth should be sufficiently short

  // for Earth based observations.  Space-based baselines will move faster, but

  // probably don't have the data rate to support full beam imaging.  An

  // alternative limit could come from |d(uvw)/dt|/|uvw|, but that is

  // guaranteed to be at least somewhat larger than this limit (replace the

  // 2.44 with 2, and remove the dish diameters and max_baseline).

  //

  // Do not raise the 0.0025 coefficient by too much, since the times used for

  // UVW calculation could be biased by as much as -0.5 * timeRes_p if more

  // than one integration fits into timeRes_p.  timeRes_p is not intended so

  // much for skipping calculations as it is for preventing antUVW_p from being

  // calculated for each row in the same integration.  The integration interval

  // may change within an MS, but ideally antUVW_p is calculated once per

  // integration (timeRes_p <~ integration) and there is no bias.

  timeRes_p = 0.0025 * 24.0 * 3600.0 / (6.283 * 2.44) *

    sqrt(smallestDiam * secondSmallestDiam) / max_baseline;