Line data Source code
1 : //# FTMachine.cc: Implementation of FTMachine class
2 : //# Copyright (C) 1997-2016
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 : #include <cmath>
28 : #include <casacore/casa/Quanta/Quantum.h>
29 : #include <casacore/casa/Quanta/UnitMap.h>
30 : #include <casacore/casa/Quanta/UnitVal.h>
31 : #include <casacore/measures/Measures/Stokes.h>
32 : #include <casacore/casa/Quanta/Euler.h>
33 : #include <casacore/casa/Quanta/RotMatrix.h>
34 : #include <casacore/measures/Measures/MFrequency.h>
35 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
36 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
37 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
38 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
39 : #include <casacore/coordinates/Coordinates/Projection.h>
40 : #include <casacore/lattices/Lattices/LatticeLocker.h>
41 : #include <casacore/ms/MeasurementSets/MSColumns.h>
42 : #include <casacore/casa/BasicSL/Constants.h>
43 : #include <synthesis/TransformMachines2/FTMachine.h>
44 : #include <synthesis/TransformMachines2/SkyJones.h>
45 : #include <synthesis/TransformMachines2/VisModelData.h>
46 : #include <synthesis/TransformMachines2/BriggsCubeWeightor.h>
47 : #include <casacore/scimath/Mathematics/RigidVector.h>
48 : #include <synthesis/TransformMachines/StokesImageUtil.h>
49 : #include <synthesis/TransformMachines2/Utils.h>
50 : #include <msvis/MSVis/VisibilityIterator2.h>
51 : #include <msvis/MSVis/VisBuffer2.h>
52 : #include <msvis/MSVis/StokesVector.h>
53 : #include <msvis/MSVis/MSUtil.h>
54 : #include <casacore/images/Images/ImageInterface.h>
55 : #include <casacore/images/Images/PagedImage.h>
56 : #include <casacore/images/Images/ImageUtilities.h>
57 : #include <casacore/casa/Containers/Block.h>
58 : #include <casacore/casa/Containers/Record.h>
59 : #include <casacore/casa/Arrays/ArrayIter.h>
60 : #include <casacore/casa/Arrays/ArrayLogical.h>
61 : #include <casacore/casa/Arrays/ArrayMath.h>
62 : #include <casacore/casa/Arrays/MatrixMath.h>
63 : #include <casacore/casa/Arrays/MaskedArray.h>
64 : #include <casacore/casa/Arrays/Array.h>
65 : #include <casacore/casa/Arrays/Vector.h>
66 : #include <casacore/casa/Arrays/Matrix.h>
67 : #include <casacore/casa/Arrays/MatrixIter.h>
68 : #include <casacore/casa/BasicSL/String.h>
69 : #include <casacore/casa/Utilities/Assert.h>
70 : #include <casacore/casa/Utilities/BinarySearch.h>
71 : #include <casacore/casa/Exceptions/Error.h>
72 : #include <casacore/scimath/Mathematics/NNGridder.h>
73 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
74 : #include <casacore/measures/Measures/UVWMachine.h>
75 :
76 : #include <casacore/casa/System/ProgressMeter.h>
77 :
78 : #include <casacore/casa/OS/Timer.h>
79 : #include <sstream>
80 : #include <iostream>
81 : #include <iomanip>
82 : using namespace casacore;
83 : namespace casa{//# CASA namespace
84 : namespace refim {//# namespace refactor imaging
85 :
86 : using namespace casacore;
87 : using namespace casa;
88 : using namespace casacore;
89 : using namespace casa::refim;
90 : using namespace casacore;
91 : using namespace casa::vi;
92 68 : FTMachine::FTMachine() : isDryRun(false), image(0), uvwMachine_p(0),
93 68 : tangentSpecified_p(false), fixMovingSource_p(false),
94 34 : ephemTableName_p(""),
95 34 : movingDirShift_p(0.0),
96 34 : distance_p(0.0), lastFieldId_p(-1),lastMSId_p(-1), romscol_p(nullptr),
97 34 : useDoubleGrid_p(false),
98 34 : freqFrameValid_p(false),
99 34 : freqInterpMethod_p(InterpolateArray1D<Double,Complex>::nearestNeighbour),
100 34 : pointingDirCol_p("DIRECTION"),
101 34 : cfStokes_p(), cfCache_p(), cfs_p(), cfwts_p(), cfs2_p(), cfwts2_p(),
102 34 : canComputeResiduals_p(false), toVis_p(true),
103 170 : numthreads_p(-1), pbLimit_p(0.05),sj_p(0), cmplxImage_p( ), vbutil_p(), phaseCenterTime_p(-1.0), doneThreadPartition_p(-1), briggsWeightor_p(nullptr), tempFileNames_p(0), ftmType_p(FTMachine::CORRECTED), avgPBReady_p(false)
104 : {
105 34 : spectralCoord_p=SpectralCoordinate();
106 34 : isPseudoI_p=false;
107 34 : spwChanSelFlag_p=0;
108 34 : polInUse_p=0;
109 34 : pop_p = new PolOuterProduct;
110 34 : ft_p=FFT2D(true);
111 34 : }
112 :
113 0 : FTMachine::FTMachine(CountedPtr<CFCache>& cfcache,CountedPtr<ConvolutionFunction>& cf):
114 0 : isDryRun(false), image(0), uvwMachine_p(0),
115 0 : tangentSpecified_p(false), fixMovingSource_p(false),
116 0 : ephemTableName_p(""),
117 0 : movingDirShift_p(0.0),
118 0 : distance_p(0.0), lastFieldId_p(-1),lastMSId_p(-1), romscol_p(nullptr),
119 0 : useDoubleGrid_p(false),
120 0 : freqFrameValid_p(false),
121 0 : freqInterpMethod_p(InterpolateArray1D<Double,Complex>::nearestNeighbour),
122 0 : pointingDirCol_p("DIRECTION"),
123 0 : cfStokes_p(), cfCache_p(cfcache), cfs_p(), cfwts_p(), cfs2_p(), cfwts2_p(),
124 0 : convFuncCtor_p(cf),canComputeResiduals_p(false), toVis_p(true), numthreads_p(-1),
125 0 : pbLimit_p(0.05),sj_p(0), cmplxImage_p( ), vbutil_p(), phaseCenterTime_p(-1.0), doneThreadPartition_p(-1), briggsWeightor_p(nullptr), tempFileNames_p(0), ftmType_p(FTMachine::CORRECTED), avgPBReady_p(false)
126 : {
127 0 : spectralCoord_p=SpectralCoordinate();
128 0 : isPseudoI_p=false;
129 0 : spwChanSelFlag_p=0;
130 0 : polInUse_p=0;
131 0 : pop_p = new PolOuterProduct;
132 0 : ft_p=FFT2D(true);
133 0 : }
134 :
135 204 : LogIO& FTMachine::logIO() {return logIO_p;};
136 :
137 : //----------------------------------------------------------------------
138 0 : FTMachine& FTMachine::operator=(const FTMachine& other)
139 : {
140 0 : if(this!=&other) {
141 0 : image=other.image;
142 : //generic selection stuff and state
143 0 : nAntenna_p=other.nAntenna_p;
144 0 : distance_p=other.distance_p;
145 0 : lastFieldId_p=other.lastFieldId_p;
146 0 : lastMSId_p=other.lastMSId_p;
147 0 : romscol_p=other.romscol_p;
148 0 : tangentSpecified_p=other.tangentSpecified_p;
149 0 : mTangent_p=other.mTangent_p;
150 0 : mImage_p=other.mImage_p;
151 0 : mFrame_p=other.mFrame_p;
152 :
153 0 : nx=other.nx;
154 0 : ny=other.ny;
155 0 : npol=other.npol;
156 0 : nchan=other.nchan;
157 0 : nvischan=other.nvischan;
158 0 : nvispol=other.nvispol;
159 0 : mLocation_p=other.mLocation_p;
160 0 : if(uvwMachine_p)
161 0 : delete uvwMachine_p;
162 0 : if(other.uvwMachine_p)
163 0 : uvwMachine_p=new casacore::UVWMachine(*other.uvwMachine_p);
164 : else
165 0 : uvwMachine_p=0;
166 0 : doUVWRotation_p=other.doUVWRotation_p;
167 : //Spectral and pol stuff
168 0 : freqInterpMethod_p=other.freqInterpMethod_p;
169 0 : spwChanSelFlag_p.resize();
170 0 : spwChanSelFlag_p=other.spwChanSelFlag_p;
171 0 : freqFrameValid_p=other.freqFrameValid_p;
172 : //selectedSpw_p.resize();
173 : //selectedSpw_p=other.selectedSpw_p;
174 0 : imageFreq_p.resize();
175 0 : imageFreq_p=other.imageFreq_p;
176 0 : lsrFreq_p.resize();
177 0 : lsrFreq_p=other.lsrFreq_p;
178 0 : interpVisFreq_p.resize();
179 0 : interpVisFreq_p=other.interpVisFreq_p;
180 : //multiChanMap_p=other.multiChanMap_p;
181 0 : chanMap.resize();
182 0 : chanMap=other.chanMap;
183 0 : polMap.resize();
184 0 : polMap=other.polMap;
185 0 : nVisChan_p.resize();
186 0 : nVisChan_p=other.nVisChan_p;
187 0 : spectralCoord_p=other.spectralCoord_p;
188 0 : visPolMap_p.resize();
189 0 : visPolMap_p=other.visPolMap_p;
190 : //doConversion_p.resize();
191 : //doConversion_p=other.doConversion_p;
192 0 : pointingDirCol_p=other.pointingDirCol_p;
193 : //moving source stuff
194 0 : movingDir_p=other.movingDir_p;
195 0 : fixMovingSource_p=other.fixMovingSource_p;
196 0 : ephemTableName_p = other.ephemTableName_p;
197 0 : firstMovingDir_p=other.firstMovingDir_p;
198 0 : movingDirShift_p=other.movingDirShift_p;
199 : //Double precision gridding for those FTMachines that can do
200 0 : useDoubleGrid_p=other.useDoubleGrid_p;
201 0 : cfStokes_p = other.cfStokes_p;
202 0 : polInUse_p = other.polInUse_p;
203 0 : cfs_p = other.cfs_p;
204 0 : cfwts_p = other.cfwts_p;
205 0 : cfs2_p = other.cfs2_p;
206 0 : cfwts2_p = other.cfwts2_p;
207 0 : canComputeResiduals_p = other.canComputeResiduals_p;
208 :
209 0 : pop_p = other.pop_p;
210 0 : toVis_p = other.toVis_p;
211 0 : spwFreqSel_p.resize();
212 0 : spwFreqSel_p = other.spwFreqSel_p;
213 0 : expandedSpwFreqSel_p = other.expandedSpwFreqSel_p;
214 0 : expandedSpwConjFreqSel_p = other.expandedSpwConjFreqSel_p;
215 0 : cmplxImage_p=other.cmplxImage_p;
216 0 : vbutil_p=other.vbutil_p;
217 0 : numthreads_p=other.numthreads_p;
218 0 : pbLimit_p=other.pbLimit_p;
219 0 : convFuncCtor_p = other.convFuncCtor_p;
220 0 : sj_p.resize();
221 0 : sj_p=other.sj_p;
222 0 : isDryRun=other.isDryRun;
223 0 : phaseCenterTime_p=other.phaseCenterTime_p;
224 0 : doneThreadPartition_p=other.doneThreadPartition_p;
225 0 : xsect_p=other.xsect_p;
226 0 : ysect_p=other.ysect_p;
227 0 : nxsect_p=other.nxsect_p;
228 0 : nysect_p=other.nysect_p;
229 0 : obsvelconv_p=other.obsvelconv_p;
230 0 : mtype_p=other.mtype_p;
231 0 : briggsWeightor_p=other.briggsWeightor_p;
232 0 : ft_p=other.ft_p;
233 0 : ftmType_p = other.ftmType_p;
234 0 : avgPBReady_p = other.avgPBReady_p;
235 : };
236 0 : return *this;
237 : };
238 :
239 0 : FTMachine* FTMachine::cloneFTM(){
240 0 : Record rec;
241 0 : String err;
242 0 : if(!(this->toRecord(err, rec)))
243 0 : throw(AipsError("Error in cloning FTMachine"));
244 0 : FTMachine* retval=VisModelData::NEW_FT(rec);
245 0 : if(retval)
246 0 : retval->briggsWeightor_p=briggsWeightor_p;
247 0 : return retval;
248 0 : }
249 :
250 : //----------------------------------------------------------------------
251 0 : Bool FTMachine::changed(const vi::VisBuffer2&) {
252 0 : return false;
253 : }
254 : //----------------------------------------------------------------------
255 0 : FTMachine::FTMachine(const FTMachine& other)
256 : {
257 0 : operator=(other);
258 0 : }
259 :
260 0 : Bool FTMachine::doublePrecGrid(){
261 0 : return useDoubleGrid_p;
262 : }
263 :
264 0 : void FTMachine::reset(){
265 : //ft_p=FFT2D(true);
266 0 : }
267 :
268 : //----------------------------------------------------------------------
269 34 : void FTMachine::initPolInfo(const vi::VisBuffer2& vb)
270 : {
271 : //
272 : // Need to figure out where to compute the following arrays/ints
273 : // in the re-factored code.
274 : // ----------------------------------------------------------------
275 : {
276 34 : polInUse_p = 0;
277 34 : uInt N=0;
278 102 : for(uInt i=0;i<polMap.nelements();i++) if (polMap(i) > -1) polInUse_p++;
279 34 : cfStokes_p.resize(polInUse_p);
280 102 : for(uInt i=0;i<polMap.nelements();i++)
281 68 : if (polMap(i) > -1) {cfStokes_p(N) = vb.correlationTypes()(i);N++;}
282 : }
283 34 : }
284 : //----------------------------------------------------------------------
285 34 : void FTMachine::initMaps(const vi::VisBuffer2& vb) {
286 34 : logIO() << LogOrigin("FTMachine", "initMaps") << LogIO::NORMAL;
287 :
288 34 : AlwaysAssert(image, AipsError);
289 :
290 : // Set the frame for the UVWMachine
291 34 : if(vb.isAttached()){
292 : //mFrame_p=MeasFrame(MEpoch(Quantity(vb.time()(0), "s"), MSColumns(vb.ms()).timeMeas()(0).getRef()), mLocation_p);
293 34 : if(vbutil_p.null())
294 17 : vbutil_p=new VisBufferUtil(vb);
295 34 : romscol_p=new MSColumns(vb.ms());
296 68 : Unit epochUnit=(romscol_p->time()).keywordSet().asArrayString("QuantumUnits")(IPosition(1,0));
297 34 : if(!mFrame_p.epoch())
298 17 : mFrame_p.set(MEpoch(Quantity(vb.time()(0), epochUnit), (romscol_p->timeMeas())(0).getRef()));
299 : else
300 17 : mFrame_p.resetEpoch(MEpoch(Quantity(vb.time()(0), epochUnit), (romscol_p->timeMeas())(0).getRef()));
301 34 : if(!mFrame_p.position())
302 17 : mFrame_p.set(mLocation_p);
303 : else
304 17 : mFrame_p.resetPosition(mLocation_p);
305 34 : if(!mFrame_p.direction())
306 17 : mFrame_p.set(vbutil_p->getEphemDir(vb, phaseCenterTime_p));
307 : else
308 17 : mFrame_p.resetDirection(vbutil_p->getEphemDir(vb, phaseCenterTime_p));
309 34 : }
310 : else{
311 0 : throw(AipsError("Cannot define some frame as no Visiter/MS is attached"));
312 : }
313 : //////TESTOOOO
314 : ///setMovingSource("COMET", "des_deedee_ephem2.tab");
315 : ///////////////////////////////////////////
316 : // First get the CoordinateSystem for the image and then find
317 : // the DirectionCoordinate
318 34 : casacore::CoordinateSystem coords=image->coordinates();
319 34 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
320 34 : AlwaysAssert(directionIndex>=0, AipsError);
321 : DirectionCoordinate
322 34 : directionCoord=coords.directionCoordinate(directionIndex);
323 34 : Int spectralIndex=coords.findCoordinate(Coordinate::SPECTRAL);
324 34 : AlwaysAssert(spectralIndex>-1, AipsError);
325 34 : spectralCoord_p=coords.spectralCoordinate(spectralIndex);
326 :
327 : // get the first position of moving source
328 34 : if(fixMovingSource_p){
329 : //cerr << "obsinfo time " << coords.obsInfo().obsDate() << " epoch used in frame " << MEpoch((mFrame_p.epoch())) << endl;
330 : ///Darn vb.time()(0) may not be the earliest time due to sort issues...
331 : //so lets try to use the same
332 : ///time as SynthesisIUtilMethods::buildCoordinateSystemCore is using
333 : //mFrame_p.resetEpoch(romscol_p->timeMeas()(0));
334 0 : mFrame_p.resetEpoch(coords.obsInfo().obsDate());
335 : //Double firstTime=romscol_p->time()(0);
336 :
337 0 : Double firstTime=coords.obsInfo().obsDate().get("s").getValue();
338 : //First convert to HA-DEC or AZEL for parallax correction
339 0 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
340 0 : MDirection tmphadec;
341 0 : if (upcase(movingDir_p.getRefString()).contains("APP")) {
342 : //cerr << "phaseCenterTime_p " << phaseCenterTime_p << endl;
343 0 : tmphadec = MDirection::Convert((vbutil_p->getEphemDir(vb, phaseCenterTime_p > 0.0 ? phaseCenterTime_p : firstTime)), outref1)();
344 0 : MeasComet mcomet(Path((romscol_p->field()).ephemPath(vb.fieldId()(0))).absoluteName());
345 0 : if (mFrame_p.comet())
346 0 : mFrame_p.resetComet(mcomet);
347 : else
348 0 : mFrame_p.set(mcomet);
349 0 : } else if (upcase(movingDir_p.getRefString()).contains("COMET")) {
350 0 : MeasComet mcomet(Path(ephemTableName_p).absoluteName());
351 0 : if (mFrame_p.comet())
352 0 : mFrame_p.resetComet(mcomet);
353 : else
354 0 : mFrame_p.set(mcomet);
355 0 : tmphadec = MDirection::Convert(MDirection(MDirection::COMET), outref1)();
356 0 : } else {
357 0 : tmphadec = MDirection::Convert(movingDir_p, outref1)();
358 : }
359 0 : MDirection::Ref outref(directionCoord.directionType(), mFrame_p);
360 0 : firstMovingDir_p=MDirection::Convert(tmphadec, outref)();
361 : ////////////////////
362 : /*ostringstream ss;
363 : Unit epochUnit=(romscol_p->time()).keywordSet().asArrayString("QuantumUnits")(IPosition(1,0));
364 : MEpoch(Quantity(vb.time()(0), epochUnit), (romscol_p->timeMeas())(0).getRef()).print(ss);
365 : cerr << std::setprecision(15) << "First time " << ss.str() << "field id " << vb.fieldId()(0) << endl;
366 : ss.clear();
367 : firstMovingDir_p.print(ss);
368 : cerr << "firstdir " << ss.str() << " " << firstMovingDir_p.toString() << endl;
369 : */
370 : //////////////
371 :
372 0 : if(spectralCoord_p.frequencySystem(False)==MFrequency::REST){
373 : ///We want the data frequency to be shifted to the SOURCE frame
374 : ///which is labelled REST as we have never defined the SOURCE frame didn't we
375 0 : initSourceFreqConv();
376 : }
377 : ///TESTOO
378 : ///waiting for CAS-11060
379 : //firstMovingDir_p=MDirection::Convert(vbutil_p->getPhaseCenter(vb, phaseCenterTime_p), outref)();
380 : ////////////////////
381 0 : }
382 :
383 :
384 : // Now we need MDirection of the image phase center. This is
385 : // what we define it to be. So we define it to be the
386 : // center pixel. So we have to do the conversion here.
387 : // This is independent of padding since we just want to know
388 : // what the world coordinates are for the phase center
389 : // pixel
390 : {
391 34 : Vector<Double> pixelPhaseCenter(2);
392 34 : pixelPhaseCenter(0) = Double( image->shape()(0) / 2 );
393 34 : pixelPhaseCenter(1) = Double( image->shape()(1) / 2 );
394 34 : directionCoord.toWorld(mImage_p, pixelPhaseCenter);
395 34 : }
396 :
397 : // Get the object distance in meters
398 34 : Record info(image->miscInfo());
399 34 : if(info.isDefined("distance")) {
400 0 : info.get("distance", distance_p);
401 0 : if(abs(distance_p)>0.0) {
402 0 : logIO() << "Distance to object is set to " << distance_p/1000.0
403 0 : << "km: applying focus correction" << LogIO::POST;
404 : }
405 : }
406 :
407 : // Set up the UVWMachine.
408 34 : initUVWMachine(vb);
409 :
410 34 : lastFieldId_p=-1;
411 :
412 34 : lastMSId_p=vb.msId();
413 :
414 : // Set up maps
415 :
416 :
417 :
418 : //Store the image/grid channels freq values
419 : {
420 34 : Int chanNumbre=image->shape()(3);
421 34 : Vector<Double> pixindex(chanNumbre);
422 34 : imageFreq_p.resize(chanNumbre);
423 34 : Vector<Double> tempStorFreq(chanNumbre);
424 34 : indgen(pixindex);
425 : // pixindex=pixindex+1.0;
426 68 : for (Int ll=0; ll< chanNumbre; ++ll){
427 34 : if( !spectralCoord_p.toWorld(tempStorFreq(ll), pixindex(ll))){
428 0 : logIO() << "Cannot get imageFreq " << LogIO::EXCEPTION;
429 :
430 : }
431 : }
432 34 : convertArray(imageFreq_p,tempStorFreq);
433 34 : }
434 : //Destroy any conversion layer Freq coord if freqframe is not valid
435 34 : if(!freqFrameValid_p){
436 0 : MFrequency::Types imageFreqType=spectralCoord_p.frequencySystem();
437 0 : spectralCoord_p.setFrequencySystem(imageFreqType);
438 0 : spectralCoord_p.setReferenceConversion(imageFreqType,
439 0 : MEpoch(Quantity(vb.time()(0), "s")),
440 0 : mLocation_p,
441 0 : mImage_p);
442 : }
443 :
444 : // Channel map: do this properly by looking up the frequencies
445 : // If a visibility channel does not map onto an image
446 : // pixel then we set the corresponding chanMap to -1.
447 : // This means that put and get must always check for this
448 : // value (see e.g. GridFT)
449 :
450 34 : nvischan = vb.getFrequencies(0).nelements();
451 34 : interpVisFreq_p.resize();
452 34 : interpVisFreq_p=vb.getFrequencies(0);
453 :
454 : // Polarization map
455 34 : visPolMap_p.resize();
456 34 : polMap.resize();
457 :
458 : //As matchChannel calls matchPol ...it has to be called after making sure
459 : //polMap and visPolMap are zero size to force a polMap matching
460 34 : chanMap.resize();
461 34 : matchChannel(vb);
462 : //chanMap=multiChanMap_p[vb.spectralWindows()(0)];
463 34 : if(chanMap.nelements() == 0)
464 0 : chanMap=Vector<Int>(vb.getFrequencies(0).nelements(), -1);
465 :
466 : {
467 : //logIO() << LogIO::DEBUGGING << "Channel Map: " << chanMap << LogIO::POST;
468 : }
469 : // Should never get here
470 34 : if(max(chanMap)>=nchan||min(chanMap)<-2) {
471 0 : logIO() << "Illegal Channel Map: " << chanMap << LogIO::EXCEPTION;
472 : }
473 :
474 :
475 34 : initPolInfo(vb);
476 34 : Vector<Int> intpolmap(visPolMap_p.nelements());
477 102 : for (uInt kk=0; kk < intpolmap.nelements(); ++kk){
478 68 : intpolmap[kk]=Int(visPolMap_p[kk]);
479 : }
480 34 : pop_p->initCFMaps(intpolmap, polMap);
481 :
482 :
483 :
484 :
485 :
486 :
487 :
488 34 : }
489 :
490 34 : void FTMachine::initUVWMachine(const vi::VisBuffer2& vb) {
491 : // Decide if uvwrotation is not necessary, if phasecenter and
492 : // image center are with in one pixel distance; Save some
493 : // computation time especially for spectral cubes.
494 34 : casacore::CoordinateSystem coords=image->coordinates();
495 34 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
496 34 : AlwaysAssert(directionIndex>=0, AipsError);
497 34 : auto const directionCoord=coords.directionCoordinate(directionIndex);
498 68 : Vector<Double> equal= (mImage_p.getAngle()-
499 102 : vbutil_p->getPhaseCenter(vb, phaseCenterTime_p).getAngle()).getValue();
500 102 : if((abs(equal(0)) < abs(directionCoord.increment()(0)))
501 102 : && (abs(equal(1)) < abs(directionCoord.increment()(1)))){
502 34 : doUVWRotation_p=false;
503 : } else {
504 0 : doUVWRotation_p=true;
505 : }
506 :
507 34 : if(uvwMachine_p) delete uvwMachine_p; uvwMachine_p=0;
508 34 : String observatory;
509 34 : if(vb.isAttached())
510 34 : observatory=(vb.subtableColumns().observation()).telescopeName()(0);
511 : else
512 0 : throw(AipsError("Cannot define frame because of no access to OBSERVATION table"));
513 68 : if(observatory.contains("ATCA") || observatory.contains("DRAO")
514 68 : || observatory.contains("WSRT")){
515 0 : uvwMachine_p=new casacore::UVWMachine(mImage_p, vbutil_p->getPhaseCenter(vb, phaseCenterTime_p), mFrame_p,
516 0 : true, false);
517 : } else {
518 68 : uvwMachine_p=new casacore::UVWMachine(mImage_p, vbutil_p->getPhaseCenter(vb, phaseCenterTime_p), mFrame_p,
519 34 : false, tangentSpecified_p);
520 : }
521 34 : AlwaysAssert(uvwMachine_p, AipsError);
522 :
523 34 : phaseShifter_p=new UVWMachine(*uvwMachine_p);
524 34 : }
525 :
526 34 : void FTMachine::initBriggsWeightor(vi::VisibilityIterator2& vi){
527 : ///Lastly initialized Briggs cube weighting scheme
528 34 : if(!briggsWeightor_p.null()){
529 0 : String error;
530 0 : Record rec;
531 0 : AlwaysAssert(image, AipsError);
532 0 : if(!toRecord(error, rec))
533 0 : throw (AipsError("Could not initialize BriggsWeightor"));
534 0 : String wgtcolname=briggsWeightor_p->initImgWeightCol(vi, *image, rec);
535 0 : tempFileNames_p.resize(tempFileNames_p.nelements()+1, True);
536 0 : tempFileNames_p[tempFileNames_p.nelements()-1]=wgtcolname;
537 :
538 0 : }
539 34 : }
540 :
541 34 : FTMachine::~FTMachine()
542 : {
543 34 : if(uvwMachine_p) delete uvwMachine_p; uvwMachine_p=0;
544 34 : }
545 :
546 :
547 0 : void FTMachine::initSourceFreqConv(){
548 0 : MRadialVelocity::Types refvel=MRadialVelocity::GEO;
549 0 : if(mFrame_p.comet()){
550 : //Has a ephem table
551 0 : if(((mFrame_p.comet())->getTopo().getLength("km").getValue()) > 1.0e-3){
552 0 : refvel=MRadialVelocity::TOPO;
553 : }
554 :
555 :
556 : }
557 : else{
558 : //using a canned DE-200 or 405 source
559 0 : MDirection::Types planetType=MDirection::castType(movingDir_p.getRef().getType());
560 0 : mtype_p=MeasTable::BARYEARTH;
561 0 : if(planetType >=MDirection::MERCURY && planetType <MDirection::COMET){
562 : //Damn these enums are not in the same order
563 0 : switch(planetType){
564 0 : case MDirection::MERCURY :
565 0 : mtype_p=MeasTable::MERCURY;
566 0 : break;
567 0 : case MDirection::VENUS :
568 0 : mtype_p=MeasTable::VENUS;
569 0 : break;
570 0 : case MDirection::MARS :
571 0 : mtype_p=MeasTable::MARS;
572 0 : break;
573 0 : case MDirection::JUPITER :
574 0 : mtype_p=MeasTable::JUPITER;
575 0 : break;
576 0 : case MDirection::SATURN :
577 0 : mtype_p=MeasTable::SATURN;
578 0 : break;
579 0 : case MDirection::URANUS :
580 0 : mtype_p=MeasTable::URANUS;
581 0 : break;
582 0 : case MDirection::NEPTUNE :
583 0 : mtype_p=MeasTable::NEPTUNE;
584 0 : break;
585 0 : case MDirection::PLUTO :
586 0 : mtype_p=MeasTable::PLUTO;
587 0 : break;
588 0 : case MDirection::MOON :
589 0 : mtype_p=MeasTable::MOON;
590 0 : break;
591 0 : case MDirection::SUN :
592 0 : mtype_p=MeasTable::SUN;
593 0 : break;
594 0 : default:
595 0 : throw(AipsError("Cannot translate to known major solar system object"));
596 : }
597 :
598 : }
599 :
600 : }
601 0 : obsvelconv_p=MRadialVelocity::Convert (MRadialVelocity(MVRadialVelocity(0.0),
602 0 : MRadialVelocity::Ref(MRadialVelocity::TOPO, mFrame_p)),
603 0 : MRadialVelocity::Ref(refvel));
604 :
605 0 : }
606 :
607 :
608 51 : Long FTMachine::estimateRAM(const CountedPtr<SIImageStore>& imstor){
609 : //not set up yet
610 51 : if(!image && !imstor)
611 0 : return -1;
612 51 : Long npixels=0;
613 51 : if(image)
614 0 : npixels=((image->shape()).product())/1024;
615 : else{
616 51 : if((imstor->getShape()).product() !=0)
617 51 : npixels=(imstor->getShape()).product()/1024;
618 : }
619 51 : if(npixels==0) npixels=1; //1 kPixels is minimum then
620 51 : Long factor=sizeof(Complex);
621 51 : if(!toVis_p && useDoubleGrid_p)
622 0 : factor=sizeof(DComplex);
623 51 : return (npixels*factor);
624 : }
625 :
626 0 : void FTMachine::shiftFreqToSource(Vector<Double>& freqs){
627 0 : MDoppler dopshift;
628 0 : MEpoch ep(mFrame_p.epoch());
629 0 : if(mFrame_p.comet()){
630 : ////Will use UT for now for ephem tables as it is not clear that they are being
631 : ///filled with TDB as intended in MeasComet.h
632 0 : MEpoch::Convert toUT(ep, MEpoch::UT);
633 0 : MVRadialVelocity cometvel;
634 0 : (*mFrame_p.comet()).getRadVel(cometvel, toUT(ep).get("d").getValue());
635 : //cerr << std::setprecision(10) << "UT " << toUT(ep).get("d").getValue() << " cometvel " << cometvel.get("km/s").getValue("km/s") << endl;
636 :
637 : //cerr << "pos " << MPosition(mFrame_p.position()) << " obsevatory vel " << obsvelconv_p().get("km/s").getValue("km/s") << endl;
638 0 : dopshift=MDoppler(Quantity(-cometvel.get("km/s").getValue("km/s")+obsvelconv_p().get("km/s").getValue("km/s") , "km/s"), MDoppler::RELATIVISTIC);
639 :
640 0 : }
641 : else{
642 0 : Vector<Double> planetparam;
643 0 : Vector<Double> earthparam;
644 0 : MEpoch::Convert toTDB(ep, MEpoch::TDB);
645 0 : earthparam=MeasTable::Planetary(MeasTable::EARTH, toTDB(ep).get("d").getValue());
646 0 : planetparam=MeasTable::Planetary(mtype_p, toTDB(ep).get("d").getValue());
647 : //GEOcentric param
648 0 : planetparam=planetparam-earthparam;
649 0 : Vector<Double> unitdirvec(3);
650 0 : Double dist=sqrt(planetparam(0)*planetparam(0)+planetparam(1)*planetparam(1)+planetparam(2)*planetparam(2));
651 0 : unitdirvec(0)=planetparam(0)/dist;
652 0 : unitdirvec(1)=planetparam(1)/dist;
653 0 : unitdirvec(2)=planetparam(2)/dist;
654 0 : Quantity planetradvel(planetparam(3)*unitdirvec(0)+planetparam(4)*unitdirvec(1)+planetparam(5)*unitdirvec(2), "AU/d");
655 0 : dopshift=MDoppler(Quantity(-planetradvel.getValue("km/s")+obsvelconv_p().get("km/s").getValue("km/s") , "km/s"), MDoppler::RELATIVISTIC);
656 :
657 0 : }
658 :
659 0 : Vector<Double> newfreqs=dopshift.shiftFrequency(freqs);
660 0 : freqs=newfreqs;
661 0 : }
662 :
663 6120 : Bool FTMachine::interpolateFrequencyTogrid(const vi::VisBuffer2& vb,
664 : const Matrix<Float>& wt,
665 : Cube<Complex>& data,
666 : Cube<Int>& flags,
667 : Matrix<Float>& weight,
668 : FTMachine::Type type){
669 6120 : Cube<Complex> origdata;
670 6120 : Cube<Bool> modflagCube;
671 : // Read flags from the vb.
672 6120 : setSpectralFlag(vb,modflagCube);
673 6120 : Vector<Double> visFreq(vb.getFrequencies(0).nelements());
674 : //if(doConversion_p[vb.spectralWindows()[0]]){
675 6120 : if(freqFrameValid_p){
676 6120 : visFreq.resize(lsrFreq_p.shape());
677 6120 : convertArray(visFreq, lsrFreq_p);
678 : }
679 : else{
680 0 : convertArray(visFreq, vb.getFrequencies(0));
681 0 : lsrFreq_p.resize();
682 0 : lsrFreq_p=vb.getFrequencies(0);
683 : }
684 6120 : if(type==FTMachine::MODEL){
685 0 : origdata.reference(vb.visCubeModel());
686 : }
687 6120 : else if(type==FTMachine::CORRECTED){
688 2340 : origdata.reference(vb.visCubeCorrected());
689 : }
690 3780 : else if(type==FTMachine::OBSERVED){
691 720 : origdata.reference(vb.visCube());
692 : }
693 3060 : else if(type==FTMachine::PSF){
694 : // make sure its a size 0 data ...psf
695 : //so avoid reading any data from disk
696 3060 : origdata.resize();
697 :
698 : }
699 : else{
700 0 : throw(AipsError("Don't know which column is being regridded"));
701 : }
702 6120 : if((imageFreq_p.nelements()==1) || (freqInterpMethod_p== InterpolateArray1D<Double, Complex>::nearestNeighbour) || (vb.nChannels()==1)){
703 6120 : data.reference(origdata);
704 : // do something here for apply flag based on spw chan sels
705 : // e.g.
706 :
707 :
708 6120 : flags.resize(modflagCube.shape());
709 6120 : flags=0;
710 : //flags(vb.flagCube())=true;
711 :
712 6120 : flags(modflagCube)=true;
713 :
714 6120 : weight.reference(wt);
715 6120 : interpVisFreq_p.resize();
716 6120 : interpVisFreq_p=lsrFreq_p;
717 : //cerr << "INTERPTOGRID " << interpVisFreq_p.nelements() << " vb.nchan " << vb.nChannels() << endl;
718 6120 : return false;
719 : }
720 :
721 0 : Cube<Bool>flag;
722 :
723 : //okay at this stage we have at least 2 channels
724 0 : Double width=fabs(imageFreq_p[1]-imageFreq_p[0])/fabs(visFreq[1]-visFreq[0]);
725 : //if width is smaller than number of points needed for interpolation ...do it directly
726 : //
727 : // If image chan width is more than twice the data chan width, make a new list of
728 : // data frequencies on which to interpolate. This new list is sync'd with the starting image chan
729 : // and have the same width as the data chans.
730 : /*if(((width >2.0) && (freqInterpMethod_p==InterpolateArray1D<Double, Complex>::linear)) ||
731 : ((width >4.0) && (freqInterpMethod_p !=InterpolateArray1D<Double, Complex>::linear))){
732 : */
733 0 : if(width > 1.0){
734 0 : Double minVF=min(visFreq);
735 0 : Double maxVF=max(visFreq);
736 0 : Double minIF=min(imageFreq_p);
737 0 : Double maxIF=max(imageFreq_p);
738 0 : if( ((minIF-fabs(imageFreq_p[1]-imageFreq_p[0])/2.0) > maxVF) ||
739 0 : ((maxIF+fabs(imageFreq_p[1]-imageFreq_p[0])/2.0) < minVF)){
740 : //This function should not have been called with image
741 : //being out of bound of data...but still
742 0 : interpVisFreq_p.resize(imageFreq_p.nelements());
743 0 : interpVisFreq_p=imageFreq_p;
744 0 : chanMap.resize(interpVisFreq_p.nelements());
745 0 : chanMap.set(-1);
746 : }
747 : else{ // Make a new list of frequencies.
748 : Bool found;
749 0 : uInt where=0;
750 : //Double interpwidth=visFreq[1]-visFreq[0];
751 0 : Double interpwidth=copysign(fabs(imageFreq_p[1]-imageFreq_p[0])/floor(width), visFreq[1]-visFreq[0]);
752 : //if(name() != "GridFT")
753 : // cerr << "width " << width << " interpwidth " << interpwidth << endl;
754 0 : if(minIF < minVF){ // Need to find the first image-channel with data in it
755 0 : where=binarySearchBrackets(found, imageFreq_p, minVF, imageFreq_p.nelements());
756 0 : if(where != imageFreq_p.nelements()){
757 0 : minIF=imageFreq_p[where];
758 : }
759 : }
760 :
761 0 : if(maxIF > maxVF){
762 0 : where=binarySearchBrackets(found, imageFreq_p, maxVF, imageFreq_p.nelements());
763 0 : if(where!= imageFreq_p.nelements()){
764 0 : maxIF=imageFreq_p[where];
765 : }
766 :
767 : }
768 :
769 : // This new list of frequencies starts at the first image channel minus half image channel.
770 : // It ends at the last image channel plus half image channel.
771 0 : Int ninterpchan=(Int)ceil((maxIF-minIF+fabs(imageFreq_p[1]-imageFreq_p[0]))/fabs(interpwidth))+2;
772 0 : chanMap.resize(ninterpchan);
773 0 : chanMap.set(-1);
774 0 : interpVisFreq_p.resize(ninterpchan);
775 0 : interpVisFreq_p[0]=(interpwidth > 0) ? minIF : maxIF;
776 0 : if(freqInterpMethod_p==InterpolateArray1D<Double, Complex>::linear)
777 0 : interpVisFreq_p[0]-=interpwidth;
778 0 : if(freqInterpMethod_p==InterpolateArray1D<Double, Complex>::cubic)
779 0 : interpVisFreq_p[0]-=2.0*interpwidth;
780 0 : Double startedge=abs(imageFreq_p[1]-imageFreq_p[0])/2.0 -abs(interpwidth)/2.0;
781 0 : interpVisFreq_p[0] =(interpwidth >0) ? (interpVisFreq_p[0]-startedge):(interpVisFreq_p[0]+startedge);
782 :
783 0 : for (Int k=1; k < ninterpchan; ++k){
784 0 : interpVisFreq_p[k] = interpVisFreq_p[k-1]+ interpwidth;
785 : }
786 0 : Double halfdiff=fabs((imageFreq_p[1]-imageFreq_p[0])/2.0);
787 0 : for (Int k=0; k < ninterpchan; ++k){
788 : ///chanmap with width
789 : // Double nearestchanval = interpVisFreq_p[k]- (imageFreq_p[1]-imageFreq_p[0])/2.0;
790 : //where=binarySearchBrackets(found, imageFreq_p, nearestchanval, imageFreq_p.nelements());
791 0 : Int which=-1;
792 0 : for (Int j=0; j< Int(imageFreq_p.nelements()); ++j){
793 : //cerr << (imageFreq_p[j]-halfdiff) << " " << (imageFreq_p[j]+halfdiff) << " val " << interpVisFreq_p[k] << endl;
794 0 : if( (interpVisFreq_p[k] >= (imageFreq_p[j]-halfdiff)) && (interpVisFreq_p[k] < (imageFreq_p[j]+halfdiff)))
795 0 : which=j;
796 : }
797 0 : if((which > -1) && (which < Int(imageFreq_p.nelements()))){
798 0 : chanMap[k]=which;
799 : }
800 : else{
801 : //if(name() != "GridFT")
802 : // cerr << "MISSED it " << interpVisFreq_p[k] << endl;
803 : }
804 :
805 :
806 : }
807 : // if(name() != "GridFT")
808 : // cerr << std::setprecision(10) << "chanMap " << chanMap << endl; //" interpvisfreq " << interpVisFreq_p << " orig " << visFreq << endl;
809 :
810 : }// By now, we have a new list of frequencies, synchronized with image channels, but with data chan widths.
811 : }// end of ' if (we have to make new frequencies) '
812 : else{
813 : // Interpolate directly onto output image frequencies.
814 0 : interpVisFreq_p.resize(imageFreq_p.nelements());
815 0 : convertArray(interpVisFreq_p, imageFreq_p);
816 0 : chanMap.resize(interpVisFreq_p.nelements());
817 0 : indgen(chanMap);
818 : }
819 0 : if(type != FTMachine::PSF){ // Interpolating the data
820 : //Need to get new interpolate functions that interpolate explicitly on the 2nd axis
821 : //2 swap of axes needed
822 0 : Cube<Complex> flipdata;
823 0 : Cube<Bool> flipflag;
824 :
825 : // Interpolate the data.
826 : // Input flags are from the previous step ( setSpectralFlag ).
827 : // Output flags contain info about channels that could not be interpolated
828 : // (for example, linear interp with only one data point)
829 0 : swapyz(flipflag,modflagCube);
830 0 : swapyz(flipdata,origdata);
831 : InterpolateArray1D<Double,Complex>::
832 0 : interpolate(data,flag,interpVisFreq_p,visFreq,flipdata,flipflag,freqInterpMethod_p, False, False);
833 0 : flipdata.resize();
834 0 : swapyz(flipdata,data);
835 0 : data.resize();
836 0 : data.reference(flipdata);
837 0 : flipflag.resize();
838 0 : swapyz(flipflag,flag);
839 0 : flag.resize();
840 0 : flag.reference(flipflag);
841 : // Note : 'flag' will get augmented with the flags coming out of weight interpolation
842 0 : }
843 : else
844 : { // get the flag array to the correct shape.
845 : // This will get filled at the end of weight-interpolation.
846 0 : flag.resize(vb.nCorrelations(), interpVisFreq_p.nelements(), vb.nRows());
847 0 : flag.set(false);
848 : }
849 : // Now, interpolate the weights also.
850 : // (1) Read in the flags from the vb ( setSpectralFlags -> modflagCube )
851 : // (2) Collapse the flags along the polarization dimension to match shape of weight.
852 : //If BriggsWeightor is used weight is already interpolated so we can bypass this
853 0 : InterpolateArray1D<casacore::Double,casacore::Complex>::InterpolationMethod weightinterp=freqInterpMethod_p;
854 :
855 0 : if(!briggsWeightor_p.null()){
856 0 : weightinterp= InterpolateArray1D<casacore::Double,casacore::Complex>::nearestNeighbour;
857 : }
858 : //InterpolateArray1D<casacore::Double,casacore::Complex>::InterpolationMethod weightinterp=InterpolateArray1D<casacore::Double,casacore::Complex>::nearestNeighbour;
859 0 : Matrix<Bool> chanflag(wt.shape());
860 0 : AlwaysAssert( chanflag.shape()[0]==modflagCube.shape()[1], AipsError);
861 0 : AlwaysAssert( chanflag.shape()[1]==modflagCube.shape()[2], AipsError);
862 0 : chanflag=false;
863 0 : for(uInt pol=0;pol<modflagCube.shape()[0];pol++)
864 0 : chanflag = chanflag | modflagCube.yzPlane(pol);
865 :
866 : // (3) Interpolate the weights.
867 : // Input flags are the collapsed vb flags : 'chanflag'
868 : // Output flags are in tempoutputflag
869 : // - contains info about channels that couldn't be interpolated.
870 0 : Matrix<Float> flipweight;
871 0 : flipweight=transpose(wt);
872 0 : Matrix<Bool> flipchanflag;
873 0 : flipchanflag=transpose(chanflag);
874 0 : Matrix<Bool> tempoutputflag;
875 : InterpolateArray1D<Double,Float>::
876 0 : interpolate(weight,tempoutputflag, interpVisFreq_p, visFreq,flipweight,flipchanflag,weightinterp, False, False);
877 0 : flipweight.resize();
878 0 : flipweight=transpose(weight);
879 0 : weight.resize();
880 0 : weight.reference(flipweight);
881 0 : flipchanflag.resize();
882 0 : flipchanflag=transpose(tempoutputflag);
883 0 : tempoutputflag.resize();
884 0 : tempoutputflag.reference(flipchanflag);
885 :
886 : // (4) Now, fill these flags back into the flag cube
887 : // so that they get USED while gridding the PSF (and data)
888 : // Taking the OR of the flags that came out of data-interpolation
889 : // and weight-interpolation, in case they're different.
890 : // Expanding flags across polarization. This will destroy any
891 : // pol-dependent flags for imaging, but msvis::VisImagingWeight
892 : // uses the OR of flags across polarization anyway
893 : // so we don't lose anything.
894 :
895 0 : AlwaysAssert( tempoutputflag.shape()[0]==flag.shape()[1], AipsError);
896 0 : AlwaysAssert( tempoutputflag.shape()[1]==flag.shape()[2], AipsError);
897 0 : for(uInt pol=0;pol<flag.shape()[0];pol++)
898 0 : flag.yzPlane(pol) = tempoutputflag | flag.yzPlane(pol);
899 :
900 : // Fill the output array of image-channel flags.
901 0 : flags.resize(flag.shape());
902 0 : flags=0;
903 0 : flags(flag)=true;
904 :
905 0 : return true;
906 6120 : }
907 :
908 0 : void FTMachine::getInterpolateArrays(const vi::VisBuffer2& vb,
909 : Cube<Complex>& data, Cube<Int>& flags){
910 :
911 0 : Vector<Double> visFreq(vb.getFrequencies(0).nelements());
912 :
913 : //if(doConversion_p[vb.spectralWindows()[0]]){
914 0 : if(freqFrameValid_p){
915 0 : convertArray(visFreq, lsrFreq_p);
916 : }
917 : else{
918 0 : convertArray(visFreq, vb.getFrequencies(0));
919 : }
920 :
921 :
922 :
923 0 : if((imageFreq_p.nelements()==1) || (freqInterpMethod_p== InterpolateArray1D<Double, Complex>::nearestNeighbour)|| (vb.nChannels()==1) ){
924 0 : Cube<Bool> modflagCube;
925 0 : setSpectralFlag(vb,modflagCube);
926 :
927 0 : data.reference(vb.visCubeModel());
928 : //flags.resize(vb.flagCube().shape());
929 0 : flags.resize(modflagCube.shape());
930 0 : flags=0;
931 : //flags(vb.flagCube())=true;
932 0 : flags(modflagCube)=true;
933 0 : interpVisFreq_p.resize();
934 0 : interpVisFreq_p=visFreq;
935 0 : return;
936 0 : }
937 :
938 0 : data.resize(vb.nCorrelations(), imageFreq_p.nelements(), vb.nRows());
939 0 : flags.resize(vb.nCorrelations(), imageFreq_p.nelements(), vb.nRows());
940 0 : data.set(Complex(0.0,0.0));
941 0 : flags.set(0);
942 : //no need to degrid channels that does map over this vb
943 0 : Int maxchan=max(chanMap);
944 0 : for (uInt k =0 ; k < chanMap.nelements() ; ++k){
945 0 : if(chanMap(k)==-1)
946 0 : chanMap(k)=maxchan;
947 : }
948 0 : Int minchan=min(chanMap);
949 0 : if(minchan==maxchan)
950 0 : minchan=-1;
951 :
952 :
953 0 : for(Int k = 0; k < minchan; ++k)
954 0 : flags.xzPlane(k).set(1);
955 :
956 0 : for(uInt k = maxchan + 1; k < imageFreq_p.nelements(); ++k)
957 0 : flags.xzPlane(k).set(1);
958 :
959 0 : interpVisFreq_p.resize(imageFreq_p.nelements());
960 0 : convertArray(interpVisFreq_p, imageFreq_p);
961 0 : chanMap.resize(imageFreq_p.nelements());
962 0 : indgen(chanMap);
963 0 : }
964 :
965 0 : Bool FTMachine::interpolateFrequencyFromgrid(vi::VisBuffer2& vb,
966 : Cube<Complex>& data,
967 : FTMachine::Type type){
968 :
969 : Cube<Complex> *origdata;
970 0 : Vector<Double> visFreq(vb.getFrequencies(0).nelements());
971 :
972 : //if(doConversion_p[vb.spectralWindows()[0]]){
973 0 : if(freqFrameValid_p){
974 0 : convertArray(visFreq, lsrFreq_p);
975 : }
976 : else{
977 0 : convertArray(visFreq, vb.getFrequencies(0));
978 : }
979 :
980 0 : if(type==FTMachine::MODEL){
981 0 : origdata=const_cast <Cube<Complex>* > (&(vb.visCubeModel()));
982 : }
983 0 : else if(type==FTMachine::CORRECTED){
984 0 : origdata=const_cast<Cube<Complex>* >(&(vb.visCubeCorrected()));
985 : }
986 : else{
987 0 : origdata=const_cast<Cube<Complex>* >(&(vb.visCube()));
988 : }
989 : //
990 : // If visibility data (vb) has only one channel, or the image cube
991 : // has only one channel, resort to nearestNeighbour interpolation.
992 : // Honour user selection of nearestNeighbour.
993 : //
994 :
995 0 : Double width=fabs(imageFreq_p[1]-imageFreq_p[0])/fabs(visFreq[1]-visFreq[0]);
996 :
997 0 : if((imageFreq_p.nelements()==1) ||
998 0 : (vb.nChannels()==1) ||
999 0 : (freqInterpMethod_p== InterpolateArray1D<Double, Complex>::nearestNeighbour) ){
1000 0 : interpVisFreq_p=visFreq;
1001 : //cerr << "INTERPFROMGRID " << interpVisFreq_p << " vb.nchan " << vb.nChannels() << endl;
1002 0 : origdata->reference(data);
1003 0 : interpVisFreq_p=visFreq;
1004 0 : return false;
1005 : }
1006 :
1007 : //Need to get new interpolate functions that interpolate explicitly on the 2nd axis
1008 : //2 swap of axes needed
1009 0 : Cube<Complex> flipgrid;
1010 0 : flipgrid.resize();
1011 0 : swapyz(flipgrid,data);
1012 0 : Vector<Double> newImFreq;
1013 0 : newImFreq=imageFreq_p;
1014 :
1015 : //cerr << "width " << width << endl;
1016 : /* if(((width >2.0) && (freqInterpMethod_p==InterpolateArray1D<Double, Complex>::linear)) ||
1017 : ((width >4.0) && (freqInterpMethod_p !=InterpolateArray1D<Double, Complex>::linear))){*/
1018 0 : if(width > 1.0){
1019 0 : Int newNchan=Int(std::floor(width))*imageFreq_p.nelements();
1020 0 : newImFreq.resize(newNchan);
1021 0 : Double newIncr= (imageFreq_p[1]-imageFreq_p[0])/std::floor(width);
1022 0 : Double newStart=imageFreq_p[0]-(imageFreq_p[1]-imageFreq_p[0])/2.0+newIncr/2.0;
1023 0 : Cube<Complex> newflipgrid(flipgrid.shape()[0], flipgrid.shape()[1], newNchan);
1024 :
1025 0 : for (Int k=0; k < newNchan; ++k){
1026 0 : newImFreq[k]=newStart+k*newIncr;
1027 0 : Int oldchan=k/Int(std::floor(width));
1028 0 : newflipgrid.xyPlane(k)=flipgrid.xyPlane(oldchan);
1029 :
1030 : }
1031 : //cerr << std::setprecision(12) << "newfreq " << newImFreq << endl;
1032 : //cerr << "oldfreq " << imageFreq_p << endl;
1033 : //InterpolateArray1D<Double,Complex>::
1034 : //interpolate(newflipgrid,newImFreq, imageFreq_p, flipgrid, InterpolateArray1D<Double, Complex>::nearestNeighbour);
1035 0 : flipgrid.resize();
1036 0 : flipgrid.reference(newflipgrid);
1037 :
1038 0 : }
1039 0 : Cube<Complex> flipdata((origdata->shape())(0),(origdata->shape())(2),
1040 0 : (origdata->shape())(1)) ;
1041 0 : flipdata.set(Complex(0.0));
1042 :
1043 : ///TESTOO
1044 : //Cube<Bool> inflag(flipgrid.shape());
1045 : //inflag.set(False);
1046 : //Cube<Bool> outflag(flipdata.shape());
1047 : //InterpolateArray1D<Double,Complex>::
1048 : // interpolate(flipdata,outflag,visFreq,newImFreq,flipgrid,inflag,freqInterpMethod_p, False, True);
1049 :
1050 : //cerr << "OUTFLAG " << anyEQ(True, outflag) << " chanmap " << chanMap << endl;
1051 : /////End TESTOO
1052 : InterpolateArray1D<Double,Complex>::
1053 0 : interpolate(flipdata,visFreq, newImFreq, flipgrid,freqInterpMethod_p);
1054 :
1055 :
1056 :
1057 0 : Cube<Bool> copyOfFlag;
1058 : //Vector<Int> mychanmap=multiChanMap_p[vb.spectralWindows()[0]];
1059 0 : matchChannel(vb);
1060 0 : copyOfFlag.assign(vb.flagCube());
1061 0 : for (uInt k=0; k< chanMap.nelements(); ++ k)
1062 0 : if(chanMap(k) == -1)
1063 0 : copyOfFlag.xzPlane(k).set(true);
1064 0 : flipgrid.resize();
1065 0 : swapyz(flipgrid, copyOfFlag, flipdata);
1066 : //swapyz(flipgrid,flipdata);
1067 0 : vb.setVisCubeModel(flipgrid);
1068 :
1069 0 : return true;
1070 0 : }
1071 :
1072 :
1073 0 : void FTMachine::girarUVW(Matrix<Double>& uvw, Vector<Double>& dphase,
1074 : const vi::VisBuffer2& vb)
1075 : {
1076 :
1077 :
1078 :
1079 : //the uvw rotation is done for common tangent reprojection or if the
1080 : //image center is different from the phasecenter
1081 : // UVrotation is false only if field never changes
1082 0 : if(lastMSId_p != vb.msId())
1083 0 : romscol_p=new MSColumns(vb.ms());
1084 0 : if((vb.fieldId()(0)!=lastFieldId_p) || (vb.msId()!=lastMSId_p)){
1085 0 : doUVWRotation_p=true;
1086 : }
1087 : else{
1088 : //if above failed it still can be changing if polynome phasecenter or ephem
1089 :
1090 0 : if( (vb.subtableColumns().field().numPoly()(lastFieldId_p) >0) || (! (vb.subtableColumns().field().ephemerisId().isNull()) && (vb.subtableColumns().field().ephemerisId()(lastFieldId_p) > -1)))
1091 0 : doUVWRotation_p=True;
1092 : }
1093 0 : if(doUVWRotation_p || fixMovingSource_p){
1094 :
1095 0 : mFrame_p.epoch() != 0 ?
1096 0 : mFrame_p.resetEpoch(MEpoch(Quantity(vb.time()(0), "s"))):
1097 0 : mFrame_p.set(mLocation_p, MEpoch(Quantity(vb.time()(0), "s"), (romscol_p->timeMeas())(0).getRef()));
1098 : MDirection::Types outType;
1099 0 : MDirection::getType(outType, mImage_p.getRefString());
1100 0 : MDirection phasecenter=MDirection::Convert(vbutil_p->getPhaseCenter(vb, phaseCenterTime_p), MDirection::Ref(outType, mFrame_p))();
1101 0 : MDirection inFieldPhaseCenter=phasecenter;
1102 :
1103 0 : if(fixMovingSource_p){
1104 :
1105 : //First convert to HA-DEC or AZEL for parallax correction
1106 0 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
1107 0 : MDirection tmphadec;
1108 0 : if(upcase(movingDir_p.getRefString()).contains("APP")){
1109 0 : tmphadec=MDirection::Convert((vbutil_p->getEphemDir(vb, phaseCenterTime_p)), outref1)();
1110 : }
1111 : else{
1112 0 : tmphadec=MDirection::Convert(movingDir_p, outref1)();
1113 : }
1114 0 : MDirection::Ref outref(mImage_p.getRef().getType(), mFrame_p);
1115 0 : MDirection sourcenow=MDirection::Convert(tmphadec, outref)();
1116 : //cerr << "Rotating to fixed moving source " << MVDirection(phasecenter.getAngle()-firstMovingDir_p.getAngle()+sourcenow.getAngle()) << endl;
1117 : //phasecenter.set(MVDirection(phasecenter.getAngle()+firstMovingDir_p.getAngle()-sourcenow.getAngle()));
1118 0 : movingDirShift_p=MVDirection(sourcenow.getAngle()-firstMovingDir_p.getAngle());
1119 : // cerr << "shift " << movingDirShift_p.getAngle() << endl;
1120 0 : inFieldPhaseCenter.shift(movingDirShift_p, False);
1121 0 : }
1122 :
1123 :
1124 : // Set up the UVWMachine only if the field id has changed. If
1125 : // the tangent plane is specified then we need a UVWMachine that
1126 : // will reproject to that plane iso the image plane
1127 0 : if(doUVWRotation_p || fixMovingSource_p) {
1128 :
1129 0 : String observatory=(vb.subtableColumns().observation()).telescopeName()(0);
1130 0 : if(uvwMachine_p) delete uvwMachine_p; uvwMachine_p=0;
1131 0 : if(observatory.contains("ATCA") || observatory.contains("WSRT")){
1132 : //Tangent specified is being wrongly used...it should be for a
1133 : //Use the safest way for now.
1134 0 : uvwMachine_p=new UVWMachine(inFieldPhaseCenter, vbutil_p->getPhaseCenter(vb, phaseCenterTime_p), mFrame_p,
1135 0 : true, false);
1136 0 : phaseShifter_p=new UVWMachine(mImage_p, phasecenter, mFrame_p,
1137 0 : true, false);
1138 : }
1139 : else{
1140 0 : uvwMachine_p=new UVWMachine(inFieldPhaseCenter, vbutil_p->getPhaseCenter(vb, phaseCenterTime_p), mFrame_p,
1141 0 : false, false);
1142 0 : phaseShifter_p=new UVWMachine(mImage_p, phasecenter, mFrame_p,
1143 0 : false, false);
1144 : }
1145 0 : }
1146 :
1147 0 : lastFieldId_p=vb.fieldId()(0);
1148 0 : lastMSId_p=vb.msId();
1149 :
1150 :
1151 0 : AlwaysAssert(uvwMachine_p, AipsError);
1152 :
1153 : // Always force a recalculation
1154 0 : uvwMachine_p->reCalculate();
1155 0 : phaseShifter_p->reCalculate();
1156 :
1157 : // Now do the conversions
1158 0 : uInt nrows=dphase.nelements();
1159 0 : Vector<Double> thisRow(3);
1160 0 : thisRow=0.0;
1161 : //CoordinateSystem csys=image->coordinates();
1162 : //DirectionCoordinate dc=csys.directionCoordinate(0);
1163 : //Vector<Double> thePix(2);
1164 : //dc.toPixel(thePix, phasecenter);
1165 : //cerr << "field id " << vb.fieldId() << " the Pix " << thePix << endl;
1166 : //Vector<Float> scale(2);
1167 : //scale(0)=dc.increment()(0);
1168 : //scale(1)=dc.increment()(1);
1169 0 : for (uInt irow=0; irow<nrows;++irow) {
1170 0 : thisRow.assign(uvw.column(irow));
1171 : //cerr << " uvw " << thisRow ;
1172 : // This is for frame change
1173 0 : uvwMachine_p->convertUVW(dphase(irow), thisRow);
1174 : // This is for correlator phase center change
1175 0 : MVPosition rotphase=phaseShifter_p->rotationPhase() ;
1176 : //cerr << " rotPhase " << rotphase << " oldphase "<< rotphase*(uvw.column(irow)) << " newphase " << (rotphase)*thisRow ;
1177 : // cerr << " phase " << dphase(irow) << " new uvw " << uvw.column(irow);
1178 : //dphase(irow)+= (thePix(0)-nx/2.0)*thisRow(0)*scale(0)+(thePix(1)-ny/2.0)*thisRow(1)*scale(1);
1179 : //Double pixphase=(thePix(0)-nx/2.0)*uvw.column(irow)(0)*scale(0)+(thePix(1)-ny/2.0)*uvw.column(irow)(1)*scale(1);
1180 : //Double pixphase2=(thePix(0)-nx/2.0)*thisRow(0)*scale(0)+(thePix(1)-ny/2.0)*thisRow(1)*scale(1);
1181 : //cerr << " pixphase " << pixphase << " pixphase2 " << pixphase2<< endl;
1182 : //dphase(irow)=pixphase;
1183 0 : RotMatrix rotMat=phaseShifter_p->rotationUVW();
1184 0 : uvw.column(irow)(0)=thisRow(0)*rotMat(0,0)+thisRow(1)*rotMat(1,0);
1185 0 : uvw.column(irow)(1)=thisRow(1)*rotMat(1,1)+thisRow(0)*rotMat(0,1);
1186 0 : uvw.column(irow)(2)=thisRow(0)*rotMat(0,2)+thisRow(1)*rotMat(1,2)+thisRow(2)*rotMat(2,2);
1187 0 : dphase(irow)+= rotphase(0)*uvw.column(irow)(0)+rotphase(1)*uvw.column(irow)(1);
1188 0 : }
1189 :
1190 :
1191 0 : }
1192 0 : }
1193 :
1194 :
1195 6120 : void FTMachine::rotateUVW(Matrix<Double>& uvw, Vector<Double>& dphase,
1196 : const vi::VisBuffer2& vb)
1197 : {
1198 :
1199 6120 : if(lastMSId_p != vb.msId())
1200 0 : romscol_p=new MSColumns(vb.ms());
1201 : //the uvw rotation is done for common tangent reprojection or if the
1202 : //image center is different from the phasecenter
1203 : // UVrotation is false only if field never changes
1204 :
1205 6120 : if((vb.fieldId()(0)!=lastFieldId_p) || (vb.msId()!=lastMSId_p)){
1206 34 : doUVWRotation_p=true;
1207 :
1208 : }
1209 : else{
1210 : //if above failed it still can be changing if polynome phasecenter or ephem
1211 6086 : if( (vb.subtableColumns().field().numPoly()(lastFieldId_p) >0) || (! (vb.subtableColumns().field().ephemerisId().isNull()) &&(vb.subtableColumns().field().ephemerisId()(lastFieldId_p) > -1)))
1212 0 : doUVWRotation_p=True;
1213 :
1214 : }
1215 6120 : if(doUVWRotation_p || tangentSpecified_p || fixMovingSource_p){
1216 6120 : ok();
1217 :
1218 6120 : mFrame_p.epoch() != 0 ?
1219 12240 : mFrame_p.resetEpoch(MEpoch(Quantity(vb.time()(0), "s"))):
1220 :
1221 6120 : mFrame_p.set(mLocation_p, MEpoch(Quantity(vb.time()(0), "s"), (romscol_p->timeMeas())(0).getRef()));
1222 :
1223 6120 : MDirection phasecenter=mImage_p;
1224 6120 : if(fixMovingSource_p){
1225 :
1226 : //First convert to HA-DEC or AZEL for parallax correction
1227 0 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
1228 0 : MDirection tmphadec;
1229 0 : if(upcase(movingDir_p.getRefString()).contains("APP")){
1230 0 : tmphadec=MDirection::Convert((vbutil_p->getEphemDir(vb, phaseCenterTime_p)), outref1)();
1231 : }
1232 : else{
1233 0 : tmphadec=MDirection::Convert(movingDir_p, outref1)();
1234 : }
1235 : ////TESTOO waiting for CAS-11060
1236 : //MDirection tmphadec=MDirection::Convert((vbutil_p->getPhaseCenter(vb, phaseCenterTime_p)), outref1)();
1237 : /////////
1238 0 : MDirection::Ref outref(mImage_p.getRef().getType(), mFrame_p);
1239 0 : MDirection sourcenow=MDirection::Convert(tmphadec, outref)();
1240 : //cerr << "Rotating to fixed moving source " << MVDirection(phasecenter.getAngle()-firstMovingDir_p.getAngle()+sourcenow.getAngle()) << endl;
1241 : //phasecenter.set(MVDirection(phasecenter.getAngle()+firstMovingDir_p.getAngle()-sourcenow.getAngle()));
1242 0 : movingDirShift_p=MVDirection(sourcenow.getAngle()-firstMovingDir_p.getAngle());
1243 0 : phasecenter.shift(movingDirShift_p, False);
1244 : //cerr << sourcenow.toString() <<" "<<(vbutil_p->getPhaseCenter(vb, phaseCenterTime_p)).toString() << " difference " << firstMovingDir_p.getAngle() - sourcenow.getAngle() << endl;
1245 0 : }
1246 :
1247 :
1248 : // Set up the UVWMachine only if the field id has changed. If
1249 : // the tangent plane is specified then we need a UVWMachine that
1250 : // will reproject to that plane iso the image plane
1251 6120 : if(doUVWRotation_p || fixMovingSource_p) {
1252 :
1253 6120 : String observatory=(vb.subtableColumns().observation()).telescopeName()(0);
1254 6120 : if(uvwMachine_p) delete uvwMachine_p; uvwMachine_p=0;
1255 6120 : if(observatory.contains("ATCA") || observatory.contains("WSRT")){
1256 : //Tangent specified is being wrongly used...it should be for a
1257 : //Use the safest way for now.
1258 0 : uvwMachine_p=new UVWMachine(phasecenter, vbutil_p->getPhaseCenter(vb, phaseCenterTime_p), mFrame_p,
1259 0 : true, false);
1260 : }
1261 : else{
1262 12240 : uvwMachine_p=new UVWMachine(phasecenter, vbutil_p->getPhaseCenter(vb, phaseCenterTime_p), mFrame_p,
1263 6120 : false,tangentSpecified_p);
1264 : }
1265 6120 : }
1266 :
1267 6120 : lastFieldId_p=vb.fieldId()(0);
1268 6120 : lastMSId_p=vb.msId();
1269 :
1270 :
1271 6120 : AlwaysAssert(uvwMachine_p, AipsError);
1272 :
1273 : // Always force a recalculation
1274 6120 : uvwMachine_p->reCalculate();
1275 :
1276 : // Now do the conversions
1277 6120 : uInt nrows=dphase.nelements();
1278 6120 : Vector<Double> thisRow(3);
1279 6120 : thisRow=0.0;
1280 : uInt irow;
1281 : //#pragma omp parallel default(shared) private(irow,thisRow)
1282 : {
1283 : //#pragma omp for
1284 740520 : for (irow=0; irow<nrows;++irow) {
1285 734400 : thisRow.reference(uvw.column(irow));
1286 734400 : convUVW(dphase(irow), thisRow);
1287 : }
1288 :
1289 : }//end pragma
1290 6120 : }
1291 6120 : }
1292 734400 : void FTMachine::convUVW(Double& dphase, Vector<Double>& thisrow){
1293 : //for (uInt i=0;i<3;i++) thisRow(i)=uvw(i,row);
1294 734400 : uvwMachine_p->convertUVW(dphase, thisrow);
1295 : //for (uint i=0;i<3;i++) uvw(i,row)=thisRow(i);
1296 :
1297 734400 : }
1298 :
1299 :
1300 0 : void FTMachine::locateuvw(const Double*& uvw, const Double*& dphase,
1301 : const Double*& freq, const Int& nvchan,
1302 : const Double*& scale, const Double*& offset, const Int& sampling, Int*& loc, Int*& off, Complex*& phasor, const Int& row, const bool& doW){
1303 :
1304 0 : Int rowoff=row*nvchan;
1305 : Double phase;
1306 : Double pos;
1307 0 : Int nel= doW ? 3 : 2;
1308 0 : for (Int f=0; f<nvchan; ++f){
1309 0 : for (Int k=0; k <2; ++k){
1310 0 : pos=(scale[k])*uvw[3*row+k]*(freq[f])/C::c+((offset[k])+1.0);
1311 0 : loc[(rowoff+f)*nel+k]=std::lround(pos);
1312 0 : off[(rowoff+f)*nel+k]=std::lround((Double(loc[(rowoff+f)*nel+k])-pos)*Double(sampling));
1313 : //off[(rowoff+f)*2+k]=(loc[(rowoff+f)*2+k]-pos(k))*sampling;
1314 : }
1315 0 : phase=-Double(2.0)*M_PI*dphase[row]*(freq[f])/C::c;
1316 0 : phasor[rowoff+f]=Complex(cos(phase), sin(phase));
1317 :
1318 : ///This is for W-Projection
1319 0 : if(doW){
1320 0 : pos=sqrt(fabs(scale[2]*uvw[3*row+2]*(freq[f])/C::c))+offset[2]+1.0;
1321 0 : loc[(rowoff+f)*nel+2]=std::lround(pos);
1322 0 : off[(rowoff+f)*nel+2]=0;
1323 : }
1324 : }
1325 :
1326 :
1327 :
1328 :
1329 0 : }
1330 :
1331 0 : void FTMachine::setnumthreads(Int num){
1332 0 : numthreads_p=num;
1333 0 : }
1334 0 : Int FTMachine::getnumthreads(){
1335 0 : return numthreads_p;
1336 : }
1337 :
1338 : //
1339 : // Refocus the array on a point at finite distance
1340 : //
1341 : // Refocus the array on a point at finite distance
1342 : //
1343 6120 : void FTMachine::refocus(Matrix<Double>& uvw, const Vector<Int>& ant1,
1344 : const Vector<Int>& ant2,
1345 : Vector<Double>& dphase, const vi::VisBuffer2& vb)
1346 : {
1347 :
1348 6120 : ok();
1349 :
1350 6120 : if(abs(distance_p)>0.0) {
1351 :
1352 0 : nAntenna_p=max(vb.antenna2())+1;
1353 :
1354 : // Positions of antennas
1355 0 : Matrix<Double> antPos(3,nAntenna_p);
1356 0 : antPos=0.0;
1357 0 : Vector<Int> nAntPos(nAntenna_p);
1358 0 : nAntPos=0;
1359 :
1360 0 : uInt aref = min(ant1);
1361 :
1362 : // Now find the antenna locations: for this we just reference to a common
1363 : // point. We ignore the time variation within this buffer.
1364 0 : uInt nrows=dphase.nelements();
1365 0 : for (uInt row=0;row<nrows;row++) {
1366 0 : uInt a1=ant1(row);
1367 0 : uInt a2=ant2(row);
1368 0 : for(uInt dim=0;dim<3;dim++) {
1369 0 : antPos(dim, a1)+=uvw(dim, row);
1370 0 : antPos(dim, a2)-=uvw(dim, row);
1371 : }
1372 0 : nAntPos(a1)+=1;
1373 0 : nAntPos(a2)+=1;
1374 : }
1375 :
1376 : // Now remove the reference location
1377 0 : Vector<Double> center(3);
1378 0 : for(uInt dim=0;dim<3;dim++) {
1379 0 : center(dim) = antPos(dim,aref)/nAntPos(aref);
1380 : }
1381 :
1382 : // Now normalize
1383 0 : for (uInt ant=0; ant<nAntenna_p; ant++) {
1384 0 : if(nAntPos(ant)>0) {
1385 0 : for(uInt dim=0;dim<3;dim++) {
1386 0 : antPos(dim,ant)/=nAntPos(ant);
1387 0 : antPos(dim,ant)-=center(dim);
1388 : }
1389 : }
1390 : }
1391 :
1392 : // Now calculate the offset needed to focus the array,
1393 : // including the w term. This must have the correct asymptotic
1394 : // form so that at infinity no net change occurs
1395 0 : for (uInt row=0;row<nrows;row++) {
1396 0 : uInt a1=ant1(row);
1397 0 : uInt a2=ant2(row);
1398 :
1399 0 : Double d1=distance_p*distance_p-2*distance_p*antPos(2,a1);
1400 0 : Double d2=distance_p*distance_p-2*distance_p*antPos(2,a2);
1401 0 : for(uInt dim=0;dim<3;dim++) {
1402 0 : d1+=antPos(dim,a1)*antPos(dim,a1);
1403 0 : d2+=antPos(dim,a2)*antPos(dim,a2);
1404 : }
1405 0 : d1=sqrt(d1);
1406 0 : d2=sqrt(d2);
1407 0 : for(uInt dim=0;dim<2;dim++) {
1408 0 : dphase(row)-=(antPos(dim,a1)*antPos(dim,a1)-antPos(dim,a2)*antPos(dim,a2))/(2*distance_p);
1409 : }
1410 0 : uvw(0,row)=distance_p*(antPos(0,a1)/d1-antPos(0,a2)/d2);
1411 0 : uvw(1,row)=distance_p*(antPos(1,a1)/d1-antPos(1,a2)/d2);
1412 0 : uvw(2,row)=distance_p*(antPos(2,a1)/d1-antPos(2,a2)/d2)+dphase(row);
1413 : }
1414 0 : }
1415 6120 : }
1416 :
1417 0 : void FTMachine::ok() {
1418 0 : AlwaysAssert(image, AipsError);
1419 0 : AlwaysAssert(uvwMachine_p, AipsError);
1420 0 : }
1421 :
1422 0 : Bool FTMachine::toRecord(String& error, RecordInterface& outRecord,
1423 : Bool withImage, const String diskimage) {
1424 : // Save the FTMachine to a Record
1425 : //
1426 0 : outRecord.define("name", this->name());
1427 0 : if(withImage){
1428 0 : if(image==nullptr)
1429 0 : throw(AipsError("Programmer error: saving to record without proper initialization"));
1430 0 : CoordinateSystem cs=image->coordinates();
1431 0 : DirectionCoordinate dircoord=cs.directionCoordinate(cs.findCoordinate(Coordinate::DIRECTION));
1432 0 : dircoord.setReferenceValue(mImage_p.getAngle().getValue());
1433 0 : if(diskimage != ""){
1434 : try{
1435 0 : PagedImage<Complex> imCopy(TiledShape(toVis_p ? griddedData.shape(): image->shape()), image->coordinates(), diskimage);
1436 0 : toVis_p ? imCopy.put(griddedData) : imCopy.copyData(*image);
1437 0 : ImageUtilities::copyMiscellaneous(imCopy, *image);
1438 0 : Vector<Double> pixcen(2);
1439 0 : pixcen(0)=Double(imCopy.shape()(0)/2); pixcen(1)=Double(imCopy.shape()(1)/2);
1440 0 : dircoord.setReferencePixel(pixcen);
1441 0 : cs.replaceCoordinate(dircoord, cs.findCoordinate(Coordinate::DIRECTION));
1442 0 : imCopy.setCoordinateInfo(cs);
1443 0 : }
1444 0 : catch(...){
1445 0 : throw(AipsError(String("Failed to save model image "+diskimage+String(" to disk"))));
1446 0 : }
1447 0 : outRecord.define("diskimage", diskimage);
1448 :
1449 : }
1450 : else{
1451 0 : Record imrec;
1452 0 : Vector<Double> pixcen(2);
1453 0 : pixcen(0)=Double(griddedData.shape()(0)/2); pixcen(1)=Double(griddedData.shape()(1)/2);
1454 0 : dircoord.setReferencePixel(pixcen);
1455 0 : cs.replaceCoordinate(dircoord, cs.findCoordinate(Coordinate::DIRECTION));
1456 0 : TempImage<Complex> imCopy(griddedData.shape(), cs);
1457 0 : imCopy.put(griddedData) ;
1458 0 : ImageUtilities::copyMiscellaneous(imCopy, *image);
1459 0 : if(imCopy.toRecord(error, imrec))
1460 0 : outRecord.defineRecord("image", imrec);
1461 0 : }
1462 0 : }
1463 0 : outRecord.define("nantenna", nAntenna_p);
1464 0 : outRecord.define("distance", distance_p);
1465 0 : outRecord.define("lastfieldid", lastFieldId_p);
1466 0 : outRecord.define("lastmsid", lastMSId_p);
1467 0 : outRecord.define("tangentspecified", tangentSpecified_p);
1468 0 : saveMeasure(outRecord, String("mtangent_rec"), error, mTangent_p);
1469 0 : saveMeasure(outRecord, "mimage_rec", error, mImage_p);
1470 : //mFrame_p not necessary to save as it is generated from mLocation_p
1471 0 : outRecord.define("nx", nx);
1472 0 : outRecord.define("ny", ny);
1473 0 : outRecord.define("npol", npol);
1474 0 : outRecord.define("nchan", nchan);
1475 0 : outRecord.define("nvischan", nvischan);
1476 0 : outRecord.define("nvispol", nvispol);
1477 : //no need to save uvwMachine_p
1478 0 : outRecord.define("douvwrotation", doUVWRotation_p);
1479 0 : outRecord.define("freqinterpmethod", static_cast<Int>(freqInterpMethod_p));
1480 0 : outRecord.define("spwchanselflag", spwChanSelFlag_p);
1481 0 : outRecord.define("freqframevalid", freqFrameValid_p);
1482 : //outRecord.define("selectedspw", selectedSpw_p);
1483 0 : outRecord.define("imagefreq", imageFreq_p);
1484 0 : outRecord.define("lsrfreq", lsrFreq_p);
1485 0 : outRecord.define("interpvisfreq", interpVisFreq_p);
1486 0 : Record multichmaprec;
1487 : //for (uInt k=0; k < multiChanMap_p.nelements(); ++k)
1488 : // multichmaprec.define(k, multiChanMap_p[k]);
1489 : //outRecord.defineRecord("multichanmaprec", multichmaprec);
1490 0 : outRecord.define("chanmap", chanMap);
1491 0 : outRecord.define("polmap", polMap);
1492 0 : outRecord.define("nvischanmulti", nVisChan_p);
1493 : //save moving source related variables
1494 0 : storeMovingSourceState(error, outRecord);
1495 : //outRecord.define("doconversion", doConversion_p);
1496 0 : outRecord.define("pointingdircol", pointingDirCol_p);
1497 0 : outRecord.define("usedoublegrid", useDoubleGrid_p);
1498 0 : outRecord.define("cfstokes", cfStokes_p);
1499 0 : outRecord.define("polinuse", polInUse_p);
1500 0 : outRecord.define("tovis", toVis_p);
1501 0 : outRecord.define("sumweight", sumWeight);
1502 0 : outRecord.define("numthreads", numthreads_p);
1503 0 : outRecord.define("phasecentertime", phaseCenterTime_p);
1504 : //Need to serialized sj_p...the user has to set the sj_p after recovering from record
1505 0 : return true;
1506 0 : };
1507 :
1508 0 : Bool FTMachine::saveMeasure(RecordInterface& rec, const String& name, String& err, const Measure& meas){
1509 0 : Record tmprec;
1510 0 : MeasureHolder mh(meas);
1511 0 : if(mh.toRecord(err, tmprec)){
1512 0 : rec.defineRecord(name, tmprec);
1513 0 : return true;
1514 : }
1515 0 : return false;
1516 0 : }
1517 :
1518 0 : Bool FTMachine::fromRecord(String& error,
1519 : const RecordInterface& inRecord
1520 : ) {
1521 : // Restore an FTMachine from a Record
1522 : //
1523 0 : uvwMachine_p=0; //when null it is reconstructed from mImage_p and mFrame_p
1524 : //mFrame_p is not necessary as it is generated in initMaps from mLocation_p
1525 0 : inRecord.get("nx", nx);
1526 0 : inRecord.get("ny", ny);
1527 0 : inRecord.get("npol", npol);
1528 0 : inRecord.get("nchan", nchan);
1529 0 : inRecord.get("nvischan", nvischan);
1530 0 : inRecord.get("nvispol", nvispol);
1531 0 : cmplxImage_p=NULL;
1532 0 : inRecord.get("tovis", toVis_p);
1533 0 : if(inRecord.isDefined("image")){
1534 0 : cmplxImage_p=new TempImage<Complex>();
1535 0 : image=&(*cmplxImage_p);
1536 :
1537 0 : const Record rec=inRecord.asRecord("image");
1538 0 : if(!cmplxImage_p->fromRecord(error, rec))
1539 0 : return false;
1540 :
1541 0 : }
1542 0 : else if(inRecord.isDefined("diskimage")){
1543 0 : String theDiskImage;
1544 0 : inRecord.get("diskimage", theDiskImage);
1545 : try{
1546 0 : File eldir(theDiskImage);
1547 0 : if(! eldir.exists()){
1548 0 : String subPathname[30];
1549 0 : String sep = "/";
1550 0 : uInt nw = split(theDiskImage, subPathname, 20, sep);
1551 0 : String theposs=(subPathname[nw-1]);
1552 0 : Bool isExistant=File(theposs).exists();
1553 0 : if(isExistant)
1554 0 : theDiskImage=theposs;
1555 0 : for (uInt i=nw-2 ; i>0; --i){
1556 0 : theposs=subPathname[i]+"/"+theposs;
1557 0 : File newEldir(theposs);
1558 0 : if(newEldir.exists()){
1559 0 : isExistant=true;
1560 0 : theDiskImage=theposs;
1561 : }
1562 0 : }
1563 0 : if(!isExistant)
1564 0 : throw(AipsError("Could not locate mage"));
1565 0 : }
1566 0 : cmplxImage_p=new PagedImage<Complex> (theDiskImage);
1567 0 : image=&(*cmplxImage_p);
1568 0 : }
1569 0 : catch(...){
1570 0 : throw(AipsError(String("Failure to load ")+theDiskImage+String(" image from disk")));
1571 0 : }
1572 0 : }
1573 0 : if(toVis_p && !cmplxImage_p.null()) {
1574 0 : griddedData.resize(image->shape());
1575 0 : griddedData=image->get();
1576 : }
1577 0 : else if(!toVis_p){
1578 0 : IPosition gridShape(4, nx, ny, npol, nchan);
1579 0 : griddedData.resize(gridShape);
1580 0 : griddedData=Complex(0.0);
1581 0 : }
1582 :
1583 0 : nAntenna_p=inRecord.asuInt("nantenna");
1584 0 : distance_p=inRecord.asDouble("distance");
1585 0 : lastFieldId_p=inRecord.asInt("lastfieldid");
1586 0 : lastMSId_p=inRecord.asInt("lastmsid");
1587 0 : inRecord.get("tangentspecified", tangentSpecified_p);
1588 0 : { const Record rec=inRecord.asRecord("mtangent_rec");
1589 0 : MeasureHolder mh;
1590 0 : if(!mh.fromRecord(error, rec))
1591 0 : return false;
1592 0 : mTangent_p=mh.asMDirection();
1593 0 : }
1594 0 : { const Record rec=inRecord.asRecord("mimage_rec");
1595 0 : MeasureHolder mh;
1596 0 : if(!mh.fromRecord(error, rec))
1597 0 : return false;
1598 0 : mImage_p=mh.asMDirection();
1599 0 : }
1600 :
1601 :
1602 :
1603 0 : inRecord.get("douvwrotation", doUVWRotation_p);
1604 :
1605 : //inRecord.get("spwchanselflag", spwChanSelFlag_p);
1606 : //We won't respect the chanselflag as the vister may have different selections
1607 0 : spwChanSelFlag_p.resize();
1608 0 : inRecord.get("freqframevalid", freqFrameValid_p);
1609 : //inRecord.get("selectedspw", selectedSpw_p);
1610 0 : inRecord.get("imagefreq", imageFreq_p);
1611 0 : inRecord.get("lsrfreq", lsrFreq_p);
1612 0 : inRecord.get("interpvisfreq", interpVisFreq_p);
1613 : //const Record multichmaprec=inRecord.asRecord("multichanmaprec");
1614 : //multiChanMap_p.resize(multichmaprec.nfields(), true, false);
1615 : //for (uInt k=0; k < multichmaprec.nfields(); ++k)
1616 : // multichmaprec.get(k, multiChanMap_p[k]);
1617 0 : inRecord.get("chanmap", chanMap);
1618 0 : inRecord.get("polmap", polMap);
1619 0 : inRecord.get("nvischanmulti", nVisChan_p);
1620 : //inRecord.get("doconversion", doConversion_p);
1621 0 : inRecord.get("pointingdircol", pointingDirCol_p);
1622 :
1623 :
1624 0 : inRecord.get("usedoublegrid", useDoubleGrid_p);
1625 0 : inRecord.get("cfstokes", cfStokes_p);
1626 0 : inRecord.get("polinuse", polInUse_p);
1627 :
1628 :
1629 0 : inRecord.get("sumweight", sumWeight);
1630 0 : if(toVis_p){
1631 0 : freqInterpMethod_p=InterpolateArray1D<Double, Complex>::nearestNeighbour;
1632 : }
1633 : else{
1634 : Int tmpInt;
1635 0 : inRecord.get("freqinterpmethod", tmpInt);
1636 0 : freqInterpMethod_p=static_cast<InterpolateArray1D<Double, Complex >::InterpolationMethod>(tmpInt);
1637 : }
1638 0 : inRecord.get("numthreads", numthreads_p);
1639 0 : inRecord.get("phasecentertime", phaseCenterTime_p);
1640 : ///No need to store this...recalculate thread partion because environment
1641 : ///may have changed.
1642 0 : doneThreadPartition_p=-1;
1643 0 : vbutil_p=nullptr;
1644 0 : briggsWeightor_p=nullptr;
1645 0 : ft_p=FFT2D(true);
1646 0 : if(!recoverMovingSourceState(error, inRecord))
1647 0 : return False;
1648 0 : return true;
1649 : };
1650 0 : Bool FTMachine::storeMovingSourceState(String& error, RecordInterface& outRecord){
1651 :
1652 0 : Bool retval=True;
1653 0 : retval=retval && saveMeasure(outRecord, "mlocation_rec", error, mLocation_p);
1654 0 : spectralCoord_p.save(outRecord, "spectralcoord");
1655 0 : retval=retval && saveMeasure(outRecord, "movingdir_rec", error, movingDir_p);
1656 0 : outRecord.define("fixmovingsource", fixMovingSource_p);
1657 0 : retval=retval && saveMeasure(outRecord, "firstmovingdir_rec", error, firstMovingDir_p);
1658 0 : movingDirShift_p=MVDirection(0.0);
1659 0 : if( mFrame_p.comet()){
1660 0 : String ephemTab=MeasComet(*(mFrame_p.comet())).getTablePath();
1661 0 : outRecord.define("ephemeristable",ephemTab);
1662 0 : }
1663 0 : return retval;
1664 : }
1665 0 : Bool FTMachine::recoverMovingSourceState(String& error, const RecordInterface& inRecord){
1666 0 : Bool retval=True;
1667 0 : inRecord.get("fixmovingsource", fixMovingSource_p);
1668 0 : { const Record rec=inRecord.asRecord("firstmovingdir_rec");
1669 0 : MeasureHolder mh;
1670 0 : if(!mh.fromRecord(error, rec))
1671 0 : return false;
1672 0 : firstMovingDir_p=mh.asMDirection();
1673 0 : }
1674 0 : { const Record rec=inRecord.asRecord("movingdir_rec");
1675 0 : MeasureHolder mh;
1676 0 : if(!mh.fromRecord(error, rec))
1677 0 : return false;
1678 0 : movingDir_p=mh.asMDirection();
1679 0 : }
1680 0 : { const Record rec=inRecord.asRecord("mlocation_rec");
1681 0 : MeasureHolder mh;
1682 0 : if(!mh.fromRecord(error, rec))
1683 0 : return false;
1684 0 : mLocation_p=mh.asMPosition();
1685 0 : }
1686 0 : SpectralCoordinate *tmpSpec=SpectralCoordinate::restore(inRecord, "spectralcoord");
1687 0 : if(tmpSpec){
1688 0 : spectralCoord_p=*tmpSpec;
1689 0 : delete tmpSpec;
1690 : }
1691 0 : visPolMap_p.resize();
1692 0 : if(inRecord.isDefined("ephemeristable")){
1693 0 : String ephemtab;
1694 0 : inRecord.get("ephemeristable", ephemtab);
1695 0 : MeasComet laComet;
1696 0 : if(Table::isReadable(ephemtab, False)){
1697 0 : Table laTable(ephemtab);
1698 0 : Path leSentier(ephemtab);
1699 0 : laComet=MeasComet(laTable, leSentier.absoluteName());
1700 0 : }
1701 : else{
1702 0 : laComet= MeasComet(ephemtab);
1703 : }
1704 0 : if(!mFrame_p.comet())
1705 0 : mFrame_p.set(laComet);
1706 : else
1707 0 : mFrame_p.resetComet(laComet);
1708 0 : }
1709 :
1710 0 : return retval;
1711 : }
1712 :
1713 :
1714 6120 : void FTMachine::getImagingWeight(Matrix<Float>& imwgt, const vi::VisBuffer2& vb){
1715 : //cerr << "BRIGGSweightor " << briggsWeightor_p.null() << " or " << !briggsWeoght_p << endl;
1716 6120 : if(briggsWeightor_p.null()){
1717 6120 : imwgt=vb.imagingWeight();
1718 : }
1719 : else{
1720 0 : briggsWeightor_p->weightUniform(imwgt, vb);
1721 : }
1722 :
1723 6120 : }
1724 : // Make a plain straightforward honest-to-FSM image. This returns
1725 : // a complex image, without conversion to Stokes. The representation
1726 : // is that required for the visibilities.
1727 : //----------------------------------------------------------------------
1728 0 : void FTMachine::makeImage(FTMachine::Type type,
1729 : vi::VisibilityIterator2& vi,
1730 : ImageInterface<Complex>& theImage,
1731 : Matrix<Float>& weight) {
1732 :
1733 :
1734 0 : logIO() << LogOrigin("FTMachine", "makeImage0") << LogIO::NORMAL;
1735 :
1736 : // Loop over all visibilities and pixels
1737 0 : vi::VisBuffer2* vb=vi.getVisBuffer();
1738 :
1739 : // Initialize put (i.e. transform to Sky) for this model
1740 0 : vi.origin();
1741 :
1742 0 : if(vb->polarizationFrame()==MSIter::Linear) {
1743 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
1744 : }
1745 : else {
1746 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
1747 : }
1748 :
1749 0 : initializeToSky(theImage,weight,*vb);
1750 : //This call is a NOP for all weighting schemes except for cube-briggs-perchanweightdensity
1751 0 : initBriggsWeightor(vi);
1752 0 : Bool useCorrected= !(MSColumns(vi.ms()).correctedData().isNull());
1753 0 : if((type==FTMachine::CORRECTED) && (!useCorrected))
1754 0 : type=FTMachine::OBSERVED;
1755 0 : Bool normalize=true;
1756 0 : if(type==FTMachine::COVERAGE)
1757 0 : normalize=false;
1758 :
1759 : // Loop over the visibilities, putting VisBuffers
1760 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
1761 0 : for (vi.origin(); vi.more(); vi.next()) {
1762 :
1763 0 : switch(type) {
1764 0 : case FTMachine::RESIDUAL:
1765 0 : vb->setVisCube(vb->visCubeCorrected());
1766 0 : vb->setVisCube(vb->visCube()-vb->visCubeModel());
1767 0 : put(*vb, -1, false);
1768 0 : break;
1769 0 : case FTMachine::MODEL:
1770 0 : put(*vb, -1, false, FTMachine::MODEL);
1771 0 : break;
1772 0 : case FTMachine::CORRECTED:
1773 0 : put(*vb, -1, false, FTMachine::CORRECTED);
1774 0 : break;
1775 0 : case FTMachine::PSF:
1776 0 : vb->setVisCube(Complex(1.0,0.0));
1777 0 : put(*vb, -1, true, FTMachine::PSF);
1778 0 : break;
1779 0 : case FTMachine::COVERAGE:
1780 0 : vb->setVisCube(Complex(1.0));
1781 0 : put(*vb, -1, true, FTMachine::COVERAGE);
1782 0 : break;
1783 0 : case FTMachine::OBSERVED:
1784 : default:
1785 0 : put(*vb, -1, false, FTMachine::OBSERVED);
1786 0 : break;
1787 : }
1788 : }
1789 : }
1790 0 : finalizeToSky();
1791 : // Normalize by dividing out weights, etc.
1792 0 : getImage(weight, normalize);
1793 0 : }
1794 :
1795 :
1796 :
1797 :
1798 34 : Bool FTMachine::setFrameValidity(Bool validFrame){
1799 :
1800 34 : freqFrameValid_p=validFrame;
1801 34 : return true;
1802 : }
1803 :
1804 :
1805 0 : Vector<Int> FTMachine::channelMap(const vi::VisBuffer2& vb){
1806 0 : matchChannel(vb);
1807 0 : return chanMap;
1808 : }
1809 6154 : Bool FTMachine::matchChannel(const vi::VisBuffer2& vb){
1810 : //Int spw=vb.spectralWindows()[0];
1811 6154 : nvischan = vb.nChannels();
1812 6154 : if(lastMSId_p != vb.msId() || romscol_p.null())
1813 0 : romscol_p = new MSColumns(vb.ms());
1814 : //Try to avoid a bug in visiter2 than once in a while gets nchan more than what is in ms
1815 6154 : Int nchaninms = romscol_p->spectralWindow().numChan()(vb.spectralWindows()(0));
1816 6154 : if(nvischan > nchaninms){
1817 0 : nvischan = nchaninms;}
1818 : //////////////////
1819 6154 : chanMap.resize(nvischan);
1820 6154 : chanMap.set(-1);
1821 :
1822 6154 : Vector<Double> lsrFreq(0);
1823 :
1824 : //cerr << "doConve " << spw << " " << doConversion_p[spw] << " freqframeval " << freqFrameValid_p << endl;
1825 : //cerr <<"valid frame " << freqFrameValid_p << " polmap "<< polMap << endl;
1826 : //cerr << "spectral coord system " << spectralCoord_p.frequencySystem(False) << endl;
1827 6154 : if (freqFrameValid_p &&spectralCoord_p.frequencySystem(False)!=MFrequency::REST ) {
1828 6154 : lsrFreq=vb.getFrequencies(0,MFrequency::LSRK);
1829 : }
1830 : else {
1831 0 : lsrFreq=vb.getFrequencies(0);
1832 : }
1833 6154 : if (spectralCoord_p.frequencySystem(False)==MFrequency::REST && fixMovingSource_p) {
1834 :
1835 0 : if (lastMSId_p != vb.msId()) {
1836 0 : romscol_p=new MSColumns(vb.ms());
1837 : //if ms changed ...reset ephem table
1838 0 : if (upcase(movingDir_p.getRefString()).contains("APP")) {
1839 0 : MeasComet mcomet(Path((romscol_p->field()).ephemPath(vb.fieldId()(0))).absoluteName());
1840 0 : mFrame_p.resetComet(mcomet);
1841 0 : }
1842 : }
1843 0 : MEpoch e0 (Quantity(vb.time()(0), "s"));
1844 0 : mFrame_p.epoch() ? mFrame_p.resetEpoch(e0)
1845 0 : : mFrame_p.set(e0);
1846 :
1847 0 : const auto ephemDirection = vbutil_p->getEphemDir(vb, phaseCenterTime_p);
1848 0 : mFrame_p.direction() ? mFrame_p.resetDirection(ephemDirection)
1849 0 : : mFrame_p.set(ephemDirection);
1850 :
1851 0 : shiftFreqToSource(lsrFreq);
1852 0 : }
1853 : //cerr << "lsrFreq " << lsrFreq.shape() << " nvischan " << nvischan << endl;
1854 : // if(doConversion_p.nelements() < uInt(spw+1))
1855 : // doConversion_p.resize(spw+1, true);
1856 : // doConversion_p[spw]=freqFrameValid_p;
1857 :
1858 6154 : if(lsrFreq.nelements() ==0){
1859 0 : matchPol(vb);
1860 0 : return false;
1861 : }
1862 6154 : lsrFreq_p.resize(lsrFreq.nelements());
1863 6154 : lsrFreq_p=lsrFreq;
1864 6154 : Vector<Double> c(1);
1865 6154 : c=0.0;
1866 6154 : Vector<Double> f(1);
1867 6154 : Int nFound=0;
1868 :
1869 : Double minFreq;
1870 : Double maxFreq;
1871 6154 : spectralCoord_p.toWorld(minFreq, Double(0));
1872 6154 : spectralCoord_p.toWorld(maxFreq, Double(nchan));
1873 6154 : if(maxFreq < minFreq){
1874 0 : f(0)=minFreq;
1875 0 : minFreq=maxFreq;
1876 0 : maxFreq=f(0);
1877 : }
1878 :
1879 :
1880 : //cout.precision(10);
1881 12308 : for (Int chan=0;chan<nvischan;chan++) {
1882 6154 : f(0)=lsrFreq[chan];
1883 6154 : if(spectralCoord_p.toPixel(c, f)) {
1884 6154 : Int pixel=Int(floor(c(0)+0.5)); // round to chan freq at chan center
1885 : //cerr << " chan " << chan << " f " << f(0) << " pixel "<< c(0) << " " << pixel << endl;
1886 : /////////////
1887 : //c(0)=pixel;
1888 : //spectralCoord_p.toWorld(f, c);
1889 : //cerr << "f1 " << f(0) << " pixel "<< c(0) << " " << pixel << endl;
1890 : ////////////////
1891 6154 : if(pixel>-1&&pixel<nchan) {
1892 6154 : chanMap(chan)=pixel;
1893 6154 : nFound++;
1894 6154 : if(nvischan>1&&(chan==0||chan==nvischan-1)) {
1895 : /*logIO() << LogIO::DEBUGGING
1896 : << "Selected visibility channel : " << chan+1
1897 : << " has frequency "
1898 : << MFrequency(Quantity(f(0), "Hz")).get("GHz").getValue()
1899 : << " GHz and maps to image pixel " << pixel+1 << LogIO::POST;
1900 : */
1901 : }
1902 : }
1903 : else{
1904 :
1905 0 : if(nvischan > 1){
1906 0 : Double fwidth=lsrFreq[1]-lsrFreq[0];
1907 0 : Double limit=0;
1908 : //Double where=c(0)*fabs(spectralCoord_p.increment()(0));
1909 0 : if( freqInterpMethod_p==InterpolateArray1D<Double,Complex>::linear)
1910 0 : limit=2;
1911 0 : else if( freqInterpMethod_p==InterpolateArray1D<Double,Complex>::cubic || freqInterpMethod_p==InterpolateArray1D<Double,Complex>::spline)
1912 0 : limit=4;
1913 : //cerr<< "where " << where << " pixel " << pixel << " fwidth " << fwidth << endl;
1914 : /*
1915 : if(((pixel<0) && (where >= (0-limit*fabs(fwidth)))) )
1916 : chanMap(chan)=-2;
1917 : if((pixel>=nchan) ) {
1918 : where=f(0);
1919 : Double fend;
1920 : spectralCoord_p.toWorld(fend, Double(nchan));
1921 : if( ( (fwidth >0) &&where < (fend+limit*fwidth)) || ( (fwidth <0) &&where > (fend+limit*fwidth)) )
1922 : chanMap(chan)=-2;
1923 : }
1924 : */
1925 :
1926 0 : if((f(0) < (maxFreq + limit*fabs(fwidth))) && (f(0) >(maxFreq-0.5*fabs(fwidth)))){
1927 0 : chanMap(chan)=-2;
1928 : }
1929 0 : if((f(0) < minFreq+0.5*fabs(fwidth)) && (f(0) > (minFreq-limit*fabs(fwidth)))){
1930 0 : chanMap(chan)=-2;
1931 : }
1932 : }
1933 :
1934 :
1935 : }
1936 : }
1937 : }
1938 : //cerr << "chanmap " << chanMap << endl;
1939 : /* if(multiChanMap_p.nelements() < uInt(spw+1))
1940 : multiChanMap_p.resize(spw+1);
1941 : multiChanMap_p[spw].resize();
1942 : multiChanMap_p[spw]=chanMap;
1943 : */
1944 :
1945 6154 : if(nFound==0) {
1946 : /*
1947 : logIO() << "Visibility channels in spw " << spw+1
1948 : << " of ms " << vb.msId() << " is not being used "
1949 : << LogIO::WARN << LogIO::POST;
1950 : */
1951 0 : matchPol(vb); ///sometimes the polmap is needed even if chanmap failed
1952 0 : return false;
1953 : }
1954 :
1955 6154 : return matchPol(vb);
1956 :
1957 :
1958 :
1959 6154 : }
1960 :
1961 6154 : Bool FTMachine::matchPol(const vi::VisBuffer2& vb){
1962 6154 : Vector<Stokes::StokesTypes> visPolMap(vb.getCorrelationTypesSelected());
1963 6154 : if((polMap.nelements() > 0) &&(visPolMap.nelements() == visPolMap_p.nelements()) &&allEQ(visPolMap, visPolMap_p))
1964 6120 : return True;
1965 34 : Int stokesIndex=image->coordinates().findCoordinate(Coordinate::STOKES);
1966 34 : AlwaysAssert(stokesIndex>-1, AipsError);
1967 34 : StokesCoordinate stokesCoord=image->coordinates().stokesCoordinate(stokesIndex);
1968 :
1969 :
1970 34 : visPolMap_p.resize();
1971 34 : visPolMap_p=visPolMap;
1972 34 : nvispol=visPolMap.nelements();
1973 34 : AlwaysAssert(nvispol>0, AipsError);
1974 34 : polMap.resize(nvispol);
1975 34 : polMap=-1;
1976 34 : Int pol=0;
1977 34 : Bool found=false;
1978 : // First we try matching Stokes in the visibilities to
1979 : // Stokes in the image that we are gridding into.
1980 102 : for (pol=0;pol<nvispol;pol++) {
1981 68 : Int p=0;
1982 68 : if(stokesCoord.toPixel(p, Stokes::type(visPolMap(pol)))) {
1983 0 : AlwaysAssert(p<npol, AipsError);
1984 0 : polMap(pol)=p;
1985 0 : found=true;
1986 : }
1987 : }
1988 : // If this fails then perhaps we were looking to grid I
1989 : // directly. If so then we need to check that the parallel
1990 : // hands are present in the visibilities.
1991 34 : if(!found) {
1992 34 : Int p=0;
1993 34 : if(stokesCoord.toPixel(p, Stokes::I)) {
1994 34 : polMap=-1;
1995 34 : if(vb.polarizationFrame()==MSIter::Linear) {
1996 34 : p=0;
1997 102 : for (pol=0;pol<nvispol;pol++) {
1998 68 : if(Stokes::type(visPolMap(pol))==Stokes::XX)
1999 34 : {polMap(pol)=0;p++;found=true;};
2000 68 : if(Stokes::type(visPolMap(pol))==Stokes::YY)
2001 34 : {polMap(pol)=0;p++;found=true;};
2002 : }
2003 : }
2004 : else {
2005 0 : p=0;
2006 0 : for (pol=0;pol<nvispol;pol++) {
2007 0 : if(Stokes::type(visPolMap(pol))==Stokes::LL)
2008 0 : {polMap(pol)=0;p++;found=true;};
2009 0 : if(Stokes::type(visPolMap(pol))==Stokes::RR)
2010 0 : {polMap(pol)=0;p++;found=true;};
2011 : }
2012 : }
2013 34 : if(!found) {
2014 : logIO() << "Cannot find polarization map: visibility polarizations = "
2015 0 : << visPolMap << LogIO::EXCEPTION;
2016 : }
2017 : else {
2018 :
2019 : //logIO() << LogIO::DEBUGGING << "Transforming I only" << LogIO::POST;
2020 : }
2021 : };
2022 : }
2023 34 : return True;
2024 6154 : }
2025 :
2026 68 : Vector<String> FTMachine::cleanupTempFiles(const String& mess){
2027 68 : briggsWeightor_p=nullptr;
2028 68 : for(uInt k=0; k < tempFileNames_p.nelements(); ++k){
2029 0 : if(Table::isReadable(tempFileNames_p[k])){
2030 0 : if(mess.size()==0){
2031 : try{
2032 0 : Table::deleteTable(tempFileNames_p[k]);
2033 : }
2034 0 : catch(AipsError &x){
2035 0 : logIO() << LogOrigin("FTMachine", "cleanupTempFiles") << LogIO::NORMAL;
2036 0 : logIO() << LogIO::WARN<< "YOU may have to delete the temporary file " << tempFileNames_p[k] << " because " << x.getMesg() << LogIO::POST;
2037 :
2038 0 : }
2039 : }
2040 : else{
2041 0 : logIO() << LogOrigin("FTMachine", "cleanupTempFiles") << LogIO::NORMAL;
2042 0 : logIO() << "YOU have to delete the temporary file " << tempFileNames_p[k] << " because " << mess << LogIO::DEBUG1 << LogIO::POST;
2043 : }
2044 : }
2045 : }
2046 68 : return tempFileNames_p;
2047 : }
2048 6120 : void FTMachine::gridOk(Int convSupport){
2049 :
2050 6120 : if (nx <= 2*convSupport) {
2051 0 : logIO_p
2052 : << "number of pixels "
2053 : << nx << " on x axis is smaller that the gridding support "
2054 : << 2*convSupport << " Please use a larger value "
2055 0 : << LogIO::EXCEPTION;
2056 : }
2057 :
2058 6120 : if (ny <= 2*convSupport) {
2059 0 : logIO_p
2060 : << "number of pixels "
2061 : << ny << " on y axis is smaller that the gridding support "
2062 : << 2*convSupport << " Please use a larger value "
2063 0 : << LogIO::EXCEPTION;
2064 : }
2065 :
2066 6120 : }
2067 :
2068 0 : void FTMachine::setLocation(const MPosition& loc){
2069 :
2070 0 : mLocation_p=loc;
2071 :
2072 0 : }
2073 :
2074 0 : MPosition& FTMachine::getLocation(){
2075 :
2076 0 : return mLocation_p;
2077 : }
2078 :
2079 :
2080 0 : void FTMachine::setMovingSource(const String& sname, const String& ephtab){
2081 0 : String sourcename=sname;
2082 0 : String ephemtab=ephtab;
2083 : //if a table is given as sourcename...assume ephemerides
2084 0 : if(Table::isReadable(sourcename, False)){
2085 0 : sourcename="COMET";
2086 0 : ephemtab=sname;
2087 0 : ephemTableName_p = sname;
2088 : }
2089 : ///Special case
2090 0 : if(upcase(sourcename)=="TRACKFIELD"){
2091 : //if(name().contains("MosaicFT"))
2092 : // throw(AipsError("Cannot use field phasecenter to track moving source in a Mosaic"));
2093 0 : fixMovingSource_p=True;
2094 0 : movingDir_p=MDirection(Quantity(0.0,"deg"), Quantity(90.0, "deg"));
2095 0 : movingDir_p.setRefString("APP");
2096 : ///Setting it to APP with movingDir_p==True => should use the phasecenter to track
2097 : ///Discussion in CAS-9004 where users are too lazy to give an ephemtable.
2098 0 : return;
2099 : }
2100 :
2101 : MDirection::Types refType;
2102 0 : Bool isValid = MDirection::getType(refType, sourcename);
2103 0 : if(!isValid)
2104 0 : throw(AipsError("Cannot recognize moving source "+sourcename));
2105 0 : if(refType < MDirection::N_Types || refType > MDirection:: N_Planets )
2106 0 : throw(AipsError(sourcename+" is not type of source we can track"));
2107 0 : if(refType==MDirection::COMET){
2108 0 : MeasComet laComet;
2109 0 : if(Table::isReadable(ephemtab, False)){
2110 0 : Table laTable(ephemtab);
2111 0 : Path leSentier(ephemtab);
2112 0 : laComet=MeasComet(laTable, leSentier.absoluteName());
2113 0 : }
2114 : else{
2115 0 : laComet= MeasComet(ephemtab);
2116 : }
2117 0 : if(!mFrame_p.comet())
2118 0 : mFrame_p.set(laComet);
2119 : else
2120 0 : mFrame_p.resetComet(laComet);
2121 0 : }
2122 0 : fixMovingSource_p=true;
2123 0 : movingDir_p=MDirection(Quantity(0.0,"deg"), Quantity(90.0, "deg"));
2124 0 : movingDir_p.setRefString(sourcename);
2125 : // cerr << "movingdir " << movingDir_p.toString() << endl;
2126 0 : }
2127 :
2128 :
2129 0 : void FTMachine::setMovingSource(const MDirection& mdir){
2130 :
2131 0 : fixMovingSource_p=true;
2132 0 : movingDir_p=mdir;
2133 :
2134 0 : }
2135 :
2136 34 : void FTMachine::setFreqInterpolation(const String& method){
2137 :
2138 34 : String meth=method;
2139 34 : meth.downcase();
2140 34 : if(meth.contains("linear")){
2141 34 : freqInterpMethod_p=InterpolateArray1D<Double,Complex>::linear;
2142 : }
2143 0 : else if(meth.contains("splin")){
2144 0 : freqInterpMethod_p=InterpolateArray1D<Double,Complex>::spline;
2145 : }
2146 0 : else if(meth.contains("cub")){
2147 0 : freqInterpMethod_p=InterpolateArray1D<Double,Complex>::cubic;
2148 : }
2149 : else{
2150 0 : freqInterpMethod_p=InterpolateArray1D<Double,Complex>::nearestNeighbour;
2151 : }
2152 :
2153 34 : }
2154 0 : void FTMachine::setFreqInterpolation(const InterpolateArray1D<Double,Complex>::InterpolationMethod type){
2155 0 : freqInterpMethod_p=type;
2156 0 : }
2157 :
2158 : // helper function to swap the y and z axes of a Cube
2159 0 : void FTMachine::swapyz(Cube<Complex>& out, const Cube<Complex>& in)
2160 : {
2161 0 : IPosition inShape=in.shape();
2162 0 : uInt nxx=inShape(0),nyy=inShape(2),nzz=inShape(1);
2163 : //resize breaks references...so out better have the right shape
2164 : //if references is not to be broken
2165 0 : if(out.nelements()==0)
2166 0 : out.resize(nxx,nyy,nzz);
2167 : Bool deleteIn,deleteOut;
2168 0 : const Complex* pin = in.getStorage(deleteIn);
2169 0 : Complex* pout = out.getStorage(deleteOut);
2170 0 : uInt i=0, zOffset=0;
2171 0 : for (uInt iz=0; iz<nzz; ++iz, zOffset+=nxx) {
2172 0 : Int yOffset=zOffset;
2173 0 : for (uInt iy=0; iy<nyy; ++iy, yOffset+=nxx*nzz) {
2174 0 : for (uInt ix=0; ix<nxx; ++ix){
2175 0 : pout[i++] = pin[ix+yOffset];
2176 : }
2177 : }
2178 : }
2179 0 : out.putStorage(pout,deleteOut);
2180 0 : in.freeStorage(pin,deleteIn);
2181 0 : }
2182 :
2183 0 : void FTMachine::swapyz(Cube<Complex>& out, const Cube<Bool>& outFlag, const Cube<Complex>& in)
2184 : {
2185 0 : IPosition inShape=in.shape();
2186 0 : uInt nxx=inShape(0),nyy=inShape(2),nzz=inShape(1);
2187 : //resize breaks references...so out better have the right shape
2188 : //if references is not to be broken
2189 0 : if(out.nelements()==0)
2190 0 : out.resize(nxx,nyy,nzz);
2191 : Bool deleteIn,deleteOut, delFlag;
2192 0 : const Complex* pin = in.getStorage(deleteIn);
2193 0 : const Bool* poutflag= outFlag.getStorage(delFlag);
2194 0 : Complex* pout = out.getStorage(deleteOut);
2195 0 : uInt i=0, zOffset=0;
2196 0 : for (uInt iz=0; iz<nzz; ++iz, zOffset+=nxx) {
2197 0 : Int yOffset=zOffset;
2198 0 : for (uInt iy=0; iy<nyy; ++iy, yOffset+=nxx*nzz) {
2199 0 : for (uInt ix=0; ix<nxx; ++ix){
2200 0 : if(!poutflag[i])
2201 0 : pout[i] = pin[ix+yOffset];
2202 0 : ++i;
2203 : }
2204 : }
2205 : }
2206 0 : out.putStorage(pout,deleteOut);
2207 0 : in.freeStorage(pin,deleteIn);
2208 0 : outFlag.freeStorage(poutflag, delFlag);
2209 0 : }
2210 :
2211 : // helper function to swap the y and z axes of a Cube
2212 0 : void FTMachine::swapyz(Cube<Bool>& out, const Cube<Bool>& in)
2213 : {
2214 0 : IPosition inShape=in.shape();
2215 0 : uInt nxx=inShape(0),nyy=inShape(2),nzz=inShape(1);
2216 0 : if(out.nelements()==0)
2217 0 : out.resize(nxx,nyy,nzz);
2218 : Bool deleteIn,deleteOut;
2219 0 : const Bool* pin = in.getStorage(deleteIn);
2220 0 : Bool* pout = out.getStorage(deleteOut);
2221 0 : uInt i=0, zOffset=0;
2222 0 : for (uInt iz=0; iz<nzz; iz++, zOffset+=nxx) {
2223 0 : Int yOffset=zOffset;
2224 0 : for (uInt iy=0; iy<nyy; iy++, yOffset+=nxx*nzz) {
2225 0 : for (uInt ix=0; ix<nxx; ix++) pout[i++] = pin[ix+yOffset];
2226 : }
2227 : }
2228 0 : out.putStorage(pout,deleteOut);
2229 0 : in.freeStorage(pin,deleteIn);
2230 0 : }
2231 :
2232 0 : void FTMachine::setPointingDirColumn(const String& column){
2233 0 : pointingDirCol_p=column;
2234 0 : pointingDirCol_p.upcase();
2235 0 : if( (pointingDirCol_p != "DIRECTION") &&(pointingDirCol_p != "TARGET") && (pointingDirCol_p != "ENCODER") && (pointingDirCol_p != "POINTING_OFFSET") && (pointingDirCol_p != "SOURCE_OFFSET")){
2236 :
2237 : //basically at this stage you don't know what you're doing...so you get the default
2238 :
2239 0 : pointingDirCol_p="DIRECTION";
2240 :
2241 : }
2242 0 : }
2243 :
2244 0 : String FTMachine::getPointingDirColumnInUse(){
2245 :
2246 0 : return pointingDirCol_p;
2247 :
2248 : }
2249 :
2250 0 : void FTMachine::setSpwChanSelection(const Cube<Int>& spwchansels) {
2251 0 : spwChanSelFlag_p.resize();
2252 0 : spwChanSelFlag_p=spwchansels;
2253 0 : }
2254 :
2255 0 : void FTMachine::setSpwFreqSelection(const Matrix<Double>& spwFreqs)
2256 : {
2257 0 : spwFreqSel_p.assign(spwFreqs);
2258 0 : SynthesisUtils::expandFreqSelection(spwFreqs,expandedSpwFreqSel_p, expandedSpwConjFreqSel_p);
2259 0 : }
2260 :
2261 6120 : void FTMachine::setSpectralFlag(const vi::VisBuffer2& vb, Cube<Bool>& modflagcube){
2262 :
2263 6120 : modflagcube.resize(vb.flagCube().shape());
2264 6120 : modflagcube=vb.flagCube();
2265 6120 : if(!isPseudoI_p){
2266 6120 : ArrayIterator<Bool> from(modflagcube, IPosition(1, 0));
2267 740520 : while(!from.pastEnd()){
2268 734400 : if(anyTrue(from.array()))
2269 0 : from.array().set(True);
2270 734400 : from.next();
2271 : }
2272 6120 : }
2273 6120 : uInt nchan = vb.nChannels();
2274 6120 : uInt msid = vb.msId();
2275 6120 : uInt selspw = vb.spectralWindows()[0];
2276 6120 : Bool spwFlagIsSet=( (spwChanSelFlag_p.nelements() > 0) && (spwChanSelFlag_p.shape()(1) > selspw) &&
2277 6120 : (spwChanSelFlag_p.shape()(0) > msid) &&
2278 0 : (spwChanSelFlag_p.shape()(2) >=nchan));
2279 : //cerr << "spwFlagIsSet " << spwFlagIsSet << endl;
2280 12240 : for (uInt i=0;i<nchan;i++) {
2281 : //Flag those channels that did not get selected...
2282 : //respect the flags from vb if selected or
2283 : //if spwChanSelFlag is wrong shape
2284 6120 : if ((spwFlagIsSet) && (spwChanSelFlag_p(msid,selspw,i)!=1)) {
2285 0 : modflagcube.xzPlane(i).set(true);
2286 : }
2287 : }
2288 6120 : }
2289 :
2290 6120 : Matrix<Double> FTMachine::negateUV(const vi::VisBuffer2& vb){
2291 6120 : Matrix<Double> uvw(vb.uvw().shape());
2292 740520 : for (rownr_t i=0;i< vb.nRows() ; ++i) {
2293 2203200 : for (Int idim=0;idim<2; ++idim) uvw(idim,i)=-vb.uvw()(idim, i);
2294 734400 : uvw(2,i)=vb.uvw()(2,i);
2295 : }
2296 6120 : return uvw;
2297 0 : }
2298 : //-----------------------------------------------------------------------------------------------------------------
2299 : //------------ Vectorized versions of initializeToVis, initializeToSky, finalizeToSky
2300 : //------------ that are called from CubeSkyEquation.
2301 : //------------ They call getImage,getWeightImage, which are implemented in all FTMs
2302 : //------------ Also, Correlation / Stokes conversions and gS/ggS normalizations.
2303 :
2304 :
2305 0 : void FTMachine::setSkyJones(Vector<CountedPtr<casa::refim::SkyJones> >& sj){
2306 0 : sj_p.resize();
2307 0 : sj_p=sj;
2308 0 : cout << "FTM : Set Sky Jones of length " << sj_p.nelements() << " and types ";
2309 0 : for( uInt i=0;i<sj_p.nelements();i++) cout << sj_p[i]->type() << " ";
2310 0 : cout << endl;
2311 0 : }
2312 : // Convert complex correlation planes to float Stokes planes
2313 34 : void FTMachine::correlationToStokes(ImageInterface<Complex>& compImage,
2314 : ImageInterface<Float>& resImage,
2315 : const Bool dopsf)
2316 : {
2317 : // Convert correlation image to IQUV format
2318 34 : AlwaysAssert(compImage.shape()[0]==resImage.shape()[0], AipsError);
2319 34 : AlwaysAssert(compImage.shape()[1]==resImage.shape()[1], AipsError);
2320 34 : AlwaysAssert(compImage.shape()[3]==resImage.shape()[3], AipsError);
2321 :
2322 34 : if(dopsf)
2323 : {
2324 : // For the PSF, choose only those stokes planes that have a valid PSF
2325 17 : StokesImageUtil::ToStokesPSF(resImage,compImage);
2326 : }
2327 : else
2328 : {
2329 17 : StokesImageUtil::To(resImage,compImage);
2330 : }
2331 34 : };
2332 :
2333 : // Convert float Stokes planes to complex correlation planes
2334 0 : void FTMachine::stokesToCorrelation(ImageInterface<Float>& modelImage,
2335 : ImageInterface<Complex>& compImage)
2336 : {
2337 : /*
2338 : StokesCoordinate stcomp=compImage.coordinates().stokesCoordinate(compImage.coordinates().findCoordinate(Coordinate::STOKES));
2339 : StokesCoordinate stfloat = modelImage.coordinates().stokesCoordinate(modelImage.coordinates().findCoordinate(Coordinate::STOKES));
2340 :
2341 : cout << "Stokes types : complex : " << stcomp.stokes() << " float : " << stfloat.stokes() << endl;
2342 : cout << "Shapes : complex : " << compImage.shape() << " float : " << modelImage.shape() << endl;
2343 : */
2344 :
2345 : //Pol axis need not be same
2346 0 : AlwaysAssert(modelImage.shape()[0]==compImage.shape()[0], AipsError);
2347 0 : AlwaysAssert(modelImage.shape()[1]==compImage.shape()[1], AipsError);
2348 0 : AlwaysAssert(modelImage.shape()[3]==compImage.shape()[3], AipsError);
2349 : // Convert from Stokes to Complex
2350 0 : StokesImageUtil::From(compImage, modelImage);
2351 0 : };
2352 :
2353 : //------------------------------------------------------------------------------------------------------------------
2354 :
2355 0 : void FTMachine::normalizeImage(ImageInterface<Float>& skyImage,
2356 : Matrix<Float>& sumOfWts,
2357 : ImageInterface<Float>& sensitivityImage,
2358 : Bool dopsf, Float pblimit, Int normtype)
2359 : {
2360 :
2361 : //Normalize the sky Image
2362 0 : Int nXX=(skyImage).shape()(0);
2363 0 : Int nYY=(skyImage).shape()(1);
2364 0 : Int npola= (skyImage).shape()(2);
2365 0 : Int nchana= (skyImage).shape()(3);
2366 :
2367 0 : IPosition pcentre(4,nXX/2,nYY/2,0,0);
2368 : // IPosition psource(4,nXX/2+22,nYY/2,0,0);
2369 :
2370 : // storeImg(String("norm_resimage.im") , skyImage);
2371 : // storeImg(String("norm_sensitivity.im"), sensitivityImage);
2372 :
2373 : ///// cout << "FTM::norm : pblimit : " << pblimit << endl;
2374 :
2375 : // Note : This is needed because initial prediction has no info about sumwt.
2376 : // Not a clean solution. // ForSB -- if you see a better way, go for it.
2377 0 : if(sumOfWts.shape() != IPosition(2,npola,nchana))
2378 : {
2379 0 : cout << "Empty SumWt.... resizing and filling with 1.0" << endl;
2380 0 : sumOfWts.resize(IPosition(2,npola,nchana));
2381 0 : sumOfWts=1.0;
2382 : }
2383 :
2384 : // if(dopsf) cout << "*** FTM::normalizeImage : Image Center : " << skyImage.getAt(pcentre) << " and weightImage : " << sensitivityImage.getAt(pcentre) << " SumWt : " << sumOfWts[0,0];
2385 : // else cout << "*** FTM::normalizeImage : Source Loc : " << skyImage.getAt(psource) << " and weightImage : " << sensitivityImage.getAt(psource) << " SumWt : " << sumOfWts[0,0];
2386 :
2387 :
2388 :
2389 0 : IPosition blc(4,nXX, nYY, npola, nchana);
2390 0 : IPosition trc(4, nXX, nYY, npola, nchana);
2391 0 : blc(0)=0; blc(1)=0; trc(0)=nXX-1; trc(1)=nYY-1;
2392 : //max weights per plane
2393 0 : for (Int pol=0; pol < npola; ++pol){
2394 0 : for (Int chan=0; chan < nchana ; ++chan){
2395 :
2396 0 : blc(2)=pol; trc(2)=pol;
2397 0 : blc(3)=chan; trc(3)=chan;
2398 0 : Slicer sl(blc, trc, Slicer::endIsLast);
2399 0 : SubImage<Float> subSkyImage(skyImage, sl, false);
2400 0 : SubImage<Float> subSensitivityImage(sensitivityImage, sl, false);
2401 0 : SubImage<Float> subOutput(skyImage, sl, true);
2402 0 : Float sumWt = sumOfWts(pol,chan);
2403 0 : if(sumWt > 0.0){
2404 0 : switch(normtype)
2405 : {
2406 0 : case 0: // only sum Of Weights - FTM only (ForSB)
2407 0 : subOutput.copyData( (LatticeExpr<Float>)
2408 0 : ((dopsf?1.0:-1.0)*subSkyImage/(sumWt)));
2409 0 : break;
2410 :
2411 0 : case 1: // only sensitivityImage Ic/avgPB (ForSB)
2412 0 : subOutput.copyData( (LatticeExpr<Float>)
2413 0 : (iif(subSensitivityImage > (pblimit),
2414 0 : (subSkyImage/(subSensitivityImage)),
2415 : (subSkyImage))));
2416 : // 0.0)));
2417 0 : break;
2418 :
2419 0 : case 2: // sum of Weights and sensitivityImage IGridded/(SoW*avgPB) and PSF --> Id (ForSB)
2420 0 : subOutput.copyData( (LatticeExpr<Float>)
2421 0 : (iif(subSensitivityImage > (pblimit),
2422 0 : ((dopsf?1.0:-1.0)*subSkyImage/(subSensitivityImage*sumWt)),
2423 : //((dopsf?1.0:-1.0)*subSkyImage))));
2424 : 0.0)));
2425 0 : break;
2426 :
2427 0 : case 3: // MULTIPLY by the sensitivityImage avgPB
2428 0 : subOutput.copyData( (LatticeExpr<Float>) (subSkyImage * subSensitivityImage) );
2429 0 : break;
2430 :
2431 0 : case 4: // DIVIDE by sqrt of sensitivityImage
2432 0 : subOutput.copyData( (LatticeExpr<Float>)
2433 0 : (iif((subSensitivityImage) > (pblimit),
2434 0 : (subSkyImage/(sqrt(subSensitivityImage))),
2435 : (subSkyImage))));
2436 : //0.0)));
2437 0 : break;
2438 :
2439 0 : case 5: // MULTIPLY by sqrt of sensitivityImage
2440 0 : subOutput.copyData( (LatticeExpr<Float>)
2441 0 : (iif((subSensitivityImage) > (pblimit),
2442 0 : (subSkyImage * (sqrt(subSensitivityImage))),
2443 : (subSkyImage))));
2444 :
2445 0 : break;
2446 :
2447 0 : case 6: // divide by non normalized sensitivity image
2448 : {
2449 0 : Float elpblimit=max( subSensitivityImage).getFloat() * pblimit;
2450 0 : subOutput.copyData( (LatticeExpr<Float>)
2451 0 : (iif(subSensitivityImage > (elpblimit),
2452 0 : ((dopsf?1.0:-1.0)*subSkyImage/(subSensitivityImage)),
2453 : 0.0)));
2454 : }
2455 0 : break;
2456 0 : default:
2457 0 : throw(AipsError("Unrecognized norm-type in FTM::normalizeImage : "+String::toString(normtype)));
2458 : }
2459 : }
2460 : else{
2461 0 : subOutput.set(0.0);
2462 : }
2463 0 : }
2464 : }
2465 :
2466 : //if(dopsf) cout << " Normalized (" << normtype << ") Image Center : " << skyImage.getAt(pcentre) << endl;
2467 : // else cout << " Normalized (" << normtype << ") Source Loc : " << skyImage.getAt(psource) << endl;
2468 :
2469 0 : }
2470 :
2471 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2472 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2473 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2474 : ///// For use with the new framework
2475 : ///// (Sorry about these copies, but need to keep old system working)
2476 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2477 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2478 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2479 :
2480 : // Vectorized InitializeToVis
2481 0 : void FTMachine::initializeToVisNew(const VisBuffer2& vb,
2482 : CountedPtr<SIImageStore> imstore)
2483 : {
2484 0 : AlwaysAssert(imstore->getNTaylorTerms(false)==1, AipsError);
2485 :
2486 0 : Matrix<Float> tempWts;
2487 :
2488 0 : if(!(imstore->forwardGrid()).get())
2489 0 : throw(AipsError("FTMAchine::InitializeToVisNew error imagestore has no valid grid initialized"));
2490 : // Convert from Stokes planes to Correlation planes
2491 0 : LatticeLocker lock1 (*(imstore->model()), FileLocker::Read);
2492 0 : stokesToCorrelation(*(imstore->model()), *(imstore->forwardGrid()));
2493 :
2494 0 : if(vb.polarizationFrame()==MSIter::Linear) {
2495 0 : StokesImageUtil::changeCStokesRep(*(imstore->forwardGrid()),
2496 : StokesImageUtil::LINEAR);
2497 : }
2498 : else {
2499 0 : StokesImageUtil::changeCStokesRep(*(imstore->forwardGrid()),
2500 : StokesImageUtil::CIRCULAR);
2501 : }
2502 :
2503 : //------------------------------------------------------------------------------------
2504 : // Image Mosaic only : Multiply the input model with the Primary Beam
2505 0 : if(sj_p.nelements() >0 ){
2506 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
2507 0 : (sj_p(k))->apply(*(imstore->forwardGrid()), *(imstore->forwardGrid()), vb, 0, true);
2508 : }
2509 : }
2510 : //------------------------------------------------------------------------------------
2511 :
2512 : // Call initializeToVis
2513 0 : initializeToVis(*(imstore->forwardGrid()), vb); // Pure virtual
2514 :
2515 0 : };
2516 :
2517 : // Vectorized finalizeToVis is not implemented because it does nothing and is never called.
2518 :
2519 : // Vectorized InitializeToSky
2520 34 : void FTMachine::initializeToSkyNew(const Bool dopsf,
2521 : const VisBuffer2& vb,
2522 : CountedPtr<SIImageStore> imstore)
2523 :
2524 : {
2525 34 : AlwaysAssert(imstore->getNTaylorTerms(false)==1, AipsError);
2526 :
2527 : // Make the relevant float grid.
2528 : // This is needed mainly for facetting (to set facet shapes), but is harmless for non-facetting.
2529 34 : if( dopsf ) { imstore->psf(); } else { imstore->residual(); }
2530 :
2531 : // Initialize the complex grid (i.e. tell FTMachine what array to use internally)
2532 34 : Matrix<Float> sumWeight;
2533 34 : if(!(imstore->backwardGrid()).get())
2534 0 : throw(AipsError("FTMAchine::InitializeToSkyNew error imagestore has no valid grid initialized"));
2535 34 : initializeToSky(*(imstore->backwardGrid()) , sumWeight , vb);
2536 :
2537 34 : };
2538 :
2539 : // Vectorized finalizeToSky
2540 34 : void FTMachine::finalizeToSkyNew(Bool dopsf,
2541 : const VisBuffer2& vb,
2542 : CountedPtr<SIImageStore> imstore )
2543 : {
2544 : // Check vector lengths.
2545 34 : AlwaysAssert( imstore->getNTaylorTerms(false)==1, AipsError);
2546 :
2547 34 : Matrix<Float> sumWeights;
2548 34 : finalizeToSky();
2549 :
2550 : //------------------------------------------------------------------------------------
2551 : // Straightforward case. No extra primary beams. No image mosaic
2552 34 : if(sj_p.nelements() == 0 )
2553 : {
2554 : // cerr << "TYPEID " << typeid( *(imstore->psf())).name() << " " << typeid(typeid( *(imstore->psf())).name()).name() << endl;
2555 34 : shared_ptr<ImageInterface<Float> > theim=dopsf ? imstore->psf() : imstore->residual();
2556 : //static_cast<decltype(imstore->residual())>(theim)->lock();
2557 34 : { LatticeLocker lock1 (*theim, FileLocker::Write);
2558 34 : correlationToStokes( getImage(sumWeights, false) , *theim, dopsf);
2559 34 : }
2560 34 : theim->unlock();
2561 :
2562 34 : if( (useWeightImage() && dopsf) || isSD() ) {
2563 :
2564 0 : LatticeLocker lock1 (*(imstore->weight()), FileLocker::Write);
2565 0 : getWeightImage( *(imstore->weight()) , sumWeights);
2566 0 : imstore->weight()->unlock();
2567 : //cerr << "FTMachine getweight " << max(imstore->weight()->get()) << endl;
2568 :
2569 :
2570 : // Fill weight image only once, during PSF generation. Remember.... it is normalized only once
2571 : // during PSF generation.
2572 0 : }
2573 :
2574 : // Take sumWeights from corrToStokes here....
2575 34 : LatticeLocker lock1 (*(imstore->sumwt()), FileLocker::Write);
2576 34 : Bool donesumwt=(max(imstore->sumwt()->get()) > 0.0);
2577 34 : if(!donesumwt){
2578 34 : Matrix<Float> sumWeightStokes( (imstore->sumwt())->shape()[2], (imstore->sumwt())->shape()[3] );
2579 17 : CoordinateSystem incoord=image->coordinates();
2580 17 : CoordinateSystem outcoord=imstore->sumwt()->coordinates();
2581 17 : StokesImageUtil::ToStokesSumWt(sumWeightStokes, sumWeights, outcoord, incoord);
2582 :
2583 :
2584 17 : Array<Float> sumWtArr(IPosition(4,1,1,sumWeights.shape()[0], sumWeights.shape()[1]));
2585 :
2586 17 : IPosition blc(4, 0, 0, 0, 0);
2587 17 : IPosition trc(4, 0, 0, sumWeightStokes.shape()[0]-1, sumWeightStokes.shape()[1]-1);
2588 17 : sumWtArr(blc, trc).reform(sumWeightStokes.shape())=sumWeightStokes;
2589 :
2590 : //StokesImageUtil::ToStokesSumWt( sumWeightStokes, sumWeights );
2591 17 : AlwaysAssert( ( (imstore->sumwt())->shape()[2] == sumWeightStokes.shape()[0] ) &&
2592 : ((imstore->sumwt())->shape()[3] == sumWeightStokes.shape()[1] ) , AipsError );
2593 :
2594 17 : (imstore->sumwt())->put( sumWeightStokes.reform((imstore->sumwt())->shape()) );
2595 17 : }
2596 34 : imstore->sumwt()->unlock();
2597 :
2598 34 : }
2599 : //------------------------------------------------------------------------------------
2600 : // Image Mosaic only : Multiply the residual, and weight image by the PB.
2601 : else
2602 : {
2603 :
2604 : // Take the FT of the gridded values. Writes into backwardGrid().
2605 0 : getImage(sumWeights, false);
2606 :
2607 : // Multiply complex image grid by PB.
2608 0 : if( !dopsf )
2609 : {
2610 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
2611 0 : (sj_p(k))->apply(*(imstore->backwardGrid()), *(imstore->backwardGrid()), vb, 0, true);
2612 : }
2613 : }
2614 :
2615 : // Convert from correlation to Stokes onto a new temporary grid.
2616 0 : SubImage<Float> targetImage( ( dopsf ? *(imstore->psf()) : *(imstore->residual()) ) , true);
2617 0 : TempImage<Float> temp( targetImage.shape(), targetImage.coordinates() );
2618 0 : correlationToStokes( *(imstore->backwardGrid()), temp, false);
2619 :
2620 : // Add the temporary Stokes image to the residual or PSF, whichever is being made.
2621 0 : LatticeExpr<Float> addToRes( targetImage + temp );
2622 0 : targetImage.copyData(addToRes);
2623 :
2624 : // Now, do the same with the weight image and sumwt ( only on the first pass )
2625 0 : if( dopsf )
2626 : {
2627 0 : SubImage<Float> weightImage( *(imstore->weight()) , true);
2628 0 : TempImage<Float> temp(weightImage.shape(), weightImage.coordinates());
2629 0 : getWeightImage(temp, sumWeights);
2630 :
2631 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
2632 0 : (sj_p(k))->applySquare(temp,temp, vb, -1);
2633 : }
2634 :
2635 0 : LatticeExpr<Float> addToWgt( weightImage + temp );
2636 0 : weightImage.copyData(addToWgt);
2637 :
2638 0 : AlwaysAssert( ( (imstore->sumwt())->shape()[2] == sumWeights.shape()[0] ) &&
2639 : ((imstore->sumwt())->shape()[3] == sumWeights.shape()[1] ) , AipsError );
2640 :
2641 0 : SubImage<Float> sumwtImage( *(imstore->sumwt()) , true);
2642 0 : TempImage<Float> temp2(sumwtImage.shape(), sumwtImage.coordinates());
2643 0 : temp2.put( sumWeights.reform(sumwtImage.shape()) );
2644 0 : LatticeExpr<Float> addToWgt2( sumwtImage + temp2 );
2645 0 : sumwtImage.copyData(addToWgt2);
2646 :
2647 : //cout << "In finalizeGridCoreMos : sumwt : " << sumwtImage.get() << endl;
2648 :
2649 0 : }
2650 :
2651 0 : }
2652 : //------------------------------------------------------------------------------------
2653 :
2654 :
2655 :
2656 68 : return;
2657 34 : };
2658 :
2659 :
2660 : /////------------------------------------------------
2661 0 : void FTMachine::finalizeToWeightImage(const VisBuffer2& vb,
2662 : CountedPtr<SIImageStore> imstore )
2663 : {
2664 : // Check vector lengths.
2665 0 : AlwaysAssert( imstore->getNTaylorTerms(false)==1, AipsError);
2666 :
2667 0 : Matrix<Float> sumWeights;
2668 :
2669 : //------------------------------------------------------------------------------------
2670 : // Straightforward case. No extra primary beams. No image mosaic
2671 0 : if(sj_p.nelements() == 0 )
2672 : {
2673 :
2674 :
2675 0 : if( useWeightImage() ) {
2676 : //if( name().contains("Mosaic") ){
2677 : {
2678 0 : finalizeToSky();
2679 : }
2680 0 : LatticeLocker lock1 (*(imstore->weight()), FileLocker::Write);
2681 0 : getWeightImage( *(imstore->weight()) , sumWeights);
2682 0 : imstore->weight()->unlock();
2683 :
2684 : // Fill weight image only once, during PSF generation. Remember.... it is normalized only once
2685 : // during PSF generation.
2686 0 : }
2687 0 : if(sumWeights.nelements() >0){
2688 : // Take sumWeights from corrToStokes here....
2689 0 : LatticeLocker lock1 (*(imstore->sumwt()), FileLocker::Write);
2690 0 : Matrix<Float> sumWeightStokes( (imstore->sumwt())->shape()[2], (imstore->sumwt())->shape()[3] );
2691 0 : StokesImageUtil::ToStokesSumWt( sumWeightStokes, sumWeights );
2692 :
2693 0 : AlwaysAssert( ( (imstore->sumwt())->shape()[2] == sumWeightStokes.shape()[0] ) &&
2694 : ((imstore->sumwt())->shape()[3] == sumWeightStokes.shape()[1] ) , AipsError );
2695 :
2696 0 : (imstore->sumwt())->put( sumWeightStokes.reform((imstore->sumwt())->shape()) );
2697 0 : imstore->sumwt()->unlock();
2698 0 : }
2699 :
2700 : }
2701 : //------------------------------------------------------------------------------------
2702 : // Image Mosaic only : Multiply the residual, and weight image by the PB.
2703 : else
2704 : {
2705 :
2706 : // Now, do the same with the weight image and sumwt ( only on the first pass )
2707 : {
2708 0 : SubImage<Float> weightImage( *(imstore->weight()) , true);
2709 0 : TempImage<Float> temp(weightImage.shape(), weightImage.coordinates());
2710 0 : getWeightImage(temp, sumWeights);
2711 :
2712 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
2713 0 : (sj_p(k))->applySquare(temp,temp, vb, -1);
2714 : }
2715 :
2716 0 : LatticeExpr<Float> addToWgt( weightImage + temp );
2717 0 : weightImage.copyData(addToWgt);
2718 :
2719 0 : AlwaysAssert( ( (imstore->sumwt())->shape()[2] == sumWeights.shape()[0] ) &&
2720 : ((imstore->sumwt())->shape()[3] == sumWeights.shape()[1] ) , AipsError );
2721 :
2722 0 : SubImage<Float> sumwtImage( *(imstore->sumwt()) , true);
2723 0 : TempImage<Float> temp2(sumwtImage.shape(), sumwtImage.coordinates());
2724 0 : temp2.put( sumWeights.reform(sumwtImage.shape()) );
2725 0 : LatticeExpr<Float> addToWgt2( sumwtImage + temp2 );
2726 0 : sumwtImage.copyData(addToWgt2);
2727 :
2728 : //cout << "In finalizeGridCoreMos : sumwt : " << sumwtImage.get() << endl;
2729 :
2730 0 : }
2731 :
2732 : }
2733 : //------------------------------------------------------------------------------------
2734 :
2735 :
2736 :
2737 0 : return;
2738 0 : };
2739 :
2740 :
2741 :
2742 :
2743 : /////-----------------------------------------------
2744 0 : Bool FTMachine::changedSkyJonesLogic(const vi::VisBuffer2& vb, Bool& firstRow, Bool& internalRow)
2745 : {
2746 0 : firstRow=false;
2747 0 : internalRow=false;
2748 :
2749 0 : if( sj_p.nelements()==0 )
2750 0 : {throw(AipsError("Internal Error : Checking changedSkyJones, but it is not yet set."));}
2751 :
2752 0 : CountedPtr<SkyJones> ej = sj_p[0];
2753 0 : if(ej.null())
2754 0 : return false;
2755 0 : if(ej->changed(vb,0))
2756 0 : firstRow=true;
2757 0 : Int row2temp=0;
2758 0 : if(ej->changedBuffer(vb,0,row2temp)) {
2759 0 : internalRow=true;
2760 : }
2761 0 : return (firstRow || internalRow) ;
2762 0 : }
2763 :
2764 0 : std::shared_ptr<std::complex<double>> FTMachine::getGridPtr(size_t& size) const
2765 : {
2766 0 : size = 0;
2767 0 : return std::shared_ptr<std::complex<double>>();
2768 : }
2769 :
2770 0 : std::shared_ptr<double> FTMachine::getSumWeightsPtr(size_t& size) const
2771 : {
2772 0 : size = 0;
2773 0 : return std::shared_ptr<double>();
2774 : }
2775 :
2776 0 : void FTMachine::setCFCache(CountedPtr<CFCache>& /*cfc*/, const Bool /*loadCFC*/)
2777 : {
2778 0 : throw(AipsError("FTMachine::setCFCache() directly called!"));
2779 : }
2780 :
2781 :
2782 :
2783 1088 : void FTMachine::findGridSector(const Int& nxp, const Int& nyp, const Int& ixsub, const Int& iysub, const Int& minx, const Int& miny, const Int& icounter, Int& x0, Int& y0, Int& nxsub, Int& nysub, const Bool linear){
2784 : /* Vector<Int> ord(36);
2785 : ord(0)=14;
2786 : ord(1)=15;
2787 : ord(2)=20;
2788 : ord(3)=21;ord(4)=13;
2789 : ord(5)=16;ord(6)=19;ord(7)=22;ord(8)=8;ord(9)=9;
2790 : ord(10)=26;ord(11)=27;ord(12)=25;ord(13)=28;ord(14)=7;
2791 : ord(15)=10;ord(16)=32;ord(17)=33;ord(18)=2;ord(19)=3;
2792 : ord(20)=18;ord(21)=23;ord(22)=12;ord(23)=17;ord(24)=1;
2793 : ord(25)=4;ord(26)=6;ord(27)=11;ord(28)=24;ord(29)=29;
2794 : ord(30)=31;ord(31)=34;ord(32)=0;ord(33)=5;ord(34)=30;
2795 : ord(35)=35;
2796 : */
2797 : /*
2798 : Int ix= (icounter+1)-((icounter)/ixsub)*ixsub;
2799 : Int iy=(icounter)/ixsub+1;
2800 : y0=(nyp/iysub)*(iy-1)+1+miny;
2801 : nysub=nyp/iysub;
2802 : if( iy == iysub) {
2803 : nysub=nyp-(nyp/iysub)*(iy-1);
2804 : }
2805 : x0=(nxp/ixsub)*(ix-1)+1+minx;
2806 : nxsub=nxp/ixsub;
2807 : if( ix == ixsub){
2808 : nxsub=nxp-(nxp/ixsub)*(ix-1);
2809 : }
2810 : */
2811 :
2812 :
2813 : {
2814 1088 : Int elrow=icounter/ixsub;
2815 1088 : Int elcol=(icounter-elrow*ixsub);
2816 : //cerr << "row "<< elrow << " col " << elcol << endl;
2817 : //nxsub=Int(floor(((ceil(fabs(float(2*elcol+1-ixsub)/2.0))-1.0)*5 +1)*nxp/36.0 + 0.5));
2818 1088 : Float factor=0;
2819 1088 : if(ixsub > 1){
2820 5440 : for (Int k=0; k < ixsub/2; ++k)
2821 4352 : factor= linear ? factor+(k+1): factor+sqrt(Float(k+1));
2822 : //factor= linear ? factor+(k+1): factor+(k+1)*(k+1)*(k+1);
2823 1088 : factor *= 2.0;
2824 1088 : if(linear)
2825 1088 : nxsub=Int(floor((ceil(fabs(float(2*elcol+1-ixsub)/2.0))/factor)*nxp + 0.5));
2826 : else
2827 : //nxsub=Int(floor((ceil(fabs(float(2*elcol+1-ixsub)/2.0))*ceil(fabs(float(2*elcol+1-ixsub)/2.0))*ceil(fabs(float(2*elcol+1-ixsub)/2.0))/factor)*nxp + 0.5));
2828 0 : nxsub=Int(floor((sqrt(ceil(fabs(float(2*elcol+1-ixsub)/2.0)))/factor)*nxp + 0.5));
2829 : }
2830 : else{
2831 0 : nxsub=nxp;
2832 : }
2833 : //cerr << nxp << " col " << elcol << " nxsub " << nxsub << endl;
2834 1088 : x0=minx;
2835 1088 : elcol-=1;
2836 4896 : while(elcol >= 0){
2837 : //x0+=Int(floor(((ceil(fabs(float(2*elcol+1-ixsub)/2.0))-1.0)*5 +1)*nxp/36.0+0.5));
2838 :
2839 3808 : if(linear)
2840 3808 : x0+=Int(floor((ceil(fabs(float(2*elcol+1-ixsub)/2.0))/factor)*nxp + 0.5));
2841 : else
2842 : //x0+=Int(floor((ceil(fabs(float(2*elcol+1-ixsub)/2.0))*ceil(fabs(float(2*elcol+1-ixsub)/2.0))*ceil(fabs(float(2*elcol+1-ixsub)/2.0))/factor)*nxp + 0.5));
2843 0 : x0+=Int(floor((sqrt(ceil(fabs(float(2*elcol+1-ixsub)/2.0)))/factor)*nxp + 0.5));
2844 3808 : elcol-=1;
2845 : }
2846 1088 : factor=0;
2847 1088 : if(iysub >1){
2848 5440 : for (Int k=0; k < iysub/2; ++k)
2849 : //factor=linear ? factor+(k+1): factor+(k+1)*(k+1)*(k+1);
2850 4352 : factor= linear ? factor+(k+1): factor+sqrt(Float(k+1));
2851 1088 : factor *= 2.0;
2852 : //nysub=Int(floor(((ceil(fabs(float(2*elrow+1-iysub)/2.0))-1.0)*5 +1)*nyp/36.0+0.5));
2853 1088 : if(linear)
2854 1088 : nysub=Int(floor((ceil(fabs(float(2*elrow+1-iysub)/2.0))/factor)*nyp + 0.5));
2855 : else
2856 0 : nysub=Int(floor((sqrt(ceil(fabs(float(2*elrow+1-iysub)/2.0)))/factor)*nyp + 0.5));
2857 : }
2858 : else{
2859 0 : nysub=nyp;
2860 : }
2861 : //nysub=Int(floor((ceil(fabs(float(2*elrow+1-iysub)/2.0))*ceil(fabs(float(2*elrow+1-iysub)/2.0))*ceil(fabs(float(2*elrow+1-iysub)/2.0))/factor)*nyp + 0.5));
2862 1088 : y0=miny;
2863 1088 : elrow-=1;
2864 :
2865 4896 : while(elrow >=0){
2866 : //y0+=Int(floor(((ceil(fabs(float(2*elrow+1-iysub)/2.0))-1.0)*5 +1)*nyp/36.0+0.5));
2867 3808 : if(linear)
2868 3808 : y0+=Int(floor((ceil(fabs(float(2*elrow+1-iysub)/2.0))/factor)*nyp + 0.5));
2869 : else
2870 0 : y0+=Int(floor((sqrt(ceil(fabs(float(2*elrow+1-iysub)/2.0)))/factor)*nyp + 0.5));
2871 : //y0+=Int(floor((ceil(fabs(float(2*elrow+1-iysub)/2.0))*ceil(fabs(float(2*elrow+1-iysub)/2.0))*ceil(fabs(float(2*elrow+1-iysub)/2.0))/factor)*nyp + 0.5));
2872 3808 : elrow-=1;
2873 : }
2874 : }
2875 1088 : if( (y0+nysub) >= nyp){
2876 136 : nysub = nyp-y0-1;
2877 : }
2878 1088 : if( (x0+nxsub) >= nxp){
2879 136 : nxsub = nxp-x0-1;
2880 : }
2881 1088 : y0+=1;
2882 1088 : x0+=1;
2883 : //cerr << icounter << " x0, y0 " << x0 << " " << y0 << " ixsub, iysub " << nxsub << " " << nysub << endl;
2884 1088 : if(doneThreadPartition_p < 0)
2885 17 : doneThreadPartition_p=1;
2886 :
2887 1088 : }
2888 :
2889 3060 : void FTMachine::tweakGridSector(const Int& nx, const Int& ny, const Int& ixsub, const Int& iysub){
2890 : //if(doneThreadPartition_p)
2891 : // return;
2892 3060 : Vector<Int> x0, y0, nxsub, nysub;
2893 3060 : Vector<Float> xcut(nx/2);
2894 3060 : Vector<Float> ycut(ny/2);
2895 3060 : if(griddedData2.nelements() >0 ){
2896 : //cerr << "shapes " << xcut.shape() << " gd " << amplitude(griddedData2(IPosition(4, 0, ny/2-1, 0, 0), IPosition(4, nx/2-1, ny/2-1, 0,0))).shape() << endl;
2897 3060 : convertArray(xcut, amplitude(griddedData2(IPosition(4, 0, ny/2-1, 0, 0), IPosition(4, nx/2-1, ny/2-1, 0,0))).reform(xcut.shape()));
2898 3060 : convertArray(ycut, amplitude(griddedData2(IPosition(4, nx/2-1, 0, 0, 0), IPosition(4, nx/2-1, ny/2-1, 0,0))).reform(ycut.shape()));
2899 : }
2900 : else{
2901 0 : xcut=amplitude(griddedData(IPosition(4, 0, ny/2-1, 0, 0), IPosition(4, nx/2-1, ny/2-1, 0,0)));
2902 0 : ycut=amplitude(griddedData(IPosition(4, nx/2-1, 0, 0, 0), IPosition(4, nx/2-1, ny/2-1, 0,0)));
2903 : }
2904 : //cerr << griddedData2.shape() << " " << griddedData.shape() << endl;
2905 3060 : Vector<Float> cumSumX(nx/2, 0);
2906 : //Vector<Float> cumSumX2(nx/2,0);
2907 3060 : cumSumX(0)=xcut(0);
2908 : //cumSumX2(0)=cumSumX(0)*cumSumX(0);
2909 3060 : Float sumX=sum(xcut);
2910 3060 : if(sumX==0.0)
2911 0 : return;
2912 : //cerr << "sumX " << sumX << endl;
2913 : //sumX *= sumX;
2914 3060 : x0.resize(ixsub);
2915 3060 : x0=nx/2-1;
2916 3060 : nxsub.resize(ixsub);
2917 3060 : nxsub=0;
2918 3060 : x0(0)=0;
2919 3060 : Int counter=1;
2920 244800 : for (Int k=1; k < nx/2; ++k){
2921 241740 : cumSumX(k)=cumSumX(k-1)+xcut(k);
2922 : //cumSumX2(k)=cumSumX(k)*cumSumX(k);
2923 241740 : Float nextEdge=sumX/(Float(ixsub/2)*Float(ixsub/2))*Float(counter)*Float(counter);
2924 241740 : if(cumSumX(k-1) < nextEdge && cumSumX(k) >= nextEdge){
2925 12240 : x0(counter)=k;
2926 : //cerr << counter << " " << k << " diff " << x0(counter)-x0(counter-1) << endl;
2927 12240 : nxsub(counter-1)=x0(counter)-x0(counter-1);
2928 12240 : ++counter;
2929 : }
2930 : }
2931 :
2932 3060 : x0(ixsub/2)=nx/2-1;
2933 3060 : nxsub(ixsub/2)=nxsub(ixsub/2-1)=nx/2-1-x0(ixsub/2-1);
2934 12240 : for(Int k=ixsub/2+1; k < ixsub; ++k){
2935 9180 : x0(k)=x0(k-1)+ nxsub(ixsub-k);
2936 9180 : nxsub(k)=nxsub(ixsub-k-1);
2937 : }
2938 3060 : nxsub(ixsub-1)+=1;
2939 :
2940 3060 : Vector<Float> cumSumY(ny/2, 0);
2941 : //Vector<Float> cumSumY2(ny/2,0);
2942 3060 : cumSumY(0)=ycut(0);
2943 : //cumSumY2(0)=cumSumY(0)*cumSumY(0);
2944 3060 : Float sumY=sum(ycut);
2945 3060 : if(sumY==0.0)
2946 0 : return;
2947 : //sumY *=sumY;
2948 3060 : y0.resize(iysub);
2949 3060 : y0=ny/2-1;
2950 3060 : nysub.resize(iysub);
2951 3060 : nysub=0;
2952 3060 : y0(0)=0;
2953 3060 : counter=1;
2954 244800 : for (Int k=1; k < ny/2; ++k){
2955 241740 : cumSumY(k)=cumSumY(k-1)+ycut(k);
2956 : //cumSumY2(k)=cumSumY(k)*cumSumY(k);
2957 241740 : Float nextEdge=sumY/(Float(iysub/2)*Float(iysub/2))*Float(counter)*Float(counter);
2958 241740 : if(cumSumY(k-1) < nextEdge && cumSumY(k) >= nextEdge){
2959 12240 : y0(counter)=k;
2960 12240 : nysub(counter-1)=y0(counter)-y0(counter-1);
2961 12240 : ++counter;
2962 : }
2963 : }
2964 :
2965 3060 : y0(ixsub/2)=ny/2-1;
2966 3060 : nysub(iysub/2)=nysub(iysub/2-1)=ny/2-1-y0(iysub/2-1);
2967 12240 : for(Int k=iysub/2+1; k < iysub; ++k){
2968 9180 : y0(k)=y0(k-1)+ nysub(iysub-k);
2969 9180 : nysub(k)=nysub(iysub-k-1);
2970 : }
2971 3060 : nysub(iysub-1)+=1;
2972 :
2973 3060 : if(anyEQ(nxsub, 0) || anyEQ(nysub, 0))
2974 0 : return;
2975 : //cerr << " x0 " << x0 << " nxsub " << nxsub << endl;
2976 : //cerr << " y0 " << y0 << " nysub " << nysub << endl;
2977 3060 : x0+=1;
2978 3060 : y0+=1;
2979 3060 : xsect_p.resize(ixsub*iysub);
2980 3060 : ysect_p.resize(ixsub*iysub);
2981 3060 : nxsect_p.resize(ixsub*iysub);
2982 3060 : nysect_p.resize(ixsub*iysub);
2983 27540 : for (Int iy=0; iy < iysub; ++iy){
2984 220320 : for (Int ix=0; ix< ixsub; ++ix){
2985 :
2986 195840 : xsect_p(iy*ixsub+ix)=x0[ix];
2987 195840 : ysect_p(iy*ixsub+ix)=y0[iy];
2988 195840 : nxsect_p(iy*ixsub+ix)=nxsub[ix];
2989 195840 : nysect_p(iy*ixsub+ix)=nysub[iy];
2990 : }
2991 : }
2992 :
2993 3060 : ++doneThreadPartition_p;
2994 :
2995 3060 : }
2996 :
2997 :
2998 : /*
2999 : /// Move to individual FTMs............ make it pure virtual.
3000 : Bool FTMachine::useWeightImage()
3001 : {
3002 : if( name() == "GridFT" || name() == "WProjectFT" )
3003 : { return false; }
3004 : else
3005 : { return true; }
3006 : }
3007 : */
3008 :
3009 : }//# namespace refim ends
3010 : }//namespace CASA ends
3011 :
|
