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