Line data Source code
1 : //# VisBufferUtil.cc: VisBuffer Utilities
2 : //# Copyright (C) 1996,1997,2001
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: aips2-request@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 :
29 : #include <casacore/casa/aips.h>
30 :
31 : #include <casacore/casa/Utilities/Assert.h>
32 : #include <casacore/casa/Exceptions/Error.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/MatrixMath.h>
35 : #include <casacore/casa/Arrays/Array.h>
36 : #include <casacore/casa/Arrays/Vector.h>
37 : #include <casacore/casa/Arrays/Matrix.h>
38 : #include <casacore/casa/Arrays/Cube.h>
39 : #include <casacore/casa/Utilities/Sort.h>
40 : #include <casacore/casa/OS/Timer.h>
41 : #include <casacore/casa/OS/Path.h>
42 : #include <casacore/measures/Measures/UVWMachine.h>
43 : #include <casacore/measures/Measures/MeasTable.h>
44 : #include <casacore/ms/MSSel/MSSelectionTools.h>
45 : #include <casacore/ms/MeasurementSets/MSPointingColumns.h>
46 : #include <msvis/MSVis/VisBufferUtil.h>
47 : #include <msvis/MSVis/StokesVector.h>
48 : #include <msvis/MSVis/ViImplementation2.h>
49 : #include <msvis/MSVis/VisibilityIterator.h>
50 : #include <msvis/MSVis/VisibilityIterator2.h>
51 : #include <msvis/MSVis/VisBuffer.h>
52 : #include <casacore/ms/MeasurementSets/MSColumns.h>
53 : #include <iostream>
54 : #include <fstream>
55 : #include <iomanip>
56 : #ifdef _OPENMP
57 : #include <omp.h>
58 : #endif
59 : using namespace std;
60 :
61 : using namespace casacore;
62 : namespace casa { //# NAMESPACE CASA - BEGIN
63 :
64 : // <summary>
65 : // </summary>
66 :
67 : // <reviewed reviewer="" date="" tests="tMEGI" demos="">
68 :
69 : // <prerequisite>
70 : // </prerequisite>
71 : //
72 : // <etymology>
73 : // </etymology>
74 : //
75 : // <synopsis>
76 : // </synopsis>
77 : //
78 : // <example>
79 : // <srcblock>
80 : // </srcblock>
81 : // </example>
82 : //
83 : // <motivation>
84 : // </motivation>
85 : //
86 : // <todo asof="">
87 : // </todo>
88 :
89 :
90 0 : VisBufferUtil::VisBufferUtil(): oldMSId_p(-1), oldPCMSId_p(-1), timeAntIndex_p(0), cachedPointingDir_p(0), cachedPhaseCenter_p(0){};
91 :
92 :
93 : // Construct from a VisBuffer (sets frame info)
94 0 : VisBufferUtil::VisBufferUtil(const VisBuffer& vb): oldMSId_p(-1), oldPCMSId_p(-1), timeAntIndex_p(0), cachedPointingDir_p(0), cachedPhaseCenter_p(0) {
95 :
96 : // The nominal epoch
97 0 : MEpoch ep=vb.msColumns().timeMeas()(0);
98 :
99 : // The nominal position
100 0 : String observatory;
101 0 : MPosition pos;
102 0 : if (vb.msColumns().observation().nrow() > 0) {
103 0 : observatory = vb.msColumns().observation().telescopeName()
104 0 : (vb.msColumns().observationId()(0));
105 : }
106 0 : if (observatory.length() == 0 ||
107 0 : !MeasTable::Observatory(pos,observatory)) {
108 : // unknown observatory, use first antenna
109 0 : pos=vb.msColumns().antenna().positionMeas()(0);
110 : }
111 :
112 : // The nominal direction
113 0 : MDirection dir=vb.phaseCenter();
114 :
115 : // The nominal MeasFrame
116 0 : mframe_=MeasFrame(ep, pos, dir);
117 :
118 0 : }
119 :
120 : // Construct from a VisBuffer (sets frame info)
121 0 : VisBufferUtil::VisBufferUtil(const vi::VisBuffer2& vb): oldMSId_p(-1), oldPCMSId_p(-1), timeAntIndex_p(0), cachedPointingDir_p(0), cachedPhaseCenter_p(0) {
122 0 : if(!vb.isAttached())
123 0 : ThrowCc("Programmer Error: used a detached Visbuffer when it should not have been so");
124 0 : MSColumns msc(vb.ms());
125 : // The nominal epoch
126 0 : MEpoch ep=msc.timeMeas()(0);
127 :
128 : // The nominal position
129 0 : String observatory;
130 0 : MPosition pos;
131 0 : if (msc.observation().nrow() > 0) {
132 0 : observatory = msc.observation().telescopeName()
133 0 : (msc.observationId()(0));
134 : }
135 0 : if (observatory.length() == 0 ||
136 0 : !MeasTable::Observatory(pos,observatory)) {
137 : // unknown observatory, use first antenna
138 0 : pos=msc.antenna().positionMeas()(0);
139 : }
140 :
141 : // The nominal direction
142 0 : MDirection dir=vb.phaseCenter();
143 :
144 : // The nominal MeasFrame
145 0 : mframe_=MeasFrame(ep, pos, dir);
146 :
147 0 : }
148 0 : VisBufferUtil::VisBufferUtil(const vi::VisibilityIterator2& iter): oldMSId_p(-1) {
149 :
150 0 : MSColumns msc(iter.ms());
151 : // The nominal epoch
152 0 : MEpoch ep=msc.timeMeas()(0);
153 :
154 : // The nominal position
155 0 : String observatory;
156 0 : MPosition pos;
157 0 : if (msc.observation().nrow() > 0) {
158 0 : observatory = msc.observation().telescopeName()
159 0 : (msc.observationId()(0));
160 : }
161 0 : if (observatory.length() == 0 ||
162 0 : !MeasTable::Observatory(pos,observatory)) {
163 : // unknown observatory, use first antenna
164 0 : pos=msc.antenna().positionMeas()(0);
165 : }
166 :
167 : // The nominal direction
168 : //MDirection dir=iter.phaseCenter();
169 0 : MDirection dir=msc.field().phaseDirMeasCol()(0)(IPosition(1,0));
170 : // The nominal MeasFrame
171 0 : mframe_=MeasFrame(ep, pos, dir);
172 :
173 0 : }
174 0 : VisBufferUtil::VisBufferUtil(const MeasFrame& mframe): oldMSId_p(-1) {
175 0 : mframe_=mframe;
176 :
177 0 : }
178 0 : Bool VisBufferUtil::rotateUVW(const vi::VisBuffer2&vb, const MDirection& desiredDir,
179 : Matrix<Double>& uvw, Vector<Double>& dphase){
180 :
181 0 : Bool retval=true;
182 0 : mframe_.resetEpoch(vb.time()(0));
183 0 : UVWMachine uvwMachine(desiredDir, vb.phaseCenter(), mframe_,
184 0 : false, false);
185 0 : retval = !uvwMachine.isNOP();
186 0 : dphase.resize(vb.nRows());
187 0 : dphase.set(0.0);
188 0 : if(uvw.nelements() ==0)
189 0 : uvw=vb.uvw();
190 0 : for (rownr_t row=0; row< vb.nRows(); ++row){
191 0 : Vector<Double> eluvw(uvw.column(row));
192 0 : uvwMachine.convertUVW(dphase(row), eluvw);
193 : }
194 :
195 0 : return retval;
196 : }
197 :
198 : // Set the visibility buffer for a PSF
199 0 : void VisBufferUtil::makePSFVisBuffer(VisBuffer& vb) {
200 0 : CStokesVector coh(Complex(1.0), Complex(0.0), Complex(0.0), Complex(1.0));
201 0 : vb.correctedVisibility()=coh;
202 0 : }
203 :
204 :
205 0 : Bool VisBufferUtil::interpolateFrequency(Cube<Complex>& data,
206 : Cube<Bool>& flag,
207 : const VisBuffer& vb,
208 : const Vector<Float>& outFreqGrid,
209 : const MS::PredefinedColumns whichCol,
210 : const MFrequency::Types freqFrame,
211 : const InterpolateArray1D<Float,Complex>::InterpolationMethod interpMethod){
212 :
213 0 : Cube<Complex> origdata;
214 : // Convert the visibility frequency to the frame requested
215 0 : Vector<Double> visFreqD;
216 0 : convertFrequency(visFreqD, vb, freqFrame);
217 : //convert it to Float
218 0 : Vector<Float> visFreq(visFreqD.nelements());
219 0 : convertArray(visFreq, visFreqD);
220 :
221 : //Assign which column is to be regridded to origdata
222 0 : if(whichCol==MS::MODEL_DATA){
223 0 : origdata.reference(vb.modelVisCube());
224 : }
225 0 : else if(whichCol==MS::CORRECTED_DATA){
226 0 : origdata.reference(vb.correctedVisCube());
227 : }
228 0 : else if(whichCol==MS::DATA){
229 0 : origdata.reference(vb.visCube());
230 : }
231 : else{
232 0 : throw(AipsError("Don't know which column is being regridded"));
233 : }
234 0 : Cube<Complex> flipdata;
235 0 : Cube<Bool> flipflag;
236 : //The interpolator interpolates on the 3rd axis only...so need to flip the axes (y,z)
237 0 : swapyz(flipflag,vb.flagCube());
238 0 : swapyz(flipdata,origdata);
239 :
240 : //interpolate the data and the flag to the output frequency grid
241 : InterpolateArray1D<Float,Complex>::
242 0 : interpolate(data,flag, outFreqGrid,visFreq,flipdata,flipflag,interpMethod);
243 0 : flipdata.resize();
244 : //reflip the data and flag to be in the same order as in Visbuffer output
245 0 : swapyz(flipdata,data);
246 0 : data.resize();
247 0 : data.reference(flipdata);
248 0 : flipflag.resize();
249 0 : swapyz(flipflag,flag);
250 0 : flag.resize();
251 0 : flag.reference(flipflag);
252 :
253 0 : return true;
254 :
255 : }
256 0 : void VisBufferUtil::getFreqRange(Double& freqMin, Double& freqMax, vi::VisibilityIterator2& vi, MFrequency::Types freqFrame){
257 0 : vi.originChunks();
258 0 : vi.origin();
259 :
260 0 : Double freqEnd=0.0;
261 0 : Double freqStart=C::dbl_max;
262 0 : vi::VisBuffer2* vb=vi.getVisBuffer();
263 0 : for (vi.originChunks(); vi.moreChunks();vi.nextChunk())
264 : {
265 0 : for (vi.origin(); vi.more();vi.next())
266 : {
267 : Double localmax, localmin;
268 0 : IPosition localmaxpos(1,0);
269 0 : IPosition localminpos(1,0);
270 0 : Vector<Double> freqs=vb->getFrequencies(0, freqFrame);
271 0 : if(freqs.nelements() ==0){
272 0 : throw(AipsError("Frequency selection error" ));
273 : }
274 0 : minMax(localmin,localmax,localminpos, localmaxpos, freqs);
275 : //localmax=max(freqs);
276 : //localmin=min(freqs);
277 : //freqEnd=max(freqEnd, localmax);
278 : //freqStart=min(freqStart, localmin);
279 :
280 0 : Int nfreq = freqs.nelements();
281 0 : Vector<Int> curspws = vb->spectralWindows();
282 : // as the vb row 0 is used for getFrequencies, the same row 0 is used here
283 0 : Vector<Double> chanWidths = vi.subtableColumns().spectralWindow().chanWidth()(curspws[0]);
284 : // freqs are in channel center freq so add the half the width to the values to return the edge frequencies
285 0 : if (nfreq==1) {
286 0 : freqEnd=max(freqEnd, localmax+fabs(chanWidths[0]/2.0));
287 0 : freqStart=min(freqStart, localmin-fabs(chanWidths[0]/2.0));
288 : }
289 : else {
290 0 : freqEnd=max(freqEnd, localmax+fabs(chanWidths[localmaxpos[0]]/2.0));
291 0 : freqStart=min(freqStart, localmin-fabs(chanWidths[localminpos[0]]/2.0));
292 : }
293 :
294 : }
295 : }
296 0 : freqMin=freqStart;
297 0 : freqMax=freqEnd;
298 0 : }
299 :
300 0 : Bool VisBufferUtil::getFreqRangeFromRange(casacore::Double& outfreqMin, casacore::Double& outfreqMax, const casacore::MFrequency::Types inFreqFrame, const casacore::Double infreqMin, const casacore::Double infreqMax, vi::VisibilityIterator2& vi, casacore::MFrequency::Types outFreqFrame){
301 :
302 :
303 0 : if(inFreqFrame==outFreqFrame){
304 0 : outfreqMin=infreqMin;
305 0 : outfreqMax=infreqMax;
306 0 : return True;
307 : }
308 :
309 :
310 0 : vi.originChunks();
311 0 : vi.origin();
312 : try{
313 0 : outfreqMin=C::dbl_max;
314 0 : outfreqMax=0;
315 0 : vi::VisBuffer2* vb=vi.getVisBuffer();
316 0 : MSColumns msc(vi.ms());
317 : // The nominal epoch
318 0 : MEpoch ep=msc.timeMeas()(0);
319 :
320 : // The nominal position
321 0 : String observatory;
322 0 : MPosition pos;
323 0 : if (msc.observation().nrow() > 0) {
324 0 : observatory = msc.observation().telescopeName()
325 0 : (msc.observationId()(0));
326 : }
327 0 : if (observatory.length() == 0 ||
328 0 : !MeasTable::Observatory(pos,observatory)) {
329 : // unknown observatory, use first antenna
330 0 : pos=msc.antenna().positionMeas()(0);
331 : }
332 :
333 : // The nominal direction
334 0 : MDirection dir=vb->phaseCenter();
335 0 : MeasFrame mFrame(ep, pos, dir);
336 : // The conversion engine:
337 : MFrequency::Convert toNewFrame(inFreqFrame,
338 0 : MFrequency::Ref(outFreqFrame, mFrame));
339 :
340 0 : for (vi.originChunks(); vi.moreChunks();vi.nextChunk())
341 : {
342 0 : for (vi.origin(); vi.more();vi.next()){
343 : //assuming time is fixed in visbuffer
344 0 : mFrame.resetEpoch(vb->time()(0)/86400.0);
345 :
346 : // Reset the direction (ASSUMES phaseCenter is constant in the VisBuffer)
347 0 : mFrame.resetDirection(vb->phaseCenter());
348 0 : Double temp=toNewFrame(infreqMin).getValue().getValue();
349 0 : if(temp < outfreqMin)
350 0 : outfreqMin = temp;
351 :
352 0 : temp=toNewFrame(infreqMax).getValue().getValue();
353 0 : if(temp > outfreqMax)
354 0 : outfreqMax = temp;
355 : }
356 : }
357 : }
358 0 : catch(...){
359 : //Could not do a conversion
360 0 : return False;
361 :
362 : }
363 : //cerr << "min " << outfreqMin << " max " << outfreqMax << endl;
364 0 : return True;
365 :
366 : }
367 :
368 0 : void VisBufferUtil::convertFrequency(Vector<Double>& outFreq,
369 : const VisBuffer& vb,
370 : const MFrequency::Types freqFrame){
371 0 : Int spw=vb.spectralWindow();
372 0 : MFrequency::Types obsMFreqType=(MFrequency::Types)(vb.msColumns().spectralWindow().measFreqRef()(spw));
373 :
374 : // The input frequencies (a reference)
375 0 : Vector<Double> inFreq(vb.frequency());
376 :
377 : // The output frequencies
378 0 : outFreq.resize(inFreq.nelements());
379 :
380 0 : MFrequency::Types newMFreqType=freqFrame;
381 0 : if (freqFrame==MFrequency::N_Types)
382 : // Opt out of conversion
383 0 : newMFreqType=obsMFreqType;
384 :
385 :
386 : // Only convert if the requested frame differs from observed frame
387 0 : if(obsMFreqType != newMFreqType){
388 :
389 : // Setting epoch to the first in this iteration
390 : // MEpoch ep=vb.msColumns().timeMeas()(0);
391 : // MEpoch ep(MVEpoch(vb.time()(0)/86400.0),MEpoch::UTC);
392 : // cout << "Time = " << ep.getValue() << endl;
393 :
394 : // Reset the timestamp (ASSUMES TIME is constant in the VisBuffer)
395 0 : mframe_.resetEpoch(vb.time()(0)/86400.0);
396 :
397 : // Reset the direction (ASSUMES phaseCenter is constant in the VisBuffer)
398 0 : mframe_.resetDirection(vb.msColumns().field().phaseDirMeasCol()(vb.fieldId())(IPosition(1,0)));
399 :
400 : // cout << "Frame = " << mframe_ << endl;
401 :
402 : // The conversion engine:
403 : MFrequency::Convert toNewFrame(obsMFreqType,
404 0 : MFrequency::Ref(newMFreqType, mframe_));
405 :
406 : // Do the conversion
407 0 : for (uInt k=0; k< inFreq.nelements(); ++k)
408 0 : outFreq(k)=toNewFrame(inFreq(k)).getValue().getValue();
409 :
410 : }
411 : else{
412 : // The requested frame is the same as the observed frame
413 0 : outFreq=inFreq;
414 : }
415 :
416 0 : }
417 :
418 0 : void VisBufferUtil::convertFrequency(Vector<Double>& outFreq,
419 : const vi::VisBuffer2& vb,
420 : const MFrequency::Types freqFrame){
421 0 : Int spw=vb.spectralWindows()(0);
422 0 : MFrequency::Types obsMFreqType=(MFrequency::Types)(MSColumns(vb.ms()).spectralWindow().measFreqRef()(spw));
423 :
424 :
425 :
426 : // The input frequencies
427 0 : Vector<Int> chanNums=vb.getChannelNumbers(0);
428 :
429 0 : Vector<Double> inFreq(chanNums.nelements());
430 0 : Vector<Double> spwfreqs=MSColumns(vb.ms()).spectralWindow().chanFreq().get(spw);
431 0 : for (uInt k=0; k < chanNums.nelements(); ++k){
432 :
433 0 : inFreq[k]=spwfreqs[chanNums[k]];
434 : }
435 :
436 : // The output frequencies
437 0 : outFreq.resize(inFreq.nelements());
438 :
439 0 : MFrequency::Types newMFreqType=freqFrame;
440 0 : if (freqFrame==MFrequency::N_Types)
441 : // Opt out of conversion
442 0 : newMFreqType=obsMFreqType;
443 :
444 :
445 : // Only convert if the requested frame differs from observed frame
446 0 : if(obsMFreqType != newMFreqType){
447 :
448 : // Setting epoch to the first in this iteration
449 : // MEpoch ep=vb.msColumns().timeMeas()(0);
450 : // MEpoch ep(MVEpoch(vb.time()(0)/86400.0),MEpoch::UTC);
451 : // cout << "Time = " << ep.getValue() << endl;
452 :
453 : // Reset the timestamp (ASSUMES TIME is constant in the VisBuffer)
454 0 : mframe_.resetEpoch(vb.time()(0)/86400.0);
455 :
456 : // Reset the direction (ASSUMES phaseCenter is constant in the VisBuffer)
457 0 : mframe_.resetDirection(vb.phaseCenter());
458 :
459 : // cout << "Frame = " << mframe_ << endl;
460 :
461 : // The conversion engine:
462 : MFrequency::Convert toNewFrame(obsMFreqType,
463 0 : MFrequency::Ref(newMFreqType, mframe_));
464 :
465 : // Do the conversion
466 0 : for (uInt k=0; k< inFreq.nelements(); ++k)
467 0 : outFreq(k)=toNewFrame(inFreq(k)).getValue().getValue();
468 :
469 : }
470 : else{
471 : // The requested frame is the same as the observed frame
472 0 : outFreq=inFreq;
473 : }
474 :
475 : //cerr << std::setprecision(9) << " infreq " << inFreq[152] << " " << outFreq[152] << " vb freq " << vb.getFrequencies(0, freqFrame)[152] << endl;
476 :
477 0 : }
478 :
479 0 : void VisBufferUtil::toVelocity(Vector<Double>& outVel,
480 : const VisBuffer& vb,
481 : const MFrequency::Types freqFrame,
482 : const MVFrequency restFreq,
483 : const MDoppler::Types veldef){
484 :
485 : // The input frequencies (a reference)
486 0 : Vector<Double> inFreq(vb.frequency());
487 :
488 : // The output velocities
489 0 : outVel.resize(inFreq.nelements());
490 :
491 : // Reset the timestamp (ASSUMES TIME is constant in the VisBuffer)
492 0 : mframe_.resetEpoch(vb.time()(0)/86400.0);
493 :
494 : // Reset the direction (ASSUMES phaseCenter is constant in the VisBuffer)
495 : //mframe_.resetDirection(vb.phaseCenter());
496 0 : mframe_.resetDirection(vb.msColumns().field().phaseDirMeasCol()(vb.fieldId())(IPosition(1,0)));
497 :
498 : // The frequency conversion engine:
499 0 : Int spw=vb.spectralWindow();
500 0 : MFrequency::Types obsMFreqType=(MFrequency::Types)(vb.msColumns().spectralWindow().measFreqRef()(spw));
501 :
502 0 : MFrequency::Types newMFreqType=freqFrame;
503 0 : if (freqFrame==MFrequency::N_Types)
504 : // Don't convert frame
505 0 : newMFreqType=obsMFreqType;
506 :
507 : MFrequency::Convert toNewFrame(obsMFreqType,
508 0 : MFrequency::Ref(newMFreqType, mframe_));
509 :
510 : // The velocity conversion engine:
511 0 : MDoppler::Ref dum1(MDoppler::RELATIVISTIC);
512 0 : MDoppler::Ref dum2(veldef);
513 0 : MDoppler::Convert dopConv(dum1, dum2);
514 :
515 : // Cope with unspecified rest freq
516 0 : MVFrequency rf=restFreq;
517 0 : if (restFreq.getValue()<=0.0)
518 0 : rf=toNewFrame(inFreq(vb.nChannel()/2)).getValue();
519 :
520 : // Do the conversions
521 0 : for (uInt k=0; k< inFreq.nelements(); ++k){
522 : //cerr <<"old freq " << toNewFrame(inFreq(k)).getValue().get().getValue() << endl;
523 0 : MDoppler eh = toNewFrame(inFreq(k)).toDoppler(rf);
524 0 : MDoppler eh2 = dopConv(eh);
525 0 : outVel(k)=eh2.getValue().get().getValue();
526 : }
527 :
528 0 : }
529 :
530 0 : void VisBufferUtil::toVelocity(Vector<Double>& outVel,
531 : const vi::VisBuffer2& vb,
532 : const MFrequency::Types freqFrame,
533 : const MVFrequency restFreq,
534 : const MDoppler::Types veldef, const Int row){
535 :
536 0 : Vector<Double> inFreq;
537 0 : inFreq=vb.getFrequencies(row, freqFrame);
538 : // cerr << "Freqs " << inFreq << endl;
539 : // The output velocities
540 0 : outVel.resize(inFreq.nelements());
541 : // The velocity conversion engine:
542 0 : MDoppler::Ref dum1(MDoppler::RELATIVISTIC);
543 0 : MDoppler::Ref dum2(veldef);
544 0 : MDoppler::Convert dopConv(dum1, dum2);
545 :
546 : // Cope with unspecified rest freq
547 0 : MVFrequency rf=restFreq;
548 0 : if (restFreq.getValue()<=0.0)
549 0 : rf=inFreq(inFreq.nelements()/2);
550 :
551 : // Do the conversions
552 0 : for (uInt k=0; k< inFreq.nelements(); ++k){
553 0 : MDoppler eh = MFrequency(Quantity(inFreq(k), "Hz"), freqFrame).toDoppler(rf);
554 0 : MDoppler eh2 = dopConv(eh);
555 0 : outVel(k)=eh2.getValue().get().getValue();
556 : }
557 :
558 :
559 0 : }
560 0 : void VisBufferUtil::toVelocity(Vector<Double>& outVel,
561 : const vi::VisBuffer2& vb,
562 : const vi::VisibilityIterator2& iter,
563 : const MFrequency::Types freqFrame,
564 : const MVFrequency restFreq,
565 : const MDoppler::Types veldef, const Int row){
566 :
567 : // The input frequencies (a reference)
568 0 : Vector<Double> inFreq(vb.getFrequencies(row));
569 0 : MSColumns msc(iter.ms());
570 :
571 0 : MEpoch ep(Quantity(vb.time()(row)/86400.0, "d"), msc.timeMeas()(0).getRef());
572 0 : MDirection dir(msc.field().phaseDirMeasCol()(vb.fieldId()(row))(IPosition(1,0)));
573 0 : Int spw=vb.spectralWindows()(row);
574 0 : MFrequency::Types obsMFreqType=(MFrequency::Types)(msc.spectralWindow().measFreqRef()(spw));
575 0 : toVelocity(outVel, freqFrame, inFreq, obsMFreqType, ep, dir, restFreq, veldef);
576 0 : }
577 0 : void VisBufferUtil::toVelocity(Vector<Double>& outVel,
578 : const MFrequency::Types outFreqFrame,
579 : const Vector<Double>& inFreq,
580 : const MFrequency::Types inFreqFrame,
581 : const MEpoch& ep,
582 : const MDirection& dir,
583 : const MVFrequency restFreq,
584 : const MDoppler::Types veldef){
585 :
586 :
587 :
588 : // The output velocities
589 0 : outVel.resize(inFreq.nelements());
590 :
591 : // Reset the timestamp
592 0 : mframe_.resetEpoch(ep);
593 :
594 : // Reset the direction
595 0 : mframe_.resetDirection(dir);
596 :
597 : // The frequency conversion engine:
598 :
599 0 : MFrequency::Types newMFreqType=outFreqFrame;
600 0 : if (outFreqFrame==MFrequency::N_Types)
601 : // Don't convert frame
602 0 : newMFreqType=inFreqFrame;
603 :
604 : MFrequency::Convert toNewFrame(inFreqFrame,
605 0 : MFrequency::Ref(newMFreqType, mframe_));
606 :
607 : // The velocity conversion engine:
608 0 : MDoppler::Ref dum1(MDoppler::RELATIVISTIC);
609 0 : MDoppler::Ref dum2(veldef);
610 0 : MDoppler::Convert dopConv(dum1, dum2);
611 :
612 : // Cope with unspecified rest freq
613 0 : MVFrequency rf=restFreq;
614 0 : if (restFreq.getValue()<=0.0)
615 0 : rf=toNewFrame(inFreq(inFreq.nelements()/2)).getValue();
616 :
617 : // Do the conversions
618 0 : for (uInt k=0; k< inFreq.nelements(); ++k){
619 0 : MDoppler eh = toNewFrame(inFreq(k)).toDoppler(rf);
620 0 : MDoppler eh2 = dopConv(eh);
621 0 : outVel(k)=eh2.getValue().get().getValue();
622 : }
623 :
624 :
625 :
626 :
627 0 : }
628 :
629 0 : MDirection VisBufferUtil::getPointingDir(const VisBuffer& vb, const Int antid, const Int vbrow){
630 :
631 0 : Timer tim;
632 0 : tim.mark();
633 0 : const MSColumns& msc=vb.msColumns();
634 : //cerr << "oldMSId_p " << oldMSId_p << " vb " << vb.msId() << endl;
635 0 : if(vb.msId() < 0)
636 0 : throw(AipsError("VisBuffer is not attached to an ms so cannot get pointing "));
637 0 : if(oldMSId_p != vb.msId()){
638 0 : oldMSId_p=vb.msId();
639 0 : if(timeAntIndex_p.shape()(0) < (oldMSId_p+1)){
640 0 : timeAntIndex_p.resize(oldMSId_p+1, true);
641 0 : cachedPointingDir_p.resize(oldMSId_p+1, true);
642 : }
643 0 : if( timeAntIndex_p[oldMSId_p].empty()){
644 0 : Vector<Double> tOrig;
645 0 : msc.time().getColumn(tOrig);
646 0 : Vector<Double> t;
647 0 : rejectConsecutive(tOrig, t);
648 0 : Vector<uInt> uniqIndx;
649 0 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, t, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
650 0 : uInt nAnt=msc.antenna().nrow();
651 0 : const MSPointingColumns& mspc=msc.pointing();
652 0 : Vector<Double> tUniq(nTimes);
653 0 : for (uInt k=0; k <nTimes; ++k){
654 0 : tUniq[k]= t[uniqIndx[k]];
655 : }
656 : Bool tstor, timcolstor, intcolstor, antcolstor;
657 0 : Double * tuniqptr=tUniq.getStorage(tstor);
658 0 : Int cshape=cachedPointingDir_p[oldMSId_p].shape()[0];
659 0 : cachedPointingDir_p[oldMSId_p].resize(cshape+nTimes*nAnt, True);
660 0 : Vector<Double> timecol;
661 0 : Vector<Double> intervalcol;
662 0 : Vector<Int> antcol;
663 0 : mspc.time().getColumn(timecol, True);
664 0 : mspc.interval().getColumn(intervalcol, True);
665 0 : mspc.antennaId().getColumn(antcol, True);
666 0 : Double *tcolptr=timecol.getStorage(timcolstor);
667 0 : Double *intcolptr=intervalcol.getStorage(intcolstor);
668 0 : Int * antcolptr=antcol.getStorage(antcolstor);
669 0 : Int npointrow=msc.pointing().nrow();
670 0 : #pragma omp parallel for firstprivate(nTimes, tuniqptr, tcolptr, antcolptr, intcolptr, npointrow), shared(mspc)
671 : for (uInt a=0; a < nAnt; ++a){
672 :
673 : //Double wtime1=omp_get_wtime();
674 : Vector<ssize_t> indices;
675 : Vector<MDirection> theDirs(nTimes);
676 : pointingIndex(tcolptr, antcolptr, intcolptr, npointrow, a, nTimes, tuniqptr, indices);
677 :
678 : #pragma omp critical
679 : {
680 : for (uInt k=0; k <nTimes; ++k){
681 :
682 :
683 : // std::ostringstream oss;
684 : // oss.precision(13);
685 : // oss << tuniqptr[k] << "_" << a;
686 : // String key=oss.str();
687 : std::pair<double, int> key=make_pair(t[uniqIndx[k]],a);
688 : timeAntIndex_p[oldMSId_p][key]=indices[k] > -1 ? cshape : -1;
689 :
690 : if(indices[k] >-1){
691 :
692 : cachedPointingDir_p[oldMSId_p][cshape]=mspc.directionMeas(indices[k]);
693 : cshape+=1;
694 : }
695 :
696 :
697 : }
698 : }//end critical
699 :
700 :
701 : }
702 :
703 0 : cachedPointingDir_p[oldMSId_p].resize(cshape, True);
704 : }
705 :
706 : }
707 :
708 : /////
709 : // String index=String::toString(vb.time()(vbrow))+String("_")+String::toString(antid);
710 : // std::ostringstream oss;
711 : // oss.precision(13);
712 : // oss << vb.time()(vbrow) << "_" << antid ;
713 : // String index=oss.str();
714 : //cerr << "index "<< index << endl;
715 : ////////////
716 : /*
717 : for (auto it = timeAntIndex_p[oldMSId_p].begin(); it != timeAntIndex_p[oldMSId_p].end(); ++it)
718 : {
719 : cerr << (*it).first << " --> " << (*it).second << endl;
720 : }
721 : */
722 : /////////////
723 0 : std::pair<double, int> index=make_pair(vb.time()(vbrow), antid);
724 0 : Int rowincache=timeAntIndex_p[oldMSId_p][index];
725 : ///////TESTOO
726 : /* if(rowincache>=0){
727 : cerr << "msid " << oldMSId_p << " key "<< index << " index " << rowincache<< " " << cachedPointingDir_p[oldMSId_p][rowincache] << endl;
728 : }*/
729 : /////////////
730 : //tim.show("retrieved cache");
731 0 : if(rowincache <0)
732 0 : return vb.phaseCenter();
733 0 : return cachedPointingDir_p[oldMSId_p][rowincache];
734 :
735 : }
736 0 : MDirection VisBufferUtil::getPointingDir(const vi::VisBuffer2& vb, const Int antid, const Int vbrow, const MDirection::Types dirframe, const Bool usePointing){
737 :
738 : //Double wtime0=omp_get_wtime();
739 0 : Int rowincache=-1;
740 0 : if(usePointing){
741 0 : if(oldMSId_p != vb.msId()){
742 0 : MSColumns msc(vb.ms());
743 0 : oldMSId_p=vb.msId();
744 0 : if(timeAntIndex_p.shape()(0) < (oldMSId_p+1)){
745 0 : timeAntIndex_p.resize(oldMSId_p+1, true);
746 0 : cachedPointingDir_p.resize(oldMSId_p+1, true);
747 : }
748 0 : MEpoch::Types timeType=MEpoch::castType(msc.timeMeas()(0).getRef().getType());
749 0 : Unit timeUnit(msc.timeMeas().measDesc().getUnits()(0).getName());
750 0 : MPosition pos;
751 0 : String tel;
752 0 : if (vb.subtableColumns().observation().nrow() > 0) {
753 0 : tel =vb.subtableColumns().observation().telescopeName()(msc.observationId()(0));
754 : }
755 0 : if (tel.length() == 0 || !tel.contains("VLA") ||
756 0 : !MeasTable::Observatory(pos,tel)) {
757 : // unknown observatory, use first antenna
758 0 : Int ant1=vb.antenna1()(0);
759 0 : pos=vb.subtableColumns().antenna().positionMeas()(ant1);
760 : }
761 0 : if( timeAntIndex_p[oldMSId_p].empty()){
762 0 : Vector<Double> tOrig;
763 0 : msc.time().getColumn(tOrig);
764 0 : Vector<Double> t;
765 0 : rejectConsecutive(tOrig, t);
766 0 : Vector<uInt> uniqIndx;
767 0 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, t, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
768 0 : uInt nAnt=msc.antenna().nrow();
769 0 : const MSPointingColumns& mspc=msc.pointing();
770 0 : Vector<Double> tUniq(nTimes);
771 0 : for (uInt k=0; k <nTimes; ++k){
772 0 : tUniq[k]= t[uniqIndx[k]];
773 : }
774 : Bool tstor, timcolstor, intcolstor, antcolstor;
775 0 : Double * tuniqptr=tUniq.getStorage(tstor);
776 0 : Int cshape=cachedPointingDir_p[oldMSId_p].shape()[0];
777 0 : cachedPointingDir_p[oldMSId_p].resize(cshape+nTimes*nAnt, True);
778 0 : Vector<Double> timecol;
779 0 : Vector<Double> intervalcol;
780 0 : Vector<Int> antcol;
781 0 : mspc.time().getColumn(timecol, True);
782 0 : mspc.interval().getColumn(intervalcol, True);
783 0 : mspc.antennaId().getColumn(antcol, True);
784 0 : Double *tcolptr=timecol.getStorage(timcolstor);
785 0 : Double *intcolptr=intervalcol.getStorage(intcolstor);
786 0 : Int * antcolptr=antcol.getStorage(antcolstor);
787 0 : Int npointrow=vb.ms().pointing().nrow();
788 : //ofstream myfile;
789 : //myfile.open ("POINTING.txt", ios::trunc);
790 0 : #pragma omp parallel for firstprivate(nTimes, tuniqptr, tcolptr, antcolptr, intcolptr, npointrow, timeType, timeUnit, pos, dirframe), shared(mspc)
791 : for (uInt a=0; a < nAnt; ++a){
792 :
793 : //Double wtime1=omp_get_wtime();
794 : Vector<ssize_t> indices;
795 : //Vector<MDirection> theDirs(nTimes);
796 : pointingIndex(tcolptr, antcolptr, intcolptr, npointrow, a, nTimes, tuniqptr, indices);
797 :
798 : #pragma omp critical
799 : {
800 : MEpoch timenow(Quantity(tuniqptr[0], timeUnit),timeType);
801 : MeasFrame mframe(timenow, pos);
802 : MDirection::Convert cvt(MDirection(), MDirection::Ref(dirframe, mframe));
803 : for (uInt k=0; k <nTimes; ++k){
804 :
805 :
806 : /*std::ostringstream oss;
807 : oss.precision(13);
808 : oss << tuniqptr[k] << "_" << a;
809 : String key=oss.str();
810 : */
811 : // myfile <<std::setprecision(13) << tuniqptr[k] << "_" << a << " index "<< indices[k] << "\n";
812 : pair<double, int> key=make_pair(tuniqptr[k],a);
813 : timeAntIndex_p[oldMSId_p][key]=indices[k] > -1 ? cshape : -1;
814 :
815 : if(indices[k] >-1){
816 : timenow=MEpoch(Quantity(tuniqptr[k], timeUnit),timeType);
817 : mframe.resetEpoch(timenow);
818 : cachedPointingDir_p[oldMSId_p][cshape]=cvt(mspc.directionMeas(indices[k]));
819 :
820 : cshape+=1;
821 : }
822 :
823 :
824 : }
825 : }//end critical
826 :
827 :
828 : }
829 0 : cachedPointingDir_p[oldMSId_p].resize(cshape, True);
830 : }
831 :
832 : }
833 :
834 : /////
835 : // String index=String::toString(vb.time()(vbrow))+String("_")+String::toString(antid);
836 : /* std::ostringstream oss;
837 : oss.precision(13);
838 : oss << vb.time()(vbrow) << "_" << antid ;
839 : String index=oss.str();
840 : */
841 0 : pair<double, int> index=make_pair(vb.time()(vbrow),antid);
842 0 : rowincache=timeAntIndex_p[oldMSId_p].at(index);
843 :
844 : //tim.show("retrieved cache");
845 : }///if usepointing
846 0 : if(rowincache <0)
847 0 : return getPhaseCenter(vb);
848 0 : return cachedPointingDir_p[oldMSId_p][rowincache];
849 :
850 :
851 :
852 : }
853 :
854 0 : void VisBufferUtil::pointingIndex(Double*& timecol, Int*& antcol, Double*& intervalcol, const rownr_t nrow, const Int ant, const Int ntimes, Double*& ptime, Vector<ssize_t>& indices){
855 :
856 0 : indices.resize(ntimes);
857 :
858 0 : indices.set(-1);
859 :
860 0 : for(Int pt=0; pt < ntimes; ++pt){
861 : //cerr << " " << guessRow ;
862 :
863 : /*for (Int k=0; k< 2; ++k){
864 : Int start=guessRow;
865 : Int end=nrow;
866 : if(k==1){
867 : start=0;
868 : end=guessRow;
869 : }
870 : */
871 0 : Double nearval=1e99;
872 0 : ssize_t nearestIndx=-1;
873 0 : ssize_t start=0;
874 0 : ssize_t end=nrow;
875 0 : for (ssize_t i=start; i<end; i++) {
876 0 : if(intervalcol[i]<=0.0 && ant==antcol[i]){
877 0 : if(abs(timecol[i]-ptime[pt]) < nearval){
878 0 : nearestIndx=i;
879 0 : nearval=abs(timecol[i]-ptime[pt]);
880 : }
881 : }
882 : }
883 0 : indices[pt]=nearestIndx;
884 :
885 :
886 0 : for (ssize_t i=start; i<end; i++) {
887 0 : if(ant == antcol[i]){
888 0 : Double halfInt=0.0;
889 0 : Bool done=False;
890 0 : if(intervalcol[i]<=0.0){
891 : // if(abs(timecol[inx[i]]-ptime[pt]) < nearval){
892 : // nearestIndx=inx[i];
893 : // }
894 0 : ssize_t counter=0;
895 0 : ssize_t adder=1;
896 0 : done=False;
897 : // while(!( (timecol[i+counter]!=timecol[i]))){
898 0 : while(!done){
899 0 : counter=counter+adder;
900 0 : if((ssize_t)nrow <= i+counter){
901 0 : adder=-1;
902 0 : counter=0;
903 : }
904 : ////Could not find another point (interval is infinite) hence only 1 valid point
905 0 : if( (i+counter) < 0){
906 0 : cerr << "HIT BREAK " << endl;
907 0 : indices[pt]=i;
908 0 : break;
909 : }
910 0 : if( (antcol[i+counter]==ant && timecol[i+counter] != timecol[i]) ){
911 0 : done=True;
912 : }
913 : }
914 :
915 : //if(ant==12 && abs(timecol[i+counter]-timecol[i]) > 10.0){
916 : //cerr << "i " << i << " counter " << counter << " done " << done << " adder " << adder << "ant count "<< antcol[i+counter] << " diff " << abs(timecol[i+counter]-timecol[i]) << endl;
917 : // }
918 0 : halfInt = abs(timecol[i+counter]-timecol[i])/2.0;
919 : }
920 : else{
921 0 : halfInt = intervalcol[i]/2.0;
922 : }
923 0 : if (halfInt>0.0) {
924 :
925 0 : if ((timecol[i] >= (ptime[pt] - halfInt)) && (timecol[i] <= (ptime[pt] + halfInt))) {
926 : ////TESTOO
927 : //if(ant==12){
928 : // cerr << "timecol " << timecol[i] << " halfInt " << halfInt << " TEST " << timecol[nearestIndx] << " inx " << i << " " << nearestIndx << endl;
929 : //}
930 0 : indices[pt]=abs(timecol[i]-ptime[pt]) < nearval ? i : nearestIndx;
931 : ////////TESTOO
932 0 : if(indices[pt] > 4688000){
933 0 : cerr << indices[pt] << " timecol " << timecol[i] << " halfInt " << halfInt << " TEST " << timecol[nearestIndx] << " nearval " << nearval << " inx " << i << " " << nearestIndx << endl;
934 :
935 : }
936 : ///////////////////
937 0 : break;
938 : }
939 : } else {
940 : // valid for all times (we should also handle interval<0 -> timestamps)
941 0 : cerr << "JUMPY " << i << " ant " << ant << " halfint " << halfInt << " done "<< done << endl;
942 0 : indices[pt]=i;
943 0 : break;
944 : }
945 :
946 : }//if ant
947 : }//start end
948 :
949 : //}//k
950 :
951 : }//pt
952 :
953 : //cerr << "ant " << ant << " indices " << indices << endl;
954 0 : }
955 :
956 0 : MDirection VisBufferUtil::getPhaseCenter(const vi::VisBuffer2& vb, const Double timeo){
957 : //Timer tim;
958 :
959 0 : Double timeph = timeo > 0.0 ? timeo : vb.time()(0);
960 : //MDirection outdir;
961 0 : if(oldPCMSId_p != vb.msId()){
962 0 : MSColumns msc(vb.ms());
963 : //tim.mark();
964 : //cerr << "MSID: "<< oldPCMSId_p << " " << vb.msId() << endl;
965 0 : oldPCMSId_p=vb.msId();
966 0 : if(cachedPhaseCenter_p.shape()(0) < (oldPCMSId_p+1)){
967 0 : cachedPhaseCenter_p.resize(oldPCMSId_p+1, true);
968 0 : cachedPhaseCenter_p[oldPCMSId_p]=map<Double, MDirection>();
969 : }
970 0 : if( cachedPhaseCenter_p[oldPCMSId_p].empty()){
971 0 : Vector<Double> tOrig;
972 0 : msc.time().getColumn(tOrig);
973 0 : Vector<Int> fieldId;
974 0 : msc.fieldId().getColumn(fieldId);
975 0 : Vector<Double> t;
976 0 : Vector<Int> origindx;
977 0 : rejectConsecutive(tOrig, t, origindx);
978 0 : Vector<uInt> uniqIndx;
979 :
980 0 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, t, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
981 : //cerr << "ntimes " << nTimes << " " << uniqIndx << "\n origInx " << origindx << endl;
982 0 : const MSFieldColumns& msfc=msc.field();
983 0 : for (uInt k=0; k <nTimes; ++k){
984 : //cerr << t[uniqIndx[k]] << " " << fieldId[origindx[uniqIndx[k]]] << endl;
985 : //cerr << msfc.phaseDirMeas(fieldId[origindx[uniqIndx[k]]], t[uniqIndx[k]]) << endl;
986 : //cerr << "size " << cachedPhaseCenter_p[oldPCMSId_p].size() << endl;
987 0 : String ephemIfAny=msfc.ephemPath(fieldId[origindx[uniqIndx[k]]]);
988 0 : if(ephemIfAny=="" || !Table::isReadable(ephemIfAny, False)){
989 0 : (cachedPhaseCenter_p[oldPCMSId_p])[t[uniqIndx[k]]]=msfc.phaseDirMeas(fieldId[origindx[uniqIndx[k]]], t[uniqIndx[k]]);
990 : }
991 : else{
992 0 : Vector<MDirection> refDir(msfc.referenceDirMeasCol()(fieldId[origindx[uniqIndx[k]]]));
993 0 : (cachedPhaseCenter_p[oldPCMSId_p])[t[uniqIndx[k]]]=getEphemBasedPhaseDir(vb, ephemIfAny, refDir(0), t[uniqIndx[k]]);
994 : }
995 :
996 : }
997 :
998 :
999 : }
1000 : //tim.show("After caching all phasecenters");
1001 : }
1002 : //tim.mark();
1003 0 : MDirection retval;
1004 0 : auto it=cachedPhaseCenter_p[oldPCMSId_p].find(timeph);
1005 0 : if(it != cachedPhaseCenter_p[oldPCMSId_p].end()){
1006 0 : retval=it->second;
1007 : }
1008 : else{
1009 0 : auto upp= cachedPhaseCenter_p[oldPCMSId_p].upper_bound(timeph);
1010 0 : auto low= cachedPhaseCenter_p[oldPCMSId_p].lower_bound(timeph);
1011 0 : if (upp==cachedPhaseCenter_p[oldPCMSId_p].begin())
1012 0 : retval=(cachedPhaseCenter_p[oldPCMSId_p].begin())->second;
1013 0 : else if(low==cachedPhaseCenter_p[oldPCMSId_p].end()){
1014 0 : --low;
1015 0 : retval=low->second;
1016 : }
1017 : else{
1018 0 : if(fabs(timeph - (low->first)) < fabs(timeph - (upp->first))){
1019 0 : retval=low->second;
1020 : }
1021 : else{
1022 0 : retval=upp->second;
1023 : }
1024 :
1025 : }
1026 : }
1027 : //tim.show("retrieved cache");
1028 : //cerr << std::setprecision(12) <<"msid " << oldPCMSId_p << " time "<< timeph << " val " << retval.toString() << endl;
1029 :
1030 0 : return retval;
1031 :
1032 :
1033 :
1034 : }
1035 :
1036 :
1037 0 : MDirection VisBufferUtil::getEphemBasedPhaseDir(const vi::VisBuffer2& vb, const String& ephemPath, const MDirection&refDir, const Double t){
1038 0 : MEpoch ep(Quantity(t, "s"), vb.getVi()->getImpl()->getEpoch().getRef());
1039 0 : mframe_.resetEpoch(ep);
1040 0 : if(!Table::isReadable(ephemPath, False))
1041 0 : return refDir;
1042 0 : MeasComet mcomet(Path(ephemPath).absoluteName());
1043 0 : mframe_.set(mcomet);
1044 0 : MDirection::Ref outref1(MDirection::AZEL, mframe_);
1045 0 : MDirection tmpazel=MDirection::Convert(MDirection(MDirection::COMET), outref1)();
1046 0 : MDirection::Types outtype=(MDirection::Types) refDir.getRef().getType();
1047 0 : MDirection::Ref outref(outtype, mframe_);
1048 0 : MDirection outdir=MDirection::Convert(tmpazel, outref)();
1049 0 : MVDirection mvoutdir(outdir.getAngle());
1050 0 : MVDirection mvrefdir(refDir.getAngle());
1051 : //copying what ROMSFieldColumns::extractDirMeas does
1052 0 : mvoutdir.shift(mvrefdir.getAngle(), True);
1053 0 : return MDirection(mvoutdir, outtype);
1054 : }
1055 :
1056 0 : MDirection VisBufferUtil::getEphemDir(const vi::VisBuffer2& vb,
1057 : const Double timeo){
1058 :
1059 0 : Double timeEphem = timeo > 0.0 ? timeo : vb.time()(0);
1060 0 : const MSFieldColumns &msfc = vb.subtableColumns().field();
1061 0 : Int fieldId=vb.fieldId()(0);
1062 0 : MDirection refDir(Quantity(0, "deg"), Quantity(0, "deg"),(MDirection::Types)msfc.ephemerisDirMeas(fieldId, timeEphem).getRef().getType());
1063 : //Now do the parallax correction
1064 0 : return getEphemBasedPhaseDir(vb, msfc.ephemPath(fieldId), refDir, timeEphem);
1065 :
1066 :
1067 : }
1068 :
1069 : //utility to reject consecutive similar value for sorting
1070 0 : void VisBufferUtil::rejectConsecutive(const Vector<Double>& t, Vector<Double>& retval){
1071 0 : uInt n=t.nelements();
1072 0 : if(n >0){
1073 0 : retval.resize(n);
1074 0 : retval[0]=t[0];
1075 : }
1076 : else
1077 0 : return;
1078 0 : Int prev=0;
1079 0 : for (uInt k=1; k < n; ++k){
1080 0 : if(t[k] != retval(prev)){
1081 0 : ++prev;
1082 :
1083 0 : retval[prev]=t[k];
1084 : }
1085 : }
1086 0 : retval.resize(prev+1, true);
1087 :
1088 : }
1089 0 : void VisBufferUtil::rejectConsecutive(const Vector<Double>& t, Vector<Double>& retval, Vector<Int>& indx){
1090 0 : uInt n=t.nelements();
1091 0 : if(n >0){
1092 0 : retval.resize(n);
1093 0 : indx.resize(n);
1094 0 : retval[0]=t[0];
1095 0 : indx[0]=0;
1096 : }
1097 : else
1098 0 : return;
1099 0 : Int prev=0;
1100 0 : for (uInt k=1; k < n; ++k){
1101 0 : if(t[k] != retval(prev)){
1102 0 : ++prev;
1103 : //retval.resize(prev+1, true);
1104 0 : retval[prev]=t[k];
1105 : //indx.resize(prev+1, true);
1106 0 : indx[prev]=k;
1107 : }
1108 : }
1109 0 : retval.resize(prev+1, true);
1110 0 : indx.resize(prev+1, true);
1111 :
1112 : }
1113 :
1114 : // helper function to swap the y and z axes of a Cube
1115 0 : void VisBufferUtil::swapyz(Cube<Complex>& out, const Cube<Complex>& in)
1116 : {
1117 0 : IPosition inShape=in.shape();
1118 0 : uInt nxx=inShape(0),nyy=inShape(2),nzz=inShape(1);
1119 : //resize breaks references...so out better have the right shape
1120 : //if references is not to be broken
1121 0 : if(out.nelements()==0)
1122 0 : out.resize(nxx,nyy,nzz);
1123 : Bool deleteIn,deleteOut;
1124 0 : const Complex* pin = in.getStorage(deleteIn);
1125 0 : Complex* pout = out.getStorage(deleteOut);
1126 0 : uInt i=0, zOffset=0;
1127 0 : for (uInt iz=0; iz<nzz; ++iz, zOffset+=nxx) {
1128 0 : Int yOffset=zOffset;
1129 0 : for (uInt iy=0; iy<nyy; ++iy, yOffset+=nxx*nzz) {
1130 0 : for (uInt ix=0; ix<nxx; ++ix){
1131 0 : pout[i++] = pin[ix+yOffset];
1132 : }
1133 : }
1134 : }
1135 0 : out.putStorage(pout,deleteOut);
1136 0 : in.freeStorage(pin,deleteIn);
1137 0 : }
1138 :
1139 : // helper function to swap the y and z axes of a Cube
1140 0 : void VisBufferUtil::swapyz(Cube<Bool>& out, const Cube<Bool>& in)
1141 : {
1142 0 : IPosition inShape=in.shape();
1143 0 : uInt nxx=inShape(0),nyy=inShape(2),nzz=inShape(1);
1144 0 : if(out.nelements()==0)
1145 0 : out.resize(nxx,nyy,nzz);
1146 : Bool deleteIn,deleteOut;
1147 0 : const Bool* pin = in.getStorage(deleteIn);
1148 0 : Bool* pout = out.getStorage(deleteOut);
1149 0 : uInt i=0, zOffset=0;
1150 0 : for (uInt iz=0; iz<nzz; iz++, zOffset+=nxx) {
1151 0 : Int yOffset=zOffset;
1152 0 : for (uInt iy=0; iy<nyy; iy++, yOffset+=nxx*nzz) {
1153 0 : for (uInt ix=0; ix<nxx; ix++) pout[i++] = pin[ix+yOffset];
1154 : }
1155 : }
1156 0 : out.putStorage(pout,deleteOut);
1157 0 : in.freeStorage(pin,deleteIn);
1158 0 : }
1159 :
1160 :
1161 :
1162 :
1163 : } //# NAMESPACE CASA - END
1164 :
|
