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         209 : SDPosInterpolator::SDPosInterpolator(
      55             :   const MSPointingColumns& pointingColumns,
      56             :   const String& columnName,
      57         209 :   const size_t nant){
      58         209 :   setup(pointingColumns, columnName, nant);
      59         209 : }
      60         254 : SDPosInterpolator::SDPosInterpolator(
      61             :   const Vector<Vector<Double> >& time,
      62         254 :   const Vector<Vector<Vector<Double> > >& dir) {
      63         254 :   setup(time, dir);
      64         254 : }
      65         463 : SDPosInterpolator::~SDPosInterpolator() {}
      66             : 
      67         254 : 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         254 :   Int nant = time.nelements();
      72         254 :   Vector<uInt> nPointingData(nant);
      73         254 :   nPointingData = 0;
      74         538 :   for (Int iant = 0; iant < nant; ++iant) {
      75         284 :     nPointingData(iant) = time(iant).nelements();
      76             :   }
      77             : 
      78             :   //(2)setup spline coefficients for each antenna ID
      79         254 :   timePointing.resize(nant);
      80         254 :   dirPointing.resize(nant);
      81         254 :   splineCoeff.resize(nant);
      82         254 :   doSplineInterpolation.resize(nant);
      83         254 :   doSplineInterpolation = false;
      84         538 :   for (Int i = 0; i < nant; ++i) {
      85         284 :     if (nPointingData(i) < 4) continue;
      86             :     
      87         284 :     doSplineInterpolation(i) = true;
      88         284 :     timePointing(i).resize(nPointingData(i));
      89         284 :     dirPointing(i).resize(nPointingData(i));
      90         284 :     splineCoeff(i).resize(nPointingData(i) - 1);
      91     1111753 :     for (uInt j = 0; j < dirPointing(i).nelements(); ++j) {
      92     1111469 :       dirPointing(i)(j).resize(2);
      93             :     }
      94     1111469 :     for (uInt j = 0; j < splineCoeff(i).nelements(); ++j) {
      95     1111185 :       splineCoeff(i)(j).resize(2);
      96     1111185 :       splineCoeff(i)(j)(0).resize(4); // x
      97     1111185 :       splineCoeff(i)(j)(1).resize(4); // y
      98             :     }
      99             :     
     100         284 :     Int npoi = nPointingData(i);
     101     1111753 :     for (Int j = 0; j < npoi; ++j) {
     102     1111469 :       timePointing(i)(j) = time(i)(j);
     103     3334407 :       for (Int k = 0; k < 2; ++k) {
     104     2222938 :         dirPointing(i)(j)(k) = dir(i)(j)(k);
     105             :       }
     106             :     }
     107             :       
     108         284 :     calcSplineCoeff(timePointing(i), dirPointing(i), splineCoeff(i));
     109             :   }
     110             : 
     111             :   //(3) keep time range
     112         254 :   timeRangeStart.resize(nant);
     113         254 :   timeRangeEnd.resize(nant);
     114         538 :   for (Int iant = 0; iant < nant; ++iant) {
     115         284 :     timeRangeStart(iant) = timePointing(iant)(0);
     116         284 :     timeRangeEnd(iant)   = timePointing(iant)(timePointing(iant).nelements()-1);
     117             :   }
     118         254 : }
     119             : 
     120         209 : 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         209 :   };
     129         209 :   std::function<Vector<Double>(Int)> get_direction;
     130             : 
     131             :   //(0)check POINTING table and set function to obtain direction data
     132         209 :   if (pointingDirCol_p == "TARGET") {
     133           0 :     get_direction = [&](Int idx){
     134           0 :       return act_mspc.targetMeas(idx).getAngle("rad").getValue();
     135           0 :     };
     136         209 :   } 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         209 :   } 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         209 :   } 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       25058 :     get_direction = [&](Int idx){
     153       50116 :       return act_mspc.directionMeas(idx).getAngle("rad").getValue();
     154         209 :     };
     155             :   }
     156             : 
     157             :   //(1)get number of pointing data for each antennaID
     158         209 :   Vector<uInt> nPointingData(nant, 0);
     159         209 :   auto pointingRows = static_cast<size_t>(act_mspc.nrow());
     160       25267 :   for (size_t row = 0; row < pointingRows ; ++row) {
     161       25058 :     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         209 :   timePointing.resize(nant);
     172         209 :   dirPointing.resize(nant);
     173         209 :   splineCoeff.resize(nant);
     174         209 :   doSplineInterpolation.resize(nant);
     175         209 :   doSplineInterpolation = false;
     176         814 :   for (uInt i = 0; i < nant; ++i) {
     177         605 :     if (nPointingData(i) < 4) continue;
     178             : 
     179         605 :     doSplineInterpolation(i) = true;
     180         605 :     timePointing(i).resize(nPointingData(i));
     181         605 :     dirPointing(i).resize(nPointingData(i));
     182         605 :     splineCoeff(i).resize(nPointingData(i) - 1);
     183       25663 :     for (uInt j = 0; j < dirPointing(i).nelements(); ++j) {
     184       25058 :       dirPointing(i)(j).resize(2);
     185             :     }
     186       25058 :     for (uInt j = 0; j < splineCoeff(i).nelements(); ++j) {
     187       24453 :       splineCoeff(i)(j).resize(2);
     188       24453 :       splineCoeff(i)(j)(0).resize(4); // x
     189       24453 :       splineCoeff(i)(j)(1).resize(4); // y
     190             :     }
     191             : 
     192             :     //set ptime array etc. need for spline calculation...
     193         605 :     size_t tidx = 0;
     194       75559 :     for (size_t j = 0; j < pointingRows; ++j) {
     195       74954 :       if (act_mspc.antennaId()(j) != i) continue;
     196             : 
     197       25058 :       timePointing(i)(tidx) = act_mspc.time()(j);
     198       25058 :       dirPointing(i)(tidx) = get_direction(j);
     199       25058 :       tidx++;
     200             :     }
     201             : 
     202         605 :     calcSplineCoeff(timePointing(i), dirPointing(i), splineCoeff(i));
     203             :   }
     204             : 
     205             :   //(3) keep time range
     206         209 :   timeRangeStart.resize(nant);
     207         209 :   timeRangeEnd.resize(nant);
     208         814 :   for (size_t iant = 0; iant < nant; ++iant) {
     209         605 :     timeRangeStart(iant) = timePointing(iant)(0);
     210         605 :     timeRangeEnd(iant)   = timePointing(iant)(timePointing(iant).nelements()-1);
     211             :   }
     212         209 : }
     213             : 
     214         889 : void SDPosInterpolator::calcSplineCoeff(const Vector<Double>& time,
     215             :                                         const Vector<Vector<Double> >& dir,
     216             :                                         Vector<Vector<Vector<Double> > >& coeff) {
     217         889 :   Vector<Double> h, vx, vy;
     218         889 :   Vector<Double> a;
     219         889 :   Vector<Double> c;
     220         889 :   Vector<Double> alpha, beta, gamma;
     221         889 :   Vector<Double> wx, wy;
     222         889 :   Vector<Double> ux, uy;
     223             : 
     224         889 :   Int const num_data = time.nelements();
     225         889 :   h.resize(num_data-1);
     226         889 :   vx.resize(num_data-1);
     227         889 :   vy.resize(num_data-1);
     228         889 :   a.resize(num_data-1);
     229         889 :   c.resize(num_data-1);
     230         889 :   alpha.resize(num_data-1);
     231         889 :   beta.resize(num_data-1);
     232         889 :   gamma.resize(num_data-1);
     233         889 :   wx.resize(num_data-1);
     234         889 :   wy.resize(num_data-1);
     235         889 :   ux.resize(num_data);
     236         889 :   uy.resize(num_data);
     237             : 
     238         889 :   h(0) = time(1) - time(0);
     239     1135638 :   for (Int i = 1; i < num_data-1; ++i) {
     240     1134749 :     h(i) = time(i+1) - time(i);
     241     1134749 :     vx(i) = 6.0*((dir(i+1)(0)-dir(i)(0))/h(i) - (dir(i)(0)-dir(i-1)(0))/h(i-1));
     242     1134749 :     vy(i) = 6.0*((dir(i+1)(1)-dir(i)(1))/h(i) - (dir(i)(1)-dir(i-1)(1))/h(i-1));
     243     1134749 :     a(i) = 2.0*(time(i+1) - time(i-1));
     244     1134749 :     c(i) = h(i);
     245     1134749 :     gamma(i) = c(i);
     246             :   }
     247         889 :   alpha(2) = c(1)/a(1);
     248     1133860 :   for (Int i = 3; i < num_data-1; ++i) {
     249     1132971 :     alpha(i) = c(i-1)/(a(i-1) - alpha(i-1)*c(i-2));
     250             :   }
     251         889 :   beta(1) = a(1);
     252     1133860 :   for (Int i = 2; i < num_data-2; ++i) {
     253     1132971 :     beta(i) = c(i)/alpha(i+1);
     254             :   }
     255         889 :   beta(num_data-2) = a(num_data-2) - alpha(num_data-2) * c(num_data-3);
     256         889 :   wx(0) = 0.0;
     257         889 :   wx(1) = vx(1);
     258         889 :   wy(0) = 0.0;
     259         889 :   wy(1) = vy(1);
     260     1134749 :   for (Int i = 2; i < num_data-1; ++i) {
     261     1133860 :     wx(i) = vx(i) - alpha(i)*wx(i-1);
     262     1133860 :     wy(i) = vy(i) - alpha(i)*wy(i-1);
     263             :   }
     264         889 :   ux(num_data-1) = 0.0;
     265         889 :   uy(num_data-1) = 0.0;
     266     1135638 :   for (Int i = num_data-2; i >= 1; --i) {
     267     1134749 :     ux(i) = (wx(i) - gamma(i)*ux(i+1))/beta(i);
     268     1134749 :     uy(i) = (wy(i) - gamma(i)*uy(i+1))/beta(i);
     269             :   }
     270         889 :   ux(0) = 0.0;
     271         889 :   uy(0) = 0.0;
     272             : 
     273     1136527 :   for (Int i = 0; i < num_data-1; ++i) {
     274     1135638 :     coeff(i)(0)(0) = dir(i)(0);
     275     1135638 :     coeff(i)(1)(0) = dir(i)(1);
     276     1135638 :     const auto dt = time(i+1)-time(i);
     277     1135638 :     coeff(i)(0)(1) = (dir(i+1)(0)-dir(i)(0))/dt - dt*(2.0*ux(i)+ux(i+1))/6.0;
     278     1135638 :     coeff(i)(1)(1) = (dir(i+1)(1)-dir(i)(1))/dt - dt*(2.0*uy(i)+uy(i+1))/6.0;
     279     1135638 :     coeff(i)(0)(2) = ux(i)/2.0;
     280     1135638 :     coeff(i)(1)(2) = uy(i)/2.0;
     281     1135638 :     coeff(i)(0)(3) = (ux(i+1)-ux(i))/dt/6.0;
     282     1135638 :     coeff(i)(1)(3) = (uy(i+1)-uy(i))/dt/6.0;
     283             :   }
     284         889 : }
     285             : 
     286       24621 : MDirection SDPosInterpolator::interpolateDirectionMeasSpline(
     287             :             const MSPointingColumns& mspc,
     288             :             const Double& time,
     289             :             const Int& index,
     290             :             const Int& antid) {
     291       24621 :   Int lastIndex = timePointing(antid).nelements() - 1;
     292       24621 :   Int aindex = lastIndex;
     293      351601 :   for (uInt i = 0; i < timePointing(antid).nelements(); ++i) {
     294      351529 :     if (time < timePointing(antid)(i)) {
     295       24549 :       aindex = i-1;
     296       24549 :       break;
     297             :     }
     298             :   }
     299       24621 :   if (aindex < 0) aindex = 0;
     300       24621 :   if (lastIndex <= aindex) aindex = lastIndex - 1;
     301             : 
     302       24621 :   auto const &coeff = splineCoeff(antid)(aindex);
     303       24621 :   Double dt = time - timePointing(antid)(aindex);
     304       24621 :   Vector<Double> newdir(2);
     305             :   // Why don't we use Horner's method here ?
     306       49242 :   newdir(0) =
     307       24621 :     coeff(0)(0) + coeff(0)(1)*dt + coeff(0)(2)*dt*dt + coeff(0)(3)*dt*dt*dt;
     308       49242 :   newdir(1) =
     309       24621 :     coeff(1)(0) + coeff(1)(1)*dt + coeff(1)(2)*dt*dt + coeff(1)(3)*dt*dt*dt;
     310             : 
     311       24621 :   Quantity rDirLon(newdir(0), "rad");
     312       24621 :   Quantity rDirLat(newdir(1), "rad");
     313       24621 :   auto const &directionMeasColumn = mspc.directionMeasCol();
     314       24621 :   MDirection::Ref rf(directionMeasColumn.measDesc().getRefCode());
     315       24621 :   if (directionMeasColumn.measDesc().isRefCodeVariable()) {
     316           0 :     rf = mspc.directionMeas(index).getRef();
     317             :   }
     318             : 
     319       49242 :   return MDirection(rDirLon, rDirLat, rf);
     320       24621 : }
     321             : 
     322         254 : Vector<Vector<Vector<Vector<Double> > > > SDPosInterpolator::getSplineCoeff() {
     323         254 :   return splineCoeff;
     324             : }
     325             : 
     326       24583 : Bool SDPosInterpolator::inTimeRange(const Double& time, const Int& antid) {
     327       24583 :   Bool inrange = false;
     328       24583 :   if ((timeRangeStart(antid) <= time) && (time <= timeRangeEnd(antid))) {
     329       24393 :     inrange = true;
     330             :   }
     331       24583 :   return inrange;
     332             : }
     333             : 
     334             : } //#End casa namespace