Line data Source code
1 : //# CTIter.cc: Implementation of CTIter.h
2 : //# Copyright (C) 2011
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //----------------------------------------------------------------------------
27 :
28 : #include <casacore/measures/Measures/Muvw.h>
29 : #include <casacore/measures/Measures/MCBaseline.h>
30 : #include <msvis/MSVis/ViImplementation2.h>
31 : #include <synthesis/CalTables/CTIter.h>
32 : #include <casacore/tables/Tables/ScalarColumn.h>
33 : #include <casacore/casa/Arrays.h>
34 : #include <casacore/casa/OS/Timer.h>
35 :
36 : using namespace casacore;
37 : namespace casa { //# NAMESPACE CASA - BEGIN
38 :
39 : //----------------------------------------------------------------------------
40 :
41 18 : ROCTIter::ROCTIter(NewCalTable tab, const Block<String>& sortcol) :
42 18 : sortCols_(sortcol.begin( ),sortcol.end( )),
43 18 : singleSpw_(false),
44 18 : parentNCT_(tab),
45 18 : calCol_(NULL),
46 18 : ti_(NULL),
47 18 : inct_(NULL),
48 18 : iROCTMainCols_(NULL),
49 18 : init_epoch_(false),
50 36 : init_uvw_(false)
51 : {
52 18 : calCol_=new ROCTColumns(tab);
53 18 : ti_=new TableIterator(tab,sortcol);
54 : // Attach initial accessors:
55 18 : attach();
56 :
57 : // If SPW a sort column, then
58 18 : singleSpw_=anyEQ(sortCols_,String("SPECTRAL_WINDOW_ID"));
59 :
60 : /*
61 : cout << "singleSpw_ = " << boolalpha << singleSpw_ << endl;
62 : cout << "inct_->nrow() = " << inct_->nrow() << endl;
63 : cout << "this->nrow() = " << this->nrow() << endl;
64 : cout << "iROCTMainCols_->spwId() = " << iROCTMainCols_->spwId().getColumn() << endl;
65 : cout << "iROCTMainCols_->spwId()(0) = " << iROCTMainCols_->spwId()(0) << endl;
66 : cout << "thisSpw() = " << this->thisSpw() << endl;
67 :
68 : cout << "done." << endl;
69 : */
70 18 : };
71 :
72 : //----------------------------------------------------------------------------
73 :
74 18 : ROCTIter::~ROCTIter()
75 : {
76 18 : if (calCol_!=NULL) delete calCol_;
77 18 : if (ti_!=NULL) delete ti_;
78 18 : if (iROCTMainCols_!=NULL) delete iROCTMainCols_;
79 18 : if (inct_!=NULL) delete inct_;
80 18 : };
81 :
82 549 : void ROCTIter::next() {
83 : // Advance the TableIterator
84 549 : ti_->next();
85 :
86 : // attach accessors to new iteration
87 549 : this->attach();
88 549 : };
89 :
90 0 : void ROCTIter::next0() {
91 : // Advance the TableIterator
92 0 : ti_->next();
93 0 : };
94 :
95 0 : void ROCTIter::setCTColumns(const NewCalTable& tab) {
96 : // Set subtable columns from another table
97 0 : delete calCol_;
98 0 : calCol_ = new ROCTColumns(tab);
99 0 : }
100 :
101 0 : Double ROCTIter::thisTime() const { return iROCTMainCols_->time()(0); };
102 0 : Vector<Double> ROCTIter::time() const { return iROCTMainCols_->time().getColumn(); };
103 0 : void ROCTIter::time(Vector<Double>& v) const { iROCTMainCols_->time().getColumn(v); };
104 :
105 0 : Int ROCTIter::thisField() const { return iROCTMainCols_->fieldId()(0); };
106 0 : Vector<Int> ROCTIter::field() const { return iROCTMainCols_->fieldId().getColumn(); };
107 0 : void ROCTIter::field(Vector<Int>& v) const { iROCTMainCols_->fieldId().getColumn(v); };
108 :
109 409 : Int ROCTIter::thisSpw() const { return iROCTMainCols_->spwId()(0); };
110 0 : Vector<Int> ROCTIter::spw() const { return iROCTMainCols_->spwId().getColumn(); };
111 0 : void ROCTIter::spw(Vector<Int>& v) const { iROCTMainCols_->spwId().getColumn(v); };
112 :
113 0 : Int ROCTIter::thisScan() const { return iROCTMainCols_->scanNo()(0); };
114 120 : Vector<Int> ROCTIter::scan() const { return iROCTMainCols_->scanNo().getColumn(); };
115 0 : void ROCTIter::scan(Vector<Int>& v) const { iROCTMainCols_->scanNo().getColumn(v); };
116 :
117 120 : Int ROCTIter::thisObs() const { return iROCTMainCols_->obsId()(0); };
118 120 : Vector<Int> ROCTIter::obs() const { return iROCTMainCols_->obsId().getColumn(); };
119 0 : void ROCTIter::obs(Vector<Int>& v) const { iROCTMainCols_->obsId().getColumn(v); };
120 :
121 409 : Int ROCTIter::thisAntenna1() const { return iROCTMainCols_->antenna1()(0); };
122 0 : Vector<Int> ROCTIter::antenna1() const { return iROCTMainCols_->antenna1().getColumn(); };
123 0 : void ROCTIter::antenna1(Vector<Int>& v) const { iROCTMainCols_->antenna1().getColumn(v); };
124 0 : Int ROCTIter::thisAntenna2() const { return iROCTMainCols_->antenna2()(0); };
125 0 : Vector<Int> ROCTIter::antenna2() const { return iROCTMainCols_->antenna2().getColumn(); };
126 0 : void ROCTIter::antenna2(Vector<Int>& v) const { iROCTMainCols_->antenna2().getColumn(v); };
127 :
128 0 : Double ROCTIter::thisInterval() const { return iROCTMainCols_->interval()(0); };
129 0 : Vector<Double> ROCTIter::interval() const { return iROCTMainCols_->interval().getColumn(); };
130 0 : void ROCTIter::interval(Vector<Double>& v) const { iROCTMainCols_->interval().getColumn(v); };
131 :
132 0 : Cube<Complex> ROCTIter::cparam() const { return iROCTMainCols_->cparam().getColumn(); };
133 0 : void ROCTIter::cparam(Cube<Complex>& c) const { iROCTMainCols_->cparam().getColumn(c); };
134 0 : Cube<Float> ROCTIter::fparam() const { return iROCTMainCols_->fparam().getColumn(); };
135 0 : void ROCTIter::fparam(Cube<Float>& f) const { iROCTMainCols_->fparam().getColumn(f); };
136 :
137 : // Complex as Float
138 0 : Cube<Float> ROCTIter::casfparam(String what) const { return iROCTMainCols_->fparamArray(what); };
139 0 : void ROCTIter::casfparam(Cube<Float>& f,String what) const { iROCTMainCols_->fparamArray(f,what); };
140 :
141 0 : Cube<Float> ROCTIter::paramErr() const { return iROCTMainCols_->paramerr().getColumn(); };
142 0 : void ROCTIter::paramErr(Cube<Float>& c) const { iROCTMainCols_->paramerr().getColumn(c); };
143 :
144 0 : Cube<Float> ROCTIter::snr() const { return iROCTMainCols_->snr().getColumn(); };
145 0 : void ROCTIter::snr(Cube<Float>& c) const { iROCTMainCols_->snr().getColumn(c); };
146 0 : Cube<Float> ROCTIter::wt() const { return iROCTMainCols_->weight().getColumn(); };
147 0 : void ROCTIter::wt(Cube<Float>& c) const { iROCTMainCols_->weight().getColumn(c); };
148 :
149 240 : Cube<Bool> ROCTIter::flag() const { return iROCTMainCols_->flag().getColumn(); };
150 0 : void ROCTIter::flag(Cube<Bool>& c) const { iROCTMainCols_->flag().getColumn(c); };
151 :
152 0 : Vector<Int> ROCTIter::chan() const {
153 0 : Vector<Int> chans;
154 0 : this->chan(chans);
155 0 : return chans;
156 0 : }
157 :
158 0 : Int ROCTIter::nchan() const {
159 0 : if (singleSpw_)
160 0 : return calCol_->spectralWindow().numChan()(this->thisSpw());
161 : else
162 : // more than one spw per iteration...
163 0 : throw(AipsError("Please sort by spw."));
164 : }
165 :
166 0 : void ROCTIter::chan(Vector<Int>& v) const {
167 0 : if (singleSpw_) {
168 0 : v.resize(calCol_->spectralWindow().numChan()(this->thisSpw()));
169 : // TBD: observe channel selection here:
170 0 : indgen(v);
171 : }
172 : else
173 : // more than one spw per iteration...
174 0 : throw(AipsError("Please sort by spw."));
175 0 : }
176 :
177 0 : Vector<Double> ROCTIter::freq() const {
178 0 : Vector<Double> freqs;
179 0 : this->freq(freqs);
180 0 : return freqs;
181 0 : }
182 :
183 0 : void ROCTIter::freq(Vector<Double>& v) const {
184 0 : if (singleSpw_) {
185 0 : v.resize();
186 0 : calCol_->spectralWindow().chanFreq().get(this->thisSpw(),v);
187 : }
188 : else
189 : // more than one spw per iteration...
190 0 : throw(AipsError("Please sort by spw."));
191 0 : }
192 :
193 0 : int ROCTIter::freqFrame(int spwId) const {
194 0 : int frame = calCol_->spectralWindow().measFreqRef()(spwId);
195 0 : return frame;
196 : }
197 :
198 570 : void ROCTIter::attach() {
199 : // Attach accessors for current iteration:
200 570 : if (iROCTMainCols_!=NULL) delete iROCTMainCols_;
201 570 : if (inct_!=NULL) delete inct_;
202 :
203 570 : inct_= new NewCalTable(ti_->table());
204 570 : iROCTMainCols_ = new ROCTMainColumns(*inct_);
205 570 : }
206 :
207 0 : casacore::MDirection ROCTIter::azel0(casacore::Double time) {
208 0 : if (!init_epoch_) {
209 0 : initEpoch();
210 : }
211 :
212 : try {
213 0 : updatePhaseCenter();
214 0 : } catch (const casacore::AipsError& err) {
215 0 : throw(casacore::AipsError("azel: " + err.getMesg()));
216 0 : }
217 :
218 0 : casacore::MSDerivedValues msd;
219 0 : msd.setAntennas(calCol_->antenna());
220 0 : msd.setFieldCenter(phaseCenter_);
221 :
222 0 : casacore::MDirection azel;
223 0 : vi::ViImplementation2::azel0Calculate(time, msd, azel, epoch_);
224 0 : return azel;
225 0 : }
226 :
227 0 : casacore::Double ROCTIter::hourang(casacore::Double time) {
228 0 : if (!init_epoch_) {
229 0 : initEpoch();
230 : }
231 :
232 : try {
233 0 : updatePhaseCenter();
234 0 : } catch (const casacore::AipsError& err) {
235 0 : throw(casacore::AipsError("hourang: " + err.getMesg()));
236 0 : }
237 :
238 0 : casacore::MSDerivedValues msd;
239 0 : msd.setAntennas(calCol_->antenna());
240 0 : msd.setFieldCenter(phaseCenter_);
241 :
242 0 : return vi::ViImplementation2::hourangCalculate(time, msd, epoch_);
243 0 : }
244 :
245 0 : casacore::Float ROCTIter::parang0(casacore::Double time) {
246 0 : if (!init_epoch_) {
247 0 : initEpoch();
248 : }
249 :
250 : try {
251 0 : updatePhaseCenter();
252 0 : } catch (const casacore::AipsError& err) {
253 0 : throw(casacore::AipsError("parang: " + err.getMesg()));
254 0 : }
255 :
256 0 : casacore::MSDerivedValues msd;
257 0 : msd.setAntennas(calCol_->antenna());
258 0 : msd.setFieldCenter(phaseCenter_);
259 :
260 0 : return vi::ViImplementation2::parang0Calculate(time, msd, epoch_);
261 0 : }
262 :
263 0 : casacore::Matrix<casacore::Double> ROCTIter::uvw() {
264 0 : casacore::Vector<casacore::Int> ant2 = antenna2();
265 0 : if (ant2(0) == -1) {
266 0 : throw(AipsError("UVW axes cannot be calculated for antenna-based calibration tables."));
267 : }
268 :
269 0 : updateAntennaUVW();
270 :
271 0 : casacore::Vector<casacore::Int> ant1 = antenna1();
272 0 : auto nbaseline = ant1.size();
273 0 : casacore::Matrix<casacore::Double> uvw(3, nbaseline);
274 :
275 0 : for (uInt i = 0; i < nbaseline; ++i) {
276 0 : uvw.column(i) = antennaUVW_[ant2(i)] - antennaUVW_[ant1(i)];
277 : }
278 :
279 0 : return uvw;
280 0 : }
281 :
282 0 : void ROCTIter::updatePhaseCenter() {
283 : // Set MSDerivedValues phase center for field
284 0 : if ((thisTime() == 0) || (thisField() < 0)) {
285 0 : throw(AipsError("cannot calculate this value with no valid timestamp or field ID."));
286 : }
287 :
288 0 : phaseCenter_ = calCol_->field().phaseDirMeas(thisField(), thisTime());
289 0 : }
290 :
291 0 : void ROCTIter::initEpoch() {
292 0 : epoch_ = iROCTMainCols_->timeMeas()(0);
293 0 : init_epoch_ = true;
294 0 : }
295 :
296 0 : void ROCTIter::initUVW() {
297 : // Calculate relative positions of antennas
298 0 : nAnt_ = calCol_->antenna().nrow();
299 0 : auto antPosMeas = calCol_->antenna().positionMeas();
300 0 : refAntPos_ = antPosMeas(0); // use first antenna for reference
301 0 : auto refAntPosValue = refAntPos_.getValue();
302 :
303 : // Set up baseline and uvw types
304 0 : MBaseline::getType(baseline_type_, MPosition::showType(refAntPos_.getRef().getType()));
305 :
306 0 : mvbaselines_.resize(nAnt_);
307 :
308 0 : for (Int ant = 0; ant < nAnt_; ++ant) {
309 : // MVBaselines are basically xyz Vectors, not Measures
310 0 : mvbaselines_[ant] = MVBaseline(refAntPosValue, antPosMeas(ant).getValue());
311 : }
312 :
313 0 : init_uvw_ = true;
314 0 : }
315 :
316 0 : void ROCTIter::updateAntennaUVW() {
317 : // Set antennaUVW_ for current iteration
318 0 : if (!init_uvw_) {
319 0 : initUVW();
320 : }
321 :
322 0 : updatePhaseCenter();
323 :
324 0 : MeasFrame measFrame(refAntPos_, epoch_, phaseCenter_);
325 :
326 : // Antenna frame
327 0 : MBaseline::Ref baselineRef(baseline_type_);
328 0 : MVBaseline mvbaseline;
329 0 : MBaseline baselineMeas;
330 0 : baselineMeas.set(mvbaseline, baselineRef);
331 0 : baselineMeas.getRefPtr()->set(measFrame);
332 :
333 : // Conversion engine to phasedir type
334 : casacore::MBaseline::Types phasedir_type;
335 0 : MBaseline::getType(phasedir_type, MDirection::showType(phaseCenter_.getRef().getType()));
336 0 : MBaseline::Ref uvwRef(phasedir_type);
337 0 : MBaseline::Convert baselineConv(baselineMeas, uvwRef);
338 :
339 : // WSRT convention: phase opposite to VLA (l increases toward increasing RA)
340 0 : Int obsId = thisObs();
341 0 : if (obsId < 0) obsId = 0;
342 0 : casacore::String telescope = calCol_->observation().telescopeName()(obsId);
343 0 : bool wsrtConvention = (telescope == "WSRT");
344 :
345 0 : antennaUVW_.resize(nAnt_);
346 :
347 0 : for (int i = 0; i < nAnt_; ++i) {
348 0 : baselineMeas.set(mvbaselines_[i], baselineRef);
349 0 : MBaseline baselineOutFrame = baselineConv(baselineMeas);
350 0 : MVuvw uvwOutFrame(baselineOutFrame.getValue(), phaseCenter_.getValue());
351 :
352 0 : if (wsrtConvention) {
353 0 : antennaUVW_[i] = -uvwOutFrame.getValue();
354 : } else {
355 0 : antennaUVW_[i] = uvwOutFrame.getValue();
356 : }
357 0 : }
358 0 : }
359 :
360 :
361 : // CTIter
362 :
363 1 : CTIter::CTIter(NewCalTable tab, const Block<String>& sortcol) :
364 : ROCTIter(tab,sortcol),
365 1 : irwnct_(NULL),
366 1 : iRWCTMainCols_(NULL)
367 : {
368 : // Attach first iteration
369 : // TBD: this unnecessarily redoes the ROCTIter attach...
370 1 : attach();
371 1 : }
372 :
373 1 : CTIter::~CTIter() {
374 1 : if (iRWCTMainCols_!=NULL) delete iRWCTMainCols_;
375 1 : if (irwnct_!=NULL) delete irwnct_;
376 1 : }
377 :
378 : // Set fieldid
379 0 : void CTIter::setfield(Int fieldid) {
380 0 : iRWCTMainCols_->fieldId().fillColumn(fieldid);
381 0 : }
382 :
383 : // Set scan number
384 0 : void CTIter::setscan(Int scan) {
385 0 : iRWCTMainCols_->scanNo().fillColumn(scan);
386 0 : }
387 :
388 : // Set obsid
389 0 : void CTIter::setobs(Int obs) {
390 0 : iRWCTMainCols_->obsId().fillColumn(obs);
391 0 : }
392 :
393 : // Set antenna2 (e.g., used for setting refant)
394 0 : void CTIter::setantenna2(const Vector<Int>& a2) {
395 0 : iRWCTMainCols_->antenna2().putColumn(a2);
396 0 : }
397 :
398 120 : void CTIter::setflag(const Cube<Bool>& fl) {
399 120 : iRWCTMainCols_->flag().putColumn(fl);
400 120 : }
401 :
402 0 : void CTIter::setfparam(const Cube<Float>& f) {
403 0 : iRWCTMainCols_->fparam().putColumn(f);
404 0 : };
405 :
406 0 : void CTIter::setcparam(const Cube<Complex>& c) {
407 0 : iRWCTMainCols_->cparam().putColumn(c);
408 0 : };
409 :
410 242 : void CTIter::attach() {
411 :
412 : // Attach readonly access
413 242 : ROCTIter::attach();
414 :
415 : // Attach writable access
416 242 : if (iRWCTMainCols_!=NULL) delete iRWCTMainCols_;
417 242 : if (irwnct_!=NULL) delete irwnct_;
418 242 : irwnct_= new NewCalTable(this->table());
419 242 : iRWCTMainCols_ = new CTMainColumns(*irwnct_);
420 242 : }
421 :
422 :
423 :
424 : } //# NAMESPACE CASA - END
|
