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 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 :
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 0 : parseConfiguration(configuration);
51 :
52 0 : initialize();
53 :
54 0 : return;
55 0 : }
56 :
57 : /**
58 : * Get max valid FIELD IDs for this iterator. MSv2 uses the FIELD
59 : * table row index as FIELD ID.
60 : */
61 0 : rownr_t PhaseShiftingTVI::getMaxMSFieldID() const
62 : {
63 0 : const auto &fieldsTable = getVii()->fieldSubtablecols();
64 0 : 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 0 : void PhaseShiftingTVI::parsePhasecenter(const Record &config)
75 : {
76 0 : auto exists = config.fieldNumber ("phasecenter");
77 0 : if (exists < 0) {
78 0 : return;
79 : }
80 :
81 : // phasecenter can be given as a string or as a dict (per-field centers)
82 0 : auto isStr = false;
83 0 : casacore::String phaseCenterStr;
84 : try {
85 0 : config.get(exists, phaseCenterStr);
86 0 : isStr = true;
87 0 : logger_p << LogIO::NORMAL << LogOrigin("PhaseShiftingTVI", __FUNCTION__)
88 : << "Interpreting phasecenter as a single string, to be applied"
89 0 : " to all input (selected) fields."<< LogIO::POST;
90 0 : } 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 0 : logger_p << LogIO::NORMAL << LogOrigin("PhaseShiftingTVI", __FUNCTION__)
94 : << "Interpreting phasecenter as a dictionary of field->center."
95 0 : << "Using following phase centers:" << LogIO::POST;
96 0 : }
97 :
98 0 : if (isStr) {
99 0 : const auto &phaseCenterDir = checkPhaseCenterStr(phaseCenterStr);
100 0 : phaseCenterSpec_p.insert({-1, phaseCenterDir});
101 0 : } else {
102 0 : parsePhasecenterDict(config);
103 : }
104 0 : }
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 0 : void PhaseShiftingTVI::parsePhasecenterDict(const Record &config)
113 : {
114 0 : const auto &inPhasecenter = config.asRecord("phasecenter");
115 0 : if (inPhasecenter.empty()) {
116 0 : throw AipsError("The dictionary 'phasecenter' is empty.");
117 : }
118 :
119 0 : const auto maxMSField = getMaxMSFieldID();
120 0 : std::set<unsigned int> fieldsSeen;
121 : // Go through items in the input dict/record
122 0 : for (unsigned int rid=0; rid < inPhasecenter.nfields(); ++rid) {
123 0 : const std::string fieldStr = inPhasecenter.name(RecordFieldId(rid));
124 0 : const auto fid = std::stoi(fieldStr);
125 0 : if (fid < 0 || static_cast<unsigned int>(fid) > maxMSField) {
126 0 : throw AipsError("Wrong field ID given: " + std::to_string(fid) +
127 0 : ". This MeasurementSet has field IDs between 0 and " +
128 0 : std::to_string(maxMSField));
129 : }
130 0 : if (not fieldsSeen.insert(fid).second) {
131 0 : throw AipsError("Field " + std::to_string(fid) + " is given multiple times");
132 : }
133 :
134 0 : std::string center;
135 : try {
136 0 : 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 0 : const auto strFieldID = std::to_string(fid);
143 0 : const auto &phaseCenterDir = checkPhaseCenterStr(center, strFieldID);
144 0 : phaseCenterSpec_p.insert({fid, phaseCenterDir});
145 0 : }
146 0 : }
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 0 : MDirection PhaseShiftingTVI::checkPhaseCenterStr(const String &phasecenter,
158 : const string &fieldInfo)
159 : {
160 : // casaMDirection requires a variant
161 0 : casac::variant phaseCenterVar(phasecenter);
162 0 : casacore::MDirection phaseCenterDir;
163 0 : if(!casaMDirection(phaseCenterVar, phaseCenterDir)) {
164 0 : throw AipsError("Cannot interpret phase center string as a direction object: "
165 0 : + phasecenter);
166 : } else {
167 0 : const auto myFrame = phaseCenterDir.getRefString();
168 : MDirection::Types mdtype;
169 0 : MDirection::getType(mdtype, myFrame);
170 0 : 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 0 : 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 0 : auto msg = "Phase center '" + std::string(phasecenter) +
189 0 : "' successfully parsed";
190 0 : if (not fieldInfo.empty()) {
191 0 : msg += " for field " + fieldInfo;
192 : }
193 0 : logger_p << LogIO::NORMAL << LogOrigin("PhaseShiftingTVI", __FUNCTION__)
194 0 : << msg << LogIO::POST;
195 0 : }
196 :
197 0 : return phaseCenterDir;
198 0 : }
199 :
200 : // -----------------------------------------------------------------------
201 : //
202 : // -----------------------------------------------------------------------
203 0 : void PhaseShiftingTVI::parseConfiguration(const Record &configuration)
204 : {
205 0 : int exists = -1;
206 :
207 0 : exists = -1;
208 0 : exists = configuration.fieldNumber ("XpcOffset");
209 0 : if (exists >= 0)
210 : {
211 0 : configuration.get (exists, dx_p);
212 : }
213 :
214 0 : exists = -1;
215 0 : exists = configuration.fieldNumber ("YpcOffset");
216 0 : if (exists >= 0)
217 : {
218 0 : configuration.get (exists, dy_p);
219 : }
220 :
221 0 : 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 0 : parsePhasecenter(configuration);
229 0 : wideFieldMode_p = true;
230 0 : }
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 0 : PhaseShiftingTVI::findConvertedPhaseCenter() const
243 : {
244 0 : const auto *vb = getVii()->getVisBuffer();
245 :
246 0 : auto centerIt = phaseCenterSpec_p.find(-1);
247 0 : if (centerIt == phaseCenterSpec_p.end()) {
248 0 : auto fieldID = vb->fieldId()[0];
249 0 : centerIt = phaseCenterSpec_p.find(fieldID);
250 0 : if (centerIt == phaseCenterSpec_p.end()) {
251 0 : return {false, MDirection()};
252 : }
253 : }
254 0 : auto convertedPhaseCenter = centerIt->second;
255 :
256 0 : if (convertedPhaseCenter.getRefString() != vb->phaseCenter().getRefString()) {
257 : MDirection::Types mdtype;
258 0 : MDirection::getType(mdtype, vb->phaseCenter().getRefString());
259 0 : convertedPhaseCenter = MDirection::Convert(convertedPhaseCenter, mdtype)();
260 : }
261 :
262 0 : return {true, convertedPhaseCenter};
263 0 : }
264 :
265 : // -----------------------------------------------------------------------
266 : //
267 : // -----------------------------------------------------------------------
268 0 : void PhaseShiftingTVI::initialize()
269 : {
270 : // Populate nchan input-output maps
271 : Int spw;
272 0 : uInt spw_idx = 0;
273 0 : map<Int,vector<Int> >::iterator iter;
274 0 : for (iter=spwInpChanIdxMap_p.begin();iter!=spwInpChanIdxMap_p.end();iter++)
275 : {
276 0 : spw = iter->first;
277 0 : spwOutChanNumMap_p[spw] = spwInpChanIdxMap_p[spw].size();
278 :
279 0 : 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 0 : if (wideFieldMode_p)
285 : {
286 0 : const auto selectedInputMsCols = std::make_unique<MSColumns>(getVii()->ms());
287 0 : observatoryPosition_p = selectedInputMsCols->antenna().positionMeas()(0);
288 0 : referenceTime_p = selectedInputMsCols->timeMeas()(0);
289 0 : referenceTimeUnits_p = selectedInputMsCols->timeQuant()(0).getUnit();
290 0 : }
291 0 : }
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 0 : void PhaseShiftingTVI::shiftUVWPhases()
298 : {
299 : // Get input VisBuffer
300 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
301 :
302 0 : bool doShift = false;
303 0 : MDirection convertedPhaseCenter;
304 0 : std::tie(doShift, convertedPhaseCenter) = findConvertedPhaseCenter();
305 0 : if (not doShift) {
306 0 : phaseShift_p.resize(0, false);
307 0 : newUVW_p = vb->uvw();
308 0 : return;
309 : }
310 :
311 : // Initialize epoch corresponding to current buffer
312 : // with time reference to the first row in the MS
313 0 : MEpoch epoch(Quantity(vb->time()(0),referenceTimeUnits_p),referenceTime_p.getRef());
314 0 : MeasFrame refFrame(epoch,observatoryPosition_p);
315 0 : UVWMachine uvwMachine(convertedPhaseCenter, vb->phaseCenter(), refFrame,false,false);
316 : // Initialize phase array and uvw matrix
317 0 : phaseShift_p.resize(vb->nRows(),false);
318 0 : newUVW_p.resize(vb->uvw().shape(),false);
319 :
320 : // Obtain phase shift and new uvw coordinates
321 0 : Vector<Double> dummy(3,0.0);
322 0 : double phase2radPerHz = -2.0 * M_PI / C::c;
323 0 : 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 0 : 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 0 : dummy(0) = -1*dummy(0);
331 0 : dummy(1) = -1*dummy(1);
332 :
333 : // Transform uvw coordinates and obtain corresponding phase shift
334 0 : 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 0 : dummy(0) = -1*dummy(0);
339 0 : dummy(1) = -1*dummy(1);
340 0 : newUVW_p.column(row) = dummy;
341 : // Convert phase shift to radian/Hz
342 0 : phaseShift_p(row) = phase2radPerHz*phaseShift_p(row);
343 : }
344 0 : }
345 :
346 : // -----------------------------------------------------------------------
347 : //
348 : // -----------------------------------------------------------------------
349 0 : void PhaseShiftingTVI::origin()
350 : {
351 : // Drive underlying ViImplementation2
352 0 : getVii()->origin();
353 :
354 : // CAS-12706 Add support for shifting across large offset/angles
355 0 : if (wideFieldMode_p) {
356 0 : shiftUVWPhases();
357 : }
358 :
359 : // Define the shapes in the VB2, patch provided by cgarcia in CAS-12706
360 0 : configureShapes();
361 :
362 : // Synchronize own VisBuffer
363 0 : configureNewSubchunk();
364 0 : }
365 :
366 : // -----------------------------------------------------------------------
367 : //
368 : // -----------------------------------------------------------------------
369 0 : void PhaseShiftingTVI::next()
370 : {
371 : // Drive underlying ViImplementation2
372 0 : getVii()->next();
373 :
374 : // CAS-12706 Add support for shifting across large offset/angles
375 0 : if (wideFieldMode_p) {
376 0 : shiftUVWPhases();
377 : }
378 :
379 : // Define the shapes in the VB2, patch provided by cgarcia in CAS-12706
380 0 : configureShapes();
381 :
382 : // Synchronize own VisBuffer
383 0 : configureNewSubchunk();
384 0 : }
385 :
386 :
387 : // -----------------------------------------------------------------------
388 : //
389 : // -----------------------------------------------------------------------
390 0 : void PhaseShiftingTVI::visibilityObserved (Cube<Complex> & vis) const
391 : {
392 : // Get input VisBuffer
393 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
394 0 : Matrix<Double> uvw = vb->uvw();
395 0 : Vector<Double> frequencies = vb->getFrequencies(0);
396 :
397 : // Reshape output data before passing it to the DataCubeHolder
398 0 : vis.resize(getVisBuffer()->getShape(),false);
399 :
400 : // Gather input data
401 0 : DataCubeMap inputData;
402 0 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCube());
403 0 : inputData.add(MS::DATA,inputVisCubeHolder);
404 :
405 : // Gather output data
406 0 : DataCubeMap outputData;
407 0 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
408 0 : outputData.add(MS::DATA,outputVisCubeHolder);
409 :
410 0 : if (wideFieldMode_p)
411 : {
412 : // Configure Transformation Engine
413 0 : WideFieldPhaseShiftingTransformEngine<Complex> transformer(phaseShift_p,&uvw,&frequencies,&inputData,&outputData);
414 :
415 : // Transform data
416 0 : transformFreqAxis2(vb->getShape(),transformer);
417 0 : }
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 0 : }
427 :
428 : // -----------------------------------------------------------------------
429 : //
430 : // -----------------------------------------------------------------------
431 0 : void PhaseShiftingTVI::visibilityCorrected (Cube<Complex> & vis) const
432 : {
433 : // Get input VisBuffer
434 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
435 0 : Matrix<Double> uvw = vb->uvw();
436 0 : Vector<Double> frequencies = vb->getFrequencies(0);
437 :
438 : // Reshape output data before passing it to the DataCubeHolder
439 0 : vis.resize(getVisBuffer()->getShape(),false);
440 :
441 : // Gather input data
442 0 : DataCubeMap inputData;
443 0 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCubeCorrected());
444 0 : inputData.add(MS::DATA,inputVisCubeHolder);
445 :
446 : // Gather output data
447 0 : DataCubeMap outputData;
448 0 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
449 0 : outputData.add(MS::DATA,outputVisCubeHolder);
450 :
451 0 : if (wideFieldMode_p)
452 : {
453 : // Configure Transformation Engine
454 0 : WideFieldPhaseShiftingTransformEngine<Complex> transformer(phaseShift_p,&uvw,&frequencies,&inputData,&outputData);
455 :
456 : // Transform data
457 0 : transformFreqAxis2(vb->getShape(),transformer);
458 0 : }
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 0 : }
468 :
469 : // -----------------------------------------------------------------------
470 : //
471 : // -----------------------------------------------------------------------
472 0 : void PhaseShiftingTVI::visibilityModel (Cube<Complex> & vis) const
473 : {
474 : // Get input VisBuffer
475 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
476 0 : Matrix<Double> uvw = vb->uvw();
477 0 : Vector<Double> frequencies = vb->getFrequencies(0);
478 :
479 : // Reshape output data before passing it to the DataCubeHolder
480 0 : vis.resize(getVisBuffer()->getShape(),false);
481 :
482 : // Gather input data
483 0 : DataCubeMap inputData;
484 0 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCubeModel());
485 0 : inputData.add(MS::DATA,inputVisCubeHolder);
486 :
487 : // Gather output data
488 0 : DataCubeMap outputData;
489 0 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
490 0 : outputData.add(MS::DATA,outputVisCubeHolder);
491 :
492 0 : if (wideFieldMode_p)
493 : {
494 : // Configure Transformation Engine
495 0 : WideFieldPhaseShiftingTransformEngine<Complex> transformer(phaseShift_p,&uvw,&frequencies,&inputData,&outputData);
496 :
497 : // Transform data
498 0 : transformFreqAxis2(vb->getShape(),transformer);
499 0 : }
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 0 : }
509 :
510 : // -----------------------------------------------------------------------
511 : //
512 : // -----------------------------------------------------------------------
513 0 : void PhaseShiftingTVI::uvw (casacore::Matrix<double> & uvw) const
514 : {
515 0 : if (wideFieldMode_p) {
516 0 : uvw.resize(newUVW_p.shape(),false);
517 0 : uvw = newUVW_p;
518 : }
519 : else {
520 0 : getVii()->uvw (uvw);
521 : }
522 0 : }
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 0 : PhaseShiftingTVILayerFactory::PhaseShiftingTVILayerFactory(Record &configuration) :
559 : ViiLayerFactory(),
560 0 : configuration_p(configuration)
561 0 : {}
562 :
563 : ViImplementation2*
564 0 : PhaseShiftingTVILayerFactory::createInstance(ViImplementation2* vii0) const
565 : {
566 : // Make the PhaseShiftingTVi2, using supplied ViImplementation2, and return it
567 0 : ViImplementation2 *vii = new PhaseShiftingTVI(vii0,configuration_p);
568 0 : 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*(M_PI / 180.0 / 3600.0);
591 0 : dy_p = dy*(M_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 * M_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 0 : 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 0 : FreqAxisTransformEngine2<T>(inputData,outputData)
643 : {
644 0 : uvw_p = uvw;
645 0 : frequencies_p = frequencies;
646 0 : phaseShift_p = phaseShift;
647 0 : }
648 :
649 : // -----------------------------------------------------------------------
650 : //
651 : // -----------------------------------------------------------------------
652 0 : template<class T> void WideFieldPhaseShiftingTransformEngine<T>::transform( )
653 : {
654 0 : transformCore(inputData_p,outputData_p);
655 0 : }
656 :
657 : // -----------------------------------------------------------------------
658 : //
659 : // -----------------------------------------------------------------------
660 0 : template<class T> void WideFieldPhaseShiftingTransformEngine<T>::transformCore( DataCubeMap *inputData,
661 : DataCubeMap *outputData)
662 : {
663 : // Get input/output data
664 0 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
665 0 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
666 :
667 0 : if (phaseShift_p.shape() == 0) {
668 : // no shift, bypass data as 'passthrough' field
669 0 : outputVector = inputVector;
670 0 : return;
671 : }
672 :
673 : // Main loop
674 : Double phase_i;
675 0 : Complex factor;
676 0 : for (uInt chan_i=0;chan_i<inputVector.size();chan_i++)
677 : {
678 0 : phase_i = phaseShift_p(rowIndex_p) * (*frequencies_p)(chan_i);
679 0 : factor = Complex(cos(phase_i), sin(phase_i));
680 0 : outputVector(chan_i) = factor*inputVector(chan_i);
681 : }
682 : }
683 :
684 :
685 : } //# NAMESPACE VI - END
686 :
687 : } //# NAMESPACE CASA - END
688 :
689 :
|
