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: aips2-request@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 204 : SDPosInterpolator::SDPosInterpolator(
41 : const vi::VisBuffer2& vb,
42 204 : const String& pointingDirCol_p) {
43 204 : const auto & pointingColumns = vb.subtableColumns().pointing();
44 204 : const auto nant = static_cast<size_t>(vb.subtableColumns().antenna().nrow());
45 204 : setup(pointingColumns, pointingDirCol_p, nant);
46 204 : }
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 5 : SDPosInterpolator::SDPosInterpolator(
55 : const MSPointingColumns& pointingColumns,
56 : const String& columnName,
57 5 : const size_t nant){
58 5 : setup(pointingColumns, columnName, nant);
59 5 : }
60 254 : SDPosInterpolator::SDPosInterpolator(const Vector<Vector<Double> >& time,
61 254 : const Vector<Vector<Vector<Double> > >& dir) {
62 254 : setup(time, dir);
63 254 : }
64 463 : SDPosInterpolator::~SDPosInterpolator() {}
65 :
66 254 : void SDPosInterpolator::setup(const Vector<Vector<Double> >& time,
67 : const Vector<Vector<Vector<Double> > >& dir) {
68 : //(1)get number of pointing data for each antennaID
69 254 : Int nant = time.nelements();
70 508 : Vector<uInt> nPointingData(nant);
71 254 : nPointingData = 0;
72 538 : for (Int iant = 0; iant < nant; ++iant) {
73 284 : nPointingData(iant) = time(iant).nelements();
74 : }
75 :
76 : //(2)setup spline coefficients for each antenna ID
77 254 : timePointing.resize(nant);
78 254 : dirPointing.resize(nant);
79 254 : splineCoeff.resize(nant);
80 254 : doSplineInterpolation.resize(nant);
81 254 : doSplineInterpolation = false;
82 538 : for (Int i = 0; i < nant; ++i) {
83 284 : if (nPointingData(i) < 4) continue;
84 :
85 284 : doSplineInterpolation(i) = true;
86 284 : timePointing(i).resize(nPointingData(i));
87 284 : dirPointing(i).resize(nPointingData(i));
88 284 : splineCoeff(i).resize(nPointingData(i) - 1);
89 1111753 : for (uInt j = 0; j < dirPointing(i).nelements(); ++j) {
90 1111469 : dirPointing(i)(j).resize(2);
91 : }
92 1111469 : for (uInt j = 0; j < splineCoeff(i).nelements(); ++j) {
93 1111185 : splineCoeff(i)(j).resize(2);
94 1111185 : splineCoeff(i)(j)(0).resize(4); // x
95 1111185 : splineCoeff(i)(j)(1).resize(4); // y
96 : }
97 :
98 284 : Int npoi = nPointingData(i);
99 1111753 : for (Int j = 0; j < npoi; ++j) {
100 1111469 : timePointing(i)(j) = time(i)(j);
101 3334407 : for (Int k = 0; k < 2; ++k) {
102 2222938 : dirPointing(i)(j)(k) = dir(i)(j)(k);
103 : }
104 : }
105 :
106 284 : calcSplineCoeff(timePointing(i), dirPointing(i), splineCoeff(i));
107 : }
108 :
109 : //(3) keep time range
110 254 : timeRangeStart.resize(nant);
111 254 : timeRangeEnd.resize(nant);
112 538 : for (Int iant = 0; iant < nant; ++iant) {
113 284 : timeRangeStart(iant) = timePointing(iant)(0);
114 284 : timeRangeEnd(iant) = timePointing(iant)(timePointing(iant).nelements()-1);
115 : }
116 254 : }
117 :
118 209 : void SDPosInterpolator::setup(
119 : const MSPointingColumns& act_mspc,
120 : const String& pointingDirCol_p,
121 : size_t nant) {
122 0 : auto check_col = [&](Bool isnull){
123 0 : if (isnull) {
124 0 : cerr << "No " << pointingDirCol_p << " column in POINTING table" << endl;
125 : }
126 209 : };
127 418 : std::function<Vector<Double>(Int)> get_direction;
128 :
129 : //(0)check POINTING table and set function to obtain direction data
130 209 : if (pointingDirCol_p == "TARGET") {
131 0 : get_direction = [&](Int idx){
132 0 : return act_mspc.targetMeas(idx).getAngle("rad").getValue();
133 0 : };
134 209 : } else if (pointingDirCol_p == "POINTING_OFFSET") {
135 0 : check_col(act_mspc.pointingOffsetMeasCol().isNull());
136 0 : get_direction = [&](Int idx){
137 0 : return act_mspc.pointingOffsetMeas(idx).getAngle("rad").getValue();
138 0 : };
139 209 : } else if (pointingDirCol_p == "SOURCE_OFFSET") {
140 0 : check_col(act_mspc.sourceOffsetMeasCol().isNull());
141 0 : get_direction = [&](Int idx){
142 0 : return act_mspc.sourceOffsetMeas(idx).getAngle("rad").getValue();
143 0 : };
144 209 : } else if (pointingDirCol_p == "ENCODER") {
145 0 : check_col(act_mspc.encoderMeas().isNull());
146 0 : get_direction = [&](Int idx){
147 0 : return act_mspc.encoderMeas()(idx).getAngle("rad").getValue();
148 0 : };
149 : } else {
150 25058 : get_direction = [&](Int idx){
151 25058 : return act_mspc.directionMeas(idx).getAngle("rad").getValue();
152 209 : };
153 : }
154 :
155 : //(1)get number of pointing data for each antennaID
156 418 : Vector<uInt> nPointingData(nant, 0);
157 209 : auto pointingRows = static_cast<size_t>(act_mspc.nrow());
158 25267 : for (size_t row = 0; row < pointingRows ; ++row) {
159 25058 : nPointingData(act_mspc.antennaId()(row)) += 1;
160 : }
161 :
162 : //(2)setup spline coefficients for each antenna ID that
163 : // appear in the main table (spectral data) if there
164 : // are enough number of pointing data (4 or more).
165 : // in case there exists antenna ID for which not enough
166 : // (i.e., 1, 2 or 3) pointing data are given, linear
167 : // interpolation is applied for that antenna ID as
168 : // previously done.
169 209 : timePointing.resize(nant);
170 209 : dirPointing.resize(nant);
171 209 : splineCoeff.resize(nant);
172 209 : doSplineInterpolation.resize(nant);
173 209 : doSplineInterpolation = false;
174 814 : for (uInt i = 0; i < nant; ++i) {
175 605 : if (nPointingData(i) < 4) continue;
176 :
177 605 : doSplineInterpolation(i) = true;
178 605 : timePointing(i).resize(nPointingData(i));
179 605 : dirPointing(i).resize(nPointingData(i));
180 605 : splineCoeff(i).resize(nPointingData(i) - 1);
181 25663 : for (uInt j = 0; j < dirPointing(i).nelements(); ++j) {
182 25058 : dirPointing(i)(j).resize(2);
183 : }
184 25058 : for (uInt j = 0; j < splineCoeff(i).nelements(); ++j) {
185 24453 : splineCoeff(i)(j).resize(2);
186 24453 : splineCoeff(i)(j)(0).resize(4); // x
187 24453 : splineCoeff(i)(j)(1).resize(4); // y
188 : }
189 :
190 : //set ptime array etc. need for spline calculation...
191 605 : size_t tidx = 0;
192 75559 : for (size_t j = 0; j < pointingRows; ++j) {
193 74954 : if (act_mspc.antennaId()(j) != i) continue;
194 :
195 25058 : timePointing(i)(tidx) = act_mspc.time()(j);
196 25058 : dirPointing(i)(tidx) = get_direction(j);
197 25058 : tidx++;
198 : }
199 :
200 605 : calcSplineCoeff(timePointing(i), dirPointing(i), splineCoeff(i));
201 : }
202 :
203 : //(3) keep time range
204 209 : timeRangeStart.resize(nant);
205 209 : timeRangeEnd.resize(nant);
206 814 : for (size_t iant = 0; iant < nant; ++iant) {
207 605 : timeRangeStart(iant) = timePointing(iant)(0);
208 605 : timeRangeEnd(iant) = timePointing(iant)(timePointing(iant).nelements()-1);
209 : }
210 209 : }
211 :
212 889 : void SDPosInterpolator::calcSplineCoeff(const Vector<Double>& time,
213 : const Vector<Vector<Double> >& dir,
214 : Vector<Vector<Vector<Double> > >& coeff) {
215 1778 : Vector<Double> h, vx, vy;
216 1778 : Vector<Double> a;
217 1778 : Vector<Double> c;
218 1778 : Vector<Double> alpha, beta, gamma;
219 1778 : Vector<Double> wx, wy;
220 1778 : Vector<Double> ux, uy;
221 :
222 889 : Int const num_data = time.nelements();
223 889 : h.resize(num_data-1);
224 889 : vx.resize(num_data-1);
225 889 : vy.resize(num_data-1);
226 889 : a.resize(num_data-1);
227 889 : c.resize(num_data-1);
228 889 : alpha.resize(num_data-1);
229 889 : beta.resize(num_data-1);
230 889 : gamma.resize(num_data-1);
231 889 : wx.resize(num_data-1);
232 889 : wy.resize(num_data-1);
233 889 : ux.resize(num_data);
234 889 : uy.resize(num_data);
235 :
236 889 : h(0) = time(1) - time(0);
237 1135638 : for (Int i = 1; i < num_data-1; ++i) {
238 1134749 : h(i) = time(i+1) - time(i);
239 1134749 : vx(i) = 6.0*((dir(i+1)(0)-dir(i)(0))/h(i) - (dir(i)(0)-dir(i-1)(0))/h(i-1));
240 1134749 : vy(i) = 6.0*((dir(i+1)(1)-dir(i)(1))/h(i) - (dir(i)(1)-dir(i-1)(1))/h(i-1));
241 1134749 : a(i) = 2.0*(time(i+1) - time(i-1));
242 1134749 : c(i) = h(i);
243 1134749 : gamma(i) = c(i);
244 : }
245 889 : alpha(2) = c(1)/a(1);
246 1133860 : for (Int i = 3; i < num_data-1; ++i) {
247 1132971 : alpha(i) = c(i-1)/(a(i-1) - alpha(i-1)*c(i-2));
248 : }
249 889 : beta(1) = a(1);
250 1133860 : for (Int i = 2; i < num_data-2; ++i) {
251 1132971 : beta(i) = c(i)/alpha(i+1);
252 : }
253 889 : beta(num_data-2) = a(num_data-2) - alpha(num_data-2) * c(num_data-3);
254 889 : wx(0) = 0.0;
255 889 : wx(1) = vx(1);
256 889 : wy(0) = 0.0;
257 889 : wy(1) = vy(1);
258 1134749 : for (Int i = 2; i < num_data-1; ++i) {
259 1133860 : wx(i) = vx(i) - alpha(i)*wx(i-1);
260 1133860 : wy(i) = vy(i) - alpha(i)*wy(i-1);
261 : }
262 889 : ux(num_data-1) = 0.0;
263 889 : uy(num_data-1) = 0.0;
264 1135638 : for (Int i = num_data-2; i >= 1; --i) {
265 1134749 : ux(i) = (wx(i) - gamma(i)*ux(i+1))/beta(i);
266 1134749 : uy(i) = (wy(i) - gamma(i)*uy(i+1))/beta(i);
267 : }
268 889 : ux(0) = 0.0;
269 889 : uy(0) = 0.0;
270 :
271 1136527 : for (Int i = 0; i < num_data-1; ++i) {
272 1135638 : coeff(i)(0)(0) = dir(i)(0);
273 1135638 : coeff(i)(1)(0) = dir(i)(1);
274 1135638 : coeff(i)(0)(1) = (dir(i+1)(0)-dir(i)(0))/(time(i+1)-time(i)) - (time(i+1)-time(i))*(2.0*ux(i)+ux(i+1))/6.0;
275 1135638 : coeff(i)(1)(1) = (dir(i+1)(1)-dir(i)(1))/(time(i+1)-time(i)) - (time(i+1)-time(i))*(2.0*uy(i)+uy(i+1))/6.0;
276 1135638 : coeff(i)(0)(2) = ux(i)/2.0;
277 1135638 : coeff(i)(1)(2) = uy(i)/2.0;
278 1135638 : coeff(i)(0)(3) = (ux(i+1)-ux(i))/(time(i+1)-time(i))/6.0;
279 1135638 : coeff(i)(1)(3) = (uy(i+1)-uy(i))/(time(i+1)-time(i))/6.0;
280 : }
281 889 : }
282 :
283 24621 : MDirection SDPosInterpolator::interpolateDirectionMeasSpline(const MSPointingColumns& mspc,
284 : const Double& time,
285 : const Int& index,
286 : const Int& antid) {
287 24621 : Int lastIndex = timePointing(antid).nelements() - 1;
288 24621 : Int aindex = lastIndex;
289 351601 : for (uInt i = 0; i < timePointing(antid).nelements(); ++i) {
290 351529 : if (time < timePointing(antid)(i)) {
291 24549 : aindex = i-1;
292 24549 : break;
293 : }
294 : }
295 24621 : if (aindex < 0) aindex = 0;
296 24621 : if (lastIndex <= aindex) aindex = lastIndex - 1;
297 :
298 24621 : auto const &coeff = splineCoeff(antid)(aindex);
299 24621 : Double dt = time - timePointing(antid)(aindex);
300 49242 : Vector<Double> newdir(2);
301 24621 : newdir(0) = coeff(0)(0) + coeff(0)(1)*dt + coeff(0)(2)*dt*dt + coeff(0)(3)*dt*dt*dt;
302 24621 : newdir(1) = coeff(1)(0) + coeff(1)(1)*dt + coeff(1)(2)*dt*dt + coeff(1)(3)*dt*dt*dt;
303 :
304 49242 : Quantity rDirLon(newdir(0), "rad");
305 49242 : Quantity rDirLat(newdir(1), "rad");
306 24621 : auto const &directionMeasColumn = mspc.directionMeasCol();
307 49242 : MDirection::Ref rf(directionMeasColumn.measDesc().getRefCode());
308 24621 : if (directionMeasColumn.measDesc().isRefCodeVariable()) {
309 0 : rf = mspc.directionMeas(index).getRef();
310 : }
311 :
312 49242 : return MDirection(rDirLon, rDirLat, rf);
313 : }
314 :
315 254 : Vector<Vector<Vector<Vector<Double> > > > SDPosInterpolator::getSplineCoeff() {
316 254 : return splineCoeff;
317 : }
318 :
319 24583 : Bool SDPosInterpolator::inTimeRange(const Double& time, const Int& antid) {
320 24583 : Bool inrange = false;
321 24583 : if ((timeRangeStart(antid) <= time) && (time <= timeRangeEnd(antid))) {
322 24393 : inrange = true;
323 : }
324 24583 : return inrange;
325 : }
326 :
327 : } //#End casa namespace
|
