source: trunk/src/STCalTsys.cpp@ 2736

Last change on this file since 2736 was 2720, checked in by Takeshi Nakazato, 12 years ago

New Development: Yes

JIRA Issue: Yes CAS-4770 and its sub-tickets

Ready for Test: Yes

Interface Changes: 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...

First version of applycal for single dish calibration.
Added new classes for the operation (STApplyCal, Calibrator, PSAlmaCalibrator,
Locator, BisectionLocator, Interpolator1D, NearestInterpolator1D).
Also, modified existing classes to fit with implementation of applycal.


File size: 3.8 KB
Line 
1//
2// C++ Implementation: STCalTsys
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
13#include <vector>
14#include "STSelector.h"
15#include "STCalTsys.h"
16#include "RowAccumulator.h"
17#include "STIdxIter.h"
18#include "STDefs.h"
19#include <atnf/PKSIO/SrcType.h>
20
21#include <casa/Arrays/ArrayMath.h>
22
23using namespace std;
24using namespace casa;
25
26namespace asap {
27 STCalTsys::STCalTsys(CountedPtr<Scantable> &s, vector<int> &iflist)
28 : STCalibration(s),
29 iflist_(iflist)
30{
31 applytable_ = new STCalTsysTable(*s);
32}
33
34void STCalTsys::calibrate()
35{
36 STSelector selOrg = scantable_->getSelection();
37 STSelector sel;
38 sel.setIFs(iflist_);
39 scantable_->setSelection(sel);
40
41 fillCalTable();
42
43 scantable_->setSelection(selOrg);
44}
45
46void STCalTsys::fillCalTable()
47{
48 RowAccumulator acc(W_TINT);
49
50 vector<string> cols(3);
51 cols[0] = "IFNO";
52 cols[1] = "POLNO";
53 cols[2] = "BEAMNO";
54 STIdxIterAcc iter(scantable_, cols);
55
56 ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(scantable_->table(), "TIME");
57 Vector<Double> timeSec = tcol->getColumn() * 86400.0;
58 tcol->attach(scantable_->table(), "INTERVAL");
59 Vector<Double> intervalSec = tcol->getColumn();
60 delete tcol;
61 ROScalarColumn<Float> *ecol = new ROScalarColumn<Float>(scantable_->table(), "ELEVATION");
62 Vector<Float> elevation = ecol->getColumn();
63 delete ecol;
64
65 ROArrayColumn<Float> specCol(scantable_->table(), "TSYS");
66 ROArrayColumn<uChar> flagCol(scantable_->table(), "FLAGTRA");
67 ROScalarColumn<uInt> freqidCol(scantable_->table(), "FREQ_ID");
68
69 // dummy Tsys: the following process doesn't need Tsys but RowAccumulator
70 // requires to set it with spectral data
71 Vector<Float> tsys(1, 1.0);
72
73 Double timeCen = 0.0;
74 Float elCen = 0.0;
75 uInt count = 0;
76
77 while(!iter.pastEnd()) {
78 Vector<uInt> rows = iter.getRows(SHARE);
79 Vector<uInt> current = iter.current();
80 uInt len = rows.nelements();
81 if (len == 0) {
82 iter.next();
83 continue;
84 }
85
86 uInt nchan = scantable_->nchan(scantable_->getIF(rows[0]));
87 Vector<uChar> flag(nchan);
88 Vector<Bool> bflag(nchan);
89 Vector<Float> spec(nchan);
90
91 Vector<Double> timeSep(len);
92 for (uInt i = 0; i < len-1; i++) {
93 timeSep[i] = timeSec[rows[i+1]] - timeSec[rows[i]] ;
94 }
95 Double tMedian = median(timeSep(IPosition(1,0), IPosition(1,len-2)));
96 timeSep[len-1] = tMedian * 10000.0 ; // any large value
97
98 uInt irow ;
99 uInt jrow ;
100 for (uInt i = 0; i < len; i++) {
101 irow = rows[i];
102 jrow = (i < len-1) ? rows[i+1] : rows[i];
103 // accumulate data
104 flagCol.get(irow, flag);
105 convertArray(bflag, flag);
106 specCol.get(irow, spec);
107 if ( !allEQ(bflag,True) )
108 acc.add( spec, !bflag, tsys, intervalSec[irow], timeSec[irow] ) ;
109 timeCen += timeSec[irow];
110 elCen += elevation[irow];
111 count++;
112
113 // check time gap
114 double gap = 2.0 * timeSep[i] / (intervalSec[jrow] + intervalSec[irow]);
115 if ( gap > 5.0 ) {
116 if ( acc.state() ) {
117 acc.replaceNaN() ;
118// const Vector<Bool> &msk = acc.getMask();
119// convertArray(flag, !msk);
120// for (uInt k = 0; k < nchan; ++k) {
121// uChar userFlag = 1 << 7;
122// if (msk[k]==True) userFlag = 0 << 7;
123// flag(k) = userFlag;
124// }
125 timeCen /= (Double)count * 86400.0; // sec->day
126 elCen /= (Float)count;
127 STCalTsysTable *p = dynamic_cast<STCalTsysTable *>(&(*applytable_));
128 p->appenddata(0, 0, current[2], current[0], current[1],
129 freqidCol(irow), timeCen, elCen, acc.getSpectrum());
130 }
131 acc.reset() ;
132 timeCen = 0.0;
133 elCen = 0.0;
134 count = 0;
135 }
136 }
137
138 iter.next() ;
139 }
140}
141
142}
Note: See TracBrowser for help on using the repository browser.