Commits

Pam Harris authored and Ville Suoranta committed 4a3b5d6e8cb Merge
Pull request #564: CAS-13450 add plotms caltable axes

Merge in CASA/casa6 from CAS-13450 to master * commit '7e55aa32f8daae5eb85f4efdc7dfaec6a6b799a8': Update plotms master wheel version Update build.conf Revert code style changes Update plotms branch wheel with test_plotms fix Update plotms branch wheel in build.conf Add errors for invalid plotms caltable axes Add uvw and azel related axes to CTIter Set plotms branch wheel in build.conf Fix new axes for averaged cal table Calculate uvw for cal tables for plotms axes Add new caltable axes az, el, ha, pa
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