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 88 : PhaseShiftingTVI::PhaseShiftingTVI( ViImplementation2 * inputVii,
38 88 : const Record &configuration):
39 88 : FreqAxisTVI (inputVii)
40 : {
41 88 : dx_p = 0;
42 88 : dy_p = 0;
43 :
44 : // CAS-12706 Zero initialization for wide-field phase shifting algorithm
45 88 : 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 88 : parseConfiguration(configuration);
51 :
52 87 : initialize();
53 :
54 87 : 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 6 : rownr_t PhaseShiftingTVI::getMaxMSFieldID() const
62 : {
63 6 : const auto &fieldsTable = getVii()->fieldSubtablecols();
64 6 : 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 88 : void PhaseShiftingTVI::parsePhasecenter(const Record &config)
75 : {
76 88 : auto exists = config.fieldNumber ("phasecenter");
77 88 : if (exists < 0) {
78 0 : return;
79 : }
80 :
81 : // phasecenter can be given as a string or as a dict (per-field centers)
82 88 : auto isStr = false;
83 88 : casacore::String phaseCenterStr;
84 : try {
85 94 : 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 6 : } 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 12 : logger_p << LogIO::NORMAL << LogOrigin("PhaseShiftingTVI", __FUNCTION__)
94 : << "Interpreting phasecenter as a dictionary of field->center."
95 12 : << "Using following phase centers:" << LogIO::POST;
96 6 : }
97 :
98 88 : if (isStr) {
99 82 : const auto &phaseCenterDir = checkPhaseCenterStr(phaseCenterStr);
100 82 : phaseCenterSpec_p.insert({-1, phaseCenterDir});
101 82 : } else {
102 6 : parsePhasecenterDict(config);
103 : }
104 88 : }
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 6 : void PhaseShiftingTVI::parsePhasecenterDict(const Record &config)
113 : {
114 6 : const auto &inPhasecenter = config.asRecord("phasecenter");
115 6 : if (inPhasecenter.empty()) {
116 0 : throw AipsError("The dictionary 'phasecenter' is empty.");
117 : }
118 :
119 6 : const auto maxMSField = getMaxMSFieldID();
120 6 : std::set<unsigned int> fieldsSeen;
121 : // Go through items in the input dict/record
122 15 : for (unsigned int rid=0; rid < inPhasecenter.nfields(); ++rid) {
123 20 : const std::string fieldStr = inPhasecenter.name(RecordFieldId(rid));
124 10 : const auto fid = std::stoi(fieldStr);
125 10 : 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 9 : if (not fieldsSeen.insert(fid).second) {
131 0 : throw AipsError("Field " + std::to_string(fid) + " is given multiple times");
132 : }
133 :
134 9 : std::string center;
135 : try {
136 9 : 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 9 : const auto strFieldID = std::to_string(fid);
143 9 : const auto &phaseCenterDir = checkPhaseCenterStr(center, strFieldID);
144 9 : phaseCenterSpec_p.insert({fid, phaseCenterDir});
145 10 : }
146 6 : }
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 91 : MDirection PhaseShiftingTVI::checkPhaseCenterStr(const String &phasecenter,
158 : const string &fieldInfo)
159 : {
160 : // casaMDirection requires a variant
161 91 : casac::variant phaseCenterVar(phasecenter);
162 91 : casacore::MDirection phaseCenterDir;
163 91 : if(!casaMDirection(phaseCenterVar, phaseCenterDir)) {
164 0 : throw AipsError("Cannot interpret phase center string as a direction object: "
165 0 : + phasecenter);
166 : } else {
167 91 : const auto myFrame = phaseCenterDir.getRefString();
168 : MDirection::Types mdtype;
169 91 : MDirection::getType(mdtype, myFrame);
170 91 : 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 91 : 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 182 : auto msg = "Phase center '" + std::string(phasecenter) +
189 91 : "' successfully parsed";
190 91 : if (not fieldInfo.empty()) {
191 9 : msg += " for field " + fieldInfo;
192 : }
193 182 : logger_p << LogIO::NORMAL << LogOrigin("PhaseShiftingTVI", __FUNCTION__)
194 182 : << msg << LogIO::POST;
195 91 : }
196 :
197 182 : return phaseCenterDir;
198 91 : }
199 :
200 : // -----------------------------------------------------------------------
201 : //
202 : // -----------------------------------------------------------------------
203 88 : void PhaseShiftingTVI::parseConfiguration(const Record &configuration)
204 : {
205 88 : int exists = -1;
206 :
207 88 : exists = -1;
208 88 : exists = configuration.fieldNumber ("XpcOffset");
209 88 : if (exists >= 0)
210 : {
211 0 : configuration.get (exists, dx_p);
212 : }
213 :
214 88 : exists = -1;
215 88 : exists = configuration.fieldNumber ("YpcOffset");
216 88 : if (exists >= 0)
217 : {
218 0 : configuration.get (exists, dy_p);
219 : }
220 :
221 88 : 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 88 : parsePhasecenter(configuration);
229 87 : wideFieldMode_p = true;
230 87 : }
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 5185 : PhaseShiftingTVI::findConvertedPhaseCenter() const
243 : {
244 5185 : const auto *vb = getVii()->getVisBuffer();
245 :
246 5185 : auto centerIt = phaseCenterSpec_p.find(-1);
247 5185 : if (centerIt == phaseCenterSpec_p.end()) {
248 2220 : auto fieldID = vb->fieldId()[0];
249 2220 : centerIt = phaseCenterSpec_p.find(fieldID);
250 2220 : if (centerIt == phaseCenterSpec_p.end()) {
251 1580 : return {false, MDirection()};
252 : }
253 : }
254 4395 : auto convertedPhaseCenter = centerIt->second;
255 :
256 4395 : 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 4395 : return {true, convertedPhaseCenter};
263 4395 : }
264 :
265 : // -----------------------------------------------------------------------
266 : //
267 : // -----------------------------------------------------------------------
268 87 : void PhaseShiftingTVI::initialize()
269 : {
270 : // Populate nchan input-output maps
271 : Int spw;
272 87 : uInt spw_idx = 0;
273 87 : map<Int,vector<Int> >::iterator iter;
274 181 : for (iter=spwInpChanIdxMap_p.begin();iter!=spwInpChanIdxMap_p.end();iter++)
275 : {
276 94 : spw = iter->first;
277 94 : spwOutChanNumMap_p[spw] = spwInpChanIdxMap_p[spw].size();
278 :
279 94 : 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 87 : if (wideFieldMode_p)
285 : {
286 87 : const auto selectedInputMsCols = std::make_unique<MSColumns>(getVii()->ms());
287 87 : observatoryPosition_p = selectedInputMsCols->antenna().positionMeas()(0);
288 87 : referenceTime_p = selectedInputMsCols->timeMeas()(0);
289 87 : referenceTimeUnits_p = selectedInputMsCols->timeQuant()(0).getUnit();
290 87 : }
291 87 : }
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 5185 : void PhaseShiftingTVI::shiftUVWPhases()
298 : {
299 : // Get input VisBuffer
300 5185 : VisBuffer2 *vb = getVii()->getVisBuffer();
301 :
302 5185 : bool doShift = false;
303 5185 : MDirection convertedPhaseCenter;
304 5185 : std::tie(doShift, convertedPhaseCenter) = findConvertedPhaseCenter();
305 5185 : if (not doShift) {
306 790 : phaseShift_p.resize(0, false);
307 : // reinit to avoid shape mismatch, CAS-14633
308 790 : auto temp = vb->uvw();
309 790 : newUVW_p.resize(temp.shape(), false);
310 790 : newUVW_p = temp;
311 790 : return;
312 790 : }
313 :
314 : // Initialize epoch corresponding to current buffer
315 : // with time reference to the first row in the MS
316 8790 : MEpoch epoch(Quantity(vb->time()(0),referenceTimeUnits_p),referenceTime_p.getRef());
317 4395 : MeasFrame refFrame(epoch,observatoryPosition_p);
318 4395 : UVWMachine uvwMachine(convertedPhaseCenter, vb->phaseCenter(), refFrame,false,false);
319 : // Initialize phase array and uvw matrix
320 4395 : phaseShift_p.resize(vb->nRows(),false);
321 4395 : newUVW_p.resize(vb->uvw().shape(),false);
322 :
323 : // Obtain phase shift and new uvw coordinates
324 4395 : Vector<Double> dummy(3,0.0);
325 4395 : double phase2radPerHz = -2.0 * M_PI / C::c;
326 770463 : for (rownr_t row=0; row<vb->nRows(); row++)
327 : {
328 : // Copy current uvw coordinates so that they are not modified
329 : // Note: Columns in uvw correspond to rows in the main table/VisBuffer!
330 766068 : dummy = vb->uvw().column(row);
331 :
332 : // Have to change (u,v,w) to (-u,-v,w) because is the convention used by uvwMachine
333 766068 : dummy(0) = -1*dummy(0);
334 766068 : dummy(1) = -1*dummy(1);
335 :
336 : // Transform uvw coordinates and obtain corresponding phase shift
337 766068 : uvwMachine.convertUVW(phaseShift_p(row), dummy);
338 :
339 : // Store new uvw coordinates
340 : // Have to change back (-u,-v,w) to (u,v,w) because is the convention used by the MS
341 766068 : dummy(0) = -1*dummy(0);
342 766068 : dummy(1) = -1*dummy(1);
343 766068 : newUVW_p.column(row) = dummy;
344 : // Convert phase shift to radian/Hz
345 766068 : phaseShift_p(row) = phase2radPerHz*phaseShift_p(row);
346 : }
347 5185 : }
348 :
349 : // -----------------------------------------------------------------------
350 : //
351 : // -----------------------------------------------------------------------
352 542 : void PhaseShiftingTVI::origin()
353 : {
354 : // Drive underlying ViImplementation2
355 542 : getVii()->origin();
356 :
357 : // CAS-12706 Add support for shifting across large offset/angles
358 542 : if (wideFieldMode_p) {
359 542 : shiftUVWPhases();
360 : }
361 :
362 : // Define the shapes in the VB2, patch provided by cgarcia in CAS-12706
363 542 : configureShapes();
364 :
365 : // Synchronize own VisBuffer
366 542 : configureNewSubchunk();
367 542 : }
368 :
369 : // -----------------------------------------------------------------------
370 : //
371 : // -----------------------------------------------------------------------
372 4643 : void PhaseShiftingTVI::next()
373 : {
374 : // Drive underlying ViImplementation2
375 4643 : getVii()->next();
376 :
377 : // CAS-12706 Add support for shifting across large offset/angles
378 4643 : if (wideFieldMode_p) {
379 4643 : shiftUVWPhases();
380 : }
381 :
382 : // Define the shapes in the VB2, patch provided by cgarcia in CAS-12706
383 4643 : configureShapes();
384 :
385 : // Synchronize own VisBuffer
386 4643 : configureNewSubchunk();
387 4643 : }
388 :
389 :
390 : // -----------------------------------------------------------------------
391 : //
392 : // -----------------------------------------------------------------------
393 4643 : void PhaseShiftingTVI::visibilityObserved (Cube<Complex> & vis) const
394 : {
395 : // Get input VisBuffer
396 4643 : VisBuffer2 *vb = getVii()->getVisBuffer();
397 4643 : Matrix<Double> uvw = vb->uvw();
398 4643 : Vector<Double> frequencies = vb->getFrequencies(0);
399 :
400 : // Reshape output data before passing it to the DataCubeHolder
401 4643 : vis.resize(getVisBuffer()->getShape(),false);
402 :
403 : // Gather input data
404 4643 : DataCubeMap inputData;
405 4643 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCube());
406 4643 : inputData.add(MS::DATA,inputVisCubeHolder);
407 :
408 : // Gather output data
409 4643 : DataCubeMap outputData;
410 4643 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
411 4643 : outputData.add(MS::DATA,outputVisCubeHolder);
412 :
413 4643 : if (wideFieldMode_p)
414 : {
415 : // Configure Transformation Engine
416 4643 : WideFieldPhaseShiftingTransformEngine<Complex> transformer(phaseShift_p,&uvw,&frequencies,&inputData,&outputData);
417 :
418 : // Transform data
419 4643 : transformFreqAxis2(vb->getShape(),transformer);
420 4643 : }
421 : else
422 : {
423 : // Configure Transformation Engine
424 0 : PhaseShiftingTransformEngine<Complex> transformer(dx_p,dy_p,&uvw,&frequencies,&inputData,&outputData);
425 :
426 : // Transform data
427 0 : transformFreqAxis2(vb->getShape(),transformer);
428 : }
429 4643 : }
430 :
431 : // -----------------------------------------------------------------------
432 : //
433 : // -----------------------------------------------------------------------
434 1757 : void PhaseShiftingTVI::visibilityCorrected (Cube<Complex> & vis) const
435 : {
436 : // Get input VisBuffer
437 1757 : VisBuffer2 *vb = getVii()->getVisBuffer();
438 1757 : Matrix<Double> uvw = vb->uvw();
439 1757 : Vector<Double> frequencies = vb->getFrequencies(0);
440 :
441 : // Reshape output data before passing it to the DataCubeHolder
442 1757 : vis.resize(getVisBuffer()->getShape(),false);
443 :
444 : // Gather input data
445 1757 : DataCubeMap inputData;
446 1757 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCubeCorrected());
447 1757 : inputData.add(MS::DATA,inputVisCubeHolder);
448 :
449 : // Gather output data
450 1757 : DataCubeMap outputData;
451 1757 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
452 1757 : outputData.add(MS::DATA,outputVisCubeHolder);
453 :
454 1757 : if (wideFieldMode_p)
455 : {
456 : // Configure Transformation Engine
457 1757 : WideFieldPhaseShiftingTransformEngine<Complex> transformer(phaseShift_p,&uvw,&frequencies,&inputData,&outputData);
458 :
459 : // Transform data
460 1757 : transformFreqAxis2(vb->getShape(),transformer);
461 1757 : }
462 : else
463 : {
464 : // Configure Transformation Engine
465 0 : PhaseShiftingTransformEngine<Complex> transformer(dx_p,dy_p,&uvw,&frequencies,&inputData,&outputData);
466 :
467 : // Transform data
468 0 : transformFreqAxis2(vb->getShape(),transformer);
469 : }
470 1757 : }
471 :
472 : // -----------------------------------------------------------------------
473 : //
474 : // -----------------------------------------------------------------------
475 1757 : void PhaseShiftingTVI::visibilityModel (Cube<Complex> & vis) const
476 : {
477 : // Get input VisBuffer
478 1757 : VisBuffer2 *vb = getVii()->getVisBuffer();
479 1757 : Matrix<Double> uvw = vb->uvw();
480 1757 : Vector<Double> frequencies = vb->getFrequencies(0);
481 :
482 : // Reshape output data before passing it to the DataCubeHolder
483 1757 : vis.resize(getVisBuffer()->getShape(),false);
484 :
485 : // Gather input data
486 1757 : DataCubeMap inputData;
487 1757 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCubeModel());
488 1757 : inputData.add(MS::DATA,inputVisCubeHolder);
489 :
490 : // Gather output data
491 1757 : DataCubeMap outputData;
492 1757 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
493 1757 : outputData.add(MS::DATA,outputVisCubeHolder);
494 :
495 1757 : if (wideFieldMode_p)
496 : {
497 : // Configure Transformation Engine
498 1757 : WideFieldPhaseShiftingTransformEngine<Complex> transformer(phaseShift_p,&uvw,&frequencies,&inputData,&outputData);
499 :
500 : // Transform data
501 1757 : transformFreqAxis2(vb->getShape(),transformer);
502 1757 : }
503 : else
504 : {
505 : // Configure Transformation Engine
506 0 : PhaseShiftingTransformEngine<Complex> transformer(dx_p,dy_p,&uvw,&frequencies,&inputData,&outputData);
507 :
508 : // Transform data
509 0 : transformFreqAxis2(vb->getShape(),transformer);
510 : }
511 1757 : }
512 :
513 : // -----------------------------------------------------------------------
514 : //
515 : // -----------------------------------------------------------------------
516 4643 : void PhaseShiftingTVI::uvw (casacore::Matrix<double> & uvw) const
517 : {
518 4643 : if (wideFieldMode_p) {
519 4643 : uvw.resize(newUVW_p.shape(),false);
520 4643 : uvw = newUVW_p;
521 : }
522 : else {
523 0 : getVii()->uvw (uvw);
524 : }
525 4643 : }
526 :
527 : //////////////////////////////////////////////////////////////////////////
528 : // PhaseShiftingTVIFactory class
529 : //////////////////////////////////////////////////////////////////////////
530 :
531 : // -----------------------------------------------------------------------
532 : //
533 : // -----------------------------------------------------------------------
534 0 : PhaseShiftingTVIFactory::PhaseShiftingTVIFactory ( Record &configuration,
535 0 : ViImplementation2 *inputVii)
536 : {
537 0 : inputVii_p = inputVii;
538 0 : configuration_p = configuration;
539 0 : }
540 :
541 : // -----------------------------------------------------------------------
542 : //
543 : // -----------------------------------------------------------------------
544 0 : vi::ViImplementation2 * PhaseShiftingTVIFactory::createVi(VisibilityIterator2 *) const
545 : {
546 0 : return new PhaseShiftingTVI(inputVii_p,configuration_p);
547 : }
548 :
549 : // -----------------------------------------------------------------------
550 : //
551 : // -----------------------------------------------------------------------
552 0 : vi::ViImplementation2 * PhaseShiftingTVIFactory::createVi() const
553 : {
554 0 : return new PhaseShiftingTVI(inputVii_p,configuration_p);
555 : }
556 :
557 : //////////////////////////////////////////////////////////////////////////
558 : // PhaseShiftingTVILayerFactory class
559 : //////////////////////////////////////////////////////////////////////////
560 :
561 88 : PhaseShiftingTVILayerFactory::PhaseShiftingTVILayerFactory(Record &configuration) :
562 : ViiLayerFactory(),
563 88 : configuration_p(configuration)
564 88 : {}
565 :
566 : ViImplementation2*
567 88 : PhaseShiftingTVILayerFactory::createInstance(ViImplementation2* vii0) const
568 : {
569 : // Make the PhaseShiftingTVi2, using supplied ViImplementation2, and return it
570 88 : ViImplementation2 *vii = new PhaseShiftingTVI(vii0,configuration_p);
571 87 : return vii;
572 : }
573 :
574 : //////////////////////////////////////////////////////////////////////////
575 : // PhaseShiftingTransformEngine class
576 : //////////////////////////////////////////////////////////////////////////
577 :
578 : // -----------------------------------------------------------------------
579 : //
580 : // -----------------------------------------------------------------------
581 0 : template<class T> PhaseShiftingTransformEngine<T>::PhaseShiftingTransformEngine(
582 : Double dx, Double dy,
583 : Matrix<Double> *uvw,
584 : Vector<Double> *frequencies,
585 : DataCubeMap *inputData,
586 : DataCubeMap *outputData):
587 0 : FreqAxisTransformEngine2<T>(inputData,outputData)
588 : {
589 0 : uvw_p = uvw;
590 0 : frequencies_p = frequencies;
591 :
592 : // Offsets in radians (input is arcsec)
593 0 : dx_p = dx*(M_PI / 180.0 / 3600.0);
594 0 : dy_p = dy*(M_PI / 180.0 / 3600.0);
595 0 : }
596 :
597 : // -----------------------------------------------------------------------
598 : //
599 : // -----------------------------------------------------------------------
600 0 : template<class T> void PhaseShiftingTransformEngine<T>::transform( )
601 : {
602 0 : transformCore(inputData_p,outputData_p);
603 0 : }
604 :
605 : // -----------------------------------------------------------------------
606 : //
607 : // -----------------------------------------------------------------------
608 0 : template<class T> void PhaseShiftingTransformEngine<T>::transformCore( DataCubeMap *inputData,
609 : DataCubeMap *outputData)
610 : {
611 : // Get input/output data
612 0 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
613 0 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
614 :
615 : // Extra path as fraction of U and V in m
616 0 : Double phase = dx_p*(*uvw_p)(0,rowIndex_p) + dy_p*(*uvw_p)(1,rowIndex_p);
617 :
618 : // In radian/Hz
619 0 : phase *= -2.0 * M_PI / C::c;
620 :
621 : // Main loop
622 : Double phase_i;
623 0 : Complex factor;
624 0 : for (uInt chan_i=0;chan_i<inputVector.size();chan_i++)
625 : {
626 0 : phase_i = phase * (*frequencies_p)(chan_i);
627 0 : factor = Complex(cos(phase_i), sin(phase_i));
628 0 : outputVector(chan_i) = factor*inputVector(chan_i);
629 : }
630 0 : }
631 :
632 : //////////////////////////////////////////////////////////////////////////
633 : // WideFieldPhaseShiftingTransformEngine class
634 : //////////////////////////////////////////////////////////////////////////
635 :
636 : // -----------------------------------------------------------------------
637 : //
638 : // -----------------------------------------------------------------------
639 8157 : template<class T> WideFieldPhaseShiftingTransformEngine<T>::WideFieldPhaseShiftingTransformEngine(
640 : const Vector<Double> &phaseShift,
641 : Matrix<Double> *uvw,
642 : Vector<Double> *frequencies,
643 : DataCubeMap *inputData,
644 : DataCubeMap *outputData):
645 8157 : FreqAxisTransformEngine2<T>(inputData,outputData)
646 : {
647 8157 : uvw_p = uvw;
648 8157 : frequencies_p = frequencies;
649 8157 : phaseShift_p = phaseShift;
650 8157 : }
651 :
652 : // -----------------------------------------------------------------------
653 : //
654 : // -----------------------------------------------------------------------
655 3944138 : template<class T> void WideFieldPhaseShiftingTransformEngine<T>::transform( )
656 : {
657 3944138 : transformCore(inputData_p,outputData_p);
658 3944138 : }
659 :
660 : // -----------------------------------------------------------------------
661 : //
662 : // -----------------------------------------------------------------------
663 3944138 : template<class T> void WideFieldPhaseShiftingTransformEngine<T>::transformCore( DataCubeMap *inputData,
664 : DataCubeMap *outputData)
665 : {
666 : // Get input/output data
667 3944138 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
668 3944138 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
669 :
670 3944138 : if (phaseShift_p.shape() == 0) {
671 : // no shift, bypass data as 'passthrough' field
672 38324 : outputVector = inputVector;
673 38324 : return;
674 : }
675 :
676 : // Main loop
677 : Double phase_i;
678 3905814 : Complex factor;
679 11948290 : for (uInt chan_i=0;chan_i<inputVector.size();chan_i++)
680 : {
681 8042476 : phase_i = phaseShift_p(rowIndex_p) * (*frequencies_p)(chan_i);
682 8042476 : factor = Complex(cos(phase_i), sin(phase_i));
683 8042476 : outputVector(chan_i) = factor*inputVector(chan_i);
684 : }
685 : }
686 :
687 :
688 : } //# NAMESPACE VI - END
689 :
690 : } //# NAMESPACE CASA - END
691 :
692 :
|
