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