Line data Source code
1 : //# MSUtil.cc: Some MS specific Utilities
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 adressed 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 : //# $Id$
28 : #include <casacore/casa/Utilities/Sort.h>
29 : #include <casacore/measures/Measures/MeasTable.h>
30 : #include <casacore/ms/MeasurementSets/MSColumns.h>
31 : #include <casacore/ms/MSSel/MSSpwIndex.h>
32 : #include <casacore/ms/MSSel/MSDataDescIndex.h>
33 : #include <msvis/MSVis/MSUtil.h>
34 : #include <casacore/casa/Arrays/ArrayMath.h>
35 : #include <casacore/casa/OS/Path.h>
36 : #include <casacore/casa/Utilities/GenSort.h>
37 : #include <iomanip>
38 : using namespace casacore;
39 : namespace casa { //# NAMESPACE CASA - BEGIN
40 :
41 0 : MSUtil::MSUtil(){};
42 7 : void MSUtil::getSpwInFreqRange(Vector<Int>& spw, Vector<Int>& start,
43 : Vector<Int>& nchan,
44 : const MeasurementSet& ms,
45 : const Double freqStart,
46 : const Double freqEnd,
47 : const Double freqStep,
48 : const MFrequency::Types freqframe,
49 : const Int fieldId)
50 : {
51 7 : spw.resize();
52 7 : start.resize();
53 7 : nchan.resize();
54 7 : Vector<Double> t;
55 7 : ScalarColumn<Double> (ms,MS::columnName(MS::TIME)).getColumn(t);
56 : //Vector<Int> ddId;
57 : //Vector<Int> fldId;
58 :
59 7 : MSFieldColumns fieldCol(ms.field());
60 7 : MSDataDescColumns ddCol(ms.dataDescription());
61 7 : MSSpWindowColumns spwCol(ms.spectralWindow());
62 7 : ROScalarMeasColumn<MEpoch> timeCol(ms, MS::columnName(MS::TIME));
63 7 : Vector<uInt> uniqIndx;
64 7 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, t, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
65 :
66 7 : t.resize(0);
67 : //ScalarColumn<Int> (ms,MS::columnName(MS::DATA_DESC_ID)).getColumn(ddId);
68 : //ScalarColumn<Int> (ms,MS::columnName(MS::FIELD_ID)).getColumn(fldId);
69 7 : ScalarColumn<Int> ddId(ms,MS::columnName(MS::DATA_DESC_ID));
70 7 : ScalarColumn<Int> fldId(ms,MS::columnName(MS::FIELD_ID));
71 : //now need to do the conversion to data frame from requested frame
72 : //Get the epoch mesasures of the first row
73 7 : MEpoch ep;
74 7 : timeCol.get(0, ep);
75 7 : String observatory;
76 7 : MPosition obsPos;
77 : /////observatory position
78 7 : MSColumns msc(ms);
79 7 : if (ms.observation().nrow() > 0) {
80 7 : observatory = msc.observation().telescopeName()(msc.observationId()(0));
81 : }
82 14 : if (observatory.length() == 0 ||
83 7 : !MeasTable::Observatory(obsPos,observatory)) {
84 : // unknown observatory, use first antenna
85 0 : obsPos=msc.antenna().positionMeas()(0);
86 : }
87 : //////
88 7 : Int oldDD=ddId(0);
89 7 : Int newDD=oldDD;
90 : //For now we will assume that the field is not moving very far from polynome 0
91 7 : MDirection dir =fieldCol.phaseDirMeas(fieldId);
92 7 : MFrequency::Types obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(ddCol.spectralWindowId()(ddId(0))));
93 : //cout << "nTimes " << nTimes << endl;
94 : //cout << " obsframe " << obsMFreqType << " reqFrame " << freqframe << endl;
95 7 : MeasFrame frame(ep, obsPos, dir);
96 : MFrequency::Convert toObs(freqframe,
97 7 : MFrequency::Ref(obsMFreqType, frame));
98 7 : Double freqEndMax=freqEnd;
99 7 : Double freqStartMin=freqStart;
100 7 : if(freqframe != obsMFreqType){
101 7 : freqEndMax=0.0;
102 7 : freqStartMin=C::dbl_max;
103 : }
104 3367 : for (uInt j=0; j< nTimes; ++j){
105 3360 : if(fldId(uniqIndx[j]) ==fieldId){
106 3360 : timeCol.get(uniqIndx[j], ep);
107 3360 : newDD=ddId(uniqIndx[j]);
108 3360 : if(oldDD != newDD){
109 0 : oldDD=newDD;
110 0 : if(spwCol.measFreqRef()(ddCol.spectralWindowId()(newDD)) != obsMFreqType){
111 0 : obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(ddCol.spectralWindowId()(newDD)));
112 0 : toObs.setOut(MFrequency::Ref(obsMFreqType, frame));
113 : }
114 : }
115 3360 : if(obsMFreqType != freqframe){
116 3360 : frame.resetEpoch(ep);
117 3360 : Double freqTmp=toObs(Quantity(freqStart, "Hz")).get("Hz").getValue();
118 3360 : freqStartMin=(freqStartMin > freqTmp) ? freqTmp : freqStartMin;
119 3360 : freqTmp=toObs(Quantity(freqEnd, "Hz")).get("Hz").getValue();
120 3360 : freqEndMax=(freqEndMax < freqTmp) ? freqTmp : freqEndMax;
121 : }
122 : }
123 : }
124 :
125 : //cout << "freqStartMin " << freqStartMin << " freqEndMax " << freqEndMax << endl;
126 7 : MSSpwIndex spwIn(ms.spectralWindow());
127 7 : spwIn.matchFrequencyRange(freqStartMin-0.5*freqStep, freqEndMax+0.5*freqStep, spw, start, nchan);
128 :
129 :
130 :
131 7 : }
132 :
133 0 : void MSUtil::getSpwInSourceFreqRange(Vector<Int>& spw, Vector<Int>& start,
134 : Vector<Int>& nchan,
135 : const MeasurementSet& ms,
136 : const Double freqStart,
137 : const Double freqEnd,
138 : const Double freqStep,
139 : const String& ephemtab,
140 : const Int fieldId)
141 : {
142 0 : spw.resize();
143 0 : start.resize();
144 0 : nchan.resize();
145 0 : Vector<Double> t;
146 0 : ScalarColumn<Double> (ms,MS::columnName(MS::TIME)).getColumn(t);
147 : //Vector<Int> ddId;
148 : //Vector<Int> fldId;
149 :
150 0 : MSFieldColumns fieldCol(ms.field());
151 0 : MSDataDescColumns ddCol(ms.dataDescription());
152 0 : MSSpWindowColumns spwCol(ms.spectralWindow());
153 0 : ROScalarMeasColumn<MEpoch> timeCol(ms, MS::columnName(MS::TIME));
154 0 : Vector<uInt> uniqIndx;
155 0 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, t, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
156 :
157 0 : t.resize(0);
158 : //ScalarColumn<Int> (ms,MS::columnName(MS::DATA_DESC_ID)).getColumn(ddId);
159 : //ScalarColumn<Int> (ms,MS::columnName(MS::FIELD_ID)).getColumn(fldId);
160 0 : ScalarColumn<Int> ddId(ms,MS::columnName(MS::DATA_DESC_ID));
161 0 : ScalarColumn<Int> fldId(ms,MS::columnName(MS::FIELD_ID));
162 : //now need to do the conversion to data frame from requested frame
163 : //Get the epoch mesasures of the first row
164 0 : MEpoch ep;
165 0 : timeCol.get(0, ep);
166 0 : String observatory;
167 0 : MPosition obsPos;
168 : /////observatory position
169 0 : MSColumns msc(ms);
170 0 : if (ms.observation().nrow() > 0) {
171 0 : observatory = msc.observation().telescopeName()(msc.observationId()(0));
172 : }
173 0 : if (observatory.length() == 0 ||
174 0 : !MeasTable::Observatory(obsPos,observatory)) {
175 : // unknown observatory, use first antenna
176 0 : obsPos=msc.antenna().positionMeas()(0);
177 : }
178 : //////
179 : // Int oldDD=ddId(0); // unused/commented out below in the oldDD!=newDD if
180 : // Int newDD=oldDD; // unused/commented out below in the oldDD!=newDD if
181 : //For now we will assume that the field is not moving very far from polynome 0
182 0 : MDirection dir =fieldCol.phaseDirMeas(fieldId);
183 0 : MFrequency::Types obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(ddCol.spectralWindowId()(ddId(0))));
184 0 : if( obsMFreqType != MFrequency::TOPO)
185 0 : throw(AipsError("No dealing with non topo data for moving source yet"));
186 : //cout << "nTimes " << nTimes << endl;
187 : //cout << " obsframe " << obsMFreqType << " reqFrame " << freqframe << endl;
188 0 : MeasFrame frame(ep, obsPos, dir);
189 : // MFrequency::Convert toObs(freqframe,MFrequency::Ref(obsMFreqType, frame);
190 0 : MDoppler toObs;
191 0 : MDoppler toSource;
192 0 : setupSourceObsVelSystem(ephemtab, ms, fieldId, toSource, toObs,frame);
193 0 : Double freqEndMax=0.0;
194 0 : Double freqStartMin=C::dbl_max;
195 :
196 0 : for (uInt j=0; j< nTimes; ++j){
197 0 : if(fldId(uniqIndx[j]) ==fieldId){
198 0 : timeCol.get(uniqIndx[j], ep);
199 : // newDD=ddId(uniqIndx[j]); // unused below
200 : /*if(oldDD != newDD){
201 : oldDD=newDD;
202 : if(spwCol.measFreqRef()(ddCol.spectralWindowId()(newDD)) != obsMFreqType){
203 : obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(ddCol.spectralWindowId()(newDD)));
204 : toObs.setOut(MFrequency::Ref(obsMFreqType, frame));
205 : }
206 : }
207 : */
208 : //if(obsMFreqType != freqframe){
209 0 : frame.resetEpoch(ep);
210 0 : Vector<Double> freqTmp(2);
211 0 : freqTmp[0]=freqStart;
212 0 : freqTmp[1]=freqEnd;
213 0 : Vector<Double> newFreqs=toObs.shiftFrequency(freqTmp);
214 : //Double freqTmp=toObs(Quantity(freqStart, "Hz")).get("Hz").getValue();
215 0 : freqStartMin=(freqStartMin > newFreqs[0]) ? newFreqs[0] : freqStartMin;
216 : //freqTmp=toObs(Quantity(freqEnd, "Hz")).get("Hz").getValue();
217 0 : freqEndMax=(freqEndMax < newFreqs[1]) ? newFreqs[1] : freqEndMax;
218 : //}
219 0 : }
220 : }
221 :
222 : //cout << "freqStartMin " << freqStartMin << " freqEndMax " << freqEndMax << endl;
223 0 : MSSpwIndex spwIn(ms.spectralWindow());
224 0 : spwIn.matchFrequencyRange(freqStartMin-0.5*freqStep, freqEndMax+0.5*freqStep, spw, start, nchan);
225 :
226 :
227 :
228 0 : }
229 0 : void MSUtil:: setupSourceObsVelSystem(const String& ephemTable, const MeasurementSet& ms, const Int& fieldid, MDoppler& toSource, MDoppler& toObs, MeasFrame& mFrame){
230 0 : String ephemtab("");
231 0 : const MSColumns mscol(ms);
232 0 : if(Table::isReadable(ephemTable)){
233 0 : ephemtab=ephemTable;
234 : }
235 0 : else if(ephemTable=="TRACKFIELD"){
236 0 : ephemtab=(mscol.field()).ephemPath(fieldid);
237 :
238 : }
239 0 : MRadialVelocity::Types refvel=MRadialVelocity::GEO;
240 0 : MEpoch ep(mFrame.epoch());
241 0 : if(ephemtab != ""){
242 :
243 0 : MeasComet mcomet(Path(ephemtab).absoluteName());
244 0 : if(mFrame.comet())
245 0 : mFrame.resetComet(mcomet);
246 : else
247 0 : mFrame.set(mcomet);
248 0 : if(mcomet.getTopo().getLength("km").getValue() > 1.0e-3){
249 0 : refvel=MRadialVelocity::TOPO;
250 : }
251 : ////Will use UT for now for ephem tables as it is not clear that they are being
252 : ///filled with TDB as intended in MeasComet.h
253 0 : MEpoch::Convert toUT(ep, MEpoch::UT);
254 0 : MVRadialVelocity cometvel;
255 0 : mcomet.getRadVel(cometvel, toUT(ep).get("d").getValue());
256 0 : MRadialVelocity::Convert obsvelconv(MRadialVelocity(MVRadialVelocity(0.0),
257 0 : MRadialVelocity::Ref(MRadialVelocity::TOPO, mFrame)), MRadialVelocity::Ref(refvel));
258 0 : toSource=MDoppler(Quantity(-cometvel.get("km/s").getValue("km/s")+obsvelconv().get("km/s").getValue("km/s") , "km/s"), MDoppler::RELATIVISTIC);
259 0 : toObs=MDoppler(Quantity(cometvel.get("km/s").getValue("km/s")-obsvelconv().get("km/s").getValue("km/s") , "km/s"), MDoppler::RELATIVISTIC);
260 :
261 :
262 0 : }
263 : else{//Must be a DE-200 canned source that measures know
264 0 : ephemtab=upcase(ephemTable);
265 0 : MeasTable::Types mtype=MeasTable::BARYEARTH;
266 : MDirection::Types planettype;
267 0 : if(!MDirection::getType(planettype, ephemtab))
268 0 : throw(AipsError("Did not understand sourcename as a known solar system object"));
269 0 : switch(planettype){
270 0 : case MDirection::MERCURY :
271 0 : mtype=MeasTable::MERCURY;
272 0 : break;
273 0 : case MDirection::VENUS :
274 0 : mtype=MeasTable::VENUS;
275 0 : break;
276 0 : case MDirection::MARS :
277 0 : mtype=MeasTable::MARS;
278 0 : break;
279 0 : case MDirection::JUPITER :
280 0 : mtype=MeasTable::JUPITER;
281 0 : break;
282 0 : case MDirection::SATURN :
283 0 : mtype=MeasTable::SATURN;
284 0 : break;
285 0 : case MDirection::URANUS :
286 0 : mtype=MeasTable::URANUS;
287 0 : break;
288 0 : case MDirection::NEPTUNE :
289 0 : mtype=MeasTable::NEPTUNE;
290 0 : break;
291 0 : case MDirection::PLUTO :
292 0 : mtype=MeasTable::PLUTO;
293 0 : break;
294 0 : case MDirection::MOON :
295 0 : mtype=MeasTable::MOON;
296 0 : break;
297 0 : case MDirection::SUN :
298 0 : mtype=MeasTable::SUN;
299 0 : break;
300 0 : default:
301 0 : throw(AipsError("Cannot translate to known major solar system object"));
302 : }
303 :
304 0 : Vector<Double> planetparam;
305 0 : Vector<Double> earthparam;
306 0 : MEpoch::Convert toTDB(ep, MEpoch::TDB);
307 0 : earthparam=MeasTable::Planetary(MeasTable::EARTH, toTDB(ep).get("d").getValue());
308 0 : planetparam=MeasTable::Planetary(mtype, toTDB(ep).get("d").getValue());
309 : //GEOcentric param
310 0 : planetparam=planetparam-earthparam;
311 0 : Vector<Double> unitdirvec(3);
312 0 : Double dist=sqrt(planetparam(0)*planetparam(0)+planetparam(1)*planetparam(1)+planetparam(2)*planetparam(2));
313 0 : unitdirvec(0)=planetparam(0)/dist;
314 0 : unitdirvec(1)=planetparam(1)/dist;
315 0 : unitdirvec(2)=planetparam(2)/dist;
316 0 : MRadialVelocity::Convert obsvelconv(MRadialVelocity(MVRadialVelocity(0.0),
317 0 : MRadialVelocity::Ref(MRadialVelocity::TOPO, mFrame)), MRadialVelocity::Ref(refvel));
318 0 : Quantity planetradvel(planetparam(3)*unitdirvec(0)+planetparam(4)*unitdirvec(1)+planetparam(5)*unitdirvec(2), "AU/d");
319 0 : toSource=MDoppler(Quantity(-planetradvel.getValue("km/s")+obsvelconv().get("km/s").getValue("km/s") , "km/s"), MDoppler::RELATIVISTIC);
320 0 : toObs=MDoppler(Quantity(planetradvel.getValue("km/s")-obsvelconv().get("km/s").getValue("km/s") , "km/s"), MDoppler::RELATIVISTIC);
321 : }
322 0 : }
323 :
324 :
325 0 : void MSUtil::getSpwInFreqRangeAllFields(Vector<Int>& outspw, Vector<Int>& outstart,
326 : Vector<Int>& outnchan,
327 : const MeasurementSet& ms,
328 : const Double freqStart,
329 : const Double freqEnd,
330 : const Double freqStep,
331 : const MFrequency::Types freqframe){
332 0 : Vector<Int> fldId;
333 0 : ScalarColumn<Int> (ms, MS::columnName (MS::FIELD_ID)).getColumn (fldId);
334 0 : const Int option = Sort::HeapSort | Sort::NoDuplicates;
335 0 : const Sort::Order order = Sort::Ascending;
336 :
337 0 : Int nfields = GenSort<Int>::sort (fldId, order, option);
338 :
339 0 : fldId.resize(nfields, true);
340 0 : outspw.resize();
341 0 : outstart.resize();
342 0 : outnchan.resize();
343 0 : for (Int k=0; k < nfields; ++k){
344 0 : Vector<Int> locspw, locstart, locnchan;
345 0 : MSUtil::getSpwInFreqRange(locspw, locstart, locnchan, ms, freqStart, freqEnd, freqStep, freqframe, fldId[k]);
346 0 : for (Int j=0; j< locspw.shape()(0); ++j ){
347 0 : Bool hasthisspw=false;
348 0 : for (Int i=0; i< outspw.shape()(0); ++i){
349 0 : if(outspw[i]==locspw[j]){
350 0 : hasthisspw=true;
351 0 : if(locstart[j] < outstart[i]){
352 0 : Int endchan=outstart[i]+outnchan[i]-1;
353 0 : outstart[i]=locstart[j];
354 0 : if(endchan < (locstart[j]+ locnchan[j]-1))
355 0 : endchan=locstart[j]+ locnchan[j]-1;
356 0 : outnchan[i]=endchan+outstart[i]+1;
357 : }
358 : }
359 : }
360 0 : if(!hasthisspw){
361 0 : uInt nout=outspw.nelements();
362 0 : outspw.resize(nout+1, true);
363 0 : outnchan.resize(nout+1, true);
364 0 : outstart.resize(nout+1, true);
365 0 : outspw[nout]=locspw[j];
366 0 : outstart[nout]=locstart[j];
367 0 : outnchan[nout]=locnchan[j];
368 :
369 :
370 :
371 : }
372 : }
373 :
374 0 : }
375 :
376 :
377 0 : }
378 17 : void MSUtil::getChannelRangeFromFreqRange(Int& start,
379 : Int& nchan,
380 : const MeasurementSet& ms,
381 : const Int spw,
382 : const Double freqStart,
383 : const Double freqEnd,
384 : const Double freqStep,
385 : const MFrequency::Types freqframe)
386 : {
387 17 : start=-1;
388 17 : nchan=-1;
389 17 : MSFieldColumns fieldCol(ms.field());
390 17 : MSDataDescColumns ddCol(ms.dataDescription());
391 17 : Vector<Int> dataDescSel=MSDataDescIndex(ms.dataDescription()).matchSpwId(spw);
392 : //cerr << "dataDescSel " << dataDescSel << endl;
393 17 : if(dataDescSel.nelements()==0)
394 0 : return;
395 17 : Vector<Double> t;
396 17 : ScalarColumn<Double> (ms,MS::columnName(MS::TIME)).getColumn(t);
397 17 : Vector<Int> ddId;
398 17 : ScalarColumn<Int> (ms,MS::columnName(MS::DATA_DESC_ID)).getColumn(ddId);
399 17 : MSSpWindowColumns spwCol(ms.spectralWindow());
400 17 : Vector<Double> ddIdD(ddId.shape());
401 17 : convertArray(ddIdD, ddId);
402 17 : ddIdD+= 1.0; //no zero id
403 : //we have to do this as one can have multiple dd for the same time.
404 17 : t*=ddIdD;
405 : //t(fldId != fieldId)=-1.0;
406 17 : Vector<Double> elt;
407 17 : Vector<Int> elindx;
408 : //rejecting the large blocks of same time for all baselines
409 : //this speeds up by a lot GenSort::sort
410 17 : rejectConsecutive(t, elt, elindx);
411 :
412 17 : ROScalarMeasColumn<MEpoch> timeCol(ms, MS::columnName(MS::TIME));
413 17 : Vector<uInt> uniqIndx;
414 17 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, elt, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
415 :
416 17 : t.resize(0);
417 :
418 : //ScalarColumn<Int> ddId(ms,MS::columnName(MS::DATA_DESC_ID));
419 17 : ScalarColumn<Int> fldId(ms,MS::columnName(MS::FIELD_ID));
420 : //now need to do the conversion to data frame from requested frame
421 : //Get the epoch mesasures of the first row
422 17 : MEpoch ep;
423 17 : timeCol.get(0, ep);
424 17 : String observatory;
425 17 : MPosition obsPos;
426 : /////observatory position
427 17 : MSColumns msc(ms);
428 17 : if (ms.observation().nrow() > 0) {
429 17 : observatory = msc.observation().telescopeName()(msc.observationId()(0));
430 : }
431 34 : if (observatory.length() == 0 ||
432 17 : !MeasTable::Observatory(obsPos,observatory)) {
433 : // unknown observatory, use first antenna
434 0 : obsPos=msc.antenna().positionMeas()(0);
435 : }
436 : //////
437 17 : Int oldDD=dataDescSel(0);
438 17 : Int newDD=oldDD;
439 : //For now we will assume that the field is not moving very far from polynome 0
440 17 : MDirection dir =fieldCol.phaseDirMeas(0);
441 17 : MFrequency::Types obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(ddCol.spectralWindowId()(dataDescSel(0))));
442 17 : MeasFrame frame(ep, obsPos, dir);
443 : MFrequency::Convert toObs(freqframe,
444 17 : MFrequency::Ref(obsMFreqType, frame));
445 17 : Double freqEndMax=freqEnd+0.5*fabs(freqStep);
446 17 : Double freqStartMin=freqStart-0.5*fabs(freqStep);
447 17 : if(freqframe != obsMFreqType){
448 17 : freqEndMax=0.0;
449 17 : freqStartMin=C::dbl_max;
450 : }
451 8177 : for (uInt j=0; j< nTimes; ++j) {
452 8160 : if(anyEQ(dataDescSel, ddId(elindx[uniqIndx[j]]))){
453 8160 : timeCol.get(elindx[uniqIndx[j]], ep);
454 8160 : newDD=ddId(elindx[uniqIndx[j]]);
455 8160 : if(oldDD != newDD) {
456 :
457 0 : oldDD=newDD;
458 0 : if(spwCol.measFreqRef()(ddCol.spectralWindowId()(newDD)) != obsMFreqType) {
459 0 : obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(ddCol.spectralWindowId()(newDD)));
460 0 : toObs.setOut(MFrequency::Ref(obsMFreqType, frame));
461 : }
462 : }
463 8160 : if(obsMFreqType != freqframe) {
464 8160 : dir=fieldCol.phaseDirMeas(fldId(elindx[uniqIndx[j]]));
465 8160 : frame.resetEpoch(ep);
466 8160 : frame.resetDirection(dir);
467 8160 : Double freqTmp=toObs(Quantity(freqStart, "Hz")).get("Hz").getValue();
468 8160 : freqStartMin=(freqStartMin > freqTmp) ? freqTmp : freqStartMin;
469 8160 : freqTmp=toObs(Quantity(freqEnd, "Hz")).get("Hz").getValue();
470 8160 : freqEndMax=(freqEndMax < freqTmp) ? freqTmp : freqEndMax;
471 : }
472 : }
473 : }
474 :
475 :
476 17 : MSSpwIndex spwIn(ms.spectralWindow());
477 17 : Vector<Int> spws;
478 17 : Vector<Int> starts;
479 17 : Vector<Int> nchans;
480 17 : if(!spwIn.matchFrequencyRange(freqStartMin-0.5*freqStep, freqEndMax+0.5*freqStep, spws, starts, nchans)){
481 0 : return;
482 : }
483 34 : for (uInt k=0; k < spws.nelements(); ++k){
484 17 : if(spws[k]==spw){
485 17 : start=starts[k];
486 17 : nchan=nchans[k];
487 : }
488 :
489 : }
490 17 : }
491 76 : Bool MSUtil::getFreqRangeInSpw( Double& freqStart,
492 : Double& freqEnd, const Vector<Int>& spw, const Vector<Int>& start,
493 : const Vector<Int>& nchan,
494 : const MeasurementSet& ms,
495 : const MFrequency::Types freqframe,
496 : const Int fieldId, const Bool fromEdge){
497 76 : Vector<Int> fields(1, fieldId);
498 76 : return MSUtil::getFreqRangeInSpw( freqStart, freqEnd, spw, start,
499 152 : nchan,ms, freqframe,fields, fromEdge);
500 :
501 :
502 76 : }
503 1883 : Bool MSUtil::getFreqRangeInSpw( Double& freqStart,
504 : Double& freqEnd, const Vector<Int>& spw, const Vector<Int>& start,
505 : const Vector<Int>& nchan,
506 : const MeasurementSet& ms,
507 : const MFrequency::Types freqframe,
508 : const Bool fromEdge){
509 1883 : Vector<Int> fields(0);
510 1883 : return MSUtil::getFreqRangeInSpw( freqStart, freqEnd, spw, start,
511 3766 : nchan,ms, freqframe,fields, fromEdge, True);
512 :
513 :
514 1883 : }
515 4339 : Bool MSUtil::getFreqRangeInSpw( Double& freqStart,
516 : Double& freqEnd, const Vector<Int>& spw, const Vector<Int>& start,
517 : const Vector<Int>& nchan,
518 : const MeasurementSet& ms,
519 : const MFrequency::Types freqframe,
520 : const Vector<Int>& fieldIds, const Bool fromEdge, const Bool useFieldsInMS){
521 :
522 4339 : Bool retval=False;
523 4339 : freqStart=C::dbl_max;
524 4339 : freqEnd=0.0;
525 4339 : Vector<Double> t;
526 4339 : ScalarColumn<Double> (ms,MS::columnName(MS::TIME)).getColumn(t);
527 4339 : Vector<Int> ddId;
528 4339 : Vector<Int> fldId;
529 4339 : ScalarColumn<Int> (ms,MS::columnName(MS::DATA_DESC_ID)).getColumn(ddId);
530 4339 : ScalarColumn<Int> (ms,MS::columnName(MS::FIELD_ID)).getColumn(fldId);
531 4339 : MSFieldColumns fieldCol(ms.field());
532 4339 : MSDataDescColumns ddCol(ms.dataDescription());
533 4339 : MSSpWindowColumns spwCol(ms.spectralWindow());
534 4339 : ROScalarMeasColumn<MEpoch> timeCol(ms, MS::columnName(MS::TIME));
535 4339 : Vector<Double> ddIdD(ddId.shape());
536 4339 : convertArray(ddIdD, ddId);
537 4339 : ddIdD+= 1.0; //no zero id
538 : //we have to do this as one can have multiple dd for the same time.
539 4339 : t*=ddIdD;
540 : //t(fldId != fieldId)=-1.0;
541 4339 : Vector<Double> elt;
542 4339 : Vector<Int> elindx;
543 : //rejecting the large blocks of same time for all baselines
544 : //this speeds up by a lot GenSort::sort
545 4339 : rejectConsecutive(t, elt, elindx);
546 4339 : Vector<uInt> uniqIndx;
547 :
548 4339 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, elt, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
549 :
550 4339 : MDirection dir;
551 4339 : if(useFieldsInMS)
552 1883 : dir=fieldCol.phaseDirMeas(fldId[0]);
553 : else
554 2456 : dir=fieldCol.phaseDirMeas(fieldIds[0]);
555 4339 : MSDataDescIndex mddin(ms.dataDescription());
556 4339 : MFrequency::Types obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(0));
557 4339 : MEpoch ep;
558 4339 : timeCol.get(0, ep);
559 4339 : String observatory;
560 4339 : MPosition obsPos;
561 : /////observatory position
562 4339 : MSColumns msc(ms);
563 4339 : if (ms.observation().nrow() > 0) {
564 4339 : observatory = msc.observation().telescopeName()(msc.observationId()(0));
565 : }
566 8678 : if (observatory.length() == 0 ||
567 4339 : !MeasTable::Observatory(obsPos,observatory)) {
568 : // unknown observatory, use first antenna
569 0 : obsPos=msc.antenna().positionMeas()(0);
570 : }
571 : //////
572 4339 : MeasFrame frame(ep, obsPos, dir);
573 :
574 :
575 9246 : for (uInt ispw =0 ; ispw < spw.nelements() ; ++ispw){
576 4907 : if(nchan[ispw]>0 && start[ispw] >-1){
577 4907 : Double freqStartObs=C::dbl_max;
578 4907 : Double freqEndObs=0.0;
579 4907 : Vector<Double> chanfreq=spwCol.chanFreq()(spw[ispw]);
580 4907 : Vector<Double> chanwid=spwCol.chanWidth()(spw[ispw]);
581 4907 : Vector<Int> ddOfSpw=mddin.matchSpwId(spw[ispw]);
582 193394 : for (Int ichan=start[ispw]; ichan<start[ispw]+nchan[ispw]; ++ichan){
583 188487 : if(fromEdge){
584 116323 : if(freqStartObs > (chanfreq[ichan]-fabs(chanwid[ichan])/2.0)) freqStartObs=chanfreq[ichan]-fabs(chanwid[ichan])/2.0;
585 116323 : if(freqEndObs < (chanfreq[ichan]+fabs(chanwid[ichan])/2.0)) freqEndObs=chanfreq[ichan]+fabs(chanwid[ichan])/2.0;
586 : }
587 : else{
588 72164 : if(freqStartObs > (chanfreq[ichan])) freqStartObs=chanfreq[ichan];
589 72164 : if(freqEndObs < (chanfreq[ichan])) freqEndObs=chanfreq[ichan];
590 : }
591 : }
592 4907 : obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(spw[ispw]));
593 4907 : if((obsMFreqType==MFrequency::REST) || (obsMFreqType==freqframe && obsMFreqType != MFrequency::TOPO)){
594 527 : if(freqStart > freqStartObs) freqStart=freqStartObs;
595 527 : if(freqEnd < freqStartObs) freqEnd=freqStartObs;
596 527 : if(freqStart > freqEndObs) freqStart=freqEndObs;
597 527 : if(freqEnd < freqEndObs) freqEnd=freqEndObs;
598 527 : retval=True;
599 : }
600 : else{
601 : MFrequency::Convert toframe(obsMFreqType,
602 4380 : MFrequency::Ref(freqframe, frame));
603 1067963 : for (uInt j=0; j< nTimes; ++j){
604 1063583 : if((useFieldsInMS || anyEQ(fieldIds, fldId[elindx[uniqIndx[j]]])) && anyEQ(ddOfSpw, ddId[elindx[uniqIndx[j]]])){
605 929273 : timeCol.get(elindx[uniqIndx[j]], ep);
606 929273 : dir=fieldCol.phaseDirMeas(fldId[elindx[uniqIndx[j]]]);
607 929273 : frame.resetEpoch(ep);
608 929273 : frame.resetDirection(dir);
609 929273 : Double freqTmp=toframe(Quantity(freqStartObs, "Hz")).get("Hz").getValue();
610 929273 : if(freqStart > freqTmp) freqStart=freqTmp;
611 929273 : if(freqEnd < freqTmp) freqEnd=freqTmp;
612 929273 : freqTmp=toframe(Quantity(freqEndObs, "Hz")).get("Hz").getValue();
613 929273 : if(freqStart > freqTmp) freqStart=freqTmp;
614 929273 : if(freqEnd < freqTmp) freqEnd=freqTmp;
615 929273 : retval=True;
616 : }
617 : }
618 4380 : }
619 4907 : }
620 : }
621 4339 : return retval;
622 4339 : }
623 :
624 :
625 5 : Bool MSUtil::getFreqRangeAndRefFreqShift( Double& freqStart,
626 : Double& freqEnd, Quantity& sysvel, const MEpoch& refEp, const Vector<Int>& spw, const Vector<Int>& start,
627 : const Vector<Int>& nchan,
628 : const MeasurementSet& ms,
629 : const String& ephemPath, const MDirection& trackDir,
630 : const Bool fromEdge){
631 :
632 5 : casacore::MDirection::Types planetType=MDirection::castType(trackDir.getRef().getType());
633 5 : if( (! Table::isReadable(ephemPath)) && ( (planetType <= MDirection::N_Types) || (planetType >= MDirection::COMET)))
634 0 : throw(AipsError("getFreqRange in SOURCE frame has to have a valid ephemeris table or major solar system object defined"));
635 5 : Bool isephem=False;
636 5 : MeasComet mcomet;
637 5 : MeasTable::Types mtype=MeasTable::BARYEARTH;
638 5 : if(Table::isReadable(ephemPath)){
639 5 : mcomet=MeasComet(Path(ephemPath).absoluteName());
640 5 : isephem=True;
641 : }
642 : else{
643 :
644 0 : if(planetType >=MDirection::MERCURY && planetType < MDirection::COMET){
645 : //Damn these enums are not in the same order
646 0 : switch(planetType){
647 0 : case MDirection::MERCURY :
648 0 : mtype=MeasTable::MERCURY;
649 0 : break;
650 0 : case MDirection::VENUS :
651 0 : mtype=MeasTable::VENUS;
652 0 : break;
653 0 : case MDirection::MARS :
654 0 : mtype=MeasTable::MARS;
655 0 : break;
656 0 : case MDirection::JUPITER :
657 0 : mtype=MeasTable::JUPITER;
658 0 : break;
659 0 : case MDirection::SATURN :
660 0 : mtype=MeasTable::SATURN;
661 0 : break;
662 0 : case MDirection::URANUS :
663 0 : mtype=MeasTable::URANUS;
664 0 : break;
665 0 : case MDirection::NEPTUNE :
666 0 : mtype=MeasTable::NEPTUNE;
667 0 : break;
668 0 : case MDirection::PLUTO :
669 0 : mtype=MeasTable::PLUTO;
670 0 : break;
671 0 : case MDirection::MOON :
672 0 : mtype=MeasTable::MOON;
673 0 : break;
674 0 : case MDirection::SUN :
675 0 : mtype=MeasTable::SUN;
676 0 : break;
677 0 : default:
678 0 : throw(AipsError("Cannot translate to known major solar system object"));
679 : }
680 : }
681 :
682 : }
683 :
684 :
685 10 : Vector<Double> planetparam;
686 10 : Vector<Double> earthparam;
687 :
688 :
689 :
690 :
691 :
692 :
693 5 : Bool retval=False;
694 5 : freqStart=C::dbl_max;
695 5 : freqEnd=0.0;
696 10 : Vector<Double> t;
697 5 : ScalarColumn<Double> (ms,MS::columnName(MS::TIME)).getColumn(t);
698 10 : Vector<Int> ddId;
699 10 : Vector<Int> fldId;
700 5 : ScalarColumn<Int> (ms,MS::columnName(MS::DATA_DESC_ID)).getColumn(ddId);
701 5 : ScalarColumn<Int> (ms,MS::columnName(MS::FIELD_ID)).getColumn(fldId);
702 10 : MSFieldColumns fieldCol(ms.field());
703 10 : MSDataDescColumns ddCol(ms.dataDescription());
704 10 : MSSpWindowColumns spwCol(ms.spectralWindow());
705 10 : ROScalarMeasColumn<MEpoch> timeCol(ms, MS::columnName(MS::TIME));
706 10 : Vector<Double> ddIdD(ddId.shape());
707 5 : convertArray(ddIdD, ddId);
708 5 : ddIdD+= 1.0; //no zero id
709 : //we have to do this as one can have multiple dd for the same time.
710 5 : t*=ddIdD;
711 : //t(fldId != fieldId)=-1.0;
712 10 : Vector<Double> elt;
713 10 : Vector<Int> elindx;
714 : //rejecting the large blocks of same time for all baselines
715 : //this speeds up by a lot GenSort::sort
716 5 : rejectConsecutive(t, elt, elindx);
717 10 : Vector<uInt> uniqIndx;
718 :
719 5 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, elt, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
720 :
721 10 : MDirection dir;
722 5 : dir=fieldCol.phaseDirMeas(fldId[0]);
723 10 : MSDataDescIndex mddin(ms.dataDescription());
724 10 : MEpoch ep;
725 5 : timeCol.get(0, ep);
726 10 : String observatory;
727 10 : MPosition obsPos;
728 : /////observatory position
729 10 : MSColumns msc(ms);
730 5 : if (ms.observation().nrow() > 0) {
731 5 : observatory = msc.observation().telescopeName()(msc.observationId()(0));
732 : }
733 10 : if (observatory.length() == 0 ||
734 5 : !MeasTable::Observatory(obsPos,observatory)) {
735 : // unknown observatory, use first antenna
736 0 : obsPos=msc.antenna().positionMeas()(0);
737 : }
738 : //////
739 : //cerr << "obspos " << obsPos << endl;
740 10 : MeasFrame mframe(ep, obsPos, dir);
741 : ////Will use UT for now for ephem tables as it is not clear that they are being
742 : ///filled with TDB as intended in MeasComet.h
743 10 : MEpoch::Convert toUT(ep, MEpoch::UT);
744 10 : MEpoch::Convert toTDB(ep, MEpoch::TDB);
745 5 : MRadialVelocity::Types refvel=MRadialVelocity::GEO;
746 5 : if(isephem){
747 : //cerr << "dist " << mcomet.getTopo().getLength("km") << endl;
748 5 : if(mcomet.getTopo().getLength("km").getValue() > 1.0e-3){
749 2 : refvel=MRadialVelocity::TOPO;
750 : }
751 : }
752 10 : MRadialVelocity::Convert obsvelconv (MRadialVelocity(MVRadialVelocity(0.0),
753 10 : MRadialVelocity::Ref(MRadialVelocity::TOPO, mframe)),
754 15 : MRadialVelocity::Ref(refvel));
755 5 : MVRadialVelocity cometvel;
756 :
757 10 : for (uInt ispw =0 ; ispw < spw.nelements() ; ++ispw){
758 5 : if(nchan[ispw]>0 && start[ispw] >-1){
759 5 : Double freqStartObs=C::dbl_max;
760 5 : Double freqEndObs=0.0;
761 5 : Vector<Double> chanfreq=spwCol.chanFreq()(spw[ispw]);
762 5 : Vector<Double> chanwid=spwCol.chanWidth()(spw[ispw]);
763 5 : Vector<Int> ddOfSpw=mddin.matchSpwId(spw[ispw]);
764 21 : for (Int ichan=start[ispw]; ichan<start[ispw]+nchan[ispw]; ++ichan){
765 16 : if(fromEdge){
766 16 : if(freqStartObs > (chanfreq[ichan]-fabs(chanwid[ichan])/2.0)) freqStartObs=chanfreq[ichan]-fabs(chanwid[ichan])/2.0;
767 16 : if(freqEndObs < (chanfreq[ichan]+fabs(chanwid[ichan])/2.0)) freqEndObs=chanfreq[ichan]+fabs(chanwid[ichan])/2.0;
768 : }
769 : else{
770 0 : if(freqStartObs > (chanfreq[ichan])) freqStartObs=chanfreq[ichan];
771 0 : if(freqEndObs < (chanfreq[ichan])) freqEndObs=chanfreq[ichan];
772 : }
773 : }
774 :
775 5 : MFrequency::Types obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(spw[ispw]));
776 5 : if(obsMFreqType==MFrequency::REST)
777 0 : throw(AipsError("cannot do Source frame conversion from REST"));
778 : MFrequency::Convert toTOPO(obsMFreqType,
779 5 : MFrequency::Ref(MFrequency::TOPO, mframe));
780 5 : Double diffepoch=1e37;
781 5 : sysvel=Quantity(0.0,"m/s");
782 :
783 1782 : for (uInt j=0; j< nTimes; ++j){
784 1777 : if(anyEQ(ddOfSpw, ddId[elindx[uniqIndx[j]]])){
785 1777 : timeCol.get(elindx[uniqIndx[j]], ep);
786 1777 : dir=fieldCol.phaseDirMeas(fldId[elindx[uniqIndx[j]]], ep.get("s").getValue());
787 1777 : mframe.resetEpoch(ep);
788 1777 : mframe.resetDirection(dir);
789 1777 : if(obsMFreqType != MFrequency::TOPO){
790 0 : freqStartObs=toTOPO(Quantity(freqStartObs, "Hz")).get("Hz").getValue();
791 0 : freqEndObs=toTOPO(Quantity(freqEndObs, "Hz")).get("Hz").getValue();
792 : }
793 1777 : MDoppler mdop;
794 1777 : if(isephem){
795 1777 : mcomet.getRadVel(cometvel, toUT(ep).get("d").getValue());
796 1777 : mdop=MDoppler(Quantity(-cometvel.get("km/s").getValue("km/s")+obsvelconv().get("km/s").getValue("km/s") , "km/s"), MDoppler::RELATIVISTIC);
797 : // cerr << std::setprecision(10) << toUT(ep).get("d").getValue() << " fieldid " << fldId[elindx[uniqIndx[j]]] << " cometvel " << cometvel.get("km/s").getValue("km/s") << " obsvel " << obsvelconv().get("km/s").getValue("km/s") << endl;
798 : }
799 : else{
800 0 : earthparam=MeasTable::Planetary(MeasTable::EARTH, toTDB(ep).get("d").getValue());
801 0 : planetparam=MeasTable::Planetary(mtype, toTDB(ep).get("d").getValue());
802 : //GEOcentric param
803 0 : planetparam=planetparam-earthparam;
804 0 : Vector<Double> unitdirvec(3);
805 0 : Double dist=sqrt(planetparam(0)*planetparam(0)+planetparam(1)*planetparam(1)+planetparam(2)*planetparam(2));
806 0 : unitdirvec(0)=planetparam(0)/dist;
807 0 : unitdirvec(1)=planetparam(1)/dist;
808 0 : unitdirvec(2)=planetparam(2)/dist;
809 0 : Quantity planetradvel(planetparam(3)*unitdirvec(0)+planetparam(4)*unitdirvec(1)+planetparam(5)*unitdirvec(2), "AU/d");
810 0 : mdop=MDoppler(Quantity(-planetradvel.getValue("km/s")+obsvelconv().get("km/s").getValue("km/s") , "km/s"), MDoppler::RELATIVISTIC);
811 :
812 0 : }
813 1777 : Vector<Double> range(2); range(0)=freqStartObs; range(1)=freqEndObs;
814 1777 : range=mdop.shiftFrequency(range);
815 1777 : if(diffepoch > fabs(ep.get("s").getValue("s")-refEp.get("s").getValue("s"))){
816 5 : diffepoch= fabs(ep.get("s").getValue("s")-refEp.get("s").getValue("s"));
817 5 : sysvel=mdop.get("km/s");
818 : //cerr << std::setprecision(10) << "shifts " << range(0)-freqStartObs << " " << range(1)-freqEndObs << endl;
819 : }
820 :
821 :
822 1777 : if(freqStart > range[0]) freqStart=range[0];
823 1777 : if(freqEnd < range[0]) freqEnd=range[0];
824 :
825 1777 : if(freqStart > range[1]) freqStart=range[1];
826 1777 : if(freqEnd < range[1]) freqEnd=range[1];
827 1777 : retval=True;
828 1777 : }
829 : }
830 5 : }
831 :
832 : }
833 10 : return retval;
834 5 : }
835 4361 : void MSUtil::rejectConsecutive(const Vector<Double>& t, Vector<Double>& retval, Vector<Int>& indx){
836 4361 : uInt n=t.nelements();
837 4361 : if(n >0){
838 4361 : retval.resize(n);
839 4361 : indx.resize(n);
840 4361 : retval[0]=t[0];
841 4361 : indx[0]=0;
842 : }
843 : else
844 0 : return;
845 4361 : Int prev=0;
846 282287347 : for (uInt k=1; k < n; ++k){
847 282282986 : if(t[k] != retval(prev)){
848 1642956 : ++prev;
849 : //retval.resize(prev+1, true);
850 1642956 : retval[prev]=t[k];
851 : //indx.resize(prev+1, true);
852 1642956 : indx[prev]=k;
853 : }
854 : }
855 4361 : retval.resize(prev+1, true);
856 4361 : indx.resize(prev+1, true);
857 :
858 : }
859 :
860 0 : Vector<String> MSUtil::getSpectralFrames(Vector<MFrequency::Types>& types, const MeasurementSet& ms)
861 : {
862 0 : Vector<String> retval;
863 0 : Vector<Int> typesAsInt=MSSpWindowColumns(ms.spectralWindow()).measFreqRef().getColumn();
864 0 : if(ms.nrow()==Table(ms.getPartNames()).nrow()){
865 0 : types.resize(typesAsInt.nelements());
866 0 : for (uInt k=0; k < types.nelements(); ++k)
867 0 : types[k]=MFrequency::castType(typesAsInt[k]);
868 : }
869 : else{
870 0 : Vector<Int> ddId;
871 0 : ScalarColumn<Int> (ms,MS::columnName(MS::DATA_DESC_ID)).getColumn(ddId);
872 0 : Vector<uInt> uniqIndx;
873 0 : uInt nTimes=GenSort<Int>::sort (ddId, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
874 0 : ddId.resize(nTimes, true);
875 0 : Vector<Int> spwids(nTimes);
876 0 : Vector<Int> spwInDD=MSDataDescColumns(ms.dataDescription()).spectralWindowId().getColumn();
877 0 : for (uInt k=0; k < nTimes; ++k)
878 0 : spwids[k]=spwInDD[ddId[k]];
879 :
880 0 : nTimes=GenSort<Int>::sort (spwids, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
881 0 : spwids.resize(nTimes, true);
882 0 : types.resize(nTimes);
883 0 : for(uInt k=0; k <nTimes; ++k)
884 0 : types[k]=MFrequency::castType(typesAsInt[spwids[k]]);
885 0 : }
886 :
887 0 : retval.resize(types.nelements());
888 0 : for (uInt k=0; k < types.nelements(); ++k){
889 0 : retval[k]=MFrequency::showType(types[k]);
890 : }
891 :
892 :
893 0 : return retval;
894 :
895 0 : }
896 42 : void MSUtil::getIndexCombination(const MSColumns& mscol, Matrix<Int>& retval2){
897 42 : Vector<Vector<Int> >retval;
898 42 : Vector<Int> state = mscol.stateId().getColumn();
899 42 : Vector<Int> scan=mscol.scanNumber().getColumn();
900 42 : Vector<Double> t=mscol.time().getColumn();
901 42 : Vector<Int> fldid=mscol.fieldId().getColumn();
902 42 : Vector<Int> ddId=mscol.dataDescId().getColumn();
903 42 : Vector<Int> spwid=mscol.dataDescription().spectralWindowId().getColumn();
904 42 : Vector<uInt> uniqIndx;
905 : {
906 42 : Vector<Double> ddIdD(ddId.shape());
907 42 : convertArray(ddIdD, ddId);
908 42 : ddIdD+= 1.0; //no zero id
909 : //we have to do this as one can have multiple dd for the same time.
910 42 : t*=ddIdD;
911 42 : }
912 42 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, t, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
913 42 : Vector<Int> comb(4);
914 :
915 11733 : for (uInt k=0; k < nTimes; ++k){
916 11691 : comb(0)=fldid[uniqIndx[k]];
917 11691 : comb(1)=spwid[ddId[uniqIndx[k]]];
918 11691 : comb(2)=scan(uniqIndx[k]);
919 11691 : comb(3)=state(uniqIndx[k]);
920 11691 : Bool hasComb=False;
921 11691 : if(retval.nelements() >0){
922 30493 : for (uInt j=0; j < retval.nelements(); ++j){
923 18844 : if(allEQ(retval[j], comb))
924 11579 : hasComb=True;
925 : }
926 :
927 : }
928 11691 : if(!hasComb){
929 112 : retval.resize(retval.nelements()+1, True);
930 112 : retval[retval.nelements()-1]=comb;
931 : }
932 : }
933 42 : retval2.resize(retval.nelements(),4);
934 154 : for (uInt j=0; j < retval.nelements(); ++j)
935 112 : retval2.row(j)=retval[j];
936 :
937 42 : }
938 :
939 :
940 : } //# NAMESPACE CASA - END
|
