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