Line data    Source code

       1             : //# SDPosInterpolator.cc: Implementation of SDPosInterpolator class
       2             : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: casa-feedback@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : //# $Id$
      27             : 
      28             : #include <synthesis/Utilities/SDPosInterpolator.h>
      29             : 
      30             : using namespace casacore;
      31             : namespace casa {
      32             : 
      33           0 : SDPosInterpolator::SDPosInterpolator(
      34             :   const VisBuffer& vb,
      35           0 :   const String& pointingDirCol_p) {
      36           0 :   const auto & pointingColumns = vb.msColumns().pointing();
      37           0 :   const auto nant = static_cast<size_t>(vb.msColumns().antenna().nrow());
      38           0 :   setup(pointingColumns, pointingDirCol_p, nant);
      39           0 : }
      40           0 : SDPosInterpolator::SDPosInterpolator(
      41             :   const vi::VisBuffer2& vb,
      42           0 :   const String& pointingDirCol_p) {
      43           0 :   const auto & pointingColumns = vb.subtableColumns().pointing();
      44           0 :   const auto nant = static_cast<size_t>(vb.subtableColumns().antenna().nrow());
      45           0 :   setup(pointingColumns, pointingDirCol_p, nant);
      46           0 : }
      47           0 : SDPosInterpolator::SDPosInterpolator(
      48             :   const MSPointing& pointingTable,
      49             :   const String& columnName,
      50           0 :   const size_t nant){
      51           0 :   MSPointingColumns pointingColumns{pointingTable};
      52           0 :   setup(pointingColumns, columnName, nant);
      53           0 : }
      54           0 : SDPosInterpolator::SDPosInterpolator(
      55             :   const MSPointingColumns& pointingColumns,
      56             :   const String& columnName,
      57           0 :   const size_t nant){
      58           0 :   setup(pointingColumns, columnName, nant);
      59           0 : }
      60           0 : SDPosInterpolator::SDPosInterpolator(
      61             :   const Vector<Vector<Double> >& time,
      62           0 :   const Vector<Vector<Vector<Double> > >& dir) {
      63           0 :   setup(time, dir);
      64           0 : }
      65           0 : SDPosInterpolator::~SDPosInterpolator() {}
      66             : 
      67           0 : void SDPosInterpolator::setup(
      68             :       const Vector<Vector<Double> >& time,
      69             :       const Vector<Vector<Vector<Double> > >& dir) {
      70             :   //(1)get number of pointing data for each antennaID
      71           0 :   Int nant = time.nelements();
      72           0 :   Vector<uInt> nPointingData(nant);
      73           0 :   nPointingData = 0;
      74           0 :   for (Int iant = 0; iant < nant; ++iant) {
      75           0 :     nPointingData(iant) = time(iant).nelements();
      76             :   }
      77             : 
      78             :   //(2)setup spline coefficients for each antenna ID
      79           0 :   timePointing.resize(nant);
      80           0 :   dirPointing.resize(nant);
      81           0 :   splineCoeff.resize(nant);
      82           0 :   doSplineInterpolation.resize(nant);
      83           0 :   doSplineInterpolation = false;
      84           0 :   for (Int i = 0; i < nant; ++i) {
      85           0 :     if (nPointingData(i) < 4) continue;
      86             :     
      87           0 :     doSplineInterpolation(i) = true;
      88           0 :     timePointing(i).resize(nPointingData(i));
      89           0 :     dirPointing(i).resize(nPointingData(i));
      90           0 :     splineCoeff(i).resize(nPointingData(i) - 1);
      91           0 :     for (uInt j = 0; j < dirPointing(i).nelements(); ++j) {
      92           0 :       dirPointing(i)(j).resize(2);
      93             :     }
      94           0 :     for (uInt j = 0; j < splineCoeff(i).nelements(); ++j) {
      95           0 :       splineCoeff(i)(j).resize(2);
      96           0 :       splineCoeff(i)(j)(0).resize(4); // x
      97           0 :       splineCoeff(i)(j)(1).resize(4); // y
      98             :     }
      99             :     
     100           0 :     Int npoi = nPointingData(i);
     101           0 :     for (Int j = 0; j < npoi; ++j) {
     102           0 :       timePointing(i)(j) = time(i)(j);
     103           0 :       for (Int k = 0; k < 2; ++k) {
     104           0 :         dirPointing(i)(j)(k) = dir(i)(j)(k);
     105             :       }
     106             :     }
     107             :       
     108           0 :     calcSplineCoeff(timePointing(i), dirPointing(i), splineCoeff(i));
     109             :   }
     110             : 
     111             :   //(3) keep time range
     112           0 :   timeRangeStart.resize(nant);
     113           0 :   timeRangeEnd.resize(nant);
     114           0 :   for (Int iant = 0; iant < nant; ++iant) {
     115           0 :     timeRangeStart(iant) = timePointing(iant)(0);
     116           0 :     timeRangeEnd(iant)   = timePointing(iant)(timePointing(iant).nelements()-1);
     117             :   }
     118           0 : }
     119             : 
     120           0 : void SDPosInterpolator::setup(
     121             :   const MSPointingColumns& act_mspc,
     122             :   const String& pointingDirCol_p,
     123             :   size_t nant) {
     124           0 :   auto check_col = [&](Bool isnull){
     125           0 :     if (isnull) {
     126           0 :       cerr << "No " << pointingDirCol_p << " column in POINTING table" << endl;
     127             :     }
     128           0 :   };
     129           0 :   std::function<Vector<Double>(Int)> get_direction;
     130             : 
     131             :   //(0)check POINTING table and set function to obtain direction data
     132           0 :   if (pointingDirCol_p == "TARGET") {
     133           0 :     get_direction = [&](Int idx){
     134           0 :       return act_mspc.targetMeas(idx).getAngle("rad").getValue();
     135           0 :     };
     136           0 :   } else if (pointingDirCol_p == "POINTING_OFFSET") {
     137           0 :     check_col(act_mspc.pointingOffsetMeasCol().isNull());
     138           0 :     get_direction = [&](Int idx){
     139           0 :       return act_mspc.pointingOffsetMeas(idx).getAngle("rad").getValue();
     140           0 :     };
     141           0 :   } else if (pointingDirCol_p == "SOURCE_OFFSET") {
     142           0 :     check_col(act_mspc.sourceOffsetMeasCol().isNull());
     143           0 :     get_direction = [&](Int idx){
     144           0 :       return act_mspc.sourceOffsetMeas(idx).getAngle("rad").getValue();
     145           0 :     };
     146           0 :   } else if (pointingDirCol_p == "ENCODER") {
     147           0 :     check_col(act_mspc.encoderMeas().isNull());
     148           0 :     get_direction = [&](Int idx){
     149           0 :       return act_mspc.encoderMeas()(idx).getAngle("rad").getValue();
     150           0 :     };
     151             :   } else {
     152           0 :     get_direction = [&](Int idx){
     153           0 :       return act_mspc.directionMeas(idx).getAngle("rad").getValue();
     154           0 :     };
     155             :   }
     156             : 
     157             :   //(1)get number of pointing data for each antennaID
     158           0 :   Vector<uInt> nPointingData(nant, 0);
     159           0 :   auto pointingRows = static_cast<size_t>(act_mspc.nrow());
     160           0 :   for (size_t row = 0; row < pointingRows ; ++row) {
     161           0 :     nPointingData(act_mspc.antennaId()(row)) += 1;
     162             :   }
     163             : 
     164             :   //(2)setup spline coefficients for each antenna ID that
     165             :   //   appear in the main table (spectral data) if there
     166             :   //   are enough number of pointing data (4 or more).
     167             :   //   in case there exists antenna ID for which not enough
     168             :   //   (i.e., 1, 2 or 3) pointing data are given, linear
     169             :   //   interpolation is applied for that antenna ID as
     170             :   //   previously done.
     171           0 :   timePointing.resize(nant);
     172           0 :   dirPointing.resize(nant);
     173           0 :   splineCoeff.resize(nant);
     174           0 :   doSplineInterpolation.resize(nant);
     175           0 :   doSplineInterpolation = false;
     176           0 :   for (uInt i = 0; i < nant; ++i) {
     177           0 :     if (nPointingData(i) < 4) continue;
     178             : 
     179           0 :     doSplineInterpolation(i) = true;
     180           0 :     timePointing(i).resize(nPointingData(i));
     181           0 :     dirPointing(i).resize(nPointingData(i));
     182           0 :     splineCoeff(i).resize(nPointingData(i) - 1);
     183           0 :     for (uInt j = 0; j < dirPointing(i).nelements(); ++j) {
     184           0 :       dirPointing(i)(j).resize(2);
     185             :     }
     186           0 :     for (uInt j = 0; j < splineCoeff(i).nelements(); ++j) {
     187           0 :       splineCoeff(i)(j).resize(2);
     188           0 :       splineCoeff(i)(j)(0).resize(4); // x
     189           0 :       splineCoeff(i)(j)(1).resize(4); // y
     190             :     }
     191             : 
     192             :     //set ptime array etc. need for spline calculation...
     193           0 :     size_t tidx = 0;
     194           0 :     for (size_t j = 0; j < pointingRows; ++j) {
     195           0 :       if (act_mspc.antennaId()(j) != i) continue;
     196             : 
     197           0 :       timePointing(i)(tidx) = act_mspc.time()(j);
     198           0 :       dirPointing(i)(tidx) = get_direction(j);
     199           0 :       tidx++;
     200             :     }
     201             : 
     202           0 :     calcSplineCoeff(timePointing(i), dirPointing(i), splineCoeff(i));
     203             :   }
     204             : 
     205             :   //(3) keep time range
     206           0 :   timeRangeStart.resize(nant);
     207           0 :   timeRangeEnd.resize(nant);
     208           0 :   for (size_t iant = 0; iant < nant; ++iant) {
     209           0 :     timeRangeStart(iant) = timePointing(iant)(0);
     210           0 :     timeRangeEnd(iant)   = timePointing(iant)(timePointing(iant).nelements()-1);
     211             :   }
     212           0 : }
     213             : 
     214           0 : void SDPosInterpolator::calcSplineCoeff(const Vector<Double>& time,
     215             :                                         const Vector<Vector<Double> >& dir,
     216             :                                         Vector<Vector<Vector<Double> > >& coeff) {
     217           0 :   Vector<Double> h, vx, vy;
     218           0 :   Vector<Double> a;
     219           0 :   Vector<Double> c;
     220           0 :   Vector<Double> alpha, beta, gamma;
     221           0 :   Vector<Double> wx, wy;
     222           0 :   Vector<Double> ux, uy;
     223             : 
     224           0 :   Int const num_data = time.nelements();
     225           0 :   h.resize(num_data-1);
     226           0 :   vx.resize(num_data-1);
     227           0 :   vy.resize(num_data-1);
     228           0 :   a.resize(num_data-1);
     229           0 :   c.resize(num_data-1);
     230           0 :   alpha.resize(num_data-1);
     231           0 :   beta.resize(num_data-1);
     232           0 :   gamma.resize(num_data-1);
     233           0 :   wx.resize(num_data-1);
     234           0 :   wy.resize(num_data-1);
     235           0 :   ux.resize(num_data);
     236           0 :   uy.resize(num_data);
     237             : 
     238           0 :   h(0) = time(1) - time(0);
     239           0 :   for (Int i = 1; i < num_data-1; ++i) {
     240           0 :     h(i) = time(i+1) - time(i);
     241           0 :     vx(i) = 6.0*((dir(i+1)(0)-dir(i)(0))/h(i) - (dir(i)(0)-dir(i-1)(0))/h(i-1));
     242           0 :     vy(i) = 6.0*((dir(i+1)(1)-dir(i)(1))/h(i) - (dir(i)(1)-dir(i-1)(1))/h(i-1));
     243           0 :     a(i) = 2.0*(time(i+1) - time(i-1));
     244           0 :     c(i) = h(i);
     245           0 :     gamma(i) = c(i);
     246             :   }
     247           0 :   alpha(2) = c(1)/a(1);
     248           0 :   for (Int i = 3; i < num_data-1; ++i) {
     249           0 :     alpha(i) = c(i-1)/(a(i-1) - alpha(i-1)*c(i-2));
     250             :   }
     251           0 :   beta(1) = a(1);
     252           0 :   for (Int i = 2; i < num_data-2; ++i) {
     253           0 :     beta(i) = c(i)/alpha(i+1);
     254             :   }
     255           0 :   beta(num_data-2) = a(num_data-2) - alpha(num_data-2) * c(num_data-3);
     256           0 :   wx(0) = 0.0;
     257           0 :   wx(1) = vx(1);
     258           0 :   wy(0) = 0.0;
     259           0 :   wy(1) = vy(1);
     260           0 :   for (Int i = 2; i < num_data-1; ++i) {
     261           0 :     wx(i) = vx(i) - alpha(i)*wx(i-1);
     262           0 :     wy(i) = vy(i) - alpha(i)*wy(i-1);
     263             :   }
     264           0 :   ux(num_data-1) = 0.0;
     265           0 :   uy(num_data-1) = 0.0;
     266           0 :   for (Int i = num_data-2; i >= 1; --i) {
     267           0 :     ux(i) = (wx(i) - gamma(i)*ux(i+1))/beta(i);
     268           0 :     uy(i) = (wy(i) - gamma(i)*uy(i+1))/beta(i);
     269             :   }
     270           0 :   ux(0) = 0.0;
     271           0 :   uy(0) = 0.0;
     272             : 
     273           0 :   for (Int i = 0; i < num_data-1; ++i) {
     274           0 :     coeff(i)(0)(0) = dir(i)(0);
     275           0 :     coeff(i)(1)(0) = dir(i)(1);
     276           0 :     const auto dt = time(i+1)-time(i);
     277           0 :     coeff(i)(0)(1) = (dir(i+1)(0)-dir(i)(0))/dt - dt*(2.0*ux(i)+ux(i+1))/6.0;
     278           0 :     coeff(i)(1)(1) = (dir(i+1)(1)-dir(i)(1))/dt - dt*(2.0*uy(i)+uy(i+1))/6.0;
     279           0 :     coeff(i)(0)(2) = ux(i)/2.0;
     280           0 :     coeff(i)(1)(2) = uy(i)/2.0;
     281           0 :     coeff(i)(0)(3) = (ux(i+1)-ux(i))/dt/6.0;
     282           0 :     coeff(i)(1)(3) = (uy(i+1)-uy(i))/dt/6.0;
     283             :   }
     284           0 : }
     285             : 
     286           0 : MDirection SDPosInterpolator::interpolateDirectionMeasSpline(
     287             :             const MSPointingColumns& mspc,
     288             :             const Double& time,
     289             :             const Int& index,
     290             :             const Int& antid) {
     291           0 :   Int lastIndex = timePointing(antid).nelements() - 1;
     292           0 :   Int aindex = lastIndex;
     293           0 :   for (uInt i = 0; i < timePointing(antid).nelements(); ++i) {
     294           0 :     if (time < timePointing(antid)(i)) {
     295           0 :       aindex = i-1;
     296           0 :       break;
     297             :     }
     298             :   }
     299           0 :   if (aindex < 0) aindex = 0;
     300           0 :   if (lastIndex <= aindex) aindex = lastIndex - 1;
     301             : 
     302           0 :   auto const &coeff = splineCoeff(antid)(aindex);
     303           0 :   Double dt = time - timePointing(antid)(aindex);
     304           0 :   Vector<Double> newdir(2);
     305             :   // Why don't we use Horner's method here ?
     306           0 :   newdir(0) =
     307           0 :     coeff(0)(0) + coeff(0)(1)*dt + coeff(0)(2)*dt*dt + coeff(0)(3)*dt*dt*dt;
     308           0 :   newdir(1) =
     309           0 :     coeff(1)(0) + coeff(1)(1)*dt + coeff(1)(2)*dt*dt + coeff(1)(3)*dt*dt*dt;
     310             : 
     311           0 :   Quantity rDirLon(newdir(0), "rad");
     312           0 :   Quantity rDirLat(newdir(1), "rad");
     313           0 :   auto const &directionMeasColumn = mspc.directionMeasCol();
     314           0 :   MDirection::Ref rf(directionMeasColumn.measDesc().getRefCode());
     315           0 :   if (directionMeasColumn.measDesc().isRefCodeVariable()) {
     316           0 :     rf = mspc.directionMeas(index).getRef();
     317             :   }
     318             : 
     319           0 :   return MDirection(rDirLon, rDirLat, rf);
     320           0 : }
     321             : 
     322           0 : Vector<Vector<Vector<Vector<Double> > > > SDPosInterpolator::getSplineCoeff() {
     323           0 :   return splineCoeff;
     324             : }
     325             : 
     326           0 : Bool SDPosInterpolator::inTimeRange(const Double& time, const Int& antid) {
     327           0 :   Bool inrange = false;
     328           0 :   if ((timeRangeStart(antid) <= time) && (time <= timeRangeEnd(antid))) {
     329           0 :     inrange = true;
     330             :   }
     331           0 :   return inrange;
     332             : }
     333             : 
     334             : } //#End casa namespace