Commits

Pam Harris authored 3c2e662404a Merge
Merge branch 'master' into CAS-13724
No tags

casatools/src/code/synthesis/CalTables/CTIter.cc

Modified
16 16 //# along with this library; if not, write to the Free Software Foundation,
17 17 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 18 //#
19 19 //# Correspondence concerning AIPS++ should be addressed as follows:
20 20 //# Internet email: aips2-request@nrao.edu.
21 21 //# Postal address: AIPS++ Project Office
22 22 //# National Radio Astronomy Observatory
23 23 //# 520 Edgemont Road
24 24 //# Charlottesville, VA 22903-2475 USA
25 25 //#
26 -//# $Id: NewCalTable.cc 15602 2011-07-14 00:03:34Z tak.tsutsumi $
27 26 //----------------------------------------------------------------------------
28 27
28 +#include <measures/Measures/Muvw.h>
29 +#include <measures/Measures/MCBaseline.h>
30 +#include <msvis/MSVis/ViImplementation2.h>
29 31 #include <synthesis/CalTables/CTIter.h>
30 32 #include <tables/Tables/ScalarColumn.h>
31 33 #include <casa/Arrays.h>
32 34 #include <casa/OS/Timer.h>
33 35
34 36 using namespace casacore;
35 37 namespace casa { //# NAMESPACE CASA - BEGIN
36 38
37 39 //----------------------------------------------------------------------------
38 40
39 41 ROCTIter::ROCTIter(NewCalTable tab, const Block<String>& sortcol) :
40 42 sortCols_(sortcol.begin( ),sortcol.end( )),
41 43 singleSpw_(false),
42 44 parentNCT_(tab),
43 - calCol_(tab),
45 + calCol_(NULL),
44 46 ti_(NULL),
45 47 inct_(NULL),
46 - iROCTMainCols_(NULL)
48 + iROCTMainCols_(NULL),
49 + init_epoch_(false),
50 + init_uvw_(false)
47 51 {
48 -
52 + calCol_=new ROCTColumns(tab);
49 53 ti_=new TableIterator(tab,sortcol);
50 -
51 54 // Attach initial accessors:
52 55 attach();
53 56
54 - // cout << "CTIter sort columns (" << sortCols_.nelements() << "): " << sortCols_ << endl;
55 -
56 57 // If SPW a sort column, then
57 58 singleSpw_=anyEQ(sortCols_,String("SPECTRAL_WINDOW_ID"));
58 -
59 +
59 60 /*
60 61 cout << "singleSpw_ = " << boolalpha << singleSpw_ << endl;
61 62 cout << "inct_->nrow() = " << inct_->nrow() << endl;
62 63 cout << "this->nrow() = " << this->nrow() << endl;
63 64 cout << "iROCTMainCols_->spwId() = " << iROCTMainCols_->spwId().getColumn() << endl;
64 65 cout << "iROCTMainCols_->spwId()(0) = " << iROCTMainCols_->spwId()(0) << endl;
65 66 cout << "thisSpw() = " << this->thisSpw() << endl;
66 67
67 68 cout << "done." << endl;
68 69 */
69 70 };
70 71
71 72 //----------------------------------------------------------------------------
72 73
73 74 ROCTIter::~ROCTIter()
74 75 {
76 + if (calCol_!=NULL) delete calCol_;
75 77 if (ti_!=NULL) delete ti_;
76 78 if (iROCTMainCols_!=NULL) delete iROCTMainCols_;
77 79 if (inct_!=NULL) delete inct_;
78 80 };
79 81
80 82 void ROCTIter::next() {
81 83 // Advance the TableIterator
82 84 ti_->next();
83 85
84 86 // attach accessors to new iteration
85 87 this->attach();
86 -
87 88 };
88 89
89 90 void ROCTIter::next0() {
90 91 // Advance the TableIterator
91 92 ti_->next();
92 93 };
93 94
95 +void ROCTIter::setCTColumns(const NewCalTable& tab) {
96 + // Set subtable columns from another table
97 + delete calCol_;
98 + calCol_ = new ROCTColumns(tab);
99 +}
100 +
94 101 Double ROCTIter::thisTime() const { return iROCTMainCols_->time()(0); };
95 102 Vector<Double> ROCTIter::time() const { return iROCTMainCols_->time().getColumn(); };
96 103 void ROCTIter::time(Vector<Double>& v) const { iROCTMainCols_->time().getColumn(v); };
97 104
98 105 Int ROCTIter::thisField() const { return iROCTMainCols_->fieldId()(0); };
99 106 Vector<Int> ROCTIter::field() const { return iROCTMainCols_->fieldId().getColumn(); };
100 107 void ROCTIter::field(Vector<Int>& v) const { iROCTMainCols_->fieldId().getColumn(v); };
101 108
102 109 Int ROCTIter::thisSpw() const { return iROCTMainCols_->spwId()(0); };
103 110 Vector<Int> ROCTIter::spw() const { return iROCTMainCols_->spwId().getColumn(); };
143 150 void ROCTIter::flag(Cube<Bool>& c) const { iROCTMainCols_->flag().getColumn(c); };
144 151
145 152 Vector<Int> ROCTIter::chan() const {
146 153 Vector<Int> chans;
147 154 this->chan(chans);
148 155 return chans;
149 156 }
150 157
151 158 Int ROCTIter::nchan() const {
152 159 if (singleSpw_)
153 - return calCol_.spectralWindow().numChan()(this->thisSpw());
160 + return calCol_->spectralWindow().numChan()(this->thisSpw());
154 161 else
155 162 // more than one spw per iteration...
156 163 throw(AipsError("Please sort by spw."));
157 164 }
158 165
159 166 void ROCTIter::chan(Vector<Int>& v) const {
160 167 if (singleSpw_) {
161 - v.resize(calCol_.spectralWindow().numChan()(this->thisSpw()));
168 + v.resize(calCol_->spectralWindow().numChan()(this->thisSpw()));
162 169 // TBD: observe channel selection here:
163 170 indgen(v);
164 171 }
165 172 else
166 173 // more than one spw per iteration...
167 174 throw(AipsError("Please sort by spw."));
168 175 }
169 176
170 177 Vector<Double> ROCTIter::freq() const {
171 178 Vector<Double> freqs;
172 179 this->freq(freqs);
173 180 return freqs;
174 181 }
175 182
176 183 void ROCTIter::freq(Vector<Double>& v) const {
177 184 if (singleSpw_) {
178 185 v.resize();
179 - calCol_.spectralWindow().chanFreq().get(this->thisSpw(),v);
186 + calCol_->spectralWindow().chanFreq().get(this->thisSpw(),v);
180 187 }
181 188 else
182 189 // more than one spw per iteration...
183 190 throw(AipsError("Please sort by spw."));
184 191 }
185 192
186 193 int ROCTIter::freqFrame(int spwId) const {
187 - int frame = calCol_.spectralWindow().measFreqRef()(spwId);
194 + int frame = calCol_->spectralWindow().measFreqRef()(spwId);
188 195 return frame;
189 196 }
190 197
191 198 void ROCTIter::attach() {
192 - // Attach accessors:
199 + // Attach accessors for current iteration:
193 200 if (iROCTMainCols_!=NULL) delete iROCTMainCols_;
194 201 if (inct_!=NULL) delete inct_;
202 +
195 203 inct_= new NewCalTable(ti_->table());
196 204 iROCTMainCols_ = new ROCTMainColumns(*inct_);
197 205 }
198 206
207 +casacore::MDirection ROCTIter::azel0(casacore::Double time) {
208 + if (!init_epoch_) {
209 + initEpoch();
210 + }
211 +
212 + try {
213 + updatePhaseCenter();
214 + } catch (const casacore::AipsError& err) {
215 + throw(casacore::AipsError("azel: " + err.getMesg()));
216 + }
217 +
218 + casacore::MSDerivedValues msd;
219 + msd.setAntennas(calCol_->antenna());
220 + msd.setFieldCenter(phaseCenter_);
221 +
222 + casacore::MDirection azel;
223 + vi::ViImplementation2::azel0Calculate(time, msd, azel, epoch_);
224 + return azel;
225 +}
226 +
227 +casacore::Double ROCTIter::hourang(casacore::Double time) {
228 + if (!init_epoch_) {
229 + initEpoch();
230 + }
231 +
232 + try {
233 + updatePhaseCenter();
234 + } catch (const casacore::AipsError& err) {
235 + throw(casacore::AipsError("hourang: " + err.getMesg()));
236 + }
237 +
238 + casacore::MSDerivedValues msd;
239 + msd.setAntennas(calCol_->antenna());
240 + msd.setFieldCenter(phaseCenter_);
241 +
242 + return vi::ViImplementation2::hourangCalculate(time, msd, epoch_);
243 +}
244 +
245 +casacore::Float ROCTIter::parang0(casacore::Double time) {
246 + if (!init_epoch_) {
247 + initEpoch();
248 + }
249 +
250 + try {
251 + updatePhaseCenter();
252 + } catch (const casacore::AipsError& err) {
253 + throw(casacore::AipsError("parang: " + err.getMesg()));
254 + }
255 +
256 + casacore::MSDerivedValues msd;
257 + msd.setAntennas(calCol_->antenna());
258 + msd.setFieldCenter(phaseCenter_);
259 +
260 + return vi::ViImplementation2::parang0Calculate(time, msd, epoch_);
261 +}
262 +
263 +casacore::Matrix<casacore::Double> ROCTIter::uvw() {
264 + casacore::Vector<casacore::Int> ant2 = antenna2();
265 + if (ant2(0) == -1) {
266 + throw(AipsError("UVW axes cannot be calculated for antenna-based calibration tables."));
267 + }
268 +
269 + updateAntennaUVW();
270 +
271 + casacore::Vector<casacore::Int> ant1 = antenna1();
272 + auto nbaseline = ant1.size();
273 + casacore::Matrix<casacore::Double> uvw(3, nbaseline);
274 +
275 + for (uInt i = 0; i < nbaseline; ++i) {
276 + uvw.column(i) = antennaUVW_[ant2(i)] - antennaUVW_[ant1(i)];
277 + }
278 +
279 + return uvw;
280 +}
281 +
282 +void ROCTIter::updatePhaseCenter() {
283 + // Set MSDerivedValues phase center for field
284 + if ((thisTime() == 0) || (thisField() < 0)) {
285 + throw(AipsError("cannot calculate this value with no valid timestamp or field ID."));
286 + }
287 +
288 + phaseCenter_ = calCol_->field().phaseDirMeas(thisField(), thisTime());
289 +}
290 +
291 +void ROCTIter::initEpoch() {
292 + epoch_ = iROCTMainCols_->timeMeas()(0);
293 + init_epoch_ = true;
294 +}
295 +
296 +void ROCTIter::initUVW() {
297 + // Calculate relative positions of antennas
298 + nAnt_ = calCol_->antenna().nrow();
299 + auto antPosMeas = calCol_->antenna().positionMeas();
300 + refAntPos_ = antPosMeas(0); // use first antenna for reference
301 + auto refAntPosValue = refAntPos_.getValue();
302 +
303 + // Set up baseline and uvw types
304 + MBaseline::getType(baseline_type_, MPosition::showType(refAntPos_.getRef().getType()));
305 +
306 + mvbaselines_.resize(nAnt_);
307 +
308 + for (Int ant = 0; ant < nAnt_; ++ant) {
309 + // MVBaselines are basically xyz Vectors, not Measures
310 + mvbaselines_[ant] = MVBaseline(refAntPosValue, antPosMeas(ant).getValue());
311 + }
312 +
313 + init_uvw_ = true;
314 +}
315 +
316 +void ROCTIter::updateAntennaUVW() {
317 + // Set antennaUVW_ for current iteration
318 + if (!init_uvw_) {
319 + initUVW();
320 + }
321 +
322 + updatePhaseCenter();
323 +
324 + MeasFrame measFrame(refAntPos_, epoch_, phaseCenter_);
325 +
326 + // Antenna frame
327 + MBaseline::Ref baselineRef(baseline_type_);
328 + MVBaseline mvbaseline;
329 + MBaseline baselineMeas;
330 + baselineMeas.set(mvbaseline, baselineRef);
331 + baselineMeas.getRefPtr()->set(measFrame);
332 +
333 + // Conversion engine to phasedir type
334 + casacore::MBaseline::Types phasedir_type;
335 + MBaseline::getType(phasedir_type, MDirection::showType(phaseCenter_.getRef().getType()));
336 + MBaseline::Ref uvwRef(phasedir_type);
337 + MBaseline::Convert baselineConv(baselineMeas, uvwRef);
338 +
339 + // WSRT convention: phase opposite to VLA (l increases toward increasing RA)
340 + Int obsId = thisObs();
341 + if (obsId < 0) obsId = 0;
342 + casacore::String telescope = calCol_->observation().telescopeName()(obsId);
343 + bool wsrtConvention = (telescope == "WSRT");
344 +
345 + antennaUVW_.resize(nAnt_);
346 +
347 + for (int i = 0; i < nAnt_; ++i) {
348 + baselineMeas.set(mvbaselines_[i], baselineRef);
349 + MBaseline baselineOutFrame = baselineConv(baselineMeas);
350 + MVuvw uvwOutFrame(baselineOutFrame.getValue(), phaseCenter_.getValue());
351 +
352 + if (wsrtConvention) {
353 + antennaUVW_[i] = -uvwOutFrame.getValue();
354 + } else {
355 + antennaUVW_[i] = uvwOutFrame.getValue();
356 + }
357 + }
358 +}
359 +
360 +
361 +// CTIter
199 362
200 363 CTIter::CTIter(NewCalTable tab, const Block<String>& sortcol) :
201 364 ROCTIter(tab,sortcol),
202 365 irwnct_(NULL),
203 366 iRWCTMainCols_(NULL)
204 367 {
205 368 // Attach first iteration
206 369 // TBD: this unnecessarily redoes the ROCTIter attach...
207 370 attach();
208 371 }
209 372
210 373 CTIter::~CTIter() {
211 374 if (iRWCTMainCols_!=NULL) delete iRWCTMainCols_;
212 375 if (irwnct_!=NULL) delete irwnct_;
213 376 }
214 377
215 -
216 378 // Set fieldid
217 379 void CTIter::setfield(Int fieldid) {
218 380 iRWCTMainCols_->fieldId().fillColumn(fieldid);
219 381 }
220 382
221 383 // Set scan number
222 384 void CTIter::setscan(Int scan) {
223 385 iRWCTMainCols_->scanNo().fillColumn(scan);
224 386 }
225 387

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut