Line data Source code
1 : //# PointingDirectionCalculator.cc: Implementation of PointingDirectionCalculator.h
2 : //# All helper functions of imager moved here for readability
3 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: casa-feedback@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 : #include <cassert>
29 : #include <sstream>
30 : #include <iostream>
31 : #include <iomanip>
32 :
33 : #include <synthesis/Utilities/PointingDirectionCalculator.h>
34 : #include <memory> // for unique_ptr<>
35 : #include <utility> // for std::pair
36 :
37 : #include <casacore/casa/aipstype.h>
38 : #include <casacore/casa/Exceptions/Error.h>
39 : #include <casacore/casa/IO/ArrayIO.h>
40 : #include <casacore/casa/Arrays/ArrayMath.h>
41 : #include <casacore/casa/BasicSL/String.h>
42 : #include <casacore/casa/Quanta/Quantum.h>
43 : #include <casacore/casa/Containers/Block.h>
44 : #include <casacore/casa/Utilities/BinarySearch.h>
45 : #include <casacore/casa/Logging/LogIO.h>
46 : #include <casacore/tables/TaQL/ExprNode.h>
47 : #include <casacore/ms/MSSel/MSSelection.h>
48 : #include <casacore/measures/Measures/MPosition.h>
49 : #include <casacore/measures/Measures/MEpoch.h>
50 : #include <casacore/measures/Measures/MDirection.h>
51 : // CAS-8418 NEW //
52 : #include <synthesis/Utilities/SDPosInterpolator.h>
53 :
54 : using namespace casacore;
55 : using namespace casacore;
56 : using namespace std;
57 :
58 : // Debug Message Handling
59 : // if DIRECTIONCALC_DEBUG is defined, the macro debuglog and
60 : // debugpost point standard output stream (std::cout and
61 : // std::endl so that debug messages are sent to standard
62 : // output. Otherwise, these macros basically does nothing.
63 : // "Do nothing" behavior is implemented in NullLogger
64 : // and its associating << operator below.
65 : //
66 : // Usage:
67 : // Similar to standard output stream.
68 : //
69 : // debuglog << "Any message" << any_value << debugpost;
70 : //
71 :
72 : // #define DIRECTIONCALC_DEBUG
73 :
74 : namespace {
75 : struct NullLogger {
76 : };
77 :
78 : template<class T>
79 22613264 : inline NullLogger &operator<<(NullLogger &logger, T /*value*/) {
80 22613264 : return logger;
81 : }
82 :
83 : #ifndef DIRECTIONCALC_DEBUG
84 : NullLogger nulllogger;
85 : #endif
86 : }
87 :
88 : #ifdef DIRECTIONCALC_DEBUG
89 : #define debuglog cout << "PointingDirectionCalculator::DEBUG "
90 : #define debugpost endl
91 : #else
92 : #define debuglog nulllogger
93 : #define debugpost 0
94 : #endif
95 :
96 : namespace {
97 : #define ARRAY_DIRECTION(ColumnName) \
98 : inline MDirection ColumnName ## Accessor(MSPointingColumns &pointingColumns, uInt rownr) { \
99 : return pointingColumns.ColumnName ## Meas(rownr); \
100 : }
101 :
102 : #define SCALAR_DIRECTION(ColumnName) \
103 : inline MDirection ColumnName ## Accessor(MSPointingColumns &pointingColumns, uInt rownr) { \
104 : return pointingColumns.ColumnName ## Meas()(rownr); \
105 : }
106 :
107 1581666 : ARRAY_DIRECTION(direction)
108 0 : ARRAY_DIRECTION(target)
109 0 : ARRAY_DIRECTION(pointingOffset)
110 0 : ARRAY_DIRECTION(sourceOffset)
111 0 : SCALAR_DIRECTION(encoder)
112 :
113 : // working function for moving source correction
114 : // convertToAzel must be configured with moving source direction and
115 : // proper reference frame. Also, convertToCelestial must refer proper
116 : // reference frame.
117 30445 : inline void performMovingSourceCorrection(
118 : CountedPtr<MDirection::Convert> &convertToAzel,
119 : CountedPtr<MDirection::Convert> &convertToCelestial,
120 : Vector<Double> &direction) {
121 : // moving source handling
122 : // If moving source is specified, output direction list is always
123 : // offset from reference position of moving source
124 :
125 30445 : debuglog << "MovingSourceCorrection <Working>." << debugpost;
126 :
127 : // DEBUG (CAS-11818) //
128 30445 : assert( convertToCelestial != nullptr );
129 30445 : assert( convertToAzel != nullptr );
130 :
131 30445 : MDirection srcAzel = (*convertToAzel)();
132 30445 : MDirection srcDirection = (*convertToCelestial)(srcAzel);
133 60890 : Vector<Double> srcDirectionVal = srcDirection.getAngle("rad").getValue();
134 30445 : direction -= srcDirectionVal;
135 30445 : }
136 :
137 436679 : inline void skipMovingSourceCorrection(
138 : CountedPtr<MDirection::Convert> &/*convertToAzel*/,
139 : CountedPtr<MDirection::Convert> &/*convertToCelestial*/,
140 : Vector<Double> &/*direction*/) {
141 :
142 436679 : debuglog << "MovingSourceCorrection <NO ACTION>" << debugpost;
143 :
144 : // do nothing
145 436679 : }
146 : } // anonymous namespace
147 :
148 : using namespace casacore;
149 : namespace casa {
150 :
151 : //+
152 : // Original Constructor
153 : //-
154 304 : PointingDirectionCalculator::PointingDirectionCalculator(
155 304 : MeasurementSet const &ms) :
156 304 : originalMS_(new MeasurementSet(ms)), selectedMS_(), pointingTable_(),
157 304 : pointingColumns_(), timeColumn_(), intervalColumn_(), antennaColumn_(),
158 304 : directionColumnName_(), accessor_(NULL), antennaPosition_(), referenceEpoch_(),
159 304 : referenceFrame_(referenceEpoch_, antennaPosition_),
160 304 : directionConvert_(NULL), directionType_(MDirection::J2000), movingSource_(NULL),
161 304 : movingSourceConvert_(NULL), movingSourceCorrection_(NULL),
162 304 : antennaBoundary_(), numAntennaBoundary_(0), pointingTimeUTC_(), lastTimeStamp_(-1.0),
163 304 : lastAntennaIndex_(-1), pointingTableIndexCache_(0),
164 304 : shape_(PointingDirectionCalculator::COLUMN_MAJOR),
165 :
166 304 : /*CAS-8418*/ useSplineInterpolation_(true), ////!! Set, when Spline is used. ////
167 304 : /*CAS-8418*/ currSpline_(nullptr),
168 304 : /*CAS-8418*/ splineObj_(PointingDirectionCalculator::PtColID::nItems),
169 304 : /*CAS-8418*/ initializeReady_(PointingDirectionCalculator::PtColID::nItems,false),// Spline initialization Ready
170 304 : /*CAS-8418*/ coefficientReady_(PointingDirectionCalculator::PtColID::nItems,false), // Spline Coefficient Ready
171 304 : /*CAS-8418*/ accessorId_(DIRECTION) // specify default accessor ID
172 : {
173 : // -- original code -- //
174 304 : accessor_ = directionAccessor;
175 :
176 304 : Block<String> sortColumns(2);
177 304 : sortColumns[0] = "ANTENNA1";
178 304 : sortColumns[1] = "TIME";
179 :
180 304 : selectedMS_ = new MeasurementSet(originalMS_->sort(sortColumns));
181 304 : debuglog << "Calling init()" << debugpost;
182 304 : init();
183 :
184 : // set default output direction reference frame
185 304 : debuglog << "Calling setFrame(J2000)" << debugpost;
186 304 : setFrame("J2000");
187 :
188 : // set default direction column name
189 304 : setDirectionColumn("DIRECTION");
190 304 : }
191 :
192 :
193 304 : PointingDirectionCalculator::~PointingDirectionCalculator()
194 : {
195 : // Nothing at the moment. //
196 304 : }
197 :
198 571 : void PointingDirectionCalculator::init() {
199 : // attach column
200 571 : timeColumn_.attach(*selectedMS_, "TIME");
201 571 : intervalColumn_.attach(*selectedMS_, "INTERVAL");
202 571 : antennaColumn_.attach(*selectedMS_, "ANTENNA1");
203 :
204 : // initial setup
205 571 : debuglog << "inspectAntenna" << debugpost;
206 571 : inspectAntenna();
207 571 : debuglog << "done" << debugpost;
208 :
209 571 : resetAntennaPosition(antennaColumn_(0));
210 571 : }
211 :
212 293 : void PointingDirectionCalculator::selectData(String const &antenna,
213 : String const &spw, String const &field, String const &time,
214 : String const &scan, String const &feed, String const &intent,
215 : String const &observation, String const &uvrange,
216 : String const &msselect) {
217 : // table selection
218 293 : MSSelection thisSelection;
219 293 : thisSelection.setAntennaExpr(antenna);
220 293 : thisSelection.setSpwExpr(spw);
221 293 : thisSelection.setFieldExpr(field);
222 293 : thisSelection.setTimeExpr(time);
223 293 : thisSelection.setScanExpr(scan);
224 293 : thisSelection.setStateExpr(intent);
225 293 : thisSelection.setObservationExpr(observation);
226 293 : thisSelection.setUvDistExpr(uvrange);
227 293 : thisSelection.setTaQLExpr(msselect);
228 :
229 293 : TableExprNode exprNode = thisSelection.getTEN(&(*originalMS_));
230 :
231 : // sort by ANTENNA1 and TIME for performance reason
232 269 : Block<String> sortColumns(2);
233 269 : sortColumns[0] = "ANTENNA1";
234 269 : sortColumns[1] = "TIME";
235 269 : if (exprNode.isNull()) {
236 0 : debuglog << "NULL selection" << debugpost;
237 0 : selectedMS_ = new MeasurementSet(originalMS_->sort(sortColumns));
238 : } else {
239 269 : debuglog << "Sort of selection" << debugpost;
240 269 : MeasurementSet tmp = (*originalMS_)(exprNode);
241 269 : selectedMS_ = new MeasurementSet(tmp.sort(sortColumns));
242 269 : }
243 269 : debuglog << "selectedMS_->nrow() = " << selectedMS_->nrow() << debugpost;
244 269 : if (selectedMS_->nrow() == 0) {
245 2 : stringstream ss;
246 2 : ss << "Selected MS is empty for given selection: " << endl;
247 2 : if (!antenna.empty()) {
248 2 : ss << "\tantenna \"" << antenna << "\"" << endl;
249 : }
250 2 : if (!spw.empty()) {
251 2 : ss << "\tspw \"" << spw << "\"" << endl;
252 : }
253 2 : if (!field.empty()) {
254 2 : ss << "\tfield \"" << field << "\"" << endl;
255 : }
256 2 : if (!time.empty()) {
257 0 : ss << "\ttime \"" << time << "\"" << endl;
258 : }
259 2 : if (!scan.empty()) {
260 0 : ss << "\tscan \"" << scan << "\"" << endl;
261 : }
262 2 : if (!feed.empty()) {
263 0 : ss << "\tfeed \"" << feed << "\"" << endl;
264 : }
265 2 : if (!intent.empty()) {
266 0 : ss << "\tintent \"" << intent << "\"" << endl;
267 : }
268 2 : if (!observation.empty()) {
269 0 : ss << "\tobservation \"" << observation << "\"" << endl;
270 : }
271 2 : if (!uvrange.empty()) {
272 0 : ss << "\tuvrange \"" << uvrange << "\"" << endl;
273 : }
274 2 : if (!msselect.empty()) {
275 0 : ss << "\tmsselect \"" << msselect << "\"" << endl;
276 : }
277 :
278 2 : throw AipsError(ss.str());
279 2 : }
280 :
281 267 : init();
282 :
283 267 : debuglog << "done selectdata" << debugpost;
284 297 : }
285 :
286 608 : void PointingDirectionCalculator::configureMovingSourceCorrection() {
287 : #if 0 // Causes Crash //
288 : if (!movingSource_.null() || directionColumnName_.contains("OFFSET")) {
289 : #else // FIXED //
290 608 : if ( !movingSource_.null() && !directionColumnName_.contains("OFFSET") )
291 : {
292 : #endif
293 9 : debuglog << "configureMovingSourceCorrection::Perfrom." << debugpost;
294 9 : movingSourceCorrection_ = performMovingSourceCorrection;
295 : } else {
296 599 : debuglog << "configureMovingSourceCorrection::Skip." << debugpost;
297 599 : movingSourceCorrection_ = skipMovingSourceCorrection;
298 : }
299 608 : }
300 :
301 :
302 599 : void PointingDirectionCalculator::setDirectionColumn(String const &columnName) {
303 599 : String columnNameUpcase = columnName;
304 599 : columnNameUpcase.upcase();
305 599 : if (!(originalMS_->pointing().tableDesc().isColumn(columnNameUpcase))) {
306 0 : stringstream ss;
307 0 : ss << "Column \"" << columnNameUpcase
308 0 : << "\" doesn't exist in POINTING table.";
309 0 : throw AipsError(ss.str());
310 0 : }
311 :
312 599 : directionColumnName_ = columnNameUpcase;
313 :
314 : //+
315 : // New code reuired by CAS-8418
316 : // - When setDirectionColumn is called ,
317 : // - new spline coeffcient table corresoinding to specified Direction column is generated.
318 : // Once generated and from next time, lately created Obj. is pointed and used.
319 : // - from init(), this function is called.
320 : // - See indetal in;
321 : // initializeSplinefromPointingColumn(*originalMS_, accessorId_ );
322 : //-
323 :
324 599 : if (directionColumnName_ == "DIRECTION") {
325 599 : accessor_ = directionAccessor;
326 599 : accessorId_ = DIRECTION;
327 599 : initializeSplinefromPointingColumn(*originalMS_, accessorId_ );
328 :
329 0 : } else if (directionColumnName_ == "TARGET") {
330 0 : accessor_ = targetAccessor;
331 0 : accessorId_ = TARGET;
332 0 : initializeSplinefromPointingColumn(*originalMS_, accessorId_ );
333 :
334 0 : } else if (directionColumnName_ == "POINTING_OFFSET") {
335 0 : accessor_ = pointingOffsetAccessor;
336 0 : accessorId_ = POINTING_OFFSET;
337 0 : initializeSplinefromPointingColumn(*originalMS_, accessorId_ );
338 :
339 0 : } else if (directionColumnName_ == "SOURCE_OFFSET") {
340 0 : accessor_ = sourceOffsetAccessor;
341 0 : accessorId_ = SOURCE_OFFSET;
342 0 : initializeSplinefromPointingColumn(*originalMS_, accessorId_ );
343 :
344 0 : } else if (directionColumnName_ == "ENCODER") {
345 0 : accessor_ = encoderAccessor;
346 0 : accessorId_ = ENCODER;
347 0 : initializeSplinefromPointingColumn(*originalMS_, accessorId_ );
348 :
349 : } else {
350 0 : stringstream ss;
351 0 : ss << "Column \"" << columnNameUpcase << "\" is not supported.";
352 0 : throw AipsError(ss.str());
353 0 : }
354 :
355 : //+
356 : // CAS:8418::Limmited service,
357 : // when indeicated by flag,force to use traditional Linear Interpolation.
358 : //-
359 :
360 599 : if (getCurrentSplineObj()->isCoefficientReady() == false )
361 : {
362 100 : useSplineInterpolation_ = false;
363 : }
364 :
365 599 : debuglog << "initializeSplinefromPointingColumn, Normal End." << debugpost;
366 :
367 : // ---org code ---
368 599 : configureMovingSourceCorrection();
369 599 : }
370 :
371 599 : void PointingDirectionCalculator::setFrame(String const frameType) {
372 599 : Bool status = MDirection::getType(directionType_, frameType);
373 599 : if (!status) {
374 0 : LogIO os(LogOrigin("PointingDirectionCalculator", "setFrame", WHERE));
375 : os << LogIO::WARN << "Conversion of frame string \"" << frameType
376 : << "\" into direction type enum failed. Use J2000."
377 0 : << LogIO::POST;
378 0 : directionType_ = MDirection::J2000;
379 0 : }
380 :
381 : // create conversion engine
382 :
383 : // Accessor
384 599 : MDirection nominalInputMeasure = accessor_(*pointingColumns_, 0);
385 :
386 : // RefFrame
387 599 : MDirection::Ref outReference(directionType_, referenceFrame_);
388 :
389 : // Conversion
390 : directionConvert_ = new MDirection::Convert(nominalInputMeasure,
391 599 : outReference);
392 : // Epoch
393 599 : const MEpoch *e = dynamic_cast<const MEpoch *>(referenceFrame_.epoch());
394 : const MPosition *p =
395 599 : dynamic_cast<const MPosition *>(referenceFrame_.position());
396 : debuglog << "Conversion Setup: Epoch "
397 1797 : << e->get("s").getValue() << " " << e->getRefString() << " Position "
398 2396 : << p->get("m").getValue() << " " << p->getRefString()
399 3594 : << debugpost;
400 599 : }
401 :
402 304 : void PointingDirectionCalculator::setDirectionListMatrixShape(
403 : PointingDirectionCalculator::MatrixShape const shape) {
404 304 : shape_ = shape;
405 304 : }
406 :
407 9 : void PointingDirectionCalculator::setMovingSource(String const sourceName) {
408 18 : MDirection sourceDirection(Quantity(0.0, "deg"), Quantity(90.0, "deg"));
409 9 : sourceDirection.setRefString(sourceName);
410 9 : setMovingSource(sourceDirection);
411 9 : }
412 :
413 9 : void PointingDirectionCalculator::setMovingSource(
414 : MDirection const &sourceDirection) {
415 9 : movingSource_ = dynamic_cast<MDirection *>(sourceDirection.clone());
416 :
417 : // create conversion engine for moving source
418 9 : MDirection::Ref refAzel(MDirection::AZEL, referenceFrame_);
419 9 : movingSourceConvert_ = new MDirection::Convert(*movingSource_, refAzel);
420 :
421 9 : configureMovingSourceCorrection();
422 9 : }
423 :
424 16 : void PointingDirectionCalculator::unsetMovingSource() {
425 16 : if (!movingSource_.null()) {
426 0 : movingSource_ = nullptr;
427 : }
428 16 : }
429 :
430 :
431 304 : Matrix<Double> PointingDirectionCalculator::getDirection() {
432 304 : assert(!selectedMS_.null());
433 304 : uInt const nrow = selectedMS_->nrow();
434 304 : debuglog << "selectedMS_->nrow() = " << nrow << debugpost;
435 304 : Vector<Double> outDirectionFlattened(2 * nrow);
436 : // column major data offset and increment for outDirectionFlattened,
437 : // and output matrix shape
438 304 : uInt offset = nrow;
439 304 : uInt increment = 1;
440 : // matrix shape: number of rows is nrow and number of columns is 2
441 304 : IPosition outShape(2, nrow, 2);
442 304 : if (shape_ == PointingDirectionCalculator::ROW_MAJOR) {
443 : // column major specific offset, increment and output shape
444 18 : offset = 1;
445 18 : increment = 2;
446 : // matrix shape: number of rows is 2 and number of columns is nrow
447 18 : outShape = IPosition(2, 2, nrow);
448 : }
449 :
450 609 : for (uInt i = 0; i < numAntennaBoundary_ - 1; ++i) {
451 305 : uInt start = antennaBoundary_[i];
452 305 : uInt end = antennaBoundary_[i + 1];
453 305 : uInt currentAntenna = antennaColumn_(start);
454 :
455 305 : resetAntennaPosition(currentAntenna);
456 :
457 : debuglog << "antenna " << currentAntenna << " start " << start
458 305 : << " end " << end << debugpost;
459 305 : uInt const nrowPointing = pointingTimeUTC_.nelements();
460 305 : debuglog << "nrowPointing = " << nrowPointing << debugpost;
461 : debuglog << "pointingTimeUTC = " << min(pointingTimeUTC_) << "~"
462 305 : << max(pointingTimeUTC_) << debugpost;
463 :
464 467429 : for (uInt j = start; j < end; ++j) {
465 467124 : debuglog << "start index " << j << debugpost;
466 :
467 : // doGetDirection call //
468 :
469 467124 : Vector<Double> direction = doGetDirection(j,currentAntenna); // CAS-8418 args werre changed.
470 :
471 : debuglog << "index for lat: " << (j * increment)
472 : << " (cf. outDirectionFlattened.nelements()="
473 467124 : << outDirectionFlattened.nelements() << ")" << debugpost;
474 467124 : debuglog << "index for lon: " << (offset + j * increment)
475 467124 : << debugpost;
476 467124 : outDirectionFlattened[j * increment] = direction[0];
477 467124 : outDirectionFlattened[offset + j * increment] = direction[1];
478 467124 : }
479 305 : debuglog << "done antenna " << currentAntenna << debugpost;
480 : }
481 304 : debuglog << "done getDirection" << debugpost;
482 608 : return Matrix < Double > (outShape, outDirectionFlattened.data());
483 304 : }
484 :
485 : //**************************************************
486 : // CAS-8418
487 : // doGetDirection(irow)
488 : // - Spline Interpolation is executed, except
489 : // number of data is insufficient
490 : // - Interpolation path in the module was separated.
491 : //***************************************************
492 467124 : Vector<Double> PointingDirectionCalculator::doGetDirection(uInt irow, uInt antID) {
493 467124 : debuglog << "doGetDirection(" << irow << "," << antID << ")" << debugpost;
494 :
495 : // sec:1 Linear , Spline common//
496 467124 : Double currentTime = timeColumn_.convert(irow, MEpoch::UTC).get("s").getValue();
497 467124 : resetTime(currentTime);
498 :
499 : // search and interpolate if necessary
500 : Bool exactMatch;
501 467124 : uInt const nrowPointing = pointingTimeUTC_.nelements();
502 : // pointingTableIndexCache_ is not so effective in terms of performance
503 : // simple binary search may be enough,
504 467124 : Int index = binarySearch(exactMatch, pointingTimeUTC_, currentTime,
505 : nrowPointing, 0);
506 467124 : debuglog << "binarySearch result " << index << debugpost;
507 : debuglog << "Time " << setprecision(16) << currentTime << " idx=" << index
508 467124 : << debugpost;
509 : //+
510 : // Check section on series of pointing data
511 : // by 'Binary match' wth Pointing Table
512 : // -
513 467124 : MDirection direction;
514 467124 : assert(accessor_ != NULL);
515 467124 : if (exactMatch) {
516 360808 : debuglog << "exact match" << debugpost;
517 360808 : direction = accessor_(*pointingColumns_, index);
518 106316 : } else if (index <= 0) {
519 120 : debuglog << "take 0th row" << debugpost;
520 120 : direction = accessor_(*pointingColumns_, 0);
521 106196 : } else if (index > (Int) (nrowPointing - 1)) {
522 72 : debuglog << "take final row" << debugpost;
523 72 : direction = accessor_(*pointingColumns_, nrowPointing - 1);
524 : } else {
525 106124 : debuglog << "linear interpolation " << debugpost;
526 :
527 : // commonly used result buffer .. //
528 106124 : Vector<Double> interpolated(2);
529 :
530 : //*
531 : // CAS-8418 Time Gap
532 : //*
533 106124 : bool fgsw = getCurrentSplineObj()->isTimeGap(antID, index);
534 106124 : if((fgsw==false) && // Not specially Marled. (CAS-8418)
535 103728 : (useSplineInterpolation_) ) // SPLINE //
536 : {
537 : //+
538 : // CAS-8418:: Spline Interpolation section.
539 : //-
540 :
541 103728 : Double t0 = pointingTimeUTC_[index - 1];
542 103728 : Double dtime = (currentTime - t0) ;
543 :
544 : // determin section on 'uIndex'
545 : uInt uIndex;
546 103728 : if( index >=1 ){
547 103728 : uIndex = index-1;
548 0 : } else if (index > (Int)(nrowPointing-1) ) {
549 0 : uIndex = nrowPointing-1;
550 : } else {
551 0 : stringstream ss; ss << "BUGCHECK, never come here. " << endl;
552 0 : throw AipsError( ss.str() );
553 0 : }
554 :
555 : //+
556 : // Execute Interpolation
557 : //-
558 :
559 103728 : if(getCurrentSplineObj() == nullptr){
560 0 : stringstream ss; ss << "FAITAL ERROR: Invalid Current Spline Object pointer.";
561 0 : throw AipsError(ss.str());
562 0 : }
563 103728 : interpolated = getCurrentSplineObj()-> calculate(uIndex, dtime, antID );
564 :
565 : // obtain refType1 (original copied)//
566 :
567 103728 : MDirection dir1 = accessor_(*pointingColumns_, index - 1);
568 103728 : String dirRef1 = dir1.getRefString();
569 : MDirection::Types refType1;
570 103728 : MDirection::getType(refType1, dirRef1);
571 :
572 : // LINEAR/SPLINE common
573 : // Convert the interpolated diretion from MDirection to Vector //
574 103728 : direction = MDirection(Quantum<Vector<Double> >(interpolated, "rad"),refType1);
575 103728 : }
576 : else // LINEAR (basically same as original code)//
577 : {
578 2396 : Double t0 = pointingTimeUTC_[index - 1];
579 2396 : Double t1 = pointingTimeUTC_[index];
580 2396 : Double dt = t1 - t0;
581 :
582 : debuglog << "Interpolate between " << setprecision(16) << index - 1
583 : << " (" << t0 << ") and " << index << " (" << t1 << ")"
584 2396 : << debugpost;
585 :
586 2396 : MDirection dir1 = accessor_(*pointingColumns_, index - 1);
587 2396 : MDirection dir2 = accessor_(*pointingColumns_, index);
588 :
589 : // obtain Ref type //
590 2396 : String dirRef1 = dir1.getRefString();
591 2396 : String dirRef2 = dir2.getRefString();
592 :
593 : MDirection::Types refType1, refType2;
594 2396 : MDirection::getType(refType1, dirRef1);
595 2396 : MDirection::getType(refType2, dirRef2);
596 :
597 : debuglog << "dirRef1 = " << dirRef1 << " ("
598 2396 : << MDirection::showType(refType1) << ")" << debugpost;
599 :
600 2396 : if (dirRef1 != dirRef2) {
601 0 : MeasFrame referenceFrameLocal((pointingColumns_->timeMeas())(index),
602 0 : *(referenceFrame_.position()));
603 0 : dir2 = MDirection::Convert(dir2,
604 0 : MDirection::Ref(refType1, referenceFrameLocal))();
605 0 : }
606 :
607 4792 : Vector<Double> dirVal1 = dir1.getAngle("rad").getValue();
608 4792 : Vector<Double> dirVal2 = dir2.getAngle("rad").getValue();
609 2396 : Vector<Double> scanRate = dirVal2 - dirVal1;
610 :
611 2396 : interpolated = dirVal1 + scanRate * (currentTime - t0) / dt;
612 :
613 : // LINEAR/SPLINE common
614 : // Convert the interpolated diretion from MDirection to Vector //
615 2396 : direction = MDirection(Quantum<Vector<Double> >(interpolated, "rad"),refType1);
616 :
617 2396 : }
618 106124 : }
619 :
620 : // CAS-8418:: Ending Section,following section is same as original //
621 : debuglog << "direction = "
622 934248 : << direction.getAngle("rad").getValue() << " (unit rad reference frame "
623 934248 : << direction.getRefString()
624 1868496 : << ")" << debugpost;
625 467124 : Vector<Double> outVal(2);
626 467124 : if (direction.getRefString() == MDirection::showType(directionType_)) {
627 360384 : outVal = direction.getAngle("rad").getValue();
628 : } else {
629 106740 : MDirection converted = (*directionConvert_)(direction);
630 106740 : outVal = converted.getAngle("rad").getValue();
631 : debuglog << "converted = " << outVal << "(unit rad reference frame "
632 106740 : << converted.getRefString() << ")" << debugpost;
633 106740 : }
634 :
635 : // moving source correction
636 467124 : assert(movingSourceCorrection_ != NULL);
637 467124 : movingSourceCorrection_(movingSourceConvert_, directionConvert_, outVal);
638 :
639 934248 : return outVal;
640 467124 : }
641 :
642 0 : Vector<Double> PointingDirectionCalculator::getDirection(uInt i) {
643 0 : if (i >= selectedMS_->nrow()) {
644 0 : stringstream ss;
645 0 : ss << "Out of range row index: " << i << " (nrow for selected MS "
646 0 : << getNrowForSelectedMS() << ")" << endl;
647 0 : throw AipsError(ss.str());
648 0 : }
649 0 : debuglog << "start row " << i << debugpost;
650 0 : Int currentAntennaIndex = antennaColumn_(i);
651 : debuglog << "currentAntennaIndex = " << currentAntennaIndex
652 0 : << " lastAntennaIndex_ = " << lastAntennaIndex_ << debugpost;
653 : Double currentTime =
654 0 : timeColumn_.convert(i, MEpoch::UTC).get("s").getValue();
655 0 : resetAntennaPosition(currentAntennaIndex);
656 : debuglog << "currentTime = " << currentTime << " lastTimeStamp_ = "
657 0 : << lastTimeStamp_ << debugpost;
658 0 : if (currentTime != lastTimeStamp_) {
659 0 : resetTime(i);
660 : }
661 0 : debuglog << "doGetDirection" << debugpost;
662 :
663 0 : Vector<Double> direction = doGetDirection(i,currentAntennaIndex); // CAS-8418 args were changed
664 :
665 0 : return direction;
666 : }
667 :
668 0 : Vector<rownr_t> PointingDirectionCalculator::getRowId() {
669 0 : return selectedMS_->rowNumbers();
670 : }
671 :
672 9 : Vector<rownr_t> PointingDirectionCalculator::getRowIdForOriginalMS() {
673 9 : return selectedMS_->rowNumbers(*originalMS_, True);
674 : }
675 :
676 0 : rownr_t PointingDirectionCalculator::getRowId(uInt i) {
677 0 : return selectedMS_->rowNumbers()[i];
678 : }
679 :
680 :
681 571 : void PointingDirectionCalculator::inspectAntenna() {
682 : // selectedMS_ must be sorted by ["ANTENNA1", "TIME"]
683 571 : antennaBoundary_.resize(selectedMS_->antenna().nrow() + 1);
684 571 : antennaBoundary_ = -1;
685 571 : uInt count = 0; // CAS-8418 :: Int ->iUnt //
686 571 : antennaBoundary_[count] = 0;
687 :
688 571 : ++count;
689 :
690 571 : Vector<Int> antennaList = antennaColumn_.getColumn();
691 571 : uInt nrow = antennaList.nelements();
692 571 : Int lastAnt = antennaList[0];
693 :
694 827318 : for (uInt i = 0; i < nrow; ++i) {
695 826747 : if (antennaList[i] != lastAnt) {
696 6 : antennaBoundary_[count] = i;
697 6 : ++count;
698 6 : lastAnt = antennaList[i];
699 : }
700 : }
701 : //+
702 : // CAS-8418 access violation check
703 : //-
704 571 : assert( count <= selectedMS_->antenna().nrow() );
705 : // end add
706 :
707 571 : antennaBoundary_[count] = nrow;
708 571 : ++count;
709 571 : numAntennaBoundary_ = count;
710 571 : debuglog << "antennaBoundary_=" << antennaBoundary_ << debugpost;
711 571 : }
712 :
713 305 : void PointingDirectionCalculator::initPointingTable(Int const antennaId) {
714 306 : if (!pointingTable_.null() && !pointingColumns_.null()
715 1 : && pointingTable_->nrow() > 0
716 306 : && pointingColumns_->antennaId()(0) == antennaId) {
717 : // no need to update
718 0 : return;
719 : }
720 305 : debuglog << "update pointing table for antenna " << antennaId << debugpost;
721 305 : MSPointing original = selectedMS_->pointing();
722 610 : MSPointing selected = original(original.col("ANTENNA_ID") == antennaId);
723 305 : if (selected.nrow() == 0) {
724 : debuglog << "no rows for antenna " << antennaId << " try -1"
725 0 : << debugpost;
726 : // try ANTENNA_ID == -1
727 0 : selected = original(original.col("ANTENNA_ID") == -1);
728 :
729 : #if 0 // follwing assert() invalidate throw AipsError() //
730 : assert(selected.nrow() > 0);
731 : #endif
732 0 : if (selected.nrow() == 0) {
733 0 : stringstream ss;
734 0 : ss << "Internal Error: POINTING table has no entry for antenna "
735 0 : << antennaId << "." << endl;
736 0 : throw AipsError(ss.str());
737 0 : }
738 : }
739 305 : debuglog << "selected pointing rows " << selected.nrow() << debugpost;
740 305 : pointingTable_ = new MSPointing(selected.sort("TIME"));
741 :
742 : // attach columns
743 305 : pointingColumns_ = new MSPointingColumns(*pointingTable_);
744 :
745 : // initialize pointingTimeUTC_
746 305 : uInt const nrowPointing = pointingTable_->nrow();
747 305 : pointingTimeUTC_.resize(nrowPointing);
748 : ScalarMeasColumn<MEpoch> pointingTimeColumn =
749 305 : pointingColumns_->timeMeas();
750 960529 : for (uInt i = 0; i < nrowPointing; ++i) {
751 960224 : MEpoch e = pointingTimeColumn(i);
752 960224 : if (e.getRefString() == MEpoch::showType(MEpoch::UTC)) {
753 960224 : pointingTimeUTC_[i] = e.get("s").getValue();
754 : } else {
755 0 : pointingTimeUTC_[i] =
756 0 : MEpoch::Convert(e, MEpoch::UTC)().get("s").getValue();
757 : }
758 960224 : }
759 :
760 : // reset index cache for pointing table
761 305 : pointingTableIndexCache_ = 0;
762 :
763 305 : debuglog << "done initPointingTable" << debugpost;
764 305 : }
765 :
766 876 : void PointingDirectionCalculator::resetAntennaPosition(Int const antennaId) {
767 876 : MSAntenna antennaTable = selectedMS_->antenna();
768 876 : uInt nrow = antennaTable.nrow();
769 876 : if (antennaId < 0 || (Int) nrow <= antennaId) {
770 0 : stringstream ss;
771 0 : ss << "Internal Error: Invalid ANTENNA_ID is specified (" << antennaId
772 0 : << ")." << endl;
773 0 : throw AipsError(ss.str());
774 876 : } else if (antennaId != lastAntennaIndex_ || lastAntennaIndex_ == -1) {
775 : ScalarMeasColumn < MPosition
776 305 : > antennaPositionColumn(antennaTable, "POSITION");
777 305 : antennaPosition_ = antennaPositionColumn(antennaId);
778 : debuglog << "antenna position: "
779 305 : << antennaPosition_.getRefString() << " "
780 610 : << setprecision(16) << antennaPosition_.get("m").getValue() << debugpost;
781 305 : referenceFrame_.resetPosition(antennaPosition_);
782 :
783 305 : initPointingTable(antennaId);
784 :
785 305 : lastAntennaIndex_ = antennaId;
786 305 : }
787 876 : }
788 :
789 467124 : void PointingDirectionCalculator::resetTime(Double const timestamp) {
790 467124 : debuglog << "resetTime(Double " << timestamp << ")" << debugpost;
791 : debuglog << "lastTimeStamp_ = " << lastTimeStamp_ << " timestamp = "
792 467124 : << timestamp << debugpost;
793 467124 : if (timestamp != lastTimeStamp_ || lastTimeStamp_ < 0.0) {
794 467124 : referenceEpoch_ = MEpoch(Quantity(timestamp, "s"), MEpoch::UTC);
795 467124 : referenceFrame_.resetEpoch(referenceEpoch_);
796 :
797 467124 : lastTimeStamp_ = timestamp;
798 : }
799 467124 : }
800 :
801 : //*********************************************************
802 : // CAS-8418::
803 : // Extended dunctions in PointingDirectionCalculator
804 : // for initializing Spline Interpolation
805 : //*********************************************************
806 :
807 : //+
808 : // Direction Columns and
809 : // AccessorId and accessor_ (function pointer)
810 : //-
811 :
812 : std::vector<string> dirColList
813 : = { "DIRECTION","TARGET","POINTING_OFFSET","SOURCE_OFFSET","ENCODER" };
814 :
815 : std::vector<PointingDirectionCalculator::ACCESSOR> accList
816 : = { directionAccessor, targetAccessor,pointingOffsetAccessor,
817 : sourceOffsetAccessor, encoderAccessor };
818 :
819 : // Column checck in Pointing Table //
820 304 : bool PointingDirectionCalculator::checkColumn(MeasurementSet const &ms,String const &columnName )
821 : {
822 304 : String columnNameUpcase = columnName;
823 304 : columnNameUpcase.upcase();
824 304 : if (true == (ms.pointing().tableDesc().isColumn(columnNameUpcase))) return true;
825 0 : else return false;
826 304 : }
827 : //+
828 : // Activate Spline Interpolation
829 : // - check specified Column if exist
830 : // - create Spline object from specified Direction column
831 : // - prepare Coeffient table for calulation.
832 : //-
833 599 : bool PointingDirectionCalculator::initializeSplinefromPointingColumn(MeasurementSet const &ms,
834 : PointingDirectionCalculator::PtColID DirColNo )
835 : {
836 599 : debuglog << "initializeSplinefromPointingColumn, columNo=" << DirColNo << debugpost;
837 :
838 599 : String colName = dirColList[DirColNo] ;
839 599 : PointingDirectionCalculator::ACCESSOR acc = accList[DirColNo] ;
840 :
841 : //+
842 : // Column Range check
843 : //-
844 599 : if( DirColNo >PointingDirectionCalculator::PtColID::nItems )
845 : {
846 0 : stringstream ss;
847 0 : ss << "Bugcheck. No column on Pointing Table." << endl;
848 0 : throw AipsError(ss.str());
849 : return false; // Bad Param //
850 0 : }
851 :
852 : //+
853 : // CASE 1: Spline Object is already available.
854 : //-
855 599 : if( initializeReady_[DirColNo] == true )
856 : {
857 295 : debuglog << "initializeSplinefromPointingColumn, Normal,already active." << debugpost;
858 :
859 : // SWITCH Master pointer //
860 295 : currSpline_ = splineObj_[DirColNo].get();
861 295 : assert(currSpline_ !=nullptr);
862 :
863 295 : return true; // Servece already OK //
864 : }
865 : //+
866 : // CASE 2: New Direction Colomn, initialize Spline Obj.
867 : //-
868 304 : if(checkColumn(ms, colName))
869 : {
870 304 : debuglog << "Spline Obj:: attempt to construct by " << colName.c_str() << debugpost;
871 :
872 : // Temporary Obj. //
873 304 : unique_ptr<SplineInterpolation> spTemp( new SplineInterpolation(ms,acc));
874 :
875 : // Spline Available (N>4)
876 304 : coefficientReady_ [DirColNo] = spTemp-> isCoefficientReady();
877 :
878 : // move to Spline obj. //
879 304 : assert(splineObj_[DirColNo] == nullptr);
880 304 : splineObj_[DirColNo] = std::move(spTemp);
881 :
882 : // copy to Master pointer, if this is desired. //
883 304 : currSpline_ = splineObj_[DirColNo].get();
884 :
885 : // Obj. available //
886 304 : initializeReady_[DirColNo] = true;
887 :
888 304 : return true;
889 304 : }
890 :
891 : //+
892 : // Initialize Faiure.
893 : //-
894 :
895 0 : stringstream ss;
896 0 : ss << "FAILED:: No spline obj, atempted to make. No column on Pointing Table." << endl;
897 0 : throw AipsError(ss.str());
898 :
899 599 : }
900 :
901 : //+
902 : // CAS-8418::
903 : // Exporting Internal Status/Object Service
904 : //-
905 :
906 : // Export COEFF: returns Coefficitent Table of the current Spline-Object
907 0 : PointingDirectionCalculator::COEFF PointingDirectionCalculator::exportCoeff()
908 : {
909 0 : SplineInterpolation *sp = this->getCurrentSplineObj();
910 0 : return (sp->getCoeff());
911 : }
912 :
913 : // Export status of COEFF: returns true if COEFF is available.
914 0 : bool PointingDirectionCalculator::isCoefficientReady()
915 : {
916 0 : SplineInterpolation *sp = this->getCurrentSplineObj();
917 0 : return (sp->isCoefficientReady() );
918 : }
919 :
920 : //***************************************************
921 : // CAS-8418:
922 : // Antenna Boundary (for Pointing Table ) methods
923 : // - create antenna information on Pointing Table.
924 : //***************************************************
925 :
926 : //+
927 : // Antenna Boundary Class
928 : //-
929 : class AntennaBoundary {
930 : public:
931 : AntennaBoundary(casacore::MeasurementSet const &ms) ;
932 304 : ~AntennaBoundary() { };
933 :
934 : std::pair<casacore::uInt, casacore::uInt> getAntennaBoundary( casacore::uInt n );
935 :
936 304 : casacore::uInt getNumOfAntenna() {return numAntennaBoundary_ - 1;}
937 :
938 304 : casacore::MSPointing getPointingHandle() { return hPointing_; };
939 :
940 : private:
941 : // AntennaBoundary on Pointing Tablle
942 : casacore::Vector<casacore::uInt> antennaBoundary_;
943 : casacore::uInt numAntennaBoundary_;
944 :
945 : // Pointing Table handle
946 : casacore::MSPointing hPointing_;
947 :
948 : };
949 :
950 : // Constructor //
951 304 : AntennaBoundary::AntennaBoundary(MeasurementSet const &ms)
952 : {
953 : // Antenna Boundary body //
954 304 : antennaBoundary_.resize(ms.antenna().nrow() + 1);
955 304 : antennaBoundary_ = -1;
956 :
957 304 : Int count = 0;
958 304 : antennaBoundary_[count] = 0;
959 304 : ++count;
960 :
961 : // Pointing Table Handle and the Columns Handle//
962 : // with sorting by AntennaID and Time
963 :
964 304 : MSPointing hPointing_org = ms.pointing();
965 :
966 : // Sort keys //
967 304 : Block<String> sortColumns(2);
968 304 : sortColumns[0]="ANTENNA_ID";
969 304 : sortColumns[1]="TIME";
970 :
971 : // Pointing Table handle //
972 304 : MSPointing hPointingTmp(hPointing_org.sort(sortColumns));
973 304 : hPointing_ = hPointingTmp;
974 :
975 : // Column Handle //
976 : std::unique_ptr<casacore::MSPointingColumns>
977 304 : columnPointing( new casacore::MSPointingColumns( hPointing_ ));
978 :
979 : // Antenna List //
980 :
981 304 : ScalarColumn<casacore::Int> antennaColumn = columnPointing->antennaId();
982 304 : Vector<Int> antennaList = antennaColumn.getColumn();
983 :
984 304 : uInt nrow = antennaList.nelements();
985 304 : Int lastAnt = antennaList[0];
986 :
987 1111851 : for (uInt i = 0; i < nrow; ++i)
988 : {
989 1111547 : if (antennaList[i] > lastAnt)
990 : {
991 30 : antennaBoundary_[count] = i;
992 30 : count++;
993 30 : lastAnt = antennaList[i];
994 : }
995 1111517 : else if (antennaList[i] < lastAnt )
996 : {
997 0 : stringstream ss;
998 0 : ss << "Bugcheck. Bad sort in creating antenna list." << endl;
999 0 : throw AipsError(ss.str());
1000 0 : }
1001 : }
1002 :
1003 304 : antennaBoundary_[count] = nrow;
1004 304 : ++count;
1005 304 : numAntennaBoundary_ = count;
1006 :
1007 :
1008 304 : }
1009 : // getAntenaBoundary(start, end) //
1010 334 : std::pair<casacore::uInt, casacore::uInt> AntennaBoundary::getAntennaBoundary( casacore::uInt n )
1011 : {
1012 334 : std::pair<casacore::uInt, casacore::uInt> pos(antennaBoundary_[n],antennaBoundary_[n+1]);
1013 334 : return pos;
1014 : }
1015 :
1016 : //***************************************************
1017 : // CAS-8418:
1018 : // Spline Inerpolation methods
1019 : //***************************************************
1020 :
1021 : // constructor (for each accessor) //
1022 304 : SplineInterpolation::SplineInterpolation(MeasurementSet const &ms,
1023 304 : PointingDirectionCalculator::ACCESSOR accessor )
1024 : {
1025 304 : stsCofficientReady = false;
1026 304 : init(ms, accessor);
1027 304 : }
1028 :
1029 : // initialize //
1030 304 : void SplineInterpolation::init(MeasurementSet const &ms,
1031 : PointingDirectionCalculator::ACCESSOR const my_accessor)
1032 : {
1033 : // Antenna Bounday //
1034 :
1035 304 : AntennaBoundary antb(ms);
1036 304 : uInt numAnt = antb.getNumOfAntenna();
1037 :
1038 : // prepere MS handle from selectedMS_
1039 304 : MSPointing hPoint = antb.getPointingHandle();
1040 : std::unique_ptr<casacore::MSPointingColumns>
1041 304 : columnPointing( new casacore::MSPointingColumns( hPoint ));
1042 :
1043 : // Resize (top level) //
1044 304 : tmp_time. resize(numAnt);
1045 304 : tmp_dir. resize(numAnt);
1046 :
1047 : // CAS-8418 Time Gap //
1048 304 : tmp_dtime. resize(numAnt);
1049 304 : tmp_timegap. resize(numAnt);
1050 :
1051 : // Column handle (only time,direction are needed, others are reserved) //
1052 :
1053 304 : ScalarColumn<Double> pointingTime = columnPointing ->time();
1054 304 : ScalarColumn<Double> pointingInterval = columnPointing ->interval();
1055 :
1056 : // Following columns are accessed by 'accessor_' //
1057 :
1058 304 : ArrayColumn<Double> pointingDirection = columnPointing ->direction();
1059 304 : ArrayColumn<Double> pointingTarget = columnPointing ->target();
1060 304 : ArrayColumn<Double> pointingPointingOffset = columnPointing ->pointingOffset();
1061 304 : ArrayColumn<Double> pointingSourceOffset = columnPointing ->sourceOffset();
1062 304 : ArrayColumn<Double> pointingencoder = columnPointing ->encoder();
1063 :
1064 638 : for(uInt ant=0; ant <numAnt; ant++)
1065 : {
1066 : // Antenna Bounday Pos(start,end) //
1067 334 : std::pair<uInt,uInt> pos = antb.getAntennaBoundary(ant);
1068 334 : uInt startPos = pos.first;
1069 334 : uInt endPos = pos.second;
1070 :
1071 : // define size of each antenna
1072 334 : uInt size = endPos - startPos;
1073 334 : tmp_dir [ant]. resize(size);
1074 334 : tmp_time[ant]. resize(size);
1075 :
1076 : // CAS-8418 Time Gap //
1077 334 : tmp_dtime[ant]. resize(size);
1078 334 : tmp_timegap[ant]. resize(size);
1079 :
1080 : // set up //
1081 334 : std::fill(tmp_timegap[ant].begin(), tmp_timegap[ant].end(), false); // DEFAULT = false //
1082 334 : Double prv_time = pointingTime.get(startPos); // For making diff time
1083 : // for each row //
1084 1111881 : for (uInt row = startPos; row < endPos; row++)
1085 : {
1086 1111547 : uInt index = row - startPos;
1087 :
1088 : // resizei (for Dir) //
1089 1111547 : tmp_dir[ant][index].resize(2);
1090 :
1091 1111547 : Double time = pointingTime.get(row);
1092 1111547 : MDirection dir = my_accessor(*columnPointing, row);
1093 2223094 : Vector<Double> dirVal = dir.getAngle("rad").getValue();
1094 :
1095 : // set on Vector //
1096 1111547 : tmp_time[ant][index] = time;
1097 1111547 : tmp_dir [ant][index] = dirVal;
1098 :
1099 : // CAS-8418 Time Gap //
1100 :
1101 : /* Pointing Interval */
1102 1111547 : Double p_interval = pointingInterval. get(row);
1103 1111547 : p_interval = round(10000.0 * p_interval)/10000.0; // Round at 0.0001 order
1104 :
1105 : /* Measured Interval */
1106 1111547 : Double dd = time - prv_time;
1107 1111547 : dd = round(10000.0 * dd)/10000.0; // Round at 0.0001 order
1108 :
1109 1111547 : tmp_dtime[ant][index] = dd; // record backward diff //
1110 1111547 : prv_time = time;
1111 :
1112 : /* Mark the position (force to exec LINEAR) */
1113 1111547 : if(size>=7)
1114 : {
1115 1111469 : if(dd > p_interval){
1116 7042 : if((index < size-3 )&&(index >=3)) {
1117 7032 : tmp_timegap [ant][index+3] =true; // forward
1118 7032 : tmp_timegap [ant][index+2] =true;
1119 7032 : tmp_timegap [ant][index+1] =true;
1120 :
1121 7032 : tmp_timegap [ant][index ] =true; // center (detected point)
1122 :
1123 7032 : tmp_timegap [ant][index-1] =true;
1124 7032 : tmp_timegap [ant][index-2] =true;
1125 7032 : tmp_timegap [ant][index-3] =true; // backward
1126 :
1127 : }
1128 : //*
1129 : // In case, marking on edges;
1130 : // (note: time diff(=dd) is backward difference)
1131 : //*
1132 10 : else if (index == 1) {
1133 2 : tmp_timegap [ant][index+3] =true; // forward
1134 2 : tmp_timegap [ant][index+2] =true;
1135 2 : tmp_timegap [ant][index+1] =true;
1136 2 : tmp_timegap [ant][index ] =true; // center (detected point)
1137 2 : tmp_timegap [ant][index-1] =true;
1138 : }
1139 8 : else if (index == 2) {
1140 2 : tmp_timegap [ant][index+3] =true; // forward
1141 2 : tmp_timegap [ant][index+2] =true;
1142 2 : tmp_timegap [ant][index+1] =true;
1143 2 : tmp_timegap [ant][index ] =true; // center (detected point)
1144 2 : tmp_timegap [ant][index-1] =true;
1145 2 : tmp_timegap [ant][index-2] =true;
1146 : }
1147 6 : else if (index == size-3) {
1148 2 : tmp_timegap [ant][index+2] =true;
1149 2 : tmp_timegap [ant][index+1] =true;
1150 2 : tmp_timegap [ant][index ] =true; // center (detected point)
1151 2 : tmp_timegap [ant][index-1] =true;
1152 2 : tmp_timegap [ant][index-2] =true;
1153 2 : tmp_timegap [ant][index-3] =true; // backward
1154 : }
1155 4 : else if (index == size-2) {
1156 2 : tmp_timegap [ant][index+1] =true;
1157 2 : tmp_timegap [ant][index ] =true; // center (detected point)
1158 2 : tmp_timegap [ant][index-1] =true;
1159 2 : tmp_timegap [ant][index-2] =true;
1160 2 : tmp_timegap [ant][index-3] =true; // backward
1161 : }
1162 2 : else if (index == size-1) {
1163 2 : tmp_timegap [ant][index ] =true; // center (detected point)
1164 2 : tmp_timegap [ant][index-1] =true;
1165 2 : tmp_timegap [ant][index-2] =true;
1166 2 : tmp_timegap [ant][index-3] =true; // backward
1167 : }
1168 : }
1169 : }
1170 1111547 : }
1171 : }
1172 :
1173 : //+
1174 : // Minimum Condition Inspection
1175 : // (N >=4)
1176 : //-
1177 :
1178 588 : for(uInt ant=0; ant <numAnt; ant++)
1179 : {
1180 334 : uInt s_time =tmp_time[ant].size();
1181 334 : uInt s_dir =tmp_dir [ant].size();
1182 334 : if ((s_time < 4)||(s_dir < 4))
1183 : {
1184 : // Warning .. //
1185 100 : LogIO os(LogOrigin("SplineInterpolation", "init()", WHERE));
1186 : os << LogIO::WARN << "INSUFFICIENT NUMBER OF POINTING DATA, must be ge. 4 "
1187 50 : << "Alternatively, Linear Interpolation will be used. " << LogIO::POST;
1188 :
1189 50 : stsCofficientReady = false; // initially in-usable ..
1190 50 : return;
1191 50 : }
1192 : }
1193 :
1194 : //+
1195 : // SDPosInterpolator Objct
1196 : // - create Coefficient Table -
1197 : //-
1198 254 : SDPosInterpolator sdp (tmp_time, tmp_dir);
1199 :
1200 : // Obtain Coeff (copy object) //
1201 254 : coeff_ = sdp.getSplineCoeff();
1202 :
1203 : // Table Active ..
1204 254 : stsCofficientReady = true;
1205 :
1206 : // In case COEFF inspection needed, locate dump here. //
1207 :
1208 754 : }
1209 :
1210 : //+
1211 : // Interpolation Calculation
1212 : //-
1213 103728 : casacore::Vector<casacore::Double> SplineInterpolation::calculate(uInt index,
1214 : Double dt,
1215 : uInt antID )
1216 : {
1217 103728 : debuglog << "SplineInterpolation::calculate()" << debugpost;
1218 :
1219 : // Error check //
1220 103728 : uInt arraySize = coeff_[antID].size();
1221 103728 : if( index >= arraySize)
1222 : {
1223 : // Exception handling is to be here... //
1224 0 : stringstream ss;
1225 0 : ss << "Bugcheck. Requested Index is too large." << endl;
1226 0 : throw AipsError(ss.str());
1227 0 : }
1228 : // Coefficient //
1229 :
1230 103728 : auto pCoeff_1 =coeff_[antID][index][0];
1231 103728 : auto pCoeff_2 =coeff_[antID][index][1];
1232 :
1233 103728 : Double a0 = pCoeff_1[0];
1234 103728 : Double a1 = pCoeff_1[1];
1235 103728 : Double a2 = pCoeff_1[2];
1236 103728 : Double a3 = pCoeff_1[3];
1237 :
1238 103728 : Double b0 = pCoeff_2[0];
1239 103728 : Double b1 = pCoeff_2[1];
1240 103728 : Double b2 = pCoeff_2[2];
1241 103728 : Double b3 = pCoeff_2[3];
1242 :
1243 : // Spline Calc //
1244 :
1245 103728 : Double Xs = (((0* dt + a3)*dt + a2)*dt + a1)*dt + a0;
1246 103728 : Double Ys = (((0* dt + b3)*dt + b2)*dt + b1)*dt + b0;
1247 :
1248 : // Return //
1249 :
1250 103728 : Vector<Double> outval(2);
1251 103728 : outval[0] = Xs;
1252 103728 : outval[1] = Ys;
1253 :
1254 103728 : debuglog << "SplineInterpolation::calculate() Normal return." << debugpost;
1255 :
1256 207456 : return outval;
1257 :
1258 103728 : }
1259 :
1260 :
1261 : } //# NAMESPACE CASA - END
|
