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-2024, 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 87 : PhaseShiftingTVI::PhaseShiftingTVI( ViImplementation2 * inputVii,
38 87 : const Record &configuration):
39 87 : FreqAxisTVI (inputVii)
40 : {
41 87 : dx_p = 0;
42 87 : dy_p = 0;
43 :
44 : // CAS-12706 Zero initialization for wide-field phase shifting algorithm
45 87 : wideFieldMode_p = false;
46 :
47 : // Parse and check configuration parameters
48 : // Note: if a constructor finishes by throwing an exception, the memory
49 : // associated with the object itself is cleaned up — there is no memory leak.
50 87 : parseConfiguration(configuration);
51 :
52 86 : initialize();
53 :
54 86 : return;
55 7 : }
56 :
57 : /**
58 : * Get max valid FIELD IDs for this iterator. MSv2 uses the FIELD
59 : * table row index as FIELD ID.
60 : */
61 5 : rownr_t PhaseShiftingTVI::getMaxMSFieldID() const
62 : {
63 5 : const auto &fieldsTable = getVii()->fieldSubtablecols();
64 5 : return fieldsTable.nrow() - 1;
65 : }
66 :
67 : /**
68 : * Parses the phase center parameter, whether given as a single string or
69 : * as a dictionary/record of per-field strings.
70 : * Populates the phaseCenterSpec_p map.
71 : *
72 : * @param config TVI configuration object
73 : */
74 87 : void PhaseShiftingTVI::parsePhasecenter(const Record &config)
75 : {
76 87 : auto exists = config.fieldNumber ("phasecenter");
77 87 : if (exists < 0) {
78 0 : return;
79 : }
80 :
81 : // phasecenter can be given as a string or as a dict (per-field centers)
82 87 : auto isStr = false;
83 87 : casacore::String phaseCenterStr;
84 : try {
85 92 : config.get(exists, phaseCenterStr);
86 82 : isStr = true;
87 164 : logger_p << LogIO::NORMAL << LogOrigin("PhaseShiftingTVI", __FUNCTION__)
88 : << "Interpreting phasecenter as a single string, to be applied"
89 164 : " to all input (selected) fields."<< LogIO::POST;
90 5 : } catch(const AipsError &exc) {
91 : // ignore "RecordRep::get_pointer - incorrect data type String used...
92 : // for field phasecenter with type Record" and similar, try record type
93 10 : logger_p << LogIO::NORMAL << LogOrigin("PhaseShiftingTVI", __FUNCTION__)
94 : << "Interpreting phasecenter as a dictionary of field->center."
95 10 : << "Using following phase centers:" << LogIO::POST;
96 5 : }
97 :
98 87 : if (isStr) {
99 82 : const auto &phaseCenterDir = checkPhaseCenterStr(phaseCenterStr);
100 82 : phaseCenterSpec_p.insert({-1, phaseCenterDir});
101 82 : } else {
102 5 : parsePhasecenterDict(config);
103 : }
104 87 : }
105 :
106 : /**
107 : * Parses the phase center strings in the items of a dict/record.
108 : * Populates the phaseCenterSpec_p map.
109 : *
110 : * @param config TVI configuration object
111 : */
112 5 : void PhaseShiftingTVI::parsePhasecenterDict(const Record &config)
113 : {
114 5 : const auto &inPhasecenter = config.asRecord("phasecenter");
115 5 : if (inPhasecenter.empty()) {
116 0 : throw AipsError("The dictionary 'phasecenter' is empty.");
117 : }
118 :
119 5 : const auto maxMSField = getMaxMSFieldID();
120 5 : std::set<unsigned int> fieldsSeen;
121 : // Go through items in the input dict/record
122 13 : for (unsigned int rid=0; rid < inPhasecenter.nfields(); ++rid) {
123 18 : const std::string fieldStr = inPhasecenter.name(RecordFieldId(rid));
124 9 : const auto fid = std::stoi(fieldStr);
125 9 : if (fid < 0 || static_cast<unsigned int>(fid) > maxMSField) {
126 3 : throw AipsError("Wrong field ID given: " + std::to_string(fid) +
127 2 : ". This MeasurementSet has field IDs between 0 and " +
128 2 : std::to_string(maxMSField));
129 : }
130 8 : if (not fieldsSeen.insert(fid).second) {
131 0 : throw AipsError("Field " + std::to_string(fid) + " is given multiple times");
132 : }
133 :
134 8 : std::string center;
135 : try {
136 8 : center = inPhasecenter.asString(RecordFieldId(rid));
137 0 : } catch (const AipsError &exc) {
138 0 : throw AipsError("For field " + std::to_string(fid) + ", cannot interpret "
139 0 : "phasecenter value as a string: " +
140 0 : std::string(exc.getMesg()));
141 0 : }
142 8 : const auto strFieldID = std::to_string(fid);
143 8 : const auto &phaseCenterDir = checkPhaseCenterStr(center, strFieldID);
144 8 : phaseCenterSpec_p.insert({fid, phaseCenterDir});
145 9 : }
146 5 : }
147 :
148 : /**
149 : * Checks that a phase center string can be correctly converted to a
150 : * MDirection object, and checks that the ref frame is supported.
151 : *
152 : * @param phasecenter center string as given in the 'phasecenter' param
153 : * @param fieldInfo more info to print about the field
154 : *
155 : * @returns direction as an MDirection
156 : */
157 90 : MDirection PhaseShiftingTVI::checkPhaseCenterStr(const String &phasecenter,
158 : const string &fieldInfo)
159 : {
160 : // casaMDirection requires a variant
161 90 : casac::variant phaseCenterVar(phasecenter);
162 90 : casacore::MDirection phaseCenterDir;
163 90 : if(!casaMDirection(phaseCenterVar, phaseCenterDir)) {
164 0 : throw AipsError("Cannot interpret phase center string as a direction object: "
165 0 : + phasecenter);
166 : } else {
167 90 : const auto myFrame = phaseCenterDir.getRefString();
168 : MDirection::Types mdtype;
169 90 : MDirection::getType(mdtype, myFrame);
170 90 : ThrowIf(
171 : mdtype == MDirection::HADEC || mdtype == MDirection::AZEL
172 : || mdtype == MDirection::AZELSW || mdtype == MDirection::AZELNE
173 : || mdtype == MDirection::AZELGEO || mdtype == MDirection::AZELSWGEO
174 : || mdtype == MDirection::MECLIPTIC || mdtype == MDirection::TECLIPTIC
175 : || mdtype == MDirection::TOPO,
176 : myFrame + " is a time dependent reference frame and so is not supported"
177 : );
178 90 : ThrowIf(
179 : mdtype == MDirection::MERCURY || mdtype == MDirection::VENUS
180 : || mdtype == MDirection::MARS || mdtype == MDirection::JUPITER
181 : || mdtype == MDirection::SATURN || mdtype == MDirection::URANUS
182 : || mdtype == MDirection::NEPTUNE || mdtype == MDirection::PLUTO
183 : || mdtype == MDirection::SUN || mdtype == MDirection::MOON
184 : || mdtype == MDirection::COMET,
185 : myFrame + " denotes an ephemeris object and so is not supported"
186 : );
187 :
188 180 : auto msg = "Phase center '" + std::string(phasecenter) +
189 90 : "' successfully parsed";
190 90 : if (not fieldInfo.empty()) {
191 8 : msg += " for field " + fieldInfo;
192 : }
193 180 : logger_p << LogIO::NORMAL << LogOrigin("PhaseShiftingTVI", __FUNCTION__)
194 180 : << msg << LogIO::POST;
195 90 : }
196 :
197 180 : return phaseCenterDir;
198 90 : }
199 :
200 : // -----------------------------------------------------------------------
201 : //
202 : // -----------------------------------------------------------------------
203 87 : void PhaseShiftingTVI::parseConfiguration(const Record &configuration)
204 : {
205 87 : int exists = -1;
206 :
207 87 : exists = -1;
208 87 : exists = configuration.fieldNumber ("XpcOffset");
209 87 : if (exists >= 0)
210 : {
211 0 : configuration.get (exists, dx_p);
212 : }
213 :
214 87 : exists = -1;
215 87 : exists = configuration.fieldNumber ("YpcOffset");
216 87 : if (exists >= 0)
217 : {
218 0 : configuration.get (exists, dy_p);
219 : }
220 :
221 87 : if (abs(dx_p) > 0 or abs(dy_p) > 0)
222 : {
223 0 : logger_p << LogIO::NORMAL << LogOrigin("PhaseShiftingTVI", __FUNCTION__)
224 0 : << "Phase shift is dx="<< dx_p << " dy=" << dy_p << LogIO::POST;
225 : }
226 :
227 : // CAS-12706 Add support for shifting across large offset/angles
228 87 : parsePhasecenter(configuration);
229 86 : wideFieldMode_p = true;
230 86 : }
231 :
232 : /**
233 : * Finds the phase center (converted as MDirection object) for the current
234 : * field, converting the ref frame if needed to match the current Vis Buffer
235 : * ref frame.
236 : *
237 : * @returns whether any shift should be applied (if an output phasecenter is
238 : * defined for this field + phase center to shift to, for current
239 : * field
240 : */
241 : std::pair<bool, MDirection>
242 5181 : PhaseShiftingTVI::findConvertedPhaseCenter() const
243 : {
244 5181 : const auto *vb = getVii()->getVisBuffer();
245 :
246 5181 : auto centerIt = phaseCenterSpec_p.find(-1);
247 5181 : if (centerIt == phaseCenterSpec_p.end()) {
248 2216 : auto fieldID = vb->fieldId()[0];
249 2216 : centerIt = phaseCenterSpec_p.find(fieldID);
250 2216 : if (centerIt == phaseCenterSpec_p.end()) {
251 1576 : return {false, MDirection()};
252 : }
253 : }
254 4393 : auto convertedPhaseCenter = centerIt->second;
255 :
256 4393 : if (convertedPhaseCenter.getRefString() != vb->phaseCenter().getRefString()) {
257 : MDirection::Types mdtype;
258 1915 : MDirection::getType(mdtype, vb->phaseCenter().getRefString());
259 1915 : convertedPhaseCenter = MDirection::Convert(convertedPhaseCenter, mdtype)();
260 : }
261 :
262 4393 : return {true, convertedPhaseCenter};
263 4393 : }
264 :
265 : // -----------------------------------------------------------------------
266 : //
267 : // -----------------------------------------------------------------------
268 86 : void PhaseShiftingTVI::initialize()
269 : {
270 : // Populate nchan input-output maps
271 : Int spw;
272 86 : uInt spw_idx = 0;
273 86 : map<Int,vector<Int> >::iterator iter;
274 179 : for (iter=spwInpChanIdxMap_p.begin();iter!=spwInpChanIdxMap_p.end();iter++)
275 : {
276 93 : spw = iter->first;
277 93 : spwOutChanNumMap_p[spw] = spwInpChanIdxMap_p[spw].size();
278 :
279 93 : spw_idx++;
280 : }
281 :
282 : // CAS-12706 Add support for shifting across large offset/angles
283 : // Access observatory position and observation start (reference) time.
284 86 : if (wideFieldMode_p)
285 : {
286 86 : const auto selectedInputMsCols = std::make_unique<MSColumns>(getVii()->ms());
287 86 : observatoryPosition_p = selectedInputMsCols->antenna().positionMeas()(0);
288 86 : referenceTime_p = selectedInputMsCols->timeMeas()(0);
289 86 : referenceTimeUnits_p = selectedInputMsCols->timeQuant()(0).getUnit();
290 86 : }
291 86 : }
292 :
293 : // -----------------------------------------------------------------------
294 : // CAS-12706 Recalculate UVW and Phases according to wide-field phase shifting
295 : // algorithm. This method should do the same as ComponentFTMachine::rotateUVW
296 : // -----------------------------------------------------------------------
297 5181 : void PhaseShiftingTVI::shiftUVWPhases()
298 : {
299 : // Get input VisBuffer
300 5181 : VisBuffer2 *vb = getVii()->getVisBuffer();
301 :
302 5181 : bool doShift = false;
303 5181 : MDirection convertedPhaseCenter;
304 5181 : std::tie(doShift, convertedPhaseCenter) = findConvertedPhaseCenter();
305 5181 : if (not doShift) {
306 788 : phaseShift_p.resize(0, false);
307 788 : newUVW_p = vb->uvw();
308 788 : return;
309 : }
310 :
311 : // Initialize epoch corresponding to current buffer
312 : // with time reference to the first row in the MS
313 8786 : MEpoch epoch(Quantity(vb->time()(0),referenceTimeUnits_p),referenceTime_p.getRef());
314 4393 : MeasFrame refFrame(epoch,observatoryPosition_p);
315 4393 : UVWMachine uvwMachine(convertedPhaseCenter, vb->phaseCenter(), refFrame,false,false);
316 : // Initialize phase array and uvw matrix
317 4393 : phaseShift_p.resize(vb->nRows(),false);
318 4393 : newUVW_p.resize(vb->uvw().shape(),false);
319 :
320 : // Obtain phase shift and new uvw coordinates
321 4393 : Vector<Double> dummy(3,0.0);
322 4393 : double phase2radPerHz = -2.0 * C::pi / C::c;
323 770351 : for (rownr_t row=0; row<vb->nRows(); row++)
324 : {
325 : // Copy current uvw coordinates so that they are not modified
326 : // Note: Columns in uvw correspond to rows in the main table/VisBuffer!
327 765958 : dummy = vb->uvw().column(row);
328 :
329 : // Have to change (u,v,w) to (-u,-v,w) because is the convention used by uvwMachine
330 765958 : dummy(0) = -1*dummy(0);
331 765958 : dummy(1) = -1*dummy(1);
332 :
333 : // Transform uvw coordinates and obtain corresponding phase shift
334 765958 : uvwMachine.convertUVW(phaseShift_p(row), dummy);
335 :
336 : // Store new uvw coordinates
337 : // Have to change back (-u,-v,w) to (u,v,w) because is the convention used by the MS
338 765958 : dummy(0) = -1*dummy(0);
339 765958 : dummy(1) = -1*dummy(1);
340 765958 : newUVW_p.column(row) = dummy;
341 : // Convert phase shift to radian/Hz
342 765958 : phaseShift_p(row) = phase2radPerHz*phaseShift_p(row);
343 : }
344 5181 : }
345 :
346 : // -----------------------------------------------------------------------
347 : //
348 : // -----------------------------------------------------------------------
349 540 : void PhaseShiftingTVI::origin()
350 : {
351 : // Drive underlying ViImplementation2
352 540 : getVii()->origin();
353 :
354 : // CAS-12706 Add support for shifting across large offset/angles
355 540 : if (wideFieldMode_p) {
356 540 : shiftUVWPhases();
357 : }
358 :
359 : // Define the shapes in the VB2, patch provided by cgarcia in CAS-12706
360 540 : configureShapes();
361 :
362 : // Synchronize own VisBuffer
363 540 : configureNewSubchunk();
364 540 : }
365 :
366 : // -----------------------------------------------------------------------
367 : //
368 : // -----------------------------------------------------------------------
369 4641 : void PhaseShiftingTVI::next()
370 : {
371 : // Drive underlying ViImplementation2
372 4641 : getVii()->next();
373 :
374 : // CAS-12706 Add support for shifting across large offset/angles
375 4641 : if (wideFieldMode_p) {
376 4641 : shiftUVWPhases();
377 : }
378 :
379 : // Define the shapes in the VB2, patch provided by cgarcia in CAS-12706
380 4641 : configureShapes();
381 :
382 : // Synchronize own VisBuffer
383 4641 : configureNewSubchunk();
384 4641 : }
385 :
386 :
387 : // -----------------------------------------------------------------------
388 : //
389 : // -----------------------------------------------------------------------
390 4641 : void PhaseShiftingTVI::visibilityObserved (Cube<Complex> & vis) const
391 : {
392 : // Get input VisBuffer
393 4641 : VisBuffer2 *vb = getVii()->getVisBuffer();
394 4641 : Matrix<Double> uvw = vb->uvw();
395 4641 : Vector<Double> frequencies = vb->getFrequencies(0);
396 :
397 : // Reshape output data before passing it to the DataCubeHolder
398 4641 : vis.resize(getVisBuffer()->getShape(),false);
399 :
400 : // Gather input data
401 4641 : DataCubeMap inputData;
402 4641 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCube());
403 4641 : inputData.add(MS::DATA,inputVisCubeHolder);
404 :
405 : // Gather output data
406 4641 : DataCubeMap outputData;
407 4641 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
408 4641 : outputData.add(MS::DATA,outputVisCubeHolder);
409 :
410 4641 : if (wideFieldMode_p)
411 : {
412 : // Configure Transformation Engine
413 4641 : WideFieldPhaseShiftingTransformEngine<Complex> transformer(phaseShift_p,&uvw,&frequencies,&inputData,&outputData);
414 :
415 : // Transform data
416 4641 : transformFreqAxis2(vb->getShape(),transformer);
417 4641 : }
418 : else
419 : {
420 : // Configure Transformation Engine
421 0 : PhaseShiftingTransformEngine<Complex> transformer(dx_p,dy_p,&uvw,&frequencies,&inputData,&outputData);
422 :
423 : // Transform data
424 0 : transformFreqAxis2(vb->getShape(),transformer);
425 : }
426 4641 : }
427 :
428 : // -----------------------------------------------------------------------
429 : //
430 : // -----------------------------------------------------------------------
431 1757 : void PhaseShiftingTVI::visibilityCorrected (Cube<Complex> & vis) const
432 : {
433 : // Get input VisBuffer
434 1757 : VisBuffer2 *vb = getVii()->getVisBuffer();
435 1757 : Matrix<Double> uvw = vb->uvw();
436 1757 : Vector<Double> frequencies = vb->getFrequencies(0);
437 :
438 : // Reshape output data before passing it to the DataCubeHolder
439 1757 : vis.resize(getVisBuffer()->getShape(),false);
440 :
441 : // Gather input data
442 1757 : DataCubeMap inputData;
443 1757 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCubeCorrected());
444 1757 : inputData.add(MS::DATA,inputVisCubeHolder);
445 :
446 : // Gather output data
447 1757 : DataCubeMap outputData;
448 1757 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
449 1757 : outputData.add(MS::DATA,outputVisCubeHolder);
450 :
451 1757 : if (wideFieldMode_p)
452 : {
453 : // Configure Transformation Engine
454 1757 : WideFieldPhaseShiftingTransformEngine<Complex> transformer(phaseShift_p,&uvw,&frequencies,&inputData,&outputData);
455 :
456 : // Transform data
457 1757 : transformFreqAxis2(vb->getShape(),transformer);
458 1757 : }
459 : else
460 : {
461 : // Configure Transformation Engine
462 0 : PhaseShiftingTransformEngine<Complex> transformer(dx_p,dy_p,&uvw,&frequencies,&inputData,&outputData);
463 :
464 : // Transform data
465 0 : transformFreqAxis2(vb->getShape(),transformer);
466 : }
467 1757 : }
468 :
469 : // -----------------------------------------------------------------------
470 : //
471 : // -----------------------------------------------------------------------
472 1757 : void PhaseShiftingTVI::visibilityModel (Cube<Complex> & vis) const
473 : {
474 : // Get input VisBuffer
475 1757 : VisBuffer2 *vb = getVii()->getVisBuffer();
476 1757 : Matrix<Double> uvw = vb->uvw();
477 1757 : Vector<Double> frequencies = vb->getFrequencies(0);
478 :
479 : // Reshape output data before passing it to the DataCubeHolder
480 1757 : vis.resize(getVisBuffer()->getShape(),false);
481 :
482 : // Gather input data
483 1757 : DataCubeMap inputData;
484 1757 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCubeModel());
485 1757 : inputData.add(MS::DATA,inputVisCubeHolder);
486 :
487 : // Gather output data
488 1757 : DataCubeMap outputData;
489 1757 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
490 1757 : outputData.add(MS::DATA,outputVisCubeHolder);
491 :
492 1757 : if (wideFieldMode_p)
493 : {
494 : // Configure Transformation Engine
495 1757 : WideFieldPhaseShiftingTransformEngine<Complex> transformer(phaseShift_p,&uvw,&frequencies,&inputData,&outputData);
496 :
497 : // Transform data
498 1757 : transformFreqAxis2(vb->getShape(),transformer);
499 1757 : }
500 : else
501 : {
502 : // Configure Transformation Engine
503 0 : PhaseShiftingTransformEngine<Complex> transformer(dx_p,dy_p,&uvw,&frequencies,&inputData,&outputData);
504 :
505 : // Transform data
506 0 : transformFreqAxis2(vb->getShape(),transformer);
507 : }
508 1757 : }
509 :
510 : // -----------------------------------------------------------------------
511 : //
512 : // -----------------------------------------------------------------------
513 4641 : void PhaseShiftingTVI::uvw (casacore::Matrix<double> & uvw) const
514 : {
515 4641 : if (wideFieldMode_p) {
516 4641 : uvw.resize(newUVW_p.shape(),false);
517 4641 : uvw = newUVW_p;
518 : }
519 : else {
520 0 : getVii()->uvw (uvw);
521 : }
522 4641 : }
523 :
524 : //////////////////////////////////////////////////////////////////////////
525 : // PhaseShiftingTVIFactory class
526 : //////////////////////////////////////////////////////////////////////////
527 :
528 : // -----------------------------------------------------------------------
529 : //
530 : // -----------------------------------------------------------------------
531 0 : PhaseShiftingTVIFactory::PhaseShiftingTVIFactory ( Record &configuration,
532 0 : ViImplementation2 *inputVii)
533 : {
534 0 : inputVii_p = inputVii;
535 0 : configuration_p = configuration;
536 0 : }
537 :
538 : // -----------------------------------------------------------------------
539 : //
540 : // -----------------------------------------------------------------------
541 0 : vi::ViImplementation2 * PhaseShiftingTVIFactory::createVi(VisibilityIterator2 *) const
542 : {
543 0 : return new PhaseShiftingTVI(inputVii_p,configuration_p);
544 : }
545 :
546 : // -----------------------------------------------------------------------
547 : //
548 : // -----------------------------------------------------------------------
549 0 : vi::ViImplementation2 * PhaseShiftingTVIFactory::createVi() const
550 : {
551 0 : return new PhaseShiftingTVI(inputVii_p,configuration_p);
552 : }
553 :
554 : //////////////////////////////////////////////////////////////////////////
555 : // PhaseShiftingTVILayerFactory class
556 : //////////////////////////////////////////////////////////////////////////
557 :
558 87 : PhaseShiftingTVILayerFactory::PhaseShiftingTVILayerFactory(Record &configuration) :
559 : ViiLayerFactory(),
560 87 : configuration_p(configuration)
561 87 : {}
562 :
563 : ViImplementation2*
564 87 : PhaseShiftingTVILayerFactory::createInstance(ViImplementation2* vii0) const
565 : {
566 : // Make the PhaseShiftingTVi2, using supplied ViImplementation2, and return it
567 87 : ViImplementation2 *vii = new PhaseShiftingTVI(vii0,configuration_p);
568 86 : return vii;
569 : }
570 :
571 : //////////////////////////////////////////////////////////////////////////
572 : // PhaseShiftingTransformEngine class
573 : //////////////////////////////////////////////////////////////////////////
574 :
575 : // -----------------------------------------------------------------------
576 : //
577 : // -----------------------------------------------------------------------
578 0 : template<class T> PhaseShiftingTransformEngine<T>::PhaseShiftingTransformEngine(
579 : Double dx, Double dy,
580 : Matrix<Double> *uvw,
581 : Vector<Double> *frequencies,
582 : DataCubeMap *inputData,
583 : DataCubeMap *outputData):
584 0 : FreqAxisTransformEngine2<T>(inputData,outputData)
585 : {
586 0 : uvw_p = uvw;
587 0 : frequencies_p = frequencies;
588 :
589 : // Offsets in radians (input is arcsec)
590 0 : dx_p = dx*(C::pi / 180.0 / 3600.0);
591 0 : dy_p = dy*(C::pi / 180.0 / 3600.0);
592 0 : }
593 :
594 : // -----------------------------------------------------------------------
595 : //
596 : // -----------------------------------------------------------------------
597 0 : template<class T> void PhaseShiftingTransformEngine<T>::transform( )
598 : {
599 0 : transformCore(inputData_p,outputData_p);
600 0 : }
601 :
602 : // -----------------------------------------------------------------------
603 : //
604 : // -----------------------------------------------------------------------
605 0 : template<class T> void PhaseShiftingTransformEngine<T>::transformCore( DataCubeMap *inputData,
606 : DataCubeMap *outputData)
607 : {
608 : // Get input/output data
609 0 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
610 0 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
611 :
612 : // Extra path as fraction of U and V in m
613 0 : Double phase = dx_p*(*uvw_p)(0,rowIndex_p) + dy_p*(*uvw_p)(1,rowIndex_p);
614 :
615 : // In radian/Hz
616 0 : phase *= -2.0 * C::pi / C::c;
617 :
618 : // Main loop
619 : Double phase_i;
620 0 : Complex factor;
621 0 : for (uInt chan_i=0;chan_i<inputVector.size();chan_i++)
622 : {
623 0 : phase_i = phase * (*frequencies_p)(chan_i);
624 0 : factor = Complex(cos(phase_i), sin(phase_i));
625 0 : outputVector(chan_i) = factor*inputVector(chan_i);
626 : }
627 0 : }
628 :
629 : //////////////////////////////////////////////////////////////////////////
630 : // WideFieldPhaseShiftingTransformEngine class
631 : //////////////////////////////////////////////////////////////////////////
632 :
633 : // -----------------------------------------------------------------------
634 : //
635 : // -----------------------------------------------------------------------
636 8155 : template<class T> WideFieldPhaseShiftingTransformEngine<T>::WideFieldPhaseShiftingTransformEngine(
637 : const Vector<Double> &phaseShift,
638 : Matrix<Double> *uvw,
639 : Vector<Double> *frequencies,
640 : DataCubeMap *inputData,
641 : DataCubeMap *outputData):
642 8155 : FreqAxisTransformEngine2<T>(inputData,outputData)
643 : {
644 8155 : uvw_p = uvw;
645 8155 : frequencies_p = frequencies;
646 8155 : phaseShift_p = phaseShift;
647 8155 : }
648 :
649 : // -----------------------------------------------------------------------
650 : //
651 : // -----------------------------------------------------------------------
652 3944008 : template<class T> void WideFieldPhaseShiftingTransformEngine<T>::transform( )
653 : {
654 3944008 : transformCore(inputData_p,outputData_p);
655 3944008 : }
656 :
657 : // -----------------------------------------------------------------------
658 : //
659 : // -----------------------------------------------------------------------
660 3944008 : template<class T> void WideFieldPhaseShiftingTransformEngine<T>::transformCore( DataCubeMap *inputData,
661 : DataCubeMap *outputData)
662 : {
663 : // Get input/output data
664 3944008 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
665 3944008 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
666 :
667 3944008 : if (phaseShift_p.shape() == 0) {
668 : // no shift, bypass data as 'passthrough' field
669 38304 : outputVector = inputVector;
670 38304 : return;
671 : }
672 :
673 : // Main loop
674 : Double phase_i;
675 3905704 : Complex factor;
676 11835540 : for (uInt chan_i=0;chan_i<inputVector.size();chan_i++)
677 : {
678 7929836 : phase_i = phaseShift_p(rowIndex_p) * (*frequencies_p)(chan_i);
679 7929836 : factor = Complex(cos(phase_i), sin(phase_i));
680 7929836 : outputVector(chan_i) = factor*inputVector(chan_i);
681 : }
682 : }
683 :
684 :
685 : } //# NAMESPACE VI - END
686 :
687 : } //# NAMESPACE CASA - END
688 :
689 :
|
