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 0 : 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 0 : spw.resize();
52 0 : start.resize();
53 0 : nchan.resize();
54 0 : Vector<Double> t;
55 0 : ScalarColumn<Double> (ms,MS::columnName(MS::TIME)).getColumn(t);
56 : //Vector<Int> ddId;
57 : //Vector<Int> fldId;
58 :
59 0 : MSFieldColumns fieldCol(ms.field());
60 0 : MSDataDescColumns ddCol(ms.dataDescription());
61 0 : MSSpWindowColumns spwCol(ms.spectralWindow());
62 0 : ROScalarMeasColumn<MEpoch> timeCol(ms, MS::columnName(MS::TIME));
63 0 : Vector<uInt> uniqIndx;
64 0 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, t, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
65 :
66 0 : 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 0 : ScalarColumn<Int> ddId(ms,MS::columnName(MS::DATA_DESC_ID));
70 0 : 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 0 : MEpoch ep;
74 0 : timeCol.get(0, ep);
75 0 : String observatory;
76 0 : MPosition obsPos;
77 : /////observatory position
78 0 : MSColumns msc(ms);
79 0 : if (ms.observation().nrow() > 0) {
80 0 : observatory = msc.observation().telescopeName()(msc.observationId()(0));
81 : }
82 0 : if (observatory.length() == 0 ||
83 0 : !MeasTable::Observatory(obsPos,observatory)) {
84 : // unknown observatory, use first antenna
85 0 : obsPos=msc.antenna().positionMeas()(0);
86 : }
87 : //////
88 0 : Int oldDD=ddId(0);
89 0 : Int newDD=oldDD;
90 : //For now we will assume that the field is not moving very far from polynome 0
91 0 : MDirection dir =fieldCol.phaseDirMeas(fieldId);
92 0 : MFrequency::Types obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(ddCol.spectralWindowId()(ddId(0))));
93 : //cout << "nTimes " << nTimes << endl;
94 : //cout << " obsframe " << obsMFreqType << " reqFrame " << freqframe << endl;
95 0 : MeasFrame frame(ep, obsPos, dir);
96 : MFrequency::Convert toObs(freqframe,
97 0 : MFrequency::Ref(obsMFreqType, frame));
98 0 : Double freqEndMax=freqEnd;
99 0 : Double freqStartMin=freqStart;
100 0 : if(freqframe != obsMFreqType){
101 0 : freqEndMax=0.0;
102 0 : freqStartMin=C::dbl_max;
103 : }
104 0 : for (uInt j=0; j< nTimes; ++j){
105 0 : if(fldId(uniqIndx[j]) ==fieldId){
106 0 : timeCol.get(uniqIndx[j], ep);
107 0 : newDD=ddId(uniqIndx[j]);
108 0 : 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 0 : if(obsMFreqType != freqframe){
116 0 : frame.resetEpoch(ep);
117 0 : Double freqTmp=toObs(Quantity(freqStart, "Hz")).get("Hz").getValue();
118 0 : freqStartMin=(freqStartMin > freqTmp) ? freqTmp : freqStartMin;
119 0 : freqTmp=toObs(Quantity(freqEnd, "Hz")).get("Hz").getValue();
120 0 : freqEndMax=(freqEndMax < freqTmp) ? freqTmp : freqEndMax;
121 : }
122 : }
123 : }
124 :
125 : //cout << "freqStartMin " << freqStartMin << " freqEndMax " << freqEndMax << endl;
126 0 : MSSpwIndex spwIn(ms.spectralWindow());
127 0 : spwIn.matchFrequencyRange(freqStartMin-0.5*freqStep, freqEndMax+0.5*freqStep, spw, start, nchan);
128 :
129 :
130 :
131 0 : }
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 0 : 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 0 : start=-1;
388 0 : nchan=-1;
389 0 : MSFieldColumns fieldCol(ms.field());
390 0 : MSDataDescColumns ddCol(ms.dataDescription());
391 0 : Vector<Int> dataDescSel=MSDataDescIndex(ms.dataDescription()).matchSpwId(spw);
392 : //cerr << "dataDescSel " << dataDescSel << endl;
393 0 : if(dataDescSel.nelements()==0)
394 0 : return;
395 0 : Vector<Double> t;
396 0 : ScalarColumn<Double> (ms,MS::columnName(MS::TIME)).getColumn(t);
397 0 : Vector<Int> ddId;
398 0 : ScalarColumn<Int> (ms,MS::columnName(MS::DATA_DESC_ID)).getColumn(ddId);
399 0 : MSSpWindowColumns spwCol(ms.spectralWindow());
400 0 : Vector<Double> ddIdD(ddId.shape());
401 0 : convertArray(ddIdD, ddId);
402 0 : ddIdD+= 1.0; //no zero id
403 : //we have to do this as one can have multiple dd for the same time.
404 0 : t*=ddIdD;
405 : //t(fldId != fieldId)=-1.0;
406 0 : Vector<Double> elt;
407 0 : Vector<Int> elindx;
408 : //rejecting the large blocks of same time for all baselines
409 : //this speeds up by a lot GenSort::sort
410 0 : rejectConsecutive(t, elt, elindx);
411 :
412 0 : ROScalarMeasColumn<MEpoch> timeCol(ms, MS::columnName(MS::TIME));
413 0 : Vector<uInt> uniqIndx;
414 0 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, elt, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
415 :
416 0 : t.resize(0);
417 :
418 : //ScalarColumn<Int> ddId(ms,MS::columnName(MS::DATA_DESC_ID));
419 0 : 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 0 : MEpoch ep;
423 0 : timeCol.get(0, ep);
424 0 : String observatory;
425 0 : MPosition obsPos;
426 : /////observatory position
427 0 : MSColumns msc(ms);
428 0 : if (ms.observation().nrow() > 0) {
429 0 : observatory = msc.observation().telescopeName()(msc.observationId()(0));
430 : }
431 0 : if (observatory.length() == 0 ||
432 0 : !MeasTable::Observatory(obsPos,observatory)) {
433 : // unknown observatory, use first antenna
434 0 : obsPos=msc.antenna().positionMeas()(0);
435 : }
436 : //////
437 0 : Int oldDD=dataDescSel(0);
438 0 : Int newDD=oldDD;
439 : //For now we will assume that the field is not moving very far from polynome 0
440 0 : MDirection dir =fieldCol.phaseDirMeas(0);
441 0 : MFrequency::Types obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(ddCol.spectralWindowId()(dataDescSel(0))));
442 0 : MeasFrame frame(ep, obsPos, dir);
443 : MFrequency::Convert toObs(freqframe,
444 0 : MFrequency::Ref(obsMFreqType, frame));
445 0 : Double freqEndMax=freqEnd+0.5*fabs(freqStep);
446 0 : Double freqStartMin=freqStart-0.5*fabs(freqStep);
447 0 : if(freqframe != obsMFreqType){
448 0 : freqEndMax=0.0;
449 0 : freqStartMin=C::dbl_max;
450 : }
451 0 : for (uInt j=0; j< nTimes; ++j) {
452 0 : if(anyEQ(dataDescSel, ddId(elindx[uniqIndx[j]]))){
453 0 : timeCol.get(elindx[uniqIndx[j]], ep);
454 0 : newDD=ddId(elindx[uniqIndx[j]]);
455 0 : 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 0 : if(obsMFreqType != freqframe) {
464 0 : dir=fieldCol.phaseDirMeas(fldId(elindx[uniqIndx[j]]));
465 0 : frame.resetEpoch(ep);
466 0 : frame.resetDirection(dir);
467 0 : Double freqTmp=toObs(Quantity(freqStart, "Hz")).get("Hz").getValue();
468 0 : freqStartMin=(freqStartMin > freqTmp) ? freqTmp : freqStartMin;
469 0 : freqTmp=toObs(Quantity(freqEnd, "Hz")).get("Hz").getValue();
470 0 : freqEndMax=(freqEndMax < freqTmp) ? freqTmp : freqEndMax;
471 : }
472 : }
473 : }
474 :
475 :
476 0 : MSSpwIndex spwIn(ms.spectralWindow());
477 0 : Vector<Int> spws;
478 0 : Vector<Int> starts;
479 0 : Vector<Int> nchans;
480 0 : if(!spwIn.matchFrequencyRange(freqStartMin-0.5*freqStep, freqEndMax+0.5*freqStep, spws, starts, nchans)){
481 0 : return;
482 : }
483 0 : for (uInt k=0; k < spws.nelements(); ++k){
484 0 : if(spws[k]==spw){
485 0 : start=starts[k];
486 0 : nchan=nchans[k];
487 : }
488 :
489 : }
490 0 : }
491 0 : 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 0 : Vector<Int> fields(1, fieldId);
498 0 : return MSUtil::getFreqRangeInSpw( freqStart, freqEnd, spw, start,
499 0 : nchan,ms, freqframe,fields, fromEdge);
500 :
501 :
502 0 : }
503 0 : 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 0 : Vector<Int> fields(0);
510 0 : return MSUtil::getFreqRangeInSpw( freqStart, freqEnd, spw, start,
511 0 : nchan,ms, freqframe,fields, fromEdge, True);
512 :
513 :
514 0 : }
515 0 : 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 0 : Bool retval=False;
523 0 : freqStart=C::dbl_max;
524 0 : freqEnd=0.0;
525 0 : Vector<Double> t;
526 0 : ScalarColumn<Double> (ms,MS::columnName(MS::TIME)).getColumn(t);
527 0 : Vector<Int> ddId;
528 0 : Vector<Int> fldId;
529 0 : ScalarColumn<Int> (ms,MS::columnName(MS::DATA_DESC_ID)).getColumn(ddId);
530 0 : ScalarColumn<Int> (ms,MS::columnName(MS::FIELD_ID)).getColumn(fldId);
531 0 : MSFieldColumns fieldCol(ms.field());
532 0 : MSDataDescColumns ddCol(ms.dataDescription());
533 0 : MSSpWindowColumns spwCol(ms.spectralWindow());
534 0 : ROScalarMeasColumn<MEpoch> timeCol(ms, MS::columnName(MS::TIME));
535 0 : Vector<Double> ddIdD(ddId.shape());
536 0 : convertArray(ddIdD, ddId);
537 0 : ddIdD+= 1.0; //no zero id
538 : //we have to do this as one can have multiple dd for the same time.
539 0 : t*=ddIdD;
540 : //t(fldId != fieldId)=-1.0;
541 0 : Vector<Double> elt;
542 0 : Vector<Int> elindx;
543 : //rejecting the large blocks of same time for all baselines
544 : //this speeds up by a lot GenSort::sort
545 0 : rejectConsecutive(t, elt, elindx);
546 0 : Vector<uInt> uniqIndx;
547 :
548 0 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, elt, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
549 :
550 0 : MDirection dir;
551 0 : if(useFieldsInMS)
552 0 : dir=fieldCol.phaseDirMeas(fldId[0]);
553 : else
554 0 : dir=fieldCol.phaseDirMeas(fieldIds[0]);
555 0 : MSDataDescIndex mddin(ms.dataDescription());
556 0 : MFrequency::Types obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(0));
557 0 : MEpoch ep;
558 0 : timeCol.get(0, ep);
559 0 : String observatory;
560 0 : MPosition obsPos;
561 : /////observatory position
562 0 : MSColumns msc(ms);
563 0 : if (ms.observation().nrow() > 0) {
564 0 : observatory = msc.observation().telescopeName()(msc.observationId()(0));
565 : }
566 0 : if (observatory.length() == 0 ||
567 0 : !MeasTable::Observatory(obsPos,observatory)) {
568 : // unknown observatory, use first antenna
569 0 : obsPos=msc.antenna().positionMeas()(0);
570 : }
571 : //////
572 0 : MeasFrame frame(ep, obsPos, dir);
573 :
574 :
575 0 : for (uInt ispw =0 ; ispw < spw.nelements() ; ++ispw){
576 0 : if(nchan[ispw]>0 && start[ispw] >-1){
577 0 : Double freqStartObs=C::dbl_max;
578 0 : Double freqEndObs=0.0;
579 0 : Vector<Double> chanfreq=spwCol.chanFreq()(spw[ispw]);
580 0 : Vector<Double> chanwid=spwCol.chanWidth()(spw[ispw]);
581 0 : Vector<Int> ddOfSpw=mddin.matchSpwId(spw[ispw]);
582 0 : for (Int ichan=start[ispw]; ichan<start[ispw]+nchan[ispw]; ++ichan){
583 0 : if(fromEdge){
584 0 : if(freqStartObs > (chanfreq[ichan]-fabs(chanwid[ichan])/2.0)) freqStartObs=chanfreq[ichan]-fabs(chanwid[ichan])/2.0;
585 0 : if(freqEndObs < (chanfreq[ichan]+fabs(chanwid[ichan])/2.0)) freqEndObs=chanfreq[ichan]+fabs(chanwid[ichan])/2.0;
586 : }
587 : else{
588 0 : if(freqStartObs > (chanfreq[ichan])) freqStartObs=chanfreq[ichan];
589 0 : if(freqEndObs < (chanfreq[ichan])) freqEndObs=chanfreq[ichan];
590 : }
591 : }
592 0 : obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(spw[ispw]));
593 0 : if((obsMFreqType==MFrequency::REST) || (obsMFreqType==freqframe && obsMFreqType != MFrequency::TOPO)){
594 0 : if(freqStart > freqStartObs) freqStart=freqStartObs;
595 0 : if(freqEnd < freqStartObs) freqEnd=freqStartObs;
596 0 : if(freqStart > freqEndObs) freqStart=freqEndObs;
597 0 : if(freqEnd < freqEndObs) freqEnd=freqEndObs;
598 0 : retval=True;
599 : }
600 : else{
601 : MFrequency::Convert toframe(obsMFreqType,
602 0 : MFrequency::Ref(freqframe, frame));
603 0 : for (uInt j=0; j< nTimes; ++j){
604 0 : if((useFieldsInMS || anyEQ(fieldIds, fldId[elindx[uniqIndx[j]]])) && anyEQ(ddOfSpw, ddId[elindx[uniqIndx[j]]])){
605 0 : timeCol.get(elindx[uniqIndx[j]], ep);
606 0 : dir=fieldCol.phaseDirMeas(fldId[elindx[uniqIndx[j]]]);
607 0 : frame.resetEpoch(ep);
608 0 : frame.resetDirection(dir);
609 0 : Double freqTmp=toframe(Quantity(freqStartObs, "Hz")).get("Hz").getValue();
610 0 : if(freqStart > freqTmp) freqStart=freqTmp;
611 0 : if(freqEnd < freqTmp) freqEnd=freqTmp;
612 0 : freqTmp=toframe(Quantity(freqEndObs, "Hz")).get("Hz").getValue();
613 0 : if(freqStart > freqTmp) freqStart=freqTmp;
614 0 : if(freqEnd < freqTmp) freqEnd=freqTmp;
615 0 : retval=True;
616 : }
617 : }
618 0 : }
619 0 : }
620 : }
621 0 : return retval;
622 0 : }
623 :
624 :
625 0 : 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 0 : casacore::MDirection::Types planetType=MDirection::castType(trackDir.getRef().getType());
633 0 : 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 0 : Bool isephem=False;
636 0 : MeasComet mcomet;
637 0 : MeasTable::Types mtype=MeasTable::BARYEARTH;
638 0 : if(Table::isReadable(ephemPath)){
639 0 : mcomet=MeasComet(Path(ephemPath).absoluteName());
640 0 : 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 0 : Vector<Double> planetparam;
686 0 : Vector<Double> earthparam;
687 :
688 :
689 :
690 :
691 :
692 :
693 0 : Bool retval=False;
694 0 : freqStart=C::dbl_max;
695 0 : freqEnd=0.0;
696 0 : Vector<Double> t;
697 0 : ScalarColumn<Double> (ms,MS::columnName(MS::TIME)).getColumn(t);
698 0 : Vector<Int> ddId;
699 0 : Vector<Int> fldId;
700 0 : ScalarColumn<Int> (ms,MS::columnName(MS::DATA_DESC_ID)).getColumn(ddId);
701 0 : ScalarColumn<Int> (ms,MS::columnName(MS::FIELD_ID)).getColumn(fldId);
702 0 : MSFieldColumns fieldCol(ms.field());
703 0 : MSDataDescColumns ddCol(ms.dataDescription());
704 0 : MSSpWindowColumns spwCol(ms.spectralWindow());
705 0 : ROScalarMeasColumn<MEpoch> timeCol(ms, MS::columnName(MS::TIME));
706 0 : Vector<Double> ddIdD(ddId.shape());
707 0 : convertArray(ddIdD, ddId);
708 0 : ddIdD+= 1.0; //no zero id
709 : //we have to do this as one can have multiple dd for the same time.
710 0 : t*=ddIdD;
711 : //t(fldId != fieldId)=-1.0;
712 0 : Vector<Double> elt;
713 0 : Vector<Int> elindx;
714 : //rejecting the large blocks of same time for all baselines
715 : //this speeds up by a lot GenSort::sort
716 0 : rejectConsecutive(t, elt, elindx);
717 0 : Vector<uInt> uniqIndx;
718 :
719 0 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, elt, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
720 :
721 0 : MDirection dir;
722 0 : dir=fieldCol.phaseDirMeas(fldId[0]);
723 0 : MSDataDescIndex mddin(ms.dataDescription());
724 0 : MEpoch ep;
725 0 : timeCol.get(0, ep);
726 0 : String observatory;
727 0 : MPosition obsPos;
728 : /////observatory position
729 0 : MSColumns msc(ms);
730 0 : if (ms.observation().nrow() > 0) {
731 0 : observatory = msc.observation().telescopeName()(msc.observationId()(0));
732 : }
733 0 : if (observatory.length() == 0 ||
734 0 : !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 0 : 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 0 : MEpoch::Convert toUT(ep, MEpoch::UT);
744 0 : MEpoch::Convert toTDB(ep, MEpoch::TDB);
745 0 : MRadialVelocity::Types refvel=MRadialVelocity::GEO;
746 0 : if(isephem){
747 : //cerr << "dist " << mcomet.getTopo().getLength("km") << endl;
748 0 : if(mcomet.getTopo().getLength("km").getValue() > 1.0e-3){
749 0 : refvel=MRadialVelocity::TOPO;
750 : }
751 : }
752 0 : MRadialVelocity::Convert obsvelconv (MRadialVelocity(MVRadialVelocity(0.0),
753 0 : MRadialVelocity::Ref(MRadialVelocity::TOPO, mframe)),
754 0 : MRadialVelocity::Ref(refvel));
755 0 : MVRadialVelocity cometvel;
756 :
757 0 : for (uInt ispw =0 ; ispw < spw.nelements() ; ++ispw){
758 0 : if(nchan[ispw]>0 && start[ispw] >-1){
759 0 : Double freqStartObs=C::dbl_max;
760 0 : Double freqEndObs=0.0;
761 0 : Vector<Double> chanfreq=spwCol.chanFreq()(spw[ispw]);
762 0 : Vector<Double> chanwid=spwCol.chanWidth()(spw[ispw]);
763 0 : Vector<Int> ddOfSpw=mddin.matchSpwId(spw[ispw]);
764 0 : for (Int ichan=start[ispw]; ichan<start[ispw]+nchan[ispw]; ++ichan){
765 0 : if(fromEdge){
766 0 : if(freqStartObs > (chanfreq[ichan]-fabs(chanwid[ichan])/2.0)) freqStartObs=chanfreq[ichan]-fabs(chanwid[ichan])/2.0;
767 0 : 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 0 : MFrequency::Types obsMFreqType= (MFrequency::Types) (spwCol.measFreqRef()(spw[ispw]));
776 0 : if(obsMFreqType==MFrequency::REST)
777 0 : throw(AipsError("cannot do Source frame conversion from REST"));
778 : MFrequency::Convert toTOPO(obsMFreqType,
779 0 : MFrequency::Ref(MFrequency::TOPO, mframe));
780 0 : Double diffepoch=1e37;
781 0 : sysvel=Quantity(0.0,"m/s");
782 :
783 0 : for (uInt j=0; j< nTimes; ++j){
784 0 : if(anyEQ(ddOfSpw, ddId[elindx[uniqIndx[j]]])){
785 0 : timeCol.get(elindx[uniqIndx[j]], ep);
786 0 : dir=fieldCol.phaseDirMeas(fldId[elindx[uniqIndx[j]]], ep.get("s").getValue());
787 0 : mframe.resetEpoch(ep);
788 0 : mframe.resetDirection(dir);
789 0 : 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 0 : MDoppler mdop;
794 0 : if(isephem){
795 0 : mcomet.getRadVel(cometvel, toUT(ep).get("d").getValue());
796 0 : 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 0 : Vector<Double> range(2); range(0)=freqStartObs; range(1)=freqEndObs;
814 0 : range=mdop.shiftFrequency(range);
815 0 : if(diffepoch > fabs(ep.get("s").getValue("s")-refEp.get("s").getValue("s"))){
816 0 : diffepoch= fabs(ep.get("s").getValue("s")-refEp.get("s").getValue("s"));
817 0 : sysvel=mdop.get("km/s");
818 : //cerr << std::setprecision(10) << "shifts " << range(0)-freqStartObs << " " << range(1)-freqEndObs << endl;
819 : }
820 :
821 :
822 0 : if(freqStart > range[0]) freqStart=range[0];
823 0 : if(freqEnd < range[0]) freqEnd=range[0];
824 :
825 0 : if(freqStart > range[1]) freqStart=range[1];
826 0 : if(freqEnd < range[1]) freqEnd=range[1];
827 0 : retval=True;
828 0 : }
829 : }
830 0 : }
831 :
832 : }
833 0 : return retval;
834 0 : }
835 0 : void MSUtil::rejectConsecutive(const Vector<Double>& t, Vector<Double>& retval, Vector<Int>& indx){
836 0 : uInt n=t.nelements();
837 0 : if(n >0){
838 0 : retval.resize(n);
839 0 : indx.resize(n);
840 0 : retval[0]=t[0];
841 0 : indx[0]=0;
842 : }
843 : else
844 0 : return;
845 0 : Int prev=0;
846 0 : for (uInt k=1; k < n; ++k){
847 0 : if(t[k] != retval(prev)){
848 0 : ++prev;
849 : //retval.resize(prev+1, true);
850 0 : retval[prev]=t[k];
851 : //indx.resize(prev+1, true);
852 0 : indx[prev]=k;
853 : }
854 : }
855 0 : retval.resize(prev+1, true);
856 0 : 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 0 : void MSUtil::getIndexCombination(const MSColumns& mscol, Matrix<Int>& retval2){
897 0 : Vector<Vector<Int> >retval;
898 0 : Vector<Int> state = mscol.stateId().getColumn();
899 0 : Vector<Int> scan=mscol.scanNumber().getColumn();
900 0 : Vector<Double> t=mscol.time().getColumn();
901 0 : Vector<Int> fldid=mscol.fieldId().getColumn();
902 0 : Vector<Int> ddId=mscol.dataDescId().getColumn();
903 0 : Vector<Int> spwid=mscol.dataDescription().spectralWindowId().getColumn();
904 0 : Vector<uInt> uniqIndx;
905 : {
906 0 : Vector<Double> ddIdD(ddId.shape());
907 0 : convertArray(ddIdD, ddId);
908 0 : ddIdD+= 1.0; //no zero id
909 : //we have to do this as one can have multiple dd for the same time.
910 0 : t*=ddIdD;
911 0 : }
912 0 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, t, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
913 0 : Vector<Int> comb(4);
914 :
915 0 : for (uInt k=0; k < nTimes; ++k){
916 0 : comb(0)=fldid[uniqIndx[k]];
917 0 : comb(1)=spwid[ddId[uniqIndx[k]]];
918 0 : comb(2)=scan(uniqIndx[k]);
919 0 : comb(3)=state(uniqIndx[k]);
920 0 : Bool hasComb=False;
921 0 : if(retval.nelements() >0){
922 0 : for (uInt j=0; j < retval.nelements(); ++j){
923 0 : if(allEQ(retval[j], comb))
924 0 : hasComb=True;
925 : }
926 :
927 : }
928 0 : if(!hasComb){
929 0 : retval.resize(retval.nelements()+1, True);
930 0 : retval[retval.nelements()-1]=comb;
931 : }
932 : }
933 0 : retval2.resize(retval.nelements(),4);
934 0 : for (uInt j=0; j < retval.nelements(); ++j)
935 0 : retval2.row(j)=retval[j];
936 :
937 0 : }
938 :
939 :
940 : } //# NAMESPACE CASA - END
|
