source: branches/casa-release-4_3-test02/src/STApplyCal.cpp@ 3036

Last change on this file since 3036 was 2965, checked in by Takeshi Nakazato, 11 years ago

New Development: No

JIRA Issue: Yes CAS-6585, CAS-6571

Ready for Test: Yes

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Refactoring STApplyCal::doapply.
File size: 21.6 KB
Line 
1//
2// C++ Implementation: STApplyCal
3//
4// Description:
5//
6//
7// Author: Takeshi Nakazato <takeshi.nakazato@nao.ac.jp> (C) 2012
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <assert.h>
13
14#include <casa/Arrays/Array.h>
15#include <casa/Arrays/Vector.h>
16#include <casa/Arrays/Matrix.h>
17#include <casa/Arrays/ArrayIO.h>
18#include <casa/Arrays/ArrayMath.h>
19#include <casa/BasicSL/String.h>
20#include <casa/Logging/LogIO.h>
21#include <casa/Exceptions/Error.h>
22#include <casa/Utilities/CountedPtr.h>
23#include <casa/Utilities/Sort.h>
24#include <casa/Utilities/Assert.h>
25#include <tables/Tables/Table.h>
26
27#include "Scantable.h"
28#include "STApplyCal.h"
29#include "STApplyTable.h"
30#include "STCalTsysTable.h"
31#include "STCalSkyTable.h"
32#include "STCalEnum.h"
33#include "STIdxIter.h"
34#include "Calibrator.h"
35#include "PSAlmaCalibrator.h"
36#include "Interpolator1D.h"
37#include "NearestInterpolator1D.h"
38#include "BufferedLinearInterpolator1D.h"
39#include "PolynomialInterpolator1D.h"
40#include "CubicSplineInterpolator1D.h"
41#include <atnf/PKSIO/SrcType.h>
42
43
44using namespace casa;
45using namespace std;
46
47namespace {
48template<class Accessor, class Type>
49class AccessorInterface
50{
51public:
52 typedef Type TableType;
53 static void GetSortedData(const vector<TableType *> &tablelist,
54 const Vector<uInt> &tableIndex,
55 const uInt nrow,
56 const uInt nchan,
57 Vector<Double> &time,
58 Matrix<Float> &data,
59 Matrix<uChar> &flag)
60 {
61 Vector<Double> timeUnsorted;
62 Matrix<Float> dataUnsorted;
63 Matrix<uChar> flagUnsorted;
64 GetFromTable(tablelist, tableIndex, nrow, nchan,
65 timeUnsorted, dataUnsorted, flagUnsorted);
66 SortData(timeUnsorted, dataUnsorted, flagUnsorted,
67 time, data, flag);
68 }
69private:
70 static void GetFromTable(const vector<TableType *> &tableList,
71 const Vector<uInt> &tableIndex,
72 const uInt nrow,
73 const uInt nchan,
74 Vector<Double> &time,
75 Matrix<Float> &data,
76 Matrix<uChar> &flag)
77 {
78 time.resize(nrow, False);
79 const IPosition shape(2, nrow, nchan);
80 data.resize(shape, False);
81 flag.resize(shape, False);
82 uInt rowIndex = 0;
83 for (uInt i = 0 ; i < tableIndex.nelements(); i++) {
84 TableType *p = tableList[tableIndex[i]];
85 Vector<Double> t = Accessor::GetTime(p);
86 Matrix<Float> dt = Accessor::GetData(p);
87 Matrix<uChar> fl = Accessor::GetFlag(p);
88 for (uInt j = 0; j < t.nelements(); j++) {
89 time[rowIndex] = t[j];
90 data.row(rowIndex) = dt.column(j);
91 flag.row(rowIndex) = fl.column(j);
92 rowIndex++;
93 }
94 }
95 }
96
97 static Vector<uInt> IndexSort(const Vector<Double> &t)
98 {
99 Sort sort;
100 sort.sortKey(&t[0], TpDouble, 0, Sort::Ascending);
101 Vector<uInt> idx;
102 sort.sort(idx, t.nelements(), Sort::QuickSort|Sort::NoDuplicates);
103 return idx;
104 }
105
106 static void SortData(const Vector<Double> &key, const Matrix<Float> &data,
107 const Matrix<uChar> &flag,
108 Vector<Double> &sortedKey, Matrix<Float> &sortedData,
109 Matrix<uChar> &sortedFlag)
110 {
111 Vector<uInt> sortIndex = IndexSort(key);
112 uInt len = sortIndex.nelements();
113 IPosition shape = data.shape();
114 shape[0] = len;
115 Int64 nelements = shape.product();
116 sortedKey.takeStorage(IPosition(1, len), new Double[len], TAKE_OVER);
117 sortedData.takeStorage(shape, new Float[nelements], TAKE_OVER);
118 sortedFlag.takeStorage(shape, new uChar[nelements], TAKE_OVER);
119 for (uInt i = 0 ; i < len; i++) {
120 sortedKey[i] = key[sortIndex[i]];
121 }
122 for (uInt i = 0; i < len; ++i) {
123 sortedData.row(i) = data.row(sortIndex[i]);
124 sortedFlag.row(i) = flag.row(sortIndex[i]);
125 }
126 }
127
128};
129
130class SkyTableAccessor : public AccessorInterface<SkyTableAccessor, asap::STCalSkyTable>
131{
132public:
133 static Vector<Double> GetTime(const TableType *t) {return t->getTime();}
134 static Matrix<Float> GetData(const TableType *t) {return t->getSpectra();}
135 static Matrix<uChar> GetFlag(const TableType *t) {return t->getFlagtra();}
136};
137
138class TsysTableAccessor : public AccessorInterface<TsysTableAccessor, asap::STCalTsysTable>
139{
140public:
141 static Vector<Double> GetTime(const TableType *t) {return t->getTime();}
142 static Matrix<Float> GetData(const TableType *t) {return t->getTsys();}
143 static Matrix<uChar> GetFlag(const TableType *t) {return t->getFlagtra();}
144};
145
146inline uInt setupWorkingData(const uInt n, const Double *xin, const Float *yin,
147 const uChar *f, Double *xout, Float *yout)
148{
149 uInt nValid = 0;
150 for (uInt i = 0; i < n; ++i) {
151 if (f[i] == 0) {
152 xout[nValid] = xin[i];
153 yout[nValid] = yin[i];
154 nValid++;
155 }
156 }
157 return nValid;
158}
159
160template<class InterpolationHelperImpl>
161class InterpolationHelperInterface
162{
163public:
164 static void Interpolate(const Double xref, const uInt nx, const uInt ny,
165 Double *xin, Float *yin, uChar *fin,
166 asap::Interpolator1D<Double, Float> *interpolator,
167 Double *xwork, Float *ywork,
168 Float *yout, uChar *fout)
169 {
170 for (uInt i = 0; i < ny; i++) {
171 Float *tmpY = &(yin[i * nx]);
172 uInt wnrow = setupWorkingData(nx, xin, tmpY, &(fin[i * nx]), xwork, ywork);
173 if (wnrow > 0) {
174 // any valid reference data
175 InterpolationHelperImpl::ProcessValid(xref, i, interpolator,
176 xwork, ywork, wnrow,
177 yout, fout);
178 }
179 else {
180 // no valid reference data
181 InterpolationHelperImpl::ProcessInvalid(xref, i, interpolator,
182 xin, tmpY, nx,
183 yout, fout);
184 }
185 }
186 }
187};
188
189class SkyInterpolationHelper : public InterpolationHelperInterface<SkyInterpolationHelper>
190{
191public:
192 static void ProcessValid(const Double xref, const uInt index,
193 asap::Interpolator1D<Double, Float> *interpolator,
194 Double *xwork, Float *ywork,
195 const uInt wnrow, Float *yout, uChar *fout)
196 {
197 interpolator->setData(xwork, ywork, wnrow);
198 yout[index] = interpolator->interpolate(xref);
199 }
200 static void ProcessInvalid(const Double xref, const uInt index,
201 asap::Interpolator1D<Double, Float> *interpolator,
202 Double *xwork, Float *ywork,
203 const uInt wnrow, Float *yout, uChar *fout)
204 {
205 // interpolate data regardless of flag
206 ProcessValid(xref, index, interpolator, xwork, ywork, wnrow, yout, fout);
207 // flag this channel for calibrated data
208 fout[index] = 1 << 7; // user flag
209 }
210};
211
212class TsysInterpolationHelper : public InterpolationHelperInterface<TsysInterpolationHelper>
213{
214public:
215 static void ProcessValid(const Double xref, const uInt index,
216 asap::Interpolator1D<Double, Float> *interpolator,
217 Double *xwork, Float *ywork,
218 const uInt wnrow, Float *yout, uChar *fout)
219 {
220 interpolator->setData(xwork, ywork, wnrow);
221 yout[index] = interpolator->interpolate(xref);
222 fout[index] = 0;
223 }
224 static void ProcessInvalid(const Double xref, const uInt index,
225 asap::Interpolator1D<Double, Float> *interpolator,
226 Double *xwork, Float *ywork,
227 const uInt wnrow, Float *yout, uChar *fout)
228 {
229 fout[index] = 1 << 7; // user flag
230 }
231};
232}
233
234namespace asap {
235STApplyCal::STApplyCal()
236{
237 init();
238}
239
240STApplyCal::STApplyCal(CountedPtr<Scantable> target)
241 : target_(target)
242{
243 init();
244}
245
246STApplyCal::~STApplyCal()
247{
248}
249
250void STApplyCal::init()
251{
252 caltype_ = STCalEnum::NoType;
253 doTsys_ = False;
254 iTime_ = STCalEnum::DefaultInterpolation;
255 iFreq_ = STCalEnum::DefaultInterpolation;
256}
257
258void STApplyCal::reset()
259{
260 // call init
261 init();
262
263 // clear apply tables
264 // do not delete object here
265 skytable_.resize(0);
266 tsystable_.resize(0);
267
268 // clear mapping for Tsys transfer
269 spwmap_.clear();
270
271 // reset selector
272 sel_.reset();
273
274 // delete interpolators
275 interpolatorT_ = 0;
276 interpolatorS_ = 0;
277 interpolatorF_ = 0;
278
279 // clear working scantable
280 work_ = 0;
281
282 // clear calibrator
283 calibrator_ = 0;
284}
285
286void STApplyCal::completeReset()
287{
288 reset();
289 target_ = 0;
290}
291
292void STApplyCal::setTarget(CountedPtr<Scantable> target)
293{
294 target_ = target;
295}
296
297void STApplyCal::setTarget(const String &name)
298{
299 // always create PlainTable
300 target_ = new Scantable(name, Table::Plain);
301}
302
303void STApplyCal::push(STCalSkyTable *table)
304{
305 os_.origin(LogOrigin("STApplyCal","push",WHERE));
306 skytable_.push_back(table);
307 STCalEnum::CalType caltype = STApplyTable::getCalType(table);
308 os_ << "caltype=" << caltype << LogIO::POST;
309 if (caltype_ == STCalEnum::NoType ||
310 caltype_ == STCalEnum::DefaultType ||
311 caltype_ == STCalEnum::CalTsys) {
312 caltype_ = caltype;
313 }
314 os_ << "caltype_=" << caltype_ << LogIO::POST;
315}
316
317void STApplyCal::push(STCalTsysTable *table)
318{
319 tsystable_.push_back(table);
320 doTsys_ = True;
321}
322
323void STApplyCal::setTimeInterpolation(STCalEnum::InterpolationType itype, Int order)
324{
325 iTime_ = itype;
326 order_ = order;
327}
328
329void STApplyCal::setFrequencyInterpolation(STCalEnum::InterpolationType itype, Int order)
330{
331 iFreq_ = itype;
332 order_ = order;
333}
334
335void STApplyCal::setTsysTransfer(uInt from, Vector<uInt> to)
336{
337 os_.origin(LogOrigin("STApplyCal","setTsysTransfer",WHERE));
338 os_ << "from=" << from << ", to=" << to << LogIO::POST;
339 map<uInt, Vector<uInt> >::iterator i = spwmap_.find(from);
340 if (i == spwmap_.end()) {
341 spwmap_.insert(pair<uInt, Vector<uInt> >(from, to));
342 }
343 else {
344 Vector<uInt> toNew = i->second;
345 spwmap_.erase(i);
346 uInt k = toNew.nelements();
347 toNew.resize(k+to.nelements(), True);
348 for (uInt i = 0; i < to.nelements(); i++)
349 toNew[i+k] = to[i];
350 spwmap_.insert(pair<uInt, Vector<uInt> >(from, toNew));
351 }
352}
353
354void STApplyCal::apply(Bool insitu, Bool filltsys)
355{
356 os_.origin(LogOrigin("STApplyCal","apply",WHERE));
357
358 //assert(!target_.null());
359 assert_<AipsError>(!target_.null(),"You have to set target scantable first.");
360
361 // calibrator
362 if (caltype_ == STCalEnum::CalPSAlma)
363 calibrator_ = new PSAlmaCalibrator();
364
365 // interpolator
366 initInterpolator();
367
368 // select data
369 sel_.reset();
370 sel_ = target_->getSelection();
371 if (caltype_ == STCalEnum::CalPSAlma ||
372 caltype_ == STCalEnum::CalPS) {
373 sel_.setTypes(vector<int>(1,(int)SrcType::PSON));
374 }
375 target_->setSelection(sel_);
376
377 //os_ << "sel_.print()=" << sel_.print() << LogIO::POST;
378
379 // working data
380 if (insitu) {
381 os_.origin(LogOrigin("STApplyCal","apply",WHERE));
382 os_ << "Overwrite input scantable" << LogIO::POST;
383 work_ = target_;
384 }
385 else {
386 os_.origin(LogOrigin("STApplyCal","apply",WHERE));
387 os_ << "Create output scantable from input" << LogIO::POST;
388 work_ = new Scantable(*target_, false);
389 }
390
391 //os_ << "work_->nrow()=" << work_->nrow() << LogIO::POST;
392
393 // list of apply tables for sky calibration
394 Vector<uInt> skycalList(skytable_.size());
395 uInt numSkyCal = 0;
396
397 // list of apply tables for Tsys calibration
398 for (uInt i = 0 ; i < skytable_.size(); i++) {
399 STCalEnum::CalType caltype = STApplyTable::getCalType(skytable_[i]);
400 if (caltype == caltype_) {
401 skycalList[numSkyCal] = i;
402 numSkyCal++;
403 }
404 }
405 skycalList.resize(numSkyCal, True);
406
407
408 vector<string> cols( 3 ) ;
409 cols[0] = "BEAMNO" ;
410 cols[1] = "POLNO" ;
411 cols[2] = "IFNO" ;
412 CountedPtr<STIdxIter2> iter = new STIdxIter2(work_, cols) ;
413 double start = mathutil::gettimeofday_sec();
414 os_ << LogIO::DEBUGGING << "start iterative doapply: " << start << LogIO::POST;
415 while (!iter->pastEnd()) {
416 Record ids = iter->currentValue();
417 Vector<uInt> rows = iter->getRows(SHARE);
418 if (rows.nelements() > 0)
419 doapply(ids.asuInt("BEAMNO"), ids.asuInt("IFNO"), ids.asuInt("POLNO"), rows, skycalList, filltsys);
420 iter->next();
421 }
422 double end = mathutil::gettimeofday_sec();
423 os_ << LogIO::DEBUGGING << "end iterative doapply: " << end << LogIO::POST;
424 os_ << LogIO::DEBUGGING << "elapsed time for doapply: " << end - start << " sec" << LogIO::POST;
425
426 target_->unsetSelection();
427}
428
429void STApplyCal::doapply(uInt beamno, uInt ifno, uInt polno,
430 Vector<uInt> &rows,
431 Vector<uInt> &skylist,
432 Bool filltsys)
433{
434 os_.origin(LogOrigin("STApplyCal","doapply",WHERE));
435 Bool doTsys = doTsys_;
436
437 STSelector sel;
438 vector<int> id(1);
439 id[0] = beamno;
440 sel.setBeams(id);
441 id[0] = ifno;
442 sel.setIFs(id);
443 id[0] = polno;
444 sel.setPolarizations(id);
445
446 // apply selection to apply tables
447 uInt nrowSkyTotal = 0;
448 uInt nrowTsysTotal = 0;
449 for (uInt i = 0; i < skylist.nelements(); i++) {
450 skytable_[skylist[i]]->setSelection(sel);
451 nrowSkyTotal += skytable_[skylist[i]]->nrow();
452 os_ << "nrowSkyTotal=" << nrowSkyTotal << LogIO::POST;
453 }
454
455 // Skip IFNO without sky data
456 if (nrowSkyTotal == 0)
457 return;
458
459 uInt nchanTsys = 0;
460 Vector<Double> ftsys;
461 uInt tsysifno = getIFForTsys(ifno);
462 os_ << "tsysifno=" << (Int)tsysifno << LogIO::POST;
463 if (tsystable_.size() == 0) {
464 os_.origin(LogOrigin("STApplyTable", "doapply", WHERE));
465 os_ << "No Tsys tables are given. Skip Tsys calibratoin." << LogIO::POST;
466 doTsys = False;
467 }
468 else if (tsysifno == (uInt)-1) {
469 os_.origin(LogOrigin("STApplyTable", "doapply", WHERE));
470 os_ << "No corresponding Tsys for IFNO " << ifno << ". Skip Tsys calibration" << LogIO::POST;
471 doTsys = False;
472 }
473 else {
474 id[0] = (int)tsysifno;
475 sel.setIFs(id);
476 for (uInt i = 0; i < tsystable_.size() ; i++) {
477 tsystable_[i]->setSelection(sel);
478 uInt nrowThisTsys = tsystable_[i]->nrow();
479 nrowTsysTotal += nrowThisTsys;
480 if (nrowThisTsys > 0 && nchanTsys == 0) {
481 nchanTsys = tsystable_[i]->nchan(tsysifno);
482 ftsys = tsystable_[i]->getBaseFrequency(0);
483 }
484 }
485 }
486
487 uInt nchanSp = skytable_[skylist[0]]->nchan(ifno);
488 uInt nrowSky = nrowSkyTotal;
489 Vector<Double> timeSky;
490 Matrix<Float> spoff;
491 Matrix<uChar> flagoff;
492 SkyTableAccessor::GetSortedData(skytable_, skylist,
493 nrowSkyTotal, nchanSp,
494 timeSky, spoff, flagoff);
495 nrowSky = timeSky.nelements();
496
497 uInt nrowTsys = nrowTsysTotal;
498 Vector<Double> timeTsys;
499 Matrix<Float> tsys;
500 Matrix<uChar> flagtsys;
501 if (doTsys) {
502 //os_ << "doTsys" << LogIO::POST;
503 Vector<uInt> tsyslist(tsystable_.size());
504 indgen(tsyslist);
505 TsysTableAccessor::GetSortedData(tsystable_, tsyslist,
506 nrowTsysTotal, nchanTsys,
507 timeTsys, tsys, flagtsys);
508 nrowTsys = timeTsys.nelements();
509 }
510
511 Table tab = work_->table();
512 ArrayColumn<Float> spCol(tab, "SPECTRA");
513 ArrayColumn<uChar> flCol(tab, "FLAGTRA");
514 ArrayColumn<Float> tsysCol(tab, "TSYS");
515 ScalarColumn<Double> timeCol(tab, "TIME");
516
517 // Array for scaling factor (aka Tsys)
518 Vector<Float> iTsys(IPosition(1, nchanSp), new Float[nchanSp], TAKE_OVER);
519 // Array for Tsys interpolation
520 // This is empty array and is never referenced if doTsys == false
521 // (i.e. nchanTsys == 0)
522 Vector<Float> iTsysT(IPosition(1, nchanTsys), new Float[nchanTsys], TAKE_OVER);
523
524 // Array for interpolated off spectrum
525 Vector<Float> iOff(IPosition(1, nchanSp), new Float[nchanSp], TAKE_OVER);
526
527 // working array for interpolation with flags
528 uInt arraySize = max(max(nrowTsys, nchanTsys), nrowSky);
529 Vector<Double> xwork(IPosition(1, arraySize), new Double[arraySize], TAKE_OVER);
530 Vector<Float> ywork(IPosition(1, arraySize), new Float[arraySize], TAKE_OVER);
531 Vector<uChar> fwork(IPosition(1, nchanTsys), new uChar[nchanTsys], TAKE_OVER);
532
533 // data array
534 Vector<Float> on(IPosition(1, nchanSp), new Float[nchanSp], TAKE_OVER);
535 Vector<uChar> flag(on.shape(), new uChar[on.shape().product()], TAKE_OVER);
536
537 for (uInt i = 0; i < rows.nelements(); i++) {
538 //os_ << "start i = " << i << " (row = " << rows[i] << ")" << LogIO::POST;
539 uInt irow = rows[i];
540
541 // target spectral data
542 spCol.get(irow, on);
543 flCol.get(irow, flag);
544 //os_ << "on=" << on[0] << LogIO::POST;
545 calibrator_->setSource(on);
546
547 // interpolation
548 Double t0 = timeCol(irow);
549 Double *xwork_p = xwork.data();
550 Float *ywork_p = ywork.data();
551 SkyInterpolationHelper::Interpolate(t0, nrowSky, nchanSp,
552 timeSky.data(), spoff.data(),
553 flagoff.data(), &(*interpolatorS_),
554 xwork_p, ywork_p,
555 iOff.data(), flag.data());
556 calibrator_->setReference(iOff);
557
558 if (doTsys) {
559 // Tsys correction
560 // interpolation on time axis
561 TsysInterpolationHelper::Interpolate(t0, nrowTsys, nchanTsys,
562 timeTsys.data(), tsys.data(),
563 flagtsys.data(), &(*interpolatorT_),
564 xwork_p, ywork_p,
565 iTsysT.data(), fwork.data());
566 uChar *fwork_p = fwork.data();
567 if (nchanSp == 1) {
568 // take average
569 iTsys[0] = mean(iTsysT);
570 }
571 else {
572 // interpolation on frequency axis
573 Vector<Double> fsp = getBaseFrequency(rows[i]);
574 uInt wnchan = setupWorkingData(nchanTsys, ftsys.data(), iTsysT.data(),
575 fwork_p, xwork_p, ywork_p);
576 interpolatorF_->setData(xwork_p, ywork_p, wnchan);
577 for (uInt ichan = 0; ichan < nchanSp; ichan++) {
578 iTsys[ichan] = interpolatorF_->interpolate(fsp[ichan]);
579 }
580 }
581 }
582 else {
583 Vector<Float> tsysInRow = tsysCol(irow);
584 if (tsysInRow.nelements() == 1) {
585 iTsys = tsysInRow[0];
586 }
587 else {
588 for (uInt ichan = 0; ichan < tsysInRow.nelements(); ++ichan)
589 iTsys[ichan] = tsysInRow[ichan];
590 }
591 }
592 //os_ << "iTsys=" << iTsys[0] << LogIO::POST;
593 calibrator_->setScaler(iTsys);
594
595 // do calibration
596 calibrator_->calibrate();
597
598 // update table
599 //os_ << "calibrated=" << calibrator_->getCalibrated()[0] << LogIO::POST;
600 spCol.put(irow, calibrator_->getCalibrated());
601 flCol.put(irow, flag);
602 if (filltsys)
603 tsysCol.put(irow, iTsys);
604 }
605
606
607 // reset selection on apply tables
608 for (uInt i = 0; i < skylist.nelements(); i++)
609 skytable_[i]->unsetSelection();
610 for (uInt i = 0; i < tsystable_.size(); i++)
611 tsystable_[i]->unsetSelection();
612
613
614 // reset interpolator
615 interpolatorS_->reset();
616 interpolatorF_->reset();
617 interpolatorT_->reset();
618}
619
620uInt STApplyCal::getIFForTsys(uInt to)
621{
622 for (map<casa::uInt, Vector<uInt> >::iterator i = spwmap_.begin();
623 i != spwmap_.end(); i++) {
624 Vector<uInt> tolist = i->second;
625 os_ << "from=" << i->first << ": tolist=" << tolist << LogIO::POST;
626 for (uInt j = 0; j < tolist.nelements(); j++) {
627 if (tolist[j] == to)
628 return i->first;
629 }
630 }
631 return (uInt)-1;
632}
633
634void STApplyCal::save(const String &name)
635{
636 //assert(!work_.null());
637 assert_<AipsError>(!work_.null(),"You have to execute apply method first.");
638
639 work_->setSelection(sel_);
640 work_->makePersistent(name);
641 work_->unsetSelection();
642}
643
644Vector<Double> STApplyCal::getBaseFrequency(uInt whichrow)
645{
646 //assert(whichrow <= (uInt)work_->nrow());
647 assert_<AipsError>(whichrow <= (uInt)work_->nrow(),"row index out of range.");
648 ROTableColumn col(work_->table(), "IFNO");
649 uInt ifno = col.asuInt(whichrow);
650 col.attach(work_->table(), "FREQ_ID");
651 uInt freqid = col.asuInt(whichrow);
652 uInt nc = work_->nchan(ifno);
653 STFrequencies ftab = work_->frequencies();
654 Double rp, rf, inc;
655 ftab.getEntry(rp, rf, inc, freqid);
656 Vector<Double> r(nc);
657 indgen(r, rf-rp*inc, inc);
658 return r;
659}
660
661void STApplyCal::initInterpolator()
662{
663 os_.origin(LogOrigin("STApplyCal","initInterpolator",WHERE));
664 int order = (order_ > 0) ? order_ : 1;
665 switch (iTime_) {
666 case STCalEnum::NearestInterpolation:
667 {
668 os_ << "use NearestInterpolator in time axis" << LogIO::POST;
669 interpolatorS_ = new NearestInterpolator1D<Double, Float>();
670 interpolatorT_ = new NearestInterpolator1D<Double, Float>();
671 break;
672 }
673 case STCalEnum::LinearInterpolation:
674 {
675 os_ << "use BufferedLinearInterpolator in time axis" << LogIO::POST;
676 interpolatorS_ = new BufferedLinearInterpolator1D<Double, Float>();
677 interpolatorT_ = new BufferedLinearInterpolator1D<Double, Float>();
678 break;
679 }
680 case STCalEnum::CubicSplineInterpolation:
681 {
682 os_ << "use CubicSplineInterpolator in time axis" << LogIO::POST;
683 interpolatorS_ = new CubicSplineInterpolator1D<Double, Float>();
684 interpolatorT_ = new CubicSplineInterpolator1D<Double, Float>();
685 break;
686 }
687 case STCalEnum::PolynomialInterpolation:
688 {
689 os_ << "use PolynomialInterpolator in time axis" << LogIO::POST;
690 if (order == 0) {
691 interpolatorS_ = new NearestInterpolator1D<Double, Float>();
692 interpolatorT_ = new NearestInterpolator1D<Double, Float>();
693 }
694 else {
695 interpolatorS_ = new PolynomialInterpolator1D<Double, Float>();
696 interpolatorT_ = new PolynomialInterpolator1D<Double, Float>();
697 interpolatorS_->setOrder(order);
698 interpolatorT_->setOrder(order);
699 }
700 break;
701 }
702 default:
703 {
704 os_ << "use BufferedLinearInterpolator in time axis" << LogIO::POST;
705 interpolatorS_ = new BufferedLinearInterpolator1D<Double, Float>();
706 interpolatorT_ = new BufferedLinearInterpolator1D<Double, Float>();
707 break;
708 }
709 }
710
711 switch (iFreq_) {
712 case STCalEnum::NearestInterpolation:
713 {
714 os_ << "use NearestInterpolator in frequency axis" << LogIO::POST;
715 interpolatorF_ = new NearestInterpolator1D<Double, Float>();
716 break;
717 }
718 case STCalEnum::LinearInterpolation:
719 {
720 os_ << "use BufferedLinearInterpolator in frequency axis" << LogIO::POST;
721 interpolatorF_ = new BufferedLinearInterpolator1D<Double, Float>();
722 break;
723 }
724 case STCalEnum::CubicSplineInterpolation:
725 {
726 os_ << "use CubicSplineInterpolator in frequency axis" << LogIO::POST;
727 interpolatorF_ = new CubicSplineInterpolator1D<Double, Float>();
728 break;
729 }
730 case STCalEnum::PolynomialInterpolation:
731 {
732 os_ << "use PolynomialInterpolator in frequency axis" << LogIO::POST;
733 if (order == 0) {
734 interpolatorF_ = new NearestInterpolator1D<Double, Float>();
735 }
736 else {
737 interpolatorF_ = new PolynomialInterpolator1D<Double, Float>();
738 interpolatorF_->setOrder(order);
739 }
740 break;
741 }
742 default:
743 {
744 os_ << "use LinearInterpolator in frequency axis" << LogIO::POST;
745 interpolatorF_ = new BufferedLinearInterpolator1D<Double, Float>();
746 break;
747 }
748 }
749}
750
751}
Note: See TracBrowser for help on using the repository browser.