Line data    Source code

       1             : #include <cassert>
       2             : #include <memory>
       3             : #include <iostream>
       4             : 
       5             : #include <casacore/casa/Arrays/ArrayMath.h>
       6             : #include <casacore/casa/Arrays/Vector.h>
       7             : #include <casacore/casa/Utilities/Assert.h>
       8             : #include <casacore/casa/Logging.h>
       9             : #include <casacore/measures/Measures/MDirection.h>
      10             : #include <casacore/measures/TableMeasures/ScalarMeasColumn.h>
      11             : #include <casacore/ms/MSSel/MSSelection.h>
      12             : #include <casacore/tables/TaQL/ExprNode.h>
      13             : 
      14             : #include <synthesis/Utilities/PointingDirectionCalculator.h>
      15             : #include <synthesis/Utilities/SingleDishBeamUtil.h>
      16             : 
      17             : using namespace std;
      18             : using namespace casacore;
      19             : 
      20             : #define _ORIGIN LogOrigin("SingleDishBeamUtil", __func__, WHERE)
      21             : 
      22             : namespace casa { //# NAMESPACE CASA - BEGIN
      23             : 
      24          12 : SingleDishBeamUtil::SingleDishBeamUtil(const MeasurementSet &ms,
      25             :                                        const String &referenceFrame,
      26             :                                        const String &movingSource,
      27             :                                        const String &pointingColumn,
      28          12 :                                        const String &antenna)
      29          12 :        : referenceFrame_(referenceFrame), movingSource_(movingSource),
      30          24 :          pointingColumn_(pointingColumn), antSel_(antenna)
      31             : {
      32          12 :   ms_ = new MeasurementSet(ms);
      33          12 :   directionUnit_ = Unit("rad");
      34          12 : }
      35             : 
      36          12 : Bool SingleDishBeamUtil::getMapPointings(Matrix<Double> &pointingList) {
      37             :     try {
      38          12 :         PointingDirectionCalculator calc(*ms_);
      39             : 
      40          12 :         calc.setDirectionColumn(pointingColumn_);
      41          12 :         calc.selectData(antSel_);
      42          12 :         calc.setFrame(referenceFrame_);
      43          12 :         MDirection::Types refType = MDirection::J2000; // any non planet value
      44          12 :         Bool status = False;
      45          12 :         status = MDirection::getType(refType, movingSource_);
      46          24 :         Bool doMovingSourceCorrection = (status == True &&
      47          12 :                 MDirection::N_Types < refType &&
      48           0 :                 refType < MDirection::N_Planets);
      49          12 :         if (doMovingSourceCorrection) {
      50           0 :             calc.setMovingSource(movingSource_);
      51             :         }
      52          12 :         calc.setDirectionListMatrixShape(PointingDirectionCalculator::COLUMN_MAJOR);
      53             : 
      54          12 :         pointingList = calc.getDirection();
      55          12 :         directionUnit_ = Unit("rad");
      56          12 :         Vector<Double> longitude = pointingList.column(0);
      57          12 :         Vector<Double> latitude = pointingList.column(1);
      58             : 
      59          12 :         if (longitude.size() < 2) return True; // no need for boundary check.
      60             : 
      61             :         // Diagnose if longitude values are divided by periodic boundary surface
      62             :         // (+-pi or 0, 2pi)
      63             :         // In this case, mean of longitude should be around 0 (+-pi) or pi (0,2pi)
      64             :         // and stddev of longitude array be around pi.
      65          12 :         Double longitudeMean = mean(longitude);
      66          12 :         Double longitudeStddev = stddev(longitude);
      67          12 :         if (longitudeStddev > 2.0 / 3.0 * C::pi) {
      68             :             // likely to be the case
      69           0 :             if (abs(longitudeMean) < 0.5 * C::pi) {
      70             :                 // periodic boundary surface would be +-pi
      71           0 :                 for (size_t i = 0; i < longitude.nelements(); ++i) {
      72           0 :                     if (longitude[i] < 0.0) {
      73           0 :                         longitude[i] += C::_2pi;
      74             :                     }
      75             :                 }
      76             :             }
      77           0 :             else if (abs(longitudeMean - C::pi) < 0.5 * C::pi ) {
      78             :                 // periodic boundary surface would be 0,2pi
      79           0 :                 for (size_t i = 0; i < longitude.nelements(); ++i) {
      80           0 :                     if (longitude[i] < C::pi) {
      81           0 :                         longitude[i] += C::_2pi;
      82             :                     }
      83             :                 }
      84             :             }
      85             :         }
      86          12 :     }
      87           0 :     catch (AipsError &e) {
      88           0 :         LogIO os(LogOrigin("Imager", "getMapPointings", WHERE));
      89           0 :         os << LogIO::SEVERE << "Failed due to the rror \"" << e.getMesg() << "\"" << LogIO::POST;
      90           0 :         return False;
      91           0 :     }
      92           0 :     catch (...) {
      93           0 :         LogIO os(LogOrigin("Imager", "getMapPointings", WHERE));
      94           0 :         os << LogIO::SEVERE << "Failed due to unknown error" << LogIO::POST;
      95           0 :         throw;
      96             :         return False;
      97           0 :     }
      98          12 :     return True;
      99             : }
     100             : 
     101          12 :   Bool SingleDishBeamUtil::getPointingSamplingRaster(Quantum<Vector<Double>> &sampling, Quantity &positionAngle) {
     102          24 :     LogIO os(_ORIGIN);
     103          12 :     os << "calculating sampling interval assuming raster scan." << LogIO::POST;
     104             :     try {
     105          12 :       Vector<Double> samplingVal(2, 0.0);
     106             :       // Get time sorted pointing (sort by ANTENNA, TIME)
     107          12 :       Matrix<Double> pointingList;
     108          12 :       ThrowIf (!getMapPointings(pointingList), "Failed to get map pointings");
     109          12 :       os << "got " << pointingList.column(0).size() << " pointings of " << antSel_  << LogIO::POST;
     110             :       // Get timestamps
     111          12 :       Block<String> sortColumns(2);
     112          12 :       sortColumns[0] = "ANTENNA1";
     113          12 :       sortColumns[1] = "TIME";
     114          12 :       MSSelection thisSelection;
     115          12 :       thisSelection.setAntennaExpr(antSel_);
     116          12 :       TableExprNode exprNode = thisSelection.getTEN(&(*ms_));
     117          12 :       if (exprNode.isNull()) {
     118           0 :         throw(AipsError("Invalid antenna selection"));
     119             :       }
     120          12 :       MeasurementSet tmp = (*ms_)(exprNode);
     121          12 :       CountedPtr<MeasurementSet> sortedMS = new MeasurementSet(tmp.sort(sortColumns));
     122          12 :       AlwaysAssert(sortedMS->nrow() == pointingList.column(0).size(), AipsError);
     123          12 :       ScalarMeasColumn<MEpoch> timeColumn_;
     124          12 :       timeColumn_.attach(*sortedMS, "TIME");
     125             :       // Get time and pointings of unique time stamp per antenna
     126          12 :       Vector<uInt> uniqueAntTimeIdx(sortedMS->nrow());
     127          12 :       Vector<Double> uniqueAntTimes(sortedMS->nrow());
     128          12 :       Double currentTime = timeColumn_.convert(0, MEpoch::UTC).get("s").getValue();
     129          12 :       uInt itime = 0;
     130             :       // Initial time stamp
     131          12 :       uniqueAntTimeIdx(itime) = 0;
     132          12 :       uniqueAntTimes(itime) = currentTime;
     133          12 :       ++itime;
     134       11264 :       for (uInt i = 1; i < sortedMS->nrow(); ++i) {
     135       11252 :         Double nextTime = timeColumn_.convert(i, MEpoch::UTC).get("s").getValue();
     136       11252 :         if (abs(nextTime-currentTime) > 1.e-14*currentTime) {
     137        9208 :           uniqueAntTimeIdx(itime) = i;
     138        9208 :           uniqueAntTimes(itime) = nextTime;
     139        9208 :           ++itime;
     140             :         }
     141       11252 :         currentTime = nextTime;
     142             :       }
     143          12 :       if (itime != uniqueAntTimeIdx.size()) {
     144           4 :         uniqueAntTimeIdx.resize(itime, True);
     145           4 :         uniqueAntTimes.resize(itime, True); 
     146             :       }
     147          12 :       os << LogIO::DEBUG1 << uniqueAntTimeIdx.size() << " unique time stamps"  << LogIO::POST;
     148          12 :       if (uniqueAntTimes.size() == 1) {
     149           0 :         samplingVal = 0.0;
     150           0 :         sampling = Quantum< Vector<Double> >(samplingVal, directionUnit_);
     151           0 :         positionAngle = Quantity(0.0, "rad");
     152             :         os << LogIO::NORMAL
     153             :            << "Got only one pointing. exiting without calculating sampling interval."
     154           0 :            << LogIO::POST;
     155           0 :         return True;
     156             :       }
     157          12 :       Vector<Double> longitude(uniqueAntTimes.size());
     158          12 :       Vector<Double> latitude(longitude.size());
     159             :       {
     160          12 :         Vector<Double> all_longitude = pointingList.column(0); // time sorted pointing list
     161          12 :         Vector<Double> all_latitude = pointingList.column(1);
     162        9232 :         for (uint i = 0; i < uniqueAntTimeIdx.size(); ++i) {
     163        9220 :           longitude(i) = all_longitude[uniqueAntTimeIdx[i]];
     164        9220 :           latitude(i) = all_latitude[uniqueAntTimeIdx[i]];
     165             :         }
     166          12 :       }
     167             :       // calculate pointing interval assuming RASTER
     168          12 :       Vector<Double> delta_lon(longitude.size()-1);
     169          12 :       Vector<Double> delta_lat(delta_lon.size());
     170             :       Double min_lat, max_lat;
     171          12 :       minMax(min_lat, max_lat, latitude);
     172          12 :       Double center_lat = 0.5*(min_lat+max_lat);
     173        9220 :       for (size_t i=0; i < delta_lon.size(); ++i) {
     174        9208 :         Double dlon = (longitude[i+1]-longitude[i])*cos(center_lat);
     175        9208 :         delta_lon[i] = abs(dlon) > 1.e-8 ? dlon : 1.e-8;
     176        9208 :         delta_lat[i] = latitude[i+1]-latitude[i];
     177             :       }
     178             :       { // sampling interval along scan row
     179          12 :         Vector<Double> delta2(delta_lon.size());
     180          12 :         delta2 = square(delta_lon) + square(delta_lat);
     181          12 :         samplingVal(0) = sqrt(median(delta2));
     182          24 :         os << LogIO::DEBUG1 << "sampling interval along scan: " << samplingVal(0)
     183          12 :            << " " << directionUnit_.getName() << LogIO::POST;
     184          12 :       }
     185             :       { // position angle
     186          12 :         Vector<Double> delta_tan(delta_lon.size());
     187          12 :         delta_tan = delta_lat/delta_lon;
     188          12 :         Double positionAngleVal = atan(median(delta_tan));
     189          12 :         positionAngleVal = std::isfinite(positionAngleVal) ? positionAngleVal: 0.5*C::pi;
     190          12 :         positionAngle = Quantity(positionAngleVal, "rad");
     191          12 :         os << LogIO::DEBUG1 << "position angle of scan direction: " << positionAngle << LogIO::POST;
     192          12 :       }
     193             :       { // sampling interval orthogonal to scan row
     194          12 :         vector<uInt> gap_idx(0);
     195          12 :         uInt numAntGap = 0;
     196          12 :         os << LogIO::DEBUG1 << "start analysing raster pattern by time gap " << LogIO::POST;
     197             :         {// detect raster gap by time interval
     198          12 :           Vector<Double> deltaTime(uniqueAntTimes.size()-1);
     199          12 :           Vector<Double> positiveTimeGap(deltaTime.size());
     200          12 :           uInt itime = 0;
     201        9220 :           for (uInt i = 0; i < deltaTime.size(); ++i) {
     202        9208 :             deltaTime[i] = uniqueAntTimes[i+1] - uniqueAntTimes[i];
     203        9208 :             if (deltaTime[i] > 0.0) {
     204        9208 :               positiveTimeGap[itime] = deltaTime[i];
     205        9208 :               ++itime;
     206             :             }
     207             :           }
     208          12 :           positiveTimeGap.resize(itime, True);
     209          12 :           Double medianInterval5 = Double(5)*median(positiveTimeGap);
     210          12 :           os << LogIO::DEBUG1 << "Gap interval threshold = " << medianInterval5  << LogIO::POST;
     211        9220 :           for (uInt i = 0; i < deltaTime.size(); ++i) {
     212        9208 :             if (deltaTime[i] > medianInterval5) {// raster row gap
     213          24 :               gap_idx.push_back(i);
     214        9184 :             } else if (deltaTime[i] < 0.0 || i == deltaTime.size()-1) { // antenna gap
     215          12 :               gap_idx.push_back(i);
     216          12 :               ++numAntGap;
     217             :             }
     218             :           }
     219          12 :         }
     220          12 :         if (gap_idx.size()==numAntGap) {// no gap detected.
     221           4 :           os << LogIO::NORMAL << "No time gap found in scans. The scan pattern may not be RASTER. Median sampling interval will be returned." << LogIO::POST;
     222           4 :           samplingVal(1) = 0.0;
     223           4 :           sampling = Quantum< Vector<Double> >(samplingVal, directionUnit_);
     224           4 :           os << "sampling interval: " << sampling << ", pa: " << positionAngle  << LogIO::POST;
     225           4 :           return True;
     226             :         }
     227           8 :         os << LogIO::DEBUG1 << gap_idx.size() << " raster rows detected" << LogIO::POST;
     228             :         //os << LogIO::DEBUGGING << "gap idx = " << Vector<uInt>(gap_idx) << LogIO::POST;
     229             :         // a unit vector of orthogonal direction
     230           8 :         Vector<Double> orth_vec(2);
     231           8 :         orth_vec(0) = cos(positionAngle.getValue("rad")+0.5*C::pi);
     232           8 :         orth_vec(1) = sin(positionAngle.getValue("rad")+0.5*C::pi);
     233           8 :         os << LogIO::DEBUG1 << "orthogonal vector = " << orth_vec << LogIO::POST;
     234             :         // median lon, lat in each raster row
     235           8 :         Vector<Double> typical_lon(gap_idx.size());
     236           8 :         Vector<Double> typical_lat(gap_idx.size());
     237          40 :         for (uInt i = 0; i < gap_idx.size(); ++i) {
     238             :           uInt start_idx, num_dump;
     239          32 :           if (i==0) start_idx = 0;
     240          24 :           else start_idx = gap_idx[i-1] + 1;
     241          32 :           num_dump = gap_idx[i]-start_idx + 1;
     242          32 :           typical_lon(i) = median(longitude(Slice(start_idx, num_dump)));
     243          32 :           typical_lat(i) = median(latitude(Slice(start_idx, num_dump)));
     244             :         }
     245             :         os << LogIO::DEBUG1 << "Typical longitude of scan row (first 10 at max) = "
     246           8 :            << typical_lon(Slice(0, min(typical_lon.size(), 10))) << LogIO::POST;
     247             :         os << LogIO::DEBUG1 << "Typical latitude of scan row (first 10 at max) = "
     248           8 :            << typical_lat(Slice(0, min(typical_lat.size(), 10))) << LogIO::POST;
     249           8 :         Vector<Double> orth_dist(typical_lon.size()-1);
     250          32 :         for (uInt i = 0; i < typical_lon.size()-1; ++i) {
     251             :           Double delta_row_lon, delta_row_lat;
     252          24 :           delta_row_lon = typical_lon(i) - typical_lon(i+1);
     253          24 :           delta_row_lat = typical_lat(i) - typical_lat(i+1);
     254             :           // product of orthogonal vector and row gap vector
     255          48 :           orth_dist(i) = abs(delta_row_lon*orth_vec(0) + 
     256          24 :                              delta_row_lat*orth_vec(1));
     257             :         }
     258           8 :         samplingVal(1) = median(orth_dist);
     259          16 :         os << LogIO::DEBUG1 << "sampling interval between scan row: " << samplingVal(1)
     260           8 :            << " " << directionUnit_.getName() << LogIO::POST;
     261          12 :       }
     262           8 :       sampling = Quantum<Vector<Double>>(samplingVal, directionUnit_);
     263          64 :     }
     264           0 :     catch (AipsError &e) {
     265           0 :         os << LogIO::SEVERE << "Failed due to the rror \"" << e.getMesg() << "\"" << LogIO::POST;
     266           0 :         return False;
     267           0 :     }
     268           0 :     catch (...) {
     269           0 :         os << LogIO::SEVERE << "Failed due to unknown error" << LogIO::POST;
     270           0 :         throw;
     271             :         return False;
     272           0 :     }
     273             :     os << LogIO::NORMAL << "sampling interval: " << sampling
     274           8 :        << ", pa: " << positionAngle  << LogIO::POST;
     275           8 :     return True;
     276          12 : }
     277             : 
     278             : } //# NAMESPACE CASA - END