Commits
24 24 | //# 520 Edgemont Road |
25 25 | //# Charlottesville, VA 22903-2475 USA |
26 26 | //# |
27 27 | //# $Id$ |
28 28 | |
29 29 | |
30 30 | |
31 31 | |
32 32 | |
33 33 | |
34 - | // for unique_ptr<> |
34 + | // for unique_ptr<> |
35 35 | // for std::pair |
36 36 | |
37 37 | |
38 38 | |
39 39 | |
40 40 | |
41 41 | |
42 42 | |
43 43 | |
44 44 | |
61 61 | // std::endl so that debug messages are sent to standard |
62 62 | // output. Otherwise, these macros basically does nothing. |
63 63 | // "Do nothing" behavior is implemented in NullLogger |
64 64 | // and its associating << operator below. |
65 65 | // |
66 66 | // Usage: |
67 67 | // Similar to standard output stream. |
68 68 | // |
69 69 | // debuglog << "Any message" << any_value << debugpost; |
70 70 | // |
71 - | |
71 + | |
72 72 | // #define DIRECTIONCALC_DEBUG |
73 73 | |
74 74 | namespace { |
75 75 | struct NullLogger { |
76 76 | }; |
77 77 | |
78 78 | template<class T> |
79 79 | inline NullLogger &operator<<(NullLogger &logger, T /*value*/) { |
80 80 | return logger; |
81 81 | } |
120 120 | Vector<Double> &direction) { |
121 121 | // moving source handling |
122 122 | // If moving source is specified, output direction list is always |
123 123 | // offset from reference position of moving source |
124 124 | |
125 125 | debuglog << "MovingSourceCorrection <Working>." << debugpost; |
126 126 | |
127 127 | // DEBUG (CAS-11818) // |
128 128 | assert( convertToCelestial != nullptr ); |
129 129 | assert( convertToAzel != nullptr ); |
130 - | |
130 + | |
131 131 | MDirection srcAzel = (*convertToAzel)(); |
132 132 | MDirection srcDirection = (*convertToCelestial)(srcAzel); |
133 133 | Vector<Double> srcDirectionVal = srcDirection.getAngle("rad").getValue(); |
134 134 | direction -= srcDirectionVal; |
135 135 | } |
136 136 | |
137 137 | inline void skipMovingSourceCorrection( |
138 138 | CountedPtr<MDirection::Convert> &/*convertToAzel*/, |
139 139 | CountedPtr<MDirection::Convert> &/*convertToCelestial*/, |
140 140 | Vector<Double> &/*direction*/) { |
143 143 | |
144 144 | // do nothing |
145 145 | } |
146 146 | } // anonymous namespace |
147 147 | |
148 148 | using namespace casacore; |
149 149 | namespace casa { |
150 150 | |
151 151 | //+ |
152 152 | // Original Constructor |
153 - | //- |
153 + | //- |
154 154 | PointingDirectionCalculator::PointingDirectionCalculator( |
155 155 | MeasurementSet const &ms) : |
156 - | originalMS_(new MeasurementSet(ms)), selectedMS_(), pointingTable_(), |
157 - | pointingColumns_(), timeColumn_(), intervalColumn_(), antennaColumn_(), |
156 + | originalMS_(new MeasurementSet(ms)), selectedMS_(), pointingTable_(), |
157 + | pointingColumns_(), timeColumn_(), intervalColumn_(), antennaColumn_(), |
158 158 | directionColumnName_(), accessor_(NULL), antennaPosition_(), referenceEpoch_(), |
159 - | referenceFrame_(referenceEpoch_, antennaPosition_), |
160 - | directionConvert_(NULL), directionType_(MDirection::J2000), movingSource_(NULL), |
161 - | movingSourceConvert_(NULL), movingSourceCorrection_(NULL), |
159 + | referenceFrame_(referenceEpoch_, antennaPosition_), |
160 + | directionConvert_(NULL), directionType_(MDirection::J2000), movingSource_(NULL), |
161 + | movingSourceConvert_(NULL), movingSourceCorrection_(NULL), |
162 162 | antennaBoundary_(), numAntennaBoundary_(0), pointingTimeUTC_(), lastTimeStamp_(-1.0), |
163 - | lastAntennaIndex_(-1), pointingTableIndexCache_(0), |
163 + | lastAntennaIndex_(-1), pointingTableIndexCache_(0), |
164 164 | shape_(PointingDirectionCalculator::COLUMN_MAJOR), |
165 165 | |
166 166 | /*CAS-8418*/ useSplineInterpolation_(true), ////!! Set, when Spline is used. //// |
167 - | /*CAS-8418*/ currSpline_(nullptr), |
167 + | /*CAS-8418*/ currSpline_(nullptr), |
168 168 | /*CAS-8418*/ splineObj_(PointingDirectionCalculator::PtColID::nItems), |
169 169 | /*CAS-8418*/ initializeReady_(PointingDirectionCalculator::PtColID::nItems,false),// Spline initialization Ready |
170 170 | /*CAS-8418*/ coefficientReady_(PointingDirectionCalculator::PtColID::nItems,false), // Spline Coefficient Ready |
171 171 | /*CAS-8418*/ accessorId_(DIRECTION) // specify default accessor ID |
172 - | { |
172 + | { |
173 173 | // -- original code -- // |
174 174 | accessor_ = directionAccessor; |
175 175 | |
176 176 | Block<String> sortColumns(2); |
177 177 | sortColumns[0] = "ANTENNA1"; |
178 178 | sortColumns[1] = "TIME"; |
179 179 | |
180 180 | selectedMS_ = new MeasurementSet(originalMS_->sort(sortColumns)); |
181 181 | debuglog << "Calling init()" << debugpost; |
182 182 | init(); |
304 304 | columnNameUpcase.upcase(); |
305 305 | if (!(originalMS_->pointing().tableDesc().isColumn(columnNameUpcase))) { |
306 306 | stringstream ss; |
307 307 | ss << "Column \"" << columnNameUpcase |
308 308 | << "\" doesn't exist in POINTING table."; |
309 309 | throw AipsError(ss.str()); |
310 310 | } |
311 311 | |
312 312 | directionColumnName_ = columnNameUpcase; |
313 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. |
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 318 | // Once generated and from next time, lately created Obj. is pointed and used. |
319 319 | // - from init(), this function is called. |
320 320 | // - See indetal in; |
321 - | // initializeSplinefromPointingColumn(*originalMS_, accessorId_ ); |
321 + | // initializeSplinefromPointingColumn(*originalMS_, accessorId_ ); |
322 322 | //- |
323 323 | |
324 324 | if (directionColumnName_ == "DIRECTION") { |
325 325 | accessor_ = directionAccessor; |
326 326 | accessorId_ = DIRECTION; |
327 327 | initializeSplinefromPointingColumn(*originalMS_, accessorId_ ); |
328 328 | |
329 329 | } else if (directionColumnName_ == "TARGET") { |
330 330 | accessor_ = targetAccessor; |
331 331 | accessorId_ = TARGET; |
332 332 | initializeSplinefromPointingColumn(*originalMS_, accessorId_ ); |
333 - | |
333 + | |
334 334 | } else if (directionColumnName_ == "POINTING_OFFSET") { |
335 335 | accessor_ = pointingOffsetAccessor; |
336 336 | accessorId_ = POINTING_OFFSET; |
337 337 | initializeSplinefromPointingColumn(*originalMS_, accessorId_ ); |
338 338 | |
339 339 | } else if (directionColumnName_ == "SOURCE_OFFSET") { |
340 340 | accessor_ = sourceOffsetAccessor; |
341 341 | accessorId_ = SOURCE_OFFSET; |
342 342 | initializeSplinefromPointingColumn(*originalMS_, accessorId_ ); |
343 343 | |
347 347 | initializeSplinefromPointingColumn(*originalMS_, accessorId_ ); |
348 348 | |
349 349 | } else { |
350 350 | stringstream ss; |
351 351 | ss << "Column \"" << columnNameUpcase << "\" is not supported."; |
352 352 | throw AipsError(ss.str()); |
353 353 | } |
354 354 | |
355 355 | //+ |
356 356 | // CAS:8418::Limmited service, |
357 - | // when indeicated by flag,force to use traditional Linear Interpolation. |
357 + | // when indeicated by flag,force to use traditional Linear Interpolation. |
358 358 | //- |
359 359 | |
360 360 | if (getCurrentSplineObj()->isCoefficientReady() == false ) |
361 361 | { |
362 362 | useSplineInterpolation_ = false; |
363 363 | } |
364 364 | |
365 365 | debuglog << "initializeSplinefromPointingColumn, Normal End." << debugpost; |
366 366 | |
367 367 | // ---org code --- |
373 373 | if (!status) { |
374 374 | LogIO os(LogOrigin("PointingDirectionCalculator", "setFrame", WHERE)); |
375 375 | os << LogIO::WARN << "Conversion of frame string \"" << frameType |
376 376 | << "\" into direction type enum failed. Use J2000." |
377 377 | << LogIO::POST; |
378 378 | directionType_ = MDirection::J2000; |
379 379 | } |
380 380 | |
381 381 | // create conversion engine |
382 382 | |
383 - | // Accessor |
383 + | // Accessor |
384 384 | MDirection nominalInputMeasure = accessor_(*pointingColumns_, 0); |
385 385 | |
386 386 | // RefFrame |
387 387 | MDirection::Ref outReference(directionType_, referenceFrame_); |
388 388 | |
389 - | // Conversion |
389 + | // Conversion |
390 390 | directionConvert_ = new MDirection::Convert(nominalInputMeasure, |
391 391 | outReference); |
392 - | // Epoch |
392 + | // Epoch |
393 393 | const MEpoch *e = dynamic_cast<const MEpoch *>(referenceFrame_.epoch()); |
394 394 | const MPosition *p = |
395 395 | dynamic_cast<const MPosition *>(referenceFrame_.position()); |
396 396 | debuglog << "Conversion Setup: Epoch " |
397 397 | << e->get("s").getValue() << " " << e->getRefString() << " Position " |
398 398 | << p->get("m").getValue() << " " << p->getRefString() |
399 399 | << debugpost; |
400 400 | } |
401 401 | |
402 402 | void PointingDirectionCalculator::setDirectionListMatrixShape( |
458 458 | << " end " << end << debugpost; |
459 459 | uInt const nrowPointing = pointingTimeUTC_.nelements(); |
460 460 | debuglog << "nrowPointing = " << nrowPointing << debugpost; |
461 461 | debuglog << "pointingTimeUTC = " << min(pointingTimeUTC_) << "~" |
462 462 | << max(pointingTimeUTC_) << debugpost; |
463 463 | |
464 464 | for (uInt j = start; j < end; ++j) { |
465 465 | debuglog << "start index " << j << debugpost; |
466 466 | |
467 467 | // doGetDirection call // |
468 - | |
469 - | Vector<Double> direction = doGetDirection(j,currentAntenna); // CAS-8418 args werre changed. |
468 + | |
469 + | Vector<Double> direction = doGetDirection(j,currentAntenna); // CAS-8418 args werre changed. |
470 470 | |
471 471 | debuglog << "index for lat: " << (j * increment) |
472 472 | << " (cf. outDirectionFlattened.nelements()=" |
473 473 | << outDirectionFlattened.nelements() << ")" << debugpost; |
474 474 | debuglog << "index for lon: " << (offset + j * increment) |
475 475 | << debugpost; |
476 476 | outDirectionFlattened[j * increment] = direction[0]; |
477 477 | outDirectionFlattened[offset + j * increment] = direction[1]; |
478 478 | } |
479 479 | debuglog << "done antenna " << currentAntenna << debugpost; |
480 480 | } |
481 481 | debuglog << "done getDirection" << debugpost; |
482 482 | return Matrix < Double > (outShape, outDirectionFlattened.data()); |
483 483 | } |
484 484 | |
485 485 | //************************************************** |
486 486 | // CAS-8418 |
487 487 | // doGetDirection(irow) |
488 488 | // - Spline Interpolation is executed, except |
489 489 | // number of data is insufficient |
490 - | // - Interpolation path in the module was separated. |
490 + | // - Interpolation path in the module was separated. |
491 491 | //*************************************************** |
492 492 | Vector<Double> PointingDirectionCalculator::doGetDirection(uInt irow, uInt antID) { |
493 493 | debuglog << "doGetDirection(" << irow << "," << antID << ")" << debugpost; |
494 494 | |
495 - | // sec:1 Linear , Spline common// |
495 + | // sec:1 Linear , Spline common// |
496 496 | Double currentTime = timeColumn_.convert(irow, MEpoch::UTC).get("s").getValue(); |
497 497 | resetTime(currentTime); |
498 498 | |
499 499 | // search and interpolate if necessary |
500 500 | Bool exactMatch; |
501 501 | uInt const nrowPointing = pointingTimeUTC_.nelements(); |
502 502 | // pointingTableIndexCache_ is not so effective in terms of performance |
503 503 | // simple binary search may be enough, |
504 504 | Int index = binarySearch(exactMatch, pointingTimeUTC_, currentTime, |
505 505 | nrowPointing, 0); |
506 506 | debuglog << "binarySearch result " << index << debugpost; |
507 507 | debuglog << "Time " << setprecision(16) << currentTime << " idx=" << index |
508 508 | << debugpost; |
509 509 | //+ |
510 - | // Check section on series of pointing data |
511 - | // by 'Binary match' wth Pointing Table |
510 + | // Check section on series of pointing data |
511 + | // by 'Binary match' wth Pointing Table |
512 512 | // - |
513 513 | MDirection direction; |
514 514 | assert(accessor_ != NULL); |
515 515 | if (exactMatch) { |
516 516 | debuglog << "exact match" << debugpost; |
517 517 | direction = accessor_(*pointingColumns_, index); |
518 518 | } else if (index <= 0) { |
519 519 | debuglog << "take 0th row" << debugpost; |
520 520 | direction = accessor_(*pointingColumns_, 0); |
521 521 | } else if (index > (Int) (nrowPointing - 1)) { |
522 522 | debuglog << "take final row" << debugpost; |
523 523 | direction = accessor_(*pointingColumns_, nrowPointing - 1); |
524 524 | } else { |
525 525 | debuglog << "linear interpolation " << debugpost; |
526 526 | |
527 527 | // commonly used result buffer .. // |
528 528 | Vector<Double> interpolated(2); |
529 529 | |
530 530 | //* |
531 - | // CAS-8418 Time Gap |
531 + | // CAS-8418 Time Gap |
532 532 | //* |
533 533 | bool fgsw = getCurrentSplineObj()->isTimeGap(antID, index); |
534 534 | if((fgsw==false) && // Not specially Marled. (CAS-8418) |
535 535 | (useSplineInterpolation_) ) // SPLINE // |
536 536 | { |
537 537 | //+ |
538 538 | // CAS-8418:: Spline Interpolation section. |
539 539 | //- |
540 540 | |
541 541 | Double t0 = pointingTimeUTC_[index - 1]; |
542 542 | Double dtime = (currentTime - t0) ; |
543 - | |
543 + | |
544 544 | // determin section on 'uIndex' |
545 545 | uInt uIndex; |
546 546 | if( index >=1 ){ |
547 547 | uIndex = index-1; |
548 548 | } else if (index > (Int)(nrowPointing-1) ) { |
549 549 | uIndex = nrowPointing-1; |
550 550 | } else { |
551 551 | stringstream ss; ss << "BUGCHECK, never come here. " << endl; |
552 552 | throw AipsError( ss.str() ); |
553 553 | } |
554 554 | |
555 555 | //+ |
556 - | // Execute Interpolation |
556 + | // Execute Interpolation |
557 557 | //- |
558 558 | |
559 559 | if(getCurrentSplineObj() == nullptr){ |
560 560 | stringstream ss; ss << "FAITAL ERROR: Invalid Current Spline Object pointer."; |
561 561 | throw AipsError(ss.str()); |
562 562 | } |
563 563 | interpolated = getCurrentSplineObj()-> calculate(uIndex, dtime, antID ); |
564 564 | |
565 565 | // obtain refType1 (original copied)// |
566 566 | |
567 567 | MDirection dir1 = accessor_(*pointingColumns_, index - 1); |
568 568 | String dirRef1 = dir1.getRefString(); |
569 569 | MDirection::Types refType1; |
570 570 | MDirection::getType(refType1, dirRef1); |
571 571 | |
572 - | // LINEAR/SPLINE common |
572 + | // LINEAR/SPLINE common |
573 573 | // Convert the interpolated diretion from MDirection to Vector // |
574 574 | direction = MDirection(Quantum<Vector<Double> >(interpolated, "rad"),refType1); |
575 575 | } |
576 576 | else // LINEAR (basically same as original code)// |
577 577 | { |
578 578 | Double t0 = pointingTimeUTC_[index - 1]; |
579 579 | Double t1 = pointingTimeUTC_[index]; |
580 580 | Double dt = t1 - t0; |
581 581 | |
582 582 | debuglog << "Interpolate between " << setprecision(16) << index - 1 |
583 583 | << " (" << t0 << ") and " << index << " (" << t1 << ")" |
584 584 | << debugpost; |
585 585 | |
586 586 | MDirection dir1 = accessor_(*pointingColumns_, index - 1); |
587 587 | MDirection dir2 = accessor_(*pointingColumns_, index); |
588 588 | |
589 589 | // obtain Ref type // |
590 590 | String dirRef1 = dir1.getRefString(); |
591 591 | String dirRef2 = dir2.getRefString(); |
592 - | |
592 + | |
593 593 | MDirection::Types refType1, refType2; |
594 594 | MDirection::getType(refType1, dirRef1); |
595 595 | MDirection::getType(refType2, dirRef2); |
596 - | |
596 + | |
597 597 | debuglog << "dirRef1 = " << dirRef1 << " (" |
598 598 | << MDirection::showType(refType1) << ")" << debugpost; |
599 599 | |
600 600 | if (dirRef1 != dirRef2) { |
601 601 | MeasFrame referenceFrameLocal((pointingColumns_->timeMeas())(index), |
602 602 | *(referenceFrame_.position())); |
603 603 | dir2 = MDirection::Convert(dir2, |
604 604 | MDirection::Ref(refType1, referenceFrameLocal))(); |
605 605 | } |
606 - | |
606 + | |
607 607 | Vector<Double> dirVal1 = dir1.getAngle("rad").getValue(); |
608 608 | Vector<Double> dirVal2 = dir2.getAngle("rad").getValue(); |
609 609 | Vector<Double> scanRate = dirVal2 - dirVal1; |
610 - | |
610 + | |
611 611 | interpolated = dirVal1 + scanRate * (currentTime - t0) / dt; |
612 612 | |
613 - | // LINEAR/SPLINE common |
613 + | // LINEAR/SPLINE common |
614 614 | // Convert the interpolated diretion from MDirection to Vector // |
615 615 | direction = MDirection(Quantum<Vector<Double> >(interpolated, "rad"),refType1); |
616 - | |
616 + | |
617 617 | } |
618 618 | } |
619 619 | |
620 620 | // CAS-8418:: Ending Section,following section is same as original // |
621 621 | debuglog << "direction = " |
622 622 | << direction.getAngle("rad").getValue() << " (unit rad reference frame " |
623 623 | << direction.getRefString() |
624 624 | << ")" << debugpost; |
625 625 | Vector<Double> outVal(2); |
626 626 | if (direction.getRefString() == MDirection::showType(directionType_)) { |
653 653 | Double currentTime = |
654 654 | timeColumn_.convert(i, MEpoch::UTC).get("s").getValue(); |
655 655 | resetAntennaPosition(currentAntennaIndex); |
656 656 | debuglog << "currentTime = " << currentTime << " lastTimeStamp_ = " |
657 657 | << lastTimeStamp_ << debugpost; |
658 658 | if (currentTime != lastTimeStamp_) { |
659 659 | resetTime(i); |
660 660 | } |
661 661 | debuglog << "doGetDirection" << debugpost; |
662 662 | |
663 - | Vector<Double> direction = doGetDirection(i,currentAntennaIndex); // CAS-8418 args were changed |
663 + | Vector<Double> direction = doGetDirection(i,currentAntennaIndex); // CAS-8418 args were changed |
664 664 | |
665 665 | return direction; |
666 666 | } |
667 667 | |
668 668 | Vector<uInt> PointingDirectionCalculator::getRowId() { |
669 669 | return selectedMS_->rowNumbers(); |
670 670 | } |
671 671 | |
672 672 | Vector<uInt> PointingDirectionCalculator::getRowIdForOriginalMS() { |
673 673 | return selectedMS_->rowNumbers(*originalMS_, True); |
695 695 | if (antennaList[i] != lastAnt) { |
696 696 | antennaBoundary_[count] = i; |
697 697 | ++count; |
698 698 | lastAnt = antennaList[i]; |
699 699 | } |
700 700 | } |
701 701 | //+ |
702 702 | // CAS-8418 access violation check |
703 703 | //- |
704 704 | assert( count <= selectedMS_->antenna().nrow() ); |
705 - | // end add |
705 + | // end add |
706 706 | |
707 707 | antennaBoundary_[count] = nrow; |
708 708 | ++count; |
709 709 | numAntennaBoundary_ = count; |
710 710 | debuglog << "antennaBoundary_=" << antennaBoundary_ << debugpost; |
711 711 | } |
712 712 | |
713 713 | void PointingDirectionCalculator::initPointingTable(Int const antennaId) { |
714 714 | if (!pointingTable_.null() && !pointingColumns_.null() |
715 715 | && pointingTable_->nrow() > 0 |
718 718 | return; |
719 719 | } |
720 720 | debuglog << "update pointing table for antenna " << antennaId << debugpost; |
721 721 | MSPointing original = selectedMS_->pointing(); |
722 722 | MSPointing selected = original(original.col("ANTENNA_ID") == antennaId); |
723 723 | if (selected.nrow() == 0) { |
724 724 | debuglog << "no rows for antenna " << antennaId << " try -1" |
725 725 | << debugpost; |
726 726 | // try ANTENNA_ID == -1 |
727 727 | selected = original(original.col("ANTENNA_ID") == -1); |
728 - | |
728 + | |
729 729 | // follwing assert() invalidate throw AipsError() // |
730 730 | assert(selected.nrow() > 0); |
731 - | |
731 + | |
732 732 | if (selected.nrow() == 0) { |
733 733 | stringstream ss; |
734 734 | ss << "Internal Error: POINTING table has no entry for antenna " |
735 735 | << antennaId << "." << endl; |
736 736 | throw AipsError(ss.str()); |
737 737 | } |
738 738 | } |
739 739 | debuglog << "selected pointing rows " << selected.nrow() << debugpost; |
740 740 | pointingTable_ = new MSPointing(selected.sort("TIME")); |
741 741 | |
798 798 | } |
799 799 | } |
800 800 | |
801 801 | //********************************************************* |
802 802 | // CAS-8418:: |
803 803 | // Extended dunctions in PointingDirectionCalculator |
804 804 | // for initializing Spline Interpolation |
805 805 | //********************************************************* |
806 806 | |
807 807 | //+ |
808 - | // Direction Columns and |
809 - | // AccessorId and accessor_ (function pointer) |
808 + | // Direction Columns and |
809 + | // AccessorId and accessor_ (function pointer) |
810 810 | //- |
811 811 | |
812 812 | std::vector<string> dirColList |
813 813 | = { "DIRECTION","TARGET","POINTING_OFFSET","SOURCE_OFFSET","ENCODER" }; |
814 814 | |
815 815 | std::vector<PointingDirectionCalculator::ACCESSOR> accList |
816 816 | = { directionAccessor, targetAccessor,pointingOffsetAccessor, |
817 817 | sourceOffsetAccessor, encoderAccessor }; |
818 818 | |
819 819 | // Column checck in Pointing Table // |
820 820 | bool PointingDirectionCalculator::checkColumn(MeasurementSet const &ms,String const &columnName ) |
821 821 | { |
822 822 | String columnNameUpcase = columnName; |
823 823 | columnNameUpcase.upcase(); |
824 824 | if (true == (ms.pointing().tableDesc().isColumn(columnNameUpcase))) return true; |
825 825 | else return false; |
826 826 | } |
827 827 | //+ |
828 828 | // Activate Spline Interpolation |
829 829 | // - check specified Column if exist |
830 830 | // - create Spline object from specified Direction column |
831 - | // - prepare Coeffient table for calulation. |
831 + | // - prepare Coeffient table for calulation. |
832 832 | //- |
833 833 | bool PointingDirectionCalculator::initializeSplinefromPointingColumn(MeasurementSet const &ms, |
834 834 | PointingDirectionCalculator::PtColID DirColNo ) |
835 835 | { |
836 836 | debuglog << "initializeSplinefromPointingColumn, columNo=" << DirColNo << debugpost; |
837 837 | |
838 838 | String colName = dirColList[DirColNo] ; |
839 839 | PointingDirectionCalculator::ACCESSOR acc = accList[DirColNo] ; |
840 840 | |
841 841 | //+ |
842 - | // Column Range check |
842 + | // Column Range check |
843 843 | //- |
844 844 | if( DirColNo >PointingDirectionCalculator::PtColID::nItems ) |
845 845 | { |
846 846 | stringstream ss; |
847 847 | ss << "Bugcheck. No column on Pointing Table." << endl; |
848 848 | throw AipsError(ss.str()); |
849 849 | return false; // Bad Param // |
850 850 | } |
851 851 | |
852 852 | //+ |
853 853 | // CASE 1: Spline Object is already available. |
854 854 | //- |
855 855 | if( initializeReady_[DirColNo] == true ) |
856 856 | { |
857 857 | debuglog << "initializeSplinefromPointingColumn, Normal,already active." << debugpost; |
858 858 | |
859 - | // SWITCH Master pointer // |
859 + | // SWITCH Master pointer // |
860 860 | currSpline_ = splineObj_[DirColNo].get(); |
861 861 | assert(currSpline_ !=nullptr); |
862 862 | |
863 863 | return true; // Servece already OK // |
864 864 | } |
865 865 | //+ |
866 - | // CASE 2: New Direction Colomn, initialize Spline Obj. |
866 + | // CASE 2: New Direction Colomn, initialize Spline Obj. |
867 867 | //- |
868 868 | if(checkColumn(ms, colName)) |
869 869 | { |
870 870 | debuglog << "Spline Obj:: attempt to construct by " << colName.c_str() << debugpost; |
871 871 | |
872 872 | // Temporary Obj. // |
873 873 | unique_ptr<SplineInterpolation> spTemp( new SplineInterpolation(ms,acc)); |
874 874 | |
875 875 | // Spline Available (N>4) |
876 876 | coefficientReady_ [DirColNo] = spTemp-> isCoefficientReady(); |
877 877 | |
878 878 | // move to Spline obj. // |
879 879 | assert(splineObj_[DirColNo] == false); // MUST BE CLEAN (see C++11 for this statement) // |
880 880 | splineObj_[DirColNo] = std::move(spTemp); |
881 881 | |
882 - | // copy to Master pointer, if this is desired. // |
882 + | // copy to Master pointer, if this is desired. // |
883 883 | currSpline_ = splineObj_[DirColNo].get(); |
884 884 | |
885 885 | // Obj. available // |
886 886 | initializeReady_[DirColNo] = true; |
887 887 | |
888 888 | return true; |
889 - | } |
890 - | |
889 + | } |
890 + | |
891 891 | //+ |
892 892 | // Initialize Faiure. |
893 893 | //- |
894 - | |
894 + | |
895 895 | stringstream ss; |
896 896 | ss << "FAILED:: No spline obj, atempted to make. No column on Pointing Table." << endl; |
897 897 | throw AipsError(ss.str()); |
898 - | |
898 + | |
899 899 | } |
900 900 | |
901 901 | //+ |
902 902 | // CAS-8418:: |
903 - | // Exporting Internal Status/Object Service |
904 - | //- |
903 + | // Exporting Internal Status/Object Service |
904 + | //- |
905 905 | |
906 906 | // Export COEFF: returns Coefficitent Table of the current Spline-Object |
907 907 | PointingDirectionCalculator::COEFF PointingDirectionCalculator::exportCoeff() |
908 908 | { |
909 909 | SplineInterpolation *sp = this->getCurrentSplineObj(); |
910 910 | return (sp->getCoeff()); |
911 911 | } |
912 912 | |
913 913 | // Export status of COEFF: returns true if COEFF is available. |
914 914 | bool PointingDirectionCalculator::isCoefficientReady() |
915 915 | { |
916 916 | SplineInterpolation *sp = this->getCurrentSplineObj(); |
917 917 | return (sp->isCoefficientReady() ); |
918 918 | } |
919 919 | |
920 920 | //*************************************************** |
921 - | // CAS-8418: |
921 + | // CAS-8418: |
922 922 | // Antenna Boundary (for Pointing Table ) methods |
923 923 | // - create antenna information on Pointing Table. |
924 924 | //*************************************************** |
925 925 | |
926 926 | //+ |
927 927 | // Antenna Boundary Class |
928 928 | //- |
929 929 | class AntennaBoundary { |
930 930 | public: |
931 931 | AntennaBoundary(casacore::MeasurementSet const &ms) ; |
932 932 | ~AntennaBoundary() { }; |
933 933 | |
934 934 | std::pair<casacore::uInt, casacore::uInt> getAntennaBoundary( casacore::uInt n ); |
935 - | |
936 - | casacore::uInt getNumOfAntenna() {return numAntennaBoundary_ - 1;} |
935 + | |
936 + | casacore::uInt getNumOfAntenna() {return numAntennaBoundary_ - 1;} |
937 937 | |
938 938 | casacore::MSPointing getPointingHandle() { return hPointing_; }; |
939 939 | |
940 940 | private: |
941 - | // AntennaBoundary on Pointing Tablle |
941 + | // AntennaBoundary on Pointing Tablle |
942 942 | casacore::Vector<casacore::uInt> antennaBoundary_; |
943 943 | casacore::uInt numAntennaBoundary_; |
944 - | |
944 + | |
945 945 | // Pointing Table handle |
946 - | casacore::MSPointing hPointing_; |
947 - | |
946 + | casacore::MSPointing hPointing_; |
947 + | |
948 948 | }; |
949 949 | |
950 950 | // Constructor // |
951 - | AntennaBoundary::AntennaBoundary(MeasurementSet const &ms) |
951 + | AntennaBoundary::AntennaBoundary(MeasurementSet const &ms) |
952 952 | { |
953 953 | // Antenna Boundary body // |
954 954 | antennaBoundary_.resize(ms.antenna().nrow() + 1); |
955 955 | antennaBoundary_ = -1; |
956 956 | |
957 957 | Int count = 0; |
958 958 | antennaBoundary_[count] = 0; |
959 959 | ++count; |
960 960 | |
961 961 | // Pointing Table Handle and the Columns Handle// |
962 - | // with sorting by AntennaID and Time |
962 + | // with sorting by AntennaID and Time |
963 963 | |
964 964 | MSPointing hPointing_org = ms.pointing(); |
965 965 | |
966 966 | // Sort keys // |
967 967 | Block<String> sortColumns(2); |
968 968 | sortColumns[0]="ANTENNA_ID"; |
969 969 | sortColumns[1]="TIME"; |
970 970 | |
971 971 | // Pointing Table handle // |
972 972 | MSPointing hPointingTmp(hPointing_org.sort(sortColumns)); |
977 977 | columnPointing( new casacore::MSPointingColumns( hPointing_ )); |
978 978 | |
979 979 | // Antenna List // |
980 980 | |
981 981 | ScalarColumn<casacore::Int> antennaColumn = columnPointing->antennaId(); |
982 982 | Vector<Int> antennaList = antennaColumn.getColumn(); |
983 983 | |
984 984 | uInt nrow = antennaList.nelements(); |
985 985 | Int lastAnt = antennaList[0]; |
986 986 | |
987 - | for (uInt i = 0; i < nrow; ++i) |
988 - | { |
989 - | if (antennaList[i] > lastAnt) |
990 - | { |
991 - | antennaBoundary_[count] = i; |
987 + | for (uInt i = 0; i < nrow; ++i) |
988 + | { |
989 + | if (antennaList[i] > lastAnt) |
990 + | { |
991 + | antennaBoundary_[count] = i; |
992 992 | count++; |
993 993 | lastAnt = antennaList[i]; |
994 - | } |
994 + | } |
995 995 | else if (antennaList[i] < lastAnt ) |
996 - | { |
996 + | { |
997 997 | stringstream ss; |
998 998 | ss << "Bugcheck. Bad sort in creating antenna list." << endl; |
999 999 | throw AipsError(ss.str()); |
1000 - | } |
1000 + | } |
1001 1001 | } |
1002 - | |
1003 - | //+ |
1004 - | // CAS-8418 access violation check |
1005 - | //- |
1006 - | assert( count <= selectedMS_->antenna().nrow() ); |
1007 - | // end add |
1008 1002 | |
1009 1003 | antennaBoundary_[count] = nrow; |
1010 1004 | ++count; |
1011 1005 | numAntennaBoundary_ = count; |
1012 1006 | |
1013 1007 | |
1014 1008 | } |
1015 1009 | // getAntenaBoundary(start, end) // |
1016 1010 | std::pair<casacore::uInt, casacore::uInt> AntennaBoundary::getAntennaBoundary( casacore::uInt n ) |
1017 1011 | { |
1018 1012 | std::pair<casacore::uInt, casacore::uInt> pos(antennaBoundary_[n],antennaBoundary_[n+1]); |
1019 1013 | return pos; |
1020 1014 | } |
1021 1015 | |
1022 1016 | //*************************************************** |
1023 - | // CAS-8418: |
1017 + | // CAS-8418: |
1024 1018 | // Spline Inerpolation methods |
1025 1019 | //*************************************************** |
1026 1020 | |
1027 1021 | // constructor (for each accessor) // |
1028 - | SplineInterpolation::SplineInterpolation(MeasurementSet const &ms, |
1029 - | PointingDirectionCalculator::ACCESSOR accessor ) |
1022 + | SplineInterpolation::SplineInterpolation(MeasurementSet const &ms, |
1023 + | PointingDirectionCalculator::ACCESSOR accessor ) |
1030 1024 | { |
1031 1025 | stsCofficientReady = false; |
1032 1026 | init(ms, accessor); |
1033 1027 | } |
1034 1028 | |
1035 1029 | // initialize // |
1036 - | void SplineInterpolation::init(MeasurementSet const &ms, |
1030 + | void SplineInterpolation::init(MeasurementSet const &ms, |
1037 1031 | PointingDirectionCalculator::ACCESSOR const my_accessor) |
1038 1032 | { |
1039 1033 | // Antenna Bounday // |
1040 1034 | |
1041 1035 | AntennaBoundary antb(ms); |
1042 1036 | uInt numAnt = antb.getNumOfAntenna(); |
1043 1037 | |
1044 1038 | // prepere MS handle from selectedMS_ |
1045 1039 | MSPointing hPoint = antb.getPointingHandle(); |
1046 1040 | std::unique_ptr<casacore::MSPointingColumns> |
1047 1041 | columnPointing( new casacore::MSPointingColumns( hPoint )); |
1048 1042 | |
1049 1043 | // Resize (top level) // |
1050 1044 | tmp_time. resize(numAnt); |
1051 1045 | tmp_dir. resize(numAnt); |
1052 1046 | |
1053 1047 | // CAS-8418 Time Gap // |
1054 1048 | tmp_dtime. resize(numAnt); |
1055 1049 | tmp_timegap. resize(numAnt); |
1056 - | |
1050 + | |
1057 1051 | // Column handle (only time,direction are needed, others are reserved) // |
1058 - | |
1052 + | |
1059 1053 | ScalarColumn<Double> pointingTime = columnPointing ->time(); |
1060 1054 | ScalarColumn<Double> pointingInterval = columnPointing ->interval(); |
1061 1055 | |
1062 1056 | // Following columns are accessed by 'accessor_' // |
1063 - | |
1057 + | |
1064 1058 | ArrayColumn<Double> pointingDirection = columnPointing ->direction(); |
1065 - | ArrayColumn<Double> pointingTarget = columnPointing ->target(); |
1059 + | ArrayColumn<Double> pointingTarget = columnPointing ->target(); |
1066 1060 | ArrayColumn<Double> pointingPointingOffset = columnPointing ->pointingOffset(); |
1067 1061 | ArrayColumn<Double> pointingSourceOffset = columnPointing ->sourceOffset(); |
1068 1062 | ArrayColumn<Double> pointingencoder = columnPointing ->encoder(); |
1069 - | |
1063 + | |
1070 1064 | for(uInt ant=0; ant <numAnt; ant++) |
1071 1065 | { |
1072 1066 | // Antenna Bounday Pos(start,end) // |
1073 1067 | std::pair<uInt,uInt> pos = antb.getAntennaBoundary(ant); |
1074 1068 | uInt startPos = pos.first; |
1075 1069 | uInt endPos = pos.second; |
1076 1070 | |
1077 1071 | // define size of each antenna |
1078 - | uInt size = endPos - startPos; |
1072 + | uInt size = endPos - startPos; |
1079 1073 | tmp_dir [ant]. resize(size); |
1080 1074 | tmp_time[ant]. resize(size); |
1081 1075 | |
1082 1076 | // CAS-8418 Time Gap // |
1083 1077 | tmp_dtime[ant]. resize(size); |
1084 1078 | tmp_timegap[ant]. resize(size); |
1085 - | |
1079 + | |
1086 1080 | // set up // |
1087 - | std::fill(tmp_timegap[ant].begin(), tmp_timegap[ant].end(), false); // DEFAULT = false // |
1081 + | std::fill(tmp_timegap[ant].begin(), tmp_timegap[ant].end(), false); // DEFAULT = false // |
1088 1082 | Double prv_time = pointingTime.get(startPos); // For making diff time |
1089 - | // for each row // |
1090 - | for (uInt row = startPos; row < endPos; row++) |
1083 + | // for each row // |
1084 + | for (uInt row = startPos; row < endPos; row++) |
1091 1085 | { |
1092 1086 | uInt index = row - startPos; |
1093 1087 | |
1094 1088 | // resizei (for Dir) // |
1095 1089 | tmp_dir[ant][index].resize(2); |
1096 1090 | |
1097 1091 | Double time = pointingTime.get(row); |
1098 1092 | MDirection dir = my_accessor(*columnPointing, row); |
1099 1093 | Vector<Double> dirVal = dir.getAngle("rad").getValue(); |
1100 1094 | |
1104 1098 | |
1105 1099 | // CAS-8418 Time Gap // |
1106 1100 | |
1107 1101 | /* Pointing Interval */ |
1108 1102 | Double p_interval = pointingInterval. get(row); |
1109 1103 | p_interval = round(10000.0 * p_interval)/10000.0; // Round at 0.0001 order |
1110 1104 | |
1111 1105 | /* Measured Interval */ |
1112 1106 | Double dd = time - prv_time; |
1113 1107 | dd = round(10000.0 * dd)/10000.0; // Round at 0.0001 order |
1114 - | |
1115 - | tmp_dtime[ant][index] = dd; // record backward diff // |
1108 + | |
1109 + | tmp_dtime[ant][index] = dd; // record backward diff // |
1116 1110 | prv_time = time; |
1117 1111 | |
1118 1112 | /* Mark the position (force to exec LINEAR) */ |
1119 1113 | if(size>=7) |
1120 - | { |
1114 + | { |
1121 1115 | if(dd > p_interval){ |
1122 1116 | if((index < size-3 )&&(index >=3)) { |
1123 - | tmp_timegap [ant][index+3] =true; // forward |
1117 + | tmp_timegap [ant][index+3] =true; // forward |
1124 1118 | tmp_timegap [ant][index+2] =true; |
1125 1119 | tmp_timegap [ant][index+1] =true; |
1126 1120 | |
1127 1121 | tmp_timegap [ant][index ] =true; // center (detected point) |
1128 1122 | |
1129 1123 | tmp_timegap [ant][index-1] =true; |
1130 1124 | tmp_timegap [ant][index-2] =true; |
1131 - | tmp_timegap [ant][index-3] =true; // backward |
1132 - | |
1125 + | tmp_timegap [ant][index-3] =true; // backward |
1126 + | |
1133 1127 | } |
1134 1128 | //* |
1135 1129 | // In case, marking on edges; |
1136 1130 | // (note: time diff(=dd) is backward difference) |
1137 1131 | //* |
1138 1132 | else if (index == 1) { |
1139 - | tmp_timegap [ant][index+3] =true; // forward |
1133 + | tmp_timegap [ant][index+3] =true; // forward |
1140 1134 | tmp_timegap [ant][index+2] =true; |
1141 1135 | tmp_timegap [ant][index+1] =true; |
1142 1136 | tmp_timegap [ant][index ] =true; // center (detected point) |
1143 1137 | tmp_timegap [ant][index-1] =true; |
1144 1138 | } |
1145 1139 | else if (index == 2) { |
1146 - | tmp_timegap [ant][index+3] =true; // forward |
1140 + | tmp_timegap [ant][index+3] =true; // forward |
1147 1141 | tmp_timegap [ant][index+2] =true; |
1148 1142 | tmp_timegap [ant][index+1] =true; |
1149 1143 | tmp_timegap [ant][index ] =true; // center (detected point) |
1150 1144 | tmp_timegap [ant][index-1] =true; |
1151 1145 | tmp_timegap [ant][index-2] =true; |
1152 1146 | } |
1153 1147 | else if (index == size-3) { |
1154 1148 | tmp_timegap [ant][index+2] =true; |
1155 1149 | tmp_timegap [ant][index+1] =true; |
1156 1150 | tmp_timegap [ant][index ] =true; // center (detected point) |
1157 1151 | tmp_timegap [ant][index-1] =true; |
1158 1152 | tmp_timegap [ant][index-2] =true; |
1159 - | tmp_timegap [ant][index-3] =true; // backward |
1153 + | tmp_timegap [ant][index-3] =true; // backward |
1160 1154 | } |
1161 1155 | else if (index == size-2) { |
1162 1156 | tmp_timegap [ant][index+1] =true; |
1163 1157 | tmp_timegap [ant][index ] =true; // center (detected point) |
1164 1158 | tmp_timegap [ant][index-1] =true; |
1165 1159 | tmp_timegap [ant][index-2] =true; |
1166 - | tmp_timegap [ant][index-3] =true; // backward |
1160 + | tmp_timegap [ant][index-3] =true; // backward |
1167 1161 | } |
1168 1162 | else if (index == size-1) { |
1169 1163 | tmp_timegap [ant][index ] =true; // center (detected point) |
1170 1164 | tmp_timegap [ant][index-1] =true; |
1171 1165 | tmp_timegap [ant][index-2] =true; |
1172 - | tmp_timegap [ant][index-3] =true; // backward |
1166 + | tmp_timegap [ant][index-3] =true; // backward |
1173 1167 | } |
1174 1168 | } |
1175 - | } |
1169 + | } |
1176 1170 | } |
1177 1171 | } |
1178 1172 | |
1179 1173 | //+ |
1180 1174 | // Minimum Condition Inspection |
1181 - | // (N >=4) |
1175 + | // (N >=4) |
1182 1176 | //- |
1183 1177 | |
1184 1178 | for(uInt ant=0; ant <numAnt; ant++) |
1185 1179 | { |
1186 1180 | uInt s_time =tmp_time[ant].size(); |
1187 1181 | uInt s_dir =tmp_dir [ant].size(); |
1188 1182 | if ((s_time < 4)||(s_dir < 4)) |
1189 1183 | { |
1190 1184 | // Warning .. // |
1191 1185 | LogIO os(LogOrigin("SplineInterpolation", "init()", WHERE)); |
1192 - | os << LogIO::WARN << "INSUFFICIENT NUMBER OF POINTING DATA, must be ge. 4 " |
1186 + | os << LogIO::WARN << "INSUFFICIENT NUMBER OF POINTING DATA, must be ge. 4 " |
1193 1187 | << "Alternatively, Linear Interpolation will be used. " << LogIO::POST; |
1194 1188 | |
1195 1189 | stsCofficientReady = false; // initially in-usable .. |
1196 1190 | return; |
1197 1191 | } |
1198 1192 | } |
1199 - | |
1193 + | |
1200 1194 | //+ |
1201 - | // SDPosInterpolator Objct |
1202 - | // - create Coefficient Table - |
1195 + | // SDPosInterpolator Objct |
1196 + | // - create Coefficient Table - |
1203 1197 | //- |
1204 1198 | SDPosInterpolator sdp (tmp_time, tmp_dir); |
1205 1199 | |
1206 1200 | // Obtain Coeff (copy object) // |
1207 1201 | coeff_ = sdp.getSplineCoeff(); |
1208 1202 | |
1209 1203 | // Table Active .. |
1210 1204 | stsCofficientReady = true; |
1211 1205 | |
1212 1206 | // In case COEFF inspection needed, locate dump here. // |
1221 1215 | uInt antID ) |
1222 1216 | { |
1223 1217 | debuglog << "SplineInterpolation::calculate()" << debugpost; |
1224 1218 | |
1225 1219 | // Error check // |
1226 1220 | uInt arraySize = coeff_[antID].size(); |
1227 1221 | if( index >= arraySize) |
1228 1222 | { |
1229 1223 | // Exception handling is to be here... // |
1230 1224 | stringstream ss; |
1231 - | ss << "Bugcheck. Requested Index is too large." << endl; |
1225 + | ss << "Bugcheck. Requested Index is too large." << endl; |
1232 1226 | throw AipsError(ss.str()); |
1233 1227 | } |
1234 1228 | // Coefficient // |
1235 1229 | |
1236 1230 | auto pCoeff_1 =coeff_[antID][index][0]; |
1237 1231 | auto pCoeff_2 =coeff_[antID][index][1]; |
1238 1232 | |
1239 1233 | Double a0 = pCoeff_1[0]; |
1240 1234 | Double a1 = pCoeff_1[1]; |
1241 1235 | Double a2 = pCoeff_1[2]; |
1246 1240 | Double b2 = pCoeff_2[2]; |
1247 1241 | Double b3 = pCoeff_2[3]; |
1248 1242 | |
1249 1243 | // Spline Calc // |
1250 1244 | |
1251 1245 | Double Xs = (((0* dt + a3)*dt + a2)*dt + a1)*dt + a0; |
1252 1246 | Double Ys = (((0* dt + b3)*dt + b2)*dt + b1)*dt + b0; |
1253 1247 | |
1254 1248 | // Return // |
1255 1249 | |
1256 - | Vector<Double> outval(2); |
1250 + | Vector<Double> outval(2); |
1257 1251 | outval[0] = Xs; |
1258 1252 | outval[1] = Ys; |
1259 1253 | |
1260 1254 | debuglog << "SplineInterpolation::calculate() Normal return." << debugpost; |
1261 1255 | |
1262 1256 | return outval; |
1263 1257 | |
1264 1258 | } |
1265 1259 | |
1266 - | |
1260 + | |
1267 1261 | } //# NAMESPACE CASA - END |