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         258 : SingleDishBeamUtil::SingleDishBeamUtil(const MeasurementSet &ms,
      25             :                                        const String &referenceFrame,
      26             :                                        const String &movingSource,
      27             :                                        const String &pointingColumn,
      28         258 :                                        const String &antenna)
      29         258 :        : referenceFrame_(referenceFrame), movingSource_(movingSource),
      30         516 :          pointingColumn_(pointingColumn), antSel_(antenna)
      31             : {
      32         258 :   ms_ = new MeasurementSet(ms);
      33         258 :   directionUnit_ = Unit("rad");
      34         258 : }
      35             : 
      36         258 : Bool SingleDishBeamUtil::getMapPointings(Matrix<Double> &pointingList) {
      37             :     try {
      38         258 :         PointingDirectionCalculator calc(*ms_);
      39             : 
      40         258 :         calc.setDirectionColumn(pointingColumn_);
      41         258 :         calc.selectData(antSel_);
      42         258 :         calc.setFrame(referenceFrame_);
      43         258 :         MDirection::Types refType = MDirection::J2000; // any non planet value
      44         258 :         Bool status = False;
      45         258 :         status = MDirection::getType(refType, movingSource_);
      46         519 :         Bool doMovingSourceCorrection = (status == True &&
      47         261 :                 MDirection::N_Types < refType &&
      48           3 :                 refType < MDirection::N_Planets);
      49         258 :         if (doMovingSourceCorrection) {
      50           3 :             calc.setMovingSource(movingSource_);
      51             :         }
      52         258 :         calc.setDirectionListMatrixShape(PointingDirectionCalculator::COLUMN_MAJOR);
      53             : 
      54         258 :         pointingList = calc.getDirection();
      55         258 :         directionUnit_ = Unit("rad");
      56         258 :         Vector<Double> longitude = pointingList.column(0);
      57         258 :         Vector<Double> latitude = pointingList.column(1);
      58             : 
      59         258 :         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         222 :         Double longitudeMean = mean(longitude);
      66         222 :         Double longitudeStddev = stddev(longitude);
      67         222 :         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         330 :     }
      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         222 :     return True;
      99             : }
     100             : 
     101         258 :   Bool SingleDishBeamUtil::getPointingSamplingRaster(Quantum<Vector<Double>> &sampling, Quantity &positionAngle) {
     102         516 :     LogIO os(_ORIGIN);
     103         258 :     os << "calculating sampling interval assuming raster scan." << LogIO::POST;
     104             :     try {
     105         258 :       Vector<Double> samplingVal(2, 0.0);
     106             :       // Get time sorted pointing (sort by ANTENNA, TIME)
     107         258 :       Matrix<Double> pointingList;
     108         258 :       ThrowIf (!getMapPointings(pointingList), "Failed to get map pointings");
     109         258 :       os << "got " << pointingList.column(0).size() << " pointings of " << antSel_  << LogIO::POST;
     110             :       // Get timestamps
     111         258 :       Block<String> sortColumns(2);
     112         258 :       sortColumns[0] = "ANTENNA1";
     113         258 :       sortColumns[1] = "TIME";
     114         258 :       MSSelection thisSelection;
     115         258 :       thisSelection.setAntennaExpr(antSel_);
     116         258 :       TableExprNode exprNode = thisSelection.getTEN(&(*ms_));
     117         258 :       if (exprNode.isNull()) {
     118           0 :         throw(AipsError("Invalid antenna selection"));
     119             :       }
     120         258 :       MeasurementSet tmp = (*ms_)(exprNode);
     121         258 :       CountedPtr<MeasurementSet> sortedMS = new MeasurementSet(tmp.sort(sortColumns));
     122         258 :       AlwaysAssert(sortedMS->nrow() == pointingList.column(0).size(), AipsError);
     123         258 :       ScalarMeasColumn<MEpoch> timeColumn_;
     124         258 :       timeColumn_.attach(*sortedMS, "TIME");
     125             :       // Get time and pointings of unique time stamp per antenna
     126         258 :       Vector<uInt> uniqueAntTimeIdx(sortedMS->nrow());
     127         258 :       Vector<Double> uniqueAntTimes(sortedMS->nrow());
     128         258 :       Double currentTime = timeColumn_.convert(0, MEpoch::UTC).get("s").getValue();
     129         258 :       uInt itime = 0;
     130             :       // Initial time stamp
     131         258 :       uniqueAntTimeIdx(itime) = 0;
     132         258 :       uniqueAntTimes(itime) = currentTime;
     133         258 :       ++itime;
     134      344658 :       for (uInt i = 1; i < sortedMS->nrow(); ++i) {
     135      344400 :         Double nextTime = timeColumn_.convert(i, MEpoch::UTC).get("s").getValue();
     136      344400 :         if (abs(nextTime-currentTime) > 1.e-14*currentTime) {
     137      330068 :           uniqueAntTimeIdx(itime) = i;
     138      330068 :           uniqueAntTimes(itime) = nextTime;
     139      330068 :           ++itime;
     140             :         }
     141      344400 :         currentTime = nextTime;
     142             :       }
     143         258 :       if (itime != uniqueAntTimeIdx.size()) {
     144          34 :         uniqueAntTimeIdx.resize(itime, True);
     145          34 :         uniqueAntTimes.resize(itime, True); 
     146             :       }
     147         258 :       os << LogIO::DEBUG1 << uniqueAntTimeIdx.size() << " unique time stamps"  << LogIO::POST;
     148         258 :       if (uniqueAntTimes.size() == 1) {
     149          36 :         samplingVal = 0.0;
     150          36 :         sampling = Quantum< Vector<Double> >(samplingVal, directionUnit_);
     151          36 :         positionAngle = Quantity(0.0, "rad");
     152             :         os << LogIO::NORMAL
     153             :            << "Got only one pointing. exiting without calculating sampling interval."
     154          36 :            << LogIO::POST;
     155          36 :         return True;
     156             :       }
     157         222 :       Vector<Double> longitude(uniqueAntTimes.size());
     158         222 :       Vector<Double> latitude(longitude.size());
     159             :       {
     160         222 :         Vector<Double> all_longitude = pointingList.column(0); // time sorted pointing list
     161         222 :         Vector<Double> all_latitude = pointingList.column(1);
     162      330512 :         for (uint i = 0; i < uniqueAntTimeIdx.size(); ++i) {
     163      330290 :           longitude(i) = all_longitude[uniqueAntTimeIdx[i]];
     164      330290 :           latitude(i) = all_latitude[uniqueAntTimeIdx[i]];
     165             :         }
     166         222 :       }
     167             :       // calculate pointing interval assuming RASTER
     168         222 :       Vector<Double> delta_lon(longitude.size()-1);
     169         222 :       Vector<Double> delta_lat(delta_lon.size());
     170             :       Double min_lat, max_lat;
     171         222 :       minMax(min_lat, max_lat, latitude);
     172         222 :       Double center_lat = 0.5*(min_lat+max_lat);
     173      330290 :       for (size_t i=0; i < delta_lon.size(); ++i) {
     174      330068 :         Double dlon = (longitude[i+1]-longitude[i])*cos(center_lat);
     175      330068 :         delta_lon[i] = abs(dlon) > 1.e-8 ? dlon : 1.e-8;
     176      330068 :         delta_lat[i] = latitude[i+1]-latitude[i];
     177             :       }
     178             :       { // sampling interval along scan row
     179         222 :         Vector<Double> delta2(delta_lon.size());
     180         222 :         delta2 = square(delta_lon) + square(delta_lat);
     181         222 :         samplingVal(0) = sqrt(median(delta2));
     182         444 :         os << LogIO::DEBUG1 << "sampling interval along scan: " << samplingVal(0)
     183         222 :            << " " << directionUnit_.getName() << LogIO::POST;
     184         222 :       }
     185             :       { // position angle
     186         222 :         Vector<Double> delta_tan(delta_lon.size());
     187         222 :         delta_tan = delta_lat/delta_lon;
     188         222 :         Double positionAngleVal = atan(median(delta_tan));
     189         222 :         positionAngleVal = std::isfinite(positionAngleVal) ? positionAngleVal: 0.5*C::pi;
     190         222 :         positionAngle = Quantity(positionAngleVal, "rad");
     191         222 :         os << LogIO::DEBUG1 << "position angle of scan direction: " << positionAngle << LogIO::POST;
     192         222 :       }
     193             :       { // sampling interval orthogonal to scan row
     194         222 :         vector<uInt> gap_idx(0);
     195         222 :         uInt numAntGap = 0;
     196         222 :         os << LogIO::DEBUG1 << "start analysing raster pattern by time gap " << LogIO::POST;
     197             :         {// detect raster gap by time interval
     198         222 :           Vector<Double> deltaTime(uniqueAntTimes.size()-1);
     199         222 :           Vector<Double> positiveTimeGap(deltaTime.size());
     200         222 :           uInt itime = 0;
     201      330290 :           for (uInt i = 0; i < deltaTime.size(); ++i) {
     202      330068 :             deltaTime[i] = uniqueAntTimes[i+1] - uniqueAntTimes[i];
     203      330068 :             if (deltaTime[i] > 0.0) {
     204      330068 :               positiveTimeGap[itime] = deltaTime[i];
     205      330068 :               ++itime;
     206             :             }
     207             :           }
     208         222 :           positiveTimeGap.resize(itime, True);
     209         222 :           Double medianInterval5 = Double(5)*median(positiveTimeGap);
     210         222 :           os << LogIO::DEBUG1 << "Gap interval threshold = " << medianInterval5  << LogIO::POST;
     211      330290 :           for (uInt i = 0; i < deltaTime.size(); ++i) {
     212      330068 :             if (deltaTime[i] > medianInterval5) {// raster row gap
     213        3730 :               gap_idx.push_back(i);
     214      326338 :             } else if (deltaTime[i] < 0.0 || i == deltaTime.size()-1) { // antenna gap
     215         222 :               gap_idx.push_back(i);
     216         222 :               ++numAntGap;
     217             :             }
     218             :           }
     219         222 :         }
     220         222 :         if (gap_idx.size()==numAntGap) {// no gap detected.
     221         105 :           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         105 :           samplingVal(1) = 0.0;
     223         105 :           sampling = Quantum< Vector<Double> >(samplingVal, directionUnit_);
     224         105 :           os << "sampling interval: " << sampling << ", pa: " << positionAngle  << LogIO::POST;
     225         105 :           return True;
     226             :         }
     227         117 :         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         117 :         Vector<Double> orth_vec(2);
     231         117 :         orth_vec(0) = cos(positionAngle.getValue("rad")+0.5*C::pi);
     232         117 :         orth_vec(1) = sin(positionAngle.getValue("rad")+0.5*C::pi);
     233         117 :         os << LogIO::DEBUG1 << "orthogonal vector = " << orth_vec << LogIO::POST;
     234             :         // median lon, lat in each raster row
     235         117 :         Vector<Double> typical_lon(gap_idx.size());
     236         117 :         Vector<Double> typical_lat(gap_idx.size());
     237        3964 :         for (uInt i = 0; i < gap_idx.size(); ++i) {
     238             :           uInt start_idx, num_dump;
     239        3847 :           if (i==0) start_idx = 0;
     240        3730 :           else start_idx = gap_idx[i-1] + 1;
     241        3847 :           num_dump = gap_idx[i]-start_idx + 1;
     242        3847 :           typical_lon(i) = median(longitude(Slice(start_idx, num_dump)));
     243        3847 :           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         117 :            << 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         117 :            << typical_lat(Slice(0, min(typical_lat.size(), 10))) << LogIO::POST;
     249         117 :         Vector<Double> orth_dist(typical_lon.size()-1);
     250        3847 :         for (uInt i = 0; i < typical_lon.size()-1; ++i) {
     251             :           Double delta_row_lon, delta_row_lat;
     252        3730 :           delta_row_lon = typical_lon(i) - typical_lon(i+1);
     253        3730 :           delta_row_lat = typical_lat(i) - typical_lat(i+1);
     254             :           // product of orthogonal vector and row gap vector
     255        7460 :           orth_dist(i) = abs(delta_row_lon*orth_vec(0) + 
     256        3730 :                              delta_row_lat*orth_vec(1));
     257             :         }
     258         117 :         samplingVal(1) = median(orth_dist);
     259         234 :         os << LogIO::DEBUG1 << "sampling interval between scan row: " << samplingVal(1)
     260         117 :            << " " << directionUnit_.getName() << LogIO::POST;
     261         222 :       }
     262         117 :       sampling = Quantum<Vector<Double>>(samplingVal, directionUnit_);
     263        1947 :     }
     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         117 :        << ", pa: " << positionAngle  << LogIO::POST;
     275         117 :     return True;
     276         258 : }
     277             : 
     278             : } //# NAMESPACE CASA - END