Line data Source code
1 : //# PhaseShiftingTVI.h: This file contains the implementation of the PhaseShiftingTVI class.
2 : //#
3 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
4 : //# Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All rights reserved.
5 : //# Copyright (C) European Southern Observatory, 2011, All rights reserved.
6 : //#
7 : //# This library is free software; you can redistribute it and/or
8 : //# modify it under the terms of the GNU Lesser General Public
9 : //# License as published by the Free software Foundation; either
10 : //# version 2.1 of the License, or (at your option) any later version.
11 : //#
12 : //# This library is distributed in the hope that it will be useful,
13 : //# but WITHOUT ANY WARRANTY, without even the implied warranty of
14 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 : //# Lesser General Public License for more details.
16 : //#
17 : //# You should have received a copy of the GNU Lesser General Public
18 : //# License along with this library; if not, write to the Free Software
19 : //# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
20 : //# MA 02111-1307 USA
21 : //# $Id: $
22 :
23 : #include <mstransform/TVI/PhaseShiftingTVI.h>
24 :
25 : using namespace casacore;
26 : namespace casa { //# NAMESPACE CASA - BEGIN
27 :
28 : namespace vi { //# NAMESPACE VI - BEGIN
29 :
30 : //////////////////////////////////////////////////////////////////////////
31 : // PhaseShiftingTVI class
32 : //////////////////////////////////////////////////////////////////////////
33 :
34 : // -----------------------------------------------------------------------
35 : //
36 : // -----------------------------------------------------------------------
37 0 : PhaseShiftingTVI::PhaseShiftingTVI( ViImplementation2 * inputVii,
38 0 : const Record &configuration):
39 0 : FreqAxisTVI (inputVii)
40 : {
41 0 : dx_p = 0;
42 0 : dy_p = 0;
43 :
44 : // CAS-12706 Zero initialization for wide-field phase shifting algorithm
45 0 : wideFieldMode_p = false;
46 0 : phaseCenterName_p = "";
47 0 : selectedInputMsCols_p = NULL;
48 :
49 : // Parse and check configuration parameters
50 : // Note: if a constructor finishes by throwing an exception, the memory
51 : // associated with the object itself is cleaned up — there is no memory leak.
52 0 : if (not parseConfiguration(configuration))
53 : {
54 0 : throw AipsError("Error parsing PhaseShiftingTVI configuration");
55 : }
56 :
57 0 : initialize();
58 :
59 0 : return;
60 : }
61 :
62 : // -----------------------------------------------------------------------
63 : //
64 : // -----------------------------------------------------------------------
65 0 : Bool PhaseShiftingTVI::parseConfiguration(const Record &configuration)
66 : {
67 0 : int exists = -1;
68 0 : Bool ret = true;
69 :
70 0 : exists = -1;
71 0 : exists = configuration.fieldNumber ("XpcOffset");
72 0 : if (exists >= 0)
73 : {
74 0 : configuration.get (exists, dx_p);
75 : }
76 :
77 0 : exists = -1;
78 0 : exists = configuration.fieldNumber ("YpcOffset");
79 0 : if (exists >= 0)
80 : {
81 0 : configuration.get (exists, dy_p);
82 : }
83 :
84 0 : if (abs(dx_p) > 0 or abs(dy_p) > 0)
85 : {
86 0 : logger_p << LogIO::NORMAL << LogOrigin("PhaseShiftingTVI", __FUNCTION__)
87 0 : << "Phase shift is dx="<< dx_p << " dy=" << dy_p << LogIO::POST;
88 : }
89 :
90 : // CAS-12706 Add support for shifting across large offset/angles
91 0 : exists = -1;
92 0 : exists = configuration.fieldNumber ("phasecenter");
93 0 : if (exists >= 0)
94 : {
95 0 : configuration.get (exists, phaseCenterName_p);
96 : // casaMDirection requires a variant
97 0 : casac::variant phaseCenterVar(phaseCenterName_p);
98 :
99 0 : if(!casaMDirection(phaseCenterVar, phaseCenter_p))
100 : {
101 0 : logger_p << LogIO::SEVERE << LogOrigin("PhaseShiftingTVI", __FUNCTION__)
102 0 : << "Cannot interpret phase center " << phaseCenterName_p << LogIO::POST;
103 0 : ret = false;
104 : }
105 : else
106 : {
107 : MDirection::Types mdtype;
108 0 : const auto myFrame = phaseCenter_p.getRefString();
109 0 : MDirection::getType(mdtype, myFrame);
110 0 : ThrowIf(
111 : mdtype == MDirection::HADEC || mdtype == MDirection::AZEL
112 : || mdtype == MDirection::AZELSW || mdtype == MDirection::AZELNE
113 : || mdtype == MDirection::AZELGEO || mdtype == MDirection::AZELSWGEO
114 : || mdtype == MDirection::MECLIPTIC || mdtype == MDirection::TECLIPTIC
115 : || mdtype == MDirection::TOPO,
116 : myFrame + " is a time dependent reference frame and so is not supported"
117 : );
118 0 : ThrowIf(
119 : mdtype == MDirection::MERCURY || mdtype == MDirection::VENUS
120 : || mdtype == MDirection::MARS || mdtype == MDirection::JUPITER
121 : || mdtype == MDirection::SATURN || mdtype == MDirection::URANUS
122 : || mdtype == MDirection::NEPTUNE || mdtype == MDirection::PLUTO
123 : || mdtype == MDirection::SUN || mdtype == MDirection::MOON
124 : || mdtype == MDirection::COMET,
125 : myFrame + " denotes an ephemeris object and so is not supported"
126 : );
127 0 : wideFieldMode_p = true;
128 0 : logger_p << LogIO::NORMAL << LogOrigin("PhaseShiftingTVI", __FUNCTION__)
129 0 : << "Phase center " << phaseCenterName_p << " successfully parsed"<< LogIO::POST;
130 : }
131 : }
132 :
133 0 : return ret;
134 : }
135 :
136 : // -----------------------------------------------------------------------
137 : //
138 : // -----------------------------------------------------------------------
139 0 : void PhaseShiftingTVI::initialize()
140 : {
141 : // Populate nchan input-output maps
142 : Int spw;
143 0 : uInt spw_idx = 0;
144 0 : map<Int,vector<Int> >::iterator iter;
145 0 : for(iter=spwInpChanIdxMap_p.begin();iter!=spwInpChanIdxMap_p.end();iter++)
146 : {
147 0 : spw = iter->first;
148 0 : spwOutChanNumMap_p[spw] = spwInpChanIdxMap_p[spw].size();
149 :
150 0 : spw_idx++;
151 : }
152 :
153 : // CAS-12706 Add support for shifting across large offset/angles
154 : // Access observatory position and observation start (reference) time.
155 0 : if (wideFieldMode_p)
156 : {
157 0 : selectedInputMsCols_p = new MSColumns(getVii()->ms());
158 0 : observatoryPosition_p = selectedInputMsCols_p->antenna().positionMeas()(0);
159 0 : referenceTime_p = selectedInputMsCols_p->timeMeas()(0);
160 0 : referenceTimeUnits_p = selectedInputMsCols_p->timeQuant()(0).getUnit();
161 : }
162 :
163 0 : return;
164 : }
165 :
166 : // -----------------------------------------------------------------------
167 : // CAS-12706 Recalculate UVW and Phases according to wide-field phase shifting
168 : // algorithm. This method should do the same as ComponentFTMachine::rotateUVW
169 : // -----------------------------------------------------------------------
170 0 : void PhaseShiftingTVI::shiftUVWPhases()
171 : {
172 : // Get input VisBuffer
173 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
174 :
175 0 : auto convertedPhaseCenter = phaseCenter_p;
176 0 : if (phaseCenter_p.getRefString() != vb->phaseCenter().getRefString()) {
177 : MDirection::Types mdtype;
178 0 : MDirection::getType(mdtype, vb->phaseCenter().getRefString());
179 0 : convertedPhaseCenter = MDirection::Convert(phaseCenter_p, mdtype)();
180 : }
181 :
182 : // Initialize epoch corresponding to current buffer
183 : // with time reference to the first row in the MS
184 0 : MEpoch epoch(Quantity(vb->time()(0),referenceTimeUnits_p),referenceTime_p.getRef());
185 0 : MeasFrame refFrame(epoch,observatoryPosition_p);
186 0 : UVWMachine uvwMachine(convertedPhaseCenter, vb->phaseCenter(), refFrame,false,false);
187 : // Initialize phase array and uvw matrix
188 0 : phaseShift_p.resize(vb->nRows(),false);
189 0 : newUVW_p.resize(vb->uvw().shape(),false);
190 :
191 : // Obtain phase shift and new uvw coordinates
192 0 : Vector<Double> dummy(3,0.0);
193 0 : double phase2radPerHz = -2.0 * C::pi / C::c;
194 0 : for (uInt row=0;row<vb->nRows();row++)
195 : {
196 : // Copy current uvw coordinates so that they are not modified
197 : // Note: Columns in uvw correspond to rows in the main table/VisBuffer!
198 0 : dummy = vb->uvw().column(row);
199 :
200 : // Have to change (u,v,w) to (-u,-v,w) because is the convention used by uvwMachine
201 0 : dummy(0) = -1*dummy(0);
202 0 : dummy(1) = -1*dummy(1);
203 :
204 : // Transform uvw coordinates and obtain corresponding phase shift
205 0 : uvwMachine.convertUVW(phaseShift_p(row), dummy);
206 :
207 : // Store new uvw coordinates
208 : // Have to change back (-u,-v,w) to (u,v,w) because is the convention used by the MS
209 0 : dummy(0) = -1*dummy(0);
210 0 : dummy(1) = -1*dummy(1);
211 0 : newUVW_p.column(row) = dummy;
212 : // Convert phase shift to radian/Hz
213 0 : phaseShift_p(row) = phase2radPerHz*phaseShift_p(row);
214 : }
215 0 : return;
216 : }
217 :
218 : // -----------------------------------------------------------------------
219 : //
220 : // -----------------------------------------------------------------------
221 0 : void PhaseShiftingTVI::origin()
222 : {
223 : // Drive underlying ViImplementation2
224 0 : getVii()->origin();
225 :
226 : // CAS-12706 Add support for shifting across large offset/angles
227 0 : if (wideFieldMode_p) shiftUVWPhases();
228 :
229 : // Define the shapes in the VB2, patch provided by cgarcia in CAS-12706
230 0 : configureShapes();
231 :
232 : // Synchronize own VisBuffer
233 0 : configureNewSubchunk();
234 :
235 0 : return;
236 : }
237 :
238 : // -----------------------------------------------------------------------
239 : //
240 : // -----------------------------------------------------------------------
241 0 : void PhaseShiftingTVI::next()
242 : {
243 : // Drive underlying ViImplementation2
244 0 : getVii()->next();
245 :
246 : // CAS-12706 Add support for shifting across large offset/angles
247 0 : if (wideFieldMode_p) shiftUVWPhases();
248 :
249 : // Define the shapes in the VB2, patch provided by cgarcia in CAS-12706
250 0 : configureShapes();
251 :
252 : // Synchronize own VisBuffer
253 0 : configureNewSubchunk();
254 :
255 0 : return;
256 : }
257 :
258 :
259 : // -----------------------------------------------------------------------
260 : //
261 : // -----------------------------------------------------------------------
262 0 : void PhaseShiftingTVI::visibilityObserved (Cube<Complex> & vis) const
263 : {
264 : // Get input VisBuffer
265 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
266 0 : Matrix<Double> uvw = vb->uvw();
267 0 : Vector<Double> frequencies = vb->getFrequencies(0);
268 :
269 : // Reshape output data before passing it to the DataCubeHolder
270 0 : vis.resize(getVisBuffer()->getShape(),false);
271 :
272 : // Gather input data
273 0 : DataCubeMap inputData;
274 0 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCube());
275 0 : inputData.add(MS::DATA,inputVisCubeHolder);
276 :
277 : // Gather output data
278 0 : DataCubeMap outputData;
279 0 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
280 0 : outputData.add(MS::DATA,outputVisCubeHolder);
281 :
282 0 : if (wideFieldMode_p)
283 : {
284 : // Configure Transformation Engine
285 0 : WideFieldPhaseShiftingTransformEngine<Complex> transformer(phaseShift_p,&uvw,&frequencies,&inputData,&outputData);
286 :
287 : // Transform data
288 0 : transformFreqAxis2(vb->getShape(),transformer);
289 : }
290 : else
291 : {
292 : // Configure Transformation Engine
293 0 : PhaseShiftingTransformEngine<Complex> transformer(dx_p,dy_p,&uvw,&frequencies,&inputData,&outputData);
294 :
295 : // Transform data
296 0 : transformFreqAxis2(vb->getShape(),transformer);
297 : }
298 :
299 0 : return;
300 : }
301 :
302 : // -----------------------------------------------------------------------
303 : //
304 : // -----------------------------------------------------------------------
305 0 : void PhaseShiftingTVI::visibilityCorrected (Cube<Complex> & vis) const
306 : {
307 : // Get input VisBuffer
308 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
309 0 : Matrix<Double> uvw = vb->uvw();
310 0 : Vector<Double> frequencies = vb->getFrequencies(0);
311 :
312 : // Reshape output data before passing it to the DataCubeHolder
313 0 : vis.resize(getVisBuffer()->getShape(),false);
314 :
315 : // Gather input data
316 0 : DataCubeMap inputData;
317 0 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCubeCorrected());
318 0 : inputData.add(MS::DATA,inputVisCubeHolder);
319 :
320 : // Gather output data
321 0 : DataCubeMap outputData;
322 0 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
323 0 : outputData.add(MS::DATA,outputVisCubeHolder);
324 :
325 0 : if (wideFieldMode_p)
326 : {
327 : // Configure Transformation Engine
328 0 : WideFieldPhaseShiftingTransformEngine<Complex> transformer(phaseShift_p,&uvw,&frequencies,&inputData,&outputData);
329 :
330 : // Transform data
331 0 : transformFreqAxis2(vb->getShape(),transformer);
332 : }
333 : else
334 : {
335 : // Configure Transformation Engine
336 0 : PhaseShiftingTransformEngine<Complex> transformer(dx_p,dy_p,&uvw,&frequencies,&inputData,&outputData);
337 :
338 : // Transform data
339 0 : transformFreqAxis2(vb->getShape(),transformer);
340 : }
341 :
342 0 : return;
343 : }
344 :
345 : // -----------------------------------------------------------------------
346 : //
347 : // -----------------------------------------------------------------------
348 0 : void PhaseShiftingTVI::visibilityModel (Cube<Complex> & vis) const
349 : {
350 : // Get input VisBuffer
351 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
352 0 : Matrix<Double> uvw = vb->uvw();
353 0 : Vector<Double> frequencies = vb->getFrequencies(0);
354 :
355 : // Reshape output data before passing it to the DataCubeHolder
356 0 : vis.resize(getVisBuffer()->getShape(),false);
357 :
358 : // Gather input data
359 0 : DataCubeMap inputData;
360 0 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCubeModel());
361 0 : inputData.add(MS::DATA,inputVisCubeHolder);
362 :
363 : // Gather output data
364 0 : DataCubeMap outputData;
365 0 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
366 0 : outputData.add(MS::DATA,outputVisCubeHolder);
367 :
368 0 : if (wideFieldMode_p)
369 : {
370 : // Configure Transformation Engine
371 0 : WideFieldPhaseShiftingTransformEngine<Complex> transformer(phaseShift_p,&uvw,&frequencies,&inputData,&outputData);
372 :
373 : // Transform data
374 0 : transformFreqAxis2(vb->getShape(),transformer);
375 : }
376 : else
377 : {
378 : // Configure Transformation Engine
379 0 : PhaseShiftingTransformEngine<Complex> transformer(dx_p,dy_p,&uvw,&frequencies,&inputData,&outputData);
380 :
381 : // Transform data
382 0 : transformFreqAxis2(vb->getShape(),transformer);
383 : }
384 :
385 0 : return;
386 : }
387 :
388 : // -----------------------------------------------------------------------
389 : //
390 : // -----------------------------------------------------------------------
391 0 : void PhaseShiftingTVI::uvw (casacore::Matrix<double> & uvw) const
392 : {
393 0 : if (wideFieldMode_p)
394 : {
395 0 : uvw.resize(newUVW_p.shape(),false);
396 0 : uvw = newUVW_p;
397 : }
398 : else
399 : {
400 0 : getVii()->uvw (uvw);
401 : }
402 :
403 0 : return;
404 : }
405 :
406 : //////////////////////////////////////////////////////////////////////////
407 : // PhaseShiftingTVIFactory class
408 : //////////////////////////////////////////////////////////////////////////
409 :
410 : // -----------------------------------------------------------------------
411 : //
412 : // -----------------------------------------------------------------------
413 0 : PhaseShiftingTVIFactory::PhaseShiftingTVIFactory ( Record &configuration,
414 0 : ViImplementation2 *inputVii)
415 : {
416 0 : inputVii_p = inputVii;
417 0 : configuration_p = configuration;
418 0 : }
419 :
420 : // -----------------------------------------------------------------------
421 : //
422 : // -----------------------------------------------------------------------
423 0 : vi::ViImplementation2 * PhaseShiftingTVIFactory::createVi(VisibilityIterator2 *) const
424 : {
425 0 : return new PhaseShiftingTVI(inputVii_p,configuration_p);
426 : }
427 :
428 : // -----------------------------------------------------------------------
429 : //
430 : // -----------------------------------------------------------------------
431 0 : vi::ViImplementation2 * PhaseShiftingTVIFactory::createVi() const
432 : {
433 0 : return new PhaseShiftingTVI(inputVii_p,configuration_p);
434 : }
435 :
436 : //////////////////////////////////////////////////////////////////////////
437 : // PhaseShiftingTVILayerFactory class
438 : //////////////////////////////////////////////////////////////////////////
439 :
440 0 : PhaseShiftingTVILayerFactory::PhaseShiftingTVILayerFactory(Record &configuration) :
441 : ViiLayerFactory(),
442 0 : configuration_p(configuration)
443 0 : {}
444 :
445 : ViImplementation2*
446 0 : PhaseShiftingTVILayerFactory::createInstance(ViImplementation2* vii0) const
447 : {
448 : // Make the PhaseShiftingTVi2, using supplied ViImplementation2, and return it
449 0 : ViImplementation2 *vii = new PhaseShiftingTVI(vii0,configuration_p);
450 0 : return vii;
451 : }
452 :
453 : //////////////////////////////////////////////////////////////////////////
454 : // PhaseShiftingTransformEngine class
455 : //////////////////////////////////////////////////////////////////////////
456 :
457 : // -----------------------------------------------------------------------
458 : //
459 : // -----------------------------------------------------------------------
460 0 : template<class T> PhaseShiftingTransformEngine<T>::PhaseShiftingTransformEngine(
461 : Double dx, Double dy,
462 : Matrix<Double> *uvw,
463 : Vector<Double> *frequencies,
464 : DataCubeMap *inputData,
465 : DataCubeMap *outputData):
466 0 : FreqAxisTransformEngine2<T>(inputData,outputData)
467 : {
468 0 : uvw_p = uvw;
469 0 : frequencies_p = frequencies;
470 :
471 : // Offsets in radians (input is arcsec)
472 0 : dx_p = dx*(C::pi / 180.0 / 3600.0);
473 0 : dy_p = dy*(C::pi / 180.0 / 3600.0);
474 0 : }
475 :
476 : // -----------------------------------------------------------------------
477 : //
478 : // -----------------------------------------------------------------------
479 0 : template<class T> void PhaseShiftingTransformEngine<T>::transform( )
480 : {
481 0 : transformCore(inputData_p,outputData_p);
482 0 : }
483 :
484 : // -----------------------------------------------------------------------
485 : //
486 : // -----------------------------------------------------------------------
487 0 : template<class T> void PhaseShiftingTransformEngine<T>::transformCore( DataCubeMap *inputData,
488 : DataCubeMap *outputData)
489 : {
490 : // Get input/output data
491 0 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
492 0 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
493 :
494 : // Extra path as fraction of U and V in m
495 0 : Double phase = dx_p*(*uvw_p)(0,rowIndex_p) + dy_p*(*uvw_p)(1,rowIndex_p);
496 :
497 : // In radian/Hz
498 0 : phase *= -2.0 * C::pi / C::c;
499 :
500 : // Main loop
501 : Double phase_i;
502 0 : Complex factor;
503 0 : for (uInt chan_i=0;chan_i<inputVector.size();chan_i++)
504 : {
505 0 : phase_i = phase * (*frequencies_p)(chan_i);
506 0 : factor = Complex(cos(phase_i), sin(phase_i));
507 0 : outputVector(chan_i) = factor*inputVector(chan_i);
508 : }
509 0 : }
510 :
511 : //////////////////////////////////////////////////////////////////////////
512 : // WideFieldPhaseShiftingTransformEngine class
513 : //////////////////////////////////////////////////////////////////////////
514 :
515 : // -----------------------------------------------------------------------
516 : //
517 : // -----------------------------------------------------------------------
518 0 : template<class T> WideFieldPhaseShiftingTransformEngine<T>::WideFieldPhaseShiftingTransformEngine(
519 : const Vector<Double> &phaseShift,
520 : Matrix<Double> *uvw,
521 : Vector<Double> *frequencies,
522 : DataCubeMap *inputData,
523 : DataCubeMap *outputData):
524 0 : FreqAxisTransformEngine2<T>(inputData,outputData)
525 : {
526 0 : uvw_p = uvw;
527 0 : frequencies_p = frequencies;
528 0 : phaseShift_p = phaseShift;
529 0 : }
530 :
531 : // -----------------------------------------------------------------------
532 : //
533 : // -----------------------------------------------------------------------
534 0 : template<class T> void WideFieldPhaseShiftingTransformEngine<T>::transform( )
535 : {
536 0 : transformCore(inputData_p,outputData_p);
537 0 : }
538 :
539 : // -----------------------------------------------------------------------
540 : //
541 : // -----------------------------------------------------------------------
542 0 : template<class T> void WideFieldPhaseShiftingTransformEngine<T>::transformCore( DataCubeMap *inputData,
543 : DataCubeMap *outputData)
544 : {
545 : // Get input/output data
546 0 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
547 0 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
548 :
549 : // Main loop
550 : Double phase_i;
551 0 : Complex factor;
552 0 : for (uInt chan_i=0;chan_i<inputVector.size();chan_i++)
553 : {
554 0 : phase_i = phaseShift_p(rowIndex_p) * (*frequencies_p)(chan_i);
555 0 : factor = Complex(cos(phase_i), sin(phase_i));
556 0 : outputVector(chan_i) = factor*inputVector(chan_i);
557 : }
558 0 : }
559 :
560 :
561 : } //# NAMESPACE VI - END
562 :
563 : } //# NAMESPACE CASA - END
564 :
565 :
|
