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          21 : SDPosInterpolator::SDPosInterpolator(
      61             :   const Vector<Vector<Double> >& time,
      62          21 :   const Vector<Vector<Vector<Double> > >& dir) {
      63          21 :   setup(time, dir);
      64          21 : }
      65          21 : SDPosInterpolator::~SDPosInterpolator() {}
      66             : 
      67          21 : 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          21 :   Int nant = time.nelements();
      72          21 :   Vector<uInt> nPointingData(nant);
      73          21 :   nPointingData = 0;
      74          42 :   for (Int iant = 0; iant < nant; ++iant) {
      75          21 :     nPointingData(iant) = time(iant).nelements();
      76             :   }
      77             : 
      78             :   //(2)setup spline coefficients for each antenna ID
      79          21 :   timePointing.resize(nant);
      80          21 :   dirPointing.resize(nant);
      81          21 :   splineCoeff.resize(nant);
      82          21 :   doSplineInterpolation.resize(nant);
      83          21 :   doSplineInterpolation = false;
      84          42 :   for (Int i = 0; i < nant; ++i) {
      85          21 :     if (nPointingData(i) < 4) continue;
      86             :     
      87          21 :     doSplineInterpolation(i) = true;
      88          21 :     timePointing(i).resize(nPointingData(i));
      89          21 :     dirPointing(i).resize(nPointingData(i));
      90          21 :     splineCoeff(i).resize(nPointingData(i) - 1);
      91       20501 :     for (uInt j = 0; j < dirPointing(i).nelements(); ++j) {
      92       20480 :       dirPointing(i)(j).resize(2);
      93             :     }
      94       20480 :     for (uInt j = 0; j < splineCoeff(i).nelements(); ++j) {
      95       20459 :       splineCoeff(i)(j).resize(2);
      96       20459 :       splineCoeff(i)(j)(0).resize(4); // x
      97       20459 :       splineCoeff(i)(j)(1).resize(4); // y
      98             :     }
      99             :     
     100          21 :     Int npoi = nPointingData(i);
     101       20501 :     for (Int j = 0; j < npoi; ++j) {
     102       20480 :       timePointing(i)(j) = time(i)(j);
     103       61440 :       for (Int k = 0; k < 2; ++k) {
     104       40960 :         dirPointing(i)(j)(k) = dir(i)(j)(k);
     105             :       }
     106             :     }
     107             :       
     108          21 :     calcSplineCoeff(timePointing(i), dirPointing(i), splineCoeff(i));
     109             :   }
     110             : 
     111             :   //(3) keep time range
     112          21 :   timeRangeStart.resize(nant);
     113          21 :   timeRangeEnd.resize(nant);
     114          42 :   for (Int iant = 0; iant < nant; ++iant) {
     115          21 :     timeRangeStart(iant) = timePointing(iant)(0);
     116          21 :     timeRangeEnd(iant)   = timePointing(iant)(timePointing(iant).nelements()-1);
     117             :   }
     118          21 : }
     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          21 : void SDPosInterpolator::calcSplineCoeff(const Vector<Double>& time,
     215             :                                         const Vector<Vector<Double> >& dir,
     216             :                                         Vector<Vector<Vector<Double> > >& coeff) {
     217          21 :   Vector<Double> h, vx, vy;
     218          21 :   Vector<Double> a;
     219          21 :   Vector<Double> c;
     220          21 :   Vector<Double> alpha, beta, gamma;
     221          21 :   Vector<Double> wx, wy;
     222          21 :   Vector<Double> ux, uy;
     223             : 
     224          21 :   Int const num_data = time.nelements();
     225          21 :   h.resize(num_data-1);
     226          21 :   vx.resize(num_data-1);
     227          21 :   vy.resize(num_data-1);
     228          21 :   a.resize(num_data-1);
     229          21 :   c.resize(num_data-1);
     230          21 :   alpha.resize(num_data-1);
     231          21 :   beta.resize(num_data-1);
     232          21 :   gamma.resize(num_data-1);
     233          21 :   wx.resize(num_data-1);
     234          21 :   wy.resize(num_data-1);
     235          21 :   ux.resize(num_data);
     236          21 :   uy.resize(num_data);
     237             : 
     238          21 :   h(0) = time(1) - time(0);
     239       20459 :   for (Int i = 1; i < num_data-1; ++i) {
     240       20438 :     h(i) = time(i+1) - time(i);
     241       20438 :     vx(i) = 6.0*((dir(i+1)(0)-dir(i)(0))/h(i) - (dir(i)(0)-dir(i-1)(0))/h(i-1));
     242       20438 :     vy(i) = 6.0*((dir(i+1)(1)-dir(i)(1))/h(i) - (dir(i)(1)-dir(i-1)(1))/h(i-1));
     243       20438 :     a(i) = 2.0*(time(i+1) - time(i-1));
     244       20438 :     c(i) = h(i);
     245       20438 :     gamma(i) = c(i);
     246             :   }
     247          21 :   alpha(2) = c(1)/a(1);
     248       20417 :   for (Int i = 3; i < num_data-1; ++i) {
     249       20396 :     alpha(i) = c(i-1)/(a(i-1) - alpha(i-1)*c(i-2));
     250             :   }
     251          21 :   beta(1) = a(1);
     252       20417 :   for (Int i = 2; i < num_data-2; ++i) {
     253       20396 :     beta(i) = c(i)/alpha(i+1);
     254             :   }
     255          21 :   beta(num_data-2) = a(num_data-2) - alpha(num_data-2) * c(num_data-3);
     256          21 :   wx(0) = 0.0;
     257          21 :   wx(1) = vx(1);
     258          21 :   wy(0) = 0.0;
     259          21 :   wy(1) = vy(1);
     260       20438 :   for (Int i = 2; i < num_data-1; ++i) {
     261       20417 :     wx(i) = vx(i) - alpha(i)*wx(i-1);
     262       20417 :     wy(i) = vy(i) - alpha(i)*wy(i-1);
     263             :   }
     264          21 :   ux(num_data-1) = 0.0;
     265          21 :   uy(num_data-1) = 0.0;
     266       20459 :   for (Int i = num_data-2; i >= 1; --i) {
     267       20438 :     ux(i) = (wx(i) - gamma(i)*ux(i+1))/beta(i);
     268       20438 :     uy(i) = (wy(i) - gamma(i)*uy(i+1))/beta(i);
     269             :   }
     270          21 :   ux(0) = 0.0;
     271          21 :   uy(0) = 0.0;
     272             : 
     273       20480 :   for (Int i = 0; i < num_data-1; ++i) {
     274       20459 :     coeff(i)(0)(0) = dir(i)(0);
     275       20459 :     coeff(i)(1)(0) = dir(i)(1);
     276       20459 :     const auto dt = time(i+1)-time(i);
     277       20459 :     coeff(i)(0)(1) = (dir(i+1)(0)-dir(i)(0))/dt - dt*(2.0*ux(i)+ux(i+1))/6.0;
     278       20459 :     coeff(i)(1)(1) = (dir(i+1)(1)-dir(i)(1))/dt - dt*(2.0*uy(i)+uy(i+1))/6.0;
     279       20459 :     coeff(i)(0)(2) = ux(i)/2.0;
     280       20459 :     coeff(i)(1)(2) = uy(i)/2.0;
     281       20459 :     coeff(i)(0)(3) = (ux(i+1)-ux(i))/dt/6.0;
     282       20459 :     coeff(i)(1)(3) = (uy(i+1)-uy(i))/dt/6.0;
     283             :   }
     284          21 : }
     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          21 : Vector<Vector<Vector<Vector<Double> > > > SDPosInterpolator::getSplineCoeff() {
     323          21 :   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