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